@cindex root finding
@cindex zero finding
@cindex finding roots
@cindex finding zeros
@cindex roots
@cindex solving a nonlinear equation
@cindex nonlinear equation, solutions of

This chapter describes routines for finding roots of arbitrary
one-dimensional functions.  The library provides low level components
for a variety of iterative solvers and convergence tests.  These can be
combined by the user to achieve the desired solution, with full access
to the intermediate steps of the iteration.  Each class of methods uses
the same framework, so that you can switch between solvers at runtime
without needing to recompile your program.  Each instance of a solver
keeps track of its own state, allowing the solvers to be used in
multi-threaded programs.

The header file @file{gsl_roots.h} contains prototypes for the root
finding functions and related declarations.

@menu
* Root Finding Overview::       
* Root Finding Caveats::        
* Initializing the Solver::     
* Providing the function to solve::  
* Search Bounds and Guesses::   
* Root Finding Iteration::      
* Search Stopping Parameters::  
* Root Bracketing Algorithms::  
* Root Finding Algorithms using Derivatives::  
* Root Finding Examples::       
* Root Finding References and Further Reading::  
@end menu

@node Root Finding Overview
@section Overview
@cindex root finding, overview

One-dimensional root finding algorithms can be divided into two classes,
@dfn{root bracketing} and @dfn{root polishing}.  Algorithms which proceed
by bracketing a root are guaranteed to converge.  Bracketing algorithms
begin with a bounded region known to contain a root.  The size of this
bounded region is reduced, iteratively, until it encloses the root to a
desired tolerance.  This provides a rigorous error estimate for the
location of the root.

The technique of @dfn{root polishing} attempts to improve an initial
guess to the root.  These algorithms converge only if started ``close
enough'' to a root, and sacrifice a rigorous error bound for speed.  By
approximating the behavior of a function in the vicinity of a root they
attempt to find a higher order improvement of an initial guess.  When the
behavior of the function is compatible with the algorithm and a good
initial guess is available a polishing algorithm can provide rapid
convergence.

In GSL both types of algorithm are available in similar frameworks.  The
user provides a high-level driver for the algorithms, and the library
provides the individual functions necessary for each of the steps.
There are three main phases of the iteration.  The steps are,

@itemize @bullet
@item
initialize solver state, @var{s}, for algorithm @var{T}

@item
update @var{s} using the iteration @var{T}

@item
test @var{s} for convergence, and repeat iteration if necessary
@end itemize

@noindent
The state for bracketing solvers is held in a @code{gsl_root_fsolver}
struct.  The updating procedure uses only function evaluations (not
derivatives).  The state for root polishing solvers is held in a
@code{gsl_root_fdfsolver} struct.  The updates require both the function
and its derivative (hence the name @code{fdf}) to be supplied by the
user.

@node Root Finding Caveats
@section Caveats
@cindex root finding, caveats

Note that root finding functions can only search for one root at a time.
When there are several roots in the search area, the first root to be
found will be returned; however it is difficult to predict which of the
roots this will be. @emph{In most cases, no error will be reported if
you try to find a root in an area where there is more than one.}

Care must be taken when a function may have a multiple root (such as 
@c{$f(x) = (x-x_0)^2$}
@math{f(x) = (x-x_0)^2} or 
@c{$f(x) = (x-x_0)^3$}
@math{f(x) = (x-x_0)^3}).  
It is not possible to use root-bracketing algorithms on
even-multiplicity roots.  For these algorithms the initial interval must
contain a zero-crossing, where the function is negative at one end of
the interval and positive at the other end.  Roots with even-multiplicity
do not cross zero, but only touch it instantaneously.  Algorithms based
on root bracketing will still work for odd-multiplicity roots
(e.g. cubic, quintic, @dots{}). 
Root polishing algorithms generally work with higher multiplicity roots,
but at a reduced rate of convergence.  In these cases the @dfn{Steffenson
algorithm} can be used to accelerate the convergence of multiple roots.

