Discrete Lagrangian Mechanics

I'm tired of thinking about one set of technical thoughts, so it's time to take a break and think about a different set of technical thoughts. These are technical thoughts, but since trying to typeset mathematics in HTML is generally about as pleasant as self-performed dental surgery, I will try to keep myself to a minimum of notation. Perhaps I will also achieve a minimum level of incomprehensibility in the process?

There are three main frameworks for writing the equations of classical mechanics. Students who passed high school physics are familiar with Newton's formulation, the oldest of the three. In Newton's equations of motion for a particle, the central equation is

F=ma,

where the force law F may depend on the particle position, but classically does not depend on the velocity (F is a function of v if the particle experiences friction).

Less familiar are the Lagrangian and Hamiltonian formulations. In Lagrange's formulation, also known as the principle of least action, the equations of motion are related to a scalar function called the action. The action S is the integral over [0,tfinal] of the Lagrangian function L(q,v), where q is position and v is velocity. The particle follows the path q(t) which minimizes the action [1]. To find such a path, we use the technique taught in elementary calculus courses: differentiate the function with respect to the free variable, then set the derivative to zero. The only complication -- which is not so much of a complication, really -- is that we are differentiating with respect to q(t), a variable which is itself a function. Taking derivatives with respect to functions is a topic in the fancily titled topic of the calculus of variations. Anyhow, if we write down the stationarity condition (derivative = 0), we find that the minimizing function q satisfies the Euler-Lagrange equation:

D₁ L(q,v) - d/dt (D₂ L(q,v)) = 0

where D_i indicates partial differentiation in the ith position.

The Lagrangian is usually written as L(q,v) = T(v)-U(v), where T(v) = 1/2 mv² is the kinetic energy and U(v) is the potential energy. Subsituting this form into the Euler-Lagrange equation, we have

-dU/dv = ma,

which is precisely Newton's law with F=-dU/dv. That is, the Newtonian formulation and the Lagrangian formulation are equivalent; at issue, then, is which is more convenient for analysis.

To solve differential equations on a computer, it is necessary to approximate them by finite systems of difference equations. There are different ways to approach this approximation, but one of the most intuitive is simply to take the original equations of motions and replace derivatives by differences. For example, we might follow Euler and write

q'(t) = (q(t+h)-q(t))/h + O(h).

As h becomes small, the divided difference becomes a good approximation to the derivative. This scheme will give good results if we approximate t at sufficiently closely spaced points (i.e. if h is small), and if the resulting recurrences satisfy certain algebraic properties (the recurrences need to be stable). Different approximations to q'(t) lead to different numerical integration schemes, each of which has its own advantages and disadvantages.

An alternative approach to discretizing the equations of motion is to use a discrete Lagrangian formulation. Instead of having a Lagrangian function that depends on q and v, we now have a function that depends on the value at consecutive time steps: L(q_k,q_k+1). And rather than having an action integral, we have an action sum. The minimization procedure is the same, however, and yields the discrete Euler-Lagrange equation:

D₁ L(q_k,q_k+1) + D₂ L(q_k-1,q_k) = 0

Unsurprisingly, if we choose the discrete Lagrangian to approximate the continuous Lagrangian, the numerical methods that we obtain from solving the discrete Euler-Lagrange equations are often the same as the methods we obtain from directly approximating Newton's equations of motion.

So why would we choose to use a discrete Lagrangian formulation rather than discretizing Newton's law? It is because the discrete Lagrangian formulation can guide our choice of numerical methods. We define the momentum of the continuous system to be

p = D₂ L(q,v).

Those who have been in a course on classical mechanics may vaguely remember that this definition is the first step in setting up a Legendre transformation, by which we can change from a Lagrangian formulation to a Hamiltonian formulation. For the discrete Lagrangian formulation, though, there are two natural definitions for the momentum:

p⁺_k = D₂ L(q_k-1,_k)
p^-_k = -D₂ L(q_k,_k+1).

This is one of the things that makes discrete analogues to continuous problems tricky: in general, well-defined classical concepts like velocity and momentum have multiple reasonable analogues in the discrete world.

So we have two reasonable definitions for momentum in the discrete world, and similarly we end up with two copies of other related objects (e.g. there are two discrete canonical one-forms). But, by an almost magical stroke, there is only one canonical two-form associated with the discrete equations. Thus, the discrete Euler-Lagrange equations for a conservative system mimic the continuous Euler-Lagrange equations in that they conserve symplectic structure (i.e. they conserve the canonical two-form). To this day, I have no good intuitive explanation of what symplectic forms are and why they're important; nevertheless, they are important, because with the exact conservation of symplectic structure under discretization comes better approximate conservation of better-known quantities like energy and overall momentum. This behavior markedly differs from the behavior of many direct discretizations of Newton's law in which the energy of the discrete system tends to either grow or decay over time.

So why am I thinking of this now? Some time ago, I was challenged by one of my professors, W. Kahan, to show that a numerical scheme for integrating a particular type of differential equation -- a matrix Riccati equation -- would maintain a symmetry property present in the original differential equation. He had in mind something enormously complicated involving the group-theoretic structure of the problem, which he had once proved but never trusted, and which was lost in the morass of papers in his office. I just attacked the problem with algebra guided by intuition, and to our mutual surprise, found a proof. My proof is short (one page), and relies on a set of coordinate changes to go from unknowns which are morally symmetric -- symmetric but for the error imposed by discretization -- to unknowns which are symmetric even in the discrete equations. A few months later (or earlier?), I gave a brief summary of discrete Lagrangian mechanics in a computational mechanics seminar that I was attending. While rummaging through old files today, I found both my notes on discrete Lagrangian mechanics and my proof of the symmetry in Kahan's integrator. It seems to me that there is a philosophical connection betweek Kahan's integrator and the integrators derived from the discrete Euler-Lagrange equation: in both cases, the discretization is set up in such a way that some interesting properties can't help but be conserved. I wish I'd thought a little longer and made that connection more precise; alas, that thinking will have to wait for another day.

1. Technically, the action need not be minimized, but need only be rendered stationary. This distinction is important in the study of dynamics of constrained systems, where the action is typically not minimized.