差分法の整合性・安定性・収束：Lax の同等定理と von Neumann 解析

0. この記事の要点

波や熱の広がりを計算機で追うには、時間と空間を細かい格子に区切り、微分を引き算と割り算に置き換えます。ところが、格子をどれだけ細かくしても、置き換え方によっては計算結果が厳密解に近づかず、かえって発散してしまいます。この記事は、「どんな置き換えなら格子を細かくするほど正しい答えに近づくのか」を、整合性と安定性という 2 つの性質で見分ける方法を示します。読み終えると、与えられた差分スキームが使ってよいものかどうかを、計算で判定できるようになります。

• 差分スキーム $u^{n+1} = C(\Delta t)u^n$ の誤差は $e^{n+1} = C(\Delta t)e^n + \Delta t\,\tau^n$ という漸化式を満たします。$\tau^n$ は厳密解をスキームに代入したときの残差（局所打ち切り誤差）です。冪 $C(\Delta t)^n$ が一様に有界（安定）なら、誤差は打ち切り誤差の $KT$ 倍以下に収まります。
• Lax の同等定理: 適切な線形初期値問題に対する整合的なスキームについて、収束することと安定であることは同値です。「安定 $\Rightarrow$ 収束」は誤差の漸化式から、「収束 $\Rightarrow$ 安定」は一様有界性原理（作用素の列が各点で有界なら、ノルムも有界になるという定理）から出ます。
• 定数係数のスカラーのスキームは、半離散 Fourier 変換で増幅因子 $g(\theta)$ の掛け算になります。$\ell^2$ で安定であることと $\max_\theta|g| \le 1 + C\Delta t$ は同値です（von Neumann）。連立のスキームでは、増幅行列のスペクトル半径（固有値の絶対値の最大）についての同じ条件は必要条件にとどまります。
• 陽的スキームが収束するには、数値的な依存領域が方程式の依存領域を含まなければなりません（CFL 条件）。移流方程式 $u_t + au_x = 0$ の 3 点スキームでは $|a|\Delta t/h \le 1$ です。これは十分条件ではなく、前進時間・中心空間のスキームは CFL 条件を満たしても不安定です。
• 風上差分と Lax–Friedrichs は 1 次精度で、誤差の先頭項は拡散（数値粘性。計算のせいで生じる見かけの拡散）です。Lax–Wendroff は 2 次精度で、誤差の先頭項は 3 階微分（分散）なので不連続の近くで振動します。単調性を保つ線形スキームは高々 1 次精度です（Godunov）。

この記事で使う事前知識

• Banach 空間: ノルムが定まり、コーシー列がすべて収束する線形空間。→ Banach 空間
• スペクトル半径: 行列の固有値の絶対値の最大。→ スペクトル半径
• ユニタリ写像: 内積（したがってノルム）を保つ、全射の線形写像。→ ユニタリ写像
• 正規行列: 自分の随伴行列と交換できる行列。ユニタリ行列に対角化でき、冪のノルムがスペクトル半径の冪になります。→ 正規行列
• Fourier 係数: 周期関数を三角関数の和に展開したときの各項の係数。→ Fourier 係数

1. 動機

偏微分方程式 $u_t = Lu$ を計算機で解くには、微分を差分商に置き換えます。置き換えたときの誤差は Taylor の定理で分かります。たとえば前進差分 $(u(x+h) - u(x))/h$ は $u_x$ を $O(h)$ の誤差で近似します。それなら、どの差分商を使っても、$h$ と $\Delta t$ を小さくすれば解は厳密解に近づくはずだ、と考えたくなります。

実際にはそうなりません。移流方程式 $u_t + u_x = 0$ を、時間は前進差分、空間は中心差分で離散化したスキームは、$\Delta t$ と $h$ をどれほど小さくしても、格子比 $\Delta t/h$ を固定するかぎり発散します（§9 の数値実験で、誤差が $10^{69}$ に達します）。差分商の 1 回ごとの誤差は小さくても、時間ステップを $T/\Delta t$ 回くり返すあいだに誤差が増幅されるからです。

したがって、スキームには 2 つの性質が要ります。1 回のステップで方程式を正しく近似していること（整合性）と、ステップをくり返しても誤差が増幅されないこと（安定性）です。この 2 つがあれば収束し、しかも整合的なスキームでは安定性は収束のために必要でもある、というのが Lax の同等定理です。線形方程式 $A\boldsymbol{x} = \boldsymbol{b}$ の数値解法での「前進誤差は条件数と後退誤差の積で抑えられる」という関係（行列のノルムと線形方程式の条件数）の、発展方程式版にあたります。打ち切り誤差が後退誤差、安定性の定数 $K$ が条件数の役です。

この章では、まず抽象的な枠組みで同等定理を証明し、次に安定性を調べる道具（von Neumann 解析）を作り、移流方程式の代表的な 4 つのスキームに適用します。この枠組みは、次章以降の拡散方程式・有限要素法・スペクトル法でも同じ形で使います。

2. 差分商と打ち切り誤差

定義 2.1（格子と格子関数）↗

格子幅 $h > 0$ に対し、$x_j = jh$（$j \in \mathbb{Z}$）を格子点、集合 $h\mathbb{Z} = \{x_j\}$ を格子という。格子上の複素数値の関数 $v = (v_j)_{j\in\mathbb{Z}}$ を格子関数という。

$$\|v\|_h = \Bigl(h\sum_{j\in\mathbb{Z}}|v_j|^2\Bigr)^{1/2},\qquad \|v\|_\infty = \sup_{j}|v_j|$$

を離散 $\ell^2$ ノルム、最大値ノルムといい、$\|v\|_h < \infty$ である格子関数の全体を $\ell^2(h\mathbb{Z})$ と書く。周期 $1$ の問題では $h = 1/N$（$N$ は正の整数）とし、$v_{j+N} = v_j$ を満たす格子関数を $\mathbb{C}^N$ の元と同一視して、和を $j = 0, \ldots, N-1$ でとる。時間の刻み $\Delta t > 0$ に対し $t_n = n\Delta t$ とし、時刻 $t_n$ の格子関数を $u^n = (u^n_j)_j$ と書く。

重み $h$ を付けるのは、$v_j = f(x_j)$ のとき $\|v\|_h^2$ が $\int|f|^2dx$ の Riemann 和になり、$h$ を変えてもノルムの大きさが比べられるようにするためです。

定義 2.2（前進・後退・中心差分）↗

関数 $f$ と $h > 0$ に対し

$$D_+f(x) = \frac{f(x+h) - f(x)}{h},\quad D_-f(x) = \frac{f(x) - f(x-h)}{h},\quad D_0f(x) = \frac{f(x+h) - f(x-h)}{2h}$$

を前進差分・後退差分・中心差分といい、$D_+D_-f(x) = \dfrac{f(x+h) - 2f(x) + f(x-h)}{h^2}$ を 2 階中心差分という。差分商と、それが近似する導関数との差を、その差分商の打ち切り誤差という。

補題 2.3（差分商の打ち切り誤差）✓

$f$ は $[x-h, x+h]$ を含む開区間で定義された実数値関数とし、$M_k = \sup_{[x-h,x+h]}|f^{(k)}|$ とおく。

1. $f$ が $C^2$ 級なら $|D_\pm f(x) - f'(x)| \le \dfrac{h}{2}M_2$。
2. $f$ が $C^3$ 級なら $|D_0f(x) - f'(x)| \le \dfrac{h^2}{6}M_3$。
3. $f$ が $C^4$ 級なら $|D_+D_-f(x) - f''(x)| \le \dfrac{h^2}{12}M_4$。

証明（補題 2.3）

定理 5.3［平均値の定理とテイラーの定理］ を点 $x$ のまわりで使います。

1. $n = 1$ として、$x$ と $x \pm h$ の間の $c_\pm$ があって $f(x \pm h) = f(x) \pm hf'(x) + \frac{h^2}{2}f''(c_\pm)$ です。したがって $D_\pm f(x) - f'(x) = \pm\frac{h}{2}f''(c_\pm)$ で、絶対値は $\frac{h}{2}M_2$ 以下です。

2. $n = 2$ として $f(x\pm h) = f(x) \pm hf'(x) + \frac{h^2}{2}f''(x) \pm \frac{h^3}{6}f'''(c_\pm)$ です。差をとると偶数次の項が消えて

$$f(x+h) - f(x-h) = 2hf'(x) + \frac{h^3}{6}\bigl(f'''(c_+) + f'''(c_-)\bigr)$$

となり、$2h$ で割って $|D_0f(x) - f'(x)| \le \frac{h^2}{12}\cdot 2M_3 = \frac{h^2}{6}M_3$ です。

3. $n = 3$ として $f(x\pm h) = f(x) \pm hf'(x) + \frac{h^2}{2}f''(x) \pm \frac{h^3}{6}f'''(x) + \frac{h^4}{24}f^{(4)}(c_\pm)$ です。和をとると奇数次の項が消えて

$$f(x+h) - 2f(x) + f(x-h) = h^2f''(x) + \frac{h^4}{24}\bigl(f^{(4)}(c_+) + f^{(4)}(c_-)\bigr)$$

となり、$h^2$ で割って $|D_+D_-f(x) - f''(x)| \le \frac{h^2}{24}\cdot 2M_4 = \frac{h^2}{12}M_4$ です。

∎

中心差分が 2 次なのは、対称性で偶数次（あるいは奇数次）の項が打ち消し合うからです。浮動小数点演算では $h$ を小さくしすぎると桁落ちで誤差が増えますが、その釣り合いは Richardson 補外の章 で扱っています。この章では丸め誤差は扱わず、打ち切り誤差とその積み上がり方だけを調べます。

3. 適切な初期値問題と差分スキーム

この節は、次の節の同等定理を述べるための準備で、「解が一意に存在し、初期値に連続に依存する」という当たり前の条件を、抽象的な言葉で言い直します。熱方程式や移流方程式のような具体的な方程式は、どれも次の定義の $E(t)$（初期値を時間 $t$ だけ先へ進める作用素。$E(t)u_0$ が時刻 $t$ の解）を与える例になっています。$X$ は解の置かれる関数の空間（Banach 空間）で、$\mathcal{B}(X)$ は $X$ から $X$ への有界な線形写像の全体です。

定義 3.1（適切な線形初期値問題）↗

$X$ を Banach 空間とする。$X$ 上の強連続半群 $E(t)$（$t \ge 0$。定義 4.1［1 径数ユニタリ群と強連続半群］）が与えられたとき、その生成子を $L$ として、初期値問題

$$\frac{du}{dt} = Lu\quad(0 \le t \le T),\qquad u(0) = u_0$$

は適切であるといい、$E(t)$ を解作用素、$u(t) = E(t)u_0$（$u_0 \in X$）を初期値 $u_0$ の解という。

命題 4.2［1 径数ユニタリ群と強連続半群］ の 1 から、定数 $M \ge 1$、$\omega \in \mathbb{R}$ があって $\|E(t)\| \le Me^{\omega t}$ なので、$0 \le t \le T$ で $\|E(t)\| \le K_E := Me^{|\omega|T}$ です。同じ命題の 2 から、$u_0 \in D(L)$ なら $u(t) = E(t)u_0$ は実際に微分可能で $u' = Lu$ を満たし、$u_0 \notin D(L)$ のときの $E(t)u_0$ は、そのような解の極限（一般化された解）です。「適切」とは、解が一意に存在して初期値に連続に依存する、という Hadamard の条件を半群の言葉で述べたものです。

定義 3.2（1 段の線形差分スキーム）↗

