把网格佩克莱数调到2，解掉到了-0.33 — 可最小化的能量消失的地方

一根杆上，四种方法给出了同一个答案

有限元讲义第1章把同一个杆件问题解了四遍：直接刚度法、最小总势能法、加权残差法，以及伽辽金法。四遍都得到同一个5×5刚度矩阵。讲义里这样写着：无论用哪一种，结果差别不大，都能得到准确值。

但这句话附带一个条件，而这个条件藏在杆件问题内部，不会浮到表面。流动方程不满足它。要量的是条件破裂时哪一项保证消失，以及稳定化还得回来什么、还不得回来什么。把 advection a 从0往上推，粉点绕一圈的功就从0.000变成恰好 $2\pi a$。

原先不横穿灰色等能量线、径直向内滑的黄球，会在 a 变大后画着螺旋往上爬。把 diffusion s 推到底，一圈的功仍然不动。

要有可最小化的量，矩阵必须对称

最小总势能法把总势能对节点未知量求导并令其为零。离散形式是这样的。

$\boldsymbol{\phi}$是节点未知量向量，$\mathbf{K}$是刚度矩阵，$\mathbf{f}$是载荷向量。现在对第$i$个分量求导。

注意出现的是$K_{ij} + K_{ji}$。要让这个式子变成$\mathbf{K}\boldsymbol{\phi} = \mathbf{f}$，就必须有$K_{ij} = K_{ji}$。不对称时，最小化实际求解的是对称化矩阵$\tfrac{1}{2}(\mathbf{K}+\mathbf{K}^{T})$，而那已经不是原来的方程。

同一处还有一个更根本的说法。残差场$\mathbf{r}(\boldsymbol{\phi}) = \mathbf{f} - \mathbf{K}\boldsymbol{\phi}$若要是某个标量函数的梯度，它的雅可比矩阵必须对称。不对称时，势函数根本不存在。检验是否存在最快的办法是沿闭合路径走一圈：梯度场绕一圈做的功必然为零。上面把 a 调高之后，小球逆着等能量线往上爬，就是这一圈不再回到零的地方。

对流项破坏的量恰好是$u$#

写出一维对流扩散方程$u\,\phi' = \epsilon\,\phi''$的弱形式，两项性格的差别就显出来了。

$w$是试函数，$\epsilon$是扩散系数，$u$是对流速度。第一项交换$w$和$\phi$后不变——它是对称的。后面那一项分部积分后符号翻转。对于在边界上为零的试函数，$\int w\,u\,\phi'\,dx = -\int \phi\,u\,w'\,dx$，所以对流项是纯粹的反对称成分。

用线性单元组装后，这个结构原样留在系数里。单元尺寸为$h$时，内部节点那一行是这样。

扩散往两侧邻居放同样的值，对流则一侧放$+u/2$，另一侧放$-u/2$。所以偏离对称的量恰好是一个数。

网格加密到什么程度，这个值都是$u$，因为$h$根本没进来。只要有流动，最小总势能原理就不会回来。刚度能折成6×6的Voigt矩阵，只发生在矩阵仍然对称的时候。对流在台阶附近把二阶格式震起来时，处方在通量限制器的TVD对流那一侧。

用Python在同一套网格上立起三种方法#

10个单元，$u = 1$，边界条件$\phi(0)=0$与$\phi(1)=1$。用一个人工扩散系数$\beta$生成三种格式：$\epsilon_{\text{eff}} = \epsilon + \beta\,u\,h/2$，$\beta = 0$是纯伽辽金，$\beta = 1$是完全迎风，$\beta = \coth(Pe_h) - 1/Pe_h$是SUPG。

import numpy as np
 
def assemble_ad(n, vel, eps, beta):
    """线性单元上的一维对流扩散。beta 为人工扩散系数。"""
    h = 1.0 / n
    eps_eff = eps + beta * vel * h / 2.0
    kd = (eps_eff / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])   # 扩散：对称
    ka = (vel / 2.0) * np.array([[-1.0, 1.0], [-1.0, 1.0]])     # 对流：反对称
    K = np.zeros((n + 1, n + 1))
    for e in range(n):
        K[e:e + 2, e:e + 2] += kd + ka
    return K
 
def solve_bvp(K):
    n = K.shape[0] - 1
    A, b = K.copy(), np.zeros(n + 1)
    A[0, :], A[0, 0], b[0] = 0.0, 1.0, 0.0
    A[n, :], A[n, n], b[n] = 0.0, 1.0, 1.0
    return np.linalg.solve(A, b)
 
