1Lagrangian Mechanics

The purpose of mechanics is to describe how bodies change their position in space with “time.” I should load my conscience with grave sins against the sacred spirit of lucidity were I to formulate the aims of mechanics in this way, without serious reflection and detailed explanations. Let us proceed to disclose these sins.

Albert Einstein, Relativity, the Special and General Theory [16], p. 9

The subject of this book is motion and the mathematical tools used to describe it.

Centuries of careful observations of the motions of the planets revealed regularities in those motions, allowing accurate predictions of phenomena such as eclipses and conjunctions. The effort to formulate these regularities and ultimately to understand them led to the development of mathematics and to the discovery that mathematics could be effectively used to describe aspects of the physical world. That mathematics can be used to describe natural phenomena is a remarkable fact.

A pin thrown by a juggler takes a rather predictable path and rotates in a rather predictable way. In fact, the skill of juggling depends crucially on this predictability. It is also a remarkable discovery that the same mathematical tools used to describe the motions of the planets can be used to describe the motion of the juggling pin.

Classical mechanics describes the motion of a system of particles, subject to forces describing their interactions. Complex physical objects, such as juggling pins, can be modeled as myriad particles with fixed spatial relationships maintained by stiff forces of interaction.

There are many conceivable ways a system could move that never occur. We can imagine that the juggling pin might pause in midair or go fourteen times around the head of the juggler before being caught, but these motions do not happen. How can we distinguish motions of a system that can actually occur from other conceivable motions? Perhaps we can invent some mathematical function that allows us to distinguish realizable motions from among all conceivable motions.

The motion of a system can be described by giving the position of every piece of the system at each moment. Such a description of the motion of the system is called a configuration path; the configuration path specifies the configuration as a function of time. The juggling pin rotates as it flies through the air; the configuration of the juggling pin is specified by giving the position and orientation of the pin. The motion of the juggling pin is specified by giving the position and orientation of the pin as a function of time.

The path-distinguishing function that we seek takes a configuration path as an input and produces some output. We want this function to have some characteristic behavior when its input is a realizable path. For example, the output could be a number, and we could try to arrange that this number be zero only on realizable paths. Newton's equations of motion are of this form; at each moment Newton's differential equations must be satisfied.

However, there is an alternate strategy that provides more insight and power: we could look for a path-distinguishing function that has a minimum on the realizable paths—on nearby unrealizable paths the value of the function is higher than it is on the realizable path. This is the variational strategy: for each physical system we invent a path-distinguishing function that distinguishes realizable motions of the system by having a stationary point for each realizable path.¹ For a great variety of systems realizable motions of the system can be formulated in terms of a variational principle.²

Mechanics, as invented by Newton and others of his era, describes the motion of a system in terms of the positions, velocities, and accelerations of each of the particles in the system. In contrast to the Newtonian formulation of mechanics, the variational formulation of mechanics describes the motion of a system in terms of aggregate quantities that are associated with the motion of the system as a whole.

In the Newtonian formulation the forces can often be written as derivatives of the potential energy of the system. The motion of the system is determined by considering how the individual component particles respond to these forces. The Newtonian formulation of the equations of motion is intrinsically a particle-by-particle description.

In the variational formulation the equations of motion are formulated in terms of the difference of the kinetic energy and the potential energy. The potential energy is a number that is characteristic of the arrangement of the particles in the system; the kinetic energy is a number that is determined by the velocities of the particles in the system. Neither the potential energy nor the kinetic energy depends on how those positions and velocities are specified. The difference is characteristic of the system as a whole and does not depend on the details of how the system is specified. So we are free to choose ways of describing the system that are easy to work with; we are liberated from the particle-by-particle description inherent in the Newtonian formulation.

The variational formulation has numerous advantages over the Newtonian formulation. The equations of motion for those parameters that describe the state of the system are derived in the same way regardless of the choice of those parameters: the method of formulation does not depend on the choice of coordinate system. If there are positional constraints among the particles of a system the Newtonian formulation requires that we consider the forces maintaining these constraints, whereas in the variational formulation the constraints can be built into the coordinates. The variational formulation reveals the association of conservation laws with symmetries. The variational formulation provides a framework for placing any particular motion of a system in the context of all possible motions of the system. We pursue the variational formulation because of these advantages.

1.1   Configuration Spaces

Let us consider mechanical systems that can be thought of as composed of constituent point particles, with mass and position, but with no internal structure.³ Extended bodies may be thought of as composed of a large number of these constituent particles with specific spatial relationships among them. Extended bodies maintain their shape because of spatial constraints among the constituent particles. Specifying the position of all the constituent particles of a system specifies the configuration of the system. The existence of constraints among parts of the system, such as those that determine the shape of an extended body, means that the constituent particles cannot assume all possible positions. The set of all configurations of the system that can be assumed is called the configuration space of the system. The dimension of the configuration space is the smallest number of parameters that have to be given to completely specify a configuration. The dimension of the configuration space is also called the number of degrees of freedom of the system.⁴

For a single unconstrained particle it takes three parameters to specify the configuration; a point particle has a three-dimensional configuration space. If we are dealing with a system with more than one point particle, the configuration space is more complicated. If there are k separate particles we need 3k parameters to describe the possible configurations. If there are constraints among the parts of a system the configuration is restricted to a lower-dimensional space. For example, a system consisting of two point particles constrained to move in three dimensions so that the distance between the particles remains fixed has a five-dimensional configuration space: thus with three numbers we can fix the position of one particle, and with two others we can give the position of the other particle relative to the first.

Consider a juggling pin. The configuration of the pin is specified if we give the positions of the atoms making up the pin. However, there exist more economical descriptions of the configuration. In the idealization that the juggling pin is truly rigid, the distances among all the atoms of the pin remain constant. So we can specify the configuration of the pin by giving the position of a single atom and the orientation of the pin. Using the constraints, the positions of all the other constituents of the pin can be determined from this information. The dimension of the configuration space of the juggling pin is six: the minimum number of parameters that specify the position in space is three, and the minimum number of parameters that specify an orientation is also three.

As a system evolves with time, the constituent particles move subject to the constraints. The motion of each constituent particle is specified by describing the changing configuration. Thus, the motion of the system may be described as evolving along a path in configuration space. The configuration path may be specified by a function, the configuration-path function, which gives the configuration of the system at any time.

Exercise 1.1: Degrees of freedom

For each of the mechanical systems described below, give the number of degrees of freedom of the configuration space.

a. Three juggling pins.

b. A spherical pendulum, consisting of a point mass (the pendulum bob) hanging from a rigid massless rod attached to a fixed support point. The pendulum bob may move in any direction subject to the constraint imposed by the rigid rod. The point mass is subject to the uniform force of gravity.

c. A spherical double pendulum, consisting of one point mass hanging from a rigid massless rod attached to a second point mass hanging from a second massless rod attached to a fixed support point. The point masses are subject to the uniform force of gravity.

d. A point mass sliding without friction on a rigid curved wire.

e. A top consisting of a rigid axisymmetric body with one point on the symmetry axis of the body attached to a fixed support, subject to a uniform gravitational force.

f. The same as e, but not axisymmetric.

1.2   Generalized Coordinates

In order to be able to talk about specific configurations we need to have a set of parameters that label the configurations. The parameters used to specify the configuration of the system are called the generalized coordinates. Consider an unconstrained free particle. The configuration of the particle is specified by giving its position. This requires three parameters. The unconstrained particle has three degrees of freedom. One way to specify the position of a particle is to specify its rectangular coordinates relative to some chosen coordinate axes. The rectangular components of the position are generalized coordinates for an unconstrained particle. Or consider an ideal planar double pendulum: a point mass constrained to be a given distance from a fixed point by a rigid rod, with a second mass constrained to be at a given distance from the first mass by another rigid rod, all confined to a vertical plane. The configuration is specified if the orientation of the two rods is given. This requires at least two parameters; the planar double pendulum has two degrees of freedom. One way to specify the orientation of each rod is to specify the angle it makes with a vertical plumb line. These two angles are generalized coordinates for the planar double pendulum.

The number of coordinates need not be the same as the dimension of the configuration space, though there must be at least that many. We may choose to work with more parameters than necessary, but then the parameters will be subject to constraints that restrict the system to possible configurations, that is, to elements of the configuration space.

For the planar double pendulum described above, the two angle coordinates are enough to specify the configuration. We could also take as generalized coordinates the rectangular coordinates of each of the masses in the plane, relative to some chosen coordinate axes. These are also fine coordinates, but we would have to explicitly keep in mind the constraints that limit the possible configurations to the actual geometry of the system. Sets of coordinates with the same dimension as the configuration space are easier to work with because we do not have to deal with explicit constraints among the coordinates. So for the time being we will consider only formulations where the number of configuration coordinates is equal to the number of degrees of freedom; later we will learn how to handle systems with redundant coordinates and explicit constraints.

