Universiti Teknologi PETRONAS

MDB 3053
Numerical Methods

NUMERICAL INTEGRATION

Reading Assignment: Chapter 21

1
 Lesson Outcomes

By the end of the lesson, students should be able

1. To perform numerical integration using Newton-

Cotes Integration Formulas

2. To analyze errors associated with Newton-

Cotes Integration Formulas

Chap 21/2

2
 NUMERICAL INTEGRATION

• Calculus is the mathematics of change. Calculus is an essential

tool for engineers as we deal with mechanical systems and
plant processes that change.

• Standing in the heart of calculus are the mathematical

concepts of integration and differentiation

• Integration is analogous to the process of summation.

b b
I   f ( x)dx   f n ( x)dx
a a

f n ( x)  a0  a1 x    an 1 x n 1  an x n
3

3
 NUMERICAL INTEGRATION

• For simple function; we have closed-formed analytical

integration. But for complex function, it is impossible to
obtain analytical expression.
𝑥 𝑛+1 𝑥 2
‫׬‬ 𝑥 𝑛 𝑑𝑥 = easy! 𝑏𝑢𝑡 න 𝑒 𝑑𝑥 ℎ𝑎𝑠 𝑛𝑜 𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐𝑎𝑙 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛!
𝑛+1

• Newton-Cotes Integration formulas are one of the common

numerical integration methods. Strategy here is to replace
complex function with a discrete polynomial function.

• 2 types of Newton-Cotes Integration Schemes

• Trapezoidal Rule
• Simpson’s Rule 3

4
 TRAPEZOIDAL RULE
• Trapezoidal Rule is the simplest numerical integration
schemes

• Suppose we have a function f(x) and we want to find the

area below the curve bound by a and b, as shown.

• The simplest way is a straight line and the area below the
straight line is approximated by a trapezoidal segment.

f(b)
Trapezoidal Rule
f(a)
f (a)  f (b)
I  (b  a)
2
a b

5
 ERROR In TRAPEZOIDAL RULE
• When we approximate the integral using a trapezoidal, the
error may be substantial:

(b  a )3
Approximate Error: Ea   f (x )
12
x lies somewhere in the interval from a to b
where f(x) = average value of the second derivative
Error

6
 EXAMPLE
Use Trapezoidal Rule to compute the integral I of the f(x)
function f(x) from a = 0 to b =0.8. Find errors Et, Ea.

Exact solution: I = 1.640533

f (a)  f (b)
I  (b  a)
2
(b  a )3
Ea   f (x )
12

7
 SOLUTION

0 .8
I  f ( x)dx
0 Exact solution: I = 1.640533

Trapezoidal rule: f (a)  f (b) (b  a )3

I  (b  a) Ea   f (x )
2 12

Error:

(huge error!)

8
 cont.

0 .8
I  f ( x)dx
0
Exact I = 1.640533

(b  a )3
Approximate Error: Ea   f (x )
12

where f(x) = average value of the second derivative

Approximate Error:
Comment: Discrepancy exists between Et and Ea due to less accuracy of Trapezoidal Rule 9
 CLASS ACTIVITY
Use Trapezoidal Rule to compute the integral of
f ( x ) = x 5 + 2x + 4 from a = 1 to b = 2. Calculate
the errors: Et Ea . Exact I = 1.64053
Ans: I = 23.5, Ea= 6, Ea= 6.25

f (a)  f (b)
I  (b  a)
2
(b  a )3
Ea   f (x )
12

10
MULTIPLE APPLICATION Trapezoidal Rule
• Is there any way we could improve the Trapezoidal Rule?
Yes, by using Multiple Application Trapezoidal Rule
• Area can be divided into multiple n segments:

Repeat f ( x0 )  f ( x1 ) f ( x1 )  f ( x2 ) f ( xn1 )  f ( xn )
n times I  h 2
h
2
 h
2
h n 1