While it is not absolutely required that @math{f} have a root within the
search region, numerical root finding functions should not be used
haphazardly to check for the @emph{existence} of roots.  There are better
ways to do this.  Because it is easy to create situations where numerical
root finders can fail, it is a bad idea to throw a root finder at a
function you do not know much about.  In general it is best to examine
the function visually by plotting before searching for a root.

@node Initializing the Solver
@section Initializing the Solver

@deftypefun {gsl_root_fsolver *} gsl_root_fsolver_alloc (const gsl_root_fsolver_type * @var{T})
This function returns a pointer to a newly allocated instance of a
solver of type @var{T}.  For example, the following code creates an
instance of a bisection solver,

@example
const gsl_root_fsolver_type * T 
  = gsl_root_fsolver_bisection;
gsl_root_fsolver * s 
  = gsl_root_fsolver_alloc (T);
@end example

If there is insufficient memory to create the solver then the function
returns a null pointer and the error handler is invoked with an error
code of @code{GSL_ENOMEM}.
@end deftypefun

@deftypefun {gsl_root_fdfsolver *} gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * @var{T})
This function returns a pointer to a newly allocated instance of a
derivative-based solver of type @var{T}.  For example, the following
code creates an instance of a Newton-Raphson solver,

@example
const gsl_root_fdfsolver_type * T 
  = gsl_root_fdfsolver_newton;
gsl_root_fdfsolver * s 
  = gsl_root_fdfsolver_alloc (T);
@end example

If there is insufficient memory to create the solver then the function
returns a null pointer and the error handler is invoked with an error
code of @code{GSL_ENOMEM}.
@end deftypefun


@deftypefun int gsl_root_fsolver_set (gsl_root_fsolver * @var{s}, gsl_function * @var{f}, double @var{x_lower}, double @var{x_upper})
This function initializes, or reinitializes, an existing solver @var{s}
to use the function @var{f} and the initial search interval
[@var{x_lower}, @var{x_upper}].
@end deftypefun

@deftypefun int gsl_root_fdfsolver_set (gsl_root_fdfsolver * @var{s}, gsl_function_fdf * @var{fdf}, double @var{root})
This function initializes, or reinitializes, an existing solver @var{s}
to use the function and derivative @var{fdf} and the initial guess
@var{root}.
@end deftypefun

@deftypefun void gsl_root_fsolver_free (gsl_root_fsolver * @var{s})
@deftypefunx void gsl_root_fdfsolver_free (gsl_root_fdfsolver * @var{s})
These functions free all the memory associated with the solver @var{s}.
@end deftypefun

@deftypefun {const char *} gsl_root_fsolver_name (const gsl_root_fsolver * @var{s})
@deftypefunx {const char *} gsl_root_fdfsolver_name (const gsl_root_fdfsolver * @var{s})
These functions return a pointer to the name of the solver.  For example,

@example
printf ("s is a '%s' solver\n",
        gsl_root_fsolver_name (s));
@end example

@noindent
would print something like @code{s is a 'bisection' solver}.
@end deftypefun

@node Providing the function to solve
@section Providing the function to solve
@cindex root finding, providing a function to solve

You must provide a continuous function of one variable for the root
finders to operate on, and, sometimes, its first derivative.  In order
to allow for general parameters the functions are defined by the
following data types:

@deftp {Data Type} gsl_function 
This data type defines a general function with parameters. 

@table @code
@item double (* function) (double @var{x}, void * @var{params})
this function should return the value
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} for argument @var{x} and parameters @var{params}

@item void * params
a pointer to the parameters of the function
@end table
@end deftp

Here is an example for the general quadratic function,
@tex
\beforedisplay
$$
f(x) = a x^2 + b x + c
$$
\afterdisplay
@end tex
@ifinfo

@example
f(x) = a x^2 + b x + c
@end example

