John Peach

Mind-sized STEM ideas and experiments, beyond the textbook.

Julia-sets-in-Julia

Home

keywords: Julia sets, Julia language, complex numbers, iterated functions, fractals

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)} f(x)=Q(x)P(x)

where P ( x ) P(x) P(x) and Q ( x ) Q(x) Q(x) are polynomial functions. Q ( x ) Q(x) Q(x) can't be zero, but it can be 1 1 1, meaning that f ( x ) = P ( x ) 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. f(x)=x+1.

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

f ( x ) = { x 2 if m o d     ( x , 2 ) = 0 3 x + 1 if m o d     ( 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} f(x)={2x3x+1ifmod(x,2)=0ifmod(x,2)=1.

In other words, if x x x is even, divide it by 2 2 2, but if it's odd, multiply by 3 3 3 and add 1 1 1. Any number that's a power of 2 2 2 immediately collapses down to 1 1 1, for example, x 0 = 8 ⇒ f ( 8 ) = 4 ⇒ f ( 4 ) = 2 ⇒ f ( 2 ) = 1 x_0 = 8 \Rightarrow f(8) = 4 \Rightarrow f(4) = 2 \Rightarrow f(2) = 1 x0=8f(8)=4f(4)=2f(2)=1. Then you get stuck in the loop 4 , 2 , 1 , 4 , 2 , 1 , … 4,2,1,4,2,1, \ldots 4,2,1,4,2,1, If you start with 3 3 3 you get the sequence { 3 , 10 , 5 , 16 , 8 , 4 , 2 , 1 } \{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 1 1 1? Nobody knows for sure, but it's thought that every sequence eventually winds down to 1 1 1. 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 ) ( x − 3 ) (x+2)(x-3) (x+2)(x3) using the FOIL method, you multiply the First two terms, x × x = x 2 x \times x = x^2 x×x=x2, then the Outer terms, x × − 3 = − 3 x x \times -3 = -3x x×3=3x, the Inner ones, 2 × x = 2 x 2 \times x = 2x 2×x=2x, and finally the Last two, 2 × − 3 = − 6 2 \times -3 = -6 2×3=6 to get

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

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}) (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 x x x such that x × x = x 2 = − 1 x \times x = x^2 = -1 x×x=x2=1. But, in the 17th century, René Descartes named the solution to − 1 \sqrt{-1} 1 as i i i, 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 + 3 i 4 + 3i 4+3i which is 4 4 4 parts real and 3 3 3 parts imaginary. You can plot complex numbers as points, using the x − x- xaxis for the real part, and the y − y- yaxis for the imaginary part. Here's a plot of z 1 = 4 + 3 i z_1 = 4 + 3i z1=4+3i:

The Complex Plane

Suppose you wanted to multiply z 1 z_1 z1 by another complex number z 2 = 2 + 5 i z_2 = 2 + 5i z2=2+5i. Just use the FOIL method but with the complex numbers,

z 1 × z 2 = ( 4 + 3 i ) × ( 2 + 5 i ) = 4 ⋅ 2 + 4 ⋅ 5 i + 3 i ⋅ 2 + 3 i ⋅ 5 i = 8 + 20 i + 6 i + 15 i 2 = ( 8 − 15 ) + 26 i = − 7 + 26 i . \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} z1×z2=(4+3i)×(2+5i)=42+45i+3i2+3i5i=8+20i+6i+15i2=(815)+26i=7+26i.

Remember that i 2 = − 1 i^2 = -1 i2=1, so multiplying the two imaginary parts together gives the negative of the product, ignoring the i i i'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 z 1 z_1 z1 is 4 4 4, which is drawn as the vector u u u, and the imaginary part is 3 3 3, drawn as the vector v v v. 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 = 4 2 + 3 2 = 25 = 5 r = \sqrt{4^2+3^2} = \sqrt{25} = 5 r=42+32 =25 =5. The sides u u u and v v v can be written in terms of r r r and θ \theta θ as

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

A Most Beautiful Equation

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

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

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

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

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

f ′ ( θ ) = d d θ e − i θ ( cos ⁡ θ + i sin ⁡ θ ) = e − i θ ( − sin ⁡ θ + i cos ⁡ θ ) − i e − i θ ( cos ⁡ θ + i sin ⁡ θ ) = e − i θ ( − sin ⁡ θ + i cos ⁡ θ ) − e − i θ ( − sin ⁡ θ + i cos ⁡ θ ) = 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} f(θ)=dθdeiθ(cosθ+isinθ)=eiθ(sinθ+icosθ)ieiθ(cosθ+isinθ)=eiθ(sinθ+icosθ)eiθ(sinθ+icosθ)=0.

Since the derivative of f ( θ ) = 0 f(\theta) = 0 f(θ)=0 for every θ \theta θ, then f ( θ ) f(\theta) f(θ) 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 x − x- xaxis 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 θ=0, where f ( θ ) = 1 f(\theta) = 1 f(θ)=1 since e − i ⋅ 0 = 1 e^{-i \cdot 0} = 1 ei0=1, cos ⁡ 0 = 1 \cos 0 = 1 cos0=1, and sin ⁡ 0 = 0 \sin 0 = 0 sin0=0. In other words, f ( θ ) = 1 f(\theta) = 1 f(θ)=1 for every θ \theta θ, which completes the proof.

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

Angles in the Complex Plane

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

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

Physicist Richard Feynman said that this is "the most remarkable formula in mathematics" because it contains e e e, the base of natural logarithms, the imaginary number i i i, π \pi π, the constant relating the radius of a circle to its circumference and area, 1 1 1 the first natural number, and 0 0 0 the center of the mathematical universe.

