Julia Sets in Julia

Subtitle: Completely invariant iterations of the holomorphic function

Gaston Julia

Gaston Julia was born on the 3rd of February, 1893 in Sidi bel Abbès, Algeria to Delorès Delavent and Joseph Julia, a mechanic. Gaston became interested in music and mathematics at an early age, and excelled in both, settling on the violin as his favorite instrument.

In 1901, the family moved to Oran on the coast of Algeria where he enrolled in the Lycée, eventually graduating with distinction in the baccalaureate examinations of science, modern languages, philosophy, and mathematics. He left Algeria for Paris and received a scholarship to the École Normale Supériore in 1911. In 1914 Germany declared war on France and Julia was called up for service with the 57th Infantry Regiment at Libourne.

Deployed to the Western Front, the following is a report of the action, and Julia’s injury that cost him his nose.

January 25, 1915, showed complete contempt for danger. Under an extremely violent bombardment, he succeeded despite his youth (22 years) to give a real example to his men. Struck by a bullet in the middle of his face causing a terrible injury, he could no longer speak but wrote on a ticket that he would not be evacuated. He only went to the ambulance when the attack had been driven back. It was the first time this officer had come under fire.

For the remainder of his life, Julia wore a leather strap to cover the missing nose.

gaston-julia

Julia married Marianne Chausson in 1918. She was the nurse who attended to him after his injury, and the daughter of composer Ernest Chausson. He had completed his Ph.D. in mathematics the year before, under Émile Picard, Henri Lebesgue, and Pierre Humbert.

At age 25, he published a 199-page manuscript, Mémoire sur l’iteration des fonctions rationelles, (Memoir on the iteration of rational functions) which brought him instant fame in mathematical circles and earned him the Grand Prix des Sciences Mathématiques of the French Academy of Sciences.

An iterated function is applied repeatedly to the previous solution, and a rational function is one like

f(x)=P(x)Q(x)f(x) = \frac{P(x)}{Q(x)}

where P(x)P(x) and Q(x)Q(x) are polynomial functions. Q(x)Q(x) can’t be zero, but it can be 11, meaning that f(x)=P(x)f(x) = P(x).

Function Iterations

Iterating a function is pretty simple. Suppose

f(x)=x+1.f(x) = x + 1.

If you start with x=0x = 0, then f(0)=1f(0) = 1. Now put 11 into the function to get f(1)=2f(1) = 2. This goes on forever, 1,2,3,1,2,3, \ldots A more interesting function is