$\Delta t_0 > 0$ とする。各 $0 < \Delta t \le \Delta t_0$ に有界線形作用素 $C(\Delta t) \in \mathcal{B}(X)$ を対応させる族で、任意の $0 < \delta \le \Delta t_0$ について $\sup_{\delta \le \Delta t \le \Delta t_0}\|C(\Delta t)\| < \infty$ を満たすものを 1 段の線形差分スキームという。$u^0 = u_0$、$u^{n+1} = C(\Delta t)u^n$、すなわち $u^n = C(\Delta t)^nu_0$ を、時刻 $t_n = n\Delta t$ での近似解とする。空間の格子幅 $h$ は $\Delta t$ の関数 $h = h(\Delta t)$（$\Delta t \to 0$ で $h \to 0$）として与えられているものとする。

$h$ を $\Delta t$ の関数とするのは、極限を 1 つのパラメータでとるためです。双曲型の方程式では $\Delta t/h = \lambda$ を、放物型では $\Delta t/h^2 = \mu$ を一定にするのが典型です。陰的なスキーム $B_1u^{n+1} = B_0u^n$ は $C = B_1^{-1}B_0$ としてこの形に入ります。

定義 3.3（局所打ち切り誤差・整合性・精度の次数）↗

初期値 $u_0 \in X$ の解 $u(t) = E(t)u_0$ に対し

$$\tau(t, \Delta t) = \frac{u(t + \Delta t) - C(\Delta t)u(t)}{\Delta t}\qquad(0 \le t \le T - \Delta t)$$

をスキームの局所打ち切り誤差という。$X$ の稠密な部分空間 $D$ があって、すべての $u_0 \in D$ について

$$\sup_{0 \le t \le T - \Delta t}\|\tau(t, \Delta t)\| \to 0\qquad(\Delta t \to 0)$$

となるとき、スキームは整合的であるという。初期値のあるクラスについて、定数 $c$（$u_0$ による）があって $\sup_t\|\tau(t,\Delta t)\| \le c(\Delta t^p + h^q)$ となるとき、スキームはそのクラスについて時間 $p$ 次・空間 $q$ 次の精度をもつという。

整合性の定義で稠密な部分空間 $D$ の上だけを要求するのは、滑らかな初期値では打ち切り誤差を Taylor の定理で評価できるからです。滑らかでない初期値は、$D$ の元で近似して扱います（これは §4 の定理の証明の後半で使います）。

$\tau$ は「厳密解を 1 ステップ進めたもの」と「厳密解にスキームを 1 回適用したもの」の差を $\Delta t$ で割った量です。$\Delta t$ で割るのは、$u_t = Lu$ の残差と同じ単位にするためで、差分商の打ち切り誤差（補題 2.3）がそのまま現れます。

定義 3.4（安定性と収束）↗

スキーム $C(\Delta t)$ が安定であるとは、定数 $K$ があって

$$\|C(\Delta t)^n\| \le K\qquad(0 < \Delta t \le \Delta t_0,\ n \ge 0,\ n\Delta t \le T)$$

が成り立つことをいう。スキームが収束するとは、任意の $u_0 \in X$、任意の $t \in [0, T]$、および $\Delta t_j \to 0$、$n_j\Delta t_j \le T$、$n_j\Delta t_j \to t$ を満たす任意の列 $(\Delta t_j, n_j)$ について

$$\|C(\Delta t_j)^{n_j}u_0 - E(t)u_0\| \to 0\qquad(j \to \infty)$$

となることをいう。

安定性は、スキームだけの性質で、方程式にも解にもよりません。また、$\|C(\Delta t)\| \le 1$ を要求しているのではなく、$n \le T/\Delta t$ 回の積が $\Delta t$ によらない定数で抑えられることを要求しています。たとえば $\|C(\Delta t)\| \le 1 + c\Delta t$ なら $\|C^n\| \le (1 + c\Delta t)^n \le e^{cn\Delta t} \le e^{cT}$ で安定です。

4. Lax の同等定理

誤差の積み上がり方を、空間によらない形の補題にしておきます。$X$ 上のスキームにも、格子関数の空間 $\ell^2(h\mathbb{Z})$ や $\mathbb{C}^N$ の上のスキームにもそのまま使えます。

補題 4.1（誤差の漸化式）✓

$Y$ をノルム空間、$C \in \mathcal{B}(Y)$、$\Delta t > 0$、$N \ge 1$ とし、$\|C^n\| \le K$（$0 \le n \le N$）とする。$Y$ の元の列 $w^0, w^1, \ldots, w^N$ に対し $\tau^n = (w^{n+1} - Cw^n)/\Delta t$（$0 \le n < N$）とおく。このとき、任意の $u^0 \in Y$ と $u^n = C^nu^0$ について、$e^n = w^n - u^n$ は

$$e^{n+1} = Ce^n + \Delta t\,\tau^n,\qquad \|e^n\| \le K\|e^0\| + K\,n\Delta t\max_{0 \le k < n}\|\tau^k\|\qquad(0 \le n \le N)$$

を満たす。

証明（補題 4.1）

$\tau^n$ の定義から $w^{n+1} = Cw^n + \Delta t\,\tau^n$ で、$u^{n+1} = Cu^n$ を引くと $e^{n+1} = Ce^n + \Delta t\,\tau^n$ です。これを $n$ についてくり返すと

$$e^n = C^ne^0 + \Delta t\sum_{k=0}^{n-1}C^{n-1-k}\tau^k$$

です（$n = 0$ では和が空で成り立ち、$n$ で成り立てば $e^{n+1} = C^{n+1}e^0 + \Delta t\sum_{k=0}^{n-1}C^{n-k}\tau^k + \Delta t\,\tau^n$ で $n+1$ でも成り立ちます）。$0 \le n - 1 - k \le N$ なので $\|C^{n-1-k}\| \le K$ で、三角不等式から $\|e^n\| \le K\|e^0\| + \Delta t\cdot n\cdot K\max_k\|\tau^k\|$ です。

∎

$w^n$ として厳密解の時刻 $t_n$ での値（格子上のスキームなら格子点での値）をとると、$\tau^n$ は局所打ち切り誤差です。この補題は「大域誤差は、初期誤差と、$T/\Delta t$ 回ぶんの局所誤差 $\Delta t\,\tau$ の和の $K$ 倍以下」と読めます。1 ステップの誤差 $\Delta t\,\tau$ が $O(\Delta t^{p+1})$ でも、積み上げると $O(\Delta t^p)$ になります。

次の定理の証明は 2 段に分かれます。まず滑らかな初期値では、上の補題の右辺がそのまま誤差の上界になります。次に一般の初期値は、滑らかな $v$ で近似して、誤差を 4 つの項に分けます。$u_0$ と $v$ の差に由来する 2 項は、安定性と解作用素の有界性のおかげで $\varepsilon$ に比例して小さくなります。$v$ についてのスキームの誤差は滑らかな場合の結果で小さくなり、残る 1 項は解の時間連続性で小さくなります。

定理 4.2（整合性と安定性から収束が従う）✓

適切な線形初期値問題（定義 3.1）に対するスキーム $C(\Delta t)$ が整合的で安定ならば、収束する。さらに、初期値 $u_0$ について $\sup_t\|\tau(t,\Delta t)\| \le c(\Delta t^p + h^q)$ ならば、安定性の定数 $K$ について

$$\|u(t_n) - C(\Delta t)^nu_0\| \le KTc\,(\Delta t^p + h^q)\qquad(n\Delta t \le T)$$

が成り立つ。

証明（定理 4.2）

滑らかな初期値. $u_0 \in X$ を任意にとり、$u(t) = E(t)u_0$、$w^n = u(t_n)$、$N = \lfloor T/\Delta t\rfloor$ とします。$\tau^n = (w^{n+1} - C(\Delta t)w^n)/\Delta t = \tau(t_n, \Delta t)$ で、安定性から $\|C(\Delta t)^n\| \le K$（$n \le N$）です。補題 4.1 を $u^0 = u_0$（$e^0 = 0$）で使うと

$$\|u(t_n) - C(\Delta t)^nu_0\| \le K\,n\Delta t\sup_{0\le t\le T-\Delta t}\|\tau(t,\Delta t)\| \le KT\sup_t\|\tau(t,\Delta t)\|$$

です。$\sup_t\|\tau\| \le c(\Delta t^p + h^q)$ なら後半の評価が得られます。また $u_0 \in D$（整合性の定義の稠密な部分空間）なら、右辺は $\Delta t \to 0$ で $0$ に収束します。

一般の初期値. $u_0 \in X$、$t \in [0,T]$、列 $(\Delta t_j, n_j)$ を収束の定義のとおりにとり、$s_j = n_j\Delta t_j$ とします。$\varepsilon > 0$ に対し、$D$ の稠密性から $\|u_0 - v\| < \varepsilon$ となる $v \in D$ をとります。

$$C(\Delta t_j)^{n_j}u_0 - E(t)u_0 = C(\Delta t_j)^{n_j}(u_0 - v) + \bigl(C(\Delta t_j)^{n_j}v - E(s_j)v\bigr) + \bigl(E(s_j)v - E(t)v\bigr) + E(t)(v - u_0)$$

と分けます。第 1 項のノルムは安定性から $K\varepsilon$ 以下、第 4 項は $K_E\varepsilon$ 以下です。第 2 項は前半から $KT\sup_t\|\tau_v(t,\Delta t_j)\|$ 以下（$\tau_v$ は $v$ の解についての局所打ち切り誤差）で、$v \in D$ なので $j \to \infty$ で $0$ に収束します。第 3 項は、$t \mapsto E(t)v$ の連続性（強連続半群の定義）と $s_j \to t$ から $0$ に収束します。したがって

$$\limsup_{j\to\infty}\|C(\Delta t_j)^{n_j}u_0 - E(t)u_0\| \le (K + K_E)\varepsilon$$

で、$\varepsilon$ は任意なので左辺は $0$ です。

∎

一般の初期値では収束はしますが、速さは保証されません。次数つきの評価には解の滑らかさが要ります（§9 で、角のある初期値では観測される次数が落ちることを見ます）。

逆向きの定理は、対偶で示します。つまり、安定でないスキームには、作用素のノルムの増え方が速すぎて、ある初期値で解が暴れてしまう例が必ず作れる、という流れです。「ノルムは大きいが各初期値では大丈夫」という状況を許さないのが、一様有界性原理です。

定理 4.3（収束するスキームは安定である）✓

適切な線形初期値問題に対するスキーム $C(\Delta t)$ が収束するならば、安定である。

証明（定理 4.3）

安定でないと仮定します。すると各 $j \ge 1$ について、$0 < \Delta t_j \le \Delta t_0$、$n_j\Delta t_j \le T$ で $\|C(\Delta t_j)^{n_j}\| > j$ となる組 $(\Delta t_j, n_j)$ があります。$T_j = C(\Delta t_j)^{n_j} \in \mathcal{B}(X)$ とおくと $\sup_j\|T_j\| = \infty$ です。$X$ は Banach 空間なので、定理 5.1［Banach 空間の基本定理］ の対偶から、$\sup_j\|T_ju_0\| = \infty$ となる $u_0 \in X$ があります。部分列をとって $\|T_ju_0\| \to \infty$ としてよく、さらに有界な数列 $\Delta t_j \in (0, \Delta t_0]$、$n_j\Delta t_j \in [0, T]$ から収束する部分列をとって、$\Delta t_j \to \delta \in [0, \Delta t_0]$、$n_j\Delta t_j \to t \in [0, T]$ としてよいです。

$\delta = 0$ の場合. 収束の定義から $T_ju_0 = C(\Delta t_j)^{n_j}u_0 \to E(t)u_0$ なので、$\|T_ju_0\|$ は有界です。これは $\|T_ju_0\| \to \infty$ に反します。

$\delta > 0$ の場合. 十分大きい $j$ で $\Delta t_j \ge \delta/2$ なので、$n_j \le T/\Delta t_j \le 2T/\delta$ です。スキームの定義（定義 3.2）から $m := \sup_{\delta/2 \le \Delta t \le \Delta t_0}\|C(\Delta t)\| < \infty$ で、$\|T_ju_0\| \le \|C(\Delta t_j)\|^{n_j}\|u_0\| \le \max(1, m)^{2T/\delta}\|u_0\|$ です。これも $\|T_ju_0\| \to \infty$ に反します。

