少取一个积分点，解就发散了 — DG 的求积下限与泰勒基

DG 的体积项是次数为 $2p-1$ 的多项式。高斯 $n$ 点公式精确到 $2n-1$，所以下限是 $n = p$；越过这条线，损失的不是一阶精度，而是整个格式。

少取一个积分点，答案整个消失了

我曾把一份间断伽辽金（DG, Discontinuous Galerkin — 在每个单元内放一个独立多项式、再用面通量把它们缝起来的高阶方法）代码的单元积分从高斯三点减到两点。按账面算，每个单元的积分开销少了三分之一。跑完一看，L2 误差在小数点后十三位都完全一致。于是我又贪心地降到一点。这次不是精度掉了一阶，而是解连一圈都没转完就发散了。

分界线不在网格尺寸上，也不在 CFL 数上。它在被积函数的多项式次数上。本文要处理的是：这条线究竟在哪里，为什么在那里，以及在任意网格上要守住它该如何选取基函数。依据是一个一维 DG-P2 求解器和一次质量矩阵条件数计算。

把 DG order p 调高，保持绿色的高斯点下限一起升高；把 Gauss points n 降到那条线以下，标识变红。

DG order p2Gauss points n2test mode i2

u_h shape0.60true 0.0000 · quad 0.0000

Set p = 2 and drag Gauss points from 3 down to 2: the badge stays green and the error bar stays empty, because the integrand only has degree 3. Drag to 1 and it turns red — no shape of u_h will bring it back. The minimum is n = p, and it is a property of the integrand’s degree, not of how fine the mesh is.

$n \ge p$ 时标识保持绿色，误差条保持为空，摇 u_h shape 也一样。 $n$ 再降一档，误差条就被填上。

Q1. DG 属于有限元还是有限体积#

两者都是。只看一个单元，它是有限元；只看单元边界，它是有限体积。

把守恒方程乘以试验函数 $\phi_i$ ，在单元 $\Omega_e$ 上积分并分部积分，得到：

\int_{\Omega_e} \frac{\partial U_h}{\partial t}\,\phi_i \, d\Omega - \int_{\Omega_e} \mathbf{F}(U_h)\cdot\nabla\phi_i \, d\Omega + \oint_{\partial\Omega_e} \phi_i\, \hat{\mathbf{F}}(U_h^-, U_h^+)\cdot\mathbf{n}\, dS = 0

其中 $U_h$ 是单元内的近似解， $\mathbf{F}$ 是对流通量， $\hat{\mathbf{F}}$ 是由两侧迹值 $U_h^-, U_h^+$ 构造的数值通量， $\mathbf{n}$ 是面法向。粘性通量和源项各自再加一项，但结构不变。

关键在于这个式子分成了两块。体积积分在单元内部闭合，只有面积分才与邻居通信，而进入其中的是黎曼求解器给出的单值通量。有限体积法用一个单元平均做的事，DG 只是改用多个多项式系数来做。单元里放多项式、只在面上用数值通量缝起来，这副骨架和非结构网格上的间断伽辽金离散是同一套。丢掉通量差分，激波速度会一起算错。

把近似解写成基函数的线性组合，

U_h(\mathbf{x}, t) = \sum_{k=1}^{K} U_k(t)\, b_k(\mathbf{x})

时间项就变成质量矩阵 $M_{ik} = \int_{\Omega_e} b_i b_k \, d\Omega$ 。每个单元一个 $K \times K$ 的小矩阵，且不与邻居耦合，因此可以逐单元预先求逆并存起来。这正是 DG 易于并行的重要原因。

Q2. 积分要精确到几次#

数一数体积积分被积函数的次数，答案自然出来。

设使用 $p$ 次多项式空间。 $U_h$ 是 $p$ 次，试验函数 $\phi_i$ 最高也是 $p$ 次，所以 $\nabla\phi_i$ 是 $p-1$ 次。若通量是线性的，乘积的次数为：

\deg\left(\mathbf{F}(U_h)\cdot\nabla\phi_i\right) = p + (p-1) = 2p - 1

高斯–勒让德 $n$ 点公式精确到 $2n-1$ 次。两式一拼，下限就出来了：

2n - 1 \ge 2p - 1 \quad \Longrightarrow \quad n \ge p

原始资料里"至少要做 $2K-1$ 阶的积分，否则精度阶会下降"说的就是这件事。网格再细，这个不等式也不会变，因为多项式次数与单元尺寸无关。

有两处需要留意。质量矩阵的被积函数是 $b_i b_k$ ，次数为 $2p$ ，它自己的下限高一档，是 $n \ge p+1$ 。另外，若通量非线性， $\mathbf{F}(U_h)$ 根本就不是多项式。这正是 Cockburn 与 Shu 建议体积用 $2p$ 次、面用 $2p+1$ 次的原因。工程中曲面单元还要乘上雅可比，因此余量留得更多。

