Xref: utzoo sci.math:3018 comp.graphics:1987
Path: utzoo!mnetor!uunet!husc6!tut.cis.ohio-state.edu!uwmcsd1!ig!agate!jiff!greg
From: greg@jiff.berkeley.edu (Greg)
Newsgroups: sci.math,comp.graphics
Subject: Re: Solution of quartic equation
Message-ID: <7662@agate.BERKELEY.EDU>
Date: 14 Mar 88 23:34:30 GMT
References: <1656@thorin.cs.unc.edu>
Sender: usenet@agate.BERKELEY.EDU
Reply-To: greg@jiff.UUCP (Greg)
Organization: U. C. Berkeley
Lines: 63
Keywords: Sturm chains, Newton's method

In article <1656@thorin.cs.unc.edu> turner@unc.cs.unc.edu (Douglass Turner) writes:
>Hello, does anyone know of a good, robust, method to find all roots of a 
>quartic equation?

There are various plans of attack.  I assume, since you are ray-tracing,
that you want the real roots, or perhaps even the real positive roots.
I will describe one reasonable method below.

>I have heard that direct solution is prone to numerical instability, so
>an iterative technique is called for.

The direct solution is prone to programmer's instability.  It's very
messy and doesn't generalize.  The most painful step is taking the
cube root of a complex number, which cannot be expressed in terms of
roots of real numbers.  In any case, the numerous roots in the formula
are themselves computed by iterative algorithms in the math package,
so you aren't really saving much by using the exact formula.

My algorithm is based on Newton's method and Sturm chains.  A Sturm
chain is an easy way to find out how many roots a polynomial has
in an interval.  You take the polynomial p(x), which we may call
p0(x), and you take its derivative, p1(x) = p'(x).  Now let pn(x)
be the remainder obtained when you divide p{n-2}(x) by p{n-1}(x).
Quit when you obtain a constant polynomial (the degree goes down
by at least one every time).  The sequence of polynomials p0,...,pk
is the Sturm chain.  Now let S(a) be the number of sign changes
in the sequence p0(a),p1(a),...,pk(a), i.e. the number of times
that pi(a) is positive and p{i+1}(a) is negative, or vice versa.
Then the number of roots in the interval [a,b] is S(a)-S(b).

One can use Sturm chains directly in an enlightened divide-and-conquer
approach, but if you want to zero in on the solution quickly it's best to
switch to Newton's method.  The trick with Newton's method is to avoid
the inflection points.  They are the roots of the second derivative,
p''(x).

Find the real roots of p''(x), say y1,...,ym, by recursion.  In your
case p(x) is a quartic, so you can use the quadratic formula for p''.
Pick y0 to be smaller than all of the roots, and y{m+1} to be larger
than all of the roots.  Now use Sturm chains to find the number of
roots in the intervals [y0,y1], [y1,y2],...,[ym,y{m+1}].  Each interval
has at most two roots.  Use Newton's method in each interval that has
roots in it.  If an interval has two roots, Newton's method will
converge to the two roots starting from each end of the interval.  For
an interval with one root, check if p(x) is concave up or concave down
in that interval (the concavity of course alternates with successive
intervals).  If it is concave up, Newton's method will converge
starting from the endpoint at which p(x) has a positive value.
Otherwise start from the other endpoint.

Since you are ray tracing, you may often only want the first
positive root.  You can discard and truncate the intervals appropriately,
and then progress from left to right, stopping at the first root you
find.

Sturm chains will fail if p(x) and p'(x) share roots, which happens
exactly when p(x) has multiple roots.  Exact numerical coincidences are
a delicate matter in floating point arithmetic, if the first constant
polynomial in the Sturm chain is the zero function, then the previous
polynomial in the chain is the greatest common divisor of p(x) and p'(x),
and you can just divide out p(x) by this and start over.
--
Greg