Comprehensive Notes on Numerical

Methods

From Lecture, Tutorial, and Example Sheets

Transcribed and Compiled for Quiz 2

A detailed compilation of topics, including:

• Polynomial and Spline Interpolation
• Least Squares Approximation
• Numerical Differentiation (all formulas, derivations, and error anal-
ysis)
• Numerical Integration (Newton-Cotes, Composite Rules, and error
analysis)

November 10, 2025

ii
Contents

1 Polynomial Interpolation 1
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Lagrange Interpolating Polynomial . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Newton’s Divided Difference Formula . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 Spline Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 Least Squares Approximation 7

2.1 Linear Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Non-linear Data Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 Numerical Differentiation 9
3.1 General Approach: Differentiating Polynomials . . . . . . . . . . . . . . . . . . . 9
3.1.1 Formulas from n = 1 Polynomial (P1 ) . . . . . . . . . . . . . . . . . . . . 9
3.1.2 Formulas from n = 2 Polynomial (P2 ) . . . . . . . . . . . . . . . . . . . . 9
3.2 Taylor Series Derivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2.1 First Derivative Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2.2 Higher-Order First Derivative Formulas . . . . . . . . . . . . . . . . . . . 11
3.2.3 Second Derivative Formulas . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Total Error and Optimal h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.4 Examples and Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Numerical Integration (Quadrature) 15

4.1 Introduction to Newton-Cotes Formulas . . . . . . . . . . . . . . . . . . . . . . 15
4.2 Composite Quadrature Rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

iii
Contents Contents

iv
Polynomial Interpolation

Chapter 1

Polynomial Interpolation

1.1 Introduction
Definition 1.1.1 (Interpolation). Given a set of n + 1 distinct data points (xi , yi ) for i =
0, 1, . . . , n, where yi = f (xi ), the goal of interpolation is to find a polynomial Pn (x) of degree
at most n such that:
Pn (xi ) = yi = f (xi ) for all i = 0, 1, . . . , n
This Pn (x) is called the interpolating polynomial.

Theorem 1.1.2 (Weierstrass Approximation Theorem). Any continuous function f (x) on a

closed interval [a, b] can be approximated by a polynomial to any desired degree of accuracy.

Theorem 1.1.3 (Existence and Uniqueness). There exists a unique polynomial Pn (x) of degree
at most n that interpolates the n + 1 distinct points (x0 , y0 ), . . . , (xn , yn ).

Theorem 1.1.4 (Interpolation Error). Let f ∈ C n+1 [a, b] and let Pn (x) be the unique interpo-
lating polynomial for f (x) at n + 1 distinct nodes x0 , . . . , xn ∈ [a, b]. Then for any x ∈ [a, b],
there exists a ξ(x) ∈ (a, b) such that:
f (n+1) (ξ(x))
f (x) − Pn (x) = (x − x0 )(x − x1 ) . . . (x − xn )
(n + 1)!

1.2 Lagrange Interpolating Polynomial

The Lagrange interpolating polynomial is constructed as a linear combination of basis polyno-
mials Ln,k (x).

Definition 1.2.1 (Lagrange Basis Polynomials). For n + 1 nodes x0 , . . . , xn , the k-th Lagrange
basis polynomial (for k = 0, . . . , n) is defined as:
n
(x − xi )
Ln,k (x) =
Y

i=0,i6=k (xk − xi )

(x − x0 ) . . . (x − xk−1 )(x − xk+1 ) . . . (x − xn )

Ln,k (x) =
(xk − x0 ) . . . (xk − xk−1 )(xk − xk+1 ) . . . (xk − xn )

1
1.3 Newton’s Divided Difference Formula Polynomial Interpolation

1 if k = j
These basis polynomials have the property: Ln,k (xj ) = δkj = .
0 if k =
6 j

Definition 1.2.2 (Lagrange Interpolating Polynomial). The interpolating polynomial Pn (x) is

given by:
n
Pn (x) = f (xk )Ln,k (x)
X

k=0

Example 1.2.3 (Quadratic Lagrange Polynomial, n = 2). Given (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )):

P2 (x) = f (x0 )L2,0 (x) + f (x1 )L2,1 (x) + f (x2 )L2,2 (x)
(x − x1 )(x − x2 ) (x − x0 )(x − x2 )
= f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 )
+ f (x2 )
(x2 − x0 )(x2 − x1 )

Problem 1.2.4. For f (x) = cos x, let x0 = 0.3, x1 = 0.6, x2 = 0.9. Construct an interpolating
polynomial of degree at most 2 to approximate f (0.45).

Solution 1. First, find the function values: f (x0 ) = cos(0.3) = 0.95534 f (x1 ) = cos(0.6) =
0.82534 f (x2 ) = cos(0.9) = 0.62161

The polynomial P2 (x) is (using the formula from the example above):
(0.45 − 0.6)(0.45 − 0.9) (0.45 − 0.3)(0.45 − 0.9)
P2 (0.45) = (0.95534) + (0.82534)
(0.3 − 0.6)(0.3 − 0.9) (0.6 − 0.3)(0.6 − 0.9)
(0.45 − 0.3)(0.45 − 0.6)
+ (0.62161)
(0.9 − 0.3)(0.9 − 0.6)
(−0.15)(−0.45) (0.15)(−0.45) (0.15)(−0.15)
= (0.95534) + (0.82534) + (0.62161)
(−0.3)(−0.6) (0.3)(−0.3) (0.6)(0.3)
0.0675 −0.0675 −0.0225
= (0.95534) + (0.82534) + (0.62161)
0.18 −0.09 0.18
= (0.95534)(0.375) + (0.82534)(0.75) + (0.62161)(−0.125)
= 0.35825 + 0.61900 − 0.07770
= 0.89955

The true value is f (0.45) = cos(0.45) = 0.90045.

1.3 Newton’s Divided Difference Formula

Definition 1.3.1 (Divided Differences). The divided differences are defined recursively:

• Zeroth DD: f [xi ] = f (xi )

f [xj ]−f [xi ]
• First DD: f [xi , xj ] = xj −xi

f [xj ,xk ]−f [xi ,xj ]

• Second DD: f [xi , xj , xk ] = xk −xi

2
Polynomial Interpolation 1.3 Newton’s Divided Difference Formula

• k-th DD: f [x0 , . . . , xk ] = f [x1 ,...,xk ]−f [x0 ,...,xk−1 ]

xk −x0

Definition 1.3.2 (Newton’s Divided Difference Formula). The interpolating polynomial is

given by:

Pn (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 ) + . . .

+ f [x0 , . . . , xn ](x − x0 ) . . . (x − xn−1 )

This can be written as:

Pn (x) = Pn−1 (x) + f [x0 , . . . , n](x − x0 ) . . . (x − xn−1 )

This form is useful for adding new data points.

Problem 1.3.3. Construct the Newton’s forward interpolating polynomial for the following
data:
x 0 1 2 3 4
f (x) 0 1 8 27 64

Solution 2. We construct the divided difference table:

xi f [xi ] 1st DD 2nd DD 3rd DD 4th DD

0 0
1
1 1 3
7 1
2 8 6 0
19 1
3 27 9
37
4 64

Using the forward formula (top diagonal):

P4 (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )

+ f [x0 , x1 , x2 , x3 ](x − x0 )(x − x1 )(x − x2 )
+ f [x0 , . . . , x4 ](x − x0 )(x − x1 )(x − x2 )(x − x3 )
= 0 + 1(x − 0) + 3(x − 0)(x − 1) + 1(x − 0)(x − 1)(x − 2) + 0(x − 0) . . .
= x + 3x(x − 1) + x(x − 1)(x − 2)
= x + 3x2 − 3x + (x2 − x)(x − 2)
= 3x2 − 2x + x3 − 2x2 − x2 + 2x
= x3

Problem 1.3.4. The cubic Newton’s forward polynomial P3 (x) = 2 − (x + 1) + x(x + 1) −

2x(x + 1)(x − 1) interpolates the first 4 points in the table. By adding the point (3, 10), find
the new polynomial P4 (x) and use it to approximate f (0.5).

3
1.4 Spline Interpolation Polynomial Interpolation

Solution 3. The data from P3 (x) corresponds to x0 = −1, x1 = 0, x2 = 1, x3 = 2. We add

x4 = 3. We can use the recursive property:

P4 (x) = P3 (x) + f [x0 , x1 , x2 , x3 , x4 ](x − x0 )(x − x1 )(x − x2 )(x − x3 )

P4 (x) = P3 (x) + C(x + 1)(x − 0)(x − 1)(x − 2)

We need to find C = f [−1, 0, 1, 2, 3]. We build the DD table:

xi f [xi ] 1st 2nd 3rd 4th

-1 2
-1
0 1 1
1 -2
1 2 -3 1.625
-9 6
2 -7 17
17
3 10

The 4th DD is 1.625. So, C = 1.625.

P4 (x) = P3 (x) + 1.625(x + 1)x(x − 1)(x − 2)

P4 (0.5) = P3 (0.5) + 1.625(0.5 + 1)(0.5)(0.5 − 1)(0.5 − 2)
P3 (0.5) = 2 − (0.5 + 1) + (0.5)(0.5 + 1) − 2(0.5)(0.5 + 1)(0.5 − 1)
= 2 − 1.5 + 0.75 − 2(0.5)(1.5)(−0.5)
= 2 − 1.5 + 0.75 + 0.75 = 2.0
P4 (0.5) = 2.0 + 1.625(1.5)(0.5)(−0.5)(−1.5)
= 2.0 + 1.625(1.125)
= 2.0 + 1.828125 = 3.828125

1.4 Spline Interpolation

Splines are piecewise polynomials used to avoid Runge’s phenomenon (high-degree oscillations).

Definition 1.4.1 (Cubic Spline). A function S(x) is a cubic spline interpolant for the data
(x0 , f (x0 )), . . . , (xn , f (xn )) if it satisfies:

1. S(x) is a cubic polynomial, Sj (x), on each subinterval [xj , xj+1 ].

2. S(xj ) = f (xj ) for all j = 0, . . . , n (Interpolation).

3. Sj (xj+1 ) = Sj+1 (xj+1 ) (Continuity of S).

4. Sj0 (xj+1 ) = Sj+1

0
(xj+1 ) (Continuity of S 0 ).

5. Sj00 (xj+1 ) = Sj+1

00
(xj+1 ) (Continuity of S 00 ).

4
Polynomial Interpolation 1.4 Spline Interpolation

This leads to 4n unknowns (4 coefficients for each of n intervals) but only 4n − 2 equations.
We need 2 more "boundary conditions".