どちらの場合も矛盾するので、スキームは安定です。

∎

この証明は整合性を使っていません。収束の定義が「すべての初期値について」であることが効いています。滑らかな初期値についてだけ収束する不安定なスキームは存在しますが、丸め誤差は滑らかではないので、実際の計算では役に立ちません。

定理 4.4（Lax の同等定理）✓

適切な線形初期値問題に対する整合的な 1 段の線形差分スキームについて、安定であることと収束することは同値である。

証明（定理 4.4）

安定なら収束することは 定理 4.2、収束するなら安定であることは 定理 4.3 です。

∎

注意 4.5（格子の上で使うとき）

実際のスキームは格子関数に作用するので、空間が $h$ ごとに変わります。そのときは、厳密解の格子点での値 $w^n_j = u(x_j, t_n)$ をとって 補題 4.1 を空間 $\ell^2(h\mathbb{Z})$ や $\mathbb{C}^N$ で使えば、定理 4.2 の後半と同じ評価 $\|w^n - u^n\| \le KT\max_k\|\tau^k\|$ が得られます。以下の具体的なスキームでは、この形で使います。同等定理の「収束 $\Rightarrow$ 安定」を格子の上で述べるには、格子関数を $\mathbb{R}$ 上の関数に延長する作用素を決める必要があります。Richtmyer–Morton（参考文献）がこの扱いを詳しく述べています。

5. von Neumann の安定性解析

安定性の定義は冪 $C^n$ のノルムについての条件で、そのままでは確かめにくいものです。係数が格子点によらないスキームでは、Fourier 変換が $C$ を掛け算に変えるので、ノルムが計算できます。

補題 5.1（半離散 Fourier 変換）✓

$v \in \ell^2(h\mathbb{Z})$ に対し

$$\hat v(\xi) = h\sum_{j\in\mathbb{Z}}v_je^{-ix_j\xi}\qquad\Bigl(-\frac{\pi}{h} \le \xi \le \frac{\pi}{h}\Bigr)$$

とおく（和は $L^2(-\pi/h, \pi/h)$ で収束する）。このとき次が成り立つ。

1. $\displaystyle\int_{-\pi/h}^{\pi/h}|\hat v(\xi)|^2d\xi = 2\pi\|v\|_h^2$ で、$v \mapsto (2\pi)^{-1/2}\hat v$ は $\ell^2(h\mathbb{Z})$ から $L^2(-\pi/h, \pi/h)$ の上へのユニタリ写像である。
2. 平行移動 $(Sv)_j = v_{j+1}$ について $\widehat{Sv}(\xi) = e^{ih\xi}\hat v(\xi)$。

証明（補題 5.1）

定理 5.1［フーリエ級数］ から、$F \mapsto (\hat F(n))_{n\in\mathbb{Z}}$ は $L^2(\mathbb{T})$（内積は $\frac{1}{2\pi}\int_{-\pi}^{\pi}\bar FG\,d\theta$）から $\ell^2(\mathbb{Z})$ の上へのユニタリ写像で、$S_NF \to F$ が $L^2$ で成り立ちます。Fourier 係数は 定義 2.1［フーリエ級数］ のものです。

1. $v \in \ell^2(h\mathbb{Z})$ に対し、上への写像であることから $\hat F(n) = v_{-n}$（$n \in \mathbb{Z}$）となる $F \in L^2(\mathbb{T})$ がただ 1 つあり、$F = \lim_N\sum_{|n|\le N}v_{-n}e^{in\theta} = \sum_jv_je^{-ij\theta}$（$L^2$ 収束）、$\frac{1}{2\pi}\int_{-\pi}^{\pi}|F|^2d\theta = \sum_j|v_j|^2$ です。$x_j\xi = j\cdot h\xi$ なので $\hat v(\xi) = hF(h\xi)$ で、変数変換 $\theta = h\xi$ により

$$\int_{-\pi/h}^{\pi/h}|\hat v(\xi)|^2d\xi = h^2\int_{-\pi}^{\pi}|F(\theta)|^2\frac{d\theta}{h} = 2\pi h\sum_j|v_j|^2 = 2\pi\|v\|_h^2$$

です。上への写像であることを示します。$G \in L^2(-\pi/h, \pi/h)$ に対し $F(\theta) = G(\theta/h)/h$ は $L^2(\mathbb{T})$ に属し、$v_j = \hat F(-j)$ とおけば $v \in \ell^2(h\mathbb{Z})$ で $\hat v = hF(h\,\cdot) = G$ です。ノルムを保つ上への線形写像なのでユニタリです。

2. $Sv$ に対応する関数を $F_S$ とすると $\hat F_S(n) = (Sv)_{-n} = v_{-n+1} = \hat F(n-1)$ です。一方 $e^{i\theta}F$ の Fourier 係数は $\frac{1}{2\pi}\int F(\theta)e^{i\theta}e^{-in\theta}d\theta = \hat F(n-1)$ です。Fourier 係数が一致するので $F_S = e^{i\theta}F$ で、$\widehat{Sv}(\xi) = hF_S(h\xi) = e^{ih\xi}\hat v(\xi)$ です。

∎

この補題の中身は、格子上の数列を周期 $2\pi$ の関数の Fourier 係数とみなせば、Fourier 級数の Parseval の等式と同じだ、ということです。変数変換が現れるのは、格子幅 $h$ に合わせて $\theta = h\xi$ と目盛りを取り直すためです。したがって、$v \mapsto \hat v$ は長さを変えない対応（ユニタリ写像）で、スキームの冪のノルムを波数ごとの掛け算で計算できます。

波数 $\xi$ が $[-\pi/h, \pi/h]$ に限られるのは、格子の上では $e^{ix_j\xi}$ と $e^{ix_j(\xi + 2\pi/h)}$ が区別できないからです。$\theta = h\xi \in [-\pi, \pi]$ は 1 格子あたりの位相の進みで、$\theta = \pm\pi$ が格子で表せる最も細かい振動 $v_j = (-1)^j$ にあたります。

定義 5.2（定数係数のスキームと増幅因子）↗

整数 $l, r \ge 0$ と、格子点 $j$ によらない複素数 $a_k = a_k(\Delta t, h)$（$-l \le k \le r$）で

$$(Cv)_j = \sum_{k=-l}^{r}a_kv_{j+k}$$

と書けるスキームを定数係数の陽的スキームといい、$g(\theta) = \sum_{k=-l}^{r}a_ke^{ik\theta}$（$\theta \in [-\pi, \pi]$）をその増幅因子（シンボル）という。$g$ も $\Delta t$ と $h$ によるので、必要なときは $g(\theta, \Delta t)$ と書く。

$v_j = e^{ij\theta}$ を代入すると $(Cv)_j = \sum_ka_ke^{i(j+k)\theta} = g(\theta)v_j$ なので、$g(\theta)$ は波数 $\theta/h$ の波が 1 ステップで何倍になるかを表します。補題 5.1 の 2 から $C = \sum_ka_kS^k$ について $\widehat{Cv}(\xi) = g(h\xi)\hat v(\xi)$ です。

補題 5.3（掛け算作用素のノルム）✓

定数係数の陽的スキーム $C$ は $\ell^2(h\mathbb{Z})$ 上の有界作用素で、すべての $n \ge 0$ について $\|C^n\| = \bigl(\max_{\theta\in[-\pi,\pi]}|g(\theta)|\bigr)^n$ である。周期 $1$ の格子（$h = 1/N$、空間 $\mathbb{C}^N$、添字は $N$ を法とする）では、$\theta_m = 2\pi m/N$ として $\|C^n\| = \bigl(\max_{0\le m<N}|g(\theta_m)|\bigr)^n \le \bigl(\max_\theta|g(\theta)|\bigr)^n$ である。

証明（補題 5.3）

$G = \max_\theta|g(\theta)|$ とおきます（$g$ は連続で $[-\pi,\pi]$ はコンパクトなので最大値があります）。

$\ell^2(h\mathbb{Z})$ の場合. $\widehat{C^nv} = g(h\xi)^n\hat v$ なので、補題 5.1 の 1 から

$$2\pi\|C^nv\|_h^2 = \int_{-\pi/h}^{\pi/h}|g(h\xi)|^{2n}|\hat v(\xi)|^2d\xi \le G^{2n}\cdot 2\pi\|v\|_h^2$$

で、$\|C^n\| \le G^n$ です。逆向きを示します。$|g(\theta_0)| = G$ となる $\theta_0$ をとり、$\varepsilon > 0$ に対し、$g$ の連続性から $|\theta - \theta_0| < \delta$ なら $|g(\theta)| > G - \varepsilon$ となる $\delta > 0$ をとります。$I = \{\xi \in [-\pi/h, \pi/h] : |h\xi - \theta_0| < \delta\}$ は長さが正の区間で、その定義関数 $1_I$ は $L^2(-\pi/h,\pi/h)$ に属するので、補題 5.1 の 1 から $\hat v = 1_I$ となる $v \in \ell^2(h\mathbb{Z})$ があります。このとき

$$2\pi\|C^nv\|_h^2 = \int_I|g(h\xi)|^{2n}d\xi \ge (G - \varepsilon)^{2n}\int_I d\xi = (G-\varepsilon)^{2n}\cdot2\pi\|v\|_h^2$$

（$G - \varepsilon < 0$ のときは自明）なので $\|C^n\| \ge (G - \varepsilon)^n$ で、$\varepsilon \to 0$ として $\|C^n\| \ge G^n$ です。

周期格子の場合. $\varphi^{(m)}_j = e^{ij\theta_m}$（$0 \le m < N$）とおくと、$m \ne m'$ のとき $\omega = e^{i(\theta_m - \theta_{m'})} \ne 1$、$\omega^N = 1$ なので

