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.

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:

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:

• Observe that much of the image is a solid black colour. The black region is the slowest to compute since it has to iterate all the way to imax for every pixel. We can avoid this by randomly sampling some pixels uniformly. Then when actually computing the image, if all the randomly sampled pixels in the surrounding area are the same colour, chances are that this pixel will be that colour too… and thus the expensive evaluation of that pixel can be avoided. I actually used repeated jittered grid sampling for this instead of uniform sampling.
• Multithreading. I made my program run on 8 threads (for my Core i7 2600k). Apparently 8 threads runs faster than 4 threads even though the CPU only has 4 physical cores (but has hyperthreading). Anyway the program is embarrassingly parallel since each pixel can be evaluated independently and performance scales linearly to number of cores. I was tempted to implement this on the GPU but it is too much work.
• In the derivatives function, change calls to powl to just multiplication. Store sine and cosine values in a lookup table (provided it is large enough to be precise).

Ultimately it is still a compromise between accuracy (smaller step size) and speed (larger step size). Anyway the 2400 × 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.