Course Title : Numerical Methods

Course Code : CUMT213

Lecturer : MR MAKIWA. S. (0777841246)

Department : Mathematics and Statistics

: Chinhoyi University of Technology

: Private Bag 7724, Chinhoyi

Course Outline

Purpose of the course

• To introduce the standard methods and algorithms for the numerical solution of algebraic
equations, numerical linear algebra as well as deferential equations.

Aims of the course

• The course aims at numerically solving all kinds of mathematical problems which arise
from practical applications and can be modeled by different mathematical equations or
inequalities, for example, linear or non-linear differential equations and integral
equations.

Objectives (Intended Learning Outcomes)

• To develop a good understanding of the mathematical basis of the various methods used
for the numerical solution of scientific problems.

Course Content
Numerical analysis provides a powerful means of studying various fundamental methods which
can be used to solve real world problem. It is intended that, on completing the course, students
will be able to apply numerical methods to problems which they meet elsewhere in mathematics
and other subjects in their degree studies. The course consists of six parts as follows:

1. Errors in numerical methods

• Absolute and relative errors.

1
 • Basic sources of errors.

2. Roots of non-linear equations

• Stability and convergence

• Bisection method
• Newtons method
• Secant method
• Fixed point method

3. Numerical solution of a system of liner systems

• Ill conditioned systems

• Gaussian elimination with pivoting
• Gaussian elimination with scaled partial pivoting
• L-U decomposition
• Jacobi process
• Gauss-Seidel method

4. Curve fitting, and polynomial Interpolation

• Lagrange’s interpolation polynomial

• Divided difference table
• Newton’s interpolation polynomial
• Hermite interpolation

5. Numerical differentiation and integration

• Two and Three point forward difference formula(s)

• Trapezoidal rule
• Simpson’s rule
• Interpolatory quadrature
• Gaussian quadrature

6. Numerical solutions of ordinary differential equations

• Taylor series method

• Euler’s method
• Second order Taylor method
• Runge-Kutta methods.

2
Assessment
• There will be one final exam, worth 60% of your total assessment for this unit. The
remaining 40% is made up of three in-class tests.

References
• George Nakos and David Joyner. (1998). Linear Algebra with applications. An
international Thomson Publishing Company

• Burden, Richard L., and Faires, Douglas J. (1993). Numerical Analysis. 5th Ed. Boston:
PWS KENT Publishing Company.

• Curtis, Gerald F, and Wheately, Patrick O. (1994). Applied Numerical Analysis 5th Ed.
Addison-wesley Publishing Company, Reading.

• T. Akai (1996). Numerical Methods for Engineers. John Wiley & sons
• Atkinson, k. E. (1989). An introduction to Numerical Analysis. 2nd Ed. John Wiley
and Sons, New York.

• Haggarty, F. B (1989). Fundamentals of mathematical Analysis. Addison-Wesley

Publishing Company, Wokingham, England.

• Schwarz, H.R. (1989). Numerical Analysis. John Wiley and Sons, Zurich.

• J.H. Mathews, (1992). Numerical Methods for Mathematics, Science and Engineering,
Prentice Hall International Editions.

• Cheney and Kincaid (1991). Numerical Mathematics and Computing. Brooks, Cole
Publishing Company.

3
1 Introduction
1.1 What is numerical analysis ?

• It is a way of solving highly complicated mathematical problems on a computer.

• It is also known as a technique widely used by scientists and engineers to solve their
problems.

1.2 Two issues of numerical analysis

1. How to compute.

• This corresponds to the algorithmic aspects.

2. How accurate is it?

• This corresponds to the error analysis aspect.

1.3 Advantages of numerical analysis

• We can obtain numerical answers of the problems that have no “analytical’ solutions.

• The techniques do not require complex approach but only the basic mathematical oper-
ations, that is, addition, subtraction, multiplication and division.

1.4 Disadvantages of numerical analysis

• Solutions are always numerical.

• Results from numerical analysis are approximations.

1.5 Errors in numerical methods

1.6 Absolute and relative errors

Definition 1. A number a∗ is said to be an approximate number to the exact number a, if it

differs slightly from a.

• This approximate number, a∗ is called a minor approximation if it is less than a and a

major approximation if its is greater than a
√
Example: 1.73 is a major approximation of 3 whereas 1.74 is a major approximation of the
same number.

Definition 2. The difference between the exact number a and the approximate number a∗ is
known as the error

4
 • The error is often denoted by ∆a , thus ∆a = a − a∗

• ∆a < 0 if a < a∗ and ∆a > 0 for a > a∗

• The sign of the error may not be known and hence it is important that we use the
absolute error

Definition 3. If a∗ is an approximation to a, then the absolute error is given by ∆a = |a − a∗ |

• a∗ − ∆a is called minor of approximation.

• a∗ + ∆a is called the major approximation.

Definition 4. If a∗ is an approximation to a, then the relative error is given by

true value − approximate value ∆a

δ= =
true value a

Example

The weight of 1 dm3 of water at 0o C is given as p = 999.847 ± 0.001 gf (gf=gram force).

Determine the relative error of the result of weighing the water.

Solution

• From this question relative error gives an indication of how good a measurement is
relative to the size of the item being measured.

• Hence it is calculated as
absolute error 0.001
δ = = × 100%
value of the item measured 999.847

= 1 × 10−4 .

Example

• Let the value of π be 3.14159265358 and its approximation be 3.14. Compute the absolute
and the relative error of such an approximation.
Solution

– Absolute error = |3.14159265358 − 3.14| = 0.0015926535898.

– Relative error
0.0015926535898
= 0.000506957382897.
3.14159265358

5
1.7 Sources of Errors

Source of error can be categorised into five groups

1. Error of the method

• In formulating conditions that govern the actual phenomena using mathematical

statements, one has to give allowances for errors.
• Generally, the tendency is to give conditions that would simplify the problem at
hand, so that a solution can be found, though approximately.

2. Residual error

• Functions involved in mathematical formulae are sometimes expressed in the form

of an infinite sequences or series e.g

x2
ex = 1 + x + + ···
2!

• In practical situations we consider only a limited number of terms of a sequence

stopping at some term which is then taken to be an approximation to the required
solution.
• Clearly, this kind of termination of the process gives rise to an error.

3. Initial error

• Some errors arise due to the use of approximate numerical parameters whose values
are determined either experimentally or through infinite processes.
• When performing computations with machines (calculator, computer, etc) these
errors are carried forward into the final data leading output.

4. Round-off error

• If the decimal system of numeration or some other positional system is used, there
maybe an infinite number of digits to the right of the decimal point, e.g we may
2
have a non-terminating repeating decimal like 3 = 0.666666.
• More often, we round-off to a specified degree of accuracy and this leads to round-off
error.

5. Operational error

• In some cases we errors may arise when one fails to correctly use the machine (calcu-
lator, computer, etc) or when approximate numbers are used during computations
and this leads to operational error.

6
2 Roots of non-linear equations
2.1 Stability

• Stability in numerical analysis refers to the trend of error change in the iterative scheme.

• It is related to the concept of convergence.

• An iterative scheme is said to be stable if the errors decrease as the iteration progresses,
leading to convergence.

• An iterative scheme is said to be unstable if the errors get larger and larger as the
iteration progresses and this leads to divergence.

2.2 Convergence

A numerical model is convergent if and only if a sequence of model solutions with increasingly
refined solution domains approaches a fixed value.

2.3 Consistency

A numerical model is consistent only if this sequence converges to the solution of the continuous
equations which govern the physical phenomenon being modeled.

2.4 Root finding

• Calculating the roots of an equation

f (x) = 0 (1)

is a common problem in applied mathematics.

• In this chapter we will explore on iterative methods used to solve (1).

• Generally, the strategy is to design an iterative process of the form

xn+1 = g(xn ) (2)

with some starting point x0 .

• So that the numerical solution as

xn → x∗ , as n → ∞ (3)

where f (x∗ ) = 0.

7
2.5 Bisection Method

Definition 5. The simplest numerical procedure for finding a root is to repeatedly halve the in-
terval [a, b], keeping the half for which f (x) changes sign. This procedure is called the Bisection
method, and is guaranteed to converge to a root, say α.

• We are looking for a root of a function f (x) which we assume is continuous on the interval
[a, b].

• If f (x) satisfies f (a)f (b) < 0. then bisection method can be used to find x∗ .

• The bisection method consists of the following steps:

a+b
1. Define c = 2 .

2. if b − c ≤ , then accept c as the root and stop.

3. If step 2 is not satisfied, then choose a new interval from [a, c] or [c, b].
Return to step 1 until step until the desired result is achieved.

Example

Use bisection method to find the positive root of f (x) = x2 − 1 in the interval 0, 23 .
 

Solution

• Test for the existence of a root, that is, f (a) · f (b) < 0.

Table 1: The results for f (x) = x2 − 1.

n a b c f (a) f(b) f(c)
1 0 0.75 1.5 - - +
2 0.75 1.125 1.5 - + +
3 0.75 0.9375 1.125 - - +
4 0.9375 1.03125 1.125 - + +
5 0.9375 1.03125 0.984375 - + -

• As n → ∞, x∗ → 1 which is the exact solution.

Example 2

The root of ex − 2 = 0 is known to exist in [0, 2]. Find an approximate value of the root within
a tolerance level of  = 0.01 using the bisection method.

8
 Table 2: The results for f (x) = ex − 2.
n a b c f (a) f(b) f(c) Error bound
1 0 2 1 - + + 1
2 0 1 0.5 - + - 0.5
3 0.5 1 0.75 - + + 0.25
4 0.5 0.75 0.625 - + - 0.125
5 0.625 0.75 0.6875 - + + 0.0625
6 0.6875 0.75 0.7186 - + + 0.03125
7 0.6875 0.7186 0.7031 - + + 0.015625
8 0.6875 0.7031 0.6953 - + + 0.0078125

• Observe that the error for the 8th iteration is 0.0078125 ≈ 0.001.

• Therefore the root in this case is 0.6953

2.5.1 Convergence of the Bisection method

• Let a1 = a and b1 = b and [an , bn ] (n ≥ 1) are the successive intervals in the bisection
process.

• Clearly

a1 ≤ a2 ≤ · · · ≤ b1 = b

and

b1 ≥ b2 ≥ · · · ≥ a1 = a

• Observe that the sequence {an } is monotonic increasing and bounded above and the
sequence {bn } is monotonic decreasing and bounded below.

• Hence both sequence converge.

• We also know that every iteration shrinks the length of the interval by a half, i.e.,
1
bn+1 − an+1 = (bn − an ), n ≥ 1.
2

• which means that

b0 − a)
bn − an = (4)
2−1

• Since the sequences are monotone and bounded, they should converge, i.e.,

lim an = lim bn = α
n→∞ n→∞

9
 • This implies that both sequences converge to the same value or point (denoted by α.)

• Thus
b−a
lim bn − lim an = lim = 0. (5)
n→∞ n→∞ n→∞ 2n−1

• since f (an ) · f (bn ) ≤ 0, it follows that [f (α)]2 ≤ 0 and this implies that f (α) = 0.

• Hence, both an and bn converges to a root of f (x) = 0

• Let us apply the bisection method to the interval [an , bn ]

bn + an
• We have cn = ,
2
• Thus α ∈ [an , cn ] or α ∈ [cn , bn ]

• It follows that

1
|α − cn | ≤ cn − an = bn − cn = (bn − an ) (6)
2

• This is the error bound for cn defined in step 2 of bisection method.

• Combining (4) and (6) gives

 
1 b−a 1
|α − cn | ≤ n−1
= n (b − a). (7)
2 2 2

• This shows that the iterates cn → α as n → ∞

• To observe how many iterations will be necessary to achieve the given tolerance level,
that is,

|α − cn | ≤ 

• This will be satisfied if

1
(b − a) ≤ .
2n

Example 3

The root of ex −2 = 0 is known to exist in [0, 2]. Use bisection method to determine the number
of iterations required to achieve an approximate value with a desired tolerance level of 0.01.

10
Solution

1
(b − a) ≤ 
2n

2(1−n) ≤ 0.01

n ≥ 7.64386

• Thus approximately 8 iterations are required to achieve the result within the desired level
of accuracy.

Advantages

• The method is guaranteed to work if f (x) is continuous in [a, b] and a zero actually exists.

• A specific accuracy of iterations is known in advance.

Disadvantages

• It requires an interval

• The convergence of interval having is very slow

• Multiple zeros within the chosen interval can cause problems

2.6 Newton’s method

• Newton’s method is a relatively simple, practical, and widely used root finding method.

• The method is based on the linear approximation of the function.

• We assume that f (x) has at least one (real) root and we denote it by α.

• we start with an initial guess for the location of the root, say x0 .

• We then let l(x) be the tangent line of f (x) at x0 , that is.,

l(x) − f (x0 ) = f 0 (x0 )(x − x0 )

• The intersection of l(x) with the x − axis serves as the next estimate of the root.

• We denote this point by x1 and write

