Section 2: Numerical Differentiation

1. Introduction and motivation

2. Review of differentiation
3. Finite difference method
4. Differentiation with Built-in ‘Symbolic Math’
in Matlab
5. Theoretical Examples
6. Matlab Examples
Textbook sections: 2.2, 2.6, 2.7, 8.1 to 8.4, 8.7, 8.9

ME 3060 Section 2: Numerical differentiation 1

1. Introduction and motivation
• When performing experiments the measured data is typically
discrete -> only a finite number of data points are available
• Derivatives are often needed for application of physical laws or for
data analysis
A simple problem from the Fluids Lab (ME-3160): Heat transfer
through a fin
• Temperature is measured
Thermocouple W ells
Heater
at N discrete points
𝑥𝑖 , i = 1,2, … , N with a
thermo-couple
• We need to find the heat-
x dx
transfer rate at several
points to confirm a
x = 0
proposed physical law
Cross-sectional Area, A

ME 3060 Section 2: Numerical differentiation 2

 1. Introduction and motivation
From the measurements we have

Calculate heat transfer q [W]

𝑑𝑇
𝑞 = −k ⋅ 𝐴 ⋅
𝑑𝑥

𝑑𝑇
How do we get if we don’t know
𝑑𝑥
T(x) and only have data at some
points?

• Solution 1: Fit an analytical curve through data points and then

calculate the derivative (good for noisy data, done in section 7)

𝑑𝑇 𝑇 𝑥𝑖+1 −𝑇(𝑥𝑖 )
• Finite-difference approximation of the derivative: ≈
𝑑𝑥 𝑥𝑖+1 −𝑥𝑖
ME 3060 Section 2: Numerical differentiation 3
1. Introduction and motivation
A second important application: numerical solution of ODEs!!!
𝑑𝑦
= −1.2 𝑦 + 7𝑒 −0.3𝑥 = 𝑓 𝑦, 𝑥 , y x1 = 0 = y1 = 3
𝑑𝑥
Find a numerical solution in the interval 𝑥 ∈ 𝑎 = 0, 𝑏 = 2.5 at N-discrete points:
𝑥1 = 0, 𝑥2 = 𝑥1 + ℎ, 𝑥3 = 𝑥1 + 2ℎ, … 𝑥𝑖 − 𝑥1 + 𝑖 − 1 ℎ, …
𝑏−𝑎
𝑥𝑁 = 2.5 = 𝑥1 + 𝑁 − 1 ℎ ---> h=
𝑁−1
Approximate derivative with a “finite difference”
𝑑𝑦 𝑦 𝑥𝑖+1 − 𝑦 𝑥𝑖
ቚ ≈ = −1.2 𝑦(𝑥𝑖 ) + 7𝑒 −0.3𝑥𝑖
𝑑𝑥 𝑥𝑖 ℎ
−0.3𝑥𝑖
Explicit numerical solution formula: 𝑦 𝑥𝑖+1 = 𝑦 𝑥𝑖 − ℎ 1.2 𝑦(𝑥𝑖 ) + 7𝑒

N=51
h = 0.05

ME 3060 Section 2: Numerical differentiation 4

2. Review of differentiation
Limit of a function: lim 𝑓(𝑥) = 𝑓 𝑎 = 𝐿
𝑥→𝑎

Definition: If f(x) is a function defined on an open interval containing a, and L

is a real number, then for each number 𝜀 > 0 there exists a number 𝛿 > 0
such that if 0 < 𝑥 − 𝑎 < 𝛿 then 𝑓 𝑥 − 𝐿 < 𝜀.

Continuity of a function: A function f(x) is said to be continuous at

x = a if the following three conditions are satisfied
i. f(a) exists A function is said to be continuous on an
ii. lim 𝑓(𝑥) exist open interval (a,b) if it is continuous at
𝑥→𝑎
each point in the interval.
iii. lim 𝑓 𝑥 = 𝑓(𝑎)
𝑥→𝑎

ME 3060 Section 2: Numerical differentiation 5

2. Review of differentiation
Intermediate value theorem:
If f(x) is continuous on [a,b] (closed interval)
and M is any number between f(a) and f(b)
then there exists at least one number 𝑐 ∈
[𝑎, 𝑏] such that f(c) = M.

The derivative: describes the rate of change of a function or,

equivalently, its slope. Formally, it is defined as a limit:

𝑑𝑓(𝑥) ′
𝑓 𝑥 − 𝑓(𝑎)
ቚ = 𝑓 𝑎 = lim
𝑑𝑥 𝑥=𝑎 𝑥→𝑎 𝑥−𝑎
Or, more general
𝑓 𝑥 + ℎ − 𝑓(𝑥)
𝑓′ 𝑥 = lim
ℎ→0 ℎ
ME 3060 Section 2: Numerical differentiation 6
 2. Review of differentiation
𝑓 𝑥+ℎ −𝑓(𝑥)
If the limit lim exists, then f(x) is said to be differentiable at x
ℎ→0 ℎ

f(x) differentiable f(x) continuous

Higher order derivatives: are obtained by taking

successive derivatives of derivatives:
𝑑𝑓 𝑑𝑓 𝑑𝑓
𝑑 2
𝑑 𝑓 ȁ𝑥+ℎ − ȁ𝑥
𝑑𝑥 𝑑𝑥 𝑑𝑥
= 2 = lim
𝑑𝑥 𝑑𝑥 ℎ→0 ℎ
A function f(x) is called a smooth function if all of it’s derivatives
are continuous!

ME 3060 Section 2: Numerical differentiation 7

2. Review of differentiation
Differentiation rules:
Product rule: for f(x) and g(x) differentiable,
𝑑
𝑓 𝑥 ⋅ 𝑔(𝑥) = 𝑓 ′ 𝑔 + 𝑔′f
𝑑𝑥
Quotient rule: for f(x) and g(x) differentiable, 𝑔(𝑥) ≠ 0

𝑑 𝑓 𝑥 𝑓′𝑔 − 𝑔′f
=
𝑑𝑥 𝑔 𝑥 𝑔2
Chain rule: for f(x) and g(x) differentiable
𝑑
𝑑 𝑑𝑓 𝑑𝑔 Example: exp(sin 𝑥 2 ) =
𝑓 𝑔(𝑥) = 𝑑𝑥
𝑑𝑥 𝑑𝑔 𝑑𝑥
𝑑 exp(sin 𝑥 2 ) 𝑑(sin 𝑥 2 )
2
×
𝑑(sin 𝑥 ) 𝑑𝑥

= exp sin 𝑥 2 cos 𝑥 2 2𝑥

ME 3060 Section 2: Numerical differentiation 8
2. Review of differentiation
Mean value theorem
If f(x) is continuous on the closed interval [a,b] and differentiable on the
open interval (a,b), then there exists a number 𝑐 ∈ 𝑎, 𝑏 such that

𝑓 𝑏 − 𝑓(𝑎) 𝑑𝑓
= ቚ = 𝑓′ 𝑐
𝑏−𝑎 𝑑𝑥 𝑥=𝑐
Important theorem for numerical analysis because
it allows us to determine bounds of the numerical
error for our chosen differentiation
approximation!

Partial derivatives: For functions depending on more than one variable

