Nash - Compact Numerical Methods for Computers (523163), страница 28
Текст из файла (страница 28)
We need to setskipped here because the while loop below tests this variable.}{main loop}while (count<=limit) and (skipped<((n*(n-1)) div 2) ) do{This is a safety check to avoid indefinite execution of the algorithm.The figure used for limit here is arbitrary. If the program terminatesby exceeding the sweep limit, the eigenvalues and eigenvectors computedmay still be good to the limitations of the machine in use, thoughresiduals should be calculated and other tests made.
Such tests arealways useful, though they are time- and space-consuming.}begincount := count+1; {to count sweeps -- STEP 1}write(‘sweep’,count,‘ ’);skipped := 0; {to count rotations skipped during the sweep.}for i := 1 to (n-1) do {STEP 2}begin {STEP 3}for j := (i+1) to n do {STEP 4}beginrotn := true; {to indicate that we carry out a rotation unlesscalculations show it unnecessary}p := 0.5*(A[i,j]+A[j,i]); {An average of the off-diagonal elementsis used because the right and left multiplications by therotation matrices are non-commutative in most cases due todifferences between floating-point and exact arithmetic.}q := A[i,i]-A[j,j]; {Note: this may suffer from digit cancellationwhen nearly equal eigenvalues exist. This cancellation is notnormally a problem, but may cause the algorithm to perform morework than necessary when the off-diagonal elements are verysmall.}130Compact numerical methods for computersAlgorithm 14. A Jacobi algorithm for eigensolutions of a real symmetric matrix (cont.)t := sqrt(4.0*p*p+q*q);if t=0.0 then {STEP 5}begin {STEP 11 -- If t is zero, no rotation is needed.}rotn := false; {to indicate no rotation is to be performed.}endelsebegin {t>0.0}if q>=0.0 then {STEP 6}begin {rotation for eigenvalue approximations already in order}{STEP 7 -- test for small rotation}oki := (abs(A[i,i])=abs(A[i,i])+l00.0*abs(p));okj := (abs(A[j,j])=abs(A[j,j])+l00.0*abs(p));if oki and okj then rotn := falseelse rotn := true;{This test for a small rotation uses an arbitrary factor of100 for scaling the off-diagonal elements.
It is chosen toensure “small but not very small” rotations are performed.}if rotn thenbegin {STEP 8}c := sqrt((t+q)/(2.0*t)); s := p/(t*c);end;end {if q>=0.0}elsebegin {q<0.0 -- always rotate to bring eigenvalues into order}rotn := true; {STEP 9}s := sqrt((t-q)/(2,0*t));if p<0.0 then s := -s;c := p/(t*s);end; {STEP 10}if 1.0=(1.0+abs(s)) then rotn := false; {test for small angle}end; {if t=0.0}if rotn then {STEP 11 -- rotate if necessary}begin {STEP 12}for k := 1 to n dobeginq := A[i,k]; A[i,k] := c*q+s*A[j,k]; A[j,k] := -s*q+c*A[j,k];end; {left multiplication of matrix A}{STEP 13}for k := l to n dobegin {right multiplication of A and V}q := A[k,i]; A[k,i] := c*q+s*A[k,j]; A[k,j] := -s*q+c*A[k,j];{STEP 14 -- can be omitted if eigenvectors not needed}q := V[k,i]; V[k,i] := c*q+s*V[k,j]; V[k,j] := -s*q+c*V[k,j];end; {loop on k for right multiplication of matrices A and V}end {rotation carried out}else{STEP 11 -- count skipped rotations}skipped := skipped+1; {to count the skipped rotations}end; {loop on j} {STEP 15}end; {loop on i.
This is also the end of the sweep. -- STEP 16}writeln(‘ ’,skipped,‘ / ’,n*(n-1) div 2,‘ rotations skipped’);end; {while -- main loop}end; {alg14.pas = evJacobi -- STEP 17}Real symmetric matrices131Example 10.2. Application of the Jacobi algorithm in celestial mechanicsIt is appropriate to illustrate the use of algorithm 14 by Jacobi’s (1846) ownexample. This arises in the study of orbital perturbations of the planets tocompute corrections to some of the parameters of the solar system. Unfortunatelyat the time Jacobi was writing his paper, Neptune had not been discovered.Leverrier reported calculations suggesting the existence of this planet on 31August 1846, to l’Académie des Sciences in Paris, and Galle in Berlin confirmedthis hypothesis by observation less than three weeks later on 18 September.
Thederivation of the eigenproblem in this particular case is lengthy and irrelevant tothe present illustration, so we will begin with Jacobi’s equations V. These give anon-symmetric matrix à which can be symmetrised by a diagonal transformationresulting in Jacobi’s equations VIII, where the off-diagonal elements are expressed in terms of their common logarithms to which 10 has been added. I decided towork with the non-symmetric form and symmetrised it by means ofAij = Aji = (à ij à ji ) ½ .The output from a Hewlett-Packard 9830 (machine precision = 1E–11) is givenbelow, and includes the logarithmic elements which in every case approximatedvery closely Jacobi’s equations VIII. For comparison, he computed eigenvalues–2·2584562,–3·7151584,–5·2986987,–7·5740431,–17·1524687,–17·8632192 and –22·4267712 after 10 rotations. At this point the largestoff-diagonal element (which is marked as being negative) had a logarithm(8·8528628–10), which has the approximate antilog –7·124E–2.
Jacobi used as acomputing system his student Ludwig Seidel, apparently operating in eight-digitdecimal arithmetic!ENHJCB JACOBI WITH ORDERING MAR 5 75ORDER= 7INPUT JACOBI'S MATRIX:ROW 1-5.5098820.4229081.8700868.81400E-030.1487113.90800E-034.50000E-05:ROW 20.287865-11.8116545.71190.0587170.7280880.0187882.24000E-04ROW 3 :0.0490994.308033-12.9706871.6890870.2293260.042585.04OOOE-04ROW 4 :6.23500E-030.2698511.397369-17.5962075.3040380.1253461.45100E-03ROW 5 :2.23100E-057.09480E-042.18227E-03l.l2462E-03-7.4890414.8154540.035319ROW 6:1.45000E-064.52200E-051.35880E-046.56500E-0511.893979-18.585410.232241ROW 7 :6.00000E-081.94000E-065.79000E-060.3138292.73000E-060.835482-2.325935SYMMETRIZE A(I,J)=A(J,I):=SQR(A(I,J)*A(J,I))LOG(S) GIVEN FOR COMPARISON WITH JACOBI'S TABLE VIIIS=A( 1,2)= 0.733711324LOG10(S)+l0= 9.865525222S=A( 1LOG10(S)+l0= 9.158659274)= 0.144098438,3132S=A(s=A(S=A(S=A(S=A(S=A(S=A(S=A(S=A(S=A(S=A(S=A(S=A(s=A(S=A(S=A(S=A(S=A(S=A(Compact numerical methods for computers1111222223333444556)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=,4, 5,6,7, 3, 4,5, 6,7,4, 5, 6, 7,5, 6,7, 6,7, 77.41318E-031.82147E-037.52768E-051.64317E-064.9605497370.1258762930.0227280429.21734E-042.08461E-050.5660857210.0607127982.40536E-035.40200E-050.0772335892.86862E-036.29383E-057.5680188130.1052811780.440491969LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=LOG10(S)+10=7.8700047527.2604213325.8766612794.21568188210.69552989.0999439458.3565620156.9646055555.3190248749.75288228.783280257.3811805985.732554558.8878062157.4576726055.79891503010.87898229.0223507369.643937995MATRIXROW 1:-5.509882 0.733711324 0.144098438 7.41318E-03 1.82147E-031.64317E-067.52768E-05:ROW 20.733711324 -11.811654 4.960549737 0.125876293 0.0227280429.21734E-04 2.08461E-05:ROW 30.144098438 4.960549737 -12.970687 0.566085721 0.0607127982.40536E-03 5.40200E-05:ROW 47.41318E-03 0.125876293 9.566085721 -17.5962070.0772335892.86362E-03 6.29383E-05:ROW 51.82147E-030.022728042 0.0607127980.077233589-7.4890417.568018813 0.105281178:ROW 67.52768E-05 9.21734E-04 2.40536E-032.86862E-037.568018813-18.585410.440491969ROW 7:1.64317E-06 2.08461E-05 5.40200E-056.29383E-050.1052811780.440491969 -2.325935NASH JACOBI ALG.0ROTATIONSROTATIONS00ROTATIONS5ROTATIONS19ROTATIONS21ROTATIONSCONVERGEDEIGENVALUE 14.90537E-040.107934624EIGENVALUE 26.13203E-030.438925002EIGENVALUE 3-0.9548354055.12620E-03EIGENVALUE 40.2953047477.55197E-03EIGENVALUE 5-0.0277000006.88542E-0314 DEC 13/77SKIPPEDSKIPPEDSKIPPEDSKIPPEDSKIPPEDSKIPPEDVECTOR:=-2.2584175961.72184E-031.37576E-030.978469440VECTOR:=-3.7136435880.0118612380.010411191-0.205670808VECTOR:=-5.298872615-0.240328992-0.174066634-1.07833E-03=-7.574719192VECTOR:-0.704148469-0.644023478-8.51399E-04VECTOR:=-17.15255764-0.5841833720.560202581-2.12488E-049.85037E-040.1759022565.53618E-030.874486723-0.0109897469.16066E-03-0.0449156260.011341429-0.5865879051.65695E-03Real symmetric matricesEIGENVALUE 60.0167129514.72676E-03EIGENVALUE 73.56838E-050.891931555=-17.86329687-0.363977437-1.36341E-04=-22.42730849-3.42339E-04-0.017179128133VECTOR:0.462076406-0.8085194486.627116-04VECTOR:2.46894E-036.41285E-03-0.451791590MAXIMUM ABSOLUTE RESIDUAL= 9.51794E-10MAXIMUM ABSOLUTE INNER PRODUCT= 1.11326E-1110.5.
A BRIEF COMPARISON OF METHODS FOR THEEIGENPROBLEM OF A REAL SYMMETRIC MATRIXThe programmer who wishes to solve the eigenproblem of a real symmetricmatrix must first choose which algorithm he will use. The literature presents a vastarray of methods of which Wilkinson (1965) and Wilkinson and Reinsch (1971)give a fairly large sample. The method that is chosen will depend on the size ofthe matrix, its structure and whether or not eigenvectors are required.
Supposethat all eigenvalues and eigenvectors are wanted for a matrix whose order is lessthan 20 with no particular structure. There are two families of algorithm whichare candidates to solve this problem. The first reduces the square matrix to atridiagonal form of which the eigenproblem is then solved by means of one ofseveral variants of either the QR algorithm or bisection. (Bisection will find onlyeigenvalues; inverse iteration is used to obtain the eigenvectors.) The eigenvectorsof the tridiagonal matrix must then be back-transformed to find those of theoriginal matrix. Because this family of methods—tridiagonalisation, solution ofreduced form, back-transformation—requires fairly complicated codes from thestandpoint of this work, they will not be discussed in detail.
For matrices of ordergreater than 10, however, the Householder tridiagonalisation with the QL algorithm as described by Wilkinson and Reinsch (1971) is probably the mostefficient method for the solution of the complete eigenproblem of a real symmetric matrix. (Some studies by Dr Maurice Cox of the UK National PhysicalLaboratory show that Givens’ tridiagonalisation, if carefully coded, usually involves less work than that of Householder.) It is probably worthwhile for the userwith many eigenproblems to solve of order greater than 10 to implement such amethod.The other family of methods is based on the Jacobi algorithm already discussed.Wilkinson, in Wilkinson and Reinsch (1971, pp 192-3), describes Rutishauser’svariant of this, called jacobi :‘The method of Jacobi is the most elegant of those developed for solving thecomplete eigenproblem.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
