Chapter 5 - Quadrature

The Goal: Adaptive Quadrature

Historically in mathematics, quadrature refers to the

act of trying to find a square with the same area as a
given circle.
In mathematical computing, quadrature refers to the
numerical approximation of definite integrals.
Let f (x) be a real-valued function of a real variable,
defined on a finite interval a ≤ x ≤ b.
We seek to compute the value of the definite integral
Z b
f (x) dx.
a

Note: This is a real number.

In line with its historical meaning, “quadrature” might
conjure up the vision of plotting the function on
graph paper and counting the little boxes (and their
fractions!) that lie underneath the curve.

1
Observations:

• The use of smaller boxes gives a more accurate

answer (but requires more work!).

• It is possible to get away with larger boxes in some

regions (where f (x) doesn’t change much) and not
make a large error.

Adaptive quadrature is the act of carefully selecting

the points where f (x) is sampled.

Of course, we want a computer to do this automatically

for us!

So, we will have to figure out strategies that can be

programmed.

In general, we also want to evaluate the function

at as few points as possible while approximating
the integral to within some specified accuracy; i.e.,
adaptive quadrature is meant to give accurate answers
in an efficient manner.

2
The basis for adaptive quadrature is the following
fundamental additive property of the definite integral.

If c is any point between1 a and b, then

Z b Z c Z b
f (x) dx = f (x) dx + f (x) dx.
a a c

The idea of adaptive quadrature is that if we can

approximate each of the two integrals on the right to
within a specified tolerance, then the sum gives us the
desired result.

If not, we can recursively apply the additive property

to each of the intervals [a, c] and [c, b]; e.g., if the error
is too large on [a, c] the we subdivide it into smaller
subintervals in order to reduce the overall error.

1
Actually, the ensuing equation usually is true even if c 6∈ [a, b] (unless
f (x) is misbehaved), but such values of c are not used by adaptive quadrature.

3
We expect the quadrature strategy to adapt to the
integrand automatically, partitioning the interval into
subintervals with fine spacing where the integrand is
varying rapidly and coarse spacing where the integrand
is varying slowly.

However, before we can learn the details of

adaptive quadrature, we need to understand the basic
quadrature rules.

4
 Basic Quadrature Rules

We begin by deriving of two basic yet useful quadrature

rules: the midpoint rule and the trapezoidal rule.

Both of these rules take their names from the

approximations used to perform the quadrature.

Let h = b − a be the length of the integration interval.

The midpoint rule M approximates the integral by the

area of a rectangle whose base has length h and whose
height is the value of f (x) at the midpoint:
 
a+b
M = hf .
2

5
The trapezoidal rule T approximates the integral by
the area of a trapezoid with base h and sides equal to
the values of f (x) at the two end points.
 
f (a) + f (b)
T =h .
2

What makes one quadrature rule more accurate than

another?

We would like to be able to quantify this somehow.

This leads us to the following definition:

The order of a quadrature rule is the degree of

the lowest-degree polynomial that the rule does not
integrate exactly.

6
An equivalent interpretation of order (and accuracy) of
a quadrature rule is as follows:

If a quadrature rule of order p is used to integrate

a smooth function over a small interval of length h,
then a Taylor series analysis shows that the error is
proportional to hp.

This means that if h is decreased by a factor of 2, the

error decreases by a factor of 2p.

So, rules with high order converge more quickly than

those with low order.

The midpoint rule and the trapezoidal rule are both

exact for constant functions and linear functions of x.

This should be obvious for the trapezoidal rule:

The trapezoidal rule forms a linear interpolant between

(a, f (a)) and (b, f (b)) and integrates the interpolant
exactly to define the rule.

7
If the underlying function is linear (so that the
interpolant and the underlying function are identical),
it should not be surprising that the trapezoidal rule
gives the exact result.

The order is a little more subtle for the midpoint rule!

Already there is the notion of a free lunch here: the

midpoint rule does a constant (degree-0) interpolation,
yet it can integrate linear (degree-1) functions exactly!

Coming back to reality, neither of them is exact for a

quadratic in x, so they both have order 2.

There are two quadrature rules that you may have

seen but that we did not cover: the so-called left-hand
rectangle rule and the right-hand rectangle rule.

