使用Runge-Kutta法绘制2D天体轨道出现直线问题求助

看起来你在实现RK4求解二体问题时踩了几个关键的坑，导致轨道无法正常弯曲，最终变成直线。咱们一步步拆解代码里的问题：

1. 错误的坐标变换逻辑 #

你在循环里直接修改了x和y的原始值：

x = bodyx - x
y = bodyy - y

这会把当前位置直接替换成相对于中心天体的偏移量，但每次循环都会重复这个操作，相当于坐标被反复翻转，完全破坏了位置的累积更新逻辑。正确的做法是用临时变量存储相对位置，而非覆盖原位置变量：

rx = bodyx - x
ry = bodyy - y

之后所有计算加速度的地方都用rx和ry即可。

2. 复用加速度函数导致y方向加速度错误 #

你没有单独定义y方向的加速度函数，而是直接用derVx(y,x)来计算y方向加速度，这完全不符合万有引力的公式。万有引力在y方向的加速度应该是：

def derVy(x,y):
    return -(G*M*y)/((x**2 + y**2)**(3/2))

之前的写法会把x和y的位置搞反，导致y方向的加速度计算错误，轨道自然无法产生弯曲。

3. k4阶段的参数混淆错误 #

在计算k4[3]时，你写错了第二个参数：

k4[3]=derVx(y+step*k3[1],vx+step*k3[0])

这里第二个参数应该是x + step*k3[0]（位置x的更新值），而不是vx + step*k3[0]（速度x的更新值），属于典型的变量混淆低级错误。

4. 未完整更新所有状态变量 #

你的代码最后只写了x=timeste...，明显没有完成所有状态变量（x、y、vx、vy）的RK4更新。RK4需要对每个状态变量（位置x/y、速度vx/vy）分别应用更新公式，不能只更新x。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as mpl

def derX(vx):
    return vx

def derY(vy):
    return vy

def derVx(x, y):
    return -(G*M*x)/((x**2 + y**2)**(3/2))

def derVy(x, y):
    return -(G*M*y)/((x**2 + y**2)**(3/2))

def timestep(x, k1, k2, k3, k4):
    return x + (step/6)*(k1 + 2*k2 + 2*k3 + k4)

G = 6.67408E-11  # m^3/kg s^2
M = 5.972E24     # kg, mass of Earth
step = 100       # seconds

# 初始条件（位置m，速度m/s）
x = 4596194
y = 4596194
vx = -6646
vy = 6646

t = 0
T = 3600 * 24  # 延长模拟时间到24小时，便于观察完整轨道
bodyx = 0      # 将中心天体放在原点更直观
bodyy = 0

tarray = []
xarray = []
yarray = []
vxarray = []
vyarray = []

while t < T:
    tarray.append(t)
    xarray.append(x)
    yarray.append(y)
    vxarray.append(vx)
    vyarray.append(vy)
    
    # 计算相对中心天体的位置
    rx = bodyx - x
    ry = bodyy - y
    
    k1 = np.zeros(4)
    k2 = np.zeros(4)
    k3 = np.zeros(4)
    k4 = np.zeros(4)
    
    # k1阶段
    k1[0] = derX(vx)
    k1[1] = derY(vy)
    k1[2] = derVx(rx, ry)
    k1[3] = derVy(rx, ry)
    
    # k2阶段
    k2[0] = derX(vx + (step/2)*k1[2])
    k2[1] = derY(vy + (step/2)*k1[3])
    k2[2] = derVx(rx + (step/2)*k1[0], ry + (step/2)*k1[1])
    k2[3] = derVy(rx + (step/2)*k1[0], ry + (step/2)*k1[1])
    
    # k3阶段
    k3[0] = derX(vx + (step/2)*k2[2])
    k3[1] = derY(vy + (step/2)*k2[3])
    k3[2] = derVx(rx + (step/2)*k2[0], ry + (step/2)*k2[1])
    k3[3] = derVy(rx + (step/2)*k2[0], ry + (step/2)*k2[1])
    
    # k4阶段
    k4[0] = derX(vx + step*k3[2])
    k4[1] = derY(vy + step*k3[3])
    k4[2] = derVx(rx + step*k3[0], ry + step*k3[1])
    k4[3] = derVy(rx + step*k3[0], ry + step*k3[1])
    
    # 更新所有状态变量
    x = timestep(x, k1[0], k2[0], k3[0], k4[0])
    y = timestep(y, k1[1], k2[1], k3[1], k4[1])
    vx = timestep(vx, k1[2], k2[2], k3[2], k4[2])
    vy = timestep(vy, k1[3], k2[3], k3[3], k4[3])
    
    t += step

# 绘制轨道
mpl.figure(figsize=(8,8))
mpl.plot(xarray, yarray, label='Orbit')
mpl.scatter(bodyx, bodyy, color='red', label='Earth')
mpl.xlabel('X Position (m)')
mpl.ylabel('Y Position (m)')
mpl.legend()
mpl.axis('equal')
mpl.show()

额外说明 #

• 我把中心天体位置改成了原点（bodyx=0, bodyy=0），轨道形态会更直观；
• 延长了模拟时间到24小时（T=3600*24），原1小时的模拟时间可能看不到明显的轨道弯曲；
• 确保所有状态变量都通过RK4公式正确更新，运行后就能得到正常的椭圆轨道了。

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