Linear System of Equations

a11 x1  a12 x 2  a13 x 3  b1

a21 x1  a22 x 2  a23 x 3  b2
a31 x1  a32 x 2  a33 x 3  b3

 a11 a12 a13   x1   b1 

a    
 21 a 22 a 23  x 2   b 2  Ax=b
a 31 a 32 a 33   x 3  b 3 

n*n Column Column

Matrix vector vector

Gauss Elimination method

• Forward elimination
– Starting with the first row, add or
subtract multiples of that row to
eliminate the first coefficient from
the second row and beyond.
– Continue this process with the
second row to remove the second
coefficient from the third row and
beyond.
– Stop when an upper triangular
matrix remains.
• Back substitution
– Starting with the last row, solve for
the unknown, then substitute that
value into the next highest row.
– Because of the upper-triangular
nature of the matrix, each row will
contain only one more unknown.

1
 Gauss Elimination method
Number of flops for Gauss Elimination method:

n Total flops Forward Back 2n3/3 % due to

Elimination substitution Elimination
10 805 705 100 667 87.58%

100 681550 671550 104 666667 98.53%

1000 6.68*108 6.67*108 106 6.67*108 99.85%

• As the system gets larger, the cost of computing increases

- Flops increases nearly 3 orders of magnitude for every order of
increase in the number of equations
• Most of the efforts is incurred in the forward elimination steps. 3

Gauss Elimination method

clc; close all; format long e;
% A= coefficient matrix % b=right hand side vector
% x=solution vector
A=[10 -7 0; -3 2 6; 5 -1 5]; b=[7;4;6];
[m,n]=size(A);
if m~=n
error('Matrix A must be square');
end
Nb=n+1;
Aug=[A b]; % Augmented matrix, combines A and b
%----forward elimilation-------
for k=1:n-1 % from row 1 to n-1
if abs(A(k,k))<eps
error('Zero Pivot encountered')
end
for j=k+1:n
multiplier=Aug(j,k)/Aug(k,k);
Aug(j,k:Nb)=Aug(j,k:Nb)-multiplier*Aug(k,k:Nb);
end
end 4

2
 Gauss Elimination method

%----Back substitution---------
x=zeros(n,1);
x(n)=Aug(n,Nb)/Aug(n,n);

for k=n-1:-1:1
for j=k+1:n
Aug(k,Nb)=Aug(k,Nb)-Aug(k,j)*x(j);
% computes b(k)-sum(akj*xj)
end

x(k)=Aug(k,Nb)/Aug(k,k);
end

Cost of computing
Example 1:
Estimate the time required to carry out back substitution on a system of
500 equations in 500 unknowns, on a computer where elimination takes
1 second.

For Elimination
approx. flops, fe= 2n3/3=2*5003/3
time needed, te=1 sec
For Back substitution
approx. flops, fb=n2=5002
time needed, tb=?

tb/te=fb/fe
tb=te*fb/fe=3/(2*500)=0.003 sec

Back substitution time=0.003 sec

Total computation time=1.003 sec

3
Cost of computing
Example 2:
Assume that your computer can solve a 200*200 linear system Ax=b in
1 second by Gauss elimination method. Estimate the time required to
solve four systems of 500 equations in 500 unknowns with the same
coefficient matrix, using LU factorization method.

For Gauss Elimination

approx. flops, fg= 2n3/3=2*2003/3
time needed, tg=1 sec
For LU
approx. flops, fl=2n3/3+2kn2=2*5003/3+2*4*5002
time needed,
tl=tg*fl/fg=(2*500^3/3+2*4*5002)/(2*2003/3)=16 sec

If we do the same 4 problems by Gauss Elimination:

tg=4*5003/2003=62.5 sec

Time save=46.5 sec  300% !!!

7

Forward and backward error

Definition
Let xc be an approximate solution of the linear system Ax=b.

The residual is the vector r=b-Axc

The backward error = ||b-Axc||

The forward error = ||xexact-xc||

x exact  x c 

x exact 
Error magnification factor=
b  Ax c 

b 
8

4
Condition number

Definition
The condition number of a square matrix A, cond(A), is the maximum
possible error magnification factor for Ax=b, over all right-hand side b.

Theorem
The condition number of a n*n matrix A is cond(A)=||A||*||A-1||