These rules approximate the area under f (x) by

rectangles with height f (a) or f (b), respectively.

It can be shown that their orders are only 1, and so

they are not good methods to use in practice.

8
The accuracy of the midpoint and trapezoidal rules
can be compared by examining their behavior on the
simple definite integral
Z 1
2 1
x dx = .
0 3

The midpoint rule gives

 2
1 1
M =1 = .
2 4

The trapezoidal rule gives

02 + 1 2 1
T =1 = .
2 2

So, we see that the error in M is 1/12, whereas the

error in T is −1/6.

9
The errors have opposite signs and, perhaps
surprisingly, the midpoint rule is twice as accurate
as the trapezoidal rule.

This turns out to not be just coincidence!

When integrating smooth f (x) over short intervals

[a, b], M is roughly twice as accurate as T , and the
errors have opposite signs.

We can extrapolate on this knowledge as follows:

Knowing these error estimates allows us to combine

these two methods and get a rule that is usually more
accurate than either one separately!

Let the exact value of the integral be I.

If the error in T were exactly −2 times the error in M ,

then
I − T = −2(I − M ).

Solving this equation for I would then give us the exact

value of the integral!

10
We settle for creating a new rule S that better
approximates I:

2 1
S = M + T.
3 3
Indeed, S is usually a more accurate approximation
than either M or T alone.

This rule is known as Simpson’s rule.

Classically, Simpson’s rule is derived by using a

parabola to interpolate f (x) at the two endpoints, a
and b, and the midpoint, c = (a + b)/2 and integrating
the parabola exactly.

This yields

h
S = (f (a) + 4f (c) + f (b)).
6

It turns out that S also integrates cubics exactly

(another free lunch!), but not quartics, so its order
is 4.

11
So far, we have fit one interpolant to the entire interval
[a, b] and integrated it to obtain a quadrature rule.

One way to derive quadrature rules of higher and

higher order is to fit polynomial interpolants of higher
and higher degree to the integrand f (x) and integrate
these interpolants exactly.

However, as we know, high-degree polynomial

interpolation based on equally spaced abscissae
generically results in oscillatory interpolants that do
not reflect the behaviour of the underlying function.

This would degrade the accuracy of the quadrature

rules based on them!

With equally spaced abscissae, a better strategy is to

use lower-order piecewise polynomial interpolation.

For example, apply Simpson’s rule on the two halves

of the interval, [a, c] and [c, b].

Let d and e be the midpoints of these two subintervals;

i.e., d = (a + c)/2 and e = (c + b)/2.

12
Applying Simpson’s rule to each subinterval, we obtain
a quadrature rule over [a, b]:

h
S2 = (f (a) + 4f (d) + 2f (c) + 4f (e) + f (b)).
12

This is an example of a composite quadrature rule.

As before, the two approximations S and S2 can be

combined to get an even more accurate approximation.

Both rules are of order four, but the step size for S2 is
half the step size for S, so the error in S2 is roughly
24 = 16 times smaller than that of S.

Thus, a new, more accurate rule Q can be obtained by

solving
Q − S = 16(Q − S2).
The result is

(S2 − S)
Q = S2 + .
15
13
As an exercise on the mock final asks you to express Q
as a weighted combination of the five function values
f (a) through f (e) and to establish that its order is six.

The rule is known as Weddle’s rule or the sixth-order

closed Newton–Cotes rule or the first step of Romberg
integration.

We will simply call it the extrapolated Simpson’s rule

because it uses Simpson’s rule for two different values
of h and then extrapolates toward h = 0.

Romberg integration generalizes this to the use of an

arbitrary number of different values of h.

14
 An Adaptive Quadrature Strategy

Composite quadrature rules with error estimates can

be used to produce adaptive quadrature methods.
In general, it is inefficient to use a uniform mesh on an
integrand f (x) that varies rapidly in some regions but
slowly in others.
A uniform mesh has a larger error where the integrand
varies rapidly; so the mesh has to be fine there.
Now if the mesh is uniform, we would have to use a
very fine mesh everywhere.
This can be extremely inefficient in regions where f (x)
varies slowly and a coarser mesh would be adequate.
Ideally, we would like to equidistribute the error over
the whole mesh; i.e., we would like each piece of the
composite integration to contribute the same amount
of error to the overall error.
Then we would be using a fine or coarse mesh where
appropriate.

