home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / scripts / polynomial / roots.m < prev   
Text File  |  1996-09-28  |  2KB  |  62 lines
  1. # Copyright (C) 1995 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 the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # 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, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function r = roots (v)
  20.  
  21. # usage: roots (v)
  22. #
  23. # For a vector v with n components, return the roots of the
  24. # polynomial v(1) * z^(n-1) + ... + v(n-1) * z + v(n).
  25.   
  26. # Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993.
  27.  
  28.   [nr, nc] = size (v);
  29.   if (nr <= 1 && nc <= 1)
  30.     r = [];
  31.     return;
  32.   elseif (! ((nr == 1 && nc > 1) || (nc == 1 && nr > 1)))
  33.     usage ("roots (v), where v is a nonzero vector");
  34.   endif
  35.  
  36.   n = nr + nc - 1;
  37.   v = reshape (v, 1, n);
  38.  
  39. # If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the
  40. # leading k zeros and n - k - l roots of the polynomial are zero.  
  41.  
  42.   f = find (v);
  43.   m = max (size (f));
  44.   if (m > 0)
  45.     v = v (f(1):f(m));
  46.     l = max (size (v));
  47.     if (l > 1)
  48.       A = diag (ones (1, l-2), -1);
  49.       A (1, :) = -v (2:l) ./ v (1);
  50.       r = eig (A);    
  51.       if (f(m) < n)
  52.     r = [r; zeros (n - f (m), 1)];
  53.       endif
  54.     else
  55.       r = zeros (n - f(m), 1);
  56.     endif
  57.   else
  58.     usage ("roots (v), where v is a nonzero vector");
  59.   endif
  60.   
  61. endfunction
  62.