## Sunday, June 28, 2015

### Orbits II (SciPy odeint)

Last night I was talking to a friend who asked me why I used python for science, and not something standard like FORTRAN or C. I said something like "they took all the most commonly used FORTRAN and C functions in the world, and wrapped them up in python for us" or something dreamy-eyed like that. While NumPy has been crafted with care for python users, a lot of SciPy is really standard, widely recognized packages, all just ready for us to invoke with a few lines of python.

Personally, I found it more satisfying to write and run and debug and use the Runge-Kutta algorithms (RK4, RK45) first, before I went ahead and just used the imensely powerful scipy.inegrate package odeint. There is a lot going on under the hood in odeint, and this is great. You can read this tutorial and the reference documentaiton.

Just for example, with dynamic step sizes (can save orders of magnitude of time) how do you get the position at fixed, predetermined times? With the explicit RK scripts in the previous post, I have to include tests and provisional one-step calculations. However, odeint takes care of all of this for us, either by stopping, backing-up, or interpolating between steps, whichever is numerically prudent. However, I still use explicit Runge-Kutta here and there, especially if SciPy is not available for some reason, or if I'm doing "weird" stuff.

If you have read the first Post (Orbits part 1), then this will look even better.

Anyway, here it is!

And this is the python, pretty nice don't you think?

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 from scipy.integrate import odeint as ODEint X0 = [r_aphelion, 0.0, 0.0, 0.0, 0.5*v_aphelion, 0.0] def func_ODE(X, t): f = np.zeros(6) # x, y, z, vx, vy, vz f[:3] = X[3:] # dx/dt = v one_over_r_cubed = ((X[:3]**2).sum())**-1.5 f[3:] = - Gm * X[:3] * one_over_r_cubed return f times = np.linspace(0.0, t_year, nstep+1) # stopping points to return data X, output_dict = ODEint(func_ODE, X0, timess, full_output=True) stepsizes = output_dict['hu'] 

How does it work?  We're solving 6 first-order differential equations based on a state vector:

X = [x, y, z, vx, vy, vz]

(dx)/(dt) = v_x,  (dy)/(dt) = v_y,  (dz)/(dt) = v_z,    is    f[:3] = X[3:],  and

(dv_x)/(dt) = a_x,  (dv_y)/(dt) = a_y,  (dv_z)/(dt) = z_z,    is    f[3:] = -Gm * X[:3] / r**3.

The algorithm will report back to you the calculated values of the state vector at the points in time you specify in times.  You can get a lot of information about what happened "under the hood" by setting full_output = True and examining the goodies returned in output_dict.