Problem Statement:

We are interested in computing the roots of the following polynomial using the bisection method:

f(x) = x⁹ - 18x⁸ + 144x⁷ - 672x⁶ + 2016x⁵ - 4032x⁴ + 5376x³ - 4608x² + 2304x - 512

Polynomial Evaluation using Horner's Rule:

//-------------------------------------------------------------------------- 
// CREATE THE ARRAY OF COEFFICIENT. 
// REMEMBER TO ORDER FROM LOWER TO HIGHER DEGREE. 
//-------------------------------------------------------------------------- 
#define DEGREE 9 
FLOAT coeffs[DEGREE+1] = {-512,2304,-4608,5376,-4032,2016,-672,144,-18,1};

//-------------------------------------------------------------------------- 
// Polynomial Evaluation routine : given a set of n+1 coefficients for a 
// degree n polynomial, the degree n of the polynomial, and a value x to 
// evaluate for. NOTE: a[0] is the degree 0 coefficient, etc. 
//-------------------------------------------------------------------------- 
FLOAT P(FLOAT x, FLOAT *a, int n) 
{ 
  FLOAT p=a[n]; 
  for (n--; n>=0; n--) p=(x*p)+a[n]; 
  return(p); 
} 

Bisection Method for Root-Finding:

#define FLOAT float   // (or double)
#define TOL 0.000001

//--------------------------------------------------------------------------
// Bisection method for root finding: given a bracketed range (x lower to
// x upper), a desired error tolerance, and a function to evaluate).
//--------------------------------------------------------------------------
FLOAT Bisect(FLOAT xl, FLOAT xu)
{
  if ( (xu-xl) <= (2*TOL) )  return( (xl+xu)*0.5 );

  FLOAT mid = (xl+xu)*0.5;
  FLOAT Pmid = P(mid,coeffs,DEGREE);
  FLOAT Plow = P(xl,coeffs,DEGREE);
  if ( (Plow*Pmid) < 0)
    return( Bisect(xl,mid) );
  else
    return( Bisect(mid,xu) );
}

The Main Driver evaluating the Bisection routine 
for different bracketing intervals (50x50 samples):
  given A in {1.96,2.00} and B in {2.00,2.04},
  the intervals are for all {A,B] (in discrete steps).

void main()
{
  FLOAT al=1.96, au=2.00, astep=0.0008,
        bl=2.00, bu=2.04, bstep=0.0008;
  FLOAT A,B;

  printf("dummy ");
  for (B=bl; B<=bu; B+=bstep) printf("%f ", B);
//  for (B=bl; B<=bu; B+=bstep) printf("%1.16lf ", B);  // FOR DOUBLE
  printf("\n");

  for (A=al; A<=au; A+=astep)
  {
    printf("%f ", A);
//    printf("%1.16lf ", A);                            // FOR DOUBLE
    for (B=bl; B<=bu; B+=bstep)
      printf("%f ", Bisect(A,B));
//      printf("%1.16lf ", Bisect(A,B));                // FOR DOUBLE
    printf("\n");
  }
}

Plot of the Polynomial Near the Actual Root:

It is important to note that the erratic behavior of the function exists with both single and double precision, but on different scales. Using single precision, the effective "window" of error is on the scale of +/-10^-3 and with double precision, it is +/-10^-11. What this suggests is that even though the overall magnitude of deviation is greatly reduced using double precision it will not really improve the performance of the bisection method since it relies on bracketing the root between interval endpoints whose evaluations are of opposite signs. Both of these plots show equally erratic sign changes and will therefore cause inaccurate convergence.

Plotted Results of the Bisection Routine over varying intervals:

As predicted earlier, both plots show about equal behavior of the bisection root-finding method regardless of whether or not we use single or double precision arithmetic. However, the double precision plot shows a slightly larger flat region emanating from the front-left of the graph. Nevertheless, the region is hardly worth noting since it will not help the root-finder.

Improving the Accuracy of the Algorithm:

So what do these graphs show about how to improve the accuracy of the bisection method? Not much. It seems that there is no reasonable way to find the root of such an erratic and "flat" function (without knowing about the nature of the function in advance). Other methods would also surely fail since it is difficult to determine if the "zero" obtained is the real root or just an erratic jump around zero. Methods based on the evaluation of the derivative (such as the Newton-Raphson method) would also fail miserably for such a "flat" function. The derivative - or "slope" of the tangent - of the function approaches (and crosses) zero in this flat area causing the root-finder to diverge to infinity. The bisection method is guaranteed to converge, but it is unlikely to converge to the actual roots.

Why does this function behave so erratically? Because this long sequence of terms must sum to zero, the values along the way must at some point be of similar value thus causing loss of most significant digits during the addition and subtraction operations; this is called catastrophic cancellation. Upon inspection of the polynomial, the error seems to surely result from the sequences of multiplications and and exponentials; however, these operations receive values between 1.96 and 2.04 and are unlikely to be a major source of error, despite the repeated multiplication of the fractional parts up to nine times. On the other hand, the addition and subtraction of similar values can have a devastating effect on the accuracy.

Here are some possibilities to improve the accuracy that were tried and unsuccessful:

• Scale the function in y for each iteration of the Horner's Method in hopes of shrinking this interval of erratic behavior. The interval does shrink, but the magnitude of the erratic spikes is also scaled causing a negating effect on improving accuracy.
• Try different stopping criteria. We currently use the interval "width" as a stopping criterion for the bisection method. Perhaps trying interval endpoint evaluation compared against the tolerance values would give better results. This too doesn't help since the interval could be quite wide when the iteration is terminated. The root would just be chosen as the midpoint and there is definitely no guarantee that this erratic interval is symmetric about the actual root. Simply put, this alteration would just try to find the endpoints of the erratic interval, terminate when this is reached, and simply use the midpoint of the region as the root.
• I suggested symbolic computation as a method for finding this root to Andrew Thall, and after he thought about it for a few minutes, he agreed that it seems to provide a method for solving such degenerate cases. Basically (with details provided by Andrew), the technique involves assuming that the function has rational roots, guessing a few trial root values, and then verifying with synthetic division. Each stage can be done algorithmically.
• Rewrite the function in an alternate form. The polynomial in question is the expanded form of (x-2)⁹. If we plot the function using this form with only single precision we obtain the following:



The erratic behavior vanishes!! The reasons are quite clear since we only have one chance of catastrophic cancellation (one subtraction of two similar values near the root), but this will not matter since the worst that can happen is that near the root the result of the subtraction will be very small (most significant digits are cancelled) and after being raised to ninth power the result will surely be zero. Raising these fractional values to a high power, causes the values to converge to zero completely eliminating any chance of erratic behavior. This requires some clarification. The primary difference between this form and the original form is that we are taking a possibly small value and raising it to a high power to reach zero versus summing a series of terms in hopes of hitting a target zero. Raising fractions to powers will clearly reach zero much more "effectively".

The following statements are not quite accurate, but I leave them here for illustration of a possible pitfall:

Does this help though? Not really since the bisection method will still terminate quite far from the actual root. Even though the graph looks much smoother, it probably isn't much more accurate since there are now MANY values of zero across the entire region that was previously erratic. However, like before, this ill-conditioned region will be reached and, since the stopping criterion is based on interval width, the "fall-through" path will be taken and one half of the interval will be taken recursively - this is not necessarily the correct interval.

Why is this not true? What is the pitfall? The table of values used for the graph were taken from the printed output of the function evaluation program which uses "printf" for display. The problem resides in the use of a fixed-point format specifier ("%f") whereas we should be using the full precision of the float by using exponential specifier ("%e"). The former truncates the result to a fixed number of decimal digits causing "chopping" of the significant digits out to 6 or 8 decimal places. So examining this output, we find that there is indeed only ONE zero - the bisection would, in fact, work fine!

Nevertheless, I left the paragraph here not only to show a possible pitfall, but also to show that potential of the type of error previously desribed is still there for higher powers. The use of higher precision floats just pushes the error farther out-of-reach.