we can define partial derivatives
𝜕𝑓(𝑥, 𝑦) 𝑓 𝑥 + ℎ, 𝑦 − 𝑓(𝑥, 𝑦)
= 𝑓𝑥 = lim
𝜕𝑥 ℎ→0 ℎ
ME 3060 Section 2: Numerical differentiation 9
2. Review of differentiation
Higher order partial derivatives can be defined similarly to the derivative
𝜕 2 𝑓 𝑥, 𝑦 𝑓𝑥 𝑥 + ℎ, 𝑦 − 𝑓𝑥 (𝑥, 𝑦)
2
= 𝑓𝑥𝑥 = lim
𝜕𝑥 ℎ→0 ℎ
𝜕2𝑓 𝑓𝑥 𝑥, 𝑦 + ℎ − 𝑓𝑥 (𝑥, 𝑦)
= 𝑓𝑥𝑦 = lim
𝜕𝑥 𝜕𝑦 ℎ→0 ℎ
Total differential
Example: For a function of three variables f(x,y,z) with continuous partial
derivatives the total differential is
𝜕𝑓 𝜕𝑓 𝜕𝑓
𝑑𝑓 = 𝑑𝑥 + 𝑑𝑦 + 𝑑𝑧
𝜕𝑥 𝜕𝑦 𝜕𝑧
The total differential allows us to find e.g. the time rate of change of a
function f that depends on three variables x(t), y(t), z(t) that all depend
on time:
𝑑𝑓 𝜕𝑓 𝑑𝑥 𝜕𝑓 𝑑𝑦 𝜕𝑓 𝑑𝑧
= + +
𝑑𝑡 𝜕𝑥 𝑑𝑡 𝜕𝑦 𝑑𝑡 𝜕𝑧 𝑑𝑡
ME 3060 Section 2: Numerical differentiation 10
2. Review of differentiation
Taylor series expansion
Assume that we are given the value of a function at a point 𝑥0 i.e 𝑓(𝑥0 ) and all
the values of all the derivatives at that point i.e 𝑓 ′ 𝑥0 , 𝑓 ′′ 𝑥0 , 𝑓 ′′′ 𝑥0 ,… Can
we extrapolate how the function f(x) looks in a neighborhood around 𝑥0 ?

Yes, using the Taylor series expansion

𝑥 − 𝑥0 2
𝑓 𝑥 = 𝑓 𝑥0 + 𝑓 ′ 𝑥0 𝑥 − 𝑥0 + 𝑓 ′′ 𝑥0 +
2!
𝑥 − 𝑥0 3 𝑥 − 𝑥0 4
𝑓 ′′′ 𝑥0 + 𝑓 (4) 𝑥0 +⋯
3! 4!
∞
𝑛
(𝑛)
𝑥 − 𝑥0
𝑓 𝑥 = ෍ 𝑓 (𝑥0 )
𝑛!
𝑛=0

ME 3060 Section 2: Numerical differentiation 11

2. Review of differentiation
Often we will want to express the Taylor series at an arbitrary point x in a
neighborhood 𝑥 + Δ𝑥 i.e we want to have 𝑓 𝑥 + Δ𝑥 = ⋯ This can simply be
achieved by using the previous definition
2
𝑥 − 𝑥 0
𝑓 𝑥 = 𝑓 𝑥0 + 𝑓 ′ 𝑥0 𝑥 − 𝑥0 + 𝑓 ′′ 𝑥0 +⋯
2!
and a change of variable 𝑥0 → 𝑥, 𝑥 → x + Δ𝑥 and thus
𝑥 − 𝑥0 → 𝑥 + Δ𝑥 − 𝑥 = Δ𝑥
∞ Δ𝑥 𝑛
𝑓 𝑥 + Δ𝑥 = ෍ 𝑓 (𝑛) 𝑥
𝑛=0 𝑛!

To use the Taylor series as a tool to approximate a function we obviously have to

truncate the series at some point, e.g. use only the first three terms:
Δ𝑥 2
𝑓 𝑥 + Δ𝑥 ≈ 𝑓 𝑥 + 𝑓 ′ 𝑥 Δ𝑥 + 𝑓 ′′ 𝑥
2!
Can we estimate the error that we commit by truncating the
infinite series?

ME 3060 Section 2: Numerical differentiation 12

 2. Review of differentiation
Taylor’s theorem:
If f and its first n derivatives f’, f’’, …, f(n) are continuous on [a,b] and f(n) is
differentiable on (a,b), then there exists a number 𝜉 between x and 𝑥 + Δ𝑥 such
that
Δ𝑥 2 Δ𝑥 𝑛 Δ𝑥 𝑛+1
𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 ′ 𝑥 Δ𝑥 + 𝑓 ′′ 𝑥 + ⋯+ 𝑓 𝑛 𝑥 + 𝑓 𝑛+1 (𝜉)
2! 𝑛! (𝑛 + 1)!

Remainder 𝑅𝑛
The remainder is 𝑅𝑛 ~Δ𝑥 𝑛+1 and thus the truncation error is:

𝐸 𝑇𝑅 = 𝑓𝑡𝑟𝑢𝑒 − 𝑓𝑎𝑝𝑝𝑟𝑜𝑥 ~ Δ𝑥 𝑛+1

We say the truncation error is of order Δ𝑥 𝑛+1 or 𝒪(Δ𝑥 𝑛+1 )

ME 3060 Section 2: Numerical differentiation 13

3. Finite difference method
Fin heat transfer example: Find heat transfer rate 𝑞(𝑥𝑖 ) given discrete
measurements 𝑇 𝑥𝑖 at the points 𝑥1 , 𝑥2 , … , 𝑥𝑁
𝑑𝑇
𝑞 = −k ⋅ 𝐴 ⋅
𝑑𝑥
Recall definition of derivative (infinitesimal
slope ≈
𝑑𝑇
ȁ difference)
𝑑𝑥 𝑥4
𝑑𝑇 𝑇 𝑥 + ℎ − 𝑇(𝑥)
= lim
𝑑𝑥 ℎ→0 ℎ
Finite difference
𝑥4 𝑥5
𝑑𝑇 𝑇 𝑥 + ℎ − 𝑇(𝑥)
≈
In the example we have 𝑑𝑥 ℎ
𝑇 𝑥 = 𝑇(𝑥4 )
𝑑𝑇 𝑇 𝑥5 − 𝑇(𝑥4 )
𝑇 𝑥 + ℎ = 𝑇(𝑥5 ) ቚ ≈
𝑑𝑥 𝑥4 𝑥5 − 𝑥4
ℎ = 𝑥5 − 𝑥4

ME 3060 Section 2: Numerical differentiation 14

3. Finite difference method
Fin heat transfer example: Find heat transfer rate 𝑞(𝑥𝑖 ) given discrete
measurements 𝑇 𝑥𝑖 at the points 𝑥1 , 𝑥2 , … , 𝑥𝑁
Using finite difference:
𝑇 𝑥𝑖+1 − 𝑇(𝑥𝑖 )
𝑞 𝑥𝑖 ≈ −k ⋅ 𝐴 ⋅
𝑥𝑖+1 − 𝑥𝑖
Works for i = 1,2,…,N-1 not for N!!

Small amount of noise in the T –

measurements cause “non-smooth”
behavior of the finite difference
𝑇 𝑥𝑖+1 −𝑇(𝑥𝑖 )
approximation 𝑥𝑖+1 −𝑥𝑖
Noise get’s amplified when
differentiating!

ME 3060 Section 2: Numerical differentiation 15

3. Finite difference method
Two-point finite difference methods

Forward difference (general definition)

