python - odeint <-> interp1d 相互作用中可能存在的错误？

Question

我是 Python 和 scipy 的新手，是从 MATLAB 转过来的。我正在快速测试 scipy.integrate 中的 odeint 函数，并遇到了这个潜在的错误。考虑以下代码段:

from scipy import *
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from pylab import *

# ODE system with forcing function u(t)
def sis(x,t,u):
    return [x[1], u(t)]

# Solution time span
t = linspace(0, 10, 1e3)

# Profile for forcing function u(t)
acel = lambda t: 3*(t<2)-3*(t>8)

# Interpolator for acelerator
acel_interp = interp1d(t, acel(t), bounds_error=False, fill_value=0)    

# ODE integration with u(t) = acel, a lambda function
x_1 = odeint(sis, [0,0], t, args=(acel,) )            # Correct result
# ODE integration with u(t) = acel_interp, an interpolator
x_2 = odeint(sis, [0,0], t, args=(acel_interp,) )     # Incorrect result

我画了一个图来说明这两个结果的差异，click here .

至少对我而言，您如何看待这种毫无根据的结果差异？我在 Python 2.6.6 之上使用 NumPy 1.5.0 版和 SciPy 0.8.0 版

Answer 1

这不是错误。问题在于您已将 bound_error 变为 False 并用零填充这些值。如果您在原始代码中将 bound_error 设置为 True，您可以看到您超出了插值的界限，因此在积分中置入了零(因此产生了与评估时不同的值函数在范围之外的那些点上，就像您对 x_1 的 lambda 所做的那样)。

尝试以下操作，您会发现一切正常。基本上，我只是扩展了 t 以覆盖足够大的值范围，以覆盖您正在使用插值的范围。

from scipy import *
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from pylab import *

# ODE system with forcing function u(t)
def sis(x,t,u):
    return [x[1], u(t)]

# Solution time span
t = linspace(0, 10, 1e3)
t_interp = linspace(0,20,2e3)

# Profile for forcing function u(t)
acel = lambda t: 3*(t<2)-3*(t>8)

# Interpolator for acelerator
acel_interp = interp1d(t_interp, acel(t_interp))    

# ODE integration with u(t) = acel, a lambda function
x_1 = odeint(sis, [0,0], t, args=(acel,) )            
# ODE integration with u(t) = acel_interp, an interpolator
x_2 = odeint(sis, [0,0], t, args=(acel_interp,) )