f(x)={x2ifmod(x,2)=03x+1ifmod(x,2)=1.f(x) = \begin{cases} \frac{x}{2} & \text{if} \mod(x,2) = 0 \\ 3x+1 & \text{if} \mod(x,2) = 1. \end{cases}

In other words, if xx is even, divide it by 22, but if it’s odd, multiply by 33 and add 11. Any number that’s a power of 22 immediately collapses down to 11, for example, x0=8f(8)=4f(4)=2f(2)=1x_0 = 8 \Rightarrow f(8) = 4 \Rightarrow f(4) = 2 \Rightarrow f(2) = 1. Then you get stuck in the loop 4,2,1,4,2,1,4,2,1,4,2,1, \ldots If you start with 33 you get the sequence {3,10,5,16,8,4,2,1}\{3,10,5,16,8,4,2,1\}.

Can you find a starting number that doesn’t wind up at 11? Nobody knows for sure, but it’s thought that every sequence eventually winds down to 11. This is known as the Collatz conjecture, named for mathematician Lothar Collatz who first proposed the conjecture in 1937. It’s also called the Hailstone Problem because of the way the number sequence jumps up and down.

Curses, FOILed Again!

Remember the FOIL method from high school algebra? If you want to expand (x+2)(x3)(x+2)(x-3) using the FOIL method, you multiply the First two terms, x×x=x2x \times x = x^2, then the Outer terms, x×3=3xx \times -3 = -3x, the Inner ones, 2×x=2x2 \times x = 2x, and finally the Last two, 2×3=62 \times -3 = -6 to get

(x+2)(x3)=x2x6.(x+2)(x-3) = x^2-x-6.

curses-foiled-again

Snidely Whiplash often said, “Curses, foiled again!” on the Rocky and Bullwinkle show, but he may not have been referring to algebra. Here’s a picture of Rocky, Bullwinkle, and Captain Peachfuzz. Because reasons.

rocky-bullwinkle-peachfuzz

What’s the square root of negative one (1)(\sqrt{-1})? For a long time, mathematicians thought there was no solution to this problem because when you multiply two positive numbers together you get a positive number, and when you multiply two negative numbers together you also get a positive number, so there can’t be an xx such that x×x=x2=1x \times x = x^2 = -1. But, in the 17th century, René Descartes named the solution to 1\sqrt{-1} as ii, an imaginary, fictitious, and useless number.

Since then, we’ve found plenty of uses for these otherwise useless numbers. Combine them with real numbers, the kind we’re used to, and you get complex numbers like 4+3i4 + 3i which is 44 parts real and 33 parts imaginary. You can plot complex numbers as points, using the xx-axis for the real part, and the yy-axis for the imaginary part. Here’s a plot of z1=4+3iz_1 = 4 + 3i:

complex-number

Suppose you wanted to multiply z1z_1 by another complex number z2=2+5iz_2 = 2 + 5i. Just use the FOIL method but with the complex numbers,

z1×z2=(4+3i)×(2+5i)=42+45i+3i2+3i5i=8+20i+6i+15i2=(815)+26i=7+26i.\begin{aligned} z_1 \times z_2 &= (4 + 3i) \times (2+5i) = 4 \cdot 2 + 4 \cdot 5i + 3i \cdot 2 + 3i \cdot 5i \\ &= 8 + 20i + 6i + 15i^2 \\ &= (8-15)+26i \\ &= -7 + 26i. \end{aligned}

Remember that i2=1i^2 = -1, so multiplying the two imaginary parts together gives the negative of the product, ignoring the ii’s. It might make more sense to do the first and last together to get the new real part, and then the inner and outer parts for the new imaginary. You could call it the FLIO method.

In the plot, the real part of z1z_1 is 44, which is drawn as the vector uu, and the imaginary part is 33, drawn as the vector vv. These two vectors form the adjacent and opposite sides of a right triangle with angle θ\theta. Using the Pythagorean Theorem, we can calculate the length of the hypotenuse, r=42+32=25=5r = \sqrt{4^2+3^2} = \sqrt{25} = 5. The sides uu and vv can be written in terms of rr and θ\theta as

u=rcosθv=rsinθ.\begin{aligned} u &= r \cos \theta \\ v &= r \sin \theta. \end{aligned}

A Most Beautiful Equation

Leonard Euler (see our previous post, “Seven Bridges for Seven Truckers”) discovered an amazing property of complex numbers,

eiθ=cosθ+isinθ,e^{i \theta} = \cos \theta + i \sin \theta,

where ee is the base of the natural logarithm. The proof is pretty simple. Define f(θ)f(\theta) as

f(θ)=cosθ+isinθeiθ=eiθ(cosθ+isinθ).f(\theta) = \frac{\cos \theta + i \sin \theta}{e^{i\theta}} = e^{-i \theta} (\cos \theta + i \sin \theta).

The reason for choosing this particular function is that we will show that f(θ)=1f(\theta) = 1 for every θ\theta, meaning that the numerator and denominator are the same everywhere. Using the product rule of differentiation,

f(θ)=ddθeiθ(cosθ+isinθ)=eiθ(sinθ+icosθ)ieiθ(cosθ+isinθ)=eiθ(sinθ+icosθ)eiθ(sinθ+icosθ)=0.\begin{aligned} f'(\theta) &= \frac{d}{d \theta} e^{-i \theta} (\cos \theta + i \sin \theta) \\ &= e^{-i \theta} (-\sin \theta + i \cos \theta) - ie^{-i \theta} (\cos \theta + i \sin \theta) \\ &= e^{-i \theta} (-\sin \theta + i \cos \theta) - e^{-i \theta} (-\sin \theta + i \cos \theta) = 0. \\ \end{aligned}

Since the derivative of f(θ)=0f(\theta) = 0 for every θ\theta, then f(θ)f(\theta) is a constant. The derivative of a function gives the slope of the function at every point, so a derivative of zero means the function has zero slope, in other words, it’s a straight line parallel to the xx-axis at some constant distance from the axis.

We can calculate the value of this constant at any θ\theta, and the easiest one is θ=0\theta = 0, where f(θ)=1f(\theta) = 1 since ei0=1e^{-i \cdot 0} = 1, cos0=1\cos 0 = 1, and sin0=0\sin 0 = 0. In other words, f(θ)=1f(\theta) = 1 for every θ\theta, which completes the proof.

The angle θ\theta is measured from the xx-axis, so θ=0\theta = 0 is along the positive real number line, θ=π2\theta = \frac{\pi}{2} is the positive imaginary axis, θ=π\theta = \pi represents negative real numbers and θ=3π2\theta = \frac{3\pi}{2} points along the negative imaginary axis.

e-i-pi

Looking at z2z_2, we see that when θ=π\theta = \pi, eiθ=eiπ=1e^{i \theta} = e^{i \pi} = -1, which gives the most beautiful equation,

eiπ+1=0.\LARGE e^{i \pi} + 1 = 0.

Physicist Richard Feynman said that this is “the most remarkable formula in mathematics” because it contains ee, the base of natural logarithms, the imaginary number ii, π\pi, the constant relating the radius of a circle to its circumference and area, 11 the first natural number, and 00 the center of the mathematical universe.

Iterations on a Complex Number

Julia wondered what would happen if he iterated a function like f(z)=z2+cf(z) = z^2 +c where both zz and cc could be complex. If z=0z = 0 and c=1ic = 1-i, then the first few iterations become

z=0z=1iz=13iz=77iz=1+97i\begin{aligned} z &= 0 \\ z &= 1-i \\ z &= 1-3i \\ z &= 7-7i \\ z &= 1+97i \\ \end{aligned}

If you start with a smaller cc, say c=0.79+0.15ic = -0.79 + 0.15i, the values of zz jump around, but they stay relatively small. Another way to look at the iterations is to plot several steps of the iterations at different starting points z0z_0 and additive constant cc.

Since the first iterate squares the initial value of z0z_0 and adds cc, then z0=1z_0 = 1 and z0=1z_0 = -1 will meet at c+1c+1 after the first step. The same thing happens for z0=±iz_0 = \pm i.

In this plot, z0=(0.1+eiθ)/2z_0 = (0.1+e^{i\theta})/2 for θ=30,120,210,300\theta = 30,120,210,300 degrees, and c=0.1+0.1ic = -0.1+0.1i.

several-iterates

Adding 0.10.1 shifts the starting points a little bit to prevent overlapping trajectories, and dividing by 22 is equivalent to setting r=12r=\frac{1}{2} in the uu and vv vectors except for the slight shift of the starting point. If you set z0=0.1+eiθz_0 = 0.1+e^{i\theta}, eventually all of the trajectories escape towards infinity.

several-iterates-escape

What is it that keeps some points from getting away, while others escape? The norm of a vector is the length of the vector and as we saw above, the length is easily calculated using the Pythagorean theorem. A shorthand notation for the norm is a pair of vertical lines around the vector or complex number,

norm(z)=z=u2+v2.\text{norm}(z) = |z| = \sqrt{u^2+v^2}.

After a single iteration, zz becomes z2+cz^2 + c, and the length or norm is z2+c|z^2+c|. One way to check if the iterates are tending towards infinity is to calculate the length of each iterate. In the case θ=30°=0.52 radians\theta = 30 \degree = 0.52 \text{ radians} (without dividing by 22), the norm of the iterates is

z=[1,1.046,1.235,1.428,2.153,4.508,20.177,407.001,165649.888,2.744e10].z = [1, 1.046, 1.235, 1.428, 2.153, 4.508, 20.177, 407.001, 165649.888, 2.744e10].

Using the same starting point, and now dividing by 22, the iterates are

z=[0.5,0.317,0.231,0.091,0.134,0.129,0.128,0.129,0.128,0.128].z = [0.5, 0.317, 0.231, 0.091, 0.134, 0.129, 0.128, 0.129, 0.128, 0.128].

Not only do the iterates remain small, but they get closer to 0.1280.128 (after rounding).

What happens if you write z=reiθz = re^{i \theta} and want to find z2z^2? This is very easy because

z2=reiθ×reiθ=r2×(eiθ)2=r2e2iθ.\begin{aligned} z^2 &= re^{i \theta} \times re^{i \theta} \\ &= r^2 \times (e^{i \theta})^2 \\ &= r^2 e^{2i \theta}. \end{aligned}

How long is the new vector z2z^2? The e2iθe^{2i \theta} part measures the angle so it doesn’t contribute to the length, which means that the length of z2z^2 is simply r2r^2.

If we choose cc such that c<1|c| < 1, then z2+c|z^2+c| will be greater than 11 if z2>2|z^2| > 2. If r>1r > 1 then r2r^2 will grow at each step. If you think of z2z^2 and cc as vectors, then the smallest value for rr is when z2z^2 and cc point in opposite directions.

r

In most cases, z2z^2 and cc won’t point in opposite directions, so rr will be even larger and the trajectory will escape faster, but we can check to see if z2>2z^2 > 2 and stop iterating when we get to that value.

Vectors are added together by connecting the point of uu to the tail of vv to get z=u+vz = u+v. Subtracting vv from uu turns the direction of vv around.

vector-addition

Using the FOIL method, we can calculate z2z^2 for combinations of positive and negative vectors uu and vv.

z=u+viz2=(u2v2)+2uviz=uviz2=(u2v2)2uviz=u+viz2=(u2v2)2uviz=uviz2=(u2v2)+2uvi.\begin{aligned} z = u + vi \Rightarrow z^2 = (u^2-v^2) + 2uvi\\ z = u - vi \Rightarrow z^2 = (u^2-v^2) - 2uvi\\ z = -u + vi \Rightarrow z^2 = (u^2-v^2) - 2uvi\\ z = -u - vi \Rightarrow z^2 = (u^2-v^2) + 2uvi.\\ \end{aligned}

In the picture above, you can see that the red zz vectors have the same z2z^2, as do the green vectors. After a single iteration, the new value of z=u+vz=u+v is the same as the new value of z=uvz = -u-v, just as the green vectors coincide. If u+vu+v results in an escaping trajectory, then so does uv-u-v.

Julia Sets

Gaston Julia wanted to find the boundary between the points that escaped and those that didn’t. The boundary is now called the Julia Set. The zz points are taken from all points in the complex plane, and the set changes depending on the initial choice of the parameter cc.

julia-set-c=-0.1+0.1im

In this image, c=0.7269+0.1889ic = -0.7269 + 0.1889i, and the complex numbers range over [1.6,1.6]×[1.6,1.6]i[-1.6,1.6] \times [-1.6,1.6]i. The image is composed of 1281 points in each dimension, and for each point zz, the number of iterations needed to exceed z>2|z| > 2 is stored in a matrix. The black regions are where z|z| starts greater than 22, and the colored points are counts. On the right side, the color bar indicates the number of iterations, up to 1000, before the norm of zz becomes large.

Another example shows that for c=0.20509091+0.71591ic = -0.20509091 + 0.71591i, the points either are clearly in the set or outside. At the top of this post is a plot of the Julia set for values of cc from c=0.29609091+0.62491ic = −0.29609091+0.62491i to c=0.20509091+0.71591ic = -0.20509091 + 0.71591i, taken from Benjamin Badger’s Julia Sets page.

julia-set-c=-0.20509091 + 0.71591i

The code (on Github) to generate these images is very straightforward. An iterCounts array is initialized to zero, and for each zz, the value of zn=zn12+c|z_n| = |z_{n-1}^2+c| is computed. The iteration is repeated until either zn>2|z_n| > 2 or the number of iterations reaches 1000. The array iterCounts is displayed as a heatmap to make the image.

The code is written in the Julia language because Julia is easy to learn and extremely fast. In fact, Julia is so fast that I didn’t bother to take advantage of symmetry.

In the figure above, there are 1.6 million points and each point requires one multiplication and one addition to get zn=zn12+cz_n = z_{n-1}^2 + c, and two more multiplications and an addition followed by a square root to get the norm of zz. Even if only a quarter of the points got to the limit of 1000 iterations, somewhere near one billion calculations were required. In Julia, the whole process completed in just a little over 10 seconds even with inefficient code.

function plotJulia(c,xRng = [-1.6,1.6], yRng = [-1.6, 1.6])
# Array of iteration counts
nx = 1281
ny = 1281
iterCounts = zeros(nx,ny)
# Locations
xLocs = range(minimum(xRng),stop = maximum(xRng),length = nx)
yLocs = range(minimum(yRng),stop = maximum(yRng),length = ny)
# Loop over x,y locations iterating until maxIter reached, or |z|² > 4
maxIters = 1000
for j = 1:nx
for k = 1:ny
z = xLocs[j] + yLocs[k]im
while iterCounts[j,k] < maxIters && norm(z) < 2
z = z^2 + c
iterCounts[j,k] += 1
end
end
end

# Display Julia set as a heatmap
display(heatmap(iterCounts'))

# Return array of iteratations
return iterCounts

end

Gaston Julia never saw the images that we can now quickly produce in the Julia language, but he was able to work out many of the details and provide great insight into the structure of Julia sets.


References