我正在做一个太阳系模拟。我正在使用 Ruth 的三阶辛积分器来避免能量漂移问题（我在 RK4 中遇到过），但是行星很快离开轨道，并且能量绝不是守恒的（就像 RK4 一样）。

我将其应用于 N 体问题，如下所示：

(KE=1/2mv^2)

我已经在 Fortran 2008 中实现了这一点，其中 x、a、v、p 和 m 都是长度为 30 的向量，它们包含 x、y、z 位置、x、y、z 加速度、x、y、z 速度， x,y,z 动量和 m,m,m 分别代表太阳系中 10 个独立的天体（行星 + 太阳 + 冥王星）。

每个物体上的加速度计算为每个物体上每个物体的 x,y,z 的 a=GM/(r^2) 之和。

这是代码的集成部分：

!----------Looping Through Time-----------
do while(t<365.250000d0) ! Length of simulation in days
    !----------Calculating Values-----------
    call calc_acc(masses,x,a)
    p1=p+(7.0d0/24.0d0)*h*m*a
    x1=x+(2.0d0/3.0d0)*h*p1/m
    call calc_acc(masses,x1,a)
    p2=p1+(3.0d0/4.0d0)*h*m*a
    x2=x1-(2.0d0/3.0d0)*h*p2/m
    call calc_acc(masses,x2,a)
    p=p2-(1.0d0/24.0d0)*h*m*a
    x=x2+h*p/m
    v=p/m
    t=t+h
    !----------Saving Values-----------
    do bodnum=1,10,1
        write((100+bodnum),*) t, x((1+3*(bodnum-1)):(3+3*(bodnum-1))), v((1+3*(bodnum-1)):(3+3*(bodnum-1)))
        write((200+bodnum),*) x((1+3*(bodnum-1))), x((2+3*(bodnum-1))), x((3+3*(bodnum-1)))
    end do 
end do 

完整的程序可以在这里找到。

请告诉我我做错了什么。