UNIVERSITY OF MINES AND TECHNOLOGY

FACULTY OF ENGINEERING

DEPARTMENT OF MATHEMATICAL SCIENCES

LECTURE NOTES ON

NUMERICAL METHODS (PE/NG/RP 258)

Course Instructors

Dr M. V. CRANKSON and K. Agyarko

June, 2022

#WEDEY4U
 UMaT, Department of Mathematics

CHAPTER 1
Arithmetic and Error Analysis

1.0 Introduction
Definition — Numerical
Method
– A numerical method is a complete and unambiguous set of procedures for the solution of a

problem, together with computable error estimates.

– any method that uses only four basic arithmetic operations: + , – , ÷ , ×

1.1 Nature of numerical problems

Solving mathematical equations is an important requirement for various branches of science. The field
of numerical analysis explores the techniques that give approximate solutions to such problems with
the desired accuracy.

1.2 Errors, basic definitions

Whenever numbers with a finite quantity of digits are used in computations and whenever real
measurements are performed, one has to acknowledge that errors will be made. This is not grave by
it-self. What matters is having an estimate of their size and knowing that, as calculations are made,
they can propagate. In the end, the best one can do is to bound the absolute error, that is, to know a
reasonable value (the bound) which is larger than the error, in order to assess, with certainty, how far
the real value can be from the computed one.

1.3 Sources of error

Error appears in different ways. On one hand, any measurement is subject to it (this is why any
measuring device is sold with an estimated margin of precision); this is intrinsic to Nature and one can
only take it into account and try to assess its magnitude (give a bound). On the other hand,
computations performed in finite precision arithmetic both propagate these errors and give rise to new
ones precisely because the quantity of available digits is finite.

Numerical Methods Page 1

#WEDEY4U
 UMaT, Department of Mathematics

The following are some of the sources of error:

• Measurement error, already mentioned. This is unavoidable.
• Truncation error: happens whenever a number (datum or the result of a computation) has
more digits than available and some of them must be “forgotten”.
• Rounding error: takes place when a number is rounded to a specified precision.
• Loss of significance (also called cancellation error): this appears when an operation on two
numbers increases relative error substantially more than the absolute error. For example,
when numbers of very similar magnitude are subtracted. The digits lost to truncation or
rounding are too relevant and the relative error is huge. A classical example is the instability
of the solution of the quadratic equation.
• Accumulation error: which appears when accumulating (with additions, essentially) small
errors of the same sign a lot of times.

All the above errors may take place when working with finite arithmetic. Numerically computed
solutions are subject to certain errors. Mainly there are three types of errors. They are inherent errors, truncation
errors and errors due to rounding.

1. Inherent errors or experimental errors arise due to the assumptions made in the mathematical
modeling of problem. It can also arise when the data is obtained from certain physical
measurements of the parameters of the problem. i.e., errors arising from measurements.
2. Truncation errors are those errors corresponding to the fact that a finite (or infinite) sequence
of computational steps necessary to produce an exact result is “truncated” prematurely after a
certain number of steps.
3. Round of errors are errors arising from the process of rounding off during computation.
These are also called chopping, i.e. discarding all decimals from some decimals on.

1.4 Error in Numerical Computation

Due to errors that we have just discussed, it can be seen that our numerical result is an approximate
value of the (sometimes unknown) exact result, except for the rare case where the exact answer is
sufficiently simple rational number.

Numerical Methods Page 2

#WEDEY4U
 UMaT, Department of Mathematics

If is an approximate value of a quantity whose exact value is a, then the difference is called
the absolute error of or, briefly, the error of . Hence, , i.e.
Approximate value = True value + Error

For example, if is an approximation to a 10.5, then the error is 0.02. The relative error, r , of is defined
by

Error

r
a True value

For example, consider the value of 2( 1.414213...) up to four decimal places, then
2 1.4142 Error .

Error 1.4142 1.41421 0.00001,

Taking 1.4.1421 as true or exact value. Hence, the relative error is

We note that
.

if .

We may also introduce the quantity and call it the correction, thus, , i.e.
is much less than

True value = Approximate value + Correction

Error bound for is a number such that i.e., .

1.5 Number representations

1.5.1 Floating-point representation

Most digital computers have two ways of representing numbers, called fixed point and floating point.
In a fixed point system the numbers are represented by a fixed number of decimal places e.g.

Numerical Methods Page 3

#WEDEY4U
 UMaT, Department of Mathematics

62.358, 0.013, 1.000.

In a floating point system the numbers are represented with a fixed number of significant digits, for
example also written as 0.6238E03 , or more simply
0.6238 03, 0.1714-13 , .

Definition – floating

The representation of a number f as

where
the sign, s, is either +1 or – 1,
the mantissa, m, satisfies ,
the exponent, x, is an integer

is called the floating point representation of f in base M.

Example 1

The decimal number f 8.75 , written in its floating-point representation in base 10 is ,

and has s 1, m 0.875 and x 1. Its binary equivalent is

where .

Definition — Significant
digit

Numerical Methods Page 4

#WEDEY4U
 UMaT, Department of Mathematics

Significant digit of a number c is any given digit of c, except possibly for zeros to the left of the first

nonzero digit that serve only to fix the position of the decimal point. (Thus, any other zero is a

significant digit of c). For example, each of the number 1360, 1.360, 0.01360 has 4 significant

digits.

Number of digits starting from the first non-zero on the left side
Example 2

Number of Number of
significant digits significant digits
2345000 7 5 1
2.345000 7 5.0 2
0.023450 5 5.000 4
0.02345 4

1.5.2 Roundoff (Machine) Error

In a computer the mantissa m of a number x is represented using a fixed number of digits. This means
that there is a limit to the precision of numbers represented in a computer. The computer either
truncates or rounds the mantissa after operations that produce extra mantissa digits. The absolute
error in storing a number on computer in base 10 using d digits can be calculated as follows. Let mˆ
be the truncated mantissa, then the magnitude of the absolute error is

Roundoff error can be particularly problematic when a number is subtracted from an almost equal
number (subtractive cancellation). Careful algebraic rearrangement of the expressions being evaluated
may help in some cases.

Numerical Methods Page 5

#WEDEY4U
 UMaT, Department of Mathematics

Example 3

The expression

However, when evaluated on my calculator the answer given is 590, 000!! Clearly the calculator's
answer is catastrophically incorrect.

Round off rule to discard the (k + 1)th and all subsequent decimals

(a) Rounding down If the number at (k + 1)th decimal to be discarded is less than half a unit in the
k th place, leave the k th decimal unchanged. For example, rounding of 8.43 to 1 decimal gives
8.4 and rounding of 6.281 to 2 decimal places gives 6.28.

(b) Rounding up If the number at (k + 1)th decimal to be discarded is greater than half a unit in
the k th place, add 1 to the k th decimal. For example, rounding of 8.48 to 1 decimal gives 8.5
and rounding of 6.277 to 2 decimals gives 6.28.

(c) If it is exactly half a unit, round off to the nearest even decimal. For example, rounding off
8.45 and 8.55 to 1 decimal gives 8.4 and 8.6 respectively. Rounding off 6.265 and 6.275 to 2
decimals gives 6.27 and 6.28 respectively.

1.6 Truncation Error

This type of error occurs when a computer is unable to evaluate explicitly a given quantity and instead
uses an approximation.
For example, an approximation to sin x might be calculated using the first three terms of the Taylor

series expansion,

where

Numerical Methods Page 6

#WEDEY4U
 UMaT, Department of Mathematics

If x is in the interval [ 0.2, 0.2] then the truncation error for this operation is at most

0.27 9

3.10
7!

EXERCISES (Try them)

1. Define the following terms

(a) Truncated Error (b) Round-off Error (c) Inherent Error (d) Relative Error
(e) Absolute Error (f) Floating Point
2. Determine the absolute and relative errors involved if x = 2/3 is represented in normalised
decimal form with 6 digits by
(a) round-off (b) truncation.
3. Given that 5 digit chopping is used for arithmetic calculations involving x and y where x = 1/3
and y = 5/7. Determine the absolute and relative errors involved.
4. If x = 3.536, determine the absolute error and relative error when (a) x is rounded (b) x is
truncated to two decimal digits.
5. If the number x = 57.46235 is rounded off to four significant figures, find the absolute error,
22
relative error and the percentage relative error. 6. If the approximate value of = 7

is 3.14, determine the absolute error, relative error

and relative percentage error.

7. Determine the true error and true percentage relative error for each case.

Numerical Methods Page 7

#WEDEY4U
 UMaT, Department of Mathematics

(a) If the measured length of a track is approximated by 9999 cm and the true value is 10,000
cm
(b) If the measured length of a track is approximated by 9 cm and the true value is 10 cm.
8. Find the roots of the following equations using 4 significant figures in the calculation.

9. Convert the decimal number (which is in the base 10) 81.5 to its binary form (of base 2).
10. Convert the binary number 1010.101 to its decimal form.

CHAPTER 2 Linear Equations

2.0 Introduction
In two dimensions a line in a rectangular xy-coordinate system can be represented by an equation of
the form
ax by c a b+ = ( , not both 0)

and in three dimensions a plane in a rectangular xyz-coordinate system can be represented by an

equation of the form
ax by cz d a b c+ + = ( , , not all 0)

These are examples of “linear equations”, the first being a linear equation in the variables x and y and
the second a linear equation in the variables x, y, and z. More generally, we define a linear equation in
the n variables x x1 2, , ,x to be one thatn can be expressed in the form

ax a x1 1 + 2 2 + +a xn n =b (2.1)

Where a a1, 2,...,an and b are constants and the a s' are not all zero. In the special case where b=0,
Equation (2.1) has the form

Numerical Methods Page 8

#WEDEY4U
 UMaT, Department of Mathematics

ax a x1 1 + 2 2 + +a xn n = 0

which is called a homogeneous linear equation.

Example 1
Observe that a linear equation does not involve any products or roots of variables. All variables occur
only to the first power and do not appear, for example, as arguments of trigonometric, logarithmic, or
exponential functions. The following are linear equations:
1. x+ =3 7y

2. x y− + =−3z 1

3. x1 − + + =2x2 3x3 x4 0

The following are not linear equations

1. x+3y2 =4

2. sin x y+ = 0

3. 3 2x+ − =y xy 5

4. x1 +2x2 + =x3 1

2.1 System of Linear Equations

A finite set of linear equations is called a system of linear equations or, more briefly, a linear system.
The variables are called unknowns. A general linear system of m equations in the n unknowns x x1 2, ,
, xn can be written as

122 + ax1n n = b1 a x11 1 + a x +

211 222 + ax2 n n = b2
ax+ax+
axn11 + axm 22 + + axnn n = bn

Numerical Methods Page 9

#WEDEY4U
 UMaT, Department of Mathematics

A solution of a linear system in n unknowns x x1 2, , ,x is na sequence of n numbers s s1 2, , ,s for which

n

the substitution x1 =s1, x2 =s2, ..., xn =sn

makes each equation a true statement. More generally, a solution x1 =s1, x2 =s2, ..., xn =sn of a linear

system in n unknowns can be written as

(s s1 2, ,...,sn)

