基于四阶Runge-Kutta法求解Lane-Emden方程的C++代码错误排查

你的代码核心问题在于RK4步骤中计算中间状态时遗漏了步长的乘法，导致数值积分完全偏离正确结果，以下是具体错误分析和修正后的代码：

错误点分析 #

1. RK4中间状态计算遗漏步长：
   RK4方法中，中间点的状态值应该是当前状态 + 导数k * 步长dx * 0.5（对应k2、k3）或当前状态 + 导数k * 步长dx（对应k4）。你的代码中仅用了导数k * 0.5或导数k，漏掉了步长dxi的乘法，这会导致状态增量的量级完全错误。

2. 函数声明与定义参数名不匹配：
   dYdx的声明是double dYdx(double x, double Y);，但定义中第二个参数是Z，虽然编译不会报错，但会增加代码理解的混淆性，建议统一参数名。

3. 冗余参数neq：
   main函数中neq=1，但实际求解的是2个一阶方程，且RK4Step中直接循环处理2个状态量，该参数未被使用，可删除以避免误解。

修正后的代码 #

#include <iostream>
#include <cmath>
#include <iomanip>
#include <fstream>
#define Nmax 400

using namespace std;

void RK4Step(double xi, double *Y, double (*dYdx)(double, double), double (*dZdx)(double, double, double), double dxi);

double dYdx(double x, double Z);

double dZdx(double xi, double Y, double Z);

int main(){
    cout << setprecision(6);

    double Y[] = {1., 0.};
    double xi0 = 0.1;
    double xif = 20.;
    double dxi = (xif - xi0)/Nmax;

    ofstream fdata;
    fdata.open("lane_emden.dat");

    fdata << xi0 << " " << Y[0] << endl;
    cout << "xi: " << scientific << xi0 << "; theta: " << Y[0] << endl;

    for(int i=0; i<Nmax; i++){
        RK4Step(xi0, Y, dYdx, dZdx, dxi);
        xi0 += dxi;

        fdata << xi0 << " " << Y[0] << endl;
        cout << "xi: " << scientific << xi0 << "; theta: " << Y[0] << "; dtheta/dxi: " << Y[1] << endl;
    }
    fdata.close();

    return 0;
}

void RK4Step(double xi, double *Y, double (*dYdx)(double, double), double (*dZdx)(double, double, double), double dxi){
    double k1[2], k2[2], k3[2], k4[2];
    
    // 计算k1：当前点的导数
    k1[0] = dYdx(xi, Y[1]);
    k1[1] = dZdx(xi, Y[0], Y[1]);

    // 计算k2：xi + dxi/2处的导数，使用中间状态Y + k1*dxi/2
    double mid_Y_k2 = Y[0] + k1[0] * dxi * 0.5;
    double mid_Z_k2 = Y[1] + k1[1] * dxi * 0.5;
    k2[0] = dYdx(xi + dxi*0.5, mid_Z_k2);
    k2[1] = dZdx(xi + dxi*0.5, mid_Y_k2, mid_Z_k2);

    // 计算k3：xi + dxi/2处的导数，使用中间状态Y + k2*dxi/2
    double mid_Y_k3 = Y[0] + k2[0] * dxi * 0.5;
    double mid_Z_k3 = Y[1] + k2[1] * dxi * 0.5;
    k3[0] = dYdx(xi + dxi*0.5, mid_Z_k3);
    k3[1] = dZdx(xi + dxi*0.5, mid_Y_k3, mid_Z_k3);

    // 计算k4：xi + dxi处的导数，使用中间状态Y + k3*dxi
    double mid_Y_k4 = Y[0] + k3[0] * dxi;
    double mid_Z_k4 = Y[1] + k3[1] * dxi;
    k4[0] = dYdx(xi + dxi, mid_Z_k4);
    k4[1] = dZdx(xi + dxi, mid_Y_k4, mid_Z_k4);

    // 更新状态量
    for(int i=0; i<2; i++){
        Y[i] += (dxi/6.)*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
    }
}

double dYdx(double x, double Z){
    return Z;
}

double dZdx(double xi, double Y, double Z){
    return -1. - (2.*Z)/xi;
}

验证说明 #

n=0时Lane-Emden方程的解析解为$\theta(\xi) = 1 - \frac{\xi^2}{6}$，你可以将计算结果与该解析解对比，验证数值解的正确性。例如当$\xi=2$时，解析解$\theta=1 - 4/6=1/3≈0.333333$，修正后的代码计算结果应接近该值。

内容的提问来源于stack exchange，提问作者Andrea Barone