Tools for Analyzing Discrete Dynamical Systems:
Complex Dynamics
Written by Richard A. Holmgren
last modified July 20, 1996

These programs were developed for use with the text  A First Course in Discrete Dynamical Systems, 2nd edition by Richard A. Holmgren,  which was published by Springer-Verlag in 1996.   Unfortunately, the programs I have written for complex dynamics are not as fast as one would like.  High-resolution pictures often require an hour or more to generate, though in most cases reasonably good pictures can be produced in under five minutes.  Users requiring faster results should note that computational time would probably be decreased by at least an order of magnitude if MathLink  were used to do the numerical calculations. The use of MathLink is described in the reference Mathematica by Wolfram and the text Applied Mathematica by Shaw and Tigg.  Shaw and Tigg also have placed a notebook with sample code for creating Mandelbrot Sets on this website, as has Terry Robb.  (The notebook reference numbers are 0205-366 and 0204-107.  The Mandelbrot set code is in Chapter 19 of the text by Shaw and Tigg.)   
	 Suggestions for improving this notebook are gratefully accepted.  Please contact Professor Holmgren at rholmgre@alleg.edu with suggestions and comments.

Julia Sets

Description

The first block of code defines the compiled function julia[x,y, c,m].  julia[...] determines whether or not the point x + I y, leaves the circle of radius |c| + 1 before m iterations of q(z) = z2 + c have been completed.  If so, julia[...] returns the value 1.  If not, julia[...] returns the value 0.  Recall that a complex number is in the Julia set of c  if and only if the orbit of the number under iteration of  q(z) = z2 + c is bounded by  |c| + 1.  Thus points for which julia[...]  returns the value 1 are probably not in the Julia set.  (Note however the caveat in the next paragraph.)	
	We then use julia  to sketch julia sets for c = -1 and c = I.  Notice that the Julia set for 
c = I is nowhere dense and repelling.  Thus, if we start at any point (even one in the Julia set) and iterate q(z) = z2 + I, then rounding will soon cause the iterate to leave the Julia set and grow larger than  2.  Hence, in graphing the Julia set for c = I we only use 12 iterates.  On the other hand, we seem to get a better picture of the Julia set for c = -1 by using a larger number of iterates.  It is instructiveto ponder whether or not this is indeed the case.

Code for julia[...]

Clear[julia,x,y,c,m,z,n,k,bound];
julia = 
	Compile[{x,y,{c, _Complex},{m, _Integer}},
	Module[{z = x+I y, n = 0, k = 0},
		bound = Abs[c] + 1;
		While[n <= m && k == 0,
			If[Abs[z] > bound , k = 1,z = z^2 + c;n++]];
			k],
	{{z,_Complex},{n,_Integer},{k,_Integer}}];

Julia set for c = -1

Timing[
DensityPlot[julia[x,y,-1,100],
	{x,-2,2},{y,-2,2},
	PlotPoints -> 100, Mesh -> False]
	]

Timing[
DensityPlot[julia[x,y,-1,12],
	{x,-2,2},{y,-2,2},
	PlotPoints -> 100, Mesh -> False]
	]

Julia set for c = I

Timing[
DensityPlot[Re[julia[x,y,I,12]],
	{x,-2,2},{y,-2,2},
	PlotPoints ->100, Mesh -> False]
	]

The Mandelbrot Set

Description

The following code will draw the Mandelbrot set.  The first block of code defines the compiled function mandel[...] which takes three variables, x, y, and m.  The first two are assumed to be real numbers and the third should be an integer.  Mandel determines whether or not the point 0 leaves the circle with center 0 and radius 2 within m iterations of the complex quadratic function q(z) = z2 + x + I y.  If it does, then mandel[...] returns the value 1.  If  it  remains in the circle for m iterations, then mandel[...] returns the value 0.  
	The second block of code uses mandel[...] and DensityPlot[...] to draw the Mandelbrot set.  We assume that points for which mandel[...] returns 0 are in the Mandelbrot set.  
	A more interesting picture can be created by creating a function which returns different values, depending on how many iterations it takes before the initial value escapes the circle with center 0 and radius 2.   For example below we write code for the function mandelLevels[...], which returns the value 4 if the initial point escapes within 10 iterations the function. If the point escapes in 10-19 iterations mandelLevels[...] returns 3, if it escapes in 20-29 iterations mandelLevels[...] returns 2, and  if it escapes in 30-39 iterations mandelLevels[...] returns 1.  If the point doesn't escape within 40 iteraions, then mandelLevels[...] returns 0.  Note that  unlike mandel[...], mandelLevels[...] only takes two variables, the real and imaginary parts of the initail point.

The Code for mandel[...]

Clear[mandel,x,y,m,w,c,n,k];
mandel = 
	Compile[{{x,_Real},{y,_Real},{m, _Integer}},
		Module[{w = 0*I,c = x + I y,n = 0, k = 0},
			w = x + I y;
			While[n <= m && k == 0,
				If[Abs[w] > 2, k=1, w = w^2 + c;n++]];
			k],
		{{k,_Integer},{n,_Integer},{w,_Complex},{c,_Complex}}
		];

The Code for mandelLevels[...]

Clear[mandelLevels,x,y,w,c,n,k];
mandelLevels = 
	Compile[{{x,_Real},{y,_Real}},
		Module[{w = 0*I,c = x + I y,n = 0, k = 0},
			w = x + I y;
			While[n < 10 && k == 0,
				If[Abs[w] > 2, k=4, w = w^2 + c;n++]];
			While[n < 20 && k == 0,
				If[Abs[w] > 2, k=3, w = w^2 + c;n++]];
			While[n < 30 && k == 0,
				If[Abs[w] > 2, k=2, w = w^2 + c;n++]];
			While[n < 40 && k == 0,
				If[Abs[w] > 2, k=1, w = w^2 + c;n++]];
			k],
		{{k,_Integer},{n,_Integer},{w,_Complex},{c,_Complex}}
		];

Drawing the Mandelbrot Set

The role of the variables defined in the first 6 lines should be readily apparent.  Care should be taken when choosing the grid size.  The time necessary for creating the image is proportional to the  square of the grid size.

minRealPart = -2	;
maxRealPart = 2		;
minImPart = -2		;
maxImPart =	2		;

maxIterations = 50	;

gridSize = 100 ;

Timing[DensityPlot[mandel[x,y,maxIterations],
	{x,minRealPart,maxRealPart},{y,minImPart,maxImPart},
	PlotPoints -> gridSize, Mesh -> False]]

Drawing the Mandelbrot Set with mandelLevels[...]

The role of the variables defined in the first 5 lines should be readily apparent.  Care should be taken when choosing the grid size.  The time necessary for creating the image is proportional to the  square of the grid size.  In this graphic we exmanine a portion of the boundary of the Mandelbrot set.

minRealPart = -1.3	;
maxRealPart = -1		;
minImPart = .1		;
maxImPart =	.4		;

gridSize = 100 ;

Timing[DensityPlot[mandelLevels[x,y],
	{x,minRealPart,maxRealPart},{y,minImPart,maxImPart},
	PlotPoints -> gridSize, Mesh -> False]]

Newton's Method for Complex-valued Functions

Description

There are two blocks of code.  The first defines the function newton[...], which takes two variables, a complex number and an integer.  newton[...] will iterate starting at the given complex number and using newton's method for p(z) = z3 -1 for up to the given integer number of iterations.  If during that time the iterations are within .25 units of one of the roots of p(z), then newton[...] will return an integer between 1 and 3, the choice depending on which root is being approached.  If after iteration the inital value has not settled on a root, then newton returns the integer 4.  This code can easily be modified for other polynomials.  Note that in the While[...] loop we test to see if the iteration has reached 0, since that is not in the domain of Newton's method for p(z).   Similar tests would need to be constructed for other polynomials.  Notice also that to write a program like this we must first find a good approximation of at least some of the roots.
	The second block of code uses newton[...] and DensityPlot[...]  to create an image of the basins of attraction.

The Code for newton[...]

Clear[newton,z,m,k,w,n];
newton = 
	Compile[{{z, _Complex},{m, _Integer}},
		Module[{w = z, n = 0, k = 4},
			While[n <= m && k == 4 && w != 0.0 + 0.0 I,
				If[Abs[w - (1.0+0. I)]<=.25,k = 1,
					If[Abs[w - (-.5+.866I)]<=.25,k = 2,
						If[Abs[w - (-.5-.866I)]<=.25,k = 3]]];
				w = (2 w^3 + 1)/(3 w^2);
				n++];
			k],
		{{w,_Complex},{n,_Integer},{k,_Integer}}
		];

Drawing the basins of attraction

Timing[
DensityPlot[newton[ x + I y,50],
	{x,-2,2},{y,-2,2},
	PlotPoints -> 40, Mesh -> False]
]