格子を8倍細かくしても2.4674は動かなかった — 1807年のフーリエ級数が残した二つの定数

格子を細かくするとL2誤差は減りますが、鋸歯モードの減衰は常にπ²/4倍遅いままです。細分化では直せない誤差が別にあります。

格子を8倍細かくしたのに一つの数だけ動かなかった#

1次元の熱伝導を格子64点で解きました。続いて128、256、512と上げていきます。L2誤差は予想どおり順調に下がりました。

ところが表に並べて出していた数が一つだけ動きません。2.4674。四つの格子すべてで小数第4位まで同じでした。

格子間隔 $\Delta x$ を8分の1にしたのに、なぜある誤差はそのままなのか。この記事ではその数の正体を追います。答えは $\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}$ を入れると、2階微分が $-k^2 e^{ikx}$ になります。つまり $e^{ikx}$ はラプラシアンの固有関数で、固有値は $-k^2$ です。

そのため各モードは互いに干渉せず、単独で $e^{-\alpha k^2 t}$ として減っていきます。偏微分方程式がモードごとの常微分方程式に分解されるわけです。

指数の中の $k^2$ が減衰の速さを決めます。波長が半分なら減衰は4倍速い。短い波ほど先に消えます。拡散がすることは、この一箇所に集まっています。

格子は $-k^2$ の代わりに何を返すのか#

同じ計算を格子の上で行います。2階微分を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$ です。1格子上がって1格子下がる鋸歯の形で、このとき $\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):
    """中心差分による2階微分の実効波数の2乗（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$ まで押すと、格子側の波は厳密解より遅く減衰します。6次でもナイキスト比は0.6124付近で、減衰は1.633倍遅いままです。

1波長に4格子以上入る波（ $\theta/\pi \le 0.25$ ）は、次数を上げればほぼ完璧になります。鋸歯はそうなりません。6次に上げても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です。スキームは約束どおり2次です。

同じスキームに矩形波を与えると1.00に落ちます。不連続があるとフーリエ係数が $1/m$ でしか減らず、この遅い尾がエイリアシングとなって戻ってきます。収束次数を決めるのはスキームの次数ではなく解の滑らかさです。階段の脇の振動を止める仕事は、flux limiterで2次TVDを組む話と地続きです。

右端の比は8行のどの格子でも2.4674です。L2誤差は8分の1になったのに、この値は小数点以下まで同じです。

この二つを混ぜてはいけません。L2が下がるのは鋸歯モードの振幅が小さくなるからで、鋸歯の扱いが上手くなったからではありません。格子を細かくすることは問題をよりよく解くことではなく、問題を移すことです。

格子64では滑らかなコサインのL2が5.3683e-04です。8倍細かい512では8.3907e-06になります。比は約64倍で、対数を取ると次数は2.00のままです。矩形波は同じ区間で2.3383e-02から2.9168e-03へ落ちます。8倍、次数1.00です。鋸歯の減衰時間は格子64で厳密値2.474e-05、格子値6.104e-05です。512では3.865e-07と9.537e-07です。比は2.4674のまま、小数第4位まで動きません。ナイキストの厳密な $k^2$ は9.86960で、2次ステンシルが返すのは4.0000だけです。 $4/\pi^2 = 0.4053$ です。 $\Delta x$ が約分された式は、格子を細かくしても同じ数字を打ちます。

Gibbs側も同じ種類です。項16個の超過は9.01421%です。512個にすると8.94905%です。差は0.065%pほどです。 $\theta/\pi = 0.125$ の行で2次はすでに0.9872なので、1波長に格子16個が入れば1%以内です。6次はナイキストで6.0444を返し、比は0.6124、減衰は1.633倍遅いままです。次数を買って得るのは低波数の桁であり、鋸歯の運命ではありません。

格子64の矩形波L2、2.3383e-02を、滑らかな側の5.3683e-04で割ると約44倍です。同じ格子、同じスキームで、初期条件だけが違います。512では2.9168e-03対8.3907e-06、約348倍です。細かくするほど不連続側が相対的に大きくなります。次数1と次数2が格子比の1乗ぶん開くからです。512/64 = 8なので、誤差の比は約8倍広がります。44掛ける8は352で、表の348と合います。

6次でも0.6124である理由、そしてフィルタ#

次数を上げるのは低波数でだけ得をする商売です。表の $\theta/\pi = 0.125$ の行で6次は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と出て安心したことがあるなら、実効波数を一度描いてみる価値があります。スキームが鋸歯について何を告白するかは、そこにしか書かれていません。