% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Numerical integration}

\section{Introduction}
{\em Numerical integration} constitutes a broad family of algorithms to
compute a numerical approximation to a definite (Riemann) integral.

Generally, the integral is approximated by a weighted sum of function
values within the domain of integration,
	\begin{equation}\label{q:e:1}
\int_a^b f(x)dx \approx \sum_{i=1}^{n} w_i f(x_i) \:.
	\end{equation}
Expression~(\ref{q:e:1}) is often referred to as {\em
quadrature}\index{quadrature} ({\em cubature}\index{cubature} for
multidimensional integrals) or {\em rule}\index{rule}.	The abscissas
$x_i$ (also called {\em nodes}\index{nodes}) and the weights $w_i$
of a quadrature are usually optimized---using one of a large number of
different strategies---to suit a particular class of integration problems.

The best quadrature algorithm for a given problem depends on
several factors, in particular on the integrand.  Different classes
of integrands generally require different quadratures for the most
effective calculation.

A popular numerical integration library is QUADPACK~\cite{quadpack}.
It includes general purpose routines---like QAGS, basen on an adaptive
Gauss–Kronrod quadrature with acceleration---as well as a number of
specialized routines.  The GNU scientific library~\cite{gsl} (GSL)
implements most of the QUADPACK routines and in addition includes a
modern general-purpose adaptive routine CQUAD based on Clenshaw-Curtis
quadratures~\cite{gonnet}.

In the following we shall consider some of the popular numerical
integration algorithms.

\section{Rectangle and trapezium rules}
In mathematics, the {\em Reimann integral}\index{Riemann integral}
is generally defined in terms of {\em Riemann sums}\index{Riemann
sum}~\cite{riemann}. If the integration interval $[a,b]$ is partitioned
into $n$ subintervals,
	\begin{equation}\label{q:r:e:p}
a=t_0<t_1<t_2<\dots<t_n=b \;.
	\end{equation}
the Riemann sum is defined as
	\begin{equation}\label{q:e:rs}
\sum_{i=1}^n f(x_i)\Delta x_i \;,
	\end{equation}
where $x_i\in[t_{i-1},t_i]$ and $\Delta x_i=t_{i}-t_{i-1}$.  Geometrically
a Riemann sum can be interpreted as the area of a collection of adjucent
rectangles with widths $\Delta x_i$ and heights $f(x_i)$.

The Rieman integral is defined as the limit of a Riemann sum as the
{\em mesh}\index{mesh}---the length of the largest subinterval---of the
partition approaches zero. Specifically, the number denoted as
	\begin{equation}
\int_a^bf(x)dx
	\end{equation}
is called the Riemann integral, if for any $\epsilon>0$ there exists
$\delta>0$ such that for any partition~(\ref{q:r:e:p}) with
$\max\Delta x_i<\delta$ we have
	\begin{equation}
\left| \sum_{i=1}^n f(x_i)\Delta x_i - \int_a^bf(x)dx \right| <
\epsilon  \;.
	\end{equation}
A definite integral can be interpreted as the net signed area bounded
by the graph of the integrand.

Now, the $n$-point {\em rectangle quadrature}\index{rectangle rule}
is simply the Riemann sum~(\ref{q:e:rs}),
	\begin{equation}\label{q:e:rec}
\int_a^bf(x)dx \approx \sum_{i=1}^n f(x_i)\Delta x_i \:,
	\end{equation}
where the node $x_i$ is often (but not always) taken at the middle
of the corresponding subinterval, $x_i=t_{i-1}+\frac12\Delta x_i$,
and the subintervals are often (but not always) chosen equal, $\Delta
x_i=(b-a)/n$.  Geometrically the $n$-point rectangle rule is an
approximation to the integral given by the area of a collection of $n$
adjucent equal rectangles whose heights are determined by the values of the
function (at the middle of the rectangle).

