深圳公布最新出行政策seo优化案例
最近涉及到使用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 jexinteger,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) :: rworkneq = 3y(1) = 1.D0y(2) = 0.D0y(3) = 0.D0t = 0.D0tout = .4D0itol = 2rtol = 1.D-4atol(1) = 1.D-6atol(2) = 1.D-10atol(3) = 1.D-6itask = 1istate = 1iopt = 0lrw = 58liw = 23mf = 21do iout = 1,12call 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 ) thenwrite (6,99020) istate99020 format (///' Error halt.. ISTATE =',i3)stop 1elsetout = tout*10.D0endifenddowrite (6,99030) iwork(11),iwork(12),iwork(13)99030 format (/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4)end program dlsode_exsubroutine 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) :: YdotYdot(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 fexsubroutine jex(Neq,T,Y,Ml,Mu,Pd,Nrpd)
implicit noneinteger,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) :: NrpdPd(1,1) = -.04D0Pd(1,2) = 1.D4*Y(3)Pd(1,3) = 1.D4*Y(2)Pd(2,1) = .04D0Pd(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进行迭代计算。
istate是用于输入和输出以指定计算状态的索引，要注意的是如果istate选择2，或3需要在第一次循环中等于1，初始化，到了第二次循环开始才赋值为2或3。