I   f ( x0 )  2 f ( xi )  f ( xn ) where h  b  a a  x0 b  xn
2 i 1  n
start point intermediate end point 11
 MULTIPLE For n segments
APPLICATION
TRAPEZOIDAL RULE

ba
interval, h 
n
where a  x0 b  xn
Repeat n times applying
Trapezoidal rule

Error is summed up as:

 f (xi)  nf 
(b  a) 3
Ea   2
f 
12n
interval

12
 EXAMPLE
Use the 2-segment trapezoidal rule to estimate the integral of

from a = 0 to b = 0.8. Exact I = 1.640533 Ans: I = 1.0688, Ea = 0.64

b  a 0.8  0 h n 1

h
n

2
 0.4 I   f ( x0 )  2 f ( xi )  f ( xn )
2 i 1 
(b  a ) 3
3 points: 0 0.4 0.8 Ea   2
f 
12n
2 segments

13
 **CLASS ACTIVITY
Use the 4-segment trapezoidal rule to estimate the integral of

from a = 0 to b = 0.8. Exact I = 1.640533 Ans: I = 1.4848, Ea = 0.16

h n 1

I   f ( x0 )  2 f ( xi )  f ( xn )
2 i 1 
ba (b  a ) 3
h Ea   2
f 
n 12n

14
 SOLUTION
Use the 4-segment trapezoidal rule to estimate the integral of

from a = 0 to b = 0.8. Exact I = 1.640533 Ans: I = 1.4848, Ea = 0.16

ba
h
n h n 1

0.8 0 I   f ( x0 )  2 f ( xi )  f ( xn )
h  0.2 2 i 1 
4
f (0)  0.2

0.2
0.2  2(1.288  2.456  3.464)  0.232
2
f (0.2) 1.288 1.4848
f (0.4)  2.456 Et = 1.640533 – 1.4848 = 0.1557
f (0.6)  3.464
f (0.8)  0.232
(b  a)3 0.83
Ea   2
f    2
(60)  0.16
12n 12(4)
We are getting closer to exact answer, Et  Ea 15
 Can we improve Trapezoidal Rule?

• can improve by sub-dividing the area

into smaller segments  Multiple
Application Trapezoidal Rule.

OR

• can improve accuracy by replacing

the straight line with a curve 
Simpson’s Rule

16
 SIMPSON’S RULE
• Simpson’s Rule is a numerical integration method of more
accurate estimate of an integral by using higher-order
polynomials to connect the points (at least 3 points).

• Types of Simpson’s Rule

• Simpson’s 1/3 Rule → uses 3 points (even segment no.)
• Simpson’s 5/8 Rule → uses 4 points

Simpson’s 1/3 Rule Simpson’s 3/8 Rule

17
 SIMPSON’S 1/3 RULE
• Simpson’s 1/3 rule results when a second-order Lagrange
interpolating polynomial is used → uses 3 points

• More accurate than Trapezoidal rule because it uses 3 points

b b
I   f ( x)dx   f 2 ( x)dx
a a

2nd order Lagrange Interpolation

a  x0 b  x2
 ( x  x1 )( x  x2 ) ( x  x0 )( x  x2 ) ( x  x0 )( x  x1 ) 
x2

I   f ( x0 )  f ( x1 )  f ( x2 )dx
x0 
( x0  x1 )( x0  x2 ) ( x1  x0 )( x1  x2 ) ( x2  x0 )( x2  x1 ) 

ba
Simpson' s 1/3 Rule : I 
h
 f ( x0 )  4 f ( x1 )  f ( x2 ) h
3 2

h5 ( 4)
• Approximated Error, Ea   f (x ) a x b
90
18
 EXAMPLE
Use Simpson’s 1/3 rule to estimate the integral I
of f (x)  0.2  25x  200x2  675x3  900x4  400x5
from a = 0 to b = 0.8. Find errors Et, Ea Exact I = 1.640533

