hopp@nbs-amrf.UUCP (Ted Hopp) (02/02/85)
Newsgroups: net.math Subject: Re: real roots of polynomial References: <7942@brl-tgr.ARPA> From: gwyn@brl-tgr.ARPA (Doug Gwyn <gwyn>) >I would appreciate constructive suggestions concerning fast, stable >methods for determining the real zeroes of a polynomial with real >coefficients. I know of stable but slow methods and of fast methods >that occasionally run amok. The polynomial in question is of degree >four, if that helps; more general solutions might be useful also. >For purposes of discussion, denote the polynomial: > P(x) == a[0] + a[1]*x + a[2]*x^2 + a[3]*x^3 + a[4]*x^4 >where a[i] are all real. What is wanted is one or more real values >of x such that P(x) = 0. Polynomials of degree 4 or less have algebraic solutions. Any good handbook like CRC will show you how to do it. The solution below came from my Schaums's Outline Series "Mathematical Handbook of Formulas and Tables" by Murray Spiegel, 1968. Unfortunately, there are some problems you have to watch out for in trying to implement code for this solution. If two (real) roots are very close, or if a complex root is very close to the real line, numerical accuracy problems rear their nasty heads. It is easy for code to decide that there are no real solutions when there are some, or vice versa. (The first time I naively used code for the algorithm below to solve (x-1)^4 = 0 on a Vax, it decided that all four roots were complex. A typical root was (1 + i*10*-8).) One possibility is to use the solution below to find points known to be close to roots and use an iterative procedure that works in the complex plane, starting with those as a guess. Since the guess will be off only by the numerical errors in computing the formulas below, this should be extremely fast for any reasonable iteration algorithm. Algebraic solution for degree-4 polynomials First normalize P(x) so the coefficient of x^4 is 1: P(x) == x^4 + a*x^3 + b*x^2 + c*x + d Let y1 be a real root of the cubic equation * y^3 - b*y^2 + (ac-4d)y +(4bd-c^2-a^2d) = 0 The roots of P(x) are the four roots of ** z^2 + {a +- sqrt(a^2-4b_4y1)}*z/2 + {y1 +- sqrt(y1^2-4d} = 0 If all roots of (*) are real, pick y1 to be the root that gives all real coefficients in (**) Algebraic solution for degree-3 polynomials To solve a cubic equation x^3 + a*x^2 + b*x + c = 0: Let Q = (3b-a^2)/9 R = (9ab - 27c -2a^3)/54 S^3 = R+sqrt(Q^3+R^2) T^3 = R-sqrt(Q^3+R^2) The roots are given by: x1 = S + T - a/3 x2 = -(S+T)/2 -a/3 + i(S-T)*sqrt(3)/2 x3 = -(S+T)/2 -a/3 - i(S-T)*sqrt(3)/2 (i = sqrt(-1)) Let D == Q^3 + R^2. If a, b, c are real, then i. one root is real and two complex conjugates if D>0 ii. all roots real and at least two are equal if D=0 iii. all roots are real and unequal if D<0. If D<0, let cos(h) = -R/sqrt(Q^3); then the roots are given by: x1 = 2*sqrt(-Q)*cos(h/3) x2 = 2*sqrt(-Q)*cos((h+pi)/3) x3 = 2*sqrt(-Q)*cos((h+2*pi)/3) -- Ted Hopp {seismo,umcp-cs}!nbs-amrf!hopp