# Why NDEigensystem does not show the minimum eigenvalue?

I would like to find eigenvalues and eigenfunctions using NDEigensystem. The following system is considered: $$H[u(\rho,z)]=-\frac{1}{2}\Delta u(\rho,z)+Vu(\rho,z)$$ where $$V=-\frac{1}{\sqrt{\rho^2+z^2}}-\frac{3e^{-3\sqrt{\rho^2+z^2}}+e^{-7\sqrt{\rho^2+z^2}}}{\sqrt{\rho^2+z^2}}$$

Below are two codes, the first code (1) is to find eigenvalues and eigenfunctions using NDEigensystem. The second code (2) is finding eigenvalues and eigenfunctions by matrix method. In this case I can estimate the value of the minimum eigenvalue using minimization, because this approach uses parameters dependent basis functions.

(1) With the use of NDEigensystem, the minimum eigenvalue is -0.329556.
In the both codes I renamed $$\rho≡r$$

ClearAll["Global*"]

rmax = 30;
zmax = 30;

V[r_, z_] := -(1/
Sqrt[r^2 + z^2]) - (3 (E^(-3 Sqrt[r^2 + z^2]) +
E^(-7 Sqrt[r^2 + z^2])))/Sqrt[r^2 + z^2]

H = Simplify[-1/2*
Laplacian[u[r, z], {r, \[Theta], z}, "Cylindrical"] +
V[r, z]*u[r, z]];

{vals, funs} =
NDEigensystem[{H + u[r, z]},
u[r, z], {r, 0, rmax}, {z, -zmax, zmax}, 100,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}},
"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];

Sort[vals - 1]

(*{-0.329556, -0.129353, -0.0996647, -0.0570438, -0.0555859, \
-0.0478791, -0.0338522, -0.0332227, -0.0321545, -0.0306885, \
-0.0228879, -0.0221102, -0.0195824, -0.0178194, -0.0134078, \
-0.0131479, -0.00638181, -0.00461305, -0.00226629, 0.000615921, \
0.00246069, 0.00713658, 0.0121566, 0.0150791, 0.0164991, 0.0195553, \
0.0209365, 0.0263934, 0.0285584, 0.0335753, 0.0355955, 0.0369134, \
0.0418356, 0.045981, 0.0508021, 0.051836, 0.0567986, 0.0568928, \
0.059726, 0.0648629, 0.0671563, 0.0732491, 0.0753428, 0.0755007, \
0.0794506, 0.0828858, 0.0836577, 0.0935551, 0.0980983, 0.0986071, \
0.0988963, 0.103791, 0.105917, 0.110704, 0.116505, 0.117394, 0.12031, \
0.12446, 0.125365, 0.129507, 0.13085, 0.13276, 0.137021, 0.144945, \
0.146664, 0.148891, 0.15654, 0.159095, 0.160471, 0.162447, 0.168228, \
0.170975, 0.173943, 0.176716, 0.179559, 0.183812, 0.184116, 0.190538, \
0.192467, 0.194558, 0.197421, 0.199266, 0.207298, 0.207793, 0.213011, \
0.216477, 0.21986, 0.220111, 0.229101, 0.231261, 0.235627, 0.240813, \
0.241493, 0.24344, 0.247215, 0.24818, 0.250832, 0.255031, 0.258738, \
0.264149}*)



(2) Now let's move on to the matrix method and consider Gaussian functions as basis functions.

The basis set: $$\psi_{ij}=e^{-b_{j} z^2}e^{-a_{i} \rho^2}$$ , where $$a_{i}=a_1 * qa^{i-1}$$, $$i=1, 2 ,3,..., i_{max}$$ and $$b_{j}=b_1 * qb^{j-1}$$, $$j=1, 2 ,3,..., j_{max}$$ are geometrical progressions ($$a_1$$, $$b_1$$, $$qa$$, $$qb$$ are parameters).

In order to estimate the minimum eigenvalue, it is sufficient to take a small basis set. The basis set of just two functions (imax=1, jmax=2) already gives a minimum eigenvalue much lower (vals min=-4.00166) than that obtained using the NDEigensystem (vals min=-0.329556).

Why are the values so different? It seems to me that NDEigensystem doesn't output the smallest eigenvalue, but starts from the second one.

ClearAll["Global*"]

(*basis set*)

imax = 1; jmax = 2;

geoseq[init_, r_, n_] := init*r^(n - 1);

Psi[a1_, b1_, qa_, qb_, r_, z_, i_, j_] :=
Exp[-geoseq[b1, qb, j]*z^2]*Exp[-geoseq[a1, qa, i]*r^2];

(*\[Minus]1/2 \[CapitalDelta]*)

Kk[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] =
Integrate[
FullSimplify[
Psi[a1, b1, qa, qb, r, z, i2, j2]*
Laplacian[Psi[a1, b1, qa, qb, r, z, i1, j1], {r, \[Theta], z},
"Cylindrical"]] r, {r, 0, \[Infinity]}, {z, -Infinity,
Infinity},
Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0,
i2 > 0, j2 > 0}]
-((Sqrt[\[Pi]] (b1 qa (qa^i1 + qa^i2) qb^(j1 + j2) +
2 a1 qa^(i1 + i2) qb (qb^j1 + qb^j2)))/(
a1 (qa^i1 + qa^i2)^2 Sqrt[b1 qb] (qb^j1 + qb^j2)^(3/2)))

