龙格-库塔RK4不比弗雷特好吗？

Question

我只是在游戏中测试几个轨道动力学的积分方案。 我在这里采用了RK4的恒定自适应步骤 http://www.physics.buffalo.edu/phy410-505/2011/topic2/app1/index.html

我将其与简单的verlet积分（和euler积分）进行了比较，但是它的性能非常差。看来RK4的步幅不变并不比verlet好。带自适应步长的RK4更好，但不是很多。我想知道我是不是做错了什么？或者从哪个意义上说RK4比verlet优越得多？

这种想法是，每RK4步计算4倍的力，但每verlet步仅计算1倍的力。所以为了获得同样的性能，我可以为verlet设置4倍小的时间步长。用4x小的时间步长verlet比用常数步长的RK4更精确，几乎可以与用附加步长的RK4相媲美。

请参见图片： https://lh4.googleusercontent.com/-I4wWQYV6o4g/UW5pK93WPVI/AAAAAAAAA7I/PHSsp2nEjx0/s800/kepler.png

10T表示10个轨道周期，下面的数字48968792048966是所需的力评估数

python代码（使用pylab）如下：

from pylab import * 
import math

G_m1_plus_m2 = 4 * math.pi**2

ForceEvals = 0

def getForce(x,y):
    global ForceEvals
    ForceEvals += 1
    r = math.sqrt( x**2 + y**2 )
    A = - G_m1_plus_m2 / r**3
    return x*A,y*A

def equations(trv):
    x  = trv[0]; y  = trv[1]; vx = trv[2]; vy = trv[3];
    ax,ay = getForce(x,y)
    flow = array([ vx, vy, ax, ay ])
    return flow

def SimpleStep( x, dt, flow ):
    x += dt*flow(x)

def verletStep1( x, dt, flow ):
    ax,ay = getForce(x[0],x[1])
    vx   = x[2] + dt*ax; vy   = x[3] + dt*ay; 
    x[0]+= vx*dt;        x[1]+= vy*dt;
    x[2] = vx;        x[3] = vy;

def RK4_step(x, dt, flow):    # replaces x(t) by x(t + dt)
    k1 = dt * flow(x);     
    x_temp = x + k1 / 2;   k2 = dt * flow(x_temp)
    x_temp = x + k2 / 2;   k3 = dt * flow(x_temp)
    x_temp = x + k3    ;   k4 = dt * flow(x_temp)
    x += (k1 + 2*k2 + 2*k3 + k4) / 6

def RK4_adaptive_step(x, dt, flow, accuracy=1e-6):  # from Numerical Recipes
    SAFETY = 0.9; PGROW = -0.2; PSHRINK = -0.25;  ERRCON = 1.89E-4; TINY = 1.0E-30
    scale = flow(x)
    scale = abs(x) + abs(scale * dt) + TINY
    while True:
        x_half = x.copy();  RK4_step(x_half, dt/2, flow); RK4_step(x_half, dt/2, flow)
        x_full = x.copy();  RK4_step(x_full, dt  , flow)
        Delta = (x_half - x_full)
        error = max( abs(Delta[:] / scale[:]) ) / accuracy
        if error <= 1:
            break;
        dt_temp = SAFETY * dt * error**PSHRINK
        if dt >= 0:
            dt = max(dt_temp, 0.1 * dt)
        else:
            dt = min(dt_temp, 0.1 * dt)
        if abs(dt) == 0.0:
            raise OverflowError("step size underflow")
    if error > ERRCON:
        dt *= SAFETY * error**PGROW
    else:
        dt *= 5
    x[:] = x_half[:] + Delta[:] / 15
    return dt    

def integrate( trv0, dt, F, t_max, method='RK4', accuracy=1e-6 ):
    global ForceEvals
    ForceEvals = 0
    trv = trv0.copy()
    step = 0
    t = 0
    print "integrating with method: ",method," ... "
    while True:
        if method=='RK4adaptive':
            dt = RK4_adaptive_step(trv, dt, equations, accuracy)
        elif method=='RK4':
            RK4_step(trv, dt, equations)
        elif method=='Euler':
            SimpleStep(trv, dt, equations)
        elif method=='Verlet':
            verletStep1(trv, dt, equations)
        step += 1
        t+=dt
        F[:,step] = trv[:]
        if t > t_max:
            break
    print " step = ", step


# ============ MAIN PROGRAM BODY =========================

