5
$\begingroup$

I tried to solve coupled nonlinear differential equations from this paper https://sci-hub.hkvisa.net/10.1016/j.jcp.2019.109058 to see eigenvalues and eigenfunctions in different dimensions. At first, I tried to solve it in 1D following the algorithm mentioned on pages 11 and 12 of the paper in the matrix method. But it did not work. Then I tried NDEigensystem method.It did not work as well. Here is the problem statement enter image description here enter image description here

I like to generate the figures of u and v in different dimensions as enter image description here

Another algorithm to solve this kind of equation is given in this paper https://arxiv.org/pdf/1806.01244.pdf but in FORTRAN.

(*program in matrix method*)
\[Phi][n_, x_] := 
  1/Sqrt[2^2 n!] \[Pi]^(-1/4) Exp[-x^2/2] HermiteH[n, x];

b = 10; a = -10; m = 5;
h = (b - a)/m; xgrid = Table[a + j h, {j, 1, m - 1}];

v[x_] := 0.5 \[Omega] x^2;
\[Beta] = 1; \[Omega] = 0.1; \[Mu] = 0.1;
d2dx[n_] := 
  SparseArray[{Band[{1, 1}] -> -2, Band[{1, 2}] -> 1, 
     Band[{2, 1}] -> 1}, {n, n}]/(1/h^2);
nonlin = SparseArray[{i_, i_} :> 
    2 \[Beta] Abs[\[Phi][1, xgrid[[i]]]]^2 - \[Mu], 
   Length[xgrid] {1, 1}];
pot = SparseArray[{i_, i_} :> v[xgrid[[i]]], Length[xgrid] {1, 1}];
Ah[n_] := d2dx[n] + nonlin + pot;
Bh = SparseArray[{i_, i_} :> - \[Beta] Abs[\[Phi][1, xgrid[[i]]]]^2, 
   Length[xgrid] {1, 1}];
BDG[n_] := MatrixForm[ArrayFlatten[{{Ah[n], Bh}, {-Bh, -Ah[n]}}]];
{vals, funcs} = Eigensystem[N[BDG[m - 1]], 4, Method -> "Arnoldi"]


(*NDEigensystem method*)
\[Omega] = 0.1; \[Beta] = 100;
{gvals, gfuns} = 
  NDEigensystem[-0.5 Laplacian[\[Phi][x], {x}] + 
    0.5*\[Omega]^2 x^2 \[Phi][x] + \[Beta]  Abs[\[Phi][x]]^2 \[Phi][
      x]  , \[Phi][x], {x, -10, 10}, 4, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {MaxCellMeasure -> 0.05}}}}];

\[Mu]g = gval[[1]]; \[Phi]g = gfuns[[1]];

{vals, funs} = 
 NDEigensystem[{ -0.5 Laplacian[vl[x], {x}] + 
    0.5 \[Omega]^2 x^2 vl[x] + 
    2 \[Beta] Abs[\[Phi]g[x]]^2 vl[x] - \[Mu]g vl[
      x] - \[Beta] \[Phi]g[x]^2 ul[x], -0.5 Laplacian[ul[x], {x}] + 
    0.5 \[Omega]^2 x^2  ul[x] + 
    2 \[Beta] Abs[\[Phi]g[x]]^2 ul[x] - \[Mu]g vl[
      x] - \[Beta] \[Phi]g[x]^2 ul[x], 
   DirichletCondition[ul[x] == 0, x == 10 || x == -10], 
   DirichletCondition[vl[x] == 0, x == 10 || x == -10]}, {ul[x], 
   vl[x]}, {x, -10, 10}, 4, 
  Method -> {"PDEDiscretization" -> "MethodOfLines"}]
