集成.ode设置t0值超出我的数据范围

Question

我想求出d y/dt=-2y+data（t），在t=0..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: A value in x_new is above the interpolation range.
odepack.error: Error occurred while calling the Python function named func

我的插值范围在0.0到3.0之间。 在func中打印t0的值时，我意识到{}实际上有时高于我的插值范围：3.07634612585，3.0203768998，3.00638459329。。。这就是linear_interpolation(t0)引发ValueError异常的原因。在

我有几个问题：

• integrate.ode如何使t0变化？为什么它使t0超过了插值范围的下确界（3.0）？

• 尽管有这些错误，integrate.ode仍返回一个似乎包含正确值的数组。那么，我是不是应该抓住并忽略这些错误呢？不管微分方程，范围和初始条件，我应该忽略它们吗？

• 如果我不应该忽视这些错误，那么最好的避免方法是什么？2建议：

  • 在interp1d中，我可以设置bounds_error=False和{}，因为插值范围之外的t0似乎接近{}：

    ^{pr2}$

    但是首先，我想确定的是，对于任何其他的func和任何其他{}，t0将始终与{}保持闭合。例如，如果integrate.ode在我的插值范围以下选择一个t0，填充值仍然是data[-1]，这是不正确的。也许知道integrate.ode是如何使t0变化的，这将有助于我确定这一点（见我的第一个问题）。

  • 在func中，我可以将linear_interpolation调用括在try/except块中，当我捕捉到一个ValueError时，我会回忆起linear_interpolation，但是{}被截断了：

    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
    

    至少，如果integrate.ode使t0>；=4.0或t0<；=-1.0时，线性插值仍然会引发异常。然后我就可以被提醒不连贯的行为。但它不是真正的可读性和截断似乎有点武断现在。

也许我只是想了想这些错误。请告诉我。在

Answer 1

odeint解算器在超过上次请求时间的时间值时计算函数是正常的。大多数ODE解算器都是这样工作的，它们采用内部时间步长，其大小由误差控制算法确定，然后在用户要求的时间内使用自己的插值来评估解。有些解算器（例如日晷库中的CVODE解算器）允许您指定时间的硬边界，超过该范围，解算器将不允许计算您的方程，但odeint没有这样的选项。在

如果您不介意从scipy.integrate.odeint切换到scipy.integrate.ode，那么"dopri5"和{}解算器在请求的时间之外不会计算您的函数。两个注意事项：

• ode解算器对定义微分方程的参数顺序使用不同的约定。在ode解算器中，t是第一个参数。（是的，我知道，发牢骚，发牢骚…）
• "dopri5"和{}解算器适用于非刚性系统。如果你的问题是stiff，他们仍然应该给出正确的答案，但是他们所做的工作将比僵硬的解算器做的多得多。在

下面是一个脚本，它演示了如何求解示例。为了强调参数的变化，我将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中的最后一个时间步很小，那么这可能不是一个足够长的扩展来避免这个问题，所以在下面，我使用了一个更通用的扩展。在

示例：

^{pr2}$

如果简单地使用最后两个数据点来线性扩展插值器太粗糙了，那么您将不得不使用其他方法来稍微超出赋予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],))