Code for Electronic Supplement

The JuliaSet Package

(* :Name: JuliaSet` *)
(* :Author: Mark McClure, April 1998 *)  (* some additions by Stan Wagon, August, 2007 *)
(* :Copyright: Copyright Mark McClure, 1998 *)
(* :Package Version: 1.1 *)
(* :Mathematica Version: 2.2 *)
(* :Summary:
	Generates images of Julia sets using several
	variants of the inverse iteration algorithm. 
*)

BeginPackage["JuliaSet`"];


InitialIterations::usage = "InitialIterations is an option to JuliaModified that sets the number of times the inverse iteration method is run.";


ComplexTicks::usage = "ComplexTicks[s] turns the complex number (or list) s into specifications for complex tick marks.";

TraceQ::usage = "TraceQ is an option to JuliaModified that prints the iteration number.";
TraceStep::usage = "";
JuliaSet::usage = "JuliaSet is a package containing several
implementations of inverse iteration algorithms for the
generation of Julia sets.";

Julia::usage = "Julia[formula,var] generates the Julia set of 
the rational function defined by formula, using the modified 
inverse iteration algorithm.";
 
JuliaSimple::usage = "JuliaSimple[c,depth] generates 
2^depth points of the Julia set for z^2 + c using a simple 
deterministic algorithm.";
 	
JuliaModified::usage = "JuliaModified[c] generates the Julia set 
for z^2 + c using a modified determinstic algorithm.";

JuliaStochastic::usage = "JuliaStochastic[c,numpoints] generates 
numpoints points of the Julia set for z^2 + c using a stochastic
algorithm.";

FilledJuliaSet::usage = "FilledJuliaSet[c, numpoints] generates a filled Julia set for z^2 + c using an escape time algorithm.";

Resolution::usage = "Resolution is a positive integer valued 
Option of Julia and JuliaModified which controls the level of 
detail of the generated set.";

Bound::usage = "Bound is a positive real valued Option of 
Julia which specifies the maximum absolute vale of points
plotted in the Julia set of a rational function.";

ColorQ::usage = "ColorQ is an option to FilledJuliaSet that chooses between a colored or a gray color scheme.";

Iterations::usage = "Iterations is an option to FilledJuliaSet that sets the number of iterations to try.";


OneBias::usage = "OneBias is an Option of JuliaStochastic which 
determines the probability of choosing the positive inverse of 
z^2 + c.  OneBias must be a real number stricly between 0 and 1.";

JuliaSetOnlyQ::usage = "JuliaSetOnlyQ is an option to FilledJuliaSet that, when True, causes only the black Julia set to appear.";
 	
Begin["`Private`"];

Attributes[ComplexTicks] = Listable;

ComplexTicks[s_?NumericQ] := 
  {s, s I} /. Thread[{-1. I, 0. I, 1. I} ->  {-I, 0, I}];


(* fR is like FilterOptions *) 

fR[sym_, opts___] := Sequence @@ FilterRules[{opts}, First/@Options[sym]];

(* Set the Options *)

(* For JuliaSimple *)
JuliaSimpleOptions = Join[{AspectRatio -> Automatic, Axes -> False,
		PlotStyle -> {PointSize[0.0005], Black}}, Options[ListPlot]];
keywords = Complement[Union[First /@ JuliaSimpleOptions], 
	{DisplayFunction, DefaultFont}];
vals = keywords /.  JuliaSimpleOptions;
special = {DisplayFunction :> $DisplayFunction, DefaultFont :> $DefaultFont};
Options[JuliaSimple] = Union[Apply[Rule,Transpose[{keywords,vals}],{1}],
	special];

(* For JuliaModified *)
Options[JuliaStochastic] = Join[{OneBias -> .5}, Options[JuliaSimple]];
Options[JuliaModified] = Join[{Resolution -> 200, TraceQ->False,TraceStep->20, InitialIterations->10}, Options[JuliaSimple]];
Options[Julia]    = Join[{Bound -> 4}, Options[JuliaModified]];


(* The functions *)

(* Error Messages *)
JuliaSet::argu = "`1` called with `2`; `3` expected.";
JuliaSet::badNumber = "Number expected at position 1 in `1`.";
JuliaSet::badInt = "Positive Integer expected at position 2 in `1`.";
JuliaSet::badOpt = "Options expected as optional arguments in `1`.";
JuliaSet::optx = "Unknown Option `1`.";
Julia::badFormula = "Rational function of degree at least two 
expected at position 1 in `1`.";
Julia::noInverses = "Can't find the the inverse functions.";
OneBias::outOfRange = "OneBias must be stricly between 0 and 1. Using 0.5.";
Resolution::outOfRange = "Resolution must be a positive integer. Using 200.";
Bound::outOfRange = "Bound must be a positive number.  Using 4.";


JuliaSimple[c_, depth_, opts___] := Module[
	{invImage, points, valid = First /@ Options[JuliaSimple]},
	
	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		Flatten[{opts}]
	];
	
	invImage := 
		( N[{1,-1} Sqrt[#-c]]& /@ # // Flatten ) &;
	points = {Re[#], Im[#]}& /@ Nest[invImage, {1}, depth];
	ListPlot[points, 
		fR[ListPlot,opts], 
		fR[ListPlot,Sequence @@ Options[JuliaSimple]]
	]
] /;  
  ( NumberQ[N[c]] || Message[JuliaSet::badNumber, JuliaSimple, c] ) &&
	( (IntegerQ[depth] && depth > 0 ) ||
	 Message[JuliaSet::badInt, JuliaSimple, depth] ) &&
	( And @@ Map[OptionQ, {opts}] ||
	  	Message[JuliaSet::badOpt,JuliaSimple] )

(* Special cases for JuliaSimple *)
JuliaSimple[]   := Message[JuliaSet::argu, "JuliaSimple", 
	"0 arguments", "2 arguments are"]
JuliaSimple[c_] := Message[JuliaSet::argu, "JuliaSimple", 
	"1 argument",  "2 arguments are"]


JuliaModified[c_, opts___] := Module[
	{cc,invImage,reducedInvImage,pointsSoFar,res, 
	valid = First /@ Options[JuliaModified]},

	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		{opts}];
	{res, trQ, ninit,step} = {Resolution, TraceQ, InitialIterations, TraceStep} /. {opts} /. Options[JuliaModified];
  cc=ninit;
        If[!(IntegerQ[res] && res>0),
		{Message[Resolution::outOfRange];
		res = Resolution /. Options[JuliaModified]}];
	
	invImage := 
		( N[{1,-1} Floor[N[res Sqrt[#-c]]]/res]&  /@ # // Flatten )&;

	reducedInvImage[points_] := Module[{newPoints},
		newPoints = Complement[invImage[points], pointsSoFar];
		pointsSoFar = Join[newPoints, pointsSoFar];
		If[trQ && Divisible[cc,step], Print[StringForm["`` levels, `` points, `` new points",cc,Length[pointsSoFar],Length[newPoints]]]];
		newPoints];  (* Mark used Union here *) 
    pointsSoFar = Nest[invImage, {1}, ninit];
FixedPoint[(cc++;reducedInvImage[#])&, pointsSoFar];

If[trQ,Print[StringForm["`` levels and `` points",cc,Length[pointsSoFar]]]];
	ListPlot[{Re[#], Im[#]}& /@ pointsSoFar, 
		fR[ListPlot,opts], 
		fR[ListPlot,Sequence @@ Options[JuliaModified]]
	]] /;  
	(NumberQ[c] || Message[JuliaSet::badNumber, JuliaModified, c]) &&
    ( And @@ Map[OptionQ, {opts}] ||
	  	Message[JuliaSet::badOpt,JuliaModified] )

(* Special case for JuliaModified *)
JuliaModified[] := Message[JuliaSet::argu, "JuliaModified", 
	"0 arguments", "1 argument is"]
	
	

Julia[formula_, var_Symbol, opts___] := Module[
	{z0,inverses,inverseRules,funcs,image,
	reducedImage,pointsSoFar, res,bound,
	valid = First /@ Options[Julia]},

	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		{opts}
	];

	res = Resolution /. {opts} /. Options[Julia];
	If[!(IntegerQ[res] && res>0),
		{Message[Resolution::outOfRange];
		res = Resolution /. Options[Julia]}];
	bound = Bound /. {opts} /. Options[Julia];
	If[!(bound>0),
		{Message[Bound::outOfRange];
		bound = Bound /. Options[Julia]}];
	
	inverseRules = NSolve[#==formula,var];
	If[And @@ ListQ /@ inverseRules, If[Length[inverseRules] > 1,
	
		inverses = var /. inverseRules;
		funcs = Function[anInverse,N[Floor[N[anInverse res]]/res] &] 
			/@ inverses;
		image = Flatten[Through[funcs[#]]& /@ #,1]&;
		If[PolynomialQ[formula,var],
			reducedImage[points_] := Module[{newPoints},
				newPoints = Complement[image[points], pointsSoFar];
				pointsSoFar = (* Union *) Join [newPoints, pointsSoFar];
				newPoints],
			reducedImage[points_] := Module[{newPoints},
				newPoints = Complement[image[points], pointsSoFar];
				newPoints = Select[newPoints,(N[Abs[#]]<=bound)&];
				pointsSoFar = Union[newPoints, pointsSoFar];
				newPoints]
		];
		z0 = N[Nest[funcs[[1]][funcs[[2]][#]]&, N[Pi],10]];
		pointsSoFar = Select[Nest[image,{z0},3],(N[Abs[#]]<=bound)&];
		FixedPoint[reducedImage, pointsSoFar];
		ListPlot[{Re[#],Im[#]}& /@ pointsSoFar, 
			fR[ListPlot,opts], 
			fR[ListPlot,Sequence @@ Options[JuliaModified]]
		],
		Message[Julia::badFormula, Julia]], Message[Julia::noInverses] 
	]
] /; (NumberQ[N[formula /. var -> Pi]]                &&
	PolynomialQ[Numerator[Together[formula]], var]    &&
	PolynomialQ[Denominator[Together[formula]], var]) || 
	Message[Julia::badFormula, Julia]

(* Special cases for Julia *)
Julia[]   := Message[JuliaSet::argu, "Julia", 
	"0 arguments", "2 arguments are"]
Julia[c_] := Message[JuliaSet::argu, "Julia", 
	"1 argument",  "2 arguments are"]



JuliaStochastic[c_, numpoints_, opts___] := Module[
	{inv1,inv2,points,z0,oneBias,chooser, 
	valid = First /@ Options[JuliaStochastic]},
	
	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		{opts}
	];

	oneBias = OneBias /. {opts} /. Options[JuliaStochastic];
	If[!(0<oneBias && oneBias<1),
		{Message[OneBias::outOfRange];
		oneBias = OneBias /. Options[JuliaStochastic]}];
	
	inv1 = N[Sqrt[#-c]]&; inv2 = N[-Sqrt[#-c]]&;
	z0 = Nest[inv1[inv2[#]]&, 1, 5];
	points = ComposeList[Table[chooser = Random[];
		If[chooser<oneBias, inv1, inv2], 
		{numpoints}], z0];
	ListPlot[{Re[#],Im[#]}& /@ points,  
		fR[ListPlot,opts], 
		fR[ListPlot,Sequence @@ Options[JuliaStochastic]]
	]
] /;  
	( NumberQ[c] || Message[JuliaSet::badNumber, JuliaStochastic, c] ) &&
	  ( (IntegerQ[numpoints] && numpoints > 0 ) ||
	  	Message[JuliaSet::badInt, JuliaStochastic, numpoints] ) &&
	  ( And @@ Map[OptionQ, {opts}] ||
	  	Message[JuliaSet::badOpt,JuliaStochastic] )
	  	
(* Special cases for JuliaStochastic *)
JuliaStochastic[]   := Message[JuliaSet::argu, "JuliaStochastic", 
	"0 arguments", "2 arguments are"]
JuliaStochastic[c_] := Message[JuliaSet::argu, "JuliaStochastic", 
	"1 argument",  "2 arguments are"]

      
JuliaIterationsCom = Compile[{{z, _Complex}, {c, _Complex}, {iters, _Integer}}, 
   Module[{i, s}, i = 1; s = z; While[i < iters && Abs[s] < 2, i++; s = s^2 + c]; i]];

Options[FilledJuliaSet] = {Iterations->100, PlotPoints->100, JuliaSetOnlyQ->False, ColorQ->True};

FilledJuliaSet[c_, {x0_,x1_},{y0_,y1_}, opts___] :=  (
  {iters,pp, JSonlyQ, colQ} = {Iterations,PlotPoints,JuliaSetOnlyQ, ColorQ} /. {opts} /. Options[FilledJuliaSet];
  DensityPlot[
    JuliaIterationsCom[x + y I, c, iters], {x, x0, x1}, {y,y0,y1}, Evaluate[fR[DensityPlot, opts]],PlotPoints -> pp, PlotRange->All,AspectRatio->Automatic,
    ColorFunction -> If[colQ, (If[# == 1, Black, If[JSonlyQ,White,Hue[0.9 #^(3/4)]]] &),(If[# == 1, Black, GrayLevel[ #^0.25]] &)]])
      
      
End[];  (* End Private Context *)

(*Protect[JuliaSet, Julia, JuliaModified, JuliaSimple, 
	JuliaStochastic, Bound, Resolution, OneBias]*)

EndPackage[];

Complex Dynamics: 
Julia Sets and the Mandelbrot Set

by Mark McClure, University of North Carolina at Asheville

Code for Chapter Opening Art

im1=FilledJuliaSet[0.32+0.043 I,{-1.5,1.5},{-1.3,1.3},FrameTicks->{Range[-1.5, 1.5, 1.5], ComplexTicks[Range[-1,1,.5]], None, None}, ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&),Frame->False];

im2=FilledJuliaSet[cRabbit=-0.123+0.745 I,{-1.4,1.4},{-1.3,1.3},PlotPoints->100,AspectRatio->Automatic,Frame->False, ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&)];

JuliaIterationsCompiledQM=Compile[{{z,_Complex},{r,_Complex},{iters,_Integer}},Module[{i,s},i=1;s=z;While[i<iters&&Abs[s]<2,i++;s=r s (1-s)];i]];
im3=DensityPlot[JuliaIterationsCompiledQM[x+y I,3,100],{x,-0.1,1.1},{y,-0.4,0.4},PlotPoints->100,ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&), AspectRatio->Automatic, Frame->False];

JuliaIterationsCompiled=Compile[{{z,_Complex},{c,_Complex},{iters,_Integer}},Module[{i,s},i=1;s=z;While[i<iters&&Abs[s]<2,i++;s=(s^2)+c];i]];MandelbrotIterations[c_, bail_] := JuliaIterationsCompiled[0,c,bail]

im4=DensityPlot[MandelbrotIterations[x+I y,100],{x,-2,0.6},{y,-1.3,1.3}, AspectRatio->Automatic,PlotPoints->100, Frame->False,ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&)];

GraphicsGrid[{{im2, im1},{im3, im4}}]

.  Complex Dynamics

f[z_] := (z^2);
(Subscript[z, 0])=0.9 + 0.15 I;
orbit=NestList[f,(Subscript[z, 0]),10]

ComplexTicks[s_?NumericQ]:={s,s I}/.Thread[{-1. I,0. I,1. I}->{-I,0,I}];ListPlot[{Re[#],Im[#]}& /@ orbit, Frame->True,FrameTicks->{Automatic, ComplexTicks[Range[0,.4,.1]],None,None}]

ListLinePlot[{Re[#],Im[#]}& /@ orbit, Frame->True,FrameTicks->{Automatic, complexTicks[Range[0,.4,.1]],None,None}]

α=((2π k)/(2^m))  implies  (f^((m)))((Subscript[z, 0]))=((((e^((2π i k)/(2^m)))))^(2^m))=(e^(2π k i))=1.

(Subscript[z, 0])=(E^((2 π I)/15));
orbit=NestList[f,(Subscript[z, 0]),12]

ListPlot[({Re[#],Im[#]}&)/@orbit,AspectRatio->Automatic,Frame->True,FrameTicks->{Automatic, complexTicks[Range[-.2,1,.2]],None,None}]

(Subscript[z, 0])=(E^I);
orbit=NestList[f,(Subscript[z, 0]),800];
ListPlot[({Re[#],Im[#]}&)/@orbit,AspectRatio->Automatic]

ListLinePlot[({Re[#],Im[#]}&)/@orbit,AspectRatio->Automatic]

ComplexTrajectory[c_,z_,{n0_,n1_},opts___]:=Graphics[data=({Re[#1],Im[#1]}&)/@NestList[((#1)^2)+c&,Nest[((#1)^2)+c&,N[z],n0],n1];{Thickness[0.003],Line[data],PointSize[0.015],Point[data],PointSize[Large],Point[data[[1]]]},opts,Frame->True,PlotRange->All]

Clear[c,n0,n1,z];n0[_] =0;n1[_]=25;
c[1]=c[2]=-0.12256117+0.74486177 I;z[1]=1.2-0.55 I;n1[1]=25;pr[1]={{-1,1.4},{-1,1}};
z[2]=-0.32+0.5 I;n1[2]=60;pr[2]={{-0.9,0.2},{-0.2,0.9}};
c[3]=c[4]=0.32+0.043 I;z[3]=-0.36-0.64 I;n1[3]=150;pr[3]={{-0.8,1},{-0.8,1}};
z[4]=z[3];n0[4]=300;n1[4]=20;pr[4]={{-0.8,1},{-0.8,1}};
c[5]=c[6]=-0.390541-0.586788 I;z[5]=0.3;pr[5]={{-0.8,0.4},{-0.8,0.3}};
z[6]=-0.2-0.3I;n1[6]=25;pr[6]={{-0.8,0.4},{-0.8,0.3}};
GraphicsGrid[Partition[Table[ComplexTrajectory[c[i],z[i],{n0[i],n1[i]},PlotRange->pr[i]],{i,6}],2]]

.  Julia Sets and Inverse Iteration

Simple inverse iteration

c=0.;
invImage[complexPoints_]:=Flatten[({1,-1} (Sqrt[#-c])&)/@complexPoints];
invImage[{1}]

invImage[%]

depth =10; 
points ={Re[#], Im[#]}&/@ Nest[invImage,{1},depth]; 
Graphics[{PointSize[Tiny],Point[points]}]

JuliaSimple[-0.123+0.745 I,14]

Improving the inverse iteration algorithm

JSet=JuliaModified[-0.123+0.745 I,Resolution->300]

histData=Tally[Cases[JuliaSimple[-0.123+0.745 I,14],Point[z_]:>z,∞][[1]],Norm[#1-#2]<0.01&];
Graphics3D[{Thin,(Line[{Flatten[{#[[1]],0}],Flatten[#]}]&)/@histData},BoxRatios->{1,1,(1/3)},PlotRange->All, ViewPoint->{1,-2.3, 1}]

GraphicsRow[{JuliaSimple[0.68 I,12],JuliaModified[0.68 I]}]

GraphicsRow[{JuliaSimple[0.24,12],JuliaSimple[0.26,12]}]

GraphicsRow[{JuliaModified[0.24, Resolution->400],JuliaModified[0.26,Resolution->400]}]

c=-0.123+0.745 I;res=300;
invImage[points_]:=Flatten[((Floor[N[{1,-1} res (Sqrt[#-c])]])/res)&/@points];

reducedInvImage[points_]:=Module[{newPoints},
newPoints=Complement[invImage[points],pointsSoFar];
pointsSoFar=Join[newPoints,pointsSoFar];
newPoints];

Take[(pointsSoFar=Nest[invImage,{1},10]),6]

FixedPoint[reducedInvImage,pointsSoFar];

Graphics[{PointSize[Tiny], Point[{Re[#],Im[#]}&/@pointsSoFar]}]

GraphicsRow[{JuliaModified[-0.77+0.22 I,Resolution->50],
JuliaModified[-0.77+0.22 I,Resolution->300]}]

JuliaModified[0.68 I,TraceQ->True,Resolution->300];

Generalizing the algorithm

f[z_]=(z^3)+z+0.6 I;res=200;
inverses=z/.NSolve[f[z]==#,z];
funcs=Function[anInverse,N[((Floor[anInverse res])/res)]&]/@inverses;

invImage[points_]:=Flatten[(Through[funcs[#]]&)/@points,1];
invImage[{1}]

invImage[%]

reducedImage[points_]:=Module[{newPoints},
newPoints=Complement[invImage[points],pointsSoFar];
pointsSoFar=Join[newPoints,pointsSoFar];newPoints];
pointsSoFar=Nest[invImage,{1.},10];
FixedPoint[reducedImage,pointsSoFar];
Graphics[{PointSize[Tiny],Point[({Re[#],Im[#]}&)/@pointsSoFar]}]

GraphicsRow[{Julia[(z^3)-I,z],Julia[(1/((z^2)-1)),z]}]

bound=4;
reducedImage[points_] := Module[{newPoints},
newPoints =Complement[image[points], pointsSoFar];
newPoints =Select[newPoints,(N[Abs[#]]≤bound)&];
pointsSoFar =Join[newPoints, pointsSoFar];
newPoints];

Julia[((2 z)/3)+(1/(3 (z^2))),z,Bound->12]

.  Escape Time Algorithms and the Mandelbrot Set

Escape time algorithms for Julia sets

JuliaIterations[z_,c_,nmax_:100]:=Length[FixedPointList[((#)^2)+c&,z,nmax,SameTest->(Abs[#]>2.0&)]];

DensityPlot[JuliaIterations[x+y I,-0.123+0.745 I],{x,-1.4,1.4},{y,-1.3,1.3},PlotPoints->100,Mesh->False,AspectRatio->Automatic,Frame->False,ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&)]

JuliaIterationsCompiled=Compile[{{z,_Complex},{c,_Complex},{iters,_Integer}},Module[{i,s},i=1;s=z;While[i<iters&&Abs[s]<2,i++;s=(s^2)+c];i]];

DensityPlot[JuliaIterationsCompiled[x+y I,-0.77+0.22 I,100],{x,-1.6,1.6},{y,-1.2,1.2},PlotPoints->120,ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&), FrameTicks->{Automatic, ComplexTicks[Range[-1,1,.5]], None, None},PlotRange->All, AspectRatio->Automatic]

FilledJuliaSet[0.32+0.043 I,{-1.5,1.5},{-1.3,1.3},FrameTicks->{Range[-1.5, 1.5, 1.5], ComplexTicks[Range[-1,1,.5]], None, None}]

FilledJuliaSet[-0.390541-0.586788 I,{-1.6,1.6},{-1.2,1.2}, 
FrameTicks->{Range[-1.5, 1.5, 1.5], ComplexTicks[Range[-1,1,.5]], None, None}]

JuliaIterationsCompiledQM=Compile[{{z,_Complex},{r,_Complex},{iters,_Integer}},Module[{i,s},i=1;s=z;While[i<iters&&Abs[s]<2,i++;s=r s (1-s)];i]];
DensityPlot[JuliaIterationsCompiledQM[x+y I,3,100],{x,-0.1,1.1},{y,-0.4,0.4},PlotPoints->100,ColorFunction->(If[#1==1,Black,Hue[0.9 ((#1)^(3/4))]]&),FrameTicks->{Range[0,1,0.2],ComplexTicks[Range[-0.2,0.2,0.2]],None,None},AspectRatio->Automatic]

The Mandelbrot set

MandelbrotIterations[c_,bail_]:=JuliaIterationsCompiled[0,c,bail]

DensityPlot[MandelbrotIterations[x+I y,100],{x,-2,0.6},{y,-1.3,1.3}, FrameTicks->{Automatic, ComplexTicks[ Range[-1,1,0.5]], None, None},PlotPoints->100, ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&)]

DensityPlot[MandelbrotIterations[x+I y,100],{x,-0.65,-0.45},{y,0.5,0.71}, FrameTicks->{Automatic,ComplexTicks[ Range[0.5,0.7,0.1]],None,None},PlotPoints->100, ColorFunction->(If[#==1,Black,Hue[0.9 ((#)^(3/4))]]&)]

Boundary scanning

ListConvolve[(({{(Subscript[k, 11]), (Subscript[k, 12])}, {(Subscript[k, 21]), (Subscript[k, 22])}})),(({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[a, 23])}, {(Subscript[a, 31]), (Subscript[a, 32]), (Subscript[a, 33])}}))]//MatrixForm

(({{1, 1, 1}, {1, -8, 1}, {1, 1, 1}})).

kernel={{1,1,1},{1,-8,1},{1,1,1}};
data=Table[If[Abs[Norm[{x,y}]-Floor[Norm[{x,y}]]]<1/3,0,1],{x,-5,5,10/200},{y,-5,5,10/200}];
convolvedData=ListConvolve[kernel,data];
ArrayPlot[convolvedData,ColorFunction->(Hue[(((1-#))^(1/2))]&)]

data=Reverse[Table[MandelbrotIterations[a+b I,100],{b,-1.3,1.3,(2.6/501)},{a,-2,0.6,(2.6/501)}]];
convolvedData=ListConvolve[kernel,data];
MandelbrotBW=ArrayPlot[Abs[convolvedData],Frame->False,ImageSize->250,DataRange->{{-2,0.6},{-1.3,1.3}},ColorFunction->(GrayLevel[(((1-#))^5)]&)]

A dynamic look at Julia sets

Manipulate[Grid[{{MandelbrotBW,
Dynamic[JuliaModified[pt[[1]]+pt[[2]] I,Resolution -> 100,
ImageSize -> 250,PlotRange -> {{-1.8,1.8},{-1.8,1.8}}]]},
{Text[Dynamic[If[pt[[2]]≥0,
StringForm["c = `` + `` I",PaddedForm[pt[[1]],{5,4}],PaddedForm[pt[[2]],{5,4}]],
StringForm["c = `` - `` I",PaddedForm[pt[[1]],{5,4}],PaddedForm[-pt[[2]],{5,4}]]]]], SpanFromLeft}
}],{{pt,{-0.75,0}},{-2,-1.3},{0.6,1.3},Locator}]