c - Runge Kutta 评估不适用于 1 以外的步长

Question

我一直在将 4 个链接微分方程的欧拉方法实现转换为四阶龙格库塔实现。我相当确信我的一般方法是正确的，并且我已经了解如何应用 RK4，但我可能已经有 6 年没有做过任何半严肃的数学了，所以我可能会错过某物。当我使用步长 1 时，我的 RK4 计算给出了合理的输出，但如果我使用任何其他步长，则很快会崩溃到零。我希望有一双新的眼睛能够很快发现我做错了什么。请不要发布完整的解决方案 - 我更喜欢指出我可能犯的错误 - 无论是代码还是对 RK4 的理解，因为我希望自己能够理解这一点。

这是我的欧拉实现。效果很好

// 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);
}

Answer 1

没有理由在 matsu_eulerStep 和 matsu_rkStep 的参数中使用 int 作为 step 的类型>。您应该像其他参数一样使用 double。