|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: quadrature
|
|
Packit |
67cb25 |
single: numerical integration (quadrature)
|
|
Packit |
67cb25 |
single: integration, numerical (quadrature)
|
|
Packit |
67cb25 |
single: QUADPACK
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*********************
|
|
Packit |
67cb25 |
Numerical Integration
|
|
Packit |
67cb25 |
*********************
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: include.rst
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This chapter describes routines for performing numerical integration
|
|
Packit |
67cb25 |
(quadrature) of a function in one dimension. There are routines for
|
|
Packit |
67cb25 |
adaptive and non-adaptive integration of general functions, with
|
|
Packit |
67cb25 |
specialised routines for specific cases. These include integration over
|
|
Packit |
67cb25 |
infinite and semi-infinite ranges, singular integrals, including
|
|
Packit |
67cb25 |
logarithmic singularities, computation of Cauchy principal values and
|
|
Packit |
67cb25 |
oscillatory integrals. The library reimplements the algorithms used in
|
|
Packit |
67cb25 |
|quadpack|, a numerical integration package written by Piessens,
|
|
Packit |
67cb25 |
de Doncker-Kapenga, Ueberhuber and Kahaner. Fortran code for |quadpack| is
|
|
Packit |
67cb25 |
available on Netlib. Also included are non-adaptive, fixed-order
|
|
Packit |
67cb25 |
Gauss-Legendre integration routines with high precision coefficients, as
|
|
Packit |
67cb25 |
well as fixed-order quadrature rules for a variety of weighting functions
|
|
Packit |
67cb25 |
from IQPACK.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this chapter are declared in the header file
|
|
Packit |
67cb25 |
:file:`gsl_integration.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Introduction
|
|
Packit |
67cb25 |
============
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Each algorithm computes an approximation to a definite integral of the
|
|
Packit |
67cb25 |
form,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: I = \int_a^b f(x) w(x) dx
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`w(x)` is a weight function (for general integrands :math:`w(x) = 1`).
|
|
Packit |
67cb25 |
The user provides absolute and relative error bounds
|
|
Packit |
67cb25 |
:math:`(epsabs, epsrel)` which specify the following accuracy requirement,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: |RESULT - I| \leq \max{(epsabs, epsrel |I|)}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: |RESULT - I| <= max(epsabs, epsrel |I|)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`RESULT` is the numerical approximation obtained by the
|
|
Packit |
67cb25 |
algorithm. The algorithms attempt to estimate the absolute error
|
|
Packit |
67cb25 |
:math:`ABSERR = |RESULT - I|` in such a way that the following inequality
|
|
Packit |
67cb25 |
holds,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: |RESULT - I| \leq ABSERR \leq \max{(epsabs, epsrel |I|)}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: |RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In short, the routines return the first approximation
|
|
Packit |
67cb25 |
which has an absolute error smaller than
|
|
Packit |
67cb25 |
:math:`epsabs` or a relative error smaller than :math:`epsrel`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that this is an *either-or* constraint,
|
|
Packit |
67cb25 |
not simultaneous. To compute to a specified absolute error, set
|
|
Packit |
67cb25 |
:math:`epsrel` to zero. To compute to a specified relative error,
|
|
Packit |
67cb25 |
set :math:`epsabs` to zero.
|
|
Packit |
67cb25 |
The routines will fail to converge if the error bounds are too
|
|
Packit |
67cb25 |
stringent, but always return the best approximation obtained up to
|
|
Packit |
67cb25 |
that stage.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithms in |quadpack| use a naming convention based on the
|
|
Packit |
67cb25 |
following letters::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Q - quadrature routine
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
N - non-adaptive integrator
|
|
Packit |
67cb25 |
A - adaptive integrator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
G - general integrand (user-defined)
|
|
Packit |
67cb25 |
W - weight function with integrand
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
S - singularities can be more readily integrated
|
|
Packit |
67cb25 |
P - points of special difficulty can be supplied
|
|
Packit |
67cb25 |
I - infinite range of integration
|
|
Packit |
67cb25 |
O - oscillatory weight function, cos or sin
|
|
Packit |
67cb25 |
F - Fourier integral
|
|
Packit |
67cb25 |
C - Cauchy principal value
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithms are built on pairs of quadrature rules, a higher order
|
|
Packit |
67cb25 |
rule and a lower order rule. The higher order rule is used to compute
|
|
Packit |
67cb25 |
the best approximation to an integral over a small range. The
|
|
Packit |
67cb25 |
difference between the results of the higher order rule and the lower
|
|
Packit |
67cb25 |
order rule gives an estimate of the error in the approximation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Gauss-Kronrod quadrature
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Integrands without weight functions
|
|
Packit |
67cb25 |
-----------------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithms for general functions (without a weight function) are
|
|
Packit |
67cb25 |
based on Gauss-Kronrod rules.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A Gauss-Kronrod rule begins with a classical Gaussian quadrature rule of
|
|
Packit |
67cb25 |
order :math:`m`. This is extended with additional points between each of
|
|
Packit |
67cb25 |
the abscissae to give a higher order Kronrod rule of order :math:`2m + 1`.
|
|
Packit |
67cb25 |
The Kronrod rule is efficient because it reuses existing function
|
|
Packit |
67cb25 |
evaluations from the Gaussian rule.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The higher order Kronrod rule is used as the best approximation to the
|
|
Packit |
67cb25 |
integral, and the difference between the two rules is used as an
|
|
Packit |
67cb25 |
estimate of the error in the approximation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Integrands with weight functions
|
|
Packit |
67cb25 |
--------------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: Clenshaw-Curtis quadrature
|
|
Packit |
67cb25 |
single: Modified Clenshaw-Curtis quadrature
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For integrands with weight functions the algorithms use Clenshaw-Curtis
|
|
Packit |
67cb25 |
quadrature rules.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A Clenshaw-Curtis rule begins with an :math:`n`-th order Chebyshev
|
|
Packit |
67cb25 |
polynomial approximation to the integrand. This polynomial can be
|
|
Packit |
67cb25 |
integrated exactly to give an approximation to the integral of the
|
|
Packit |
67cb25 |
original function. The Chebyshev expansion can be extended to higher
|
|
Packit |
67cb25 |
orders to improve the approximation and provide an estimate of the
|
|
Packit |
67cb25 |
error.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Integrands with singular weight functions
|
|
Packit |
67cb25 |
-----------------------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The presence of singularities (or other behavior) in the integrand can
|
|
Packit |
67cb25 |
cause slow convergence in the Chebyshev approximation. The modified
|
|
Packit |
67cb25 |
Clenshaw-Curtis rules used in |quadpack| separate out several common
|
|
Packit |
67cb25 |
weight functions which cause slow convergence.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These weight functions are integrated analytically against the Chebyshev
|
|
Packit |
67cb25 |
polynomials to precompute *modified Chebyshev moments*. Combining
|
|
Packit |
67cb25 |
the moments with the Chebyshev approximation to the function gives the
|
|
Packit |
67cb25 |
desired integral. The use of analytic integration for the singular part
|
|
Packit |
67cb25 |
of the function allows exact cancellations and substantially improves
|
|
Packit |
67cb25 |
the overall convergence behavior of the integration.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QNG non-adaptive Gauss-Kronrod integration
|
|
Packit |
67cb25 |
==========================================
|
|
Packit |
67cb25 |
.. index:: QNG quadrature algorithm
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The QNG algorithm is a non-adaptive procedure which uses fixed
|
|
Packit |
67cb25 |
Gauss-Kronrod-Patterson abscissae to sample the integrand at a maximum of 87
|
|
Packit |
67cb25 |
points. It is provided for fast integration of smooth functions.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qng (const gsl_function * f, double a, double b, double epsabs, double epsrel, double * result, double * abserr, size_t * neval)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and
|
|
Packit |
67cb25 |
87-point integration rules in succession until an estimate of the
|
|
Packit |
67cb25 |
integral of :math:`f` over :math:`(a,b)` is achieved within the desired
|
|
Packit |
67cb25 |
absolute and relative error limits, :data:`epsabs` and :data:`epsrel`. The
|
|
Packit |
67cb25 |
function returns the final approximation, :data:`result`, an estimate of
|
|
Packit |
67cb25 |
the absolute error, :data:`abserr` and the number of function evaluations
|
|
Packit |
67cb25 |
used, :data:`neval`. The Gauss-Kronrod rules are designed in such a way
|
|
Packit |
67cb25 |
that each rule uses all the results of its predecessors, in order to
|
|
Packit |
67cb25 |
minimize the total number of function evaluations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAG adaptive integration
|
|
Packit |
67cb25 |
========================
|
|
Packit |
67cb25 |
.. index:: QAG quadrature algorithm
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The QAG algorithm is a simple adaptive integration procedure. The
|
|
Packit |
67cb25 |
integration region is divided into subintervals, and on each iteration
|
|
Packit |
67cb25 |
the subinterval with the largest estimated error is bisected. This
|
|
Packit |
67cb25 |
reduces the overall error rapidly, as the subintervals become
|
|
Packit |
67cb25 |
concentrated around local difficulties in the integrand. These
|
|
Packit |
67cb25 |
subintervals are managed by the following struct,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_integration_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This workspace handles the memory for the subinterval ranges, results and error
|
|
Packit |
67cb25 |
estimates.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: gsl_integration_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_workspace * gsl_integration_workspace_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace sufficient to hold :data:`n` double
|
|
Packit |
67cb25 |
precision intervals, their integration results and error estimates.
|
|
Packit |
67cb25 |
One workspace may be used multiple times as all necessary reinitialization
|
|
Packit |
67cb25 |
is performed automatically by the integration routines.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_workspace_free (gsl_integration_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qag (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies an integration rule adaptively until an estimate
|
|
Packit |
67cb25 |
of the integral of :math:`f` over :math:`(a,b)` is achieved within the
|
|
Packit |
67cb25 |
desired absolute and relative error limits, :data:`epsabs` and
|
|
Packit |
67cb25 |
:data:`epsrel`. The function returns the final approximation,
|
|
Packit |
67cb25 |
:data:`result`, and an estimate of the absolute error, :data:`abserr`. The
|
|
Packit |
67cb25 |
integration rule is determined by the value of :data:`key`, which should
|
|
Packit |
67cb25 |
be chosen from the following symbolic names,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* .. macro:: GSL_INTEG_GAUSS15 (key = 1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* .. macro:: GSL_INTEG_GAUSS21 (key = 2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* .. macro:: GSL_INTEG_GAUSS31 (key = 3)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* .. macro:: GSL_INTEG_GAUSS41 (key = 4)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* .. macro:: GSL_INTEG_GAUSS51 (key = 5)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* .. macro:: GSL_INTEG_GAUSS61 (key = 6)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
corresponding to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod
|
|
Packit |
67cb25 |
rules. The higher-order rules give better accuracy for smooth functions,
|
|
Packit |
67cb25 |
while lower-order rules save time when the function contains local
|
|
Packit |
67cb25 |
difficulties, such as discontinuities.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
On each iteration the adaptive integration strategy bisects the interval
|
|
Packit |
67cb25 |
with the largest error estimate. The subintervals and their results are
|
|
Packit |
67cb25 |
stored in the memory provided by :data:`workspace`. The maximum number of
|
|
Packit |
67cb25 |
subintervals is given by :data:`limit`, which may not exceed the allocated
|
|
Packit |
67cb25 |
size of the workspace.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAGS adaptive integration with singularities
|
|
Packit |
67cb25 |
============================================
|
|
Packit |
67cb25 |
.. index:: QAGS quadrature algorithm
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The presence of an integrable singularity in the integration region
|
|
Packit |
67cb25 |
causes an adaptive routine to concentrate new subintervals around the
|
|
Packit |
67cb25 |
singularity. As the subintervals decrease in size the successive
|
|
Packit |
67cb25 |
approximations to the integral converge in a limiting fashion. This
|
|
Packit |
67cb25 |
approach to the limit can be accelerated using an extrapolation
|
|
Packit |
67cb25 |
procedure. The QAGS algorithm combines adaptive bisection with the Wynn
|
|
Packit |
67cb25 |
epsilon-algorithm to speed up the integration of many types of
|
|
Packit |
67cb25 |
integrable singularities.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qags (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Gauss-Kronrod 21-point integration rule
|
|
Packit |
67cb25 |
adaptively until an estimate of the integral of :math:`f` over
|
|
Packit |
67cb25 |
:math:`(a,b)` is achieved within the desired absolute and relative error
|
|
Packit |
67cb25 |
limits, :data:`epsabs` and :data:`epsrel`. The results are extrapolated
|
|
Packit |
67cb25 |
using the epsilon-algorithm, which accelerates the convergence of the
|
|
Packit |
67cb25 |
integral in the presence of discontinuities and integrable
|
|
Packit |
67cb25 |
singularities. The function returns the final approximation from the
|
|
Packit |
67cb25 |
extrapolation, :data:`result`, and an estimate of the absolute error,
|
|
Packit |
67cb25 |
:data:`abserr`. The subintervals and their results are stored in the
|
|
Packit |
67cb25 |
memory provided by :data:`workspace`. The maximum number of subintervals
|
|
Packit |
67cb25 |
is given by :data:`limit`, which may not exceed the allocated size of the
|
|
Packit |
67cb25 |
workspace.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAGP adaptive integration with known singular points
|
|
Packit |
67cb25 |
====================================================
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: QAGP quadrature algorithm
|
|
Packit |
67cb25 |
single: singular points, specifying positions in quadrature
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qagp (const gsl_function * f, double * pts, size_t npts, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the adaptive integration algorithm QAGS taking
|
|
Packit |
67cb25 |
account of the user-supplied locations of singular points. The array
|
|
Packit |
67cb25 |
:data:`pts` of length :data:`npts` should contain the endpoints of the
|
|
Packit |
67cb25 |
integration ranges defined by the integration region and locations of
|
|
Packit |
67cb25 |
the singularities. For example, to integrate over the region
|
|
Packit |
67cb25 |
:math:`(a,b)` with break-points at :math:`x_1, x_2, x_3` (where
|
|
Packit |
67cb25 |
:math:`a < x_1 < x_2 < x_3 < b`) the following :data:`pts` array should be used::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
pts[0] = a
|
|
Packit |
67cb25 |
pts[1] = x_1
|
|
Packit |
67cb25 |
pts[2] = x_2
|
|
Packit |
67cb25 |
pts[3] = x_3
|
|
Packit |
67cb25 |
pts[4] = b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :data:`npts` = 5.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If you know the locations of the singular points in the integration
|
|
Packit |
67cb25 |
region then this routine will be faster than :func:`gsl_integration_qags`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAGI adaptive integration on infinite intervals
|
|
Packit |
67cb25 |
===============================================
|
|
Packit |
67cb25 |
.. index:: QAGI quadrature algorithm
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qagi (gsl_function * f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the integral of the function :data:`f` over the
|
|
Packit |
67cb25 |
infinite interval :math:`(-\infty,+\infty)`. The integral is mapped onto the
|
|
Packit |
67cb25 |
semi-open interval :math:`(0,1]` using the transformation :math:`x = (1-t)/t`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \int_{-\infty}^{+\infty} dx f(x) = \int_0^1 dt (f((1-t)/t) + f(-(1-t)/t))/t^2.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
It is then integrated using the QAGS algorithm. The normal 21-point
|
|
Packit |
67cb25 |
Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, because the
|
|
Packit |
67cb25 |
transformation can generate an integrable singularity at the origin. In
|
|
Packit |
67cb25 |
this case a lower-order rule is more efficient.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qagiu (gsl_function * f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the integral of the function :data:`f` over the
|
|
Packit |
67cb25 |
semi-infinite interval :math:`(a,+\infty)`. The integral is mapped onto the
|
|
Packit |
67cb25 |
semi-open interval :math:`(0,1]` using the transformation :math:`x = a + (1-t)/t`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \int_{a}^{+\infty} dx f(x) = \int_0^1 dt f(a + (1-t)/t)/t^2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and then integrated using the QAGS algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qagil (gsl_function * f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the integral of the function :data:`f` over the
|
|
Packit |
67cb25 |
semi-infinite interval :math:`(-\infty,b)`. The integral is mapped onto the
|
|
Packit |
67cb25 |
semi-open interval :math:`(0,1]` using the transformation :math:`x = b - (1-t)/t`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \int_{-\infty}^{b} dx f(x) = \int_0^1 dt f(b - (1-t)/t)/t^2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and then integrated using the QAGS algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAWC adaptive integration for Cauchy principal values
|
|
Packit |
67cb25 |
=====================================================
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: QAWC quadrature algorithm
|
|
Packit |
67cb25 |
single: Cauchy principal value, by numerical quadrature
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qawc (gsl_function * f, double a, double b, double c, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the Cauchy principal value of the integral of
|
|
Packit |
67cb25 |
:math:`f` over :math:`(a,b)`, with a singularity at :data:`c`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
I = \int_a^b dx\, {f(x) \over x - c}
|
|
Packit |
67cb25 |
= \lim_{\epsilon \to 0}
|
|
Packit |
67cb25 |
\left\{
|
|
Packit |
67cb25 |
\int_a^{c-\epsilon} dx\, {f(x) \over x - c}
|
|
Packit |
67cb25 |
+
|
|
Packit |
67cb25 |
\int_{c+\epsilon}^b dx\, {f(x) \over x - c}
|
|
Packit |
67cb25 |
\right\}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: I = \int_a^b dx f(x) / (x - c)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The adaptive bisection algorithm of QAG is used, with modifications to
|
|
Packit |
67cb25 |
ensure that subdivisions do not occur at the singular point :math:`x = c`.
|
|
Packit |
67cb25 |
When a subinterval contains the point :math:`x = c` or is close to
|
|
Packit |
67cb25 |
it then a special 25-point modified Clenshaw-Curtis rule is used to control
|
|
Packit |
67cb25 |
the singularity. Further away from the
|
|
Packit |
67cb25 |
singularity the algorithm uses an ordinary 15-point Gauss-Kronrod
|
|
Packit |
67cb25 |
integration rule.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAWS adaptive integration for singular functions
|
|
Packit |
67cb25 |
================================================
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: QAWS quadrature algorithm
|
|
Packit |
67cb25 |
single: singular functions, numerical integration of
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The QAWS algorithm is designed for integrands with algebraic-logarithmic
|
|
Packit |
67cb25 |
singularities at the end-points of an integration region. In order to
|
|
Packit |
67cb25 |
work efficiently the algorithm requires a precomputed table of
|
|
Packit |
67cb25 |
Chebyshev moments.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_integration_qaws_table
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This structure contains precomputed quantities for the QAWS algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_qaws_table * gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates space for a :type:`gsl_integration_qaws_table`
|
|
Packit |
67cb25 |
struct describing a singular weight function
|
|
Packit |
67cb25 |
:math:`w(x)` with the parameters :math:`(\alpha, \beta, \mu, \nu)`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: w(x) = (x - a)^\alpha (b - x)^\beta \log^\mu (x - a) \log^\nu (b - x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`\alpha > -1`, :math:`\beta > -1`, and :math:`\mu = 0, 1`,
|
|
Packit |
67cb25 |
:math:`\nu = 0, 1`. The weight function can take four different forms
|
|
Packit |
67cb25 |
depending on the values of :math:`\mu` and :math:`\nu`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
============================================================ =================
|
|
Packit |
67cb25 |
Weight function :math:`w(x)` :math:`(\mu,\nu)`
|
|
Packit |
67cb25 |
============================================================ =================
|
|
Packit |
67cb25 |
:math:`(x - a)^\alpha (b - x)^\beta` :math:`(0,0)`
|
|
Packit |
67cb25 |
:math:`(x - a)^\alpha (b - x)^\beta \log{(x-a)}` :math:`(1,0)`
|
|
Packit |
67cb25 |
:math:`(x - a)^\alpha (b - x)^\beta \log{(b-x)}` :math:`(0,1)`
|
|
Packit |
67cb25 |
:math:`(x - a)^\alpha (b - x)^\beta \log{(x-a)} \log{(b-x)}` :math:`(1,1)`
|
|
Packit |
67cb25 |
============================================================ =================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The singular points :math:`(a,b)` do not have to be specified until the
|
|
Packit |
67cb25 |
integral is computed, where they are the endpoints of the integration
|
|
Packit |
67cb25 |
range.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The function returns a pointer to the newly allocated table
|
|
Packit |
67cb25 |
:type:`gsl_integration_qaws_table` if no errors were detected, and 0 in
|
|
Packit |
67cb25 |
the case of error.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qaws_table_set (gsl_integration_qaws_table * t, double alpha, double beta, int mu, int nu)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function modifies the parameters :math:`(\alpha, \beta, \mu, \nu)` of
|
|
Packit |
67cb25 |
an existing :type:`gsl_integration_qaws_table` struct :data:`t`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_qaws_table_free (gsl_integration_qaws_table * t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees all the memory associated with the
|
|
Packit |
67cb25 |
:type:`gsl_integration_qaws_table` struct :data:`t`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qaws (gsl_function * f, const double a, const double b, gsl_integration_qaws_table * t, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the integral of the function :math:`f(x)` over the
|
|
Packit |
67cb25 |
interval :math:`(a,b)` with the singular weight function
|
|
Packit |
67cb25 |
:math:`(x-a)^\alpha (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x)`. The parameters
|
|
Packit |
67cb25 |
of the weight function :math:`(\alpha, \beta, \mu, \nu)` are taken from the
|
|
Packit |
67cb25 |
table :data:`t`. The integral is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: I = \int_a^b dx f(x) (x - a)^\alpha (b - x)^\beta \log^\mu (x - a) \log^\nu (b - x).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The adaptive bisection algorithm of QAG is used. When a subinterval
|
|
Packit |
67cb25 |
contains one of the endpoints then a special 25-point modified
|
|
Packit |
67cb25 |
Clenshaw-Curtis rule is used to control the singularities. For
|
|
Packit |
67cb25 |
subintervals which do not include the endpoints an ordinary 15-point
|
|
Packit |
67cb25 |
Gauss-Kronrod integration rule is used.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAWO adaptive integration for oscillatory functions
|
|
Packit |
67cb25 |
===================================================
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: QAWO quadrature algorithm
|
|
Packit |
67cb25 |
single: oscillatory functions, numerical integration of
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The QAWO algorithm is designed for integrands with an oscillatory
|
|
Packit |
67cb25 |
factor, :math:`\sin(\omega x)` or :math:`\cos(\omega x)`. In order to
|
|
Packit |
67cb25 |
work efficiently the algorithm requires a table of Chebyshev moments
|
|
Packit |
67cb25 |
which must be pre-computed with calls to the functions below.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: gsl_integration_qawo_table
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_qawo_table * gsl_integration_qawo_table_alloc (double omega, double L, enum gsl_integration_qawo_enum sine, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates space for a :type:`gsl_integration_qawo_table`
|
|
Packit |
67cb25 |
struct and its associated workspace describing a sine or cosine weight
|
|
Packit |
67cb25 |
function :math:`w(x)` with the parameters :math:`(\omega, L)`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w(x) =
|
|
Packit |
67cb25 |
\left\{
|
|
Packit |
67cb25 |
\begin{array}{c}
|
|
Packit |
67cb25 |
\sin{(\omega x)} \\
|
|
Packit |
67cb25 |
\cos{(\omega x)} \\
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right\}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
| w(x) = sin(\omega x)
|
|
Packit |
67cb25 |
| w(x) = cos(\omega x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The parameter :data:`L` must be the length of the interval over which the
|
|
Packit |
67cb25 |
function will be integrated :math:`L = b - a`. The choice of sine or
|
|
Packit |
67cb25 |
cosine is made with the parameter :data:`sine` which should be chosen from
|
|
Packit |
67cb25 |
one of the two following symbolic values:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. macro:: GSL_INTEG_COSINE
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. macro:: GSL_INTEG_SINE
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The :type:`gsl_integration_qawo_table` is a table of the trigonometric
|
|
Packit |
67cb25 |
coefficients required in the integration process. The parameter :data:`n`
|
|
Packit |
67cb25 |
determines the number of levels of coefficients that are computed. Each
|
|
Packit |
67cb25 |
level corresponds to one bisection of the interval :math:`L`, so that
|
|
Packit |
67cb25 |
:data:`n` levels are sufficient for subintervals down to the length
|
|
Packit |
67cb25 |
:math:`L/2^n`. The integration routine :func:`gsl_integration_qawo`
|
|
Packit |
67cb25 |
returns the error :code:`GSL_ETABLE` if the number of levels is
|
|
Packit |
67cb25 |
insufficient for the requested accuracy.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qawo_table_set (gsl_integration_qawo_table * t, double omega, double L, enum gsl_integration_qawo_enum sine)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function changes the parameters :data:`omega`, :data:`L` and :data:`sine`
|
|
Packit |
67cb25 |
of the existing workspace :data:`t`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t, double L)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allows the length parameter :data:`L` of the workspace
|
|
Packit |
67cb25 |
:data:`t` to be changed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_qawo_table_free (gsl_integration_qawo_table * t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees all the memory associated with the workspace :data:`t`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qawo (gsl_function * f, const double a, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace * workspace, gsl_integration_qawo_table * wf, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function uses an adaptive algorithm to compute the integral of
|
|
Packit |
67cb25 |
:math:`f` over :math:`(a,b)` with the weight function
|
|
Packit |
67cb25 |
:math:`\sin(\omega x)` or :math:`\cos(\omega x)` defined
|
|
Packit |
67cb25 |
by the table :data:`wf`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
I = \int_a^b dx f(x)
|
|
Packit |
67cb25 |
\left\{
|
|
Packit |
67cb25 |
\begin{array}{c}
|
|
Packit |
67cb25 |
\sin{(\omega x)} \\
|
|
Packit |
67cb25 |
\cos{(\omega x)}
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right\}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
| I = \int_a^b dx f(x) sin(omega x)
|
|
Packit |
67cb25 |
| I = \int_a^b dx f(x) cos(omega x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The results are extrapolated using the epsilon-algorithm to accelerate
|
|
Packit |
67cb25 |
the convergence of the integral. The function returns the final
|
|
Packit |
67cb25 |
approximation from the extrapolation, :data:`result`, and an estimate of
|
|
Packit |
67cb25 |
the absolute error, :data:`abserr`. The subintervals and their results are
|
|
Packit |
67cb25 |
stored in the memory provided by :data:`workspace`. The maximum number of
|
|
Packit |
67cb25 |
subintervals is given by :data:`limit`, which may not exceed the allocated
|
|
Packit |
67cb25 |
size of the workspace.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Those subintervals with "large" widths :math:`d` where :math:`d\omega > 4` are
|
|
Packit |
67cb25 |
computed using a 25-point Clenshaw-Curtis integration rule, which handles the
|
|
Packit |
67cb25 |
oscillatory behavior. Subintervals with a "small" widths where
|
|
Packit |
67cb25 |
:math:`d\omega < 4` are computed using a 15-point Gauss-Kronrod integration.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QAWF adaptive integration for Fourier integrals
|
|
Packit |
67cb25 |
===============================================
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: QAWF quadrature algorithm
|
|
Packit |
67cb25 |
single: Fourier integrals, numerical
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_qawf (gsl_function * f, const double a, const double epsabs, const size_t limit, gsl_integration_workspace * workspace, gsl_integration_workspace * cycle_workspace, gsl_integration_qawo_table * wf, double * result, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function attempts to compute a Fourier integral of the function
|
|
Packit |
67cb25 |
:data:`f` over the semi-infinite interval :math:`[a,+\infty)`
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
I = \int_a^{+\infty} dx f(x)
|
|
Packit |
67cb25 |
\left\{
|
|
Packit |
67cb25 |
\begin{array}{c}
|
|
Packit |
67cb25 |
\sin{(\omega x)} \\
|
|
Packit |
67cb25 |
\cos{(\omega x)}
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right\}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
I = \int_a^{+\infty} dx f(x) sin(omega x)
|
|
Packit |
67cb25 |
I = \int_a^{+\infty} dx f(x) cos(omega x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The parameter :math:`\omega` and choice of :math:`\sin` or :math:`\cos` is
|
|
Packit |
67cb25 |
taken from the table :data:`wf` (the length :data:`L` can take any value,
|
|
Packit |
67cb25 |
since it is overridden by this function to a value appropriate for the
|
|
Packit |
67cb25 |
Fourier integration). The integral is computed using the QAWO algorithm
|
|
Packit |
67cb25 |
over each of the subintervals,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
C_1 &= [a,a+c] \\
|
|
Packit |
67cb25 |
C_2 &= [a+c,a+2c] \\
|
|
Packit |
67cb25 |
\dots &= \dots \\
|
|
Packit |
67cb25 |
C_k &= [a+(k-1)c,a+kc]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
C_1 = [a, a + c]
|
|
Packit |
67cb25 |
C_2 = [a + c, a + 2 c]
|
|
Packit |
67cb25 |
... = ...
|
|
Packit |
67cb25 |
C_k = [a + (k-1) c, a + k c]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
:math:`c = (2 floor(|\omega|) + 1) \pi/|\omega|`. The width :math:`c` is
|
|
Packit |
67cb25 |
chosen to cover an odd number of periods so that the contributions from
|
|
Packit |
67cb25 |
the intervals alternate in sign and are monotonically decreasing when
|
|
Packit |
67cb25 |
:data:`f` is positive and monotonically decreasing. The sum of this
|
|
Packit |
67cb25 |
sequence of contributions is accelerated using the epsilon-algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function works to an overall absolute tolerance of
|
|
Packit |
67cb25 |
:data:`abserr`. The following strategy is used: on each interval
|
|
Packit |
67cb25 |
:math:`C_k` the algorithm tries to achieve the tolerance
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: TOL_k = u_k abserr
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
:math:`u_k = (1 - p)p^{k-1}` and :math:`p = 9/10`.
|
|
Packit |
67cb25 |
The sum of the geometric series of contributions from each interval
|
|
Packit |
67cb25 |
gives an overall tolerance of :data:`abserr`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If the integration of a subinterval leads to difficulties then the
|
|
Packit |
67cb25 |
accuracy requirement for subsequent intervals is relaxed,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: TOL_k = u_k \max{(abserr, \max_{i
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: TOL_k = u_k max(abserr, max_{i
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`E_k` is the estimated error on the interval :math:`C_k`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The subintervals and their results are stored in the memory provided by
|
|
Packit |
67cb25 |
:data:`workspace`. The maximum number of subintervals is given by
|
|
Packit |
67cb25 |
:data:`limit`, which may not exceed the allocated size of the workspace.
|
|
Packit |
67cb25 |
The integration over each subinterval uses the memory provided by
|
|
Packit |
67cb25 |
:data:`cycle_workspace` as workspace for the QAWO algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: cquad, doubly-adaptive integration
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
CQUAD doubly-adaptive integration
|
|
Packit |
67cb25 |
=================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|cquad| is a new doubly-adaptive general-purpose quadrature
|
|
Packit |
67cb25 |
routine which can handle most types of singularities,
|
|
Packit |
67cb25 |
non-numerical function values such as :code:`Inf` or :code:`NaN`,
|
|
Packit |
67cb25 |
as well as some divergent integrals. It generally requires more
|
|
Packit |
67cb25 |
function evaluations than the integration routines in
|
|
Packit |
67cb25 |
|quadpack|, yet fails less often for difficult integrands.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The underlying algorithm uses a doubly-adaptive scheme in which
|
|
Packit |
67cb25 |
Clenshaw-Curtis quadrature rules of increasing degree are used
|
|
Packit |
67cb25 |
to compute the integral in each interval. The :math:`L_2`-norm of
|
|
Packit |
67cb25 |
the difference between the underlying interpolatory polynomials
|
|
Packit |
67cb25 |
of two successive rules is used as an error estimate. The
|
|
Packit |
67cb25 |
interval is subdivided if the difference between two successive
|
|
Packit |
67cb25 |
rules is too large or a rule of maximum degree has been reached.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: gsl_integration_cquad_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_cquad_workspace * gsl_integration_cquad_workspace_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace sufficient to hold the data for
|
|
Packit |
67cb25 |
:data:`n` intervals. The number :data:`n` is not the maximum number of
|
|
Packit |
67cb25 |
intervals that will be evaluated. If the workspace is full, intervals
|
|
Packit |
67cb25 |
with smaller error estimates will be discarded. A minimum of 3
|
|
Packit |
67cb25 |
intervals is required and for most functions, a workspace of size 100
|
|
Packit |
67cb25 |
is sufficient.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_cquad_workspace_free (gsl_integration_cquad_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_cquad (const gsl_function * f, double a, double b, double epsabs, double epsrel, gsl_integration_cquad_workspace * workspace, double * result, double * abserr, size_t * nevals)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the integral of :math:`f` over :math:`(a,b)`
|
|
Packit |
67cb25 |
within the desired absolute and relative error limits, :data:`epsabs`
|
|
Packit |
67cb25 |
and :data:`epsrel` using the |cquad| algorithm. The function returns
|
|
Packit |
67cb25 |
the final approximation, :data:`result`, an estimate of the absolute
|
|
Packit |
67cb25 |
error, :data:`abserr`, and the number of function evaluations required,
|
|
Packit |
67cb25 |
:data:`nevals`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The |cquad| algorithm divides the integration region into
|
|
Packit |
67cb25 |
subintervals, and in each iteration, the subinterval with the largest
|
|
Packit |
67cb25 |
estimated error is processed. The algorithm uses Clenshaw-Curits
|
|
Packit |
67cb25 |
quadrature rules of degree 4, 8, 16 and 32 over 5, 9, 17 and 33 nodes
|
|
Packit |
67cb25 |
respectively. Each interval is initialized with the lowest-degree
|
|
Packit |
67cb25 |
rule. When an interval is processed, the next-higher degree rule is
|
|
Packit |
67cb25 |
evaluated and an error estimate is computed based on the
|
|
Packit |
67cb25 |
:math:`L_2`-norm of the difference between the underlying interpolating
|
|
Packit |
67cb25 |
polynomials of both rules. If the highest-degree rule has already been
|
|
Packit |
67cb25 |
used, or the interpolatory polynomials differ significantly, the
|
|
Packit |
67cb25 |
interval is bisected.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The subintervals and their results are stored in the memory
|
|
Packit |
67cb25 |
provided by :data:`workspace`. If the error estimate or the number of
|
|
Packit |
67cb25 |
function evaluations is not needed, the pointers :data:`abserr` and :data:`nevals`
|
|
Packit |
67cb25 |
can be set to :code:`NULL`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Romberg integration
|
|
Packit |
67cb25 |
===================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The Romberg integration method estimates the definite integral
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: I = \int_a^b f(x) dx
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
by applying Richardson extrapolation on the trapezoidal rule, using
|
|
Packit |
67cb25 |
equally spaced points with spacing
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: h_k = (b - a) 2^{-k}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for :math:`k = 1, \dots, n`. For each :math:`k`, Richardson extrapolation
|
|
Packit |
67cb25 |
is used :math:`k-1` times on previous approximations to improve the order
|
|
Packit |
67cb25 |
of accuracy as much as possible. Romberg integration typically works
|
|
Packit |
67cb25 |
well (and converges quickly) for smooth integrands with no singularities in
|
|
Packit |
67cb25 |
the interval or at the end points.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_romberg_workspace * gsl_integration_romberg_alloc(const size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for Romberg integration, specifying
|
|
Packit |
67cb25 |
a maximum of :math:`n` iterations, or divisions of the interval. Since
|
|
Packit |
67cb25 |
the number of divisions is :math:`2^n + 1`, :math:`n` can be kept relatively
|
|
Packit |
67cb25 |
small (i.e. :math:`10` or :math:`20`). It is capped at a maximum value of
|
|
Packit |
67cb25 |
:math:`30` to prevent overflow. The size of the workspace is :math:`O(2n)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_romberg_free(gsl_integration_romberg_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_romberg(const gsl_function * f, const double a, const double b, const double epsabs, const double epsrel, double * result, size_t * neval, gsl_integration_romberg_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function integrates :math:`f(x)`, specified by :data:`f`, from :data:`a` to
|
|
Packit |
67cb25 |
:data:`b`, storing the answer in :data:`result`. At each step in the iteration,
|
|
Packit |
67cb25 |
convergence is tested by checking:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: | I_k - I_{k-1} | \le \textrm{max} \left( epsabs, epsrel \times |I_k| \right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`I_k` is the current approximation and :math:`I_{k-1}` is the approximation
|
|
Packit |
67cb25 |
of the previous iteration. If the method does not converge within the previously
|
|
Packit |
67cb25 |
specified :math:`n` iterations, the function stores the best current estimate in
|
|
Packit |
67cb25 |
:data:`result` and returns :macro:`GSL_EMAXITER`. If the method converges, the function
|
|
Packit |
67cb25 |
returns :macro:`GSL_SUCCESS`. The total number of function evaluations is returned
|
|
Packit |
67cb25 |
in :data:`neval`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Gauss-Legendre integration
|
|
Packit |
67cb25 |
==========================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The fixed-order Gauss-Legendre integration routines are provided for fast
|
|
Packit |
67cb25 |
integration of smooth functions with known polynomial order. The
|
|
Packit |
67cb25 |
:math:`n`-point Gauss-Legendre rule is exact for polynomials of order
|
|
Packit |
67cb25 |
:math:`2n - 1` or less. For example, these rules are useful when integrating
|
|
Packit |
67cb25 |
basis functions to form mass matrices for the Galerkin method. Unlike other
|
|
Packit |
67cb25 |
numerical integration routines within the library, these routines do not accept
|
|
Packit |
67cb25 |
absolute or relative error bounds.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: gsl_integration_glfixed_table
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_glfixed_table * gsl_integration_glfixed_table_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function determines the Gauss-Legendre abscissae and weights necessary for
|
|
Packit |
67cb25 |
an :math:`n`-point fixed order integration scheme. If possible, high precision
|
|
Packit |
67cb25 |
precomputed coefficients are used. If precomputed weights are not available,
|
|
Packit |
67cb25 |
lower precision coefficients are computed on the fly.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_integration_glfixed (const gsl_function * f, double a, double b, const gsl_integration_glfixed_table * t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Gauss-Legendre integration rule
|
|
Packit |
67cb25 |
contained in table :data:`t` and returns the result.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_glfixed_point (double a, double b, size_t i, double * xi, double * wi, const gsl_integration_glfixed_table * t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For :data:`i` in :math:`[0, \dots, n - 1]`, this function obtains the
|
|
Packit |
67cb25 |
:data:`i`-th Gauss-Legendre point :data:`xi` and weight :data:`wi` on the interval
|
|
Packit |
67cb25 |
[:data:`a`, :data:`b`]. The points and weights are ordered by increasing point
|
|
Packit |
67cb25 |
value. A function :math:`f` may be integrated on [:data:`a`, :data:`b`] by summing
|
|
Packit |
67cb25 |
:math:`wi * f(xi)` over :data:`i`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: gsl_integration_glfixed_table
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_glfixed_table_free (gsl_integration_glfixed_table * t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the table :data:`t`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Fixed point quadratures
|
|
Packit |
67cb25 |
=======================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: interpolating quadrature
|
|
Packit |
67cb25 |
single: quadrature, interpolating
|
|
Packit |
67cb25 |
single: quadrature, fixed point
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The routines in this section approximate an integral by the sum
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \int_a^b w(x) f(x) dx = \sum_{i=1}^n w_i f(x_i)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`f(x)` is the function to be integrated and :math:`w(x)` is
|
|
Packit |
67cb25 |
a weighting function. The :math:`n` weights :math:`w_i` and nodes :math:`x_i` are carefully chosen
|
|
Packit |
67cb25 |
so that the result is exact when :math:`f(x)` is a polynomial of degree :math:`2n - 1`
|
|
Packit |
67cb25 |
or less. Once the user chooses the order :math:`n` and weighting function :math:`w(x)`, the
|
|
Packit |
67cb25 |
weights :math:`w_i` and nodes :math:`x_i` can be precomputed and used to efficiently evaluate
|
|
Packit |
67cb25 |
integrals for any number of functions :math:`f(x)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This method works best when :math:`f(x)` is well approximated by a polynomial on the interval
|
|
Packit |
67cb25 |
:math:`(a,b)`, and so is not suitable for functions with singularities.
|
|
Packit |
67cb25 |
Since the user specifies ahead of time how many quadrature nodes will be used, these
|
|
Packit |
67cb25 |
routines do not accept absolute or relative error bounds. The table below lists
|
|
Packit |
67cb25 |
the weighting functions currently supported.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _tab_fixed-quadratures:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
================ ======================== =========================================== =======================================================
|
|
Packit |
67cb25 |
Name Interval Weighting function :math:`w(x)` Constraints
|
|
Packit |
67cb25 |
================ ======================== =========================================== =======================================================
|
|
Packit |
67cb25 |
Legendre :math:`(a,b)` :math:`1` :math:`b > a`
|
|
Packit |
67cb25 |
Chebyshev Type 1 :math:`(a,b)` :math:`1 / \sqrt{(b - x) (x - a)}` :math:`b > a`
|
|
Packit |
67cb25 |
Gegenbauer :math:`(a,b)` :math:`((b - x) (x - a))^{\alpha}` :math:`\alpha > -1, b > a`
|
|
Packit |
67cb25 |
Jacobi :math:`(a,b)` :math:`(b - x)^{\alpha} (x - a)^{\beta}` :math:`\alpha,\beta > -1, b > a`
|
|
Packit |
67cb25 |
Laguerre :math:`(a,\infty)` :math:`(x-a)^{\alpha} \exp{( -b (x - a) )}` :math:`\alpha > -1, b > 0`
|
|
Packit |
67cb25 |
Hermite :math:`(-\infty,\infty)` :math:`|x-a|^{\alpha} \exp{( -b (x-a)^2 )}` :math:`\alpha > -1, b > 0`
|
|
Packit |
67cb25 |
Exponential :math:`(a,b)` :math:`|x - (a + b)/2|^{\alpha}` :math:`\alpha > -1, b > a`
|
|
Packit |
67cb25 |
Rational :math:`(a,\infty)` :math:`(x - a)^{\alpha} (x + b)^{\beta}` :math:`\alpha > -1, \alpha + \beta + 2n < 0, a + b > 0`
|
|
Packit |
67cb25 |
Chebyshev Type 2 :math:`(a,b)` :math:`\sqrt{(b - x) (x - a)}` :math:`b > a`
|
|
Packit |
67cb25 |
================ ======================== =========================================== =======================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The fixed point quadrature routines use the following workspace to store the nodes and weights,
|
|
Packit |
67cb25 |
as well as additional variables for intermediate calculations:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_integration_fixed_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This workspace is used for fixed point quadrature rules and looks like this::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n; /* number of nodes/weights */
|
|
Packit |
67cb25 |
double *weights; /* quadrature weights */
|
|
Packit |
67cb25 |
double *x; /* quadrature nodes */
|
|
Packit |
67cb25 |
double *diag; /* diagonal of Jacobi matrix */
|
|
Packit |
67cb25 |
double *subdiag; /* subdiagonal of Jacobi matrix */
|
|
Packit |
67cb25 |
const gsl_integration_fixed_type * type;
|
|
Packit |
67cb25 |
} gsl_integration_fixed_workspace;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_integration_fixed_workspace * gsl_integration_fixed_alloc(const gsl_integration_fixed_type * T, const size_t n, const double a, const double b, const double alpha, const double beta)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for computing integrals with interpolating quadratures using :data:`n`
|
|
Packit |
67cb25 |
quadrature nodes. The parameters :data:`a`, :data:`b`, :data:`alpha`, and :data:`beta` specify the integration
|
|
Packit |
67cb25 |
interval and/or weighting function for the various quadrature types. See the :ref:`table <tab_fixed-quadratures>` above
|
|
Packit |
67cb25 |
for constraints on these parameters. The size of the workspace is :math:`O(4n)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_integration_fixed_type
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The type of quadrature used is specified by :data:`T` which can be set to the following choices:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_legendre
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Legendre quadrature integration. The parameters :data:`alpha` and
|
|
Packit |
67cb25 |
:data:`beta` are ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_chebyshev
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Chebyshev type 1 quadrature integration. The parameters :data:`alpha` and
|
|
Packit |
67cb25 |
:data:`beta` are ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_gegenbauer
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Gegenbauer quadrature integration. The parameter :data:`beta` is ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_jacobi
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Jacobi quadrature integration.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_laguerre
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Laguerre quadrature integration. The parameter :data:`beta` is ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_hermite
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Hermite quadrature integration. The parameter :data:`beta` is ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_exponential
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies exponential quadrature integration. The parameter :data:`beta` is ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_rational
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies rational quadrature integration.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_integration_fixed_chebyshev2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies Chebyshev type 2 quadrature integration. The parameters :data:`alpha` and
|
|
Packit |
67cb25 |
:data:`beta` are ignored for this type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_integration_fixed_free(gsl_integration_fixed_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory assocated with the workspace :data:`w`
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: size_t gsl_integration_fixed_n(const gsl_integration_fixed_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns the number of quadrature nodes and weights.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double * gsl_integration_fixed_nodes(const gsl_integration_fixed_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a pointer to an array of size :data:`n` containing the quadrature nodes :math:`x_i`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double * gsl_integration_fixed_weights(const gsl_integration_fixed_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a pointer to an array of size :data:`n` containing the quadrature weights :math:`w_i`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_integration_fixed(const gsl_function * func, double * result, const gsl_integration_fixed_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function integrates the function :math:`f(x)` provided in :data:`func` using previously
|
|
Packit |
67cb25 |
computed fixed quadrature rules. The integral is approximated as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \sum_{i=1}^n w_i f(x_i)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`w_i` are the quadrature weights and :math:`x_i` are the quadrature nodes computed
|
|
Packit |
67cb25 |
previously by :func:`gsl_integration_fixed_alloc`. The sum is stored in :data:`result` on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Error codes
|
|
Packit |
67cb25 |
===========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In addition to the standard error codes for invalid arguments the
|
|
Packit |
67cb25 |
functions can return the following values,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
===================== ============================================================================================================
|
|
Packit |
67cb25 |
:macro:`GSL_EMAXITER` the maximum number of subdivisions was exceeded.
|
|
Packit |
67cb25 |
:macro:`GSL_EROUND` cannot reach tolerance because of roundoff error, or roundoff error was detected in the extrapolation table.
|
|
Packit |
67cb25 |
:macro:`GSL_ESING` a non-integrable singularity or other bad integrand behavior was found in the integration interval.
|
|
Packit |
67cb25 |
:macro:`GSL_EDIVERGE` the integral is divergent, or too slowly convergent to be integrated numerically.
|
|
Packit |
67cb25 |
:macro:`GSL_EDOM` error in the values of the input arguments
|
|
Packit |
67cb25 |
===================== ============================================================================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Examples
|
|
Packit |
67cb25 |
========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Adaptive integration example
|
|
Packit |
67cb25 |
----------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The integrator :code:`QAGS` will handle a large class of definite
|
|
Packit |
67cb25 |
integrals. For example, consider the following integral, which has an
|
|
Packit |
67cb25 |
algebraic-logarithmic singularity at the origin,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \int_0^1 x^{-1/2} \log(x) dx = -4
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The program below computes this integral to a relative accuracy bound of
|
|
Packit |
67cb25 |
:code:`1e-7`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/integration.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The results below show that the desired accuracy is achieved after 8
|
|
Packit |
67cb25 |
subdivisions.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/integration.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In fact, the extrapolation procedure used by :code:`QAGS` produces an
|
|
Packit |
67cb25 |
accuracy of almost twice as many digits. The error estimate returned by
|
|
Packit |
67cb25 |
the extrapolation procedure is larger than the actual error, giving a
|
|
Packit |
67cb25 |
margin of safety of one order of magnitude.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Fixed-point quadrature example
|
|
Packit |
67cb25 |
------------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In this example, we use a fixed-point quadrature rule to integrate the
|
|
Packit |
67cb25 |
integral
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
\int_{-\infty}^{\infty} e^{-x^2} \left( x^m + 1 \right) dx =
|
|
Packit |
67cb25 |
\left\{
|
|
Packit |
67cb25 |
\begin{array}{cc}
|
|
Packit |
67cb25 |
\sqrt{\pi} + \Gamma{\left( \frac{m+1}{2} \right)}, & m \textrm{ even} \\
|
|
Packit |
67cb25 |
\sqrt{\pi}, & m \textrm{ odd}
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for integer :math:`m`. Consulting our :ref:`table <tab_fixed-quadratures>` of fixed point quadratures,
|
|
Packit |
67cb25 |
we see that this integral can be evaluated with a Hermite quadrature rule,
|
|
Packit |
67cb25 |
setting :math:`\alpha = 0, a = 0, b = 1`. Since we are integrating a polynomial
|
|
Packit |
67cb25 |
of degree :math:`m`, we need to choose the number of nodes :math:`n \ge (m+1)/2`
|
|
Packit |
67cb25 |
to achieve the best results.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
First we will try integrating for :math:`m = 10, n = 5`, which does not satisfy
|
|
Packit |
67cb25 |
our criteria above::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
$ ./integration2 10 5
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The output is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/integration2a.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
So, we find a large error. Now we try integrating for :math:`m = 10, n = 6` which
|
|
Packit |
67cb25 |
does satisfy the criteria above::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
$ ./integration2 10 6
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The output is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/integration2b.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The program is given below.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/integration2.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
References and Further Reading
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following book is the definitive reference for |quadpack|, and was
|
|
Packit |
67cb25 |
written by the original authors. It provides descriptions of the
|
|
Packit |
67cb25 |
algorithms, program listings, test programs and examples. It also
|
|
Packit |
67cb25 |
includes useful advice on numerical integration and many references to
|
|
Packit |
67cb25 |
the numerical integration literature used in developing |quadpack|.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* R. Piessens, E. de Doncker-Kapenga, C.W. Ueberhuber, D.K. Kahaner.
|
|
Packit |
67cb25 |
|quadpack| A subroutine package for automatic integration
|
|
Packit |
67cb25 |
Springer Verlag, 1983.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The |cquad| integration algorithm is described in the following paper:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. Gonnet, "Increasing the Reliability of Adaptive Quadrature Using
|
|
Packit |
67cb25 |
Explicit Interpolants", ACM Transactions on Mathematical
|
|
Packit |
67cb25 |
Software, Volume 37 (2010), Issue 3, Article 26.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The fixed-point quadrature routines are based on IQPACK, described in the
|
|
Packit |
67cb25 |
following papers:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* S. Elhay, J. Kautsky, Algorithm 655: IQPACK, FORTRAN Subroutines for the
|
|
Packit |
67cb25 |
Weights of Interpolatory Quadrature, ACM Transactions on Mathematical Software,
|
|
Packit |
67cb25 |
Volume 13, Number 4, December 1987, pages 399-415.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* J. Kautsky, S. Elhay, Calculation of the Weights of Interpolatory Quadratures,
|
|
Packit |
67cb25 |
Numerische Mathematik, Volume 40, Number 3, October 1982, pages 407-422.
|