NUMERICAL ANALYSIS I

f (x )
xn+1 = xn − f 0(xn )
n

[Link]
Department of Mathematics
Makerere University

D.W-Ddumba numerical analysis i Page 1 of 133

Page 2 of 133 D.W-Ddumba numerical analysis i
Contents

1 Introduction 5
1.1 What is Numerical Analysis? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Two issues of Numerical Analysis: . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Advantages and Disadvantages of Numerical Analysis: . . . . . . . . . . . . . . . 7
1.3.1 Hard or impossible classical solutions . . . . . . . . . . . . . . . . . . . . 7
1.3.2 Tabulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Important Notes: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Numerical Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.1 Sources of Errors: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.2 Types of Errors: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Numerical Integration 11
2.1 Manual Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Trapezoidal/Trapezium rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.1 Composite Trapezoidal Rule . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Simpson’s rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.1 Composite Simpson’s rule . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.2 Program-FORTRAN (an alternative to MAPLE) . . . . . . . . . . . . . 21
2.4 Mid-Point Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.1 Error in midpoint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3 Solution to Non-Linear Equations 27

3.1 Newton Raphson’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.1.1 Derivation of the Newton’s method . . . . . . . . . . . . . . . . . . . . . 28
3.1.2 General Newton’s algorithm for extracting roots of positive numbers . . . 41
3.1.3 Using Newton’s general formula for roots in finding reciprocals of numbers 43
3.1.4 Advantages of Newton Raphson’s method . . . . . . . . . . . . . . . . . 44
3.1.5 Some limitations of the Newton Raphson’s method . . . . . . . . . . . . 44
3.2 Bisection Method (Interval Halving) . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2.3 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2.4 Explanation on the bisection method . . . . . . . . . . . . . . . . . . . . 46
3.2.5 Advantages of the Bisection method . . . . . . . . . . . . . . . . . . . . . 49
3.2.6 Disadvantages of the Bisection method . . . . . . . . . . . . . . . . . . . 49
3.3 Secant Method (Chords Method) . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.3.1 Derivation of the Secant method . . . . . . . . . . . . . . . . . . . . . . . 56
3.3.2 Advantages and Disadvantages of the Secant method . . . . . . . . . . . 62
3.4 The Regula-Falsi method (Method of false position) . . . . . . . . . . . . . . . . 63
3.4.1 Geometric representation and derivation of the Regula falsi algorithm . . 63

3
 CONTENTS

3.4.2 Order of Convergence of the Regula algorithm . . . . . . . . . . . . . . . 73

3.4.3 Advantages and disadvantages of the regula falsi algorithm . . . . . . . . 74
3.5 Successive Substitution (Fixed Point Method) . . . . . . . . . . . . . . . . . . . 76
3.5.1 Background knowledge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.5.2 Successive Substitutions . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.5.3 Convergence criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4 Interpolation 83
4.1 Review- Linear interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3 History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.4 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.5 Lagrange interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.5.1 Alternatively . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.5.2 Examples of interpolating polynomials . . . . . . . . . . . . . . . . . . . 88
4.5.3 Error analysis in Lagrange’s interpolation . . . . . . . . . . . . . . . . . 92
4.5.4 Rounding errors in Lagrange polynomials . . . . . . . . . . . . . . . . . . 94

5 Numerical Differentiation 101

5.1 Why Numerical techniques for finding derivatives . . . . . . . . . . . . . . . . . 101
5.2 Analytic definition of a derivative as compared to a numerical definition . . . . . 101
5.3 Forward difference approximation . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.3.1 Analytical derivation of the forward difference approximation . . . . . . 104
5.4 Backward difference approximation . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.4.1 Analytical derivation of the backward difference approximation . . . . . . 106
5.5 The Central difference approximation . . . . . . . . . . . . . . . . . . . . . . . . 108
5.5.1 Analytical derivation of the central difference approximation . . . . . . . 109
5.6 Comparision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.7 The second derivative approximation . . . . . . . . . . . . . . . . . . . . . . . . 111
5.7.1 Error analysis in numerical differentiation . . . . . . . . . . . . . . . . . 112

6 Finite difference operators 115

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.2 Finite differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.3 Finite difference operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.3.1 The forward difference operator . . . . . . . . . . . . . . . . . . . . . . . 116
6.3.2 The backward difference operator . . . . . . . . . . . . . . . . . . . . . . 117
6.3.3 The shift operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.3.4 The central difference operator . . . . . . . . . . . . . . . . . . . . . . . 120
6.3.5 Averaging operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.4 Finite difference tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

7 Sample Questions 125

8 Further Reading 133

Page 4 of 133 D.W-Ddumba numerical analysis i

Chapter 1

Introduction

Most real mathematical problems do not have analytical solutions. However, they do have
real solutions. In order to obtain these solutions we must use other methods such as graphical
representations, or numerical analysis. Numerical analysis is the mathematical method that
uses numerical approximations to obtain numerical answers to the problem. Numerical analysis
also considers the accuracy of an approximation, and when the approximation is good enough.
Numerical answers are useful because we use numbers to build our world, not with the exact
π
analytical solution, such as √e27
The ever-increasing advances in computer technology has enabled many in science and engi-
neering to apply numerical methods to simulate physical phenomena. Numerical methods are
often divided into elementary ones such as finding the root of an equation, integrating a func-
tion or solving a linear system of equations to intensive ones like the finite element method.
Intensive methods are often needed for the solution of practical problems and they often re-
quire the systematic application of a range of elementary methods, often thousands or millions
of times over. In the development of numerical methods, simplifications need to be made to
progress towards a solution: for example general functions may need to be approximated by
polynomials and computers cannot generally represent numbers exactly anyway. As a result,
numerical methods do not usually give the exact answer to a given problem, or they can only
tend towards a solution getting closer and closer with each iteration. Numerical methods are
generally only useful when they are implemented on computer using a computer program-
ming language .
The study of the behavior of numerical methods is called numerical analysis. This is a mathe-
matical subject that considers the modeling of the error in the processing of numerical methods
and the subsequent re-design of methods.
Numerical analysis involves the study of methods of computing numerical data. In many prob-
lems this implies producing a sequence of approximations; thus the questions involve the rate
of convergence, the accuracy (or even validity) of the answer, and the completeness of the re-
sponse. (With many problems it is difficult to decide from a program’s termination whether
other solutions exist.) Since many problems across mathematics can be reduced to linear alge-
bra, this too is studied numerically; here there are significant problems with the amount of time
necessary to process the initial data. Numerical solutions to differential equations require the
determination not of a few numbers but of an entire function; in particular, convergence must
be judged by some global criterion. Other topics include numerical simulation, optimization,
and graphical analysis, and the development of robust working code.
Numerical linear algebra topics: Solutions of linear systems AX = B, eigenvalues and eigenvec-
tors, matrix factorizations. Calculus topics: numerical differentiation and integration, interpo-

5
lation, solutions of nonlinear equations f (x) = 0. Statistical topics: polynomial approximation,
curve fitting.
Further information on the elementary methods can be found in books on numerical methods
or books on numerical analysis. Dedicated text books can be found on each of the intensive
methods. Details of available books can be accessed through [Link] .

Need help understanding numerical methods?

(1) What is the use of numerical methods in real life application?

(2) Need a brief explanation of numerical methods

(3) Fixed Point Iteration, Linear Interpolation and Newton-Raphson Method, what are the
differences to their uses?

Best Answer

1. um, everywhere? From a cash machine, to calculating how much chemicals to put to
produce laundry detergent, to construction of buildings and bridges.

2. The ever-increasing advances in computer technology has enabled many in science and
engineering to apply numerical methods to simulate physical phenomena. Numerical
methods are often divided into elementary ones such as finding the root of an equation,
integrating a function or solving a linear system of equations to intensive ones like the
finite element method. Intensive methods are often needed for the solution of practical
problems and they often require the systematic application of a range of elementary
methods, often thousands or millions of times over. In the development of numerical
methods, simplifications need to be made to progress towards a solution: for example
general functions may need to be approximated by polynomials and computers cannot
generally represent numbers exactly anyway. As a result, numerical methods do not
usually give the exact answer to a given problem, or they can only tend towards a solution
getting closer and closer with each iteration. Numerical methods are generally only useful
when they are implemented on computer using a computer programming language .

3. Visit these sites: [Link]

[Link]
[Link]

Other answers

2. Numerical Methods refers to procedures to find approximate solutions when exact solu-
tions cannot be found in a straightforward manner.

3. Linear interpolation assumes that if two points on a graph are given, any point in between
them can be found by connecting the original two points by a straight line.
Newton Raphson is a method to find approximate solutions to an equation through an
iterative process where each calculated value is used as the starting point for the next
calculated value. NRM requires that you can evaluate the first derivative of that equation.

Page 6 of 133 D.W-Ddumba numerical analysis i

 1.1. WHAT IS NUMERICAL ANALYSIS?

1.1 What is Numerical Analysis?

- It is a way to do highly complicated mathematics 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:

- How to compute? This corresponds to algorithmic aspects;

- How accurate is it? That corresponds to error analysis aspects.

1.3 Advantages and Disadvantages of Numerical Analy-

sis:
Although numerical solutions are just an approximation, an estimate of an exact solution
(numerical values are ”near” solutions or often referred to as inexact solutions) due to

(i) truncation errors,

(ii) round off errors,

(iii) chop off errors,

Other than the truncation errors involved in numerical methods, numerical algorithms are
computationally involved, calling for powerful computing abilities and machines. The schemes
are usually multi steps and multi terms.
However, numerical techniques are widely applied in solving real world problem.

1.3.1 Hard or impossible classical solutions

It can obtain numerical answers of the problems that have very hard analytical techniques or
sometimes impossible to solve (no ”exact or classical” solution).

Example 1.3.1 Finding the classical solution to an integral

ˆ
2
ex dx

is quite very hard if not impossible.

Example 1.3.2 Solving a simple looking first order ordinary differential equation

dy
= x2 + y 2
dx
analytically is not easy at all.

D.W-Ddumba numerical analysis i Page 7 of 133

 1.4. IMPORTANT NOTES:

1.3.2 Tabulated data

Usually we only have tabulated data when faced with complex real world models, and do not
have an explicitly defined function.

Example 1.3.3 To compute the area under the curve defined by the table below
x 1 3 5 7
f (x) -3 5 21 45

will be analytically impossible since we do not know the explicit definition of f (x) that was
used to generate the table above. To compute
ˆ 7
f (x)dx
1

will not be possible.

The table above was actually generated by f (x) = x2 − 4 which we did not know prior, and
too hard to interpolate.

1.4 Important Notes:

- Numerical analysis solution is always numerical.

- Results from numerical analysis is an approximation.

1.5 Numerical Errors

When we get into the real world from an ideal world and finite to infinite, errors arise.

1.5.1 Sources of Errors:

- Mathematical problems involving quantities of infinite precision.

- Numerical methods bridge the precision gap by putting errors under firm control.

- Computer can only handle quantities of finite precision.

1.5.2 Types of Errors:

- Truncation error (finite speed and time)
Truncation error is a consequence of doing only a finite number of steps in a calculation
that would require an infinite number of steps to do exactly. A simple example of a
calculation that will be affected by truncation error is the evaluation of an infinite sum
using the NSum function. The computer certainly isn’t going to compute values for all
of the terms in an infinite sum. The terms that are left out lead to truncation error.
Truncation Error: The essence of any numerical method is that it is approximate this
usually occurs because of truncation, e.g., cos x ∼
2
= 1 − x2 or terminating an infinite
sequence of operations after a finite number have been performed.
It is not possible by numerical techniques alone to get an accurate estimate of the size of

Page 8 of 133 D.W-Ddumba numerical analysis i

 1.5. NUMERICAL ERRORS

the truncation error in a result. It is possible for any purely numerical algorithm, including
the algorithms used by numerical functions in Mathematica, to produce incorrect results,
and to do so without warning. The only way to be certain that results from functions
like NIntegrate and NDSolve are correct is to do some independent analysis, possibly
including detailed investigation of the algorithm that was used to do the calculation.
Such investigations are an important part of the field of numerical analysis.
For example, after using Taylor expansion

x x2 x3
ex = 1 + + + + ···
1! 2! 3!
x2 x4 x6 x8
cos x = 1 + + + + + ···
2! 4! 6! 8!
x3 x5 x7 x9
sin x = x − + − + + ···
3! 5! 7! 9!

You could realize that there are many terms you have truncated off in the expansion,
thats why · · ·

- Round off errors

Numbers can be stored only up to a fixed finite number of digits: Additional digits may be
rounded or chopped. Rounding error is sometimes characterized by ξmachine , the largest
(positive) number that the machine (computer) cannot distinguish between 1 and 1+
ξmachine .
Roundoff error, or representation error, is the error associated with the fact that the
computer keeps only a finite number of digits in calculations with inexact numbers. Since
it is not√possible (except in special cases) to represent all of the digits in numbers like 1/3
or π or 2, computers store only the first few digits in numerical approximations of these
numbers. In typical situations, the computer will store only the first 16 decimal digits or
the first 53 binary digits. The remaining digits are discarded. The discarded digits lead
to errors in the result. One of the more conspicuous symptoms of roundoff error is the
appearance of tiny non-zero numbers in results that would otherwise be zero.

Although the ability to reduce the effects of roundoff error by raising the precision of a
calculation is certainly very useful, it is far from a universal solution to all problems with
numerical error.
1.8625 to three decimal places, it becomes 1.863
1.8625 to two decimal places, it becomes 1.86
1.8625 to one decimal place, it become 1.9

- Human Errors Such as Computing tools/machines, Mathematical equation/model,

propagated error.

D.W-Ddumba numerical analysis i Page 9 of 133

1.5. NUMERICAL ERRORS

Page 10 of 133 D.W-Ddumba numerical analysis i

Chapter 2

Numerical Integration

There are two main reasons for you to need to do numerical integration: analytical integration
may be impossible or infeasible, or you may wish to integrate tabulated data rather than known
functions. In this section we outline the main approaches to numerical integration. Which is
preferable depends in part on the results required, and in part on the function or data to be
integrated.
That is
This will be useful when we cannot find an elementary antiderivative for f (x) or if the function
is defined using data obtained from some experiment.
Numerical integration is the numerical approximation of the integral of a function. For a
function of one variable, it amounts to finding the area under the graph of the function. That
is finding I where ˆ b
I= f (x) dx
a
Methods generally replace the integral by a weighted sum of n weights and n function evalua-
tions, so that
n
X
I= Wi f (xi ) dx
i=1

For a function of two variables it is equivalent to finding an approximation to the volume under
the surface. Numerical integration is often also referred to as quadrature or sometimes cubature
for functions of two or more variables. Returning to the one variable case, numerical integration
involves finding the approximation to an integral of a function f(x) through its evaluation at
a set of discrete points. There are two distinct approaches to this. Firstly methods like the
trapezium rule or Simpson’s rule determine the integral through evaluating f(x) at regularly
spaced points. These are generally referred to as Newton-Cotes formulae.
Alternative methods termed Gaussian Quadrature methods have arisen that select irregularly-
placed evaluation points, chosen to determine the integral as accurately as possible with a given
set of points.
Gaussian Quadrature methods are important as they often lead to very efficient methods. In
numerical integration the efficiency of the method relates to the accuracy obtained with respect
to the number of evaluations of the function f (x). In intensive methods such as the boundary
element method integrations may need to be performed millions of times so the efficiency of
the methods needs to be considered sometimes.
In general, care must be taken to match the numerical integration method to the expected

11
 2.1. MANUAL METHOD

nature of the function f (x). For example it may be known that f (x) is regular. On the other
hand f (x) may be singular or oscillatory and will then need special treatment. Often a special
method called a product integration method can be developed for the integration of functions
of the form f (x) = w(x)g(x) where w(x) is a pre-set function and the function and g(x) is
known to be a relatively nice function.
There are books devoted to numerical integration. Numerical integration is a basic numerical
method and the subject is generally covered in books on numerical methods or numerical
analysis (see below). Numerical methods for carrying out numerical integration can often be
easily programmed and can also be found in general numerical libraries.

2.1 Manual Method

If you were to perform the integration by hand, one approach is to superimpose a grid on
a graph of the function to be integrated, and simply count the squares, counting only those
covered by 50% or more of the function. Provided the grid is sufficiently fine, a reasonably
accurate estimate may be obtained.

Page 12 of 133 D.W-Ddumba numerical analysis i

 2.2. TRAPEZOIDAL/TRAPEZIUM RULE

2.2 Trapezoidal/Trapezium rule

2.2.1 Composite Trapezoidal Rule
To derive the rule, we use the following figure

Figure 2.1: Illustration of composite Trapezoidal rule

(b−a)
We divide [a, b] into n equal intervals of width h = n
Let
ˆ b ˆ x1
I= f (x)dx = f (x)dx+
a a
ˆ x2 ˆ xi+1 ˆ xn
f (x)dx + . . . + f (x)dx + . . . + f (x)dx.
x1 xi xn−1

Applying the trapezium equation , we get,

h
I= [f0 + 2f1 + 2f2 + . . . + 2fn−1 + fn ] + Etrunc ,
2
where

h3 00
Etrunc = − [f (c1 ) + f 00 (c2 ) + . . . + f 00 (cn )]
12
M h3
|Etrunc | ≤
12n2
such that the second derivative f 00 is continuous on [a, b] and that |f 00 (x)| < M for all x in [a, b].

Example 2.2.1 Approximate the integral

ˆ 1
I= x2 dx
0

using the composite trapezoidal rule with step length h = 0.2

Solution
Since

D.W-Ddumba numerical analysis i Page 13 of 133

 2.2. TRAPEZOIDAL/TRAPEZIUM RULE

ˆ 1
h
I= f (x)dx ' {f0 + 2f1 + 2f2 + 2f3 + 2f4 + f5 }
0 2
Partitioning the interval [0, 1], we get,

f0 = f (x0 ) = f (0) = 02 = 0
f1 = f (x1 ) = f (0.2) = (0.2)2 = 0.04
f2 = f (x2 ) = f (0.4) = (0.4)2 = 0.16
f3 = f (x3 ) = f (0.6) = (0.6)2 = 0.36
f4 = f (x4 ) = f (0.8) = (0.8)2 = 0.64
f5 = f (x5 ) = f (1) = 12 = 1

Therefore
0.2
I' {0 + 2[0.04 + 0.16 + 0.36 + 0.64] + 1}
2
= 0.1{0 + 2[1.2] + 1} = 0.340.
Absolute error committed is
0.340 − 0.333 ' 0.00667
We have noted that the error obtained is much smaller than that obtained with the pure
trapezoidal rule in the previous example. In fact the smaller the error, i.e the better the ap-
proximation.

The truncation error is

M h3 2(0.2)3
|Etrunc | ≤ 2
= 2
since f 00 (x) = 2, |f 00 (x)| < 2 on [0, 1]
12n (12)5

Example 2.2.2 Compute the numerical value of the definite the integral
ˆ 2
x2
I= e− 2 dx
−2

by the composite trapezoidal rule with h = 1.0. The exact value of the integral I to 4 decimal
places is 2.3925.
Solution
Since h = 1.0, then ,

−(1)2
 
1.0 − (−2)2 (−1)2
− 2
(0)2
− 2 + 12 e− 2
−(2)2
− 2
I = e 2 + 2e + 2e +e
2
= 0.5(0.13534 + 2 × 0.60653 + 2 × 1.00000 + 2 × 0.60653 + 0.73534)
= 2.3484

with error of 0.0441.

Example 2.2.3 Estimate ˆ 2

1
dx
1 x+1
using

Page 14 of 133 D.W-Ddumba numerical analysis i

 2.2. TRAPEZOIDAL/TRAPEZIUM RULE

(i) Trapezium rule with one subinterval

 
h 1 1 1
[f0 + f1 ] = + = 0.416
2 2 2 3

(ii) Composite trapezium rule scheme with h = 0.2

h
' {f (1.0) + 2f (1.2) + 2f (1.4) + 2f (1.6) + 2f (1.8) + f (2.0)}
2
0.2
' {0.5 + 2(0.45) + 2(0.42) + 2(0.38) + 2(0.36) + 0.33}
2
' 0.4057

Example 2.2.4 Evaluate

ˆ 3
(2x + 3) dx
0

by the Trapezium rule with four intervals (5 ordinates).

h
' {f (0) + 2f (0.75) + 2f (1.5) + 2f (2.25) + f (3.0)}
2
0.75
' {3 + 2(4.5) + 2(6) + 2(7.5) + 9}
2
' 18

Example 2.2.5 Evaluating

ˆ 3
sin x dx
1

by the Trapezium rule with 100 points, it gives the answer as 1.5302437923
But can you think of a way of programme this easily??

Example 2.2.6 Evaluating

ˆ 3
sin x2 dx
0

by the Trapezium rule

n Sum of areas of trapezoids

4 0.43358
8 0.70404
16 0.75723
32 0.76954
64 0.77256
128 0.77331
256 0.77350
512 0.77355
1024 0.77356
2048 0.77356

0.77356 appears to be a reasonable estimate of our integral.

D.W-Ddumba numerical analysis i Page 15 of 133

 2.2. TRAPEZOIDAL/TRAPEZIUM RULE

Example 2.2.7 Using the Trapezoidal algorithm, approximate

ˆ 1
√
2x + 1 dx
0

with h = 0.25

√ x 0 0.25 0.5 0.75 1

2x + 1 1 1.22 1.41 1.58 1.73
1
A≈ [1 + 1.73 + 2(1.22 + 1.41 + 1.58)] = 1.39
2(0.25)
(i) Repeat the problem with h = 0.2

(ii) Repeat the problem with h = 0.1

(iii) Compute the analytical, classical, or exact solution to the definite integral above

(iv) State the effect of h on the numerical solution

Exercise 2.1
´1
1. Compute the approximate value of 0 (x2 + 1)−1 dx by using the trapezoidal rule with
ten subintervals. Then compare with the actual value of the integral. Determine the
truncation error bound and compare with the actual error.
´5
2. If the trapezoidal rule is used to compute 2 sin xdx with h = 0.01, obtain an upper
bound error.
´2
3. How large must n be if the trapezoidal rule is to estimate 0 exp{−x2 }dx with an error
not exceeding 10−6 ?
´1 2
4. Consider the integral 0 sin(π x2 )dx. Suppose that we wish to integrate numerically with
error < 10−5 . What interval width h is needed if the trapezoidal rule is to be used?
´3
5. Approximate 1 x1 dx by the trapezoidal rule with an error of utmost 0.1.

Page 16 of 133 D.W-Ddumba numerical analysis i

 2.3. SIMPSON’S RULE

2.3 Simpson’s rule

We can obtain the Simpson’s rule by various ways. One of the most popular ways is from
Lagrange’s quadratic interpolation polynomial. The Simpson’s rule approximates the area
under the curve y = f (x) from x0 to x2 by the area under a parabolic curve. The figure below
illustrates parabolic fitting on to the curve y = f (x) from x0 to x2 .
Interpolating f (x) by a Lagrange polynomial of degree 2 i.e. P2 (x) then

f (x) = P2 (x) + Etrunc (x)

So ˆ ˆ ˆ
x2 x2 x2
f (x)dx = P2 (x)dx + Etrunc (x)dx (2.1)
x0 x0 x0

Results
Summing up all the three cases, equation (2.1) becomes
ˆ x2
h
P2 (x)dx = [f0 + 4f1 + f2 ] + Etrunc (x) Thus
x0 3
ˆ x2
h
P2 (x)dx = [f0 + 4f1 + f2 ] (2.2)
x0 3
Relation equation (2.2) is the Simpson’s rule for approximating the integral. The integral for
the error in equation (2.1), becomes
ˆ x2 ˆ x2
(x − x0 )(x − x1 )(x − x2 ) (4)
Etrunc (x)dx = f (c(x))dx
x0 x0 3!

This can be shown (with difficulty!) to be

M h5
(2.3)
180n4
That the fourth derivative f 4 is continuous on [a, b] and that |f 4 (x)| < M for all x in [a, b].

´Example
1 2
2.3.1 Use Simpson’s method (or with n = 1) to approximate the integral I =
0
x dx
Solution ˆ 1
h
Since I = f (x )dx ' [f (x0 ) + 4f (x1 ) + f (x2 )]
0 3

1 (1 − 0 ) 1 (b − a)
But x0 = 0 , x1 = , x2 = 1 , h = = = and n = 1.
2 (2 )(1 ) 2 2n
Therefore
1 1 1 1 1
I ' [f (0) + 4f ( ) + f (1)] = [02 + 4.( )2 + 12 ] = = 0.33
6 2 6 2 3
But the exact value of the integral is 13 = 0.33. It should not surprise you that the Simpson”
rule has generated the exact value of the integral. In fact the general result is that for f (x) a
polynomial of degree two or less, the Simpson” rule will always generate the exact value of the
integral. This will later be stated as a theorem.

D.W-Ddumba numerical analysis i Page 17 of 133

 2.3. SIMPSON’S RULE

Figure 2.2: Illustration of the composite Simpson’s rule

2.3.1 Composite Simpson’s rule

Let’s consider the Figure (2.2),
(b−a)
We divide the interval [a, b] into 2n equal intervals of width h = 2n
.
Thus the integral,
ˆ b
I= f (x)dx
a

becomes,
ˆ b ˆ x2 ˆ x4 ˆ 2n
I= f (x)dx = f (x)dx + f (x)dx + . . . + f (x)dx
a a x2 x2n−2

h
= {[f0 + 4f1 + f2 ] + [f2 + 4f3 + f4 ]+
3
[f4 + 4f5 + f6 ] + . . . + [f2n−2 + 4f2n−1 + f2n ]}
h
= [f0 + 4(f1 + f3 + . . . + f2n−1 ) + 2(f2 + f4 + . . . + f2n−2 ) + fn ] + Etrunc
3
Where the truncation error,

h5 (4)
Etrunc = − [f (c1 ) + f (4) (c2 ) + . . . + f (4) (cn )]
90
h5 (4)
=− f (c2 ) × n
90
(b − a)h4 (4)
=− f (cf )
180
where a ≤ cf ≤ b.

Example 2.3.2 Use the Simpson’s rule to compute the integral

ˆ 2
x2
I= e− 2 dx
−2

Page 18 of 133 D.W-Ddumba numerical analysis i

 2.3. SIMPSON’S RULE

Using step size h = 1.0. Recall the exact value of I to 4 decimal places is 2.3925.
Solution
Using the Simpson’s rule, we have,
 
