1 @cindex solving nonlinear systems of equations
2 @cindex nonlinear systems of equations, solution of
3 @cindex systems of equations, nonlinear
5 This chapter describes functions for multidimensional root-finding
6 (solving nonlinear systems with @math{n} equations in @math{n}
7 unknowns). The library provides low level components for a variety of
8 iterative solvers and convergence tests. These can be combined by the
9 user to achieve the desired solution, with full access to the
10 intermediate steps of the iteration. Each class of methods uses the
11 same framework, so that you can switch between solvers at runtime
12 without needing to recompile your program. Each instance of a solver
13 keeps track of its own state, allowing the solvers to be used in
14 multi-threaded programs. The solvers are based on the original Fortran
17 The header file @file{gsl_multiroots.h} contains prototypes for the
18 multidimensional root finding functions and related declarations.
21 * Overview of Multidimensional Root Finding::
22 * Initializing the Multidimensional Solver::
23 * Providing the multidimensional system of equations to solve::
24 * Iteration of the multidimensional solver::
25 * Search Stopping Parameters for the multidimensional solver::
26 * Algorithms using Derivatives::
27 * Algorithms without Derivatives::
28 * Example programs for Multidimensional Root finding::
29 * References and Further Reading for Multidimensional Root Finding::
32 @node Overview of Multidimensional Root Finding
34 @cindex multidimensional root finding, overview
36 The problem of multidimensional root finding requires the simultaneous
37 solution of @math{n} equations, @math{f_i}, in @math{n} variables,
42 f_i (x_1, \dots, x_n) = 0 \qquad\hbox{for}~i = 1 \dots n.
49 f_i (x_1, ..., x_n) = 0 for i = 1 ... n.
54 In general there are no bracketing methods available for @math{n}
55 dimensional systems, and no way of knowing whether any solutions
56 exist. All algorithms proceed from an initial guess using a variant of
61 x \to x' = x - J^{-1} f(x)
68 x -> x' = x - J^@{-1@} f(x)
73 where @math{x}, @math{f} are vector quantities and @math{J} is the
74 Jacobian matrix @c{$J_{ij} = \partial f_i / \partial x_j$}
75 @math{J_@{ij@} = d f_i / d x_j}.
76 Additional strategies can be used to enlarge the region of
77 convergence. These include requiring a decrease in the norm @math{|f|} on
78 each step proposed by Newton's method, or taking steepest-descent steps in
79 the direction of the negative gradient of @math{|f|}.
81 Several root-finding algorithms are available within a single framework.
82 The user provides a high-level driver for the algorithms, and the
83 library provides the individual functions necessary for each of the
84 steps. There are three main phases of the iteration. The steps are,
88 initialize solver state, @var{s}, for algorithm @var{T}
91 update @var{s} using the iteration @var{T}
94 test @var{s} for convergence, and repeat iteration if necessary
98 The evaluation of the Jacobian matrix can be problematic, either because
99 programming the derivatives is intractable or because computation of the
100 @math{n^2} terms of the matrix becomes too expensive. For these reasons
101 the algorithms provided by the library are divided into two classes according
102 to whether the derivatives are available or not.
104 The state for solvers with an analytic Jacobian matrix is held in a
105 @code{gsl_multiroot_fdfsolver} struct. The updating procedure requires
106 both the function and its derivatives to be supplied by the user.
108 The state for solvers which do not use an analytic Jacobian matrix is
109 held in a @code{gsl_multiroot_fsolver} struct. The updating procedure
110 uses only function evaluations (not derivatives). The algorithms
111 estimate the matrix @math{J} or @c{$J^{-1}$}
112 @math{J^@{-1@}} by approximate methods.
114 @node Initializing the Multidimensional Solver
115 @section Initializing the Solver
117 The following functions initialize a multidimensional solver, either
118 with or without derivatives. The solver itself depends only on the
119 dimension of the problem and the algorithm and can be reused for
122 @deftypefun {gsl_multiroot_fsolver *} gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * @var{T}, size_t @var{n})
123 This function returns a pointer to a newly allocated instance of a
124 solver of type @var{T} for a system of @var{n} dimensions.
125 For example, the following code creates an instance of a hybrid solver,
126 to solve a 3-dimensional system of equations.
129 const gsl_multiroot_fsolver_type * T
130 = gsl_multiroot_fsolver_hybrid;
131 gsl_multiroot_fsolver * s
132 = gsl_multiroot_fsolver_alloc (T, 3);
136 If there is insufficient memory to create the solver then the function
137 returns a null pointer and the error handler is invoked with an error
138 code of @code{GSL_ENOMEM}.
141 @deftypefun {gsl_multiroot_fdfsolver *} gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * @var{T}, size_t @var{n})
142 This function returns a pointer to a newly allocated instance of a
143 derivative solver of type @var{T} for a system of @var{n} dimensions.
144 For example, the following code creates an instance of a Newton-Raphson solver,
145 for a 2-dimensional system of equations.
148 const gsl_multiroot_fdfsolver_type * T
149 = gsl_multiroot_fdfsolver_newton;
150 gsl_multiroot_fdfsolver * s =
151 gsl_multiroot_fdfsolver_alloc (T, 2);
155 If there is insufficient memory to create the solver then the function
156 returns a null pointer and the error handler is invoked with an error
157 code of @code{GSL_ENOMEM}.
160 @deftypefun int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * @var{s}, gsl_multiroot_function * @var{f}, const gsl_vector * @var{x})
161 @deftypefunx int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * @var{s}, gsl_multiroot_function_fdf * @var{fdf}, const gsl_vector * @var{x})
162 These functions set, or reset, an existing solver @var{s} to use the
163 function @var{f} or function and derivative @var{fdf}, and the initial
164 guess @var{x}. Note that the initial position is copied from @var{x}, this
165 argument is not modified by subsequent iterations.
168 @deftypefun void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * @var{s})
169 @deftypefunx void gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * @var{s})
170 These functions free all the memory associated with the solver @var{s}.
173 @deftypefun {const char *} gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * @var{s})
174 @deftypefunx {const char *} gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * @var{s})
175 These functions return a pointer to the name of the solver. For example,
178 printf ("s is a '%s' solver\n",
179 gsl_multiroot_fdfsolver_name (s));
183 would print something like @code{s is a 'newton' solver}.
186 @node Providing the multidimensional system of equations to solve
187 @section Providing the function to solve
188 @cindex multidimensional root finding, providing a function to solve
190 You must provide @math{n} functions of @math{n} variables for the root
191 finders to operate on. In order to allow for general parameters the
192 functions are defined by the following data types:
194 @deftp {Data Type} gsl_multiroot_function
195 This data type defines a general system of functions with parameters.
198 @item int (* f) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f})
199 this function should store the vector result
200 @c{$f(x,\hbox{\it params})$}
201 @math{f(x,params)} in @var{f} for argument @var{x} and parameters @var{params},
202 returning an appropriate error code if the function cannot be computed.
205 the dimension of the system, i.e. the number of components of the
206 vectors @var{x} and @var{f}.
209 a pointer to the parameters of the function.
214 Here is an example using Powell's test function,
218 f_1(x) = A x_0 x_1 - 1,
219 f_2(x) = \exp(-x_0) + \exp(-x_1) - (1 + 1/A)
226 f_1(x) = A x_0 x_1 - 1,
227 f_2(x) = exp(-x_0) + exp(-x_1) - (1 + 1/A)
232 with @math{A = 10^4}. The following code defines a
233 @code{gsl_multiroot_function} system @code{F} which you could pass to a
237 struct powell_params @{ double A; @};
240 powell (gsl_vector * x, void * p, gsl_vector * f) @{
241 struct powell_params * params
242 = *(struct powell_params *)p;
243 const double A = (params->A);
244 const double x0 = gsl_vector_get(x,0);
245 const double x1 = gsl_vector_get(x,1);
247 gsl_vector_set (f, 0, A * x0 * x1 - 1);
248 gsl_vector_set (f, 1, (exp(-x0) + exp(-x1)
253 gsl_multiroot_function F;
254 struct powell_params params = @{ 10000.0 @};
261 @deftp {Data Type} gsl_multiroot_function_fdf
262 This data type defines a general system of functions with parameters and
263 the corresponding Jacobian matrix of derivatives,
266 @item int (* f) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f})
267 this function should store the vector result
268 @c{$f(x,\hbox{\it params})$}
269 @math{f(x,params)} in @var{f} for argument @var{x} and parameters @var{params},
270 returning an appropriate error code if the function cannot be computed.
272 @item int (* df) (const gsl_vector * @var{x}, void * @var{params}, gsl_matrix * @var{J})
273 this function should store the @var{n}-by-@var{n} matrix result
274 @c{$J_{ij} = \partial f_i(x,\hbox{\it params}) / \partial x_j$}
275 @math{J_ij = d f_i(x,params) / d x_j} in @var{J} for argument @var{x}
276 and parameters @var{params}, returning an appropriate error code if the
277 function cannot be computed.
279 @item int (* fdf) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f}, gsl_matrix * @var{J})
280 This function should set the values of the @var{f} and @var{J} as above,
281 for arguments @var{x} and parameters @var{params}. This function
282 provides an optimization of the separate functions for @math{f(x)} and
283 @math{J(x)}---it is always faster to compute the function and its
284 derivative at the same time.
287 the dimension of the system, i.e. the number of components of the
288 vectors @var{x} and @var{f}.
291 a pointer to the parameters of the function.
296 The example of Powell's test function defined above can be extended to
297 include analytic derivatives using the following code,
301 powell_df (gsl_vector * x, void * p, gsl_matrix * J)
303 struct powell_params * params
304 = *(struct powell_params *)p;
305 const double A = (params->A);
306 const double x0 = gsl_vector_get(x,0);
307 const double x1 = gsl_vector_get(x,1);
308 gsl_matrix_set (J, 0, 0, A * x1);
309 gsl_matrix_set (J, 0, 1, A * x0);
310 gsl_matrix_set (J, 1, 0, -exp(-x0));
311 gsl_matrix_set (J, 1, 1, -exp(-x1));
316 powell_fdf (gsl_vector * x, void * p,
317 gsl_matrix * f, gsl_matrix * J) @{
318 struct powell_params * params
319 = *(struct powell_params *)p;
320 const double A = (params->A);
321 const double x0 = gsl_vector_get(x,0);
322 const double x1 = gsl_vector_get(x,1);
324 const double u0 = exp(-x0);
325 const double u1 = exp(-x1);
327 gsl_vector_set (f, 0, A * x0 * x1 - 1);
328 gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A));
330 gsl_matrix_set (J, 0, 0, A * x1);
331 gsl_matrix_set (J, 0, 1, A * x0);
332 gsl_matrix_set (J, 1, 0, -u0);
333 gsl_matrix_set (J, 1, 1, -u1);
337 gsl_multiroot_function_fdf FDF;
341 FDF.fdf = &powell_fdf;
347 Note that the function @code{powell_fdf} is able to reuse existing terms
348 from the function when calculating the Jacobian, thus saving time.
350 @node Iteration of the multidimensional solver
353 The following functions drive the iteration of each algorithm. Each
354 function performs one iteration to update the state of any solver of the
355 corresponding type. The same functions work for all solvers so that
356 different methods can be substituted at runtime without modifications to
359 @deftypefun int gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * @var{s})
360 @deftypefunx int gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * @var{s})
361 These functions perform a single iteration of the solver @var{s}. If the
362 iteration encounters an unexpected problem then an error code will be
367 the iteration encountered a singular point where the function or its
368 derivative evaluated to @code{Inf} or @code{NaN}.
371 the iteration is not making any progress, preventing the algorithm from
376 The solver maintains a current best estimate of the root @code{s->x}
377 and its function value @code{s->f} at all times. This information can
378 be accessed with the following auxiliary functions,
380 @deftypefun {gsl_vector *} gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * @var{s})
381 @deftypefunx {gsl_vector *} gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * @var{s})
382 These functions return the current estimate of the root for the solver @var{s}, given by @code{s->x}.
385 @deftypefun {gsl_vector *} gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * @var{s})
386 @deftypefunx {gsl_vector *} gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * @var{s})
387 These functions return the function value @math{f(x)} at the current
388 estimate of the root for the solver @var{s}, given by @code{s->f}.
391 @deftypefun {gsl_vector *} gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * @var{s})
392 @deftypefunx {gsl_vector *} gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * @var{s})
393 These functions return the last step @math{dx} taken by the solver
394 @var{s}, given by @code{s->dx}.
397 @node Search Stopping Parameters for the multidimensional solver
398 @section Search Stopping Parameters
399 @cindex root finding, stopping parameters
401 A root finding procedure should stop when one of the following conditions is
406 A multidimensional root has been found to within the user-specified precision.
409 A user-specified maximum number of iterations has been reached.
412 An error has occurred.
416 The handling of these conditions is under user control. The functions
417 below allow the user to test the precision of the current result in
418 several standard ways.
420 @deftypefun int gsl_multiroot_test_delta (const gsl_vector * @var{dx}, const gsl_vector * @var{x}, double @var{epsabs}, double @var{epsrel})
422 This function tests for the convergence of the sequence by comparing the
423 last step @var{dx} with the absolute error @var{epsabs} and relative
424 error @var{epsrel} to the current position @var{x}. The test returns
425 @code{GSL_SUCCESS} if the following condition is achieved,
429 |dx_i| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, |x_i|
436 |dx_i| < epsabs + epsrel |x_i|
441 for each component of @var{x} and returns @code{GSL_CONTINUE} otherwise.
444 @cindex residual, in nonlinear systems of equations
445 @deftypefun int gsl_multiroot_test_residual (const gsl_vector * @var{f}, double @var{epsabs})
446 This function tests the residual value @var{f} against the absolute
447 error bound @var{epsabs}. The test returns @code{GSL_SUCCESS} if the
448 following condition is achieved,
452 \sum_i |f_i| < \hbox{\it epsabs}
459 \sum_i |f_i| < epsabs
464 and returns @code{GSL_CONTINUE} otherwise. This criterion is suitable
465 for situations where the precise location of the root, @math{x}, is
466 unimportant provided a value can be found where the residual is small
470 @comment ============================================================
472 @node Algorithms using Derivatives
473 @section Algorithms using Derivatives
475 The root finding algorithms described in this section make use of both
476 the function and its derivative. They require an initial guess for the
477 location of the root, but there is no absolute guarantee of
478 convergence---the function must be suitable for this technique and the
479 initial guess must be sufficiently close to the root for it to work.
480 When the conditions are satisfied then convergence is quadratic.
483 @comment ============================================================
484 @cindex HYBRID algorithms for nonlinear systems
485 @deffn {Derivative Solver} gsl_multiroot_fdfsolver_hybridsj
486 @cindex HYBRIDSJ algorithm
487 @cindex MINPACK, minimization algorithms
488 This is a modified version of Powell's Hybrid method as implemented in
489 the @sc{hybrj} algorithm in @sc{minpack}. Minpack was written by Jorge
490 J. Mor@'e, Burton S. Garbow and Kenneth E. Hillstrom. The Hybrid
491 algorithm retains the fast convergence of Newton's method but will also
492 reduce the residual when Newton's method is unreliable.
494 The algorithm uses a generalized trust region to keep each step under
495 control. In order to be accepted a proposed new position @math{x'} must
496 satisfy the condition @math{|D (x' - x)| < \delta}, where @math{D} is a
497 diagonal scaling matrix and @math{\delta} is the size of the trust
498 region. The components of @math{D} are computed internally, using the
499 column norms of the Jacobian to estimate the sensitivity of the residual
500 to each component of @math{x}. This improves the behavior of the
501 algorithm for badly scaled functions.
503 On each iteration the algorithm first determines the standard Newton
504 step by solving the system @math{J dx = - f}. If this step falls inside
505 the trust region it is used as a trial step in the next stage. If not,
506 the algorithm uses the linear combination of the Newton and gradient
507 directions which is predicted to minimize the norm of the function while
508 staying inside the trust region,
512 dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2.
519 dx = - \alpha J^@{-1@} f(x) - \beta \nabla |f(x)|^2.
524 This combination of Newton and gradient directions is referred to as a
527 The proposed step is now tested by evaluating the function at the
528 resulting point, @math{x'}. If the step reduces the norm of the function
529 sufficiently then it is accepted and size of the trust region is
530 increased. If the proposed step fails to improve the solution then the
531 size of the trust region is decreased and another trial step is
534 The speed of the algorithm is increased by computing the changes to the
535 Jacobian approximately, using a rank-1 update. If two successive
536 attempts fail to reduce the residual then the full Jacobian is
537 recomputed. The algorithm also monitors the progress of the solution
538 and returns an error if several steps fail to make any improvement,
542 the iteration is not making any progress, preventing the algorithm from
546 re-evaluations of the Jacobian indicate that the iteration is not
547 making any progress, preventing the algorithm from continuing.
552 @deffn {Derivative Solver} gsl_multiroot_fdfsolver_hybridj
553 @cindex HYBRIDJ algorithm
554 This algorithm is an unscaled version of @code{hybridsj}. The steps are
555 controlled by a spherical trust region @math{|x' - x| < \delta}, instead
556 of a generalized region. This can be useful if the generalized region
557 estimated by @code{hybridsj} is inappropriate.
561 @deffn {Derivative Solver} gsl_multiroot_fdfsolver_newton
562 @cindex Newton's method for systems of nonlinear equations
564 Newton's Method is the standard root-polishing algorithm. The algorithm
565 begins with an initial guess for the location of the solution. On each
566 iteration a linear approximation to the function @math{F} is used to
567 estimate the step which will zero all the components of the residual.
568 The iteration is defined by the following sequence,
572 x \to x' = x - J^{-1} f(x)
579 x -> x' = x - J^@{-1@} f(x)
584 where the Jacobian matrix @math{J} is computed from the derivative
585 functions provided by @var{f}. The step @math{dx} is obtained by solving
602 using LU decomposition.
605 @comment ============================================================
607 @deffn {Derivative Solver} gsl_multiroot_fdfsolver_gnewton
608 @cindex Modified Newton's method for nonlinear systems
609 @cindex Newton algorithm, globally convergent
610 This is a modified version of Newton's method which attempts to improve
611 global convergence by requiring every step to reduce the Euclidean norm
612 of the residual, @math{|f(x)|}. If the Newton step leads to an increase
613 in the norm then a reduced step of relative size,
617 t = (\sqrt{1 + 6 r} - 1) / (3 r)
624 t = (\sqrt(1 + 6 r) - 1) / (3 r)
629 is proposed, with @math{r} being the ratio of norms
630 @math{|f(x')|^2/|f(x)|^2}. This procedure is repeated until a suitable step
634 @comment ============================================================
636 @node Algorithms without Derivatives
637 @section Algorithms without Derivatives
639 The algorithms described in this section do not require any derivative
640 information to be supplied by the user. Any derivatives needed are
641 approximated by finite differences. Note that if the
642 finite-differencing step size chosen by these routines is inappropriate,
643 an explicit user-supplied numerical derivative can always be used with
644 the algorithms described in the previous section.
646 @deffn {Solver} gsl_multiroot_fsolver_hybrids
647 @cindex HYBRIDS algorithm, scaled without derivatives
648 This is a version of the Hybrid algorithm which replaces calls to the
649 Jacobian function by its finite difference approximation. The finite
650 difference approximation is computed using @code{gsl_multiroots_fdjac}
651 with a relative step size of @code{GSL_SQRT_DBL_EPSILON}. Note that
652 this step size will not be suitable for all problems.
655 @deffn {Solver} gsl_multiroot_fsolver_hybrid
656 @cindex HYBRID algorithm, unscaled without derivatives
657 This is a finite difference version of the Hybrid algorithm without
661 @comment ============================================================
663 @deffn {Solver} gsl_multiroot_fsolver_dnewton
665 @cindex Discrete Newton algorithm for multidimensional roots
666 @cindex Newton algorithm, discrete
668 The @dfn{discrete Newton algorithm} is the simplest method of solving a
669 multidimensional system. It uses the Newton iteration
673 x \to x - J^{-1} f(x)
680 x -> x - J^@{-1@} f(x)
685 where the Jacobian matrix @math{J} is approximated by taking finite
686 differences of the function @var{f}. The approximation scheme used by
687 this implementation is,
691 J_{ij} = (f_i(x + \delta_j) - f_i(x)) / \delta_j
698 J_@{ij@} = (f_i(x + \delta_j) - f_i(x)) / \delta_j
703 where @math{\delta_j} is a step of size @math{\sqrt\epsilon |x_j|} with
704 @math{\epsilon} being the machine precision
705 (@c{$\epsilon \approx 2.22 \times 10^{-16}$}
706 @math{\epsilon \approx 2.22 \times 10^-16}).
707 The order of convergence of Newton's algorithm is quadratic, but the
708 finite differences require @math{n^2} function evaluations on each
709 iteration. The algorithm may become unstable if the finite differences
710 are not a good approximation to the true derivatives.
713 @comment ============================================================
715 @deffn {Solver} gsl_multiroot_fsolver_broyden
716 @cindex Broyden algorithm for multidimensional roots
717 @cindex multidimensional root finding, Broyden algorithm
719 The @dfn{Broyden algorithm} is a version of the discrete Newton
720 algorithm which attempts to avoids the expensive update of the Jacobian
721 matrix on each iteration. The changes to the Jacobian are also
722 approximated, using a rank-1 update,
726 J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df
733 J^@{-1@} \to J^@{-1@} - (J^@{-1@} df - dx) dx^T J^@{-1@} / dx^T J^@{-1@} df
738 where the vectors @math{dx} and @math{df} are the changes in @math{x}
739 and @math{f}. On the first iteration the inverse Jacobian is estimated
740 using finite differences, as in the discrete Newton algorithm.
742 This approximation gives a fast update but is unreliable if the changes
743 are not small, and the estimate of the inverse Jacobian becomes worse as
744 time passes. The algorithm has a tendency to become unstable unless it
745 starts close to the root. The Jacobian is refreshed if this instability
746 is detected (consult the source for details).
748 This algorithm is included only for demonstration purposes, and is not
749 recommended for serious use.
752 @comment ============================================================
755 @node Example programs for Multidimensional Root finding
758 The multidimensional solvers are used in a similar way to the
759 one-dimensional root finding algorithms. This first example
760 demonstrates the @code{hybrids} scaled-hybrid algorithm, which does not
761 require derivatives. The program solves the Rosenbrock system of equations,
765 f_1 (x, y) = a (1 - x),~
766 f_2 (x, y) = b (y - x^2)
773 f_1 (x, y) = a (1 - x)
774 f_2 (x, y) = b (y - x^2)
779 with @math{a = 1, b = 10}. The solution of this system lies at
780 @math{(x,y) = (1,1)} in a narrow valley.
782 The first stage of the program is to define the system of equations,
787 #include <gsl/gsl_vector.h>
788 #include <gsl/gsl_multiroots.h>
797 rosenbrock_f (const gsl_vector * x, void *params,
800 double a = ((struct rparams *) params)->a;
801 double b = ((struct rparams *) params)->b;
803 const double x0 = gsl_vector_get (x, 0);
804 const double x1 = gsl_vector_get (x, 1);
806 const double y0 = a * (1 - x0);
807 const double y1 = b * (x1 - x0 * x0);
809 gsl_vector_set (f, 0, y0);
810 gsl_vector_set (f, 1, y1);
817 The main program begins by creating the function object @code{f}, with
818 the arguments @code{(x,y)} and parameters @code{(a,b)}. The solver
819 @code{s} is initialized to use this function, with the @code{hybrids}
826 const gsl_multiroot_fsolver_type *T;
827 gsl_multiroot_fsolver *s;
833 struct rparams p = @{1.0, 10.0@};
834 gsl_multiroot_function f = @{&rosenbrock_f, n, &p@};
836 double x_init[2] = @{-10.0, -5.0@};
837 gsl_vector *x = gsl_vector_alloc (n);
839 gsl_vector_set (x, 0, x_init[0]);
840 gsl_vector_set (x, 1, x_init[1]);
842 T = gsl_multiroot_fsolver_hybrids;
843 s = gsl_multiroot_fsolver_alloc (T, 2);
844 gsl_multiroot_fsolver_set (s, &f, x);
846 print_state (iter, s);
851 status = gsl_multiroot_fsolver_iterate (s);
853 print_state (iter, s);
855 if (status) /* check if solver is stuck */
859 gsl_multiroot_test_residual (s->f, 1e-7);
861 while (status == GSL_CONTINUE && iter < 1000);
863 printf ("status = %s\n", gsl_strerror (status));
865 gsl_multiroot_fsolver_free (s);
872 Note that it is important to check the return status of each solver
873 step, in case the algorithm becomes stuck. If an error condition is
874 detected, indicating that the algorithm cannot proceed, then the error
875 can be reported to the user, a new starting point chosen or a different
878 The intermediate state of the solution is displayed by the following
879 function. The solver state contains the vector @code{s->x} which is the
880 current position, and the vector @code{s->f} with corresponding function
885 print_state (size_t iter, gsl_multiroot_fsolver * s)
887 printf ("iter = %3u x = % .3f % .3f "
888 "f(x) = % .3e % .3e\n",
890 gsl_vector_get (s->x, 0),
891 gsl_vector_get (s->x, 1),
892 gsl_vector_get (s->f, 0),
893 gsl_vector_get (s->f, 1));
898 Here are the results of running the program. The algorithm is started at
899 @math{(-10,-5)} far from the solution. Since the solution is hidden in
900 a narrow valley the earliest steps follow the gradient of the function
901 downhill, in an attempt to reduce the large value of the residual. Once
902 the root has been approximately located, on iteration 8, the Newton
903 behavior takes over and convergence is very rapid.
906 iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
907 iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
908 iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
909 iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
910 iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
911 iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
912 iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
913 iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
914 iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
915 iter = 9 x = 1.000 0.878 f(x) = 1.268e-10 -1.218e+00
916 iter = 10 x = 1.000 0.989 f(x) = 1.124e-11 -1.080e-01
917 iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00
922 Note that the algorithm does not update the location on every
923 iteration. Some iterations are used to adjust the trust-region
924 parameter, after trying a step which was found to be divergent, or to
925 recompute the Jacobian, when poor convergence behavior is detected.
927 The next example program adds derivative information, in order to
928 accelerate the solution. There are two derivative functions
929 @code{rosenbrock_df} and @code{rosenbrock_fdf}. The latter computes both
930 the function and its derivative simultaneously. This allows the
931 optimization of any common terms. For simplicity we substitute calls to
932 the separate @code{f} and @code{df} functions at this point in the code
937 rosenbrock_df (const gsl_vector * x, void *params,
940 const double a = ((struct rparams *) params)->a;
941 const double b = ((struct rparams *) params)->b;
943 const double x0 = gsl_vector_get (x, 0);
945 const double df00 = -a;
946 const double df01 = 0;
947 const double df10 = -2 * b * x0;
948 const double df11 = b;
950 gsl_matrix_set (J, 0, 0, df00);
951 gsl_matrix_set (J, 0, 1, df01);
952 gsl_matrix_set (J, 1, 0, df10);
953 gsl_matrix_set (J, 1, 1, df11);
959 rosenbrock_fdf (const gsl_vector * x, void *params,
960 gsl_vector * f, gsl_matrix * J)
962 rosenbrock_f (x, params, f);
963 rosenbrock_df (x, params, J);
970 The main program now makes calls to the corresponding @code{fdfsolver}
971 versions of the functions,
977 const gsl_multiroot_fdfsolver_type *T;
978 gsl_multiroot_fdfsolver *s;
984 struct rparams p = @{1.0, 10.0@};
985 gsl_multiroot_function_fdf f = @{&rosenbrock_f,
990 double x_init[2] = @{-10.0, -5.0@};
991 gsl_vector *x = gsl_vector_alloc (n);
993 gsl_vector_set (x, 0, x_init[0]);
994 gsl_vector_set (x, 1, x_init[1]);
996 T = gsl_multiroot_fdfsolver_gnewton;
997 s = gsl_multiroot_fdfsolver_alloc (T, n);
998 gsl_multiroot_fdfsolver_set (s, &f, x);
1000 print_state (iter, s);
1006 status = gsl_multiroot_fdfsolver_iterate (s);
1008 print_state (iter, s);
1013 status = gsl_multiroot_test_residual (s->f, 1e-7);
1015 while (status == GSL_CONTINUE && iter < 1000);
1017 printf ("status = %s\n", gsl_strerror (status));
1019 gsl_multiroot_fdfsolver_free (s);
1020 gsl_vector_free (x);
1026 The addition of derivative information to the @code{hybrids} solver does
1027 not make any significant difference to its behavior, since it able to
1028 approximate the Jacobian numerically with sufficient accuracy. To
1029 illustrate the behavior of a different derivative solver we switch to
1030 @code{gnewton}. This is a traditional Newton solver with the constraint
1031 that it scales back its step if the full step would lead ``uphill''. Here
1032 is the output for the @code{gnewton} algorithm,
1035 iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
1036 iter = 1 x = -4.231 -65.317 f(x) = 5.231e+00 -8.321e+02
1037 iter = 2 x = 1.000 -26.358 f(x) = -8.882e-16 -2.736e+02
1038 iter = 3 x = 1.000 1.000 f(x) = -2.220e-16 -4.441e-15
1043 The convergence is much more rapid, but takes a wide excursion out to
1044 the point @math{(-4.23,-65.3)}. This could cause the algorithm to go
1045 astray in a realistic application. The hybrid algorithm follows the
1046 downhill path to the solution more reliably.
1048 @node References and Further Reading for Multidimensional Root Finding
1049 @section References and Further Reading
1051 The original version of the Hybrid method is described in the following
1056 M.J.D. Powell, ``A Hybrid Method for Nonlinear Equations'' (Chap 6, p
1057 87--114) and ``A Fortran Subroutine for Solving systems of Nonlinear
1058 Algebraic Equations'' (Chap 7, p 115--161), in @cite{Numerical Methods for
1059 Nonlinear Algebraic Equations}, P. Rabinowitz, editor. Gordon and
1064 The following papers are also relevant to the algorithms described in
1069 J.J. Mor@'e, M.Y. Cosnard, ``Numerical Solution of Nonlinear Equations'',
1070 @cite{ACM Transactions on Mathematical Software}, Vol 5, No 1, (1979), p 64--85
1073 C.G. Broyden, ``A Class of Methods for Solving Nonlinear
1074 Simultaneous Equations'', @cite{Mathematics of Computation}, Vol 19 (1965),
1078 J.J. Mor@'e, B.S. Garbow, K.E. Hillstrom, ``Testing Unconstrained
1079 Optimization Software'', ACM Transactions on Mathematical Software, Vol
1080 7, No 1 (1981), p 17--41