|
559
|
1 function [r, p, k, e] = residue(b,a,toler) |
|
904
|
2 |
|
|
3 # [r p k e] = residue(b,a) |
|
|
4 # If b and a are vectors of polynomial coefficients, then residue |
|
|
5 # calculates the partial fraction expansion corresponding to the |
|
|
6 # ratio of the two polynomials. The vector r contains the residue |
|
|
7 # terms, p contains the pole values, k contains the coefficients of |
|
|
8 # a direct polynomial term (if it exists) and e is a vector containing |
|
|
9 # the powers of the denominators in the partial fraction terms. |
|
|
10 # Assuming b and a represent polynomials P(s) and Q(s) we have: |
|
559
|
11 # |
|
904
|
12 # P(s) M r(m) N |
|
|
13 # ---- = # ------------- + # k(n)*s^(N-n) |
|
|
14 # Q(s) m=1 (s-p(m))^e(m) n=1 |
|
559
|
15 # |
|
904
|
16 # (# represents summation) where M is the number of poles (the length of |
|
|
17 # the r, p, and e vectors) and N is the length of the k vector. |
|
559
|
18 # |
|
904
|
19 # [r p k e] = residue(b,a,tol) |
|
|
20 # This form of the function call may be used to set a tolerance value. |
|
|
21 # The default value is 0.001. The tolerance value is used to determine |
|
|
22 # whether poles with small imaginary components are declared real. It is |
|
|
23 # also used to determine if two poles are distinct. If the ratio of the |
|
|
24 # imaginary part of a pole to the real part is less than tol, the imaginary |
|
|
25 # part is discarded. If two poles are farther apart than tol they are |
|
|
26 # distinct. |
|
559
|
27 # |
|
904
|
28 # Example: |
|
|
29 # b = [1 1 1]; |
|
|
30 # a = [1 -5 8 -4]; |
|
559
|
31 # |
|
904
|
32 # [r p k e] = residue(b,a) |
|
559
|
33 # |
|
904
|
34 # returns |
|
559
|
35 # |
|
904
|
36 # r = [-2 7 3]; p = [2 2 1]; k = []; e = [1 2 1]; |
|
559
|
37 # |
|
904
|
38 # which implies the following partial fraction expansion |
|
559
|
39 # |
|
904
|
40 # s^2 + s + 1 -2 7 3 |
|
|
41 # ------------------- = ----- + ------- + ----- |
|
|
42 # s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) |
|
559
|
43 # |
|
904
|
44 # SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg |
|
559
|
45 |
|
|
46 # Author: |
|
|
47 # Tony Richardson |
|
|
48 # amr@mpl.ucsd.edu |
|
|
49 # June 1994 |
|
|
50 |
|
|
51 # Here's the method used to find the residues. |
|
|
52 # The partial fraction expansion can be written as: |
|
|
53 # |
|
904
|
54 # |
|
559
|
55 # P(s) D M(k) A(k,m) |
|
|
56 # ---- = # # ------------- |
|
|
57 # Q(s) k=1 m=1 (s - pr(k))^m |
|
|
58 # |
|
|
59 # (# is used to represent a summation) where D is the number of |
|
|
60 # distinct roots, pr(k) is the kth distinct root, M(k) is the |
|
|
61 # multiplicity of the root, and A(k,m) is the residue cooresponding |
|
|
62 # to the kth distinct root with multiplicity m. For example, |
|
|
63 # |
|
|
64 # s^2 A(1,1) A(2,1) A(2,2) |
|
|
65 # ------------------- = ------ + ------ + ------- |
|
|
66 # s^3 + 4s^2 + 5s + 2 (s+2) (s+1) (s+1)^2 |
|
|
67 # |
|
|
68 # In this case there are two distinct roots (D=2 and pr = [-2 -1]), |
|
|
69 # the first root has multiplicity one and the second multiplicity |
|
|
70 # two (M = [1 2]) The residues are actually stored in vector format as |
|
|
71 # r = [ A(1,1) A(2,1) A(2,2) ]. |
|
|
72 # |
|
|
73 # We then multiply both sides by Q(s). Continuing the example: |
|
|
74 # |
|
|
75 # s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2) |
|
|
76 # |
|
|
77 # or |
|
|
78 # |
|
|
79 # s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2) |
|
|
80 # |
|
|
81 # The coefficients of the polynomials on the right are stored |
|
|
82 # in a row vector called rhs, while the coefficients of the |
|
|
83 # polynomial on the left is stored in a row vector called lhs. |
|
|
84 # If the multiplicity of any root is greater than one we'll |
|
|
85 # also need derivatives of this equation of order up to the |
|
|
86 # maximum multiplicity minus one. The derivative coefficients |
|
|
87 # are stored in successive rows of lhs and rhs. |
|
|
88 # |
|
|
89 # For our example lhs and rhs would be: |
|
|
90 # |
|
|
91 # | 1 0 0 | |
|
|
92 # lhs = | | |
|
|
93 # | 0 2 0 | |
|
|
94 # |
|
|
95 # | 1 2 1 1 3 2 0 1 2 | |
|
|
96 # rhs = | | |
|
|
97 # | 0 2 2 0 2 3 0 0 1 | |
|
|
98 # |
|
|
99 # We then form a vector B and a matrix A obtained by evaluating the |
|
|
100 # polynomials in lhs and rhs at the pole values. If a pole has a |
|
|
101 # multiplicity greater than one we also evaluate the derivative |
|
|
102 # polynomials (successive rows) at the pole value. |
|
|
103 # |
|
|
104 # For our example we would have |
|
|
105 # |
|
|
106 # | 4| | 1 0 0 | | r(1) | |
|
|
107 # | 1| = | 0 0 1 | * | r(2) | |
|
|
108 # |-2| | 0 1 1 | | r(3) | |
|
|
109 # |
|
|
110 # We then solve for the residues using matrix division. |
|
|
111 |
|
|
112 if(nargin < 2 || nargin > 3) |
|
904
|
113 usage ("residue(b,a[,toler])"); |
|
559
|
114 endif |
|
|
115 |
|
|
116 if (nargin == 2) |
|
|
117 # Set the default tolerance level |
|
|
118 toler = .001; |
|
|
119 endif |
|
|
120 |
|
|
121 # Make sure both polynomials are in reduced form. |
|
|
122 a = polyreduce(a); |
|
|
123 b = polyreduce(b); |
|
|
124 |
|
|
125 b = b/a(1); |
|
|
126 a = a/a(1); |
|
|
127 |
|
|
128 la = length(a); |
|
|
129 lb = length(b); |
|
|
130 |
|
|
131 # Handle special cases here. |
|
|
132 if(la == 0 || lb == 0) |
|
|
133 k = r = p = e = []; |
|
|
134 return; |
|
|
135 elseif (la == 1) |
|
|
136 k = b/a; |
|
|
137 r = p = e = []; |
|
|
138 return; |
|
|
139 endif |
|
|
140 |
|
|
141 # Find the poles. |
|
|
142 p = roots(a); |
|
|
143 lp = length(p); |
|
|
144 |
|
|
145 # Determine if the poles are (effectively) real. |
|
|
146 index = find(abs(imag(p) ./ real(p)) < toler); |
|
|
147 if (length(index) != 0) |
|
|
148 p(index) = real(p(index)); |
|
|
149 endif |
|
|
150 |
|
|
151 # Find the direct term if there is one. |
|
|
152 if(lb>=la) |
|
|
153 # Also returns the reduced numerator. |
|
|
154 [k, b] = deconv(b,a); |
|
|
155 lb = length(b); |
|
|
156 else |
|
|
157 k = []; |
|
|
158 endif |
|
|
159 |
|
|
160 if(lp == 1) |
|
|
161 r = polyval(b,p); |
|
|
162 e = 1; |
|
|
163 return; |
|
|
164 endif |
|
|
165 |
|
|
166 |
|
|
167 # We need to determine the number and multiplicity of the roots. |
|
|
168 # D is the number of distinct roots. |
|
|
169 # M is a vector of length D containing the multiplicity of each root. |
|
|
170 # pr is a vector of length D containing only the distinct roots. |
|
|
171 # e is a vector of length lp which indicates the power in the partial |
|
|
172 # fraction expansion of each term in p. |
|
|
173 |
|
|
174 # Set initial values. We'll remove elements from pr as we find |
|
|
175 # multiplicities. We'll shorten M afterwards. |
|
|
176 e = ones(lp,1); |
|
|
177 M = zeros(lp,1); |
|
|
178 pr = p; |
|
|
179 D = 1; M(1) = 1; |
|
|
180 |
|
|
181 old_p_index = 1; new_p_index = 2; |
|
|
182 M_index = 1; pr_index = 2; |
|
|
183 while(new_p_index<=lp) |
|
|
184 if(abs(p(new_p_index) - p(old_p_index)) < toler) |
|
|
185 # We've found a multiple pole. |
|
|
186 M(M_index) = M(M_index) + 1; |
|
|
187 e(new_p_index) = e(new_p_index-1) + 1; |
|
|
188 # Remove the pole from pr. |
|
|
189 pr(pr_index) = []; |
|
|
190 else |
|
|
191 # It's a different pole. |
|
|
192 D++; M_index++; M(M_index) = 1; |
|
|
193 old_p_index = new_p_index; pr_index++; |
|
|
194 endif |
|
|
195 new_p_index++; |
|
|
196 endwhile |
|
|
197 |
|
|
198 # Shorten M to it's proper length |
|
|
199 M = M(1:D); |
|
|
200 |
|
|
201 # Now set up the polynomial matrices. |
|
|
202 |
|
|
203 MM = max(M); |
|
|
204 # Left hand side polynomial |
|
|
205 lhs = zeros(MM,lb); |
|
|
206 rhs = zeros(MM,lp*lp); |
|
|
207 lhs(1,:) = b; |
|
|
208 rhi = 1; |
|
|
209 dpi = 1; |
|
|
210 mpi = 1; |
|
|
211 while(dpi<=D) |
|
|
212 for ind = 1:M(dpi) |
|
|
213 if(mpi>1 && (mpi+ind)<=lp) |
|
|
214 cp = [p(1:mpi-1); p(mpi+ind:lp)]; |
|
|
215 elseif (mpi==1) |
|
|
216 cp = p(mpi+ind:lp); |
|
|
217 else |
|
|
218 cp = p(1:mpi-1); |
|
|
219 endif |
|
|
220 rhs(1,rhi:rhi+lp-1) = prepad(poly(cp),lp); |
|
|
221 rhi = rhi + lp; |
|
|
222 endfor |
|
|
223 mpi = mpi + M(dpi); |
|
|
224 dpi++; |
|
|
225 endwhile |
|
|
226 if(MM > 1) |
|
|
227 for index = 2:MM |
|
|
228 lhs(index,:) = prepad(polyderiv(lhs(index-1,:)),lb); |
|
|
229 ind = 1; |
|
|
230 for rhi = 1:lp |
|
|
231 cp = rhs(index-1,ind:ind+lp-1); |
|
|
232 rhs(index,ind:ind+lp-1) = prepad(polyderiv(cp),lp); |
|
|
233 ind = ind + lp; |
|
|
234 endfor |
|
|
235 endfor |
|
|
236 endif |
|
|
237 |
|
|
238 # Now lhs contains the numerator polynomial and as many derivatives as are |
|
|
239 # required. rhs is a matrix of polynomials, the first row contains the |
|
|
240 # corresponding polynomial for each residue and successive rows are |
|
|
241 # derivatives. |
|
|
242 |
|
|
243 # Now we need to evaluate the first row of lhs and rhs at each distinct |
|
|
244 # pole value. If there are multiple poles we will also need to evaluate |
|
|
245 # the derivatives at the pole value also. |
|
|
246 |
|
|
247 B = zeros(lp,1); |
|
|
248 A = zeros(lp,lp); |
|
|
249 |
|
|
250 dpi = 1; |
|
|
251 row = 1; |
|
|
252 while(dpi<=D) |
|
|
253 for mi = 1:M(dpi) |
|
|
254 B(row) = polyval(lhs(mi,:),pr(dpi)); |
|
|
255 ci = 1; |
|
|
256 for col = 1:lp |
|
|
257 cp = rhs(mi,ci:ci+lp-1); |
|
|
258 A(row,col) = polyval(cp,pr(dpi)); |
|
|
259 ci = ci + lp; |
|
|
260 endfor |
|
|
261 row++; |
|
|
262 endfor |
|
|
263 dpi++; |
|
|
264 endwhile |
|
|
265 |
|
|
266 # Solve for the residues. |
|
|
267 r = A\B; |
|
|
268 |
|
|
269 endfunction |
