c2-9 (779467), страница 2
Текст из файла (страница 2)
(See §5.3 for more on syntheticdivision.) Since each such division is only a process of order N , the total procedure isof order N 2 .You should be warned that Vandermonde systems are notoriously ill-conditioned, bytheir very nature. (As an aside anticipating §5.8, the reason is the same as that which makesChebyshev fitting so impressively accurate: there exist high-order polynomials that are verygood uniform fits to zero.
Hence roundoff error can introduce rather substantial coefficientsof the leading terms of these polynomials.) It is a good idea always to compute Vandermondeproblems in double precision.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The alternative identification of rows and columns leads to the set of equations 11···1w1q1 x1x2···xN w2 q2 2 x22···x2N · w3 = q3 (2.8.2) x1 ··· ······−1−1−1xNwNqNxN· · · xN12NWrite this out and you will see that it relates to the problem of moments: Given the valuesof N points xi , find the unknown weights wi , assigned so as to match the given valuesqj of the first N moments. (For more on this problem, consult [3].) The routine given inthis section solves (2.8.2).92Chapter 2.Solution of Linear Algebraic EquationsThe routine for (2.8.2) which follows is due to G.B.
Rybicki.#include "nrutil.h"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 foundfor (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 turnxx=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 MatricesAn 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 constantalong the (upper-left to lower-right) diagonals of the matrix: R0 R1 R2 ···RN −2RN −1R−1R0R1R−2R−1R0RN −3RN −2RN −4RN −3··················R−(N −2)R−(N −3)R−(N −4)R0R1R−(N −1) R−(N −2) R−(N −3) R−1R0(2.8.8)The linear Toeplitz problem can thus be written asNXRi−j xj = yi(i = 1, .
. . , N )(2.8.9)j=1where 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 analgorithm for fast solution of the symmetric Toeplitz problem, by a bordering method, that is,Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void vander(double x[], double w[], double q[], int n)Pk−1Solves the Vandermonde linear system Nwi = qk (k = 1, . . . , N ). Input consists ofi=1 xithe 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;2.8 Vandermonde Matrices and Toeplitz Matrices93a recursive procedure that solves the M -dimensional Toeplitz problemMX(M )Ri−j xj= yi(i = 1, . . .
, M )(2.8.10)j=1(M )MX(M )Ri−j xj= yii = 1, . . . , M(2.8.11)j=1becomesMX(M +1)Ri−j xj(M +1)+ Ri−(M +1) xM +1 = yii = 1, . . . , M + 1(2.8.12)i = 1, . . . , M(2.8.13)j=1By eliminating yi we findMX(M )Ri−jxj(M +1)− xj!= Ri−(M +1)(M +1)xM +1j=1or by letting i → M + 1 − i and j → M + 1 − j,MX(M )Rj−i Gj= R−i(2.8.14)j=1where(M )(M )≡Gj(M +1)xM +1−j − xM +1−j(M +1)xM +1(2.8.15)To put this another way,(M +1)(M )(M +1)(M )xM +1−j = xM +1−j − xM +1 Gjj = 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(M +1)(M +1)order M + 1 quantity xM +1 , then all of the other xjwill follow.
Fortunately, the(M +1)quantity xM +1follows from equation (2.8.12) with i = M + 1,MX(M +1)RM +1−j xj(M +1)+ R0 xM +1 = yM +1(2.8.17)j=1(M +1)For the unknown order M + 1 quantities xjquantities in G since(M )(M )GM +1−j =xjwe can substitute the previous order(M +1)− xj(2.8.18)(M +1)xM +1The result of this operation isPM(M +1)j=1xM +1 = PM(M )RM +1−j xj(M )− yM +1j=1 RM +1−j GM +1−j − R0(2.8.19)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).in turn for M = 1, 2, . . . until M = N , the desired result, is finally reached. The vector xjis the result at the M th 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 thatthe method generalizes to the nonsymmetric case seems to be less well known. At some riskof 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 solutionx(M ) changes in this way:94Chapter 2.Solution of Linear Algebraic EquationsThe only remaining problem is to develop a recursion relation for G. Before we dothat, however, we should point out that there are actually two distinct sets of solutions to theoriginal linear problem for a nonsymmetric matrix, namely right-hand solutions (which wehave been discussing) and left-hand solutions zi .
The formalism for the left-hand solutionsdiffers only in that we deal with the equationsMX(M )= yii = 1, . . . , M(2.8.20)Then, the same sequence of operations on this set leads toMX(M )Ri−j Hj= Ri(2.8.21)j=1where(M )(M )Hj≡(M +1)zM +1−j − zM +1−j(2.8.22)(M +1)zM +1(compare with 2.8.14 – 2.8.15). The reason for mentioning the left-hand solutions now isthat, by equation (2.8.21), the Hj satisfy exactly the same equation as the xj except forthe substitution yi → Ri on the right-hand side.
Therefore we can quickly deduce fromequation (2.8.19) thatPM(M )− RM +1j=1 RM +1−j Hj(M +1)HM +1 = PM(2.8.23)(M )j=1 RM +1−j GM +1−j − R0By the same token, G satisfies the same equation as z, except for the substitution yi → R−i .This givesPM(M )− R−M −1j=1 Rj−M −1 Gj(M +1)GM +1 = PM(2.8.24)(M )j=1 Rj−M −1 HM +1−j − R0The same “morphism” also turns equation (2.8.16), and its partner for z, into the final equations(M +1)= Gj(M +1)= HjGjHj(M )− GM +1 HM +1−j(M +1)(M )− HM +1 GM +1−j(M +1)(M )(M )(2.8.25)Now, starting with the initial values(1)x1 = y1 /R0(1)G1 = R−1 /R0(1)H1= R1 /R0(2.8.26)we can recurse away. At each stage M we use equations (2.8.23) and (2.8.24) to find(M +1)(M +1)HM +1 , GM +1 , and then equation (2.8.25) to find the other components of H (M +1), G(M +1) .From there the vectors x(M +1) and/or z (M +1) are easily calculated.The program below does this.
It incorporates the second equation in (2.8.25) in the form(M +1)(M )(M +1)(M )HM +1−j = HM +1−j − HM +1 Gj(2.8.27)so that the computation can be done “in place.”Notice that the above algorithm fails if R0 = 0. In fact, because the bordering methoddoes not allow pivoting, the algorithm will fail if any of the diagonal principal minors of theoriginal Toeplitz matrix vanish. (Compare with discussion of the tridiagonal algorithm in§2.4.) If the algorithm fails, your matrix is not necessarily singular — you might just haveto solve your problem by a slower and more general algorithm such as LU decompositionwith pivoting.The routine that implements equations (2.8.23)–(2.8.27) is also due to Rybicki.
Notethat the routine’s r[n+j] is equal to Rj above, so that subscripts on the r array vary from1 to 2N − 1.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Rj−i zjj=12.8 Vandermonde Matrices and Toeplitz Matrices95#include "nrutil.h"#define FREERETURN {free_vector(h,1,n);free_vector(g,1,n);return;}if (r[n] == 0.0) nrerror("toeplz-1 singular principal minor");g=vector(1,n);h=vector(1,n);x[1]=y[1]/r[n];Initialize for the recursion.if (n == 1) FREERETURNg[1]=r[n-1]/r[n];h[1]=r[n+1]/r[n];for (m=1;m<=n;m++) {Main loop over the recursion.m1=m+1;sxn = -y[m1];Compute numerator and denominator for x,sd = -r[n];for (j=1;j<=m;j++) {sxn += r[n+m1-j]*x[j];sd += r[n+m1-j]*g[m-j+1];}if (sd == 0.0) nrerror("toeplz-2 singular principal minor");x[m1]=sxn/sd;whence x.for (j=1;j<=m;j++) x[j] -= x[m1]*g[m-j+1];if (m1 == n) FREERETURNsgn = -r[n-m1];Compute numerator and denominator for G and H,shn = -r[n+m1];sgd = -r[n];for (j=1;j<=m;j++) {sgn += r[n+j-m1]*g[j];shn += r[n+m1-j]*h[j];sgd += r[n+j-m1]*h[m-j+1];}if (sd == 0.0 || sgd == 0.0) nrerror("toeplz-3 singular principal minor");g[m1]=sgn/sgd;whence G and H.h[m1]=shn/sd;k=m;m2=(m+1) >> 1;pp=g[m1];qq=h[m1];for (j=1;j<=m2;j++) {pt1=g[j];pt2=g[k];qt1=h[j];qt2=h[k];g[j]=pt1-pp*qt2;g[k]=pt2-pp*qt1;h[j]=qt1-qq*pt2;h[k--]=qt2-qq*pt1;}}Back for another recurrence.nrerror("toeplz - should not arrive here!");}If you are in the business of solving very large Toeplitz systems, you should find out aboutso-called “new, fast” algorithms, which require only on the order of N (log N )2 operations,compared to N 2 for Levinson’s method.
These methods are too complicated to include here.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void toeplz(float r[], float x[], float y[], int n)PSolves the Toeplitz system Nj=1 R(N +i−j) xj = yi (i = 1, .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
