RK4算法求解行星轨道微分方程所得图像异常问题排查

问题概述 #

我正在构建太阳系内行星间火箭飞行的2D模型，先以行星轨道模拟为基础。采用无量纲单位（令G*M_Sun/A.U**2 = 1），通过RK4方法求解如下微分方程组：

x' = vx
vx'= -x/(x**2+y**2)**1.5
y' = vy
vy'= -y/(x**2+y**2)**1.5

将方程表示为向量形式U' = f(U, t)（其中U = [x, y, vx, vy]，f = [vx, vy , -x/(x^2+y^2)^1.5, -y/(x^2+y^2)^1.5]），初始条件U0 = [x_0, 0, 0, vy_0]。

但自行实现的RK4代码运行后未得到预期的椭圆轨道，反而出现异常图像；而使用scipy.integrate.solve_ivp求解时能得到近乎正确的结果，需排查RK4实现中的错误。

原始错误代码 #

import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt

a1 = 1
eps1 = 0.0167 #Earth parameters
r_01 = (1-eps1)*(a1) 
v_01 = 1*((1+eps1)/(1-eps1))**0.5

def runge_kutta(x, t, func, dt):
    k1 = func(x, t)
    k2 = func(x+k1/2, t+dt/2)
    k3 = func(x+k2/2, t+dt/2)
    k4 = func(x+k3, t+dt)
    return dt*(k1+2*k2+2*k3+k4)/6

def f(U, t): #U = [x, y, vx, vy]
    res = np.zeros(4)
    res[0] = U[2]
    res[1] = U[3]
    res[2] = -U[0]/((U[0]**2+U[1]**2)**1.5)
    res[3] = -U[1]/((U[0]**2+U[1]**2)**1.5)
    return res

U_0 = np.zeros(4)
U_0[0]+=r_01
U_0[3]+=v_01

x=[]
y=[]
t=0
dt=0.1
while t<=2*np.pi:
    x.append(U_0[0])
    y.append(U_0[1])
    res = runge_kutta(U_0, t, f, dt)
    print(res)
    U_0 += res
    t += dt 

x = np.array(x)
y = np.array(y)
fig1, ax1 = plt.subplots(layout="tight")
ax1.scatter(x, y)
plt.show()

错误定位与修复 #

核心错误：RK4步长计算逻辑错误 #

RK4的标准公式中，k1/k2/k3/k4是状态量的增量，需由导数（func的返回值）乘以步长dt得到；而你的实现中：

1. 直接将导数赋值给k1/k2/k3/k4，未乘以dt，导致状态量更新时用导数直接累加，维度不匹配
2. 最终返回增量时又额外乘以dt，相当于对步长做了二次缩放，完全偏离了RK4的计算逻辑

修复后的RK4函数 #

def runge_kutta(x, t, func, dt):
    k1 = dt * func(x, t)
    k2 = dt * func(x + k1/2, t + dt/2)
    k3 = dt * func(x + k2/2, t + dt/2)
    k4 = dt * func(x + k3, t + dt)
    return (k1 + 2*k2 + 2*k3 + k4)/6

修复原理 #

• k1：当前状态下的导数乘以步长，得到初始增量
• k2：用初始增量的一半更新状态，计算该中间状态的导数并乘以步长
• k3：用k2的一半更新状态，计算对应导数并乘以步长
• k4：用k3更新状态，计算终点导数并乘以步长
• 最终增量为四个k值的加权平均，符合RK4的精度要求

验证结果 #

替换修复后的runge_kutta函数后，运行代码将得到与scipy.integrate.solve_ivp一致的椭圆轨道，符合开普勒运动的预期结果。

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