Surface-of- section for a system of coupled oscillators

For the following system of non linear coupled ODEs:

\begin{align} & \ddot{x}=-x^{1/3}+k_2(y-x)^{1/3} \\ & \ddot{y}=-k_2(y-x)^{1/3} \end{align}

and the corresponding potential function:

\begin{equation} V(x,y)=\frac{3}{4}\left( x^{4/3}+k_2(y-x)^{4/3} \right) \end{equation}

I would like to produce some surfaces of section for various values of $k_2 \in \mathbb{R}$ which controls the coupling term.

For the sections, I have the following piece of code:

SurfaceOfSection[{x0_, {y0_, dy0_}, e_}, tmax_] :=
Module[
{dx0 =
Sqrt[2 e -dy0^2 - (2/(μ + 1)) (Abs[x0]^(μ + 1) + 2*Abs[y0 - x0]^(μ + 1))],
x, y, t},
If[# == {}, {}, First[#]] & @
Last[
Reap[
NDSolve[
{x''[t] == -CubeRoot[x[t]] + 2 CubeRoot[(y[t] - x[t])],
y''[t] == -2 CubeRoot[y[t] - x[t]],
y == y0, y' == dy0, x == x0, x' == dx0},
{x, y}, {t, 0, tmax},
Method ->
{EventLocator, "Event" -> x[t],
"EventAction" :> Sow[{y[t], y'[t]}]}]
Assumptions :> {Element[x, Reals], Element[y, Reals]}]]

where I chose in this case μ=1/3 and $k_2=2$, on the plane $y, \dot{y}$. As it can be seen, I have solved the energy equation with respect to $\dot{x}$:

\begin{equation} \dot{x}=\sqrt{2E-\dot{y}^2-\frac{3}{2}\left( |x|^{4/3}+k_2|(y-x)|^{4/3} \right)} \end{equation}

where I had to use Abs[x] and Abs[y - x] as well the CubeRoot[x[t]]^(μ), since Mathematica by default returns the principal value of a number raised to a rational exponent:

N[(-1)^(1/3)] = 0.5 + 0.866025 I

and not the real root which is needed in this case for the sections. (in fact I got an error message when I was trying to iterate without the Abs[]).

While I was able to (slowly) get some plots in the following way

x = Array[0 &, 151];
Do[
x[[i]] =
Join @@
Table[SurfaceOfSection[{0.001, {0, 0.01*(i - 1)}, 1.2}, 3000], {1}],
{i,151}]
l = Catenate[Table[x[[i]], {i, 151}]];
With[{x = .0001},
InternalDeactivateMessages[
ListPlot[l,
PlotStyle -> {PointSize[.001], Red},
AspectRatio -> GoldenRatio,
AxesLabel -> TraditionalForm /@ {y[t], Overscript[y, .][t]},
ImageSize -> 500, Frame -> True, Axes -> False]] I was wondering the following:

• How is it possible to acquire the correct result when evaluating x_0^(μ+1) and (y_0-x_0)^(μ+1) without using the absolute value of it? By doing so I am only getting back the positive corresponding point on the Sections which creates a false behavior on the dynamics of the ODEs. There must be a way to tell Mathematica to return the real root of these quantities.
• Is it possible to acquire the resulting plots faster? In order to get a satisfying result I had to use tmax=3000 which took like half an hour of calculations on my laptop.
• Try Surd, which gives the real root of any real quantity. – bbgodfrey May 26 '17 at 19:09
• Use Surd[-1, 3]. – user21 May 26 '17 at 19:10
• The correct way appears to be Surd[x0^(Numerator@(mu + 1)), Denominator@(mu + 1)], because the second argument must be an integer. – bbgodfrey May 26 '17 at 19:30
• NDSolve doesn't pay attention to assumptions, so just remove  Assumptions :> {Element[x, Reals], Element[y, Reals]` from your code. – m_goldberg May 26 '17 at 20:05
• By the way, "WhenEvent" is preferred over "EventLocator" in newer versions of Mathematica. – bbgodfrey May 26 '17 at 20:17

SurfaceOfSection can be simplified,

SurfaceOfSection[{x0_, {y0_, dy0_}, e_}, tmax_] :=
Reap[NDSolve[{x''[t] == -CubeRoot[x[t]] + 2 CubeRoot[(y[t] - x[t])],
y''[t] == -2 CubeRoot[y[t] - x[t]], y == y0, y' == dy0, x == x0,
x' == Sqrt[2 e - dy0^2 - (2/(4/3)) (Surd[x0^4, 3] + 2*Surd[(y0 - x0)^4, 3])],
WhenEvent[x[t] == 0, Sow[{y[t], y'[t]}]]}, {}, {t, 0, tmax}]][[2, 1]]

but the simplifications reduce runtime by only about 5%, mostly from replacing {x, y} by { }. As recommended in the comments, using Surd effectively addresses the first issue in the question. Then, reduction by a factor of three or more is achieved with ParallelTable (on a four-processor PC)

l = Join @@ ParallelTable[
SurfaceOfSection[{0.001, {0, 0.01*(i - 1)}, 1.2}, 3000], {i, 151}];

after which

ListPlot[l, PlotStyle -> {PointSize[.001], Red}, AspectRatio -> GoldenRatio,
ImageSize -> 500, Frame -> True, Axes -> False]

reproduces the plot in the question. For completeness Length@l is 113793.

Minor Generalization

In response to comments below, the earlier computation can be generalized as follows.

SurfaceOfSection[μ_, κ_, {x0_, {y0_, dy0_}, e_}, tmax_] :=
Reap[NDSolve[{x''[t] == -Surd[x[t]^(Numerator@μ), Denominator@μ] +
κ*Surd[(y[t] - x[t])^(Numerator@μ), Denominator@μ],
y''[t] == -κ*Surd[(y[t] - x[t])^(Numerator@μ), Denominator@μ],
y == y0, y' == dy0, x == x0,
x' == Sqrt[2 e - dy0^2 - (2/(μ + 1)) (Surd[x0^(Numerator@(μ + 1)),
Denominator@(μ + 1)] + κ*Surd[(y0 - x0)^(Numerator@(μ + 1)), Denominator@(μ + 1)])],
WhenEvent[x[t] == 0, Sow[{y[t], y'[t]}]]}, {}, {t, 0, tmax}]][[2, 1]]

Then, the results above are reproduced by

l = Join @@ ParallelTable[
SurfaceOfSection[1/3, 2, {0.001, {0, 0.01*(i - 1)}, 1.2}, 3000], {i, 151}];
• Seems like a new laptop would be of good use to me. :) Thank you @bbgodfrey! – Mitscaype May 27 '17 at 18:34
• Btw I tried to run and got the following error message Transpose::nmtx : The first two levels of {} cannot be transposed.. Is it something wrong that am I doing? – Mitscaype May 27 '17 at 18:51
• @Mitscaype What were you running? What version of Mathematica do you have? With respect to a new PC, if you buy one with four processors, I recommend you also obtain 16GB of memory. A solid state disk (in addition to a standard disk) is convenient but not necessary. If you do purchase a PC with SSD, be sure it is at least 200GB. – bbgodfrey May 27 '17 at 18:56
• @Mitscaype How are you inserting these values into SurfaceOfSection? Using {...} /. {μ = 1/3,κ = 2} does not work. However, With[{μ = 1/3, κ = 2}, ...] does. Perhaps better would be to add μ and κ as additional arguments of SurfaceOfSection. By the way, your expressions for x''[t] and y''[t] can be simplified. – bbgodfrey May 28 '17 at 17:31
• @Mitscaype Easy to do. Omit Join from the definition of l and Red from PlotStyle, and the set of points from each IC will have color chosen by ListPlot. However, after enough ICs, ListPlot will repeat colors. To prevent that, you must define your own table of PlotStyle values. – bbgodfrey May 31 '17 at 11:50