c2-4 (779462)
Текст из файла
50Chapter 2.Solution of Linear Algebraic EquationsA quick-and-dirty way to solve complex systems is to take the real and imaginaryparts of (2.3.16), givingA·x−C·y=b(2.3.17)C·x+A·y=d(2.3.18)and then solved with ludcmp and lubksb in their present forms. This scheme is a factor of2 inefficient in storage, since A and C are stored twice. It is also a factor of 2 inefficient intime, since the complex multiplies in a complexified version of the routines would each use4 real multiplies, while the solution of a 2N × 2N problem involves 8 times the work ofan N × N one.
If you can tolerate these factor-of-two inefficiencies, then equation (2.3.18)is an easy way to proceed.CITED REFERENCES AND FURTHER READING:Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), Chapter 4.Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.).Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), §3.3, and p. 50.Forsythe, G.E., and Moler, C.B.
1967, Computer Solution of Linear Algebraic Systems (Englewood Cliffs, NJ: Prentice-Hall), Chapters 9, 16, and 18.Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New York: Wiley).Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§4.2.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed.
(New York:McGraw-Hill), §9.11.Horn, R.A., and Johnson, C.R. 1985, Matrix Analysis (Cambridge: Cambridge University Press).2.4 Tridiagonal and Band Diagonal Systemsof EquationsThe special case of a system of linear equations that is tridiagonal, that is, hasnonzero elements only on the diagonal plus or minus one column, is one that occursfrequently. Also common are systems that are band diagonal, with nonzero elementsonly along a few diagonal lines adjacent to the main diagonal (above and below).For tridiagonal sets, the procedures of LU decomposition, forward- and backsubstitution each take only O(N ) operations, and the whole solution can be encodedvery concisely.
The resulting routine tridag is one that we will use in later chapters.Naturally, one does not reserve storage for the full N × N matrix, but only forthe nonzero components, stored as three vectors. The set of equations to be solved is u1r1b 1 c1 0 · · · u2 r2 a2 b 2 c2 · · · ··· · · · · = · · · (2.4.1) uN−1rN−1· · · aN−1 bN−1 cN−1···0aNbNuNrNSample 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).which can be written as a 2N × 2N set of real equations, A −Cxb·=C Ayd2.4 Tridiagonal and Band Diagonal Systems of Equations51Notice that a1 and cN are undefined and are not referenced by the routine that follows.#include "nrutil.h"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 u2trivially 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 beu[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 evenwhen the underlying matrix is nonsingular: A zero pivot can be encountered even fora nonsingular matrix. In practice, this is not something to lose sleep about. The kindsof problems that lead to tridiagonal linear sets usually have additional propertieswhich 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 encountera zero pivot.It is possible to construct special examples in which the lack of pivoting in thealgorithm causes numerical instability.
In practice, however, such instability is almostnever encountered — unlike the general matrix problem where pivoting is essential.The tridiagonal algorithm is the rare case of an algorithm that, in practice, ismore robust than theory says it should be. Of course, should you ever encounter aproblem for which tridag fails, you can instead use the more general method forband diagonal systems, now described (routines bandec and banbks).Some other matrix forms consisting of tridiagonal with a small number ofadditional elements (e.g., upper right and lower left corners) also allow rapidsolution; see §2.7.Band Diagonal SystemsWhere tridiagonal systems have nonzero elements only on the diagonal plus or minusone, band diagonal systems are slightly more general and have (say) m1 ≥ 0 nonzero elementsimmediately to the left of (below) the diagonal and m2 ≥ 0 nonzero elements immediately toits right (above it).
Of course, this is only a useful classification if m1 and m2 are both N .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 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;52Chapter 2.Solution of Linear Algebraic EquationsIn that case, the solution of the linear system by LU decomposition can be accomplishedmuch faster, and in much less storage, than for the general N × N case.The precise definition of a band diagonal matrix with elements aij is thataij = 0j > i + m2whenori > j + m1(2.4.3)3490000112300005657000058930000938200002440000064(2.4.4)which has N = 7, m1 = 2, and m2 = 1, is stored compactly as the 7 × 4 matrix,x xx 49 23 57 93 82 43168344155926x(2.4.5)Here x denotes elements that are wasted space in the compact format; these will not bereferenced by any manipulations and can have arbitrary values. Notice that the diagonalof the original matrix appears in column m1 + 1, with subdiagonal elements to its left,superdiagonal elements to its right.The simplest manipulation of a band diagonal matrix, stored compactly, is to multiplyit by a vector to its right.
Although this is algorithmically trivial, you might want to studythe following routine carefully, as an example of how to pull nonzero elements aij out of thecompact storage format in an orderly fashion.#include "nrutil.h"void banmul(float **a, unsigned long n, int m1, int m2, float x[], float b[])Matrix multiply b = A · x, where A is band diagonal with m1 rows below the diagonal and m2rows above. The input vector x and output vector b are stored as x[1..n] and b[1..n],respectively. The array a[1..n][1..m1+m2+1] stores A as follows: The diagonal elementsare in a[1..n][m1+1].
Subdiagonal elements are in a[j ..n][1..m1] (with j > 1 appropriate to the number of elements on each subdiagonal). Superdiagonal elements are ina[1..j ][m1+2..m1+m2+1] with j < n appropriate to the number of elements on each superdiagonal.{unsigned long i,j,k,tmploop;for (i=1;i<=n;i++) {k=i-m1-1;tmploop=LMIN(m1+m2+1,n-k);b[i]=0.0;for (j=LMAX(1,1-k);j<=tmploop;j++) b[i] += a[i][j]*x[j+k];}}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).Band diagonal matrices are stored and manipulated in a so-called compact form, which resultsif the matrix is tilted 45◦ clockwise, so that its nonzero elements lie in a long, narrowmatrix with m1 + 1 + m2 columns and N rows. This is best illustrated by an example:The band diagonal matrix2.4 Tridiagonal and Band Diagonal Systems of Equations53#include <math.h>#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}#define TINY 1.0e-20void bandec(float **a, unsigned long n, int m1, int m2, float **al,unsigned long indx[], float *d)Given an n × n band diagonal matrix A with m1 subdiagonal rows and m2 superdiagonal rows,compactly stored in the array a[1..n][1..m1+m2+1] as described in the comment for routinebanmul, this routine constructs an LU decomposition of a rowwise permutation of A.
The uppertriangular matrix replaces a, while the lower triangular matrix is returned in al[1..n][1..m1].indx[1..n] is an output vector which records the row permutation effected by the partialpivoting; d is output as ±1 depending on whether the number of row interchanges was evenor odd, respectively. This routine is used in combination with banbks to solve band-diagonalsets of equations.{unsigned long i,j,k,l;int mm;float dum;mm=m1+m2+1;l=m1;for (i=1;i<=m1;i++) {Rearrange the storage a bit.for (j=m1+2-i;j<=mm;j++) a[i][j-l]=a[i][j];l--;for (j=mm-l;j<=mm;j++) a[i][j]=0.0;}*d=1.0;l=m1;for (k=1;k<=n;k++) {For each row...dum=a[k][1];i=k;if (l < n) l++;for (j=k+1;j<=l;j++) {Find the pivot element.if (fabs(a[j][1]) > fabs(dum)) {dum=a[j][1];i=j;}}indx[k]=i;if (dum == 0.0) a[k][1]=TINY;Matrix is algorithmically singular, but proceed anyway with TINY pivot (desirable insome applications).if (i != k) {Interchange rows.*d = -(*d);for (j=1;j<=mm;j++) SWAP(a[k][j],a[i][j])}for (i=k+1;i<=l;i++) {Do the elimination.dum=a[i][1]/a[k][1];al[k][i-k]=dum;for (j=2;j<=mm;j++) a[i][j-1]=a[i][j]-dum*a[k][j];a[i][mm]=0.0;}}}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).It is not possible to store the LU decomposition of a band diagonal matrix A quiteas compactly as the compact form of A itself. The decomposition (essentially by Crout’smethod, see §2.3) produces additional nonzero “fill-ins.” One straightforward storage schemeis to return the upper triangular factor (U ) in the same space that A previously occupied, andto return the lower triangular factor (L) in a separate compact matrix of size N × m1 .
Thediagonal elements of U (whose product, times d = ±1, gives the determinant) are returnedin the first column of A’s storage space.The following routine, bandec, is the band-diagonal analog of ludcmp in §2.3:54Chapter 2.Solution of Linear Algebraic Equations#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}void banbks(float **a, unsigned long n, int m1, int m2, float **al,unsigned long indx[], float b[])Given the arrays a, al, and indx as returned from bandec, and given a right-hand side vectorb[1..n], solves the band diagonal linear equations A · x = b.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
