使用scipy solve_ivp求解一阶运动ODE仅得到少量点的问题排查

你遇到的问题由3处核心错误导致，逐一修正即可得到预期椭圆轨迹：

• 微分方程逻辑错误：motion函数内没有从状态变量Z中提取当前位置x、y，且加速度计算错误叠加了速度分量。万有引力加速度仅与位置相关，不需要乘速度vx/vy。
• solve_ivp参数传值错误：t_span参数仅接受「起始时间, 结束时间」的二元序列，你传入的全量时间数组需要赋值给t_eval参数，才能让求解器返回对应时间点的结果，这是你仅得到4个点的直接原因。
• 模拟参数冗余：设置的1000单位模拟时间过长，会产生大量冗余数据，也容易累积数值误差，可适当缩短时间范围。

修正后完整代码如下：

import math
import matplotlib.pyplot as plt 
import numpy as np 
import scipy.integrate

gim = 4*(math.pi**2)
x0 = 1 # x初始位置
y0 = 0 # y初始位置
vx0 = 0 # x方向初始速度
vy0 = 1.1 * 2 * math.pi # y方向初始速度
initial = [x0, y0, vx0, vy0] # 系统初始状态
# 缩短模拟时间到5个周期单位，步长0.01足够绘制光滑轨迹
time = np.arange(0, 5, 0.01)

def motion(t, Z): 
    x, y = Z[0], Z[1]
    dx = Z[2] # vx
    dy = Z[3] # vy
    r_cubed = (x**2 + y**2) ** (3/2)
    dvx = -gim * x / r_cubed
    dvy = -gim * y / r_cubed
    return [dx, dy, dvx, dvy]

# 修正参数：t_span传时间上下限，t_eval传要输出的时间点数组
sol = scipy.integrate.solve_ivp(motion, t_span=(time[0], time[-1]), y0=initial, method='RK45', t_eval=time)
plt.plot(sol.y[0], sol.y[1], label="Scipy RK45 solution")
plt.axis('equal') # 保证坐标轴比例一致，避免椭圆被拉伸
plt.legend()
plt.show()

补充说明：新增了plt.axis('equal')保证横纵坐标轴比例一致，避免绘制出的椭圆被拉伸变形。

内容的提问来源于stack exchange，提问作者Laura V.