Section 4.

1 Numerical Differentiation

1
Motivation.
𝜕𝜕𝜕𝜕 1 2 2 𝜕𝜕2 𝑓𝑓 𝜕𝜕𝜕𝜕
• Consider to solve Black-Scholes equation + 𝜎𝜎 𝑆𝑆 + 𝑟𝑟𝑟𝑟 −
𝜕𝜕𝜕𝜕 2 𝜕𝜕𝑆𝑆 2 𝜕𝜕𝜕𝜕
𝑟𝑟𝑟𝑟 = 0. Here 𝑓𝑓 is the price of a derivative security, 𝑡𝑡 is time, 𝑆𝑆 is the
varying price of the underlying asset, 𝑟𝑟 is the risk-free interest rate,
and 𝜎𝜎 is the market volatility.
𝜕𝜕𝜕𝜕 𝜕𝜕2 𝑢𝑢 𝜕𝜕2 𝑢𝑢
• The heat equation of a plate: = 𝑘𝑘 � 2 + �. Here 𝑘𝑘 is the heat-
𝜕𝜕𝜕𝜕 𝜕𝜕𝑥𝑥 𝜕𝜕𝑦𝑦 2
diffusivity coefficient.
Goal: Compute accurate approximation to the derivative(s) of a function.
𝑓𝑓(𝑥𝑥0 +ℎ)−𝑓𝑓(𝑥𝑥0 )
The derivative of 𝑓𝑓 at 𝑥𝑥0 is: 𝑓𝑓 ′ (𝑥𝑥0 ) = lim .
ℎ→0 ℎ
𝑓𝑓(𝑥𝑥0 +ℎ)−𝑓𝑓(𝑥𝑥0 )
Obviously, when ℎ is small, is a “good” approximation to
ℎ
𝑓𝑓 ′ (𝑥𝑥0 ).

2
What is the error of approximation?
Big idea:
Build an interpolating polynomial to approximate 𝑓𝑓(𝑥𝑥), then use the
derivative of the interpolating polynomial as the approximation of the
𝑓𝑓 ′ (𝑥𝑥0 ).

Example: Consider to approximate 𝑓𝑓 ′ (𝑥𝑥0 ) using two points 𝑥𝑥0

and 𝑥𝑥0 + ℎ.

3
Example 4.4.1 Use forward difference formula with ℎ = 0.1 to
approximate the derivative of 𝑓𝑓(𝑥𝑥) = ln(𝑥𝑥) at 𝑥𝑥0 = 1.8. Determine the
bound of the approximation error.

𝑓𝑓(𝑥𝑥0 +ℎ)−𝑓𝑓(𝑥𝑥0 )
Forward-difference: 𝑓𝑓 ′ (𝑥𝑥0 ) ≈ when ℎ > 0.
ℎ

𝑓𝑓(𝑥𝑥0 +ℎ)−𝑓𝑓(𝑥𝑥0 )
Backward-difference: 𝑓𝑓 ′ (𝑥𝑥0 ) ≈ when ℎ < 0.
ℎ

4
 1st derivative approximation (obtained by Lagrange interpolation)
