The fundamental theorem of algebra says that a polynomial of degreee
**n** with complex coefficients (which includes real coefficients as
a special case, of course) has exactly **n** complex roots, counting
multiplicities. Sometimes one wants to know just what these roots
are. More often, one just wants to know something about the location of
the roots in the complex plane (are they all in , for
example).

Many people only * think* they want roots of polynomials. For
example, to solve an eigenvalue problem they first compute the
characteristic polynomial, (which is a bad idea for computing eigenvalues
but is the first method taught --- often the * only * method taught) and
then ask the hard computational question `What are the roots of this
polynomial?', never dreaming that the computational question `What are
the eigenvalues of this matrix?' is far easier. Similarly using polynomial
roots to specify equilibria of differential equations turns an easy
computational problem into a hard one.

Nevertheless, sometimes you do want polynomial roots, and, after all,
the computations can be interesting. The numerical approximation of
polynomial roots is pursued farther in Part III, but for now we note
that the routine `realroot`

provides guaranteed bounds of
arbitrary tightness for real roots of polynomials. For complex
roots we will do something else later. The routines `sturm`

and `sturmseq`

let us count the number of roots in any given
real interval, also.

> readlib(realroot): > with(linalg): Warning: new definition for norm Warning: new definition for trace > A := hilbert(5): > p := charpoly(A,x); # This polynomial has real roots. 5 563 4 735781 3 852401 2 61501 p := x - --- x + ------- x - --------- x + ----------- x - 1/266716800000 315 2116800 222264000 53343360000 > st := time(): > realroot(p,Float(1,-10)); # This command gets the real roots to 10 places. 28243 56487 656911 5255289 [[----------, -----------], [----------, -----------], 8589934592 17179869184 2147483648 17179869184 195979213 97989607 895647649 3582590597 [-----------, ----------], [----------, -----------], 17179869184 8589934592 4294967296 17179869184 26921725877 13460862939 [-----------, -----------]] 17179869184 8589934592 > time_taken := time() - st; time_taken := 7.000 > evalf(""); -5 -5 [[.3287917934*10 , .3287976142*10 ], [.0003058980219, .0003058980801], [.01140749158, .01140749163], [.2085342186, .2085342187], [1.567050691, 1.567050691]]The routine

`realroot`

returns a list of pairs of dyadic rationals
(rationals whose denominators are pure powers of two --- computing with
these is more efficient than arbitrary rationals). Each pair of rationals
represents an interval on the real line containing at least one root of
the polynomial. In the above example there were five separate intervals
returned, each of width less than , so we have successfully
located all roots to ten decimal places absolutely.
Hereafter we regard the problem of the computation of arbitrarily accurate approximations or guaranteed bounds for roots of polynomials as a `kitchen problem' --- that is to say, a solved problem, one that we can carry out at any time we like. This isn't quite true in Maple --- there are high-degree polynomials that give trouble and we will see these later --- but for the purpose of the following discussion assume that this is so.

If we can compute roots at any time, then in particular we can compute
them `later' and work with a symbolic representation of the roots.
The way mathematicians do this is to use a phrase something like
`Let be a root of the polynomial . Then ...' and
compute formulae in terms of . This is possible in Maple
as well, and often very useful. The syntax is similar to the
above mathematical phrase, as is the meaning. Maple uses the
`RootOf`

construct, as follows.

> z := RootOf(x^5+x+1,x); 5 z := RootOf(_Z + _Z + 1)

> z^5; 5 5 RootOf(_Z + _Z + 1) > simplify("); 5 - RootOf(_Z + _Z + 1) - 1Use the definition to simplify expressions containing the RootOf.

> simplify(1/z); 5 4 - RootOf(_Z + _Z + 1) - 1Those expressions are ugly to look at, so we can use

`alias`

to clean up the presentation.
> alias(alpha=RootOf(x^5+x+1,x)); I, alpha > z; alpha > simplify(z^5); - alpha - 1 > simplify(1/z); 4 - alpha - 1 > simplify((z^3+z)/(1-z-z^5)); 3 1/2 alpha + 1/2 alpha > simplify(alpha^3/(1-alpha)); 4 3 2 1/3 alpha + 1/3 alpha - 2/3 alpha - 2/3 alpha - 1/3The Maple

`RootOf`

is best thought of as `RootOf(p,x)`

to `-RootOf(p,x)`

gives zero.
This means that Maple assumes that every single occurrence of a given
`RootOf`

refers to the * same* unspecified root. Usually this
is the case, but sometimes it isn't and there are bugs due to this
interpretation known in Maple. As of this writing there is no
satisfactory solution to this bug --- we need an acceptable way of
telling the difference between different unspecified roots.

We have already seen some examples of the uses for `RootOf`

.
We saw it used to represent an algebraic function for the purpose
of integration. In fact, `RootOf`

can be used to represent algebraic functions more generally than
radical expressions can. It can also be used to give concise formulae
for eigenvectors in terms of eigenvalues, as follows.

> with(linalg): Warning: new definition for norm Warning: new definition for trace > A := hilbert(3): > p := charpoly(A,x); 3 23 2 127 p := x - ---- x + --- x - 1/2160 15 720 > alias(lambda=RootOf(p,x)); I, lambda > eye := n -> array(1..n,1..n,identity); eye := n -> array(1 .. n, 1 .. n, identity) > B := evalm(A - lambda*eye(3)); [ 1 - lambda 1/2 1/3 ] [ ] B := [ 1/2 1/3 - lambda 1/4 ] [ ] [ 1/3 1/4 1/5 - lambda ] > z := vector(3,0); z := [ 0, 0, 0 ] > eigenvector_corresponding_to_lambda := linsolve(B,z); eigenvector_corresponding_to_lambda := [ _t[1], 154 2 - --- _t[1] + 42 _t[1] lambda - 80/3 _t[1] lambda , 27 2 50/9 _t[1] - 60 _t[1] lambda + 40 _t[1] lambda ] > > > > subs(_t[1]=1,"); 154 2 2 [ 1, - --- + 42 lambda - 80/3 lambda , 50/9 - 60 lambda + 40 lambda ] 27In the example above, we have found all the eigenvectors of the matrix

The following example shows that algebraic functions and perturbation
expansions can also be investigated using `RootOf`

.

> p := x^5 + e*x - 1; 5 p := x + e x - 1 > alias(z=RootOf(p,x)); I, z > alias(z0 = RootOf(x^5-1,x)); I, z, z0 > diff(z,e); 3 4 2 3 e e z z e z - 320 ------------- + 64 ------------- - 625 ------------- - 500 ------------- 5 5 5 5 256 e + 3125 256 e + 3125 256 e + 3125 256 e + 3125 2 4 e z - 400 ------------- 5 256 e + 3125 > series(z,e); 2 3 2 4 3 21 5 6 z0 - 1/5 z0 e - 1/25 z0 e - 1/125 z0 e + ----- z0 e + O(e ) 15625 > convert(",polynom); 2 3 2 4 3 21 5 z0 - 1/5 z0 e - 1/25 z0 e - 1/125 z0 e + ----- z0 e 15625 > subs(e=1/100,"); 156250000000021 2 3 4 --------------- z0 - 1/500 z0 - 1/250000 z0 - 1/125000000 z0 156250000000000This gives approximate formulae for the roots of . These formulae are valid no matter which root we take.

> > evalf(subs(z0=1,")); .9979959920

** Exercises**

- Use
`realroot`

to compute tight bounds on the real roots of and compare with the series solution given above. - Try to integrate
**z**with respect to**e**, where .

Wed Feb 28 18:25:56 EST 1996