# Fractal basins of attraction in a Magnetic Pendulum

I am trying to write a Mathematica program that realizes a graphical approximation of the basins of attraction in a Magnetic pendulum subject to friction and gravity, in which the three magnets are disposed on the vertices of an equilateral triangle. This system is chaotic and has very interesting properties. The basins of attraction look something like this:

A code I wrote can produce a $400 \times 400$ image such as this (Caution: no fancy ColorFunctions involved) in about two hours. The computation seems to be extremely slow.

Is there any way of having a better rendering, say, full HD 1920x1080 resolution, for the basins of attraction of a magnetic pendulum as the one mentioned that can be run in a farely quick time on a common machine?

## Code

Here is the code I used to produce the above image. I set the position of the magnets and define the Lagrange equations

X1 = 1; X2 = -(1/2); X3 = -(1/2); Y1 = 0; Y2 = Sqrt[3]/2; Y3 = -(Sqrt[3]/2);
X[1] = X1; X[2] = X2; X[3] = X3; Y[1] = Y1; Y[2] = Y2; Y[3] = Y3;
Eqs[k_,c_,h_]:={
x''[t]+k x'[t]+c x[t]-Sum[(X[i]-x[t])/(h^2+(X[i]-x[t])^2+(Y[i]-y[t])^2)^(3/2),{i,3}]==0,
y''[t]+k y'[t]+c y[t]-Sum[(Y[i]-y[t])/(h^2+(X[i]-x[t])^2+(Y[i]-y[t])^2)^(3/2),{i,3}]==0
}


I define a function that numerically integrates the equations up until $t=100$.

Sol[k_, c_, h_, xo_, yo_] :=
NDSolve[
Flatten[{Evaluate[Eqs[k, c, h]],
x'[0] == 0, y'[0]== 0, x[0] == xo, y[0] == yo}],
{x, y}, {t, 99.5, 100.5}, Method -> "Adams"
];


I define a function tt that gives a value between $\frac13, \frac 23, 1$ based on magnet proximity at time $100$ for fixed $k,c,h$ (in this case $.15$,$.2$,$.2$) and a function k that evaluates tt on a grid.

