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));
        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;
          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 outputs H as the identity matrix, which will be fed into tqli for the eigenvectors
//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.

  • $\begingroup$ 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. $\endgroup$ – davidhigh Mar 20 '20 at 19:43
  • $\begingroup$ @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. $\endgroup$ – user8384493 Mar 20 '20 at 23:16
  • $\begingroup$ 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. $\endgroup$ – Alone Programmer Mar 21 '20 at 0:30
  • $\begingroup$ @AloneProgrammer Educational purposes. $\endgroup$ – 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.
  • $\begingroup$ I already know it doesn't work. I stumped as to why. $\endgroup$ – 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.

  • $\begingroup$ 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. $\endgroup$ – user8384493 Mar 20 '20 at 23:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.