6

私は、ODEとSciPyの常微分方程式関数を使用して火星への衛星アプローチをシミュレートするように依頼された大学のプロジェクトを持っています。

2次のODEを2つの1次のODEにすることで、2Dでシミュレートすることができます。ただし、コードがSI単位を使用しているため、数秒で実行され、Pythonのlinspace制限は、1つの完全な軌道をシミュレートすらしないため、時間制限にとらわれています。

変数と定数を時間とキロメートルに変換しようとしましたが、コードでエラーが発生し続けます。

私はこの方法に従いました：

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]）の2番目のエラー：

2番目のタイプのエラー

誰かがコードを改善する方法について何か提案があれば、私は非常に感謝します。

どうもありがとうございます！

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]ます。

2階常微分方程式を解こうとしているようです

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

ここで、rhatは半径方向の単位ベクトルです。

xとy座標を別々の2次常微分方程式に分離しているように見えます。

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

初期条件x0=0およびy0=20056の場合。

これは非常に問題があります。問題の中にはx0 = 0、x''が爆発するというものがあります。の元の2次常微分方程式にr''はこの問題はありません。分母は、がゼロから遠く離れているx0 = 0ためy0 = 20056に爆発しません。r0 = (x**2+y**2)**(1/2)

r''結論： ODEを2つのODEに分割する方法は正しくx''ありませんy''。

r''ODEを解くための別の方法を探してみてください。

ヒント：

  • 「状態」ベクトルがである場合はどうなりますz = [x, y, x', y']か？
  • 、、、 の1次常微分z'方程式xを書き留めてもらえますか？yx'y'
  • 1回の呼び出しで解決できますintegrate.odeintか？
于 2012-12-29T21:43:04.497 に答える