f2-2

2.2GaussianEliminationwithBacksubstitution 33Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).2.2 GaussianEliminationwithBacksubstitution The usefulness of Gaussian elimination with backsubstitution is primarily pedagogical. It stands between full elimination schemes such as Gauss-Jordan, and triangular decomposition schemes such as will be discussed in the next section. Gaussian elimination reduces a matrix not all the way to the identity matrix, but onlyhalfway,toamatrixwhosecomponentsonthediagonalandabove(say)remain nontrivial. Let us now see what advantages accrue. Suppose that in doing Gauss-Jordan elimination, as described in §2.1, we at eachstagesubtractawayrowsonly belowthethen-currentpivotelement. When a22 is the pivot element, forexample,we dividethe secondrow by its value (as before), but now use the pivot row to zero only a32and a42, not a12(see equation 2.1.1). Suppose,also,thatwedoonlypartialpivoting,neverinterchangingcolumns,sothat the 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 reduced equation that looks like this (in the case of a single right-handside vector):  a/prime 11 a/prime12 a/prime13 a/prime14 0 a/prime 22 a/prime23 a/prime24 00 a/prime 33 a/prime34 000 a/prime 44 · x1 x2 x3 x4 = b/prime 1 b/prime2 b/prime 3 b/prime 4  (2.2.1 ) Here the primes signify that the a’s and b’s do not have their original numerical values, but have been modified by all the row operations in the elimination to this point. The procedure up to this point is termed Gaussian elimination . Backsubstitution But how do we solve for the x’s? The last x(x4in this example) is already isolated, namely x4=b/prime 4/a/prime44(2.2.2 ) With the last xknown we can move to the penultimate x, x3=1 a/prime 33[b/prime 3−x4a/prime 34]( 2.2.3 ) and then proceed with the xbefore that one. The typical step is xi=1 a/prime ii b/prime i−N/summationdisplay j=i+1a/prime ijxj  (2.2.4 ) The procedure defined by equation (2.2.4) is called backsubstitution . The com- bination of Gaussian elimination and backsubstitution yields a solution to the set of equations. 34 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).TheadvantageofGaussianeliminationandbacksubstitutionoverGauss-Jordan elimination is simply that the former is faster in raw operations count: Theinnermost loops of Gauss-Jordan elimination, each containing one subtraction and one multiplication, are executed N 3and N2Mtimes (where there are Nequations and Munknowns). The corresponding loops in Gaussian elimination are executed only1 3N3times (only half the matrix is reduced, and the increasing numbers of predictable zeros reduce the count to one-third), and1 2N2Mtimes, respectively. Each backsubstitutionof a right-handside is1 2N2executionsof a similar loop (one multiplication plus one subtraction). For M/lessmuchN(only a few right-hand sides) Gaussian elimination thus has about a factor three advantage over Gauss-Jordan.(We couldreducethis advantagetoa factor1.5by notcomputingtheinversematrix as part of the Gauss-Jordan scheme.) For computing the inverse matrix (which we can view as the case of M =N right-hand sides, namely the Nunit vectors which are the columns of the identity matrix),Gaussianeliminationandbacksubstitutionatfirstglancerequire 1 3N3(matrix reduction) +1 2N3(right-hand side manipulations) +1 2N3(Nbacksubstitutions) =4 3N3loopexecutions,whichismorethanthe N3forGauss-Jordan. However,the unit vectors are quite special in containing all zeros except for one element. If this is taken intoaccount,the right-sidemanipulationscanbereducedto only1 6N3loop executions,and,for matrixinversion,thetwo methodshave identicalefficiencies. BothGaussianeliminationandGauss-Jordaneliminationsharethedisadvantage that all right-handsides must beknownin advance. The LUdecompositionmethod in the next section does not share that deficiency, and also has an equally small operations count, both for solution with any number of right-hand sides, and formatrix inversion. For this reason we will not implement the method of Gaussian elimination 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. Isaacson, 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: Addison- Wesley), §2.2.1. Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). 2.3 LU Decomposition and Its Applications Suppose we are able to write the matrix Aas a product of two matrices, L·U=A (2.3.1 ) whereLislower triangular (has elements only on the diagonal and below) and U isupper triangular (has elements only on the diagonal and above). For the case of