# Find a stable 3-cycle of the sine map

I am trying to find an stable 3-cycle for the sine map $$g(x) = a \sin(\pi x)$$ but I do not know exactly how to use mathematica to do so. Is there anyone familiar with this?

$$x$$ lies as before between 0 and 1 but the real parameter $$a$$ is positive but not greater than 1. Use the bifurcation diagram to locate the 3-cycle

• Welcome to Mathematica.SE! There are a lot of posts on this site with answers explaining how to make bifurcation diagrams. Maybe start by searching for those and seeing if one of them helps you! Here's one. – march May 1 '19 at 18:37
• is there any other ways to do this more simply? – Guest1928 Apr 20 at 16:55

This is how you define the function $$g$$ in Mathematica:

g[a_, x_] = a*Sin[π*x];


Making a color-plot of the function $$g(g(g(x)))-x$$ in the $$a$$-$$x$$ plane and indicating the zero-contours, we see the lobes of 3-cycles:

DensityPlot[g[a, g[a, g[a, x]]] - x, {a, 0, 1}, {x, 0, 1},
MeshFunctions -> {#3 &}, Mesh -> {{0}}, PlotPoints -> 100] Now we know where they all are, find a 3-cycle for a specific value of $$a$$:

With[{a = 0.95},
x3 = x /. FindRoot[g[a, g[a, g[a, x]]] == x, {x, 0.2}];
{{x3, g[a, x3], g[a, g[a, x3]], g[a, g[a, g[a, x3]]]},
Abs[D[g[a, g[a, g[a, x]]], x]] /. x -> x3}]


{{0.201558, 0.562152, 0.931948, 0.201558}, 4.06318}

Oops, this is an unstable 3-cycle as the derivative (last number) is larger than 1 in magnitude (thanks @march!). Try again:

With[{a = 0.94},
x3 = x /. FindRoot[g[a, g[a, g[a, x]]] == x, {x, 0.15}];
{{x3, g[a, x3], g[a, g[a, x3]], g[a, g[a, g[a, x3]]]},
Abs[D[g[a, g[a, g[a, x]]], x]] /. x -> x3}]


{{0.176489, 0.494893, 0.939879, 0.176489}, 0.345007}

This time it worked: we found a stable 3-cycle $$0.176489, 0.494893, 0.939879$$ at $$a=0.94$$.

Let's make a new plot where the cycle contours are only shown in the stable regions (i.e., only stable 3-cycles):

g3[a_, x_] = Nest[g[a, #] &, x, 3];
dg3[a_, x_] = D[g3[a, x], x];
DensityPlot[g3[a, x] - x, {a, 0, 1}, {x, 0, 1},
MeshFunctions -> {#3 &}, Mesh -> {{0}}, MeshStyle -> Red,
PlotPoints -> 100,
RegionFunction -> Function[{a, x, f}, Abs[dg3[a, x]] <= 1]] Zooming in on one of the regions of stable 3-cycles: The two black dots that mark the limits of the stable 3-cycle band are found with

Q1 = {a, x} /. FindRoot[{g3[a, x] == x, dg3[a, x] == 1}, {{a, 0.94}, {x, 0.5}}]


{0.937818, 0.5152}

Q2 = {a, x} /. FindRoot[{g3[a, x] == x, dg3[a, x] == -1}, {{a, 0.94}, {x, 0.5}}]


{0.942488, 0.485444}

The region of stable 3-cycles is therefore $$a\in[0.937818, 0.942488]$$.

• I'm a little suspicious of your results here. I made the bifurcation diagram, and I think that 0.2 is in the region where there is exactly one fixed point. Furthermore the 3-cycles occupy a region in a-space. Something like $0.938 < a < 0.942$, roughly. – march May 1 '19 at 18:55
• @march you're right I didn't do a stability analysis, so what I wrote is too simplistic. – Roman May 1 '19 at 19:07
• This is cool! I wasn’t aware that you could do this kind of analysis. – march May 1 '19 at 19:31
• @Roman The colored background looks awesome, but I guess the only relevant bits are the zero contours, right? I also like the unstable 3-cycle on the other side of the stable one! – Chris K May 1 '19 at 19:50
• @ChrisK yes the colors are just to be fancy, they're otherwise useless. You could do a ContourPlot instead of a DensityPlot to go easier on the eyes if you want. – Roman May 1 '19 at 21:16

Here's a more brute-force approach to generate the whole bifurcation diagram using Nest and NestList.

g[x_] := a*Sin[π*x];

warmup = 200; (* # warmup iterations *)
pts = 200; (* # points to save *)
ic = 0.01; (* initial condition *)

ListPlot[Flatten[Table[
Transpose[{Table[a, pts],
NestList[g, Nest[g, ic, warmup], pts - 1]}],
{a, 0, 1, 0.001}], 1],
PlotStyle -> {Black, Opacity[0.2], PointSize[0.001]},
AxesLabel -> {a, x}] As @Roman found, the 3-cycle looks like it lives around a=0.94. We can zoom in to see it better:

warmup = 500;
pts = 500;

ListPlot[Flatten[Table[
Transpose[{Table[a, pts],
NestList[g, Nest[g, ic, warmup], pts - 1]}],
{a, 0.93, 0.95, 0.00001}], 1],
PlotStyle -> {Black, Opacity[0.1], PointSize[0.001]},
AxesLabel -> {a, x}] 