一个带有下标变量的数组

Question

我在使用主程序中的子例程的结果时遇到问题。我写了这段代码：

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

  call Scale_Factor(h,n,t,a)

  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

子例程求解一个微分方程，并为每个a生成t。我需要在主程序中调用子例程的结果，并在主程序中使用a(t)。我想将a(t)定义为一个数组，但是由于t是实数，所以我不能这样做。谁能帮我吗？

Answer 1

您有一个系统

R'(t)=f(a(t))
a'(t)=g(a(t))

是半耦合的。要将两个功能集成在一起，请使用主形式的耦合RK4方法

rk4step(R,a,h)
    k1a = h*g(a)
    k1R = h*f(a)
    k2a = h*g(a+k1a/2)
    k2R = h*f(a+k1a/2)
    k3a = h*g(a+k2a/2)
    k3R = h*f(a+k2a/2)
    k4a = h*g(a+k3a)
    k4R = h*f(a+k3a)
    a += (k1a+2*k2a+2*k3a+k4a)/6
    R += (k1R+2*k2R+2*k3R+k4R)/6
    return R,a

更好的方法是使用向量值状态和函数，以避免重复相似的步骤。