Partial differential equations: modern theory

To find the equilibrium temperature $u$ in a region $\Omega$ with prescribed boundary values, nineteenth-century physicists minimised the Dirichlet energy $E(u) = \frac12 \int_\Omega |\nabla u|^2$: a minimiser satisfies Laplace's equation $\Delta u = 0$. Riemann (1851) built his theory of Riemann surfaces on this "Dirichlet principle". In 1870 Weierstrass pointed out that an energy bounded below need not attain its infimum: he exhibited a similar functional whose infimum over smooth functions is never achieved. Existence was not a formality, and the principle fell from favour until Hilbert (1900) rescued it.

The rescue, completed by Sobolev, Leray, Schwartz and others in the 1930s and 1940s, is the modern theory. Its strategy has two steps. First, enlarge the space of candidates to a complete space — a Sobolev space of functions with square-integrable weak derivatives — so that minimising sequences have limits and existence follows from functional analysis. Second, prove regularity: show that the weak solution is in fact smooth. Many problems stop at the first step, because their natural solutions are not smooth: shock waves in gas dynamics, the bending of a plate with a corner, the turbulent velocity fields of the Navier–Stokes equations.

If $u \in C^2(\overline\Omega)$ solves $-\Delta u = f$ with $u = 0$ on $\partial\Omega$, multiply by any smooth $v$ vanishing near $\partial\Omega$ and integrate by parts:

This weak formulation makes sense when $u$ has only one derivative, and only in $L^2$. Conversely, a $C^2$ function satisfying it for all $v$ solves the equation classically. So we may define a solution by the integral identity, look for it in the space where both sides are continuous in $u$ and $v$, and recognise the left side as an inner product. The equation becomes: find $u$ with $\langle u, v \rangle_{\text{energy}} = \ell(v)$ for all $v$ — and the Riesz representation theorem solves that instantly. Every quantitative estimate (stability, convergence of numerical schemes) comes from this Hilbert-space picture, with coercivity (here the Poincaré inequality) as the key input.

The same move — test against smooth functions, put derivatives on the test functions — defines derivatives of objects that are not functions at all: the Dirac delta, the derivative of a jump. That is the theory of distributions.

Let $\Omega \subseteq \R^n$ be open. Multi-indices $\alpha = (\alpha_1, \ldots, \alpha_n)$ index derivatives $\partial^\alpha = \partial_1^{\alpha_1} \cdots \partial_n^{\alpha_n}$, of order $|\alpha| = \sum \alpha_i$.

$H^1(\Omega)$ is a Hilbert space: if $(u_k)$ is Cauchy in $H^1$, then $u_k \to u$ and $\partial_i u_k \to g_i$ in $L^2$ by Riesz–Fischer, and passing to the limit in $\int u_k \partial_i \varphi = -\int \partial_i u_k \varphi$ shows $g_i = \partial_i u$. Its dual space $H^{-1}(\Omega) = (H^1_0(\Omega))^*$ contains every $L^2$ function and also distributions such as $\partial_i f$ with $f \in L^2$.

If $B$ is also symmetric, it is an inner product equivalent to the original one and Lax–Milgram is just the Riesz theorem; moreover $u$ is the unique minimiser of $J(v) = \frac12 B(v, v) - \ell(v)$, since $J(u + w) = J(u) + \frac12 B(w, w)$. That is Dirichlet's principle, now with a proof. The force of Lax–Milgram is that it also covers non-symmetric forms, such as those of convection–diffusion equations $-\Delta u + b \cdot \nabla u = f$.

The inequality fails in $H^1(\Omega)$ (constants have $\nabla u = 0$) and on unbounded domains such as $\R^n$: the boundary condition and the boundedness both matter.

The same argument handles $-\operatorname{div}(a(x) \nabla u) + b(x) \cdot \nabla u + c(x) u = f$ whenever the coefficients are bounded and the form is coercive (e.g. $a \ge a_0 > 0$, $b = 0$, $c \ge 0$), with no smoothness of $a$ at all — the setting of composite materials, where classical solutions do not exist. Combined with the Rellich–Kondrachov theorem (the embedding $H^1_0(\Omega) \hookrightarrow L^2(\Omega)$ is compact for bounded $\Omega$), it makes the solution operator $f \mapsto u$ a compact self-adjoint operator on $L^2(\Omega)$, and the spectral theorem of functional analysis yields an orthonormal basis of Dirichlet eigenfunctions $-\Delta \varphi_k = \lambda_k \varphi_k$ with $0 < \lambda_1 \le \lambda_2 \le \cdots \to \infty$.

