Variational Calculus

The calculus of variations is concerned with finding functions that extremize (minimize or maximize) functionals — mappings from a space of functions to the real numbers. This framework lies at the heart of classical mechanics (Hamilton's principle), optics (Fermat's principle), general relativity (the Einstein-Hilbert action), and quantum field theory (path integrals). Here we develop the core machinery from first principles, derive the Euler-Lagrange equation, and apply it to celebrated problems including the brachistochrone and geodesics.

1. Functionals and Their Extrema

Ordinary calculus deals with functions $f: \mathbb{R} \to \mathbb{R}$ and finding points $x^*$ where$f'(x^*) = 0$. The calculus of variations generalizes this to functionals, which are mappings from a space of functions to the reals.

Definition. A functional $J[y]$ assigns a real number to each function$y(x)$ in some admissible class. The canonical form is:

$$J[y] = \int_a^b F(x, y, y')\, dx$$

where $F$ is a given function of three arguments called the Lagrangian (or integrand), and $y' = dy/dx$. The square brackets in $J[y]$ are conventional notation emphasizing that$J$ depends on the entire function $y(x)$, not just on its value at a point.

Examples of Functionals

Arc length: $J[y] = \int_a^b \sqrt{1 + (y')^2}\, dx$ with $F = \sqrt{1 + (y')^2}$.
Surface of revolution area: $J[y] = 2\pi \int_a^b y\sqrt{1 + (y')^2}\, dx$.
Action in mechanics: $S[q] = \int_{t_1}^{t_2} L(q, \dot{q}, t)\, dt$ where $L = T - V$ is the Lagrangian.
Optical path length: $J[y] = \int_a^b n(x,y)\sqrt{1 + (y')^2}\, dx$ where $n$ is the refractive index.

Boundary Conditions

We typically seek extrema among functions satisfying fixed boundary conditions:

$$y(a) = y_a, \qquad y(b) = y_b$$

The admissible class consists of all sufficiently smooth functions satisfying these endpoint conditions. A function $y^*(x)$ is an extremal if $J[y^*]$ is stationary with respect to all admissible variations. To make this precise, we introduce the concept of the first variation.

2. The Euler-Lagrange Equation: Full Derivation

Let $y^*(x)$ be the extremal function. Consider a one-parameter family of comparison functions:

$$Y(x) = y^*(x) + \epsilon\, \eta(x)$$

where $\epsilon$ is a small parameter and $\eta(x)$ is an arbitrary smooth function satisfying$\eta(a) = \eta(b) = 0$ (so the boundary conditions are preserved). The functional becomes a function of $\epsilon$:

$$\Phi(\epsilon) = J[y^* + \epsilon\eta] = \int_a^b F(x,\, y^* + \epsilon\eta,\, y^{*\prime} + \epsilon\eta')\, dx$$

For $y^*$ to be an extremal, we need $\Phi'(0) = 0$ for all admissible $\eta$. Differentiating under the integral sign:

$$\Phi'(\epsilon) = \int_a^b \left[\frac{\partial F}{\partial y}\,\eta + \frac{\partial F}{\partial y'}\,\eta'\right] dx$$

Setting $\epsilon = 0$ gives the first variation:

$$\delta J = \Phi'(0) = \int_a^b \left[\frac{\partial F}{\partial y}\,\eta + \frac{\partial F}{\partial y'}\,\eta'\right] dx = 0$$

Integration by Parts

The key step is to integrate the second term by parts. Let $u = \partial F/\partial y'$ and $dv = \eta'\, dx$:

$$\int_a^b \frac{\partial F}{\partial y'}\,\eta'\, dx = \left[\frac{\partial F}{\partial y'}\,\eta\right]_a^b - \int_a^b \frac{d}{dx}\frac{\partial F}{\partial y'}\,\eta\, dx$$

The boundary term vanishes because $\eta(a) = \eta(b) = 0$. Therefore:

$$\delta J = \int_a^b \left[\frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'}\right]\eta(x)\, dx = 0$$

The Fundamental Lemma

We invoke the fundamental lemma of the calculus of variations: if$\int_a^b g(x)\,\eta(x)\,dx = 0$ for all smooth $\eta$ vanishing at the endpoints, and $g$ is continuous, then $g(x) \equiv 0$ on $[a,b]$.

Applying the lemma yields the Euler-Lagrange equation:

$$\boxed{\frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'} = 0}$$

This is a second-order ODE for $y(x)$ (since the total derivative $d/dx$ of $\partial F/\partial y'$generally involves $y''$). It is the necessary condition for $y$ to be an extremal.

Beltrami Identity (When F Does Not Depend on x Explicitly)

If $F = F(y, y')$ (no explicit $x$-dependence), there is a first integral. Compute:

$$\frac{d}{dx}\left[F - y'\frac{\partial F}{\partial y'}\right] = \frac{\partial F}{\partial y}y' + \frac{\partial F}{\partial y'}y'' - y''\frac{\partial F}{\partial y'} - y'\frac{d}{dx}\frac{\partial F}{\partial y'} = y'\left[\frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'}\right] = 0$$

where we used the Euler-Lagrange equation in the last step. Thus:

$$\boxed{F - y'\frac{\partial F}{\partial y'} = C \quad \text{(Beltrami identity)}}$$

This reduces the problem from a second-order to a first-order ODE, which is a significant simplification. In mechanics, when the Lagrangian has no explicit time dependence, the Beltrami identity gives conservation of energy.

3. The Brachistochrone Problem

The brachistochrone ("shortest time" in Greek) was posed by Johann Bernoulli in 1696: find the curve connecting two points along which a bead slides (frictionlessly, under gravity) in the least time.

Setting Up the Functional

Let the bead start from rest at the origin and slide to point $(x_1, y_1)$ with $y$ measured downward. Conservation of energy gives the speed:

$$v = \sqrt{2gy}$$

The time to traverse an infinitesimal arc length $ds = \sqrt{dx^2 + dy^2}$ is $dt = ds/v$. Using$y$ as the independent variable and writing $x = x(y)$:

$$T[x] = \int_0^{y_1} \frac{\sqrt{1 + (x')^2}}{\sqrt{2gy}}\, dy$$

Here $x' = dx/dy$ and the integrand is $F(y, x, x') = \sqrt{(1 + x'^2)/(2gy)}$. Since $F$ does not depend explicitly on $x$, the Euler-Lagrange equation gives a first integral (the Beltrami identity adapted for the absent variable):

$$\frac{\partial F}{\partial x'} = \frac{x'}{\sqrt{2gy}\sqrt{1+x'^2}} = \text{const}$$

Solution: The Cycloid

After algebraic manipulation, setting the constant as $1/\sqrt{2g \cdot 2a}$, we find:

$$x'^2 = \frac{2a - y}{y} \quad \Longrightarrow \quad dx = \sqrt{\frac{2a - y}{y}}\, dy$$

The substitution $y = a(1 - \cos\theta)$ yields $dy = a\sin\theta\, d\theta$ and:

$$dx = \sqrt{\frac{a(1+\cos\theta)}{a(1-\cos\theta)}}\, a\sin\theta\, d\theta = a(1 - \cos\theta)\, d\theta$$

where we used the half-angle identities. Integrating:

$$\boxed{x = a(\theta - \sin\theta), \qquad y = a(1 - \cos\theta)}$$

This is the parametric equation of a cycloid — the curve traced by a point on the rim of a circle of radius $a$ rolling along a horizontal line. The parameter $a$ is determined by the endpoint condition $(x_1, y_1)$.

4. Geodesics

A geodesic is a curve of extremal length on a surface or in a Riemannian manifold. In general relativity, freely falling particles follow geodesics of spacetime.

Geodesics in a General Metric

On a surface with metric $ds^2 = g_{ij}\, dx^i\, dx^j$, the arc length functional is:

$$L[\gamma] = \int_{\lambda_1}^{\lambda_2} \sqrt{g_{ij}\frac{dx^i}{d\lambda}\frac{dx^j}{d\lambda}}\, d\lambda$$

Applying the Euler-Lagrange equation to this functional yields the geodesic equation:

$$\boxed{\frac{d^2 x^\mu}{d\lambda^2} + \Gamma^\mu_{\alpha\beta}\frac{dx^\alpha}{d\lambda}\frac{dx^\beta}{d\lambda} = 0}$$

where $\Gamma^\mu_{\alpha\beta}$ are the Christoffel symbols derived from the metric.

Geodesics on a Sphere

On a sphere of radius $R$, the metric is $ds^2 = R^2(d\theta^2 + \sin^2\theta\, d\phi^2)$. Using $\theta$ as the independent variable, the arc length becomes:

$$L[\phi] = R\int_{\theta_1}^{\theta_2} \sqrt{1 + \sin^2\theta\, (\phi')^2}\, d\theta$$

Since $F$ does not depend on $\phi$ explicitly, the Euler-Lagrange equation gives:

$$\frac{\partial F}{\partial \phi'} = \frac{\sin^2\theta\, \phi'}{\sqrt{1 + \sin^2\theta\, (\phi')^2}} = c$$

This is Clairaut's relation. Solving for $\phi'$ and integrating gives great circles as the geodesics on a sphere, confirming that the shortest path between two points on a sphere is an arc of a great circle.

5. Isoperimetric Problems and Lagrange Multipliers

The classical isoperimetric problem asks: among all closed curves of a given perimeter, which encloses the maximum area? The answer is a circle. More generally, we want to extremize a functional subject to constraints.

Constrained Variational Problems

Suppose we wish to extremize $J[y] = \int_a^b F(x, y, y')\, dx$ subject to the constraint:

$$K[y] = \int_a^b G(x, y, y')\, dx = \ell \quad \text{(given constant)}$$

By analogy with finite-dimensional optimization, we introduce a Lagrange multiplier $\lambda$and form the augmented functional:

$$\tilde{J}[y] = J[y] + \lambda\, K[y] = \int_a^b \left[F + \lambda\, G\right] dx$$

The extremal satisfies the Euler-Lagrange equation for the combined integrand $\tilde{F} = F + \lambda G$:

$$\boxed{\frac{\partial}{\partial y}(F + \lambda G) - \frac{d}{dx}\frac{\partial}{\partial y'}(F + \lambda G) = 0}$$

Example: Maximum Area for Fixed Perimeter

For a curve $y(x)$ from $(0,0)$ to $(L,0)$, the enclosed area is$A[y] = \int_0^L y\, dx$ and the arc length constraint is$\int_0^L \sqrt{1 + y'^2}\, dx = \ell$. The augmented integrand is:

$$\tilde{F} = y + \lambda\sqrt{1 + y'^2}$$

The Euler-Lagrange equation yields:

$$1 - \lambda\frac{d}{dx}\frac{y'}{\sqrt{1 + y'^2}} = 0 \quad \Longrightarrow \quad \frac{y'}{\sqrt{1+y'^2}} = \frac{x - c_1}{\lambda}$$

Solving gives $(x - c_1)^2 + (y - c_2)^2 = \lambda^2$, which is a circle of radius $|\lambda|$. This confirms that among all curves of fixed length, the circle encloses the greatest area.

Multiple Constraints

With $m$ constraints $K_j[y] = \ell_j$ for $j = 1, \ldots, m$, we introduce multipliers$\lambda_1, \ldots, \lambda_m$ and extremize:

$$\tilde{J}[y] = \int_a^b \left[F + \sum_{j=1}^{m}\lambda_j\, G_j\right] dx$$

This technique extends naturally to higher-dimensional problems and is the foundation of the method used in quantum mechanics for constrained optimization (e.g., the Rayleigh-Ritz method for finding bound-state energies).

6. Hamilton's Principle

Hamilton's principle (the principle of stationary action) states that the physical trajectory of a mechanical system between times $t_1$ and $t_2$ is the one that makes the action stationary:

$$\boxed{\delta S = \delta \int_{t_1}^{t_2} L(q_i, \dot{q}_i, t)\, dt = 0}$$

where $L = T - V$ is the Lagrangian, $q_i$ are generalized coordinates, and $\dot{q}_i = dq_i/dt$.

Derivation of Lagrange's Equations from Hamilton's Principle

Applying the Euler-Lagrange equation to the action with $F = L$, $y \to q_i$, $y' \to \dot{q}_i$, $x \to t$:

$$\frac{\partial L}{\partial q_i} - \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} = 0, \qquad i = 1, \ldots, n$$

These are the Lagrange equations of motion. For a particle with$L = \frac{1}{2}m\dot{x}^2 - V(x)$:

$$\frac{\partial L}{\partial x} = -\frac{dV}{dx}, \qquad \frac{\partial L}{\partial \dot{x}} = m\dot{x}$$

$$-\frac{dV}{dx} - \frac{d}{dt}(m\dot{x}) = 0 \quad \Longrightarrow \quad m\ddot{x} = -\frac{dV}{dx} = F$$

recovering Newton's second law.

The Hamiltonian and Legendre Transform

The Beltrami identity applied to the action (when $L$ has no explicit time dependence) gives:

$$H = \sum_i \dot{q}_i \frac{\partial L}{\partial \dot{q}_i} - L = \text{const}$$

This is the Hamiltonian, and its conservation is the energy conservation theorem. Defining the conjugate momenta $p_i = \partial L / \partial \dot{q}_i$, we write$H = H(q_i, p_i, t)$ via the Legendre transform. Hamilton's equations follow:

$$\boxed{\dot{q}_i = \frac{\partial H}{\partial p_i}, \qquad \dot{p}_i = -\frac{\partial H}{\partial q_i}}$$

7. Noether's Theorem from the Variational Perspective

Emmy Noether's theorem (1918) establishes a deep connection between symmetries and conservation laws. We derive it within the variational framework.

Statement

Theorem (Noether): If the action $S = \int L\, dt$ is invariant under a continuous one-parameter family of transformations $q_i \to q_i + \epsilon\, \xi_i(q, t)$,$t \to t + \epsilon\, \tau(q, t)$, then the quantity

$$\boxed{I = \sum_i \frac{\partial L}{\partial \dot{q}_i}\xi_i - \left(\sum_i \frac{\partial L}{\partial \dot{q}_i}\dot{q}_i - L\right)\tau}$$

is conserved along physical trajectories ($dI/dt = 0$).

Derivation

Consider an infinitesimal transformation with parameter $\epsilon$:

$$\bar{q}_i = q_i + \epsilon\, \xi_i, \qquad \bar{t} = t + \epsilon\, \tau$$

Invariance of the action means $\delta S = 0$ to first order in $\epsilon$. The change in the action has two contributions: the change in the integrand and the change in the limits. After careful computation (using the chain rule and the Leibniz rule for the endpoint variation):

$$\delta S = \epsilon\int_{t_1}^{t_2}\left[\sum_i\left(\frac{\partial L}{\partial q_i} - \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i}\right)(\xi_i - \dot{q}_i\tau) + \frac{d}{dt}\left(\sum_i \frac{\partial L}{\partial \dot{q}_i}\xi_i - H\tau\right)\right] dt = 0$$

The first group of terms vanishes by the Euler-Lagrange equations (since we evaluate on physical trajectories). Since $\delta S = 0$ and $t_1, t_2$ are arbitrary, the integrand of the remaining total derivative must vanish:

$$\frac{d}{dt}\left[\sum_i p_i\, \xi_i - H\, \tau\right] = 0$$

Applications

Time translation ($\tau = 1, \xi_i = 0$): Conserved quantity is $-H = -(T+V)$, giving energy conservation.
Spatial translation ($\tau = 0, \xi_j = \delta_{ij}$): Conserved quantity is $p_i$, giving momentum conservation.
Rotation ($\tau = 0$, $\xi_i$ generates rotation): Conserved quantity is a component of angular momentum, giving angular momentum conservation.

Noether's theorem has far-reaching consequences in quantum field theory, where gauge symmetries lead to conserved currents (electric charge, color charge, etc.) and the Ward-Takahashi identities.

8. Extensions of the Euler-Lagrange Framework

Several Dependent Variables

For a functional depending on $n$ functions $y_1(x), \ldots, y_n(x)$:

$$J[y_1, \ldots, y_n] = \int_a^b F(x, y_1, \ldots, y_n, y_1', \ldots, y_n')\, dx$$

we obtain $n$ coupled Euler-Lagrange equations:

$$\frac{\partial F}{\partial y_i} - \frac{d}{dx}\frac{\partial F}{\partial y_i'} = 0, \qquad i = 1, \ldots, n$$

Higher Derivatives

If the integrand depends on higher derivatives $y''$, the functional$J = \int_a^b F(x, y, y', y'')\, dx$ leads to the Euler-Poisson equation:

$$\frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'} + \frac{d^2}{dx^2}\frac{\partial F}{\partial y''} = 0$$

This arises in elasticity theory (Euler-Bernoulli beam equation) and in certain higher-derivative field theories.

Functionals of Several Independent Variables

For field theories, we need functionals of the form$J[\phi] = \int \mathcal{L}(\phi, \partial_\mu\phi, x^\mu)\, d^4x$ where $\phi$ is a field. The corresponding Euler-Lagrange equation is:

$$\boxed{\frac{\partial \mathcal{L}}{\partial \phi} - \partial_\mu \frac{\partial \mathcal{L}}{\partial(\partial_\mu \phi)} = 0}$$

This is the equation of motion for a classical field. For the Klein-Gordon Lagrangian$\mathcal{L} = \frac{1}{2}(\partial_\mu\phi)(\partial^\mu\phi) - \frac{1}{2}m^2\phi^2$, the Euler-Lagrange equation gives $(\Box + m^2)\phi = 0$, the Klein-Gordon equation.

9. Python Simulation: Brachistochrone and Geodesics

Below we verify the variational calculus results computationally. We compare the brachistochrone (cycloid) against straight-line and circular-arc paths, compute travel times numerically, and visualize geodesics on a sphere.

Python

script.py214 lines

#!/usr/bin/env python3
"""
variational_calculus_sim.py
Brachistochrone verification and geodesics on a sphere.  numpy only.
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ─── Part 1: Brachistochrone (Cycloid) ──────────────────────────────
g = 9.81

def cycloid(a, theta_arr):
    """Parametric cycloid: x = a(theta - sin(theta)), y = a(1 - cos(theta))"""
    x = a * (theta_arr - np.sin(theta_arr))
    y = a * (1.0 - np.cos(theta_arr))
    return x, y

def travel_time(x, y, g=9.81):
    """Compute travel time along a curve y(x) from rest under gravity.
    y measured downward (positive). Uses trapezoidal integration."""
    dx = np.diff(x)
    dy = np.diff(y)
    ds = np.sqrt(dx**2 + dy**2)
    y_mid = 0.5 * (y[:-1] + y[1:])
    y_mid = np.maximum(y_mid, 1e-12)  # avoid division by zero
    v_mid = np.sqrt(2.0 * g * y_mid)
    dt = ds / v_mid
    return np.sum(dt)

# Target point
x1, y1 = 2.0, 1.0

# Find cycloid parameter a by bisection
def endpoint_error(a, x1, y1, N=2000):
    theta_max_guess = np.pi
    for _ in range(100):
        theta_test = np.linspace(0, theta_max_guess, N)
        xc, yc = cycloid(a, theta_test)
        idx = np.argmin(np.abs(yc - y1))
        if idx == 0:
            theta_max_guess *= 1.5
            continue
        return xc[idx] - x1
    return 1e10

a_low, a_high = 0.1, 10.0
for _ in range(80):
    a_mid = 0.5 * (a_low + a_high)
    if endpoint_error(a_mid, x1, y1) > 0:
        a_high = a_mid
    else:
        a_low = a_mid
a_opt = 0.5 * (a_low + a_high)

# Generate cycloid
N = 5000
theta_arr = np.linspace(0, np.pi, N)
xc, yc = cycloid(a_opt, theta_arr)
mask = yc <= y1 * 1.01
xc, yc = xc[mask], yc[mask]
# Interpolate to exactly hit endpoint
xc_interp = np.linspace(0, x1, N)
yc_interp = np.interp(xc_interp, xc, yc)

# Straight line
xs = np.linspace(0, x1, N)
ys = np.linspace(0, y1, N)

# Circular arc (passing through origin and endpoint)
t_param = np.linspace(0, 1, N)
# Quadratic: y = 4*y1*t*(1-... no, let us use a parabola
xp = x1 * t_param
yp = 4.0 * y1 * t_param * (1.0 - 0.5 * t_param)  # parabolic arc

T_cycloid = travel_time(xc_interp, yc_interp)
T_straight = travel_time(xs, ys)
T_parabola = travel_time(xp, yp)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot brachistochrone comparison
ax = axes[0]
ax.plot(xc_interp, yc_interp, 'b-', lw=2.5, label=f'Cycloid  T={T_cycloid:.4f} s')
ax.plot(xs, ys, 'r--', lw=2, label=f'Straight T={T_straight:.4f} s')
ax.plot(xp, yp, 'g-.', lw=2, label=f'Parabola T={T_parabola:.4f} s')
ax.set_xlabel('x (m)', fontsize=12)
ax.set_ylabel('y (m, downward)', fontsize=12)
ax.set_title('Brachistochrone Problem', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.invert_yaxis()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')

# ─── Part 2: Geodesics on a Sphere ─────────────────────────────────
ax2 = axes[1]

# Draw sphere wireframe
u_sphere = np.linspace(0, 2*np.pi, 60)
v_sphere = np.linspace(0, np.pi, 30)

# We will draw the equatorial and meridional lines
for v in np.linspace(0, np.pi, 7):
    x_circ = np.sin(v) * np.cos(u_sphere)
    y_circ = np.sin(v) * np.sin(u_sphere)
    ax2.plot(x_circ, y_circ, 'gray', alpha=0.15, lw=0.5)

for u in np.linspace(0, np.pi, 7):
    x_merid = np.sin(v_sphere) * np.cos(u)
    y_merid = np.cos(v_sphere)
    ax2.plot(x_merid, y_merid, 'gray', alpha=0.15, lw=0.5)

# Great circle geodesic between two points
theta1, phi1 = np.pi/6, 0.0
theta2, phi2 = 2*np.pi/3, np.pi/3

# Convert to Cartesian
def sph_to_cart(theta, phi):
    return np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)

p1 = np.array(sph_to_cart(theta1, phi1))
p2 = np.array(sph_to_cart(theta2, phi2))

# Great circle via slerp
dot = np.clip(np.dot(p1, p2), -1, 1)
omega = np.arccos(dot)
t_gc = np.linspace(0, 1, 200)
pts = np.array([(np.sin((1-t)*omega)/np.sin(omega))*p1 +
                (np.sin(t*omega)/np.sin(omega))*p2 for t in t_gc])

# Project to 2D (stereographic-like: just use x, z for visualization)
ax2.plot(pts[:, 0], pts[:, 2], 'b-', lw=3, label='Geodesic (great circle)')

# Non-geodesic path (linear interpolation in spherical coords)
theta_lin = theta1 + t_gc * (theta2 - theta1)
phi_lin = phi1 + t_gc * (phi2 - phi1)
x_lin = np.sin(theta_lin) * np.cos(phi_lin)
z_lin = np.cos(theta_lin)
ax2.plot(x_lin, z_lin, 'r--', lw=2, label='Non-geodesic path')

# Mark endpoints
ax2.plot([p1[0]], [p1[2]], 'ko', ms=8, zorder=5)
ax2.plot([p2[0]], [p2[2]], 'ko', ms=8, zorder=5)
ax2.annotate('P1', (p1[0]+0.05, p1[2]+0.05), fontsize=12, color='white')
ax2.annotate('P2', (p2[0]+0.05, p2[2]+0.05), fontsize=12, color='white')

# Unit circle outline
circ_t = np.linspace(0, 2*np.pi, 200)
ax2.plot(np.cos(circ_t), np.sin(circ_t), 'gray', lw=1, alpha=0.5)

ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('z', fontsize=12)
ax2.set_title('Geodesics on Unit Sphere (x-z projection)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10, loc='lower left')
ax2.set_aspect('equal')
ax2.grid(True, alpha=0.3)
ax2.set_facecolor('#0a0a2e')

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight',
            facecolor='#0a0a2e', edgecolor='none')
plt.close()

# ─── Part 3: Numerical Verification ────────────────────────────────
print("=" * 65)
print("VARIATIONAL CALCULUS: NUMERICAL RESULTS")
print("=" * 65)

print(f"\nBrachistochrone Problem: endpoint ({x1}, {y1})")
print(f"  Cycloid parameter a = {a_opt:.6f}")
print(f"  Travel time (cycloid):   {T_cycloid:.6f} s")
print(f"  Travel time (straight):  {T_straight:.6f} s")
print(f"  Travel time (parabola):  {T_parabola:.6f} s")
print(f"  Cycloid is fastest: {T_cycloid < T_straight and T_cycloid < T_parabola}")

# Arc lengths
def arc_length_sphere(pts):
    diffs = np.diff(pts, axis=0)
    return np.sum(np.sqrt(np.sum(diffs**2, axis=1)))

L_geodesic = arc_length_sphere(pts)
pts_nongeo = np.column_stack([x_lin, np.sin(theta_lin)*np.sin(phi_lin), z_lin])
L_nongeo = arc_length_sphere(pts_nongeo)

print(f"\nGeodesics on Unit Sphere:")
print(f"  Arc length (great circle): {L_geodesic:.6f}")
print(f"  Arc length (non-geodesic): {L_nongeo:.6f}")
print(f"  Great circle is shorter:   {L_geodesic < L_nongeo}")
print(f"  Theoretical great-circle:  {omega:.6f}")
print(f"  Numerical error:           {abs(L_geodesic - omega):.2e}")

# Verify Euler-Lagrange for simple case: shortest path in plane
print(f"\nEuler-Lagrange verification (flat plane):")
print(f"  F = sqrt(1 + y'^2), extremal should be straight line")
print(f"  For y = mx + b: y'' = 0")
print(f"  E-L: d/dx[y'/sqrt(1+y'^2)] = 0 => y'/sqrt(1+y'^2) = const")
print(f"  This is satisfied iff y' = const, i.e., straight line. QED")

# Numerical E-L residual for cycloid
print(f"\nEuler-Lagrange residual for brachistochrone:")
dx_num = xc_interp[1] - xc_interp[0]
yp_num = np.gradient(yc_interp, dx_num)
ypp_num = np.gradient(yp_num, dx_num)
# F = sqrt(1+y'^2) / sqrt(2gy), E-L residual
F_y = -0.5 * np.sqrt(1 + yp_num**2) / (np.sqrt(2*g) * yc_interp**1.5 + 1e-15)
dF_yp = yp_num / (np.sqrt(2*g * yc_interp + 1e-15) * np.sqrt(1 + yp_num**2))
d_dF_yp = np.gradient(dF_yp, dx_num)
residual = F_y - d_dF_yp
# Skip endpoints (numerical issues)
interior = residual[100:-100]
print(f"  Max |residual| (interior): {np.max(np.abs(interior)):.6e}")
print(f"  Mean |residual|:           {np.mean(np.abs(interior)):.6e}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

10. Summary and Connections

The calculus of variations provides a unifying framework that permeates theoretical physics:

Concept	Key Result	Physics Application
Euler-Lagrange equation	$\partial_y F - \frac{d}{dx}\partial_{y'} F = 0$	Equations of motion
Beltrami identity	$F - y' F_{y'} = C$	Energy conservation
Brachistochrone	Cycloid	Optimal transport
Geodesic equation	$\ddot{x}^\mu + \Gamma^\mu_{\alpha\beta}\dot{x}^\alpha\dot{x}^\beta = 0$	General relativity
Isoperimetric / multipliers	Augmented Euler-Lagrange	Constrained optimization
Hamilton's principle	$\delta S = 0$	All of classical mechanics
Noether's theorem	Symmetry $\Leftrightarrow$ conservation law	QFT, gauge theories

The variational perspective reveals that the fundamental laws of physics are not arbitrary differential equations but rather consequences of extremizing action functionals. This insight extends to quantum mechanics through Feynman's path integral formulation, where a particle explores all possible paths weighted by$e^{iS/\hbar}$, and the classical trajectory emerges in the limit $\hbar \to 0$ via stationary phase.

← Part IV Overview Integral Equations →