home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / roots < prev    next >
Encoding:
Text File  |  1995-11-15  |  1.8 KB  |  71 lines
  1. //-------------------------------------------------------------------------
  2. //
  3. // roots
  4. //
  5. // Syntax: r=roots(a)
  6. //
  7. // This routine finds the roots of the polynomial given by the
  8. // coefficients in the input vector a. If a has N components,
  9. // the polynomial is given by,
  10. //
  11. //  polynomial = a[1]*X^(N-1) + a[2]*X^(N-2) + ... + a[N-1]*X + a[N]
  12. //
  13. // Ref: Golub, G., and Van Loan, C., Matrix Computations, Second Edition,
  14. //      The John Hopkins University Press, Baltimore MA, 1989, p. 369.
  15. //      (Note: The method in Golub and Van Loan was modified to be
  16. //             compatible with the Transfer function notation - which
  17. //             is also compatible with the MATLAB version).
  18. //
  19. // Copyright (C), by Jeffrey B. Layton, 1994-95
  20. // Version JBL 950105
  21. // Bill Lash -- Bug Fix for initializing n
  22. //-------------------------------------------------------------------------
  23.  
  24. rfile isvec
  25.  
  26. roots = function(a)
  27. {
  28.    local(n,inz,nnz,A,c,r)
  29.  
  30. // Check if input is a vector
  31.  
  32.    if (!isvec(a)) {
  33.        error("ROOTS: Input must be a Vector.");
  34.    }
  35.  
  36. // Make sure input is a row vector (if it is a row vector it is a
  37. // little easier to create companion matrix).
  38.  
  39.    if (a.nr > 1) {
  40.        A=a.';    // Make sure it's a row vector
  41.    else
  42.        A=a;
  43.    }
  44.  
  45. // Strip leading zeros.  Strip trailing zeros, but keep them as roots
  46. // at zero (store in roots vector r).
  47.  
  48.    n=max(size(A));
  49.  
  50.    inz=find(abs(A));
  51.    nnz=max(size(inz));
  52.    if (nnz != 0) {
  53.        A=A[inz[1]:inz[nnz]];
  54.        r=zeros(n-inz[nnz],1);
  55.    }
  56.  
  57. // Compute the polynomial roots by forming the Companion Matrix
  58. // and then finding the eigenvalues of it.
  59.  
  60.    n=max(size(A));
  61.    c=diag(ones(1,n-2),-1);
  62.    if (n > 1) {
  63.        c[1;]=-A[2:n] ./ A[1];
  64.    }
  65. // Form the roots by adding the eigenvalues of the companion matrix
  66. // to the zero roots previously found.
  67.    r=[r;eig(c).val];
  68.  
  69.    return r
  70. };
  71.