Plotting Arnold tongues

Question

I would like to understand how can I plot Arnold tongues (see figure below) with help of Wolfram Mathematica.

From this question, I have borrowed the function that computes the winding number of the circle map. My own version (in order to be consistent with the figure presented above) is

WindingNumber = 
  Compile[{n, η, ϵ, ϕ0}, (Nest[# + η + ϵ*Sin[#] &, ϕ0, n] - ϕ0)/(2*Pi*n)];

Here, η is the parameter, ϕ0 is the initial value and ϵ is the parameter, too. I know that each Arnold tongue corresponds to the rational value of WindingNumber. So, my idea is to check that for a given values of ϵ and η is rational or not. To do it, I need the Wolfram Mathematica analog of RationalQ function (see discussion here). As I understand, there is no analog of RationalQ function, but in this question, I see the naive realization:

TrueQ@Element[x, Rationals]

where (as I understand) one checks that x belongs to Rationals (or not). So, I have tried to create a grid of parameters η and ϵ with help of Table. For each point of the grid, I check is WindingNumber rational or not with the help of the function written above. If true, I mark this point by 1, if false, I mark by 0, so I write

ArnoldTongues = Table[{ϵ, 2*Pi*i, 
    If[TrueQ@Element[WindingNumber [100, 2*Pi*i, 1.0, 0], Rationals], 1, 
     0]}, {i, 0, 1, 0.01}, {ϵ, 0, 1, 0.01}];

As a result, I have the array that is suitable for ListDensityPlot.

However, I obtain that the winding number is always irrational for every point on my grid, which is impossible. What am I doing wrong?

Answer 1

(* Set compitation target, use "WVM" if no C compiler is avaliable *)
ClearAll[target] ;
target = "C" ;
ClearAll[windingNumber];
windingNumber = Compile[
    {{n,_Integer},{eta,_Real},{epsilon,_Real},{phi,_Real}},
    (Nest[# + eta + epsilon*Sin[#] &, phi, n] - phi)/(2*Pi*n),
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    CompilationTarget -> target,
    RuntimeOptions -> "Speed"
];

Generate wnums on a rectangular eta-epsilon grid:

(* Set number of iterations and grid size *)
n = 10^3;
m = 5*10^3;

(* Generate grid *)
eta = 2*Pi*Subdivide[0.0, 1.0, m - 1] ;
epsilon = Subdivide[0.0, 1.0, m - 1] ;

(* Set initial condition *)
phi = 0.0 ;

(* Compute winding numbers *)
wnums = Reverse[Table[windingNumber[n, eta, eps, phi], {eps, epsilon}]] ;

Plot all values:

(* Plot all winding numbers *)
ArrayPlot[wnums, DataRange -> {MinMax[eta], MinMax[epsilon]}, FrameTicks -> {{True, None}, {True, None}}, ImageSize -> 900, AspectRatio -> 1/4, PlotRangePadding -> None]

Filter and plot tongues:

(* Select close to a rational with a given tolerance *)
tol = 10.^-3;
rationals = {0, 1/5, 1/4, 1/3, 1/2, 2/3, 3/4, 4/5, 1} ;
ClearAll[select] ;
select[tol_, list_, rational_] := N[UnitStep[Clip[Abs[list - rational], {0.0, tol}, {-1.0, -1.0}]]] ; 
tongues = Fold[Plus, Map[select[tol, wnums, #] &, rationals]] ;
ArrayPlot[tongues, DataRange -> {MinMax[eta], MinMax[epsilon]}, FrameTicks -> {{True, None}, {True, None}}, ImageSize -> 900, AspectRatio -> 1/4, PlotRangePadding -> None]