222TCE100

ADVANCED NUMERICAL
METHODS
Introduction
• The limitations of analytical methods have led the engineers and scientists
to evolve graphical and numerical methods. The graphical methods ,
though simple, give results to a low degree of accuracy.

• Numerical methods can, however be derived which are more accurate.

• With the advent of high-speed digital computers and increasing demand

for numerical answers to various problems, numerical techniques have
become indispensable tools in the hands of Engineers.
Simultaneous linear system of Equations
• System of simultaneous linear equations occur in solving problems in wide variety
of disciplines including mathematics, statistics, physical, biological and social
sciences, engineering and business.
• They arise directly in solving real world problems and they also occur as a part of
the solution process for other problems. Eg: solving a system of simultaneous
non-linear equations.
• Numerical solutions of boundary value problems and initial boundary value
problems for differential equations are rich source of linear systems
 Solution of System of Linear
Equations

Direct Indirect
(Exact) (Iterative)

➢ Gauss Siedel
➢ Gauss Elimination
➢ Jacobi’s
➢ Gauss Jordan
➢ Relaxation
➢ LU Decomposition
➢
➢ Cramer’s rule
➢
➢ Matrix Inversion
• A system of m linear equations in n unknows x1,x2,x3,……xn has the
general form
a11x1 + a12x2 + ,……a1nxn = b1
a21x1 + a22x2 + ,……a2nxn = b2
.
.
am1x1 + am2x2 + ,……amnxn = bm
Can be written compactly as a matrix equation Ax = b
Where A =
• A, x, b are respectively called the coefficient matrix, solution vector,
and right-hand side column vector
• The orders of the matrices A, x, b are mxn, nx1 & mx1 respectively
• The matrix (A|b) = is

called the Augmented Matrix of the given system of Equations and

has m rows and (n+ 1) columns, i.e. of the order m x (n+1)
Elementary row operations
• The system of linear equations may be most conveniently described
in matrix form. We introduce 3 elementary row operations on general
rectangular matrices. They are
(i) Interchange of rows
(ii) Multiplication of a row by a non zero scalar
(iii) Addition of a non-zero multiple of one row to another row

We do not normally use column elementary operations in solving the

linear system of equations
 DIRECT METHODS
Gauss Elimination method Gauss Jordan method
Unknows are eliminated successively and the system is Modification of Gauss Elimination method
reduced to an upper triangular system from which
unknows are found by back substitution

General method well adopted for computer operations Elimination of unknows performed by row
transformations ultimately reducing the system to a
diagonal matrix form, i.e. each equation involving only
one unknown. From these equations, unknowns can
be obtained readily by inspection

Thus, back substitution is avoided at the cost of

additional calculations
Gauss Elimination method Gauss Jordan method
Consider the eqns: a11x+a12y+a13z = b1, a21x+a22y+a23z = b2, a31x+a32y+a33z = b3
In matrix form: Ax = b

The system is reduced to an upper triangular system by The system is reduced to diagonal form
row operations
Gauss Elimination method Gauss Jordan method

Values of x,y,z by back substitution x, y, z by direct inspection

a'11x = b1’
a‘’22y = b2’’
a‘’’33z =b3’’’
Gauss Elimination method Gauss Jordan method

Note: Note:

For a system of 10 eqns, the number of

• AX = B multiplications required for the Gauss Jordan
• Transform the coeff matrix A to upper triangular method is about 500 whereas for the Gauss
matrix by elementary row transformations only
elimination method, we need only 333
multiplications. This shows that Gauss Jordan
• The method will fail if any one of the pivots a11, a22’ method appears to be easier, but required 50%
or a33’’ becomes zero. In such cases, we rewrite the more operations than Gauss Elimination method. As
such, Gauss elimination method is preferred for
eqn in a diff order so that pivots are non zero
large systems
Gauss Elimination method Gauss Jordan method

Note: Note:
• We require max (m-1) stages of elimination to reduce • Whenever a row of zeros appear to the left of the
augmentation line and a non zero number appear to
given augmented matrix to equivalent row Echelon the right- system is inconsistent and has no solution
form. But this process may terminate at an early stage.
Eg:

