Friday, February 3, 2012

better late than never: Mandelbrot Set


Although I have a math and computer science background, this (above) is my first attempt to draw the Mandelbrot set. It seems overdue since plotting the set is such a well established computer-math project that it's almost cliche (there is so much online about this, but this Math Munch article has some nice pointers). But if I am just finally getting around to this, then it's not too late for you too.

So, if you haven't written a little program to draw the set yet, I strongly recommend it. There are lots of nice things to think about as you explore this, and very little of it has much to do with fractals - although having that strange inkblot-like image appear at the end is the carrot, or cauliflower, that will motivate you along.

I'll share my little program and explanations below. You can find or create much more slick and efficient code if you try. For example, you don't really need to know the magnitude of the complex points (the square of the magnitude is enough), and I'm sure my point generation isn't the best.

what are you plotting?
I would guess that most students on the pre-calc train tend to think of plotting points almost exclusively in terms of single variable functions of real numbers. This is not one of those kinds of plots. Instead of (x,y) indicating y = (x), we are plotting points belonging to C, and making them one colour if they belong to a particular set M, and a different colour if they don't.

To figure out whether a point is in M, for every point that you consider, you need to construct a particular sequence - if that sequence stays bounded, then it is in M, but if the sequence is unbounded it isn't in M. The further you go in calculating your sequence, the more sure you are that your points are actually in the set.

some complex number stuff
You don't need to write much code to implement the complex number operations that you need for this, but I'm partial to encapsulating this sort of thing in a general purpose class like this Processing example:

class CPoint {
 float xvalue;
 float yvalue;
 float pointsize = 0.2;
 CPoint(float x, float y){
  xvalue = x;
  yvalue = y;
 } 
 CPoint sum(CPoint p){  
  return new CPoint(xvalue +p.xvalue, yvalue+ p.yvalue);
 }
 CPoint prod(CPoint p){
  float newx = xvalue*(p.xvalue) - yvalue*(p.yvalue);
  float newy = xvalue*(p.yvalue) + yvalue*(p.xvalue); 
  return new CPoint(newx, newy);
 }
 CPoint squared(){
   return this.prod(this);
 } 
 float magnit(){
   float part = pow(xvalue,2) + pow(yvalue,2);
   return sqrt(part);
 }
 void display(float xshift, float yshift, float zoom) {
    fill(255,255);
    ellipse(xvalue*zoom + xshift, yvalue*zoom  + yshift, pointsize, pointsize);
  }  
}

The key item is the complex multiplication - this is what distinguishes a point in C from one in R^2.

I have a like/dislike relationship with Processing - I like how quickly things can be created and that it makes  nice pictures without much effort, but I dislike how I can't seem to help breaking fundamental programming rules when using it (I end up using global variables, and always break model-view separation) - likely a personal problem, rather than a problem with Processing itself.

determining if a point is in the set M
You are going to be finding a sequence of points based on a special rule - if applying the special rule repeatedly causes the points to get too big, then they are out of the set. In addition to knowing the special rule, you also need to know "how big is too big" and "how many times will I apply the rule."

The calculation is encapsulated in this other Processing class Map. Can you figure out what the rule is, and how you specify the "how big" and "how many" parts of the calculation?

class Map {
 CPoint c;
 CPoint first;
 CPoint current;

  Map(CPoint initial, CPoint cValue){
   first = initial;
   current = initial;
   c = cValue;
  }
    
  void iterate(){
    current = current.squared().sum(c);
  }
  
  boolean iterate(int iterations, int bound){    
    for(int i=0; i< iterations; i++){
      this.iterate();    
      if(current.magnit()>bound) break;
    }
    return (current.magnit() < bound);
  }
  
  void display(float xshift, float yshift, float zoom){
   c.display(xshift, yshift, zoom);
  }    
}

If you only calculate for a few iterations, some points that should not be in the set get included - after 5 iterations, you get something that looks like a skate.


Having a somewhat less precise plot gives you funkier looking pictures than the crisper image that you get with more iterations. For this reason, people often include points not quite in the set in their pictures of the Mandelbrot set, and often use colour to show at which iteration a given point failed the membership test. These fuzzy pictures are really good examples of fuzzy sets - where instead of set membership being true or false, it is a range of values indicating by what level the point failed to be in the actual set. The picture below was generated using 20 iterations.



plotting the set
The last thing you need to do is to set up your window, and generate your points. Here's the main Processing file that I used for this (no zoom or panning in this one - exercise left to the reader):

//various magic numbers
int windowX = 600;
int windowY = 500;
int iterations = 200;
int zoom = 250;
int disk = 2;
int numberPoints = 100000;
int randomRange = 100;
float centerX = -0.5;
float centerY = 0;

//init
void setup() {
  size(windowX, windowY);
  noStroke();
  smooth();
  background(0);
}

void draw() 
{
  loop();
  CPoint zero = new CPoint(0,0);
  float x;
  float y;
  int signx = -1;
  int signy = -1;
  Map newMap;
  CPoint cPoint;  
  for( int i=0; i < numberPoints; i++){
    cPoint = randomCPoint(1.5, randomRange);
    newMap = new Map(zero,cPoint);
    if(newMap.iterate(iterations, disk)){
        newMap.display(windowX/2 - centerX*zoom, windowY/2 - centerY*zoom, zoom );
    }  
  }  
}

CPoint randomCPoint(float bound, float depth){
  float x;
  float y;
  int signx = -1;
  int signy = -1;
  x = bound*random(depth)/depth;
  y = bound*random(depth)/depth;
  if(random(100)<50) signx *= (-1);
  if(random(100)<50) signy *= (-1);
  return new CPoint(signx*x,signy*y); 
 }

If you implement it like this, the set will slowly emerge as more and more points are tested. This crisper version below used 200 iterations (maybe a bit excessive).



Update: A Javascript version of this is posted here.