python - Runge–Kutta RK4 不比 Verlet 好？

Question

我只是在测试游戏中轨道动力学的几种集成方案。我在这里以持续和适应性的步骤参加了 RK4 http://www.physics.buffalo.edu/phy410-505/2011/topic2/app1/index.html

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

我们的想法是，每个 RK4 步对 Force 进行 4 次评估，但每个 verlet 步仅评估 1 次。因此，为了获得相同的性能，我可以将 verlet 的 time_step 设置为小 4 倍。 verlet 的时间步长减少了 4 倍，比具有恒定步长的 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 系列求解器非常擅长以相当高的精度解决几乎所有问题，如果是自适应方法，性能也很好。但是，它们不是 symplectic ，这很快就意味着他们在问题中不守恒。

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

你的问题是节能；经过任意次数的旋转后，行星体的总能量(动能 + 势能)应该与初始条件下的能量相同。这对于 Verlet 积分器来说是正确的(至少作为时间窗口平均值)，但它不会与来自 RK 系列的积分器一起使用 - 随着时间的推移，RK 求解器将建立一个不会减少的错误，由于能量迷失在数值积分中。

要验证这一点，请尝试在每个时间步保存系统中的总能量，并绘制它(您可能需要旋转十多次才能注意到差异)。 RK方法的能量稳定下降，而Verlet方法的能量围绕恒定能量振荡。

如果这正是您需要解决的问题，我还推荐开普勒方程，它以解析的方式解决了这个问题。即使对于具有许多行星、卫星等的相当复杂的系统，行星间的相互作用也微不足道，以至于您可以对每个天体及其旋转中心单独使用开普勒方程，而不会损失太多精度。但是，如果您正在编写游戏，您可能真的对其他一些问题感兴趣，而这只是一个学习示例 - 在这种情况下，请阅读各种求解器系列的属性，然后选择一个适合的你的问题。