How to create a version of FixedPointList for asymptotically periodic functions?

Use built in hash find duplicate after rounding to precision

reapList::usage = "remove unneeded generality from Reap, do correct thing upon zero Sow";
SetAttributes[reapList, HoldAll];
reapList[r_] := Join @@ Reap[r][[2]]

periodicList[f_, x0_, maxLength_, precision_] := Block[
  {s, r, list0, x, xx, i, wholeList},
  r[y_] := Round[y, 10.^(-precision)];
  s[_] = False;
  s[r@x0] = (i = 1);
  x = x0;
  wholeList = False;
  list0 = reapList[
    Sow[x0];
    While[True,
     xx = r@Sow[x = f@x];
     If[s[xx] =!= False, Break[]];
     s[xx] = i;
     If[i === maxLength, wholeList = True; Break[]];
     i++
    ]
   ];
  If[wholeList, list0, list0[[s[xx] + 1 ;;]]]
 ]

Pitifully small test suite

periodicList[#^2 - 1.754 &, -1.754, 10000, 3]
(* {1.32252, -0.00495143, -1.75398} *)

periodicList[# + 2 &, 1, 3, 3]
(* {1, 3, 5, 7} *)

I am too lazy to build up test suite, would appreciate help from lazyweb

Only worthwhile for expensive f

My idea is to make a new function g that builds up a list of values, then use FixedPoint on g with a SameTest that looks for repeat values. First, a helper (since MemberQ wasn't checking for numerical equality):

NMemberQ[list_List, num_?NumericQ, thres_: $MachineEpsilon] := 
  If[Select[list, Abs[# - num] < thres &, 1] == {}, False, True];

Then the main function:

Options[FixedPeriodList] = {Threshold -> $MachineEpsilon};

FixedPeriodList[func_, init_, max_: 1000, opts : OptionsPattern[]] := 
  Module[{threshold, g},
   threshold = OptionValue[Threshold];
   g[list_] := Append[list, func[list[[-1]]]];
   Return[FixedPoint[g, {init}, max, SameTest -> (NMemberQ[#1, #2[[-1]], threshold] &)]]
  ];

Testing it out:

f[z_] := z^2 + x;
x=-1.1; (* a 2-cycle *)
FixedPeriodList[f,1,100]

(* {1, -0.1, -1.09, 0.0881, -1.09224, ... , 0.091608, -1.09161, 0.091608, -1.09161} *)
(* ^^ Length=71 ^^ *)

Here's an analog to FixedPoint that just returns the cycle:

Options[FixedPeriod] = {Threshold -> $MachineEpsilon};

FixedPeriod[func_, init_, max_: 1000, opts : OptionsPattern[]] := 
  Module[{threshold, g, fplist, firstrep},
   threshold = OptionValue[Threshold];
   fplist = FixedPeriodList[func, init, max, Threshold -> threshold];
   firstrep = FirstPosition[fplist, n_ /; Abs[n - fplist[[-1]]] < threshold][[1]];
   Return[fplist[[firstrep + 1 ;;]]];
  ];

Testing:

FixedPeriod[f, 1, 1000]
(* {0.091608, -1.09161} *)

x = -1.3;
FixedPeriod[f, 1, 1000]
(* {0.389019, -1.14866, 0.0194303, -1.29962} *)

Feedback & improvements welcome!

As is evident from the plot in the question, any function designed to avoid unnecessarily computing data points should focus on 1, 2, and 4 cycle solutions, because higher cycles have extremely narrow bands in x. This can be accomplished by

nodups[x_] := Module[{t1, t2}, t1 = Nest[#^2 + x &, x, 999]; 
    t2 = DeleteDuplicatesBy[NestList[#^2 + x &, t1, 5], Round[#, .00000000001] &]; 
    If[Length@t2 < 5, t2, DeleteDuplicates@NestList[#^2 + x &, t1, 100]]]

Then,

ListPlot[Catenate[Table[{x, #} & /@ nodups[x], {x, -2, .25, .0001}]], 
    PlotStyle -> PointSize[Tiny], PlotRange -> All]

yields the same plot shown in the question but in 9.7 sec on my PC, as opposed to 10.5 sec for the code in the question. Incidentally, replacing DeleteDuplicatesBy with DeleteDuplicates is less effective at preventing the computation of duplicate values but reduces run time to 9.4 sec.

In general, we should not expect to reduce run time much by making the computation of data points more efficient, because the majority of the run time for this problem goes to plotting points. (See 107640 for more discussion of timing.) Moreover, code developed by users is unlikely to be anywhere near as efficient as Mathematica routines, like NestList and DeleteDuplicates.

Addendum

A plot of the number of data points vs x may be interesting.

ListLogPlot[Table[Length@nodups[x], {x, -2, .25, .0001}], 
    PlotStyle -> PointSize[Tiny], DataRange -> {-2, .25}]

As explained above, nodups works well for 1, 2, and 4 cycle portions of the return map. A perhaps unexpected feature of the plot is the spikes at bifurcations between 1, 2, and 4 cycles. Evidently, the return map has not reached steady state there even after 1000 iterations.

Second Addendum

FixedPointList also can be replaced by NestWhileList, which is more elegant than nodup but a bit slower.

ListPlot[Catenate[Table[{x, #} & /@ NestWhileList[#^2 + x &, Evaluate[Nest
    [#^2 + x &, x, 999]], Abs[#1 - #5] > 0.0000000001 &, 5, 100], {x, -2, .25, .0001}]], 
    PlotStyle -> PointSize[Tiny], PlotRange -> All]

The similar replacement of Nest by NestWhile is prohibitively slow.