Kx[a1_, b1_, qa_, qb_,
c_List?MatrixQ] := -1/2 2 Pi Sum[
c[[i1, j1]] c[[i2, j2]] Kk[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1,
imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}]

(*V(r,z)*)

V[r_, z_] := -(1/
Sqrt[r^2 + z^2]) - (3 (E^(-7 Sqrt[r^2 + z^2]) +
E^(-3 Sqrt[r^2 + z^2])))/Sqrt[r^2 + z^2];

VP[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] :=
NIntegrate[
Psi[a1, b1, qa, qb, r, z, i2, j2]*V[r, z]*
Psi[a1, b1, qa, qb, r, z, i1, j1]*r, {r,
0, \[Infinity]}, {z, -Infinity, Infinity}];

Vx[a1_, b1_, qa_, qb_, c_List?MatrixQ] :=
2 Pi Sum[c[[i1, j1]] c[[i2, j2]] VP[a1, b1, qa, qb, i1, j1, i2,
j2], {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];

int[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] =
Integrate[
Psi[a1, b1, qa, qb, r, z, i2, j2]*
Psi[a1, b1, qa, qb, r, z, i1, j1] r, {r,
0, \[Infinity]}, {z, -Infinity, Infinity},
Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0,
i2 > 0, j2 > 0}];

norm[a1_, b1_, qa_, qb_, c_List?MatrixQ] :=
2 Pi Sum[c[[i1, j1]] c[[i2, j2]] int[a1, b1, qa, qb, i1, j1, i2,
j2], {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];

H[a1_, b1_, qa_, qb_,
c_List?MatrixQ] := (Kx[a1, b1, qa, qb, c] + Vx[a1, b1, qa, qb, c])/
norm[a1, b1, qa, qb, c];

(*finding the minimum eigenvalue*)

FindMinimum[{valsmin =
H[a1, b1, qa,
qb, {{c11, c12}}], (0.1 < a1 < 10) && (0.1 < b1 < 10) && (1.1 <
qa < 10) && (1.1 < qb < 10) && (0.1 < c11 < 10) && (0.1 < c12 <
10)}, {{a1, 1.22}, {b1, 0.44}, {qa, 1.40}, {qb, 1.81}, {c11,
1}, {c12, 1}}, AccuracyGoal -> 3, PrecisionGoal -> 3,
EvaluationMonitor :> {Print["a1=", a1, "   b1=", b1, "   qa=", qa,
"   с=", {{c11, c12}}, "   vals min=", valsmin]}]

(*{-4.00166, {a1 -> 8.21271, b1 -> 4.56641, qa -> 3.79284,
qb -> 5.76842, c11 -> 5.0677, c12 -> 3.29988}}*)

• Using Method -> {"Eigensystem" -> {"Direct"}} yields a lowest eigenvalue of -5.61306, so it may be that the problem lies in some details of the Arnoldi solver. Jun 9, 2023 at 20:54
• @MamMam Please, check, how your code works for Coulomb potential V[r_, z_] := -(1/ Sqrt[r^2 + z^2]). Jun 10, 2023 at 3:08
• @MamMam Thank you. This is good. On the next step please use V[r_, z_] := -(1/ Sqrt[r^2 + z^2]) - eps (3 (E^(-3 Sqrt[r^2 + z^2]) + E^(-7 Sqrt[r^2 + z^2])))/Sqrt[r^2 + z^2] and plot Min[vals] for {eps,0,1/3,1/30}. Let compare it with the same plot for your code with a, b. Jun 10, 2023 at 7:44
• @MamMam Perfect! On the next step please compute first and second code on {eps, 0, 2/3, 1/30}. You will see that NDEigensystem can't compute ground state at eps>=7/15. Instead it compute second state. Jun 10, 2023 at 9:41
• @MamMam To compute up to eps=0.7 use code UT[r_] := -(1/r) - eps (3 (E^(-3 r) + E^(-7 r)))/r; RR = Simplify[-1/2*Laplacian[u[r], {r, \[Theta], \[Phi]}, "Spherical"] + UT[r]*u[r]]; sol[x_] := NDEigensystem[{RR /. eps -> x}, u[r], {r, 0, 12}, 10, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> .001}}}]; lst = Table[{x, Min[sol[x][[1]]]}, {x, 0, 1, 1/30}] Jun 10, 2023 at 10:30

For this case can be used second code from my answer here for the Slater type orbitals. First we define potential in the form

VP1[r_] := -1/r - eps (3 (E^(-3 r) + E^(-7 r)))/r;


where 0<=eps<=1 is parameter. Code is given by

ClearAll["Global*"]
nmax = 5; var[n_] := Table[c[n, l], {l, nmax}];
Psi[r_, n_] := r^(n - 1) Exp[-r];
(*kinetic energy*)
Kk[r_, n1_, n2_] :=
FullSimplify[
Psi[r, n2]*
Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n_] := -1/2 Sum[
c[n, l1] c[n, l2] Integrate[
Kk[r, l1, l2]*r^2, {r, 0, \[Infinity]}], {l1, nmax}, {l2, nmax}];
(*potential energy*)

VP1[r_] := -1/r - eps (3 (E^(-3 r) + E^(-7 r)))/r;
Px1[n1_, n2_] :=
Integrate[Psi[r, n2]*VP1[r]*Psi[r, n1]*r^2, {r, 0, \[Infinity]}];
Px[n_] := Sum[c[n, l1] c[n, l2] Px1[l1, l2], {l1, nmax}, {l2, nmax}];

norm[n_] = {Sum[
c[n, l1] c[n, l2] NIntegrate[
Psi[r, l1] Psi[r, l2] r^2, {r, 0, \[Infinity]}], {l1,
nmax}, {l2, nmax}] == 1};

sol[n_, x_] :=
NMinimize[{Kx[n] + Px[n] /. eps -> x, norm[n]}, var[n]];


With this code we compute data

lst5 = Table[{x, sol[1, x][[1]]}, {x, 0, 1, .1}]


Data lst5 we can compare with 2D model same as above

rmax = 19;
zmax = 6;

V[r_, z_] := -(1/Sqrt[r^2 + z^2]) -
eps 3 (E^(-3 Sqrt[r^2 + z^2]) + E^(-7 Sqrt[r^2 + z^2]))/
Sqrt[r^2 + z^2]

H = Simplify[-1/2*
Laplacian[u[r, z], {r, \[Theta], z}, "Cylindrical"] +
V[r, z]*u[r, z]];

sol2[x_] :=
NDEigensystem[{H /. eps -> x},
u[r, z], {r, 0, rmax}, {z, -zmax, zmax}, 60];

lst2 = Table[{x, sol2[x][[1]] // Min}, {x, 0, 1, 1/30}]


Note, that to compute the ground state we can transform 2D model in 1D model as follows

UT[r_] := -(1/r) - eps (3 (E^(-3 r) + E^(-7 r)))/r;
RR = Simplify[-1/2*
Laplacian[u[r], {r, \[Theta], \[Phi]}, "Spherical"] + UT[r]*u[r]];

sol1[x_] :=
NDEigensystem[{RR /. eps -> x}, u[r], {r, 0, 12}, 10,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> .001}}}];
lst1 = Table[{x, Min[Sort[sol1[x][[1]]]]}, {x, 0, 1, 1/30}]


Finally we compare ground state computed with NMinimize (solid line) and with NDEigensystem in 1D (red points), and in 2D (blue 0)

Show[ListLinePlot[lst5, Frame -> True],
ListPlot[lst1, PlotStyle -> Red],
ListPlot[lst2, PlotStyle -> Blue, PlotMarkers -> "0"]]


From this picture we see similar jump from the ground state to the first excited state in 1D and 2D model as well. The reason is that Arnoldi method used as default is unstable. We can try option Method -> {"Eigensystem" -> "Direct"}. But with this method we compute ground state with lower accuracy,

UT[r_] := -(1/r) -
eps (3 (E^(-3 r) + E^(-7 r)))/r; oper = -1/2/r^2 D[r^2 D[u[r], r],
r] + UT[r] u[r];

sol[x_] :=
NDEigenvalues[{oper /. eps -> x}, u[r], {r, 0, 12}, 20,
Method -> {"Eigensystem" -> "Direct"}]; lstd =
Table[{x, Min[sol[x]]}, {x, 0, 1, 1/30}]


Visualization lstd(black dashed line) with data shown above

• Thanks a lot for the detailed description! Could you please write down the expression that uses this method Method -> {"Eigensystem" -> "Direct"} in full. Because when I write like this {vals, funs} = NDEigensystem[{H + u[r, z]}, u[r, z], {r, 0, rmax}, {z, -zmax, zmax}, 100, Method -> {"Eigensystem" -> "Direct"}], Mathematica outputs the following eigenvalues {-0.222278 + 0. I, -0.125203 + 0. I, -0.0790437 + 0. I, -0.0555639 + 0. I, -0.0555584 + 0. I, -0.0413612 + 0. I. That is, the ground state is also skipped Jun 11, 2023 at 7:53
• @MamMam Please, see update to my answer. Jun 11, 2023 at 10:26
• @MamMam Actually this expression 1/r^2 D[r^2 D[u[r], r], r] is the Laplacian in spherical coordinates. So there is no any typo in my expression oper. Please pay attention that 4 methods give us -0.5 at eps=0. Jun 11, 2023 at 11:14
• @MamMam Actually in my answer there are 4 models and one of them is 2D model in cylindrical coordinates - blue 0 markers. As I understand, you try to test your model with a, b it is why you using cylindrical coordinates. But your potential depends on r in spherical coordinates. Therefore we can compute the ground state in the spherical coordinates. Jun 11, 2023 at 11:58
• @MamMam Please pay attention that model for the Slater type orbitals is the best one in this case, and your model with a, b is also good. But models with NDEigensystem` look very poor. Jun 11, 2023 at 12:49