15
 Rb
The basic problem is to approximate a
f (x) dx to
within some specified tolerance  > 0.

Suppose we wish to base our adaptive method on

Simpson’s rule.

Step 1: Apply Simpson’s rule with h = b − a.

Z b
f (x) dx ≈ Sh(a, b),
a

where
   
h a+b
Sh(a, b) = f (a) + 4f + f (b) .
6 2

Step 2: Determine an estimate of the accuracy of the

approximation.

One way to do this is to halve the mesh and re-solve

the problem; i.e., compute Sh/2(a, b).

16
Recalling Weddle’s rule (translated to the present
notation),

Q(a, b) = Sh/2(a, b) + Eh/2,

where
Sh/2(a, b) − Sh(a, b)
Eh/2 = .
15

This can be viewed as an approximation Sh/2(a, b) to

the definite integral, plus an error estimate Eh/2.

In the case of Weddle’s rule, we add the error estimate

to the approximation to (hopefully) obtain an even
better approximation.

In the case of adaptive quadrature, we use Eh/2

to decide whether the approximation Sh/2(a, b) is
accurate enough.
Rb
If so, Sh/2(a, b) should approximate a f (x) dx to
within .

17
If not, we repeat this process on the subintervals
[a, (a + b)/2] and [(a + b)/2, b] individually.

However, this time we will compare the error estimates

to /2.

This way, if each of the 2 subintervals has an error of at

most /2, the total error will be at most /2 + /2 = .

Again, if any subinterval reports an error of more than

/2, it is halved again, and its error is recomputed and
compared to /4, etc.

We keep going until the error estimate on every

subinterval does not exceed the appropriate fraction
of .

18
Two important notes:

1. The text by Moler argues that this strategy is too

conservative, leading to too much accuracy.
In fact, some textbooks advocate simply using the
same accuracy tolerance  on every subinterval.

2. This method is not foolproof!

Integrands with discontinuities, singularities, or
misbehaved higher derivatives can cause this
method to never terminate.
However, it tends to work well on many problems,
roughly because for each subdivision, the accuracy
increases by a factor of 16, but the tolerance is only
restricted by a factor of 2.

Example: Evaluate
Z 3
100 10
2
sin dx.
1 x x

See quadDemo.m

19
 Specifying Integrands

Matlab’s built-in global adaptive quadrature routine

is called integral.

The basic syntax is q = integral(f,a,b), which

integrates f over the interval a to b.

There are a few different ways of specifying f.

For simple formulas, we can use anonymous functions.

For example, to compute

Z 1
dx
√ ,
1+x 4
0

we can use

f = @(x) 1./sqrt(1+x.^4)
q = integral(f,0,1)

20
For more complicated functions, we can create a file
(for the above example, say called integrand.m)
function f = integrand(x)
f = 1./sqrt(1+x.^4);
end
Then the statement
q = integral(@integrand,0,1)
uses a function handle to compute the integral.
The integral routine can sometimes handle
integrable singularities at endpoints automatically.
For example, if we want to compute
Z π
sin x
dx,
0 x

we can type
q = integral(@(x) sin(x)./x, 0, 1)

21
However, integral is not likely to work for integrands
that have singularities stronger than about x−1/2 as
x → 0+, even if mathematically the integral is well
defined.

For example, if we try

q = integral(@(x) x.^(-0.9), 0, 1)
the output is
Warning: Infinite or Not-a-Number value encountered.
> In funfun/private/integralCalc>iterateScalarValued at 349
In funfun/private/integralCalc>vadapt at 132
In funfun/private/integralCalc at 75
In integral at 88

ans =

Inf

even though the exact answer is of course 10.

22
Definite integrals that depend on parameters are also
encountered frequently.

An example is the beta function, defined by

Z 1
β(z, w) = tz−1(1 − t)w−1 dt.
0

Matlab already has a beta function, but we can use

this example to illustrate how to handle parameters.
Create an anonymous function with three arguments,
F = @(t,z,w) t.^(z-1).*(1-t).^(w-1);