The interpolation points are given as:
(𝑥𝑥0 , 𝑓𝑓(𝑥𝑥0 ))
(𝑥𝑥1 , 𝑓𝑓(𝑥𝑥1 ))
(𝑥𝑥2 , 𝑓𝑓(𝑥𝑥2 ))
…
(𝑥𝑥𝑁𝑁 , 𝑓𝑓(𝑥𝑥𝑁𝑁 ))
By Lagrange Interpolation Theorem (Thm 3.3):
(𝑥𝑥−𝑥𝑥0 )⋯(𝑥𝑥−𝑥𝑥𝑁𝑁 ) (𝑁𝑁+1)
𝑓𝑓(𝑥𝑥) = ∑𝑛𝑛𝑘𝑘=0 𝑓𝑓(𝑥𝑥𝑘𝑘 )𝐿𝐿𝑁𝑁,𝑘𝑘 (𝑥𝑥) + (𝑁𝑁+1)!
𝑓𝑓 �𝜉𝜉 (𝑥𝑥)� (1)
Take 1st derivative for Eq. (1):
𝑛𝑛 (𝑁𝑁+1) ( )
(𝑥𝑥 − 𝑥𝑥0 ) ⋯ (𝑥𝑥 − 𝑥𝑥𝑁𝑁 ) 𝑑𝑑 �𝑓𝑓 �𝜉𝜉 𝑥𝑥 ��
′ ′
𝑓𝑓 (𝑥𝑥) = � 𝑓𝑓(𝑥𝑥𝑘𝑘 )𝐿𝐿 𝑁𝑁,𝑘𝑘 (𝑥𝑥) + � �
(𝑁𝑁 + 1)! 𝑑𝑑𝑑𝑑
𝑘𝑘=0
1 𝑑𝑑�(𝑥𝑥 − 𝑥𝑥0 ) ⋯ (𝑥𝑥 − 𝑥𝑥𝑁𝑁 )� (𝑁𝑁+1)
+ � � 𝑓𝑓 �𝜉𝜉 (𝑥𝑥)�
(𝑁𝑁 + 1)! 𝑑𝑑𝑥𝑥
5
Set 𝑥𝑥 = 𝑥𝑥𝑗𝑗 , with 𝑥𝑥𝑗𝑗 being x-coordinate of one of interpolation nodes. 𝑗𝑗 =
0, … , 𝑁𝑁.
𝑓𝑓 (𝑁𝑁+1) �𝜉𝜉�𝑥𝑥𝑗𝑗 ��
𝑓𝑓 ′ �𝑥𝑥𝑗𝑗 � = ∑𝑛𝑛𝑘𝑘=0 𝑓𝑓 (𝑥𝑥𝑘𝑘 )𝐿𝐿′ 𝑁𝑁,𝑘𝑘 �𝑥𝑥𝑗𝑗 �+ ∏𝑁𝑁
𝑘𝑘=0;(𝑥𝑥𝑗𝑗 − 𝑥𝑥𝑘𝑘 )
(𝑁𝑁+1)!
𝑘𝑘≠𝑗𝑗

---- (N+1)-point formula with error to approximate 𝑓𝑓 ′ (𝑥𝑥𝑗𝑗 ) (4.2)

𝑓𝑓 (𝑁𝑁+1) �𝜉𝜉�𝑥𝑥𝑗𝑗 ��
The error of (N+1)-point formula is ∏𝑁𝑁
𝑘𝑘=0;(𝑥𝑥𝑗𝑗 − 𝑥𝑥𝑘𝑘 ).
(𝑁𝑁+1)!
𝑘𝑘≠𝑗𝑗

Remark: 𝑓𝑓 ′ �𝑥𝑥𝑗𝑗 � ≈ ∑𝑛𝑛𝑘𝑘=0 𝑓𝑓 (𝑥𝑥𝑘𝑘 )𝐿𝐿′ 𝑁𝑁,𝑘𝑘 �𝑥𝑥𝑗𝑗 �

6
Example. Derive the three-point formula with error to approximate 𝑓𝑓 ′ (𝑥𝑥𝑗𝑗 ).
Let interpolation nodes be (𝑥𝑥0 , 𝑓𝑓(𝑥𝑥0 )), (𝑥𝑥1 , 𝑓𝑓(𝑥𝑥1 )) and (𝑥𝑥2 , 𝑓𝑓(𝑥𝑥2 )).
′
2𝑥𝑥𝑗𝑗 − 𝑥𝑥1 − 𝑥𝑥2 2𝑥𝑥𝑗𝑗 − 𝑥𝑥0 − 𝑥𝑥2
𝑓𝑓 �𝑥𝑥𝑗𝑗 � = 𝑓𝑓 (𝑥𝑥0 ) � � + 𝑓𝑓(𝑥𝑥1 ) � �
(𝑥𝑥0 − 𝑥𝑥1 )(𝑥𝑥0 − 𝑥𝑥2 ) (𝑥𝑥1 − 𝑥𝑥0 )(𝑥𝑥1 − 𝑥𝑥2 )
2
2𝑥𝑥𝑗𝑗 − 𝑥𝑥0 − 𝑥𝑥1 𝑓𝑓 (3) �𝜉𝜉�𝑥𝑥𝑗𝑗 ��
+ 𝑓𝑓 (𝑥𝑥2 ) � �+ �(𝑥𝑥𝑗𝑗 − 𝑥𝑥𝑘𝑘 )
(𝑥𝑥2 − 𝑥𝑥0 )(𝑥𝑥2 − 𝑥𝑥1 ) 6
𝑘𝑘=0;
𝑘𝑘≠𝑗𝑗

7
Mostly used three-point formula (see Figure 1)
Let 𝑥𝑥0 , 𝑥𝑥1 , and 𝑥𝑥2 be equally spaced and the grid spacing be ℎ.
Thus 𝑥𝑥1 = 𝑥𝑥0 + ℎ; and 𝑥𝑥2 = 𝑥𝑥0 + 2ℎ.

