This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Numerical Methods in Engineering
ENGR 391 System of linear equations Faculty of Engineering and Computer Sciences
Concordia University LU decomposition
Recall from last lecture • We have seen how we can decompose
[A]=[L][U]
• Useful to solve systems of equations Ax=b: Ax=LUx=b
Write z=Ux => Lz=b
Step 1: Decompose the matrix A by LU decomposition: A=LU
Step 2: Compute z by forward substitution using Lz=b
Step 3: Compute x by backward substitution using Ux=z. Problem
• Since LU decomposition is related to Gauss
elimination
same pitfalls
Division by zero
0 1 1 1 1 1 1 1 1 How can we obtain the LU decomposition of the
above matrix? Partial pivoting
• Remember that we introduced partial pivoting
strategy in the Gauss algorithm to avoid
numerical problems
• How can we implement partial pivoting in the
A=LU decomposition ?
• This can be done with the help of permutation
matrices. Permutation matrix
Definition:
A permutation matrix P is a matrix consisting of all zeros, except for a single 1 in
every row and column.
Equivalently, a permutation matrix is created by applying arbitrary row exchanges for
the identity matrix. 10 01 01 10 001 001 010 100 010 100 010 100 001 001 100 010 100 010 100 010 001 001 Theorem:
Let P be a permutation matrix formed by a particular set of row exchanges
applied to the identity matrix. Then, for any matrix A, PA is the matrix obtained
applying exactly the same set of row exchange to A. Example: 100ab c ab c 001d e f g h i 010g h i d e f PA=LU decomposition
Example: find the PA=LU decomposition of following matrix:
21
A 44 5
4 13 1 5 1
1
4
1
2 Solution:
01021
00144
10013 4
1 0 044 1 002 2 100 8 1
2 4 Solving equations with the PA=LU
decomposition
• Step 1: Decompose the matrix A by
PA=LU decomposition: PA=LU
• Step 2: Compute z by forward substitution
using Lz=Pb
• Step 3: Compute x by backward
substitution using Ux=z. An example
2 x 1y 5 z 5
4x 4 y 4z 0
x 3y z
21
44
13
[A] 5 6 x 5 4y 0 1 z
{x} 6
= {b} Step 1: PA=LU
01021
00144
10013 5
4
1 1
1
4
1
2 0 044 1 002 2 100 8 1
2 4 Step 2: Compute z with Lz=Pb
1
1
4
1
2 0
1
1
2 0 z1
0 z2
z3
1 0105
0010
1006 z1
z2 0
6 z3 8 0
6
5 Step 3: Compute x with Ux=z
44 4x 0 02 2 y 6 00 8 z 8 x 1 y 2 z 1 Matrix inverse
If a matrix [A] is square, there is another matrix [A]1, called the
inverse of [A]; such as:
[A][A]1=[A]1[A]=[I] ; Identity matrix 10
AA 1 A 1 A 01 0
0 0 0 1 0
00 Matrix inverse
To compute the inverse matrix, the first column of [A]1 is
obtained by solving the problem (for 3 3 matrix):
Ax b 1
0
0 ; second column: Ax ; and third column: A x b b 0
1
0
0
0
1 The best way to implement such a calculation is to use LU
decomposition. An other way to read the definition
• Each column j of A1 is the solution of :
0 0
Ax bj 1
0 0 j Computing A1
Without Partial Pivoting
• Step 1: Decompose A=LU
• Step 2: Solve all equations Ax=bj
 Compute z by forward substitution using Lz=bj
 Compute x by backward substitution using Ux=z
• Step 3: Write A1
With Partial Pivoting
• Step 1: Decompose PA=LU
• Step 2: Solve all equations PAx=Pbj
 Compute z by forward substitution using Lz=Pbj
 Compute x by backward substitution using Ux=z
