IE2106 Engineering Mathematics I A.

Introduction

Numerical Methods 1. What are Numerical Methods? 5

2. Sources of Errors 9
by 3. Absolute & Relative Errors 13
4. Error Propagation 17
Dr Patricia Wong

Office: S1-B1b-58, Tel: 67904219 B. Solution of Nonlinear Equations

Email: ejywong@[Link]
1. Solving by Iteration 19
2. Fixed-Point Iteration Method 26
3. Newton-Raphson Method 53
PJY Wong 1 PJY Wong 2

4. Speed of Convergence 67 E. Numerical Methods for Differential Equations

1. Euler Method 125

C. Polynomial Interpolation 2. Improved Euler Method (Heun’s Method) 140
3. Runge-Kutta Method 151
1. Lagrange Interpolation 73
4. Systems & Higher Order Equations 160
2. Newton’s Divided Difference Interpolation 88

D. Least Squares Approximation

1. What is Least Squares Approximation? 102

2. Normal Equations 104

PJY Wong 3 PJY Wong 4

 A. Introduction • Analytical methods – What we have learnt so far, for
1. What are Numerical Methods? example,

ax2 + bx + c = 0
Numerical methods are methods for solving problems using √
−b ± b2 − 4ac
computers. Given a mathematical problem, there are two ⇒ x= .
2a
ways to get the solution to the problem.
It gives a solution in closed form. However, in most
• Analytical methods problems we do not know the analytical methods, for
example, solve the following for x
• Numerical methods

PJY Wong 5 PJY Wong

O sin x − x3 = 0.

numerical 6

method

• Numerical methods – Using computer to generate Any value of t would determine the value of y(t).
‘approximate’ solution to the required precision.
On the other hand, numerical method gives y(t) at specific
values of t, for example,
Numerical solution to a problem is quite different from the
analytical solution. As an example, consider the problem
y(0.1) = 1.1787, y(0.2) = 1.3105, y(0.3) = 1.3900, ···.
d2 y
+ 4y = 0, y(0) = 1, y ′(π) = 2. Note that these values are only approximations of the true
dt2
values given by (1).
The analytical solution is

y(t) = sin 2t + cos 2t. (1)

PJY Wong 7 PJY Wong 8
-
cut
off '
2. Sources of Errors Another example, the Taylor series
7
h2
• Truncation Error – This is the error made when a limiting y(x + h) = y(x) + hy ′(x) + y ′′(x) + · · ·
2!
process is truncated before one has reached the limiting
value. For example, the infinite series it truncated to give the approximation

x3 x5 x7 x9 y(x + h) ≈ y(x) + hy ′(x).

sin x = x − + − + − ···
3! 5! 7! 9!
is truncated to give the approximation Here, the truncation error decreases with decreasing step
size h.
x3 x5
sin x ≈ x − + .
3! 5!
PJY Wong
wa 9 PJY Wong 10

first 3 terms

• Round-off Error – This is the error caused by the limited

accuracy of the calculator/computer used to perform the
numerical calculations. It can be reduced by using higher
precision. Decreasing step size h results in larger number
of arithmetic operations and hence higher roundoff error
accumulated.

PJY Wong 11 PJY Wong 12

 3. Absolute & Relative Errors *• A number which is correct to O
n decimal places has an
absolute error
• Let ã be an approximate value of the exact value a. O.
|ϵ| < 0.5 × 10−n (4)
-

• The absolute error in ã is • A number which is correct to n significant digits has a

relative error
ϵ = a − ã. (2) |r| < 0.5 × 101−n. (5)
-

• The relative error in ã is absolute error

ϵ 6
a − ã
r= = . (3)
a a
PJY Wong 13 PJY Wong 14

Example 3.1. Consider a = 13 and ã = 0.333. Clearly, ã is

-
-
O Example 3.2. Consider a = 23.494 and ã = 23.491. Clearly,
correct to 3 decimal places and 3 significant digits. ã is correct to 2 decimal places and 4 significant digits.
- -

|ϵ| = 0.00033... < 0.5 × 10−3 = 0.0005 |ϵ| = 0.003 < 0.5 × 10−2 = 0.005

0.00033... 0.003
|r| = ≈ 0.001 < 0.5 × 101−3 = 0.005 |r| = ≈ 0.00013 < 0.5 × 101−4 = 0.0005
a a

PJY Wong 15 PJY Wong 16

 4. Error Propagation (b) In multiplication and division, the relative error of the result
is bounded approximately by the SUM of the relative errors
Error propagation refers to how errors propagate at each of the terms, i.e., if the relative errors of x̃ and ỹ are
arithmetic operation during computation and affect accuracy. respectively rx and ry , then
Theorem 4.1.
|rx·y | ≤ |rx| + |ry |, |rx/y | ≤ |rx| + |ry |.
(a) In addition and subtraction, the absolute error of the result
is bounded by the SUM of the absolute errors of the terms,
i.e., if the absolute errors of x̃ and ỹ are respectively ϵx and
ϵy , then |ϵx±y | ≤ |ϵx| + |ϵy |.

PJY Wong 17 PJY Wong 18

B. Solution of Nonlinear Equations To find an interval that contains a root of (1):

1. Solving by Iteration Lemma 1.1. Suppose f is continuous in the interval [a, b] and
-

one
number the

-
Consider the equation ave
one number f (a) · f (b) < 0.

O f (x) = 0. (1)
Then, f has at least
- -
a zero in [a, b].

O ∗
We say that x is a zero of f (or a root of equation (1)) if
-
Summary:
f (x∗) = 0.
20
f (a) · f (b) < 0 =⇒ f (x) = 0 has a root in [a, b]
--

PJY Wong 19 PJY Wong 20

 Initial approximations to x∗ can be obtained by sketching values are close together, or a combination of both. The
the function f. We are concerned with numerical methods following theorem will be useful.
for finding x∗ to a prescribed accuracy. The methods
used are iterative, i.e., they give a sequence of numbers Theorem 1.1. (Error estimate) Let x̃ be an approximate root
x0, x1, x2, · · · which converges to x∗ (hopefully). and x∗ a root of (1) and both are located in the interval [a, b].
Suppose f is continuously differentiable on [a, b] and
Important questions to ask:
|f ′(x)| ≥ m > 0 for all x ∈ [a, b]. (2)
(a) Under what conditions will the method converge?
Then,
(b) When to stop the iterative process? The stopping condition |f (x̃)|
|x̃ − x∗| ≤ . (3)
could be that f (x) is close to zero, or that successive x m
PJY Wong 21 PJY Wong 22

Summary: approx root = x̃, exact root = x∗, both are in

[a, b]. Suppose that
Proof. From the mean value theorem, we have
|f ′(x)| ≥ m > 0 for all x ∈ [a, b].
f (x̃) − f (x∗) = (x̃ − x∗)f ′(c) Then,
|f (x̃)|
∗
absolute error = |x̃ − x∗| ≤ .
where c lies between x̃ and x and hence c ∈ (a, b). Since m
f (x∗) = 0, it follows that Application: if we want x̃ to be correct to n decimal places,
just check
|f (x̃)−f (x∗)| = |f (x̃)| = |x̃−x∗|·|f ′(c)| ≥ |x̃−x∗|·m. |f (x̃)|

-
absolute error = |x̃ − x∗| ≤ < 0.5 × 10−n
m

holds
PJY Wong 23 PJY Wong If inequality , 24

means it is correct to

n decimal place
.
 How to get m? 2. Fixed-Point Iteration Method
|f ′(x)| ≥ m > 0 for all x ∈ [a, b].
• First, we transform the equation f (x) = 0 algebraically into
* To get m, there are 2 ways: the form
-

(i) Find the minimum of |f ′(x)| for x ∈ [a, b], take rewrite x = g(x). (4)
-

D m = min |f ′(x)|. x∈[a,b]

to
• Next, choose an x0 and compute x1 = g(x0), x2 = g(x1),
- -

(ii) Find a lower bound of |f ′(x)| for x ∈ [a, b], take and in general
m = this lower bound.
-
e .

S .

x = g(x ), n = 0, 1, 2, · · · . (5)
(f'()(
n+1 n
= x+1 , x z(i , b]
bounds :
*
PJY Wong
lower
25 PJY Wong will it converse to 26

If(2)/ >1 for all 32t[1 3]

,
can have
infinite lower
If '(x))70 . 9
bounds :
-

bourd)
1 Subsc =

