Fractals with Octave: More on the Classic Sets

The story so far

• Classic Mandelbrot and Julia

In the second instalment of this series on fractals with Octave, I will generate more images by changing some of the parameters of the different functions presented in the first article. But before doing that, I'll briefly explain how to save the generated images and say a word on colour maps.

Saving images

The images that the functions described in the previous post generate are colour mapped images. Each element of a result matrix is an integer value that corresponds to a colour index. Octave has a very easy way to save that type of images using the well named saveimage function. To save the images generated, we just need to do:

octave:1> saveimage(file name, matrix, type, colour map);

The file name parameter is the name of the image file we want to save the image to. Give it a .ppm extension, we'll see why in a second. The matrix parameter is the matrix generated by the mandelbrot or julia functions presented in the previous article. The type parameter is the type of the image: we will use the PPM type as it is supported natively by Octave so that attribute should just be ppm. The colour map attribute is a colour map to use in the file. Octave has a default colour map that contains 64 colours from dark blue to light blue, then from yellow to red and white, similar to a contour map. It can be invoked by the colormap command. Putting all this together, creating a classic Mandelbrot image and saving it can be done this way:

octave:1> M=mandelbrot(-2.1+1.05i,0.7-1.05i,320,64);
octave:2> saveimage("mandelbrot.ppm", M, "ppm", colormap);

The PPM file format is understood natively by Linux but if you want to transform it into a more widely used format such as GIF or PNG, you can use ImageMagick or The GIMP. Avoid the JPEG format as it's not well suited to that type of computer generated images.

That's all well and good but what if you want your generated images to use a different colour map from the default? Let's create one. Octave understands a colour as a row vector that has 3 elements corresponding to the red, green and blue components of a colour. Each element is a real value between 0.0 and 1.0. So for example bright red would be [1, 0, 0] and 50% grey would be [0.5, 0.5, 0.5]. A colour map is a 3 column matrix where each row in the column is a colour. As an example, the following is a 5 colour map that defines a gradient from bright red to bright green:

   1.00000   0.00000   0.00000
   0.75000   0.25000   0.00000
   0.50000   0.50000   0.00000
   0.25000   0.75000   0.00000
   0.00000   1.00000   0.00000

We can build such a gradient by using the linspace command to evenly spread values. We will need to do that for each red, green and blue component and then combine all three with the cat command. Each vector produced by linspace will represent a column in the colour map matrix. However linspace produces row vectors so we need to transpose them before concatenating them using the ' operator. The following will produce a 64 value gradient from dark blue to white:

mycolmap=cat(2,
linspace(0, 1, 64)',
linspace(0, 1, 64)',
linspace(0.2, 1, 64)'
);

The resulting colour map can then be used to save the image. Note that for best results, the number of colours in the colour map should be the same as the number of iterations used when producing the fractal image.

Zooming into the Mandelbrot set

With the preamble out of the way, let's create some more fractal images. You can now save them for posterity! Some great images can be obtained by zooming somewhere on the edges of the Mandelbrot set. Let's zoom to the top centre of the set:

octave:1> Mcz=mandelbrot(-0.27+i,0.05+0.76i,320,64);
octave:1> imagesc(Mcz)

Octave should open the Gnuplot window and display the following image:

Mandelbrot Set, Top Centre Zoom

This image shows a number of artefacts that are similar to the complete set. Zooming further on those would show them in more details. And as this is a fractal set, you can zoom as much as you want and discover more interesting structures.

Mandelbrot with a non-null z₀

In the first article, we built the mandelbrot function based on the assumption that the initial value of the z_n series, z₀, was 0 for all values of c. It doesn't have to be. So let's create a modified mandelbrot function that can set a different value for z₀ and call it mandelbrotz0. Copy the mandebrot.m file described in the previous article into a mandelbrotz0.m and modify it accordingly:

mandelbrotz0.m

function M=mandelbrotz0(cmin,cmax,hpx,niter,z0)
  vpx=round(hpx*abs(imag(cmax-cmin)/real(cmax-cmin)));
  z=zeros(vpx,hpx).+z0;
  [cRe,cIm]=meshgrid(linspace(real(cmin),real(cmax),hpx),
                     linspace(imag(cmin),imag(cmax),vpx));
  c=cRe+i*cIm;
  M=mjcore(z,c,niter);
endfunction

Then, let's run this function with z₀=0.2+0.2i:

octave:1> Mz0=mandelbrotz0(-2.1+1.05i,0.7-1.05i,320,64,0.2+0.2i);
octave:2> imagesc(Mz0)

Octave should refresh the Gnuplot window and display the following image:

Mandelbrot Set, z₀=0.2+0.2i

Quirky! It's as if the Mandelbrot set was exploding like a ruptured balloon! Other values produce interesting results and the further they are from 0, the more distorted the set is.

Zooming into the Julia set

After the experiments with the Mandelbrot set above, let's explore the Julia set a bit more than we did last time, using c=-0.1+0.9i. Let's calculate the whole set for this value:

octave:1> Jc2=julia(-1.6+1.2i,1.6-1.2i,320,64,-0.1+0.9i);
octave:2> imagesc(Jc2)

Which produces the following image:

Julia Set, c=-0.1+0.9i

Then, let's zoom 10 times to a point slightly off centre:

octave:1> Jc2z=julia(-0.16+0.24i,0.16,320,64,-0.1+0.9i);
octave:2> imagesc(Jc2z)

Julia Set, c=-0.1+0.9i, zoom x10

And again:

octave:1> Jc2zz=julia(-0.016+0.072i,0.016+0.048i,320,64,-0.1+0.9i);
octave:2> imagesc(Jc2zz)

Julia Set, c=-0.1+0.9i, zoom x100

Cool pictures, don't you think?

Bootnote

Having created a slightly more generic mandelbrotz0 function above, we can now re-write and simplify the original mandelbrot function by calling mandelbrotz0:

mandelbrot.m

function M=mandelbrot(cmin,cmax,hpx,niter)
  M=mandelbrotz0(cmin,cmax,hpx,niter,0);
endfunction

Next

In the next article of this series, Other Polynomials, we look at what happens when we use a different complex series than the standard one.