multiroots.texi

开放gsl矩阵运算

@cindex solving nonlinear systems of equations
@cindex nonlinear systems of equations, solution of
@cindex systems of equations, nonlinear

This chapter describes functions for multidimensional root-finding
(solving nonlinear systems with @math{n} equations in @math{n}
unknowns).  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 solvers are based on the original Fortran
library @sc{minpack}.

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

@menu
* Overview of Multidimensional Root Finding::  
* Initializing the Multidimensional Solver::  
* Providing the multidimensional system of equations to solve::  
* Iteration of the multidimensional solver::  
* Search Stopping Parameters for the multidimensional solver::  
* Algorithms using Derivatives::  
* Algorithms without Derivatives::  
* Example programs for Multidimensional Root finding::  
* References and Further Reading for Multidimensional Root Finding::  
@end menu

@node Overview of Multidimensional Root Finding
@section Overview
@cindex multidimensional root finding, overview

The problem of multidimensional root finding requires the simultaneous
solution of @math{n} equations, @math{f_i}, in @math{n} variables,
@math{x_i},

@tex
\beforedisplay
$$
f_i (x_1, \dots, x_n) = 0 \qquad\hbox{for}~i = 1 \dots n.
$$
\afterdisplay
@end tex
@ifinfo
@example
f_i (x_1, ..., x_n) = 0    for i = 1 ... n.
@end example
@end ifinfo
@noindent
In general there are no bracketing methods available for @math{n}
dimensional systems, and no way of knowing whether any solutions
exist.  All algorithms proceed from an initial guess using a variant of
the Newton iteration,

@tex
\beforedisplay
$$
x \to x' = x - J^{-1} f(x)
$$
\afterdisplay
@end tex
@ifinfo
@example
x -> x' = x - J^@{-1@} f(x)
@end example
@end ifinfo
@noindent
where @math{x}, @math{f} are vector quantities and @math{J} is the
Jacobian matrix @c{$J_{ij} = \partial f_i / \partial x_j$} 
@math{J_@{ij@} = d f_i / d x_j}.
Additional strategies can be used to enlarge the region of
convergence.  These include requiring a decrease in the norm @math{|f|} on
each step proposed by Newton's method, or taking steepest-descent steps in
the direction of the negative gradient of @math{|f|}.

Several root-finding algorithms are available within a single framework.
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 evaluation of the Jacobian matrix can be problematic, either because
programming the derivatives is intractable or because computation of the
@math{n^2} terms of the matrix becomes too expensive.  For these reasons
the algorithms provided by the library are divided into two classes according
to whether the derivatives are available or not.

The state for solvers with an analytic Jacobian matrix is held in a
@code{gsl_multiroot_fdfsolver} struct.  The updating procedure requires
both the function and its derivatives to be supplied by the user.

The state for solvers which do not use an analytic Jacobian matrix is
held in a @code{gsl_multiroot_fsolver} struct.  The updating procedure
uses only function evaluations (not derivatives).  The algorithms
estimate the matrix @math{J} or @c{$J^{-1}$} 
@math{J^@{-1@}} by approximate methods.

@node Initializing the Multidimensional Solver
@section Initializing the Solver

The following functions initialize a multidimensional solver, either
with or without derivatives.  The solver itself depends only on the
dimension of the problem and the algorithm and can be reused for
different problems.

@deftypefun {gsl_multiroot_fsolver *} gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * @var{T}, size_t @var{n})
This function returns a pointer to a a newly allocated instance of a
solver of type @var{T} for a system of @var{n} dimensions.
For example, the following code creates an instance of a hybrid solver, 
to solve a 3-dimensional system of equations.

@example
const gsl_multiroot_fsolver_type * T 
    = gsl_multiroot_fsolver_hybrid;
gsl_multiroot_fsolver * s 
    = gsl_multiroot_fsolver_alloc (T, 3);
@end example
@noindent
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_multiroot_fdfsolver *} gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * @var{T}, size_t @var{n})
This function returns a pointer to a a newly allocated instance of a
derivative solver of type @var{T} for a system of @var{n} dimensions.
For example, the following code creates an instance of a Newton-Raphson solver,
for a 2-dimensional system of equations.

@example
const gsl_multiroot_fdfsolver_type * T 
    = gsl_multiroot_fdfsolver_newton;
gsl_multiroot_fdfsolver * s = 
    gsl_multiroot_fdfsolver_alloc (T, 2);
@end example
@noindent
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_multiroot_fsolver_set (gsl_multiroot_fsolver * @var{s}, gsl_multiroot_function * @var{f}, gsl_vector * @var{x})
This function sets, or resets, an existing solver @var{s} to use the
function @var{f} and the initial guess @var{x}.
@end deftypefun

@deftypefun int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * @var{s}, gsl_function_fdf * @var{fdf}, gsl_vector * @var{x})
This function sets, or resets, an existing solver @var{s} to use the
function and derivative @var{fdf} and the initial guess @var{x}.
@end deftypefun

@deftypefun void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * @var{s})
@deftypefunx void gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * @var{s})
These functions free all the memory associated with the solver @var{s}.
@end deftypefun

@deftypefun {const char *} gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * @var{s})
@deftypefunx {const char *} gsl_multiroot_fdfsolver_name (const gsl_multiroot_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_multiroot_fdfsolver_name (s));
@end example
@noindent
would print something like @code{s is a 'newton' solver}.
@end deftypefun

@node Providing the multidimensional system of equations to solve
@section Providing the function to solve
@cindex multidimensional root finding, providing a function to solve

You must provide @math{n} functions of @math{n} variables for the root
finders to operate on.  In order to allow for general parameters the
functions are defined by the following data types:

@deftp {Data Type} gsl_multiroot_function 
This data type defines a general system of functions with parameters.

@table @code
@item int (* f) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f})
this function should store the vector result
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} in @var{f} for argument @var{x} and parameters @var{params},
returning an appropriate error code if the function cannot be computed.

@item size_t @var{n}
the dimension of the system, i.e. the number of components of the
vectors @var{x} and @var{f}.

@item void * @var{params}
a pointer to the parameters of the function.
@end table
@end deftp

@noindent
Here is an example using Powell's test function,

@tex
\beforedisplay
$$
f_1(x) = A x_0 x_1 - 1,
f_2(x) = \exp(-x_0) + \exp(-x_1) - (1 + 1/A)
$$
\afterdisplay
@end tex
@ifinfo
@example
f_1(x) = A x_0 x_1 - 1,
f_2(x) = exp(-x_0) + exp(-x_1) - (1 + 1/A)
@end example
@end ifinfo
@noindent
with @math{A = 10^4}.  The following code defines a
@code{gsl_multiroot_function} system @code{F} which you could pass to a
solver:

@example
struct powell_params @{ double A; @};

int
powell (gsl_vector * x, void * p, gsl_vector * f) @{
   struct powell_params * params 
     = *(struct powell_params *)p;
   double A = (params->A);
   double x0 = gsl_vector_get(x,0);
   double x1 = gsl_vector_get(x,1);

   gsl_vector_set (f, 0, A * x0 * x1 - 1)
   gsl_vector_set (f, 1, (exp(-x0) + exp(-x1) 
                          - (1.0 + 1.0/A)))
   return GSL_SUCCESS
@}

gsl_multiroot_function F;
struct powell_params params = @{ 10000.0 @};

F.function = &powell;
F.n = 2;
F.params = &params;
@end example

@deftp {Data Type} gsl_multiroot_function_fdf
This data type defines a general system of functions with parameters and
the corresponding Jacobian matrix of derivatives,

@table @code
@item int (* f) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f})
this function should store the vector result
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} in @var{f} for argument @var{x} and parameters @var{params},
returning an appropriate error code if the function cannot be computed.

@item int (* df) (const gsl_vector * @var{x}, void * @var{params}, gsl_matrix * @var{J})
this function should store the @var{n}-by-@var{n} matrix result
@c{$J_{ij} = \partial f_i(x,\hbox{\it params}) / \partial x_j$}
@math{J_ij = d f_i(x,params) / d x_j} in @var{J} for argument @var{x} 
and parameters @var{params}, returning an appropriate error code if the
function cannot be computed.

@item int (* fdf) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f}, gsl_matrix * @var{J})
This function should set the values of the @var{f} and @var{J} as above,
for arguments @var{x} and parameters @var{params}.  This function provides
an optimization of the separate functions for @math{f(x)} and @math{J(x)} -- 
it is always faster to compute the function and its derivative at the
same time. 

@item size_t @var{n}
the dimension of the system, i.e. the number of components of the
vectors @var{x} and @var{f}.

@item void * @var{params}
a pointer to the parameters of the function.
@end table
@end deftp

@noindent
The example of Powell's test function defined above can be extended to
include analytic derivatives using the following code,

@example
int
powell_df (gsl_vector * x, void * p, gsl_matrix * J) 
@{
   struct powell_params * params 
     = *(struct powell_params *)p;
   double A = (params->A);
   double x0 = gsl_vector_get(x,0);
   double x1 = gsl_vector_get(x,1);
   gsl_vector_set (J, 0, 0, A * x1)
   gsl_vector_set (J, 0, 1, A * x0)
   gsl_vector_set (J, 1, 0, -exp(-x0))
   gsl_vector_set (J, 1, 1, -exp(-x1))
   return GSL_SUCCESS
@}

int
powell_fdf (gsl_vector * x, void * p, 
            gsl_matrix * f, gsl_matrix * J) @{
   struct powell_params * params 
     = *(struct powell_params *)p;
   double A = (params->A);
   double x0 = gsl_vector_get(x,0);
   double x1 = gsl_vector_get(x,1);

   double u0 = exp(-x0);
   double u1 = exp(-x1);

   gsl_vector_set (f, 0, A * x0 * x1 - 1)
   gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A))

   gsl_vector_set (J, 0, 0, A * x1)
   gsl_vector_set (J, 0, 1, A * x0)
   gsl_vector_set (J, 1, 0, -u0)
   gsl_vector_set (J, 1, 1, -u1)
   return GSL_SUCCESS
@}

gsl_multiroot_function_fdf FDF;

FDF.f = &powell_f;
FDF.df = &powell_df;
FDF.fdf = &powell_fdf;
FDF.n = 2;
FDF.params = 0;
@end example
@noindent
Note that the function @code{powell_fdf} is able to reuse existing terms
from the function when calculating the Jacobian, thus saving time.

@node Iteration of the multidimensional solver
@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