A matrix is in row echelon form (REF) when it meets the following criteria:
1. All non-zero rows are above any rows of all zeros,
2. The leading entry (or pivot) of each non-zero row is to the right of the leading entry of the previous row, and
3. All entries in a column below a leading entry are zeros.
This form is useful for solving systems of linear equations and understanding linear transformations.
Gauss Elimination method Gauss Jordan method

Note: Note:
• Complete row of zeros → infinite solution

▪ If r<m & one or more of elements br+1*, br+2*,…bm* are

not zero → No solution
▪ If m ≥ n & r = n & br+1*, br+2*,…bm* are zeros, → Unique
• Otherwise → Unique solution
sol
▪ r<n, br+1*, br+2*,…bm* are all zeros → Infinite solution
• Apply Gauss-Jordan method to solve the equations shown below:

x+y+z = 9
2x-3y+4z = 13
3x+4y+5z = 40
• Apply Gauss-Jordan method to solve the following simultaneous
equations shown below:

10x+y+z = 12
2x+10y+z = 13
x+y+5z = 7
LU Decomposition/Factorisation/Triangularization

▪ LU decomposition or factorization is a fundamental technique in linear algebra.

It's used to solve systems of linear equations, invert matrices, and calculate

determinant

▪ It is a technique for breaking down a square matrix into two triangular matrices,

an upper triangular matrix and a lower triangular matrix.

LU Decomposition/Factorisation/Triangularization
• The lower triangular matrix, L, has ones on its diagonal.
The upper triangular matrix, U, has numbers above its diagonal
The product of L and U is the original matrix.

Doolittle’s method Crout’s method

Working rule – Doolittle’s method -4 STEPS
Consider the following system of equations
a11x + a12y+ a13z = b1
a21x + a22y+ a23z = b2
a31x + a32y+ a33z = b3
STEP 1
[A][X] = [B] -------(1)
Where [A] =
Working rule – Doolittle’s method -4 STEPS
STEP 2
• Let [A] = [L][U]---- (2) [L] = Lower Triangular Matrix
[U] =Upper Triangular Matrix
STEP 3
Sub (2) in (1) → [L] [U] [X] = [B] ---(3)
Multiply & find [L] & [U]
STEP 4
Put [U] [X] = [Y] ---(4)
Working rule – Doolittle’s method -4 STEPS

(3) Becomes [L] [Y] = [B] Where [Y] =

Solve to find [Y] & then put [Y] in (4) & solve to find [X]
• For Crout’s method, steps are same

• But,
• Solve the following Eqns by LU decomposition method

• x+y+z = 1, 4x+3y-z = 6, 3x+5y+3z = 4

 Indirect/Iterative methods
• A linear system Ax=b that occur in many applications can have a very
large order
• For such systems, direct methods may be too expensive either in
computation time or computer memory requirements or both.
Moreover, accumulation of round-off errors can prevent the solution
from being accurate.
• Thus, such linear systems are usually solved with iteration methods
 Indirect/Iterative methods
• The direct methods yield exact solutions after a certain amount of
fixed computations. On the other hand, an iterative method is that in
which we start from an approximation to the true solution and obtain
better and better approximations from a computation cycle, repeated
as often as may be necessary for achieving a desired accuracy.
• Thus, in an iterative method, the amount of computation depends
upon the degree of accuracy required.
• Iteration is a self-correcting process and any error made at any stage
of computation gets automatically corrected in the subsequent steps
 Jacobi’s Iteration / Method of simultaneous
replacements
• Consider the equations, a11x + a12y+ a13z = b1
a21x + a22y+ a23z = b2
a31x + a32y+ a33z = b3
If a11, a22, a33 are large as compared to other coefficients, solve for x,y,z.

• The system can be written as

 Jacobi’s Iteration / Method of simultaneous
replacements
Start with an initial approximation x0, y0, z0 for the values of x,y,z
respectively. Sub these on the RHS, the first approximations are given by
 Jacobi’s Iteration / Method of simultaneous
replacements
Sub the values of x,y,z on the RHS, the second approximations are given by

This process can be repeated until the difference between 2 consecutive

