z2 = -4.2169 9999 9999 9 E17 (1.68E38, 3.72E39)
z3 = 1085.8557 4101 61 (-5E10, 4.9E10)
where the numbers in parentheses are y1(x-e) and y1(x+e). This method returns the correct first and
third roots, but the second root fails the sign test.
If you save your polynomials in coefficient list form and evaluate them with polyEval(), it is easy to use
the list reversal method in tip [3.4] to reverse the coefficients. Bhuvanesh Bhatt's MathTools package
also include a list reversal function, reverse().
The following function prooti() can be used to find the roots with the inversion method.
prooti(c)
Func
©(coefList) Poly roots by inversion
©11jun02/dburkett@infinet.com
local ï
seq(c[ï],ï,dim(c),1,⁻1)→c © reverse coefficients
return 1/(zeros(polyeval(c,ä),ä)) © solve for roots and invert
EndFunc
Method 2: Laguerre's algorithm
Another alternative is to use a completely different algorithm to find the polynomial roots. solve() may
recognize a polynomial and so use a specific algorithm to find its roots, in any event, it doesn't work for
our example, and Laguerre's algorithm does. The function shown below, polyroot(), combines the code
for the programs laguer() and zroots() from Numerical Recipes in FORTRAN. Call polyroot() with the
polynomial coefficients as a list:
polyroot({an, ... a0})
and the roots are returned as a list. You must set the Complex Format mode (with [MODE]) to
RECTANGULAR before running polyroot(), since it uses complex arithmetic even for real roots. Also
set the Exact/Approx mode to APPROX, or polyroot() will take an extremely long time.
polyroot() may return the string "polyroot too many its" if a solution does not converge in 80 iterations
or less, for each root. The call for our example is
polyroot({1,4.217E17,-3.981E20,-6.494E22})
which returns these roots in about 7 seconds:
{-4.2170000000001,1085.855741 0161,-141.81969643465}
All of these roots pass the sign test, and the coefficients are reconstructed to, at worst, 5 least
significant digits.
Some comments on the operation of polyroot() follow this code listing.
polyroot(a)
Func
©({an,...,a0}) return roots of polynomial
6 - 119
To Next Page
Loading...
Need help?
General