1.0 − (−2)2 −
(−1)2
−
(0)2
−
(1)2
−
(−2)2
I' e 2 + 4e 2 + 2e 2 + 4e 2 + e 2 = 2.3743
3
The error committed is 0.0182.
Note 2.3.1 We note that is error is much smaller than that obtained when using Trapezoidal
rule in the previous lecture though same step size is used.
Note 2.3.2 Simpson’s rule only work for even number of subintervals.
Note 2.3.3 For f (x) a trignometric function, we use (Your calculator should be in) Radians.

Example 2.3.3 Determine the fairly accurate solution to the definite integral
ˆ 2
1
dx
1 x

using Simpson’s algorithm with n = 6

ˆ 2      
1 1 1 3 3 6 2 6
dx = 1+ +2 + +4 + +
1 x 18 2 4 5 7 3 4
14411
= ≈ 0.6931697
20790

Example 2.3.4 Approximate the definite integral

ˆ 3
dx
dx
2 x+1
using Simpson’s scheme with n = 4
ˆ 3
dx 0.25
dx = [(f (2) + f (3)) + 2 (f (2.5)) + 4 (f (2.25) + f (2.75))]
2 x+1 3
0.25
= [(0.33333 + 0.25) + 2 (0.2857142) + 4 (0.3076923 + 0.2666667)]
3
≈ 0.2876831

(i) Compute its exact value and hence the absolute error committed by the Simpson’s tech-
nique with n = 4.
ˆ 3
dx
dx = [ln(x + 1)]32 = 0.287682
2 x + 1
⇒ error = 0.287682 − 0.2876831
≈ 0.00036

(ii) Repeat the problem above using h = 0.1

(iii) Comment on the effect of h.

D.W-Ddumba numerical analysis i Page 19 of 133

 2.3. SIMPSON’S RULE

(iv) Compute the inexact value of the problem using the Trapezoidal rule.

(v) Which of the two numerical techniques (the Trapezoidal or the Simpson) is more accurate
and superior?

Example 2.3.5 It is required to obtain

ˆ 2
2
ex dx
0

exact to 4 decimal places. What should h be for Simpson’s rule.

Solution
Since the error term is
(b − a) 4 (4)
− h f (cf )
180
But
2
f (x) = ex
therefore
2
f 0 (x) = 2xex
2 2 2 2
f 00 (x) = 2(ex + 2xex ) = 2ex (1 + 2x2 ) = ex (2 + 4x2 )
2 2
therefore, f 000 (x) = 2xex (2 + 4x2 ) + 8xex
2 2
= ex (4x + 8x3 + 8x) = ex (12x + 8x3 )
2
and f (iv) (x) = 8ex (2x4 + 5x2 + 1) < 424e4
2h4
So h! · 424e4 < (0.5)10−4 = 0.057
180
Say choose h = 0.05.

Exercise 2.2 Compute an approximate value of

ˆ 1
(x2 + 1)−1
0

Using Composite Simpson’s rule with

(i) h = 0.1,

(ii) h = 0.5.

Compare with the actual value of the integral in each case. Next, determine the truncation
error bound and compare with the actual error.

Exercise 2.3 If the Simpson’s rule is used to compute

ˆ 5
sin xdx
2

with h = 0.75, obtain an upper bound on the error.

Page 20 of 133 D.W-Ddumba numerical analysis i

 2.3. SIMPSON’S RULE

Exercise 2.4 Establish the composite Simpson’s rule over (n − 1) even subintervals

ˆ b 2
(n−1)
2
(n−3)

h X X
f (x)dx = [(f (a) + f (b)) + 4 f (a + (2i − 1)h) + 2 f (a + 2ih)] + Etrunc
a 3 i=1 i=1

(b−a)
where, h = (n−1)
and
(b − a) 4 (4)
Etrunc = − h f (c)
180
for some c ∈ [a, b]

Exercise 2.5 Consider the integral

ˆ 1
πx2
sin( )dx.
0 2

Suppose that we wish to integrate numerically with error < 10−5 . What interval width h is
needed if the Simpson’s rule is to be used?

Exercise 2.6 Compute ˆ 2

(x3 + 1)dx
0
1
by using h = 4
and compare with the exact value of the integral.

Example 2.3.6 Using Maple; use Simpson’s Rule with n = 100 to approximate the integral
ˆ 1
1
2
dx
0 1+x

> with(student);
1
> simpson( 1+x 2 , x = 0..1, 100);

50
! 49
!
1 1 X 1 1 X 1
+ +
200 75 1 1 2 150 1+ 1 2
i


i=1 1+ 50
i − 100 i=1 2500

> evalf(”);
.7853981634
> evalf(Pi/4);
.7853981635
When trying for the trapezium rule, instead of ”simpson”, you replace it with ”trapezoid”
in step II

2.3.2 Program-FORTRAN (an alternative to MAPLE)

Note that this program is written for clarity rather than speed. The number of function
evaluations actually computed may be approximately halved for the Trapezium rule and reduced
by one third for Simpson’s rule if the compound formulations are used. Note also that this
example is included for illustrative purposes only. No knowledge of Fortran or any other
programming language is required in this course
REAL*8 FUNCTION TrapeziumRule(x0,x1,nx)
C=====parameters

D.W-Ddumba numerical analysis i Page 21 of 133

 2.3. SIMPSON’S RULE

INTEGER*4 nx
REAL*8 x0,x1
C=====functions
REAL*8 f
C=====local variables
INTEGER*4 i
REAL*8 dx,xa,xb,fa,fb,Sum
dx = (x1 - x0)/DFLOAT(nx)
Sum = 0.0
DO i=0,nx-1
xa = x0 + DFLOAT(i)*dx
xb = x0 + DFLOAT(i+1)*dx
fa = f(xa)
fb = f(xb)
Sum = Sum + fa + fb
ENDDO
Sum = Sum * dx / 2.0
TrapeziumRule = Sum
RETURN
END

Assignment 2.3.1 Evaluate

ˆ 2
2−1 1
x2 cos xdx, f (x) = x2 cos x with h = = . [a, b] = [1, 2]
1 6 6

using (i)Trapezium rule, (ii) Simpson’s rule

Example 2.3.7

(a) (i) State one reason to justify numerical techniques of integration.

The analytic integral can be very complicated, and some times impossible. For
´ 3 −x 2
example 0 e dx, but also The function f (x) might not be known, but can be
defined on some discrete points.
(ii) One of the commonly used numerical methods of integrations is the Simpson’s rule,
state the rules.
IS = h3 [f0 + 4(f1 + f3 + . . . + f2n−1 ) + 2(f2 + f4 + . . . + f2n−2 ) + fn ] + Etrunc

(b) Derive Trapezium’s rule for the integration of a function f (x) between a and b,

h
IT = [f0 + 2f1 + 2f2 + . . . + 2fn−1 + fn ] + Etrunc
2
State its truncation error.
(b−a)
We divide [a, b] into n equal intervals of width h = n
Let
ˆ b ˆ x1
I= f (x)dx = f (x)dx+
a a
ˆ x2 ˆ xi+1 ˆ xn
f (x)dx + . . . + f (x)dx + . . . + f (x)dx.
x1 xi xn−1

Page 22 of 133 D.W-Ddumba numerical analysis i

 2.3. SIMPSON’S RULE

Figure 2.3: Illusttration of Trapezoidal rule

Applying the trapezium equation , we get,

h
I= [f0 + 2f1 + 2f2 + . . . + 2fn−1 + fn ] + Etrunc ,
2
where

h3 00
Etrunc = − [f (c1 ) + f 00 (c2 ) + . . . + f 00 (cn )]
12
M h3
|Etrunc | ≤
12n2
such that the second derivative f is continuous on [a, b] and that |f 00 (x)| < M for all x
00

in [a, b].
(c) Numerically approximate ˆ 2  √ 
2 + cos(2 x ) dx
0
by using the Trapezium rule with
1
(i) n = 4 ⇒ h = 2

h
I ≈ [f0 + 2f1 + 2f2 + . . . + 2fn−1 + fn ]
2
1 h √   √  h √   √   √ ii
2 + cos(2 0) + 2 + cos(2 2) + 2 2 + cos(2 0.5) + 2 + cos(2 1) + 2 + cos(2 1.5)
4

1h √ h √ √ ii
≈ 3 + 2 + cos(2 2) + 2 (2 + cos 2) + (2 + cos 2) + (2 + cos 6)
4
1h √ h √ √ ii
≈ 5 + cos(2 2) + 2 (2 + cos 2) + (2 + cos 2) + (2 + cos 6)
4
1
= [13.98841913]
4
≈ 3.4971
1
(ii) n = 8 ⇒ h = 4

h
I = [f0 + 2f1 + 2f2 + . . . + 2fn−1 + fn ]
2 √   √   √   √ 
= 2 + cos(2 0) + 2 + cos(2 2) + 2 2 + cos(2 0.25) + 2 2 + 2 cos(2 0.5)
 √   √   √   √ 
+ 2 2 + cos(2 0.75) + 2 2 + cos(2 1) + 2 2 + cos(2 1.25) + 2 2 + cos(2 1.5) +
≈ 3.46928

D.W-Ddumba numerical analysis i Page 23 of 133

 2.3. SIMPSON’S RULE

(d) Considering your results in part (c) above, state any two reasons on how to reduce the
errors in numerical integration.

(e) (i) Solve the integral in part (c) above using the Simpson’s rule with n = 4.
h
I = [f0 + 4(f1 + f3 + . . . + f2n−1 ) + 2(f2 + f4 + . . . + f2n−2 ) + fn ] + Etrunc
3
1 h √   √  h √   √ i
≈ 2 + cos(2 0) + 2 + cos(2 2) + 2 2 + cos(2 0.5) + 2 + cos(2 1.5)
6 
√ i
+ 4 2 + cos(2 1)
1h √ h √ √ i i
≈ 3 + 2 + cos(2 2) + 2 (2 + cos 2) + (2 + cos 6) + 4(2 + cos 2)
6
1h √ h √ √ ii
≈ 5 + cos(2 2) + 2 (2 + cos 2) + (2 + cos 6) + 4(2 + cos 2)
6
≈ 3.46008250981

(ii) Show that using the Simpson’s rule with n = 8, the integral in part (c) above is
I ≈ 3.460002979.
h
I = [f0 + 4(f1 + f3 + . . . + f2n−1 ) + 2(f2 + f4 + . . . + f2n−2 ) + fn ] + Etrunc
3
1 h √   √  n √   √ 
≈ 2 + cos(2 0) + 2 + cos(2 2) + 2 2 + cos(2 0.5) + 2 + cos(2 1)
12 
√ o n √   √ 
+ 2 + cos(2 1.5) + 4 2 + cos(2 0.25) + 2 + cos(2 0.75)
 √   √ oi
+ 2 + cos(2 1.25) + 2 + cos(2 1.75)
≈ 3.460002979

(f) Which of the two techniques of integration is superior.

It is the Simpson’s rule that is more superior.

Exercise 2.7 Compute ˆ 2

1
dx
1 x
with n = 10 using
(i) Trapezium rule
0.6938

(ii) Simpson’s rule

0.69315
´1 4
Exercise 2.8 Use the simpson rule to solve 0 1+x 2 dx with n = 10

3.14159

´2 1
Exercise 2.9 Show that the error (absolute error ) in the trapezium rule to compute 1 x
dx is
2
' 0.0017
600

Page 24 of 133 D.W-Ddumba numerical analysis i

 2.4. MID-POINT RULE

2.4 Mid-Point Rule

If f is continous on [a,b] and f (x) ≥ 0∀x ∈ (a, b), we partition [a, b] into n subintervals of equal
lengths.
If mi is the midpoint of the ith interval then
ˆ b n
X
f (x)dx = (∆x)f (mi )
a i=1

´2 1
Example 2.4.1 Find the integral using the midpoint rule. 1 x
dx with n = 10
2−1 1
∆x = 10
= 10
Sometimes called step width.
1 1
f (x) = x
⇒ f (mi ) = mi

ˆ 2 n  
1 X 1 1 1 1 1 1 1
⇒ dx = ∆x f (mi ) = + + + + + ........... + = 0.69264
1 x i=1
10 1.15 1.25 1.35 1.45 1.55 1.95

Exercise 2.10 Use the Riemann and check the answer.

2.4.1 Error in midpoint

If ||f ”(x)|| < m ∀ x ∈ [a, b] then
M (b − a)3
ξm ≤
24n2
Exercise 2.11 Using n = 4 and all three rules to approximate the value of the following
integral ˆ 2
2
ex dx
0

(i) Mid-Point rule 14.48561253

(ii) Trapezium rule 20.64455905

(iii) Simpson’s rule 17.35362645

Maple gives the exact solution as

ˆ 2
2
ex dx = 16.45262776
0

D.W-Ddumba numerical analysis i Page 25 of 133

2.4. MID-POINT RULE

Page 26 of 133 D.W-Ddumba numerical analysis i

Chapter 3

Solution to Non-Linear Equations

Motivation
Given a second degree polynomial, we can accurately compute it’s root, that is
√
2 −b ± b2 − 4ac
ax + bx + c = 0 ⇒ x =
2a
But what if ax5 + bx4 + cx3 + dx2 + ex + f = 0 ⇒ x =?, sin x + x = 0 ⇒ x =?

Problem Description
• Given a non-linear equation f (x) = 0, find a x? such that f (x? ) = 0. Thus, x? is a root
of f (x) = 0.

• Galois theory in math tells us that only polynomials of degree ≤ 4 can be solved with
close forms using +, −, ×, ÷ and taking roots.

• General non-linear equations can be solved with iterative methods.

• Basically, we try to guess the location of a root, and approximate it iteratively.

• Unfortunately, this process can go wrong, leading to another root or even diverge.

Methods to be discussed
• There are two types of methods, bracketing and open. The bracketing methods require
an interval that is known to contain a root, while the open method does not.

• Commonly seen bracketing methods include the bisection method and the regula falsi
method, and the open methods are Newtons method, the secant method, and the successive
substitution.

27
 3.1. NEWTON RAPHSON’S METHOD

3.1 Newton Raphson’s method

3.1.1 Derivation of the Newton’s method
If we are given a non-linear equation f (x) = 0 and we are to apply the Newton Raphson’s
method, we linear approximate the graph of y = f (x) by a straight line passing through the
point (x0 , f0 ) and tangential to the graph of y = f (x). Take the slope of this line to be p.
Geometrically this is given in figure below.

Figure 3.1: Geometrical representation of a tangent to a curve at a point

The equation of the line with slope p and passing through the point (x0 , f0 ) is

y − f0
=p (3.1)
x − x0

However, we know that p is the slope of the tangent to y = f (x) at (x0 , f0 ). This is given by;

p = f 0 (x0 ) = f00 (3.2)

Substituting equation (3.2) in equation (3.1) we get

y − f0
= f00
x − x0

y − f0 = (x − x0 )f00 (3.3)
From figure 3.1, line of equation 13.3 cuts the x-axis at the point (x1 , 0) i.e when x = x1 and
y = 0.
Substituting in equation (3.3) we get,

0 − f0 = (x1 − x0 )f00

Making x1 the subject, we get,

f0
x1 = x0 − or
f00
f (x0 )
x1 = x0 − (3.4)
f 0 x0

Page 28 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

equation (3.4) is actually the Newton’s method for obtaining the next iterate x1 from the
previous iterate x0 . The equation (3.4) is generalized and written;

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

since the linear approximation of the curve is done at each of the iterates xn , xn+1 , xn+1 , . . . as
reflected in figure above.

D.W-Ddumba numerical analysis i Page 29 of 133

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.1 Use Newton Raphson’s method to find the root of

x2 − 3 = 0 on [1, 2]

f (xn ) = x2n − 3
therefore f 0 (xn ) = 2xn
f (xn )
But Raphson’s formula is xn+1 = xn − 0
f (xn )

Substituting in the Raphson’s formula we get

(x2n − 3)
xn+1 = xn −
2xn
2xn − x2n + 3
2
=
2x
 2 n
xn + 3
=
2xn

Taking the initial guess/approximation as x0 = 2, but you could also consider x0 = 1 you come
up with the same answer.

x20 + 3 22 + 3
⇒ x1 = = = 1.7500000
2x0 2(2)
2
x +3 (1.75)2 + 3
x2 = 1 = = 1.7321000
2x1 2(1.75)
2
x +3 (1.7321)2 + 3
x3 = 2 = = 1.7320508
2x2 2(1.7321)
x2 + 3 (1.7320508)2 + 3
x4 = 3 = = 1.7320508
2x3 2(1.7320508)

Thus the root is 1.7320508

Page 30 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.2 Use Newton-Raphson Method to find the only real root of the equation
x3 − x − 1 = 0
correct to 9 decimal places.
Since f (1) = −1 and f (2) = 5, the function has a root in the interval [1, 2] since the function
changes sign between [1, 2]. Let us make an initial guess x0 = 1.5.
f (x) = x3 − x − 1 ⇒ f (xn ) = x3n − xn − 1
f 0 (x) = 3x2 − 1 ⇒ f 0 (xn ) = 3x2n − 1

f (xn )
xn+1 = xn −
f 0 (xn )
(x3 − xn − 1)
xn+1 = xn − n 2
(3xn − 1)
3
2xn + 1
xn+1 =
3x2n − 1
2x30 + 1 2(1.5)3 + 1
x1 = = = 1.347826087
3x20 − 1 3(1.5)2 − 1
2x31 + 1 2(1.347826087)3 + 1
x2 = = = 1.325200399
3x21 − 1 3(1.347826087)2 − 1
2x32 + 1 2(1.325200399)3 + 1
x3 = = = 1.324718174
3x22 − 1 3(1.325200399)2 − 1
2x33 + 1 2(1.324718174)3 + 1
x4 = = = 1.324717957
3x23 − 1 3(1.324718174)2 − 1
2x34 + 1 2(1.324717957)3 + 1
x5 = = = 1.324717957
3x24 − 1 3(1.324717957)2 − 1
The approximated root of
x3 − x − 1 = 0
to 9 decimal places is 1.324717957
Remark 3.1.1 The more the decimal places the better the approximation.

Example 3.1.3 Repeat Example 3.1.2

x3 − x − 1 = 0
using 4 decimal places.
f (xn ) (x3n − xn − 1) 2x3n + 1
xn+1 = xn − = x n − =
f 0 (xn ) (3x2n − 1) 3x2n − 1
2x30 + 1 2(1.5)3 + 1
x1 = = = 1.3478
3x20 − 1 3(1.5)2 − 1
2x31 + 1 2(1.3478)3 + 1
x2 = = = 1.3252
3x21 − 1 3(1.3478)2 − 1
2x32 + 1 2(1.3252)3 + 1
x3 = = = 1.3247
3x22 − 1 3(1.3252)2 − 1
2x33 + 1 2(1.3247)3 + 1
x4 = = = 1.3247
3x23 − 1 3(1.3247)2 − 1

D.W-Ddumba numerical analysis i Page 31 of 133

 3.1. NEWTON RAPHSON’S METHOD

The approximated root of

x3 − x − 1 = 0
to 4 decimal places is 1.3247

Example 3.1.4 Repeat Example 3.1.2

x3 − x − 1 = 0

using 1 decimal places.

f (xn ) (x3n − xn − 1) 2x3n + 1

xn+1 = xn − = x n − =
f 0 (xn ) (3x2n − 1) 3x2n − 1
2x30 + 1 2(1.5)3 + 1
x1 = = = 1.3
3x20 − 1 3(1.5)2 − 1
2x31 + 1 2(1.3)3 + 1
x2 = = = 1.3
3x21 − 1 3(1.3)2 − 1

The approximated root of

x3 − x − 1 = 0
to 1 decimal places is 1.3

Remark 3.1.2 The fewer the decimal places the faster the iteration scheme converges.

Page 32 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.5 Find the root of the function

y = e−x − x

in the vicinity of x = 0.5 correct to 4 decimal places.

f (x) = e−x − x ⇒ f (xn ) = e−xn − xn

f 0 (x) = −e−x − 1 ⇒ f 0 (xn ) = −e−xn − 1

f (xn )
xn+1 = xn −
f 0 (xn )
(e−xn − xn ) −xn e−xn − xn − e−xn + xn (xn + 1) e−xn
xn+1 = xn − = =
(−e−xn − 1) (−e−xn − 1) 1 + e−xn
(x0 + 1) e−x0 ([0.5] + 1) e−0.5
x1 = = = 0.5663
1 + e−x0 1 + e−0.5
(x1 + 1) e−x1 ([0.5663] + 1) e−0.5663
x2 = = = 0.5671
1 + e−x1 1 + e−0.5663
(x2 + 1) e−x2 ([0.5671] + 1) e−0.5671
x3 = = = 0.5671
1 + e−x2 1 + e−0.5671
The approximated root of
e−x − x = 0
to 4 decimal places is 0.5671

Example 3.1.6 Pitfall - no real root

Newtons method will fail (the sequence {xn } does not converge) on solving

f (x) = x2 − 4x + 5 = 0

since f does not have any real root.

D.W-Ddumba numerical analysis i Page 33 of 133

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.7 Find the root of

f (x) = 3x + sin x − ex = 0

in the interval [0, 1] numerically by the famous Newton Raphson’s method

Solution
Since f (0) = −1 < 0 and f (1) = 3 + sin(1) − e > 0, so there is a real root in [0, 1]. Using

f (xn ) (3xn + sin xn − exn )

xn+1 = xn − = x n −
f 0 (xn ) (3 + cos xn − exn )
3xn + xn cos xn − xn exn − 3xn − sin xn + exn
=
3 + cos xn − exn
xn [cos xn − e ] − sin xn + exn
xn
=
3 + cos xn − exn
to have the iterations
x0 [cos x0 − ex0 ] − sin x0 + ex0 0 [cos 0 − e0 ] − sin 0 + e0
x1 = = = 0.33333
3 + cos x0 − ex0 3 + cos 0 − e0
x1 [cos x1 − ex1 ] − sin x1 + ex1
x2 =
3 + cos x1 − ex1
0.33333 [cos 0.33333 − e0.33333 ] − sin 0.33333 + e0.33333
= = 0.36017
3 + cos 0.33333 − e0.33333
x2 [cos x2 − ex2 ] − sin x2 + ex2
x3 =
3 + cos x2 − ex2
0.36017 [cos 0.36017 − e0.36017 ] − sin 0.36017 + e0.36017
= = 0.36042
3 + cos 0.36017 − e0.36017
x3 [cos x3 − ex3 ] − sin x3 + ex3
x4 =
3 + cos x3 − ex3
0.36042 [cos 0.36042 − e0.36042 ] − sin 0.36042 + e0.36042
= = 0.36042
3 + cos 0.36042 − e0.36042
The approximated root of
3x + sin x − ex = 0
to 5 decimal places is 0.36042

Page 34 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.8 Pitfall - Alternating or oscillating solutions

Use Newtons Method to find the only real root of the equation

f (x) = ex − 2x = 0

with an initial guess of x0 = 1.

f (x) = ex − 2x ⇒ f (xn ) = exn − 2xn

f 0 (x) = ex − 2 ⇒ f 0 (xn ) = exn − 2

f (xn ) exn − 2xn exn (xn − 1)

xn+1 = xn − = x n − =
f 0 (xn ) exn − 2 e xn − 2
ex0 (x0 − 1) e1 (1 − 1)
x1 = = =0
ex0 − 2 e1 − 2
ex1 (x1 − 1) e0 (0 − 1)
x2 = = 0 =1
ex1 − 2 e −2
ex2 (x2 − 1) e1 (1 − 1)
x3 = = =0
ex2 − 2 e1 − 2
ex3 (x3 − 1) e0 (0 − 1)
x4 = = =1
ex3 − 2 e0 − 2
ex4 (x4 − 1) e1 (1 − 1)
x5 = = 1 =0
ex4 − 2 e −2
ex5 (x5 − 1) e0 (0 − 1)
x6 = = =1
ex5 − 2 e0 − 2
..
.
..
.

Its an alternating solutions, the solutions oscillate near the local maxima or local minima.
In other words, Newton’s Method fails to produce a solution. Why is this? Because there is no
solution to be found!
Mathematicians are often very happy when, after a great deal of work, they are just able to
say that a solution to a problem exists. This is because once they know it exists, there might
be some nice method, such as Newton’s Method, to actually compute the solution.

D.W-Ddumba numerical analysis i Page 35 of 133

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.9 Pitfall - Diverging sequence

When Newtons method is applied to

f (x) = xe−x

with x0 = 2, it produces

n xn f (xn )

0 2.00000 2.70671 × 10−1

1 4.00000 7.32626 × 10−2

2 5.33333 2.57491 × 10−3

3 6.56410 9.25597 × 10−4

.. .. ..
. . .

20 24.96488 3.59105 × 10−10

21 26.00660 1.31995 × 10−10

22 27.04659 4.85206 × 10−11

◦ The sequence {xn } diverges to ∞ slowly;

◦ However, f (xn ) goes to zero rapidly as xn gets larger in a finite precision environment,
and could be mistaken as a zero of f .

Definition 3.1.1 Convergence of Newtons Method

If the iterations are getting closer and closer to the correct answer the method is said to converge.
That is if |xn+1 − xn | → 0
However, Newtons method will not converge if

? If f 0 (xn ) = 0 for some n

? If lim xn does not exist

n→∞

Page 36 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.10 Pitfall - Interval size

If the numerical When Newtons method applied to

f (x) = x3 − x − 3

The algorithm given by

f (xn ) x3n − xn − 3
xn+1 = xn − = x n −
f 0 (xn ) 3x2n − 1
2x3n + 3
xn+1 =
3x2n − 1

with x0 = 0, it produces

x1 = −3.00000
x2 = −1.96154
x3 = −1.14718
x4 = −0.00658
x5 = −3.00039
x6 = −1.96182
x7 = −1.14743

The sequence will not converge. But if the algorithm starts with x0 = 2, then it produces

x1 = 1.7272727
x2 = 1.6736912
x3 = 1.6717026
x4 = 1.6716999
x5 = 1.6716999

The sequence converges to the root 1.6716999 correct to seven decimal places.
This examples illustrates again that the starting point x0 must be close enough to the zero of
f.

D.W-Ddumba numerical analysis i Page 37 of 133

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.11 Pitfall - Choice of x0 and f 0 (x) = 0

Use the Newton Method to find a non-zero solution of

