Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 11
Текст из файла (страница 11)
NowMultiplication bymeans to multiply row n - 2 by mi,n-2 and add it to row i fori = n - 1, n. In the special case oftimesthis clearly results in2.2 MATRIX FACTORIZATION47Repetition of this argument shows thatwhich is a lower triangular matrix that we shall call L. Finally, then, we see thatA = LU, where the L and U arise in an simple way during elimination.
Because thediagonal elements of L are all ones, we do not need to store them. The matrix L isformed a column at a time and the elements can be written in the space occupied byelements of A that are set to zero. As the scheme was described in the precedingsection, the elements of U are written over the elements of A as they are formed. Oneof the virtues of describing elimination in terms of a matrix factorization is that it isclear how to handle more than one vector b in solving Ax = LUx = b. For any givenb we first solveLy = band thenUx=y.This yields the desired x, for substitution shows thatL y = L ( U x) = (LU)x = Ax = b.Forward substitution to solve the lower triangular system Ly = b, ory 1 = b1m2,1y1 +y2 = b2is justy1 = b1y2 = b2 - m2,1y1yn = bn - mn,1y1 - mn,2y2 - ··· - mn,n-1yn-l.Back substitution is used to solve Ux = y, oru11 + u12x2 + ··· + ulnxn = ylu22x2 + ··· + u2nxn = y2unnxn = yn,48CHAPTER 2SYSTEMS OF LINEAR EQUATIONSbut now the order is xn,xn-l,.
. . ,x1 :xn = yn/un,nx n - l = (y n - 1 - u n - l , n x n )/u n - l , n - lx1 = ( y1 - u12x2 - u13x3 - ··· - ul,nxn)/u1,1.There is another important decomposition of A that arises naturally in a discussionof least squares fitting of data. The reader should turn to the advanced texts cited fordetails, but a little perspective is useful. If the Gram-Schmidt process is used to form aset of orthonormal vectors from the columns of A, a decomposition A = QR is obtained,where Q is an orthogonal matrix and R is an upper triangular matrix. An orthogonalmatrix Q is one for which Q-1 = QT, so to solve A x = QRx = b, all we have to dois form Rx = QTb and solve it by backward substitution. The Gram-Schmidt processin its classic form is not numerically stable, but there is a modification that is.
Amore popular way to obtain a QR decomposition stably is by means of Householdertransformations. Error bounds for this way of solving systems of linear equations aremuch better than those for Gaussian elimination because there is no growth factor.However, the method is about twice as expensive. In another section dealing withmatrices with special structure we take up important circumstances that favor Gaussianelimination over QR decomposition. Because the accuracy of Gaussian elimination isalmost always satisfactory in practice, it is preferred except in special circumstances.One exception is in the solution of least squares problems where the QR decompositionis especially convenient and problems are often very ill-conditioned.EXERCISES(c) Exercise 2.
1f.2.5 Find the L and U in the LU decomposition (no pivoting) for the coefficient matrices in(a) Exercise 2.1a;(b) Exercise 2.1b;2.6 Find an LU decomposition for the singular coefficientmatrix in Exercise 2.1d. Is the decomposition unique?2.3 ACCURACYThere are two main sources of error in the computed solution z of the linear systemAx = b. The data A and b may not be measured exactly, and even if they are, errors aregenerally made in representing them as floating point numbers. Roundoff errors occurin the elimination and forward/backward substitution algorithms.
It seems obviousthat we should study the errore = x - z,but it turns out that a different way of approaching the issue of accuracy is illuminating.A backward error analysis views z as the exact solution of a perturbed problem(A +A ) z = b + b.2.3 ACCURACY 49If the perturbations,andare comparable to the measurement errors or the roundoff in the entries of A or b, then it is reasonable to say that z is about as good a solutionas we might hope to get.A BACKWARD ERROR ANALYSISA floating point error analysis of a simple system of linear equations will be illuminating. Suppose that the systemu11x1 + u12x2 = b1u22x2 = b2has arisen directly or as the result of applying Gaussian elimination to a more generalsystem.
In our development of elimination we discussed how a small pivot, here u11and u22, could be dangerous both for its direct effects and because it might lead tolarge elements in the upper triangular matrix, here u12. Analyzing the error in thissimple case will help us to understand this. Backward substitution in exact arithmeticproduces the true solution asIn floating point arithmetic,Computation of the other component involves several steps. First we computethenand finallyIn a backward error analysis, we express the solutioncomputed in floatingpoint arithmetic as the solution in exact arithmetic of a perturbed problem:50CHAPTER 2SYSTEMS OF LINEAR EQUATIONSAs the notation suggests, this can be done without perturbing the right-hand side ofthe equation.
The equationwill be valid if we defineSimilarly, the equationwill be valid if we defineWith these definitions we have expressed the computed solution of the given problem as the exact solution of a problem with perturbed matrix. It is seen that none ofthe coefficients of the matrix is perturbed by more than about two units of roundoff.This analysis tells us that the backward substitution algorithm is sure to produce agood result in the sense that the computed solution is the exact solution of a problemclose to the one given. However, that is not the same as saying that the computedsolution is close to the true solution.
A forward error analysis bounds directly thedifference between the computed and true solutions.Our basic assumption about floating point arithmetic is that a single operation iscarried out with a relative error bounded by the unit of roundoff u, so we haveSubstitution of the expressions developed earlier and a little manipulation shows thatwhereThis implies that2.3 ACCURACY 51According to this bound, the relative error is generally small.
A large relative erroris possible only when |x2u12| >> |x1u11|. If the solution is such that both componentsare of comparable size, a large relative error is possible only when the pivot u11 issmall and/or the entry u12 in the upper triangular matrix is large. Large relative errorsare more likely when |x2| >> |x1|. The denominator can be written in the formxl u 11 = b 1 - x2 u 12 ,showing that the relative error can be large when the numerator is large and the denominator is small because of cancellation.ROUNDOFF ANALYSISA natural way to measure the quality of an approximate solution z is by how well itsatisfies the equation.
A virtue of this is that it is easy to compute the residualr = b - Az .In this measure, a good solution z has a small residual. Because of cancellation (seeExample 1.10), if we should want an accurate residual for a good solution, it will benecessary to compute it in higher precision arithmetic, and this may not be available.The residual provides afor the backward error analysis, namely,= -r.The residual r is connected to the error e byr = b - Az = Ax - Az = A(x - z) = Ae-1or e = A r.
A small residual r, hence a smallmay be perfectly satisfactory fromthe point of view of backward error analysis even when the corresponding error e isnot small.Example 2.8.the systemTo illustrate the distinction between the two points of view, consider(2.11)We carry out the elimination process using three-digit chopped decimal arithmetic.After the first step we haveIt then follows thatz1 = (0.200 - 0.547z 2 )/0.747 = 0.267,52CHAPTER 2SYSTEMS OF LINEAR EQUATIONSso the computed solution isThe exact solution to (2.11) is easily found to be xl = 1 and x2 = -1. Therefore theerror (in exact arithmetic) isIn contrast, the residual (in exact arithmetic) isr = b - AzThis says that z is the exact solution of Az = b +where b1 = 0.200 is perturbed to0.199449 and b2 is perturbed to 0.166341. Thus, z is the solution of a problem veryclose to the one posed, even though it differs considerably from the solution x of theoriginal problem.nThe fundamental difficulty in Example 2.8 is that the matrix in the system (2.11)is nearly singular.
In fact, the first equation is, to within roundoff error, 1.2 times thesecond. If we examine the elimination process we see that z2 was computed from twoquantities that were themselves on the order of roundoff error. Carrying more digits inour arithmetic would have produced a totally different z2. The error in z2 propagatesto an error in z1. This accounts for the computed solution being in error. Why thenare the residuals small? Regardless of z2, the number z1 was computed to make theresidual for the first equation as nearly zero as possible in the arithmetic being used.The residual for the second equation should also be small because the system is closeto singular: the first equation is approximately a multiple of the second.
In Section 2.2we observed that any matrix A could have its rows interchanged to obtain a matrixPA, which can be decomposed as the product of a lower triangular matrix L and anupper triangular matrix U. For simplicity we ignore the permutation matrix P in whatfollows. An error analysis of elimination using floating point arithmetic shows that Land U are computed with errorsandrespectively. Then A is not exactly equalto the product (L +) (U + ). Letbe defined so thatthat is,2.3 ACCURACY 53We might reasonably hope to compute L with errors M that are small relative to L, andthe same for U. However, the expression for AA shows that the sizes of L and U playimportant roles in how well A is represented by the computed factors.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