1 ℎ2 (3)
1. 𝑓𝑓 ′ (𝑥𝑥0 ) = [−3𝑓𝑓(𝑥𝑥0 ) + 4𝑓𝑓(𝑥𝑥1 ) − 𝑓𝑓 (𝑥𝑥2 )]
+ 𝑓𝑓 �𝜉𝜉 (𝑥𝑥0 )�
2ℎ 3
(three-point endpoint formula with error) (4.4)
1 ℎ2 (3)
2. 𝑓𝑓 ′ (𝑥𝑥1 ) = [−𝑓𝑓(𝑥𝑥0 ) + 𝑓𝑓(𝑥𝑥2 )] + 𝑓𝑓 �𝜉𝜉 (𝑥𝑥1 )�
2ℎ 6
(three-point midpoint formula with error) (4.5)
1 ℎ2 (3)
3. 𝑓𝑓 ′ (𝑥𝑥2 ) = [𝑓𝑓(𝑥𝑥0 ) − 4𝑓𝑓(𝑥𝑥1 ) + 3𝑓𝑓(𝑥𝑥2 )] +𝑓𝑓 �𝜉𝜉 (𝑥𝑥2 )�
2ℎ 3
(three-point endpoint formula with error) (4.4.1)

8
Remark: Eq. (4.4) in textbook is:
2
1 ℎ
𝑓𝑓 ′ (𝑥𝑥0 ) = [−3𝑓𝑓(𝑥𝑥0 ) + 4𝑓𝑓(𝑥𝑥0 + ℎ) − 𝑓𝑓(𝑥𝑥0 + 2ℎ)] + 𝑓𝑓 (3) �𝜉𝜉 (𝑥𝑥0 )�
2ℎ 3
ℎ can be both positive and negative

9
Mostly used five-point formula
1.Five-point midpoint formula

𝑥𝑥0 − 2ℎ 𝑥𝑥0 − ℎ 𝑥𝑥0 𝑥𝑥0 + ℎ 𝑥𝑥0 + 2ℎ

1
𝑓𝑓 ′ (𝑥𝑥0 ) = [𝑓𝑓(𝑥𝑥0 − 2ℎ) − 8𝑓𝑓(𝑥𝑥0 − ℎ) + 8𝑓𝑓(𝑥𝑥0 + ℎ) − 𝑓𝑓(𝑥𝑥0 + 2ℎ)]
12ℎ
ℎ4 (5)
+ 𝑓𝑓 (𝜉𝜉 ) (𝟒𝟒. 𝟔𝟔)
30

2.Five-point endpoint formula

𝑥𝑥0 𝑥𝑥0 + ℎ 𝑥𝑥0 + 2ℎ 𝑥𝑥0 + 3ℎ 𝑥𝑥0 + 4ℎ

′(
1
𝑓𝑓 𝑥𝑥0 ) = [−25𝑓𝑓(𝑥𝑥0 ) + 48𝑓𝑓(𝑥𝑥0 + ℎ) − 36𝑓𝑓(𝑥𝑥0 + 2ℎ)
12ℎ
ℎ4 (5)
+ 16𝑓𝑓(𝑥𝑥0 + 3ℎ) − 3𝑓𝑓(𝑥𝑥0 + 4ℎ)] + 𝑓𝑓 (𝜉𝜉 ) (𝟒𝟒. 𝟕𝟕)
5 10
Example 4.1.2 Values for 𝑓𝑓(𝑥𝑥) = 𝑥𝑥𝑒𝑒 𝑥𝑥 are given in the following table.
Use all applicable 3-point and 5-point formulas to approximate 𝑓𝑓′(2.0).
𝑥𝑥 1.8 1.9 2.0 2.1 2.2
𝑓𝑓(𝑥𝑥) 10.889365 12.703199 14.778112 17.148957 19.855030

11
 2nd derivative approximation (obtained by Taylor polynomial)

𝑥𝑥0 − ℎ 𝑥𝑥0 𝑥𝑥0 + ℎ

Approximate 𝑓𝑓(𝑥𝑥0 + ℎ) by expansion about 𝑥𝑥0 :

