FALL 2020

CE 305
NUMERICAL METHODS FOR ENGINEERS

MIDDLE EAST TECHNICAL UNIVERSITY

CIVIL ENGINEERING DEPARTMENT
These lecture notes have been prepared by multiple lecturers between 2007 and 2020 for the undergraduate
course CE 305 “Numerical Methods for Engineers” in METU Department of Civil Engineering. We acknowledge all
of these efforts of contributors Yalın Arıcı, Aysegul Askan, Alp Caner, Serdar Göktepe, Engin Karaesmen,
Shaghayegh Karimzadeh, Elçin Kentel, Mete Koken, Onur Pekcan, Kartal Toker and Ahmet Cevdet Yalçıner. The
initial typing of the original handwritings was performed by Eren Yağmur. The final typesetting was prepared by
Shaghayegh Karimzadeh. All rights are reserved by METU Department of Civil Engineering. Copying or distribution
without permission is forbidden.
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

CHAPTER 6: NUMERICAL INTEGRATION

Motivation: We need numerical integration in case analytical integration is not possible or complex

• Some functions cannot be integrated

• Some integrals are complex
• The function is unknown (i.e. it is in the form of data points)

6.1. TRAPEZOIDAL RULE AND ITS COMPOSITE METHOD

6.1.1. Trapezoidal Rule

f(x)

Real f(x) (may or may

not be known)

Approx. f(x)

Error

𝑏
∫𝑎 𝑓(𝑥)𝑑𝑥
x
a b

𝑏 𝑓(𝑏)+𝑓(𝑎)
∫𝑎 𝑓(𝑥)𝑑𝑥 ≈ 𝐴𝑟𝑒𝑎 𝑜𝑓 𝑡𝑟𝑎𝑝𝑒𝑧𝑜𝑖𝑑𝑎𝑙 = (𝑏 − 𝑎) 2

(𝑏−𝑎)3
𝜀𝐴𝑇 = − 𝑓 ′′ (𝜉) where 𝑎 < 𝜉 < 𝑏
12

This is not good if a and b are far apart.

If the function is linear the error is zero.

6-1
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

6.1.2. Composite Trapezoidal Rule

f(x)

𝐴1 𝐴2 … 𝐴𝑛

x
𝑥0 𝑥1 𝑥2 … 𝑥𝑛−1 𝑥𝑛

𝑥𝑛
∫𝑥 𝑓(𝑥)𝑑𝑥 = 𝑠𝑢𝑚 𝑜𝑓 𝑛 𝑡𝑟𝑎𝑝𝑒𝑧𝑜𝑖𝑑 𝑠𝑙𝑖𝑐𝑒𝑠
0

𝑓(𝑥𝑖 )+𝑓(𝑥𝑖−1 )
= ∑𝑛𝑖=1 𝐴𝑖 = ∑𝑛𝑖=1 [(𝑥𝑖 − 𝑥𝑖−1 ) ]
2

If intervals of x are constant, then:

𝑥𝑛 −𝑥0
∀ 𝑖, 𝑥𝑖 − 𝑥𝑖−1 = ℎ =
𝑛

𝑥 ℎ
∫𝑥 𝑓(𝑥)𝑑𝑥 = 2 [𝑓(𝑥0 ) + 2 ∑𝑛−1
𝑛
0
𝑖=1 𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑛 )]

(𝑥𝑖 −𝑥𝑖−1 )3
Error of 𝑖𝑡ℎ trapezoid=− 𝑓 " (𝜉) where 𝑥𝑖−1 < 𝜉 < 𝑥𝑖
12

