This chaotic strange attractor represent the final resting positions for a magnetic pendulum suspended over some magnets (shown as black dots). It kind of looks like mixed paint. The 2D animation shows what happens as you decrease the damping factor. The 3D animation was shaded by path length.
(* runtime: 25 seconds, increase n for higher resolution *)
n = 40; h = 0.25; g = 0.2; mu = 0.07; zlist = {Sqrt[3] + I, -Sqrt[3] + I, -2I};
image = Table[z2 = z[25] /. NDSolve[{z''[t] == Plus @@ ((zlist - z[t])/(h^2 + Abs[zlist - z[t]]^2)^1.5) - g z[t] - mu z'[t], z[0] == x + I y, z'[0] == 0}, z, {t, 0, 25}, MaxSteps -> 200000][[1]]; r = Abs[z2 - zlist]; i = Position[r, Min[r]][[1, 1]]; Hue[i/3], {y, -5.0, 5.0, 10.0/n}, {x, -5.0,5.0, 10.0/n}];
Show[Graphics[RasterArray[image]], AspectRatio -> 1];
The picture on the left shows another version with five magnets. See also my physics pendulums.
Here is some Mathematica code to numerically solve this using the Beeman integration scheme with the predictor-corrector modification:
(* runtime: 45 seconds, increase n for higher resolution *)
n = 40; tmax = 25; dt = 0.1; h = 0.25; g = 0.2; mu = 0.07; zlist = {Sqrt[3] + I, -Sqrt[3] + I, -2I};
image = Table[z = x + I y; v = a = a1 = 0; Do[z += v dt + (4a - a1)dt^2/6; vpredict = v + (3a - a1)dt/2; a2 = Plus @@ ((zlist - z)/(h^2 + Abs[zlist - z]^2)^1.5) - g z - mu vpredict; v += (2a2 + 5a - a1)dt/6; a1 =a; a = a2, {t, 0, tmax, dt}]; r = Abs[z - zlist]; Hue[Position[r, Min[r]][[1, 1]]/3], {y, -5.0, 5.0, 10.0/n}, {x, -5.0, 5.0, 10.0/n}];
Show[Graphics[RasterArray[image]], AspectRatio -> 1];
Links
- Magnetic Pendulum – C++ program by Ingo Berg, uses Beeman integration scheme, here is a nice big picture and here is an animation
- Mathematica code – by Urijah Kaplan
- Butterfly Effect – systems that are highly sensitive to initial conditions
- homemade magnetic pendulum – by David Williamson
RSS feed 
Paul Nylander 
Hi Paul, nice to see you have opened your own blog. I’m always impressed when i see your work. As an addition i’d like to point all interested readers to a youtube video that explains the magnetic pendulum simulation very nicely: http://www.youtube.com/watch?v=Qe5Enm96MFQ
Regards,
Ingo
my comment it all looks realy interesting but what is it all for is there any use behind your
logig or are the colours the important part for you:if you like science then take a peak at
my foto albums on face book science is also my favourite:stay healthy and prosper:
Beautiful! Thanks for the code. I am making a real life magnetic pendulum right now. I have also adapted the code to run under a recent version of Mathematica. I have adapted it to produce increasing levels of zoom. I was expecting to find new patterns at every increased zoom level, like you would in a fractal. This would imply that no matter how precise you determine initial conditions, it will never be precise enough to determine the outcome. However, this is not what I found.
Why?
Here’s the code. Varying the “max” would determine where you are zooming in the positive quadrant.
n=30; h=0.25; g=0.2; mu=0.07; zlist={Sqrt[3]+I,-Sqrt[3]+I,-2I};
image[xmin_,xmax_,ymin_,ymax_] :=
ParallelTable[z2=z[25]/.NDSolve[{z”[t]==Plus@@((zlist-z[t])/(h^2+Abs[zlist-z[t]]^2)^1.5)-g z[t]-mu z'[t],z[0]==x+I y,z'[0]==0},z,{t,0,25},MaxSteps->20000][[1]];
r=Abs[z2-zlist];
i=Position[r,Min[r]][[1,1]];
i,
{y,ymin,ymax,(ymax-ymin)/n},{x,xmin,xmax,(xmax-ymin)/n}]
max=3.2;
min=Table[max-10^-i,{i,0,3}]//N;
images=image[#,max, #, max]&/@min;
colors={1->{1,1,1},2->{0,0,0},3->{1,0,0}};
Graphics[Raster[#/.colors],ImageSize->Small]&/@images
Really good blog, thank you very much for your effort in writing this post.