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

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

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

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

結論から言えば三人とも正しかったのです。答えが分かれるのは級数を有限項で打ち切るときだけです。この記事では二つの表現を同じコードに並べます。初期条件の滑らかさが誤差の種類をどう変えるかを見ます。スペクトル法や高次スキームで出会う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 では角が最後まで鋭いままです。移流は何も鈍らせません。

ベルヌーイは同じ答えを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 %

三つのブロックは別々の話をしています。滑らかな鐘形は32項で倍精度の底に届きます。折れた三角形は $N$ を4倍するたびに誤差が4分の1になります。1次収束です。跳びがある場合だけオーバーシュートが消えず、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$ でしかないかを確認します。一つでも当たれば、その振動はバグではありません。表現方式の代価です。

検証ケースを選ぶときも同じ基準が働きます。滑らかな解なら設計次数がそのまま出ます。跳びが入った瞬間に最大ノルムの収束は消え、 $L^1$ だけが残ります。 1747年の三人はこの区別を言葉として持っていませんでした。私たちはそれをノルムと呼びます。