[SciPy-User] bug in signal.lsim2

Ryan Krauss ryanlists at gmail.com
Wed Feb 17 08:39:48 EST 2010 

I like it alot.  It goes way beyond solving my problem to adding new
functionality, improving the documentation, and creating some well
thought out tests.  Nice work Warren.

I didn't see how to download the original version of the patch that
created test_ltisys.py (maybe it should already exist on my computer
if I had a more recent version of scipy).  So, patch yelled at me
about not finding the file.  So, I didn't actually run Warren's tests;
I just ran a couple of simple cases I made up.  But it worked
perfectly.

Thanks again,

Ryan

On Sun, Feb 14, 2010 at 8:43 PM,  <josef.pktd at gmail.com> wrote:
> On Sun, Feb 14, 2010 at 8:02 PM, Warren Weckesser
> <warren.weckesser at enthought.com> wrote:
>> I added a patch to ticket #1112.  In addition to adding the suggested
>> change to pass keyword args to odeint, I also added a new funciton,
>> impulse2(), that computes the impulse response by using odeint.  See
>> this thread for details:
>>
>>    http://mail.scipy.org/pipermail/scipy-user/2009-November/023416.html
>
> Thanks, the test examples look nice and serve also as some
> documentation for ltisys. The lsim2 test05 is interesting, I was
> convinced last year (and there were some previous threads) that ltisys
> cannot handle multi-input systems. I always got shape errors when I
> tried and without documentation it's difficult to figure out what is
> supposed to work and what not. (I will look at it more closely when I
> have more time.)
>
> I think the changes are reasonable including turning u and t into
> keywords instead of required arguments. Maybe Ryan or someone who is
> using this functions can comment on this.
>
> Just one improvement to the tests, assert_almost_equal takes a decimal
> argument. If you know the precision of your tests, then it would be
> useful to increase it from the default decimal=7.
>
> Josef
>
>
>>
>> Warren
>>
>>
>> Ryan Krauss wrote:
>>> Yeah, **kwds is a little harder to understand and slightly less user
>>> friendly, but it is easier to maintain and less work for the patch
>>> writer.  And it is more future proof if the integrator ever changes.
>>>
>>> On Tue, Feb 9, 2010 at 12:17 PM, Warren Weckesser
>>> <warren.weckesser at enthought.com> wrote:
>>>
>>>> josef.pktd at gmail.com wrote:
>>>>
>>>>> On Wed, Feb 3, 2010 at 9:00 AM, Ryan Krauss <ryanlists at gmail.com> wrote:
>>>>>
>>>>>
>>>>>> FYI, I am quite happy with passing in an hmax value.  I basically
>>>>>> copied and pasted lsim2 from signal.ltisys and adapted it just a
>>>>>> little to make it a method of my derived class.  Then I added the hmas
>>>>>> kwarg that gets passed to odeint.
>>>>>>
>>>>>> Is there any reason not to allow the user to pass in a kwargs to lsim2
>>>>>> that gets passed to odeint?
>>>>>>
>>>>>>
>>>>> I don't see a reason why we cannot add a **kwargs, it should be
>>>>> completely backwards compatible.
>>>>> Can you file a ticket and add your adjusted version or a patch? And
>>>>> even better, add your original example as a test case?
>>>>>
>>>>>
>>>>>
>>>> Josef,
>>>>
>>>> I just created ticket #1112 for this.  Unless Ryan wants to adapt his
>>>> change to lsim2, I can make a patch this week to implement the enhancement.
>>>>
>>>> Warren
>>>>
>>>>
>>>>
>>>>> Josef
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>> On Fri, Jan 29, 2010 at 6:44 AM, Ryan Krauss <ryanlists at gmail.com> wrote:
>>>>>>
>>>>>>
>>>>>>> Thanks to Warren and Josef for their time and thoughts.  I feel like I
>>>>>>> now understand the underlying problem and have some good options to
>>>>>>> solve my short term issues (I assigned the project last night and they
>>>>>>> need to be able to start working on it immediately).  I actually use a
>>>>>>> TransferFunction class that derives from ltisys.  I could override its
>>>>>>> lsim2 method to try out some of these solutions quickly and fairly
>>>>>>> easily.
>>>>>>>
>>>>>>> Ryan
>>>>>>>
>>>>>>> On Thu, Jan 28, 2010 at 10:00 PM,  <josef.pktd at gmail.com> wrote:
>>>>>>>
>>>>>>>
>>>>>>>> On Thu, Jan 28, 2010 at 10:33 PM, Warren Weckesser
>>>>>>>> <warren.weckesser at enthought.com> wrote:
>>>>>>>>
>>>>>>>>
>>>>>>>>> josef.pktd at gmail.com wrote:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>> On Thu, Jan 28, 2010 at 8:50 PM, Warren Weckesser
>>>>>>>>>> <warren.weckesser at enthought.com> wrote:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> 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.)
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>> It's something what I suspected. I don't know much about odeint, but
>>>>>>>>>> do you think it would be useful to let lsim2 pass through some
>>>>>>>>>> parameters to odeint?
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>> Sounds useful to me.  A simple implementation is an optional keyword
>>>>>>>>> argument that is a dict of odeint arguments.   But this would almost
>>>>>>>>> certainly break if lsim2 were ever reimplemented with a different
>>>>>>>>> solver.  So perhaps it should allow a common set of ODE solver
>>>>>>>>> parameters (e.g. absolute and relative error tolerances, max and min
>>>>>>>>> step sizes, others?).
>>>>>>>>>
>>>>>>>>> Perhaps this should wait until after the ODE solver redesign that is
>>>>>>>>> occasionally discussed:
>>>>>>>>>    http://projects.scipy.org/scipy/wiki/OdeintRedesign
>>>>>>>>> Then the solver itself could be an optional argument to lsim2.
>>>>>>>>>
>>>>>>>>>
>>>>>>>> I was just thinking of adding to the argument list a **kwds argument
>>>>>>>> that is directly passed on to whatever ODE solver is used. This should
>>>>>>>> be pretty flexible for any changes and be backwards compatible.
>>>>>>>>
>>>>>>>> I've seen and used it in a similar way for calls to optimization
>>>>>>>> routines, e.g. also optimize.curve_fit, does it. What are actually
>>>>>>>> valid keywords would depend on which function is called.
>>>>>>>>
>>>>>>>> (But I'm not a user of lsim, I'm just stealing some ideas from lti and
>>>>>>>> friends for time series analysis.)
>>>>>>>>
>>>>>>>> Josef
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>> Warren
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>> Josef
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> 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
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>> from pylab import *
>>>>>>>>>>> from scipy import signal
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> g = 100.0
>>>>>>>>>>> p = 15.0
>>>>>>>>>>> G = signal.ltisys.lti(g, [1,p,0])
>>>>>>>>>>>
>>>>>>>>>>> t = arange(0, 1.0, 0.002)
>>>>>>>>>>> N = len(t)
>>>>>>>>>>>
>>>>>>>>>>> # u for the whole interval (not used in lsim2, only for plotting later).
>>>>>>>>>>> amp = 50.0
>>>>>>>>>>> u = zeros(N)
>>>>>>>>>>> k1 = 50
>>>>>>>>>>> k2 = 100
>>>>>>>>>>> u[k1:k2] = amp
>>>>>>>>>>>
>>>>>>>>>>> # Create input functions for each smooth interval. (This could be simpler,
>>>>>>>>>>> since u
>>>>>>>>>>> # is constant on each interval.)
>>>>>>>>>>> a = float(k1)/N
>>>>>>>>>>> b = float(k2)/N
>>>>>>>>>>> T1 = linspace(0, a, 201)
>>>>>>>>>>> u1 = zeros_like(T1)
>>>>>>>>>>> T2 = linspace(a, b, 201)
>>>>>>>>>>> u2 = amp*ones_like(T2)
>>>>>>>>>>> T3 = linspace(b, 1.0, 201)
>>>>>>>>>>> u3 = zeros_like(T3)
>>>>>>>>>>>
>>>>>>>>>>> # Solve on each interval; use the final value of one solution as the
>>>>>>>>>>> starting
>>>>>>>>>>> # point of the next solution.
>>>>>>>>>>> # (We could skip the first calculation, since we know the solution will be
>>>>>>>>>>> 0.)
>>>>>>>>>>> (t1, y1, x1) = signal.lsim2(G,u1,T1)
>>>>>>>>>>> (t2, y2, x2) = signal.lsim2(G, u2, T2, X0=x1[-1])
>>>>>>>>>>> (t3, y3, x3) = signal.lsim2(G, u3, T3, X0=x2[-1])
>>>>>>>>>>>
>>>>>>>>>>> figure(1)
>>>>>>>>>>> clf()
>>>>>>>>>>> plot(t, u, 'k', linewidth=3)
>>>>>>>>>>> plot(t1, y1, 'y', linewidth=3)
>>>>>>>>>>> plot(t2, y2, 'b', linewidth=3)
>>>>>>>>>>> plot(t3, y3, 'g', linewidth=3)
>>>>>>>>>>>
>>>>>>>>>>> show()
>>>>>>>>>>>
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> SciPy-User mailing list
>>>>>>>>>>> SciPy-User at scipy.org
>>>>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>> _______________________________________________
>>>>>>>>>> SciPy-User mailing list
>>>>>>>>>> SciPy-User at scipy.org
>>>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> SciPy-User mailing list
>>>>>>>>> SciPy-User at scipy.org
>>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> SciPy-User mailing list
>>>>>>>> SciPy-User at scipy.org
>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>> _______________________________________________
>>>>>> SciPy-User mailing list
>>>>>> SciPy-User at scipy.org
>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>
>>>>>>
>>>>>>
>>>>> _______________________________________________
>>>>> SciPy-User mailing list
>>>>> SciPy-User at scipy.org
>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>
>>>>>
>>>> _______________________________________________
>>>> SciPy-User mailing list
>>>> SciPy-User at scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>
>>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>

More information about the SciPy-User mailing list