ode_solver_comparison.py

import odespy
import numpy as np
dt = np.float64

mu = dt(1.0)/dt(82.45)
lamb = dt(1.0) - mu

def gprim_1(t,z):
    return z[1]

def gprim_2(t,z):
    r1 = np.sqrt((z[0]+mu)**2 + z[2]**2)
    r2 = np.sqrt((z[0]-lamb)**2 + z[2]**2)
    return z[0] + 2*z[3] - (lamb*(z[0]+mu))/(r1**3) - (mu*(z[0]-lamb))/(r2**3)

def gprim_3(t,z):
    return z[3]

def gprim_4(t,z):
    r1 = np.sqrt((z[0]+mu)**2 + z[2]**2)
    r2 = np.sqrt((z[0]-lamb)**2 + z[2]**2)
    return z[2] - 2*z[1] - (lamb*z[2])/(r1**3) - (mu*z[2])/(r2**3)

def gprim(z,t):
    return np.array([gprim_1(t,z), gprim_2(t,z), gprim_3(t,z), gprim_4(t,z)], dtype=dt)


def cnst(z):
    r1 = math.sqrt((z[0]+mu)**2 + z[2]**2)
    r2 = math.sqrt((z[0]-lamb)**2 + z[2]**2)
    return 0.5*(z[1]**2 + z[3]**2 - z[0]**2 - z[2]**2) - lamb/r1 - mu/r2

solver = odespy.RK3(gprim, rtol=0.0)
init = [1.2, 0, 0, -1.04935750983]
solver.set_initial_condition(init)
t_points = np.linspace(0, 7, int(7/1e-4))
u, t = solver.solve(t_points)
print abs(cnst(init) - cnst(u[-1]))

solver = odespy.AdamsBashMoulton3(gprim, rtol=0.0)
init = [1.2, 0, 0, -1.04935750983]
solver.set_initial_condition(init)
t_points = np.linspace(0, 7, int(7/1e-4))
u, t = solver.solve(t_points)
print abs(cnst(init) - cnst(u[-1]))