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!
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:
z = [x, y, x', y']
?z'
in terms of x
, y
,
x'
and y'
?integrate.odeint
?