x = 2 sin x (3.6)

Let f (x) = x − 2 sin x. Then f 0 (x) = 1 − 2 cos x, and the Newton-Raphson iteration is

f (xn ) xn − 2 sin xn 2 (sin xn − xn cos xn )

xn+1 = xn − 0
= xn − =
f (xn ) 1 − 2 cos xn 1 − 2 cos xn

Let x0 = 1.1. The next six estimates, to 3 decimal places, are:

x1 = 8.453 x3 = 203.384 x5 = −87.471

x2 = 5.256 x4 = 118.019 x6 = −203.637

Things dont look good, and they get worse. It turns out that x35 < −64000000. We could be
stubborn and soldier on. Miracles happen-but not often. (One happens here, around n = 212.)

The trouble was caused by the choice of x0 , let us consider x0 = 1.5. Here are the next six
estimates, to 19 decimal places - to indicate how fast the convergence is

x1 = 2.0765582006304348291 x4 = 1.8954942764727706570

x2 = 1.9105066156590806258 x5 = 1.8954942670339809987

x3 = 1.8956220029878460925 x6 = 1.8954942670339809471

The next iterate x7 agrees with x6 in the first 19 decimal places, indeed in the first 32, and the
true root is equal to x6 to 32 decimal places.

Remark 3.1.3 Note that choosing

π
x0 = ≈ 1.0472
3
leads to immediate disaster, since then the denominator f 0 (x) = 1 − 2 cos x0 = 0 and therefore
x1 does not exist.

Comment 3.1.1 The remedy to this is to rewrite the Equation (3.6) into other forms say
2 1 sin x 1
− = 0 or =
x sin x x 2
works nicely - to have f 0 (x) harder to be too small or almost zero for ∀ x ∈ (0, π) - as we shall
see in the section of successive substitution, Section 3.5

Page 38 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.12 Consider the problem of finding the positive number x with cos x = x3 We
can rephrase that as finding the zero of

f (x) = cos x − x3 ⇒ f 0 (x) = − sin x − 3x2

Since cos x ≤ 1 for all x and x3 > 1 for x > 1, we know that our zero lies in [0, 1]. We try a
starting value of x0 = 0.5.

f (xn ) cos xn − x3n xn sin xn + cos xn + 2x3n

xn+1 = xn − = x n − =
f 0 (xn ) − sin xn − 3x2n sin xn + 3x2n
x0 sin x0 + cos x0 + 2x30 (0.5) sin(0.5) + cos(0.5) + 2(0.5)3
x1 = = = 1.1121416
sin x0 + 3x20 sin(0.5) + 3(0.5)2
x1 sin x1 + cos x1 + 2x31
x2 =
sin x1 + 3x21
(1.1121416) sin(1.1121416) + cos(1.112141637097) + 2(1.1121416)3
= = 0.9096727
sin(1.1121416) + 3(1.1121416)2
x2 sin x2 + cos x2 + 2x32
x3 =
sin x2 + 3x22
(0.9096727) sin(0.9096727) + cos(0.9096727) + 2(0.9096727)3
= = 0.8672638
sin(0.9096727) + 3(0.9096727)2
x3 sin x3 + cos x3 + 2x33
x4 =
sin x3 + 3x23
(0.8672638) sin(0.8672638) + cos(0.8672638) + 2(0.8672638)3
= = 0.8654771
sin(0.8672638) + 3(0.8672638)2
x4 sin x4 + cos x4 + 2x34
x5 =
sin x4 + 3x24
(0.8654771) sin(0.8654771) + cos(0.8654771) + 2(0.8654771)3
= = 0.8654740
sin(0.8654771) + 3(0.8654771)2
x5 sin x5 + cos x5 + 2x35
x6 =
sin x5 + 3x25
(0.8654740) sin(0.8654740) + cos(0.8654740) + 2(0.8654740)3
= = 0.8654740
sin(0.8654740) + 3(0.8654740)2

D.W-Ddumba numerical analysis i Page 39 of 133

 3.1. NEWTON RAPHSON’S METHOD

Example 3.1.13 Use Newtons Method to estimate the point of intersection of

2
y = e−x and y = x.
2
−x2 xn − e−xn
f (x) = x − e , ⇒ xn+1 = xn −
1 + 2xn e−x2n
f (xn ) f (xn )
n xn f (xn ) f 0 (xn ) f 0 (xn )
xn − f 0 (xn )

1 0.500000 -0.27880 1.77880 -0.15673 0.65673

2 0.65673 0.00706 1.85331 0.00381 0.65292
3 0.65292 0.00000 1.85261 0.00000 0.65292
4 0.65292

The solution exists and it is 0.65292

Exercise 3.1 Let f (x) = x2 − a. Show that the Newton Method leads to the recurrence
 
1 a
xn+1 = xn +
2 xn

Exercise 3.2 Newtons equation y 3 − 2y − 5 = 0 has a root near y = 2. Starting with y0 = 2,

compute y1 , y2 , and y3 , the next three Newton-Raphson estimates for the root.

Exercise 3.3 Find the root of the equation

x2 + x − 1 = 0

in the interval [0, 1], giving your answer correct to 4 decimal places.

Exercise 3.4 Show that the cubic equation

2x3 + 3x2 − 3x − 5 = 0

has a real root in the interval [1, 2]. Approximate this root correct to five decimal places using
Newton Raphson’s method.

Exercise 3.5 Use Newton’s method to approximate the root of the equation

g(x) = x3 − 2 sin x

on [0.5, 2].

Page 40 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

3.1.2 General Newton’s algorithm for extracting roots of positive

numbers
Suppose that our interest is to find the rth root of a real positive number A. If x is the value
of this root, then x is related to A by the equation,
1
A r = x or xr = A or xr − A = 0

Let f (x) = xr − A, then x is the root of the nonlinear equation,

f (x) = xr − A = 0, ⇒ f (xn ) = xrn − A, ⇒ f 0 (xn ) = rxr−1

n

Applying the Newton Raphson’s method, we have;

f (xn )
xn+1 = xn −
f 0 (xn )
(xr − A)
xn+1 = xn − n r−1
rxn
rxn − xrn + A
r
=
rxr−1
n
(r − 1)xrn + A
=
rxr−1
n

 
1 A
xn+1 = (r − 1)xn + r−1 (3.7)
r xn
Equation (3.7) is a general formula from which we can obtain quadratically convergent iterative
processes for finding approximations to arbitrary roots of numbers.

Note 3.1.1 The root, or answer interested in is the x. And also what is f (x) is the function
f (x) = 0

Example 3.1.14 When r = 2, the square root of a number A we have

√
x= A ⇒ x2 = A

Thus for r = 2 in the general formula, we get,

 
1 A
xn+1 = xn +
2 xn

Which is the Newton’s square root algorithm for extracting roots of positive numbers.

Example 3.1.15 Use Newton’s square root algorithm to find the square root of 5 correct to
six decimal places.
1
5 2 = x ⇒ x2 − 5 = 0 ⇒ f (x) = x2 − 5 = 0
Or substituting in the general formula r = 2 and A = 5 we get,
 
1 A
xn+1 = xn +
2 xn

D.W-Ddumba numerical analysis i Page 41 of 133

 3.1. NEWTON RAPHSON’S METHOD

Starting with x0 = 2.0 we get,

   
1 5 1 5
x1 = x0 + = 2.000000 + = 2.250000
2 x0 2 2.000000
   
1 5 1 5
x2 = x1 + = 2.250000 + = 2.236111
2 x1 2 2.250000
   
1 5 1 5
x3 = x2 + = 2.236111 + = 2.236068
2 x2 2 2.236111
   
1 5 1 5
x4 = x3 + = 2.236068 + = 2.236068
2 x3 2 2.236068

The value of x2 = 2.236111 is correct only to one decimal place since it agrees with the previous
iterate x1 = 2.250000 only in one decimal place.
However, x3 = 2.236068 is correct to three decimal places since it is in agreement with the
previous iterate x2 in exactly three places of decimal.
But x4 = 2.236068 is exactly the same as x3 . In fact they are exactly the same up to nine
decimal places.
This means that x4 = 2.236068 is the value of the root correct to nine decimal places. Thus
x4 , must also be correct up to six decimal places.
Hence, the value of the root that you state as being correct to six decimal places or nine decimal
places is x4 = 2.236068. Compare with the value obtained from calculator.

Page 42 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

3.1.3 Using Newton’s general formula for roots in finding reciprocals

of numbers
The reciprocal of a number A

A−1 = x ⇒ x−1 = A
1
⇒ =A
x
⇒ f (x) = Ax − 1 = 0 wrong choice, f 0 (x) = c
1 1
⇒ f (x) = − A = 0 ⇒ f 0 (x) = − 2
x x
1
f (xn ) −A
xn+1 = xn − 0 ⇒ xn+1 = xn − xn 1
f (xn ) − x2
n

⇒ xn+1 = xn (2 − Axn )

Alternatively, if we have that r = −1 then x−1 = A (positive number), this means x = 1

A
(reciprocal of A) with r = −1 in the general formula in equation (3.7) then we get,

xn+1 = xn (2 − Axn )

This formula is quadratically convergent and can suitably be applied to calculate the reciprocal
of numbers.

Example 3.1.16 Use Newton’s reciprocal algorithm to find the reciprocal of 3.

xn+1 = xn (2 − Axn ) ⇒ xn+1 = xn (2 − 3xn )

Let x0 = 0.5.

x1 = x0 [2 − 3(x0 )] = 0.50000[2 − 3(0.50000)] = 0.25000

x2 = x1 [2 − 3(x1 )] = 0.25000[2 − 3(0.25000)] = 0.31250
x3 = x2 [2 − 3(x2 )] = 0.31250[2 − 3(0.31250)] = 0.33203
x4 = x3 [2 − 3(x3 )] = 0.33203[2 − 3(0.33203)] = 0.33333
x5 = x4 [2 − 3(x4 )] = 0.33333[2 − 3(0.33333)] = 0.33333
1
Thus x5 = 0.33333 is the value of the reciprocal of 3 i.e 3
correct to five decimal places.

Exercise 3.6 Use the Newton Raphson method to estimate the reciprocal of 18 correct to 8
decimal places.

D.W-Ddumba numerical analysis i Page 43 of 133

 3.1. NEWTON RAPHSON’S METHOD

3.1.4 Advantages of Newton Raphson’s method

• A faster method for converging on a single root of a function is the Newton- Raphson
method.

• Perhaps it is the most widely used method of all locating formulas.

• Because the convergence rate is high, no worry for initial guess, interval size, and number
of decimal places required.

• Requires only one guess

3.1.5 Some limitations of the Newton Raphson’s method

Good though it is, the method has some limitations.

? The Newton Raphson’s method may also fail if f (x) has a point of inflection in the neigh-
borhood of the root.

? Newtons method is an extremely powerful technique, but it has a major weakness: the
need to know the value of the derivative of f at each approximation.

? Frequently, f 0 (x) is far more difficult and needs more arithmetic operations to calculate
than f (x).

? If in the immediate neighborhood of a root of f (x), f 0 (x) vanishes or is very small, the
Newton Raphson’s method will not converge. The reason for this failure is; since f 0 (x) is
very small, the quantity, − ff0(x n)
(xn )
becomes very large. In case it is zero - the denominator,
then function is not defined. The consequence is that we are thrown away from the root
we are approximating.

Exercise 3.7 Use Newton’s square root algorithm to find the square root of 2 correct to 6
decimal places.

Exercise 3.8 Use cube root Newton’s algorithm to find the cube root of 7 correct to four
decimal places.

Exercise 3.9 Use Newton’s reciprocal algorithm to find

(a) the reciprocal of the square root of 2.

(b) The reciprocal of the cube root of 4.

Exercise 3.10 State the advantages and disadvantages of the Newton’s method for nonlinear
equations.

Page 44 of 133 D.W-Ddumba numerical analysis i

 3.1. NEWTON RAPHSON’S METHOD

Exercise 3.11

(a) Define the Newton-Raphson method formula for finding the root of a non-linear equation
f (x) = 0

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

(b) The convergence of the Newton-Raphson method technique highly depends on the initial
guess. Discuss.
Yes, when the initial guess is in the interval given, the iterations converge faster than
otherwise.

(c) Use Newton Raphson method to estimate one of the solutions of x2 − 4 = 0 using x0 = 6
to 2 decimal places.

xn xn+1
x0 = 6 x1 = 3.33
x1 = 3.33 x2 = 2.27
x2 = 2.27 x3 = 2.01
x3 = 2.01 x4 = 2.00
x4 = 2.00 x5 = 2.00

(d) Newton’s Raphson’s method is one of the popular schemes for solving a non-linear equation
f (x) = 0. Prove that the Newton Raphson’s method for finding the square root of a positive
number A is given by,  
1 A
xn+1 = xn +
2 xn
√
Use the scheme above to approximate the square root of 5( 5) to three decimal places
with x0 = 2.

xn xn+1
x0 = 2 x1 = 2.25
x1 = 2.25 x2 = 2.236
x2 = 2.236 x3 = 2.2.236
x3 = 2.236 x4 = 2.236

Exercise 3.12 With any simple example, write short notes on the ”NRM” for solving non-
linear equations.

D.W-Ddumba numerical analysis i Page 45 of 133

 3.2. BISECTION METHOD (INTERVAL HALVING)

3.2 Bisection Method (Interval Halving)

3.2.1 Background
The bisection method is one of the bracketing methods for finding roots of equations.

3.2.2 Implementation
Given a function f(x) and an interval which might contain a root, perform a predetermined
number of iterations using the bisection method.

3.2.3 Limitations
Investigate the result of applying the bisection method over an interval where there is a discon-
tinuity. Apply the bisection method for a function using an interval where there are distinct
roots. Apply the bisection method over a ”large” interval.

3.2.4 Explanation on the bisection method

The bisection method takes a similar geometrical approach with the Regula falsi algorithm.
You need two initial guesses x0 and x1 to the root x? of the nonlinear equation f (x) = 0, such
that
f (x0 )f (x1 ) < 0
The next approximation is obtained by getting the arithmetic mean of the previous two. That
is
xn + xn−1
xn+1 = (3.8)
2
However, the pair xn , xn+1 to be used to get xr+2 must satisfy the condition

f (xn )f (xn−1 ) < 0 (3.9)

Masenge (1989) called this method a trivial simplification of the regula falsi.
In mathematics, the bisection method is a root-finding algorithm which works by repeatedly
dividing an interval in half and then selecting the subinterval in which the root exists.

Page 46 of 133 D.W-Ddumba numerical analysis i

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.1 Use the Bisection method to estimate the root of x2 = 3

The non linear function x2 − 3 = 0 ⇒ f (x) = x2 − 3
Let

x0 = 1 ⇒ f0 = −2.000 < 0
x1 = 2 ⇒ f1 = 1.000 > 0

For the Bisection scheme,

xn + xn−1
xn+1 =
2
x1 + x0 1+2
x2 = = = 1.500 ⇒ f2 = −0.75 < 0
2 2
x2 + x1 1.500 + 2.000
x3 = = = 1.750 ⇒ f3 = 0.062 > 0
2 2
x3 + x2 1.750 + 1.500
x4 = = = 1.625 ⇒ f4 = −0.359 < 0
2 2
x4 + x3 1.625 + 1.750
x5 = = = 1.688 ⇒ f5 = −0.151 < 0
2 2
x5 + x3 1.688 + 1.750
x6 = = = 1.719 ⇒ f6 = −0.045 < 0
2 2
x6 + x3 1.719 + 1.750
x7 = = = 1.735 ⇒ f7 = 0.010 > 0
2 2
x7 + x6 1.735 + 1.719
x8 = = = 1.727 ⇒ f8 = −0.017 < 0
2 2
x8 + x7 1.727 + 1.735
x9 = = = 1.731 ⇒ f9 = −0.004 < 0
2 2
x9 + x7 1.731 + 1.735
x10 = = = 1.733 ⇒ f10 = 0.003 > 0
2 2
x10 + x9 1.733 + 1.731
x11 = = = 1.732 ⇒ f11 = 0.000
2 2

Since f11 = 0.000, then x11 = 1.732 it is the root or zero of the function x2 − 3 = 0, that is the
approximated square root of 3 correct to three decimal places.

D.W-Ddumba numerical analysis i Page 47 of 133

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.2 Compute the numerical root of

3x + sin x − ex = 0

using the Bisection method on the interval [0, 0.5] correct to 4 decimal places.
Let

x0 = 0 ⇒ f0 = −1.0000 < 0
x1 = 0.5 ⇒ f1 = 0.3307 > 0

For the Bisection scheme,

xn + xn−1
xn+1 =
2
x1 + x0 0 + 0.5
x2 = = = 0.2500 ⇒ f2 = −0.2866 < 0
2 2
x2 + x1 0.2500 + 0.5000
x3 = = = 0.3750 ⇒ f3 = 0.0363 > 0
2 2
x3 + x2 0.3750 + 0.2500
x4 = = = 0.3125 ⇒ f4 = −0.1219 < 0
2 2
x4 + x3 0.3125 + 0.3750
x5 = = = 0.3438 ⇒ f5 = −0.0418 < 0
2 2
x5 + x3 0.3438 + 0.3750
x6 = = = 0.3594 ⇒ f6 = −0.0026 < 0
2 2
x6 + x3 0.3594 + 0.3750
x7 = = = 0.3672 ⇒ f7 = 0.0169 > 0
2 2
x7 + x6 0.3672 + 0.3594
x8 = = = 0.3633 ⇒ f8 = 0.0072 > 0
2 2
x8 + x6 0.3633 + 0.3594
x9 = = = 0.3614 ⇒ f9 = 0.0024 > 0
2 2
x9 + x6 0.3614 + 0.3594
x10 = = = 0.3604 ⇒ f10 = −0.0001 < 0
2 2
x10 + x9 0.3604 + 0.3614
x11 = = = 0.3609 ⇒ f11 = 0.0012 > 0
2 2
x11 + x10 0.3609 + 0.3604
x12 = = = 0.3607 ⇒ f12 = 0.0007 > 0
2 2
x12 + x10 0.3607 + 0.3604
x13 = = = 0.3606 ⇒ f13 = 0.0004 > 0
2 2
x13 + x10 0.3606 + 0.3604
x14 = = = 0.3605 ⇒ f14 = 0.0002 > 0
2 2
x14 + x10 0.3605 + 0.3604
x15 = = = 0.3604 ⇒ f15 = −0.0001 < 1
2 2
x15 + x14 0.3604 + 0.3605
x16 = = = 0.3604
2 2
So one of the roots of
3x + sin x − ex = 0
is approximately
0.3604
The exact solution is in fact 0.36044

Remark 3.2.1 The Bisection method takes long to converge.

Page 48 of 133 D.W-Ddumba numerical analysis i

 3.2. BISECTION METHOD (INTERVAL HALVING)

3.2.5 Advantages of the Bisection method

• The method is simple.

• The method is always convergent.

3.2.6 Disadvantages of the Bisection method

• It requires the values of a = x0 and b = x1 .

• The convergence of interval halving is very slow (slow at converging to the root x? ). For
example, a simple non linear equation

x3 − x − 2

will converge to about 1.521 after 15-fifteen iterations.

• The method fails in case of approximating a double root or a root of even multiplicity.

• A faster method for converging on a single root of a function is the Newton-Raphson

iteration Method.

D.W-Ddumba numerical analysis i Page 49 of 133

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.3 Using the Bisection algorithm, find the root of

cos x − xex = 0

correct to 3 decimal places on the interval 0 ≤ x ≤ 1

Let x0 = 0 ⇒ f0 = 1.000 > 0 and x1 = 1 ⇒ f1 = −2.178 < 0, for the Bisection scheme,
xn + xn−1
xn+1 =
2
x1 + x0 0+1
x2 = = = 0.500 ⇒ f2 = 0.053 > 0
2 2
x2 + x1 0.500 + 1.000
x3 = = = 0.750 ⇒ f3 = −0.856 < 0
2 2
x3 + x2 0.750 + 0.500
x4 = = = 0.625 ⇒ f4 = −0.357 < 0
2 2
x4 + x2 0.625 + 0.500
x5 = = = 0.562 ⇒ f5 = −0.140 < 0
2 2
x5 + x2 0.562 + 0.500
x6 = = = 0.531 ⇒ f6 = −0.041 < 0
2 2
x6 + x2 0.531 + 0.500
x7 = = = 0.516 ⇒ f7 = 0.005 > 0
2 2
x7 + x6 0.516 + 0.531
x8 = = = 0.524 ⇒ f8 = −0.019 < 0
2 2
x8 + x6 0.524 + 0.516
x9 = = = 0.520 ⇒ f9 = −0.007 < 0
2 2
x9 + x6 0.520 + 0.516
x10 = = = 0.518 ⇒ f10 = −0.001 < 0
2 2
x10 + x6 0.518 + 0.516
x11 = = = 0.517 ⇒ f11 = 0.002 > 0
2 2
x11 + x10 0.517 + 0.518
x12 = = = 0.518 ⇒ f12 = −0.001 < 0
2 2
x12 + x11 0.518 + 0.517
x13 = = = 0.518
2 2
The scheme converges at x = 0.518. Analytically, the root is 0.5177

Page 50 of 133 D.W-Ddumba numerical analysis i

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.4 Estimate the zero to the non linear equation

f (x) = x3 − 5x2 − 2x + 10

using numerical Bisection method, if the graphical methods found the real root between x = 1
and x = 3:
Let

x0 = 1 ⇒ f0 = 4.0000 > 0
x1 = 3 ⇒ f1 = −14.0000 < 0

By interval halving,
xn + xn−1
xn+1 =
2
x1 + x0 3+1
x2 = = = 2.0000 ⇒ f2 = −6.0000 < 0
2 2
x2 + x0 2.0000 + 1.0000
x3 = = = 1.5000 ⇒ f3 = −0.8750 < 0
2 2
x3 + x0 1.5000 + 1.0000
x4 = = = 1.2500 ⇒ f4 = 1.6406 > 0
2 2
x4 + x3 1.2500 + 1.5000
x5 = = = 1.3750 ⇒ f5 = 0.3965 > 0
2 2
x5 + x3 1.3750 + 1.5000
x6 = = = 1.4375 ⇒ f6 = −0.2366 < 0
2 2
x6 + x5 1.4375 + 1.3750
x7 = = = 1.4063 ⇒ f7 = 0.0802 > 0
2 2
....
..

It is evident that the functional values f (xi ) are approaching zero as the number of iterations
is increased.
After more
√ six iteration the approximated root of 1.40625 compares favorably with the exact
value of 2

D.W-Ddumba numerical analysis i Page 51 of 133

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.5 The following polynomial has a root within the interval 3.75 ≤ x ≤ 5.00

f (x) = x3 − x2 − 10x − 8 = 0

If a tolerance of 0.01(1%) is required, find this root using bisection method.

The Bisection algorithm is given by
xn + xn−1 x s + xe
xn+1 = alternatively xm =
2 2
where xm is the mid point, xs starting point of the two, and xe end point of the two

Iteration
i xs xm xe f (xs ) f (xm ) f (xe ) f (xs )f (xm ) f (xm )f (xe ) error d

1 3.7500 4.3750 5.0000 -6.8281 12.8496 42.0000 - + –

2 3.7500 4.0625 4.3750 -6.8281 1.9182 12.8496 - + 0.31250

3 3.7500 3.9063 4.0625 -6.8281 -2.7166 1.9182 + - 0.15625

4 3.9063 .9844 4.0625 -2.7166 -0.4661 1.9182 + - 0.07813

5 3.9844 4.0234 4.0625 -0.4661 0.7092 1.9182 - + 0.03906

6 3.9844 4.0039 4.0234 -0.4661 0.1174 0.7092 - + 0.01953

7 3.9844 3.9941 4.0039 -0.4661 -0.1754 0.1174 + - 0.00977

8 3.9941 3.9990 4.0039 -0.1754 -0.0293 0.1174 + - 0.00488

9 3.9990 4.0015 4.0039 -0.0293 0.0440 0.1174 - + 0.00244

10 3.9990 4.0002 4.0015 -0.0293 0.0073 0.0440 - + 0.00122

11 3.9990 3.9996 4.0002 -0.0293 -0.0110 0.0073 + - 0.00061

12 3.9996 3.9999 4.0002 -0.0110 -0.0018 0.0073 + - 0.00031

13 3.9999 4.0001 4.0002 -0.0018 0.0027 0.0073 - + 0.00015

14 3.9999 4.0000 4.0001 -0.0018 0.0005 0.0027 - + 0.00008

15 3.9999 4.0000 4.0000 -0.0018 -0.0007 0.0005 + - 0.00004

Page 52 of 133 D.W-Ddumba numerical analysis i

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.6 Consider finding the root of

f (x) = e−x (3.2 sin x − 0.5 cos x)

on the interval [3, 4], this time with step = 0.001, abs = 0.001.

Iteration
i xs xe f (xs ) f (xe ) xm f (xm ) abs

1 3.0000 4.0000 0.047127 -0.038372 3.5000 -0.019757 0.5

2 3.0000 3.5000 0.047127 -0.019757 3.2500 0.0058479 0.25

3 3.2500 3.5000 0.0058479 -0.019757 3.3750 -0.0086808 0.125

4 3.2500 3.3750 0.0058479 -0.0086808 3.3125 -0.0018773 0.0625

5 3.2500 3.3125 0.0058479 -0.0018773 3.2812 0.0018739 0.0313

6 3.2812 3.3125 0.0018739 -0.0018773 3.2968 -0.000024791 0.0156

7 3.2812 3.2968 0.0018739 -0.000024791 3.2890 0.00091736 0.0078

8 3.2890 3.2968 0.00091736 -0.000024791 3.2929 0.00044352 0.0039

9 3.2929 3.2968 0.00044352 -0.000024791 3.2948 0.00021466 0.002

10 3.2948 3.2968 0.00021466 -0.000024791 3.2958 0.000094077 0.001

11 3.2958 3.2968 0.000094077 -0.000024791 3.2963 0.000034799 0.0005

Remark 3.2.2 Thus, after the 11th iteration, we note that the final interval, [3.2958, 3.2968]
has a width less than 0.001 and |f (3.2968)| = 0.000025 < 0.001 and therefore we chose xe =
3.2968 to be our approximation of the root.
Although it is xm = 3.2963 in the last interval, but it is not the better approximation, since
scheme has not yet converged and yet |f (3.2963)| = 0.000035 a much more deviation from the
exact value 0.00000 compared to f (3.2968)

