<p>My first recommendation is to drastically increase the step size. Both <code>odeint</code> and <code>solve_ivp</code> 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.</p>
<p>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.</p>
<p>Regarding <code>rtol</code> and <code>atol</code>, the first (<code>rtol</code>) addresses relative accuracy. Setting <code>rtol</code> 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 (<code>atol</code>) addresses absolute accuracy. The intent is to let values between -1.0 and +1.0 be a bit sloppy. Setting <code>atol</code> to 1e-4 is perhaps a bit too sloppy. I suggest you leave <code>atol</code> alone as well. Note: the relative tolerance (<code>rtol</code>) dominates for integrated quantities whose absolute values are well above one while the absolute tolerance (<code>atol</code>) dominates for integrated quantities whose absolute values are well below one.</p>

CC BY-SA 4.0
2023-03-14T15:42:33.800Z
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 !
Lunthang Peter
CC BY-SA 4.0
2023-03-15T00:14:33.607Z
You don't want to change the tolerance like that. It's a bad idea as it results in getting completely garbage results quickly.
David Hammen
CC BY-SA 4.0
2023-03-15T21:49:24.633Z
(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)](https://scicomp.stackexchange.com/q/42613/17869)
uhoh
CC BY-SA 4.0
2023-03-15T21:53:45.617Z
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
CC BY-SA 4.0
2023-03-16T11:18:03.813Z
@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.
Lunthang Peter
CC BY-SA 4.0
2023-03-14T13:44:45.917Z
2023-03-14T13:44:45.917Z
David Hammen
<p>I would suggest taking a look at JIT acceleration, e.g. with the <a href="https://numba.pydata.org/" rel="nofollow noreferrer">Numba</a> package, <a href="https://cython.org/" rel="nofollow noreferrer">Cython</a> 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. <a href="https://github.com/hgrecco/numbakit-ode" rel="nofollow noreferrer">numbakit-ode</a>.</p>
<p>Beware when increasing the tolerance or the time step, as your solution may become quite wrong...</p>

CC BY-SA 4.0
2023-03-23T17:20:03.820Z
2023-03-23T22:38:24.840Z
2023-03-23T22:38:24.840Z
Laurent90
<p>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= (3600<em>24</em>30)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:</p>
<pre><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
ε = 8.85e-12              # Farad/metre        
β = (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                      #
    Ω = 1.7e-4                           # rad/second
    μ_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])
</code></pre>
<h1>initial conditions</h1>
<pre><code>p0 = np.array([1.11*R,0.0,90.0*(pi/180),0.0*(pi/180), 0.0*(pi/180), 0.0206*(pi/180)])
</code></pre>
<h1>time window</h1>
<pre><code>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
</code></pre>
<h1>Plot Solutions</h1>
<pre><code>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()
</code></pre>
<h1>spherical to cartesian</h1>
<pre><code>x1= r * np.sin(θ)* np.cos(ϕ)
y1= r * np.sin(θ)* np.sin(ϕ)
z1= r * np.cos(θ)

ax = plt.axes(projection=&quot;3d&quot;)
ax.plot3D(x1,y1,z1, color=&quot;red&quot;,linewidth='0.5')
ax.set_xlabel('X-Axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
plt.show()
</code></pre>
<p><strong>OSERVATIONS AFTER CHANGES MADE AS PER SUGGESIONS</strong>
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) &lt;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.</p>
<p><a href="https://i.sstatic.net/nwe7P.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/nwe7P.png" alt="For, t = np.arange(0,86400,100) rtol=None,atol=None" /></a></p>
<p><a href="https://i.sstatic.net/4Oc8e.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/4Oc8e.png" alt="t = np.arange(0,86400,100) rtol=1e-4,atol=1e-4" /></a></p>
<p><a href="https://i.sstatic.net/gM0g6.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gM0g6.png" alt="rtol=1e-3,atol=1e-3" /></a></p>
<p><a href="https://i.sstatic.net/ldlsR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ldlsR.png" alt="rtol=1e-3,atol=1e-3" /></a></p>
<p><a href="https://i.sstatic.net/vsHDw.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/vsHDw.png" alt="t = np.linspace(t_span[0], t_span1, 100)
sol = solve_ivp(odes,t_span ,p0,t_eval= t,method= &quot;DOP853&quot;,dense_output=True, vectorised= False, atol=1e-4,rtol=1e" /></a></p>
<p>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).</p>

