The system of differential equation that this post hinges on is the following:

\(x_{n + 1} = a_0 + a_1x_n + a_2x_n^2 + a_3x_ny_n + a_4y_n + a_5y_n^2\)

\(y_{n + 1} = a_6 + a_7 x_n + a_8 x_n^2 + a_9x_ny_n + a_{10}y_n + a_{11}y_n^2\)

The code looks super simple on Fronkonstin's blog. It consists of the differential equation function definitions, a set of coefficients, a function that generates the steps iteratively, and a plotter. These pieces map to python well

The Equations

def attractor(x, y, z):
    a = z[0] + z[1] * x + z[2] * x**2 + z[3] * x * y + z[4] * y + z[5] * y**2
    b = z[6] + z[7] * x + z[8] * x**2 + z[9] * x * y + z[10] * y + z[11] * y**2
    return [a, b]

The Coefficients

z1 = [1.0, -0.1, -0.2,  1.0,  0.3,  0.6,  0.0,  0.2, -0.6, -0.4, -0.6,  0.6]
z2 = [-1.1, -1.0,  0.4, -1.2, -0.7,  0.0, -0.7,  0.9,  0.3,  1.1, -0.2,  0.4]
z3 = [-0.3,  0.7,  0.7,  0.6,  0.0, -1.1,  0.2, -0.6, -0.1, -0.1,  0.4, -0.7]
z4 = [-1.2, -0.6, -0.5,  0.1, -0.7,  0.2, -0.9,  0.9,  0.1, -0.3, -0.9,  0.3]
z5 = [1.0, 0, -1.4, 0, 0.3, 0, 0, 1, 0, 0, 0, 0]

The function that generates the steps iteratively

def make_galaxy(z_input, iters=1000):
    output = [[0, 0]]
    for i in range(1, iters):
        new_vals = attractor(*output[i - 1], z = z_input)
        output.append(new_vals)
    df = pd.DataFrame(output)
    df.columns = ['x', 'y']
    return df

the plotter

For the plotter, I'm going to use the color palette that I generated for this blog in this post [Color Palette Generation]​:

colors = ['#a252f2', '#ae6bf2', '#f9bf76', '#d1313d', '#e5625c', '#8eb2c5', '#4dc9e3']
colors_extended = colors + ['#4dc9e3', '#1f3f6c', '#a96e14', '#a129df', '#00ff9f']

def plot_galaxy(df, color_palette=colors_extended, save_location=None):
    color_vals = [choice(color_palette) for _ in range(len(df))]
    _, ax = plt.subplots()
    ax.scatter(df['x'], df['y'], c=color_vals, alpha=np.random.uniform(0.5, 0.9), s=3)
    plt.tick_params(left = False, right = False, labelleft = False,
    labelbottom = False, bottom = False)
    if save_location:
        ax.set_axis_off()
        plt.savefig(save_location,
                    bbox_inches='tight',
                    pad_inches=0,
                    transparent=True)
    else:
        plt.show()

Now we have all the pieces, and can generate some plots!

for z in [z1, z2, z3, z4, z5]:
    plot_galaxy(make_galaxy(z))

Z1 Coefficients: [1.0, -0.1, -0.2, 1.0, 0.3, 0.6, 0.0, 0.2, -0.6, -0.4, -0.6, 0.6]

Z2 Coefficients: [-1.1, -1.0, 0.4, -1.2, -0.7, 0.0, -0.7, 0.9, 0.3, 1.1, -0.2, 0.4]

Z3 Coefficients: [-0.3, 0.7, 0.7, 0.6, 0.0, -1.1, 0.2, -0.6, -0.1, -0.1, 0.4, -0.7]

Z4 Coefficients: [-1.2, -0.6, -0.5, 0.1, -0.7, 0.2, -0.9, 0.9, 0.1, -0.3, -0.9, 0.3]

Z5 Coefficients: [1.0, 0, -1.4, 0, 0.3, 0, 0, 1, 0, 0, 0, 0]

Z5 is the Hénon map, a special case of our more general map.

After generating the plots on the blog, I hoped to begin tweaking the coefficients and generating some new pictures.

Fronkenstin gives the words of encouragement:

Changing the vector of parameters you can obtain other galaxies. Do you want to try?

Chaotic systems must have these properties:

• it must be sensitive to initial conditions
• it must be topologically transitive
• it must have dense periodic orbits