which is called an ordered n-tuple. With this notation it is understood that all variables appear in the
same order in each equation. If n = 2, then the n-tuple is called an ordered pair, and if n = 3, then it
is called an ordered triple.
In general, we say that a linear system is consistent if it has at least one solution and inconsistent if it
has no solutions. Thus, a consistent linear system of two equations in two unknowns has either one
solution or infinitely many solutions—there are no other possibilities.

2.2 Direct Methods for Solving Systems of Linear Algebraic Equations

In direct methods, a finite number of arithmetic operations lead to an “exact” (within round-off errors)
solution. Examples of such direct methods include Gauss elimination, Gauss-Jordan elimination, the
matrix inverse method, and LU factorization. The average number of operations to solve a system of
linear equations for these methods is n3. In many applications where linear

Numerical Methods Page 10

#WEDEY4U
 UMaT, Department of Mathematics

systems appear, one needs to solve for many different vectors b. The systems of equations are
used to represent physical problems that involve the interaction of various properties. The variables in
the system represent the properties being studied, and the equations describe the interaction between
the variables.

A system of linear equations has the following form

ax
111 + ax
122 + + ax
1n n = b1
ax211 ax
+ 222 + ax
+ 2n n = b2
(2.2 )
axn11 + axm 22 + + axnn n = bn
where are unknown variables, are the coefficients, and

are the nonhomogeneous terms. The first subscript i identifies the row of the
equations and the second subscript j identifies the column of the system of equations. The system

(2.2) can also be written as a matrix equation

2.3
or in a compact form as

where A is a square matrix of order n and b is a vector in . We always assume unless explicitly
stated that A is non-singular, so that the system is consistent and has a unique solution.

Methods for solving linear systems are normally taught in mathematical classes and courses of linear
algebra. The standard syllabus includes the substitution method, Cramer’s rule, and the inverse matrix.

Numerical Methods Page 11

#WEDEY4U
 UMaT, Department of Mathematics

2.1 Crammer’s Rule

Cramer’s rule provides a method for solving a system of linear algebraic equations for which the
associated matrix problem has a coefficient matrix, which is nonsingular. It is of no use if this
criterion is not met. So, let us assume the n-by-n matrix A is nonsingular, that b is a known vector in ,
and that we wish to solve the equation for an unknown (unique) vector

. Cramer’s rule requires the construction of matrices , where each is built from the
original A and b. These are constructed as follows: the jth column of A is replaced by b to form Aj.

Example 1

Construction of A1, A2, A3, where A is 3-by-3

Example 2

Use Cramer’s rule to solve the system of equations

Numerical Methods Page 12

#WEDEY4U
 UMaT, Department of Mathematics

Solution:
Here, A and b are given by

Thus,

Numerical Methods Page 13

#WEDEY4U
 UMaT, Department of Mathematics

2.2 Inverse Matrix

The inverse of a matrix is analogous to the reciprocal of a scalar. We already encountered

the inverses of matrices corresponding to elementary row operations. We begin with the formal
definition.
Definition — Inverse Matrix
Let I be the n-by-n identity matrix. Let A be an n-by-n matrix. A n-by-n matrix B, is called an
inverse of A if and only if AB = I and BA = I. B is often written as A-1.

The requirements for obtaining a unique inverse of a matrix are:

1. The matrix is a square matrix
2. The determinant of the matrix is not zero (the matrix is non-singular) The inverse of
a matrix is also defined by the relationship:
A A−1 = I
Step

1. Augment the n-by n matrix A with the n-by-n identity matrix In : A I| n

2. Pivot on matrix elements in positions i = j (i,j=1, 2, 3, … n) continuing through n,n in that
order, with the goal of creating a copy of the identity matrix (In) in the left portion of the
augmented matrix. If one of the pivoting elements is zero, then first interchange it’s row with
a lower row. If no such interchange produces a non-zero pivot element, then the matrix A has

no inverse. A I| n In | A−1

3. When step [2] above is done, the right half of the latest augmented matrix will be the desired

inverse, A−1 , write it separately, and you're done, as in the example below.

Example 1

Use inverse matrix to solve the system of equations

Numerical Methods Page 14

#WEDEY4U
 UMaT, Department of Mathematics

Solution:

2.3 Gaussian Elimination

Consider the following system of linear simultaneous equations:

Gauss elimination is a popular technique for solving simultaneous linear algebraic equations. It
reduces the coefficient matrix into an upper triangular matrix through a sequence of operations carried
out on the matrix. The vector b is also modified in the process. The solution vector (x) is obtained
from a backward substitution procedure.
Suppose we subject this system to the system of following operations:
1. Multiplication of one equation by a non-zero constant.
2. Addition of a multiple of one equation to another equation.
3. Interchange of two equations.
If the sequence of operations produces the new system Ax b = , then both the systems Ax b= and
Ax b = are equivalent. In particular, then A is invertible if A is invertible. In Gaussian elimination
method, we adopt this and the elimination process is based on this theorem.

Numerical Methods Page 15

#WEDEY4U
 UMaT, Department of Mathematics

In Gaussian elimination method, the unknowns are eliminated such that the elimination process leads
to an upper triangular system and the unknowns are obtained by back substitution. It is assumed a11
0 . The method can be described by the following steps:
Step 1: Eliminate x1 from the second and third equations.

Step 2: Eliminate x2 from the new third equation using the new second equation by assuming the new

a22 0.

Step 3: To find x1, x2 and x3, we apply back substitution. Let

denote any n-by-n matrix, the elements are called diagonal elements of A. if the

matrix A is of the form

we say that AU is upper triangular matrix and if A is of the form

we say AL is lower triangular matrix.

If the system Ax b= is such that A is upper triangular, then the system takes the form

Numerical Methods Page 16

#WEDEY4U
 UMaT, Department of Mathematics

In this case if aii 0 for any i, then the system may be solved easily as follows. First, solve the last

equation for xn ; this gives

bn
xn =
.
ann

Next solve the last but one equation for xn−1 ; this gives

n = bn-1 - an-1,nxn
x -1
an-1,n-1

Proceeding in this manner we may solve the ith equation for by

2.3.1 The Algorithm for the Gaussian Elimination

Let consider the matrix and decide to reduce A to an equivalent upper

triangular matrix. This means we want to manipulate all sub-diagonal elements to zero. To eliminate

a21, find a factor f which when multiplied onto the first row and added to the second, term by term,

make a21 zero. The second row becauses:

Numerical Methods Page 17

#WEDEY4U
 UMaT, Department of Mathematics

To develop an algorithm for the elimination (1) must be generalized. The first stage is an

expression for f . From the first equation in (1), since a21 = 0,

Notice the pattern of the subscript in (1); the column subscripts only vary from equation to equation.
These equations can be summarized by:

 a2¢k= a2k + fa1k , k =1, 2, 3,, n (3)

At each stage of the elimination, the diagonal element and the row it is in play a central role. These
elements are called pivot elements and their rows, pivot rows. If the subscript i refers to the pivot row
and j to the elimination row (the row being operated on), then (2) and (3) become:
a
f =- ji

aii

a¢=a + faik
jk jk

Example 1

Solve the system

to 3 decimal places, by
(a) reducing to an equivalent upper triangular system.
(b) reducing to an equivalent lower triangular system.

Numerical Methods Page 18

#WEDEY4U
 UMaT, Department of Mathematics

Solution:
(a) The augmented matrix of the system is

f1 =- 3.02.1 a21¢ =3.0+æçè- 3.02.1ö÷ø ´ 2.1=0; a22¢ =2.5+æçè-

æ
3.02.1ö÷ø ´ -4.5=8.92857; è 2.1÷ø ´ -2.0 =7.15714; a24¢ =3.21+ çè-

3.0 ö æ 3.0ö
2.1 ÷ø ´19.07 =-24.03286; a23¢ =4.3+ ç-

f2 =- æçè -2.16.0ö÷ø

a31¢ =-6+æçè 6.02.1ö÷ø ´ 2.1=0; a32¢ =3.5+æçè 6.02.1ö÷ø ´ -

æ 6.0ö
4.5=-9.35714; a33¢ =2.5+ ç

æ 6.0 ö
è 2.1÷ø ´ -2.0 =-3.21429; a34¢ =- 18.25+çè 2.1 ÷ø ´19.07
=36.23571;

The above reduction produces the equivalent matrix

f3 =- æçè -8.926579.35714ö÷ø =1.04800

Numerical Methods Page 19

#WEDEY4U
 UMaT, Department of Mathematics

a32¢ =-9.35714 +(1.04800) ´(8.92857) =0

a33¢ =-3.21429 +(1.04800) ´(7.15714) =4.28639

a34¢ =36.23571+(1.04800) ´(-24.03286)

=11.04927

Solving by back substitution, we have

x3 = = 2.57776

From the second row: 8.92857x2 +7.15714x3 =−24.03286

x2 =-=-4.75801

From the first row: 2.1x1 −4.5x2 −2.0x3 =19.07 x1 =1.34023 Therefore to 3

decimal places we have:

x1 =1.340; x2 =−4.758; x3 = 2.578

(b)

f1 =-æçè -2.52.0ö÷ø =0.8

a14 =19.07 +0.8(-18.25) =4.47; a13 =-2.0+0.8(2.5) =0; a12

=-4.5+0.8(3.5) =-1.7; a11 =2.1+0.8(-6.0) =-2.7;

Numerical Methods Page 20

#WEDEY4U
 UMaT, Department of Mathematics

f2 =-æçè 4.52.5ö÷ø =-1.72

a24 =3.21- 1.72(-18.25) =34.6; a23 =4.3- 1.72(2.5) =0; a21

a22 =2.5- 1.72(3.5) =-3.52; =3.0- 1.72(-6.0) =13.32;

The new matrix is

f3 =- æçè 3.521.7 ö÷ø =- 0.48295

a14 =4.47 - 0.48295(34.6) =- 12.24007; a21 =- 17 - 0.48295(- 3.52) =0

a11 =- 2.7 - 0.48295(13.32) =- 9.13289

The equivalent lower triangular matrix is

Solving by forward substitution we have

x1 ==1.34022

From the second row we obtain: 13.32x1 - 3.52x2 =34.6 Þ x2 =-4.75803

Finally, from the third row we have: -6x1 +3.5x2 +2.5x3 =-18.25 Þ x3 =2.57777 Therefore to 3
decimal places we have:

Numerical Methods Page 21

#WEDEY4U
 UMaT, Department of Mathematics

x1 =1.340; x2 =-4.758; x3 =2.578

2.3.2 Partial and Complete Pivoting

If at any stage during the reduction of a system to an equivalent upper or lower triangular system it

aji is found that the pivot aii is too

small to quotients of the form are inaccurate, we may examine aii

the coefficients of aji for j >i for the largest number. If the coefficient is aki , interchange the ith equation

with the kth equation below so thataki becomes the new pivot and the reduction continues.

This modified procedure aimed at avoiding small is called Gaussian elimination with partial
pivoting.
If at the ith stage of the reduction the pivot aii turns out to be too small, we may examine all the

coefficients in the system and pick the largest one. If this is apq then interchange the ith row with the

pth row and the ith column with the qth column. This places the largest coefficient available at this stage

in the position (i, i). The reduction is then continued. This modified process is called Gaussian

elimination with complete pivoting.

The extra work involved in Gaussian elimination with complete pivoting does not provide solutions
any more accurate than the solution obtained by using Gaussian elimination with partial pivoting.
Hence, elimination with partial pivoting is always preferred.
Example 2

Solve the system

Solution:
The augmented matrix of the system is

Numerical Methods Page 22

#WEDEY4U
 UMaT, Department of Mathematics

a11 = 0.01 is to be the first pivot element, but it is too small compared with a21 = 2.1and a31 = 3.2.
We therefore replace 0.01 with 3.2. This requires that we interchange row 1 and row 3. This
gives

f1 =- æçè 3.22.1ö÷ø =- 0.656250

a22 =- 0.02 - 0.656250(- 0.01) =- 0.013438; a23 =0.2 - 0.656250(4.2) =- 2.556250

a24 =1.1- 0.656250(2.1) =- 0.278125

f2 =- æçè 0.013.2 ö÷ø =- 0.003125

a32 =2.0 - 0.003125(- 0.01) =2.000031; a33 =3.0 - 0.003125(4.2) =2.986875

a34 =0.1- 0.003125(2.1) =0.093438

The equivalent matrix is

The next pivot element should be a22 =-0.013438, which is small compared with a32 = 2.000031.
We therefore interchange row 2 and row 3 to make a32 = 2.000031 the pivot element.
This gives

Numerical Methods Page 23

#WEDEY4U
 UMaT, Department of Mathematics

æ- 0.013438 ö = 0.06719
f3 = - ç
è 2.000031 ÷ø

a33 =- 2.55650 +0.06719(2.986875) =- 2.536431;

a34 =- 0.278125 +0.06719(0.093438) =- 0.277497

The equivalent lower triangular matrix is

Solving by back substitution, we have

- 2.536431x3 =-0277497 Þ x3 =0.109405

2.000031x2 +2.986875(0.109405) =0.093438 Þ x2 =-0.116669

3.2x1 - (0.01)(-0.116669)+4.2(0.109405) =2.1 Þ x1 =0.512291

Therefore to 4 decimal places we have: x1 =0.5123; x2 =-0.1167; x3 =0.1094

2.4 Gauss Jordan Method

The method is based on the idea of reducing the given system of equations Ax = b, to a diagonal
system of equations Ix = d, where I is the identity matrix, using elementary row operations. We know
that the solutions of both the systems are identical. This reduced system gives the solution vector x.
This reduction is equivalent to finding the solution as x A 1b.

In this case, a system of 3 equations in 3 unknowns

Numerical Methods Page 24

#WEDEY4U
 UMaT, Department of Mathematics

is written as

After some linear transformations, we obtain the 3-by-3 system as

To obtain the system as given in (2), first we augment the matrices given is (1) as,

and after some elementary operations, it is written as,

This helps us to write the given system as given in (2). Then it is easy to get the solution of the system
as x1 d1, x2 d2 and x3 d3.

Elimination procedure: The first step is same as in Gauss elimination method, which is, we make the
elements below the first pivot in the augmented matrix as zeros, using the elementary row
transformations. From the second step onwards, we make the elements below and above the pivots as
zeros using the elementary row transformations. Lastly, we divide each row by its pivot so that the
final matrix is of the form (3). Partial pivoting can also be used in the solution. We may also make the
pivots as 1 before performing the elimination.
Numerical Methods Page 25

#WEDEY4U
 UMaT, Department of Mathematics

Example 1

Solve the following system of equation

using the Gauss-Jordan method without partial pivoting

Solution:

We have the matrix form as

Then the augmented matrix is

To do the eliminations follow the operations, R2 = R2 – 4R1, and R3 = R3 – 3R1. This gives,

Then, R1 = R1 + R2 and R3 = R3 + 2R2, this gives

Again, R1 = R1 – (4/10)R3, and R2 = R2 – (5/10)R3 which gives,

Numerical Methods Page 26

#WEDEY4U
 UMaT, Department of Mathematics

Now, making the pivotal as 1, R2 = (-R2) and R3 = (R3/(-10)), we get

Hence,

Therefore, the solution of the system is,

2.5 LU decomposition
In linear algebra, LU decomposition is also called LU-factorization. The nonsingular matrix A has an
LU-factorization if it can be expressed as the product of a lower-triangular matrix L and an upper
triangular matrix U, e.g.

It turns out that this factorization (when it exists) is not unique. If L has ones on it’s diagonal, then it
is called a Doolittle factorization. If U has ones on its diagonal, then it is called a Crout
T
factorization. When U = L it is called
a Cholesky decomposition.

2.5.1 Using LU factorization to solve a system

Numerical Methods Page 27

#WEDEY4U
 UMaT, Department of Mathematics

Assuming Ax = b can be factorized into

LU x = b Then

let y = Ux and use forward substitution to solve

Ly = b

for y. Then use backward substitution to solve

Ux = y

for x. This has the advantage that once A = LU is known we only have to do a forward and back
substitution again for a different righthand side b.

2.5.2 Doolittle’s Method

Multiplying out

Comparing coefficients implies that

Sometimes it is necessary to use pivoting.

Numerical Methods Page 28

#WEDEY4U
 UMaT, Department of Mathematics

Example 1

Use LU factorization to solve

apply the Doolittle’s method.

Solution:

By the Doolittle’s method implies

Comparing coefficients we find

Now let y = Ux ⇒ LUx = b becomes Ly = b, which is solved by forward substitution (i.e. start at
top of matrix):

Numerical Methods Page 29

#WEDEY4U
 UMaT, Department of Mathematics

Then solve Ux = y for x by back substitution:

2.6 Thomas Algorithm for Solving Tridiagonal Systems

Consider the system of linear simultaneous algebraic equations given by
Ax b=
T T

where A is a tridiagonal matrix, x = x x1, 2, , x n

and b b b= , ,
1 2 , bn . Hence, we consider a
4 × 4 tridiagonal system of equations given by

a12 a13 0 0 x1 b1

a21 a22 a23 0 x2 =

a32 b2
a41 a33 x3 b3
0 a31

a42 x4 b4
0 0
The equation can be written as
a x12 1 + a x13 2 = b1
a x21 1 + a x22 2 + a x23 3 = b2

a x31 2 + a x32 3 + a x33 4 =

b3 a x41 3 + a x42 4
= b4
The system of equations is solved using Thomas Algorithm which is described in three steps as shown
below:

Step 1: Set y1 = a12 and compute

Numerical Methods Page 30

#WEDEY4U
 UMaT, Department of Mathematics

a ai1 (i−13) ,n
yi = ai2 − i = 2, 3,
yi−1

b1 and compute Step

2: Set z1 =
a12
bi − azi1 i −1
zi = i = 2,3, ,n
yi
axi 3 i +1
Step 3: xi =−
zi i =−
n 1, n − 2, ,1,where xn = zn
yi
Example 1

Solve the following system of equation by Thomas Algorithm

3x1 − x2 = 5

2x1 −3x2 + 2x3 = 5 x2 +

2x3 + 5x4 =10 x3 −

x4 =1

Solution:
Let

3 −1 0 0 x1 5 Here
2 −3 2 0
x2 =
1 2
0 0 1 5

0 5 x3 10

−1 x4
1
a2, a3, a4 = 2,1,1

b b b b1, 2, 3, 4 = 3, 3, 2,− −1

Numerical Methods Page 31

#WEDEY4U
 UMaT, Department of Mathematics

c c1, 2, c3 =− 1, 2, 5
Step 1: Set y1 =b1 and compute

,n
aci i−1 i = 2,3, yi =
i−1
−bi y
y1 = 3

a c2 1 3 2(−1) =− 7
i = 2, y2 = b2 − =− −
y1 3 3

a c3 2 2 1 2 = 20
i = 3, y3 = b3 − y2 = − − 7 7

a c4 3 1 1.5 =− 55
i = 4, y4 = b4 − =− −
y3 20 20
7

,n − a zi
d1 5, zi = di i−1
i = 2, 3,
Step 2: Set z1 = =
b1 3 yi
5

z = d2 − a z2 1 = 5− 2 3 =− 5

i = 2, 2 y2 −7 7

−−

i = 3, z3 = d3 − a z3 2 = 101 75 = 75
20
y3 7 20
75

Numerical Methods Page 32

#WEDEY4U
 UMaT, Department of Mathematics

i = 4, z4 = 4 −ya z4 3 = 1−155
20 =1 d
4 −
20
cxi i+1
; i = −n 1, n ,1; xn = z
− 2, n Step 3: Set xi = −zi
yi

x4 = z4 =1
c x3 4
75 5 1 2
i = 3, x3 = z3 − = − =
20 20
y3 7
cx 5 2
i = 2, x2 = z2 − 2 3y =− −7 72 =1
2 −
3
cx 5 (−1) 1
i =1, x
1 = z1 − 1 2y1 = −3 3 =2

EXERCISES (Try all)

1. Using the Gaussian elimination method, solve the system of equations [A] {x} = {b} where
1 1 1 1 3
2 −1 3 0 3
2 0 b =
A = 3
0 2 and
0
1
1
−1
0
2. Solve the following set of equations by Gauss-Jordan method.

2x1 + x2 − 3x3 =11

Numerical Methods Page 33
#WEDEY4U
 UMaT, Department of Mathematics

4x1 − 2x2 + 3x3 = 8

6
−2x1 + 2x2 − x3 = −

1 3 6

3. Express the matrix A = 2 8 16 as the product LU of a lower triangular matrix L

and
5 21 45

0 0 u13
1 u11
1
upper triangular matrix U , where L = l21 u23
0 and U= 0 u12
l32
l31 1 0 u22 u33
0
By writing Ux = y and Ly =b, solve the system of equation

x+ + =3y 6z 17
2x+ +8y 16z= 42
5x+ 21y+ 45z= 91
4. Solve each of the following systems by Gaussian elimination with partial pivoting.

4x1 − x2 − x3 −3x4 = 4
x
i. 1 − 4x2 + x3 + x4 = 0

x1 + x2 − 4x3 + x4 =−5
3x1 + x2 + x3 − 4x4 =−8

4.331x1 − 0.132x2 + 0.029x3 + 0.02x4 = 3.210 ii. 1

+ 4.397x2 − 0.332x3 − 0.115x4 =−1.001
−1.031x
0.415x1 −1.447x2 −5.137x3 + 2.014x4 = 7.394
1.974x1 − 2.106x2 + 0.847x3 − 7.130x4 =−5.241

Numerical Methods Page 34

#WEDEY4U
 UMaT, Department of Mathematics

2.66x1 +1.06x2 +1.09x3 = 0.60 iii.

1.06x1 + 2.66x2 +1.09x3 = 2.26

0.24x1 +1.24x2 + 2.78x3 = −1.13

5. Apply Gauss elimination method to solve the equations:

7. Solve the following system Gauss Jordan Method:

i.

ii.

Numerical Methods Page 35

#WEDEY4U
 UMaT, Department of Mathematics

