Calculus of Variations

General Idea

In calculus, we take derivatives of functions of a few variables. In calculus of variations, we take derivatives of functions of functions (also called functionals). The main purpose of this is to formulate the optimization problem over function spaces and derive its optimality condition (KKT Condition).

These optimality conditions are often differential equations, called Euler-Lagrange Equations, and most physical equations arise as Euler-Lagrange equations of some optimization problem. That way, we can “design one optimal functional” instead of “modeling forces”.

Example 1

Let

$$M=\{y:[0,T]\to \mathbb{R}\mid y(0)=y_0,y(T)=y_T\}$$

be the manifold of all possible valid trajectories of a particle in one dimension, from $y_0$ to $y_T$ in time $T$. A single “point” on this manifold is an entire function $y(t)$. Because it takes infinitely many values (one for each $t$), we can think of $M$ as an infinite-dimensional manifold.

Let $\mathcal{E}:M\to \mathbb{R}$ be the covector where

$$\mathcal{E}(y):={\int }_{t=0}^T\left(\frac{1}{2}{y}^{\prime }(t{)}^2+\cos (y(t))\right)dt$$

The problem is, what is the optimality condition for

$$\underset{y\in M}{\min }\mathcal{E}(y)$$

If we fixed some $p\in M$, the tangent space $T_{y}M$ represents all valid perturbations we can apply to $y$ without leaving the Manifold. This perturbation is $\mathring{y}$. So, for perturbed path $y(t)+\epsilon \mathring{y}(t)\in M$, it must satisfy the fixed boundary conditions. In particular,

$$T_{y}M=\{\mathring{y}:[0,T]\to \mathbb{R}\mid \mathring{y}(0)=0,\mathring{y}(T)=0\}$$

Visually, we have

where the blue curve at the bottom represents one example of a perturbation $\mathring{y}$, and the purple curves represent some examples of $y$.

The optimality condition is

$$d{\mathcal{E}}_y[[\mathring{y}]]=0$$

for all $\mathring{y}\in T_{y}M$, just like the first derivative test for optimality in regular calculus. Geometrically, this means the covector itself is zero. I.e. if we were at the bottom of some “energy bowl”, no matter which direction $\mathring{y}$ we perturb, the energy will not change. Indeed,

$$\begin{array}{rl}0& =d{\mathcal{E}}_y[[\mathring{y}]]\\ & =\frac{d}{d\epsilon }{|}_{\epsilon =0}\mathcal{E}(y+\epsilon \mathring{y})\\ & =\frac{d}{d\epsilon }{|}_{\epsilon =0}{\int }_{t=0}^T\left(\frac{1}{2}(y^{\prime }(t)+\epsilon {\mathring{y}}^{\prime }(t){)}^2+\cos (y(t)+\epsilon \mathring{y}(t))\right)dt\\ & ={\int }_{t=0}^T\left(y^{\prime }{\mathring{y}}^{\prime }-\sin (y)\mathring{y}\right)dt\\ & =-{\int }_{t=0}^T\left(y^{\prime \prime }+\sin (y)\right)\mathring{y}dt\text{via integration by parts}\\ & ⟹0=y^{\prime \prime }(t)+\sin (y(t))\end{array}$$

which is precisely the physical equation of motion for a pendulum.

Physics is just optimization over function spaces.

Euler-Lagrange Equation

We can derive the optimality condition for a more general functional. Suppose we have the same conditions as above but with

$$\mathcal{E}(y):={\int }_{t=0}^{T}L(t,y(t),y^{\prime }(t))dt$$

for some integrand $L:[0,T]\times {\mathbb{R}}^m\times {\mathbb{R}}^m\to \mathbb{R}$. The directional derivative is then

$$\begin{array}{rl}d{\mathcal{E}}_{\mathbf{y}}[[\mathring{\mathbf{y}}]]& ={\frac{d}{d\epsilon }|}_{\epsilon =0}\mathcal{E}(\mathbf{y}+\epsilon \mathring{\mathbf{y}})\\ & ={\int }_0^T{\frac{\partial }{\partial \epsilon }|}_{\epsilon =0}L(t,\mathbf{y}(t)+\epsilon \mathring{\mathbf{y}}(t),{\mathbf{y}}^{\prime }(t)+\epsilon {\mathring{\mathbf{y}}}^{\prime }(t))dt\\ & ={\int }_0^T⟨{\frac{\partial L}{\partial \mathbf{y}}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}|\mathring{\mathbf{y}}(t)⟩dt+{\int }_0^T⟨{\frac{\partial L}{\partial {\mathbf{y}}^{\prime }}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}|{\mathring{\mathbf{y}}}^{\prime }(t)⟩dt\\ & ={\int }_0^T⟨{\frac{\partial L}{\partial \mathbf{y}}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}|\mathring{\mathbf{y}}(t)⟩dt-{\int }_0^T⟨\frac{d}{dt}\left({\frac{\partial L}{\partial {\mathbf{y}}^{\prime }}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}\right)|\mathring{\mathbf{y}}(t)⟩dt\\ & ={\int }_0^T⟨{\frac{\partial L}{\partial \mathbf{y}}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}-\frac{d}{dt}\left({\frac{\partial L}{\partial {\mathbf{y}}^{\prime }}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}\right)|\mathring{\mathbf{y}}(t)⟩dt\end{array}$$

