# Eigenvalue problem

my question is about solving an eigenvalue problem of the Helmholtz equation using sinc approximation
$$\nabla^2E + V (x) = \lambda E$$ and $$V(x)= X^2 / 2$$

I have a problem in calculating the eigenvalues of this system because the eigen values if this system is known to be $$\lambda$$=$$\{1/2 , 3/2 , 5/2\}$$

my code is :

Nn = 8
h = π/Sqrt[2 Nn]
h1 = h

exactdelta2 = h1^2 D[S1[ϕ1[x], j, h1], { ϕ1[x], 2}]

ex5 = exactdelta2 /. x -> ψ1[x] /. ϕ1[ψ1[x_]] :> x /. x -> k h1

ex6 = ex5 /. S1 -> Function[{x, k, h}, Sinc[π (x - k h)/h]]

nabla = Table[ex6, {k, -Nn, Nn}, {j, -Nn, Nn}] /.Indeterminate :> Limit[ex6, k -> j]

vx[x_] = x^2/2 /. x -> k h

vx1 = Table[vx[x], {k, -Nn, Nn}]

vxmat = DiagonalMatrix[vx1]

A = nabla + vxmat // N

lamda = Eigenvalues[A]


Weak method, the solution converges to the exact one at Nn >= 200. It takes time. Figure 1 shows the solution for Nn=100, 200 in comparison with the exact solution.

Nn = 200;
h = \[Pi]/Sqrt[2 Nn];
h1 = h;

exactdelta2 = h1^2 D[S1[\[Phi]1[x], j, h1], {\[Phi]1[x], 2}];

ex5 = (exactdelta2 /. x -> \[Psi]1[x]) /. \[Phi]1[\[Psi]1[x_]] :>
x /. x -> k h1;
f[x_, k_, h_] :=
If[Abs[x - k h] > 0, Sin[\[Pi]*(x - k h)/h]/(\[Pi]*(x - k h)/h), 1]
ex6 = ex5 /. S1 -> f;

nabla = Table[ex6, {k, -Nn, Nn}, {j, -Nn, Nn}];

vx[x_] := -x^2/2;

vx1 = Table[vx[x] /. x -> k h, {k, -Nn, Nn}];

vxmat = DiagonalMatrix[vx1];

A = nabla + vxmat // N;

lamda = Reverse[Eigenvalues[A, -15]]/h1*Sqrt[2]

(*Out[]= {0.61922, -1.38078, 2.61922, -3.38078, 4.61922, -5.38078, \
6.61922, -7.38078, 8.61922, -9.38078, 10.6192, -11.3808, 12.6192, \
-13.3808, 14.6192}*)

Show[
ListPlot[{Table[1/2 + i, {i, 0, 14}], Abs[lamda]},
PlotMarkers -> Automatic, PlotLegends -> {"Exact", "Numeric"},
PlotLabel -> Row[{"Nn = ", Nn}]], Plot[x - .5, {x, 0, 15}]]


The standard method is much more accurate and faster.

 h = 1/10; {vals, func} =
NDEigensystem[{-h^2*Laplacian[u[x], {x}] + x^2/2*u[x],
DirichletCondition[u[x] == 0, True]}, u[x], {x, -3, 3}, 15,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {MaxCellMeasure -> 0.05}}}}];

vals/h/Sqrt[2]

(*Out[]= {0.5, 1.5, 2.50001, 3.50003, 4.50005, 5.50009, 6.50015, \
7.50023, 8.50033, 9.50046, 10.5006, 11.5008, 12.501, 13.5013, 14.5016}*)

Show[
ListPlot[{Table[1/2 + i, {i, 0, 14}], vals/h/Sqrt[2]},
PlotLegends -> {"Exact", "Numeric"}, PlotMarkers -> Automatic],
Plot[x - .5, {x, 0, 15}]]


The author's code slows down when calculating nabla. To reduce the time by 2-3 times, we applied Compile[] and ParallelTable[]. Figure 3 shows the numerical solution for Nn=111, 137 in comparison with the exact solution.

Nn = 137;
h = \[Pi]/Sqrt[2 Nn];
h1 = h;

exactdelta2 = h1^2 D[S1[\[Phi]1[x], j, h1], {\[Phi]1[x], 2}];

ex5 = (exactdelta2 /. x -> \[Psi]1[x]) /. \[Phi]1[\[Psi]1[x_]] :>
x /. x -> k h1;

cf = Compile[ {{x, _Real}, {k, _Integer}, {h, _Real}},
If[Abs[x - k h] > 0, Sinc[\[Pi]*(x - k h)/h], 1],
CompilationTarget -> "C"];

ex6 = ex5 /. S1 -> cf;
nabla = ParallelTable[N[ex6], {k, -Nn, Nn}, {j, -Nn, Nn}];

vx[x_] := -x^2/2;

vx1 = Table[vx[x] /. x -> k h, {k, -Nn, Nn}];

vxmat = DiagonalMatrix[vx1];

A = nabla + vxmat;

