home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 4
/
DATAFILE_PDCD4.iso
/
languages
/
rlab1_23a
/
CTB
/
roots
< prev
next >
Wrap
Text File
|
1995-11-14
|
2KB
|
71 lines
//-------------------------------------------------------------------------
//
// roots
//
// Syntax: r=roots(a)
//
// This routine finds the roots of the polynomial given by the
// coefficients in the input vector a. If a has N components,
// the polynomial is given by,
//
// polynomial = a[1]*X^(N-1) + a[2]*X^(N-2) + ... + a[N-1]*X + a[N]
//
// Ref: Golub, G., and Van Loan, C., Matrix Computations, Second Edition,
// The John Hopkins University Press, Baltimore MA, 1989, p. 369.
// (Note: The method in Golub and Van Loan was modified to be
// compatible with the Transfer function notation - which
// is also compatible with the MATLAB version).
//
// Copyright (C), by Jeffrey B. Layton, 1994-95
// Version JBL 950105
// Bill Lash -- Bug Fix for initializing n
//-------------------------------------------------------------------------
rfile isvec
roots = function(a)
{
local(n,inz,nnz,A,c,r)
// Check if input is a vector
if (!isvec(a)) {
error("ROOTS: Input must be a Vector.");
}
// Make sure input is a row vector (if it is a row vector it is a
// little easier to create companion matrix).
if (a.nr > 1) {
A=a.'; // Make sure it's a row vector
else
A=a;
}
// Strip leading zeros. Strip trailing zeros, but keep them as roots
// at zero (store in roots vector r).
n=max(size(A));
inz=find(abs(A));
nnz=max(size(inz));
if (nnz != 0) {
A=A[inz[1]:inz[nnz]];
r=zeros(n-inz[nnz],1);
}
// Compute the polynomial roots by forming the Companion Matrix
// and then finding the eigenvalues of it.
n=max(size(A));
c=diag(ones(1,n-2),-1);
if (n > 1) {
c[1;]=-A[2:n] ./ A[1];
}
// Form the roots by adding the eigenvalues of the companion matrix
// to the zero roots previously found.
r=[r;eig(c).val];
return r
};