𝑑𝑓 𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖 𝑓𝑖+1 − 𝑓𝑖
ቚ ≈ =
𝑑𝑥 𝑥=𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖

Backward difference
𝑑𝑓 𝑓 𝑥𝑖 − 𝑓(𝑥𝑖−1 )
ቚ ≈
𝑑𝑥 𝑥=𝑥𝑖 𝑥𝑖 − 𝑥𝑖−1

Central difference
𝑑𝑓 𝑓 𝑥𝑖+1 − 𝑓(𝑥𝑖−1 )
ቚ ≈
𝑑𝑥 𝑥=𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖−1

ME 3060 Section 2: Numerical differentiation 16

3. Finite difference method
Fin heat transfer example:
Find heat transfer rate 𝑞(𝑥𝑖 ) given discrete
measurements 𝑇 𝑥𝑖 at the points 𝑥1 , 𝑥2 , … , 𝑥𝑁
“Optimal solution” using two-point finite difference:

First point i = 1 ---> forward difference

𝑇 𝑥2 − 𝑇(𝑥1 )
𝑞(𝑥1 ) = −k ⋅ 𝐴 ⋅
𝑥2 − 𝑥1
Points i = 2,N-1 ---> central difference
𝑇 𝑥𝑖+1 − 𝑇(𝑥𝑖−1 )
𝑞(𝑥𝑖 ) = −k ⋅ 𝐴 ⋅
𝑥𝑖+1 − 𝑥𝑖−1
Last point i = N ---> backward difference

𝑇 𝑥𝑁 − 𝑇(𝑥𝑁−1 )
𝑞(𝑥𝑁 ) = −k ⋅ 𝐴 ⋅
𝑥𝑁 − 𝑥𝑁−1

ME 3060 Section 2: Numerical differentiation 17

3. Finite difference method
Comparing numerical (FD) and analytical differentiation for: 𝒇 𝒙 = 𝒙𝟑
Find the numerical error of the forward (f), backward (b),
and central (c) difference scheme at 𝑥2 using the points:
a) 𝑥1 = 2, 𝒙𝟐 = 𝟑, 𝑥3 = 4
b) 𝑥1 = 2.75, 𝒙𝟐 = 𝟑, 𝑥3 = 3.25
The point spacing in a) is Δ𝑥𝑎 = 𝑥𝑖+1 − 𝑥𝑖 = 1 and
in b) Δ𝑥𝑏 = 0.25
---> expect much smaller error in b)
′
Analytical solution: 𝑓 ′ 𝑥 = 3𝑥 2 → 𝑓𝑒𝑥𝑎𝑐𝑡 𝑥2 = 3 = 27
Now let’s find 𝑓 ′ for case a):
Forward difference
𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖 𝑓 4 −𝑓 3
f𝑓′ (𝑥2 = 3) ≈ = = 37
𝑥𝑖+1 − 𝑥𝑖 4−3
27 − 37
Total relative error 𝐸𝑓𝑟𝑒𝑙 = ⋅ 100% = 37.04%
27
ME 3060 Section 2: Numerical differentiation 18
 3. Finite difference method
Backward difference
𝑓 𝑥𝑖 − 𝑓 𝑥𝑖−1 𝑓 3 −𝑓 2
f𝑏′
𝑥2 = 3 ≈ =
𝑥𝑖 − 𝑥𝑖−1 3−2
= 19
𝑟𝑒𝑙 27 − 19
𝐸𝑏 = ⋅ 100% = 29.63%
27
Central difference
𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖−1 𝑓 4 −𝑓 2
f𝑐′
𝑥2 = 3 ≈ = = 28
𝑥𝑖+1 − 𝑥𝑖−1 4−2
𝑟𝑒𝑙
27 − 28
𝐸𝑐 = ⋅ 100% = 3.704%
27
Forward 𝐸𝑓𝑟𝑒𝑙 [%] Backward 𝐸𝒃𝑟𝑒𝑙 [%] Central 𝐸𝒄𝑟𝑒𝑙 [%]
a) Δ𝑥 = 1 37.04 29.63 3.704
Central FD with
b) Δ𝑥 = 0.25 8.565 8.102 0.2315
small Δ𝑥 has
smallest error!
ME 3060 Section 2: Numerical differentiation 19
 3. Finite difference method
Error analysis of two-point FD methods:
Recall Taylor’s theorem
Δ𝑥 2 Δ𝑥 𝑛 Δ𝑥 𝑛+1
𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 ′ 𝑥 Δ𝑥 + 𝑓 ′′ 𝑥 + ⋯+ 𝑓 𝑛 𝑥 +𝑓 𝑛+1 (𝜉)
2! 𝑛! (𝑛 + 1)!
where 𝜉 is a number between 𝑥 and 𝑥 + Δ𝑥; assume 𝜟𝒙 = 𝒙𝒊+𝟏 − 𝒙𝒊 = 𝒄𝒐𝒏𝒔𝒕
Forward difference error (two-point):
Δ𝑥 2
𝑓 𝑥𝑖 + Δ𝑥 = 𝑓(𝑥𝑖+1 ) = 𝑓 𝑥𝑖 + 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝜉
2
𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖 𝑓 ′′ 𝜉
𝑓 ′ 𝑥𝑖 = − Δ𝑥
Δ𝑥 2
′′
′
𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖 𝑓 𝜉
ȁErrorȁ = ȁ𝑓 𝑥𝑖 − ȁ=ȁ Δ𝑥ȁ ⇒ 𝒪(Δ𝑥)
Δ𝑥 2

ME 3060 Section 2: Numerical differentiation 20

3. Finite difference method
Error analysis of two-point FD methods:
Forward difference error (two-point):

𝒪(Δ𝑥) means the “scheme” is first order accurate ---> if we half Δ𝑥

the error will decrease by half
Order of accuracy is typically shown in a log-log plot of Error vs Δ𝑥:

𝐸 ~Δx → 𝐸 = c ⋅ Δ𝑥 1
For some unknown
constant c

𝑙𝑛 𝐸 = ln c + 1 ⋅ ln(Δ𝑥)

𝑓 𝑥 = 𝑥 3 used as example, Matlab script

FD_Analysis.m posted on course webpage

ME 3060 Section 2: Numerical differentiation 21

 3. Finite difference method 𝜟𝒙𝒌 𝚫𝒙𝒏+𝟏
Taylor: 𝒇 𝒙 + 𝚫𝒙 = σ𝒏𝒌=𝟎 𝒇(𝒌) 𝒙 +𝒇 𝒏+𝟏
(𝝃)
Error analysis of two point FD methods 𝒌! (𝒏+𝟏)!

Backward difference error (two-point):

Δ𝑥 2
𝑓 𝑥𝑖 − Δ𝑥 = 𝑓 𝑥𝑖−1 = 𝑓 𝑥𝑖 − 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝜉
2
𝑓 𝑥𝑖 − 𝑓 𝑥𝑖−1 𝑓 ′′ 𝜉
𝑓 ′ 𝑥𝑖 = + Δ𝑥
Δ𝑥 2
′ 𝑓 𝑥𝑖 −𝑓 𝑥𝑖−1 𝑓 ′′ 𝜉
|Errorȁ = ȁ𝑓 𝑥𝑖 − ȁ= ȁ Δ𝑥ȁ ⇒ 𝒪(Δ𝑥)
Δ𝑥 2

Backward difference scheme

is first order accurate

ME 3060 Section 2: Numerical differentiation 22

 3. Finite difference method