eigs = Eigensystem[A];
lambd = Reverse[Select[Abs[First[eigs]]/h1*Sqrt[2], # < 30 &]]
(*{0.485734, 1.51427, 2.48573, 3.51427, 4.48573, 5.51427, 6.48573, \
7.51427, 8.48573, 9.51427, 10.4857, 11.5143, 12.4857, 13.5143, \
14.4857, 15.5143, 16.4857, 17.5143, 18.4857, 19.5143, 20.4857, \
21.5143, 22.4857, 23.5143, 24.4857, 26.4858, 28.4854}*)
Show[ListPlot[{Table[1/2 + i, {i, 0, 29}], Abs[lambd]},
PlotMarkers -> Automatic, PlotLegends -> {"Exact", "Numeric"},
PlotLabel -> Row[{"Nn = ", Nn}]], Plot[x - .5, {x, 0, 30}]]


It was possible to improve the algorithm, improve accuracy and speed

exactdelta2 = h1^2 D[S1[\[Phi]1[x], j, h1], {\[Phi]1[x], 2}];

ex5 = exactdelta2 /. x -> \[Psi]1[x] /. \[Phi]1[\[Psi]1[x_]] :> x /.
x -> k h1;

ex6 = ex5 /. S1 -> Function[{x, k, h}, Sinc[\[Pi] (x - k h)/h]];

Nn = 137;
h = \[Pi]/Sqrt[2 Nn];
h1 = h; nabla =
ParallelTable[
If[k != j, N[ex6], -(\[Pi]^2/3)], {k, -Nn, Nn}, {j, -Nn, Nn}];
vx[x_] = -x^2/2 /. x -> k h;

vx1 = Table[vx[x], {k, -Nn, Nn}];

vxmat = DiagonalMatrix[vx1];

A = nabla + vxmat;
eigs = Eigensystem[A];
lambd = Reverse[Select[Abs[First[eigs]]/h/Sqrt[2], # < 30 &]]
Show[ListPlot[{Table[1/2 + i, {i, 0, 29}], Abs[lambd]},
PlotMarkers -> Automatic, PlotLegends -> {"Exact", "Numeric"},
PlotLabel -> Row[{"Nn = ", Nn}]], Plot[x - .5, {x, 0, 30}]]
(*{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, \
13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, \
24.5, 25.5, 26.4998, 27.5004, 28.4986, 29.5033}*)


Find the eigenvalues of the Laplacian on the segment -1<=x<=1

exactdelta2 = h1^2 D[S1[\[Phi]1[x], j, h1], {\[Phi]1[x], 2}];

ex5 = exactdelta2 /. x -> \[Psi]1[x] /. \[Phi]1[\[Psi]1[x_]] :> x /.
x -> k h1;

ex6 = ex5 /. S1 -> Function[{x, k, h}, Sinc[\[Pi] (x - k h)/h]];
Nn = 137;
h = 2/(2*Nn + 1);
h1 = h; nabla =
ParallelTable[
If[k != j, N[ex6], N[-(\[Pi]^2/3)]], {k, -Nn, Nn}, {j, -Nn, Nn}];
vx[x_] = 0;

vx1 = Table[0, {k, -Nn, Nn}];

vxmat = DiagonalMatrix[vx1];

A = nabla + vxmat;

eigs = Eigensystem[A];
lambd = Reverse[Select[Abs[First[eigs]]/h^2, # < 2300 &]]
Show[ListPlot[{Table[i^2*Pi^2/(h*(2*Nn + 1))^2, {i, 1, 29}],
Abs[lambd]}, PlotMarkers -> Automatic,
PlotLegends -> {"Exact", "Numeric"},
PlotLabel -> Row[{"Nn = ", Nn}]],
Plot[x^2*Pi^2/(h*(2*Nn + 1))^2, {x, 0, 30}]]
(*{2.45382, 9.81528, 22.0844, 39.2611, 61.3455, 88.3375, 120.237, \
157.044, 198.759, 245.382, 296.912, 353.35, 414.695, 480.949, \
552.109, 628.178, 709.154, 795.038, 885.829, 981.528, 1082.13, \
1187.65, 1298.07, 1413.4, 1533.64, 1658.78, 1788.84, 1923.8, 2063.66, \
2208.44}*)


Compare the numerical and exact eigenvalues of the Laplacian on the segment

Table[lambd[[i]]/(i^2*Pi^2/(2)^2), {i, 1, Length[lambd]}]

(*Out[]= {0.994495, 0.994495, 0.994495, 0.994495, 0.994495, 0.994495, \
0.994495, 0.994495, 0.994495, 0.994495, 0.994496, 0.994496, 0.994496, \
0.994496, 0.994496, 0.994496, 0.994496, 0.994496, 0.994496, 0.994496, \
0.994496, 0.994496, 0.994496, 0.994496, 0.994496, 0.994496, 0.994497, \
0.994497, 0.994497, 0.994497}*)


Find the length of the segment that corresponds to lambd[[i]]

dl = y /.
FindRoot[lambd[[1]]/(1^2*Pi^2/(h*(2*Nn + 1 + y))^2) == 1, {y, 1}]
(*0.760036*)


Consequently

L = h*(2*Nn + 1 + dl)

(*Out[]= 2.00553*)


Check that the eigenvalues correspond to L

Table[lambd[[i]]/(i^2*Pi^2/L^2), {i, 1, Length[lambd]}]

(*Out[]= {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., \
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}*)

• Thank you so much it helped me, but I can not use the standard method because i use sinc approximation to test its convergence on many models of the Schrodinger equation. Commented Jun 17, 2019 at 14:50
• @MohamedRaafat I added faster code. Commented Jun 17, 2019 at 15:36
• thank you so much it helped me figure the difference between the bounded and unbounded intervals for the sturm-liouville problem Commented Jun 19, 2019 at 12:15
• the only problem is that when i try to evaluate the laplacian equation (zero potential) the eigen values are not correct because they are not equal to the theoritical eigen values=n^2 pi^2/L^2 Commented Jun 19, 2019 at 12:43
• @MohamedRaafat It was possible to improve the algorithm, improve accuracy and speed. With this code you can solve the problem with vx[x_] =0 Commented Jun 19, 2019 at 14:50