Download Gaussian Elimination - Computational Methods, Notes | MCEN 3030 and more Study notes Mechanical Engineering in PDF only on Docsity!
Lectures 5 to 7
Direct methods for solving Ax = b
Notation: A ∈ R
m×m and x, b ∈ R
m
A =
a 11 a 12 · · · · · · a 1 m
a 21 a 22 a 2 m
. .
. a 33
am 1 am 2 · · · · · · amm
, b =
b 1
b 2
. . .
bm− 1
bm
, x =
x 1
x 2
. . .
xm− 1
xm
Ax = b represents a system of linear equations to solve for x:
a 11 x 1 + a 12 x 2 + · · · + a 1 mxm = b 1
a 21 x 1 + a 22 x 2 + · · · + a 2 mxm = b 2
. . .
am 1 x 1 + am 2 x 2 + · · · + ammxm = bm
Gaussian elimination
In the above system it is possible to eliminate all the ak 1 x 1 by multiplying
the first line by
ak 1 a 11
and subtracting it to the k-th line. Doing so will lead to
the modified system:
a 11 x 1 + a 12 x 2 + · · · + a 1 mxm = b 1
0 + a
(1) 22 x^2 +^ · · ·^ +^ a
(1) 2 mxm^ =^ b
(1) 2 . . .
0 + a
(1) m 2 x^2 +^ · · ·^ +^ a
(1) mmxm =^ b
(1) m
The modified values a
(1) and b
(1) are resulting from the subtraction
a
(1) ij =^ aij^ −^ a^1 j
ai 1
a 11
and
b
(1) k =^ bk^ −^ b^1
ak 1
a 11
where (i, j) > 1 and k > 1
The step of eliminating the ak 1 for k > 1 in the system can be written as
a matrix operation called a Gaussian transformation:
M 1 =
−a 21 /a 11 1 0 0
−a 31 /a 11 0 1
−am 1 /a 11 0 · · · · · · 1
or equivalently as
M 1 = I + m 1 e
T 1 (3)
where I is the m × m identity matrix, e
T 1 = (1,^0 ,^0 ,^0 ,^ · · ·^ ,^ 0) and^ m
T 1 =
(0, −a 21 /a 11 , −a 31 /a 11 , · · · , −am 1 /a 11 ).
Thus,
M 1 A =
a 11 a 12 · · · a 1 m
0 a
(1) 22 · · ·^ a
(1) 2 m . . .
0 a
(1) m 2 · · ·^ a
(1) mm
and
M 1 b =
b 1
b
(1) 2
b
(1) 2 . . .
b
(1) m
The second elimination step is then written as :
M 2 M 1 Ax = M 2 M 1 b (4)
where the Gauss elimination matrix is
M 2 =
0 −a
(1) 32 /a
(1) 22 1
0 −a
(1) m 2 /a
(1) 22 0 · · ·^1
k-th positions (1 ≤ l ≤ k). Carrying the procedure until the last block is
reached brings to the upper triangular system
U x = Mm− 1 Mm− 2 · · · M 2 M 1 Ax = Mm− 1 Mm− 2 · · · M 2 M 1 b =
b (6)
or simply U x = ˆb. Explicitly the matrix U is
U =
a 11 a 12 · · · · · · a 1 m
0 a
(1) 22 · · ·^ · · ·^ a
(1) 2 m . .
. 0 a
(2) 33 · · ·^ a
(2) 3 m . . .
0 0 · · · · · · a
(m−1) mm
It is easy to see that xm = b
(m−1) m /a
(m−1) mm since the vector
ˆb is
ˆb =
b 1
b
(1) 2 . . .
b
(m−1) m
In a similar fashion, since we now know xm, it is possible to obtain xm− 1 :
xm− 1 =
b
(m−2) m− 1 −^ a
(m−2) m− 1 ,mxm
a
(m−2) m− 1 ,m− 1
The multiplications of A by the matrices Mk transforming A to U is the
(Gauss) elimination step. Solving for x with the matrix U is the back(ward)substitution
step. It is possible to write a mathematical formula for the backsubstitution:
xm−k =
b
(m−k−1) m−k
m j=m−k+
a
(m−k−1) m−k,j
xj
a
(m−k−1) m−k,m−k
The pseudo code performing the transformation of A into U is
given: a(,), b() and m
do k=1,...,m-
do i=k+1,...,m
fac = a(i,k)/a(k,k)
do j=k+1,...,m
a(i,j) = a(i,j) - fac * a(k,j)
end do
b(i) = b(i) - fac * b(k)
end do
end do
the output of the above code returns U inside A (the upper triangular part of
A is replaced by U ) and ˆb. The pseudo code performing the backsubstitution
is simply
do k=m,m-1,...,
x(k) = b(k)
do i=k+1,...,m
x(k) = x(k)-a(k,i) * x(i)
end do
x(k) = x(k)/a(k,k)
end do
Up until now we have assumed that the divisions by diagonal elements during
the elimination process was defined (all non-zeros). In the case of a zero
diagonal element, one needs to pivot lines or rows in the matrix. The subject
will be discussed later.
LU factorization
Gaussian elimination is not very attractive when the system needs to be
solved for multiple right hand sides (rhs). This is obvious since b is trans-
formed to ˆb using the Gaussian transformation matrices Mk. This change of
b is also depicted in the Gaussian elimination pseudo-code. In what follows
we will show that in fact the inverse of the combined Mk putting A in upper
triangular form is a lower triangular matrix that we will call L. Generating L
requires only a simple modification to the Gaussian elimination pseudo-code.
Thus the goal of this section is to transform A into LU or in other words
construct the identity A = LU. This is quite useful if the following steps are
considered:
The Gauss elimination pseudo code can be modified to create the L and U
matrices. Without any pivoting strategy, we are going to call this algorithm
naive LU :
given: a(,) and m
set: L(,) = 0
do k=1,...,m-
do i=k+1,...,m
fac = a(i,k)/a(k,k)
L(i,k) = fac
a(i,k) = 0
do j=k+1,...,m
a(i,j) = a(i,j) - fac * a(k,j)
end do
end do
end do
do k=1,...,m
L(k,k)=
end do
at completion, A is transformed into U with zeros in its lower part and L is
created.
Pivoting strategies
The pivot is the element of the matrix with whom we perform divisions in
both the Gaussian elimination and the LU factorization. A zero pivot obvi-
ously pose problem. All hopes are not lost since if A is regular then a non
zero pivot must exist in the same column.
- Partial pivoting: before eliminating values in column k, find the largest
|a
(k−1) ik |^ for^ i^ =^ k, k^ + 1, ..., m. The largest value becomes the pivot instead of
the naive a
(k−1) kk
- Full pivoting: before eliminating values in column k, find the largest
max(|a
(k−1) ik
|, |a
(k−1) ki
|) for i = k, k + 1, ..., m. The largest value becomes the
pivot instead of the naive a
(k−1) kk.
We will interest ourselves to the first approach only since it involves
only row interchanges while the second performs both row and column inter-
changes.
LU with partial pivoting:
given: a(,) and m
set: L(,) = 0
do k=1,...,m
P(k)=k
end do
do k=1,...,m-
q = findpivot(a(:,k),k)
tmp = P(k)
P(k) = P(q)
P(q) = tmp
Exchange: row k with row q in a(,)
Exchange: row k with row q in L(,)
do i=k+1,...,m
fac = a(i,k)/a(k,k)
L(i,k) = fac
a(i,k) = 0
do j=k+1,...,m
a(i,j) = a(i,j) - fac * a(k,j)
end do
end do
end do
do k=1,...,m
L(k,k)=
end do
The pseudo code performing the forward and backward substitution is simply
given: P, L, U, b and m
do k=1,...,m
bnew(k)=b(P(k))
end do
do k=1,2,...,m
y(k) = bnew(k)
Examples
- Gaussian elimination with partial pivoting
A =
the pivot in the first column is 6:
P 1 =
P 1 A =
the Gaussian transformation is
M 1 =
we get
M 1 P 1 A =
The second pivot is 8 in the second column
P 2 =
thus
P 2 M 1 P 1 A =
the second Gauss transformation is
M 2 =
and finally we obtain:
M 2 P 2 M 1 P 1 A =
Note that if given a b then ˆb = M 2 P 2 M 1 P 1 b.
- LU with partial pivoting Using the same matrix but not considering the
ones in the Mk
the permuted inverse Gaussian transformation is (minus identity)
P 1 (M
− 1 1 −^ I) =
proceed as before (last step) no reordering
P 1 (M
− 1 1 −^ I) + (M^
− 1 2 −^ I) =
adding the diagonal and multiplying by U on the right:
With P = P 2 P 1
multiplying A with P
