Computational Physics II
WS 2013/14
Ass.-Prof. Peter Puschnig
Institut f¨ur Physik, Fachbereich Theoretische Physik
Karl-Franzens-Universit¨at Graz
Universit¨atsplatz 5, A-8010 Graz
http://physik.uni-graz.at/~pep
Graz, 14. J¨anner 2014
Prepare-se para as provas
Estude fácil! Tem muito documento disponível na Docsity
Ganhe pontos para baixar
Ganhe pontos ajudando outros esrudantes ou compre um plano Premium
Guias e Dicas
Prepare-se para as provas
Estude fácil! Tem muito documento disponível na Docsity
Encontrar documentos
Prepare-se para as provas com trabalhos de outros alunos como você, aqui na Docsity
Encontra documentos específicos para os exames da tua universidade
Encontra documentos específicos para os exames da tua universidade
Videoaulas
Prepare-se com as videoaulas e exercícios resolvidos criados a partir da grade da sua Universidade
Quiz
Responda perguntas de provas passadas e avalie sua preparação.
Ganhe pontos para baixar
Ganhe pontos ajudando outros esrudantes ou compre um plano Premium
Física Computacional
Tipologia: Notas de estudo
2015
1 / 124
Esta página não é visível na pré-visualização
Não perca as partes importantes!
Documentos relacionados
Pré-visualização parcial do texto
Baixe Física Computacional e outras Notas de estudo em PDF para Física, somente na Docsity!
Computational Physics II
WS 2013/
Ass.-Prof. Peter Puschnig
Institut f¨ur Physik, Fachbereich Theoretische Physik Karl-Franzens-Universit¨at Graz Universit¨atsplatz 5, A-8010 Graz [email protected] http://physik.uni-graz.at/~pep
Graz, 14. J¨anner 2014
ii
iv
Contents
- 1 Linear Systems of Equations
- 1.1 Introduction
- Exercise 1: Matrix Multiplications
- 1.2 The LU Decomposition
- 1.2.1 LU decomposition for general matrices
- Exercise 2: Implementation of the LU-decomposition
- 1.2.2 Pivoting
- 1.2.3 Error estimation and condition number
- Exercise 3: LU-decomposition with LAPACK
- 1.2.4 Linear Least Squares Fits
- Project 1: Force Constant Matrix of a Molecule
- 1.2.5 LU decomposition for tridiagonal matrices
- 1.3 Iterative Methods
- 1.3.1 Construction of iterative methods
- 1.3.2 The Jacobi Method
- 1.3.3 Gauss-Seidel and Successive Over-Relaxation Method
- Exercise 4: Matrix inversion for tridiagonal matrices
- Project 2: Finite difference solution of the stationary heat equation
- 1.3.4 Conjugate Gradient Method
- Exercise 5: Conjugate gradient method
- 1.1 Introduction
- 2 Eigenvalues and Eigenvectors of Matrices
- 2.1 Introduction
- 2.2 Subspace Methods
- 2.2.1 Power iteration (Von Mises Method)
- Exercise 6: Von Mises Method
- 2.2.2 Simultaneous vector iterations
- 2.2.3 Lanczos algorithm
- 2.3 Transformation Methods
- 2.3.1 Jacobi-Method
- Exercise 7: Jacobi Method
- 2.3.2 The QR algorithm
- 2.4 Applications in Physics
- 2.4.1 Normal modes of vibration
- 2.4.2 One-dimensional Boundary Value Problems
- Project 3: Eigenvalues of the stationary Schr¨odinger equation
- 2.4.3 Secular equation of Schr¨odinger equation
- Exercise 8: The Band Structure of fcc-Al
- 3 Partial Differential Equations
- 3.1 Classification of PDEs
- 3.1.1 Discriminant of a quadratic form
- 3.1.2 Boundary vs. initial value problem
- 3.2 Boundary Value Problems
- 3.2.1 Finite differencing in ≥ 2 dimensions
- Exercise 9: Solution of Poisson’s Equation in 2D
- 3.3 Initial Value Problems
- 3.3.1 von Neumann stability analysis
- 3.3.2 Heat equation
- 3.3.3 Time dependent Schr¨odinger equation
- Exercise 10: Time-dependent Schr¨odinger equation in 1D
- 3.3.4 Initial value problems for > 1 space dimensions
- Project 4: Time-dependent Schr¨odinger equation in 2D
- 3.3.5 Consistency, Order, Stability, and Convergence
- Exercise 11: Finite difference schemes for the wave equation
- 3.4 Beyond Finite Differencing
- 3.4.1 Finite Elements Method
- 3.4.2 Finite Volume Method
- 3.1 Classification of PDEs
- A Solutions to Exercises
- Solution to Exercise 1: Matrix Multiplications
- Solution to Exercise 2: LU-decomposition
- Solution to Exercise 3: LU-decomposition with LAPACK
- Solution to Exercise 4: Matrix inversion for tridiagonal matrices
- Solution to Exercise 5: Conjugate gradient method
- Solution to Exercise 6: Von Mises Method
- Solution to Exercise 7: Jacobi Method
- Solution to Exercise 8: The Band Structure of fcc-Al
- Solution to Exercise 9: Solution of Poisson’s Equation in 2D
- Solution to Exercise 10: Time-dependent Schr¨odinger equation in 1D
- Solution to Exercise 11: Finite difference scheme for the wave equation
- B Solutions to Projects
- Solution to Project 1: Force Constant Matrix of a Molecule
- Solution to Project 2: Finite difference solution of the stationary heat equation
- Solution to Project 3: Eigenvalues of the stationary Schr¨odinger equation
- Solution to Project 4: Time-dependent Schr¨odinger equation in 2D
- Bibliography
viii
minant of the coefficient matrix A does not vanish
det A 6 = 0,
i.e., the problem is not singular, the solution can be obtained by applying the well-known Cramer’s Rule. However, this method turns out to be rather impracticable for implementing it in a computer code for n & 4. Moreover, its is computationally inefficient for large matrices. In this chapter, we will therefore discuss various methods which can be implemented in an efficient way and also work in a satisfactory manner for very large problem sizes n. Generally we can distinguish between direct and iterative methods. Direct methods lead in principle to the exact solution, while iterative methods only lead to an approximate solution. However, due to inevitable rounding errors in the numerical treatment on a computer, the usage of an iterative method may prove superior over a direct method particular for very large or ill-conditioned coefficient matrices. In this lecture we will learn about one direct method, the so-called LU - decomposition according to Doolittle and Crout (see Section 1.2) which is the most widely used direct method. In Sec. 1.3 we will outline various iterative procedures including the Gauss- Seidel method (GS), the successive overrelaxation (SOR) method, and the conjugate gradient (CG) approach. The presentation of the above two classes of methods within this lecture notes will be kept to a mini- mum. For further reading, the standard text book ’Numerical Recipes’ by Press et al. is recommended [1]. A detailed description of various numerical methods can also be found in the book by T¨ornig and Spellucci [2]. A description of the LU -decomposition and the Gauss-Seidel method can also be found in the recent book by Stickler and Schachinger [3], and also in the lecture notes of Sormann [4]. Before we start with the description of the LU -decomposition, it will be a good exercise to look at how a matrix-matrix product is implemented Exercise 1. Throughout this lecture notes, Fortran^1 will be used to present some coding examples, though other programming languages such as C, C++ or Python may be used as well. Thus the goal of the first exercise will be to compute the matrix product C = A · B which may be implemented by a simple triple loop
1 do i = 1 , n 2 do j = 1 ,m 3 do k = 1 , l 4 C( i , j ) = C( i , j ) + A( i , k ) ∗B( k , j ) 5 end do 6 end do 7 end do
Apart from this easy-to-implement solution, we will also learn to use highly-optimized numerical libraries. For linear algebra operations, these are the so-called BLAS (Basic Linear Algebra Subpro- (^1) An online tutorial for modular programming with Fortran90 can for instance be found here [5]
2
grams) [6] and the LAPACK (Linear Algebra PACKage) [7]. When using Intel compilers (ifort for Fortran or ifc for C), a particularly simple and efficient way to use the suite of BLAS and LAPACK routines is to make use of Intel’s Math-Kernel-Library [8].
1.2 The LU Decomposition
The direct solution method called LU -decomposition, first discussed by Doolittle and Crout, is based on the Gauss’s elimination method [1, 3, 4]. This in turn is possible due to the following property of a system of linear equations: The solution of a linear system of equation is not altered if a linear combination of equations is added to another equation. Now, Doolittle and Crout have shown that any matrix A can be decomposed into the product of a lower and upper triangular matrix L and U , respectively: A = L · U , (1.4)
where
L =
l 21 1 0 · · · 0 l 31 l 32 1 · · · 0 ... ... ln 1 ln 2 ln 3 · · · 1
and U =
u 11 u 12 u 13 · · · u 1 n 0 u 22 u 23 · · · u 2 n 0 0 u 33 · · · u 3 n ... ... 0 0 0 · · · unn
Once such a decomposition has been found, it can be used to obtain the solution vector x by simple forward- and backward substitution. We can write
A · x = L · U · x = L · y = b,
where we have introduced the new vector y as U · x = y, and see that due the lower triangle form of L, the auxiliary vector y is obtained from forward substitution, that is
yi = bi −
∑^ i−^1 j=
lij yj , i = 1, 2 , · · · , n. (1.6)
Using this result for y, the solution vector x can be obtained via backward substitution, i.e. starting with the index i = n, due to the upper triangle form of the matrix U :
xi = (^) u^1 ii
[
yi −
∑^ n
j=i+
uij xj
]
, i = n, n − 1 , · · · , 1. (1.7)
Of course now the crucial point is how the decomposition of the matrix A into the triangular matrices L and U can be accomplished. This will be shown in the next two subsections, first for a general matrix, and second for a so-called tridiagonal matrix, a form which is often encountered when solving differential equations by the finite difference approach.
5
1.2.1 LU decomposition for general matrices
Here, we simply state the surprisingly simple formulae by Doolittle and Crout for achieving the LU - decomposition of a general matrix A. More details and a brief derivation of the equations below can be found in Ref. [1] and an illustration of the algorithm can be found here: LU.nb.
uij = aij −
∑^ i−^1
k=
likukj , i = 1, · · · , j − 1 (1.8)
γij = aij −
∑^ j−^1 k=
likukj , i = j, · · · , n (1.9)
ujj = γjj (1.10) lij =
γij γjj^ ,^ i^ =^ j^ + 1,^ · · ·^ , n^ (1.11)
Note that the evaluation of the above equations proceeds column-wise. Also note that in the code example below, we have used the fact that both matrices L and U can be stored in only one two- dimensional array LU due to the special shape of the matrices. 1 do j = 1 , n! l o o p o v e r columns 2 do i = 1 , j − 1! Eq. ( 1. 8 ) 3 s = A( i , j ) 4 do k = 1 , i − 1 5 s = s − LU( i , k ) ∗LU( k , j ) 6 end do 7 LU( i , j ) = s 8 end do 9 do i = j , n! Eq. ( 1. 9 ) 10 s = a ( i , j ) 11 do k = 1 , j − 1 12 s = s − LU( i , k ) ∗LU( k , j ) 13 end do 14 g ( i , j ) = s 15 end do 16 LU( j , j ) = g ( j , j )! Eq. ( 1. 1 0 ) 17 i f ( g ( j , j ). eq. 0. 0 d0 ) then 18 g ( j , j ) = 1. 0 d− 30 19 end i f 20 i f ( j. l t. n ) then! Eq. ( 1. 1 1 ) 21 do i = j +1,n 22 LU( i , j ) = g ( i , j ) /g ( j , j ) 23 end do 24 end i f 25 end do! end l o o p o v e r columns
6
1.2.2 Pivoting
As has been demonstrated in the previous exercise, round-off errors in the simple procedure for LU - decomposition may severely depend on the order of the equation rows. Thus, without a proper ordering or permutations in the matrix, the factorization may fail. For instance, if a 11 happens to be zero, then the factorization demands a 11 = l 11 u 11. But since a 11 = 0, then at least one of l 11 and u 11 has to be zero, which would imply either being L or U singular. This, however, is impossible if A is nonsingular. The problem arises only from the specific procedure. It can be removed by simply reordering the rows of A so that the first element of the permuted matrix is nonzero. This is called pivoting. Using such a LU factorization with partial pivoting, that is with row permutations, analogous problems in subsequent factorization steps can be removed as well.
When analyzing the origin of possible round-off errors (or indeed the complete failure of the simple algorithm without pivoting), it turns out that a simple recipe to overcome these problems is to require that the absolute values of lij are kept as small as possible. When looking at Eq. 1.11, this in turn means that the γii values must be large in absolute value. According to Eq. 1.9, this can be achieved by searching for the row in which the largest |γij | appears and swapping this row with the j−row. Only then, Eqs. 1.10 and 1.11 are evaluated. It is clear that such a strategy, partial pivoting, also avoids the above mentioned problem that a diagonal element γjj may turn zero.
If it turns out that all γij for i = j, j + 1, · · · , n are zero simultaneously, then the coefficient matrix Aij is indeed singular and the system of equations can not be solved with the present method.
1.2.3 Error estimation and condition number
While pivoting greatly improves the numerical stability of the LU -factorization, matrices may be ill- conditioned, meaning that the solution vector is very sensitive to numerical errors in the matrix and vectors elements of A and b. Knowing this property about a given matrix is of course very important. In the subsection, we outline the general idea of a condition number which is a measure for how much the output value of a function (solution vector x) can change for a small change in the input argument (inhomogeneous vector b).
Quite generally, the condition number is a means of estimating the accuracy of a result in a given calculation, in particular, this number tells us how accurate we can expect the solution vector x when solving the system of linear equations A · x = b. Assume that there is an error in representing the vector b, call it δb and otherwise the solution is given to absolute accuracy. That is we have to solve A · (x + δx) = b + δb for which we formally get the solution
x + δx = A−^1 · (b + δb) =⇒ δx = A−^1 · δb. (1.12)
8
Using the operator norm of the inverse matrix ‖A−^1 ‖ as defined below, we can estimate the error δx in the solution vector. ‖δx‖ ≤ ‖A−^1 ‖ · ‖δb‖ (1.13)
In order words, the operator norm for the inverse matrix ‖A−^1 ‖ is the so-called condition number for the absolute error in the solution vector. If we want to have a relative error estimate, that is ‖ ‖δxx‖‖ , we can start with the ratio of the relative errors of the solution and the inhomogeneous vectors
‖δx‖ ‖x‖ ‖δb‖ ‖b‖
= ‖ ‖δxx‖‖ · ‖‖δbb‖‖ ≤ ‖A
− (^1) ‖ · ‖b‖ ‖x‖ ≤ ‖A
− 1 ‖ · ‖A‖. (1.14)
Here, we have used the fact that for consistent operator norms, the inequalities ‖δx‖ ≤ ‖A−^1 ‖ · ‖δb‖ and ‖b‖ ≤ ‖A‖ · ‖x‖ hold. Thus, we can finally write for the relative error in the solution vector
‖δx‖ ‖x‖
≤ ‖A−^1 ‖ · ‖A‖ · ‖δb‖ ‖b‖
= cond(A)‖δb‖ ‖b‖
where we have defined the condition number of the matrix A as
cond(A) = ‖A−^1 ‖ · ‖A‖. (1.16)
For calculating the matrix norm, the Euclidean norm is rather difficult to calculate since it would require the full eigenvalue spectrum. Instead, we can use as operator norm either the maximum absolute column sum (1-norm) or the maximum absolute row sum (∞-norm) of the matrix which are easy to compute [9]
‖A‖ 1 = max 1 ≤j≤n
{ (^) n ∑ i=
|aij |
, ‖A‖∞ = max 1 ≤i≤n
{ (^) n ∑ j=
|aij |
1.2.4 Linear Least Squares Fits
In numerous physics applications a set of corresponding data points (xk; yk) are measured or calculated and a set of certain parameters {αj } characterizing a function f (xk; {αj }) are to be determined in such a way that χ^2 =
k
wk [yk − f (xk; {αj })]^2 −→ min (1.18)
This is referred to as a least squares fit problem where wk are weights indicating the relevance of a certain data point (xk; yk) [3]. We distinguish between two different cases: (i) the function f (xk; {αj }) is a linear function of the parameters {αj }, and, (ii), the function f (xk; {αj }) is nonlinear in its parameters {αj }. Here, we restrict ourselves to the linear fit problem since in this case the requirement of the minimization of the least squares sum χ^2 leads to a system of linear equations which can be solved by the method presented in th previous sections.
In the linear case, we can write the test function f (xk; {αj }) as
f (x; {αj }) =
j
αj ϕj (x), (1.19)
where ϕj (x) are linearly independent basis functions. Insertion into Eq. 1.18 and requiring that ∂χ ∂α^2 j = 0 leads to a system of linear equations for the unknown fit parameters αj
M · α = β. (1.20)
Here, the coefficient matrix M and the inhomogeneous vector β are given by the following expression [3]
Mij =
∑^ N
k=
wkϕi(xk)ϕj (xk), (1.21)
βi =
∑^ N
k=
wkykϕi(xk). (1.22)
Note that k runs over all N data points while the size of the matrix Mij is given by the number of fit parameters n.
Project 1
Force Constant Matrix of a Molecule
In this project we will apply our knowledge about how to numerically solve a linear system of equations to a specific linear least squares problem: We want to determine the force constant matrix for atomic displacements in molecule with n vibrational degrees of freedom. The total energy E as a function of the atomic displacement vector u = (u 1 , u 2 , · · · , un) can be written as a Taylor series
E(u) = E( 0 ) +
∑^ n
i=
∂E
∂ui
u= 0
· ui +^12
∑^ n
i=
∑^ n
j=
∂^2 E
∂ui∂uj
u= 0
· uiuj + · · · (1.23)
Taking into account that at the equilibrium geometry u = 0 the forces vanish, i.e., the linear term in the Taylor expansion is zero, and the i−th force component is given by
Fi = −
∂E
∂ui^ =^ −
∑^ n j=
φij uj with φij =
∂^2 E
∂ui∂uj
u= 0
The matrix elements φij of the force constant matrix can be determined by a fit to N data sets (u(k)), F (k)). For instance, consider an H 2 O molecule assuming all atoms in the xy-plane. There are n = 6 degrees of freedom and correspondingly for each data set (k) there will be six coordinates and six forces. The structure of the data file to be read in is illustrated for the first two data sets^2 1 C o n f i g u r a t i o n # 000 2 x y Fx Fy 3 0. 5 9 5 6 5 0 0 0. 7 6 4 9 2 6 1 0. 0 0 0 4 1 7 0 −0. 4 0. 5 9 5 6 5 0 0 −0.7649261 0. 0 0 0 4 1 7 0 0. 0 0 3 0 6 1 0 5 0. 0 0 0 0 0 0 0 0. 0 0 0 0 0 0 0 −0.0008340 0. 0 0 0 0 0 0 0 6 C o n f i g u r a t i o n # 001 7 x y Fx Fy 8 0. 5 9 2 8 8 6 5 0. 7 6 1 0 1 4 4 0. 1 6 0 7 8 7 0 0. 1 8 9 7 4 3 0 9 0. 5 9 9 1 0 0 2 −0.7695293 −0.1226700 0. 1 7 4 9 8 5 0 10 0. 0 0 1 4 4 0 2 −0.0008279 −0.0381170 −0. 11...
Now, the linear fit problem can be defined by demanding
χ^2 =
∑^ N
k=
F (k)^ − F (u(k); φij )
−→ min ⇐⇒ ∂χ
2 ∂φij
(^2) Coordinates are given in ˚A, and forces are in eV/˚A. The data has been obtained from ab-initio quantum-mechanical calculations within the framework of density functional theory.
12