# Why is this Mandelbrot set's implementation infeasible: takes a massive amount of time to do?

The Mandelbrot set is defined by complex numbers such as $z=z^2+c$ where $z_0=0$ for the initial point and $c\in\mathbb C$. The numbers grow very fast in the iteration.

z = 0; n = 0; l = {0};
While[n < 9, c = 1 + I; l = Join[l, {z}]; z = z^2 + c; n++];l


If n is very large, the numbers become too large and impossible to calculate in practical time limits. I don't know what it would look after long time but doubting whether it would look like here.

What is wrong with this implementation? Why does it take so long time to calculate?

-
somewhat relevant discussion in chat –  acl Mar 20 '13 at 10:35
It takes a long time to calculate because you are doing exact arithmetic with very large numbers. You will want to use machine precision for computing fractals in any reasonable time. More importantly, why are you computing the sequence for an escaping point? I think you need to read more and code less until you have a basic understanding of what you're trying to do. –  Simon Woods Mar 20 '13 at 10:45
@SimonWoods is there some command to make Mathematica use only Machine-precision numbers? –  hhh Mar 20 '13 at 11:15
A useful introduction wikihow.com/Plot-the-Mandelbrot-Set-By-Hand –  cormullion Mar 20 '13 at 11:21
Add a decimal point to the numbers you use, or use N. –  Sjoerd C. de Vries Mar 20 '13 at 11:37

What is wrong: a) you're using exact arithmetic. b) You keep iterating even if the point seems to be escaping.

Try this