0 − f (x0 ) = f 0 (x0 )(x1 − x0 )

• which means that

f (x0 )
x1 = x0 −
f 0 (x0 )

11
 • Since x1 is expected to be an improvement over x0 as an estimate of α, we repeat the
procedure with x1 as a guess to get:
f (x1 )
x2 = x1 −
f 0 (x1 )

• Repeating this process, we obtain a sequence of iterates, x1 , x2 , x3 , .... hopefully ap-

proaching the root α.

• The iteration formula

f (xn )
xn+1 = xn − (8)
f 0 (xn )
is referred to as the Newton’s method or Newton-Raphson, for solving f (x) = 0.

Example 4

Use Newton-Raphson’s method to approximate the positive root of f (x) = x2 − 1.

Solution

• Observe that

f (xn ) = x2n − 1, f 0 (xn ) = 2xn

• Let x0 = 0.2, f (0.2) = −0.96, f 0 (0.2) = 0.4

• Thus:
[−0.96]
x1 = 0.2 − = 2.6
0.4

• Repeating the process up to x5 gives x5 = 1.000004486540379.

• The Newton’s method may converge slowly at first.

• However, as the iterated progress the speed of convergence increases.

2.6.1 Speed of converge of Newton’s Method

• A Taylor series expansion for f (x + h) is given by

1 1
f (x + h) = f (x) + hf 0 (x) + h2 f 00 (x) + h3 f 000 (x) + · · · .
2 3!

• If α is root of f (x), then f (α) = 0

• We can write α = xn + en , where xn and en denote the estimate and the error after n
iterations

12
 • Thus,

f (α) = f (xn + en )
1 1
= f (xn ) + en f 0 (xn ) + e2n f 00 (x) + e3n f 000 (x) + · · ·
2 3!
• If we assume that the third and higher order terms are negligible we have
1
f (α) = f (xn ) + en f 0 (xn ) + e2n f 00 (cn ) (9)
2
where cn is an unknown point between α and xn .

• Since f (α) = 0, if we divide (9) by f 0 (xn ) we have

f (xn ) 00
2 1 f (cn )
0= + α − x n + (α − x n ) (10)
f 0 (xn ) 2 f 0 (xn )

• But from Netwon’s method we have

f (xn )
− xn = −xn+1 .
f 0 (xn )

• Solving for α − xn+1 , we have

1 f 00 (cn )
α − xn+1 = −(α − xn )2 (11)
2 f 0 (xn )

• We know that the error after n iterations is en = α − xn and the error after (n + 1)
iterations is en+1 = α − xn+1

• Thus
1 f 00 (cn )
en+1 = e2n ⇒ en+1 ≤ |A||en |2 (12)
2 f 0 (xn )

• Equation (12) implies that the error in xn+1 is near proportional to the square of the
error in xn .

• When the initial error is sufficiently small, this shows that the error in the succeeding
iterates will decrease very rapidly.

Advantages of Newton’s method

• Starting point x0 can be arbitrary

• The convergence is fast compared to other root finding methods such as bisection method.

Disadvantages of Newton’s method

• The method requires that f 0 (x) exists

• Divergence may occur

13
2.7 Secant Method

• It is a modified version of Newton-Raphson’s method in which

f (xn ) − f (xn−1 )
f 0 (xn ) =
xn − xn−1

• Thus,
xn − xn−1
xn+1 = xn − f (xn ). (13)
f (xn ) − f (xn−1 )

Example 5

Use secant’s method to solve f (x) = x2 − 1, starting with x0 = 0 and x1 = 0.75

Solution

x1 − x0 4
x2 = x1 − f (x1 ) = .
f (x1 ) − f (x0 ) 3

x3 = 0.96, x4 = 0.994186, x5 = 1.00012, x6 =????????

2.7.1 Converge of Secant method

• Let the error after (n + 1) iteration be:

en+1 = α − xn+1

xn − xn−1
= α − xn + f (xn )
f (xn ) − f (xn−1 )
 
xn − xn−1
= en + f (xn )
f (xn ) − f (xn−1 )
 
(xn − α) + (α − xn−1 )
= en + f (xn )
f (xn ) − f (xn−1 )
 
en − en−1
= en − f (xn )
f (xn ) − f (xn−1 )
en−1 f (xn ) − en f (xn−1 )
=
f (xn ) − f (xn−1 )
  
xn − xn−1 en−1 f (xn ) − en f (xn−1 )
=
f (xn ) − f (xn−1 ) xn − xn−1
  
xn − xn−1 en−1 f (xn ) − en f (xn−1 )
=
f (xn ) − f (xn−1 ) xn − xn−1
 
  f (xn ) − f (xn−1 )
xn − xn−1  en en−1 
= en en−1
f (xn ) − f (xn−1 ) xn − xn−1

14
• Now by Taylor’s Theorem we have (xn = α − en )
1
f (xn ) = f (α) − f 0 (α)en + f 00 (α)e2n − O(e3n )
2
1
= 0 − f 0 (α)en + f 00 (α)e2n − O(e3n )
2

• Recall that xn = α − en , and f (α) = 0.

• Thus,
f (xn ) 1
= −f 0 (α) + en f 00 (α) − O(e2n )
en 2

• Similarly

f (xn−1 ) 1
= −f 0 (α) + en−1 f 00 (α) − O(e2n−1 )
en−1 2

• Therefore
f (xn ) f (xn−1 ) 1
− = f 00 (α)(en − en−1 ) − O(e2n ) − O(e2n−1 )
en en−1 2

• Now
 " 1 00 (α)(e − e
#
xn − xn−1 2 f n n−1 )
en+1 ≈ en en−1 .
f (xn ) − f (xn−1 ) xn − xn−1

• Now

en − en−1 = (α − xn ) − (α − xn−1 ) = xn − xn−1

and for xn and xn−1 sufficiently close to α

xn − xn−1 1
≈ 0 .
f (xn ) − f (xn−1 ) f (α)

• So
1 f 00 (α)
en+1 ≈ en en−1 = Cen en−1 . (14)
2 f 0 (α)

• In order to determine the order of convergence, we now suppose an asymptomatic rela-

tionship of the form

|en+1 | ≈ A|en |r

• This relationship requires

|en | ≈ A|en−1 |r ⇒ |en−1 | = (A−1 |en |)1/r

15
 • In order to be consistent with (14)

|en+1 | = A|en |r = C|en en−1 | = C|en |(A−1 |en |)1/r

• Which reduces to
1 1
A|en |r = CA− r |en |1+ r

• or

1 1
A1+ r C −1 = |en |1−r+ r . (15)

• The left-hand-side of (15) is non-zero constant while the right-hand-side of (15) tends to
zero as n → ∞ (assuming, of course, that the method converges). This is possible only
if
1
0=1−r+ .
r
•
√
2 1 5
r −r−1=0⇒r = ± .
2 2
• Taking the positive root (otherwise, the error terms asymptotically diverge), we find

r ≈ 1.62 < 2.

• Thus, the rate of convergence of the secant method is super-linear, but not quadratic.

Advantages of secant method

• Converges faster compared to bisection method

• The method does not require the derivative

Disadvantage(s)

In some cases the method may fail if we start too far from the solution.

2.8 Fixed point method

• In this method, one writes f (x) = 0 in the form x = g(x) so that any solution of x = g(x)
is a solution of f (x) = 0.

• Hence the method is also known as the x = g(x) method.

• The function g(x) is also called an iteration function.

• Once g(x) is chosen, then we carry out the iterations ( starting from an initial guess x0 )

xn+1 = g(xn ), n = 0, 1, 2, · · ·

16
Example

Given f (x) = x2 −2x−3 = 0. Use the fixed point method to determine the roots of the function.

Solution

Use direct calculations one can easily see that the roots of f (x) are x = −1 and x = 3.

Applying the fixed point method

Arrangement 1: x2 = 2x + 3

√
• Thus xn+1 = 2xn + 3

• Choosing x0 = 4 we have
p √
• x1 = 2(4) + 3 = 11 = 3.31662

• x2 = 3.10375, x3 = 3.03439, x4 = 3.01144

• From the above iterations it appears that we will converge to x = 3, which is another
root of f (x).

3
Arrangement 2: xn+1 =
xn − 2
• Using x0 = 4 we have

• x1 = 1.5 x2 = −6, x3 = −0.375, x4 = −1.263158,

• x5 = −0.919355, x6 = −1.02762, x7 = −0.99087,

• x8 = −1.00305

• It appears that we will converged to x = −1, which is another root of f (x).

• We also observe that the convergence is oscillatory rather than monotonic, as in the first
arrangement.

Example
x2n − 3
Arrangement 3: xn+1 =
2

• Choosing x0 = 4 we have

• x1 = 6.5, x2 = 19.625, x3 = 191.070

• It appears that we are diverging

17
2.8.1 Convergence Theorem

• Consider a function f (x) and suppose it has a zero on the interval [a, b]

• Also. consider the iteration scheme

xn+1 = g(xn )

derived using fixed-point method.

• Then this scheme converges, i.e.,

xn → x∗ , as, n → ∞

if the following condition(s) are satisfied:

1. |g 0 (x)| < K for all x ∈ [a, b], for all 0 < K < 1

2.8.2 Condition for a unique fixed point

Theorem 2.1. Let g(x) be a continuous and differentiable function that maps its domain onto
a subset of itself, i.e., g(x) : [a, b] 7→ S ⊂ [a, b]. Suppose further that there exists some positive
constant 0 < K < 1 such that

|g 0 (x)| ≤ K (16)

for all x ∈ [a, b]. Then g(x) has a unique fixed point p ∈ [a, b].

Proof

• Assume that g(x) is continuous and bounded on the real numbers and has a fixed point.

• Let the fixed point be p so that

g(p) = p. (17)

• Suppose that a second fixed point q 6= p exists.

• Since q is also a fixed point we have

q = g(q). (18)

• By the Mean Value theorem, there exists some number c ∈ [min(p, q), max(p, q)] such
that
g(p) − g(q)
g 0 (c) = (19)
p−q

18
 • By equation (16), |g 0 (c)| ≤ K we have

g(p) − g(q)
≤ K. (20)
p−q

• Thus

|g(p) − g(q)| ≤ K|p − q| < |p − q| (21)

because K < 1.

• But by equations (17) and (18) we have

|g(p) − g(q)| = |p − q| (22)

• and therefore

|p − q| < |p − q| (23)

• Since p 6= q we know that |p − q| =

6 0 hence we can cancel it both sides of the equation
(23) to get 1 < 1, which is a contradiction. hence our original assumption p 6= q must be
wrong. Thus the fixed point is unique.

Example

Show that
1 x
g(x) = π + sin
2 2
has a unique fixed point.

Solution

• Observe that g(x) is continuous and differentiable, and that

 
1 1
Range of g = π − , π + ⊂ (−∞, ∞) = Domain of g.
2 2

• Since g(x) is continuous and differentiable it follows that

1 x 1
|g 0 (x)| = cos ≤ < 1.
4 2 4

• Therefore, we conclude that a unique fixed point exists.

19
2.8.3 Order of convergence of Fixed point method

• Suppose the iteration xn+1 = g(xn ) converges to x∗ where x∗ is a fixed point of the
function g(x).

• Thus, f (x∗ ) = 0

• Consider

xn+1 = g(xn ) Estimate (24)

x∗ = g(x∗ ) Exact (25)

• Thus the error is after (n + 1) iterations given by

xn+1 − x∗ = g(xn ) − g(x∗ ) (26)

• In general, the error at the nth iteration is given by

e n = xn − x∗

• We can write xn as

xn = x∗ + xn − x∗
= x∗ + en

• Thus,

xn+1 − x∗ = g(xn ) − g(x∗ )

= g(x∗ + en ) − g(x∗ )
= en+1

• In general, the error at the (n + 1)th iteration is given by

en+1 = xn+1 − x∗
= g(x∗ + en ) − g(x∗ ) (27)

• Since we assume convergence it follows that

lim xn = x∗ , lim en = 0
n→∞ n→∞

• In order to deduce how fast en converges to zero, we carry out the Taylor series expansion
of g(x∗ + en ) around x∗ , i.e.,
1 00 ∗ 2
g(x∗ + en ) = g(x∗ ) + g 0 (x∗ )en + g (x )en + .... (28)
2!

20
 • Substituting (28) into (27) we have

en+1 = g(x∗ + en ) − g(x∗ )

1
= g 0 (x∗ )en + g 00 (x∗ )e2n + ...... (29)
2!

• For small en , the first non-zero term on the right hand side of (29) is the dominant term
(or the leading term)

• This term determines the order of convergence.

• If g 0 (x∗ ) 6= 0, then from (29) we have

|en+1 | ≈ |g 0 (x∗ )||en | (30)

• Equation (30) is called Linear (or First order) convergence.

• if g 0 (x∗ ) = 0 but g 00 (x∗ ) 6= 0, then we have

1 00 ∗
|en+1 | ≈ g (x ) |en |2
2!

• This is called quadratic convergence.