approximations are negligible
Note: In the absence of any better estimate for x0, y0, z0,these may be taken as 0
• Solve by Jacobi’s iteration method, the equations

20x + y – 2z = 17
3x + 20y – z = -18
2x -3y +20z = 25
 Gauss Seidel/ Method of successive
replacements
• Modification of Jacobi’s iteration method
• More rapidly convergent than Jacobi’s method
Consider the equations, a11x + a12y+ a13z = b1
a21x + a22y+ a23z = b2
a31x + a32y+ a33z = b3

The system can be written as

 Gauss Seidel/ Method of successive
replacements
Start with an initial approximation x0, y0, z0 for the values of x,y,z which
may be taken as zeros.
Sub y=y0, z= z0,

Now, sub x= x1, z=z0,

Next, sub x=x1, y=y1 and so on

 Gauss Seidel/ Method of successive
replacements
• i.e. as soon as a new approximation for an unknown is found, it is
immediately used in the next step

• This process of iteration is repeated until the values of x,y,z are obtained
to a desired degree of accuracy.
 Gauss Seidel/ Method of successive
replacements
Note:
• Since the most recent approximations of the unknowns are used while
proceeding to the next step, the convergence in the Gauss Seidel
method is twice as fast as in Jacobi’s method
• Jacobi and Gauss Seidel methods converge for any choice of the initial
approximations if in each equation the absolute value of the largest
coefficient is greater than the sum of the absolute values of the
remaining coefficients
• Apply Gauss Seidel iteration method to solve the equations:

20x + y – 2z = 17
3x + 20y – z = -18
2x -3y +20z = 25
 Sparse & Dense matrices
• In numerical analysis and scientific computing, a sparse matrix or
sparse array is a matrix in which most of the elements are zeros.
• There is no definition strict definition regarding the proportion of zero
value elements for a matrix to qualify as sparse, but a common
criterion is that the number of non zero elements is roughly equal to
the number of rows or columns
• By contrast, if most of the elements are non zero the matrix is
considered as dense
• The number of non zero valued elements divided by the total number
of elements is referred to as sparsity of the matrix
 Sparse & Dense matrices

Advantages of using a sparse matrix instead of a simple matrix

(i) Storage – There are lesser non-zero elements than zeros and thus
lesser memory used to store only those elements
(ii) Computing time – Can be saved by logically designing a data
structure, traversing only non-zero elements
 Sparse & Dense matrices
TYPES OF SPARSE MATRICES

There are different variations of sparse matrices which depends upon

the nature of the spacing of the matrices
Based on these properties sparse matrices can be
(i) Regular sparse matrices
(ii) Irregular/Non regular sparse matrices
 Sparse & Dense matrices
REGULAR SPARSE MATRICES

• A regular sparse matrix is the square matrix with a well-defined

sparsity pattern, i.e. non zero elements occur in a well-defined
pattern
• The various types of regular sparse matrices are
❑ Lower triangular regular sparse matrix
❑ Upper triangular regular sparse matrix
❑ Tridiagonal regular sparse matrix
 Sparse & Dense matrices
LOWER TRIANGULAR REGULAR SPARSE MATRICES

• A lower triangular sparse matrix is one where all elements above the
main diagonal are zeros
• Maximum number of non 0 elements of a lower triangular regular
sparse matrix of m rows =
 Sparse & Dense matrices
LOWER TRIANGULAR REGULAR SPARSE MATRICES

• It can be stored in a 1-D dimensional array row by row

Eg: B = { A11, A21, A22, A31, A32, A33, A41, A42, A43, A44, A51, A52, A53, A54, A55}
 Sparse & Dense matrices
UPPER TRIANGULAR REGULAR SPARSE MATRICES

• All elements below the main diagonal are zeros

• Maximum number of non 0 elements of an upper triangular regular
sparse matrix of m rows =
 Sparse & Dense matrices
UPPER TRIANGULAR REGULAR SPARSE MATRICES

• It can be also be stored in a 1-D dimensional array column by column

Eg: B = { A11, A21, A22, A13, A23, A33, A14, A24, A34, A44, A15, A25, A35, A45, A55}
 Sparse & Dense matrices