Q3. 降到一点，究竟坏在哪里#

直接测比争论快。在周期区域 $[0, 2\pi]$ 上用 DG-P2 求解 $u_t + u_x = 0$ ：勒让德基，时间推进用 SSP-RK3，面通量取迎风。质量矩阵按解析式给出，这样体积积分的点数就成了唯一的变量。

from math import pi, sin, exp, log, sqrt, ceil
 
GAUSS = {                                   # Gauss-Legendre on [-1,1]: exact to degree 2n-1
    1: ([0.0], [2.0]),
    2: ([-0.5773502691896257, 0.5773502691896257], [1.0, 1.0]),
    3: ([-0.7745966692414834, 0.0, 0.7745966692414834], [5/9, 8/9, 5/9]),
    6: ([-0.9324695142031521, -0.6612093864662645, -0.2386191860831969,
          0.2386191860831969,  0.6612093864662645,  0.9324695142031521],
        [0.1713244923791704, 0.3607615730481386, 0.4679139345726910,
         0.4679139345726910, 0.3607615730481386, 0.1713244923791704]),
}
PHI  = [lambda s: 1.0, lambda s: s,   lambda s: 1.5*s*s - 0.5]   # 勒让德模态, p = 2
DPHI = [lambda s: 0.0, lambda s: 1.0, lambda s: 3.0*s]
K = 3
 
def dg_rhs(U, h, nq):
    """u_t + u_x = 0 的 DG 半离散残差，面上取迎风通量。"""
    xq, wq = GAUSS[nq]
    N = len(U)
    uR = [sum(U[j][i]*PHI[i](1.0) for i in range(K)) for j in range(N)]   # 右迹值
    R  = []
    for j in range(N):
        fR = uR[j]                       # a = 1 > 0，所以面取左单元的值
        fL = uR[j-1]
        row = []
        for i in range(K):
            vol = 0.0
            for xk, wk in zip(xq, wq):
                uh = sum(U[j][m]*PHI[m](xk) for m in range(K))
                vol += wk*DPHI[i](xk)*uh
            surf = PHI[i](1.0)*fR - PHI[i](-1.0)*fL
            row.append((vol - surf)*(2*i+1)/h)                # M_ii = h/(2i+1)
        R.append(row)
    return R
 
def run_dg(N, nq, T=1.0, cfl=0.05):
    h  = 2*pi/N
    xc = [h*(j + 0.5) for j in range(N)]
    xg, wg = GAUSS[6]
    u0 = lambda x: exp(sin(x))
    U  = [[(2*i+1)/2*sum(w*PHI[i](s)*u0(xc[j] + h/2*s) for s, w in zip(xg, wg))
           for i in range(K)] for j in range(N)]
    nt = int(ceil(T/(cfl*h/5))); dt = T/nt
    for _ in range(nt):                                        # SSP-RK3
        R0 = dg_rhs(U, h, nq)
        U1 = [[U[j][i] + dt*R0[j][i] for i in range(K)] for j in range(N)]
        R1 = dg_rhs(U1, h, nq)
        U2 = [[0.75*U[j][i] + 0.25*(U1[j][i] + dt*R1[j][i]) for i in range(K)] for j in range(N)]
        R2 = dg_rhs(U2, h, nq)
        U  = [[(U[j][i] + 2*(U2[j][i] + dt*R2[j][i]))/3 for i in range(K)] for j in range(N)]
    e2 = 0.0
    for j in range(N):
        for s, w in zip(xg, wg):
            uh = sum(U[j][i]*PHI[i](s) for i in range(K))
            e2 += w*(uh - u0(xc[j] + h/2*s - T))**2*h/2
    return sqrt(e2)
 
print("nq  exact-to-deg |   N=10       N=20       N=40    | order")
for nq in (1, 2, 3):
    e = [run_dg(N, nq) for N in (10, 20, 40)]
    print(f" {nq}       {2*nq-1}      | {e[0]:.3e}  {e[1]:.3e}  {e[2]:.3e} |  {log(e[1]/e[2], 2):.2f}")

nq  exact-to-deg |   N=10       N=20       N=40    | order
 1       1      | 1.170e+01  1.302e+01  1.186e+01 |  0.13
 2       3      | 5.989e-03  7.369e-04  9.211e-05 |  3.00
 3       5      | 5.989e-03  7.369e-04  9.211e-05 |  3.00

逐行来读。 $n=2$ 与 $n=3$ 在三套网格上打印出的位数完全相同。实际上它们在第十三位有效数字才分开，那点差别是舍入误差。原因是被积函数只有三次，两点公式给出的已经是精确值，多加点数没有任何收益。

