Lagrange multipliers

This is a post that builds on my previous notes about the calculus of variations. This week I’m going to teach about modeling constraint forces using the Lagrangian approach, and I wanted to put these notes down for my current students and my future self.

Multiple variables

If you have a system with multiple variables, it’s pretty straightforward to convince yourself that the action is now the integral of a function of several functions and their slopes. Each variable will have a best path that minimizes the functional, with each best being surrounded by independent deviations that must all make the integral worse (larger, if you’re looking for a minimum). You can parametrize those deviations (the separate shapes of the deviations) with the same \alpha that we used before. When you minimize the integral according to \alpha you get

\frac{dJ}{d\alpha}=\int_{x_1}^{x_2}\left[\left(\frac{\partial f}{\partial y}-\frac{d}{dx}\frac{\partial f}{\partial y'}\right)\eta_1(x)+\left(\frac{\partial f}{\partial z}-\frac{d}{dx}\frac{\partial f}{\partial z'}\right)\eta_2(x)\right]\,dx

Now, if the deviations are truly independent, then the only way to make that be zero is to have the two things in the parentheses to be zero independently. Basically, you just get an Euler equation for both y and z.

Constraints

Ok, but what if the deviations aren’t independent (\eta_1 and \eta_2)? In other words, what if there are some additional constraints on y and z such that they need to conspire to keep the whole system on some sort of curve? We call that an additional constraint. It’s a pain because the trick we did above (setting both parentheses to zero) won’t work.

There’s two major ways of dealing with such a constraint. The first is to reparametrize your system, recognizing that you don’t really have two degrees of freedom, but rather one. If you do that, you reduce the problem to one dimension and you’re off and running. This is actually one of the big selling points of the Lagrangian approach to mechanics problems. The constraints in that case are usually constraining forces (like normal or tension forces) and it’s nice to be able to ignore them, especially since they usually don’t do any work (they force the constraint, which means they act perpendicular to the constraining surface, which the particle never leaves, hence no work).

The second way to deal with a constraint is what this post is all about. We want to be explicit about how the two \eta‘s are related to each other (if one changes, the other has to change in such a way as to make sure the constraint is still met). Maybe we can express one in terms of the other, being able to factor one out of the above equation, leaving behind something that has to be zero. That’s the game plan.

Buckle up, this gets a little nasty, but it’s worth it. In the end we’ll get a way of leveraging the constraint rules to be able to obtain not only the motion of the particle we care about, but also information about the forces ensuring the constraint in the first place. That latter information is lost in the first approach mentioned above.

First we have to have a way to mathematically talk about the constraint. Here’s the approach I’ll take: We’ll assume there’s a combination of y’s and z’s that must always be a constant: g(y, z)=constant. Note that often people like to put that as g(y,z)=0 by redefining g to have “-constant” in the definition. For how I’m going to use it, that’s not strictly necessary.

Next we need to force g not to change, even if \alpha changes. That’s easy, we’ll just say that dg/d \alpha is zero:

\frac{dg}{d\alpha}=\frac{\partial g}{\partial y}\frac{\partial y}{\partial \alpha}+\frac{\partial g}{\partial z}\frac{\partial z}{\partial \alpha}=\frac{\partial g}{\partial y}\eta_1(x)+\frac{\partial g}{\partial y}\eta_2(x)=0

where I’ve made use of some knowledge about the \eta‘s. Cool, I can solve that for \eta_2 in terms of \eta_1:

\eta_2(x)=-\frac{\frac{\partial g}{\partial y}}{\frac{\partial g}{\partial z}}\eta_1

Alright, now we’re in business. We can plug that into the first equation above and then be able to factor out \eta_1 out of the mix. Whatever’s leftover in the integrand has to be zero since \eta_1 can be anything. That leaves us with:

\left(\frac{\partial f}{\partial y}-\frac{d}{dx}\frac{\partial f}{\partial \dot{y}}\right)-\left(\frac{\partial f}{\partial z}-\frac{d}{dx}\frac{\partial f}{\partial \dot{z}}\right)\frac{\frac{\partial g}{\partial y}}{\frac{\partial g}{\partial z}}

(sorry for the ugly pixelization, wordpress didn’t want to render it for some reason)

That’s pretty ugly, but, with some manipulation, we get a cool separation of variables:

\left(\frac{\partial f}{\partial y}-\frac{d}{dx}\frac{\partial f}{\partial \dot{y}}\right)\frac{1}{\frac{\partial g}{\partial y}}=\left(\frac{\partial f}{\partial z}-\frac{d}{dx}\frac{\partial f}{\partial \dot{z}}\right)\frac{1}{\frac{\partial g}{\partial z}}=\lambda(x)

where I’ve set both sides to a function that’s only a function of x (not y or z). Classic separation of variables approach.

Ok, then, we seem to have the newly modified Euler equation. For each variable we now have:

\frac{\partial f}{\partial y}-\frac{d}{dx}\frac{\partial f}{\partial y'}-\lambda(x)\frac{\partial g}{\partial y}=0

Oof, done. Now, what can we do with it?

Simple pendulum

Let’s apply this to a simple pendulum and see what we can do. First, we have to shift from the typical calculus of variations variables (y(x) and z(x)) to typical Lagrangian variables (x(t) and y(t)). The Euler-Lagrange (with constraint) equations will be:

\frac{\partial L}{\partial x}-\frac{d}{dx}\frac{\partial L}{\partial x'}-\lambda(t)\frac{\partial g}{\partial x}=0

and

\frac{\partial L}{\partial y}-\frac{d}{dx}\frac{\partial L}{\partial y'}-\lambda(t)\frac{\partial g}{\partial y}=0

where L is the Lagrangian (Kinetic energy – potential energy).

Ok, for a pendulum, we have x and y with a constraint that says the distance from the tether point (a convenient origin) should be the length of the string, a (I would love to use “ell” but it looks too much like a 1). So, a good constraint function would be:

g(x,y)=\sqrt{x^2+y^2}

Plugging and chugging gives us the following equations of motion:

m x''(t)=-\lambda(t) \frac{x}{\sqrt{x^2+y^2}}

and

m y''(t)=-mg-\lambda(t) \frac{y}{\sqrt{x^2+y^2}}

(please take note of the nasty double use of g in that derivation. Sorry about that).

If you stare at those long enough, you’ll see that \lambda(t) is really just the negative of the tension force of the string. That’s cool, and representative of what \lambda often turns out to be. Don’t be too bugged by the sign of the constraint force. Often we care most about when it goes to zero.

So, what did we learn? If you build in the constraints, you get equations of motion for more variables than the degrees of freedom of the system for the benefit of the ability to solve for the constraint force. Cool!

Typical analytical use

My friend Rhett Allain has a great post about one of the most famous uses of lagrange multipliers. There he talks about ice sliding off a bowl, trying to figure out when the ice leaves contact with the spherical surface of the (upturned) bowl, and also trying to figure out things about the constraint force by modeling it as a partial spring.

I have a comment there talking about how I used to teach this. I would walk students through the analytical argument, showing that you can find the angle where the constraint force (or \lambda(t)) goes to zero. But then I tried to just put the equations of motion into Mathematica to numerically integrate the solutions for x(t), y(t), and \lambda(t), and I found that you can’t start with the initial conditions that are recommended in the typical analytical approach, namely having the ice start right at the top of the bowl. Mathematica dutifully shows that the ice will just sit there if you start it there. You have to nudge it off just a little to make any interesting progress.

But really my problem with the typical analytical approach is that there’s such a small number of “perfect” problems that provide such a unique end. Just like the rest of the Lagrangian stuff I’ve done with students, I really wanted to show that this approach will work, even for hard problems that don’t have analytical results.

Mathematica break through

Unfortunately, over the last few years of toying with this in Mathematica, I’ve been stymied. What I’ve tried to do is to put the constraint Euler-Lagrange equations, the initial conditions (that meet the constraint), and the constraint equations into an NDSolve command, hoping to numerically integrate the variables and the constraint forces. However, the most common problem was this:

NDSolve::nderr: Error test failure at t == 0.013181091203435103`; unable to continue.

Arg, that’s incredibly frustrating. As far as I could ever figure, it would numerically integrate for basically one time step (note the early time when it failed) and decide that the constraint wasn’t being enforced anymore. I dug around a lot to try to figure out what was going on, but I never got any satisfaction.

But this year, as I was preparing the resources for my students for this material, I decided to take one more crack at it, and I’m glad I did, because I finally solved the problem!

I coded up a simple pendulum, giving an E-L equation for x and y, a constraint that said \sqrt{x^2+y^2}-1 had to be zero, and some initial conditions that had the pendulum up at some angle. Mathematica failed as usual, so I decided to look more carefully at the equations. I took out a sheet of paper, and I came to realize that the constraint equation is used by making x”(t) and y”(t) have a relationship, not so much using the direct relationship between x and y. So, I tried something. Instead of putting in the constraint, I defined the constraint as a function (namely g(x,y)=\sqrt{x^2+y^2}-1 and I told NDSolve that the second time derivative of g with respect to t should be zero. This gave NDSolve direct information about how x” and y” are related. And lo and behold . . . IT WORKED! Since then I’ve been using that same trick, even for problems with multiple constraints, and I’m very satisfied with the results.

Ice on a Mathematica bowl

Here’s how I can use Mathematica to do the ice on an upturned bowl:iceonblockmma

From which I can plot the constraint force:iceblockconstraintplot

It starts with a magnitude of 9.8, which is to be expected, since it needs to balance gravity near the top of the bowl. The most interesting feature of the curve is where it crosses zero. The bowl can only have a negative constraint force (it can only push, it can’t pull), so when it goes to zero, that’s when the ice leaves the surface. Using FindRoot[\lambda[t]/.sol, {t, 1.5}] we can find the time when that happens, and we can find the angle (arctan(x/y)) when that happens (~48 degrees).

Some oddities

If you look carefully at how I’m using the constraint, you’ll notice that it is never used without taking a derivative of it. That means that any constants involved in the definition of g are not used. This bothered me at first, so I looked to see what happened if I gave initial conditions that didn’t match g. In the case of the ice block, the result was that the ice stayed on a circle with whatever initial radius I gave it in the form of the initial conditions. Cool, I guess. Basically it means that you establish the functional form of the constraint to get the derivatives right (ie, the square root in the ice block case, showing the relationship between x and y), but you use initial conditions to set the  particulars.

Bead on a wire

Ok, let’s flex our muscles a little. What about a bead on an oddly shaped wire:beadonwiremma

Note that the only changes from the previous code is to change g and to change the initial conditions. Here’s the result:beadonwiremovie

Here’s a plot of the normal force as a function of time:beadonwireconstraintplot

Since the mass is set to 1, you can see there are times when the wire has to provide almost 30 g-forces to maintain its shape. Cool.

Time-dependent constraints

What about a particle on a floor that starts to tilt?tiltingplanemma

producing thistiltingplanemovie

Multiple constraints

Finally, here’s an example of a multi-pendulum, having constraints on the distance between each successive mass. When there are multiple constraints, the extension of the Euler Lagrange equation is

\frac{\partial L}{\partial y}-\frac{d}{dx}\frac{\partial L}{\partial y'}-\sum_i^N\lambda_i(t)\frac{\partial g_i}{\partial y}=0

where there are N constraints. For each one, you have to specify both g and provide for a different \lambda.

Here’s the code for a n-pendulum, constraining each length to be 1:npendulummma

which makes this animation:npendulummovie

And allows us to see the tensions necessary to pull this off:npendulumconstraintplot

Note in this plot that there’s a couple of points when some of the constraint forces change sign. At that point the strings are being asked to push instead of pull. In other words, this wouldn’t work with string. You’d have to do it with light-weight bars or something that can both push and pull.

Can you tell I’m happy I finally got this to work in Mathematica? I’m looking forward to working with my students on this on Wednesday!

Advertisements

About Andy "SuperFly" Rundquist

Professor of physics at Hamline University in St. Paul, MN
This entry was posted in mathematica, physics, teaching. Bookmark the permalink.

5 Responses to Lagrange multipliers

  1. Wes says:

    Thanks for this! I’ve been stuck with the same issue and the answer is so simple! Do you understand exactly why this works though? Is NDSolve flawed? Or perhaps the standard treatment of Lagrangian Mechanics with constraints is? What gives?

    • Andy "SuperFly" Rundquist says:

      I think that the notion that the constraint has to be satisfied doesn’t give mathematica enough freedom to adjust the variables. Telling it that the constraint doesn’t change in time seems to do the trick. So mathematically it’s basically the same, it just took me a while to convince mathematica to give it a shot.

  2. I like the valuable information you provide in your articles.
    I’ll bookmark your weblog and check again here regularly. I am quite
    certain I’ll learn a lot of new stuff right here!
    Best of luck for the next!

  3. unnamed says:

    wow! super!

  4. Pingback: Lagrange multipliers revisited | SuperFly Physics

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s