or create an M-file with a name like my beta.m,

function f = my beta(t,z,w)
f = t.^(z-1).*(1-t).^(w-1)

As with all functions, the order of the arguments is

important!

23
The functions used with quadrature routines must have
the variable of integration as the first argument.

Values for the parameters are then given as extra

arguments.

For example, to compute β(8/3, 10/3), you could do

z = 8/3;
w = 10/3;
atol = 1e-10;
rtol = 1e-6;
and then use
Q = integral(@t f(t,z,w),0,1,’AbsTol’, atol, ’RelTol’, rtol);

or
Q = integral(@(t) my beta,0,1,’AbsTol’, atol, ’RelTol’, rtol);

The “function functions” in Matlab usually expect

the first argument to be in vectorized form.

24
In other words, the function you input to the function
function should be able to accept a vector argument;
i.e., it should be written in such a way as to be able to
process a bunch of inputs simultaneously .

Usually, this means that the mathematical expressions

should be specified with Matlab array (.) notation.

For example,
sin x
f (x) =
1 + x2
should be input as

sin(x)./(1 + x.^2)

Without the (.) operators, the code

sin(x)/(1 + x^2)

calls for linear algebraic vector operations that are not

appropriate here.

The Matlab function vectorize automatically

transforms a scalar expression into a form that can
be used as input to function functions.

25
So

vectorize(’sin(x)/(1+x^2)’)

produces

sin(x)./(1+x.^2)

Many of the function functions in Matlab require the

specification of an interval of the x-axis.

Mathematically, we have two notations, a ≤ x ≤ b or

x ∈ [a, b].

With Matlab, we also have two possibilities.

The endpoints can be given as two separate arguments,

a and b, or they can be combined into one vector
argument, [a,b].

The quadrature functions use two separate arguments.

As we have seen, the root-finder fzero uses a single

argument, but this can be either a single starting point
or a two-element vector specifying the interval.

26
 Integrating Discrete Data

So far, we have been concerned with computing an

approximation to the definite integral of a specified
function f (x).

We have assumed that we can evaluate f (x) at any

point in a given interval.

However, in many situations, the function is only

known at a discrete set of points, say (xk , yk ),
k = 1, 2, 3, . . . , n. Assume the xk are sorted in
increasing order, with

a = x1 < x2 < · · · < xn−1 < xn = b.

Rb
When approximating a f (x) dx, it is now not possible
to sample f (x) wherever we want.

This rules out the use of adaptive quadrature methods!

27
The most obvious approach is to do piecewise linear
interpolation through the data and integrate the
interpolant exactly.

i.e., use the composite trapezoidal rule

n−1
X yk + yk+1
T = hk ,
2
k=1

where hk = xk+1 − xk .

The composite trapezoidal rule can be implemented in

Matlab in one line:

T = sum(diff(x).*(y(1:end-1)+y(2:end))/2)

The Matlab function trapz also provides an

implementation.

28
For example, consider the data set

x = 1:6
y = [6 8 11 7 5 2]

For these data, the trapezoidal rule gives

T = 35.

The trapezoidal rule is often satisfactory in practice,

and more complicated methods may not be necessary.

Nevertheless, methods based on higher-order

interpolation can give other estimates of the integral.

Whether or not they are better is impossible to decide

without further knowledge of the data!

29
Recall that both the spline and pchip interpolants
are based on piecewise cubic Hermite interpolation:

3hs2 − 2s3 h3 − 3hs2 + 2s3

P (x) = yk+1 + yk
h3 h3
s2(s − h) s(s − h)2
+ dk+1 + dk ,
h2 h2

where xk ≤ x ≤ xk+1, s = x − xk , and h = hk .

Recall also that this is a cubic polynomial in s (and

hence in x) that satisfies 4 interpolation conditions, 2
on function values and 2 on derivative values:

P (xk ) = yk , P (xk+1) = yk+1,

P 0(xk ) = dk , P 0(xk+1) = dk+1.

The slopes dk are computed in spline or pchip.

30
As an exercise, show that
xk+1
2 dk+1 − dk
Z
yk+1 + yk
P (x) dx = hk − hk .
xk 2 12