r_aphelion   = 1
eccentricity = 0.95
a = r_aphelion / (1 + eccentricity)
T = a**1.5
vy0 = math.sqrt(G_m1_plus_m2 * (2 / r_aphelion - 1 / a))
print " Semimajor axis a = ", a, " AU"
print " Period T = ", T, " yr"
print " v_y(0) = ", vy0, " AU/yr"
dt       = 0.0003
accuracy = 0.0001

#                 x        y     vx  vy
trv0 = array([ r_aphelion, 0,    0, vy0 ])             

def testMethod( trv0, dt, fT, n, method, style ):
    print " "
    F = zeros((4,n));
    integrate(trv0, dt, F, T*fT, method, accuracy);
    print "Periods ",fT," ForceEvals ",  ForceEvals
    plot(F[0],F[1], style ,label=method+" "+str(fT)+"T "+str(  ForceEvals ) ); 

testMethod( trv0, dt, 10, 20000  , 'RK4', '-' )
testMethod( trv0, dt, 10, 10000  , 'RK4adaptive', 'o-' )
testMethod( trv0, dt/4, 10, 100000, 'Verlet', '-' )
#testMethod( trv0, dt/160, 2, 1000000, 'Euler', '-' )

legend();
axis("equal")
savefig("kepler.png")
show();

Answer 1

好的，最后，我使用了自适应Runge-Kutta-Fehlberg（RKF45）。 有趣的是，当我要求更高的精度（最优值是1E-9）时，它会更快（需要更少的步骤），因为在较低的精度（<；1E-6）下，解是不稳定的，并且许多迭代被讨论的步骤浪费（步骤太长且不精确）。当我要求更高的精度（1E-12）时，它需要更多的步骤，因为时间步骤更短。

由于圆轨道的精度可以定为（1e-5）和高达3倍的速度增益，然而，当我需要模拟高度偏心（椭圆轨道）时，我宁愿保持安全的1e-9精度。

如果有人能解决类似的问题 http://www.openprocessing.org/sketch/96977 它还显示了模拟一个时间单位轨迹长度所需的力计算量

Answer 2

我不知道我是否要回答你的具体问题，但这里是我的想法。

你定义了一个非常简单的力模型。在这种情况下，保存一些步骤可能不会提高性能，因为在RK4中计算新步骤可能需要更长的时间。如果力模型比较复杂，采用自适应步长的RK4可以节省很多时间。从你的曲线图上我认为小红帽也偏离了正确的解，一个重复的椭圆。

对于轨道力学，你也可以尝试RK7（8）自适应积分器，Adams-Bashforth多步或高斯-杰克逊方法。这里有一篇文章展示了其中一些方法：http://drum.lib.umd.edu/bitstream/1903/2202/7/2004-berry-healy-jas.pdf

最后，如果你的力模型总是一个简单的中心力，就像在这个例子中，看看开普勒方程。解决它是精确的，快速的，你可以跳到一个任意的时间。

Answer 3

我知道这个问题现在已经很老了，但这与这些方法中的一种比另一种更“优越”无关，也与你对它们的编程无关——它们只擅长不同的事情。（所以不，这个答案并不是关于代码。甚至关于编程。更多的是关于数学，真的…）

Runge-Kutta解算器家族非常擅长以相当高的精度解决几乎所有问题，在自适应方法的情况下，性能也很好。然而，它们不是symplectic，这很快就意味着它们在问题中不节省能量。

另一方面，为了使解的振荡最小化，Verlet方法可能需要比RK方法小得多的步长，但是该方法是辛的。

你的问题是能量守恒；在任意数转之后，行星体的总能量（动能+势能）应该与初始条件相同。对于Verlet积分器（至少作为时间窗平均值），这是正确的，但对于RK族积分器，情况并非如此-随着时间的推移，RK解算器将建立一个错误，该错误不会因为数值积分中的能量损失而减少。

要验证这一点，请尝试在每个时间步节省系统中的总能量，并绘制出图（您可能需要旋转10圈以上才能注意到差异）。RK方法将有稳定的能量递减，而Verlet方法将给出围绕恒定能量的振荡。

如果这正是你需要解决的问题，我也推荐开普勒方程，它解析地解决这个问题。即使对于有许多行星、卫星等的相当复杂的系统，行星际的相互作用是如此的微不足道，以至于你可以对每个物体使用开普勒方程，并且它的旋转中心是独立的，而不会有太多的精度损失。然而，如果你正在编写一个游戏，你可能真的对其他问题感兴趣，这只是一个学习的例子-在这种情况下，阅读各种解算器族的属性，并选择一个适合你的问题。