diff options
Diffstat (limited to 'gsl-1.9/doc/roots.texi')
| -rw-r--r-- | gsl-1.9/doc/roots.texi | 928 |
1 files changed, 928 insertions, 0 deletions
diff --git a/gsl-1.9/doc/roots.texi b/gsl-1.9/doc/roots.texi new file mode 100644 index 0000000..fdcefe1 --- /dev/null +++ b/gsl-1.9/doc/roots.texi @@ -0,0 +1,928 @@ +@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 + |