min If'()) = 1 +) = 2
(largest lower

x[ [ 3]
,

Note that we may get several different forms of (4), Theorem 2.1. (Convergence of fixed-point iteration)
and the behaviour of corresponding iterative sequences of Let x∗ be a solution of x = g(x) and suppose that g is
x0, x1, · · · may differ accordingly, in particular, with respect continuously differentiable in some interval J containing x∗.
to divergence/convergence and how fast it converges. A Then, if |g ′(x)| ≤ L < 1 in J , the iteration process defined
sufficient condition for convergence is given by the following by (5) converges for any x0 in J .
theorem.
Summary: J is an interval containing a root x∗. If
|g ′(x)| ≤ L < 1 for all x ∈ J ,
then the iteration process xn+1 = g(xn) converges for any x0
in J .

PJY Wong 27 PJY Wong 28

 How to get L? The Steps of Fixed-Point Iteration Method
|g ′(x)| ≤ L for all x ∈ J = [a, b]. Solve f (x) = 0.

* To get L, there are 2 ways: Step 1. Find an interval [a, b] that contains a root of
f (x) = 0. (Note: f (a)f (b) < 0)
(i) Find the maximum of |g ′(x)| for x ∈ [a, b], take
′ Step 2. Change f (x) = 0 to the form x = g(x).
D L = max |g (x)|. x∈[a,b]
* There are many different g(x) possible. How to make sure that the g(x)
′ chosen is such that xn+1 = g(xn) converges to a solution of f (x) = 0?
(ii) Find an upper bound of |g (x)| for x ∈ [a, b], take

e
L = this upper bound.
C ′
* Theorem 2.1. If |g (x)| ≤ L < 1 in [a, b], then x n+1 = g(xn) converges
.
S .
for any x0 in [a, b].
(g'x)) = 3 -

j ,
x[[1 2) ,
check if converge
PJY Wong 29 PJY Wong 30
upper bounds :

↓ findanhaveer 194)[3 for all (2ET1 2) ,

bounds
1g' ( >) / < 3 01
.

:
-
Max 1g'ki)
Submi = 2 3 (smallest upper bound
xE[1 2)
.

Step 3. Check that Step 4. Start the iteration process xn+1 = g(xn) with any
|g ′(x)| ≤ L < 1 for all x ∈ [a, b]. x0 ∈ [a, b].
Step 5. When to stop? Suppose we need the answer to be
* To get L, there are 2 ways:
correct up to 3 decimal places. There are 2 ways:
D (i) Find the maximum of |g ′(x)| for x ∈ [a, b], take (i) Stop when the first 3 decimals of xn stop varying.
L = max |g ′(x)|.
x∈[a,b] (Not so good ...)
(ii) Find an upper bound of |g ′(x)| for x ∈ [a, b], take *(ii) Use Theorem 1.1. First, obtain m such that
L = this upper bound. |f ′(x)| ≥ m > 0 for all x ∈ [a, b].

PJY Wong 31 PJY Wong 32

 D
Then, check to see whether J
|f (xn)|
absolute error = |xn − x∗| ≤ m
< 0.5 × 10−3. redo : until

* LC)

If yes, stop, this xn is correct up to 3 decimal places.

converge?
if
( <I
* To get m, there are 2 ways:
D (i) Find the minimum of |f ′(x)| for x ∈ [a, b], take
m = min |f ′(x)|.
x∈[a,b]

(ii) Find a lower bound of |f ′(x)| for x ∈ [a, b], take

m = this lower bound.

PJY Wong
IN Radians 33 PJY Wong 34

check for courgence

M
1
Example 2.1. Set up an iteration process for f (x) = x − 2 Step 3. With J = [0, 1], let us find L. biggest x
=

3x + 1 = 0. (True solutions: x = 2.618034 and 0.381966) Ex ! !

tre
! 2x ! 2x 2
L= max |g1′ (x)| = max !! !! = max = < 1.
Solution Case (a). Step 1. Since f (0) = 1 > 0 and x∈J -
x∈J 3 x∈J 3 3
f (1) = −1 < 0, take J = [a, b] = [0, 1]. f has a zero in J .
-ve Hence, by Theorem 2.1 the iteration process defined by (6)
f (x) = 0
Step 2. Rewrite the given equation as converges for any x0 ∈ J .
x -
3x + 1 = 0

x = g1(x) = 13 (x2 + 1) x + 1
3x =

[
gives the iterative scheme x = j
(x + 1)
<= (2)
g ,

1
xn+1 = (x2n + 1). (6)
3
PJY Wong 35 PJY Wong 36
 E((5 + 1) =
5(1+ 1) 5 = = 0 . 0667

((d)
=

(2
=
g
,

j
Check
0 66066067
Step 4. (Stopping criterion) Suppose we require the solution converges to the smaller solution. .

to be correct to 2 decimal places. With J = [0,O

-

1], let us find

-
~
n xn |f (xn)|/m < 0 . 5x152X 2
3(x, )H)
remove modulas,
m. -

-
xi
0 1
# If kil I

- -
′
m = min |f (x)| = min |3 − 2x| = min(3 − 2x) = 1. x
=
92(x) 1 0.66666667 0.56 m
0
36
1) +
.

6666667
x∈J x∈J x∈J
5(
- +
Subsc
=
1
=

2 0.48148148 0.21
5(0 --
.

Hence, the stopping criterion is =

= 0
48148148 3 0.41060814
.
0.06
|f (xn)| 4 0.38953301 0.017
= |f (xn)| < 0.5 × 10−2. (7)
m 5 0.38391199 0.00435 < 0.5 × 10−2
Step 5. Choose x0 = 1 ∈ J . Using (6) and (7), we get xn x5 = 0.38 is correct to 2 d.p.
-

↓
jeto 1
PJY Wong 37 PJY Wong 38
,

f(x) x
1
=
-
3x +

1)
g(x) 5(x"
= +

If, however, we choose x0 = 3 ̸∈ J , the sequence is now

I
x0 = 3.000, x1 = 3.333, x2 = 4.037, x3 = 5.766,

which seems to diverge!

···

&

"
x= 1

PJY Wong 39 PJY Wong

32=
0 .

38391199 40
 f(x) 32 3x + 1
-

<" >(2) + 10 ve
L :
converge
f(z)
- = - -

f(a)f(b)20 f (3) = 3" -

3(3) + =
130
tre
I :
diverge
-Case (b). Step 1. Since f (2) < 0 and f (3) > 0, take Step 3. With J = [2, 3], let us find L.
subsi = 2
J = [a, b] = [2, 3]. f has a zero in J . 5 t [ 2 , 3] ! !

O
!1! 1 1
Step 2. Rewrite the given equation as L = max |g2(x)| = max !! 2 !! = max 2 = < 1.
′
x∈J x∈J x x∈J x 4
Make &"
converge
<" 3x + 1 =
1
0)
-

Hence, by Theorem 2.1 the iteration process defined by (8)

x = g2(x) = 3 − subject
x x = 3x
- 1
converges for any x0 ∈ J .
+
giving
1 x =
3
-

j x
,

xn+1 = 3 − . (8) provided

xn !!
s is non-zero

interval [2 , 3]

PJY Wong 41 PJY Wong 42

Step 4. (Stopping criterion) Suppose we require the solution Step 5. Choose x0 = 3 ∈ J . Using (8) and (9), we get xn
to be correct to 3 decimal places. With J = [2, 3], let us find converges to the larger solution.
m.
n Godp
x
n |f (xn)|/m

Jo
m = min |f ′(x)| = min |3 − 2x| = min(2x − 3) = 1. 0 3
x∈J x∈J x∈J
1 2.66666667 0.11
Hence, the stopping criterion is 2 2.62500000 0.015
|f (xn)| 3 2.61904762 0.0023
= |f (xn)| < 0.5 × 10−3. (9)
m 4 2.61818182 0.00033 < 0.5 × 10−3

x4 = 2.618 is correct to 3 d.p.

PJY Wong 43 PJY Wong 44
 Example 2.2. Find a solution of f (x) = x3 + x − 1 = 0 by
iteration. 0
3
> +x -
2

Solution Rewrite the equation as 3)(x+ 1) 1 x)

I
=
x= 1 -

+
x=
>2+ /

"sesion
1
x = g1(x) =
1 + x2

so that
1
xn+1 = . (10)
1 + x2n
PJY Wong 45 PJY Wong 46

The derivative maxIg'l 2|x|

|g1′ (x)| = <1
(1 + x2)2

for any x (see fig), so we have convergence for any x0.

(Steps 1–3 considered done)

PJY Wong 47 PJY Wong 48

 If / ~
can pick any co

Step 4. (Stopping criterion) Suppose we require the solution Step 5. Choose -

x0 = 0.7. Using (10) and (11) give

-
to be correct to 2 decimal places. Let us find m. Clearly,
n Ostp
xn |f (xn)|/m

|f ′(x)| = |3x2 + 1| = 3x2 + 1 ≥ 1 for any x. 0 0.7

1 0.67114094 0.027
Hence, m = 1. 2 0.68945064 0.017
The stopping criterion is 3 0.67780886 0.011
4 0.68520143 0.0069
|f (xn)|
0.00437 < 0.5 × 10−2 ~
−2
= |f (xn)| < 0.5 × 10 . (11) 5 0.68050311
m
x5 = 0.68 is correct to 2 d.p.

PJY Wong 49 PJY Wong 50

The equation can also be written as

of
converge
to pet x = g2(x) = 1 − x3

i intersection so that
|g2′ (x)| = 3x2.

-
Since |g2′ (x)| is greater than 1 near the solution x∗ =

·
0.682328..., we cannot guarantee convergence for the
iterative scheme
xn+1 = 1 − x3n.

PJY Wong 51 PJY Wong 52

 3. Newton-Raphson Method

A
This method is commonly used. Advantages include

TR

TR
simplicity and fast convergence. The underlying idea is that
the graph of f is approximated by suitable tangents.
EX

EX
PJY Wong 53 PJY Wong 54

Step 1. tan β1 = f ′(x0) = f (x 0 ) If it happens that f ′(xn) = 0 for some n, then try another
A
x0−x1
f (x 0 ) x0. Newton-Raphson method converges provided x0 is close
⇒ x1 = x0 − f ′ (x 0 ) enough to x∗. From the figure below we see that if x0 is
TR

Step 2. tan β2 = f ′(x1) = f (x 1 ) chosen from approximately the interval [a, b], then the method
x1−x2
f (x 1 )
should converge.
⇒ x2 = x1 − f ′ (x 1 )
So find a, b such that f (a)f (b) < 0, and pick x0 ∈ [a, b].
Step (n + 1). Newton-Raphson method
EX

f (xn)
xn+1 = xn − , n = 0, 1, 2, · · · . (12)
f ′(xn)

PJY Wong 55 PJY Wong 56

 The Steps of Newton-Raphson Method

A
Solve f (x) = 0.

TR

TR
Step 1. Find an interval [a, b] that contains a root of
f (x) = 0. (Note: f (a)f (b) < 0)
Step 2. Pick x0 ∈ [a, b] that is close to the root. Start the
iteration process
f (xn)
EX

EX
xn+1 = xn − ′ , n = 0, 1, 2, · · · .
f (xn)
Step 3. When to stop? Same stopping criterion as fixed-
point iteration method (see step 5).

PJY Wong 57 PJY Wong 58

Example 3.1. Find a root of the equation x3 + x − 1 = 0 in

A
the interval [0, 1].

Solution We have f (x) = x3 + x − 1 and f ′(x) = 3x2 + 1.

TR
Since f (1) > 0 and f (0) < 0, the equation has a root in
[a, b] = [0, 1].
From (12) we have
x3n + xn − 1

EX
xn+1 = xn − .
3x2n + 1
(Steps 1–2 considered done)
E

PJY Wong 59 PJY Wong 60

(Stopping criterion) Suppose we want the solution to correct Pick x0 = 1 ∈ [a, b], we obtain

A
to 4 decimal places. For x ∈ [0, 1],
n xn f (xn) f ′(xn) |f (xn)|/m

TR
0 1
′ 2
|f (x)| = |3x + 1| ≥ 1 ≡ m. 1 0.75 ··· ··· 0.17
2 0.68604651 ··· ··· 0.009
Hence, the stopping criterion is 3 0.68233958 ··· ··· 0.000028 < 0.5 × 10−4
EX
Hence, x3 = 0.6823 is correct to 4 decimal places.
|f (xn)|
= |f (xn)| < 0.5 × 10−4.
m
PJY Wong 61 PJY Wong 62

Example 3.2. Find the positive solution of 2 sin x = x. (Stopping criterion) Suppose we want the solution to correct

A
$ %
to 4 decimal places. For x ∈ π2 , 2 ,
Solution Setting f (x) = x − 2 sin x, we have f ′(x) = 1 −
" #
2 cos x. Since f π2 < 0 and f (2) > 0, the equation has a π

TR
$ % |f ′(x)| = |1 − 2 cos x| ≥ 1 − 2 cos = 1 ≡ m.
root in [a, b] = π2 , 2 . 2

From (12) we have (See fig)

xn − 2 sin xn Hence, the stopping criterion is

xn+1 = xn − .

EX
1 − 2 cos xn |f (xn)|
= |f (xn)| < 0.5 × 10−4.
(Steps 1–2 considered done) m

PJY Wong 63 PJY Wong 64

 Pick x0 = 2 ∈ [a, b], we obtain

A
n xn f (xn) f ′(xn) |f (xn)|/m

TR
0 2
1 1.90099559 ··· ··· 0.009
2 1.89551165 ··· ··· 0.0000285 < 0.5 × 10−4

EX
Hence, x3 = 1.8955 is correct to 4 decimal places.

PJY Wong 65 PJY Wong 66

4. Speed of Convergence Suppose that g is sufficiently differentiable, then the Taylor

A

A
formula gives
The quality of an iteration method can be judged by its
TR

TR
speed of convergence. xn+1 = g(xn)
= g(x∗ + xn − x∗)
Let xn+1 = g(xn) define an iteration method. Let xn
= g(x∗) + g ′(x∗)(xn − x∗) + 12 g ′′(x∗)(xn − x∗)2 + · · ·
approximate a solution x∗ of x = g(x). Then,
= g(x∗) − g ′(x∗)ϵn + 12 g ′′(x∗)ϵ2n + · · · . (13)
EX

EX
xn = x∗ − ϵn,
The exponent of ϵn in the first non-vanishing term after
where ϵn is the error of xn. g(x∗) is called the order of the iteration process defined by
g. The order measures the speed of convergence.
PJY Wong 67 PJY Wong 68
 Subtract x∗ (= g(x∗)) from both sides of equation (13), of convergence. Thus,

A
we have
ϵn+1 ≈ g ′(x∗)ϵn in the case of first order

TR

TR
1
xn+1 − x∗ = −g ′(x∗)ϵn + g ′′(x∗)ϵ2n + · · ·
2 1
ϵn+1 ≈ − g ′′(x∗)ϵ2n in the case of second order.
2
1
− ϵn+1 = −g ′(x∗)ϵn + g ′′(x∗)ϵ2n + · · ·
⇒ (14)
2
EX

EX
where ϵn+1 is the error of xn+1.
On the RHS of (14), the remaining expression approximately
equals its first nonzero term because |ϵn| is small in the case
PJY Wong 69 PJY Wong 70

Example 4.1. Consider Newton-Raphson method. We have Note that

A
f ′′(x∗)
g ′′(x∗) =
f (x) f ′(x∗)
g(x) = x − which will not be zero in general. Hence,
f ′(x)

TR
and by differentiation, 1 f ′′(x∗)
ϵn+1 ≈ − g ′′(x∗)ϵ2n = − ′ ∗ ϵ2n.
2 2f (x )
f ′(x)2 − f (x)f ′′(x) f (x)f ′′(x)
g ′(x) = 1 − = .
f ′(x)2 f ′(x)2

EX
Since f (x∗) = 0, we get g ′(x∗) = 0. Thus Newton-Raphson’s
method is at least of second order.

PJY Wong 71 PJY Wong 72

 /
+ Pp(x)
xX Is degre of s

X one less than

> no .
of pts

&
C. Polynomial Interpolation of x between the points x0, x1, · · · , xn?
1. Lagrange Interpolation In polynomial interpolation, we find a polynomial pn(x)
of degree n (or less) that assumes the values f0, f1, · · · , fn
Suppose we are given a set of (n + 1) points x0, x1, · · · , xn at x0, x1, · · · , xn, i.e.,
and corresponding values f0, f1, · · · , fn of a function f at
pn(x) = anxn + an−1xn−1 + · · · + a1x + a0 (1)
these points.
x x0 x1 · · · xn fulfills
pn(x0) = f0, pn(x1) = f1, · · · , pn(xn) = fn
f (x) f0 f1 · · · fn

Base on the data given, how to construct an approximation or equivalently

to f so that we can find approximate values of f for values pn(xi) = fi, 0 ≤ i ≤ n. (2)
PJY Wong 73 PJY Wong 74

Theorem 1.1. If the points x0, x1, · · · , xn are distinct, then Note that (2) yields a system of (n+1) equations in (n+1)
there exists a unique polynomial pn(x) of degree n (or less) unknowns an, an−1, · · · , a0. By solving the system, we can
passing through the points (x0, f0), (x1, f1), · · · , (xn, fn). obtain the constants an, an−1, · · · , a0 and hence pn(x).

The system of equations (2) can be solved in principle.

The polynomial pn is called the interpolating polynomial
But practically it is not so easy. Therefore, instead of writing
of f with respect to the nodes x0, x1, · · · , xn .
pn(x) as in (1), we express it in terms of Lagrange polynomials
l0(x), l1(x), · · · , ln(x).
Using pn to get approximate values of f for x’s between
x0 and xn is interpolation, whereas outside this interval is What are Lagrange polynomials?
extrapolation.

PJY Wong 75 PJY Wong 76

What are Lagrange polynomials l0(x), l1(x), · · · , ln(x)? From (4), we can see that

A
Each of l0(x), l1(x), · · · , ln(x) is of degree n, and pn(x0) = f0, pn(x1) = f1, · · · , pn(xn) = fn.

TR

TR
⎧ ⎫
⎨ 1, if i = j ⎬ Therefore, pn(x) given in (4) is indeed the interpolating
lj (xi) = 0 ≤ i, j ≤ n. (3)
⎩ 0, if i ̸= j ⎭ polynomial of f with respect to the nodes x0, x1, · · · , xn.
Since pn(x) given in (4) is in terms of Lagrange
The interpolating polynomial pn(x) can be expressed as polynomials, it is called the Lagrange interpolating
EX

EX
,n polynomial of f .
pn(x) = f0l0(x) + f1l1(x) + · · · + fnln(x) = k=0 fklk(x).
(4)

PJY Wong 77 PJY Wong 78

What is the expression of Lagrange polynomial lj (x), 0 ≤ Theorem 1.2. (Error estimate) Suppose f is (n + 1) times
A
j ≤ n? continuously differentiable. Then,

(x−x0)(x−x1)···(x−xj−1)(x−xj +1)···(x−xn)
zoi
TR

lj (x) = (xj −x0)(xj −x1)···(xj −xj−1)(xj −xj +1)···(xj −xn) f (n+1) n

(ξ) -
f (x) − pn(x) = (x − xk) (6)
(n + 1)! k=0
n
- (x−xk ) (5)
= (xj −xk )
. where
k=0
EX

k ̸= j
min{x, x0, x1, · · · , xn} ≤ ξ ≤ max{x, x0, x1, · · · , xn}.
(7)
Moreover, from (6) we have the error inequality
PJY Wong 79 PJY Wong 80
 ! n !
!- 1 ! Summary
! !

A
|f (x) − pn(x)| ≤ ! (x − x k !×
)
-
-

(n + 1)!
-
-

! ! -
-

k=0 1. Given (n + 1) points (xi, fi), i = 0, 1, · · · , n, the Lagrange

|f (n+1)(ξ)|.

TR
max (8) interpolating polynomial pn(x) can be expressed as (4)
min{x,x0,x1,···,xn}≤ξ≤max{x,x0,x1,···,xn}
=
-
-

,n
pn(x) = f0l0(x) + f1l1(x) + · · · + fnln(x) = k=0 fklk(x)
** Think 1.1. If f is a polynomial of degree n, then what can
you observe from (6)? where the Lagrange polynomials lj (x) have the expression

EX
(5).
** Think 1.2. What do you think of the error |f (x) −
pn(x)| if x ∈ [x0, xn] (interpolation) and if x ̸∈ [x0, xn] 2. The error inequality is given by (8).
(extrapolation)?
PJY Wong 81 PJY Wong 82

Example 1.1. Given The interpolating polynomial p1(x) is

A

A
ln 9.0 = 2.1972 and ln 9.5 = 2.2513,
p1(x) = f0l0(x)+f1l1(x) = 2.1972×2(9.5−x)+2.2513×2(x−9).
find the interpolating polynomial of ln x. Hence, compute an
TR

TR
approximate value of ln 9.2.
Further,
Solution From the data, we have
ln 9.2 ≈ p1(9.2) = 2.1972 · 2(9.5 − 9.2) + 2.2513 · 2(9.2 − 9)
n = 1, x0 = 9.0, x1 = 9.5, f0 = 2.1972, f1 = 2.2513.
= 2.2188.
EX

EX
We write down the Lagrange polynomials
(x−x1)
l0(x) = (x0−x1)
= x− 9.5
−0.5
= 2(9.5 − x),
(x−x0)
l1(x) = (x1−x0) = x− 9.0
0.5
= 2(x − 9).

PJY Wong 83 PJY Wong 84

Example 1.2. Given the data of Example 1.1 plus one more We write down the Lagrange polynomials

A
point l0(x) = (x−9.5)(x−11)
,
(9.0−9.5)(9.0−11.0)
ln 9.0 = 2.1972, ln 9.5 = 2.2513, ln 11.0 = 2.3979, (x−9.0)(x−11)
l1(x) = ,

TR

TR
(9.5−9.0)(9.5−11.0)
compute an approximate value of ln 9.2 and obtain an error (x−9.0)(x−9.5)
l2(x) = (11.0−9.0)(11.0−9.5)
.
estimate.
The interpolating polynomial p2(x) is
Solution From the date, we have
EX

EX
n = 2, x0 = 9.0, x1 = 9.5, x2 = 11.0, p2(x) = f0l0(x) + f1l1(x) + f2l2(x).

f0 = 2.1972, f1 = 2.2513, f2 = 2.3979. Hence,

ln 9.2 ≈ p2(9.2) = 2.2192.

PJY Wong 85 PJY Wong 86

Here, f (x) = ln x, f (n+1)(x) = f (3)(x) = 2

, 2. Newton’s Divided Difference Interpolation
A
x3

2 2 2 Suppose we are given a set of (n + 1) points x0, x1, · · · , xn

max |f (3)(ξ)| = max = = .
ξ∈[9,11] ξ 3 93 729
TR

ξ∈[9,11] and corresponding values f0, f1, · · · , fn of a function f at

Using the error inequality (8) we have these points.
x x0 x1 · · · xn
1 2
| ln 9.2 − p2(9.2)| ≤ |(9.2 − 9)(9.2 − 9.5)(9.2 − 11)| f (x) f0 f1 · · · fn
3! 729
EX

Once again, we shall find the interpolating polynomial pn(x)

= 0.4938 × 10−4 < 0.5 × 10−4. of f, this time, via divided differences.

Therefore, p2(9.2) = 2.2192 is correct to 4 decimal places. What are divided differences?

PJY Wong 87 PJY Wong 88

• First order divided difference of f Derive the interpolating polynomial of f
-
f (z1)−f (z0) f (x0)−f (x)
f [z0, z1] = z1−z0 First, f [x, x0] =
i
x0−x

• Second order divided difference of f ⇒ f (x) = f (x0) + (x − x0)f [x, x0]. (9)
f (z1,z2)−f (z0,z1)
f [z0, z1, z2] = z2−z0
f [x0,x1]−f [x,x0]
Next, f [x, x0, x1] = x1−x
• Third order divided difference of f
f (z1,z2,z3)−f (z0,z1,z2)
f [z0, z1, z2, z3] = z3−z0
⇒ f [x, x0] = f [x0, x1] + (x − x1)f [x, x0, x1]. (10)

• kth order divided difference of f Substitute (10) into (9) gives

f (z1,z2,···,zk )−f (z0,z1,···,zk−1)
f [z0, z1, · · · , zk] = zk −z0
f (x) = f (x0)+(x−x0)f [x0, x1]+(x−x0)(x−x1)f [x, x0, x1].
PJY Wong 89 PJY Wong 90

Continuing in this way, we get where

Remainder
f (x) Rn(x) = (x − x0)(x − x1) · · · (x − xn)f [x, x0, x1, · · · , xn]

-
= f (x0) + (x − x0)f [x0, x1] (12)
+(x − x0)(x − x1)f [x0, x1, x2] + · · · and
+(x − x0)(x − x1) · · · (x − xn−1)f [x0, x1, · · · , xn]
Interpolating Polynomail
pn(x)
+ (x − x0)(x − x1) · · · (x − xn)f [x, x0, x1, · · · , xn] = f (x0) + (x − x0)f [x0, x1]
-

= pn(x) + Rn(x). +(x − x0)(x − x1)f [x0, x1, x2] + · · ·

(11) -

-
+(x − x0)(x − x1) · · · (x − xn−1)f [x0, x1, · · · , xn].
-
(13)
Divided
PJY Wong 91 PJY Wong

& differences
92
 S
It can be shown that pn(xi) = f (xi), i = 0, 1, · · · , n,

A
therefore pn(x) given in (13) is indeed the interpolating
polynomial of f with respect to the nodes x0, x1, · · · , xn.

TR
Since pn(x) given in (13) is in terms of divided Lagrange vs Newton
differences, it is called the Newton’s divided difference
interpolating polynomial of f .

EX
Remark 2.1. Is (13) different from (4)? NO, they are the
same pn(x) expressed in different forms. After simplifying,
the final expression is the same.

PJY Wong 93 PJY Wong 94

Remark 2.2. Suppose we add one more point (xn+1, fn+1) Remark 2.3. The computations of the finite differences can
to the original (n + 1) points (xi, fi), i = 0, 1, · · · , n. Then, be simplified by using the following table.
in (13) the first (n + 1) terms remain the same and we just
need to add one more term
(x − x0)(x − x1) · · · (x − xn)f [x0, x1, · · · , xn+1].
This is of great advantage over Lagrange’s representation (4)
where we need to re-calculate all the l0(x), l1(x), · · · , ln(x), in
addition to ln+1(x). Thus, if we seek better approximation by
adding more points, we can still use our previous calculations
in Newton’s representation (13).

PJY Wong 95 PJY Wong 96

 E
-coordinates
-coordinates
in

I
,
i
Summary
xi f [·] f [·, ·] f [·, ·, ·] f [·, ·, ·, ·]

1. Given (n+1) points (xi, fi), i = 0, 1, · · · , n, the Newton’s

x0 Q f [x0]
divided difference interpolating polynomial pn(x) can be
-
②-Q -
x,
- Po
expressed as (13)
x1 ②f [x1] f [x0, x1]
-③ ②
x
-

3,
-
6 x To
-
-
8
pn(x)

-
-

x2 ③f [x2] f [x1, x2] f [x0, x1, x2] = f (x0) + (x − x0)f [x0, x1]
No
-B ⑰ -6
-

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

L
-

xy xz -

Jy-
3,

x3 ④f [x3] f [x2, x3] f [x1, x2, x3] f [x0, x1, x2, x3] +(x − x0)(x − x1) · · · (x − xn−1)f [x0, x1, · · · , xn].

PJY Wong 97 PJY Wong 98

2. The error inequality is given by (8). Solution The divided difference table is as follows.

3. Newton’s representation of pn(x) is better than Lagrange’s i xi f [·] f [·, ·] f [·, ·, ·] f [·, ·, ·, ·]
representation when new data are added. 0 −1 6

·1
1 0 1 0 -(
-
1) −5
Example 2.1. With a divided difference table, find the third
degree interpolating polynomial p3(x) passing through the
-
-
2
3
2
5
3
66 s 21
2

5 4
- (1) 2
* =
1
3
I
points (−1, 6), (0, 1), (2, 3), (5, 66).
↓ and order

T . 1-15
pts it order eg s -
21
y-coordinates
-

2
ald
x-coordinates di
VI

PJY Wong 99 PJY Wong

difference 100
-3
66 = 2)
2 ]
. -

-
2
5 -
 From (13), the interpolating polynomial is D. Least Squares Approximation
1. What is Least Squares Approximation?
p3(x) = f (x0) + (x − x0)f [x0, x1]
+ (x − x0)(x − x1)f [x0, x1, x2] It is not always advisable to construct interpolation
+ (x − x0)(x − x1)(x − x2)f [x0, x1, x2, x3] polynomial passing through all the points (xi, yi), i =
= 6 + (x + 1)(−5) + (x + 1)(x − 0)(2) 0, 1, · · · , n, in particular when n is large, or when the data is
+ (x + 1)(x − 0)(x − 2) 13 likely to contain error. In such a case, to find a best fit curve
through a set of points is the method of least squares.
= 6 − 5(x + 1) + 2x(x + 1) + 13 x(x + 1)(x − 2).

PJY Wong 101 PJY Wong 102

Here, we shall construct a polynomial pm of degree at most In other words, we plan to minimize the sum of the square of
A

A
m (where m < n), errors at the points (xi, yi), i = 0, 1, · · · , n. Hence, the term
‘least squares’. The polynomial pm is called the least-square
TR

TR
m
.
m m−1
polynomial of degree m.
pm(x) = amx + am−1x + · · · + a1 x + a0 = aj xj ,
j =0
(1)
2. Normal Equations
such that the following sum is minimized To minimize the sum of the square of errors S (see (2)) at
EX

EX
the points (xi, yi), i = 0, 1, · · · , n, we differentiate S w.r.t.
⎛ ⎞2 ∂S
n
. n
. m
. ak, k = 0, 1, · · · , m, and set ∂a = 0.
aj xji − yi⎠ .
k
S= [pm(xi) − yi]2 = ⎝ (2)
i=0 i=0 j =0

PJY Wong 103 PJY Wong 104

 ⎛ ⎞
n
. m
. or equivalently

A
∂S
=2 ⎝ aj xji − yi⎠ xki = 0, 0 ≤ k ≤ m. (3) ⎧
∂ak ⎪ s0 a 0 + s1 a 1 + · · · + sm a m = v 0
i=0 j =0 ⎪
⎪
⎪
⎪

TR
Using the notations ⎨ s a + s a + ··· + s a =v
n n 1 0 2 1 m+1 m 1
. . (5)
sk = xki, vk = yixki, ⎪
⎪ ·········
⎪
⎪
i=0 i=0 ⎪
⎩ s a +s
m 0 m+1a1 + · · · + s2mam = vm.
the system (3) can be written as
These equations are called normal equations.
EX
m
.
a j sj + k − v k = 0
j =0 (4)
,m Note: s0 = n + 1 = the no. of points given
⇒ j =0 aj sj +k = vk , 0 ≤ k ≤ m

PJY Wong 105 PJY Wong 106

There are two things to be shown: (i) To show that the system (4) has a unique solution
A

A
(i) show that the system (4) (or equivalently (5)) has a unique This is equivalent to showing that the homogeneous system
TR

TR
solution a∗j , j = 0, 1, · · · , m;
m
.
(ii) show that this solution indeed minimizes S. aj sj +k = 0, 0 ≤ k ≤ m (6)
j =0

Once these are proven, then we can conclude that

EX

EX
has only the trivial solution, i.e., aj = 0, j = 0, 1, · · · , m.
There are altogether (m + 1) unknowns aj , j = 0, 1, · · · , m
in (m + 1) normal equations, solving the equations for aj ’s Suppose that NOT all aj ’s are zero, say ak ̸= 0, for some
will give the least-square polynomial pm (see (1)). k ∈ {0, 1, · · · , m}. Multiply the kth equation of (6) by ak
PJY Wong 107 PJY Wong 108
and sum all equations, to obtain Equation (7) implies that

A
⎛ ⎞ 4 n 5 m
.
m m m .
m
. . . . j +k akxki = 0 = pm(xi)

TR

TR
⎝ a k a j sj + k ⎠ = ak aj xi k=0
k=0 j =0 k=0 j =0 i=0
4 5⎛ m ⎞
n
. m
. . for all 0 ≤ i ≤ n, i.e., the polynomial pm has (n + 1) zeros.
= akxk ⎝
i aj xji ⎠ Since n > m, this is impossible! Hence, the homogeneous
i=0 k=0 j =0
system has only the trivial solution.
EX

EX
n
4 m
52
. .
= akxki =0 (7)
i=0 k=0

PJY Wong 109 PJY Wong 110

(ii) To show that the unique solution of system (4) Then, from (2)
A

A
⎡ ⎤2
indeed minimizes S n
!
⎣
m
! ∗ j
m
! ∗ j
S = aj xi − yi + ( a j − aj ) x i ⎦

a∗j ,
TR

TR
i=0
Let 0 ≤ j ≤ m be the unique solution of system (4). ⎧j⎛
=0 j =0
⎞2 ⎡ ⎤2
!n ⎪⎨ ! m !m
⎝ ∗ j ∗ j
= a j x i − yi ⎠ + ⎣ ( aj − aj ) x i ⎦
⎪
i=0 ⎩ j =0 j =0
⎛ ⎞⎡ ⎤⎫
!m !m ⎬
∗ j ∗ j
+ 2⎝ a j x i − yi ⎠ ⎣ ( aj − a j ) x i ⎦
⎭
j =0 j =0
⎛ ⎞2 ⎛ ⎞1 2
EX

EX
!n m
! !n !m !m
⎝ ∗ j ⎠ ⎝ ∗ j ⎠ ∗ k
≥ aj xi − yi +2 aj xi − yi ( ak − a k ) x i
i=0 j =0 i=0 j =0 k=0
⎛ ⎞2 ⎛ ⎞
n
! m
! m
! n
! m
!
⎝ ∗ j ∗ ∗ j k
= aj x i − yi ⎠ + 2 (ak − ak ) ⎝ aj x i − yi ⎠ x i .
i=0 j =0 k=0 i=0 j =0
(8)
PJY Wong 111 PJY Wong 112
 Note that the second term in (8) is zero because of (3). Summary

A
Hence, (8) gives
1. Given (n + 1) points (xi, yi), i = 0, 1, · · · , n. The least-
⎛ ⎞2

TR

TR
n
. m
. square polynomial pm(x) of degree at most m (where
S≥ ⎝ a∗j xji − yi⎠ , m < n) is (1)
i=0 j =0

pm(x) = amxm + am−1xm−1 + · · · + a1x + a0.

i.e., a∗j , 0 ≤ j ≤ m minimizes S.
EX

EX
2. The coefficients a0, a1, · · · , am are obtained by solving (m+
1) normal equations (5)

PJY Wong 113 PJY Wong 114

⎧ Example 2.1. Use the method of least squares, fit a straight

A
⎪
⎪ s0 a 0 + s1 a 1 + · · · + sm a m = v 0
⎪
⎪ line to the following points:
⎪
⎨ s a + s a + ··· + s
1 0 2 1 m+1am = v1 (−1.0, 1.000), (−0.1, 1.099), (0.2, 0.808), (1.0, 1.000).

TR
⎪
⎪ ·········
⎪
⎪
⎪
⎩ s a +s Solution Here m = 1. The least-square straight line is
m 0 m+1a1 + · · · + s2mam = vm.
Here, y = a0 + a1x.
n
. n
.
sk = xki, vk = yixki.
The normal equations are

EX
i=0 i=0
⎧
(s0 = the no. of points given) ⎨ s a +s a =v
0 0 1 1 0
⎩ s1 a 0 + s2 a 1 = v 1

PJY Wong 115 PJY Wong 116

We find the least-square straight line is

A
. .
s0 = 4, s1 = xi = 0.1, s2 = x2i = 2.05, y = a0 + a1x = 0.9773 − 0.0224x.

TR

TR
. .
v0 = yi = 3.907, v1 = xiyi = 0.0517.
Example 2.2. Fit a parabola (degree 2 polynomial) to the data
Hence, the normal equations become
x 1 2 3 4 5 6 7 8 9
⎧
⎨ 4a + 0.10a = 3.9070
EX

EX
0 1 y 2 6 7 8 10 11 11 10 9
⎩ 0.1a0 + 2.05a1 = 0.0517

from which we obtain a0 = 0.9773, a1 = −0.0224. Hence,

PJY Wong 117 PJY Wong 118

. .
Solution Usual method: Here, m = 2. The least-square sk = xki, vk = yixki.
A

A
parabola is i i

y = a 0 + a 1 x + a 2 x2 . Find the coefficients sk and vk of the system (9), then

TR

TR
substitute into (9) and solve for a0, a1, a2.
The normal equations are
(Do this as an exercise)
⎧
⎪
⎪
⎨ s0 a 0 + s 1 a 1 + s 2 a 2 = v 0
⎪
s1 a 0 + s2 a 1 + s3 a 2 = v 1 (9)
EX

EX
⎪
⎪
⎪
⎩ s2 a 0 + s3 a 1 + s4 a 2 = v 2

where
PJY Wong 119 PJY Wong 120
Substitution method: We use the substitution
u w uw u2 u2 w u3 u4

A
u = x − 5, w = y − 11.
−4 −9
The least-square parabola is

TR

TR
−3 −5
−2 −4
w = a 0 + a 1 u + a 2 u2 .
−1 −3
The normal equations are as in (9), where 0 −1
. .
ui k , w i ui k .
EX

EX
sk = vk = 1 0
i i
2 0
We shall calculate the coefficients sk and vk before solving
the system (9). 3 −1
4 −2
PJY Wong 121 PJY Wong 122

0 −25 51 60 −249 0 708

Thus, Hence, the required parabola is

A

A
s0 = 9, s1 = 0, v0 = −25, v1 = 51,
s2 = 60, v2 = −249, s3 = 0, s4 = 708. w = a 0 + a 1 u + a 2 u2
TR

TR
Substituting these values into the normal equations (9) gives
⎧ ⇒ y − 11 = a0 + a1(x − 5) + a2(x − 5)2
⎪
⎪
⎨ 9a0 + 0 + 60a2 = −25
⎪
0 + 60a1 + 0 = 51 ⇒ y = −1 + 3.55x − 0.27x2.
⎪
⎪
⎪
EX

EX
⎩ 60a0 + 0 + 708a2 = −249

which we solve to get

a0 = − 230
231
51
, a1 = − 60 , a2 = − 247
924
.
PJY Wong 123 PJY Wong 124
 E. Numerical Methods for Differential Equations • The order of an ODE is the order of the highest derivative
1. Euler Method in it.

1.1. Basic Terminology • The general solution of a nth order ODE contains n
arbitrary constants. These constants can be determined
We only deal with ordinary differential equations (ODE).
if n initial/boundary conditions are given.
• An ODE is a differential equation in which an unknown
function y (dependent variable) depends on 1 independent • Initial conditions (IC) – conditions specified at the same
variable x. For example value of x. For example

y ′ = cos x, y ′′ + 4y = 0, y (3) + 2y ′′ + y = x. y(0) = 1, y ′(0) = 2; y(2) = 5, y ′(2) = 0.3, y ′′(2) = 1.

PJY Wong ① ② 125 PJY Wong 126

• boundary conditions (BC) – conditions specified at the satisfy BC. For example
boundary points of a solution interval. For example
y ′′ − y = 0, y(0) = 5, y ′(1) = −3.
y(0) = 1, y(π) = 2; y(0) = 1, y ′(π) = 3.

1.2. Initial Value Problems

• Initial value problem (IVP) – solution of ODE has to satisfy
IC. For example It is sufficient to focus on solving first order ordinary
differential equations (ODE), since every mth order equation
y ′′ − y = 0, y(0) = 5, y ′(0) = 1. can be written as a system of m first order equations.

• Boundary value problem (BVP) – solution of ODE has to

PJY Wong 127 PJY Wong 128
 Consider the first order initial value problem (IVP) The theoretical solution at xn is y(xn), and the
approximate solution at xn will be denoted by yn,
y ′ = f (x, y), y(x0) = y0. (1)
y(xn) ≈ yn.
To solve (1) numerically, we discretize the problem by taking
the grid points (or mesh points) x0, x1, x2, · · · with step
size h, 1.3. Euler Method

xn = x0 + nh, n = 0, 1, 2, · · · . Let y(x) be the exact solution of (1). Expand y(x + h)

using Taylor series
We shall seek the solution y at x0, x1, x2, · · · .

PJY Wong 129 PJY Wong 130

This gives the following recursive formula, known as the

2 Euler method,
h
y(x + h) = y(x) + hy ′(x) + y ′′(x) + · · · .
2!
yn+1 = yn + hf (xn, yn), n = 0, 1, · · · . (2)
In subscript form,

h2 Euler method is of order 1.

y(xn+1) = y(xn) + hy ′(xn) + y ′′(xn) + · · · .
2!
This method is the easiest of all the methods we discuss.
If h is sufficiently small, we can do a first order approximation Geometrically, Euler method amounts to approximating the
y(xn+1) ≈ y(xn) + hy ′(xn) curve of y(x) by a polygon whose first side is tangent to the
curve at x0.
= y(xn) + hf (xn, yn).

PJY Wong 131 PJY Wong 132

 Example 1.1. Use the Euler method with step size h = 0.1 to
solve the IVP at x = 0.3.
y ′ = sin xy, y(0) = 1

Solution Here, f (x, y) = sin xy, x0 = 0, y0 = 1. We have

yn+1 = yn + h sin xnyn = yn + 0.1 sin xnyn, n = 0, 1, · · ·
which gives
y(0.1) ≈ y1 = y0 + 0.1 sin x0y0 = 1 + 0.1(0) = 1
y(0.2) ≈ y2 = y1 + 0.1 sin x1y1 = 1 + 0.1 sin 0.1 =
1.009983342
PJY Wong 133 PJY Wong 134

y(0.3) ≈ y3 = y2 + 0.1 sin x2y2 = 1.030045922. n xn yn exact values error

Example 1.2. Apply the Euler method to the following IVP (8 d.p.) y(xn) y(xn) − yn
0 0.0 0.000 0.000 0.000
y ′ = x + y, y(0) = 0.
1 0.2 0.000 0.021 0.021
Choose h = 0.2 and compute y1, y2, · · · , y5. 2 0.4 0.040 0.092 0.052
3 0.6 0.128 0.222 0.094
Solution Noting that f (x, y) = x + y, x0 = y0 = 0, we have
4 0.8 0.274 0.426 0.152
yn+1 = yn + 0.2(xn + yn), n = 0, 1, · · · . 5 1.0 0.489 0.718 0.229

The computation is recorded in the following table.

PJY Wong 135 PJY Wong 136
In fact, using integrating factor, the IVP can be solved to get error ϵln+1 at xn+1. We have

y(x) = ex − x − 1, h2 h3
ϵln+1 = y ′′(xn) + y ′′′(xn) + · · ·
2! 3!
from which we obtain the values of the exact solution y(xn).
h2
Notice that the error increases as the computational process ≈ y ′′(xn) if h is sufficiently small
2!
is continued to higher values of x.

1.4. Error Analysis for Euler Method h2

≈ f ′(xn, yn) (3a)
2!
The truncation of the Taylor series for the step from xn to or
xn+1 causes a per step truncation error or local truncation ϵln+1 = O(h2). (3b)
PJY Wong 137 PJY Wong 138

The total truncation error or global truncation error ϵgn for 2. Improved Euler Method (Heun’s Method)
n steps is
.n
Idea behind Improved Euler Method
ϵg n = ϵlk+1 (4a)
k=1 From the IVP (1): y ′ = f (x, y), y(x0) = y0, it is clear that
-

= n × O(h ) 2 6 x1 6 x1
′
y (x)dx = f (x, y(x))dx
x0 x0
xn − x 0 2 6
= × O(h ) x1
h ⇒ y(x1) − y(x0) = f (x, y(x))dx
= O(h). (4b) To x0
6 x1
Hence, Euler method is of first order. ⇒ y(x1) = y0 + f (x, y(x))dx
x0

PJY Wong 139 PJY Wong

~
this integral x
140
Now a numerical integration method can be used to
7x
approximate the integral x01 f (x, y(x))dx. -

↓
PJY Wong 141 PJY Wong 142

Let us use trapezoidal rule to approximate the integral Using Euler method gives
7 x1
x0
f (x, y(x))dx. Then,
6 x1 y(x1) ≈ y0 + hf (x0, y0).
x 1 − x0
f (x, y(x))dx ≈ [f (x0, y(x0)) + f (x1, y(x1))] -

x0 2
h It follows that
= [f (x0, y0) + f (x1, y(x1))].
2
6 x1
h
Notice that the expression above involves the yet to be f (x, y(x))dx ≈ [f (x0, y0) + f (x1, y(x1))]
x0 2 -

determined y(x1). Hence, we need to use an approximation

of y(x1). h
≈ [f (x0, y0) + f (x1, y0 + hf (x0, y0))]
2 -

PJY Wong 143 PJY Wong 144

and therefore In each step of the method, we compute first the auxiliary
6 x1 value
y(x1) = y0 + f (x, y(x))dx ∗
yn
x0 +1 = yn + hf (xn, yn) (5)
and then the new value
h
≈ y0 + [f (x0, y0) + f (x1, y0 + hf (x0, y0))]. 1
2 yn+1 = yn + ∗
h[f (xn, yn) + f (xn+1, yn +1)]. (6)
2
In general, we have the improved Euler method
The improved Euler method is a predictor-corrector
h
yn+1 = yn + [f (xn, yn) + f (xn+1, yn + hf (xn, yn))]. method, because in each step we first predict a value by
2
(5) and then correct it by (6). The improved Euler method is
a second order method.
PJY Wong 145 PJY Wong 146

Y f (x y]
=
, Y (3) =
Yo
For simplicity, we denote Example 2.1. Apply the improved Euler method to Example
1.2, choosing again h = 0.2.
a = hf (xn, yn), f (32 y) y ′ = x + y,
, = y(0) = 0. %= 0 CLo
0 =

(7)
∗ Solution We have n (Xn + Yn)
b = hf (xn+1, yn +1) = hf (xn+1, yn + a)
a = hf (xn, yn) = 0.2(xn + yn),
so that the improved Euler method can be written as b = hf (xn+1, yn + a) = 0.2(xn+1 + yn + a)
↓ (Sn+ + Ynta)
1 so that
yn+1 = yn + (a + b). (8)
2 1
yn+1 = yn + (a + b), n = 0, 1, · · · .
2 -

PJY Wong 147 PJY Wong 148

 2(xn + Yn) b = 0 .
2(xm + Yn +a)
a= 0 .

Yntl
=
Yn + = Catb) ,
n = 0, 1
,
2 ...
x0 =
Yo = 0

Step D Step 2 Step 3

2.1. Error Analysis for Improved Euler Method
How ???
n xn yn (8 d.p.) a (8 d.p.) b (8 d.p.) error
0 0.0 0.0000 0 · ·+· yo)
( <Lo 0 2 .
(PatYota)
··· 0.0000 Using Taylor series, the local truncation error at xn+1 can

-
~Y+Yeta)
, = YotE (a+b)
Y

1 0.2 0.0200 1) 0.0014 be shown to be

Yz = Y +
,
=(a+b) 3
2 0.4 0.0884 0.0034 ϵln+1 ≈ − h12 f ′′(xn, yn) if h is sufficiently small. (9)
3 0.6 0.2158 0.0063 Hence,
ϵln+1 = O(h3).
4 0.8 0.4153 0.0102
The global truncation error for n steps is thus (see (4))
5 1.0 0.7027 0.0156
Note that now the error in each step is greatly reduced, ϵgn = O(h2). (10)
compared with that in Example 1.2. The improved Euler method is therefore of second order.

PJY Wong 149 PJY Wong 150

3. Runge-Kutta Method The four auxiliary quantities are defined as

A still more accurate method of great practical importance a = hf (xn, yn)

" #
is the Runge-Kutta method. In particular we would b = hf xn + 12 h, yn + 12 a
" # (11)
examine the fourth order Runge-Kutta method. The local c = hf xn + 12 h, yn + 12 b
truncation error can be shown to be O(h5) and the global d = hf (xn + h, yn + c)
truncation error is O(h4).
and the value of yn+1 is
In this method, four auxiliary quantities a, b, c, d are first
1
computed before getting the new value yn+1. yn+1 = yn + (a + 2b + 2c + d). (12)
6

PJY Wong 151 PJY Wong 152

Example 3.1. Apply the fourth order Runge-Kutta method to and
the IVP yn+1 = yn + 16 (a + 2b + 2c + d).
′
& (hY) y = x + y,
=
y(0) = 0.
The table of computation is as follows: yn, a, b, c, d are in 8 d.p.

Choose h = 0.2 as before and compute five steps.

n xn yn a b c d error
Solution We have f (x, y) = x + y. Hence, 0 0 0 >
-
-
-
0 0.02 0.022
+
0.0444 0
1 0.2 0.021400 0.04428 0.068708 0.071151 0.098510 3E − 16
a = hf (xn, yn) = 0.2(xn + yn) o
2 0.4 0.091818 0.098364 0.128200 0.131184 0.164600 7E − 16
" # O o
b = hf xn + 12 h, yn + 12 a = 0.2(xn + 0.1 + yn + 0.5a) 3 0.6 0.222107 0.164451 0.200863 0.204507 0.245322 11E − 16
" # 4 0.8 0.425521 0.245104 0.289614 0.294065 0.343917 20E − 16
c = hf xn + 12 h, yn + 12 b = 0.2(xn + 0.1 + yn + 0.5b)
5 1.0 0.718251 31E − 16
d = hf (xn + h, yn + c) = 0.2(xn + 0.2 + yn + c)

PJY Wong 153 PJY Wong 154

Step 1 n = 0, x0 = 0, y0 = 0. Step 2 n = 1, x1 = 0.2, y1 = 0.0214

a = 0.2(x0 + y0) = 0 a = 0.2(x1 + y1) = 0.04428

b = 0.2(x0 + 0.1 + y0 + 0.5a) b = 0.2(x1 + 0.1 + y1 + 0.5a) = 0.068708
= 0.2(0 + 0.1 + 0 + 0.5 × 0) = 0.02 c = 0.2(x1 + 0.1 + y1 + 0.5b) = 0.071151
c = 0.2(x0 + 0.1 + y0 + 0.5b) d = 0.2(x1 + 0.2 + y1 + c) = 0.09851
= 0.2(0 + 0.1 + 0 + 0.5 × 0.02) = 0.022 y2 = y1 + 16 (a + 2b + 2c + d) = 0.091818.
d = 0.2(x0 + 0.2 + y0 + c)
= 0.2(0 + 0.2 + 0 + 0.022) = 0.0444
y1 = y0 + 16 (a + 2b + 2c + d) = 0.0214.

PJY Wong 155 PJY Wong 156

 Y' = f(x y)Y(x0)
,
=
Yo

Example 3.2. Apply the fourth order Runge-Kutta method to Solution We have f (x, y) = (y − x − 1)2 + 2. Hence,
the IVP

y ′ = (y − x − 1)2 + 2, y(0) = 1 a = hf (xn, yn) = 0.2[(yn − xn − 1)2 + 2]

" #
b = hf xn + 12 h, yn + 12 a
for 0 ≤ x ≤ 0.4 with h = 0.2.
= 0.2[(yn + 0.5a − xn − 0.5h − 1)2 + 2]
" #
f (x , x) = <Y -
x -

1)" + 2 c = hf xn + 12 h, yn + 12 b

0
= 0.2[(yn + 0.5b − xn − 0.5h − 1)2 + 2]
2Lo =

Yo = d = hf (xn + h, yn + c)
= 0.2[(yn + c − xn − h − 1)2 + 2]

PJY Wong 157 PJY Wong 158

and 4. Systems & Higher Order Equations

1
yn+1 = yn + (a + 2b + 2c + d).
6 The various methods discussed so far apply to ODE of first
order. These methods can be extended to deal with ODE
n xn yn a b c d
of higher order. This is so because an mth order differential
0 0 1 ? ? ? ? equation can always be written as a system of m first order
1 0.2 1.402707341 ? ? ? ? equations. We would start off the discussion here with the
2 0.4 1.822788917 numerical solution of a system of first order ODE.

As an exercise, fill in all the ‘?’ above.

PJY Wong 159 PJY Wong 160

4.1. System of First Order ODE
We shall denote

ũ = (u1, u2, · · · , um).

Consider the first order system

ỹ ′ = f˜(x, ỹ), ỹ(x0) = ỹ0. (13)

PJY Wong 161 PJY Wong 161

~
↓

PJY Wong 161 PJY Wong 16

 Note that system (13) actually means The various methods discussed earlier like the Euler,
⎧ the improved Euler and the Runge-Kutta methods are all
⎪ ′
⎪ y1 = f1(x, y1, y2, · · · , ym)
⎪
⎪ applicable to system (13). These methods for single first
⎪
⎪ ′
⎪
⎨ y2 = f2(x, y1, y2, · · · , ym)
⎪ order equations can be extended to system (13) by writing
⎪
······ vector functions ỹ, f˜, ... instead of scalar functions y, f, ...,
⎪
⎪
⎪
⎪ ym′
= fm(x, y1, y2, · · · , ym) whereas x remains a scalar variable, i.e.,
⎪
⎪
⎪
⎩ y (x ) = y , y (x ) = y , · · · , y (x ) = y
1 0 1,0 2 0 2,0 m 0 m,0. y → ỹ, f → f˜ All he
(14) a → ã, b → b̃ Vector
A solution y˜∗ of the system (14) is c → c̃, d → d˜
x → x, h→h
y˜∗ = (y1∗, y2∗, · · · , ym
∗
).
PJY Wong 162 PJY Wong 163

