(*********************************************************************** Mathematica-Compatible Notebook This notebook can be used on any computer system with Mathematica 4.0, MathReader 4.0, or any compatible application. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. ***********************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 34927, 1127]*) (*NotebookOutlinePosition[ 35966, 1160]*) (* CellTagsIndexPosition[ 35922, 1156]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Inverse Iteration Algorithms for Julia Sets", "Title", PageWidth->PaperWidth], Cell["by Mark McClure ", "Author", PageWidth->PaperWidth], Cell["\<\ Inverse iteration algorithms are extremely fast methods to generate \ images of Julia sets. While they are fairly simple to understand and \ implement, a little thought can improve their performance dramatically. In \ this article, we discuss several variations of the basic inverse iteration \ scheme.\ \>", "Abstract", PageWidth->PaperWidth], Cell["\<\ \ \>", "Text"], Cell[CellGroupData[{ Cell["Introduction", "Section", PageWidth->PaperWidth], Cell[TextData[{ "Julia sets form a beautiful and varied class of fractals. Spectacular \ color images of these sets may be generated using the well known escape time \ algorithm (see [Peitgen and Richter 1989]). In this article, however, we \ discuss several variations of the inverse iteration algorithm. Inverse \ iteration algorithms are extremely fast, broadly applicable, and easily \ understood. Furthermore, an understanding of inverse iteration illuminates \ the connection between Julia sets of quadratic functions, more general \ rational functions and, even, self-similar sets.\n We begin with a brief \ look at the theory of Julia sets. An intuitive understanding of the theory \ will explain the idea behind inverse iteration and why it works. We then \ look at the implementation of inverse iteration for certain quadratic \ functions, carefully refine the technique, and generalize it to more \ arbitrary functions. Finally, we discuss a stochastic algorithm, which is \ reminiscent of the chaos game.\n All the code in this article is \ encapsulated in the package ", StyleBox["JuliaSet", "Input"], ", whose use will be illustrated throughout." }], "Text", PageWidth->PaperWidth] }, Closed]], Cell[CellGroupData[{ Cell["The Theory of Fatou and Julia", "Section", PageWidth->PaperWidth], Cell[TextData[{ "Julia Sets arise in the study of complex dynamics. This study was \ initiated in the late teens through the seminal works of Fatou and Julia \ [Fatou 1919; Julia 1918]. Many expositions at various levels have appeared \ in recent years. [Beardon1991] and [Carleson and Gamelin 1993] are both \ outstanding texts, the first being more advanced. [Devaney 1992] is an \ excellent choice for those with a calculus background. It, also, has \ information on algorithms, including a version of the inverse iteration \ algorithm (sec. 16.6).\n In complex dynamics, we study the iteration of a \ function mapping the complex plane into itself, ", Cell[BoxData[ \(TraditionalForm\`\(\ f : \(\[DoubleStruckCapitalC]\ \[LongRightArrow]\ \[DoubleStruckCapitalC] . \)\)\)]], " The function, ", Cell[BoxData[ \(TraditionalForm\`\(f, \)\)]], " is usually assumed to be rational. The complex plane, ", Cell[BoxData[ \(TraditionalForm\`\(\[DoubleStruckCapitalC], \)\)]], " may be thought of as the state space and the application of ", Cell[BoxData[ \(TraditionalForm\`f\)]], " represents the passage of one unit of time. If ", Cell[BoxData[ \(TraditionalForm\`z\_0\)]], " is an initial point, then the position of ", Cell[BoxData[ \(TraditionalForm\`z\_0\)]], " after the passage of ", Cell[BoxData[ \(TraditionalForm\`n\)]], " units of time is ", Cell[BoxData[ \(TraditionalForm\`\(f\^n\)(z\_0)\)]], " and the orbit of ", Cell[BoxData[ \(TraditionalForm\`z\_0\)]], " is defined to be ", Cell[BoxData[ \(TraditionalForm\`\({\(f\^n\)(z\_0)}\_n . \)\)]], " \n In this study, two subsets of ", Cell[BoxData[ \(TraditionalForm\`\[DoubleStruckCapitalC]\)]], " naturally arise - the Fatou set and the Julia set. The ", StyleBox["Fatou set", FontSlant->"Italic"], ", ", Cell[BoxData[ \(TraditionalForm\`F\)]], ", may be defined to be the largest open set on which the set of iterates \ of ", Cell[BoxData[ \(TraditionalForm\`f\)]], ", ", Cell[BoxData[ \(TraditionalForm\`\({f\^n}\_n, \)\)]], " form a normal family. Intuitively, this may be thought of as the largest \ open set on which the dynamics of ", Cell[BoxData[ \(TraditionalForm\`f\)]], " are relatively tame in the sense that points close to one another have \ similar long term behavior. The ", StyleBox["Julia set", FontSlant->"Italic"], ", ", Cell[BoxData[ \(TraditionalForm\`\(J, \)\)]], " is defined to be the complement of the Fatou set and the dynamics of ", Cell[BoxData[ \(TraditionalForm\`f\)]], " are quite chaotic on ", Cell[BoxData[ \(TraditionalForm\`J\)]], ". \n For example, suppose ", Cell[BoxData[ \(TraditionalForm\`f(z)\ = \ z\^2. \)]], " Then, ", Cell[BoxData[ \(TraditionalForm\`J\ = \ {z : \(| z | \) = 1}\)]], ", the unit circle. To see this, note that if ", Cell[BoxData[ \(TraditionalForm\`z\_0\)]], " is in the interior of the unit circle, then ", Cell[BoxData[ \(TraditionalForm\`\(\(f\^n\)(z\_0)\)\ \[LongRightArrow]\ 0\)]], " as ", Cell[BoxData[ \(TraditionalForm\`\(n\ \[LongRightArrow]\ \[Infinity] . \)\)]], " Thus, the entire interior of the unit circle is attracted to the fixed \ point ", Cell[BoxData[ \(TraditionalForm\`0\)]], ". Similarly, if ", Cell[BoxData[ \(TraditionalForm\`z\_0\)]], " is in the exterior of the unit circle, then ", Cell[BoxData[ \(TraditionalForm \`\(\(f\^n\)(z\_0)\)\ \[LongRightArrow]\ \[Infinity]\)]], " as ", Cell[BoxData[ \(TraditionalForm\`n\ \[LongRightArrow]\ \[Infinity]\)]], ". Again, the dynamics of any point outside the unit circle is similar to \ that of any other point, outside the unit circle. On the other hand, the \ dynamics on the unit circle itself are quite complex. Note that any point ", Cell[BoxData[ \(TraditionalForm\`z\)]], " on the unit circle may be written ", Cell[BoxData[ \(TraditionalForm\`z\ = \ \(e\^\(\[Alpha]\ \[Pi]\ i\) . \)\)]], " It may be shown that if ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\)]], " is rational, then the orbit of ", Cell[BoxData[ \(TraditionalForm\`z\)]], " consists of finitely many points, while the orbit of ", Cell[BoxData[ \(TraditionalForm\`z\)]], " is dense in the unit circle for irrational ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\)]], " [Beardon 1991, sec. 1.3]. Thus we may find distinct points in the unit \ circle, arbitrarily close to one another, with distinctly different long term \ behavior.\n We will begin by developing algorithms to generate Julia sets \ for functions of the form ", Cell[BoxData[ \(TraditionalForm\`f\_c = z\^2 + \(c . \)\)]], " We denote the Julia Set of ", Cell[BoxData[ \(TraditionalForm\`f\_c\)]], " by ", Cell[BoxData[ \(TraditionalForm\`J\_\(c . \)\)]], " This seemingly small subset of the rational functions is much less \ restrictive than it may first appear. It may be shown that the dynamical \ behavior of any quadratic function is exhibited by ", Cell[BoxData[ \(TraditionalForm\`f\_c\)]], " for exactly one ", Cell[BoxData[ \(TraditionalForm\`c\)]], " [Carleson and Gamelin 1993, p. 123]. Thus, the Julia set of any \ quadratic will be closely related (homeomorphic) to ", Cell[BoxData[ \(TraditionalForm\`J\_c\)]], " for some ", Cell[BoxData[ \(TraditionalForm\`\(c . \)\)]], " Furthermore, our techniques will easily generalize to more arbitrary \ functions." }], "Text", PageWidth->PaperWidth] }, Closed]], Cell[CellGroupData[{ Cell[" A Simple Deterministic Inverse Iteration Algorithm", "Section", PageWidth->PaperWidth], Cell[TextData[{ "Note that in the preceeding example, points in ", Cell[BoxData[ \(TraditionalForm\`F\)]], " were repelled from ", Cell[BoxData[ \(TraditionalForm\`J\)]], " under the dynamics of ", Cell[BoxData[ \(TraditionalForm\`\(f.\)\)]], " This is true in general. In fact, ", Cell[BoxData[ \(TraditionalForm\`J\)]], " may be characterized as the closure of the set of repelling periodic \ points of ", Cell[BoxData[ \(TraditionalForm\`f\)]], " [Beardon 1991, Thm. 6.9.2] Thus, the Julia set of ", Cell[BoxData[ \(TraditionalForm\`f\)]], " may be thought of as a dynamical repeller for ", Cell[BoxData[ \(TraditionalForm\`\(f.\)\)]], " This observation lies at the heart of all inverse iteration algorithms. \ Let's use it to construct a first algorithm to generate the Julia set, ", Cell[BoxData[ \(TraditionalForm\`\(J\_c, \)\)]], " for the function ", Cell[BoxData[ \(TraditionalForm\`\(f\_c\)(z)\ = \ z\^\(\ 2\)\ + \ \(c.\)\)]], " Since ", Cell[BoxData[ \(TraditionalForm\`J\_c\)]], " is a repeller for ", Cell[BoxData[ \(TraditionalForm\`\(f\_c, \)\)]], " it should, also, be an attractor for an inverse of ", Cell[BoxData[ \(TraditionalForm\`\(f\_c.\)\)]], " Of course, ", Cell[BoxData[ \(TraditionalForm\`f\_c\)]], " has two inverses ", Cell[BoxData[ \(TraditionalForm\`f\_1\%\(-1\)\ = \ \@\(z - c\)\)]], " and ", Cell[BoxData[ \(TraditionalForm\`f\_2\%\(-1\)\ = \ \(-\(\@\(z - c\).\)\)\)]], " Thus, if ", Cell[BoxData[ \(TraditionalForm\`{z\_0}\)]], " is an initial point, then ", Cell[BoxData[ \(TraditionalForm \`{\(f\_1\%\(-1\)\)(z\_0), \ \(f\_2\%\(-1\)\)(z\_0)}\)]], " are two points that will be closer to the Julia set ", Cell[BoxData[ \(TraditionalForm\`\(J\_c, \)\)]], " ", Cell[BoxData[ \(TraditionalForm \`{\(f\_1\%\(-1\)\)(\(f\_1\%\(-1\)\)(z\_0)), \ \(f\_1\%\(-1\)\)(\(f\_2\%\(-1\)\)(z\_0)), \ \(f\_2\%\(-1\)\)(\(f\_1\%\(-1\)\)(z\_0)), \ \(f\_2\%\(-1\)\)(\(f\_2\%\(-1\)\)(z\_0))}\)]], " are four points that will be even closer to ", Cell[BoxData[ \(TraditionalForm\`\(J\_c, \)\)]], " and so on.\n To apply the above idea, we need a ", StyleBox["Mathematica", FontSlant->"Italic"], " function, ", StyleBox["invImage", "Input"], ", that accepts a list of complex numbers and returns the inverse image. \ Our first attempt will generate ", Cell[BoxData[ \(TraditionalForm\`J\_c\)]], " for ", Cell[BoxData[ \(TraditionalForm\`\(c = 0, \)\)]], " which we know to be the unit circle." }], "Text", PageWidth->PaperWidth], Cell["\<\ c = 0.; invImage[complexPoints_] := \tFlatten[({1, -1} Sqrt[#1 - c] & ) /@ \tcomplexPoints]; \ \>", "Input", CellLabel->"In[1]:=", PageWidth->PaperWidth], Cell["invImage[{1}]", "Input", CellLabel->"In[2]:=", PageWidth->PaperWidth], Cell["invImage[%]", "Input", CellLabel->"In[3]:=", PageWidth->PaperWidth], Cell[TextData[{ "Now, we'll nest ", StyleBox["invImage", "Input"], " several times and plot the points after converting them to ordered \ pairs." }], "Text", PageWidth->PaperWidth], Cell["\<\ depth = 10; points = ({Re[#1], Im[#1]} & ) /@ Nest[invImage, \t{1},depth]; ListPlot[points, AspectRatio -> Automatic, \tAxes -> False, \tPlotStyle -> {AbsolutePointSize[0.4]}]\ \>", "Input", CellLabel->"In[4]:=", PageWidth->PaperWidth], Cell[TextData[{ " The package function which encapsulates these commands is ", StyleBox["JuliaSimple", "Input"], ". Note that the second argument to ", StyleBox["JuliaSimple", "Input"], " represents the depth rather than the total number of points. Thus, the \ following command generates ", Cell[BoxData[ \(TraditionalForm\`2\^12\)]], " points of an interesting Julia set." }], "Text", PageWidth->PaperWidth], Cell["Needs[\"JuliaSet`\"]", "Input", CellLabel->"In[5]:=", PageWidth->PaperWidth], Cell["\<\ JuliaSimple[-0.123 + 0.745 I, 12, \tPlotLabel -> \"c = -0.123 + 0.745 I\"]\ \>", "Input", CellLabel->"In[6]:=", PageWidth->PaperWidth] }, Closed]], Cell[CellGroupData[{ Cell["Improving The Simple Deterministic Algorithm", "Section", PageWidth->PaperWidth], Cell[TextData[{ StyleBox["JuliaSimple", "Input"], " is a short, fast, and easily understood algorithm. However, our first \ interesting image above already indicates a weakness in the algorithm. Some \ parts of the image seem more detailed than others. This is because some \ parts of the Julia set are more attractive than others under the action of ", StyleBox["invImage", "Input"], ". In this section, we'll describe the function ", StyleBox["JuliaModified", "Input"], " which fixes the problem. Here is an example which compares the \ algorithms." }], "Text", PageWidth->PaperWidth], Cell["\<\ Show[GraphicsArray[{ \t{JuliaSimple[0.68 I, 12, \t\tDisplayFunction -> Identity, \t\tPlotLabel -> \"Simple\"]}, \t{JuliaModified[0.68 I, \t\tDisplayFunction -> Identity, \t\tPlotLabel -> \"Modified\"]}}], \tPlotLabel -> \"c = 0.68 I\"]\ \>", "Input", CellLabel->"In[7]:=", PageWidth->PaperWidth], Cell[TextData[{ " The simple algorithm clearly misses a great amount of detail. A \ naive solution would be to increase the depth. However, the required depth \ is frequently greater than ", Cell[BoxData[ \(TraditionalForm\`\(50, \)\)]], " resulting in more than ", Cell[BoxData[ \(TraditionalForm\`2\^50\)]], " points in the image.\n Here is a solution. After a certain depth, \ we'll keep track of all of the points that have been plotted. Points close \ together, as measured by the variable ", StyleBox["res", "Input"], ", will be treated the same. With each iteration, we'll apply ", StyleBox["invImage", "Input"], " and discard points that have already been plotted. As a result, the \ program will be able to run to a much larger depth since the length of the \ list of points no longer grows as ", Cell[BoxData[ \(TraditionalForm\`\(2\^depth.\)\)]], " \n First, we need to rewrite ", StyleBox["invImage", "Input"], " to treat points close together similarly. We'll regenerate ", Cell[BoxData[ \(TraditionalForm\`J\_c\)]], " for ", Cell[BoxData[ \(TraditionalForm\`c = \(-0.123\)\ + \ 0.745 \( I.\)\)]] }], "Text", PageWidth->PaperWidth], Cell["\<\ c = -0.123 + 0.745 I; res = 200; invImage[points_] := Flatten[(Floor[{1, -1} \tres Sqrt[#1 - c]]/res & ) /@ \tpoints]; \ \>", "Input", CellLabel->"In[8]:=", PageWidth->PaperWidth], Cell[TextData[{ "The function ", StyleBox["reducedInvImage", "Input"], " will accept a list of points, apply ", StyleBox["invImage", "Input"], ", and return only those points that don't appear in the auxiliary \ variable, ", StyleBox["pointsSoFar", "Input"], ". As a side effect, it will update ", StyleBox["pointsSoFar", "Input"], "." }], "Text", PageWidth->PaperWidth], Cell["\<\ reducedInvImage[points_] := Module[{newPoints}, \tnewPoints = Complement[invImage[points], \t\tpointsSoFar]; \tpointsSoFar = Union[newPoints, pointsSoFar]; \tnewPoints]; \ \>", "Input", CellLabel->"In[9]:=", PageWidth->PaperWidth], Cell[TextData[{ "Next, we'll iterate ", StyleBox["invImage", "Input"], " several times from an arbitrary starting value to obtain some points \ close to ", Cell[BoxData[ \(TraditionalForm\`\(J\_c.\)\)]], " We'll store the result in ", StyleBox["pointsSoFar", "Input"], "." }], "Text", PageWidth->PaperWidth], Cell["pointsSoFar = Nest[invImage, {1}, 5]; ", "Input", CellLabel->"In[10]:=", PageWidth->PaperWidth], Cell[TextData[{ "Now, we'll iterate ", StyleBox["reducedInvImage", "Input"], ", starting with ", StyleBox["pointsSoFar", "Input"], ", until it returns the empty list." }], "Text", PageWidth->PaperWidth], Cell["FixedPoint[reducedInvImage, pointsSoFar]; ", "Input", CellLabel->"In[11]:=", PageWidth->PaperWidth], Cell[TextData[{ "Finally, we ", StyleBox["ListPlot", "Input"], " the points in ", StyleBox["pointsSoFar", "Input"], ", after converting the complex numbers to ordered pairs." }], "Text", PageWidth->PaperWidth], Cell["\<\ ListPlot[({Re[#1], Im[#1]} & ) /@ pointsSoFar, \tAspectRatio -> Automatic, \tAxes -> False, \tPlotStyle -> {AbsolutePointSize[0.4]}, \tPlotLabel -> \"c = -0.123 + 0.745 I, Modified\"]\ \>", "Input", CellLabel->"In[12]:=", PageWidth->PaperWidth], Cell[TextData[{ " Note that the level of detail is controlled by the option ", StyleBox["Resolution", "Input"], ". The default is ", StyleBox["Resolution\[Rule]200", "Input"], ". We may get a quick, less detailed image by using a smaller value for ", StyleBox["Resolution", "Input"], ". Here is an illustration of the effect of ", StyleBox["Resolution", "Input"], "." }], "Text", PageWidth->PaperWidth], Cell["\<\ Show[GraphicsArray[{ \t{JuliaModified[-0.77 + 0.22 I, \t\tResolution -> 40, \t\tPlotLabel -> \"Resolution \[Rule] 40\", \t\tDisplayFunction -> Identity]}, {JuliaModified[-0.77 + 0.22*I, \t\tPlotLabel -> \"Resolution \[Rule] 200\", \t\tDisplayFunction -> Identity]}}], \tPlotLabel -> \"c = -0.77 + 0.22 I\"]\ \>", "Input", CellLabel->"In[13]:=", PageWidth->PaperWidth], Cell[CellGroupData[{ Cell["Comparing The Algorithms", "Subsection", PageWidth->PaperWidth], Cell[TextData[{ "The above images show much greater detail than was possible with ", StyleBox["JuliaSimple", "Input"], ". We can measure how large ", StyleBox["depth", "Input", FontWeight->"Bold"], " would have to be in ", StyleBox["JuliaSimple", "Input"], " to achieve the same level of detail by replacing ", StyleBox["FixedPoint", "Input"], " with ", StyleBox["FixedPointList", "Input"], " and measuring its length. Note that this computation will depend on the \ complexity of the Julia set. We'll measure the depth for a fairly \ complicated Julia set." }], "Text", PageWidth->PaperWidth], Cell["\<\ c = 0.68 I; res = 200; invImage[points_] := Flatten[(Floor[{1, -1} \tres Sqrt[#1 - c]]/res & ) /@ \t\tpoints]; reducedInvImage[points_] := Module[{newPoints}, \tnewPoints = Complement[invImage[points], \t\tpointsSoFar]; \tpointsSoFar = Union[newPoints, pointsSoFar]; \tnewPoints]; pointsSoFar = Nest[invImage, {1}, 5]; Length[FixedPointList[reducedInvImage, pointsSoFar]]\ \>", "Input", CellLabel->"In[14]:=", PageWidth->PaperWidth], Cell[TextData[{ "Thus, ", StyleBox["JuliaSimple", "Input"], " would need a depth of 77+5 to plot this particular Julia set resulting in \ a list of ", Cell[BoxData[ \(TraditionalForm\`2\^82\)]], " points. ", StyleBox["JuliaModified", "Input"], " used barely 20,000 points to plot much more detail, as the following \ computation shows." }], "Text", PageWidth->PaperWidth], Cell["Length[pointsSoFar]", "Input", CellLabel->"In[15]:=", PageWidth->PaperWidth] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["The Mandelbrot Set", "Section", PageWidth->PaperWidth], Cell[TextData[{ "Julia sets for quadratic functions of the form ", Cell[BoxData[ \(TraditionalForm\`z\^2 + c\)]], " are among the most important, because of their relationship with the \ Mandelbrot set. The Mandelbrot set is defined to be the set of all ", Cell[BoxData[ \(TraditionalForm\`c\)]], " values that lead to connected Julia sets, ", Cell[BoxData[ \(TraditionalForm\`J\_c\)]], ". While it is not the purpose of this article to discuss the Mandelbrot \ set in detail, it is nice to have an image of it to assist in choosing \ interesting values of ", Cell[BoxData[ \(TraditionalForm\`\(c.\)\)]], " The following code generates an image of the Mandelbrot set. It is a \ minor modification of the code found in [Dickau 1997]. Points near the \ boundary frequently lead to interesting Julia sets." }], "Text", PageWidth->PaperWidth], Cell["\<\ MandelbrotFunction = Compile[{{c, _Complex}}, \tLength[FixedPointList[#1^2 + c & , c, 100, \t\tSameTest -> (Abs[#2] > 2. & )]]]; DensityPlot[MandelbrotFunction[x + y*I], \t{x, -2, 0.6}, {y, -1.3, 1.3}, \tMesh -> False, AspectRatio -> Automatic, \tPlotPoints -> 300, \tColorFunction -> (If[#1 == 1, \t\tRGBColor[0, 0, 0], Hue[0.9*#1]] & )]; \ \>", "Input", CellLabel->"In[16]:=", PageWidth->PaperWidth] }, Closed]], Cell[CellGroupData[{ Cell["Generalizing the Deterministic Algorithm", "Section", PageWidth->PaperWidth], Cell[TextData[{ "The algorithms described above apply only to functions of the specific \ form ", Cell[BoxData[ \(TraditionalForm\`\(\(f\_c\)(z) = z\^2 + c, \)\)]], " but the theory is much more broadly applicable. We should be able to \ apply the same ideas to any rational function whose inverses are known. \ Let's apply these ideas to plot the Julia set of ", Cell[BoxData[ \(TraditionalForm\`f(z)\ = z\^3 + z + .6 \( I.\)\)]] }], "Text", PageWidth->PaperWidth], Cell["\<\ f[z_] = z^3 + z + 0.6*I; res = 200; inverses = z /. NSolve[f[z] == #1, z]; funcs = (Function[anInverse, N[Floor[anInverse*res]/res] & ]) /@ inverses; \ \ \>", "Input", CellLabel->"In[17]:=", PageWidth->PaperWidth], Cell[TextData[{ "A list of the inverses of ", Cell[BoxData[ \(TraditionalForm\`f(z)\)]], " is now held in ", StyleBox["funcs", "Input"], ". We need to turn this into ", StyleBox["invImage", "Input"], "." }], "Text", PageWidth->PaperWidth], Cell["\<\ invImage[points_] := \tFlatten[(Through[funcs[#1]] & ) /@ points, 1]; \tinvImage[{1}]\ \>", "Input", CellLabel->"In[18]:=", PageWidth->PaperWidth], Cell["invImage[%]", "Input", CellLabel->"In[19]:=", PageWidth->PaperWidth], Cell["Now the algorithm proceeds as before.", "Text", PageWidth->PaperWidth], Cell["\<\ reducedImage[points_] := Module[{newPoints}, \tnewPoints = Complement[invImage[points], \t\tpointsSoFar]; \tpointsSoFar = \t\tUnion[newPoints, pointsSoFar]; newPoints]; pointsSoFar = Nest[invImage, {1.}, 5]; FixedPoint[reducedImage, pointsSoFar]; ListPlot[({Re[#1], Im[#1]} & ) /@ pointsSoFar, \tAspectRatio -> Automatic, \tAxes -> False, \tPlotStyle -> {AbsolutePointSize[0.4]}]\ \>", "Input", CellLabel->"In[20]:=", PageWidth->PaperWidth], Cell[TextData[{ " These commands are encapsulated in the package function ", StyleBox["Julia", "Input"], ". Here are two more examples." }], "Text", PageWidth->PaperWidth], Cell[BoxData[ \(Show[ GraphicsArray[\n\t{{Julia[z\^3 - I, z, \n\t\tPlotLabel\ \[Rule] \ \*"\"\\"", \n\t\tDisplayFunction\ \[Rule] \ Identity]}, \n\t{Julia[ 1\/\(z\^2 - 1\), z, \n\t\tPlotLabel\ \[Rule] \ \*"\"\\"", \n\t\tDisplayFunction\ \[Rule] \ Identity]}}]]\)], "Input", CellLabel->"In[21]:=", PageWidth->PaperWidth], Cell[TextData[{ " Note that the second image in the preceeding example is the Julia set \ of a rational function. We, actually, need to make one final adjustment to \ generate Julia sets of rational functions. Our method for pruning points \ depends upon the fact that the Julia set is bounded. While this assumption \ is valid for polynomials, the Julia set of a rational function need not be \ bounded. In this case, the function ", StyleBox["reducedImage", "Input"], " should throw out points that are larger than some specified bound, in \ addition to points that have already been plotted. Here is the modified \ version of ", StyleBox["reducedImage", "Input"], "." }], "Text", PageWidth->PaperWidth], Cell["\<\ bound = 4; reducedImage[points_] := Module[{newPoints}, \tnewPoints = \t\tComplement[image[points], pointsSoFar]; \tnewPoints = \t\tSelect[newPoints, N[Abs[#1]] <= bound & ]; \tpointsSoFar = \t\tUnion[newPoints, pointsSoFar]; \tnewPoints]\ \>", "Input", CellLabel->"In[22]:=", PageWidth->PaperWidth], Cell[TextData[{ "The package function ", StyleBox["Julia", "Input"], " uses this version of ", StyleBox["reducedImage", "Input"], " when called with a rational function. The value of ", StyleBox["bound", "Input"], " is controlled by the option ", StyleBox["Bound", "Input"], ". The default is ", StyleBox["Bound\[Rule]4.", "Input"], " We illustrate the use of ", StyleBox["Bound", "Input"], " by plotting the Julia set of the rational function used to find the roots \ of ", Cell[BoxData[ \(TraditionalForm\`z\^3 = 1\)]], " with Newton's method." }], "Text", PageWidth->PaperWidth], Cell[BoxData[ \(Julia[\(2 z\)\/3 + 1\/\(3 z\^2\), z, \n\tBound\ \[Rule] \ 12, \n\t PlotLabel\ \[Rule] \ \*"\"\\""]\)], "Input", CellLabel->"In[23]:="] }, Closed]], Cell[CellGroupData[{ Cell["A Stochastic Algorithm", "Section", PageWidth->PaperWidth], Cell[TextData[{ "This article would not be complete without a discussion of the random \ inverse iteration algorithm. This is probably the simplest, and fastest, way \ to generate the Julia set of a quadratic function. The random inverse \ iteration algorithm, also, highlights the similarities between Julia sets and \ self-similar sets, because of its resemblance to the chaos game [Devaney \ 1992, sec. 14.1]. \n\nSuppose we would like to generate 1000 points of the \ Julia set for ", Cell[BoxData[ \(TraditionalForm\`c = 0.\)]] }], "Text", PageWidth->PaperWidth], Cell["c = 0.; numPoints = 1000; ", "Input", CellLabel->"In[24]:=", PageWidth->PaperWidth], Cell[TextData[{ "As we've seen, ", Cell[BoxData[ \(TraditionalForm\`J\_c\)]], " is attractive under the action of both inverses of ", Cell[BoxData[ \(TraditionalForm\`\(f\_c.\)\)]] }], "Text", PageWidth->PaperWidth], Cell["inv1 = Sqrt[#1 - c] & ; inv2 = -Sqrt[#1 - c] & ; ", "Input", CellLabel->"In[25]:=", PageWidth->PaperWidth], Cell[TextData[{ "We choose the number 1 as an arbitrary starting point. We then choose one \ of the inverses ", StyleBox["inv1", "Input"], " or ", StyleBox["inv2", "Input"], " randomly and apply it to the starting point 1. Continuing, we iterate \ this procedure, choosing a random inverse and applying it to the previous \ point. This generates a list of points converging to ", Cell[BoxData[ \(TraditionalForm\`\(J\_c.\)\)]], " We drop the first few points, since they might not be close to ", Cell[BoxData[ \(TraditionalForm\`\(J\_c, \)\)]], " and plot the rest. We implement this idea in ", StyleBox["Mathematica", FontSlant->"Italic"], " by using ", StyleBox["Table", "Input"], " to generate a list of length ", StyleBox["numPoints,", "Input"], " where each entry is either ", StyleBox["inv1", "Input"], " or ", StyleBox["inv2", "Input"], " with equal probability. We then use ", StyleBox["ComposeList", "Input"], " to build the list of points." }], "Text", PageWidth->PaperWidth], Cell["\<\ points = Drop[ComposeList[Table[chooser = Random[]; \tIf[chooser < 0.5, inv1, \t\tinv2], {numPoints}], 1], 10]; \ \>", "Input", CellLabel->"In[26]:=", PageWidth->PaperWidth], Cell["\<\ Finally, we plot these points, after converting them to ordered \ pairs.\ \>", "Text", PageWidth->PaperWidth], Cell["\<\ ListPlot[({Re[#1], Im[#1]} & ) /@ points, \tAspectRatio -> Automatic, \tAxes -> False, \tPlotStyle -> {AbsolutePointSize[0.4]}]\ \>", "Input", CellLabel->"In[27]:=", PageWidth->PaperWidth], Cell[TextData[{ " These commands are encapsulated in the package function ", StyleBox["JuliaStochastic", "Input"], ". Unfortunately, images generated in this manner frequently display the \ same nonuniform distribution we saw with ", StyleBox["JuliaSimple", "Input"], ". The random inverse iteration algorithm may be improved, for some Julia \ sets, by skewing the probability of choosing one inverse over the other. \ This accomplished in our function ", StyleBox["JuliaStochastic", "Input"], " with the option ", StyleBox["OneBias", "Input"], ", which is a real number, strictly between 0 and 1, indicating the \ probability of choosing ", StyleBox["inv1", "Input"], ". The probability of choosing ", StyleBox["inv2", "Input"], " is then, necessarily, ", StyleBox["1-OneBias", "Input"], ". The following example illustrates the effect of ", StyleBox["OneBias", "Input"], "." }], "Text", PageWidth->PaperWidth], Cell["\<\ Show[GraphicsArray[{ \t{JuliaStochastic[-1, 8000, \t\tPlotLabel -> \"OneBias \[Rule] 0.5\", \t\tDisplayFunction -> Identity]}, {JuliaStochastic[-1, 8000, \t\tOneBias -> 0.28, \t\tPlotLabel -> \"OneBias \[Rule] 0.28\", \t\tDisplayFunction -> Identity]}}], \tPlotLabel -> \"c = -1\"]\ \>", "Input", CellLabel->"In[28]:=", PageWidth->PaperWidth], Cell[TextData[{ " In general, the inverse which is more contractive about its fixed \ point should be assigned a lower probability. The contractivity of a \ function about its fixed point is measured by the value of its derivative. \ Thus for ", Cell[BoxData[ \(TraditionalForm\`\(c = \(-1\), \)\)]], " as above we have the following." }], "Text", PageWidth->PaperWidth], Cell["\<\ c = -1.; inv1 = Sqrt[#1 - c] & ; inv2 = -Sqrt[#1 - c] & ; {inv1FixedPoint, inv2FixedPoint} = Chop[{FixedPoint[inv1, 1, 1000], \tFixedPoint[inv2, 1, 1000]}]\ \>", "Input", CellLabel->"In[29]:=", PageWidth->PaperWidth], Cell["With the fixed points, we may now measure the contractivity.", "Text", PageWidth->PaperWidth], Cell["\<\ {Abs[D[inv1[x], x]] /. \tx -> inv1FixedPoint, Abs[D[inv2[x], x]] /. \tx -> inv2FixedPoint}\ \>", "Input", CellLabel->"In[30]:=", PageWidth->PaperWidth], Cell[TextData[{ "This suggests we should choose ", StyleBox["OneBias", "Input"], " to be less than .5, since ", StyleBox["inv1", "Input"], " is more contractive than ", StyleBox["inv2", "Input"], ". Unfortunately, many Julia sets have much more detail than the random \ inverse algorithm will generate for any choice of ", StyleBox["OneBias", "Input"], ". For example, suppose ", Cell[BoxData[ \(TraditionalForm\`c = 0.68\ \(I.\)\)]], " " }], "Text", PageWidth->PaperWidth], Cell["JuliaStochastic[0.68 I, 40000, OneBias -> 0.23]", "Input", CellLabel->"In[31]:=", PageWidth->PaperWidth], Cell[TextData[{ " Experimentation suggests that a value of 0.23 is close to optimal for \ ", StyleBox["OneBias", "Input"], ". However, the modified inverse iterated algorithm shows much more detail \ using far fewer points. Therefore, we will not pursue the random algorithm \ further." }], "Text", PageWidth->PaperWidth] }, Closed]], Cell[CellGroupData[{ Cell["References", "Section", PageWidth->PaperWidth], Cell[TextData[{ "Beardon, A.F. 1991. ", StyleBox["Iteration of Rational Functions. ", FontSlant->"Italic"], "Springer-Verlag New York Graduates Texts in Mathematics 132." }], "Reference", PageWidth->PaperWidth], Cell[TextData[{ "Carleson, L. and Gamelin, T.W. 1993. ", StyleBox["Complex Dynamics", FontSlant->"Italic"], ". Springer-Verlag New York." }], "Reference", PageWidth->PaperWidth], Cell[TextData[{ "Devaney, R.L. 1992. ", StyleBox["A First Course in Chaotic Dynamical Systems. ", FontSlant->"Italic"], "Addison-Wesley Reading, Mass." }], "Reference", PageWidth->PaperWidth], Cell[TextData[{ "Dickau, R.M. 1997. Compilation of iterative and list operations in \ \"Tricks of the Trade.\" ", StyleBox["The Mathematica Journal.", FontSlant->"Italic"], " 7(1): 14-15." }], "Reference", PageWidth->PaperWidth], Cell[TextData[{ "Fatou, M.P. 1919. Sur les equations fonctionelles. ", StyleBox[" Bull. Soc. Math. France.", FontSlant->"Italic"], " 47:161-271." }], "Reference", PageWidth->PaperWidth], Cell[TextData[{ "Julia, G. 1918. Memoire sur l'iteration des fonctions rationelles. ", StyleBox["Math. Pures Appl.,", FontSlant->"Italic"], " 8:47-245." }], "Reference", PageWidth->PaperWidth], Cell[TextData[{ "Peitgen, H.O. and Richter, P. 1989. ", StyleBox["The Beauty of Fractals.", FontSlant->"Italic"], " Springer-Verlag Heidelberg." }], "Reference", PageWidth->PaperWidth] }, Closed]], Cell[CellGroupData[{ Cell["About the Author", "Section"], Cell[TextData[{ "Mark McClure is an assistant professor of mathematics at the University of \ North Carolina at Asheville. He received his Ph.D. in mathematics from the \ Ohio State University under the direction of Gerald Edgar. His primary \ research interest is the study of fractal dimensions of various sets arising \ in real analysis. He was introduced to ", StyleBox["Mathematica", FontSlant->"Italic"], " through the Calculus&", StyleBox["Mathematica", FontSlant->"Italic"], " program at OSU." }], "Text", PageWidth->PaperWidth] }, Closed]], Cell["\<\ Mark McClure Department of Mathematics University of North Carolina at Asheville Asheville, North Carolina 28804 mcmcclur@bulldog.unca.edu http://www.unca.edu/~mcmcclur/\ \>", "Address", PageWidth->PaperWidth] }, Open ]] }, FrontEndVersion->"4.0 for Macintosh", ScreenRectangle->{{0, 832}, {0, 604}}, AutoGeneratedPackage->None, ScreenStyleEnvironment->"Working", PrintingStyleEnvironment->"Working", WindowToolbars->"EditBar", WindowSize->{630, 560}, WindowMargins->{{38, Automatic}, {6, Automatic}}, CellLabelAutoDelete->False, StyleDefinitions -> "ArticleClassic.nb", MacintoshSystemPageSetup->"\<\ 00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001 0000I00000400`<300000BL?00400@0000000000000006P801T1T00000000000 00000000000000000000000000000000\>" ] (*********************************************************************** Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. ***********************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1739, 51, 85, 1, 105, "Title"], Cell[1827, 54, 60, 1, 40, "Author"], Cell[1890, 57, 357, 7, 82, "Abstract"], Cell[2250, 66, 33, 11, 150, "Text"], Cell[CellGroupData[{ Cell[2308, 81, 56, 1, 57, "Section"], Cell[2367, 84, 1218, 19, 193, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[3622, 108, 73, 1, 31, "Section"], Cell[3698, 111, 5724, 153, 470, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[9459, 269, 95, 1, 31, "Section"], Cell[9557, 272, 2730, 82, 203, "Text"], Cell[12290, 356, 169, 7, 73, "Input"], Cell[12462, 365, 79, 2, 28, "Input"], Cell[12544, 369, 77, 2, 28, "Input"], Cell[12624, 373, 189, 6, 25, "Text"], Cell[12816, 381, 252, 9, 103, "Input"], Cell[13071, 392, 437, 11, 55, "Text"], Cell[13511, 405, 86, 2, 28, "Input"], Cell[13600, 409, 149, 5, 43, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[13786, 419, 88, 1, 31, "Section"], Cell[13877, 422, 605, 12, 82, "Text"], Cell[14485, 436, 313, 11, 133, "Input"], Cell[14801, 449, 1229, 29, 156, "Text"], Cell[16033, 480, 195, 7, 73, "Input"], Cell[16231, 489, 391, 12, 55, "Text"], Cell[16625, 503, 246, 8, 88, "Input"], Cell[16874, 513, 330, 11, 40, "Text"], Cell[17207, 526, 105, 2, 28, "Input"], Cell[17315, 530, 214, 7, 39, "Text"], Cell[17532, 539, 109, 2, 28, "Input"], Cell[17644, 543, 221, 7, 39, "Text"], Cell[17868, 552, 262, 8, 88, "Input"], Cell[18133, 562, 433, 12, 55, "Text"], Cell[18569, 576, 392, 12, 148, "Input"], Cell[CellGroupData[{ Cell[18986, 592, 71, 1, 42, "Subsection"], Cell[19060, 595, 624, 16, 83, "Text"], Cell[19687, 613, 455, 14, 178, "Input"], Cell[20145, 629, 394, 12, 54, "Text"], Cell[20542, 643, 86, 2, 28, "Input"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[20677, 651, 62, 1, 31, "Section"], Cell[20742, 654, 887, 20, 95, "Text"], Cell[21632, 676, 422, 12, 148, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[22091, 693, 84, 1, 31, "Section"], Cell[22178, 696, 490, 11, 68, "Text"], Cell[22671, 709, 229, 7, 73, "Input"], Cell[22903, 718, 261, 10, 25, "Text"], Cell[23167, 730, 162, 6, 58, "Input"], Cell[23332, 738, 78, 2, 28, "Input"], Cell[23413, 742, 78, 1, 24, "Text"], Cell[23494, 745, 462, 14, 178, "Input"], Cell[23959, 761, 184, 5, 39, "Text"], Cell[24146, 768, 473, 10, 174, "Input"], Cell[24622, 780, 726, 14, 111, "Text"], Cell[25351, 796, 318, 12, 148, "Input"], Cell[25672, 810, 622, 19, 70, "Text"], Cell[26297, 831, 244, 4, 94, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[26578, 840, 66, 1, 31, "Section"], Cell[26647, 843, 582, 11, 94, "Text"], Cell[27232, 856, 93, 2, 28, "Input"], Cell[27328, 860, 237, 8, 24, "Text"], Cell[27568, 870, 116, 2, 28, "Input"], Cell[27687, 874, 1051, 29, 113, "Text"], Cell[28741, 905, 189, 6, 58, "Input"], Cell[28933, 913, 121, 4, 24, "Text"], Cell[29057, 919, 205, 7, 73, "Input"], Cell[29265, 928, 956, 23, 113, "Text"], Cell[30224, 953, 367, 12, 148, "Input"], Cell[30594, 967, 390, 9, 52, "Text"], Cell[30987, 978, 234, 8, 103, "Input"], Cell[31224, 988, 101, 1, 24, "Text"], Cell[31328, 991, 168, 7, 73, "Input"], Cell[31499, 1000, 506, 15, 55, "Text"], Cell[32008, 1017, 114, 2, 28, "Input"], Cell[32125, 1021, 335, 8, 53, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[32497, 1034, 54, 1, 31, "Section"], Cell[32554, 1037, 225, 6, 30, "Reference"], Cell[32782, 1045, 192, 6, 18, "Reference"], Cell[32977, 1053, 203, 6, 30, "Reference"], Cell[33183, 1061, 242, 7, 32, "Reference"], Cell[33428, 1070, 200, 6, 18, "Reference"], Cell[33631, 1078, 207, 6, 30, "Reference"], Cell[33841, 1086, 199, 6, 18, "Reference"] }, Closed]], Cell[CellGroupData[{ Cell[34077, 1097, 35, 0, 31, "Section"], Cell[34115, 1099, 559, 13, 82, "Text"] }, Closed]], Cell[34689, 1115, 222, 9, 118, "Address"] }, Open ]] } ] *) (*********************************************************************** End of Mathematica Notebook file. ***********************************************************************)