(𝑥𝑖 −𝑥𝑖−1 )3
Error of composite trapezoid=∑𝑛𝑖=1 [− 𝑓 " (𝜉𝑖 )] where 𝑥𝑖−1 < 𝜉𝑖 < 𝑥𝑖
12

𝑥𝑛 −𝑥0 ℎ3 𝑛 ℎ3 𝑛 "
For constant intervals of ℎ = 𝜀𝐴𝑇 = − ∑ 𝑓 " (𝜉𝑖 ) ≈− 𝑓 𝑎𝑣
𝑛 12 𝑖=1 12

(𝑥𝑛 −𝑥0 )3
or alternatively 𝜀𝐴𝑇 = − 𝑓 " 𝑎𝑣
12𝑛2

∑𝑛 "
𝑖=1 𝑓 (𝜉𝑖 )
where 𝑓 " 𝑎𝑣 = that is the average of 𝑓 " over [𝑥0 , 𝑥𝑛 ]
𝑛

6-2
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

6.2. SIMPSON'S RULE AND ITS COMPOSITE METHOD

6.2.1. Simpson’s 1/3 Rule

𝑓(𝑥)

Real f(x) (may or may not be known)

Approx. f(x): 𝑓̃2 (𝑥)

2nd order parabola

𝑥2
∫ 𝑓(𝑥)𝑑𝑥
𝑥0
𝑥
𝑥0 𝑥1 𝑥2

h h

𝑥0 +𝑥2
From 3 points with equal intervals in x (i.e. 𝑥1 = )
2

f(x) is approximated by a parabola

Is based on second-order polynomial (parabola) interpolation using 3 points with equal intervals in x
𝑥 +𝑥 𝑥 −𝑥
(i.e. 𝑥1 = 0 2 𝑜𝑟 ℎ = 2 1):
2 2
𝑥2 𝑥2
𝐼 = ∫ 𝑓 (𝑥)𝑑𝑥 = ∫ 𝑓2 (𝑥)𝑑𝑥
𝑥0 𝑥0
𝑥2
(𝑥 − 𝑥1 )(𝑥 − 𝑥2 ) (𝑥 − 𝑥0 )(𝑥 − 𝑥2 )
=∫ [ 𝑓(𝑥0 ) + 𝑓(𝑥1 )
𝑥0 (𝑥0 − 𝑥1 )(𝑥0 − 𝑥2 ) (𝑥1 − 𝑥0 )(𝑥1 − 𝑥2 )
(𝑥 − 𝑥0 )(𝑥 − 𝑥1 ) 1
+ 𝑓(𝑥2 )] 𝑑𝑥 ≅ ℎ[𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )]
(𝑥2 − 𝑥0 )(𝑥2 − 𝑥1 ) 3

𝑥 1
Then, integral of this parabola is, exactly: 𝐼 = ∫𝑥 2 𝑓̃2 (𝑥) = ℎ[𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )]
0 3

In other words:
𝑥2
1
𝐼 = ∫ 𝑓 (𝑥)𝑑𝑥 ≅ ℎ[𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )]
𝑥0 3

(𝑥2 −𝑥0 )5 ℎ5 (4) 𝑑4 𝑓

𝜀𝐴𝑇 = − 𝑓 (4) (𝜉) = − 𝑓 (𝜉) where 𝑥0 < 𝜉 < 𝑥2 and 𝑓 (4) (𝜉) =
2880 90 𝑑𝑥 4

6-3
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

Error of Simpson’s 1/3 rule is proportional to 𝑓 (4). This means Simpson can get the integral of a
third-degree polynomial (cubic polynomial) exactly, as 𝑓 (4) will be 0.

6.2.2. Composite Simpson’s 1/3 Rule

*Intervals must be equal

*If the integral involves n+1 data points from x0 to xn: The number n must be even (odd # of data
points)
𝑥 𝑛
*∫𝑥 𝑛 𝑓 (𝑥)𝑑𝑥 = 𝑠𝑢𝑚 𝑜𝑓 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙𝑠 𝑐𝑜𝑟𝑟𝑒𝑠𝑝𝑜𝑛𝑑𝑖𝑛𝑔 𝑡𝑜 𝑝𝑎𝑟𝑎𝑏𝑜𝑙𝑎𝑠
0 2
f(x)

n+1= 7 data points, n=6, n/2=3 parabolas

𝒉 = (𝒙𝟔 − 𝒙𝟎 )/𝟔
parabola 2

parabola 3
parabola 1

𝐴1 𝐴2 𝐴3

𝑥0 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6

h h

Note: In this Chapter

• “h” is always the distance between 2 successive data points.

• Each slice is composed of 2h distances from xi-1 to xi+1.
• 1 Simpson parabola passes through each slice.

For example, for the above figure, applying Simpson’s 1/3 Rule for each of the three slices results in
the following:
𝑥
6 ℎ
∫𝑥 𝑓(𝑥)𝑑𝑥 = 3 [𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 2𝑓(𝑥2 ) + 4𝑓(𝑥3 ) + 2𝑓(𝑥4 ) + 4𝑓(𝑥5 ) + 𝑓(𝑥6 )]
0