An $n$-point {\em trapezium rule}\index{trapezium rule} uses instead
a collection of trapezia fitted under the graph,
	\begin{equation}
\int_a^bf(x)dx \approx
\sum_{i=1}^n \frac{f(t_{i-1})+f(t_i)}{2}\Delta x_i \;.
	\end{equation}
Importantly, the trapezium rule is the average of two Riemann sums,
	\begin{equation}\label{q:e:2s}
\sum_{i=1}^n \frac{f(t_{i-1})+f(t_i)}{2}\Delta x_i =
\frac12 \sum_{i=1}^n f(t_{i-1}) \Delta x_i
+ \frac12 \sum_{i=1}^n f(t_i) \Delta x_i \;.
	\end{equation}

Rectangle and trapezium quadratures both have the important feature of
closely following the very mathematical definition of the integral as
the limit of the Riemann sums.  Therefore---disregarding the round-off
errors---these two rules cannot fail if the integral exists.

For certain partitions of the interval the rectangle and trapezium
rules coincide.  For example, for the nodes
	\begin{equation}\label{q:eq:pa}
x_i = a+(b-a)\frac{i-\frac12}{n} \;,\;\;\; i=1,\dots,n \;
	\end{equation}
both rules give the same quadrature with equal weights, $w_i=(b-a)/n$,
	\begin{equation}\label{q:eq:rt}
\int_a^b f(x)dx \approx \frac{b-a}{n}\sum_{i=1}^n
f\left(a+(b-a)\frac{i-\frac12}{n}\right) \;.
	\end{equation}

Rectangle and trapezium quadratures are rarely used on their own---because
of the slow convergence---but they often serve as the basis for more
advanced quadratures, for example adaptive quadratures and variable
transformation quadratures considered below.

\section{Quadratures with regularly spaced abscissas}

A quadrature~(\ref{q:e:1}) with $n$ predefined nodes $x_i$ has $n$ free
parameters: the weights $w_i$.	A set of $n$ parameters can generally
be tuned to satisfy $n$ conditions. The archetypal set of conditions
in quadratures is that the quadrature integrates exactly a set of $n$
functions,
	\begin{equation}\label{q:eq:fis}
\left\{\phi_1(x),\dots,\phi_n(x)\right\} \;.
	\end{equation}
This leads to a set of $n$ equations,
	\begin{equation}\label{q:eq:w}
\sum_{i=1}^n w_i \phi_k(x_i)=I_k \;\Big|_{k=1,\dots,n}\;,
	\end{equation}
where the integrals
	\begin{equation}
I_k\doteq\int_a^b\phi_k(x)dx
	\end{equation}
are assumed to be known. Equations~(\ref{q:eq:w}) are linear in $w_i$
and can be easily solved.

Since integration is a linear operation, the quadrature will then also
integrate exactly any linear combination of functions~(\ref{q:eq:fis}).

A popular choice for predefined nodes is a {\em closed} set---that is,
including the end-points of the interval---of evenly spaced abscissas,
	\begin{equation}
x_i=a+\frac{i-1}{n-1}(b-a) \;\Big|_{i=1,\dots,n}\;.
	\end{equation}
However, in practice it often happens that the integrand has an integrable
sinfularity at one or both ends of the interval.  In this case one can
choose an {\em open set} of equidistant nodes,
	\begin{equation}
x_i=a+\frac{i-\frac12}{n}(b-a) \;\Big|_{i=1,\dots,n}\;.
	\end{equation}

The set of functions to be integrated exactly is generally chosen to
suite the properties of the integrands at hand: the integrands must be
well represented by linear combinations of the chosen functions.

\subsection{Classical quadratures}

Suppose the integrand can be well represented by the first few terms of
its Taylor series,
	\begin{equation}\label{q:eq:tay}
f(x)=\sum_{k=0}^\infty {f^{(k)}(a)\over k!}(x-a)^k \:,
	\end{equation}