in particular,

1. We expanded the definition of directional derivative
2. We applied the multivariate chain rule since $\epsilon$ appears in both $\mathbf{y}$ and ${\mathbf{y}}^{\prime }$
3. We took the partial derivative of $L$ with respect to $\epsilon$. The use of $⟨,⟩$ is just to denote the dot product between the gradient and the perturbation. Here, it is the inner product in ${\mathbb{R}}^m$.
4. We applied integration by parts to move the derivative from ${\mathring{\mathbf{y}}}^{\prime }$ to $\frac{\partial L}{\partial {\mathbf{y}}^{\prime }}$. The boundary term vanishes since $\mathring{\mathbf{y}}(0)=\mathring{\mathbf{y}}(T)=0$.
5. We factored out $\mathring{\mathbf{y}}$ via linearity of the inner product.

Theorem (Euler-Lagrange Equation)

The optimality condition $d{\mathcal{E}}_{\mathbf{y}}[[\mathring{\mathbf{y}}]]=0$ for all $\mathring{\mathbf{y}}\in T_{\mathbf{y}}M$ is equivalent to the following system of ordinary differential equations:

$$\frac{\partial \mathcal{E}}{\partial \mathbf{y}}(t)={\frac{\partial L}{\partial \mathbf{y}}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}-\frac{d}{dt}\left({\frac{\partial L}{\partial {\mathbf{y}}^{\prime }}}_{(t,\mathbf{y}(t),{\mathbf{y}}^{\prime }(t))}\right)=0$$

called the Euler-Lagrange Equation.

Applications in Physics

The Euler-Lagrange equations are foundational to physics. For detailed derivations of Newton’s laws, conservation of energy, and rotational dynamics from variational principles, see Lagrangian Mechanics.

Example 2 (Hanging Chain)

What is the shape of a hanging chain on two poles? It is the curve with the lowest potential energy with fixed total length. Let the shape be the function graph

$$y=f(x)$$

First, we want to describe the physical geometry of the curve. Any section of the chain forms a right triangle microscopically. With the Pythagorean theorem, we have that the length of the chain

$$ds=\sqrt{d{x}^2+d{y}^2}=\sqrt{1+{\left(\frac{dy}{dx}\right)}^2}dx=\sqrt{1+f^{\prime 2}}dx$$

The total length is found by integrating $ds$ along the curve

$$\mathcal{L}(f)={\int }_a^{b}ds={\int }_a^b\sqrt{1+f^{\prime 2}}dx$$

The total gravitational potential energy is

$$\mathcal{E}(f)={\int }_a^{b}f\sqrt{1+f^{\prime 2}}dx$$

We must minimize $\mathcal{E}$. However, without any constraints, the chain would drop straight down into infinity (to achieve negative infinite potential energy). So we must constrain it to have a fixed length $\ell$. Thus, the constraint is $\mathcal{L}(f)-\ell =0$. The KKT Condition for minimizing $\mathcal{E}$ with constraint $\mathcal{L}-\ell =0$ says that there exists some $\lambda$ such that

$$d\mathcal{E}+\lambda d\mathcal{L}=0⟺d\tilde{\mathcal{E}}=0\text{ where }\tilde{\mathcal{E}}=\mathcal{E}+\lambda \mathcal{L}={\int }_a^b(f+\lambda )\sqrt{1+f^{\prime 2}}dx$$

We want to minimize $\tilde{\mathcal{E}}$. By applying Noether’s theorem of time independence (in this case, “time” is just the $x$-axis), we have that

$$\frac{\partial L}{\partial f^{\prime }}{f}^{\prime }-L=c$$

implying that

$$(f+\lambda )\left(\frac{f^{\prime 2}}{\sqrt{1+f^{\prime 2}}}-\sqrt{1+f^{\prime 2}}\right)=c$$

such that

