The double compound pendulum is a fascinating physical system that exhibits chaotic behavior. In a previous article, we solved and animated the double compound pendulum. In this follow-up article, we explore the fractal structure of the double compound pendulum.

We have seen that the double compound pendulum can flip over, and whether it does so depends sensitively on the initial conditions. In this article, we calculate the time until the first flip-over as a function of the initial conditions. We use a different solver function which is better suited for situations like this and integrate the equations of motion until the first flip-over occurs for a wide range of initial conditions.

When we plot the results, an intriguing structure will appear with the typical signs of a fractal, most notably self-similarity on many scales. A word of caution if you want to reproduce the results: it is quite computationally expensive. Creating the plots in this article required more than 100 CPU hours! But if you look at the pictures, I think you agree that it was worth the effort.

The following figure is a reminder of the nomenclature that we are using:

The following code snippet ist just for your convenience. It is the right-hand side of the equations of motion and was derived in yesterday’s article.

def f(t, y):
    dy0_dt = 6.0 / (m * L ** 2) * (2 * y[2] - 3 * np.cos(y[0] - y[1]) * y[3]) / (16 - 9 * np.cos(y[0] - y[1]) ** 2)
    dy1_dt = 6.0 / (m * L ** 2) * (8 * y[3] - 3 * np.cos(y[0] - y[1]) * y[2]) / (16 - 9 * np.cos(y[0] - y[1]) ** 2)
    dy2_dt = -0.5 * m * L ** 2 * (dy0_dt * dy1_dt * np.sin(y[0] - y[1]) + 3 * (g / L) * np.sin(y[0]))
    dy3_dt = -0.5 * m * L ** 2 * (-dy0_dt * dy1_dt * np.sin(y[0] - y[1]) + (g / L) * np.sin(y[1]))

    return [dy0_dt, dy1_dt, dy2_dt, dy3_dt]

In principle, we could solve this just like yesterday, but there we calculated the trajectories for a given time interval. However, this time, we only want to calculate until the first flip-over occurs and we don’t know in advance when that happens. So we have to change the code a bit, and in particular, we use a different solver function which is better suited for situations like this:

mport numpy as np
from scipy.integrate import ode
import time
import sys

n1, n2 = [int(a) for a in sys.argv[1:3]]
tag = sys.argv[3]
print(n1, n2, tag)

θ_1s = np.linspace(0.5 * np.pi, 0.75 * np.pi, n1)
θ_2s = np.linspace(0.5 * np.pi, 0.75 * np.pi, n2)


def compute_portrait(θ_1s, θ_2s, dt=0.1, max_t=50):
    results = []
    for θ_1 in θ_1s:
        row = []
        for θ_2 in θ_2s:
            r = ode(f)
            y_init = [θ_1, θ_2, 0, 0]
            y = r.set_initial_value(y_init, 0)

            while r.successful() and r.y[1] < np.pi and r.t < max_t:
                r.integrate(r.t + dt)

            row.append(r.t)
        results.append(row)
    return np.array(results)


start = time.time()
θ_1s, θ_2s = select_quadrant(0, 1, θ_1s, θ_2s)
results = compute_portrait(θ_1s, θ_2s)
stop = time.time()
print(stop - start)
np.save("chaos.out.{}.{}x{}".format(tag, n1, n2), results)

Here we are integrating the equations of motion until the first flip-over occurs. We do this for a wide range of initial conditions, where both θ_1 and θ_2 range from π/2 to 3π/4, as an example.

The code above is not meant to run interactively in a Jupyter notebook but as a standalone Python script. The reason is that the computations we are going to do will be quite time consuming. In fact, the results that you see in this article have consumed more than 100 CPU hours, so I ran them on the cluster that I am glad to have access to.

Let’s run this with a resolution of 40 points both for θ_1 and θ_2 and plot the resulting time of first flip-over as a function of the initial angles:

import matplotlib.pyplot as plt


def make_plot(n, tag=""):
    fig, ax = plt.subplots(figsize=(10, 10))
    if tag:
        results = np.load(f"chaos.out.{tag}.{n}x{n}.npy")
    else:
        results = np.load(f"chaos.out.{n}x{n}.npy")
    # plt.axis("off")
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_xlabel(r"initial $\theta_1$")
    ax.set_ylabel(r"initial $\theta_2$")
    ax.imshow(results, origin="lower")

make_plot(40)

This is still a quite coarse picture but we can already guess that there is some interesting structure hiding. The dark areas are the ones for which the time to first flip-over is low, wherease the bright points are those initial conditions where there was no flip-over for a long time or even none until the maximum time computed. Let’s increase the resolution to 100x100:

make_plot(100)

Now we start to see that even in the bright regions, there are dark regions again! Increase the resolution to 250x250:

make_plot(250)

Wow, the dark spots in the bright regions have turned into smaller dark regions which have about the same shape as the big dark region! That kind of self-similarity on different scales is typical for fractal structures. Now let’s go to extremes and calculate for a resolution of 1000x1000. This takes about 25 CPU hours for this picture alone.

make_plot(1000)

Isn’t that beautiful? Just for fun, let’s redo the calculations only for the lower right quadrant:

def select_quadrant(qx, qy, θ_1s, θ_2s):
    n1, n2 = len(θ_1s) // 2, len(θ_2s) // 2
    if qx == 0:
        θ_1s = np.linspace(θ_1s[0], θ_1s[n1], n1)
    else:
        θ_1s = np.linspace(θ_1s[n1], θ_1s[-1], n1)
    if qy == 0:
        θ_2s = np.linspace(θ_2s[0], θ_2s[n2], n2)
    else:
        θ_2s = np.linspace(θ_2s[n2], θ_2s[-1], n2)
    return θ_1s, θ_2s

make_plot(1000, "quad")

And again, we see the same pattern repeated on a still smaller scale.

Conclusion

In conclusion, we have explored the fractal structure of the double compound pendulum, a fascinating physical system that exhibits chaotic behavior. By calculating the time until the first flip-over as a function of the initial conditions, we discovered an intriguing structure with self-similarity on many scales, typical of fractals. It is fascinating to see how a simple physical system can exhibit such complex behavior, and how it can be explored using numerical methods. The fractal structure of the double compound pendulum is a beautiful example of the power of mathematics and computation in understanding the world around us.