where $f^{(k)}$ is the $k$-th derivative of the integrand.  This is often
the case for analytic---that is, infinitely differentiable---functions.
For such integrands one can obviously choose polynomials
	\begin{equation}\label{q:eq:ps}
\{1,x,x^2,\dots,x^{n-1}\}
	\end{equation}
as the set of functions to be integrated exactly.

This leads to the so called {\em classical quadratures}\index{classical
quadrature}: quadratures with regularly spaced abscissas and {\em
polynomials} as exactly integrable functions.

An $n$-point classical quadrature integrates exactly the first $n$ terms
of the function's Taylor expansion~(\ref{q:eq:tay}).  The $x^n$ order
term will not be integrated exactly and will lead to an error of the
quadrature.  Thus the error $E_n$ of the $n$-point classical quadrature
is on the order of the integral of the $x^n$ term in~(\ref{q:eq:tay}),
	\begin{equation}
E_n \approx \int_a^b {f^{(n)}(a)\over n!}(x-a)^n dx ={f^{(n)}(a)\over
(n+1)!}h^{n+1} \propto h^{n+1}\:,
	\end{equation}
where $h=b-a$ is the length of the integration interval.  A quadrature
with the error of the order $h^{n+1}$ is often called a {\em degree}-$n$
quadrature\index{degree-$n$ quadrature}.

If the integrand is smooth enough and the length $h$ is small enough
a classical quadrature with not so large $n$ can provide a good
approximation for the integral.  However, for large $n$ the weights of
classical quadratures tend to have alternating signs, which leads to large
round-off errors, which in turn negates the potentially higher accuracy of
the quadrature. Again, if the integrand violates the assumption of Taylor
expansion---for example by having an integrable singularity inside the
integration interval---the higher order quadratures may perform poorly.

Classical quadratures are mostly of historical interest nowadays.
Alternative methods---such as quadratures with optimized abscissas,
adaptive, and variable transformation quadratures---are more stable and
accurate and are normally preferred to classical quadratures.

\begin{table}
\caption{Maxima script to calculate analytically the weights of an
$n$-point classical quadrature with predefined abscissas in the interval
$[0,1]$.}
\label{q:t:nc}
{\small
\begin{verbatim}
n: 8; xs: makelist((i-1)/(n-1),i,1,n); /* nodes: adapt to your needs */
ws: makelist(concat(w,i),i,1,n);
ps: append([1],makelist(x^i,i,1,n-1)); /* polynomials */
fs: makelist(buildq([i:i,ps:ps],lambda([x],ps[i])),i,1,n);
integ01: lambda([f],integrate(f(x),x,0,1));
Is: maplist(integ01,fs); /* calculate the integrals */
eq: lambda([f],lreduce("+",maplist(f,xs)*ws));
eqs: maplist(eq,fs)-Is; /* build equations */
solve(eqs,ws); /* solve for the weights */
\end{verbatim}}
\hrule
\end{table}

Classical quadratures with equally spaced abscissas---both closed
and open sets---are generally referred to as {\em Newton-Cotes
quadratures}\index{Newton-Cotes quadrature}.  An interested reader can
generate Newton–-Cotes quadratures of any degree $n$ using the Maxima
script in Table~(\ref{q:t:nc}).

\section{Quadratures with optimized abscissas}

In quadratures with optimized abscissas not only the weights $w_i$ but also
the abscissas $x_i$ are chosen optimally.  The number of free parameters
is thus $2n$ and one can choose a set of $2n$ functions,
	\begin{equation}
\{\phi_1(x),\dots,\phi_{2n}(x)\} \;,
	\end{equation}
to be integrated exactly. This gives a system of $2n$ equations, linear
in $w_i$ and non-linear in $x_i$,
	\begin{equation}\label{q:eq:wx}
\sum_{i=1}^n w_i \phi_k(x_i)=I_k \:\Big|_{k=1,\dots,2n}\;,
	\end{equation}
where again
	\begin{equation}
I_k\doteq\int_a^b \phi_k(x)dx \:.
	\end{equation}
The weights and abscissas of the quadrature can be determined by solving
this system of equations\footnote{Here is, for example, an $n=2$
quadrature with optimized abscissas,
	\begin{equation}
\int_{-1}^{1}f(x)dx\approx
f\left(-\sqrt{\tfrac{1}{3}}\right)+f\left(+\sqrt{\tfrac{1}{3}}\right)\:.
	\end{equation}}.

Although quadratures with optimized abcissas are generally of much
higher order, $2n-1$ compared to $n-1$ for non-optimal abscissas, the
optimal points generally can not be reused at the next iteration in an
adaptive algorithm.

\subsection{Gauss quadratures}

Gauss quadratures deal with a slightly more general form of integrals,
	\begin{equation}
\int_a^b \omega(x) f(x) dx \;,
	\end{equation}
where $\omega(x)$ is a positive weight function. For $\omega(x)=1$ the
problem is the same as considered above.  Popular choices of the
weight function include $\omega(x)=(1-x^2)^{\pm 1/2}$, $\exp(-x)$,
$\exp(-x^2)$ and others.  The idea is to represent the integrand as a
product $\omega(x) f(x)$ such that all the difficulties go into the
weight function $\omega(x)$ while the remaining factor $f(x)$ is smooth
and well represented by polynomials.

An $N$-point {\em Gauss quadrature} is a quadrature with optimized
abcissas,
	\begin{equation}
\int_a^b \omega(x) f(x) dx \approx \sum_{i=1}^N w_i f(x_i) \;,
	\end{equation}
which integrates exactly a set of $2N$ polynomials of the orders
$1,\dots,2N-1$ with the given weight $\omega(x)$.

\subsubsection{Fundamental theorem}
There is a theorem stating that there exists a set of polynomials
$p_n(x)$, orthogonal on the interval $[a,b]$ with the weight function
$\omega(x)$,
	\begin{equation}
\int_a^b\omega(x) p_n(x)p_k(x) \propto \delta_{nk}\;.
	\end{equation}

Now, one can prove that the optimal nodes for the $N$-point Gauss
quadrature are the roots of the polynomial $p_N(x)$,
	\begin{equation}
p_N(x_i)=0 \;.
	\end{equation}

The idea behind the proof is to consider the integral
	\begin{equation}
\int_a^b\omega(x)q(x)p_N(x)dx = 0 \;,
	\end{equation}
where $q(x)$ is an arbitrary polynomial of degree less than $N$. The
quadrature should represent this integral exactly,
	\begin{equation}
\sum_{i=1}^N q(x_i)p_N(x_i) = 0 \;.
	\end{equation}
Apparently this is only possible if $x_i$ are the roots of $p_N$
$\blacksquare$.
\subsubsection{Calculation of nodes and weights}
A neat algorithm---usually refered to as Golub-Welsch~\cite{Golub-Welsch}
algorithm---for calculation of the nodes and weights of a Gauss quadrature
is based on the symmetric form of the three-term reccurence relation
for orthogonal polynomials,
	\begin{equation}\label{q:eq:3r}
x p_{n-1}(x) = \beta_n p_n(x)+\alpha_n p_{n-1}(x)
+\beta_{n-1}p_{n-2}(x)\;,
	\end{equation}
where $p_{-1}(x)\doteq 0$, $p_1(x)\doteq 1$, and $n=1,\dots,N$. This
reccurence relation can be written in the matrix form,
	\begin{equation}\label{q:eq:rm}
x \ve{p}(x)=\ma{J} \ve{p}(x) + \beta_N p_N(x) \ve{e}_N \;,
	\end{equation}
