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!

Lederer

You will need to show the exact error you are getting. – Boatload

I have just edited the original post. Thanks! – Lederer

Does it matter that deriv_x(xinit,0) is not defined. – Occlude

unrelated to your "divide by zero" error, try: ipython notebook --pylab or ipython qtconsole --pylab to get a nicer interface. – Trinomial

I changed t to >t=linspace(0.01,24.0,100) and >xinit=array([0.001,5251]) and not I am getting another. – Lederer

How your equations look like? – Precontract

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?

Miniskirt

Thank you very much for the reply! I believe that is exactly the source of the error. However I modified t and xinit as seen on the second image on the original post and I still get an error. Have you got any suggestion on how to solve this issue? Sorry for the rather obvious questions, I am very new to Python and programming in general and my course has not been very well taught... Thanks again! – Lederer