Chapters: 3, 4 and 6

Solving Linear Systems

A system of n simultaneous, non-homogenous, algebraic, linear equations is given in
the form:

Ax b
A : is an n n coefficient matrix.
b : is an n 1 right-hand side vector.
x : is an n 1 vector of unknowns.
A11
A
21

An1

A12 A1n x1 b1
A22 A2 n x2 b2
=


An 2 Ann xn bn

Ann

xn1

bn1

First: Exact Methods

1)

Cramer's Rule:

This method is one of the least efficient for solving a large number of nonhomogeneous linear equations.
In general, given:

A x b
ij

xi

detA j

i j 1,2,3,, n

Aj

det A
A
where,
A : is the determinant of the coefficient matrix A .
A j : is the determinant of the coefficient matrix A with its j th column replaced by

the bi vector.

Example:
Use Cramer's Rule to solve the following system of equations:
2x1 3x2 5
x1 x2 5
Solution:

Ax b
2 3 x1 5
1 1 x 5

2
5 3
A
5 1
x1 1
4
2 3
A
1 1
2 5
A
1 5
x2 2
1
2 3
A
1 1

x 4
1
x2 1
Elementary Row Operations:
Interchange rows k and j .
(2) Multiply row k by 0 .
(3) Add times row j to row k .
(1)

2)

Gauss-Jordan Elimination Method:

We can solve Ax b using Gauss-Jordan elimination method.

Algorithm:
1. Form the augmented matrix: C A, bn( n1)
2. For j 1 to n do
{
a) Compute the pivot index "p": n p j such that:
n

C pj max Cij
i j

b) If C pj 0 , exit (no solution).

c) If p j , interchange rows p and j .
d) Divide row j by the pivot C jj .

e) For each i j , subtract Cij times row j from row i . Ri = Ri Cij R j

3. Partition C as C I , x .

Example:
Solve the following system of equations using Gauss-Jordan Elimination method:

1 1 0 x1 2
2 2 1 x 1

2
0
1 2 x3 6
Solution:
2
1 1 0
C 2 2 1 1

0
1 2 6 34

j 1

3 p 1

p max Cij 2
i j

C p1 2 2 0
p 2 j 1

interchange rows 2 and 1

2 2 1 1
1 1 0
2

0
1 2 6

R1
2

1 1 0.5 0.5
1 1 0
2

0 1 2 6

1 1 0.5 0.5
0 0 0.5 1.5

0 1
2
6

row2 - C21.row1
row3 - C31.row1

j2
32

3 p 2

j 3

R2
1

1 1 0.5 0.5
0 1
2
6

0 0 0.5 1.5

1 0 1.5 6.5
0 1 2
6

0 0 0.5 1.5

row1 - C12.row2
row3 - C32.row2

p3
C32 1

interchange rows 3 and 2

1 1 0.5 0.5
0 1
2
6

0 0 0.5 1.5

R3
0 .5

3 p 3

p3

C33 0.5 0
1 0 1.5 6.5
0 1 2
6

0 0
1
3

1 0 1.5 6.5
0 1 2
6

0 0
1
3

R1 C13 R3
R2 C23 R3

1 0 0 2
0 1 0 0

0 0 1 3

I
2 x1
x 0 x2

3 x3

3)

Matrix Inverse Method:

Ax b

x A 1b

1
- First, find Anxn

- Then, x A b
1

Algorithm:
1. Form the augmented matrix: C A, I n2n
2. For j 1 to n do
{
a) Compute the pivot index n p j such that:
n

C pj max Cij
i j

b) If C pj 0 , exit (no inverse).

c) If p j , interchange rows p and j .
d) Divide row j by the pivot C jj .
e) For each i j , subtract Cij times row j from row i .

3. Partition C as C I , A1 .

Example:
Consider the following 3 3 matrix:
1 1 0
A 2 0 4 Find A 1

