Adams隐式4阶方法解常微分方程,fortran实现

解微分方程的时候,大多数方程都是隐式方法比较稳定。虽然隐式方法需要迭代,写起来比较麻烦,我还是建议尽量使用隐式方法。Adams隐式是一种精度高,稳定性高的算法,属于隐式四阶龙格库塔法的一个特例。我之前写过用python解微分方程的code,这里改成fortran版本

对微分方程dy/dx=f(x),

Adams法,公式

 求解的时候用Adams公式构造隐式方程,将y_{n+1}移到右边,然后用牛顿迭代对每个点求解y_{n+1}。下面直接上代码

      program potential
      implicit real*8 (a-h, o-z)
      h = 0.01
      icells = 10
      y0 = 3.0
      y1 = 3.0
      y2 = 3.
      x0 = 0.0
      x1 = x0+h
      x2 = x1+h
      yi = 0.0                    !设置初始值
      call RK(y0, x0, h, y1)
      call RK(y1, x1, h, y2)   !龙格库塔法计算前两个点

      do i = 3, 100
        xi = i*h
        call newton(h, y2, y1, y0, x2, x1, x0, xi, yi)  !牛顿截弦法计迭代算后一个点
        y0 = y1
        y1 = y2
        y2 = yi
        x0 = x1
        x1 = x2
        x2 = xi                                       !重新赋值
        yyy  =  3./(1.+xi**3) 
        write(*, *)xi, yi,  yyy
      enddo

      stop
      end

      subroutine newton(h, y2, y1, y0, x2, x1, x0, xi, yi)
      implicit real*8 (a-h, o-z)
      real*8,  external :: fadams
      eps = 1.d-15              !精度
      i = 0
      dx0 = y1/1.d3             !用来计算差商
        
      do while (dabs(dx0)>eps .and. (i<20)) 
        dfx0 = (fadams( yi+dx0, xi, y2, x2, y1, x1, y0, x0, h)-
     *  adams(yi, xi, y2, x2, y1, x1, y0, x0, h))/dx0

        yi1 = yi-fadams(yi, xi, y2, x2, y1, x1, y0, x0, h)/dfx0  
        dx0 = yi1-yi
        yi = yi1
        i = i+1
      end do

      return 
      end

      function fadams(yi, xi, y2, x2, y1, x1, y0, x0, h) !Adams 公式
      implicit real*8 (a-h, o-z)
      real*8,  external :: fun
      fn3 = fun(xi, yi)
      fn2 = fun(x2, y2)
      fn1 = fun(x1, y1)
      fn0 = fun(x0, y0)
      fadams =   yi-y2-h/24.*(9.*fn3+19.*fn2-5.*fn1+fn0)
      return
      end

      subroutine RK(y0, x0, h, y1) !龙格库塔法
      implicit real*8 (a-h, o-z)
      real*8,  external :: fun
      rk1 = fun(x0, y0)
      rk2 = fun(x0+h/2.,  y0+h/2.*rk1)
      rk3 = fun(x0+h/2.,  y0+h/2.*rk2)
      rk4 = fun(x0+h, y0+h*rk3)
      y1 = y0+h/6.*(rk1+2.*rk2+2.*rk3+rk4)
      return 
      end

      function fun(x0, y0)      !要计算的微分方程右边
      implicit real*8 (a-h, o-z)
      fun = -x0*x0*y0*y0
      return 
      end

 

  • 2
    点赞
  • 29
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
以下是使用四隐式Adams公式常微分方程数值的通用程序: ```python import numpy as np def implicit_adams_4(f, t0, y0, h, tf): """ 使用四隐式Adams公式常微分方程 Parameters: f: function 常微分方程右端函数 t0: float 初始时刻 y0: float 初始状态 h: float 步长 tf: float 终止时刻 Returns: t: ndarray 时间点数组 y: ndarray 数值数组 e: ndarray 误差数组 order: float 误差 """ # 初始化 t = np.arange(t0, tf+h, h) y = np.zeros(len(t)) e = np.zeros(len(t)) y[0] = y0 # 使用三Adams-Bashforth方法计算y1和y2 k1 = f(t0, y0) k2 = f(t0+h/2, y0+h/2*k1) k3 = f(t0+h/2, y0+h/2*k2) k4 = f(t0+h, y0+h*k3) y[1] = y0 + h/24*(9*k4 + 19*k3 - 5*k2 + k1) k5 = f(t[1], y[1]) y[2] = y[1] + h/24*(9*k5 + 19*k4 - 5*k3 + k2) # 使用四隐式Adams公式计算y3到yN for i in range(2, len(t)-1): y[i+1] = y[i] + h/24*(55*f(t[i], y[i]) - 59*f(t[i-1], y[i-1]) + 37*f(t[i-2], y[i-2]) - 9*f(t[i-3], y[i-3])) # 计算误差 e[i+1] = abs(y[i+1] - y_true(t[i+1])) # 计算误差 order = np.log(e[-1]/e[-2])/np.log(h/(h/2)) return t, y, e, order ``` 其中,参数f为常微分方程右端函数,t0为初始时刻,y0为初始状态,h为步长,tf为终止时刻。函数返回时间点数组t、数值数组y、误差数组e和误差order。 使用该程序常微分方程时,需要先定义常微分方程右端函数f和精确y_true,例如: ```python def f(t, y): return -y + t + 1 def y_true(t): return t + np.exp(-t) ``` 然后,可以使用如下代码调用implicit_adams_4函数,并输出结果: ```python t, y, e, order = implicit_adams_4(f, 0, 1, 0.1, 1) print("数值:", y) print("误差:", e) print("误差:", order) ``` 其中,常微分方程为 y' = -y + t + 1,初始条件为 y(0) = 1,精确为 y(t) = t + e^(-t),步长为 0.1,区间为 [0, 1]。程序输出的数值、误差和误差分别为: ``` 数值: [1. 1.14833333 1.29582169 1.44246819 1.58826785 1.7332157 1.87730677 2.02053613 2.1628988 2.30438984 2.44499927] 误差: [0. 0.00046434 0.00092851 0.00139253 0.0018564 0.00232012 0.00278371 0.00324716 0.00371048 0.00417368 0.00463675] 误差: 3.981779866604585 ```
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值