8
$\begingroup$

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

bifurcation

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?

$\endgroup$
  • 2
    $\begingroup$ +1 just because "you've had fingers wagged at you" $\endgroup$ – Dr. belisarius Feb 19 '16 at 1:20
  • $\begingroup$ I probably should release that package soon, if you need this functionality... anyway, "MMa jocks scoff at Do loops!" is not necessarily true; sometimes loops are the clearest way to express an idea. $\endgroup$ – J. M. is away Feb 19 '16 at 2:45
  • $\begingroup$ @J.M. What package is that? Sounds intriguing $\endgroup$ – Chris K Feb 21 '16 at 1:37
  • $\begingroup$ @Chris, I wrote a package, but I'm not yet done with the docs. I'll post here once it's done. $\endgroup$ – J. M. is away Feb 21 '16 at 1:43
5
$\begingroup$

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}]

enter image description here

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.

$\endgroup$
  • $\begingroup$ Yes, the Answer you gave me on my first question about this accomplish most of the efficiency gain, but then I got interested in squeezing harder, because I need to do some of these calculations even if I don't plot the output. $\endgroup$ – Jerry Guern Feb 19 '16 at 5:11
2
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Thank you for this. Your coding is a few grade levels above mine, so I had to look up some things to understand what you are doing. I have two questions. Was there some reasons you wanted to use the existing FixedPoint function? And isn't it really slow to build up a long list through repeated use of Append? $\endgroup$ – Jerry Guern Feb 19 '16 at 3:43
  • $\begingroup$ It seemed like a good idea at the time, but testing against your diagram shows that it is in fact, much slower than @bbgodfrey's solution. $\endgroup$ – Chris K Feb 19 '16 at 13:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.