home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 4
/
DATAFILE_PDCD4.iso
/
languages
/
rlab1_23a
/
CTB
/
poly
< prev
next >
Wrap
Text File
|
1995-11-14
|
2KB
|
78 lines
//----------------------------------------------------------------------
//
// poly
//
// Syntax: c=poly(x)
//
// This routine evaluates the characteristic polynomial contained
// in x. If x is an N by N matrix, then poly(x) will return a row
// vector of length N+1 whose elements are the coefficients of the
// characteristic polynomial (det(lambda*eye(N,N)-a)).
//
// If x is a vector, then poly(x) will return a vector whose elements
// are the coefficients of the polynomial whose roots are the elements
// of x.
//
// For example,
//
// > rfile poly
// > A=[0,4,5,6];
// > poly(A)
// 1 -15 74 -120 0
//
//
// > rfile poly
// > A=[1,2,3;4,5,6;7,8,0];
// > poly(A)
// 1 -6 -72 -27
//
// In general, the routines poly and roots are the inverse of each
// other (But only for vector inputs).
//
// Ref: MATLAB Reference Guide, The Mathworks Inc. August 1992.
//
// Copyright (C), by Jeffrey B. Layton, 1994-95
// Version JBL 950105
//----------------------------------------------------------------------
poly = function(x)
{
local(n,z,i,c)
// Check for the empty case
if (max(size(x)) == 0) {
c=1;
return c
}
// find eigenvalues of x
if (x.nc == x.nr) {
z=eig(x).val;
else if (x.nr != 1) {
z=x.';
else
z=x;
}}
// strip away infinities and NaN
z=z[finite(z).*[1:z.nc]];
n=max([z.nr,z.nc]);
// Initialize c (characteristic polynomial)
c=zeros(1,n+1);
c[1]=1;
// Expand recursion formula
for (i in 1:n) {
c[2:(i+1)]=c[2:(i+1)]-z[i].*c[1:i];
}
// If x was real, then make c real.
if (x.type == "real") {
c=real(c);
}
return c
};