Remark 3.2.3 The Bisection scheme has not yet converged, but iteration terminated because
of the required error achieved.

D.W-Ddumba numerical analysis i Page 53 of 133

 3.2. BISECTION METHOD (INTERVAL HALVING)

Example 3.2.7 The root of

ex − 2 = 0
is known to exist in [0, 2]. Use 8 iterations to find an approximate value of the root (or find an
approximate value of the root to within a tolerance of )

Iteration
i xs xm xe f (xs ) f (xe ) f (xm )

1 0.0000 1.0000000 2.000000 -1.0000 0.7183 5.3891

2 0.0000 0.5000000 1.000000 -1.0000 -0.3513 0.7183
3 0.5000 0.7500000 1.000000 -0.3513 0.1170 0.7183
4 0.5000 0.6250000 0.750000 -0.3513 -0.1318 0.1170
5 0.6250 0.6875000 0.750000 -0.1318 -0.0113 0.1170
6 0.6875 0.7187500 0.750000 -0.0113 0.0519 0.1170
7 0.6875 0.7031250 0.718750 -0.0113 0.0201 0.0519
8 0.6875 0.6953125 0.703125 -0.0113 0.0043 0.0201

Exercise 3.13 Will the Bisection Method applied to f (x) = tan x and initial interval [a, b] =
[1, 2] converge to a root? Why or why not? To which value, if any, will the Bisection Method
converge?
Exercise 3.14 If g(x) = cos x − 2x, and [a1 , b1 ] = [0, 1], use the Bisection Method to compute
x3 . Show your work.
Exercise 3.15 Consider the equations
(a) x5 + x = 1
(b) sin x = 6x + 5
(c) ln x + x2 = 3
Apply two steps of the Bisection Method to find an approximate root within 1/8 of the true
root.
Exercise 3.16 Approximate to 2 decimal places the roots of the following equations using
the bisection method.
(a) x2 = 3
(b) x3 = 2
(c) x4 = 2

Exercise 3.17 The function h(x) = x sin x occurs in the study of damped forced oscillation.
Find the value of x that lies in the interval [0, 2] where the function takes on the value h(x) = 1.
Use interval bisection.
Exercise 3.18 If a = 0.1 and b = 1.0, how many steps of the bisection method are needed to
determine a root in this interval with an error of at most 21 × 10−8 ?
Exercise 3.19 Consider obtaining the root of:
ex + 1 + sin x
f (x) = .
(x − 2)
Show that f (1.9) < 0, f (2.1) > 0 and use the bisection method to obtain the root.

Page 54 of 133 D.W-Ddumba numerical analysis i

 3.2. BISECTION METHOD (INTERVAL HALVING)

Exercise 3.20 Find the real root of the equation

x3 − x2 − x + 1 = 0

using the bisection algorithm.

Exercise 3.21 The bisection method generates intervals [a0 , b0 ], [a1 , b1 ], and so on, which of
these inequalities are true for the root r that is being calculated?

(a) |r − an | ≤ 2|r − bn |

(b) |r − an | ≤ 2−n−1 (b0 − a0 )

(c) |r − bn | ≤ 2−n−1 (b0 − a0 )

(d) 0 ≤ r − an 2−n (b0 − a0 )

(e) |r − 21 (an + bn )| ≤ 2−n−2 (b0 − a0 )

Example 3.2.8 Find all the real solutions to the cubic equation

x3 + 4x2 − 10 = 0

in the interval [1, 2].

Example 3.2.9 Use Newton’s method to find the roots of the cubic polynomial x3 −3x+2 = 0
in the interval

(a) [0, 2]

(b) [−3, −1]

D.W-Ddumba numerical analysis i Page 55 of 133

 3.3. SECANT METHOD (CHORDS METHOD)

3.3 Secant Method (Chords Method)

The Secant method needs two points near the root before the algorithm can be applied. Thus
it is of the form
xr+1 = f (xr , xr−1 ).

3.3.1 Derivation of the Secant method

We linear approximate the graph of y = f (x) in figure (10.1), by a chord passing through the
points A and B. The equation of this chord is,

The equation of the Chord is,

[f (x0 ) − f (x1 )](x − x1 )
y = f (x1 ) =
x0 − x1
This Chord cuts the x-axis at x2 i.e.
f (x0 ) − f (x1 )
−f (x1 ) = (x2 − x1 )
(x0 − x1 )

f (x1 )(x1 − x0 )
giving, x2 = x1 −
f (x1 ) − f (x0 )
or in general we have that, xn+1 = g(xn , xn−1 )
That is
f (xn ) [xn − xn−1 ]
xn+1 = xn −
f (xn ) − f (xn−1 )
xn [f (xn ) − f (xn−1 )] − f (xn ) [xn − xn−1 ]
xn+1 =
f (xn ) − f (xn−1 )
xn f (xn ) − xn f (xn−1 ) − xn f (xn ) + xn−1 f (xn )
xn+1 =
f (xn ) − f (xn−1 )
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )

xn−1 f (xn ) − xn f (xn−1 )

xn+1 = (3.10)
f (xn ) − f (xn−1 )
The error in the (n + 1)th iterate is related to the error in the nth iterate en by the relation
en+1 ' Aekn where k ' 1.618 . . . and A is a constant. This relation suggests that the method
has order of convergence 1.618.

Page 56 of 133 D.W-Ddumba numerical analysis i

 3.3. SECANT METHOD (CHORDS METHOD)

Example 3.3.1 Use the Secant method to find the root near 2 of the equation
x3 − 2x − 5 = 0
Start the iteration with x0 = 1.9, ⇒ f (x0 ) = −1.941 and x1 = 2.0 ⇒ f (x1 ) = −1.000.
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (1.9)(−1.000) − (2.0)(−1.941)
x2 = = = 2.10627, f2 = 0.13166
f (x1 ) − f (x0 ) −1.000 −− 1.941
x1 f (x2 ) − x2 f (x1 ) (2.0)(0.13166) − (2.10627)(−1.000)
x3 = = = 2.09391, f3 = −0.00716
f (x2 ) − f (x1 ) 0.13166 −− 1.000
x2 f (x3 ) − x3 f (x2 ) (2.10627)(−0.00716) − (2.09391)(0.13166)
x4 = = = 2.09455, f4 = −0.00002
f (x3 ) − f (x2 ) −0.00716 − 0.13166
x3 f (x4 ) − x4 f (x3 ) (2.09391)(−0.00002) − (2.09455)(−0.00716)
x5 = = = 2.09455
f (x4 ) − f (x3 ) −0.00002 −− 0.00716
Thus since x4 and x5 are identical to 5 decimal places, so x5 = 2.09455 is the value of the root
correct to five decimal places.
Example 3.3.2 Estimate the root of
x4 − x − 10 = 0
with the initial guess as 1.0 and 2.0 using the numerical Secant scheme. Take solutions to 5
decimal places.
Let
x0 = 1.0 ⇒ f0 = −10.0
x1 = 2.0 ⇒ f1 = 4.0
For the Secant rule
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (1.0)(4.0) − (2.0)(−10.0)
x2 = = = 1.71429
f (x1 ) − f (x0 ) 4.0 −− 10.0
⇒ f2 = −3.07780
x1 f (x2 ) − x2 f (x1 ) (2.0)(−3.07780) − (1.71429)(4.0)
x3 = = = 1.83853
f (x2 ) − f (x1 ) −3.07780 − 4.0
⇒ f3 = −0.41283
x2 f (x3 ) − x3 f (x2 ) (1.71429)(−0.41283) − (1.83853)(−3.07780)
x4 = = = 1.85778
f (x3 ) − f (x2 ) −0.41283 −− 3.07780
⇒ f4 = 0.05401
x3 f (x4 ) − x4 f (x3 ) (1.83853)(0.05401) − (1.85778)(−0.41283)
x5 = = = 1.85555
f (x4 ) − f (x3 ) 0.05401 −− 0.41283
⇒ f5 = −0.00085
x4 f (x5 ) − x5 f (x4 ) (1.85778)(−0.00085) − (1.85555)(0.05401)
x6 = = = 1.85558
f (x5 ) − f (x4 ) −0.00085 − 0.05401
⇒ f6 = −0.00011
x5 f (x6 ) − x6 f (x5 ) (1.85555)(−0.00011) − (1.85558)(−0.00085)
x7 = = = 1.85558
f (x6 ) − f (x5 ) −0.00011 −− 0.00085
So the iterative process converges at 1.85558

D.W-Ddumba numerical analysis i Page 57 of 133

 3.3. SECANT METHOD (CHORDS METHOD)

Note 3.3.1 If the function f (x) has a trigonometric term, the calculators better be in radians.

Example 3.3.3 Approximate the root of

1
x − sin x − =0
2

Let the initial guess be 1.0 and 2.0 using the Secant algorithm.
Let x0 = 1.0 ⇒ f0 = −0.34147 and x1 = 2.0 ⇒ f1 = 0.59070

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (1.0)(0.59070) − (2.0)(−0.34147)
x2 = = = 1.36632
f (x1 ) − f (x0 ) 0.59070 −− 0.34147
⇒ f2 = −0.11285
x1 f (x2 ) − x2 f (x1 ) (2.0)(−0.11285) − (1.36632)(0.59070)
x3 = = = 1.46796
f (x2 ) − f (x1 ) −0.11285 − 0.59070
⇒ f3 = −0.02676
x2 f (x3 ) − x3 f (x2 ) (1.36632)(−0.02676) − (1.46796)(−0.11285)
x4 = = = 1.49955
f (x3 ) − f (x2 ) −0.02676 −− 0.11285
⇒ f4 = 0.00209
x3 f (x4 ) − x4 f (x3 ) (1.46796)(0.00209) − (1.49955)(−0.02676)
x5 = = = 1.49726
f (x4 ) − f (x3 ) 0.00209 −− 0.02676
⇒ f5 = −0.00004
x4 f (x5 ) − x5 f (x4 ) (1.49955)(−0.00004) − (1.49726)(0.00209)
x6 = = = 1.49730
f (x5 ) − f (x4 ) −0.00004 − 0.00209
⇒ f6 = 0.00000
x5 f (x6 ) − x6 f (x5 ) (1.49726)(0.00000) − (1.49730)(−0.00004)
x7 = = = 1.49730
f (x6 ) − f (x5 ) −0.00000 −− 0.00004

So the iterative process converges at 1.49730, and it is almost exact.

Example 3.3.4 Let the initial guess be −2.0 and −1.0, find the root of

(x2 + 5x + 2)e−x + 1 = 0

using the Secant technique approximated to 5 decimal places.

Page 58 of 133 D.W-Ddumba numerical analysis i

 3.3. SECANT METHOD (CHORDS METHOD)

Let x0 = −2.0 ⇒ f0 = −28.55622 and x1 = −1.0 ⇒ f1 = −4.43656

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (−2.0)(−4.43656) − (−1.0)(−28.55622)
x2 = = = −0.81606
f (x1 ) − f (x0 ) −4.43656 −− 28.55622
⇒ f2 = −2.19865
x1 f (x2 ) − x2 f (x1 ) (−1.0)(−2.19865) − (−0.81606)(−4.43656)
x3 = = = −0.63535
f (x2 ) − f (x1 ) −2.19865 −− 4.43656
⇒ f3 = −0.45933
x2 f (x3 ) − x3 f (x2 ) (−0.81606)(−0.45933) − (−0.63535)(−2.19865)
x4 = = = −0.58763
f (x3 ) − f (x2 ) −0.45933 −− 2.19865
⇒ f4 = −0.06695
x3 f (x4 ) − x4 f (x3 ) (−0.63535)(−0.06695) − (−0.58763)(−0.45933)
x5 = = = −0.57949
f (x4 ) − f (x3 ) −0.06695 −− 0.45933
⇒ f5 = −0.00260
x4 f (x5 ) − x5 f (x4 ) (−0.58763)(−0.00260) − (−0.57949)(−0.06695)
x6 = = = −0.57916
f (x5 ) − f (x4 ) −0.00260 −− 0.06695
⇒ f6 = −0.00001
x5 f (x6 ) − x6 f (x5 ) (−0.57949)(−0.00001) − (−0.57916)(−0.00260)
x7 = = = −0.57916
f (x6 ) − f (x5 ) −0.00001 −− 0.00260

The approximated zero of

(x2 + 5x + 2)e−x + 1 = 0
is
−0.57916
to five decimal places.

D.W-Ddumba numerical analysis i Page 59 of 133

 3.3. SECANT METHOD (CHORDS METHOD)

Example 3.3.5 Apply the Secant method to show that for the non linear equation

cos x − xex = 0

the approximated roots are given by

x0 x1 x2 x3 x4 x5 x6 x7 x8

1.0 2.0 0.83267 0.72878 0.56240 0.52478 0.51801 0.51776 0.51776

So the iterative process converges at 0.51776

Example 3.3.6 With initial conditions as x0 = 1.0 and x1 = 2.0, iterate with Secant method
to show that for

x − e−x = 0

the inexact solutions approximated to 5 decimal places are

x0 x1 x2 x3 x4 x5 x6

1.0 2.0 0.48714 0.58378 0.56739 0.56714 0.56714

So one of the roots of f (x) = x − e−x is approximately 0.56714

Example 3.3.7 Find the zero of the non linear equation

e−x = 3 log10 x

by the Secant method to 5 decimal places with the initial guess as 1.0 and 2.0
x0 x1 x2 x3 x4 x5 x6 x7

1.0 2.0 1.32394 1.22325 1.24759 1.24683 1.24682 1.24682

So one of the roots of e−x = 3 log10 x is approximately 1.24682.

Hint: The function f (x) used is

f (x) = e−x − 3 log10 x = 0

Page 60 of 133 D.W-Ddumba numerical analysis i

 3.3. SECANT METHOD (CHORDS METHOD)

Example 3.3.8 Use the Secant method to find a solution to x = cos x , and compare the
approximations with those given by Newtons method with x0 = π/4.
For the Secant method we need two initial approximations. Suppose we use x0 = 0.5 and
x1 = π/4.

n xn−1 xn xn+1 |xn+1 − xn |

1 0.500000000 0.785398163 0.736384139 0.0490140246

2 0.785398163 0.736384139 0.739058139 0.0026740004

3 0.736384139 0.739085149 0.739085149 0.0000270101

4 0.739058139 0.739085149 0.739085133 0.0000000161

The Newton’s method for x = cos x with x0 = π/4 is given by

n xn f (xn ) f 0 (xn ) xn+1 |xn+1 − xn |

0 0.78539816 -0.078291 -1.707107 0.73953613 0.04586203

1 0.73953613 -0.000755 -1.673945 0.73908518 0.00045096

2 0.73908518 -0.000000 -1.673612 0.73908513 0.00000004

3 0.73908513 -0.000000 -1.673612 0.73908513 0.00000000

Observations:

• For Newton’s method, an excellent approximation is obtained with n = 2

• Because of the agreement of x2 and x3 we could reasonably expect this result to be

accurate to the places listed.

• Comparing results, we see that the Secant Method approximation x4 is accurate to the
tenth decimal place, whereas Newtons method obtained this accuracy by x2 .

• Here, the convergence of the Secant method is much faster than functional iteration but
slightly slower than Newtons method.

D.W-Ddumba numerical analysis i Page 61 of 133

 3.3. SECANT METHOD (CHORDS METHOD)

3.3.2 Advantages and Disadvantages of the Secant method

The method;

(i) can work for double roots.

(ii) has order of convergence of 1.618.

(iii) is not always convergent.

The above are advantages or disadvantages depending on the comparison technique in question.

Remark 3.3.1 The Secant method and Newtons method are often used to refine an answer
obtained by another technique (such as the Bisection Method).

Exercise 3.22 Find the real root of

f (x) = x3 + x2 − 3x − 3 = 0

using Secant technique correct to 2 decimal places.

Exercise 3.23 Show that there is a root of the equation

f (x) = 3x − sin x − ex = 0

in the interval (0, 1). Estimate this root to 2 decimal places using the Secant method.

Exercise 3.24 Find the roots of the following equations, using the methods of Secant.

(a) ex = cos x

(b) x3 − 2x + 1 = 0

(c) sin 2x − ex − 1 = 0

(d) ln(x − 1) = x2

Exercise 3.25 Use the Secant method to find the real root of the equation

x3 + 2x2 − x + 5 = 0.

Exercise 3.26 Consider obtaining the root of;

ex + sin x
f (x) =
(x − 2)

Show that f (1.9) < 0, f (2.1) > 0 and use the Secant method to obtain the root.

Page 62 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

3.4 The Regula-Falsi method (Method of false position)

The regula falsi algorithm uses two points near the root before the algorithm is applied (a
similar geometric approach like the Secant method-the regula falsi method is an associate of
the Secant method) with the exception that,

f (xr )f (xr−1 ) < 0

at each stage of the algorithm. This is termed as root bracketing, like in the Bisection method.
Thus to start the method, you need two points x0 = a0 and x1 = b0 near the root such that

f (a0 )f (b1 ) < 0

That is f (a0 ) and f (b0 ) are of opposite signs, which implies by the intermediate value theorem
that the function f has a root in the interval [a0 , b0 ], assuming continuity of the function f .
The method proceeds by producing a sequence of shrinking intervals [ak , bk ] that all contain a
root of f .

3.4.1 Geometric representation and derivation of the Regula falsi

algorithm

Figure 3.2: The first two iterations of the false position method

From Figure (3.2) above, we note that the produce f (x0 )(f (x1 ) < 0, which is in conformity
with the regula falsi. The equation of the Chord CD is,
[f (x0 ) − f (x1 )](x − x1 )
y = f (x1 ) =
x0 − x1
This Chord cuts the x-axis at x2 i.e.
f (x0 ) − f (x1 )
−f (x1 ) = (x2 − x1 )
(x0 − x1 )
f (x1 )(x1 − x0 )
giving, x2 = x1 −
f (x1 ) − f (x0 )

D.W-Ddumba numerical analysis i Page 63 of 133

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

or in general we have that, xn+1 = g(xn , xn−1 )

That is

f (xn ) [xn − xn−1 ]

xn+1 = xn −
f (xn ) − f (xn−1 )
xn [f (xn ) − f (xn−1 )] − f (xn ) [xn − xn−1 ]
xn+1 =
f (xn ) − f (xn−1 )
xn f (xn ) − xn f (xn−1 ) − xn f (xn ) + xn−1 f (xn )
xn+1 =
f (xn ) − f (xn−1 )
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )

xn−1 f (xn ) − xn f (xn−1 )

xn+1 = (3.11)
f (xn ) − f (xn−1 )
Equation (3.11) is the popular Regula false/falsi position method with condition that at each
stage of the algorithm,

f (xn )f (xn−1 ) < 0 (3.12)

Page 64 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.1 Find the root between (2, 3) of x3 − 2x − 5 = 0, by using regula falsi method.
Approximate values to 3 decimal places.
Let us take x0 = 2 and x1 = 3.

f (x) = x3 − 2x − 5
f (x0 ) = f (2) = 23 − 2(2) − 5 = −1 < 0 (negative)
f (x1 ) = f (3) = 33 − 2(3) − 5 = 16 > 0 (positive)

The first approximation to root is x2 and is given by

xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) 2f (3) − (3)f (2) 2(16) − (3)(−1) 35
x2 = = = = = 2.059
f (x1 ) − f (x0 ) f (3) − f (2) 16 − (−1) 17
⇒ f (2.059) = −0.389 < 0

The root lies between x2 = 2.059 and x1 = 3

Taking x2 = 2.059 and x1 = 3. we have the second approximation to the root given by
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )
x1 f (x2 ) − x2 f (x1 ) (3)f (2.059) − 2.059f (3) (3)(−0.389) − 2.059(16)
x3 = = = = 2.081
f (x2 ) − f (x1 ) f (2.059) − f (3) (−0.389) − 16
⇒ f (2.081) = −0.150 < 0

The root lies between x3 = 2.081 and x1 = 3

Taking x3 = 2.081 and x1 = 3. we have the third approximation to the root given by
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )
x1 f (x3 ) − x3 f (x1 ) (3)f (2.081) − 2.081f (3) (3)(−0.150) − 2.081(16)
x4 = = = = 2.090
f (x3 ) − f (x1 ) f (2.081) − f (3) (−0.15) − 16
⇒ f (2.090) = −0.051 < 0

The root lies between x4 = 2.090 and x1 = 3

Taking x4 = 2.090 and x1 = 3. we have the fourth approximation to the root given by
xn−1 f (xn ) − xn f (xn−1 )
xn+1 =
f (xn ) − f (xn−1 )
x1 f (x4 ) − x4 f (x1 ) (3)f (2.090) − 2.090f (3) (3)(−0.051) − (2.090)(16)
x5 = = = = 2.091
f (x4 ) − f (x1 ) f (2.090) − f (3) (−0.051) − 16
If we need a root to 3 decimal places, we could still continue with the iterations (the root lies
between x5 = 2.091 and x1 = 3), but if to 2 decimal places, the required root is 2.09
Example 3.4.2 Approximate a root of

3x + sin x − ex = 0

using the regula falsi scheme of solving non linear equations f (x) = 0 to 3 decimal places on
the interval [0, 0.5].

D.W-Ddumba numerical analysis i Page 65 of 133

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

If we sketch the function f (x) it’s clear that there is a root between 0 and 0.5 and also another
root between 1.5 and 2.0. Now let us consider the function f (x) in the interval [0, 0.5] where
f (0) × f (0.5) is less than zero and use the regula-falsi scheme to obtain the zero of f (x) = 0.
Let x0 = 0.0 ⇒ f0 = −1.000 < 0 and x1 = −0.5 ⇒ f1 = 0.331 > 0

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (0.0)(0.331) − (0.5)(−1.000)
x2 = = = 0.376, ⇒ f2 = 0.039 > 0
f (x1 ) − f (x0 ) 0.331 −− 1.000
x0 f (x2 ) − x2 f (x0 ) (0.0)(0.039) − (0.376)(−1.000)
x3 = = = 0.362, ⇒ f3 = 0.004 > 0
f (x2 ) − f (x0 ) 0.039 −− 1.000
x0 f (x3 ) − x3 f (x0 ) (0.0)(0.004) − (0.362)(−1.000)
x4 = = = 0.361, ⇒ f4 = 0.001 > 0
f (x3 ) − f (x0 ) 0.004 −− 1.000
x0 f (x4 ) − x4 f (x0 ) (0.0)(0.001) − (0.361)(−1.000)
x5 = = = 0.360, ⇒ f5 = −0.001 < 0
f (x4 ) − f (x0 ) 0.001 −− 1.000
x4 f (x5 ) − x5 f (x4 ) (0.361)(−0.001) − (0.360)(0.001)
x5 = = = 0.360
f (x5 ) − f (x4 ) −0.001 − 0.001

So one of the roots of

3x + sin x − ex = 0
is approximately 0.360.

Page 66 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.3 Estimate the zero of the equation

 
x
x cos = 0
x−2

with the initial guess of x0 = 1 and x1 = 1.5 to 3 decimal places using the popular regula falsi
technique.
Let

x0 = 1.0 ⇒ f0 = 0.540 > 0

x1 = 1.5 ⇒ f1 = −1.485 < 0

For the regula falsi

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (1.0)(−1.485) − (1.5)(0.540)
x2 = = = 1.133, ⇒ f2 = 0.296 > 0
f (x1 ) − f (x0 ) −1.485 − 0.540
x1 f (x2 ) − x2 f (x1 ) (1.5)(0.296) − (1.133)(−1.485)
x3 = = = 1.194, ⇒ f3 = 0.107 > 0
f (x2 ) − f (x1 ) 0.296 −− 1.485
x1 f (x3 ) − x3 f (x1 ) (1.5)(0.107) − (1.194)(−1.485)
x4 = = = 1.214, ⇒ f4 = 0.032 > 0
f (x3 ) − f (x1 ) 0.107 −− 1.485
x1 f (x4 ) − x4 f (x1 ) (1.5)(0.032) − (1.214)(−1.485)
x5 = = = 1.220, ⇒ f5 = 0.008 > 0
f (x4 ) − f (x1 ) 0.032 −− 1.485
x1 f (x5 ) − x5 f (x1 ) (1.5)(0.008) − (1.220)(−1.485)
x6 = = = 1.222, ⇒ f6 = 0.000 > 0
f (x5 ) − f (x1 ) 0.008 −− 1.485
x1 f (x6 ) − x6 f (x1 ) (1.5)(0.000) − (1.222)(−1.485)
x7 = = = 1.222
f (x6 ) − f (x1 ) 0.000 −− 1.485

So one of the roots of  

x
x cos =0
x−2
is approximately
1.222
to 3 decimal places.

D.W-Ddumba numerical analysis i Page 67 of 133

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.4 Use the Regula-Falsi method to compute a real root of the equation

x3 − 9x + 1 = 0,

(a) if the root lies between 2 and 4

(b) if the root lies between 2 and 3.

Comment on the results.

Let x0 = 2.0 ⇒ f0 = −9 < 0 and x1 = 4.0 ⇒ f1 = 29 > 0, For the regula falsi

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 )
x2 = = 2.4736, ⇒ f2 = −6.12644 < 0
f (x1 ) − f (x0 )
x1 f (x2 ) − x2 f (x1 )
x3 = = 2.73989, ⇒ f3 = −3.090707 < 0
f (x2 ) − f (x1 )
x1 f (x3 ) − x3 f (x1 )
x4 = = 2.86125, ⇒ f4 = −1.326868 < 0
f (x3 ) − f (x1 )
..
.

Let x0 = 2.0 ⇒ f0 = −9 < 0 and x1 = 3.0 ⇒ f1 = 1 > 0, For the regula falsi

x0 f (x1 ) − x1 f (x0 )
x2 = = 2.90000, ⇒ f2 = −0.7110 < 0
f (x1 ) − f (x0 )
x1 f (x2 ) − x2 f (x1 )
x3 = = 2.94156, ⇒ f3 = −0.0207 < 0
f (x2 ) − f (x1 )
x1 f (x3 ) − x3 f (x1 )
x4 = = 2.94275, ⇒ f4 = −0.0011896 < 0
f (x3 ) − f (x1 )
..
.

We observe that the value of the root as a third approximation is evidently different in both
the cases, while the value of x4 , when the interval considered is (2, 3), is closer to the root.
Important observation: The initial interval (x0 , x1 ) in which the root of the equation lies should
be sufficiently small.

