Poincaré section of the particle motion near the black hole horizon

Question

I am trying to recreate the Poincaré sections given in Figure 2 of this paper.

For this purpose, I use the following two methods:

Method1: Using ResourceFunction["ClickPoincarePlot2D"]

    f = 2  k  x[t];

lagrangian = -Sqrt[f - x'[t]^2/f - y'[t]^2] - \[Omega]^2/
    2  ((x[t] - xc)^2 + y[t]^2);

px1 = D[lagrangian, x'[t]]

py1 = D[lagrangian, y'[t]]

testham = (px1  x'[t] + py1  y'[t]) - lagrangian

k = 1/2; \[Omega] = 10; xc = 1;

hamiltonian[x_, px_, y_, py_] = 
  testham /. x[t] -> x /. x'[t] -> px /. y[t] -> y /. y'[t] -> py;


H := Function[{S}, (2  k  S[[1]])/Sqrt[
   2  k  S[[1]] - S[[2]]^2/(2  k  S[[1]]) - S[[4]]^2] + 
   1/2  (S[[1]]^2 + S[[3]]^2 - 2  S[[1]]  xc + xc^2)  \[Omega]^2]
    
    {a1, b1, c1, d1} = {D[hamiltonian[x[t], px[t], y[t], py[t]], px[t]], 
      D[hamiltonian[x[t], px[t], y[t], py[t]], 
       py[t]], -D[hamiltonian[x[t], px[t], y[t], py[t]], x[t]], -D[
        hamiltonian[x[t], px[t], y[t], py[t]], y[t]]}; abcd = {x'[t] == 
       a1, px'[t] == c1, y'[t] == b1, py'[t] == d1};
    
    cross = y[t]; recover = {x[t], px[t]};
    
    ResourceFunction[
      "ClickPoincarePlot2D"][abcd, H, 15, t, 10000, cross, recover, \
    {PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
      PlotHighlighting -> None, PlotRange -> {{0, 2}, {-20, 20}}}]

Method2: Using algorithm of this answer by @Alex Trounev

psect[{x0_, px0_, y0_, py0_}] := 
 Reap[NDSolve[{abcd, x[0] == x0 + 10^-10, y[0] == y0, px[0] == px0, 
     py[0] == py0, WhenEvent[x[t] == 0, Sow[{y[t], py[t]}]]}, {}, {t, 
     0, 1000}, MaxSteps -> \[Infinity]]][[-1, 1]]

inivalues[E0_, n_, pxmin_, pxmax_] := 
  Module[{x, px, X, y, P, py, Y}, px = RandomReal[{pxmin, pxmax}, n];
   X = Table[0, {n}]; Y = Table[0, {n}]; P = Table[0, {n}];
   Do[{X[[i]], P[[i]], Y[[i]]} = {x, py, y} /. 
       NMinimize[(hamiltonian[x, px[[i]], y, py] - E0)^2, {x, py, 
          y}][[2]];, {i, n}]; Transpose[{X, px, Y, P}]];

iniv = inivalues[15, 100, -0.1, 0.1];

H0[{x_, px_, y_, py_}] = -Sqrt[f - px^2/f - py^2] - \[Omega]^2/
    2 ((x - xc)^2 + y^2) + px^2 + py^2;

abcdata = Map[psect, iniv];

ListPlot[abcdata, ImageSize -> Large, Background -> White, 
 Frame -> True]

I tried to alter the cross and recover variables in both methods but to no avail! The first method gives me an empty plot while the second one provides me with errors!

Any help in this regard would be truly beneficial!

ResourceFunction["ClickPoincarePlot2D"] des not work for me. Even the examples are not working. MMA version 14 — Daniel Huber, Commented Jan 21 at 12:00
The examples are working on MMA version 13. It is developed by @E. Chan-López, who is also a regular user here, he told about his tool in this answer: mathematica.stackexchange.com/a/290237/78049 — codebpr, Commented Jan 21 at 12:39
I'm working on an update and solving a significant issue that I observed together with Alex Trounev in order to generalize the routine so that it can be applied to more complicated systems. I estimate that I will be able to send an update by the end of February. I appreciate the feedback. — E. Chan-López, Commented Jan 21 at 17:45
Suddenly it stopped working for my MMA version 13 as well! Seems like they are working on it. But what about the second method by @Alex Trounev ? Why am I getting errors for that? — codebpr, Commented Jan 22 at 5:01
@codebpr I'm going to make an interactive version using ClickPane and post it as an answer. Once I have updated the resource function and the new version is accepted in WFR, then I will repost your example using the resource function. — E. Chan-López, Commented Jan 22 at 19:40

E. Chan-López · Accepted Answer · 2024-01-24 03:32:24Z

To detect more islands, it is better to use the interactive method. Here, I provide an example with the result for the case in which the energy is equal to 40.

PoincareSections[{x0_, px0_, y0_, py0_}, e_, tmax_] := 
 Block[{pdy0, x, px, y, py},
  pdy0 = 
   Sqrt[2500 + e^2 - 100*e*(-1 + x0)^2 + 
      x0*(-10001 - px0^2*x0 + 2500*x0*(6 + (-4 + x0)*x0))]/Sqrt[x0];
  If[Head[N[pdy0]] === Complex, Nothing,
   If[# == {}, {}, First[#]] &@Last[Reap[NDSolve[{
        x'[t] == px[t] x[t] Sqrt[x[t]/(1 + py[t]^2 + px[t]^2 x[t])],
        px'[t] == 
         100 - 100 x[t] - 
          px[t]^2 Sqrt[x[t]/(
           1 + py[t]^2 + px[t]^2 x[t])] - ((1 + py[t]^2) Sqrt[x[t]/(
           1 + py[t]^2 + px[t]^2 x[t])])/(2 x[t]),
        y'[t] == py[t] Sqrt[x[t]/(1 + py[t]^2 + px[t]^2 x[t])],
        py'[t] == -100 y[t],
        x[0] == x0,
        px[0] == px0,
        y[0] == y0,
        py[0] == pdy0,
        WhenEvent[y[t], 
         Sow[{x[t], px[t]}, "LocationMethod" -> "EvenLocator"]]
        },
       {x, px, y, py}, {t, 0, tmax}, MaxSteps -> ∞, 
       MaxStepSize -> 0.1
       ]]]]]