Hence Z b
P (x) dx = T − D,
a
where T is the trapezoidal rule, and

n−1
X dk+1 − dk
D= h2k ,
12
k=1

is a higher-order correction.

Note: If the xk are equally spaced with spacing h, then

D simplifies to

dn − d1
D = h2 .
12
31
For the data shown earlier, the area obtained by trapz
(piecewise linear interpolation) is 35.00.

By spline interpolation, the area is 35.25.

By the pchip, the area is 35.41667.

The integration process tends to average out the

variation in the interpolants, so even though the
interpolants might have rather different shapes, the
approximations to the integral are often quite close to
each other.

32
 High-Powered Tricks and Theory

• The method of undetermined coefficients

Quadrature formulas take the form

Z b n
X
f (x) dx ≈ wif (xi);
a i=1

i.e., the definite integral is approximated by a weighted

sum of function values.

We have seen how quadrature rules can be derived

based on interpolating a set of function values by a
polynomial or piecewise polynomial of a given degree
and integrating the interpolant exactly.

An alternative derivation to explicitly forming

and integrating an interpolant is the method of
undetermined coefficients.

33
Suppose the nodes xi have been chosen; e.g., they
could be n equally spaced points in [a, b], including the
endpoints.

The idea is that we treat the weights wi as unknowns

and solve for them by forcing the quadrature rule to
integrate the usual polynomial basis functions exactly.

If we have n nodes, we have n weights; so, we

can require the quadrature rule to exactly integrate
polynomials up to degree n − 1.

→ a system of n linear equations in n unknowns.

Example: Derive a 3-point quadrature rule that uses

the nodes x1 = a, x2 = (a + b)/2, and x3 = b.

The quadrature rule is

Z b  
a+b
f (x) dx ≈ w1f (a) + w2f + w3f (b).
a 2

We want the rule to be exact for f (x) = 1, x, x2.

34
(Then by the linearity property of integration, the rule
will be exact for all polynomials of up to degree 2.)

f (x) = 1 :
Z b
w1 · 1 + w2 · 1 + w3 · 1 = 1 dx = b − a
a

f (x) = x :

b
b2 − a2
Z
w1 · a + w2 · w2 + w3 · b = x dx =
a 2

f (x) = x2 :
 2 b 3 3
−
Z
2 a+b b a
w1 · a + w2 · + w3 · b2 = x2 dx =
2 a 3

35
Solution:

b−a 2(b − a) b−a

w1 = , w2 = , w3 = .
6 3 6

This is of course Simpson’s rule.

See undetCoeffs.m

Try Simpson’s rule on f (x) = x3.

Are you surprised by the result?

36
• Gaussian quadrature

In the quadrature rules that we have studied so far, we

have always assumed that the integrand is evaluated at
the end points and at equally spaced points in between.

There are disadvantages to this!

1. Insisting that the nodes be equally spaced

removes degrees of freedom from the method of
undetermined coefficients.
→ If we remove this restriction, we can derive
higher-order quadrature rules for a given number of
function evaluations!

2. Integrands are often singular (and yet integrable) at

the endpoints.
→ Evaluating f (x) at x = a or x = b is problematic,
even though the integral itself is well defined.

37
Recall that when we fixed the xi, we chose the wi to
integrate exactly the functions f (x) = 1, x, x2, . . ..

So for n nodes, we expect to integrate exactly

polynomials of degree n − 1.

The key idea behind Gaussian quadrature is to also

treat the xi as unknowns.

Now we have 2n degrees of freedom, and we expect to

integrate exactly polynomials up to degree 2n − 1.

Example: Derive the Gaussian quadrature rule with

n = 2 and [a, b] = [−1, 1].

We are to determine x1, x2, w1, w2 such that

Z 1
f (x) dx = w1f (x1) + w2f (x2)
−1

holds exactly for all polynomials of degree 3 or less.

38
 Z 1
f (x) = 1 : w1 · 1 + w2 · 1 = 1 dx = 2
−1
Z 1
f (x) = x : w1 · x1 + w2 · x2 = x dx = 0
−1
Z 1
2 2
f (x) = x : w1 · x21 + w2 · x22 = x dx = 2
−1 3
Z 1
f (x) = x3 : w1 · x31 + w2 · x32 = x3 dx = 0
−1

