<a href="https://colab.research.google.com/github/BachiLi/A-Tour-of-Computer-Animation/blob/main/A_Tour_of_Computer_Animation_2_Lagrangian_Mechanics_and_Pendulums.ipynb" target="_parent"></a>

Last time we briefly discussed the Newtonian mechanics, where the position of each object is governed by the forces applied on them ($F=ma$). This is a very general model and covers all phenomonon in classical mechanics. However, in some scenarios, Newtonian mechanics can be a bit cumbersome for deriving the equations of motion. For example, what if our object is a pendulum affected by gravity (see the figure below)? The gravity force pulls the object downwards, but the pendulum wire applies another force to make the object goes towards the tangent direction. Deriving the pendulum wire force is a little bit annoying. Can we do better?



The trick is to reformulate Newtonian mechanics. Instead of considering the local impact of forces to an object, we consider the whole *trajectory* of the motion. Consider the following scenario: given an initial position $x_0$ at time $t_0$, and an end position at $x_1$ at time $t_1$. There are infinitely many trajectories an object can take to go from $x_0$ to $x_1$. In reality, can we predict which trajectory the object will take?



The principle for selecting the actual trajectory is called the [stationary action principle](https://en.wikipedia.org/wiki/Stationary_Action_Principle). Feynman explained this concept much better than me in his [lecture](https://www.feynmanlectures.caltech.edu/II_19.html), but I'll still try to explain it here. 

Before I start, note that the stationary action principle is asking a different problem compared to the Newtonian mechanics we introduced last time: in Newtonian mechanics, we are given an initial position and an initial velocity, and are asked to solve an [*initial value problem*](https://en.wikipedia.org/wiki/Initial_value_problem). Here we are asked to solve a [*boundary value problem*](https://en.wikipedia.org/wiki/Boundary_value_problem): given the endpoints (without the velocity), can we solve the intermediate motions? However, this point of view makes the stationary action principle a historically controversial law, philosophically: does an object somehow know where it is going in advance, then it decides the path in between? 
It turns out, mathematically, as I will show later, these two laws lead to the exact same motions in classical mechanics. So we will just use the stationary action principle as a convienent mathematical tool for deriving motions. We will leave the philosophical debate as an excersice to the readers. : )



The stationary action principle thinks in terms of energies instead of forces. The idea is that in a physical system, there is a magical quantity called energy, which is a sum of multiple components, that is kept as a constant. In classical mechanics, the energy is the sum of two components: kinetic energy $K$ and potential energy $U$. The kinetic energy represents how fast the object is moving, and the potential energy represents how much *potential* does the object have to become even faster. In classical mechanics, all motions are trade-offs between the kinetic energy and potential energy: when an object becomes faster, it is converting potential energy into kinetic energy, and vice versa. The kinetic energy $K$ is defined as:

$$K = \frac{1}{2} m \dot{x}^2$$

where $m$ is the mass of the object (so far we've assumed the mass to be 1). The potential energy $U$ is defined as the integration of force $F$:

$$\nabla_{x} U = -F$$

For example, in the case of a stone falling to the earth, the potential energy is $U = 9.8 m x$ where $x$ is the height (recall that $F = -9.8 m$).

A magical conjecture is that the sum of kinetic energy and potential energy is a constant over time: $K + U = \text{const}$, if the force $F$ is [*conservative*](https://en.wikipedia.org/wiki/Conservative_force). A conservative force means that the force is only a function of the position, and does not depend on time or velocity. 

We can verify this in the stone falling case: we know the motion is $x = -\frac{9.8}{2} t^2 + \dot{x}(0)t + x(0)$. So $K + U = -\frac{1}{2} m (-9.8t + \dot{x}(0))^2 + 9.8 m (-\frac{9.8}{2} t^2 + \dot{x}(0)t + x(0)) = -\frac{1}{2} m \dot{x}(0)^2 + 9.8 m x(0)$ (notice that all the $t$ and $t^2$ terms cancel out). 

This is actually pretty easy to prove for the general case using $F=ma$: we can look at the time derivative of $K + U$:

$$\frac{d}{dt}(K+U) = m\dot{x}\ddot{x} + \frac{dU}{dx} \dot{x} = \dot{x} \left( m \ddot{x} - F \right) = 0$$

(note that we use the conservative force assumption in the first equation: $\frac{dU}{dt} = \frac{dU}{dx}\dot{x}$ only if $U$ does not directly depend on time)

Instead of the summation $K+U$, the stationary action principle looks at the difference of them $L=K-U$. The difference is called the *Lagrangian* (probably because [Lagrange](https://en.wikipedia.org/wiki/Joseph-Louis_Lagrange) was the first person who tried to do this. By the way, the sum $H=K+U$ is called the Hamiltonian and is crucial to the [*Hamiltonian mechanics*](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) -- yet another formulation of classical mechanics. We may or may not get to Hamiltonian mechanics in the future). 

We define an *action* $S$ which is an integral over the Lagrangian $L$ over time:

$$S(x) = \int_{t=t_0}^{t_1} L(x, \dot{x}, t) = \int_{t=t_0}^{t_1} K - U$$

Note that the domain of $S$ is the set of all possible paths $x$. This makes $S$ a *functional*: it is a meta function that takes a function and outputs a real number.

What does the action $S$ mean intuitively? At each point of time, we are measuring the amount of movement $K$ plus the amount of integrated force (recall $\nabla_{x} U = -F$) applied on the object. This is why this is called an action: it is measuring how much the object is moving, plus how much force is applied. 

The stationary action principle claims that the motions in the nature will be the local extrema of the action $S(x)$ (most of the time it will be the global minimum -- nature is lazy and will minimize actions). 

To find the extrema of the action $S$, we need to first define the derivative of $S$ -- remember that $S$ is a *functional* -- what does it mean to differentiate a function whose input is a function? For the rigorous math the readers should consult [calculus of variation](https://en.wikipedia.org/wiki/Calculus_of_variations). We can roughly define the derivative $\delta S$ as the changes of $S$ when we apply a perturbation $\delta_x$ of $x$:

$$\delta S(x) = \frac{S(x + \delta_x) - S(x)}{|\delta_x|}$$

So the stationary action principle can be mathematically stated as $\delta S(x) = 0$: we want to find trajectories $x$ that satisfy this equation. Let's verify this using our stone falling to earth example again. We have $L(x, \dot{x}, t) = \frac{1}{2}m\dot{x}^2 - 9.8 m x$. We want to figure out the motion between time $t_0 = 0, t_1 = 1$, and we know $x(0)$ and $x(1)$. To solve for the motion, let's *guess* that it's going to be a qudratic function over time: $x(t) = At^2 + Bt + C$. Our goal is to solve for the parameters $A, B, $ and $C$. Our boundary constraints tells us $x(0) = C$ and $x(1) = A + B + x(0)$. The action integral is

$$S(A, B) = \int_{t=t_0}^{t_1} \frac{1}{2} m (2At+B)^2 - 9.8m(At^2+Bt+x(0))$$

With a little bit of algebra we can actually solve the action integral in closed form:

$$S(A, B) = m \left( (2A^2 - 9.8A) t^2 + (2AB - 9.8B) t + \frac{1}{2}B^2 - 9.8 x(0) \right)$$

To solve for the extrema of $S(A, B)$ under our constraints $A+B = x(1) - x(0)$, we use a trick called the [Lagrange multiplier](https://en.wikipedia.org/wiki/Lagrange_multiplier) (yes it's the same Lagrange): we solve for the extrema of $G(A, B, \lambda) = S(A, B) + \lambda (A+B+x(0)-x(1))$, where $\lambda$ is a free variable. Setting the derivative of $G$ to zero we get three equations with three unknowns ($A, B, \lambda$). Solving for the system we get $A = -\frac{9.8}{2}$ and $B=x(1)-x(0)-\frac{9.8}{2}$.

Compare the motion equation $x(t) = -\frac{9.8}{2}t^2 + (x(1)-x(0)-\frac{9.8}{2})t + x(0)$ to the equation from our last chapter $x(t) -\frac{9.8}{2}t^2 + \dot{x}(0)t + x(0)$, we know that if we set our initial velocity $\dot{x}(0)$ as $x(1)-x(0)-\frac{9.8}{2}$, we will arrive at $x(1)$ at time $1$.

**Exercise**: we can state a different version of stationary action principle assuming we only know $x(0)$ and $\dot{x}(0)$, instead of $x(0)$ and $x(1)$. Does this different version predict the same motion?

In practice, we don't always want to assume a parametric form of our motion trajectory. What can we do in that case?

It turns out that the least action principle is a very strong constraint: even though we are imposing a constraint on an integral over some time, it ends up constraining each point in the time (intuitively, since $x$ is an infinite-dimensional function, taking its derivative and setting it to zero gives us infinitely many constraints).

The constraint is called the [Euler-Lagrange equation](https://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation) (Euler proposed a first version of this, and Lagrange refined it to the current form), and can be written as follows:

$$\frac{\partial L}{\partial x} = \frac{d}{dt} \frac{\partial L}{\partial \dot{x}}$$

For example, if $L(x, \dot{x}, t) = \frac{1}{2}m\dot{x}^2 - 9.8 m x$, then the corresponding Euler-Langrange equation is

$$
\begin{split}
\frac{\partial{L}}{\partial x} &= -9.8m \\
\frac{\partial L}{\partial \dot{x}} &= m\dot{x} \\
\frac{d}{dt} \frac{\partial L}{\partial \dot{x}} &= m\ddot{x} \\
\frac{\partial L}{\partial x} &= \frac{d}{dt} \frac{\partial L}{\partial \dot{x}} \rightarrow \ddot{x} = -9.8 
\end{split}
$$

We recover the equation of motion from Newtonian mechanics we had from our last chapter! In fact, the Euler-Lagrange equation is usually just a different way of writing $F=ma$. If we look at it closely, since the kinetic energy is always $\frac{1}{2}m\dot{x}^2$, $\frac{\partial L}{\partial x} = -\frac{\partial U}{\partial x} = F$, and if the potential energy $U$ does not depend on the velocity $\dot{x}$, $\frac{d}{dt} \frac{\partial L}{\partial \dot{x}} = m \ddot{x}$.

The derivation of the Euler-Lagrange equation from the stationary action principle can be found online: I recommend the proof from [Preetum Nakkiran](https://preetum.nakkiran.org/lagrange.html) which is geometric and very intuitive. [Wikipedia](https://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation#Statement) also has the standard proof using integration by parts. I also like the derivation from [Jason Gross](https://web.mit.edu/~jgross/Public/least_action/Principle%20of%20Least%20Action%20with%20Derivation.html).

The Euler-Lagrange equation connects the stationary action principle and Newtonian mechanics. Remember we set up the problem formulation as a *boundary value problem*. However, Newtonian mechanics predicts an *initial value problem*. Interestingly, since the stationary action principle puts a constraint on every single time instance, **no matter where the particle is going eventually, they need to follow the stationary action principle and Euler-Lagrange equation all the time**. 

So, what have we gained over all this crazy math? Let's go back to the example at the beginning, where we have a pendulum. The trick is to use something called the *generalized coordinates*. Instead of specifying the Cartesian position $x$ of the pendulum in 2D, we use a coordinate $q$ based on the angle of the pendulum. We have the relation $x(q) = (r \sin(q), r \cos(q))$ where $r$ is the length of the wire.




Recall that, earlier, we parametrize the trajectory of a stone falling to earth using a quadratic $At^2+Bt+C$, and the sationary principle helps us solve the quadratic parameters. We are doing to same thing here: we "parametrize" the trajector $x$ using the angle $q$, and solve for q instead. Our Lagrangian is now:

$$L(q, \dot{q}, t) = \frac{1}{2} m \dot{x}(q)^2 - 9.8 m (r - x(q)) = \frac{1}{2} m r^2 \dot{q}^2 - 9.8m r (1 - \cos(q))$$

(since $x$ is a vector, $\dot{x}(q) = -r \sin(q) \dot{q}, r \cos(q) \dot{q}$, $\dot{x}(q)^2 = \dot{x}(q) \cdot \dot{x}(q) = r^2 \dot{q}^2$)

The cool thing is that the Euler-Lagrange equation still holds:

$$
\begin{split}
\frac{\partial L}{\partial q} &= \frac{d}{dt} \frac{\partial L}{\partial \dot{q}} \\
\frac{\partial L}{\partial q} &= - 9.8mr \sin(q) \\
\frac{\partial L}{\partial \dot{q}} &= m r^2 \dot{q} \\
\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} &= mr^2 \ddot{q} \\
- 9.8mr \sin(q) &= mr^2 \ddot{q} \\
&\downarrow \\
\ddot{q} &= \frac{-9.8}{r} \sin(q) \\
\end{split}
$$

We get a simple second-order ODE! In fact you can derive the same thing using Newtonian mechanics (see [here](https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendula.html)), however it requires you to factor out the forces and constraints -- something that demands a bit of cleverness. Lagrangian mechanics allows us to do things in a very mechanical, automatic way: it's all basic algebra.

Interestingly, despite its simplicity, this ODE does not have a known closed form. Let's try to visualize it using forward Euler method (finally some code!!):


```python
import math
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML

fps = 20

q0 = math.pi / 4
v0 = 0.0
r = 20.0
ts = 0.01
def pendulum_forward_euler(t, q, v):
    ct = 0
    while True:
        nq = q + ts * v
        nv = v - ts * (9.8 / r) * math.sin(q)
        q = nq
        v = nv
        ct += ts
        if ct >= t:
            break
    return ct, q, v

def visualize(ode_solver):
    fig = plt.figure(figsize=(8, 4))
    ax = plt.axes(xlim=(-25, 25), ylim=(-25, 0))
    line, = ax.plot([], [], '-', lw=2)
    point, = ax.plot([], [], 'g.', ms=20)

    t = 0
    q = q0
    v = v0
    def animate(i):
        nonlocal t, q, v
        dt, q, v = ode_solver(1.0/fps, q, v)
        x = r * math.sin(q)
        y = -r * math.cos(q)
        t += dt
        line.set_data([0, x], [0, y])
        point.set_data([x], [y])
        return point,

    plt.close()
    return animation.FuncAnimation(fig, animate, frames=400, interval=fps, blit=True)

anim = visualize(pendulum_forward_euler)
HTML(anim.to_html5_video())
```








**Exercise**: can you modify the pendulum to handle the collision with the ceiling?

We'll close this chapter by deriving the Lagrangian mechanics for a double pendulum. The derivation is pretty messy, so we will use sympy to help us deriving the equations.

The first mass on the pendulum is parametrized by $x_0 = (r_0 \sin(q_0), -r_0\cos(q_0))$, the second mass is parametrized by $x_1 = x_0 + (r_1 \sin(q_1), -r_1 \cos(q_1))$. 



Our Lagrangian is:

$$L(q,\dot{q},t) = \frac{1}{2} m_0 r_0^2 \dot{q}_0^2 + \frac{1}{2} m_1 \left( r_0^2 \dot{q}_0^2 + 2 r_0 r_1 \dot{q}_0 \dot{q}_1 \cos(q_0 - q_1) + r_1^2 \dot{q}_1^2 \right) - g \left(m_0 r_0 (- \cos(q_0)) + m_1 \left( r_0 (- \cos(q_0)) + r_1 (- \cos(q_1)) \right) \right)$$

(we finally replace the 9.8 magic constant with g)

Let's define it in sympy:


```python
import sympy

sympy.init_printing() # printing latex

t, m0, r0, m1, r1, g = sympy.symbols('t, m0, r0, m1, r1, g')
q0 = sympy.Function('q0')(t)
q1 = sympy.Function('q1')(t)
v0 = q0.diff(t) # \dot{q}_0
v1 = q1.diff(t) # \dot{q}_1
half = sympy.Rational(1, 2)

# Lagrangian
K = half * m0 * r0 * r0 * v0 * v0 + \
     half * m1 * (r0 * r0 * v0 * v0 + 2 * r0 * r1 * v0 * v1 * sympy.cos(q0 - q1) + r1 * r1 * v1 * v1)
U = g * (m0 * r0 * (- sympy.cos(q0)) + m1 * (r0 * (- sympy.cos(q0)) + r1 * (- sympy.cos(q1))))
L = sympy.simplify(K - U)

L
```

Next we use sympy to derive the Euler-Lagrange equation (sympy has a function [euler_equations](https://docs.sympy.org/latest/modules/calculus/index.html) for this purpose, but we'll do it manually just for this time):


```python
dLdq0 = L.diff(q0).simplify()
dLdq1 = L.diff(q1).simplify()
dLdv0 = L.diff(v0).simplify()
dLdv1 = L.diff(v1).simplify()
dLdv0dt = dLdv0.diff(t).simplify()
dLdv1dt = dLdv1.diff(t).simplify()

eq0 = sympy.simplify(dLdq0 - dLdv0dt) # this should be 0
eq1 = sympy.simplify(dLdq1 - dLdv1dt) # this should be 0
```


```python
eq0
```

sympy's equation solver didn't do what I wanted. So we'll have to manually arrange the equations.


```python
# simplify eq0 manually
eq0 = -eq0 / r0
eq0.args
```


```python
eq0 = sympy.simplify(eq0.args[0] + eq0.args[1]) + sympy.simplify(eq0.args[2] + eq0.args[3]) + eq0.args[4] + eq0.args[5]
eq0.args
```


```python
eq1
```


```python
# simplify eq1 manually
eq1 = -eq1 / (m1 * r1)
eq1.args
```


```python
# eliminate a1
a0 = (eq0 - eq1 * m1 * sympy.cos(q0 - q1)).expand()
a0.args
```


```python
# manually simplify a0
a0 = sympy.simplify(a0.args[0] + a0.args[1]) + sympy.simplify(a0.args[2] + a0.args[3] + a0.args[6]) + a0.args[4] + a0.args[5] + a0.args[7]
a0.args
```


```python
# r0 (m0 - m1 cos^2(q0 - q1) + m1) = r0 (m0 + m1 sin^2(q0 - q1))
s = sympy.sin(q0 - q1)
# minus sign for moving d^2 dq0/dt^2 to the other side of the equation
a0solve = -(a0.args[0] + sum(a0.args[2:])) / (r0 * (m0 + m1 * s * s))
a0solve
```


```python
# eliminate a0
a1 = (eq0 * sympy.cos(q0 - q1) - eq1 * (m0 + m1)).expand()
a1.args
```


```python
# manually simplify a1
a1 = sympy.simplify(a1.args[0] + a1.args[1]) + sympy.simplify(a1.args[2] + a1.args[3] + a1.args[8]) + sympy.simplify(a1.args[4] + a1.args[5]) + sympy.simplify(a1.args[6] + a1.args[7]) + a1.args[9]
a1.args
```


```python
# r1 (-m0 + m1 cos^2(q0 - q1) - m1) = -r1 (m0 + m1 sin^2(q0 - q1))
# a minus sign for the denominator, another minus sign for moving d^2 dq1/dt^2 to the other side of the equation
a1solve = sum(a1.args[1:]) / (r1 * (m0 + m1 * s * s))
a1solve
```


```python
sympy.Eq(q0.diff(t).diff(t), a0solve)
```


```python
sympy.Eq(q1.diff(t).diff(t), a1solve)
```

Above is our velocity update formula. Let's try to implement this!


```python
import math
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
fps = 20
g = 9.8
q0_init = math.pi
v0_init = 0.0
q1_init = -math.pi / 4
v1_init = 0.0
r0 = 20.0
r1 = 20.0
m0 = 1
m1 = 1
ts = 0.001
def double_pendulum_forward_euler(t, q0, v0, q1, v1):
    ct = 0
    while True:
        nq0 = q0 + ts * v0
        nq1 = q1 + ts * v1
        s0 = math.sin(q0)
        s1 = math.sin(q1)
        sd = math.sin(q0 - q1)
        cd = math.cos(q0 - q1)
        nv0 = v0 + ts * (g * m1 * s1 * cd - g * (m0 + m1) * s0 - m1 * r0 * sd * cd * v0 * v0 - m1 * r1 * sd * v1 * v1) / (r0 * (m0 + m1 * sd * sd))
        nv1 = v1 + ts * (g * (m0 + m1) * s0 * cd - g * (m0 + m1) * s1 + m1 * r1 * sd * cd * v1 * v1 + r0 * (m0 + m1) * sd * v0 * v0) / (r1 * (m0 + m1 * sd * sd))

        q0 = nq0
        q1 = nq1
        v0 = nv0
        v1 = nv1
        ct += ts
        if ct >= t:
            break
    return ct, q0, v0, q1, v1

def visualize(ode_solver):
    fig = plt.figure(figsize=(8, 8))
    ax = plt.axes(xlim=(-55, 55), ylim=(-55, 55))
    line, = ax.plot([], [], '-', lw=2)
    point, = ax.plot([], [], 'g.', ms=20)

    t = 0
    q0 = q0_init
    v0 = v0_init
    q1 = q1_init
    v1 = v1_init
    def animate(i):
        nonlocal t, q0, v0, q1, v1
        dt, q0, v0, q1, v1 = ode_solver(1.0/fps, q0, v0, q1, v1)
        x0 = r0 * math.sin(q0)
        y0 = -r0 * math.cos(q0)
        x1 = x0 + r1 * math.sin(q1)
        y1 = y0 - r1 * math.cos(q1)
        t += dt
        line.set_data([0, x0, x1], [0, y0, y1])
        point.set_data([x0, x1], [y0, y1])
        return point,

    plt.close()
    return animation.FuncAnimation(fig, animate, frames=800, interval=fps, blit=True)

anim = visualize(double_pendulum_forward_euler)
HTML(anim.to_html5_video())
```








**Exercise**: can you implement the collision between the masses and the wires?

This actually already allows us to do some fun character animations! Imagine the double pendulum above being an arm of a character. You can setup an articulated character and derive its dynamics yourself using the Lagrangian mechanics.

In summary, we have learned a very powerful way, called Lagrangian mechanics, to derive equations of motion. The usual process to set up a computer animation, is to specify our animation in a generalized coordinates $q$ and their mapping to the Cartesian coordinates $x$. We then derive our Lagrangian $L$: the kinetic energy $K$ is always $\frac{1}{2} m v^2$, the potential energy $U$ depends on the interactions you want to have between the objects and user inputs. We then solve the sationary action principle or the Euler-Lagrange equation to compute the animation at certain time. This process is very mechanical and mostly involving just doing some algebra, and importantly, it enables us to design physically-based animation with great flexibility!

In the next chapter, we will analyze the discretization we've been doing to the differential equation: Are there better way to do it other than the forward Euler method? (yes) What are the different desirable criteria for a discretization?

b980620e727eaf80f8a7e48c5b684559afdab416
ipynb
Jupyter Notebook
A_Tour_of_Computer_Animation_2_Lagrangian_Mechanics_and_Pendulums.ipynb
BachiLi/A-Tour-of-Computer-Animation
c739e333a5519b1ec379b8a867536de65233d0ea
MIT
2021-06-15T16:07:28.000Z
2021-08-12T02:01:30.000Z
A_Tour_of_Computer_Animation_2_Lagrangian_Mechanics_and_Pendulums.ipynb
BachiLi/A-Tour-of-Computer-Animation
c739e333a5519b1ec379b8a867536de65233d0ea
MIT
A_Tour_of_Computer_Animation_2_Lagrangian_Mechanics_and_Pendulums.ipynb
BachiLi/A-Tour-of-Computer-Animation
c739e333a5519b1ec379b8a867536de65233d0ea
MIT
Qwen/Qwen-72B
1. YES
2. YES