Heath - Scientific Computing (523150), страница 15

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 15)

We first choose an elementary elimination matrixM1 according to the recipe given in Section 2.2.2, with the first diagonal entry a11 as pivot,so that, when premultiplied by M1 , the first column of A becomes zero below the first row.Of course, all of the remaining columns of A, as well as the right-hand-side vector b, arealso multiplied by M1 , so the new system becomes M1 Ax = M1 b, but by our previousdiscussion the solution is unchanged.Next we use the second diagonal entry as pivot to determine a second elementaryelimination matrix M2 that annihilates all of the entries of the second column of thenew matrix, M1 A, below the second row. Again, M2 must be applied to the entirematrix and right-hand-side vector, so that we obtain the further modified linear systemM2 M1 Ax = M2 M1 b. Note that the first column of the matrix M1 A is not affected byM2 because all of its entries are zero in the relevant rows.

This process is continued for eachsuccessive column until all of the subdiagonal entries of the matrix have been annihilated,so that the linear systemM Ax = Mn−1 · · · M1 Ax = Mn−1 · · · M1 b = M bis upper triangular and can be solved by back-substitution to obtain the solution to theoriginal linear system Ax = b.2.2. SOLVING LINEAR SYSTEMS43The process we have just described is known as Gaussian elimination. It is also knownas LU factorization or LU decomposition because it decomposes the matrix A into a productof a unit lower triangular matrix, L, and an upper triangular matrix, U . To see this, recallthat the product Lk Lj is unit lower triangular if k < j, so that−1L = M −1 = (Mn−1 · · · M1 )−1 = M1−1 · · · Mn−1= L1 · · · Ln−1is unit lower triangular. We have already seen that, by design, the matrix U = M A isupper triangular.

Therefore, we have expressed A as a productA = LU ,where L is unit lower triangular and U is upper triangular. Given such a factorization,the linear system Ax = b can then be written as LU x = b and hence can be solved byfirst solving the lower triangular system Ly = b by forward-substitution, then the uppertriangular system U x = y by back-substitution. Note that the intermediate solution yis the same as the transformed right-hand-side vector, M b, in the previous formulation.Thus, Gaussian elimination and LU factorization are simply two ways of expressing thesame solution process.Example 2.6 Gaussian Elimination. We illustrate Gaussian elimination by solving thelinear system2x1 + 4x2 − 2x3 = 2,4x1 + 9x2 − 3x3 = 8,−2x1 − 3x2 + 7x3 = 10,or in matrix notation2Ax =  4−2   4 −2x129 −3   x2  =  8  = b.−37x310To annihilate the subdiagonalfirst row from the second row,1M1 A = −21entries of the first column of A, we subtract two times theand add the first row to the third row: 0 024 −22 4 −21 0 49 −3  =  0 11,0 1−2 −370 15   1 0 022M1 b =  −2 1 0   8  =  4  .1 0 11012Now to annihilate the subdiagonal entry of the second column of M1 A, we subtract thesecond row from the third row: 10 02 4 −22 4 −2M2 M1 A =  01 00 11 = 0 11,0 −1 10 150 0444CHAPTER 2.

SYSTEMS OF LINEAR EQUATIONS1M2 M1 b =  00   0 0221 0 4 = 4.−1 1128We have therefore reduced the original system to the equivalent upper triangular system   2 4 −22x10 11   x2  =  4  ,x30 048which can now be solved by back-substitution (as in Example 2.4) to obtain x = [ −1To write out the LU factorization explicitly, we have 1 0 01 0 01 0 0L = L1 L2 =  2 1 0   0 1 0  =  2 1 0  ,−1 0 10 1 1−1 1 1so that2A= 4−22.2.4 4 −21 09 −3  =  2 1−37−1 102001041022 ]T .−21  = LU .4PivotingThere is one obvious problem with the Gaussian elimination process as we have described it,as well as another, somewhat more subtle, problem. The obvious potential difficulty is thatthe process breaks down if the leading diagonal entry of the remaining unreduced portionof the matrix is zero at any stage, as computing the multipliers mi for a given columnrequires division by the diagonal entry in that column.

The solution to this problem isalmost equally obvious: if the diagonal entry is zero at stage k, then interchange row k ofthe system with some subsequent row whose entry in column k is nonzero. As we knowfrom Example 2.2, such an interchange does not alter the solution to the system. With anonzero diagonal entry as pivot, the process can then proceed as usual.But what if there is no nonzero entry on or below the diagonal in column k? Then thereis nothing to do at this stage, since all the entries to be annihilated are already zero, andwe can simply move on to the next column (i.e., Mk = I). Note that this step leaves azero on the diagonal, and hence the resulting upper triangular matrix U is singular, butthe LU factorization can still be completed. It does mean, however, that the subsequentback-substitution process will fail, since it requires a division by each diagonal entry of U ,but this is not surprising because the original matrix must have been singular anyway.