One solution is

1 1
x1 = − √ , x2 = √ , w1 = 1, w2 = 1.
3 3

Thus, the Gaussian quadrature rule takes the form

Z 1    
1 1
f (x) dx ≈ f − √ +f √ .
−1 3 3

39
It turns out that the abscissae for these maximum-
order (2n) quadrature rules are the roots of the nth
Legendre polynomial.
Hence, these methods are also sometimes called Gauss–
Legendre quadrature rules.
Knowledge of the location of the xi makes the
computation of the wi much easier!
→ A linear system of size n can be solved for the wi,
as opposed solving to a nonlinear system of size 2n for
the xi and wi simultaneously.
Note: The Legendre polynomials are defined on the
interval [−1, 1].
So in general, the definite integral defined on [a, b] is
transformed to [−1, 1] to make the evaluations more
convenient.
This means we change the variable of integration by
means of

2x − b − a 1
t= ⇐⇒ x = [(b − a)t + (b + a)].
b−a 2
40
Thus,

b 1  
b−a
Z Z
1
f (x) dx = f [(b − a)t + (b + a)] dt.
a −1 2 2

−x2
R1
Example: Approximate 0
e dx by the 2-point
Gaussian quadrature rule.

First, we transform the interval of integration:

1
a = 0, b = 1 =⇒ x = (t + 1).
2

Thus
1
1 − 0 1 −{ 1 (t+1)}2
Z Z
2
e−x dx = e 2 dt
0 2 −1
 n  o2 n  o2 
1 − 12 − √13 +1 − 1 √1 +1
≈ e +e 2 3
2
= 0.746595.

41
Note: The exact answer is 0.746824.

This approximation is more accurate than that obtained

using (non-composite) Simpson’s rule and with one
fewer function evaluation!

See gaussQuad.m

42
• Infinite intervals of integration

The problem we consider here is

Z b
f (x) dx,
a

where now a = −∞ and/or b = ∞.

There are 3 standard ways to approach such problems.

1. Replace the infinite limit(s) by finite value(s).

In order to do this sensibly, an error estimate of the
neglected “tail(s)” is required.

2. Change variables so that the new interval is finite;

e.g.,
 
1
x = − log t = log ⇐⇒ t = e−x,
t
x = 0 =⇒ t = 1, x = ∞ =⇒ t = 0,

43
 or

t x
x= ⇐⇒ t = ,
1−t x+1
x = 0 =⇒ t = 0, x = ∞ =⇒ t = 1.

3. Apply a quadrature rule like Gauss–Laguerre (for

[0, ∞)) or Gauss–Hermite (for (−∞, ∞)).

We briefly describe Gauss–Laguerre quadrature.

Gauss–Laguerre quadrature rules are designed to
approximate definite integrals
Z ∞ n
X
e−xf (x) dx ≈ wif (xi)
0 i=1

by choosing the quadrature weights wi and nodes xi

such that the quadrature rule is exact for f (x) that are
polynomials up to and including degree 2n − 1.
We can derive the 2-stage Gauss–Laguerre quadrature
rule using the method of undetermined coefficients.

44
 Z ∞
f (x) = 1 : w1 · 1 + w2 · 1 = e−x(1) dx = 1
0
Z ∞
f (x) = x : w1 · x1 + w2 · x2 = e−xx dx = 1
0
Z ∞
f (x) = x2 : w1 · x21 + w2 · x22 = e−xx2 dx = 2
0
Z ∞
f (x) = x3 : w1 · x31 + w2 · x32 = e−xx3 dx = 6
0

Solution:
√ √
√ √ 2+ 2 2− 2
x1 = 2− 2, x2 = 2+ 2, w1 = , w2 = .
4 4

45
Implementation note: When solving
Z ∞
f˜(x) dx,
0

we re-write the integral as

Z ∞
e−x[ex f˜(x)] dx,
0

and use the quadrature rule in the form

Z ∞ n
X
f˜(x) dx ≈ wiexi f˜(xi);
0 i=1

i.e., we identify f (x) with exf˜(x).

