Marcos Duarte

Laboratory of Biomechanics and Motor Control (http://demotu.org/)

Federal University of ABC, Brazil

"The theoretical development of the laws of motion of bodies is a problem of such interest and importance that it has engaged the attention of all the most eminent mathematicians since the invention of dynamics as a mathematical science by Galileo, and especially since the wonderful extension which was given to that science by Newton.

Among the successors of those illustrious men, Lagrange has perhaps done more than any other analyst to give extent and harmony to such deductive researches, by showing that the most varied consequences respecting the motions of systems of bodies may be derived from one radical formula; the beauty of the methods so suiting the dignity of the results as to make of his great work a kind of scientific poem."

— Hamilton, 1834 (apud Taylor, 2005).

The Lagrangian mechanics is a formulation of classical mechanics where the equations of motion are obtained from the kinetic and potential energy of the system (scalar quantities) represented in generalized coordinates instead of using Newton's laws of motion to deduce the equations of motion from the forces on the system (vector quantities) represented in Cartesian coordinates. Lagrangian mechanics was introduced by Joseph-Louis Lagrange in the late 18th century.

Let's deduce the Lagrange equations, but first let's review the basics of the Newtonian approach.

One can describe the motion of a particle by specifying its position with respect to a frame of reference in the three-dimensional space as a function of time:

\begin{equation} x(t),\, y(t),\, z(t) \quad \equiv \quad x_i(t) \quad i=1,\dotsc,3 \label{eq1} \end{equation}A system of $N$ particles will require $3N$ equations to describe their motion.

The basic problem in classical mechanics is to find ways to determine functions such as these, also known as equations of motion, capable of specifying the position of objects over time, for any mechanical situation. Assuming as known the meaning of $x_i(t)$, one can define the components of velocity, $v_i$, and acceleration, $a_i$, at time $t$, as:

\begin{equation} \begin{array}{rcl} v_i(t) = \dfrac{\mathrm d x_i(t)}{\mathrm d t} = \dot{x}_i(t) \\ a_i(t) = \dfrac{\mathrm d^2 x_i(t)}{\mathrm d t^2} = \dot{v}_i(t) \label{eq3} \end{array} \end{equation}Where we used the Newton's notation for differentiation (also called the dot notation), a dot over the dependent variable. Of note, Joseph Louis Lagrange introduced the prime mark to denote a derivative: $x'(t)$. Read more about the different notations for differentiation at Wikipedia.

The Newton's laws of motion laid the foundation for classical mechanics. They describe the relationship between the motion of a body and the possible forces acting upon it.

In Newtonian mechanics, the body's linear momentum is defined as:

\begin{equation} \mathbf{p} = m\mathbf{v} \label{eq5} \end{equation}If the mass of the body is constant, remember that Newton's second law can be expressed by:

\begin{equation} \mathbf{F} = \frac{\mathrm d \mathbf{p}}{\mathrm d t}=\frac{\mathrm d \big(m\mathbf{v}\big)}{\mathrm d t} = m\mathbf{a} \label{eq6} \end{equation}Using Newton's second law, to determine the position of the body we will have to solve the following second order ordinary differential equation:

\begin{equation} \frac{\mathrm d^2 x_i(t)}{\mathrm d t^2} = \frac{\mathbf{F}}{m} \label{eqa7} \end{equation}Which has the general solution:

\begin{equation} \mathbf{x}(t) = \int\!\bigg(\int\frac{\mathbf{F}}{m} \mathrm{d}t\bigg)\mathrm{d}t \label{eq8} \end{equation}A related physical quantity is the mechanical energy, which is the sum of kinetic and potential energies.

The kinetic energy, $T$ of a particle is given by:

The kinetic energy of a particle can be expressed in terms of its linear momentum:

\begin{equation} T = \frac{p^2}{2m} \label{eq10} \end{equation}And for a given coordinate of the particle's motion, its linear momentum can be obtained from its kinetic energy by:

\begin{equation} p_i = \frac{\partial T}{\partial v_i} \label{eq11} \end{equation}The potential energy, $V$ is the stored energy of a particle and its formulation is dependent on the force acting on the particle. For example, for a conservative force dependent solely on the particle position, such as due to the gravitational field near the Earth surface or due to a linear spring, force and potential energy are related by:

\begin{equation} \mathbf{F} = - \frac{\partial \mathbf{V}}{\partial x} \label{eq12} \end{equation}The Lagrangian mechanics can be formulated independent of the Newtonian mechanics and Cartesian coordinates; in fact Joseph-Louis Lagrange this new formalism from the principle of least action.

For simplicity, let's first deduce the Lagrange's equation for a particle in Cartesian Coordinates and from Newton's second law.

Because we want to deduce the laws of motion based on the mechanical energy of the particle, one can see that the time derivative of the expression for the linear momentum as a function of the kinetic energy, cf. Eq. (\ref{eq11}), is equal to the force acting on the particle and we can substitute the force in Newton's second law by this term:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\bigg(\frac{\partial T}{\partial \dot x}\bigg) = m\ddot x \label{eq13} \end{equation}We saw that a conservative force can also be expressed in terms of the potential energy of the particle, cf. Eq. (\ref{eq12}); substituting the right side of the equation above by this expression, we have:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\bigg(\frac{\partial T}{\partial \dot x}\bigg) = -\frac{\partial V}{\partial x} \label{eq14} \end{equation}Using the fact that:

\begin{equation} \frac{\partial T}{\partial x} = 0 \quad and \quad \frac{\partial V}{\partial \dot x} = 0 \label{eq15} \end{equation}We can write:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\bigg(\frac{\partial (T-V)}{\partial \dot x}\bigg) - \frac{\partial (T-V)}{\partial x} = 0 \label{eq16} \end{equation}Defining the Lagrange or Lagrangian function, $\mathcal{L}$, as the difference between the kinetic and potential energy in the system:

\begin{equation} \mathcal{L}(x,\dot x, t) = T - V \label{eq17} \end{equation}We have the Lagrange's equation in Cartesian Coordinates for a conservative force acting on a particle:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\bigg(\frac{\partial \mathcal{L}}{\partial \dot x}\bigg) - \frac{\partial \mathcal{L}}{\partial x} = 0 \label{eq18} \end{equation}Once all derivatives of the Lagrangian function are calculated, this equation will be the equation of motion for the particle. If there are $N$ independent particles in a three-dimensional space, there will be $3N$ equations for the system.

The set of equations above for a system are known as Euler–Lagrange equations, or Lagrange's equations of the second kind.

Let's see some simple examples of the Lagrange's equation in Cartesian Coordinates just to consolidate what was deduced above. The real application of Lagrangian mechanics is in generalized coordinates, what will see after these examples.

Let's deduce the equation of motion using the Lagrangian mechanics for a particle with mass $m$ moving in the 3D space under the influence of a conservative force.

The Lagrangian $(\mathcal{L} = T - V)$ of the particle is:

The equations of motion for the particle are found by applying the Lagrange's equation for each coordinate.

For the x coordinate:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\left( {\frac{\partial \mathcal{L}}{\partial \dot{x}}} \right) = \frac{\partial \mathcal{L}}{\partial x } \end{equation}And the derivatives of $\mathcal{L}$ are given by:

\begin{equation} \begin{array}{rcl} &\dfrac{\partial \mathcal{L}}{\partial x} &=& -\dfrac{\partial \mathbf{V}}{\partial x} \\ &\dfrac{\partial \mathcal{L}}{\partial \dot{x}} &=& m\dot{x} \\ &\dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{x}}} \right) &=& m\ddot{x} \end{array} \end{equation}hence:

