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更优越？

我们认为Force的每个RK4步骤评估为4x,但每个verlet步骤只评估1x.所以为了获得相同的性能,我可以为verlet设置更小的time_step 4x.4x较小的时间步长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求解器会产生一个不会因为能量而减少的误差迷失在数值整合中.

要验证这一点,请尝试在每个时间步骤中保存系统中的总能量,然后绘制它(您可能需要做十次以上的旋转才能注意到差异).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

最后，如果您的力模型始终是一个简单的中心力（如本例所示），请查看开普勒方程。求解它精确、快速，并且可以跳到任意时间。