Nash - Compact Numerical Methods for Computers (523163), страница 27
Текст из файла (страница 27)
This versioncalculates only a positive shift.}shift:=0.0;for i:=1 to n dobegint:=A[i,i];for j:=1 to n doif i<>j then t:=t-abs(A[ij]);if t<shift then shift:=t; {looking for lowest bound to eigenvalue}end; {loop over rows}shift:=-shift; {change sign, since current value < 0 if useful}if shift0.0 then shift:=0.0;writeln(‘Adding a shift of ’,shift,‘ to diagonal of matrix.’);for i:=1 to n dobeginfor j:=1 to n dobeginW[i,j]:=A[i,j]; {copy matrix to working array}if i=j then W[i,i]:=A[i,i]+shift; {adding shift in process}if initev then {initialize eigenvector matrix}beginif i=j then W[i+n,i]:=0.0elsebeginW[i+n,j]:=0.0;end,end; {eigenvector initialization}end; {loop on j}end; {loop on i}NashSVD(n, n, W, Z); {call alg01 to do the work}for i:=1 to n dobeginZ[i]:=sqrt(Z[i])-shift; {to adjust eigenvalues}for j:=1 to n doV[i,j]:=W[n+i,j]; {extract eivenvectors}end; {loop on i}end; {alg13.pas == evsvd}Real symmetric matrices125Example 10.1.
Principal axes of a cubeConsider a cube of uniform density, mass m and edge length a situated so thatthree of its edges which meet at a vertex form the coordinate axes x, y, z. Themoments and products of inertia (see, for instance, Synge and Griffith 1959,pp 282-93) for the cube in this coordinate frame areIxx = I yy = I zz = 2ma 2/3Ixy = I xz = I yz = -ma2/4.These can be formed into the moment-of-inertia tensorwhere all the elements have been measured in units 12/(ma2 ).Algorithm 13 can be used to find the eigensolutions of this matrix. Theeigenvalues then give the principal moments of inertia and the eigenvectors givethe principal axes to which these moments apply. Algorithm 13 when used on aData General NOVA operating in 23-bit binary arithmetic findsI 1 =2ma2/12I 2 =11ma2/12I3 =11ma2/12v l = (0·57735, 0·577351, 0·57735)Tv 2 = (0·707107, –0·707107, –4·33488E–8)Tv 3 = (0·408248, 0·408248, –0·816496)T.(The maximum absolute residual was 3·8147E–6, the maximum inner product4·4226E–7.) The last two principal moments of inertia are the same ordegenerate.
Thus any linear combination of v2 and v 3 will give a new vectorwhich forms a new set of principal axes with v 1 andwhich is orthogonal toI 1 = 2ma2/12I2 =11ma2/12I 3 =11ma 2/12Indeed algorithm 14 on the same system foundv 1 = (0·57735, 0·57735, 0·57735)Tv 2 = (0·732793, –0·678262, –5·45312E–2)Tv 3 = (0·360111, 0·454562, –0·814674)Twith a maximum absolute residual of 7·62939E–6 and a maximum inner product2·38419E–7.126Compact numerical methods for computers10.3.
THE JACOBI ALGORITHM FOR THE EIGENSOLUTIONS OF AREAL SYMMETRIC MATRIXEquation (10.1) immediately leads toV T AV = E(10.31)(using V in place of X). The fact that a real symmetric matrix can be diagonalisedby its eigenvectors gives rise to a number of approaches to the algebraiceigenvalue problem for such matrices. One of the earliest of these was suggestedby Jacobi (1846). This proposes the formation of the sequence of matricesA( 0 ) = AA (k+1) =(V(k) ) T A ( k )V ( k )(10.32)where the V(k) are the plane rotations introduced in §3.3.
The limit of thesequence is a diagonal matrix under some conditions on the angles of rotation.Each rotation is chosen to set one off-diagonal element of the matrix A (k) to zero.In general an element made zero by one rotation will be made non-zero byanother so that a series of sweeps through the off-diagonal elements are needed toreduce the matrix to diagonal form. Note that the rotations in equation (10.32)preserve symmetry, so that there are n(n -1)/2 rotations in one sweep if A is oforder n.Consider now the effect of a single rotation, equation (3.11), in the ij plane.Then for m i, j(10.33)(10.34)while(10.35)(10.36)(10.37)By allowing(10.38)and(10.39)the angle calculation defined by equations (3.22)-(3.27) will causezero.
By lettingto be(10.40)be the measure of the non-diagonal character of A(k) in a similar fashion to thenon-orthogonality measure (3.17), it is straightforward (if a little tedious) to show127Real symmetric matricesthat(10.41)(k+1)(k)so that each rotation causes Ato be ‘more diagonal’ than A .Specification of the order in which the off-diagonal elements are to be madezero defines a particular variant of the Jacobi algorithm. One of the simplest is totake the elements in row-wise fashion: (1, 2), (1, 3), .
. . , (1, n), (2, 3), (2, 4), . . . ,(2, n), . . . , ( n – 1, n ). Such cyclic Jacobi algorithms have not been proved toconverge in general, except in the case where the angle of rotation φ isconstrained so that- π/4< φ < π/ 4 .(10.42)Similarly to the orthogonalisation algorithm of §3.3, the difficulty lies indemonstrating that the ordering of the diagonal elements is stable. For φrestricted as in (10.42), Forsythe and Henrici (1960) have carried out thecomplicated proof, and most authors (Wilkinson 1965, Rutishauser 1966,Schwarz et al 1973, Ralston 1965) discuss algorithms using angle calculationswhich satisfy the restriction (10.42). In fact, among the variety of texts availableon numerical linear algebra, the author has noted only one (Fröberg 1965) whichuses the calculation based on equations (10.38), (10.39) and (3.22)-(3.27).
Theadvantage of the calculation given here is that the eigenvalues, which areproduced as the diagonal elements of the matrix that is the limit of the sequenceA(k), are ordered from most positive to most negative. Most applications whichrequire eigensolutions also require them ordered in some way and the orderingthat the current process yields is the one I have almost invariably been required toproduce.
No extra program code is therefore needed to order the eigenvalues andeigenvectors and there may be some savings in execution time as well.We have already noted in chapter 3 the research of Booker (1985) relating to theconvergence of the cyclic Jacobi method with the ordering angle calculation.Jacobi (1846) did not use a cyclic pattern, but chose to zero the largestoff-diagonal element in the current matrix. This process has generally beenconsidered inappropriate for automatic computation due to the time required forthe search before each plane rotation. However, for comparison purposes Imodified a BASIC precursor of algorithm 14. The changes made were only as many aswere needed to obtain the zeroing of the largest off-diagonal element in thepresent matrix, and no attempt was made to remove computations present inthe algorithm which were used to provide convergence test information. Theprocessor time required for the order-5 Moler matrix (appendix 1) on a DataGeneral NOVA operating in six hexadecimal digit arithmetic was 15·7 secondsfor algorithm 14 while the Jacobi method required 13·7 seconds.
Indeed the latterrequired only 30 rotations while algorithm 14 took 5 sweeps (up to 50 rotations).The comparison of timings on this one example may be misleading, especially asthe system uses an interpreter, that is, programs are executed by translating theBASIC at run time instead of compiling it first. (This has an analogy in thetranslation of a speech either in total or simultaneously.) However, for matriceswith only a few large off-diagonal elements, the relatively simple changes to128Compact numerical methods for computersalgorithm 14 to perform the search are probably worthwhile.
In a larger set oftimings run on both a Data General ECLIPSE and an IBM/370 model 168 in sixhexadecimal digit arithmetic, the ‘original Jacobi method was on average theslowest of five methods tested. (Note that the convergence test at STEP 17 below cannot be used.)10.4. ORGANISATION OF THE JACOBI ALGORITHMTo reduce program code length, the following procedure performs left and rightmultiplications on the matrix A(k) separately, rather than use the formulae(10.33)-(10.37).
This may cause some slight loss in accuracy of the computedresults (compared, for instance, to the procedure jacobi discussed in §10.5).Convergence testing is a relatively complicated matter. In the following algorithm, convergence is assessed by counting the number of rotations skippedduring one sweep of the matrix because the rotation angle is too small to have anyeffect on the current matrix. Rotations needed to order the diagonal elements ofthe matrix (hence the eigenvalues) are always performed. The test that the sine ofthe rotation angle is smaller than the machine precision is used to decide (at STEP10) if the rotation is to be carried out when the diagonal elements are in order.Unfortunately, if two diagonal elements are close in value, then the rotation anglemay not be small even if the off-diagonal element to be set to zero is quite small, sothat it is unnecessary to perform the rotation. Thus at STEP 7, the algorithm tests tosee if the off-diagonal element has a magnitude which when multiplied by 100 isincapable of adjusting the diagonal elements, in which case the rotation is skipped.
Asafety check is provided on the number of sweeps which are carried out since it doesnot seem possible to guarantee the convergence of the algorithm to the satisfaction ofthe above tests, Even if the algorithm terminates after 30 sweeps (my arbitrary choicefor the safety check limit) the approximations to the eigenvalues and eigenvectorsmay still be good, and it is recommended that they be checked by computingresiduals and other tests.Algorithm 14.
A Jacobi algorithm for eigensolutions of a real symmetric matrixProcedure evJacobi(n: integer; {order of matrices}var A,V : matrix; {matrix and eigenvector array}hitev: boolean); {flag to initialize eigenvectorarray to a unit matrix if TRUE}{alg14.pas ==a variant of the Jacobi algorithm for computing eigensolutions of areal symmetric matrixn is the order of the eigenproblemA is the real symmetric matrix in fullV will be rotated by the Jacobi transformations.initev is TRUE if we wish this procedure to initializeV to a unit matrix before commencing; otherwise,V is assumed to be initialized outside the procedure,e.g. for computing the eigenvectors of a generalizedeigenproblem as in alg15.pas.Copyright 1988 J.C.NashReal symmetric matrices129Algorithm 14. A Jacobi algorithm for eigensolutions of a real symmetric matrix (cont.)}{STEP 0 -- via the calling sequence of the procedure, we supply the matrixand its dimensions to the program.}varcount, i, j, k, limit, skipped : integer;c, p, q, s, t : real;ch : char;oki, okj, rotn : boolean;beginwriteln(‘alg14.pas -- eigensolutions of a real symmetric’);writeln(‘matrix via a Jacobi method’);if initev then {Do we initialize the eigenvectors to columns ofthe identity?}beginfor i := l to n dobeginfor j := 1 to n do V[i,j] := 0.0;V[i,i] := 1.0; {to set V to a unit matrix -- rotated to becomethe eigenvectors}end; {loop on i;}end; {initialize eigenvectors}count := 0;limit := 30; {an arbitrary choice following lead of Eberlein}skipped := 0; {so far no rotations have been skipped.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