$n=1$ 这一行性质完全不同。误差停在 $10^1$ 量级，网格加密四倍也不见缩小。收敛阶 0.13 不是"掉到一阶"，而是"不收敛"。欠积分（under-integration）每一步都喂给格式一个错误的体积项，这个误差会随时间放大。下面的模拟把这个过程原样呈现出来。

Gauss points n2cells N16play rate0.50

t = 0.00 · L2 0.00e+0

Watch the face jumps first: they are tiny while the rule is consistent, and they are what the upwind flux has to reconcile. Now drag Gauss points from 3 to 2 — nothing moves, the L2 readout does not budge. Drag to 1 and the parabolas tear apart within a fraction of a revolution. Adding cells only makes it happen sooner.

先确认把 Gauss points 从 3 降到 2 时 L2 误差读数纹丝不动。然后降到 1：每个单元的抛物线还没转完一圈就撕裂了。把 cells N 调大，崩得更快。

把输出的三行按次数拆开，下限就和表对上了。P2 的体积被积函数是 $2p-1=3$ 次。高斯 $n$ 点精确到 $2n-1$ ，所以 $n=2$ 覆盖 3 次， $n=3$ 覆盖 5 次。已经能积出 3 次的两点公式，和三点给出同一个积分值。因此 N=10、20、40 上的 L2 重叠在 $5.989\times 10^{-3}$ 、 $7.369\times 10^{-4}$ 、 $9.211\times 10^{-5}$ 。收敛阶是 $\log_2(7.369\times 10^{-4}/9.211\times 10^{-5})$ 。比值大约 8.00， $\log_2 8 = 3$ ，输出里的 3.00 就是这个除法。

$n=1$ 的精确次数是 1，够不到 $2p-1=3$ 。误差是 $1.170\times 10^{1}$ 、 $1.302\times 10^{1}$ 、 $1.186\times 10^{1}$ 。从 N=20 到 N=40 的比是 $1.302/1.186 \approx 1.098$ ， $\log_2 1.098 \approx 0.135$ ，和输出的 0.13 一致。网格加密四倍误差仍停在 10 附近，原因就是这个比。说到十三位都相同，指的是 $n=2$ 和 $n=3$ 之间的舍入， $n=1$ 则是另一个数量级。

条件数列也是同一种读法。 $h$ 为 1、 $10^{-1}$ 、 $10^{-2}$ 、 $10^{-3}$ 时，归一化前的值是 $7.225\times 10^{2}$ 、 $7.200\times 10^{6}$ 、 $7.200\times 10^{10}$ 、 $7.200\times 10^{14}$ 。长度每缩小十倍，条件数大约放大 $10^{4}$ 倍。 $7.200\times 10^{6}/7.225\times 10^{2} \approx 9.97\times 10^{3}$ ，所以只有第一档因为 7.225 而不是整 $10^{4}$ 。下面几档就是 $10^{4}$ 。 $p=2$ 时指数是 $2p=4$ ，与此相符。归一化那一列钉在 $7.225\times 10^{2}$ 。

Q4. 为什么偏偏用泰勒基#

前面之所以轻松，是因为只有一维。真实网格里四面体、六面体、棱柱、金字塔和多面体是混在一起的。标准有限元要把每种形状映射到参考单元，并在其上定义形函数——也就是说，那套雅可比工作，每种形状都要来一遍。而多面体压根没有参考单元。

Luo 等人提出的泰勒基跳过了映射，直接在单元中心 $\mathbf{x}_c$ 做泰勒展开：

U_h = \bar{U} + \left.\frac{\partial U}{\partial x}\right|_c (x - x_c) + \left.\frac{\partial U}{\partial y}\right|_c (y - y_c) + \left.\frac{\partial^2 U}{\partial x^2}\right|_c \frac{(x - x_c)^2}{2} + \cdots

把每一项各自的单元平均减掉，首项系数 $\bar{U}$ 就恰好是单元平均。这个性质在工程上分量很重。取 $p=0$ ，DG 与有限体积法完全重合，有限体积用的限制器可以直接搬过来。有限体积限制器进入 DG 代码走的就是这条通道。由于不区分单元形状，混合网格上一套代码就够了。

代价也有一个。若直接使用 $(x-x_c)^k$ ，质量矩阵元素按 $h^{k+l+1}$ 缩放，条件数随单元尺寸失控。边界层网格的 $h$ 大约是 $10^{-3}$ ，看看那时会怎样。

from math import factorial, sqrt
 
