龙格-库塔RK4并不比Verlet好

Question

我只是在测试游戏中轨道动力学的几个整合方案。我把RK4与常数和自适应的步骤在这里http://www.physics.buffalo.edu/phy410-505/2011/topic2/app1/index.html

，我将它与简单的verlet集成(和euler，但它的性能非常非常差)进行了比较。不变步进的RK4似乎并不比verlet好。具有自适应步进的RK4更好，但也没有那么好。我想知道我是否做错了什么?或者在什么意义上说RK4比verlet优越得多?

认为Force在每个RK4步骤中评估4次，但在每个verlet步骤中仅评估1次。因此，为了获得相同的性能，我可以将verlet的time_step设置为小4倍。当时间步长减小4倍时，verlet比具有恒定步长的RK4更精确，几乎与具有自适应步长的RK4相当。

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

10T表示10个轨道周期，下面的数字48968、7920、48966是需要进行力评估的次数

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解算器家族非常擅长以相当高的精度解决几乎任何问题，并且在自适应方法的情况下，性能也很好。然而，它们不是辛的，这就意味着它们在问题中不保存能量。

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

你的问题是节能;任意转数后，行星体的总能量(动能+势能)应与初始条件下相同。对于Verlet积分器(至少作为时间窗平均值)，这将是正确的，但对于来自RK家族的积分器，它不会-随着时间的推移，RK解算器将建立一个不减少的误差，由于数值积分中的能量损失。

为了验证这一点，尝试在每个时间步骤中节省系统中的总能量，并绘制它(您可能需要进行超过10转才能注意到差异)。RK方法的能量稳定下降，而Verlet方法将在恒定能量周围振荡。

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

Answer 2

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

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

对于轨道力学，您还可以尝试RK7(8)自适应积分器，Adams-Bashforth多步或高斯杰克逊方法。下面是一篇展示其中一些方法的论文:http://drum.lib.umd.edu/bitstream/1903/2202/7/2004-berry-healy-jas.pdf

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

Answer 3

好的，最后，我使用了自适应龙格-库塔-费伯格(RKF45)。有趣的是，当我要求更高的精度(最优是1E-9)时，它会更快(需要更少的步骤)，因为在较低的精度(<1e-6)下，解决方案是不稳定的，并且由于丢弃的步骤(步骤太长且不精确)而浪费了许多迭代。当我要求更高的精度(1E-12)时，它需要更多的步骤，因为时间步骤更短。

对于圆形轨道，精度可以降低到(1e-5)，速度增益可达3倍，然而，当我需要模拟高度偏心(椭圆轨道)时，我宁愿保持安全的1E-9精度。

有代码，如果有人会解决类似的问题http://www.openprocessing.org/sketch/96977它还显示了模拟一个时间单位的轨迹长度需要多少次力评估