2 @cindex numerical integration (quadrature)
3 @cindex integration, numerical (quadrature)
6 This chapter describes routines for performing numerical integration
7 (quadrature) of a function in one dimension. There are routines for
8 adaptive and non-adaptive integration of general functions, with
9 specialised routines for specific cases. These include integration over
10 infinite and semi-infinite ranges, singular integrals, including
11 logarithmic singularities, computation of Cauchy principal values and
12 oscillatory integrals. The library reimplements the algorithms used in
13 @sc{quadpack}, a numerical integration package written by Piessens,
14 Doncker-Kapenga, Uberhuber and Kahaner. Fortran code for @sc{quadpack} is
17 The functions described in this chapter are declared in the header file
18 @file{gsl_integration.h}.
21 * Numerical Integration Introduction::
22 * QNG non-adaptive Gauss-Kronrod integration::
23 * QAG adaptive integration::
24 * QAGS adaptive integration with singularities::
25 * QAGP adaptive integration with known singular points::
26 * QAGI adaptive integration on infinite intervals::
27 * QAWC adaptive integration for Cauchy principal values::
28 * QAWS adaptive integration for singular functions::
29 * QAWO adaptive integration for oscillatory functions::
30 * QAWF adaptive integration for Fourier integrals::
31 * Numerical integration error codes::
32 * Numerical integration examples::
33 * Numerical integration References and Further Reading::
36 @node Numerical Integration Introduction
39 Each algorithm computes an approximation to a definite integral of the
44 I = \int_a^b f(x) w(x) \,dx
51 I = \int_a^b f(x) w(x) dx
56 where @math{w(x)} is a weight function (for general integrands @math{w(x)=1}).
57 The user provides absolute and relative error bounds
58 @c{$(\hbox{\it epsabs}, \hbox{\it epsrel}\,)$}
59 @math{(epsabs, epsrel)} which specify the following accuracy requirement,
63 |\hbox{\it RESULT} - I| \leq \max(\hbox{\it epsabs}, \hbox{\it epsrel}\, |I|)
70 |RESULT - I| <= max(epsabs, epsrel |I|)
76 @c{$\hbox{\it RESULT}$}
77 @math{RESULT} is the numerical approximation obtained by the
78 algorithm. The algorithms attempt to estimate the absolute error
79 @c{$\hbox{\it ABSERR} = |\hbox{\it RESULT} - I|$}
80 @math{ABSERR = |RESULT - I|} in such a way that the following inequality
85 |\hbox{\it RESULT} - I| \leq \hbox{\it ABSERR} \leq \max(\hbox{\it epsabs}, \hbox{\it epsrel}\,|I|)
92 |RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|)
97 In short, the routines return the first approximation
98 which has an absolute error smaller than @c{$\hbox{\it epsabs}$}
99 @math{epsabs} or a relative error smaller than @c{$\hbox{\it epsrel}$}
102 Note that this is an @i{either-or} constraint,
103 not simultaneous. To compute to a specified absolute error, set @c{$\hbox{\it epsrel}$}
104 @math{epsrel} to zero. To compute to a specified relative error,
105 set @c{$\hbox{\it epsabs}$}
106 @math{epsabs} to zero.
107 The routines will fail to converge if the error bounds are too
108 stringent, but always return the best approximation obtained up to
111 The algorithms in @sc{quadpack} use a naming convention based on the
115 @code{Q} - quadrature routine
117 @code{N} - non-adaptive integrator
118 @code{A} - adaptive integrator
120 @code{G} - general integrand (user-defined)
121 @code{W} - weight function with integrand
123 @code{S} - singularities can be more readily integrated
124 @code{P} - points of special difficulty can be supplied
125 @code{I} - infinite range of integration
126 @code{O} - oscillatory weight function, cos or sin
127 @code{F} - Fourier integral
128 @code{C} - Cauchy principal value
132 The algorithms are built on pairs of quadrature rules, a higher order
133 rule and a lower order rule. The higher order rule is used to compute
134 the best approximation to an integral over a small range. The
135 difference between the results of the higher order rule and the lower
136 order rule gives an estimate of the error in the approximation.
139 * Integrands without weight functions::
140 * Integrands with weight functions::
141 * Integrands with singular weight functions::
144 @node Integrands without weight functions
145 @subsection Integrands without weight functions
146 @cindex Gauss-Kronrod quadrature
147 The algorithms for general functions (without a weight function) are
148 based on Gauss-Kronrod rules.
150 A Gauss-Kronrod rule begins with a classical Gaussian quadrature rule of
151 order @math{m}. This is extended with additional points between each of
152 the abscissae to give a higher order Kronrod rule of order @math{2m+1}.
153 The Kronrod rule is efficient because it reuses existing function
154 evaluations from the Gaussian rule.
156 The higher order Kronrod rule is used as the best approximation to the
157 integral, and the difference between the two rules is used as an
158 estimate of the error in the approximation.
160 @node Integrands with weight functions
161 @subsection Integrands with weight functions
162 @cindex Clenshaw-Curtis quadrature
163 @cindex Modified Clenshaw-Curtis quadrature
164 For integrands with weight functions the algorithms use Clenshaw-Curtis
167 A Clenshaw-Curtis rule begins with an @math{n}-th order Chebyshev
168 polynomial approximation to the integrand. This polynomial can be
169 integrated exactly to give an approximation to the integral of the
170 original function. The Chebyshev expansion can be extended to higher
171 orders to improve the approximation and provide an estimate of the
174 @node Integrands with singular weight functions
175 @subsection Integrands with singular weight functions
177 The presence of singularities (or other behavior) in the integrand can
178 cause slow convergence in the Chebyshev approximation. The modified
179 Clenshaw-Curtis rules used in @sc{quadpack} separate out several common
180 weight functions which cause slow convergence.
182 These weight functions are integrated analytically against the Chebyshev
183 polynomials to precompute @dfn{modified Chebyshev moments}. Combining
184 the moments with the Chebyshev approximation to the function gives the
185 desired integral. The use of analytic integration for the singular part
186 of the function allows exact cancellations and substantially improves
187 the overall convergence behavior of the integration.
190 @node QNG non-adaptive Gauss-Kronrod integration
191 @section QNG non-adaptive Gauss-Kronrod integration
192 @cindex QNG quadrature algorithm
194 The QNG algorithm is a non-adaptive procedure which uses fixed
195 Gauss-Kronrod-Patterson abscissae to sample the integrand at a maximum of 87
196 points. It is provided for fast integration of smooth functions.
198 @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})
200 This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and
201 87-point integration rules in succession until an estimate of the
202 integral of @math{f} over @math{(a,b)} is achieved within the desired
203 absolute and relative error limits, @var{epsabs} and @var{epsrel}. The
204 function returns the final approximation, @var{result}, an estimate of
205 the absolute error, @var{abserr} and the number of function evaluations
206 used, @var{neval}. The Gauss-Kronrod rules are designed in such a way
207 that each rule uses all the results of its predecessors, in order to
208 minimize the total number of function evaluations.
212 @node QAG adaptive integration
213 @section QAG adaptive integration
214 @cindex QAG quadrature algorithm
216 The QAG algorithm is a simple adaptive integration procedure. The
217 integration region is divided into subintervals, and on each iteration
218 the subinterval with the largest estimated error is bisected. This
219 reduces the overall error rapidly, as the subintervals become
220 concentrated around local difficulties in the integrand. These
221 subintervals are managed by a @code{gsl_integration_workspace} struct,
222 which handles the memory for the subinterval ranges, results and error
225 @deftypefun {gsl_integration_workspace *} gsl_integration_workspace_alloc (size_t @var{n})
226 This function allocates a workspace sufficient to hold @var{n} double
227 precision intervals, their integration results and error estimates.
230 @deftypefun void gsl_integration_workspace_free (gsl_integration_workspace * @var{w})
231 This function frees the memory associated with the workspace @var{w}.
234 @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})
236 This function applies an integration rule adaptively until an estimate
237 of the integral of @math{f} over @math{(a,b)} is achieved within the
238 desired absolute and relative error limits, @var{epsabs} and
239 @var{epsrel}. The function returns the final approximation,
240 @var{result}, and an estimate of the absolute error, @var{abserr}. The
241 integration rule is determined by the value of @var{key}, which should
242 be chosen from the following symbolic names,
245 GSL_INTEG_GAUSS15 (key = 1)
246 GSL_INTEG_GAUSS21 (key = 2)
247 GSL_INTEG_GAUSS31 (key = 3)
248 GSL_INTEG_GAUSS41 (key = 4)
249 GSL_INTEG_GAUSS51 (key = 5)
250 GSL_INTEG_GAUSS61 (key = 6)
254 corresponding to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod
255 rules. The higher-order rules give better accuracy for smooth functions,
256 while lower-order rules save time when the function contains local
257 difficulties, such as discontinuities.
259 On each iteration the adaptive integration strategy bisects the interval
260 with the largest error estimate. The subintervals and their results are
261 stored in the memory provided by @var{workspace}. The maximum number of
262 subintervals is given by @var{limit}, which may not exceed the allocated
263 size of the workspace.
267 @node QAGS adaptive integration with singularities
268 @section QAGS adaptive integration with singularities
269 @cindex QAGS quadrature algorithm
271 The presence of an integrable singularity in the integration region
272 causes an adaptive routine to concentrate new subintervals around the
273 singularity. As the subintervals decrease in size the successive
274 approximations to the integral converge in a limiting fashion. This
275 approach to the limit can be accelerated using an extrapolation
276 procedure. The QAGS algorithm combines adaptive bisection with the Wynn
277 epsilon-algorithm to speed up the integration of many types of
278 integrable singularities.
280 @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})
282 This function applies the Gauss-Kronrod 21-point integration rule
283 adaptively until an estimate of the integral of @math{f} over
284 @math{(a,b)} is achieved within the desired absolute and relative error
285 limits, @var{epsabs} and @var{epsrel}. The results are extrapolated
286 using the epsilon-algorithm, which accelerates the convergence of the
287 integral in the presence of discontinuities and integrable
288 singularities. The function returns the final approximation from the
289 extrapolation, @var{result}, and an estimate of the absolute error,
290 @var{abserr}. The subintervals and their results are stored in the
291 memory provided by @var{workspace}. The maximum number of subintervals
292 is given by @var{limit}, which may not exceed the allocated size of the
297 @node QAGP adaptive integration with known singular points
298 @section QAGP adaptive integration with known singular points
299 @cindex QAGP quadrature algorithm
300 @cindex singular points, specifying positions in quadrature
301 @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})
303 This function applies the adaptive integration algorithm QAGS taking
304 account of the user-supplied locations of singular points. The array
305 @var{pts} of length @var{npts} should contain the endpoints of the
306 integration ranges defined by the integration region and locations of
307 the singularities. For example, to integrate over the region
308 @math{(a,b)} with break-points at @math{x_1, x_2, x_3} (where
309 @math{a < x_1 < x_2 < x_3 < b}) the following @var{pts} array should be used
323 If you know the locations of the singular points in the integration
324 region then this routine will be faster than @code{QAGS}.
328 @node QAGI adaptive integration on infinite intervals
329 @section QAGI adaptive integration on infinite intervals
330 @cindex QAGI quadrature algorithm
332 @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})
334 This function computes the integral of the function @var{f} over the
335 infinite interval @math{(-\infty,+\infty)}. The integral is mapped onto the
336 semi-open interval @math{(0,1]} using the transformation @math{x = (1-t)/t},
340 \int_{-\infty}^{+\infty} dx \, f(x)
341 = \int_0^1 dt \, (f((1-t)/t) + f(-(1-t)/t))/t^2.
348 \int_@{-\infty@}^@{+\infty@} dx f(x) =
349 \int_0^1 dt (f((1-t)/t) + f((-1+t)/t))/t^2.
354 It is then integrated using the QAGS algorithm. The normal 21-point
355 Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, because the
356 transformation can generate an integrable singularity at the origin. In
357 this case a lower-order rule is more efficient.
360 @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})
362 This function computes the integral of the function @var{f} over the
363 semi-infinite interval @math{(a,+\infty)}. The integral is mapped onto the
364 semi-open interval @math{(0,1]} using the transformation @math{x = a + (1-t)/t},
368 \int_{a}^{+\infty} dx \, f(x)
369 = \int_0^1 dt \, f(a + (1-t)/t)/t^2
376 \int_@{a@}^@{+\infty@} dx f(x) =
377 \int_0^1 dt f(a + (1-t)/t)/t^2
382 and then integrated using the QAGS algorithm.
385 @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})
386 This function computes the integral of the function @var{f} over the
387 semi-infinite interval @math{(-\infty,b)}. The integral is mapped onto the
388 semi-open interval @math{(0,1]} using the transformation @math{x = b - (1-t)/t},
392 \int_{-\infty}^{b} dx \, f(x)
393 = \int_0^1 dt \, f(b - (1-t)/t)/t^2
400 \int_@{-\infty@}^@{b@} dx f(x) =
401 \int_0^1 dt f(b - (1-t)/t)/t^2
406 and then integrated using the QAGS algorithm.
409 @node QAWC adaptive integration for Cauchy principal values
410 @section QAWC adaptive integration for Cauchy principal values
411 @cindex QAWC quadrature algorithm
412 @cindex Cauchy principal value, by numerical quadrature
413 @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})
415 This function computes the Cauchy principal value of the integral of
416 @math{f} over @math{(a,b)}, with a singularity at @var{c},
420 I = \int_a^b dx\, {f(x) \over x - c}
421 = \lim_{\epsilon \to 0}
423 \int_a^{c-\epsilon} dx\, {f(x) \over x - c}
425 \int_{c+\epsilon}^b dx\, {f(x) \over x - c}
433 I = \int_a^b dx f(x) / (x - c)
438 The adaptive bisection algorithm of QAG is used, with modifications to
439 ensure that subdivisions do not occur at the singular point @math{x = c}.
440 When a subinterval contains the point @math{x = c} or is close to
441 it then a special 25-point modified Clenshaw-Curtis rule is used to control
442 the singularity. Further away from the
443 singularity the algorithm uses an ordinary 15-point Gauss-Kronrod
448 @node QAWS adaptive integration for singular functions
449 @section QAWS adaptive integration for singular functions
450 @cindex QAWS quadrature algorithm
451 @cindex singular functions, numerical integration of
452 The QAWS algorithm is designed for integrands with algebraic-logarithmic
453 singularities at the end-points of an integration region. In order to
454 work efficiently the algorithm requires a precomputed table of
457 @deftypefun {gsl_integration_qaws_table *} gsl_integration_qaws_table_alloc (double @var{alpha}, double @var{beta}, int @var{mu}, int @var{nu})
459 This function allocates space for a @code{gsl_integration_qaws_table}
460 struct describing a singular weight function
461 @math{W(x)} with the parameters @math{(\alpha, \beta, \mu, \nu)},
465 W(x) = (x - a)^\alpha (b - x)^\beta \log^\mu (x - a) \log^\nu (b - x)
472 W(x) = (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x)
477 where @math{\alpha > -1}, @math{\beta > -1}, and @math{\mu = 0, 1},
478 @math{\nu = 0, 1}. The weight function can take four different forms
479 depending on the values of @math{\mu} and @math{\nu},
484 W(x) = (x - a)^\alpha (b - x)^\beta
485 \hfill~ (\mu = 0, \nu = 0) \cr
486 W(x) = (x - a)^\alpha (b - x)^\beta \log(x - a)
487 \hfill~ (\mu = 1, \nu = 0) \cr
488 W(x) = (x - a)^\alpha (b - x)^\beta \log(b - x)
489 \hfill~ (\mu = 0, \nu = 1) \cr
490 W(x) = (x - a)^\alpha (b - x)^\beta \log(x - a) \log(b - x)
491 \hfill~ (\mu = 1, \nu = 1)
499 W(x) = (x-a)^alpha (b-x)^beta (mu = 0, nu = 0)
500 W(x) = (x-a)^alpha (b-x)^beta log(x-a) (mu = 1, nu = 0)
501 W(x) = (x-a)^alpha (b-x)^beta log(b-x) (mu = 0, nu = 1)
502 W(x) = (x-a)^alpha (b-x)^beta log(x-a) log(b-x) (mu = 1, nu = 1)
507 The singular points @math{(a,b)} do not have to be specified until the
508 integral is computed, where they are the endpoints of the integration
511 The function returns a pointer to the newly allocated table
512 @code{gsl_integration_qaws_table} if no errors were detected, and 0 in
516 @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})
517 This function modifies the parameters @math{(\alpha, \beta, \mu, \nu)} of
518 an existing @code{gsl_integration_qaws_table} struct @var{t}.
521 @deftypefun void gsl_integration_qaws_table_free (gsl_integration_qaws_table * @var{t})
522 This function frees all the memory associated with the
523 @code{gsl_integration_qaws_table} struct @var{t}.
526 @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})
528 This function computes the integral of the function @math{f(x)} over the
529 interval @math{(a,b)} with the singular weight function
530 @math{(x-a)^\alpha (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x)}. The parameters
531 of the weight function @math{(\alpha, \beta, \mu, \nu)} are taken from the
532 table @var{t}. The integral is,
536 I = \int_a^b dx\, f(x) (x - a)^\alpha (b - x)^\beta
537 \log^\mu (x - a) \log^\nu (b - x).
544 I = \int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x).
549 The adaptive bisection algorithm of QAG is used. When a subinterval
550 contains one of the endpoints then a special 25-point modified
551 Clenshaw-Curtis rule is used to control the singularities. For
552 subintervals which do not include the endpoints an ordinary 15-point
553 Gauss-Kronrod integration rule is used.
557 @node QAWO adaptive integration for oscillatory functions
558 @section QAWO adaptive integration for oscillatory functions
559 @cindex QAWO quadrature algorithm
560 @cindex oscillatory functions, numerical integration of
561 The QAWO algorithm is designed for integrands with an oscillatory
562 factor, @math{\sin(\omega x)} or @math{\cos(\omega x)}. In order to
563 work efficiently the algorithm requires a table of Chebyshev moments
564 which must be pre-computed with calls to the functions below.
566 @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})
568 This function allocates space for a @code{gsl_integration_qawo_table}
569 struct and its associated workspace describing a sine or cosine weight
570 function @math{W(x)} with the parameters @math{(\omega, L)},
575 W(x) & = \left\{\matrix{\sin(\omega x) \cr \cos(\omega x)} \right\}
589 The parameter @var{L} must be the length of the interval over which the
590 function will be integrated @math{L = b - a}. The choice of sine or
591 cosine is made with the parameter @var{sine} which should be chosen from
592 one of the two following symbolic values:
600 The @code{gsl_integration_qawo_table} is a table of the trigonometric
601 coefficients required in the integration process. The parameter @var{n}
602 determines the number of levels of coefficients that are computed. Each
603 level corresponds to one bisection of the interval @math{L}, so that
604 @var{n} levels are sufficient for subintervals down to the length
605 @math{L/2^n}. The integration routine @code{gsl_integration_qawo}
606 returns the error @code{GSL_ETABLE} if the number of levels is
607 insufficient for the requested accuracy.
611 @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})
612 This function changes the parameters @var{omega}, @var{L} and @var{sine}
613 of the existing workspace @var{t}.
616 @deftypefun int gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * @var{t}, double @var{L})
617 This function allows the length parameter @var{L} of the workspace
618 @var{t} to be changed.
621 @deftypefun void gsl_integration_qawo_table_free (gsl_integration_qawo_table * @var{t})
622 This function frees all the memory associated with the workspace @var{t}.
625 @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})
627 This function uses an adaptive algorithm to compute the integral of
628 @math{f} over @math{(a,b)} with the weight function
629 @math{\sin(\omega x)} or @math{\cos(\omega x)} defined
630 by the table @var{wf},
635 I & = \int_a^b dx\, f(x) \left\{ \matrix{\sin(\omega x) \cr \cos(\omega x)}\right\}
643 I = \int_a^b dx f(x) sin(omega x)
644 I = \int_a^b dx f(x) cos(omega x)
649 The results are extrapolated using the epsilon-algorithm to accelerate
650 the convergence of the integral. The function returns the final
651 approximation from the extrapolation, @var{result}, and an estimate of
652 the absolute error, @var{abserr}. The subintervals and their results are
653 stored in the memory provided by @var{workspace}. The maximum number of
654 subintervals is given by @var{limit}, which may not exceed the allocated
655 size of the workspace.
657 Those subintervals with ``large'' widths @math{d} where @math{d\omega > 4} are
658 computed using a 25-point Clenshaw-Curtis integration rule, which handles the
659 oscillatory behavior. Subintervals with a ``small'' widths where
660 @math{d\omega < 4} are computed using a 15-point Gauss-Kronrod integration.
664 @node QAWF adaptive integration for Fourier integrals
665 @section QAWF adaptive integration for Fourier integrals
666 @cindex QAWF quadrature algorithm
667 @cindex Fourier integrals, numerical
669 @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})
671 This function attempts to compute a Fourier integral of the function
672 @var{f} over the semi-infinite interval @math{[a,+\infty)}.
677 I & = \int_a^{+\infty} dx\, f(x) \left\{ \matrix{ \sin(\omega x) \cr
678 \cos(\omega x) } \right\}
686 I = \int_a^@{+\infty@} dx f(x) sin(omega x)
687 I = \int_a^@{+\infty@} dx f(x) cos(omega x)
691 The parameter @math{\omega} and choice of @math{\sin} or @math{\cos} is
692 taken from the table @var{wf} (the length @var{L} can take any value,
693 since it is overridden by this function to a value appropriate for the
694 fourier integration). The integral is computed using the QAWO algorithm
695 over each of the subintervals,
700 C_1 & = [a, a + c] \cr
701 C_2 & = [a + c, a + 2c] \cr
703 C_k & = [a + (k-1) c, a + k c]
712 C_2 = [a + c, a + 2 c]
714 C_k = [a + (k-1) c, a + k c]
720 @c{$c = (2 \,\hbox{floor}(|\omega|) + 1) \pi/|\omega|$}
721 @math{c = (2 floor(|\omega|) + 1) \pi/|\omega|}. The width @math{c} is
722 chosen to cover an odd number of periods so that the contributions from
723 the intervals alternate in sign and are monotonically decreasing when
724 @var{f} is positive and monotonically decreasing. The sum of this
725 sequence of contributions is accelerated using the epsilon-algorithm.
727 This function works to an overall absolute tolerance of
728 @var{abserr}. The following strategy is used: on each interval
729 @math{C_k} the algorithm tries to achieve the tolerance
733 TOL_k = u_k \hbox{\it abserr}
746 @c{$u_k = (1 - p)p^{k-1}$}
747 @math{u_k = (1 - p)p^@{k-1@}} and @math{p = 9/10}.
748 The sum of the geometric series of contributions from each interval
749 gives an overall tolerance of @var{abserr}.
751 If the integration of a subinterval leads to difficulties then the
752 accuracy requirement for subsequent intervals is relaxed,
756 TOL_k = u_k \max(\hbox{\it abserr}, \max_{i<k}\{E_i\})
763 TOL_k = u_k max(abserr, max_@{i<k@}@{E_i@})
768 where @math{E_k} is the estimated error on the interval @math{C_k}.
770 The subintervals and their results are stored in the memory provided by
771 @var{workspace}. The maximum number of subintervals is given by
772 @var{limit}, which may not exceed the allocated size of the workspace.
773 The integration over each subinterval uses the memory provided by
774 @var{cycle_workspace} as workspace for the QAWO algorithm.
778 @node Numerical integration error codes
781 In addition to the standard error codes for invalid arguments the
782 functions can return the following values,
786 the maximum number of subdivisions was exceeded.
788 cannot reach tolerance because of roundoff error,
789 or roundoff error was detected in the extrapolation table.
791 a non-integrable singularity or other bad integrand behavior was found
792 in the integration interval.
794 the integral is divergent, or too slowly convergent to be integrated
798 @node Numerical integration examples
801 The integrator @code{QAGS} will handle a large class of definite
802 integrals. For example, consider the following integral, which has a
803 algebraic-logarithmic singularity at the origin,
807 \int_0^1 x^{-1/2} \log(x) \,dx = -4
814 \int_0^1 x^@{-1/2@} log(x) dx = -4
819 The program below computes this integral to a relative accuracy bound of
823 @verbatiminclude examples/integration.c
827 The results below show that the desired accuracy is achieved after 8
832 @verbatiminclude examples/integration.out
836 In fact, the extrapolation procedure used by @code{QAGS} produces an
837 accuracy of almost twice as many digits. The error estimate returned by
838 the extrapolation procedure is larger than the actual error, giving a
839 margin of safety of one order of magnitude.
842 @node Numerical integration References and Further Reading
843 @section References and Further Reading
845 The following book is the definitive reference for @sc{quadpack}, and was
846 written by the original authors. It provides descriptions of the
847 algorithms, program listings, test programs and examples. It also
848 includes useful advice on numerical integration and many references to
849 the numerical integration literature used in developing @sc{quadpack}.
853 R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. Kahaner.
854 @cite{@sc{quadpack} A subroutine package for automatic integration}
855 Springer Verlag, 1983.