ℎ
𝐼 1 = [𝑓(𝑥0 ) + 4 ∑ 𝑓 (𝑥 𝑜𝑑𝑑 ) + 2 ∑ 𝑓 (𝑥 𝑒𝑣𝑒𝑛 ) + 𝑓(𝑥𝑛 )]
(𝑐𝑜𝑚𝑝𝑜𝑠𝑖𝑡𝑒 3 𝑆𝑖𝑚𝑝𝑠𝑜𝑛′ 𝑠 𝑟𝑢𝑙𝑒) 3 1:𝑛−1 2:𝑛−2

General form for n+1 data points:

6-4
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

𝑛
𝑥𝑛 ℎ 𝑛/2 −1
2
∫𝑥 𝑓(𝑥)𝑑𝑥 = 3 [𝑓(𝑥0 ) + 4 ∑𝑖=1 𝑓(𝑥2𝑖−1 ) + 2 ∑𝑖=1 𝑓(𝑥2𝑖 ) + 𝑓(𝑥𝑛 )]
0

**Error of 𝑖𝑡ℎ parabola, which includes point at 𝑥2𝑖−2 , 𝑥2𝑖−1 , 𝑥2𝑖 is:

ℎ5 (4)
𝜀𝐴𝑇 = − 𝑓 (𝜉2𝑖 ) where at 𝑥2𝑖−2 < 𝜉2𝑖 < 𝑥2𝑖
90

**Error of composite Simpson is:

𝑛
5 2
ℎ ℎ5 𝑛 (4) (𝑥𝑛 − 𝑥0 )5 (4)
𝜀𝐴𝑇 = − ∑ 𝑓 (4) (𝜉2𝑖 ) ≈ − 𝑓 𝑎𝑣 = 𝑓 𝑎𝑣
90 180 180𝑛4
𝑖=1

where 𝑓 (4)𝑎𝑣 is the average 𝑓 (4) over [𝑥0 , 𝑥𝑛 ].

Important note: In both composite methods, the theoretical error decreases with increasing number
of intervals/considering more slices (or smaller h), but with 𝑛 ↑ 𝑎𝑛𝑑 ℎ ↓, rounding errors accumulate
more.

Practical note: If you have even number of data points but want to use Composite Simpson, you can
use it up to the (𝑛)𝑠𝑡 point and add a trapezoid slice at the beginning or the end.

𝑥 𝑓(𝑥)
Ex.: −1 −0.5
0 1
1 1
2 0
3 −1
4 0
5 7

5
a) Estimate ∫−1 𝑓(𝑥)𝑑𝑥 using Comp. Trapezoid
5
b) Estimate ∫−1 𝑓(𝑥)𝑑𝑥 using Comp. Simpson
4
c) Estimate ∫−1 𝑓(𝑥)𝑑𝑥 using Comp. Simpson
5
d) The real function is 𝑓(𝑥) = 2𝑥 − 𝑥 2 . Calculate ∫−1 𝑓(𝑥)𝑑𝑥 analytically and compare it with
the results of parts a and b.

Solution
𝑥 ℎ 1
a) ∫𝑥 6 𝑓(𝑥)𝑑𝑥 = [𝑓(𝑥0 ) + 2 ∑6−1
𝑖=1 𝑓(𝑥𝑖 ) + 𝑓(𝑥6 )] = 2 [−0.5 + 2(1 + 1 + 0 − 1 + 0) + 7] =
0 2
4.25
𝑥 ℎ 6/2 6/2−1 1
b) ∫𝑥 6 𝑓(𝑥)𝑑𝑥 = [𝑓(𝑥0 ) + 4 ∑𝑖=1 𝑓(𝑥2𝑖−1 ) + 2 ∑𝑖=1 𝑓(𝑥2𝑖 ) + 𝑓(𝑥6 )] = [−0.5 +
0 3 3
4(1 + 0 + 0) + 2(1 − 1) + 7] = 3.5

6-5
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

c) 𝑆𝑖𝑚𝑝𝑠𝑜𝑛|4−1 → 𝑆𝑖𝑚𝑝𝑠𝑜𝑛|3−1 + 𝑇𝑟𝑎𝑝𝑒𝑧𝑜𝑖𝑑 |43

1 1
[−0.5 + 4(1 + 0) + 2 × 1 − 1] + (−1 + 0) = 1
3 2

