50 Chapter 2.Solution of Linear Algebraic Equations A quick-and-dirty way to solve complex systems is to take the real and imaginary parts of (2.3.16),giving A·x-C·y=b (2.3.17) C·x+A·y=d which can be written as a 2N x 2N set of real equations, (è)(（)=() (2.3.18) and then solved with ludcmp and lubksb in their present forms.This scheme is a factor of 2 inefficient in storage,since A and C are stored twice.It is also a factor of 2 inefficient in time,since the complex multiplies in a complexified version of the routines would each 81 use 4 real multiplies,while the solution of a 2N x 2N problem involves 8 times the work of an N x N one.If you can tolerate these factor-of-two inefficiencies,then equation (2.3.18) is an easy way to proceed. ICAL CITED REFERENCES AND FURTHER READING: Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),Chapter 4. RECIPES 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),83.3,and p.50. Forsythe.G.E.,and Moler,C.B.1967.Computer Solution of Linear Algebraic Systems(Engle- 代》 Press. wood Cliffs,NJ:Prentice-Hall),Chapters 9,16,and 18. Westlake,J.R.1968,A Handbook of Numerical Matrix Inversion and Solution of Linear Equations (New York:Wiley). Stoer.J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). 84.2. SCIENTIFIC Ralston,A.,and Rabinowitz,P.1978,A First Course in Numerical Analysis,2nd ed.(New York: McGraw-Hill),89.11. 6 Horn,R.A.,and Johnson,C.R.1985,Matrix Analysis(Cambridge:Cambridge University Press). 2.4 Tridiagonal and Band Diagonal Systems of Equations Numerica 10621 uctio 43106 The special case of a system of linear equations that is tridiagonal,that is,has (outside Recipes nonzero elements only on the diagonal plus or minus one column,is one that occurs frequently.Also common are systems that are band diagonal,with nonzero elements only along a few diagonal lines adjacent to the main diagonal (above and below). North Software. For tridiagonal sets,the procedures of LU decomposition,forward-and back- substitution each take only O(N)operations,and the whole solution can be encoded very concisely.The resulting routine tridag is one that we will use in later chapters Naturally,one does not reserve storage for the full N x N matrix,but only for the nonzero components,stored as three vectors.The set ofequations to be solved is [b1c10- u1 a2b2c2··· u2 T2 (2.4.1) ·aN-1bN-1 CN-1 UN-1 TN-1 . uN TN
50 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 machinereadable 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). A quick-and-dirty way to solve complex systems is to take the real and imaginary parts of (2.3.16), giving A · x − C · y = b C · x + A · y = d (2.3.17) which can be written as a 2N × 2N set of real equations, A −C C A · x y = b d (2.3.18) and then solved with ludcmp and lubksb in their present forms. This scheme is a factor of 2 inefficient in storage, since A and C are stored twice. It is also a factor of 2 inefficient in time, since the complex multiplies in a complexified version of the routines would each use 4 real multiplies, while the solution of a 2N × 2N problem involves 8 times the work of an N × N one. 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., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins 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 (Englewood Cliffs, NJ: Prentice-Hall), Chapters 9, 16, and 18. Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations (New York: Wiley). Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: 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., and Johnson, C.R. 1985, Matrix Analysis (Cambridge: Cambridge University Press). 2.4 Tridiagonal and Band Diagonal Systems of Equations The special case of a system of linear equations that is tridiagonal, that is, has nonzero elements only on the diagonal plus or minus one column, is one that occurs frequently. Also common are systems that are band diagonal, with nonzero elements only along a few diagonal lines adjacent to the main diagonal (above and below). For tridiagonal sets, the procedures of LU decomposition, forward- and backsubstitution each take only O(N) operations, and the whole solution can be encoded very concisely. The resulting routine tridag is one that we will use in later chapters. Naturally, one does not reserve storage for the full N × N matrix, but only for the nonzero components, stored as three vectors. The set of equations to be solved is b1 c1 0 ··· a2 b2 c2 ··· ··· ··· aN−1 bN−1 cN−1 ··· 0 aN bN · u1 u2 ··· uN−1 uN = r1 r2 ··· rN−1 rN (2.4.1)
2.4 Tridiagonal and Band Diagonal Systems of Equations 51 Notice that a and cy are undefined and are not referenced by the routine that follows. #include "nrutil.h" void tridag(float a[],float b[],float c[],float r[],float u[], unsigned long n) Solves for a vector u[1..n]the tridiagonal linear set given by equation (2.4.1).a[1..n], b[1..n],c[1..n],and r[1..n]are input vectors and are not modified. unsigned long j; float bet,*gami gam=vector(1,n); One vector of workspace,gam is needed. if (b[1]==0.0)nrerror("Error 1 in tridag"); If this happens then you should rewrite your equations as a set of order N-1,with u2 trivially eliminated. u[1]=r[1]/(bet=b[1]); for(j=2;j=1j-) u[j]-=gam[j+1]*u[j+1]; Backsubstitution. free_vector(gam,1,n); 子g Press. There is no pivoting in tridag.It is for this reason that tridag can fail even 9 when the underlying matrix is nonsingular:A zero pivot can be encountered even for a nonsingular matrix.In practice,this is not something to lose sleep about.The kinds SCIENTIFIC of problems that lead to tridiagonal linear sets usually have additional properties which guarantee that the algorithm in tridag will succeed.For example,if 61 loil>lajl +lcil j=1,.,N (2.4.2) (called diagonal dominance)then it can be shown that the algorithm cannot encounter a zero pivot. It is possible to construct special examples in which the lack of pivoting in the Numerica 10621 algorithm causes numerical instability.In practice,however,such instability is almost 43106 never encountered-unlike the general matrix problem where pivoting is essential. The tridiagonal algorithm is the rare case of an algorithm that,in practice,is Recipes more robust than theory says it should be.Of course,should you ever encounter a problem for which tridag fails,you can instead use the more general method for North Software. band diagonal systems,now described(routines bandec and banbks). 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 Diagonal Systems Where tridiagonal systems have nonzero elements only on the diagonal plus or minus one,band diagonal systems are slightly more general and have(say)m>0 nonzero elements immediately to the left of(below)the diagonal and m2>0 nonzero elements immediately to its right (above it).Of course,this is only a useful classification if m and m2 are bothN
2.4 Tridiagonal and Band Diagonal Systems of Equations 51 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 machinereadable 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). Notice that a1 and cN are undefined and are not referenced by the routine that follows. #include "nrutil.h" void tridag(float a[], float b[], float c[], float r[], float u[], unsigned long n) Solves for a vector u[1..n] the tridiagonal linear set given by equation (2.4.1). a[1..n], b[1..n], c[1..n], and r[1..n] are input vectors and are not modified. { unsigned long j; float bet,*gam; gam=vector(1,n); One vector of workspace, gam is needed. if (b[1] == 0.0) nrerror("Error 1 in tridag"); If this happens then you should rewrite your equations as a set of order N − 1, w ith u2 trivially eliminated. u[1]=r[1]/(bet=b[1]); for (j=2;j=1;j--) u[j] -= gam[j+1]*u[j+1]; Backsubstitution. free_vector(gam,1,n); } There is no pivoting in tridag. It is for this reason that tridag can fail even when the underlying matrix is nonsingular: A zero pivot can be encountered even for a nonsingular matrix. In practice, this is not something to lose sleep about. The kinds of problems that lead to tridiagonal linear sets usually have additional properties which guarantee that the algorithm in tridag will succeed. For example, if |bj| > |aj | + |cj | j = 1,...,N (2.4.2) (called diagonal dominance) then it can be shown that the algorithm cannot encounter a zero pivot. It is possible to construct special examples in which the lack of pivoting in the algorithm causes numerical instability. In practice, however, such instability is almost never encountered — unlike the general matrix problem where pivoting is 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 a problem for which tridag fails, you can instead use the more general method for band diagonal systems, now described (routines bandec and banbks). 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 Diagonal Systems Where tridiagonal systems have nonzero elements only on the diagonal plus or minus one, band diagonal systems are slightly more general and have (say) m1 ≥ 0 nonzero elements immediately to the left of (below) the diagonal and m2 ≥ 0 nonzero elements immediately to its right (above it). Of course, this is only a useful classification if m1 and m2 are both N
52 Chapter 2.Solution of Linear Algebraic Equations In that case,the solution of the linear system by LU decomposition can be accomplished much faster,and in much less storage,than for the general N x N case. The precise definition of a band diagonal matrix with elements is that aij=0 when j>i+m2 or i>j+m (2.4.3) Band diagonal matrices are stored and manipulated in a so-called compact form,which results if the matrix is tilted 45 clockwise,so that its nonzero elements lie in a long,narrow matrix with m+1 +m2 columns and N rows.This is best illustrated by an example: The band diagonal matrix 3 100000 ￥ 15000 0 9 6 令 0 0 0 0 9 0 0 (2.4.4) granted for 18881892 0 0 9 2 0 g 0 0 0 2 38 4 1-800 0 0 0 0 2 4 from NUMERICAL RECIPES I which has N =7,m=2,and m2 =1,is stored compactly as the 7 x 4 matrix, 3 S 4 (North 1 5 9 2 6 5 3 8 9 (2.4.5) computer Ameri Press. 9 2 ART 3 8 4 6 2 4 4 9 Program Here x denotes 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 diagonal of the original matrix appears in column m+1,with subdiagonal elements to its left, superdiagonal elements to its right. to dir 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 study the following routine carefully,as an example of how to pull nonzero elements aj out of the OF SCIENTIFIC COMPUTING(ISBN compact storage format in an orderly fashion. 198918920 #include "nrutil.h" void banmul(float **a,unsigned long n,int m1,int m2,float x[],float b[]) Numerical Recipes 10-621 43108 Matrix multiply b =A.x,where A is band diagonal with m1 rows below the diagonal and m2 rows above.The input vector x and output vector b are stored as x[1..n]and b[1..n], respectively. The array a[1..n][1..m1+m2+1]stores A as follows:The diagonal elements (outside are in a[1..n][m1+1].Subdiagonal elements are in a[j..n][1..m1](with>1 ap- 膜 propriate to the number of elements on each subdiagonal).Superdiagonal elements are in a[1..j][m1+2..m1+m2+1]with j<n appropriate to the number of elements on each su- North Software. perdiagonal. Ame unsigned long i,j,k,tmploop; for(i=1;1<=n;1+){ k=1-m1-1; tmploop=LMIN(m1+m2+1,n-k); b[i]=0.0; for (j=LMAX(1,1-k);j<=tmploop;j++)b[i]+a[i][j]*x[j+k];
52 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 machinereadable 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). In that case, the solution of the linear system by LU decomposition can be accomplished much faster, and in much less storage, than for the general N × N case. The precise definition of a band diagonal matrix with elements aij is that aij = 0 when j>i + m2 or i>j + m1 (2.4.3) Band diagonal matrices are stored and manipulated in a so-called compact form, which results if the matrix is tilted 45◦ clockwise, so that its nonzero elements lie in a long, narrow matrix with m1 +1+ m2 columns and N rows. This is best illustrated by an example: The band diagonal matrix 3100000 4150000 9265000 0358900 0079320 0003846 0000244 (2.4.4) which has N = 7, m1 = 2, and m2 = 1, is stored compactly as the 7 × 4 matrix, x x 3 1 x 415 9265 3589 7932 3846 244 x (2.4.5) Here x denotes 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 diagonal of the original matrix appears in column m1 + 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 study the following routine carefully, as an example of how to pull nonzero elements aij out of the compact storage format in an orderly fashion. #include "nrutil.h" void banmul(float **a, unsigned long n, int m1, int m2, float x[], float b[]) Matrix multiply b = A · x, w here A is band diagonal with m1 rows below the diagonal and m2 rows above. The input vector x and output vector b are stored as x[1..n] and b[1..n], respectively. The array a[1..n][1..m1+m2+1] stores A as 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 are in a[1..j][m1+2..m1+m2+1] with j < n appropriate to the number of elements on each superdiagonal. { unsigned long i,j,k,tmploop; for (i=1;i<=n;i++) { k=i-m1-1; tmploop=LMIN(m1+m2+1,n-k); b[i]=0.0; for (j=LMAX(1,1-k);j<=tmploop;j++) b[i] += a[i][j]*x[j+k]; } }
2.4 Tridiagonal and Band Diagonal Systems of Equations 53 It is not possible to store the LU decomposition of a band diagonal matrix A quite as compactly as the compact form of A itself.The decomposition (essentially by Crout's method,see 82.3)produces additional nonzero"fill-ins."One straightforward storage scheme is to return the upper triangular factor (in the same space that A previously occupied,and to return the lower triangular factor(L)in a separate compact matrix of size N x m.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 ludcmp in 82.3: #include #define SWAP(a,b){dum=(a);(a)=(b);(b)=dum;} #define TINY 1.0e-20 83 void bandec(float **a,unsigned long n,int m1,int m2,float **al unsigned long indx,float *d) 鱼 granted for 19881992 Given an n x n band diagonal matrix A with m1 subdiagonal rows and m2 superdiagonal rows, 1800 豆 compactly stored in the array a[1..n][1..m1+m2+1]as described in the comment for routine banmul,this routine constructs an LU decomposition of a rowwise permutation of A.The upper triangular matrix replaces a,while the lower triangular matrix is returned in al[1..n][1..m1] indx[1..n]is an output vector which records the row permutation effected by the partial from NUMERICAL RECIPES IN pivoting:d is output as +1 depending 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. unsigned long i,j,k,1; server computer, (North America to make one paper THE int mm; float dum; ART 是 mm=m1+m2+1; st st Programs 1=m1; for(i=1;1fabs(dum)){ dum=a[j][1]; @cambridge.org 1988-1992 by Numerical Recipes 10-:6211 43106 1=j; indx[k]=i; (outside if (dum =0.0)a[k][1]=TINY; Software. Matrix is algorithmically singular,but proceed anyway with TINY pivot (desirable in some applications). 1f(1I=k)[ Interchange rows Ame *d=-(*d); for (j=1;j<=mm;j++)SWAP(a[k][j],a[i][i]) for(1=k+1;1<=1;i+) Do the elimination dum=a[i][1]/a[k][1]; al[k][i-k]=dum; for (j=2;j<=mm;j++)a[i][j-1]=a[i][j]-dum*a[k][j]; a[i][mm]=0.0;
2.4 Tridiagonal and Band Diagonal Systems of Equations 53 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 machinereadable 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). It is not possible to store the LU decomposition of a band diagonal matrix A quite as compactly as the compact form of A itself. The decomposition (essentially by Crout’s method, see §2.3) produces additional nonzero “fill-ins.” One straightforward storage scheme is to return the upper triangular factor (U) in the same space that A previously 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 ludcmp in §2.3: #include #define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;} #define TINY 1.0e-20 void bandec(float **a, unsigned long n, int m1, int m2, float **al, unsigned long indx[], float *d) Given an n × n band diagonal matrix A with m1 subdiagonal rows and m2 superdiagonal 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 LU decomposition of a rowwise permutation of A. The upper triangular matrix replaces a, while the lower triangular matrix is returned in al[1..n][1..m1]. indx[1..n] is an output vector which records the row permutation effected by the partial pivoting; d is output as ±1 depending 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. { unsigned long i,j,k,l; int mm; float dum; mm=m1+m2+1; l=m1; for (i=1;i fabs(dum)) { dum=a[j][1]; i=j; } } indx[k]=i; if (dum == 0.0) a[k][1]=TINY; Matrix is algorithmically singular, but proceed anyway with TINY pivot (desirable in some applications). if (i != k) { Interchange rows. *d = -(*d); for (j=1;j<=mm;j++) SWAP(a[k][j],a[i][j]) } for (i=k+1;i<=l;i++) { Do the elimination. dum=a[i][1]/a[k][1]; al[k][i-k]=dum; for (j=2;j<=mm;j++) a[i][j-1]=a[i][j]-dum*a[k][j]; a[i][mm]=0.0; } } }
54 Chapter 2.Solution of Linear Algebraic Equations Some pivoting is possible within the storage limitations of bandec,and the above routine does take advantage of the opportunity.In general,when TINY is returned as a diagonal element of 0,then the original matrix (perhaps as modified by roundoff error) is in fact singular.In this regard,bandec is somewhat more robust than tridag above, which can fail algorithmically even for nonsingular matrices;bandec is thus also useful (with m=m2=1)for some ill-behaved tridiagonal systems. Once the matrix A has been decomposed,any number ofright-hand sides can be solved in turn by repeated calls to banbks,the backsubstitution routine whose analog in 82.3 is lubksb. #define SWAP(a,b){dum=(a)(a)=(b);(b)=dum; void banbks(float *a,unsigned long n,int mi,int m2,float **al, 81 nsigned1 Long indx☐，f1oatb[]) Given the arrays a,al,and 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. unsigned long i,k,1; 3 int mm; float dum; RECIPES mm=m1+m2+1; 1m1; for (k=1;k=1;i--)[ Backsubstitution. SCIENTIFIC dum=b[i]; for(k=2;k<=1;k++)dum-=a[i][k]*b[k+i-1]; b[i]=dum/a[i][1]; to dir 1f(1<mm)1++: 2 192 COMPUTING (ISBN The routines bandec and banbks are based on the Handbook routines bandet/and bansoll in [1]. Fuunrggoirioh Numerica 10.621 43106 CITED REFERENCES AND FURTHER READING: Recipes Keller,H.B.1968,Numerical Methods for Two-Point Boundary-Value Problems (Waltham,MA: Blaisdell),p.74. Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall), North Software. Example 5.4.3,p.166. Ralston,A.,and Rabinowitz,P.1978,A First Course in Numerical Analysis,2nd ed.(New York: McGraw-Hill),89.11. Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.ll of Handbook for Automatic Com- putation (New York:Springer-Verlag),Chapter 1/6.[1] Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),84.3
54 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 machinereadable 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). Some pivoting is possible within the storage limitations of bandec, and the above routine does take advantage of the opportunity. In general, when TINY is returned as a diagonal element of U, then the original matrix (perhaps as modified by roundoff error) is in fact singular. In this regard, bandec is somewhat more robust than tridag above, which can fail algorithmically even for nonsingular matrices; bandec is thus also useful (with m1 = m2 = 1) for some ill-behaved tridiagonal systems. Once the matrix A has been decomposed, any number of right-hand sides can be solved in turn by repeated calls to banbks, the backsubstitution routine whose analog in §2.3 is lubksb. #define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;} void banbks(float **a, unsigned long n, int m1, int m2, float **al, unsigned long indx[], float b[]) Given the arrays a, al, and 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. { unsigned long i,k,l; int mm; float dum; mm=m1+m2+1; l=m1; for (k=1;k=1;i--) { Backsubstitution. dum=b[i]; for (k=2;k<=l;k++) dum -= a[i][k]*b[k+i-1]; b[i]=dum/a[i][1]; if (l < mm) l++; } } The routines bandec and banbks are based on the Handbook routines bandet1 and bansol1 in [1]. CITED REFERENCES AND FURTHER READING: Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell), p. 74. 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 Computation (New York: Springer-Verlag), Chapter I/6. [1] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §4.3
2.5 Iterative Improvement of a Solution to Linear Equations 55 x+δx b+ob 83 旦 granted for 11800 (including this one) /Cambridge from NUMERICAL RECIPES IN 19881992 Figure 2.5.1.Iterative improvement of the solution to A.x b.The first guess x +6x is multiplied by A to produce b+5b.The known vector b is subtracted,giving b.The linear set with this right-hand (Nort server side is inverted,giving &x.This is subtracted from the first guess giving an improved solution x. America computer, users to make one paper UnN电.t THE ART 2.5 Iterative Improvement of a Solution to Linear Equations 9 ictly proh Progra 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,for to dir large 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 OF SCIENTIFIC COMPUTING(ISBN errors accumulate,and they are magnified to the extent that your matrix is close 1988-1992 to singular.You can easily lose two or three significant figures for matrices which (you thought)were far from singular. 10-621 If this happens to you,there is a neat trick to restore the full machine precision, 43106 called iterative improvement of the solution.The theory is very straightforward(see Numerical Recipes Figure 2.5.1):Suppose that a vector x is the exact solution of the linear set (outside A·x=b (2.5.1) Software. You don't,however,know x.You only know some slightly wrong solution x +ox, where &x is the unknownerror.When multiplied by the matrix A,your slightly wrong visit website solution gives a product slightly discrepant from the desired right-hand side b,namely machine A·(x+6x)=b+b (2.5.2) Subtracting (2.5.1)from (2.5.2)gives A·6x=6b (2.5.3)
2.5 Iterative Improvement of a Solution to Linear Equations 55 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 machinereadable 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). A A−1 δx x + δx x b b + δb δb Figure 2.5.1. Iterative improvement of the solution to A · x = b. The first guess x + δx is multiplied by A to produce b + δb. The known vector b is subtracted, giving δb. The linear set with this right-hand side is inverted, giving δx. This is subtracted from the first guess giving an improved solution x. 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, for large 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 close to singular. You can easily lose two or three significant figures for matrices which (you thought) were far from singular. If this happens to you, there is a neat trick to restore the full machine precision, called iterative improvement of the solution. The theory is very straightforward (see Figure 2.5.1): Suppose that a vector x is 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 δx is the unknown error. When multiplied by the matrix A, your slightly wrong solution gives a product slightly discrepant from the desired right-hand side 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)