# Numerical Programming using ODEINT takes more than 17 minutes

I am trying to track the trajectories of a charged particles under the influence of Gravitational and Electromagnetic effect. Computing for time points, t0=0second -tf= (36002430)sec with stepsize 0.5 more than 17 minutes. Is it normal for a programming to take such a long time or is it something to do with my laptop (core i3, 8th gen, 2.54Hz) or my code. Please someone review my code:

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi
from scipy.integrate import odeint

μ = 379312077e8           # m^3/s^2 # G= 6.67408e-11m^3/kg*s^ # M = 5.68e26 # μ=GM
R = 60268e3               # metre
g_10 = 21141e-9
Ω = 17e-5                 # rad/s = 9.74e-7

# charge by mass ratio
ρ = 1e3                   # 1gm/cm^3 Headmen 2021= 10^3 kg/m^3,
b = 1e-11                 # 1nm = 1e-3 micro
V = 100                   # Volt
β = (3*ε*V)/(ρ*b**2)     #
m = (4/3)*pi*(b**3)*ρ

def LzG(p,t):
# assigning each ODE to a vector element
r,x,θ,y,ϕ,z = p

# constants
μ = 379312077e9                   # m^3/s^2, G = 6.67408e-11m^3/kg*s^,M = 5.68e26kg
R = 60268e3                          # metre
j_2 = 1.629071e-2
g_10 = 21141e-9                      #
μ_0 = 4*pi*1e-7
B_θ = μ_0*(R/r)**3*g_10*sin(θ)
B_r = μ_0*2*(R/r)**3*g_10*cos(θ)
β = (3*ε*V)/(ρ*b**2)                # q/m = -3.46 x 103 C/kg,

# defining the ODEs
drdt = x
dxdt = r*(y**2+(z+Ω)**2*sin(θ)**2-β1*z*sin(θ)*B_θ)-(μ/r**2)*(1-(3/2)*j_2*(R/r)**2*(3*cos(θ)**2-1))
dθdt = y
dydt = (-2*x*y + r*(z+Ω)**2*sin(θ)*cos(θ)+ β1*r*z*sin(θ)*B_r)/r + (3*μ/r**2)*j_2*(R/r)**2*sin(θ)*cos(θ)
dϕdt = z
dzdt = (-2*(z+Ω)*(x*sin(θ)+r*y*cos(θ)) + β1*(x*B_θ-r*y*B_r))/(r*sin(θ))

return np.array([drdt,dxdt,dθdt,dydt,dϕdt,dzdt])


# initial conditions

p0 = np.array([1.11*R,0.0,90.0*(pi/180),0.0*(pi/180), 0.0*(pi/180), 0.0206*(pi/180)])


# time window

t = np.arange(0,2680000,0.5)# 30 days

# Solutions
p_LzG = odeint(LzG, p0,t,rtol=1e-4,atol=1e-4)

r,x,θ,y,ϕ,z = p_LzG.T


# Plot Solutions

r,x,θ,y,ϕ,z = p_LzG.T
fig,ax=plt.subplots(2,3,figsize=(10,4))

for a,s in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
a.plot(t,s); a.grid()
plt.tight_layout(); plt.show()

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

fig = plt.figure()


# spherical to cartesian

x1= r * np.sin(θ)* np.cos(ϕ)
y1= r * np.sin(θ)* np.sin(ϕ)
z1= r * np.cos(θ)

