Mandelbrot Fractal v2

I had so much fun with the original Mandelbrot program that I decided to see if I could improve it. I was amazed by how easily I was able to generate a fractal picture that looked exactly like other pictures online. I didn’t look at anyone else’s code, just went down the path and the amazing Mandelbrot just happened.

Before You Read, Play

Click here to play with it in your web browser.

Ok, now that you’re back: here are some new “features” that were implemented, and the interesting math that happened:

Zooming

Left click to zoom in on that point (70% smaller window). Right click to zoom out ( 1/0.7 = 143% bigger window).

The first problem that I had with zooming was the conflict of the virtual coordinate complex plane (-2.0<x<2.0 and -2.0<y<2.0) and the processing.org window (0<width<500 pixels and 0<height<500 pixels, with the height pixels arranged from top to bottom (instead of bottom to top).

2014-04-04_14h13_02

After much brooding, lots of print(variable_name_here)s, lots of sketching on paper, and finally changing a + to a - or - to a +, the zooming finally worked as desired. The code is more ugly than I’d normally like, but it’s generalized so you can change the resolution to whatever you’d like and it should work.
Here’s the processing code that takes a 1D pixel array (0 to 250,000) and produces each complex coordinate:

  loadPixels();
  color c;
  float windowwidth = maxwindowx - minwindowx;
  float windowheight = maxwindowy - minwindowy;
  for (int i = 0; i < pixels.length; i++)
  {
    //this code maps the one dimensional pixel array to a cx and cy in the complex plane
    float startcx = windowwidth * (float(i % width) / width) - (windowwidth/2) + centerx; 
    float startcy = windowheight * (float(floor(float(i)/width))/height) - (windowheight/2) + centery;

The interesting math to get 2D information from a 1D array is here:

i % width
floor(i / width)

The first line gets the row (the % sign means remainder, so the 1052nd pixel is 1052 % 500 = 52, so it’s on row 52 (of 500)).
The second line gets the column (divides i and width, and rounds down, so floor(1052 / 500) = 2, so it’s on column 2.

Adaptive Resolution (and Color)

The resolution is how the program calculates the color of each pixel.

Some background: the mandelbrot set is the set of all coordinates in the complex plane (x is reals, y is imaginary), c=a+bi, where if you start with the complex number z=0+0i, and iterate the following operation:
z := z^2 + c
The resultant complex number stays with a size (modulus) of less than or equal to 2. So in otherwords, if c = 0+0i, then z^2+c will always have a size of less than 2. Color this point black.

On the otherhand, if I start out with c=1+1i then I get the following sequence for z:

z_1 = (0+0i)^2+c = 1+1i  size <= 2, keep going.
z_2 = (1+1i)^2+c = 2i  size <= 2, keep going.
z_3 = (2i)^2+c = 1+3i  size > 2, stop.

So the complex point 1+1i fails after 3 steps. Color this point something dark because it fails somewhat quickly. If you make the rule that the longer a point goes before it fails, the lighter the color then you get something exactly like the above images. The coloring is that easy.

Here’s the coloring code:

if  (mandelvalue < 2)
    {
        c = color(0,100,0);
    }
    else
    {
        //the color of the pixel depends on how far the iteration got before it was out of bounds
        float ratio = 1.0*count/resolution;
        c = color(10+1.0*ratio*140,100,1.0*ratio*200);
    }

If the final z (mandelvalue) is less than 2, then color it black, otherwise color it based on how many cycles it took to exceed two (count). Note: this color system is based on 0 to 100 HSB color format.

A big problem that has to be dealt with is when do you stop checking to see if a point converges to <= 2 or diverges….

Here’s a video of the resolution increasing one step at a time. A resolution of 1 means that all points whose size is <= 2 after 1 step are colored black. A resolution of 50 means that all points whose size is <= 2 after 50 steps are colored black. If the point strays outside of a size of 2, then color the point dependent on how quickly it strayed.

Ok. So clearly more resolution is better right? The problem is that it’s starting to get computationally slow. 500×500 is 250,000 pixels and at a resolution of 50 each one of these pixels may need up to 50 sets of calculations.

If we’re able to zoom in, then how can we make a clear enough picture for the given window size? If the picture is really zoomed in, then we’re going to need a higher resolution to get a nice picture.

So what I need is a way to dynamically change the resolution based on the virtual window size. If the virtual window size has a width of 4 units (-2 to 2), then a resolution of 18 is pretty good. But if the virtual window size has a width of 0.2 (-1.5 to -1.3) then I need a higher resolution. I need an equation to set the resolution based on the window size. Enter regressions. I took some points and made a quadratic regression:

The x-axis is the windowwidth, and the y-axis is the resolution. The red graph was the first regression that I used. Not great. As the window got smaller and smaller, the increase in the resolution necessary didn’t keep up. So I went to a log regression (duh, vertical asymptote at x=0), and it picked a resolution with more accuracy.
Here’s the resulting code:

resolution = int(27.71-5.156*log(windowwidth));

Varieties of Mandelbrot

z^2+c is only one possibility.
You can try z^n+c as well. Here’s a video walk through of z^2+c, z^3+c, z^4+c, … , z^10 + c. You can try this yourself in the real program by clicking the 2, 3, 4, 5, … , 9, and 0 (10) keys:


Here’s the code for the z^2 + c and the z^3 + c using binomial expansion that was done by hand (and simplification of powers of i).

  if (power < 2.01 && power > 1.99)
  {
    fx = x*x - y*y + startcx;
    fy = 2*x*y + startcy;
  }
  else if (power < 3.01 && power > 2.99)
  {
    fx = x*x*x-3*x*y*y + startcx;
    fy = 3*x*x*y-y*y*y + startcy;
  }

But as I was going to wolframalpha to expand (a+bi)^4, I remembered DeMoivre’s Theorem. Duh. Here’s the general code:

    float r = sqrt(x*x + y*y);
    float theta = 0;
    r = pow(r,power);
    if (x > 0)
    {
      theta = atan(y/x);
    }
    else if (x<0)
    {
      theta = atan(y/x);
      theta = theta + PI;
    }
    else
    {
      theta = 0;
    }
    theta = power * theta;
    fx = r*cos(theta) + startcx;
    fy = r*sin(theta) + startcy;

This will work for any integer n for z^n+c. Great.
However this method is noticeably slower than the specific solutions. Bummer.

Why Stick to Integer Powers?

I (now) know that DeMoivre's Theorem doesn't "work" for non-integer powers, but it sort of works.
Here's a video of the fractal walking from z^2+c to z^3+c by steps up in the power of 0.05. On the live version, you can try this yourself by clicking the "Q" key to go up by 0.05, and the "A" key to go down by 0.05. Notice the vertical line right down the middle? That's due to converting a rectangular point to polar and the program is trying to get pi/2 out of a theta of arctan(some_number / some_number_thats_really_really_small). See the trouble?

Float Fail

Eventually we zoom too far down because the numbers are too small to accurately represent with a 32-bit float. Here's a picture of what happens when it eventually fails. At the bottom of the picture, I've printed out the resolution and the window width.
Screen Shot 2014-04-05 at 12.33.00 PM
A float only can store about 7 (decimal) digits of precision, so when we're squaring the same number 134 times, the errors will eventually catch up to us.

Mandelbrot Fractal

I was inspired by a My Favorite presentation from a student on the Mandelbrot fractal. I played around for a half hour at school without a ton of success, but I sorted out the sticking points on the way home, and finished it up tonight. Amazingly simple (the code may not look simple, but it’s all mine and most of it is setup. The actual fractal code is ~15 lines).

Screen Shot 2014-03-27 at 8.41.55 PM

Link to live openprocessing demo with source code.

Here’s the code, stripped of the color fluff.

float fx = 0;
float fy = 0;
int count;

int resolution = 18;

void setup() { 
    background(0);
    size(500, 500);
    loadPixels();
    color c;
    for (int i = 0; i < pixels.length; i++)
    {
        float startcx = 3.0 * (float(i % width) / width) - 2; //This provides a virtual x window from -2 to 1.
        float startcy = 3.0 * (float(floor(float(i)/width))/height) - 1.5; //This provides a virtual y window of -1.5 to 1.5.
        float mandelvalue = mandel(startcx,startcy);
        if  (mandelvalue < 2)  c = color(0);
        else  c = color(255);    
        pixels[i] = c;
    }
    updatePixels();
}

void draw() 
{
}

float mandel(float startcx, float startcy)
{
    float xm = 0;
    float ym = 0;
    count = 0;
    fx = 0;
    fy = 0;
    while ((count < resolution) && (fx*fx+fy*fy < 4))
    {
        func(xm,ym,startcx,startcy);
        xm = fx;
        ym = fy;
        count++;
    }
    return sqrt(fx*fx + fy*fy);
}

void func(float x, float y, float startcx, float startcy)
{
    // this squares the complex number x+yi, and then adds the constant
    fx = x*x - y*y + startcx;
    fy = 2*x*y + startcy;
}

Chaos Game v2

I prepared a small demonstration of the Chaos game for some math teachers. First we used transparencies, markers, dice, and rulers, but humans are mistake-prone and slow.
20140319-101745.jpg

To The Computer!

I decided to experiment with the rules of the game to see where it’d go.

Human Error

What happens when you bring in a random amount of error based on the distance traveled? Can you predict the result? The further the distance, the more error in the calculation.

Random Square

Random Hexagon

Weighted Die

What happens when you use a weighted die? Can you tell which vertex is being chosen more often than the others?

Weighted Die Triangle

Weighted Die Hexagon

Generalized Chaos Game

Live Demo.

Screen Shot 2014-03-18 at 8.46.47 PM

By this time, I was getting frustrated at using different processing files just because I wanted to change the number of vertices. I had originally “cheated” to find the coordinates of each point by using a website to calculate the rectangular coordinates of a hexagon. I knew there must be an easier way, but last summer I was more focused on the creation of the images. I figured it out this time.

Here’s the code:

//create points 
  float h1 = random(0.95);
  for (float theta = 0; theta < TWO_PI; theta += TWO_PI / numOfPoints)
  {
    //graph the points using polar form (rotate theta), then convert back to rectangular
    int xposition = int(radius * cos(theta))+width/2; //Add half width to move to 1st quad
    int yposition = int(radius * sin(theta))+width/2;

    h1 = h1 + 0.618033988749895; //Golden ratio!!
    h1 = h1 % 1;

    color c = color(h1,0.8,0.95);
    sp.add(new Spoint(xposition,yposition,c));
  }

For those who are code-adverse, instead of trying to calculate the rectangular coordinates of a regular n-gon directly, I calculate the polar coordinates and then convert back to rectangular. Duh! Notice the theta going from 0 to TWO_PI, and going up by (TWO_PI / numOfPoints), and then going back to rectangular by using x = r*cos (theta) and y = r*sin(theta).

This code snippet also highlights a neat trick to make the random colors that were spread out nicely (idea from here). The problem is that using random to get a decimal from 0 to 1 doesn’t tend to spread the color values out, there are lots of repeat colors. But if you chose the first color randomly, and then rotate by the golden ratio (add the conjugate:  0.618033988749895), then you get colors that start off random, but are nicely spread out. Here’s an example of the problem. The top graph is twenty random numbers, and the bottom graph is one random number spread out to 20 points with help from the golden ratio.

Screen Shot 2014-03-18 at 8.42.22 PM

Code.

Some Interesting Results

Go to the generalized chaos game, and reduce the number of vertices to 4, and increase the distance by two clicks (60% of the way to the next random choice). You should get something like this:

2014-03-19_11h29_00

Do the white lines occur at a geometric sequence? Why??

Now go to the hexagon, and increase the distance two clicks again, you should get something like this:

2014-03-19_11h33_37

What shapes are the overlapping colors (and the primary colors)? Why?

Lastly, some fun patterns happen when you continue to decrease the distance. Eventually, each point will actually move away from the randomly chosen point. Here are two such versions.

2014-03-19_11h37_11
2014-03-19_11h36_51

Happy Accidents

While trying to make the generalized chaos program, I ran into the following by accident:

Unique Quad v1:

2014-03-19_11h24_38

Unique Quad v2

2014-03-19_11h26_46

Now It’s Your Turn

Find me something cool. Let me know about it. Enjoy!

Buffon’s Needles simulation in Processing

Buffon’s Needle is a famous way to (slowly) estimate \pi.
Here’s a processing.org program to calculate \tau (to keep the math-hipster hatred of \pi-day at a critical point and concave up).

Link to live simulation and code. All variables are easy to change, size of window, length of needle, spacing of lines, etc.

2014-03-14_12h23_26

Neat program to write. Funnily enough, the code wasn’t working in the javascript version of processing because javascript didn’t understand the constant TAU. Ha.

note: I got the graph idea at the bottom idea from this sketch, but I never looked at the actual code.

Power Series on Desmos

This is a crosspost from my Photo 180 blog.

Power Series work in AP Calculus BC.

\sum_{n=0}^{\infty}(-1)^{n}\frac{x^{2n+1}}{(2n+1)!}=\frac{x}{1!}-\frac{x^3}{3!}+\frac{x^5}{5!}-...

Process: Since it’s a infinite series, look at partial sums to get an idea what this graph looks like.
So look at
y_1=\frac{x}{1!}
y_2=\frac{x}{1!}-\frac{x^3}{3!}
y_2=\frac{x}{1!}-\frac{x^3}{3!}+\frac{x^5}{5!}

Perfect time to use technology.

Texas Instruments Method

Go to y1. Enter in y=x.
Graph.
Wait 3-5 seconds.
Go to y1. Subtract a \frac{x^3}{3!}.
Graph.
Wait 3-5 seconds.
Go to y1. Add on a \frac{x^5}{5!}.
Graph.
Wait 3-5 seconds.
Go to y1. Subtract a \frac{x^7}{7!}.

2014-03-05_13h53_31 2014-03-05_13h53_44

REALLY crappy resolution. Awful zoom system. Where’s the factorial sign? Hopefully you remember what the previous graphs looked like. Lots of waiting. Ugh.

Desmos: version 1

Graph y=x
(NO WAIT STEP)
Subtract a \frac{x^3}{3!}.
(STILL NO WAIT STEP)
Add on a \frac{x^5}{5!}.
Subtract a \frac{x^7}{7!}.
2014-03-05_13h55_46
Great!

Desmos: version 2

Students teaching teachers: Have one your students find this out for himself, and remark that they can enter in the entire series, but he’s having trouble finding the infinity sign. SLIDERS!
Link.
Whoa. (That seemed like it should be much harder to type in. Took me 5 minutes to type in the original latex code at the beginning of this post. Took me < 1 minute to actually enter the series into desmos. I had to help 1 student out of 15. That’s it. Wow easy.)
2014-03-05_13h58_56

Watch the power series create sin(x) step by step by moving a slider. LIVE!

We live in good times.

 

Two James Tanton Questions

It’s midterm week at school, and James Tanton threw out two interesting questions in two days. I spent a little time programming “solutions” to these problems (not solutions, just verifications for an infinitesimally small portion of the natural numbers).

Problem One:

Here’s my processing.org code for this problem. And here’s the output of the code, each time the sequence gets longer, it prints out the new “max” sequence length. 2014-01-28_08h29_27

I didn’t use any of processing.org’s graphics but I had the prime function optimized, so it was quick work.

Problem Two:

Here’s the python code for the “solution”. And here’s the last six lines of the output of the code. 2014-01-28_08h31_32 I checked all numbers under 1,000,000, and all the sequences were finite (they stopped at a multiple of 13). The starting number whose sequence ended in the largest multiple of 13 was 964,665, and the multiple of 13 had 384 digits (BIG NUMBER! There are only ~10^80 particles in the entire universe). Fun stuff.

[Edit: 1/28/2014 9:08am] Ok, ok. James asked for a proof for the second problem. Here you go :-) 2014-01-28 09-05-38

Prime Matrix – Processing

Saw this tweet yesterday:

Pretty cool. Josh Giesbrecht also did some great work with Mathematica to replicate the image.

I also replicated the image with processing.org. The (live) code can be found here.

screen-0000

(update 1-16-2014 11:38, image above is now lossless and will stand up to zooming.)

Four and Five – Cartalk Puzzler

From cartalk.

Ray: Get a piece of paper and write the number four, leave a little space, and write the number five. What common mathematical symbol, when placed between the numbers four and five, will result in a number that is greater than four but less than six?

Tom: It has to be a mathematical symbol? It can’t be, like, the word “or”?

Ray: No, it’s got to be something that’s commonly used in mathematics. You’ve used it many, many times – maybe even today.

Pebbling a Checkerboard Game (Or Chessboard)

I happened upon this tweet when I got to work this morning:

So I watched the fantastic numberphile video Pebbling a Chessboard. I wanted to get the kids to play the game and I didn’t have any checkers in the classroom, so I decided to program the game in Processing. Thankfully the programming went smoothly, and I finished it in time to have the Pre-Calculus class try to beat the game (it’s the day before winter break, and we just finished sequences and series). Fantastic timing.

Rules:

  • Goal: remove all checkers from the green prison.
  • If you click on a checker, and there is space to the right and above, then that checker will disappear, and two clones will appear to the right and above.
  • Theoretically this game board extends to infinity to the right and above.

The game is hosted by openprocessing.org or hosted by recursiveprocess.com. It should work in any web browser, including smartphones. If you’d like to increase the checkerboard size, just go to the openprocessing.org site and “tweak” the code. The size of the board is set on the first two lines.

Enjoy!