PHY5210 W08

Lagrangian Approach: The Double Pendulum

Lagrangian Approach for Two Carts and Three Springs

Already done in lecture.

The Double Pendulum

U1 = m1gL1(1 - cosφ1)
U2 = m2g[L1(1 - cosφ1 + L2(1 - cosφ2)]

In the small angle approximation, the total potential energy is

U(φ1, φ2) = ½(m1 + m2)gL1φ1² + ½m2gL2φ2².

For the first mass, the position is (x1, y1) = L1(sinφ1, 1-cosφ1), and kinetic energy is T1 = ½ m1L1²φ1dot².

The position of the second mass is (x2, y2) = (x1, y1) + L2(sinφ2, 1-cosφ2). The square of the velocity is v2² = x2dot² + y2dot² = L1²φ1dot² + L2²φ2dot² + L1L2φ1dotφ2dot(cosφ1cosφ2 + sinφ1sinφ2). The final trigonometric terms can be written as cos(φ1 - φ2) ≈ 1 for small angles.

Putting these together, we can write the kinetic energy as

T = ½(m1 + m2)L1²φ1dot² + ½m2L2²φ2dot² + m2L1L2φ1dotφ2dot.
And the Lagrangian is T - U. Lagrange's equations yield, for φ1
(m1 + m2)L1²φ1ddot + m2L1L2φ2ddot = -(m1 + m2)gL1φ1
and for φ2
m2L1L2φ1ddot + m2L2²φ2ddot = -m2gL1φ2.
These can be written in our standard matrix form
Mφddot = -Kφ


φ = [φ1, φ2]


M = [((m1 + m2)L1², m2L1L2), (m2L1L2, m2L2²)]


K = [((m1 + m2)gL1, 0), (0, m2gL1)].

Try the same trick again to solve this equation: try the solution φ = aeiωt, and get the characteristic equation (K - ω²M)a = 0 which must be solved for ω and a.

The Double Pendulum

Mφddot = -Kφ


φ = [φ1, φ2]


M = [((m1 + m2)L1², m2L1L2), (m2L1L2, m2L2²)]


K = [((m1 + m2)gL1, 0), (0, m2gL1)].

Try the same trick again to solve this equation: try the solution φ = aeiωt, and get the characteristic equation (K - ω²M)a = 0 which must be solved for ω and a.

Equal Lengths and Masses

As a simple example, consider the case when the masses and lengths are equal. This simplifies K and M to

M = mL²[(2, 1), (1, 1)]


K = mL²[(2ω0², 0), (0, ω0²)]

where ω0 = √(g/L) is the natural frequency of a pendulum of length L. The determinant is

m²L4[2(ω0² - ω²)² - ω4] = m²L44 - 4ω0²ω² + 2ω0²] = 0

This has solutions

ω1² = (2 - √2)ω0²  and  ω2² = (2 + √2)ω0²

The normal modes are found by substituting back into the characteristic equation. For ω1 we get

(K - ω²M)a = mL²ω0²(√2 - 1)[(2, -√2), (-√2, 1)][a1, a2] = 0,

implying that √2 a1 = a2. For ω2 we get

(K - ω²M)a = mL²ω0²(√2 + 1)[(2, √2), (√2, 1)][a1, a2] = 0,

implying that -√2 a1 = a2.

In the first mode, the masses oscillate together, but the amplitude of the lower mass is √2 times larger than the amplitude of the upper mass. In the second mode, the masses oscillate 180° out of phase, again with the lower mass swinging with an amplitude √2 times greater than the upper mass.

The General Case

Consider the general case of a system with n degrees of freedom oscillating with small amplitude about a point of stable equilibrium. Assume the system is holonomic so that the n degrees of freedom translats to n generalized coordinates q1, ..., qn. It is convenient to represent these as the column matrix q = [q1, ..., qn]. (Recall my note about the difference between tensors and matrices. The same applies between "true" vectors and column matrices. We often call q a vector, meaning it is a column matrix, not a "true" vector.)

We assume that the system is conservative and has a potential energy function

U(q1, ..., qn) = U(q).

The potential energy is given by

T = ½ ∑a mava²

where the sum is over N particles in the system. The N velocities are related to derivatives of the generalized coordinates. Begin by writing the position of particle a as a function of the generalized coordinates, ra = ra(q1, ..., qn). The time derivative of ra is found using the chain rule

va = dra/dt = (∂ra/∂q1)q1dot + ... + (∂ra/∂qn)qndot = ∑i(∂ra/∂qi)qidot

The square of the velocity is the square of the above expression. We want to be careful to keep track of all the terms. Just writing down the product we have

va² = vava = ∑i,j(∂ra/∂qi)⋅(∂ra/∂qj)qidot qjdot

where we must take the dot product of the vector quantities in parentheses.

Now that we have the most general equations for this situation, we want to specialize to the case of small oscillations about a position of stable equilibrium. Small oscillations imply that any generalized coordinate or its derivatives are of order ε where ε represents a small quantity. We want to keep only the lowest order terms in ε, and we will see that the first order terms in ε are zero, so the lowest order terms are proportional to ε², that is, terms like qi², qiqj, and qidot qjdot.

Let's begin with the potential energy. We assume there's a position of stable equilibrium, and for convenience, adjust the generalized coordinates so that this occurs at q = 0 (this notation means that every qi is zero). Next, we can expand the potential and kinetic energies about this point in a Taylor series. The kinetic energy can depend on the qi through the partial derivative terms. Since each term in the sum for the kinetic energy contains a factor of qidot qjdot, we take only the constant term in the Taylor series expansion of the partial derivatives, so that

T ≈ ½∑a mai,j(∂ra/∂qi)0⋅(∂ra/∂qj)0qidot qjdot

The subscript zero on the partial derivatives indicates that they are evaluated at q = 0. These terms are therefore constants. It is convenient to change the order of the sums

T ≈ ½ ∑i,j [∑a ma(∂ra/∂qi)0⋅(∂ra/∂qj)0]qidot qjdot

Then note that the sum over a just adds up constant partial derivative terms, and these can be moved out of sight by defining the matrix M with components

Mij = ∑a ma(∂ra/∂qi)0⋅(∂ra/∂qj)0

allowing us to write the kinetic energy as

T ≈ ½ ∑i,j Mijqidot qjdot = ½ qdotT M qdot

where I've explicitly shown the matrix expression for the record, though you'll normally use the sum notation.

The Taylor series expansion of the potential energy is a bit more straightforward. First, note that, at a point of stable equilibrium, all the forces are zero, that is, -∂U/∂qi = 0. Since we're expanding around a point of stable equilibrium, the term that is first order in qi is zero, and we must go to the second order terms, yielding

U ≈ U(0) + ½ ∑i,j (∂²U/∂qi∂qj)0 qi qj = U(0) + ½ ∑i,j Kij qi qj = ½ qT K q,

where I've introduced the matrix K with components Kij = (∂²U/∂qi∂qj)0.

With these approximations for the kinetic and potential energies, the Lagrangian is

L = T - U = ½ ∑i,j Mijqidot qjdot - ½ ∑i,j Kij qi qj = ½ qdotT M qdot - ½ qT K q.

The Equation of Motion

We're now ready to look at the equations of motion. There are n of them, corresponding to the n generalized coordinates, each of the form

(d/dt)(∂L/∂qkdot) = ∂L/∂qk,

for k = 1 to n. The double sums present a small complication to the derivatives. Here's a way to get the correct result. First, note that the matrices M and K are symmetric, that is, Mij = Mji and Kij = Kji. Second, notice that there is a term for i=k and one for j=k, so that

L/∂qk = - ½ ∑j Kkj qj - ½ ∑i Kik qi = -∑i Kik qi,

where I've used the fact that i and j are just dummy indices for the summation, and changed j in the first sum to i, then used the symmetry property of K to combine the two terms. Likewise, the other derivative yields

L/∂qkdot = ½ ∑j Mkj qjdot + ½ ∑i Mik qidot = ∑i Mik qidot.
With this simplification the kth equation of motion reads
i Mik qiddot = -∑i Kik qi.

In matrix format this looks like

Mqddot = -Kq

As we did with the pair of oscillators, try a solution of the form q = aeiωt, yielding the eigenvalue equation

(K - ω²M)a = 0

which has a non-trivial solution only if the characteristic equation is satisfied

det(K - ω²M) = 0.

This is an nth order equation in ω². Aside from a few special cases, it must generally be solved numerically. This type of problem is rather common in science and engineering, so there are numerous computer programs available for solving such systems.

© 2008 Robert Harr