Newton Excel Bach, not (just) an Excel Blog

… and a better solution for quadratics.

Having written a better solver for high order polynomials, that raises a few other issues.  The error in the roots found by the RPolyJT was checked by simply substituting each root as the x value in: ax^n + bx^(n-1) + cx^(n-2) + …+e = 0.  The problem is that finding the difference between two similar large numbers significantly reduces the precision of the result.  A much better approach is known as the Horner Scheme, which uses Synthetic Division to arrive at the value of the polynomial without subtraction of the values raised to high powers.

This method has now been added to the Polynomial spreadsheet with the following changes:

• RPolyJT has been modified to call the new function EvalErrorHC() to evaluate the error in the roots found.
• EvalErrorHC may also be used as a UDF from the spreadsheet.
• EvalPolyHC will evaluate any polynomial for a list of values of x (real or complex), which may take any value (i.e. not necessarily estimated roots).

For similar reasons, the standard quadratic formula will not give accurate results when b^2 and 4ac are large and of similar magnitude.  Better accuracy is given by a method included in the Fortran code used as the basis of the RPolyJT function.  The Quadratic function in the Polynomial spreadsheet has now been revised to use this method.

The revised spreadsheet may be downloaded (including full open source code) from:

Polynomial.zip

Results of the new RPolyJT and EvalErrorHC functions are shown in the screen shot below:

Click for full size view

The Quartic spreadsheet presented here previously (most recently here) uses algebraic techniques to provide an “exact” solution to polynomial equations of up to fouth order.  For polynomials of higher order there is no general algebraic solution, and numerical methods are required.  Also under some circumstances the “exact” solution, when evaluated with floating point numbers of limited precision, produces results that are very far from exact.  In these circumstances a well written iterative numerical procedure can give much better results.

One such procedure is the Jenkins-Traub algorithm, which apart from the advantage of one author having an excellent family name, is also said to be both fast and stable.  The Wikipedia article says that it “has been described as practically a standard in black-box polynomial root-finders”.

I have added a VBA user defined function (UDF) implementation of the Jenkins-Traub method to the Quartic spreadsheet, which has been re-named Polynomial.xlsb, and is available for free download, including open-source code. 

Download Polynomial.zip

The VBA code is a translation of Fortran77 code at the NetLib library.  Similar code in Fortran90 format can be found at Alan Miller’s Fortran Software.  Both of these sites are invaluable sources of open source scientific Fortran code.  The chosen routine provides all real and complex solutions to any polynomial with real coefficients.  Fortran routines are also available for polynomials with complex coefficients, but I have not yet translated these.

The screen-shots below show output from the Jenkins-Traub UDF (RPolyJT) for a quartic polynomial, followed by the results from the Quartic UDF for the same input.  Note that the Jenkins-Traub result is close to machine precision, but in this case the Quartic output has errors in the third decimal place.

RPolyJT output for quartic polynomial

Quartic output for the same polynomial

Output from RPolyJT for higher order polynomials is shown in the screenshots below, where the function has been used to solve the equation : 2.x^60 – 1 = 0.  The roots are 60 equally spaced points in the complex plane, with a modulus of 0.5^(1/60).  Note that the input is the range A26:A86; the cells containing zero coefficients could have been left blank, with the same results.  The top and bottom of the input and output results are shown in the two screenshots below:

Top of input and output for 60th order poynomial

Bottom of input and output for 60th order polynomial