CHAPTER 3
Iterative Methods for Solving Linear Systems
3.0 Introduction
The methods discussed in the previous section belong to the direct methods for solving systems of
linear equations; these are methods that yield solutions after an amount of computations that can be
specified in advance.

Numerical Methods Page 36

#WEDEY4U
 UMaT, Department of Mathematics

In this section, we discuss indirect or iterative methods in which we start from an initial value and
obtain better and better approximations from a computational cycle repeated as often as may be
necessary, for achieving a required accuracy, so that the amount of arithmetic depends upon the
accuracy required.

Consider a linear system of n linear equations in n

unknowns of the form

in which the diagonal elements aii do not vanish. Now the system (1) can be written as

Numerical Methods Page 37

#WEDEY4U
 UMaT, Department of Mathematics

3.1 Jacobi’s Iteration Method

Jacobi’s iteration method, also called the method of simultaneous displacements, is as follows:

,
Step 2: Similarly,  x1(2),x2(2) ,xn(2)are evaluated by just replacing xr(0) in the right hand sides

equations in (3) by xr(1).

Step n 1: In general, if x1(n),x2(n),,x(n) are a system of nth approximations, then the next
 n approximation is given by the

formula

Numerical Methods Page 38

#WEDEY4U
 UMaT, Department of Mathematics

The system in (4) can also be briefly described as follows:

æ ö

(k+1) = 1 ç - n (k ) ÷ ( =1,2, , =0,1,2, )

x aii çb åj=1 a x ÷ i n k i i ij j

 è j¹i ø
A sufficient condition for obtaining a solution by Jacobi’s iteration method is the diagonal dominance,
>
i.e., aii åa, n
ij i =1,2,,n
j=1
 j¹ i

i.e., in each row of A the modulus of the diagonal element exceeds the sum of the off diagonal

elements and also the diagonal elements aii ¹ 0. If any diagonal element is 0, the equations can always

be re-arranged to satisfy this condition.

Remark: The solution to any iterative procedure is terminated when a predetermined tolerance is

( ) ()
reduced by all xi’s. This condition is usually set as xi k+1 - xi k £Tol .

Example 1
Solve the following system of equation by Jacobi’s iteration method
4x1 + x2 + x3 + x4 =1 + x3
− 2x +5x2 x1 − x4 =0 + 4x3 + x4
− x2 =−1

x1 + x2 − 2x3 +5x4 =0

Correct your answers to five (5) decimal places

Solution:

Numerical Methods Page 39

#WEDEY4U
 UMaT, Department of
Mathematics

Clearly, the coefficient matrix is;

1 1 −2 5
4 1 1 1
5 1 𝑤ℎ𝑖𝑐ℎ 𝑖𝑠 𝑠𝑡𝑟𝑐𝑡𝑙𝑦 𝑑𝑖𝑎𝑔𝑜𝑛𝑎𝑙𝑙𝑦 𝑑𝑜𝑚𝑖𝑛𝑎𝑛𝑡.
4 −1 𝐻𝑒𝑛𝑐𝑒 𝐽𝑎𝑐𝑜𝑏𝑖′𝑠 𝑚𝑒𝑡ℎ𝑜𝑑 𝑤𝑖𝑙𝑙 𝑐𝑜𝑛𝑣𝑒𝑟𝑔𝑒 𝑓𝑜𝑟
The coefficient matrix is A −1
− 2 1 𝑎𝑙𝑙 𝑖𝑛𝑖𝑡𝑖𝑎𝑙
= 1
𝑔𝑢𝑒𝑠𝑠 𝑥(0).

Let 𝑥(0) = [0,0,0,0] 𝑇 , |𝑥 (𝑘+1) − 𝑥 (𝑘)| ≤ 0.00005

Solving for 𝑥𝑖 in each of the equations, we obtain;

𝑥1(𝑘+1) = 14 (1 − 𝑥2(𝑘) − 𝑥3(𝑘) − 𝑥4(𝑘))

𝑥2(𝑘+1) = 15 (0 + 2𝑥1(𝑘) + 𝑥2(𝑘) − 𝑥4(𝑘))

𝑥3(𝑘+1) = 15 (−1 − 𝑥1(𝑘) + 𝑥2(𝑘) − 𝑥4(𝑘))

𝑥4(𝑘+1) = 15 (0 − 𝑥1(𝑘) − 𝑥2(𝑘) + 2𝑥3(𝑘))

Hence 𝑥1(1) = 14 (1 − 0) == 0.25

𝑥2(1) = 0

𝑥3(1) =
−

𝑥4(1) = 0

∴ 𝑋(1) = (0.25 , 0 , −0.25 , 0) 𝑇

Numerical Methods Page 40

#WEDEY4U
 UMaT, Department of Mathematics

𝑥1(2) = 14 (1 − 0 + + 0)
=

𝑥2(2) = 15 (2() + + 0) =

𝑥3(2) = 14 (−1 −+
0 − 0) = −

𝑥4(2) = 15 (−− 0 + 2() ) = −

∴ 𝑋(2) = (0.3125 , 0.15 , −0.3125 , −0.15) 𝑇

1

𝑥1(3) =4 𝑥2(4) = 15 (2 (
(1 −
++) =

𝑥2(3) = 15
(2 ( ) +−) =

− ++) = −
𝑥3(3) = 14
(−1
−+ 2
(−) ) = −
𝑥4(3) = 15
(− ∴ 𝑋(3) = (0.328125 , 0.157500 , −0.253125 , −0.217500) 𝑇

++) =

𝑥1(4) = 14 ) +−) =
(1 −

𝑥3(4) = 14 (−1 − + + )=−

Numerical Methods Page 41

#WEDEY4U
 UMaT, Department of Mathematics

𝑥4(4) = 15 (− − + 2(− ) )=−

∴ 𝑋(4) = (0.328281 , 0.138375 , −0.238281 , −0.198375) 𝑇

𝑥1(5) = 14 (1 − 0.138375 + 0.238281 + 0.198375) = 0.324570

𝑥2(5) = 15 (2(0.328281) + 0.238281 − 0.198375) = 0.139294

𝑥3(5) = 14 (−1 − 0.328281 + 0.138375 + 0.198375) = −0.247883

𝑥4(5) = 15 (−0.328281 − 0.138375 − 2(−0.238281)) = −0.188644

∴ 𝑋(5) = (0.324570 , 0.139294 , − 0.247883 , −0.188644 ) 𝑇

𝑥1(6) = 14 (1 − 0.139294 + 0.247883 + 0.188644) = 0.324308

𝑥2(6) = 15 (2(0.324495) + 0.247883 − 0.188644) = 0.141676

𝑥3(6) = 14 (−1 − 0.324495 + 0.139294 + 0.188644) = −0.247677

𝑥4(6) = 15 (−0.324495 − 0.139294 + 2(−0.247883)) = −0.191926

∴ 𝑋(6) = ( 0.324308, 0.141676 , −0.247677 , −0.191926) 𝑇

𝑥1(7) = 14 (−0.141676 + 0.247677 + 0.191926) = 0.324852

𝑥2(7) = 15 (2(0.324308) + 0.247677 − 0.191926 ) = 0.141140

𝑥3(7) = 14 (−1 − 0.324308 + 0.141676 + 0.191926 ) = −0.247677

𝑥4(7) = 15 (− 0.324308 − 0.141676 + 2(−0.247677)) = −0.192860

Numerical Methods Page 42

#WEDEY4U
 UMaT, Department of Mathematics

∴ 𝑋(7) = (0.324852 , 0.141140 , −0.247677 , −0.192860) 𝑇

𝑥1(8) = 14 (1 − 0.141770 + 0.247677 − 0.192860) = 0.324842

𝑥2(8) = 15 (2(0.324852) + 0.247677 − 0.191926 ) = 0.140904

𝑥3(8) = 14 (−1 − 0.324308 + 0.141676 + 0.192860 ) = −0.247706

𝑥4(8) = 15 (− 0.324852 − 0.141170 + 2(−0.247677)) = −0.192275

∴ 𝑋(8) = (0.324842, 0.140904 , −0.247706, −0.192275 ) 𝑇

𝑥1(9) = 14 (1 − 0.140904 + 0.247706 − 0.192275) = 0.324769

𝑥2(9) = 15 (2(0.324842) + 0.247706 − 0.192275 ) = 0.141023

𝑥3(9) = 14 (−1 − 0.324842 + 0.140904 + 0.192275 ) = −0.249716

𝑥4(9) = 15 (− 0.324542 − 0.140904 − 2(0.247706)) = −0.192231

∴ 𝑋(9) = (0.324769, 0.141023, −0.247916 , −0.192231 ) 𝑇

𝑥1(10) = 14 (1 − 0.141023 + 0.247916 + 0.192231) = 0.324781

𝑥2(10) = 15 (2(0.324769) + 0.247916 − 0.192231 ) = 0.141044

𝑥3(10) = 14 (−1 − 0.324769 + 0.41023 + 0.192231 ) = −0.247853

𝑥4(10) = 15 (− 0.324769 − 0.141023 − 2(−0.247916)) = −0.192325

∴ 𝑋(10) = ( 0.324781 , 0.141044 , −0.247853, −0.192325 ) 𝑇

Numerical Methods Page 43

#WEDEY4U
 UMaT, Department of Mathematics

3.2 Gauss Seidel Iteration Method

A simple modification to Jacobi’s iteration method is given by Gauss-Seidel method. It is also known
as the Liebmann method or the method of successive displacement.

Step 1 (Gauss-Seidel method): Determination of first approximation using

.

Step n 1: In general, if  x1(n),x2(n),,xn(n) are a system of nth approximations, then the next

approximation is given by the formula

Equation (6) can briefly described as follows:

+1) 1 æ - åi-1 a x( +1) - ån a x( )ö

Numerical Methods Page 44

#WEDEY4U
 UMaT, Department of Mathematics

