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];
MatrixPlot@ham
HermitianMatrixQ@ham
But still Eigensystem
takes a really long time to solve for its eigenvectors
{time, {energies, states}} =
Timing@Eigensystem[ham,
Method -> {"FEAST", MaxIterations -> 1000,
Interval -> {minV, -minV}}];