# 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, 2013 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. Mar 20, 2013 at 10:45
• @SimonWoods is there some command to make Mathematica use only Machine-precision numbers?
– hhh
Mar 20, 2013 at 11:15
• A useful introduction wikihow.com/Plot-the-Mandelbrot-Set-By-Hand Mar 20, 2013 at 11:21
• Add a decimal point to the numbers you use, or use N. Mar 20, 2013 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.,
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, 2013 at 11:30
• @hhh: "0" is an exact number. "0." is a machine-precision number, and calculations with machine-precision numbers are much faster. Mar 20, 2013 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, 2013 at 11:50
• After seeing the chat you referred to, +1 for patience. Mar 20, 2013 at 13:17
• @MarkMcClure Indeed (and more, if you include the previous day)!
– rm -rf
Mar 20, 2013 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! Mar 20, 2013 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. Mar 20, 2013 at 14:58
• @halirutan What was the OP asking? Mar 20, 2013 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). Mar 20, 2013 at 15:03
• Sweet! Can it dance? Perhaps like the eyes? :D
– rm -rf
Mar 20, 2013 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}];

{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.

• Put RuntimeOptions->"Speed" in my code and watch yourself ;-) Mar 20, 2013 at 15:28
• @chyanog On my machine, the timing is 0.169407 for your method and 0.163131 for my one with option "Speed". Mar 20, 2013 at 15:49
• Here chyanog's is about three times faster, RuntimeOptions -> "Speed" included: i.stack.imgur.com/RtHDG.png OS X, clang-425.0.27 Mar 20, 2013 at 16:39
• @chyanog Actually both run in the same amount of time. The mistake I made was to use your data with halirutan's function. Maybe you made the same mistake. Both run in about 0.17 seconds here. Mar 20, 2013 at 16:47

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, 2013 at 16:19
• @acl Thanks! I can live without the upvotes; I tend to be conservative with mine, anyway. Mar 20, 2013 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, 2013 at 17:32
• Where can you find the notebooks for the book Mathematica in Action by Stan Wagon?
– hhh
Mar 20, 2013 at 19:33
• @hhh They come with the book! An antiquated concept, I know. Mar 20, 2013 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.)

ClearAll[mandel];
mandel =
Compile[{{A, _Real, 1}, {maxIt, _Real}},
Module[{i = 0, re = 0., im = 0., temp, x = A[], y = A[]},
While[i++ < maxIt && re^2 + im^2 <= 4,
temp = re^2 - im^2 + x;
im = 2 re im + y;
re = temp;
];
{Log[i], Sqrt[re^2 + im^2]}
],
CompilationTarget -> "C", RuntimeOptions -> "Speed", RuntimeAttributes -> {Listable}
];

t1 = AbsoluteTime[];
data = mandel[Table[{i, j}, {j, -2, 2, 4/1023.}, {i, -2, 2, 4/1023.}], 100];
Image@Rescale@data
AbsoluteTime[] - t1 