Fortran 微分方程求解 --ODEPACK

这篇具有很好参考价值的文章主要介绍了Fortran 微分方程求解 --ODEPACK。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

最近涉及到使用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函数:写的就是你的微分方程组

另外,

Fortran 微分方程求解 --ODEPACK,Fortran,算法

Fortran 微分方程求解 --ODEPACK,Fortran,算法

 rwork,iwork也是两个一维数组,大小如图所示。

以及,

lrw = 22 +  9*NEQ + NEQ**2
liw = 20 + NEQ

整体使用的逻辑就是先设置t值,然后设置循环,tout不断累加,下次循环就使用上次计算得到的新y值以及tout进行迭代计算。

istate是用于输入和输出以指定计算状态的索引,要注意的是如果istate选择2,或3需要在第一次循环中等于1,初始化,到了第二次循环开始才赋值为2或3。文章来源地址https://www.toymoban.com/news/detail-675688.html

到了这里,关于Fortran 微分方程求解 --ODEPACK的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • MATLAB求解偏微分方程【PDE和差分法】

    目录 前言  1.用差分法求解 显示差分 其他方程举例 : r是什么 2.PDETOOL 3.pdepe函数 示例:热方程 代码:   在我们处理一些公式时,常常会有偏微分方程出现,所以我今天整理了一下求解偏微分方程的常用方法,希望有所帮助 在1979年复旦大学学者的一篇论文里,谈到了偏微分

    2024年02月04日
    浏览(49)
  • 40.利用欧拉法求解微分方程组(matlab程序)

    1. 简述        求解微分方程的时候,如果不能将求出结果的表达式,则可以对利用数值积分对微分方程求解,获取数值解。欧拉方法是最简单的一种数值解法。前面介绍过MATLAB实例讲解欧拉法求解微分方程,今天实例讲解欧拉法求解一阶微分方程组。 本文理论部分来自知乎

    2024年02月14日
    浏览(57)
  • matlab使用教程(27)—微分代数方程(DAE)求解

            微分代数方程是一类微分方程,其中一个或多个因变量导数未出现在方程中。方程中出现的未包含其导数的变量称为代数变量,代数变量的存在意味着您不能将这些方程记为显式形式 y ′ = f t , y 。相反,您可以解算下列形式的 DAE:         • ode15s 和 ode23t

    2024年02月11日
    浏览(46)
  • Simulink基础【1】-弹簧-阻尼模型的常微分方程求解

    Simulink是Matlab软件的框图设计环境,可用于各种动态系统的建模、分析与仿真过程。如:导航制导、通讯、电子、机械、热力学等诸多领域。这些系统在数学角度描述上涉及连续、离散、非线性、时变等用解析方法难以求解的系统,因而采用Simulink进行建模与仿真是指导这些系

    2023年04月08日
    浏览(47)
  • matlab使用教程(28)—微分方程(ODE)求解常见问题

            本博客说明如何将 ODE 解约束为非负解。施加非负约束不一定总是可有可无,在某些情况下,由于方程的物理解释或解性质的原因,可能有必要施加非负约束。仅在必要时对解施加此约束,例如不这样做积分就会失败或者解将不适用的情况。         如果解的

    2024年02月11日
    浏览(54)
  • Python和Julia TensorFlow科学计算常微分方程求解器

    常微分方程(ODE)可用于描述动态系统。 从某种程度上来说,我们生活在一个动态系统中,窗外的天气从黎明到黄昏都在变化,我们体内发生的新陈代谢也是一个动态系统,因为随着时间的推移,成千上万的反应和分子被合成和降解。 更正式地说,如果我们定义一组变量,

    2024年02月01日
    浏览(34)
  • 【数学建模\MATLAB】掌握用Matlab求解微分方程问题

    d y d t = a y 2 cfrac{dy}{dt}=ay^2 d t d y ​ = a y 2 结果 d 3 y d t 3 = b y cfrac{d^3y}{dt^3}=by d t 3 d 3 y ​ = b y 结果 x 2 + y + ( x − 2 y ) y ′ = 0 x^2+y+(x-2y)yprime=0 x 2 + y + ( x − 2 y ) y ′ = 0 结果: d 2 y d t 2 = a y , 初 始 条 件 为 y ( 0 ) = 5 , y ′ ( 0 ) = 1 cfrac{d^2y}{dt^2}=ay, text初始条件为y(0)=5,yprime(

    2024年02月03日
    浏览(43)
  • 欧拉法与梯形法求解微分方程【含matlab源代码】

    本文介绍两种入门级求解微分方程的方法 —— 梯形法与欧拉法。 将上述方程组改写成matlab语言: 一、欧拉法 1.1 向前欧拉公式 1.2 向后欧拉公式     欧拉法求解源代码       欧拉法结果图   二、梯形法    梯形法求解源代码       梯形法结果图    感谢 Miracle 向公众

    2024年02月15日
    浏览(46)
  • 微分方程+传染病模型(指数传播、SI、SIS、SIR模型)+MATLAB求解

    本文为北海的数模课程学习笔记,课程出自微信公众号:数学建模BOOM。 求赞!求收藏!求关注! 微分方程结合传染病模型(如指数传播、SI、SIS、SIR模型)提供了一种用数学公式描述疾病传播动态的方法,有助于理解和预测疾病在人群中的传播路径和速度。 目录 指数传播模

    2024年02月04日
    浏览(57)
  • 0702可分类变量的微分方程-微分方程

    本节至第四节我们学习的都是一阶微分方程 ​ y ′ = f ( x , y ) y^{\\\'}=f(x,y) y ′ = f ( x , y ) (2-1) 一阶微分方程对称形式 p ( x , y ) d x + Q ( x , y ) d y = 0 ( 2 − 2 ) p(x,y)dx+Q(x,y)dy=0qquad (2-2) p ( x , y ) d x + Q ( x , y ) d y = 0 ( 2 − 2 ) 若以x为自变量,y为因变量,则 d y d x = − P ( x , y ) Q (

    2024年02月04日
    浏览(50)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包