Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 13
Текст из файла (страница 13)
The special case of a change in only one entry of b is(2.25)and(2.26)This says that the relative error of a solution component xi will be sensitive to therelative error in a data component bp whenever the factoris large. The resultsusing norms told us that “large” components in the inverse matrix indicate possiblesensitivity to errors in the data.
This result goes much further. It brings to our attentionthat “small” components in the solution vector indicate possible sensitivity to errors inthe data. More to the point, it shows which components are sensitive and how much.2.3 ACCURACY 59What if A is perturbed? For simplicity we study the effect on a solution componentxi when only a single entry apq of A is altered. In component form Ax = b isTaking the partial derivative with respect to apq leads toandIn terms of the vectorsw = (w i), where wi = 0 for ip, wp= -xq,this is a system of equations Av = w. Then v = A-1 w, or in component form,We conclude that for “small” perturbationsperturbed to xi +whereof apq, the solution component xi isIn terms of relative changes, this is(2.27)This approximate result shows much more detail than the bound (2.22).
In particular,it is clear that if there is a solution component xq that is large compared to a componentxi , then xi can be sensitive to perturbations in column q of A.ESTIMATING CONDITIONAlthough there is far more information in the equality (2.26) and the approximation (2.27) than in the condition number inequality (2.22), they require knowledgeof A-1, which is relatively expensive to compute.
An efficient algorithm for the calculation of A-1 requires roughly three times the work needed for elimination. Tocompute cond(A) exactly requires A-l, but for our purposes an estimate for ||A-1||60CHAPTER 2SYSTEMS OF LINEAR EQUATIONSwould suffice. An adequate estimate can be obtained with little more work than a forward/backward substitution. The calculation of ||A|| using the maximum norm is easy.To estimate ||A-1|| the basic observation is that if Ay = d, then from (2.16)||y|| = ||A-1 d|| < ||A-1|| ||d||.For any vector y we can form d and then obtain a lower bound for ||A-l|| from||A-1|| > ||y||/||d||Using a factorization of A, Cline, Moler, Stewart, and Wilkinson [3] construct y andd that result in a “large” ratio ||y||/||d||.
A lower bound for ||A-1|| is obtained in anycase, and often the algorithm results in a value comparable to ||A-1||. The code Factordiscussed below uses this approach to estimating cond(A). LAPACK [l] refines thisestimate by an iterative procedure that involves repeated solutions of linear systemsinvolving A.The idea behind the condition estimate explains a common way that ill-conditioning is revealed in a computation.
Let y be the computed solution of Ax = d. In theinequalitylet us approximate ||x|| by ||y||. Suppose we find that the right-hand side of the inequality is then “large.” If the computed solution is comparable to the size of the truesolution, this says that the problem is ill-conditioned.
If the computed solution is noteven comparable in size to the true solution, it is very inaccurate. Either way, a largevalue of this quantity is a warning. Often a problem is scaled naturally so that ||A||and ||d|| are about the same size, and in this common situation a large value of thequantity corresponds to a large computed solution.
With this in mind, if you shouldcompute a solution that seems “large,” you should question the conditioning of theproblem and/or the accuracy of the solution.Often a problem is scaled so that ||A|| and ||d|| are about 1. If this is the case and yshould turn out to be large, the inequality shows that the matrix must be ill-conditioned.Example 2.13. An example of Wilkinson [15] illustrates this and makes the pointthat a matrix can be very ill-conditioned without any small pivots arising in the elimination process.
Because the matrix of the set of 100 equationsis upper triangular, the pivot elements are on the diagonal, and obviously none is“small.” Back substitution shows thatxl = l/(0.600 × 0.599 × ··· × 0.502 × 0.501) > (0.6)-100 > 1022,2.4 ROUTINES FACTOR AND SOLVE61nhence this matrix is very ill-conditioned.The purpose of the condition number estimator in Factor is not to get an accuratevalue for the condition number, rather to recognize when it is “large.” The inequality(2.22) is a worst-case bound.
If the condition estimate of Factor indicates that theproblem is ill-conditioned, you might well feel justified in computing A-1 so that youcan use the more informative (2.26) and (2.27) to assess the effects of perturbations tothe data.EXERCISES(d) Use exact arithmetic to calculate residuals foreach solution. Which method seems better [comparewith part (c)]? Do the residuals indicate this?2.7 For the linear system in Exercise 1.1, letIn exact arithmetic, calculate the residuals r = b - A yand s = b - Az. Does the better approximation havethe smaller residual?2.8 Consider the systemx1 + x2 = 2(e) Using the formula (2.21), compute A-1 for thissystem, and cond(A).2.9 Assume that the computed solution to a nonsingularlinear system is(-10.4631, 0.00318429, 3.79144, -0.000422790)10x1 + 1018x2 = 10 + 1018and the condition number is 1200.Do not use more than 15-digit decimal arithmetic inthe computations of parts (a) and (b).
This will, forexample, result in 10 + 1018 becoming just 10.(a) Solve using Gaussian elimination with partial pivoting.(a) What is the uncertainty (±?) in each componentof the computed solution? Assume exact data and aunit roundoff of 10-6.(b) Repeat (a) but for the case whenandare each 10-5.(b) Divide each row by its largest |aij| and then useGaussian elimination with partial pivoting.2.10 On a machine with unit roundoff 10-l7 with b exactbuthow large a condition num(c) Solve by hand using any method and exact arithber can we tolerate if we wantmetic.2.4 ROUTINES FACTOR AND SOLVEIn this section we describe routines Factor and Solve that can be used to solve thesystem Ax = b.
Routine Factor performs Gaussian elimination with partial pivotingon the matrix A and estimates its condition number. It saves the multipliers and thepivot information to be used by Solve to obtain the solution x for any right-hand sideb. A typical call to Factor in FORTRAN isCALL FACTOR (A,MAXROW,NEQ,COND,PVTIDX,FLAG,TEMP)while a typical function evaluation of Factor isflag = Factor(a, neq, cond, pivot_index);62CHAPTER 2SYSTEMS OF LINEAR EQUATIONSin the C++ version. The parameter cond must be passed by reference since its value isset by the function Factor; in C the address of cond must be explicitly passed so thatits call looks likeflag = Factor(a, neq, &cond, pivot_index);Input variables in FORTRAN are A, the array containing the coefficient matrix A;MAXROW, the declared row dimension of A in the program that calls Factor; andNEQ, the number of equations to be solved.
Output variables are A, containing theupper triangular matrix U in positions aij, i < j, and the lower triangular matrix L inpositions aij, i > j (as long as FLAG is zero); FLAG, an integer variable that indicateswhether or not zero pivots were encountered; COND, an estimate of the conditionnumber of A in the maximum norm; PVTIDX, an array that records the row interchanges; and in FORTRAN 77 we need TEMP, an array used for temporary storage.In the C and C++ versions corresponding variables are a for A, neq for NEQ,cond for COND, and pivotindex for PVTIDX.
Note that in the C and C++ versions,(1) instead of declaring the matrix A to have two indices, we use a pointer to a vectorconsisting of the rows of A; (2) there is no need to reserve space for the temporaryindex TEMP because this allocation can be made dynamically when needed (as can bedone for a and for pivotindex); (3) the output flag is the return variable for the functionFactor.
Because arrays are typically indexed starting with zero in C and C++, we havemodified the algorithm accordingly. Some further comments are in order about thevariables for FactorMAXROW and NEQ (or neq). If, for example, the array A is declaredas A(10,10) in the calling program and we are solving three equationsin three unknowns, then MAXROW would have the value 10 in theFORTRAN version and NEQ (or neq in the C version) the value 3.COND.
COND is a lower bound for cond(A) in the maximum normand is often a good approximation to cond(A). If fl(COND + 1) =COND, the matrix A is singular to working precision. Because weare working with arithmetic of limited precision, exact singularity isdifficult to detect. In particular, the occurrence of a zero pivot does notnecessarily mean that the matrix is singular nor does a singular matrixnecessarily produce a zero pivot (see Exercise 2.16). When FLAG isnonzero, the output COND is meaningless.PVTIDX (or pivotindex).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
