c2-8 (779466)

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла

90Chapter 2.Solution of Linear Algebraic Equations2.8 Vandermonde Matrices and ToeplitzMatricesVandermonde MatricesA Vandermonde matrix of size N × N is completely determined by N arbitrarynumbers x1 , x2 , . . . , xN , in terms of which its N 2 components are the integer powersxj−1, i, j = 1, . . . , N . Evidently there are two possible such forms, depending on whetheriwe view the i’s as rows, j’s as columns, or vice versa. In the former case, we get a linearsystem of equations that looks like this,11. ..1x1x21···x2...xNx22...x2N······−1xN1  −1  xN2·.. .  −1xNNc1c2...cN  =  y1y2 .. .

yN(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.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 §2.4 the case of a tridiagonal matrix was treated specially, because thatparticular 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 typesexist, it is important to know about them. Your computational savings, should youever happen to be working on a problem that involves the right kind of particulartype, can be enormous.This section treats two special types of matrices that can be solved in of orderN 2 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 problemshaving to do with the fitting of polynomials, the reconstruction of distributions fromtheir moments, and also other contexts.

In this book, for example, a Vandermondeproblem crops up in §3.5. Matrices of the second type, termed Toeplitz matrices,tend to occur in problems involving deconvolution and signal processing. In thisbook, a Toeplitz problem is encountered in §13.7.These are not the only special types of matrices worth knowing about. TheHilbert 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 toinvert 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 sometimesbe used to convert new special forms into old ones. Reference [2] gives some otherspecial forms.

We have not found these additional forms to arise as frequently asthe two that we now discuss.2.8 Vandermonde Matrices and Toeplitz Matrices91The 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) =NYn=1(n6=j)Xx − 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 6= j, and has a valueof unity at x = xj .

In other words,NXPj (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 =NXAjk 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 =NXAkj 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) ≡NY(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.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.

Характеристики

Тип файла PDF

PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.

Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.

Список файлов книги