Error analysis of two point FD methods
Central difference error (two-point):
2 3
Δ𝑥 Δ𝑥
𝑓 𝑥𝑖+1 = 𝑓 𝑥𝑖 + 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝑥𝑖 + 𝑓′′′(𝜉1 )
2 6
2 3
Δ𝑥 Δ𝑥
𝑓 𝑥𝑖−1 = 𝑓 𝑥𝑖 − 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝑥𝑖 − 𝑓′′′(𝜉2 )
2 6
Δ𝑥 3
𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖−1 = 2𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′′ 𝜉1 + 𝑓′′′(𝜉2 )
6
𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖−1 𝑓 ′′′ 𝜉1 + 𝑓 ′′′ 𝜉2
𝑓 ′ 𝑥𝑖 = − Δ𝑥 2
2 Δ𝑥 12

𝑓 ′′′ 𝜉1 + 𝑓 ′′′ 𝜉2
ȁErrorȁ = ȁ Δ𝑥 2 ȁ ⇒ 𝒪(Δ𝑥 2 )
12

ME 3060 Section 2: Numerical differentiation 23

 3. Finite difference method
Error analysis of two point FD methods
Central difference error (two-point):
𝒪(Δ𝑥 2 ) => central difference is second order accurate: if we
half Δ𝑥 the error will decrease by a quarter

Why slope = 2?
𝑓 ′′′ 𝜉1 + 𝑓 ′′′ 𝜉2
𝐸= Δ𝑥 2
12
𝐸 = 𝑐 ⋅ Δ𝑥 2
𝑙𝑛 𝐸 = ln 𝑐 + 2 ⋅ ln(Δx)
For small Δ𝑥 central difference
(𝒪(Δ𝑥 2 )) is a lot more accurate
Round off errors dominate truncation than forward/backward (𝒪(Δ𝑥))
error! schemes!

ME 3060 Section 2: Numerical differentiation 24

3. Finite difference method
Back to Fin heat transfer example:
We used two-point methods
• Forward difference at 𝑥1 (first order)
• Central difference for 𝑥2 , … , 𝑥𝑁−1 (second order)
• Backward difference at 𝑥𝑁 (first order)

Now we would like to have second order approximations also at the endpoints!
Can we find e.g. a second-order accurate forward difference scheme using more
points i.e. 𝑥𝑖 , 𝑥𝑖+1 , 𝑥𝑖+2 and 𝑓(𝑥𝑖 ), 𝑓(𝑥𝑖+1 ), 𝑓(𝑥𝑖+2 )?

Yes, using Taylor series expansion at 𝑥𝑖+2 and assuming Δ𝑥 = 𝑐𝑜𝑛𝑠𝑡

2Δ𝑥 2 2Δ𝑥 3
𝑓 𝑥𝑖+2 = 𝑓 𝑥𝑖 + 𝑓 ′ 𝑥𝑖 2Δ𝑥 + 𝑓 ′′ 𝑥𝑖 + 𝑓′′′(𝜉2 )
2 6

ME 3060 Section 2: Numerical differentiation 25

 3. Finite difference method
Three-point forward difference (second order)
2 3
Δ𝑥 Δ𝑥
1 𝑓 𝑥𝑖+1 = 𝑓 𝑥𝑖 + 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝑥𝑖 + 𝑓′′′(𝜉1 )
2 6
4Δ𝑥 2 8Δ𝑥 3
2 𝑓 𝑥𝑖+2 = 𝑓 𝑥𝑖 + 2𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝑥𝑖 + 𝑓′′′(𝜉2 )
2 6
Solve for 𝑓′(𝑥𝑖 ) such that the Δ𝑥 2 terms cancel: take 4 ⋅ 1 − (2)
Δ𝑥 3
4𝑓 𝑥𝑖+1 − 𝑓(𝑥𝑖+2 ) = 3𝑓 𝑥𝑖 + 2𝑓 ′ 𝑥𝑖 Δ𝑥 + 4𝑓 ′′′ 𝜉1 − 8𝑓′′′(𝜉2 )
6
4𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖+2 − 3𝑓 𝑥𝑖 Δ𝑥 2
𝑓 ′ 𝑥𝑖 = − 4𝑓 ′′′ 𝜉1 − 8𝑓′′′(𝜉2 )
2Δ𝑥 12
Second order truncation error:
Δ𝑥 2
𝐸𝑟𝑟𝑜𝑟 = 4𝑓 ′′′ 𝜉1 − 8𝑓′′′(𝜉2 ) = 𝒪(Δ𝑥 2 )
12

ME 3060 Section 2: Numerical differentiation 26

 3. Finite difference method
Three-point backward difference (second order)
Using a similar approach, we can derive a second –order accurate three-point
backward difference method:

′
𝑓 𝑥𝑖−2 − 4𝑓 𝑥𝑖−1 + 3𝑓 𝑥𝑖
𝑓 𝑥𝑖 = + 𝒪(Δ𝑥 2 )
2Δ𝑥
Using the three-point schemes at the end
points in our fin example:

How can we further improve the accuracy

of the derivative approximation?
• Temperature measurements are
fixed at 𝑥1 , 𝑥2 , … , 𝑥𝑁 ---> cannot
decrease Δ𝑥
• Need to use higher-order FD
approximations using more points!
ME 3060 Section 2: Numerical differentiation 27
 3. Finite difference method
Fourth-order FD example: mixed five-point forward & backward FD schemes and
four-point central scheme
• First and second points use five-point forward FD based on 𝑥𝑖 , 𝑥𝑖+1 , 𝑥𝑖+2 , 𝑥𝑖+3 , 𝑥𝑖+4
• Points i = 3,N-2 use four-point central FD based on 𝑥𝑖−2 , 𝑥𝑖−1 , 𝑥𝑖+1 , 𝑥𝑖+2
• Last and second to last points use five-point backward FD schemes based on
𝑥𝑖 , 𝑥𝑖−1 , 𝑥𝑖−2 , 𝑥𝑖−3 , 𝑥𝑖−4

For “real” (somewhat noisy) data,

higher order FD approximations often
amplify the noise too much!

ME 3060 Section 2: Numerical differentiation 28

3. Finite difference method
Second derivatives: same concept of using Taylor series approximation can be
applied to derive FD formulas for any order derivative!
Example: Three-point central difference for the second derivative 𝑓′′(𝑥𝑖 ) using
𝑥𝑖−1 , 𝑥𝑖 , 𝑥𝑖+1
Δ𝑥 2 Δ𝑥 3 Δ𝑥 4
1 𝑓 𝑥𝑖+1 = 𝑓 𝑥𝑖 + 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝑥𝑖 + 𝑓 ′′′ 𝑥𝑖 + 𝑓 (4) (𝜉1 )
2 6 24
Δ𝑥 2 Δ𝑥 3 Δ𝑥 4
2 𝑓 𝑥𝑖−1 = 𝑓 𝑥𝑖 − 𝑓 ′ 𝑥𝑖 Δ𝑥 + 𝑓 ′′ 𝑥𝑖 − 𝑓 ′′′ 𝑥𝑖 + 𝑓 (4) (𝜉2 )
2 6 24
Get rid of 𝑓′(𝑥𝑖 ) by adding (1) and (2) --->
Δ𝑥 2 Δ𝑥 4
𝑓 𝑥𝑖+1 + 𝑓 𝑥𝑖−1 = 2𝑓 𝑥𝑖 + 2𝑓 ′′ 𝑥𝑖 + 𝑓 4 𝜉1 + 𝑓 (4) (𝜉2 )
2 24

