Subject: A Matlab program to generate the Mandelbrot Set

From: Shlomi Fish <[email-masked]>
Date: Fri, 18 May 2001 06:27:32 +0300 (IDT)

The Mandelbrot Set is the set of points in the complex plane for which the
iterative process z <- z^2 + c where c is the coordinate in the plane and
z is initialized to zero does not propagate to infinity. While it is easy
to determine if a point has overflowed after the N'th step (if its
absolute value is greater than 2), one cannot be completetly sure if it
will not overflow after an infinite number of steps. However, a limited
number of them (even 10 or 20) can usually determine it quite accurately.

It is a common practice to colour an image of the Mandelbrot Set with
different colours based on the number of iterations it took them to
propagate to infinity.

The day before yesterday I decided to try and write a Matlab program to
generate the Mandelbrot set. As Electrical Engineers we use Matlab a lot
, because it is very suited to what we need as engineers. However, Matlab
programs are generally interpreted, which makes their loops very slow.
Thus, we constantly need to think how to avoid loops as much as possible.

What Matlab does well is bulk operations on matrixes (or
multi-dimensional tensors in gernarl), because such routines are
hard-coded in C. So for example, we can add two matrixes, multiply them
(either matrix product or an element-by-element product), generete [
e^A(i,j) ], or [ A(i,j)^B(i,j) ], take a sub-tensor, etc. So we constantly
have to put all of our data inside one matrix, and use the Matlab
functions to further manipulate it quickly.

Thus, I thought that I should process all the pixels of the part of the
complex plane that I want to show at once, while keeping them in one big
matrix. The Matlab operation that performs a Mandelbrot iteration for
every pixel in a matrix is "value = (value .* value) + Z", where value
contains the current value and Z is a matrix containing the complex
coordinates. ".*" is an element-by-element product and "+" is naturally
matrix addition (which is done element-by-element anyhow).

I could repeat that (in a loop) for the number of steps I want, and that
way I will get the Mandelbrot values. This worked well for a small number
of steps (I tried 10), but in 20 I got disruptions in the resultant image.

What I found out was that because of the large number of iterations, some
of values overflowed to infinity, which confused the test for their
absolute value. So, I looked for a way to keep the values at bay.

Matlab also allows the operation "abs(values) > 2" which gives a matrix of
0 and 1 based on the comparison. I decided to use it as a mask to perform
two things:

1. Zero the values whose absolute value is greater than 2 so they will not
overflow.

2. Keep a mask of all the values that have already overflowed (at least
once) so far. This mask will be updated with the current mask by using a
logical or operation. (Yes, you guessed it - "or" in Matlab operates on
entire matrixes!)

That way, I got much better results for an iteration limit of 20.

Naturally, in a compiled programming language it would have been just as
fast (and probably faster) to draw the Mandelbrot set using a separate
loop for each pixel that stops just when it overflows. But because of
Matlab's tensor-oriented nature, I cannot afford to do this.

I believe it illustrates how the philosophy of a certain language affects
how one should program with it, and how it gives a new perspective and new
insights on programming in general.

In any case I included the code of the script below. Don't be alarmed by
its size because most of it is comments and default variable handling.
Those that do not have access to Matlab, can use Scilab (a free
Matlab-like replacement) and its Matlab->Scilab converter to run it.

To display it type something like that at the console:
steps=20;
A = mandel4(200,200,steps);
image(A/steps*64);
colormap(gray(64));

Regards,

	Shlomi Fish

%%%%%%%% Snip here and save as mandel.m %%%%%%%%%%%%%%%
function ret = mandel(x,y,steps,init_value)

% mandel - a function to generate the Mandelbrot Set.
% Written by Shlomi Fish, 2001
%
% This file is under the public domain.

% Assign some default values to the parameters
% nargin is an internal Matlab variable that specifies how many
% parameters were passed to the function.
if nargin < 1
    x = 100;
end

if nargin < 2
    y = 100;
end

if nargin < 3
    steps = 10;
end

if nargin < 4
    init_value = 0;
end

% Generate the coordinates in the complex plane
[X,Y]=meshgrid(-2:(4/x):2,-2:(4/y):2);
% Combine them into a matrix of complex numbers
Z=X+j*Y;
% Retrieve the dimensions of Z
dims = size(Z);
% The length in the x direction
x_len = dims(1);
% The length in the y direction
y_len = dims(2);
% value is initialized to init_value in every point of the plane
value = ones(x_len,y_len) * init_value;
% In the beginning all points are considered as part of the Mandelbrot
% set. Thus, they are initialized to zero.
ret   = zeros(x_len,y_len);
% The mask which indicates which points have already overflowed, is set
% to zero, to indicate that none have so far.
mask  = zeros(x_len,y_len);
% Perform the check "steps" times
for step=1:steps
    % For every point with a mandel value of "v" and a coordinate of "z"
    % perform  v <- (v ^ 2) + z
    %
    % .* is an element-by-element multiplication of two matrixes of the same size.
    % (regular * indicates matrix multiplication)
    value = (value .* value) + Z;
    % Retrieve the points that overflowed in this iteration
    % An overflowed point has a mandel value with an absolute value greater
    % than 2.
    current_mask = (abs(value) > 2);
    % Update the mask. We use "or" in order to avoid a situation where
    % 1 and 1 become two or so.
    mask = or(mask, current_mask);
    % Upgrade the points in the mask to a greater value in the returned
    % Mandelbrot-map. 
    ret = ret + mask;
    % Zero the points that have overflowed, so they will not propagate
    % to infinity.
    value = value .* not(current_mask);
end

% Now ret is ready for prime time so we return it.

%%%% --- SNIP HERE --- %%%%

----------------------------------------------------------------------
Shlomi Fish        [email-masked] 
Home Page:         http://t2.technion.ac.il/~shlomif/
Home E-mail:       [email-masked]

A more experienced programmer does not make less bugs. He just realizes
what went wrong more quickly.