0 2 1
Solution:
1 1 0 1 0 0
C A, I 2 0 4 0 1 0

0 2 1 0 0 1 36
1 0 0 0.8 0.1 0.4
0 1 0 0.2 0.1 0.4 I , A1

0 0 1 0.4 0.2 0.2

0.8 0.1 0.4
1
A 0.2 0.1 0.4

0.4 0.2 0.2

5

4)

Gaussian Elimination Method:

Ax b
C A, bn( n1)
D, en( n1)
D : is an upper-triangular matrix (matrix with zeros below the main diagonal)
Algorithm:
1. Form the augmented matrix: C A, bn( n1) .
2. For j 1 to n do
{
a) Compute the pivot index n p j such that:
n

C pj max Cij
i j

b) If C pj 0 , exit (no solution).

c) If p j , interchange rows p and j .

C ij
d) For each i j , subtract
C
jj

times row j from row i .

c
Ri Ri ij R j
c jj

3. Partition C as C D, e , where D is an n n and e is n 1 .

4. For j n down to 1 , compute:

xj

n
1

j D ji xi
D jj
i j 1

Example:
Solve the following linear algebraic system using Gaussian Elimination method:
1 0 2 x1 1
2 1 3 x 1

2
4 1 8 x3 2

Solution:

1 0 2 1
C 2 1 3 1

4 1 8 2
j 1:
C31 4

p3
c21
R2
c
R1
11

4 1 8 2
2 1 3 1

1 0 2 1

c31
R3
c
R1
11

1
8
2
4
0 1.5 1 2

0 0.25 0 0.5

j2
j 3
1
8
2
4
C 0 1.5
1
2

0
0
0.1667 0.8333

j 3

3
1
0
.
8333

D ji xi 5

0.1667

i 4
3
1
1

2 D23 x3 1 2 (1)5 3 2
x2
2 D ji xi
1.5
1.5
1.5
1.5
i 3
3
1
1
1
x1 2 D ji xi 2 D12 x2 D13 x3 2 (1 2) (8 5) 9
4
4
4
i 2

x3

x1 9
x 2
2
x3 5

LU Decomposition (Triangular Factorization)

Definition:
The nonsingular matrix A has a triangular factorization if it can be expressed as the
product of a lower triangular matrix L and an upper triangular matrix U.
A = LU
L: it has 1s along the diagonal.
U: it has nonzero diagonal elements.
1
m
L = 21
m31

m n1

0
1

0
0

m32 1
m n2

0
0
0

u11

0
U=
0

u1n

u nn

u12
u 22
0
0

e.g
Given,
4 3 1
A = 2 4 5
1
2 6

Convert the matrix A to the form A=LU.

Solution:
1 0 0
4 3 1
2 4 5 = m 1 0
21

1
2 6
m31 m32 1

u11
0

u13
u 23
u 33

u12
u 22
0

4 = u 11 , 3 = u12 , -1 = u13
m21 u11 = -2 m 21 = - 0.5
m21u12 + u22 = - 4 m21 = - 4 (- 0.5*3) = -2.5
m21u13 + u23 = 5 u23 = 5 (- 0.5*-1) = 4.5
m31u11 = 1 m31 = 0.25
1
2 *3
4
m31u12 +m32u22 = 2 m32 =
= -0.5
2.5

m31u13 + m32u23 + u33 = 6

u33 = 6

1
1
* 1 ( * 4.5) 8.5
4
2

4 3 1
2 4
5

1
2
6

1
= 0.5
0.25

0
1
0.5

0
0
1

4
0

3
2.5
0

1
4.5
8.5

U
8

5)

Solution of a Linear System: AX = b

Suppose that the coefficient matrix A of the linear system: AX = b has a triangular
factorization A = LU, then the solution to: AX = b LUX = b can be obtained by
defining:
Y = UX and then solving two systems:
First, solve LY = b for Y using forward substitution.
Then, solve UX = Y for X using back substitution.
e.g.
Solve the following system of equations using the triangular factorization method.
1
2

