ode integration in python versus mathematica results

Question

Edit: So I found out that NDSolve for ODE is using Runge Kutta to solve the equations. How can I use the Runge Kutta method on my python code to solve the ODE I have below?

From my post on text files with float entries, I was able to determine that python and mathematica immediately start diverging with a tolerance of 10 to the negative 6.

End Edit

For last few hours, I have been trying to figure out why my solutions in Mathematica and Python differ by 5000 something km.

I am led to believe one program has a higher error tolerance when simulating over millions of seconds in flight time.

My question is which program is more accurate, and if it isn't python, how can I adjust the precision?

With Mathematica, I am less than 10km away from L4 where as with Python I am 5947 km away.

The codes are listed below:

Python

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import linspace
from scipy.optimize import brentq

me = 5.974 * 10 ** (24)  #  mass of the earth                                     
mm = 7.348 * 10 ** (22)  #  mass of the moon                                      
G = 6.67259 * 10 ** (-20)  #  gravitational parameter                             
re = 6378.0  #  radius of the earth in km                                         
rm = 1737.0  #  radius of the moon in km                                          
r12 = 384400.0  #  distance between the CoM of the earth and moon                 
d = 300 #  distance the spacecraft is above the Earth                             
pi1 = me / (me + mm)
pi2 = mm / (me + mm)
mue = 398600.0  #  gravitational parameter of earth km^3/sec^2                    
mum = G * mm  #  grav param of the moon                                           
mu = mue + mum
omega = np.sqrt(mu / (r12 ** 3))

nu = -np.pi / 4  #  true anomaly  pick yourself                                   

xl4 = r12 / 2 - 4671  #  x location of L4                                         
yl4 = np.sqrt(3) / 2 * r12  #  y                                                  

print("The location of L4 is", xl4, yl4)

#  Solve for Jacobi's constant                                                    
def f(C):
    return (omega ** 2 * (xl4 ** 2 + yl4 ** 2) + 2 * mue / r12 + 2 * mum / r12
            + 2 * C)


c = brentq(f, -5, 0)

print("Jacobi's constant is",c)

x0 = (re + 200) * np.cos(nu) - pi2 * r12  #  x location of the satellite          
y0 = (re + 200) * np.sin(nu)  #  y location                                       

print("The satellite's initial position is", x0, y0)

vbo = (np.sqrt(omega ** 2 * (x0 ** 2 + y0 ** 2) + 2 * mue /
               np.sqrt((x0 + pi2 * r12) ** 2 + y0 ** 2) + 2 * mum /
               np.sqrt((x0 - pi1 * r12) ** 2 + y0 ** 2) + 2 * -1.21))

print("Burnout velocity is", vbo)

gamma = 0.4678 * np.pi / 180  #  flight path angle pick yourself                  

vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
#  velocity of the bo in the x direction                                          
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
#  velocity of the bo in the y direction                                          

print("The satellite's initial velocity is", vx, vy)

#  r0 = [x, y, 0]                                                                 
#  v0 = [vx, vy, 0]                                                               
u0 = [x0, y0, 0, vx, vy, 0]


def deriv(u, dt):
return [u[3],  #  dotu[0] = u[3]                                                 
        u[4],  #  dotu[1] = u[4]                                                 
        u[5],  #  dotu[2] = u[5]                                                 
        (2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
         np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
         (u[0] - pi1 * r12) /
         np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
        #  dotu[3] = that                                                        
        (-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
         np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
         np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
        #  dotu[4] = that                                                        
        0]  #  dotu[5] = 0                                                       


dt = np.linspace(0.0, 6.0 * 86400.0, 2000000.0)  #  secs to run the simulation    
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = u.T

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, color = 'r')
#  adding the Lagrange point                                                      
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = 2000 * np.outer(np.cos(phi), np.sin(theta)) + xl4
ym = 2000 * np.outer(np.sin(phi), np.sin(theta)) + yl4
zm = 2000 * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
#  adding the earth                                                               
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = 2000 * np.outer(np.cos(phi), np.sin(theta))
ym = 2000 * np.outer(np.sin(phi), np.sin(theta))
zm = 2000 * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])

plt.show()

#  The code below finds the distance between path and l4                          
my_x, my_y, my_z = (xl4, yl4, 0.0)

delta_x = x - my_x
delta_y = y - my_y
delta_z = z - my_z
distance = np.array([np.sqrt(delta_x ** 2 + delta_y ** 2 + delta_z ** 2)])

minimum = np.amin(distance)

print("Closet approach to L4 is", minimum)

Mathematica

ClearAll["Global`*"];
me = 5.974*10^(24);
mm = 7.348*10^(22);
G = 6.67259*10^(-20);
re = 6378;
rm = 1737;
r12 = 384400;

\[Pi]1 = me/(me + mm);
\[Pi]2 = mm/(me + mm);
M = me + mm;
\[Mu]1 = 398600;
\[Mu]2 = G*mm;
\[Mu] = \[Mu]1 + \[Mu]2;
\[CapitalOmega] = Sqrt[\[Mu]/r12^3];
\[Nu] = -\[Pi]/4;

xl4 = 384400/2 - 4671;
yl4 = Sqrt[3]/2*384400 // N;

Solve[\[CapitalOmega]^2*(xl4^2 + yl4^2) + 2 \[Mu]1/r12 + 
   2 \[Mu]2/r12 + 2*C == 0, C]
x = (re + 200)*Cos[\[Nu]] - \[Pi]2*r12 // N
y = (re + 200)*Sin[\[Nu]] // N

{{C -> -1.56824}}

-19.3098

-4651.35

vbo = Sqrt[\[CapitalOmega]^2*((x)^2 + (y)^2) + 
   2*\[Mu]1/Sqrt[(x + \[Pi]2*r12)^2 + (y)^2] + 
   2*\[Mu]2/Sqrt[(x - \[Pi]1*r12)^2 + (y)^2] + 2*(-1.21)]

10.8994

\[Gamma] = 0.4678*Pi/180;
vx = vbo*(Sin[\[Gamma]]*Cos[\[Nu]] - Cos[\[Gamma]]*Sin[\[Nu]]);
vy = vbo*(Sin[\[Gamma]]*Sin[\[Nu]] + Cos[\[Gamma]]*Cos[\[Nu]]);

r0 = {x, y, 0};
v0 = {vx, vy, 0}

{7.76974, 7.64389, 0}

s = NDSolve[{x1''[t] - 
      2*\[CapitalOmega]*x2'[t] - \[CapitalOmega]^2*
       x1[t] == -\[Mu]1/((Sqrt[(x1[t] + \[Pi]2*r12)^2 + 
             x2[t]^2])^3)*(x1[t] + \[Pi]2*
          r12) - \[Mu]2/((Sqrt[(x1[t] - \[Pi]1*r12)^2 + 
             x2[t]^2])^3)*(x1[t] - \[Pi]1*r12), 
    x2''[t] + 
      2*\[CapitalOmega]*x1'[t] - \[CapitalOmega]^2*
       x2[t] == -\[Mu]1/(Sqrt[(x1[t] + \[Pi]2*r12)^2 + x2[t]^2])^3*
       x2[t] - \[Mu]2/(Sqrt[(x1[t] - \[Pi]1*r12)^2 + x2[t]^2])^3*
       x2[t], x3''[t] == 0, x1[0] == r0[[1]], x1'[0] == v0[[1]], 
    x2[0] == r0[[2]], x2'[0] == v0[[2]], x3[0] == r0[[3]], 
    x3'[0] == v0[[3]]}, {x1, x2, x3}, {t, 0, 1000000}];

ParametricPlot3D[
 Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, 10*24*3600}, 
 PlotStyle -> {Red, Thick}]

g1 = ParametricPlot3D[
   Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, 5.75*3600*24}, 
   PlotStyle -> {Red}, 
   PlotRange -> {{-10000, 400000}, {-10000, 400000}}];
g2 = Graphics3D[{Blue, Opacity[0.6], Sphere[{-4671, 0, 0}, re]}];
g3 = Graphics3D[{Green, Opacity[0.6], Sphere[{379729, 0, 0}, rm]}];
g4 = Graphics3D[{Black, Sphere[{xl4, yl4, 0}, 2000]}];
Show[g2, g1, g3, g4, Boxed -> False]


(*XYdata=Flatten[Table[Evaluate[{x1[t],x2[t],x3[t]}/.s],{t,5.5*24*\
3600,5.78*24*3600,1}],1];
X1Y1data=Flatten[Table[Evaluate[{x1'[t],x2'[t],x3'[t]}/.s],{t,5.5*24*\
3600,5.78*24*3600,1}],1];
SetDirectory[NotebookDirectory[]];
Export["OrbitData.txt",XYdata,"CSV"];
Export["OrbVeloc.txt",X1Y1data,"CSV"];*)

Answer 1

If at this point your problem has reduced to just wanting to use Runge-Kutta, you can for example replace this line:

u = odeint(deriv, u0, dt)

with something like this:

#reverse the order of arguments
def deriv2(t,u):
    return deriv(u,t)

# initialize a 4th order Runge-Kutta ODE solver
solver = ode(deriv2).set_integrator('dopri5')
solver.set_initial_value(u0)
u = np.empty((len(dt), 6))
u[0,:] = u0
for ii in range(1,len(dt)):
    u[ii] = solver.integrate(dt[ii])

(+obviously replace the odeint import with ode).

Note that this is significantly slower for this type of ODE.

To use the dop853, use solver.set_integrator('dop853').

Answer 2

I re-wrote the def deriv part of the ode and it works now! So the Mathematica plot and the Python agree.

def deriv(u, dt):
    return [u[3],  #  dotu[0] = u[3]                                                 
            u[4],  #  dotu[1] = u[4]                                                 
            u[5],  #  dotu[2] = u[5]                                                 
            (2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
             np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
             (u[0] - pi1 * r12) /
             np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
            #  dotu[3] = that                                                        
            (-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
             np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
             np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
            #  dotu[4] = that                                                        
            0]  #  dotu[5] = 0      

Answer 3

One way to check the accuracy could be this:

take a time where the 2 trajectories are sufficiently different from one another (the last time/point in the trajectory for example). Use this as the new starting point and integrate back in time to the initial time (or, if you prefer, integrate forward in time, but reverse the velocities of all bodies present in the simulation). The final point in this simulation should be (close to) your initial point in the original simulation. By comparing how close you actually get to your original initial point, you can judge how good python vs mathematica solver routines are. My guess is that the Mathematica routines are much better, so all this Python thing is a waste of time.

Another way could be this: With Mathematica what you get is an interpolating function which you can symbolically derive and then numerically evaluate its derivatives. Plug these values in the ODE at different points and see if they satisfy the ODE. Do something similar in python and compare the results.