Definition 1.4.2 (Natural Cubic Spline). A cubic spline is natural if it satisfies the boundary
conditions:
S 00 (x0 ) = 0 and S 00 (xn ) = 0

Definition 1.4.3 (Clamped Cubic Spline). A cubic spline is clamped if the first derivatives
at the endpoints are specified:

S 0 (x0 ) = f 0 (x0 ) and S 0 (xn ) = f 0 (xn )

5
1.4 Spline Interpolation Polynomial Interpolation

6
Least Squares Approximation

Chapter 2

Least Squares Approximation

2.1 Linear Least Squares

When data is noisy or there are more points than the degree of the desired polynomial, inter-
polation is not ideal. Instead, we find a "best fit" line.

Definition 2.1.1 (Linear Least Squares). Given n data points (xi , yi ), i = 1, . . . , n, we want
to find the line y(x) = a0 + a1 x that "best" fits the data. We define "best" by minimizing the
sum of the squared errors (SSE).

The error for each point is ei = y(xi ) − yi = (a0 + a1 xi ) − yi . We want to minimize the total
error function E(a0 , a1 ):
n n
E(a0 , a1 ) = (y(xi ) − yi ) =2
(a0 + a1 xi − yi )2
X X

i=1 i=1

To find the minimum, we take the partial derivatives with respect to a0 and a1 and set them
to zero:
n
∂E
= 2(a0 + a1 xi − yi )(1) = 0
X
∂a0 i=1
n
∂E
= 2(a0 + a1 xi − yi )(xi ) = 0
X
∂a1 i=1

This simplifies to the Normal Equations:

n n n
1 + a1 xi =
X X X
a0 yi
i=1 i=1 i=1
n n n
x i + a1 x2i =
X X X
a0 xi y i
i=1 i=1 i=1

Let n = 1. In matrix form:

Pn
i=1
" P #" # " P #
n x a0 y
P P 2i = P i
xi xi a1 xi y i

7
2.2 Non-linear Data Linearization Least Squares Approximation

Solving for a0 and a1 (using Cramer’s rule or substitution) gives:

xi yi − ( xi )( yi )
P P P
n
a1 =
n x2i − ( xi )2
P P

x2i yi − xi xi yi
P P P P
a0 = = ȳ − a1 x̄
n x2i − ( xi )2
P P
P P
xi yi
where x̄ = n
and ȳ = n
.

Problem 2.1.2. Construct a linear least square fitting polynomial for the data (0, 6), (1, 0), (2, 0).

Solution 4. We have n = 3. Let’s create a table of sums:

xi yi x2i xi y i
0 6 0 0
1 0 1 0
2 0 4 0
=3 =6 =5 =0
P P P P

The normal equations are:

3a0 + 3a1 = 6
3a0 + 5a1 = 0

From the first equation, a0 + a1 = 2 =⇒ a0 = 2 − a1 . Substitute into the second:

3(2 − a1 ) + 5a1 =0
6 − 3a1 + 5a1 =0
6 + 2a1 =0
a1 = −3

Then, a0 = 2 − (−3) = 5. The best fit line is y(x) = 5 − 3x.

2.2 Non-linear Data Linearization

We can linearize some non-linear models to use the same method.

• Model 1: y = beax Take natural log: ln y = ln b + ax Let Y = ln y, A0 = ln b, A1 = a.

We fit the line Y = A0 + A1 x. Then a = A1 and b = eA0 .

• Model 2: y = bxa Take natural log: ln y = ln b + a ln x Let Y = ln y, X = ln x, A0 = ln b,

A1 = a. We fit the line Y = A0 + A1 X. Then a = A1 and b = eA0 .

8
Numerical Differentiation

Chapter 3

Numerical Differentiation

3.1 General Approach: Differentiating Polynomials

The core idea is to find the polynomial Pn (x) that interpolates n + 1 points, and then use its
derivative Pn0 (x) as an approximation for f 0 (x).
f 0 (xj ) ≈ Pn0 (xj )
The general error formula for an (n + 1)-point formula is:
f (n+1) (ξj ) Y
n
f 0 (xj ) − Pn0 (xj ) = (xj − xk )
(n + 1)! k=0,k6=j

3.1.1 Formulas from n = 1 Polynomial (P1 )

Given two points x0 and x1 = x0 + h. The polynomial is P1 (x) = f (x0 ) xx−x 1
0 −x1
+ f (x1 ) xx−x 0
1 −x0
. The
derivative is constant: P1 (x) = x1 −x0 =
0 f (x1 )−f (x0 ) f (x0 +h)−f (x0 )
h
.

2-Point Forward Difference (at xj = x0 ):

f (x0 + h) − f (x0 )
f 0 (x0 ) ≈ P10 (x0 ) =
h
f 00 (ξ0 )
Error term: 2!
(x0 − x1 ) = − h2 f 00 (ξ0 ).

2-Point Backward Difference (at xj = x1 ):

f (x1 ) − f (x0 ) f (x1 ) − f (x1 − h)
f 0 (x1 ) ≈ P10 (x1 ) = =
h h
f 00 (ξ1 )
Error term: 2!
(x1 − x0 ) = h2 f 00 (ξ1 ).

3.1.2 Formulas from n = 2 Polynomial (P2 )

