# Poincare section for a perturbed string

I am trying to find the Poincare section of a perturbed string in a gravitational background with a magnetic field. The plots that I am trying to reproduce are from this paper (page 11).

The code that I have used for this task is explained in detail below:

a) Firstly, I enter the values of different coefficients which are involved in the equation and input the lagrangian.

K1 = 6.90; K2 = 16.05; K3 = 9.65; K4 = 3.27; K5 = 6.55;
\[Omega]sq = -0.923; \[Omega]sq = 6.478;
lagrangian =
Sum[c[n]'[t]^2 - c[n][t]^2 \[Omega]sq[n], {n, {0, 1}}] +
K1 c[t]^3 + K2 c[t] c[t]^2 + K3 c[t] c'[t]^2 +
K4 c[t] c'[t]^2 + K5 c'[t] c[t] c'[t];


b) Since in some regions of the potential the kinetic term is negative the coefficients are replaced as under:

c[t_] := OverTilde[c][t] + \[Alpha]1*OverTilde[c][t]^2 + \[Alpha]2*OverTilde[c][t]^2;
c[t_] := OverTilde[c][t] + \[Alpha]3*OverTilde[c][t]*OverTilde[c][t];
\[Alpha]1 = -2; \[Alpha]2 = -0.5; \[Alpha]3 = -1;


c) I find the modified lagrangian and hamiltonian from the above equations:

n = Expand[lagrangian];
vars = {OverTilde[c][t], OverTilde[c][t],
Derivative[OverTilde[c]][t],
Derivative[OverTilde[c]][t]};
lagrangian =
Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1;
momentum[n_] := D[lagrangian, Derivative[OverTilde[c][n]][t]]
hamiltonian = Expand[Sum[momentum[n]*Derivative[OverTilde[c][n]][t], {n, {0, 1}}] - lagrangian];


d) I solve the modified lagrangian to find the equations of motion:


eulerLagrange[lagrangian_, vars_, dvars_] :=
Thread[Table[D[D[lagrangian, dvar], t], {dvar, dvars}] - Table[D[lagrangian, var], {var, vars}] ==
ConstantArray[0, Length[vars]]];
equationsOfMotion = eulerLagrange[lagrangian, {OverTilde[c][t], OverTilde[c][t]},
{Derivative[OverTilde[c]][t], Derivative[OverTilde[c]][t]}];


e) I plot the Poincare section using NDSolve and WhenEvent technique using the "StiffnessSwitching" method:

sol = Table[Block[{a, b, \[Chi], d}, {a, b, \[Chi], d} = {x, 0, 0.01, 0.0001};
Reap[NDSolve[{equationsOfMotion, OverTilde[c] == a, Derivative[OverTilde[c]] == b,
OverTilde[c] == \[Chi], Derivative[OverTilde[c]] == d,
WhenEvent[OverTilde[c][t] == 0 && Derivative[OverTilde[c]][t] >= 0,
Sow[{OverTilde[c][t], Derivative[OverTilde[c]][t]}]]},
{OverTilde[c][t], OverTilde[c][t]}, {t, 0, 8000}, Method -> "StiffnessSwitching"]]],
{x, -0.2, 0., 0.01}];
ListPlot[Table[Flatten[sol[[i]][], 1], {i, Length[sol]}], PlotTheme -> "Detailed",
FrameLabel -> {OverTilde[Subscript[c, 0]], OverDot[OverTilde[Subscript[c, 0]], 1]}]


f) I try to zoom in on the above plot in the desired range to observe chaos :

ListPlot[Table[Flatten[sol[[i]][], 1], {i, Length[sol]}], PlotTheme -> "Detailed",
PlotRange -> {{-0.1, 0.}, {-0.05, 0.05}}, FrameLabel -> {OverTilde[Subscript[c, 0]],
OverDot[OverTilde[Subscript[c, 0]], 1]}]


After highlighting all the steps, here are the problems which I face with the above code:

