Iteration and inverse iteration of Weierstrass elliptic ℘ functions

Here’s the code I used to generate the images of Weierstrass elliptic ℘ functions for my presentation Computational adventures in complex dynamics. For the most part, this should just run, though some images depend on definitions for previous images.

Weierstrass Nest for intro page

{g1,g2} = {5.739525,0};
f[z_] := N@WeierstrassP[z,{g1,g2}];
fInv[p_] := N@InverseWeierstrassP[p,{g1,g2}];
{w,w'} = N@WeierstrassHalfPeriods[{g1,g2}];
 bound = Abs[w + w'];
invImage[p_] :=Select[Flatten[
(({1,-1}fInv[p] + #)& /@ 
{0, 2w, 2w', 2w + 2w', 2w - 2w', -2w, 
-2w', -2w-2w', -2w + 2w'})
],
N[Abs[#]]≤bound&];
invImage[l_List] := Flatten[invImage /@ l];
sym1[z_] :=-z;
sym2[z_] := I z;

points0 = Nest[invImage, 
Random[] + Random[] I, 8]; // Timing

farOutPoints[n_, depth_:2] := Module[{points, points0, z1},
z1 = Nest[fInv, Random[] + Random[] I, 100];
points0 = Nest[invImage,z1, depth];
points =Flatten[Flatten[Join[
Table[points0 /. z_Complex -> fInv[z + 2 n w + 2 m w'],
{m,0,n}], 
Table[points0 /. z_Complex -> fInv[z + 2 n w' + 2 m w],
{m,0,n-1}]]]];
points = Join[points, sym1 /@ points];
points = Join[points, sym2 /@ points];
points = Join[points, Conjugate /@ points];
points];

morePoints = Flatten[{farOutPoints[1,6],farOutPoints[2,5],
farOutPoints[3,4], farOutPoints[4,3], farOutPoints[5,3],
Table[farOutPoints[n], {n,6,9}]}]; // Timing
joinedPoints = Join[points0,morePoints];
shiftedPoints = {#,#+2w, #+2w + 2w', # + 2w', # - 2w + 2w',
# - 2w, # - 2w - 2w', # - 2w', # + 2w - 2w'}& /@ joinedPoints;
evenMorePoints = invImage[Flatten[joinedPoints]];//Timing

goodPic = ListPlot[{Re[#],Im[#]}& /@ Join[joinedPoints, evenMorePoints,
Flatten[Table[farOutPoints[k],{k,10,20}]]],
PlotStyle -> AbsolutePointSize[.4], ImageSize->800,
AspectRatio -> Automatic, Axes->False]

r

r = Abs[w];
goodPic = ListPlot[{Re[#],Im[#]}& /@ Join[joinedPoints, evenMorePoints,
Flatten[Table[farOutPoints[k],{k,10,20}]]],
PlotStyle ->{AbsolutePointSize[.4]}, ImageSize->500,
AspectRatio -> Automatic, Axes->False, PlotRange -> {{-r,r},{-r,r}},
PlotRangePadding->None]

Export["wpTile.png", goodPic, ImageSize -> 500]

ImageAssemble[Table[goodImage,{3},{3}]]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "WeierstrassNestForIntro.png"}],goodPic,
ImageSize->{800,800}]

A simple contour plot

r = 3;
rePlot = ContourPlot[Re[WeierstrassP[x+I*y,{5.739525,0}]],{x,-r,r},{y,-r,r},
ContourShading->None, Contours->10, PlotPoints->50];
imPlot = ContourPlot[Im[WeierstrassP[x+I*y,{5.739525,0}]],{x,-r,r},{y,-r,r},
ContourShading->None, Contours->10, PlotPoints->50, ContourStyle->Blue];
contourPic = Show[{rePlot,imPlot}]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "wpContours.png"}],
contourPic,ImageSize->{800,800}]

Naive inverse iteration

{g1,g2} = {5.739525,0};
f[z_] := N@WeierstrassP[z,{g1,g2}];
fInv[p_] := N@InverseWeierstrassP[p,{g1,g2}];
simpleInvImage[p_] :=Flatten[
(({1,-1}fInv[p] + #)& /@ 
{0})
];
simpleInvImage[l_List] := Flatten[simpleInvImage /@ l];
pts = Nest[simpleInvImage, 
Random[] + Random[] I, 11];
badPic = ListPlot[{Re[#],Im[#]}& /@ pts, AspectRatio->Automatic, PlotStyle->Red]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "naive_inverse_iteration.png"}],badPic,
ImageSize->{800,800}]

Fitting them together

Needs["JLink`"];
InstallJava[];
AddToClassPath[FileNameJoin[{NotebookDirectory[],"JavaFiles"}]];
SetComplexClass["ORG.netlib.math.complex.Complex"];
WP = JavaNew["WeierstrassP"];

g3=0.0I;
{z0,g2} = {z,g2} /. FindRoot[{
WeierstrassP[z,{g2,0}]==z,
WeierstrassPPrime[z,{g2,0}]==0},
{z,1},{g2,5}];

n =10;
zMin = -1.4-1.4I;
zMax = 1.4+1.4I;
bail = 200;
tol = 0.001;
resRe = 800;
resIm = 800;
fps = {z0+0.0I};

t = AbsoluteTime[];
orbitData = WP@itCounts[zMin, zMax,g2+0.0I, g3+0.0I, n, fps, 
bail,tol, resRe,resIm];
AbsoluteTime[]-t

p = 10;
color[{cnt_,index_}] := {
1-(1-cnt/bail)^p, 1-(1-cnt/bail)^p,
1-0.5(1-cnt/bail)^p};
colors = Map[color, orbitData, {2}];

naiveEscapeCombo = Graphics[Raster[colors, {{Re[zMin],Im[zMin]},{Re[zMax],Im[zMax]}}],Frame->False,
AspectRatio -> (Im[zMax]-Im[zMin])/(Re[zMax]-Re[zMin]),Epilog->
{Red,Point[{Re[#],Im[#]}& /@ pts]}, PlotRangePadding->0];

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "NaiveInverseIterationOnEscapeTime.png"}],
naiveEscapeCombo,ImageSize->{800,800}]

Improving inverse iteration

{g1,g2} = {5.739525,0};
f[z_] := N@WeierstrassP[z,{g1,g2}];
fInv[p_] := N@InverseWeierstrassP[p,{g1,g2}];
{w,w'} = N@WeierstrassHalfPeriods[{g1,g2}];
 bound = Abs[w + w'];
invImage[p_] :=Select[Flatten[
(({1,-1}fInv[p] + #)& /@ 
{0, 2w, 2w', 2w + 2w', 2w - 2w', -2w, 
-2w', -2w-2w', -2w + 2w'})
],
N[Abs[#]]≤bound&];
invImage[l_List] := Flatten[invImage /@ l];
sym1[z_] :=-z;
sym2[z_] := I z;

points0 = Nest[invImage, 
Random[] + Random[] I, 10];

ListPlot[{Re[#],Im[#]}&/@points0, AspectRatio->Automatic]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "SlightlyImprovedInverseIteration.png"}],
%,ImageSize->{800,800}]

More  improvement

zInf1 = w+w';
infBack = Flatten@NestList[invImage, 
{zInf1, I zInf1, -zInf1, I^3 zInf1}, 2];
infBack2 = Flatten@NestList[invImage,
{2w, 2w + 2w', 2w', -2w + 2w',
 - 2w,  - 2w - 2w',  - 2w', 2w - 2w'}, 2];
infBacks = Join[infBack, infBack2];

pic = ListPlot[{Re[#],Im[#]}& /@ Join[joinedPoints, evenMorePoints,
Flatten[Table[farOutPoints[k],{k,10,20}]]],
PlotRange -> {{-1.6,1.6},{-1.6,1.6}}, PlotRangePadding->None,
PlotStyle -> AbsolutePointSize[.4], ImageSize->800,
Epilog -> {Join[{PointSize[.01], RGBColor[0,0,1]}, 
Point /@ ({Re[#],Im[#]}& /@  infBacks)]},
AspectRatio -> Automatic];

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "MoreImprovedInverseIteration.png"}],
pic,ImageSize->{800,800}]

(Subscript[g, 2])-space

radius = 9;
n = 10;
bail = 300;
tol = 0.0001;
res =300;
maxOrbitLength = 50;
t = AbsoluteTime[];
data = WP@g2Triangle[radius, n, bail, tol, res, maxOrbitLength];
AbsoluteTime[]-t
Clear[color];
p = 10; (* You can fiddle with p to change the contrast. *)
maxMaxValue = Max[Last /@ data];
color[{cnt_,orbitLength_, maxValue_}] := If[cnt<bail,
{0,0,0},{1,(1-Log[maxValue]/Log[maxMaxValue])^p,(0.7-0.5Log[maxValue]/Log[maxMaxValue])^p}];
MapIndexed[(color[#2[[1]]] = color[#])&,data];
pts = {Re[#],Im[#]}& /@ Flatten[Table[Table[
i*radius/res + j(1/2 + Sqrt[3]I/2)*radius/res,
{j,0,res-i}],{i,0,res} ]];
indices = Table[start = Sum[res+2-j,{j,1,i}]+1; fin = start+res-i;
Range[Length[pts]][[start;;fin]],{i,0,res}];
polys1 = Flatten[Table[Table[
{indices[[i,j]],indices[[i+1,j]],indices[[i,j+1]]},
{j,1,res+1-i}],{i,1,res}],1];
polys2 = Flatten[Table[Table[
{indices[[i,j]],indices[[i,j+1]],indices[[i-1,j+1]]},
{j,1,res+1-i}],{i,2,res}],1];
polys = Join[polys1,polys2];
vc = Map[color, polys,{2}];
initPolygon = GraphicsComplex[pts,
{Polygon[polys, VertexColors -> vc]}];
pastedPolygon = {initPolygon,GeometricTransformation[initPolygon,
ReflectionTransform[{-Sqrt[3]/2,1/2}]]};
Graphics[{pastedPolygon, 
GeometricTransformation[pastedPolygon,
RotationTransform[2Pi/3]],GeometricTransformation[pastedPolygon,
RotationTransform[4Pi/3]]}]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "g2Space.png"}],
%,ImageSize->{800,800}]

A  super-attractive orbit of period 2

g3=0.0I;
{z0,g2} = {z,g2} /. FindRoot[{
WeierstrassP[WeierstrassP[z,{g2,0}],{g2,0}]==z,
WeierstrassPPrime[z,{g2,0}]==0},
{z,1},{g2,3}]

n =10;
zMin = -2.2-2.2I;
zMax = 2.2+2.2I;
bail = 200;
tol = 0.001;
resRe = 800;
resIm = 800;
fps = {z0+0.0I};

t = AbsoluteTime[];
orbitData = WP@itCounts[zMin, zMax,g2+0.0I, g3+0.0I, n, fps, 
bail,tol, resRe,resIm];
AbsoluteTime[]-t

p = 10;
color[{cnt_,index_}] := {
1-(1-cnt/bail)^p, 1-(1-cnt/bail)^p,
1-0.5(1-cnt/bail)^p};
colors = Map[color, orbitData, {2}];

Graphics[Raster[colors, {{Re[zMin],Im[zMin]},{Re[zMax],Im[zMax]}}],Frame->False,
AspectRatio -> (Im[zMax]-Im[zMin])/(Re[zMax]-Re[zMin]), PlotRangePadding->0]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "p2Weierstrass.png"}],
%,ImageSize->{800,800}]

Inverse iteration for p2

g3=0.0I;
{z0,g2} = {z,g2} /. FindRoot[{
WeierstrassP[WeierstrassP[z,{g2,0}],{g2,0}]==z,
WeierstrassPPrime[z,{g2,0}]==0},
{z,1},{g2,3}]

Clear[w,w]

Flatten[Table[2i*w+2j*w',{i,-1,1},{j,-1,1}]]

invImage[p_] :=Select[Flatten[Map[({1,-1}fInv[p] + #)&,Table[2i*w+2j*w',{i,-1,1},{j,-1,1}],{2}]],
N[Abs[#]]≤bound&]

invImage[p_] :=Select[Flatten[
(({1,-1}fInv[p] + #)& /@ 
{0, 2w, 2w', 2w + 2w', 2w - 2w', -2w, 
-2w', -2w-2w', -2w + 2w'})
],
N[Abs[#]]≤bound&]

bound = 4;
res = 200;
f[z_] := N@WeierstrassP[z,{g2,0}];
fInv[p_] := N@InverseWeierstrassP[p,{g2,0}];
{w,w'} = N@WeierstrassHalfPeriods[{g2,0}];
 bound = Abs[w + w'];
invImage[p_] :=Select[Flatten[
(({1,-1}fInv[p] + #)& /@ 
{0, 2w, 2w', 2w + 2w', 2w - 2w', -2w, 
-2w', -2w-2w', -2w + 2w'})
],
N[Abs[#]]≤bound&];
invImage[l_List] := Flatten[invImage /@ l];
reducedInvImage[points_]:=Module[
{newPoints, croppedNewPoints, actuallyNewPoints,},
newPoints = invImage[points];
croppedNewPoints = Floor[res newPoints]/res;
croppedPointsSoFar = Floor[res pointsSoFar]/res;
actuallyNewPoints = Complement[croppedNewPoints,
croppedPointsSoFar];
positions = Flatten[Position[croppedNewPoints,#]& /@ 
actuallyNewPoints];
newPoints = newPoints[[positions]];
pointsSoFar=Union[newPoints,pointsSoFar];
newPoints
];

z1=Nest[fInv,Random[]+Random[] I,100];
(pointsSoFar=lastPoints=Nest[invImage,z1,2];
lastPoints=FixedPoint[reducedInvImage,lastPoints];) // Timing

points0 = pointsSoFar;
farOutPoints[n_, depth_:2] := Module[{points, points0, z1},
z1 = Nest[fInv, Random[] + Random[] I, 100];
points0 = Nest[invImage,z1, depth];
points =Flatten[Flatten[Join[
Table[points0 /. z_Complex -> fInv[z + 2 n w + 2 m w'],
{m,0,n}], 
Table[points0 /. z_Complex -> fInv[z + 2 n w' + 2 m w],
{m,0,n-1}]]]];
points = Join[points, Conjugate /@ points];
points];

morePoints = Flatten[{farOutPoints[1,6],farOutPoints[2,5],
farOutPoints[3,4], farOutPoints[4,3], farOutPoints[5,3],
Table[farOutPoints[n], {n,6,9}]}]; // Timing
joinedPoints = Join[points0,morePoints];
shiftedPoints = {#,#+2w, #+2w + 2w', # + 2w', # - 2w + 2w',
# - 2w, # - 2w - 2w', # - 2w', # + 2w - 2w'}& /@ joinedPoints;
evenMorePoints = invImage[Flatten[joinedPoints]];//Timing

allPoints=Join[joinedPoints, evenMorePoints,
Flatten[Table[farOutPoints[k],{k,10,20}]]];
allPoints = Flatten[Table[(I^k)allPoints,{k,0,3}]];

allPoints[[1]]

ListPlot[{Re[#],Im[#]}& /@ allPoints,
PlotStyle -> AbsolutePointSize[.4], ImageSize->800,
AspectRatio -> Automatic, Axes->True, PlotRangePadding->None,
PlotRange -> {{-2,2},{-2,2}}]

Export[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images", "p2InverseIteration.png"}],
%,ImageSize->{800,800}]

ListPlot[{Re[#],Im[#]}& /@ pointsSoFar,
PlotStyle -> AbsolutePointSize[.4], AspectRatio -> Automatic]

sym1[z_] :=-z;
sym2[z_] := I z;

points0 = Nest[invImage, 
Random[] + Random[] I, 12]; // Timing

farOutPoints[n_, depth_:2] := Module[{points, points0, z1},
z1 = Nest[fInv, Random[] + Random[] I, 100];
points0 = Nest[invImage,z1, depth];
points =Flatten[Flatten[Join[
Table[points0 /. z_Complex -> fInv[z + 2 n w + 2 m w'],
{m,0,n}], 
Table[points0 /. z_Complex -> fInv[z + 2 n w' + 2 m w],
{m,0,n-1}]]]];
points = Join[points, sym1 /@ points];
points = Join[points, sym2 /@ points];
points = Join[points, Conjugate /@ points];
points];

morePoints = Flatten[{farOutPoints[1,6],farOutPoints[2,5],
farOutPoints[3,4], farOutPoints[4,3], farOutPoints[5,3],
Table[farOutPoints[n], {n,6,9}]}]; // Timing
joinedPoints = Join[points0,morePoints];
shiftedPoints = {#,#+2w, #+2w + 2w', # + 2w', # - 2w + 2w',
# - 2w, # - 2w - 2w', # - 2w', # + 2w - 2w'}& /@ joinedPoints;
evenMorePoints = invImage[Flatten[joinedPoints]];//Timing

ListPlot[{Re[#],Im[#]}& /@ Join[joinedPoints, evenMorePoints,
Flatten[Table[farOutPoints[k],{k,10,20}]]],
PlotStyle -> AbsolutePointSize[.4], ImageSize->800,
AspectRatio -> Automatic, Axes->False]