Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 12
Текст из файла (страница 12)
Partial pivotingkeeps the elements of L less than or equal to 1 in magnitude. We also saw in (2.10) thatwas moderated with partial pivoting. In particular,the size of elements of U, thethey cannot exceed 2n-1 maxij |aij| for an n × n matrix. It can be shown rigorously, ontaking into account the errors of decomposition and of forward/backward substitution,that the computed solution z of Ax = b satisfies(A +)z = b,(2.12)where the entries ofare usually small.
To make precise how small these entriesare, we need a way of measuring the sizes of vectors and matrices. One way to measure the size of a vector x of n components is by its norm, which is denoted by ||x||.Several definitions of norm are common in numerical analysis. One that is likely to beAll vector norms possess many offamiliar is the Euclidean length of x,the properties of length.
The norm used in this chapter is the maximum norm(2.13)If A is an n × n matrix and x is an n -vector, then Ax is also an n-vector. A matrix normcan be defined in terms of a vector norm by(2.14)Geometrically, this says that ||A|| is the maximum relative distortion that the matrix Acreates when it multiplies a vector x 0. It is not easy to evaluate ||A|| directly from(2.14), but it can be shown that for the maximum norm (2.13)(2.15)which is easy enough to evaluate.
An important inequality connects norms of vectorsand matrices:||Ax|| < ||A|| ||x||.(2.16)For x0 this follows immediately from the definition (2.14). For x = 0 we note thatA x = 0 and that ||x|| = 0, from which the inequality is seen to hold.Example 2.9. Let x =||x|| = max [| - 1|, |2|, |3|] = 3.54CHAPTER 2SYSTEMS OF LINEAR EQUATIONSLet A =||A|| = [(|l| + |-l| + |0|), (|2|+ |-2| + |3|), (|-4| + |1| + |-1|)]= max[(2), (7), (6)] = 7.Returning to the roundoff analysis for Gaussian elimination, it can be shown rigorously [11] that the computed solution z satisfies the perturbed equation (2.12) where(2.17)As usual, u is the unit roundoff.
The factor γn depends on n and can grow as fastas 2n-1. To put this in perspective, suppose that AA arises from rounding A to formmachine numbers. Thencould be as large as u|aij| andcould be as largeasAccording to the bounds, the perturbations due to the decomposition and forward/backward substitution process are at worst a factor of γn times the error made in the initialrounding of the entries of A. If the rigorous bound 2 n-1 on γn truly reflected practice, we would have to resort to another algorithm for large n.
Fortunately, for mostproblems γn is more like 10, independent of the size of n.From this it can be concluded that Gaussian elimination practically always produces a solution z that is the exact solution of a problem close to the one posed. SinceAz - b =the residual r satisfiesThis says that the size of the residual is nearly always small relative to the sizes of Aand z. However, recall that this does not imply that the actual error e is small.For additional insight as to why Gaussian elimination tends to produce solutionswith small residuals, think of the LU factorization of A discussed in Section 2.2. Theforward substitution process used to solve the lower triangular system Ly = b successively computes y1, y2,.
. . , yn so as to make the residual zero. For example, regardlessof the errors in y1 and m2,1 the value of y2 is computed so thatm2,lyl + y2 = b2,that is, the residual of this equation is zero (in exact arithmetic) with this value of y2.The same thing happens in the back substitution process to compute xn,xn-1,. . .
,x1that satisfy Ux = y. Thus, the very nature of the process responds to errors in thedata in such a way as to yield a small residual. This is not at all true when x is-1computed by first calculating the inverse A and then forming A-l b. With a littleextra work it is possible to make Gaussian elimination stable in a very strong sense.2.3 ACCURACY 55Suppose that we have solved Ax = b to obtain an approximate solution z. We canexpect it to have some accuracy, although perhaps not all the accuracy possible in theprecision used. A little manipulation shows that the error e = x - z satisfies A e = r,where r is the residual of the approximate solution z.
We have seen that if Gaussianelimination is organized properly, it is inexpensive to solve this additional system ofequations. Of course, we do not expect to solve it exactly either, but we do expectthat the computed approximation d to the error in z will have some accuracy. If itdoes, w = z + d will approximate x better than z does. In principle this process, callediterative refinement, can be repeated to obtain an approximation to x correct in all itsdigits. The trouble in practice is that for the process to work as described, we have tohave an accurate residual, and the better the approximate solution, the more difficultthis is to obtain. Skeel [14] has shown that just one step of iterative refinement withthe residual computed in the working precision will provide a computed solution thatis very satisfactory.
This solution will have a small residual and will satisfy exactlya system of equations with each coefficient differing slightly from that of the givensystem. This is much better than the result for z that states that the perturbation ina coefficient is small compared to the norm of the whole matrix, not that it is smallcompared to the coefficient itself. So, if we are concerned about the reliability ofGaussian elimination with partial pivoting, we could save copies of the matrix andright-hand side and perform one step of iterative refinement in the working precisionto correct the result as necessary.NORM BOUNDS FOR THE ERRORIn the preceding subsection, we found that roundoff errors in the algorithm could beconsidered equivalent to errors in the data A. We now study the effect of such perturbations, as well as errors in the given data, on the error e.
For simplicity, let us firstconsider the case where only the data b is in error. Let x +be the solution ofMultiply this by A-1 and use the fact that x = A-l b to get(2.18)Norm inequalities say that(2.19)But b = Ax similarly implies ||b|| < ||A|| ||x||, hence(2.20)Inequality (2.19) says that in an absolute sense, input errorscan be amplified byas much as ||A-1|| in the solution. In contrast, (2.20) says that in a relative sense, inputerrors can be magnified by as much as ||A|| ||A-1||. The important quantity ||A|| ||A-1||denoted by cond(A) is called the56CHAPTER 2 SYSTEMS OF LINEAR EQUATIONScondition number of Acond(A) = ||A|| ||A-l||.A theorem that helps us understand the condition number isIn words this says that there is a singular matrix S that differs from A in a relativesense by the reciprocal of the condition number of A.
Put differently, if A has a “large”condition number, it is “close” to a singular matrix.The condition number clearly depends on the norm, but in this book we consideronly the maximum norm (2. 15).Example 2.10. For A, find ||A||,||A-l||, cond(A ). First,||A|| = max[(|l|+|2|),(|3|+|4|)] = max[(3),(7)] = 7.For 2 x 2 matrices, the inverse matrix is easy to work out. If A =(2.21)So, in our caseandThenExample 2.11.The matrixis much closer to being singular than the matrix in Example 2.10 since2.3 ACCURACY57and||A|| = 2, ||A-l|| = 2 x 105, cond(A) = 4 x 105.The theorem about the condition number says that there is asingular matrix that iswithin 1/cond(A) = 2.5 x 10–6 of A. Although not quite this close, the simple matrixis obviously singular andThe effect of perturbing A is more complicated because it is possible that A + ∆ A issingular.
However, if the perturbation is sufficiently small, say ||A–1|| ||∆ A|| < 1, thenit can be shown [9] that A + ∆ A is nonsingular and further that if (A + ∆A) (x + ∆x) =b + ∆ b, then we have the so-calledcondition number inequality(2.22)valid for ||A–l || ||∆ A|| < 1.Inequality (2. 17) and the related discussion say that rounding errors in the course ofGaussian elimination are equivalent to solving a perturbed system for which we usually have(2.23)where u is the unit roundoff. In some applications data errors maybe much larger thanthis, and they must be used in (2.22).Example 2.12.obtainSuppose we solve Ax = b on a machine with u = 5 x 10–11 andThen from (2.22), assuming exact data so that ||∆A||/||A||the relative error is5 x 10-10, the bound on58CHAPTER 2SYSTEMS OF LINEAR EQUATIONSOn the other hand, if the data are known to be in error, saythenWith ||x||||z||and18.6 the absolute error bound is 0.37, so this analysis says thatxl = 6.23 ± 0.37x2 = 18.62 ± 0.37.One criticism of the analysis based on norms is that it is a worst-case analysis.The condition number does not depend on b, so the inequality (2.22) must allow forthe worst choice of b andA large condition number is cause for concern, but itdoes not mean a given problem will be solved inaccurately.
Furthermore, the fact thatis small does not mean that for each i,is small.OTHER ERROR EXPRESSIONS AND APPROXIMATIONSA better understanding of the size of the error can be obtained with a more carefulanalysis. Again, it is simpler to consider first only changes in b. If the entries of A-1thenin component form isare denoted by(2.24)Hence, foris the exact formula for the relative change in xi as a function of the relative changesin the bp.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