ba
I
h
 f ( x0 )  4 f ( x1 )  f ( x2 ) h
3 points: 0 0.4 0.8 3 2
1 segment h5 ( 4)
Ea   f (x ) a x b
90
0.8
h  0.4
2
f (0)  0.2 I
0.4
0.2  4(2.456)  0.232  1.367467
3
f (0.4)  2.456
Et  |1.640533 1.367467 | 0.2703667
f (0.8)  0.232
 t 16.6%
h5 ( 4) 0.45
Ea   f (x )   (2400)  0.2730667
90 90
Comment: more accurate than 2-segment Trapezoidal Rule because it is up to third order 19
 CLASS ACTIVITY
Use a single segment Simpson’s 1/3 rule to compute the integral of
ƒ ( x)= x 5 + 2x + 4 from a = 1 to b = 2
and calculate the errors: Et , Ea Exact I = 17.5

Ans: I = 17.5625, Et = 0.0625, Ea = 0.0625

ba
3 points: 1 1.5 2 I
h
 f ( x0 )  4 f ( x1 )  f ( x2 ) h
3 2
h5 ( 4)
Ea   f (x ) a x b
90

20
 Multiple-Application of Simpson’s 1/3 Rule
• Just as the Trapezoidal rule, Simpson’s rule can be improved
by dividing the integration interval into multiple segments of
equal width.

• More accurate than Trapezoidal. However it is limited to

cases where values are equispaced.

• Further, it is limited to situations where there are an even

number of segments and odd number of points.

(b  a)  
 f  x0   4  f  xi   2  f x j  f xn 
n1 n2
I
3n  i1,3,5... j2,4,6... 
(b  a)5 ( 4)
Approximate Error Ea   4
f (x ) a x b
180n
21
Multiple-Application of Simpson’s 1/3 Rule

• Method is applicable if the number of segments is even,

i.e., 2, 4, 6, 8, 10, etc

# of segment = 10 (EVEN)

# of points = 11 (ODD)

22
 CLASS ACTIVITY
Use 2 segments of Simpson’s 1/3 rule to estimate the integral of
f (x)  0.2  25x  200x2  675x3  900x4  400x5

from a = 0 to b = 0.8. Use n = 4

Exact I = 1.64053

5 points: 0 0.2 0.4 0.6 0.8

2 segments, h = 0.2

23
 SOLUTION
Use 2 segments of Simpson’s 1/3 rule to estimate the integral of
f (x)  0.2  25x  200x2  675x3  900x4  400x5

from a = 0 to b = 0.8. Use n = 4

Exact I = 1.64053

f (0)  0.2; f (0.2) 1.288;

5 points: 0 0.2 0.4 0.6 0.8
2 segments, h = 0.2 f (0.4)  2.456; f (0.6)  3.464;
f (0.8)  0.232 h = 0.2
I  I1  I 2


0 .2
0.2  4(1.288)  2.456  0.2 2.456  4(3.464)  0.232
3 3
 1.62347
Et  1.64053 1.62347  0.017067
 t  1.04% 0.85
Ea   4
(2400)  0.017067 24
180(4)
 CONT.
OR
(b  a)  
 f  x0   4  f  xi   2  f x j  f xn 
n1 n2
I
3n  i1,3,5... j2,4,6... 

f (0)  0.2; f (0.2) 1.288;

f (0.4)  2.456; f (0.6)  3.464;
f (0.8)  0.232

(0.8  0)
I 0.2  4(1.288  3.464)  2(2.456)  0.232
3(4)
 1.62347

Et  1.640533 1.623467  0.017067

 t  1.04%
0.85
Ea   4
(2400)  0.017067
180(4) 25
 SIMPSON’S 3/8 RULE
• Simpson’s 3/8 Rule results when third-order Lagrange
interpolating polynomial is used → uses 4 points
• 3/8 rule is slightly more accurate than 1/3 rule