$$h\sum_{j=0}^{N-1}\overline{\varphi^{(m')}_j}\varphi^{(m)}_j = h\sum_{j=0}^{N-1}\omega^j = h\,\frac{\omega^N - 1}{\omega - 1} = 0$$

で、$m = m'$ のときは $hN = 1$ です。したがって $\{\varphi^{(m)}\}$ は内積 $h\sum_j\bar v_jw_j$ についての $\mathbb{C}^N$ の正規直交基底です。$e^{i(j+k)\theta_m}$ は $j + k$ について周期 $N$ なので $C\varphi^{(m)} = g(\theta_m)\varphi^{(m)}$ です。$v = \sum_mc_m\varphi^{(m)}$ と展開すると $\|C^nv\|_h^2 = \sum_m|g(\theta_m)|^{2n}|c_m|^2 \le \max_m|g(\theta_m)|^{2n}\|v\|_h^2$ で、最大を与える $m$ の $\varphi^{(m)}$ で等号になります。

∎

補題 5.4（冪の有界性から 1 ステップの評価へ）✓

$K \ge 1$、$0 < \Delta t_0 \le T$ とする。$0 < \Delta t \le \Delta t_0$ と $\rho \ge 0$ が、$n\Delta t \le T$ を満たすすべての整数 $n \ge 0$ について $\rho^n \le K$ を満たすならば、$\rho \le 1 + C_0\Delta t$ である。ここで $C_0 = (K^{2\Delta t_0/T} - 1)/\Delta t_0$ は $\Delta t$ と $\rho$ によらない。

証明（補題 5.4）

$n = \lfloor T/\Delta t\rfloor$ とすると、$T/\Delta t \ge 1$ なので $n \ge 1$ で、$y \ge 1$ のとき $\lfloor y\rfloor \ge y/2$ であること（$1 \le y < 2$ では $\lfloor y\rfloor = 1 \ge y/2$、$y \ge 2$ では $\lfloor y\rfloor \ge y - 1 \ge y/2$）から $n \ge T/(2\Delta t)$ です。$\rho^n \le K$ と $K \ge 1$ から $\rho \le K^{1/n} \le K^{2\Delta t/T} = e^{c\Delta t}$（$c = 2\ln K/T \ge 0$）です。関数 $s \mapsto e^{cs}$ は凸なので、区間 $[0, \Delta t_0]$ では両端を結ぶ弦の下にあり、$e^{c\Delta t} \le 1 + \frac{e^{c\Delta t_0} - 1}{\Delta t_0}\Delta t = 1 + C_0\Delta t$ です。

∎

定理 5.5（von Neumann の安定性条件）✓

定数係数の陽的スキーム $C(\Delta t)$（$0 < \Delta t \le \Delta t_0 \le T$、$h = h(\Delta t)$）について、次は同値である。

1. スキームは $\ell^2(h\mathbb{Z})$ で安定である。
2. 定数 $C_0 \ge 0$ があって、すべての $0 < \Delta t \le \Delta t_0$ で $\max_{\theta\in[-\pi,\pi]}|g(\theta, \Delta t)| \le 1 + C_0\Delta t$。

特に、$g$ が $\Delta t$ によらないとき（係数 $a_k$ が格子比だけで決まるとき）、安定であることと $\max_\theta|g(\theta)| \le 1$ は同値である。周期格子の上では、2 が成り立てば安定である。

証明（定理 5.5）

$G(\Delta t) = \max_\theta|g(\theta, \Delta t)|$ とおくと、補題 5.3 から $\|C(\Delta t)^n\| = G(\Delta t)^n$ です。

2 $\Rightarrow$ 1. $n\Delta t \le T$ のとき $\|C(\Delta t)^n\| \le (1 + C_0\Delta t)^n \le e^{C_0n\Delta t} \le e^{C_0T}$ です（$1 + s \le e^s$ を使いました）。周期格子でも 補題 5.3 の不等式から同じ評価が成り立ちます。

1 $\Rightarrow$ 2. 安定性の定数を $K$ とします（$n = 0$ の場合から $K \ge 1$）。$n\Delta t \le T$ のとき $G(\Delta t)^n = \|C(\Delta t)^n\| \le K$ なので、補題 5.4 から $G(\Delta t) \le 1 + C_0\Delta t$ です。

$g$ が $\Delta t$ によらないとき. $\max|g| \le 1$ なら 2 が $C_0 = 0$ で成り立ちます。逆に 2 が成り立てば、$\Delta t \to 0$ として $\max|g| \le 1$ です。

∎

$1 + C_0\Delta t$ の余裕は、$u_t = u$ のように解そのものが指数的に増える方程式を扱うために必要です。解が増えない方程式（移流・拡散）では、$|g| \le 1$ を条件にするのが普通です。

連立の方程式では、$a_k$ が $m$ 次の行列 $A_k$ になり、増幅因子は増幅行列 $G(\theta) = \sum_kA_ke^{ik\theta}$ になります。格子関数は $\mathbb{C}^m$ に値をとり、ノルムは $\|v\|_h^2 = h\sum_j|v_j|^2$（$|\cdot|$ は Euclid ノルム）、Fourier 変換は成分ごとにとります。補題 5.1 を成分ごとに使えば、$\widehat{Cv}(\xi) = G(h\xi)\hat v(\xi)$ と $\int|\hat v|^2d\xi = 2\pi\|v\|_h^2$ が成り立ちます。

補題 5.6（増幅行列による下からの評価）✓

$(Cv)_j = \sum_{k=-l}^{r}A_kv_{j+k}$（$A_k \in M_m(\mathbb{C})$）とし、$G(\theta) = \sum_kA_ke^{ik\theta}$ とおく。任意の $\theta_0 \in [-\pi,\pi]$、$|w| = 1$ である $w \in \mathbb{C}^m$、$n \ge 0$ について、$\ell^2(h\mathbb{Z};\mathbb{C}^m)$ 上の作用素ノルムは $\|C^n\| \ge |G(\theta_0)^nw|$ を満たす。

証明（補題 5.6）

$\theta \mapsto G(\theta)^nw$ は連続なので、$\varepsilon > 0$ に対し $\delta > 0$ があって、$|\theta - \theta_0| < \delta$ なら $|G(\theta)^nw| > |G(\theta_0)^nw| - \varepsilon$ です。補題 5.3 の証明と同じ区間 $I$ をとり、$\hat v = 1_Iw$ となる $v$ をとります（成分ごとに 補題 5.1 を使う）。$\widehat{C^nv}(\xi) = G(h\xi)^nw\,1_I(\xi)$ なので

$$2\pi\|C^nv\|_h^2 = \int_I|G(h\xi)^nw|^2d\xi \ge \bigl(|G(\theta_0)^nw| - \varepsilon\bigr)^2\int_Id\xi = \bigl(|G(\theta_0)^nw| - \varepsilon\bigr)^2\cdot2\pi\|v\|_h^2$$

（括弧の中が負のときは自明）です。$\varepsilon \to 0$ として結論を得ます。

∎

系 5.7（連立のスキームの必要条件）✓

補題 5.6 の連立のスキームが $\ell^2$ で安定（$0 < \Delta t \le \Delta t_0 \le T$）ならば、定数 $C_0$ があって、増幅行列のスペクトル半径（定義 5.1［行列のノルムと線形方程式の条件数］）は $\rho(G(\theta, \Delta t)) \le 1 + C_0\Delta t$（すべての $\theta$ と $\Delta t$）を満たす。

証明（系 5.7）

$\theta_0$ を固定し、$G(\theta_0)$ の固有値 $\mu$ で $|\mu| = \rho(G(\theta_0))$ となるものと、その固有ベクトル $w$（$|w| = 1$）をとります。$G(\theta_0)^nw = \mu^nw$ なので、補題 5.6 と安定性から $\rho(G(\theta_0))^n = |\mu|^n \le \|C^n\| \le K$（$n\Delta t \le T$）です。補題 5.4 から $\rho(G(\theta_0)) \le 1 + C_0\Delta t$ で、$C_0$ は $\theta_0$ によりません。

∎

スカラーの場合と違い、この条件は十分ではありません。定理 5.2［行列のノルムと線形方程式の条件数］ のとおり $\rho(G) \le \|G\|$ ですが、等号は一般には成り立たず、$\rho(G) = 1$ でも $\|G^n\|$ が $n$ とともに増えることがあるからです。

例 5.8（スペクトル半径が 1 でも不安定な連立スキーム）

方程式 $u_t = 0$、$v_t = 0$（解は初期値のまま動かない）に対するスキーム

$$u^{n+1}_j = u^n_j + (v^n_{j+1} - 2v^n_j + v^n_{j-1}),\qquad v^{n+1}_j = v^n_j$$

を、$\Delta t/h = \lambda$ を固定して考えます。滑らかな解 $(u(x), v(x))$ を代入した残差は、第 1 式で $-(v(x+h) - 2v(x) + v(x-h))/\Delta t$ です。Taylor の定理を 2 階の剰余まで使うと $v(x\pm h) = v(x) \pm hv'(x) + \frac{h^2}{2}v''(c_\pm)$ なので、残差の絶対値は $h^2\sup|v''|/\Delta t = (h/\lambda)\sup|v''|$ 以下です。したがって整合的です。増幅行列は

$$G(\theta) = \begin{pmatrix}1 & 2\cos\theta - 2\\ 0 & 1\end{pmatrix},\qquad G(\theta)^n = \begin{pmatrix}1 & -4n\sin^2(\theta/2)\\ 0 & 1\end{pmatrix}$$

で、固有値は $1$ だけなので $\rho(G(\theta)) = 1$ です。ところが $\theta_0 = \pi$、$w = (0, 1)^{\mathsf{T}}$ とすると $|G(\pi)^nw| = \sqrt{16n^2 + 1} \ge 4n$ で、補題 5.6 から $\|C^n\| \ge 4n$ です。$n = \lfloor T/\Delta t\rfloor$ は $\Delta t \to 0$ で限りなく大きくなるので、このスキームは不安定です。$G(\theta)$ が正規行列でない（対角化できない）ことが原因です。

注意 5.9（連立のスキームの安定性の判定）

増幅行列の族 $\{G(\theta, \Delta t)\}$ の冪が一様に有界であるための必要十分条件は Kreiss の行列定理が与えます（レゾルベントの評価 $\|(zI - G)^{-1}\| \le C/(|z| - 1)$ などの同値な条件）。境界のある領域でのスキームの安定性は Gustafsson–Kreiss–Sundström の理論（GKS 理論）が扱います。どちらも証明が長いので、この章では扱いません。Gustafsson–Kreiss–Oliger（参考文献）に証明があります。$G$ が正規行列（Hermite 行列やユニタリ行列など）なら $\|G^n\|_2 = \rho(G)^n$ なので、スカラーの場合と同じ判定ができます。

6. 依存領域と CFL 条件

ここからは移流方程式 $u_t + au_x = 0$（$a \in \mathbb{R}$ は定数）を扱います。

補題 6.1（移流方程式の解）✓

$u_0 \in C^1(\mathbb{R})$ とする。$u(x,t) = u_0(x - at)$ は、$u_t + au_x = 0$（$x \in \mathbb{R}$、$t \ge 0$）、$u(x, 0) = u_0(x)$ を満たすただ 1 つの $C^1$ 級の関数である。$u_0 \in C^m$ なら $\partial_t^i\partial_x^ju = (-a)^iu_0^{(i+j)}(x - at)$（$i + j \le m$）である。

証明（補題 6.1）

連鎖律から $\partial_t^i\partial_x^j[u_0(x - at)] = (-a)^iu_0^{(i+j)}(x-at)$ で、特に $u_t + au_x = -au_0' + au_0' = 0$、$u(x,0) = u_0(x)$ です。一意性を示します。$u$ を $C^1$ 級の解とし、$x_0 \in \mathbb{R}$ を固定して $\phi(t) = u(x_0 + at, t)$ とおくと、$\phi'(t) = au_x + u_t = 0$ なので $\phi(t) = \phi(0) = u_0(x_0)$ です。$x = x_0 + at$ と書き直せば $u(x,t) = u_0(x - at)$ です。

∎

直線 $x = x_0 + at$ が特性曲線で、解はその上で一定です。一般の準線形 1 階方程式での同じ構成が 定理 6.1［安定性・力学系・特性曲線法］ です。

定義 6.2（依存領域と CFL 数）↗

移流方程式の解の点 $(x,t)$ での値は、初期値の 1 点 $x - at$ での値だけで決まる。この点（の集合 $\{x - at\}$）を $(x,t)$ の依存領域という。定数係数の陽的スキーム $u^{n+1}_j = \sum_{k=-l}^{r}a_ku^n_{j+k}$ では、$u^n_j$ は $u^0_{j-nl}, \ldots, u^0_{j+nr}$ だけで決まる。区間 $[x_j - nlh,\ x_j + nrh]$ を $(x_j, t_n)$ の数値的な依存領域という。$\lambda = \Delta t/h$ を格子比、$\nu = a\Delta t/h = a\lambda$ を CFL 数（Courant 数）という。

$\lambda$ を固定すると $nh = t_n/\lambda$ なので、数値的な依存領域は $[x - lt/\lambda,\ x + rt/\lambda]$ で、格子を細かくしても変わりません。

3 点スキーム（l = r = 1）の数値的な依存領域（三角形の底辺）と特性曲線。左は |ν| ≤ 1 で、特性曲線の足が底辺の中にある。右は |ν| > 1 で、足が底辺の外に出る。右の場合、足の近くの初期値を変えても数値解は変わらないので、収束しえない。

定理 6.3（Courant–Friedrichs–Lewy の条件）✓

$a \in \mathbb{R}$、$\lambda > 0$ とし、移流方程式 $u_t + au_x = 0$ に対する定数係数の陽的スキーム $u^{n+1}_j = \sum_{k=-l}^{r}a_ku^n_{j+k}$ を、$\Delta t = 1/n$、$h = \Delta t/\lambda$、$u^0_j = u_0(x_j)$ で考える（係数 $a_k$ は任意で、$n$ によってもよい）。$-a\lambda \notin [-l, r]$、すなわち点 $(0, 1)$ の依存領域が数値的な依存領域に含まれないならば、$C^4$ 級で台がコンパクトな $u_0$ と定数 $c > 0$ があって、すべての $n \ge 1$ について

$$u^n_0 = 0,\qquad u(0, 1) \ge c,\qquad \Bigl(h\sum_j|u(x_j, 1) - u^n_j|^2\Bigr)^{1/2} \ge c$$

が成り立つ。したがって、スキームが収束するには $-l \le -a\lambda \le r$ が必要である。3 点スキーム（$l = r = 1$）では、これは $|\nu| = |a|\Delta t/h \le 1$ である。

証明（定理 6.3）

$J = [-l/\lambda,\ r/\lambda]$ とおきます。仮定から $-a \notin J$ で、$J$ は閉区間なので $d = \mathrm{dist}(-a, J) > 0$ です。$\varphi(s) = (1 - s^2)^5$（$|s| < 1$）、$\varphi(s) = 0$（$|s| \ge 1$）とおくと、$\varphi$ は $C^4$ 級（$s = \pm1$ で 4 階までの導関数が $0$ につながる）で、$|s| \le 1/2$ で $\varphi(s) \ge (3/4)^5$ です。$u_0(x) = \varphi\bigl(3(x + a)/d\bigr)$ とおくと、$u_0$ の台は $[-a - d/3,\ -a + d/3]$ で、$|x + a| \le d/6$ で $u_0(x) \ge (3/4)^5$ です。

数値解. $|x_j| \le d/6$ である $j$ を考えます。帰納法で、$u^n_j$ は $u^0_{j+k}$（$-nl \le k \le nr$）の一次結合です（$n = 0$ は自明で、$u^{n+1}_j = \sum_ka_ku^n_{j+k}$ から添字の範囲が左に $l$、右に $r$ 広がります）。$nh = 1/\lambda$ なので、これらの格子点 $x_{j+k}$ は $[x_j - l/\lambda,\ x_j + r/\lambda]$ に入り、この区間の点は $J$ からの距離が $d/6$ 以下です。$-a$ と $J$ の距離は $d$ なので、これらの点と $-a$ の距離は $5d/6 > d/3$ 以上で、$u_0$ の台の外です。したがって $u^0_{j+k} = 0$ で、$u^n_j = 0$ です。特に $u^n_0 = 0$ です。

厳密解. 補題 6.1 から $u(x, 1) = u_0(x - a)$ で、$|x| \le d/6$ なら $u(x,1) \ge (3/4)^5$ です。

$\ell^2$ での誤差. $|x_j| \le d/6$ である $j$ は $2\lfloor d/(6h)\rfloor + 1 \ge d/(6h)$ 個あるので（$y \ge 0$ のとき $2\lfloor y\rfloor + 1 \ge y$）、

$$h\sum_j|u(x_j,1) - u^n_j|^2 \ge h\cdot\frac{d}{6h}\cdot\Bigl(\frac34\Bigr)^{10} = \frac{d}{6}\Bigl(\frac34\Bigr)^{10}$$

です。$c = \min\bigl((3/4)^5,\ (3/4)^5\sqrt{d/6}\bigr)$ とすれば主張が成り立ちます。$\Delta t = 1/n \to 0$ でも誤差が $0$ に近づかないので、収束しません。

∎

CFL 条件は、スキームの係数を一切使わずに、情報の伝わる範囲だけから出る必要条件です。波動方程式 $u_{tt} = c^2u_{xx}$ では依存領域が区間 $[x - ct, x + ct]$ になり（定理 5.3［波動方程式］ の有限伝播速度）、同じ議論で $c\Delta t/h \le 1$ が必要になります。一方、熱方程式は依存領域が全空間なので、陽的スキームでは $\Delta t/h \to 0$ が必要です（次章の $\Delta t \le h^2/2$ はこれを満たします）。陰的スキームは、1 ステップで全格子点の値が結びつくので、数値的な依存領域が全空間になり、CFL 条件による制限を受けません。

CFL 条件は十分条件ではありません。次節で見るとおり、前進時間・中心空間のスキームは $|\nu| \le 1$ でも不安定です。

7. 移流方程式の 4 つのスキーム

$u_t + au_x = 0$ の代表的な 3 点スキームを、CFL 数 $\nu = a\Delta t/h$ を使って書きます。

名前	$u^{n+1}_j =$
風上差分（$a > 0$）	$u^n_j - \nu(u^n_j - u^n_{j-1})$
Lax–Friedrichs	$\frac12(u^n_{j+1} + u^n_{j-1}) - \frac{\nu}{2}(u^n_{j+1} - u^n_{j-1})$
Lax–Wendroff	$u^n_j - \frac{\nu}{2}(u^n_{j+1} - u^n_{j-1}) + \frac{\nu^2}{2}(u^n_{j+1} - 2u^n_j + u^n_{j-1})$
前進時間・中心空間（FTCS）	$u^n_j - \frac{\nu}{2}(u^n_{j+1} - u^n_{j-1})$

風上差分は、情報が来る側（$a > 0$ なら左）の差分を使います。$a < 0$ のときは前進差分 $u^n_j - \nu(u^n_{j+1} - u^n_j)$ を使い、以下の結果は $\nu$ を $|\nu|$ に置き換えて成り立ちます。Lax–Friedrichs は FTCS の $u^n_j$ を両隣の平均に置き換えたもの、Lax–Wendroff は $u(x, t + \Delta t) = u + \Delta t\,u_t + \frac{\Delta t^2}{2}u_{tt} + \cdots$ に $u_t = -au_x$、$u_{tt} = a^2u_{xx}$ を代入して中心差分で置き換えたものです。

定理 7.1（4 つのスキームの安定性と精度）✓

$a \ne 0$、格子比 $\lambda = \Delta t/h$ を固定し、$\nu = a\lambda$ とする。$u_0 \in C^4(\mathbb{R})$ は 4 階までの導関数が有界とし、$U_m = \sup|u_0^{(m)}|$、$u(x,t) = u_0(x - at)$ とする。$s = \sin(\theta/2)$ と書く。

1. 風上差分（$a > 0$）。$|g(\theta)|^2 = 1 - 4\nu(1-\nu)s^2$ で、$\ell^2$ で安定であることと $0 < \nu \le 1$ は同値。このとき $\|u^{n+1}\|_\infty \le \|u^n\|_\infty$ でもある。$|\tau| \le \frac{ah}{2}(1 + \nu)U_2$（1 次精度）。
2. Lax–Friedrichs。$|g(\theta)|^2 = 1 - (1 - \nu^2)\sin^2\theta$ で、$\ell^2$ で安定であることと $|\nu| \le 1$ は同値。このとき $\|u^{n+1}\|_\infty \le \|u^n\|_\infty$ でもある。$|\tau| \le \frac{h}{2}\bigl(a^2\lambda + \frac{1}{\lambda}\bigr)U_2 + \frac{|a|h^2}{6}U_3$（1 次精度）。
3. Lax–Wendroff。$|g(\theta)|^2 = 1 - 4\nu^2(1 - \nu^2)s^4$ で、$\ell^2$ で安定であることと $|\nu| \le 1$ は同値。$|\tau| \le \frac{|a|h^2}{6}(1 + \nu^2)U_3 + \frac{a^2\Delta t\,h^2}{24}U_4$（2 次精度）。
4. FTCS。$|g(\theta)|^2 = 1 + \nu^2\sin^2\theta$ で、すべての $\lambda > 0$ について $\ell^2$ で不安定。$|\tau| \le \frac{a^2\Delta t}{2}U_2 + \frac{|a|h^2}{6}U_3$ なので整合的ではある。

証明（定理 7.1）

増幅因子は $u^n_j = e^{ij\theta}$ を代入して求めます。係数は $\nu$ だけで決まり $\Delta t$ によらないので、定理 5.5 から、安定性は $\max_\theta|g(\theta)| \le 1$ と同値です。打ち切り誤差は、$x$ を固定して $t$ について、また $t$ を固定して $x$ について 定理 5.3［平均値の定理とテイラーの定理］ を使い、補題 6.1 の $|\partial_t^i\partial_x^ju| \le |a|^iU_{i+j}$、$u_t = -au_x$、$u_{tt} = a^2u_{xx}$ を使って評価します。$\nu/\Delta t = a/h$ に注意します。

1. 風上差分. $g = 1 - \nu + \nu e^{-i\theta}$ なので

$$|g|^2 = (1 - \nu + \nu\cos\theta)^2 + \nu^2\sin^2\theta = (1-\nu)^2 + 2\nu(1-\nu)\cos\theta + \nu^2 = 1 - 2\nu(1-\nu)(1 - \cos\theta) = 1 - 4\nu(1-\nu)s^2$$

です。$0 < \nu \le 1$ なら $\nu(1 - \nu) \ge 0$ で $|g| \le 1$、$\nu > 1$ なら $\theta = \pi$ で $|g|^2 = 1 + 4\nu(\nu - 1) > 1$ です（$a > 0$ なので $\nu > 0$）。$0 < \nu \le 1$ のとき、$u^{n+1}_j = (1-\nu)u^n_j + \nu u^n_{j-1}$ は非負の係数（和が $1$）の結合なので $|u^{n+1}_j| \le (1-\nu)\|u^n\|_\infty + \nu\|u^n\|_\infty = \|u^n\|_\infty$ です。打ち切り誤差は

$$\tau = \frac{u(x, t+\Delta t) - u(x,t)}{\Delta t} + a\,\frac{u(x,t) - u(x-h,t)}{h} = \Bigl(u_t + \frac{\Delta t}{2}u_{tt}(x, t')\Bigr) + a\Bigl(u_x - \frac{h}{2}u_{xx}(x', t)\Bigr)$$

（$t'$、$x'$ は Taylor の定理の中間の点）で、$u_t + au_x = 0$ から $|\tau| \le \frac{\Delta t}{2}a^2U_2 + \frac{ah}{2}U_2 = \frac{ah}{2}(\nu + 1)U_2$ です。

2. Lax–Friedrichs. $g = \cos\theta - i\nu\sin\theta$ なので $|g|^2 = \cos^2\theta + \nu^2\sin^2\theta = 1 - (1-\nu^2)\sin^2\theta$ です。$|\nu| \le 1$ なら $|g| \le 1$、$|\nu| > 1$ なら $\theta = \pi/2$ で $|g| = |\nu| > 1$ です。$|\nu| \le 1$ のとき $u^{n+1}_j = \frac{1-\nu}{2}u^n_{j+1} + \frac{1+\nu}{2}u^n_{j-1}$ は非負の係数（和が $1$）の結合なので $\|u^{n+1}\|_\infty \le \|u^n\|_\infty$ です。打ち切り誤差は

$$\tau = \frac{u(x,t+\Delta t) - u}{\Delta t} - \frac{h^2}{2\Delta t}D_+D_-u + aD_0u$$

です。2 階の剰余までの Taylor 展開 $u(x\pm h) = u \pm hu_x + \frac{h^2}{2}u_{xx}(x'_\pm)$ から $|D_+D_-u| \le U_2$、補題 2.3 の 2 から $|D_0u - u_x| \le \frac{h^2}{6}U_3$ なので、$|\tau| \le \frac{\Delta t}{2}a^2U_2 + \frac{h^2}{2\Delta t}U_2 + \frac{|a|h^2}{6}U_3$ で、$\Delta t = \lambda h$ を代入すると主張の形です。

3. Lax–Wendroff. $1 - \cos\theta = 2s^2$ なので $g = 1 - 2\nu^2s^2 - i\nu\sin\theta$ で、$\sin^2\theta = 4s^2(1 - s^2)$ から

$$|g|^2 = (1 - 2\nu^2s^2)^2 + 4\nu^2s^2(1 - s^2) = 1 - 4\nu^2s^2 + 4\nu^4s^4 + 4\nu^2s^2 - 4\nu^2s^4 = 1 - 4\nu^2(1-\nu^2)s^4$$

です。$|\nu| \le 1$ なら $|g| \le 1$、$|\nu| > 1$ なら $\theta = \pi$ で $|g|^2 = 1 + 4\nu^2(\nu^2 - 1) > 1$ です。打ち切り誤差は

$$\tau = \frac{u(x,t+\Delta t) - u}{\Delta t} + aD_0u - \frac{a^2\Delta t}{2}D_+D_-u$$

で、$\frac{u(x,t+\Delta t) - u}{\Delta t} = u_t + \frac{\Delta t}{2}u_{tt} + \frac{\Delta t^2}{6}u_{ttt}(x,t')$、$D_0u = u_x + R_1$、$D_+D_-u = u_{xx} + R_2$（$|R_1| \le \frac{h^2}{6}U_3$、$|R_2| \le \frac{h^2}{12}U_4$。補題 2.3）を代入すると

$$\tau = \underbrace{(u_t + au_x)}_{=0} + \frac{\Delta t}{2}\underbrace{(u_{tt} - a^2u_{xx})}_{=0} + \frac{\Delta t^2}{6}u_{ttt}(x,t') + aR_1 - \frac{a^2\Delta t}{2}R_2$$

です。$|u_{ttt}| \le |a|^3U_3$ と $|a|^3\Delta t^2 = |a|\nu^2h^2$ から $|\tau| \le \frac{|a|h^2}{6}\nu^2U_3 + \frac{|a|h^2}{6}U_3 + \frac{a^2\Delta t\,h^2}{24}U_4$ です。

4. FTCS. $g = 1 - i\nu\sin\theta$ なので $|g|^2 = 1 + \nu^2\sin^2\theta$ で、$\theta = \pi/2$ で $|g| = \sqrt{1 + \nu^2} > 1$ です（$a \ne 0$、$\lambda > 0$ なので $\nu \ne 0$）。したがって不安定です。打ち切り誤差は $\tau = \frac{u(x,t+\Delta t) - u}{\Delta t} + aD_0u = \frac{\Delta t}{2}u_{tt}(x,t') + aR_1$ で、主張の評価が成り立ちます。

∎

FTCS は整合的で CFL 条件も満たしますが不安定なので、定理 4.4 により収束しません。不安定性は格子比を固定したときのものです。$\Delta t = \mu h^2$ のように $\Delta t$ を $h^2$ に比例させれば $\nu^2 = a^2\mu\Delta t$ となって安定になります（演習 10.3）が、時間ステップの数が $h^{-2}$ に比例するので実用的ではありません。

CFL 数 ν = 0.8 のときの増幅因子の大きさ |g(θ)|（0 ≤ θ ≤ π）。破線の水平線が 1。FTCS だけが 1 を超え、θ = π/2 の細かい波が 1 ステップごとに約 1.28 倍に育つ。風上差分・Lax–Friedrichs・Lax–Wendroff は 1 以下で、その分だけ波が減衰する。Lax–Wendroff は θ が小さいところでほとんど減衰せず、θ = π 付近で強く減衰する。

矩形波を実際に運んでみると、風上差分と Lax–Wendroff の違いが目に見えます。次の動画は $\nu = 0.8$、$h = 0.01$ で 1 周期ぶん進めたもので、風上差分は角がなまり（数値粘性）、Lax–Wendroff は角は鋭いものの後ろに振動が出ます。

周期境界で矩形波を 1 周期運ぶ。破線が厳密解。風上差分（左）は 1 次精度で波形がなまり、Lax–Wendroff（右）は 2 次精度で立ち上がりは鋭いが、後ろに振動が出て 1 を超える。

系 7.2（周期境界での収束の次数）✓

$u_0$ は周期 $1$ の $C^4$ 級関数、$h = 1/N$、$\lambda$ は固定で、定理 7.1 の安定条件（風上差分は $0 < \nu \le 1$、Lax–Friedrichs と Lax–Wendroff は $|\nu| \le 1$）が成り立つとする。周期格子の上で $u^0_j = u_0(x_j)$ から計算した $u^n$ は

$$\Bigl(h\sum_{j=0}^{N-1}|u(x_j, t_n) - u^n_j|^2\Bigr)^{1/2} \le t_n\max_{x,t}|\tau(x,t)|$$

を満たす。右辺は風上差分と Lax–Friedrichs で $O(h)$、Lax–Wendroff で $O(h^2)$ である。

証明（系 7.2）

補題 5.3 と $\max_\theta|g| \le 1$ から、$\mathbb{C}^N$ 上で $\|C^n\| \le 1$ です。$w^n_j = u(x_j, t_n)$ とおくと $\tau^n_j = (w^{n+1}_j - (Cw^n)_j)/\Delta t = \tau(x_j, t_n)$ で、$hN = 1$ から $\|\tau^n\|_h \le \max|\tau|$ です。補題 4.1 を $K = 1$、$e^0 = 0$ で使えば結論を得ます。次数は 定理 7.1 の $\tau$ の評価から従います（周期関数の導関数は有界です）。

∎

例 7.3（CFL 数 1 では風上差分が厳密になる）

$\nu = 1$ のとき風上差分は $u^{n+1}_j = u^n_{j-1}$ です。$a\Delta t = h$ なので、厳密解も $u(x_j, t_{n+1}) = u_0(x_j - at_{n+1}) = u_0(x_{j-1} - at_n) = u(x_{j-1}, t_n)$ を満たし、数値解は格子点で厳密解に一致します。特性曲線がちょうど格子点を通るからです。Lax–Friedrichs と Lax–Wendroff も $\nu = 1$ では同じ $u^{n+1}_j = u^n_{j-1}$ になります。ただし、係数 $a$ が場所によって変わる問題ではすべての点で $\nu = 1$ にはできないので、この性質は当てにできません。

8. 数値的な散逸と分散、Godunov の次数の壁

定義 8.1（数値的な散逸と分散）↗

移流方程式の厳密解では、波 $e^{i\xi x}$ は 1 ステップで $e^{-ia\xi\Delta t} = e^{-i\nu\theta}$ 倍（$\theta = h\xi$）になり、振幅は変わらない。スキームの増幅因子 $g(\theta)$ について、$|g(\theta)| < 1$ であること（振幅が減ること）を数値的な散逸、$\arg g(\theta) \ne -\nu\theta$ であること（位相の進みがずれ、波数によって伝わる速さが変わること）を数値的な分散という。

定理 7.1 の $|g|^2$ を見ると、風上差分は $1 - |g|^2$ が $s^2 \approx \theta^2/4$ に比例し、Lax–Wendroff は $s^4$ に比例します。つまり風上差分は滑らかな成分まで減衰させ、Lax–Wendroff は細かい成分だけを弱く減衰させます。この違いは、スキームが「実際に」近似している方程式を調べるとはっきりします。

次の命題の考え方は、「移流方程式の近似として作ったスキームは、実は少しずれた別の方程式をより高い精度で解いている」というものです。風上差分の場合、そのずれは拡散項 $\kappa w_{xx}$ になります。拡散は山を低く広げるので、矩形波の角がなまるのは、このためです。Lax–Wendroff の場合、ずれは 3 階微分の項（波数によって速さを変える分散）になります。

命題 8.2（修正方程式：数値粘性と分散）✓

$a > 0$、$\nu = a\Delta t/h$ とする。

1. $\kappa = \frac{ah}{2}(1 - \nu)$ とし、$w$ は $w_t + aw_x = \kappa w_{xx}$ を満たす滑らかな関数で、以下に現れる導関数が有界であるとする。風上差分に $w$ を代入した残差 $\tau_w = \frac{w(x,t+\Delta t) - w}{\Delta t} + aD_-w$ は

$$|\tau_w| \le \Bigl(a\kappa\Delta t + \frac{ah^2}{6}\Bigr)\sup|w_{xxx}| + \frac{\kappa^2\Delta t}{2}\sup|w_{xxxx}| + \frac{\Delta t^2}{6}\sup|w_{ttt}|$$

を満たす。$\Delta t/h$ を固定すれば右辺は $O(h^2)$ である。 2. $u_0 \in C^5$ は 5 階までの導関数が有界とし、$u = u_0(x - at)$ とする。Lax–Wendroff の局所打ち切り誤差は

$$\tau = \frac{ah^2}{6}(1 - \nu^2)u_{xxx} + \rho,\qquad|\rho| \le \frac{a^4\Delta t^3}{24}U_4 + \frac{a^2\Delta t\,h^2}{24}U_4 + \frac{ah^4}{120}U_5$$

を満たす。

証明（命題 8.2）

1. 3 階の剰余までの Taylor 展開から

$$\frac{w(x,t+\Delta t) - w}{\Delta t} = w_t + \frac{\Delta t}{2}w_{tt} + \frac{\Delta t^2}{6}w_{ttt}(x,t'),\qquad D_-w = w_x - \frac{h}{2}w_{xx} + \frac{h^2}{6}w_{xxx}(x',t)$$

です。$w_t = -aw_x + \kappa w_{xx}$ を $t$ と $x$ で微分して

$$w_{tt} = -aw_{xt} + \kappa w_{xxt} = -a(-aw_{xx} + \kappa w_{xxx}) + \kappa(-aw_{xxx} + \kappa w_{xxxx}) = a^2w_{xx} - 2a\kappa w_{xxx} + \kappa^2w_{xxxx}$$

です。代入すると

$$\tau_w = \bigl(w_t + aw_x\bigr) + \Bigl(\frac{a^2\Delta t}{2} - \frac{ah}{2}\Bigr)w_{xx} - a\kappa\Delta t\,w_{xxx} + \frac{\kappa^2\Delta t}{2}w_{xxxx} + \frac{\Delta t^2}{6}w_{ttt}(x,t') + \frac{ah^2}{6}w_{xxx}(x',t)$$

で、$\frac{a^2\Delta t}{2} - \frac{ah}{2} = -\frac{ah}{2}(1 - \nu) = -\kappa$ なので、最初の 2 項の和は $w_t + aw_x - \kappa w_{xx} = 0$ です。残りの項の絶対値を評価すれば主張を得ます。$\kappa = O(h)$、$\Delta t = O(h)$ なので各項は $O(h^2)$ です。

2. 定理 7.1 の 3 の証明より 1 つ高い次数まで展開します。$\frac{u(x,t+\Delta t) - u}{\Delta t} = u_t + \frac{\Delta t}{2}u_{tt} + \frac{\Delta t^2}{6}u_{ttt} + \frac{\Delta t^3}{24}u_{tttt}(x,t')$ です。$D_0u$ は、5 階の剰余までの展開 $u(x \pm h) = \sum_{k=0}^{4}\frac{(\pm h)^k}{k!}\partial_x^ku \pm \frac{h^5}{120}\partial_x^5u(x'_\pm)$ の差をとって $D_0u = u_x + \frac{h^2}{6}u_{xxx} + R_1'$、$|R_1'| \le \frac{h^4}{120}U_5$ です。$D_+D_-u = u_{xx} + R_2$、$|R_2| \le \frac{h^2}{12}U_4$ は 補題 2.3 の 3 です。$u_t + au_x = 0$、$u_{tt} = a^2u_{xx}$、$u_{ttt} = -a^3u_{xxx}$ を使うと

$$\tau = -\frac{a^3\Delta t^2}{6}u_{xxx} + \frac{ah^2}{6}u_{xxx} + \frac{\Delta t^3}{24}u_{tttt}(x,t') + aR_1' - \frac{a^2\Delta t}{2}R_2$$

で、$a^3\Delta t^2 = a\nu^2h^2$ から最初の 2 項は $\frac{ah^2}{6}(1-\nu^2)u_{xxx}$ です。残りが $\rho$ で、$|u_{tttt}| \le a^4U_4$ から主張の評価を得ます。

∎

1 は「風上差分は、移流方程式を 1 次の精度で近似するが、移流拡散方程式 $w_t + aw_x = \kappa w_{xx}$ を 2 次の精度で近似する」と読めます。$\kappa = \frac{ah}{2}(1 - \nu) \ge 0$ が数値粘性で、不連続をなまらせる原因です。$\nu > 1$ では $\kappa < 0$ となり、逆向きの拡散方程式（非適切な問題）を近似することになります。これは不安定性の別の見方です。2 は、Lax–Wendroff の誤差の先頭が 3 階微分であることを示しています。$u_t + au_x = -\frac{ah^2}{6}(1-\nu^2)u_{xxx}$ は波数 $\xi$ の波を速さ $a\bigl(1 - \frac{h^2\xi^2}{6}(1 - \nu^2)\bigr)$ で運ぶ方程式で、細かい波ほど遅れます。不連続な初期値は全波数の成分を含むので、遅れた成分が不連続の後ろに振動として残ります。

定義 8.3（単調性を保つスキーム）↗

スキームが単調性を保つとは、$u^n_j \le u^n_{j+1}$（すべての $j$）ならば $u^{n+1}_j \le u^{n+1}_{j+1}$（すべての $j$）となることをいう。

単調性を保つスキームは、単調な初期値に新しい極値（振動）を作りません。

定理 8.4（Godunov の次数の壁）✓

$a \ne 0$、$\lambda = \Delta t/h$ を固定し、実数の係数 $a_k$（$-l \le k \le r$）が $h$ によらない定数係数の陽的スキーム $u^{n+1}_j = \sum_ka_ku^n_{j+k}$ を考える。

1. スキームが単調性を保つことと、すべての $k$ で $a_k \ge 0$ であることは同値である。
2. 3 階までの導関数が有界なすべての $u_0 \in C^3(\mathbb{R})$ について、移流方程式 $u_t + au_x = 0$ に対する局所打ち切り誤差が $\sup|\tau| \le ch^2$（$c$ は $u_0$ による定数）を満たす（2 次精度）ならば、$\sum_ka_k = 1$、$\sum_kka_k = -\nu$、$\sum_kk^2a_k = \nu^2$ である。
3. 1 と 2 の両方が成り立つのは、$\nu = a\lambda$ が整数で、スキームが $u^{n+1}_j = u^n_{j-\nu}$（厳密解そのもの）である場合に限る。したがって、それ以外の単調性を保つ線形スキームは 2 次精度になりえない。

証明（定理 8.4）

1. $u^{n+1}_{j+1} - u^{n+1}_j = \sum_ka_k(u^n_{j+k+1} - u^n_{j+k})$ です。$a_k \ge 0$ で $u^n$ が単調増加なら、右辺の各項は非負です。逆に、単調性を保つとし、$u^n_i = 0$（$i \le 0$）、$u^n_i = 1$（$i \ge 1$）という単調増加な格子関数をとります。$u^n_{j+k+1} - u^n_{j+k}$ は $j + k = 0$ のときだけ $1$ で他は $0$ なので、$u^{n+1}_{j+1} - u^{n+1}_j = a_{-j}$（$-l \le -j \le r$ のとき。それ以外は $0$）です。これが非負なので、すべての $k$ で $a_k \ge 0$ です。

2. $y = x - at$ とおくと、$u(x + kh, t) = u_0(y + kh)$、$u(x, t + \Delta t) = u_0(y - a\Delta t) = u_0(y - \nu h)$ です。3 階の剰余までの Taylor 展開 $u_0(y + \sigma h) = u_0(y) + \sigma hu_0'(y) + \frac{\sigma^2h^2}{2}u_0''(y) + \frac{\sigma^3h^3}{6}u_0'''(c_\sigma)$ を $\sigma = -\nu$ と $\sigma = k$ で使うと

$$\Delta t\,\tau = u_0(y - \nu h) - \sum_ka_ku_0(y + kh) = m_0u_0(y) + m_1hu_0'(y) + \frac{m_2h^2}{2}u_0''(y) + R,$$$$m_0 = 1 - \sum_ka_k,\quad m_1 = -\nu - \sum_kka_k,\quad m_2 = \nu^2 - \sum_kk^2a_k,\quad|R| \le \frac{h^3}{6}\Bigl(|\nu|^3 + \sum_k|a_k||k|^3\Bigr)U_3$$

です。仮定から $|\Delta t\,\tau| \le c\lambda h^3$ なので、$|m_0u_0(y) + m_1hu_0'(y) + \frac{m_2h^2}{2}u_0''(y)| \le c'h^3$（$c'$ は $h$ によらない）です。$u_0 = 1$ とすると $|m_0| \le c'h^3$ で、$m_0$ は $h$ によらないので $h \to 0$ として $m_0 = 0$ です。次に $u_0 = \sin x$、$y = 0$ とすると $u_0(0) = 0$、$u_0'(0) = 1$、$u_0''(0) = 0$ で $|m_1|h \le c'h^3$、したがって $m_1 = 0$ です。最後に $u_0 = \cos x$、$y = 0$ とすると $u_0'(0) = 0$、$u_0''(0) = -1$ で、$m_0 = 0$ と合わせて $\frac{|m_2|h^2}{2} \le c'h^3$、したがって $m_2 = 0$ です。

3. $a_k \ge 0$ と 2 の 3 つの等式が成り立つとします。Cauchy–Schwarz の不等式から

$$\nu^2 = \Bigl(\sum_kk\sqrt{a_k}\cdot\sqrt{a_k}\Bigr)^2 \le \Bigl(\sum_kk^2a_k\Bigr)\Bigl(\sum_ka_k\Bigr) = \nu^2\cdot1$$

で、等号が成り立っています。Cauchy–Schwarz の不等式の等号は、ベクトル $(k\sqrt{a_k})_k$ と $(\sqrt{a_k})_k$ が比例するときに限ります。$\sum_ka_k = 1$ から $(\sqrt{a_k})_k \ne 0$ なので、定数 $\beta$ があって $k\sqrt{a_k} = \beta\sqrt{a_k}$（すべての $k$）、すなわち $a_k \ne 0$ となる $k$ はすべて $k = \beta$ です。したがって $a_k \ne 0$ となる $k$ はただ 1 つ $k_0$ で、$a_{k_0} = 1$、$k_0 = \sum_kka_k = -\nu$ です。$\nu = -k_0$ は整数で、スキームは $u^{n+1}_j = u^n_{j-\nu}$ です。$a\Delta t = \nu h$ なので、これは厳密解の関係 $u(x_j, t_{n+1}) = u(x_j - \nu h, t_n)$ そのものです。

∎

2 次精度と「振動を作らない」は、線形スキームでは両立しません。風上差分（$a_{-1} = \nu$、$a_0 = 1 - \nu$）は単調性を保ち 1 次精度、Lax–Wendroff（$a_1 = \frac{\nu(\nu - 1)}{2} < 0$、$0 < \nu < 1$）は 2 次精度で単調性を保ちません。

注意 8.5（この先：非線形の保存則と高解像度スキーム）

Godunov の壁を越えるには、スキームを非線形にする（解の滑らかさに応じて係数を変える）しかありません。流束制限関数を使う TVD スキームや ENO/WENO スキームがそれで、滑らかなところでは 2 次以上、不連続の近くでは 1 次の風上差分に切り替わります。非線形の保存則 $u_t + f(u)_x = 0$ では、滑らかな初期値からも有限時間で衝撃波ができ、弱解の一意性にエントロピー条件が要ります。これらの理論は証明が長く、この章の枠組み（線形の同等定理）の外にあるので扱いません。LeVeque『Numerical Methods for Conservation Laws』（参考文献）が標準的な入門書です。

9. 数値実験

周期 $1$ の移流方程式 $u_t + u_x = 0$ を、$\nu = 0.8$ で時刻 $T = 1$ まで解きます。厳密解は $u_0(x - t)$ で、$T = 1$ でちょうど 1 周して初期値に戻ります。滑らかな初期値 $e^{\sin 2\pi x}$ と、角のある初期値（幅 $1/2$ の三角形）で、$\ell^2$ 誤差と、$N$ を 2 倍にしたときの観測次数 $\log_2(e_N/e_{2N})$ を表にします。

# Python 3.9.6 / NumPy 2.0.2（arm64）
import numpy as np

def step(u, nu, scheme):
    um, up = np.roll(u, 1), np.roll(u, -1)      # u_{j-1}, u_{j+1}（周期境界）
    if scheme == "upwind":
        return u - nu * (u - um)
    if scheme == "LF":
        return 0.5 * (up + um) - 0.5 * nu * (up - um)
    if scheme == "LW":
        return u - 0.5 * nu * (up - um) + 0.5 * nu**2 * (up - 2 * u + um)
    return u - 0.5 * nu * (up - um)             # FTCS

def solve(N, nu, scheme, u0, T=1.0):
    h = 1.0 / N
    n = round(T / (nu * h))                     # a = 1, dt = nu * h
    x = h * np.arange(N)
    u = u0(x)
    for _ in range(n):
        u = step(u, nu, scheme)
    err = u - u0((x - n * nu * h) % 1.0)        # 厳密解は平行移動
    return np.sqrt(h * np.sum(err**2))

smooth = lambda x: np.exp(np.sin(2 * np.pi * x))
hat = lambda x: np.maximum(0.0, 1 - 4 * np.abs(x - 0.5))
for name, u0 in [("smooth", smooth), ("hat", hat)]:
    print(name)
    for s in ["upwind", "LF", "LW", "FTCS"]:
        e = [solve(N, 0.8, s, u0) for N in (40, 80, 160, 320, 640)]
        o = [np.log2(e[i] / e[i + 1]) for i in range(4)]
        print(f"{s:7s}" + " ".join(f"{v:9.2e}" for v in e)
              + "  |" + " ".join(f"{v:6.2f}" for v in o))

smooth
upwind  9.97e-02  5.29e-02  2.73e-02  1.39e-02  7.00e-03  |  0.91   0.95   0.98   0.99
LF      1.97e-01  1.10e-01  5.90e-02  3.06e-02  1.56e-02  |  0.84   0.90   0.95   0.97
LW      1.77e-02  4.50e-03  1.13e-03  2.82e-04  7.06e-05  |  1.98   2.00   2.00   2.00
FTCS    1.27e+00  3.32e-01  4.85e+04  2.47e+26  1.61e+69  |  1.94 -17.16 -72.11 -142.23
hat
upwind  6.24e-02  3.68e-02  2.17e-02  1.29e-02  7.66e-03  |  0.76   0.76   0.75   0.75
LF      1.09e-01  6.69e-02  3.98e-02  2.37e-02  1.41e-02  |  0.71   0.75   0.75   0.75
LW      2.85e-02  1.43e-02  7.24e-03  3.68e-03  1.87e-03  |  1.00   0.98   0.98   0.98
FTCS    1.93e+03  1.11e+08  1.88e+18  1.70e+39  4.69e+81  |-15.81 -33.98 -69.62 -140.99

検算. (a) 参照解は閉じた形（平行移動）です。(b) 滑らかな初期値での観測次数は、風上差分と Lax–Friedrichs で $1$、Lax–Wendroff で $2$ に近づき、系 7.2 と合います。FTCS は 定理 7.1 の 4 のとおり発散します。$N = 640$ ではステップ数が $n = 800$ で、$\max|g|^n = (1 + 0.8^2)^{400} \approx 10^{86}$ です。丸め誤差（$10^{-16}$ 程度）に含まれる $\theta \approx \pi/2$ の成分がこの倍率で増えるとすると $10^{70}$ 程度になり、表の $1.6\times10^{69}$ と桁が合います。(c) 別の方法として、同じコードで $\nu = 1$、$N = 160$ とすると、風上差分・Lax–Friedrichs・Lax–Wendroff の誤差は $2.1\times10^{-16}$、$2.8\times10^{-16}$、$2.9\times10^{-16}$ で、例 7.3 のとおり丸めの大きさです。(d) 修正方程式の確認として、初期値を $\sin 2\pi x$、$N = 320$、$\nu = 0.8$ とすると、風上差分の解と厳密解 $\sin 2\pi(x - 1)$ の $\ell^2$ 距離は $8.67\times10^{-3}$、移流拡散方程式の解 $e^{-4\pi^2\kappa}\sin2\pi(x - 1)$（$\kappa = \frac{h}{2}(1 - \nu)$）との距離は $3.38\times10^{-5}$ で、命題 8.2 の 1 のとおり後者のほうが 2 桁以上近くなります。

例 9.1（滑らかでない初期値では次数が落ちる）

三角形の初期値は連続ですが $C^1$ 級ではないので、系 7.2 の仮定を満たしません。表では、風上差分と Lax–Friedrichs の観測次数は約 $0.75$、Lax–Wendroff は約 $1$ で、どれも収束はしています（定理 4.2 の前半は一般の初期値での収束を保証します）が、次数は公称の値より下がっています。精度 $p$ 次のスキームを滑らかさの足りない初期値に使ったときの次数は、Brenner–Thomée–Wahlbin（参考文献）が Besov 空間を使って調べています。ここでは、次数つきの評価には解の滑らかさが本当に必要であることを確認するにとどめます。

次に、矩形波を初期値にして、CFL 数を変えたときの風上差分の最大値と、$\nu = 0.5$ での全変動 $\sum_j|u_j - u_{j-1}|$ を調べます。

# Python 3.9.6 / NumPy 2.0.2（arm64）
import numpy as np

N = 100
x = np.arange(N) / N
sq = np.where((x >= 0.25) & (x < 0.75), 1.0, 0.0)       # 矩形波
tv = lambda u: np.sum(np.abs(u - np.roll(u, 1)))        # 全変動

def upwind(u, nu):
    return u - nu * (u - np.roll(u, 1))

def lw(u, nu):
    um, up = np.roll(u, 1), np.roll(u, -1)
    return u - 0.5 * nu * (up - um) + 0.5 * nu**2 * (up - 2 * u + um)

print("風上差分の最大値 max|u^n| と増幅因子の最大値の n 乗")
for nu in (0.5, 1.0, 1.1):
    u, g = sq.copy(), max(1.0, abs(1 - 2 * nu))
    row = []
    for n in range(1, 201):
        u = upwind(u, nu)
        if n in (50, 100, 200):
            row.append(f"n={n}: {np.max(np.abs(u)):9.3e} ({g**n:9.3e})")
    print(f"nu={nu}: " + "  ".join(row))

print("全変動（nu = 0.5）")
for name, f in [("upwind", upwind), ("LW", lw)]:
    u, row = sq.copy(), [tv(sq)]
    for n in range(1, 101):
        u = f(u, 0.5)
        if n in (1, 10, 100):
            row.append(tv(u))
    print(f"{name:7s}" + " ".join(f"{v:8.4f}" for v in row)
          + f"   min={u.min():8.4f} max={u.max():8.4f}")

風上差分の最大値 max|u^n| と増幅因子の最大値の n 乗
nu=0.5: n=50: 1.000e+00 (1.000e+00)  n=100: 1.000e+00 (1.000e+00)  n=200: 9.996e-01 (1.000e+00)
nu=1.0: n=50: 1.000e+00 (1.000e+00)  n=100: 1.000e+00 (1.000e+00)  n=200: 1.000e+00 (1.000e+00)
nu=1.1: n=50: 9.574e+02 (9.100e+03)  n=100: 6.020e+06 (8.282e+07)  n=200: 3.524e+14 (6.859e+15)
全変動（nu = 0.5）
upwind   2.0000   2.0000   2.0000   2.0000   min=  0.0000 max=  1.0000
LW       2.0000   2.5000   2.8801   3.4455   min= -0.2041 max=  1.2041

検算. (a) $\nu \le 1$ では最大値は $1$ を超えず、定理 7.1 の 1 の $\|u^{n+1}\|_\infty \le \|u^n\|_\infty$ と合います。(b) $\nu = 1.1$ では $\max_\theta|g| = |g(\pi)| = |1 - 2\nu| = 1.2$ で、補題 5.3 から $\|u^n\|_h \le 1.2^n\|u^0\|_h$ です。観測された最大値はこの上界より小さく、1 ステップあたりの増え方は $n = 50 \to 100$ で $(6.020\times10^6/957.4)^{1/50} = 1.191$、$n = 100 \to 200$ で $1.196$ と、$1.2$ に下から近づいています。この矩形波は $\theta = \pi$ の成分 $\sum_j(-1)^ju_j$ がちょうど $0$ なので、$\pi$ に近い波数の成分が先に育つためです。(c) 風上差分の全変動は増えず、最小値と最大値も $[0, 1]$ に収まっています（単調性を保つこと。定理 8.4 の 1）。Lax–Wendroff の全変動は増え、$[0,1]$ をはみ出す振動（$-0.20$ と $1.20$）が現れます。これが 命題 8.2 の 2 の分散による振動です。

10. 演習

演習 10.1易

$a > 0$ のとき、風下側の差分を使うスキーム $u^{n+1}_j = u^n_j - \nu(u^n_{j+1} - u^n_j)$（$\nu = a\Delta t/h > 0$）は、どんな $\nu$ でも不安定であることを、(a) 増幅因子、(b) CFL 条件の 2 通りで示せ。

解答

(a) $g(\theta) = 1 + \nu - \nu e^{i\theta}$ で、$g(\pi) = 1 + 2\nu > 1$ です。$g$ は $\Delta t$ によらないので、定理 5.5 から不安定です。

(b) このスキームは $u^n_j$ と $u^n_{j+1}$ だけを使うので $l = 0$、$r = 1$ です。定理 6.3 の必要条件は $0 \le -a\lambda \le 1$ ですが、$-a\lambda = -\nu < 0$ なので成り立ちません。数値的な依存領域 $[x, x + t/\lambda]$ は点の右側にあり、特性曲線の足 $x - at$ は左側にあります。したがって収束せず、整合的なので 定理 4.4 により不安定です（整合性は風上差分と同じ計算で $|\tau| \le \frac{ah}{2}(1+\nu)U_2$）。

演習 10.2標準

$a > 0$ とする。陰的な風上差分 $u^{n+1}_j + \nu(u^{n+1}_j - u^{n+1}_{j-1}) = u^n_j$ は、周期格子の上で、すべての $\nu > 0$ について $\|u^{n+1}\|_h \le \|u^n\|_h$ を満たすことを示せ。CFL 条件に反しないのはなぜか。

解答

左辺の作用素を $B = (1 + \nu)I - \nu S^{-1}$ と書きます（$S$ は 補題 5.1 の平行移動）。補題 5.3 の証明の基底 $\varphi^{(m)}$ について $B\varphi^{(m)} = b(\theta_m)\varphi^{(m)}$、$b(\theta) = 1 + \nu - \nu e^{-i\theta}$ です。

$$|b(\theta)|^2 = (1 + \nu - \nu\cos\theta)^2 + \nu^2\sin^2\theta = 1 + 2\nu(1+\nu)(1 - \cos\theta) = 1 + 4\nu(1+\nu)\sin^2(\theta/2) \ge 1$$

なので $B$ は正則で、$u^{n+1} = B^{-1}u^n$ は $\varphi^{(m)}$ の成分を $1/b(\theta_m)$ 倍します。$|1/b| \le 1$ と $\{\varphi^{(m)}\}$ の正規直交性から $\|u^{n+1}\|_h \le \|u^n\|_h$ です。

$B^{-1}$ は一般にすべての成分が $0$ でない行列で、$u^{n+1}_j$ はすべての $u^n_i$ に依存します。数値的な依存領域が全区間なので、どんな $\nu$ でも特性曲線の足を含みます。定理 6.3 は陽的スキームについての定理で、陰的スキームには制限を与えません。

演習 10.3標準

FTCS スキームを、$\Delta t = \mu h^2$（$\mu > 0$ は定数）として使うと $\ell^2$ で安定になることを示せ。このときの打ち切り誤差の次数と、時刻 $T$ までのステップ数を $h$ で表せ。

解答

$\nu = a\Delta t/h = a\mu h$ なので $\nu^2 = a^2\mu^2h^2 = a^2\mu\Delta t$ です。定理 7.1 の 4 から $|g(\theta)|^2 = 1 + \nu^2\sin^2\theta \le 1 + a^2\mu\Delta t$ で、$\sqrt{1 + s} \le 1 + s/2$ から $\max_\theta|g| \le 1 + \frac{a^2\mu}{2}\Delta t$ です。定理 5.5 の条件 2 が $C_0 = a^2\mu/2$ で成り立つので安定で、$\|C^n\| \le e^{a^2\mu T/2}$ です。打ち切り誤差は $|\tau| \le \frac{a^2\Delta t}{2}U_2 + \frac{|a|h^2}{6}U_3 = O(h^2)$、ステップ数は $T/\Delta t = T/(\mu h^2)$ です。格子比を固定する Lax–Wendroff は同じ $O(h^2)$ の精度をステップ数 $O(h^{-1})$ で達成するので、こちらを使う理由はありません。

演習 10.4難

$u^n_{j-1}, u^n_j, u^n_{j+1}$ だけを使う定数係数の陽的スキーム（係数は $\nu$ だけによる）で、移流方程式に対して 定理 8.4 の 2 の意味で 2 次精度をもつものは Lax–Wendroff に限ることを示せ。また、$0 < \nu < 1$ で負の係数が現れることを確かめよ。

解答

定理 8.4 の 2 から、係数は $a_{-1} + a_0 + a_1 = 1$、$a_1 - a_{-1} = -\nu$、$a_1 + a_{-1} = \nu^2$ を満たします。後の 2 式から $a_1 = \frac{\nu^2 - \nu}{2}$、$a_{-1} = \frac{\nu^2 + \nu}{2}$ で、最初の式から $a_0 = 1 - \nu^2$ です。Lax–Wendroff の式を $u^n_{j\pm1}$、$u^n_j$ について整理すると、$u^n_{j+1}$ の係数は $-\frac{\nu}{2} + \frac{\nu^2}{2}$、$u^n_{j-1}$ の係数は $\frac{\nu}{2} + \frac{\nu^2}{2}$、$u^n_j$ の係数は $1 - \nu^2$ で、一致します。$0 < \nu < 1$ では $a_1 = \frac{\nu(\nu-1)}{2} < 0$ なので、定理 8.4 の 1 から単調性を保ちません。これは同じ定理の 3 と整合しています。

参考文献

• P. D. Lax, R. D. Richtmyer, “Survey of the stability of linear finite difference equations”, Communications on Pure and Applied Mathematics 9 (1956), 267–293 — 同等定理の原論文。
• R. D. Richtmyer, K. W. Morton『Difference Methods for Initial-Value Problems』2nd ed., Interscience, 1967 — 適切性・整合性・安定性の枠組みと同等定理、von Neumann 条件。
• R. Courant, K. Friedrichs, H. Lewy, “Über die partiellen Differenzengleichungen der mathematischen Physik”, Mathematische Annalen 100 (1928), 32–74 — 依存領域の条件の原論文。
• J. C. Strikwerda『Finite Difference Schemes and Partial Differential Equations』2nd ed., SIAM, 2004 — Fourier 解析による安定性、散逸と分散、同等定理の証明。
• B. Gustafsson, H.-O. Kreiss, J. Oliger『Time Dependent Problems and Difference Methods』Wiley, 1995 — Kreiss の行列定理と境界のあるスキームの安定性。
• R. J. LeVeque『Numerical Methods for Conservation Laws』2nd ed., Birkhäuser, 1992 — 修正方程式、単調スキームと Godunov の定理、非線形の保存則。
• P. Brenner, V. Thomée, L. B. Wahlbin『Besov Spaces and Applications to Difference Methods for Initial Value Problems』Lecture Notes in Mathematics 434, Springer, 1975 — 滑らかでない初期値での収束の次数。