间断伽辽金法

间断伽辽金法（discontinuous Galerkin method，DG）是在每个网格单元内使用多项式、允许跨单元跳跃，并通过数值通量连接相邻单元的偏微分方程离散方法。它同时保留有限元的局部弱形式与有限体积法的通量守恒。单元内部可以提高多项式次数，而单元界面仍需要说明信息从哪一侧传来。

单元平均值之外再保存一个斜率

考虑周期区间上的线性输运 $u_t+a{u}_x=0$，其中 $a>0$。均匀单元宽为 $\Delta x$，在第 $i$ 个单元引入局部坐标 $\xi =2(x-x_i)/\Delta x\in [-1,1]$。一次 DG 近似写成 $$u_h{|}_{I_i}=U_i+V_i\xi .$$ $U_i$ 是单元平均值；$V_i$ 表示从中点到右端的变化量，实际空间斜率为 $2V_i/\Delta x$。左右端迹分别为 $U_i-V_i$ 和 $U_i+V_i$。两个相邻单元给同一界面的数值可以不同，这正是“间断”的含义。

图取两单元的 U=(0.8,0.3)、V=(0.3,−0.2)，所以界面左迹为 1.1，右迹为 0.5。向右输运时，界面通量采用左迹；后文另取分片常数初值演示时间导数。

弱形式与系数方程

在单元内以测试多项式 $v$ 乘方程并分部积分，得到 $${\int }_{I_i}(u_h{)}_{t}v\mathrm{d}x-a{\int }_{I_i}{u}_h{v}_x\mathrm{d}x+{\hat{f}}_{i+1/2}v(x_{i+1/2}^{-})-{\hat{f}}_{i-1/2}v(x_{i-1/2}^{+})=0.$$ 对向右输运，迎风通量是左侧迹乘以 $a$。记 $R_i=U_i+V_i$，便有 ${\hat{f}}_{i+1/2}=a{R}_i$、${\hat{f}}_{i-1/2}=a{R}_{i-1}$。用 $v=1$，体积分的第二项消失，得到 $${U_i}^{\prime }=-\frac{a}{\Delta x}(R_i-R_{i-1}).$$ 因此单元平均值仍按共享通量更新。

再用 $v=\xi$。由于 ${\int }_{-1}^1{\xi }^2\mathrm{d}\xi =2/3$、$v_x=2/\Delta x$，时间项为 $(\Delta x/3){V_i}^{\prime }$，体积项为 $-2a{U}_i$，两个边界项的和为 $a(R_i+R_{i-1})$。于是 $${V_i}^{\prime }=\frac{3a}{\Delta x}(2U_i-R_i-R_{i-1}).$$ 这两条方程完整规定了空间离散，尚未规定时间积分。更高次数会得到更大的单元局部质量矩阵和更多界面迹，而不是简单添加几个差分点。MFEM 的 DG 输运示例展示了弱形式、质量矩阵与时间求解器之间的分工。

一个可逐项计算的双单元例子

取两个周期单元、$a=\Delta x=1$，初值 $U=(1,0),V=(0,0)$。右端迹为 $R=(1,0)$，所以 $$U^{\prime }=(-1,1),V^{\prime }=(3,-3).$$ 如果只为观察初始变化，作一次 $\Delta t=0.05$ 的显式 Euler 试步，得到 $U=(0.95,0.05),V=(0.15,-0.15)$。第二单元的右端值为 $-0.1$，尽管平均值还是正的。这说明高阶多项式在不连续处可能振荡；单元平均值守恒并不保证多项式处处非负。正式计算需要相匹配的时间方法、步长限制，以及必要时的限制器。

有解变式取常数状态 $U=(2,2),V=(0,0)$。所有右端迹都等于二，两类导数均为零，常数被精确保留。若把第二条式子的上游迹误写成右邻，虽然常数测试仍可能通过，上面的非平凡例子会立即给出不同结果；两个测试应当一起核对。

守恒与半离散能量估计

将所有平均值方程乘以 $\Delta x$ 求和，周期通量逐项抵消，所以 $\sum _i\Delta x{U}_i$ 不变。进一步可以证明迎风通量具有耗散性。

定义 $E=\frac{1}{2}\sum _i{\int }_{I_i}{u}_h^2\mathrm{d}x$。在弱形式中取 $v=u_h$，并记每个单元左右迹为 $L_i,R_i$，则 $${E_i}^{\prime }=\frac{a}{2}(R_i^2-L_i^2)-a{R}_i^2+a{R}_{i-1}{L}_i.$$ 求和后，$\sum _i{R}_{i-1}^2=\sum _i{R}_i^2$，因此 $$E^{\prime }=-\frac{a}{2}\sum _i(R_{i-1}-L_i{)}^2\le 0.$$ 耗散恰好来自单元间的跳跃。若相邻迹连续，这一时刻没有界面耗散；出现跳跃时，迎风通量使离散能量下降。该证明针对连续时间的空间离散系统；任意时间步进不自动继承这个不等式。上例初始两处跳跃的平方均为一，因而 $E^{\prime }(0)=-1$。

可运行的局部算子

下面函数返回一次 DG 系数的时间导数，不自行选择时间步长。输入同长度的平均值和斜率系数，速度必须为正，边界为周期。返回两个新列表，不修改输入。$M$ 个单元只需 $O(M)$ 时间和空间；高阶、多维和变系数问题的成本还取决于局部基函数数目。

Python 3

import math


def dg_p1_rhs(averages, slopes, a, dx):
    U = [float(x) for x in averages]
    V = [float(x) for x in slopes]
    a, dx = float(a), float(dx)
    if not U or len(U) != len(V):
        raise ValueError("equal nonempty coefficient lists required")
    if not all(math.isfinite(x) for x in U + V + [a, dx]) or a <= 0 or dx <= 0:
        raise ValueError("finite coefficients and positive a, dx required")
    right = [u + v for u, v in zip(U, V)]
    dU = [-a / dx * (right[i] - right[i-1]) for i in range(len(U))]
    dV = [3 * a / dx * (2 * U[i] - right[i] - right[i-1])
          for i in range(len(U))]
    if not all(math.isfinite(x) for x in dU + dV):
        raise OverflowError("nonfinite derivative")
    return dU, dV


print(dg_p1_rhs([1, 0], [0, 0], 1, 1))  # ([-1.0, 1.0], [3.0, -3.0])
print(dg_p1_rhs([2, 2], [0, 0], 1, 1))  # zero derivatives

与有限体积法相比，这里保存的不只单元平均值，而是每个单元的一组多项式系数。与连续有限元相比，它不强制共享界面节点值，适合局部不同次数和非匹配网格；这些自由也使数值通量、惩罚项和稳定性条件成为方法的一部分。椭圆问题的内罚 DG 格式还需要另外定义通量与惩罚参数，不能用本篇的单向输运公式代替。