where $\ve{p}(x)\doteq \{p_0(x),\dots,p_{N-1}(x)\}^T$,
$\ve{e}_N=\{0,\dots,0,1\}^T$, and the tridiagonal matrix $\ma{J}$ ---
usually refered to as {\em Jacobi matrix}\index{Jacobi matrix} or {\em
Jacobi operator}\index{Jacobi operator} --- is given as
	\begin{equation}
\ma{J}=\left( \begin{array}{ccccc} \alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 &
\beta_2 \\
	& \beta_2 & \alpha_3 & \beta_3 \\
&&\ddots&\ddots& \\ &&&\beta_{N-1}&\alpha_N \end{array} \right) \,.
	\end{equation}

Substituting the roots $x_i$ of $p_N$---that is, the set $\{x_i\mid
p_N(x_i)=0\}$---into the matrix equation~(\ref{q:eq:rm}) leads to
eigenvalue problem for the Jacobi matrix,
	\begin{equation}
\ma{J} \ve{p}(x_i) = x_i \ve{p}(x_i) \;.
	\end{equation}
Thus, the nodes of an $N$-point Gauss quadrature---the roots of the
polynomial $p_N$---are the eigenvalues of the Jacobi matrix $J$ and
can be calculated by a standard diagonalization\footnote{ A symmetric
tridiagonal matrix can be diagonalized very effectively using the QR/RL
algorithm.} routine $\blacksquare$.

The weights can be obtained considering $N$ integrals,
	\begin{equation}
\int_a^b \omega(x) p_n(x) dx = \delta_{n0}\int_a^b\omega(x)dx \,,
\; n=0,\dots,N-1 \,.
	\end{equation}
Applying our quadrature gives the matrix equation,
	\begin{equation}\label{q:e:pw}
\ma{P} \ve{w} = \ve{e}_1 \int_a^b \omega(x) dx \;,
	\end{equation}
where $\ve{w}\doteq\{w_1,\dots,w_N\}^T$, $\ve{e}_1=\{1,0,\dots,0\}^T$, and
	\begin{equation}
\ma{P}\doteq\left(\begin{array}{ccc} p_0(x_1) & \dots & p_0(x_N) \\ p_1(x_1)
& \dots & p_1(x_N) \\ \dots & \dots & \dots \\ p_{N-1}(x_1) & \dots &
p_{N-1}(x_N) \end{array} \right) \;.
	\end{equation}
Equation~(\ref{q:e:pw}) is linear in $w_i$ and can be solved directly.
However, if diagonalization of the Jacobi matrix provided the normalized
eigenvectors, the weigths can be readily obtained using the following
method.

The matrix $\ma{P}$ apparently consists of non-normalized column
eigenvectors of the matrix $\ma{J}$. The eigenvectors are orthogonal
and therefore $\ma{P}^T\ma{P}$ is a diagonal matrix with positive
elements. Multiplying~(\ref{q:e:pw}) by $\ma{P}^T$ and then by
$(\ma{P}^T\ma{P})^{-1}$ from the left gives
	\begin{equation}
\ve{w} = (\ma{P}^T\ma{P})^{-1}\ma{P}^T \ve{e}_1 \int_a^b\omega(x)dx \;.
	\end{equation}
From $p_0(x)=1$ it follows that $\ma{P}^T\ve{e}_1=\{1,\dots,1\}^T$
and therefore
	\begin{equation}\label{q:e:wi}
w_i = \frac{1}{(\ma{P}^T\ma{P})_{ii}} \int_a^b\omega(x)dx \;.
	\end{equation}
Let the matrix $\ma{V}$ be the set of the normalized column eigenvectors
of the matrix $\ma{J}$. The matrix $\ma{V}$ is then connected with the
matrix $\ma{P}$ through the normalization equation,
	\begin{equation}
\ma{V} = \sqrt{(\ma{P}^T\ma{P})^{-1}}\ma{P} \;.
	\end{equation}
Therefore, again taking into account that $p_0(x)=1$,
equation~(\ref{q:e:wi}) can be written as
	\begin{equation}
