Numerical solution of systems of linear algebraic equations

1
Numerical solution of systems of linear algebraic
equations
In many problems of engineering interest, you
will need to be able to solve a set of linear
algebraic equations. The methods can be broadly
classified as
Solution of linear algebraic equations
Direct methods
Iterative methods

• Gauss-Elimination
• Gauss-Jordan
• LU Decomposition

• Gauss-Seidel
• Jacobi

2
Consider the following set of equations
where all the as and cs are constants. We
need to find all the xs such that all the above
equations are satisfied. Sometimes, I will write
the above equations using matrix notation
(2.1)
3
where
4
We can solve Eq. (2.1) by just inverting the
matrix A and then multiply C by A-1.
But this might be difficult and time consuming,
sometimes impossible to do if A is large. Can we
numerically try to solve it ?
5
Gauss-Elimination If you want to use the Gauss
elimination procedure to solve the following
system of linear equations,
you have to
6

• first augment A with C. We have to next do
  forward elimination
• to get the matrix into a form that is easily
  solvable.

7

• Eliminate x1 from the second and subsequent
  equations. For example,
• to eliminate x1 from the i th row, multiply row
  1 by ai1(1)/ a11(1) and
• subtract from row i. Hence,

for i1 and j1,2,3,n1
for i2,3,,n and j1
for i2,3,,n and j2,3,.n1
8

• Next, eliminate x2 from the third and subsequent
  equations. For example,
• to eliminate x2 from the i th row, multiply row
  2 by ai2(2)/ a22(2) and
• subtract from row i. Hence,

for i1, 2 and ji, i1, i2,,n1
for i3, 4,,n and j1,2
for i3,4,n and j3,4,.n1
9

• If you keep repeating the steps above, you will
  end up with a matrix
• that look something like

(2.2)
where
for i1,2,,k-1 and ji,i1,i2,n1
for ik,k1,,n and j1,2,k-1
for ik,k1,k2,,n and jk,k1,.n1
10
We can solve the system of equations given by Eq.
(2.2) by using back substitution. First look at
the last row of Eq. (2.2)
(2.3)
So we can get a numerical value for xn. Next,
look at the 2nd last row of (2.2)
(2.4)
Since we already know xn from Eq. (2.3), we can
also get a numerical value for xn-1 from Eq.
(2.4). Similarly, looking at the 3rd last row of
(2.2)
11
In general, the formula for getting the numerical
value for xi is
and
for in-1, n-2, n-3,.., 2,1
12

• The pitfalls of the Gauss elimination method for
  solving
• a linear system of equations are
• It is very expensive, compared to other methods
  that are
• currently available.
• The formula for the forward elimination stage
  requires that
• we divide by the diagonal components of the
  matrix A.
• So, we have to always rearrange our equations
  so to ensure
• that we are not dividing by zero. For example
  the following
• system

For this case, a110, so you will run into
problems in the first forward elimination stage.
If you re-write the above equation as
13
will eliminate the problem. The process of
interchanging rows so that you are not dividing
by zero is called partial pivoting. In general,
if the diagonal elements are very small, you
would like to pivot the rows so that you are
not dividing by a small number.
14
Example Solve the following system of 4
equations using Gauss Elimination
15
LU Decomposition Suppose we want to solve the
following linear equation
(2.5)
If we can somehow re-express A such at
(2.6)
where L and U are lower- and upper-triangular
matrices. The structure of L and U are
shown below
16
If we now substitute Eq. (2.6) into Eq. (2.5),
then we will have
(2.7)
Let
(2.8)
17
and substitute Eq. (2.8) into Eq. (2.7) to give
(2.9)
Thus we have replaced the original problem, Eq.
(2.5), by two problems given by Eq. (2.8) and Eq.
(2.9). However, Eqs. (2.8) and (2.9) are
much easier to solve because L and U are
lower- and upper-triangular matrices respectively.
This means that to solve Eq. (2.8) and (2.9),
all you have to do is do backward and forward
substitution respectively.

• A strategy of solving Eq. (2.5) using Eqs. (2.8)
  and (2.9) is the following
• Decompose A into L and U
• Solve Eq. (2.9) and obtain R by forward
  subtitution.
• Use the R obtain in the step above and
  substitute into Eq. (2.8)
• so that you can solve for X by backward
  substitution.

18
(1) decomposition
(2) forward substitution
(3) backward substitution
19
The problem now is how do we find L and R
from A ? To illustrate the procedure of
finding L and R, we will look at how it can
be done for a small 4x4 A matrix.

