Impulse/Step response troubles
Dear SciPy users, (Long mail) I am having several problems in working out the impulse response of a simple function in SciPy. I would request you to please point out whether the fault lies with SciPy or me. The question is to determine (plot) the step and ramp responses of K / (s^2 + 8s + K) for K = 7, 16 and 80 So, the following code is used, and I'll keep changing a line: <code> from scipy import * from pylab import * K = 16.0 # should be redone with 7.0, 16.0 and 80.0 r = signal.impulse(([K],[1.0,8.0,K]), T=r_[0:5:0.001]) plot(r[0], r[1], linewidth=2) show() </code> The above code plots the impulse response of the given function. Now, on to the action: 1. Running the above code with K = 7.0 and K = 16.0 gives expected results. However, running the code with K = 16.0 doesn't work; if K = 16.0, the response should be 16 * t * exp(-4t) * u(t), which it isn't. Changing K to 16 + or - .00000000001 fixes it. There surely is some problem when the value of K is such that a double root is hit. 2. For the step response, I change "impulse" in the above code to "step", things seem all right. Again, for K = 16, it doesn't work. I then switch back to impulse, but make the second (denominator) array [1.0, 8.0, K, 0], and still the same thing is observed. 3. To get the ramp response, I now have the following code in line 8 of the above code. r = signal.step(([K],[1.0,8.0,K, 0]), T=r_[0:5:0.001]) But this causes the error: Traceback (most recent call last): File "<stdin>", line 8, in ? File "/usr/lib/python2.4/site-packages/scipy/signal/ltisys.py", line 498, in step vals = lsim(sys, U, T, X0=X0) File "/usr/lib/python2.4/site-packages/scipy/signal/ltisys.py", line 403, in lsim ATm1 = linalg.inv(AT) File "/usr/lib/python2.4/site-packages/scipy/linalg/basic.py", line 306, in inv if info>0: raise LinAlgError, "singular matrix" numpy.linalg.linalg.LinAlgError: singular matrix So let's switch back to impulse: r = signal.impulse(([K],[1.0,8.0,K, 0.0, 0.0]), T=r_[0:5:0.001]) Now, it seems to work, but I am yet to verify the correctness of the plots. Naturally, for K = 16.0, it doesn't work; I have to use a value close to it. I would really appreciate help in getting this debugged. Thanks a lot! Kumar -- Kumar Appaiah, 458, Jamuna Hostel, Indian Institute of Technology Madras, Chennai - 600 036
Reply below quote. On Mon, Mar 31, 2008 at 08:52:57AM +0530, Kumar Appaiah wrote:
The question is to determine (plot) the step and ramp responses of K / (s^2 + 8s + K) for K = 7, 16 and 80
So, the following code is used, and I'll keep changing a line: <code> from scipy import * from pylab import *
K = 16.0 # should be redone with 7.0, 16.0 and 80.0
r = signal.impulse(([K],[1.0,8.0,K]), T=r_[0:5:0.001]) plot(r[0], r[1], linewidth=2) show() </code>
The above code plots the impulse response of the given function. Now, on to the action:
1. Running the above code with K = 7.0 and K = 16.0 gives expected results. However, running the code with K = 16.0 doesn't work; if K = 16.0, the response should be 16 * t * exp(-4t) * u(t), which it isn't. Changing K to 16 + or - .00000000001 fixes it. There surely is some problem when the value of K is such that a double root is hit.
Well, I've zeroed in on the issue. Running the code of signal.lti.impulse, we get here: s,v = linalg.eig(sys.A) vi = linalg.inv(v) Now, v is a matrix which has the eigen vectors, which are EQUAL in this case: print sys.A [[-0.9701425 -0.9701425 ] [ 0.24253563 0.24253563]] which is singular. However, the code goes on to invert this happily, resulting in a bad matrix and horrendous values. I would be surprised if nobody has encountered this yet (didn't find this on the Tickets). What would be the best way to handle repeated roots in the transfer function's denominator? Thanks. Kumar -- Kumar Appaiah, 458, Jamuna Hostel, Indian Institute of Technology Madras, Chennai - 600 036
I have seen the problem with repeated roots. I thought I reported it. Your second case with a ramp response could be better handled with signal.lsim using your original transfer function ([K],[1.0,8.0,K]) and u = t. On Mon, Mar 31, 2008 at 11:50 AM, Kumar Appaiah <akumar@iitm.ac.in> wrote:
Reply below quote.
On Mon, Mar 31, 2008 at 08:52:57AM +0530, Kumar Appaiah wrote:
The question is to determine (plot) the step and ramp responses of K / (s^2 + 8s + K) for K = 7, 16 and 80
So, the following code is used, and I'll keep changing a line: <code> from scipy import * from pylab import *
K = 16.0 # should be redone with 7.0, 16.0 and 80.0
r = signal.impulse(([K],[1.0,8.0,K]), T=r_[0:5:0.001]) plot(r[0], r[1], linewidth=2) show() </code>
The above code plots the impulse response of the given function. Now, on to the action:
1. Running the above code with K = 7.0 and K = 16.0 gives expected results. However, running the code with K = 16.0 doesn't work; if K = 16.0, the response should be 16 * t * exp(-4t) * u(t), which it isn't. Changing K to 16 + or - .00000000001 fixes it. There surely is some problem when the value of K is such that a double root is hit.
Well, I've zeroed in on the issue. Running the code of signal.lti.impulse, we get here:
s,v = linalg.eig(sys.A) vi = linalg.inv(v)
Now, v is a matrix which has the eigen vectors, which are EQUAL in this case:
print sys.A [[-0.9701425 -0.9701425 ] [ 0.24253563 0.24253563]]
which is singular. However, the code goes on to invert this happily, resulting in a bad matrix and horrendous values. I would be surprised if nobody has encountered this yet (didn't find this on the Tickets).
What would be the best way to handle repeated roots in the transfer function's denominator?
Thanks.
Kumar -- Kumar Appaiah, 458, Jamuna Hostel, Indian Institute of Technology Madras, Chennai - 600 036 _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
On Tue, Apr 01, 2008 at 01:16:28PM -0500, Ryan Krauss wrote:
I have seen the problem with repeated roots. I thought I reported it.
Oh, I'll look again. Even otherwise, I would request someone to please fix this (I am looking for the right method to alter the code to handle this).
Your second case with a ramp response could be better handled with signal.lsim using your original transfer function ([K],[1.0,8.0,K]) and u = t.
Many thanks. I shall try this. Kumar -- Kumar Appaiah, 458, Jamuna Hostel, Indian Institute of Technology Madras, Chennai - 600 036
participants (2)
-
Kumar Appaiah
-
Ryan Krauss