TRIDIAGONAL REGULAR SPARSE MATRIX

• All non zero elements lie on the main diagonal and on the diagonals
immediately above and below the main diagonal
TRIDIAGONAL REGULAR SPARSE MATRIX

• It can be stored in 1-D array as

B = {A11,A12,A21,A22,A23,A32,A33,A34,A43,A44,A45,A54,A55}

Total number of non zero elements in nxn square matrix is

n + (n-1) + (n-1) = 3n - 2
Irregular Sparse Matrices
• The irregular sparse mattresses are the ones that have an irregular or
unstructured pattern of occurrences of non zero elements

• Storage → By index method

→ Compressed row format
Solution of Tridiagonal systems
• Tridiagonal systems occur commonly in solving problems in many different areas of
numerical analysis

Eg: Spline function interpolation, finite difference solution of 2-point boundary value problems of 2nd order,
ordinary differential equations, reactor problems, membrane problems (reverse osmosis) etc..

• For Gauss Elimination method → No of computations α n3 (n- no of unknowns)

When n = 10 → 1000

n = 100 → 106

n=1000 → 109

• Thus, Gauss Elimination is computationally intensive when no of unknowns are large

Solution of Tridiagonal systems
• For a general purpose problem Ax = b, Gauss Elimination is one efficient method.

• However, a lot of problems of interest to engineers (chemical, mechanical etc) are not
general purpose, they have a specific structure (banded structure or sparse structure)

• Banded solvers can be used to solve them. Eg. Tridiagonal matrix algorithm (TDMA)

• Computation load can be reduced compared to a general Gauss Elimination method,

since we know in advance that some elements above and below diagonal are 0

• For system Ax= b, (where A is a tridiagonal matrix) TDMA method can be used to reduce
the number of computation from n3 to n1
Solution of Tridiagonal systems
• Tridiagonal structure is common in lot of engineering problem,

especially when we express a model for an engineering problem in

terms of a differential equation

• Consider a hot rod losing heat to surroundings. We want to find out

how temperature varies along the rod

• This results in an ODE-boundary value problem

• Now we discretize the rod into (n – 1) segments and we get n equations in

n unknowns
• On rearranging we get,

This is the tridiagonal system which occurs in a heat conduction problem

Thomas Algorithm (TDMA)
• The augmented matrix has a general form:

• The Thomas algorithm is a procedure for

solving tridiagonal system of equations

• The method of solution in the Thomas

algorithm is similar to the Gauss Elimination
method in which the system is first changed
to UPPER TRIANGULAR FORM and then
solved using back substitution
Thomas Algorithm (TDMA)
• The Thomas algorithm is much more efficient since only the non 0

elements of the matrix of coefficients are stored and only the necessary

operations are executed

• Unnecessary operations on the zero elements are eliminated

• Once the matrix of coefficients is in upper triangular form the values of

the unknowns are calculated by back substitution

THOMAS ALGORITHM
STEP 1
Row 1 is the pivot row
• Divide R1 by a11……………R1→ R1/a11
The augmented matrix becomes
THOMAS ALGORITHM
• Now use pivot row R1 to get a21 = 0
i.e. R2 → R2 – a21R1 a22’ = a22 – a21a12’
b2’ = b2 – a21b1’
THOMAS ALGORITHM
STEP 2
Row 2 is the pivot row
• Divide R2 by a22……………R2→ R2/a’22
The augmented matrix becomes
THOMAS ALGORITHM
STEP 2
Use pivot row R2 to make a32 = 0 ---- R3→R3-a32R2
The augmented matrix becomes
THOMAS ALGORITHM
• Continue the steps for (n-1) rows
• Total computational requirement = (n-1)2 + (n-1) 4

• Thus, instead of n3 computations, only n1 computations are required

in Thomas algorithm for tridiagonal matrices
• Eg: from n=10 to n=100
THOMAS ALGORITHM
• The sequence of operation finally leads to

• Back substitution gives the solution vector

• Solve by Thomas algorithm, the following system of equations
2x1 + x2 = 1
2x1 + 3x2 + x3 = 2
x2 + 4x3 + 2x4 = 3
x3 + 3x4 = 4