So the finite element error is, up to the factor $M/c$, the best approximation error of the finite element space — a question of approximation theory (for piecewise linear elements of size $h$, of order $h \|u\|_{H^2}$).

Consequences: the classical Dirichlet problem has at most one solution (the difference of two solutions is harmonic with zero boundary values), and solutions depend continuously on the boundary data in the sup norm, $\max_{\overline\Omega} |u_1 - u_2| = \max_{\partial\Omega} |u_1 - u_2|$. The same argument proves the weak maximum principle for the heat equation on space–time cylinders. In one dimension "subharmonic" means $u'' \ge 0$, i.e. convex, and the principle says a convex function on an interval lies below the larger of its endpoint values — the version formalised below.

Energy methods are robust: they extend to variable coefficients and nonlinear problems, and localised versions prove finite speed of propagation for the wave equation (a disturbance supported in a ball of radius $r$ at time $0$ stays in the ball of radius $r + t$), in contrast with the heat equation, whose kernel is positive everywhere for $t > 0$.

Elliptic regularity. Weak solutions of elliptic equations are as smooth as the data allow. If $u \in H^1(\Omega)$ is a weak solution of $-\Delta u = f$ and $f \in H^k_{\mathrm{loc}}(\Omega)$, then $u \in H^{k+2}_{\mathrm{loc}}(\Omega)$; by Sobolev embedding ($H^m_{\mathrm{loc}} \subseteq C^j$ when $m > j + n/2$), $f \in C^\infty$ implies $u \in C^\infty(\Omega)$. Weyl's lemma (1940) is the case $f = 0$: a distribution with $\Delta T = 0$ is a smooth harmonic function. The proof (Evans, §6.3) differentiates the equation weakly using Nirenberg's difference quotients $\frac{u(x + he_i) - u(x)}{h}$, which are bounded in $H^1$ by the energy estimate and therefore converge weakly to $\partial_i u \in H^1$; iterating gains one derivative at a time. Up to the boundary the same holds if $\partial\Omega$ is smooth. For nonlinear elliptic equations regularity is much harder: Hilbert's nineteenth problem, on the regularity of minimisers of variational integrals, was solved independently by De Giorgi (1957) and Nash (1958).

Well-posedness can fail. Hadamard's example: $u_n(x, y) = \frac{\sin(nx) \sinh(ny)}{n^2}$ solves $\Delta u_n = 0$ with Cauchy data $u_n(x, 0) = 0$ and $\partial_y u_n(x, 0) = \frac{\sin(nx)}{n}$. The data tend to $0$ uniformly, yet $u_n\big(\frac{\pi}{2n}, 1\big) = \frac{\sinh n}{n^2} \to \infty$. The zero data have the solution $u = 0$, so the Cauchy problem for Laplace's equation does not depend continuously on its data: it is ill-posed, although (by Cauchy–Kovalevskaya) it is uniquely solvable for analytic data. The right problems for elliptic equations are boundary value problems; Cauchy problems belong to hyperbolic equations such as the wave equation. Ill-posed problems are not useless — inverse problems in imaging are ill-posed — but they require regularisation.

Weak solutions can be the only solutions. For Burgers' equation $u_t + (u^2/2)_x = 0$, smooth initial data develop shocks in finite time, after which only weak (integral) solutions exist; an extra entropy condition restores uniqueness (Oleinik, Kružkov). And for the three-dimensional Navier–Stokes equations, Leray (1934) proved that global weak solutions exist; whether smooth data always lead to global smooth solutions is one of the Clay Millennium Prize problems, and is open.

Mathlib has the Lax–Milgram theorem in the form "the operator $A$ of the proof is a continuous linear equivalence". Combined with the Riesz representation it gives the statement of the lesson:

