[SciPy-User] bug in signal.lsim2

Ryan Krauss ryanlists at gmail.com
Wed Feb 3 09:00:36 EST 2010 

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?

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

More information about the SciPy-User mailing list