For ill-conditioned system, condition number is high

If the coefficient matrix A have t-digit precision, and cond(A)=10k

The solution vector x is accurate upto (t-k) digits.

Pivoting

• Problems arise with Gauss elimination if a coefficient along the

diagonal (Pivot element) is 0 (problem: division by 0) or close to 0
(problem: round-off error)

• One way to combat these issues is to determine the coefficient with

the largest absolute value in the column below the pivot element.
The rows can then be switched so that the largest element is the
pivot element. This is called partial pivoting.

• If the rows to the right of the pivot element are also checked and
columns switched, this is called complete pivoting.

10

5
 Gauss Elimination with partial pivoting
clc; close all; format long e;
% A= coefficient matrix % b=right hand side vector
% x=solution vector
A=[10 -7 0; -3 2 6; 5 -1 5]; b=[7;4;6];
[m,n]=size(A);
if m~=n, error('Matrix A must be square');end
Nb=n+1;
Aug=[A b]; % Augmented matrix, combines A and b
%----forward elimilation-------
for k=1:n-1 % from row 1 to n-1
%-----------partial pivoting------------
[big, big_pos]=max(abs(Aug(k:n,k)));
pivot_pos=big_pos+k-1;
if pivot_pos~=k
Aug([k pivot_pos],:)=Aug([pivot_pos k],:);
end
%--------------------------------------
for j=k+1:n
multiplier=Aug(j,k)/Aug(k,k);
Aug(j,k:Nb)=Aug(j,k:Nb)-multiplier*Aug(k,k:Nb);
end 11
end

Gauss Elimination method

%----Back substitution---------
x=zeros(n,1);
x(n)=Aug(n,Nb)/Aug(n,n);

for k=n-1:-1:1
for j=k+1:n
Aug(k,Nb)=Aug(k,Nb)-Aug(k,j)*x(j);
% computes b(k)-sum(akj*xj)
end

x(k)=Aug(k,Nb)/Aug(k,k);
end

12

6
Scaling

Scaling is the operation of adjusting the coefficients of a set of

equations so that they are all of the same order of magnitude.

• Whenever the coefficients in one column are widely different from

those in another column, scaling is beneficial.

13

Tridiagonal systems
function [xout]=tridiagonal(a,b,c,d)
N=length(b);
%---forward elimination------------
for k=2:N
multiplier=a(k)/b(k-1);
b(k)=b(k)-multiplier*c(k-1);
d(k)=d(k)-multiplier*d(k-1);
end
%-----back substitution----------
x(N)=d(N)/b(N);
for k=N-1:-1:1
x(k)=(d(k)-c(k)*x(k+1))/b(k);
end
xout=x.';
%---------------------------------
Main:
format long e; clc; clear all; close all;
a=[0 -1 -1 -1]; b=[2.04 2.04 2.04 2.04];
c=[-1 -1 -1 0]; d=[40.8 0.8 0.8 200.8];
[xout]=tridiagonal(a,b,c,d)
14

7
Jacobi method
function [xout iter_number]=Jacobi(A,b,max_iter,TOL)
N=length(b);
d=diag(A); % takes the main diagonal elements only
x=zeros(N,1); % initial guess vector
for k=1:max_iter
x=(b-(A-diag(d))*x)./d;
maxresidual=norm(b-A*x,inf);
if maxresidual<TOL
break;
end
End
iter_number=k;
xout=x;
%----------------------------------------------------
Main:
format long e; clc; clear all; close all;
A=[3 -0.1 -0.2; 0.1 7 -0.3; 0.3 -0.2 10];
b=[7.85;-19.3; 71.4];
max_iter=1000;TOL=1e-5;
[xc iter_number]=Jacobi(A,b,max_iter,TOL); 15

Gauss Seidel and Jacobi method: Graphical depiction

16

8
 Gauss Seidel method
function [xout iter_number]=GaussSeidel(A,b,max_iter,TOL)
N=length(b);
d=diag(diag(A)); % takes the main diagonal elements only
x=zeros(N,1); % initial guess vector
U=triu(A,1); % above the main diagonal
L=tril(A,-1); % below the main diagonal
for k=1:max_iter
bL=b-U*x;
for j=1:N
x(j)=(bL(j)-L(j,:)*x)./d(j,j);
end
maxresidual=norm(b-A*x,inf);
if maxresidual<TOL
break;
end
End
iter_number=k;
xout=x;
%---------------------------------------------------------
17

