计算科学 - 一个程序中的两个 RK4 方法 - 吾爱随笔录

我想通过在 Fortran 中编码来使用 RK4 解决这个积分：

R = \int 1 / a (t) d t \to d R / d t = 1 / a (t) = f (t)

$R=∫1/a(t) dt → dR/dt=1/a(t) =f(t)$

初始点：t=0（或a=0.001）且R=0

我必须通过求解另一个微分方程得到 a(t)：

d a / d t = 1 / a + 1 / a^{2} = g (a)

$da/dt=1/a+1/a^2 =g(a)$

初始点：t=0 和 a=0.001

我写了这段代码来得到一个（t）：

    PROGRAM RK4
      implicit none
      real h,t
      integer n
      read*,h,n
      call Scale_Factor(h,n,t,a)
    END PROGRAM

    !---------------------------------------------

    SUBROUTINE Scale_Factor(h,n,t,a)
      implicit none
      real t,a,k1,k2,k3,k4,h,g
      integer i,n

      t=0
      a=0.001


    Do i=1,n

       k1=h*g(a)

       k2=h*g(a+k1/2.0)

       k3=h*g(a+k2/2.0)

       k4=h*g(a+k3)

       t=t+h

       a=a+(k1+2*k2+2*k3+k4)*(1/6.0)

       write(*,*)t,a

    END DO
    END SUBROUTINE

    !-------------------------
    FUNCTION g(a)
      implicit none
      real a,g
      g=sqrt((1.0/a)+(1.0/a**2)) 
    END FUNCTION

我还有另一个类似的程序来解决第一个积分。但是我需要使用这个程序产生的 a(t) 来求解积分，我不知道如何将它们组合在一个程序中。

我写的把它们结合起来是这样的：

    Program RK4
    implicit none

    real k1,k2,k3,k4,h,t,R
    integer i,n
    real a
    read*,n,h 


    t=0
    R=0

    Do i=1,n

      k1=h*(1/a(t))

      k2=h*(1/a(t+h/2.0))

      k3=h*(1/a(t+h/2.0))

      k4=h*(1/a(t+h))

      t=t+h

      R=R+(k1+2*k2+2*k3+k4)*(1/6.0)

      write(*,*)t,R

    End Do

    end program

    !-----------------------------------------

    SUBROUTINE Scale_Factor(h,n,t,a)
      implicit none
      real t,a,k1,k2,k3,k4,h,g
      integer i,n

      t=0
      a=0.001


    Do i=1,n

       k1=h*g(a)

       k2=h*g(a+k1/2.0)

       k3=h*g(a+k2/2.0)

       k4=h*g(a+k3)

       t=t+h

       a=a+(k1+2*k2+2*k3+k4)*(1/6.0)

       write(*,*)t,a

    END DO
    END SUBROUTINE

    !-------------------------
    FUNCTION g(a)
      implicit none
      real a,g
      g=sqrt((1.0/a)+(1.0/a**2)) 
    END FUNCTION

但我知道这是不正确的。