C语言 龙格库塔评估不适用于 1 以外的步长

Question

我一直在将 4 个链接微分方程的欧拉方法实现转换为 4 阶 Runge Kutta 实现。我有理由确定我已经掌握了正确的一般方法，并且我已经了解如何应用 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。