c11-3 (779551), страница 2

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

Текст из файла (страница 2)

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).011.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrixn(11.3.12)nThe nth row of this matrix equation isαn qTn−1 + βn qTn = qTn · A(11.3.13)qTn · qm = δnm(11.3.14)Since Q is orthogonal,Thus if we postmultiply equation (11.3.13) by qn , we findβn = qTn · A · qn(11.3.15)which is known since qn is known. Then equation (11.3.13) givesαn qTn−1 = zTn−1(11.3.16)zTn−1 ≡ qTn · A − βn qTn(11.3.17)α2n = zTn−1 zn−1 ,(11.3.18)αn = |zn−1 |(11.3.19)qTn−1 = zTn−1 /αn(11.3.20)whereis known.

Thereforeorand(where αn is nonzero by hypothesis). Similarly, one can show by induction that if we knowqn , qn−1 , . . . , qn−j and the α’s, β’s, and γ’s up to level n − j, one can determine thequantities at level n − (j + 1).To apply the lemma in practice, suppose one can somehow find a tridiagonal matrixAs+1 such thatTAs+1 = Qs · As · Qs(11.3.21)Twhere Qs is orthogonal and has the same last row as QTs in the original QL algorithm.Then Qs = Qs and As+1 = As+1 .Now, in the original algorithm, from equation (11.3.11) we see that the last row of QTs(s)(s)is the same as the last row of Pn−1 .

But recall that Pn−1 is a plane rotation designed toannihilate the (n − 1, n) element of As − ks 1. A simple calculation using the expression(11.1.1) shows that it has parametersd n − ksc= p2en + (dn − ks )2(s)(s)T,−en−1s= p2en + (dn − ks )2The matrix Pn−1 · As · Pn−1 is tridiagonal with 2···× × ×× ××x(11.3.22)extra elements:× x× ×× ×(11.3.23)We must now reduce this to tridiagonal form with an orthogonal matrix whose last row isT(s)[0, 0, . .

. , 0, 1] so that the last row of Qs will stay equal to Pn−1 . This can be done bySample 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).matrix Q.

The relation B · QT = QT · A can be written  qT   qT β 1 γ111 qT2 qT2   α 2 β 2 γ2  .   . ..· .  = . ·A  .   . .  T   T αn−1 βn−1 γn−1qn−1qn−1αnβnqTqT479480Chapter 11.Eigensystemsa sequence of Householder or Givens transformations. For the special form of the matrix(11.3.23), Givens is better.

We rotate in the plane (n − 2, n − 1) to annihilate the (n − 2, n)element. [By symmetry, the (n, n − 2) element will also be zeroed.] This leaves us withtridiagonal form except for extra elements (n − 3, n − 1) and (n − 1, n − 3). We annihilatethese with a rotation in the (n − 3, n − 2) plane, and so on. Thus a sequence of n − 2Givens rotations is required. The result is that(s)(s)(s)(s)(11.3.24)where the P’s are the Givens rotations and Pn−1 is the same plane rotation as in the originalalgorithm. Then equation (11.3.21) gives the next iterate of A. Note that the shift ks entersimplicitly through the parameters (11.3.22).The following routine tqli (“Tridiagonal QL Implicit”), based algorithmicallyon the implementations in [2,3], works extremely well in practice.

The numberof iterations for the first few eigenvalues might be 4 or 5, say, but meanwhilethe off-diagonal elements in the lower right-hand corner have been reduced too.The later eigenvalues are liberated with very little work. The average number ofiterations per eigenvalue is typically 1.3 − 1.6. The operation count per iteration isO(n), with a fairly large effective coefficient, say, ∼ 20n. The total operation countfor the diagonalization is then ∼ 20n × (1.3 − 1.6)n ∼ 30n2 . If the eigenvectorsare required, the statements indicated by comments are included and there is anadditional, much larger, workload of about 3n3 operations.#include <math.h>#include "nrutil.h"void tqli(float d[], float e[], int n, float **z)QL algorithm with implicit shifts, to determine the eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix, or of a real, symmetric matrix previously reduced by tred2 §11.2.

Oninput, d[1..n] contains the diagonal elements of the tridiagonal matrix. On output, it returnsthe eigenvalues. The vector e[1..n] inputs the subdiagonal elements of the tridiagonal matrix,with e[1] arbitrary. On output e is destroyed. When finding only the eigenvalues, several linesmay be omitted, as noted in the comments.

If the eigenvectors of a tridiagonal matrix are desired, the matrix z[1..n][1..n] is input as the identity matrix. If the eigenvectors of a matrixthat has been reduced by tred2 are required, then z is input as the matrix output by tred2.In either case, the kth column of z returns the normalized eigenvector corresponding to d[k].{float pythag(float a, float b);int m,l,iter,i,k;float s,r,p,g,f,dd,c,b;for (i=2;i<=n;i++) e[i-1]=e[i];Convenient to renumber the ele[n]=0.0;ements of e.for (l=1;l<=n;l++) {iter=0;do {for (m=l;m<=n-1;m++) {Look for a single small subdidd=fabs(d[m])+fabs(d[m+1]);agonal element to splitif ((float)(fabs(e[m])+dd) == dd) break;the matrix.}if (m != l) {if (iter++ == 30) nrerror("Too many iterations in tqli");g=(d[l+1]-d[l])/(2.0*e[l]);Form shift.r=pythag(g,1.0);g=d[m]-d[l]+e[l]/(g+SIGN(r,g));This is dm − ks .s=c=1.0;p=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).TQTs = Qs = P1 · P2 · · · Pn−2 · Pn−111.4 Hermitian Matrices481}} while (m != l);}}CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 331–335. [1]Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2]Smith, B.T., et al.

1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [3]Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§6.6.6. [4]11.4 Hermitian MatricesThe complex analog of a real, symmetric matrix is a Hermitian matrix,satisfying equation (11.0.4). Jacobi transformations can be used to find eigenvaluesand eigenvectors, as also can Householder reduction to tridiagonal form followed byQL iteration.

Complex versions of the previous routines jacobi, tred2, and tqliare quite analogous to their real counterparts. For working routines, consult [1,2] .An alternative, using the routines in this book, is to convert the Hermitianproblem to a real, symmetric one: If C = A + iB is a Hermitian matrix, then then × n complex eigenvalue problem(A + iB) · (u + iv) = λ(u + iv)(11.4.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).for (i=m-1;i>=l;i--) {A plane rotation as in the origif=s*e[i];nal QL, followed by Givensb=c*e[i];rotations to restore tridiage[i+1]=(r=pythag(f,g));onal form.if (r == 0.0) {Recover from underflow.d[i+1] -= p;e[m]=0.0;break;}s=f/r;c=g/r;g=d[i+1]-p;r=(d[i]-g)*s+2.0*c*b;d[i+1]=g+(p=s*r);g=c*r-b;/* Next loop can be omitted if eigenvectors not wanted*/for (k=1;k<=n;k++) {Form eigenvectors.f=z[k][i+1];z[k][i+1]=s*z[k][i]+c*f;z[k][i]=c*z[k][i]-s*f;}}if (r == 0.0 && i >= l) continue;d[l] -= p;e[l]=g;e[m]=0.0;.

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

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