xi(r = aii èçbi j=1 ij jr j=i+1 ij jr ø÷ (i =1,2,,n r =0,1,2,)

Remark: We note the difference between Jacobi’s method and Gauss-Seidel method.

(a) Gauss Seidel Method (b) Jacobi Method

Remark

5x1 − 2x2 + 2x3 =10

Given this system of linear equations 2x1 − 4x2 − x3 = −7

3x1 − x2 + 6x3 =12

Numerical Methods Page 45

#WEDEY4U
 UMaT, Department of Mathematics

Jacobi Gauss − Seidel

2 old 2 old 10 2 old 2 old 10
x1new = x2 − x3 + x1new = x2 − x3 +
5 5 5 5 5 5
2 1 old 7 2 new 1 old 7
x2new = x1old − 4 x3 + 4 x2new = 4 x1 − 4 x3 + 4
4
3 1 old 12 3 new 1 new 12
x3new =− x1old + 6 x2 +6 x3new =− 6 x1 + 6 x2 +6
6

Note:
 This method is very useful with less work for the given systems of equation whose augmented
matrix have a large number of zero elements.

 We say a matrix is diagonally dominant if the numerical value of the leading diagonal element
in each row is greater than or equal to the sum of the numerical values of the other elements in
that row.

 For the Gauss-Seidal method to coverage quickly, the coefficient matrix must be diagonally
dominant. If it is not so, we have to rearrange the equations in such a way that the coefficient
matrix is diagonally dominant and then only we can apply Gauss-Seidal method.
Example 1

Determine if the Guass-Seidel method will converge for the system of equations
5x1 − x2 + x3 − x4 + x5 = 1
− x1 +8x2 + 2x3 + x4 −2x5 = −1
1
2x1 − x2 +6x3 − x4 − x5 =
−1
x1 + 2x2 +3x3 −9x4 −2x5 =
1
− x2 − x3 + x4
2x1 −6x5 =
And solve the system correct to two decimal places if the method converges.

Solution:

Clearly, the coefficient matrix is;

Numerical Methods Page 46

#WEDEY4U
 UMaT, Department of Mathematics

−1 1
5 1
2
−1 86 −2
A= 2 −1 3 −1 −1
−1 1
12 −1 −2
2 −1 −9 −6
1
is strictly diagonally dominant. Hence Jacobi’s method converges for arbitrary initial guesses.
Consequently, the Gauss –Seidel method also converges for arbitrary initial guesses. The algorithm
for the method is;

( k 1) 1 i1 n

a
xi + = ii bi − j−=1 a xij(jk+1) − j i = +1a xij (jk)

Gives;

𝑥1(𝑘+1) = 15 [1 + 𝑥2(𝑘) − 𝑥3(𝑘) + 𝑥4(𝑘) − 𝑥5(𝑘)]

𝑥2(𝑘+1) = 18 [−1 + 𝑥1(𝑘+1) − 2𝑥3(𝑘) + 𝑥5(𝑘)]

𝑥3(𝑘+1) = 16 [1 − 2𝑥1(𝑘+1) + 𝑥2(𝑘+1) + 𝑥4(𝑘) + 𝑥5(𝑘)]

𝑥4(𝑘+1) = 1 9 [1 + 𝑥1(𝑘+1) + 2𝑥2(𝑘+1) + 3𝑥3(𝑘+1) − 2𝑥5(𝑘)]

𝑥1(𝑘+1) = 16 [−1 + 2𝑥1(𝑘+1) − 𝑥2(𝑘+1) − 𝑥3(𝑘+1) + 𝑥4(𝑘+1)]

Using the initial guess [0,0,0,0] 𝑇 we have:

𝑥1(1) = 15 (1 + 0 − 0 + 0 − 0) = 0.20000

𝑥2(1) = 18 (−1 + .2 − 2(0) − 0 + 2(0)) = −0.10000

Numerical Methods Page 47

#WEDEY4U
 UMaT, Department of Mathematics

𝑥3(1) = 16 (1 − 2(0.2) − 0.10000 + 0 + 0) = = 0.08333

𝑥4(1) = 1 9 (−1 + (0.2) − 2(0.1) + 3(0.08333) − 2(0) = 0.13889

𝑥4(1) = 16 (−1 + 2(0.2) + 0.1 − + − 0) = − = −0.7407

∴ 𝑋(1) = ( 0.20000 , −0.10000, 0.08333, 0.113889, −0.07407 ) 𝑇

𝑥1(2) = 15 (1 − 0.1 − 0.08333 + 0.13889 + 0.07407) = 0.205930

𝑥2(2) = 18 (−1 + 0.20593 − 2(0.08333) − 0.13889 + 2(−0.07407) = −0.15597

𝑥3(2) = 16 (1 − 2(0.2093) − 0.15597 + 0.13889 − 0.07407 ) = −0.08283

𝑥4(2) = 1 9 (1 + − 0.20593 + 2(−0.15597) + 3(0.08283) − 2(0.07407)) = 0.14340

𝑥5(2) = 16 (−1 + 2(0.20593) + 0.15597 − 0.8333 + 0.14340) = −0.06193

∴ 𝑋(2) = (0.220593, −0.15597, 0.08283, 0.114340, −0.06193) 𝑇

𝑥1(3) = 15 (1 − 0.15597 − 0.08283 + 0.14340 + 0.06193) = 0.193307

𝑥2(3) = 18 (−1 + 0.19331 − 2(0.08283) − 0.14340 + 2(−0.06193)) = −0.154950 𝑥3(3) = 16

(1 − 2(0.19331) + (−0.15495) + 0.14340 − 0.06193 ) = −0.08998

𝑥4(3) = 1 9 (1 + 0.19331 + 2(−0.15495) + 3(0.08998) + 2(0.06193)) = 0.14191

𝑥5(3) = 16 (−1 + 2(0.19331) + 0.15495 − 0.08998 + 0.14191) = −0.067750

∴ 𝑋(3) = ( 0.19331, −0.15495, 0.8998 ,0.14191, −0.06775) 𝑇

Numerical Methods Page 48

#WEDEY4U
 UMaT, Department of Mathematics

𝑥1(4) = 15 (1 − 0.15495 − 0.08998 + 0.14191 + 0.06775) = 0.19295

𝑥2(4) = 18 (−1 + 0.19295 − 2(0.08998 − 0.14191 2(−0.06775) = −0.158053

𝑥3(4) = 16 (1 − 2(0.19295) − 0.15803 + 0.14191 − 0.06775) = 0.08836

𝑥4(4) = 1 9 (1 + 0.19295 − 2(0.15805) + 3(0.08836) + 2(0.06775) ) = 0.14191

𝑥5(4) = 16 (−1 + 2(0.19295) + 0.15805 − 0.08836 + 0.14194) = −0.06708

∴ 𝑋(4) = ( 0.19295, −0.15805 , 0.08836 , 0.14194, −0.06708) 𝑇

𝑥1(5) = 15 (1 − 0.15805 − 0.08836 + 0.14194 + 0.06708) = 0.19252

𝑥2(5) = 18 (−1 + 0.19252 − 2(0.08836) − 0.14194 − 2(0.06708)) = −0.15754

𝑥3(5) = 16 (1 − 2(0.19252) − 0.15754 ) = −0.08998

𝑥4(5) = 1 9 (1 + 0.19331 + 2(−0.15495) + 3(0.08998) + 2(0.06193)) = 0.14191

𝑥5(5) = 16 (−1 + 2(0.19331) + 0.15495 − 0.08998 + 0.14191) = −0.067750

∴ 𝑋(5) = ( 0.19252, −0.15754, 0.8871, 0.14197 , −0.06736) 𝑇

𝑥1(6) = 15 (1 − 0.15754 − 0.08871 + 0.14197 + 0.06736) = 0.19262

𝑥2(6) = 18 (−1 + 0.19262 − 2(0.08871) − 0.14197 + 2(−0.06736)) = −0.15769

𝑥3(6) = 16 (1 − 2(0.19262) − 0.15769 + 0.14197 − 0.06736 ) = 0.08861

𝑥4(6) = 1 9 (1 + 0.19262 − 2(−0.15769) + 3(0.08861) + 2(0.06736)) = 0.14198

Numerical Methods Page 49

#WEDEY4U
 UMaT, Department of Mathematics

𝑥5(6) = 16 (−1 + 2(0.19262) + 0.15736 − 0.08861 + 0.14198) = −0.06728

∴ 𝑋(6) = ( 0.19262, − 0.15736, 0.08861, 0.14198, − 0.06728 ) 𝑇

𝑥1(7) = 15 (1 − 0.15769 − 0.08861 + 0.14198 + 0.06728) = 0.19259

𝑥2(7) = 18 (−1 + 0.19259 − 2(0.08861) − 0.14198 − 2(0.06728)) = −0.15765

𝑥3(7) = 16 (1 − 2(0.19259) − 0.15765 + 0.14198 − 0.06728 ) = 0.08865

𝑥4(7) = 1 9 (1 + 0.19259 − 2(0.15765) + 3(0.08865) + 2(0.06728)) = 0.14198

𝑥5(7) = 16 (−1 + 2(0.19259) + 0.15765 − 0.08865 + 0.14198) = −0.06731

∴ 𝑋(7) = ( 0.19259, − 0.15765, 0.08865, − 0.06731 ) 𝑇

𝑆𝑖𝑛𝑐𝑒 |𝑥𝑖(7) − 𝑥𝑖(6)| < 0.00005 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑖, 𝑖𝑡 𝑓𝑜𝑙𝑙𝑜𝑤𝑠 𝑡ℎ𝑎𝑡 𝑥(4) 𝑖𝑠 𝑡ℎ𝑒 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑡𝑜 𝑡ℎ𝑒

𝑠𝑦𝑠𝑡𝑒𝑚 𝑡𝑜 𝑓𝑜𝑢𝑟 𝑑𝑒𝑐𝑖𝑚𝑎𝑙 𝑝𝑙𝑎𝑐𝑒𝑠 𝑓𝑜𝑟 𝑖 = 1, 2, . . . 5.

∴ 𝑋(7) = ( 0.19259, − 0.15765, 0.08865, − 0.06731 ) 𝑇

Example 2

Solve the following system of equation by Jacobi’s iteration method

-5x1 +12x3 =80

4x1 - x2 - x3 =-2

6x1 +8x2 =45

Correct your answers to four (4) decimal places

Numerical Methods Page 50

#WEDEY4U
 UMaT, Department of Mathematics

Solution:

Clearly, the coefficient matrix is not diagonally dominant!!!

Order of the equation can be important. Rearrange the equations to ensure convergence.
4x1 − x2 x1 = (x2 +
6x1 + 8x2 8
x3 − 2) 4
12
−5x1
x2 = (45
− 6x1) x3

Rearrange = (80 +
5x1)

Assume x1 = x2 = x3 = 0

x1 =(0 + −0 2) 4 =−0.5

First iteration
− x3 = − 2 4 −1 −1 x2 = 45− 6 (−0.5) 8 = 6.0

= 45 6 8 0
+ 12x3 = 80 0 x3 = 80 + 5 (−0.5) 12 = 6.4583
−5 12

x1 = − + +( 2 6 6.4582) 4 = 2.6146

Second iteration x2 = 45− 6 2.6146() 8 = 3.6641

x3 = 80 + 5 2.6146() 12 = 7.7561

x1 = − +( 23.6641+ 7.7561) 4 = 2.3550

Numerical Methods Page 51

#WEDEY4U
 UMaT, Department of Mathematics

Third iteration x2 = 45− 6 2.3550() 8 = 3.8587

x3 = 80 + 5 2.3550() 12 = 7.6479

x1 = − +( 23.8587 + 7.6479) 4 = 2.3767

Fourth iteration x2 = 45− 6 2.3767() 8 = 3.8425

x3 = 80 + 5 2.3767() 12 = 7.6569
5th: x1 = 2.3749, x2 = 3.8439, x3 = 7.6562
6th: x1 = 2.3750, x2 = 3.8437, x3 = 7.6563
7th: x1 = 2.3750, x2 = 3.8438, x3 = 7.6562

3.3 Method Of Successive Over – Relaxation (SOR)

The SOR method is similar to the Jacobi and Gauss-Seidel methods, but it uses a scaling factor to
more rapidly reduce the approximation error.

The SOR technique is one of a class of relaxation methods that compute approximations x(k) by the

formula

xi(k) =(1- w)xi(k-1) +wêébi - åi-1 aijx(jk) - ån aijx(jk-1) ùú i =1,2,

 aii ë j=1 j=i+1 û

where is the relaxation parameter.

We distinguish three cases, depending on the values of . =1 is equivalent to the Gauss Seidel

method, 1 is called under-relaxation, and 1 is called over-relaxation. This procedure is

illustrated as.

Numerical Methods
x1new = (1− )x1old + (b − a x
1 12 2old − a x13 3old − a x14 4old ) a11 Page 52
#WEDEY4U
 UMaT, Department of Mathematics

xx2newnew == ((11−− ))xx32oldold ++ ((bb 32 −−a xa x31 121 1newnew −−a xa x32 223

3newold −−a xa x24 434 4oldold ))aa2233

3

xnew = (1− )x4old + (b − a x

4 41 1new − a x42 2new − a x43 3new ) a44
4

ALTERNATIVELY (VECTOR FORM)

w
(
Ax = L+D +U x =b Þ ) (D +L)x =-wUx +wb
(D - wL)x =((1- w)D - wU)x +wb x (
= D - wL )
((1- w)D +wU)x +w D - wL b =Cx +d new -1 old (
) -1 old

Remarks:

 It can be known that the SOR method converges arbitrary initial guess, 0 <w<2. Also it can be

shown that if 0 <w<2 and A is strictly diagonally dominant, then the SOR method converges

for all initial guesses.

( )
 In order to avoid inverting the matrix D-wL , the SOR method is applied to the form:

(D - w)x (k+1)
={(1- w)D +wU} x (k)
+wb, and x(k+1) is then found by forward elimination,

( )
since D-wL is lower triangular with non-zero elements aii on its diagonal.

Example 1
Numerical Methods Page 53
#WEDEY4U
 UMaT, Department of Mathematics

Determine whether the SOR method can be used to solve the linear system
3x1 − x2 + x3 =2
x + 4x2 − x3 =0
1 + x4
− 2x1 + x2 − 6x3 − x4 =−1
x2
+ x3 − 4x4 = 0
Solve using the relaxation parameter w=1.2 and give your answer correct to 2-decimal places.

Solution:

The coefficient matrix of the system is

3 0

−1 1
1
4 −1
A= −12 −1
1 −6
1 1
0 −
4

Since A (a ) is strictly diagonally dominant, we know that for any w set for 0 <w<2, the SOR
=
ij

method converges for all initial guesses. Hence the method can be applied to solve the given system.
The method may be applied in vector form or in component form. Computing by component form,
we have:

( k 1) 1
xi + = −( )xi(k) + bi − i−1 a xij(jk+1) − n a xij(jk)
aii j=1 j i= +1

Hence with = 1.2 𝑎𝑛𝑑 𝑥(0) = (0 ,0 ,0 , 0) 𝑇. We have;

𝑥1(𝑘+1) = (1 − 1.2)𝑥1(𝑘) + (2 + 𝑥2(𝑘) − 𝑥3(𝑘) )

𝑥2(𝑘+1) = (1 − 1.2)𝑥2(𝑘) + (−𝑥1(𝑘+1) + 𝑥3(𝑘) − 𝑥4(𝑘))

Numerical Methods Page 54

#WEDEY4U
 UMaT, Department of Mathematics

𝑥3(𝑘+1) = (1 − 1.2)𝑥3(𝑘) + (2𝑥1(𝑘+1) − 𝑥2(𝑘+1) + 𝑥4(𝑘) − 1)

𝑥4(𝑘+1) = (1 − 1.2)𝑥4(𝑘) + (−𝑥2(𝑘+1) + 𝑥3(𝑘+1))

𝑥1(𝑘+1) = (−0.2)𝑥1(𝑘) + (2 + 𝑥2(𝑘) − 𝑥3(𝑘) )

𝑥2(𝑘+1) = (−0.2)𝑥2(𝑘) + (−𝑥1(𝑘+1) + 𝑥3(𝑘) − 𝑥4(𝑘))

𝑥3(𝑘+1) = (−0.2)𝑥3(𝑘) + (2𝑥1(𝑘+1) − 𝑥2(𝑘+1) + 𝑥4(𝑘) − 1)

𝑥4(𝑘+1) = (1 − 1.2)𝑥4(𝑘) + (−𝑥2(𝑘+1) + 𝑥3(𝑘+1))

𝑥1(1) = 0.4(0.2) = 0.8

𝑥2(1) = 0.3(0.8) = −0.24

𝑥3(1) = −0.2(−1 + 2(0.8 + 0.24) = −0.168

𝑥4(1) = −0.3( 0.24 + 0.168) = −0.1224

∴ 𝑥(1) = [ 0.8, −0.24 , − 0.168, − 0.1224] 𝑇 𝑥1(2) =

−0.2(0.8) + 0.4(2 − 0.24 + 0.168) = 0.6112

𝑥2(2) = (−0.2)(−0.24) − 0.3(−0.6112 − 0.168 + 1224 ) = −0.1490

𝑥3(2) = −0.2(0.168) − 0.2(−1 + 2(0.6112) + 0.1490 − 0.1224) = −0.0162

𝑥4(2) = −0.2(−0.1224) − 0.3(−0.6112 + 0.0162) = 0.2030

∴ 𝑥(2) = [ 0.6112, − 0.1490, −0.0162, 0.2030] 𝑇

𝑥1(3) = −0.2(0.6112) + 0.4(2 − 0.1490 + 0.0162) = 0.6246

𝑥2(3) = (−0.2)(−0.1490) + 0.3(−0.6246 − 0.0162 − 0.2030) = −0.2233

Numerical Methods Page 55
#WEDEY4U
 UMaT, Department of Mathematics

𝑥3(3) = −0.2(0.0162) − 0.2(−1 + 2(0.6112) + 0.1490 + 0.2030) = −0.1116

𝑥4(3) = (−0.2)(−0.0162) − 0.3(0.1490 + 0.0162) = −0.0902

∴ 𝑥(3) = [ 0.6264, −0.2233, − 0.1116, − 0.0902] 𝑇

𝑥1(4) = −0.2(0.6264) + 0.4(2 − 0.2233 + 0.1116) = 0.6304

𝑥2(4) = (−0.2)(−0.2233) + 0.3(−0.6304 − 0.1116 + 0.0902) = −0.1509

𝑥3(4) = −0.2(−0.1116) − 0.2(−1 + 2(0.6304) + 0.1509 − 0.0902) = −0.0420

𝑥4(4) = −0.2(−0.0902) − 0.3(0.1509 + 0.0420) = −0.0398

∴ 𝑥(4) = [ 0.6304, −0.1509, −0.0420, −0.0398] 𝑇

𝑥1(5) = −0.2(0.6304) + 0.4(2 − 0.1509 + 0.0420) = 0.6304

𝑥2(5) = (−0.2)(−0.1509) + 0.3(−0.6304 − 0.0420 + 0.0398) = −0.1596

𝑥3(5) = −0.2(−0.0420) − 0.2(−1 + 2(0.6304) + 0.1596 − 0.0398) = −.0677

𝑥4(5) = −0.2(−0.0398) − 0.3(0.1596 + 0.0677) = −0.0602

∴ 𝑥(5) = [ 0.6304, − 0.1596, − 0.0677, − 0.0602] 𝑇 𝑥1(6) =

−0.2(0.6304) + 0.4(2 − 0.1596 + 0.0677) = 0.6372

𝑥2(6) = (−0.2)(−0.1596) + 0.3(−0.6372 − 0.0677 + 0.0602) = −0.1615

𝑥3(6) = −0.2(−0.0677) − 0.2(−1 + 2(0.6372) + 0.1615 − 0.0602) = −.0616

𝑥4(6) = −0.2(−0.0602) − 0.3(0.1615 + 0.0616) = −0.0549

∴ 𝑥(6) = [ 0.6372, − 0.1615, −0.0616, − 0.0549] 𝑇

𝑥1(7) = −0.2(0.6372)
Numerical Methods + 0.4(2 − 0.1615 + 0.0616) = 0.6326 Page 56
#WEDEY4U
 UMaT, Department of Mathematics

𝑥2(7) = (−0.2)(−0.1615) + 0.3(−0.6326 − 0.0616 + 0.0549) = −0.1595

𝑥37 = −0.2(−0.0616) − 0.2(−1 + 2(0.6326) + 0.1595 − 0.0549) = −.0616

𝑥4(7) = −0.2(−0.0549) − 0.3(0.1595 + 0.0616) = −0.0554

∴ 𝑥(7) = [ 0.6326, − 0.1595, − 0.0616, − 0.0554] 𝑇

𝑆𝑖𝑛𝑐𝑒 |𝑥𝑖(7) − 𝑥𝑖(6)| < 0.005 𝑓𝑜𝑟 𝑖 = 1,2,3,4 𝑖𝑡 𝑓𝑜𝑙𝑙𝑜𝑤𝑠 0.005 𝑓𝑜𝑟 𝑖 = 1,2,3,4

𝑖𝑡 𝑓𝑜𝑙𝑙𝑜𝑤𝑠 𝑡ℎ𝑎𝑡 𝑐𝑢𝑟𝑟𝑒𝑛𝑡 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑥(7) 𝑖𝑠 𝑐𝑜𝑟𝑟𝑒𝑐𝑡 𝑡𝑜 2 𝑑𝑒𝑐𝑖𝑚𝑎𝑙 𝑝𝑙𝑎𝑐𝑒𝑠.

∴ 𝑥(7) = [ 0.63, − 0.16, − 0.06, − 0.05] 𝑇

EXERCISES (Try all)

1. Apply Gauss Seidel iteration method to solve:

10x +2y +z =9
2x +20y - 2z =-44
-2x +3y +10z =22
2. Apply Gauss Seidel iteration method to solve:

1.2x +2.1y +4.2z =9.9

5.3x +6.1y +4.7z =21.6
9.2x +8.3y +z =15.2

3. Apply Jacobi’s iteration method to solve:

5x - y +z =10
2x - y +z =10
x +y +5z =-1

4. Apply Jacobi’s iteration method to solve:

Numerical Methods Page 57

#WEDEY4U
 UMaT, Department of Mathematics

5x +2y +z =12 x
+4y +2z =15
x +2y +5z =20

CHAPTER 4
Numerical Methods for Non-Linear Equations in One Variable
4.0 Introduction
A numerical iteration method or simply iteration method is a mathematical procedure that
generates a sequence of improving approximate solutions for a class of problems. A specific way of
implementation of an iteration method, including the termination criteria, is called an algorithm of the
iteration method. In the problems of finding the solution of an equation an iteration method uses an
initial guess to generate successive approximations to the solution.
Since the iteration methods involve repetition of the same process many times, computers can act well
for finding solutions of equation numerically. Some of the iteration methods for finding solution of
Numerical Methods Page 58
#WEDEY4U
 UMaT, Department of Mathematics

equations involves (1) Bisection method, (2) Method of false position (Regula-falsi Method), (3)
Newton-Raphson method.
A numerical method to solve equations may be a long process in some cases. If the method leads to
value close to the exact solution, then we say that the method is convergent. Otherwise, the method is
said to be divergent.
Solution of Algebraic and Transcendental Equations
One of the most common problems encountered in engineering analysis is that given a function f x ,

find the values of x for which f x 0 . The solution (values of x) is known as the roots of the equation f

x 0 , or the zeroes of the function f x . The roots of equations may be real or complex.

In general, an equation may have any number of (real) roots, or no roots at all. For example, sinx x 0
has a single root, namely, x 0, whereas tanx x 0 has infinite number of roots
x 0, 4.493, 7.725,... .

Algebraic and Transcendental Equations

If then the equation f x 0 is called an algebraic equation.

An example is 2 x3 x 4 0 . If f x 0 contains trigonometric, or exponential or logarithmic
functions then fx 0 is called transcendental equation. Examples of transcendental
, equations
are and 7x3 log 3x 6 3ex cosx 0.

There are two types of methods available to find the roots of algebraic and transcendental equations

of the form f x 0.

1. Direct Methods: Direct methods give the exact value of the roots in a finite number of steps.
We assume here that there are no round off errors. Direct methods determine all the roots at the same
time.
2. Indirect or Iterative Methods: Indirect or iterative methods are based on the concept of
successive approximations. The general procedure is to start with one or more initial approximation
to the root and obtain a sequence of iterates xk which in the limit converges to the actual or true solution
to the root. Indirect or iterative methods determine one or two roots at a time. The indirect or iterative
methods are further divided into two categories: bracketing and open methods. The bracketing
Numerical Methods Page 59
#WEDEY4U
 UMaT, Department of Mathematics

methods require the limits between which the root lies, whereas the open methods require the initial
estimation of the solution. Bisection and False position methods are two known examples of the
bracketing methods. Among the open methods, the Newton-Raphson is most commonly used. The
most popular method for solving a non-linear equation is the Newton-Raphson method and this
method has a high rate of convergence to a solution.
In this chapter and in the coming chapters, we present the following indirect or iterative methods with
illustrative examples:

4.1 Bisection Method

The bisection method is one of the bracketing methods for finding roots of an equation. For a given a

function , guess an interval which might contain a root and perform a number of iterations,
where, in each iteration the interval containing the root is get halved. The bisection method is based
on the intermediate value theorem for continuous functions.

Intermediate Value Theorem for Continuous Functions

If is a continuous function and and have opposite signs, then at least one root lies in between a and

b. If the interval is small enough, it is likely to contain a single root. i.e., an interval must contain a

zero of a continuous function f if the product .

Geometrically, this means that if , then the curve f has to cross the x-axis at some point in between a
and b.

Numerical Methods Page 60

#WEDEY4U
 UMaT, Department of Mathematics

Algorithm: Bisection Method

Suppose we want to find the solution to the equation , where f is continuous. Given a function

continuous on an interval and satisfying .

For until termination do:

Compute .

If , accept xn as a solution and stop.

Else continue.
If , a root lies in the interval .

Set .

If , a root lies in the interval .

Set .

Then for some x in .

Test for termination.

Criterion for termination

A convenient criterion is to compute the percentage error defined by where is

the new value of xr . The computation can be terminated when becomes less than a

Numerical Methods Page 61

#WEDEY4U
 UMaT, Department of Mathematics

prescribed tolerance, say . In addition, the maximum number of iterations may also be specified in
advance.

Some other termination criteria are as follows:

• Termination after N steps (N given, fixed)

• Termination ( if given)

• Termination ( if given)

In this chapter our criterion for termination is terminate the iteration process after some finite steps.
However, we note that this is generally not advisable, as the steps may not be sufficient to get an
approximate solution.

Example 1 Choose an interval of length 1 unit which contains the root of x −sin x =1. Apply 3 steps
of the bisection method to improve this approximation to the root.
Solution Making a rough sketch of y = −x 1 and y = sin x shows a root near 2.

We could choose a =1.5, b = 2.5 and f x( ) = x − −1 sin x.

f a( ) f b( ) f ( p)
Step a b p

1 1.5 0 2.5 0 2 0
2 1.5 0 2 0

Merits of bisection method

a) The iteration using bisection method always produces a root, since the method brackets the
root between two values.
b) As iterations are conducted, the length of the interval gets halved. So one can guarantee the
convergence in case of the solution of the equation.
c) the Bisection Method is simple to program in a computer.

Demerits of bisection method

a) The convergence of the bisection method is slow as it is simply based on halving the interval.
Numerical Methods Page 62
#WEDEY4U
 UMaT, Department of Mathematics

b) Bisection method cannot be applied over an interval where there is a discontinuity.

c) Bisection method cannot be applied over an interval where the function takes always values
of the same sign.
d) The method fails to determine complex roots.

e) If one of the initial guesses a0 or b0 is closer to the exact solution, it will take larger number

of iterations to reach the root.

EXERCISES

Find a real root of the following equations by bisection method.

4.2 Regula Falsi Method (Method of False Position)

This method is also based on the intermediate value theorem. In this method also, as in bisection
method, we choose two points a and b such that and are of opposite signs (i.e.,
n n

) . Then, intermediate value theorem suggests that a zero of f lies in between an

, if
and bn f is a continuous function.

Numerical Methods Page 63

#WEDEY4U
 UMaT, Department of Mathematics

Algorithm: False Position Method

Given a function continuous on an interval and satisfying .

For until termination do:

Compute

If , accept xn as a solution and stop.

Else continue.

If , a root lies in the interval .

Set .

Numerical Methods Page 64

#WEDEY4U
 UMaT, Department of Mathematics

If , a root lies in the interval .

Set .

Then for some x in .

Example: Using regula-falsi method, find a real root of the equation, , near .

Here note that and . Hence , so by intermediate value theorem a

root lies in between 0 and 1. We search for that root by regula falsi method and we will get an
approximate root.

Then
Numerical Methods Page 65
#WEDEY4U
 UMaT, Department of Mathematics

Theoretical Exercises with Answers:

1. What is the difference between algebraic and transcendental equations?

