网格加密8倍，2.4674纹丝不动 — 1807年傅里叶级数留下的两个常数

网格加密能让L2误差下降，但锯齿模态的衰减永远慢π²/4倍。有一类误差靠加密根本修不掉。

网格加密了8倍，只有一个数字纹丝不动#

我用64个网格解了一维热传导，接着换成128、256、512。L2误差如预期般稳步下降。

可是并排打印的另一个数字始终没动。2.4674。四套网格在小数点后四位上完全一致。

网格间距 $\Delta x$ 缩小到八分之一，为什么有的误差原地不动？本文追查这个数字的来历。答案是 $\pi^2/4$ ，起点在1807年的巴黎。

1807年，三位审稿人揪住的是收敛性#

傅里叶（Joseph Fourier）在1807年向巴黎科学院提交了一篇关于固体热传导的论文，主张固体内部的温度分布可以写成正弦与余弦之和。

审稿人是Lagrange、Laplace、Legendre三位。论文被退回。理由不是物理，而是数学：没有证明级数收敛。

Lagrange亲自尝试证明，未能完成，于1813年去世。这个空缺很久以后才由傅里叶的学生Dirichlet补上。

审稿人并非苛刻。他们看到的东西确实存在。把正弦级数套在含间断的函数上，无论加多少项，总有东西留下来。

宽0.5的方波上，滑块 N 是保留的项数。往上推，峰贴向间断，高度几乎不动。

terms N = 41Fejer averagezoom near jump

max overshoot = 0.000%  · peak at x = 0.0000 (jump at x = 0.25)  · distance to jump = -0.2500

把 N 调大，峰会贴紧间断。高度停在8.95%附近。打开 Fejer average，峰消失，台阶被抹糊。

用正弦堆出方波，剩下8.9490%#

在周期区间上放一个宽0.5的方波。傅里叶系数可以手写出来。

f(x) = \frac{1}{2} + \sum_{m=1}^{\infty} \frac{2}{\pi m}\sin\!\left(\frac{m\pi}{2}\right)\cos\big(2\pi m (x - 0.5)\big)

其中 $m$ 是项号， $x$ 是位置。只有奇数 $m$ 存活。把级数截断到 $N$ 项，部分和会在间断点旁边冲过1。

这就是Gibbs现象。极限值为8.9490%，是与项数无关的常数。这不是花精度就能买掉的误差。

热方程为什么偏爱正弦级数

傅里叶选正弦不是出于偏好。在热方程

\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}

中 $\alpha$ 是热扩散率。代入 $T = e^{ikx}$ ，二阶导数变成 $-k^2 e^{ikx}$ 。也就是说 $e^{ikx}$ 是拉普拉斯算子的本征函数，本征值为 $-k^2$ 。

于是各个模态互不干扰，各自按 $e^{-\alpha k^2 t}$ 衰减。一个偏微分方程被拆成每个模态一个的常微分方程。

指数里的 $k^2$ 决定衰减有多快。波长减半，衰减快四倍。短波先被抹掉。扩散要做的事都挤在这一处。

网格给出的不是 $-k^2$ ，那是什么#

现在把同样的计算搬到网格上。把二阶导数换成二阶中心差分，代入 $T_j = e^{ikj\Delta x}$ 。

\frac{T_{j+1} - 2T_j + T_{j-1}}{\Delta x^2} = -\frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right) e^{ikj\Delta x}

$j$ 是网格点编号。离散算子同样以 $e^{ikj\Delta x}$ 为本征函数，这意味着傅里叶分解在网格上照样成立。

变的是本征值。不是 $-k^2$ 而是 $-k'^2$ ，且

k'^2 = \frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right)

这个 $k'$ 叫有效波数（modified wavenumber）。它是格式在坦白"我以为的 $k$ 是这个"。

$\theta = k\Delta x$ 小的时候，有效曲线贴着 $\theta^2$ 。越往右越分开。同一个波，网格给出的那份死得更慢。

