[Help-gsl] Root finding, false position method: excessive function evaluations

[Maxim Nikulin]

Hi,

I have noticed that the GSL implementation of the false position method 
for one-dimensional root finding may require twice more function 
evaluations than the bisection method for an almost linear function that 
should not be considered as a degenerate case.

Since the reference states:
> Its convergence is linear, but it is usually faster than bisection.

I think it worths a remark in the reference.

I suppose that it is a design flaw that user supplied tolerance is not 
available in the iteration functions. The implementation of the Brent 
method uses DBL_EPSILON, but in general user might request significantly 
lower precision. As a result it deteriorates programs performance.

The problem appears if at a certain step the edge of a bracketing 
interval almost coincides the root. Further linear interpolation gives 
the same point (either due to rounding error or withing the required 
precision). Function evaluation should be skipped for this point before 
bisection step.

The following program demonstrates the issue. The bisection method 
converges after ~50 function evaluations, the false position method 
requires ~100 function evaluations for the same precision. I would 
stress again that it is a well behaviored almost linear function.

Maxim Nikulin

/* Excessive function evaluations for an almost linear function
 * in the current implementation of the regula falsi method
 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

double fpol2 (double x, void *params)
{
   double *p = params;
   return (x*p[2] + p[1])*x + p[0];
}

struct f_count_wrap_par {
   gsl_function *F;
   int n;
};

double f_count_wrap (double x, void *params)
{
   struct f_count_wrap_par *f = params;
   ++f->n;
   return GSL_FN_EVAL(f->F, x);
}

int main (void)
{
   int i, max_iter = 50;
   double x_i, x_l = 0.5, x_h = 1.1;
   double tol_abs = 5*GSL_DBL_EPSILON;
   double tol_rel = 5*GSL_DBL_EPSILON;
   gsl_root_fsolver *s;
   double par[] = { -1, 1, -1e-10 }; /* f(x) = -1e-10*x*x + (x - 1) */
   gsl_function F = {fpol2, par}, Fwrap;
   struct f_count_wrap_par wrap_par = {&F, 0};
   int status;
   double r = -2*par[0]/(sqrt(par[1]*par[1]-4*par[0]*par[2])+par[1]);

   Fwrap.function = f_count_wrap;
   Fwrap.params = &wrap_par;

   s = gsl_root_fsolver_alloc (gsl_root_fsolver_falsepos);
   /* s = gsl_root_fsolver_alloc (gsl_root_fsolver_bisection); */
   gsl_root_fsolver_set (s, &Fwrap, x_l, x_h);
   printf ("# i   N_f   x_l-r      x_i-r      x_h-r      x_h-x_l\n");
   for (i = 0, status = GSL_CONTINUE ;
         status == GSL_CONTINUE && i < max_iter ; ++i) {
      status = gsl_root_fsolver_iterate (s);
      x_i = gsl_root_fsolver_root (s);
      x_l = gsl_root_fsolver_x_lower (s);
      x_h = gsl_root_fsolver_x_upper (s);
      if (status != GSL_SUCCESS)
         break;
      status = gsl_root_test_interval (x_l, x_h, tol_abs, tol_rel);
      printf ("%3d %5d % .3e % .3e % .3e % .3e\n", i, wrap_par.n,
            x_l - r, x_i - r, x_h - r, x_h - x_l);
   }
   if (status == GSL_SUCCESS) {
      printf ("# Converged\n");
      printf ("# f(x_l) = % .4e, f(x_h) = % .4e\n",
            GSL_FN_EVAL(&F, x_l), GSL_FN_EVAL(&F, x_h));
   } else if (i == max_iter)
      printf ("# Max iteration number exceeded\n");
   gsl_root_fsolver_free (s);
   exit (EXIT_SUCCESS);
}

Output: iteration, number of function evaluations, shifts from the root 
of the left interval edge, the current root estimation, and the right 
edge, and inally the bracketing interval length

# i   N_f   x_l-r      x_i-r      x_h-r      x_h-x_l
  0     4 -2.000e-01  5.000e-12  5.000e-12  2.000e-01
  1     6 -1.000e-01  0.000e+00  0.000e+00  1.000e-01
  2     8 -5.000e-02  0.000e+00  0.000e+00  5.000e-02
  3    10 -2.500e-02  0.000e+00  0.000e+00  2.500e-02
  4    12 -1.250e-02  0.000e+00  0.000e+00  1.250e-02
  5    14 -6.250e-03  0.000e+00  0.000e+00  6.250e-03
  6    16 -3.125e-03  0.000e+00  0.000e+00  3.125e-03
  7    18 -1.563e-03  0.000e+00  0.000e+00  1.563e-03
  8    20 -7.813e-04  0.000e+00  0.000e+00  7.813e-04
  9    22 -3.906e-04  0.000e+00  0.000e+00  3.906e-04
 10    24 -1.953e-04  0.000e+00  0.000e+00  1.953e-04
 11    26 -9.766e-05  0.000e+00  0.000e+00  9.766e-05
 12    28 -4.883e-05  0.000e+00  0.000e+00  4.883e-05
 13    30 -2.441e-05  0.000e+00  0.000e+00  2.441e-05
 14    32 -1.221e-05  0.000e+00  0.000e+00  1.221e-05
 15    34 -6.104e-06  0.000e+00  0.000e+00  6.104e-06
 16    36 -3.052e-06  0.000e+00  0.000e+00  3.052e-06
 17    38 -1.526e-06  0.000e+00  0.000e+00  1.526e-06
 18    40 -7.629e-07  0.000e+00  0.000e+00  7.629e-07
 19    42 -3.815e-07  0.000e+00  0.000e+00  3.815e-07
 20    44 -1.907e-07  0.000e+00  0.000e+00  1.907e-07
 21    46 -9.537e-08  0.000e+00  0.000e+00  9.537e-08
 22    48 -4.768e-08  0.000e+00  0.000e+00  4.768e-08
 23    50 -2.384e-08  0.000e+00  0.000e+00  2.384e-08
 24    52 -1.192e-08  0.000e+00  0.000e+00  1.192e-08
 25    54 -5.960e-09  0.000e+00  0.000e+00  5.960e-09
 26    56 -2.980e-09  0.000e+00  0.000e+00  2.980e-09
 27    58 -1.490e-09  0.000e+00  0.000e+00  1.490e-09
 28    60 -7.451e-10  0.000e+00  0.000e+00  7.451e-10
 29    62 -3.725e-10  0.000e+00  0.000e+00  3.725e-10
 30    64 -1.863e-10  0.000e+00  0.000e+00  1.863e-10
 31    66 -9.313e-11  0.000e+00  0.000e+00  9.313e-11
 32    68 -4.657e-11  0.000e+00  0.000e+00  4.657e-11
 33    70 -2.328e-11  0.000e+00  0.000e+00  2.328e-11
 34    72 -1.164e-11  0.000e+00  0.000e+00  1.164e-11
 35    74 -5.821e-12  0.000e+00  0.000e+00  5.821e-12
 36    76 -2.910e-12  0.000e+00  0.000e+00  2.910e-12
 37    78 -1.455e-12  0.000e+00  0.000e+00  1.455e-12
 38    80 -7.276e-13  0.000e+00  0.000e+00  7.276e-13
 39    82 -3.637e-13  0.000e+00  0.000e+00  3.637e-13
 40    84 -1.819e-13  0.000e+00  0.000e+00  1.819e-13
 41    86 -9.104e-14  0.000e+00  0.000e+00  9.104e-14
 42    88 -4.552e-14  0.000e+00  0.000e+00  4.552e-14
 43    90 -2.265e-14  0.000e+00  0.000e+00  2.265e-14
 44    92 -1.132e-14  0.000e+00  0.000e+00  1.132e-14
 45    94 -5.773e-15  0.000e+00  0.000e+00  5.773e-15
 46    96 -2.887e-15  0.000e+00  0.000e+00  2.887e-15
 47    98 -1.332e-15  0.000e+00  0.000e+00  1.332e-15
# Converged
# f(x_l) = -1.3240e-15, f(x_h) =  8.2399e-18