6
$\begingroup$

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}}*)
Improve this question
$\endgroup$
14
  • 2
    $\begingroup$ 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. $\endgroup$
    – Michael Seifert
    Jun 9 at 20:54
  • 1
    $\begingroup$ @MamMam Please, check, how your code works for Coulomb potential V[r_, z_] := -(1/ Sqrt[r^2 + z^2]). $\endgroup$
    – Alex Trounev
    Jun 10 at 3:08
  • 1
    $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Jun 10 at 7:44
  • 1
    $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Jun 10 at 9:41
  • 1
    $\begingroup$ @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}] $\endgroup$
    – Alex Trounev
    Jun 10 at 10:30

1 Answer 1

Reset to default
5
$\begingroup$

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"]]

Figure 1

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 Figure 2

Improve this answer
$\endgroup$
12
  • $\begingroup$ 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 $\endgroup$
    – Mam Mam
    Jun 11 at 7:53
  • 1
    $\begingroup$ @MamMam Please, see update to my answer. $\endgroup$
    – Alex Trounev
    Jun 11 at 10:26
  • 1
    $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Jun 11 at 11:14
  • 1
    $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Jun 11 at 11:58
  • 1
    $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Jun 11 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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