pythonscipydifferential-equationsnumerical-integrationorbital-mechanics

Errors while solving ODE's python


I have a university project in which we are asked to simulate a satellite approach to Mars using ODE's and SciPy's odeint function.

I manage to simulate it in 2D by making a second-order ODE into two first-order ODE's. However I am stuck in the time limitation because my code is using SI units therefore running in seconds and Python's linspace limits does not even simulate one complete orbit.

I tried converting the variables and constants to hours and kilometers but now the code keeps giving errors.

I followed this method:

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

And the code is:

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)

I don't know how to copy/paste the error code from PyLab so I took a PrintScreen of the error:

Error when running odeint for x

Second error with t=linspace(0.01,24.0,100) and xinit=array([0.001,5251]):

Second type of error

If anyone has any suggestions on how to improve the code I will be very grateful.

Thank you very much!


Solution

  • odeint(deriv_x, xinit, t)
    

    uses xinit as its initial guess for x. This value for x is used when evaluating deriv_x.

    deriv_x(xinit, t)
    

    raises a divide-by-zero error since x[0] = xinit[0] equals 0, and deriv_x divides by x[0].


    It looks like you are trying to solve the second-order ODE

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

    where rhat is the unit vector in the radial direction.

    You appear to be separating the x and y coordinates into separate second-order ODES:

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

    with initial conditions x0 = 0 and y0 = 20056.

    This is very problematic. Among the problems is that when x0 = 0, x'' blows up. The original second-order ODE for r'' does not have this problem -- the denominator does not blow up when x0 = 0 because y0 = 20056, and so r0 = (x**2+y**2)**(1/2) is far from zero.

    Conclusion: Your method of separating the r'' ODE into two ODEs for x'' and y'' is incorrect.

    Try searching for a different way to solve the r'' ODE.

    Hint: