This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Construct an approximate solution to the functional equation f(z)=α f(f(z /α)) where α=1/f(1) in the form f(z)=1+(∑_(k=1)^n)(Subscript[c, k])(z^(2k)) with f''(0)<0 [✶], [✶], [✶], [✶], [✶]. Calculate α (the second Feigenbaum constant) to 20 significant digits. Use the power series and argument reduction to plot the function f(x) in the interval [100, 100].

f[z_, cList_] := 1 + cList.z^(2 Range[Length[cList]]);

rGSeries[z_, cList_] := With[{α = 1/f[1, cList]}, f[z, cList] - 
         α f[f[z/α, cList], cList] + O[z]^(2 Length[cList] + 1)];

eqs[cList_?(VectorQ[#, NumericQ]&)] :=
Module[{δ = rGSeries[z, cList]}, 
      Table[Coefficient[δ, z, 2j], {j, Length[cList]}]]

An approximate value for (Subscript[c, 1])is -1.4.

rGSeries[z, {c1}]

Solve[CoefficientList[%, z] == 0, {c1}] // N

Adding now one (Subscript[c, k]) at a time, we can quickly calculate a 20-digit approximation for 1/f(0). To avoid large symbolic expressions of the truncated functional equation, we restrict the first argument of FindRoot to purely numeric (Subscript[c, k]).

cListT = {-3/2}; ε = 10^-6; n = 1; αOld = 0;
  
While[αNew = 1/f[1, cListT]; Abs[αNew - αOld] > 10^-12,
      n++; αOld = αNew; cList = Table[c[k], {k, n}];
      (* smallest coefficient plus some digits *)
      scale = Round[-Log[10, Abs[Last[cListT]]]] + 10; 
      (* new coefficient list *)
      cListT = Last /@ (* solve for the c_k *)
        FindRoot[eqs[cList] == 0, (* starting values *)
         Evaluate[MapIndexed[{c[#2[[1]]], (1 - ε) #1, (1 + ε) #1}&, 
                  AppendTo[cListT, Last[cListT]/10]]],
         (* use adaptive precision and goals *)         
         WorkingPrecision -> 2 scale, PrecisionGoal -> scale,
         AccuracyGoal -> {scale, Infinity}, Method -> Secant]]

N[αNew, 20]

This results agrees with the known value of α [✶], [✶].

Now we will plot the function. A direct use of the power series allows to plot the function for small x, say x≲3. By observing that in the interval [-1,1], the function f(x) is always contracting, we repeatedly apply the functional equation to bring the argument finally to a size where we can apply the power series f(x)=(α^n)(f^(((2^n))))(x/(α^n)) (here (f^((k))) stands for the kth iterate of f) [✶]. The function fCF implements a compiled version of this procedure.

(* machine number version inside the interval [-1, 1] *)
fUD[z_] = 1. - 1.5276329970363323 z^2 + 
0.1048151947874277 z^4 + 0.026705670524930787 z^6 - 
0.003527409660464297 z^8 + 0.00008160096594827505 z^10 + 
0.000025285084886512315 z^12 - 2.5563177536625283*^-6 z^14 - 
9.65122702290271*^-8 z^16 + 2.8193175723520713*^-8 z^18 - 
2.771441260107602*^-10 z^20 - 3.0292086423142963*^-10 z^22 + 
2.6739057855563045*^-11 z^24 + 9.838888060875235*^-13 z^26 - 
3.5838769501333333*^-13 z^28 + 2.063994985307743*^-14 z^30;

fCF = Compile[{z},
             Module[{α = -2.5029078750959130867, n, ζ},
                    If[Abs[z] <= 1., n = 0, (* reduce argument size *)
                       n = Ceiling[Log[-α, Abs[z]]]]; 
                    ζ = z/α^n; Do[ζ = #, {2^n}]; (α)^n ζ]]&[fUD[ζ]];

Here is a plot of f(x). It has a quite complicated structure [✶].

Plot[fCF[x], {x, -100, 100}, PlotPoints -> 200, PlotRange -> All]

We now could try to make a plot in the complex plane. Unfortunately, the last construction does not always converge because for some numbers z from the unit disk, the values f(z) become too large [✶]. Trying to use the recursion directly might not terminate. So we restrict the number of recursions in the following (to 40) and calculate the function over a domain where the recursions end quickly.

(* use larger value to obtain potentially more points *)
$RecursionLimit = 40 ;
α = 2.5029078750959130867`(* use machine arithmetic! *);
fR[z_] := -α fR[fR[z/α]] /; Abs[z] >= 1
fR[z_] := fUD[z] /; Abs[z] < 1

L = 16; pp = 240;
Off[$RecursionLimit::reclim];
data = Table[If[NumberQ[#], #, 
                (* no recursion termination *) C]&[fR[x + I y]],
             {y, -L, L, 2L/pp}, {x, -L, L, 2L/pp}];

Here is the absolute value of the Feigenbaum–Cvitanović function for the z-values where the recursion terminates within 40 iterations. It shows a very complicated behavior over the complex plane.

Show[Graphics3D[
ListPlot3D[ArcTan[Abs[data]] /. ArcTan[Abs[C]] -> -1,
           (* color cyclically with argument *)
           Map[Hue, (Arg[data] + Pi)/(2Pi), {2}], ClipFill -> None,
           PlotRange -> {0, Pi/2}, BoxRatios -> {1, 1, 0.3}, 
           Mesh -> False, DisplayFunction -> Identity,
           MeshRange -> L {{-1, 1}, {-1, 1}}]] /. 
           (* remove steep vertical polygons *)
           {p:Polygon[l_] :> If[Min[Last /@ l] < 0, {}, p],
            h_Hue :> SurfaceColor[h, h, 2.2]},
           DisplayFunction -> $DisplayFunction, Lighting -> True]

For the Feigenbaum attractor, see [✶].