In general, the configurations form a space M of some dimension n. The n-dimensional configuration space can be parameterized by choosing a coordinate function χ that maps elements of the configuration space to n-tuples of real numbers.⁵ If there is more than one dimension, the function χ is a tuple of n independent coordinate functions⁶ χⁱ, i = 0, …, n − 1, where each χⁱ is a real-valued function defined on some region of the configuration space.⁷ For a given configuration m in the configuration space M the values χⁱ(m) of the coordinate functions are the generalized coordinates of the configuration. These generalized coordinates permit us to identify points of the n-dimensional configuration space with n-tuples of real numbers.⁸ For any given configuration space, there are a great variety of ways to choose generalized coordinates. Even for a single point moving without constraints, we can choose rectangular coordinates, polar coordinates, or any other coordinate system that strikes our fancy.

The motion of the system can be described by a configuration path γ mapping time to configuration-space points. Corresponding to the configuration path is a coordinate path q = χ ∘ γ mapping time to tuples of generalized coordinates.⁹ If there is more than one degree of freedom the coordinate path is a structured object: q is a tuple of component coordinate path functions qⁱ = χⁱ ∘ γ. At each instant of time t, the values q(t) = (q⁰(t), …, qⁿ⁻¹(t)) are the generalized coordinates of a configuration.

The derivative Dq of the coordinate path q is a function¹⁰ that gives the rate of change of the configuration coordinates at a given time: Dq(t) = (Dq⁰(t), …, Dqⁿ⁻¹(t)). The rate of change of a generalized coordinate is called a generalized velocity.

Exercise 1.2: Generalized coordinates

For each of the systems in exercise 1.1, specify a system of generalized coordinates that can be used to describe the behavior of the system.

1.3   The Principle of Stationary Action

Let us suppose that for each physical system there is a path-distinguishing function that is stationary on realizable paths. We will try to deduce some of its properties.

Experience of motion

Our ordinary experience suggests that physical motion can be described by configuration paths that are continuous and smooth.¹¹ We do not see the juggling pin jump from one place to another. Nor do we see the juggling pin suddenly change the way it is moving.

Our ordinary experience suggests that the motion of physical systems does not depend upon the entire history of the system. If we enter the room after the juggling pin has been thrown into the air we cannot tell when it left the juggler's hand. The juggler could have thrown the pin from a variety of places at a variety of times with the same apparent result as we walk through the door.¹² So the motion of the pin does not depend on the details of the history.

Our ordinary experience suggests that the motion of physical systems is deterministic. In fact, a small number of parameters summarize the important aspects of the history of the system and determine its future evolution. For example, at any moment the position, velocity, orientation, and rate of change of the orientation of the juggling pin are enough to completely determine the future motion of the pin.

Realizable paths

From our experience of motion we develop certain expectations about realizable configuration paths. If a path is realizable, then any segment of the path is a realizable path segment. Conversely, a path is realizable if every segment of the path is a realizable path segment. The realizability of a path segment depends on all points of the path in the segment. The realizability of a path segment depends on every point of the path segment in the same way; no part of the path is special. The realizability of a path segment depends only on points of the path within the segment; the realizability of a path segment is a local property.

So the path-distinguishing function aggregates some local property of the system measured at each moment along the path segment. Each moment along the path must be treated in the same way. The contributions from each moment along the path segment must be combined in a way that maintains the independence of the contributions from disjoint subsegments. One method of combination that satisfies these requirements is to add up the contributions, making the path-distinguishing function an integral over the path segment of some local property of the path.¹³

So we will try to arrange that the path-distinguishing function, constructed as an integral of a local property along the path, assumes a stationary value for any realizable path. Such a path-distinguishing function is traditionally called an action for the system. We use the word “action” to be consistent with common usage. Perhaps it would be clearer to continue to call it “path-distinguishing function,” but then it would be more difficult for others to know what we were talking about.¹⁴

In order to pursue the agenda of variational mechanics, we must invent action functions that are stationary on the realizable trajectories of the systems we are studying. We will consider actions that are integrals of some local property of the configuration path at each moment. Let q = χ ∘ γ be a coordinate path in the configuration space; q(t) are the coordinates of the configuration at time t. Then the action of a segment of the path in the time interval from t₁ to t₂ is¹⁵

$$\begin{array}{ll}S[q](t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}F[q].}\hfill & (1.1)\hfill \end{array}$$

where F[q] is a function of time that measures some local property of the path. It may depend upon the value of the function q at that time and the value of any derivatives of q at that time.¹⁶

The configuration path can be locally described at a moment in terms of the coordinates, the rate of change of the coordinates, and all the higher derivatives of the coordinates at the given moment. Given this information the path can be reconstructed in some interval containing that moment.¹⁷ Local properties of paths can depend on no more than the local description of the path.

The function F measures some local property of the coordinate path q. We can decompose F [q] into two parts: a part that measures some property of a local description and a part that extracts a local description of the path from the path function. The function that measures the local property of the system depends on the particular physical system; the method of construction of a local description of a path from a path is the same for any system. We can write F[q] as a composition of these two functions:¹⁸

$$\begin{array}{ll}F[q]=L\circ \Gamma [q].\hfill & (1.2)\hfill \end{array}$$

The function Γ takes the coordinate path and produces a function of time whose value is an ordered tuple containing the time, the coordinates at that time, the rate of change of the coordinates at that time, and the values of higher derivatives of the coordinates evaluated at that time. For the path q and time t:

$$\begin{array}{ll}\Gamma [q](t)=(t,q(t),Dq(t),\dots ).\hfill & (1.3)\hfill \end{array}$$

We refer to this tuple, which includes as many derivatives as are needed, as the local tuple. The function Γ[q] depends only on the coordinate path q and its derivatives; the function Γ[q] does not depend on χ or the fact that q is made by composing χ with γ.

The function L depends on the specific details of the physical system being investigated, but does not depend on any particular configuration path. The function L computes a real-valued local property of the path. We will find that L needs only a finite number of components of the local tuple to compute this property: The path can be locally reconstructed from the full local description; that L depends on a finite number of components of the local tuple guarantees that it measures a local property.¹⁹

The advantage of this decomposition is that the local description of the path is computed by a uniform process from the configuration path, independent of the system being considered. All of the system-specific information is captured in the function L.

The function L is called a Lagrangian²⁰ for the system, and the resulting action,

$$\begin{array}{ll}S[q](t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}L\circ \Gamma [q],}\hfill & (1.4)\hfill \end{array}$$

is called the Lagrangian action. For Lagrangians that depend only on time, positions, and velocities the action can also be written

$$\begin{array}{ll}S[q](t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}L(t,q(t),Dq(t))dt.}\hfill & (1.5)\hfill \end{array}$$

Lagrangians can be found for a great variety of systems. We will see that for many systems the Lagrangian can be taken to be the difference between kinetic and potential energy. Such Lagrangians depend only on the time, the configuration, and the rate of change of the configuration. We will focus on this class of systems, but will also consider more general systems from time to time.

A realizable path of the system is to be distinguished from others by having stationary action with respect to some set of nearby unrealizable paths. Now some paths near realizable paths will also be realizable: for any motion of the juggling pin there is another that is slightly different. So when addressing the question of whether the action is stationary with respect to variations of the path we must somehow restrict the set of paths we are considering to contain only one realizable path. It will turn out that for Lagrangians that depend only on the configuration and rate of change of configuration it is enough to restrict the set of paths to those that have the same configuration at the endpoints of the path segment.

The principle of stationary action asserts that for each dynamical system we can cook up a Lagrangian such that a realizable path connecting the configurations at two times t₁ and t₂ is distinguished from all conceivable paths by the fact that the action S[q](t₁, t₂) is stationary with respect to variations of the path.²¹ For Lagrangians that depend only on the configuration and rate of change of configuration, the variations are restricted to those that preserve the configurations at t₁ and t₂.²²

Exercise 1.3: Fermat optics

Fermat observed that the laws of reflection and refraction could be accounted for by the following facts: Light travels in a straight line in any particular medium with a velocity that depends upon the medium. The path taken by a ray from a source to a destination through any sequence of media is a path of least total time, compared to neighboring paths. Show that these facts imply the laws of reflection and refraction.²³

1.4   Computing Actions

To illustrate the above ideas, and to introduce their formulation as computer programs, we consider the simplest mechanical system—a free particle moving in three dimensions. Euler and Lagrange discovered that for a free particle the time integral of the kinetic energy over the particle's actual path is smaller than the same integral along any alternative path between the same points: a free particle moves according to the principle of stationary action, provided we take the Lagrangian to be the kinetic energy. The kinetic energy for a particle of mass m and velocity $\overrightarrow{v}$ is $\frac{1}{2}m{v}^2$, where v is the magnitude of $\overrightarrow{v}$. In this case we can choose the generalized coordinates to be the ordinary rectangular coordinates.

Following Euler and Lagrange, the Lagrangian for the free particle is²⁴

$$\begin{array}{ll}L(t,x,v)=\frac{1}{2}m(v·v),\hfill & (1.6)\hfill \end{array}$$

where the formal parameter x names a tuple of components of the position with respect to a given rectangular coordinate system, and the formal parameter v names a tuple of velocity components.²⁵

We can express this formula as a procedure:

(define ((L-free-particle mass) local)
  (let ((v (velocity local)))
    (* 1/2 mass (dot-product v v))))

The definition indicates that L-free-particle is a procedure that takes mass as an argument and returns a procedure that takes a local tuple local, extracts the generalized velocity with the procedure velocity, and uses the velocity to compute the value of the Lagrangian.²⁶

Suppose we let q denote a coordinate path function that maps time to position components:²⁷

$$\begin{array}{ll}q(t)=(x(t),y(t),z(t)).\hfill & (1.7)\hfill \end{array}$$

We can make this definition²⁸

(define q
  (up (literal-function 'x)
      (literal-function 'y)
      (literal-function 'z)))

where literal-function makes a procedure that represents a function of one argument that has no known properties other than the given symbolic name. The symbol q now names a procedure of one real argument (time) that produces a tuple of three components representing the coordinates at that time. For example, we can evaluate this procedure for a symbolic time t as follows:

(q 't)

(up (x t) (y t) (z t))

The derivative of the coordinate path Dq is the function that maps time to velocity components:

$Dq(t)=(Dx(t),Dy(t),Dz(t)).$

We can make and use the derivative of a function.²⁹ For example, we can write:

((D q) 't)

(up ((D x) t) ((D y) t) ((D z) t))

The function Γ takes a coordinate path and returns a function of time that gives the local tuple (t, q(t), Dq(t), …). We implement this Γ with the procedure Gamma.³⁰ Here is what Gamma does:

((Gamma q) 't)

(up t
    (up (x t) (y t) (z t))
    (up ((D x) t) ((D y) t) ((D z) t)))

So the composition L ∘ Γ is a function of time that returns the value of the Lagrangian for this point on the path:³¹

((compose (L-free-particle 'm) (Gamma q)) 't)

(+ (* 1/2 m (expt ((D x) t) 2))
    (* 1/2 m (expt ((D y) t) 2))
    (* 1/2 m (expt ((D z) t) 2)))

The procedure show-expression simplifies the expression and uses T_EX to display the result in traditional infix form. We use this method of display to make the boxed expressions in this book. The procedure show-expression also produces the prefix form, but we usually do not show this.³²

(show-expression
  ((compose (L-free-particle 'm) (Gamma q)) 't))

According to equation (1.4) we can compute the Lagrangian action from time t₁ to time t₂ as:

(define (Lagrangian-action L q t1 t2)
  (definite-integral (compose L (Gamma q)) t1 t2))

Lagrangian-action takes as arguments a procedure L that computes the Lagrangian, a procedure q that computes a coordinate path, and starting and ending times t1 and t2. The definite-integral used here takes as arguments a function and two limits t1 and t2, and computes the definite integral of the function over the interval from t1 to t2.³³ Notice that the definition of Lagrangian-action does not depend on any particular set of coordinates or even the dimension of the configuration space. The method of computing the action from the coordinate representation of a Lagrangian and a coordinate path does not depend on the coordinate system.

We can now compute the action for the free particle along a path. For example, consider a particle moving at uniform speed along a straight line t ↦ (4t + 7, 3t + 5, 2t + 1).³⁴ We represent the path as a procedure

(define (test-path t)
  (up (+ (* 4 t) 7)
      (+ (* 3 t) 5)
      (+ (* 2 t) 1)))

For a particle of mass 3, we obtain the action between t = 0 and t = 10 as³⁵

(Lagrangian-action (L-free-particle 3.0)
                   test-path 0.0 10.0)

435

Exercise 1.4: Lagrangian actions

For a free particle an appropriate Lagrangian is³⁶

$$\begin{array}{ll}L(t,x,v)=\frac{1}{2}m{v}^2.\hfill & (1.8)\hfill \end{array}$$

Suppose that x is the constant-velocity straight-line path of a free particle, such that x_a = x(t_a) and x_b = x(t_b). Show that the action on the solution path is

$$\begin{array}{ll}\frac{m}{2}\frac{{(x_b-x_a)}^2}{t_b-t_a}.\hfill & (1.9)\hfill \end{array}$$

Paths of minimum action

We already know that the actual path of a free particle is uniform motion in a straight line. According to Euler and Lagrange, the action is smaller along a straight-line test path than along nearby paths. Let q be a straight-line test path with action S[q](t₁, t₂). Let q + ϵη be a nearby path, obtained from q by adding a path variation η scaled by the real parameter ϵ.³⁷ The action on the varied path is S[q + ϵη](t₁, t₂). Euler and Lagrange found that S[q + ϵη](t₁, t₂) > S[q](t₁, t₂) for any η that is zero at the endpoints and for any small nonzero ϵ.

Let's check this numerically by varying the test path, adding some amount of a test function that is zero at the endpoints t = t₁ and t = t₂. To make a function η that is zero at the endpoints, given a sufficiently well-behaved function ν, we can use η(t) = (t − t₁)(t − t₂)ν(t). This can be implemented:

(define ((make-eta nu t1 t2) t)
  (* (- t t1) (- t t2) (nu t)))

We can use this to compute the action for a free particle over a path varied from the given path, as a function of ϵ:³⁸

(define ((varied-free-particle-action mass q nu t1 t2) eps)
  (let ((eta (make-eta nu t1 t2)))
    (Lagrangian-action (L-free-particle mass)
                       (+ q (* eps eta))
                       t1
                       t2)))

The action for the varied path, with ν(t) = (sin t, cos t, t²) and ϵ = 0.001, is, as expected, larger than for the test path:

((varied-free-particle-action 3.0 test-path
                              (up sin cos square)
                              0.0 10.0)
 0.001)

436.29121428571153

We can numerically compute the value of ϵ for which the action is minimized. We search between, say, −2 and 1:³⁹

(minimize
  (varied-free-particle-action 3.0 test-path
                               (up sin cos square)
                               0.0 10.0)
  -2.0 1.0)

(-1.5987211554602254e-14 435.0000000000237 5)

We find exactly what is expected—that the best value for ϵ is zero,⁴⁰ and the minimum value of the action is the action along the straight path.

Finding trajectories that minimize the action

We have used the variational principle to determine if a given trajectory is realizable. We can also use the variational principle to find trajectories. Given a set of trajectories that are specified by a finite number of parameters, we can search the parameter space looking for the trajectory in the set that best approximates the real trajectory by finding one that minimizes the action. By choosing a good set of approximating functions we can get arbitrarily close to the real trajectory.⁴¹

One way to make a parametric path that has fixed endpoints is to use a polynomial that goes through the endpoints as well as a number of intermediate points. Variation of the positions of the intermediate points varies the path; the parameters of the varied path are the coordinates of the intermediate positions. The procedure make-path constructs such a path using a Lagrange interpolation polynomial. The procedure make-path is called with five arguments: (make-path t0 q0 t1 q1 qs), where q0 and q1 are the endpoints, t0 and t1 are the corresponding times, and qs is a list of intermediate points.⁴²

Having specified a parametric path, we can construct a parametric action that is just the action computed along the parametric path:

(define ((parametric-path-action Lagrangian t0 q0 t1 q1) qs)
  (let ((path (make-path t0 q0 t1 q1 qs)))
    (Lagrangian-action Lagrangian path t0 t1)))

We can find approximate solution paths by finding parameters that minimize the action. We do this minimization with a canned multidimensional minimization procedure:⁴³

(define (find-path Lagrangian t0 q0 t1 q1 n)
  (let ((initial-qs (linear-interpolants q0 q1 n)))
    (let ((minimizing-qs
            (multidimensional-minimize
              (parametric-path-action Lagrangian t0 q0 t1 q1)
              initial-qs)))
      (make-path t0 q0 t1 q1 minimizing-qs))))

The procedure multidimensional-minimize takes a procedure (in this case the value of the call to parametric-path-action) that computes the function to be minimized (in this case the action) and an initial guess for the parameters. Here we choose the initial guess to be equally spaced points on a straight line between the two endpoints, computed with linear-interpolants.

To illustrate the use of this strategy, we will find trajectories of the harmonic oscillator, with Lagrangian⁴⁴

$$\begin{array}{ll}L(t,q,v)=\frac{1}{2}m{v}^2-\frac{1}{2}k{q}^2,\hfill & (1.10)\hfill \end{array}$$

for mass m and spring constant k. This Lagrangian is implemented by⁴⁵

(define ((L-harmonic m k) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (- (* 1/2 m (square v)) (* 1/2 k (square q)))))

We can find an approximate path taken by the harmonic oscillator for m = 1 and k = 1 between q(0) = 1 and q(π/2) = 0 as follows:⁴⁶

(define q
  (find-path (L-harmonic 1.0 1.0) 0.0 1.0 :pi/2 0.0 3))

We know that the trajectories of this harmonic oscillator, for m = 1 and k = 1, are

$$\begin{array}{ll}q(t)=A\text{cos}(t+\phi )\hfill & (1.11)\hfill \end{array}$$

where the amplitude A and the phase φ are determined by the initial conditions. For the chosen endpoint conditions the solution is q(t) = cos(t). The approximate path should be an approximation to cosine over the range from 0 to π/2. Figure 1.1 shows the error in the polynomial approximation produced by this process. The maximum error in the approximation with three intermediate points is less than 1.7 × 10⁻⁴. We find, as expected, that the error in the approximation decreases as the number of intermediate points is increased. For four intermediate points it is about a factor of 15 better.

Exercise 1.5: Solution process

We can watch the progress of the minimization by modifying the procedure parametric-path-action to plot the path each time the action is computed. Try this:

(define win2 (frame 0.0 :pi/2 0.0 1.2))

(define ((parametric-path-action Lagrangian t0 q0 t1 q1)
         intermediate-qs)
  (let ((path (make-path t0 q0 t1 q1 intermediate-qs)))
    ;; display path
    (graphics-clear win2)
    (plot-function win2 path t0 t1 (/ (- t1 t0) 100))
    ;; compute action
    (Lagrangian-action Lagrangian path t0 t1)))

(find-path (L-harmonic 1.0 1.0) 0.0 1.0 :pi/2 0.0 2)

Exercise 1.6: Minimizing action

Suppose we try to obtain a path by minimizing an action for an impossible problem. For example, suppose we have a free particle and we impose endpoint conditions on the velocities as well as the positions that are inconsistent with the particle being free. Does the formalism protect itself from such an unpleasant attack? You may find it illuminating to program it and see what happens.

1.5   The Euler–Lagrange Equations

The principle of stationary action characterizes the realizable paths of systems in configuration space as those for which the action has a stationary value. In elementary calculus, we learn that the critical points of a function are the points where the derivative vanishes. In an analogous way, the paths along which the action is stationary are solutions of a system of differential equations. This system, called the Euler–Lagrange equations or just the Lagrange equations, is the link that permits us to use the principle of stationary action to compute the motions of mechanical systems, and to relate the variational and Newtonian formulations of mechanics.

Lagrange equations

We will find that if L is a Lagrangian for a system that depends on time, coordinates, and velocities, and if q is a coordinate path for which the action S[q](t₁, t₂) is stationary (with respect to any variation in the path that keeps the endpoints of the path fixed), then

$$\begin{array}{ll}D({\partial }_{2}L\circ \Gamma [q])-{\partial }_{1}L\circ \Gamma [q]=0.\hfill & (1.12)\hfill \end{array}$$

Here L is a real-valued function of a local tuple; ∂₁L and ∂₂L denote the partial derivatives of L with respect to its generalized position argument and generalized velocity argument respectively.⁴⁷ The function ∂₂L maps a local tuple to a structure whose components are the derivatives of L with respect to each component of the generalized velocity. The function Γ[q] maps time to the local tuple: Γ[q](t) = (t, q(t), Dq(t), …). Thus the compositions ∂₁L ∘ Γ[q] and ∂₂L ∘ Γ[q] are functions of one argument, time. The Lagrange equations assert that the derivative of ∂₂L ∘ Γ[q] is equal to ∂₁L ∘ Γ[q], at any time. Given a Lagrangian, the Lagrange equations form a system of ordinary differential equations that must be satisfied by realizable paths.

Lagrange's equations are traditionally written as a separate equation for each component of q:

$\begin{array}{ll}\frac{d}{dt}\frac{\partial L}{\partial {\dot{q}}^i}-\frac{\partial L}{\partial q^i}=0\hfill & i=0,\dots ,n-1.\hfill \end{array}$

In this way of writing Lagrange's equations the notation does not distinguish between L, which is a real-valued function of three variables (t, q, $$ \dot{q}$$), and L ∘ Γ[q], which is a real-valued function of one real variable t. If we do not realize this notational pun, the equations don't make sense as written—∂L/∂$$ \dot{q}$$ is a function of three variables, so we must regard the arguments q, $$ \dot{q}$$ as functions of t before taking d/dt of the expression. Similarly, ∂L/∂q is a function of three variables, which we must view as a function of t before setting it equal to $d/dt\left(\partial L/\partial \dot{q}\right)$.

A correct use of the traditional notation is more explicit:

$\frac{d}{dt}\left({\frac{\partial L\left(t,w,\dot{w}\right)}{\partial {\dot{w}}^i}|}_{\begin{array}{l}w=q\left(t\right)\hfill \\ \dot{w}=\frac{dq\left(t\right)}{dt}\hfill \end{array}}\right)-{\frac{\partial L\left(t,w,\dot{w}\right)}{\partial w^i}|}_{\begin{array}{l}w=q\left(t\right)\hfill \\ \dot{w}=\frac{dq\left(t\right)}{dt}\hfill \end{array}}=0,$

where i = 0, …, n − 1. In these equations we see that the partial derivatives of the Lagrangian function are taken, then the path and its derivative are substituted for the position and velocity arguments of the Lagrangian, resulting in an expression in terms of the time.

1.5.1 Derivation of the Lagrange Equations

We will show that the principle of stationary action implies that realizable paths satisfy the Euler–Lagrange equations.

A Direct Derivation

Let q be a realizable coordinate path from (t₁, q(t₁)) to (t₂, q(t₂)). Consider nearby paths q + ϵη where η(t₁) = η(t₂) = 0. Let

$$\begin{array}{lll}g(\mathit{ϵ})\hfill & =S[q+\mathit{ϵ}\eta ](t_1,t_2)\hfill & \hfill \\ \hfill & ={\displaystyle {\int }_{t_1}^{t_2}L(t,q(t)+\mathit{ϵ}\eta (t),Dq(t)+\mathit{ϵ}D\eta (t))dt.}\hfill & (1.13)\hfill \end{array}$$

Expanding as a power series in ϵ

$$\begin{array}{ll}g(\mathit{ϵ})=g(0)+\mathit{ϵ}Dg(0)+\cdots \hfill & (1.14)\hfill \end{array}$$

and using the chain rule we get

$$\begin{array}{lll}Dg(0)\hfill & ={\displaystyle {\int }_{t_1}^{t_2}({\partial }_{1}L(t,q(t),Dq(t))\eta (t))dt}\hfill & \hfill \\ \hfill & +{\displaystyle {\int }_{t_1}^{t_2}({\partial }_{2}L(t,q(t),Dq(t))D\eta (t))dt.}\hfill & (1.15)\hfill \end{array}$$

Integrating the second term by parts we obtain

$$\begin{array}{lll}Dg(0)\hfill & ={\displaystyle {\int }_{t_1}^{t_2}({\partial }_{1}L(t,q(t),Dq(t))\eta (t))dt}\hfill & \hfill \\ \hfill & +{\partial }_{2}L(t,q(t),Dq(t)D\eta (t))\eta (t){|}_{t_1}^{t_2}\hfill & \hfill \\ \hfill & -{\displaystyle {\int }_{t_1}^{t_2}\frac{d}{dt}({\partial }_{2}L(t,q(t),Dq(t)))\eta (t)dt}.\hfill & (1.16)\hfill \end{array}$$

The increment ΔS in the action due to the variation in the path is, to first order in ϵ, ϵDg(0). Because η is zero at the endpoints the integrated term is zero. Collecting together the other two terms, and reverting to functional notation, we find the increment to be

$$\begin{array}{ll}\Delta S=\mathit{ϵ}{\displaystyle {\int }_{t_1}^{t_2}\{{\partial }_{1}L\circ \Gamma [q]-D({\partial }_{2}L\circ \Gamma [q])\}}\eta .\hfill & (1.17)\hfill \end{array}$$

If ΔS is zero the action is stationary. We retain enough freedom in the choice of the variation that the factor in the integrand multiplying η is forced to be zero at each point along the path. We argue by contradiction: Suppose this factor were nonzero at some particular time. Then it would have to be nonzero in at least one of its components. But if we choose our η to be a bump that is nonzero only in that component in a neighborhood of that time, and zero everywhere else, then the integral will be nonzero. So we may conclude that the factor in curly brackets is identically zero and thus obtain Lagrange's equations:⁴⁸

$$\begin{array}{ll}D({\partial }_{2}L\circ \Gamma [q])-{\partial }_{1}L\circ \Gamma [q]=0.\hfill & (1.18)\hfill \end{array}$$

The Variation Operator

First we will develop tools for investigating how path-dependent functions vary as the paths are varied. We will then apply these tools to the action, to derive the Lagrange equations.

Suppose that we have a function f[q] that depends on a path q. How does the function vary as the path is varied? Let q be a coordinate path and q + ϵη be a varied path, where the function η is a path-like function that can be added to the path q, and the factor ϵ is a scale factor. We define the variation δ_ηf[q] of the function f on the path q by⁴⁹

$$\begin{array}{ll}{\delta }_{\eta }f[q]=\underset{\mathit{ϵ}\to 0}{\lim }\left(\frac{f[q+\mathit{ϵ}\eta ]-f[q]}{\mathit{ϵ}}\right).\hfill & (1.19)\hfill \end{array}$$

The variation of f is a linear approximation to the change in the function f for small variations in the path. The variation of f depends on η.

A simple example is the variation of the identity path function: I[q] = q. Applying the definition, we find

$$\begin{array}{ll}{\delta }_{\eta }I[q]=\underset{\mathit{ϵ}\to 0}{\lim }\left(\frac{(q+\mathit{ϵ}\eta )-q}{\mathit{ϵ}}\right)=\eta .\hfill & (1.20)\hfill \end{array}$$

It is traditional to write δ_ηI[q] simply as δq. Another example is the variation of the path function that returns the derivative of the path. We have⁵⁰

$$\begin{array}{ll}{\delta }_{\eta }g[q]=\underset{\mathit{ϵ}\to 0}{\lim }\left(\frac{D(q+\mathit{ϵ}\eta )-Dq}{\mathit{ϵ}}\right)=D\eta \text{with}g[q]=Dq.\hfill & (1.21)\hfill \end{array}$$

It is traditional to write δ_ηg[q] as δDq.

The variation may be represented in terms of a derivative. Let g(ϵ) = f[q + ϵη]; then

$$\begin{array}{ll}{\delta }_{\eta }f[q]=\underset{\mathit{ϵ}\to 0}{\lim }\left(\frac{g(\mathit{ϵ})-g(0)}{\mathit{ϵ}}\right)=Dg(0).\hfill & (1.22)\hfill \end{array}$$

Variations have the following derivative-like properties. For path-dependent functions f and g and constant c:

$$\begin{array}{ll}{\delta }_{\eta }(fg)[q]={\delta }_{\eta }f[q]g[q]+f[q]{\delta }_{\eta }g[q]\hfill & (1.23)\hfill \end{array}$$

$$\begin{array}{ll}{\delta }_{\eta }(f+g)[q]={\delta }_{\eta }f[q]+{\delta }_{\eta }g[q]\hfill & (1.24)\hfill \end{array}$$

$$\begin{array}{ll}{\delta }_{\eta }(cf)[q]=c{\delta }_{\eta }f[q].\hfill & (1.25)\hfill \end{array}$$

Let F be a path-independent function and g be a path-dependent function; then

$$\begin{array}{ll}{\delta }_{\eta }h[q]=(DF\circ g[q]){\delta }_{\eta }g[q]\text{with}h[q]=F\circ g[q].\hfill & (1.26)\hfill \end{array}$$

The operators D (differentiation) and δ (variation) commute in the following sense:

$$D\begin{array}{ll}{\delta }_{\eta }f[q]={\delta }_{\eta }g[q]\text{with}g[q]=D(f[q]).\hfill & (1.27)\hfill \end{array}$$

Variations also commute with integration in a similar sense.

If a path-dependent function f is stationary for a particular path q with respect to small changes in that path, then it must be stationary for a subset of those variations that results from adding small multiples of a particular function η to q. So the statement δ_ηf[q] = 0 for arbitrary η implies the function f is stationary for small variations of the path around q.

Exercise 1.7: Properties of δ

Show that δ has the properties 1.23–1.27.

Exercise 1.8: Implementation of δ

a. Suppose we have a procedure f that implements a path-dependent function: for path q and time t it has the value ((f q) t). The procedure delta computes the variation (δ_ηf)[q](t) as the value of the expression ((((delta eta) f) q) t). Complete the definition of delta:

(define (((delta eta) f) q)
    ...
)

b. Use your delta procedure to verify the properties of δ listed in exercise 1.7 for simple functions such as implemented by the procedure f:⁵¹

(define (f q)
  (compose
    (literal-function 'F
                      (-> (UP Real (UP* Real) (UP* Real)) Real))
    (Gamma q)))

This implements an n-degree-of-freedom path-dependent function that depends on the local tuple of the path at each moment. You can define a literal two-dimensional path by

(define q (literal-function 'q (-> Real (UP Real Real))))

You should compute both sides of the equalities and subtract the results. The answer should be zero.

A Derivation with the Variation Operator

The action is the integral of the Lagrangian along a path:

$$\begin{array}{ll}S[q](t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}L\circ \Gamma [q].}\hfill & (1.28)\hfill \end{array}$$

For a realizable path q the variation of the action with respect to any variation η that preserves the endpoints, η(t₁) = η(t₂) = 0, is zero:

$$\begin{array}{ll}{\delta }_{\eta }S[q](t_1,t_2)=0.\hfill & (1.29)\hfill \end{array}$$

Variation commutes with integration, so the variation of the action is

$$\begin{array}{ll}{\delta }_{\eta }S[q](t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}{\delta }_{\eta }h[q]\text{where}h[q]=L\circ \Gamma [q].}\hfill & (1.30)\hfill \end{array}$$

Using the fact that

$$\begin{array}{ll}{\delta }_{\eta }\Gamma [q](t)=(0,\eta (t),D\eta (t)),\hfill & (1.31)\hfill \end{array}$$

which follows from equations (1.20) and (1.21), and using the chain rule for variations (1.26), we get⁵²

$$\begin{array}{lll}{\delta }_{\eta }S[q](t_1,t_2)\hfill & ={\displaystyle {\int }_{t_1}^{t_2}(DL\circ \Gamma [q]){\delta }_{\eta }\Gamma [q]}\hfill & \hfill \\ \hfill & ={\displaystyle {\int }_{t_1}^{t_2}(({\partial }_{1}L\circ \Gamma [q])\eta +({\partial }_{2}L\circ \Gamma [q])D\eta ).}\hfill & (1.32)\hfill \end{array}$$

Integrating the last term of equation (1.32) by parts gives

$$\begin{array}{lll}{\delta }_{\eta }S[q](t_1,t_2)=\hfill & ({\partial }_{2}L\circ \Gamma [q])\eta {|}_{t_1}^{t_2}\hfill & \hfill \\ \hfill & +{\displaystyle {\int }_{t_1}^{t_2}\{({\partial }_{1}L\circ \Gamma [q])-D({\partial }_{2}L\circ \Gamma [q])\}\eta .}\hfill & (1.33)\hfill \end{array}$$

For our variation η we have η(t₁) = η(t₂) = 0, so the first term vanishes.

Thus the variation of the action is zero if and only if

$$\begin{array}{ll}0={\displaystyle {\int }_{t_1}^{t_2}\{({\partial }_{1}L\circ \Gamma [q])-D({\partial }_{2}L\circ \Gamma [q])\}\eta .}\hfill & (1.34)\hfill \end{array}$$

The variation of the action is zero because, by assumption, q is a realizable path. Thus (1.34) must be true for any function η that is zero at the endpoints. Since η is arbitrary, except for being zero at the endpoints, the bracketed factor of the integrand is zero. So

$$\begin{array}{ll}D({\partial }_{2}L\circ \Gamma [q])-({\partial }_{1}L\circ \Gamma [q])=0.\hfill & (1.35)\hfill \end{array}$$

This is just what we set out to obtain, the Lagrange equations.

A path satisfying Lagrange's equations is one for which the action is stationary, and the fact that the action is stationary depends only on the values of L at each point of the path (and at each point on nearby paths), not on the coordinate system we use to compute these values. So if the system's path satisfies Lagrange's equations in some particular coordinate system, it must satisfy Lagrange's equations in any coordinate system. Thus the equations of variational mechanics are derived the same way in any configuration space and any coordinate system.

Harmonic oscillator

For an example, consider the harmonic oscillator. A Lagrangian is

$$\begin{array}{ll}L(t,x,v)=\frac{1}{2}m{v}^2-\frac{1}{2}k{x}^2.\hfill & (1.36)\hfill \end{array}$$

Then

$$\begin{array}{ll}{\partial }_{1}L(t,x,v)=-kx\text{and}{\partial }_{2}L(t,x,v)=mv.\hfill & (1.37)\hfill \end{array}$$

The Lagrangian is applied to a tuple of the time, a coordinate, and a velocity. The symbols t, x, and v are arbitrary; they are used to specify formal parameters of the Lagrangian.

Now suppose we have a configuration path y, which gives the coordinate of the oscillator y(t) for each time t. The initial segment of the corresponding local tuple at time t is

$$\begin{array}{ll}\Gamma [y](t)=(t,y(t),Dy(t)).\hfill & (1.38)\hfill \end{array}$$

So

$$\begin{array}{ll}({\partial }_{1}L\circ \Gamma [y])(t)=-ky(t)\text{and}({\partial }_{2}L\circ \Gamma [y])(t)=mDy(t),\hfill & (1.39)\hfill \end{array}$$

and

$$\begin{array}{ll}D({\partial }_{2}L\circ \Gamma [y])(t)=m{D}^{2}y(t),\hfill & (1.40)\hfill \end{array}$$

so the Lagrange equation is

$$\begin{array}{ll}m{D}^{2}y(t)+ky(t)=0,\hfill & (1.41)\hfill \end{array}$$

which is the equation of motion of the harmonic oscillator.

Orbital motion

As another example, consider the two-dimensional motion of a particle of mass m orbiting a fixed center of attraction, with gravitational potential energy −μ/r, where r is the distance to the center of attraction. This is called the Kepler problem.

A Lagrangian for this problem is⁵³

$$\begin{array}{ll}L(t;\xi ,\eta ;v_{\xi },v_{\eta })=\frac{1}{2}m(v_{\xi }^2+v_{\eta }^2)+\frac{\mu }{\sqrt{{\xi }^2+{\eta }^2}},\hfill & (1.42)\hfill \end{array}$$

where ξ and η are formal parameters for rectangular coordinates of the particle, and v_ξ and v_η are formal parameters for corresponding rectangular velocity components. Then

$$\begin{array}{lll}{\partial }_{1}L(t;\xi ,\eta ;v_{\xi },v_{\eta })\hfill & =[{\partial }_{1,0}L(t;\xi ,\eta ;v_{\xi },v_{\eta }),{\partial }_{1,1}L(t;\xi ,\eta ;v_{\xi },v_{\eta })]\hfill & \hfill \\ \hfill & =\left[\frac{-\mu \xi }{{({\xi }^2+{\eta }^2)}^{3/2}},\frac{-\mu \eta }{{({\xi }^2+{\eta }^2)}^{3/2}}\right].\hfill & (1.43)\hfill \end{array}$$

Similarly,

$$\begin{array}{ll}{\partial }_{2}L(t;\xi ,\eta ;v_{\xi },v_{\eta })=[m{v}_{\xi },m{v}_{\eta }].\hfill & (1.44)\hfill \end{array}$$

Now suppose we have a configuration path q = (x, y), so that the coordinate tuple at time t is q(t) = (x(t), y(t)). The initial segment of the local tuple at time t is

$$\begin{array}{ll}\Gamma [q](t)=(t;x(t),y(t);Dx(t),Dy(t)).\hfill & (1.45)\hfill \end{array}$$

So

$$\begin{array}{lll}({\partial }_{1}L\circ \Gamma [q])(t)\hfill & =\left[\frac{-\mu x(t)}{{({(x(t))}^2+{(y(t))}^2)}^{3/2}},\frac{-\mu y(t)}{{({(x(t))}^2+{(y(t))}^2)}^{3/2}}\right]\hfill & \hfill \\ ({\partial }_{2}L\circ \Gamma [q])(t)\hfill & =[mDx(t),mDy(t)]\hfill & (1.46)\hfill \end{array}$$

and

$$\begin{array}{ll}D({\partial }_{2}L\circ \Gamma [q])(t)=[m{D}^{2}x(t),m{D}^{2}y(t)].\hfill & (1.47)\hfill \end{array}$$

The component Lagrange equations at time t are

$$\begin{array}{ll}m{D}^{2}x(t)+\frac{\mu x(t)}{{({(x(t))}^2+{(y(t))}^2)}^{3/2}}=0\hfill & \hfill \\ m{D}^{2}y(t)+\frac{\mu y(t)}{{({(x(t))}^2+{(y(t))}^2)}^{3/2}}=0.\hfill & (1.48)\hfill \end{array}$$

Exercise 1.9: Lagrange's equations

Derive the Lagrange equations for the following systems, showing all of the intermediate steps as in the harmonic oscillator and orbital motion examples.

a. An ideal planar pendulum consists of a bob of mass m connected to a pivot by a massless rod of length l subject to uniform gravitational acceleration g. A Lagrangian is $L(t,\theta ,\dot{\theta })=\frac{1}{2}m{l}^2{\dot{\theta }}^2+mgl\cos \theta$. The formal parameters of L are t, θ, and $\dot{\theta }$; θ measures the angle of the pendulum rod to a plumb line and $\dot{\theta }$ is the angular velocity of the rod.⁵⁴

b. A particle of mass m moves in a two-dimensional potential V (x, y) = (x² + y²)/2 + x²y − y³/3, where x and y are rectangular coordinates of the particle. A Lagrangian is $L(t;x,y;v_x,v_y)=\frac{1}{2}m\left(v_x^2+v_y^2\right)-V(x,y)$.

c. A Lagrangian for a particle of mass m constrained to move on a sphere of radius R is $L(t;\theta ,\phi ;\alpha ,\beta )=\frac{1}{2}m{R}^2({\alpha }^2+{(\beta \sin \theta )}^2)$. The angle θ is the colatitude of the particle and φ is the longitude; the rate of change of the colatitude is α and the rate of change of the longitude is β.

Exercise 1.10: Higher-derivative Lagrangians

Derive Lagrange's equations for Lagrangians that depend on accelerations. In particular, show that the Lagrange equations for Lagrangians of the form L(t, q, $$ \dot{q}$$, $$ \ddot{q}$$) with $$ \ddot{q}$$ terms are⁵⁵

$$\begin{array}{ll}{D}^2({\partial }_{3}L\circ \Gamma [q])-D({\partial }_{2}L\circ \Gamma [q])+{\partial }_{1}L\circ \Gamma [q])=0.\hfill & (1.49)\hfill \end{array}$$

In general, these equations, first derived by Poisson, will involve the fourth derivative of q. Note that the derivation is completely analogous to the derivation of the Lagrange equations without accelerations; it is just longer. What restrictions must we place on the variations so that the critical path satisfies a differential equation?

1.5.2 Computing Lagrange's Equations

The procedure for computing Lagrange's equations mirrors the functional expression (1.12), where the procedure Gamma implements Γ:⁵⁶

(define ((Lagrange-equations Lagrangian) q)
  (- (D (compose ((partial 2) Lagrangian) (Gamma q)))
     (compose ((partial 1) Lagrangian) (Gamma q))))

The argument of Lagrange-equations is a procedure that computes a Lagrangian. The Lagrange-equations procedure returns a procedure that when applied to a path q returns a procedure of one argument (time) that computes the left-hand side of the Lagrange equations (1.12). These residual values are zero if q is a path for which the Lagrangian action is stationary.

Observe that the Lagrange-equations procedure, like the Lagrange equations themselves, is valid for any generalized coordinate system. When we write programs to investigate particular systems, the procedures that implement the Lagrangian function and the path q will reflect the actual coordinates chosen to represent the system, but we use the same Lagrange-equations procedure in each case. This abstraction reflects the important fact that the method of derivation of Lagrange's equations from a Lagrangian is always the same; it is independent of the number of degrees of freedom, the topology of the configuration space, and the coordinate system used to describe points in the configuration space.

The free particle

Consider again the case of a free particle. The Lagrangian is implemented by the procedure L-free-particle. Rather than numerically integrating and minimizing the action, as we did in section 1.4, we can check Lagrange's equations for an arbitrary straight-line path t ↦ (at + a₀, bt + b₀, ct + c₀):

(define (test-path t)
  (up (+ (* 'a t) 'a0)
      (+ (* 'b t) 'b0)
      (+ (* 'c t) 'c0)))

(((Lagrange-equations (L-free-particle 'm))
  test-path)
 't)

(down 0 0 0)

That the residuals are zero indicates that the test path satisfies the Lagrange equations.⁵⁷

We can also apply the Lagrange-equations procedure to an arbitrary function:⁵⁸

(show-expression
  (((Lagrange-equations (L-free-particle 'm))
    (literal-function 'x))
   't))
(*(((expt D 2) x) t) m)

The result is an expression containing the arbitrary time t and mass m, so it is zero precisely when D²x = 0, which is the expected equation for a free particle.

The harmonic oscillator

Consider the harmonic oscillator again, with Lagrangian (1.10). We know that the motion of a harmonic oscillator is a sinusoid with a given amplitude, frequency, and phase:

$$\begin{array}{ll}x(t)=a\cos (\omega t+\phi ).\hfill & (1.50)\hfill \end{array}$$

Suppose we have forgotten how the constants in the solution relate to the mass m and spring constant k of the oscillator. Let's plug in the proposed solution and look at the residual:

(define (proposed-solution t)
  (* 'A (cos (+ (* 'omega t) 'phi))))

(show-expression
  (((Lagrange-equations (L-harmonic 'm 'k))
    proposed-solution)
   't))

The residual here shows that for nonzero amplitude, the only solutions allowed are ones where (k − mω²) = 0 or $\omega =\sqrt{k/m}.$

Exercise 1.11: Kepler's third law

A Lagrangian suitable for studying the relative motion of two particles, of masses m₁ and m₂, with potential energy V, is:

(define ((L-central-polar m V) local)
  (let ((q (coordinate local))
        (qdot (velocity local)))
    (let ((r (ref q 0))     (phi (ref q 1))
                            (rdot (ref qdot 0)) (phidot (ref qdot 1)))
      (- (* 1/2 m
            (+ (square rdot) (square (* r phidot))) )
         (V r)))))

The argument m is the reduced mass of the system

$$\begin{array}{ll}m=\frac{m_1{m}_2}{m_1+m_2}.\hfill & (1.51)\hfill \end{array}$$

For gravity, the potential energy function is

(define ((gravitational-energy G m1 m2) r)
  (- (/ (* G m1 m2) r)))

where r is the distance between the two particles.

Consider the simple situation of the particles in circular orbits around their common center of mass. Construct a circular orbit and plug it into the Lagrange equations. Show that the residual gives Kepler's law:

$$\begin{array}{ll}{n}^2{a}^3=G(m_1+m_2)\hfill & (1.52)\hfill \end{array}$$

where n is the angular frequency of the orbit and a is the distance between the particles.

Exercise 1.12: Lagrange's Equations

Compute Lagrange's equations for the Lagrangians in exercise 1.9 using the Lagrange-equations procedure. Additionally, use the computer to perform each of the steps in the Lagrange-equations procedure and show the intermediate results. Relate these steps to the ones you showed in the hand derivation of exercise 1.9.

Exercise 1.13: Higher-derivative Lagrangians

a. Write a procedure to compute the Lagrange equations for Lagrangians that depend upon acceleration, as in exercise 1.10. Note that Gamma can take an optional argument giving the length of the initial segment of the local tuple needed. The default length is 3, giving components of the local tuple up to and including the velocities.

b. Use your procedure to compute the Lagrange equations for the Lagrangian

$L(t,x,v,a)=-\frac{1}{2}mxa-\frac{1}{2}k{x}^2.$

Do you recognize the resulting equation of motion?

c. For more fun, write the general Lagrange equation procedure that takes a Lagrangian that depends on any number of derivatives, and the number of derivatives, to produce the required equations of motion.

1.6   How to Find Lagrangians

Lagrange's equations are a system of second-order differential equations. In order to use them to compute the evolution of a mechanical system, we must find a suitable Lagrangian for the system. There is no general way to construct a Lagrangian for every system, but there is an important class of systems for which we can identify Lagrangians in a straightforward way in terms of kinetic and potential energy. The key idea is to construct a Lagrangian L such that Lagrange's equations are Newton's equations $\overrightarrow{F}=m\overrightarrow{a}$.

Suppose our system consists of N particles indexed by α, with mass m_α and vector position ${\overrightarrow{x}}_{\alpha }(t)$. Suppose further that the forces acting on the particles can be written in terms of a gradient of a potential energy $\mathcal{V}$ that is a function of the positions of the particles and possibly time, but does not depend on the velocities. In other words, the force on particle α is ${\overrightarrow{F}}_{\alpha }=-{\overrightarrow{\nabla }}_{{\overrightarrow{x}}_{\alpha }}\mathcal{V}$, where ${\overrightarrow{\nabla }}_{{\overrightarrow{x}}_{\alpha }}\mathcal{V}$ is the gradient of $\mathcal{V}$ with respect to the position of the particle with index α. We can write Newton's equations as

$$\begin{array}{ll}D(m_{\alpha }D{\overrightarrow{x}}_{\alpha })(t)+{\overrightarrow{\nabla }}_{{\overrightarrow{x}}_{\alpha }}\mathcal{V}(t,{\overrightarrow{x}}_0(t),\dots ,{\overrightarrow{x}}_{N-1}(t))=0.\hfill & (1.53)\hfill \end{array}$$

Vectors can be represented as tuples of components of the vectors on a rectangular basis. So ${\overrightarrow{x}}_1(t)$ is represented as the tuple x₁(t). Let V be the potential energy function expressed in terms of components:

$$\begin{array}{ll}V(t;x_0(t),\dots ,x_{N-1}(t))=\mathcal{V}(t,{\overrightarrow{x}}_0(t),\dots ,{\overrightarrow{x}}_{N-1}(t)).\hfill & (1.54)\hfill \end{array}$$

Newton's equations are

$$\begin{array}{ll}D(m_{\alpha }D{x}_{\alpha })(t)+{\partial }_{1,\alpha }V(t;x_0(t),\dots ,x_{\alpha }(t),\dots ,x_{N-1}(t))=0,\hfill & (1.55)\hfill \end{array}$$

where ∂_1,αV is the partial derivative of V with respect to the x_α(t) argument slot.

To form the Lagrange equations we collect all the position components of all the particles into one tuple x(t), so x(t) = (x₀(t), …, x_N−1(t)). The Lagrange equations for the coordinate path x are

$$\begin{array}{ll}D({\partial }_{2}L\circ \Gamma [x])-{\partial }_{1}L\circ \Gamma [x]=0.\hfill & (1.56)\hfill \end{array}$$

Observe that Newton's equations (1.55) are just the components of the Lagrange equations (1.56) if we choose L to have the properties

$$\begin{array}{lll}({\partial }_{2}L\circ \Gamma [x])(t)\hfill & =[m_{0}D{x}_0(t),\dots ,m_{N-1}D{x}_{N-1}(t)]\hfill & \hfill \\ ({\partial }_{1}L\circ \Gamma [x])(t)\hfill & =[-{\partial }_{1,0}V(t,x(t)),\dots ,-{\partial }_{1,N-1}V(t,x(t))];\hfill & (1.57)\hfill \end{array}$$

here V (t, x(t)) = V (t; x₀(t), …, x_N−1(t)) and ∂_1,αV (t, x(t)) is the tuple of the components of the derivative of V with respect to the coordinates of the particle with index α, evaluated at time t and coordinates x(t). These conditions are satisfied if for every x_α and v_α

$$\begin{array}{ll}{\partial }_{2}L(t;x_0,\dots ,x_{N-1};v_0,\dots ,v_{N-1})\hfill & \hfill \\ =[m_0{v}_0,\dots ,m_{N-1}{v}_{N-1}]\hfill & (1.58)\hfill \end{array}$$

and

$$\begin{array}{ll}{\partial }_{1}L(t;x_0,\dots ,x_{N-1};v_0,\dots ,v_{N-1})\hfill & \hfill \\ =[-{\partial }_{1,0}V(t,x),\dots ,-{\partial }_{1,N-1}V(t,x)],\hfill & (1.59)\hfill \end{array}$$

where x = (x₀, …, x_N−1). One choice for L that has the required properties (1.58–1.59) is

$$\begin{array}{ll}L(t,x,v)=\frac{1}{2}{\displaystyle \sum _{\alpha }{m}_{\alpha }{v}_{\alpha }^2-V(t,x),}\hfill & (1.60)\hfill \end{array}$$

where $v_{\alpha }^2$ is the sum of the squares of the components of v_α.⁵⁹

The first term is the kinetic energy, conventionally denoted T. So this choice for the Lagrangian is L(t, x, v) = T (t, x, v)−V (t, x), the difference of the kinetic and potential energy. We will often extend the arguments of the potential energy function to include the velocities so that we can write L = T − V.⁶⁰

Hamilton's principle

Given a system of point particles for which we can identify the force as the (negative) derivative of a potential energy V that is independent of velocity, we have shown that the system evolves along a path that satisfies Lagrange's equations with L = T − V. Having identified a Lagrangian for this class of systems, we can restate the principle of stationary action in terms of energies. This statement is known as Hamilton's principle: A point-particle system for which the force is derived from a velocity-independent potential energy evolves along a path q for which the action

$S[q](t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}L\circ \Gamma [q]}$

is stationary with respect to variations of the path q that leave the endpoints fixed, where L = T − V is the difference between kinetic and potential energy.⁶¹

It might seem that we have reduced Lagrange's equations to nothing more than $\overrightarrow{F}=m\overrightarrow{a}$, and indeed, the principle is motivated by comparing the two equations for this special class of systems. However, the Lagrangian formulation of the equations of motion has an important advantage over $\overrightarrow{F}=m\overrightarrow{a}$. Our derivation used the rectangular components x_α of the positions of the constituent particles for the generalized coordinates, but if the system's path satisfies Lagrange's equations in some particular coordinate system, it must satisfy the equations in any coordinate system. Thus we see that L = T − V is suitable as a Lagrangian with any set of generalized coordinates. The equations of variational mechanics are derived the same way in any configuration space and any coordinate system. In contrast, the Newtonian formulation is based on elementary geometry: In order for $D^2\overrightarrow{x}(t)$ to be meaningful as an acceleration, $\overrightarrow{x}(t)$ must be a vector in physical space. Lagrange's equations have no such restriction on the meaning of the coordinate q. The generalized coordinates can be any parameters that conveniently describe the configurations of the system.

Constant acceleration

Consider a particle of mass m in a uniform gravitational field with acceleration g. The potential energy is mgh where h is the height of the particle. The kinetic energy is just $\frac{1}{2}m{v}^2$. A Lagrangian for the system is the difference of the kinetic and potential energies. In rectangular coordinates, with y measuring the vertical position and x measuring the horizontal position, the Lagrangian is $L(t;x,y;v_x,v_y)=\frac{1}{2}m\left(v_x^2+v_y^2\right)-mgy$. We have⁶²

(define ((L-uniform-acceleration m g) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (let ((y (ref q 1)))
      (- (* 1/2 m (square v)) (* m g y)))))

(show-expression
  (((Lagrange-equations
      (L-uniform-acceleration 'm 'g))
    (up (literal-function 'x)
        (literal-function 'y)))
   't))

This equation describes unaccelerated motion in the horizontal direction (mD²x(t) = 0) and constant acceleration in the vertical direction (mD²y(t) = −gm).

Central force field

Consider planar motion of a particle of mass m in a central force field, with an arbitrary potential energy U(r) depending only upon the distance r to the center of attraction. We will derive the Lagrange equations for this system in both rectangular coordinates and polar coordinates.

In rectangular coordinates (x, y), with origin at the center of attraction, the potential energy is $V(t;x,y)=U(\sqrt{x^2+y^2})$ and the kinetic energy is $T(t;x,y;v_x,v_y)=\frac{1}{2}m\left(v_x^2+v_y^2\right)$. A Lagrangian for the system is L = T − V:

$$\begin{array}{ll}L(t;x,y;v_x,v_y)=\frac{1}{2}m\left(v_x^2+v_y^2\right)-U(\sqrt{x^2+y^2}).\hfill & (1.61)\hfill \end{array}$$

As a procedure:

(define ((L-central-rectangular m U) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (- (* 1/2 m (square v))
       (U (sqrt (square q))))))

The Lagrange equations are

(show-expression
  (((Lagrange-equations
      (L-central-rectangular 'm (literal-function 'U)))
    (up (literal-function 'x)
        (literal-function 'y)))
   't))

We can rewrite these Lagrange equations as:

$$\begin{array}{ll}m{D}^{2}x(t)=-\frac{x(t)}{r(t)}DU(r(t))\hfill & (1.62)\hfill \end{array}$$

$$\begin{array}{ll}m{D}^{2}y(t)=-\frac{y(t)}{r(t)}DU(r(t)),\hfill & (1.63)\hfill \end{array}$$

where $r(t)=\sqrt{{(x(t))}^2+{(y(t))}^2}$. We can interpret these as follows. The particle is subject to a radially directed force with magnitude −DU(r). Newton's equations equate the force with the product of the mass and the acceleration. The two Lagrange equations are just the rectangular components of Newton's equations.

We can describe the same system in polar coordinates. The relationship between rectangular coordinates (x, y) and polar coordinates (r, φ) is

$$\begin{array}{ll}x=r\cos \phi \hfill & \hfill \\ y=r\sin \phi .\hfill & (1.64)\hfill \end{array}$$

The relationship of the generalized velocities is derived from the coordinate transformation. Consider a configuration path that is represented in both rectangular and polar coordinates. Let $\tilde{x}$ and $\tilde{y}$ be components of the rectangular coordinate path, and let $\tilde{r}$ and $\tilde{\phi }$ be components of the corresponding polar coordinate path. The rectangular components at time t are $(\tilde{x}(t),\tilde{y}(t))$ and the polar coordinates at time t are $(\tilde{r}(t),\tilde{\phi }(t))$. They are related by (1.64):

$$\begin{array}{ll}\tilde{x}(t)=\tilde{r}(t)\cos \tilde{\phi }(t)\hfill & \hfill \\ \tilde{y}(t)=\tilde{r}(t)\sin \tilde{\phi }(t).\hfill & (1.65)\hfill \end{array}$$

The rectangular velocity at time t is $(D\tilde{x}(t),D\tilde{y}(t))$. Differentiating (1.65) gives the relationship among the velocities

$$\begin{array}{lll}D\tilde{x}(t)\hfill & =D\tilde{r}(t)\cos \tilde{\phi }(t)-\tilde{r}(t)D\tilde{\phi }(t)\sin \tilde{\phi }(t)\hfill & \hfill \\ D\tilde{y}(t)\hfill & =D\tilde{r}(t)\sin \tilde{\phi }(t)+\tilde{r}(t)D\tilde{\phi }(t)\cos \tilde{\phi }(t).\hfill & (1.66)\hfill \end{array}$$

These relations are valid for any configuration path at any moment, so we can abstract them to relations among coordinate representations of an arbitrary velocity. Let v_x and v_y be the rectangular components of the velocity and $\dot{r}$ and $\dot{\phi }$ be the rate of change of r and φ. Then

$$\begin{array}{ll}{v}_x=\dot{r}\cos \phi -r\dot{\phi }\sin \phi \hfill & \hfill \\ v_y=\dot{r}\sin \phi +r\dot{\phi }\cos \phi .\hfill & (1.67)\hfill \end{array}$$

The kinetic energy is $\frac{1}{2}m(v_x^2+v_y^2)$:

$$\begin{array}{ll}T(t;r,\phi ;\dot{r},\dot{\phi })=\frac{1}{2}m\left({\dot{r}}^2+r^2{\dot{\phi }}^2\right),\hfill & (1.68)\hfill \end{array}$$

and the Lagrangian is

$$\begin{array}{ll}L(t;r,\phi ;\dot{r},\dot{\phi })=\frac{1}{2}m\left({\dot{r}}^2+r^2{\dot{\phi }}^2\right)-U(r).\hfill & (1.69)\hfill \end{array}$$

We express this Lagrangian as follows:

(define ((L-central-polar m U) local)
  (let ((q (coordinate local))
        (qdot (velocity local)))
    (let ((r (ref q 0)) (phi (ref q 1))
                        (rdot (ref qdot 0)) (phidot (ref qdot 1)))
      (- (* 1/2 m
            (+ (square rdot)
               (square (* r phidot))) )
         (U r)))))

Lagrange's equations are

(show-expression
  (((Lagrange-equations
      (L-central-polar 'm (literal-function 'U)))
    (up (literal-function 'r)
        (literal-function 'phi)))
   't))

We can interpret the first equation as saying that the product of the mass and the radial acceleration is the sum of the force due to the potential and the centrifugal force. The second equation can be interpreted as saying that the derivative of the angular momentum mr²Dφ is zero, so angular momentum is conserved.

Note that we used the same Lagrange-equations procedure for the derivation in both coordinate systems. Coordinate representations of the Lagrangian are different for different coordinate systems, and the Lagrange equations in different coordinate systems look different. Yet the same method is used to derive the Lagrange equations in any coordinate system.

Exercise 1.14: Coordinate-independence of Lagrange equations

Check that the Lagrange equations for central force motion in polar coordinates and in rectangular coordinates are equivalent. Determine the relationship among the second derivatives by substituting paths into the transformation equations and computing derivatives, then substitute these relations into the equations of motion.

1.6.1 Coordinate Transformations

The motion of a system is independent of the coordinates we use to describe it. This coordinate-free nature of the motion is apparent in the action principle. The action depends only on the value of the Lagrangian along the path and not on the particular coordinates used in the representation of the Lagrangian. We can use this property to find a Lagrangian in one coordinate system in terms of a Lagrangian in another coordinate system.

Suppose we have a mechanical system whose motion is described by a Lagrangian L that depends on time, coordinates, and velocities. And suppose we have a coordinate transformation F such that x = F (t, x′). The Lagrangian L is expressed in terms of the unprimed coordinates. We want to find a Lagrangian L′ expressed in the primed coordinates that describes the same system. One way to do this is to require that the value of the Lagrangian along any configuration path be independent of the coordinate system. If q is a path in the unprimed coordinates and q′ is the corresponding path in primed coordinates, then the Lagrangians must satisfy:

$$\begin{array}{ll}L\prime \circ \Gamma [q\prime ]=L\circ \Gamma [q].\hfill & (1.70)\hfill \end{array}$$

We have seen that the transformation from rectangular to polar coordinates implies that the generalized velocities transform in a certain way. The velocity transformation can be deduced from the requirement that a path in polar coordinates and a corresponding path in rectangular coordinates are consistent with the coordinate transformation. In general, the requirement that paths in two different coordinate systems be consistent with the coordinate transformation can be used to deduce how all of the components of the local tuple transform. Given a coordinate transformation F, let C be the corresponding function that maps local tuples in the primed coordinate system to corresponding local tuples in the unprimed coordinate system:

$$\begin{array}{ll}C\circ \Gamma [q\prime ]=\Gamma [q].\hfill & (1.71)\hfill \end{array}$$

We will deduce the general form of C below.

Given such a local-tuple transformation C, a Lagrangian L′ that satisfies equation (1.70) is

$$\begin{array}{ll}L\prime =L\circ C.\hfill & (1.72)\hfill \end{array}$$

We can see this by substituting for L′ in equation (1.70):

$$\begin{array}{ll}L\prime \circ \Gamma [q\prime ]=L\circ C\circ \Gamma [q\prime ]=L\circ \Gamma [q].\hfill & (1.73)\hfill \end{array}$$