Implement the Jacobi eigenvalue algorithm for matrix diagonalizaion (also called matrix eigenvalue decomposition or eigendecomposition).
Tasks
Implement a function, timesJ, that multiplies the given
matrix A and the Jacobi matrix J(p,q,θ) in O(n) operations,
A←AJ. Something like
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;
}
}
Implement a function, Jtimes, that multiplies the
Jacobi matrix J(p,q,θ) and the given matrix A in O(n) operations,
A←JA. Something like
void timesJ(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;
}
}
Make a function, say jacobi.cyclic, that implements
the Jacobi eigenvalue algorithm for real symmetric matrices using the
cyclic sweeps: eliminates the off-diagonal elements in a predefined
sequence which spans all off-diagonal elements, for example row after row,
and repeats the sweeps until converged. The convergence criterion could
be that the eigenvalues did not change 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=matrix.get(A,p,q);
double app=matrix.get(A,p,p);
double aqq=matrix.get(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);
Jtimes(A,p,q,-theta); // A←J^T*A*J
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,
(3 points) Hydrogen atom: s-wave radial Schrödinger equation
The s-wave radial Schrödinger equation for the Hydrogen atom reads (in units "Bohr radius" and "Hartree"),
The bound s-state wave-function satisfies this equation and the two boundary conditions,
f(r → 0) = r-r², (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),
f(∞)=0 → f(rmax)=0 .
Numerically one can represent this problem on an equidistant N-point grid with step Δr,
where the wave-function is represented by a vector fi = f(ri), and where the second derivative operation in the Hamiltonian can then be approximated using the three point finite difference formula,
Now the Schrödinger equation becomes a matrix eigenvalue equation,
or, in matrix notation,
H f = ε f ,
where the Hamiltonian matrix H=K+V is given as
( -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 )
( -1/r1 0 0 0 ... 0 )
( 0 -1/r2 0 0 ... 0 )
( 0 0 -1/r3 0 ... 0 )
V = ( ... ... ... ... ... ... )
( 0 ... 0 0 -1/rN-2 0 )
( 0 ... 0 0 0 -1/rN-1 )
Notice that the boundary condition f(0)=0,
f(rmax)=0 is in-built in the K-matrix (and the
V-matrix): indeed the (ommited) components u0
and uN+1 have been assumed to be zero. Thus the
zero-boundary conditions are in-built in the formulation of the matrix
problem.
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 eigenvalues and a matrix V whith the corresponding eigenvectors as columns, such that
∑j HijVjk = ek Vik .
Now calculate numerically (using your diagonalization routine) the lowest egenvalues and eigenfunctions of the s-wave hydrogen atom and compare with the exact results.
Build the Hamiltonian matrix
double 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++){
matrix.set(H,i,i,-2);
matrix.set(H,i,i+1,1);
matrix.set(H,i+1,i,1);
}
matrix.set(H,npoints-1,npoints-1,-2);
H*=-0.5/dr/dr;
for(int i=0;i<npoints;i++)H[i,i]+=-1/r[i];
Diagonalize the matrix using your Jacobi routine.
Investigate convergence of your energies with respect to rmax and dr.
Plot the several lowest eigenfunctions and compare with the analytical result.
(1 point) Scaling and optimization