Pictures of Julia and Mandelbrot Sets/Computer programs

From Wikibooks, open books for an open world
Jump to: navigation, search

We will show two ways of making a program that can draw a Mandelbrot set where the iterations are towards ∞: a program where all is done from the ground, but which does nothing more than drawing a single picture, and a program where you can do all you desire (zoom, alter colouring and render as file), but where these facilities are handed over to another program, namely Ultra Fractal.

A single picture[edit]

Section of the Mandelbrot set for z^2 + 0.009/z + c

Let us make a program that is as simple as possible: it draws a single picture, and does nothing more (when it has done its work it closes). It was originally written in assembly language, but here we write it in pseudo-code. It consists of two parts: the drawing procedure and the colour scale.

We draw a section of the Mandelbrot set for the family f(z) = z^2 + \frac{p}{z} + c \,, for p = 0.009. We let the size of the window be 800x600 pixels, and we imagine that the section has its lower left corner in the point with coordinates (ax, ay) and that is has width h. We have f'(z) = 2z - \frac{p}{z^2} \,, so that the finite critical point we need, is solution to the equation z^3 = p/2. We choose the real solution cri = (p/2)^{1/3}.

The drawing procedure

p = 0.009
cri = (p/2)^{1/3}
r = 1.0E200 (square of the bail-out radius)
u = log(log(r))
v = 1/log(2) (2 is the degree of the rational function)
ax = -0.7028233689263 (x-coordinate of lower left corner)
ay = 0.1142331238418 (y-coordinate of lower left corner)
h = 0.0092906805859 (width of the section)
g = h / 800
m = 850 (maximum iteration number)
thick = h * 0.0005 (thickness of the boundary)
dens = 24.9 (density of the colours)
disp = 432.0 (displacement of the colour scale)

for i = 0 to 799 do

begin
cx = ax + i * g (x-coordinate of pixel (i, j))
for j = 0 to 599 do
begin
cy = ay + j * g (y-coordinate of pixel (i, j))
x = cri
y = 0
xd = 0
yd = 0
f = 0
n = 0
while (n < m) and (f < r) do
begin
n = n + 1
x1 = x * x - y * y
y1 = 2 * x * y
f = x * x + y * y
x2 = 2 * x - p * x1 / (f * f)
y2 = 2 * y + p * y1 / (f * f)
temp = xd
xd = x2 * xd - y2 * yd + 1
yd = x2 * yd + y2 * temp
fd = xd * xd + yd * yd
x = x1 + p * x / f + cx
y = y1 - p * y / f + cy
f = x * x + y * y
end
if (n = m) or (log(f) * sqrt(f) < thick * sqrt(fd)) then
setpixel(i, 599 - j, 0)
else
begin
s = n - v * (log(log(f)) - u)
n = round(dens * s + disp) mod 720
col = paletteRGB(col1[n], col2[n], col3[n])
setpixel(i, 599 - j, col)
end
end
end

Explanation

We have set c = cx + icy, z = x + iy, z' = xd + iyd and x1 + iy1 = (x + iy)^2, and we have used that 1/z = z^-/|z|^2.

The successive calculation of the derivative z' is xd + iyd = (x2 + iy2) * (xd + iyd) + 1, where x2 + iy2 = 2 * (x + iy) - p * (x1 - iy1) / (f * f) and f = x^2 + y^2. The next point in the iteration is x + iy = (x1 + iy1) + p * (x - iy) / f + (cx + icy). The distance function is

log(√f)*√f/√fd   (= log(f)*√f/(2*√fd))

