Possible error in odeint & # 8596; interp1d interplay?

Question

I am relatively new to python and scipy, being a converter from MATLAB. I quickly checked the odeint function in scipy.integrate and came across this potential error. Consider the following snippet:

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 

I made a conspiracy that illustrates the difference between the two results, click here .

What do you make of this, at least for me, an unjustified difference in the results? I am using NumPy version 1.5.0 and SciPy version 0.8.0 on top of Python 2.6.6

Answer 1

It's not a mistake. The problem is that you turned bound_error to False and filled these values ​​with zeros. If you set bound_error to True in the source code, you will see that you exceed the boundaries of your interpolation and thus insert zeros into the integration (and thus get a different value than if you evaluated the function at these points out of range, as for lambda for x_1 ).

Try the following and you will see that everything is working correctly. Basically, I just expanded t to cover a range of values ​​large enough to cover the range in which you use interpolation.

 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,) )