• If g 0 (x∗ ) = · · · g (p−1) (x∗ ) = 0 but g (p) 6= 0, then we have

g (p) ∗
|en+1 | ≈ (x ) |en |p
p!

2.8.4 Error estimate for fixed point iteration

If fixed point iteration is terminated after n ≥ 1 steps then the error is limited by

K n |p1 − p0 |
|pn − p| ≤ . (31)
1−K
Proof: We will utilize induction.

• For n = 1 we need to prove

K
|p1 − p| ≤ |p1 − p0 | (32)
1−K

• To demonstrate (32) we use the Mean Value Theorem: there is some number c between
p0 and p such that

g(p0 ) − g(p) p1 − p
|g 0 (c)| = = ≤K (33)
p0 − p p0 − p

21
 • Hence by the triangle inequality we have

|p1 − p| ≤ K|p0 − p|
= K|p0 − p1 + p1 − p|
≤ K(|p0 − p1 | + |p1 − p|)
= K|p0 − p1 | + K|p1 − p| (34)

• Solving the last equation (34) for |p1 − p| yields (32). Thus, it holds for n = 1.

• For inductive step we assume that equation (31) holds.

• Now we need to show that it holds for (n + 1) that is

K n+1 |p1 − p0 |
|pn+1 − p| ≤ (35)
1−K

• We again use the Mean Value Theorem: there is some number c between pn and p such
that
g(pn ) − g(p) pn+1 − p
|g 0 (c)| = = ≤K (36)
pn − p pn − p

• Hence

|pn+1 − p| ≤ K|pn − p| (37)

• Substituting equation (31) on the right of (37) gives (38), that is,

K n+1 |p1 − p0 |
|pn+1 − p| ≤ . (38)
1−K

Example

Estimate the number of iterations required for fixed point iteration to converge to he fixed
point of
1 x
g(x) = π + sin
2 2
with 4 digit accuracy i.e 10−4 , using p0 = π.

Solution

To achieve an accuracy of , it is sufficient to find the smallest n such that

K n |p1 − p0 |
|pn − p| ≤ < (39)
1−K

22
Example

From (39) it follows that

 
1 (1 − K)
n> log
log K |p1 − p0 |

• To find K we calculate

0 1 x 1
|g (x)| = cos ≤
4 2 4

• so we choose K = 14 . Recall that p0 = π so that

(0.75) × 10−4
 
1
n> log
log 0.25 |p1 − π|

• To get p1 we iterate once, that is

1 π  1
p1 = π + sin =π+
2 2 2

• Thus
(0.75) × 10−4
 
1
n> log ≈ 6.3.
log 0.25 0.5

• Hence we will need at most 7 iterations.

23
3 Numerical methods for solving linear systems
Numerical methods for solving linear systems of equations can be categorized as
1. Direct methods

2. Iterative methods

3.1 Direct methods

3.1.1 Ill-conditioned systems

Definition 6. A matrix A is said to be ill-conditioned if there exists a vector b for which a small
change in the coefficient matrix A or ~b will produce large changes in the solution x = A−1 b.
Definition 7. A system of equations is said to be well-conditioned if a small change in the
coefficient matrix or a small change in the right hand side results in a small change in the
solution vector.

• System Ax = b is ill-conditioned when the coefficient matrix A is ill-conditioned.

• In this case numerical methods for computing the solution are prone to more errors.

• The Gaussian elimination method usually falters when we have an ill-conditioned system.

Example

Show that the system of equations


 x+y =2


(40)


 x+ 600
601 y = 20.
is ill-conditioned if rounding is done to three significant figures.

Solution

• From the first equation we have x = 2 − y

• substituting this into the second equation gives y = −10818

