[SciPy-User] bug in signal.lsim2
Warren Weckesser
warren.weckesser@enthought....
Sun Feb 14 19:02:53 CST 2010
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
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@enthought.com> wrote:
>
>> josef.pktd@gmail.com wrote:
>>
>>> On Wed, Feb 3, 2010 at 9:00 AM, Ryan Krauss <ryanlists@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@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@gmail.com> wrote:
>>>>>
>>>>>
>>>>>> On Thu, Jan 28, 2010 at 10:33 PM, Warren Weckesser
>>>>>> <warren.weckesser@enthought.com> wrote:
>>>>>>
>>>>>>
>>>>>>> josef.pktd@gmail.com wrote:
>>>>>>>
>>>>>>>
>>>>>>>> On Thu, Jan 28, 2010 at 8:50 PM, Warren Weckesser
>>>>>>>> <warren.weckesser@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@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@scipy.org
>>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> SciPy-User mailing list
>>>>>>>> SciPy-User@scipy.org
>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> SciPy-User mailing list
>>>>>>> SciPy-User@scipy.org
>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>> _______________________________________________
>>>>>> SciPy-User mailing list
>>>>>> SciPy-User@scipy.org
>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>
>>>>>>
>>>>>>
>>>> _______________________________________________
>>>> SciPy-User mailing list
>>>> SciPy-User@scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>
>>>>
>>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list