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](http://en.wikipedia.org/wiki/double_pendulum) to flip over, as a function of initial conditions.
pic http://upload.wikimedia.org/wikipedia/commons/8/87/Double_pendulum_flips_graph.png Double pendulum fractal : Fractal generated from a double pendulum. Image credit [Jeremy S. Heyl](http://www.phas.ubc.ca/~heyl/) from UBC. You can read more about what it means at the [image description page](http://en.wikipedia.org/wiki/File:Double_pendulum_flips_graph.png).
Inspired by this, I decided to make my own picture, but for the [swinging Atwood's machine](http://en.wikipedia.org/wiki/Swinging_Atwood%27s_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:
pic FRACTAL_2400x2400s_small.png Swinging Atwood's machine fractal : 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.
# Implementation details
This was generated by the standard [fourth order Runge-Kutta method](http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Common_fourth-order_Runge.E2.80.93Kutta_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.
~~~
lang c++
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.
~~~
lang c++
int i=0, imax = 40000;
for(int x=0; x 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](http://en.wikipedia.org/wiki/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.