roots.texi

开放gsl矩阵运算

@cindex root finding
@cindex zero finding
@cindex finding roots
@cindex finding zeros
@cindex roots
@cindex solving a non-linear equation
@cindex non-linear 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 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 go awry, 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 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 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 (* @var{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 * @var{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 (* @var{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 (* @var{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 (* @var{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 * @var{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