Mercurial > octave
annotate scripts/polynomial/residue.m @ 32632:2e484f9f1f18 stable
maint: update Octave Project Developers copyright for the new year
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Fri, 22 Dec 2023 12:08:17 -0500 |
| parents | 597f3ee61a48 |
| children |
| rev | line source |
|---|---|
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 ######################################################################## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 ## |
|
32632
2e484f9f1f18
maint: update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
31706
diff
changeset
|
3 ## Copyright (C) 1994-2024 The Octave Project Developers |
|
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 ## |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 ## See the file COPYRIGHT.md in the top-level directory of this |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 ## distribution or <https://octave.org/copyright/>. |
| 2313 | 7 ## |
| 8 ## This file is part of Octave. | |
| 9 ## | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23573
diff
changeset
|
10 ## Octave is free software: you can redistribute it and/or modify it |
| 2313 | 11 ## under the terms of the GNU General Public License as published by |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23573
diff
changeset
|
12 ## the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22609
diff
changeset
|
13 ## (at your option) any later version. |
| 2313 | 14 ## |
| 15 ## Octave is distributed in the hope that it will be useful, but | |
| 16 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22609
diff
changeset
|
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22609
diff
changeset
|
18 ## GNU General Public License for more details. |
| 2313 | 19 ## |
| 20 ## You should have received a copy of the GNU General Public License | |
| 7016 | 21 ## along with Octave; see the file COPYING. If not, see |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23573
diff
changeset
|
22 ## <https://www.gnu.org/licenses/>. |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 ## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 ######################################################################## |
| 904 | 25 |
| 3368 | 26 ## -*- texinfo -*- |
|
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
27 ## @deftypefn {} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a}) |
|
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
28 ## @deftypefnx {} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) |
|
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
29 ## @deftypefnx {} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) |
|
13127
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
30 ## The first calling form computes the partial fraction expansion for the |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
31 ## quotient of the polynomials, @var{b} and @var{a}. |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
32 ## |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
33 ## The quotient is defined as |
| 3368 | 34 ## @tex |
| 35 ## $$ | |
| 6978 | 36 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
| 3368 | 37 ## + \sum_{i=1}^N k_i s^{N-i}. |
| 38 ## $$ | |
| 39 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
40 ## @ifnottex |
| 3426 | 41 ## |
| 3368 | 42 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
43 ## @group |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
44 ## B(s) M r(m) N |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
45 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
46 ## A(s) m=1 (s-p(m))^e(m) i=1 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
47 ## @end group |
| 3368 | 48 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
49 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
50 ## @end ifnottex |
| 3368 | 51 ## @noindent |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
52 ## where @math{M} is the number of poles (the length of the @var{r}, @var{p}, |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
53 ## and @var{e}), the @var{k} vector is a polynomial of order @math{N-1} |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
54 ## representing the direct contribution, and the @var{e} vector specifies the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
55 ## multiplicity of the m-th residue's pole. |
| 3426 | 56 ## |
| 6964 | 57 ## For example, |
| 3426 | 58 ## |
| 3368 | 59 ## @example |
| 60 ## @group | |
| 6964 | 61 ## b = [1, 1, 1]; |
| 62 ## a = [1, -5, 8, -4]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
63 ## [r, p, k, e] = residue (b, a) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
64 ## @result{} r = [-2; 7; 3] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
65 ## @result{} p = [2; 2; 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
66 ## @result{} k = [](0x0) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
67 ## @result{} e = [1; 2; 1] |
| 3368 | 68 ## @end group |
| 69 ## @end example | |
| 3426 | 70 ## |
| 3368 | 71 ## @noindent |
| 6978 | 72 ## which represents the following partial fraction expansion |
| 3368 | 73 ## @tex |
| 74 ## $$ | |
| 75 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} | |
| 76 ## $$ | |
| 77 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
78 ## @ifnottex |
| 3426 | 79 ## |
| 3368 | 80 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
81 ## @group |
| 3368 | 82 ## s^2 + s + 1 -2 7 3 |
| 2311 | 83 ## ------------------- = ----- + ------- + ----- |
| 84 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
85 ## @end group |
| 3368 | 86 ## @end example |
| 6964 | 87 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
88 ## @end ifnottex |
| 6978 | 89 ## |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
90 ## The second calling form performs the inverse operation and computes the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
91 ## reconstituted quotient of polynomials, @var{b}(s)/@var{a}(s), from the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
92 ## partial fraction expansion; represented by the residues, poles, and a direct |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
93 ## polynomial specified by @var{r}, @var{p} and @var{k}, and the pole |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
94 ## multiplicity @var{e}. |
| 7011 | 95 ## |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
96 ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
97 ## is determined by the function @code{mpoles}. |
| 6978 | 98 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
99 ## For example: |
| 6964 | 100 ## |
| 101 ## @example | |
| 102 ## @group | |
| 7011 | 103 ## r = [-2; 7; 3]; |
| 104 ## p = [2; 2; 1]; | |
| 105 ## k = [1, 0]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
106 ## [b, a] = residue (r, p, k) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
107 ## @result{} b = [1, -5, 9, -3, 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
108 ## @result{} a = [1, -5, 8, -4] |
| 7011 | 109 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
110 ## where mpoles is used to determine e = [1; 2; 1] |
| 7011 | 111 ## @end group |
| 112 ## @end example | |
| 113 ## | |
| 114 ## Alternatively the multiplicity may be defined explicitly, for example, | |
| 115 ## | |
| 116 ## @example | |
| 117 ## @group | |
| 118 ## r = [7; 3; -2]; | |
| 119 ## p = [2; 1; 2]; | |
| 120 ## k = [1, 0]; | |
| 121 ## e = [2; 1; 1]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
122 ## [b, a] = residue (r, p, k, e) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
123 ## @result{} b = [1, -5, 9, -3, 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
124 ## @result{} a = [1, -5, 8, -4] |
| 6964 | 125 ## @end group |
| 126 ## @end example | |
| 127 ## | |
| 128 ## @noindent | |
| 6978 | 129 ## which represents the following partial fraction expansion |
| 6964 | 130 ## @tex |
| 131 ## $$ | |
| 6978 | 132 ## {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s = {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} |
| 6964 | 133 ## $$ |
| 134 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
135 ## @ifnottex |
| 6964 | 136 ## |
| 137 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
138 ## @group |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
139 ## -2 7 3 s^4 - 5s^3 + 9s^2 - 3s + 1 |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
140 ## ----- + ------- + ----- + s = -------------------------- |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
141 ## (s-2) (s-2)^2 (s-1) s^3 - 5s^2 + 8s - 4 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
142 ## @end group |
| 6964 | 143 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
144 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
145 ## @end ifnottex |
|
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
13963
diff
changeset
|
146 ## @seealso{mpoles, poly, roots, conv, deconv} |
| 3368 | 147 ## @end deftypefn |
| 1025 | 148 |
| 6978 | 149 function [r, p, k, e] = residue (b, a, varargin) |
| 559 | 150 |
| 7011 | 151 if (nargin < 2 || nargin > 4) |
| 6046 | 152 print_usage (); |
| 559 | 153 endif |
| 154 | |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
155 tol = .001; |
| 6964 | 156 |
| 7011 | 157 if (nargin >= 3) |
| 158 if (nargin >= 4) | |
| 159 e = varargin{2}; | |
| 160 else | |
| 161 e = []; | |
| 162 endif | |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
163 ## The inputs are the residue, pole, and direct part. |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
164 ## Solve for the corresponding numerator and denominator polynomials. |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
165 [r, p] = rresidue (b, a, varargin{1}, tol, e); |
|
17312
088d014a7fe2
Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents:
14868
diff
changeset
|
166 return; |
| 7011 | 167 endif |
| 559 | 168 |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
169 ## Make sure both polynomials are in reduced form, and scaled. |
| 1025 | 170 a = polyreduce (a); |
| 171 b = polyreduce (b); | |
| 559 | 172 |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
173 b /= a(1); |
|
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
174 a /= a(1); |
| 559 | 175 |
| 1025 | 176 la = length (a); |
| 177 lb = length (b); | |
| 559 | 178 |
| 2303 | 179 ## Handle special cases here. |
| 1025 | 180 if (la == 0 || lb == 0) |
| 559 | 181 k = r = p = e = []; |
| 182 return; | |
| 183 elseif (la == 1) | |
| 1025 | 184 k = b / a; |
| 559 | 185 r = p = e = []; |
| 186 return; | |
| 187 endif | |
| 188 | |
| 2303 | 189 ## Find the poles. |
| 1025 | 190 p = roots (a); |
| 191 lp = length (p); | |
| 559 | 192 |
| 6964 | 193 ## Sort poles so that multiplicity loop will work. |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
194 [e, idx] = mpoles (p, tol, 1); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
195 p = p(idx); |
| 559 | 196 |
| 7398 | 197 ## For each group of pole multiplicity, set the value of each |
|
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21546
diff
changeset
|
198 ## pole to the average of the group. This reduces the error in |
| 7398 | 199 ## the resulting poles. |
| 200 p_group = cumsum (e == 1); | |
| 201 for ng = 1:p_group(end) | |
| 202 m = find (p_group == ng); | |
| 203 p(m) = mean (p(m)); | |
| 204 endfor | |
| 205 | |
| 2303 | 206 ## Find the direct term if there is one. |
| 1025 | 207 if (lb >= la) |
| 6964 | 208 ## Also return the reduced numerator. |
| 1025 | 209 [k, b] = deconv (b, a); |
| 210 lb = length (b); | |
| 559 | 211 else |
| 212 k = []; | |
| 213 endif | |
| 214 | |
| 7398 | 215 ## Determine if the poles are (effectively) zero. |
| 216 small = max (abs (p)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
217 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
218 small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
219 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
220 small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
221 endif |
| 7398 | 222 p(abs (p) < small) = 0; |
| 223 | |
| 224 ## Determine if the poles are (effectively) real, or imaginary. | |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
225 idx = (abs (imag (p)) < small); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
226 p(idx) = real (p(idx)); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
227 idx = (abs (real (p)) < small); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
228 p(idx) = 1i * imag (p(idx)); |
| 7398 | 229 |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
230 ## The remainder determines the residues. The case of one pole is trivial. |
| 1025 | 231 if (lp == 1) |
| 232 r = polyval (b, p); | |
| 559 | 233 return; |
| 234 endif | |
| 235 | |
| 6964 | 236 ## Determine the order of the denominator and remaining numerator. |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
237 ## With the direct term removed, the potential order of the numerator |
| 6964 | 238 ## is one less than the order of the denominator. |
| 239 aorder = numel (a) - 1; | |
| 240 border = aorder - 1; | |
| 1025 | 241 |
| 6964 | 242 ## Construct a system of equations relating the individual |
| 243 ## contributions from each residue to the complete numerator. | |
| 244 A = zeros (border+1, border+1); | |
| 245 B = prepad (reshape (b, [numel(b), 1]), border+1, 0); | |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14552
diff
changeset
|
246 for ip = 1:numel (p) |
| 6964 | 247 ri = zeros (size (p)); |
| 248 ri(ip) = 1; | |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
249 A(:,ip) = prepad (rresidue (ri, p, [], tol), border+1, 0).'; |
| 6964 | 250 endfor |
| 559 | 251 |
| 2303 | 252 ## Solve for the residues. |
|
25378
8e30ef611002
residue.m: Use a pre-conditioner on final matrix division for numerical stability (bug #53869)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
253 ## FIXME: Use a pre-conditioner d to make A \ B work better (bug #53869). |
|
8e30ef611002
residue.m: Use a pre-conditioner on final matrix division for numerical stability (bug #53869)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
254 ## It would be better to construct A and B so they are not close to |
|
8e30ef611002
residue.m: Use a pre-conditioner on final matrix division for numerical stability (bug #53869)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
255 ## singular in the first place. |
|
8e30ef611002
residue.m: Use a pre-conditioner on final matrix division for numerical stability (bug #53869)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
256 d = max (abs (A), [], 2); |
|
8e30ef611002
residue.m: Use a pre-conditioner on final matrix division for numerical stability (bug #53869)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
257 r = (diag (d) \ A) \ (B ./ d); |
| 559 | 258 |
| 259 endfunction | |
| 6964 | 260 |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
261 ## Reconstitute the numerator and denominator polynomials |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
262 ## from the residues, poles, and direct term. |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
263 function [pnum, pden, e] = rresidue (r, p, k = [], tol = [], e = []) |
| 6964 | 264 |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
265 if (! isempty (e)) |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
266 idx = 1:numel (p); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
267 else |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
268 [e, idx] = mpoles (p, tol, 0); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
269 p = p(idx); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
270 r = r(idx); |
| 7011 | 271 endif |
| 272 | |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
273 idx = 1:numel (p); |
|
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
274 for n = idx |
| 6964 | 275 pn = [1, -p(n)]; |
|
14552
86854d032a37
maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents:
14363
diff
changeset
|
276 if (n == 1) |
| 6964 | 277 pden = pn; |
| 278 else | |
| 279 pden = conv (pden, pn); | |
| 280 endif | |
| 281 endfor | |
| 282 | |
| 283 ## D is the order of the denominator | |
| 284 ## K is the order of the direct polynomial | |
| 285 ## N is the order of the resulting numerator | |
| 286 ## pnum(1:(N+1)) is the numerator's polynomial | |
| 287 ## pden(1:(D+1)) is the denominator's polynomial | |
|
27956
2310164737b3
fix many spelling errors (bug #57613)
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
288 ## pm is the multiple pole for the nth residue |
| 6964 | 289 ## pn is the numerator contribution for the nth residue |
| 290 | |
| 291 D = numel (pden) - 1; | |
| 292 K = numel (k) - 1; | |
| 293 N = K + D; | |
| 294 pnum = zeros (1, N+1); | |
|
25547
2b9a30925a9c
residue.m: Clean up variable names, spacing, etc.
Rik <rik@octave.org>
parents:
25378
diff
changeset
|
295 for n = idx(abs (r) > 0) |
| 6964 | 296 p1 = [1, -p(n)]; |
|
23558
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
297 pn = 1; |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
298 for j = 1:n - 1 |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
299 pn = conv (pn, [1, -p(j)]); |
|
28947
b73a54c31152
maint: Use Octave-specific end keywords rather than just 'end'.
Rik <rik@octave.org>
parents:
27985
diff
changeset
|
300 endfor |
|
23558
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
301 for j = n + 1:numel (p) |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
302 pn = conv (pn, [1, -p(j)]); |
|
28947
b73a54c31152
maint: Use Octave-specific end keywords rather than just 'end'.
Rik <rik@octave.org>
parents:
27985
diff
changeset
|
303 endfor |
|
23558
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
304 for j = 1:e(n) - 1 |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
305 pn = deconv (pn, p1); |
|
28947
b73a54c31152
maint: Use Octave-specific end keywords rather than just 'end'.
Rik <rik@octave.org>
parents:
27985
diff
changeset
|
306 endfor |
| 6964 | 307 pn = r(n) * pn; |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
308 pnum += prepad (pn, N+1, 0, 2); |
| 6964 | 309 endfor |
| 310 | |
| 311 ## Add the direct term. | |
| 312 if (numel (k)) | |
|
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
313 pnum += conv (pden, k); |
| 6964 | 314 endif |
| 315 | |
| 316 pnum = polyreduce (pnum); | |
| 317 pden = polyreduce (pden); | |
| 318 | |
| 319 endfunction | |
| 6968 | 320 |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
321 |
| 6968 | 322 %!test |
| 323 %! b = [1, 1, 1]; | |
| 324 %! a = [1, -5, 8, -4]; | |
| 325 %! [r, p, k, e] = residue (b, a); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
326 %! assert (r, [-2; 7; 3], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
327 %! assert (p, [2; 2; 1], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
328 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
329 %! assert (e, [1; 2; 1]); |
| 6994 | 330 %! k = [1 0]; |
| 7011 | 331 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); |
| 332 %! a = a; | |
| 333 %! [br, ar] = residue (r, p, k); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
334 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
335 %! assert (ar, a, 1e-12); |
| 7011 | 336 %! [br, ar] = residue (r, p, k, e); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
337 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
338 %! assert (ar, a, 1e-12); |
| 6994 | 339 |
| 340 %!test | |
| 341 %! b = [1, 0, 1]; | |
| 342 %! a = [1, 0, 18, 0, 81]; | |
| 7398 | 343 %! [r, p, k, e] = residue (b, a); |
| 6994 | 344 %! r1 = [-5i; 12; +5i; 12]/54; |
| 345 %! p1 = [+3i; +3i; -3i; -3i]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
346 %! assert (r, r1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
347 %! assert (p, p1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
348 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
349 %! assert (e, [1; 2; 1; 2]); |
| 6994 | 350 %! [br, ar] = residue (r, p, k); |
|
25378
8e30ef611002
residue.m: Use a pre-conditioner on final matrix division for numerical stability (bug #53869)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
351 %! assert (br, b, 1e-12); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
352 %! assert (ar, a, 1e-12); |
| 7011 | 353 |
| 354 %!test | |
| 355 %! r = [7; 3; -2]; | |
| 356 %! p = [2; 1; 2]; | |
| 357 %! k = [1 0]; | |
| 358 %! e = [2; 1; 1]; | |
| 359 %! [b, a] = residue (r, p, k, e); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
360 %! assert (b, [1, -5, 9, -3, 1], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
361 %! assert (a, [1, -5, 8, -4], 1e-12); |
| 7011 | 362 %! [rr, pr, kr, er] = residue (b, a); |
|
29533
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
363 %! [~, m] = mpoles (rr); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
364 %! [~, n] = mpoles (r); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
365 %! assert (rr(m), r(n), 1e-12); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
366 %! assert (pr(m), p(n), 1e-12); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
367 %! assert (kr, k, 1e-12); |
|
29533
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
368 %! assert (er(m), e(n), 1e-12); |
| 7011 | 369 |
| 7188 | 370 %!test |
| 371 %! b = [1]; | |
| 372 %! a = [1, 10, 25]; | |
| 7398 | 373 %! [r, p, k, e] = residue (b, a); |
| 7188 | 374 %! r1 = [0; 1]; |
| 375 %! p1 = [-5; -5]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
376 %! assert (r, r1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
377 %! assert (p, p1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
378 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
379 %! assert (e, [1; 2]); |
| 7188 | 380 %! [br, ar] = residue (r, p, k); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
381 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
382 %! assert (ar, a, 1e-12); |
|
13128
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
383 |
|
22489
93ea313301f9
test: Add bug ids (<#####>) to BIST tests.
Rik <rik@octave.org>
parents:
22470
diff
changeset
|
384 ## The following test is due to Bernard Grung |
|
23573
1b4f4ec53b4a
use new script to tag fixed bugs in tests
John W. Eaton <jwe@octave.org>
parents:
23572
diff
changeset
|
385 %!test <*34266> |
|
13128
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
386 %! z1 = 7.0372976777e6; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
387 %! p1 = -3.1415926536e9; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
388 %! p2 = -4.9964813512e8; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
389 %! r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
390 %! r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
391 %! r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
392 %! r4 = z1/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
393 %! r = [r1; r2; r3; r4]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
394 %! p = [p1; p2; 0; 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
395 %! k = []; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
396 %! e = [1; 1; 1; 2]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
397 %! b = [1, z1]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
398 %! a = [1, -(p1 + p2), p1*p2, 0, 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
399 %! [br, ar] = residue (r, p, k, e); |
|
23544
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
400 %! assert (br, [0,0,b], 1e-7); |
|
13128
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
401 %! assert (ar, a, 1e-8); |
|
23544
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
402 |
|
23573
1b4f4ec53b4a
use new script to tag fixed bugs in tests
John W. Eaton <jwe@octave.org>
parents:
23572
diff
changeset
|
403 %!test <*49291> |
|
23544
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
404 %! rf = [1e3, 2e3, 1e3, 2e3]; |
|
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
405 %! cf = [316.2e-9, 50e-9, 31.6e-9, 5e-9]; |
|
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
406 %! [num, den] = residue (1./cf,-1./(rf.*cf),0); |
|
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
407 %! assert (numel (num), 4); |
|
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
408 %! assert (numel (den), 5); |
|
7eba3cd5b6e0
residue.m: Remove code that filters out small return values (bug #34266, bug #49291).
Rik <rik@octave.org>
parents:
23219
diff
changeset
|
409 %! assert (den(1), 1); |
|
23558
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
410 |
|
23573
1b4f4ec53b4a
use new script to tag fixed bugs in tests
John W. Eaton <jwe@octave.org>
parents:
23572
diff
changeset
|
411 %!test <*51148> |
|
23558
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
412 %! r = [1.0000e+18, 3.5714e+12, 2.2222e+11, 2.1739e+10]; |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
413 %! pin = [-1.9231e+15, -1.6234e+09, -4.1152e+07, -1.8116e+06]; |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
414 %! k = 0; |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
415 %! [p, q] = residue (r, pin, k); |
|
5faf1767f9a9
Improve accuracy of residue for inputs with very different magnitudes (bug #51148).
Marco Caliari <marco.caliari@univr.it>
parents:
23544
diff
changeset
|
416 %! assert (p(4), 4.6828e+42, -1e-5); |
|
29533
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
417 |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
418 %!test <*60384> |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
419 %! B = [1315.789473684211]; |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
420 %! A = [1, 1.100000536842105e+04, 1.703789473684211e+03, 0]; |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
421 %! poles1 = roots (A); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
422 %! [r, p, k, e] = residue (B, A); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
423 %! [B1, A1] = residue (r, p, k, e); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
424 %! assert (B, B1); |
|
6dfc06f55cd2
mpoles.m: Fix detection of pole multiplicity (bug #60384).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
425 %! assert (A, A1); |