Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 10
Текст из файла (страница 10)
Wilkinson [15] points out that there isequality in this bound for matrices of the formHowever, usually the growth is very modest. Research into this matter is surveyedin [11]. There are other ways of selecting pivot elements that lead to better boundson the error and there are other ways of solving linear systems that have still betterbounds. Some details will be mentioned later, but in practice the numerical propertiesof Gaussian elimination with partial pivoting are so good that it is the method of choiceexcept in special circumstances, and when one speaks of “Gaussian elimination” it isassumed that partial pivoting is done unless something to the contrary is said.
Gaussianelimination with partial pivoting is the basic method used in the popular computingenvironments MATLAB, Mathematica, and MATHCAD.40CHAPTER 2SYSTEMS OF LINEAR EQUATIONSExample 2.7. Using exact arithmetic and elimination with partial pivoting, solve thefollowing system:Since |2| > |1|, we interchange the first and second equations to getUsing 2 as a pivot, we eliminate the coefficients of xl in equations two and three to getSince |-2| > |l|, we interchange equations two and three,and using -2 as pivot obtainBack substitution then givesorThe algorithm for elimination is quite compact:Elimination, modification of b.for k = l,...,n - 1 begininterchange rows so that |akk| = maxk<i<n |aik|if |akk| = 0, set singularity indicator, returnfor i = k + l, .
. . , n begin2.1 GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING41t := a ik /akkfor j = k + l, . . . , n beginaij := aij - t * akjend jbi := bi - t * bkend iend kif |ann| = 0, set singularity indicator.Back substitutionfor i = n, . . . , 1 beginxi := bifor j = i + l, .
. . ,n beginxi := xi - aij * xjend jxi := xi /aiiend i.Sometimes we are interested in solving problems involving one matrix A and severalright-hand sides b. Examples are given in the exercises of problems with the righthand sides corresponding to different data sets. Also, if we should want to compute theinverse of an n × n matrix A, this can be done a column at a time. It is left as an exerciseto show that column i of A-1 is the result of solving the system of equations withcolumn i of the identity matrix as b. If we know all the right-hand sides in advance, itis clear that we can process them simultaneously. It is not always the case that they areall known in advance. The residual correction process we take up later is an example.For such problems it is important to observe that if we save the multipliers and recordhow the rows are interchanged when processing A, we can process b separately.
Tounderstand how this can be valuable, we first need to look at the costs of the variousportions of this algorithm.As a measure of work we count arithmetic operations. Since the number of additions and subtractions equals the number of multiplications, only the latter (as well asdivisions) are counted. It is easy enough to see that elimination requiresn( n - 1) (2n - 1)/6 multiplications and n( n - 1)/2 divisions.Modification of b requiresn( n - 1)/2 multiplications.Back substitution requiresn( n - 1)/2 multiplications and n divisions.For large n the multiplications dominate the cost, both because there are more of themand because they are relatively expensive. The most important point is that processingthe matrix A is the bulk of the cost of solving a system of linear equations of evenmoderate size.42CHAPTER 2SYSTEMS OF LINEAR EQUATIONSSeveral designs are seen in popular codes.
The most straightforward is to inputA and b and have the code compute the solution x and return it. It is quite easy tomodify such a code to accept input of m right-hand sides, process all the right-handsides along with A, and return all the solutions in an array. This is considerably cheaperthan solving the problems one after another because A is processed only once and thisis the most expensive part of the computation.
In detail, solving m systems with thesame A simultaneously requiresmultiplications. Solving them independently requiresmultiplications. If, for example, we wished to invert A, we would have m = n and thecost would bemultiplications. For large n, there is a considerable difference.
The most flexible design separates the two phases of the computation. By saving the information necessaryfor processing the right-hand side, systems involving the same matrix A can be solvedindependently and just as inexpensively as if all the right-hand sides were available tobegin with. This is the design found in production-grade software and in the programsof this chapter. Because it is a little more trouble to use than the simplest design, itis not unusual for libraries to have both. The computing environment MATLAB is anexample of this.EXERCISES2.1 Using elimination with partial pivoting, determine which of the following systems are singular and which are nonsingular.
For the nonsingular problems, findsolutions. Use exact arithmetic.(a)(b)2.1 GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING43(c)xl2 x13 x1+++2x24x26x2- x3+ x3- 2x3===277xl2 x1x1++x2x22x2+-x3x32x3===0-3-2xl2 x12x1++x2x2+ x3- x3- 4x3===0-3-6(d)(e)(f)2 x1 - 3x2 + 2x3 + 5x4 = 3x1 - x2 + x3 + 2x4 = 13 x1 + 2x2 + 2x3 + x4 = 0xl + x2 - 3x3 - x4 = 02.2 Four loads applied to a three-legged table yield the following system for the reactions on the legs:R l + R 2 + R 3 = 110.00R1+ R2 = 78.33R2 + R3 = 58.33.Solve for R1, R2, and R3 by hand.2.3 The following set of equations arises in analyzing loads on an A-frame:8.00RE - 1784.00 = 0.00-8.00RA + 1416.00 = 0.00Ch + Dh = 0.00Cv + D v + 223.00 = 0.00-5.18Cv - 5.18Ch + 446.00 = 0.00-5.77 Dv - 1456.00 = 0.00-5.77Bv - 852.00 = 0.00Bh + Dh = 0.00.Solve the equations by hand.2.4 Consider the linear systemxl + 1/2x2+1/3x3=11/2x1 + 1/3x2 + 1/4x3= 01/3x1 + 1/4x2 + 1/5x3= 0.44CHAPTER 2SYSTEMS OF LINEAR EQUATIONS(a) Solve the system using exact arithmetic (any method).(b) Put the system in matrix form using a two-digit decimal chopped representation.(c) Solve the system in (b) without partial pivoting [same arithmetic as (b)].(d) Solve the system in (b) with partial pivoting [same arithmetic as (b)].(e) Solve the system in (b) using exact arithmetic.2.2 MATRIX FACTORIZATIONBecause it is illuminating, more advanced books invariably study elimination by viewing it as a matrix factorization.
In this section we shall see that if no pivoting is done,the elimination algorithm factors, or decomposes, the matrix A into the product LU ofa lower triangular matrix L =where= 0 if i < j,and an upper triangular matrix U = (u ij ), whereuij = 0 if i > j.When partial pivoting is done, it is a version of A with its rows interchanged that isdecomposed. Rows can be interchanged in A by multiplication with a matrix P calleda permutation matrix. It is easy to construct P. If PA is to be the result of interchangingsome rows of A, all we need do is take P to be the result of interchanging these rowsin the identity matrix I. For example, to interchange rows 1 and 3 of the 3 × 3 matrixA in PA, we useThe entire elimination process with partial pivoting can be written as theLU factorizationPA = LU.Rather than sort out the permutations, we concentrate here on the factorization withoutpivoting to show how the “LU decomposition” arises. The remainder of this sectionprovides the details of this factorization and it may be skipped by the reader unfamiliar with linear algebra.
Looking back at the elimination described in the precedingsection, we see that ifwe could multiply row 1 byand subtract it from row i to eliminate the first unknown from row i. This is done for rows2.2 MATRIX FACTORIZATION45i = 2,3,. . . ,n. It is easily verified that when the matrixmultiplies any matrix A on the left, it has the effect of multiplying row 1 of A by mil andsubtracting the result from row i of A for i = 2,. . . ,n. As with permutation matrices,this kind of matrix is found by performing the operations on the identity matrix. Withthe multipliers mil chosen as specified, the product M1A has the formFor later use we note that multiplication by the inverse of a matrix “undoes” amultiplication by the matrix.
To “undo” multiplication by Ml, we need to multiplyrow 1 by mil and add the result to row i for i = 2,. . . , n. In this way, we see thatIt is also easy to verify directly that M1-1M1 = I. Suppose we have formedIfwe want to multiply row k of this matrix by mik =and subtract itfrom row i of the matrix for i = k + l, . . . , n. This is done by multiplying the matrix by46CHAPTER 2SYSTEMS OF LINEAR EQUATIONSThenandElimination without pivoting results after n - 1 steps inwhich is an upper triangular matrix that we shall call U. Multiplication of this equationon the left bythenresults inEarlier we saw the simple form of these inverse matrices. It is a delightful fact thattheir product is also extremely simple.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