600
• Now, if we round-off to 3. s.f we have 0.998
601
• System (40) becomes (
x+y =2
(41)
x + 0.998y = 20
• Solving this system gives y = −9000 and x = 9002

• Clearly, the change on solutions after is very big and highly unacceptable. Thus the
system is ill-conditioned.

24
3.1.2 Gaussian elimination method

Definition 8. Row-echelon form: If A is an m × n matrix, then A is said to be in echelon

form if

1. All zero rows are at the bottom of the matrix.

2. The leading entry of each non-zero row after the first occurs to the right of the leading
entry of the previous row.

3. The leading entry in any non-zero row is 1.

4. All the entries in the column below a leading 1 are zero.

• If a matrix satisfies the first two conditions, we say that it is in row echelon for (or simply
echelon form).

• If a matrix satisfies all four conditions, we say that it is in reduced row echelon form (or
just reduced echelon form).

• A reduced echelon form matrix is always in echelon form.

Example

Determine which of the following matrices are in echelon form and reduced echelon form
   
−1 0 1 2 1 −6 0 2
0 2 −6
 
M =  0 0 −3 4  , N =  0 0 1 4  , P =
   
0 −1 1
0 0 0 0 0 0 0 0

Solution

• Matrices M and N are in echelon form

• Matrix N is also in reduced echelon form

• Matrix M is not in reduce echelon form, because condition 3 fails

• Matrix P is not in echelon form, because condition 2 fails.

3.1.3 Weakness of Gaussian elimination method

• Any matrix can be reduced to row echelon form by Gaussian elimination. However, it is
important to note that we may encounter serious computational challenges when using
this method.

• For example if we divide by a number that has been rounded-off the result could carry a
significant round-off error

25
 1 1
• For instance 0.0003 ≈ 3333 while 0.00034 ≈ 2941

• We can clearly observe the discrepancy, for most calculation may be totally unacceptable

• In order to avoid this problem we used a method known as Gaussian elimination with
partial pivoting

• To achieve his we divide by the largest component in a column thereby reducing the
chances of the type of error shown above as much as possible

3.1.4 Gaussian elimination with pivoting

Definition 9. An element app that is used to eliminate xp in rows p + 1, p + 2, ....n is called a

pivot or pivotal element and the row containing this element is called the pivot row

• Generally, in applying the Gaussian elimination method to solve a system of equations if

app = 0, row p cannot be used to eliminate the elements in column p below the diagonal

• It is necessary to find a row k, where akp 6= 0 and k > p and the interchange row p and
row k so that a non-zero pivot element is obtained.

• This process is called pivoting and the criterion for deciding which row to choose is called
pivoting strategy

• A commonly used pivoting strategy is to choose row k which is such that

|akp | = max |aip |, p≤i≤n

i

• Rows k and p are then interchanges and the elimination process is carried out. This
process is called partial pivoting or maximal column pivoting method.

Example

Use Gaussian elimination with partial pivoting to solve

x1 − x2 + x3 = 1
−3x1 + 2x2 − 3x3 = −6
2x1 − 5x2 + 4x3 = 5

Solution

• Step 1: form the augmented matrix, and select the pivot column ( first column with
non-zero components). Rearrange the rows in order to omve the pivot to the top
 
1 −1 1 1
 
 −3 2 −3 −6 
 
2 −5 4 5

26
• Observe that row 2 is our pivot row, thus we interchange row 1 & 2 to have
 
−3 2 −3 −6
 
 1 −1 1 1 
 
2 −5 4 5

• Step 2: Divide the first row by the pivot element and add multiples of the new first row
to make all the other components in he pivot column zero
 
1 − 32 1 2
  R1
−1 1 1   R1 → − 3
 1

2 −5 4 5

• Using the following operations

1
R2 → R2 + R1
3
2
R3 → R3 + R1
3

we have  
1 − 32 1 2
 
 0 −1 0 −1 
 3 
0 − 11
3 2 1
Step 3: Now cover the first row and column, and perform steps 1&2 on the resulting
sub-matrix. Continue with the procedure until the matrix is in two echelon form:

• Observe that our pivot row is the last row, thus

 
1 − 32 1 2
 
 0
 − 11
3 2 1 
0 − 31 0 −1

3
• Using R2 → − 11 R2 we have
 
1 − 23 1 2
 
6 3 
 0 1
 − 11 − 11 
0 − 13 0 −1

1
• Using R3 → R3 − 11 R2 we have
 
1 − 23 1 2
 
6 3 
 0
 1 − 11 − 11 
2
0 0 − 11 − 12
11

27
 Further R3 → − 11
2 R3  
1 − 23 1 2
 
6 3 
 0
 1 − 11 − 11 
0 0 1 6
using back substitution, we have
   
x1 −2
 x2  =  3  .
   

x3 6

3.1.5 Gaussian elimination with scaled partial pivoting

• The most widely used pivoting strategy is called scaled partial pivoting

• In the k th step of elimination process one chooses as pivot equation, an equation from
the n − k + 1 equations which is such that the coefficient of xk and the absolutely largest
coefficient inn the equation is maximum in absolute value

• If there is more than one such equation, take the first of them.

• To find the pivots we use

|ajk |
, k = j = 1, 2, 3, · · · , n
max |ajk |

Example

Apply the Gaussian elimination method with scaled partial pivoting to solve the following
equations

3x1 − 4x2 + 5x3 = −1

−3x1 + 2x2 + x3 = 1
6x1 + 8x2 − x3 = 35

Solution

• To find the first pivot equation we consider the following

|a11 | 3 |a21 | 3 |a31 | 6
= , = , =
max |a1k | 5 max |a2k | 3 max |a3k | 8

• From the above calculations, the second is our pivot equation.

• Thus we interchange rows 1& 2 and eliminate x1

 
−3 2 1 1
 
 3
 −4 5 −1 

6 8 −1 35

28
 • Using the operations

R2 → R2 + R1
R3 → R3 + 2R2

 
−3 2 1 1
 
 0
 −2 6 0 
0 12 1 37

• To find the second pivot we consider

|a22 | 2 |a32 | 12
= , =
max |a2k | 6 max |a3k | 12

• Hence we select R3 as the pivot row

 
−3 2 1 1
 
 0 12 1 37  ,
  R3 → 6R3 + R2
0 −2 6 0

 
−3 2 1 1
 
 0 12 1 37 
 
0 0 37 37

• Using back-substitution we obtain

 
 
x1 2
 x2  =  3  .
   

x3 1

Diagonally dominant matrices

Definition 10. An n×n matrix is said to be diagonally dominant if, in every row, the absolute
value of the diagonal component is greater than the sum of the absolute values of the off-diagonal
components, i.e.
n
X
|ajj | > |ajk |, k 6= j, j = 1, 2, 3....n
k=1

Example

Consider the matrix  

4.1 −2.2 0.8
 
A=
 0.8 2.5 1.1 
−1.6 0.4 −5.2

29
Matrix A is diagonally dominant since

|4.1| > | − 2.2| + |0.8| = 3

|2.5| > |0.8| + |1.1| = 1.9
| − 5.2| > | − 1.6| + |0.4| = 2

Remark: No pivoting is required for diagonally dominant systems.

3.1.6 LU-decomposition

• Given a system, Ax = ~b, we express a non-singular matrix A as a product of a lower

triangular matrix L and an upper triangular matrix U , that is

A = LU (42)

• Thus system (42) can thus be presented as

LU x = b (43)

• We may write (43) as

Ly = b (44)

where y = [y1 , y2 , y3 , · · · , yn ]T and

Ux = y (45)

Theorem 3.1. Let A be a given n × n matrix. denote by Am the m × m matrix formed by the
intersection of the first m rows and m columns of A for m < n. If |Am | =
6 0, m = 1, 2, ...n−1,
the there exists a unique lower triangular matrix L = lij with lii = 1, for i = 1, 2, 3....n and an
upper triangular matrix U = uij such that LU = A.

Proof. • For m = 1, we have a11 = u11 (Induction base) which is unique

• Assume that the decomposition is true for n = m − 1 (induction hypothesis)

• For n = m we partition Am Lm and Um according to

Am−1 b Lm−1 0 Um−1 u
     
Am = , Lm = , Um = (46)
ct Amm pt 1 0 Umm
where c, u and p are column vectors with m − 1 components and b is the right side of
(44).

• We form the product, Am = Lm Um , that is,

Am−1 b Lm−1 0 Um−1 u
    
= (47)
ct Amm pt 1 0 umm

30
 • Equating the elements after doing the multiplication on the right hand side we have

Am−1 = Lm−1 Um−1

b = Lm−1 u
ct = pt Um−1
Amm = pt u + Umm

• By induction hypothesis, Lm−1 and Um−1 are uniquely determined and non-singular.

• It follows that u and p are uniquely determined from the triangular systems:

b = Lm−1 u
t
c = Um−1 p

• Finally, since Umm = Amm −pt u Lm and Um are uniquely determined, and this completes
the proof.

Remark If for some m, we have |Am | = 0, LU − decomposition may not exist.

Example

Solve the following system using the LU − decomposition method

2x + 4y − z = −5
x + y − 3z = −9
4x + y + 2z = 9

Solution

We have to write Ax = LU x = b, thus

    
2 4 −1 1 0 0 u11 u12 u13
 1 1 −3  =  l21 1 0   0 u22 u23 
    

4 1 2  l31 l32 1 0 0 u32 

u11 u12 u13
=  l21 u11 l21 u12 + u22 l21 u13 + u23
 

l31 u11 l32 u12 + l32 u22 l31 u13 + l32 u23 + u33

Equating corresponding entries gives

u11 = 2, u12 = 4, u13 = −1,

l21 u11 = 1 ⇒ l21 = 1/2,
l31 u11 = 4 ⇒ l31 = 2
l21 u12 + u22 = 1 ⇒ u22 = −1

31
 l21 u13 + u23 = 1 ⇒ l32 = 7
l21 u13 + u23 = −3 ⇒ u23 = −5/2
l31 u13 + l32 u23 + u33 = 2 ⇒ u33 = 43/2

Since Ax = LU x = Ly = b we now solve for y = [y1 , y2 , y3 ]T Thus

    
1 0 0 y1 −5
1
 2 1 0   y2  =  −9 
   

2 7 1 y3 9

We have
13 129
y1 = −5, y2 = − , y3 =
2 2
Now, we consider U x = y, that is
    
2 4 −1 x1 −5
 0 −1 − 25   x2  =  − 13
    
2 
43 129
0 0 2 x3 2

Thus

x1 = 1, x2 = −1, x3 = 3

This implies that

   
x 1
 y  =  −1 
   

z 3

3.1.7 Error analysis for linear systems

• If many physical phenomena, the elements of A and b (for Ax = b) are not known exactly,
hence there is need to investigate the corresponding uncertainty in x

• if x∗ is the computed solution, how can one estimate the error x∗ − A−1 b?

• Let x∗ be the computed solution for the linear system Ax = b, then its error is e = x−x∗ ,
which is usually not known.

• The residual error r is defined as

r = b − Ax∗ = A(x − x∗ ) = Ae

• If r = 0 it follows that e = 0 and x∗ is the exact solution.

32
Example

Consider the following system

1.01x1 + 0.99x2 = 2
0.99x1 + 1.01x2 = 2

whose unique solution is x = y = 1.

• Suppose the approximate solution is

 ∗ 
x1 1.01

x∗ = =
x∗2 1.01

• Now
1 1.01 −0.01
     
∗
e=x−x = − =
1 1.01 −0.01

• and
1.01 0.99 −0.01 −0.02
    
r = Ae = =
0.99 1.01 −0.01 −0.2

• In this case, a small residual error corresponds to a small error. However, this is not
always the case, for instance, Let
2
 
x∗ =
0
then
−1 −0.02
   
e= , and, r =
1 0.02
Remarks: Observe that in this case the residual error is still relatively small but the error
is now relatively large. This shows that the size of the residual error r of an approximate
solution x∗ is not a reliable indicator of the size of the error e, in the approximate solution.
In order to solve this problem we shall now discuss the condition number

Conditioning

Definition 11. The condition number of a matrix A is given by

k(A) = A−1 kAk

where k·k is the matrix norm

Proof. Before we embark on the proof we need to consider the properties of matrix norms:

1. For a matrix A, kAk ≥ 0.

33
 2. For a matrix A and a scalar k, kkAk = |k| kAk

3. For two matrices A and B of the same order, kA + Bk ≤ kAk + kBk

4. For two matrices A and B that can be multiplied, kABk ≤ kAk kBk

The following are the most frequently used norms:

P
1. kAk1 = maxj i |aij | = maximum column sum
P
2. kAk∞ = maxj i |aij | = maximum row sum

3. kAk2 = {max eigenvalue of AT A}1/2

Remark: The matrix norm kAk2 that corresponds to the 2-norm of a vector is related to
eigenvalue of the matrix. This norm is special in the sense that no other norm is smaller than
it, and therefore it provides the strictest measure of magnitude of a matrix. It is also called
the spectral norm.

Example

Consider the matrix

 
5 −5 −7
A =  −4 2 −4 
 

−7 −4 5
Find

1. kAk1

2. kAk∞

Solution

1. We have

|5| + | − 4| + | − 7| = 16
| − 5| + |2| + | − 4| = 11
| − 7| + | − 4| + |5| = 16

Thus, kAk1 = 16

2. kAk∞ = 17

• Let the exact problem be Ax = b and the approximate be Ax∗ = b + b∗

• The residual error is

r = b∗ = A(x − x∗ ) (48)

34
 • Taking the norm of both sides we obtain

kx − x∗ k = A−1 b∗ .

• The relative error leads to

kx − x∗ k kbk A−1 b∗ kbk A−1 kb∗ k

= ≤ (49)
kxk kxk kbk kxk kbk

• Using the fact that b = Ax it follows that

kbk = kAxk ≤ kAk kxk

• Thus (49) becomes

kx − x∗ k kb∗ k kb∗ k
≤ A−1 kAk = k(A) . (50)
kxk kbk kbk

• Since b∗ = r, (50) can be presented as

kx − x∗ k krk
≤ k(A) (51)
kxk kbk

which gives a relationship between the residual and the relative error in the solution.

Remarks

1. k(A) ≥ 1.

2. If k(A) >> 1 the system is ill-conditioned.

3. If k(A) = O(1) the system is well conditioned.

4. Small errors in b for ill-conditioned systems result in large errors in x.

Application of condition number

• The relative error in a solution vector norm is ≤ Cond(A)× relative error in the right
hand side vector norm.

• The possible relative error in the solution vector norm is ≤ Cond(A) × machine

• Hence, Cond(A) × machine should give us the number of significant digits, m that
are at least correct in our solution by finding out the largest value of m for which
Cond(A) × machine is less that 0.5 × 10m

35
Example

How many significant digits can I trust in the solution of the following system of equations?

1 2 x 2
    
= .
2 3.999 y 4

Solution

• Observe that
−3999 2000
 
−1
A = .
2000 −1000

• It can easily be verified that

kAk∞ = 5.999, kA−1 k∞ = 5999.4

• Thus k(A) = kAkkA−1 k = 5.999 × 5999.4 = 35990.

• Assuming precision with 23 bits used in the mantissa for real numbers, the machine
tolerance (machine ) is 2−23 = 0.119209 × 10−6 .

• Therefore

k(A) × machine = 35990 × 0.119209 × 10−6 = 0.4290 × 10−2

0.4290 × 10−2 ≤ 0.5 × 10−m

m ≤ 2.067

Remark Thus, two significant digits are at least correct in the solution vector.

3.2 Indirect methods

• In addition to direct methods, we also have iterative methods, where we try to approxi-
mate the solution of a system by using iterations, starting with an initial guess.

• If the successive iterations approach the solution, we say that the iteration converges.
Otherwise, we say that it diverges.

• The procedure ends when two consecutive iterations yield the same answer within a
desired accuracy

• Unlike the direct methods, the number of steps need is not known beforehand.

36
3.2.1 Jacobi iteration

• Jacobi iteration applies to a square systems, that is., systems with as many equations as
unknowns.

• Suppose we have a system with n equations in n unknowns x1 , x2 , · · · , xn such as

5x + y − z = 14
x − 5y + 2z = −9
x − 2y + 10z = −30

Step 1: Solve the ith equations of the system for xi , that is:

 xk+1 = −0.2yk + 0.2zk + 2.8


yk+1 = 0.2xk + 0.4zk + 1.8 (52)

= −0.1xk + 0.2yk − 3.0

 z
k+1

(0) (0) (0)

Step 2: Start with an initial guess, x1 , x2 , · · · , xn for the solution. In the absence of any
(0) (0) (0)
information initialize all the variables at zero, that is, x1 , x2 , · · · , xn = 0.
(k−1) (k−1) (k−1)
Step 3: Substitute the values x1 , x2 , · · · , xn obtained after the (k − 1) iteration
(k) (k) (k)
into the right side of (52). to obtain new values x1 , x2 , · · · , xn .

Step 4: Stop the process when a desired accuracy has been achieved. Usually we stop when
two consecutive iterations yield the same values up to this accuracy.

Solution

Table 3:
n x(n) y (n) z (n)
Initial guess 0.0000 0.0000 0.0000
1 2.8000 1.8000 -3.0000
2 1.8400 1.1600 -2.9200
3 1.9840 1.0000 -2.9520
4 2.0096 1.0160 -2.9984
5 1.9971 1.0026 -2.9978
6 1.9999 1.0003 -2.9992
7 2.0001 1.0003 -2.9999
8 2.0000 1.0001 -3.0000
9 2.0000 1.0000 -3.0000
10 2.0000 1.0000 -3.0000

• The iterations suggest that x = 2, y = 1 and z = −3 is the solution of the system, correct
to at least four decimal places. In fact, this is the exact solution in this case.

37
Failure of the Jacobi process

• Consider the system 

 x + 6y − 2z = 7


x − 3y + 7z = 16 (53)


 5x + y + z = 10

• If we set x0 = y0 = z0 = 0 as our initial guess.

• Applying, the Jacobi method we have


 xk+1 = −6yk + 2zk + 7


1 (54)
yk+1 = 3 [xk + 7zk − 16]

= −5xk − yk + 10

 z
k+1

• Solving using a comp (54) use computer yields

8.50136 × 1037
   
x55
 y55  =  3.0087 × 1037
   

z55 z55 × 4.4372 × 1037

• However if (53) is solved using direct methods the solution is

   
x 1
y  = 2
   

z 3

• Clearly, we can conclude that the Jacobi process has diverged in this case

• The results suggests that we should be able to tell in advance, whether or not the process
is going to converge for any given system.

Theorem 3.2. For a system Ax = b, the Jacobi iterative process will converge to the exact
solution if the coefficient matrix A is strictly diagonally dominant, that is., if for each row we
have
n
X
|ajj | > |ajk |, k 6= j, j = 1, 2, 3....n
k=1

Proof. • Consider the linear system


 a11 x + a12 y + a13 z = b1


a21 x + a22 y + a23 z = b2 (55)


 a x+a y+a z =b
31 32 33 3

• Suppose the true solutions of these equations are x, y, z and the iterative values after n
steps of the Jacobi process are xn , yn , zn

38
• Let


 xn = x + αn


yn = y + βn (56)


 z =z+γ
n n

• Further, let En = max{|αn |, |βn |, |γn |} so that Er is the largest error after n steps.

• If Er gets smaller and smaller as n → ∞, then the system converges, otherwise it diverges

• From (55) he Jacobi iterative schemes are

1
xn+1 = [b1 − a12 yn − a13 zn (57)
a11
• From (56), after (n + 1) iterations we have
1
x + αn+1 = [b1 − a12 yn − a13 zn ] (58)
a11

• Since x is a solution, from (55) we have

1
x= [b1 − a12 y − a13 z] (59)
a11

• Thus

αn+1 = xn+1 − x

a12 a13
= − βn − γn (60)
a11 a11
and it follows that
|a12 | |a13 | |a12 | |a13 |
|αn+1 | = |βn | + |γn | ≤ Er + En (61)
|a11 | |a11 | |a11 | |a11 |

• Thus
 
|a12 | |a13 |
|αn+1 | ≤ + En (62)
|a11 | |a11 |

• Applying the same techniques for the other two equations, we can show that
 
|a21 | |a23 |
|βn+1 | ≤ + En ,
|a22 | |a22 |
 
|a31 | |a32 |
|γn+1 | ≤ + En
|a33 | |a33 |

39
 • Now, if 
 |a11 | > |a12 | + |a13 |


|a22 | > |a21 | + |a23 | (63)

 |a | > |a | + |a |

33 31 32

then

|αn+1 | < En , |βn+1 | < En , |γn+1 | < En

• It follows that

En+1 < En (64)

• Thus, the process will converge if the inequalities (63) are satisfies. Therefore if the
matrix A is strictly diagonally dominant, then the Jacobi iterative process for the set of
equations Ax = b will converge.

3.2.2 Gauss-Seidel Iterative method

• The Gauss-Seidel iterative method also applies to square systems and it applies the
following steps

Step 1: Similar to the Jacobi iterative method.

Step 2: Similar to the Jacobi iterative method.

Step 3: Substitute the most recently calculated unknown into the right side of the equation
(k)
obtained in step 1 to get the new approximation, xi .

Step 4: Same as the Jacobi iterative method.

Example

Consider the system

5x + y − z = 14
x − 5y + 2z = −9
x − 2y + 10z = −30

Apply the Gauss-Seidel method to find the solution.

40
Solution

Step 1: Solve the ith equations of the system for xi , that is:

 xk+1 = −0.2yk + 0.2zk + 2.8


yk+1 = 0.2xk+1 + 0.4zk + 1.8 (65)

= −0.1x − 3.0

 z
k+1 + 0.2y
k+1 k+1

Step 2: Consider initial guess to be x0 = y0 = z0 = 0

Table 4:
n x(n) y (n) z (n)
Initial guess 0.0000 0.0000 0.0000
1 2.8000 2.3600 -2.8080
2 1.7664 1.0301 -2.9706
3 1.9999 1.0117 -2.9976
4 1.9981 1.0006 -2.9997
5 1.9999 1.0001 -3.0000
6 2.0000 1.0000 -3.0000
7 2.0000 1.0000 -3.0000

• Thus, the iterations suggest that x = 2, y = 1 and z = −3 is the solution of the

system, correct to at least four decimal places. In fact, this is the exact solution in
this case.
• From this example we can conclude that the Gauss-Seidel method converges faster
compared to the Jacobi iterative method.

41
4 Interpolation
The word ”interpolation” refers to interpolating some unknown information from a given set of
known information. The technique of interpolation is widely used as a valuable tool in science
and engineering. Mathematical applications of interpolation include derivation of computa-
tional techniques for

• Numerical differentiation

• Numerical integration

• Numerical solutions of differential equations

Problem statement

The basic idea of polynomial interpolation is that we find a polynomial which agrees with the
data from the function f of interest. Given a set of measured data, say n + 1 pairs:

(x0 , y0 ), (x1 , y1 ), · · · , (xn , yn )

the problem of interpolation is to find a function f (x) such that

f (xi ) = yi , i = 0, 1, · · · , n

• xi denotes the nodes

• f (x) is said to interpolate the data and is called an interpolation function

• f (x) is said to approximate g(x) if the data are from a function g(x)

• It is called interpolation (or extrapolation) if f (x) uses values within (or outside) the
interval [x0 , xn ]

4.1 Lagrange Interpolation

• The idea is to construct the interpolating polynomial Pn (x) in the form:

Pn (x) = L0 (x)f0 + L1 (x)f1 + · · · + Ln (x)fn (66)

where each Lk (x) for k = 0, 1, 2.....n, is a polynomial of degree at most n, called the
Lagrange polynomial.

• The Lagrange polynomials {Lk (x)} each of degree n are defined as follows

(x − x1 )(x − x2 ) · · · (x − xn )
L0 (x) = , the factor (x − x0 ) is missing
(x0 − x1 )(x0 − x2 ) · · · (x0 − xn )
(x − x0 )(x − x1 ) · · · (x − xn )
L1 (x) = , the factor (x − x1 ) is missing
(x1 − x0 )(x1 − x2 ) · · · (x1 − xn )

42
 • In general, the Lagrange polynomial is given by

(x − x0 )(x − x1 ) · · · (x − xk−1 )
Lk (x) = , the factor (x − xk−1 ) is missing (67)
(xk − x0 )(xk − xk−1 ) · · · (xk − xk−1 )

where k = 0, 1, 2, .....n.

Fitting two points

Fit the linear polynomial for two given points (x0 , f0 ) and (x1 , f1 )

P1 (x) = L0 (x)f0 + L1 (x)f1

where
x − x1 x − x0
L0 (x) = , L1 (x) = .
x0 − x1 x1 − x0

Fitting three points

Fit the Lagrange polynomial for the given points (x0 , f0 ), (x1 , f1 ) and (x2 , f2 )

P2 (x) = L0 (x)f0 + L1 (x)f1 + L2 (x)f2 (x)

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

= f0 + f1 + f2
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x0 )

Fitting (n + 1) data points

The Lagrange polynomial for fitting (n + 1) data points (x0 , f0 ), (x1 , f1 ), · · · , (xn , fn ) is given
by
n n
X X lk (x)
P − n(x) = Lk (x)fk = fk (68)
lk (xk )
k=0 k=0

where
(x − x0 )(x − x1 ) · · · (x − xn )
lk (x) =
(x − xk )

Example

Use Lagrange’s interpolation polynomial to estimate the value of f (1.2) from the following
data

x 1 1.5 2
f 0.0000 0.4055 0.6931

43
Solution

Since these are three nodes we have

P2 (x) = L0 (x)f0 + L1 (x)f1 + L2 (x)f2 (x)

(x − 1.5)(x − 2) (x − 1)(x − 2) (x − 1)(x − 1.5)

= × 0.0000 + × 0.4055 + × 0.6931
(1 − 1.5)(1 − 2) (1.5 − 1)(1.5 − 2) (2 − 1)(2 − 1.5)

For x = 1.2 we have

(1.2 − 1)(1.2 − 2) (1.2 − 1)(1.2 − 1.5)
P (1.2) = × 0.4055 + × 0.6931
−0.25 0.5

= 0.176348.

Remark: Note that f (x) in this case is f (x) = ln x. hence the interpolation error is

ln(1.2) − 0.176348 = 0.00597.

4.2 Newton’s divided difference method

The following two disadvantages of Lagrangian polynomial method lead to the development of
Newton’s interpolation. They are:

1. it involves a cumbersome arithmetic computations

2. we essentially need to start over the computations if we desire to add or subtract a point
from the data set

Divided difference

• Define the first order divided difference between two nodes (xi , xi+1 ) as

fi+1 − fi
f [1] [xi , xi+1 ] = (69)
xi+1 − xi

• The second order divided difference is given by

f [xi+1 , xi+2 ] − f [xi , xi+1 ]

f [2] [xi , xi+1 , xi+2 ] = (70)
xi+2 − xi

• The third order divided difference is given by

f [xi+1 , xi+2 , xi+3 ] − f [xi , xi+1 , xi+2 ]

f [3] [xi , xi+1 , xi+2 , xi+3 ] = (71)
xi+3 − xi

and the higher oder divided difference as

f [xi+1 , · · · , xi+m ] − f [xi , · · · , xi+m−1 ]

f [m] [xi , xi+1 , xi+2 , · · · , xi+m ] = (72)
xi+m − xi

44
Example

Construct the divided difference table for the following data

x 1 1.5 2
f (x) 0.0000 0.4055 0.6931

Solution

f1 − f0 0.4055 − 0
f [x0 , x1 ] = = = 0.8110.
x1 − x0 1.5 − 1
f2 − f1 0.6931 − 04055
f [x1 , x2 ] = = = 0.5752.
x2 − x1 2 − 1.5
f [x1 , x2 ] − f [x0 , x1 ] 0.5752 − 0.8110
f [x0 , x1 , x2 ] = = = −0.1179.
x2 − x0 2−0
The results can be presented in Table as follows
[1] [2]
i xi fi fi fi
0 1 0.000
0.8110
1 1.5 0.4055 -0.1179
0.57527
2 2 0.6931

Newton’s divided difference interpolation polynomial

Consider the nth degree polynomial written as

Pn (x) = a0 + (x − x0 )a1 + (x − x0 )(x − x1 )a2 + · · · + (x − x0 )(x − x1 ) · · · (x − xn−1 )an (73)

The idea is to find a0 , ..., an so that Pn interpolates the given data:

(x0 , y0 ), (x1 , y1 ), · · · , (xn , yn )

Example

Use Newton’s divided interpolation polynomial to estimate f (1.2) from the following data

x 1 1.5 2
f (x) 0.0000 0.4055 0.6931

45
Solution

• The first step is to construct the divided difference table as done on the previous example.

• Since we have three nodes our polynomial should be of the form

P2 = a0 + (x − x0 )a1 + (x − x0 )(x − x1 )a1

• From the divided difference Table

a0 = f0 = 0.0000,

f1 − f0
a1 = = f [1] = 0.8110,
x1 − x0

a2 = f [2] = −0.1179, (74)

• Therefore the Newton’s divided interpolation polynomial is

P2 (x) = 0.8110(x − 1) − 0.1179(x − 1)(x − 1.5) ≈ −0.1179x2 + 1.10575x − 0.98785.

• Now f (1.2) is given by

P (1.2) = 0.8110(1.2 − 1) − 0.1179(1.2 − 1)(1.2 − 1.5) = 0.169274.

• Therefore the error in this case is ln(1.2) − 0.169274 = 0.0130476.

4.3 Newton’s Interpolation with equally spaced nodes: Forward differences

• The problem of interpolation from tabulated data is considerably simplified if the values
of the function are given at equally spaced intervals of the independent variable.

• If this is the case, then the divided differences can be computed without any divisions
by the associated factors

• Assume that the given set of data is equally spaced such that xi+ − xi = h, then the first
order differences of the functions are defined as:

∆f0 = f1 − f0 , at x0
∆f1 = f2 − f1 , at x1
..
.
∆fi = fi+1 − fi , at xi .

• The second order differences of the functions are given by:

∆2 f0 = ∆(∆f0 ) = ∆f1 − ∆f0 , at x0

46
 ∆2 f1 = ∆(∆f1 ) = ∆f2 − ∆f1 , at x1
..
.
∆2 fi = ∆(∆fi ) = ∆fi+1 − ∆fi , at xi

• The nth order differences of the functions are given by:

∆n fi = ∆n−1 fi+1 − ∆n−1 fi .

• Thus given a set of measured data

(x0 , f0 ), (x1 , f1 ), · · · , (xn , fn )

in which xi+1 − xi = h, then the Newtons’ forward interpolation polynomial is given by

     
s s s
Pn (x) = f0 + ∆f1 + ∆2 f2 + · · · + ∆n fn (75)
1 2 n
where
 
x − x0 s s(s − 1) · · · (s − k + 1)
s= , and, = , k = 1, 2, · · · n
h k k!
and ∆k f0 is the k th order difference of the given data.

Example

Interpolate the following set of data using Newton’s forward difference method

x 0 1 2 3 4
f (x) 1.000 2.7183 7.3891 20.0855 54.5982

Solution

i xi fi ∆fi ∆2 fi ∆3 f i ∆4 f i
0 0 1.0000
1.7183
1 1 2.7183 2.9525
4.6708 5.0731
2 2 7.3891 8.0256 8.7176
12.6964 13.7907
3 3 20.0855 21.8163
34.5127
4 4 54.5982

• Thus the polynomial is given by

       
s s s s
Pn (x) = 1 + 1.7183 + 2.9525 + 5.0731 + 8.7176 .
1 2 3 4

47
 • Suppose that we are interested in interpolating at x = 1.5, we will have
x − x0 1.5 − 0
s= = = 1.5
h 1
and it follows that

 
s
= 1.5.
1
 
s s(s − 1) 1.5(1.5 − 1)
= = = 0.375.
2 2 2
 
s s(s − 1)(s − 2) 1.5(1.5 − 1)(1.5 − 2)
= = = −0.0625.
3 6 6
 
s s(s − 1)(s − 2)(s − 3) 1.5(1.5 − 1)(1.5 − 2)(1.5 − 3)
= = = 0.0234375.
4 24 24

• Finally we have

P4 (1.5) = 4.5716.

• Since the data set was generated using the function ex we can calculate the error

|e1.5 − 4.5716| = 0.0899.

4.4 Newton’s Interpolation with equally spaced nodes: Backward differ-

ences

• while interpolating at some value of x near the end of the difference table, it is logical to
reorder the noes so that the end-differences can be used in comparison.

• Assume that the nodes nodes are equally spaced with spacing h, the backward difference
can be applied using the following steps:

Step 1 Set x0 ≡ xk , where xk is a point near the end of the table close to x = xi
Step 2 Label the nodes above x0 as x−1 , x−2 , · · · , x−k and the corresponding diagonal
divided differences as ∆f−1 , ∆f−2 , · · · , ∆f−k .
x−x0
Step 3 Set s = h . Then the general backward Newton interpolating polynomial is
given by
     
−s −s 2 n −s
Pn (x) = f0 − ∆f−1 + ∆ f−2 + · · · + (−1) ∆n f0 (76)
1 2 n
where
 
x − x0 s s(s − 1) · · · (s − k + 1)
s= , and, = , k = 1, 2, · · · n
h k k!

and ∆k f0 is the k th order difference of the given data.

48
Example

Using table of values in the previous example (forward difference), interpolate f (3.5) using
Newton’s backward difference method.

Solution

• Since the point x = 3.5 is near the end of the table, we take x0 = 4.

• Then
x − x0 3.5 − 4 1
s= = =− .
h 1 2

• Denote x−1 = 3, x−2 = 2, x−3 = 1 and

• ∆f0 = 54.5982, ∆f−1 = 34.5127, ∆f−2 = 21.8163, ∆f−3 = 13.7907 and ∆f−4 = 8.7176.

• Further, we have
 
−s 1
= −s = .
1 2
 
−s −s(−s − 1) 1
= =− .
2 2 8
 
s −s(−s − 1)(−s − 2) 0.5(0.5 − 1)(0.5 − 2)
= = = 0.0625.
3 6 6
 
s −s(−s − 1)(−s − 2)(−s − 3) 0.5(0.5 − 1)(0.5 − 2)(0.5 − 3)
= = = −0.15625.
4 24 24

• Therefore P4 (3.5) = 32.3908.

• Error is |e3.5 − 32.3908| = 0.7247.

4.5 Hermite Interpolation

• Given (n + 1) distinct points (x0 , f0 ), (x1 , f1 ), · · · , (xn , fn ), the Lagrange of Newton in-
terpolating polynomial Pn (x) matches the function values at the nodes, that is

Pn (xi ) = fi , i = 0, 1, · · · , n

• Some applications, however, demand that the interpolant be such that not only the
functional values are matched, but also the derivatives at the nodes. The procedure of
constructing such interpolates is called Hermite interpolation

49
 Generating Hermite Interpolating polynomial

• Given {xk , fk , fk0 }nk=0

Step 1. Form the Lagrangian polynomials Lk (x) and find their derivatives L0k (x),
k = 0, 1, · · · , n.
Step 2. Compute the polynomials {Bk }nk=0 and {B̂k }nk=0 as follows

Bk (x) = [(1 − 2(x − xk ))L0k (xk )]L2k (x)

B̂k (x) = (x − xk )L2k (x)

Step 3. Form the Hermite polynomial H2n+1 :

n
X n
X
H2n+1 (x) = f (xk )Bk (x) + f 0 (xk )B̂k (x).
k=0 k=0

Example

Construct the Hermite polynomial for the following data and use it to approximate f (1.5).

x f (x) f 0 (x)
1 0 1
2 0.6931 0.5

Solution

• Step 1. Compute the Lagrangian polynomials L0 (x) and L1 (x) and compute the deriva-
tives

x − x1 x−2
L0 (x) = = =2−x
x0 − x1 −1
x − x0 x−1
L1 (x) = = =x−1
x1 − x0 1

L00 (x) = −1, L01 (x) = 1

• Step 2. Compute the polynomials B0 (x), B1 (x), B̂0 (x) and B̂1 (x) from the Lagrange
polynomials:

B0 (x) = [1 − 2(x − x0 )L00 (x)]L20 (x)

= [1 − 2(x − 1)(−1)][2 − x]2
= [1 + 2(x − 1)][2 − x]2

50
 = (2x − 1)(2 − x)2
B1 (x) = [1 − 2(x − x1 )L01 (x)]L21 (x)
= [1 − 2(x − 2)][x − 1]2
= (5 − 2x)(x − 1)2
B̂0 (x) = (x − x0 )L20 (x) = (x − 1)(2 − x)2
B̂1 (x) = (x − x1 )L21 (x) = (x − 2)(x − 1)2

• Step 3. From the Hermite polynomial H3 (x) of degree 3:

H3 (x) = f0 B0 (x) + f1 B1 (x) + f00 B̂0 (x) + f10 B̂1 (x)

= 0.6931(5 − 2x)(x − 1)2 + (x − 1)(2 − x)2 + 0.5(x − 2)(x − 1)2 .

• Therefore H3 (1.5) = 0.4091.

• Accuracy check: The above table corresponds to f (x) = ln x. Thus the error is given by

| ln(1.5) − 0.4091| = |0.4055 − 0.4091| = 0.0036.

51
5 Numerical differentiation and integration
Numerical differentiation

• Numerical differentiation is the process of calculating the value of the derivative of a

function, y = f (x), at some assigned value of x from a given set of values (xi , yi ).

• The derivation of the formulae for numerical differentiation hinges on Taylor series ex-
pansions of f (x + h) and f (x − h), for h > 0 which are
1 2 00 1
f (x − h) = f (x) − hf 0 (x) + h f (x) − h3 f 000 (x) + · · · (77)
2! 3!
0 1 2 00 1 3 000
f (x + h) = f (x) + hf (x) + h f (x) + h f (x) + · · · (78)
2 3

Two point backward difference formula

• If the series (77) is truncated at the second derivative, we have

1 2 00
f (x − h) = f (x) − hf 0 (x) + h f (c0 ) (79)
2!
where c0 ∈ [x − h, x].

• Solving (79) for f 0 (x) gives

f (x) − f (x − h)
f 0 (x) = + ET (f, h) (80)
h
where ET (f, h) = 12 hf 00 (c0 ).

• Equation (80) is called the two point backward difference formula.

Two point forward difference formula

• If the series (78) is truncated at the second derivative, we have

1 2 00
f (x + h) = f (x) + hf 0 (x) + h f (c1 ) (81)
2!
where c1 ∈ [x, x + h].

• Solving (79) for f 0 (x) gives

f (x + h) − f (x)
f 0 (x) = + ET (f, h) (82)
h
where ET (f, h) = − 12 hf 00 (c1 ).

• Equation (82) is called the two point forward difference formula.

52
Central difference formula

• Subtracting equation (77) from (78) , will eliminate the even order derivatives, hat is,
f 00 (x), f iv (x),· · · and we obtain
2 3 000
f (x + h) − f (x − h) = 2hf 0 (x) + h f (x) + · · · (83)
3!
• Solving equation (83) gives the central difference formula
f (x + h) − f (x − h)
f 0 (x) = + E( f, h)
2h
where ET (f, h) = 16 h2 f 000 (c) with c ∈ [x − h, x + h].

Example

Consider the following table of values generated with f (x) = cos x.

x 0 0.2 0.4 0.6 0.8

f (x) 1.000000 0.980067 0.921061 0.825330 0.696707

Find an approximation of f 0 (0.4) (consider h = 0.2) using

1. The two point backward difference formula.

2. The two point forward difference formula.
3. The central difference formula.

Solution

• The two point backward difference formula is given by

f (x) − f (x − h) 0.921061 − 0.980067
f 0 (x) = = = −0.295028.
h 0.2
• We know that f 0 (x) = − sin x ⇒ f 0 (0.4) = − sin(0.4) = −0.389418

• Thus, the error is | − 0.389418 − (−0.295028)| = 0.09439.

Three point difference formula for first order derivatives

• Derivative formulae may also be obtained using interpolating polynomials.

• Consider an approximation to f (x) through x0 , x1 = x0 + h and x2 = x0 + 2h, using

Lagrange’s interpolating polynomial we have
(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )
P2 (x) = f0 + f1 + f2
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x0 )
(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )
= f0 − f1 + f2 . (84)
2h2 h2 2h2

