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}}*)
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. $\endgroup$V[r_, z_] := -(1/ Sqrt[r^2 + z^2])
. $\endgroup$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 witha, b
. $\endgroup${eps, 0, 2/3, 1/30}
. You will see thatNDEigensystem
can't compute ground state ateps>=7/15
. Instead it compute second state. $\endgroup$eps=0.7
use codeUT[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}]
$\endgroup$