/-- **Lax–Milgram theorem.** Let `B` be a bounded bilinear form on a real Hilbert space `V`
which is coercive (`C ‖u‖² ≤ B(u, u)` for some `C > 0`). Then every bounded linear functional
`ℓ` is represented by `B`: there is `u` with `B(u, v) = ℓ(v)` for all `v`. This combines
Mathlib's `IsCoercive.continuousLinearEquivOfBilin` (the operator `u ↦ B(u, ·)` is an
isomorphism) with the Riesz representation `InnerProductSpace.toDual`. -/
theorem lax_milgram {V : Type*} [NormedAddCommGroup V] [InnerProductSpace ℝ V] [CompleteSpace V]
    {B : V →L[ℝ] V →L[ℝ] ℝ} (hB : IsCoercive B) (ℓ : StrongDual ℝ V) :
    ∃ u : V, ∀ v, B u v = ℓ v := by
  refine ⟨hB.continuousLinearEquivOfBilin.symm ((InnerProductSpace.toDual ℝ V).symm ℓ),
    fun v => ?_⟩
  rw [← IsCoercive.continuousLinearEquivOfBilin_apply hB,
    ContinuousLinearEquiv.apply_symm_apply, InnerProductSpace.toDual_symm_apply]

Uniqueness and the stability estimate $c\|u\| \le \|\ell\|$, proved from coercivity exactly as above:

/-- **Uniqueness of weak solutions.** If `B` is coercive and `B(u₁, v) = B(u₂, v)` for all `v`,
then `u₁ = u₂`: with `w = u₁ − u₂`, coercivity gives `C ‖w‖² ≤ B(w, w) = 0`. -/
theorem lax_milgram_unique {V : Type*} [NormedAddCommGroup V] [InnerProductSpace ℝ V]
    {B : V →L[ℝ] V →L[ℝ] ℝ} (hB : IsCoercive B) {u₁ u₂ : V} (h : ∀ v, B u₁ v = B u₂ v) :
    u₁ = u₂ := by
  obtain ⟨C, hC, hcoer⟩ := hB
  have h0 : B (u₁ - u₂) (u₁ - u₂) = 0 := by
    rw [map_sub B u₁ u₂, sub_apply, h, sub_self]
  have hw := hcoer (u₁ - u₂)
  rw [h0] at hw
  have hn : ‖u₁ - u₂‖ = 0 := by
    by_contra hne
    have hpos : 0 < ‖u₁ - u₂‖ := lt_of_le_of_ne (norm_nonneg _) (Ne.symm hne)
    have := mul_pos (mul_pos hC hpos) hpos
    linarith
  exact sub_eq_zero.mp (norm_eq_zero.mp hn)

/-- **Stability (continuous dependence on the data).** If `C ‖u‖² ≤ B(u, u)` and
`B(u, v) = ℓ(v)` for all `v`, then `C ‖u‖ ≤ ‖ℓ‖`: indeed
`C ‖u‖² ≤ B(u, u) = ℓ(u) ≤ ‖ℓ‖ ‖u‖`; divide by `‖u‖` (the case `u = 0` is trivial). For
`−Δu = f` this is the energy estimate `‖u‖_{H¹₀} ≤ C⁻¹ ‖f‖_{H⁻¹}`. -/
theorem lax_milgram_bound {V : Type*} [NormedAddCommGroup V] [InnerProductSpace ℝ V]
    {B : V →L[ℝ] V →L[ℝ] ℝ} {C : ℝ} (hcoer : ∀ u, C * ‖u‖ * ‖u‖ ≤ B u u)
    (ℓ : StrongDual ℝ V) {u : V} (hu : ∀ v, B u v = ℓ v) : C * ‖u‖ ≤ ‖ℓ‖ := by
  rcases eq_or_lt_of_le (norm_nonneg u) with h | h
  · rw [← h, mul_zero]
    exact norm_nonneg ℓ
  · have h1 : C * ‖u‖ * ‖u‖ ≤ ‖ℓ‖ * ‖u‖ :=
      calc C * ‖u‖ * ‖u‖ ≤ B u u := hcoer u
        _ = ℓ u := hu u
        _ ≤ |ℓ u| := le_abs_self _
        _ = ‖ℓ u‖ := (Real.norm_eq_abs _).symm
        _ ≤ ‖ℓ‖ * ‖u‖ := ℓ.le_opNorm u
    exact le_of_mul_le_mul_right h1 h

