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

Question

I need to efficiently generate a bunch of bifurcation diagrams. A simple familiar example illustrates the efficiency problems I'm running into:

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

This method is a bit inefficient.

FixedPointList[] does nothing unless the series z->(z^2+x) becomes 1-periodic (to within machine precision) within 999 cycles. Some numbers like -.75 and -1.25 still don't quite reach their asymptotes after 999 cycles, which is why those bifurcations look soft on my diagram above. But then for some numbers like the neighborhoods of -1, 999 is ridiculous overkill. DeleteDuplicates cleans up the 100-fold redundancy I wish I hadn't had to generate in the first place.

All these problems would be solved if I had a function similiar-ish to FixedPointList that worked on periodic functions. But I don't know how to write a function that takes a function as an input. Here's my best effort below with a hard-wired function z->(z^2+x), but it's so sloooooow it's worse than the one-line code above, but I'd bet that's a problem with my coding. I've had fingers wagged at me for using Do Loops, but I don't know how to get by without one in this case.

My code below returns a 3-element result:

{List become periodic True/False, Element in list that would match the next element in the list, the result list}

FixedPeriodList[x_, maxLength_] := Module[{ret},
    ret = 0;
    ans = ConstantArray[0, maxLength];
    ans[[1]] = x;
    Do[   (*MMa jocks scoff at Do loops!*)
       ans[[i]] = ans[[i - 1]]^2 + x;
        If[ContainsAny[Part[ans, 1 ;; i - 1], {ans[[i]]}],
           ret = {True, Position[Part[ans, 1 ;; i - 1], ans[[i]]][[1, 1]],Part[ans, 1 ;; i - 1]};
           Break[]], {i, 2, maxLength}];
    If[ret == 0, ret = {False, 1, ans}];
    ret];

a=FixedPeriodList[-1.754, 10000]

Part[#[[3]], #[[2]] ;; Length[#[[3]]]] &@ a

(* {True, 34, {-1.754, 1.32252, -0.00495143, -1.75398, 1.32243, -0.00517891, -1.75397, 1.32242, -0.00520029, -1.75397, 1.32242, -0.00520234, -1.75397, 1.32242, -0.00520254, -1.75397, 1.32242, -0.00520256, -1.75397, 1.32242, -0.00520256, -1.75397, 1.32242, -0.00520256, -1.75397, 1.32242, -0.00520256, -1.75397, 1.32242, -0.00520256, -1.75397, 1.32242, -0.00520256, -1.75397, 1.32242, -0.00520256}}

{-1.75397, 1.32242, -0.00520256} *)

How can we create a FixedPeriodList that works more efficiently? And with the function specifiable like in FixedPointList?

Answer 1

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.

Answer 2

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

Answer 3

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!