def taylor_mass(h, K, scale):
    """宽度为 h 的单元上，泰勒基 b_k = ((x-xc)/scale)^k / k! 的质量矩阵。"""
    M = [[0.0]*K for _ in range(K)]
    for i in range(K):
        for j in range(K):
            n = i + j
            if n % 2:                                   # 关于形心的奇数阶矩为零
                continue
            M[i][j] = (h/scale)**n * h / (2**n * (n+1) * factorial(i) * factorial(j))
    return M
 
def jacobi_eig(A, sweeps=60):
    """对称矩阵特征值 — 循环雅可比旋转。"""
    K = len(A); A = [row[:] for row in A]
    for _ in range(sweeps):
        for p in range(K-1):
            for q in range(p+1, K):
                if abs(A[p][q]) < 1e-300:
                    continue
                th = 0.5*(A[q][q]-A[p][p])/A[p][q]
                t  = (1 if th >= 0 else -1)/(abs(th)+sqrt(th*th+1))
                c  = 1/sqrt(t*t+1); s = t*c
                for k in range(K):
                    akp, akq = A[k][p], A[k][q]
                    A[k][p], A[k][q] = c*akp - s*akq, s*akp + c*akq
                for k in range(K):
                    apk, aqk = A[p][k], A[q][k]
                    A[p][k], A[q][k] = c*apk - s*aqk, s*apk + c*aqk
    return [A[k][k] for k in range(K)]
 
print(" h        raw Taylor      normalized")
for h in (1.0, 1e-1, 1e-2, 1e-3):
    out = []
    for scale in (1.0, h):
        ev = [abs(v) for v in jacobi_eig(taylor_mass(h, 3, scale))]
        out.append(max(ev)/min(ev))
    print(f" {h:<8.0e} {out[0]:.3e}       {out[1]:.3e}")

 h        raw Taylor      normalized
 1e+00    7.225e+02       7.225e+02
 1e-01    7.200e+06       7.225e+02
 1e-02    7.200e+10       7.225e+02
 1e-03    7.200e+14       7.225e+02

$h$ 每缩小十倍，条件数就放大 $10^4$ 倍；在 $p=2$ 时指数正是 $2p$ 。 $h = 10^{-3}$ 时达到 $7.2 \times 10^{14}$ ，几乎耗尽双精度 $10^{16}$ 的余量。右边按单元尺寸归一化的一列，则与 $h$ 无关地稳定在 722。差别全部来自单元内除以 $\Delta x$ 的那一行。到了 $p=3$ ，指数变成 6，不做归一化在任何实用网格上都无法使用。

Q5. 需要预先放进内存的是什么#

DG 代码的初始化阶段本质上就是在造表，顺序如下。

按形状给单元分类——四面体／六面体／棱柱／金字塔／多面体。
按形状给面分类——三角形／四边形／多边形。
为每种形状准备所需阶数的高斯求积规则。
在每个高斯点上计算基函数值及其梯度并存下来。

三维中 $p$ 次完全多项式空间的自由度是 $\binom{p+3}{3}$ 。

$p$	0	1	2	3	4
每单元模态数 $K$	1	4	10	20	35

原始资料里的 (1,4,10,20,35) 就是这一行，带 *3 的那组是每个模态的三个梯度分量。三维可压缩计算有五个守恒变量，所以在 $p=2$ 的六面体网格上，仅状态向量每单元就是 $5 \times 10 \times 8 = 400$ 字节，之后还要加上各高斯点处的基函数值。 $p=2$ 的六面体若体积求积取 $3^3 = 27$ 点，每单元还需再存 $27 \times 10$ 个实数。

这张表不必逐单元各存一份。参考坐标下的基函数值只要形状相同就完全一样，所以每种形状造一套即可，单元本身只需带上雅可比、形心和尺寸。只有多面体属于例外，必须自带一张表。

从 P1 升到 P2，账单落在哪里#

把 $p$ 从 1 升到 2，三维模态数从 4 变成 10，内存是 2.5 倍。体积点数下限跟着 $n \ge p$ 抬高，张量积是 $n^3$ ，点数变成八倍。稳定 CFL 大约是 $1/(2p+1)$ ，从 $p=1$ 的 $1/3$ 降到 $p=2$ 的 $1/5$ ，时间步缩成五分之三。若用泰勒基， $\Delta x^k$ 的指数变大，在 $h=10^{-3}$ 要躲开 $7.2\times 10^{14}$ 的条件数，归一化就不能省。

解在大范围内光滑时，误差按 $h^{p+1}$ 下降，提高 $p$ 比加密网格便宜。激波主导时，限制器会把这笔收益吃掉。事后怀疑一个单元，和 MOOD subcell 限制器在 DG 欧拉方程里做的事是同一个位置。唯一不划算的交易，是把点数省到 $n = p$ 以下。那条线下面，格式在解另一个方程。