←to practical programming

Homework EVD (eigenvalue decomposition)

Objective

Implement the Jacobi eigenvalue algorithm for matrix diagonalization (matrix eigenvalue decomposition).

Jacobi eigenvalue algorithm is not the fastest. However, it is probably the most stable and the most accurate. And the easiest to implement.

The EVD decomposition of a matrix is given by the vector of the eigenvalues (let's call it "w", like in numpy) and the matrix of the corresponding eigenvector-columns (let's call it "V"). There are several options for the interface to your routines:

  1. Return "(w,V)" in a tuple (like in "numpy.linalg.eig"), something like 
    public static class jacobi{
public static void timesJ(matrix A, int p, int q, double theta){/*...*/}
public static void Jtimes(matrix A, int p, int q, double theta){/*...*/}
public static (vector,matrix) cyclic(matrix M){
	matrix A=M.copy();
	matrix V=matrix.id(M.size1);
	vector w=new vector(M.size1);
	/* run Jacobi rotations on A and update V */
	/* copy diagonal elements into w */
	return (w,V);
	}
}
  2. Keep "w" and "V" in an instance of the class (let's call it "EVD)", something like 
    public class EVD{
public vector w;
public matrix V;
public static void timesJ(matrix A, int p, int q, double theta){/*...*/}
public static void Jtimes(matrix A, int p, int q, double theta){/*...*/}
public EVD(matrix M){
	matrix A=M.copy();
	V=matrix.id(M.size1);
	w=new vector(M.size1);
	/* run Jacobi rotations on A and update V */
	/* copy diagonal elements into w */
	}
}
  3. Write your eigenvalues and eigenvectors into the caller-provided vector "w" and matrix "V" (like in GSL), 
    public static class jacobi{
public static void timesJ(matrix A, int p, int q, double theta){/*...*/}
public static void Jtimes(matrix A, int p, int q, double theta){/*...*/}
public static void cyclic(matrix A, vector w, matrix V){
	/* run Jacobi rotations on A keeping diagonal elemets in w and updating V */
	}
}

Tasks

  1. (6 points) Jacobi diagonalization with cyclic sweeps

    1. Make a function, timesJ, that multiplies (in-place) the given matrix A with the Jacobi matrix J(p,q,θ) from the right, A←AJ, in O(n) operations. Something like 

      static void timesJ(matrix A, int p, int q, double theta){
	double c=cos(theta),s=sin(theta);
	for(int i=0;i<A.size1;i++){
		double aip=A[i,p],aiq=A[i,q];
		A[i,p]=c*aip-s*aiq;
		A[i,q]=s*aip+c*aiq;
		}
}

    2. Make a function, Jtimes, that multiplies (in-place) the given matrix A with the Jacobi matrix J(p,q,θ) from the left, A←JA, in O(n) operations. Something like 

      static void Jtimes(matrix A, int p, int q, double theta){
	double c=cos(theta),s=sin(theta);
	for(int j=0;j<A.size1;j++){
		double apj=A[p,j],aqj=A[q,j];
		A[p,j]= c*apj+s*aqj;
		A[q,j]=-s*apj+c*aqj;
		}
}

    3. Make a function, say cyclic, that implements the Jacobi eigenvalue algorithm for real symmetric matrices using the cyclic sweeps: eliminate the off-diagonal elements in a predefined sequence which spans all off-diagonal elements, for example row after row, repeating the sweeps until converged. The convergence criterion could be that none of the diagonal elements changed after a whole sweep. Something like 

      bool changed;
do{
	changed=false;
	for(int p=0;p<n-1;p++)
	for(int q=p+1;q<n;q++){
		double apq=A[p,q], app=A[p,p], aqq=A[q,q];
		double theta=0.5*Atan2(2*apq,aqq-app);
		double c=Cos(theta),s=Sin(theta);
		double new_app=c*c*app-2*s*c*apq+s*s*aqq;
		double new_aqq=s*s*app+2*s*c*apq+c*c*aqq;
		if(new_app!=app || new_aqq!=aqq) // do rotation
			{
			changed=true;
			timesJ(A,p,q, theta); // A←A*J 
			Jtimes(A,p,q,-theta); // A←JT*A 
			timesJ(V,p,q, theta); // V←V*J
			}
	}
}while(changed);

      Prove that your implementation works as intended: generate a random symmetric matrix A, apply your routine to perform the eigenvalue-decomposition, A=VDVT (where V is the orthogonal matrix of eigenvectors and D is the diagonal matrix with the corresponding eigenvalues), and check that VTAV==D, VDVT==A, VTV==1, VVT==1,

  2. (3 points) Hydrogen atom, s-wave radial Schrödinger equation on a grid

    Radial Schrödinger equation and boundary conditions

    The s-wave (l=0) radial Schrödinger equation for the Hydrogen atom reads (in units "Bohr radius" for length and "Hartree" for energy),

    -(1/2)f'' -(1/r)f = εf ,

    where f(r) is the reduced radial wave-function, f(r)=r*Ψ(r); ε is the energy; and primes denote the derivative over r.

    The bound s-wave wave-function satisfies the above equation and the two boundary conditions, 

    f(r → 0) = r-r² → 0, (prove this)
    
f(r → ∞) → 0 .

    These two boundary conditions can only be satisfied for certain discrete (negative) values of the energy.

    Since one cannot integrate numerically to ∞ one substitutes ∞ with a reasonably large number, rmax, such that it is much larger than the typical size of the hydrogen atom but still managable for the numerical inregrator (say, rmax = 10 Bohr radii for the ground state). The boundary conditions then become 

    f(0)=0 ,  f(rmax)=0 .
    Representation on a grid

    Numerically one can represent the above problem on an equidistant grid with step Δr,

    ri = i*Δr , i = 1...N , N = rmax/Δr - 1 ,

    where the wave-function is represented by a vector {fi = f(ri)}i=1…N , and where the second derivative in the kinetic energy operator can be approximated using the three point finite difference formula,

    f''i = (fi-1-2fi+fi+1)/Δr² .

    With the grid representation the Schrödinger equation becomes a matrix eigenvalue equation,

    ∑j=1..N-1 Hij fj = εfi , i=1...N ,

    or, in matrix notation, 

    
H f = ε f ,

    where the Hamiltonian matrix H=K+W is given as a sum of the kinetic energy matrix K, 

                ( -2   1   0   0  ...  0  )
            (  1  -2   1   0  ...  0  )
            (  0   1  -2   1  ...  0  )
K = - 1/2Δr² ( ... ... ... ... ... ... )
            (  0  ...  0   1  -2   1  )
            (  0  ...  0   0   1  -2  )
    and the potential energy matrix W, 
        ( -1/r1  0    0    0   ...     0  )
    (   0  -1/r2  0    0   ...     0  )
    (   0    0  -1/r3  0   ...     0  )
W = (  ...  ...  ...  ...  ...    ... )
    (   0   ...   0    0  -1/rN-1   0  )
    (   0   ...   0    0    0   -1/rN )
    Notice that the boundary conditions f(0)=0, f(rmax)=0 are in-built in the K-matrix and the W-matrix: indeed the (omitted) components f0 and fN+1 are set to zero. Thus the matrix formulation naturally satisfies the zero-boundary conditions.

    Solving the matrix eigenvalue equation amounts to matrix diagonalization. If you give the matrix H to your diagonalization routine it will return a vector e with the eigenvalues and a matrix V whith the corresponding eigenvectors as columns, such that

    ∑j HijVjk = Vik ek .
    Numerical calculation

    Now, calculate numerically (using your diagonalization routine) the lowest egenvalues and eigenfunctions of the s-wave states in the hydrogen atom and compare with the exact results.

    Build the Hamiltonian matrix 

    /* your main should get rmax and dr from the command line, like this */
/* mono main.exe -rmax:10 -dr:0.3 */
int npoints = (int)(rmax/dr)-1;
vector r = new vector(npoints);
for(int i=0;i<npoints;i++)r[i]=dr*(i+1);
matrix H = new matrix(npoints,npoints);
for(int i=0;i<npoints-1;i++){
   H[i,i]  =-2*(-0.5/dr/dr);
   H[i,i+1]= 1*(-0.5/dr/dr);
   H[i+1,i]= 1*(-0.5/dr/dr);
  }
H[npoints-1,npoints-1]=-2*(-0.5/dr/dr);
for(int i=0;i<npoints;i++)H[i,i]+=-1/r[i];

    Diagonalize the matrix using your Jacobi routine and obtain the eigenvalues and igenvectors.

    Convergence

    Investigate convergence of your energies with respect to rmax and Δr :

    • Fix rmax to a reasonable value and calculate ε0 for several different values of Δr ; plot the resulting curve.
    • Fix Δr to a reasonable value and calculate ε0 for several different values of rmax ; plot the resulting curve.
    Wave-functions

    Plot several lowest eigen-functions and compare with analytical results.

    The hydrogen s-wave reduced radial eigen-function f(k)(r) with the radial quantum number k is represented by the k'th eigenvector of our Hamiltonian matrix. That is, by the k'th column of the matrix V, 

    f(k)(ri) = Const×Vik ,
    where the columns of the matrix V are normalized, 
    ∑iVik² = 1 .
    The constant is found from the normalization condition for the radial functions, 
    ∫0∞ |f(k)(r)|²dr ≈ ∑i|f(k)(ri)|²Δr = ∑iConst²|Vik|²Δr = Const²Δr = 1 ,
    which gives 
    Const = 1/√(Δr) .

  3. (1 point) Scaling and optimization At your choice, do one (or several) of the following exercises,