53
 • Differentiating (84), with respect to x yields

(2x − x1 − x2 ) (2x − x0 − x2 ) (2x − x0 − x1 )

P20 (x) = 2
f0 − 2
f1 + f2 . (85)
2h h 2h2

• Evaluating (85) at x0 , x1 and x2 we obtain the three point formulae

−3f0 + 4f1 − f2
f 0 (x0 ) ≈ Forward difference formula
2h
−f0 + f2
f 0 (x1 ) ≈ , Central difference formula
2h
f0 − 4f1 + 3f2
f 0 (x2 ) ≈ , Backward difference formula.
2h

Example

Consider the following table of values generated with f (x) = cos x.

x 0 0.2 0.4 0.6 0.8

f (x) 1.000000 0.980067 0.921061 0.825330 0.696707

Find an approximation of f 0 (0.4) (consider h = 0.2) using

1. The three point backward difference formula.

2. The three point forward difference formula.

3. The central difference formula.

Solution

• The three point backward difference formula is given by

f0 − 4f1 + 3f2
f 0 (x) =
2h
f (x − 2h) − 4f (x − h) + 3f (x)
=
2h
f (0) − 4f (0.2) + 3f (0.4)
f 0 (0.4) = = −0.392712.
2(0.2)

• We know that y 0 = − sin x, so that f 0 (0.4) = − sin(0.4) = −0.3894183423086505.

