Implement the Jacobi 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(gsl_matrix* A, int p, int q, double theta){
double c=cos(theta),s=sin(theta);
for(int i=0;i<A->size1;i++){
double new_aip=c*gsl_matrix_get(A,i,p)-s*gsl_matrix_get(A,i,q);
double new_aiq=s*gsl_matrix_get(A,i,p)+c*gsl_matrix_get(A,i,q);
gsl_matrix_set(A,i,p,new_aip);
gsl_matrix_set(A,i,q,new_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 Jtimes(gsl_matrix* A, int p, int q, double theta){
double c=cos(theta),s=sin(theta);
for(int j=0;j<A->size2;j++){
double new_apj= c*gsl_matrix_get(A,p,j)+s*gsl_matrix_get(A,q,j);
double new_aqj=-s*gsl_matrix_get(A,p,i)+c*gsl_matrix_get(A,q,j);
gsl_matrix_set(A,p,j,new_apj);
gsl_matrix_set(A,q,j,new_aqj);
}
}
Make a function, jacobi_diag, 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
int changed;
do{
changed=0;
for(int p=0;p<n-1;p++)
for(int q=p+1;q<n;q++){
double apq=gsl_matrix_get(A,p,q);
double app=gsl_matrix_get(A,p,p);
double aqq=gsl_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=1;
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!=0);
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,
(3 points) Quantum particle in a box
Consider a quantum particle in a box (or a standing wave in a string). You might also want to read the articles "Quantum particle in a box" and/or "Standing wave in a string".
The wave-function ψ(x) of a particle in a box of length L satisfies the Schrödinger equation
Ĥψ = Eψ
where E is the energy of the particle and where the Hamiltonian Ĥ is given as
Ĥ = - ℏ²/2m d²/dx²
where m is the mass of the particle. The boundary condition for an infinitely deep box is that the wave-function vanishes at the borders,
ψ(0) = ψ(L) = 0 .
Let us introduce dimensionless variables
ξ=x/L , ε = 2mL²E/ℏ² ,
and the new function
u(ξ)=ψ(ξL) .
The Schrödinger equation then reduces to
ĥ u = ε u ,
where
ĥ = - d²/ξ² ,
with the boundary condition
u(0)=u(1)=0 .
Numerically one can represent this problem on an equidistant N-point grid with step s,
ξi = i·s , s=1/N , i=0...N ,
where the wave-function is represented by a vector ui = u(ξi), and where the second derivative operation in the Hamiltonian can then be approximated using the three point finite difference formula,
d²ui/dξ² = (ui-1-2ui+ui+1)/s² .
Now the Schrödinger equation becomes a matrix eigenvalue equation,
∑j=1..N-1 hij uj = εui , i=1...N-1 ,
or, in matrix notation,
h u = ε u ,
where the matrix h is given as
( -2 1 0 0 ... 0 )
( 1 -2 1 0 ... 0 )
( 0 1 -2 1 ... 0 )
h = - 1/s² ( ... ... ... ... ... ... )
( 0 ... 0 1 -2 1 )
( 0 ... 0 0 1 -2 )
Notice that—since the wave-function at the end-points is
zero—the summation goes over 1...N-1 and the vanishing
components u0=0 and uN=0 have been omitted and the boundary
condition is thus 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 egenvalues and eigenfunctions for a quantum particle in a box and compare with the exact results.
Build the Hamiltonian matrix
int n=20;
double s=1.0/(n+1);
gsl_matrix* H = gsl_matrix_alloc(n,n);
for(int i=0;i<n-1;i++){
gsl_matrix_set(H,i,i,-2);
gsl_matrix_set(H,i,i+1,1);
gsl_matrix_set(H,i+1,i,1);
}
gsl_matrix_set(H,n-1,n-1,-2);
gsl_matrix_scale(H,-1/s/s);
Diagonalize the matrix using your Jacobi routine,
gsl_matrix V = gsl_matrix_alloc(n,n); jacobi_diag(H,V);
Check that your energies are correct,
for (int k=0; k < n/3; k++){
double exact = M_PI*M_PI*(k+1)*(k+1);
double calculated = gsl_matrix_get(H,k,k);
printf("%i %g %g\n",k,calculated,exact);
}
Plot the several lowest eigenfunctions (and compare with the analytical result),
for(int k=0;k<3;k++){
printf("0 0\n");
for(int i=0;i<n;i++)
printf("%g %g\n",(i+1.0)/(n+1), gsl_matrix_get(V,i,k));
printf("1 0\n");
}
(1 point) Scaling, optimization, and comparison with GSL