Iteration of complex cubics

This notebook provides some functions to investigate the complex iteration of functions of the form

(Subscript[f, a,b])(z)=(z^3)-3(a^2)z+b.

Why?

We first claim that any cubic is conjugate to some (Subscript[f, a,b]). Thus, we are essentially studying the family of cubics. To see this suppose that

g(z)=(Subscript[c, 3])(z^3)+(Subscript[c, 2])(z^2)+(Subscript[c, 1])z+(Subscript[c, 0])

represents an arbitrary cubic. We need to find a function φ(z)=m z+d such that φ(g(z))=(Subscript[f, a,b])(φ(z)). We can set up the equations we need like so:

phi[z_]=m*z+d;
f[a_,b_][z_]= z^3-3a^2*z+b;
g[z_] = c3*z^3+c2*z^2+c1*z+c0;
eqs=Thread[CoefficientList[phi[g[z]]-f[a,b][phi[z]],z]=={0,0,0,0}]

And, now, let’s try to solve them

Solve[eqs,{a,b,d,m}]

I would interpret this to mean that φ(z)=(Sqrt[(Subscript[c, 3])])z+(Subscript[c, 2])/(3(Sqrt[(Subscript[c, 3])])) conjugates (Subscript[f, a,b]) to g when a and b are given by the last solution above. Evidently, the correspondence is not one-to-one, as in the quadratic case for (z^2)+c.

Investigating the parameter space

The initialization cells at the bottom of the notebook defines several functions for investigating the parameter space for (Subscript[f, a,b]).  The function cubicJuliaPic simply generates the Julia set.  It should be called like so:

cubicJuliaPic[a_?NumericQ, b_?NumericQ, bail_Integer, resolution_Integer]

cubicJuliaPic[-1/Sqrt[2],0,100,600]

There’s a dynamic version that accepts the same arguments but includes a locator allowing you to visualize the orbit starting from where the user clicks.

dynamicCubicJuliaPic[-1/Sqrt[2],0,100,600]

Now, the critical points of (Subscript[f, a,b])(z)=(z^3)-3(a^2)z+b are ±a.  To study the parameter space, we should partition the four dimensional space (ℂ^2) into the sets of points where

a and -a both converge to the same finite orbit

a and -a both converge but to different finite orbits

a escapes to ∞ but -a converges to a finite orbit

-a escapes to ∞ but a converges to a finite orbit

a and -a both escape to ∞

Black means that the classification was unsuccessful.  There are two functions to study the parameter space - one to generate a slice with a fixed and the other to generate a slice with b fixed.

parameterPic["a", a_?NumericQ, bMin_Complex, bMax_Complex, res_Real, graphicsOpts:OptionsPattern[]]
parameterPic["b", b_?NumericQ, aMin_Complex, aMax_Complex, res_Real, graphicsOpts:OptionsPattern[]]

For example:

parameterPic["a", 1/Sqrt[2], -2-I, 2+I, 0.005, 
Frame->True, PlotRangePadding -> None]

As with the Julia set command, there is a dynamic version that allows you to explore how Julia sets vary as the free parameter changes.

dynamicParameterPic["a", 1/Sqrt[2], -2-2I, 2+2I, 0.01]

dynamicParameterPic["a", ((I (Sqrt[59]))/10), -1-I, 1+I, 0.002]

dynamicParameterPic["b", 1/Sqrt[2], -2-I, 2+I, 0.01]

Here’s are some pictures of grids of slices to, hopefully, illustrate the 4D set.  These take about a half minute on my machine.

