Linear Equations-Direct Methods
Linear Equations-Direct Methods
Linear Equations-Direct Methods
Computational Methods
(SECE)
Yonas Y.
Computational Methods (SECE) Chapter 5 Yonas Y. 1 / 50
Outline
3 Pivoting strategies
Partial pivoting
GEPP stability
Matrices requiring no pivoting
Introduction
Kirchhoff’s laws of electrical circuits state that both the net flow of current
through each junction and the net voltage drop around each closed loop of a
circuit are zero.
Suppose we have a circuit shown in the diagram below.
Using G as a reference point, Kirchhoff’s laws imply that the currents satisfy
the following system of linear equations:
5i1 + 5i2 = V ,
i3 − i4 − i5 = 0,
2i4 − 3i5 = 0,
i1 − i2 − i3 = 0,
5i2 − 7i3 − 2i4 = 0.
Computational Methods (SECE) Chapter 5 Yonas Y. 3 / 50
Linear Systems: Direct Methods Introduction
or
Ax = b.
The matrix elements aik and the right-hand side elements bi are given.
We are looking for the unknowns xk .
Ax = b
Backward substitution
Occasionally, A has a special structure
If A is diagonal
a11
a22
A=
..
.
ann
then the linear equations are uncoupled, and the solution is
bi
xi = , i = 1, 2, . . . , n.
aii
If A is an upper triangular matrix
a11 a12 · · · a1n
. ..
a22 . .
.
A= .
,
.. ..
.
ann
Where aij = 0 for all i > j.
Computational Methods (SECE) Chapter 5 Yonas Y. 6 / 50
Gaussian elimination and backward substitution Backward substitution
Note that the successful completion of the algorithm is possible only if the
diagonal elements of A, aii , are nonzero.
In addition, data locality and data movement are crucial to the execution of
an algorithm.
Forward substitution
5x1 = 15,
x1 + 2x2 = 7,
−x1 + 3x2 + 2x3 = 5.
Gaussian elimination
Elementary row transformations are used to reduce the given linear system to an
equivalent upper triangular form.
The reason this works is that the solution to our problem Ax = b is invariant to
multiplication of a row by a constant,
subtraction of a multiple of one row from another, and
row interchanges.
x1 x2 x3 x4
a11 a12 a13 a14 b1
a21 a22 a23 a24 b2
a31 a32 a33 a34 b3
a41 a42 a43 a44 b4
x1 x2 x3 x4
a11 a12 a13 a14 b1
0 0 0 0
0 a22 a23 a24 b2
0 0 0 0
0 a32 a33 a34 b3
0 0 0 0
0 a42 a43 a44 b4
0
1 Permute rows i = 2, . . . , 4 (if necessary) such that a22 6= 0. This element is
the next pivot.
0 0 0
2 Subtract multiples li2 = ai2 /a22 of row 2 from row i, i = 3, . . . , 4.
00 0 0 0
3 Set aik = aik − li2 a2k , i, k = 3, . . . , 4.
00 0 0 0
4 Set bi = bi − li2 b2 , i = 3, ..., 4.
x1 x2 x3 x4
a11 a12 a13 a14 b1
0 0 0 0
0 a22 a23 a24 b2
00 00 00
0 0 a33 a34 b3
00 00 00
0 0 a43 a44 b4
00
1 Permute rows i = 3, . . . , 4 (if necessary) such that a33 6= 0. This is the next
pivot.
00 00 00
2 Subtract multiples l43 = a43 /a33 of row 3 from row 4.
000 00 00 00
3 Set a44 = a44 − l43 a34 .
000 00 00 00
4 Set b4 = b4 − l43 b3 .
x1 − 4x2 = −10
0.5x1 − x2 = −2
Example 5.3. Use Gaussian elimination with backward substitution to solve the
following linear system.
4x1 − x2 + x3 = 8,
2x1 + 5x2 + 2x3 = 3,
x1 + 2x2 + 4x3 = 11.
We have n = 3; hence the loops in the general algorithm collapse into defining
after two update’s for k = 1, 2, i = 2, 3, j = 2, 3:
For k = 1,
a21 a31
l21 = = 0.5, l31 = = 0.25.
a11 a11
0 0 0
a22 = 5−0.5·(−1) = 5.5, a23 = 2−0.5·(1) = 1.5, b2 = 3−0.5·(8) = −1.
0 0 0
a32 = 2−0.25·(−1) = 2.25, a33 = 4−0.25·(1) = 3.75 b3 = 11−0.25·(8) = 9.
Computational Methods (SECE) Chapter 5 Yonas Y. 17 / 50
Gaussian elimination and backward substitution Gaussian elimination
Thus
0 0
4 −1 1 8
(A |b ) = 0 5.5 1.5 −1 ,
0 2.25 3.75 9
For k = 2,
0
0 a 2.25
l32 = 32
0 = = 0.41.
a22 5.5
00 00
a33 = 3.75 − 0.41 · (1.5) = 3.135, b3 = 9 − 0.41 · (−1) = 9.41.
Thus
00 00
4 −1 1 8
(A |b ) = 0 5.5 1.5 −1 ,
0 0 3.135 9.41
LU decomposition
Gaussian elimination decomposes matrix A into a product L × U of a lower
triangular matrix L and an upper triangular matrix U.
This matrix has zeros everywhere except in the main diagonal and the first
column.
Then,
A(1) = M (1) A and b(1) = M (1) b
Likewise, to zero out the second column of A(1) below the main diagonal
A(2) = M (2) A(1) = M (2) M (1) A and b(1) = M (2) b(1) = M (2) M (1) b
where
1
1
..
M (2)
=
−l32 . .
.. ..
. .
−ln2 1
and likewise
b(n−1) = M (n−1) · · · M (2) M (1) b
(n−1) −1
Multiplying U by [M ] , then by [M (n−2) ]−1 , etc., we obtain
A = LU
where L = [M (1) ]−1 · · · [M (n−2) ]−1 [M (n−1) ]−1 .
1
1. The Gaussian elimination process for the first column yields l21 = 1 = 1,
l31 = 13 = 3, so
1 0 0 1 −1 3
M (1) = −1 1 0 , A(1) = M (1) A = 0 2 −3 .
−3 0 1 0 1 −8
2. The Gaussian elimination process for the second column yields l32 = 12 , so
1 0 0 1 −1 3
M (2) = 0 1 0 , U = A(2) = M (2) A(1) = 0 2 −3 .
0 −0.5 1 0 0 −6.5
3. Note that
1 0 0 1 0 0
[M (1) ]−1 = 1 1 0 , [M (2) ]−1 = 0 1 0 ,
3 0 1 0 0.5 1
1 0 0
[M (1) ]−1 [M (2) ]−1 = 1 1 0 = L,
3 0.5 1
and
1 0 0 1 −1 3 1 −1 3
LU = 1 1 0 0 2 −3 = 1 1 0 = A
3 0.5 1 0 0 −6.5 3 −2 1
A = LU.
Example 5.5. Determine the LU factorization for matrix A in the linear system
Ax = b, where
1 1 0 3 8
2 1 −1 1 7
A= 3 −1 −1 2
and b= 14 .
−1 2 3 −1 −7
and solve the system.
To solve
1 0 0 0 1 1 0 3 x1 8
2 1 0 0 0 −1 −1 −5 x2 7
Ax = LUx =
3
= .
4 1 0 0 0 3 13 x3 14
−1 3 0 1 0 0 0 −13 x4 −7
Computational Methods (SECE) Chapter 5 Yonas Y. 25 / 50
Gaussian elimination and backward substitution LU decomposition
We first introduce the substitution y = Ux. Then b = L(Ux) = Ly. That is,
1 0 0 0 y1 8
2 1 0 0 y2 7
Ax = LUx = 3 4 1 0 y3 = 14 .
−1 3 0 1 y4 −7
We then solve Ux = y for x, the solution of the original system; that is,
1 1 0 3 x1 8
0
−1 −1 −5 x2 = −9 .
0 0 3 13 x3 26
0 0 0 −13 x4 −26
Pivoting strategies
The process of Gaussian elimination (or LU decomposition) cannot be applied
unmodified to all nonsingular matrices.
For instance
0 1 x1 4
Ax = b ⇐⇒ =
1 1 x2 7
x1 x2
0 1 4
1 1 7
x1 x2 x1 x2
0.00035 1.2654 3.5267 0.00035 1.2654 3.5267
1.2547 1.3182 6.8541 0 -4535.0 -12636
with
l21 = 1.2547/0.00035 ≈ 3584.9
0
a22 = 1.3182 − 3584.9 × 1.2654 ≈ 1.3182 − 4536.3 ≈ −4535.0
0
b2 = 6.8541 − 3584.9 × (3.5267) ≈ 6.8541 − 12643 ≈ −12636.
Backsubstitution gives
x2 = −12636/(−4535.0) ≈ 2.7863
x1 = (3.5267 − 1.2654 × 2.7863)/0.00035
≈ (3.5267 − 3.5258)/0.00035 = 0.0009/0.00035 = 2.5714.
x1 = 2.5354, x2 = 2.7863
Partial pivoting
(k−1)
If a pivotal element akk is small in magnitude
Roundoff errors are amplified, both in subsequent stages of the Gaussian
elimination process and (even more apparently) in the backward substitution
phase.
The most common strategy for pivoting is to search for the maximal element in
modulus:
The index p of the pivot in the k-th step of GE is determined by
(k−1) (k−1)
|apk | = max |aik |.
i≥k
Let’s consider example 5.7 using Gaussian elimination with partial pivoting(GEPP)
x1 x2 x1 x2
1.2547 1.3182 6.8541 1.2547 1.3182 6.8541
0.00035 1.2654 3.5267 0 1.2650 3.5248
with
l21 = 0.00027895
x2 = 2.7864
x1 = (6.8541 − 1.3182 × 2.7864)/1.2547
≈ (6.8541 − 3.6730)/1.2547 = 3.1811/1.2547 ≈ 2.5353.
There is a deviation from the exact solution in the last digit only.
GEPP stability
Complete pivoting
In complete pivoting we look for (one of) the largest element in modulus
(k−1) (k−1)
|aqr | = max |aij |,
k≤i,j≤n
at each stage of GE and interchange both row q and column r with row i and
column j, respectively.
Theorem
If A is nonsingular and diagonally dominant then the LU factorization can be
computed without pivoting.
xT Ax > 0, ∀x 6= 0.
A = LDLT ,
Note that L is a unit lower triangular matrix and D is a diagonal matrix with
positive diagonal entries.
Then,
A = GG T ,
Then
−10−8
br = b − Ab
x= =⇒ k br k∞ = 10−8 .
10−8
However
2
x= , =⇒ k x−b
x k∞ = 1.513.
−2
Thus, the error is of the same size as the solution itself, roughly 108 times that of
the residual!
Therefore
kzk kb
x−xk k br k k br k
= ≤ k A kk A−1 k =: κ(A)
kxk kxk kbk kbk
Definition
We therefore define the condition number of the matrix A as
κ(A) =k A kk A−1 k
is called the condition number of A.
The relative error in the solution is bounded by the condition number of the
matrix A times the relative error in the residual.
(A + δA)b
x = b + δb, x = x + δx.
b
br = b − Ab x − δb.
x = (δA)b
k δA k∞ ≤ ηφ(n)gn (A),
where
φ(n) is a low order polynomial in n (cubic at most) and
η is the rounding unit.
The extra perturbations which the forward and backward substitutions yield,
in particular the bound on k δb k∞ , are significantly smaller.
k x−b
xk k br k
≤ κ(A)
kxk kbk
we finally obtain
k x−b
xk k δb k + k δA kk b
xk
≤ κ(A)
kxk kbk
End of Chapter 5
Questions?