(b  a )
Simpson' s 3/8 Rule : I   f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 )
3h
h
8 3
(b  a ) 5 ( 4 )
Approximat e error : Ea   f (x )
6480
# of segment = 3 (odd)
• Method is applicable/useful if # of points = 4 (even)
the number of segments is
multiplier of 3 e.g. 3, 6, 9…

26
 CLASS ACTIVITY
Use Simpson’s 3/8 rule to estimate the integral of

f (x)  0.2  25x  200x2  675x3  900x4  400x5

from a = 0 to b = 0.8. Exact I = 1.64053

(b  a )
I
3h
 f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 ) h
8 3
4 points: 0 4/15 8/15 4/5 (b  a ) 5 ( 4 )
Ea   f (x )
6480
 SOLUTION
Use Simpson’s 3/8 rule to estimate the integral of

f (x)  0.2  25x  200x2  675x3  900x4  400x5

from a = 0 to b = 0.8. Exact I = 1.64053

(b  a )
I
3h
 f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 ) h
8 3
4 points: 0 4/15 8/15 4/5 (b  a ) 5 ( 4 )
Ea   f (x )
6480
0.8  0
h  4/15 = 0.2667
3
f (0)  0.2
I  1.51917
f (0.2667) 1.432724
Et  |1.64053 – 1.51917|= 0.121363
f (0.5333)  3.487177
f (0.8)  0.232  t 7.4%
0.85
Ea   (2400)  0.121363
6480
 Combined Simpson’s Rules

Useful if you have

peculiar numbers
of segments.
Hence combine
1/3-rule and 3/8-
rule.

29
TRAPEZOIDAL RULE (MATLAB Code)
% Number of segments
n = 5;
% Integral limits
lowerLimit = 0;
upperLimit = 0.8;

% Discretization into equally spaced segments

x = linspace(lowerLimit,upperLimit,n+1);
% Define function to integrate
y = 0.2 + 2*x + 90*x.^2 - 120*x.^3 + 25*x.^4;

% Trapezoidal formula
A = sum(y(2:n));
M = 2*A;
I = (upperLimit-lowerLimit)*(y(1)+M+y(end))/(2*n);

disp(['Estimation of the integral is ', num2str(I)]);

Alternatively, use built-in function trapz

> x = linspace(0, 0.8, 6)
> f = x.^2 – 120*x + 4
> q = trapz(x,f)

30
 Concluding Remark
• We covers 3 types of Newton-Cotes Integration Schemes:

Trapezoidal Rule Simpson’s 1/3 Rule Simpson’s 5/8 Rule

• Simplest scheme • Needs 3 points • Needs 4 points

• Need 2 points • Even no. of • Slightly more
• not so accurate but segments accurate than 1/3
can be improved by • More accurate than rule but more
multiple-application trapezoidal rule computation
• Up to 1st-order • Up to third-order • Up to third-order
(linear) accurate accurate accurate
41

31
 Universiti Teknologi PETRONAS
MDB 3053
Numerical Methods

**Gauss-Legendre Integration Rule

Reading Assignment: Chapter 22

32
 **GAUSS-LEGENDRE RULE

The main idea of Gauss-Legendre:

n
f (x)dx   f (x )dx   wi f (x )i
b 1
a 1
i1

Gauss-Legendre Integration requires:

- Weighting functions: c0, c1
- Sampling points: x0, x1

Integral: I  c0 f ( x0 )  c1 f ( x1 )

33
 GAUSS-LEGENDRE RULE

The main idea of Gauss-Legendre:

I  c0 f ( x0 )  c1 f ( x1 )

• The integration strategy is to position any two points on a

curve to define a straight that would balance the positive and
negative area errors.

Integration by
Gauss-Legendre

34
 Method of Undetermined Coefficients
