Mercurial > hg > octave-lyh
comparison scripts/polynomial/roots.m @ 787:c5d35bb139b6
[project @ 1994-10-11 00:34:13 by jwe]
Initial revision
| author | jwe |
|---|---|
| date | Tue, 11 Oct 1994 00:34:39 +0000 |
| parents | |
| children | 3470f1e25a79 |
comparison
equal
deleted
inserted
replaced
| 786:4fcd2e68dd3b | 787:c5d35bb139b6 |
|---|---|
| 1 function r = roots (v) | |
| 2 # | |
| 3 # For a vector v with n components, return the roots of the | |
| 4 # polynomial v(1) * z^(n-1) + ... + v(n-1) * z + v(n). | |
| 5 | |
| 6 # Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993 | |
| 7 # Copyright Dept of Probability Theory and Statistics TU Wien | |
| 8 | |
| 9 [nr, nc] = size(v); | |
| 10 if !((nr == 1 && nc > 1) || (nc == 1 && nr > 1)) | |
| 11 error ("usage: roots (v), where v is a nonzero vector"); | |
| 12 endif | |
| 13 | |
| 14 n = nr + nc - 1; | |
| 15 v = reshape (v, 1, n); | |
| 16 | |
| 17 # If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the | |
| 18 # leading k zeros and n - k - l roots of the polynomial are zero. | |
| 19 f = find (v); | |
| 20 m = max (size (f)); | |
| 21 if (m > 0) | |
| 22 v = v(f(1):f(m)); | |
| 23 l = max (size (v)); | |
| 24 if (l > 1) | |
| 25 A = diag (ones (1, l-2), -1); | |
| 26 A(1,:) = -v(2:l) ./ v(1); | |
| 27 r = eig (A); | |
| 28 if (f(m) < n) | |
| 29 r = [r; zeros (n - f(m), 1)]; | |
| 30 endif | |
| 31 else | |
| 32 r = zeros (n - f(m), 1); | |
| 33 endif | |
| 34 else | |
| 35 error ("usage: roots(v), where v is a nonzero vector"); | |
| 36 endif | |
| 37 | |
| 38 endfunction |