tt = Compile[{{x1, _Real}, {y1, _Real}}, Module[{},
Final = ({x[100], y[100]} /. (Sol[0.15, .2, .2, x1, y1])[[1]]);
Distances = Map[(Final - #).(Final - #) &, {{1, 0}, {-(1/2), Sqrt[3]/2}, {-(1/2), -(Sqrt[3]/2)}}];
Magnet = Min[Distances];
Position[Distances, Magnet][[1, 1]]/3]];
k[n_, xm_, ym_, xM_, yM_] := ParallelTable[tt[xi, yi], {yi, ym, yM, Abs[yM - ym]/n}, {xi, xm, xM, Abs[xM - xm]/n}];


Finally, I rasterize the table produced by k.

G = Graphics[Raster[k[400, -2, -2, 2, 2], ColorFunction -> Hue]]


and, after a while, I obtain the previous image. I attempted using a dynamic energy control (i.e. using EvaluationMonitor to monitor the energy level of ther trajectory: if it falls in a potential hole NDSolve throws the position) but this did not increase the speed as much as I was hoping; it actually seems to slow the computation down.

• Why don't you start by sharing the code you wrote? Perhaps we can improve on that, rather than reinventing the wheel. – MarcoB Jul 18 '16 at 0:30
• You might want to look up previous work by Paul Nylander. – J. M. will be back soon Jul 18 '16 at 9:06
• There is a syntax error (Y3 = -(Sqrt3]/2);) in the first code block. When I run Sol[0.15, .2, .2, 1, 1] I get an error message saying that the system is undetermined, because of lack of spaces between c and x and y. – C. E. Jul 18 '16 at 11:57
• After modernizing Nylander's code for the current version, I managed this. – J. M. will be back soon Jul 18 '16 at 16:37
• @J.M. This seems very interesting. Can you share your code with us? How much time did the computation take? – Lonidard Jul 18 '16 at 17:46

JM commented:

If you want to try things out, use Nylander's second snippet, which is using a Beeman integrator. This looks to be faster than native NDSolve[] for this specific case.

Paul Nylander's code is here.

Below is a modified version of his code which computes all points simultaneously using the fact that all the operations in Beeman's algorithm are Listable functions in Mathematica.

The run time for the 400x400 image is around 30 seconds.

n = 400;
{tmax, dt} = {25, 0.05};
{k, c, h} = {0.15, 0.2, 0.2};
{z1, z2, z3} = N@Exp[I 2 Pi {1, 2, 3}/3];
l = 2.0;

z = DeveloperToPackedArray @ Table[x + I y, {y, -l, l, 2 l/n}, {x, -l, l, 2 l/n}];
v = a = aold = 0 z;
Do[
z += v dt + (4 a - aold) dt^2/6;
vpredict = v + (3 a - aold) dt/2;
anew = (z1 - z)/(h^2 + Abs[z1 - z]^2)^1.5 + (z2 - z)/(h^2 + Abs[z2 - z]^2)^1.5 +
(z3 - z)/(h^2 + Abs[z3 - z]^2)^1.5 - c z - k vpredict;
v += (5 anew + 8 a - aold) dt/12;
aold = a; a = anew,
{t, 0, tmax, dt}];
res = Abs[{z - z1, z - z2, z - z3}];
Image[0.2/res, Interleaving -> False]


• complexGrid from here might be useful. – J. M. will be back soon Jul 18 '16 at 23:55
• Turn the Do into Table and return z each time to see a trippy animation: cl.ly/3C163h1A3F26 – Chip Hurst Jul 19 '16 at 15:04
• @Chip, I believe Nylander did something like that in one of the animations on his website. – J. M. will be back soon Jul 19 '16 at 15:08
• @J.M. Ah-ha! I should've clicked the link! – Chip Hurst Jul 19 '16 at 15:09

I don't have any breakthrough ideas, but I am able to cut the computation time in half on my computer by optimizing the usage of NDSolve.

My version of your code looks like this:

X[1] = 1;
X[2] = -(1/2);
X[3] = -(1/2);
Y[1] = 0;
Y[2] = Sqrt[3]/2;
Y[3] = -Sqrt[3]/2;

Sol[k_, c_, h_, xo_, yo_] := NDSolve[{
x''[t] + k x'[t] + c x[t] - Sum[(X[i] - x[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i, 3}] == 0,
y''[t] + k y'[t] + c y[t] - Sum[(Y[i] - y[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i, 3}] == 0,
x'[0] == 0,
y'[0] == 0,
x[0] == xo,
y[0] == yo
}, {x, y}, {t, 99.5, 100.5}, Method -> "Adams"
];

nf = Nearest[{{1, 0}, {-0.5, Sqrt[3]/2}, {-0.5, -Sqrt[3]/2}} -> Automatic] /* First;
getBasin[x1_, y1_] := nf[{x[100], y[100]} /. Sol[0.15, .2, .2, x1, y1] // Flatten, 1]


You are essentially building your own Nearest function, but one is already built in. This approach is as performant as your code.

Some advanced usage tips for NDSolve are documented here. The idea is that before NDSolve can start to integrate the equations it needs to rewrite them, and this takes time. Instead of doing the rewriting for every single point we can do it just once and then use that for all of the points.

getStateData[k_, c_, h_, x0_, y0_] :=
First@NDSolveProcessEquations[{
x''[t] + k x'[t] + c x[t] - Sum[(X[i] - x[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i, 3}] == 0,
y''[t] + k y'[t] + c y[t] - Sum[(Y[i] - y[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i, 3}] == 0,
x'[0] == 0,
y'[0] == 0,
x[0] == x0,
y[0] == y0
}, {x, y}, t, Method -> "Adams"]

sd = getStateData[.15, .2, .2, 1, 1];

getBasin2[x0_, y0_] := Module[{state = sd, sol},
state = First@NDSolveReinitialize[state, {x[0] == x0, y[0] == y0}];
NDSolveIterate[state, 100.5];
sol = {x[100], y[100]} /. NDSolveProcessSolutions[state];
nf[sol]
]


Let's test it:

ArrayPlot[
ParallelTable[getBasin2[xpos, ypos], {xpos, -2, 2, 0.1}, {ypos, -2, 2, 0.1}],
ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}
] // AbsoluteTiming
`

This simple example took 22 seconds for me to generate, as opposed to 44.5 seconds for my first rewritten version of your code.

Your image I was able to generate in 33 minutes rather than the two hours it took for you:

Implementing a stopping condition seems cumbersome, but you can also achieve speed enhancements by lowering the amount of integration time. Integrating to 100 seems excessive, with 25 it looks the same for the 400x400 case and it takes only 11 minutes.