Mercurial > hg > octave-lyh

comparison scripts/polynomial/residue.m @ 6964:33f20a41aeea

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-06 04:31:18 by jwe]
author jwe
date Sat, 06 Oct 2007 04:31:18 +0000
parents 34f96dd5441b
children c8fc3487ed2c
comparison
equal deleted inserted replaced
6963:642f481d2d50 6964:33f20a41aeea
1 ## Copyright (C) 1996, 1997 John W. Eaton 1 ## Copyright (C) 1996, 1997 John W. Eaton
2 ## Copyright (C) 2007 Ben Abbott
2 ## 3 ##
3 ## This file is part of Octave. 4 ## This file is part of Octave.
4 ## 5 ##
5 ## Octave is free software; you can redistribute it and/or modify it 6 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by 7 ## under the terms of the GNU General Public License as published by
16 ## along with Octave; see the file COPYING. If not, write to the Free 17 ## along with Octave; see the file COPYING. If not, write to the Free
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 18 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 ## 02110-1301, USA. 19 ## 02110-1301, USA.
19 20
20 ## -*- texinfo -*- 21 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {} residue (@var{b}, @var{a}, @var{tol}) 22 ## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}] =} residue (@var{b}, @var{a})
22 ## If @var{b} and @var{a} are vectors of polynomial coefficients, then 23 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k})
23 ## residue calculates the partial fraction expansion corresponding to the 24 ## In the instance of two inputs, they are assumed to correspond to vectors
24 ## ratio of the two polynomials. 25 ## of polynomial coefficients, @var{b} and @var{a}. From these polynomial
26 ## coefficients, the function residue calculates the partial fraction
27 ## expansion corresponding to the ratio of the two polynomials.
25 ## @cindex partial fraction expansion 28 ## @cindex partial fraction expansion
26 ## 29 ##
27 ## The function @code{residue} returns @var{r}, @var{p}, @var{k}, and 30 ## In this instance, the function @code{residue} returns @var{r}, @var{p},
28 ## @var{e}, where the vector @var{r} contains the residue terms, @var{p} 31 ## and @var{k}, where the vector @var{r} contains the residue terms,
29 ## contains the pole values, @var{k} contains the coefficients of a direct 32 ## @var{p} contains the pole values, @var{k} contains the coefficients of a
30 ## polynomial term (if it exists) and @var{e} is a vector containing the 33 ## direct polynomial term (if it exists).
31 ## powers of the denominators in the partial fraction terms. 34 ##
35 ## In the instance of three inputs, the function 'residue' performs the
36 ## reciprocal task. Meaning the partial fraction expansion is reconstituted
37 ## into the corresponding pair of polynomial coefficients.
32 ## 38 ##
33 ## Assuming @var{b} and @var{a} represent polynomials 39 ## Assuming @var{b} and @var{a} represent polynomials
34 ## @iftex 40 ## @iftex
35 ## @tex 41 ## @tex
36 ## $P(s)$ and $Q(s)$ 42 ## $P(s)$ and $Q(s)$
58 ## @end ifinfo 64 ## @end ifinfo
59 ## 65 ##
60 ## @noindent 66 ## @noindent
61 ## where @math{M} is the number of poles (the length of the @var{r}, 67 ## where @math{M} is the number of poles (the length of the @var{r},
62 ## @var{p}, and @var{e} vectors) and @math{N} is the length of the 68 ## @var{p}, and @var{e} vectors) and @math{N} is the length of the
63 ## @var{k} vector. 69 ## @var{k} vector. The @var{e} vector specifies the multiplicity of the
64 ## 70 ## mth residue's pole.
65 ## The argument @var{tol} is optional, and if not specified, a default 71 ##
66 ## value of 0.001 is assumed. The tolerance value is used to determine 72 ## For example,
67 ## whether poles with small imaginary components are declared real. It is
68 ## also used to determine if two poles are distinct. If the ratio of the
69 ## imaginary part of a pole to the real part is less than @var{tol}, the
70 ## imaginary part is discarded. If two poles are farther apart than
71 ## @var{tol} they are distinct. For example,
72 ## 73 ##
73 ## @example 74 ## @example
74 ## @group 75 ## @group
75 ## b = [1, 1, 1]; 76 ## b = [1, 1, 1];
76 ## a = [1, -5, 8, -4]; 77 ## a = [1, -5, 8, -4];
77 ## [r, p, k, e] = residue (b, a); 78 ## [r, p, k] = residue (b, a);
78 ## @result{} r = [-2, 7, 3] 79 ## @result{} r = [-2, 7, 3]
79 ## @result{} p = [2, 2, 1] 80 ## @result{} p = [2, 2, 1]
80 ## @result{} k = [](0x0) 81 ## @result{} k = [](0x0)
81 ## @result{} e = [1, 2, 1]
82 ## @end group 82 ## @end group
83 ## @end example 83 ## @end example
84 ## 84 ##
85 ## @noindent 85 ## @noindent
86 ## which implies the following partial fraction expansion 86 ## which implies the following partial fraction expansion
96 ## @example 96 ## @example
97 ## s^2 + s + 1 -2 7 3 97 ## s^2 + s + 1 -2 7 3
98 ## ------------------- = ----- + ------- + ----- 98 ## ------------------- = ----- + ------- + -----
99 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) 99 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1)
100 ## @end example 100 ## @end example
101 ##
101 ## @end ifinfo 102 ## @end ifinfo
102 ## @seealso{poly, roots, conv, deconv, polyval, polyderiv, and polyinteg} 103 ## A similar, but reciprocal example, where the fraction's polynomials are
104 ## reconstituted from the residues, poles, and direct term is
105 ##
106 ## @example
107 ## @group
108 ## r = [-2, 7, 3];
109 ## p = [2, 2, 1];
110 ## k = [1 0];
111 ## [b, a] = residue (r, p, k);
112 ## @result{} b = [1, -5, 9, -3, 1]
113 ## @result{} a = [1, -5, 8, 4]
114 ## @end group
115 ## @end example
116 ##
117 ## @noindent
118 ## which implies the following partial fraction expansion
119 ## @iftex
120 ## @tex
121 ## $$
122 ## {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s
123 ## $$
124 ## @end tex
125 ## @end iftex
126 ## @ifinfo
127 ##
128 ## @example
129 ## s^4 - 5s^3 + 9s^2 - 3s + 1 -2 7 3
130 ## -------------------------- = ----- + ------- + ----- + s
131 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1)
132 ## @end example
133 ## @end ifinfo
134 ## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg}
103 ## @end deftypefn 135 ## @end deftypefn
104 136
105 ## Author: Tony Richardson <arichard@stark.cc.oh.us> 137 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
138 ## Author: Ben Abbott <bpabbott@mac.com>
106 ## Created: June 1994 139 ## Created: June 1994
107 ## Adapted-By: jwe 140 ## Adapted-By: jwe
108 141
109 function [r, p, k, e] = residue (b, a, toler) 142 function [r, p, k] = residue (b, a, varargin)
110
111 ## Here's the method used to find the residues.
112 ## The partial fraction expansion can be written as:
113 ##
114 ##
115 ## P(s) D M(k) A(k,m)
116 ## ---- = # # -------------
117 ## Q(s) k=1 m=1 (s - pr(k))^m
118 ##
119 ## (# is used to represent a summation) where D is the number of
120 ## distinct roots, pr(k) is the kth distinct root, M(k) is the
121 ## multiplicity of the root, and A(k,m) is the residue cooresponding
122 ## to the kth distinct root with multiplicity m. For example,
123 ##
124 ## s^2 A(1,1) A(2,1) A(2,2)
125 ## ------------------- = ------ + ------ + -------
126 ## s^3 + 4s^2 + 5s + 2 (s+2) (s+1) (s+1)^2
127 ##
128 ## In this case there are two distinct roots (D=2 and pr = [-2 -1]),
129 ## the first root has multiplicity one and the second multiplicity
130 ## two (M = [1 2]) The residues are actually stored in vector format as
131 ## r = [ A(1,1) A(2,1) A(2,2) ].
132 ##
133 ## We then multiply both sides by Q(s). Continuing the example:
134 ##
135 ## s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2)
136 ##
137 ## or
138 ##
139 ## s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2)
140 ##
141 ## The coefficients of the polynomials on the right are stored in a row
142 ## vector called rhs, while the coefficients of the polynomial on the
143 ## left is stored in a row vector called lhs. If the multiplicity of
144 ## any root is greater than one we'll also need derivatives of this
145 ## equation of order up to the maximum multiplicity minus one. The
146 ## derivative coefficients are stored in successive rows of lhs and
147 ## rhs.
148 ##
149 ## For our example lhs and rhs would be:
150 ##
151 ## | 1 0 0 |
152 ## lhs = | |
153 ## | 0 2 0 |
154 ##
155 ## | 1 2 1 1 3 2 0 1 2 |
156 ## rhs = | |
157 ## | 0 2 2 0 2 3 0 0 1 |
158 ##
159 ## We then form a vector B and a matrix A obtained by evaluating the
160 ## polynomials in lhs and rhs at the pole values. If a pole has a
161 ## multiplicity greater than one we also evaluate the derivative
162 ## polynomials (successive rows) at the pole value.
163 ##
164 ## For our example we would have
165 ##
166 ## | 4| | 1 0 0 | | r(1) |
167 ## | 1| = | 0 0 1 | * | r(2) |
168 ## |-2| | 0 1 1 | | r(3) |
169 ##
170 ## We then solve for the residues using matrix division.
171 143
172 if (nargin < 2 || nargin > 3) 144 if (nargin < 2 || nargin > 3)
173 print_usage (); 145 print_usage ();
174 endif 146 endif
175 147
176 if (nargin == 2) 148 toler = .001;
177 toler = .001; 149
178 endif 150 if (nargin == 3)
151 ## The inputs are the residue, pole, and direct part. Solve for the
152 ## corresponding numerator and denominator polynomials
153 [r, p] = rresidue (b, a, varargin{1}, toler);
154 return
155 end
179 156
180 ## Make sure both polynomials are in reduced form. 157 ## Make sure both polynomials are in reduced form.
181 158
182 a = polyreduce (a); 159 a = polyreduce (a);
183 b = polyreduce (b); 160 b = polyreduce (b);
204 p = roots (a); 181 p = roots (a);
205 lp = length (p); 182 lp = length (p);
206 183
207 ## Determine if the poles are (effectively) zero. 184 ## Determine if the poles are (effectively) zero.
208 185
209 p(abs (p) < toler) = 0; 186 small = max (abs (p));
210 187 small = max ([small, 1] ) * 1e-8 * (1 + numel (p))^2;
211 ## Determine if the poles are (effectively) real. 188 p(abs (p) < small) = 0;
212 189
213 index = (abs (p) >= toler & (abs (imag (p)) ./ abs (p) < toler)); 190 ## Determine if the poles are (effectively) real, or imaginary.
191
192 index = (abs (imag (p)) < small);
214 p(index) = real (p(index)); 193 p(index) = real (p(index));
194 index = (abs (real (p)) < small);
195 p(index) = 1i * imag (p(index));
196
197 ## Sort poles so that multiplicity loop will work.
198
199 [e, indx] = mpoles (p, toler, 1);
200 p = p (indx);
215 201
216 ## Find the direct term if there is one. 202 ## Find the direct term if there is one.
217 203
218 if (lb >= la) 204 if (lb >= la)
219 ## Also returns the reduced numerator. 205 ## Also return the reduced numerator.
220 [k, b] = deconv (b, a); 206 [k, b] = deconv (b, a);
221 lb = length (b); 207 lb = length (b);
222 else 208 else
223 k = []; 209 k = [];
224 endif 210 endif
225 211
226 if (lp == 1) 212 if (lp == 1)
227 r = polyval (b, p); 213 r = polyval (b, p);
228 e = 1;
229 return; 214 return;
230 endif 215 endif
231 216
232 ## We need to determine the number and multiplicity of the roots. 217 ## Determine the order of the denominator and remaining numerator.
233 ## 218 ## With the direct term removed the potential order of the numerator
234 ## D is the number of distinct roots. 219 ## is one less than the order of the denominator.
235 ## M is a vector of length D containing the multiplicity of each root. 220
236 ## pr is a vector of length D containing only the distinct roots. 221 aorder = numel (a) - 1;
237 ## e is a vector of length lp which indicates the power in the partial 222 border = aorder - 1;
238 ## fraction expansion of each term in p. 223
239 224 ## Construct a system of equations relating the individual
240 ## Set initial values. We'll remove elements from pr as we find 225 ## contributions from each residue to the complete numerator.
241 ## multiplicities. We'll shorten M afterwards. 226
242 227 A = zeros (border+1, border+1);
243 e = ones (lp, 1); 228 B = prepad (reshape (b, [numel(b), 1]), border+1, 0);
244 M = zeros (lp, 1); 229 for ip = 1:numel(p)
245 pr = p; 230 ri = zeros (size (p));
246 D = 1; 231 ri(ip) = 1;
247 M(1) = 1; 232 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).';
248 233 endfor
249 old_p_index = 1; 234
250 new_p_index = 2; 235 ## Solve for the residues.
251 M_index = 1; 236
252 pr_index = 2; 237 r = A \ B;
253 238
254 while (new_p_index <= lp) 239 endfunction
255 if (abs (p (new_p_index) - p (old_p_index)) < toler) 240
256 ## We've found a multiple pole. 241 function [pnum, pden] = rresidue (r, p, k, toler)
257 M (M_index) = M (M_index) + 1; 242
258 e (new_p_index) = e (new_p_index-1) + 1; 243 ## Reconstitute the numerator and denominator polynomials from the
259 ## Remove the pole from pr. 244 ## residues, poles, and direct term.
260 pr (pr_index) = []; 245
246 if (nargin < 2 || nargin > 4)
247 print_usage ();
248 endif
249
250 if (nargin < 4)
251 toler = [];
252 endif
253
254 if (nargin < 3)
255 k = [];
256 endif
257
258 [multp, indx] = mpoles (p, toler, 0);
259
260 p = p (indx);
261 r = r (indx);
262
263 indx = 1:numel(p);
264
265 for n = indx
266 pn = [1, -p(n)];
267 if n == 1
268 pden = pn;
261 else 269 else
262 ## It's a different pole. 270 pden = conv (pden, pn);
263 D++;
264 M_index++;
265 M (M_index) = 1;
266 old_p_index = new_p_index;
267 pr_index++;
268 endif 271 endif
269 new_p_index++; 272 endfor
270 endwhile 273
271 274 ## D is the order of the denominator
272 ## Shorten M to it's proper length 275 ## K is the order of the direct polynomial
273 276 ## N is the order of the resulting numerator
274 M = M (1:D); 277 ## pnum(1:(N+1)) is the numerator's polynomial
275 278 ## pden(1:(D+1)) is the denominator's polynomial
276 ## Now set up the polynomial matrices. 279 ## pm is the multible pole for the nth residue
277 280 ## pn is the numerator contribution for the nth residue
278 MM = max(M); 281
279 282 D = numel (pden) - 1;
280 ## Left hand side polynomial 283 K = numel (k) - 1;
281 284 N = K + D;
282 lhs = zeros (MM, lb); 285 pnum = zeros (1, N+1);
283 rhs = zeros (MM, lp*lp); 286 for n = indx(abs(r)>0)
284 lhs (1, :) = b; 287 p1 = [1, -p(n)];
285 rhi = 1; 288 for m = 1:multp(n)
286 dpi = 1; 289 if m == 1
287 mpi = 1; 290 pm = p1;
288 while (dpi <= D)
289 for ind = 1:M(dpi)
290 if (mpi > 1 && (mpi+ind) <= lp)
291 cp = [p(1:mpi-1); p(mpi+ind:lp)];
292 elseif (mpi == 1)
293 cp = p (mpi+ind:lp);
294 else 291 else
295 cp = p (1:mpi-1); 292 pm = conv (pm, p1);
296 endif 293 endif
297 rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp, 0, 2);
298 rhi = rhi + lp;
299 endfor 294 endfor
300 mpi = mpi + M (dpi); 295 pn = deconv (pden, pm);
301 dpi++; 296 pn = r(n) * pn;
302 endwhile 297 pnum = pnum + prepad ( pn, N+1, 0);
303 if (MM > 1) 298 endfor
304 for index = 2:MM 299
305 lhs (index, :) = prepad (polyderiv (lhs (index-1, :)), lb, 0, 2); 300 ## Add the direct term.
306 ind = 1; 301
307 for rhi = 1:lp 302 if (numel (k))
308 cp = rhs (index-1, ind:ind+lp-1); 303 pnum = pnum + conv (pden, k);
309 rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp, 0, 2); 304 endif
310 ind = ind + lp; 305
311 endfor 306 ## Check for leading zeros and trim the polynomial coefficients.
312 endfor 307
313 endif 308 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
314 309
315 ## Now lhs contains the numerator polynomial and as many derivatives as 310 pnum (abs (pnum) < small) = 0;
316 ## are required. rhs is a matrix of polynomials, the first row 311 pden (abs (pden) < small) = 0;
317 ## contains the corresponding polynomial for each residue and 312
318 ## successive rows are derivatives. 313 pnum = polyreduce (pnum);
319 314 pden = polyreduce (pden);
320 ## Now we need to evaluate the first row of lhs and rhs at each
321 ## distinct pole value. If there are multiple poles we will also need
322 ## to evaluate the derivatives at the pole value also.
323
324 B = zeros (lp, 1);
325 A = zeros (lp, lp);
326
327 dpi = 1;
328 row = 1;
329 while (dpi <= D)
330 for mi = 1:M(dpi)
331 B (row) = polyval (lhs (mi, :), pr (dpi));
332 ci = 1;
333 for col = 1:lp
334 cp = rhs (mi, ci:ci+lp-1);
335 A (row, col) = polyval (cp, pr(dpi));
336 ci = ci + lp;
337 endfor
338 row++;
339 endfor
340 dpi++;
341 endwhile
342
343 ## Solve for the residues.
344
345 r = A \ B;
346 315
347 endfunction 316 endfunction