Ans: An equation is called an algebraic equation if the corresponding is a

polynomial, while, is called transcendental equation if the contains

trigonometric, or exponential or logarithmic functions.

2. Why do we use numerical iterative methods for solving equations?

Ans: As analytic solutions are often either too tiresome or simply do not exist, we need to find
NumericalanMethods
approximate method of solution. This is where numerical analysis comes into the picture.
Page 66
#WEDEY4U
 UMaT, Department of Mathematics

3. Based on which principle, the bisection and regula-falsi method is developed?

Ans: These methods are based on the intermediate value theorem for continuous functions:

stated as, “If f is a continuous function and and have opposite signs, then at

least one root lies in between a and b. If the interval (a, b) is small enough, it is likely to contain
a single root.”

4. What are the advantages and disadvantages of the bracketing methods like bisection and
regulafalsi?

Ans: (i) The bisection and regula-falsi method is always convergent. Since the method brackets
the root, the method is guaranteed to converge. The main disadvantage is, if it is not possible
to bracket the roots, the methods cannot be applicable. For example, if f (x) is such that it
always takes the values with same sign, say, always positive or always negative, then we
cannot work with bisection method. Some examples of such functions are

• , which take only non-negative values and

• , which take only non-positive values.

EXERCISES

Find a real root of the following equations by false position method.

Numerical Methods Page 67

#WEDEY4U
 UMaT, Department of Mathematics

4.3 Newton – Raphson Method

The Newton-Raphson method, or Newton Method, is a powerful technique for solving equations
numerically. Like so much of the differential calculus, it is based on the simple idea of linear
approximation.

Consider , where f has continuous derivative . From the figure we can say that at ;

which means that a is a solution to the equation . In order to find