10
12

8
10

4
8

6
1

x1 21
x
2 52
x3 79

x4 82

Solution:
1. A = L U
1
2

10
12

8
10

4
21

m31
8

6
m41
1

0
1

0
0

m32

m42

m43

0
0
0

u11 u12
0
u 22

0
0

0
0

Find the values of the ms and us. Therefore,

1
2
L=
3

0
0
0

1
0
,U=
0

1
2
3

2. solve LY = b for Y using forward substitution:

1
2

0
0
0

y1
21
y

2 = 52
y3
79


82
y4

u13
u 23
u33
0

u14
u 24
u 34

u 44

y1 = 21
2y1 + y2 = 52 y2 = 10
3y1 + y2 + y3 = 79 y3 = 6
4y1 + y2 + 2y3 + y4 = 82 y4 = -24

Y=

21
10

24

3. solve UX = Y for X using back substitution:

1
0

1
2
3

x1
x
2 =
x3

x4

21
10

24

x1 + 2x2 + 4x3 + x4 = 21 x1 = 1
4x2 2x3 +2x4 = 10
x2 = 2
-2x3 +3x4 = 6
x3 = 3
- 6x4 = -24
x4 = 4
1
2
X=
3

4

is the solution of the system of the equations.

e.g.
Show that the following matrix can not be factored directly as: A = LU.
1
A = 4
2

2
8
3

6
1
5

10

ILL-Conditional Systems:
Consider the following system:
100 200 x1 100
200 401 x 100

x1 201, x2 100

Now, consider the system:

101 200 x1 100
200 400 x 100

x1 50 , x2 24.75

Clearly the solution is very sensitive to the values of A in this case. When the
solution is highly sensitive to the values of the coefficient matrix A or the righthand side vector b , the equations are said to be ill-conditioned.

Vector and Matrix Norms:

A norm of a vector is a generalization of the absolute value of a scalar. That is, a
norm is a measure of the length or magnitude of the vector.
There are several ways to define a norm of an n 1 vector x .
1) l norm

x max xi

i 1

x : is called the infinity norm of x .

2) l 2 norm

x2

xi2
i 1

l 2 norm: is called the Euclidean norm (distance from the origin).

In order to develop a measure of the sensitivity of a system, we need to consider the
notion of a matrix norm:
n
n

A max Aij
i 1 j 1

11

Example:

4 6 9 x1 1
5 8 2. x 12

2
7 3 1 x3 13
3

A max Aij max{19,15,11} 19

i 1
j 1
3

The solution of this system is:

x1 2
x 0
2
x3 1
X

max{2,0,1} 2

4 1 5

Properties of norms:
1. x 0
2.
3.
4.
5.

for x 0 .
x 0
for x 0 .
: is a scalar.
x . x
x y x y .
Ax A . x .

Condition Number:
AX = b
Consider the effect of changes in the right-hand side vector b, say from b to b+b.
Let X denotes the solution to the original system, and let X + X denotes the
solution to the perturbed system. Then, A(X + X) = b+b which reduces to A X
= b because AX = b
X = A-1 b
12

Using:
AX A . X
X A 1 . b
Using the properties of norms, the relative change in the solution can be expressed in
terms of the perturbation as follows:
X
X

K ( A)

b
b

K(A): is a scalar called the condition number of the matrix A and is defined as follows,
assuming A is nonsingular.
K(A) = A . A1
Note that K(A) is a measure of the relative sensitivity of the solution to changes in b .
Properties of the condition number:
1. When K(A) becomes large, the system is regarded as being an ill-conditioned.
e.g
100

A=
200

200
401

determine K(A).

Solution :
4.01

2
1
A max{300,601} 601

A-1 =
2

A 1 max{6.01,3} 6.01
K ( A) 601 * 6.01 3612.01 large value

Our system is ill-conditional system.