The idea was to tweak initial parameters and get some cool plots. Unfortunately, just choosing random initial parameters is not guaranteed to give us an interesting output. Only a small subset of the possible sets of coefficients will give us an interesting strange attractor output. Most of the sets of coefficients will either blow up to \(\pm \infty\), or return nothing interesting at all.

I tried out some random sets of the 12 coefficients in the hopes that I would find something, but that was inefficient and returned mostly garbage. The overwhelming majority of images were just random splashes of color on a dark background.

Here are some examples of the majority of plots produced:

Not too thrilling!

After a while of tinkering, I ended up constraining the space in which I generated the coefficients \(z\) and got some better results:

Here is the script for generating the more interesting outputs:

valid_galaxies = []
for _ in range(10000):
    # generate array
    z = [round(np.random.uniform(-1.2, 1.2), 1) for _ in range(12)]
    try:
        df = make_galaxy(z)
        # check if out of bounds
        a, b = df.max().to_list()
        if 0.9 < abs(a - b) < 1:
            valid_galaxies.append(df)
        else:
            pass
    except OverflowError:
        pass

plot_galaxy(choice(valid_galaxies))

You can tweak the values that \(|a - b|\) lies in to get different outcomes.

Most of the results look like the underwhelming sparse dots of color on a dark background. Out of 10,000 iterations, the algorithm above yielded 70 results that fit the criteria and didn't result in an overflow error.

Of those 70 results, only a handful looked interesting.

Here are some of the more interesting ones:

The next batch was ran with 0.75 < abs(a - b) < 1.5

Of the 10,000 iterations, the algorithm found 451 results.

In the book, the author uses a coding system to represent different coefficient values. It looks like this:

Then he outputs the key at the top of each plot:

In this case, he leads with an E as a code for a quadratic 2-dimensional map and the rest of the letters translate to coefficients in the key.

For example, J = -0.3, T = 0.7, S = 0.6, M = 0.0, B = -1.1, and so on until you have all 12 coefficients for this map. This becomes handy further in the book as the degree of the polynomials increases (and allows for even more complexity of output).

If I do a follow-up post, this method could be used to make more intricate strange attractors (with higher degree polynomials)

def letter_range(start, stop, step=1):
    """Yield a range of letters"""
    for ord_ in range(ord(start), ord(stop), step):
        yield chr(ord_)

def translate_to_coefs(input_str: str) -> List[float]:
    coefs = merge([{l: d / 100} for l, d in zip(letter_range('A', '['), range(-120, 140, 10))])
    # drop the E that leads the string
    return [coefs[a] for a in list(input_str)[1:]]

def plot_galaxy_with_key(key, color_palette=colors_extended, save_location=None):
    df = make_galaxy(translate_to_coefs(key))
    color_vals = [choice(color_palette) for _ in range(len(df))]
    fig, ax = plt.subplots()
    plt.title(key, loc='left', color='#f9bf76')
    ax.scatter(df['x'], df['y'], c=color_vals, alpha=np.random.uniform(0.5, 0.9), s=3)
    ax.set_facecolor('#1B1D21')  # this is the background color of my blog :D
    plt.tick_params(left = False, right = False, labelleft = False,
    labelbottom = False, bottom = False)
    if save_location:
        ax.set_axis_off()
        plt.savefig(save_location,
                    bbox_inches='tight',
                    pad_inches=0,
                    transparent=True)
    else:
        plt.show()


codenames = ["EDSYUECINGQNV", "EELXAPXMPQOBT", "EEYYMKTUMXUVC", "EZPMSGCNFRENG",
             "ENNMJRCTVVTYG", "EOUGFJKDHSAJU", "EQKOCSIDVTPGY", "EQLOIARXYGHAJ",
             "ETJUBWEDNRORR", "ETSILUNDQSIFA", "EUEBJLCDISIIQ", "EVDUOTLRBKTJD"]


for code in codenames:
    plot_galaxy_with_key(code)

Unfortunately, random choosing these values doesn't seem to be super helpful (we tried that above), but it does allow us to recreate some of the nice figures from the book with a simple codename:

What a neat title! That could be a band name.

What I want is to find a way to consistently generate bulbous (roundish, relatively self-contained) strange attractors.

How do we measure the bulbousness of a strange attractor?

Did it work? kind of. Here are the plots in descending order of the bulbous measure: