Solving Schrodinger Equation - Linear Algebra with LAPACK

Today I solve the Time-independent Schrodinger equation (TISE) using the finite difference method that I explained in previous blog and the Matrix method. Being able to solve the TISE numerically is important since only small idealized system can be solved analytically. I consider here only one dimensional case for a particle in a potential V(x). The one dimensional TISE is

-(1/2) [d^2 psi(x) /dx^2] + V(x) psi(x) = E psi(x)

Note I use natural unit m = hbar = 1.

The first step we need to do is to transform differential equation into matrix equation. To use finite difference, we evaluate the wavefunction, psi(x) i.e. the solution of the TISE, at N+1 grid points with spacing dx, {psi(i)}. Similarly the potential is evaluated at i*dx i.e. V(i). Note i=0, 1, ..., N. The computational length is L = N*dx.

Using central finite difference,

d^2 psi(x)/dx^2 = (psi(x+dx) - 2psi(x) + psi(x-dx))/dx^2
= ( psi(i+1) - 2psi(i) + psi(i-1) )/dx^2

This transforms the 1D TISE into
-( psi(i+1) - 2psi(i) + psi(i-1) )/(2dx^2) + V(i)psi(i) = E psi(i)

after simple manipulation, we have

-psi(i+1) + (2 + (2dx^2)V(i))psi(i) - psi(i-1) = (2dx^2)E psi(i)

For computational convenience we use the boundary conditions where psi(0)=psi(L)=0;

The above discretized 1D TISE and using boundary conditions, in matrix form, is given by





This is an eigen matrix problem. To solve it I use LAPACK function dstev_() to compute eigenvalues for symmetric matrix. As an example I consider a finite square well potential given by V(x) = 1 for |x-5|>1 and V(x) = 0 for |x-5|<1 and grid spacing dx=0.1.


--- eigen.c ---

#include <stdio.h>
#include <math.h>

#define N 100

int main(void)
{
int i;
double dx = 0.1;
double dxsq = dx*dx;
double hdiag[(N-1)],hoff[(N-2)]; //hold diagonal and off-diagonal
double v[(N-1)];
char jobType = 'V';
long int n = (N-1);
long int info,ldz = N-1;
double work[2*N-4];
double z[(N-1)][(N-1)];

/* initialize V(x) */
for(i=0;i<(N-1);i++){
if(fabs((i+1)*dx - 5)<1){ v[i]= 0;}
else v[i] = 1;
}

/* initialize hoff */
for(i=0;i<(N-2);i++){
hoff[i] = -1;
}
/* initialize hdiag */
for(i=0;i<(N-1);i++){
hdiag[i] = 2.0 + 2.0*dxsq*v[i];
}

dstev_(&jobType,&n,hdiag,hoff,z,&ldz,work,&info);

/* Print results for four eigen states */
printf("info = %ld\n",info);
printf("The Lowest Four Eigen Values \n");
for(i=0;i<4;i++){
printf("%ld %lf\n",i,hdiag[i]/(2.0*dxsq));
}

printf("The Lowest four Eigen Vectors\n");
printf("x groundstate 1st 2nd 3rd \n");
for(i=0;i<(N-1);i++){
printf("%lf %lf %lf %lf %lf\n",i*dx,z[0][i],z[1][i],z[2][i],z[3][i]);
}

return 0;
}



To compile:

gcc eigen.c -lm -llapack -o try