$$\begin{array}{rl}{f}^{\prime }& =\pm \frac{1}{c}\sqrt{(f+\lambda {)}^2-c^2}\\ \frac{df}{\sqrt{(f+\lambda {)}^2-c^2}}& =\pm \frac{1}{c}dx\\ \int \frac{df}{\sqrt{(f+\lambda {)}^2-c^2}}& =\pm \int \frac{1}{c}dx\\ {\cosh }^{-1}\left(\frac{f+\lambda }{c}\right)& =\pm \frac{x}{c}+c_1\\ f(x)& =c\cosh \left(\frac{x}{c}+c_2\right)-\lambda \end{array}$$

Example 3 (Mass Spring System)

Consider a graph $G=(V,E)$ where each vertex is assigned a mass $m_i$ and position $x_i$. Each edge is assigned a spring with a spring constant $k_e$. The total kinetic energy and potential energy is

$$K=\sum _{i\in V}\frac{1}{2}{m}_i|{\dot{x}}_i{|}^{2}U=\sum _{e=(i,j)\in E}\frac{1}{2}{k}_e|x_i-x_j{|}^2$$

The Euler-Lagrange equation for this gives

$$\frac{d}{dt}\frac{\partial L}{\partial {\dot{x}}_i}=-\frac{\partial L}{\partial x_i}⟹m_i{\ddot{x}}_i=-\frac{\partial U}{\partial x_i}$$

We can find the derivative with respect to one vertex position

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1.0]\n\n$

source code

One interesting thing to note is that this is similar to a feed forward neural network. We input raw vertex positions and they transform step by step into a single scalar output $U$.

1. The first transformation layer subtracts the position of the source vertex from the position of the destination vertex to get the edge vector.
2. The second transformation layer normalizes (calculates the magnitude) the edge vector to get the current physical length of the string.
3. The third transformation layer computes the local energy of each edge via the formula $\frac{1}{2}{k}_e{\ell }_e^2$.
4. The sum of the individual spring energies gives us the total potential energy of the system.

Now, like in machine learning, we do a backward pass. To find the force on a specific vertex $i$, we must compute $\frac{\partial U}{\partial {\mathbf{x}}_i}$. Because the energy was calculated through a chain of operations, we must use the multivariate chain rule to backpropagate the gradient through the computational graph.

On each edge $e$ of the graph, we can compute ${\sigma }_e=k_e(d{\mathbf{x}}_e)$, which is the stress (stiffness times strain) on that edge. To get the actual force ${\mathbf{f}}_i$ actuing on each node, the chain rule requires us to sum up the stress vectors from all the edges that touch that specific node. Indeed,

$$f_i=\text{div}(\sigma )$$

where $\text{div}$ is the graph divergence operator. Thus, the force on a node is the aggregate of the stress from all edges touching that node.

Remark

But of course, since this is just like machine learning, we can efficiently compute the forces on all nodes by matrix multiplications. Truly, we have

$$\left[\begin{array}{ccccc}{m}_1& & & & \\ & \ddots & & & \\ & & m_i& & \\ & & & \ddots & \\ & & & & m_{|\mathcal{V}|}\end{array}\right]\left[\begin{array}{c}{\ddot{\mathbf{x}}}_1\\ ⋮\\ {\ddot{\mathbf{x}}}_i\\ ⋮\\ {\ddot{\mathbf{x}}}_{|\mathcal{V}|}\end{array}\right]=\underset{\text{Graph Laplacian}}{\underbrace{-{\mathbf{d}}^{\top }\left[\begin{array}{ccccc}{k}_1& & & & \\ & \ddots & & & \\ & & k_e& & \\ & & & \ddots & \\ & & & & k_{|\mathcal{E}|}\end{array}\right]\mathbf{d}}}\left[\begin{array}{c}{\mathbf{x}}_1\\ ⋮\\ {\mathbf{x}}_i\\ ⋮\\ {\mathbf{x}}_{|\mathcal{V}|}\end{array}\right]$$

where the left side is global Newton’s Second Law $F=ma$ and the right side is the restorative spring forces. Note that $\mathbf{K}(\mathbf{dx})$ here is Hooke’s Law (forces on linear springs). The ${\mathbf{d}}^{\top }$ is actually the adjoint linear map (recall that the transpose is the adjoint of a linear map with respect to the standard inner product) acting as the discrete gradient. The multiplication of ${\mathbf{d}}^{\top }(\mathbf{K}\mathbf{dx})$ sums up all the converging spring tensions onto their shared nodes, calculating the net force on each node.

It’s particularly important to mention the Graph Laplacian shown in the diagram. The final system becomes

$$\mathbf{M}\ddot{\mathbf{x}}=-\mathbf{Lx}$$

where $\mathbf{L}$ encodes the connectivity of the graph and the stiffness of each edge.

This is a beautiful result.