Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 27
Текст из файла (страница 27)
. .. .. N −1cNyNxN(2.8.1)Performing the matrix multiplication, you will see that this equation solves for the unknowncoefficients 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 themethod that we are about to describe.2.8 Vandermonde Matrices and Toeplitz Matrices91The 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).The method of solution of both (2.8.1) and (2.8.2) is closely related to Lagrange’spolynomial 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 byPj (x) =Nn=1(n=j)x − xn=Ajk xk−1xj − xnk=1N(2.8.3)Here the meaning of the last equality is to define the components of the matrix Aij as thecoefficients 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 isspecifically designed so that it takes on a value of zero at all xi with i = j, and has a valueof unity at x = xj . In other words,NPj (xi ) = δij =Ajk xk−1i(2.8.4)k=1, whichBut (2.8.4) says that Ajk is exactly the inverse of the matrix of components xk−1iappears 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 =NAjk qk(2.8.5)k=1As for the transpose problem (2.8.1), we can use the fact that the inverse of the transposeis the transpose of the inverse, socj =NAkj yk(2.8.6)k=1The routine in §3.5 implements this.It remains to find a good way of multiplying out the monomial terms in (2.8.3), in orderto get the components of Ajk . This is essentially a bookkeeping problem, and we will let youread the routine itself to see how it can be solved.
One trick is to define a master P (x) byP (x) ≡N(x − xn )(2.8.7)n=1work out its coefficients, and then obtain the numerators and denominators of the specificPj ’s via synthetic division by the one supernumerary term. (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.92Chapter 2.Solution of Linear Algebraic EquationsThe 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)"k−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;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 asNRi−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,2.8 Vandermonde Matrices and Toeplitz Matrices93a recursive procedure that solves the M -dimensional Toeplitz problemM(M )= yiRi−j xj(i = 1, . . . , M )(2.8.10)j=1(M )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:M(M )Ri−j xj= yii = 1, . . . , M(2.8.11)j=1becomesM(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 find (M )(M +1)Mxj − xj= Ri−(M +1)Ri−j(M +1)xM +1j=1or by letting i → M + 1 − i and j → M + 1 − j,M(M )Rj−i Gj= R−i(2.8.14)j=1where(M )(M )≡Gj(M +1)xM +1−j − xM +1−j(2.8.15)(M +1)xM +1To 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,M(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 sincewe can substitute the previous order(M )(M )GM +1−j =xj(M +1)− xj(M +1)xM +1(2.8.18)The result of this operation is"M(M +1)xM +1j=1= "M(M )RM +1−j xj− yM +1(M )j=1 RM +1−j GM +1−j− R0(2.8.19)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 equationsM(M )= yiRj−i zji = 1, . . . , M(2.8.20)j=1Then, the same sequence of operations on this set leads toM(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) that"M(M )− RM +1j=1 RM +1−j Hj(M +1)HM +1 = "M(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 gives"M(M )− R−M −1j=1 Rj−M −1 Gj(M +1)GM +1 = "M(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.2.8 Vandermonde Matrices and Toeplitz Matrices95#include "nrutil.h"#define FREERETURN {free_vector(h,1,n);free_vector(g,1,n);return;}void toeplz(float r[], float x[], float y[], int n)"Solves the Toeplitz system Nj=1 R(N +i−j) xj = yi (i = 1, . . . , N ). The Toeplitz matrix neednot be symmetric. y[1..n] and r[1..2*n-1] are input arrays; x[1..n] is the output array.{int j,k,m,m1,m2;float pp,pt1,pt2,qq,qt1,qt2,sd,sgd,sgn,shn,sxn;float *g,*h;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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
