# 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:

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

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[0] == y0, y'[0] == dy0, x[0] == x0, x'[0] == 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}$:

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

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. Commented May 26, 2017 at 19:09
• Use Surd[-1, 3]. Commented May 26, 2017 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. Commented May 26, 2017 at 19:30
• NDSolve doesn't pay attention to assumptions, so just remove  Assumptions :> {Element[x, Reals], Element[y, Reals] from your code. Commented May 26, 2017 at 20:05
• By the way, "WhenEvent" is preferred over "EventLocator" in newer versions of Mathematica. Commented May 26, 2017 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[0] == y0, y'[0] == dy0, x[0] == x0,
x'[0] == 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[0] == y0, y'[0] == dy0, x[0] == x0,
x'[0] == 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! Commented May 27, 2017 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? Commented May 27, 2017 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. Commented May 27, 2017 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. Commented May 28, 2017 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. Commented May 31, 2017 at 11:50