ルンゲクッタ4次粒子の移流コードサンプル

Question

は私がやっているシミュレーションにRK4統合を実装しようとしています。機能は以下のthisサイトで、12ページの式に基づいて、3次元で力場を経由統合するRK4を使用しての私の最高の試みです。私のコードで

、粒子クラスは、本質的に、速度と位置のリストを格納し、位置（力は速度とは無関係である）所定の力を計算することができます。また、私は私の関数が長いと引き出されていることを知っている、とループのために使用して低減することができますが、私は（現在は）私がリンクされた紙に使用される構造と一致します。

私はこの方法を使用して粒子をシミュレートしようとすると、エラーが私はリープフロッグ積分法を使用する場合よりも大幅に悪化しています。したがって、RK4の実装に何か問題があると思います。 RK4が連動微分方程式を解くために使用されているときに、RK4がどのように動作するのか誤解しているかどうか教えてください。

// 4th Order Runge-Kutta 
void Update(Particle * p, double dt) { 

    double * v = p->getVel(); 
    double * pos = p->getPos(); 

    double initPos[3] = {pos[0], pos[1], pos[2]}; 
    double initVel[3] = {v[0], v[1], v[2]}; 
    double mass = 0.01; 

    double k[4][3]; // related to dv 
    double l[4][3]; // related to dr 

    p->findForce(); 

    k[0][0] = dt*p->force[0]/mass; 
    k[0][1] = dt*p->force[1]/mass; 
    k[0][2] = dt*p->force[2]/mass; 

    l[0][0] = dt*v[0]; 
    l[0][1] = dt*v[1]; 
    l[0][2] = dt*v[2]; 

    // Set position to midpoint, using l[0] 
    pos[0] = initPos[0] + l[0][0]/2; 
    pos[1] = initPos[1] + l[0][1]/2; 
    pos[2] = initPos[2] + l[0][2]/2; 

    p->findForce(); 

    k[1][0] = dt*p->force[0]/mass; 
    k[1][1] = dt*p->force[1]/mass; 
    k[1][2] = dt*p->force[2]/mass; 

    l[1][0] = dt*(v[0]+k[0][0]/2); 
    l[1][1] = dt*(v[1]+k[0][1]/2); 
    l[1][2] = dt*(v[2]+k[0][2]/2); 

    // Set position to midpoint, using l[1] 
    pos[0] = initPos[0] + l[1][0]/2; 
    pos[1] = initPos[1] + l[1][1]/2; 
    pos[2] = initPos[2] + l[1][2]/2; 

    p->findForce(); 

    k[2][0] = dt*p->force[0]/mass; 
    k[2][1] = dt*p->force[1]/mass; 
    k[2][2] = dt*p->force[2]/mass; 

    l[2][0] = dt*(v[0]+k[1][0]/2); 
    l[2][1] = dt*(v[1]+k[1][1]/2); 
    l[2][2] = dt*(v[2]+k[1][2]/2); 

    // Set position to endpoint, using l[2] 
    pos[0] = initPos[0] + l[2][0]; 
    pos[1] = initPos[1] + l[2][1]; 
    pos[2] = initPos[2] + l[2][2]; 

    p->findForce(); 

    k[3][0] = dt*p->force[0]/mass; 
    k[3][1] = dt*p->force[1]/mass; 
    k[3][2] = dt*p->force[2]/mass; 

    l[3][0] = dt*(v[0]+k[2][0]); 
    l[3][1] = dt*(v[1]+k[2][1]); 
    l[3][2] = dt*(v[2]+k[2][2]); 

    // Finalize pos and v 
    pos[0] = initPos[0] + (l[0][0] + 2*l[1][0] + 2*l[2][0] + l[3][0])/6; 
    pos[1] = initPos[1] + (l[0][1] + 2*l[1][1] + 2*l[2][1] + l[3][1])/6; 
    pos[2] = initPos[2] + (l[0][2] + 2*l[1][2] + 2*l[2][2] + l[3][2])/6; 

    v[0] = initVel[0] + (k[0][0] + 2*k[1][0] + 2*k[2][0] + k[3][0])/6; 
    v[1] = initVel[1] + (k[0][1] + 2*k[1][1] + 2*k[2][1] + k[3][1])/6; 
    v[2] = initVel[2] + (k[0][2] + 2*k[1][2] + 2*k[2][2] + k[3][2])/6; 
} 

Answer 1

一度に一つの粒子を統合することができない、これは、粒子のアンサンブルため1方法をもたらし、したがって、注文4の方法に適しているステップサイズで大きなドリフトを提示します。

すべてのパーティクルを一度に統合する必要があります。つまり、すべてのパーティクルについてステージ0を計算し、すべてのパーティクルを使用してステージ1に状態1を設定し、力と速度を計算します。つまり、ステージ1ではk状態1から、一度に全ての粒子次に一度に全ての粒子のものkベクトルを計算し、ステージ2用に状態2を計算等