1. I get an NDSolve error: "step size is effectively zero; singularity or stiff system suspected."
2. I am unable to reproduce the chaotic aspect of the Poincare plots.
3. I am not sure how to change the initial conditions with fixed Energy as mentioned in the paper. I tried to find the hamiltonian but don't understand how to use it effectively in the code.
4. The effect of a Magnetic field is to mitigate the chaotic behavior which can be done by changing the values of K variables as given in the paper, I am unable to see any effect in that regard.
5. The plot range, mainly the highlighted section differs from that given in the paper.

I have been stuck with this problem for weeks now, any help in reproducing the plots would be truly beneficial!

Acknowledgments: I would like to thank @Diffycue for his generosity in helping me to start with the code and @bmf to help me with the code in order to neglect the higher-order coefficients.

Edit 1 Thanks to @Alex Trounev 's answer, I tried reproducing the plot of a similar paper (Figure 4, page 5) and got this result: The code that I used to get the above plot taking help from @Alex Trounev is given below:

\[Omega]sq = -1.4; \[Omega]sq = 7.57;
lagrangian =
Sum[c[n]'[t]^2 - c[n][t]^2 \[Omega]sq[n], {n, {0, 1}}] +
7.11 c[t]^3 + 35.3 c[t] c[t]^2 +
4.66 c[t] c'[t]^2 + 1.32 c[t] c'[t]^2 -
7.57 c'[t] c[t] c'[t];
momentum[n_] := D[lagrangian, c[n]'[t]]
hamiltonian =
Sum[momentum[n] c[n]'[t], {n, {0, 1}}] - lagrangian // Expand;

eulerLagrange[lagrangian_, vars_, dvars_] :=
Thread[(Table[D[D[lagrangian, dvar], t], {dvar, dvars}] -
Table[D[lagrangian, var], {var, vars}]) ==
ConstantArray[0, Length@vars]];
equationsOfMotion =
eulerLagrange[lagrangian, {c[t], c[t]}, {c'[t], c'[t]}];
sol = Table[
Block[{a, b, \[Chi], d, habd},
habd = hamiltonian /. t -> 0 /. {c -> a, c' -> b,
c -> \[Chi], c' -> d}; {a, b, d} = {x, 10^-2,
10^-4}; \[Chi] = \[Chi] /.
FindRoot[habd == 9.28 10^-6, {\[Chi], -x}];
Reap[NDSolve[{equationsOfMotion, c == a, c' == b,
c == \[Chi], c' == d,
WhenEvent[c[t] == 0 && c'[t] >= 0,
Sow[{c[t], c'[t]}]]}, {c[t], c[t]}, {t, 0, 8000},
Method -> "StiffnessSwitching"]]], {x, -0.2, -0.1, 0.001}];
ListPlot[Table[Flatten[sol[[i]][], 1], {i, Length@sol}],
PlotTheme -> "Scientific",
PlotRange -> {{-0.05, 0.00}, {-0.05, 0.05}}]


There is still room for improvement. I was not able to get the proper Plot Range for the above plot as in the paper. Also, the scattered points indicating chaos differ from that given in the paper.

• Please, pay attention for action defined by (B12) and (B14) from Appendix B of Chaos of Wilson Loop from String Motion near Black Hole Horizon. It seems that you use (B14) and then normalize variables with using (B13), while they used (B12) before normalization. May 4 at 10:08
• @AlexTrounev, actually I have directly used the normalized variables from (B14) and replaced OverTilde[c] with regular [c]. It might be the source of confusion. Still, I will again verify the same. May 4 at 13:26
• It looks like they used Mathematica as well. May 4 at 16:34
• Yes, they have definitely used Mathematica. I have played with different initial conditions and now I think it all boils down with choosing right set of initial conditions to start with. May 4 at 16:47

Small modification leads to appearance of random behaver and solves step size problem

par = {6.90, 16.05, 9.65, 3.27, 6.55, -0.923, 6.478}; {K1, K2, K3, K4,
K5, \[Omega]sq, \[Omega]sq} = Rationalize[par, 10^-30];
lagrangian =
Sum[c[n]'[t]^2 - c[n][t]^2 \[Omega]sq[n], {n, {0, 1}}] +
K1 c[t]^3 + K2 c[t] c[t]^2 + K3 c[t] c'[t]^2 +
K4 c[t] c'[t]^2 + K5 c'[t] c[t] c'[t];

c[t_] :=
OverTilde[c][t] + \[Alpha]1*OverTilde[c][t]^2 + \[Alpha]2*
OverTilde[c][t]^2;
c[t_] :=
OverTilde[c][t] + \[Alpha]3*OverTilde[c][t]*OverTilde[c][t];
\[Alpha]1 = -2; \[Alpha]2 = -1/2; \[Alpha]3 = -1;

n = Expand[lagrangian];
vars = {OverTilde[c][t], OverTilde[c][t],
Derivative[OverTilde[c]][t],
Derivative[OverTilde[c]][t]};
lagrangian =
Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1;
momentum[n_] := D[lagrangian, Derivative[OverTilde[c][n]][t]]
hamiltonian =
Expand[Sum[
momentum[n]*Derivative[OverTilde[c][n]][t], {n, {0, 1}}] -
lagrangian];
(* I solve the modified lagrangian to find the equations of motion:*)
eulerLagrange[lagrangian_, vars_, dvars_] :=
Thread[Table[D[D[lagrangian, dvar], t], {dvar, dvars}] -
Table[D[lagrangian, var], {var, vars}] ==
ConstantArray[0, Length[vars]]];
equationsOfMotion =
eulerLagrange[
lagrangian, {OverTilde[c][t],
OverTilde[c][t]}, {Derivative[OverTilde[c]][t],
Derivative[OverTilde[c]][t]}];

sol = Table[Block[{a, b, \[Chi], d}, {a, b, \[Chi], d} = {x, 0, -x, 0};
Reap[NDSolve[{equationsOfMotion, OverTilde[c] == a,
Derivative[OverTilde[c]] == b,
OverTilde[c] == \[Chi],
Derivative[OverTilde[c]] == d,
WhenEvent[
OverTilde[c][t] == 0 &&
Derivative[OverTilde[c]][t] >= 0,
Sow[{OverTilde[c][t],
Derivative[OverTilde[c]][t]}]]}, {OverTilde[c][t],
OverTilde[c][t]}, {t, 0, 8000},
Method ->
"StiffnessSwitching"]]], {x, {-1.5 10^-3, -10^-3, -5 10^-4,
-10^-4, -5 10^-5}}];

ListPlot[Table[Flatten[sol[[i]][], 1], {i, Length[sol]}],
PlotTheme -> "Detailed",
PlotRange -> {{-2 10^-3, 0}, {-10^-3, 10^-3}},
FrameLabel -> {OverTilde[Subscript[c, 0]],
OverDot[OverTilde[Subscript[c, 0]]]},
PlotStyle -> PointSize[.0025]] Update 1. We can fixed Hamiltonian, for instance at hamiltonian=10^-5, in computation as follows (example for x1 with B=0.6)

par = {9.57, 19.56, 10.24, 3.33, 6.67, -1.202, 7.213}; {K1, K2, K3,
K4, K5, \[Omega]sq, \[Omega]sq} = Rationalize[par, 10^-30];
lagrangian =
Sum[c[n]'[t]^2 - c[n][t]^2 \[Omega]sq[n], {n, {0, 1}}] +
K1 c[t]^3 + K2 c[t] c[t]^2 + K3 c[t] c'[t]^2 +
K4 c[t] c'[t]^2 + K5 c'[t] c[t] c'[t];

c[t_] :=
OverTilde[c][t] + \[Alpha]1*OverTilde[c][t]^2 + \[Alpha]2*
OverTilde[c][t]^2;
c[t_] :=
OverTilde[c][t] + \[Alpha]3*OverTilde[c][t]*OverTilde[c][t];
\[Alpha]1 = -2; \[Alpha]2 = -1/2; \[Alpha]3 = -1;

n = Expand[lagrangian];
vars = {OverTilde[c][t], OverTilde[c][t],
Derivative[OverTilde[c]][t],
Derivative[OverTilde[c]][t]};
lagrangian =
Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1;
momentum[n_] := D[lagrangian, Derivative[OverTilde[c][n]][t]]
hamiltonian =
Expand[Sum[
momentum[n]*Derivative[OverTilde[c][n]][t], {n, {0, 1}}] -
lagrangian];
(* I solve the modified lagrangian to find the equations of motion:*)
eulerLagrange[lagrangian_, vars_, dvars_] :=
Thread[Table[D[D[lagrangian, dvar], t], {dvar, dvars}] -
Table[D[lagrangian, var], {var, vars}] ==
ConstantArray[0, Length[vars]]];
equationsOfMotion =
eulerLagrange[
lagrangian, {OverTilde[c][t],
OverTilde[c][t]}, {Derivative[OverTilde[c]][t],
Derivative[OverTilde[c]][t]}];

sol = Table[
Block[{a, b, \[Chi], d, habd},
habd = hamiltonian /. t -> 0 /. {OverTilde[c] -> a,
Derivative[OverTilde[c]] -> b,
OverTilde[c] -> \[Chi],
Derivative[OverTilde[c]] -> d}; {a, b, d} = {x, 0,
0}; \[Chi] = \[Chi] /. FindRoot[habd == 10^-5, {\[Chi], -x}];
Reap[NDSolve[{equationsOfMotion, OverTilde[c] == a,
Derivative[OverTilde[c]] == b,
OverTilde[c] == \[Chi],
Derivative[OverTilde[c]] == d,
WhenEvent[
OverTilde[c][t] == 0 &&
Derivative[OverTilde[c]][t] >= 0,
Sow[{OverTilde[c][t],
Derivative[OverTilde[c]][t]}]]}, {OverTilde[c][t],
OverTilde[c][t]}, {t, 0, 8000},
Method ->
"StiffnessSwitching"]]], {x, {-5 10^-4, -4 10^-4, -3 10^-4, -2 10^-4, -1.5 10^-4}}];


Visualization

ListPlot[Table[Flatten[sol[[i]][], 1], {i, Length[sol]}],
PlotTheme -> "Detailed",
PlotRange -> {{-  8 10^-4, 0}, {-4 10^-4, 4 10^-4}},
FrameLabel -> {OverTilde[Subscript[c, 0]],
OverDot[OverTilde[Subscript[c, 0]]]}, PlotStyle -> PointSize[.005]] • Thank you for such an elaborate answer! I have a doubt regarding the plot for x1 with B=0.6, it differs from the paper in regards to the location of scattered points. In the paper, near zero the orbits are scattered points which depend on the intial conditions. So, is it because of less number of initial conditions? I tried adding more values for x but got similar plot. May 3 at 6:50
• @codebpr As I understand in Figure 7 in the paper cited they shown chaotic solution for B=0.6. It is not clear why they supposed Poincare section with lines. Note, that for B=1 sections look like lines, but it is also scatter points. May 3 at 7:37
• Actually, in a similar paper(with B=0), they get the result with the orbits forming regular tori, with some tori near the saddle point of the potential collapsing to become scattered plots. Figure 4 of this paper on page 5: arxiv.org/pdf/1803.06756.pdf May 3 at 7:51
• I edited the code to add more initial conditions, I again get step size errors and I am getting this plot as output: imgur.com/a/W9FS6B3 It differs a little bit from the paper in terms of nonlinear resonances and island structures, and also the plot range is different. May 4 at 6:05
• @codebpr You did a very nice picture. The difference is due to differ scaling or may be there is a typo in the paper. May 4 at 7:01