Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
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 – 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
up vote 21 down vote accepted

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

Try this

prodOrb[c_, maxIters_: 100, escapeRadius_: 1] := 
 NestWhileList[#^2 + c &,
  Abs[#] < escapeRadius &,

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:

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

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

Mathematica graphics

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];
  {p, 0, 500, 10}


Withering mandelbrot

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)

enter image description here

Maximal silliness has now been achieved. Or has it?

share|improve this answer
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)! – R. M. Mar 20 '13 at 15:16

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

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


share|improve this answer
+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 – R. M. 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,
    {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}

Mathematica graphics

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,
    {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.}];

Mathematica graphics

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

share|improve this answer
Put RuntimeOptions->"Speed" in my code and watch yourself ;-) – halirutan Mar 20 '13 at 15:28
@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 – chyaong Mar 20 '13 at 16:11
Here chyanog's is about three times faster, RuntimeOptions -> "Speed" included: OS X, clang-425.0.27 – Szabolcs Mar 20 '13 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. – Szabolcs Mar 20 '13 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.

share|improve this answer
+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 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
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.}];
AbsoluteTime[] - t1

enter image description here

share|improve this answer

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

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.