[SciPy-User] bug in signal.lsim2
Warren Weckesser
warren.weckesser at enthought.com
Thu Jan 28 20:50:05 EST 2010
Ryan,
The problem is that the ODE solver used by lsim2 is too good. :)
It uses scipy.integrate.odeint, which in turn uses the Fortran library
LSODA. Like any good solver, LSODA is an adaptive solver--it adjusts
its step size to be as large as possible while keeping estimates of the
error bounded. For the problem you are solving, with initial condition
0, the exact solution is initially exactly 0. This is such a nice
smooth solution that the solver's step size quickly grows--so big, in
fact, that it skips right over your pulse and never sees it.
So how does it create all those intermediate points at the requested
time values? It uses interpolation between the steps that it computed
to create the solution values at the times that you requested. So using
a finer grid of time values won't help. (If lsim2 gave you a hook into
the parameters passed to odeint, you could set odeint's 'hmax' to a
value smaller than your pulse width, which would force the solver to see
the pulse. But there is no way to set that parameter from lsim2.)
The basic problem is you are passing in a discontinuous function to a
solver that expects a smooth function. A better way to solve this
problem is to explicitly account for the discontinuity. One possibility
is the attached script.
This is an excellent "learning opportunity" for your students on the
hazards of numerical computing!
Warren
Ryan Krauss wrote:
> I believe I have discovered a bug in signal.lsim2. I believe the
> short attached script illustrates the problem. I was trying to
> predict the response of a transfer function with a pure integrator:
>
> g
> G = -------------
> s(s+p)
>
> to a finite width pulse. lsim2 seems to handle the step response just
> fine, but says that the pulse response is exactly 0.0 for the entire
> time of the simulation. Obviously, this isn't the right answer.
>
> I am running scipy 0.7.0 and numpy 1.2.1 on Ubuntu 9.04, but I also
> have the same problem on Windows running 0.7.1 and 1.4.0.
>
> Thanks,
>
> Ryan
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: lsim2_solution.py
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100128/22baf859/attachment.ksh>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lsim2_solution.png
Type: image/png
Size: 13212 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100128/22baf859/attachment.png>
More information about the SciPy-User
mailing list