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

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

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

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

先说结论：三个人都对。答案分岔的地方，是把级数截断在有限项的时候。本文把两种表示放进同一段代码，看初始条件的光滑性如何改变误差的类型。谱方法和高阶格式里遇到的那个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$ 得到的函数。固定端条件不需要另行施加。奇延拓本身就在墙面处翻转了符号。没有展开，没有系数，没有频率。整个解就是把初始形状劈成两半，分别送往两侧。

在下面的模拟中亲手试一试。

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.

用 halves 按钮把两列行波关掉再打开，就能看出那条白线其实是两个半高副本的和。在 jump 形状下，注意拐角始终锋利如初。对流不会把任何东西抹平。

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

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

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 %

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

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

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

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

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 滑块推到最右端。在 jump 上，红色振荡只是变细，个头没变。右边的收敛曲线也在底部走平。同一个滑块换到 smooth 上再拖，曲线像悬崖一样掉下去。改变的只有初始条件。

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

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

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

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

先确认初始条件的光滑性

新格式出现振荡时，先看初始条件和边界数据，再去怀疑格式。初始场是不是按单元常数铺的？界面是不是当作跳跃塞进去的？入口型线在时间上是不是只有 $C^0$ ？只要中了一条，那个振荡就不是bug，而是表示方式的代价。

挑选验证算例时也是同一条标准。光滑解能干净地跑出设计精度。一旦引入跳跃，最大范数收敛就消失，只剩下 $L^1$ 。 1747年的三个人还没有描述这种区别的词汇。我们把它叫作范数的选择。