python - 解决ODE的python时出错

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

通过将一个二阶ODE变成两个一阶ODE，我设法在2D中对其进行了仿真。但是我受时间限制，因为我的代码使用的是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])的第二个错误:

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

非常感谢你!

最佳答案

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']怎么办？

您可以按照z'，x，y和x'？

您可以通过一次调用y'来解决它吗？

关于python - 解决ODE的python时出错，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/14085412/