5.2 Adaptive Integration 115
The polynomial interpolatory quadrature of ordernwith equidistant quadrature points
xk=a+khand steph= (b−a)/nis called the Newton-Cotes quadrature formula of ordern.
General expressions can be found in for example Refs. [24,25].
5.2 Adaptive Integration.
Before we proceed with more advanced methods like Gaussian quadrature, we mention
breefly how an adaptive integration method can be implemented.
The above methods are all based on a defined step length, normally provided by the user,
dividing the integration domain with a fixed number of subintervals. This is rather simple
to implement may be inefficient, in particular if the integrand varies considerably in certain
areas of the integration domain. In these areas the number offixed integration points may
not be adequate. In other regions, the integrand may vary slowly and fewer integration points
may be needed.
In order to account for such features, it may be convenient tofirst study the properties
of integrand, via for example a plot of the function to integrate. If this function oscillates
largely in some specific domain we may then opt for adding moreintegration points to that
particular domain. However, this procedure needs to be repeated for every new integrand
and lacks obviously the advantages of a more generic code.
The algorithm we present here is based on a recursive procedure and allows us to automate
an adaptive domain. The procedure is very simple to implement.
Assume that we want to compute an integral using say the trapezoidal rule. We limit our-
selves to a one-dimensional integral. Our integration domain is defined byx∈[a,b]. The algo-
rithm goes as follows
- We compute our first approximation by computing the integral for the full domain. We label
this asI(^0 ). It is obtained by calling our previously discussed functiontrapezoidal_ruleas
I0 = trapezoidal_rule(a, b, n, function); - In the next step we split the integration in two, withc= (a+b)/ 2. We compute then the two
integralsI(^1 L)andI(^1 R)
I1L = trapezoidal_rule(a, c, n, function);
and
I1R = trapezoidal_rule(c, b, n, function);With a given defined tolerance, being a small number providedby us, we estimate the
difference|I(^1 L)+I(^1 R)−I(^0 )|<tolerance. If this test is satisfied, our first approximation is
satisfactory.- If not, we can set up a recursive procedure where the integral is split into subsequent
subintervals until our tolerance is satisfied.
This recursive procedure can be easily implemented via the following function
// Simple recursive function that implements the
// adaptive integration using the trapezoidal rule
// It is convenient to define as global variables
// the tolerance and the number of recursive steps
const intmaxrecursions = 50;
const doubletolerance = 1.0E-10;
// Takes as input the integration limits, number of points, function to integrate
// and the number of steps