Finding real roots of a polynomial

[sbasak]

What is the formula for finding all real roots of a polynomial of the form

a0.x + a1.x + a2.x^2 + ... + an.x^n = 0


Thanks for any help.

[Guv]

The formula for a quadratic (ax^2 + bx + c = 0) is the following.

x1 = (-b + SquareRoot(b^2 - 4ac) / 2a

x2 = (-b + SquareRoot(b^2 - 4ac) / 2a

The formulae for third and fourth order polynomials are so complicated that nobody uses them except students who are taking certain mathematics courses. Even before hand calculators and computers, it was easier to use some successive approximation method.

Over 100 years ago, it was proven that there could be no formula for fifth & higher order polynomials.

Note that even order polynomials might not have any real roots, while odd order polynomials always have at least one real root.

A good successive approximation method is Newton-Raphson, which works as follows (Fourth order polynomial used as example).

Polynomial(x) = ax^4 + bx^3 + cx^2 + dx + e
Derivative(x) = 4ax^3 + 3bx^4 +2cx + d

NextGuess = LastGuess - Polynomial(LastGuess) / Derivative(LastGuess)

Start with a guess. Use the above to compute next guess. Use next guess to compute a better guess. Keep it up until two successive guesses are nearly equal.

If you study the formula for above derivative a bit, you will see the rule for determining the derivative of any polynomial.

If you do complex arithmetic, you can use the above to find complex as well as real roots. If you do complex arithmetic, you can use the above to find complex as well as real roots. Complex roots always come in pairs.

If (a + ib) is a root, then (a - ib) is also a root.

[wossname]

Guv, I just spotted an undocumented feature in your last post.

x1 = (-b + SquareRoot(b^2 - 4ac) / 2a

x2 = (-b - SquareRoot(b^2 - 4ac) / 2a

[riis]

Wossname, I also just spotted in this thread, even in your post

x1 = (-b + SquareRoot(b^2 - 4ac)) / 2a

x2 = (-b - SquareRoot(b^2 - 4ac)) / 2a

Nonetheless I think it was obvious those brackets should be there.

[wy125]

I wrote some code one time for finding the roots of a cubic. I could dig it up if you're interested. You'll probably have to make some slight modifications.

[Guv]

Riis and Wossname: Sloppy posts should result in a reduction of Post Count.

You folks have sharp eyes.

[sbasak]

Yes, I know about Newton Raphson method. However, I like to make a program which can give the real root (if any) of any given polynomial equation. The user will give only the equation and no guess value.

What algorithm do you recommend? If you already have made such program, I shall be glad if you post it here.

Thanks

[Guv]

As mentioned in a previous post, there is no formula for polynomials of order 5 or greater. It was proven over 100 years ago that such formulae cannot exist.

For 3rd & 4th order polynomials, there are algorithms leading to roots, real or complex. The formulae derivable from these algorithms are so messy that nobody ever uses them. I do not think they are published anywhere.

Fighting your way through the algorithms to get a solution is bad enough. If you really want to investigate, try the following site.

http://mathworld.wolfram.com/CubicEquation.html

If you check out that site, you might decide you do not want a formula for even the 3rd order polynomial.

Why not program Newton-Raphson? It is as good or better than any other method, and not difficult to program for a general polynomial. I have written programs using that method in several languages, including capabilities for getting complex roots. I have also used it many times, doing the computations on a hand calculator. My current calculator has a built in polynomial roots function.

There is a pitfall or two when using Newton Raphson, but they are not difficult to avoid.

I think there might be iterative methods for finding quadratic factors, which can then be solved using the formula for quadratic roots. Perhaps you might be able to find a method like this if you do not like Newton Raphson.

[wy125]

This is the code I wrote for getting the roots of a cubic polynomial. As I recall, I was only interested in real positive roots.


Public Function getCubicRoots(coefs() As Double) As Roots
'This function computes the roots of a cubic equation with coeficients contained in coefs():
'The first element of coefs() is the coefficient of the cubic term, the second element is the coeficient of the
'quadratic term, and so on with the last term being the constant term.
'So that if we have an equation of the form: aX^3 + bX^2 +cX + d
' coefs(0) = a
' coefs(1) = b
' coefs(2) = c
' coefs(3) = d



Dim e As Double, f As Double, g As Double
Dim p As Double, q As Double, R As Double
Dim Y1 As Double, A As Double, B As Double
Dim X As Double
Dim w As Double, theta As Double
Dim Y(2) As Double
Dim root As Roots


Dim z As Integer



'put the equation in the form X^3 + eX^2 + fX + g = 0
'to do this we divide all coeficients in the coefs() array by the cubic term ('a')

e = coefs(1) / coefs(0)
f = coefs(2) / coefs(0)
g = coefs(3) / coefs(0)


p = 1 / 3 * (3 * f - e ^ 2)

q = 1 / 27 * (27 * g - 9 * e * f + 2 * e ^ 3)
R = (p / 3) ^ 3 + (q / 2) ^ 2
If R > 0 Then
A = (-q / 2 + Sqr(R)) ^ (1 / 3)
B = (-q / 2 - Sqr(R)) ^ (1 / 3)
Y1 = A + B
X = Y1 - e / 3 'the cubic was reduced by substition, cH = y - e/3

ElseIf R < 0 Then

w = Sqr((q ^ 2 / 4) / (-p ^ 3 / 27))
'Visual Basic does not have an arc cosine function
'so the expression below is used as a substitute
theta = Atn(-w / Sqr(-w * w + 1)) + 2 * Atn(1)

'calculate the roots
For z = 0 To 2

If q > 0 Then

Y(z) = -2 * Sqr(-p / 3) * Cos(theta / 3 + 2 * PI / 3 * z)
root.root(z) = Y(z) - e / 3

ElseIf q < 0 Then
Y(z) = 2 * Sqr(-p / 3) * Cos(theta / 3 + 2 * PI / 3 * z)
root.root(z) = Y(z) - e / 3

End If
Next z


ElseIf R = 0 Then
ErrorHandler (2)

End If


getCubicRoots = root




End Function

[Guv]

Wy125: Your code looks like an implementation of the Tartaglia/Cardano solution.

Have you tested it?

What does it do for a cubic with three real roots?

What happens with only one root?

[sbasak]

Well, if I use Newton Raphson formula, I always need a guess value to start with. Now my question is - how to find the guess value?

One option might be:

x = a very large value
find f(-x) and f(x)
if f(-x) < 0 and f(x) > 0 then root lies between them,
so start with guess value of x
else
set x = x/2
repeat the above method

The problem is, once I find a root, how do I venture for the 2nd or 3rd root of the equation?

For example, how do Newton Raphson method find all real roots of following equations without having to specify a guess value for each root?

x^2-5*x+6=0 or (x-1)*(x-2)*(x-3)=0 ?

Any help?

[wy125]

Guv:

I don't know who developed the solution that I used in my code. I got it from Perry's Chemical Engineers' Handbook.

The code I posted only worked for the case when R<0, which was when there were 3 real unequal roots.

I worked on it this morning and added the functionality to handle the other cases; R = 0 (three real roots, of which at least two are equal), R > 0 (one real root and two conjugate complex roots).

I tested (not extensively) the latest version against mathCAD's polyroots function and it worked for few examples I tried for each of the three states of R.

If anyone is interested, I could post the new code that now returns all the roots (complex and real).

[Guv]

Sbasak: A program implementing Newton-Raphson for finding roots of a polynomial requires a polynomial evaluation Function. The polynomial evaluation function can be used to find initial guesses for roots.

It might be worthwhile to check a good text. Theory of Equations might be the discipline dealing with polynomial roots. A lot of theorems have been proved about roots. Some of those theorems might be useful.

If the polynomial has been adjusted so that the highest order coefficient is +1, the following are valid statements.

• The constant term is the product of the roots: As is for even order polynomials; Reverse sign for odd order polynomials. Perhaps the nth root of the constant term is a worthwhile first guess. Perhaps this is a good initial value for Delta (See below).
  
• The negative sum of the roots is the second highest order coefficient. Perhaps 1/n times the sum is a worth while first guess. Perhaps this is a good initial value for Delta (See below).
  
• You could write a For x = 0 to 10000 Step Delta . . .Next loop, which evaluates the polynomial, checking for sign changes. There is always a real root between values for which the sign changes. Similarly do a loop For x = 0 to -10000 Step -Delta. This method can miss two real roots which are very close to each other, so you might want to start with a large value for Delta, and then repeat the loop for smaller values (Say Delta/5) if you do not find any sign changes.
  
• You might want to use different algorithms for odd & even order polynomials. Odd polynomials always have one real root. Using a large positive value for the first guess is likely to result in convergence to a root.

Unfortunately, complex roots can cause Newton-Raphson to get into an infinite loop when you are trying to find real roots. It is not a high frequency problem, but an application must have logic to avoid such a loop. A simple iteration counter should do the job. Newton-Raphson should almost always converge in 10-20 iterations if it is going to converge at all. Quitting after 50-100 iterations will almost never result in missing a root. You could have an infinite loop problem starting with a particular guess, while a different initial guess will result in convergence.

Newton-Raphson will converge to complex roots if you do complex arithmetic.

BTW: My HP hand calculator can find all real & complex roots of a polynomial in a few seconds. Since it has nothing like the power and capacity of a PC, it should not be difficult to write a root finder algorithm in VB.

[sbasak]

Thanks for your elaborate reply. I too have the HP Graph calculator which calculates all roots of polynomial. For the very reason, I wondered how to do it! In face, in MathCad too, I found that they can show all roots of the equation. But they don't tell how to do it!

I have a program which can evaluate an expression eg. 2*3-5 = 1
So, I can evaluate an equation pretty easily.

I've searched on the net and even some text books on finding polynomial roots. But unfortunately most numerical analysis books only shows very basic example of finding only one root using Newton Raphson or equivalent methods!

Do you have any such program (with code) which can perform this task?

Thanks a lot.

[sbasak]

Could you please post the code which finds all 3 roots of a cubic equation?

Thanks.

[sbasak]

f(x) = a0 + a1.x + a2.x^2 + a3.x^3 + ... + an.x^n = 0

The coefficients and power of x is stored as input matrix a[n] where

by formula all roots should lie in the region of

LowerLimit = -Sqr((a(n - 2) / a(n - 1)) ^ 2 - 2 * (a(n - 3) / a(n - 1)))
UpperLimit = Sqr((a(n - 2) / a(n - 1)) ^ 2 - 2 * (a(n - 3) / a(n - 1)))

tabulate f(x) in the above zone with delta step
Note where f(x) changes sign. In all those zone, there should be a root.

Using that guess value, use Newton Raphson method to find out a root.

Venture for next root.

This method answers to my first question - how to get a suitable guess value for NewtonRaphson method to start with.

This formula, helps to get the first guess value of Newton Raphson method. Once a root is found, we can proceed for 2nd root where f(x) crosses x-axis again.

However, my question is, if I proceed in this way, why do I need to apply NewtonRaphson method at all? Since to find out where
f(x) crosses x-axis by above method, I already get the value for which f(x)=0, so I can get all the roots in this manner.

So, what do you think about it?

[sbasak]

I finally made the working model for finding real roots of polynomial. Though it's not perfect yet - but works anyway!

The file is attached.

:-)

[Guv]

Sbasak: I downloaded and tried your root finder application.

For the default cubic, it only gave about 4-5 digits of precision. VB should be able to get at least 12 digits, if not 14-15

A few equations I tried could not be solved. Once I was erroneously told that zero was a root, after getting a messge about slope being too small.

Why require the user to enter the derivative? It is easy for a program to derive.

I only tried cubics. what is the limit for order?

[sbasak]

It solves polynomial of the form:
a0 + a1*x + a2*x^2 + ... + an*x^n = 0

Input: a0,a1,a2,...,an (as string - values separated by comma)
Output: all real roots

It tries to find root zone by itself, if fails, it asks to specify a limit.

For more examples, see source code. (attached zip file)

Thanks.

[Guv]

Sbasak: I downloaded and tried your second program on

x^3 - 6x^2 + 11x -6

The coefficients were entered as: -6, 11, -6, 1

Result was: Error 48 File not found calcdllw.dll

Click on Okay button resulted in: Error finding roots.

I have VB 6.0 with SP5 running on a Windows 98SE system.

[sbasak]

I can't understand why you did get an error (file not found)!

Please make sure that after deflating the zip files, you keep them all 4 files in same folder.
Files are: polynomial.vbp, polynomial.frm, polynomial.frx and calcdllw.dll

FYI: I made the program in VB5 on Windows 98 and it perfectly worked!

for your problem, x^3 - 6x^2 + 11x -6 = 0
give input as: -6,11,-6,1
It will give 3 real roots as: 1, 2, 3 (some rounding of might happen)

Hope this will solve your problem.

If you still get dll file not found error, try changing the line

Private Declare Function expr Lib "calcdllw.dll" (ByVal x As String, ByRef y As Double) As Integer

to as

Private Declare Function expr Lib "c:\your file path\calcdllw.dll" (ByVal x As String, ByRef y As Double) As Integer

Thanks.

[Guv]

Sbasak: Made the suggested change, and the dll was found. There is probably a way to make the program use app.path to find such files, so that compiling & runnin on another system will avoid this type of problem.

The application worked, although I do not understand why it only gets 4-6 digits of precision.

[Guv]

After seeing and using the applications posted here, I decided to make a Root Finder application which comes close to doing the job the way it should be done.

The application was created using VB 6.00 (SP5) under Windows 98. I have no idea what other environments it can cope with.

I have zipped all the source files rather than making a running version of the application. Compile (with or without source tweaking) should allow to allow it to run under other environments.

If any of you download and try it, I would be interested in comments, especially if something goes wrong. Note that it will not cope with multiple roots, making it unsuitable as a commercial package.

Try x³ + 3x² + 3x + 1, which has 3 equal roots, and cannot be handled by the application. I am not sure what other anomalous situations will cause trouble.

This application is user friendly, and it will easily handle polynomials up to order twenty. It will handle higher order polynomials, but the user interface gets confusing.

It will find complex as well as real roots, storing all roots in List Boxes (One List Box for real parts and one for imaginary parts). The two List Boxes have coordinated scrolling.

Precision is pretty good, but documentation is lacking. There is a Readme file which can be opened by MS Word, WordPad, or Corel WordPerfect (the latter is my WP of choice).

There is only one feature which might be difficult to hack without the Readme file.

[Guv]

Oops: I forgot to attach the .zip file to previous post.

[Guv]

I just remembered that a multiple root is also a root of the derivative of the polynomial, casuing Newton-Raphson to get into trouble with 0/0, or close to it.

When a polynomial causes trouble, you try to find roots of the derivative. Then check those roots to see if they are also roots of the polynomial.

If a root is a triple or quadruple root, Newton-Raphson has trouble finding roots of the derivative. If this happens, try to find roots of the second derivative.

[Guv]

Just discovered a bug in my application.

It only bites if the high order coefficient is not one.

[Guv]

Spent some time cleaning up the Root finder application. Some existing bugs were fixed. It will find all roots, real and complex for polynomials of order 20 or more. Not sure of the limit.

It uses the Rnd function to determine initial guesses in the complex plain. So far, this approach has worked well. When a root is found, synthetic division is used to obtain a lower order polynomial. A linear factor is removed for a real root, while a quadratic factor is removed when a complex root is found.

While most of the roots are found using reduced polynomials, Newton-Raphson is applied to the original polynomial to improve the accuracy of each root. This avoids loss of precision due to roundoff errors introduced by synthetic division.

The application includes a polynomial evaluation option. Clicking on a root will evaluate the polynomial and the derivative at that root. This feature was used to verify root finder logic & precision, and is the only feature which might be difficult to figure out.

The sum and product of the roots are displayed as a check on the overall precision of the roots.

[Guv]

Some how the zipped application was not attached to previous post.

[DavidHooper]

I've got to say Guv, that has perhaps one of the worst interfaces I have ever seen. No offence. However, the actual operation of the program is well sweet!

Some suggestions to tidy it up:
*Could get rid of the complex coefficients option
*Could use a simpler method of entering coefficients
*Could generally tidy up the interface
*Could give it a name change

Do persevere, I think this program could be really useful.

[Guv]

DavidHooper: I would be intertested in what sort of interface you would use to enter coefficients and display roots.

[DavidHooper]

Having read back what I wrote, I apologise - it was too strongly worded.

I should have suggested generally providing a more intuitive interface.