Mysteries of the Simulated Pendulum

Last week, I spent a bunch of time using VPython to simulate a simple pendulum, which was a fun way to fritter away several hours (yes, I'm a great big nerd), and led to some fun physics. I had a little more time to kill, so I did one of the things I mentioned as a possible follow-on, which turned out to be kind of baffling, in a good way.

Last week's post was written very quickly, and thus ended up a little more jargon-y than I usually shoot for, so let me try to set the stage a little better for this one. the physical system I'm talking about is just a simple pendulum, a mass on the end of a string swinging back and forth. The idealized version of this is very simple mathematically, and gives a very regular oscillation at a frequency that just depends on the length (which is why old clocks contain pendulums). When you dig into the details, though, there's some rich physics involved in terms of forces, because the force the string needs to exert varies over the course of the swing, getting bigger toward the bottom and smaller out at the ends (this is why, if you have small kids or are just a kid at heart and thus spend time on swingsets, you feel "heavier" at the bottom of the arc of a playground swing, and lighter out at the turning point).

What I did in the simulation was rather than using physics to figure out the force needed to make the bob swing back and forth in advance, I treated the "string" on the simulated pendulum as if it were a spring, with a force that increases as it gets "stretched." There's a nice simple mathematical form for this, characterized by a number called the "spring constant," and the bigger that constant is, the bigger the force you get for a given stretch.

This is nice because it saves me having to code equations into VPython, and also because it is, in some ways, an accurate representation of what actually happens. If you make a pendulum out of a ball and a bit of string, the string doesn't know how to solve physics equations. It supplies the appropriate force by stretching a small amount and exerting a force in response. The string you typically use to make a pendulum doesn't stretch very much, so it looks very close to the "ideal" case where the equations are easy to solve, namely a string that doesn't stretch at all. And, in fact, my toy computer model does exactly that when you make the spring constant big. Physics works, hooray.

Of course, having put together the code with an arbitrary spring constant, the fun part is sticking in parameters that aren't very close to the "ideal" case, to see what happens. In this case, the spring stretches substantially as the pendulum swings, and the path it follows makes nifty patterns-- kind of a braided look in the post from last week, or the cartoon mustache kind of thing on the left in the "featured image" at the top of this post (RSS readers, click through). That happens because I'm lazy, and didn't want to work out the details of the starting conditions, so I just imagined that the spring-string was unstretched when it was released, because that's easy to do.

When that happens, you get two things going on: a back-and-forth motion that's the "normal" pendulum track, and an in-and-out motion characteristic of a spring. You can sort of see what's going on in that left-hand image above, which is a spring with a small (25 N/m) spring constant released from an unstretched position. Right as it's let go on the right side of the swing, it falls more or less straight down, with the spring stretching a lot. As the stretch gets big enough, it pulls the bob gets pulled back in, and the combination of the sideways pendulum motion and the inward spring motion flattens out the bottom of the arc, then makes the upturn at the left end, whereupon the whole thing reverses itself.

Of course, this isn't the only possibility, and in fact, it happens only because I was lazy. If I started the spring with some stretch, though, the force it provided would be bigger sooner, and the oscillation would look more like a "normal" pendulum. This is what you see on the right above: I found a value of the stretch for which this looks much simpler, and the spring hardly stretches at all.

Now, there are two ways to do find this proper prestretch value: brute force trial and error, or using physics. What I did was a combination of the two, and that's the baffling part. I tried to use physics to estimate the starting stretch, but ended up being way wrong, and finding the actual value by trial and error.

So, what's the physics guess? Well, the simple treatment of a pendulum has the string not stretch at all, which means that the initial force along the string must be big enough to counter the portion of the pull of gravity that tends to stretch the string. This leaves just the sideways component, and drives the regular oscillation. The value of that radial component of gravity is easy to calculate, and turns out to be the weight of the bob multiplied by the cosine of the release angle.

Of course, for the simulated pendulum to work right, the string needs to stretch a little, so as to allow for the increasing force you need to make the thing work. So, I reasoned correctly that the appropriate stretch wouldn't be exactly the stretch needed to cancel the stretching part of gravity, but something a little smaller. I had no idea what the exact value should be, though, and didn't see an easy way to calculate it, so what I did was to stick in an arbitrary factor in front of that extra force when calculating the stretch. So, the pre-stretch of the spring in the simulation is given by the weight times the cosine of the angle times an arbitrary factor. I varied the arbitrary factor and looked for the value that gave the minimum amplitude of the oscillation-- that is, the track that looked most like a "normal" pendulum.

So, I did this for a bunch of different factors and a bunch of different spring constants, and like a good physicist, I made a graph:

Radial oscillation amplitude vs. pre-stretch factor, for various simulated pendula.

This plots the amplitude of the in-and-out oscillation on the vertical axis, and the arbitrary factor on the horizontal axis. And it's totally bizarre.

For one thing, my initial assumption had been that the factor should be a number less than one-- the "ideal" value that would keep the spring from stretching in the first instant of the swing would leave it too short to provide the needed force in the next instant, leading to a big oscillation. I figured that the correct value would end up being a bit smaller than that-- say, 75% of the "ideal"--allowing some spring stretch and keeping everything nice.

What you see from the graph, though, is that the correct value turns out to be substantially bigger than that-- of the five spring constants I tested, the smallest value is around 2.4-- that is, 240% of the "ideal" value. My initial thought was that this was just a numerical error on my part, but I don't see where I made a mistake.

And there's another reason to think that this isn't just a coding screw-up in my definition of the starting point, which is that the value depends on the spring constant in a complicated way. If I had just managed to get the force wrong by a factor of 2.5, you would expect that factor to be the same for each different spring, but the five spring constants I tried give different values for the factor needed to get the smallest amplitude. I can graph that, too, and it looks like this:

Pre-stretch factor for minimum amplitude vs. spring constant.

The factor gets bigger as you increase the spring constant, though it seems to reach some maximum value a bit below 5. I have no idea what causes that, though. Which is why I didn't bother to fit a function to those points, because I don't even know what I should guess as the appropriate form.

The other weird thing about the first graph is an aesthetic point-- all the data points fall on V-shaped curves. On either side of the minimum value, the points are beautifully linear-- in fact, that's how I got the data for the second graph: I fit straight lines to the data to either side of the minimum, and calculated the point where those lines crossed. (I could find the uncertainty in that value, but there's a limit to how much math I'm willing to do for a blog post, and anyway, they're close to the by-hand optimized values for the 100 N/m and 400 N/m cases, where I just plugged in factors until I found the minimum amplitude to four decimal places.) That's kind of a weird shape, though, in terms of physics graphs. You don't see a lot of natural phenomena that have that kind of pointy shape to them-- a parabola, or some other kind of smooth curve would be more "normal." You kind of need to be a physicist to be bothered by that, but it's definitely a thing that jumped out at me.

So, anyway, I'm baffled by this. It's very robust, in that I get the same basic behavior even if I make changes to the simulation code, but also very odd. I'm also not at all sure how to go about calculating the right answers. The easiest way to do it would probably involve a Lagrangian, but even there, I suspect you wouldn't get an analytical solution, but would need to do some numerical solving of equations. Which might or might not be more reliable than my crude modeling of the situation.

If I were going to teach intermediate classical mechanics at any time in the near future, this might be a great way to re-learn all that Lagrangian stuff that I forgot from my undergrad classical mechanics courses. It'll be at least six years before I do that, though, so I'm not likely to put in the effort (then again, I've spent umpteen hours screwing around with VPython, so...). I might also be able to dig up somebody else's solution-- the pendulum is, after all, a very important physical system, and I'm sure it's been studied in detail by lots of people-- but I have too much else to do to go in for a long literature search. Thus, this blog post, where I'll just throw this out there for people to either recognize, or calculate themselves, or point out the embarrassing flaw in my whole approach to this question.

In the absence of a better conclusion, then, let me just say again that this is a nice demonstration of how rich the physics of even very simple and well-studied systems can be. And also, what an enormous dork I am, spending hours and hours playing around with this to no concrete purpose...

UPDATE: The ScienceBlogs uploader rejects the code "for security reasons," but here's a Dropbox link to the poorly documented VPython program I used to do these simulations.