If we multiply the rows of L by the first
column of U and compare it with A, we get
(2.10)
So the first column of L is the first column of
A
Multiply the first row of L by the columns of
U to get
Hence, the first row of U is determined
(2.11)
20
We keep alternating between getting a column of
L and a row of U. So we next get the equation
for the second column of L by multiplying the
rows of L by the second column of U
This gives the 2nd column of L to be
(2.12)
Since we already know u12, l21, l31 and l41 from
Eqs. (2.10) and (2.11)
21
Multiply the second row of L by the columns of
U to get
This allows us to obtain the 2nd row of U to be
Since Eqs. (2.10), (2.11) and (2.12) will give
you l21 u13, u14 and l22
If you keep going you can get
22
In general, the formula for getting the elements
of L and U from the elements an n x n A
matrix can be written as
(2.13)
23
The LU decomposition algorithm outlined in Eq.
(2.13) is also known as the Crouts LU
decomposition algorithm.
After obtaining the L and the U matrix, we
need to solve for R
(2.14)
by forward substitution. From the first row,
(2.15)
24
From the 2nd row,
(2.16)
we can get a numerical value for r2 because we
know r1 (from Eq. (2.15)). From the third and
fourth rows of Eq. (2.14) we can get
In general, the formula for ri for a n x n system
is
(2.17)
25
Lastly, we need to solve
Use the method of back substitution which we was
introduced in the method of Gauss elimination.
Start with the last (fourth row)
(2.18)
26
Next, use the 2nd last (third row) to give
(2.19)
You can get a value for x3 because we know x4
from Eq. (2.18). The second and first row will
give
Which can be used to solve for x2 and x1 (in that
sequence). In general, the formula for a n x n
system is
(2.20)
27
Example Solve the following system of 4
equations using LU Decomposition
28
Iterative method Gauss-Seidel
This method is the most commonly used iterative
method. For conciseness, I will illustrate the
method with the 3x3 system below
(2.21)
(2.22)
(2.23)
Eqs. (2.21), (2.22) and (2.23) can be use to
solve for x1, x2, and x3 respectively,
(2.24)
(2.25)
(2.26)
29
We can now start the iterative process. WARNING,
you have to be careful that none of your aii s
are zero !!
1) Assume initial guess values for x1, x2, and
x3 x1(0), x2(0)and x3 (0). 2) Substitute the
guess values for x2 (0), and x3 (0) into Eq.
(2.24) to obtain a better estimate for x1.
Call this value x1(1). 3) Substitute x1(1) and
x3 (0) into Eq. (2.25) to obtain a new estimate
for x2 x2 (1). 4) Use x1(1) and x2 (1) in
Eq. (2.26) to obtain a better estimate for x3 x3
(1).
30
First iteration
Second iteration
31
Bit of Matlab code to do Gauss Seidel iteration.
for i1n, x(i)C(i) for j1n,
if j i x(i)x(i)-A(i,j)x(j) end
end x(i)x(i)/A(i,i) end
One big advantage of the Gauss-Seidel method is
that it works for nonlinear systems too !
32
Example Use Gauss-seidel iterative method to
solve the following linear system of equations
Example Use Gauss-seidel iterative method to
solve the following nonlinear system of equations
33
Newtons method for nonlinear equations
(revisited)
Earlier in the course, we talked about finding
roots of two simultaneous nonlinear equations.
What if we want to find roots of 3 or more
simultaneous nonlinear equations ? The methods
that we have mentioned so far can only cope with
linear equations. What can we do for nonlinear
equations ?
For simplicity, the theory for solving a system
of 3 or more simultaneous nonlinear equations
will be developed for 3x3 system. Extension to
larger systems is straight forward. Consider
the following system of 3 simultaneous nonlinear
equations
(2.27)
34
Do a Taylor series expansion on the left hand
side and truncate at the first order terms to get
where
etc.
We can put the above system of equations into
matrix form
Commonly known as the Jacobian matrix
(2.28)
35
Equation (2.28) is a linear equation so we should
be able to use any of the methods discussed
earlier to solve it !
Hence, to solve the system of equations (2.27),
do the following
1) Using Eq. (2.27), calculate all the
derivatives that are required by Eq.
(2.27). 2) Guess the (initial) values for xi, yi
and zi. 3) Solve Eq. (2.28) to obtain (xi1-
xi), (yi1- yi) and (zi1 -zi). You can
use Gauss elimination, LU decomposition,
Gauss-Seidel etc. 4) Calculate xi1 , yi1 ,
zi1 since you know xi, yi and zi from step
2 above. 5) Repeat steps 2, 3 and 4 until
converge.
This method is the multidimensional version of
the Newton-Raphson method mentioned earlier.
36
Example Solve
Using the multi-dimensional Newton-Raphson method.