四阶Runge-Kutta求解自治微分方程出现振荡不收敛问题求助

已发现的致命错误 #

• 语法逻辑错误（导致结果异常的核心原因）：状态更新段中，第三个分量的计算漏了dt与后续括号之间的乘号，原代码写为：

  f00[3*kk+2] = f00[3*kk+2]+1.0/6.0*dt(FieldLK1[3*kk+2]+2.0*FieldLK2[3*kk+2]+2.0*FieldLK3[3*kk+2]+FieldLK4[3*kk+2]);
  

  C语言中dt(...)会被解析为调用名为dt的函数，但dt是双精度浮点数变量，要么编译失败，要么运行时产生完全不符合预期的计算结果，正确写法应该补充乘号：

  f00[3*kk+2] = f00[3*kk+2]+1.0/6.0*dt*(FieldLK1[3*kk+2]+2.0*FieldLK2[3*kk+2]+2.0*FieldLK3[3*kk+2]+FieldLK4[3*kk+2]);
  

• 函数定义错误：主函数定义int main {缺少参数括号，不符合C语言标准，正确写法为int main() {。
• 迭代次数偏差：for (int i = 1; i < 100; i++)实际仅执行99次积分步，若预期执行100次，可修改为for (int i = 0; i < 100; i++)。
• 冗余初始化：每次积分步开头对所有临时变量执行归零操作完全多余，所有临时变量后续都会被新的赋值覆盖，该部分代码可直接删除，不影响计算逻辑。

修正后完整参考代码 #

int main() {
    int numberparticles = 6;
    double dt = 0.01;
    double f00[3*numberparticles];
    double FieldLK1[3*numberparticles],FieldLK2[3*numberparticles],FieldLK3[3*numberparticles],FieldLK4[3*numberparticles];
    double ArGk1[3*numberparticles],ArGk2[3*numberparticles],ArGk3[3*numberparticles],ArGk4[3*numberparticles];
    
    // 初始条件赋值
    for (int k = 0; k < numberparticles; k++) {
        f00[3*k] = 0.1;
        f00[3*k+1] = 0.2;
        f00[3*k+2] = 0.9;
    }
    
    // RK4迭代100次
    for (int i = 0; i < 100; i++) {
        // 计算k1
        RightSideEq(f00, ArGk1);
        for (int j = 0; j < numberparticles; j++) {
            FieldLK1[3*j] = ArGk1[3*j];
            FieldLK1[3*j+1] = ArGk1[3*j+1];
            FieldLK1[3*j+2] = ArGk1[3*j+2];
            // 计算k2对应的中间状态
            ArGk2[3*j] = f00[3*j] + 0.5 * FieldLK1[3*j] * dt;
            ArGk2[3*j+1] = f00[3*j+1] + 0.5 * FieldLK1[3*j+1] * dt;
            ArGk2[3*j+2] = f00[3*j+2] + 0.5 * FieldLK1[3*j+2] * dt;
        }
        
        // 计算k2
        RightSideEq(ArGk2, FieldLK2);
        for (int j = 0; j < numberparticles; j++) {
            // 计算k3对应的中间状态
            ArGk3[3*j] = f00[3*j] + 0.5 * FieldLK2[3*j] * dt;
            ArGk3[3*j+1] = f00[3*j+1] + 0.5 * FieldLK2[3*j+1] * dt;
            ArGk3[3*j+2] = f00[3*j+2] + 0.5 * FieldLK2[3*j+2] * dt;
        }
        
        // 计算k3
        RightSideEq(ArGk3, FieldLK3);
        for (int j = 0; j < numberparticles; j++) {
            // 计算k4对应的中间状态
            ArGk4[3*j] = f00[3*j] + FieldLK3[3*j] * dt;
            ArGk4[3*j+1] = f00[3*j+1] + FieldLK3[3*j+1] * dt;
            ArGk4[3*j+2] = f00[3*j+2] + FieldLK3[3*j+2] * dt;
        }
        
        // 计算k4
        RightSideEq(ArGk4, FieldLK4);
        
        // 更新状态变量
        for (int kk = 0; kk < numberparticles; kk++) {
            f00[3*kk] += 1.0/6.0 * dt * (FieldLK1[3*kk] + 2.0*FieldLK2[3*kk] + 2.0*FieldLK3[3*kk] + FieldLK4[3*kk]);
            f00[3*kk+1] += 1.0/6.0 * dt * (FieldLK1[3*kk+1] + 2.0*FieldLK2[3*kk+1] + 2.0*FieldLK3[3*kk+1] + FieldLK4[3*kk+1]);
            f00[3*kk+2] += 1.0/6.0 * dt * (FieldLK1[3*kk+2] + 2.0*FieldLK2[3*kk+2] + 2.0*FieldLK3[3*kk+2] + FieldLK4[3*kk+2]);
        }
    }
    return 0;
}

验证建议 #

• 修正代码后，先用简单自治方程验证RK4实现正确性，例如积分df/dt = -f，将数值解与解析解f(t) = f0 * exp(-t)对比，确认误差符合四阶方法的精度预期。
• 确认RightSideEq函数的数组长度要求、参数顺序与调用时一致，避免出现数组越界、参数传反的问题。
• 若积分对象为守恒系统，长时间积分下RK4本身不是辛算法，可能出现能量漂移，若对守恒性要求高可替换为辛RK算法，但本问题中你已经用CVODE验证过右侧函数正确，优先确认代码语法错误即可。

内容的提问来源于stack exchange，提问作者F.Mark