python - 求解 ODE 的 python 时出错

Question

我有一个大学项目，我们被要求使用 ODE 和 SciPy 的 odeint 函数模拟火星的卫星方法。

我设法通过将二阶 ODE 变成两个一阶 ODE 来在二维中对其进行仿真。然而，我被时间限制困住了，因为我的代码使用的是 SI 单位，因此在几秒钟内运行，而 Python 的 linspace 限制甚至不能模拟一个完整的轨道。

我尝试将变量和常量转换为小时和公里，但现在代码一直出错。

我遵循了这个方法:

http://bulldog2.redlands.edu/facultyfolder/deweerd/tutorials/Tutorial-ODEs.pdf

代码是:

import numpy

import scipy

from scipy.integrate import odeint

def deriv_x(x,t):
    return array([ x[1], -55.3E10/(x[0])**2 ]) #55.3E10 is the value for G*M in km and hours

xinit = array([0,5251]) # this is the velocity for an orbit of period 24 hours

t=linspace(0,24.0,100) 

x=odeint(deriv_x, xinit, t)

def deriv_y(y,t):
    return array([ y[1], -55.3E10/(y[0])**2 ])

yinit = array([20056,0]) # this is the radius for an orbit of period 24 hours

t=linspace(0,24.0,100) 

y=odeint(deriv_y, yinit, t)

我不知道如何从 PyLab 复制/粘贴错误代码，所以我截取了错误的 PrintScreen:

t=linspace(0.01,24.0,100) 和 xinit=array([0.001,5251]) 的第二个错误:

如果有人对如何改进代码有任何建议，我将不胜感激。

非常感谢!

Answer 1

odeint(deriv_x, xinit, t)

使用 xinit 作为它对 x 的初始猜测。 x 的这个值在评估 deriv_x 时使用。

deriv_x(xinit, t)

引发被零除错误，因为 x[0] = xinit[0] 等于 0，并且 deriv_x 除以 x[0].

---

看起来您正在尝试求解二阶 ODE

r'' = - C rhat
      ---------
        |r|**2

其中 rhat 是径向方向的单位矢量。

您似乎将 x 和 y 坐标分离为单独的二阶 ODES:

x'' = - C             y'' = - C
      -----    and          -----
       x**2                  y**2

初始条件 x0 = 0 和 y0 = 20056。

这是非常有问题的。问题之一是当 x0 = 0 时，x'' 爆炸了。 r'' 的原始二阶 ODE 没有这个问题——当 x0 = 0 时分母不会爆炸，因为 y0 = 20056，因此 r0 = (x**2+y**2)**(1/2) 远非零。

结论:您将 r'' ODE 分为 x'' 和 y'' 两个 ODE 的方法不正确。

尝试寻找一种不同的方法来求解 r'' ODE。

提示:

• 如果您的“状态”向量是 z = [x, y, x', y'] 会怎么样？
• 你能根据 x、y 写下 z' 的一阶 ODE，x' 和 y'?
• 你能一次调用 integrate.odeint 来解决这个问题吗？