Leap Frog Integration 模拟地球绕太阳运行的问题

Question

我一直在尝试对地球绕太阳运行的系统实施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;

}

...我不断遇到同样的错误，我只得到一条垂直直线而不是轨道。

我已经尝试了几乎所有算法和不同的实现，问题还在于我需要进行踢-漂移-踢集成。我期待得到地球绕太阳的正确轨道。

Answer 1

好吧，你有一些问题 - 主要是身体上的问题。

您的代码将模拟 48 * 0.01 = 0.48 秒，但一年有 3.16e7 秒，所以它不会是一个轨道。增加步数并使时间步成为轨道时间（周长/速度）的实际分数。

px

和

py

指的是动量，所以在它们的变化中你需要的是总力，而不是总加速度：你忘了

Mearth

！更改为

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

py

也类似。

您应该在步骤结束时将

x

和

y

更新为

x2

和

y2

，而不是

x1

和

y1

。另外，您只需要每个时间步输出一次。

下图绘制了一个圆形轨道。如果你想要一个椭圆，稍微改变一下初始速度。

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

double Grav = 6.6726e-11;        // m^3/kg/s^2 (Universal gravitational constant)
double Msun = 1.989e30;          // kg
double Mearth = 5.973e24;        // kg
double Dse = 1.496e11;           // distance Sun-Earth in meters
int nsteps = 1000;


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 = ( 2 * 3.141592654 * Dse / vy0 ) / nsteps;      // fraction of a year

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

        double px = px0 - dt * Grav* Msun * Mearth * x / pow(sqrt(x*x + y*y), 3);        // Don't forget the earth!
        double py = py0 - dt * Grav* Msun * Mearth * 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 = x2;
        y0 = y2;

        px0 = px;
        py0 = py;

    }

    outputFile.close();
}



int main()
{
    LeapFrog( nsteps * 1 );                   // One orbit
}