Iterations on a Complex Number

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

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

If you start with a smaller c c c, say c = − 0.79 + 0.15 i c = -0.79 + 0.15i c=0.79+0.15i, the values of z z z 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 z 0 z_0 z0 and additive constant c c c.

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

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

several-iterates

Adding 0.1 0.1 0.1 shifts the starting points a little bit to prevent overlapping trajectories, and dividing by 2 2 2 is equivalent to setting r = 1 2 r=\frac{1}{2} r=21 in the u u u and v v v vectors except for the slight shift of the starting point. If you set z 0 = 0.1 + e i θ z_0 = 0.1+e^{i\theta} z0=0.1+eiθ, 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 ∣ = u 2 + v 2 . \text{norm}(z) = |z| = \sqrt{u^2+v^2}. norm(z)=z=u2+v2 .

After a single iteration, z z z becomes z 2 + c z^2 + c z2+c, and the length or norm is ∣ z 2 + c ∣ |z^2+c| z2+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} θ=30°=0.52 radians (without dividing by 2 2 2), 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.744 e 10 ] . 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 2 2 2, 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]. 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.128 0.128 0.128 (after rounding).

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

z 2 = r e i θ × r e i θ = r 2 × ( e i θ ) 2 = r 2 e 2 i θ . \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} z2=reiθ×reiθ=r2×(eiθ)2=r2e2iθ.

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

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

r

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

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

vector-addition

Using the FOIL method, we can calculate z 2 z^2 z2 for combinations of positive and negative vectors u u u and v v v.

z = u + v i ⇒ z 2 = ( u 2 − v 2 ) + 2 u v i z = u − v i ⇒ z 2 = ( u 2 − v 2 ) − 2 u v i z = − u + v i ⇒ z 2 = ( u 2 − v 2 ) − 2 u v i z = − u − v i ⇒ z 2 = ( u 2 − v 2 ) + 2 u v i . \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} z=u+viz2=(u2v2)+2uviz=uviz2=(u2v2)2uviz=u+viz2=(u2v2)2uviz=uviz2=(u2v2)+2uvi.

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

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 z z z points are taken from all points in the complex plane, and the set changes depending on the initial choice of the parameter c c c.

Julia Set (c = 0.7269 + 0.1889i)

In this image, c = − 0.7269 + 0.1889 i c = -0.7269 + 0.1889i c=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 [1.6,1.6]×[1.6,1.6]i. The image is composed of 1281 points in each dimension, and for each point z z z, the number of iterations needed to exceed ∣ z ∣ > 2 |z| > 2 z>2 is stored in a matrix. The black regions are where ∣ z ∣ |z| z starts greater than 2 2 2, 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 z z z becomes large.

Another example shows that for c = − 0.20509091 + 0.71591 i c = -0.20509091 + 0.71591i c=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 c c c from c = − 0.29609091 + 0.62491 i c = −0.29609091+0.62491i c=0.29609091+0.62491i to c = − 0.20509091 + 0.71591 i c = -0.20509091 + 0.71591i c=0.20509091+0.71591i, taken from Benjamin Badger's Julia Sets page.

Julia Set (c = -0.2 + 0.7i)

The code (on Github) to generate these images is very straightforward. An iterCounts array is initialized to zero, and for each z z z, the value of ∣ z n ∣ = ∣ z n − 1 2 + c ∣ |z_n| = |z_{n-1}^2+c| zn=zn12+c is computed. The iteration is repeated until either ∣ z n ∣ > 2 |z_n| > 2 zn>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 z n = z n − 1 2 + c z_n = z_{n-1}^2 + c zn=zn12+c, and two more multiplications and an addition followed by a square root to get the norm of z z z. 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

Exploring the Infinite Beauty of Julia Sets

The journey through the world of Julia sets reveals the intricate beauty and complexity that arise from simple mathematical iterations. By exploring the boundary between points that escape to infinity and those that remain bounded, we uncover stunning fractal patterns that captivate mathematicians and artists alike.

If you'd like to go deeper, consider experimenting with different values of the complex parameter c c c to generate your own unique Julia sets. The Julia programming language, with its speed and ease of use, provides an excellent platform for these explorations. You could also look at the relationship between Julia sets and the Mandelbrot set, another fascinating fractal that serves as a map of all possible Julia sets.

You might also consider investigating the applications of complex dynamics in fields such as computer graphics, chaos theory, and even natural phenomena. The interplay between mathematics and art in fractals offers endless possibilities for creative expression and scientific inquiry.

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.


Image credits


Code and Software

JuliaSetPlots.jl - Plots trajectories and filled Julia sets: plotIterates, plotJulia

  • Julia - The Julia Project as a whole is about bringing usable, scalable technical computing to a greater audience: allowing scientists and researchers to use computation more rapidly and effectively; letting businesses do harder and more interesting analyses more easily and cheaply.
  • VSCode and VSCodium - Visual Studio Code is a lightweight but powerful source code editor which runs on your desktop and is available for Windows, macOS and Linux. It comes with built-in support for JavaScript, TypeScript and Node.js and has a rich ecosystem of extensions for other languages and runtimes (such as C++, C#, Java, Python, PHP, Go, .NET). VSCodium is a community-driven, freely-licensed binary distribution of Microsoft’s editor VS Code.
  • Pluto.jl - Pluto is an environment to work with the Julia programming language. Easy to use like Python, fast like C.

References

📬 Subscribe and stay in the loop. · Email · GitHub · Forum · Facebook

© 2020 – 2025 Jan De Wilde & John Peach