Magnetic Pendulum Fractal Basin Boundaries in Mathematica

[majinsock]

Homework Statement

Hey all, this is for a laboratory. I need to determine the fractal basin boundaries for a magnetic pendulum swinging chaotically about 3 other magnets underneath. I was able to plot a single path for the pendulum in Mathematica, but now I need to expand my code to determine each endpoint for each starting point. I'm thinking that a Do loop involving lists might work, but for the life of me I can't figure out how to incorporate the functions NDSolve, List and Do. I COULD do each path individually and record the results, but I really don't think that's practical.

Homework Equations

y''[t] + r y'[t] - Ʃ(Yi-y[t])/√((Xi-x[t])^2 + (Yi-y[t])^2 + d^2)^3 + c y[t] ==0

The sum is for Y1, Y2 and Y3 (and X1, X2 and X3) which are the coordinates of the magnets. The other equation for the displacement of x[t] is the same as for y[t], so I don't see the need to write it twice.

r,c,d are all constants, where r is air resistance, d is the distance between the pendulum and the sitting magnets if it were directly over it. c is the spring constant of the pendulum.

I've attached a couple notebooks so that you can see more clearly what I'm doing.

The Attempt at a Solution

My attempt at the solution is to make a Do loop in which Mathematica uses NDSolve for each starting coordinate from {-5,-5} to {5,5} in increments of 0.1 and outputs the endpoint at around 100 seconds when the pendulum is basically circling a single magnet. I understand what I have to find out, but for the life of me I can't figure out how to make it work!

 

[Nate810]

Any help would be greatly appreciated! Here's what I have so far: Do[ {X1, Y1} = {-5, -5}; {X2, Y2} = {0, 5}; {X3, Y3} = {5, -5}; eqns = {y''[t] + r y'[t] - Sum[(Yi - y[t])/Sqrt[(Xi - x[t])^2 + (Yi - y[t])^2 + d^2]^3 + c y[t] == 0, {i, 1, 3}], x''[t] + r x'[t] - Sum[(Xi - x[t])/Sqrt[(Xi - x[t])^2 + (Yi - y[t])^2 + d^2]^3 + c x[t] == 0}, {x[0], y[0]} = {0, 0}, {x'[0], y'[0]} = {xstart, ystart}, {x[100], y[100]}]; soln = NDSolve[eqns, {x[t], y[t]}, {t, 0, 100}]; {xend, yend} = {x[100], y[100]} /. soln; Print[{xstart, ystart} -> {xend, yend}], {xstart, ystart} = {xstart + 0.1, ystart + 0.1}, {xend, yend} < {5, 5}]