• Thus, the error is 0.00329416.

• Remark Clearly from this example we observe that the three point formula give a
superior result compared to the two point formula.

54
Second order derivatives

• Adding (77) and (78) will eliminate the terms involving the odd order derivatives, that
is, f 0 (x), f 000 (x) · · ·, hence
1 4 (iv)
f (x + h) + f (x − h) = 2f (x) + h2 f 00 (x) + h f (x) + · · · (86)
12

• Solving (86) for f 0 (x) yields the central difference formula:

f (x + h) − 2f (x) + f (x − h)
f 00 (x) = + ET (f, h)
h2
2
where ET (f, h) = − h12 f (iv) (c) for some c ∈ [x − h, x + h].

Example

Given that f (x) = xex , use the three point formula with h = 0.1, 0.001, 0.001 and 0.0001 to
find approximations to f 00 (0.5). Compare the calculated value with the true value.

Solution

f (x + h) − 2f (x) + f (x − h)
f 00 (x) = .
h2
• Consider h = 0.01
f (0.51) − 2f (0.5) + f (0.49)
f 00 (0.5) =
(0.01)2
0.849298509 − 2(0.824360635) + 0.799834948
= = 4.121864950.
0.0001

• We known that f 00 (x) = ex (x + 2), so that f 00 (0.5) = 2.5e0.5 = 4.121803177.

• Error is |4.121803177 − 4.121864950| = 6.1773 × 10−5 .

h Calculated value Error

0.1 4.127988859 0.006185682
0.01 4.121864950 0.000061773
0.001 4.1279000 0.000005177
0.0001 4.121400000 0.000403177

Remark: Note that the error initially decreases as h decreases. However, it starts increasing
again if h is reduced below 0.001. This is because of rounding error is more prominent for
smaller step size.

55
Error analysis

• Let f (x) = fk + ek where ek is the rounding error in computing f (xk ) and fk is the
approximate. The error in three point central difference formula for f 00 (x) is

E(h) = rounding error + truncation error

e1 − 2e0 + e−1 h2 (iv)

= − f (c)
h2 12
• If we assume that the round-off error are bounded by some positive number e, and
|f (iv) (x)| ≤ M , then
4e M h2
|E(h)| ≤ + .
h2 12
4e
• If h is small, then the contribution of h2
due to round-off is large. However, when h is
M h2
large the contribution of 12 is large. Thus, the optimum size of h can be found by
minimizing

4e M h2
g(h) = + .
h2 12
• Thus,
8e M h
g 0 (h) = − + =0
h3 6
which yields

 1
48e 4
h=
M
Example

Determine the optimal step size h for solving f (x) = xex in the interval [0.4, 0.6].

Solution

• Observe that |f (iv) (x)| = |(x + 4)ex | ≤ 8.4 on the given interval.

• Note the for the calculation in the previous example the solutions were given correct to
9 decimal places which corresponds to the maximum round-off error e = 0.5 × 10−9 .

• Thus
 14
48 × 0.5 × 10−9

h= = 0.0073.
8.4
Remark Observe that the solution here is consistent with our observation using direct calcu-
lations.

56
Numerical Integration
Rb
• It is often either difficult or impossible to evaluate definite integral of the form a f (x)dx.

• It is well known that the definite integral may be interpreted as the area under the curve
y = f (x) for a ≤ x ≤ b and may be evaluated by subdivision.
b−a
• Le h = n , be the interval length where n is the number of strips.

• For any subinterval xj ≤ x ≤ xj+n of the interval a ≤ x ≤ b, we may approximate f (x)

by interpolating polynomial Pn (x).

• Thus, we have the approximation

Z xj+n Z xj+n
f (x)dx ≈ Pn (x)dx
xj xj

which will be a good approximation if n is chosen so that the error [f (x) − Pn (x)] in each
tabular subinterval xj+k−1 ≤ x ≤ xj+k (k = 1, 2, 3 · · · n) is sufficiently small.

Trapezoidal

• Consider the following data set

x x0 x1
f (x) f0 f1

• The Lagrange polynomial of the above data is give by

(x − x1 ) (x − x0 )
P1 (x) = f0 + f1 (87)
(x0 − x1 ) (x1 − x0 )

• Let x1 = x0 + h so that
(x − x1 ) (x − x0 )
P1 (x) = − f0 + f1 (88)
h h

• To derive the Trapezoidal rule we integrate (88), that is,

Z x1 Z x1 Z x1  
(x − x1 ) (x − x0 )
f (x)(x)dx ≈ P1 (x)dx = − f0 + f1 dx
x0 x0 x0 h h
Z hh
x x i
≈ f1 − f0 + f0 dx.
0 h h
h
= (f1 + f0 ).
2

57
 • In general, the total area under the curve can therefore be approximated by the sum of
the areas of the n trapezia, that is,
Z b Z x1 Z x2 Z xn 1
f (x)dx = f (x)dx + f (x)dx + · · · + f (x)dx
a x0 x1 xn−1

h h h h
≈ (f0 + f1 ) + (f1 + f2 ) + (f2 + f3 ) + · · · (fn−1 + fn )
2 2 2 2
h
= [f0 + 2(f1 + f2 + · · · + fn−1 ) + fn ]
2
n−1
" #
1 X 1
= h f (x0 ) + f (xi ) + f (xn ) .
2 2
i=1

Example

Use the Trapezoidal rule to evaluate

Z 1
1
I= dx
0 1+x
with h = 0.25.

Solution

• We have
4 2 4 1
f (0) = 1, f (0.25) = , f (0.5) = , f (0.75) = , f (1) = .
5 3 7 2

   
0.25 4 2 4 1
I = 1+2 + + +
2 5 3 7 2

= 0.697023809.

• We know that the analytical solution is ln(1 + x)|10 = ln 2

• Error is | ln 2 − 0.697023809| = 0.00387663.

Error of the Trapezoidal rule

• The trapezoidal rule correspond to a rather crude polynomial approximation (a straight

line) between successive points xj and xj+1 = xj + h, and hence can only be accurate for
sufficiently small h.

