最近涉及到使用Fortran对微分方程求解,我们知道MATLAB已有内置的函数,比如ode家族,ode15s,对应着不同的求解办法。通过查看odepack的官方文档,我尝试使用了dlsode求解刚性和非刚性常微分方程组。
首先是github网址:https://github.com/jacobwilliams/odepack
具体使用办法:
1.我使用的是vs2022,比较简单的用法就是把,src文件夹所有的文件复制到和项目一个文件夹即可,将M_odepack.f90文件放入到项目中,这样就可以用了。
2.在使用前要use m_odepack
3.这里以官方文档中的例子为例:
program dlsode_ex
use m_odepack
implicit none
external fex
external jex
integer,parameter :: dp=kind(0.0d0)
real(kind=dp),dimension(3) :: atol,y
integer :: iopt,iout,istate,itask,itol,liw,lrw,mf,neq
integer,dimension(23) :: iwork
real(kind=dp) :: rtol,t,tout
real(kind=dp),dimension(58) :: rwork
neq = 3
y(1) = 1.D0
y(2) = 0.D0
y(3) = 0.D0
t = 0.D0
tout = .4D0
itol = 2
rtol = 1.D-4
atol(1) = 1.D-6
atol(2) = 1.D-10
atol(3) = 1.D-6
itask = 1
istate = 1
iopt = 0
lrw = 58
liw = 23
mf = 21
do iout = 1,12
call dlsode(fex,[neq],y,t,tout,itol,[rtol],atol,itask,istate,iopt, &
& rwork,lrw,iwork,liw,jex,mf)
write (6,99010) t,y(1),y(2),y(3)
99010 format (' At t =',d12.4,' y =',3D14.6)
if ( istate<0 ) then
write (6,99020) istate
99020 format (///' Error halt.. ISTATE =',i3)
stop 1
else
tout = tout*10.D0
endif
enddo
write (6,99030) iwork(11),iwork(12),iwork(13)
99030 format (/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4)
end program dlsode_ex
subroutine fex(Neq,T,Y,Ydot)
implicit none
integer,parameter :: dp=kind(0.0d0)
integer :: Neq
real(kind=dp) :: T
real(kind=dp),intent(in),dimension(3) :: Y
real(kind=dp),intent(inout),dimension(3) :: Ydot
Ydot(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)
Ydot(3) = 3.D7*Y(2)*Y(2)
Ydot(2) = -Ydot(1) - Ydot(3)
end subroutine fex
subroutine jex(Neq,T,Y,Ml,Mu,Pd,Nrpd)
implicit none
integer,parameter :: dp=kind(0.0d0)
integer :: Neq
real(kind=dp) :: T
real(kind=dp),intent(in),dimension(3) :: Y
integer :: Ml
integer :: Mu
real(kind=dp),intent(inout),dimension(Nrpd,3) :: Pd
integer,intent(in) :: Nrpd
Pd(1,1) = -.04D0
Pd(1,2) = 1.D4*Y(3)
Pd(1,3) = 1.D4*Y(2)
Pd(2,1) = .04D0
Pd(2,3) = -Pd(1,3)
Pd(3,2) = 6.D7*Y(2)
Pd(2,2) = -Pd(1,2) - Pd(3,2)
end subroutine jex
一些变量意义具体看文档说明:https://jacobwilliams.github.io/odepack/proc/dlsode.html
其中,假设n是方程个数,
y:是初值,数组,y(n)
atol:每个方程的绝对误差,数组,atol(n)
t:输入的初始点,tout是下一个点。
mf:是求解方法,其中如果等于21,24需要使用者自己提供雅各比矩阵,如示例代码中jex函数中那样,如果等于10,22,25则不需要自己写,但是jex函数还是需要定义,就是函数框架,函数名,变量声明就可。
fex函数:写的就是你的微分方程组
另外,
rwork,iwork也是两个一维数组,大小如图所示。
以及,
lrw = 22 + 9*NEQ + NEQ**2
liw = 20 + NEQ
整体使用的逻辑就是先设置t值,然后设置循环,tout不断累加,下次循环就使用上次计算得到的新y值以及tout进行迭代计算。文章来源:https://www.toymoban.com/news/detail-675688.html
istate是用于输入和输出以指定计算状态的索引,要注意的是如果istate选择2,或3需要在第一次循环中等于1,初始化,到了第二次循环开始才赋值为2或3。文章来源地址https://www.toymoban.com/news/detail-675688.html
到了这里,关于Fortran 微分方程求解 --ODEPACK的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!