octave_packages/m/polynomial/roots.m

   1 ## Copyright (C) 1994-2012 John W. Eaton
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## 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 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} roots (@var{v})
  21 ##
  22 ## For a vector @var{v} with @math{N} components, return
  23 ## the roots of the polynomial
  24 ## @tex
  25 ## $$
  26 ## v_1 z^{N-1} + \cdots + v_{N-1} z + v_N.
  27 ## $$
  28 ## @end tex
  29 ## @ifnottex
  30 ##
  31 ## @example
  32 ## v(1) * z^(N-1) + @dots{} + v(N-1) * z + v(N)
  33 ## @end example
  34 ##
  35 ## @end ifnottex
  36 ##
  37 ## As an example, the following code finds the roots of the quadratic
  38 ## polynomial
  39 ## @tex
  40 ## $$ p(x) = x^2 - 5. $$
  41 ## @end tex
  42 ## @ifnottex
  43 ##
  44 ## @example
  45 ## p(x) = x^2 - 5.
  46 ## @end example
  47 ##
  48 ## @end ifnottex
  49 ##
  50 ## @example
  51 ## @group
  52 ## c = [1, 0, -5];
  53 ## roots (c)
  54 ## @result{}  2.2361
  55 ## @result{} -2.2361
  56 ## @end group
  57 ## @end example
  58 ##
  59 ## Note that the true result is
  60 ## @tex
  61 ## $\pm \sqrt{5}$
  62 ## @end tex
  63 ## @ifnottex
  64 ## @math{+/- sqrt(5)}
  65 ## @end ifnottex
  66 ## which is roughly
  67 ## @tex
  68 ## $\pm 2.2361$.
  69 ## @end tex
  70 ## @ifnottex
  71 ## @math{+/- 2.2361}.
  72 ## @end ifnottex
  73 ## @seealso{poly, compan, fzero}
  74 ## @end deftypefn
  75
  76 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
  77 ## Created: 24 December 1993
  78 ## Adapted-By: jwe
  79
  80 function r = roots (v)
  81
  82   if (nargin != 1 || min (size (v)) > 1)
  83     print_usage ();
  84   elseif (any (isnan(v) | isinf(v)))
  85     error ("roots: inputs must not contain Inf or NaN");
  86   endif
  87
  88   n = numel (v);
  89   v = v(:);
  90
  91   ## If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the
  92   ## leading k zeros and n - k - l roots of the polynomial are zero.
  93
  94   if (isempty (v))
  95     f = v;
  96   else
  97     f = find (v ./ max (abs (v)));
  98   endif
  99   m = numel (f);
 100
 101   if (m > 0 && n > 1)
 102     v = v(f(1):f(m));
 103     l = max (size (v));
 104     if (l > 1)
 105       A = diag (ones (1, l-2), -1);
 106       A(1,:) = -v(2:l) ./ v(1);
 107       r = eig (A);
 108       if (f(m) < n)
 109         tmp = zeros (n - f(m), 1);
 110         r = [r; tmp];
 111       endif
 112     else
 113       r = zeros (n - f(m), 1);
 114     endif
 115   else
 116     r = [];
 117   endif
 118
 119 endfunction
 120
 121 %!test
 122 %! p = [poly([3 3 3 3]), 0 0 0 0];
 123 %! r = sort (roots (p));
 124 %! assert (r, [0; 0; 0; 0; 3; 3; 3; 3], 0.001)
 125
 126 %!assert(all (all (abs (roots ([1, -6, 11, -6]) - [3; 2; 1]) < sqrt (eps))));
 127
 128 %!assert(isempty (roots ([])));
 129
 130 %!error roots ([1, 2; 3, 4]);
 131
 132 %!assert(isempty (roots (1)));
 133
 134 %!error roots ([1, 2; 3, 4]);
 135
 136 %!error roots ([1 Inf 1]);
 137
 138 %!error roots ([1 NaN 1]);
 139
 140 %!assert(roots ([1e-200, -1e200, 1]), 1e-200)
 141 %!assert(roots ([1e-200, -1e200 * 1i, 1]), -1e-200 * 1i)