𝑓 𝑥𝑖+1 − 2𝑓 𝑥𝑖 + 𝑓 𝑥𝑖−1
𝑓 ′′ 𝑥𝑖 = + 𝒪(Δ𝑥 2)
Δ𝑥 2

ME 3060 Section 2: Numerical differentiation 29

3. Finite difference method

General Taylor Expansion (for Unequally Spaced Points):

′ 2
𝑓 ′′ 𝑥𝑖 3
𝑓 ′′′ 𝑥𝑖 𝑛
𝑓 𝑛
𝜉𝑖
𝑓 𝑥𝑖+𝑘 = 𝑓 𝑥𝑖 + 𝑥𝑖+𝑘 − 𝑥𝑖 𝑓 𝑥𝑖 + 𝑥𝑖+𝑘 − 𝑥𝑖 + 𝑥𝑖+𝑘 − 𝑥𝑖 + ⋯ + 𝑥𝑖+𝑘 − 𝑥𝑖
2! 3! 𝑛!
ℎ𝑘

𝑓 ′′ 𝑖 𝑓 ′′′ 𝑖 𝑓 𝑛 𝜉𝑖
𝑓𝑖+𝑘 = 𝑓𝑖 + ℎ𝑘 𝑓 ′ 𝑖 + ℎ𝑘 2 + ℎ𝑘 3 + ⋯ + ℎ𝑘 𝑛
2 6 𝑛!

𝑓𝑖+𝑘 𝑓𝑖 𝑓 ′′ 𝑖 ′′′
2𝑓 𝑖 𝑛−1 𝑓
𝑛
𝜉𝑖
′
𝑓𝑖= − − ℎ𝑘 − ℎ𝑘 − ⋯ − ℎ𝑘 → 𝜀𝑇 𝑓 ′ = 𝑂 ℎ𝑘 𝑛−1
ℎ𝑘 ℎ𝑘 2 6 𝑛!
𝜀𝑇

2 𝑓𝑖+𝑘 2𝑓 ′ 𝑖
2 𝑓𝑖 𝑓 ′′′ 𝑖 2𝑓
4
𝑛−2 𝑓
𝑛
𝜉𝑖
𝑓 ′′
𝑖
= 2 − 2− − 2ℎ𝑘 − 2ℎ𝑘 𝑖
− ⋯ − 2ℎ𝑘 → 𝜀𝑇 𝑓 ′′ = 𝑂 ℎ𝑘 𝑛−2
ℎ𝑘 ℎ𝑘 ℎ𝑘 6 24 𝑛!
𝜀𝑇

𝑓 ′′′ 𝑖 = ⋯ → 𝜀𝑇 𝑓 ′′′ = 𝑂 ℎ𝑘 𝑛−3

In general: For estimating 𝑓 (𝑝) with order of accuracy 𝑛, the truncation error 𝜀𝑇 should be of order 𝑛 + 𝑝

ME 3060 Section 2: Numerical differentiation 30

 3. Finite difference method
Simplified Case: Taylor Expansion for Equally Spaced Points
(𝑙)
𝑓′′𝑖 𝑓′′′𝑖 𝑙 𝑓𝑖
Recall the general form from previous slide: 𝑓𝑖+𝑘 = 𝑓𝑖 + ℎ𝑘 𝑓 ′ 𝑖 + ℎ𝑘 2 2
+ ℎ𝑘 3
6
+ ⋯ = σ∞
𝑙=0 ℎ𝑘 𝑙!
where ℎ𝑘 = 𝑥𝑖+𝑘 − 𝑥𝑖

Now assume equal spacing ℎ = 𝑥𝑖+1 − 𝑥𝑖 :

′′ ′′′ 𝑥
2𝑓
𝑥𝑖 3𝑓
𝑓 𝑥𝑖+1 = 𝑓 𝑥𝑖 + ℎ = 𝑓 𝑥𝑖 + ℎ 𝑓 ′ 𝑥𝑖 + ℎ + ℎ 𝑖
+⋯
2! 3!
𝑓 ′′ 𝑥 𝑓 ′′′ 𝑥
𝑓 𝑥𝑖+2 = 𝑓 𝑥𝑖 + 2ℎ = 𝑓 𝑥𝑖 + 2ℎ 𝑓 ′ 𝑥𝑖 + 2ℎ 2 2! 𝑖 + 2ℎ 3 3! 𝑖 + ⋯
𝑓 ′′ 𝑥𝑖 𝑓 ′′′ 𝑥𝑖
𝑓 𝑥𝑖−3 = 𝑓 𝑥𝑖 − 3ℎ = 𝑓 𝑥𝑖 + −3ℎ 𝑓 ′ 𝑥𝑖 + −3ℎ 2 + −3ℎ 3 +⋯
2! 3!

Therefore in general: ℎ𝑘 = 𝑘 × ℎ where k could be positive or negative integers

∞ (𝑙)
𝑙
𝑓𝑖
𝑓𝑖+𝑘 = ෍ 𝑘ℎ 𝑓𝑜𝑟 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 𝑎𝑛𝑑 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 𝑘
𝑙!
𝑙=0

Alternatively:
∞ (𝑙)
𝑙
𝑓𝑖
෍ 𝑘ℎ 𝑖𝑓 𝑘 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒
𝑙!
𝑙=0
𝑓𝑖+𝑘 = ∞ (𝑙) ∞ 𝑙
𝑓
𝑙 𝑖 𝑙
𝑓𝑖
෍ 𝑘ℎ − ෍ 𝑘ℎ 𝑖𝑓 𝑘 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒
𝑙! 𝑙!
𝑙=0,2,4,… 𝑙=1,3,5,…

ME 3060 Section 2: Numerical differentiation 31

3. Finite difference method

Table 8-1 for equally-spaced points:

In general, with 𝑚 points we get 𝑓 ′ with order of accuracy 𝑛 = 𝑚

ME 3060 Section 2: Numerical differentiation 32

3. Finite difference method

Table 8-1 for equally-spaced points:

For the most part, with 𝑚 points we get 𝑓 ′′ with order of accuracy 𝑛 = 𝑚 − 1, except for central schemes

ME 3060 Section 2: Numerical differentiation 33

3. Finite difference method

Table 8-1 for equally-spaced points:

In general, with 𝑚 points we get 𝑓 ′′′ with order of accuracy 𝑛 = 𝑚 − 2

ME 3060 Section 2: Numerical differentiation 34

3. Finite difference method

Table 8-1 for equally-spaced points:

For the most part, with 𝑚 points we get 𝑓 (4) with order of accuracy 𝑛 = 𝑚 − 3, except for central schemes

𝑚 − 𝑝 + 2 (𝑖𝑓 𝑝 𝑖𝑠 𝑒𝑣𝑒𝑛 𝑎𝑛𝑑 𝑝𝑜𝑖𝑛𝑡𝑠 𝑎𝑟𝑒 𝑐𝑒𝑛𝑡𝑟𝑎𝑙)

In general: with 𝑚 points we get 𝑓 (𝑝) with order of accuracy 𝑛 = ቊ
𝑚−𝑝+1 (𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒)

ME 3060 Section 2: Numerical differentiation 35

3. Finite difference method
Summary:
(𝑙)
𝑓
The general Taylor expansion for unequally spaced points: 𝑓𝑖+𝑘 = σ∞
𝑙=0 ℎ𝑘 𝑙 𝑖𝑙!
where ℎ𝑘 = 𝑥𝑖+𝑘 − 𝑥𝑖