def exact_ad(x, pe):
    return (np.exp(pe * (x - 1.0)) - np.exp(-pe)) / (1.0 - np.exp(-pe))
 
def skew_ratio(K):
    return np.linalg.norm(K - K.T) / np.linalg.norm(K + K.T)
 
def loop_work(K, m=20000):
    """残差场 -K.phi 沿自由度空间单位圆绕一圈所做的功。"""
    t = np.linspace(0.0, 2.0 * np.pi, m, endpoint=False)
    path = np.stack([np.cos(t), np.sin(t)])          # 位置
    tang = np.stack([-np.sin(t), np.cos(t)])         # dl / dt
    return float(np.sum(np.sum(-(K @ path) * tang, axis=0)) * (2.0 * np.pi / m))
 
n, vel = 10, 1.0
x = np.linspace(0.0, 1.0, n + 1)
print(f"{'Pe_h':>5} {'scheme':>9} {'max|err|%':>10} {'min phi':>9} {'skew/sym':>9} {'loop W':>8}")
for pe_h in [0.5, 1.0, 2.0, 5.0]:
    eps = vel / (2.0 * n * pe_h)
    ex = exact_ad(x, vel / eps)
    for name, beta in [("Galerkin", 0.0), ("upwind", 1.0),
                       ("SUPG", 1.0 / np.tanh(pe_h) - 1.0 / pe_h)]:
        K = assemble_ad(n, vel, eps, beta)
        phi = solve_bvp(K)
        print(f"{pe_h:5.1f} {name:>9} {100 * np.max(np.abs(phi - ex)):10.2f} "
              f"{phi.min():9.4f} {skew_ratio(K):9.4f} {loop_work(K[4:6, 4:6]):8.4f}")
 
K0 = assemble_ad(n, 0.0, 0.1, 0.0)
print(f"\nvel = 0 : skew/sym = {skew_ratio(K0):.2e},  loop W = {loop_work(K0[4:6, 4:6]):.2e}")
print(f"pi * u  = {np.pi * vel:.4f}")

 Pe_h    scheme  max|err|%   min phi  skew/sym   loop W
  0.5  Galerkin       3.45    0.0000    0.2924   3.1416
  0.5    upwind      13.17    0.0000    0.1954   3.1416
  0.5      SUPG       0.00   -0.0000    0.2704   3.1416
  1.0  Galerkin      13.53    0.0000    0.5774   3.1416
  1.0    upwind      19.80    0.0000    0.2924   3.1416
  1.0      SUPG       0.00   -0.0000    0.4428   3.1416
  2.0  Galerkin      35.17   -0.3334    1.1010   3.1416
  2.0    upwind      18.17   -0.0000    0.3885   3.1416
  2.0      SUPG       0.00    0.0000    0.5572   3.1416
  5.0  Galerkin      69.61   -0.6961    2.1517   3.1416
  5.0    upwind       9.09    0.0000    0.4836   3.1416
  5.0      SUPG       0.00    0.0000    0.5773   3.1416
 
vel = 0 : skew/sym = 0.00e+00,  loop W = -4.29e-16
pi * u  = 3.1416

把 Pe_h 滑块推过1，右侧 a_E 条变红越过，节点值沉到0以下。

beta 留在零、Pe_h 放到2时，能看到表里 min phi 那种 -0.3334 附近。按下 SUPG，误差贴在0 %，粉色的反对称那一行停在3.1416。

边界值被钉在0和1，$Pe_h = 2$的伽辽金解却掉到$-0.3334$。$Pe_h = 5$时是$-0.6961$。这是物理上不可能出现的值。

在$Pe_h = 1$处，一个系数改变了符号#

把网格佩克莱数定义为$Pe_h = u h / (2\epsilon)$，上面那一行右侧的系数就能重写成：

一旦$Pe_h > 1$，非对角元素变为正值。就在这一刻，矩阵不再是M矩阵，离散极值原理——内部解不会跑出边界值所构成区间的保证——也一同消失。表中min phi到$Pe_h = 1.0$为止都是零、到2.0变负，原因就在这里。上面 a_E 变红越过的那一刻，就是这个符号变化。

稳定化并不能找回对称性

这里有一种常见的误读。说迎风格式或SUPG"恢复了稳定性"，并不意味着对称性和最小化原理回来了。看表里的loop W列：三种格式、四个$Pe_h$值，一共12行全是3.1416。这个数就是$\pi u$，里面既没有$\epsilon$也没有$\beta$。

