The calculus of variations and the Euler–Lagrange equation

Ordinary calculus finds the points where a function has a minimum, a maximum or a flat spot by setting its derivative to zero. The calculus of variations asks the same kind of question one level up: among all functions \(y(x)\) satisfying some end conditions, which one makes a given integral stationary? The answer is not a number but a function, and the condition it must satisfy turns out to be a differential equation, the Euler–Lagrange equation. This lesson derives that equation, explores its special cases and first integrals, and applies it to optics, geometry and mechanics.

Euler developed the subject in his 1744 book on curves with maximum or minimum properties, using a geometric approach in which a curve was replaced by a polygon. Lagrange, then nineteen, wrote to Euler in 1755 with a cleaner method based on varying the whole function at once, which he called the method of variations. Euler adopted Lagrange's method and gave the field its name. The derivation below is Lagrange's.

Functionals

A functional takes a function and returns a number. The most important kind for us has the form

where \(F\) is a given function of three variables and \(y(x)\) is any smooth function with fixed end values \(y(a) = y_a\), \(y(b) = y_b\). Some examples: the length of a curve in the plane is \(\int \sqrt{1 + y'^2}\,dx\); the time for light to cross a medium with refractive index \(n(x, y)\) is \(\frac1c\int n\sqrt{1 + y'^2}\,dx\); and the action of lesson 1 is \(\int L(x, \dot x, t)\,dt\), the same form with \(t\) in place of \(x\).

Deriving the Euler–Lagrange equation

Suppose \(y(x)\) makes \(J\) stationary. Compare it with a nearby function \(y(x) + \varepsilon\eta(x)\), where \(\varepsilon\) is a small number and \(\eta(x)\) is any smooth function that vanishes at both ends, \(\eta(a) = \eta(b) = 0\), so that the varied function satisfies the same end conditions. Then \(J[y + \varepsilon\eta]\) is an ordinary function of the single number \(\varepsilon\), and stationarity means its derivative at \(\varepsilon = 0\) vanishes for every choice of \(\eta\).

Differentiating under the integral sign, using the chain rule on \(F\),

