Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
2  
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
1  
A useful introduction wikihow.com/Plot-the-Mandelbrot-Set-By-Hand –  cormullion Mar 20 '13 at 11:21
1  
Add a decimal point to the numbers you use, or use N. –  Sjoerd C. de Vries Mar 20 '13 at 11:37
add comment

closed as too localized by acl, Yves Klett, rcollyer, Sjoerd C. de Vries, Mark McClure Mar 20 '13 at 17:41

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

5 Answers

up vote 15 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

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}
 ]

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];
  ArrayPlot[Abs@lst],
  {p, 0, 500, 10}
  ];

ListAnimate[tbl]

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
2  
@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
2  
@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
2  
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
show 2 more comments

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}]]

fractal

share|improve this answer
4  
+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
show 2 more comments

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}
]

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,
      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.}];
ImageAdjust[Image[fc[data]]]

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
    
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
show 3 more comments
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.}];
ImageAdjust[Image[cf[data]]]
AbsoluteTime[] - t1

enter image description here

share|improve this answer
add comment

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 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
1  
@hhh They come with the book! An antiquated concept, I know. –  Mark McClure Mar 20 '13 at 20:47
add comment

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