Next: , Previous: Root Finding Algorithms using Derivatives, Up: One dimensional Root-Finding


32.10 Examples

For any root finding algorithm we need to prepare the function to be solved. For this example we will use the general quadratic equation described earlier. We first need a header file (demo_fn.h) to define the function parameters,

     struct quadratic_params
       {
         double a, b, c;
       };
     
     double quadratic (double x, void *params);
     double quadratic_deriv (double x, void *params);
     void quadratic_fdf (double x, void *params, 
                         double *y, double *dy);

We place the function definitions in a separate file (demo_fn.c),

     double
     quadratic (double x, void *params)
     {
       struct quadratic_params *p 
         = (struct quadratic_params *) params;
     
       double a = p->a;
       double b = p->b;
       double c = p->c;
     
       return (a * x + b) * x + c;
     }
     
     double
     quadratic_deriv (double x, void *params)
     {
       struct quadratic_params *p 
         = (struct quadratic_params *) params;
     
       double a = p->a;
       double b = p->b;
       double c = p->c;
     
       return 2.0 * a * x + b;
     }
     
     void
     quadratic_fdf (double x, void *params, 
                    double *y, double *dy)
     {
       struct quadratic_params *p 
         = (struct quadratic_params *) params;
     
       double a = p->a;
       double b = p->b;
       double c = p->c;
     
       *y = (a * x + b) * x + c;
       *dy = 2.0 * a * x + b;
     }

The first program uses the function solver gsl_root_fsolver_brent for Brent's method and the general quadratic defined above to solve the following equation,

     x^2 - 5 = 0

with solution x = \sqrt 5 = 2.236068...

     #include <stdio.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_roots.h>
     
     #include "demo_fn.h"
     #include "demo_fn.c"
     
     int
     main (void)
     {
       int status;
       int iter = 0, max_iter = 100;
       const gsl_root_fsolver_type *T;
       gsl_root_fsolver *s;
       double r = 0, r_expected = sqrt (5.0);
       double x_lo = 0.0, x_hi = 5.0;
       gsl_function F;
       struct quadratic_params params = {1.0, 0.0, -5.0};
     
       F.function = &quadratic;
       F.params = &params;
     
       T = gsl_root_fsolver_brent;
       s = gsl_root_fsolver_alloc (T);
       gsl_root_fsolver_set (s, &F, x_lo, x_hi);
     
       printf ("using %s method\n", 
               gsl_root_fsolver_name (s));
     
       printf ("%5s [%9s, %9s] %9s %10s %9s\n",
               "iter", "lower", "upper", "root", 
               "err", "err(est)");
     
       do
         {
           iter++;
           status = gsl_root_fsolver_iterate (s);
           r = gsl_root_fsolver_root (s);
           x_lo = gsl_root_fsolver_x_lower (s);
           x_hi = gsl_root_fsolver_x_upper (s);
           status = gsl_root_test_interval (x_lo, x_hi,
                                            0, 0.001);
     
           if (status == GSL_SUCCESS)
             printf ("Converged:\n");
     
           printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
                   iter, x_lo, x_hi,
                   r, r - r_expected, 
                   x_hi - x_lo);
         }
       while (status == GSL_CONTINUE && iter < max_iter);
       return status;
     }

Here are the results of the iterations,

     $ ./a.out
     using brent method
      iter [    lower,     upper]      root        err  err(est)
         1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
         2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
         3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
         4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
         5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
     Converged:
         6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666

If the program is modified to use the bisection solver instead of Brent's method, by changing gsl_root_fsolver_brent to gsl_root_fsolver_bisection the slower convergence of the Bisection method can be observed,

     $ ./a.out
     using bisection method
      iter [    lower,     upper]      root        err  err(est)
         1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000
         2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000
         3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000
         4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000
         5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500
         6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250
         7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625
         8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312
         9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656
        10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828
        11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414
     Converged:
        12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207

The next program solves the same function using a derivative solver instead.

     #include <stdio.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_roots.h>
     
     #include "demo_fn.h"
     #include "demo_fn.c"
     
     int
     main (void)
     {
       int status;
       int iter = 0, max_iter = 100;
       const gsl_root_fdfsolver_type *T;
       gsl_root_fdfsolver *s;
       double x0, x = 5.0, r_expected = sqrt (5.0);
       gsl_function_fdf FDF;
       struct quadratic_params params = {1.0, 0.0, -5.0};
     
       FDF.f = &quadratic;
       FDF.df = &quadratic_deriv;
       FDF.fdf = &quadratic_fdf;
       FDF.params = &params;
     
       T = gsl_root_fdfsolver_newton;
       s = gsl_root_fdfsolver_alloc (T);
       gsl_root_fdfsolver_set (s, &FDF, x);
     
       printf ("using %s method\n", 
               gsl_root_fdfsolver_name (s));
     
       printf ("%-5s %10s %10s %10s\n",
               "iter", "root", "err", "err(est)");
       do
         {
           iter++;
           status = gsl_root_fdfsolver_iterate (s);
           x0 = x;
           x = gsl_root_fdfsolver_root (s);
           status = gsl_root_test_delta (x, x0, 0, 1e-3);
     
           if (status == GSL_SUCCESS)
             printf ("Converged:\n");
     
           printf ("%5d %10.7f %+10.7f %10.7f\n",
                   iter, x, x - r_expected, x - x0);
         }
       while (status == GSL_CONTINUE && iter < max_iter);
       return status;
     }

Here are the results for Newton's method,

     $ ./a.out
     using newton method
     iter        root        err   err(est)
         1  3.0000000 +0.7639320 -2.0000000
         2  2.3333333 +0.0972654 -0.6666667
         3  2.2380952 +0.0020273 -0.0952381
     Converged:
         4  2.2360689 +0.0000009 -0.0020263

Note that the error can be estimated more accurately by taking the difference between the current iterate and next iterate rather than the previous iterate. The other derivative solvers can be investigated by changing gsl_root_fdfsolver_newton to gsl_root_fdfsolver_secant or gsl_root_fdfsolver_steffenson.