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
On Thu, Jan 28, 2010 at 3:39 PM, Ryan Krauss <ryanlists@gmail.com> 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
When I add a small noise u2 = zeros(N) + 1e-14 or for u2[:50] = amp or for u2[50:200] = amp it seems to work. This might be a tricky bug. Josef
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Hmmm. Thanks. That solves the immediate problem. I am letting my students choose between Matlab and Python for projects in my course. This one my erode their confidence in Python/Scipy a bit. On Thu, Jan 28, 2010 at 3:05 PM, <josef.pktd@gmail.com> wrote:
On Thu, Jan 28, 2010 at 3:39 PM, Ryan Krauss <ryanlists@gmail.com> 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
When I add a small noise
u2 = zeros(N) + 1e-14
or for u2[:50] = amp
or for u2[50:200] = amp
it seems to work.
This might be a tricky bug.
Josef
_______________________________________________ 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
On Thu, Jan 28, 2010 at 6:00 PM, Ryan Krauss <ryanlists@gmail.com> wrote:
Hmmm. Thanks. That solves the immediate problem. I am letting my students choose between Matlab and Python for projects in my course. This one my erode their confidence in Python/Scipy a bit.
It's always good to have a rough idea about whether the results are correct and not trust the computer too much, whether it's matlab or scipy. But if some of your students are willing to submit bug reports or tests, then the next generation of students can be more confident that it's the bugs in their own program that might be causing problems and not the code in scipy. Josef
On Thu, Jan 28, 2010 at 3:05 PM, <josef.pktd@gmail.com> wrote:
On Thu, Jan 28, 2010 at 3:39 PM, Ryan Krauss <ryanlists@gmail.com> 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
When I add a small noise
u2 = zeros(N) + 1e-14
or for u2[:50] = amp
or for u2[50:200] = amp
it seems to work.
This might be a tricky bug.
Josef
_______________________________________________ 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
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.) 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()
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
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
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
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
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@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
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
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
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
On Tue, Feb 9, 2010 at 1: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.
Thanks, Stefan (I think) or I will look at it. For statsmodels (maximum likelihood) I was looking at a second pattern instead of **kw that keeps the options for the optimizer, or ode in this case, separate. It's just a taste difference. **kw is the more usual pattern in scipy. Using a separate keyword for the options instead, keeps the signature a bit cleaner (and I have been reading more gauss and matlab lately), but it requires slightly more typing in the call. Josef def fun1a(a, opt1=3, **kw): #options to second function as keywords print '\nfun1a a',a print 'fun1a opt1', opt1 print 'fun1a kw', kw fun2(a, **kw) return 'end fun1a' def fun1b(a, opt1=3, optfun={}): #options to second function as separate dictionary print '\nfun1b a',a print 'fun1 opt1', opt1 print 'fun1b optfun', optfun fun2(a, **optfun) return 'end fun1b' def fun2(a, tol=0, maxiter=1): print 'fun2 a',a print 'fun2 kw: tol, maxiter', tol, maxiter return 'end fun2' print fun1a('1a') print fun1a('1a', tol=5, maxiter=10) print fun1a('1a', opt1=4, tol=5, maxiter=10) my_options = dict(tol=5, maxiter=10) print fun1a('1a', opt1=4, **my_options) print fun1a('1b') print fun1b('1b', optfun=dict(tol=5, maxiter=10)) print fun1b('1b', opt1=4, optfun={ 'tol':5, 'maxiter':10}) my_options = dict(tol=5, maxiter=10) print fun1b('1b', opt1=4, optfun=my_options)
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
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
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
On Sun, Feb 14, 2010 at 8:02 PM, Warren Weckesser <warren.weckesser@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@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
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
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@gmail.com> wrote:
On Sun, Feb 14, 2010 at 8:02 PM, Warren Weckesser <warren.weckesser@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@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
_______________________________________________ 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
josef.pktd@gmail.com wrote:
On Sun, Feb 14, 2010 at 8:02 PM, Warren Weckesser <warren.weckesser@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'd like to look into this, too. I don't think it was originally designed for MIMO systems.
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.
I suspect the default precision is fine. To actually know the expected precision, one has to be able to convert the absolute and relative tolerances used by the ODE solver, which are *local* error control parameters, into a global tolerance (i.e. what is the error at t=f_final). This is a nontrivial problem. Warren
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@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
_______________________________________________ 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
josef.pktd@gmail.com wrote:
On Sun, Feb 14, 2010 at 8:02 PM, Warren Weckesser <warren.weckesser@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'd like to look into this, too. I don't think it was originally designed for MIMO systems.
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.
I suspect the default precision is fine. To actually know the expected precision, one has to be able to convert the absolute and relative tolerances used by the ODE solver, which are *local* error control parameters, into a global tolerance (i.e. what is the error at t=f_final). This is a nontrivial problem. Warren
Or, you could modify your script lsim2_problem.py to eliminate the initial interval during which the solution is exactly zero, by changing this line: u2[50:100] = amp to this: u2[:50] = amp This just translates the pulse so that it starts at t=0. Since the default initial condition used by lsim2 is zero, it will give the same solution, just translated in time. Warren Warren Weckesser 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.)
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
------------------------------------------------------------------------
------------------------------------------------------------------------
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
participants (3)
-
josef.pktd＠gmail.com -
Ryan Krauss -
Warren Weckesser