Page 68 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.5 Use Regula-Falsi method to find a real root of the equation

ln x − cos x = 0

accurate to four decimal places after three successive approximations.

Let

x0 = 1.0 ⇒ f0 = −0.540302 < 0

x1 = 2.0 ⇒ f1 = 1.109 > 0

For the regula falsi

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (1.0)(1.109) − (2.0)(−0.540302)
x2 = = = 1.3275
f (x1 ) − f (x0 ) 1.109 −− 0.540302
⇒ f2 = 0.0424 > 0
Next iteration to be for x0 and x2 since f0 f2 < 0, change signs
x0 f (x2 ) − x2 f (x0 ) (1.0)(0.0424) − (1.3275)(−0.540302)
x3 = = = 1.3037
f (x2 ) − f (x0 ) 0.0424 −− 0.540302
⇒ f3 = 0.001248 > 0
Next iteration to be for x0 and x3 since f0 f3 < 0, change signs
x0 f (x3 ) − x3 f (x0 ) (1.0)(0.001248) − (1.3037)(−0.540302)
x4 = = = 1.3030
f (x3 ) − f (x0 ) 0.001248 −− 0.540302

The required real root is 1.3030

Note 3.4.1 The iterations have not converged, but we were interested in up to the third
iteration.

Note 3.4.2 For Matlab program for the non linear equation, we write ln x as log x, where as
log x as log10 x

D.W-Ddumba numerical analysis i Page 69 of 133

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.6 Use the method of false position to solve

ex + 2−x + 2 cos x − 6 = 0 1≤x≤2

correct your approximations to 4 decimal places.

Let

x0 = 1.0 ⇒ f0 = −1.7011 < 0

x1 = 2.0 ⇒ f1 = 0.8068 > 0

For the regula falsi

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (1.0)(0.8068) − (2.0)(−1.7011)
x2 = = = 1.6783, ⇒ f2 = −0.5457 < 0
f (x1 ) − f (x0 ) 0.8068 −− 1.7011
x1 f (x2 ) − x2 f (x1 ) (2.0)(−0.5457) − (1.6783)(0.8068)
x3 = = = 1.8081, ⇒ f3 = −0.0858 < 0
f (x2 ) − f (x1 ) −0.5457 − 0.8068
x1 f (x3 ) − x3 f (x1 ) (2.0)(−0.0858) − (1.8081)(0.8068)
x4 = = = 1.8265, ⇒ f4 = −0.0118 < 0
f (x3 ) − f (x1 ) −0.0858 − 0.8068
x1 f (x4 ) − x4 f (x1 ) (2.0)(−0.0118) − (1.8265)(0.8068)
x5 = = = 1.8290, ⇒ f5 = −0.0016 < 0
f (x4 ) − f (x1 ) −0.0118 − 0.8068
x1 f (x5 ) − x5 f (x1 ) (2.0)(−0.0016) − (1.8290)(0.8068)
x6 = = = 1.8293, ⇒ f6 = −0.0003 < 0
f (x5 ) − f (x1 ) −0.0016 − 0.8068
x1 f (x6 ) − x6 f (x1 ) (2.0)(−0.0003) − (1.8294)(0.8068)
x7 = = = 1.8294, ⇒ f7 = 0.0001 > 0
f (x6 ) − f (x1 ) −0.0003 − 0.8068
x6 f (x7 ) − x7 f (x6 ) (1.8293)(0.0001) − (1.8294)(−0.0003)
x8 = = = 1.8294
f (x7 ) − f (x6 ) 0.0001 −− 0.0003

So the root of
ex + 2−x + 2 cos x − 6 = 0
is approximately 1.8294 to 4 decimal places. Analytically, the exact value would be 1.82938.

Remark 3.4.1 Look back at how x8 was computed, how and why it was x6 and x7 paired
together, and not x7 with x1 .

Page 70 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.7 Solve the equation

2x cos 2x − (x − 2)2 = 0 2≤x≤3

correct your approximations to 4 decimal places. Perform only four iterations.

Let

x0 = 2.0 ⇒ f0 = −2.6146 < 0

x1 = 3.0 ⇒ f1 = 4.7610 > 0

For the regula falsi

xn−1 f (xn ) − xn f (xn−1 )

xn+1 =
f (xn ) − f (xn−1 )
x0 f (x1 ) − x1 f (x0 ) (2.0)(4.7610) − (3.0)(−2.6146)
x2 = = = 2.3545
f (x1 ) − f (x0 ) 4.7610 −− 2.6146
⇒ f2 = −0.1416 < 0
x1 f (x2 ) − x2 f (x1 ) (3.0)(−0.1416) − (2.3545)(4.7610)
x3 = = = 2.3731
f (x2 ) − f (x1 ) −0.1416 − 4.7610
⇒ f3 = 0.0212 > 0
x2 f (x3 ) − x3 f (x2 ) (2.3545)(0.0212) − (2.3731)(−0.1416)
x4 = = = 2.3707
f (x3 ) − f (x2 ) 0.0212 −− 0.1416
⇒ f4 = 0.0001 > 0
x2 f (x4 ) − x4 f (x2 ) (2.3545)(0.0001) − (2.3707)(−0.1416)
x5 = = = 2.3707
f (x4 ) − f (x2 ) 0.0001 −− 0.1416

Remark 3.4.2 The scheme converges after the fourth iteration and converges to 2.3707.
In fact, one of the classical root of the non linear equation

2x cos 2x − (x − 2)2 = 0

is
2.37069

Remark 3.4.3 In every next iteration, the previous point must be part of it.

D.W-Ddumba numerical analysis i Page 71 of 133

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.8 Using Regula-Falsi algorithm, approximate the real root of

f (x) = x3 − 2x − 2 = 0

Now since, f (1) = −3 < 0 and f (2) = 2 > 0 and f (x) is continuous for all real x, there exists
x? ∈ (1, 2) such that f (x? ) = 0 (the intermediate value theorem).
Let x0 = 1.0 ⇒ f0 = −3.00000 < 0 and x1 = 2.0 ⇒ f1 = 2.00000 > 0,

x0 f (x1 ) − x1 f (x0 ) (1.0)(2.00000) − (2.0)(−3.00000)

x2 = = = 1.60000, ⇒ f2 = −1.10400 < 0
f (x1 ) − f (x0 ) 2.00000 −− 3.00000
x1 f (x2 ) − x2 f (x1 ) (2.0)(−1.10400) − (1.60000)(2.00000)
x3 = = = 1.74227
f (x2 ) − f (x1 ) −1.10400 − 2.00000
⇒ f3 = −0.19587 < 0
x1 f (x3 ) − x3 f (x1 ) (2.0)(−0.19587) − (1.74227)(2.00000)
x4 = = = 1.76526
f (x3 ) − f (x1 ) −0.19587 − 2.00000
⇒ f4 = −0.02972 < 0
x1 f (x4 ) − x4 f (x1 ) (2.0)(−0.02972) − (1.76526)(2.00000)
x5 = = = 1.76870
f (x4 ) − f (x1 ) −0.02972 − 2.00000
⇒ f5 = −0.00438 < 0
x1 f (x5 ) − x5 f (x1 ) (2.0)(−0.00438) − (1.76870)(2.00000)
x6 = = = 1.76921
f (x5 ) − f (x1 ) −0.00438 − 2.00000
⇒ f6 = −0.00061 < 0
x1 f (x6 ) − x6 f (x1 ) (2.0)(−0.00061) − (1.76921)(2.00000)
x7 = = = 1.76928
f (x6 ) − f (x1 ) −0.00061 − 2.00000
⇒ f7 = −0.00009 < 0
x1 f (x7 ) − x7 f (x1 ) (2.0)(−0.00009) − (1.76928)(2.00000)
x8 = = = 1.76929
f (x7 ) − f (x1 ) −0.00009 − 2.00000
⇒ f8 = −0.00002 < 0
x1 f (x8 ) − x8 f (x1 ) (2.0)(−0.00002) − (1.76929)(2.00000)
x9 = = = 1.76929
f (x8 ) − f (x1 ) −0.00002 − 2.00000

The numerical solution is therefore 1.76929 to five decimal points.

Page 72 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

Example 3.4.9 Use the method of False Position to find a solution to x = cos x, and compare
the approximations with those given by the Newtons method and the Secant Method. Let the
initial approximations be x0 = 0.5 and x1 = π4

False Position Secant Newton-Raphson

n xn xn xn
0 0.5 0.5 0.5
1 0.7853981635 0.7853981635 0.7853981635
2 0.7363841388 0.7363841388 0.7395361337
3 0.7390581392 0.7390581392 0.7390851781
4 0.7390848638 0.7390851493 0.7390851332
5 0.7390851305 0.7390851332 0.7390851332
5 0.7390851332 0.7390851332
7 0.7390851332

Note 3.4.3 Note that the False Position and Secant approximations agree through x3 and
that the method of False Position requires an additional iteration to obtain the same accuracy
as the Secant method.
Remark 3.4.4
• The added insurance of the method of False Position commonly requires more calculation
than the Secant method, . . .
• just as the simplification that the Secant method provides over Newtons method usually
comes at the expense of additional iterations.
Example 3.4.10 The function
f (x) = x2 ex − 1
has a root in the interval [0, 1] since f (0)f (1) < 0. The results from the false position and
secant methods, both started with x0 = 0 and x1 = 1, are shown in the table
Iterates False position Secant
x2 0.3679 0.3679
x3 0.5695 0.5695
x4 0.6551 0.7974
x5 0.6868 0.6855
x6 0.6978 0.7012
x7 0.7016 0.7035
It appears from these results that the secant method gives the correct result x = 0.7035 a little
more quickly.

3.4.2 Order of Convergence of the Regula algorithm

The error in the (n + 1)th iterate (denoted en+1 ) is related to the error in the nth iterate en by
the equation,
en+1 ' Aekn
where k = 1. This suggests that although the regula falsi uses the same formula as the
Secant method, the order of convergence of the regula falsi is one compared to ' 1.618 for the
Secant. Thus, the method is slower at converging to the root compared to the Secant method.
However, with the condition f (xr )f (xr−1 ) < 0 at each stage, ensures that the regula falsi is
always convergent which is not the case with Secant method.

D.W-Ddumba numerical analysis i Page 73 of 133

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

3.4.3 Advantages and disadvantages of the regula falsi algorithm

(i) The regula falsi algorithm is always convergent.

(ii) The order of convergence of the method is one.

The two basic points on the advantages and disadvantages. Whether it is an advantage
or a disadvantage it all depends on the comparison in question.
For instance in comparison with the Secant method, it is disadvantageous that the regula
falsi has order of convergence one.
While it is advantageous that it is always convergent.

Exercise 3.27 Find the approximate value of the real root of

x log10 x = 1.2

by regula falsi method

Exercise 3.28 Find the root of the

xex = 3
by regula falsi method and correct to the three decimal places

Exercise 3.29 Find a root which lies between 1 and 2 of

f (x) = x3 + 2x2 + 10x − 20

(Leonardo’s Equation) using the regula falsi method

Exercise 3.30 Use the regula false algorithm to find the root of

f (x) = x2 − 4x + 2 = 0

that lies in the interval (0, 1) and state your answer correct to three decimal places.

Exercise 3.31 Verify that x = 3 is a solution of x = g(x) where

18x
g(x) = .
(x2 + 9)

Use the regula false to approximate this root.

Exercise 3.32 Consider the equation

ex + 1 + sin x
f (x) = =0
(x − 2)

whose root you would want to find. Show that

f (1.9) < 0, f (2.1) > 0

and use the regula false algorithm to compute this root.

Exercise 3.33 Approximate to three decimal places the roots of the following equations using
the regula false algorithm.

Page 74 of 133 D.W-Ddumba numerical analysis i

 3.4. THE REGULA-FALSI METHOD (METHOD OF FALSE POSITION)

(a) x3 = 2 (c) x4 = 2

(b) x2 = 3 (d) x5 = 3

Exercise 3.34

(a) Derive the regula falsi algorithm by clearly giving its geometrical illustration.

(b) What advantages and disadvantages does the Secant method enjoy over other methods so
far considered for solving nonlinear equations.

Exercise 3.35 Use both the Falsi Position and Bisection to solve the same equation

f (d) = 2552 − 30d2 + d3 = 0

and see if there is a difference in the number of steps falsi Position takes to converge versus
Bisection.

Exercise 3.36
Find the solution of x3 + x − 4 = 0 in the interval [1, 4] with accuracy 10−3 . Apply

(a) Newton Raphson method

(b) Bisection algorithm

(c) Regular Falsi scheme

Exercise 3.37 Approximate the solution of x3 − x − 1 = 0 in the interval [1, 2] with accuracy
10−4 with the Bisection method.

Exercise 3.38 Use Newton’s method to solve the equation x3 − x − 1 = 0

Exercise 3.39 The fourth-degree polynomial

f (x) = 230x4 + 18x3 + 9x2 − 221x − 9

has two real zeros, one in [−1, 0] and the other in [0, 1]. Attempt to approximate these zeros
to within 10−6 using the Regular Falsi method.

Summary 3.1 Increased number of decimal points and/or increased size of interval given
both increase the number of iterations as they increase.
However, increased number of decimal points and/or increased size of interval given both im-
prove on the accuracy of the scheme.

D.W-Ddumba numerical analysis i Page 75 of 133

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

3.5 Successive Substitution (Fixed Point Method)

3.5.1 Background knowledge
Successive substitution is one of the iterative techniques for solving nonlinear equations. Itera-
tive techniques start with an initial value/guess x0 to the root α and then using a suitable recur-
rence relation we generate a sequence of approximations {xk }∞ k=o . If the sequence {x0 , x1 , . . .}
converges, then it does so on the required root. Iterative techniques are written in the form

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

if the next iterate xn+1 depends on the previous one xn .

or xn+1 = gr (xn , xn−1 ), n = 1, 2, . . . ,

if the next iterate depends on the previous two i.e. xn and xn−1 .

3.5.2 Successive Substitutions

In the method, we seek the roots of,
f (x) = 0. (3.13)
We try to split f (x) in the form,
f (x) = x − g(x) (3.14)
However, this splitting may not be unique. But not all the different splittings may be useful to
us. We can determine the type of splitting which is useful to a numerical analyst. Now, instead
of solving equation (3.13) we now solve x = g(x). The scheme for solving this problem is given
by the algorithm;
xn+1 = g(xn ), n = 0, 1, 2, . . .
Thus we start with a suitable value x0 , and generate the sequence of approximations

x1 = g(x0 )
x2 = g(x1 )
x3 = g(x2 )
x4 = g(x3 )
..
.
xn+1 = g(xn )
..
.
..
.

That is, the sequence is {x1 , x2 , . . . , xn , . . .}

Page 76 of 133 D.W-Ddumba numerical analysis i

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

Example 3.5.1 Find the real root of the equation

x2 − 2x − 3 = 0

in the interval [2, 4].

Splitting
f (x) = x2 − 2x − 3 = 0
in the form
f (x) = x − g(x) = 0.
(a) Can have the form
3
x = g1 (x) =
x−2
Taking the initial approximation x0 = 4, we get the iterates
3
xn+1 =
xn − 2
3 3
x1 = = = 1.500000
x0 − 2 (4.000000 − 2)
3 3
x2 = = = −6.000000
x1 − 2 (1.500000 − 2)
3 3
x3 = = = −0.375000
x2 − 2 (−6.000000 − 2)
3 3
x4 = = = −1.263158
x3 − 2 (−0.375000 − 2)
3 3
x5 = = = −0.919355
x4 − 2 (−1.263158 − 2)
3 3
x6 = = = −1.027624
x5 − 2 (−0.919355 − 2)
3 3
x7 = = = −0.990876
x6 − 2 (−1.027624 − 2)
..
.

According to the behavior of the iterates, there is no hope for convergence in the interval
[2, 4]. Hence such a rearrangement-choice of g(x)- is not good.
(b) Splitting f (x) = 0 in the form,
x2 − 3
x = g2 (x) =
2
giving the iterative scheme,
x2n − 3
xn+1 =
2
x20 − 3
x1 = = 6.500000
2
x21 − 3
x2 = = 19.625000
2
x2 − 3
x3 = 2 = 191.070312
2
..
.

D.W-Ddumba numerical analysis i Page 77 of 133

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

This shows that the iterates are obviously diverging. Hence such a rearrangement of g2 (x)
is not good too.
(c) Splitting f (x) = 0 in the form,
p
x = g3 (x) = (2x + 3)

Giving the iteration formula, p

xn+1 = (2xn + 3)
Thus with x0 = 4 we get,
p
x1 = (2x0 + 3) = 3.316625
p
x2 = (2x1 + 3) = 3.103748
p
x3 = (2x2 + 3) = 3.034386
p
x4 = (2x3 + 3) = 3.011440
p
x5 = (2x4 + 3) = 3.003811
p
x6 = (2x5 + 3) = 3.001270
p
x7 = (2x6 + 3) = 3.000423
p
x8 = (2x7 + 3) = 3.000141
p
x9 = (2x8 + 3) = 3.000047
p
x10 = (2x9 + 3) = 3.000016
p
x11 = (2x10 + 3) = 3.000005
p
x12 = (2x11 + 3) = 3.000002
p
x13 = (2x12 + 3) = 3.000001
p
x14 = (2x13 + 3) = 3.000000
p
x15 = (2x14 + 3) = 3.000000

In fact this is an arrangement which is giving a sequence of iterates which are converging
to the root α. The sequence is converging to the root x = 3.000000.

Note 3.5.1
One would actually wonder as to whether you have to keep trying the splitting until you get
one which converges to the root. We can test and know the splitting which gives a convergence
sequence of approximations before starting to compute the iterates. This wonderful criterion
is called the convergence criterion for the iterative scheme of the form

xn+1 = g(xn ).

However, before stating the criterion, lets first formerly state what is meant by an iterative
scheme
xn+1 = g(xn )
being convergent.
Definition 3.5.1
An iterative scheme/process. xn+1 = g(xn ) is convergent if,

lim |xn+1 − xn | = 0,
r→∞

otherwise we say that the scheme is divergent.

Page 78 of 133 D.W-Ddumba numerical analysis i

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

3.5.3 Convergence criterion

(i) Let the function g(x) be continuous in a small interval [a, b] containing a simple (single)
root of the function f (x).

(ii) Let also g(x) be differentiable in the open interval (a, b).

If there exists a real number L such that

• 0 ≤ |g 0 (x)| ≤ L ≤ 1 ∀ x ∈ (a, b)

• a ≤ g(x) ≤ b ∀ x ∈ (a, b)

then for an arbitrary starting value x0 taken from (a, b) the iteration formula xn+1 = g(xn ) will
converge.
The rate of convergence of the iteration will depend on the smallness of the constant L relative
to unity.

Example 3.5.2

(i) Given the function f (x) = x3 − sin x = 0, [0, 1], using successive substitution technique,
we can generate
1 sin xn
xn+1 = (sin xn ) 3 & xn+1 =
x2n
Which of the two methods converge? Why?
Since for
1 1 −2
g(x) = (sin x) 3 , g 0 (x) = (sin) 3 cos x ⇒ 0 < g 0 (x) < 1
3
1
which implies, g(x) = (sin x) 3 is the converging formula.

(ii) Use the converging formula in (ii) above, to approximate the root of

f (x) = x3 − sin x = 0

in [0, 1] to 3 decimal places.

[Hint: Let x0 = 1 and use radians.]

1
n (sin x) 3
1 1.0
2 0.944
3 0.932
4 0.929
5 0.929

Note 3.5.2
Sometimes this method is known as the fixed point method, i.e. a fixed point is a number x?
such that x? = g(x? ). Thus a root of f (x) = 0 is a fixed point of the scheme x = g(x).

Example 3.5.3 Given

f (x) = x3 − 7x + 2 = 0
in [0, 1]. Find a sequence that {xn } that converges to the root of f (x) = 0 in [0, 1].

D.W-Ddumba numerical analysis i Page 79 of 133

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

Rewrite f (x) = 0 as x = 71 (x3 + 2). Then

1 3 3
g(x) = (x3 + 2) ⇒ g 0 (x) = x2 < < 1 ∀ x ∈ [0, 1]
7 7 7
Hence by the convergence criterion, the sequence {xn } defined by

1 3
xn+1 = (x + 2)
7 n
converges to a root of x3 − 7x + 2 = 0

Example 3.5.4 Approximate

f (x) = x3 − x − 1 = 0
on (1, 2) to nine decimal places by the successive substitution method.
Note that f (−1) = −1 and f (2) = 5. Therefore by the IVT (Intermediate Value Theorem of
Calculus) a root exists on (0, 2).
Set
1
xn+1 = g(x) = (1 + x) 3
Note that
1 2 1 1
g 0 (x) = (1 + x)− 3 ⇒ 2/3
≤ g 0 (x) ≤
3 3(1 + 2) 3(1 + 1)2/3
1
⇒ 0 ≤ g 0 (x) ≤ < 1 ∀ x ∈ (1, 2)
3(22/3 )
By the convergence criterion, the sequence
1
xn+1 = (1 + xn ) 3

will converge to a fixed point on (1, 2)

Let the initial guess be x0 = 1.3

x1 = (1 + x0 )1/3 = (1 + 1.300000000)1/3 = 1.320006122

x2 = (1 + x1 )1/3 = (1 + 1.320006122)1/3 = 1.323822354
x3 = (1 + x2 )1/3 = (1 + 1.323822354)1/3 = 1.324547818
x4 = (1 + x3 )1/3 = (1 + 1.324547818)1/3 = 1.324685639
x5 = (1 + x4 )1/3 = (1 + 1.324685639)1/3 = 1.324711818
x6 = (1 + x5 )1/3 = (1 + 1.324711818)1/3 = 1.324716791
x7 = (1 + x6 )1/3 = (1 + 1.324716791)1/3 = 1.324717736
x8 = (1 + x7 )1/3 = (1 + 1.324717736)1/3 = 1.324717915
x9 = (1 + x8 )1/3 = (1 + 1.324717915)1/3 = 1.324717949
x10 = (1 + x9 )1/3 = (1 + 1.324717949)1/3 = 1.324717956
x11 = (1 + x10 )1/3 = (1 + 1.324717956)1/3 = 1.324717957
x12 = (1 + x11 )1/3 = (1 + 1.324717957)1/3 = 1.324717957
x13 = (1 + x12 )1/3 = (1 + 1.324717957)1/3 = 1.324717957

It converges to a fixed point 1.324717957

Page 80 of 133 D.W-Ddumba numerical analysis i

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

Example 3.5.5 Find the root of

f (x) = e−x − x = 0

starting with an initial guess of x0 = 0 correct to one decimal places.

We use the substitution
xn+1 = e−xn
to have the solution as
x = 0.5
after six iterations

Remark 3.5.1 Correct to one decimal place is too much an approximation, but for faster
convergence - not to have more than 14 iterations - it is the pay-off between accuracy and
computational difficulties.

D.W-Ddumba numerical analysis i Page 81 of 133

 3.5. SUCCESSIVE SUBSTITUTION (FIXED POINT METHOD)

Find the root of

f (x) = e−x − x = 0
using both the Newton-Raphson and Fixed point iteration
Fixed Point Iteration with
xn+1 = e−xn
n xn
Newton-Raphson 0 0
1 1.000000
n xn 2 0.367879
0 0 3 0.692201
1 0.500000000 4 0.500473
2 0.566311003 5 0.606244
3 0.567143165 6 0.545396
4 0.567143290 7 0.579612
8 0.560115
9 0.571143
10 0.564879

The correct solution is x = 0.56714329. We realise that even at 10th iteration, the fixed point
method has not yet converged.

Exercise 3.40 Verify for the following example that x = 3 is a solution of x = g(x).
18x
(a) g(x) = (x2 +9)

(b) g(x) = x3 − 24
2x3
(c) g(x) = (3x2 −9)

81
(d) g(x) = (x2 +18)

Starting with x0 = 3.1, calculate the first few iteration and justify theoretically the apparent
behavior.
1
Exercise 3.41 Consider the fixed iteration xk = (xk−1 ) for g(x) = 2(x − 1) 2 for x ≥ 1. Show
that only one fixed point exists (at x = 2) and that g 0 (2) = 1. Compute iterations starting
from

(a) x0 = 1.5 and

(b) x0 = 2.5 and show them on a plot of x and g(x).

Exercise 3.42 By splitting f (x) = x3 − x − 1 = 0 in the form f (x) = x − g(x) = 0 for finding
the root α ∈ [1, 2],

(i) get three different splitting and with their corresponding iterative formulae.

(ii) Test using the convergence criterion which of the splitting lead to a convergent sequence.

(iii) For the scheme leading to a convergent sequence, start with a suitable initial approxima-
tion and find the root correct to 3 decimal places.

Page 82 of 133 D.W-Ddumba numerical analysis i

Chapter 4

Interpolation

Linear interpolation is often used to approximate a value of some function f using two known
values of that function at other points.

4.1 Review- Linear interpolation

Using the similarities of triangles,
f (b) − f (a) f (c) − f (a)
=
b−a c−a
(b − a)
f (b) − f (a) = [f (c) − f (a)]
(c − a)

thus
(b − a)
f (b) = f (a) + [f (c) − f (a)] (4.1)
(c − a)
and
 
f (b) − f (a)
b = a+ (c − a) (4.2)
f (c) − f (a)
Example 4.1.1 Given the data below
Time 0 1 2 3 4

Distance 0 6 39 67 100

(i) Find the distance traveled when t = 2.3 hrs.

Then a = 2, c = 3, & b = 2.3
f (a) = 39, f (c) = 67, & f (b)?
(b − a) (2.3 − 2)
f (b) = f (a) + [f (c) − f (a)] = 39 + (67 − 39) = 47.4
(c − a) (3 − 2)

(ii) The time taken when the distance traveled is 80 miles.

Then a = 3, c = 4, & b =?
f (a) = 67, f (c) = 100, & f (b) = 80
 
f (b) − f (a) (80 − 67)
b=a+ (c − a) = 3 + (4 − 1) = 3.39394
f (c) − f (a) (100 − 67)

83
 4.1. REVIEW- LINEAR INTERPOLATION

Example 4.1.2 The bus stages along Kampala-Jinja are 10 km apart. An express bus trav-
eling between the two towns only stops at these stages except in case of an emergency when its
permitted to stop at a point between the two stages.
The fares (fee) between the first, second, third and fourth stages from Jinja are Sh 110, Sh 150, Sh 185
and Sh 200 respectively. On a certain day, a passenger paid to travel from Jinja in the bus
upto the fourth stage, but he fell sick and had to be left on a health center 33 km from Jinja.
(a) Given that he was refunded money for the distance he had not traveled, find the approxi-
mate amount of money he received.

