折角的弦没事，切断的弦跳出了9% — 1747年波动方程之争的数值结局

用什么表示解，就决定了会得到哪一类误差。行波搬运形状，模态求和在拐角处振荡。

1747年，三个人为一根弦分成了三派#

达朗贝尔在1747年写下了振动弦的方程。这是人类的第一个偏微分方程。但争论的焦点不是方程本身，而是解可以长成什么样。欧拉和丹尼尔·伯努利各自给出了不同的答案，三个人三十多年互不相让。

先说结论：三个人都对。答案分岔的地方，是把级数截断在有限项的时候。

关掉 halves，白线拆成两列半高行波。换成 jump，拐角一直锋利。

speed0.35

Turn halves off and you see only the string. Turn it on and the same picture is two rigid copies at half height sliding in opposite directions — no frequencies anywhere. Switch to jump: the corners stay razor sharp forever, because transport never smooths anything. The clock reads t·c/L = 0.00; at 2 the shape is exactly back.

弦的形状就是这两半之和。拐角不被抹平的那条是达朗贝尔解。把模态叠上去的那条，后面会留下9%。

谱方法和高阶格式里那个9%振荡，源头就在这里。

达朗贝尔什么都没有展开

对张力为 $T$ 、线密度为 $\rho$ 的弦取一小段，套用牛顿第二定律，就得到下式。

u_{tt} = c^2 u_{xx}, \qquad c = \sqrt{T/\rho}

$u$ 是弦的横向位移， $c$ 是波速。达朗贝尔注意到，令 $\xi = x - ct$ 、 $\eta = x + ct$ 之后，方程变成 $u_{\xi\eta} = 0$ 。积分两次，解就直接出来了。

u(x,t) = \tfrac{1}{2}\left[\,F(x - ct) + F(x + ct)\,\right]

$F$ 是把初始位移 $f(x)$ 关于两端做奇延拓、周期取 $2L$ 得到的函数。固定端条件不需要另行施加。奇延拓本身就在墙面处翻转了符号。没有展开，没有系数，没有频率。整个解就是把初始形状劈成两半，分别送往两侧。对流不会把任何东西抹平。

伯努利用正弦重写了同一个答案

丹尼尔·伯努利的反驳来自音乐。一根弦同时发出基音和泛音。那么解也应该是驻波的叠加。

u(x,t) = \sum_{n=1}^{\infty} b_n \sin\!\frac{n\pi x}{L}\cos\!\frac{n\pi c t}{L}, \qquad b_n = \frac{2}{L}\int_0^L f(x)\,\sin\!\frac{n\pi x}{L}\,dx

$b_n$ 是第 $n$ 个模态的振幅， $n\pi c/L$ 是它的角频率。模态频率成 $1:2:3$ 的整数比，这一点从毕达哥拉斯起就为人所知。当时的巴黎正为拉莫的和声理论争论不休，达朗贝尔本人就站在用牛顿力学替它辩护的一边。音阶的物理依据，正是从这个方程里来的。

欧拉真正挑战的是"函数是什么"

欧拉站在达朗贝尔一边，理由却不同。用手指拨弦，初始形状是一个折角三角形。在折点处 $u_{xx}$ 并不存在。达朗贝尔根本不愿承认这种曲线是解，因为在他看来只有能用一个解析式写出的曲线才算函数。欧拉则坚持，手绘的任意曲线同样可以作初始条件。

伯努利走得更远：任何曲线都能写成正弦的无穷和。在当时这看起来像是毫无根据的乐观。这件事直到1822年傅里叶才有定论。线性波只搬运形状。特征线一旦相交就会自己长出激波的非线性问题即使初值光滑也会造出跳跃。达朗贝尔在这里的直觉只对了一半。

用Python把两种表示摆在同一时刻#

准备了三种初始位移：带跳跃的方波、只有折角的三角形、光滑的钟形。在同一时刻 $tc/L = 0.15$ 比较达朗贝尔解与 $N$ 项正弦级数。

import numpy as np
 
L, c, T = 1.0, 1.0, 0.15
X = np.linspace(0.0, L, 2001)
 
def hat_profile(x, a=0.35, b=0.55):
    """跳跃: 仅在 [a,b] 上为1，其余为0 — 锤子敲打的区间"""
    return np.where((x >= a) & (x <= b), 1.0, 0.0)
 
def kink_profile(x, a=0.45):
    """折角: 连续，但斜率在 x=a 处跳变 — 拨弦的形状"""
    return np.where(x < a, x / a, (L - x) / (L - a))
 
def bell_profile(x, a=0.45, s=0.055):
    """光滑: 无穷次可微的钟形"""
    return np.exp(-((x - a) ** 2) / (2 * s ** 2))
 
def odd_extend(f0, xq):
    """在固定端所要求的奇函数、周期2L的延拓上做插值"""
    xs = np.mod(xq, 2 * L)
    sign = np.where(xs > L, -1.0, 1.0)
    xs = np.where(xs > L, 2 * L - xs, xs)
    return sign * np.interp(xs, X, f0)
 
def dalembert_wave(f0, t):
    """达朗贝尔解 — 各取一半、分别向左右搬运的两列行波之和"""
    return 0.5 * (odd_extend(f0, X - c * t) + odd_extend(f0, X + c * t))
 
def modal_wave(f0, t, n_modes):
    """伯努利解 — n_modes 个正弦模态的叠加"""
    n = np.arange(1, n_modes + 1)[:, None]
    k = n * np.pi / L
    b = 2.0 / L * np.trapezoid(f0[None, :] * np.sin(k * X[None, :]), X, axis=1)
    return (b[:, None] * np.sin(k * X[None, :]) * np.cos(k * c * t)).sum(axis=0)
 
