Generalized eigenvalue problem

1. Implement a function to solve a generalized eigenvalue problem AV=BVE (where A is a real symmetric matrix, B is a real symmetric positive definite matrix, V is the matrix of the generalized eigenvectors and E is the diagonal matrix with the generalized eigenvalues) via diagonalization of the matrix B. Specifically, the diagonalization of matrix B, given as B=QSQ^T, where S is a diagnoal matrix with positive elements and Q is an orthogonal matrix, turns the eigenproblem

   AV=BVE
   

   into

   AV=QSQTVE .
   

   This can be identically rewritten as

   S-½QTAQS-½ S½QTV = S½QTV E ,
   

   which is an ordinary real symmetric eigenvalue problem

   à X = X E ,
   

   where

   à = S-½QTAQS-½ ,
   

   and

   X = S½QTV .
   

   Diagonalization of matrix à gives directly the eigenvalues E while the original eigenvectors can be calculated as

   V = QS-½X .
   

   (Note that since S is diagonal so is S^½ and S^-½ and therefore you don't need to build this matrices explicitly and waste O(n³) operations for matrix multiplication.)

2. Find the ground state of the hydrogen atom with the Schrödinger equation

   (-1/2 d²/dr² - 1/r)u(r) = εu(r),
   

   using the variational method,

   u(r) = ∑i=1..nciφi(r)
   

   with the basis functions φ_i(r) in the form (satisfying the boundary conditions u(0)=u(∞)=0)

   φi(r) = r e-αir²

   where α_i are the variational parameters (to be found by your own minimization routine) and where the coefficients c_i should be found by solving the generalized eigenvalue problem

   H c = ε N c
   

   where H and N are matrices with the elements

   Hij = ∫0∞dr φi(r) (-1/2 d²/dr² - 1/r) φj(r) ,
   Nij = ∫0∞dr φi(r) φj(r) .
   

   The matrix elements are actually analytic (check the formulae before using),

   ∫0∞dr r e-αr² r e-βr² = 1/4 √π (α+β)-3/2 ,
   

   ∫0∞dr r e-αr² (d²/dr²) r e-βr² = -3/2 √π αβ (α+β)-5/2 ,
   

   ∫0∞dr r e-αr² (1/r) r e-βr² = 1/2 (α+β)-1 .
   

   Try n=1,2,3,... gaussians in your basis and investigate the convergence. The method will probably break down at about 5 Gaussians due to their excessive overlap.

3. Implement a function to solve the generalized eigenvalue problem via Cholesky decomposition of the matrix B.