Next, use the code with ClickPane found at: The Poincaré Sections of Yang-Mills-Higgs System:

    style = {PlotStyle -> {{AbsolutePointSize[1], GrayLevel[0], 
             Opacity[0.4]}}, AspectRatio -> 1, ImageSize -> 420};
(*Omit the PlotHighlighting option if you're using a version prior to 13.3.*)

    DynamicModule[{f = {}}, 
     Column[{Dynamic[
        MatrixForm@{MousePosition["Graphics", "Mouse not in graphics!"]}],
        Row[{ClickPane[
          Dynamic@ListPlot[f, 
            style], (AppendTo[f, 
             PoincareSections[{#[[1]], #[[2]], 0, 0}, 40, 3000]]) &], 
         Button[Style["Copy orbit list", 12], 
          CopyToClipboard[Iconize[f, "Orbit list"]], 
          ImageSize -> Automatic, Alignment -> Center, 
          Appearance -> {"Palette", "Pressed"}]}, "  "]}, 
      Alignment -> Center]]

The Poincaré section plot looks like this:

style ={PlotStyle -> {{Black, Opacity[0.5], PointSize[0.002]}}, ImageSize -> Large};

ListPlot["Orbit list", style]

Note that you can copy the list of sections by clicking on the button that says Copy orbit list. Afterward, you can save your generated image in PNG format.

Interactively:

The plots of the Poincaré sections for energies 15 and 49.5 respectively look like this:

To export your images with high quality, I recommend the following post: Obtaining better quality in ListPlot.

@E.Chan-López thank you for the interactive method. Just a small question, what does pdy0 denote in your code? Have you solved for the Hamiltonian to get it, where e is the energy of the system? — codebpr, Commented Jan 23 at 4:41
@codebpr The variable pdy0 is the initial condition that is updated every time you click on the plane (x, px) considered in WhenEvent. — E. Chan-López, Commented Jan 23 at 5:04
@E.Chan-López thanks for clarifying. One more doubt, can I export the image as a png file from this interactive plot? Or in other words, how did you get the first image in your answer? — codebpr, Commented Jan 23 at 5:17
@codebpr I hope you generate high-quality images and share the result with the community. I'll move on to placing an additional answer here when I've updated the ClickPoincarePlot2D resource function. — E. Chan-López, Commented Jan 23 at 5:44

I.M. · Accepted Answer · 2024-01-23 02:57:09Z

There are several problems with your code:

Hamiltonian construction is not correct (you need to solve for velocities in terms of momenta, not just substitute $x' \to p_x$ and $y' \to p_y$)
In the linked paper the section is in $(x, p_x)$, while your event is WhenEvent[x[t] == 0, Sow[{y[t], py[t]}]] which is wrong

Note, while solving for $p_x$ and $p_y$ I just took one of the branches. Similarly solving for $p_y$ for given $(x, p_x)$ and $y = 0$, I just assume the last solution is positive. Sections looks close, but you need to check my assumptions.

k = 1/2 ;
w = 10 ; 
c = 1 ;

lagrangian = -Sqrt[2 k x[t] - x'[t]^2/(2 k x[t]) - y'[t]^2] - w^2/2 ((x[t] - c)^2 + y[t]^2) ;

Px = D[lagrangian, x'[t]] ;
Py = D[lagrangian, y'[t]] ;
hamiltonian = ((Px x'[t] + Py y'[t]) - lagrangian) /. (Solve[{px[t] == Px, py[t] == Py}, {x'[t], y'[t]}] // Last) // FullSimplify ;

equations = {
    x'[t] == +D[hamiltonian, px[t]],
    y'[t] == +D[hamiltonian, py[t]],
    px'[t] == -D[hamiltonian, x[t]], 
    py'[t] == -D[hamiltonian, y[t]]
} // Simplify ;


ClearAll[solution] ;
solution[X_, PX_] := Reap@NDSolve[
        Flatten[
            {
                equations, 
                {x[0] == X, px[0] == PX, y[0] == 0, py[0] == initial[X, PX]},
                WhenEvent[y[t] == 0 && y'[t] > 0, Sow[{x[t], px[t]}]]
            }
        ],
        {},
        {t, 0.0, T},
        MaxSteps -> Infinity
] ;

T = 2000.0 ;

energy = 15.0 ;
constraint = (py /. Solve[energy == hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> 0, py[t] -> py}, py]) // Last ;
ClearAll[initial] ;
initial[x_?NumericQ, px_?NumericQ] := Evaluate[constraint]
ListPlot[Table[solution[X, 0.0] // Last // First, {X, 0.5, 1.5, 0.01}], PlotRange->All, PlotTheme -> "Detailed", PlotStyle -> Directive[{PointSize[Tiny], Black}], AspectRatio -> 1/GoldenRatio, ImageSize->Large, GridLines->None]

energy = 40.0 ;
constraint = (py /. Solve[energy == hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> 0, py[t] -> py}, py]) // Last ;
ClearAll[initial] ;
initial[x_?NumericQ, px_?NumericQ] := Evaluate[constraint]
ListPlot[Table[solution[X, 0.0] // Last // First, {X, 0.2, 1.5, 0.01}], PlotRange->All, PlotTheme -> "Detailed", PlotStyle -> Directive[{PointSize[Tiny], Black}], AspectRatio -> 1/GoldenRatio, ImageSize->Large, GridLines->None]

energy = 49.5 ;
constraint = (py /. Solve[energy == hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> 0, py[t] -> py}, py]) // Last ;
ClearAll[initial] ;
initial[x_?NumericQ, px_?NumericQ] := Evaluate[constraint]
ListPlot[Table[solution[X, 0.0] // Last // First, {X, 0.2, 1.5, 0.01}], PlotRange->All, PlotTheme -> "Detailed", PlotStyle -> Directive[{PointSize[Tiny], Black}], AspectRatio -> 1/GoldenRatio, ImageSize->Large, GridLines->None]

It is a nice solution (+1). Maybe it could be better to variate y[0] as well. — Alex Trounev, Commented Jan 22 at 11:38
@AlexTrounev, right, it seems they've used $y=0.1$ (see page 5) — I.M., Commented Jan 22 at 11:48
In general case it could be better to use random values as in my answer. — Alex Trounev, Commented Jan 22 at 12:08
@AlexTrounev, could you explain about using random being beneficial? From my understanding with random initials you hit more islands? — I.M., Commented Jan 22 at 17:44

Alex Trounev · Accepted Answer · 2024-01-22 12:06:08Z

As mentioned by I.M. the Hamiltonian could be derive with using Legendre transformation - see for example here. Therefore we have

k = 1/2;
w = 10;
c = 1;

lagrangian = -Sqrt[2  k  x[t] - x'[t]^2/(2  k  x[t]) - y'[t]^2] - 
   w^2/2  ((x[t] - c)^2 + y[t]^2);

Px = D[lagrangian, x'[t]];
Py = D[lagrangian, y'[t]];
hamiltonian = ((Px  x'[t] + Py  y'[t]) - 
     lagrangian) /. (Solve[{px[t] == Px, py[t] == Py}, {x'[t], 
       y'[t]}] // First) /. {x[t] -> x, px[t] -> px, y[t] -> y, 
   py[t] -> py}

Last expression we use in two forms

h[x_, px_, y_, py_] := (
   py^2 x)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + (
   px^2 x^2)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + Sqrt[
   x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
    1 + py^2 + px^2 x)] + 50  ((-1 + x)^2 + y^2);

H0[{x_, px_, y_, py_}] := (
   py^2 x)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + (
   px^2 x^2)/((1 + py^2 + px^2 x) Sqrt[
    x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
     1 + py^2 + px^2 x)]) + Sqrt[
   x - (py^2 x)/(1 + py^2 + px^2 x) - (px^2 x^2)/(
    1 + py^2 + px^2 x)] + 50  ((-1 + x)^2 + y^2);

In this case we can generate initial values using NMiimize as follows

inivalues[E0_, n_, xmin_, xmax_] := 
  Module[{X, x, px, PX, y, P, py, Y}, 
   x = RandomReal[{xmin, xmax}, n];
   PX = Table[0, {n}]; Y = Table[0, {n}]; P = Table[0, {n}];
   Do[{PX[[i]], P[[i]], Y[[i]]} = {px, py, y} /. 
       NMinimize[(h[x[[i]], px, y, py] - E0)^2, {px, py, 
          y}][[2]];, {i, n}]; Transpose[{x, PX, Y, P}]];

For E0=15 we have iniv = inivalues[15, 100, .5, 1.5];

Map[H0, iniv]

(*Out[]= {15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., 15., \
15., 15., 15., 15.}*)

Therefore the energy is saved for every run. Last portion of code is same as in my answer here, also we use {x[t],px[t]} for visualization and y[t] == 0 && y'[t] > 0 for events detection as in a paper linked

{a1, b1, c1, d1} = {D[h[x[t], px[t], y[t], py[t]], px[t]], 
  D[h[x[t], px[t], y[t], py[t]], 
   py[t]], -D[h[x[t], px[t], y[t], py[t]], x[t]], -D[
    h[x[t], px[t], y[t], py[t]], y[t]]}; abcd = {x'[t] == a1, 
  px'[t] == c1, y'[t] == b1, py'[t] == d1};
psect[{x0_, px0_, y0_, py0_}] := 
 Reap[NDSolve[{abcd, x[0] == x0, y[0] == y0, px[0] == px0, 
     py[0] == py0, 
     WhenEvent[y[t] == 0 && y'[t] > 0, Sow[{x[t], px[t]}]]}, {}, {t, 
     0, 1000}, MaxSteps -> 10^6]][[-1, 1]]

Finally we have

abcdata = Map[psect, iniv];

ListPlot[abcdata, ImageSize -> Large, Background -> Black, 
 Frame -> True]

For E0=40 we have picture

And for E0=49.6

Thank you again for providing such an elaborate solution :) I hope the problem with resources.wolframcloud.com/FunctionRepository/resources/… also gets rectified soon, it is not working for any MMA version, as of now. — codebpr, Commented Jan 22 at 12:12

codebpr · Accepted Answer · 2024-03-04 13:57:39Z

The resource function ["ClickPoincarePlot2D"] is also working now, thanks to the new update.

I have also used the ScalarPureFunction by @E. Chan-López presented here to simplify my Hamiltonian into a scalar function form.

ScalarPureFunction[f_, 
    argvars_?
     VectorQ] /; (SameQ[Length[Variables[Level[f, {-1}]]], 
      Length[Flatten[argvars]]] && Length[argvars] > 1) := 
  Module[{positions, heldPart, replacements}, 
   positions = MapIndexed[List, argvars, {-1}];
   heldPart = 
    MapThread[
     Function[
      HoldForm[Part][\[FormalP], 
       Apply[Sequence, Part[SlotSequence@1, 2, All]]]], {positions}];
   replacements = Thread[argvars -> heldPart];
   ReleaseHold[
    Fold[Function[#, #2] &, {\[FormalP]}, {f /. replacements}]]];

ScalarPureFunction[f_][argvars_?VectorQ] := 
  ScalarPureFunction[f, argvars];

k = 1/2;
w = 10;
c = 1;

lagrangian = -Sqrt[2  k  x[t] - x'[t]^2/(2  k  x[t]) - y'[t]^2] - 
   w^2/2  ((x[t] - c)^2 + y[t]^2);

Px = D[lagrangian, x'[t]];
Py = D[lagrangian, y'[t]];
hamiltonian = ((Px  x'[t] + Py  y'[t]) - 
      lagrangian) /. (Solve[{px[t] == Px, py[t] == Py}, {x'[t], 
        y'[t]}] // Last) // FullSimplify;

H := ScalarPureFunction[
  hamiltonian /. {x[t] -> x, px[t] -> px, y[t] -> y, py[t] -> py}, {x,
    px, y, py}]

{a1, b1, c1, d1} = {D[hamiltonian, px[t]], 
  D[hamiltonian, 
   py[t]], -D[hamiltonian, x[t]], -D[hamiltonian, 
    y[t]]}; hameqns = {x'[t] == a1, px'[t] == c1, y'[t] == b1, 
  py'[t] == d1};

cross = y[t]; recover = {x[t], px[t]};

E=15

ResourceFunction[
          "ClickPoincarePlot2D"][hameqns, H, 15, t, 6000, cross, recover,
        {PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
          PlotRange -> {{0, 2}, Automatic}, PlotHighlighting -> None}]

E=40

ResourceFunction
  "ClickPoincarePlot2D"][hameqns, H, 40, t, 6000, cross, recover,
{PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
  PlotRange -> {{0, 2}, Automatic}, PlotHighlighting -> None}]

E=49.5

ResourceFunction
      "ClickPoincarePlot2D"][hameqns, H, 49.5, t, 6000, cross, recover,
    {PlotStyle -> AbsolutePointSize[1], AspectRatio -> 1, 
      PlotRange -> {{0, 2}, Automatic}, PlotHighlighting -> None}]

Stack Exchange Network

Poincaré section of the particle motion near the black hole horizon

4 Answers 4

Not the answer you're looking for? Browse other questions tagged
differential-equations
equation-solving
calculus-and-analysis
physics
sow-reap
or ask your own question.

Linked

Hot Network Questions

Poincaré section of the particle motion near the black hole horizon

4 Answers 4

Not the answer you're looking for? Browse other questions tagged differential-equationsequation-solvingcalculus-and-analysisphysicssow-reap or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
differential-equations
equation-solving
calculus-and-analysis
physics
sow-reap
or ask your own question.