Distance (x) 10 20 30 40

Amount Paid 110 150 185 200

Then a = 30, c = 40, & b = 33

f (a) = 185, f (c) = 200, & f (b)?
(b − a) (200 − 185)
f (b) = f (a) + [f (c) − f (a)] = 185 + (33 − 30) = 189.5
(c − a) (40 − 30)
The journey he had moved costed him Sh 189.5, thus he was refunded
200 − 189.5 = 10.5 shillings

(b) Another person who had only Sh.165 was allowed to board a bus but would be left at a
point worth his money, how far from Jinja would he be left.
Then a = 20, c = 30, & b =?
f (a) = 150, f (c) = 185, & f (b) = 165
 
f (b) − f (a) (165 − 150)
b=a+ (c − a) = 20 + (30 − 20) = 24.286 km
f (c) − f (a) (185 − 150)

Example 4.1.3 The table below shows the values of cos x.

80◦ x0 0 10 20 30 40 50

cos 80◦ x0 0.1736 0.1708 0.1679 0.1650 0.1622 .01593

(i) Find the value of cos 80◦ 350 .

Then a = 300 , c = 400 , & b = 350
f (a) = 0.1650, f (c) = .01622, & f (b)?
(b − a) (−0.0028) 0
f (b) = f (a) + [f (c) − f (a)] = 0.1650 + 5 = .01636
(c − a) 100

(ii) Find the cos−1 0.1655.

Then a = 200 , c = 300 , & b =?
f (a) = 0.1679, f (c) = 0.1650, & f (b) = 0.1655
 
f (b) − f (a) (0.1655 − 0.1679) 0
b = a+ (c − a) = 200 + (30 − 200 ) = 200 + 8.2760 = 28.2760
f (c) − f (a) (0.1650 − 0.1679)

cos−1 0.1655 = 80◦ 28.30

Page 84 of 133 D.W-Ddumba numerical analysis i

 4.2. APPLICATION

Example 4.1.4 Use Linear interpolation to find the root of the equation x3 − x − 1 = 0 which
lies between (1, 2)
x 1 2

f (x) −1 5
Recall that a value x is a root or a solution to an equation f (x) if x satisfy it, i.e if f (x) = 0.
The question is to find the value of x at f (x) = 0.
Then a = 1, c = 2, & b =?
f (a) = −1, f (c) = 5, & f (b) = 0
(0 −− 1)
 
f (b) − f (a)
b=a+ (c − a) = 1 + (2 − 1) = 1.16667
f (c) − f (a) (5 −− 1)

Example 4.1.5 Use linear interpolation to estimate the root of the equation

x3 − 2x − 5 = 0

which lies in the interval (2, 3). 2.0588

Example 4.1.6 Use linear interpolation to estimate the root of the equation

x2 − 2 = 0

which lies in the interval (−2, −1). −1.3333

Example 4.1.7 The following data gives the distance covered by a particle for a certain period
of time.
Time (s) 0 1 2

Distance(m) 0 5 7
Estimate the time taken by a particle to cover a distance of 6 m. 1.0833 s

4.2 Application
Linear interpolation is often used to fill the gaps in a table. Suppose you have a table listing
the population of some country in 1970, 1980, 1990 and 2000, and that you want to estimate
the population in 1994. Linear interpolation gives you an easy way to do this.
The basic operation of linear interpolation between two values is so commonly used in computer
graphics that it is sometimes called a lerp in the jargon of computer graphics. The term can
be used as a verb or noun for the operation. e.g. ”Bresenham’s algorithm lerps incrementally
between the two endpoints of the line.”
Lerp operations are built into the hardware of all modern computer graphics processors. They
are often used as building blocks for more complex operations: for example, a bilinear interpo-
lation can be accomplished in three lerps. Because this operation is cheap, it’s also a good way
to implement accurate lookup tables with quick lookup for smooth functions without having
too many table entries.

D.W-Ddumba numerical analysis i Page 85 of 133

 4.3. HISTORY

4.3 History
Linear interpolation has been used since antiquity for filling the gaps in tables, often with
astronomical data. It is believed that it was used in the Seleucid Empire (last three centuries
BC) and by the Greek astronomer and mathematician Hipparchus (second century BC). A
description of linear interpolation can be found in the Almagest (second century AD) of Ptolemy.

4.4 Extensions
In demanding situations, linear interpolation is often not accurate enough (since not all points
can be approximated to be on a straight line). In that case, it can be replaced by polynomial
interpolation or spline interpolation.
Linear interpolation can also be extended to bilinear interpolation for interpolating functions
of two variables. Bilinear interpolation is often used as a crude anti-aliasing filter. Similarly,
trilinear interpolation is used to interpolate functions of three variables. Other extensions of
linear interpolation can be applied to other kinds of mesh such as triangular and tetrahedral
meshes.

4.5 Lagrange interpolation

We may know the value of a function f at a set of points x1 , x2 , . . . , xN . How do we estimate
the value of the function at any point x∗; how do we compute f (x∗). In general this is done
by developing a smooth curve though the points (x1 , f (x1 )), (x2 , f (x2 )), ...(xN , f (xN )).
Lagrange interpolation is a way to pass a polynomial of degree N − 1 through N points.
Lagrange polynomials are the interpolating polynomials that equal zero in all given points,
save one. Say, given points x1 , x2 , . . . , xN , Lagrange’s polynomial number k is the product
(x − x1 ) (x − x2 ) (x − xk−1 ) (x − xk+1 ) (x − xN )
Pk (x) = ··· ···
(xk − x1 ) (xk − x2 ) (xk − xk−1 ) (xk − xk+1 ) (xk − xN )

such that Pk (xk ) = 1 and Pk (xj ) = 0, for j different from k. In terms of Lagrange’s polynomials
the polynomial interpolation through the points (x1 , y1 ), (x2 , y2 ), . . . , (xN , yN ) could be defined
simply as

P (x) = y1 P1 (x) + y2 P2 (x) + ... + yN PN (x)

Definition 4.5.1 Polynomial interpolation is a method of constructing a function and esti-

mating values at non-tabular points between x0 and xn .

4.5.1 Alternatively
Let f (x) be a continuous function on [a, b], such that

f (xk ) = f (xk ), k = 0, 1, ..., n,

with xk ∈ [a, b]. We call the set of points {xk } tabular points (or interpolating points), while
the set of values {f (xk )} are called the tabular values of f (x). We seek for a polynomial of
degree n, such that
Pn (xk ) = f (xk ), for k = 0 , 1 , 2 , . . . n. (4.3)

Page 86 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

Such a polynomial is called Lagrange’s interpolation polynomial. However, the process of

calculating f (x), x ∈ [a, b], x 6= xk ∀ k from 0 to n, is called interpolation.
We now proceed to derive a formula for Pn (x). However, there is need to introduce some key
definitions.

Definition 4.5.2
Let {Dk : k = 0, 1, ..., n} be any set of numbers. We define the products;
n
Y
Dk = D0 .D1 ...Dn
k=0

n
Y
and Dk = D0 .D1 ...Dj −1 .Dj +1 ...Dn .
k =0 ,k 6=j

This definition is important in the following defining equation for Lagrange’s interpolating
polynomial.

Theorem 4.5.1
The Lagrange’s interpolation polynomial Pn (x) is given by,
n
X
Pn (x) = Lk (x)f (xk ) (4.4)
k=0

n
Y x − xj
where Lk (x ) =
j =0 ,j 6=k
xk − xj

and f (xk ) ≡ f (xk )

Proof
Because Pn (x) is of degree n, we may write
n
Y
Lk (x) = αk (x − xj ),
j=0,j6=k

αk constant.
For Pn (x) to satisfy equation (4.3) we must have,

Lk (xj ) = δkj ( Kroneker delta)

with 
 1, j = k
δkj = (4.5)
0, j 6= k


With the condition (1.3) we have,

1
αk = n
Q
(xk − xj )
j=0,j6=k

This completes the proof.

D.W-Ddumba numerical analysis i Page 87 of 133

 4.5. LAGRANGE INTERPOLATION

4.5.2 Examples of interpolating polynomials

(i) Linear interpolating polynomials
When n = 1, in equation (4.4) we have the polynomial as
1
X
P1 (x) = Lk (x)f (xk )
k=0
= L0 (x)f (x0 ) + L1 (x)f (x1 )
n
Y (x − xj )
But Lk (x) =
j=0 j6=k
(xk − xj )
 
x − x1
therefore L0 (x) =
x − x1
 0 
x − x0
L1 (x) =
x1 − x0

Thus,    
x − x1 x − x0
P1 (x) = f (x0 ) + f (x1 ) (4.6)
x0 − x1 x 1 − x0
Which is known as Lagrange’s interpolating polynomial of degree one, popularly known
as Linear interpolating polynomial.

Example 4.5.1 Construct a linear interpolation polynomial for the data,

x 0 1
f(x) 1.0000 2.7183

Hence interpolate f (0.5)

Solution
Let x0 = 0 and x1 = 1
therefore f (0) = 1.0000 and f (1) = 2.7183
Substituting in equation (4.6) for linear interpolation we have,

P1 (x) = (1 − x)(1.0000) + (x)(2.7183)

= 1.0000 + 1.7183x
and hence, P1 (0.5) = 1.0000 + (1.7183)(0.5)
= 1.0000 + 0.8592 = 1.8592

In fact the data in this particular example describes the graph of f (x) = ex on [0, 1].
The geometrical interpretation of the linear interpolation on [0, 1] .

Page 88 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

(ii) Quadrative interpolating polynomials

When n = 2 in equation (4.4) we get

2
X
P2 (x) = Lk (x)f (xk )
k=0
= L0 (x)f (x0 ) + L1 (x)f (x1 ) + L2 (x)f (x2 )
n
Y (x − xj )
with Lk (x) =
j=0,6=k
(xk − xj )
(x − x1 )(x − x2 )
then L0 (x) =
(x0 − x1 )(x0 − x2 )
(x − x0 )(x − x2 )
L1 (x) =
(x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 )
L2 (x) =
(x2 − x0 )(x2 − x1 )

Thus, we have,
(x − x1 )(x − x2 )
P2 (x) = f (x0 )
(x0 − x1 )(x0 − x2 )