t=AbsoluteTime[];
pics = Table[parameterPic["a", ax+ay*I, -1.5-1.5I, 1.5+1.5I,0.02],
{ay,1,0,-0.1},{ax,0,1,0.1}];//AbsoluteTiming
corner = Graphics[{Line[{{-1,-1},{2,2}}],
Text["re(a)",{1,0}],
Text["im(a)",{0,1}]},
PlotRange->{{0,1},{0,1}},
PlotRangePadding -> 1];
GraphicsGrid[Join[MapIndexed[Prepend[#,Chop[N[1+((1-#2[[1]])/(Length[pics]-1))]]]&,pics],
{Flatten[{corner,Chop[0.1(#-1)]& /@ Range[Length[First[pics]]]}]}]]

AbsoluteTime[]-t

t = AbsoluteTime[];
pics = Table[parameterPic["b", bx+by*I, -1.5-1.5I, 1.5+1.5I,0.02],
{by,1,0,-0.1},{bx,0,1,0.1}];//AbsoluteTiming
corner = Graphics[{Line[{{-1,-1},{2,2}}],
Text["re(b)",{1,0}],
Text["im(b)",{0,1}]},
PlotRange->{{0,1},{0,1}},
PlotRangePadding -> 1];
GraphicsGrid[Join[MapIndexed[Prepend[#,Chop[N[1+((1-#2[[1]])/(Length[pics]-1))]]]&,pics],
{Flatten[{corner,Chop[0.1(#-1)]& /@ Range[Length[First[pics]]]}]}]]

AbsoluteTime[]-t

Initialization

cubicCriticalEscapeColor = 
  Compile[{{a, _Complex}, {b, _Complex}, {bail, _Integer}},
   Module[{cnt, n1, n2, z, orb1, orb2, escapeRadius, escapeFlag1, 
     escapeFlag2, 
     finiteFlag1, finiteFlag2, maxOrbitDetectionLength,
     orbRep1, orbRep2, period1, period2, tolerance},
    
    tolerance = 10.0^(-8);
    maxOrbitDetectionLength = 10;
    escapeRadius = Max[Abs[b], Sqrt[Abs[9 a^2] + 2]];
    
    (* Test z=a *)
    z = a;
    cnt = 1;
    escapeFlag1 = False;
    finiteFlag1 = False;
    period1 = 0; orbRep1 = 0.0 + 0.0 I;
    orb1 = Table[200.0 + 300.0 I, {maxOrbitDetectionLength}];
    While[++cnt < bail && escapeFlag1 === False && 
      finiteFlag1 === False,
     z = z^3 - 3 a^2 z + b;
     If[Abs[z] > escapeRadius,
      escapeFlag1 = True];
     If[Min[Abs[z - orb1]] < tolerance,
      finiteFlag1 = True;
      orbRep1 = z;
      z = z^3 - 3 a^2 z + b;
      period1 = 1;
      While[Abs[z - orbRep1] >= 100 tolerance,
       z = z^3 - 3 a^2 z + b;
       period1 = period1 + 1];
      ];
     orb1[[Mod[cnt, maxOrbitDetectionLength, 1]]] = z;
     ];
    n1 = cnt;
    
    (* Test z=-a *)
    z = -a;
    cnt = 1;
    escapeFlag2 = False;
    finiteFlag2 = False;
    period2 = 0; orbRep2 = 0.0 + 0.0 I;
    orb2 = Table[200.0 + 300.0 I, {maxOrbitDetectionLength}];
    While[++cnt < bail && escapeFlag2 === False && 
      finiteFlag2 === False,
     z = z^3 - 3 a^2 z + b;
     If[Abs[z] > escapeRadius,
      escapeFlag2 = True];
     If[Min[Abs[z - orb2]] < tolerance,
      finiteFlag2 = True;
      orbRep2 = z;
      z = z^3 - 3 a^2 z + b;
      period2 = 1;
      While[Abs[z - orbRep2] >= 100 tolerance,
       z = z^3 - 3 a^2 z + b;
       period2 = period2 + 1];
      ];
     orb2[[Mod[cnt, maxOrbitDetectionLength, 1]]] = z
     ];
    n2 = cnt;
    
    (* Color *)
    Which[
     escapeFlag1 === True (* a escapes *),
     Which[
      escapeFlag2 === True (* -a also escapes *),
      (* light shade *)
      {(1 - n1/bail)^20.0, (1. - 
          n1/(2. bail) - n2/(2 bail))^20, (1. - 2/bail)^20},
      finiteFlag2 === True (* -a is attracted to a finite orbit *),
      {(1.0 - n1/bail), 0.0, 0.5 (1.0 - n2/bail)},
      True, (* n2 === bail *)
      {(1.0 - n1/bail), 0.0, 0.0}],
     finiteFlag1 === True (* a is attracted to a finite orbit *),
     Which[
      escapeFlag2 === True (* -a escapes *),
      (* light shade *)
      {0.5 (1.0 - n1/bail), 
       0.0, (1.0 - n2/bail)},
      finiteFlag2 === True (* -a is attracted to a finite orbit *),
      If[period1 > 0 && period1 === period2 &&
        
        Min[Abs[orbRep1 - 
            NestList[#^3 - 3 a^2 # + b &, orbRep2, period2]]] < 
         100 tolerance &&
        
        Min[Abs[orbRep2 - 
            NestList[#^3 - 3 a^2 # + b &, orbRep1, period1]]] < 
         100 tolerance,
       {(n1/bail)^0.5, 0.0, (n2/bail)^0.5},
       {(n1/bail)^2, ((n1 + n2)/(2 bail))^0.25, (n2/bail)^2}],
      True, (* n2 === bail *)
      {0.5 (1.0 - n1/bail), 0.0, 0.0}],
     True (* n1 === bail *),
     Which[
      escapeFlag2 === True (* -a escapes *),
      (* light shade *)
      {0.0, 0.0, (1.0 - n2/bail)},
      finiteFlag2 === True (* -a is attracted to a finite orbit *),
      {0.0, 0.0, 0.5 (1.0 - n2/bail)},
      True, (* n2 === bail *)
      {0.0, 0.0, 0.0}]
     ]
    ],
   RuntimeAttributes -> {Listable}, Parallelization -> True,
   RuntimeOptions -> "Speed", CompilationTarget -> "MVM"
   ];
parameterPic["a", a_?NumericQ, bMin_Complex, bMax_Complex, res_Real, graphicsOpts : OptionsPattern[]] :=
  Graphics[Raster[cubicCriticalEscapeColor[
     a, Table[x + I*y, {y, Im[bMin], Im[bMax], res}, {x, Re[bMin], Re[bMax], res}], 1000],
    {Re[#], Im[#]} & /@ {bMin, bMax}], graphicsOpts];
parameterPic["b", b_?NumericQ, aMin_Complex, aMax_Complex, res_Real, graphicsOpts : OptionsPattern[]] :=
  Graphics[Raster[cubicCriticalEscapeColor[
     Table[x + I*y, {y, Im[aMin], Im[aMax], res}, {x, Re[aMin], Re[aMax], res}], b, 1000],
    {Re[#], Im[#]} & /@ {aMin, aMax}], graphicsOpts];

quickPts[a_, b_] := Block[{y, z, invfuncs, invImage, invfunc1,invfunc2,invfunc3},
	invfuncs = Quiet[z /.Solve[y==z^3 - 3a^2 z + b, z], Solve::ratnz];
	invfunc1 = Function[y,Evaluate[invfuncs[[1]]]];
	invfunc2 = Function[y,Evaluate[invfuncs[[2]]]];
	invfunc3 = Function[y,Evaluate[invfuncs[[3]]]];
	invImage[pts_] := Flatten[{invfunc1/@pts, invfunc2/@pts, invfunc3/@pts}];
	Nest[invImage, {1.0},5]
];
colorFunction = If[#==1, RGBColor[0.0,0.0,0.0], ColorData["BlueGreenYellow"][Mod[#,1]]]&;
colorLookup = List @@@ Array[colorFunction[#/2047.0]&, 2048, 0];
cubicJuliaColor = With[{colorLookup = colorLookup},
	Compile[{{a, _Complex}, {b, _Complex}, {z0, _Complex}, {bail,_Integer}},
	Module[{z, cnt, escapeRadius},
		escapeRadius = Max[Abs[b], Sqrt[Abs[9(a^2)]+2]];
		z = z0;
		cnt = 0;
		While[Abs[z] < escapeRadius && ++cnt < bail,
			z = z^3 - 3a^2 z + b
		];
		colorLookup[[Round[2047*(cnt/bail)^(0.25)]+1]]
	],
	RuntimeAttributes->{Listable},Parallelization -> True,
	RuntimeOptions->"Speed", CompilationTarget -> "MVM"
]];
cubicJuliaPic[a_?NumericQ, b_?NumericQ, bail_Integer, resolution_Integer] := 
	Module[{xmin,xmax, ymin,ymax, xStep,yStep,step, colors},
	{{xmin,xmax},{ymin,ymax}} = PlotRange[ListPlot[{Re[#],Im[#]}& /@ quickPts[a,b]]];
	{xmin,xmax} = (xmax+xmin)/2 -{1,-1} 2/3(xmax-xmin);
	{ymin,ymax} = (ymax+ymin)/2 -{1,-1} 2/3(ymax-ymin);
	xStep = N[(xmax - xmin)/resolution];
	yStep = N[(ymax - ymin)/resolution];
	step = Max[xStep,yStep];
	colors = cubicJuliaColor[a, b,
		Table[x+I*y, {y,ymin,ymax,step}, {x,xmin,xmax,step}], bail];
	Graphics[Raster[colors, {{xmin,ymin},{xmax,ymax}}]]
];
dynamicCubicJuliaPic[a_?NumericQ,b_?NumericQ, bail_, resolution_Integer] := Module[
	{(* background, plotRange *)},
	dcjpBackground = cubicJuliaPic[a,b,bail,resolution];
	dcjpPlotRange = PlotRange[dcjpBackground];
	DynamicModule[{pt={Re[a],Im[a]}},
		LocatorPane[Dynamic[pt], Dynamic@Show[{dcjpBackground,
			ListPlot[{Re[#],Im[#]}& /@
				NestWhileList[((#)^3)-3*(a^2)#+b&,N[pt.{1,I}],Abs[#]<5&,1,100],
				Joined -> True, PlotStyle -> Opacity[0.9]]},
			PlotRange -> dcjpPlotRange, ImageSize -> {400,400}]
		]
	]
];
dynamicParameterPic["a", a_?NumericQ, bMin_Complex, bMax_Complex, res_Real] := Module[
	{parameterImage},
	parameterImage = parameterPic["a", a, bMin, bMax,res, Frame -> True,
		ImageSize -> 400, PlotRangePadding -> None];
	DynamicModule[{pt={0.0,0.0I}}, Row[{
	LocatorPane[Dynamic[pt],parameterImage,{Re[#],Im[#]}& /@ {bMin,bMax}],
		Dynamic@dynamicCubicJuliaPic[a, pt.{1,I}, 100, 400]
	}]]];
dynamicParameterPic["b", b_?NumericQ, aMin_Complex, aMax_Complex, res_Real] := Module[
	{parameterImage},
	parameterImage = parameterPic["b", b, aMin, aMax,res, Frame -> True,
		ImageSize -> 400, PlotRangePadding -> None];
	DynamicModule[{pt={0.0,0.0I}}, Row[{
	LocatorPane[Dynamic[pt], parameterImage, {Re[#],Im[#]}& /@ {aMin,aMax}],
		Dynamic@dynamicCubicJuliaPic[pt.{1,I},b, 100, 400]
	}]]];