| @cindex root finding |
| @cindex zero finding |
| @cindex finding roots |
| @cindex finding zeros |
| @cindex roots |
| @cindex solving a nonlinear equation |
| @cindex nonlinear equation, solutions of |
| |
| This chapter describes routines for finding roots of arbitrary |
| one-dimensional functions. The library provides low level components |
| for a variety of iterative solvers and convergence tests. These can be |
| combined by the user to achieve the desired solution, with full access |
| to the intermediate steps of the iteration. Each class of methods uses |
| the same framework, so that you can switch between solvers at runtime |
| without needing to recompile your program. Each instance of a solver |
| keeps track of its own state, allowing the solvers to be used in |
| multi-threaded programs. |
| |
| The header file @file{gsl_roots.h} contains prototypes for the root |
| finding functions and related declarations. |
| |
| @menu |
| * Root Finding Overview:: |
| * Root Finding Caveats:: |
| * Initializing the Solver:: |
| * Providing the function to solve:: |
| * Search Bounds and Guesses:: |
| * Root Finding Iteration:: |
| * Search Stopping Parameters:: |
| * Root Bracketing Algorithms:: |
| * Root Finding Algorithms using Derivatives:: |
| * Root Finding Examples:: |
| * Root Finding References and Further Reading:: |
| @end menu |
| |
| @node Root Finding Overview |
| @section Overview |
| @cindex root finding, overview |
| |
| One-dimensional root finding algorithms can be divided into two classes, |
| @dfn{root bracketing} and @dfn{root polishing}. Algorithms which proceed |
| by bracketing a root are guaranteed to converge. Bracketing algorithms |
| begin with a bounded region known to contain a root. The size of this |
| bounded region is reduced, iteratively, until it encloses the root to a |
| desired tolerance. This provides a rigorous error estimate for the |
| location of the root. |
| |
| The technique of @dfn{root polishing} attempts to improve an initial |
| guess to the root. These algorithms converge only if started ``close |
| enough'' to a root, and sacrifice a rigorous error bound for speed. By |
| approximating the behavior of a function in the vicinity of a root they |
| attempt to find a higher order improvement of an initial guess. When the |
| behavior of the function is compatible with the algorithm and a good |
| initial guess is available a polishing algorithm can provide rapid |
| convergence. |
| |
| In GSL both types of algorithm are available in similar frameworks. The |
| user provides a high-level driver for the algorithms, and the library |
| provides the individual functions necessary for each of the steps. |
| There are three main phases of the iteration. The steps are, |
| |
| @itemize @bullet |
| @item |
| initialize solver state, @var{s}, for algorithm @var{T} |
| |
| @item |
| update @var{s} using the iteration @var{T} |
| |
| @item |
| test @var{s} for convergence, and repeat iteration if necessary |
| @end itemize |
| |
| @noindent |
| The state for bracketing solvers is held in a @code{gsl_root_fsolver} |
| struct. The updating procedure uses only function evaluations (not |
| derivatives). The state for root polishing solvers is held in a |
| @code{gsl_root_fdfsolver} struct. The updates require both the function |
| and its derivative (hence the name @code{fdf}) to be supplied by the |
| user. |
| |
| @node Root Finding Caveats |
| @section Caveats |
| @cindex root finding, caveats |
| |
| Note that root finding functions can only search for one root at a time. |
| When there are several roots in the search area, the first root to be |
| found will be returned; however it is difficult to predict which of the |
| roots this will be. @emph{In most cases, no error will be reported if |
| you try to find a root in an area where there is more than one.} |
| |
| Care must be taken when a function may have a multiple root (such as |
| @c{$f(x) = (x-x_0)^2$} |
| @math{f(x) = (x-x_0)^2} or |
| @c{$f(x) = (x-x_0)^3$} |
| @math{f(x) = (x-x_0)^3}). |
| It is not possible to use root-bracketing algorithms on |
| even-multiplicity roots. For these algorithms the initial interval must |
| contain a zero-crossing, where the function is negative at one end of |
| the interval and positive at the other end. Roots with even-multiplicity |
| do not cross zero, but only touch it instantaneously. Algorithms based |
| on root bracketing will still work for odd-multiplicity roots |
| (e.g. cubic, quintic, @dots{}). |
| Root polishing algorithms generally work with higher multiplicity roots, |
| but at a reduced rate of convergence. In these cases the @dfn{Steffenson |
| algorithm} can be used to accelerate the convergence of multiple roots. |
| |
| While it is not absolutely required that @math{f} have a root within the |
| search region, numerical root finding functions should not be used |
| haphazardly to check for the @emph{existence} of roots. There are better |
| ways to do this. Because it is easy to create situations where numerical |
| root finders can fail, it is a bad idea to throw a root finder at a |
| function you do not know much about. In general it is best to examine |
| the function visually by plotting before searching for a root. |
| |
| @node Initializing the Solver |
| @section Initializing the Solver |
| |
| @deftypefun {gsl_root_fsolver *} gsl_root_fsolver_alloc (const gsl_root_fsolver_type * @var{T}) |
| This function returns a pointer to a newly allocated instance of a |
| solver of type @var{T}. For example, the following code creates an |
| instance of a bisection solver, |
| |
| @example |
| const gsl_root_fsolver_type * T |
| = gsl_root_fsolver_bisection; |
| gsl_root_fsolver * s |
| = gsl_root_fsolver_alloc (T); |
| @end example |
| |
| If there is insufficient memory to create the solver then the function |
| returns a null pointer and the error handler is invoked with an error |
| code of @code{GSL_ENOMEM}. |
| @end deftypefun |
| |
| @deftypefun {gsl_root_fdfsolver *} gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * @var{T}) |
| This function returns a pointer to a newly allocated instance of a |
| derivative-based solver of type @var{T}. For example, the following |
| code creates an instance of a Newton-Raphson solver, |
| |
| @example |
| const gsl_root_fdfsolver_type * T |
| = gsl_root_fdfsolver_newton; |
| gsl_root_fdfsolver * s |
| = gsl_root_fdfsolver_alloc (T); |
| @end example |
| |
| If there is insufficient memory to create the solver then the function |
| returns a null pointer and the error handler is invoked with an error |
| code of @code{GSL_ENOMEM}. |
| @end deftypefun |
| |
| |
| @deftypefun int gsl_root_fsolver_set (gsl_root_fsolver * @var{s}, gsl_function * @var{f}, double @var{x_lower}, double @var{x_upper}) |
| This function initializes, or reinitializes, an existing solver @var{s} |
| to use the function @var{f} and the initial search interval |
| [@var{x_lower}, @var{x_upper}]. |
| @end deftypefun |
| |
| @deftypefun int gsl_root_fdfsolver_set (gsl_root_fdfsolver * @var{s}, gsl_function_fdf * @var{fdf}, double @var{root}) |
| This function initializes, or reinitializes, an existing solver @var{s} |
| to use the function and derivative @var{fdf} and the initial guess |
| @var{root}. |
| @end deftypefun |
| |
| @deftypefun void gsl_root_fsolver_free (gsl_root_fsolver * @var{s}) |
| @deftypefunx void gsl_root_fdfsolver_free (gsl_root_fdfsolver * @var{s}) |
| These functions free all the memory associated with the solver @var{s}. |
| @end deftypefun |
| |
| @deftypefun {const char *} gsl_root_fsolver_name (const gsl_root_fsolver * @var{s}) |
| @deftypefunx {const char *} gsl_root_fdfsolver_name (const gsl_root_fdfsolver * @var{s}) |
| These functions return a pointer to the name of the solver. For example, |
| |
| @example |
| printf ("s is a '%s' solver\n", |
| gsl_root_fsolver_name (s)); |
| @end example |
| |
| @noindent |
| would print something like @code{s is a 'bisection' solver}. |
| @end deftypefun |
| |
| @node Providing the function to solve |
| @section Providing the function to solve |
| @cindex root finding, providing a function to solve |
| |
| You must provide a continuous function of one variable for the root |
| finders to operate on, and, sometimes, its first derivative. In order |
| to allow for general parameters the functions are defined by the |
| following data types: |
| |
| @deftp {Data Type} gsl_function |
| This data type defines a general function with parameters. |
| |
| @table @code |
| @item double (* function) (double @var{x}, void * @var{params}) |
| this function should return the value |
| @c{$f(x,\hbox{\it params})$} |
| @math{f(x,params)} for argument @var{x} and parameters @var{params} |
| |
| @item void * params |
| a pointer to the parameters of the function |
| @end table |
| @end deftp |
| |
| Here is an example for the general quadratic function, |
| @tex |
| \beforedisplay |
| $$ |
| f(x) = a x^2 + b x + c |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| f(x) = a x^2 + b x + c |
| @end example |
| |
| @end ifinfo |
| @noindent |
| with @math{a = 3}, @math{b = 2}, @math{c = 1}. The following code |
| defines a @code{gsl_function} @code{F} which you could pass to a root |
| finder: |
| |
| @example |
| struct my_f_params @{ double a; double b; double c; @}; |
| |
| double |
| my_f (double x, void * p) @{ |
| struct my_f_params * params |
| = (struct my_f_params *)p; |
| double a = (params->a); |
| double b = (params->b); |
| double c = (params->c); |
| |
| return (a * x + b) * x + c; |
| @} |
| |
| gsl_function F; |
| struct my_f_params params = @{ 3.0, 2.0, 1.0 @}; |
| |
| F.function = &my_f; |
| F.params = ¶ms; |
| @end example |
| |
| @noindent |
| The function @math{f(x)} can be evaluated using the following macro, |
| |
| @example |
| #define GSL_FN_EVAL(F,x) |
| (*((F)->function))(x,(F)->params) |
| @end example |
| |
| @deftp {Data Type} gsl_function_fdf |
| This data type defines a general function with parameters and its first |
| derivative. |
| |
| @table @code |
| @item double (* f) (double @var{x}, void * @var{params}) |
| this function should return the value of |
| @c{$f(x,\hbox{\it params})$} |
| @math{f(x,params)} for argument @var{x} and parameters @var{params} |
| |
| @item double (* df) (double @var{x}, void * @var{params}) |
| this function should return the value of the derivative of @var{f} with |
| respect to @var{x}, |
| @c{$f'(x,\hbox{\it params})$} |
| @math{f'(x,params)}, for argument @var{x} and parameters @var{params} |
| |
| @item void (* fdf) (double @var{x}, void * @var{params}, double * @var{f}, double * @var{d}f) |
| this function should set the values of the function @var{f} to |
| @c{$f(x,\hbox{\it params})$} |
| @math{f(x,params)} |
| and its derivative @var{df} to |
| @c{$f'(x,\hbox{\it params})$} |
| @math{f'(x,params)} |
| for argument @var{x} and parameters @var{params}. This function |
| provides an optimization of the separate functions for @math{f(x)} and |
| @math{f'(x)}---it is always faster to compute the function and its |
| derivative at the same time. |
| |
| @item void * params |
| a pointer to the parameters of the function |
| @end table |
| @end deftp |
| |
| Here is an example where |
| @c{$f(x) = \exp(2x)$} |
| @math{f(x) = 2\exp(2x)}: |
| |
| @example |
| double |
| my_f (double x, void * params) |
| @{ |
| return exp (2 * x); |
| @} |
| |
| double |
| my_df (double x, void * params) |
| @{ |
| return 2 * exp (2 * x); |
| @} |
| |
| void |
| my_fdf (double x, void * params, |
| double * f, double * df) |
| @{ |
| double t = exp (2 * x); |
| |
| *f = t; |
| *df = 2 * t; /* uses existing value */ |
| @} |
| |
| gsl_function_fdf FDF; |
| |
| FDF.f = &my_f; |
| FDF.df = &my_df; |
| FDF.fdf = &my_fdf; |
| FDF.params = 0; |
| @end example |
| |
| @noindent |
| The function @math{f(x)} can be evaluated using the following macro, |
| |
| @example |
| #define GSL_FN_FDF_EVAL_F(FDF,x) |
| (*((FDF)->f))(x,(FDF)->params) |
| @end example |
| |
| @noindent |
| The derivative @math{f'(x)} can be evaluated using the following macro, |
| |
| @example |
| #define GSL_FN_FDF_EVAL_DF(FDF,x) |
| (*((FDF)->df))(x,(FDF)->params) |
| @end example |
| |
| @noindent |
| and both the function @math{y = f(x)} and its derivative @math{dy = f'(x)} |
| can be evaluated at the same time using the following macro, |
| |
| @example |
| #define GSL_FN_FDF_EVAL_F_DF(FDF,x,y,dy) |
| (*((FDF)->fdf))(x,(FDF)->params,(y),(dy)) |
| @end example |
| |
| @noindent |
| The macro stores @math{f(x)} in its @var{y} argument and @math{f'(x)} in |
| its @var{dy} argument---both of these should be pointers to |
| @code{double}. |
| |
| @node Search Bounds and Guesses |
| @section Search Bounds and Guesses |
| @cindex root finding, search bounds |
| @cindex root finding, initial guess |
| |
| You provide either search bounds or an initial guess; this section |
| explains how search bounds and guesses work and how function arguments |
| control them. |
| |
| A guess is simply an @math{x} value which is iterated until it is within |
| the desired precision of a root. It takes the form of a @code{double}. |
| |
| Search bounds are the endpoints of a interval which is iterated until |
| the length of the interval is smaller than the requested precision. The |
| interval is defined by two values, the lower limit and the upper limit. |
| Whether the endpoints are intended to be included in the interval or not |
| depends on the context in which the interval is used. |
| |
| @node Root Finding Iteration |
| @section Iteration |
| |
| The following functions drive the iteration of each algorithm. Each |
| function performs one iteration to update the state of any solver of the |
| corresponding type. The same functions work for all solvers so that |
| different methods can be substituted at runtime without modifications to |
| the code. |
| |
| @deftypefun int gsl_root_fsolver_iterate (gsl_root_fsolver * @var{s}) |
| @deftypefunx int gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * @var{s}) |
| These functions perform a single iteration of the solver @var{s}. If the |
| iteration encounters an unexpected problem then an error code will be |
| returned, |
| |
| @table @code |
| @item GSL_EBADFUNC |
| the iteration encountered a singular point where the function or its |
| derivative evaluated to @code{Inf} or @code{NaN}. |
| |
| @item GSL_EZERODIV |
| the derivative of the function vanished at the iteration point, |
| preventing the algorithm from continuing without a division by zero. |
| @end table |
| @end deftypefun |
| |
| The solver maintains a current best estimate of the root at all |
| times. The bracketing solvers also keep track of the current best |
| interval bounding the root. This information can be accessed with the |
| following auxiliary functions, |
| |
| @deftypefun double gsl_root_fsolver_root (const gsl_root_fsolver * @var{s}) |
| @deftypefunx double gsl_root_fdfsolver_root (const gsl_root_fdfsolver * @var{s}) |
| These functions return the current estimate of the root for the solver @var{s}. |
| @end deftypefun |
| |
| @deftypefun double gsl_root_fsolver_x_lower (const gsl_root_fsolver * @var{s}) |
| @deftypefunx double gsl_root_fsolver_x_upper (const gsl_root_fsolver * @var{s}) |
| These functions return the current bracketing interval for the solver @var{s}. |
| @end deftypefun |
| |
| @node Search Stopping Parameters |
| @section Search Stopping Parameters |
| @cindex root finding, stopping parameters |
| |
| A root finding procedure should stop when one of the following conditions is |
| true: |
| |
| @itemize @bullet |
| @item |
| A root has been found to within the user-specified precision. |
| |
| @item |
| A user-specified maximum number of iterations has been reached. |
| |
| @item |
| An error has occurred. |
| @end itemize |
| |
| @noindent |
| The handling of these conditions is under user control. The functions |
| below allow the user to test the precision of the current result in |
| several standard ways. |
| |
| @deftypefun int gsl_root_test_interval (double @var{x_lower}, double @var{x_upper}, double @var{epsabs}, double @var{epsrel}) |
| This function tests for the convergence of the interval [@var{x_lower}, |
| @var{x_upper}] with absolute error @var{epsabs} and relative error |
| @var{epsrel}. The test returns @code{GSL_SUCCESS} if the following |
| condition is achieved, |
| @tex |
| \beforedisplay |
| $$ |
| |a - b| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, \min(|a|,|b|) |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| |a - b| < epsabs + epsrel min(|a|,|b|) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| when the interval @math{x = [a,b]} does not include the origin. If the |
| interval includes the origin then @math{\min(|a|,|b|)} is replaced by |
| zero (which is the minimum value of @math{|x|} over the interval). This |
| ensures that the relative error is accurately estimated for roots close |
| to the origin. |
| |
| This condition on the interval also implies that any estimate of the |
| root @math{r} in the interval satisfies the same condition with respect |
| to the true root @math{r^*}, |
| @tex |
| \beforedisplay |
| $$ |
| |r - r^*| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, r^* |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| |r - r^*| < epsabs + epsrel r^* |
| @end example |
| |
| @end ifinfo |
| @noindent |
| assuming that the true root @math{r^*} is contained within the interval. |
| @end deftypefun |
| |
| @deftypefun int gsl_root_test_delta (double @var{x1}, double @var{x0}, double @var{epsabs}, double @var{epsrel}) |
| |
| This function tests for the convergence of the sequence @dots{}, @var{x0}, |
| @var{x1} with absolute error @var{epsabs} and relative error |
| @var{epsrel}. The test returns @code{GSL_SUCCESS} if the following |
| condition is achieved, |
| @tex |
| \beforedisplay |
| $$ |
| |x_1 - x_0| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, |x_1| |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| |x_1 - x_0| < epsabs + epsrel |x_1| |
| @end example |
| |
| @end ifinfo |
| @noindent |
| and returns @code{GSL_CONTINUE} otherwise. |
| @end deftypefun |
| |
| |
| @deftypefun int gsl_root_test_residual (double @var{f}, double @var{epsabs}) |
| This function tests the residual value @var{f} against the absolute |
| error bound @var{epsabs}. The test returns @code{GSL_SUCCESS} if the |
| following condition is achieved, |
| @tex |
| \beforedisplay |
| $$ |
| |f| < \hbox{\it epsabs} |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| |f| < epsabs |
| @end example |
| |
| @end ifinfo |
| @noindent |
| and returns @code{GSL_CONTINUE} otherwise. This criterion is suitable |
| for situations where the precise location of the root, @math{x}, is |
| unimportant provided a value can be found where the residual, |
| @math{|f(x)|}, is small enough. |
| @end deftypefun |
| |
| @comment ============================================================ |
| |
| @node Root Bracketing Algorithms |
| @section Root Bracketing Algorithms |
| |
| The root bracketing algorithms described in this section require an |
| initial interval which is guaranteed to contain a root---if @math{a} |
| and @math{b} are the endpoints of the interval then @math{f(a)} must |
| differ in sign from @math{f(b)}. This ensures that the function crosses |
| zero at least once in the interval. If a valid initial interval is used |
| then these algorithm cannot fail, provided the function is well-behaved. |
| |
| Note that a bracketing algorithm cannot find roots of even degree, since |
| these do not cross the @math{x}-axis. |
| |
| @deffn {Solver} gsl_root_fsolver_bisection |
| |
| @cindex bisection algorithm for finding roots |
| @cindex root finding, bisection algorithm |
| |
| The @dfn{bisection algorithm} is the simplest method of bracketing the |
| roots of a function. It is the slowest algorithm provided by |
| the library, with linear convergence. |
| |
| On each iteration, the interval is bisected and the value of the |
| function at the midpoint is calculated. The sign of this value is used |
| to determine which half of the interval does not contain a root. That |
| half is discarded to give a new, smaller interval containing the |
| root. This procedure can be continued indefinitely until the interval is |
| sufficiently small. |
| |
| At any time the current estimate of the root is taken as the midpoint of |
| the interval. |
| |
| @comment eps file "roots-bisection.eps" |
| @comment @iftex |
| @comment @sp 1 |
| @comment @center @image{roots-bisection,3.4in} |
| |
| @comment @quotation |
| @comment Four iterations of bisection, where @math{a_n} is @math{n}th position of |
| @comment the beginning of the interval and @math{b_n} is the @math{n}th position |
| @comment of the end. The midpoint of each interval is also indicated. |
| @comment @end quotation |
| @comment @end iftex |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Solver} gsl_root_fsolver_falsepos |
| @cindex false position algorithm for finding roots |
| @cindex root finding, false position algorithm |
| |
| The @dfn{false position algorithm} is a method of finding roots based on |
| linear interpolation. Its convergence is linear, but it is usually |
| faster than bisection. |
| |
| On each iteration a line is drawn between the endpoints @math{(a,f(a))} |
| and @math{(b,f(b))} and the point where this line crosses the |
| @math{x}-axis taken as a ``midpoint''. The value of the function at |
| this point is calculated and its sign is used to determine which side of |
| the interval does not contain a root. That side is discarded to give a |
| new, smaller interval containing the root. This procedure can be |
| continued indefinitely until the interval is sufficiently small. |
| |
| The best estimate of the root is taken from the linear interpolation of |
| the interval on the current iteration. |
| |
| @comment eps file "roots-false-position.eps" |
| @comment @iftex |
| @comment @image{roots-false-position,4in} |
| @comment @quotation |
| @comment Several iterations of false position, where @math{a_n} is @math{n}th |
| @comment position of the beginning of the interval and @math{b_n} is the |
| @comment @math{n}th position of the end. |
| @comment @end quotation |
| @comment @end iftex |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Solver} gsl_root_fsolver_brent |
| @cindex Brent's method for finding roots |
| @cindex root finding, Brent's method |
| |
| The @dfn{Brent-Dekker method} (referred to here as @dfn{Brent's method}) |
| combines an interpolation strategy with the bisection algorithm. This |
| produces a fast algorithm which is still robust. |
| |
| On each iteration Brent's method approximates the function using an |
| interpolating curve. On the first iteration this is a linear |
| interpolation of the two endpoints. For subsequent iterations the |
| algorithm uses an inverse quadratic fit to the last three points, for |
| higher accuracy. The intercept of the interpolating curve with the |
| @math{x}-axis is taken as a guess for the root. If it lies within the |
| bounds of the current interval then the interpolating point is accepted, |
| and used to generate a smaller interval. If the interpolating point is |
| not accepted then the algorithm falls back to an ordinary bisection |
| step. |
| |
| The best estimate of the root is taken from the most recent |
| interpolation or bisection. |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @node Root Finding Algorithms using Derivatives |
| @section Root Finding Algorithms using Derivatives |
| |
| The root polishing algorithms described in this section require an |
| initial guess for the location of the root. There is no absolute |
| guarantee of convergence---the function must be suitable for this |
| technique and the initial guess must be sufficiently close to the root |
| for it to work. When these conditions are satisfied then convergence is |
| quadratic. |
| |
| These algorithms make use of both the function and its derivative. |
| |
| @deffn {Derivative Solver} gsl_root_fdfsolver_newton |
| @cindex Newton's method for finding roots |
| @cindex root finding, Newton's method |
| |
| Newton's Method is the standard root-polishing algorithm. The algorithm |
| begins with an initial guess for the location of the root. On each |
| iteration, a line tangent to the function @math{f} is drawn at that |
| position. The point where this line crosses the @math{x}-axis becomes |
| the new guess. The iteration is defined by the following sequence, |
| @tex |
| \beforedisplay |
| $$ |
| x_{i+1} = x_i - {f(x_i) \over f'(x_i)} |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| x_@{i+1@} = x_i - f(x_i)/f'(x_i) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| Newton's method converges quadratically for single roots, and linearly |
| for multiple roots. |
| |
| @comment eps file "roots-newtons-method.eps" |
| @comment @iftex |
| @comment @sp 1 |
| @comment @center @image{roots-newtons-method,3.4in} |
| |
| @comment @quotation |
| @comment Several iterations of Newton's Method, where @math{g_n} is the |
| @comment @math{n}th guess. |
| @comment @end quotation |
| @comment @end iftex |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Derivative Solver} gsl_root_fdfsolver_secant |
| @cindex secant method for finding roots |
| @cindex root finding, secant method |
| |
| The @dfn{secant method} is a simplified version of Newton's method which does |
| not require the computation of the derivative on every step. |
| |
| On its first iteration the algorithm begins with Newton's method, using |
| the derivative to compute a first step, |
| @tex |
| \beforedisplay |
| $$ |
| x_1 = x_0 - {f(x_0) \over f'(x_0)} |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| x_1 = x_0 - f(x_0)/f'(x_0) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| Subsequent iterations avoid the evaluation of the derivative by |
| replacing it with a numerical estimate, the slope of the line through |
| the previous two points, |
| @tex |
| \beforedisplay |
| $$ |
| x_{i+1} = x_i - {f(x_i) \over f'_{est}} |
| ~\hbox{where}~ |
| f'_{est} = {f(x_{i}) - f(x_{i-1}) \over x_i - x_{i-1}} |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| x_@{i+1@} = x_i f(x_i) / f'_@{est@} where |
| f'_@{est@} = (f(x_i) - f(x_@{i-1@})/(x_i - x_@{i-1@}) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| When the derivative does not change significantly in the vicinity of the |
| root the secant method gives a useful saving. Asymptotically the secant |
| method is faster than Newton's method whenever the cost of evaluating |
| the derivative is more than 0.44 times the cost of evaluating the |
| function itself. As with all methods of computing a numerical |
| derivative the estimate can suffer from cancellation errors if the |
| separation of the points becomes too small. |
| |
| On single roots, the method has a convergence of order @math{(1 + \sqrt |
| 5)/2} (approximately @math{1.62}). It converges linearly for multiple |
| roots. |
| |
| @comment eps file "roots-secant-method.eps" |
| @comment @iftex |
| @comment @tex |
| @comment \input epsf |
| @comment \medskip |
| @comment \centerline{\epsfxsize=5in\epsfbox{roots-secant-method.eps}} |
| @comment @end tex |
| @comment @quotation |
| @comment Several iterations of Secant Method, where @math{g_n} is the @math{n}th |
| @comment guess. |
| @comment @end quotation |
| @comment @end iftex |
| @end deffn |
| |
| @comment ============================================================ |
| |
| @deffn {Derivative Solver} gsl_root_fdfsolver_steffenson |
| @cindex Steffenson's method for finding roots |
| @cindex root finding, Steffenson's method |
| |
| The @dfn{Steffenson Method} provides the fastest convergence of all the |
| routines. It combines the basic Newton algorithm with an Aitken |
| ``delta-squared'' acceleration. If the Newton iterates are @math{x_i} |
| then the acceleration procedure generates a new sequence @math{R_i}, |
| @tex |
| \beforedisplay |
| $$ |
| R_i = x_i - {(x_{i+1} - x_i)^2 \over (x_{i+2} - 2 x_{i+1} + x_i)} |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| R_i = x_i - (x_@{i+1@} - x_i)^2 / (x_@{i+2@} - 2 x_@{i+1@} + x_@{i@}) |
| @end example |
| |
| @end ifinfo |
| @noindent |
| which converges faster than the original sequence under reasonable |
| conditions. The new sequence requires three terms before it can produce |
| its first value so the method returns accelerated values on the second |
| and subsequent iterations. On the first iteration it returns the |
| ordinary Newton estimate. The Newton iterate is also returned if the |
| denominator of the acceleration term ever becomes zero. |
| |
| As with all acceleration procedures this method can become unstable if |
| the function is not well-behaved. |
| @end deffn |
| |
| @node Root Finding Examples |
| @section Examples |
| |
| For any root finding algorithm we need to prepare the function to be |
| solved. For this example we will use the general quadratic equation |
| described earlier. We first need a header file (@file{demo_fn.h}) to |
| define the function parameters, |
| |
| @example |
| @verbatiminclude examples/demo_fn.h |
| @end example |
| |
| @noindent |
| We place the function definitions in a separate file (@file{demo_fn.c}), |
| |
| @example |
| @verbatiminclude examples/demo_fn.c |
| @end example |
| |
| @noindent |
| The first program uses the function solver @code{gsl_root_fsolver_brent} |
| for Brent's method and the general quadratic defined above to solve the |
| following equation, |
| @tex |
| \beforedisplay |
| $$ |
| x^2 - 5 = 0 |
| $$ |
| \afterdisplay |
| @end tex |
| @ifinfo |
| |
| @example |
| x^2 - 5 = 0 |
| @end example |
| |
| @end ifinfo |
| @noindent |
| with solution @math{x = \sqrt 5 = 2.236068...} |
| |
| @example |
| @verbatiminclude examples/roots.c |
| @end example |
| |
| @noindent |
| Here are the results of the iterations, |
| |
| @smallexample |
| $ ./a.out |
| using brent method |
| iter [ lower, upper] root err err(est) |
| 1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000 |
| 2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000 |
| 3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000 |
| 4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000 |
| 5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300 |
| Converged: |
| 6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666 |
| @end smallexample |
| |
| @noindent |
| If the program is modified to use the bisection solver instead of |
| Brent's method, by changing @code{gsl_root_fsolver_brent} to |
| @code{gsl_root_fsolver_bisection} the slower convergence of the |
| Bisection method can be observed, |
| |
| @smallexample |
| $ ./a.out |
| using bisection method |
| iter [ lower, upper] root err err(est) |
| 1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000 |
| 2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000 |
| 3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000 |
| 4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000 |
| 5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500 |
| 6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250 |
| 7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625 |
| 8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312 |
| 9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656 |
| 10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828 |
| 11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414 |
| Converged: |
| 12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207 |
| @end smallexample |
| |
| The next program solves the same function using a derivative solver |
| instead. |
| |
| @example |
| @verbatiminclude examples/rootnewt.c |
| @end example |
| |
| @noindent |
| Here are the results for Newton's method, |
| |
| @example |
| $ ./a.out |
| using newton method |
| iter root err err(est) |
| 1 3.0000000 +0.7639320 -2.0000000 |
| 2 2.3333333 +0.0972654 -0.6666667 |
| 3 2.2380952 +0.0020273 -0.0952381 |
| Converged: |
| 4 2.2360689 +0.0000009 -0.0020263 |
| @end example |
| |
| @noindent |
| Note that the error can be estimated more accurately by taking the |
| difference between the current iterate and next iterate rather than the |
| previous iterate. The other derivative solvers can be investigated by |
| changing @code{gsl_root_fdfsolver_newton} to |
| @code{gsl_root_fdfsolver_secant} or |
| @code{gsl_root_fdfsolver_steffenson}. |
| |
| @node Root Finding References and Further Reading |
| @section References and Further Reading |
| |
| For information on the Brent-Dekker algorithm see the following two |
| papers, |
| |
| @itemize @asis |
| @item |
| R. P. Brent, ``An algorithm with guaranteed convergence for finding a |
| zero of a function'', @cite{Computer Journal}, 14 (1971) 422--425 |
| |
| @item |
| J. C. P. Bus and T. J. Dekker, ``Two Efficient Algorithms with Guaranteed |
| Convergence for Finding a Zero of a Function'', @cite{ACM Transactions of |
| Mathematical Software}, Vol.@: 1 No.@: 4 (1975) 330--345 |
| @end itemize |
| |