I am modeling charged particles moving through an electromagnetic field and am using scipy ode. The code here is simplified, obviously, but works as an example. The problem I have is that I want to end the integration after a limit on r, not on t. So, integrate dx/dt up to the point where norm(x) > r.
I don't want to just change the function to integrate over r, however, because the position is a function of t. Can I do a definite integral over an unrelated variable or something?
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def RHS(state, t, Efield, q, mp):
ds0 = state[3]
ds1 = state[4]
ds2 = state[5]
ds3 = q/mp * Efield*state[0]
ds4 = q/mp * Efield*state[1]
ds5 = q/mp * Efield*state[2]
r = np.linalg.norm((state[0], state[1], state[2]))
# if r > 30000 then do stop integration.....?
# return the two state derivatives
return [ds0, ds1, ds2, ds3, ds4, ds5]
ts = np.arange(0.0, 10.0, 0.1)
state0 = [1.0, 2.0, 3.0, 0.0, 0.0, 0.0]
Efield=1.0
q=1.0
mp=1.0
stateFinal = odeint(RHS, state0, ts, args=(Efield, q, mp))
print(np.linalg.norm(stateFinal[-1,0:2]))
You can control the process by performing integration step by step using scipy.integrate.ode. Example:
Note that the signature of RHS will have to be changed to
def RHS(t, state, Efield, q, mp):
forode
, the independent variable comes first, unlike inodeint
.The output is the solution computed in the increments of
dt
of independent variable, up to the time when the loop ends (either becauset_max
is reached, or integrator fails, or the condition forbreak
was encountered).