Summer project: more twirling chains

There are three main origins of this post:

1. I really like playing with twirling chains.
2. Yet another example of how mixing analytical and numeric approaches in Mathematica can be cool.
3. This seems like a cool distributed project for next summer.

The twirling chains post hinted at some aspects of that problem that I think are really interesting. What stymied me back then, though, was an inability to fully model the system. Thankfully, I finally figured out a way to do it, and now I’m excited to learn more about this interesting system.

Modeling problems

When I first tried to model a system of masses on a string, I really just approached it like a multi-pendulum. I’ve done lots of that before in 2D, but the 3D was really giving me headaches. Here’s how I set up the equations of motion:

Mathematica to set up the 3D beads equations of motion

Here are a few comments that go chronologically with the code above:

• nmax is the total number of beads minus 1 since there’s always one that’s driven
• phase[t] is put inside the Cos[] and Sin[] of the driven bead
• x[n], y[n], and z[n] are the coordinates of the nth bead, defined recursively from the bead above it.
• ke[n] is the kinetic energy of the nth bead (note that I set the mass to 1 for each)
• ketotal just adds up all the kinetic energy
• petotal adds up the gravitational potential energy of the beads
• power is used to calculate air resistance (using the power formalism in Lagrangian dynamics). 0.5 is the friction coefficient
• L is the Lagrangian
• lag is a function that calculates the Euler Lagrange equation (with the friction correction)
• eoms are the equations of motion. They need to be flattened to be joined with the initial conditions later
• inits has two definitions. Mathematica will use the last one (I used to use the first one). It turns out that there’s problems at theta and phi=0 so I offset them a little.
• alleqns puts the equations of motion and the initial conditions together.
• vars enables me to tell NDSolve what to solve for (all the thetas and phis for each bead).

One of the reasons I really like using Mathematica to do this is how it can handle the very complicated analytical equations involved. As an example, here’s the equation of motion for the 4th bead in a 4-bead chain:

Equation of motion for the 4th bead in a 4 bead chain

They just get uglier as you use more and more beads. Having Mathematica do this analytical work is great, especially since you don’t have to print the results of the equations of motion before just sticking them in the differential equation solver (NDSolve).

So what’s the problem? Well, when I plugged it all in and ran it with nmax quite low, it gobbled up all the memory and crashed without showing any results. This really frustrated me, as I didn’t see what the immediate problem was. Yes, the equations of motion are more complicated than those for the 2D case, but ultimately it’s just having to take baby steps and determine the appropriate changes for the variables. I knew that, for a single (Euler-method-type) step, you have everything worked out numerically except for the derivatives of the variables (which, strictly speaking, are both the thetas and phis and their first derivatives). Those you just need to solve for, solving n equations for n unknowns, then you take a step. Solving equations like that are what differential equations solvers are really good at, and I didn’t quite understand why Mathematica was choking on this problem.

So I decided to break the task down a little and see where the logjam was. I thought I’d solve the equations of motions in a way that had the derivatives on one side, and all the (previous step) values on the right hand side. Note that that’s the form most non-Mathematica differential equations solvers want as well. Then, I figured, I could do a simple Solve[] command and see what the problem was. Well, this sounded good, but the first step, which is really just a rearrangement of the equations, gobbled up all the memory! Ah ha! NDSolve had never actually gotten around to taking steps, it was trying to find an analytical pre-solution of sorts, and it was running out of memory. I was excited to have found the problem, but crestfallen that I wouldn’t be able to model this system.

Modeling solution!

So I did what most people did: I complained to google! And, thankfully, google (more specifically StackExchange) came through for me. It turns out you can relatively easily convince Mathematica to not try to find the analytical rearrangement mentioned above. Instead, you can have it do it with numerical values at every step. All you do in NDSolve[] is add SolveDelay->True as an option! Awesome!

Pretty pictures

Ok, so now I had a model that could work. Would it show the same results that I can get with my kids’ toys at home? You bet:

10 beads twirled at a constant rate with friction

and, of course, you can run it with lots of different bead numbers and superimpose the movies:

2-15 beads superimposed

And you can investigate whether the resonance is maintained as the frequency increases (which is what the phase[] function is set to above):

accelerated rotation with 15 beads

You can see how the chain rises and bulges out, and, sort-of, maintains its Bessel-like shape as discussed here.

Next summer

After playing with this just a little, I can tell there’s a lot of richness here. What I’d like to do is put a few of my summer research students on this next year. We’ll model all kinds of things, and build some driving mechanisms using Arduino/Mathematica. One thing I’m really interested in is whether you can decouple x- and y- motions if the amplitudes are small enough, so that you can get one mode in x-z and another in y-z happening at the same time. That’s what this paper hints at, anyways. We could look at this both using Mathematica and in the lab.

I call it a distributed project because I think others could join in if they were interested. The code is ready to be distributed and run with lots of different boundary/initial conditions. And building a setup in the lab really shouldn’t be that hard. If you’re interested, let me know!