(x − x0 )(x − x2 ) (x − x0 )(x − x1 )
+ f (x1 ) + f (x2 )
(x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Which is called a quadratic interpolating polynomial. Generally they are better interpo-
lating polynomials than the linear ones.

Example 4.5.2 You are given that, f (0) = −2, f (2) = 4, and f (3) = 10. Find a
Lagrange polynomial of degree 2 that fits the data.

Solution
Since x0 = 0, x1 = 2 and x2 = 3 therefore f (x0 ) = −2, f (x1 ) = 4, f (x2 ) = 10.
But

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

(x − x1 ).(x − x2 ) (x − 2)(x − 3)
L0 = =
(x0 − x1 )(x0 − x2 ) (−2)(−3)
1 2 5
= x − x+1
6 6
(x − x0 ).(x − x2 ) (x − 0)(x − 3)
L1 = =
(x1 − x0 )(x1 − x2 ) (2 − 0)(2 − 3)
1 2 3
= − x + x
2 2
(x − x0 ).(x − x1 ) (x − 0)(x − 2)
L2 = =
(x2 − x0 )(x2 − x1 ) (3 − 0)(3 − 2)
1 2 2
= x − x
3 3
D.W-Ddumba numerical analysis i Page 89 of 133
 4.5. LAGRANGE INTERPOLATION

P2 (x) = f (x) = L0 fx0 + L1 fx1 + L2 fx2

     
1 2 5 1 2 3 1 2 2
f (x) = (−2) x − x + 1 + (4) − x + x + (10) x − x
6 6 2 2 3 3
2
= x +x−2

Thus P2 (x) can be used to interpolate f (x) at any of the non tabular points.

Page 90 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

(iii) The Cubic interpolating polynomial

When n = 3 in equation (4.4) we get,

3
X
P3 (x) = Lk (x)f (xk )
k=0
= L0 (x)f (x0 ) + L1 (x)f (x1 ) + L2 (x)f (x2 ) + L3 (x)f (x3 )
n
Y (x − xj )
and with, Lk (x) =
j=0,j6=k
(xk − xj )
(x − x1 )(x − x2 )(x − x3 )
then, L0 (x) =
(x0 − x1 )(x0 − x2 )(x0 − x3 )
(x − x0 )(x − x2 )(x − x3 )
L1 (x) =
(x1 − x0 )(x1 − x2 )(x1 − x3 )
(x − x0 )(x − x1 )(x − x3 )
L2 (x) =
(x2 − x0 )(x2 − x1 )(x2 − x3 )
(x − x0 )(x − x1 )(x − x2 )
L3 (x) =
(x3 − x0 )(x3 − x1 )(x3 − x2 )

Thus,
(x − x1 )(x − x2 )(x − x3 ) (x − x0 )(x − x2 )(x − x3 )
P3 = f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (x1 − x0 )(x1 − x2 )(x1 − x3 )
(x − x0 )(x − x1 )(x − x3 ) (x − x0 )(x − x1 )(x − x2 )
+ f (x2 ) + f (x3 )
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (x3 − x0 )(x3 − x1 )(x3 − x2 )
Construction of the cubic polynomials from available data will be tested in the text
questions at the end of this lecture. Note that the higher degree Lagrange polynomials
can be constructed with ease.

Theorem 4.5.2
Lagrange’s interpolation polynomial Pn (x) is unique.

Note 4.5.1
The uniqueness means that you cannot find any other polynomial of the same degree
which can interpolate the data.
Proof
The proof proceeds by contradiction. Let Pn (x) and Qn (x) be two different polynomials
which interpolate f (x) over the set of points {xk : k = 0, 1, ..., n} which belong to the
interval [a, b], then

Pn (xk ) = Q(xk ) = f (xk ), for k = 0 , 1 , 2 , . . . , n.

Lets define,
r(x) = P (x) − Q(x) ∀ x ∈ [a, b],
then r(x) has at most degree n. Since r(xk ) = 0, for k = 0, 1, 2, ..., n, it has n + 1 distinct
zeros in [a, b]. This contradicts the fundamental law of algebra which states that a non-
zero polynomial of degree n cannot have more than n zeros and so Pn (x) and Qn (x) are
the same polynomials.

D.W-Ddumba numerical analysis i Page 91 of 133

 4.5. LAGRANGE INTERPOLATION

Exercise 4.1
1
1. Construct a linear interpolating polynomial P1 (x) for the function f (x) = x
on the interval
[1, 2]. Use your polynomial to interpolate f (x) at x = 1.2 .

2. Given that x0 = 0, x1 = 21 and x2 = 1 for f (x) = ex . Construct a Lagrange polynomial

that agrees with f(x) at the interpolating points.

3. Find a third degree Lagrange polynomial that goes through the points (0, 0), (1, 1), (8, 2)
and (27, 3).
Use the polynomial to find q for (20, q). Also construct a linear interpolating polynomial
using only (8, 2) and (27, 3), then use the linear polynomial to estimate q. Compare the
1
estimated q s and comment on your results given that the data is of the function y = x 3 .

4. Find the interpolating polynomials going through

(i) (0, 1) and (2, 3)

(ii) (−1, a), (0, b) and (1, c)
(iii) (0, 1), (1, 0), (2, 1) and (3, 0)

5. Given the table below, use Lagrange’s interpolation polynomials of degree one, two and
three to find f (2.5)

x 2.0 2.2 2.4 2.6

f (x) 0.5102 0.5208 0.5104 0.4813

Solution

(i) For n = 1, since need to predict f (2.5) then x0 = 2.4, x2 = 2.6

x − x1 x − x0
P1 (x) = ( )f (x0 ) + ( )f (x1 )
x0 − x1 x1 − x0
x − 2.6 x − 2.4
= ( )0.5104 + ( )0.4813
2.4 − 2.6 2.6 − 2.4
= 0.49585

so for x = 2.5, f (2.5) = 0.49585

(ii) For n = 2, x0 = 2.2, x1 = 2.4, x2 = 2.6

(x − x0 )(x − x2 ) (x − x0 )(x − x1 )
= f (x1 ) + f (x2 )
(x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )
= 0.53814

6. Repeat questions in linear-interpolation (review-section) using the Lagrange linear inter-

polating polynomial. comment on the accuracy of the two techniques.

4.5.3 Error analysis in Lagrange’s interpolation

Truncation errors in Lagrange interpolation
We assume that f (xk ) ≡ f (xk ), k = 0, 1, 2, . . . , n are exact and we consider the truncation
error e(x) = f (x) − Pn (x) for x ∈ [a, b] and Pn (x) the Lagrange polynomial of degree n as
defined in earlier.

Page 92 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

Apart from the fact that e(xk ) = 0, for k = 0, 1, 2, . . . , n with xk ∈ [a, b], we can say nothing
more about e(x) for any x 6= xk .
In addition we have that f (x) has at least n+1 continuous derivatives on [a, b], then it is possible
to express e(x) in terms of f n+1 (x). We now state without proof two necessary Lemmas.

Lemma 4.5.1
Given
n
Y
qn+1 (x) = (x − xk ), x ∈ [a, b]
k=0

of degree (n + 1),
0
qn+1 (x) is of degree n such that
n
X
0
qn+1 (xj ) = (xj − xk )
k=0,k6=j

Lemma 4.5.2
Let a function g(x) be defined on [a, b]. Let {Sk : k = 0, 1, . . . , n} be a set of distinct points
each belonging to [a, b], with S0 < S1 < S2 < . . . < Sn . Suppose that:

(i) g (n) (x)(n ≥ integer) is continuous on [S0 , Sn ].

(ii) g(Sk ) = 0, k = 0, 1, . . . , n.

Then there is at least a number ξ ∈ (S0 , Sn ) such that , g (n) (ξ) = 0

Lemma 4.5.3
The expression in the truncation error e(x) is given by:
n
Y f (n+1) (ξ)
e(x) = (x − xk )
k=0
(n + 1)!

with ξ = ξ(x) and mink (x, xk ) < ξ < maxk (x, xk )

Proof

For points x = xk , k = 0, 1, 2, . . . , n, the theorem is trivially satisfied.

Suppose x 6= xk , we define
Yn
q(t) = (t − xk )
k=0

and
q(t)
g(t) = f (t) − Pn (t) − (f (x) − Pn (x)) for t, x ∈ [a, b].
q(x)
Now g(xk ) = 0, k = 0, 1, 2, . . . n Also g(x) = 0 (x 6= xk ). So for each fixed x 6= xk , g(t) has
n + 2 distinct zeros.
Since g(t) satisfies all the conditions of Lemma 1.2, we deduce that there is a number ξ such
that g (n+1) (ξ) = 0 where
min(x, xk) < ξ < max(x, xk ).
k k

D.W-Ddumba numerical analysis i Page 93 of 133

 4.5. LAGRANGE INTERPOLATION

Now
(f (x) − Pn )q (n+1) (t)
g (n+1) (t) = f (n+1) (t) − Pn(n+1) (t) −
q(x)
(f (x) − Pn (x))
= f (n+1) (t) − (n + 1)!
q(x)
Since
P n+1 (t) = 0
and
qn+1 (t) = (n + 1)!
The result of the theorem follows when t = ξ.
Example 4.5.3 Using two point linear interpolation of ex on [0, 1] therefore
x(x − 1) ξ
ex = e , ξ ∈ (0, 1).
2
But x(x − 1) is maximum or minimum at x = 21 with maximum absolute value equal to 1. So
e(x) has a maximum value equal to 8e = 0.3398 (4dp)

4.5.4 Rounding errors in Lagrange polynomials

Errors known as rounding errors are usually introduced in the functional evaluation of f (xk ).
Suppose rounding errors Ek occur in the data f (xk ) = f (xk ), k = 0, 1, 2, . . . , n respectively.
Let P (x) and P ? (x) denote, respectively, the interpolating polynomials using exact and inexact
data.
Thus,
n
X
Pn (x) = (f (xk ) + Ek )Lk (x)
k=0
n
X
Pn? (x) = f (xk )Lk (x)
k=0
n
X
so |Pn (x ) − Pn? (x )| ≤ |Ek ||Lk (x)|
k=0

Which is the rounding error bound. If data was rounded to m decimal places, then the absolute
maximum error is |Ek | ≤ 21 × 10−m for each k = 0, 1, . . . , n.
Example 4.5.4 Find the rounding error bound in linear interpolation, of ex for x ∈ [0, 1],
when data were rounded to four digits.

Solution
1 −4
|P1 (x) − P1 ∗ (x)| ≤ 10 (|L0 (x)| + |L1 (x)|)
2
Now, L0 (x) = 1 − x
and L1 (x) = x
and so on [0, 1], |L0 (x)| + |L1 (x)| = 1 − x + x = 1.
1 −4
Thus, |P1 (x ) − P1? (x )| ≤ 10 ,
2
which says that the effect of rounding errors in the data, on P1 (x) maintains the same maximum
magnitude.

Page 94 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

Example 4.5.5 The table below gives the tabulated values of the probability integral,
r ˆ x
2 −t2
I= e 2 dt
π 0
Use linear interpolation to find the value of I when x = 1.125. Estimate also the error bound
on the truncation error over [1, 1.25]
x 1.00 1.05 1.10 1.15 1.20 1.25
I 0.683 0.705 0.729 0.750 0.770 0.789
Solution
Note 4.5.2
At the tabular points x0 , x1 , ..., xn , the truncation error is zero, hence we can write
n
Y
f (x) = Pn (x) + R (x − xi ) (R constant depending on x)
i=0

Let,
n
Q
(t − xi )
i=0
F (x) = f (t) − Pn (t) − {f (x) − Pn (x)} Qn (∗)
(x − xi )
i=0

then
F (xi ) = f (xi ) − Pn (xi ) = 0, i = 0, 1, 2, . . . n
Further
F (x) = f (x) − Pn (x) − f (x) + Pn (x) = 0, ∀ x ∈ (a, b), x 6= xi ,
necessarily and i = 0, 1, 2, . . . n.
Thus, F (t) vanishes in (a, b) at n + 2 distinct points. On applying Rolle’s theorem repeatedly
we conclude that there exists c ∈ (a, b) such that, F (n+1) (c) = 0.
Differentiating equation (∗) (n + 1) times and putting t = c, gives,
n
Y f (n+1) (c)
f (x) = Pn (x) + (x − xi )
i=0
(n + 1)!

For linear interpolation, n = 1, so,

x − x1 x − x0
P1 (x) = f (x0 ) + f (x1 )
x0 − x1 x1 − x0
with
x0 = 1.10, x1 = 1.15 ⇒ f (x0 ) = .729, f (x1 ) = 0.750
thus ,
(1.125 − 1.15)(0.729) (1.125 − 1.10)(0.75)
P1 (1.125) = +
−0.05 0.05
= 0.5(0.729 + 0.750) = 0.7395

Truncation error
1
Y f (2) (c)
(x − xi )
i=0
2!

D.W-Ddumba numerical analysis i Page 95 of 133

 4.5. LAGRANGE INTERPOLATION

2
Now (x − x0 )(x − x1 ) is minimum at x = x0 +x 2
1
and is equal to (x1 −x
4
0)
in magnitude.
r √
2 −x2 2 − x2
f 0 (x) = e 2 , f 00 (x) = − xe 2
π π
√ √
2 −x2 2 2 −x2
f 000 (x) = − e 2 + xe 2
π π
and is 0 if x = ±1. Thus , √
00 2
|f (x)| ≤ √
π e
A bound on the truncation error is given by,
r
(0.05)2 2
' 0.00009
8π e

Exercise 4.2
1. Compute a bound on the truncation error for ex on [1, 1.4] when fourth degree polynomial
is used to interpolate ex .

2. Obtain error bounds for both linear and quadratic interpolation for sin hx over the interval
[1.90, 2.10]
4
3. A table of values for x 12−x is constructed for 0 ≤ x ≤ 1 in such a way that the error in
linear interpolation does not exceed  if rounding errors
√ were neglected. Show that for
uniform spacing h, then h does not exceed the value 2 2.

4. Find the rounding error bound when quadratic interpolation, of ex for x ∈ [0, 1], when
data were rounded to four digits.

Note 4.5.3 Use Lagrange interpolation to find an appropriate function passing through the
given points. Sketch a graph of this function based only on the given points and what you
think the curve must be. Compare your sketch with the graph created by graphing technology.
(a) A linear function passing through the points (−1, 3) and (2, 1).

(b) A quadratic function passing through the points (−1, 3), (0, 2), and (2, 1).

(c) A cubic function passing through the points (−1, 3), (0, 2), (1, 5), and (2, 1).

(d) A quartic (fourth degree) polynomial function passing through (−2, 4), (−1, 3), (0, 2), (1, 5),
and (2, 1).

Note 4.5.4 Finding a quadratic function that resembles other functions: By choosing three
noncolinear points on any curve we can use Lagrange’s interpolation to find a parabola that
passes through those points. For each of the following functions find a parabola that passes
through the graphs of the functions when at points with the indicated first coordinates. Use
graphing technology to draw a sketch of the function and the quadratic function you find.
Discuss how you might use the function you find to estimate the value of the given function.
(a) f (x) = x5 − 4x3 + 2; x = 0, 1, 2.
√
(b) f (x) = x; x = 0, 1, 4.

Page 96 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

(c) f (x) = 2x ; x = −1, 0, 1.

(d) f (x) = 2x ; x = 0, 1, 2.
(e) f (x) = sin π2 x ; x = 0, 1, 2.

(f) f (x) = cos π2 x ; x = −1, 0, 1 .

Describe a general procedure for finding a polynomial function of degree n that passes through
n + 1 given points with distinct first coordinates.
Note 4.5.5 Lagrange’s interpolation formula has the disadvantage that the degree of the
approximating polynomial must be chosen at the outset; an alternative approach is discussed in
the next Step. Thus, Lagrange’s formula is mainly of theoretical interest for us here; in passing,
we mention that there are some important applications of this formula beyond the scope of this
book - for example, the construction of basis functions to solve differential equations using a
spectral (discrete ordinate) method.
Note 4.5.6 Given the data below,
x -3 1 4 5 7
f(x) -28 4 28 36 52
Use linear interpolation to approximate
(a) f (3).
(b) f (x) if x = 5.
(c) x if f (x) = 12.
(d) the root/solution of f (x). f (x) = 8x − 4

Example 4.5.6 Use Lagrange’s polynomial to approximate f (2) for

x0 x1
x 1 3
f(x) 5 21
f (x) = 8x − 3
(x − x1 ) 1
L0 (x) = = (3 − x)
(x0 − x1 ) 2
(x − x0 ) 1
L1 (x) = = (x − 1)
(x1 − x0 ) 2

P1 (x) = L0 (x)f (x0 ) + L1 (x)f (x1 )

1 1
= (3 − x)(5) + (x − 1)(21)
2 2
1
= [15 − 5x + 21x − 21]
2
1
= [16x − 6]
2
≈ 8x − 3
f (2) = 8(2) − 3 = 13

D.W-Ddumba numerical analysis i Page 97 of 133

 4.5. LAGRANGE INTERPOLATION

Example 4.5.7 The Lagrange basis polynomials (n=N-1) on the data to approximate f (3.5)

x0 x1 x2
x 1 2 4
f(x) 3 2 1

(x − x1 )(x − x2 ) 1
L0 (x) = = (x − 2)(x − 4)
(x0 − x1 )(x0 − x2 ) 3
(x − x0 )(x − x2 ) 1
L1 (x) = = − (x − 1)(x − 4)
(x1 − x0 )(x1 − x2 ) 2
(x − x0 )(x − x1 ) 1
L2 (x) = = (x − 1)(x − 2)
(x2 − x0 )(x2 − x1 ) 6

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

1 1 1
= (x − 2)(x − 4)(3) − (x − 1)(x − 4)(2) + (x − 1)(x − 2)
3 2 6
1
= (x − 2)(x − 4) − (x − 1)(x − 4) + (x − 1)(x − 2)
6
1 2 3 1
= x − x+4+
6 2 3
≈ 0.1667x2 − 1.5x + 4.3333

f (3.5) ≈ 0.1667(3.5)2 − 1.5(3.5) + 4.3333 ≈ 1.125075

Note 4.5.7 Lagrange error for n = 2,

f (3) (ξ)
f (x) = P2 (x) + (x − x0 )(x − x1 )(x − x2 ), ξ ∈ [a, b], where {x0 , x1 , x2 } ⊆ [a, b]
3!

Note 4.5.8 The identity

n
X
Lk (x) = 1
k=0

(established by setting f (x) = 1) may be used as a check. Note also that with n = 1 we recover
the linear interpolation formula:

P1 (x) = L0 (x)f (x0 ) + L1 (x)f (x1 )

(x − x1 ) (x − x0 )
= f (x0 ) + f (x1 )
(x0 − x1 ) (x1 − x0 )
(x − x0 )
= f (x0 ) + [f (x1 ) − f (x0 )]
(x1 − x0 )

Example 4.5.8 Use Lagrange’s interpolation formula to find the interpolating polynomial
P3 (x) through the points (0, 3), (1, 2), (2, 7), and (4, 59) and then find the approximate value of
P3 (3).

Page 98 of 133 D.W-Ddumba numerical analysis i

 4.5. LAGRANGE INTERPOLATION

x0 x1 x2 x3
x 0 1 2 4
f(x) 3 2 7 59

Here n = 3 = 4 − 1
The Lagrange coefficients are:

(x − x1 )(x − x2 )(x − x3 ) (x − 1)(x − 2)(x − 4) 1

L0 (x) = = = − (x3 − 7x2 + 14x − 8)
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (0 − 1)(0 − 2)(0 − 4) 8
(x − x0 )(x − x2 )(x − x3 ) (x − 0)(x − 2)(x − 4) 1
L1 (x) = = = (x3 − 6x2 + 8x)
(x1 − x0 )(x1 − x2 )(x1 − x3 ) (1 − 0)(1 − 2)(1 − 4) 3
(x − x0 )(x − x1 )(x − x3 ) (x − 0)(x − 1)(x − 4) 1
L2 (x) = = = − (x3 − 5x2 + 4x)
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (2 − 0)(2 − 1)(2 − 4) 4
(x − x0 )(x − x1 )(x − x2 ) (x − 0)(x − 1)(x − 2) 1
L3 (x) = = = (x3 − 3x2 + 2x)
(x3 − x0 )(x3 − x1 )(x3 − x2 ) (4 − 0)(4 − 1)(4 − 2) 24

(The student should verify that [L0 (x) + L1 (x) + L2 (x) + L3 (x) = 1]
Hence, the required polynomial is

3 2
P3 (x) = − (x3 − 7x2 + 14x − 8) + (x3 − 6x2 + 8x)
8 3
7 3 59
− (x − 5x2 + 4x) + (x3 − 3x2 + 2x)
4 24

1
−9x3 + 63x2 − 126x + 72 + 16x3 − 96x2 + 128x − 42x3 + 210x2 − 168x + 59x3 − 177x2 + 118x
 
= 24

1 
−9x3 + 63x2 − 126x + 72 + 16x3 − 96x2 + 128x − 42x3 + 210x2 − 168x + 59x3 − 177x2 + 118x

=
24
1 
24x3 + 0x2 − 48x + 72

=
24
= x3 − 2x + 3

Consequently f (3) ≈ P (3) = 33 − 2(3) + 3 = 24, However, note that, if the explicit form of
the interpolating polynomial were not required, one would proceed to evaluate P3 (x) for some
value of x directly from the factored forms of Lk (x). Thus, in order to evaluate P3 (3), one has

(x − x1 )(x − x2 )(x − x3 ) (3 − 1)(3 − 2)(3 − 4) 1

L0 (3) = = = . etc.
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (0 − 1)(0 − 2)(0 − 4) 4

Exercise 4.3 Given that f (−2) = 46, f (−1) = 4, f (1) = 4, f (3) = 156, and f (4) = 484, use
Lagrange’s interpolation formula to estimate the value of f (0).

Example 4.5.9 Use Lagrange interpolation polynomial for the data below

x0 x1 x2 x3
x 1 2 4 5
f(x) 3 8 54 107

D.W-Ddumba numerical analysis i Page 99 of 133

 4.5. LAGRANGE INTERPOLATION

to show that P (x) = x3 − x2 + x + 2 The Lagrange coefficients are:

(x − x1 )(x − x2 )(x − x3 ) (x − 2)(x − 4)(x − 5) 1 3

L0 (x) = = =− (x − 11x2 + 38x − 40)
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (1 − 3)(1 − 4)(1 − 5) 12
(x − x0 )(x − x2 )(x − x3 ) (x − 1)(x − 4)(x − 5) 1
L1 (x) = = = (x3 − 10x2 + 29x − 20)
(x1 − x0 )(x1 − x2 )(x1 − x3 ) (2 − 1)(2 − 4)(2 − 5) 6
(x − x0 )(x − x1 )(x − x3 ) (x − 1)(x − 2)(x − 5) 1
L2 (x) = = = − (x3 − 8x2 + 17x − 10)
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (4 − 1)(4 − 2)(4 − 5) 6
(x − x0 )(x − x1 )(x − x2 ) (x − 1)(x − 2)(x − 4) 1
L3 (x) = = = (x3 − 7x2 + 14x − 8)
(x3 − x0 )(x3 − x1 )(x3 − x2 ) (5 − 1)(5 − 2)(5 − 4) 12

(The student should verify that [L0 (x) + L1 (x) + L2 (x) + L3 (x) = 1] Hence, the required
polynomial is
1 3 1
P3 (x) = − (x − 11x2 + 38x − 40)[3] + (x3 − 10x2 + 29x − 20)[8]
12 6
1 3 1
+ (x − 8x2 + 17x − 10)[54] + (x3 − 7x2 + 14x − 8)[107]
6 12

1

= 12 −3x3 + 33x2 − 114x + 120 + 16x3 − 160x2 + 464x − 320 − 108x3 + 864x2 − 1836x + 1080 + 107x3 − 749x2 + 1498x −

1 
= −3x3 + 33x2 − 114x + 120 + 16x3 − 160x2 + 464x − 320 − 108x3 + 864x2 − 1836x + 1080 + 107x
12
1 
12x3 − 12x2 + 12x + 24

=
12
= x3 − x2 + x + 2

Consequently f (3.5) ≈ P (3.5) = (3.5)3 − (3.5)2 + (3.5) + 2 = 36.125,

Page 100 of 133 D.W-Ddumba numerical analysis i

Chapter 5

Numerical Differentiation

5.1 Why Numerical techniques for finding derivatives

A large number of physical problems involve functions which are known only through experi-
mental measurements. In order to find derivatives of such functions one is forced to use methods
based on the available discrete data and the methods are purely numerical.
Sometimes, some functions with known analytical expressions, their derivatives may be so
complex and computationally involved in their evaluation that one must resort to numerical
methods, though less accurate, but have less involved expressions.
Throughout numerical weather prediction, you often need to calculate the gradient of a function
at a number of points. As we rarely know the equation which defines the function, we need to
calculate numerically an estimate of the gradient. This method requires the use of finite differ-
ence schemes. In this module we will be considering finite difference schemes to approximate
the gradient of a function in one variable. The problem can be stated more specifically in these
terms: how to estimate the gradient of a function in f (x), at a point a.
Numerical differentiation is the process of finding the numerical value of a derivative of a given
function at a given point. In general, numerical differentiation is more difficult than numer-
ical integration. This is because while numerical integration requires only good continuity
properties of the function being integrated, numerical differentiation requires more complicated
properties such as Lipschitz classes.

5.2 Analytic definition of a derivative as compared to a

numerical definition
Analytically, we define the derivative of f (x) at x = a denoted f 0 (x) as the limiting process,
f (x + h) − f (x)
f 0 (x) = lim .
h→0 h
However, a numerical approximation to f 0 (x) involves a difference such as
f (x + h) − f (x)
Φ(h) =
h
In geometrical terms, this can be represented as in Figure 5.1.
Note 5.2.1 f 0 (x) is the slope of the tangent T and, Φ(h) is the slope of the line L. As
h → 0, Φ(h) → f 0 (x) suggesting that h should be small for a good approximation. However, in
computing Φ(h) the two terms f (x) and f (x + h) are close in value. So there is likely to be a
loss of significant digits.

101
 5.3. FORWARD DIFFERENCE APPROXIMATION

Figure 5.1: Geometrical interpretation of a derivative

5.3 Forward difference approximation

The forward difference approximation for a derivative is given by,
f (x + h) − f (x)
f 0 (x) =
h
The formula can geometrically be represented as in Figure (5.2)
What the above expression actually calculates is the gradient of the line which intersects
the points x, f (x) and x + h, f (x + h) as you can see from the figure above. In the figure, the
function is shown with a red line (a curve). The blue line (a lower straight line) represents the
approximation to the gradient and the black line (the upperline) is the actual gradient of the
function at x = a .

Example 5.3.1
For the function f (x) = x2 , approximate f 0 (2) with step length
(i) h = 0.1

(ii) h = 0.01

(iii) h = 0.001.
Compare your results with the analytic/exact value of f 0 (2).
Solution
(i) Since
f (x + h) − f (x)
f (x) '
h
But f (x) = x2 and x = 2, h = 0.1

therefore f(x + h) = (2 + 0.1)2 = (2.1)2 = 4.41

Page 102 of 133 D.W-Ddumba numerical analysis i

 5.3. FORWARD DIFFERENCE APPROXIMATION

Figure 5.2: Geometrical representation of the forward difference approximation

and
f (x) = 22 = 4
therefore
4.41 − 4 0.41
f 0 (2) ' = = 4.1
0.1 0.1
However, the exact value of
f 0 (2) = (2)(2) = 4
Since f 0 (x) = 2x. This yields an error of 0.1.

(ii) Since h = 0.01, therefore

f (x + h) = (2.01)2 = 4.0401
and
f (x) = 22 = 4
therefore
4.0401 − 4 0.0401
f 0 (2) ' = = 4.01
0.01 0.01
The error committed is 0.01 i.e smaller than in part (i).

(iii) For h = 0.001,

f (x + h) = (2.001)2 = 4.004001
and
f (x) = 22 = 4
therefore
4.004001 − 4 0.004001
f 0 (2) ' = = 4.001
0.001 0.001
The error committed in this case is 0.001.

D.W-Ddumba numerical analysis i Page 103 of 133

 5.3. FORWARD DIFFERENCE APPROXIMATION

5.3.1 Analytical derivation of the forward difference approximation

We derive the forward difference approximation analytically from Taylor series expansion of
f (x + h). By Taylor series we have,
 h2 00 3

 2
f (x) + h6 f (3) (x) + . . . 1(a)



0
f (x + h) = f (x) + hf (x) +




 h2 00
2
f (c1 ) x ≤ c1 ≤ x + h 1(b)

If we rearrange equation 1(b), we get,

(x + h) − f (x) h 00
f 0 (x) = − f (c1 ) (5.1)
h 2
Equation (5.1) is the forward difference form which we can write in better notation as:

f (x + h) − f (x)
f 0 (x) ' + Truncation Error
h

Note 5.3.1

(i) Error 0(h) means that the error → 0 (goes to zero) near x = a at the rate of Ah where
A is a constant (i.e. halving h, halves the error).

(ii) The magnitude of the truncation error i.e

h
Etruc = − f 0 (c1 ).
2
Indeed | − h2 f 0 (c1 ) is the bound on the truncation error.

Page 104 of 133 D.W-Ddumba numerical analysis i

 5.4. BACKWARD DIFFERENCE APPROXIMATION

5.4 Backward difference approximation

The backward difference approximation for the derivative is given by,

f (x) − f (x − h) h 00
f 0 (x) = + f (c2 )
h 2
The formula can be geometrically represented as in Figure (5.3)

Figure 5.3: Geometrical representation of the backward difference approximation

Example 5.4.1 For the function f (x) = x2 , approximate f 0 (2) using the backward difference
approximation with step length

(i) h = 0.1,

(ii) h = 0.01,

(iii) h = 0.001.

Solution

(i) Since the backward difference approximation is,

f (x) − f (x − h)
f 0 (x) '
h
But f (x) = x2 , x = 2 and h = 0.1, therefore

f (x) = f (2) = 22 = 4

f (x − h) = f (2 − 0.1) = f (1.9) = (1.9)2 = 3.61

D.W-Ddumba numerical analysis i Page 105 of 133

 5.4. BACKWARD DIFFERENCE APPROXIMATION

therefore
4 − 3.6100 0.3900
f 0 (2) '= = = 3.900
0.1 0.1
Since the exact value of
f 0 (2) = (2)(2) = 4,
therefore the absolute error committed by using the numerical formula is 0.1. This is not
very bad.

(ii) Since h = 0.01 and

f (x) = f (2) = 22 = 4
f (x − h) = f (2 − 0.01) = f (1.99) = (1.99)2 = 3.9601,
therefore
4 − 3.9601 0.0399
f 0 (2) ' = = 3.99
0.01 0.01
Error committed is 0.01. This is better than the previous.

(iii) Since h = 0.001 and f (x) = 4

f (x − h) = f (2 − 0.001) = f (1.999) = (1.999)2 = 3.996001,

therefore
4 − 3.996001 0.003999
f 0 (2) ' = = 3.999
0.001 0.001
Error committed is 0.00001. This error is smaller than any of the previous two. Infact the zero
error means the exact value of the derivative is generated.

Note 5.4.1 We note that the smaller h is the better is the numerical approximate to the
derivative.

5.4.1 Analytical derivation of the backward difference approxima-

tion
Using the Taylor series expansion
h2 00 h3 (3)


 2
f (x) − 6
f (x) + . . . (?)



f (x − h) = f (x) − hf 0 (x) +




h2 00
x ≤ c1 ≤ x + h

2
f (c2 ) (??)
Rearranging equation (??) i.e

h2 00
f (x − h) = f (x) − hf 0 (x) + f (c2 ),
2
we get,
f (x) − f (x − h) h 00
f 0 (x) = + f (c2 )
h 2
with x − h ≤ c2 ≤ x
Which is the backward difference approximation. The approximation could also be written
as,

Page 106 of 133 D.W-Ddumba numerical analysis i

 5.4. BACKWARD DIFFERENCE APPROXIMATION

f (x) − f (x − h)
f 0 (x) = + Etrunc
h
or
f (x) − f (x − h)
f 0 (x) = + 0(h)
h

D.W-Ddumba numerical analysis i Page 107 of 133

 5.5. THE CENTRAL DIFFERENCE APPROXIMATION

5.5 The Central difference approximation

We state the central difference approximation to the first derivative as;

0 f (x + h) − f (x − h) h2 3
f (x) = + f (c)
2h 3!
The formula is geometrically represented as seen in the Figure (figure:central)

Figure 5.4: Geometrical interpretation of the central difference approximation

The slope of the line AB is the central difference approximation

f (x + h) − f (x − h)
f 0 (x) ' .
2h
From Figure 5.4, it is clear that the gradient of AB is closer to the gradient of the tangent at
x = a. That is, the two lines are almost parallel. This is expected to be a better approximation
than the forward and the backward. Indeed it is as we shall see in the following example.

Example 5.5.1
Approximate f 0 (2) for f (x) = x2 using the central difference approximation with step size

(i) h = 0.1

(ii) h = 0.01

(iii) h = 0.001

Solution

(i) Since a = 2, h = 0.1, but

f (x + h) − f (x − h)
f 0 (x) '
2h
f (2.1) − f (1.9)
=
2(0.1)
(2.1)2 − f (1.9)2
= = 4.000000000
0.2
Since the exact value is,
f 0 (2) = 2(2) = 4

Page 108 of 133 D.W-Ddumba numerical analysis i

 5.5. THE CENTRAL DIFFERENCE APPROXIMATION

therefore the error is zero to 9 places of decimals.

This could not be achieved with a forward or backward formula using same step size.
This is because the central formula is of higher order i.e 0(h2 )

(ii) with h = 0.01,

f (2.01) − f (1.99)
f 0 (q) '
2(0.01)

(2.01)2 − (1.99)2
= = 4.000000000.
0.02
This also generates a zero error.

(iii) with h = 0.001,

f (2.01) − f (1.999)
f 0 (2) '
2(0.001)

(2.001)2 − (1.999)2
= = 4.00000000
0.002

Which also generates a zero error.

5.5.1 Analytical derivation of the central difference approximation

Using Taylor series expansion of f (x + h) and f (a − h) about x = a, we have;

h2 00
0 h3 (3)
f (x + h) = f (x) + hf (x) + f (x) + f (c1 ) (5.2)
2 6

and
h2 00 h3
f (x − h) = f (x) − hf 0 (x) + f (x) − f (3) (c2 ) (5.3)
2 6

equation (5.2) minus equation (5.3) we get;

h3 (3)
f (x + h) = f (x − h) = 2f 0 (x)h + (f (c1 ) − f (3) (c2 )) (5.4)
6
Now,
f (3) (c1 ) + f (3) (2)
= f (3) (c3 ), c1 ≤ c3 ≤ c2
2
This is derived from the intermediate value theorem.
So
f (x + h) − f (x − h) h2 (3)
f 0 (x) = + f (c3 )
2h 3!
i.e, the error is of order 0(h2 ).

D.W-Ddumba numerical analysis i Page 109 of 133

 5.5. THE CENTRAL DIFFERENCE APPROXIMATION

Exercise 5.1 Derive the common approximations below

f (x + h) − f (x)

f 0 (x) ≈ (a) 

h 





f (x) − f (x − h)


f 0 (x) ≈ (b)



h 





f (x + h) − f (x − h)

f 0 (x) ≈ (c) (5.5)
2h 




−f (x + 2h) + 4f (x + h) − 3f (x)


f 0 (x) ≈

(d) 

2h 





−f (x + 2h) + 4f (x + h) − 4f (x − h) + f (x + 2h)


f 0 (x) ≈

(e) 
4h

Exercise 5.2
1. For f (x) = ex , approximate f 0 (1) using forward, backward and central difference formulae
with,

(i) h = 0.1
(ii) h = 0.01
(iii) h = 0.001

Compare the results with the analytic/exact value of f 0 (1). Comment on your results.

2. Let f (x) be given by the table below. The inherent round off error has the bound

|ek | ≤ 5 × 10−5

use the rounded values in your calculations.

x 1.100 1.190 1.199 1.200 1.201 1.210

f (x) 0.4536 0.3717 0.3633 0.3624 0.314 0.3630

(a) Find approximation for f 0 (1.2) using the central difference formula with h = 0.01,
and h = 0.001.
(b) Compare with f 0 (1.2) = − sin(1.2) = −0.932.
(c) Find the total error bound for the three cases in part (a).

3. For f (x) = cos x, Use the forward, backward and central difference formulae to approxi-
mate f 0 (0.8) using h = 0.001. Compare your results with the analytic values.[Hint: either
all in radians or all in degrees, Ans= -0.717]

4. Repeat question three using the formula

f1 − f−1
f 0 (x0 ) ' (5.6)
2h
Compare your results with analytic results. What can you say about the order of the
truncation error in (5.6) in comparison to the forward, backward and central difference
formulae.

Page 110 of 133 D.W-Ddumba numerical analysis i

 5.6. COMPARISION

5.6 Comparision
It is clear that the central difference gives a much more accurate approximation of the derivative
compared to the forward and backward differences. Central differences are useful in solving
partial differential equations. If the data values are available both in the past and in the future,
the numerical derivative should be approximated by the central difference.

5.7 The second derivative approximation

The most commonly used approximation is of a central difference form which we obtain from
the following Taylor series approximations earlier considered in the previous lecture. Consider
again the expansions,

h2 00
0 h3 (3) h4 (4)
f (x + h) = f (x) + hf (x) + f (x) + f (x) + f (c)) (5.7)
2 6 24
and

h2 00 h3 h4
f (x − h) = f (x) − hf 0 (x) + f (x) − f (3) (x) + f (4) (c2 ) (5.8)
2 6 24
addition of equations (5.7) and (5.8) gives,

h4 (4)
 
2 00 4
f (x + h) + f (x − h) = 2f (x) + h f (x) + f (c1 ) + f (c2 ) (5.9)
12
and using,
f (4) (c1 ) + f (4) (c2 )
= f (4) (c3 )
2
with c1 ≤ c3 ≤ c2 . This is derived from the intermediate value theorem and making f 00 (x) the
subject in equation (5.9), we get,

00 f (x − h) − 2f (x) + f (x + h) h2 (4)
f (x) = + f (c3 ) (5.10)
h2 6
or

f (x − h) − 2f (x) + f (x + h)
f 00 (x) = + 0(h2 ) (5.11)
h2
Equation (5.11) is a second order approximation for the second order derivative. The formula
is handy for approximating second order derivatives.

Example 5.7.1 Approximate f 00 (1) for the function f (x) = x3 with

(i) h = 0.100

(ii) h = 0.010

(iii) h = 0.001.

Solution

D.W-Ddumba numerical analysis i Page 111 of 133

 5.7. THE SECOND DERIVATIVE APPROXIMATION

(i) Since,
f (x − h) − 2f (x) + f (x + h)
f 00 (x) '
h2
But x = 1, h = 0.1, therefore
f (0.9) − 2f (1) + f (1.1)
f 00 (1) '
(0.1)2
(0.9)3 − 2(1)3 + (1.1)3
=
0.01
= 6.00000000000005

However, the exact value is f 00 (x) = 6x, therefore f 00 (1) = 6. Hence error committed is
just 5.0 × 10−14 . This is really small.
(ii) With h = 0.01
f (0.99) − 2f (1) + f (1.01)
f 00 (1) '
(0.01)2
(0.99)3 − 2 + (1.01)3
=
0.0001
= 6.00000000000000
Since the exact value is 6, the error committed is zero to 14 decimal places.
(iii) This also generates a zero error, since,
f (0.999) − 2f (1) + f (1.001)
f 00 (1) '
(0.001)2
(0.999)3 − 2 + (1.001)3
=
0.000001
= 6.00000000000000

5.7.1 Error analysis in numerical differentiation

Addition of a rounding error term
When we use a formula such as
f (x + h) − f (x)
f 0 (x) = + Etrun (h2 ) (5.12)
h
in general the difference f (x+h)−f
h
(x)
will be evaluated with a rounding error. For example if f
is evaluated to four decimal places with h = 0.1, there is a possible rounding error of
2 × 0.00005
= 0.001
0.1
in the difference
f (x + h) − f (x)
.
h
This is because when a number is evaluated to n- decimal places, the maximum absolute error
committed is 21 10−n . Thus maximum error in
1 −4
f (x + h) is 10 = 0.00005
2
Page 112 of 133 D.W-Ddumba numerical analysis i
 5.7. THE SECOND DERIVATIVE APPROXIMATION

and maximum absolute error in

1 −4
f (x) is 10 = 0.00005.
2
But the error in the difference is a sum of the absolute errors in the individual numbers. So
the maximum error in
1
f (x + h) − f (x) is 2( 10−4 ) or 2(0.00005)
2
Thus, the error in
f (x + h) − f (x) (0.00005)
is 2
h 0.1
assuming h = 0.1 is exact. Thus equation (5.12) becomes
y1 − y0
f 0 (x) = + Etruc + Eround
h
where, y1 , y0 are the rounded values of f (x + h) and f (x) respectively.

Optimum step size

We now have that:
Total error = Etrunc + Eround .
Thus there is no point in choosing h so that |Etrunc | is small. If |Eround | is large since the benefit
is swamped.

Example 5.7.2 Suppose

f (x + h) = y1 + e1
and
f (x) = y0 + e0
and |e1 |, |e0 | < e (small number). Find the optimum choice for h when using the approximation

f (x + h) − f (x)
f 0 (x) '
h
Solution
For the approximation,
f (x + h) − f (x)
f 0 (x) '
h
Maximum absolute error in f (x + h) is e and that in f (x) is also e. So

e+e 2e
|Eround | < =
h h
Now
h 00 h
|Etrunc | =
|f (c1 )| ≤ M2
2 2
00
say, where M2 = max |f (x)|, a ≤ x ≤ a + h Thus,

2e h
|Total error| ≤ + M2 .
h 2
D.W-Ddumba numerical analysis i Page 113 of 133
 5.7. THE SECOND DERIVATIVE APPROXIMATION

The optimum choice for h is

d
|total error| = 0
dh
2e M2
that is, − 2 + =0
h 2
r
e
giving h = 2
M2

Exercise 5.3

1. For f (x) = x2 + x, approximate f 00 (1) using a step length,

(i) h = 0.4
(ii) h = 0.04
(iii) h = 0.004

2. Given the function f (x) = ex , approximate f 00 (x) at x = 2 using step length h of magni-
tude.

(i) 0.1
(ii) 0.01
(iii) 0.001

3. Let f (x) = cosx. Use the formula considered in this Lecture for approximating f 00 (x),
with h = 0.01 to calculate approximations for f 0 (0.8). Compare with the true value.

4. Given that
−f2 + 16f1 − 30f0 + 16f−1 − f−2
f 00 (x0 ) '
12h2
Use the formula to approximate f 00 (1) for f (x) = ex using h = 0.5. Compare your answer
with the analytic answer and the answer obtained when using formula equation (5.11).
What do you notice and suggest?

5. Show that for the central difference formula

f (x + h) − f (x) h2 (3)
f 0 (x) = + f (c),
2h 3!
e Mh2
(i) |Total error| ≤ h
+ 6
3e 3 1
(ii) and Optimum h = ( M )
(iii) for f (x) = cosx, x = 0.8 and h = 0.0001, Show that, | Etrunc | ' 0.1 × 10−8 . And if
e = 0.5 × 10−9 , Show that, |Eround | < 0.5 × 10−5 and Optimum h = 0.0011.

Page 114 of 133 D.W-Ddumba numerical analysis i

Chapter 6

Finite difference operators

6.1 Introduction
The lecture defines and inter-relates the most common finite difference operators i.e. forward,
backward, shift, central and the averaging operator. Use of the operators to prove finite differ-
ence identities is also done.

6.2 Finite differences

We consider a function f (x) known for a certain set of equally spaced values of x such that

x = x0 + rh, r = 0, 1, . . . , n

and h > 0. This generates a set of n + 1 pairs i.e. (xr , fr ).

The xr are called pivotal points and fr = f (xr ) are the pivotal values. The pairs can be
written as;

x0 f (x0 )
x1 f (x1 )
x2 f (x2 )
· ·
· ·
· ·
xr f (xr )
We can form differences between consecutive pivotal values i.e. fr+1 − fr where

xr+1 = xr+h for r = 0, 1, 2, . . .

The differences fr+1 − fr are called the first differences. We similarly define higher order
differences.

Example 6.2.1

Masenge (1989) used in his book the example of tabulating the function f (x) = ln(1 + x)
in the interval [2.00, 2.10]. This was tabulated as,

115
 6.3. FINITE DIFFERENCE OPERATORS

Differences
xr fr 1st 2nd 3rd
2.00 1.098612
0.006645
2.02 1.105257 −0.000044
0.00601 0.00000
2.04 1.111858 −0.000044
0.006557 0.000002
2.06 1.118415 −0.000042
0.006515 −0.000001
2.08 1.124930 −0.000043
0.006472
2.10 1.31402

Definition 6.2.1
We say that differences are converging if elements of higher order differences decrease rapidly in
magnitude. Such a behavior is common from tables of well-behaved functions like polynomials.

6.3 Finite difference operators

The most common finite difference operators are; the forward difference operator, the back-
ward difference operator, the central difference operator, the averaging operator and the shift
operator.

6.3.1 The forward difference operator

The forward difference operator is denoted by 4 and defined by the difference equation;

4fr = fr+1 − fr

Thus, the 1st difference column in a difference table consists of the elements,

4f0 , 4f1 , 4f2 . . . . . .

The second differences are got by differentiating the first differences i.e

4 fr+1 − 4fr (6.1)

But we can easily show that 4 is a linear operator i.e.

4(αfr + βgr ) = α 4 fr + β 4 gr

Thus equation (6.1), becomes,

4fr+1 − 4fr = 4(fr+1 − fr ) = 4(4fr ).

If we write 42 fr for 4(4fr ), then the second differences are the quantities,

42 f 0 , 42 f 1 , 42 f 2 , . . . . . .

We similarly get third differences as consisting of the elements.

4(42 fr ) = 43 fr , r = 0, 1, 2, . . .

Page 116 of 133 D.W-Ddumba numerical analysis i

 6.3. FINITE DIFFERENCE OPERATORS

i.e 42 fr+1 − 42 fr = 42 (fr+1 − fr )

= (4fr )
= 43 fr
Generally, we have that the nth order differences consist of the elements/quantities,

4n fr = 4n−1 fr+1 − 4n−1 fr

= 4(4n−1 fr )
Thus 4k fr involves information at the pivoted points xr , xr+1 , . . . , xr+k . Thus,

4f0 involves information at x0 and x1

42 f0 involves information at x0 , x1 , andx2
43 f0 involves information at x0 , x1 , x2 and x3
and so on. In fact this is how the name forward comes about.

Differences as related to ordinates

By definition of the forward difference operator, we have:

4fr = fr+1 − fr

42 fr = 4fr+1 − 4fr
= (fr+2 − fr+1 ) − (fr+1 − fr )
= fr+2 − 2fr+1 + fr
43 fr = 42 fr+1 − 42 fr
= (fr+3 − 2fr+2 + fr+1 ) − (fr+2 − 2fr+1 + fr )
= fr+3 − 3fr+2 + 3fr+1 − fr .
Similarly we can show that,

44 fr = fr+4 − 4fr+3 + 6fr+2 − 4fr+1 + fr

Try and show this on a rough paper before you continue.

6.3.2 The backward difference operator

This is denoted by 5 and defined by the difference equation,

5fr = fr − fr+1 .

With similar argument as for 4− operator, show that

5n fr = 5n−1 fr − 5n−1 fr−1

before you continue reading.

We also note that this operator has a nice property of linearity i.e

5(αfr + βgr ) = α 5 fr + βgr

where α and β are real scalars.

D.W-Ddumba numerical analysis i Page 117 of 133

 6.3. FINITE DIFFERENCE OPERATORS

6.3.3 The shift operator

The shift operator is denoted by E and defined by

Efr = fr+1 .

Hence, by this definition we can define higher order shift operators as,

E(Efr ) = Efr+1 = fr+2 .

We usually denote E(Efr ) by E 2 fr . and,

E k fr = E(E k−1 fr ) = fr+k .

For k negative, we get what is known as the backward shift. For example E −1 fr is a
backward shift defined by,
E −1 fr = fr−1
We also note that k can be fractional.
1
Example 6.3.1 Find E 2 fr and E −5 fr .

Solution
By definitions of the forward and backward shifts i.e.

Efr = fr+1

and
E −1 fr = fr−1
we have,
1
E 2 fr = fr+ 1
2

and
E −5 fr = fr−5 .

Relation between 4 and E

From
4fr = fr+1 − fr
and
Efr = fr+1
combining the two we get,
4fr = Efr − [Link]

= (E − 1)fr
where 1 is the identity operator I, that is,

Ifr = fr .

Thus 4 = E − 1 and E = 4 + 1.

Page 118 of 133 D.W-Ddumba numerical analysis i

 6.3. FINITE DIFFERENCE OPERATORS

Example 6.3.2 Show that

4=E5.

Solution
Since
E 5 fr = E(5fr ) = E(fr − fr−1 )
= Efr − Efr−1 (because of linearity of E)
= fr+1 − fr (By definition of E).
And 4 fr = fr+1 − fr (By definition of 4)
Hence shown.

Differences in terms of pivoted values

From the relation
4k fr = (E − 1)k fr ,
We can easily show that
     
k k k k
4 fr = fr+k − fr+k−1 + fr+k−2 + . . .
0 1 2
 
k k!
Try and show this before you continue. Recall: = r!(k−r)! .
r
Example 6.3.3 Express 43 f0 in terms of pivoted values.

Solution
Since
     
k k k k
4 fr = fr+k − fr+k−1 + fr+k−2 + . . .
0 1 2

      
3 3 3 3 3
therefore 4 f0 = f3+0 − f3+0−1 + f3+0−2 − f3+0−3
0 1 2 3
       
3 3 3 3
= f3 − f2 + f1 − f0
0 1 2 3
= f3 − 3f2 + 3f1 − f0

Relation between 5 and E

From,
5fr = fr − fr−1
and
E −1 fr = fr−1 .
Combining the two, we have
5fr = fr − E −1 fr = (1 − E −1 )fr
⇒ 5 = 1 − E −1
and
E = (1 − 5)−1

D.W-Ddumba numerical analysis i Page 119 of 133

 6.3. FINITE DIFFERENCE OPERATORS

6.3.4 The central difference operator

The central difference operator is denoted by δ and defined by the equation

δfr = fr+ 1 − fr− 1 .

2 2

And
δ 2 fr = δfr+ 1 − δfr− 1
2 2

= (fr+1 − fr ) − (fr − fr−1 )

= fr+1 − 2fr + fr−1
With higher order differences defined by

δ k fr = δ k−1 fr+ 1 − δ k−1 fr− 1 .

2 2

Try to find the expressions for δ 3 fr and δ 4 fr before you continue reading.

Relation between δ and E

From,
δfr = fr+ 1 − fr− 1
2 2

and
1
E 2 fr = fr+ 1
2

and
1
E − 2 fr = fr− 1 ,
2

we get
1 1
δfr = (E 2 − E − 2 )fr
and hence conclude that,
1 1 1
δ = E 2 − E − 2 = E 2 (1 − E −1 ).

6.3.5 Averaging operator

The averaging operator is denoted by u and defined by the equation,
1
ufr = {fr+ 1 + fr− 1 }
2 2 2

Exercise 6.1 Prove the following finite difference identities.

1
(i) E − 2 4 = δ
(ii) 5 = E −1 4
1
(iii) 4 = E 2 δ
1
(iv) δ = E 2 5
1
(v) 5 = E − 2 δ

Exercise 6.2 Prove the finite difference identity

1 1
uδ = 4 E −1 + 4 .
2 2

Page 120 of 133 D.W-Ddumba numerical analysis i

 6.4. FINITE DIFFERENCE TABLES

Exercise 6.3 From the fact that 4 = (E − 1), we can write

4k fr = (E − 1)k fr .

Hence show that the forward difference can be written in terms of pivotal values by the formula,
     
k k k k
4 fr = fr+k − fr+k−1 + fr+k−2 + . . .
0 1 2

Hence express 45 f0 in terms of pivotal values.

Exercise 6.4 Prove the following identities

1
(i) 4 = uδ + 21 δ 2 = E 2 δ

(ii) 4 − 5 = 45 = δ 2

(iii) u(fr gr ) = ufr ugr + 14 δfr δgr

(iv) 4(fr gr ) = fr 4 gr + gr+1 4 fr and deduce that

4(fr2 ) = (fr + fr+1 ) 4 fr

gr 4fr −fr 4gr

(v) 5( fgrr ) = gr gr +1
and deduce that,

1 − 4 fr
4( )=
f −r fr fr+1

Exercise 6.5 Use mathematical induction to prove that,

4r fk = 5r fk+r = δ r fk+ r2 ,

6.4 Finite difference tables

The fact that the nth order differences may be expressed in terms of (n−1)th differences enables
us to construct finite difference tables as follows:

D.W-Ddumba numerical analysis i Page 121 of 133

6.4. FINITE DIFFERENCE TABLES

Table 6.1: Showing a table of forward differences.

f−3
4f−3
f−2 42 f−3
4f−2 43 f−3
f−1 42 f−2 44 f−3
4f−1 43 f−2 45 f−3
2 4
f0 4 f−1 4 f−2 45 f−3
4f0 43 f−1 45 f−2
2 4
f1 4 f0 4 f−1
3
4f1 4 f0
2
f2 4 f1
4f2
f3

Table 6.2: Showing a table of backward differences.

f−3
5f−2
f−2 52 f−1
5f−1 53 f 0
f−1 52 f 0 54 f1
3
5f0 5 f1 55 f 2
f0 52 f 1 54 f2 56 f 3
5f1 53 f 2 55 f 3
2 4
f1 5 f2 5 f3
3
5f2 5 f3
2
f2 5 f3
5f3
f3

Page 122 of 133 D.W-Ddumba numerical analysis i

 6.4. FINITE DIFFERENCE TABLES

Table 6.3: Showing a table of central differences.

f−3
δf−2 12
f−2 δ 2 f−2
δf−1 1 δ 3 f−1 1
2 2
f−1 δ 2 f−1 δ 4 f−1
δf−1 1 δ 3 f− 1 δ 5 f− 1
2 2 2
f0 δ 2 f0 δ 4 f0 δ 6 f0
3 5
δf 1 δ f1 δ f1
2 2 2
f1 δ 2 f1 δ 4 f1
δf1 1 δ 3 f1 1
2 2
f2 δ 2 f2
δf2 21
f3

Example 6.4.1

Given f (x) = x4 , form a table of finite differences, for x from x = −2 in steps of one to 5.

Solution

Table 6.4: Showing a finite difference table for f (x) = x4

xr f r 4fr 42 f r 43 fr 44 fr 45 fr
−2 16
-15
−1 1 14
−1 −12
0 0 2 24
1 12 0
1 1 14 24
15 36 0
2 16 50 24
65 60 0
3 81 110 24
175 84
4 256 194
369
5 625

D.W-Ddumba numerical analysis i Page 123 of 133

 6.4. FINITE DIFFERENCE TABLES

We note that the fourth differences are constants and fifth differences are all zeros. We state
the following theorem to support this observation.

Theorem 6.4.1
Let Pn (x) be a polynomial of degree n(n ≥ 1 integer). Then

(i) 4r Pn (x) is a polynomial of degree n − r(r ≥ n integer)

(ii) 4n Pn (x) = constant

(iii) 4n+1 Pn (x) = 0.

Thus, parts (ii) and (iii) of the theorem confirm the validity of the results of the difference table
(4.4) in which fourth order differences are all constants equal to 24 and fifth differences are all
zero.
In fact the theorem can help us check the behaviour of a function (not polynomial) on an
interval. If the nth order differences are fairly constant on an interval for a particular function,
then such a function more or less behaves like a polynomial of degree n in that interval.

Exercise 6.6 For a polynomial

P6 (x) = x6 − 6x2 + 1,

Construct a table of forward differences on the interval [4, −4] with step size 1.

Exercise 6.7 Construct a table of backward differences for the polynomial in Exer6.6.

Exercise 6.8 Construct a table of differences for the function f (x) = ex on the interval [−2, 2]
with step size 0.5.
From your results what polynomial function fairly approximate ex on this interval?

Exercise 6.9 What do you think are the uses of finite difference tables

Page 124 of 133 D.W-Ddumba numerical analysis i

Chapter 7

Sample Questions

Question 1

(a) What is Polynomial interpolation?

(b) The Lagrange’s interpolating polynomial P (x) is given by,

n
X
P (x) = Li (x)fi
i=0

where the coefficients,

n
Y (x − xj )
Li (x) = i = 0, 1, 2, · · · , n
j=0,j6=i
(xi − xj )

(i) Prove that P (x) is unique.

(ii) Given the data
x 0 2 3
f (x) -2 4 10
Find a Lagrange quadratic interpolating polynomial that fits the data and use it to
interpolate f (1.5).

125
 Question 2

(a) state one advantage and one disadvantage of the following iterative methods used in the
solution of a non-linear equation f (x) = 0 :

(i) Bisection
(ii) Regular falsi
(iii) Newton Raphson’s

(b) (i) Prove that the general Newton-Raphson iterative method for finding the rth root of
a number, A > 0 is given by,

1 A
xn+1 = {(p − 1)xn + r−1 }
p xn

(ii) From 2b(i) deduce the corresponding iteration formula when p = −2, and use it to
find an approximate value of √12 correct to two decimal places.
√
(iii) Use 2b(ii) to deduce the approximate√ value of 2. Give reasons why this approach
to finding the approximate value of 2 is better than using the iterative formula
 
1 A
xn+1 = xn +
2 xn

Question 3

(a) The derivative of a function f (x) at a point x0 is defined by the limiting process,

f (x0 + h) − f (x0 )
f 0 (x0 ) = lim ,
h→0 h
provided the limit exists. However, a numerical approximation to f 0 (x0 ) involves differ-
ences such as the forward difference approximation

f (x0 + h) − f (x0 )
f 0 (x0 ) ≈ (3.1)
h
2
(i) Using (3.1), approximate f 0 (1) for f (x) = ex using a step size h = 0.1. Compare
your answer with the exact value of the derivative. Account for the difference in the
results. Suggest one way of improving on your result so obtained.
(ii) Geometrically represent the difference approximation (3.1).
(iii) Derive the difference approximation (3.1) together with its truncation error term.
State the order of the approximation

(b) For the approximation (3.1), find expression for

(i) total error and the optimal size of h

(ii) State two outstanding problems associated with numerical differentiation of a func-
tion f (x).

Question 4

Page 126 of 133 D.W-Ddumba numerical analysis i

 (a) Using the composite form of the Trapezoidal Rule with h = 0.25, approximate the value
of I where ˆ 3
I= ex sin x dx.
1
Compare your result with the exact value of I.
(b) Determine the step-length h needed so that the approximate value of I obtained from the
composite Trapezoidal rule in 4(a) has maximum truncation error 10−4 .

Question 5

(a) (i) State the general Runge-Kutta four stage method of fourth order popular for solving
initial value problems of the form, y 0 = f (x, y), y(xn ) = yn .
(ii) State the order of the local truncation error for the scheme in part (i).
(b) Given the initial value problem y 0 = x − y, y(0) = 0 with classical solution given by, y =
e−x + x − 1. Using h = 0.1, with Runge-Kutta fourth order method, find
(i) y(0.1)
(ii) y(0.2)
(ii) y(0.3)
(iii) y(0.4)
Make a table of comparison of the classical solutions and numerical solutions with the absolute
errors committed at each of the points 0.1, 0.2, 0.3 and 0.4. Account for the errors in the
numerical solutions obtained.

Question 6

(a) One of the oldest numerical techniques for solving a differential equation
dy
= f (x, y) (7.1)
dx
with the initial condition y0 = y(x0 ) where x0 and y0 are given constants and x0 not a singular
point of the function, is by Taylor series. Develop the Taylor series method for equation (7.1).
(b) Using the Taylor series method find an expression for the solution y(x) given that,
dy
= x3 − y
dx
and the initial condition y = 1 when x = 0. Use this expression to find values of y for
(i) x = x0 + h
(ii) x = x0 + 2h with h = 0.1

Question 7

(a) Runge-Kutta methods are popular Numerical techniques for solving ordinary differential
equations. Given the equation,
dy
= f (x, y) with y(a) = c,
dx

D.W-Ddumba numerical analysis i Page 127 of 133

 (i) State the general Runge-Kutta method of order four for the equation
(ii) Hence using the Runge-Kutta fourth process, find y(1) for the equation

y 0 = x + y, y(0) = 1, with h = 1

(iii) Compare the solution y(1) with the exact/analytic solution at that point.
(iv) Suggest anyway of increasing the accuracy in the numerical solution y(1).
(v) Explain the origin of the error in the numerical solution y(1).

(b) The oldest technique for obtaining a numerical solution of an ordinary differential equation
is the Euler’s method. Derive the Euler’s method together with its truncation error for
the equation,
dy
= f (x, y) with y(a) = c
dx
Question 8

(a) Given
y 0 = f (x, y) with y(x0 ) = y0 , (7.2)
derive the Euler’s method for solving equation (2)

(b)(i) For the equation

y 0 = −y with y(0) = 1
with step length h = 0.01, find y(0.04) using Euler’s method.
Compare your results with the analytic values at the discrete points. Account for the
difference in your results.

(ii) State one advantage and one disadvantage of Euler’s method as compared to Runge-
Kutta.

(c) Given a two point Gauss-Quadrature rule,

ˆ 1
f (x) dx ' f (α) + f (−α)
−1

(i) determine the value of α

(ii) Hence or otherwise compute
ˆ 4
sin t
dt.
0 t

Question 9

(a) State the defining equations for the following finite difference operators on a function f (x).

(i) µ (ii) δ (iii) E k (iv) 5

(b) Prove that E = ehD ; where D is a differential operator and E is a shift operator acting
on a function y(x).

(c) Let (x0 , y0 ), (x1 , y1 ) · · · , (xn , yn ) be a given set of (n + 1) points. Define the divided
differences

Page 128 of 133 D.W-Ddumba numerical analysis i

 (i) [x0 , x1 ]
(ii) [x0 , x1 , x2 ]
(iii) [x0 , x1 , x2 , x3 ]
(iv) [x0 , x1 , x2 , · · · , xn ]

(d) Using Newton’s divided difference formula, fit a cubic polynomial y(x) on the data :

x 0 1 2 3
y(x) 1 0 1 10

Hence approximate y(0.5).

Question 10

(a) (i) Define the finite difference operator 4 on a function y(x)

(ii) Newton’s Forward Difference Interpretation formula (NFDIF) is popular in inter-
polation of equispaced data near the beginning of tabular values. Given the set of
(n + 1) values, (x0 , y0 ), (x1 , y1 ) · · · , (xn , yn ) of x and y, Show that.

p(p − 1) 2 p(p − 1)(p − 2) · · · (p − n + 1) n

yn (x) = y0 +p4y0 + 4 y0 +· · ·+ 4 y0 ; x = x0 +ph
2! n!

(b) Using NFDIF, show that

(i)
42 y0 43 y0 44 y0
   
dy 1
= 4y0 − + − + ···
dx x=x0 h 2 3 4
(ii)
d2 y
   
1 2 3 11 4
= 2 4 y0 − 4 y0 + 4 y0 − · · · .
dx2 x=x0 h 12

(c) Given the data

x 0 2 4 6 8
y 7 13 43 145 367
 
d2 y
compute dx2
.
x=0

Question 11
´ xn
(a) (i) State the Trapezoidal rule for approximating the integral x0
y(x) dx on [x0 , xn ].
Write down the expression for the truncation error term.
(ii) Derive the rule in part 3(a)(i).
´ 1 dx
(b) Approximate I = 0 1+x 2 using Trapezoidal rule with stepsize h = 0.2. Hence approxi-

mate the value of π . Compare the numerical value of I with it’s analytic value. Explain
how the accuracy of I is affected by the value of h.

D.W-Ddumba numerical analysis i Page 129 of 133

 Question 12

(a) (i) With a relevant sketch, describe the Bisection method for approximating the root α
of a non-linear equation f (x) = 0 on [x0 , xn ].
(ii) Use the bisection method to find an approximation to the root of x3 − x − 1 = 0 on
the interval [1, 2] correct to one decimal place.

(b) Newton’s Raphson’s method is one of the popular schemes for solving a non-linear equa-
tion f (x) = 0. Prove that the general Newton Raphson’s iterative method for finding the
pth root of a positive number B is given by
 
1 B
xn+1 = (p − 1)xn + p−1 .
p xn
Using this scheme, approximate √1 correct to 4 decimal places. Take initial guess as 0.6.
2

Q13.(a) Common iterative schemes for solving non linear equations f (x) = 0 are of the form,
xr+1 = g(xr ) for r = 0, 1, . . .

(i) State what is meant by the iterative scheme being convergent.

(ii) State one advantage and one disadvantage of Newton’s method.

(b)(i) Prove that the general Newton Raphson’s method for finding the rth root of a number
(N > 0) is given by
1 N
xn+1 = {(r − 1)xn + r−1 }
r xn
(ii) From b(i), deduce the corresponding iterative formula when r = −2, and use it to find an
approximate value of √12 correct to 3 decimal places. Use initial approximation x0 = 0.7.

(iii) State two limitations of the Newton’s method.

Q14(a)(i) Explain what is meant by polynomial interpolation.

(ii) A function f (x) passes through the points (0, 1.0000) and (1, 2.7183). Construct La-
grange’s interpolating polynomial satisfying the above data. Use the polynomial to com-
pute the approximate value of f (0.5).

(b)(i) Prove that the bound on the rounding error when using Lagrange’s polynomial Pn (x)
does not blow beyond Σnk=0 |ξk ||Lk (x)| where ξ is the rounding error in fk and Lk (x) are
the Lagrange coefficients.

(ii) Given that f (x) = ex for part a(ii) and data were rounded to four digits, show that the
effect of rounding errors in the on P1 (x) maintains the same maximum magnitude of
1
2
× 10−4

(c) Give the defining equations for the finite difference operators, µ, σ and E. Hence prove
1
the finite difference identity µσ + 12 σ 2 = E 2 σ.
´x
Q15(a) x0 f (x)dx ' h3 [f0 + 4f1 + f2 ] is the popular Simpson’s rule for approximating
an integral.

(i) State two major sources of error when the Simpson’s rule is used to approximate a definite
integral.

Page 130 of 133 D.W-Ddumba numerical analysis i

 (ii) By interpolating f (x) by a Lagrange polynomial P2 (x) of degree two,
[Link] (x) = P2 (x) + E(x), derive the Simpson’s rule.
´1
(iii) Using the simpson’s rule, approximate 0 x2 dx and comment on your result.
´2 2
(b) It is required to obtain 0 ex dx exact to 4 decimal places.
What should h be for Simpson’s rule?

Q16(a) The second order central difference approximation for f 0 (x) is given as

f (x + h) − f (x − h)
f 0 (x) ' .
2h

(i) Give the geometrical interpretation of the central difference approximation.

(ii) By the help of appropriate Taylor series, derive the central difference approximation with
its truncation error.

(iii) Using the central difference approximation, approximate f 0 (2) with h = 1.0 and h = 0.1
for f (x) = x2 . What do you notice from your results?

(b) Given that |e0 | and |e1 | are the maximum absolute errors committed when evaluating
f (x + h) and f (x − h) respectively, and |e0 |, |e1 | < e, show that for the central difference
approximation,
e M h2
(i) |Total error| ≤ h
+ 6
1
) where M =max|f 000 (x)| for a − h ≤ c ≤ a + h.
3e 3
(ii) Optimal h = ( M

D.W-Ddumba numerical analysis i Page 131 of 133

Page 132 of 133 D.W-Ddumba numerical analysis i
Chapter 8

Further Reading

1. Curtis F. G and Patrick O.W. (1994). Applied Numerical Analysis. Addison - Wesley
Publishing Company.

2. Maron M. J. (1982). Numerical Analysis. Macmillan Publishing Company.

3. Maron M.J. and Robert J.L. (1991). Numerical Analysis. Wadsworth Publishing Com-
pany.

4. Samuel D. Coute (1980). Elementary Numerical Analysis. Mcgraw hill Kogakusha, Ltd.

5. Masenge R.W.P. (1987). Fundamentals of Numerical Methods. University of Dar-es-

Salaam. Tanzania.

133