计算科学 - 我的四阶 Runge-Kutta 方法实现中的错误在哪里？ - 吾爱随笔录

我正在尝试构建一个 Runge-Kutta 代码来整合一个简单谐振子的运动方程。但是，当我运行代码时，我只看到错误的一阶改进，因为我减小了步长。什么可能导致此问题？

这里我从二阶方程 x''+x=0 开始，引入一个变量 v=x'。然后我们有一个由两个一阶方程组成的系统：

x'=v
v'=-x.

以矩阵形式写成 (x',v')^T=M(x,v)^T（T 代表转置），其中 M=(0,1;-1,0)。

然后我们使用 (ts=timestep) 在 RK 公式 x(i+1)=x(i)+h/6*(A+2B+2C+2D) 中计算 A,B,C,D

A=M(x0,v0)
B=M[(x0,v0)+ts/2*A]
C=M[(x0,v0)+ts/2*B]
D=M[(x0,v0)+ts*C]

并应用上面写的 RK 公式。但这仅提供线性改进（即，将时间步长减半可将精度提高 2 倍）。

program rk4
implicit none

integer :: i,j,k
integer :: nos

double precision, allocatable :: y(:,:)
double precision :: x0,x1,x2,x3,x4,v0,v1,v2,v3,v4
double precision :: ts

allocate(y(2,10000))

x0=1.D0 !initial position
v0=0.D0 !initial velocity
ts=.025D0 !timestep
nos=1000 !number of steps

do i=1,nos

x1=v0
v1=-x0

x2=v0+ts/(2.D0)*v1 !a
v2=-x0-ts/(2.D0)*x1

x3=v0+ts/(2.D0)*v2 !b
v3=-x0-ts/(2.D0)*x2

x4=v0+ts/(2.D0)*v3 !c
v4=-x0-ts/(2.D0)*x3

x0=x0+ts/(6.D0)*(x1+(2.D0)*x2+(2.D0)*x3+x4) !d
v0=v0+ts/(6.D0)*(v1+(2.D0)*v2+(2.D0)*v3+v4)

y(1,i)=x0
y(2,i)=v0

end do

open(1,file='rk4.dat', status='unknown')
do i=1,nos
write(1,*) y(1,i), dble(i)*ts !use to plot a graph of the solution
end do
write(*,*) cos(6.3D0)-y(1,252) !error
end