Errors while solving ODE's python

Question

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:

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

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

Thank you very much!

Answer 1

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:

• What if your "state" vector is z = [x, y, x', y']?
• Can you write down a first-order ODE for z' in terms of x, y, x' and y'?
• Can you solve it with one call to integrate.odeint?