# 2d Schrodinger Equation via matrix diagonalization in C

I am trying to solve the time-independent Schrodinger equation in two dimensions via discrete matrix diagonalization. I want the energy eigenvalues and the corresponding eigenfunctions for a given potential.

Eventually I want the solutions for the following potential (although as shown later I can't seem to get it working for any potential): $$V(x,y) = x^{4} + y^{4}$$

To discretize the problem I use the following basis functions: $$\phi_{i}(x,y) = \begin{cases} \frac{1}{h^{2}} & \frac{-h}{2}\leq x-x_{i}<\frac{h}{2}, \frac{-h}{2}\leq y-y_{i}<\frac{h}{2} \\ 0 & \text{else} \end{cases}$$ Essentially, each function occupies a single grid square with a normalized height of $$\frac{1}{h^{2}}$$ where $$h$$ is the width of a square. This is a good basis to use because if I can find the eigenvectors expressed in this basis, it should be trivial to plot them. To map the single index $$i$$ to the 2d grid I use the following scheme:

$$x_{i} = h((i\mod N) - R) \\ y_{i} = h\left(\left\lfloor\frac{i}{N}\right\rfloor - R\right) \\ N = 2R + 1$$ where R is an integer 'radius'. This is to say the domain is an $$N\times N$$ grid centered at the origin. The ordering starts at $$(-R,-R)$$ and reads until $$(R,-R)$$ before advancing to the above row at $$(-R,-R+1)$$ and so on and so forth.

The matrix for the kinetic term of the Hamiltonian will be constructed from the Laplacian matrix operator. I use $$\nabla^{2} = \begin{pmatrix} -4 & 1 & 0 & 1 \\ 1 & -4 & 1 & \ddots & \ddots \\ 0 & 1 & \ddots & \ddots & \ddots & 1\\ 1 & \ddots & \ddots & \ddots & 1 & 0 \\ & \ddots & \ddots & 1 & -4& 0 \\ & & 1 & 0 & 0 & -4 \\ \end{pmatrix}$$ Each row in the grid domain is punctuated by zeros on the sub and super-diagonals, then the pattern repeats. The potential matrix is diagonal given by $$V_{ij} = \delta_{ij} \langle i | V | j\rangle$$

Thus the matrix Hamiltonian: $$\mathbf{H} = -\frac{\hbar^{2}}{2m}\nabla^{2} + \mathbf{V}$$

Currently, $$\mathbf{H}$$ is not tridiagonal. So first I use Numerical Recipe's tred2 routine. This outputs the diagonal terms into a 1d array and the off-diagonals into another array. I then pass these two arrays into Numerical Recipe's tqli routine for the diagonalization. This routine outputs the eigenvalues by overwriting the diagonal array (d[k] for $$k$$-th eigenvalue) and outputs the full set of eigenvectors as columns into a 2d array.

The output I am getting now is completely nonsensical. My understanding was that the output array would have eigenvectors expressed in the original basis, so that I could visualize them via a heat map or the like with gnuplot. Am I missing a step, or is my code just erroneous?

Earlier I accidentally got it working for the 1d case. I was using the wrong kinetic matrix (1d case is already tridiagonal) but with the 2d potential and some of the eigenvectors I saw actually resembled realistic wavefunctions. So, if I didn't touch tqli since then then the problem isn't there.

Since I am unsure of what or where the problem is, I am going to provide the full program.

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

//L is the unit width of the domain.
//R is the integer radius.  Number of states is (1+2R)^2
//m is particle mass.
#define L 2.0
#define R 15
#define m 1.0

double V(int,double);
double x(int,double);
double y(int,double);
void tred2(double**,int,double*,double*);
void tqli(double*,double*,int,double**);
double pythag(double a,double b);

double x(int i,double h) {return h*((i%(1+2*R) - R));}
double y(int i,double h) {return h*(i/(1+2*R) - R);}

//potential function. returns the integral over the grid square indexed by i
double V(int i,double h) {
double xi=pow(x(i,h),2),yi=pow(y(i,h),2),v;
v = pow(h,5)/40.0;
v += h*h*(xi*xi+yi*yi)/2.0;
v += pow(xi,4) + pow(yi,4);
return v;}

double pythag(double a, double b) {
double A=fabs(a),B=fabs(b);
if (A>B) {return A*sqrt(1.0+pow(B/A,2));}
return (B==0.0?0.0:B*sqrt(1.0+pow(A/B,2)));
}

void tred2(double **a,int n,double *d,double *e) {
int l,k,j,i;
double scale,hh,h,g,f;

for (i=n-1;i>0;i--) {
l = i - 1;
h = scale = 0.0;
if (l>0) {
for (k=0;k<=l;k++) {scale += fabs(a[i][k]);}
if (scale==0.0) {e[i] = a[i][l];}
else {
for (k=0;k<=l;k++) {
a[i][k] /= scale;
h += a[i][k]*a[i][k];
}
f = a[i][l];
g = (f>=0.0?-sqrt(h):sqrt(h));
e[i]=scale*g;
h -= f*g;
a[i][l] = f-g;
f = 0.0;
for (j=0;j<=l;j++) {
a[j][i] = a[i][j]/h;
g = 0.0;
for (k=0;k<=j;k++) {g += a[j][k]*a[i][k];}
for (k=j+1;k<=l;k++) {g += a[k][j]*a[i][k];}
}
}
d[i] = a[i][i];
a[i][i] = 1.0;
for (j=0;j<=l;j++) {a[j][i] = a[i][j] = 0.0;}
}
}

void tqli(double *d,double *e,int n,double **z) {
int j,l,it,i,k;
double s,r,p,g,f,dd,c,b;

for (i=2;i<n;i++) {e[i-1] = e[i];}
e[n-1] = 0.0;
for (l=0;l<n;l++) {
it = 0;
do {
for (j=l;j<n-1;j++) {
dd = fabs(d[j])+fabs(d[j+1]);
if ((double)(fabs(e[j])+dd) == dd) break;
}
if (j!=l) {
if (it++ == 30) {printf("Too many iterations in tqli\n");}
g = (d[l+1]-d[l])/(2.0*e[l]);
r = pythag(g,1.0);
g = d[j]-d[l]+e[l]/(g+(g<0.0)?-r:r);
s = c = 1.0;
p = 0.0;
for (i=j-1;i>=l;i--) {
f = s*e[i];
b = c*e[i];
r = pythag(f,g);
e[i+1] = r;
if (r == 0.0) {
d[i+1] -= p;
e[j] = 0.0;
break;
}
s = f/r;
c = g/r;
g = d[i+1] - p;
r = (d[i]-g)*s + 2.0*c*b;
d[i+1] = g + (p = s*r);
g = c*r - b;
for (k=0;k<n;k++) {
f = z[k][i+1];
z[k][i+1] = s*z[k][i] + c*f;
z[k][i] = c*z[k][i] - s*f;
}
}
if (r == 0.0 && i >= l) continue;
d[l] -= p;
e[l] = g;
e[j] = 0.0;
}
} while (j != l);
}
}

void main() {
int i,j,N=1+2*R,M,a=0;
double **H,h=L/N,dd,*D,*T,q;
FILE *f;
dd = pow(h,-2);
M = N*N; //M is total number of states

D = (double *) malloc(M*sizeof(double));
T = (double *) malloc(M*sizeof(double));
H = (double **) malloc(M*sizeof(double*));
H[0] = (double *) malloc(M*M*sizeof(double));
for (i=1;i<M;i++) {
H[i] = H[i-1] + M;
}

//initialize Hamiltonian in discrete basis.  triangular loop by symmetry.
for (i=0;i<M;i++) {
*(H[i]+i) = V(i,h) + 2*dd/m;
for (j=i+1;j<M;j++) {
*(H[i]+j) = 0.0;
if (j==i+1 && i%N!=N-1) {*(H[i]+j) = -dd/2;}
else if (j==i+3) {*(H[i]+j) = -dd/2;}
*(H[j]+i) = *(H[i]+j);
}
}

//D and T are empty, but tred2 will write into them. D for diagonals, T for the off-diagonal terms
tred2(H,M,D,T);
//tred2 outputs H as the identity matrix, which will be fed into tqli for the eigenvectors
tqli(D,T,M,H);
//D[k] now contains kth eigenvalue.  H columns are eigenvectors
//you may decide how you want to visualize the output
}


The original routines can be found in the eigensystems section in http://www2.geog.ucl.ac.uk/~mdisney/teaching/GEOGG121/inversion/NumericalRecipesinC.pdf Note: they use arrays of a[1...n] so I converted the functions so they would work on a[0...n-1] indexed arrays.

I have tried many different values for size and $$h$$ but all the output is garbage. A typical output is zero in most places but with completely random valued stripes or spots without any kind of symmetry.

• If its about solving the TDSE for this particular potential, this can be better done in radial coordinates. You should then be able ro treat the angular part analytically. – davidhigh Mar 20 '20 at 19:43
• @davidhigh, Ah, that might help for verifying the solutions, but I need to solve this problem numerically, specifically through matrix diagonalization.. Besides, if I can get this to work it will also work with any arbitrary potential. – user8384493 Mar 20 '20 at 23:16
• Just a comment: Do you insist to solve it in this way for probably educational purposes or you have some freedom to choose whatever works for you? I think there are lots of more straightforward ways to do this. – Alone Programmer Mar 21 '20 at 0:30
• @AloneProgrammer Educational purposes. – user8384493 Mar 21 '20 at 17:49

1. Test your code. "I don't know if my linear algebra routines work or not" is a problem you can solve easily. It is easy to check if you have computed the correct eigenvalues or not; just check that $$VDV^{-1}=H$$, or if you don't fancy the inversion $$HV=VD$$. If you are not sure that your individual pieces work, you are walking in the dark.
2. Possibly a hot take: use a different language. With C all code takes longer to write (and read) with respect to a language that has more numerical libraries and some more syntactic sugar for faster developing/prototyping. You are not even sure if your discretization strategy is sound here; it is not the time to use the fastest-language-in-the-world for performance reasons. Just do your early development and debugging in a friendlier environment. Chances are that you won't even need to move to C for performance reasons in this problem, since 90% of your CPU time will be spent in a linear algebra library anyway.
• I already know it doesn't work. I stumped as to why. – user8384493 Mar 21 '20 at 17:48

Your whole code is not understandable to me. Try this alternative components:

int grid(int i, double a, double h)
{
return a+i*h;
}

int index(int i, int j, int N)
{
return i*N+j;
}

//potential function. returns the integral over the grid square indexed by i
double V(int i, double x0, double hx, int j, double y0, double hy)
{
double xi=grid(i,x0,hx);
double yi=grid(j,y0,hy);
return pow(xi,4) + pow(yi,4);
}


Then use the following code to initialize your Hamiltonian:

//initialize Hamiltonian in discrete basis.  no triangular loop at first!
int M2=M*M;
//set Hamiltonian to zero
for (i=0;i<M;i++)
{
for (j=0;j<M;j++)
{
H((i*M+i)*M2+(j*M+j)) = -2.0/(hx*hx) - 2.0/(hy*hy) + V(i,x0,hx,j,y0,hx);
H((i*M+(i+1))*M2+(j*M+j)) = 1.0/(hx*hx);
H((i*M+(i-1))*M2+(j*M+j)) = 1.0/(hx*hx);
H((i*M+i)*M2+(j*M+(j+1))) = 1.0/(hy*hy);
H((i*M+i))*M2+(j*M+(j-1))) = 1.0/(hy*hy);
}
}


Note that this code will probably not run, because I'm not a C guy. but it should give you the idea. Notably, it's not necessary to do the modulo-and-integer-division index hackery.

• The grid is not the problem. I have tested that part and it returns what I want, something is happening down the line in either tred2 or tlqi. – user8384493 Mar 20 '20 at 23:13