Implement the Jacobi eigenvalue algorithm for matrix diagonalizaion (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 eigenvectors (let's call it "V"). There are several options for the interface to your routines:
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);
}
}
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 */
}
}
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
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;
}
}
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;
}
}
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,
(3 points) Hydrogen atom, s-wave radial Schrödinger equation on a grid
The s-wave (l=0) radial Schrödinger equation for the Hydrogen atom reads (in units "Bohr radius" for length and "Hartree" for energy),
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),
f(∞)=0 ⇒ f(rmax)=0 .
Numerically one can represent the above problem on an equidistant grid with step Δr,
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,
With the grid representation 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-1 0 )
( 0 ... 0 0 0 -1/rN )
Notice that the boundary condition f(0)=0,
f(rmax)=0 is inbuilt in the K-matrix and the
V-matrix: indeed the (omitted) components f0
and fN+1 are set to zero. Thus the (zero-) boundary
conditions are inbuilt 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 the eigenvalues and a matrix V whith the corresponding eigenvectors as columns, such that
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;
H[i,i+1]= 1;
H[i+1,i]= 1;
}
H[npoints-1,npoints-1]=-2;
H.scale(-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.
Investigate convergence of your energies with respect to rmax and Δr :
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) .
(1 point) Scaling and optimization