Python 3.x Runge Kutta 简单轨道

Question

我正处于创建使用 Runge-Kutta 方法绘制轨道的程序的早期阶段，并且希望以 2D 形式绘制轨道，但是，无论初始条件是什么，我都会得到一条直线。我看到了一个类似的问题，但它并没有解决我的问题。为什么会这样？

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 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
x=4596194 #initial conditions in m and m/s
y=4596194
vx=-6646
vy=6646
t=0
T=3600
bodyx = 444 #stationary body position metres
bodyy = 444
tarray=[]
xarray=[]
yarray=[]
vxarray=[]
vyarray=[]
while t<T:
    k1 = np.zeros(4)
    k2 = np.zeros(4)
    k3 = np.zeros(4)
    k4 = np.zeros(4)
    tarray.append(t)
    xarray.append(x)
    yarray.append(y)
    vxarray.append(vx)
    vyarray.append(vy)
    x = bodyx - x
    y = bodyy - y
    k1[0]=derX(vx)
    k1[1]=derY(vy)
    k1[2]=derVx(x,y)
    k1[3]=derVx(y,x)


    k2[0]=derX(vx+(step/2)*k1[2])
    k2[1]=derY(vy+(step/2)*k1[3])
    k2[2]=derVx(x+(step/2)*k1[0],y+(step/2)*k1[1])
    k2[3]=derVx(y+(step/2)*k1[1],x+(step/2)*k1[0])


    k3[0]=derX(vx+(step/2)*k2[2])
    k3[1]=derY(vy+(step/2)*k2[3])
    k3[2]=derVx(x+(step/2)*k2[0],y+(step/2)*k2[1])
    k3[3]=derVx(y+(step/2)*k2[1],x+(step/2)*k2[0])


    k4[0]=derX(vx+step*k3[2])
    k4[1]=derY(vy+step*k3[3])
    k4[2]=derVx(x+step*k3[0],y+step*k3[1])
    k4[3]=derVx(y+step*k3[1],vx+step*k3[0])

    t=t+step
    x=timestep(x,k1[0],k2[0],k3[0],k4[0])
    y=timestep(x,k1[1],k2[1],k3[1],k4[1])
    vx=timestep(x,k1[2],k2[2],k3[2],k4[2])
    vy=timestep(x,k1[3],k2[3],k3[3],k4[3])

mpl.plot(xarray, yarray)

Answer 1

k4[3] 的计算中有一个虚假的v。

timestep 的调用有 x 作为参数，它应该是 y, vx, vy。

另一个错误似乎是在差异计算中

x = bodyx - x
y = bodyy - y

你也改变了绝对位置。力的方向也会反转。

改成类似

diffx = x - bodyx
diffy = y - bodyy

并在力计算中使用这些相对位置。

---

为了比较，内置程序生成scipy.integrate.odeint 和数据

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

bodyx = 444 #stationary body position metres
bodyy = 444

def system(u,t):
    x,y,vx,vy = u
    x -= bodyx
    y -= bodyy
    f = -(G*M)/((x**2 + y**2)**(1.5))
    return [ vx, vy, f*x, f*y ]

x0=4596194 #initial conditions in m and m/s
y0=4596194
vx0=-6646
vy0=6646
u0 = [ x0, y0, vx0, vy0 ]

T= np.linspace(0,3600,36+1)

sol = odeint(system, u0, T)

mpl.plot(sol[:,0], sol[:,1]); mpl.show()

提供一个很好弯曲的弓，大约是完整轨道的 1/4。