1 1
𝑓𝑓(𝑥𝑥0 + ℎ) = 𝑓𝑓(𝑥𝑥0 ) + 𝑓𝑓 ′ (𝑥𝑥0 )ℎ + 𝑓𝑓 ′′ (𝑥𝑥0 )ℎ2 + 𝑓𝑓 ′′′ (𝑥𝑥0 )ℎ3 +
2 6
1 (4)
𝑓𝑓 (𝜉𝜉1 )ℎ4 (3)
24
Approximate 𝑓𝑓(𝑥𝑥0 − ℎ) by expansion about 𝑥𝑥0 :
1 ′′ 1 ′′′
𝑓𝑓(𝑥𝑥0 − ℎ) = 𝑓𝑓(𝑥𝑥0 ) − 𝑓𝑓 ′ (𝑥𝑥0 )ℎ + 𝑓𝑓 (𝑥𝑥0 )ℎ2 − 𝑓𝑓 (𝑥𝑥0 )ℎ3 +
2 6
1 (4)
𝑓𝑓 (𝜉𝜉2 )ℎ4 (4)
24
Add Eqns. (3) and (4): 𝑓𝑓 (𝑥𝑥0 − ℎ) + 𝑓𝑓(𝑥𝑥0 + ℎ) = 2𝑓𝑓(𝑥𝑥0 ) + 𝑓𝑓 ′′ (𝑥𝑥0 )ℎ2 +
1 (4) 1 (4)
[ 𝑓𝑓 (𝜉𝜉1 )ℎ + 𝑓𝑓 (𝜉𝜉2 )ℎ4 ]
4
24 24
Thus Second derivative midpoint formula
1 ℎ 2
𝑓𝑓 ′′ (𝑥𝑥0 ) = 2 [𝑓𝑓(𝑥𝑥0 − ℎ) − 2𝑓𝑓(𝑥𝑥0 ) + 𝑓𝑓(𝑥𝑥0 + ℎ)] − 𝑓𝑓 (4) (𝜉𝜉 )
ℎ 12 12
Example 3. Values for 𝑓𝑓 (𝑥𝑥) = 𝑥𝑥𝑒𝑒 𝑥𝑥 are given in the following table. Use
second derivative approximation formula to approximate 𝑓𝑓′′(2.0).
𝑥𝑥 1.8 1.9 2.0 2.1 2.2
𝑓𝑓(𝑥𝑥) 10.889365 12.703199 14.778112 17.148957 19.855030

1
Solution: 𝑓𝑓 ′′ (2.0) ≈ (0.1)2
[𝑓𝑓(1.9) − 2𝑓𝑓 (2.0) + 𝑓𝑓(2.1)] =
Or
1
𝑓𝑓 ′′ (2.0) ≈ 2
[𝑓𝑓(1.8) − 2𝑓𝑓(2.0) + 𝑓𝑓(2.2)] =
(0.2)

13
Round-Off Error Instability
Let the round-off errors associated with 𝑓𝑓(𝑥𝑥0 + ℎ) and 𝑓𝑓(𝑥𝑥0 − ℎ) be
𝑒𝑒(𝑥𝑥0 + ℎ) and 𝑒𝑒(𝑥𝑥0 − ℎ), respectively.
Then 𝑓𝑓(𝑥𝑥0 + ℎ) = 𝑓𝑓̃(𝑥𝑥0 + ℎ) + 𝑒𝑒(𝑥𝑥0 + ℎ);
𝑓𝑓 (𝑥𝑥0 − ℎ) = 𝑓𝑓̃(𝑥𝑥0 − ℎ) + 𝑒𝑒(𝑥𝑥0 − ℎ).
Here 𝑓𝑓̃(𝑥𝑥0 + ℎ) and 𝑓𝑓̃(𝑥𝑥0 − ℎ) are actual values used by computer.
The total error of approximation using three-point midpoint formula:
̃(𝑥𝑥0 + ℎ) − 𝑓𝑓̃(𝑥𝑥0 − ℎ) 𝑒𝑒(𝑥𝑥0 + ℎ) − 𝑒𝑒(𝑥𝑥0 − ℎ) ℎ2
𝑓𝑓
𝑓𝑓 ′ (𝑥𝑥0 ) − = − 𝑓𝑓 (3) (𝜉𝜉1 )
2ℎ 2ℎ 6
Assume round-off errors are bounded by 𝜀𝜀 ≥ 0, �𝑓𝑓 (3) (𝜉𝜉1 )� ≤ 𝑀𝑀.
Then:

14
 ̃(𝑥𝑥0 + ℎ) − 𝑓𝑓̃(𝑥𝑥0 − ℎ)
𝑓𝑓 𝜀𝜀 ℎ 2
�𝑓𝑓 ′ (𝑥𝑥0 ) − � ≤ + 𝑀𝑀.
2ℎ ℎ 6
Remark:
𝜀𝜀
1. As ℎ reduces, grows;
ℎ
2.In practice, it’s rare to let ℎ be too small;
𝜀𝜀 ℎ2
3. Let the total error be + 𝑀𝑀, a minimum of the total error occurs at
ℎ 6
1
3𝜀𝜀 3
ℎ=� � .
𝑀𝑀

15