Example: Use 2-stage Gauss–Laguerre quadrature to
approximate Z ∞
ln(1 + e−x) dx.
0

See infiniteInterval.m

46
• Multiple integrals (cubature)

Multiple integrals are integrals in dimension greater

than 1.

Needless to say, the computational effort required to

compute definite integrals generally increases greatly
with increasing dimension.

Fortunately, definite integrals derived from physical

problems are usually limited to 3 or 4 dimensions.

However, applications such as mathematical finance

can have integrals in an arbitrary number of
dimensions; e.g., depending on how many different
financial products are relevant.

There are 3 general approaches for computing multiple

integrals.

1. Write the integral as an iterated integral and use

1D quadrature as with Matlab’s integral2 or
integral3 functions.
This approach works best for rectangular domains.

47
The syntax for integral2 is

q = integral2(fun,xmin,xmax,ymin,ymax)

xmin and xmax must be constants.

If ymin and ymax are constants, then the domain of

integration is a rectangle.

If the domain has a nice mathematical description,

ymin and ymax can be functions of x.

e.g., ymin=0, ymax=@(x) x, xmin=0, xmax=1

describes the triangle (0, 0), (1, 0), and (1, 1).

→ Keeping reasonable tolerances at each stage can

be tricky!

For complicated domains, one can define a rectangle

that contains the domain and use the characteristic
function χ(x, y) that takes the value 1 if (x, y) is in
the domain and 0 otherwise.

For example, the domain x2 + y 2 ≤ 1 could

be described using xmin,ymin=-1, xmax,ymax=1,
and χ = max(sgn(1 − x2 − y 2), 0).

48
2. Use a higher-dimensional quadrature rule.
This generalizes what we have been doing to higher
dimensions; e.g., a 2D quadrature rule looks like
Z bZ d n
X
f (x, y) dy dx ≈ wif (xi, yi).
a c i=1

Again, this approach works well for regular domains;

various rules are derived for other domains.

3. Use a Monte Carlo method.

This is often the method of choice to approximate
definite integrals in dimensions higher than 2 or 3.
The idea is to sample the integrand at N points
distributed randomly in the domain of integration.
These values are averaged and multiplied by the
“volume” of the domain to approximate the integral.
It can be shown that such a procedure converges
like N −1/2, which is pretty slow!

49
 For example, to gain an extra decimal of accuracy,
you will need 100 times more points!
This makes Monte Carlo uncompetitive for low-
dimensional quadrature.
However, the beauty is that the convergence rate is
independent of the number of dimensions!
Also, the software is easy to write and parallelize.
Monte Carlo methods can be made more efficient
by biasing the sampling, e.g., concentrating where
the integrand is large or rapidly varying.
→ This makes it behave like adaptive quadrature.

Example: Estimate
Z 1Z 1
1
dy dx.
0 0 1 − xy

Note: Exact value is π 2/6.

See multipleIntegral.m

50
 Shared-Memory Adaptive Quadrature

The natural progression of scientific discovery and

industrial research is to continually seek out harder
and harder problems to study.

Parallel programming and parallel computers allow us

to take better advantage of computing resources to
solve larger and more memory-intensive problems in
less time than is possible on a single serial computer.

Trends in computer architecture towards multi-

core systems make parallel programming and high-
performance computing (HPC) necessary practices
going forward.

Computational scientists who ignore HPC do so at

their own risk!

It is already difficult to do leading-edge computational

science without HPC, and this situation can only get
worse in the future.

51
 The Shared-Memory Paradigm

Parallel computers can be built in different ways.

In a shared-memory environment, each processor is able

to read from and write to the same (global) memory.

Each processor (or node) also has its own local memory
(cache) and runs its own program.

CPU CPU CPU

Cache Cache Cache

High−speed Interconnect

Shared Memory

Parallel programming is the act of writing a program

for each node, co-ordinating the work that they do.

52
We assume each node has access to the integral
function and the integrand f.

First, we must initialize any required global variables,

i.e., the integration limits a, b as well as any optional
arguments such as tolerances, parameters, etc.

When parallel programs are executed, the number of

desired processors p is known in advance.

Suppose they are numbered from 1 to p.

Each processor can find out its local ID number.

We usually view one processor (say the first one) as

