Geek Answers Handbook

Errors in solving ODE python

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

I managed to simulate it in 2D by making a second-order ODE into two first-order ODEs. However, I am stuck in a time limit because my code uses SI units, so it works in seconds, and the limits of Python Linspace do not even simulate one full orbit.

I tried to convert variables and constants to hours and kilometers, but now the code continues to give errors.

I followed this method:

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

And the code:

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 do not know how to copy / paste the error code from PyLab, so I took the PrintScreen 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 suggestions for improving the code, I will be very grateful.

Many thanks!

+6
source share
1 answer
 odeint(deriv_x, xinit, t)

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

 deriv_x(xinit, t)

causes a division by zero error, since x[0] = xinit[0] is 0, and deriv_x divides by x[0] .

It looks like you are trying to solve a second order ODE

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

where rhat is the unit vector in the radial direction.

It seems that you are 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'' explodes. The initial second-order ODE for r'' does not have this problem - the denominator does not explode when x0 = 0 , because y0 = 20056 , and therefore r0 = (x**2+y**2)**(1/2) far from zero.

Conclusion: your method of splitting an ODE r'' into two ODEs for x'' and y'' incorrect.

Try to find another way to solve r'' ODE.

Hint:

  • What if your state vector z = [x, y, x', y'] ?
  • Can you write first-order ODEs for z' in terms of x , y , x' and y' ?
  • Can you solve this problem with a single call to integrate.odeint ?
+5
source

More articles:


All Articles