I'm numerically solving a time-independent Schrödinger equation using Eigensystem's FEAST method. It takes a lot longer than I would have thought; 63 seconds to grab the first 8 eigenvectors of a very sparse $1608\times1608$ matrix.

Is this just the kind of performance I should expect in the problem of numerically diagonalizing a matrix? Even if I only want a few low eigenvalues and the matrix is very sparse? Or is there something in my implementation that is causing this to be very inefficient?

Here is the code in its entirety.

These functions are used to construct the hamiltonian matrix

laplacian[xmin_, xmax_, dx_] := 
 SparseArray[{Band[{1, 1}] -> -2, Band[{1, 2}] -> 1, 
 Band[{2, 1}] -> 1}, {Floor[(xmax - xmin)/dx], 
 Floor[(xmax - xmin)/dx]}]/dx^2 // N

potential[xmin_, xmax_, dx_, f_] := 
 SparseArray[{Band[{1, 1}] -> 
 Table[f[x], {x, xmin, xmax - dx, dx}]}, {Floor[(xmax - xmin)/dx],
 Floor[(xmax - xmin)/dx]}] // N

The actual potential I'm using is a 1d quartic potential selected so that the system will have four eigenvectors of negative energy. Sorry if this looks bad on stackexchange; it should be readable in a notebook though.

V[x_]:=-(m \[Omega]^2  x^2/2)(\[HBar]/(2 m \[Omega]))+(\[Epsilon]x^4/4)(\[HBar]/(2 m \[Omega]))^2

m=1; \[Omega]=1; \[Epsilon]=0.1; \[HBar]=1; E0=0;\[Delta]E=10;dx=0.01;
minV = Minimize[V[x],x][[1]];

Then I chop off the matrix wherever the potential is sufficiently large.

turningpoints = NSolve[{10 == V[x], x \[Element] Reals}, x]
qmax = Max[x /. turningpoints];
qmin = Min[x /. turningpoints];

Here is the hamiltonian (and I check that it's hermitian so we can use FEAST). It is very sparse.

ham = -\[HBar]^2/(2 m)*(2 m \[Omega])/\[HBar] laplacian[qmin,qmax,dx]+potential[qmin,qmax,dx,V];

But still Eigensystem takes a really long time to solve for its eigenvectors

{time, {energies, states}} = 
Method -> {"FEAST", MaxIterations -> 1000, 
  Interval -> {minV, -minV}}];

Your Answer

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

Browse other questions tagged or ask your own question.