The Mandelbrot Set in R

I was having a conversation with someone I know about weather forecasts the other day and it went something like this:

“Yes, their forecasts are really not very good. They really need to work on improving them.”

“Well, yes, I think they’re okay, considering what they’re trying to do is impossible.”

The thought that a relatively simple set of equations can give rise to infinitely complex behavior is something which has intrigue and appeal beyond those mathematically minded and academic. Fractals became a part of pop culture, and the idea of chaos is not unknown to those outside mathematical research, as it was popularized by the term “The Butterfly Effect” (though unfortunately I can’t say the movie starring Ashton Kutcher was anything to write home about).

When I was finishing high school and starting my undergraduate degree I got very interested in fractals and chaos. What really spurred this interest was James Gleick’s Chaos, which outlined the history of the “mathematics of chaos” in a highly readable, narrative way. So much so that I later went on to take a pure mathematics course in chaos during my degree, and wrote a program in MATLAB for exploring the Mandelbrot set.

So I thought I’d give it a shot in R.

Unquestionably the most well known fractal of all is the Mandelbrot fractal. Without going into laborious mathematical detail, the Mandelbrot set stems from using complex numbers, which are numbers of the form z = x + yi. The number i, the “imaginary unit”, is the special quantity such that i² = –1. As such complex numbers are unique in that multiplying two of them together can result in much different behavior than numbers on the real number line: the magnitude of their product is not guaranteed to be greater than the terms multiplied, even for two quantities with real and imaginary parts with magnitude greater than one.

The Mandelbrot set the set of numbers produced by the following:

1. Initialize a complex set of numbers in the complex plane that are all z = 0.

2. Iterate the formula

Execution

In R this is straightforward code. Pick some parameters about the space (range and resolution for x and y, number of iterations, etc.) and then iterate. Naively, the code could be executed as below:

However, if you know anything about writing good code in R, you know that for loops are bad, and it’s better to rely upon the highly optimized underpinnings of the language by using array-based functions like which and lapply.  So a more efficient version of the same code is below:

We can then plunk these into functions where the different criteria for the rendered fractal are parameters:

Let’s compare the runtime between the two shall we? For the base settings, how does the naive version compare to the vectorized one? 

The results speak for themselves: a ~65x decrease in runtime by using R indexing functions instead of for loops! This definitely speaks to the importance of writing optimized code taking advantage of R’s vectorized functions.

You can tweak the different parameters for the function to return different parts of the complex space. Below are some pretty example plots of what the script can do with different parameters:

What does all this have to do with my conversation about the weather? Why, everything, of course! It’s where chaos theory came from. Enjoy the pretty pictures. It was fun getting R to churn out some nice fractals and good to take a trip down memory lane.

The Mandelbot Set: 

code on github: