gwyn@brl-tgr.ARPA (Doug Gwyn <gwyn>) (02/01/85)
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.
hopp@nbs-amrf.UUCP (Ted Hopp) (02/02/85)
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
ashby@uiucdcsp.UUCP (02/07/85)
One reliable and robust means of computing the roots of a polynomial is to find the eigenvalues of the polynomial's companion matrix. First set up the companion matrix; see any linera algebra text for a definition. Then just compute that matrix's eigenvalues, via EISPACK for example. Good luck!
jfh@browngr.UUCP (John "Spike" Hughes) (02/11/85)
Why not use the quartic formula (which is a lot like the quadratic formula,
only more complicated)? It can be found in most books on the theory of equations,
and in Jacobson's (abstract) Algebra book, I think. It may even be in the
handbook of mathematical tables, published by the chemical rubber company
(CRC) which is on the refernce shelf in most libraries. If you want the
'C' code for it, write to me and I'll mail it to you.
-jfh{ashby@uiucdcsp.UUCP (02/12/85)
A problem with simply coding a CRC formula is that it may be numerically unstable. I believe finding the roots of the polynomial as the eigenvalues of the companion matrix is numerically reliable as well as easy to code -- if you happen to have EISPACK around.
lindley@ut-ngp.UUCP (John L. Templer) (02/13/85)
From ashby@uiucdcsp.UUCP: > One reliable and robust means of computing the roots of > a polynomial is to find the eigenvalues of the polynomial's > companion matrix. > > First set up the companion matrix; see any linera algebra > text for a definition. Then just compute that matrix's > eigenvalues, via EISPACK for example. I'm curious; does this method fufill the requirement that it be efficient? I thought finding the eigenvalues of a matrix was not all that simple a process. -- John L. Templer University of Texas at Austin {allegra,gatech,seismo!ut-sally,vortex}!ut-ngp!lindley "and they called it, yuppy love."
ashby@uiucdcsp.UUCP (02/20/85)
EISPACK is an easy-to-use and, I believe, efficient way of computing the egvals of a "moderately sized" matrix. Surely a 4x4 matrix is no problem. However, using the CRC quartic formulae may cause trouble.