McMaster University, Physics 4G03/6G03 October 8 - October 29 2020
Homework 3 :
Eigenvalue Problems, Jacobi and Lanczos methods
You should submit source and makefiles for the exercises marked with PROGRAM.
Files for this assignment are available in /home/4G03/HW3 on phys-ugrad
Due: Thursday October 29 2020, 11:30pm on avenue
☛Reading: In the first weeks we have covered material, parts of which can be found in
Chapter 2, Chapter 3 and Chapter 10 (Monte Carlo) of Pang. This week we will discus diagonalization
for dense (Jacobi) and sparse (Lanczos) matrices and continue with Linear Algebra
and matrices (Pang Chapter 5).
Exercise 1 Quantum States in a Potential Well
You may remember from Physics 3MM3 (or similar Quantum Mechanics course) that it is
possible to solve a differential equation in terms of the Legendre and Laguerre polynomials.
However, in practice it is often much more useful to solve Schr¨odinger’s equation by mapping
it onto a Matrix which we can solve numerically for its eigenvalues and eigenvectors. That’s
the goal of this exercise. You may find that basic QM like the material developed in Chapter
3 in Griffiths can be quite handy, so you may want to review that a bit. We wish to find the
eigenvalues and eigenfunctions of the following 1D Schr¨odinger equation:
H|ψi = E|ψi, H = H0 + V (1). (3)
As you can see this potential diverges at x = 0, π/2 so we may think of the particle as being
confined to this interval. The eigenstates of H0 on the same interval are simply the eigenstates
of the infinite potential well of width L. They are well known:
They also form a complete set of basis states on the Hilbert space defined on the interval
0 < x < L. We can therefore expand the solution we are looking for in the following manner:(7)
If we insert this expression in Eq. (1) we find:
If we take the inner product of Eq. (8) with the bra hψn| we obtain the matrix equation:
X∞
m=1
Hnmcm = Ecn, (9)
where
Hnm = hψn|(H0 + V )|ψmi = δnmE
0
n + hψn|V |ψmi. (10)
Note that this matrix equation is infinite dimensional ! Hence, for a practical application we
will have to make it finite-dimensional by simply truncating the basis set:
|ψi ' X
M
m=1
cm|ψmi. (11)
With M = 5 the finite matrix equation then looks like this
H11 H12 H13 H14 H15
H21 H22 H23 H24 H25
H31 H32 H33 H34 H35
H41 H42 H43 H44 H45
H51 H52 H53 H54 H55. (12)
Solving the Schr¨odinger equation has become equivalent to determining the eigenvalues and
eigenvectors of a matrix.
1.1 For the potential in Eq. (4) write out the 5 × 5 matrix Hnm that corresponds to truncating
the basis set with M = 5. You may find the following integrals useful:
dx = (−1)|m−n|π min(n, m), (14)
with m, n both integers. When using this integrals be careful with overall constants. (No
programming needed).
1.2PROGRAM Write a module that implements Jacobi’s method for diagonalizing. Note that
in the handout Eq. 11.1.11 should read 1/
√
t
2 + 1. We shall use that to study the eigenvalues
of the matrix H. You should write this as a subroutine starting something like this:
2
#include
#include
using namespace std;
typedef vector Row; // One row of the matrix
typedef vector Matrix; // Matrix: a vector of rows
void jacdiag(Matrix & A, vector & d)
.
.
Here A is the matrix to be diagonalized and d is a vector that on output contains the calculated
eigenvalues in any order. Try the simplest possible implementation first, by doing 10 sweeps
through the elements apq above the diagonal and performing a jacobi rotation on each one. Use
the following equations, assuming that |apq| > 2.2 × 10−16:
Then perform the rotation as discussed in class. In a first implementation you can try something
like this in pseudocode: (watch out I haven’t checked this)
AP=A(p,:)
A(p,:)=c*A(p,:)-s*A(q,:)
A(q,:)=s*AP +c*A(q,:)
AP=A(:,p)
A(:,p)=c*A(:,p)-s*A(:,q)
A(:,q)=s*AP +c*A(:,q)
A(p,q)=0.0_dbl
A(q,p)=0.0_dbl
After each sweep, print out the sum of the squares of the remaining elements above the diagonal.
Then try to add some of the refinements described in the handout and write a routine that only
work on the upperdiagonal of A. The matrix A is assumed symmetric so the transformations
you perform should preserve that. A simple way of doing that is to use only elements above
the diagonal (as it is done in the handout from Numerical Recipes by W. Press et al.). Check
your program on simple matrices where you know the result.
For this assignment there is no need to also calculate the eigenvectors. However, if you
feel up to it you can add the calculation of the eigenvectors to your implementation by adding
another argument Matrix & V.
1.3PROGRAM Use your implementation of Jacobi’s routine to determine the eigenvalues of the
matrix equation with M = 10, 20, 30, 40. Determine the three lowest eigenvalues this way for
the different values of M. The eigenvalues should be in units of ~
2/(2m). This method
becomes exact in the limit where M → ∞. How many digits have you reliably determined ?
Exercise 2 S = 1/2 Heisenberg Spin Chain
While Jacobi’s method is very reliable it is not very efficient for very large matrices. In that
case we have to use sparse matrix methods. We shall look at spin waves as the occur in the
Heisenberg spin chain. We will do this by exactly diagonalizing a small model system using
Lanczos method. Hopefully you will be able to study the ground-sate properties of a system
with 20 spins corresponding to a matrix of size 220 × 2
20 or 1048576 × 1048576 !!! In part, the
3
2016 nobel prize in physics was awarded to Duncan Haldane for studies of this system. I’ve
put the nobel committees overview of the scientific background on avenue.
Some materials consists of long 1-dimensional chains of atoms. These chains interact only
weakly and hence we can model the magnetic properties of such a material by a linear chain
of interacting quantum mechanical spins as shown in the figure below. On the i’th site of the
figure 1: The Heisenberg Spin Chain.
The matrices are here written in terms of the z-component of the spin, that’s why σ
z
is diagonal.
The z-component of the spin can be either up, | ↑i, or down, | ↓i. If we have two spins we
need to work in a product basis, | ↑↑i, | ↑↓i, | ↓↑i, | ↓↓i. We see that we have 4 states. If we
have 3 spins interacting we get 8 states and in the general case for L spins interacting we get
N = 2L
states. The interaction between 2 spins is usually written in terms of the Heisenberg
interaction S~i· S~j, with a little bit of work it is not difficult to see that in the above basis we
find that:
If we number the basis states 0, 1, 2, 3 we see that the action of the matrix P01, describing the
interaction between the spins on site 0 and 1, is the following:
nothing more than permuting the coefficient in the vector. It is easy to see that the coefficient
xi
is mapped onto the coefficient j of the resulting vector where j is equal to i with the bits 0,1
interchanged. In general 2S~k · S~l +12I permutes the bits k and l. If we have L interacting spin
their interactions are described by the Hamiltonian H which we can write:
Here it is implied that we are using periodic boundary conditions in the sense that across the
boundary we get S~
L−1 · S~
0.
file: hv.cc
#define BIT_SET(a,b) ((a) |= (1U<<(b)))
#define BIT_CLEAR(a,b) ((a) &= ~(1U<<(b)))
#define BIT_FLIP(a,b) ((a) ^= (1U<<(b)))
#define BIT_CHECK(a,b) ((bool)((a) & (1U<<(b))))
// Set on the condition f else clear
//bool f; // conditional flag
//unsigned int m; // the bit mask
//unsigned int w; // the word to modify: if (f) w |= m; else w &= ~m;
#define COND_BIT_SET(a,b,f) ((a) = ((a) & ~(1U<<(b))) | ((-(unsigned int)f) & (1U<<(b))))
#include
#include
#include
using namespace std;
void hv(vector& y, const vector& x,int L)
{
//for (vector::iterator it = y.begin() ; it != y.end(); ++it)
// *it=0.0;
bool b ;
unsigned int k;
for (unsigned int i=0;i