解决ODE的python时出错

Question

我有一个大学项目,要求我们使用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])的第二个错误:

如果有人对如何改进代码有任何建议,我将非常感激.

非常感谢你!

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''没有这个问题-当分母不吹了x0 = 0,因为y0 = 20056,所以r0 = (x**2+y**2)**(1/2)是远离零.

结论:你分离的方法r''ODE成两个常微分方程的x''和y''不正确.

尝试搜索解决r''ODE 的不同方法.

暗示:

• 如果您的"状态"向量是z = [x, y, x', y']什么？
• 你可以写下一阶微分方程的z'来讲x,y, x'和y'？
• 你能打个电话解决integrate.odeint吗？