2
$\begingroup$

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]
Improve this question
$\endgroup$

1 Answer 1

Reset to default
11
$\begingroup$

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

Figure 1

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

Figure 2

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

Figure 3

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}*)

fig4

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}*)

figure 5

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.}*)
Improve this answer
$\endgroup$
11
  • $\begingroup$ 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. $\endgroup$
    – Mohamed Raafat
    Commented Jun 17, 2019 at 14:50
  • $\begingroup$ @MohamedRaafat I added faster code. $\endgroup$
    – Alex Trounev
    Commented Jun 17, 2019 at 15:36
  • $\begingroup$ thank you so much it helped me figure the difference between the bounded and unbounded intervals for the sturm-liouville problem $\endgroup$
    – Mohamed Raafat
    Commented Jun 19, 2019 at 12:15
  • $\begingroup$ 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 $\endgroup$
    – Mohamed Raafat
    Commented Jun 19, 2019 at 12:43
  • $\begingroup$ @MohamedRaafat It was possible to improve the algorithm, improve accuracy and speed. With this code you can solve the problem with vx[x_] =0 $\endgroup$
    – Alex Trounev
    Commented Jun 19, 2019 at 14:50

Your Answer

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

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