Orbiting the Hénon Attractor

Note: This was originally published as an Observable notebook.

The Hénon strange attractor is a system of two simple equations that gives rise to chaotic behavior.

I first encountered Jared Tarbell's beautiful visualizations of the attractor many years ago and recently explored making my own pictures using my favorite programming language at the time, Julia.

Here's a gallery of the most interesting pictures I created along the way. The code for these is described at the bottom of the post.

Click on any picture to see it in full resolution.

Behind the Scenes

Here's the Julia code that I wrote to explore this system.

1. Introducing the attractor

To explore a chaotic attractor we need to define it in code. The Hénon map has an angle parameter so we define henon(α) to return the attractor function for that angle.

# define the attractor
function henon(α)
  c = cos(α)
  s = sin(α)
  attractor((x, y)) = (x*c - (y - x^2)*s, x*s + (y - x^2)*c)
end

We'll also need some utilities to help us along the way.

First, a function to generate random points within our specified bounds, so that we can create detailed zoom images of specific parts of the phase space.

# randomly initialize an point within our bounds, then
# iterate it a few times to remove initial noisy jumps
function init(h, (xlo, xhi), (ylo, yhi))
  init_x = xlo + (xhi - xlo)rand()
  init_y = ylo + (yhi - ylo)rand()
  p = (init_x, init_y)
  for _ in 1:100
    p = h(p)
  end
  p
end

Second, a linear scale function, similar to d3.scaleLinear, which interpolates values in the interval to produce integers from to that we can use for array indexing.

# map from continuous (lo, hi) to integer 1:n
function lerp((lo, hi), n)
  @assert lo < hi
  x -> 1 + round(Int, (x - lo) / (hi - lo) * (n - 1))
end

Finally, a simple fuzz function to randomly jitter points. This technique is responsible for the organic/analog feeling of the images, particularly in deep zooms.

function fuzz(x, amount)
  x + amount * (rand() - 0.5)
end

2. Computing orbits

Now we get to the heart of the system — it samples Hénon orbits and returns a 2D array of counts indicating how many times each pixel bin was visited by an orbit.

function orbit_counts(
  α;                                 # the henon angle parameter
  sz = 1000,                         # size of the output canvas
  n_samples = 100_000,               # number of total points we want to record on our canvas
  max_samples_per_orbit = 50_000,    # the maximum number of samples for any individual orbit
  x_bounds = (-1.0, 1.0),            # the x zoom bounds  
  y_bounds = (-1.0, 1.0),            # the y zoom bounds
  fuzz_factor = 0.0                  # a fuzz factor to jitter points before plotting
)
  # create a henon attractor with our given angle
  h = henon(α)
  
  # linear scales, to integer coordinates in A
  ix, iy = lerp(x_bounds, sz), lerp(y_bounds, sz)
  
  # output array of visit counts for each bin
  counts = zeros(sz, sz)
  
  # loop until we collect the desired number of attractor samples
  count = 0
  while count < n_samples
    # initialize a point
    p = init(h, x_bounds, y_bounds)

    for _ in 1:max_samples_per_orbit
      # iterate the attractor
      p = h(p)

      # if the orbit diverges, break out of the loop for this orbit
      if any(isinf, p) || any(isnan, p)
        break
      end

      # if the point is within bounds, increment the count in the corresponding bin
      if x_bounds[1] ≤ p[1] ≤ x_bounds[2] && y_bounds[1] ≤ p[2] ≤ y_bounds[2]
        x = clamp(ix(fuzz(p[1], fuzz_factor)), 1, sz)
        y = clamp(iy(fuzz(p[2], fuzz_factor)), 1, sz)
        counts[y, x] += 1
        count += 1
      end
    end
  end

  # return the array of counts
  counts
end

3. Making pictures

Next, let's generate six 1,000 × 1,000 arrays of orbit counts to visualize:

orbits = [
  orbit_counts(
    1.5957,
    sz = 1000,
    n_samples = 1_000_000,
    max_samples_per_orbit = 100_000,
    x_bounds = (.45, .55),
    y_bounds = (0.05, .15),
    fuzz_factor = 0.0001
  )
  for _ in 1:6
]

To make our pictures more interesting we generate multiple "layers" of counts and composite them together into the final image. This lets us take advantage of the variation between layers to make our pictures more interesting.

To talk about colors and color spaces we're going to use the Julia Images package:

using Images

We define one last utility function, after which we're ready to visualize the Hénon attractor:

# normalize an array of positive numbers to (0, 1), being
# careful not to introduce NaNs via division by zero.
normalize(A) = A ./ max(1e-6, maximum(A))

Finally, the pièce de résistance! The visualize function takes a collection of 2D orbit arrays as input and produces an image as output:

function visualize(orbits)
  # normalize each layer of orbit counts to (0, 1), then
  # concatenate them to create a 3D array of weights
  weights = cat(normalize.(orbits)..., dims=3)

  # try visualizing each pixel based on the difference between its 
  # minimum value and maximum value across all of the orbit count layers.
  # there are other approaches; this is a simple one to start with.
  A = normalize(maximum(weights, dims=3) .- minimum(weights, dims=3))
  
  # we're going to use the LCH color space.
  # lightness goes from 0-100
  L = 100(1 .- A)

  # contrast goes from 0-300ish
  C = 300A

  # hue cycles with a period of 360 degrees
  H = 360normalize(sum(weights .^ 3, dims=3)) 
  
  # create an array of colors by combining our lightness, contrast, and hue
  colors = RGBA.(LCHuv.(L, C, H), 1)
  
  # we want to return a 2D array. By this point our array is
  # width × height × 1, so we can just drop the third dimension
  dropdims(colors, dims=3)
end

# visualize the six orbit count matrices we computed
visualize(orbits)

All of the pictures in this post were made with variations of this code.

Interesting directions

That's it. If you make something cool, let me know!

Yuri Vishnevsky
July 2020