Download Gaussian Elimination and LU-Factorization: Solving Linear Systems of Equations - Prof. Lot and more Study notes Mathematics in PDF only on Docsity!
3 Gaussian elimination and LU-factorization
Many computational problems require the solution of linear systems of equations
Ax = b, (1)
where
A =
a 1 , 1 a 1 , 2... a 1 ,n− 1 a 1 ,n a 2 , 1 a 2 , 2... a 2 ,n− 1 a 2 ,n .. .
an− 1 , 1 an− 1 , 2... an− 1 ,n− 1 an− 1 ,n an, 1 an, 2... an,n− 1 an,n
and b =
b 1 b 2 .. . bn− 1 bn
denote an n × n matrix with real entries ai,j and an n-vector with real entries bj. The matrix A typically represents a model and the right-hand side b represents data to which we seek to fit the model. We will usually refer to A as the model matrix and to b as the data. In some applications the data is contaminated my errors, e.g., measurement errors, and in those cases we are interested in how much the solution is affected by these errors. This topic will be discussed in Chapter 6. The number of equations n may be large for many practical scientific applications. For instance, problems in structural analysis frequently give rise to linear systems of equations of order about 10^6. These systems are commonly solved on computers with several processors. Some applications require the solution of linear systems of equations with 10^8 or more unknowns. The solution of so large systems is a challenge even on the fastest computers available. We will discuss methods for the solution of linear systems of equations based on Gaussian elimination. The basic ideas will be presented for small systems, starting with systems of of order two.
3.1 Linear systems of equations of order two
Consider the linear system of equations
a 1 , 1 x 1 + a 1 , 2 x 2 = b 1 (2) a 2 , 1 x 1 + a 2 , 2 x 2 = b 2 (3)
Each equation represents a line in the (x 1 , x 2 )-plane. The solution of this system amounts to the familiar problem of finding the point at which two lines intersect. If these lines intersect at unique point, then the linear system of equations has a unique solution. On the other hand, if the lines are parallel then two cases have to be distinguished: i) if the lines are distinct, then the linear system of equations has no solution, and ii) if the lines coalesce, then the system has infinitely many solutions. The angle between the lines determines how accurately the intersection can be computed in the presence of errors in the data. Example 1. Consider the linear system of equations
1 · x 1 + 0 · x 2 = 1 0 · x 1 + 1 · x 2 = 2
The first equation represents the vertical line {(x 1 , x 2 ) : x 1 = 1, x 2 ∈ R} and the second equation represents the horizontal line {(x 1 , x 2 ) : x 1 ∈ R, x 2 = 2}. These lines are perpendicular and intersect at the point (1, 2); the lines are shown in Figure 1(a). Matrices, such that the solutions of the individual equations intersect
−1 0 1 2 3
−
0
1
2
3
x 1
x^2
−1 0 1 2 3
−
0
1
2
3
x 1
x^2
−1 0 1 2 3
−
0
1
2
3
x 1
x^2
Data error
(a) (b)
Figure 1: Perpendicular lines (a), and with error in our knowledge of the position of the vertical line (b).
at right angles are said to be orthogonal. The matrix in the present example is a very simple orthogonal matrix, the identity matrix. Orthogonality is an important concept that we will study in more detail in later chapters. Orthogonality can be extended to higher-dimensional systems. Figure 1(b) depicts a plot in which the first right-hand side component is perturbed by −1 (labeled as “data error” in the plot. This perturbation may, for instance, stem from a measurement error, and we are interested in how it affects the computed solution. The equation with the perturbed right-hand side component represents the horizontal line {(x 1 , x 2 ) : x 1 = 0, x 2 ∈ R}. The second equation is not changed. The lines represented by these equations now intersect at the point (0, 2). Figure 1(b) shows the new solution with a black disc. 2 In the above example the error 1 in the data, as measured by the Euclidean norm ‖ · ‖, gives rise to an error 1 in the computed solution. This is as good a result as we can hope for and it depends on that the lines represented by the linear system of equations (2)-(3) are perpendicular. When the angle between the lines is acute, a larger error in the computed solution may result. This is illustrated in the following example. Example 2. Consider the linear system of equations
x 1 + 0 · x 2 = 1 −x 1 + x 2 = 1
The lines corresponding to the two equations are depicted in Figure 3.1(a). As in the last example, Figure 3.1(b) shows the result of perturbing the first component of the data by −1. However, unlike in the last example, we see that the Euclidean-norm error in the computed solution is now larger than the error in the data. (Elementary trigonometry shows that the computed solution error is
As Examples 1 and 2 illustrate, errors in the data may be maginfied in the computed solution when the
In many applications one would like to apply the same model to more than one set of measured data vectors. One then needs to solve the system (1) for several right-hand sides, all of which may not be available at the same time. The upper triangular matrix U is independent of the right-hand side vector b in (1), however, the vector b˜ in the upper triangular system (4) is not. We would like to determine the vector b˜ for each right-hand side vector b without having to recompute U. It turns out that this is possible and computationally inexpensive. It can be achieved by storing some auxiliary quantities, which are computed during the evaluation of U from A, in an n × n lower triangular matrix with unit diagonal,
L =
ℓn− 1 , 1 ℓn− 1 , 2 · · · 1 ℓn, 1 ℓn, 2 · · · 1
Then we have A = LU. (5)
This factorization is called the LU-factorization of A. Given this LU-factorization of the matrix A, we can solve the linear system (1) in two steps. We have
LU x = b
and solve, in order,
Ly = b, (6) U x = y. (7)
The first system is solved by computing the components of y = [y 1 , y 2 ,... , yn]T^ in the order of increasing index, i.e., y 1 is computed from the first equation of (6), then y 2 is computed from the second equation of (6), and so on, until the last component yn is computed from the last equation in (6). This solution method is known as forward substitution. One can show that the vector y obtained in this manner is the same as the vector b˜ in (4).
3.3 Computation of the LU-factorization
We illustrate the computation of the triangular factor L and U with a 4 × 4 matrix. Let
A =
The first step of Gaussian elimination entails subtracting suitable multiples of the first row from rows 2 through 4, so that the (2, 1), (3, 1), and (4, 1) entries in the matrix so obtained vanish. We will store these multiples in the first column below the diagonal of a lower triangular matrix L with unit diagonal. Thus,
we subtract 2 times the first row from the second row and store the multiple 2 in the (2, 1)-entry of L. This gives the matrices
A′^ =
,^ L^ =
To proceed, we subtract suitable multiples of the second row of the matrix A′^ from rows three and four and store these multiples in the second column below the diagonal of L. Specifically, we subtract 3 times row two from row three and 4 times row two from row four and obtain
A′′^ =
,^ L^ =
Finally, we subtract row three of A′′^ from row four and record the multiple 1 in the (4, 3)-entry of L. This yields the upper triangular matrix U in (4) as well as the lower triangular matrix L,
U =
,^ L^ =
One can verify by direct computation that L and U are the desired factors in the LU-factorization (5) of A. The upper triangular matrix U is determined in the standard way by Gaussian elimination of A, and the lower triangular matrix L is obtained by recording the multiples used in the subtractions in Gaussian elimination. The follwing algorithm determines the LU-factorization of an n × n matrix A.
Algorithm 3.1: Computation of the LU-factorization
INPUT: n × n matrix A = [aj,k]nj,k=1.
OUTPUT: Lower triangular n×n matrix L = [ℓj,k]nj,k=1, upper triangular n×n matrix U = [uj,k]nj,k=1.
Let L be the n × n identity I. Let U := A.
For k = 1, 2 ,... , n − 1,
For j = k + 1, k + 2,... , n, ℓj,k := (^) uuk,kj,k. For i = k + 1, k + 2,... , n, uj,i := uj,i − ℓj,kuk,i. End i End k
End j
