我有一条穿过时空的光线,即一条曲线,由某个变量 λ 参数化。确切的轨迹,即曲线的坐标函数由一些 ODE。
现在,我的代码要满足两个要求:
结果必须在规则网格上给出,即步长 dλ 必须在任何地方都相同。
我不知道区间 λ 预先运行,即积分必须在满足某个条件后立即动态结束(基本上我正在测试坐标是否在我不知道时空属性的某个坐标区间之外)。
要求 2 显然不包括scipy.odeint()
. 这给我留下了一个类,如果我选择或算法scipy.ode
,它允许我通过solout
回调将中断条件传递给积分器。dopri5
dop853
不过,这有点复杂:由于这些(和所有其他)积分器动态调整它们的内部步长(这是不手动进行积分的一种点),如果我愿意,我不能简单地存储传递给回调的值满足要求1。(更新:事实上,这样做是没有意义的。请参阅下面的我的注释(*)。)相反,我仍然需要循环λ并scipy.ode.integrate()
为每个步骤调用方法。然而,这意味着我必须以某种方式停止循环:
from scipy.integrate import ode
def deriv(λ, x, params):
# …
return dx_dλ
end_reached = False
def solout(λ, x):
if g(λ, x):
# Stop integration
# Returning -1 will NOT make solver.successful() return False, so
# we need to inform the integration loop (see below) separately.
nonlocal end_reached
end_reached = True
return -1
# Everything's ok. Move on.
return 0
λ0 = 0
x0 = …
params = …
dλ = 0.1
solver = ode(deriv)
solver.set_integrator('dopri5', nsteps=500)
solver.set_solout(solout)
solver.set_initial_value(x0, λ0)
solver.set_f_params(params)
λ = [λ0]
x = [x0]
while True:
solver.integrate(solver.t + dλ)
if not solver.successful() or end_reached:
break
λ.append(solver.t)
x.append(solver.y)
总而言之,这是相当丑陋的。不仅逻辑基本上无处不在(例如,我无法测试它所属的循环中的条件)。但是由于solout
回调(我想),它也至少比没有回调慢一个数量级。更新:事实上,即使没有回调,dopri5
积分器也比vode
(至少对于我想使用这种方法的另一个 ODE)慢几个数量级。
我是否遗漏了什么,或者这真的是满足上述要求的唯一方法吗?我还尝试引发异常deriv()
(并随后在循环中捕获它)以停止积分器。虽然这适用于所有集成商,但它仍然会导致一堆错误被打印到标准输出。
(*) 据我了解,传递给的值deriv()
并不是一个可靠的近似值,而只是求解器返回的最终值。例如,考虑到(快速)变化的斜率,求解器可能会进行过冲以选择合理的步长,从而得出最佳估计。