CSCI3340 - Scientific Computing
Lecture 8: Quadrature
Dr. Colin Tsang
1/31
Table of Contents
1 Introduction
▶ Introduction
▶ Newton-Cotes Rules
▶ Gauss Quadrature
▶ Adaptive Quadrature
▶ Further readings
2/31
Numerical Quadrature
1 Introduction
Quadrature means the approximate evaluation of a definite integral
∫ b
I= f(x)dx
a
From calculus, we know that a Riemann sum can be used to approximate a definite integral. To form a Riemann
sum, we partition the interval [a, b] into n intervals with a = x0 < x1 < ... < xn = b and approximate
∑
n
I ≈ In = wi f(δi )
i=1
with wi = xi − xi−1 and δi ∈ [xi−1 , xi ].
Every rule for numerical quadrature has the form of the above In . The key problem is to determine the nodes δi
(height) and the weight wi (width) to obtain the best possible result with the lowest computational cost.
3/31
Numerical Quadrature
1 Introduction
The value for the integral computed by a Riemann sum can also be interpreted as follows: we approximate the
function f by a piecewise constant function g
g(x) = f(δi ), xi−1 ≤ x ≤ xi , i = 1, ..., n
and compute
∫ b ∫ b
I= f(x)dx ≈ g(x)dx
a a
In fact, any other approximation of f by a function g which can be integrated analytically will lead to a quadrature
rule.
4/31
Table of Contents
2 Newton-Cotes Rules
▶ Introduction
▶ Newton-Cotes Rules
▶ Gauss Quadrature
▶ Adaptive Quadrature
▶ Further readings
5/31
Polynomial Interpolation
2 Newton-Cotes Rules
One obvious way to approximate f is by polynomial interpolation, i.e., our Lecture 4. Suppose that we have n + 1
pairs of data, (δi , yi ) for n = 0, ..., n. Then, we can represent the polynomial using Lagrange interpolation, which we
have discussed in Lecture 4. We have
∑n
f(x) = Pn (x) = li (x)yi
i=0
where
∏
n
x − δj
li (x) =
δi − δj
j̸=i
for i = 0, ..., n. We can obtain a quadrature rule by integrating both sides
∫ b ∫ b ∫ b∑
n n (∫
∑ b )
f(x)dx = Pn (x)dx = li (x)yi dx = li (x)dx yi
a a a i=0 i=0 a
6/31
Newton-Cotes Rules
2 Newton-Cotes Rules
We can define the weight as
∫ b
wi = li (x)dx
a
Then, we have the Newton-Contes Quadrature Rules
∫ b ∑
n
I= f(x)dx ≈ wi yi
a i=0
where yi = f(δi ). Note that
• The weights in a Newton–Cotes quadrature rule are given by the integrals of the corresponding Lagrange basis
polynomials.
• Newton–Cotes rules are typically constructed using equally spaced nodes (so as the interpolation problem),
although this is not a strict requirement.
• The nodes may, in principle, be arbitrarily chosen.
7/31
Some special cases
2 Newton-Cotes Rules
Let’s assume equidistance nodes in this slide.
For n = 1, we have the Trapezoidal Rule
∫ b
b−a
I= f(x)dx ≈ [f(a) + f(b)]
a 2
b−a
For n = 2, we set h = 2 . Then, we have the Simpson’s Rule:
∫ b
h
I= f(x)dx ≈ [f(a) + 4f(a + h) + f(b)]
a 3
b−a
For n = 3, we set h = 3 . Then, we have the Simpson’s 3-8 Rule:
∫ b
3h
I= f(x)dx ≈ [f(a) + 3f(a + h) + 3f(a + 2h) + f(b)]
a 8
b−a
For n = 5, we set h = 4 . Then, we have the Boole’s Rule:
∫ b
2h
I= f(x)dx ≈ [7f(a) + 32f(a + h) + 12f(a + 2h) + 32f(a + 3h) + 7f(b)]
a 45
8/31
Order of a Quadrature Formula
2 Newton-Cotes Rules
A quadrature formula with nodes δi and weights wi is of order p, if it computes the exact integral value for all
polynomials pk up to degree k ≤ p − 1, i.e.,
∫ b ∑
n
pk (x)dx = wi pk (δi )
a i=0
for k = 0, ..., p − 1.
We have:
• Trapezoidal Rule: Order = 2
• Simpson’s Rule: Order = 4
• Simpson’s 3-8 Rule: Order = 4
• Boole’s Rule: Order = 6
9/31
Order of the Trapezoidal Rule
2 Newton-Cotes Rules
For degree 0, i.e., p0 (x) = C for some constant C, the Trapezoidal rule is equal to the exact integral
∫
b−a b−a b
(f(a) + f(b)) = (C + C) = C(b − a) = C[x]ba = p0 (x)dx
2 2 a
Next, for degree 1, i.e., p1 (x) = mx + c, the Trapezoidal rule is equal to the exact integral
∫
b−a b−a m(b2 − a2 ) m b
(f(a) + f(b)) = (ma + c + mb + c) = + c(b − a) = [x2 ]ba + c[x]ba = p1 (x)dx
2 2 2 2 a
∫b b3 −a3
For degree 2, we can let f(x) = x2 . The exact integral is a f(x)dx = 3 . However, the Trapezoidal rule gives
b−a b−a 2 b3 − a3
(f(a) + f(b)) = (a + b2 ) ̸=
2 2 3
Hence, the Trapezoidal rule is equal to the exact integral of pk (x) for k = 0, 1. Therefore, the Trapezoidal rule is
order 2.
10/31
Order of a Quadrature Formula
2 Newton-Cotes Rules
Another way to find the order is to consider the monomials, i.e., 1, x, x2 ,...,xm .
Theorem (Order by Monomials)
If a quadrature rule is exact for f(x) = xk for k = 0, ..., m, i.e.,
∫ b ∑
n
k
x dx = wi δik
a i=0
then it has an order of m + 1.
11/31
Order of a Quadrature Formula
2 Newton-Cotes Rules
Proof.
Any degree k polynomial can be written as pk (x) = a0 + a1 x + ... + ak xk , for k = 0, ..., m. Thus, we have
∫ b ∫ b
pk (x)dx = (a0 + a1 x + ... + ak xk )dx
a a
∫ b ∫ b ∫ b
= a0 1dx + a1 xdx + ... + ak xk dx
a a a
∑
n ∑
n ∑
n
= a0 wi + a1 wi δi + ... + ak wi δik (By assumption)
i=0 i=0 i=0
∑
n
= wi (a0 + a1 δi + ... + ak δik )
i=0
∑
n
= wi pk (δi ) (for k = 0, ..., m)
i=0
Thus, the quadrature rule is exact for pk (x), where k = 0, ..., m. Thus, by definition, it has the order m + 1.
12/31
Minimum Order of the Newton-Cotes Rule
2 Newton-Cotes Rules
Theorem (Minimum Order of the Newton-Cotes Rule)
If the Newton-Cotes Rule has n-nodes, then the order of it must be at least n.
Proof.
Consider the Newton-Cotes Rules with n-nodes. We have
∑n n (∫ b
∑ ) ∫ b∑
n ∫ b
wi f(δi ) = li (x)dx f(δi ) = li (x)f(δi )dx ≈ f(x)dx
i=1 i=1 a a i=1 a
∫b
The weights wi = a li (x)dx is coming from the fact that we are using Lagrange interpolation with n nodes to
approximate f(x). However, if f(x) is a polynomial of degree ≤ n − 1, then the Lagrange interpolation with n nodes
will be exact. Thus, the Newton-Cotes Rule with n nodes is exact for all polynomial up to degree n − 1. So, it has
order n.
13/31
Error Analysis
2 Newton-Cotes Rules
• In general, it is impossible to know the exact error of a quadrature rule for an arbitrary function.
• Unless we use infinitely many nodes—which is impossible in practice—the rule may miss features like narrow
peaks that contribute significantly to the integral.
• Nevertheless, the order of a quadrature formula provides useful information on its accuracy.
• Any sufficiently smooth function f(x) can be approximated by a polynomial using its Taylor expansion.
• If the quadrature rule is exact for all polynomials up to degree p − 1, the error depends on the p-th and higher
derivatives of f(x).
• Thus, a higher order p generally leads to a smaller error.
14/31
Table of Contents
3 Gauss Quadrature
▶ Introduction
▶ Newton-Cotes Rules
▶ Gauss Quadrature
▶ Adaptive Quadrature
▶ Further readings
15/31
Order Bound
3 Gauss Quadrature
Theorem (Order Bound)
For any quadrature rule with n nodes,
∫ b ∑
n
I= f(x)dx ≈ wi f(δi )
a i=1
the order of it is bounded by 2n, i.e., an n-point quadrature rule can exactly integrate a polynomial of degree 2n − 1
at most.
16/31
Order Bound
3 Gauss Quadrature
Proof.
If there exists a polynomial of degree 2n such that the quadrature rule cannot integrate it exactly, then the proof is
complete. Consider a polynomial with the nodes as roots
∏
n
Q(x) = (x − δi )2
i=1
The degree of Q is 2n. Since Q(x) > 0 for any x ∈ [a, b], we have
∫ b
Q(x)dx > 0
a
However, the n-point quadrature rule will give us
∑
n ∑
n
wi f(δi ) = wi Q(δi ) = 0
i=1 i=1
Hence, the order of any n-point quadrature is at most 2n.
17/31
Positivity of Weights
3 Gauss Quadrature
There is another very desirable property for a high order quadrature rule. If a quadrature has a very high order
(2n − 1 or 2n), their weights are positive.
Theorem (Positivity of Weights)
If an n-nodes quadrature rule
∫ b ∑
n
I= f(x)dx ≈ wi f(δi )
a i=1
is of order p ≥ 2n − 1, then the weights wi of the Newton-Cotes Rules, for i = 1, ..., n, are positive.
18/31
Positivity of Weights
3 Gauss Quadrature
Proof.
Consider the Lagrange polynomials
∏
n
x − δj
lk (x) =
δk − δj
j̸=k
for k = 1, ..., n. The degree is n − 1. Thus, l2k (x) is of degree 2n − 2. If the order is p ≥ 2n − 1, then the quadrature
rule is exact for any polynomial of degree less than 2n − 2. Hence, we have
∑
n ∫ b
wi l2k (δi ) = l2k (x)dx
i=1 a
Since∑lk (δi ) = 0 for any i ̸= k and lk (δk ) = 1, all terms on the LHS will vanish except i = k. Thus, we
have ni=1 wi l2k (δi ) = wk l2k (δk ) = wk (1)2 = wk . Hence, we have
∫ b
wk = l2k (x)dx > 0
a
for k = 1, ..., n.
19/31
From [a, b] to [−1, 1]
3 Gauss Quadrature
Consider a change of variable
b−a a+b
x= t+
2 2
Hence, we have
2x − b − a
t=
b−a
b−a = −1 when x = a, and t = b−a = 1 when x = b. Therefore, we can always
This means that t = 2a−b−a 2b−b−a
transform a general finite interval [a, b] to [−1, 1] because
∫ b ∫
b−a 1 b−a a+b
f(x)dx = f( t+ )dt
a 2 −1 2 2
WLOG, consider an integral on the interval [−1, 1] with a quadrature rule
∫ 1 ∑
n
I= f(x)dx ≈ wi f(δi )
−1 i=1
For a given n, we want to find wi and δi such that the order is as high as possible.
20/31
Gauss-Legendre Rule
3 Gauss Quadrature
∫1
Suppose that n = 2. We have a quadrature rule −1 f(x)dx ≈ w1 f(δ1 ) + w2 f(δ2 ). From the Order Bound theorem, the
highest possible order is 4. Hence, we want to satisfy four equations:
∫ 1
w1 δ1k + w2 δ2k = xk dx , for k = 0, 1, 2, 3
−1
Hence, we have a system of nonlinear equations
w1 + w2 = 2
w δ + w δ = 0
1 1 2 2
w1 δ1 + w2 δ22 = 23
2
w δ 3 + w δ 3 = 0
1 1 2 2
We have learned how to solve such a system in Lecture 3. The solution is given by
−1 1
{w1 = 1, w2 = 1, δ1 = √ , δ2 = √ }
3 3
Hence, we have the Gauss-Legendre Rule for n = 2
∫ 1
−1 1
f(x)dx ≈ f( √ ) + f( √ )
−1 3 3
21/31
Gauss-Legendre Rule
3 Gauss Quadrature
For n = 3, we can obtain 6 equations with 6 unknowns following the same method. Once can compute that the
corresponding quadrature rule is given by
∫ 1 √ √
5 3 8 5 3
f(x)dx = f(− ) + f(0) + f( )
−1 9 5 9 9 5
The nonlinear system becomes very hard to solve when n is large.
22/31
Gauss Quadrature
3 Gauss Quadrature
Definition: If qn (x) is a degree n polynomial such that
∫ b
pk (x)qn (x)dx = 0
a
for any polynomial pk (x) with degree k ≤ n − 1, then we say that qn (x) is an orthogonal polynomial on the interval
[a, b].
Theorem (Gauss Quadrature)
If an n-node quadrature rule satisfies the following conditions:
• The nodes δi are the roots of Qn (x), an orthogonal polynomial on the interval [−1, 1].
∫1
• The weights wi is obtained by the Newton-Cotes Rules, i.e., wi = −1 li (x)dx.
Then, it is named the Gauss Quadrature. It has an order of 2n, which is the maximum possible order.
23/31
Gauss Quadrature
3 Gauss Quadrature
Proof.
For any polynomial with degree 2n − 1, divide it by Qn (x) and then integrate both sides
∫ 1 ∫ 1 ∫ 1
P2n−1 (x) = Hn−1 (x)Qn (x) + Rn−1 (x)
−1 −1 −1
where Hn−1 (x) with degree n − 1 is the quotient and Rn−1 (x) with degree n − 1 is the remainder. The error of our
n-nodes quadrature rule is given by
∫ 1 ∑
n ∫ 1 ∑
n ∫ 1 ∑
n
E(n) = P2n−1 (x)− wi P2n−1 (δi ) = Hn−1 (x)Qn (x)− wi Hn−1 (δi )Qn (δi )+ Rn−1 (x)− wi Rn−1 (δi )
−1 i=1 −1 i=1 −1 i=1
• Since Qn (x) is orthogonal, the first term vanishes.
• Since δi are roots of Qn (x), the second term vanishes.
• By the minimum order of the Newton-Cotes Rule theorem, our n-nodes quadrature rule has an order of at least
n. As the remainder Rn−1 (x) is of degree n − 1, the last two terms also vanish.
Thus, the error is zero for any polynomial of degree 2n − 1.
24/31
Orthogonal Polynomials
3 Gauss Quadrature
There is one remaining problem: How to find the orthogonal polynomial Qn (x)?
• For Gauss quadrature on [−1, 1], the standard orthogonal polynomials are Legendre polynomials Qn (x).
• Qn (x) can be computed with the recursion formula:
(n + 1)Qn+1 (x) = (2n + 1)xQn (x) − nQn−1 (x)
with Q0 (x) = 1, Q1 (x) = x.
• For small n:
1 1
Q2 (x) = (3x2 − 1), Q3 (x) = (5x3 − 3x)
2 2
• Keep in mind that this formula only works on [−1, 1].
• The proof of this recursion formula will be skipped in this course. You can read Chapter 9.3.2 of our textbook if
you are interested.
25/31
Gauss Quadrature
3 Gauss Quadrature
∫b
Here is an algorithm to apply the Gauss Quadrature to approximate an integral a f(x)dx with the maximum possible
order:
1. Change the interval to [−1, 1] by substitution.
2. Find Qn (x) by (n + 1)Qn+1 (x) = (2n + 1)xQn (x) − nQn−1 (x) on [−1, 1].
3. Solve Qn (x) = 0 for x ∈ [−1, 1]. Take the roots as δi .
∫1
4. Compute the weights by wi = −1 li (x)dx, where li (x) is the Lagrange basis polynomial for the nodes.
∑
5. The Gauss Quadrature is given by ni=1 wi f(δi )
26/31
Table of Contents
4 Adaptive Quadrature
▶ Introduction
▶ Newton-Cotes Rules
▶ Gauss Quadrature
▶ Adaptive Quadrature
▶ Further readings
27/31
Adaptive nodes
4 Adaptive Quadrature
Key idea: We want to have more nodes in the “difficult” interval and less nodes in the “easy” interval.
1. Compute a rough approximation I1 (e.g., Simpson’s rule) on [a, b].
a+b
2. Split [a, b] at midpoint m = 2 .
3. Compute I2 = Simpson’s rule on [a, m] plus on [m, b].
4. If |I2 − I1 | is small enough, accept I2 .
5. Otherwise, recurse on subintervals [a, m] and [m, b].
28/31
Example
4 Adaptive Quadrature
∫1 1
Consider 0 1+100(1−x)2 dx with the Adaptive Simpson Rule.
29/31
Table of Contents
5 Further readings
▶ Introduction
▶ Newton-Cotes Rules
▶ Gauss Quadrature
▶ Adaptive Quadrature
▶ Further readings
30/31
Further readings
5 Further readings
More details and examples can be found in CHAPTER 9 of this reference book:
• Gander, Walter, Martin J. Gander, and Felix Kwok. Scientific computing-An introduction using Maple and
MATLAB. Vol. 11. Springer Science & Business, 2014
31/31