Integrated.ode将t0值设置为超出我的数据范围 [英] integrate.ode sets t0 values outside of my data range

问题描述

对于y(t = 0)= 1，我想在t = 0..3之间求解ODE dy/dt = -2y + data(t).

I would like to solve the ODE dy/dt = -2y + data(t), between 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)

当我运行这段代码时，我得到了一些类似这样的错误:

When I run this code, I get several errors like this:

ValueError:x_new中的值超出插值范围.
odepack.error:调用名为Python的函数时发生错误 功能

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的值时，我意识到t0有时实际上超出了我的插值范围:3.07634612585、3.0203768998、3.000638459329，...这就是linear_interpolation(t0)引发ValueError异常的原因.

My interpolation range is between 0.0 and 3.0. Printing the value of t0 in func, I realized that t0 is actually sometimes above my interpolation range: 3.07634612585, 3.0203768998, 3.00638459329, ... It's why linear_interpolation(t0) raises the ValueError exceptions.

我有几个问题:

• integrate.ode如何使t0变化?为什么t0超出插值范围的最小值(3.0)?

• how does integrate.ode makes t0 vary? Why does it make t0 exceed the infimum (3.0) of my interpolation range?

仍会返回一个似乎包含正确值的数组.因此，我是否应该捕获并忽略这些错误?无论微分方程，t范围和初始条件如何，我都应该忽略它们吗?

in spite of these errors, integrate.ode returns an array which seems to contain correct value. So, should I just catch and ignore these errors? Should I ignore them whatever the differential equation(s), the t range and the initial condition(s)?

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

if I shouldn't ignore these errors, what is the best way to avoid them? 2 suggestions:

在interp1d中的

• ，我可以设置bounds_error=False和fill_value=data[-1]，因为插值范围之外的t0似乎与t[-1]接近:

• in interp1d, I could set bounds_error=False and fill_value=data[-1] since the t0 outside of my interpolation range seem to be closed to 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有所不同将有助于我确定这一点(请参阅第一个问题).

But first I would like to be sure that with any other func and any other data the t0 will always remain closed to t[-1]. For example, if integrate.ode chooses a t0 below my interpolation range, the fill_value would be still data[-1], which would not be correct. Maybe to know how integrate.ode makes t0 vary would help me to be sure of that (see my first question).

中，我可以将linear_interpolation调用包含在try/except块中，并且当我捕获到ValueError时，我会回想起linear_interpolation，但将t0截断:

in func, I could enclose the linear_interpolation call in a try/except block, and, when I catch a ValueError, I recall linear_interpolation but with t0 truncated:

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，则至少此解决方案允许linear_interpolation仍然引发异常.然后，我会收到有关不连贯行为的警报.但这不是真正可读的，到现在，截断在我看来似乎有些武断.

At least this solution permits linear_interpolation to still raise an exception if integrate.ode makes t0 >= 4.0 or t0 <= -1.0. I can then be alerted of incoherent behavior. But it is not really readable and the truncation seems to me a little arbitrary by now.

也许我只是对这些错误有过多的考虑.请让我知道.

Maybe I'm just over-thinking about these errors. Please let me know.

推荐答案

odeint求解器在超过最后请求时间的时间值上评估函数是正常的.大多数ODE求解器都是以这种方式工作的-他们采用内部时间步长，其大小由其错误控制算法确定，然后使用自己的插值法在用户要求的时间评估解决方案.某些求解器(例如Sundials库中的CVODE求解器)允许您指定时间的硬限制，超过该时间则不允许求解器评估您的方程式，但是odeint没有这种选择.

It is normal for the odeint solver to evaluate your function at time values past the last requested time. Most ODE solvers work this way--they take internal time steps with sizes determined by their error control algorithm, and then use their own interpolation to evaluate the solution at the times requested by the user. Some solvers (e.g. the CVODE solver in the Sundials library) allow you to specify a hard bound on the time, beyond which the solver is not allowed to evaluate your equations, but odeint does not have such an option.

如果您不介意从scipy.integrate.odeint切换到scipy.integrate.ode，则看起来"dopri5"和"dop853"求解器不会在超出请求时间的时间内评估您的功能.两个警告:

If you don't mind switching from scipy.integrate.odeint to scipy.integrate.ode, it looks like the "dopri5" and "dop853" solvers don't evaluate your function at times beyond the requested time. Two caveats:

• ode求解器对定义微分方程的参数顺序使用不同的约定.在ode求解器中，t是第一个参数. (是的，我知道，抱怨，抱怨...)
• "dopri5"和"dop853"求解器适用于非刚性系统.如果您的问题是僵硬，他们仍然应该给出正确的答案，但是他们会做更多的工作比僵硬的求解器要好.

• The ode solvers use a different convention for the order of the arguments that define the differential equation. In the ode solvers, t is the first argument. (Yeah, I know, grumble, grumble...)
• The "dopri5" and "dop853" solvers are for non-stiff systems. If your problem is stiff, they should still give correct answers, but they will do a lot more work than a stiff solver would do.

这是一个脚本，显示如何解决您的示例.为了强调参数的变化，我将func重命名为rhs.

Here's a script that shows how to solve your example. To emphasize the change in the arguments, I renamed func to 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.

The rest of this answer looks at how you could continue to use odeint.

如果您只对线性插值感兴趣，则可以使用数据的最后两个点简单地线性扩展数据.扩展data数组的一种简单方法是将值2*data[-1] - data[-2]附加到数组的末尾，并对t数组执行相同的操作.如果t中的最后一个时间步很小，则该扩展名可能不会足够长以避免发生该问题，因此在下面，我将使用更通用的扩展名.

If you are only interested in linear interpolation, you could simply extend your data linearly, using the last two points of the data. A simple way to extend the data array is to append the value 2*data[-1] - data[-2] to the end of the array, and do the same for the t array. If the last time step in t is small, this might not be a sufficiently long extension to avoid the problem, so in the following, I've used a more general extension.

示例:

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值.

If simply using the last two data points to extend the interpolator linearly is too crude, then you'll have to use some other method to extrapolate a little beyond the final t value given to odeint.

另一种选择是将最终的t值包含为func的参数，并显式处理大于func中的t值.像这样的东西，其中extrapolation是您必须弄清楚的:

Another alternative is to include the final t value as an argument to func, and explicitly handle t values larger than it in func. Something like this, where extrapolation is something you'll have to figure out:

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],))

这篇关于Integrated.ode将t0值设置为超出我的数据范围的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！