# Solving Coupled nonlinear Differential Equations For Eigenvalue And Eigen functions

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

I like to generate the figures of u and v in different dimensions as

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"}]
• 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. Commented Nov 19, 2022 at 7:22
• 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. Commented Nov 19, 2022 at 7:22
• 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. Commented Nov 24, 2022 at 5:54
• 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? Commented Nov 26, 2022 at 8:23
• 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, ...$. Commented Nov 26, 2022 at 14:31

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

Using LDG method we can compute several excitations, but problem is correct on $$-\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]

• 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. Commented Nov 29, 2022 at 6:47
• 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? Commented Nov 29, 2022 at 6:53
• How can I get @AlexTrounev Other Eigenvalues and EigenFunctions? Commented Nov 29, 2022 at 6:54
• @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. Commented Nov 29, 2022 at 9:05