Date: February 19, 2025
Topics: Discrete Lagrangian mechanics, discrete integration by parts, Noether constraints for stencils, polynomial reproduction via symmetry
This lecture completes the bridge between classical mechanics and numerical methods by applying the discrete action principle to spatial discretization of PDEs. We previously developed the Lagrangian formulation in continuous space-time; now we discretize the spatial domain while keeping time continuous, deriving finite difference stencils from variational principles.
The key innovation is recognizing that discrete stencils must respect discrete analogs of continuous symmetries. Just as Noether's theorem connects continuous symmetries to conserved quantities, discrete symmetries impose algebraic constraints on stencil coefficients. For the wave equation, translation invariance requires stencils to sum to zero, while polynomial reproduction (the discrete analog of Taylor consistency) imposes additional constraints on higher-order terms.
We introduce discrete integration by parts using shift operators, which elegantly mirrors the continuous integration by parts used to derive Euler-Lagrange equations. This leads to a systematic framework: discretize the Lagrangian, apply discrete variational calculus, and obtain stencil weights from polynomial reproduction constraints. The resulting stencils automatically inherit symmetry properties from the continuous Lagrangian.
This approach unifies our previous ad-hoc stencil constructions (centered differences, compact schemes) under a single variational principle. It also reveals a deep connection: searching for minimal-norm stencils with polynomial reproduction (Lecture 6) is equivalent to applying the discrete action principle with appropriate symmetry constraints.
Pedagogical significance: This material demonstrates how mathematical physics principles (symmetries, variational principles) translate directly into numerical algorithms. Understanding this connection is essential for designing structure-preserving numerical methods and physics-informed neural networks that respect conservation laws at the discrete level.
Recall the action for the 1D wave equation:
$$S = \int_0^T \int_0^L \left[\frac{1}{2}(\partial_t u)^2 - \frac{c^2}{2}(\partial_x u)^2\right] dx\, dt$$where:
Define piecewise constant extension of a grid function $u_i(t)$ on nodes $x_i = ih$ for $i = 1, \ldots, N$:
$$\phi_i(x) = \mathbb{I}_{(x_i - h/2, x_i + h/2)}(x) = \begin{cases} 1 & \text{if } x \in (x_i - h/2, x_i + h/2) \\ 0 & \text{else} \end{cases}$$Discrete representation:
$$u(x, t) = \sum_{i=1}^N u_i(t) \phi_i(x)$$Key observation:
$$\int_0^L u(x, t)\, dx = \sum_{i=1}^N h u_i(t)$$Each basis function $\phi_i$ contributes area $h$ (trapezoid rule).
Selecting $q_i = u_i$ as generalized coordinates:
$$L = \sum_{i=1}^N \left[\frac{1}{2}\dot{q}_i^2 h - \frac{c^2}{2}(D_h q_i)^2 h\right]$$where $D_h q_i$ is the spatial finite difference stencil.
Important: In light of Noether's theorem, we analyze symmetries:
(1) Time translation invariance:
Time only enters via $\dot{q}_i$, so the Lagrangian is time-translation invariant.
Conserved quantity (Hamiltonian):
$$H = \partial_{\dot{q}} L \cdot \dot{q} - L$$where the conjugate momentum is:
$$p_i = \frac{\partial L}{\partial \dot{q}_i} = \dot{q}_i h$$(2) Spatial translation invariance:
$q \to q + \delta q$ (shift all grid values by a constant) is a symmetry only if $D_h \delta q = 0$ for constant $\delta q$.
Conserved quantity (total momentum):
$$\text{const.} = \sum_{i=1}^N p_i = \sum_{i=1}^N h\dot{q}_i$$This is the total discrete momentum (integral of velocity).
Critical constraint: Translation invariance requires $D_h(\text{const.}) = 0$.
Lemma: Consider two periodic grid functions $u_i, v_i$ satisfying $u_0 = u_N, v_0 = v_N$.
Then:
$$\sum_{i=1}^N x_i y_{i+\alpha} = \sum_{i=1}^N x_{i-\alpha} y_i$$Compact notation using shift operators:
Define $E^\alpha u_i = u_{i+\alpha}$ (forward shift by $\alpha$). Then:
$$\langle x, E^\alpha y \rangle = \langle E^{-\alpha} x, y \rangle$$where $\langle \cdot, \cdot \rangle$ denotes the discrete inner product $\sum_i x_i y_i$.
Consider a stencil operator:
$$D_h u_i = \sum_{k=-m}^m c_k E^k u_i$$Definition: The adjoint operator is:
$$D_h^* u_i = \sum_{k=-m}^m c_{-k} E^k u_i$$(flip stencil coefficients: $c_k \to c_{-k}$).
Theorem:
$$\langle D_h^* u, \delta u \rangle = \langle u, D_h \delta u \rangle$$This is the discrete analog of integration by parts!
Continuous:
$$F[u] = \frac{1}{2}\int |\nabla u|^2\, dx$$Step 1: Take variation:
$$(\delta_u F, \delta u) = \int \nabla u \cdot \nabla \delta u\, dx$$Step 2: Integration by parts:
$$(\delta_u F, \delta u) = -\int (\nabla^2 u) \delta u\, dx = -(\nabla^2 u, \delta u)$$Result:
$$\delta_u F = -\nabla^2 u$$Discrete functional:
$$F_h[u] = \frac{1}{2}\int \left(\sum_i D_0 u_i \phi_i(x)\right)^2 dx = \sum_i \frac{h}{2}(D_0 u_i)^2$$where $D_0 u_i = \frac{u_{i+1} - u_{i-1}}{2h}$ (centered difference).
Step 1: Take variation:
$$\langle \delta_{u_i} F_h[u], \delta u_i \rangle = \langle D_0 u, D_0 \delta u \rangle$$Step 2: Apply discrete integration by parts:
$$\langle D_0 u, D_0 \delta u \rangle = \langle D_0^* D_0 u, \delta u \rangle$$Result:
$$\delta_{u_i} F_h = D_0^* D_0 u$$This is the discrete Laplacian!
Variation of $S_1$ (kinetic term):
$$\delta_{q_i} S_1 = -h\ddot{q}_i$$Variation of $S_2$ (potential term):
$$\delta_{q_i} S_2 = -c^2 h \cdot D_h^* D_h q_i$$Euler-Lagrange equation:
$$\ddot{q}_i = -c^2 D_h^* D_h q_i$$Remark: While tedious, no thought was needed to reach this point — just careful application of calculus!
We now have two constraints for a good stencil:
(1) Noether constraint (translation invariance):
$$D_h(\text{const.}) = 0 \quad \Rightarrow \quad \sum_k c_k = 0$$Derivation:
$$D_h \delta q = \sum_k c_k E^k \delta q = \delta q \sum_k c_k$$If $\delta q = \text{const.}$, then $D_h \delta q = 0$ requires $\sum_k c_k = 0$.
(2) Polynomial reproduction:
Does $D_h^* D_h q$ give a stable prediction of $\nabla^2 u$?
Recall from earlier lectures: this is true if $D_h^* D_h p = \nabla^2 p$ for any quadratic polynomial $p$.
General stencil:
$$D_h u_i = \frac{1}{h}\sum_j c_j u_{i+j}$$Composed operator:
$$D_h^* D_h u_i = \frac{1}{h^2} \sum_{j,k} c_j c_{-k} u_{i+j+k}$$Set:
$$c_{-1} = \alpha, \quad c_0 = \beta, \quad c_1 = \gamma$$Build table of products $c_j c_{-k}$ for indices $i+j+k$:
| $k \backslash j$ | -1 | 0 | 1 |
|---|---|---|---|
| -1 | $\alpha^2$ (i-2) | $\alpha\beta$ (i-1) | $\alpha\gamma$ (i) |
| 0 | $\beta\alpha$ (i-1) | $\beta^2$ (i) | $\beta\gamma$ (i+1) |
| 1 | $\gamma\alpha$ (i) | $\gamma\beta$ (i+1) | $\gamma^2$ (i+2) |
Sum along diagonals:
$$D_h^* D_h u_i = \frac{1}{h^2}\left[\alpha\gamma u_{i-2} + \beta(\alpha + \gamma) u_{i-1} + (\alpha^2 + \beta^2 + \gamma^2) u_i + \beta(\alpha + \gamma) u_{i+1} + \alpha\gamma u_{i+2}\right]$$(1) Noether constraint (constant reproduction):
$$\alpha + \beta + \gamma = 0$$(2) Constant reproduction (apply to $u_i = 1$):
$$2\alpha\gamma + 2\beta(\alpha + \gamma) + \alpha^2 + \beta^2 + \gamma^2 = 0$$(3) Linear reproduction (apply to $u_i = ih$):
This is automatically satisfied (terms cancel).
(4) Quadratic reproduction (apply to $u_i = (ih)^2$, want $\nabla^2 u = 2$):
$$8\alpha\gamma + 2\beta(\alpha + \gamma) = -2$$Solving (1)–(4) is non-unique. Let's try some simplifications:
Set $\gamma = 0$:
From (1): $\alpha = -\beta$
From (4): $-2\beta^2 = -2 \Rightarrow \beta = \pm 1$
Taking $\beta = 1$:
$$\alpha = -1, \quad \beta = 2, \quad \gamma = -1$$ $$D_h u = \frac{1}{h^2}(-u_{i-1} + 2u_i - u_{i+1})$$We recover the 3-point Laplacian $D_+ D_- = D_- D_+$! ✓
Set $\alpha = \gamma$:
From (1): $\alpha = -\beta/2$
From (4): $8\alpha^2 + 2\beta \cdot 2\alpha = -2$
Substitute $\beta = -2\alpha$:
$$8\alpha^2 - 8\alpha^2 = -2 \quad \Rightarrow \quad 0 = -2$$No real-valued solutions! Symmetric 3-point stencils cannot satisfy polynomial reproduction up to degree 2.
Set $\alpha = -\gamma$:
From (1): $\beta = 0$
From (4): $-8\alpha^2 = -2 \Rightarrow \alpha^2 = 1/4 \Rightarrow \alpha = \pm 1/2$
Taking $\alpha = 1/2$:
$$\alpha = \frac{1}{2}, \quad \beta = 0, \quad \gamma = -\frac{1}{2}$$ $$D_h u = \frac{1}{h^2}\left[-\frac{1}{4}u_{i-2} + \frac{1}{2}u_i - \frac{1}{4}u_{i+2}\right]$$This is a 5-point compact stencil using only nodes $i-2, i, i+2$.
This lecture unified variational principles with numerical analysis to derive finite difference stencils from symmetry and polynomial reproduction:
Key Takeaway: Finite difference stencils are not arbitrary — they emerge naturally from applying the discrete action principle with appropriate symmetry constraints. The variational framework guarantees that discrete schemes inherit conservation laws from continuous physics. Noether's theorem provides the algebraic constraints (e.g., $\sum c_k = 0$), while polynomial reproduction ensures accuracy. This unifies our previous ad-hoc constructions under a single principle and provides a blueprint for designing novel structure-preserving schemes.
Connection to machine learning: When learning stencils from data, we should encode these symmetry constraints (translation invariance, polynomial reproduction) as architectural inductive biases. This reduces the hypothesis space, improves sample efficiency, and guarantees physically meaningful solutions.