Generating Julia Sets

Retired for over 2 years, I've devoted most of my time to photography, fell walking, and renovating my long neglected home. It's now time to pay some attention to my other passions: mathematics and computing.

I decided to indulge myself with the visual side of mathematics by generating some fractal images. This was something I'd tried to do back in the 1980's, but the speed of PC's in those days made the exercise prohibitive.

Although these images are complex, the maths and programming isn't. ( I can only assume that those people who choose to use ready written program like Apophysis, either don't understand the maths or can't program - I'm for doing it all myself).

Take for instance the Julia sets (one example is given below).
Julia Set for z2+ 0.32 + 0.043i

The Julia sets live in the complex plane and are crucial to the understanding of iterations of polynomials like z2 + c.

For the construction of Julia sets we fix c, and choose some initial value for z = z0, say. Then we repeat the process and calculate z1 = z02 + c, z2 = z12 + c, and so on.

In other words, for an arbitrary but fixed value of c, we generate a sequence of complex numbers {z0, z1, z2,......, zn, zn + 1, ......} for each pixel. That's 16002 = 2,560,000 pixels in the images shown here.

Each sequence has one of two properties, either:
  1. the sequence remains bounded, i.e. you can draw a circle of (predefined), radius M around z0, and all members of the sequence will lie within its boundary, or
  2. the sequence becomes unbounded, i.e. after a number of iterations the members of the sequence move outside the boundary of the circle.
The collection of points which lead to the first kind of behaviour is called the prisoner set for c, while the sequence of points which lead to the second kind of behaviour is called the escape set for c.

By writing a simple program we can investigate the behaviour of a whole range of starting values for z, and determine which starting values belong to the prisoner and escape sets.

In order to illustrate graphically if a point is in an escape or prisoner set we need to assign a colour to each point. Now if we define K to be the number of colours we can use, we can assign a colour number, that lies in the range 0 to K-1, to each starting point. That is, for each starting value, if the kth iteration (k < K), escapes the circle then assign the starting value the colour k. But if after K iterations all members of the sequence are still within the circle then assign the starting value the colour 0. Thus all prisoner values will have colour 0 and all escape values will have a colour in the range 1 to K-1.

In the fractals shown here, all pixels coloured Red represent prisoner values and all the other pixels (including those coloured a lighter shade of red, represent escape values.

zoomed in at 3 times magnification
So how do we plot the values? Well if you remember we said that z is a complex number, i.e. z is of the form z = x + iy, where the numbers x and y are real numbers and i is the special number equal to the square root of -1, (so i2 = 1).
We call x the Real part of z and y the imaginary part of z.
i is called an imaginary number. Like infinity it cannot be calculated it is a mathematical definition, a convenience that enables mathematicians to investigate complex mathematical problems - like fractals.

This means we can represent any value z as a pixel in the (x,y) plane with a colour k. Here is the Algorithm I used:

Decompose the complex numbers z and c into their real and imaginary parts,
i.e. z = x + iy and c = p + iq, then (using the definition i2 = 1)
z2 + c = (x+iy)2 + p + iq
= x2 + 2ixy + i2y2 +p + iq
(but i2 = -1, so substituting -1 for i2 and rearranging we have)
=x2 - y2 + p + i(2xy + q).
In other words the real part of z2 + c is x2 - y2 + p, and the imaginary part (the bit multiplied by i )is 2xy + q.
So in order to evaluate the sequence of complex values zk = zk-12 + c we evaluate:
xk = xk-12 - yk-12 + p, 
yk = 2xk-1yk-1 + q 
........... Eqns 1

Now; assume the image we are going to produce has a graphical resolution of a times b pixels, (1600 x 1600 say). Let there be K+1 colours that can be displayed simultaneously, numbered 0 through K, (K = 512 say).
Step 0:
Choose a parameter c = p + iq, (0.32 + 0.043i say)
Choose xmin = ymin = -1.5, xmax = ymax = 1.5 say,
and define M (=100, say. i.e. our boundary circle will have a radius of 10).
Set dx = (xmax - xmin)/(a-1), and
dy = (ymax - ymin)/(b-1).
Then for all pixels (nx, ny) with nx = 0, ....., a-1 and
ny = 0, ......, b-1
we go through the following routine:
Step 1: Set 
x0 = xmin + nxdx
y0 = ymin + nydy
k = 0
Step 2: (iteration step):
Calculate (xk + 1, yk + 1) from (xk, yk) using Eqns 1.
Increase the counter k by 1 (k = k + 1).

Step 3: (evaluation step):
Calculate r = xk2 + yk2
  1. If r > M then choose colour k and go to Step 4.
  2. If k = K then choose colour 0 and go to Step 4.
  3. r ≤ M and k < K. Go back to Step 2.
Step 4:  Assign colour k to the pixel (nx, ny) and go to the next pixel, starting with Step 1.

Remarks: Varying the values of xmin, xmax, ... in Step 0 you may produce blow ups. For example: the first image was produced using values -1.5 to 1.5, while the second image was produced by substituting values -0.5 to 0.5.

You can experiment by programming this algorithm in any language, Basic, Java, C, etc. but I used a favourite package of many mathematical programmers called Mathematica. Below is the actual code I used.

The boundary of the prisoner set is simultaneously the boundary of the escape set, and that joint boundary is the Julia set for the given c (and M).

No comments:

Post a Comment


© All rights reserved

Search This Blog