This call tests the roots returned by solve():
proot_t1({1,4.217E17,-3.981 E20,-6.494E22},{-4.217E17,-141.81969643465,0})
which returns {true,true,false}, indicating that the first two roots pass the sign test and the last root fails.
Note the proot_t1() also tests for f(x) = 0 for each root, and returns true in that case.
proot_t1() can only be used to test real roots of polynomials with real coefficients. The method used to
perturb the roots for the test does not make sense in the complex plane.
Another method to check polynomial solutions is to expand the roots into a polynomial and compare
the resulting coefficients with the original coefficients. One way to do this is
expand(product(x-{zn, ... z0}))
where {zn, ... z0} is a list of the roots and x is a symbolic variable. Applying this method to the roots
above with
expand(product(x-{0,-4.217E17,-141 .81969643465}))
gives
x
3
+ 4.217E17 $ x
2
+ 5.9805365986492e19 $ x
The two highest-order coefficients are the same, but the original constant term has disappeared and
the coefficient for x has the wrong sign, let alone the right value. By now you should be starting to
suspect that solve() has failed us.
You might be tempted to try factor() or cFactor() to find the roots, but they fail in the same way as
solve(), returning (x+0) as one of the factors.
The following sections show a few methods to find the correct result for this root:
Method 1: Invert polynomial coefficients (routine prooti())
Method 2: Laguerre's algorithm (routine polyroot())
Method 3: Eigenvalue methods (routine proots())
These methods have trade-offs in speed, accuracy and code size, as well as their abilities to find all
the real and complex roots.
Method 1: Invert polynomial coefficients
An interesting property of polynomials is that the roots map onto the reciprocals of the roots of the
polynomial with the coefficients reversed, that is, if
and
f
1
(
x
)
= a
3
x
3
+ a
2
x
2
+ a
1
x + a
0
f
2
(
x
)
= a
0
x
3
+ a
1
x
2
+ a
2
x + a
3
and the solutions to f
2
(x) are z
1
, z
2
and z
3
, then the solutions to f
1
(x) are 1/z
1
, 1/z
2
and 1/z
3
. Sometimes
solving this 'reversed' polynomial gives better solutions than the original polynomial. Using this method
and zeros() with our example polynomial, we get
z1 = -141.8196 9643 465 (7E9, -3E9)
6 - 118