ax = plt.axes(projection="3d")
ax.plot3D(x1,y1,z1, color="red",linewidth='0.5')
ax.set_xlabel('X-Axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
plt.show()


OSERVATIONS AFTER CHANGES MADE AS PER SUGGESIONS 1.Tolerance rather than timestep reduces the computing time to seconds 2.With the tolerance to rtol=1e-3 and atol=1e-3 it takes about 10sec. 3. The solutions/graphs with (rtol=1e-3 and atol=1e-3) and ( rtol=1e-4 and atol=1e-4 and lower values) are different ! 4. With tolerance 1e-2 and above the solution suddenly falls off to Zero, with default (None) No Solution. 5. With regard to the time stepsize, strangely ODEINTS seems to accepts only certain values (in this case) <1, 10, 36 and 100 6. solution of ODEINT and solve_ivp are different ! NB: The equations of motion are in a co-Rotating frame.     I have accepted the answers since the question was related to computing time and has been significantly reduces with the help of the David Hammen and Uhoh. However, I still doubt the robustness of ODEINT and Solve_ivp with regard to complicated coupled Non linear non- homogenous equations such as this and for long duration simulations (without re-interpolating the internal results).

• You're using python, which is ten to a hundred times slower than a compiled language. You're also using a rather small step size. 72 million steps is going to take quite a bit of time, even in a compiled language, let alone python. You have two inconsistent definitions of mu that differ by a factor of ten. The first one is correct (for Saturn). The overriding value in the derivative function (which is the value that is used) is incorrect. Even stranger yet, your ODE looks ... strange. Mar 12 at 10:29
• A much better approach is to take larger (much, much larger) step sizes for integration and interpolate if needed. Python scipy's odeint and solve_ivp are fairly good and fairly fast, but only if you give them nice step sizes. Adaptive numerical integrators perform much better, both in terms of accuracy and CPU usage, with large step sizes. Mar 12 at 10:34
• @DavidHammen ODENT chooses and varies its own step size, based on the value of rtol (relative tolerance) and atol (absolute tolerance) So I guess right now atol=1E-04 meters (or m/sec) is requesting 100 micron accuracy per integration step!
– uhoh
Mar 12 at 17:03
• Also, as you can see on the documentation page I linked to in my comment to DH, they recommend that new code use scipy.integrate.solve_ivp ODEINT is still supported, but solve_ivp is "better" and will be better supported in the future.
– uhoh
Mar 12 at 17:05
• @uhoh You appear to have a misunderstanding of atol. The intent is to address integrated values that become less than one, where rtol is not a good specifier. Keep in mind that the default values for rtol and atol are 1.49e-8. odeint cannot vary its step size if the user-requested step size is overly small. That the factor of ten error in Saturn's gravitational parameter (mu) might well have resulted in a highly elliptical orbit, and maybe odeint tried hard to zero in on the periapsis time. However, I suspect the key culprit was the request for a minimum of 72 million steps. Mar 12 at 18:05

My first recommendation is to drastically increase the step size. Both odeint and solve_ivp use adaptive integrators. Keep in mind that the integrator will take at least one step per the user-desired step size. Your step size (half a second, or about 35000 steps per orbit) is so small that it precludes that adaptation. You might as well be using the Euler method (which is incredibly lousy) for integration at that small of a step size. Even generating the plots for that much collected data will take a good amount of CPU time.

A hundred or so steps per orbit (which means a step size of two to three minutes or so rather than half of a second) should be more than enough for visualization purposes, and this much larger step size will enable an adaptive integrator to come much closer to finding the optimal internal step size. Your 17 minute run time should drop by a factor of ten, minimum, maybe by a factor of one hundred -- and your results will be much more accurate.

Regarding rtol and atol, the first (rtol) addresses relative accuracy. Setting rtol to 1e-4 means you are going after four decimal digits of precision. I suggest leaving that alone; use the default. The default is to pursue about seven or eight decimal places of relative accuracy. You might be able to do much better. Don't make it worse! The second (atol) addresses absolute accuracy. The intent is to let values between -1.0 and +1.0 be a bit sloppy. Setting atol to 1e-4 is perhaps a bit too sloppy. I suggest you leave atol alone as well. Note: the relative tolerance (rtol) dominates for integrated quantities whose absolute values are well above one while the absolute tolerance (atol) dominates for integrated quantities whose absolute values are well below one.

• Tolerance rather than timestep reduces the computing time to seconds (1)With the tolerance to rtol=1e-3 and atol=1e-3 it takes about 10sec. (2) The solutions/graphs with (rtol=1e-3 and atol=1e-3) and ( rtol=1e-4 and atol=1e-4 and lower values) are different ! (3) With tolerance 1e-2 and above the solution suddenly falls off to Zero, with default (None) No Solution. (4) With regard to the time stepsize, strangely ODEINTS seems to accepts only certain values (in this case) <1, 10, 36 and 100 (5) solution of ODEINT and solve_ivp are different ! Mar 14 at 15:42
• You don't want to change the tolerance like that. It's a bad idea as it results in getting completely garbage results quickly. Mar 15 at 0:14
• (reposting this comment because I forgot to include the quote the first time) @DavidHammen "...the integrator will take at least one step per the user-desired step size" is not really right. The number of function calls is essentially independent of the size of the time array in the call. The extra time spent at small "requested step size" is likely to be the time spent in post-interpolation. See Understanding SciPy's odeint's step sizes, evaluation times number of function calls and total process time (context is a question in AstronomySE)
– uhoh
Mar 15 at 21:49
• Qhether the time array size is 11 or 1000001, there's about 500 function calls. I think best way to proceed is to use solve_ivp instead of odeint, request a coarse spacing for output, find the regions of interest and then having set dense_output=True in the call, the interpolation coefficients will be exposed and can be used to get results in higher density in those regions of interest, or at a later time (like how Horizons works).
– uhoh
Mar 15 at 21:53
• @DavidHammen, the intent here is not about manipulating and getting quick results. But rather, the vulnerability and reliability of the integrator (odeint) with change of the tolerance arguments. Mar 16 at 11:18

I would suggest taking a look at JIT acceleration, e.g. with the Numba package, Cython or others, which help compile the Python code more efficiently. You also have other integration libraries for Python which can be much faster, e.g. numbakit-ode.

Beware when increasing the tolerance or the time step, as your solution may become quite wrong...