折れた弦は無事、切れた弦は9%跳ねた — 1747年の波動方程式論争の数値的な決着

解を何で表すかが誤差の種類を決めます。進行波は形を運び、モードの和は角で振動します。

1747年、一本の弦をめぐって三人が割れた#

ダランベールは1747年に振動する弦の方程式を書きました。人類最初の偏微分方程式です。ところが論争になったのは方程式そのものではありませんでした。解がどんな形をしてよいか、でした。オイラーとダニエル・ベルヌーイがそれぞれ別の答えを出し、三人は30年以上譲りませんでした。

結論から言えば三人とも正しかったのです。答えが分かれるのは級数を有限項で打ち切るときだけです。

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$ に延長したものです。固定端条件を別に課すことはありません。奇延長そのものが壁で符号を反転させてくれるからです。展開も係数も周波数もありません。初期形状を半分に割って両側へ運ぶ、それだけです。移流は何も鈍らせません。

ベルヌーイは同じ答えをsineで書き直した#

ダニエル・ベルヌーイの反論は音楽から来ました。弦は基音と倍音を同時に鳴らします。ならば解も定在波の重ね合わせであるはずだ、というわけです。

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}$ は存在しません。ダランベールはそもそもそんな曲線を解として認めたがりませんでした。一つの解析的な式で書ける曲線だけが関数だと考えていたからです。オイラーは、手で描いた任意の曲線も初期条件になりうると反論しました。

ベルヌーイはさらに踏み込みます。どんな曲線もsineの無限和で書ける、という主張でした。当時は根拠のない楽観に見えました。決着は1822年のフーリエまで待つことになります。線形の波は形を運ぶだけです。特性線が重なると衝撃波が生まれる非線形は、初期形状が滑らかでも自分で跳びを作ります。ダランベールの直観は、ここでは半分だけ正しかったのです。

Pythonで二つの表現を同じ時刻に並べた#

三つの初期変位を用意しました。跳びのある矩形、折れだけの三角形、滑らかな鐘形です。同じ時刻 $tc/L = 0.15$ でダランベール解と $N$ 項のsine級数を比べます。

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):
    """ベルヌーイ解 — sineモード 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倍するたびに誤差が4分の1になります。1次収束です。跳びがある場合だけオーバーシュートが消えず、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 倍のとき誤差が約 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年の三人は、この区別をノルムという言葉では持っていませんでした。