APPLIED NUMERICAL METHODS

SEMM/SKMM3023

Numerical Integration

Department of Aeronautics, Automotive, & Ocean Engineering

School of Mechanical Engineering, Faculty of Engineering,
Universiti Teknologi Malaysia

1
 Why numerical integration?

■ Many engineering application requires the evaluation of the integral:

𝑏
𝐼 = න 𝑓 𝑥 𝑑𝑥
𝑎
■ Compared to differentiation, integration is more difficult to solve and
often, analytical methods may not exist.
■ Hence, the practical applications of numerical methods for integration
is actually quite common.

2
 Why numerical integration?

■ Numerical integration is an approximation of definite integrals and

used in cases where:
1. f(x) is a complicated continuous function that is difficult or impossible
to integrate in closed form
2. f(x) is not known in analytical form; it may be known in tabular form
involving discrete numbers.
3. Limits of integration may be infinite
4. f(x) may be discontinuous or may become infinite at some point in
range a ≤ x ≤ b

3
 Why numerical integration?

■ When you integrate a function, what do you get?

4
 Numerical Integration

■ In this chapter we will learn four types of numerical integration

methods:

1. Trapezoid rule
2. Simpson’s 1/3rd rule Newton-Cotes Formulas
3. Simpson’s 3/8th rule
4. Gauss quadrature

5
 Trapezoid Rule

■ Consider the following function f(x) in the figure below. The area under
the curve can be calculated by dividing the curve into several equal
segments and approximating the area for each segment as a trapezoid

1
𝐼1 = (𝑥1 − 𝑥0 )(𝑓0 + 𝑓1 )
2
1
𝐼1 = ℎ(𝑓0 + 𝑓1 )
2
1
𝐼2 = ℎ(𝑓1 + 𝑓2 )
2
1
𝐼3 = ℎ(𝑓2 + 𝑓3 )
2

6
 Trapezoid Rule

■ In general, the integral using the trapezoid rule can be expressed as:
𝑏 𝑥𝑛
1
𝐼 = න 𝑓 𝑥 𝑑𝑥 = න 𝑓 𝑥 𝑑𝑥 ≈ ℎ 𝑓0 + 2𝑓1 + 2𝑓2 + 2𝑓3 + ⋯ 𝑓𝑛
𝑎 𝑥0 2
𝑛−1
1
≈ ℎ 𝑓0 + 𝑓𝑛 + 2 ෍ 𝑓𝑖
2
𝑖=1

■ This is also known as the Composite Trapezoid Rule

■ The accuracy of the trapezoid rule increases with decreasing step size, h.
■ In some cases, the trapezoid rule is also used for discrete data. However,
the step size may not be uniform in this case.

7
 In MATLAB, there is a function
Trapezoid Rule called ‘trapz’ which takes in input of
x and f(x) data and calculates the
integral using the trapezoid rule!

■ Example: Evaluate the following integral using the trapezoid rule and a
step size of h = 0.1. n = (xn-x0)/h
1.6
2
𝐼 = න 𝑒 𝑥 𝑑𝑥
1 n is number of segment

x0 x1 x2 x3 x4 x5 x6 = xn
x 1 1.1 1.2 1.3 1.4 1.5 1.6
f(x) 2.7183 3.3535 4.2207 5.4195 7.0993 9.4877 12.9358

𝐼
1
≈ (0.1)ሾ2.7183 + 2 3.3535 + 4.2207 + 5.4195 + 7.0993 + 9.4877
2
+ 12.9358ሿ ≈

8
 Simpson’s 1/3rd Rule

■ Previously in the trapezoid rule, a linear line was used to connect two
discrete points. The Simpson’s 1/3rd rule tries to improve that by
connecting three discrete points with a second order polynomial
(quadratic) equation instead. The general form of equation is:

𝑓 𝑥 = 𝑎2 𝑥 2 + 𝑎1 𝑥 + 𝑎0

9
 Simpson’s 1/3rd Rule

■ Constants a0, a1 and a2 can be determined by making the approximating

polynomial pass through three consecutive points (xi-1, fi-1),(xi, fi) and
(xi+1, fi+1). This means that the coefficients can be solved using the
collocation polynomial method you have learnt before!
■ Lets consider xi at the origin location (xi = 0) with a uniform step size of
h. This means that xi-1 = -h and xi+1 = h. Applying the collocation
polynomial method:

10
 Simpson’s 1/3rd Rule

■ Hence, we get:

■ The integral of a second order polynomial between xi-1 and xi+1 is:

11
 Simpson’s 1/3rd Rule

■ Substituting the coefficient terms, we finally get

■ The term ‘1/3’ in Simpson’s one-third rule refers to the factor of 1/3.
■ For multistage application, we need to divide the range a ≤ x ≤ b into n
segments of equal width h = (b−a) / n , with n being an EVEN number so
that it can be applied for groups of two segments.

12
 Simpson’s 1/3rd Rule

■ The general equation is thus:

13
 Simpson’s 1/3rd Rule

■ Example: Evaluate the following integral using the Simpson’s 1/3rd rule
and a step size of h = 0.1.
1.6
2
𝐼 = න 𝑒 𝑥 𝑑𝑥
1

Check if this is possible first. Remember, n must be an even number.

𝑏 − 𝑎 1.6 − 1
𝑛= = = 6 → 𝑒𝑣𝑒𝑛 𝑛𝑢𝑚𝑏𝑒𝑟
ℎ 0.1

So this means that the problem can be solve using Simpson’s 1/3rd rule.

14
 Simpson’s 1/3rd Rule

1.6
2
𝐼 = න 𝑒 𝑥 𝑑𝑥
1

x0 x1 x2 x3 x4 x5 x6
x 1 1.1 1.2 1.3 1.4 1.5 1.6
f(x) 2.7183 3.3535 4.2207 5.4195 7.0993 9.4877 12.9358

𝐼
1
≈ (0.1)ሾ2.7183 + 4 3.3535 + 5.4195 + 9.4877 + 2 4.2207 + 7.0993
3
+ 12.9358ሿ ≈

15
 Simpson’s 3/8th Rule

■ Simpson’s 3/8th rule aims to further improve the approximation by

fitting a third order polynomial (cubic) equation through four data
points instead. The general form of equation is:

𝑓 𝑥 = 𝑎3 𝑥 3 + 𝑎2 𝑥 2 + 𝑎1 𝑥 + 𝑎0

■ Again, the collocation polynomial method is used to determine the

unknown coefficients of the cubic equation

16
 Simpson’s 3/8th Rule

■ Consider four consecutive points (xi-1, fi-1),(xi, fi), (xi+1, fi+1) and (xi+2, fi+2).
Again, assume that xi is at the origin location (xi = 0) with a uniform step
size of h. This means that xi-1 = -h, xi+1 = h and xi+2 = 2h. Applying the
collocation polynomial method:

17
 Simpson’s 3/8th Rule

■ Solving for the four unknow coefficients we get:

18
 Simpson’s 3/8th Rule

■ The integral of a third order polynomial between xi-1 and xi+2 is:

19
 Simpson’s 3/8th Rule

■ Substituting in the expressions for the unknown coefficients results in:

■ The term ‘3/8’ in Simpson’s three-eighth rule refers to the factor of 3/8.

20
 Simpson’s 3/8th Rule

■ The general equation is thus:

■ The number of segments n must be a multiple of 3 so that can be applied

for groups of three segments.

21
 Simpson’s 3/8th Rule

■ Example: Evaluate the following integral using the Simpson’s 3/8th rule
and a step size of h = 0.1.
1.6
2
𝐼 = න 𝑒 𝑥 𝑑𝑥
1

Check if this is possible first. Remember, n must be a multiple of 3.

𝑏 − 𝑎 1.6 − 1
𝑛= = = 6 → 𝑚𝑢𝑙𝑡𝑖𝑝𝑙𝑒 𝑜𝑓 3
ℎ 0.1

So this means that the problem can be solve using Simpson’s 3/8th rule.

22
 Simpson’s 3/8th Rule

1.6
2
𝐼 = න 𝑒 𝑥 𝑑𝑥
1

x0 x1 x2 x3 x4 x5 x6
x 1 1.1 1.2 1.3 1.4 1.5 1.6
f(x) 2.7183 3.3535 4.2207 5.4195 7.0993 9.4877 12.9358

𝐼
3
≈ (0.1)ሾ2.7183 + 3 3.3535 + 4.2207 + 7.0993 + 9.4877 + 2 5.4195
8
+ 12.9358ሿ ≈

23
 Gauss Quadrature

■ Consider the case where the trapezoidal rule is used to calculate the
area under a straight line which connects the function values at the ends
of the integration intervals is

where a and b are limits of integration and (b − a) is the integration

interval.

24
 Gauss Quadrature

■ Because in this case the trapezoidal rule must pass through end points,
there are cases where the formula results in large error

25
 Gauss Quadrature

■ Suppose we take away the constraint of fixed end points and evaluate
the area under a straight line joining any two points on the curve. Then
we could define a straight line that would balance the positive and
negative errors as shown in the figure below.

26
 Gauss Quadrature

■ The domain of integration for this rule is conventionally taken as [−1,1]

to simplify the mathematics and make formulation as general as
possible. A change of variable is used to translate the original [a, b]
integration limits into [−1,1] form.

■ where wi are unknown coefficient also called weights, xi are the

transformed values of xd where xd are called Gauss points at which
integrand is evaluated.

27
Gauss Quadrature

𝑏−𝑎
𝑑𝑥 = 𝑑𝑥𝑑
2

28
 Gauss Quadrature

■ Hence, the transformation is:

1
𝑏−𝑎 𝑏 + 𝑎 + 𝑏 − 𝑎 𝑥𝑑
𝐼= න𝑓 𝑑𝑥𝑑
2 2
−1

