Возможная ошибка в odeint <-> interp1d interplay?

Question

Я относительно новичок в python и scipy, будучи конвертированным из MATLAB.Я делал быструю проверку функции odeint в scipy.integrate и наткнулся на эту потенциальную ошибку.Рассмотрим следующий фрагмент:

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, вы можете увидеть, что вы выходите за пределы интерполяции и, таким образом, ставите в интеграцию нули (и, таким образом, получаете другое значение, чем если бы вы оценивали функцию в этих точках вне диапазона, как вы делаете с лямбдой для 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,) )