I  c0 f ( x0 )  c1 f ( x1 )
Consider area under constant line y=1
To integrate function y=1:
area, I = (ba)
I  c0 f (a)  c1 f (b)
(ba ) / 2

 c0 c1   1 dx  (b  a)
( b  a ) / 2
       (1)

Consider area under linear line y=x

To integrate: I  c0 f (a)  c1 f (b) area, I = 0
(ba ) / 2
 ba ba
I  c0   
 1c    x dx
 2   2  ( b  a ) / 2
 ba ba
 c0   
 1c 0
 2  2
35
solving : c0  c1          (2)
 Method of Undetermined Coefficients

From (1) and (2), we get c0 c1 b  a and c0  c1

ba
solving : c0  c1 
2
Rewrite integral :
ba ba
I f (a)  f (b) which is equivalent to
2 2 Trapezoidal rule!

• Next, we extend this Gauss-Legendre Rule to any area

I  c0 f ( x0 )  c1 f ( x1 )
need to find 4 unknowns!
c0, c1, x0, x1 are all unknowns
36
 Integration by Gauss-Legendre
Function:
1
Constant c0 f ( x0 )  c1 f ( x1 )   1 dx  2
1 To derive 2-point Gauss-Legendre formula
1
Linear c0 f ( x0 )  c1 f ( x1 )   x dx  0
1
Solved simultaneously to find
1 4 unknowns: c0 , c1, x0, x1
2
Quadratic c0 f ( x0 )  c1 f ( x1 )   x 2 dx 
1
3
1
Cubic c0 f ( x0 )  c1 f ( x1 )   x 3 dx  0
1

SOLUTION: c  c  1
Coefficients: 0 1

1
x0    0.5773503
3 I  c0 f ( x0 )  c1 f ( x1 )
1
x1   0.5773503
3
2-point Gauss-
 1   1  yields an integral estimate that is
Legendre I  f  f  third order accurate
formula  3  3 Chap22/37
 2-point Gauss-Legendre Rule
 1   1 
I  f   f   Error up to 3rd order (which is
 3  3 similar to Simpson’s 1/3 rule)

Integration limits
from 1 to +1

1 1
2 sampling x0   x0  
points: 3 3
38
simplest
GAUSS-LEGENDRE RULE
Gauss-Legendre Rule from 2-points to 5-point
Rules

sampling pts

39
EXAMPLE

assume: x = a + bxd
0 = a + b(-1)
0.8 = a + b(1)
solve: a = 0.4, b = 0.4
rewrite: x = 0.4 + 0.4xd

Chap21/36

40
 *CLASS ACTIVITY

• Use 2-point Gauss Legendre formula to find the integral of:

STEP 1: convert the integration limit to [1,+1]

i.e. coordinates transform from [a,b]=[0,3] to [ 1,1]
Assume: x = at + b find a and b?

Rewrite: x = 3t/2 + 3/2 derivative: dx = 1.5 dt

41
 SOLUTION

STEP 2: Substitute into integral

STEP 3: use 2-point Gauss-Legendre formula

I  c0 f (x0 )  c1 f (x1 )
c0  c1 1
1 1
x0   , x1 
3 3
 1  1   1  1
I  f  f   f    1.9174, f    2.7520
 3  3  3  3
Hence: I = 1.9174 + 2.7520 = 4.6694
error t = 0.06%

42
Universiti Teknologi PETRONAS
MDB3053
Numerical Methods

DIFFERENTIATION

Chapter 23 in Textbook

39

43
 NUMERICAL DIFFERENTIATION

• Taylor Series (infinite terms):

f ( xi ) f n ( xi )
forward f ( xi 1 )  f ( xi )  f ( xi )( xi 1  xi )  xi 1  xi    
2
( xi 1  xi ) n  Rn
point 2! n!

• Consider first-order approximation (chop after two terms)

1st order derivative
f ( xi 1 )  f ( xi )  f ( xi )( xi 1  xi )  

