Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 16

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

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

1969, Applied Numerical Methods (New York:Wiley), Example 5.2, p. 282.Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York:McGraw-Hill), Program B-2, p. 298.Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New York: Wiley).Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:McGraw-Hill), §9.3–1.2.2 Gaussian Elimination with BacksubstitutionThe usefulness of Gaussian elimination with backsubstitution is primarilypedagogical.

It stands between full elimination schemes such as Gauss-Jordan, andtriangular decomposition schemes such as will be discussed in the next section.Gaussian elimination reduces a matrix not all the way to the identity matrix, butonly halfway, to a matrix whose components on the diagonal and above (say) remainnontrivial. Let us now see what advantages accrue.Suppose that in doing Gauss-Jordan elimination, as described in §2.1, we ateach stage subtract away rows only below the then-current pivot element. When a22is the pivot element, for example, we divide the second row by its value (as before),but now use the pivot row to zero only a32 and a42 , not a12 (see equation 2.1.1).Suppose, also, that we do only partial pivoting, never interchanging columns, so thatthe order of the unknowns never needs to be modified.Then, when we have done this for all the pivots, we will be left with a reducedequation that looks like this (in the case of a single right-hand side vector):     x1b1a11 a12 a13 a14 0 a22 a23 a24   x2   b2 (2.2.1)·  =  00 a33 a34x3b3000 a44x4b4Here the primes signify that the a’s and b’s do not have their original numericalvalues, but have been modified by all the row operations in the elimination to thispoint.

The procedure up to this point is termed Gaussian elimination.42Chapter 2.Solution of Linear Algebraic EquationsBacksubstitutionBut how do we solve for the x’s? The last x (x4 in this example) is alreadyisolated, namelyx4 = b4 /a44(2.2.2)With the last x known we can move to the penultimate x,x3 =1 [b − x4 a34 ]a33 3and then proceed with the x before that one. The typical step isN1aij xj xi = bi −aii(2.2.3)(2.2.4)j=i+1The procedure defined by equation (2.2.4) is called backsubstitution. The combination of Gaussian elimination and backsubstitution yields a solution to the setof equations.The advantage of Gaussian elimination and backsubstitution over Gauss-Jordanelimination is simply that the former is faster in raw operations count: Theinnermost loops of Gauss-Jordan elimination, each containing one subtraction andone multiplication, are executed N 3 and N 2 M times (where there are N equationsand M unknowns).

The corresponding loops in Gaussian elimination are executedonly 13 N 3 times (only half the matrix is reduced, and the increasing numbers ofpredictable zeros reduce the count to one-third), and 12 N 2 M times, respectively.Each backsubstitution of a right-hand side is 12 N 2 executions of a similar loop (onemultiplication plus one subtraction). For M N (only a few right-hand sides)Gaussian elimination thus has about a factor three advantage over Gauss-Jordan.(We could reduce this advantage to a factor 1.5 by not computing the inverse matrixas part of the Gauss-Jordan scheme.)For computing the inverse matrix (which we can view as the case of M = Nright-hand sides, namely the N unit vectors which are the columns of the identitymatrix), Gaussian elimination and backsubstitution at first glance require 13 N 3 (matrixreduction) + 12 N 3 (right-hand side manipulations) + 12 N 3 (N backsubstitutions)= 43 N 3 loop executions, which is more than the N 3 for Gauss-Jordan.

However, theunit vectors are quite special in containing all zeros except for one element. If thisis taken into account, the right-side manipulations can be reduced to only 16 N 3 loopexecutions, and, for matrix inversion, the two methods have identical efficiencies.Both Gaussian elimination and Gauss-Jordan elimination share the disadvantagethat all right-hand sides must be known in advance.

The LU decomposition methodin the next section does not share that deficiency, and also has an equally smalloperations count, both for solution with any number of right-hand sides, and formatrix inversion. For this reason we will not implement the method of Gaussianelimination as a routine.CITED REFERENCES AND FURTHER READING:Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:McGraw-Hill), §9.3–1.432.3 LU Decomposition and Its ApplicationsIsaacson, E., and Keller, H.B.

1966, Analysis of Numerical Methods (New York: Wiley), §2.1.Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), §2.2.1.Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New York: Wiley).2.3 LU Decomposition and Its ApplicationsSuppose we are able to write the matrix A as a product of two matrices,L·U=A(2.3.1)where L is lower triangular (has elements only on the diagonal and below) and Uis upper triangular (has elements only on the diagonal and above). For the case ofa 4 × 4 matrix A, for example, equation (2.3.1) would look like this:  α11 α21α31α410α22α32α4200α33α430β110   0· 00α440β12β2200β13β23β330β14β24 β34β44a11a=  21a31a41a12a22a32a42a13a23a33a43a14a24 a34a44(2.3.2)We can use a decomposition such as (2.3.1) to solve the linear setA · x = (L · U) · x = L · (U · x) = b(2.3.3)by first solving for the vector y such thatL·y=b(2.3.4)U·x=y(2.3.5)and then solvingWhat is the advantage of breaking up one linear set into two successive ones?The advantage is that the solution of a triangular set of equations is quite trivial, aswe have already seen in §2.2 (equation 2.2.4).

