c11-3 (779551)
Текст из файла
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix475}} elsee[i]=a[i][l];d[i]=h;}/* Next statement can be omitted if eigenvectors not wanted */d[1]=0.0;e[1]=0.0;/* Contents of this loop can be omitted if eigenvectors notwanted except for statement d[i]=a[i][i]; */for (i=1;i<=n;i++) {Begin accumulation of transformation mal=i-1;trices.if (d[i]) {This block skipped when i=1.for (j=1;j<=l;j++) {g=0.0;for (k=1;k<=l;k++)Use u and u/H stored in a to form P·Q.g += a[i][k]*a[k][j];for (k=1;k<=l;k++)a[k][j] -= g*a[k][i];}}d[i]=a[i][i];This statement remains.a[i][i]=1.0;Reset row and column of a to identityfor (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;matrix for next iteration.}}CITED REFERENCES AND FURTHER READING:Golub, G.H., and Van Loan, C.F.
1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §5.1. [1]Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag).Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag).
[2]11.3 Eigenvalues and Eigenvectors of aTridiagonal MatrixEvaluation of the Characteristic PolynomialOnce our original, real, symmetric matrix has been reduced to tridiagonal form,one possible way to determine its eigenvalues is to find the roots of the characteristicpolynomial pn (λ) directly. The characteristic polynomial of a tridiagonal matrix canbe evaluated for any trial value of λ by an efficient recursion relation (see [1], forexample). The polynomials of lower degree produced during the recurrence form aSample 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).f += e[j]*a[i][j];}hh=f/(h+h);Form K, equation (11.2.11).for (j=1;j<=l;j++) {Form q and store in e overwriting p.f=a[i][j];e[j]=g=e[j]-hh*f;for (k=1;k<=j;k++)Reduce a, equation (11.2.13).a[j][k] -= (f*e[k]+g*a[i][k]);}476Chapter 11.EigensystemsThe QR and QL AlgorithmsThe basic idea behind the QR algorithm is that any real matrix can bedecomposed in the formA= Q·R(11.3.1)where Q is orthogonal and R is upper triangular.
For a general matrix, thedecomposition is constructed by applying Householder transformations to annihilatesuccessive columns of A below the diagonal (see §2.10).Now consider the matrix formed by writing the factors in (11.3.1) in theopposite order:A0 = R · Q(11.3.2)Since Q is orthogonal, equation (11.3.1) gives R = QT · A. Thus equation (11.3.2)becomesA0 = QT · A · Q(11.3.3)We see that A0 is an orthogonal transformation of A.You can verify that a QR transformation preserves the following properties ofa matrix: symmetry, tridiagonal form, and Hessenberg form (to be defined in §11.5).There is nothing special about choosing one of the factors of A to be uppertriangular; one could equally well make it lower triangular.
This is called the QLalgorithm, sinceA= Q·L(11.3.4)where L is lower triangular. (The standard, but confusing, nomenclature R and Lstands for whether the right or left of the matrix is nonzero.)Recall that in the Householder reduction to tridiagonal form in §11.2, we startedin the nth (last) column of the original matrix. To minimize roundoff, we thenexhorted you to put the biggest elements of the matrix in the lower right-handcorner, if you can.
If we now wish to diagonalize the resulting tridiagonal matrix,the QL algorithm will have smaller roundoff than the QR algorithm, so we shalluse QL henceforth.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).Sturmian sequence that can be used to localize the eigenvalues to intervals on thereal axis. A root-finding method such as bisection or Newton’s method can thenbe employed to refine the intervals.
The corresponding eigenvectors can then befound by inverse iteration (see §11.7).Procedures based on these ideas can be found in [2,3] . If, however, morethan a small fraction of all the eigenvalues and eigenvectors are required, then thefactorization method next considered is much more efficient.47711.3 Eigenvalues and Eigenvectors of a Tridiagonal MatrixThe QL algorithm consists of a sequence of orthogonal transformations:As = Qs · LsAs+1 = Ls · Qs(= QTs · As · Qs )(11.3.5)(s)aij∼λiλjs(11.3.6)Although λi < λj , convergence can be slow if λi is close to λj .
Convergence canbe accelerated by the technique of shifting: If k is any constant, then A − k1 haseigenvalues λi − k. If we decomposeAs − ks 1 = Qs · Lsso that(11.3.7)As+1 = Ls · Qs + ks 1= QTs · As · Qs(11.3.8)then the convergence is determined by the ratioλi − ksλj − ks(11.3.9)The idea is to choose the shift ks at each stage to maximize the rate ofconvergence. A good choice for the shift initially would be ks close to λ1 , thesmallest eigenvalue. Then the first row of off-diagonal elements would tend rapidlyto zero. However, λ1 is not usually known a priori. A very effective strategy inpractice (although there is no proof that it is optimal) is to compute the eigenvaluesSample 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 following (nonobvious!) theorem is the basis of the algorithm for a generalmatrix A: (i) If A has eigenvalues of different absolute value |λi |, then As → [lowertriangular form] as s → ∞.
The eigenvalues appear on the diagonal in increasingorder of absolute magnitude. (ii) If A has an eigenvalue |λi | of multiplicity p,As → [lower triangular form] as s → ∞, except for a diagonal block matrixof order p, whose eigenvalues → λi . The proof of this theorem is fairly lengthy;see, for example, [4].The workload in the QL algorithm is O(n3 ) per iteration for a general matrix,which is prohibitive.
However, the workload is only O(n) per iteration for atridiagonal matrix and O(n2 ) for a Hessenberg matrix, which makes it highlyefficient on these forms.In this section we are concerned only with the case where A is a real, symmetric,tridiagonal matrix. All the eigenvalues λi are thus real. According to the theorem,if any λi has a multiplicity p, then there must be at least p − 1 zeros on thesub- and superdiagonal.
Thus the matrix can be split into submatrices that can bediagonalized separately, and the complication of diagonal blocks that can arise inthe general case is irrelevant.In the proof of the theorem quoted above, one finds that in general a superdiagonal element converges to zero like478Chapter 11.Eigensystemsof the leading 2 × 2 diagonal submatrix of A. Then set ks equal to the eigenvaluecloser to a11 .More generally, suppose you have already found r − 1 eigenvalues of A. Thenyou can deflate the matrix by crossing out the first r − 1 rows and columns, leaving. ..A=. ..0·········drererdr+100······0dn−1en−1.. .
0 en−1 dn(11.3.10)Choose ks equal to the eigenvalue of the leading 2 × 2 submatrix that is closer to dr .One can show that the convergence of the algorithm with this strategy is generallycubic (and at worst quadratic for degenerate eigenvalues). This rapid convergenceis what makes the algorithm so attractive.Note that with shifting, the eigenvalues no longer necessarily appear on thediagonal in order of increasing absolute magnitude. The routine eigsrt (§11.1)can be used if required.As we mentioned earlier, the QL decomposition of a general matrix is effectedby a sequence of Householder transformations. For a tridiagonal matrix,however, it ismore efficient to use plane rotations Ppq .
One uses the sequence P12 , P23 , . . . , Pn−1,nto annihilate the elements a12 , a23, . . . , an−1,n. By symmetry, the subdiagonalelements a21 , a32 , . . . , an,n−1 will be annihilated too. Thus each Qs is a productof plane rotations:(s)(s)(s)QTs = P1 · P2 · · · Pn−1(11.3.11)where Pi annihilates ai,i+1 . Note that it is QT in equation (11.3.11), not Q, becausewe defined L = QT · A.QL Algorithm with Implicit ShiftsThe algorithm as described so far can be very successful. However, whenthe elements of A differ widely in order of magnitude, subtracting a large ksfrom the diagonal elements can lead to loss of accuracy for the small eigenvalues.This difficulty is avoided by the QL algorithm with implicit shifts.
The implicitQL algorithm is mathematically equivalent to the original QL algorithm, but thecomputation does not require ks 1 to be actually subtracted from A.The algorithm is based on the following lemma: If A is a symmetric nonsingular matrixand B = QT · A · Q, where Q is orthogonal and B is tridiagonal with positive off-diagonalelements, then Q and B are fully determined when the last row of QT is specified. Proof:Let qTi denote the ith row vector of the matrix QT . Then qi is the ith column vector of theSample 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-файл и есть ли нужная программа для его просмотра.














