44 Chapter 2.Solution of Linear Algebraic Equations Equations(2.3.6)and(2.3.7)total (for each right-hand side b)N2 executions of an inner loop containing one multiply and one add.If we have N right-hand sides which are the unit column vectors (which is the case when we are inverting a matrix),then taking into account the leading zeros reduces the total execution count of (2.3.6)from N3 to N3,while (2.3.7)is unchanged at N3 Notice that,once we have the LU decomposition of A,we can solve with as many right-hand sides as we then care to,one at a time.This is a distinct advantage over the methods of $2.1 and $2.2. Performing the LU Decomposition How then can we solve for L and U,given A?First,we write out the i,ith component of equation(2.3.1)or(2.3.2).That component always is a sum beginning with ICAL 01B1j十···=ai对 RECIPES The number of terms in the sum depends,however,on whether i or j is the smaller number.We have.in fact.the three cases. 9 i<j: 011j+02B21+··+aiP对=ai (2.3.8) i=j:ai1B1j+ai2B2j+...+aiiBii=aij (2.3.9) 9 i>j: a1B1j+a2P2j+…+afi=a (2.3.10) 色 9 Equations(2.3.8)-(2.3.10)total N2 equations for the N2+N unknown a's and B's(the diagonal being represented twice).Since the number of unknowns is greater than the number ofequations,we are invited to specify N of the unknowns arbitrarily and then try to solve for the others.In fact,as we shall see,it is always possible to take ai≡1i=1,.,N (2.3.11) A surprising procedure,now,is Crout's algorithm,which quite trivially solves the set of N2+N equations(2.3.8)-(2.3.11)for all the a's and B's by just arranging 兰丝母 10621 the equations in a certain order!That order is as follows: Numerica .Set ai =1,i =1,...,N (equation 2.3.11). 43126 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 Bij,namely i1 North =a1 (2.3.12) k= (Wheni=1 in 2.3.12 the summation term is taken to mean zero.)Second. fori=j+1,j+2,...,N use (2.3.10)to solve for aii,namely 2.3.13) Be sure to do both procedures before going on to the next j.44 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Equations (2.3.6) and (2.3.7) total (for each right-hand side b) N 2 executions of an inner loop containing one multiply and one add. If we have N right-hand sides which are the unit column vectors (which is the case when we are inverting a matrix), then taking into account the leading zeros reduces the total execution count of (2.3.6) from 1 2N3 to 1 6N3, while (2.3.7) is unchanged at 1 2N3. Notice that, once we have the LU decomposition of A, we can solve with as many right-hand sides as we then care to, one at a time. This is a distinct advantage over the methods of §2.1 and §2.2. Performing the LU Decomposition How then can we solve for L and U, given A? First, we write out the i, jth component of equation (2.3.1) or (2.3.2). That component always is a sum beginning with αi1β1j + ··· = aij The number of terms in the sum depends, however, on whether i or j is the smaller number. We have, in fact, the three cases, i<j : αi1β1j + αi2β2j + ··· + αiiβij = aij (2.3.8) i = j : αi1β1j + αi2β2j + ··· + αiiβjj = aij (2.3.9) i>j : αi1β1j + αi2β2j + ··· + αijβjj = aij (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 greater than the number of equations, we are invited to specify N of the unknowns arbitrarily and then try to solve for the others. In fact, as we shall see, it is always possible to take αii ≡ 1 i = 1,...,N (2.3.11) A surprising procedure, now, is Crout’s algorithm, which quite trivially solves the set of N2 + N equations (2.3.8)–(2.3.11) for all the α’s and β’s by just arranging the 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 k=1 αikβkj . (2.3.12) (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 , namely αij = 1 βjj  aij − j−1 k=1 αikβkj . (2.3.13) Be sure to do both procedures before going on to the next j