4.2. Euler Method for Systems Note that

Recall that the Euler method is ỹn = (y1,n, y2,n, · · · , ym,n)

ỹn+1 = (y1,n+1, y2,n+1, · · · , ym,n+1)
yn+1 = yn + hf (xn, yn). (15)

Generalizing the method to system (13), we have

-v W

ỹn+1 = ỹn + hf˜(xn, ỹn).

~ 2 change to
vecfur

(16)
Notation of
• - which dependent variable
• - at which grid point
y•,•

exceptx & h e.g. y1,n ≈ y1(xn), y2,n+1 ≈ y2(xn+1).

Y,
dependent variable

at
cn

PJY Wong 164 PJY Wong Y, 165

n approximate , at In
 The Euler method (16): ỹn+1 = ỹn + hf˜(xn, ỹn) in 4.3. Improved Euler Method for Systems
component form is
Recall that the improved Euler method is
⎧
⎪
⎪ y1,n+1 = y1,n + hf1(xn, y1,n, y2,n, · · · , ym,n) 1
⎪
⎪
⎪
⎪ yn+1 = yn + (a + b) (18)
⎪ y2,n+1 = y2,n + hf2(xn, y1,n, y2,n, · · · , ym,n)
⎪ 2
⎨
······ where
⎪
⎪
⎪
⎪ a = hf (xn, yn),
⎪
⎪ ym,n+1 = ym,n + hfm(xn, y1,n, y2,n, · · · , ym,n)
⎪
⎪ (19)
⎩ y (x ) = y , y (x ) = y , · · · , y (x ) = y b = hf (xn+1, yn + a).
1 0 1,0 2 0 2,0 m 0 m,0.
(17) Generalizing the method to system (13) gives

PJY Wong 166 PJY Wong 167

&h
except

n
ỹn+1 = ỹn + 12 (ã + b̃), (20)

ã = hf˜(xn, ỹn),
(21)
b̃ = hf˜(xn+1, ỹn + ã).

In component form, we have for i = 1, 2, · · · , m,

yi,n+1 = yi,n + 12 (ai + bi), (22)

ai = hfi(xn, y1,n, y2,n, · · · , ym,n),

bi = hfi(xn+1, y1,n + a1, y2,n + a2, · · · , ym,n + am).
(23) *
PJY Wong 168 PJY Wong 168
4.4. Runge-Kutta Method for Systems Generalizing the method to system (13) gives

Recall that the Runge-Kutta method is 1 ˜

ỹn+1 = ỹn + (ã + 2b̃ + 2c̃ + d) (26)
1 6
yn+1 = yn + (a + 2b + 2c + d) (24)
6

d)
where
where
a = hf (xn, yn),
-except ã = hf˜(xn, ỹn),
" #
" # b̃ = hf˜ xn + 12 h, ỹn + 12 ã ,
8 9 (27)
b = hf xn + 12 h, yn + 12 a ,
" # (25) c̃ = hf˜ xn + 12 h, ỹn + 12 b̃ ,
c = hf xn + 12 h, yn + 12 b ,
d˜ = hf˜(xn + h, ỹn + c̃).
d = hf (xn + h, yn + c).

