roots

roots of a polynomial

Syntax

x = roots(p)
x = roots(p, 'e')

Arguments

p

a polynomial with real or complex coefficients, or a vector of the real or complex polynomial coefficients in decreasing power order.

'e', algo

a character: the algorithm to be used (default "f", for "fast"):

"f":	The Jenkins-Traub method is used. The polynomial must be real and of degree < 100; otherwise, an error is yielded.
"e":	eigenvalues of the companion matrix are returned.

Description

This function returns in the complex vector x the roots of the polynomial p.

The "f" option corresponds to the fast RPOLY algorithm, based on Jenkins-Traub method.

The "e" option corresponds to method based on the eigenvalues of the companion matrix.

For real polynomials of degree <=100, users may consider the "f" option, which might be faster in some cases. On the other hand, some specific polynomials are known to be able to make this option to fail. For instance, p=poly([1.e300,1.e0,1.e-300],'x'); provokes infinite looping of roots(p,"f")

Examples

In the following examples, we compute roots of polynomials.

// Roots given a real polynomial
p = poly([1 2 3],"x")
roots(p)
// Roots, given the real coefficients
p = [3 2 1];
roots(p)
// The roots of a complex polynomial
p = poly([0,10,1+%i,1-%i],'x');
roots(p)
// The roots of the polynomial of a matrix
A = rand(3,3);
p = poly(A,'x')
roots(p)
spec(A)

The polynomial representation can have a significant impact on the roots. In the following example, suggested by Wilkinson in the 60s and presented by Moler, we consider a diagonal matrix with diagonal entries equal to 1, 2, ..., 20. The eigenvalues are obviously equal to 1, 2, ..., 20. If we compute the associated characteristic polynomial and compute its roots, we can see that the eigenvalues are significantly different from the expected ones. This implies that just representing the coefficients as IEEE doubles changes the roots.

A = diag(1:20);
p = poly(A,'x')
roots(p)

The "f" option produces an error if the polynomial is complex or if the degree is greater than 100.

// The following case produces an error.
p = %i+%s;
roots(p,"f")
// The following case produces an error.
p = ones(101,1);
roots(p,"f")

The following script is a simple way of checking that the companion matrix gives the same result as the "e" option. It explicitly uses the companion matrix to compute the roots. There is a small step to reverse the coefficients of the polynomial ; indeed, "roots" expects the coefficients in decreasing degree order, while "poly" expects the coefficients in increasing degree order.

v= [1.12119799 0 3.512D+13 32 3.275D+27 0 1.117D+41 4.952D+27 1.722D+54 0 1.224D+67 0 3.262D+79 ];
r1 = roots(v,"e"); // With "e" option
dv = size(v,"*");
p = poly(v(dv:-1:1),"x","coeff"); // Reversing v's coefficients
A = companion(p);
r2 = spec(A); // With the companion matrix
max(abs(r1-r2))

References

The RPOLY algorithm is described in "Algorithm 493: Zeros of a Real Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189

Jenkins, M. A. and Traub, J. F. (1970), A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration, SIAM J. Numer. Anal., 7(1970), 545-566.

Jenkins, M. A. and Traub, J. F. (1970), Principles for Testing Polynomial Zerofinding Programs. ACM TOMS 1, 1 (March 1975), pp. 26-34

Report an issue
<< residu	Polynomials	rowcompr >>