■ Notice that the x in f(x) has been transformed as so:

𝑏 + 𝑎 + 𝑏 − 𝑎 𝑥𝑑
𝑥𝑖 =
2

29
 Gauss Quadrature

■ Where the transformed function is usually written as

1 1
𝑏−𝑎 𝑏 + 𝑎 + 𝑏 − 𝑎 𝑥𝑑 𝑏−𝑎
න𝑓 𝑑𝑥𝑑 = න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑
2 2 2
−1 −1

■ Once we have transformed the integral into the form of:

𝑏 1
𝑏−𝑎
න 𝑓 𝑥 𝑑𝑥 = න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑
2
𝑎 −1

30
 Gauss Quadrature

■ We can focus on the right-hand side of the equation:

1 𝑛

න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑 ≈ ෍ 𝑤𝑖 𝑓 𝑥𝑖
−1 𝑖=1
1 𝑛
𝑏−𝑎 𝑏−𝑎
න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑 ≈ ෍ 𝑤𝑖 𝑓 𝑥𝑖
2 2
−1 𝑖=1

■ Again, 𝑥𝑑 are known as the Gauss points.

31
 Gauss Legendre Quadrature

■ Consider the case where n = 2 for the right-hand side equation

■ For the case of n = 2, the formula should give exact value for third order
polynomials (cubic) and below. Hence, we extend this reasoning by
assuming that it fits the integrals of f (x) = 1, x, x2, x3 functions.

32
 Gauss Legendre Quadrature

■ Evaluation of integrals for f (x) = 1, x, x2, x3 equations:

33
 Gauss Legendre Quadrature

■ Since limits of integration, −1 and +1, are symmetric about x = 0, we

expect x1 and x2 also to be symmetric about x = 0.
■ By setting x2 = −x1, we obtain the following w1 = w2 = 1
■ Once we get w1 and w2 , we need to determine x2 and x1. Referring to the
previous slide, we see that for f (x) = x and x3, the solutions are zero.
This would lead to an infinite solution for x2 and x1. For the case where
f(x) = 1, the equation is independent of x2 and x1. Hence, we need to
solve using the equation f(x) = x2.
2 2
2 2 2
2 2
2
𝑤1 𝑥1 + 𝑤2 𝑥2 = → 𝑥1 + 𝑥2 = → 2𝑥1 = 𝑥1 and 𝑥2 are
3 3 3
the Gauss
1 1 1 points (𝑥𝑑 )
𝑥12 = → 𝑥1 = ,𝑥 = − 𝑜𝑟 𝑣𝑖𝑐𝑒 𝑣𝑒𝑟𝑠𝑎
3 3 2 3

34
 Gauss Legendre Quadrature

■ For higher number of n points, the gauss quadrature formula can be

expressed as the following

■ where n = number of points. Values for w’s and x’s for up to and
including the five-point formula are summarized in the next two slides.

35
Gauss Legendre Quadrature

36
Gauss Legendre Quadrature
𝑥𝑑

37
 Gauss Legendre Quadrature

■ Example:

Use the Gauss legendre quadrature with n = 3

0.8 1 1
𝑏−𝑎
𝐼 = න 𝑓 𝑥 𝑑𝑥 = න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑 = 0.4 න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑
2
0 −1 −1
1 𝑛

න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑 ≈ ෍ 𝑤𝑖 𝑓 𝑥𝑖
−1 𝑖=1

𝑏 + 𝑎 + 𝑏 − 𝑎 𝑥𝑑
𝑥𝑖 = → 𝑥𝑖 = 0.4 + 0.4𝑥𝑑
2

38
 Gauss Legendre Quadrature

■ Example:

Use the Gauss legendre quadrature with n = 3

𝑥1 = 0.4 + 0.4 − 3/5 = 0.090161332

𝑓 𝑥1 = 1.265855804
𝑥2 = 0.4 + 0.4 0 = 0.4
𝑓 𝑥2 = 2.4560
𝑥3 = 0.4 + 0.4 3/5 = 0.709838667
𝑓 𝑥3 = 2.1869442

39
 Gauss Legendre Quadrature

■ Example:

Use the Gauss legendre quadrature with n = 3

1 𝑛=3

0.4 න 𝑓𝑑 𝑥𝑑 𝑑𝑥𝑑 ≈ 0.4 ෍ 𝑤𝑖 𝑓 𝑥𝑖 = 0.4ሾ𝑤1 𝑓 𝑥1 + 𝑤2 𝑓 𝑥2 + 𝑤3 𝑓 𝑥3 ሿ

−1 𝑖=1
𝑛=3
5 8 5
0.4 ෍ 𝑤𝑖 𝑓 𝑥𝑖 = 0.4 1.265855804 + 2.4560 + 2.1869442
9 9 9
𝑖=1
𝑛=3

0.4 ෍ 𝑤𝑖 𝑓 𝑥𝑖 = 1.64053
𝑖=1

40
Thank you

41