无法让RK4解决Python中轨道物体的位置 [英] Cannot get RK4 to solve for position of orbiting body in Python

问题描述

我正在尝试使用更大的物体不移动的理想化方法来解决绕更大的物体运行的物体的位置.我正在尝试使用python中的四阶Runge-Kutta来求解笛卡尔坐标中的位置.

I am trying to solve for the position of a body orbiting a much more massive body, using the idealization that the much more massive body doesn't move. I am trying to solve for the position in cartesian coordinates using 4th order Runge-Kutta in python.

这是我的代码:

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到109.的整数的x数组和从0.到99的整数的y数组.

Originally, no matter which way I tweaked the code, I was getting an error on my for loop, either "object of type 'float' has no len()" (I didn't understand what float python could be referring to), or "setting an array element with a sequence" (I also didn't understand what sequence it meant). I've managed to get rid of the errors, but my results are just wrong. I get vx and vy arrays of 10s, an x array of integers from 10. to 109., and a y array of integers from 0. to 99.

我怀疑fx(x，y，t)和fy(x，y，t)或我将runge-kutta函数编码为fx和fy的方式存在问题，因为我已经使用过其他功能使用相同的runge-kutta代码，效果很好.

I suspect there are issues with fx(x,y,t) and fy(x,y,t) or with the way I have coded the runge-kutta functions to go with fx and fy, because I've used the same runge-kutta code for other functions and it works fine.

我非常感谢您提供帮助，以弄清我的代码为何无法正常工作.谢谢.

I greatly appreciate any help in figuring out why my code isn't working. Thank you.

推荐答案

物理

牛顿定律为您提供二阶ODE u''=F(u)和u=[x,y].使用v=[x',y']您可以获得一阶系统

Physics

The Newton law gives you a second order ODE u''=F(u) with u=[x,y]. Using v=[x',y'] you get the first order system

u' = v
v' = F(u)

是4维的，必须使用4维状态求解.唯一可用的减小方法是使用开普勒定律，该定律允许将系统减小到标量阶为角度的一个ODE.但这不是这里的任务.

which is 4-dimensional and has to be solved using a 4 dimensional state. The only reduction available is to use the Kepler laws which allow to reduce the system to a scalar order one ODE for the angle. But that is not the task here.

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

You correctly implemented the Euler method to calculate values in the last loop of your code. That it may look un-physical can be because the Euler method continuously increases the orbit, as it moves to the outside of convex trajectories following the tangent. In your implementation this outward spiral can be seen for G=100.

可以通过选择较小的步长(例如dt=0.001)来减小这种影响.

This can be reduced in effect by choosing a smaller step size, such as dt=0.001.

您应该选择积分时间作为整个轨道的良好部分，以获得可观的结果，使用上述参数，您将获得大约2个循环，这是很好的.

You should select the integration time to be a good part of a full orbit to get a presentable result, with above parameters you get about 2 loops, which is good.

您犯了几个错误.您以某种方式失去了速度，位置更新应基于速度.

You made several errors. Somehow you lost the velocities, the position updates should be based on the velocities.

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

Then you should have halted at fx(x + .5*kx1, y + .5*kx1, t + .5*dt) to reconsider your approach as that is inconsistent with any naming convention. The consistent, correct variant is

fx(x + .5*kx1, y + .5*ky1, t + .5*dt) 

这表明您无法解耦已耦合系统的集成，因为您需要y更新以及x更新.此外，函数值是加速度，因此更新了速度.位置更新使用当前状态的速度.因此，该步骤应从

which shows that you can not decouple the integration of a coupled system, as you need the y updates alongside the x updates. Further, the function values are the accelerations, thus update the velocities. The position updates use the velocities of the current state. Thus the step should start as

 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)
 ky2 = dt * fy(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
 my2 = dt * (vy + 0.5*ky1)

等

但是，正如您所看到的，这已经开始变得笨拙.将状态组装成向量，并对系统方程使用向量值函数

However, as you see, this already starts to become unwieldy. Assemble the state into a vector and use a vector valued function for the system equations

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步骤的菜谱实现

Then you can use a cook-book implementation of the Euler or Runge-Kutta step

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

并将它们组合成一个集成循环

and combine them into an integration loop

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

并根据您的给定问题调用它们

and invoke them with your given problem

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)

这篇关于无法让RK4解决Python中轨道物体的位置的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！