\begin{equation} m\ddot{x} = -\frac{\partial \mathbf{V}}{\partial x} \end{equation}and similarly for the $y$ and $z$ coordinates.

For instance, if the conservative force is due to the gravitational field near Earth's surface $(\mathbf{V}=[0, mgy, 0])$, the Lagrange's equations (the equations of motion) are:

\begin{equation} \begin{array}{rcl} m\ddot{x} &=& -\dfrac{\partial (0)}{\partial x} = 0 \\ m\ddot{y} &=& -\dfrac{\partial (mgy)}{\partial y} = -mg \\ m\ddot{z} &=& -\dfrac{\partial (0)}{\partial z} = 0 \end{array} \end{equation}Consider a system with a mass $m$ attached to an ideal spring (massless) with spring constant $k$ at the horizontal direction $x$. If the system is perturbed (a force is momentarily applied to the mass), we know the mass will oscillate around the rest position of the spring. Let's deduce the equation of motion for this system using the Lagrangian mechanics.

The Lagrangian $(\mathcal{L} = T - V)$ of the system is:

The derivatives of $L$ are given by:

\begin{equation} \begin{array}{rcl} &\dfrac{\partial \mathcal{L}}{\partial x} &=& -kx \\ &\dfrac{\partial \mathcal{L}}{\partial \dot{x}} &=& m\dot{x} \\ &\dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{x}}} \right) &=& m\ddot{x} \end{array} \end{equation}And the Lagrange's equation (the equation of motion) is:

\begin{equation} m\ddot{x} + kx = 0 \end{equation}The direct application of Newton's laws to mechanical systems results in a set of equations of motion in terms of Cartesian coordinates of each of the particles that make up the system. In many cases, this is not the most convenient coordinate system to solve the problem or describe the movement of the system. For example, in problems involving many particles, it may be convenient to choose a system that includes the coordinate of the center of mass. Another example is a serial chain of rigid links, such as a member of the human body or from a robot manipulator, it may be simpler to describe the positions of each link by the angles between links.

Coordinate systems such as these are referred as generalized coordinates. Generalized coordinates uniquely specify the positions of the particles in a system. Although there may be several generalized coordinates to describe a system, usually a judicious choice of generalized coordinates provides the minimum number of independent coordinates that define the configuration of a system (which is the number of degrees of freedom of the system), turning the problem simpler to solve.

Being a little more technical, according to Wikipedia):

"In classical mechanics, the parameters that define the configuration of a system are called generalized coordinates, and the vector space defined by these coordinates is called the configuration space of the physical system. It is often the case that these parameters satisfy mathematical constraints, such that the set of actual configurations of the system is a manifold in the space of generalized coordinates. This manifold is called the configuration manifold of the system."

In problems where it is desired to use generalized coordinates, one can write Newton's equations of motion in terms of Cartesian coordinates and then transform them into generalized coordinates. However, it would be desirable and convenient to have a general method that would directly establish the equations of motion in terms of a set of convenient generalized coordinates. In addition, general methods for writing, and perhaps solving, the equations of motion in terms of any coordinate system would also be desirable. The Lagrangian mechanics is such a method.

We have deduced the Lagrange's equation in Cartesian Coordinates from the Newton's law just because it was a simple form of getting to the final equations. But, by no means the Lagrangian Mechanics should be viewed as a consequence of Newton's laws and specific to Cartesian Coordinates. The Lagrangian Mechanics could be deduced completely independent of Newton's law.

The Lagrange's equation can be expressed in terms of generalized coordinates what makes the Lagrangian formalism even more powerful. In fact, we will have the same equation as we deduced before (the only explicit difference will be that we will use $q_i$ instead of the Cartesian coordinate).

See this notebook for a deduction of the Lagrange's equation in generalized coordinates.

Defining the Lagrange or Lagrangian function of a system with $N$ generalized coordinates:

