6

ゲームの軌道ダイナミクスのいくつかの統合スキームをテストしているところです。ここで一定の適応ステップでRK4を取りました http://www.physics.buffalo.edu/phy410-505/2011/topic2/app1/index.html

そして、それを単純なベルレット統合(およびオイラーですが、パフォーマンスが非常に悪い)と比較しました。一定のステップを備えたRK4がベルレよりも優れているとは思えません。アダプティブステップを備えたRK4の方が優れていますが、それほどではありません。私は何か間違ったことをしているのだろうか?それとも、どのような意味で RK4 がベルレよりもはるかに優れていると言われていますか?

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();
4

3 に答える 3

13

この質問は今ではかなり古いことを知っていますが、これは実際には、これらのメソッドのいずれかが他のメソッドよりも「優れている」ことや、それらのプログラミングとは何の関係もありません。(だからいいえ、この答えは実際にはコードに関するものではありません。プログラミングに関するものでさえありません。それは数学に関するものです、本当に...)

Runge-Kutta ファミリーのソルバーは、ほぼすべての問題を非常に優れた精度で攻撃するのに優れており、適応法の場合はパフォーマンスも優れています。ただし、それらはシンプレクティックではありません。つまり、問題でエネルギーを保存しないということです。

一方、Verlet 法では、解の振動を最小限に抑えるために RK 法よりもはるかに小さいステップ サイズが必要になる場合がありますが、この方法はシンプレクティックです。

あなたの問題はエネルギー保存です。任意の回数の回転の後、惑星体の総エネルギー (運動 + ポテンシャル) は、初期条件の場合と同じになるはずです。これは、Verlet 積分器では (少なくとも時間ウィンドウ平均として) 当てはまりますが、RK ファミリの積分器では当てはまりません。時間の経過とともに、RK ソルバーは、エネルギーのために減少しないエラーを蓄積します。数値積分で負けました。

これを確認するには、各タイム ステップでシステムの総エネルギーを保存してプロットします (違いに気付くには、10 回以上の回転が必要になる場合があります)。RK 法ではエネルギーが着実に減少しますが、Verlet 法では一定のエネルギーを中心に振動します。

これがまさにあなたが解決する必要がある問題である場合は、この問題を分析的に解決するケプラー方程式もお勧めします。多くの惑星、月などを含むかなり複雑なシステムでも、惑星間の相互作用は非常に重要ではないため、精度をあまり損なうことなく、各ボディとその回転中心に個別にケプラー方程式を使用できます。ただし、ゲームを作成している場合は、他の問題に本当に興味があるかもしれません。これは学習するための単なる例にすぎません。その場合は、さまざまなソルバー ファミリのプロパティを読んで、適切なものを選択してください。あなたの問題。

于 2013-08-14T01:43:57.107 に答える
4

あなたの特定の質問に答えるかどうかはわかりませんが、ここに私の考えがあります。

非常に単純な力モデルを定義しました。この場合、いくつかのステップを保存しても、RK4 での新しいステップの計算に時間がかかる可能性があるため、パフォーマンスが向上しない可能性があります。力モデルがより複雑な場合、適応ステップを備えた RK4 を使用すると、時間を大幅に節約できます。あなたのプロットから、ベルレも正しい解である繰り返し楕円から逸脱していると思います。

軌道力学については、RK7(8) アダプティブ インテグレーター、Adams-Bashforth マルチステップ、または Gauss Jackson メソッドを試すこともできます。これらの方法のいくつかを示す論文は次のとおりです

最後に、この例のように、力モデルが常に単純な中心力である場合は、ケプラーの方程式を見てください。それを解くことは正確で速く、任意の時間にジャンプできます。

于 2013-04-18T20:43:00.390 に答える
2

OK、最後に、適応ルンゲ クッタ フェルベルク (RKF45) を使用しました。興味深いことに、より高い精度 (最適は 1E-9 ) を要求すると、より高速です (必要なステップが少なくなります)。これは、精度が低い (<1e-6) 場合、解が不安定になり、破棄されたステップによって多くの反復が無駄になるためです (長くて不正確な手順)。さらに高い精度 ( 1E-12 ) を求めると、タイム ステップが短くなるため、より多くのステップが必要になります。

円軌道の場合、精度は ( 1e-5 ) まで下げることができ、最大 3 倍の速度が得られますが、非常に偏心した (楕円軌道) をシミュレートする必要がある間は、安全な 1E-9 精度を維持することを好みます。

誰かが同様の問題を解決する場合のコードがあり ます http://www.openprocessing.org/sketch/96977 また、軌跡の長さの1時間単位をシミュレートするために必要な力の評価の数も示しています

于 2013-04-22T13:08:53.960 に答える