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<=n;j++){ Decomposition and forward substitution. gam[j]=c[j-1]/bet; bet=b[j]-a[j]*gam[j]; 1f(bet==0.0) nrerror("Error 2 in tridag"); Algorithm fails;see be- u[j]=(r[j]-a[j]*u[j-1])/bet; low. RECIPES I for(j=(n-1);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 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). 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<=n;j++) { Decomposition and forward substitution. gam[j]=c[j-1]/bet; bet=b[j]-a[j]*gam[j]; if (bet == 0.0) nrerror("Error 2 in tridag"); Algorithm fails; see be￾u[j]=(r[j]-a[j]*u[j-1])/bet; low. } for (j=(n-1);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