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
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值