Thus, equation (2.3.4) can be solvedby forward substitution as follows,y1 =yi =b1α11i−11 αij yj bi −αiij=1(2.3.6)i = 2, 3, . . . , Nwhile (2.3.5) can then be solved by backsubstitution exactly as in equations (2.2.2)–(2.2.4),yNxN =βNNN(2.3.7)1 xi =βij xj yi −i = N − 1, N − 2, .

. . , 1βiij=i+144Chapter 2.Solution of Linear Algebraic EquationsEquations (2.3.6) and (2.3.7) total (for each right-hand side b) N 2 executionsof an inner loop containing one multiply and one add. If we have N right-handsides which are the unit column vectors (which is the case when we are inverting amatrix), then taking into account the leading zeros reduces the total execution countof (2.3.6) from 12 N 3 to 16 N 3 , while (2.3.7) is unchanged at 12 N 3 .Notice that, once we have the LU decomposition of A, we can solve with asmany right-hand sides as we then care to, one at a time. This is a distinct advantageover the methods of §2.1 and §2.2.Performing the LU DecompositionHow then can we solve for L and U, given A? First, we write out thei, jth component of equation (2.3.1) or (2.3.2). That component always is a sumbeginning withαi1 β1j + · · · = aijThe number of terms in the sum depends, however, on whether i or j is the smallernumber.

We have, in fact, the three cases,i<j:i=j:i>j:αi1 β1j + αi2 β2j + · · · + αii βij = aijαi1 β1j + αi2 β2j + · · · + αii βjj = aijαi1 β1j + αi2 β2j + · · · + αij βjj = aij(2.3.8)(2.3.9)(2.3.10)Equations (2.3.8)–(2.3.10) total N 2 equations for the N 2 + N unknown α’s andβ’s (the diagonal being represented twice). Since the number of unknowns is greaterthan the number of equations, we are invited to specify N of the unknowns arbitrarilyand then try to solve for the others.

In fact, as we shall see, it is always possible to takeαii ≡ 1i = 1, . . . , N(2.3.11)A surprising procedure, now, is Crout’s algorithm, which quite trivially solvesthe set of N 2 + N equations (2.3.8)–(2.3.11) for all the α’s and β’s by just arrangingthe equations in a certain order! That order is as follows:• Set αii = 1, i = 1, . . . , N (equation 2.3.11).• For each j = 1, 2, 3, . .

. , N do these two procedures: First, for i =1, 2, . . ., j, use (2.3.8), (2.3.9), and (2.3.11) to solve for βij , namelyβij = aij −i−1αik βkj .(2.3.12)k=1(When i = 1 in 2.3.12 the summation term is taken to mean zero.) Second,for i = j + 1, j + 2, . . . , N use (2.3.10) to solve for αij , namely1αij =βjjaij −j−1αik βkj.k=1Be sure to do both procedures before going on to the next j.(2.3.13)452.3 LU Decomposition and Its Applicationsacegietc.xbddifaghjonalsubdelemeniatsgonaxleleetc.mentsFigure 2.3.1.

Crout’s algorithm for LU decomposition of a matrix. Elements of the original matrix aremodified in the order indicated by lower case letters: a, b, c, etc. Shaded boxes show the previouslymodified elements that are used in modifying two typical elements, each indicated by an “x”.If you work through a few iterations of the above procedure, you will see thatthe α’s and β’s that occur on the right-hand side of equations (2.3.12) and (2.3.13)are already determined by the time they are needed.

You will also see that every aijis used only once and never again. This means that the corresponding αij or βij canbe stored in the location that the a used to occupy: the decomposition is “in place.”[The diagonal unity elements αii (equation 2.3.11) are not stored at all.] In brief,Crout’s method fills in the combined matrix of α’s and β’s,β11 α21α31α41β12β22α32α42β13β23β33α43β14β24 β34β44(2.3.14)by columns from left to right, and within each column from top to bottom (seeFigure 2.3.1).What about pivoting? Pivoting (i.e., selection of a salubrious pivot elementfor the division in equation 2.3.13) is absolutely essential for the stability of Crout’s46Chapter 2.Solution of Linear Algebraic Equationsmethod.

Only partial pivoting (interchange of rows) can be implemented efficiently.However this is enough to make the method stable. This means, incidentally, thatwe don’t actually decompose the matrix A into LU form, but rather we decomposea rowwise permutation of A. (If we keep track of what that permutation is, thisdecomposition is just as useful as the original one would have been.)Pivoting is slightly subtle in Crout’s algorithm. The key point to notice is thatequation (2.3.12) in the case of i = j (its final application) is exactly the same asequation (2.3.13) except for the division in the latter equation; in both cases theupper limit of the sum is k = j − 1 (= i − 1).

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

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