积分.Ode在我的数据范围之外设置了0个值

Question

我想解ODE dy/dt = -2y + data(t)， t=0。3、for 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:调用Python函数时发生错误func

我的插值范围在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范围和初始条件如何，我都应该忽略它们吗?

• 如果我不应该忽略这些错误，那么避免它们的最佳方法是什么?2 .建议:

  • 在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低于我的插值范围，fill_value仍然是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
    

    至少这个解决方案允许线性插值仍然引发异常，如果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],))