| @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 |
| |
| |
| |
| |
| |
| |