Leap Frog积分模拟地日轨道异常：仅输出竖直线问题排查

问题描述 #

尝试实现kick-drift-kick型Leap Frog积分模拟地球绕太阳的轨道系统，但运行后得到一条垂直直线而非预期的椭圆轨道，原代码如下：

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;


#define Grav    6.6726e-11             // m^3/kg/s^2
#define Msun    1.989e30               // kg
#define Mearth  5.973e24               // kg
#define Dse     1.496e11               // distance Sun-Earth in meters



void LeapFrog(int n) {

    // Initial Conditions
    double x0 = Dse;
    double y0 = 0;

    double vx0 = 0;
    double vy0 = sqrt(Grav*Msun/Dse);

    double px0 = Mearth*vx0;
    double py0 = Mearth*vy0;

    // Delta time
    double dt = 0.01;

    //we create a file for saving the coordinates
    std::ofstream outputFile("LF.dat");
    outputFile << x0 << " " << y0 << std::endl;


    for (int i = 0; i<n; i++) {

        double x = x0 + dt*px0/(2*Mearth);
        double y = y0 + dt*py0/(2*Mearth);

        outputFile << x << " " << y << std::endl;

        double px = px0 - dt* Grav* Msun * x/pow(sqrt(x*x + y*y), 3);
        double py = py0 - dt* Grav* Msun* y/pow(sqrt(x*x + y*y), 3);

        double x2 = x + dt*px/(2*Mearth);
        double y2 = y + dt*py/(2*Mearth);

        
        outputFile << x2 << " " << y2 << std::endl;

        //initialize the values for the following step
        x0 = x;
        y0 = y;

        px0 = px;
        py0 = py;

    }

    outputFile.close();
}



int main() {

    LeapFrog(48);

    return 0;

}

问题分析 #

1. 时间步长不合理：原代码dt=0.01（秒），模拟48步仅0.48秒，地球轨道周期约3.15×10^7秒，这么短时间内地球的位移微乎其微，大量密集输出的点看起来像直线。
2. Leap Frog实现逻辑错误：标准Kick-Drift-Kick流程应为：
   • 半步冲量：用当前位置计算力，更新动量到t+dt/2时刻
   • 全步漂移：用更新后的动量，将位置从t时刻推进到t+dt时刻
   • 半步冲量：用新位置计算力，将动量从t+dt/2更新到t+dt时刻
     原代码的状态迭代和输出逻辑混乱，错误地在一次循环中输出两个中间点，且未用最终位置更新下一次迭代的初始状态。
3. 状态更新错误：循环结束时将x0设为中间位置x而非最终位置x2，导致轨道计算完全偏离预期。

修正后的代码 #

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

#define Grav    6.6726e-11             // m^3/kg/s^2
#define Msun    1.989e30               // kg
#define Mearth  5.973e24               // kg
#define Dse     1.496e11               // distance Sun-Earth in meters
#define DAY_SEC 86400.0                // 一天的秒数

void LeapFrog(int n_steps) {
    // 初始条件：地球在x轴正方向，y=0；y方向初速度为圆周运动速度
    double x = Dse;
    double y = 0.0;
    double px = Mearth * 0.0;           // x方向初始动量为0
    double py = Mearth * sqrt(Grav * Msun / Dse); // y方向初始动量

    double dt = DAY_SEC;                // 以一天为时间步长
    ofstream outputFile("LF_correct.dat");
    outputFile << x << " " << y << endl; // 输出初始位置

    for (int i = 0; i < n_steps; ++i) {
        // 第一步：半步冲量（计算t时刻的力，更新动量到t+dt/2）
        double r = sqrt(x*x + y*y);
        double fx = -Grav * Msun * Mearth * x / (r*r*r); // 万有引力的x分量
        double fy = -Grav * Msun * Mearth * y / (r*r*r); // 万有引力的y分量
        double px_half = px + 0.5 * dt * fx;
        double py_half = py + 0.5 * dt * fy;

        // 第二步：全步漂移（用t+dt/2的动量更新位置到t+dt）
        double x_new = x + dt * px_half / Mearth;
        double y_new = y + dt * py_half / Mearth;

        // 第三步：半步冲量（计算t+dt时刻的力，更新动量到t+dt）
        double r_new = sqrt(x_new*x_new + y_new*y_new);
        double fx_new = -Grav * Msun * Mearth * x_new / (r_new*r_new*r_new);
        double fy_new = -Grav * Msun * Mearth * y_new / (r_new*r_new*r_new);
        double px_new = px_half + 0.5 * dt * fx_new;
        double py_new = py_half + 0.5 * dt * fy_new;

        // 输出新位置
        outputFile << x_new << " " << y_new << endl;

        // 更新状态为下一次迭代做准备
        x = x_new;
        y = y_new;
        px = px_new;
        py = py_new;
    }

    outputFile.close();
}

int main() {
    LeapFrog(365); // 模拟365天，接近一年的轨道
    return 0;
}

修正说明 #

• 时间步长调整：改用DAY_SEC（86400秒，即一天）作为步长，模拟365步对应一年，能清晰看到椭圆轨道。
• 标准Kick-Drift-Kick实现：严格按照半步冲量→全步漂移→半步冲量的流程迭代，确保积分精度。
• 状态更新与输出：每次迭代仅输出t+dt时刻的位置，用最终状态更新下一次迭代的初始值，保证轨道计算正确。

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