f2-4

42 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).and then solved with ludcmpandlubksbin their present forms. This scheme is a factor of 2 inefficient in storage, since AandCare stored twice. It is also a factor of 2 inefficient in time, since the complex multiplies in a complexified version of the routines would eachuse 4 real multiplies, while the solution of a 2N×2Nproblem involves 8 times the work of anN×None. If you can tolerate these factor-of-two inefficiencies, then equation (2.3.18) is an easy way to proceed. CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), Chapter 4. Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.). Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §3.3, and p. 50. Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Engle- wood Cliffs, NJ: Prentice-Hall), Chapters 9, 16, and 18. Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §4.2. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §9.11. Horn,R.A.,andJohnson,C.R.1985, MatrixAnalysis (Cambridge:CambridgeUniversityPress). 2.4 Tridiagonal and Band Diagonal Systems of Equations The special case of a system of linear equations that is tridiagonal , that is, has nonzeroelementsonlyonthediagonalplusorminusonecolumn,is onethat occursfrequently. Alsocommonaresystemsthatare banddiagonal ,withnonzeroelements onlyalonga fewdiagonallines adjacent to the maindiagonal(aboveandbelow). For tridiagonal sets, the procedures of LUdecomposition, forward- and back- substitutioneachtakeonly O(N)operations,andthewholesolutioncanbeencoded veryconcisely. Theresultingroutine tridagisonethatwewilluseinlaterchapters. Naturally,one does not reservestorage forthe full N×Nmatrix, but only for thenonzerocomponents,storedasthreevectors. Thesetofequationstobesolvedis  b 1c1 0··· a2b2c2······ ··· a N−1bN−1cN−1 ··· 0 aN bN · u 1 u2 ··· uN−1 uN = r 1 r2 ··· rN−1 rN (2.4.1 ) Noticethat a 1andcNareundefinedandarenotreferencedbytheroutinethatfollows. 2.4TridiagonalandBandDiagonalSystems ofEquations 43Sample 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).SUBROUTINE tridag(a,b,c,r,u,n) INTEGER n,NMAX REAL a(n),b(n),c(n),r(n),u(n) PARAMETER (NMAX=500) Solves for a vector u(1:n) of length nthe tridiagonal linear set given by equation (2.4.1). a(1:n) ,b(1:n) ,c(1:n) ,a n d r(1:n) are input vectors and are not modified. Parameter: NMAX is the maximum expected value of n. INTEGER jREAL bet,gam(NMAX) One vector of workspace, gamis needed. if(b(1).eq.0.)pause ’tridag: rewrite equations’ If this happens then you should rewrite your equations as a set of order N−1,w i t h u 2 trivially eliminated. bet=b(1) u(1)=r(1)/bet do11j=2,n Decomposition and forward substitution. gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.)pause ’tridag failed’ Algorithm fails; see below. u(j)=(r(j)-a(j)*u(j-1))/bet enddo 11 do12j=n-1,1,-1 Backsubstitution. u(j)=u(j)-gam(j+1)*u(j+1) enddo 12 returnEND There is no pivoting in tridag. It is for this reason that tridagcan fail (pause) even when the underlying matrix is nonsingular: A zero pivot can be encounteredevenfora nonsingularmatrix. In practice,this is notsomethingto lose sleep about. The kinds of problems that lead to tridiagonal linear sets usually have additional properties which guarantee that the algorithm in tridagwill succeed. For example, if |bj|>|aj|+|cj| j=1,...,N (2.4.2 ) (calleddiagonaldominance )thenitcanbeshownthatthealgorithmcannotencounter a zero pivot. It is possible to construct special examples in which the lack of pivoting in the algorithmcausesnumericalinstability. Inpractice,however,suchinstabilityisalmost neverencountered— unlikethegeneralmatrixproblemwherepivotingis essential. The tridiagonal algorithm is the rare case of an algorithm that, in practice, is more robust than theory says it should be. Of course, should you ever encounter aproblem for which tridagfails, you can instead use the more general method for band diagonal systems, now described (routines bandecandbanbks). Some other matrix forms consisting of tridiagonal with a small number of additional elements (e.g., upper right and lower left corners) also allow rapid solution; see §2.7. Band DiagonalSystems Where tridiagonal systems have nonzero elements only on the diagonal plus or minus one,banddiagonalsystemsareslightlymoregeneralandhave(say) m1≥0nonzeroelements immediately tothe leftof (below) thediagonal and m2≥0nonzero elements immediatelyto itsright(above it). Ofcourse, thisisonly ausefulclassification if m1and m2areboth /lessmuchN. In that case, the solution of the linear system by LUdecomposition can be accomplished much faster, and in much less storage, than for the general N×Ncase. 44 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).The precise definition of a band diagonal matrix with elements aijis that aij=0when j>i +m2or i>j +m1 (2.4.3 ) Banddiagonalmatricesarestoredandmanipulatedinaso-calledcompactform,whichresults if the matrix is tilted 45◦clockwise, so that its nonzero elements lie in a long, narrow matrix with m1+1+ m2columns and Nrows. This is best illustrated by an example: The band diagonal matrix  3100000 415000092650000358900007932000038460000244 (2.4.4 ) which has N=7,m 1=2, and m2=1, is stored compactly as the 7×4matrix,  xx 31 x 415 9265358979323846244 x (2.4.5 ) Here xdenotes elements that are wasted space in the compact format; these will not be referenced by any manipulations and can have arbitrary values. Notice that the diagonalof the original matrix appears in column m 1+1, with subdiagonal elements to its left, superdiagonal elements to its right. The simplest manipulation of a band diagonal matrix, stored compactly, is to multiply it by a vector to its right. Although this is algorithmically trivial, you might want to studythe following routine carefully, as an example of how to pull nonzero elements a ijout of the compact storage format in an orderly fashion. Notice that, as always, the logical and physicaldimensions of a two-dimensional array can be different. Our convention is to pass N,m 1, m2, and the physicaldimensions np≥Nandmp≥m1+1+ m2. SUBROUTINE banmul(a,n,m1,m2,np,mp,x,b) INTEGER m1,m2,mp,n,np REAL a(np,mp),b(n),x(n) Matrix multiply b=A·x,w h e r e Ais band diagonal with m1rows below the diagonal andm2rows above. The input vector xand output vector bare stored as x(1:n) and b(1:n) , respectively. The array a(1:n,1:m1+m2+1) storesAas follows:The diagonal elements are in a(1:n,m1+1) . Subdiagonal elements are in a(j:n,1:m1) (with j> 1 appropriate to the number of elements on each subdiagonal). Superdiagonal elements arein a(1: j,m1+2:m1+m2+1) with j< nappropriate to the number of elements on each superdiagonal. INTEGER i,j,kdo 12i=1,n b(i)=0. k=i-m1-1 do11j=max(1,1-k),min(m1+m2+1,n-k) b(i)=b(i)+a(i,j)*x(j+k) enddo 11 enddo 12 returnEND 2.4TridiagonalandBandDiagonalSystems ofEquations 45Sample 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).It is not possible to store the LUdecomposition of a band diagonal matrix Aquite as compactly as the compact form of Aitself. The decomposition (essentially by Crout’s method, see §2.3)produces additional nonzero “fill-ins.” Onestraightforward storagescheme is to return the upper triangular factor ( U) in the same space that Apreviously occupied, and to return the lower triangular factor ( L) in a separate compact matrix of size N×m1. The diagonal elements of U(whose product, times d=±1, gives the determinant) are returned in the first column of A’s storage space. The following routine, bandec, is the band-diagonal analog of ludcmpin§2.3: SUBROUTINE bandec(a,n,m1,m2,np,mp,al,mpl,indx,d) INTEGER m1,m2,mp,mpl,n,np,indx(n)REAL d,a(np,mp),al(np,mpl),TINY PARAMETER (TINY=1.e-20) Given an n×nband diagonal matrix Awith m1subdiagonal rows and m2superdiagonal rows, compactly stored in the array a(1:n,1:m1+m2+1) as described in the comment for routine banmul , this routine constructs an LUdecomposition of a rowwise permutation ofA. The upper triangular matrix replaces a, while the lower triangular matrix is returned inal(1:n,1:m1) .indx(1:n) is an output vector which records the row permutation effected by the partial pivoting; dis output as ±1depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with banbks to solve band-diagonal sets of equations. INTEGER i,j,k,l,mmREAL dum mm=m1+m2+1 if(mm.gt.mp.or.m1.gt.mpl.or.n.gt.np) pause ’bad args in bandec’l=m1do 13i=1,m1 Rearrange the storage a bit. do11j=m1+2-i,mm a(i,j-l)=a(i,j) enddo 11 l=l-1do 12j=mm-l,mm a(i,j)=0. enddo 12 enddo 13 d=1.l=m1 do 18k=1,n For each row... dum=a(k,1)i=kif(l.lt.n)l=l+1 do 14j=k+1,l Find the pivot element. if(abs(a(j,1)).gt.abs(dum))then dum=a(j,1)i=j endif enddo 14 indx(k)=i if(dum.eq.0.) a(k,1)=TINY Matrix is algorithmically singular, but proceed anyway with TINY pivot (desirable in some applications). if(i.ne.k)then Interchange rows. d=-d do15j=1,mm dum=a(k,j) a(k,j)=a(i,j) a(i,j)=dum enddo 15 endif do17i=k+1,l Do the elimination. dum=a(i,1)/a(k,1)al(k,i-k)=dum do 16j=2,mm 46 Chapter2. Solutionof LinearAlgebraicEquationsSample 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).a(i,j-1)=a(i,j)-dum*a(k,j) enddo 16 a(i,mm)=0. enddo 17 enddo 18 returnEND Some pivoting is possible within the storage limitations of bandec, and the above routine does take advantage of the opportunity. In general, when TINYis returned as a diagonal element of U, then the original matrix (perhaps as modified by roundoff error) is in fact singular. In this regard, bandecis somewhat more robust than tridagabove, whichcanfailalgorithmicallyevenfornonsingular matrices; bandecisthusalsouseful(with m1=m2=1) for some ill-behaved tridiagonal systems. Oncethematrix Ahasbeendecomposed,anynumberofright-handsidescanbesolvedin turnbyrepeatedcallsto banbks,thebacksubstitutionroutinewhoseanalogin §2.3is lubksb. SUBROUTINE banbks(a,n,m1,m2,np,mp,al,mpl,indx,b) INTEGER m1,m2,mp,mpl,n,np,indx(n)REAL a(np,mp),al(np,mpl),b(n) Given the arrays a,al,a n d indx as returned from bandec , and given a right-hand side vector b(1:n) , solves the band diagonal linear equations A·x=b. The solution vector x overwrites b(1:n) . The other input arrays are not modified, and can be left in place for successive calls with different right-hand sides. INTEGER i,k,l,mm REAL dummm=m1+m2+1if(mm.gt.mp.or.m1.gt.mpl.or.n.gt.np) pause ’bad args in banbks’ l=m1 do 12k=1,n Forward substitution, unscrambling the permuted rows as we go. i=indx(k) if(i.ne.k)then dum=b(k)b(k)=b(i)b(i)=dum endif if(l.lt.n)l=l+1do 11i=k+1,l b(i)=b(i)-al(k,i-k)*b(k) enddo 11 enddo 12 l=1do 14i=n,1,-1 Backsubstitution. dum=b(i)do 13k=2,l dum=dum-a(i,k)*b(k+i-1) enddo 13 b(i)=dum/a(i,1) if(l.lt.mm) l=l+1 enddo 14 return END The routines bandecandbanbksare based on the Handbook routines bandet1and bansol1in[1]. CITED REFERENCES AND FURTHER READING: Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell), p. 74. 2.5IterativeImprovementofaSolutiontoLinearEquations 47Sample 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).Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Example 5.4.3, p. 166. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §9.11. Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag), Chapter I/6. [1] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §4.3. 2.5 Iterative Improvement of a Solution to Linear Equations Obviously it is not easy to obtain greater precision for the solution of a linear set than the precision of your computer’s floating-point word. Unfortunately, forlarge sets of linear equations, it is not always easy to obtain precision equal to, or even comparable to, the computer’s limit. In direct methods of solution, roundoff errors accumulate, and they are magnified to the extent that your matrix is closeto singular. You can easily lose two or three significant figures for matrices which (you thought) were farfrom singular. Ifthishappenstoyou,thereisa neattricktorestorethefullmachineprecision, callediterative improvement of thesolution. Thetheoryis verystraightforward(see Figure 2.5.1): Suppose that a vector xis the exact solution of the linear set A·x=b (2.5.1 ) You don’t, however, know x. You only know some slightly wrong solution x+δx, where δxistheunknownerror. Whenmultipliedbythematrix A,yourslightlywrong solutiongivesaproductslightlydiscrepantfromthedesiredright-handside b,namely A·(x+δx)=b+δb (2.5.2 ) Subtracting (2.5.1) from (2.5.2) gives A·δx=δb (2.5.3 ) But (2.5.2)can also be solved,trivially,for δb. Substitutingthis into (2.5.3)gives A·δx=A·(x+δx)−b (2.5.4 ) In this equation, the whole right-hand side is known, since x+δxis the wrong solution that you want to improve. It is essential to calculate the right-hand sidein double precision, since there will be a lot of cancellation in the subtraction of b. Then, we need only solve (2.5.4)for the error δx, then subtract this from the wrong solution to get an improved solution. An important extra benefit occurs if we obtained the original solution by LU decomposition. Inthiscase we alreadyhavethe LUdecomposedformof A, andall we need do to solve (2.5.4)is compute the right-handside and backsubstitute! The code to do all this is concise and straightforward: