Mercurial > hg > octave-lyh

comparison scripts/polynomial/roots.m @ 2325:b5568c31ee2c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-07-15 22:20:21 by jwe]
author jwe
date Mon, 15 Jul 1996 22:20:21 +0000
parents 5ca126254d15
children 8b262e771614
comparison
equal deleted inserted replaced
2324:fdc6e2f81333 2325:b5568c31ee2c
29 function r = roots (v) 29 function r = roots (v)
30 30
31 if (min (size (v)) > 1 || nargin != 1) 31 if (min (size (v)) > 1 || nargin != 1)
32 usage ("roots (v), where v is a vector"); 32 usage ("roots (v), where v is a vector");
33 endif 33 endif
34 34
35 n = length (v); 35 n = length (v);
36 v = reshape (v, 1, n); 36 v = reshape (v, 1, n);
37 37
38 ## If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the 38 ## If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the
39 ## leading k zeros and n - k - l roots of the polynomial are zero. 39 ## leading k zeros and n - k - l roots of the polynomial are zero.
40 40
41 f = find (v); 41 f = find (v);
42 m = max (size (f)); 42 m = max (size (f));
43 43
44 if (m > 0 && n > 1) 44 if (m > 0 && n > 1)
45 v = v(f(1):f(m)); 45 v = v(f(1):f(m));
46 l = max (size (v)); 46 l = max (size (v));
47 if (l > 1) 47 if (l > 1)
48 A = diag (ones (1, l-2), -1); 48 A = diag (ones (1, l-2), -1);
49 A(1,:) = -v(2:l) ./ v(1); 49 A(1,:) = -v(2:l) ./ v(1);
50 r = eig (A); 50 r = eig (A);
51 if (f(m) < n) 51 if (f(m) < n)
52 tmp = zeros (n - f(m), 1); 52 tmp = zeros (n - f(m), 1);
53 r = [r; tmp]; 53 r = [r; tmp];
54 endif 54 endif
55 else 55 else
56 r = zeros (n - f(m), 1); 56 r = zeros (n - f(m), 1);
57 endif 57 endif
58 else 58 else
59 r = []; 59 r = [];
60 endif 60 endif
61 61
62 endfunction 62 endfunction