@end ifinfo
@noindent
with @math{a = 3}, @math{b = 2}, @math{c = 1}.  The following code
defines a @code{gsl_function} @code{F} which you could pass to a root
finder:

@example
struct my_f_params @{ double a; double b; double c; @};

double
my_f (double x, void * p) @{
   struct my_f_params * params 
     = (struct my_f_params *)p;
   double a = (params->a);
   double b = (params->b);
   double c = (params->c);

   return  (a * x + b) * x + c;
@}

gsl_function F;
struct my_f_params params = @{ 3.0, 2.0, 1.0 @};

F.function = &my_f;
F.params = &params;
@end example

@noindent
The function @math{f(x)} can be evaluated using the following macro,

@example
#define GSL_FN_EVAL(F,x) 
    (*((F)->function))(x,(F)->params)
@end example

@deftp {Data Type} gsl_function_fdf
This data type defines a general function with parameters and its first
derivative.

@table @code
@item double (* f) (double @var{x}, void * @var{params})
this function should return the value of
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} for argument @var{x} and parameters @var{params}

@item double (* df) (double @var{x}, void * @var{params})
this function should return the value of the derivative of @var{f} with
respect to @var{x},
@c{$f'(x,\hbox{\it params})$}
@math{f'(x,params)}, for argument @var{x} and parameters @var{params}

@item void (* fdf) (double @var{x}, void * @var{params}, double * @var{f}, double * @var{d}f)
this function should set the values of the function @var{f} to 
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)}
and its derivative @var{df} to
@c{$f'(x,\hbox{\it params})$}
@math{f'(x,params)} 
for argument @var{x} and parameters @var{params}.  This function
provides an optimization of the separate functions for @math{f(x)} and
@math{f'(x)}---it is always faster to compute the function and its
derivative at the same time.

@item void * params
a pointer to the parameters of the function
@end table
@end deftp

Here is an example where 
@c{$f(x) = \exp(2x)$}
@math{f(x) = 2\exp(2x)}:

@example
double
my_f (double x, void * params)
@{
   return exp (2 * x);
@}

double
my_df (double x, void * params)
@{
   return 2 * exp (2 * x);
@}

void
my_fdf (double x, void * params, 
        double * f, double * df)
@{
   double t = exp (2 * x);

   *f = t;
   *df = 2 * t;   /* uses existing value */
@}

gsl_function_fdf FDF;

FDF.f = &my_f;
FDF.df = &my_df;
FDF.fdf = &my_fdf;
FDF.params = 0;
@end example

@noindent
The function @math{f(x)} can be evaluated using the following macro,

@example
#define GSL_FN_FDF_EVAL_F(FDF,x) 
     (*((FDF)->f))(x,(FDF)->params)
@end example

@noindent
The derivative @math{f'(x)} can be evaluated using the following macro,

@example
#define GSL_FN_FDF_EVAL_DF(FDF,x) 
     (*((FDF)->df))(x,(FDF)->params)
@end example

@noindent
and both the function @math{y = f(x)} and its derivative @math{dy = f'(x)}
can be evaluated at the same time using the following macro,

@example
#define GSL_FN_FDF_EVAL_F_DF(FDF,x,y,dy) 
     (*((FDF)->fdf))(x,(FDF)->params,(y),(dy))
@end example

@noindent
The macro stores @math{f(x)} in its @var{y} argument and @math{f'(x)} in
its @var{dy} argument---both of these should be pointers to
@code{double}.

@node Search Bounds and Guesses
@section Search Bounds and Guesses
@cindex root finding, search bounds
@cindex root finding, initial guess

You provide either search bounds or an initial guess; this section
explains how search bounds and guesses work and how function arguments
control them.

A guess is simply an @math{x} value which is iterated until it is within
the desired precision of a root.  It takes the form of a @code{double}.

Search bounds are the endpoints of a interval which is iterated until
the length of the interval is smaller than the requested precision.  The
interval is defined by two values, the lower limit and the upper limit.
Whether the endpoints are intended to be included in the interval or not
depends on the context in which the interval is used.

@node Root Finding Iteration
@section Iteration

The following functions drive the iteration of each algorithm.  Each
function performs one iteration to update the state of any solver of the
corresponding type.  The same functions work for all solvers so that
different methods can be substituted at runtime without modifications to
the code.

@deftypefun int gsl_root_fsolver_iterate (gsl_root_fsolver * @var{s})
@deftypefunx int gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * @var{s})
These functions perform a single iteration of the solver @var{s}.  If the
iteration encounters an unexpected problem then an error code will be
returned,

@table @code
@item GSL_EBADFUNC
the iteration encountered a singular point where the function or its
derivative evaluated to @code{Inf} or @code{NaN}.

@item GSL_EZERODIV
the derivative of the function vanished at the iteration point,
preventing the algorithm from continuing without a division by zero.
@end table
@end deftypefun

The solver maintains a current best estimate of the root at all
times.  The bracketing solvers also keep track of the current best
interval bounding the root.  This information can be accessed with the
following auxiliary functions,

@deftypefun double gsl_root_fsolver_root (const gsl_root_fsolver * @var{s})
@deftypefunx double gsl_root_fdfsolver_root (const gsl_root_fdfsolver * @var{s})
These functions return the current estimate of the root for the solver @var{s}.
@end deftypefun

@deftypefun double gsl_root_fsolver_x_lower (const gsl_root_fsolver * @var{s})
@deftypefunx double gsl_root_fsolver_x_upper (const gsl_root_fsolver * @var{s})
These functions return the current bracketing interval for the solver @var{s}.
@end deftypefun

@node Search Stopping Parameters
@section Search Stopping Parameters
@cindex root finding, stopping parameters

A root finding procedure should stop when one of the following conditions is
true:

@itemize @bullet
@item
A root has been found to within the user-specified precision.

@item
A user-specified maximum number of iterations has been reached.

@item
An error has occurred.
@end itemize

@noindent
The handling of these conditions is under user control.  The functions
below allow the user to test the precision of the current result in
several standard ways.

@deftypefun int gsl_root_test_interval (double @var{x_lower}, double @var{x_upper}, double @var{epsabs}, double @var{epsrel})
This function tests for the convergence of the interval [@var{x_lower},
@var{x_upper}] with absolute error @var{epsabs} and relative error
@var{epsrel}.  The test returns @code{GSL_SUCCESS} if the following
condition is achieved,
@tex
\beforedisplay
$$
|a - b| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, \min(|a|,|b|)
$$
\afterdisplay
@end tex
@ifinfo

@example
|a - b| < epsabs + epsrel min(|a|,|b|) 
@end example

@end ifinfo
@noindent
when the interval @math{x = [a,b]} does not include the origin.  If the
interval includes the origin then @math{\min(|a|,|b|)} is replaced by
zero (which is the minimum value of @math{|x|} over the interval).  This
ensures that the relative error is accurately estimated for roots close
to the origin.

This condition on the interval also implies that any estimate of the
root @math{r} in the interval satisfies the same condition with respect
to the true root @math{r^*},
@tex
\beforedisplay
$$
|r - r^*| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, r^*
$$
\afterdisplay
@end tex
@ifinfo

@example
|r - r^*| < epsabs + epsrel r^*
@end example

@end ifinfo
@noindent
assuming that the true root @math{r^*} is contained within the interval.
@end deftypefun

@deftypefun int gsl_root_test_delta (double @var{x1}, double @var{x0}, double @var{epsabs}, double @var{epsrel})

This function tests for the convergence of the sequence @dots{}, @var{x0},
@var{x1} with absolute error @var{epsabs} and relative error
@var{epsrel}.  The test returns @code{GSL_SUCCESS} if the following
condition is achieved,
@tex
\beforedisplay
$$
|x_1 - x_0| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, |x_1|
$$
\afterdisplay
@end tex
@ifinfo

@example
|x_1 - x_0| < epsabs + epsrel |x_1|
@end example

@end ifinfo
@noindent
and returns @code{GSL_CONTINUE} otherwise.
@end deftypefun


@deftypefun int gsl_root_test_residual (double @var{f}, double @var{epsabs})
This function tests the residual value @var{f} against the absolute
error bound @var{epsabs}.  The test returns @code{GSL_SUCCESS} if the
following condition is achieved,
@tex
\beforedisplay
$$
|f| < \hbox{\it epsabs}
$$
\afterdisplay
@end tex
@ifinfo

@example
|f| < epsabs
@end example

@end ifinfo
@noindent
and returns @code{GSL_CONTINUE} otherwise.  This criterion is suitable
for situations where the precise location of the root, @math{x}, is
unimportant provided a value can be found where the residual,
@math{|f(x)|}, is small enough.
@end deftypefun

@comment ============================================================

@node Root Bracketing Algorithms
@section Root Bracketing Algorithms

The root bracketing algorithms described in this section require an
initial interval which is guaranteed to contain a root---if @math{a}
and @math{b} are the endpoints of the interval then @math{f(a)} must
differ in sign from @math{f(b)}.  This ensures that the function crosses
zero at least once in the interval.  If a valid initial interval is used
then these algorithm cannot fail, provided the function is well-behaved.

Note that a bracketing algorithm cannot find roots of even degree, since
these do not cross the @math{x}-axis.

@deffn {Solver} gsl_root_fsolver_bisection

@cindex bisection algorithm for finding roots
@cindex root finding, bisection algorithm

The @dfn{bisection algorithm} is the simplest method of bracketing the
roots of a function.   It is the slowest algorithm provided by
the library, with linear convergence.

On each iteration, the interval is bisected and the value of the
function at the midpoint is calculated.  The sign of this value is used
to determine which half of the interval does not contain a root.  That
half is discarded to give a new, smaller interval containing the
root.  This procedure can be continued indefinitely until the interval is
sufficiently small.

At any time the current estimate of the root is taken as the midpoint of
the interval.

@comment eps file "roots-bisection.eps"
@comment @iftex
@comment @sp 1
@comment @center @image{roots-bisection,3.4in}

@comment @quotation
@comment Four iterations of bisection, where @math{a_n} is @math{n}th position of
@comment the beginning of the interval and @math{b_n} is the @math{n}th position
@comment of the end.  The midpoint of each interval is also indicated.
@comment @end quotation
@comment @end iftex
@end deffn

@comment ============================================================

@deffn {Solver} gsl_root_fsolver_falsepos
@cindex false position algorithm for finding roots
@cindex root finding, false position algorithm

The @dfn{false position algorithm} is a method of finding roots based on
linear interpolation.  Its convergence is linear, but it is usually
faster than bisection.

On each iteration a line is drawn between the endpoints @math{(a,f(a))}
and @math{(b,f(b))} and the point where this line crosses the
@math{x}-axis taken as a ``midpoint''.  The value of the function at
this point is calculated and its sign is used to determine which side of
the interval does not contain a root.  That side is discarded to give a
new, smaller interval containing the root.  This procedure can be
continued indefinitely until the interval is sufficiently small.

The best estimate of the root is taken from the linear interpolation of
the interval on the current iteration.

@comment eps file "roots-false-position.eps"
@comment @iftex
@comment @image{roots-false-position,4in}
@comment @quotation
@comment Several iterations of false position, where @math{a_n} is @math{n}th
@comment position of the beginning of the interval and @math{b_n} is the
@comment @math{n}th position of the end.
@comment @end quotation
@comment @end iftex
@end deffn

@comment ============================================================

@deffn {Solver} gsl_root_fsolver_brent
@cindex Brent's method for finding roots
@cindex root finding, Brent's method

The @dfn{Brent-Dekker method} (referred to here as @dfn{Brent's method})
combines an interpolation strategy with the bisection algorithm.  This
produces a fast algorithm which is still robust.

On each iteration Brent's method approximates the function using an
interpolating curve.  On the first iteration this is a linear
interpolation of the two endpoints.  For subsequent iterations the
algorithm uses an inverse quadratic fit to the last three points, for
higher accuracy.  The intercept of the interpolating curve with the
@math{x}-axis is taken as a guess for the root.  If it lies within the
bounds of the current interval then the interpolating point is accepted,
and used to generate a smaller interval.  If the interpolating point is
not accepted then the algorithm falls back to an ordinary bisection
step.

The best estimate of the root is taken from the most recent
interpolation or bisection.
@end deffn

@comment ============================================================

@node Root Finding Algorithms using Derivatives
@section Root Finding Algorithms using Derivatives

The root polishing algorithms described in this section require an
initial guess for the location of the root.  There is no absolute
guarantee of convergence---the function must be suitable for this
technique and the initial guess must be sufficiently close to the root
for it to work.  When these conditions are satisfied then convergence is
quadratic.

These algorithms make use of both the function and its derivative. 

@deffn {Derivative Solver} gsl_root_fdfsolver_newton
@cindex Newton's method for finding roots
@cindex root finding, Newton's method

Newton's Method is the standard root-polishing algorithm.  The algorithm
begins with an initial guess for the location of the root.  On each
iteration, a line tangent to the function @math{f} is drawn at that
position.  The point where this line crosses the @math{x}-axis becomes
the new guess.  The iteration is defined by the following sequence,
@tex
\beforedisplay
$$
x_{i+1} = x_i - {f(x_i) \over f'(x_i)}
$$
\afterdisplay
@end tex
@ifinfo

@example
x_@{i+1@} = x_i - f(x_i)/f'(x_i)
@end example

@end ifinfo
@noindent
Newton's method converges quadratically for single roots, and linearly
for multiple roots.

@comment eps file "roots-newtons-method.eps"
@comment @iftex
@comment @sp 1
@comment @center @image{roots-newtons-method,3.4in}

@comment @quotation
@comment Several iterations of Newton's Method, where @math{g_n} is the
@comment @math{n}th guess.
@comment @end quotation
@comment @end iftex
@end deffn

@comment ============================================================

@deffn {Derivative Solver} gsl_root_fdfsolver_secant
@cindex secant method for finding roots
@cindex root finding, secant method

The @dfn{secant method} is a simplified version of Newton's method which does
not require the computation of the derivative on every step.

On its first iteration the algorithm begins with Newton's method, using
the derivative to compute a first step,
@tex
\beforedisplay
$$
x_1 = x_0 - {f(x_0) \over f'(x_0)}
$$
\afterdisplay
@end tex
@ifinfo

@example
x_1 = x_0 - f(x_0)/f'(x_0)
@end example

@end ifinfo
@noindent
Subsequent iterations avoid the evaluation of the derivative by
replacing it with a numerical estimate, the slope of the line through
the previous two points,
@tex
\beforedisplay
$$
x_{i+1} = x_i - {f(x_i) \over f'_{est}}
 ~\hbox{where}~
 f'_{est} =  {f(x_{i}) - f(x_{i-1}) \over x_i - x_{i-1}}
$$
\afterdisplay
@end tex
@ifinfo

@example
x_@{i+1@} = x_i f(x_i) / f'_@{est@} where
 f'_@{est@} = (f(x_i) - f(x_@{i-1@})/(x_i - x_@{i-1@})
@end example

@end ifinfo
@noindent
When the derivative does not change significantly in the vicinity of the
root the secant method gives a useful saving.  Asymptotically the secant
method is faster than Newton's method whenever the cost of evaluating
the derivative is more than 0.44 times the cost of evaluating the
function itself.  As with all methods of computing a numerical
derivative the estimate can suffer from cancellation errors if the
separation of the points becomes too small.

On single roots, the method has a convergence of order @math{(1 + \sqrt
5)/2} (approximately @math{1.62}).  It converges linearly for multiple
roots.  

@comment eps file "roots-secant-method.eps"
@comment @iftex
@comment @tex
@comment \input epsf
@comment \medskip
@comment \centerline{\epsfxsize=5in\epsfbox{roots-secant-method.eps}}
@comment @end tex
@comment @quotation
@comment Several iterations of Secant Method, where @math{g_n} is the @math{n}th
@comment guess.
@comment @end quotation
@comment @end iftex
@end deffn

@comment ============================================================

@deffn {Derivative Solver} gsl_root_fdfsolver_steffenson
@cindex Steffenson's method for finding roots
@cindex root finding, Steffenson's method

The @dfn{Steffenson Method} provides the fastest convergence of all the
routines.  It combines the basic Newton algorithm with an Aitken
``delta-squared'' acceleration.  If the Newton iterates are @math{x_i}
then the acceleration procedure generates a new sequence @math{R_i},
@tex
\beforedisplay
$$
R_i = x_i - {(x_{i+1} - x_i)^2 \over (x_{i+2} - 2 x_{i+1} + x_i)}
$$
\afterdisplay
@end tex
@ifinfo

@example
R_i = x_i - (x_@{i+1@} - x_i)^2 / (x_@{i+2@} - 2 x_@{i+1@} + x_@{i@})
@end example

@end ifinfo
@noindent
which converges faster than the original sequence under reasonable
conditions.  The new sequence requires three terms before it can produce
its first value so the method returns accelerated values on the second
and subsequent iterations.  On the first iteration it returns the
ordinary Newton estimate.  The Newton iterate is also returned if the
denominator of the acceleration term ever becomes zero.

As with all acceleration procedures this method can become unstable if
the function is not well-behaved. 
@end deffn

@node Root Finding Examples
@section 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 (@file{demo_fn.h}) to
define the function parameters,

@example
@verbatiminclude examples/demo_fn.h
@end example

@noindent
We place the function definitions in a separate file (@file{demo_fn.c}),

@example
@verbatiminclude examples/demo_fn.c
@end example

@noindent
The first program uses the function solver @code{gsl_root_fsolver_brent}
for Brent's method and the general quadratic defined above to solve the
following equation,
@tex
\beforedisplay
$$
x^2 - 5 = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
x^2 - 5 = 0
@end example

@end ifinfo
@noindent
with solution @math{x = \sqrt 5 = 2.236068...}

@example
@verbatiminclude examples/roots.c
@end example

@noindent
Here are the results of the iterations,

@smallexample
$ ./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
@end smallexample

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

@smallexample
$ ./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
@end smallexample

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

@example
@verbatiminclude examples/rootnewt.c
@end example

@noindent
Here are the results for Newton's method,

@example
$ ./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
@end example

@noindent
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 @code{gsl_root_fdfsolver_newton} to
@code{gsl_root_fdfsolver_secant} or
@code{gsl_root_fdfsolver_steffenson}.

@node Root Finding References and Further Reading
@section References and Further Reading

For information on the Brent-Dekker algorithm see the following two
papers,

@itemize @asis
@item
R. P. Brent, ``An algorithm with guaranteed convergence for finding a
zero of a function'', @cite{Computer Journal}, 14 (1971) 422--425

@item
J. C. P. Bus and T. J. Dekker, ``Two Efficient Algorithms with Guaranteed
Convergence for Finding a Zero of a Function'', @cite{ACM Transactions of
Mathematical Software}, Vol.@: 1 No.@: 4 (1975) 330--345
@end itemize