Swinging Atwood’s machine fractal

Here’s an interesting fact: most chaotic systems can generate a fancy picture that may resemble a fractal. You can choose almost any aspect of the system to plot, and the result will probably be quite interesting. For example, here is a picture of the time taken for a double pendulum to flip over, as a function of initial conditions.

Double pendulum fractal
FIGURE 1 Fractal generated from a double pendulum. Image credit Jeremy S. Heyl from UBC. You can read more about what it means at the image description page.

Inspired by this, I decided to make my own picture, but for the swinging Atwood’s machine. The swinging Atwood’s machine is simply two masses connected by a string on a pulley, where the lighter mass is allowed to swing. Surprisingly such a simple system is chaotic too! The shape of the trajectory of the system can be uniquely defined by just the ratio between the masses and the initial angle of the swinging mass. The initial radius does not matter, since it essentially just scales the entire trajectory without affecting its shape.

So, I plotted the time taken for the swinging Atwood’s machine to flip over, as a function of its initial angle (horizontal axis, ranging from 0 on the left to pi on the right), and mass ratio (vertical axis, ranging linearly from 1 at the top to 3.4 on the bottom). The system is said to have “flipped over” if it crosses the theta = 0 line, i.e. the positive y axis. Here’s the picture:

Swinging Atwood's machine fractal
FIGURE 2 Fractal generated from a swinging Atwood’s machine.

Unfortunately it doesn’t look quite as fancy as the double pendulum one. Maybe it is because of my choice of colour scheme. Here, light blue means it flips over quickly; dark blue means it takes a long time to flip over; black means it never flips over. Notice the some jagged things around the edges. This is due to error in the RK4 numeric scheme I’m using. Because near the regions of sharp transition, the trajectory gets close to the origin, the radius becomes small and the equations of motion becomes close to singular. Naturally, computers (or maybe just me) are not good at dealing with singularities so it becomes imprecise.

1 Implementation details

This was generated by the standard fourth order Runge-Kutta method. In fact, for this problem, a symplectic integrator would be more suitable since that would exactly conserve the total amount of energy. However, I was too lazy to implement that.

inline void derivatives(long double *y, long double *dy, long double mu, long double g) {
    // evaluate the equations of motion
    dy[0] = y[2]/(1+mu); //r
    dy[1] = y[3]/powl(y[0],2); //theta
    dy[2] = powl(y[3],2)/powl(y[0],3) - g*(mu-cosl(y[1])); //pr
    dy[3] = -g*y[0]*sinl(y[1]); //ptheta

void rk4(long double *x, long double mu, long double g) {
    /* fourth order Runge Kutta integrator */
    long double dx[4], yt[4], dym[4], dyt[4];
    derivatives(x, dx, mu, g); // k1
    for(int i=0; i<4; i++) {
        yt[i] = x[i] + h/2 * dx[i];
    derivatives(yt, dyt, mu, g); // k2
    for(int i=0; i<4; i++) {
        yt[i] = x[i] + h/2 * dyt[i];
    derivatives(yt, dym, mu, g); // k3
    for(int i=0; i<4; i++) {
        yt[i] = x[i] + h * dym[i];
        dym[i] += dyt[i];
    derivatives(yt, dyt, mu, g); // k4
    for(int i=0; i<4; i++) {
        x[i] += h/6*(dx[i]+dyt[i]+2*dym[i]);

So, basically, for each pixel, we set the initial conditions and run rk4 up to some set maximum limit, and count how many iterations has passed before the trajectory crosses over the y axis.

int i=0, imax = 40000;
for(int x=0; x<xsize; x++) {
    for(int y=0; y<ysize; y++) {
        long double q[4];
        q[0] = 500; // initial radius
        q[1] = y/(long double)ysize * PI; // initial angle
        q[2] = 0; // initial radial momentum
        q[3] = 0; // initial angular momentum

        long double T = q[1], Told = q[1];
        for(i=0; i<imax; i++) {
            T = q[1];
            if( Told < PI && Told > PI/2 && T>PI ||
                Told > PI && T < PI && T > PI/2 || q[0] <= 0) break;
            Told = q[1];
            rk4(q, mu, g);
        data[x][y] = i;

This is a very straightforward algorithm, but it is also highly inefficient. Here were some things I did to make the program run faster:

Ultimately it is still a compromise between accuracy (smaller step size) and speed (larger step size). Anyway the 2400 &times; 2400 picture still took the following amount of time to compute:

    real    2052m6.549s
    user    9538m28.860s
    sys     6m33.584s

Yet I am still not satisfied with the level of accuracy. Oh well. The picture didn’t look all that interesting anyway.