2. > K(A) 1
3. K(I) = 1
4. K(CA) = K(A)

C is a scalar.

5. If K(A) is large , then A is close to singular matrix.

13

e.g.
4.1

Given, b = , x=
9.7
0
0.34

4.11

if b is changed to bnew = and x is changed to xnew = , determine K(A).

9.7
0.97
X

K(A)

X
b
b

0.66

X =

0.97
X 0.97

b 0.01 ,

0.01

b =

X 1
b 9.7

0.97
K(A) 1
940.9
0.01/ 9.7

K(A) 940.9

Second: Iterative (Numerical) Methods

1)
2)
3)

Jacobi method.
Gauss-Seidel method.
Relaxation methods:
a. Successive Under-Relaxation.
b. Successive Over-Relaxation.

These are methods by which an approximation to the solution of a system of

linear algebraic equations may be obtained. The iterative procedure may not
always yield a solution.
A sufficient condition for a solution to be found is that the absolute value of the
diagonal coefficient in any equation must be greater than the sum of the absolute
values of all other coefficients appearing in that equation.
Akk Akj
: n k 1

j k

14

Example:
a11 x1 a12 x2 a13 x3 b1
a21 x1 a22 x2 a23 x3 b2
a31 x1 a32 x2 a33 x3 b3
Test:
a11 a12 a13

a22 a21 a23

a33 a31 a32
To start the iterative scheme, an initial guess x 0 must be chosen. For example, if
nothing is known about the solution, the initial guess x 0 0 is as good as any.
Typically, the iterations are repeated until the norm of the difference between
successive estimates is sufficiently small, that is:
x ( k 1) x ( k )

" : tolerance"

1) Jacobi Method:
Jacobi iteration can be written as a set of scalar equations expressed directly in
terms of the components of A and b ,
1
(k )
xi( k 1)
:n i 1
bi Aij x j
Aii
j i

Example:
Given the following system of equations:
5 0 2 x1 7
3 5
1 . x2 2

0 3 4 x3 4
Use Jacobi iterative method to find x1 , x2 and x3 .
Do two iterations, let the initial guess x 0 0 .
5 0 2

531
4 0 3
15

1
x1( k 1) (7 2 x3( k ) )
5
1
x2( k 1) (2 3x1( k ) x3( k ) )
5
1
x3( k 1) (4 3x2( k ) )
4
1st iteration (k=0)
(1)
1

0 7
1
(0)
(7 2 x3 ) 1.4 :
5
5

(1)
2

0
0 2
1
( 0)
( 0)
(2 3x1 x3 ) 0.4
5
5

x1(1) is the value of x1 in the first iteration.

1
4
x3(1) (4 3x2( 0) ) 1.0
4
4
x (1) x ( 0 ) ?

1.4 0 1.4
x x 0.4 0 0.4

1.0 0 1.0
x (1) x ( 0 ) 1.4
(1)

(0)

2nd iteration (k=1)

1
1
x1( 2) (7 2 x3(1) ) [7 (2 1.0)] 1.0
5
5
1
1
x2( 2) (2 3x1(1) x3(1) ) [2 (3 1.4) (1.0)] 0.24
5
5
1
1
x3( 2) (4 3x2(1) ) [4 (3 0.4)] 0.7
4
4
x ( 2 ) x (1) ?

16

 1.0 1.4 0.4

x ( 2 ) x (1) 0.24 0.4 0.64

0.7 1.0 0.3

x ( 2 ) x (1) 0.64

Jacobi Estimates:

x1

0
1
2
3
4
5
6
7
8
9

0.0000000
1.4000000
1.0000000
1.1200000
0.9280000
0.9820000
0.9892000
1.0156599
1.0048599
0.9995950

x2
0.0000000
0.4000000
- 0.2400000
- 0.0600000
- 0.0360000
0.0522000
0.0162000
- 0.0013500
- 0.0118259
- 0.0027135

