4

作为 MATLAB 的转换者，我对 python 和 scipy 比较陌生。我正在对 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

我制作了一个图表来说明两种结果的差异，请单击此处。

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

4

1 回答 1

3

这不是错误。问题在于您已转向bound_errorFalse 并用零填充这些值。如果您bound_error在原始代码中设置为 True ，您可以看到您超出了插值的范围，因此在积分中输入了零（因此产生的值与您在范围就像你对 lambda 所做的那样x_1。

尝试以下操作，您会发现一切正常。基本上，我只是扩展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,) )
于 2011-04-13T14:14:26.293 回答