CC BY-SA 4.0
2023-03-12T09:55:35.157Z
have you tried with  scipy.integrate.solve_ivp  which is said to be better?
James K
CC BY-SA 4.0
2023-03-12T09:59:13.280Z
@James K, Thank you for the suggestion. I shall change the tolerance to 0.01 and get back to you. Yes, I have tried scipy.integrate.solve_ivp before ODEINT but gave me similar troubles.
Lunthang Peter
CC BY-SA 4.0
2023-03-12T10:29:03.300Z
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.
David Hammen
CC BY-SA 4.0
2023-03-12T10:34:53.287Z
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.
David Hammen
CC BY-SA 4.0
2023-03-12T17:03:41.017Z
@DavidHammen [ODENT](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) 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
CC BY-SA 4.0
2023-03-12T17:03:51.107Z
Hi @LunthangPeter see my comment above (and the documentation!) can you remove `atol=1e-4` from your call to odeint, or set `atol=None` (which is the default anyway) and see if that speeds it up? In my experience with NumPy and SciPy it's almost always better to go with the defaults first, and not set parameters you don't need to. Also if you set `full_output=True` you will get a python dictionary-like object that has lots of goodies, and it will tell you how many function calls have been made. I think it's making zillions more than it needs to because of your 1E-04 kilometer specificlation
uhoh
CC BY-SA 4.0
2023-03-12T17:05:53.587Z
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`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) ODEINT is still supported, but solve_ivp is "better" and will be better supported in the future.
uhoh
CC BY-SA 4.0
2023-03-12T17:58:37.320Z
@uhoh If one forces the step size to be small, which is exactly what this user did, it will use that overly small step size rather than optimal step size. The user asked for 72 million steps, and I suspect that that is exactly how many steps were taken. A good integrator does not need that many steps for 30 days in Saturn orbit.
David Hammen
CC BY-SA 4.0
2023-03-12T18:05:04.857Z
@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.
David Hammen
CC BY-SA 4.0
2023-03-12T19:25:23.777Z
I concur with @DavidHammen: I think OP is trying to do too much here, too quickly. Science computing isn't "hack this complex problem into code and instantly reproduce paper". You'd do good to first simplify your code to test it. 1.) No forces: Do you get straight lines, as expected? 2.) Only gravity with $\vec v_0=\vec 0$: can you reproduce the free-fall velocities, do the units all check out? 3.) Only EM, no gravity, do you get correct gyroradii etc.    and of course all that with a number of timesteps that gives you results in a few seconds..
AtmosphericPrisonEscape
CC BY-SA 4.0
2023-03-12T23:10:35.833Z
@DavidHammen No, I am not convinced I have misunderstood the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) "The solver will control the vector, e, of estimated local errors in y, according to an inequality of the form max-norm of (e / ewt) <= 1, where ewt is a vector of positive error weights computed as ewt = rtol * abs(y) + atol." If we worked in some kind of reduced units (with say 1 year = 1 AU = 1) which one really should try to do, then 1.49E-8 would make sense as a default to get started. Once my coffee kicks in I'll run some tests.
uhoh
CC BY-SA 4.0
2023-03-12T23:20:19.190Z
@DavidHammen again I think, for the first situation ever, you're not quite right! :-) The step size is internal and not visible unless the user does as I suggest and sets `full_output=True`. The line `t = np.arange(0,2680000,0.5)# 30 days` specifies the *interpolation grid*, something that's done after the integration is complete. If they then do as I suggest and use `scipy.integrate.solve_ivp` they will then be able to set `dense_output=True` and expose the raw results directly from the integrator, including the interpolation coefficients, in case they'd like to interpolate differently later.
uhoh
CC BY-SA 4.0
2023-03-12T23:22:25.733Z
@DavidHammen ODEINT is what I grew up on. Except for a hand-rolled RK45 once (for fun and learning) it's the only numerical integrator I've ever used, until just recently when I switched to the (very similar cousin) solve_ivp.
uhoh
CC BY-SA 4.0
2023-03-13T06:59:59.670Z
@DavidHammen and uhoh, Thank you both for your edifying discussion. The computing times reduces to seconds by changing the tolerance to 1e-3. Tolerance from  1e-4 and below including None (defaults)  it takes about 17minutes. With tolerance 1e-2 and above the solution fall off to Zero, which is not the case with smaller values of tolerance. The same is with solve_ivp. Also with solve_ivp the graph is so dense that it is difficult to decipher the physics.
Lunthang Peter
CC BY-SA 4.0
2023-03-13T07:00:06.997Z
Another interesting fact which I cant understand isthat with smaller time step (0.5s) the instability point ( the point in the graph) where the particle leaves the orbit happens around time 72,00000. But with larger step size (36s) the instability point is obtained at 2711750s. It appears that the instaibility depends on the choice of the timestep and not on the physics !.
Lunthang Peter
CC BY-SA 4.0
2023-03-13T07:19:43.323Z
@DavidHamme, would you please specify " The overriding value in the derivative function" which is incorrect? the mu error has been noted and corrected.
Lunthang Peter
CC BY-SA 4.0
2023-03-13T07:29:11.063Z
@LunthangPeter now I'm confused; per my best understanding, odeint does not allow you to set the integration step size, that's internal to the integrator. The times you specify are only for re-interpolating the internal results to a new set of time points you're interested in. The only way I see to affect step size is through the parameters h0, hmax hmin, and mxstep, and that's affect, not set. Make your call `p_LzG, infodict = odeint(..., full_output=True)` The infodict object is now a dictionary-like thing you can query for internal steps, e.g. `infodict['hu']` returns an array of step sizes
uhoh
CC BY-SA 4.0
2023-03-13T07:59:52.617Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/144559/discussion-between-lunthang-peter-and-uhoh).
Lunthang Peter
CC BY-SA 4.0
2023-03-15T01:58:19.810Z
@LunthangPeter have a look at [Understanding SciPy's odeint's step sizes, evaluation times number of function calls and total process time (context is a question in AstronomySE)](https://scicomp.stackexchange.com/q/42613/17869)
uhoh
<astrophysics><planetary-science><modeling><plasma-physics><saturn-rings>
Numerical Programming using ODEINT takes more than 17 minutes
Lunthang Peter
# Numerical Programming using ODEINT takes more than 17 minutes

(Asked by: username_0 on 2023-03-12)

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= (3600*24*30)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
ε = 8.85e-12              # Farad/metre        
β = (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                      #
    Ω = 1.7e-4                           # rad/second
    μ_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.

[![For, t = np.arange(0,86400,100) rtol=None,atol=None](https://i.sstatic.net/nwe7P.png)](https://i.sstatic.net/nwe7P.png)

[![t = np.arange(0,86400,100) rtol=1e-4,atol=1e-4](https://i.sstatic.net/4Oc8e.png)](https://i.sstatic.net/4Oc8e.png)

[![rtol=1e-3,atol=1e-3](https://i.sstatic.net/gM0g6.png)](https://i.sstatic.net/gM0g6.png)

[![rtol=1e-3,atol=1e-3](https://i.sstatic.net/ldlsR.png)](https://i.sstatic.net/ldlsR.png)

[![t = np.linspace(t_span[0], t_span1, 100)
sol = solve_ivp(odes,t_span ,p0,t_eval= t,method= "DOP853",dense_output=True, vectorised= False, atol=1e-4,rtol=1e](https://i.sstatic.net/vsHDw.png)](https://i.sstatic.net/vsHDw.png)

I have accepted the answers since the question was related to computing time and has been significantly reduces with the help of the username_2 username_2 and username_3. 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).

**Tags:** astrophysics, planetary-science, modeling, plasma-physics, saturn-rings

> **Comments:**
> **username_1**: have you tried with  scipy.integrate.solve_ivp  which is said to be better?  
> **username_0**: @username_1, Thank you for the suggestion. I shall change the tolerance to 0.01 and get back to you. Yes, I have tried scipy.integrate.solve_ivp before ODEINT but gave me similar troubles.  
> **username_2**: 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.  
> **username_2**: 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.  
> **username_3**: @username_2 [ODENT](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) 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!*  
> **username_3**: Hi @username_0 see my comment above (and the documentation!) can you remove `atol=1e-4` from your call to odeint, or set `atol=None` (which is the default anyway) and see if that speeds it up? In my experience with NumPy and SciPy it's almost always better to go with the defaults first, and not set parameters you don't need to. Also if you set `full_output=True` you will get a python dictionary-like object that has lots of goodies, and it will tell you how many function calls have been made. I think it's making zillions more than it needs to because of your 1E-04 kilometer specificlation  
> **username_3**: 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`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) ODEINT is still supported, but solve_ivp is "better" and will be better supported in the future.  
> **username_2**: @username_3 If one forces the step size to be small, which is exactly what this user did, it will use that overly small step size rather than optimal step size. The user asked for 72 million steps, and I suspect that that is exactly how many steps were taken. A good integrator does not need that many steps for 30 days in Saturn orbit.  
> **username_2**: @username_3 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.  
> **username_4**: I concur with @username_2: I think OP is trying to do too much here, too quickly. Science computing isn't "hack this complex problem into code and instantly reproduce paper". You'd do good to first simplify your code to test it. 1.) No forces: Do you get straight lines, as expected? 2.) Only gravity with $\vec v_0=\vec 0$: can you reproduce the free-fall velocities, do the units all check out? 3.) Only EM, no gravity, do you get correct gyroradii etc.    and of course all that with a number of timesteps that gives you results in a few seconds..  
> **username_3**: @username_2 No, I am not convinced I have misunderstood the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) "The solver will control the vector, e, of estimated local errors in y, according to an inequality of the form max-norm of (e / ewt) <= 1, where ewt is a vector of positive error weights computed as ewt = rtol * abs(y) + atol." If we worked in some kind of reduced units (with say 1 year = 1 AU = 1) which one really should try to do, then 1.49E-8 would make sense as a default to get started. Once my coffee kicks in I'll run some tests.  
> **username_3**: @username_2 again I think, for the first situation ever, you're not quite right! :-) The step size is internal and not visible unless the user does as I suggest and sets `full_output=True`. The line `t = np.arange(0,2680000,0.5)# 30 days` specifies the *interpolation grid*, something that's done after the integration is complete. If they then do as I suggest and use `scipy.integrate.solve_ivp` they will then be able to set `dense_output=True` and expose the raw results directly from the integrator, including the interpolation coefficients, in case they'd like to interpolate differently later.  
> **username_3**: @username_2 ODEINT is what I grew up on. Except for a hand-rolled RK45 once (for fun and learning) it's the only numerical integrator I've ever used, until just recently when I switched to the (very similar cousin) solve_ivp.  
> **username_0**: @username_2 and username_3, Thank you both for your edifying discussion. The computing times reduces to seconds by changing the tolerance to 1e-3. Tolerance from  1e-4 and below including None (defaults)  it takes about 17minutes. With tolerance 1e-2 and above the solution fall off to Zero, which is not the case with smaller values of tolerance. The same is with solve_ivp. Also with solve_ivp the graph is so dense that it is difficult to decipher the physics.  
> **username_0**: Another interesting fact which I cant understand isthat with smaller time step (0.5s) the instability point ( the point in the graph) where the particle leaves the orbit happens around time 72,00000. But with larger step size (36s) the instability point is obtained at 2711750s. It appears that the instaibility depends on the choice of the timestep and not on the physics !.  
> **username_0**: @DavidHamme, would you please specify " The overriding value in the derivative function" which is incorrect? the mu error has been noted and corrected.  
> **username_3**: @username_0 now I'm confused; per my best understanding, odeint does not allow you to set the integration step size, that's internal to the integrator. The times you specify are only for re-interpolating the internal results to a new set of time points you're interested in. The only way I see to affect step size is through the parameters h0, hmax hmin, and mxstep, and that's affect, not set. Make your call `p_LzG, infodict = odeint(..., full_output=True)` The infodict object is now a dictionary-like thing you can query for internal steps, e.g. `infodict['hu']` returns an array of step sizes  
> **username_0**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/144559/discussion-between-username_0-username_0-and-username_3).  
> **username_3**: @username_0 have a look at [Understanding SciPy's odeint's step sizes, evaluation times number of function calls and total process time (context is a question in AstronomySE)](https://scicomp.stackexchange.com/q/42613/17869)  

---

## 2 answers

### Answer by username_2 on 2023-03-14. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_0**: 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 !  
> **username_2**: You don't want to change the tolerance like that. It's a bad idea as it results in getting completely garbage results quickly.  
> **username_3**: (reposting this comment because I forgot to include the quote the first time) @username_2 "...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)](https://scicomp.stackexchange.com/q/42613/17869)  
> **username_3**: 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).  
> **username_0**: @username_2, 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.  

---

### Answer by username_5 on 2023-03-23. Score: 1

I would suggest taking a look at JIT acceleration, e.g. with the [Numba](https://numba.pydata.org/) package, [Cython](https://cython.org/) 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](https://github.com/hgrecco/numbakit-ode).

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


---