\begin{equation} \mathcal{L} \equiv \mathcal{L}(q_1,\dotsc,q_{N} ,\dot{q}_1,\dotsc,\dot{q}_{N} ) = T - V \label{eq46} \end{equation}We have the Lagrange's equation:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\left( {\frac{\partial \mathcal{L}}{\partial \dot{q}_i }} \right)-\frac{\partial \mathcal{L}}{\partial q_i } = Q_{NCi} \quad i=1,\dotsc,N \label{eq47} \end{equation}Where $Q_{NCi}$ is the generalized force due to a non-conservative force, any force that can't be expressed in terms of a potential.

Once all derivatives of the Lagrangian function are calculated, this equation will be the equation of motion for each particle. If there are $N$ generalized coordinates to define the configuration of a system, there will be $N$ equations for the system.

The set of equations above for a system are known as Euler–Lagrange equations, or Lagrange's equations of the second kind.

Consider a pendulum with a massless rod of length $d$ and a mass $m$ at the extremity swinging in a plane forming the angle $\theta$ with vertical, $g=10 m/s^2$ and the origin of the coordinate system at the point of the pendulum suspension.

Although the pendulum moves at the plane, it only has one degree of freedom, which can be described by the angle $\theta$, the generalized coordinate. It is not difficult to find the equation of motion using Newton's law, but let's find it using the Lagrangian mechanics.

The kinetic energy is:

\begin{equation} T = \frac{1}{2}md^2\dot\theta^2 \end{equation}And the potential energy is:

\begin{equation} V = -mgd\cos\theta \end{equation}The Lagrangian function is:

\begin{equation} \mathcal{L}(\theta, \dot\theta, t) = \frac{1}{2}md^2\dot\theta^2(t) + mgd\cos\theta(t) \end{equation}And the derivatives are given by:

\begin{equation} \begin{array}{rcl} &\dfrac{\partial \mathcal{L}}{\partial \theta} &=& -mgd\sin\theta \\ &\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}} &=& md^2\dot{\theta} \\ &\dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}}} \right) &=& md^2\ddot{\theta} \end{array} \end{equation}Finally, the Lagrange's equation (the equation of motion) is:

\begin{equation} md^2\ddot\theta + mgd\sin\theta = 0 \end{equation}A classical approach to solve the equation of motion for the simple pendulum is to consider the motion for small angles where $\sin\theta \approx \theta$ and the differential equation is linearized to $d\ddot\theta + g\theta = 0$. This equation has an analytical solution of the type $\theta(t) = A \sin(\omega t + \phi)$, where $\omega = \sqrt{g/d}$ and $A$ and $\phi$ are constants related to the initial position and velocity.

Let's solve the differential equation for the pendulum numerically using the Euler’s method.

Remember that we have to (1) transform the second-order ODE into two coupled first-order ODEs, (2) approximate the derivative of each variable by its discrete first order difference and (3) write equation to calculate the variable in a recursive way, updating its value with an equation based on the first order difference.

In Python:

In [1]:

```
import numpy as np
def pendulum_euler(T, y0, v0, h):
"""
Two coupled first-order ODEs for the pendulum (y=theta and v is a new variable):
dydt = v
dvdt = -g/d * y
Two equations to update the values of the variables based on first-order difference:
y[n+1] = y[n] + h*v[n]
v[n+1] = v[n] + h*dvdt[n]
"""
N = int(np.ceil(T/h))
y, v = np.zeros(N), np.zeros(N)
y[0], v[0] = y0, v0
d = 2 # length of the pendulum in m
g = 10 # acceleration of gravity in m/s2
for i in range(N-1):
y[i+1] = y[i] + h*v[i]
v[i+1] = v[i] + h*(-g/d*y[i])
t = np.linspace(0, T, N, endpoint=False)
return t, y, v
```

In [2]:

```
%matplotlib notebook
import matplotlib.pyplot as plt
def plot(t, y, v):
"""
Plot data
"""
fig, ax = plt.subplots(1, 1, figsize=(8, 3))
ax.plot(t, y, 'b', label='Position')
ax.plot(t, v, 'r', label='Velocity')
ax.legend()
plt.show()
```

In [3]:

```
T, y0, v0, h = 10, 45*np.pi/180, 0, .001
t, theta, vtheta = pendulum_euler(T, y0, v0, h)
plot(t, theta, vtheta)
```

Consider a simple pendulum with massless rod of length $d$ and mass $m$ at the extremity of the rod forming an angle $\theta$ with the vertical direction under the action of gravity and $g=10 m/s^2$. The pendulum swings freely from a cart with mass $M$ that moves at the horizontal direction pushed by a force $F_x$.

Let's use the Lagrangian mechanics to derive the equations of motion for the system.

From the figure on the right, because of the constraints introduced by the constant length of the rod and the motion the cart can perform, good generalized coordinates to describe the configuration of the system are $x$ and $\theta$. Let's use these coordinates.

The positions of the cart (c) and of the pendulum tip (p) are:

$ x_c = x $

$ y_c = 0 $

$ x_p = x + d \sin \theta $

$ y_p = -d \cos \theta $

The velocities of the cart and of the pendulum are:

$ \dot{x}_c = \dot{x} $

$ \dot{y}_c = 0 $

$ \dot{x}_p = \dot{x} + d \dot{\theta} \cos \theta $

$ \dot{y}_p = d \dot{\theta} \sin \theta $

The kinetic energy of the system is:

\begin{equation} T = \frac{1}{2} M \big(\dot x_c^2 + \dot y_c^2 \big) + \frac{1}{2}m\,\big( \dot x_p^2 + \dot y_p^2 \big) \end{equation}And the potential energy of the system is:

\begin{equation} V = M g y_c + m g y_p \end{equation}The Lagrangian function is:

\begin{equation}\begin{array}{rcl} \mathcal{L} &=& \frac{1}{2} M \dot x^2 + \frac{1}{2}m\,\bigg[ \big(\dot{x} + d \dot{\theta} \cos\theta \big)^2 + \big(d \dot{\theta} \sin \theta \big)^2 \bigg] + m g d \cos \theta \\ &=& \frac{1}{2} (M+m) \dot x^2 + m\dot{x}d\dot{\theta}\cos\theta + \frac{1}{2}md^2\dot{\theta}^2 + m g d \cos \theta \end{array}\end{equation}The derivatives w.r.t. $x$ are:

\begin{equation}\begin{array}{rcl} \dfrac{\partial \mathcal{L}}{\partial x} &=& 0 \\ \dfrac{\partial \mathcal{L}}{\partial \dot{x}} &=& (M+m) \dot{x} + m d \dot{\theta} \cos \theta \\ \dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{x}}} \right) &=& (M+m)\ddot{x} + m d \ddot{\theta} \cos \theta - m d \dot{\theta}^2 \sin \theta \end{array}\end{equation}The derivatives w.r.t. $\theta$ are:

\begin{equation}\begin{array}{rcl} \dfrac{\partial \mathcal{L}}{\partial \theta} &=& -m\dot{x}d\dot{\theta}\sin\theta - mgd\sin\theta \\ \dfrac{\partial \mathcal{L}}{\partial \dot{\theta}} &=& m\dot{x}d\cos\theta + md^2\dot{\theta} \\ \dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}}} \right) &=& m\ddot{x}d\cos\theta - m\dot{x}d\dot{\theta}\sin\theta + md^2\ddot{\theta} \end{array}\end{equation}Finally, the Lagrange's equations (the equations of motion) are:

\begin{equation}\begin{array}{rcl} (M+m)\ddot{x} + md\big(\ddot{\theta} \cos\theta - \dot{\theta}^2 \sin\theta\big) = F_x \\ \ddot{\theta} + \dfrac{\ddot{x}}{d}\cos\theta + \dfrac{g}{d}\sin\theta = 0 \end{array}\end{equation}Consider a double pendulum (one pendulum attached to another) with massless rods of length $d_1$ and $d_2$ and masses $m_1$ and $m_2$ at the extremities of each rod swinging in a plane forming the angles $\theta_1$ and $\theta_2$ with vertical and $g=10 m/s^2$.

This case could be solved using Newtonian mechanics, but it's not simple (e.g., see this link). Instead, let's use the Lagrangian mechanics to derive the equations of motion for each pendulum.

The system has two degrees of freedom and we need two generalized coordinates ($\theta_1, \theta_2$) to describe the system's configuration.

The position of masses $m_1$ and $m_2$ are:

$x_1 = d_1\sin\theta_1$

$y_1 = -d_1\cos\theta_1$

$x_2 = d_1\sin\theta_1 + d_2\sin\theta_2$

$y_2 = -d_1\cos\theta_1 - d_2\cos\theta_2$

The kinetic and potential energies of the system are:

$ T = \frac{1}{2}m_1(\dot x_1^2 + \dot y_1^2) + \frac{1}{2}m_2(\dot x_2^2 + \dot y_2^2) $

$ V = m_1gy_1 + m_2gy_2 $

Let's use Sympy to help us; in fact we could solve this problem entirely in Sympy, see Lagrange’s method in Sympy, but for now let's do just the derivatives.

Let's import Sympy libraries and define some variables:

In [4]:

```
from sympy import Symbol, symbols, cos, sin, Matrix, simplify, Eq, latex
from sympy.physics.mechanics import dynamicsymbols, mlatex, init_vprinting
init_vprinting()
from IPython.display import display, Math
t = Symbol('t')
d1, d2, m1, m2, g = symbols('d1 d2 m1 m2 g', positive=True)
a1, a2 = dynamicsymbols('theta1 theta2')
```

The positions and velocities of masses $m_1$ and $m_2$ are:

In [5]:

```
x1, y1 = d1*sin(a1), -d1*cos(a1)
x2, y2 = d1*sin(a1) + d2*sin(a2), -d1*cos(a1) - d2*cos(a2)
x1d, y1d = x1.diff(t), y1.diff(t)
x2d, y2d = x2.diff(t), y2.diff(t)
display(Math(r'x_1=' + mlatex(x1) + r'\quad \text{and} \quad \dot{x}_1=' + mlatex(x1d)))
display(Math(r'y_1=' + mlatex(y1) + r'\quad \text{and} \quad \dot{y}_1=' + mlatex(y1d)))
display(Math(r'x_2=' + mlatex(x2) + r'\quad \text{and} \quad \dot{x}_2=' + mlatex(x2d)))
display(Math(r'y_2=' + mlatex(y2) + r'\quad \text{and} \quad \dot{y}_2=' + mlatex(y2d)))
```

The kinetic and potential energies of the system are:

In [6]:

```
T = m1*(x1d**2 + y1d**2)/2 + m2*(x2d**2 + y2d**2)/2
V = m1*g*y1 + m2*g*y2
display(Math(r'T=' + mlatex(simplify(T))))
display(Math(r'V=' + mlatex(simplify(V))))
```

The Lagrangian function is:

In [7]:

```
L = T - V
display(Math(r'\mathcal{L}=' + mlatex(simplify(L))))
```

And the derivatives are (let's write a function to automate this process):

In [8]:

```
def lagrange_terms(L, *q, show=True):
"""Calculate terms of Lagrange equations given the Lagrangian and q's.
"""
Lterms = []
for qi in q:
dLdqi = simplify(L.diff(qi))
Lterms.append(dLdqi)
dLdqdi = simplify(L.diff(qi.diff(t)))
Lterms.append(dLdqdi)
dtdLdqdi = simplify(dLdqdi.diff(t))
Lterms.append(dtdLdqdi)
if show:
display(Math(r'w.r.t.\;%s:'%latex(qi.func)))
display(Math(r'\dfrac{\partial\mathcal{L}}{\partial %s}='
%latex(qi.func) + mlatex(dLdqi)))
display(Math(r'\dfrac{\partial\mathcal{L}}{\partial\dot{%s}}='
%latex(qi.func) + mlatex(dLdqdi)))
display(Math(r'\dfrac{\mathrm d}{\mathrm{dt}}\left({\dfrac{'+
r'\partial\mathcal{L}}{\partial\dot{%s}}}\right)='
%latex(qi.func) + mlatex(dtdLdqdi)))
return Lterms
```

In [9]:

```
Lterms = lagrange_terms(L, a1, a2)
```

Finally, the Lagrange's equations (the equations of motion) are:

In [10]:

```
for i in range(int(len(Lterms)/3)):
display(Eq(simplify(Lterms[3*i+2]-Lterms[3*i]), 0))
```

The motion of a double pendulum is very interesting; most of times it presents a chaotic behavior, see for example https://www.myphysicslab.com/pendulum/double-pendulum-en.html.

In order to solve numerically the ODEs for the double pendulum we will transform each equation above into two first ODEs. As the two variables $\ddot{\theta_1}$ and $\ddot{\theta_2}$ appear in both equations, first we have to rearrange the equations to find expressions for $\ddot{\theta_1}$ and $\ddot{\theta_2}$. Second, we define two variables and then we can use these equations as we did for the single pendulum.

But we should avoid to use the Euler's method because of the non-negligible error in the numerical integration in this case; more accurate methods such as Runge-Kutta should be employed.

Consider the double compound pendulum shown on the the right with length $d$ and mass $m$ of each rod swinging in a plane forming the angles $\theta_1$ and $\theta_2$ with vertical and $g=10 m/s^2$.

The system has two degrees of freedom and we need two generalized coordinates ($\theta_1, \theta_2$) to describe the system's configuration.

Let's use the Lagrangian mechanics to derive the equations of motion for each pendulum.

To calculate the potential and kinetic energy of the system, we will need to calculate the position and velocity of each pendulum. Now each pendulum is a rod with distributed mass and we will have to calculate the moment of rotational inertia of the rod. In this case, the kinetic energy of each pendulum will be given as the kinetic energy due to rotation of the pendulum plus the kinetic energy due to the speed of the center of mass of the pendulum, such that the total kinetic energy of the system is:

$ T = \overbrace{\underbrace{\frac{1}{2}I_{cm}\dot\theta_1^2}_{\text{rotation}} + \underbrace{\frac{1}{2}m(\dot x_{1,cm}^2 + \dot y_{1,cm}^2)}_{\text{translation}}}^{\text{pendulum 1}} + \overbrace{\underbrace{\frac{1}{2}I_{cm}\dot\theta_2^2}_{\text{rotation}} + \underbrace{\frac{1}{2}m(\dot x_{2,cm}^2 + \dot y_{2,cm}^2)}_{\text{translation}}}^{\text{pendulum 2}} $

And the potential energy of the system is:

$ V = mg\big(y_{1,cm} + y_{2,cm}\big) $

Let's use Sympy once again.

The position and velocity of the center of mass of the rods $1$ and $2$ are:

In [11]:

```
d, m, g = symbols('d m g', positive=True)
a1, a2 = dynamicsymbols('theta1 theta2')
I = m*d*d/12
x1, y1 = d*sin(a1)/2, -d*cos(a1)/2
x2, y2 = d*sin(a1) + d*sin(a2)/2, -d*cos(a1) - d*cos(a2)/2
x1d, y1d = x1.diff(t), y1.diff(t)
x2d, y2d = x2.diff(t), y2.diff(t)
display(Math(r'x_1=' + mlatex(x1) + r'\quad \text{and} \quad \dot{x}_1=' + mlatex(x1d)))
display(Math(r'y_1=' + mlatex(y1) + r'\quad \text{and} \quad \dot{y}_1=' + mlatex(y1d)))
display(Math(r'x_2=' + mlatex(x2) + r'\quad \text{and} \quad \dot{x}_2=' + mlatex(x2d)))
display(Math(r'y_2=' + mlatex(y2) + r'\quad \text{and} \quad \dot{y}_2=' + mlatex(y2d)))
```

The kinetic and potential energies of the system are:

In [12]:

```
T = I/2*(a1.diff(t))**2 + m*(x1d**2+y1d**2)/2 + I/2*(a2.diff(t))**2 + m*(x2d**2+y2d**2)/2
V = m*g*y1 + m*g*y2
display(Math(r'T=' + mlatex(simplify(T))))
display(Math(r'V=' + mlatex(simplify(V))))
```

The Lagrangian function is:

In [13]:

```
L = T - V
display(Math(r'\mathcal{L}=' + mlatex(simplify(L))))
```

And the derivatives are (let's write a function to automate this process):

In [14]:

```
Lterms = lagrange_terms(L, a1, a2)
```

Finally, the Lagrange's equations (the equations of motion) are:

In [15]:

```
for i in range(int(len(Lterms)/3)):
display(Eq(simplify(Lterms[3*i+2]-Lterms[3*i]), 0))
```

Consider a system composed by two masses and two springs attached in series with massless springs under gravity and a force on $m_2$ and $g=10 m/s^2$ as shown in the figure.

The system has two degrees of freedom and we need two generalized coordinates to describe the system's configuration, for example, ${y_1, y_2}$, the positions of masses $m_1, m_2$ w.r.t. the ceiling (the origin).

(We could also have used ${z_1, z_2}$, the position of mass $m_1$ w.r.t. the ceiling and the position of mass $m_2$ w.r.t. the mass $m_1$. But the first set of coordinates would be the only choice if we were solving the problem using Newtonian mechanics because we can only apply Newton's second law to an inertial frame of reference.)

The kinetic energy of the system is:

$ T = \frac{1}{2}m_1\dot y_1^2 + \frac{1}{2}m_2\dot y_2^2 $

The potential energy of the system is:

$ V = \frac{1}{2}k_1 (y_1-\ell_1)^2 + \frac{1}{2}k_2 ((y_2-y_1)-\ell_2)^2 - m_1gy_1 - m_2g y_2 $

Where $\ell_1, \ell_2$ are the resting lengths (constants) of the two springs because the elastic potential energy is proportional to the deformation of the spring. For simplicity, let's ignore the resting position of each spring, so:

$ V = \frac{1}{2}k_1 y_1^2 + \frac{1}{2}k_2 (y_2-y_1)^2 - m_1gy_1 - m_2g y_2 $

Sympy is our friend:

In [16]:

```
from sympy import Symbol, symbols, cos, sin, Matrix, simplify, Eq, latex
from sympy.physics.mechanics import dynamicsymbols, mlatex, init_vprinting
init_vprinting()
from IPython.display import display, Math
t = Symbol('t')
m1, m2, g, k1, k2 = symbols('m1 m2 g k1 k2', positive=True)
y1, y2, F = dynamicsymbols('y1 y2 F')
```

The Lagrangian function is:

In [17]:

```
y1d, y2d = y1.diff(t), y2.diff(t)
T = (m1*y1d**2)/2 + (m2*y2d**2)/2
V = (k1*y1**2)/2 + (k2*(y2-y1)**2)/2 - m1*g*y1 - m2*g*y2
#display(Math(r'T=' + mlatex(simplify(T))))
#display(Math(r'V=' + mlatex(simplify(V))))
L = T - V
display(Math(r'\mathcal{L}=' + mlatex(simplify(L))))
```

And the derivatives are (using the function we wrote before to automate this process):

In [18]:

```
Lterms = lagrange_terms(L, y1, y2)
```

Finally, the Lagrange's equations (the equations of motion) are:

In [19]:

```
display(Eq(simplify(Lterms[2]-Lterms[0]), 0))
display(Eq(simplify(Lterms[5]-Lterms[3]), F))
```

**Same problem, but with the other set of coordinates**

Using ${z_1, z_2}$ as the position of mass $m_1$ w.r.t. the ceiling and the position of mass $m_2$ w.r.t. the mass $m_1$, the solution is:

In [20]:

```
z1, z2 = dynamicsymbols('z1 z2')
z1d, z2d = z1.diff(t), z2.diff(t)
T = (m1*z1d**2)/2 + (m2*(z1d + z2d)**2)/2
V = (k1*z1**2)/2 + (k2*z2**2)/2 - m1*g*z1 - m2*g*(z1+z2)
display(Math(r'T=' + mlatex(simplify(T))))
display(Math(r'V=' + mlatex(simplify(V))))
L = T - V
display(Math(r'\mathcal{L}=' + mlatex(simplify(L))))
```

In [21]:

```
Lterms = lagrange_terms(L, z1, z2)
```

Finally, the Lagrange's equations (the equations of motion) are:

In [22]:

```
display(Eq(simplify(Lterms[2]-Lterms[0]), F))
display(Eq(simplify(Lterms[5]-Lterms[3]), F))
```

The solutions using the two sets of coordinates seem different; the reader is invited to verify that in fact they are the same (remember that $y_1 = z_1,\, y_2 = z_1+z_2,\, \ddot{y}_2 = \ddot{z}_1+\ddot{z}_2$).

**Same problem, but considering the spring resting length**

And here is the solution considering the rest length of the spring:

In [23]:

```
t = Symbol('t')
m1, m2, g, k1, k2 = symbols('m1 m2 g k1 k2', positive=True)
y1, y2, l1, l2, F = dynamicsymbols('y1 y2 ell1 ell2 F')
y1d, y2d = y1.diff(t), y2.diff(t)
T = (m1*y1d**2)/2 + (m2*y2d**2)/2
V = (k1*(y1-l1)**2)/2 + (k2*((y2-y1)-l2)**2)/2 - m1*g*y1 - m2*g*y2
L = T - V
display(Math(r'\mathcal{L}=' + mlatex(simplify(L))))
Lterms = lagrange_terms(L, y1, y2, show=False)
display(Eq(simplify(Lterms[2]-Lterms[0]), 0))
display(Eq(simplify(Lterms[5]-Lterms[3]), F))
```

The dissipation energy of a non-conservative system with a non-conservative force (e.g., the viscous force from a damper, which is proportional to velocity) can be expressed as:

\begin{equation} D_i = \frac{1}{2}C \, \dot{q}_i^2 \label{dissipation} \end{equation}And the Lagrange's equation can be extended to include such non-conservative force in the following way:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\left( {\frac{\partial \mathcal{L}}{\partial \dot{q}_i }} \right)-\frac{\partial \mathcal{L}}{\partial q_i } + \frac{\partial D_i}{\partial \dot{q}_i }= 0 \label{lagrange_dissip} \end{equation}Consider a mass-spring-damper system with an external force acting on the mass.

The massless spring has a stiffness coefficient $k$ and length at rest $x_0$.

The massless damper has a damping coefficient $b$.

For simplicity, consider that the system starts at the resting position of the spring ($x=0$ at $x_0$).

The system has one degree of freedom and we need only one generalized coordinate ($x$) to describe the system's configuration.

Let's use the Lagrangian mechanics to derive the equations of motion for the system.

The kinetic energy of the system is:

\begin{equation} T = \frac{1}{2} m \dot x^2 \end{equation}The potential energy of the system is:

\begin{equation} V = \frac{1}{2} k x^2 \end{equation}The Lagrangian function is:

\begin{equation} \mathcal{L} = \frac{1}{2} m \dot x^2 - \frac{1}{2} k x^2 \end{equation}The dissipation energy of the system is:

\begin{equation} D = \frac{1}{2} b \dot x^2 \end{equation}Calculating all the terms in the Lagrange's equation for a dissipative process (cf. Eq. (\ref{lagrange_dissip})), the classical equation for a mass-spring-damper system can be found:

\begin{equation} m\ddot{x} + b\dot{x} + kx = F(t) \end{equation}Consider a mass-spring-damper system under the action of the gravitational force ($g=10 m/s^2$) and an external force acting on the mass.

The massless spring has a stiffness coefficient $k$ and length at rest $y_0$.

The massless damper has a damping coefficient $b$.

The gravitational force acts downwards and it is negative (see figure).
For simplicity, consider that the system starts at the resting position of the spring ($y=0$ at $y_0$).

The system has one degree of freedom and we need only one generalized coordinate ($y$) to describe the system's configuration.

Let's use the Lagrangian mechanics to derive the equations of motion for the system.

The kinetic energy of the system is:

\begin{equation} T = \frac{1}{2} m \dot y^2 \end{equation}The potential energy of the system is:

\begin{equation} V = \frac{1}{2} k y^2 + m g y \end{equation}The Lagrangian function is:

\begin{equation} \mathcal{L} = \frac{1}{2} m \dot y^2 - \frac{1}{2} k y^2 - m g y \end{equation}The dissipation energy of the system is:

\begin{equation} D = \frac{1}{2} b \dot y^2 \end{equation}The derivatives of the Lagrangian w.r.t. $y$ and $t$ are:

\begin{equation}\begin{array}{rcl} \dfrac{\partial \mathcal{L}}{\partial y} &=& -ky - mg \\ \dfrac{\partial \mathcal{L}}{\partial \dot{y}} &=& m \dot{y} \\ \dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{y}}} \right) &=& m\ddot{y} \end{array}\end{equation}The derivative of the dissipation energy w.r.t. $\dot{y}$ is:

\begin{equation} \frac{\partial D_i}{\partial \dot{y}_i } = b \dot y \end{equation}Substituting all these terms in the Lagrange's equation:

\begin{equation} \frac{\mathrm d }{\mathrm d t}\left( {\frac{\partial \mathcal{L}}{\partial \dot{q}_i }} \right)-\frac{\partial \mathcal{L}}{\partial q_i } + \frac{\partial D_i}{\partial \dot{q}_i } = Q_{NCi} \label{lagrange_dissip2} \end{equation}Results in:

\begin{equation} m\ddot{y} + b\dot{y} + ky + mg = F_0 \cos(\omega t) \end{equation}Let's solve the differential equation for the pendulum numerically using the Euler’s method.

Remember that we have to (1) transform the second-order ODE into two coupled first-order ODEs, (2) approximate the derivative of each variable by its discrete first order difference and (3) write equation to calculate the variable in a recursive way, updating its value with an equation based on the first order difference.

In Python:

In [24]:

```
import numpy as np
def msdg_euler(T, y0, v0, h):
"""
Two coupled first-order ODEs for the pendulum (v is a new variable):
dydt = v
dvdt = (F0*np.cos(omega*t) - b*v - k*y - m*g)/m
Two equations to update the values of the variables based on first-order difference:
y[n+1] = y[n] + h*v[n]
v[n+1] = v[n] + h*dvdt[n]
"""
N = int(np.ceil(T/h))
y, v = np.zeros(N), np.zeros(N)
y[0], v[0] = y0, v0
m = 1 # mass, kg
k = 100 # spring coefficient, N/m
b = 2 # damping coefficient, N/m/s
F0 = 1 # external force, N
w = 1 # angular frequency, Hz
g = 10 # acceleration of gravity, m/s2
t = np.linspace(0, T, N, endpoint=False)
F = F0*np.cos(w*t)
for i in range(N-1):
y[i+1] = y[i] + h*v[i]
v[i+1] = v[i] + h*((F[i] - b*v[i] - k*y[i] - m*g)/m)
return t, y, v
```

In [25]:

```
%matplotlib notebook
import matplotlib.pyplot as plt
def plot(t, y, v):
"""
Plot data
"""
fig, ax1 = plt.subplots(1, 1, figsize=(8, 3))
ax1.set_title('Simulation of mass-spring-damper system under gravity and external force')
ax1.plot(t, y, 'b', label='Position')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Position (m)', color='b')
ax1.tick_params('y', colors='b')
ax2 = ax1.twinx()
ax2.plot(t, v, 'r', label='Velocity')
ax2.set_ylabel('Velocity (m/s)', color='r')
ax2.tick_params('y', colors='r')
plt.tight_layout()
plt.show()
```

In [26]:

```
T, y0, v0, h = 10, .1, 0, .01
t, y, v = msdg_euler(T, y0, v0, h)
plot(t, y, v)
```

The fact the Lagrangian formalism uses generalized coordinates means that in a system with constraints we typically have fewer coordinates (in turn, fewer equations of motion) and we don't need to worry about forces of constraint that we would have to consider in the Newtonian formalism.

However, when we do want to determine a force of constraint, using Lagrangian formalism in fact will be disadvantageous! Let's see now one way of determining a force of constraint using Lagrangian formalism. The trick is to postpone the consideration that there is a constraint in the system; this will increase the number of generalized coordinates but will allow the determination of a force of constraint.
Let's exemplify this approach determining the tension at the rod in the simple pendulum under the influence of gravity we saw earlier.

Consider a pendulum with a massless rod of length $d$ and a mass $m$ at the extremity swinging in a plane forming the angle $\theta$ with vertical and $g=10 m/s^2$.

Although the pendulum moves at the plane, it only has one degree of freedom, which can be described by the angle $\theta$, the generalized coordinate. But because we want to determine the force of constraint tension at the rod, let's also consider for now the variable $r$ for the 'varying' length of the rod (instead of the constant $d$).

In this case, the kinetic energy of the system will be:

\begin{equation} T = \frac{1}{2}mr^2\dot\theta^2 + \frac{1}{2}m\dot r^2 \end{equation}And for the potential energy we will also have to consider the constraining potential, $V_r(r(t))$:

\begin{equation} V = -mgr\cos\theta + V_r(r(t)) \end{equation}The Lagrangian function is:

\begin{equation} \mathcal{L}(\theta, \dot\theta, t) = \frac{1}{2}m(\dot r^2(t) + r^2(t)\,\dot\theta^2(t)) + mgr(t)\cos\theta(t) - V_r(r(t)) \end{equation}The derivatives w.r.t. $\theta$ are:

\begin{equation} \begin{array}{rcl} &\dfrac{\partial \mathcal{L}}{\partial \theta} &=& -mgr\sin\theta \\ &\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}} &=& mr^2\dot{\theta} \\ &\dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}}} \right) &=& 2mr\dot{r}\dot{\theta} + mr^2\ddot{\theta} \end{array} \end{equation}The derivatives w.r.t. $r$ are:

\begin{equation} \begin{array}{rcl} &\dfrac{\partial \mathcal{L}}{\partial r} &=& mr \dot\theta^2 + mg\cos\theta - \dot{V}_r(r) \\ &\dfrac{\partial \mathcal{L}}{\partial \dot{r}} &=& m\dot r \\ &\dfrac{\mathrm d }{\mathrm d t}\left( {\dfrac{\partial \mathcal{L}}{\partial \dot{r}}} \right) &=& m\ddot{r} \end{array} \end{equation}The Lagrange's equations (the equations of motion) are:

\begin{equation} \begin{array}{rcl} &2mr\dot{r}\dot{\theta} + mr^2\ddot{\theta} + mgr\sin\theta &=& 0 \\ &m\ddot{r} - mr \dot\theta^2 - mg\cos\theta + \dot{V}_r(r) &=& 0 \\ \end{array} \end{equation}Now, we will apply the constraint condition, $r(t)=d$. This means that $\dot{r}=\ddot{r}=0$.

With this constraint applied, the first Lagrange's equation is the equation for the simple pendulum:

\begin{equation} md^2\ddot{\theta} + mgd\sin\theta = 0 \end{equation}The second equation yields:

\begin{equation} -\dfrac{\mathrm d V_r}{\mathrm d r}\bigg{\rvert}_{r=d} = - md \dot\theta^2 - mg\cos\theta \end{equation}But the tension force, $F_T$, is by definition equal to the gradient of the constraining potential, so:

\begin{equation} F_T = - md \dot\theta^2 - mg\cos\theta \end{equation}As expected, the tension at the rod is proportional to the centripetal and the gravitational forces.

The Lagrangian mechanics does not constitute a new theory in classical mechanics; the results of a Lagrangian or Newtonian analysis must be the same for any mechanical system, only the method used to obtain the results is different.

We are accustomed to think of mechanical systems in terms of vector quantities such as force, velocity, angular momentum, torque, etc., but in the Lagrangian formalism the equations of motion are obtained entirely in terms of the kinetic and potential energies (scalar operations) in the configuration space. Another important aspect of the force vs. energy analogy is that in situations where it is not possible to make explicit all the forces acting on the body, it is still possible to obtain expressions for the kinetic and potential energies.

In fact, the concept of force does not enter into Lagrangian mechanics. This is an important property of the method. Since energy is a scalar quantity, the Lagrangian function for a system is invariant for coordinate transformations. Therefore, it is possible to move from a certain configuration space (in which the equations of motion can be somewhat complicated) to a space that can be chosen to allow maximum simplification of the problem.

There is an elegant form to display the equations of motion using generalized coordinates $q_i$ grouping the terms proportional to common quantities in matrices, see for example, Craig (2005, page 180), Pandy (2001), and Zatsiorsky (2002, page 383):

\begin{equation} \quad M(q_i)\ddot{q_i} + C(q_i,\dot{q}_i) + G(q_i) = Q(q_i,\dot{q}_i) \label{} \end{equation}Where, for a system with $N$ generalized coordinates:

- $M$ is the inertia matrix ($NxN$);
- $\ddot{q_i}$ is a matrix ($Nx1$ of generalized accelerations);
- $C$ is a matrix ($Nx1$) of centripetal and Coriolis generalized forces;
- $G$ is a matrix ($Nx1$) of gravitational generalized forces;
- $Q$ is a matrix ($Nx1$) of external generalized forces.

The reader is invited to express the equations of motion from the examples above in this form.

- Derive the Lagrange's equation (the equation of motion) for a mass-spring system where the spring is attached to the ceiling and the mass in hanging in the vertical.
- Derive the Lagrange's equation for an inverted pendulum in the vertical.
- Derive the Lagrange's equation for the following system:
- Derive the Lagrange's equation for a spring pendulum, a simple pendulum where a mass $m$ is attached to a massless spring with spring constant $k$ and length at rest $d_0$.
- Derive the Lagrange's equation for the system shown below.
- Derive the Lagrange's equation for the following Atwood machine (consider that $m_1 > m_2$, i.e., the pulley will rotate counter-clockwise, and that moving down is in the positive direction):
- Write computer programs (in Python!) to solve numerically the equations of motion from the problems above.

- Marion JB (1970) Classical Dynamics of particles and systems, 2nd ed., Academic Press.
- Synge JL (1949) Principles of Mechanics, 2nd ed., McGraw-hill.
- Taylor J (2005) Classical Mechanics. University Science Books.