being the co-ordinating processor (or server); the rest
are workers (or clients).

The idea is that we will get each processor to compute

part of the integral, and then one processor can sum
up the partial results and report the answer.

We will also allocate an array Qp of dimension p in

global memory for the partial results and one space Q
for the final answer.

53
For concreteness, suppose we are trying to compute2
Z 1  
1 log(x)
Q= cos dx
0 x x

using p = 3 processors.

1 2 3

Cache Cache Cache

High−speed Interconnect

a=0, b=1, p=3, Qp(1:3), Q

2
This is the first problem from the SIAM 100-Digit Challenge.

54
 Static Scheduling

A natural way to partition (or schedule) the work

among the p processors is to assign

a+ih
b−a
Z
Qp(i) = f (x) dx, h := ,
a+(i−1)h p

to processor i, i = 1, 2, . . . , p.
The plan is for processor i to work on its part of Q
then write its result to Qp(i).
When all processors have finished, processor 1 sums
the entries of Qp and reports the final answer Q.
This schedule is called “static” because the assignment
of the tasks is fixed before the computation begins.
Because the time taken to solve the problem is
governed by the time it takes the slowest processor to
finish its task, this strategy works well if the workload
can be divided evenly among the processors.

55
Unfortunately, this is often not the case in quadrature
problems.

(If it were, adaptive quadrature would not be so

important!)

One way to mitigate this problem is to subdivide the

interval [a, b] into rp subintervals and assign r different
subintervals to each processor.

For example, if r = 5 and p = 3, processor 1 would be

assigned intervals 1, 4, 7, 10, 13, processor 2 would be
assigned intervals 2, 5, 8, 11, 14, etc.

The idea is that

R bthe processors would share the parts
of [a, b] where a f (x) dx is tough.

Clearly, this becomes more likely for large r.

But if r is too large, some intervals may be computed

“too accurately”!

Hence there is a tradeoff in the choice of r.

It seems that imposing a static scheduling on an

adaptive quadrature method may be ineffective.

56
 Dynamic Scheduling

In parallel computing, we strive to keep all processors

busy doing useful work all the time.

Such a computation is said to be load-balanced.

An idle processor is a wasted resource!

We would like a way for processors with easy

assignments to help out with the hard parts of the
integration when they are done.

An effective way to do this is instead of fixing the

assignment of the rp intervals, we allocate assignments
on a first-come, first-served basis.

For example, we begin by assigning interval i to

processor i for i = 1, 2, . . . , p.

Then interval p + 1 is assigned to whichever processor

is done its initial assignment first, interval p + 2 is
assigned to the first processor done its assignment
after that, etc.

57
This process continues until all rp intervals have been
assigned.

In this way, processors that have easy assignments and

are done quickly can return to the “assignment pool”
(or queue) to get a new assignment.

There is still no guarantee of good performance (e.g.,

interval rp turns out to be extremely difficult), but
it is likely that with some foresight extremely poor
performance can be avoided.

Care must be taken to ensure more than one processor

is not assigned the same job!

This means there may be idling while processors

determine which interval is to be handled next.

A more esoteric way to accomplish load balancing is to

have the roles of clients and servers be dynamic, so a
client with a subinterval that is too difficult can act as
the server for that subinterval.

58
 Example of Parallel Quadrature

We have already seen multi-threaded parallelism in

Matlab using vector operations.

The Parallel Computing Toolbox (PCT) allows us to

access multiple cores (or GPUs) on a single computer.

(To run Matlab on clusters of computers (or GPUs)

requires the Matlab Distributed Computing Server.)

The PCT is available and configured for us with up to

512 cores on [Link].

When running Matlab in parallel, the cores are known

as workers in a pool.

The worker pool can be opened automatically by

issuing certain Matlab commands, or it can be
started manually.

The default setting for the pool size on vlab is 8.

59
We now re-solve the two-dimensional quadrature
problem that we solved with the Monte Carlo method:

Z 1Z 1
1
dy dx,
0 0 1 − xy

with exact value π 2/6.

The basic strategy is to parallelize the for loops using

Matlab’s parfor command to allow each iteration
of the loops to run concurrently rather than serially.

See

parallelQuad.m

60