for the last z-value, here f = x^2 + y^2 and fd = xd^2 + yd^2. The number in the interval [0, 1[ to be subtracted from the iteration number (in order to get the real iteration number), is log(log(|z|)/log(r))/log(2) = v * (log(log(f)) - u), for the last z, here v = 1/log(2) and u = log(log(r)).

paletteRGB(r, g, b) is the integer indexing the colour with RGB-values (r, g, b), 0 gives black. col, col2 and col3 are the arrays from 0 to 719 of integers from 0 to 255 constructed in the next section.

The colour scale

A colour is immediately (in the computer) given by its composition of the three primary colours red, green and blue, and their shares are measured in whole numbers between 0 and 255. This triple of numbers is the set of RGB values of the colour. The colours correspond accordingly to the entire points in a cube with side length 256, and we construct our cyclic colour scala by going regularly along an ellipse in this cube. An ellipse in the space is determined by:

a centre with coordinates (a, b, c)
a major axe rmaj and a minor axe rmin
two angles angle g and h determining the direction of the plane of the ellipse
an angle q determining the direction of the ellipse in this plane
u = cos(g * pi / 180) * cos(h * pi / 180) {(u,v,w) is a unit vector orthogonal to the plane}
v = cos(g * pi / 180) * sin(h * pi / 180)
w = sin(g * pi / 180)
x = -a {(x,y,z) is a vector in the plane}
y = -b
z = a * u / w + b * v / w
e = sqrt(sqr(x) + sqr(y) + sqr(z))
x1 = x / e {the unit vector corresponding to (x,y,z)}
y1 = y / e
z1 = z / e
x2 = v * z1 - w * y1 {the unit vector in the plane orthogonal to (x1,y1,z1)}
y2 = w * x1 - u * z1
z2 = u * y1 - v * x1
e = cos(q * pi / 180)
f = sin(q * pi / 180)
x1 = e * x1 + f * x2 {the unit vector in the plane in the direction of the angle q}
y1 = e * y1 + f * y2
z1 = e * z1 + f * z2
x2 = v * z1 - w * y1 {the unit vector in the plane orthogonal to this direction}
y2 = w * x1 - u * z1
z2 = u * y1 - v * x1
for i = 0 to 719 do
begin
e = rmaj * cos(i * pi / 360)
f = rmin * sin(i * pi / 360)
col1[i] = round(a + e * x1 + f * x2) {the three coordinates of the point on the ellipse}
col2[i] = round(b + e * y1 + f * y2)
col3[i] = round(c + e * z1 + f * z2)
end

In the picture we have: (a, b, c) = (176, 176, 176), rmaj = rmin = 88, g = 146, h = -32, q = 0.

Ultra Fractal[edit]

Ultra Fractal is a program for drawing fractals where the user to a great extent can write his own formula programs and where the main program performs all that is common for all the fractal procedures, that is, the procedure of going from pixel to pixel, zooming, colouring, production of a large picture as a file. Your effort is reduced to a minimum, and apart from pictures that cannot be drawn regularly from pixel to pixel - attractors and landscapes, for instance - you can do all you desire: non-complex functions, quaternions, ... You can operate directly with complex numbers.

First you should download a free trial-version of the program and learn how it works. Then you should see the article Make attractive pictures of Mandelbrot and Julia sets with Ultra Fractal, and download the (free) programs from this site to be run with Ultra Fractal.

For the drawing two programs have to be combined: a formula program and a colouring program. In the formula program is a section called loop, and when this loop has done its work, the number of iterations and the last value of the complex number z are transferred to the colouring program, which calculates the colour from these two numbers. However, Ultra fractal is designed to "all" types of fractals, and unfortunately the Julia and Mandelbrot sets do not quite fit into this form: for iterations towards cycles, the iteration has to be continued in order to find the order and the attraction of the cycle, and it is not the last value of z in the sequence we do use for the colour. Therefore we replace this loop by a fictive loop (which only runs one time) and perform all the operations in the section init. This implies that we cannot use the default colouring program of Ultra Fractal, we must write a new colouring program.

Another thing you must be aware of, is that in Ultra Fractal the norm |z| of the complex number z is the square of its length: |z| = x^2 + y^2.

Section of the Mandelbrot set for z^2/2 + 0.26*z^4 + c

We will show a program that draws the Mandelbrot set for the family f(z) + c, where f(z) is a polynomial whose first two terms are zero, so that it begins with z^2 or a larger power of z, because we then can take 0 as the finite critical point.

The two programs can be copied and inserted in an empty formula and colouring document, respectively. We have inserted the polynomial z^2/2 + p*z^4 + c of degree 4, where p is a real parameter. The picture shows a section of the Mandelbrot set for p = 0.26.

The formula program

Mandelbrot {
global:
float p = 0.26  ;parameter
float deg = 4  ;degree of the polynomial
float v = 1 / log(deg)
float g = 10 * log(10)
float r = exp(g)  ;square of the radius of the bail-out circle
u = log(g)
float tb = sqr(@thick / (1000 * #magn))  ;for the thickness of the boundary
float h = 1 / (1500 * #magn * @width)  ;a very small real number
init:
complex z = 0  ;critical point
complex zd = 0  ;the sequence of the derivatives
complex z1 = 0
float w = 0
int n = 0
while n < #maxit && |z| < r
n = n + 1
z1 = z^2/2 + p*z^4
zd = ((z+h)^2/2 + p*(z+h)^4 - z1) * zd / h + 1
z = z1 + #pixel
endwhile
if n == #maxit || sqr(log(|z|)) * |z| < tb * |zd|
w = -1
else
w = n - v * (log(log(|z|)) - u)
endif
 ;begin fictive loop
z = 0
n = 0
loop:
n = n + 1
z = z + #pixel
if n == 1
z = w
endif
bailout:
n < 1
 ;end fictive loop
default:
title = "Mandelbrot"
maxiter = 100
param thick
caption = "boundary"
default = 1.0
endparam
param width
caption = "width"
default = 640
endparam
}

Explanation

We have set the square of the radius r of the bail-out circle to 10^{10}, and set the thickness of the boundary to "thick/(1000 * #magn)", so that it becomes thinner when we zoom in.

We have calculated the derivative f(z) by (f(z+h) - f(z))/h for a small number h, namely "h = 1/(1500*#magn*width)", where "width" is the width of the picture. The default value is the width of the window in pixels (here set to 640), but if you draw a large picture, you should enter the width of this.

The successive calculation of the derivative z'_{k+1} = f'(z_k)z'_k + 1 is "zd = (f(z+h) - f(z)) * zd / h + 1". The next iteration z_{k+1} = f(z_k) + c is "z = f(z) + #pixel". The square of the distance estimation log|z_k| * |z_k| / |z'_k| is "sqr(log(|z|)) * |z| / |zd|". The number to be subtracted from the iteration number log(log(|z_k|)/log(r))/log(deg) is "v * (log(log(|z|)) - u)", where "v = 1/log(deg)" and "u = log(log(r))".

The final complex number z (to be transferred to the colouring program), which here actually is real, is set to -1, when the point belongs to the Mandelbrot set or to the boundary, otherwise it is set to the real iteration number. It is transferred to the colouring program, which we have called Gradient, as #z, and is here set equal to the real number s. When s is negative the colour is set to "#solid", otherwise we multiply s by a number "dens" determining the density, and add to this number a number "disp" determining the displacement of the scale, and as we want this displacement to be a per cent, we divide the result by 100: "u = (dens * s + disp) / 100". In Ultra Fractal the colours of the cyclic colour scales are indexed by the real numbers in the interval [0, 1[, and we let this number, "#index", be the non-integral part of the number u: "#index = u - trunc(u)".

The colouring program

Gradient {
final:
float s = real(#z)
float u = 0
if s < 0
#solid = true
else
u = (@dens * s + @disp) / 100
#index = u - trunc(u)
endif
default:
title = "Gradient"
param disp
caption = "displace"
default = 0
endparam
param dens
caption = "density"
default = 1.0
endparam
}