2.7 Sparse Linear Systems 89 void dsprstx(double sa[],unsigned long ija[],double x[],double b[], unsigned long n); These are double versions of sprsax and sprstx. if (itrnsp)dsprstx(sa,ija,x,r,n); else dsprsax(sa,ija,x,r,n); 2 extern unsigned long ija☐； extern double sa☐； The matrix is stored somewhere void asolve(unsigned long n,double b[],double x[],int itrnsp) 83g 19881992 unsigned long i; for(i=1;i<=n;i+)x[i]=(sa[i]I=0.0？b[i]/sa[i]:b[i]); 三之h会州 1-800 The matrix A is the diagonal part of A,stored in the first n elements of sa.Since the transpose matrix has the same diagonal,the flag itrnsp is not used. from NUMERICAL RECIPES I America computer, Press. CITED REFERENCES AND FURTHER READING: Tewarson,R.P.1973,Sparse Matrices (New York:Academic Press).[1] Jacobs,D.A.H.(ed.)1977,The State of the Art in Numerical Analysis (London:Academic Press). 9 Chapter 1.3 (by J.K.Reid).[2] George,A.,and Liu,J.W.H.1981,Computer Solution of Large Sparse Positive Definite Systems SCIENTIFIC (Englewood Cliffs,NJ:Prentice-Hall).[3] NAG Fortran Library (Numerical Algorithms Group,256 Banbury Road,Oxford OX27DE,U.K.). [4 IMSL Math/Library Users Manual (IMSL Inc..2500 CityWest Boulevard,Houston TX 77042).[5] Eisenstat,S.C.,Gursky,M.C.,Schultz,M.H.,and Sherman,A.H.1977,Yale Sparse Matrix Pack- age,Technical Reports 112 and 114(Yale University Department of Computer Science).[6] Knuth.D.E.1968.Fundamental Algorithms,vol.1 of The Art of Computer Programming(Reading. MA:Addison-Wesley),82.2.6.[7] Kincaid,D.R.,Respess,J.R.,Young,D.M.,and Grimes,R.G.1982,ACM Transactions on Math- ematical Software,vol.8,pp.302-322.[8] Numerical Recipes 10621 43106 PCGPAK User's Guide (New Haven:Scientific Computing Associates,Inc.).[9] Bentley.J.1986,Programming Pearls(Reading,MA:Addison-Wesley).$9.[10] (outside Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press).Chapters 4 and 10,particularly $10.2-10.3.[11] North Software. Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). Chapter 8.[12] Baker,L.1991,More C Tools for Scientists and Engineers (New York:McGraw-Hill).[13] Fletcher,R.1976,in Numerical Analysis Dundee 1975,Lecture Notes in Mathematics,vol.506 A.Dold and B Eckmann,eds.(Berlin:Springer-Verlag),pp.73-89.[14] Saad,Y.,and Schulz,M.1986,SIAM Journal on Scientific and Statistical Computing.vol.7, pp.856-869.[15] Bunch.J.R.,and Rose.D.J.(eds.)1976.Sparse Matrix Computations (New York:Academic Press). Duff,I.S.,and Stewart,G.W.(eds.)1979,Sparse Matrix Proceedings 1978 (Philadelphia S.IA.M.)
2.7 Sparse Linear Systems 89 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). void dsprstx(double sa[], unsigned long ija[], double x[], double b[], unsigned long n); These are double versions of sprsax and sprstx. if (itrnsp) dsprstx(sa,ija,x,r,n); else dsprsax(sa,ija,x,r,n); } extern unsigned long ija[]; extern double sa[]; The matrix is stored somewhere. void asolve(unsigned long n, double b[], double x[], int itrnsp) { unsigned long i; for(i=1;i<=n;i++) x[i]=(sa[i] != 0.0 ? b[i]/sa[i] : b[i]); The matrix A is the diagonal part of A, stored in the first n elements of sa. Since the transpose matrix has the same diagonal, the flag itrnsp is not used. } CITED REFERENCES AND FURTHER READING: Tewarson, R.P. 1973, Sparse Matrices (New York: Academic Press). [1] Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: Academic Press), Chapter I.3 (by J.K. Reid). [2] George, A., and Liu, J.W.H. 1981, Computer Solution of Large Sparse Positive Definite Systems (Englewood Cliffs, NJ: Prentice-Hall). [3] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.). [4] IMSL Math/Library Users Manual (IMSL Inc., 2500 CityWest Boulevard, Houston TX 77042). [5] Eisenstat, S.C., Gursky, M.C., Schultz, M.H., and Sherman, A.H. 1977, Yale Sparse Matrix Package, Technical Reports 112 and 114 (Yale University Department of Computer Science). [6] Knuth, D.E. 1968, Fundamental Algorithms, vol. 1 of The Art of Computer Programming (Reading, MA: Addison-Wesley), §2.2.6. [7] Kincaid, D.R., Respess, J.R., Young, D.M., and Grimes, R.G. 1982, ACM Transactions on Mathematical Software, vol. 8, pp. 302–322. [8] PCGPAK User’s Guide (New Haven: Scientific Computing Associates, Inc.). [9] Bentley, J. 1986, Programming Pearls (Reading, MA: Addison-Wesley), §9. [10] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), Chapters 4 and 10, particularly §§10.2–10.3. [11] Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), Chapter 8. [12] Baker, L. 1991, More C Tools for Scientists and Engineers (New York: McGraw-Hill). [13] Fletcher, R. 1976, in Numerical Analysis Dundee 1975, Lecture Notes in Mathematics, vol. 506, A. Dold and B Eckmann, eds. (Berlin: Springer-Verlag), pp. 73–89. [14] Saad, Y., and Schulz, M. 1986, SIAM Journal on Scientific and Statistical Computing, vol. 7, pp. 856–869. [15] Bunch, J.R., and Rose, D.J. (eds.) 1976, Sparse Matrix Computations (New York: Academic Press). Duff, I.S., and Stewart, G.W. (eds.) 1979, Sparse Matrix Proceedings 1978 (Philadelphia: S.I.A.M.)
90 Chapter 2.Solution of Linear Algebraic Equations 2.8 Vandermonde Matrices and Toeplitz Matrices In 82.4 the case of a tridiagonal matrix was treated specially,because that particular type of linear system admits a solution in only of order N operations, rather than of order N3 for the general linear problem.When such particular types exist,it is important to know about them.Your computational savings,should you ever happen to be working on a problem that involves the right kind of particular type,can be enormous. This section treats two special types of matrices that can be solved in of order N2 operations,not as good as tridiagonal,but a lot better than the general case. ted for (Other than the operations count,these two types having nothing in common.) Matrices of the first type,termed Vandermonde matrices,occur in some problems having to do with the fitting of polynomials,the reconstruction of distributions from their moments.and also other contexts.In this book.for example.a Vandermonde RECIPES I problem crops up in 83.5.Matrices of the second type,termed Toeplit=matrices, 2 tend to occur in problems involving deconvolution and signal processing.In this book,a Toeplitz problem is encountered in $13.7. These are not the only special types of matrices worth knowing about.The Press. Hilbert matrices,whose components are of the form aij=1/(i+j-1),i,j= 1,...,N can be inverted by an exact integer algorithm,and are very difficult to invert in any other way,since they are notoriously ill-conditioned (see [1]for details) 9 The Sherman-Morrison and Woodbury formulas,discussed in $2.7,can sometimes be used to convert new special forms into old ones.Reference [2]gives some other OF SCIENTIFIC( special forms.We have not found these additional forms to arise as frequently as 6 the two that we now discuss. 1988-19920 COMPUTING (ISBN Vandermonde Matrices 10-521 A Vandermonde matrix of size N x N is completely determined by N arbitrary Numerical Recipes 43106 numbers 1,22,...,N,in terms of which its N2 components are the integer powers 1...N.Evidently there are two possible such forms,depending on whether (outside we view the i's as rows,i's as columns,or vice versa.In the former case,we get a linear system of equations that looks like this, North Software. 2 2 4+ Y2 (2.8.1) : x一 CN N Performing the matrix multiplication,you will see that this equation solves for the unknown coefficients c;which fit a polynomial to the N pairs of abscissas and ordinates (xj,y). Precisely this problem will arise in 83.5,and the routine given there will solve (2.8.1)by the method that we are about to describe
90 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). 2.8 Vandermonde Matrices and Toeplitz Matrices In §2.4 the case of a tridiagonal matrix was treated specially, because that particular type of linear system admits a solution in only of order N operations, rather than of order N 3 for the general linear problem. When such particular types exist, it is important to know about them. Your computational savings, should you ever happen to be working on a problem that involves the right kind of particular type, can be enormous. This section treats two special types of matrices that can be solved in of order N2 operations, not as good as tridiagonal, but a lot better than the general case. (Other than the operations count, these two types having nothing in common.) Matrices of the first type, termed Vandermonde matrices, occur in some problems having to do with the fitting of polynomials, the reconstruction of distributions from their moments, and also other contexts. In this book, for example, a Vandermonde problem crops up in §3.5. Matrices of the second type, termed Toeplitz matrices, tend to occur in problems involving deconvolution and signal processing. In this book, a Toeplitz problem is encountered in §13.7. These are not the only special types of matrices worth knowing about. The Hilbert matrices, whose components are of the form aij = 1/(i + j − 1), i, j = 1,...,N can be inverted by an exact integer algorithm, and are very difficult to invert in any other way, since they are notoriously ill-conditioned (see [1] for details). The Sherman-Morrison and Woodbury formulas, discussed in §2.7, can sometimes be used to convert new special forms into old ones. Reference [2] gives some other special forms. We have not found these additional forms to arise as frequently as the two that we now discuss. Vandermonde Matrices A Vandermonde matrix of size N × N is completely determined by N arbitrary numbers x1, x2,...,xN , in terms of which its N2 components are the integer powers xj−1 i , i, j = 1,...,N. Evidently there are two possible such forms, depending on whether we view the i’s as rows, j’s as columns, or vice versa. In the former case, we get a linear system of equations that looks like this, 1 x1 x2 1 ··· xN−1 1 1 x2 x2 2 ··· xN−1 2 . . . . . . . . . . . . 1 xN x2 N ··· xN−1 N · c1 c2 . . . cN = y1 y2 . . . yN (2.8.1) Performing the matrix multiplication, you will see that this equation solves for the unknown coefficients ci which fit a polynomial to the N pairs of abscissas and ordinates (xj , yj ). Precisely this problem will arise in §3.5, and the routine given there will solve (2.8.1) by the method that we are about to describe
2.8 Vandermonde Matrices and Toeplitz Matrices 91 The alternative identification of rows and columns leads to the set of equations 1 1 1 W1 工2 2 名 N W3 3 (2.8.2) 21 4。 IN-I N Write this out and you will see that it relates to the problem of moments:Given the values of N points ri,find the unknown weights wi,assigned so as to match the given values gi of the first N moments.(For more on this problem,consult [3).)The routine given in this section solves (2.8.2). EOEE 8 The method of solution of both (2.8.1)and (2.8.2)is closely related to Lagrange's polynomial interpolation formula,which we will not formally meet until $3.I below.Notwith- standing,the following derivation should be comprehensible: Let Pi()be the polynomial of degree N-1 defined by N ICAL P(x)= E一tn Ij -In =∑Ax-1 (2.8.3) k三1 RECIPES Here the meaning of the last equality is to define the components of the matrix As as the coefficients that arise when the product is multiplied out and like terms collected. The polynomial Pi()is a function of generally.But you will notice that it is specifically designed so that it takes on a value of zero at all i with ij,and has a value of unity at x =j.In other words, Press. Programs B(e)==∑Ax-1 (2.8.4) k=1 But(2.8.4)says that A is exactly the inverse of the matrix of components,which appears in (2.8.2),with the subscript as the column index.Therefore the solution of(2.8.2) IENTIFIC is just that matrix inverse times the right-hand side, N =∑A9 (2.8.5) k As for the transpose problem(2.8.1),we can use the fact that the inverse of the transpose is the transpose of the inverse,so Numer (2.8.6) The routine in 83.5 implements this. It remains to find a good way of multiplying out the monomial terms in (2.8.3),in order to get the components of Ajk.This is essentially a bookkeeping problem,and we will let you read the routine itself to see how it can be solved.One trick is to define a master P()by North Software. N P(x)≡ (-2n) (2.8.7) work out its coefficients,and then obtain the numerators and denominators of the specific P,'s viasynthetic division by the one supernumerary term.(See 5.3 for more on synthetic division.) Since each such division is only a process of order N,the total procedure is of order N2 You should be warned that Vandermonde systems are notoriously ill-conditioned,by their very nature.(As an aside anticipating $5.8,the reason is the same as that which makes Chebyshev fitting so impressively accurate:there exist high-order polynomials that are very good uniform fits to zero.Hence roundoff error can introduce rather substantial coefficients of the leading terms of these polynomials.It is a good idea always to compute Vandermonde problems in double precision
2.8 Vandermonde Matrices and Toeplitz Matrices 91 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). The alternative identification of rows and columns leads to the set of equations 1 1 ··· 1 x1 x2 ··· xN x2 1 x2 2 ··· x2 N ··· xN−1 1 xN−1 2 ··· xN−1 N · w1 w2 w3 ··· wN = q1 q2 q3 ··· qN (2.8.2) Write this out and you will see that it relates to the problem of moments: Given the values of N points xi, find the unknown weights wi, assigned so as to match the given values qj of the first N moments. (For more on this problem, consult [3].) The routine given in this section solves (2.8.2). The method of solution of both (2.8.1) and (2.8.2) is closely related to Lagrange’s polynomial interpolation formula, which we will not formally meet until §3.1 below. Notwithstanding, the following derivation should be comprehensible: Let Pj (x) be the polynomial of degree N − 1 defined by Pj (x) = N n=1 (n=j) x − xn xj − xn = N k=1 Ajkxk−1 (2.8.3) Here the meaning of the last equality is to define the components of the matrix Aij as the coefficients that arise when the product is multiplied out and like terms collected. The polynomial Pj (x) is a function of x generally. But you will notice that it is specifically designed so that it takes on a value of zero at all xi with i = j, and has a value of unity at x = xj . In other words, Pj (xi) = δij = N k=1 Ajkxk−1 i (2.8.4) But (2.8.4) says that Ajk is exactly the inverse of the matrix of components xk−1 i , which appears in (2.8.2), with the subscript as the column index. Therefore the solution of (2.8.2) is just that matrix inverse times the right-hand side, wj = N k=1 Ajkqk (2.8.5) As for the transpose problem (2.8.1), we can use the fact that the inverse of the transpose is the transpose of the inverse, so cj = N k=1 Akj yk (2.8.6) The routine in §3.5 implements this. It remains to find a good way of multiplying out the monomial terms in (2.8.3), in order to get the components of Ajk. This is essentially a bookkeeping problem, and we will let you read the routine itself to see how it can be solved. One trick is to define a master P(x) by P(x) ≡ N n=1 (x − xn) (2.8.7) work out its coefficients, and then obtain the numerators and denominators of the specific Pj ’s via synthetic division by the one supernumerary term. (See §5.3 for more on synthetic division.) Since each such division is only a process of order N, the total procedure is of order N2. You should be warned that Vandermonde systems are notoriously ill-conditioned, by their very nature. (As an aside anticipating §5.8, the reason is the same as that which makes Chebyshev fitting so impressively accurate: there exist high-order polynomials that are very good uniform fits to zero. Hence roundoff error can introduce rather substantial coefficients of the leading terms of these polynomials.) It is a good idea always to compute Vandermonde problems in double precision
92 Chapter 2.Solution of Linear Algebraic Equations The routine for (2.8.2)which follows is due to G.B.Rybicki. #include "nrutil.h" void vander(double x[],double w[],double g[],int n) Solves the Vandermonde linear system(1..N).Input consists of the vectors x[1..n]and q[1..n];the vector w[1..n]is output. f int i,j,k; double b,s,t,xx; double米c; c=dvector(1,n); if(n=1)[1]=g[1] 83 else 鱼 granted for 19881992 for(i=1;i<=n;i++)c[i]=0.0: Initialize array. c[n]=-x[1]; Coefficients of the master polynomial are found 1800 for(1=2;i<n;i++)( by recursion. xx=-x[1]; for(j=(n+1-1);j<=(n-1);j++)c[j]+=xx*c[j+1]: to any Cambridge c[n]+xx; from NUMERICAL RECIPES IN 2 for(i=1;i<=n;i++){ Each subfactor in turn (North to make Xx=x[1]; t=b=1.0; America 州bMe se one paper UnN电.t THE s=q[n]; for(k=n;k>=2;k--)[ is synthetically divided, ART b=c[k]+xx*b; s+=qk-1]*b; matrix-multiplied by the right-hand side, t=xx*t+b; st st Programs [1]=s/t; and supplied with a denominator 2 free_dvector(c,1,n); to dir 18881920 OF SCIENTIFIC COMPUTING(ISBN Toeplitz Matrices v@cam 12-521 An N x N Toeplitz matrix is specified by giving 2N-1 numbers R&,k =-N+ 1,...,-1,0,1,...,N-1.Those numbers are then emplaced as matrix elements constant Numerical Recipes 43108 along the (upper-left to lower-right)diagonals of the matrix: Ro R-1 R-2 R-(N-2) R-N-1) Ro R-1 R-(N-3) (outside R-(N-2) R2 Ri Ro R-(N-4) R-(N-3) (2.8.8) North Software. 。。 RN-2 N-3 RN-4 Ro R-1 Ame .RN-1 RN-2 RN-3 Ri Ro The linear Toeplitz problem can thus be written as Ri-jxj =Vi (i=1,,N) (2.8.9) j=】 where the j's,j=1,...,N,are the unknowns to be solved for. The Toeplitz matrix is symmetric if Rk=R-k for all k.Levinson [4]developed an algorithm for fast solution of the symmetric Toeplitz problem,by a bordering method,that is
92 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). The routine for (2.8.2) which follows is due to G.B. Rybicki. #include "nrutil.h" void vander(double x[], double w[], double q[], int n) Solves the Vandermonde linear system N i=1 xk−1 i wi = qk (k = 1,...,N). Input consists of the vectors x[1..n] and q[1..n]; the vector w[1..n] is output. { int i,j,k; double b,s,t,xx; double *c; c=dvector(1,n); if (n == 1) w[1]=q[1]; else { for (i=1;i<=n;i++) c[i]=0.0; Initialize array. c[n] = -x[1]; Coefficients of the master polynomial are found for (i=2;i<=n;i++) { by recursion. xx = -x[i]; for (j=(n+1-i);j<=(n-1);j++) c[j] += xx*c[j+1]; c[n] += xx; } for (i=1;i<=n;i++) { Each subfactor in turn xx=x[i]; t=b=1.0; s=q[n]; for (k=n;k>=2;k--) { is synthetically divided, b=c[k]+xx*b; s += q[k-1]*b; matrix-multiplied by the right-hand side, t=xx*t+b; } w[i]=s/t; and supplied with a denominator. } } free_dvector(c,1,n); } Toeplitz Matrices An N × N Toeplitz matrix is specified by giving 2N − 1 numbers Rk, k = −N + 1,..., −1, 0, 1,...,N − 1. Those numbers are then emplaced as matrix elements constant along the (upper-left to lower-right) diagonals of the matrix: R0 R−1 R−2 ··· R−(N−2) R−(N−1) R1 R0 R−1 ··· R−(N−3) R−(N−2) R2 R1 R0 ··· R−(N−4) R−(N−3) ··· ··· RN−2 RN−3 RN−4 ··· R0 R−1 RN−1 RN−2 RN−3 ··· R1 R0 (2.8.8) The linear Toeplitz problem can thus be written as N j=1 Ri−jxj = yi (i = 1,...,N) (2.8.9) where the xj ’s, j = 1,...,N, are the unknowns to be solved for. The Toeplitz matrix is symmetric if Rk = R−k for all k. Levinson [4] developed an algorithm for fast solution of the symmetric Toeplitz problem, by a bordering method, that is
2.8 Vandermonde Matrices and Toeplitz Matrices 93 a recursive procedure that solves the M-dimensional Toeplitz problem M R-= (i=1,,M0 (2.8.10) j=1 intum for M1,2..until M-N,the desired result,is finally reached.The vector) is the result at the Mth stage,and becomes the desired answer only when N is reached. Levinson's method is well documented in standard texts(e.g.,[5]).The useful fact that the method generalizes to the nonsymmetric case seems to be less well known.At some risk of excessive detail,we therefore give a derivation here,due to G.B.Rybicki. In following a recursion from step M to step M+1 we find that our developing solution (M)changes in this way: ∑R-0= i=1,,M (2.8.11) j=1 becomes M RM)+R-+)=i=1....,M+1 (2.8.12) A8/ j=1 By eliminating y we find ∑R（ C(MD)-2(M+) =R-(M+1) i=1,,M (2.8.13) America Press. = or by letting i M +1-i and jM+1-j, Programs ∑R-G)=R- (2.8.14) SCIENTIFIC j=1 where 6 c0=, (2.8.15) COMPUTING (ISBN To put this another way, ,=牌--G写” j=1,,M (2.8.16) Thus,if we can use recursion to find the order M quantities(M)and G()and the single Numerical 10621 orderM+1 quantity,then all of the other)will follow.Fortunately,the 43106 quantity)follows from equation (28.12)withiM+1. (outside Recipes M RM+1-写+)+RoxH)=M+1 (2.8.17) North Software. j=1 For the unknown orderM+1 quantities we can substitute the previous order quantities in G since c2-,= 0-x+ M+1) (2.8.18) IM+1 The result of this operation is = ∑Ru+1-0-1 (M) (2.8.19) ∑1R+1-,G21-,-R
2.8 Vandermonde Matrices and Toeplitz Matrices 93 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 recursive procedure that solves the M-dimensional Toeplitz problem M j=1 Ri−jx(M) j = yi (i = 1,...,M) (2.8.10) in turn for M = 1, 2,... until M = N, the desired result, is finally reached. The vector x(M) j is the result at the Mth stage, and becomes the desired answer only when N is reached. Levinson’s method is well documented in standard texts (e.g., [5]). The useful fact that the method generalizes to the nonsymmetric case seems to be less well known. At some risk of excessive detail, we therefore give a derivation here, due to G.B. Rybicki. In following a recursion from step M to step M + 1 we find that our developing solution x(M) changes in this way: M j=1 Ri−jx(M) j = yi i = 1,...,M (2.8.11) becomes M j=1 Ri−jx(M+1) j + Ri−(M+1)x(M+1) M+1 = yi i = 1,...,M +1 (2.8.12) By eliminating yi we find M j=1 Ri−j x(M) j − x(M+1) j x(M+1) M+1 = Ri−(M+1) i = 1,...,M (2.8.13) or by letting i → M + 1 − i and j → M + 1 − j, M j=1 Rj−iG(M) j = R−i (2.8.14) where G(M) j ≡ x(M) M+1−j − x(M+1) M+1−j x(M+1) M+1 (2.8.15) To put this another way, x(M+1) M+1−j = x(M) M+1−j − x(M+1) M+1 G(M) j j = 1,...,M (2.8.16) Thus, if we can use recursion to find the order M quantities x(M) and G(M) and the single order M + 1 quantity x(M+1) M+1 , then all of the other x(M+1) j will follow. Fortunately, the quantity x(M+1) M+1 follows from equation (2.8.12) with i = M + 1, M j=1 RM+1−jx(M+1) j + R0x(M+1) M+1 = yM+1 (2.8.17) For the unknown order M + 1 quantities x(M+1) j we can substitute the previous order quantities in G since G(M) M+1−j = x(M) j − x(M+1) j x(M+1) M+1 (2.8.18) The result of this operation is x(M+1) M+1 = M j=1 RM+1−jx(M) j − yM+1 M j=1 RM+1−jG(M) M+1−j − R0 (2.8.19)