6

我有一个大学项目，要求我们使用 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：

为 x 运行 odeint 时出错

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

第二种错误

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

非常感谢！

4

1 回答 1

5 
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''爆炸。的原始二阶 ODE不存在这个问题——因为, 所以远非零时，r''分母不会爆炸。x0 = 0y0 = 20056r0 = (x**2+y**2)**(1/2)

结论：您将r''ODE 分成两个 ODE 的x''方法y''是不正确的。

尝试寻找解决r''ODE 的不同方法。

暗示：

  • 如果您的“状态”向量是z = [x, y, x', y']怎么办？
  • z'你能写出关于x、y和 x'的一阶常微分方程y'吗？
  • 你可以通过一个电话解决它integrate.odeint吗？
于 2012-12-29T21:43:04.497 回答