x3
0.0000000
- 1.0000000
- 0.7000000
- 1.1799999
- 1.0450000
- 1.0270000
- 0.9608500
- 0.9878500
- 1.0010124
- 1.0088694

2) Gauss-Seidel Method:
The convergence rate of the Jacobi method can be improved using Gauss-Seidel
method. In general, it is better to take advantage of the most recent information
when updating an estimate of a solution as follows:
( k 1)
i

i 1
n
1
( k 1)
(k )

bi Aij x j Aij x j
Aii
j 1
j i 1

17

:n i 1

Example:
Use Gauss-Seidel iterative method to find x1 , x2 and x3 for the same system in the
previous example:
5 0 2 x1 7
3 5
1 . x2 2

0 3 4 x3 4
Do two iterations, and let the initial guess x 0 0 .
Solution:
1
x1( k 1) (7 2 x3( k ) )
5
1
x2( k 1) (2 3x1( k 1) x3( k ) )
5
1
x3( k 1) (4 3x2( k 1) )
4
1st iteration (k=0)
1
7
x1(1) (7 2 x3( 0 ) ) 1.4
5
5
1
1
x2(1) (2 3x1(1) x3( 0) ) [2 (3 1.4) 0] 0.44
5
5
1
1
x3(1) (4 3x2(1) ) [4 (3 0.44)] 1.33
4
4
x (1) x ( 0 ) ?

1.4 0 1.4
x (1) x ( 0 ) 0.44 0 0.44

1.33 0 1.33
x (1) x ( 0 ) 1.4
2nd iteration (k=1)
1
1
x1( 2) (7 2 x3(1) ) [7 (2 1.33)] 0.868
5
5
1
1
x2( 2) (2 3x1( 2) x3(1) ) [2 (3 0.868) (1.33)] 0.1452
5
5
1
1
x3( 2) (4 3x2( 2) ) [4 (3 0.1452)] 0.8911
4
4
18

x ( 2 ) x (1) ?

0.868 1.4 0.532

x ( 2 ) x (1) 0.1452 0.44 0.5852

0.8911 1.33 0.4389

x ( 2 ) x (1) 0.5852

Gauss-Seidel Estimates:
k

0
1
2
3
4
5
6
7
8
9

x1
0.0000000
1.4000000
0.8680000
1.0435599
0.9856252
1.0047437
0.9984345
1.0005165
0.9998295
1.0000563

x2
0.0000000
- 0.4400000
0.1452000
- 0.0479159
0.0158123
- 0.0052180
0.0017220
- 0.0005682
0.0001875
- 0.0000619

19

x3
0.0000000
- 1.3299999
- 0.8911000
- 1.0359370
- 0.9881408
- 1.0039135
- 0.9987085
- 1.0004262
- 0.9998593
- 1.0000464

Jacobi vs Gauss-Seidel Estimates:

x1
k

0
1
2
3
4
5
6
7
8
9

x2

x3

Jacobi

Gauss-Seidel

Jacobi

Gauss-Seidel

Jacobi

Gauss-Seidel

0.0000000

0.0000000

0.0000000

0.0000000

0.0000000

0.0000000

1.4000000

1.4000000

0.4000000

- 0.4400000

- 1.0000000

- 1.3299999

1.0000000

0.8680000

- 0.2400000

0.1452000

- 0.7000000

- 0.8911000

1.1200000

1.0435599

- 0.0600000 - 0.0479159

- 1.1799999

- 1.0359370

0.9280000

0.9856252

- 0.0360000

0.0158123

- 1.0450000

- 0.9881408

0.9820000

1.0047437

0.0522000

- 0.0052180

- 1.0270000

- 1.0039135

0.9892000

0.9984345

0.0162000

0.0017220

- 0.9608500

- 0.9987085

1.0156599

1.0005165

- 0.0013500 - 0.0005682

- 0.9878500

- 1.0004262

1.0048599

0.9998295

- 0.0118259

0.0001875