Given x0 , x1 = x0 + h, x2 = x0 + 2h. The derivative of P2 (x) is:
2x − x1 − x2 2x − x0 − x2 2x − x0 − x1
P20 (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

9
3.2 Taylor Series Derivations Numerical Differentiation

3-Point Endpoint Formula (at xj = x0 ): Substitute x = x0 , x1 = x0 + h, x2 = x0 + 2h into

P20 (x):

2x0 − (x0 + h) − (x0 + 2h) 2x0 − x0 − (x0 + 2h) 2x0 − x0 − (x0 +

f 0 (x0 ) ≈ P20 (x0 ) = f (x0 ) + f (x1 ) + f (x2 )
(−h)(−2h) (h)(−h) (2h)(h)
−3h −2h −h
= f (x0 ) 2 + f (x1 ) 2 + f (x2 ) 2
2h −h 2h
1
= [−3f (x0 ) + 4f (x1 ) − f (x2 )]
2h
h2 000
The error is 3
f (ξ0 ).

3-Point Midpoint Formula (at xj = x1 ): Substitute x = x1 into P20 (x):

2x1 − x1 − x2 2x1 − x0 − x2 2x1 − x0 − x1

f 0 (x1 ) ≈ P20 (x1 ) = f (x0 ) + f (x1 ) + f (x2 )
(−h)(−2h) (h)(−h) (2h)(h)
x1 − x2 (x1 − x0 ) + (x1 − x2 ) x1 − x0
= f (x0 ) + f (x 1 ) + f (x 2 )
2h2 −h2 2h2
−h h + (−h) h
= f (x0 ) 2 + f (x1 ) + f (x2 ) 2
2h −h 2 2h
1 f (x2 ) − f (x0 )
= [−f (x0 ) + 0 · f (x1 ) + f (x2 )] =
2h 2h
f (x1 + h) − f (x1 − h)
f 0 (x1 ) ≈
2h
2
The error is − h6 f 000 (ξ1 ).

3.2 Taylor Series Derivations

3.2.1 First Derivative Formulas
Definition 3.2.1 (Forward Difference). From Taylor expansion f (x + h) = f (x) + hf 0 (x) +
h2 00
2
f (ξ):
f (x + h) − f (x) h 00
f 0 (x) = − f (ξ)
h 2
The forward difference formula is f 0 (x) ≈ f (x+h)−f (x)
h
. The truncation error is E(f ) =
− h2 f 00 (ξ) for ξ ∈ (x, x + h). This is order O(h).
h2 00
Definition 3.2.2 (Backward Difference). From f (x − h) = f (x) − hf 0 (x) + 2
f (ξ):

f (x) − f (x − h) h 00
f 0 (x) = + f (ξ)
h 2

The backward difference formula is f 0 (x) ≈ f (x)−f (x−h)

h
. The truncation error is E(f ) =
f (ξ) for ξ ∈ (x − h, x). This is order O(h).
h 00
2

Definition 3.2.3 (Central Difference). Subtract the Taylor expansions: f (x + h) = f (x) +

2 3 2 3
hf 0 (x) + h2 f 00 (x) + h6 f 000 (ξ1 ) f (x − h) = f (x) − hf 0 (x) + h2 f 00 (x) − h6 f 000 (ξ2 )

10
Numerical Differentiation 3.2 Taylor Series Derivations

h3 f 000 (ξ1 )+f 000 (ξ2 )

f (x + h) − f (x − h) = 2hf 0 (x) + 6
[f 000 (ξ1 ) + f 000 (ξ2 )] By IVT, 2
= f 000 (ξ) for some
ξ ∈ (x − h, x + h).
f (x + h) − f (x − h) h2 000
f 0 (x) = − f (ξ)
2h 6
The central difference formula (or 3-point midpoint) is f 0 (x) ≈ f (x+h)−f (x−h)
2h
. The trunca-
2
tion error is E(f ) = − h6 f 000 (ξ). This is order O(h2 ).

3.2.2 Higher-Order First Derivative Formulas

Definition 3.2.4 (Five-Point Midpoint Formula).
1
f 0 (x) ≈ [f (x − 2h) − 8f (x − h) + 8f (x + h) − f (x + 2h)]
12h

Error Derivation. Expand f (x ± h) and f (x ± 2h) using Taylor series up to h5 :

h2 00 h3 h4 h5 (5)
f (x + h) = f (x) + hf 0 (x) + f (x) + f 000 (x) + f (4) (x) + f (ξ1 )
2 6 24 120
h2 h3 h4 h5 (5)
f (x − h) = f (x) − hf 0 (x) + f 00 (x) − f 000 (x) + f (4) (x) − f (ξ2 )
2 6 24 120
4h2 00 8h3 000 16h4 (4) 32h5 (5)
f (x + 2h) = f (x) + 2hf 0 (x) + f (x) + f (x) + f (x) + f (ξ3 )
2 6 24 120
4h2 00 8h3 000 16h4 (4) 32h5 (5)
f (x − 2h) = f (x) − 2hf 0 (x) + f (x) − f (x) + f (x) − f (ξ4 )
2 6 24 120
Now, compute the combination f (x − 2h) − 8f (x − h) + 8f (x + h) − f (x + 2h):

• f (x) terms: (1 − 8 + 8 − 1)f (x) = 0

• f 0 (x) terms: (−2h − 8(−h) + 8(h) − (2h))f 0 (x) = (−2 + 8 + 8 − 2)hf 0 (x) = 12hf 0 (x)
2 2
• f 00 (x) terms: (2h2 − 8( h2 ) + 8( h2 ) − 2h2 )f 00 (x) = (2 − 4 + 4 − 2)h2 f 00 (x) = 0
3 3 3 8h3
• f 000 (x) terms: (− 8h6 − 8(− h6 ) + 8( h6 ) − 6
)f 000 (x) = (− 43 + 43 + 43 − 43 )h3 f 000 (x) = 0
4 4 4 16h4
• f (4) (x) terms: ( 16h
24
− 8( h24 ) + 8( h24 ) − 24
)f (4) (x) =0

• f (5) terms: ≈ (− 120

32 1
−8(− 120 )+8( 120
1
)− 120
32
)h5 f (5) (ξ) = h f (ξ)
−32+8+8−32 5 (5)
120
= − 120 h f (ξ) =
48 5 (5)

− 25 h5 f (5) (ξ)

So, f (x − 2h) − 8f (x − h) + 8f (x + h) − f (x + 2h) = 12hf 0 (x) − 25 h5 f (5) (ξ) Rearranging for

f 0 (x):
1 1 2 5 (5)
 
f (x) =
0
[f (x − 2h) − · · · − f (x + 2h)] + h f (ξ)
12h 12h 5
2h5 (5) 4
The error is E(f ) = 60h f (ξ) = h30 f (5) (ξ). This is O(h4 ).

Definition 3.2.5 (Five-Point Endpoint Formula).

1
f 0 (x) ≈ [−25f (x) + 48f (x + h) − 36f (x + 2h) + 16f (x + 3h) − 3f (x + 4h)]
12h
h4 (5)
The error is E(f ) = 5
f (ξ). This is O(h4 ).

11
3.3 Total Error and Optimal h Numerical Differentiation

3.2.3 Second Derivative Formulas

Definition 3.2.6 (3-Point Midpoint Formula for f 00 ). We can approximate f 00 (x) by adding
the Taylor expansions for f (x + h) and f (x − h):
h2 00 h3 h4
f (x + h) = f (x) + hf 0 (x) + f (x) + f 000 (x) + f (4) (ξ1 )
2 6 24
2 3
h h h4
f (x − h) = f (x) − hf 0 (x) + f 00 (x) − f 000 (x) + f (4) (ξ2 )
2 6 24
h4
f (x + h) + f (x − h) = 2f (x) + h2 f 00 (x) + 24
[f (4) (ξ1 ) + f (4) (ξ2 )] Rearranging for f 00 (x):
f (x + h) − 2f (x) + f (x − h) h2 (4)
f 00 (x) = − f (ξ)
h2 12
2
The formula is f 00 (x) ≈ f (x+h)−2f (x)+f (x−h)
h2
. The truncation error is E(f ) = − h12 f (4) (ξ). This is
O(h2 ).

Definition 3.2.7 (3-Point Endpoint Formula for f 00 ). We can also find f 00 (x) by applying the
forward difference formula to f 0 (x):
f 0 (x + h) − f 0 (x)
f 00 (x) ≈
h
1 f (x + 2h) − f (x + h) f (x + h) − f (x)
" ! !#
≈ −
h h h
f (x + 2h) − 2f (x + h) + f (x)
=
h2
The error for this formula is O(h).

3.3 Total Error and Optimal h

Numerical differentiation is unstable because of roundoff error.

Definition 3.3.1 (Total Error). The total error is the sum of Truncation Error and Round-
off Error. Let f˜(x) be the computed value of f (x), and |f (x) − f˜(x)| ≤ .

Example 3.3.2 (Optimal h for Forward Difference). The total error e(h) is:
f˜(x + h) − f˜(x)
e(h) = − f 0 (x)
h
f˜(x + h) − f (x + h) + f (x) − f˜(x) f (x + h) − f (x)
= + − f 0 (x)
h h
|ex+h | + |ex |
!
h
≤ + f 0 (x) + f 00 (ξ) − f 0 (x)
h 2
2 hM
≤ + (where |f 00 (ξ)| ≤ M )
h 2
To find the h that minimizes e(h) = 2
h
+ hM
2
:
2 M 4 
r
e (h) = − 2 +
0
= 0 =⇒ h2 = =⇒ hoptimal = 2
h 2 M M

12
Numerical Differentiation 3.4 Examples and Problems

f˜(x+h)−f˜(x−h)
Example 3.3.3 (Optimal h for Central Difference). The total error e(h) is: e(h) = 2h
− f 0 (x) ≤
h2 000 h2 M
2
2h
+ 6
f (ξ) e(h) ≤ 
h
+ 6
(where |f 000 (ξ)| ≤ M ) To find the h that minimizes e(h):

 2hM hM  3
e0 (h) = − + = 0 =⇒ = 2 =⇒ h3 =
h 2 6 3 h M
s
3
hoptimal = 3

3.4 Examples and Problems

Problem 3.4.1. Let f (x) = ln x. (i) Use the 2-point endpoint (forward) and 3-point midpoint
(central) formulas to approximate f 0 (3) using h = 0.1. (ii) Find the absolute error for each
case. (iii) Find the error bound for each case.

Solution 5. Exact derivative: f 0 (x) = 1/x, so f 0 (3) = 1/3 ≈ 0.333333... Derivatives for
bounds: f 00 (x) = −1/x2 , f 000 (x) = 2/x3 , f (4) (x) = −6/x4 .

(i) Approximations

• 2-Point Forward: f 0 (3) ≈ f (3.1)−f (3)

0.1
= ln(3.1)−ln(3)
0.1
= 1.13140−1.09861
0.1
= 0.3279

• 3-Point Midpoint: f 0 (3) ≈ f (3.1)−f (2.9)

2(0.1)
= ln(3.1)−ln(2.9)
0.2
= 1.13140−1.06471
0.2
= 0.33345

(ii) Absolute Errors

• 2-Point: |1/3 − 0.3279| ≈ |0.33333 − 0.3279| ≈ 0.00543

• 3-Point: |1/3 − 0.33345| ≈ |0.33333 − 0.33345| ≈ 0.00012

(iii) Error Bounds

• 2-Point (Forward): Error = |− h2 f 00 (ξ)|, ξ ∈ (3, 3.1) M = maxξ∈[3,3.1] |f 00 (ξ)| = |f 00 (3)| =

1/9. Bound ≤ 0.1
2
· 19 ≈ 0.00556
2
• 3-Point (Midpoint): Error = | − h6 f 000 (ξ)|, ξ ∈ (2.9, 3.1) M = maxξ∈[2.9,3.1] |f 000 (ξ)| =
2
|f 000 (2.9)| = 2/(2.9)3 ≈ 0.0821. Bound ≤ (0.1)
6
· (0.0821) ≈ 0.000137

Problem 3.4.2. Given f (x) = ln x at x0 = 2, x1 = 2.2, x2 = 2.6. Find an approximate value

of f 0 (2) using the method of quadratic interpolation.

Solution 6. We use the derivative of the P2 (x) Lagrange polynomial, evaluated at x = x0 = 2.

The general formula for P20 (x) is:
2x − x1 − x2 2x − x0 − x2 2x − x0 − x1
P20 (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Evaluate at x = x0 = 2:
4 − 2.2 − 2.6 4 − 2 − 2.6 4 − 2 − 2.2
f 0 (2) ≈ P20 (2) = f (2) +f (2.2) +f (2.6)
(2 − 2.2)(2 − 2.6) (2.2 − 2)(2.2 − 2.6) (2.6 − 2)(2.6 − 2.2)

13
3.4 Examples and Problems Numerical Differentiation

−0.8 −0.6 −0.2

P20 (2) = f (2) + f (2.2) + f (2.6)
(−0.2)(−0.6) (0.2)(−0.4) (0.6)(0.4)
−0.8 −0.6 −0.2
     
P20 (2) = (ln 2) + (ln 2.2) + (ln 2.6)
0.12 −0.08 0.24
P2 (2) = (0.69315)(−6.6667) + (0.78846)(7.5) + (0.95551)(−0.8333)
0

P20 (2) = −4.6210 + 5.9134 − 0.7962 = 0.4962

|0.5−0.4962|
Relative Error: Exact f 0 (2) = 1/2 = 0.5. Rel. Error = |0.5|
= 0.0076.

Problem 3.4.3. Use the second derivative 3-point midpoint formula to compute f 00 (3) with
h = 0.1 where f (x) = ln x. Find the absolute error and error bound.

Solution 7. Formula: f 00 (3) ≈ f (2.9)−2f (3)+f (3.1)

(0.1)2

ln(2.9) − 2 ln(3) + ln(3.1) 1.064710 − 2(1.098612) + 1.131402

f 00 (3) ≈ = = −0.1112
0.01 0.01
Exact value: f 00 (x) = −1/x2 , so f 00 (3) = −1/9 = −0.111111... Absolute Error = | − 0.111111 −
(−0.1112)| ≈ 8.9 × 10−5 (Note: Using higher precision f 00 (3) ≈ −0.1111728, Abs. Error
≈ 6.2 × 10−5 )
2
Error Bound: E ≤ h12 maxξ∈[2.9,3.1] |f (4) (ξ)| f (4) (x) = −6/x4 . The max magnitude is at x = 2.9.
2
M = |f (4) (2.9)| = 6/(2.9)4 ≈ 0.0847 Bound ≤ (0.1) 12
(0.0847) ≈ 7.06 × 10−5

Problem 3.4.4. Given the following data, compute f 0 (2.0) using all applicable 2-point, 3-point,
and 5-point formulas.

x 1.8 1.9 2.0 2.1 2.2

f (x) 10.89 12.70 14.77 17.15 19.86

Solution 8. We want f 0 (2.0). The data is equally spaced with h = 0.1.

• 2-Point Forward (h = 0.1): f 0 (2.0) ≈ f (2.1)−f (2.0)

0.1
= 17.15−14.77
0.1
= 23.8

• 2-Point Backward (h = 0.1): f 0 (2.0) ≈ f (2.0)−f (1.9)

0.1
= 14.77−12.70
0.1
= 20.7

• 3-Point Midpoint (h = 0.1): (Most accurate 3-point) f 0 (2.0) ≈ f (2.1)−f (1.9)

2(0.1)
= 17.15−12.70
0.2
=
22.25

• 3-Point Endpoint (h = 0.1, forward): f 0 (2.0) ≈ 2(0.1) 1

[−3f (2.0) + 4f (2.1) − f (2.2)]
= 0.2 [−3(14.77) + 4(17.15) − 19.86] = 0.2 [−44.31 + 68.6 − 19.86] = 4.43
1 1
0.2
= 22.15

• 3-Point Endpoint (h = 0.1, backward): f 0 (2.0) ≈ 2(0.1) 1

[3f (2.0) − 4f (1.9) + f (1.8)]
= 0.2 [3(14.77) − 4(12.70) + 10.89] = 0.2 [44.31 − 50.8 + 10.89] = 4.4
1 1
0.2
= 22.0

• 5-Point Midpoint (h = 0.1): (Most accurate overall) f 0 (2.0) ≈ 12(0.1) 1

[f (1.8)−8f (1.9)+
8f (2.1) − f (2.2)] = 1.2 [10.89 − 8(12.70) + 8(17.15) − 19.86] = 1.2 [10.89 − 101.6 + 137.2 −
1 1

19.86] = 26.63
1.2
≈ 22.19

The 3-point and 5-point formulas all give results around 22.0-22.25, which are likely more
accurate than the 2-point formulas.

14
Numerical Integration (Quadrature)

Chapter 4

Numerical Integration (Quadrature)

4.1 Introduction to Newton-Cotes Formulas

The strategy is to approximate the integrand f (x) with an interpolating polynomial Pn (x) and
then integrate the polynomial.
Z b Z b
f (x)dx ≈ Pn (x)dx
a a

Using the Lagrange form Pn (x) = f (xk )Ln,k (x):

Pn
k=0
Z b Z n
bX n Z b
f (x)dx ≈ f (xk )Ln,k (x)dx = f (xk ) Ln,k (x)dx
X
a a k=0 k=0 |a {z }
weights wk
Z b n
f (x)dx ≈ wk f (xk )
X
a k=0
These are the Newton-Cotes Formulas. The nodes xk are equally spaced.

Definition 4.1.1 (Rectangle Rule, n = 0). P0 (x) = f (x0 ). Let x0 = a.

Z b Z b
f (x)dx ≈ f (a)dx = f (a)(b − a)
a a

f 0 (η)(b−a)2
Error: E = f 0 (ξ)(x − a)dx = f 0 (η) a (x − a)dx = .
Rb Rb
a 2

Definition 4.1.2 (Midpoint Rule, n = 0). Let x0 = a+b

2
.

a+b a+b
! !
Z b Z b
f (x)dx ≈ f dx = f (b − a)
a a 2 2
f 00 (η)(b−a)3
Error: E = 24
. This rule is more accurate.

Definition 4.1.3 (Trapezoidal Rule, n = 1). Use nodes x0 = a, x1 = b. Let h = b − a.

P1 (x) = f (a) x−b
a−b
+ f (b) x−a
b−a
Z b Z b h
f (x)dx ≈ P1 (x)dx = [f (a) + f (b)]
a a 2

15
4.1 Introduction to Newton-Cotes Formulas Numerical Integration (Quadrature)

f 00 (ξ)
Error: f (x) − P1 (x) = 2
(x − a)(x − b).
Z b f 00 (ξ(x))
E= (x − a)(x − b)dx
a 2
Since (x − a)(x − b) ≤ 0 on [a, h b], by i Weighted Mean Value Theorem for Integrals: E =
f 00 (η) R b f 00 (η) (b−a)3
2 a (x − a)(x − b)dx = 2
− 6

h3 00
E = − f (η) for η ∈ (a, b)
12
The rule is exact for polynomials of degree ≤ 1.

Definition 4.1.4 (Simpson’s 1/3 Rule, n = 2). Use nodes x0 = a, x1 = a + h, x2 = a + 2h = b.

Let h = b−a
2
.
Z b Z x2
h
f (x)dx = f (x)dx ≈ [f (x0 ) + 4f (x1 ) + f (x2 )]
a x0 3

Error Derivation for Simpson’s Rule. Instead of integrating P2 (x), we use Taylor expansion
2 3 4
around x1 . f (x) = f (x1 ) + (x − x1 )f 0 (x1 ) + (x−x
2
1)
f 00 (x1 ) + (x−x
6
1)
f 000 (x1 ) + (x−x
24
1)
f (4) (ξ)
Integrate term by term from x0 = x1 − h to x2 = x1 + h:
#x1 +h #x1 +h
(x − x1 )2 (x − x1 )3
" "
Z x1 +h
f (x)dx = [f (x1 ) · x]xx11 +h + f (x1 ) 0
+ f (x1 )00
+ ...
x1 −h
−h
2 x1 −h
6 x1 −h

• 1st term: f (x1 )((x1 + h) − (x1 − h)) = 2hf (x1 )

f 0 (x1 )
• 2nd term: 2
((h)2 − (−h)2 ) = 0
f 00 (x1 ) f 00 (x1 ) h3 00
• 3rd term: 6
((h)3 − (−h)3 ) = 6
(2h3 ) = 3
f (x1 )

• 4th term: (Also integrates to 0)

(x−x1 )4 (4) f (4) (η) R x1 +h f (4) (η) (x−x1 )5 x1 +h

h i
• 5th term (Error): f (ξ)dx x1 −h (x−x1 ) dx = =
R x1 +h 4
x1 −h 24
≈ 24 24 5 x1 −h
f (4) (η) 2h5 (4) h5 (4)
120
(h5 − (−h)5 ) = 120
f (η) = 60
f (η)
h3 00 h5 (4)
So, x0 f (x)dx = 2hf (x1 )+ 3 f (x1 )+ 60 f (η1 ). Now, use the 2nd derivative formula f 00 (x1 ) ≈
R x2

h2 (4)
f (x0 )−2f (x1 )+f (x2 )
h2
− 12 f (η2 ). Substitute this f 00 (x1 ) back into the integral expansion:

h3 f (x0 ) − 2f (x1 ) + f (x2 ) h2 (4)

" #
Z x2 h5 (4)
f (x)dx ≈ 2hf (x1 ) + − f (η2 ) + f (η1 )
x0 3 h2 12 60
Z x2 h h5 h5
f (x)dx ≈ [6f (x1 ) + f (x0 ) − 2f (x1 ) + f (x2 )] − f (4) (η2 ) + f (4) (η1 )
x0 3 36 60
Z x2 h h5
f (x)dx = [f (x0 ) + 4f (x1 ) + f (x2 )] − f (4) (η)
x0 3 90
5
The truncation error is E(f ) = − h90 f (4) (η).

16
Numerical Integration (Quadrature) 4.2 Composite Quadrature Rules

Definition 4.1.5 (Degree of Precision). The degree of precision (DOP) is the largest
integer n such that the formula is exact for all polynomials of degree ≤ n (i.e., f (x) = xk for
all k = 0, . . . , n).

• Trapezoidal Rule: DOP = 1.

• Simpson’s 1/3 Rule: DOP = 3.

Example 4.1.6 (Degree of Precision for Simpson’s Rule). Test f (x) = xk on [0, 2] with h = 1.
x0 = 0, x1 = 1, x2 = 2. Rule: 0 f (x)dx ≈ 13 [f (0) + 4f (1) + f (2)].
R2

• f (x) = 1: 1dx = 2. Rule: 13 [1 + 4(1) + 1] = = 2. (Exact)

R2 6
0 3

2
• f (x) = x: xdx = [ x2 ]20 = 2. Rule: 13 [0 + 4(1) + 2] = = 2. (Exact)
R2 6
0 3

3
• f (x) = x2 : x2 dx = [ x3 ]20 = 83 . Rule: 13 [0 + 4(1)2 + 22 ] = 83 . (Exact)
R2
0

4
• f (x) = x3 : x3 dx = [ x4 ]20 = 4. Rule: 13 [0 + 4(1)3 + 23 ] = = 4. (Exact)
R2 12
0 3

5
• f (x) = x4 : x4 dx = [ x5 ]20 = = 6.4. Rule: 13 [0 + 4(1)4 + 24 ] = ≈ 6.67. (NOT exact)
R2 32 20
0 5 3

The rule is exact for degrees 0, 1, 2, 3, but not 4. Thus, the Degree of Precision is 3.

4.2 Composite Quadrature Rules

To improve accuracy, we divide [a, b] into n subintervals.

Definition 4.2.1 (Composite Trapezoidal Rule). Divide [a, b] into n intervals with h = (b−a)/n
and nodes xj = a + jh.

b n−1
X Z xj+1 n−1
Z
h
f (x)dx = f (x)dx ≈ [f (xj ) + f (xj+1 )]
X
a j=0 xj j=0 2
h
= [(f (x0 ) + f (x1 )) + (f (x1 ) + f (x2 )) + · · · + (f (xn−1 ) + f (xn ))]
2
n−1
h
= [f (x0 ) + 2 f (xj ) + f (xn )]
X
2 j=1

− h12 f 00 (ηj ) = − h12 (nf¯00 ) =

 3
 3
Error: The total error is the sum of n individual errors: Etotal =
Pn
j=1
2
− h12 (nh)f 00 (η)

(b − a)h2 00
Etotal = − f (η)
12

The error is O(h2 ).

Definition 4.2.2 (Composite Simpson’s Rule). Divide [a, b] into n intervals (where n must be

17
4.2 Composite Quadrature Rules Numerical Integration (Quadrature)

even). Let h = (b − a)/n. Apply Simpson’s 1/3 rule on n/2 pairs of intervals [x2j , x2j+2 ].

b n/2−1 Z x n/2−1
Z 2j+2 h
f (x)dx = f (x)dx ≈ [f (x2j ) + 4f (x2j+1 ) + f (x2j+2 )]
X X
a j=0 x2j j=0 3
h
= [(f (x0 ) + 4f (x1 ) + f (x2 )) + (f (x2 ) + 4f (x3 ) + f (x4 )) + . . . ]
3
n/2 n/2−1
h
= [f (x0 ) + 4 f (x2j−1 ) + 2 f (x2j ) + f (xn )]
X X
3 j=1 j=1

The formula separates odd (multiply by 4) and even (multiply by 2) indices. Error: Etotal =
(ηj ) = − h90 ( n2 f ¯(4) ) = − 180
Pn/2  h5 (4)  5 h4
j=1 − 90 f (nh)f (4) (η)

(b − a)h4 (4)
Etotal = − f (η)
180
The error is O(h4 ).

Remark 4.2.3 (Stability). Total error for Composite Simpson’s Rule: e(h) = Roundoff +
4M 4M
Truncation |e(h)| ≤ ni=0 |wi i | + (b−a)h ≈ (b − a) + (b−a)h Unlike differentiation, as
P
180 180
h → 0, the roundoff error is stable and the truncation error goes to 0. Numerical integration is
a stable process.

18