5
$\begingroup$

Note: This question is for those chaotic systems where the order of maximal Lyapunov exponent is around $10^{-3}$. So it's different from normal chaotic systems like Hénon–Heiles system, Yang-Mills-Higgs System, Double Pendulum or the famous Lorenz Attractor.

The related problem: I am trying to reproduce my results for the Poincaré section published here (Specifically one of the sub-figures in Figure $8$ on page $15$ given below), which I plotted by making use of an algorithm developed by @Alex Trounev. I am trying to plot the same this time by a different interactive algorithm of @E. Chan-López.

poincare-string

The code being used is as follows:

lagrangian = (8 x[t]^2)/3 + (4 x[t]^3)/3 - 6 y[t]^2 + 148/7 x[t] y[t]^2 + Derivative[1][x][t]^2 + 1/3 y[t] Derivative[1][x][t] Derivative[1][y][t] + Derivative[1][y][t]^2 + 1/5 x[t] Derivative[1][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; {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}; constraint = FullSimplify[ Last[(py0 /. Solve[e == (hamiltonian /. {x[t] -> x0, px[t] -> px0, y[t] -> 0, py[t] -> py0}), py0])]] PoincareSections[{x0_, px0_, y0_, py0_}, e_, tmax_] := Block[{constraint, x, px, y, py}, constraint = ( I Sqrt[5 + x0] Sqrt[-12 e + 3 px0^2 - 16 x0^2 (2 + x0)])/Sqrt[15]; If[Head[N[constraint]] === Complex, Nothing, If[# == {}, {}, First[#]] &@ Last[Reap[ NDSolve[{hameqns, x[0] == x0, px[0] == px0, y[0] == y0, py[0] == constraint, WhenEvent[y[t], Sow[{x[t], px[t]}, "LocationMethod" -> "EvenLocator"]]}, {x, px, y, py}, {t, 0, tmax}, MaxSteps -> ∞, MaxStepSize -> 0.1]]]]] style = {PlotStyle -> {{AbsolutePointSize[1], Opacity[0.4]}}, AspectRatio -> 1, ImageSize -> 420, PlotHighlighting -> None}; DynamicModule[{f = {}}, Column[{Dynamic[ MatrixForm@{MousePosition["Graphics", "Mouse not in graphics!"]}], Row[{ClickPane[ Dynamic@ListPlot[f, style], (AppendTo[f, PoincareSections[{#[[1]], #[[2]], 0, 0}, 10^-5, 3000]]) &], Button[Style["Copy orbit list", 12], CopyToClipboard[Iconize[f, "Orbit list"]], ImageSize -> Automatic, Alignment -> Center, Appearance -> {"Palette", "Pressed"}]}, " "]}, Alignment -> Center]]

Now, the above algorithm doesn't work simply because the Head of the constraint turns out to be complex. But if we calculate it, technically, it's not a complex number. A simple example can explain the reason:

The issue:

a = N[constraint] /. x0 -> -1/1000 /. px0 -> 1/1000 /. e -> 10^-5 (*-0.007046375030231264`+0.` I*) Head[a] (*Complex*)

The number is Real, but the Head turns out to be Complex only because it is in the form of a+Ib.

If I comment out the condition If[Head[N[constraint]] === Complex in the code, it works, but the system hangs after 3 or 4 mouse clicks. How to rectify this issue?

Edit: Issue resolved Thanks to @BobHanlon, the issue is solved.

 PoincareSections[{x0_, px0_, y0_, py0_}, e_, tmax_] := Block[{constraint, x, px, y, py}, constraint = N[(I Sqrt[ 5 + x0] Sqrt[-12 e + 3 px0^2 - 16 x0^2 (2 + x0)])/ Sqrt[15]] // Chop; If[Head[N[constraint]] === Complex, Nothing, If[# == {}, {}, First[#]] &@ Last[Reap[ NDSolve[{hameqns, x[0] == x0, px[0] == px0, y[0] == y0, py[0] == constraint, WhenEvent[y[t], Sow[{x[t], px[t]}, "LocationMethod" -> "EvenLocator"]]}, {x, px, y, py}, {t, 0, tmax}, MaxSteps -> ∞, MaxStepSize -> 0.1]]]]] style = {PlotStyle -> {{PointSize[.005], Opacity[0.9]}}, AspectRatio -> 1, ImageSize -> 420, PlotHighlighting -> None}; DynamicModule[{f = {}}, Column[{Dynamic[ MatrixForm@{MousePosition["Graphics", "Mouse not in graphics!"]}], Row[{ClickPane[ Dynamic@ListPlot[f, style], (AppendTo[f, PoincareSections[{#[[1]], #[[2]], 0, 0}, 10^-5, 3000]]) &], Button[Style["Copy orbit list", 12], CopyToClipboard[Iconize[f, "Orbit list"]], ImageSize -> Automatic, Alignment -> Center, Appearance -> {"Palette", "Pressed"}]}, " "]}, Alignment -> Center]]

Choosing the best NDSolve Method

Rectified code used for a close up:

par = {9.07, 11.73, 5.80, 2.18, 4.36, -2.645, 6.023}; {K1, K2, K3, K4, K5, ωsq[0], ωsq[1]} = Rationalize[par, 10^-30]; lagrangianold = Sum[c[n]'[t]^2 - c[n][t]^2 ωsq[n], {n, {0, 1}}] + K1 c[0][t]^3 + K2 c[0][t] c[1][t]^2 + K3 c[0][t] c[0]'[t]^2 + K4 c[0][t] c[1]'[t]^2 + K5 c[0]'[t] c[1][t] c[1]'[t]; c[0][t_] := OverTilde[c][0][t] + α1*OverTilde[c][0][t]^2 + α2* OverTilde[c][1][t]^2; c[1][t_] := OverTilde[c][1][t] + α3*OverTilde[c][0][t]*OverTilde[c][1][t]; α1 = -29/20; α2 = -1/2; α3 = -1; n = Expand[lagrangianold]; vars = {OverTilde[c][0][t], OverTilde[c][1][t], Derivative[1][OverTilde[c][0]][t], Derivative[1][OverTilde[c][1]][t]}; lagrangianew = Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1; lagrangian = Simplify[ lagrangianew /. OverTilde[c][0] -> x /. OverTilde[c][1] -> y]; 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; {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}; constraint = FullSimplify[ Last[(py0 /. Solve[e == (hamiltonian /. {x[t] -> x0, px[t] -> px0, y[t] -> 0, py[t] -> py0}), py0])]] PoincareSections[{x0_, px0_, y0_, py0_}, e_, tmax_] := Block[{constraint, x, px, y, py}, constraint = Chop[N[1/50 Sqrt[-5 - (9 x0)/10] Sqrt[-2000 e + 500 px0^2 - x0^2 (5290 + 2799 x0)]]]; If[Head[N[constraint]] === Complex, Nothing, If[# == {}, {}, First[#]] &@ Last[Reap[ NDSolve[{hameqns, x[0] == x0, px[0] == px0, y[0] == y0, py[0] == constraint, WhenEvent[y[t], Sow[{x[t], px[t]}, "LocationMethod" -> "EvenLocator"]]}, {x, px, y, py}, {t, 0, tmax}, MaxSteps -> ∞, MaxStepSize -> 0.1, Method -> "StiffnessSwitching"]]]]] style = {PlotStyle -> {{PointSize[.005], Opacity[0.9]}}, AspectRatio -> 1, ImageSize -> 300, PlotRange -> {{-0.0001, 0}, {-0.0001, 0.0001}}, PlotHighlighting -> None}; Off[NDSolve::ndsz] DynamicModule[{f = {}}, Column[{Dynamic[ MatrixForm@{MousePosition["Graphics", "Mouse not in graphics!"]}], Row[{ClickPane[ Dynamic@ListPlot[f, style], (AppendTo[f, PoincareSections[{#[[1]], #[[2]], 0, 0}, 10^-5, 2000]]) &], Button[Style["Copy orbit list", 12], CopyToClipboard[Iconize[f, "Orbit list"]], ImageSize -> Automatic, Alignment -> Center, Appearance -> {"Palette", "Pressed"}]}, " "]}, Alignment -> Center]]

I have one minor problem now for a close-up of the above section: I get an NDSolve warning message, which I turn off. Still, the plot I get differs from the original when I use Method -> "ImplicitRungeKutta" instead of Automatic or some other like "StiffnessSwitching".

Here is a result found using two different methods for comparison:

comparison

Why is it that the results are scaled in the y-axis, and which method should be considered the best in such a case? The methods that I used are also not easily seen in the documentation. Also, any help in making the code more elegant(like avoiding warnings) would be genuinely beneficial!

Improve this question
$\endgroup$
11
  • 4
    $\begingroup$ In the definition of PoincareSections change the definition of constraint to read constraint = N[(I Sqrt[5 + x0] Sqrt[-12 e + 3 px0^2 - 16 x0^2 (2 + x0)])/ Sqrt[15]] // Chop; The Chop will remove imaginary artifacts. $\endgroup$ Commented Mar 5, 2024 at 5:25
  • $\begingroup$ @BobHanlon, thanks a lot. It works now. Can you see the updated question for a minor query? $\endgroup$ Commented Mar 5, 2024 at 7:50
  • $\begingroup$ @codebpr It could be better to reproduce known results, maybe from your paper as well. $\endgroup$ Commented Mar 5, 2024 at 14:20
  • $\begingroup$ @AlexTrounev should I ask a new question for that or just update this one to reflect the comparison for reproduced results? $\endgroup$ Commented Mar 5, 2024 at 15:19
  • 1
    $\begingroup$ @I.M. thanks for the useful response. I was just wondering why the Wolfram guys have not kept them in the Options section of NDSolve but hidden inside the Tutorials? Also, one query? Why is the result rescaled when I use Hamilton equations of motion instead of Euler- Lagrange? $\endgroup$ Commented Mar 6, 2024 at 6:23

1 Answer 1

Reset to default
6
+100
$\begingroup$

We don't know what method is most efficient, since to answer this question we need too much test performed. We can only suggest code without some NDSolve messages concerning solution behaver and WhenEvent. It seems impossible to exclude all messages. We can't avoid, for example, message which appears when we put mouse on an axis

General::munfl: -(1./9.49657*10^307) is too small to represent as a normalized machine number; precision may be lost.

First, note, that for chaotic system research we should exclude any Automatic option, like automatic step size. Second, note, that all implicit scheme working in this case like averaging, and therefore depresses chaotic mode. Third, note, that WhenEvent with some DetectionMethod make code unstable. Taken into account all these notes we proposed next options

par = {9.07, 11.73, 5.80, 2.18, 4.36, -2.645, 6.023}; {K1, K2, K3, K4, K5, \[Omega]sq[0], \[Omega]sq[1]} = Rationalize[par, 10^-30]; lagrangianold = Sum[c[n]'[t]^2 - c[n][t]^2 \[Omega]sq[n], {n, {0, 1}}] + K1 c[0][t]^3 + K2 c[0][t] c[1][t]^2 + K3 c[0][t] c[0]'[t]^2 + K4 c[0][t] c[1]'[t]^2 + K5 c[0]'[t] c[1][t] c[1]'[t]; c[0][t_] := OverTilde[c][0][t] + \[Alpha]1*OverTilde[c][0][t]^2 + \[Alpha]2* OverTilde[c][1][t]^2; c[1][t_] := OverTilde[c][1][t] + \[Alpha]3*OverTilde[c][0][t]*OverTilde[c][1][t]; \[Alpha]1 = -29/20; \[Alpha]2 = -1/2; \[Alpha]3 = -1; n = Expand[lagrangianold]; vars = {OverTilde[c][0][t], OverTilde[c][1][t], Derivative[1][OverTilde[c][0]][t], Derivative[1][OverTilde[c][1]][t]}; lagrangianew = Normal[Series[n /. Thread[vars -> m*vars], {m, 0, 3}]] /. m -> 1; lagrangian = Simplify[ lagrangianew /. OverTilde[c][0] -> x /. OverTilde[c][1] -> y]; 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; {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}; constraint = FullSimplify[ Last[(py0 /. Solve[e == (hamiltonian /. {x[t] -> x0, px[t] -> px0, y[t] -> 0, py[t] -> py0}), py0])]]; PoincareSections[{x0_, px0_, y0_, py0_}, e_, tmax_] := Block[{constraint, x, px, y, py}, constraint = Chop[N[1/ 50 Sqrt[-5 - (9 x0)/10] Sqrt[-2000 e + 500 px0^2 - x0^2 (5290 + 2799 x0)]]]; If[Head[N[constraint]] === Complex, Nothing, If[# == {}, {}, First[#]] &@ Last[Reap[ NDSolve[{hameqns, x[0] == x0, px[0] == px0, y[0] == y0, py[0] == constraint, WhenEvent[y[t], Sow[{x[t], px[t]}], "DetectionMethod" -> "Interpolation"]}, {x, px, y, py}, {t, 0, tmax}, MaxSteps -> \[Infinity], StartingStepSize -> 0.1, Method -> {"FixedStep", Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 8}}]]]]] style = {PlotStyle -> {{PointSize[.005], Opacity[0.5]}}, AspectRatio -> 1, ImageSize -> 500, PlotRange -> {{-0.0001, 0}, {-0.0001, 0.0001}}, PlotHighlighting -> None, Frame -> True, Background -> Black}; (*Off[NDSolve::ndsz]*) DynamicModule[{f = {}}, Column[{Dynamic[ MatrixForm@{MousePosition["Graphics", "Mouse not in graphics!"]}], Row[{ClickPane[ Dynamic@ListPlot[f, style], (AppendTo[f, PoincareSections[{#[[1]], #[[2]], 0, 0}, 10^-5, 2000]]) &], Button[Style["Copy orbit list", 12], CopyToClipboard[Iconize[f, "Orbit list"]], ImageSize -> Automatic, Alignment -> Center, Appearance -> {"Palette", "Pressed"}]}, " "]}, Alignment -> Center]]

Figure 1

Using this code we can performed research with step size decreasing.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Many thanks for this! Just a small query, when I tried testing with SPRK method, I don't get any points. Here's what I used: Method -> {"SymplecticPartitionedRungeKutta", "DifferenceOrder" -> 8, "PositionVariables" -> {x[t], y[t]}} $\endgroup$ Commented Mar 8, 2024 at 4:19
  • 1
    $\begingroup$ @codebpr Hamiltonian hamiltonian is not separable. Therefore SPRK has been replaced in this case by Automatic method. In turn Automatic method is optimized for variable step size that leads to "StiffnessSwitching" and automatic Stop. $\endgroup$ Commented Mar 8, 2024 at 12:39

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.