[SciPy-User] bug in signal.lsim2

josef.pktd@gmai... josef.pktd@gmai...
Thu Jan 28 19:57:31 CST 2010


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?

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


More information about the SciPy-User mailing list