Solving higher order polynomials

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

About these ads
This entry was posted in Excel, Maths, Newton, UDFs, VBA and tagged , , , , , , . Bookmark the permalink.

6 Responses to Solving higher order polynomials

  1. Sergio says:

    Great! I started using your code right away. I have a nonlinear system of equation which can be reduced to a quartic equation in a single unknown, and both Quartic and PolyRJT gave accurate and fast results. Since the overall computational time is similar, in my case, and you warn that PolyRJT can be more accurate, I’m using it instead of Quartic, but both worked really fine. I’m extremely grateful, and booked your website for future reference.

  2. Pingback: Ping Back Ping Backs | Newton Excel Bach, not (just) an Excel Blog

  3. Pingback: Daily Download 19: Solving polynomials | Newton Excel Bach, not (just) an Excel Blog

  4. Pingback: what Microsoft think VBA is good for … | Newton Excel Bach, not (just) an Excel Blog

  5. Kostas says:

    Great function! Solved my problem. Many many many thanks

  6. Pingback: Solving Quadratic, Cubic, Quartic and higher order equations; examples | Newton Excel Bach, not (just) an Excel Blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s