integrate.ode设置t0值超出数据范围

Question

我想解决这个微分方程 dy/dt = -2y + data(t)，在 t=0 到 t=3 之间，初始条件是 y(t=0)=1。

我写了以下代码：

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

t = np.linspace(0, 3, 4)

data = [1, 2, 3, 4]

linear_interpolation = interp1d(t, data)

def func(y, t0):
    print 't0', t0
    return -2*y + linear_interpolation(t0)

soln = odeint(func, 1, t)

当我运行这段代码时，出现了几个错误，像这样：

ValueError: x_new 中的一个值超出了插值范围。
odepack.error: 调用名为 func 的 Python 函数时发生错误

我的插值范围是 0.0 到 3.0。打印 func 中的 t0 的值时，我发现 t0 有时确实超出了我的插值范围，比如：3.07634612585, 3.0203768998, 3.00638459329，等等。这就是为什么 linear_interpolation(t0) 会引发 ValueError 异常的原因。

我有几个问题：

• integrate.ode 是怎么让 t0 变化的？为什么它会让 t0 超过我插值范围的上限（3.0）？

• 尽管出现了这些错误，integrate.ode 还是返回了一个看起来正确的值的数组。那么，我应该只是捕获并忽略这些错误吗？无论微分方程、t 范围和初始条件是什么，我都应该忽略它们吗？

• 如果我不应该忽略这些错误，避免它们的最佳方法是什么？我有两个建议：

  • 在 interp1d 中，我可以设置 bounds_error=False 和 fill_value=data[-1]，因为超出我插值范围的 t0 看起来接近 t[-1]：

    linear_interpolation = interp1d(t, data, bounds_error=False, fill_value=data[-1])
    

    但首先我想确保无论其他 func 和其他 data 是什么，t0 都会始终接近 t[-1]。例如，如果 integrate.ode 选择了一个低于我插值范围的 t0，那么填充值仍然是 data[-1]，这就不正确了。也许了解 integrate.ode 是如何让 t0 变化的会帮助我确认这一点（见我的第一个问题）。

  • 在 func 中，我可以把 linear_interpolation 的调用放在一个 try/except 块里，当我捕获到 ValueError 时，我重新调用 linear_interpolation，但用截断后的 t0：

    def func(y, t0):    
        try:
            interpolated_value = linear_interpolation(t0)
        except ValueError:
            interpolated_value = linear_interpolation(int(t0)) # truncate t0
        return -2*y + interpolated_value
    

    至少这个解决方案允许 linear_interpolation 在 integrate.ode 让 t0 >= 4.0 或 t0 <= -1.0 时仍然引发异常。这样我就能被提醒到不一致的行为。但这并不是很清晰，而且截断现在看起来有点随意。

也许我只是对这些错误想得太多了。请告诉我你的看法。

Answer 1

在使用odeint这个求解器时，看到它在你请求的最后时间之后评估你的函数是很正常的。大多数常微分方程（ODE）求解器都是这样工作的——它们会根据内部的错误控制算法，选择一些时间步长，然后用自己的方法来计算用户请求的时间点的解。有些求解器（比如Sundials库中的CVODE求解器）允许你设定一个时间的上限，超过这个时间求解器就不再评估你的方程，但odeint没有这个选项。

如果你不介意从scipy.integrate.odeint切换到scipy.integrate.ode，那么"dopri5"和"dop853"这两个求解器就不会在请求的时间之后评估你的函数。不过有两个注意事项：

• ode求解器在定义微分方程的参数顺序上有不同的约定。在ode求解器中，t是第一个参数。（是的，我知道，这让人有点不爽...）
• "dopri5"和"dop853"求解器适用于非刚性系统。如果你的问题是刚性的，它们仍然会给出正确的答案，但会比刚性求解器做更多的计算。

这里有一个脚本，展示了如何解决你的例子。为了强调参数的变化，我把func改名为rhs。

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


t = np.linspace(0, 3, 4)
data = [1, 2, 3, 4]
linear_interpolation = interp1d(t, data)

def rhs(t, y):
    """The "right-hand side" of the differential equation."""
    #print 't', t
    return -2*y + linear_interpolation(t)


# Initial condition
y0 = 1

solver = ode(rhs).set_integrator("dop853")
solver.set_initial_value(y0)

k = 0
soln = [y0]
while solver.successful() and solver.t < t[-1]:
    k += 1
    solver.integrate(t[k])
    soln.append(solver.y)

# Convert the list to a numpy array.
soln = np.array(soln)

接下来这部分内容会讲讲如何继续使用odeint。

如果你只对线性插值感兴趣，可以简单地用数据的最后两个点进行线性延伸。扩展data数组的一个简单方法是把值2*data[-1] - data[-2]添加到数组的末尾，t数组也可以这样做。如果t中的最后一个时间步长很小，这样的延伸可能不够长，无法避免问题，因此在下面我使用了更通用的延伸方法。

示例：

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

t = np.linspace(0, 3, 4)

data = [1, 2, 3, 4]

# Slope of the last segment.
m = (data[-1] - data[-2]) / (t[-1] - t[-2])
# Amount of time by which to extend the interpolation.
dt = 3.0
# Extended final time.
t_ext = t[-1] + dt
# Extended final data value.
data_ext = data[-1] + m*dt

# Extended arrays.
extended_t = np.append(t, t_ext)
extended_data = np.append(data, data_ext)

linear_interpolation = interp1d(extended_t, extended_data)

def func(y, t0):
    print 't0', t0
    return -2*y + linear_interpolation(t0)

soln = odeint(func, 1, t)

如果仅仅用最后两个数据点来线性延伸插值线太粗糙，那么你就需要使用其他方法来在给odeint的最后t值之外进行外推。

另一种选择是将最后的t值作为参数传递给func，并在func中明确处理大于这个值的t值。可以像这样做，其中extrapolation是你需要自己搞清楚的内容：

def func(y, t0, tmax):
    if t0 > tmax:
        f = -2*y + extrapolation(t0)
    else:
        f = -2*y + linear_interpolation(t0)
    return f

soln = odeint(func, 1, t, args=(t[-1],))