计算物理ComputationalPhysics计算物理 (6).pdf

《计算物理ComputationalPhysics计算物理 (6).pdf》由会员分享，可在线阅读，更多相关《计算物理ComputationalPhysics计算物理 (6).pdf（40页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Computational physiCsNumerical methods for matrices Matrices in physics Basic matrix operations Linear equation systems Eigenvalue problemsMany problems in physics can be formulated in a matrix form Q:Where is matrix?A:Matrices are everywhere in physics!An eigenvalue problem given in the form of a p

2、artial differential equation can be rewritten as a matrix problem.A boundary-value problem after discretization is essentially a linear algebra problem.The vibrational spectrum of a molecule with n vibrational degrees of freedom can be also sloved using matrix methods.Vibration of a Moleculenjijiijn

3、qqAqqqU1,2121),(Elastic energyKinetic energynjijiijnqqMqqqT1,2121),(Lagrange equationL=T-U0iiqLdtdqLnnnnnnnnnnxxMMMMxxAAAA1111111111The roots of this secular equation,k with k=1,2,.,n,give all the possible vibrational angular frequencies wk2=k of the molecule.2 w MxAxThe determinant0MALet:0)(1njjijj

4、ijqMqA tijjexqw0)(12jnjijijxMAwKirchhoff Equations for Circuits The unbalanced Wheatstone bridge,0)(,0)(,v3322133211032211irrririririrrririri riraaaaxxsvR iaaaaxxsrrrrrrrrrrrrr323121RvR1-ii1i2i3i1-i2i2-i3i1-i3Ohms lawMatrix Quantum Mechanics In fact,any operators can be written as matrices.For examp

5、le:The spin-1/2 electron.Bases:spin up:|;spin down:|Choosing a set of bases:orthogonal,complete,unitary.Then the Hamiltonian can be written as a matrix with the elements as:jiijHThe Schrodinger equationkkk 1001 00 0110321zyxiiBasic matrix operations An n x m matrix A is defined through its elements

6、Aij with the row index i=1,2,.,n and the column index j=1,2,.,m.It is called a square matrix if n=m.A variable array x with elements x1,x2,.,xn arranged into a column is viewed as an n 1 matrix,or an n-element column matrix.A typical set of linear algebraic equations is given byfor i=1,2,.,n,where x

7、j are the unknowns to be solved,Aij are the given coefficients,and bi are the given constants.njijijbxA1bAxThe standard matrix multiplication:The summation over k requires the number of columns of the first matrix to be the same as the number of rows of the second matrix.Otherwise,the product does n

8、ot exist.The inverse of a square matrix A(written as A-1)is defined bywhere I is a unit matrix with the elements Iij=dij.ABC kkjikijBACIAAAA-1-1 for any j=1,2,.,n,where|Rij|is the determinant of the residual matrix Rij of A with its ith row and jth column removed.Cij=(-1)i+j|Rij|is called a cofactor

9、 of Aij.The determinant of an n x n matrix A is defined as:ijijnijiRA1)1(ADeterminantInverse of A In principle,the inverse of A can be obtained throughIf a matrix has an inverse or nonzero determinant,it is called a nonsingular matrix.Otherwise,it is a singular matrix.AAjiijC1-The trace of a matrix

10、A is the sum of all its diagonal elements,written asThe transpose of a matrix A(written as AT)has elements with the row and column indices of A interchanged,that is,We call A an orthogonal matrix if AT=A-1.The complex conjugate of AT is called the Hermitian operation of A(written as A)with Aij=A*ji.

11、We call A a Hermitian matrix if A=A and a unitary matrix if A=A-1.niiiAA1Tr jiijAA Tcode example To calculate the determinant,trace,transpose,and inverse of a square matrix.5.1.Matrix.cppTriangular matrixA matrix is called an upper-triangular(lower-triangular)matrix if the elements below(above)the d

12、iagonal are all zero.Ax=b.|A|!=0 and b!=0 -a unique solution.Method:Gaussian elimination-to transform the original matrix A to a triangular matrix The original matrix A=A(0)A(j):the resultant matrix after j matrix operations.Similar notation is used for the transformed b:b(j)Linear equationGaussian

13、elimination for linear equation systems CoefficientmatrixA(0)&b(0)intermediatematrixA(j)&b(j)Gaussianeliminationtriangularmatrixroots of equationsbackward substitutionsThe inverse and the determinant of a matrix can also be obtained in such a manner.）（）（）（）（）（）（）（）（0000000012222111211nnnnnAAAAAAAAmu

14、ltiply the first equation by-Ai1(0)/A11(0)add it to the ith equation for i 1）（）（）（）（）（）（）（）（）（）（）（）（0000000000001222211111111121nnnniniiAAAAAAAAAAAA)1()1()1()1()1()1(0022211211nnnnAAAAAAmultiply the second equation by-Ai2(1)/A22(1)1()1()1()1()1()1()1()1(002222211211nnininAAAAAAAAadd it to the ith eq

15、uation for i 2)2()2()2()2()2()2(00022211211nnnnAAAAAAcontinued with the third,fourth,.,and(n-1)th equations)1()1()1()1()1()1(000000022211211nnnnnnnnnnAAAAAAThen the coefficient matrix becomes an upper-triangular matrix A(n-1).A linear equation set with an upper-triangularcoefficient matrix can easil

16、y be solved with backward substitutions.Because all the diagonal elements are used in the denominators,the scheme would fail if any of them happened to be zero or a very small quantity.This problem can be circumvented in most cases by interchanging the rows to have the elements used for divisions be

17、ing the ones with largest magnitudes possible.This is the so-called pivoting procedure.This procedure will not change the solutions of the linear equation set.The partial-pivoting scheme,which searches for the pivoting element only from the remaining elements of the given column.We first search for

18、the element with the largest magnitude from|Ai1(0)|for i=1,2,.,n.Assuming that the element obtained is Ak11(0),we then interchange the first row and the k1th row and eliminate the first element of each row except the first row.Similarly,we can search for the second pivoting element with the largest

19、magnitude from|Ai2(1)|for i=2,3,.,n.Solution of a linear equation set After the Gaussian elimination,the solution of a linear equation set is then obtained through backward substitutions withfor i=n-1,n-2,.,1,starting with)1()1()1()1()1()1(000000022211211nnnnnnnnnnAAAAAADeterminant of a matrix The d

20、eterminant of the original matrix can easily be obtained after we have transformed it into a triangular matrix through Gaussian elimination.The partial-pivoting Gaussian elimination does not change the value of the determinant only its sign,which can be fixed with the knowledge of the order of pivot

21、ing elements.For a triangular matrix,the determinant is given by the product of all its diagonal elements.Inverse of a matrix The inverse of a matrix A can be obtained by solving the linear equation set.AA-1=I.I:a unit matrix Let A-1ij=xij,for i,j=1,2,.,n,then xij is the solution of the equation Axj

22、=bj,with bij=ij.The solution of above equation corresponds to each column of the unit matrix forms the corresponding column of A1.Code example5.2.GaussianElimination.cppHomework Solve the unbalanced Wheatstone bridge.i1i2i3i1-i2i2-i3i1-i3Requirements:Input the values of resistors.Output the effectiv

23、e resistane of the circuit.Matrix Eigenvaluewhere x and are an eigenvector and its corresponding eigenvalue of the matrix,respectively.Determined from the secular equation:|A-I|=0An nxn matrix has a total of n eigenvalues.xxA The eigenvalue problem is quite general in physics,which can come from man

24、y different problems.For example,the Lagrange equation for the vibrational modes of a large molecule,or the quantum mechanics Hamiltonian.In many problems in physics and related fields,the matrix in question is Hermitian with A=A.Eigenvalues of a Hermitian matrixThe simplicity of the Hermitian eigen

25、value problem is due to three important properties:1.the eigenvalues are all real;2.the eigenvectors can be orthonormal;3.can be transformed into a diagonal matrix with the same set of eigenvalues under a similarity transformation of a unitary matrix.Similarity transformation:B=P-1AP.Unitary matrix:

26、PP=I,P=P-1;means the conjugate transpose The eigenvalue problem of an n x n complex Hermitian matrix is equivalent to that of a 2n x 2n real symmetric matrix.Separate its real part from its imaginary part with A=B+iC.B is a real symmetric matrix:Bij=Bji C is a real skew symmetric matrix:Cij=-Cji.Dec

27、ompose the eigenvector z in a similar fashion:z=x+iy.The original eigenvalue problem becomes(B+iC)(x+iy)=(x+iy),which is equivalent toit becomes a real symmetric eigenvalue problem with the same set of eigenvalues that have an overall double degeneracy.Therefore we need to solve only the real symmet

28、ric eigenvalue problem if the matrix is Hermitian.Diagonalization For a diagonal matrix,the eigenvalues are the diagonal elements.Diagonalize a matrix -DiagonalizationTo obtain the eigenvalues of a Hermitian matrix=to transform it to a diagonal matrix by similarity transformationsLinear Algebra PACK

29、age LAPACK(Linear Algebra Package)is a software library for numerical linear algebra,provided by Univ.of Tennessee;Univ.of California,Berkeley;Univ.of Colorado,Denver;and NAG Ltd.It provides routines for solving systems of linear equations and linear least squares,eigenvalue problems,and singular va

30、lue decomposition.The routines handle both real and complex matrices in both single and double precision.LAPACK was originally written in FORTRAN 77,but moved to Fortran 90 in version 3.2(2008).The latest version is 3.9.0(Nov.2019).LAPACK is licensed under a three-clause BSD style license,a permissi

31、ve free software license with few restrictions.Website lib.org/lapackMore implements Lapack is also included(and may be optimized for specialized platforms)in MKLIntel,ACMLAMD,Matlab,Mathematica,.It can be included in commercial software packages(and has been).Netlib only ask that proper credit be g

32、iven to the authors.Code example 5.3.Diagonalization.cpp 9.6.Diagonalization.cpphttp:/lib.org/lapack/lapacke.htmlHydrogen molecule H=K+UR=0.53 AHydrogen molecule H=K+U=K1+U1+K2+U2+K12=H0+K12+r=0.74 AHydrogen molecule+r=0.74 AE=V+tE=V-tantibonding state:(|1-|2)/sqrt(2)Reaction energy:4.73 eVt=2.365 eVbonding state:(|1+|2)/sqrt(2)Hydrogen pyramid 5.3.H4.cppHome work One dimensional hydrogen chain Eigen-energies?Surface(terminal)energy?