For estimating 𝑓 (𝑝) with order of accuracy 𝑛, the truncation error 𝜀𝑇 should
be of order 𝑛 + 𝑝

Using 𝑚 points, we get 𝑓 (𝑝) with order of accuracy:

𝑚 − 𝑝 + 2 (𝑖𝑓 𝑝 𝑖𝑠 𝑒𝑣𝑒𝑛 𝑎𝑛𝑑 𝑝𝑜𝑖𝑛𝑡𝑠 𝑎𝑟𝑒 𝑐𝑒𝑛𝑡𝑟𝑎𝑙)

𝑛=ቊ
𝑚−𝑝+1 (𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒)

ME 3060 Section 2: Numerical differentiation 36

4. Differentiation with Built-in ‘Symbolic Math’ in Matlab
• Matlab provides powerful symbolic math functionality (it uses the
“Maple” kernel). Note: Symbolic Math Toolbox is required.
• Symbolic variables are defined by ‘syms x f fprime’, where “syms”
tells Matlab that you are defining symbolic variables x, f and
fprime.
• You can do symbolic math in the command window or within a
script file
• Example: running the following script

Get rid of 𝑓′(𝑥𝑖 ) by adding (1) and (2) --->

Δ𝑥 2 Δ𝑥 4
𝑓 𝑥𝑖+1 + 𝑓 𝑥𝑖−1 = 2𝑓 𝑥𝑖 + 2𝑓 ′′ 𝑥𝑖 + 𝑓 4 𝜉1 + 𝑓 (4) (𝜉2 )
2 24
will return:
f = x + x^2*exp(2*x^3)
g = cos(3*x^2)*(5*x + 5*x^2*exp(2*x^3))
ans = 5*x*cos(3*x^2)*(x*exp(2*x^3) + 1)

ME 3060 Section 2: Numerical differentiation 37

4. Differentiation with Built-in ‘Symbolic Math’ in Matlab
Useful commands are (here ‘F’ stands for any symbolic expression):

• ‘simplify(F)’: simplifies the symbolic expression if possible

• ‘collect(F,var)’: rewrites the symbolic expression in powers of var

• ‘diff(F,var)’: differentiates F with respect to var

• ‘int(F,var)’: indefinite integral (antiderivative) of F with respect to

Get rid of 𝑓′(𝑥𝑖 ) by adding (1) and (2) --->
var
2 4
Δ𝑥 Δ𝑥
𝑓 𝑥𝑖+1 + 𝑓 𝑥𝑖−1 = 2𝑓 𝑥𝑖 + 2𝑓 ′′ 𝑥𝑖 + 𝑓 4 𝜉1 + 𝑓 (4) (𝜉2 )
• ‘int(F,var,a,b)’: definite integral from2a to b of F with respect to24var

ME 3060 Section 2: Numerical differentiation 38

 Example 1
Find a finite difference formula for the first derivative 𝑓 ′ 𝑥𝑖 using 𝑓 𝑥𝑖+2 and
𝑓 𝑥𝑖−2 , assuming constant spacing. Show that the formula is 2nd order accurate.
′ ′′ Δ𝑥 2 ′′′ Δ𝑥 3
Taylor: 𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 𝑥 Δ𝑥 + 𝑓 𝑥 2!
+𝑓 𝑥 3!
+⋯
Solution:
Since we want a second-order accurate formula, we need to write down the Taylor expansions up
to the third derivative:
2Δ𝑥 2 2Δ𝑥 3
𝑓𝑖+2 = 𝑓𝑖 + 𝑓𝑖 ′ 2Δ𝑥 + 𝑓𝑖 ′′ + 𝑓𝑖 ′′′
2 6
−2Δ𝑥 2 −2Δ𝑥 3
𝑓𝑖−2 = 𝑓𝑖 + 𝑓𝑖 ′ −2Δ𝑥 + 𝑓𝑖 ′′ + 𝑓𝑖 ′′′
2 6
Now we need to eliminate the second derivative by subtracting the send equation from the first
equation:
8
𝑓𝑖+2 − 𝑓𝑖−2 = 4Δ𝑥𝑓𝑖′ + Δ𝑥 3 𝑓𝑖 ′′′
3
Solve for 𝑓𝑖′ :
𝑓𝑖+2 − 𝑓𝑖−2 2 2
𝑓𝑖′ = − Δ𝑥 𝑓𝑖 ′′′
4Δ𝑥 3
𝑂(Δ𝑥 2 )

ME 3060 Section 2: Numerical differentiation 39

 Example 2

Δ𝑥 2 Δ𝑥 3
Taylor: 𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 ′ 𝑥 Δ𝑥 + 𝑓 ′′ 𝑥 + 𝑓 ′′′ 𝑥 +⋯
2! 3!
Solution:

ME 3060 Section 2: Numerical differentiation 40

 Example 3

′ ′′ Δ𝑥 2 ′′′ Δ𝑥 3
Taylor: 𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 𝑥 Δ𝑥 + 𝑓 𝑥 +𝑓 𝑥 +⋯
2! 3!
Solution:

ME 3060 Section 2: Numerical differentiation 41

 Example 4

Solution:

ME 3060 Section 2: Numerical differentiation 42

 Example 4-Continued

ME 3060 Section 2: Numerical differentiation 43

 Example 5
Given the equally spaced points 𝑥𝑖+1 and 𝑥𝑖−1 , use Taylor series expansions to develop a 2nd order accurate finite difference
formula to evaluate the second derivative 𝑓′′(𝑥) at the point 𝑥 = 𝑥𝑖 .
Solution:
𝑛 = 2; 𝑝 = 2 ⇒ 𝜀𝑇 = 𝑂 𝑛 + 𝑝 = 𝑂 4 → the Taylor expansions should be written up to the 4th derivatives
For point (𝑖 + 1):

′
𝑓 ′′ 𝑖
2 3
𝑓 ′′′ 𝑖 4
𝑓 (4) 𝑖
𝑓𝑖+1 = 𝑓𝑖 + ℎ𝑓 𝑖 + ℎ +ℎ +ℎ +⋯
2 6 24
For point (𝑖 − 1):
𝑓 ′′ 𝑖 𝑓 ′′′ 𝑖 𝑓 (4) 𝑖
𝑓𝑖−1 = 𝑓𝑖 − ℎ𝑓 ′ 𝑖 + ℎ2
− ℎ3 + ℎ4 +⋯
2 6 24
Now, we want to find a formula for 𝑓 ′′ , so we need to eliminate all the smaller derivatives, i.e., 𝑓 ′ in this case. That can be done by
adding the two expansions:

2 ′′
𝑓 (4) 𝑖4
𝑓𝑖+1 + 𝑓𝑖−1 = 2𝑓𝑖 + 0 + ℎ 𝑓 𝑖 + 0 + ℎ +⋯
12
As we can see, the next derivative (𝑓 ′′′ in this case) also cancels out when the points are centrally located around 𝑥𝑖 . Re-arrange to find
𝑓 ′′ 𝑖 we get:

′′
𝑓𝑖+1 + 𝑓𝑖−1 − 2𝑓𝑖 2
𝑓 (4) 𝑖
𝑓 = +ℎ + ⋯ → 𝜀𝑇 = 𝑂(ℎ2 )
𝑖 ℎ2 12
𝜀𝑇

ME 3060 Section 2: Numerical differentiation 44

 Example 6