PJY Wong 169 PJY Wong 170

In component form, we have for i = 1, 2, · · · , m,

1
yi,n+1 = yi,n + (ai + 2bi + 2ci + di) (28)
6
where

ai = hfi(xn, y1,n, y2,n, · · · , ym,n),

" #
bi = hfi xn + 12 h, y1,n + 12 a1, y2,n + 12 a2, · · · , ym,n + 12 am ,
" #
ci = hfi xn + 12 h, y1,n + 12 b1, y2,n + 12 b2, · · · , ym,n + 12 bm ,
di = hfi(xn + h, y1,n + c1, y2,n + c2, · · · , ym,n + cm).
(29)

PJY Wong 171 PJY Wong 17

 *

eq (14)
↓z
J

Example 4.1. Solve the system Use Runge-Kutta method (28), (29). Here, m = 2,
the 32d y1 → y, y2 → z.
dy
dx O
= 2z − 1, y(0) = 0.5; 2 times
argurement From f1(x, y, z) = 2z − 1, we get
&

-
dz ~ a1 = 0.1(2zn − 1), If (Xm Yn En)
= 2y + 1, z(0) = 3.5 , ,

dx hf(xLutzh Yn + Ear Ent >

-

, ,

b1 = 0.1[2(zn + 0.5a2) − 1] = 0.1(2zn + a2 − 1),

to obtain approximate values of y(0.2) and z(0.2) to 3 decimal
places. Use h = 0.1. c1 = 0.1[2(zn + 0.5b2) − 1] = 0.1(2zn + b2 − 1),
d1 = 0.1[2(zn + c2) − 1]. = 0 .
1 (27n + (2
-

1)
Solution We have f1(x, y, z) = 2z − 1, f2(x, y, z) = 2y + 1,

-
x0 = 0, y0 = 0.5, z0 = 3.5.

PJY Wong 172 PJY Wong 173

h step size
,

From f2(x, y, z) = 2y + 1, we get n xn yn zn a1 a2

a2 = 0.1(2yn + 1), hfz (iniyn h[2Yn + ] 0 0 0.5 3.5 0.6 0.2

I
,
En) =

b2 = 0.1[2(yn + 0.5a1) + 1] = 0.1(2yn + a1 + 1), 1 0.1 1.124 3.762 0.6524 0.3248

c2 = 0.1[2(yn + 0.5b1) + 1] = 0.1(2yn + b1 + 1), 2 0.2 1.813 4.154
d2 = 0.1[2(yn + c1) + 1].
b1 b2 c1 c2 d1 d2
Also, 0.62 0.26 0.626 0.262 0.6524 0.3253
1
yn+1 = yn + 6
(a1 + 2b1 + 2c1 + d1), 0.6849 0.3900 0.6914 0.3933 0.7311 0.4631
1 Sub n= 0
zn+1 = zn + 6
(a2 + 2b2 + 2c2 + d2).
Hence, y(0.2) ≈ 1.813 and z(0.2) ≈ 4.154.
The table of computation is below: yn, zn, ai, bi, ci, di are in 8 d.p. likewise for subsequent
Yn & Zr
PJY Wong 174 PJY Wong 175
4.5. Higher Order ODE Set y1 = y, y2 = y ′, y3 = y ′′, · · · , ym = y (m−1).

Given an mth order differential equation, it is first rewritten Then, (30) becomes a system of m first order equations
as a system of m first order differential equations by ⎧
⎪
⎪ y1′ = y2
appropriate substitutions. The resultant system is then solved ⎪
⎪
⎪
⎪
⎪
⎪ y2′ = y3
by any of the methods discussed in Sections 4.2–4.4. ⎪
⎪
⎪
⎨ ······
Consider the mth order IVP ⎪ ′ 1) 1)
y(m
2 +
⎪ ym− 1 = ym
-

y(m
-

⎧ ⎪
⎪ Ym +
= =
You =

⎨ y (m) = f (x, y, y ′, y ′′, · · · , y (m−1)) ⎪

⎪ ′ (m)
⎪
⎪ ym = f (x, y1, y2, · · · , ym) Ym' y(m 1) y f(x y
-
1 +
Y, y(m
&
1)
⎪
-

=
= =
, , . . . .
⎪
⎪ f (x =
Y, Y Ym)
S

⎩ y(x0) = K1, y ′(x0) = K2, · · · , y (m−1)(x0) = Km. ⎩ y1(x0) = K1, y2(x0) = K2, , · · · , ym(x0) = Km. ,
,
. . .,

(30) (31)
PJY Wong 176 PJY Wong 177

We can solve (31) by the methods discussed in previous Solution First convert the second order ODE into a system of
sections. 2 first order ODE. Let let Y Y Y =

,
=

Y = Y Yz
Example 4.2. Consider the IVP
=

- y1 = y, y2 = y ′ .
Yi =
y" =
-

y = -

Y,
′′ ′
y = −y, y(0) = 0, y (0) = 1. (32) Then, (32) is equivalent to Y 10) =
,
0
,
"(0) =
1
11
- -
-
Y Ll

⎧ 1 o

Ye , o
Using -
h = 0.1, find an approximate value of y(0.2) to 3 ⎪
⎪ ′
- - ⎨ y 1 = y2
⎪
decimal places. y2′ = −y1 (33)
-
⎪
⎪
⎪
⎩ y1(0) = 0, y2(0) = 1.

PJY Wong 178 PJY Wong 179

 E
8dp 8dp
Here, n xn yn = y1,n ′
yn = y2,n
f1(x, y1, y2) = y2, f2(x, y1, y2) = −y1, 0 0 0 Y, 0
1 Y2 , 0

E
(
x0 = 0, y1,0 = 0, y2,0 = 1. 1 0.1 0.1 .
Y
1
=
Y 0thY2
,
,
0

1 Ye 1
=
Ye 0-h Y,
,
0
,

Method 1 Solve (33) by the Euler method (17). We have 2 0.2 0.2 Y , z
=
Y ,,, thYo,
0.99 Ye ,
2
=
Ye ,
1
he Y
, ,

Hence, we have y(0.2) ≈ 0.2, y ′(0.2) ≈ 0.99.

yn+1 = y1,n+1 = y1,n + hf1(xn, y1,n, y2,n) Y1 ,
2 Y 2 , 2

= y1,n + hy2,n, Method 2 Solve (33) by the Runge-Kuttar method (28), (29).
Here,
′
yn +1 = y2,n+1 = y2,n + hf2(xn, y1,n, y2,n) yn+1 = y1,n+1 = y1,n + 16 (a1 + 2b1 + 2c1 + d1),
= y2,n − hy1,n. 4 ′
yn 1
+1 = y2,n+1 = y2,n + 6 (a2 + 2b2 + 2c2 + d2),

PJY Wong 180 PJY Wong 181

same

where n xn yn = y1,n ′
yn = y2,n a1 a2
0 0 0 1 0.1 0
f1(x, y1, y2) = y2, f2(x, y1, y2) = −y1, 1 0.1 0.09985 0.9950 0.0995 -0.009985
a1 = 0.1y2,n, a2 = −0.1y1,n, 2 0.2 0.1987
b1 = 0.1(y2,n + 0.5a2), b2 = −0.1(y1,n + 0.5a1), b1 b2 c1 c2 d1 d2
c1 = 0.1(y2,n + 0.5b2), c2 = −0.1(y1,n + 0.5b1), 0.1 -0.005 0.0998 -0.005 0.0995 -0.00998
d1 = 0.1(y2,n + c2), d2 = −0.1(y1,n + c1) 0.09900 -0.01496 0.09875 -0.01494 0.09801 -0.01986
Hence, we have y(0.2) ≈ 0.199.
The table of computation is shown below.

PJY Wong 182 PJY Wong 183

Method 3 Solve (33) by the improved Euler method (20),
(21).

Do this as an exercise.

PJY Wong 184