在Python中无法得到RK4来求解轨道物体的位置

Question

我试图用一个更大质量的物体不运动的理想化方法，求出一个物体绕着一个更大的物体运行的位置。我试图用python中的四阶Runge-Kutta来求解笛卡尔坐标系中的位置。在

这是我的代码：

dt = .1
t = np.arange(0,10,dt)

vx = np.zeros(len(t))
vy = np.zeros(len(t))
x = np.zeros(len(t))
y = np.zeros(len(t))

vx[0] = 10 #initial x velocity
vy[0] = 10 #initial y velocity
x[0] = 10 #initial x position
y[0] = 0 #initial y position

M = 20

def fx(x,y,t): #x acceleration
     return -G*M*x/((x**2+y**2)**(3/2))

def fy(x,y,t): #y acceleration
     return -G*M*y/((x**2+y**2)**(3/2))

def rkx(x,y,t,dt): #runge-kutta for x

     kx1 = dt * fx(x,y,t)
     mx1 = dt * x
     kx2 = dt * fx(x + .5*kx1, y + .5*kx1, t + .5*dt)
     mx2 = dt * (x + kx1/2)
     kx3 = dt * fx(x + .5*kx2, y + .5*kx2, t + .5*dt)
     mx3 = dt * (x + kx2/2)
     kx4 = dt * fx(x + kx3, y + x3, t + dt)
     mx4 = dt * (x + kx3)

     return (kx1 + 2*kx2 + 2*kx3 + kx4)/6
     return (mx1 + 2*mx2 + 2*mx3 + mx4)/6

 def rky(x,y,t,dt): #runge-kutta for y

     ky1 = dt * fy(x,y,t)
     my1 = dt * y
     ky2 = dt * fy(x + .5*ky1, y + .5*ky1, t + .5*dt)
     my2 = dt * (y + ky1/2)
     ky3 = dt * fy(x + .5*ky2, y + .5*ky2, t + .5*dt)
     my3 = dt * (y + ky2/2)
     ky4 = dt * fy(x + ky3, y + ky3, t + dt)
     my4 = dt * (y + ky3)

     return (ky1 + 2*ky2 + 2*ky3 + ky4)/6
     return (my1 + 2*my2 + 2*my3 + my4)/6

for n in range(1,len(t)): #solve using RK4 functions
    vx[n] = vx[n-1] + fx(x[n-1],y[n-1],t[n-1])*dt
    vy[n] = vy[n-1] + fy(x[n-1],y[n-1],t[n-1])*dt
    x[n] = x[n-1] + vx[n-1]*dt
    y[n] = y[n-1] + vy[n-1]*dt

最初，无论我用哪种方式修改代码，我的for循环都会出错，要么是“类型为‘float’的对象没有len（）”（我不明白float python可能指的是什么），要么是“用序列设置数组元素”（我也不明白它的意思是什么序列）。我已经设法排除了错误，但我的结果是错误的。我得到了10的vx和vy数组，10的x数组。和一个从0到109的整数的y数组。到99。在

我怀疑fx（x，y，t）和fy（x，y，t）或者我把runge-kutta函数编码成fx和fy的方式有问题，因为我对其他函数使用了相同的runge-kutta代码，而且它工作得很好。在

我非常感谢任何帮助我弄清楚我的代码为什么不能工作的帮助。谢谢您。在

Answer 1

您没有在任何地方使用rkx，rky函数！ 在函数定义的末尾有两个return应该使用 return [(kx1 + 2*kx2 + 2*kx3 + kx4)/6, (mx1 + 2*mx2 + 2*mx3 + mx4)/6]（如@eapetcho所指出）。另外，你对Runge Kutta的实现我也不清楚。在

你有dv/dt，所以你求解v，然后相应地更新r。 在

for n in range(1,len(t)): #solve using RK4 functions
    vx[n] = vx[n-1] + rkx(vx[n-1],vy[n-1],t[n-1])*dt
    vy[n] = vy[n-1] + rky(vx[n-1],vy[n-1],t[n-1])*dt
    x[n] = x[n-1] + vx[n-1]*dt
    y[n] = y[n-1] + vy[n-1]*dt

