diff options
Diffstat (limited to 'gsl-1.9/doc/integration.texi')
| -rw-r--r-- | gsl-1.9/doc/integration.texi | 854 |
1 files changed, 854 insertions, 0 deletions
diff --git a/gsl-1.9/doc/integration.texi b/gsl-1.9/doc/integration.texi new file mode 100644 index 0000000..62afb3e --- /dev/null +++ b/gsl-1.9/doc/integration.texi @@ -0,0 +1,854 @@ +@cindex quadrature +@cindex numerical integration (quadrature) +@cindex integration, numerical (quadrature) +@cindex QUADPACK + +This chapter describes routines for performing numerical integration +(quadrature) of a function in one dimension. There are routines for +adaptive and non-adaptive integration of general functions, with +specialised routines for specific cases. These include integration over +infinite and semi-infinite ranges, singular integrals, including +logarithmic singularities, computation of Cauchy principal values and +oscillatory integrals. The library reimplements the algorithms used in +@sc{quadpack}, a numerical integration package written by Piessens, +Doncker-Kapenga, Uberhuber and Kahaner. Fortran code for @sc{quadpack} is +available on Netlib. + +The functions described in this chapter are declared in the header file +@file{gsl_integration.h}. + +@menu +* Numerical Integration Introduction:: +* QNG non-adaptive Gauss-Kronrod integration:: +* QAG adaptive integration:: +* QAGS adaptive integration with singularities:: +* QAGP adaptive integration with known singular points:: +* QAGI adaptive integration on infinite intervals:: +* QAWC adaptive integration for Cauchy principal values:: +* QAWS adaptive integration for singular functions:: +* QAWO adaptive integration for oscillatory functions:: +* QAWF adaptive integration for Fourier integrals:: +* Numerical integration error codes:: +* Numerical integration examples:: +* Numerical integration References and Further Reading:: +@end menu + +@node Numerical Integration Introduction +@section Introduction + +Each algorithm computes an approximation to a definite integral of the +form, +@tex +\beforedisplay +$$ +I = \int_a^b f(x) w(x) \,dx +$$ +\afterdisplay +@end tex +@ifinfo + +@example +I = \int_a^b f(x) w(x) dx +@end example + +@end ifinfo +@noindent +where @math{w(x)} is a weight function (for general integrands @math{w(x)=1}). +The user provides absolute and relative error bounds +@c{$(\hbox{\it epsabs}, \hbox{\it epsrel}\,)$} +@math{(epsabs, epsrel)} which specify the following accuracy requirement, +@tex +\beforedisplay +$$ +|\hbox{\it RESULT} - I| \leq \max(\hbox{\it epsabs}, \hbox{\it epsrel}\, |I|) +$$ +\afterdisplay +@end tex +@ifinfo + +@example +|RESULT - I| <= max(epsabs, epsrel |I|) +@end example + +@end ifinfo +@noindent +where +@c{$\hbox{\it RESULT}$} +@math{RESULT} is the numerical approximation obtained by the +algorithm. The algorithms attempt to estimate the absolute error +@c{$\hbox{\it ABSERR} = |\hbox{\it RESULT} - I|$} +@math{ABSERR = |RESULT - I|} in such a way that the following inequality +holds, +@tex +\beforedisplay +$$ +|\hbox{\it RESULT} - I| \leq \hbox{\it ABSERR} \leq \max(\hbox{\it epsabs}, \hbox{\it epsrel}\,|I|) +$$ +\afterdisplay +@end tex +@ifinfo + +@example +|RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|) +@end example + +@end ifinfo +@noindent +The routines will fail to converge if the error bounds are too +stringent, but always return the best approximation obtained up to that +stage. + +The algorithms in @sc{quadpack} use a naming convention based on the +following letters, + +@display +@code{Q} - quadrature routine + +@code{N} - non-adaptive integrator +@code{A} - adaptive integrator + +@code{G} - general integrand (user-defined) +@code{W} - weight function with integrand + +@code{S} - singularities can be more readily integrated +@code{P} - points of special difficulty can be supplied +@code{I} - infinite range of integration +@code{O} - oscillatory weight function, cos or sin +@code{F} - Fourier integral +@code{C} - Cauchy principal value +@end display + +@noindent +The algorithms are built on pairs of quadrature rules, a higher order +rule and a lower order rule. The higher order rule is used to compute +the best approximation to an integral over a small range. The +difference between the results of the higher order rule and the lower +order rule gives an estimate of the error in the approximation. + +@menu +* Integrands without weight functions:: +* Integrands with weight functions:: +* Integrands with singular weight functions:: +@end menu + +@node Integrands without weight functions +@subsection Integrands without weight functions +@cindex Gauss-Kronrod quadrature +The algorithms for general functions (without a weight function) are +based on Gauss-Kronrod rules. + +A Gauss-Kronrod rule begins with a classical Gaussian quadrature rule of +order @math{m}. This is extended with additional points between each of +the abscissae to give a higher order Kronrod rule of order @math{2m+1}. +The Kronrod rule is efficient because it reuses existing function +evaluations from the Gaussian rule. + +The higher order Kronrod rule is used as the best approximation to the +integral, and the difference between the two rules is used as an +estimate of the error in the approximation. + +@node Integrands with weight functions +@subsection Integrands with weight functions +@cindex Clenshaw-Curtis quadrature +@cindex Modified Clenshaw-Curtis quadrature +For integrands with weight functions the algorithms use Clenshaw-Curtis +quadrature rules. + +A Clenshaw-Curtis rule begins with an @math{n}-th order Chebyshev +polynomial approximation to the integrand. This polynomial can be +integrated exactly to give an approximation to the integral of the +original function. The Chebyshev expansion can be extended to higher +orders to improve the approximation and provide an estimate of the +error. + +@node Integrands with singular weight functions +@subsection Integrands with singular weight functions + +The presence of singularities (or other behavior) in the integrand can +cause slow convergence in the Chebyshev approximation. The modified +Clenshaw-Curtis rules used in @sc{quadpack} separate out several common +weight functions which cause slow convergence. + +These weight functions are integrated analytically against the Chebyshev +polynomials to precompute @dfn{modified Chebyshev moments}. Combining +the moments with the Chebyshev approximation to the function gives the +desired integral. The use of analytic integration for the singular part +of the function allows exact cancellations and substantially improves +the overall convergence behavior of the integration. + + +@node QNG non-adaptive Gauss-Kronrod integration +@section QNG non-adaptive Gauss-Kronrod integration +@cindex QNG quadrature algorithm + +The QNG algorithm is a non-adaptive procedure which uses fixed +Gauss-Kronrod abscissae to sample the integrand at a maximum of 87 +points. It is provided for fast integration of smooth functions. + +@deftypefun int gsl_integration_qng (const gsl_function * @var{f}, double @var{a}, double @var{b}, double @var{epsabs}, double @var{epsrel}, double * @var{result}, double * @var{abserr}, size_t * @var{neval}) + +This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and +87-point integration rules in succession until an estimate of the +integral of @math{f} over @math{(a,b)} is achieved within the desired +absolute and relative error limits, @var{epsabs} and @var{epsrel}. The +function returns the final approximation, @var{result}, an estimate of +the absolute error, @var{abserr} and the number of function evaluations +used, @var{neval}. The Gauss-Kronrod rules are designed in such a way +that each rule uses all the results of its predecessors, in order to +minimize the total number of function evaluations. +@end deftypefun + + +@node QAG adaptive integration +@section QAG adaptive integration +@cindex QAG quadrature algorithm + +The QAG algorithm is a simple adaptive integration procedure. The +integration region is divided into subintervals, and on each iteration +the subinterval with the largest estimated error is bisected. This +reduces the overall error rapidly, as the subintervals become +concentrated around local difficulties in the integrand. These +subintervals are managed by a @code{gsl_integration_workspace} struct, +which handles the memory for the subinterval ranges, results and error +estimates. + +@deftypefun {gsl_integration_workspace *} gsl_integration_workspace_alloc (size_t @var{n}) +This function allocates a workspace sufficient to hold @var{n} double +precision intervals, their integration results and error estimates. +@end deftypefun + +@deftypefun void gsl_integration_workspace_free (gsl_integration_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +@end deftypefun + +@deftypefun int gsl_integration_qag (const gsl_function * @var{f}, double @var{a}, double @var{b}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, int @var{key}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function applies an integration rule adaptively until an estimate +of the integral of @math{f} over @math{(a,b)} is achieved within the +desired absolute and relative error limits, @var{epsabs} and +@var{epsrel}. The function returns the final approximation, +@var{result}, and an estimate of the absolute error, @var{abserr}. The +integration rule is determined by the value of @var{key}, which should +be chosen from the following symbolic names, + +@example +GSL_INTEG_GAUSS15 (key = 1) +GSL_INTEG_GAUSS21 (key = 2) +GSL_INTEG_GAUSS31 (key = 3) +GSL_INTEG_GAUSS41 (key = 4) +GSL_INTEG_GAUSS51 (key = 5) +GSL_INTEG_GAUSS61 (key = 6) +@end example + +@noindent +corresponding to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod +rules. The higher-order rules give better accuracy for smooth functions, +while lower-order rules save time when the function contains local +difficulties, such as discontinuities. + +On each iteration the adaptive integration strategy bisects the interval +with the largest error estimate. The subintervals and their results are +stored in the memory provided by @var{workspace}. The maximum number of +subintervals is given by @var{limit}, which may not exceed the allocated +size of the workspace. +@end deftypefun + + +@node QAGS adaptive integration with singularities +@section QAGS adaptive integration with singularities +@cindex QAGS quadrature algorithm + +The presence of an integrable singularity in the integration region +causes an adaptive routine to concentrate new subintervals around the +singularity. As the subintervals decrease in size the successive +approximations to the integral converge in a limiting fashion. This +approach to the limit can be accelerated using an extrapolation +procedure. The QAGS algorithm combines adaptive bisection with the Wynn +epsilon-algorithm to speed up the integration of many types of +integrable singularities. + +@deftypefun int gsl_integration_qags (const gsl_function * @var{f}, double @var{a}, double @var{b}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function applies the Gauss-Kronrod 21-point integration rule +adaptively until an estimate of the integral of @math{f} over +@math{(a,b)} is achieved within the desired absolute and relative error +limits, @var{epsabs} and @var{epsrel}. The results are extrapolated +using the epsilon-algorithm, which accelerates the convergence of the +integral in the presence of discontinuities and integrable +singularities. The function returns the final approximation from the +extrapolation, @var{result}, and an estimate of the absolute error, +@var{abserr}. The subintervals and their results are stored in the +memory provided by @var{workspace}. The maximum number of subintervals +is given by @var{limit}, which may not exceed the allocated size of the +workspace. + +@end deftypefun + +@node QAGP adaptive integration with known singular points +@section QAGP adaptive integration with known singular points +@cindex QAGP quadrature algorithm +@cindex singular points, specifying positions in quadrature +@deftypefun int gsl_integration_qagp (const gsl_function * @var{f}, double * @var{pts}, size_t @var{npts}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function applies the adaptive integration algorithm QAGS taking +account of the user-supplied locations of singular points. The array +@var{pts} of length @var{npts} should contain the endpoints of the +integration ranges defined by the integration region and locations of +the singularities. For example, to integrate over the region +@math{(a,b)} with break-points at @math{x_1, x_2, x_3} (where +@math{a < x_1 < x_2 < x_3 < b}) the following @var{pts} array should be used + +@example +pts[0] = a +pts[1] = x_1 +pts[2] = x_2 +pts[3] = x_3 +pts[4] = b +@end example + +@noindent +with @var{npts} = 5. + +@noindent +If you know the locations of the singular points in the integration +region then this routine will be faster than @code{QAGS}. + +@end deftypefun + +@node QAGI adaptive integration on infinite intervals +@section QAGI adaptive integration on infinite intervals +@cindex QAGI quadrature algorithm + +@deftypefun int gsl_integration_qagi (gsl_function * @var{f}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function computes the integral of the function @var{f} over the +infinite interval @math{(-\infty,+\infty)}. The integral is mapped onto the +semi-open interval @math{(0,1]} using the transformation @math{x = (1-t)/t}, +@tex +\beforedisplay +$$ +\int_{-\infty}^{+\infty} dx \, f(x) + = \int_0^1 dt \, (f((1-t)/t) + f(-(1-t)/t))/t^2. +$$ +\afterdisplay +@end tex +@ifinfo + +@example +\int_@{-\infty@}^@{+\infty@} dx f(x) = + \int_0^1 dt (f((1-t)/t) + f((-1+t)/t))/t^2. +@end example + +@end ifinfo +@noindent +It is then integrated using the QAGS algorithm. The normal 21-point +Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, because the +transformation can generate an integrable singularity at the origin. In +this case a lower-order rule is more efficient. +@end deftypefun + +@deftypefun int gsl_integration_qagiu (gsl_function * @var{f}, double @var{a}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function computes the integral of the function @var{f} over the +semi-infinite interval @math{(a,+\infty)}. The integral is mapped onto the +semi-open interval @math{(0,1]} using the transformation @math{x = a + (1-t)/t}, +@tex +\beforedisplay +$$ +\int_{a}^{+\infty} dx \, f(x) + = \int_0^1 dt \, f(a + (1-t)/t)/t^2 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +\int_@{a@}^@{+\infty@} dx f(x) = + \int_0^1 dt f(a + (1-t)/t)/t^2 +@end example + +@end ifinfo +@noindent +and then integrated using the QAGS algorithm. +@end deftypefun + +@deftypefun int gsl_integration_qagil (gsl_function * @var{f}, double @var{b}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) +This function computes the integral of the function @var{f} over the +semi-infinite interval @math{(-\infty,b)}. The integral is mapped onto the +semi-open interval @math{(0,1]} using the transformation @math{x = b - (1-t)/t}, +@tex +\beforedisplay +$$ +\int_{-\infty}^{b} dx \, f(x) + = \int_0^1 dt \, f(b - (1-t)/t)/t^2 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +\int_@{+\infty@}^@{b@} dx f(x) = + \int_0^1 dt f(b - (1-t)/t)/t^2 +@end example + +@end ifinfo +@noindent +and then integrated using the QAGS algorithm. +@end deftypefun + +@node QAWC adaptive integration for Cauchy principal values +@section QAWC adaptive integration for Cauchy principal values +@cindex QAWC quadrature algorithm +@cindex Cauchy principal value, by numerical quadrature +@deftypefun int gsl_integration_qawc (gsl_function * @var{f}, double @var{a}, double @var{b}, double @var{c}, double @var{epsabs}, double @var{epsrel}, size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function computes the Cauchy principal value of the integral of +@math{f} over @math{(a,b)}, with a singularity at @var{c}, +@tex +\beforedisplay +$$ +I = \int_a^b dx\, {f(x) \over x - c} + = \lim_{\epsilon \to 0} +\left\{ +\int_a^{c-\epsilon} dx\, {f(x) \over x - c} ++ +\int_{c+\epsilon}^b dx\, {f(x) \over x - c} +\right\} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +I = \int_a^b dx f(x) / (x - c) +@end example + +@end ifinfo +@noindent +The adaptive bisection algorithm of QAG is used, with modifications to +ensure that subdivisions do not occur at the singular point @math{x = c}. +When a subinterval contains the point @math{x = c} or is close to +it then a special 25-point modified Clenshaw-Curtis rule is used to control +the singularity. Further away from the +singularity the algorithm uses an ordinary 15-point Gauss-Kronrod +integration rule. + +@end deftypefun + +@node QAWS adaptive integration for singular functions +@section QAWS adaptive integration for singular functions +@cindex QAWS quadrature algorithm +@cindex singular functions, numerical integration of +The QAWS algorithm is designed for integrands with algebraic-logarithmic +singularities at the end-points of an integration region. In order to +work efficiently the algorithm requires a precomputed table of +Chebyshev moments. + +@deftypefun {gsl_integration_qaws_table *} gsl_integration_qaws_table_alloc (double @var{alpha}, double @var{beta}, int @var{mu}, int @var{nu}) + +This function allocates space for a @code{gsl_integration_qaws_table} +struct describing a singular weight function +@math{W(x)} with the parameters @math{(\alpha, \beta, \mu, \nu)}, +@tex +\beforedisplay +$$ +W(x) = (x - a)^\alpha (b - x)^\beta \log^\mu (x - a) \log^\nu (b - x) +$$ +\afterdisplay +@end tex +@ifinfo + +@example +W(x) = (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x) +@end example + +@end ifinfo +@noindent +where @math{\alpha > -1}, @math{\beta > -1}, and @math{\mu = 0, 1}, +@math{\nu = 0, 1}. The weight function can take four different forms +depending on the values of @math{\mu} and @math{\nu}, +@tex +\beforedisplay +$$ +\matrix{ +W(x) = (x - a)^\alpha (b - x)^\beta + \hfill~ (\mu = 0, \nu = 0) \cr +W(x) = (x - a)^\alpha (b - x)^\beta \log(x - a) + \hfill~ (\mu = 1, \nu = 0) \cr +W(x) = (x - a)^\alpha (b - x)^\beta \log(b - x) + \hfill~ (\mu = 0, \nu = 1) \cr +W(x) = (x - a)^\alpha (b - x)^\beta \log(x - a) \log(b - x) + \hfill~ (\mu = 1, \nu = 1) +} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +W(x) = (x-a)^alpha (b-x)^beta (mu = 0, nu = 0) +W(x) = (x-a)^alpha (b-x)^beta log(x-a) (mu = 1, nu = 0) +W(x) = (x-a)^alpha (b-x)^beta log(b-x) (mu = 0, nu = 1) +W(x) = (x-a)^alpha (b-x)^beta log(x-a) log(b-x) (mu = 1, nu = 1) +@end example + +@end ifinfo +@noindent +The singular points @math{(a,b)} do not have to be specified until the +integral is computed, where they are the endpoints of the integration +range. + +The function returns a pointer to the newly allocated table +@code{gsl_integration_qaws_table} if no errors were detected, and 0 in +the case of error. +@end deftypefun + +@deftypefun int gsl_integration_qaws_table_set (gsl_integration_qaws_table * @var{t}, double @var{alpha}, double @var{beta}, int @var{mu}, int @var{nu}) +This function modifies the parameters @math{(\alpha, \beta, \mu, \nu)} of +an existing @code{gsl_integration_qaws_table} struct @var{t}. +@end deftypefun + +@deftypefun void gsl_integration_qaws_table_free (gsl_integration_qaws_table * @var{t}) +This function frees all the memory associated with the +@code{gsl_integration_qaws_table} struct @var{t}. +@end deftypefun + +@deftypefun int gsl_integration_qaws (gsl_function * @var{f}, const double @var{a}, const double @var{b}, gsl_integration_qaws_table * @var{t}, const double @var{epsabs}, const double @var{epsrel}, const size_t @var{limit}, gsl_integration_workspace * @var{workspace}, double * @var{result}, double * @var{abserr}) + +This function computes the integral of the function @math{f(x)} over the +interval @math{(a,b)} with the singular weight function +@math{(x-a)^\alpha (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x)}. The parameters +of the weight function @math{(\alpha, \beta, \mu, \nu)} are taken from the +table @var{t}. The integral is, +@tex +\beforedisplay +$$ +I = \int_a^b dx\, f(x) (x - a)^\alpha (b - x)^\beta + \log^\mu (x - a) \log^\nu (b - x). +$$ +\afterdisplay +@end tex +@ifinfo + +@example +I = \int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x). +@end example + +@end ifinfo +@noindent +The adaptive bisection algorithm of QAG is used. When a subinterval +contains one of the endpoints then a special 25-point modified +Clenshaw-Curtis rule is used to control the singularities. For +subintervals which do not include the endpoints an ordinary 15-point +Gauss-Kronrod integration rule is used. + +@end deftypefun + +@node QAWO adaptive integration for oscillatory functions +@section QAWO adaptive integration for oscillatory functions +@cindex QAWO quadrature algorithm +@cindex oscillatory functions, numerical integration of +The QAWO algorithm is designed for integrands with an oscillatory +factor, @math{\sin(\omega x)} or @math{\cos(\omega x)}. In order to +work efficiently the algorithm requires a table of Chebyshev moments +which must be pre-computed with calls to the functions below. + +@deftypefun {gsl_integration_qawo_table *} gsl_integration_qawo_table_alloc (double @var{omega}, double @var{L}, enum gsl_integration_qawo_enum @var{sine}, size_t @var{n}) + +This function allocates space for a @code{gsl_integration_qawo_table} +struct and its associated workspace describing a sine or cosine weight +function @math{W(x)} with the parameters @math{(\omega, L)}, +@tex +\beforedisplay +$$ +\eqalign{ +W(x) & = \left\{\matrix{\sin(\omega x) \cr \cos(\omega x)} \right\} +} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +W(x) = sin(omega x) +W(x) = cos(omega x) +@end example + +@end ifinfo +@noindent +The parameter @var{L} must be the length of the interval over which the +function will be integrated @math{L = b - a}. The choice of sine or +cosine is made with the parameter @var{sine} which should be chosen from +one of the two following symbolic values: + +@example +GSL_INTEG_COSINE +GSL_INTEG_SINE +@end example + +@noindent +The @code{gsl_integration_qawo_table} is a table of the trigonometric +coefficients required in the integration process. The parameter @var{n} +determines the number of levels of coefficients that are computed. Each +level corresponds to one bisection of the interval @math{L}, so that +@var{n} levels are sufficient for subintervals down to the length +@math{L/2^n}. The integration routine @code{gsl_integration_qawo} +returns the error @code{GSL_ETABLE} if the number of levels is +insufficient for the requested accuracy. + +@end deftypefun + +@deftypefun int gsl_integration_qawo_table_set (gsl_integration_qawo_table * @var{t}, double @var{omega}, double @var{L}, enum gsl_integration_qawo_enum @var{sine}) +This function changes the parameters @var{omega}, @var{L} and @var{sine} +of the existing workspace @var{t}. +@end deftypefun + +@deftypefun int gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * @var{t}, double @var{L}) +This function allows the length parameter @var{L} of the workspace +@var{t} to be changed. +@end deftypefun + +@deftypefun void gsl_integration_qawo_table_free (gsl_integration_qawo_table * @var{t}) +This function frees all the memory associated with the workspace @var{t}. +@end deftypefun + +@deftypefun int gsl_integration_qawo (gsl_function * @var{f}, const double @var{a}, const double @var{epsabs}, const double @var{epsrel}, const size_t @var{limit}, gsl_integration_workspace * @var{workspace}, gsl_integration_qawo_table * @var{wf}, double * @var{result}, double * @var{abserr}) + +This function uses an adaptive algorithm to compute the integral of +@math{f} over @math{(a,b)} with the weight function +@math{\sin(\omega x)} or @math{\cos(\omega x)} defined +by the table @var{wf}, +@tex +\beforedisplay +$$ +\eqalign{ +I & = \int_a^b dx\, f(x) \left\{ \matrix{\sin(\omega x) \cr \cos(\omega x)}\right\} +} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +I = \int_a^b dx f(x) sin(omega x) +I = \int_a^b dx f(x) cos(omega x) +@end example + +@end ifinfo +@noindent +The results are extrapolated using the epsilon-algorithm to accelerate +the convergence of the integral. The function returns the final +approximation from the extrapolation, @var{result}, and an estimate of +the absolute error, @var{abserr}. The subintervals and their results are +stored in the memory provided by @var{workspace}. The maximum number of +subintervals is given by @var{limit}, which may not exceed the allocated +size of the workspace. + +Those subintervals with ``large'' widths @math{d} where @math{d\omega > 4} are +computed using a 25-point Clenshaw-Curtis integration rule, which handles the +oscillatory behavior. Subintervals with a ``small'' widths where +@math{d\omega < 4} are computed using a 15-point Gauss-Kronrod integration. + +@end deftypefun + +@node QAWF adaptive integration for Fourier integrals +@section QAWF adaptive integration for Fourier integrals +@cindex QAWF quadrature algorithm +@cindex Fourier integrals, numerical + +@deftypefun int gsl_integration_qawf (gsl_function * @var{f}, const double @var{a}, const double @var{epsabs}, const size_t @var{limit}, gsl_integration_workspace * @var{workspace}, gsl_integration_workspace * @var{cycle_workspace}, gsl_integration_qawo_table * @var{wf}, double * @var{result}, double * @var{abserr}) + +This function attempts to compute a Fourier integral of the function +@var{f} over the semi-infinite interval @math{[a,+\infty)}. +@tex +\beforedisplay +$$ +\eqalign{ +I & = \int_a^{+\infty} dx\, f(x) \left\{ \matrix{ \sin(\omega x) \cr + \cos(\omega x) } \right\} +} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +I = \int_a^@{+\infty@} dx f(x) sin(omega x) +I = \int_a^@{+\infty@} dx f(x) cos(omega x) +@end example +@end ifinfo + +The parameter @math{\omega} and choice of @math{\sin} or @math{\cos} is +taken from the table @var{wf} (the length @var{L} can take any value, +since it is overridden by this function to a value appropriate for the +fourier integration). The integral is computed using the QAWO algorithm +over each of the subintervals, +@tex +\beforedisplay +$$ +\eqalign{ +C_1 & = [a, a + c] \cr +C_2 & = [a + c, a + 2c] \cr +\dots & = \dots \cr +C_k & = [a + (k-1) c, a + k c] +} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +C_1 = [a, a + c] +C_2 = [a + c, a + 2 c] +... = ... +C_k = [a + (k-1) c, a + k c] +@end example + +@end ifinfo +@noindent +where +@c{$c = (2 \,\hbox{floor}(|\omega|) + 1) \pi/|\omega|$} +@math{c = (2 floor(|\omega|) + 1) \pi/|\omega|}. The width @math{c} is +chosen to cover an odd number of periods so that the contributions from +the intervals alternate in sign and are monotonically decreasing when +@var{f} is positive and monotonically decreasing. The sum of this +sequence of contributions is accelerated using the epsilon-algorithm. + +This function works to an overall absolute tolerance of +@var{abserr}. The following strategy is used: on each interval +@math{C_k} the algorithm tries to achieve the tolerance +@tex +\beforedisplay +$$ +TOL_k = u_k \hbox{\it abserr} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +TOL_k = u_k abserr +@end example + +@end ifinfo +@noindent +where +@c{$u_k = (1 - p)p^{k-1}$} +@math{u_k = (1 - p)p^@{k-1@}} and @math{p = 9/10}. +The sum of the geometric series of contributions from each interval +gives an overall tolerance of @var{abserr}. + +If the integration of a subinterval leads to difficulties then the +accuracy requirement for subsequent intervals is relaxed, +@tex +\beforedisplay +$$ +TOL_k = u_k \max(\hbox{\it abserr}, \max_{i<k}\{E_i\}) +$$ +\afterdisplay +@end tex +@ifinfo + +@example +TOL_k = u_k max(abserr, max_@{i<k@}@{E_i@}) +@end example + +@end ifinfo +@noindent +where @math{E_k} is the estimated error on the interval @math{C_k}. + +The subintervals and their results are stored in the memory provided by +@var{workspace}. The maximum number of subintervals is given by +@var{limit}, which may not exceed the allocated size of the workspace. +The integration over each subinterval uses the memory provided by +@var{cycle_workspace} as workspace for the QAWO algorithm. + +@end deftypefun + +@node Numerical integration error codes +@section Error codes + +In addition to the standard error codes for invalid arguments the +functions can return the following values, + +@table @code +@item GSL_EMAXITER +the maximum number of subdivisions was exceeded. +@item GSL_EROUND +cannot reach tolerance because of roundoff error, +or roundoff error was detected in the extrapolation table. +@item GSL_ESING +a non-integrable singularity or other bad integrand behavior was found +in the integration interval. +@item GSL_EDIVERGE +the integral is divergent, or too slowly convergent to be integrated +numerically. +@end table + +@node Numerical integration examples +@section Examples + +The integrator @code{QAGS} will handle a large class of definite +integrals. For example, consider the following integral, which has a +algebraic-logarithmic singularity at the origin, +@tex +\beforedisplay +$$ +\int_0^1 x^{-1/2} \log(x) \,dx = -4 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +\int_0^1 x^@{-1/2@} log(x) dx = -4 +@end example + +@end ifinfo +@noindent +The program below computes this integral to a relative accuracy bound of +@code{1e-7}. + +@example +@verbatiminclude examples/integration.c +@end example + +@noindent +The results below show that the desired accuracy is achieved after 8 +subdivisions. + +@example +$ ./a.out +@verbatiminclude examples/integration.out +@end example + +@noindent +In fact, the extrapolation procedure used by @code{QAGS} produces an +accuracy of almost twice as many digits. The error estimate returned by +the extrapolation procedure is larger than the actual error, giving a +margin of safety of one order of magnitude. + + +@node Numerical integration References and Further Reading +@section References and Further Reading + +The following book is the definitive reference for @sc{quadpack}, and was +written by the original authors. It provides descriptions of the +algorithms, program listings, test programs and examples. It also +includes useful advice on numerical integration and many references to +the numerical integration literature used in developing @sc{quadpack}. + +@itemize @asis +@item +R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. Kahaner. +@cite{@sc{quadpack} A subroutine package for automatic integration} +Springer Verlag, 1983. +@end itemize + +@noindent + + + + + + |