• Approximate (upper) bound the error may be derived from the Taylor expansion as
follows
t2 00
f (xr + t) = f (xr ) + tf 0 (xr ) + f (xr ) + · · ·
2!

58
 = a0 + a1 t + a2 t2 + a3 t3 + · · ·
1 m 1 m
where am = m! f (xr ) = m! fr

• Now if we let x = xr + t, then

Z xr +t Z h
f (x)dx = f (xr + t)dt
xr 0
Z h
= [a0 + a1 t + a2 t2 + a3 t3 + · · ·]dt
0

1 1 1
= a0 h + a1 h2 + a2 h3 + a3 h4 (89)
2 3 4
• In the trapezoidal rule, we approximate the integral by
Z xr+1
h
f (x)dx = [fr + fr+1 ] (90)
xr 2

• Since fr = a0 and fr+1 = f (xr + h) we have

h h
[fr + fr+1 ] = [a0 + (a0 + a1 h + a2 h2 + · · ·)]
2 2
1 1
= a0 h + a1 h2 + a2 h3 + · · · (91)
2 2
• If we compare (89) and (90), ignoring higher order terms we observe that the error in
this cases is
1 1
E = a2 h3 − a2 h3
3 2
1
= − a 2 h3
6
1 3 00
= − h fr . (92)
12
• Equation (92) represents he error in one strip. Therefore the error in n strips is given by
nh3 00
En = − f (93)
12 r
b−a
• But n = h , thus equation (93) can be rewritten as
(b − a) 2
En = h max |f 00 (x)|.
12 a≤x≤b

Example
R1 2
If the Trapezoidal rule is to be used to compute 0 e−x dx with an error of at most 0.5 × 10−4 ,
how many points should be used?

59
Solution
2
• We know that f 00 (x) = (4x2 − 2)e−x , so that |f 00 (x)| ≤ 2

• Thus
1 3 00
− h fr < 0.5 × 10−4
12

h < 0.01732.

b−a 1
• But h = n = n ⇒ n ≥ 59.

Simpson’s rule

• If it is undesirable to increasingly subdivide an interval a ≤ x ≤ b in order to get

increasingly accurate quadrature, the alternative is to use an approximating polynomial
of higher degree.

• In our previous discussion we observed that the trapezoidal rule utilizes a linear approx-
imation (two nodes). So trying a quadratic interpolating polynomial (3 nodes) yields
Simpson’s rule.

• Consider the following table

x x0 x1 x2
f (x) f0 f1 f2

• Let x1 = x0 + h and x2 = x1 + h

• The Lagrange’s interpolating polynomial for the above data set is

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

P2 (x) = 2
f0 − 2
f1 + f2 (94)
2h h 2h2
• Now

Z x2 Z x0 +2h
f (x)dx = P2 (x)dx
x0 x0

h
= (f0 + 4f1 + f2 ).
3
• In general, the total area under the curve can therefore be approximated by the sum of
the areas of the n double strips, that is,
Z b Z x2 Z x4 Z xn 1
f (x)dx = f (x)dx + f (x)dx + · · · + f (x)dx
a x0 x2 xn−2

60
 h h h
≈ (f0 + 4f1 + f2 ) + (f2 + 4f3 + f4 ) + · · · + (fn−2 + 4fn−1 + fn )
3 3 3
h
= [f0 + 4(f1 + f3 + f5 · · · + fn−1 ) + 2(f2 + f4 + f6 + · · · + fn−2 ) + fn ]
3
n n−1
" #
h X X
= f (x0 ) + 4 f (x2i−1 ) + 2 f (x2i ) + f (xn ) .
3
i=1 i=1

Example

Use the Simpson’s rule to evaluate

Z 1
1
I= dx
0 1+x

with h = 0.25.

Solution

• We have
4 2 4 1
f (0) = 1, f (0.25) = , f (0.5) = , f (0.75) = , f (1) = .
5 3 7 2

       
0.25 4 2 2 4 1
I = 1+4 + + +4 +
3 5 3 3 7 2

= 0.693254.

• We know that the analytical solution is ln(1 + x)|10 = ln 2

• Error is | ln 2 − 0.697023809| = 0.000106788.

• Remark Observe that the Simpson gives a better result compared to the Trapezoidal
rule.

Error of the Simpson’s rule

• Instead of considering the single strip from xr to xr+1 , we consider two strips from xr−1
to xr+1 .

• The Taylor series expansion of f (xr + t) is given by

f (xr + t) = a0 + a1 t + a2 t2 + a3 t3 + · · · (95)

61
• The Taylor series expansion of f (xr − t)

f (xr − t) = a0 − a1 t + a2 t2 − a3 t3 + · · · (96)

1 m 1 m
where am = m! f (xr ) = m! fr

• Then
Z xr+1 Z h
f (x)dx = f (xr + t)dt
xr−1 −h

= [a0 + a1 t + a2 t2 + a3 t3 + · · ·]h−h

2 2
= 2a0 h + a2 h3 + a4 h5 + · · · (97)
3 5

• Note that fr = f (xr ) = a0 . Adding the Taylor series expansions of f (xr −t) and f (xr +t),
we have

fr−1 + fr+1 = 2a0 + 2a2 t2 + 2a4 t4 + · · ·

= 2a0 + 2a2 h2 + 2a4 h4 + · · ·

1 2 2
h[fr−1 + 4fr + fr+1 ] = 2a0 h + a2 h3 + a4 h5 + · · · (98)
3 3 3

• Comparing (97) and (98) the error of the Simpson rule is given by
 
2 2
E = − a4 h5
5 3
4
= − a4 h5
15
1 5 (iv)
= − h fr
90

• Ignoring the higher-order terms, we conclude that the approximate bound on this error
Rb
is estimating a f (x)dx by Simpson’s rule (with n2 subintervals of width 2h) is
 n   h5  (b − a)h4
max |f (iv) (x)| = max |f (iv) (x)|.
2 90 a≤x≤b 180 a≤x≤b

62
Gaussian quadrature formulas

• The numerical integration procedures previously discussed (namely the trapezoidal rile
and Simpson’s rule) involve equally spaced values of the argument, and from these nu-
merical formulas we observed that
Z b Z b
f (x)dx ≈ Pn (x)dx = w0 f0 + w1 f1 + · · · wn fn . (99)
a a

Here xi ∈ [a, b] are called nodes, and wi0 s are the weights. For example, on the trapezoidal
rule we have
b−a h
x0 = a, x1 = b, w0 = w1 = = .
2 2

• From this discussion, we can deduce the Interpolatory Quadrature as follows

Z b
f (x)dx ≈ Qn (f ) + Rn (f )
a
n
X
= wi f (xi ) + Rn (f ), (100)
r=1

with
Z b
wr = li (x)dx
a

• Equation (100) is called the n− point quadrature rule, Qn (f ) is the quadrature sum and
Rn (f ) is the error term.

Example

Determine the quadrature formula of the form

n
X
Qn (f ) = wi f (xi )
i=1

given the interval [−2, 2] and the nodes −1, 0, 1. Hence, evaluate the integral
Z 2
I= ex dx.
−2

Solution

We have 3 nodes, thus the polynomial takes the form

P2 (x) = L0 f (x0 ) + L1 f (x1 ) + L2 f (x2 )

63
 with
x(x − 1) x(x + 1)
L0 (x) = , L1 (x) = −(x − 1)(x + 1), L2 = .
2 2

Now, we compute the weights, w1 , w2 , w3

Z 2 −2
x(x − 1)
Z
8
w1 = L0 (x)dx = dx =
−2 2 2 3
Z 2 Z 2
4
w2 = L1 (x)dx = −(x − 1)(x + 1)dx = −
−2 −2 3
Z 2 Z 2
x(x + 1) 8
w3 = L2 (x)dx = dx =
−2 −2 2 3

Thus,
8 4 8 8 4 8
Qn (f ) = f (x0 ) − f (x1 ) + f (x1 ) = f (−1) − f (0) + f (1). (101)
3 3 3 3 3 3

Hence
Z 2
8 4 8
ex dx ≈ e−1 − e0 + e1 = 6.89643.
−2 3 3 3

The error for such an approximation is given by

(e2 − e−2 ) − 6.89643 = 0.357291.

Gaussian quadrature rule (The method of undetermined coefficients)

• Form the derivation of the interpolatory quadrature, we note that we fix the points xi
and we adjust the weights wi .

• If now we allow the nodes xi also to be adjusted, we gain higher degree of freedom, and
higher accuracy as well.

• One chooses the nodes xi and the weights wi (i = 0, 1, 2, · · · , n) such that equation (99)
gives the exact values for polynomial functions f (x) of highest possible degree m

Z b Z b
f (x)dx ≈ Pn (x)dx = w0 f0 + w1 f1 + · · · wn fn , if, f ∈ Pm (102)
a a

• Here, m is called the degree of precision. To fix the idea, we consider the interval [−1, 1],
for n = 1, thus we have have two nodes and two weights, that is
Z b
f (x)dx ≈ w0 f (x0 ) + w1 f (x1 ). (103)
a

64
 • If (99) gives the exact value for polynomials of degree m the it must be exact for functions
x0 = 1, x1 , x2 , x3 , · · · , xm . Thus, for, equation (103) we need 4 equations, since we 4
unknowns, therefore we must have m = 2n − 1 = 3.

• Thus, polynomials of degree 2n−1 is the highest degree of precision which can be obtained
using n such nodes.

• If the nodes are given and the weights are unknown the highest degree polynomial is
given by n − 1.

Example

Find a quadrature rule of the form

Z 1  
1
f (x)dx ≈ A0 f (0) + A1 f + A2 f (1)
0 2

Solution

• Since we have 3 nodes the quadrature should be exact for all polynomials of degree ≤ 2,
that is, we need to consider, x0 , x1 , x2 . Thus

 R1

 f (x) = 1 ⇒ 0 1dx = 1 = A0 + A1 + A2




 R1
f (x) = x ⇒ 0 xdx = 21 = 12 A1 + A2





 R1
f (x) = x2 ⇒ 0 x2 dx = 13 = 14 A1 + A2 .


Solving the system

 A + A + A = 1,
 0 1 2




1

1
A1 + A2 = , (104)
 2 2


 1A + A = 1



1 2
4 3
gives
1 2 1
A0 = , A1 = , A2 = (105)
6 3 6
Therefore the desired quadrature rule takes the form
Z 1    
1 1
f (x)dx ≈ f (0) + 4f + f (1) .
0 6 2

Given f (x) we can approximate the solution of the integral.

65
Exercise

Find a quadrature rule of the form

Z 1
f (x)dx ≈ A0 f (x0 ) + A1 f (x1 ).
−1

6 Numerical solutions of ordinary differential equations

• A differential equation is a equation in which an unknown function appears together with
one or more of its derivatives

• A solution of a differential equation is a specific function which satisfies that equation.

• A number of problems in science and technology can be presented in the form of ordinary
differential equations.

• The analytical methods of solving differential equations are applicable to a limited class
of equations. Quite often, differential equations modelling real world problems, do not
belong to any of these familiar types and as such one is obliged to resort to numerical
methods for a solution.

Taylor series method

• Consider an initial value problem (IVP)

dy
= f (x, y), y(a) = y0 , x ∈ [a, b]. (106)
dx
• We seek to find a solution of equation (106) by some numerical methods.

• If we subdivide the interval [a, b] into n subintervals and select mesh points, we have

xk+1 = xk + h, k = 0, 1, 2, · · · , n − 1

where
b−a
h= , x0 = a, and xn = b.
n
• Using Taylor’s theorem, we have

y(xk+1 ) = y(xk + h)
1 2 00 1
= y(xk ) + hy 0 (xk ) + h y (xk ) + · · · + hn y (n) (xk ) + · · · (107)
2! n!
• If the step size h is chosen small enough, we may neglect higher order terms of (107),
that is,

y(xk+1 ) = y(xk + h)
1 2 00 1
= y(xk ) + hy 0 (xk ) + h y (xk ) + · · · + hn y (n) (xk ) (108)
2! n!

66
 • Equation (108) is called the Taylor’s formula of order n and it can be used to compute
the numerical solution of the IVP (106).

• In equation (108), the error is of the order O(hn+1 ).

Example

Use Taylor’s method of order 2 to solve the following initial value problem

y 0 = y + ex , y(0) = 1, x ∈ [0, 0.003], h = 0.001.

Solution

• Taylor’s method of order 2 is given by

yk+1 = y(xk ) + hy 0 (xk ) + hy 00 (xk ).

• We know that

y 0 (xk ) = y(xk ) + exk .

y 00 (xk ) = y 0 (xk ) + exk
= [y(xk ) + exk ] + exk
= y(xk ) + 2exk .

• Thus,
1
y(xk+1 ) = y(xk ) + hy[y(xk ) + exk ] + h2 [y(xk ) + 2exk ]
2
 
1
= 1 + h + h2 yk + (h + h2 )exk .
2

• For k = 0 we have
 
1
y1 = 1 + 0.001 + (0.001) y0 + (0.001 + (0.001)2 )ex0
2
2
= 1.0010005y0 + 0.001001ex0
= 1.001(1) + 0.001001(1)
= 1.0020015.

