Возможная ошибка во взаимодействии odeint ‹-› interp1d?

Я относительно новичок в 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.


python scipy interpolation ode numerical-integration
person Luis E.    schedule 13.04.2011    source источник
comment
Поздравляю с переходом с Matlab. Я сделал шаг около 18 месяцев назад. Замечательно иметь доступ ко всем другим возможностям Python и иметь возможность делиться своим кодом с людьми, которые не заплатили тысячи за лицензию Matlab.   -  person Carl F.    schedule 17.04.2011

Ответы (1)


arrow_upward
3
arrow_downward

Это не ошибка. Проблема в том, что вы изменили 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,) )
person JoshAdel    schedule 13.04.2011
comment
Вы опередили меня на 13 минут :) Когда я ловлю первое исключение, созданное u(t) при оценке sis(), значение для t равно 17.8811987021. Тем не менее, odeint() велят интегрировать от 0 до 10. Как это работает? - person lafras; 13.04.2011
comment
Посмотрите на x_2,idct = odeint(sis, [0,0], t, args=(acel_interp,) ,full_output=True); idct['tcur']. Какая бы схема интеграции ни использовалась, она требует оценок за пределами диапазона t. Как указано в документации: 'tcur' - вектор со значением t, достигнутым для каждого временного шага. (всегда будет не меньше времени ввода). - person JoshAdel; 13.04.2011
comment
Спасибо за вашу мудрость! Я принял ваш ответ, так как он работает правильно. Однако я до сих пор не понимаю, почему он не работал с исходным кодом, поскольку начальный интервал времени (массив t) включает все соответствующие изменения в acel(t) (последний переход происходит в t = 8, а интервал времени установлен равным т = 10). Не могли бы вы подробнее рассказать об этом? - person Luis E.; 15.04.2011
comment
@Luis: вам нужно найти реализацию odeint, основанную на алгоритме lsoda из ODEPACK, чтобы точно ответить на этот вопрос. - person JoshAdel; 15.04.2011