scipy odeint 中的当前迭代

Question

我正在使用 Scipy 的 odeint (scipy.integrate.odeint) 为我求解一些 ODE，一切都运行良好。但是，我现在想在我的计算中包含另一组时间相关的数据，即对于 t = [0, 1, 2, 3] 我有数据 z = [0.1, 0.2, 0.25, 0.22] 包含在计算中。我可以将向量作为参数传递，但这会为我提供每个时间步长的整个向量。有没有一种有效的方法来获取计算的当前步骤（迭代器）？这样我就可以在第 i 个时间步获得 z[i]。请注意，z 的长度为 t，并且两者都可以包含数千个元素。

谢谢

一个非常简单的例子：

import numpy as np
from scipy.integrate import odeint

def func(y, t, z):
    # I'd like to get the i-th element
    # of z, corresponding to t[i]
    return y+z[i]

result = odeint(func, [0], t, (z,))

Answer 1

此问题的解决方法是使用更通用的 scipy.integrate.ode 函数。此函数内置了多个集成方案，您可以更好地控制每次迭代期间发生的事情。请参阅以下示例：

import numpy as np
from scipy.integrate import ode

def func(t, y, z):
    return y+z

t = np.linspace(0, 1.0, 100)
dt = t[1]-t[0]
z = np.random.rand(100)
output = np.empty_like(t)
r = ode(func).set_integrator("dop853")
r.set_initial_value(0, 0).set_f_params(z[0])

for i in xrange(len(t)):
    r.set_f_params(z[i])
    r.integrate(r.t+dt)
    output[i] = r.y

在每次迭代中，z 的求解器值会相应更新。

Answer 2

对于 ode 求解器，您可以使用 interpolation 作为时变输入。在这种情况下：

import numpy as np
from scipy.integrate import odeint
from scipy.interpolate import interp1d

t_arr = np.array([0,1,2,3])
z_arr = np.array([0.1, 0.2, 0.25, 0.22])
finterp = interp1d(t_arr,z_arr,fill_value='extrapolate') #create interpolation function

def func(y,t,z):
    print(t)
    zt = finterp(t) # call interpolation at time t
    return y+zt

result = odeint(func, [0], t_arr, (z_arr,))

但是，odeint 中的求解器可能会请求超过 t_arr 的时间点（在您的情况下为 3.028），因此您必须指定超过 t=3 的 z 值或允许使用 [=17= 进行外推]（如上图）。选择一种插值时要小心，以反映您期望的行为。

scipy odeint 中的当前迭代

Current iteration in scipy odeint

python

scipy