## Friday, 27 March 2009

### Fractals with Octave: Other Polynomials

#### The story so far

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:

zn+1=zn2+c

In the previous articles, I showed what happened when you varied the initial values for z0 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 z0 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 z0=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
endfor
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 zn+1=zn2+c (i.e. the function f=z2+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)
```
```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)
```

#### 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)
```
```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)
```

Note that for a series defined by the equation zn+1=znp+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)
```

#### 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)
```

#### More complex polynomial

Using a series based on a more complex polynomial, such as zn+1=zn2+zn+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 zn+1=zn2+zn+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.