def overshoot_pct(u, exact):
    """相对精确解振幅的最大超出量 (%)"""
    return 100.0 * (u.max() - exact.max()) / (exact.max() - exact.min())
 
for name, f0 in (("jump  ", hat_profile(X)),
                 ("kink  ", kink_profile(X)),
                 ("smooth", bell_profile(X))):
    exact = dalembert_wave(f0, T)
    print(f"[{name}]  t*c/L = {T}")
    for n_modes in (8, 32, 128, 512):
        u = modal_wave(f0, T, n_modes)
        print(f"   N={n_modes:4d}   max|modal - dAlembert| = {np.abs(u - exact).max():.5f}"
              f"   overshoot = {overshoot_pct(u, exact):+6.2f} %")

[jump  ]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.30480   overshoot = +21.01 %
   N=  32   max|modal - dAlembert| = 0.26454   overshoot = +12.07 %
   N= 128   max|modal - dAlembert| = 0.23847   overshoot =  +9.86 %
   N= 512   max|modal - dAlembert| = 0.18738   overshoot =  +8.94 %
[kink  ]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.02073   overshoot =  -0.23 %
   N=  32   max|modal - dAlembert| = 0.00639   overshoot =  -0.19 %
   N= 128   max|modal - dAlembert| = 0.00157   overshoot =  -0.03 %
   N= 512   max|modal - dAlembert| = 0.00038   overshoot =  -0.01 %
[smooth]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.04435   overshoot =  -6.25 %
   N=  32   max|modal - dAlembert| = 0.00000   overshoot =  -0.00 %
   N= 128   max|modal - dAlembert| = 0.00000   overshoot =  +0.00 %
   N= 512   max|modal - dAlembert| = 0.00000   overshoot =  -0.00 %

modes N24speed0.25

Drag modes N to the right on jump: the red ringing narrows but its height parks near 9% (now 0.00%), and the right-hand curve flattens out. Switch to kink and then smooth with N untouched — the same slider that bought nothing now buys everything, and the error curve turns into a cliff (max error 0.0000).

把 modes N 从8拖到512。jump 的过冲停在9%附近，红色振荡只变细。smooth 到32项误差就灭了。kink 每把项数乘4，误差大约降到1/4。

三组数字讲了三个不同的故事。光滑钟形在32项就触到双精度的底。折角三角形每把 $N$ 提高4倍，误差就降到四分之一，这是一阶收敛。只有带跳跃的那一组，过冲不消失，而是停在9%附近。

跳跃的最大误差从 $N = 8$ 的 0.30480 只降到 $N = 512$ 的 0.18738。项数加到 64 倍，误差只缩到大约 1.6 倍。过冲从 +21.01% 落到 +8.94%，和理论值 8.95% 差 0.01 个百分点。同一段折角给出 0.02073、0.00639、0.00157、0.00038。比值大约是 3.24、4.07、4.13，接近 $N$ 乘 4 时误差降到 1/4 的一阶。钟形在 $N = 8$ 的误差是 0.04435，过冲 −6.25%，到 $N = 32$ 已是 0.00000。 $tc/L = 0.15$ ，行波只走了长度的 15%。截断级数的疤，在这个时刻已经坐在拐角上。

9%不会因为增加项数而变小#

这个数值就是吉布斯现象。1848年威尔布拉汉姆先发现， 1899年吉布斯重新确认，名字留了下来。理论值是跳跃高度的8.95%。上面512项给出的8.94%，就是这个常数。

过冲只是变窄，高度并不下降。项数增加时，振荡区间的宽度按 $1/N$ 收缩。所以用积分或 $L^2$ 范数来量，它是收敛的。用最大值来量则不收敛。这两种范数的差别，在工程上表现为负密度。改变的只有初始条件。

系数的衰减率解释了全部差别。有跳跃时 $b_n \sim n^{-1}$ ，只有折角时 $b_n \sim n^{-2}$ ，光滑时比任何幂次都衰减得快。被截掉的尾巴就是误差，尾巴越厚，切口留下的疤越大。

1755年，同一个人遇上非线性方程#

欧拉在1755年第一次把流体运动写成偏微分方程，也就是欧拉方程。与波动方程不同，这里特征线的斜率依赖于解本身。无论初始条件多么光滑，特征线一旦相交，有限时间内就会生成不连续。这与超声速流为什么不知道上游发生了什么是同一种结构。

所以在可压缩计算中，挑一个光滑的初始条件毫无用处。激波会自己造出跳跃，从那一刻起高阶格式就回到了1747年的问题。 von Neumann在1950年故意把激波抹开，针对的正是这种振荡。TVD限制器与WENO权重只在激波附近降阶，因为抹掉那9%的唯一办法，就是在局部把光滑性找回来。

用最大范数去量，阶数在跳跃处消失

新格式出现振荡时，不必先去翻限制器。先看初始场是不是单元常数，界面是不是当跳跃塞进去的，入口在时间上是不是只有 $C^0$ 。512 项过冲仍停在 +8.94% 的数据，振荡是表示方式的代价。

验证算例也按这一刀来选。钟形在 $N = 32$ 给出误差 0.00000，设计精度原样出来。跳跃一进来，最大范数收敛就消失，只剩下 $L^1$ 。1747 年的三个人还没有用范数这个词来分开这两种误差。