I am trying to solve Equation number (1.2) numerically in MATHEMATICA. This equation is solved in the papers https://arxiv.org/pdf/2205.05193.pdf, https://arxiv.org/pdf/2202.13264.pdf, and https://arxiv.org/pdf/2005.05718.pdf. This equation is quite similar to I am trying to solve nonlinear Schrödinger equation with dipolar interaction which @AlexTrounev solved. I like to find the figures posted below and energy especially kinetic and dipolar energy. This paper might be a little help https://arxiv.org/pdf/1506.03283.pdf as Equation number (1.3) goes to infinity at the center.
Here I paste my failed attempt to modify the code from @AlexTrounev
Needs["NumericalDifferentialEquationAnalysis`"];
np = 4; xg = GaussianQuadratureWeights[np, -L, L]; xpoints =
xg[[All, 1]];
xweights = xg[[All, 2]];
yg = GaussianQuadratureWeights[np, -L, L]; ypoints = yg[[All, 1]];
yweights = yg[[All, 2]];
zg = GaussianQuadratureWeights[np, -L, L]; zpoints = zg[[All, 1]];
zweights = zg[[All, 2]];
sol[0][x_, y_, z_, t_] :=
E^(- (x^2 + y^2 + z^2)/2);
\[Omega]x = 45; \[Omega]y = 45; \[Omega]z = 133; gmean\[Omega] = (\
\[Omega]x \[Omega]y \[Omega]z)^(
1/3); \[Gamma] = \[Omega]x/gmean\[Omega]; \[Lambda] = \
\[Omega]y/gmean\[Omega]; \[Nu] = \[Omega]z/gmean\[Omega];
a0 = 5.29*10^-5; a = 85 a0; add = 131 a0; \[Epsilon]dd = add/a;
\[Gamma]QF =
128/3 Sqrt[\[Pi] a^5] (1 + 1.5 \[Epsilon]dd^2); angle = ((
3 Cos[\[Phi]]^2 - 1)/2); \[Phi] = 0;
Nat = 10^6; L = 10; dt = 1/100; nt = 100; T =
Table[i dt, {i, 0, 1001}];
nkx[0] = Table[
Sum[xweights[[i]] Exp[
I xpoints[[j]] xpoints[[i]]] Abs[
sol[0][xpoints[[i]], 0, 0, 0]]^2, {i, Length[xg]}], {j,
Length[xg]}];
nky[0] = Table[
Sum[yweights[[i]] Exp[
I ypoints[[j]] ypoints[[i]]] Abs[
sol[0][0, ypoints[[i]], 0, 0]]^2, {i, Length[yg]}], {j,
Length[yg]}];
nkz[0] = Table[
Sum[zweights[[i]] Exp[
I zpoints[[j]] zpoints[[i]]] Abs[
sol[0][0, 0, zpoints[[i]], 0]]^2, {i, Length[zg]}], {j,
Length[zg]}];
finalnk[0] = nkx[0] nky[0] nkz[0];
intn[0] =
1/(6 \[Pi]^2)
angle Table[{xpoints[[j]], ypoints[[k]], zpoints[[m]],
Re[Sum[xweights[[j]] yweights[[k]] zweights[[
m]] Exp[-I xpoints[[j]] xpoints[[i]]] Exp[-I ypoints[[
k]] ypoints[[i]]] Exp[-I zpoints[[m]] zpoints[[i]]] ((
3 zpoints[[m]]^2)/(
xpoints[[j]]^2 + ypoints[[k]]^2 + zpoints[[m]]^2) -
1) finalnk[0][[i]], {i, Length[xg]}]]}, {j, Length[xg]}, {k,
Length[yg]}, {m, Length[zg]}];
Vdd[0] = Interpolation[
Join[{{{-L, -L, -L, intn[0][[1, 1, 1, 4]]}}},
intn[0], {{{L, L, L, intn[0][[np, np, np, 4]]}}}]]
Do[sol[s] =
NDSolveValue[{-I D[\[Psi][x, y, z, t], t] -
1/2 Laplacian[\[Psi][x, y, z, t], {x, y, z}] +
1/2 (\[Gamma]^2 x^2 + \[Nu]^2 y^2 + \[Lambda]^2 z^2) \[Psi][x,
y, z, t] +
4 \[Pi] a Nat Abs[\[Psi][x, y, z, t]]^2 \[Psi][x, y, z, t] +
3 add Nat Vdd[s - 1][x, y, z] \[Psi][x, y, z,
t] + \[Gamma]QF Nat^(3/2)
Abs[\[Psi][x, y, z, t]]^3 \[Psi][x, y, z, t] == 0,
\[Psi][x, y, z, T[[s]]] ==
sol[s - 1][x, y, z, T[[s]]], \[Psi][L, y, z, t] ==
0, \[Psi][-L, y, z, t] == 0,
\[Psi][x, L, z, t] == 0, \[Psi][x, -L, z, t] ==
0, \[Psi][x, y, L, t] == 0, \[Psi][x, y, -L, t] ==
0}, \[Psi], {t, T[[s]], T[[s + 1]]}, {x, -L, L}, {y, -L,
L}, {z, -L, L},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 40, "MaxPoints" -> np,
"DifferenceOrder" -> "Pseudospectral"}}] // Quiet;
nkx[s] =
Table[Sum[
xweights[[i]] Exp[
I xpoints[[j]] xpoints[[i]]] Abs[
sol[s][xpoints[[i]], 0, 0, T[[s + 1]]]]^2, {i,
Length[xg]}], {j, Length[xg]}];
nky[s] =
Table[Sum[
yweights[[i]] Exp[
I ypoints[[j]] ypoints[[i]]] Abs[
sol[s][0, ypoints[[i]], 0, T[[s + 1]]]]^2, {i,
Length[yg]}], {j, Length[yg]}];
nkz[s] =
Table[Sum[
zweights[[i]] Exp[
I zpoints[[j]] zpoints[[i]]] Abs[
sol[s][0, 0, zpoints[[i]], T[[s + 1]]]]^2, {i,
Length[zg]}], {j, Length[zg]}];
finalnk[s] = nkx[s] nky[s] nkz[s];
intn[s] =
1/(6 \[Pi]^2)
angle Table[{xpoints[[j]], ypoints[[k]], zpoints[[m]],
Re[Sum[xweights[[j]] yweights[[k]] zweights[[
m]] Exp[-I xpoints[[j]] xpoints[[i]]] Exp[-I ypoints[[
k]] ypoints[[i]]] Exp[-I zpoints[[m]] zpoints[[i]]] ((
3 zpoints[[m]]^2)/(
xpoints[[j]]^2 + ypoints[[k]]^2 + zpoints[[m]]^2) -
1) finalnk[s][[i]], {i, Length[xg]}]]}, {j, Length[xg]}, {k,
Length[yg]}, {m, Length[zg]}];
Vdd[s] =
Interpolation[
Join[{{{-L, -L, -L, intn[s][[1, 1, 1, 4]]}}},
intn[s], {{{L, L, L, intn[s][[np, np, np, 4]]}}}]];, {s, 1, 15}]
DensityPlot[
Abs[sol[10][x, y, 0, T[[10 + 1]]]]^2, {x, -L, L}, {y, -L, L},
ColorFunction -> "BlueGreenYellow", PlotPoints -> 200, Frame -> True,
FrameTicks -> Automatic, LabelStyle -> {24, Bold, Large, Black},
FrameLabel -> {{Style["y", FontFamily -> "Times New Roman",
FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30],
None}, {Style["x", FontFamily -> "Times New Roman",
FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30],
None}}]