w_i = (V_{1i})^2 \int_a^b\omega(x)dx \;\blacksquare.
	\end{equation}

\subsubsection{Example: Gauss-Legendre quadrature}
Gauss-Legendre quadrature deals with the weight $\omega(x)=1$ on the
interval $[-1,1]$.  The associated polynomials are Legendre polynomials
${\mathcal P}_n(x)$, hence the name. Their reccurence relation is
usually given as
	\begin{equation}
(2n-1)x{\mathcal P}_{n-1}(x)=
n{\mathcal P}_n(x)+(n-1){\mathcal P}_{n-2}(x) \,.
	\end{equation}
Rescaling the polynomials (preserving $p_0(x)=1$) as
	\begin{equation}
\sqrt{2n+1}{\mathcal P}_n(x) = p_n(x)
	\end{equation}
reduces this reccurence relation to the symmetric form~(\ref{q:eq:3r}),
	\begin{equation}
xp_{n-1}(x)=
\frac12\frac{1}{\sqrt{1-(2n)^{-2}}}p_n(x)
+ \frac12\frac{1}{\sqrt{1-(2(n-1))^{-2}}}p_{n-2}(x) \;.
	\end{equation}
Correspondingly, the coefficients in the matrix $\ma{J}$ are
	\begin{equation}\label{q:e:glb}
\alpha_n=0\;,\;\;\;
\left\{\beta_n =
\frac12\frac{1}{\sqrt{1-(2n)^{-2}}}\;\Big|\;n=1,\dots,N-1\right\} \;.
	\end{equation}

The problem of finding the nodes and the weights of the $N$-point
Gauss-Legendre quadrature is thus reduced to the eigenvalue problem for
the Jacobi matrix with coefficients~(\ref{q:e:glb}).

	\begin{table}
	\caption{An Octave function which calculates the nodes and weights of the
$N$-point Gauss-Legendre quadrature and then integrates a given function.}
	\label{q:t:o}
	\begin{lstlisting}[language=Octave]
function Q = gauss_legendre(f,a,b,N)
beta = .5./sqrt(1-(2*(1:N-1)).^(-2)); % reccurence relation
J = diag(beta,1) + diag(beta,-1);     % Jacobi matrix
[V,D] = eig(J);                       % diagonalization of J
x = diag(D); [x,i] = sort(x);         % sorted nodes
w = V(1,i).^2*2;                      % weights
Q = w*f((a+b)/2+(b-a)/2*x)*(b-a)/2;   % integral
endfunction;
	\end{lstlisting}
	\end{table}

As an illustration of this algorithm Table~(\ref{q:t:o}) shows an Octave
function which calculates the nodes and the weights of the $N$-point
Gauss-Legendre quadrature and then integrates a given function.

% \section{Chebyshev quadratures (Clenshaw-Curtis)}
% The weights are equal, the abscissas are chosen such that the first
% $n+1$ Chebyshev polynomials $\{T_0(x),T_1(x),\dots,T_{n-1}(x)\}$ are
% integrated exactly. The latter are better behaving and generally provide
% better accuracy. The abscissas are the roots of Chebyshev polynomials,
% which are easily computed.

\subsection{Gauss-Kronrod quadratures}
Generally, the error of a numerical integration is estimated by comparing
the results from two rules of different orders.  However, for ordinary
Gauss quadratures the nodes for two rules of different orders almost
never coinside.  This means that one can not reuse the points of the
lower order rule when calculating the hihger order rule.

{\em Gauss-Kronrod} algorithm~\cite{kronrod} remedies this inefficiency.
The points inherited from the lower order rule are reused in the higher
order rule as predefined nodes (with $n$ weights as free parameters),
and then $m$ more optimal points are added ($m$ abscissas and $m$
weights as free parameters). The order of the method is $n+2m-1$.
The lower order rule becomes {\em embedded}\index{embedded rule}---that
is, it uses a subset of the nodes---into the higher order rule.  On the next
iteration the procedure is repeated.