Integrate the second term by parts. The boundary term \(\big[\frac{\partial F}{\partial y'}\eta\big]_a^b\) vanishes because \(\eta\) is zero at both ends, leaving

This must hold for every smooth \(\eta\) that vanishes at the ends. If the bracket were nonzero at some point, say positive, we could choose \(\eta\) to be a small positive bump around that point and zero elsewhere, making the integral positive. So the bracket must vanish everywhere. This step is the fundamental lemma of the calculus of variations, and the result is the Euler–Lagrange equation:

It is a second-order differential equation for \(y(x)\), because \(\frac{d}{dx}\) acting on \(\partial F/\partial y'\) brings in \(y''\). Its two constants of integration are fixed by the two end conditions. A solution of the Euler–Lagrange equation is called an extremal; whether it is a minimum, a maximum or a saddle needs further analysis involving the second variation.

The simplest examples

Shortest path in the plane. With \(F = \sqrt{1 + y'^2}\), \(\partial F/\partial y = 0\), so the equation says \(\frac{d}{dx}\frac{y'}{\sqrt{1 + y'^2}} = 0\). Hence \(y'/\sqrt{1 + y'^2}\) is constant, so \(y'\) is constant and the extremal is a straight line. The calculus of variations confirms what everyone believes, and the same calculation on a curved surface will give its geodesics in lesson 3.

Mechanics. With \(t\) in place of \(x\), \(q\) in place of \(y\) and \(F = L(q, \dot q, t)\), the equation becomes \(\frac{d}{dt}\frac{\partial L}{\partial\dot q} - \frac{\partial L}{\partial q} = 0\), Lagrange's equation. With several coordinates \(q_1, \dots, q_n\) we vary each one independently, and we get one Euler–Lagrange equation for each.

A quadratic functional. Take \(J = \int_0^1 (y'^2 + y^2)\,dx\) with \(y(0) = 0\), \(y(1) = 1\). Then \(\partial F/\partial y' = 2y'\) and \(\partial F/\partial y = 2y\), so \(y'' = y\), with solution \(y = \sinh x/\sinh 1\). Substituting back, the minimum value is \(J = \coth 1 = 1.3130\). The program below finds this numerically by discretizing the integral and relaxing the interior values, just as lesson 1 did for the action.

First integrals

Two special cases of the Euler–Lagrange equation can be integrated once immediately, and both appear constantly.

Case 1: \(F\) does not depend on \(y\). Then \(\partial F/\partial y = 0\) and the equation says \(\partial F/\partial y'\) is constant. In mechanics, a coordinate that does not appear in the Lagrangian is called cyclic or ignorable, and its conjugate momentum \(p = \partial L/\partial\dot q\) is conserved. The straight-line example above was of this kind.

Case 2: \(F\) does not depend explicitly on \(x\). Then the Beltrami identity holds:

To prove it, differentiate the left side with respect to \(x\) along a solution. Using the chain rule, \(\frac{dF}{dx} = \frac{\partial F}{\partial y}y' + \frac{\partial F}{\partial y'}y''\) (no \(\partial F/\partial x\) term). The derivative of the second term is \(y''\frac{\partial F}{\partial y'} + y'\frac{d}{dx}\frac{\partial F}{\partial y'}\). Subtracting, the \(y''\) terms cancel and what remains is \(y'\big(\frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'}\big) = 0\) by the Euler–Lagrange equation. In mechanics, with \(F = L\) and \(x = t\), the Beltrami identity says that \(\dot q\,\partial L/\partial\dot q - L\) is conserved when \(L\) has no explicit time dependence. For \(L = \tfrac12 m\dot q^2 - V(q)\) this quantity is \(\tfrac12 m\dot q^2 + V\): the energy. Conservation of energy is the Beltrami identity applied to the action.

Fermat's principle and Snell's law

Light travels at speed \(c/n\) in a medium of refractive index \(n\). Fermat's principle of 1662 states that the path taken by light between two points makes the travel time \(\frac1c\int n\,ds\) stationary. In a uniform medium the path is a straight line. When light crosses a flat boundary between media with indices \(n_1\) and \(n_2\), the path consists of two straight segments, and only the crossing point \(x\) on the boundary is free. With source at height \(h_1\) above the surface and target at depth \(h_2\) below, horizontal separation \(D\), the time is

Setting \(dT/dx = 0\) gives \(n_1\frac{x}{\sqrt{x^2 + h_1^2}} = n_2\frac{D - x}{\sqrt{(D - x)^2 + h_2^2}}\), which is \(n_1\sin\theta_1 = n_2\sin\theta_2\) with angles measured from the normal: Snell's law, found by Ibn Sahl around 984 and by Willebrord Snellius in 1621. The program checks this by searching numerically for the quickest crossing point.

When \(n\) varies continuously with height, \(n = n(y)\), the travel time has the form \(\int n(y)\sqrt{1 + y'^2}\,dx\), with no explicit \(x\). The Beltrami identity gives \(n(y)/\sqrt{1 + y'^2} = \text{constant}\), which is \(n\cos\alpha = \text{constant}\) with \(\alpha\) the angle of the ray to the horizontal: Snell's law applied layer by layer. This explains mirages, the bending of starlight near the horizon, and the way sound bends in the atmosphere on a cold night.

import math

# 1. Minimize J[y] = integral_0^1 (y'^2 + y^2) dx with y(0) = 0, y(1) = 1.
#    Euler-Lagrange: y'' = y, so y = sinh(x)/sinh(1).
N = 50; h = 1.0 / N
y = [i * h for i in range(N + 1)]            # straight-line first guess
def J(y):
    return sum(((y[i + 1] - y[i]) / h) ** 2 * h + (0.5 * (y[i] + y[i + 1])) ** 2 * h for i in range(N))
print(f"straight line: J = {J(y):.6f}")
for sweep in range(3000):                     # Gauss-Seidel on the discrete Euler-Lagrange equation
    for i in range(1, N):
        y[i] = (y[i + 1] + y[i - 1]) / (2 + h * h)
exact = [math.sinh(i * h) / math.sinh(1) for i in range(N + 1)]
print(f"after relaxation: J = {J(y):.6f}, exact minimum = {1 / math.tanh(1):.6f}")
print(f"y(0.5) = {y[N // 2]:.6f}, exact sinh(0.5)/sinh(1) = {exact[N // 2]:.6f}")

# 2. Fermat's principle: light crosses from air (n = 1.00) to water (n = 1.33).
#    Find the crossing point that minimizes travel time and test Snell's law.
n1, n2 = 1.00, 1.33
A, B = (0.0, 1.0), (1.0, -1.0)              # source above the surface y = 0, target below
time = lambda x: n1 * math.hypot(x - A[0], A[1]) + n2 * math.hypot(B[0] - x, B[1])
lo, hi = 0.0, 1.0
for _ in range(200):                         # golden-section search
    a = hi - 0.618034 * (hi - lo); b = lo + 0.618034 * (hi - lo)
    lo, hi = (lo, b) if time(a) < time(b) else (a, hi)
x = 0.5 * (lo + hi)
s1 = (x - A[0]) / math.hypot(x - A[0], A[1]); s2 = (B[0] - x) / math.hypot(B[0] - x, B[1])
print(f"crossing point x = {x:.6f}; n1 sin(theta1) = {n1 * s1:.6f}, n2 sin(theta2) = {n2 * s2:.6f}")

# 3. The Beltrami identity: if F does not depend on x, then F - y' dF/dy' is constant on the extremal.
#    Check it for the minimal surface of revolution F = y sqrt(1 + y'^2), whose extremal is y = c cosh(x/c).
c = 0.8
for xx in (-0.6, -0.2, 0.0, 0.3, 0.7):
    yy = c * math.cosh(xx / c); yp = math.sinh(xx / c)
    F = yy * math.sqrt(1 + yp * yp); Fp = yy * yp / math.sqrt(1 + yp * yp)
    print(f"x = {xx:+.1f}: F - y' dF/dy' = {F - yp * Fp:.6f}")

straight line: J = 1.333300
after relaxation: J = 1.313011, exact minimum = 1.313035
y(0.5) = 0.443411, exact sinh(0.5)/sinh(1) = 0.443409
crossing point x = 0.587867; n1 sin(theta1) = 0.506784, n2 sin(theta2) = 0.506784
x = -0.6: F - y' dF/dy' = 0.800000
x = -0.2: F - y' dF/dy' = 0.800000
x = +0.0: F - y' dF/dy' = 0.800000
x = +0.3: F - y' dF/dy' = 0.800000
x = +0.7: F - y' dF/dy' = 0.800000

The relaxed function reaches \(J = 1.31301\), within \(3\times10^{-5}\) of the exact \(\coth 1\), and its midpoint value 0.443411 matches \(\sinh 0.5/\sinh 1\) to six figures; the tiny gap comes from replacing the integral by a sum on 50 intervals. The straight line, which satisfies the end conditions but not the Euler–Lagrange equation, gives the larger value 1.3333. The quickest crossing point for light gives \(n_1\sin\theta_1\) and \(n_2\sin\theta_2\) equal to all six printed digits. Finally, for the minimal surface of revolution, whose integrand \(y\sqrt{1 + y'^2}\) has no explicit \(x\), the Beltrami expression evaluates to the same constant, 0.800000, at every point of the extremal \(y = 0.8\cosh(x/0.8)\).

Several functions and several variables

If the functional depends on several functions \(y_1(x), \dots, y_n(x)\), each can be varied independently, and the same argument gives one Euler–Lagrange equation per function. This is the multi-coordinate case of mechanics. If the unknown is a function of several variables, such as the shape \(u(x, y)\) of a soap film or the displacement of a drum, the functional is a multiple integral and integration by parts in each variable gives a partial differential equation. For example, minimizing \(\iint (u_x^2 + u_y^2)\,dx\,dy\) with fixed boundary values gives Laplace's equation \(u_{xx} + u_{yy} = 0\). This is how field theories, from elasticity to electromagnetism, are written in Lagrangian form; lesson 12 sketches the idea.

Constraints on the whole function

Sometimes the competing functions must satisfy an integral constraint, such as a fixed total length. The classic case is Dido's problem: enclose the largest area with a rope of given length. The method is a Lagrange multiplier, as in ordinary constrained optimization (MATH-372 used the same idea for finite-dimensional problems). To make \(J[y] = \int F\,dx\) stationary subject to \(K[y] = \int G\,dx = \ell\), make \(\int (F - \lambda G)\,dx\) stationary with \(\lambda\) a constant to be found from the constraint. For the hanging chain in lesson 3, \(F\) is the potential energy and \(G\) the length.