这是我的代码版本 在

^{pr2}$

I suspect there are issues with fx(x,y,t) and fy(x,y,t)

在这种情况下，我刚刚检查了fx=3和{}的代码，得到了一个很好的轨迹。在

这是ry与rx图：

Answer 2

物理

牛顿定律给你一个二阶常微分方程u''=F(u)，其中u=[x,y]。使用v=[x',y']得到一阶系统

u' = v
v' = F(u)

它是四维的，必须用四维状态来求解。唯一可用的简化方法是使用开普勒定律，它允许将系统缩减到一个标量级的角度。但这不是这里的任务。在

欧拉法

您正确地实现了Euler方法来计算代码的最后一个循环中的值。它看起来不是物理的可能是因为Euler方法不断地增加轨道，因为它沿着切线移动到凸轨迹的外部。在您的实现中，可以看到G=100的这种向外螺旋。在

通过选择较小的步长，例如dt=0.001，可以有效地减少这种情况。在

你应该选择积分时间作为一个完整轨道的一个很好的部分，以得到一个体面的结果，与上述参数，你得到大约2个循环，这是好的。在

RK4实施

你犯了几个错误。不知怎么的你失去了速度，位置更新应该基于速度。在

那么您应该在fx(x + .5*kx1, y + .5*kx1, t + .5*dt)处停下来重新考虑您的方法，因为这与任何命名约定都不一致。一致、正确的变体是

^{pr2}$

这表明您不能解耦耦合系统的集成，因为您需要在x更新的同时进行y更新。此外，函数值是加速度，从而更新速度。位置更新使用当前状态的速度。因此，步骤应该从

 kx1 = dt * fx(x,y,t) # vx update
 mx1 = dt * vx        # x update
 ky1 = dt * fy(x,y,t) # vy update
 my1 = dt * vy        # y update

 kx2 = dt * fx(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
 mx2 = dt * (vx + 0.5*kx1/2)
 ky2 = dt * fy(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
 my2 = dt * (vy + 0.5*ky1/2)

等等

然而，正如您所见，这已经开始变得难以处理。将状态组合成一个向量，并对系统方程使用向量值函数

M, G = 20, 100
def orbitsys(u):
     x,y,vx,vy = u
     r = np.hypot(x,y)
     f = G*M/r**3
     return np.array([vx, vy, -f*x, -f*y]);

然后可以使用Euler或Runge-Kutta步骤的烹饪书实现

def Eulerstep(f,u,dt): return u+dt*f(u)

def RK4step(f,u,dt):
    k1 = dt*f(u)
    k2 = dt*f(u+0.5*k1)
    k3 = dt*f(u+0.5*k2)
    k4 = dt*f(u+k3)
    return u + (k1+2*k2+2*k3+k4)/6

并将它们组合成一个集成回路

def Eulerintegrate(f, y0, tspan):
    y = np.zeros([len(tspan),len(y0)])
    y[0,:]=y0
    for k in range(1, len(tspan)):
        y[k,:] = Eulerstep(f, y[k-1], tspan[k]-tspan[k-1])
    return y


def RK4integrate(f, y0, tspan):
    y = np.zeros([len(tspan),len(y0)])
    y[0,:]=y0
    for k in range(1, len(tspan)):
        y[k,:] = RK4step(f, y[k-1], tspan[k]-tspan[k-1])
    return y

用你给定的问题来调用它们

dt = .1
t = np.arange(0,10,dt)
y0 = np.array([10, 0.0, 10, 10])

sol_euler = Eulerintegrate(orbitsys, y0, t)
x,y,vx,vy = sol_euler.T
plt.plot(x,y)

sol_RK4 = RK4integrate(orbitsys, y0, t)
x,y,vx,vy = sol_RK4.T
plt.plot(x,y)