ClearAll@prodOrb;
prodOrb[c_, maxIters_: 100, escapeRadius_: 1] :=
NestWhileList[#^2 + c &,
0.,
Abs[#] < escapeRadius &,
1,
maxIters
]

prodOrb[0. + 10. I]
prodOrb[0. + .1 I]


(if you don't need the entire list but only the final point, replace NestWhileList by NestWhile).

Here, I use approximate numbers by using 0. rather than 0. See this tutorial for more.

EDIT: Since we're doing interactive manipulation:

ClearAll[mnd];
mnd = Compile[{{maxiter, _Integer}, {zinit, _Complex}, {dt, _Real}},
Module[{z, c, iters},
Table[
z = zinit;
c = cr + I*ci;
iters = 0.;
While[(iters < maxiter) && (Abs@z < 2),
iters++;
z = z^2 + c
];
Sqrt[iters/maxiter],
{cr, -2, 2, dt}, {ci, -2, 2, dt}
]
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];

Manipulate[
lst = mnd[100, {1., 1.*I}.p/500, .01];
ArrayPlot[Abs@lst],
{{p, {250, 250}}, Locator}
]


Clicking around changes the fractal. Note the magic numbers sprinkled throughout the code. Why? Because ListContourPlot is way too slow, so that using the coords of the clicked point ended up being too much of a waste of time (and my coffee break is over).

EDIT2: So much for the break being over. Here we have the Mandelbrot set being blown away by strong winds:

tbl = Table[
lst = mnd[100, (1 + 1.*I)*p/500, .01];
ArrayPlot[Abs@lst],
{p, 0, 500, 10}
];

ListAnimate[tbl]


And see, here it is, sliding off the table while being melted:

tbl2 = Table[
mnd[100, (1 + 1.*I)*p/500, .05] // Abs //
ListPlot3D[#, PlotRange -> {0, 1},
ColorFunction -> "BlueGreenYellow", Axes -> False,
Boxed -> False,
ViewVertical -> {0, (p/500), Sqrt[1 - (p/500)^2]}] &,
{p, 0, 500, 25}
];


(this is very slow, because ListPlot3D is very slow)

Maximal silliness has now been achieved. Or has it?

-
why do you the dot "0." in the term "NestWhileList[#^2+c &,0.,Abs[\$]<escapeRadius&," and not "0"? –  hhh Mar 20 '13 at 11:30
@hhh: "0" is an exact number. "0." is a machine-precision number, and calculations with machine-precision numbers are much faster. –  nikie Mar 20 '13 at 11:37
@nikie is correct; press F1 to bring up the help centre, then paste tutorial/ExactAndApproximateResults into the top bar to see this tutorial. –  acl Mar 20 '13 at 11:50
After seeing the chat you referred to, +1 for patience. –  Mark McClure Mar 20 '13 at 13:17
@MarkMcClure Indeed (and more, if you include the previous day)! –  rm -rf Mar 20 '13 at 15:16

Just a bit of fun with @acl's code:

ArrayPlot[Table[
NestWhile[#^2 - (0. - 1 I) & , r + i I, Abs[#] < 2.0 &, 1, 10],
{r, -2, 2, 0.005},
{i, -2, 2, 0.005}]]


-
+1 I swear I didn't even know that ArrayPlot accepted complex values! –  Mark McClure Mar 20 '13 at 13:22
Hehe.. look at what you implemented and what the OP was asking :-) And don't take my implementation as reference, because I was only looking at yours and did it as wrong. Let's see who realizes it. –  halirutan Mar 20 '13 at 14:58
@halirutan What was the OP asking? –  Mark McClure Mar 20 '13 at 15:02
@halirutan well, i wasn't really answering the OP - just thought I'd post a nice picture I made while tinkering with acl's code (a cast silver lightning bolt wrapped in an infinite amount of black velvet). –  cormullion Mar 20 '13 at 15:03
Sweet! Can it dance? Perhaps like the eyes? :D –  rm -rf Mar 20 '13 at 15:17

With a compiled version you get it so fast, that you can manipulate it in real time.

fc = Compile[{{in, _Complex, 0}, {c, _Complex, 0}},
Module[{iter = 0, max = 10, z = in},
While[iter++ < max,
If[Abs[z = z^2 + c] > 2.0,
Break[]
]
];
{Abs[z], iter}
], CompilationTarget -> "C", Parallelization -> True,
RuntimeAttributes -> {Listable}
];

data = Table[i + I j, {j, -2, 2, .01}, {i, -2, 2, .01}];

Manipulate[ImageAdjust[Image[fc[data, cx + I cy]]],
{cx, -.5, .5},
{cy, -.5, .5}
]


And since my iteration does calculate a Julia set instead of the Mandelbrot set, let's fix this

fc = Compile[{{c, _Complex, 0}},
Module[{iter = 0, max = 100, z = 0 + 0 I},
While[iter++ < max,
If[Abs[z = z^2 + c] > 2.0,
Break[]
]
];
{Log[iter], Abs[z]}
], CompilationTarget -> "C", Parallelization -> True,
RuntimeAttributes -> {Listable}
];

data = Table[i + I j, {j, -2, 2, 4/1023.}, {i, -2, 2, 4/1023.}];


To calculate this 1024x1024 image, it took only 0.2 seconds on my machine here. I wouldn't consider this slow.

-
My version is a little faster:) –  chyaong Mar 20 '13 at 15:23
Put RuntimeOptions->"Speed" in my code and watch yourself ;-) –  halirutan Mar 20 '13 at 15:28
I tried this, my code about three times faster. –  chyaong Mar 20 '13 at 15:39
@chyanog On my machine, the timing is 0.169407 for your method and 0.163131 for my one with option "Speed". –  halirutan Mar 20 '13 at 15:49
Oh, Something strange about this, the reason might be that the difference C Compiler. See mathematica-project.googlecode.com/files/2013-03-21_000830.jpg –  chyaong Mar 20 '13 at 16:11

As @acl mentioned in chat, your question really indicates that you should read some fundamental sources. Two that I'd recommend are:

Note that the second reference, which I happen to have written largely avoids Compile and the speediest of speedy tricks presented by the other answers here in favor of simpler code. This is probably better, if your objective is to learn the basics of the material.

-
+1 I agree that avoiding C code is important here. I think all the other answers (including mine, which started out as a bit of a joke and now essentially parodies how mma.se works) are just noise, and do not help hhh to learn what they need to learn. Yours is probably the most useful answer. Lacking graphics as it does, though, it'll get less upvotes! –  acl Mar 20 '13 at 16:19
@acl Thanks! I can live without the upvotes; I tend to be conservative with mine, anyway. –  Mark McClure Mar 20 '13 at 16:33
they are practically useful in that the visibility of upvoted answers increases. However, in this particular case, the question really is framed in too specific terms (it's really 1. how do I use floating-point numbers 2. what's wrong with my algorithm to obtain a mandelbrot set), so unlikely to be useful to future visitors with similar questions. Well maybe it'll be useful to people who are new to complex numbers and fractals, to whom your answer is the only useful one... –  acl Mar 20 '13 at 17:32
Where can you find the notebooks for the book Mathematica in Action by Stan Wagon? –  hhh Mar 20 '13 at 19:33
@hhh They come with the book! An antiquated concept, I know. –  Mark McClure Mar 20 '13 at 20:47

In Mathematica 10.0:

    MandelbrotSetPlot[{-2 - I, 1 + I}]


(Admittedly, this doesn't address what's wrong with the code in the original question, but if you just want to get a Mandelbrot set plot, surely a built-in function is likely to be reasonably efficient.)

-
cf = Compile[{{a, _Real, 1}},
Module[{z = 0 I, i = 0, max = 100},
While[i++ < max && Abs[z] <= 2,
z = z^2 + a[[1]] + a[[2]]*I]; {Log[i], Abs[z]}],
RuntimeAttributes -> {Listable}, CompilationTarget -> "C",
RuntimeOptions -> "Speed", Parallelization -> True];

t1 = AbsoluteTime[];
data = Table[{i, j}, {j, -2, 2, 4/1023.}, {i, -2, 2, 4/1023.}];