2
$\begingroup$

I am trying to solve Eq.(3) from (https://www.sciencedirect.com/science/article/pii/S0378437113002161 or https://sci-hub.se/https://doi.org/10.1016/j.physa.2013.03.013 ). $$ i\frac{\partial\phi}{\partial t}=(-\frac{\partial^2}{\partial x^2}+\frac{x^2}{4}-\frac{|\phi|^2}{x^2}-2i\xi\frac{|\phi|^4}{x^4}+i\frac{\Gamma}{2})\phi $$ $\phi(0)=\phi(\infty)=0$

By switching off the nonconservative terms, a stable solution is obtained. By turning on the nonconservative terms, they saw the dynamical collapse. They have used $\Gamma=0.15;\Delta x=0.004;\Delta t=0.001;l=(0:40);\xi=0.001$. My goal is to obtain the phase plot in fig.(1) between $X(t)=\int_0^\infty{x^2|\phi(x,t)|^2dx}$ vs $\frac{dX(t)}{dt}$.

enter image description here

I solved the equation using imaginary time propagation method without the nonconservative terms.

Needs["NDSolve`FEM`"]; l = 40; finalt = 1; dt = 1/1000; \[CapitalGamma] = 0.15; \[Zeta] = 0.001; g = 1; pse[x_, t_] := Exp[-((x - 0.3)/0.05)^2]; pnumber = NIntegrate[Abs[pse[x, 0]]^2, {x, 0, l}] pot[x_] := x^2/4; Psi[0][x_] := pse[x, 0]; k[0] = 1; Do[Psi[i] = Quiet@NDSolveValue[{(\[Psi][x] - k[i - 1] Psi[i - 1][x])/dt == Laplacian[\[Psi][x], {x}] + g Abs[k[i - 1] Psi[i - 1][x]]^2 \[Psi][x]/x^2 - pot[x] \[Psi][x], DirichletCondition[\[Psi][x] == 0, True]}, \[Psi], {x, 0, l}, Method -> {"FiniteElement", InterpolationOrder -> {\[Psi] -> 2}, "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}]; num[i] = NIntegrate[Abs[Psi[i][x]]^2, {x, 0, l}]; k[i] = Sqrt[pnumber/num[i]], {i, 1, finalt/dt}] // AbsoluteTiming

The solution

Plot[Evaluate[ Table[Abs[k[1000] Psi[1000][x]]^2, {i, 100, 200, 20}]], {x, 0, 5}, PlotRange -> All]

The obtained solution

Then using the solution and switching on the nonconservative terms The dynamics can be seen in Real time

Psi1[0][x_] := Abs[k[1000] Psi[1000][x]]; Do[Psi1[i] = Quiet@NDSolveValue[{I (\[Phi][x] - Psi1[i - 1][x])/ dt == -Laplacian[\[Phi][x], {x}] - g Abs[Psi1[i - 1][x]]^2 \[Phi][x]/x^2 + pot[x] \[Phi][x] - 2 I \[Zeta] Abs[Psi1[i - 1][x]]^4 \[Phi][x]/x^4 + I \[CapitalGamma] \[Phi][x]/2, DirichletCondition[\[Phi][x] == 0, True]}, \[Phi], {x, 0, l}, Method -> {"FiniteElement", InterpolationOrder -> {\[Phi] -> 2}, "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}];, {i, 1, 1000}] // AbsoluteTiming

The dynamics

data = Flatten[ Evaluate@Table[{i, x, Abs[Psi1[i][x]]^2}, {i, 1, 2000}, {x, 0, 40}], 1]; ListDensityPlot[data, ColorFunction -> "SolarColors", FrameLabel -> Automatic, AspectRatio -> 1/2, PlotLegends -> Automatic]

enter image description here

The calculation for phase plot

Needs["NumericalDifferentialEquationAnalysis`"]; gg = GaussianQuadratureWeights[1000, 0, 40]; point = gg[[All, 1]]; weight = gg[[All, 2]]; rms = Table[{t, Sqrt[weight . Table[Evaluate[x^2 Abs[Psi1[t][x]]^2], {x, point}]]}, {t, 0, 1000}]; x\[Tau] = Interpolation[rms, InterpolationOrder -> 3]; dx\[Tau][t_] = D[x\[Tau][t], {t, 1}];

Variation of $X(t)$ and $\frac{dX}{dt}$ with $t$

Plot[{x\[Tau][t], dx\[Tau][t]}, {t, 0, 2000}, PlotRange -> All, PlotLegends -> "Expressions"]

enter image description here

The phase plot

 ParametricPlot[{x\[Tau][t], dx\[Tau][t]}, {t, 0, 2000}, PlotRange -> All]
Improve this question
$\endgroup$
16
  • 3
    $\begingroup$ There is no question to answer here. $\endgroup$ Commented Oct 21, 2024 at 3:23
  • $\begingroup$ @AlexTrounev I would like to obtain the phase plot attached above and know the process of obtaining the Lyapunov exponent I am not getting. Is my method of obtaining the solution without the nonconservative terms and seeing the dynamics switching on them correct? It is kind of the reverse of quench dynamics. $\endgroup$ Commented Oct 21, 2024 at 5:42
  • 1
    $\begingroup$ Can you upload the picture you got? $\endgroup$ Commented Oct 22, 2024 at 9:29
  • $\begingroup$ @AlexTrounev I could not upload the parametric phase plot as it's not at all giving anything as $dX/dt$ is coming very small which I uploaded here also the dynamics I am getting. $\endgroup$ Commented Oct 22, 2024 at 10:43
  • $\begingroup$ 1. Seems that you're using the method here for discretization in $t$ direction, but as pointed out therein, your original attempt fails probably because you don't code the equation correctly, you probably need to re-visit this part. 2. How is the initial condition decided? I've checked the paper and the cited papers, but this doesn't seem to be explicitly mentioned. $\endgroup$ Commented Oct 23, 2024 at 2:36

1 Answer 1

Reset to default
2
+50
$\begingroup$

From the data present in the paper in Figure 1, it follows that the number of particles is not conserved and can increase or decrease. While from the model itself this does not follow. For any initial data, the number of particles decreases in this model. On the other hand, the model with particles number stabilizations shows chaotic behaver

Needs["NDSolve`FEM`"]; l = 40; finalt = 30; dt = 1/100; \[CapitalGamma] = 0.15; \[Zeta] = 0.001; g = 1; pse[x_, t_] := 7 x^2 Exp[-((x - 1))^2]; pnumber = NIntegrate[Abs[pse[x, 0]]^2, {x, 0, l}]; pot[x_] := x^2/4; Psi1[0][x_] := pse[x, 0]; k[0] = 1; Do[Psi1[i] = Quiet@NDSolveValue[{I (\[Phi][x] - k[i - 1] Psi1[i - 1][x])/ dt == -Laplacian[\[Phi][x], {x}] - g Abs[k[i - 1] Psi1[i - 1][x]]^2 \[Phi][x]/x^2 + pot[x] \[Phi][x] - 2 I \[Zeta] Abs[k[i - 1] Psi1[i - 1][x]]^4 \[Phi][x]/ x^4 + I \[CapitalGamma] \[Phi][x]/2, DirichletCondition[\[Phi][x] == 0, True]}, \[Phi], {x, 0, l}, Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}]; numb[i] = NIntegrate[Abs[Psi1[i][x]]^2, {x, 0, l}, AccuracyGoal -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"]; k[i] = Sqrt[pnumber/numb[i]];, {i, 1, 3000}] // AbsoluteTiming X = Sqrt@ Table[NIntegrate[x^2 Abs[Psi1[i][x]]^2, {x, 0, l}, AccuracyGoal -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"], {i, 3000}]; dX = Differences[X]/dt; ListPlot[Transpose[{Drop[X, 1], dX}], PlotRange -> All, ColorFunction -> Hue, Background -> Black, Frame -> True]

Figure 1

Update 1. Note that for pse[x_, t_] := .544915 x^2 Exp[-((x - 1))^2] we have pnumber=1, and the corresponding motion looks very regular Figure 2

But even for pnumber=3 we have chaotic patterns Figure 3

Improve this answer
$\endgroup$
8
  • $\begingroup$ Could you give me some insights on how their phase plot can be exactly reproduced.@AlexTrounev $\endgroup$ Commented Oct 25, 2024 at 12:21
  • $\begingroup$ Thanks for the answare@AlexTrounev. Can I take the obtained solution from imaginary time method without the non-conservative terms and use it as initial solution for real time propagation with non-conservative terms as I posted above? Is there any another process of obtaining Lyapunov exponent other than ( Eq.[7] fig[2]) from the paper. $\endgroup$ Commented Oct 25, 2024 at 12:23
  • $\begingroup$ In this paper @AlexTrounev sci-hub.se/https://doi.org/10.1016/j.chaos.2021.111193 they saw the same dynamics. Do you think to generate fig 1 I need to include particle stabilizer k[I] or without it will suffice $\endgroup$ Commented Oct 25, 2024 at 12:23
  • $\begingroup$ 1. For my opinion data in Figure 1 cannot be reproduced with implicit step used in your code and above. 2. We don't need preprocessing to get chaotic behaver. Hence, your code with imaginary time is useless. Nevertheless, we need some initial data to reproduce Figure 1. But we don't know how function $\Phi (x,0)$ used in the paper looks like. It is clear that pnumber>>1, and in our example we take pnumber=165. But it is unclear what pnumber they took to get data in Figure 1. $\endgroup$ Commented Oct 25, 2024 at 14:50
  • $\begingroup$ Thanks, @AlexTrounev for clarification. $\endgroup$ Commented Oct 25, 2024 at 17:31

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.