• Step 3: Write A1 An example
1
A 2 1 2 1 2 31 1 Step 1: A=LU
Without Partial Pivoting 1 2 1 1 0 2 1 2 2 1 00
7
10
3 31 1 3 01 2
3
0 1
0
2 Step 2: Solve all Ax=bj
1 2 1x 1 1 2 1x 0 2 1 2y 0 2 1 2y 1 z 0 z 0 31 1 31 1 2 1x 0 2 1 2y 0 z 1 31 1 1 Step 2: Solve all Ax=bj
1
2 2
1 31 1x
2y
1 1
0 z 1 0 z1 1 z1 2 0 0 1 0 z2
7
1 z3
3 0
0 z2
z3 3 1 2 1x 0
0 3
0 0y
2z 1
2
5
3 x
y
z 1
2
2
3
5
6 1
2
5
3 Step 2: Solve all Ax=bj
1
2 2
1 31 1x
2y
1 0
1 z 1
0
0 1 2
3
0 0 z1 0 z1 0 2 0 0 1 0 z2
7
1 z3
3 1
0 z2
z3 1
7
3 3 1x
0 y 2z 0 x 1
7
3 y
z 1
2
1
3
7
6 Step 2: Solve all Ax=bj
1
2 2
1 31 1x
2y
1 0
0 z 1
0
0 1 2
3
0 0 z1 0 2
1 0 1 0 z2
7
1 z3
3 0
1 3 1x 0 x y 0 y 2z 1 z 0 z1
z2
z3 1
2
0
1
2 0
0
1 Step 3: Write A1
x
y
z 1
2
2
3
5
6 1
2
1
3
7
6 x
y
z 1 2 1 2 1 2 31 1 1 1
2
2
3
5
6 x
y
z 1
2
1
3
7
6 1
2
0
1
2 1
2
0
1
2 GaussSeidel Method
• An iterative method
Approach
• The system Ax=b is reshaped by solving the first
equation for x1, the second equation for x2, and
the third for x3, …and nth equation for xn. GaussSeidel Method
For a 3x3 system x1 b1 a12 x2 a13 x3
a11 x2 b2 a21 x1 a23 x3
a22 x3 b3 a31 x1 a32 x2
a33 A simple way to obtain initial guesses for all x’s to start the solution process is
to assume that they are zero.
These zeros can be substituted into x1 equation to calculate a new x1=b1/a11.
Then we substitute this new value of x1 along with the previous guess of zero
for x3 into the second equation to find x2. The process is repeated for the third
equation to calculate a new estimate for x3 GaussSeidel Method
• New x1 is obtained by using the new estimates of x2 and
x3. The procedure is repeated until the convergence
criterion is satisfied: a ,i xij xij
xi j 1 100% s For all i, where j and j1 are the present and previous
iterations. Fig. 11.4 Gauss Seidel Jacobi method p.291 k
In GaussSeidel method, the values x1k 1 , x 2 1 ,..., xik 1
k
k
k
computed in the current iteration as well as xi 2 , xi 3 ,..., x n k1
are used in finding the value xi 1 This implies that always the most recent approximations
are used during the computation. The general expression is:
k1
i x 1
bi
aii i1
j1 aij x k1
j NEW n
ji1 aij xik ; i=1,2,...,n; k=1,2,3, …..
OLD In Jacobi iteration method, all the new values are computed using the
values at the previous iteration. This implies that both the present and
the previous set of values have to be stored.
GaussSeidel method will improve the storage requirement as well as
the convergence. Convergence Criterion for GaussSeidel Method
• • The GaussSeidel method has two fundamental problems as many
iterative method:
– It is sometimes nonconvergent, and
– If it converges, converges very slowly.
Recalling that sufficient conditions for convergence of two linear
equations, u(x,y) and v(x,y) are u
x u
y 1 v
x v
y 1 • Similarly, in case of two simultaneous equations,
the GaussSeidel algorithm can be expressed as
u ( x1 , x2 ) b1
a11 a12
x2
a11 v( x1 , x2 ) b2
a22 a21
x1
a22 u
x1
v
x1 u
x2 0
a21
a22 v
x2 a12
a11
0 • Substitution into convergence criterion of two linear
equations yield:
a12
a11 a21
a22 1 1 • In other words, the absolute values of the slopes
must be less than unity for convergence:
a11 a12 a22 a21 For n equations :
n aii ai , j
j1
ji Figure 11.5 ...
View Full
Document
 Winter '09
 HOIDICK
 Computer Science

Click to edit the document details