$\pi^2/4$ 出现在奈奎斯特极限的什么位置#

网格能表达的最短波长是 $2\Delta x$ ：上一格、下一格，一副锯齿。此时 $\theta = k\Delta x = \pi$ 。

代入上式：

\frac{k'^2}{k^2} = \frac{4\sin^2(\pi/2)}{\pi^2} = \frac{4}{\pi^2} = 0.4053

倒数是 $\pi^2/4 = 2.4674$ ，正是前面纹丝不动的那个数。

要紧的是 $\Delta x$ 被约掉了。这个比值不依赖网格间距。无论加密到多细，"该网格的锯齿"永远慢2.4674倍才死掉。

提高阶数也差不多。用真实模板验证一下。

import numpy as np
 
STENCILS = {
    2: [1.0, -2.0, 1.0],
    4: [-1/12, 4/3, -5/2, 4/3, -1/12],
    6: [1/90, -3/20, 3/2, -49/18, 3/2, -3/20, 1/90],
}
 
def modified_k2(order, theta):
    """中心差分二阶导数的有效波数平方（以 dx=1 为准）。精确值为 theta**2。"""
    c = STENCILS[order]
    h = len(c) // 2
    sym = sum(c[h + j] * np.cos(j * theta) for j in range(-h, h + 1))
    return -sym
 
theta = np.array([np.pi/8, np.pi/4, np.pi/2, 3*np.pi/4, np.pi])
print("theta/pi   exact k^2    order2    order4    order6    (ratio)")
for th in theta:
    ex = th**2
    row = [modified_k2(o, th) for o in (2, 4, 6)]
    print(f"{th/np.pi:6.3f}  {ex:10.5f}  " + "  ".join(f"{v:7.4f}" for v in row)
          + "   " + " ".join(f"{v/ex:6.4f}" for v in row))
 
print()
for o in (2, 4, 6):
    r = modified_k2(o, np.pi) / np.pi**2
    print(f"order {o}: Nyquist ratio = {r:.4f}   decay is {1/r:.3f}x slower")

theta/pi   exact k^2    order2    order4    order6    (ratio)
 0.125     0.15421   0.1522   0.1542   0.1542   0.9872 0.9997 1.0000
 0.250     0.61685   0.5858   0.6144   0.6166   0.9496 0.9960 0.9996
 0.500     2.46740   2.0000   2.3333   2.4222   0.8106 0.9457 0.9817
 0.750     5.55165   3.4142   4.3856   4.8278   0.6150 0.7900 0.8696
 1.000     9.86960   4.0000   5.3333   6.0444   0.4053 0.5404 0.6124
 
order 2: Nyquist ratio = 0.4053   decay is 2.467x slower
order 4: Nyquist ratio = 0.5404   decay is 1.851x slower
order 6: Nyquist ratio = 0.6124   decay is 1.633x slower

旋钮把表里的 $\theta$ 和差分阶数照原样转一遍。阶数升高时，低波数贴上 $\theta^2$ ，锯齿比只挪一点。

theta/pi = 1.00

k'^2 / k^2 = 0.4053  · grid decays 2.467x slower  · a*t/dx^2 = 0.000

把 $\theta$ 推到 $\pi$ ，网格那份比精确解死得慢。六阶的奈奎斯特比仍在0.6124附近，衰减慢1.633倍。

一个波长里放得下四格以上的波（ $\theta/\pi \le 0.25$ ），提高阶数后几乎完美。锯齿不是这样。升到六阶仍然是0.6124。

把同一个方波压到四套网格上

现在真的去解。为了不混进时间积分误差，时间方向用解析积分，只留下空间离散的差别。

import numpy as np
 
L, ALPHA, TEND = 1.0, 1.0, 2.0e-3
 
def tophat_coeffs(m):
    """周期 [0,1) 上 0.25 <= x < 0.75 方波的复傅里叶系数。"""
    m = np.asarray(m, dtype=float)
    out = np.zeros(m.shape, dtype=complex)
    out[m == 0] = 0.5
    mm = m[m != 0]
    out[m != 0] = (np.exp(-2j*np.pi*mm*0.75) - np.exp(-2j*np.pi*mm*0.25)) / (-2j*np.pi*mm)
    return out
 
def gibbs_overshoot(nterm):
    """2N+1 项部分和的最大超出量(%)。只在跳跃点附近密集采样。"""
    x = np.linspace(0.25, 0.25 + 4.0/nterm, 40001)
    m = np.arange(-nterm, nterm + 1)
    c = tophat_coeffs(m)
    s = (c[None, :] * np.exp(2j*np.pi*m[None, :]*x[:, None])).sum(axis=1).real
    return (s.max() - 1.0) * 100.0
 
def heat_error(ngrid, smooth):
    """热方程在时间方向精确积分，只保留空间离散化的差异。"""
    dx = L / ngrid
    x = np.arange(ngrid) * dx
    m = np.fft.fftfreq(ngrid, d=1.0/ngrid)
    k = 2*np.pi*m
    kmod2 = (4.0/dx**2) * np.sin(k*dx/2.0)**2          # 网格给出的 k^2
    mref = np.arange(-8192, 8193)
    if smooth:
        t0 = np.cos(2*np.pi*x) + 0.3*np.cos(6*np.pi*x)
        ref = (np.cos(2*np.pi*x)*np.exp(-ALPHA*(2*np.pi)**2*TEND)
               + 0.3*np.cos(6*np.pi*x)*np.exp(-ALPHA*(6*np.pi)**2*TEND))
    else:
        t0 = ((x >= 0.25) & (x < 0.75)).astype(float)
        c = tophat_coeffs(mref)*np.exp(-ALPHA*(2*np.pi*mref)**2*TEND)
        ref = (c[None, :]*np.exp(2j*np.pi*mref[None, :]*x[:, None])).sum(axis=1).real
    num = np.fft.ifft(np.fft.fft(t0)*np.exp(-ALPHA*kmod2*TEND)).real
    nyq = ngrid//2
    return np.sqrt(np.mean((num-ref)**2)), 1.0/(ALPHA*k[nyq]**2), 1.0/(ALPHA*kmod2[nyq])
 
print("N_terms   Gibbs overshoot[%]")
for n in (16, 32, 64, 128, 256, 512):
    print(f"{n:6d}    {gibbs_overshoot(n):9.5f}")
 
for smooth in (True, False):
    print(f"\ninitial condition: {'smooth cosine' if smooth else 'square wave (discontinuous)'}")
    print("grid      L2 error     order     sawtooth decay time exact/grid    ratio")
    prev = None
    for ng in (64, 128, 256, 512):
        l2, t_ex, t_gr = heat_error(ng, smooth)
        order = "  -  " if prev is None else f"{np.log2(prev/l2):5.2f}"
        print(f"{ng:4d}   {l2:.4e}    {order}     {t_ex:.3e} / {t_gr:.3e}   {t_gr/t_ex:.4f}")
        prev = l2

N_terms   Gibbs overshoot[%]
    16      9.01421
    32      8.96527
    64      8.95306
   128      8.95000
   256      8.94924
   512      8.94905
 
initial condition: smooth cosine
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   5.3683e-04      -       2.474e-05 / 6.104e-05   2.4674
 128   1.3424e-04     2.00     6.184e-06 / 1.526e-05   2.4674
 256   3.3562e-05     2.00     1.546e-06 / 3.815e-06   2.4674
 512   8.3907e-06     2.00     3.865e-07 / 9.537e-07   2.4674
 
initial condition: square wave (discontinuous)
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   2.3383e-02      -       2.474e-05 / 6.104e-05   2.4674
 128   1.1673e-02     1.00     6.184e-06 / 1.526e-05   2.4674
 256   5.8342e-03     1.00     1.546e-06 / 3.815e-06   2.4674
 512   2.9168e-03     1.00     3.865e-07 / 9.537e-07   2.4674

加密会减小的误差，和加密也不减的误差

光滑初始条件下收敛阶正好是2.00。格式如约是二阶的。

把方波交给同一个格式，收敛阶掉到1.00。存在间断时傅里叶系数只按 $1/m$ 衰减，这条慢尾巴以混叠的形式回来。决定收敛阶的是解的光滑性，不是格式的阶数。压掉台阶旁边的振荡，和用flux limiter搭出二阶TVD是同一类选择。

最右列的比值在八行里都是2.4674。L2误差降了8倍，这个值连小数位都一样。

这两件事不能混在一起。L2下降是因为锯齿模态的幅值变小了，不是因为格式处理锯齿的本事变强了。加密网格不是把问题解得更好，而是把问题挪了个地方。

64个网格上，光滑余弦的L2是5.3683e-04。加密到512，也就是八倍，变成8.3907e-06。相除大约64倍，取对数仍是2.00阶。方波在同一段网格上从2.3383e-02降到2.9168e-03，八倍，阶数1.00。锯齿衰减时间在64个网格上，精确值2.474e-05，网格值6.104e-05。512上是3.865e-07和9.537e-07。比值停在2.4674，小数第四位都不动。奈奎斯特处精确的 $k^2$ 是9.86960，二阶模板只返回4.0000。 $4/\pi^2 = 0.4053$ 。 $\Delta x$ 被约掉之后，再加密也只能打出同一个数。

Gibbs也是同一类。16项时超出是9.01421%。加到512项是8.94905%。差大约0.065个百分点。 $\theta/\pi = 0.125$ 那一行，二阶已经是0.9872，一个波长里放进16个格子就进了1%以内。六阶在奈奎斯特给出6.0444，比值0.6124，衰减慢1.633倍。买阶数买到的是低波数的位数，不是锯齿的命运。

把64个网格上方波的L2 2.3383e-02除以光滑侧的5.3683e-04，大约44倍。同一套网格、同一个格式，只换了初值。512上是2.9168e-03对8.3907e-06，大约348倍。越加密，间断问题相对越大。一阶和二阶之间隔着网格比的一次方。512/64 = 8，误差比大约再拉开8倍。44乘8是352，和表里的348对得上。

为什么六阶还是0.6124，以及滤波#

提高阶数只在低波数上划算。表中 $\theta/\pi = 0.125$ 那一行六阶给出1.0000，而在奈奎斯特处是0.6124。任何有限模板都处理不好锯齿，因为 $2\Delta x$ 的波在网格上连求导本身都是含糊的。

于是工程实践分成两路。一路是不让锯齿产生。交错网格和Rhie–Chow插值干的就是这件事。只加密误差大的格子，那个格子上的锯齿比仍是 $\pi^2/4$ 。让网格自己长出来的四叉树也抹不掉这个比。

另一路是把已经产生的锯齿显式抹掉。谱滤波、人工粘性、WENO的权重都属于这一类。上面viz里的 Fejer average 是其中最简单的形式，代价也在同一张图里：Gibbs峰没了，台阶糊了。

两条路都不免费。能选的不是误差消不消失，而是把它放在哪里。

Dirichlet在1829年填上的那一格#

Lagrange的反对不是形式逻辑。他怀疑的是级数遇到间断时会发生什么，而那个位置上确实有8.9490%。

Dirichlet在1829年证明了条件收敛，把那一格填上：逐点收敛，但不一致收敛。Gibbs峰从来不是反例，它就是"不一致"的样子。

今天跑网格收敛测试时，我们在重复同一个问题。这个误差属于加密就会减小的那种，还是像 $\pi^2/4$ 那样网格间距被约掉、怎么都不消失的那种？收敛阶那一行只回答前一种。

如果你曾因为收敛阶打出2.00而放心，那么把有效波数画一次是值得的。格式对锯齿的坦白，只写在那里。