原因很简单。$\beta$放大的只是$\epsilon_{\text{eff}}$，而它属于矩阵的对称部分。反对称部分$u$没被碰过。skew/sym比值在$Pe_h = 5$时从2.1517降到0.4836，也不是反对称变小了，而是分母上的对称部分变大了。

$Pe_h = 0.5$ 时伽辽金的 max|err|% 是3.45，min phi 是0.0000，skew/sym 是0.2924。$Pe_h = 1.0$ 时误差升到13.53，最小值仍是0.0000。由 $K_{i,i+1} = (\epsilon/h)(Pe_h - 1)$，$Pe_h = 1$ 时右侧系数为零，一过1就变成正数。$Pe_h = 2$ 时伽辽金最小值为 -0.3334，误差35.17 %，skew/sym 为1.1010。同一行迎风格式最小值 -0.0000，误差18.17 %，skew/sym 0.3885。SUPG误差0.00，skew/sym 0.5572。三行的 loop W 全是3.1416。$u = 1$ 时 $\pi u$ 就是3.1416，所以 $\beta$ 取0、取1、或取 $\coth(Pe_h) - 1/Pe_h$，绕圈的功都一样。$Pe_h = 5$ 时伽辽金最小值 -0.6961，误差69.61 %，skew/sym 2.1517。迎风把误差降到9.09 %，把 skew/sym 降到0.4836。变小的是分母上的对称部分。速度为零的对照是 skew/sym 0.00e+00，loop W -4.29e-16。10个单元、$u = 1$ 时 $h = 0.1$，由 $Pe_h = uh/(2\epsilon)$，$Pe_h = 2$ 的 $\epsilon$ 是0.025。用 $\beta$ 把这个 $\epsilon$ 放大，$K_{i,i+1} - K_{i+1,i} = u$ 仍然是1。

所以稳定化真正买到的是更弱的保证。放弃最小化原理（能量范数下的最佳逼近），换来M矩阵性质和离散极值原理。代价用精度支付。$Pe_h = 2$时迎风格式消掉了振荡，但误差是18.17 %。而SUPG只加恰好需要的那点$\beta$，在节点上是精确的。

这个值在$Pe_h \to 0$时趋于0，在$Pe_h \to \infty$时趋于1：扩散占主导就关掉稳定化，对流占主导就走向完全迎风。不过节点精确性是一维常系数问题的特权。二维需要只沿流线方向加人工扩散的原始SUPG形式，那时"精确"也就没有了。

有限体积法在同一处使用的另一个名字

上面的计算用的是有限元的语言，但结论与离散方法无关。有限体积法里中心差分的对流项造出完全相同的模板，并在完全相同的$Pe_h = 1$处改变系数符号。为此登场的是一阶迎风，它带来的数值扩散正等于$\beta = 1$时的$u h / 2$。

两个世界分道扬镳的地方是守恒性。有限体积的迎风是改动面通量，总量守恒得以保持。这和守恒型与原始型在激波速度上分开，是同一种分岔。有限元的人工扩散是往刚度矩阵里加项，因此它到底守恒了什么要另外核对。

在加权残差这一族里换个选项，付的是另一种代价。最小二乘法最小化残差范数，因此总能造出对称正定矩阵——最小化原理回来了。代价是条件数被平方，而且在线性单元上，单元内部的二阶导数为零，扩散项会整个消失。基函数与积分规则会悄悄改变格式的性质。

$Pe_h$ 到了2，先买有界性再谈加密#

讲义里"用哪种方法都一样"，只在自伴算子里面成立。扩散、弹性、势流在这个范围里。对流一进来，$K_{i,i+1} - K_{i+1,i} = u$，对称就破了，loop W 不论 $\epsilon$ 怎么改都停在3.1416，也就是 $\pi u$。这时CG许诺的能量范数最佳逼近已经不在。非对称矩阵交给Krylov这一侧，从Jacobi到BiCGSTAB更老实。

min phi 为负，就是单元佩克莱数过了1。$Pe_h = 2$ 时伽辽金是 -0.3334，误差35.17 %。迎风消掉振荡，留下18.17 %的误差。SUPG的节点误差是0.00，skew/sym 只是从1.1010降到0.5572，反对称本身没动。三维里把 $Pe_h$ 从5降到1，单元数变成125倍。扛不住这个倍数时，买到的是有界性，不是精度。