- 1.0010124

- 0.9998593

0.9995950

1.0000563

- 0.0027135 - 0.0000619

- 1.0088694

- 1.0000464

3) Successive Relaxation (SR) Method:

The convergence rate can sometimes be increased by introducing a relaxation factor
0,
i 1
n

( k 1)
xi( k 1) (1 ) xi( k )
b

A
x

:n i 1
i ij j
Aij x (jk )
Aii
j 1
j i 1

If: 1 , SR method reduces to the Gauss-Seidel method.

If: 1 0
Successive under relaxation (SUR).
If: 2 1 Successive over relaxation (SOR).
Example:
Use Successive Relaxation method to find x1 , x2 and x3 in the previous example
5 0 2 x1 7
3 5
1 . x2 2

0 3 4 x3 4
Do two iterations, let the initial guess x 0 0 and 1.1.
20

Solution:
3
1.1
(7 0 Aij x (jk ) )
5
j 2
(k )
0.22(7 2x3 )
1.54 0.44x3( k )

x1( k 1) 0.1x1( k )
0.1x1( k )
0.1x1( k )

1
3
1.1
(2 A2 j x (jk 1) A2 j x (jk ) )
5
j 1
j 3
(k )
( k 1)
(k )
0.1x2 0.22(2 3x1 x3 )
0.1x2( k ) 0.44 0.66x1( k 1) 0.22x3( k )

x2( k 1) 0.1x2( k )

( k 1)
3

2
1.1
0.1x
(4 A3 j x (jk 1) 0)
4
j 1
(k )
0.1x3 0.275(4 0 3x2( k 1) )
0.1x3( k ) 1.1 0.825x2( k 1)
(k )
3

x1( k 1) 0.1x1( k ) 0.44x3( k ) 1.54

x2( k 1) 0.1x2( k ) 0.66x1( k 1) 0.22x3( k ) 0.44
x3( k 1) 0.1x3( k ) 0.825x2( k 1) 1.1

1st iteration (k=0)

x1(1) 0 0 1.54 1.54

x2(1) 0 0.66x1(1) 0 0.44 (0.66 1.54) 0.44 0.5764
x3(1) 0 0.825x2(1) 1.1 (0.825 0.5764) 1.1 1.5755
1.54
x (1) x ( 0 ) 0.5764

1.5755
x (1) x ( 0 ) 1.5755

21

2nd iteration (k=1)

x1( 2) 0.1x1(1) 0.44x3(1) 1.54 0.6928
x2( 2) 0.1x2(1) 0.66x1( 2) 0.22x3(1) 0.44 0.387
x3( 2) 0.1x3(1) 0.825x2( 2) 1.1 0.6232

0.6928 1.54 0.8472

x ( 2 ) x (1) 0.387 0.5764 0.9634

0.6232 1.5755 0.9523

x ( 2 ) x (1) 0.9634

Solving a system of m-equations in n-unknowns with m > n

(more equations than unknowns) in the least square sence
(Normal equations)
Solve: Amn Xn1 = bm1
Such systems can be solved using Normal equations as follows:
AT AX = AT b
Note:
A Tnm Amn results in a square matrix of size nn
A Tnm bm1 results in a vector of size n1
Now, ATAX = AT b
Can be solved using the previous exact methods (Gauss-Jordan
elimination, inverse method, Gaussian method, etc)

22

e.g.
Find the least squares solution to the inconsistent system (more
equations than unknowns):
x1 + 4x2 = -2
x1 + 2x2 = 6
2x1 + 3x2 = 1
Solution:
1
1

4
2
3

ATA =

6
12

A=

6
12

12
29

2
, b 6
1

12
29

6
AT b
7

x1
6
x 7

2

Solving the system above,

X 1 3
X
X 2 1

least squares solution

This solution is NOT an exact solution but rather an approximate

that minimizes the distance between AX and b (i.e. it minimizes
the length of the error vector: e = AX b).

23