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