Fractals with Octave: Other Polynomials

The story so far

• Classic Mandelbrot and Julia
• More on the Classic Sets

In the third instalment of this series on fractals with Octave, I want to change the polynomial series used to generate the Mandelbrot and Julia sets. So far, I've used the classic series, defined by:

z_n+1=z_n²+c

In the previous articles, I showed what happened when you varied the initial values for z₀ and c but why not go further and use a different polynomial altogether?

Code review

Before jumping into the thick of it, let's have a look at the code produced so far. We created 4 functions:

function M=mjcore(z,c,niter)

The core function shared by the mandelbrot and julia functions;

function M=mandelbrot(cmin,cmax,hpx,niter)

The function that produces a colour map image as per the Mandelbrot algorithm;

function M=mandelbrotz0(cmin,cmax,hpx,niter,z0)

A variation of the mandelbrot function that takes an extra argument in order to specify an initial z₀ value;

function M=julia(zmin,zmax,hpx,niter,c)

The function that produces a colour map image as per the Julia algorithm.

There is very little difference between the mandelbrot and mandelbrotz0 functions: mandelbrot just uses a default value of 0 for z0. In fact, we could merge both functions into one using the default argument construct supported by Octave. Here's a new version of the mandelbrot.m file that uses it:

mandelbrot.m

function M=mandelbrot(cmin,cmax,hpx,niter,z0=0)
  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

Having done this, we can now create the classic Mandelbrot image as we did previously, using four arguments:

octave-3.0.1:1> Mc=mandelbrot(-2.1+1.05i,0.7-1.05i,640,64);

Create the Mandelbrot image with z₀=0.2+0.2i by adding a fifth argument:

octave-3.0.1:2> Mz0=mandelbrot(-2.1+1.05i,0.7-1.05i,640,64,0.2+0.2i);

And we can delete the mandelbrotz0.m file. Note that you could also explicitly tell Octave to use the default for the fifth argument by specifying a colon (:) in its place, so you could produce the classic Mandelbrot this way:

octave-3.0.1:3> Mc=mandelbrot(-2.1+1.05i,0.7-1.05i,640,64,:);

This is not very useful in this context but we will see a use for it later.

Function handles and anonymous functions

In order to be able to change the polynomial used to calculate the fractal image, we need to be able to specify to the mjcore function what function to use when calculating the series. We could create a whole library of mjcore derivatives but that would quickly get tedious. Instead, we can use Octave's function handles. A function handle is a fancy way to say that in Octave, any variable can be a function. Therefore, we can just pass the function as an argument of the mjcore function. And while we're at it, we'll also pass r as an argument because, as we saw in the first article, the value for r is dependant on the function used in the series. So here is a modified version of that function:

mjcore.m

function M=mjcore(z,c,niter,f,r)
  M=zeros(length(z(:,1)),length(z(1,:)));
  for s=1:niter
    mask=abs(z)<r;
    M(mask)=M(mask)+1;
    z(mask)=f(z(mask),c(mask));
  endfor
  M(mask)=0;
endfunction

That's all good and well but we've just moved the problem further into the mandelbrot and julia functions: those functions will have to specify parameters accordingly and use sensible defaults so that if we omit the values for f and r, it will calculate the classic sets using the series z_n+1=z_n²+c (i.e. the function f=z²+c) and r=2. To do that, we are going to use the default argument construct explained above and combine it with an anonymous function. Using anonymous functions, we can define a variable f that is a handle on a function that implements the classic algorithm:

octave-3.0.1:4> f=@(z,c) z.^2.+c;

As shown above, we can use exactly the same syntax to specify default arguments to the mandelbrot and julia functions. Here are the new versions of those two functions:

mandelbrot.m

function M=mandelbrot(cmin,cmax,hpx,niter,
                      f=@(z,c) z.^2.+c,r=2,
                      z0=0)
  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,f,r);
endfunction

julia.m

function M=julia(zmin,zmax,hpx,niter,c,
                 f=@(z,c) z.^2.+c,r=2
                 )
  vpx=round(hpx*abs(imag(zmax-zmin)/real(zmax-zmin)));
  [zRe,zIm]=meshgrid(linspace(real(zmin),real(zmax),hpx),
                     linspace(imag(zmin),imag(zmax),vpx));
  z=zRe+i*zIm;
  cc=zeros(vpx,hpx)+c;
  M=mjcore(z,cc,niter,f,r);
endfunction

Cubic sets

We can now produce cubic Mandelbrot and Julia sets by setting the power to 3 instead of 2:

octave-3.0.1:5> Mn3=mandelbrot(-1.2+1.6i,1.2-1.6i,240,64,
> @(z,c) z.^3.+c);
octave-3.0.1:6> imagesc(Mn3)

Cubic Mandelbrot Set

octave-3.0.1:7> Jn3=julia(-1.6+1.2i,1.6-1.2i,320,64,0.4,
> @(z,c) z.^3.+c);
octave-3.0.1:8> imagesc(Jn3)

Cubic Julia Set for c=0.4

Quartic sets

Increasing the power to 4, we get quartic sets:

octave-3.0.1:9> Mn4=mandelbrot(-1.2+1.6i,1.2-1.6i,320,64,
> @(z,c) z.^4.+c);
octave-3.0.1:10> imagesc(Mn4)

Quartic Mandelbrot Set

octave-3.0.1:11> Jn4=julia(-1.6+1.2i,1.6-1.2i,320,64,-1,
> @(z,c) z.^4.+c);
octave-3.0.1:12> imagesc(Jn4)

Quartic Julia Set for c=-1

Note that for a series defined by the equation z_n+1=z_n^p+c, the resulting Mandelbrot set will show p-1 secondary "blobs" while the Julia set will show p secondary "blobs".

Decimal power

The power doesn't have to be an integer. Using a decimal value, such as 2.5, generates the following Mandelbrot image, that is halfway between the normal and cubic sets:

octave-3.0.1:13> Mn2_5=mandelbrot(-1.6+1.2i,1.6-1.2i,320,64,
> @(z,c) z.^2.5.+c);
octave-3.0.1:14> imagesc(Mn2_5)

Mandelbrot Set for p=2.5

Complex power

Using a complex power, such as 2+0.1i, generates a skewed set:

octave-3.0.1:15> Mn2_01i=mandelbrot(-2.2+1.2i,1.0-1.2i,320,64,
> @(z,c) z.^(2+.1i).+c);
octave-3.0.1:16> imagesc(Mn2_01i)

Mandelbrot Set for p=2+0.1i

More complex polynomial

Using a series based on a more complex polynomial, such as z_n+1=z_n²+z_n+c, generates the following set:

octave-3.0.1:17> Mp=mandelbrot(-2+1.05i,0.8-1.05i,320,64,
> @(z,c) z.^2.+z.+c);
octave-3.0.1:18> imagesc(Mp)

Mandelbrot Set for series z_n+1=z_n²+z_n+c

Bootnote

Having modified the mandelbrot function to include the f and r arguments before the z0 argument, we now need to use the colon argument for f and r to produce a classic Mandelbrot with a different value of z0:

octave-3.0.1:2> Mz0=mandelbrot(-2.1+1.05i,0.7-1.05i,640,64,:,:,0.2+0.2i);

I could have left the z0 argument in fifth position to avoid having to do this but I felt that the f and r arguments were more likely to be provided so are better specified first in the function's signature.

Next

In the next article in this series, Trigonometric and Exponential Functions, we'll look at non-polynomial series.