Pictures of Julia and Mandelbrot Sets/Computer programs

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

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

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

float 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

}