| @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 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 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_multiroot_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 n |
| the dimension of the system, i.e. the number of components of the |
| vectors @var{x} and @var{f}. |
| |
| @item void * 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; |
| const double A = (params->A); |
| const double x0 = gsl_vector_get(x,0); |
| const 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.f = &powell; |
| F.n = 2; |
| F.params = ¶ms; |
| @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 n |
| the dimension of the system, i.e. the number of components of the |
| vectors @var{x} and @var{f}. |
| |
| @item void * 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; |
| const double A = (params->A); |
| const double x0 = gsl_vector_get(x,0); |
| const double x1 = gsl_vector_get(x,1); |
| gsl_matrix_set (J, 0, 0, A * x1); |
| gsl_matrix_set (J, 0, 1, A * x0); |
| gsl_matrix_set (J, 1, 0, -exp(-x0)); |
| gsl_matrix_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; |
| const double A = (params->A); |
| const double x0 = gsl_vector_get(x,0); |
| const double x1 = gsl_vector_get(x,1); |
| |
| const double u0 = exp(-x0); |
| const double u1 = exp(-x1); |
| |
| gsl_vector_set (f, 0, A * x0 * x1 - 1); |
| gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A)); |
| |
| gsl_matrix_set (J, 0, 0, A * x1); |
| gsl_matrix_set (J, 0, 1, A * x0); |
| gsl_matrix_set (J, 1, 0, -u0); |
| gsl_matrix_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 |
| 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_multiroot_fsolver_iterate (gsl_multiroot_fsolver * @var{s}) |
| @deftypefunx int gsl_multiroot_fdfsolver_iterate (gsl_multiroot_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_ENOPROG |
| the iteration is not making any progress, preventing the algorithm from |
| continuing. |
| @end table |
| @end deftypefun |
| |
| The solver maintains a current best estimate of the root at all times. |
| This information can be accessed with the following auxiliary functions, |
| |
| @deftypefun {gsl_vector *} gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * @var{s}) |
| @deftypefunx {gsl_vector *} gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * @var{s}) |
| These functions return the current estimate of the root for the solver @var{s}. |
| @end deftypefun |
| |
| @deftypefun {gsl_vector *} gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * @var{s}) |
| @deftypefunx {gsl_vector *} gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * @var{s}) |
| These functions return the function value @math{f(x)} at the current |
| estimate of the root for the solver @var{s}. |
| @end deftypefun |
| |
| @deftypefun {gsl_vector *} gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * @var{s}) |
| @deftypefunx {gsl_vector *} gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * @var{s}) |
| These functions return the last step @math{dx} taken by the solver |
| @var{s}. |
| @end deftypefun |
| |
| @node Search Stopping Parameters for the multidimensional solver |
| @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 multidimensional 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_multiroot_test_delta (const gsl_vector * @var{dx}, const gsl_vector * @var{x}, double @var{epsabs}, double @var{epsrel}) |
| |
| This function tests for the convergence of the sequence by comparing the |
| last step @var{dx} with the absolute error @var{epsabs} and relative |
| error @var{epsrel} to the current position @var{x}. The test returns |
| @code{GSL_SUCCESS} if the following condition is achieved, |
| @tex |
| \beforedisplay |
| $$ |
| |dx_i| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, |x_i| |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| |dx_i| < epsabs + epsrel |x_i| |
| @end example |
| |
| @end ifinfo |
| @noindent |
| for each component of @var{x} and returns @code{GSL_CONTINUE} otherwise. |
| @end deftypefun |
| |
| @cindex residual, in nonlinear systems of equations |
| @deftypefun int gsl_multiroot_test_residual (const gsl_vector * @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 |
| $$ |
| \sum_i |f_i| < \hbox{\it epsabs} |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| \sum_i |f_i| < 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 is small |
| enough. |
| @end deftypefun |
| |
| @comment ============================================================ |
| |
| @node Algorithms using Derivatives |
| @section Algorithms using Derivatives |
| |
| The root finding algorithms described in this section make use of both |
| the function and its derivative. They require an initial guess for the |
| location of the root, but 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 the conditions are satisfied then convergence is quadratic. |
| |
| |
| @comment ============================================================ |
| @cindex HYBRID algorithms for nonlinear systems |
| @deffn {Derivative Solver} gsl_multiroot_fdfsolver_hybridsj |
| @cindex HYBRIDSJ algorithm |
| @cindex MINPACK, minimization algorithms |
| This is a modified version of Powell's Hybrid method as implemented in |
| the @sc{hybrj} algorithm in @sc{minpack}. Minpack was written by Jorge |
| J. Mor@'e, Burton S. Garbow and Kenneth E. Hillstrom. The Hybrid |
| algorithm retains the fast convergence of Newton's method but will also |
| reduce the residual when Newton's method is unreliable. |
| |
| The algorithm uses a generalized trust region to keep each step under |
| control. In order to be accepted a proposed new position @math{x'} must |
| satisfy the condition @math{|D (x' - x)| < \delta}, where @math{D} is a |
| diagonal scaling matrix and @math{\delta} is the size of the trust |
| region. The components of @math{D} are computed internally, using the |
| column norms of the Jacobian to estimate the sensitivity of the residual |
| to each component of @math{x}. This improves the behavior of the |
| algorithm for badly scaled functions. |
| |
| On each iteration the algorithm first determines the standard Newton |
| step by solving the system @math{J dx = - f}. If this step falls inside |
| the trust region it is used as a trial step in the next stage. If not, |
| the algorithm uses the linear combination of the Newton and gradient |
| directions which is predicted to minimize the norm of the function while |
| staying inside the trust region, |
| @tex |
| \beforedisplay |
| $$ |
| dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2. |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| dx = - \alpha J^@{-1@} f(x) - \beta \nabla |f(x)|^2. |
| @end example |
| |
| @end ifinfo |
| @noindent |
| This combination of Newton and gradient directions is referred to as a |
| @dfn{dogleg step}. |
| |
| The proposed step is now tested by evaluating the function at the |
| resulting point, @math{x'}. If the step reduces the norm of the function |
| sufficiently then it is accepted and size of the trust region is |
| increased. If the proposed step fails to improve the solution then the |
| size of the trust region is decreased and another trial step is |
| computed. |
| |
| The speed of the algorithm is increased by computing the changes to the |
| Jacobian approximately, using a rank-1 update. If two successive |
| attempts fail to reduce the residual then the full Jacobian is |
| recomputed. The algorithm also monitors the progress of the solution |
| and returns an error if several steps fail to make any improvement, |
| |
| @table @code |
| @item GSL_ENOPROG |
| the iteration is not making any progress, preventing the algorithm from |
| continuing. |
| |
| @item GSL_ENOPROGJ |
| re-evaluations of the Jacobian indicate that the iteration is not |
| making any progress, preventing the algorithm from continuing. |
| @end table |
| |
| @end deffn |
| |
| @deffn {Derivative Solver} gsl_multiroot_fdfsolver_hybridj |
| @cindex HYBRIDJ algorithm |
| This algorithm is an unscaled version of @code{hybridsj}. The steps are |
| controlled by a spherical trust region @math{|x' - x| < \delta}, instead |
| of a generalized region. This can be useful if the generalized region |
| estimated by @code{hybridsj} is inappropriate. |
| @end deffn |
| |
| |
| @deffn {Derivative Solver} gsl_multiroot_fdfsolver_newton |
| @cindex Newton's method for systems of nonlinear equations |
| |
| Newton's Method is the standard root-polishing algorithm. The algorithm |
| begins with an initial guess for the location of the solution. On each |
| iteration a linear approximation to the function @math{F} is used to |
| estimate the step which will zero all the components of the residual. |
| The iteration is defined by the following sequence, |
| @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 the Jacobian matrix @math{J} is computed from the derivative |
| functions provided by @var{f}. The step @math{dx} is obtained by solving |
| the linear system, |
| @tex |
| \beforedisplay |
| $$ |
| J \,dx = - f(x) |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| J dx = - f(x) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| using LU decomposition. |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Derivative Solver} gsl_multiroot_fdfsolver_gnewton |
| @cindex Modified Newton's method for nonlinear systems |
| @cindex Newton algorithm, globally convergent |
| This is a modified version of Newton's method which attempts to improve |
| global convergence by requiring every step to reduce the Euclidean norm |
| of the residual, @math{|f(x)|}. If the Newton step leads to an increase |
| in the norm then a reduced step of relative size, |
| @tex |
| \beforedisplay |
| $$ |
| t = (\sqrt{1 + 6 r} - 1) / (3 r) |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| t = (\sqrt(1 + 6 r) - 1) / (3 r) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| is proposed, with @math{r} being the ratio of norms |
| @math{|f(x')|^2/|f(x)|^2}. This procedure is repeated until a suitable step |
| size is found. |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @node Algorithms without Derivatives |
| @section Algorithms without Derivatives |
| |
| The algorithms described in this section do not require any derivative |
| information to be supplied by the user. Any derivatives needed are |
| approximated by finite differences. Note that if the |
| finite-differencing step size chosen by these routines is inappropriate, |
| an explicit user-supplied numerical derivative can always be used with |
| the algorithms described in the previous section. |
| |
| @deffn {Solver} gsl_multiroot_fsolver_hybrids |
| @cindex HYBRIDS algorithm, scaled without derivatives |
| This is a version of the Hybrid algorithm which replaces calls to the |
| Jacobian function by its finite difference approximation. The finite |
| difference approximation is computed using @code{gsl_multiroots_fdjac} |
| with a relative step size of @code{GSL_SQRT_DBL_EPSILON}. Note that |
| this step size will not be suitable for all problems. |
| @end deffn |
| |
| @deffn {Solver} gsl_multiroot_fsolver_hybrid |
| @cindex HYBRID algorithm, unscaled without derivatives |
| This is a finite difference version of the Hybrid algorithm without |
| internal scaling. |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Solver} gsl_multiroot_fsolver_dnewton |
| |
| @cindex Discrete Newton algorithm for multidimensional roots |
| @cindex Newton algorithm, discrete |
| |
| The @dfn{discrete Newton algorithm} is the simplest method of solving a |
| multidimensional system. It uses the Newton iteration |
| @tex |
| \beforedisplay |
| $$ |
| x \to x - J^{-1} f(x) |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| x -> x - J^@{-1@} f(x) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| where the Jacobian matrix @math{J} is approximated by taking finite |
| differences of the function @var{f}. The approximation scheme used by |
| this implementation is, |
| @tex |
| \beforedisplay |
| $$ |
| J_{ij} = (f_i(x + \delta_j) - f_i(x)) / \delta_j |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| J_@{ij@} = (f_i(x + \delta_j) - f_i(x)) / \delta_j |
| @end example |
| |
| @end ifinfo |
| @noindent |
| where @math{\delta_j} is a step of size @math{\sqrt\epsilon |x_j|} with |
| @math{\epsilon} being the machine precision |
| (@c{$\epsilon \approx 2.22 \times 10^{-16}$} |
| @math{\epsilon \approx 2.22 \times 10^-16}). |
| The order of convergence of Newton's algorithm is quadratic, but the |
| finite differences require @math{n^2} function evaluations on each |
| iteration. The algorithm may become unstable if the finite differences |
| are not a good approximation to the true derivatives. |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Solver} gsl_multiroot_fsolver_broyden |
| @cindex Broyden algorithm for multidimensional roots |
| @cindex multidimensional root finding, Broyden algorithm |
| |
| The @dfn{Broyden algorithm} is a version of the discrete Newton |
| algorithm which attempts to avoids the expensive update of the Jacobian |
| matrix on each iteration. The changes to the Jacobian are also |
| approximated, using a rank-1 update, |
| @tex |
| \beforedisplay |
| $$ |
| J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| J^@{-1@} \to J^@{-1@} - (J^@{-1@} df - dx) dx^T J^@{-1@} / dx^T J^@{-1@} df |
| @end example |
| |
| @end ifinfo |
| @noindent |
| where the vectors @math{dx} and @math{df} are the changes in @math{x} |
| and @math{f}. On the first iteration the inverse Jacobian is estimated |
| using finite differences, as in the discrete Newton algorithm. |
| |
| This approximation gives a fast update but is unreliable if the changes |
| are not small, and the estimate of the inverse Jacobian becomes worse as |
| time passes. The algorithm has a tendency to become unstable unless it |
| starts close to the root. The Jacobian is refreshed if this instability |
| is detected (consult the source for details). |
| |
| This algorithm is included only for demonstration purposes, and is not |
| recommended for serious use. |
| @end deffn |
| |
| @comment ============================================================ |
| |
| |
| @node Example programs for Multidimensional Root finding |
| @section Examples |
| |
| The multidimensional solvers are used in a similar way to the |
| one-dimensional root finding algorithms. This first example |
| demonstrates the @code{hybrids} scaled-hybrid algorithm, which does not |
| require derivatives. The program solves the Rosenbrock system of equations, |
| @tex |
| \beforedisplay |
| $$ |
| f_1 (x, y) = a (1 - x),~ |
| f_2 (x, y) = b (y - x^2) |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| f_1 (x, y) = a (1 - x) |
| f_2 (x, y) = b (y - x^2) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| with @math{a = 1, b = 10}. The solution of this system lies at |
| @math{(x,y) = (1,1)} in a narrow valley. |
| |
| The first stage of the program is to define the system of equations, |
| |
| @example |
| #include <stdlib.h> |
| #include <stdio.h> |
| #include <gsl/gsl_vector.h> |
| #include <gsl/gsl_multiroots.h> |
| |
| struct rparams |
| @{ |
| double a; |
| double b; |
| @}; |
| |
| int |
| rosenbrock_f (const gsl_vector * x, void *params, |
| gsl_vector * f) |
| @{ |
| double a = ((struct rparams *) params)->a; |
| double b = ((struct rparams *) params)->b; |
| |
| const double x0 = gsl_vector_get (x, 0); |
| const double x1 = gsl_vector_get (x, 1); |
| |
| const double y0 = a * (1 - x0); |
| const double y1 = b * (x1 - x0 * x0); |
| |
| gsl_vector_set (f, 0, y0); |
| gsl_vector_set (f, 1, y1); |
| |
| return GSL_SUCCESS; |
| @} |
| @end example |
| |
| @noindent |
| The main program begins by creating the function object @code{f}, with |
| the arguments @code{(x,y)} and parameters @code{(a,b)}. The solver |
| @code{s} is initialized to use this function, with the @code{hybrids} |
| method. |
| |
| @example |
| int |
| main (void) |
| @{ |
| const gsl_multiroot_fsolver_type *T; |
| gsl_multiroot_fsolver *s; |
| |
| int status; |
| size_t i, iter = 0; |
| |
| const size_t n = 2; |
| struct rparams p = @{1.0, 10.0@}; |
| gsl_multiroot_function f = @{&rosenbrock_f, n, &p@}; |
| |
| double x_init[2] = @{-10.0, -5.0@}; |
| gsl_vector *x = gsl_vector_alloc (n); |
| |
| gsl_vector_set (x, 0, x_init[0]); |
| gsl_vector_set (x, 1, x_init[1]); |
| |
| T = gsl_multiroot_fsolver_hybrids; |
| s = gsl_multiroot_fsolver_alloc (T, 2); |
| gsl_multiroot_fsolver_set (s, &f, x); |
| |
| print_state (iter, s); |
| |
| do |
| @{ |
| iter++; |
| status = gsl_multiroot_fsolver_iterate (s); |
| |
| print_state (iter, s); |
| |
| if (status) /* check if solver is stuck */ |
| break; |
| |
| status = |
| gsl_multiroot_test_residual (s->f, 1e-7); |
| @} |
| while (status == GSL_CONTINUE && iter < 1000); |
| |
| printf ("status = %s\n", gsl_strerror (status)); |
| |
| gsl_multiroot_fsolver_free (s); |
| gsl_vector_free (x); |
| return 0; |
| @} |
| @end example |
| |
| @noindent |
| Note that it is important to check the return status of each solver |
| step, in case the algorithm becomes stuck. If an error condition is |
| detected, indicating that the algorithm cannot proceed, then the error |
| can be reported to the user, a new starting point chosen or a different |
| algorithm used. |
| |
| The intermediate state of the solution is displayed by the following |
| function. The solver state contains the vector @code{s->x} which is the |
| current position, and the vector @code{s->f} with corresponding function |
| values. |
| |
| @example |
| int |
| print_state (size_t iter, gsl_multiroot_fsolver * s) |
| @{ |
| printf ("iter = %3u x = % .3f % .3f " |
| "f(x) = % .3e % .3e\n", |
| iter, |
| gsl_vector_get (s->x, 0), |
| gsl_vector_get (s->x, 1), |
| gsl_vector_get (s->f, 0), |
| gsl_vector_get (s->f, 1)); |
| @} |
| @end example |
| |
| @noindent |
| Here are the results of running the program. The algorithm is started at |
| @math{(-10,-5)} far from the solution. Since the solution is hidden in |
| a narrow valley the earliest steps follow the gradient of the function |
| downhill, in an attempt to reduce the large value of the residual. Once |
| the root has been approximately located, on iteration 8, the Newton |
| behavior takes over and convergence is very rapid. |
| |
| @smallexample |
| iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 |
| iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 |
| iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 |
| iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 |
| iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 |
| iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01 |
| iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01 |
| iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00 |
| iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00 |
| iter = 9 x = 1.000 0.878 f(x) = 1.268e-10 -1.218e+00 |
| iter = 10 x = 1.000 0.989 f(x) = 1.124e-11 -1.080e-01 |
| iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00 |
| status = success |
| @end smallexample |
| |
| @noindent |
| Note that the algorithm does not update the location on every |
| iteration. Some iterations are used to adjust the trust-region |
| parameter, after trying a step which was found to be divergent, or to |
| recompute the Jacobian, when poor convergence behavior is detected. |
| |
| The next example program adds derivative information, in order to |
| accelerate the solution. There are two derivative functions |
| @code{rosenbrock_df} and @code{rosenbrock_fdf}. The latter computes both |
| the function and its derivative simultaneously. This allows the |
| optimization of any common terms. For simplicity we substitute calls to |
| the separate @code{f} and @code{df} functions at this point in the code |
| below. |
| |
| @example |
| int |
| rosenbrock_df (const gsl_vector * x, void *params, |
| gsl_matrix * J) |
| @{ |
| const double a = ((struct rparams *) params)->a; |
| const double b = ((struct rparams *) params)->b; |
| |
| const double x0 = gsl_vector_get (x, 0); |
| |
| const double df00 = -a; |
| const double df01 = 0; |
| const double df10 = -2 * b * x0; |
| const double df11 = b; |
| |
| gsl_matrix_set (J, 0, 0, df00); |
| gsl_matrix_set (J, 0, 1, df01); |
| gsl_matrix_set (J, 1, 0, df10); |
| gsl_matrix_set (J, 1, 1, df11); |
| |
| return GSL_SUCCESS; |
| @} |
| |
| int |
| rosenbrock_fdf (const gsl_vector * x, void *params, |
| gsl_vector * f, gsl_matrix * J) |
| @{ |
| rosenbrock_f (x, params, f); |
| rosenbrock_df (x, params, J); |
| |
| return GSL_SUCCESS; |
| @} |
| @end example |
| |
| @noindent |
| The main program now makes calls to the corresponding @code{fdfsolver} |
| versions of the functions, |
| |
| @example |
| int |
| main (void) |
| @{ |
| const gsl_multiroot_fdfsolver_type *T; |
| gsl_multiroot_fdfsolver *s; |
| |
| int status; |
| size_t i, iter = 0; |
| |
| const size_t n = 2; |
| struct rparams p = @{1.0, 10.0@}; |
| gsl_multiroot_function_fdf f = @{&rosenbrock_f, |
| &rosenbrock_df, |
| &rosenbrock_fdf, |
| n, &p@}; |
| |
| double x_init[2] = @{-10.0, -5.0@}; |
| gsl_vector *x = gsl_vector_alloc (n); |
| |
| gsl_vector_set (x, 0, x_init[0]); |
| gsl_vector_set (x, 1, x_init[1]); |
| |
| T = gsl_multiroot_fdfsolver_gnewton; |
| s = gsl_multiroot_fdfsolver_alloc (T, n); |
| gsl_multiroot_fdfsolver_set (s, &f, x); |
| |
| print_state (iter, s); |
| |
| do |
| @{ |
| iter++; |
| |
| status = gsl_multiroot_fdfsolver_iterate (s); |
| |
| print_state (iter, s); |
| |
| if (status) |
| break; |
| |
| status = gsl_multiroot_test_residual (s->f, 1e-7); |
| @} |
| while (status == GSL_CONTINUE && iter < 1000); |
| |
| printf ("status = %s\n", gsl_strerror (status)); |
| |
| gsl_multiroot_fdfsolver_free (s); |
| gsl_vector_free (x); |
| return 0; |
| @} |
| @end example |
| |
| @noindent |
| The addition of derivative information to the @code{hybrids} solver does |
| not make any significant difference to its behavior, since it able to |
| approximate the Jacobian numerically with sufficient accuracy. To |
| illustrate the behavior of a different derivative solver we switch to |
| @code{gnewton}. This is a traditional Newton solver with the constraint |
| that it scales back its step if the full step would lead ``uphill''. Here |
| is the output for the @code{gnewton} algorithm, |
| |
| @smallexample |
| iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 |
| iter = 1 x = -4.231 -65.317 f(x) = 5.231e+00 -8.321e+02 |
| iter = 2 x = 1.000 -26.358 f(x) = -8.882e-16 -2.736e+02 |
| iter = 3 x = 1.000 1.000 f(x) = -2.220e-16 -4.441e-15 |
| status = success |
| @end smallexample |
| |
| @noindent |
| The convergence is much more rapid, but takes a wide excursion out to |
| the point @math{(-4.23,-65.3)}. This could cause the algorithm to go |
| astray in a realistic application. The hybrid algorithm follows the |
| downhill path to the solution more reliably. |
| |
| @node References and Further Reading for Multidimensional Root Finding |
| @section References and Further Reading |
| |
| The original version of the Hybrid method is described in the following |
| articles by Powell, |
| |
| @itemize @asis |
| @item |
| M.J.D. Powell, ``A Hybrid Method for Nonlinear Equations'' (Chap 6, p |
| 87--114) and ``A Fortran Subroutine for Solving systems of Nonlinear |
| Algebraic Equations'' (Chap 7, p 115--161), in @cite{Numerical Methods for |
| Nonlinear Algebraic Equations}, P. Rabinowitz, editor. Gordon and |
| Breach, 1970. |
| @end itemize |
| |
| @noindent |
| The following papers are also relevant to the algorithms described in |
| this section, |
| |
| @itemize @asis |
| @item |
| J.J. Mor@'e, M.Y. Cosnard, ``Numerical Solution of Nonlinear Equations'', |
| @cite{ACM Transactions on Mathematical Software}, Vol 5, No 1, (1979), p 64--85 |
| |
| @item |
| C.G. Broyden, ``A Class of Methods for Solving Nonlinear |
| Simultaneous Equations'', @cite{Mathematics of Computation}, Vol 19 (1965), |
| p 577--593 |
| |
| @item |
| J.J. Mor@'e, B.S. Garbow, K.E. Hillstrom, ``Testing Unconstrained |
| Optimization Software'', ACM Transactions on Mathematical Software, Vol |
| 7, No 1 (1981), p 17--41 |
| @end itemize |
| |