• For k = 1 we have

y2 = 1.001y1 + 0.001001ex1
= 1.0010005(1.0020015) + 0.001001e0.001
= 1.004006004001417.

67
 • For k = 2 we have

y3 = 1.0010005y1 + 0.001001ex1
= 1.0060135160117556.

• We can easily show that the analytical solution is

y(x) = ex (1 + x).

• Thus, the error is given by

x Estimate Exact Error

0.001 1.0020015 1.00200150066875 6.66875 × 10−10
0.002 1.004006004001417 1.0040060053366684 1.33525 × 10−9
0.003 1.0060135160117556 1.0060135180168872 2.00513 × 10−9

Euler’s method

• The Taylor’s method of order 1 is known as the Euler’s method, that is

y(xk + h) = y(xk ) + hy 0 (xk ) (109)

• Recall that y 0 (x) = f (x, y), thus (109) becomes

y(xk + h) = y(xk ) + hf (xk , yk ).

Example

Use Euler’s method to solve the following initial value problem

y 0 = y + ex , y(0) = 1, x ∈ [0, 0.003], h = 0.001.

Solution

• f (xk , yk ) = yx + exk so the Euler’s method is

yk+1 = yk + h(yk + exk )

= hexk + (1 + h)yk .

• Thus

x Estimate Exact Error

0.001 1.001 1.00200150066875 0.0010015006668750992
0.002 1.0030020005001663 1.0040060053366684 0.0010040048365020926
0.003 1.0050080070051697 1.0060135180168872 0.0010055110117175303

68
Modified Euler’s method

• Consider the initial value problem y 0 = f (x, y), x ∈ [a, b] with y(a) = y0 and y(b) = yn .
b−a
• If we subdivide the interval [a, b] into n subintervals of equal length h and set h = n ,
so that xk+1 = xk + h for k = 0, 1, 2, 3, · · · .

• Integrating y 0 = f (x, y) over [xk , xk+1 ] we get

Z xk+1 Z xk+1
0
y dx = f (x, y)dx
xk xk

h
y(xk+1 ) − y(xk ) = [f (xk , yk ) + f (xk+1 , yk+1 )]
2
h
= y(xk ) + [f (xk , yk ) + f (xk+1 , yk+1 )]. (110)
2

• Equation(110) is called the modified Euler’s method or the Euler-Cauchy’s method.

• It is also known as the predictor-corrector method.

Example

Use Euler-Cauchy’s method to solve the following initial value problem

y 0 = −2x − y, y(0) = −1, x ∈ [0, 0.2], h = 0.1.

Solution

• For k = 0 we have
h
y1 = y0 + [f (x0 , y0 ) + f (x1 , y1 )]
2

• We know that f (xk , yk ) = −2xk − yk so that

f (x0 , y0 ) = −2x0 − y0 ⇒ f (0, −1) = −2(0) − (−1) = 1

• Since we are not aware of y1 for to use on f (x1 , y1 ) we will first use Euler’s method to
find y1 , that is;

yk+1 = yk + hf (xk , yk )
y1 = y0 + hf (x0 , y0 )
= −1 + 0.1(1) = −0.9.

• Now

f (x1 , y1 ) = f (0.1, −0.9) = −2(0.1) − (−0.9) = 0.7

69
 • Thus
0.1
y1 = −1 + [1 + 0.7] = −0.915.
2
• For k = 1 we have
h
y2 = y1 + [f (x1 , y1 ) + f (x2 , y2 )]
2
•

f (x1 , y1 ) = f (0.1, −0.915) = −2(0.1) − (−0.915) = 0.715

y2 = y1 + hf (x1 , y1 ) = −0.86 + 0.2[−2(0.2) − (−0.86)] = −0.8435.
f (x2 , y2 ) = f (0.2, −0.8435) = −2(0.2) − (−0.8435) = 0.4435.

• Thus
0.1
y2 = −0.915 + [0.715 + 0.4435] = −0.857075.
2
• It can easily be verified that the analytical solution is

y(x) = 2(1 − x) − 3e−x

• Thus

x Estimate Exact Error

0.1 -0.915 -0.9145122541078787 0.00048774589212130515
0.2 -0.857075 -0.8561922592339456 0.0008827407660544351

Runge-Kutta methods

Runge-Kutta method of order 1

• The Taylor’s method of order 1 is also known as the Runge-kutta method of order 1 and
is given by

yk+1 = yk + hf (xk , yk ).

Runge-Kutta method of order 2

• It is given by
1
yk+1 = yk + h[k1 + k2 ]
2
where

k1 = f (xk , yk )
k2 = f (xk+1 , yk+1 )
= f (xk + h, yk + hk1 ).

Remark: Here we observe that Runge-Kutta method of order 2 (RK-2) is similar to the
modfied Euler’s method.

70
Runge-Kutta method of order 3

• It is given by
h
yk+1 = yk + [k1 + 4k2 + k3 ]
6
where

k1 = f (x
 k , yk ) 
h h
k2 = f xk + , yk + k1
2 2
k3 = f (xk + h, yk − hk1 + 2hk2 ).

Example

Use Runge-Kutta’s method of order 3 (RK-3) to solve the following initial value problem

y 0 = −2x − y, y(0) = 1, x ∈ [0, 0.2], h = 0.1.

Solution

• We have f (xk , yk ) = −2xk − yk so that

 0 , y0 ) = f (0, −1) 
k1 = f (x = −2(0) − (−1) = 1,
h h
k2 = f x0 + , y0 + k1 = f (0.05, −0.95) = −2(0.05) − (−0.95) = 0.85
2 2
k3 = f (x0 + h, y0 − hk1 + 2hk2 ) = f (0.1, −0.93) = −2(0.1) − (−0.93) = 0.73

• Thus

0.1
y1 = −1 + [1 + 4(0.85) + 0.73] = −0.9145.
6

• For k = 1

 1 , y1 ) = f (0.1, −0.9145)
k1 = f (x  = −2(0.1) − (−0.9145) = 0.7145,
h h
k2 = f x1 + , y1 + k1 = f (0.15, −0.878775) = −2(0.15) − (−0.878775) = 0.578775.
2 2
k3 = f (x1 + h, y1 − hk1 + 2hk2 ) = f (0.2, −0.870195) = −2(0.1) − (−0.93) = 0.470195

• Thus

0.1
y2 = −0.9145 + [0.7145 + 4(0.578775) + 0.470195] = −0.8561700833333333.
6

71
 • Therefore

x Estimate Exact Error

0.1 -0.9145 -0.9145122541078787 0.000012254107878750808
0.2 -0.8561700833333333 -0.8561922592339456 0.0000221759006122646

Runge-Kutta method of order 4 (RK-4)

• It is given by
h
yk+1 = yk + [k1 + 2k2 + 2k3 + k4 ]
6
where

k1 = f (x
 k , yk ) 
h h
k2 = f xk + , yk + k1
 2 2 
h h
k3 = f xk + , yk + k2
2 2
k4 = f (xk + h, yk + hk3 ).

Example

Use Runge-Kutta’s method of order 3 (RK-3) to solve the following initial value problem

y 0 = −2x − y, y(0) = 1, x ∈ [0, 0.2], h = 0.1.

Solution

• We have f (xk , yk ) = −2xk − yk so that

 0 , y0 ) = f (0, −1) 
k1 = f (x = −2(0) − (−1) = 1,
h h
k2 = f x0 + , y0 + k1 = f (0.05, 0.95) = 0.85,
 2 2 
h h
k3 = f x0 + , y0 + k2 = f (0.05, −0.9575) = 0.8575,
2 2
k4 = f (x0 + h, y0 + hk3 ) = f (0.1, −0.91425) = 0.71425.

• Thus
0.1
y1 = −1 + [1 + 2(0.85) + 2(0.8575) + 0.71425] = −0.9145125.
6

• Therefore the error is 2.45892 × 10−7 .

72
Second order differential equations

• Majority of real life problems can be presented as a system of first order differential
equations, that is;

y10 = f1 (x, y1 , y2 , · · · , yp ),
y20 = f2 (x, y1 , y2 , · · · , yp ),
.. .
. = ..
yp0 = fp (x, y1 , y2 , · · · , yp ).

• If a real world problem is modeled as a single higher order differential equation, it is

possible for one to reduce the differential equation to a system of first order differential
equations.

Example

Reduce the second order equation

y 00 + 3y 0 + 2y = ex , y(0) − 1, y 0 (0) = 2

to a system of first order equations.

Solution

• Let u = y and v = y 0 so that u0 = y 0 = v and v 0 = y 00

• Thus

u0 = v u(0) = 1
0
v = −2u − 3v + ex , v(0) = 2

• Once a higher-order system is reduced to a system of first-order equations, we can then

used numerical techniques discussed earlier in this chapter to obtain solutions.

Example

Consider the initial value problem

y 00 − y = x, x ∈ [0, 1], y(0) = 0, y 0 (0) = 1.

(i) Determine the analytical solution of the differential equation.

(ii) Reduce the above differential equation to a system of first order differential equations.

(iii) Use Euler’s method with h = 0.05 to estimate y(0.1).

(iv) Compute the error of y(0.1) and y 0 (0.1).

73
Solution

(i) It can easily be verified that the analytical solution is

y(x) = ex − e−x − x.

(ii) Let u = y and v = y 0 so that u0 = y 0 = v and v 0 = y 00

u0 = v u(0) = 0
0
v = u + x, v(0) = 1

(iii) Let f (xn , un , vn ) = vn and g(xn , un , vn ) = un + xn . Thus,

un+1 = un + hf (xn , un , vn )
vn+1 = vn + hg(xn , un , vn )

• For n = 0, we have

u1 = u0 + hf (x0 , u0 , v0 )
= 0 + 0.05f (0, 0, 1)
= 0.05
v1 = v0 + hg(x0 , u0 , v0 )
= 1 + 0.05g(0, 0, 1)
= 1

• For n = 1, we have

u2 = u1 + hf (x1 , u1 , v1 )
= 0.05 + 0.05f (0.05, 0.05, 1)
= 0.05 + 0.05(1) = 0.1
v2 = v1 + hg(x1 , u1 , v1 )
= 1 + 0.05g(0.05, 0.05, 1)
= 1 + 0.05(0.1)
= 1.005.

• The analytical solution is y(0.1) = 0.100334.

• Therefore the error is 0.0003335.

• Observe that y 0 (x) = ex + e−x − 1, so that y 0 (0.1) = 1.0100083361116072.

• Hence the error is 0.00500834.

74
Example

Consider the initial value problem

y 00 − y = x, x ∈ [0, 1], y(0) = 0, y 0 (0) = 1.

(i) Use RK-4 method with h = 0.05 to estimate y(0.1).

(ii) Compute the error of y(0.1) and y 0 (0.1).

Solution

• We need to reduce our differential equation to a system of first order differential equations
as done in the previous example. Thus,

u0 = v u(0) = 0
v 0 = u + x, v(0) = 1

• Let

f (xn , un , vn ) = vn
g(xn , un , vn ) = un + xn

• To apply the Rk-4 method we have

h
un+1 = un + [k1 + 2k2 + 2k3 + k4 ]
6
h
vn+1 = vn + [l1 + 2l2 + 2l3 + l4 ]
6
where

k1 = f (xn , un , vn ),

l1 = g(xn , un , vn ),
 
h h h
k2 = f xn + , un + k1 , vn + l1 ,
2 2 2
 
h h h
l2 = g xn + , un + k1 , vn + l1 ,
2 2 2
 
h h h
k3 = f xn + , un + k2 , vn + l2 ,
2 2 2
 
h h h
l3 = g xn + , un + k2 , vn + l2 ,
2 2 2

75
 k4 = f (xn + h, un + hk3 , vn + hl3 ) ,

l4 = g (xn + h, un + hk3 , vn + hl3 ) .

• Now

k1 = f (x0 , u0 , v0 ) = f (0, 0, 1) = 1,
l1 = g(x0 , u0 , v0 ) = g(0, 0, 1) = u0 + x0 = 0.

k2 = f (0.025, 0.025, 1) = 1,
l2 = g(0.025, 0.025, 1) = 0.05.

k3 = f (0.025, 0.025, 1.00125) = 1.00125,

l3 = g(0.025, 0.025, 1.00125) = 0.05.

k4 = f (0.05, 0.0500625, 1.0025) = 1.0025,

l4 = g(0.05, 0.025, 1.0025) = 0.05 + 0.0500625 = 0.1000625

• Thus
0.05
u1 = 0 + [1 + 2(1) + 2(1.00125) + 1.0025] = 0.0500417.
6
0.05
v1 = 1 + [0 + 2(0.05) + 2(0.05) + 0.1000625] = 1.0025005208333333
6

• Recall that the exact solution is y(0.05) = e0.05 − e−0.05 − 0.05 = 0.0500416718753101.

• Therefore the error is 5.20864 × 10−9 .

76