[SciPy-User] bug in signal.lsim2
Warren Weckesser
warren.weckesser@enthought....
Thu Jan 28 21:33:09 CST 2010
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.
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
>
More information about the SciPy-User
mailing list