SOR method
function [xout iter_number]=SOR(A,b,relax,max_iter,TOL)
N=length(b); d=diag(diag(A));
x=zeros(N,1); % initial guess vector
U=triu(A,1); % above the main diagonal
L=tril(A,-1); % below the main diagonal
for k=1:max_iter
bL=relax*(b-U*x)+(1-relax)*d*x;
for j=1:N
x(j)=(bL(j)-relax*L(j,:)*x)./d(j,j);
end
maxresidual=norm(b-A*x,inf);
if maxresidual<TOL
break;
end
end
iter_number=k;
xout=x;

18

9
 Sparse matrix computations
 3 1 0 0 0 0 0 0 0 0 0 2 
1

 1 3  1 0 0 0 0 0 0 0 1
0 
 2
 0 1 3 1 0 0 0 0 0 1
0 0
 
2

 0 0  1 3  1 0 0 0 1
2 0 0 0
 0 0 0 1 3 1 0 1
0 0 0 0
 2

 0 0 0 0 1 3 1 0 0 0 0 0 
AS  
0 0 0 0 0 1 3 1 0 0 0 0 
 
0 0 0 0 1
2 0 1 3 1 0 0 0 
0 0 0 1
0 0 0 1 3 1 0 0 
 2 
0 0 1
2 0 0 0 0 0 1 3 1 0 
 
0
1
2 0 0 0 0 0 0 0  1 3  1
 12 0 0 0 0 0 0 0 0 0  1 3 

• No row of A has more than 4 non-zero entries.

• Since fewer than 4n of the n2 potential entries are non zero.
We may call this matrix Sprase. 19

Sparse matrix computations

We want to solve a system of equations defined by AS for n=105
What are the options?

Treating the coefficient matrix A as a full matrix means

 storing n2=1010 elements each as a double precision floating point number.

One double precision floating point number requires 8 bytes (64 bit) of storage.

So, n2=1010 elements will require 8*1010 bytes

= 80 gigabytes (approx.) of RAM
Not only the size is enemy, but so is time.

The number of flops by Gauss Elimination will be order of n3  1015

If a PC runs on the order of a few GHz (109 cycle per second), an upper
bound of floating point operations per second is around 108.
Therefore, time required of Gauss Elimination=1015/108=107 seconds.
Note: 1 year =3*107 seconds.
20

10
 Sparse matrix computations
• So, it is clear that Gauss Elimination method for this problem is not an
overnight computation.

On the other hand,

• One step of an iterative method will require 2*4n=8*105 (approx.) operations.

• If the solution needs 100 iterations in Jacobi method, total operations=108

• Time required in modern PC for Jacobi method will be

=1 second (approx.) or less !!!

21

Rank Deficiency
Kronecker-Capelli theorem:

A system has a SOLUTION if and only if its coeffiecient matrix (denoted by A)

and its augmented matrix (denoted by Ab) have the same Rank.

Let, Rank of A matrix=R

Number of equation=n

If R=n,  the system has unique solution (full Rank system)

If R<n,  the system has an infinite number of solution (Rank-deficient system)

• We can express the system in terms of (n-R) variables

• under-determined system (number of equation<number of unknowns)

22

11
Rank Deficiency
Example 1: Solution:

A=[ 1 -2 3 ; RA=2
2 4 -1 ; RAb=3
-1 -14 11 ];
Warning: Matrix is singular to
b=[ 5; working precision.
7;
xc =
2];
NaN
Ab=[A b]; % Augmented matrix
-Inf
RA=rank(A);
RAb=rank(Ab); -Inf
xc=A\b

23

Rank Deficiency
Example 2: Solution:
RA=2
A=[ 1 -2 3 ;
2 4 -1 ; RAb=2
-1 -14 11 ];
Warning: Matrix is singular to
working precision.
b=[ 5;
xc1 =
7;
NaN
1];
NaN
NaN
Ab=[A b]; % Augmented matrix
RA=rank(A);
xc2 =
RAb=rank(Ab);
2.131455399061033e+000
1.107981220657277e+000
xc1=A\b
1.694835680751174e+000
xc2=pinv(A)*b

24

12