Patterson~\cite{patterson} has tabulated nodes and weigths for several
sequences of embedded Gauss-Kronrod rules.


\section{Adaptive quadratures}

Higher order quadratures suffer from round-off errors as the weights $w_i$
generally have alternating signs. Again, using high order polynomials is
dangerous as they typically oscillate wildly and may lead to Runge's
phenomenon. Therefore, if the error of the quadrature is yet too
large for a quadrature with sufficiently large~$n$, the best strategy
is to subdivide the interval in two and then use the quadrature on
the half-intervals.  Indeed, if the error is of the order $h^k$, the
subdivision would lead to reduced error, $2\left(h/2\right)^k< h^k$,
if $k>1$.

An {\em adaptive quadrature} is an algorithm where the integration
interval is subdivided into adaptively refined subintervals until the
given accuracy goal is reached.

Adaptive algorithms are usually built on pairs of quadrature rules
-- a higher order rule,
	\begin{equation}
Q = \sum_i w_i f(x_i),
	\end{equation}
where $w_i$ are the weights of the higher order rule and $Q$ is the
higher order estimate of the integral, and a lower order rule,
	\begin{equation}
q = \sum_i v_i f(x_i),
	\end{equation}
where $v_i$ are the weights of the lower order rule and $q$ is the the
lower order estimate of the integral. The difference between the higher
order rule and the lower order rule gives an estimate of the error,
	\begin{equation}
\delta Q = |Q-q| \:.
	\end{equation}
The integration result is accepted, if the error $\delta Q$ is smaller
than tolerance,
	\begin{equation}
\delta Q < \delta + \epsilon|Q|\:,
	\end{equation}
where $\delta$ is the absolute accuracy goal and $\epsilon$ is the
relative accuracy goal of the integration.

If the error estimate is larger than tolerance, the interval is
subdivided into two half-intervals and the procedure applies recursively
to subintervals with the same relative accuracy goal $\epsilon$ and
rescaled absolute accuracy goal $\delta/\sqrt{2}$.

The points $x_i$ are usually chosen such that the two quadratures use
the same points, and that the points can be reused in the subsequent
recursive steps.  The reuse of the function evaluations made at the
previous step of adaptive integration is very important for the efficiency
of the algorithm. The equally-spaced abscissas naturally provide for
such a reuse.

As an example, Table~\ref{q:t:aq} shows an implementation of
the described algorithm using
	\begin{eqnarray}
x_i &=& \left\{\frac16,\frac26,\frac46,\frac56\right\}
\, \text{(easily reusable points)}\,, \\ w_i &=&
\left\{\frac26,\frac16,\frac16,\frac26\right\} \, \text{(trapezium
rule)}\,, \\ v_i &=& \left\{\frac14,\frac14,\frac14,\frac14\right\} \,
\text{(rectangle rule)} \,.
	\end{eqnarray}
During recursion the function values at the points \#2 and \#3 are inherited
from the previous step and need not be recalculated.

The points and weights are cited for the rescaled integration interval
$[0,1]$.  The transformation of the points and weights to the original
interval $[a,b]$ is given as
	\begin{eqnarray}
x_i &\rightarrow& a+(b-a)x_i \:, \nonumber \\
w_i &\rightarrow& (b-a)w_i \:.
	\end{eqnarray}

This implementation calculates directly the Riemann sums and can therefore
deal with integrable singularities, although rather inefficiently.

More efficient adaptive routines keep track of the subdivisions of the
interval and the local errors~\cite{gonnet}.  This allows detection
of singularities and switching in their vicinity to specifically tuned
quadratures.  It also allows better estimates of local and global errors.

\begin{table}[thb]
\caption{Recursive adaptive integrator in C}\label{q:t:aq}
\lstinputlisting{c/integration/adapt.c}
\end{table}

