为何RK4算法在地球轨道模拟中无法更新位置与速度值？

以下是代码中导致变量无法更新的核心问题及修正方案：

关键错误点 #

• 整数除法导致增量完全失效：(1/6)是整数运算，结果为0，所有RK4的加权平均增量都变成0，变量自然无法更新。必须改为1.0/6（浮点数除法）。
• 微分方程定义完全错误：
  • 位置微分：dx/dt等于当前x方向速度vx，dy/dt等于当前y方向速度vy，但你写成了xVel*t和yVel*t（这是匀速位移公式，不是微分）。
  • 速度微分：引力加速度公式为a_x = -G*M*x / r³、a_y = -G*M*y / r³，不需要乘t，且dvy()函数里错误地将分子写成了x，应为y。
• RK4步骤中的变量混用与拼写错误：
  • k1x = dt(x,y,t); 缺少乘号，正确写法是dt * dx(x,y,t)。
  • 计算kx时最后一项误写为k，正确应为k4x。
  • 更新x的中间值时，错误用x的k值更新y（比如y + k1x/2），应对应使用y的k值（y + k1y/2），同理其他k2/k3项也要对应各自维度的增量。
  • 计算k1vy时错误调用dvx()，应调用dvy()。
• 变量拼写错误：intiial_vx应为initial_vx，避免初始速度值异常。

修正后的完整代码 #

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define dt 86400 // 1 day in seconds

const double G = 6.67e-11;
const double au = 1.496e11;
const double M = 1.99e30;

double vx(double x, double y);
double vy(double x, double y);
double dx(double x, double y);
double dy(double x, double y);
double dvx(double x, double y);
double dvy(double x, double y);

int main(){
    double initial_x = au;
    double initial_y = 0;
    double initial_vx = vx(initial_x, initial_y);
    double initial_vy = vy(initial_x, initial_y);
    double t = 0;

    double x = initial_x;
    double y = initial_y;
    double vx_val = initial_vx;
    double vy_val = initial_vy;

    for(int i=0;i<365;i++){
        // 计算位置x的RK4增量
        double k1x = dt * dx(x, y);  
        double k2x = dt * dx(x + k1x/2, y + k1y/2);
        double k3x = dt * dx(x + k2x/2, y + k2y/2);
        double k4x = dt * dx(x + k3x, y + k3y);
        double kx = (1.0/6) * (k1x + 2*k2x + 2*k3x + k4x);

        // 计算位置y的RK4增量
        double k1y = dt * dy(x, y);
        double k2y = dt * dy(x + k1x/2, y + k1y/2);
        double k3y = dt * dy(x + k2x/2, y + k2y/2);
        double k4y = dt * dy(x + k3x, y + k3y);
        double ky = (1.0/6) * (k1y + 2*k2y + 2*k3y + k4y);

        // 计算速度vx的RK4增量
        double k1vx = dt * dvx(x, y);
        double k2vx = dt * dvx(x + k1x/2, y + k1y/2);
        double k3vx = dt * dvx(x + k2x/2, y + k2y/2);
        double k4vx = dt * dvx(x + k3x, y + k3y);
        double kvx = (1.0/6) * (k1vx + 2*k2vx + 2*k3vx + k4vx);

        // 计算速度vy的RK4增量
        double k1vy = dt * dvy(x, y);
        double k2vy = dt * dvy(x + k1x/2, y + k1y/2);
        double k3vy = dt * dvy(x + k2x/2, y + k2y/2);
        double k4vy = dt * dvy(x + k3x, y + k3y);
        double kvy = (1.0/6) * (k1vy + 2*k2vy + 2*k3vy + k4vy);

        // 更新变量
        x += kx;
        y += ky;
        vx_val += kvx;
        vy_val += kvy;

        printf("%.3e\t%.3e\t%.3e\t%.3e\n", x, y, vx_val, vy_val);
        t += dt;
    }

    return 0;
}

// x方向速度（初始圆周运动速度）
double vx(double x, double y)
{
    double r = sqrt(x*x + y*y);
    // 初始位置y=0时，vx=0，避免atan(0/x)的问题
    if (y == 0) return 0.0;
    double theta = atan2(y, x);
    return sqrt(G*M / r) * sin(theta);
}

// y方向速度（初始圆周运动速度）
double vy(double x, double y) 
{
    double r = sqrt(x*x + y*y);
    double theta = atan2(y, x);
    return sqrt(G*M / r) * cos(theta);
}

// dx/dt = vx
double dx(double x, double y)
{
    return vx(x, y);
}

// dy/dt = vy
double dy(double x, double y)
{
    return vy(x, y);
}

// dvx/dt = 引力加速度x分量
double dvx(double x, double y)
{
    double r = sqrt(x*x + y*y);
    return (-G*M*x) / pow(r, 3);
}

// dvy/dt = 引力加速度y分量
double dvy(double x, double y)
{
    double r = sqrt(x*x + y*y);
    return (-G*M*y) / pow(r, 3);
}

额外说明：修正中把变量名vx、vy改为vx_val、vy_val，避免和函数名冲突；用atan2(y,x)替代atan(y/x)，解决x为0时的定义域问题。

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