Date: 6/20/96 at 15:36:45 From: Ramon Handel Subject: Polynomial Roots Dear Doctor Math: I am Ramon van Handel, a high school student in Amsterdam, the Netherlands. I have been interested lately in finding the roots of polynomials. I have found several methods to find their roots in books. First of all, for polynomials of the second degree: AX^2+BX+C = 0 -B +/- SQRT(D) X = ---------------- 2A D = B^2-4AC And Newton Raphson's method: f(x) x (n+1) = ----- f'(x) but neither of these methods is suitable for my needs. The first works only up to the 2nd degree, and the second doesn't give all the roots, just the one that is closest to your first approximation. Is there a reliable method to do both? Yours truly, Ramon van Handel RVHANDEL@DDS.NL
Date: 6/20/96 at 17:39:10 From: Doctor Anthony Subject: Re: Polynomial Roots Solving cubics and quartics is a couple of orders of magnitude more difficult than solving a quadratic. The 5th degree equation cannot be solved by any normal analytical process. In the case of a cubic, if one real root (x1) is known, then dividing by (x-x1) will reduce the cubic to a quadratic which can be solved in the usual way. For the general cubic you must first do a transformation to reduce the equation to the form z^3 + 3Hz + G = 0. (1) This follows by putting z = ax+b or x = (z-b)/a in the equation ax^3 + bx^2 + cx + d = 0. We then find the roots of the z equation as follows: Let z = m^(1/3) + n^(1/3) and by cubing we obtain z^3 - 3m^(1/3)n^(1/3) - (m+n) = 0 Compare with (1) and we see that m^(1/3)n^(1/3) = -H and m + n = -G Therefore m and n are roots of t^2 + Gt - H^3 = 0 We may take m = (1/2)(-G+sqt(G^2 + 4H^3)) and if Q denotes any of the three values of m^(1/3) then the three values of m^(1/3) are Q, wQ, w^2(Q) where w = imaginary cube root of 1. Also m^(1/3)n^(1/3) = -H and so corresponding values of n^(1/3) are -H/Q, -w^2H/Q, -wH/Q Hence values of z i.e. of ax+b are Q-H/Q, wQ-w^2H/Q, w^2(Q)-wH/Q This is known as Cardan's Solution, though it was originally given by Tartaglia. He unwisely told Cardan, who promptly published it as his own. Note that if G^2 + 4H^3 < 0 then m is imaginary and this method is not satisfactory. It means the three roots of the z equation are all real, and a trig. solution is preferred, using the fact that cos(3A) = 4cos^3(A) - 3cos(A) This gives cos^3(A) - (3/4)cos(A) - (1/4)cos(3A) = 0 Let z = qcos(A) so that our z equation becomes cos^3(A) + (3H/q^2)cos(A) + G/q^3 = 0 Comparing the two cubics in z and cos(A) we see that q = 2*sqt(-H) and cos(3A) = -4G/q^3 = -G/[2sqt(-H^3)] Having found cos(3A) we can find 3A and then A. Then z = qcos(A), qcos(2pi/3 + A) and qcos(2pi/3 - A) The above working should be enough to discourage you from analytical solutions of cubics. Quartics are even worse, and quintics are, as I mentioned earlier, impossible in algebraic terms using the coefficients of the equation. Incidentally, you misquoted the Newton Raphson formula. It should read x(n+1) = x(n) - f(x(n))/f'(x(n)) If you do a curve sketch or evaluate f(r) where r = (say) -5, -4, ...3, 4, 5 or whatever seems appropriate and check for sign changes in f(r), then you can bracket the positions of the real roots. As you can see, polynomials are not easy, and computer solutions are the only practical method for equations of high degree. -Doctor Anthony, The Math Forum Check out our web site! http://mathforum.org/dr.math/
Date: 01/28/2003 at 01:58:29 From: Andrew Korsak Subject: Polynomial Roots There exists a very reliable algorithm I use to find roots of any polynomial, even with complex coefficients and multiple roots. I first read about it back in the 1960's in the CACM journal. It is so simple that it seems unbelivable, but it works! Actually, my colleagues at Stanford Research Institute (now renamed SRI International) in Menlo Park, California, discovered that this algorithm was first introduced by Weierstrass in 1903 as the second phase of a two phase proof of the "fundamental theorem of algebra", i.e. that every n-th degree polynomial has n roots if one allows complex numbers. We published a note about this in the SIAM Review Vol. 18, No. 3, July 1976, p.501. Here is the algorithm: (It is assumed that the x^n coefficient is factored out and is 1.) Start with distinct trial roots x_i, i = 1....n. A good logical starting set of roots lie on the unit circle in the complex plane: x_k = e^(2pi*i *k/n) where i = sqrt(-1). At each iteration, let x'_i be the next set of trial roots. The iteration is: For i = 1 to n x'_i = x_i - P(x_i) divided by: (product for j = 1..n, j not = i, of:)(x_i - x_j) Next n The algorithm works regardless of whether the trial roots are updated "on the fly", i.e. x_i is replaced by the new x'_i before going on to i+1, or whether the x_i's are replaced by the x'_i's before proceeding to the next iteration. Our SIAM article asked the world of mathematics to see if anyone can explain why this works, but people submitted proofs only for n=2. Mike Green, one of my co-authors at SRI at the time, told me recently that he actually proved it for n=3 and maybe he thought he had done it for n=4 as well, but I have never heard any more results in this regard. Weirestrass' proof involves an initial step using typical epsilon delta arguments, leading to an initial set of trial roots from which the above algorithm is guaranteed to converge. The interesting thing is that this "phase two" always works, without that "phase one" of Weierstrass. Theoretically, it seems that one could arrive at a polynomial and trial roots for which the iteration process would "lock up" on some set of trial roots and never move from there, or it might just "thrash around" and never converge, but none of us experimenting with it could ever find a case of non- convergence. It is known that the "block update" form (as opposed to the "on the fly" update form earlier described) will "blow up" instantly when two trial roots become equal. Presumably there is some nonlinearly distorted analogy to this for the "on the fly" version. For example, in the quadratic case, this happens if the trial roots are equally spaced away from the center of the right bisector of the true roots in the complex plane. However, such a situation amounts to a set of measure zero in the n-dimensional complex space of trial points. The "bottom line" is that numerically this is never a problem -- any reasonable practical implementation of the algorithm will obviously kick in a fudge factor to move away from a trial roots position where numbers get too large to multiply or too small to divide by. This is what has been found but needs to be further verified empirically, and of course it would be nice to prove theoretically that the algorithm can never arrive at a situation where, for example, some set of trial roots arrived at will never change, i.e. they are permuted among each other but any small numerical perturbation will just bring you back to that same set of wrong trial roots. To put it in more precise mathematical terms, the question is: can one find a polynomial and a set of trial roots forming an open set in complex n-space such that the algorithm will cycle somehow among trial roots within that open set, never converging to the true roots outside of that open set? I did a lot of experimenting using MathCad. I also wrote a program using FORTH (Win32Forth at forth.org/compilers) which plots the progress of the trial roots as they simultaneously converge by this method. I can provide this to any interested parties. Andrew J. Korsak, Ph.D. Mathematician and Software/Firmware Engineer (retired)
Search the Dr. Math Library:
Ask Dr. MathTM
© 1994- The Math Forum at NCTM. All rights reserved.