\section{Variable transformation quadratures}
The idea behind {\em variable transformation quadratures} is to apply the
given quad\-rature---either with optimimized or regularly spaced
nodes---not to the original integral, but to a variable transformed
integral~\cite{MoriDE},
	\begin{equation}
\int_a^b f(x)dx = \int_{t_a}^{t_b} f\big(g(t)\big)g'(t)dt
\approx \sum_{i=1}^N w_i f\big(g(t_i)\big)g'(t_i) \,,
	\end{equation}
where the transformation $x=g(t)$ is chosen such that the transformed
integral better suits the given quadrature. Here $g'$ denotes the
derivative and $[t_a, t_b]$ is the corresponding interval in the new
variable.

For example, the Gauss-Legendre quadrature assumes the integrand can be
well represented with polynomials and performs poorly on integrals with
integrable singularities like
	\begin{equation}
I=\int_0^1 \frac{1}{2\sqrt x}dx \,.
	\end{equation}
However, a simple varibale transformation $x=t^2$ removes the singularity,
	\begin{equation}
I=\int_0^1 dt \,,
	\end{equation}
and the Gauss-Legendre quadrature for the transformed integral gives
exact result.  The price is that the transformed quadrature performs
less effectively on smooth functions.

Some of the popular variable transformation quadratures are
Clenshaw-Curtis~\cite{Clenshaw-Curtis}, based on the transformation
	\begin{equation}
\int_{-1}^1 f(x)dx = \int_0^\pi f(\cos\theta)\sin\theta d\theta
%\approx\sum_{i=1}^{N} \frac{\pi}{N}
%f(\cos\frac{2i-1}{2N}\pi)\sin\frac{2i-1}{2N}\pi
\,,
	\end{equation}
and ``tanh-sinh'' quadrature~\cite{MoriDE}, based on the transformation
	\begin{equation}
\int_{-1}^1 f(x)dx = \int_{-\infty}^{\infty}
f\left(\tanh\left(\frac{\pi}{2}\sinh(t)\right)\right)
\frac{\pi}{2}\frac{\cosh(t)}{\cosh^2\left(\frac{\pi}{2}\sinh(t)\right)}
dt\,.
	\end{equation}

Generally, the equally spaced trapezium rule is used after the
transformation.


\section{Infinite intervals}

One way to calculate an integral over infinite interval is to transform
it by a variable sustitution into an integral over a finite interval. The
latter can then be evaluated by ordinary integration methods.
Table~\ref{q:t:inf} lists several of such transformation.

\begin{table}[htb]
\caption{Variable transformations reducing infinite interval integrals
into integrals over finite intervals.}
\label{q:t:inf}
	\begin{eqnarray}
\int_{-\infty}^{+\infty}f(x)dx &=&
\int_{-1}^{+1}f\left(\frac{t}{1-t^2}\right)\frac{1+t^2}{(1-t^2)^2}dt
\;,\\
\int_{-\infty}^{+\infty} f(x) dx &=& \int_0^1
\left(f\left(\frac{1-t}{t}\right) + f\left(-\frac{1-t}{t}\right)\right)
{dt\over t^2}
\;,\\
\int_{a}^{+\infty}f(x)dx &=&
\int_{0}^{1}f\left(a+\frac{t}{1-t}\right)\frac{1}{(1-t)^2}dt
\;,\\
\int_{a}^{+\infty} f(x) dx &=&
\int_0^1 f\left(a + \frac{1-t}{t}\right) \frac{dt}{t^2}
\;,\\
\int_{-\infty}^{a}f(x)dx &=&
\int_{-1}^{0}f\left(a-\frac{t}{1+t}\right)\frac{-1}{(1+t)^2}dt
\;,\\
\int_{-\infty}^{a} f(x) dx &=&
\int_0^1 f\left(a - \frac{1-t}{t}\right)
\frac{dt}{t^2}
\;.
	\end{eqnarray}
\hrule
	\end{table}