Given three unequally spaced points (𝑥𝑖 , 𝑦𝑖 ), (𝑥𝑖+1 , 𝑦𝑖+1 ), and (𝑥𝑖−1 , 𝑦𝑖−1 ), use Taylor series expansions to
develop a finite difference formula to evaluate the first derivative 𝑦′(𝑥) at the point 𝑥 = 𝑥𝑖 . After obtaining the
difference formula, show that when the spacing between these points is equal it reduces to the three-point
forward difference formula from Table 8-1:

(𝑙)
𝑓′′𝑖 𝑓′′′𝑖 𝑙 𝑓𝑖
Recall: 𝑓𝑖+𝑘 = 𝑓𝑖 + ℎ𝑘 𝑓 ′ 𝑖 + ℎ𝑘 2 2
+ ℎ𝑘 3 6
+ ⋯ = σ∞
𝑙=0 ℎ𝑘 𝑙!
where ℎ𝑘 = 𝑥𝑖+𝑘 − 𝑥𝑖

Solution
𝑚 = 2; 𝑝 = 1 ⇒ 𝑛 = 𝑚 − 𝑃 + 1 = 2 ⇒ 𝜀𝑇 = 𝑂 𝑛 + 𝑝 = 𝑂 3 → Taylor expansions should be written up to 𝑦 ′′′

ME 3060 Section 2: Numerical differentiation 45

 Example 6-Continued
Given three unequally spaced points (𝑥𝑖 , 𝑦𝑖 ), (𝑥𝑖+1 , 𝑦𝑖+1 ), and (𝑥𝑖−1 , 𝑦𝑖−1 ), use Taylor series expansions to
develop a finite difference formula to evaluate the first derivative 𝑦′(𝑥) at the point 𝑥 = 𝑥𝑖 . After obtaining the
difference formula, show that when the spacing between these points is equal it reduces to the three-point
forward difference formula from Table 8-1:

ME 3060 Section 2: Numerical differentiation 46

 Example 6-Continued
Given three unequally spaced points (𝑥𝑖 , 𝑦𝑖 ), (𝑥𝑖+1 , 𝑦𝑖+1 ), and (𝑥𝑖−1 , 𝑦𝑖−1 ), use Taylor series expansions to
develop a finite difference formula to evaluate the first derivative 𝑦′(𝑥) at the point 𝑥 = 𝑥𝑖 . After obtaining the
difference formula, show that when the spacing between these points is equal it reduces to the three-point
forward difference formula from Table 8-1:

ME 3060 Section 2: Numerical differentiation 47

 Matlab Example 1
• In a script file called ‘SymbolicMath.m’, define symbolic variables ‘x
a f1 f2 f3 f3p f3pi f3piab f3poly’
4
• Set 𝑓1 = 1 + 2 𝑎𝑥 2 + 3sin(2𝑥 2 ) and 𝑓2 = 𝑥 3

𝑓1 ×𝑓2
• Set 𝑓3 = 5
−3
𝑥

• Differentiate 𝑓3 to find 𝑓3 ′ = 𝑓3𝑝 and check it

• Integrate 𝑓3𝑝 to get 𝑓3𝑝𝑖 and verify your result (it should be same as
𝑓3 )
• Integrate 𝑓3𝑝 from xa=0 to xb=1 to get 𝑓3𝑝𝑖𝑎𝑏

• Use collect to get a “polynomial” expression for 𝑓3

ME 3060 Section 2: Numerical differentiation 48

 Matlab Example 2
Part (1) Create a Matlab ‘function’ that calculate the two-point finite difference
approximation based on the equations on slide 16 (mixed forward, central,
backward). You may review and build on the ‘script’ I have provided:
“FD_Analysis.m”. Then follow these steps to write your function:
• The function gets the vectors x and f as input and should return a two-point FD
approximation of the first derivative f’ at all points. So the function header
should be ‘function dfdx = fd_mixed12(x,f)’ .
• Initialize the return argument dfdx: since ‘x’ and ‘f’ are vectors of the same
length ‘N = length(x)’, we can initialize dfdx as a vector with the same number of
elements by writing: ‘dfdx=zeros(N,1)’ .
• Next, we will calculate the derivative for the first point, i.e., we use index ‘i = 1’
in the general formula for the forward difference approximation below:
‘Math notation for
𝑑𝑓 𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖 𝑓𝑖+1 − 𝑓𝑖
ቚ ≈ = discrete values 𝑓1 , 𝑓2 , … for
𝑑𝑥 𝑥=𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖 a vector ‘f’

‘Math notation for an arbitrary analytical

function f that we can evaluate at a point 𝑥𝑖

ME 3060 Section 2: Numerical differentiation 49

 Matlab Example 2
• The first expression on the rhs above is a more ‘mathematical’ form, the second
is actually quite close to Matlab language. The correct line of code should thus
be: ‘dfdx(1) = (f(2)-f(1)) / (x(2)-x(1))’ .
• For the derivative at the last point (‘i=N’) we have to use the backward
difference formula:
𝑑𝑓 𝑓 𝑥𝑖 − 𝑓(𝑥𝑖−1 )
ቚ ≈
𝑑𝑥 𝑥=𝑥𝑖 𝑥𝑖 − 𝑥𝑖−1

• In Matlab, this will be: ‘dfdx(N) = (f(N)-f(N-1))/(x(N)-x(N-1))’

• For all interior points i=2,3,…,N-1, we can use the central difference
approximation:
𝑑𝑓 𝑓 𝑥𝑖+1 − 𝑓(𝑥𝑖−1 )
ቚ ≈
𝑑𝑥 𝑥=𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖−1
• In Matlab, you should use a for loop to calculate all the points:
for i = 2:N-1
dfdx(i) = ….
end
ME 3060 Section 2: Numerical differentiation 50
 Matlab Example 2
Part (2) Open the script “FinHeatTransfer.m”: it provides the basic variables that
you will need to calculate the Fin heat transfer rate for the example we did in class.
Executing the script will create the temperature plot. The script and the resulting
plot are shown below.
%% FinHeatTransfer.m :
clc; clear; close all;
%% Given information
Tr = 20.8; %[C]
k = 60.5; % [W / m K] for steel
% Height, thickness and area of the fin
H = 0.0381; t = 0.00650; A = H*t;
% positions
x = [0;0.025;0.05;0.075;0.1;0.125;0.15;0.175;0.2;0.225;0.25;0.275;0.3;0.325;0.35;0.375;0.4;0.425;0.45];
% Temperature measurements
T = [165.6;132.;107.;89.0;73.2;61.4;52.2;45.8;40.4;35.8;32.8;30.2;28.4;27.0;26.0;25.2;24.6;24.2;23.6];
%% Plot the temperature vs. x positions
% Define a font size
v = 20;
% Define axis line width
alw = 2.0;
% Define plot line and marker widths
lw = 2; mw = 8.;
figure(1)
plot(x,T,'ko','MarkerSize',mw,'LineWidth',lw);
% Set the axis labels (uses LaTex notation)
xlabel('x [m]');
ylabel('T [C]');
% Set the font size and axis line width
set(gca,'FontSize',v);
set(gca,'LineWidth',alw);

ME 3060 Section 2: Numerical differentiation 51

 Matlab Example 2
Now add lines to the script to calculate and plot the heat transfer rate using the
function we wrote in Part (1), ‘fd_mixed12’. Follow these steps:

• Copy the script “FinHeatTransfer.m” and save it as a new script and call it
“FinHeatTransferSolution.m”.
𝑑𝑇
• Calculate the heat transfer rate as: 𝑞 = −k ⋅ 𝐴 ⋅ where the derivative is
calculated by calling the ‘fd_mixed12’ function. 𝑑𝑥

• Plot the heat transfer rate with a black line and black circles. The horizontal
axis label should be ‘x [m]’ and the vertical label should be ‘q [W]’. The
legend for the plotted data should be ‘mixed first/second order accurate’.

ME 3060 Section 2: Numerical differentiation 52

 Matlab Example 2
Part (3) Create a function that calculates the fully second order finite difference
approximation using three-point methods at the endpoints (slides 26 and 27).
Follow these steps:
• Copy the function fd_mixed12(x,f) and save it as a new function called
fd_second(x,f).
• Assume Δx=constant and calculate it as: ‘dx = x(2) -x(1)’.
• Modify the first and the last point calculations by implementing the
second-order three-point forward and backward FD approximations:
4𝑓 𝑥𝑖+1 − 𝑓 𝑥𝑖+2 − 3𝑓 𝑥𝑖 𝑓 𝑥𝑖−2 − 4𝑓 𝑥𝑖−1 + 3𝑓 𝑥𝑖
𝑓′ 𝑥𝑖 = 𝑓 ′ 𝑥𝑖 =
2Δ𝑥 2Δ𝑥

Part (4) Add lines to the “FinHeatTransferSolution.m” script to calculate the heat
transfer rate using fd_second(x,f) function. Then plot the two heat transfer rates in
one plot for comparison. Use a black line for the mixed accuracy and a red line with
red circles for the ‘fully second order accurate’ data series.

ME 3060 Section 2: Numerical differentiation 53

 Matlab Example 3

ME 3060 Section 2: Numerical differentiation 54

 Example 1
Find a finite difference formula for the first derivative 𝑓 ′ 𝑥𝑖 using 𝑓 𝑥𝑖+2 and
𝑓 𝑥𝑖−2 , assuming constant spacing. Show that the formula is 2nd order accurate.
′ ′′ Δ𝑥 2 ′′′ Δ𝑥 3
Taylor: 𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 𝑥 Δ𝑥 + 𝑓 𝑥 2!
+𝑓 𝑥 3!
+⋯

ME 3060 Section 2: Numerical differentiation 55

 Example 2

Δ𝑥 2 Δ𝑥 3
Taylor: 𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 ′ 𝑥 Δ𝑥 + 𝑓 ′′ 𝑥 + 𝑓 ′′′ 𝑥 +⋯
2! 3!

ME 3060 Section 2: Numerical differentiation 56

 Example 3

′ ′′ Δ𝑥 2 ′′′ Δ𝑥 3
Taylor: 𝑓 𝑥 + Δ𝑥 = 𝑓 𝑥 + 𝑓 𝑥 Δ𝑥 + 𝑓 𝑥 +𝑓 𝑥 +⋯
2! 3!

ME 3060 Section 2: Numerical differentiation 57

 Example 4

ME 3060 Section 2: Numerical differentiation 58

 Example 4-Continued

ME 3060 Section 2: Numerical differentiation 59

 Matlab Example 1 - Solution
%% SymbolicMath.m f3 = x^3*(3*sin(2*x^2) +
% Define symbolic variables 2*a*x^2 + 1)
syms x a f1 f2 f3 f3p f3pi f3piab f3poly;

% Define the functions

f1 = 1 + 2*(sqrt(a)*x)^2 + 3*sin(2*x^2);
f3p = 3*x^2*(3*sin(2*x^2) +
2*a*x^2 + 1) + x^3*(4*a*x +
f2 = x^(4/3);
12*x*cos(2*x^2))
f3 = f1 * f2 / x^(-5/3)

% Differentiate to find f3_prime

f3p = diff(f3,x)
f3pi = x^3*(3*sin(2*x^2) +
2*a*x^2 + 1)
% Integrate f3_prime
f3pi = int(f3p,x)

% Integrate f3_prime from a=0 to b=1

f3piab = 2*a + 3*sin(2) + 1
f3piab = int(f3p,0,1)

% Get a polynomial expression for f3 f3poly = 2*a*x^5 +

f3poly=collect(f3,x)
(3*sin(2*x^2) + 1)*x^3

ME 3060 Section 2: Numerical differentiation 60

 Matlab Example 2 – Solution: Part (1)
function dfdx = fd_mixed12(x,f) % Calculates first derivative
using a mixed first/second order approach

% Find the length of the vectors

N = length(x);

% Intialize the output vector

dfdx = zeros(N,1);

% First point -> forward difference

dfdx(1) = (f(2) - f(1))/(x(2)-x(1));

% Last point -> backward difference

dfdx(N) = (f(N) - f(N-1))/(x(N)-x(N-1));

% In between -> central difference (second order)

for i = 2:N-1
dfdx(i) = (f(i+1) - f(i-1))/(x(i+1)-x(i-1));
end
end % End the function

ME 3060 Section 2: Numerical differentiation 61

 Matlab Example 2 – Solution: Part (2)
%% Mixed two-point forward, central, backward difference approximation

q_mixed = -k*A*fd_mixed12(x,T);

figure(2)
plot(x,q_mixed,'ko-','LineWidth',lw);
xlabel('x [m]');
ylabel('q [W]');
set(gca,'FontSize',v);
set(gca,'LineWidth',alw);
legend('mixed first/second order accurate')

ME 3060 Section 2: Numerical differentiation 62

 Matlab Example 2 – Solution: Part (3)
function dfdx = fd_second(x,f)

% Find the length of the vectors

N = length(f);

% Assume dx = const
dx = x(2) -x(1);

% Define array to hold derivatives

dfdx = zeros(N,1);

% First point -> three-point forward difference (second order)

dfdx(1) = (4*f(2) - f(3) - 3*f(1))/(2*dx);

% Last point -> three point backward difference (second order)

dfdx(N) = (f(N-2) - 4*f(N-1) + 3*f(N))/(2*dx);

% In between -> two point central difference (second order)

for i = 2:N-1
dfdx(i) = (f(i+1) - f(i-1))/(2*dx);
end

end

ME 3060 Section 2: Numerical differentiation 63

 Matlab Example 2 – Solution: Part (4)
%% Mixed three-point forward & backward and two
point central
% difference approximation (fully second order
accurate)

q_sec = -k*A*fd_second(x,T);

figure(3)
plot(x,q_sec,'r-
o','MarkerSize',mw,'LineWidth',lw);
hold on;
plot(x,q_mixed,'k-','LineWidth',lw);
hold off;
xlabel('x [m]');
ylabel('q [W]');
set(gca,'FontSize',v);
set(gca,'LineWidth',alw);
legend('fully second order accurate','mixed
first/second order')

ME 3060 Section 2: Numerical differentiation 64

 Matlab Example 3 – Solution
function dfx=DiffAnaly(Fun,xi)
h=0.05*xi;
dfx=(Fun(xi+h)-Fun(xi-h))/(2*h);
end

clc;clear;
Fun_a=@(x) exp(x)*log(x);
dF_a=DiffAnaly(Fun_a,2)

Fun_b=@(x) (x^2+sqrt(x))*cos(x)/sin(x);
dF_b=DiffAnaly(Fun_b,2)

dF_a =

8.8371

dF_b =

-8.6124
ME 3060 Section 2: Numerical differentiation 65