• Re-arrange to solve for first-order derivative. This

become the Numerical Differentiation scheme:

f ( xi 1 )  f ( xi )
f ( xi )  Error is first order O(h).
( xi 1  xi )

44
 NUMERICAL DIFFERENTIATION

• 3 types of finite-divided difference schemes:

Forward, Backward and Centered

Forward Backward Centered

forward pt backward pt
forward pt base pt base pt backward pt

f ( xi 1 )  f ( xi ) f ( xi )  f ( xi 1 ) f ( x )  f ( xi 1 )  f ( xi 1 )
f ' ( xi )  f ( xi )  i
2h
h h
Truncation error = O(h) Truncation error = O(h) Truncation error = O(h2)

• Up to 1st order • Up to 1st order • Up to 2nd order

(linear) accurate (linear) accurate accurate

45
 NUMERICAL DIFFERENTIATION
• If we retain second-derivative term in Taylor series, the new
differential scheme becomes more accurate.

• Example: Forward scheme by chopping after 2nd order term

f ( xi ) 2
f ( xi 1 )  f ( xi )  f ( xi )h  h 
2
f ( xi 1 )  f ( xi ) f ( xi )
f ( xi )   h  O(h 2 )
h 2
2nd derivative  f ( xi )  f ( xi  2 )  2 f ( xi 1 )  f ( xi )  O(h)
2
error O(h) h
f ( xi 1 )  f ( xi ) f ( xi  2 )  2 f ( xi 1 )  f ( xi )
f ( xi )   2
h  O ( h 2
)
h 2h
1st derivative  f ( x )   f ( xi  2 )  4 f ( xi 1 )  3 f ( xi )  O(h 2 )
i
(improved accuracy) 2h
error O(h2)

46
 NUMERICAL DIFFERENTIATION

• If we retain second-derivative term in Taylor series, the

new differential scheme becomes more accurate:
Forward Backward Centered

 f ( xi  2 )  4 f ( xi 1 )  3 f ( xi ) 3 f ( xi )  4 f ( xi 1 )  f ( xi 2 ) f ( xi 1 )  f ( xi 1 )
f ' ( xi )  f ' ( xi )  f ( xi ) 
2h 2h 2h

Truncation error = O(h2) Truncation error = O(h2) Truncation error = O(h2)

• 2nd order accurate • 2nd order accurate • 2nd order accurate

Centered

 f ( xi 2 )  8 f ( xi 1 )  8 f ( xi 1 )  f ( xi 2 )
f ' ( xi ) 
12h

Truncation error = O(h6)

• 6th order accurate
47
Forward
From first-
to forth-
order
derivative

48
Backward
From first- to forth-order
derivative

Chap22/45

49
Centered
From first- to forth-
order derivative

50
 CLASS ACTIVITY
• Table shows the distance traveled versus time for a rocket
t (s) 0 25 50 75 100
y (km) 0 32 58 78 92
Estimate the velocity (km/s) and acceleration (km/s2) for each point by using
finite-difference approximations with second-order error O(h2)

51
 SOLUTION
t (s) 0 25 50 75 100
y (km) 0 32 58 78 92
Calculate v (km/s) and a (km/s2) for each point.

at t = 0. Use forward difference scheme:

velocity at t=0:
Acceleration at t=0,

at t = 25s Use centered difference scheme:

velocity at t=25 s
Acceleration at t=25 s
Table of solution

Continue the centered diff. scheme for

t = 50, 75 s
at t = 100 s, use backward diff. scheme
52
 SUMMARY
• Different types of Numerical Integration

• Newton-Cotes: Trapezoidal Rule, Simpson’s 1/3 Rule,

Simpson’s 3/8 Rule
• Gauss-Legendre Rule

• Numerical Differentiation
• Forward, Backward and Centred
• First order and second order

 We have covered Chapter 21, 22, 23 in textbook

53