Improve this question
$\endgroup$
7
  • $\begingroup$ I am afraid that complex problems of this kind are ill-suited for this forum. Please run single lines at a time, read every message that shows up and try to understand what the problem is. If you cannot solve it, find a maximally simple example that gives that error, and then that may give a question suitable for this forum. Main problem with your {gvals, gfuns} = NDEigensystem[... is that the expression is not linear in \[Phi]. Also, stick to default methods unless you know what you are doing. $\endgroup$
    – user293787
    Nov 19, 2022 at 7:22
  • $\begingroup$ I did notice something strange, a simple example is f[_,_]:=Abort[]; NDEigensystem[f''[x],f[x],{x,-1,1},1] which leads to an abort in V12.3, which I find surprising (anyone?). For this reason, I suggest you drop the \[Phi][n_,x_]:=... definition in your code, do not forget to clear old definitions or restart your kernel, and perhaps use another name for it that does not interfere with the function \[Phi] in your NDEigensystem. $\endgroup$
    – user293787
    Nov 19, 2022 at 7:22
  • $\begingroup$ There is no chance to use NDEigensystem since problem is nonlinear any way due to constraint $\int {(u^2-v^2)dx=1}$. I think that we can try to solve nonlinear problem with using NMinimize to optimize energy. $\endgroup$
    – Alex Trounev
    Nov 24, 2022 at 5:54
  • $\begingroup$ If not with the NDEigensystem @AlexTrounev then can it be solved through the Matrix method as mentioned in the papers sci-hub.hkvisa.net/10.1016/j.jcp.2019.109058 (pages 11 and 12) and arxiv.org/pdf/1806.01244.pdf in Mathematica? $\endgroup$
    – Argha Debnath
    Nov 26, 2022 at 8:23
  • $\begingroup$ I don't understand methods described in the paper "Numerical methods for Bogoliubov-de Gennes excitations of Bose-Einstein condensates". Why they used normalization $\int {(u^2-v^2)}dx=1$, while there is a ground state with the same normalization and $u, v$ are small excitations over the ground state? Please, pay attention that in a linear case for $\beta=0, \omega_x=1$ they compute energy for oscillator as $\omega_i = 1, 2, 3, 4, ...$ (see Example 1, Case 1), but it should be $\omega_i = 1/2, 3/2, 5/2, 7/2, ... $. $\endgroup$
    – Alex Trounev
    Nov 26, 2022 at 14:31

1 Answer 1

Reset to default
4
$\begingroup$

As initial step to compute ground state we can use LDG method as follows

Clear["Global`*"]

Get["NumericalDifferentialEquationAnalysis`"];

L = 10; \[Beta] = 100; \[Omega] = 1; muTF = 
 1/2 (3 \[Omega] \[Beta]/2)^(2/3);
f = {2 (- mu p[t] + 1/2*\[Omega]^2 t^2 p[t] + \[Beta] p[t]^2 p[t]), 
  v[t]}; ini = {0, b1}; bc1 = 0;
LDGODEs[M0_, nn_, ns_, np_, f_, ini_, tmax_, tmin_] := 
 Module[{dx = (tmax - tmin)/nn,   A = Array[a, {M0 + 1, nn, 2}]}, 
  xl = Table[tmin + l*dx, {l, 0, nn}]; UT[m_, t_] := EulerE[m, t];
  psi[m_, k_, t_] := 
   Piecewise[{{UT[
       m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])], 
      xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
  g = Table[
    GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
  dp = Table[
    D[UT[m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])], 
     t], {k, nn}];
  rul = {p[t] -> 
     Sum[a[m, k, 2] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, nn}], 
    v[t] -> Sum[
      a[m, k, 1] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, nn}]};
  eq = Flatten[
    Table[-Sum[
         a[i + 1, n, 
           ks] (Table[(psi[i, n, t] If[j == 0, 0, 
                dp[[n]] /. m -> j]), {t, g[[n]][[All, 1]]}] . 
            g[[n]][[All, 2]]), {i, 0, 
          M0}] - (Table[(f[[ks]] /. rul) psi[j, n, t], {t, 
           g[[n]][[All, 1]]}] . g[[n]][[All, 2]]), {j, 0, M0}, {n, 1, 
       nn}, {ks, 1, 2}] + 
     Table[Sum[
        a[i + 1, n, 
          ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0, 
         M0}] - Sum[
        If[n < 2, ini[[ks]]/(M0 + 1), 
          a[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n, 
          xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks, 
       2 ns}]];
  eq1 = (p[t] - bc1) /. rul /. t -> tmax;
  eqn = Join[eq, {eq1}];
  var = Join[Flatten[A], {b1, mu}]; 
  energy = 
   Sum[Table[((1/2 v[t]^2 + 
          1/2 \[Omega]^2 t^2 p[t]^2 + \[Beta]/2 p[t]^4) /. rul), {t, 
       g[[n]][[All, 1]]}] . g[[n]][[All, 2]], {n, 1, nn}]; 
  norm = Sum[
     Table[(( p[t]^2) /. rul), {t, g[[n]][[All, 1]]}] . 
      g[[n]][[All, 2]], {n, 1, nn}] == 1/2; 
  con = Join[{norm}, Table[eqn[[i]] == 0, {i, Length[eqn]}]];
  sol1 = NMinimize[{energy, con}, var]; sol1]

Solution

M0 = 2; nn = 10; ns = 1; np = 5; tmax = L; tmin = 0; 
 ldgsol = 
  LDGODEs[M0, nn, ns, np, f, ini, tmax, tmin]; // AbsoluteTiming

Please note, that we minimize energy to compute $\mu _g$, Finally we have

mu /. sol1[[2]]

Out[]= 14.1341

This value we can compare to theoretical value muTF = 14.1155. Wave function we should extend on $-L\le x\le L$ as

lst = Table[{t, 
    Evaluate[
     Sum[a[i + 1, k, 2] psi[i, k, Abs[t]], {i, 0, M0}, {k, nn}] /. 
      sol1[[2]]]}, {t, -L, L, .0101}];

\[Phi]g = Interpolation[Join[lst, {{10, 0}}], InterpolationOrder -> 4];

Compare this function with theoretical one (red line)

Show[Plot[\[Phi]g[x], {x, -L, L}], 
 Plot[If[muTF > x^2/2, Sqrt[(muTF - x^2/2)/100], 0], {x, -10, 10}, 
  PlotStyle -> {Red, Dashed}]]

Figure 1

Using LDG method we can compute several excitations, but problem is correct on $-\infty <x<\infty$ only, while on limited interval we have some mixture of states. For example

 ns = 2;μg=mu/.sol1[[2]]; x1[t_] = Table[Symbol["x1" <> ToString[i]][t], {i, 1, ns}]; 
v1[t_] = 
 Table[Symbol["v1" <> ToString[i]][t], {i, 1, ns}]; bc = {0, 0}; f1 = 
 Join[{2 (0.5 \[Omega]^2 x^2 vl[x] + 
        2 \[Beta] Abs[\[Phi]g[x]]^2 vl[x] - \[Mu]g vl[
          x] - \[Beta] \[Phi]g[x]^2 ul[x] + mu1 vl[x]), 
     2 (0.5 \[Omega]^2 x^2 ul[x] + 
        2 \[Beta] Abs[\[Phi]g[x]]^2 ul[x] - \[Mu]g ul[
          x] - \[Beta] \[Phi]g[x]^2 vl[x] - mu1 ul[x])} /. {ul[x] -> 
      x11[t], vl[x] -> x12[t], \[Phi]g[x] -> 
      If[Abs[t] <= L, \[Phi]g[t], 0]} /. x -> t, v1[t]]; 
LDGODE[M0_, nn_, ns_, np_, f_, ini_, tmax_, tmin_] := 
 Module[{dx = (tmax - tmin)/nn,   A = Array[a, {M0 + 1, nn, 2 ns}]}, 
  xl = Table[tmin + l*dx, {l, 0, nn}]; 
  UT[m_, t_] := BernoulliB[m, t];
  psi[m_, k_, t_] := 
   Piecewise[{{UT[
       m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])], 
      xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
  g = Table[
    GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
  dp = Table[
    D[UT[m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])], 
     t], {k, nn}];
  rul1 = 
   Join[Table[
     x1[t][[i]] -> 
      Sum[a[m, k, i + ns] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, 
        nn}], {i, ns}], 
    Table[v1[t][[i]] -> 
      Sum[a[m, k, i] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, 
        nn}], {i, 1, ns}]];
  eq = Flatten[
    Table[-Sum[
         a[i + 1, n, 
           ks] (Table[(psi[i, n, t] If[j == 0, 0, 
                dp[[n]] /. m -> j]), {t, g[[n]][[All, 1]]}] . 
            g[[n]][[All, 2]]), {i, 0, 
          M0}] - (Table[(f[[ks]] /. rul1) psi[j, n, t], {t, 
           g[[n]][[All, 1]]}] . g[[n]][[All, 2]]), {j, 0, M0}, {n, 1, 
       nn}, {ks, 2 ns}] + 
     Table[Sum[
        a[i + 1, n, 
          ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0, 
         M0}] - Sum[
        If[n < 2, ini[[ks]]/(M0 + 1), 
          a[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n, 
          xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks, 
       2 ns}]];
  eqb = {(x1[t][[1]] - bc[[1]]), x1[t][[2]] - bc[[2]]} /. rul1 /. 
    t -> tmax;
  
  var = Join[Flatten[A], {b1, b2, mu1}]; 
  norm1 = {Sum[
      Table[(( x11[t]^2 - x12[t]^2) /. rul1), {t, 
         g[[n]][[All, 1]]}] . g[[n]][[All, 2]], {n, 1, nn}] == 1/2}; 
  eqn = Join[eq, eqb];
  sol2 = 
   FindRoot[Join[Table[eqn[[i]] == 0, {i, Length[eqn]}], norm1], 
    Table[{var[[i]], 1/10}, {i, Length[var]}]]; sol2];

Solution

ini = Join[{b1, b2}, {0, 0}]; M0 = 4; nn = 40; ns = 2; np = 5; 
 tmax = 9.9863; tmin = 0; 
 ldgsol2 = 
  LDGODE[M0, nn, ns, np, f1, ini, tmax, tmin]; // AbsoluteTiming

It should be first eigenvalue $\omega_1=1$, but we have instead

mu1 /. sol2

Out[]= 1.00091

Nevertheless, we have norm 2 norm1[[1, 1]] /. sol2 =1. Theoretical solution for eigenfunction is given by

u0[x_] := 1/Sqrt[2] (x \[Phi]g[x] - \[Phi]g'[x]); 
v0[x_] := 1/Sqrt[2] (x \[Phi]g[x] + \[Phi]g'[x])

Numerical solution is

pl1 = Table[
  Plot[Evaluate[
    Sum[a[i + 1, k, s] psi[i, k, t], {i, 0, M0}, {k, nn}] /. 
     sol2], {t, tmin, tmax}, PlotRange -> All, PlotStyle -> {Red}, 
   Exclusions -> None], {s, ns + 1, 2 ns}]

Compare to theory they look like

pl = {Plot[u0[x], {x, 0, Sqrt[2 muTF] - .07}], 
  Plot[v0[x], {x, 0, Sqrt[2 muTF] - .07}]}

Show[pl1, pl]

Figure 2

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks, @AlexTrounev for this solution. When I ran this program. In the first part where you minimize energy to compute μg I am getting -0.163079 instead of 14.1341 and consequently, the resultant \phi_g profile is also different than posted one. In the second program, μg is not defined in f1. $\endgroup$
    – Argha Debnath
    Nov 29, 2022 at 6:47
  • $\begingroup$ I have many questions regarding the program @AlexTrounev. How f is defined mainly second part (v[t]) in {}. How the definition of \psi is working? Can you provide in short the algorithm in the eq part? $\endgroup$
    – Argha Debnath
    Nov 29, 2022 at 6:53
  • $\begingroup$ How can I get @AlexTrounev Other Eigenvalues and EigenFunctions? $\endgroup$
    – Argha Debnath
    Nov 29, 2022 at 6:54
  • $\begingroup$ @ArghaDebnath 1) I run first part of code and evaluate mu /. sol1[[2]] = 14.1341. Plot also looks same as above. Could you run again with a fresh kernel? Please, don't run together first and second part. Run first part, compute $\mu g$ as μg=mu/.sol1[[2]], and then run second part of code. 2) we solve first order ODEs. 3) With LDG method we can compute several modes. $\endgroup$
    – Alex Trounev
    Nov 29, 2022 at 9:05

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.