the value of a, we start with any arbitrary point xi. From figure we can see that, the tangent to the
curve f at (with slope touches the x-axis at .

Now, .

As , the above simplifies to a Newton-Raphson formula, for

Numerical Methods Page 68

#WEDEY4U
 UMaT, Department of Mathematics

The refinement on the value of the root xi is terminated by any of the following conditions.

(i) Termination after a pre-fixed number of steps

(ii) After i iterations where, for a given , or

(iii) After i iterations, . where for a given

Termination after a fixed number of steps is not advisable, because a fine approximation cannot be
ensured by a fixed number of steps.

Algorithm: Newton-Raphson Method

The steps of the Newton-Raphson method to find the root of an equation are

1. Evaluate

2. Use an initial guess of the root, xi, to estimate the new value of the root, , as

3. Find the absolute relative as approximate error

4. Compare the absolute relative approximate error with the pre-specified relative error

tolerance, . If then go to Step 2, else stop the algorithm. Also, check if the
number of iterations has exceeded the maximum number of iterations allowed. If so, one
needs to terminate the algorithm and notify the user.

Numerical Methods Page 69

#WEDEY4U
 UMaT, Department of Mathematics

The method can be used for both algebraic and transcendental equations, and it also works when
coefficients or roots are complex. It should be noted, however, that in the case of an algebraic equation
with real coefficients, a complex root cannot be reached with a real starting value.

Example. Let us find an approximation to 5 to ten decimal places.

Note that 5 is an irrational number. Therefore the sequence of decimals which defines 5 will not stop.
Clearly 5 is the only zero of on the interval [1, 3].

Solution. Let be the successive approximations obtained through Newton's method. We have Let us

start this process by taking

EXERCISES (Try all)

5. Approximate the real root to two four decimal places of

6. Approximate to four decimal 3

3 places

7. Find a positive root of the equation correct to 4 places of decimals.

(Choose )
Numerical Methods Page 70
#WEDEY4U
 UMaT, Department of Mathematics

8. Explain how to determine the square root of a real number by N R method and using it

determine 3 correct to three decimal places.

9. Find the value of 2 correct to four decimals places using Newton Raphson method.
10. Use the Newton-Raphson method, with 3 as starting point, to find a fraction that is within

of 10 .

3
11. Design Newton iteration for the cube root. Calculate 7 , starting from and

performing 3 steps.

12. Calculate 7 by Newton’s iteration, starting from and calculating x1, x2, x3.

Compare the results with the value

th
13. Design a Newton’s iteration for computing k root of a positive number c.
14. Find all real solutions of the following equations by Newton’s iteration method.

4.4 Secant Method

We have seen that the Newton-Raphson method requires the evaluation of derivatives of the function
and this is not always possible, particularly in the case of functions arising in practical problems. In
the secant method, the derivative at xn is approximated by the formula

which can be written as

where . Hence, the Newton-Raphson formula becomes

Numerical Methods Page 71

#WEDEY4U
 UMaT, Department of Mathematics

It should be noted that this formula requires two initial approximations to the root.

We have
Example: Find a real root of the equation using secant method. Let the two initial
approximations be given by and .

Also,

Again

EXERCISES (Try all)

1. Determine the real root of the equation using the secant method. Compare your result
with the true value of x 0.567143…
2. Use the secant method to determine the root, lying between 5 and 8, of the equation
.
3. Use the secant method to solve the following equations, with the starting values given, correct
to 3 decimal places:

Numerical Methods (i) cosx =0.6x x-1 =1 and x0 =1.1 Page 72

#WEDEY4U
 UMaT, Department of Mathematics

(ii) 2cosx+5x- 1=0 x-1 =-0.1 and x0 =0

(iii) x4 +x- 3=0 x-1 =1.3 and x0 =1.5

(iv) x4 - 2 =0 x-1 =1.3 and x0 =1.5

Theoretical Exercises with Answers:

1. What is the difference between bracketing and open method?

Ans: For finding roots of a nonlinear equation , bracketing method requires two
guesses which contain the exact root. But in open method initial guess of the root is needed
without any condition of bracketing for starting the iterative process to find the solution of an
equation.
2. What is the importance of Secant method over Newton-Raphson method?

Ans: Newton-Raphson method requires the evaluation of derivatives of the function and this
is not always possible, particularly in the case of functions arising in practical problems. In
such situations Secant method helps to solve the equation with an approximation to the
derivative.

Numerical Methods Page 73

#WEDEY4U
 UMaT, Department of Mathematics

CHAPTER 5
Solution of Non-Linear Systems of Equations
5.0 Introduction
In this chapter we will study the solution of systems of nonlinear equations. As opposed to linear
equations, no explicit solution techniques are, in general, available for nonlinear equations, and hence
their solution completely relies on iterative method. Nonlinear equations are assumed to be those
which contain powers or products of the variables and/or transcendental functions such as x y3, −2, xy,

x y , 2( y z− )2, sin ,x e2x is called nonlinear.

A system of n equations in n unknowns, x1, x2, , x is called

n nonlinear if one or more of the equations is

nonlinear. By bringing all non zero terms to the left side of all equations any nonlinear n n system

can be put in the general form.

Numerical Methods Page 74

#WEDEY4U
 UMaT, Department of Mathematics

We may abbreviate (i) as f x( ) =0.

ë
A vector  é x1, x2,, xnùû that satisfies
fx
( )= 0 will be called a root of the nonlinear system(i).
We discuss two methods for solving such systems.

5.1 Newton-Raphson Method

Let Z = f x y( , ) and Z = g x y( , ) be two 2 2 nonlinear equations. The graphs of the equations

Z = f x y( , ) and Z = g x y( , ) are surfaces S f and Sg whose intersection is a curve C in the xyz –

plane. Setting f x y( , ) and g x y( , )to zero corresponds to intersecting S f and Sg with the xy –

plane (Z = 0). So a root x of f x y( , ) = 0 , g x y( , ) = 0 corresponds to a point at which the curve C

meets the xy – plane.

Suppose Xk = (xk , yk ) is a current approximation of Xk = x y, ( ) . Let f denote the tangent plane to S

f at the point pf (xk , g x( k )) and g denote the tangent plane to Sg at the point pg (xk , g x( k )).
Numerical Methods Page 75
#WEDEY4U
 UMaT, Department of Mathematics

Then f and g will intersect in a line L. If xk is sufficiently close to x , then L will lie close to the

curve C, and so the point X k+1 , where L meets the xy – plane should be closer to the desired root x

then xk . The z – coordinate of the point (x y z, , ) on the plane f and g can be found using the tangent

plane equation obtained in calculus:

f : z = f x y(
(ii)
k , k )+ x
g : z = g x y( k , k )+ g x y( k , k )dx + g x y( k , f x y( k , )dx
k
k )dy
x y + y f x

y( k , k )dy
Setting z to zero (ii) yields the system

x f x y( , ) y f x y( , ) dx f x( k )
=−
(iii)

dy g x( k )
x g x y( , ) y g x y( , )

Where xk = (xk , yk )

We may represent (iii) compactly as

Jdx = − f x( k ) (iv)

Numerical Methods Page 76

#WEDEY4U
 UMaT, Department of Mathematics

dx = −J−1 f x( k ) (v)

The solution of (iii) , which we denote by dxk is what must be added to xk to get the point xk+1 where

L meets the xy – plane.

Thus xk+1 = +xk dxk
xk+1 = xk − J−1 f x( k ) (vi)

If we write J−1 as f (xk )−1 then (vi) becomes

f x(
xk+1 = xk − k

) (vii)
f (xk )

Thus (vi) is the Newton-Raphson formula for solving 2 2 nonlinear system. We call the matrix J as
the Jacobian matrix at X = (x y, ) .
Example 1

Use Newton-Raphson to solve the system yex = 2 and x2 + =y 4 which has unique root located near
−0.8, 3.8 T
.

Solution:

Since Newton-Raphson can only be used to solve f x( ) = 0 , we must first convert the given equations
to an equivalent system of the form

f1(x y, ): yex − 2 = 0, f2 (x y, ): x2 + y − 4 = 0

yex − 2 f x( yex ex
J f x( ) = 2x
)= 2 y4 = 0, x +
1

−
If
Numerical x0 = −0.8, 3.8
we takeMethods T
Page 77
#WEDEY4U
 UMaT, Department of Mathematics

We now use

For x1

x1 = −3.80.8 − 1.707450−1.6 0.4493291 −1

−0.2925500.44

−0.8 −0.246927

= 3.8 − 0.044917

−0.553073
=
3.755083

−0.553073 2.159847

0.57518 −1
0.159847 x2 = 3.755083 −

−1.100146 1 0.060973

−0.553073 0.044626
= 3.755083 − 0.110333

−0.597699
=
3.644750

x3 = −3.6447500.597699 − −2.0048891.1953980.5500761 −1

0.0048890.001994

−0.597699 0.001424
= 3.644750 − 0.003697

Numerical Methods −0.599123 Page 78

#WEDEY4U
 UMaT, Department of Mathematics

=
3.641053
−0.599123 2.000005 0.549293 −1

0.000005 x4 = 3.641053 − −1.198246 1

0.000001

−0.599123 0.000002
= 3.641053 − 0.000003

−0.599125
=
3.641050
T
Since x4 − x3 0.000005 , the solution is −0.599125,3.641050 .

Example 2

Solve the system

x = 2In y (1)
y
= −xy 1 (2)
T
Using the Newton – Raphson iterative scheme till xk +1 − xk 0.00005 ; using x0 = 1.1,1.7 .

Solution:

Writing the equations in the form f x( ) = 0 we obtain

x−2Iny = 0, y− + =xy 1 0
1 −2
Thus f x( )= yx−−xy2In J f x( ) = y
−y 1−
y+1 ,
Numerical Methods x Page 79
#WEDEY4U
 UMaT, Department of Mathematics

Using Xk+1 = −Xk J−1f x( k ) to obtain x1 and x2

−1

1.1 1 −1.176471 0.038743 x1

= 1.7 − −1.7 −0.1 0.83

1.563141
=
2.126601

1.563141 1 −0.940418 −1
0.054091
x2 = 2.126601 − −2.126601

1.478762 −0.563141

= 2.094396 −0.197576

1.478762 1
x3 = 2.094396 − −2.094396 −0.954929 −1
0.000232

1.477670
= −0.478762
2.093496
−0.002717
1.477670 1
x4 = 2.093496 − −2.126601
−0.955340 −1

1.477670
−0.000001
=
2.093495
−0.477670

−0.0000023
Since x4 − x3 0.000005 , x4 is the solution. The solution is 1.477670, 2.093495 T
.

Numerical Methods Page 80

#WEDEY4U
 UMaT, Department of Mathematics

5.2 Gauss Seidel Method

A nonlinear system if equations may be worked on in the following way so as to apply the Gauss
Seidel iterative scheme to find the solutions.

1. Put each equation in the system in the form f x( ) = 0 , obtain f (x) and f (x0), where x0 is

the initial solution.

2. If on the jth column of f (x0 ), the ith element is strictly dominant, then make the element
variable the subject of the ith equation.
3. Form iterative schemes from the equations in (2), and iterate.
Example 1

Solve the system

x + In y + cosz = 0 (1)
yIn(4x) − 2z − ez +1= 0 (2)
4x + ysin z − =1 0 (3)
T

(
0.00005 ; using x0 = 1 3 , 1 3 , 6
).
ii i
Using the Gauss Seidel iterative scheme till xk( +)1 − xk( )

Solution:

The equations are all in the form f x( ) =0

1 1 1
2 x y − 2sin z

(x) = In(4x) −(2 + ez )

y
f x

4 sin z ycosz
Numerical Methods Page 81
#WEDEY4U
 UMaT, Department of Mathematics

0.866025 3 −0.25000

f (x0) = 1 0.287082 −3.688092

4 0.500000 0.288675

Now, the element a31 of the first column dominates elements in that column, we therefore make x, the
subject of the equation (3) we obtain

x = 14(1− ysin z) (a)

The second column is dominated by the equation (1) entry 3, a12 . We therefore make y the subject of
equation (1), to obtain
1

−2
(cosz+ x) (b) y e=
The entry dominates column 3 entries. We make z the subject of the equation (2)

z= (1+ yIn(4x)−e ) z
(c)
NOTE: The variable we make subject can also appear on the rhs as seen in (c) We

now form the iterative schemes from (a), (b) and (c):

1 xn+1 = (1− yn

sin(zn ))
4
1

+ (
−2 coszn+ xn+1 )
yn 1 = e
1 zn

(
zn+1 = 2 1+ y Inn+1 (4xn+1) − e )

Numerical Methods Page 82

#WEDEY4U
 UMaT, Department of Mathematics

1 1( 1 sin 6 ) = 0.208333
x1 = (1− y0 sin z0) = 1− 3

4 4
−1 +

yn+1 = e 2(cos 6 xn+1 ) = 0.410882

+ 1( ) = −0.381502
6

zn 1 = 1+ 0.410882In(4 0.208333) − e 2

Following in like manner, we obtain

x2 = 0.288244
y2 = 0.367535
z2 = 0.184741

x(2) = 0.288244, 0.367535, 0.184741 T

x(3) = 0.233122, 0.377449, −

0.114646 T
x(4) = 0.260795, 0.365168,

0.061877 T
x(5) = 0.244355, 0.370328,
T (6)
−0.036145 x = 0.253346, 0.366775,

0.020188 T
x(7) = 0.248149, 0.368600,

− 0.011566 T
x(8) = 0.251066,

0.367500, 0.006531 T
x(9) = 0.249400,

0.368104, − 0.003719 T
x(10)

= 0.250342, 0.367755, 0.002107 T

x(11)

= 0.249806, 0.367951, − 0.001197 T

x(12) = 0.250110, 0.367839, 0.000679 T

x(13) = 0.249938, 0.367902,

Numerical Methods Page 83
#WEDEY4U
 UMaT, Department of Mathematics

−0.000386 T
x(14) = 0.250035,

0.367866, 0.000219 T
x(15) = 0.249980,

0.367887, − 0.000124 T
x(16)

= 0.250011, 0.367875, 0.000071 T

x(17)

= 0.249994, 0.367882, − 0.000040 T

x(18) = 0.250004, 0.367878, 0.000023 T

x(19) = 0.249998, 0.367880, −

0.000013 T
x(20) = 0.250001, 0.367879,

0.000007 T
x(21) = 0.249999, 0.367880,

− 0.000004 T

The solution set is

0.249999, 0.367880, −0.000004 T

Numerical Methods Page 84

#WEDEY4U