方程解的错误图，通过Python中的RK4算法求解

Question

我正在尝试制作一个火箭从太阳系中的一个行星飞往另一个行星的二维模型。在模拟太空飞行之前，我必须制作行星本身的模型。因此，对于每个行星，我必须通过 RK4 方法求解微分方程组（使用速度和时间的测量单位，使 G*M_Sun/A.U**2 = 1 以减少舍入误差）：

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

由于我不想写4次类似的表达式，所以我使用4-向量来制作方程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}

但是，这段代码（在我原来的基础上重写了一点，成为“独立的”，如果你的眼睛流血了，抱歉……）

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()

没有像我预期的那样画出椭圆，但是

使用 scipy.integrate.solve_ivp 函数求解方程可以得到（实际上）正确的椭圆轨道。所以，我认为这个错误是在我的 RK4 函数中，并且绝对是非常愚蠢的错误，因为该函数非常简单......但我找不到它。寻求帮助！

Answer 1

你的RK4功能确实有问题。

如果 X 的变化是 k = dt*gradient，那么您应该将 dt 包含在 k 的原始表达式中（并从最后的加权和中将其删除）。因此：

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

然后你得到