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.
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.
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: