MacKinnon - Computational Physics, страница 8

Naı̈vely we would expect faster convergence.This is equivalent to rewriting the Jacobi equation (3.14) asQQS#TXQS#T Ã" & e T $D LuBUXQ [v(3.15)This Gauss–Seidel method has better convergence than the Jacobi method, but only marginally so.As before we consider the example of the solution of the 1D Poisson’s equation. As C programs thebasic structure might be something likeint i;const int N = ??; //incomplete codedouble x[N], b[N];while ( ... ) //incomplete code{for ( i = 0; i < N; i++ )x[i] = (b[i] - x[i-1] - x[i+1]) * 0.5;}Matrix Algebra31Note that only one array is now required to store X, whereas the Jacobi method needed 2.The time required for each iteration is proportional to for each column of B, assuming A is sparse.The number of iterations required depends on the details of the problem, on the quality of the initial guessfor X, and on the accuracy of the required solution.3.6 Matrix Eigenvalue ProblemsNo attempt will be made here to describe the detailed algorithms used for matrix diagonalisation.

A fulldiscussion can be found in the book by Wilkinson and Reinsch (Wilkinson 1964). The routines found inthe common libraries, such as NAG (Numerical Algorithms Group n.d.) or LaPack (Lapack NumericalLibrary n.d.) are almost all based on the algorithms in this book.3.6.1 Schrödinger’s equation$ Ý 0 ½I" & º [(3.16)0The Laplacian, Ý , can be represented in discrete form as in the case of Laplace’s or Poisson’s (section 3.4)€€ € €€equations. For example, in 1D (3.16) becomes$ 3 egT " f O Ky½ & 3 $ 3 S#T º 3(3.17)which can in turn be written in terms of a tridiagonal matrix H as3 º 3 [(3.18)HIn dimensionless form the time–independent Schrödinger equation can be written as4333rAn alternative and more common procedure is to represent the eigenfunction in terms of a linear combination of basis functions so that we have5 5365rr(3.19)" &i … œ ¶ " &\[5The basis functions are usually chosen for convenience and as some approximate analytical solution of theproblem.

Thus in chemistry it is common to choose theto be known atomic orbitals. In solid statephysics often plane waves are chosen.Inserting (3.19) into (3.16) gives55 5555rrr(3.20)¶… œ Ç $ Ý 0 ½I" & È ¶ " &;hº … œ ¶ " &\[Multiplying this by one of the ¶ ’s and integrating gives… 5 > r¶ 7! " r& Ç $ Ý 0 ½I" r& È ¶ 5 " r& œ 5 hº … 5 > r¶ 7! " r& ¶ 5 " r& œ 5 [We now define 2 matrices50>"&$I½"&"&ÈH 89 7r ¶ 7! r Ç Ýr ¶ r55>ß"&"&7S8r ¶ 7! r ¶ r(3.21)5…5 5…5 ߺ7 œœ º SaHa(3.22)so that the whole problem can be written concisely5 5 as59 7(3.23)(3.24)(3.25)55which has the form of the generalised eigenvalue problem (section 3.6.4).

Often the ’s are chosen to be7 and the matrix S is eliminated from the problem.orthogonal so that 7ß O¶Matrix Algebra323.6.2 General PrinciplesThe usual form of the eigenvalue problem is written x(3.26)where A is a square matrix x is an eigenvector and is an eigenvalue. Sometimes the eigenvalue and¥ matrix usually has distincteigenvector are called latent root and latent vector respectively.

An Â:eigenvalue/eigenvector pairs . The full solution of the eigenvalue problem can then be written in the form U; AU ;(3.27)U < A U<(3.28)where is a diagonal matrix of eigenvalues and U ; (U< ) are matrices whose columns (rows) are theeT [correspondingeigenvectors. U < and U ; are the left and right handed eigenvectors respectively, and U <Ax2U;3For Hermitian matrices U< and U ; are unitary and are therefore Hermitian transposes of each other:U< U; .For Real Symmetric matrices U< and U ; are also real.

Real unitary matrices are sometimes calledorthogonal.3.6.3 Full DiagonalisationRoutines are available to diagonalise real symmetric, Hermitian, tridiagonal and general matrices.In the first 2 cases this is usually a 2 step process in which the matrix is first tridiagonalised (transformedto tridiagonal form) and then passed to a routine for diagonalising a tridiagonal matrix. Routines areavailable which find only the eigenvalues or both eigenvalues and eigenvectors. The former are usuallymuch faster than the latter.Usually the eigenvalues of a Hermitian matrix are returned sorted into ascending order, but this is notalways the case (check the description of the routine). Also the eigenvectors are usually normalised tounity.For non–Hermitian matrices only the right–handed eigenvectors are returned and are not normalised.In fact it is not always clear what normalisation means in the general case.Some older FORTRAN and all C (not C++) and PASCAL routines for complex matrices store the realand imaginary parts as separate arrays.

The eigenvalues and eigenvectors may also be returned in this form.This is due to 2 factsThe original routines from Wilkinson and Reinsch (Wilkinson 1964) were written in ALGOL, whichhad no complex type.Many FORTRAN compilers (even recent ones) handle complex numbers very inefficiently, in thatthe use a function even for complex addition rather than inline code. In C++ it is worthwhile checkingthe complex header file to see how this is implemented.3.6.4 The Generalised Eigenvalue ProblemA common generalisation of the simple eigenvalue problem involves 2 matricesAx Bx[egT(3.29)This can easily be transformed into a simple eigenvalue problem by multiplying both sides by the inverseof either A or B.

This has the disadvantage however that if both matrices are Hermitian B A is not, andthe advantages of the symmetry are lost, together, possibly, with some important physics.2 Foran exception consider the matrix = >@?BA@?BC/D .U E may be singular. In such cases the F vectors do not span an F3 Occasionally,dimensional space (e.g. 2 are parallel)Matrix Algebra33There is actually a more efficient way of handling the transformation.

Using Cholesky factorisation anLU decomposition of a positive definite matrix can be carried out such thatBLL(3.30)which can be interpreted as a sort of square root of B. Using this we can transform the problem into theformp L egT A ° L ± eT t · L x¸ B Ay · L x¸ y[(3.31)(3.32)Most libraries contain routines for solving the generalised eigenvalue problem for Hermitian and RealSymmetric matrices using Cholesky Factorisation followed by a standard routine. Problem 6 contains asimple and informative example.3.6.5 Partial DiagonalisationOften only a fraction of the eigenvalues are required, sometimes only the largest or smallest.

Generally ifmore than 10% are required there is nothing to be gained by not calculating all of them, as the algorithms forpartial diagonalisation are much less efficient per eigenvalue than those for full diagonalisation. Routinesare available to calculate thelargest or smallest eigenvalues and also all eigenvalues in a particularrange.Sturm Sequence— — The Sturm sequence is a very nice algorithm found in most libraries. It finds all the eigenvalues in agiven range, minmax , and the corresponding eigenvectors. It is also able to find the number ofsuch eigenvalues very quickly and will return a message if insufficient storage has been allocated for theeigenvectors.

It does require a tridiagonalisation beforehand and is often combined with the Lanczos (seesection 3.6.6) algorithm, to deal with sparse matrices.3.6.6 Sparse Matrices and the Lanczos AlgorithmNone of the above diagonalisation procedures make any use of sparsity. A very useful algorithm fortridiagonalising a sparse matrix is the Lanczos algorithm.This algorithm was developed into a very useful tool by the Cambridge group of Volker Heine (including (e.g.) Roger Haydock, Mike Kelly, John Pendry) in the late 60’s and early 70’s.A suite of programs based on the Lanczos algorithm can be obtained by anonymous ftp from theHENSA4 archive.3.7 Problems€ differenced€ in€ 2 dimensions€ as€Ý 0 ¶ ¯ $ Ü 7€ | is  usually€΀ ¶ egT\á 5¶ S#T^á 5¶ á eT ¶ á S#T $ ׶ ß € €ß € ¯$ O € KÑ0 7 | and Ò : â : : f : [[ [ : $ .whereÜQ ù Qj‰Q HG ù where j ÎS€ ¶The difference equation can be written as a matrix equation äG S ßand.

. Assume the boundaries are at zero potential, i.e. ¶ iffWriteA for Ò or â ¯down$ ortheÒ matrix .or â1. Poisson’s equation4 ftp://unix.hensa.ac.uk/pub/netlib/lanczos/Oscillations of a Crane€ T„ƒ¶‚€ ƒ  € ƒ¶ ‚0 ¶ ‚ ҁ€ Qƒ¶‚΀34ßPTRT and ß ¶ after 8 iterations.2. The equation in question 1 can be solved by the Gauss–Seidel method. Ifotherwise, findandiffor all and , whereis the value ofâ3. Work out the factors, L and U, in the LU decomposition of the matrix,L MNJIKLAf [ ×Hence,Solve the simultaneous equations,AxeT " & b:for a variety of right hand sides, b.Evaluate det AegT "Find A &å4.

Show that the Jacobi method (3.14) is stable as long as, , the eigenvalue of largest modulus ofDL U is less than unity.5. Find both eigenvalues and the corresponding left and right handed eigenvectors of the matrix| $ [Œ ΀ € by an equation of motion in the form6. The vibrational modes of a certain molecule are  described/ * 0 K P O K : KÒwhere / and are the mass and the displacement respectively of the th atom and the real symmetric matrix K describes the interactions between the atoms.Show that this problem can be represented inthe form of a generalised eigenvalue problem: Axö Bx in which the matrix B is positivedefinite.T„’KÃ" / & 0YK , show how to transform this into a simple eigenBBy considering the transformationvalue problem in which the matrix is real symmetric.7.

Write down the types of matrices which occur in the following problems:(a) A simple discretisation (as in question1) of Poisson’s equation in 1 dimension.(b) The same but in more than 1D.(c) A simple discretisation (as above) of Schrödinger’s equation in 3 dimensions.(d) Schrödinger’s equation for a molecule written in terms of atomic basis functions.(e) Schrödinger’s equation for a crystal at a general k point in the Brillouin zone.3.8 Project — Oscillations of a CraneA crane manufacturer wants to know how his newly designed crane is going to behave in a high wind. Hewants to know the oscillation frequencies of a steel wire with and without masses on the end.

This projectis to write a program which calculates the modes of oscillation and their frequencies.Phonons in a Quasicrystal353.8.1 AnalysisThe difference equations for the oscillating wire can be derived by dividing it into N segments (10 or 20should be sufficient) each of which can be considered rigid. This would be the case if the wire consisted ofa series of rods connected together.Letbe the displacement of the bottom of the th segment. Letbe the angle it makes with thevertical.

Let Q be the tension in the th segment. Assume that the mass of the wire is located at the jointsof the segments,at each joint.whereis the length of each segment, and is the massper unit length of the wire. The equation of motion for the mass at the bottom of the nth segment is,KÑQQ8 Q8/ O‰­Ü ‰O ­/ 0Y K 0 Q Q QS#T µ )E¹ "   QS#T&k$ Q Q µ )E¹ "   Q1& [µ "   Q&iÃ"<KÑQo$}KÑQeT &R7 O‰­ .Assume small oscillations so that )E¹/ 0 K10 Q u Q QS#T K QS#T $Š" Q QS#T Q Q &„K Q Q Q K QeT v [O‰­/ 7f/Ü(3.33)(3.34)The mass associated with the end of the wire will be onlysince there is no contribution from the(N+1)th segment.

Consequently, the equation of motion for this point is,RQ Q (3.35)whereis any mass carried by the crane. In addition the displacement of the top of the wire is zero, sothatin the equation for.The modes of oscillation are calculated by seeking solutions of the form. Substituting this into the equations of motion gives,SQQQQ(3.36)R(3.37)QK " 0T / & 11 0K 0 L$8 µ )E¹ "   &iL$ "]K $VK ge T & 7 ‰O ­ :K Q "<%&i Q J!g"AÒ * % &$ / * 0 QÃ" QS#T^QS#TU$l" Q QS=TJ&CnQ QQeT &17 O‰­T$w" 0 / & * 0 "A egTŸ$ø R& 7 O‰­ [The specification of the equations is completed by noting that, from the equilibrium conditions,"R 0T / & x : Q Q Q QS#T Tx / :Q (3.38)xwhere is the acceleration due to gravity.The equations can be organised in the form,L$ M* 0 y :(3.39)Ay<"[[[&T: 0 : : , A is a symmetricwhere y is the column vector of displacements,$ 0 tridiagonal matrix andM is a diagonal matrix.