edit: It's been five years, has SciPy.integrate.odeint learned to stop yet?
The script below integrates magnetic field lines around closed paths and stops when it returns to original value within some tolerance, using Runge-Kutta RK4 in Python. I would like to use SciPy.integrate.odeint, but I can not see how I can tell it to stop when the path is approximately closed.
Of course odeint may be much faster than integrating in Python, I could just let it go around blindly and look for closure in the results, but in the future I'll do much larger problems.
Is there a way that I can implement a "OK that's close enough - you can stop now!" method into odeint? Or should I just integrate for a while, check, integrate more, check...
This discussion seems relevant, and seems to suggest that "you can't from within SciPy" might be the answer.
Note: I usually use RK45 (Runge-Kutta-Fehlberg) which is more accurate at a given steop size to speed it up, but I kept it simple here. It also makes variable step size possible.
Update: But sometimes I need fixed step size. I've found that Scipy.integrate.ode does provide a testing/stopping method ode.solout(t, y) but doesn't seem to have the ability to evaluate at fixed points of t. odeint allows evaluation at fixed points of t, but doesn't seem to have a testing/stopping method.
def rk4Bds_stops(x, h, n, F, fclose=0.1):
h_over_two, h_over_six = h/2.0, h/6.0
watching = False
distance_max = 0.0
distance_old = -1.0
i = 0
while i < n and not (watching and greater):
k1 = F( x[i] )
k2 = F( x[i] + k1*h_over_two)
k3 = F( x[i] + k2*h_over_two)
k4 = F( x[i] + k3*h )
x[i+1] = x[i] + h_over_six * (k1 + 2.*(k2 + k3) + k4)
distance = np.sqrt(((x[i+1] - x[0])**2).sum())
distance_max = max(distance, distance_max)
getting_closer = distance < distance_old
if getting_closer and distance < fclose*distance_max:
watching = True
greater = distance > distance_old
distance_old = distance
i += 1
return i
def get_BrBztanVec(rz):
Brz = np.zeros(2)
B_zero = 0.5 * i * mu0 / a
zz = rz[1] - h
alpha = rz[0] / a
beta = zz / a
gamma = zz / rz[0]
Q = ((1.0 + alpha)**2 + beta**2)
k = np.sqrt(4. * alpha / Q)
C1 = 1.0 / (pi * np.sqrt(Q))
C2 = gamma / (pi * np.sqrt(Q))
C3 = (1.0 - alpha**2 - beta**2) / (Q - 4.0*alpha)
C4 = (1.0 + alpha**2 + beta**2) / (Q - 4.0*alpha)
E, K = spe.ellipe(k**2), spe.ellipk(k**2)
Brz[0] += B_zero * C2 * (C4*E - K)
Brz[1] += B_zero * C1 * (C3*E + K)
Bmag = np.sqrt((Brz**2).sum())
return Brz/Bmag
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as spe
from scipy.integrate import odeint as ODEint
pi = np.pi
mu0 = 4.0 * pi * 1.0E-07
i = 1.0 # amperes
a = 1.0 # meters
h = 0.0 # meters
ds = 0.04 # step distance (meters)
r_list, z_list, n_list = [], [], []
dr_list, dz_list = [], []
r_try = np.linspace(0.15, 0.95, 17)
x = np.zeros((1000, 2))
nsteps = 500
for rt in r_try:
x[:] = np.nan
x[0] = np.array([rt, 0.0])
n = rk4Bds_stops(x, ds, nsteps, get_BrBztanVec)
n_list.append(n)
r, z = x[:n+1].T.copy() # make a copy is necessary
dr, dz = r[1:] - r[:-1], z[1:] - z[:-1]
r_list.append(r)
z_list.append(z)
dr_list.append(dr)
dz_list.append(dz)
plt.figure(figsize=[14, 8])
fs = 20
plt.subplot(2,3,1)
for r in r_list:
plt.plot(r)
plt.title("r", fontsize=fs)
plt.subplot(2,3,2)
for z in z_list:
plt.plot(z)
plt.title("z", fontsize=fs)
plt.subplot(2,3,3)
for r, z in zip(r_list, z_list):
plt.plot(r, z)
plt.title("r, z", fontsize=fs)
plt.subplot(2,3,4)
for dr, dz in zip(dr_list, dz_list):
plt.plot(dr, dz)
plt.title("dr, dz", fontsize=fs)
plt.subplot(2, 3, 5)
plt.plot(n_list)
plt.title("n", fontsize=fs)
plt.show()
