This post got its spark when I read this

John had also mentioned the list in last week’s Global Physics Department meeting, and I decided to check out the list of great questions. My problem is, I got stuck on the first one. It asks what pattern three snails would leave if they start at the vertices of an equilateral triangle, and then head toward the snail at the next vertex.

I thought it was a cool problem, though I don’t actually think it’s a physics problem, since there’s attraction, but no Newton’s third law. It’s really a math problem, but I don’t have any problem with that. In fact, after playing with it for a while, now I’m really intrigued, and I think I need your help in figuring out some of the subtleties.

## Equations of motion

Having said that bit about not being physics, it might surprise you to hear I used it in a physics class I’m teaching right now. We’re learning about equations of motion; how to determine them, what they mean, why we care. I thought this might be a cool problem to tackle, so we ran with it.

Each snail has its velocity direction determined by the relative position of another snail. For now, we set the speed to be a constant. We decided that the equation of motion for snail number 1 is given by:

,

with similar equations for the other snails. What else do we need? Well, it’s a set of first-order differential equations so, unlike, say, a Newton’s law problem, we only need initial positions, not initial velocities. That’s easy, because the problem states they start at the vertices of the triangle.

Cool, we ought to be able to plug those six differential equations into Mathematica and have some fun. Six, you say? Well, the equation is a vector equation with two components, and we have three equations like it. Ok, here’s the result:

## Up a notch

Ok, that’s cool, but what if we step it up a notch to see other shapes, or more snails? Here’s 20 snails randomly placed:

Cool, huh? They go into a circle. If I let the simulation go longer, that circle shrinks down to a point. First question: is that point the center of mass? Something else? I’m not sure.

## Next notch

Ok, that’s fun, but then things got weird. I wondered what would happen in 3D. I figured they’d go to a point again, but would they go into a circle? Take a look:

Ok, here’s my big question: They do tend to end in a single plane, what does that plane represent about the initial placements of the snails?

I think I get why they end in a single plane. Basically they’re dragging each other around, always trying to catch up with the next guy. If they aren’t in a plane, as the next snail gets dragged over the corner that takes it out of the plane, it’ll cut the corner and flatten it. I think this continues until you’re in a plane. I just can’t decide how that plane is determined. Thoughts?

By the way, the pics seem cooler to me if you change the speed from being constant to being proportional to the separation between the snails. Here’s an example with 20 snails:

So, to recap, here are my questions:

- Do they tend to the center of mass of the initial distribution (in 2D or 3D)?
- Does the limit plane in 3D represent something of the original distribution, and, if so, what?

Here’s the Mathematica for the 3D stuff. Quite easy to shift it to 2D if you’re interested:

Don’t laugh too much at that Thread/@Thread command. I was in a hurry, and it worked. I know I bang this gong a lot, but I just can’t help noting how easy this was to extend to an arbitrary number of particles. You set “n” and then all the equations of motion are crafted and integrated. I’ve run it with 100 and, even on my puny netbook, it only takes 1 second to do the integration (it takes 10 to do the graphics with that many lines, but on my work computer even that’s pretty fast).

If you make velocity proportional to the displacement (like you do at the end) the dimensions decouple so there is no difference between 3-D and three separate 1-D problems. Better yet since it’s now a linear problem you can find the matrix that corresponds to differentiation and calculate eigenvalues which all have non-positive real parts (you need one neutrally stable mode for the trivial solutions where the snails all start in the same place) which correspond to decaying or neutrally stable modes. Some theorem from linear algebra (http://mathworld.wolfram.com/GershgorinCircleTheorem.html) says you are stuck with in a radius 1 circle centered at -1+0i in the complex plane if you choose the constant of proportionality to be unity (the circle and the origin scale with your choice of proportionality constant so the stability is the same) and all your eigenvalues must have a non-positive real part.

I’m with you about the linear problem, but I’m not sure I get what the eigenvalues will get me. It sounds a lot like diagonalizing the inertia tensor, which would lead to preferred directions.

The eigenvectors describe the spatial distribution, it’s easier to think about for each 1-D problem seperately. One (the only?) zero mode is where all the x values are the same or [1,1,1,1,…] which gives you translational invariance (removing this mode is equivalent to setting the initial center of mass at the origin). After long times only the two (one if the eigenvalue is real, as in the N=2 case) eigenmodes have any appreciable magnitude, those with the lest negative real part, and we can acurately describe the motion in that dimension with them. Repeating the process for y and z is easy since they have the same matrix form and therefore the same eigenvalues and eigenvectors for their derivative.

The plane of the final spiralis determined by the relative amplitudes and phases of the slowest decaying modes in the x,y and z components. The motion is in each co-ordinate is a decaying sinusoid of the same frequency so in three dimensions it’s a spiral. I can’t see apriori why there can’t be two modes with similar decay rates, in fact since the eigenvalues are bounded in a small part of the complex plane the spiral should take longer to become aperent for many snails since the eigenvalues will fill the space (or possibly become degenerate but I don’t think the determinant is zero for any N).

P.S. if you assert that your particles are equidistant on a circle of radius R you can actually find analytic decay rates and angular velocities for the fully non-linear problem

d R / d t = speed * sin( pi /N )

d theta / d t = (speed/R) *cos ( pi / N )

So in the limit of large N the spiral is very slow to decay and in the limit of small N the configuration rotates more slowly, though the behaviour at small R has some issues 😉