The bridge between classical and weak derivatives — integration by parts against a test function vanishing at the endpoints:

/-- **Classical derivatives are weak derivatives.** If `u` is differentiable on `[a, b]` with
integrable derivative `u'`, then for every test function `φ` vanishing at `a` and `b`,
`∫ u φ' = − ∫ u' φ`. This identity is the *definition* of "`u'` is the weak derivative of `u`",
so the weak notion extends the classical one. Proof: integration by parts (Mathlib's
`intervalIntegral.integral_mul_deriv_eq_deriv_mul`); the boundary terms vanish. -/
theorem classical_is_weak {u u' φ φ' : ℝ → ℝ} {a b : ℝ}
    (hu : ∀ x ∈ uIcc a b, HasDerivAt u (u' x) x) (hφ : ∀ x ∈ uIcc a b, HasDerivAt φ (φ' x) x)
    (hu' : IntervalIntegrable u' MeasureTheory.volume a b)
    (hφ' : IntervalIntegrable φ' MeasureTheory.volume a b) (ha : φ a = 0) (hb : φ b = 0) :
    ∫ x in a..b, u x * φ' x = -∫ x in a..b, u' x * φ x := by
  rw [intervalIntegral.integral_mul_deriv_eq_deriv_mul hu hφ hu' hφ', ha, hb]
  ring

The one-dimensional maximum principle (subharmonic = convex), and Mathlib's mean value property for harmonic functions in the plane, from which the strong maximum principle follows:

/-- **Weak maximum principle in one dimension.** If `u'' ≥ 0` on `(a, b)` (`u` is
*subharmonic*: `−Δu ≤ 0`) and `u` is continuous on `[a, b]`, then `u ≤ max(u(a), u(b))` on
`[a, b]`: the maximum is attained on the boundary. In one dimension subharmonic means convex,
so this follows from Mathlib's `convexOn_of_deriv2_nonneg` and `ConvexOn.le_on_segment`. -/
theorem max_principle_1d {u : ℝ → ℝ} {a b : ℝ} (hab : a ≤ b)
    (hc : ContinuousOn u (Icc a b)) (hd : DifferentiableOn ℝ u (Ioo a b))
    (hd2 : DifferentiableOn ℝ (deriv u) (Ioo a b))
    (hsub : ∀ x ∈ Ioo a b, 0 ≤ deriv^[2] u x) {x : ℝ} (hx : x ∈ Icc a b) :
    u x ≤ max (u a) (u b) := by
  have hconv : ConvexOn ℝ (Icc a b) u :=
    convexOn_of_deriv2_nonneg (convex_Icc a b) hc (by rwa [interior_Icc])
      (by rwa [interior_Icc]) (by rwa [interior_Icc])
  exact hconv.le_on_segment (left_mem_Icc.2 hab) (right_mem_Icc.2 hab)
    (by rwa [segment_eq_Icc hab])

/-- **Mean value property.** A function harmonic on a neighbourhood of the closed disc of radius
`R` about `c` equals its average over the circle of radius `R` about `c`. (The maximum
principle for harmonic functions in the plane follows from it.) This is Mathlib's
`HarmonicOnNhd.circleAverage_eq`. -/
theorem mean_value_property {f : ℂ → ℝ} {c : ℂ} {R : ℝ}
    (hf : InnerProductSpace.HarmonicOnNhd f (Metric.closedBall c |R|)) :
    Real.circleAverage f c R = f c :=
  hf.circleAverage_eq

The blow-up in Hadamard's example — the quantitative heart of ill-posedness:

/-- **Hadamard's example.** `uₙ(x, y) = sin(n x) sinh(n y) / n²` solves `Δuₙ = 0` with Cauchy
data `uₙ(x, 0) = 0` and `∂_y uₙ(x, 0) = sin(n x)/n`, which tends to `0` uniformly. Yet at
`x = π/(2n)`, `y = 1` the solution equals `sinh(n)/n²`, which tends to `+∞`. Proof:
`sinh n ≥ (eⁿ − 1)/2` and `eⁿ/n² → ∞` (Mathlib's `Real.tendsto_exp_div_pow_atTop`). So the
solution does not depend continuously on the data. -/
theorem hadamard_blowup :
    Tendsto (fun n : ℕ => Real.sinh n / (n : ℝ) ^ 2) atTop atTop := by
  have h1 : Tendsto (fun x : ℝ => Real.exp x / x ^ 2 / 2 - 1 / 2) atTop atTop :=
    tendsto_atTop_add_const_right _ _
      ((Real.tendsto_exp_div_pow_atTop 2).atTop_div_const (by norm_num))
  refine tendsto_atTop_mono' atTop ?_ (h1.comp tendsto_natCast_atTop_atTop)
  filter_upwards [eventually_ge_atTop 1] with n hn
  have hn' : (1 : ℝ) ≤ n := by exact_mod_cast hn
  have hexp : Real.exp (-(n : ℝ)) ≤ 1 := Real.exp_le_one_iff.mpr (by linarith)
  have key : Real.exp (-(n : ℝ)) / (2 * (n : ℝ) ^ 2) ≤ 1 / 2 := by
    rw [div_le_div_iff₀ (by positivity) (by norm_num)]
    nlinarith
  have hsplit : Real.exp n / (n : ℝ) ^ 2 / 2 = Real.exp n / (2 * (n : ℝ) ^ 2) := by
    rw [div_div, mul_comm]
  show Real.exp (n : ℝ) / (n : ℝ) ^ 2 / 2 - 1 / 2 ≤ Real.sinh n / (n : ℝ) ^ 2
  rw [hsplit]
  calc Real.exp n / (2 * (n : ℝ) ^ 2) - 1 / 2
      ≤ Real.exp n / (2 * (n : ℝ) ^ 2) - Real.exp (-(n : ℝ)) / (2 * (n : ℝ) ^ 2) := by linarith
    _ = Real.sinh n / (n : ℝ) ^ 2 := by rw [Real.sinh_eq]; ring

What is and is not formalised. Mathlib does not yet have Sobolev spaces $H^1_0(\Omega)$ on domains, the Poincaré inequality, Rellich–Kondrachov or elliptic regularity (it does have distributions in development, Schwartz functions and tempered distributions, and harmonic functions). So our formalisation covers the abstract Hilbert-space core of the weak theory — which is all that Lax–Milgram and Céa's lemma use — and one-dimensional or planar instances of the other ideas. The application to $-\Delta u = f$ on a domain is proved in the text, not in Lean. We also do not formalise that Hadamard's $u_n$ are harmonic (a two-line computation with $\sin'' = -\sin$, $\sinh'' = \sinh$).

The experiment solves the weak form of $-u'' = \pi^2 \sin(\pi x)$ on $(0, 1)$ with zero boundary values, using continuous piecewise-linear "hat" functions on $N$ equal elements. The Galerkin equations are the tridiagonal system $\frac1h \operatorname{tridiag}(-1, 2, -1) U = F$. As $h = 1/N$ halves, the energy error halves and the $L^2$ error quarters, as Céa's lemma plus interpolation estimates (and a duality argument for $L^2$) predict. In one dimension something special happens: the nodal values are exact, because the Green's function of $-d^2/dx^2$ with pole at a node is itself piecewise linear. Indeed the inverse of the stiffness matrix is the Green's function sampled at the nodes, and its positivity is a discrete maximum principle.

"""Galerkin finite elements for -u'' = f on (0, 1), u(0) = u(1) = 0.

The weak formulation asks for u in H^1_0 with int u'v' = int f v for all v.
Galerkin's method imposes this only for v in the space of continuous
piecewise-linear "hat" functions on a uniform mesh of width h = 1/N, giving
the tridiagonal system (1/h) tridiag(-1, 2, -1) U = F. For f = pi^2 sin(pi x)
(exact solution sin(pi x)) we measure the errors: the energy error falls like
h, the L^2 error like h^2, and in one dimension the nodal values are exact.
We also check the discrete maximum principle: the inverse of the stiffness
matrix has positive entries, so f >= 0 gives U >= 0.
"""

import math

import numpy as np

GAUSS_X, GAUSS_W = np.polynomial.legendre.leggauss(6)

def stiffness(n):
    """The (n-1) x (n-1) matrix int phi_i' phi_j' = (1/h) tridiag(-1, 2, -1)."""
    if n < 2:
        raise ValueError("need at least 2 elements")
    h = 1 / n
    main = 2 * np.ones(n - 1)
    off = -np.ones(n - 2)
    return (np.diag(main) + np.diag(off, 1) + np.diag(off, -1)) / h

def load(n):
    """F_i = int f phi_i for f = pi^2 sin(pi x), computed exactly."""
    h = 1 / n
    nodes = np.arange(1, n) * h
    return 2 * (1 - math.cos(math.pi * h)) * np.sin(math.pi * nodes) / h

def solve(n):
    """Nodal values (including the zero boundary values) of the Galerkin
    solution with n elements."""
    interior = np.linalg.solve(stiffness(n), load(n))
    return np.concatenate([[0.0], interior, [0.0]])

def errors(n):
    """(energy error, L^2 error) of the Galerkin solution, by Gauss
    quadrature on each element; the exact solution is sin(pi x)."""
    values = solve(n)
    h = 1 / n
    energy = l2 = 0.0
    for i in range(n):
        a = i * h
        x = a + h * (GAUSS_X + 1) / 2
        w = GAUSS_W * h / 2
        slope = (values[i + 1] - values[i]) / h
        uh = values[i] + slope * (x - a)
        energy += np.sum(w * (math.pi * np.cos(math.pi * x) - slope) ** 2)
        l2 += np.sum(w * (np.sin(math.pi * x) - uh) ** 2)
    return math.sqrt(energy), math.sqrt(l2)

def nodal_error(n):
    """Largest error of the Galerkin solution at the mesh nodes."""
    nodes = np.arange(n + 1) / n
    return float(np.max(np.abs(solve(n) - np.sin(math.pi * nodes))))

def main():
    print("-u'' = pi^2 sin(pi x) on (0, 1), u(0) = u(1) = 0; exact u = sin(pi x)")
    print(f"{'N':>4} {'energy err':>11} {'ratio':>6} {'L2 err':>11} {'ratio':>6}"
          f" {'nodal err':>10}")
    previous = None
    for n in [4, 8, 16, 32, 64]:
        e1, e0 = errors(n)
        r1 = f"{previous[0] / e1:6.2f}" if previous else "     -"
        r0 = f"{previous[1] / e0:6.2f}" if previous else "     -"
        print(f"{n:>4} {e1:>11.3e} {r1} {e0:>11.3e} {r0} {nodal_error(n):>10.1e}")
        previous = (e1, e0)
    print("Halving h halves the energy error (rate 1) and quarters the L^2")
    print("error (rate 2); in 1D the nodal values are exact up to rounding.")
    print()
    inverse = np.linalg.inv(stiffness(16))
    print("Discrete maximum principle, N = 16:")
    print(f"  smallest entry of the inverse stiffness matrix: {inverse.min():.4f}")
    print("  all entries positive, so f >= 0 forces U >= 0")
    print()
    print("Epistemic status: numerical evidence for one smooth right-hand side")
    print("and five meshes. That the energy error is O(h) follows from Cea's")
    print("lemma (Lax-Milgram stability) plus interpolation estimates, the L^2")
    print("rate from a duality argument, and positivity from the M-matrix")
    print("structure; this table illustrates those theorems but proves none.")

if __name__ == "__main__":
    main()

This is evidence for one problem on a few meshes. The theorems — Lax–Milgram for existence and stability of both $u$ and $u_h$, Céa for the error — are what make the method trustworthy for problems where no exact solution is known.

The heat kernel $\frac{1}{\sqrt{4\pi t}} e^{-x^2/4t}$ is the solution of $u_t = u_{xx}$ starting from a unit point mass, the Dirac delta. Move $t$ towards $0$ and watch it concentrate (integral always $1$); move it up and watch it spread. For every $t > 0$ it is positive everywhere — infinite speed of propagation — and infinitely smooth: the heat equation smooths instantly.

Hadamard's example, along the line $y = 1$: the Cauchy data $\partial_y u_n(x, 0) = \sin(nx)/n$ shrink as $n$ grows, while the solution $u_n(x, 1) = \sin(nx)\sinh(n)/n^2$ explodes.