0

私は、4 つのリンクされた微分方程式のオイラー法実装を 4 次ルンゲ クッタ実装に変換してきました。私は一般的なアプローチが正しいと合理的に確信しており、RK4の適用方法を理解していますが、おそらく6年間、中途半端な深刻な数学を行っていないため、何かが欠けている可能性があります. 私の RK4 計算は、ステップ サイズ 1 を使用すると適切な出力が得られますが、他のステップ サイズを使用するとすぐにゼロになります。新鮮な目が私が間違ったことをすぐに見つけられることを願っています. 完全な解決策を投稿しないでください-私が作った可能性のあるエラーへのポインタを好むでしょう-それがコードであろうとRK4の理解であろうと、私はこれを自分で理解できるようにしたいからです。

これが私の Euler の実装です。うまく機能します

// these defs are used in both Euler and RK4 versions
#define POSPART(X)  (X > 0.0 ? X : 0.0)
#define NEGPART(X)  (X < 0.0 ? X : 0.0)
#define STEP_DIVISION 10.0

double matsu_calc_nextVal(double in, double t1, double t2,
                              double c, double b, double g,
                              matsu_state *state)
    {
        double step = 1.0 / STEP_DIVISION;
        double posX1, posX2, posIn, negIn, dx1, dx2, dv1, dv2;
        double t1recip = 1.0/ t1;
        double t2recip = 1.0/ t2;


        for(int i=0; i<STEP_DIVISION; i++){
            posX1 = POSPART(state->x1);   
            posX2 = POSPART(state->x2);
            posIn = POSPART(in);
            negIn = NEGPART(in);
            dx1 = (c - state->x1 - (b*(state->v1)) - (posX2*g) - posIn) * t1recip;
            dx2 = (c - state->x2 - (b*(state->v2)) - (posX1*g) - negIn) * t1recip;
            dv1 = (posX1 - state->v1) * t2recip;
            dv2 = (posX2 - state->v2) * t2recip;
            state->x1 += dx1*step;         state->x2 += dx2*step;
            state->v1 += dv1*step;         state->v2 += dv2*step;
        }
        return POSPART(state->x1) - POSPART(state->x2);
    }

これは私のRK4実装です

// calculates derivative
matsu_state matsu_calc_deriv(double in, double t1recip, double t2recip,
                          double c, double b, double g,
                          matsu_state state)
{
    double posX1 = POSPART(state.x1); 
    double posX2 = POSPART(state.x2);
    double posIn = POSPART(in);
    double negIn = NEGPART(in);

    matsu_state deriv;
    deriv.x1 = (c - state.x1 - (b*(state.v1)) - (g*posX2) - posIn) * t1recip;
    deriv.x2 = (c - state.x2 - (b*(state.v2)) - (g*posX1) - negIn) * t1recip;
    deriv.v1 = (posX1 - state.v1) * t2recip;
    deriv.v2 = (posX2 - state.v2) * t2recip;

    return deriv;
}

//helper function for moving the state forward by derivative*step
matsu_state matsu_eulerStep(matsu_state init, matsu_state deriv, int step)
{
    matsu_state newVal;
    newVal.x1 = init.x1 + (deriv.x1*step);
    newVal.x2 = init.x2 + (deriv.x2*step);
    newVal.v1 = init.v1 + (deriv.v1*step);
    newVal.v2 = init.v2 + (deriv.v2*step);
    return newVal;
}

// single RK4 step
void matsu_rkStep (double in, double t1recip, double t2recip,
                          double c, double b, double g,
                          matsu_state *state, int step){
  matsu_state k1, k2, k3, k4;

  k1=matsu_calc_deriv(in, t1recip, t2recip, c, b, g, 
                                        *state);
  k2=matsu_calc_deriv(in, t1recip, t2recip, c, b, g, 
                                        matsu_eulerStep(*state, k1,step*0.5));
  k3=matsu_calc_deriv(in, t1recip, t2recip, c, b, g, 
                                        matsu_eulerStep(*state, k2,step*0.5));
  k4=matsu_calc_deriv(in, t1recip, t2recip, c, b, g, 
                                        matsu_eulerStep(*state, k3,step));

  state->x1 = state->x1 + (k1.x1 + 2*(k2.x1+k3.x1) +k4.x1)*(1.0/6.0)*step; 
  state->x2 = state->x2 + (k1.x2 + 2*(k2.x2+k3.x2) +k4.x2)*(1.0/6.0)*step; 
  state->v1 = state->v1 + (k1.v1 + 2*(k2.v1+k3.v1) +k4.v1)*(1.0/6.0)*step; 
  state->v2 = state->v2 + (k1.v2 + 2*(k2.v2+k3.v2) +k4.v2)*(1.0/6.0)*step; 

}

// wrapper to loop Rk4 step enough times to move forward by 1
double matsu_calc_nextVal_RK(double in, double t1, double t2,
                          double c, double b, double g,
                          matsu_state *state)
{
    double step = 1.0 / STEP_DIVISION;
    double t1recip = 1.0/t1;
    double t2recip = 1.0/t2;

    for(int i=0; i<STEP_DIVISION; i++){
        matsu_rkStep(in, t1recip, t2recip, c, b, g, state, step);
    }
    return POSPART(state->x1) - POSPART(state->x2);
}
4

1 に答える 1

1

およびの引数で のint型として使用する必要がある理由はありません。他のパラメーターと同じように使用する必要があります。stepmatsu_eulerStepmatsu_rkStepdouble

于 2016-03-17T01:54:39.497 に答える