【发布时间】：2016-01-09 10:18:08
【问题描述】：
编辑：五年了，SciPy.integrate.odeint 学会停止了吗？
下面的脚本使用 Python 中的 Runge-Kutta RK4 集成闭合路径周围的磁场线，并在它返回到某个容差内的原始值时停止。我想使用SciPy.integrate.odeint，但我不知道如何告诉它在路径接近关闭时停止。
当然odeint 可能比在 Python 中集成要快得多，我可以让它盲目地四处寻找结果，但将来我会做更大的问题。
有没有一种方法可以在 odeint 中实现“OK 足够接近 - 你现在可以停下来了！”方法？或者我应该只是整合一段时间，检查，整合更多，检查......
这个discussion 似乎是相关的，并且似乎暗示“你不能来自 SciPy”可能是答案。
注意：我通常使用 RK45 (Runge-Kutta-Fehlberg)，它在给定的台阶尺寸下更准确，以加快速度，但我在这里保持简单。它还使可变步长成为可能。
更新：但有时我需要固定步长。我发现Scipy.integrate.ode 确实提供了一种测试/停止方法ode.solout(t, y)，但似乎没有能力在t 的固定点进行评估。 odeint 允许在t 的固定点进行评估，但似乎没有测试/停止方法。
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()
【问题讨论】：
-
您现在拥有
scipy.integrate.solve_ivp，其中一种方法是 RK45（Dormand-Prince），它具有事件处理（不够灵活）和给定时间列表中的密集输出/输出插值。使用 LSODA 方法，您可以使用与 odeint 相同的积分器方法，因此从某种意义上说，它已经学会了停止。 -
@LutzLehmann 谢谢！对于“如何让 SciPy.integrate... 停止？”，这听起来确实是更好且可以接受的答案
标签： python scipy numerical-methods ode odeint