Amore insidious problem is that in floating-point arithmetic we may not get an exact zero,but only a very small diagonal entry, which brings us to the more subtle point.In principle, any nonzero value will do as the pivot for computing the multipliers, but inpractice the choice should be made with some care to minimize error. When the remainingportion of the matrix is multiplied by the resulting elementary elimination matrix, weshould try to limit the growth of the entries of the transformed matrix in order not to2.2. SOLVING LINEAR SYSTEMS45amplify rounding errors. For this reason, it is desirable for the multipliers not to exceed 1in magnitude. This requirement can be met by choosing the entry of largest magnitude onor below the diagonal as pivot. Such a policy is called partial pivoting, and it is essential inpractice for a numerically stable implementation of Gaussian elimination for general linearsystems.The row interchanges required by partial pivoting slightly complicate the formal description of LU factorization given earlier.

In particular, each elementary elimination matrix Mkis preceded by a permutation matrix Pk that interchanges rows to bring the entry of largestmagnitude into the diagonal pivot position. We still have M A = U , where U is uppertriangular, but nowM = Mn−1 Pn−1 · · · M1 P1 .M −1 is still triangular in the general sense defined earlier, but because of the permutations,M −1 is not necessarily lower triangular, though we still denote it by L.

Thus, “LU”factorization no longer literally means “lower times upper” triangular, but it is still equallyuseful for solving linear systems by successive substitution.We note that the permutation matrixP = Pn−1 · · · P1permutes the rows of A into the order determined by partial pivoting. An alternativeinterpretation, therefore, is to think of partial pivoting as a way of determining a rowordering for the system under which no pivoting would be required for numerical stability.Thus, we obtain the factorizationP A = LU ,where now L really is lower triangular.

To solve the linear system Ax = b, we first solvethe lower triangular system Ly = P b by forward-substitution, then the upper triangularsystem U x = y by back-substitution.The name “partial” pivoting comes from the fact that only the current column issearched for a suitable pivot.

A more exhaustive pivoting strategy is complete pivoting,in which the entire remaining unreduced submatrix is searched for the largest entry, whichis then permuted into the diagonal pivot position. Note that this requires interchangingcolumns as well as rows, and hence it leads to a factorization of the formP AQ = LU ,where L is unit lower triangular, U is upper triangular, and P and Q are permutationmatrices that reorder the rows and columns, respectively, of A. To solve the linear systemAx = b, we first solve the lower triangular system Ly = P b by forward-substitution,then the upper triangular system U z = y by back-substitution, and finally we permutethe solution components to obtain x = Qz.

Although the numerical stability of completepivoting is theoretically superior, it requires a much more expensive pivot search thanpartial pivoting. Because the numerical stability of partial pivoting is more than adequatein practice, it is almost universally used in solving linear systems by Gaussian elimination.Since pivot selection depends on the magnitudes of individual matrix entries, the particular choice obviously depends on the scaling of the matrix. A diagonal scaling of the46CHAPTER 2. SYSTEMS OF LINEAR EQUATIONSmatrix (recall Example 2.3) may result in a different sequence of pivots. For example, anynonzero entry in a given column can be made the largest in magnitude simply by giving thatrow a sufficiently heavy weighting.

This does not mean that an arbitrary pivot sequenceis acceptable, however: a badly skewed scaling can result in an inherently sensitive systemand a correspondingly inaccurate solution. A well-formulated problem should have appropriately commensurate units for measuring the unknown variables (column scaling), and aweighting of the individual equations that properly reflects their relative importance (rowscaling).

It should also account for the relative accuracy of the input data. Under thesecircumstances, the pivoting procedure will usually produce a solution that is as accurate asthe problem warrants (see Section 2.4).Example 2.7 Pivoting. Here are some examples to illustrate the necessity of pivoting,both in theory and practice, for a stable implementation of Gaussian elimination. We firstobserve that the need for pivoting has nothing to do with whether the matrix is singular ornearly singular.

For example, the matrix0 1A=1 0is nonsingular yet has no LU factorization unless we interchange rows, whereas the singularmatrix1 1A=1 1does have an LU factorization.In practice, using finite-precision arithmetic, we must avoid not only zero pivots butalso small pivots in order to prevent unacceptable error growth, as shown in the followingexample. Let 1A=,1 1where is a positive number smaller than the unit roundoff mach in a given floating-pointsystem. If we do not interchange rows, then the pivot is and the resulting multiplier is−1/, so that we get the elimination matrix10M=,−1/ 1and hence 1 011L=and U ==1/ 10 1 − 1/0 −1/in floating-point arithmetic.

Характеристики

Список файлов книги