[Tutor] Problems using odeint

[Oscar Benjamin]

On Tue, Jan 14, 2014 at 06:40:03PM +0000, Grace Roberts wrote:
> Hi,
> I'm a python beginner, currently using scipy's 'odeint' to compute to a set of ODEs (obtained by splitting Newton's law of gravity ma=-Gm1m2r/r^3 in two ordinary differentials). The idea is that I'm modelling the orbit of a satellite as it approaches Mars.
> I have two problems:-For certain initial conditions the programme displays impossible orbits, showing the satellite making immediate sharp turns or jumping up and down in velocity. The exact same code can also end up displaying different graphs when run multiple times. For example when initial conditions are set at:  xx0=[1000.,10000.,1000.,10000.].-Often when run an error message appears saying: "Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information." You can see that as part of the odeint function I tried to stop the excess work by increasing the maxstep and I also tried to run with full_output =1. It's quite possible I haven't done this right. If I have, what else can I do to stop this message from appearing?
> I've attached a quick print-screen of the code.
> Any help is much appreciated.

Others have already said this but I'd like to repeat it: it would be a lot
easier for me to help with something like this if you put your code in the
email.

There are two problems with your code: one physical and one programmatic.

The physical problem is that your initial conditions are actually inside the
planet. Your code assumes that Mars is centred at the (0, 0) point. The radius
of Mars is ~3e6m but your initial condition is at (1000, 1000)m which is very
near the centre of Mars. The gravitational field at that point in Mars is
massive using the formula that you have shown. The large numbers will almost
certainly cause problems with the numerical accuracy of the integrator which
is probably why you see the warning messages.

In reality the gravitational field that close to the centre of Mars will be
very weak because of the shell theorem but your formula assumes a point mass
distribution so it doesn't take account of that:
http://en.wikipedia.org/wiki/Shell_theorem#Inside_a_shell
So the orbits you see there would be unphysical even if it wasn't for the
numerical accuracy problems.

If you change your initial condition to the following then your code will
produce sensible results:

xx0 = [
    4e6, # xpos - m - 1000 km above the surface of Mars along the x-axis
    0,   # xvel - m/s
    0,   # ypos - m
    2000,# yvel - m/s - initial velocity perpendicular to displacement
    ]

The second problem is programmatic. In your code you are treating your state
vector as [xpos, xvel, ypos, yvel] but at the end when you extract the
timeseries for each variable you have instead assumed a different ordering:
[xpos, ypos, xvel, yvel]. You should change the lines at the end to:

xx = soln[:, 0]
vx = soln[:, 1]
xy = soln[:, 2]
vy = soln[:, 3]

With those two changes I think you'll see a nice elliptic orbit.

I would suggest assigning the variable indices to names and then using those
in your code. So for example at the top of your file you could write:

# Variable indices
XPOS, YPOS, XVEL, YVEL = range(4)

Then your derivative function becomes more readable:

def f(x, t):
    r = ((x[XPOS]**2) + (x[YPOS]**2))**0.5
    dxdt = np.zeros_like(x)
    dxdt[XPOS] = x[XVEL]
    dxdt[YPOS] = x[YVEL]
    dxdt[XVEL] = G*Mm * x[XPOS]/r**3
    dxdt[YVEL] = G*Mm * x[YPOS]/r**3
    return dxdt

And at the end you can just write:

xx = soln[:, XPOS]
xy = soln[:, YPOS]
vx = soln[:, XVEL]
vy = soln[:, YVEL]

This way you'll be less likely to make mistakes with the indices.


Oscar