odeint与interp1d的交互中可能存在错误？

Question

我刚开始学习Python和SciPy，之前是用MATLAB的。我在测试SciPy库里的odeint函数时，发现了一个可能的bug。看看下面这段代码：

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

我做了一个图，展示了这两个结果的差异，点这里查看.

你觉得这个结果之间的差异有什么问题吗？对我来说，这种差异似乎没有理由。我使用的是NumPy 1.5.0版本和SciPy 0.8.0版本，Python版本是2.6.6。

Answer 1

这不是一个错误。问题在于你把 bound_error 设置成了 False，并且用零填充了那些值。如果你在原来的代码中把 bound_error 设置为 True，你会发现你超出了插值的范围，因此在积分时填入了零（这样得到的值就和你在那些超出范围的点上用 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,) )