5
5 2𝑥 𝑥3 25 53 2−1 (−1)3
d) ∫−1(2𝑥 − 𝑥 2 )𝑑𝑥 = ( − )| =( − )−( − ) = 3.4449
𝑙𝑛2 3 −1 𝑙𝑛2 3 𝑙𝑛2 3

Simpson is much closer to the true answer. Why? See the function 𝑓(𝑥) = 2𝑥 − 𝑥 2 below

1.5

0.5

0
-2 -1 0 1 2 3 4 5
-0.5

-1

-1.5

6.3. GAUSS QUADRATURE

Modern numerical integration approaches mostly rely on Gauss quadrature method. The particular
formulas described herein are Gauss Legendre formulas. The analytical form of the function to be
integrated must be known.

y
f(x)
1
𝑤𝑖

x
-1 -1/√3 1/√3 1

𝑥𝑖

1
∫−1 𝑓(𝑥)𝑑𝑥 ≈ ∑𝑖 𝑤𝑖 𝑓(𝑥𝑖 ) (See table in the handout for higher order quadrature)

6-6
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

𝑤𝑖 f(x)
5/9

8/9

5/9

x
-1 0 1
3 3
−√ √
5 5
𝑥𝑖

22𝑛+3 [(𝑛+1)!]4
Error: (2𝑛+3)[(2𝑛+2)!]3
𝑓 (2𝑛+2) (𝜉) −1<𝜉 <1

where n: is the number of Gauss-Quadrature points minus 1

Thus:

2 point-Gauss-quadrature: Exact up to cubic polynomial (where error includes 4th order derivative)

3 point-Gauss-quadrature: exact up to 5th-order polynomial (where error includes 6th order

derivative)

*Finally, the only issue is the boundaries of the integral. Transform the function so that it fits the
quadrature:
𝑏
∫𝑎 𝑔(𝑡)𝑑𝑡 is the original question (Before starting your GQ solution, replace the independent
variable such as x,z,y with t’s in the given function):
𝑏 𝑏
Ex1: ∫𝑎 𝑔(𝑥)𝑑𝑥 replace x’s wıth t’s: ∫𝑎 𝑔(𝑡)𝑑𝑡

𝑏 𝑏
Ex2: ∫𝑎 𝑔(𝑧)𝑑𝑧 replace z’s with t’s: ∫𝑎 𝑔(𝑡)𝑑𝑡

We need to transform 𝑡 ∈ [𝑎, 𝑏] into 𝑥 ∈ [−1,1]

𝑏+𝑎 𝑏−𝑎
t=c0+c1x and when t=a x=-1 and when t=b x=1 will yield 𝑐0 = and 𝑐1 =
2 2

6-7
 METU - Department of Civil Engineering

Course Code and Title: CE 305 Numerical Methods for Engineers (3-1) 3

𝑏+𝑎 𝑏−𝑎 𝑏−𝑎

This can be done by: 𝑡= + 𝑥 𝑎𝑛𝑑 𝑑𝑡 = 𝑑𝑥
2 2 2

𝑏 1 𝑏+𝑎 𝑏−𝑎 𝑏−𝑎 1

∫𝑎 𝑔(𝑡)𝑑𝑡 = ∫−1 𝑔 ( + 𝑥) 𝑑𝑥 = ∫−1 𝑓(𝑥)𝑑𝑥 ≈ ∑𝑖 𝑤𝑖 𝑓(𝑥𝑖 )
2 2 2

𝑓(𝑥)
5
Ex: ∫−1(2𝑥 − 𝑥 2 )𝑑𝑥 =? Use 3-point Gauss quadrature rule to integrate

Solution:

5 5−1 5−(−1)
∫−1(2𝑡 − 𝑡 2 )𝑑𝑡 𝑡=
2
+
2
𝑥 = 2 + 3𝑥 𝑑𝑡 = 3𝑑𝑥

1
∫−1[2(2+3𝑥) − (2 + 3𝑥)2 ]. 3𝑑𝑥 =?

𝑓(𝑥) = 3[2(2+3𝑥) − (2 + 3𝑥)2 ]

2
5 (2−3√3⁄5) 8
𝐼 = . 3. [2 − (2 − 3√3⁄5) ] + . 3. [2(2+3(0)) − (2 + 3(0))2 ]
9 9
2
5 (2+3√3⁄5)
+ . 3. [2 − (2 + 3√3⁄5) ] = 3.37466
9

6-8