# 楕円型・放物型方程式の差分法：離散最大値原理と θ 法

> Poisson 方程式の 5 点差分について離散最大値原理から最大値ノルムの安定性と O(h²) の収束を証明し、熱方程式の θ 法（陽的 Euler・陰的 Euler・Crank–Nicolson）の安定条件、最大値原理、誤差の次数、高波数の減衰を導く。
> https://rikai.mugen-giken.com/computer-science/pde-discretization/elliptic-parabolic-difference

## 0. この記事の要点

金属板の温度分布や、棒を伝わる熱のように、連続的に広がる量をコンピュータで扱うには、場所と時間を格子の点に区切って近似します。この記事で調べるのは、そうして作った近似が本当に元の方程式の解に近づくのか、近づくために時間の刻み幅をどう選べばよいのか、という素朴な疑問です。読み終えると、板の定常な温度を求める計算と、熱の広がりを時間順に追う計算のそれぞれで、近似の誤差が格子幅に応じてどれだけ小さくなるか、どの刻み幅なら計算が壊れないかを説明できるようになります。

- Poisson 方程式 $-\Delta u = f$ の 5 点差分は、連続の場合と同じ形の**最大値原理**を満たします。そこから、境界で $0$ の格子関数について $\|u\|_\infty \le \frac18\|\Delta_hu\|_\infty$（安定性）が出て、打ち切り誤差（差分近似が元の微分とずれる大きさ）$O(h^2)$（[整合性](/computer-science/pde-discretization/finite-difference-convergence#def-consistency)）と合わせて $\|u - u_h\|_\infty \le \frac{h^2}{96}\bigl(\|\partial_x^4u\|_\infty + \|\partial_y^4u\|_\infty\bigr)$ が従います。
- 離散 Laplace 作用素の固有値は閉じた形で分かり、最小固有値は $2\pi^2$ に近い値です。これは $\ell^2$ ノルムでの安定性と、連立方程式を反復法で解くときの[条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number) $O(h^{-2})$ を与えます。
- 熱方程式 $u_t = u_{xx}$ の **θ 法**は、$\theta \ge 1/2$ なら無条件に、$\theta < 1/2$ なら格子比（時間刻みと格子幅の 2 乗の比）$\mu = \Delta t/h^2 \le \frac{1}{2(1 - 2\theta)}$ のときに限り $\ell^2$ で安定です。陽的 Euler（$\theta = 0$）では $\Delta t \le h^2/2$ です。
- $\mu(1 - \theta) \le 1/2$ なら θ 法は離散最大値原理を満たします。Crank–Nicolson（$\theta = 1/2$）は $\ell^2$ では無条件に安定ですが、最大値原理が保証されるのは $\mu \le 1$ のときです。
- 誤差は $\theta \ne 1/2$ で $O(\Delta t + h^2)$、$\theta = 1/2$ で $O(\Delta t^2 + h^2)$ です。Crank–Nicolson は高波数の成分をほとんど減衰させない（1 ステップで波の大きさが何倍になるかを表す増幅因子が $-1$ に近い）ので、角や跳びのある初期値では振動が残ります。最初の数ステップを陰的 Euler にすると消えます。

<Aside type="tip" title="この記事で使う事前知識">

- **整合性**: 格子幅を 0 に近づけると、差分近似と元の微分方程式とのずれ（打ち切り誤差）が 0 に近づく、という性質。→ [整合性](/computer-science/pde-discretization/finite-difference-convergence#def-consistency)
- **中心差分**: $u(x+h)$ と $u(x-h)$ を使って $x$ での微分を近似する差分。→ [中心差分](/computer-science/pde-discretization/finite-difference-convergence#def-difference-quotients)
- **調和関数・劣調和関数**: $\Delta u = 0$ を満たす関数と、$-\Delta u \le 0$ を満たす関数。最大値原理は後者の最大が境界でとられるという主張です。→ [調和関数](/mathematics/pde/laplace-equation#def-harmonic)
- **条件数**: 連立方程式の解が、データの小さな変化にどれだけ敏感かを表す数。→ [条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number)
- **正定値**: 零でないどんなベクトルに対しても、$x^{\top}Ax > 0$ となる対称行列の性質。→ [正定値](/mathematics/calculus/multivariable-differentiation#def-definite)

</Aside>

## 1. 動機

前章では、発展方程式の差分スキームについて「整合性と安定性から収束が従う」ことを示しました。この章では、同じ枠組みを 2 つの基本的な方程式に適用します。時間によらない Poisson 方程式（楕円型）と、熱方程式（放物型）です。

定常問題には時間ステップがないので、安定性の意味が変わります。差分方程式 $-\Delta_hu_h = f$ は連立一次方程式で、安定性とは「逆行列のノルムが格子幅によらず有界」ということです。誤差 $e = u - u_h$ は打ち切り誤差 $\tau$ を右辺とする同じ差分方程式 $-\Delta_he = \tau$ を満たすので、$\|e\| \le \|(-\Delta_h)^{-1}\|\|\tau\|$ となります。逆行列のノルムを抑える道具が、最大値ノルムでは**離散最大値原理**、$\ell^2$ ノルムでは固有値です。

熱方程式では、もう 1 つ問題が現れます。陽的な方法の安定条件 $\Delta t \le h^2/2$ は、格子を 2 倍細かくすると時間ステップを 4 倍細かくすることを要求します。これを避けるには陰的な方法を使いますが、陰的な方法の中でも、高波数の成分の扱いに違いがあります。θ 法は、陽的 Euler・陰的 Euler・Crank–Nicolson を 1 つのパラメータでつなぐ族で、この違いを 1 つの式で調べられます。

方程式そのものの性質（[Laplace 方程式](/mathematics/pde/laplace-equation) の最大値原理、[熱方程式](/mathematics/pde/heat-equation) の最大値原理と平滑化）は既存の章で証明されています。この章で示すのは、差分スキームがそれらの性質をどこまで受け継ぐか、です。

## 2. 5 点差分と離散最大値原理

<Definition id="def-discrete-laplacian" title="格子領域と離散 Laplace 作用素">
$d \in \{1, 2\}$、$h > 0$ とし、$\boldsymbol{e}_1, \ldots, \boldsymbol{e}_d$ を $\mathbb{R}^d$ の標準基底とする。格子 $h\mathbb{Z}^d$ の点 $P$ の**隣接点**を $P \pm h\boldsymbol{e}_k$（$k = 1, \ldots, d$）とし、その集合を $N(P)$ と書く（$2d$ 個）。$h\mathbb{Z}^d$ の空でない有限部分集合 $\Omega_h$ を**内部の格子点**の集合とし、
$$
\Gamma_h = \{Q \notin \Omega_h : Q \in N(P)\ \text{となる}\ P \in \Omega_h\ \text{がある}\},\qquad \bar\Omega_h = \Omega_h \cup \Gamma_h
$$
を**境界の格子点**の集合、**格子領域**という。$\bar\Omega_h$ 上の実数値関数 $u$ に対し
$$
\Delta_hu(P) = \frac{1}{h^2}\sum_{Q \in N(P)}\bigl(u(Q) - u(P)\bigr)\qquad(P \in \Omega_h)
$$
を**離散 Laplace 作用素**という。$d = 1$ では $\Delta_hu(x) = \frac{u(x+h) - 2u(x) + u(x-h)}{h^2}$（**3 点差分**）、$d = 2$ では
$$
\Delta_hu(x,y) = \frac{u(x+h,y) + u(x-h,y) + u(x,y+h) + u(x,y-h) - 4u(x,y)}{h^2}
$$
（**5 点差分**）である。$\Omega_h$ のどの 2 点も、$\Omega_h$ の点だけを通って隣接点をたどる列で結べるとき、$\Omega_h$ は**連結**であるという。
</Definition>

以下では主に、$h = 1/(m+1)$（$m \ge 1$）として、単位区間 $\Omega_h = \{ih : 1 \le i \le m\}$ と単位正方形 $\Omega_h = \{(ih, jh) : 1 \le i, j \le m\}$ を扱います。単位正方形では $\Gamma_h$ は 4 辺の上の格子点（4 隅を除く）です。4 隅は内部の点の隣接点ではないので、5 点差分には現れません。

$\Delta_hu(P) \ge 0$ は「$u(P)$ が隣接点での値の平均以下」ということです。調和関数の平均値の性質の離散版で、最大値原理はここから出ます。

<Theorem id="thm-discrete-maximum-principle" title="離散最大値原理">
$u$ を格子領域 $\bar\Omega_h$ 上の実数値関数とし、すべての $P \in \Omega_h$ で $-\Delta_hu(P) \le 0$ とする。

1. $\max_{\bar\Omega_h}u = \max_{\Gamma_h}u$。
2. $\Omega_h$ が連結で、$u$ が $\Omega_h$ の点で最大値 $\max_{\bar\Omega_h}u$ をとるならば、$u$ は $\bar\Omega_h$ で定数である。

$-\Delta_hu \ge 0$ のときは、最大を最小に置き換えた主張が成り立つ。
</Theorem>

<Proof of="thm-discrete-maximum-principle">
$M = \max_{\bar\Omega_h}u$ とおきます（$\bar\Omega_h$ は有限集合なので最大値があります）。

**鍵になる事実.** $P \in \Omega_h$ で $u(P) = M$ ならば、すべての $Q \in N(P)$ で $u(Q) = M$ です。実際、$\Delta_hu(P) \ge 0$ から $\sum_{Q\in N(P)}(u(Q) - M) \ge 0$ ですが、各項は $u(Q) - M \le 0$ なので、すべての項が $0$ でなければなりません。

**1.** $u(P) = M$ となる $P \in \bar\Omega_h$ をとります。$P \in \Gamma_h$ なら示すことはありません。$P \in \Omega_h$ のとき、点列 $P_k = P + kh\boldsymbol{e}_1$（$k = 0, 1, 2, \ldots$）を考えます。$\Omega_h$ は有限集合なので、$P_k \notin \Omega_h$ となる最小の $k \ge 1$ があります。$P_0, \ldots, P_{k-1} \in \Omega_h$ で、鍵になる事実をくり返し使うと $u(P_1) = \cdots = u(P_k) = M$ です。$P_k$ は $P_{k-1} \in \Omega_h$ の隣接点で $\Omega_h$ に属さないので $P_k \in \Gamma_h$ です。したがって $\max_{\Gamma_h}u \ge M$ で、逆向きの不等式は自明です。

**2.** $u(P) = M$、$P \in \Omega_h$ とします。任意の $P' \in \Omega_h$ に対し、連結性から $\Omega_h$ の中の隣接点の列 $P = R_0, R_1, \ldots, R_s = P'$ があり、鍵になる事実を順に使うと $u(R_1) = \cdots = u(P') = M$ です。したがって $u$ は $\Omega_h$ で $M$ で、$\Gamma_h$ の各点は $\Omega_h$ のある点の隣接点なので、再び鍵になる事実から $\Gamma_h$ でも $M$ です。

最小についての主張は $-u$ に適用すれば得られます。
</Proof>

連続の場合の弱最大値原理（<Ref to="mathematics/pde/laplace-equation#cor-weak-maximum" />）では、[劣調和関数](/mathematics/pde/laplace-equation#def-harmonic) $-\Delta u \le 0$ の最大が境界でとられます。その証明は、劣調和関数の平均値の性質（球の中心での値は球の上の平均以下）から、内部の最大点のまわりで $u$ が定数になることを示すものでした。離散版では「$u(P)$ が隣接点での値の平均以下」がその役をし、上の証明の「鍵になる事実」が同じ議論になっています。

<Corollary id="cor-discrete-dirichlet-unique" title="離散 Dirichlet 問題の一意可解性">
格子領域 $\bar\Omega_h$ と、$\Omega_h$ 上の関数 $f$、$\Gamma_h$ 上の関数 $g$ に対し、$-\Delta_hu_h = f$（$\Omega_h$）、$u_h = g$（$\Gamma_h$）を満たす $u_h$ がただ 1 つ存在する。
</Corollary>

<Proof of="cor-discrete-dirichlet-unique">
$\Gamma_h$ での値を代入すると、未知数（$\Omega_h$ での値）と方程式がどちらも $\#\Omega_h$ 個の連立一次方程式です。$f = 0$、$g = 0$ のとき、<Ref to="thm-discrete-maximum-principle" /> の 1 を $u_h$ と $-u_h$ に使うと $\max u_h \le 0$、$\min u_h \ge 0$ なので $u_h = 0$ です。同次方程式の解が $0$ だけなので係数行列は正則で、任意の $f$、$g$ について解がただ 1 つあります。
</Proof>

これは連続の場合の一意性（<Ref to="mathematics/pde/laplace-equation#cor-dirichlet-uniqueness" />）に対応します。有限次元なので、一意性から存在も従います。

<Theorem id="thm-five-point-stability" title="最大値ノルムでの安定性">
$\bar\Omega_h$ を単位区間（$d = 1$）または単位正方形（$d = 2$）の格子領域（$h = 1/(m+1)$）とする。$\Gamma_h$ で $u = 0$ である格子関数 $u$ は
$$
\max_{\bar\Omega_h}|u| \le \frac18\max_{\Omega_h}|\Delta_hu|
$$
を満たす。
</Theorem>

<Proof of="thm-five-point-stability">
$d = 2$ の場合を書きます。比較関数
$$
w(x,y) = \frac14\Bigl(\bigl(x - \tfrac12\bigr)^2 + \bigl(y - \tfrac12\bigr)^2\Bigr)
$$
をとります。2 次関数 $p(x) = (x - \frac12)^2$ では $p(x+h) - 2p(x) + p(x-h) = 2h^2$ が正確に成り立つので、$\Delta_hw = \frac14(2 + 2) = 1$ です。また $[0,1]^2$ で $0 \le w \le \frac14(\frac14 + \frac14) = \frac18$ です。

$F = \max_{\Omega_h}|\Delta_hu|$ とし、$v = u + Fw$ とおくと、$\Delta_hv = \Delta_hu + F \ge 0$ です。<Ref to="thm-discrete-maximum-principle" /> の 1 から、$\Gamma_h$ で $u = 0$ であることを使って
$$
\max_{\bar\Omega_h}v = \max_{\Gamma_h}v = F\max_{\Gamma_h}w \le \frac{F}{8}
$$
です。$w \ge 0$ なので $u \le v \le F/8$ です。$-u$ も同じ仮定を満たすので $-u \le F/8$ で、合わせて $|u| \le F/8$ です。

$d = 1$ では $w(x) = \frac12(x - \frac12)^2$ とすれば、$\Delta_hw = 1$、$0 \le w \le \frac18$ で、同じ議論が成り立ちます。
</Proof>

この定理は、$-\Delta_h$ を行列とみたとき、最大値ノルムから誘導されるノルムで $\|(-\Delta_h)^{-1}\|_\infty \le 1/8$ であることを述べています。$h$ によらないことが要点で、これが定常問題での「安定性」です。

<Theorem id="thm-five-point-convergence" title="5 点差分の収束">
$\Omega$ を単位区間 $(0,1)$ または単位正方形 $(0,1)^2$、$u \in C^4(\bar\Omega)$ とし、$f = -\Delta u$、$g = u|_{\partial\Omega}$ とおく。$u_h$ を $-\Delta_hu_h = f$（$\Omega_h$）、$u_h = g$（$\Gamma_h$）の解とすると
$$
\max_{\bar\Omega_h}|u - u_h| \le \frac{h^2}{96}\Bigl(\max_{\bar\Omega}|\partial_x^4u| + \max_{\bar\Omega}|\partial_y^4u|\Bigr)
$$
が成り立つ（$d = 1$ では右辺は $\frac{h^2}{96}\max|u^{(4)}|$）。
</Theorem>

<Proof of="thm-five-point-convergence">
$d = 2$ で書きます。$M_x = \max|\partial_x^4u|$、$M_y = \max|\partial_y^4u|$ とします。

**整合性.** $P = (x,y) \in \Omega_h$ で、$y$ を固定して $x$ の関数とみた 2 階中心差分に <Ref to="computer-science/pde-discretization/finite-difference-convergence#lem-difference-truncation" /> の 3 を使うと、$\bigl|\frac{u(x+h,y) - 2u(x,y) + u(x-h,y)}{h^2} - u_{xx}(x,y)\bigr| \le \frac{h^2}{12}M_x$ で、$y$ 方向も同様です（この補題の証明は閉区間 $[x-h, x+h]$ での $C^4$ 級だけを使うので、境界に隣り合う点でも使えます）。したがって $\tau = \Delta_hu - \Delta u$ は $|\tau(P)| \le \frac{h^2}{12}(M_x + M_y)$ を満たします。

**誤差の方程式.** $e = u - u_h$ とおくと、$\Gamma_h$ で $e = g - g = 0$ で、$\Omega_h$ で
$$
\Delta_he = \Delta_hu - \Delta_hu_h = (\Delta u + \tau) + f = \tau
$$
です（$\Delta u = -f$）。

**安定性.** <Ref to="thm-five-point-stability" /> から $\max|e| \le \frac18\max|\tau| \le \frac{h^2}{96}(M_x + M_y)$ です。$d = 1$ も同じです。
</Proof>

証明の 3 段（整合性・誤差の方程式・安定性）は、前章の <Ref to="computer-science/pde-discretization/finite-difference-convergence#thm-lax-sufficiency" /> と同じ構造です。$u \in C^4$ の仮定は、右辺 $f$ が滑らかで領域に角がなければ満たされますが、正方形では角での整合条件が要ります。解の滑らかさが足りないと次数が下がることを §6 の数値実験で見ます。

<Remark id="rem-curved-boundary" title="一般の領域">
<Ref to="thm-discrete-maximum-principle" /> は任意の格子領域で成り立つので、比較関数さえ作れば安定性の評価は一般の有界領域に拡張できます（$\Omega$ が半径 $R$ の球に含まれるなら、$w = (R^2 - |x - x_0|^2)/(2d)$ で定数 $R^2/(2d)$）。曲がった境界では格子点が境界に乗らないので、境界の近くで不等間隔の差分（Shortley–Weller の近似）を使います。その場合も最大値原理による証明がそのまま通り、$O(h^2)$ の収束が得られます。Hackbusch（参考文献）に詳しい扱いがあります。
</Remark>

## 3. 離散 Laplace 作用素の固有値と ℓ² での安定性

この節の目標は、離散 Laplace 作用素の固有値と固有ベクトルを具体的に書き下し、最小の固有値が $h$ によらず連続の場合の値の近くにあることを確かめることです。次の命題の 3 の不等式の連鎖は、$\lambda_1 = \frac{4}{h^2}\sin^2\frac{\pi h}{2}$ を、$\sin x$ を $x$ とその 3 次の項で上下から挟んで $\pi^2$ に近いと示すための評価です。

境界で $0$ の格子関数を $\Omega_h$ 上の関数と同一視し、単位正方形（$d = 2$）では内積とノルムを
$$
\langle v, w\rangle_h = h^2\sum_{i,j=1}^{m}v_{ij}w_{ij},\qquad \|v\|_h = \langle v, v\rangle_h^{1/2}
$$
で定めます（$v_{ij} = v(ih, jh)$）。単位区間では $\langle v, w\rangle_h = h\sum_iv_iw_i$ です。$m^2h^2 < 1$ なので $\|v\|_h \le \|v\|_\infty$ です。

<Proposition id="prop-discrete-laplacian-spectrum" title="離散 Laplace 作用素のスペクトル">
$h = 1/(m+1)$、$\lambda_k = \dfrac{4}{h^2}\sin^2\dfrac{k\pi h}{2}$（$1 \le k \le m$）とおく。

1. 単位区間では、$\varphi_k(ih) = \sin(ik\pi h)$ が $-\Delta_h\varphi_k = \lambda_k\varphi_k$ を満たし、$\{\varphi_k\}_{k=1}^m$ は $\langle\cdot,\cdot\rangle_h$ について直交する基底である。
2. 単位正方形では、$\varphi_{kl}(ih, jh) = \sin(ik\pi h)\sin(jl\pi h)$（$1 \le k, l \le m$）が $-\Delta_h\varphi_{kl} = (\lambda_k + \lambda_l)\varphi_{kl}$ を満たし、$\{\varphi_{kl}\}$ は直交する基底である。したがって境界で $0$ の格子関数の空間の上で $-\Delta_h$ は対称で正定値であり、固有値は $\lambda_k + \lambda_l$ の全体である。
3. $4 \le \pi^2\bigl(1 - \frac{\pi^2h^2}{12}\bigr) \le \lambda_1 \le \pi^2$ である。したがって、境界で $0$ の格子関数 $v$ について、$d = 1, 2$ で $\|v\|_h \le \dfrac{1}{d\lambda_1}\|\Delta_hv\|_h \le \dfrac{1}{4d}\|\Delta_hv\|_h$ と $\dfrac{1}{d\lambda_1} \le \dfrac{1}{d\pi^2}\Bigl(1 - \dfrac{\pi^2h^2}{12}\Bigr)^{-1}$ が成り立つ。
4. <Ref to="thm-five-point-convergence" /> の設定で $\|u - u_h\|_h \le \frac{h^2}{48d}\bigl(\max|\partial_x^4u| + \max|\partial_y^4u|\bigr)$（$d = 1$ では括弧の中は $\max|u^{(4)}|$）。
</Proposition>

<Proof of="prop-discrete-laplacian-spectrum">
**1.** 境界で $0$ の格子関数 $v = (v_1, \ldots, v_m)$ に対し、$-h^2\Delta_hv$ は、対角成分が $2$、副対角成分が $-1$ の三重対角行列 $T$ を掛けることです（$v_0 = v_{m+1} = 0$）。<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#lem-tridiagonal-toeplitz-eigen" /> を $\alpha = 2$、$\beta = -1$、$n = m$、$\theta_k = k\pi/(m+1) = k\pi h$ で使うと、$T\varphi_k = (2 - 2\cos k\pi h)\varphi_k = 4\sin^2\frac{k\pi h}{2}\varphi_k$ で、固有値は相異なります（<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#cor-laplacian-condition" /> の形）。$T$ は対称なので、$k \ne k'$ のとき $\lambda_k\langle\varphi_k, \varphi_{k'}\rangle_h = \langle-\Delta_h\varphi_k, \varphi_{k'}\rangle_h = \langle\varphi_k, -\Delta_h\varphi_{k'}\rangle_h = \lambda_{k'}\langle\varphi_k, \varphi_{k'}\rangle_h$ と $\lambda_k \ne \lambda_{k'}$ から $\langle\varphi_k, \varphi_{k'}\rangle_h = 0$ です。$0$ でない互いに直交する $m$ 個のベクトルなので基底です。

**2.** 5 点差分は $x$ 方向の 3 点差分と $y$ 方向の 3 点差分の和です。$\varphi_{kl}$ は $x$ について $\varphi_k$ の定数倍、$y$ について $\varphi_l$ の定数倍で、辺の上では $\sin0 = \sin k\pi = 0$ から $0$ なので、1 から $-\Delta_h\varphi_{kl} = \lambda_k\varphi_{kl} + \lambda_l\varphi_{kl}$ です。内積は
$$
\langle\varphi_{kl}, \varphi_{k'l'}\rangle_h = \Bigl(h\sum_i\varphi_k(ih)\varphi_{k'}(ih)\Bigr)\Bigl(h\sum_j\varphi_l(jh)\varphi_{l'}(jh)\Bigr)
$$
と積に分かれ、1 から $(k,l) \ne (k',l')$ のとき $0$ です。$0$ でない互いに直交する $m^2$ 個のベクトルなので、$m^2$ 次元の空間の基底です。直交する固有ベクトルの基底をもち固有値が実数なので $-\Delta_h$ は対称で、固有値はすべて正なので正定値です。行列の言葉では、自然な順序で $-h^2\Delta_h = I \otimes T + T \otimes I$ です。

**3.** $x = \pi h/2 \in (0, \pi/4]$ とします。$\sin x \le x$ から $\lambda_1 \le \frac{4}{h^2}x^2 = \pi^2$ です。$\sin$ の 3 次までの Taylor 展開の剰余は $\frac{x^5}{120}\cos c \ge 0$（$0 < c < x \le \pi/4$）なので $\sin x \ge x - \frac{x^3}{6}$ で、
$$
\lambda_1 \ge \frac{4}{h^2}x^2\Bigl(1 - \frac{x^2}{6}\Bigr)^2 = \pi^2\Bigl(1 - \frac{\pi^2h^2}{24}\Bigr)^2 \ge \pi^2\Bigl(1 - \frac{\pi^2h^2}{12}\Bigr)
$$
です（$(1-s)^2 \ge 1 - 2s$）。$h \le 1/2$ なので右辺は $\pi^2(1 - \pi^2/48) > 9.86 \times 0.79 > 4$ です。$v = \sum c_{kl}\varphi_{kl}$ と展開すると、直交性から $\|\Delta_hv\|_h^2 = \sum(\lambda_k + \lambda_l)^2c_{kl}^2\|\varphi_{kl}\|_h^2 \ge (2\lambda_1)^2\|v\|_h^2$ で、$d = 2$ の主張を得ます。$d = 1$ も同じです。

**4.** <Ref to="thm-five-point-convergence" /> の証明のとおり $\Delta_he = \tau$、$e|_{\Gamma_h} = 0$、$\max|\tau| \le \frac{h^2}{12}(M_x + M_y)$ です。3 と $\|\tau\|_h \le \max|\tau|$ から $\|e\|_h \le \frac{1}{4d}\cdot\frac{h^2}{12}(M_x + M_y)$ です。
</Proof>

連続の問題では、単位正方形の Dirichlet ラプラシアンの最小固有値は $2\pi^2$（固有関数 $\sin\pi x\sin\pi y$）です。3 は、離散化しても最小固有値が $2\lambda_1 = 2\pi^2(1 + O(h^2))$ にとどまることを示しています。固有値の離散化誤差そのものは、このテーマの最後の章で扱います。

最大固有値は $2\lambda_m = \frac{8}{h^2}\cos^2\frac{\pi h}{2}$ なので、条件数は $\lambda_m/\lambda_1 = \cot^2\frac{\pi h}{2} \approx \frac{4}{\pi^2h^2}$ で、1 次元と同じです（<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#cor-laplacian-condition" />）。条件数が大きいほど、反復法は解に近づくのが遅くなります。このため、連立方程式 $-\Delta_hu_h = f$ を反復法で解くと、格子を細かくするほど収束が遅くなります。Jacobi 法と Gauss–Seidel 法では反復回数が $h^{-2}$ に、最適な SOR 法（Gauss–Seidel 法に緩和係数を入れて加速した反復法）では $h^{-1}$ に比例し（<Ref to="mathematics/numerical-linear-algebra/stationary-iterative-methods#prop-model-problem-rates" />）、共役勾配法では $\sqrt{\kappa} \approx \frac{2}{\pi h}$ に比例します（<Ref to="mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient#thm-cg-chebyshev-bound" />）。反復回数が $h$ によらない方法（多重格子法）は Hackbusch（参考文献）を見てください。

## 4. 熱方程式の θ 法

熱方程式の初期値境界値問題
$$
u_t = u_{xx}\quad(0 < x < 1,\ 0 < t \le T),\qquad u(0,t) = u(1,t) = 0,\qquad u(x,0) = u_0(x)
$$
を考えます。$h = 1/(m+1)$、$x_j = jh$ とし、$A$ を $-\Delta_h$（境界で $0$）の行列、すなわち対角成分が $2/h^2$、副対角成分が $-1/h^2$ の $m$ 次三重対角行列とします。

<Definition id="def-theta-method" title="θ 法">
$0 \le \theta \le 1$、$\Delta t > 0$ とする。$u^n \in \mathbb{R}^m$ から $u^{n+1}$ を
$$
\frac{u^{n+1} - u^n}{\Delta t} = -A\bigl(\theta u^{n+1} + (1-\theta)u^n\bigr),\qquad\text{すなわち}\qquad(I + \theta\Delta tA)u^{n+1} = \bigl(I - (1-\theta)\Delta tA\bigr)u^n
$$
で定める方法を **θ 法**という。$\theta = 0$ を**陽的 Euler 法**、$\theta = 1$ を**陰的 Euler 法**、$\theta = 1/2$ を **Crank–Nicolson 法**という。$\mu = \Delta t/h^2$ を**格子比**という。$B_\theta = I + \theta\Delta tA$、$E_\theta = I - (1-\theta)\Delta tA$、$C_\theta = B_\theta^{-1}E_\theta$ と書く。
</Definition>

$A$ は[正定値](/mathematics/calculus/multivariable-differentiation#def-definite)（<Ref to="prop-discrete-laplacian-spectrum" />）なので $B_\theta$ は正則です。成分で書くと、$j$ 番目の式は
$$
-\theta\mu u^{n+1}_{j-1} + (1 + 2\theta\mu)u^{n+1}_j - \theta\mu u^{n+1}_{j+1} = (1-\theta)\mu u^n_{j-1} + \bigl(1 - 2(1-\theta)\mu\bigr)u^n_j + (1-\theta)\mu u^n_{j+1}
$$
です（$u_0 = u_{m+1} = 0$）。$\theta > 0$ では毎ステップ三重対角の連立方程式を解きます。

次の定理のねらいは、1 ステップで各固有成分が $g_\theta(\sigma_k)$ 倍になる（これが増幅因子です）と見て、その絶対値が 1 を超えないための条件を調べることです。不安定を示す 3 の後半は、$|g_\theta|$ が 1 を超える範囲に最も細かい波 $k = m$ が入ってしまうことを、補助の数 $\sigma'$ を使って「$h$ を小さくすれば確実に入る」と言い直したものです。増幅因子の様子は次の図でも確かめられます。

<Figure caption="増幅因子 g(σ) の曲線。縦軸の 1 と -1 の線を超えない範囲が安定。陽的 Euler（θ = 0）は σ > 2 で -1 を下回る。Crank–Nicolson（θ = 1/2）は -1 に近づくだけで超えず、陰的 Euler（θ = 1）は 0 に近づく。">
<svg viewBox="0 0 500 270" width="100%" role="img" aria-label="横軸 σ、縦軸 g(σ) のグラフ。陽的 Euler は直線的に下がって σ = 2 で -1 を横切る。Crank–Nicolson は -1 に漸近する。陰的 Euler は 0 に漸近する。">
<g stroke="currentColor" fill="none" stroke-width="1">
<line x1="50" y1="30" x2="470" y2="30" stroke-dasharray="4 3" opacity="0.5"/>
<line x1="50" y1="190" x2="470" y2="190" stroke-dasharray="4 3" opacity="0.5"/>
<line x1="50" y1="110" x2="470" y2="110"/>
<line x1="50" y1="30" x2="50" y2="230"/>
</g>
<g fill="currentColor" font-size="12">
<text x="30" y="34">1</text>
<text x="30" y="114">0</text>
<text x="24" y="194">-1</text>
<text x="472" y="114">σ</text>
<text x="56" y="24">g(σ)</text>
<text x="148" y="126">2</text>
<text x="248" y="126">4</text>
<text x="348" y="126">6</text>
<text x="448" y="126">8</text>
</g>
<polyline fill="none" stroke="var(--sl-color-accent)" stroke-width="2" points="50,30 175,230"/>
<polyline fill="none" stroke="currentColor" stroke-width="2" points="50,30 62,47.8 75,62 88,73.6 100,83.3 125,98.6 150,110 175,118.9 200,126 250,136.7 300,144.3 350,150 400,154.4 450,158"/>
<polyline fill="none" stroke="currentColor" stroke-width="2" stroke-dasharray="6 3" points="50,30 62,46 75,56.7 88,64.3 100,70 125,78 150,83.3 175,87.1 200,90 250,94 300,96.7 350,98.6 400,100 450,101.1"/>
<g fill="currentColor" font-size="12">
<text x="180" y="226" fill="var(--sl-color-accent)">陽的 Euler</text>
<text x="300" y="170">Crank–Nicolson</text>
<text x="300" y="88">陰的 Euler</text>
</g>
</svg>
</Figure>

<Theorem id="thm-theta-stability" title="θ 法の ℓ² 安定性">
$0 \le \theta \le 1$ とし、$g_\theta(\sigma) = \dfrac{1 - (1-\theta)\sigma}{1 + \theta\sigma}$、$\sigma_k = \Delta t\lambda_k = 4\mu\sin^2\dfrac{k\pi h}{2}$ とおく。

1. $C_\theta$ の固有値は $g_\theta(\sigma_k)$（$1 \le k \le m$）で、$\|\cdot\|_h$ から誘導されるノルムについて $\|C_\theta^n\| = \max_k|g_\theta(\sigma_k)|^n$。
2. $\frac12 \le \theta \le 1$ ならば、すべての $\mu > 0$ と $h$ で $\|C_\theta\| \le 1$（無条件安定）。
3. $0 \le \theta < \frac12$ とする。$\mu \le \dfrac{1}{2(1-2\theta)}$ ならばすべての $h$ で $\|C_\theta\| \le 1$。$\mu > \dfrac{1}{2(1-2\theta)}$ ならば、$\mu$ を固定して $h \to 0$ としたとき $\sup\{\|C_\theta^n\| : n\Delta t \le T\} \to \infty$ で、安定でない。
</Theorem>

<Proof of="thm-theta-stability">
**1.** <Ref to="prop-discrete-laplacian-spectrum" /> の 1 から $A\varphi_k = \lambda_k\varphi_k$ なので、$B_\theta\varphi_k = (1 + \theta\sigma_k)\varphi_k$、$E_\theta\varphi_k = (1 - (1-\theta)\sigma_k)\varphi_k$ で、$C_\theta\varphi_k = g_\theta(\sigma_k)\varphi_k$ です。$\{\varphi_k\}$ は直交する基底なので、$v = \sum_kc_k\varphi_k$ について $\|C_\theta^nv\|_h^2 = \sum_k|g_\theta(\sigma_k)|^{2n}c_k^2\|\varphi_k\|_h^2 \le \max_k|g_\theta(\sigma_k)|^{2n}\|v\|_h^2$ で、最大を与える $k$ の $\varphi_k$ で等号になります。

**増幅因子の性質.** $\sigma > 0$ で $1 + \theta\sigma > 0$ です。$g_\theta(\sigma) \le 1$ は $1 - (1-\theta)\sigma \le 1 + \theta\sigma$、すなわち $-\sigma \le 0$ と同値で、常に成り立ちます。$g_\theta(\sigma) \ge -1$ は $1 - (1-\theta)\sigma \ge -1 - \theta\sigma$、すなわち
$$
(1 - 2\theta)\sigma \le 2
$$
と同値です。また $g_\theta'(\sigma) = \dfrac{-(1-\theta)(1+\theta\sigma) - \theta(1 - (1-\theta)\sigma)}{(1+\theta\sigma)^2} = \dfrac{-1}{(1+\theta\sigma)^2} < 0$ なので $g_\theta$ は減少関数です。

**2.** $\theta \ge \frac12$ なら $(1-2\theta)\sigma \le 0 \le 2$ なので、すべての $\sigma > 0$ で $|g_\theta(\sigma)| \le 1$ です。

**3.** $\theta < \frac12$ とし、$\sigma^{*} = \frac{2}{1-2\theta}$ とおくと、$|g_\theta(\sigma)| \le 1$ は $\sigma \le \sigma^{*}$ と同値です。$\sigma_k < 4\mu$ なので、$4\mu \le \sigma^{*}$、すなわち $\mu \le \frac{1}{2(1-2\theta)}$ ならば $\|C_\theta\| \le 1$ です。逆に $4\mu > \sigma^{*}$ とし、$\sigma^{*} < \sigma' < 4\mu$ となる $\sigma'$ をとって $q = -g_\theta(\sigma')$ とおくと、$g_\theta$ は減少で $g_\theta(\sigma^{*}) = -1$ なので $q > 1$ です。$\sigma_m = 4\mu\sin^2\frac{m\pi h}{2} = 4\mu\cos^2\frac{\pi h}{2} \to 4\mu$（$h \to 0$）なので、十分小さい $h$ で $\sigma_m \ge \sigma'$、したがって $g_\theta(\sigma_m) \le -q$ です。1 から $\|C_\theta^n\| \ge q^n$ で、$n = \lfloor T/\Delta t\rfloor = \lfloor T/(\mu h^2)\rfloor \to \infty$ なので有界ではありません。
</Proof>

不安定になるのは最も細かい波 $k = m$（格子点ごとに符号が変わる成分）です。陽的 Euler では $\mu \le 1/2$、すなわち $\Delta t \le h^2/2$ が条件です。全空間の格子 $h\mathbb{Z}$ では、前章の <Ref to="computer-science/pde-discretization/finite-difference-convergence#thm-von-neumann" /> の増幅因子が $g_\theta(4\mu\sin^2(\theta'/2))$（$\theta'$ は波数）になり、同じ条件が出ます。有界な区間では、Fourier 変換の代わりに固有ベクトル $\varphi_k$ での展開が同じ役をします。

<Theorem id="thm-theta-maximum" title="θ 法の離散最大値原理">
$0 \le \theta \le 1$、$\mu(1 - \theta) \le \frac12$ とする。θ 法の解は、すべての $j$ について
$$
\min\bigl(0, \min_iu^n_i\bigr) \le u^{n+1}_j \le \max\bigl(0, \max_iu^n_i\bigr)
$$
を満たす。特に $\|u^{n+1}\|_\infty \le \|u^n\|_\infty$ で、$u^n \ge 0$ なら $u^{n+1} \ge 0$ である。また、$\mu$ によらず、$B_\theta v = w$ ならば $\|v\|_\infty \le \|w\|_\infty$ である。
</Theorem>

<Proof of="thm-theta-maximum">
$M = \max(0, \max_iu^n_i)$ とおき、$w = E_\theta u^n$ とします。

**陽的な部分.** $w_j = (1-\theta)\mu u^n_{j-1} + (1 - 2(1-\theta)\mu)u^n_j + (1-\theta)\mu u^n_{j+1}$ の 3 つの係数は、仮定 $2(1-\theta)\mu \le 1$ から非負で、和は $1$ です。$u^n_i \le M$（境界の $u^n_0 = u^n_{m+1} = 0$ も $M \ge 0$ から $M$ 以下）なので $w_j \le M$ です。

**陰的な部分.** $B_\theta u^{n+1} = w$ の $j$ 番目の式は $(1 + 2\theta\mu)u^{n+1}_j = w_j + \theta\mu(u^{n+1}_{j-1} + u^{n+1}_{j+1})$ です。$M' = \max_iu^{n+1}_i$ とし、$M' \le 0$ なら $M' \le M$ です。$M' > 0$ のとき、$u^{n+1}_j = M'$ となる $j$ をとると、隣の値は（境界の $0$ も含めて）$M'$ 以下なので
$$
(1 + 2\theta\mu)M' \le w_j + 2\theta\mu M' \le M + 2\theta\mu M'
$$
で、$M' \le M$ です。どちらの場合も $u^{n+1}_j \le M$ です。下からの評価は $-u^n$、$-u^{n+1}$ に同じ議論を使います。$\|u^{n+1}\|_\infty \le \|u^n\|_\infty$ と非負性はここから従います。

**最後の主張.** $B_\theta v = w$ とし、$|v_j| = \|v\|_\infty$ となる $j$ をとると、$(1 + 2\theta\mu)|v_j| \le |w_j| + \theta\mu(|v_{j-1}| + |v_{j+1}|) \le \|w\|_\infty + 2\theta\mu|v_j|$ から $\|v\|_\infty \le \|w\|_\infty$ です。
</Proof>

これは連続の場合の放物型の最大値原理（<Ref to="mathematics/pde/heat-equation#thm-parabolic-maximum" />。最大は初期時刻か側面の境界でとられる）の離散版です。$B_\theta$ は非対角成分が $0$ 以下の正定値対称行列、つまり Stieltjes 行列（<Ref to="mathematics/numerical-linear-algebra/gmres-and-preconditioning#def-stieltjes" />）で、上の証明の「陰的な部分」は、その逆行列の成分が非負であることを直接示したものです。条件 $\mu(1-\theta) \le 1/2$ は、陽的 Euler で $\mu \le 1/2$（$\ell^2$ の安定条件と同じ）、Crank–Nicolson で $\mu \le 1$、陰的 Euler では無条件です。

<Example id="ex-cn-positivity" title="Crank–Nicolson の最大値原理は μ > 1 で破れうる">
$m = 1$（$h = 1/2$、未知数は中点の値 1 つ）では $A = 2/h^2$、$\Delta tA = 2\mu$ で、Crank–Nicolson は $u^{n+1} = \frac{1 - \mu}{1 + \mu}u^n$ です。$\mu > 1$ なら正の値が 1 ステップで負になり、<Ref to="thm-theta-maximum" /> の結論が成り立ちません。したがって、すべての格子で成り立つ条件としては $\mu \le 1$ を緩められません。一方、$|u^{n+1}| \le |u^n|$ は保たれています。格子が細かいときの様子を、$C_{1/2}$ の成分の最小値と行和のノルム（$\|C\|_\infty$。<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#prop-induced-norm-formulas" /> の 2）で調べます。

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

def cn_matrix(m, mu):                           # Crank–Nicolson の C = B^{-1} E
    T = 2 * np.eye(m) - np.eye(m, k=1) - np.eye(m, k=-1)
    return np.linalg.solve(np.eye(m) + 0.5 * mu * T, np.eye(m) - 0.5 * mu * T)

for m in (1, 2, 50):
    for mu in (1.0, 1.1, 1.2, 1.5, 1.6):
        C = cn_matrix(m, mu)
        print(f"m={m:2d} mu={mu:3.1f}  最小成分={C.min():8.4f}"
              f"  行和ノルム={np.abs(C).sum(axis=1).max():6.4f}")
```

```text
m= 1 mu=1.0  最小成分=  0.0000  行和ノルム=0.0000
m= 1 mu=1.1  最小成分= -0.0476  行和ノルム=0.0476
m= 1 mu=1.2  最小成分= -0.0909  行和ノルム=0.0909
m= 1 mu=1.5  最小成分= -0.2000  行和ノルム=0.2000
m= 1 mu=1.6  最小成分= -0.2308  行和ノルム=0.2308
m= 2 mu=1.0  最小成分=  0.0667  行和ノルム=0.3333
m= 2 mu=1.1  最小成分=  0.0225  行和ノルム=0.2903
m= 2 mu=1.2  最小成分= -0.0179  行和ノルム=0.2857
m= 2 mu=1.5  最小成分= -0.1209  行和ノルム=0.3846
m= 2 mu=1.6  最小成分= -0.1503  行和ノルム=0.4118
m=50 mu=1.0  最小成分=  0.0000  行和ノルム=1.0000
m=50 mu=1.1  最小成分=  0.0000  行和ノルム=1.0000
m=50 mu=1.2  最小成分= -0.0109  行和ノルム=1.0000
m=50 mu=1.5  最小成分= -0.1111  行和ノルム=1.0000
m=50 mu=1.6  最小成分= -0.1397  行和ノルム=1.0482
```

$m = 1$ の行は $\frac{1-\mu}{1+\mu}$（$\mu = 1.1$ で $-0.0476$）と一致します。$m = 2$ では $C$ の対角成分は $\frac12\bigl(g(\mu) + g(3\mu)\bigr)$（$g = g_{1/2}$。固有値 $\sigma = \mu, 3\mu$ から）で、これが $0$ になるのは $\mu = 2/\sqrt3 \approx 1.155$ です。表の符号の変わり目と合います。$m = 50$ でも $\mu = 1.2$ で負の成分が現れ、非負性は $\mu \le 1$ を少し超えたところで失われます（$m = 50$、$\mu = 1.1$ の最小成分は $1.3\times10^{-27}$ で正です）。つまり、$C_{1/2}$ の成分が非負である $\mu$ の範囲は格子点の数 $m$ によって変わり、$m = 1$ では $\mu \le 1$ ちょうど、$m = 2$ では $\mu \le 2/\sqrt3$ で、$m = 50$ でも $\mu = 1.1$ では非負です。<Ref to="thm-theta-maximum" /> が $\mu(1-\theta) \le \frac12$ を十分条件としてだけ述べているのはこのためです。一方、行和のノルムはこの表では $\mu = 1.5$ まで $1$ 以下です。Crank–Nicolson が最大値ノルムで縮小的であるための格子比の限界を調べた研究に Kraaijevanger の論文（参考文献）があります。ここでは、$\ell^2$ での無条件安定性が最大値原理を意味しないことを確認するにとどめます。
</Example>

次の定理のねらいは、熱方程式の厳密解を格子点で測った値 $w^n$ と θ 法の計算値 $u^n$ の差が、$\Delta t$ と $h$ のどんな冪で小さくなるかを調べることです。証明は、厳密解を差分式に代入したときに残るずれ（打ち切り誤差 $\tau$）を Taylor 展開で評価し、それを安定性で足し合わせる、という前章と同じ流れです。長い剰余項は、最終的に $\bar\tau$ の 3 つの項（$\Delta t$ の 1 次、$\Delta t^2$、$h^2$）にまとまります。

<Theorem id="thm-theta-convergence" title="θ 法の収束の次数">
$u$ を熱方程式の初期値境界値問題の解とし、$u$、$u_t$、$u_{tt}$、$u_{ttt}$、$\partial_x^ku$（$k \le 4$）と $u_{xxt}$、$u_{xxtt}$ が $[0,1]\times[0,T]$ で連続とする。$M_k = \max|\partial_t^ku|$ とおく。$u^0_j = u_0(x_j)$ から θ 法で計算した $u^n$ と $w^n_j = u(x_j, t_n)$ について、次が成り立つ（$t_n \le T$）。

1. $\ell^2$ の安定条件（$\theta \ge \frac12$、または $\mu \le \frac{1}{2(1-2\theta)}$）のもとで $\|w^n - u^n\|_h \le t_n\bar\tau$。ここで $\bar\tau = \bigl|\theta - \tfrac12\bigr|\Delta t\,M_2 + \dfrac{\Delta t^2}{6}M_3 + \dfrac{h^2}{12}M_2$。
2. $\mu(1-\theta) \le \frac12$ のもとでは、同じ評価が $\|w^n - u^n\|_\infty$ についても成り立つ。
3. $\theta = 0$、$\mu = \frac16$ で、さらに $\partial_x^6u$ が連続ならば、$\|w^n - u^n\|_\infty \le \dfrac{t_nh^4}{135}M_3$。
</Theorem>

<Proof of="thm-theta-convergence">
**打ち切り誤差.** $t^{*} = t_n + \Delta t/2$ とし、格子点 $x = x_j$ で
$$
\tau^n_j = \frac{u(x,t_{n+1}) - u(x,t_n)}{\Delta t} - \theta\,D_+D_-u(x,t_{n+1}) - (1-\theta)D_+D_-u(x,t_n)
$$
とおきます（$D_+D_-$ は $x$ についての 2 階中心差分）。<Ref to="computer-science/pde-discretization/finite-difference-convergence#lem-difference-truncation" /> の 2 を、$t^{*}$ を中心に刻み $\Delta t/2$ で使うと、第 1 項は $u_t(x,t^{*}) + R_a$、$|R_a| \le \frac{(\Delta t/2)^2}{6}M_3 = \frac{\Delta t^2}{24}M_3$ です。同じ補題の 3 から $D_+D_-u(x,t) = u_{xx}(x,t) + R_b(t)$、$|R_b| \le \frac{h^2}{12}\max|\partial_x^4u|$ です。$u_{xx}$ を $t$ について $t^{*}$ のまわりで 2 階の剰余まで展開すると（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-taylor" />）
$$
u_{xx}\bigl(x, t^{*} \pm \tfrac{\Delta t}{2}\bigr) = u_{xx}(x,t^{*}) \pm \frac{\Delta t}{2}u_{xxt}(x,t^{*}) + \frac{\Delta t^2}{8}u_{xxtt}(x,c_\pm)
$$
なので、$\theta u_{xx}(t_{n+1}) + (1-\theta)u_{xx}(t_n) = u_{xx}(t^{*}) + (\theta - \frac12)\Delta t\,u_{xxt}(t^{*}) + R_c$、$|R_c| \le \frac{\Delta t^2}{8}\max|u_{xxtt}|$ です。以上から
$$
\tau^n_j = \underbrace{u_t - u_{xx}}_{=0} - \bigl(\theta - \tfrac12\bigr)\Delta t\,u_{xxt}(x,t^{*}) + R_a - R_c - \theta R_b(t_{n+1}) - (1-\theta)R_b(t_n)
$$
です。$u_{xx} = u_t$ を $t$ で微分して $u_{xxt} = u_{tt}$、$u_{xxtt} = u_{ttt}$、$x$ で 2 回微分して $\partial_x^4u = \partial_x^2u_t = \partial_tu_{xx} = u_{tt}$ なので
$$
|\tau^n_j| \le \bigl|\theta - \tfrac12\bigr|\Delta t\,M_2 + \Bigl(\frac{1}{24} + \frac18\Bigr)\Delta t^2M_3 + \frac{h^2}{12}M_2 =: \bar\tau
$$
です。

**誤差の方程式.** $\tau^n$ の定義は $B_\theta w^{n+1} - E_\theta w^n = \Delta t\,\tau^n$ と書けます（境界条件 $u(0,t) = u(1,t) = 0$ から、$D_+D_-$ は $-A$ を掛けることに一致します）。したがって $w^{n+1} = C_\theta w^n + \Delta t\,B_\theta^{-1}\tau^n$ です。

**1.** $B_\theta^{-1}$ の固有値は $1/(1 + \theta\sigma_k) \in (0, 1]$ なので $\|B_\theta^{-1}\tau^n\|_h \le \|\tau^n\|_h \le \max_j|\tau^n_j| \le \bar\tau$ です。<Ref to="thm-theta-stability" /> から $\|C_\theta^n\| \le 1$ なので、<Ref to="computer-science/pde-discretization/finite-difference-convergence#lem-error-accumulation" /> を $K = 1$、$e^0 = 0$、局所誤差 $B_\theta^{-1}\tau^n$ で使って $\|w^n - u^n\|_h \le t_n\bar\tau$ です。

**2.** <Ref to="thm-theta-maximum" /> から、最大値ノルムで $\|C_\theta v\|_\infty \le \|v\|_\infty$、$\|B_\theta^{-1}\tau^n\|_\infty \le \|\tau^n\|_\infty \le \bar\tau$ です。同じ補題を最大値ノルムで使えば結論を得ます。

**3.** $\theta = 0$ のときは $t_n$ のまわりで展開します。$\frac{u(t_{n+1}) - u(t_n)}{\Delta t} = u_t + \frac{\Delta t}{2}u_{tt} + \frac{\Delta t^2}{6}u_{ttt}(x,c)$ です。6 階の剰余までの展開 $u(x\pm h) = \sum_{k=0}^{5}\frac{(\pm h)^k}{k!}\partial_x^ku + \frac{h^6}{720}\partial_x^6u(c_\pm)$ の和から $D_+D_-u = u_{xx} + \frac{h^2}{12}\partial_x^4u + R$、$|R| \le \frac{h^4}{360}\max|\partial_x^6u|$ です。$\partial_x^4u = u_{tt}$、$\partial_x^6u = u_{ttt}$ なので
$$
\tau^n_j = \Bigl(\frac{\Delta t}{2} - \frac{h^2}{12}\Bigr)u_{tt} + \frac{\Delta t^2}{6}u_{ttt}(x,c) - R
$$
で、$\Delta t = h^2/6$ なら第 1 項が消えて $|\tau^n_j| \le \bigl(\frac{1}{216} + \frac{1}{360}\bigr)h^4M_3 = \frac{h^4}{135}M_3$ です。$\mu = \frac16$ は $\mu(1-\theta) \le \frac12$ を満たすので、2 と同じ議論で結論を得ます。
</Proof>

Crank–Nicolson が 2 次なのは、$t^{*} = t_n + \Delta t/2$ について時間の差分が中心差分になり、空間の項も $t_n$ と $t_{n+1}$ の平均で $t^{*}$ での値を 2 次の精度で近似するからです。$\mu = 1/6$ の陽的 Euler では、時間と空間の誤差の先頭項が $u_{xxxx} = u_{tt}$ のおかげで打ち消し合います。熱方程式に特有の偶然で、低階の項や変数係数があると成り立ちません。

## 5. 高波数の減衰、エネルギー法、演算量

熱方程式の厳密解では、固有関数 $\sin k\pi x$ の成分は $\Delta t$ のあいだに $e^{-k^2\pi^2\Delta t}$ 倍になり、$k$ が大きいほど速く減衰します。初期値に角や跳びがあっても、$t > 0$ では解が滑らかになるのはこのためです（<Ref to="mathematics/pde/heat-equation#thm-heat-semigroup" /> の 4 は全空間での同じ性質です）。差分スキームがこの減衰を再現するかどうかは、$\sigma \to \infty$ での増幅因子 $g_\theta(\sigma)$ の振る舞いで決まります。

<Proposition id="prop-cn-damping" title="Crank–Nicolson と陰的 Euler の高波数での増幅因子">
1. Crank–Nicolson の増幅因子は $g_{1/2}(\sigma) = -1 + \dfrac{4}{2 + \sigma}$ で、$\sigma \to \infty$ で $-1$ に近づく。$\sigma \ge 2$ のとき $|g_{1/2}(\sigma)|^n \ge 1 - \dfrac{4n}{2 + \sigma}$。特に $h \le \frac12$ のとき、最も細かい波 $k = m$ について $|g_{1/2}(\sigma_m)|^n \ge 1 - \dfrac{2n}{1 + \mu}$。
2. 陰的 Euler の増幅因子は $g_1(\sigma) = \dfrac{1}{1 + \sigma} \in (0, 1)$ で、$\sigma \to \infty$ で $0$ に近づく。
3. 最初の $n_0$ ステップを陰的 Euler、その後を Crank–Nicolson とする方法では、固有成分 $\varphi_k$ は $n$ ステップ後に $(1 + \sigma_k)^{-n_0}g_{1/2}(\sigma_k)^{n - n_0}$ 倍になる。<Ref to="thm-theta-convergence" /> の仮定のもとで、誤差は
$$
\|w^n - u^n\|_h \le n_0\Delta t\Bigl(\frac{\Delta t}{2}M_2 + \frac{\Delta t^2}{6}M_3 + \frac{h^2}{12}M_2\Bigr) + t_n\Bigl(\frac{\Delta t^2}{6}M_3 + \frac{h^2}{12}M_2\Bigr)
$$
を満たし、$n_0$ を固定すれば $O(\Delta t^2 + h^2)$ のままである。
</Proposition>

<Proof of="prop-cn-damping">
**1.** $g_{1/2}(\sigma) = \frac{1 - \sigma/2}{1 + \sigma/2} = \frac{2 - \sigma}{2 + \sigma} = -1 + \frac{4}{2+\sigma}$ です。$\sigma \ge 2$ なら $|g_{1/2}(\sigma)| = 1 - \frac{4}{2+\sigma} \in [0, 1)$ で、Bernoulli の不等式 $(1 - s)^n \ge 1 - ns$（$0 \le s \le 1$）から $|g_{1/2}(\sigma)|^n \ge 1 - \frac{4n}{2+\sigma}$ です。$h \le \frac12$ なら $\sigma_m = 4\mu\cos^2\frac{\pi h}{2} \ge 4\mu\cos^2\frac{\pi}{4} = 2\mu$ です。$s \mapsto 1 - \frac{4n}{2+s}$ は増加なので、$\sigma_m \ge 2$ のとき $|g_{1/2}(\sigma_m)|^n \ge 1 - \frac{4n}{2 + 2\mu} = 1 - \frac{2n}{1+\mu}$ です。$\sigma_m < 2$ のときは $\mu < 1$ で右辺は $n \ge 1$ について負なので、不等式は自明です。

**2.** $g_1(\sigma) = \frac{1 - 0}{1 + \sigma}$ です。

**3.** <Ref to="thm-theta-stability" /> の 1 のとおり $C_1\varphi_k = g_1(\sigma_k)\varphi_k$、$C_{1/2}\varphi_k = g_{1/2}(\sigma_k)\varphi_k$ なので、倍率は積です。誤差について、ステップ $k$ で使う θ を $\theta_k$（$k < n_0$ で $1$、それ以外は $\frac12$）とすると、<Ref to="thm-theta-convergence" /> の証明から $e^{k+1} = C_{\theta_k}e^k + \Delta t\,B_{\theta_k}^{-1}\tau^k$、$\|C_{\theta_k}\| \le 1$、$\|B_{\theta_k}^{-1}\| \le 1$ なので $\|e^{k+1}\|_h \le \|e^k\|_h + \Delta t\|\tau^k\|_h$ です。$e^0 = 0$ から $\|e^n\|_h \le \Delta t\sum_{k<n}\|\tau^k\|_h$ で、$\tau^k$ の評価（$\theta = 1$ で $|\theta - \frac12| = \frac12$、$\theta = \frac12$ で $0$）を代入すれば主張を得ます。第 1 項は $n_0\Delta t\cdot O(\Delta t + h^2) = O(\Delta t^2 + \Delta t\,h^2)$ です。
</Proof>

Crank–Nicolson を $\Delta t = ch$ のように $h$ に比例させて使うと $\mu = c/h$ は大きく、1 から、最も細かい成分が半分になるまでに少なくとも $(1+\mu)/4$ ステップかかります（$1 - \frac{2n}{1+\mu} \le \frac12$ となるのは $n \ge \frac{1+\mu}{4}$ のときです）。しかも符号が毎ステップ反転します。滑らかな初期値ではこの成分がもともと小さいので問題になりませんが、跳びのある初期値では振動として見えます（§6）。陰的 Euler は高波数を 1 ステップで $1/(1+\sigma)$ 倍にするので、最初の 2 ステップ程度を陰的 Euler にすれば、3 のとおり 2 次精度を保ったまま高波数を落とせます。常微分方程式の数値解法の言葉では、陰的 Euler は L 安定、Crank–Nicolson は A 安定だが L 安定でない、と表現されます（L 安定は、A 安定に加えて、非常に速く減衰する成分を 1 ステップで 0 の近くまで落とせる性質です）。

固有値が閉じた形で分からない問題でも、次のエネルギー法なら安定性を示せます。

<Proposition id="prop-energy-estimate" title="エネルギー法による安定性と変数係数への拡張">
$a$ は $[0,1]$ 上の関数で $0 < a_0 \le a(x) \le a_1$ とし、$a_{j+1/2} = a((j + \frac12)h)$（$0 \le j \le m$）とおく。境界で $0$（$v_0 = v_{m+1} = 0$）の格子関数に作用する
$$
(A_av)_j = -\frac{1}{h^2}\Bigl(a_{j+1/2}(v_{j+1} - v_j) - a_{j-1/2}(v_j - v_{j-1})\Bigr)\qquad(1 \le j \le m)
$$
（$-(au_x)_x$ の差分近似）について、次が成り立つ。

1. （部分和分）$\langle A_av, w\rangle_h = h\displaystyle\sum_{j=0}^{m}a_{j+1/2}\frac{v_{j+1} - v_j}{h}\cdot\frac{w_{j+1} - w_j}{h}$。したがって $A_a$ は対称で正定値。
2. $A_a$ の固有値は $\Lambda = 4a_1/h^2$ 以下で、$\|A_aw\|_h^2 \le \Lambda\langle A_aw, w\rangle_h$。
3. θ 法 $\frac{v^{n+1} - v^n}{\Delta t} = -A_a(\theta v^{n+1} + (1-\theta)v^n)$ の解は、$w = \theta v^{n+1} + (1-\theta)v^n$ として
$$
\|v^{n+1}\|_h^2 - \|v^n\|_h^2 = -2\Delta t\langle A_aw, w\rangle_h - (2\theta - 1)\Delta t^2\|A_aw\|_h^2
$$
を満たす。$\theta \ge \frac12$、または $\theta < \frac12$ かつ $a_1\mu \le \frac{1}{2(1-2\theta)}$ ならば $\|v^{n+1}\|_h \le \|v^n\|_h$。
</Proposition>

<Proof of="prop-energy-estimate">
**1.** $p_j = a_{j+1/2}(v_{j+1} - v_j)$（$0 \le j \le m$）とおくと $(A_av)_j = -(p_j - p_{j-1})/h^2$ です。$w_0 = w_{m+1} = 0$ を使って
$$
\sum_{j=1}^{m}(p_j - p_{j-1})w_j = \sum_{j=0}^{m}p_jw_j - \sum_{j=0}^{m}p_jw_{j+1} = -\sum_{j=0}^{m}p_j(w_{j+1} - w_j)
$$
です（第 1 の和に $p_0w_0 = 0$ を、第 2 の和に $p_mw_{m+1} = 0$ を足し、第 2 の和の添字を 1 つずらしました）。したがって $\langle A_av, w\rangle_h = h\sum_{j=1}^m(A_av)_jw_j = \frac1h\sum_{j=0}^{m}p_j(w_{j+1} - w_j)$ で、主張の式です。右辺は $v$ と $w$ について対称で、$w = v$ のとき $a_0h\sum_j\bigl(\frac{v_{j+1} - v_j}{h}\bigr)^2$ 以上です。これが $0$ なら $v_{j+1} = v_j$（すべての $j$）で、$v_0 = 0$ から $v = 0$ です。したがって正定値です。

**2.** $A_a$ の第 $j$ 行は、対角成分が $(a_{j+1/2} + a_{j-1/2})/h^2$、非対角成分の絶対値の和がそれ以下です。<Ref to="mathematics/linear-algebra/matrix-inequalities#thm-gershgorin" /> から、固有値は $2(a_{j+1/2} + a_{j-1/2})/h^2 \le 4a_1/h^2 = \Lambda$ 以下です。$A_a$ は実対称なので、<Ref to="mathematics/linear-algebra/spectral-theorem#cor-real-symmetric" /> から固有ベクトルからなる正規直交基底 $\{q_k\}$（固有値 $0 < \kappa_k \le \Lambda$）があります。$w = \sum_kc_kq_k$ とすると、$\|A_aw\|^2 = \sum_k\kappa_k^2c_k^2 \le \Lambda\sum_k\kappa_kc_k^2 = \Lambda\langle A_aw, w\rangle$（Euclid 内積で。両辺に $h$ を掛ければ $\langle\cdot,\cdot\rangle_h$ での式）です。

**3.** $I + \theta\Delta tA_a$ は正定値なので $v^{n+1}$ は一意に定まります。$w = \frac12(v^{n+1} + v^n) + (\theta - \frac12)(v^{n+1} - v^n)$ と書けるので
$$
\langle v^{n+1} - v^n, w\rangle_h = \tfrac12\bigl(\|v^{n+1}\|_h^2 - \|v^n\|_h^2\bigr) + \bigl(\theta - \tfrac12\bigr)\|v^{n+1} - v^n\|_h^2
$$
です。スキームから $v^{n+1} - v^n = -\Delta tA_aw$ なので、左辺は $-\Delta t\langle A_aw, w\rangle_h$、右辺の最後の項は $(\theta - \frac12)\Delta t^2\|A_aw\|_h^2$ で、2 倍して整理すると主張の等式です。$\theta \ge \frac12$ なら右辺は $0$ 以下です。$\theta < \frac12$ のとき、2 から
$$
-(2\theta - 1)\Delta t^2\|A_aw\|_h^2 \le (1 - 2\theta)\Delta t^2\Lambda\langle A_aw, w\rangle_h
$$
なので、$(1-2\theta)\Delta t\Lambda \le 2$、すなわち $a_1\mu \le \frac{1}{2(1-2\theta)}$ ならば右辺全体が $0$ 以下です。
</Proof>

$a = 1$ のとき $A_a = A$ で、3 は <Ref to="thm-theta-stability" /> の十分条件を固有値の式なしで再現します。$\langle A_av, v\rangle_h$ は連続の問題のエネルギー $\int a|u_x|^2dx$ の離散版で、$\theta \ge 1/2$ では「$\|v^n\|_h^2$ が減る速さは散逸 $2\Delta t\langle A_aw, w\rangle_h$ 以上」という、熱方程式のエネルギー等式 $\frac{d}{dt}\|u\|^2 = -2\int|u_x|^2dx$ の離散版になっています。部分和分は部分積分の離散版です。

<Proposition id="prop-explicit-stiffness" title="陽的な方法と陰的な方法の演算量">
熱方程式を時刻 $T$ まで解くとする。

1. 陽的 Euler 法は、安定であるために $\Delta t \le h^2/2$ が必要で、ステップ数は $2T/h^2$ 以上である。1 ステップは $m$ 個の成分の更新なので、演算量は $h^{-3}$ に比例する。
2. $\theta > 0$ の θ 法の 1 ステップは、$m$ 次の正定値対称な三重対角行列 $B_\theta$ を係数とする連立方程式で、Cholesky 分解（最初に 1 回、演算量 $6m$ 以下）と前進・後退代入（演算量は $m$ に比例）で解ける。演算量は $m\cdot T/\Delta t$、すなわち $h^{-1}\Delta t^{-1}$ に比例する。
3. 誤差を $O(h^2)$ にするのに、陽的 Euler は $h^{-3}$ に比例する演算量を要し、Crank–Nicolson は $\Delta t$ を $h$ に比例させて $h^{-2}$ に比例する演算量で足りる。
</Proposition>

<Proof of="prop-explicit-stiffness">
**1.** <Ref to="thm-theta-stability" /> の 3（$\theta = 0$）から、$\mu > \frac12$ を固定すると安定でないので $\Delta t \le h^2/2$ が必要です。ステップ数は $T/\Delta t \ge 2T/h^2$、1 ステップは各 $j$ で定数回の演算（$m = 1/h - 1$ 個）です。

**2.** $B_\theta = I + \theta\Delta tA$ は対称で、固有値 $1 + \theta\sigma_k > 0$ なので正定値、帯幅は上下とも $1$ です。<Ref to="mathematics/numerical-linear-algebra/cholesky-and-sparse-direct#thm-band-preserved" /> の 2 を $p = 1$、$n = m$ で使うと、Cholesky 因子も上帯幅 $1$ で、分解の演算量は $m(1+1)(2+1) = 6m$ 以下です。$B_\theta$ は $n$ によらないので分解は 1 回で済み、各ステップは右辺 $E_\theta u^n$ の計算（$m$ に比例）と、帯幅 $1$ の三角行列 2 つによる代入（各行で定数回の演算）です。

**3.** <Ref to="thm-theta-convergence" /> から、陽的 Euler の誤差は $\Delta t \le h^2/2$ のもとで $O(\Delta t + h^2) = O(h^2)$、Crank–Nicolson の誤差は $O(\Delta t^2 + h^2)$ で、$\Delta t = ch$ とすれば $O(h^2)$ です。1 と 2 に代入します。
</Proof>

安定性のために必要な時間刻みが、精度のために必要な時間刻みよりずっと小さい問題を**硬い**問題といいます。拡散を含む方程式は、格子を細かくするほど硬くなります。2 次元では $B_\theta$ が三重対角でなくなり、1 ステップが 5 点差分の連立方程式になります。これを 1 次元の三重対角の方程式の列に分ける方法が交互方向陰解法（ADI）で、Strikwerda（参考文献）に解説があります。

## 6. 数値実験

まず、単位正方形の Poisson 方程式を 5 点差分で解きます。滑らかな解 $u = \sin\pi x\sin\pi y$ と、中心からの距離 $r$ の関数 $u = r^{3/2}$（$f = -\frac94r^{-1/2}$ は中心で非有界で、$u$ は $C^2$ 級でさえない）を比べます。$m$ を偶数にして、中心 $(\frac12, \frac12)$ が格子点に乗らないようにしています。

```python
# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1（arm64）
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve

def poisson(m, u, f):
    """単位正方形で -Δu = f、境界値は u。内部 m×m 点の 5 点差分。"""
    h = 1.0 / (m + 1)
    x = h * np.arange(1, m + 1)
    X, Y = np.meshgrid(x, x, indexing="ij")
    T = sp.diags([-1, 2, -1], [-1, 0, 1], shape=(m, m))
    I = sp.identity(m)
    A = (sp.kron(I, T) + sp.kron(T, I)).tocsr() / h**2
    b = f(X, Y)
    b[0, :] += u(0.0, x) / h**2                 # 境界値を右辺へ
    b[-1, :] += u(1.0, x) / h**2
    b[:, 0] += u(x, 0.0) / h**2
    b[:, -1] += u(x, 1.0) / h**2
    uh = spsolve(A, b.ravel()).reshape(m, m)
    return np.max(np.abs(uh - u(X, Y)))

pi = np.pi
r = lambda x, y: np.sqrt((x - 0.5)**2 + (y - 0.5)**2)
u1 = lambda x, y: np.sin(pi * x) * np.sin(pi * y)
f1 = lambda x, y: 2 * pi**2 * u1(x, y)
u2 = lambda x, y: r(x, y)**1.5                  # Δ r^α = α² r^(α-2)
f2 = lambda x, y: -2.25 * r(x, y)**(-0.5)
for name, u, f in [("sin*sin", u1, f1), ("r^1.5", u2, f2)]:
    e = [poisson(m, u, f) for m in (16, 32, 64, 128, 256)]
    o = [np.log2(e[i] / e[i + 1]) for i in range(4)]
    print(f"{name:8s}" + " ".join(f"{v:9.2e}" for v in e)
          + "  |" + " ".join(f"{v:5.2f}" for v in o))
print("上界 (h^2/96) * 2π^4（m = 16）:", f"{(1 / 17)**2 / 96 * 2 * pi**4:.2e}")
```

```text
sin*sin  2.83e-03  7.54e-04  1.95e-04  4.94e-05  1.25e-05  | 1.91  1.95  1.98  1.99
r^1.5    3.12e-03  1.54e-03  7.07e-04  3.09e-04  1.30e-04  | 1.02  1.12  1.19  1.24
上界 (h^2/96) * 2π^4（m = 16）: 7.02e-03
```

**検算.** (a) 参照解は製造解（$u$ を決めて $f = -\Delta u$ を計算したもの）です。2 次元で $\Delta r^\alpha = \alpha^2r^{\alpha-2}$ であることは、極座標の $\Delta = \partial_r^2 + r^{-1}\partial_r$ から $\alpha(\alpha - 1) + \alpha = \alpha^2$ です。(b) 滑らかな解では[観測次数](/computer-science/research-computing/verification-of-computations#def-observed-order)（$h$ を半分にしたときの誤差の比 $e(h)/e(h/2)$ の $\log_2$。表の `|` の右の列）が $2$ に近づき、$m = 16$ の誤差 $2.83\times10^{-3}$ は <Ref to="thm-five-point-convergence" /> の上界 $7.02\times10^{-3}$（$\partial_x^4u$ と $\partial_y^4u$ の最大はどちらも $\pi^4$）の内側です。(c) 別の方法として、$\sin\pi x\sin\pi y$ は離散 Laplace 作用素の固有ベクトル $\varphi_{11}$ なので、離散解は閉じた形 $u_h = \frac{2\pi^2}{2\lambda_1}\varphi_{11}$ で求まり、誤差の最大は $(\frac{\pi^2}{\lambda_1} - 1)\max\varphi_{11}$ です。$m = 16$（$h = 1/17$）では $\lambda_1 = 4\cdot17^2\sin^2(\pi/34) = 9.8415$ で $\frac{\pi^2}{\lambda_1} - 1 = 2.85\times10^{-3}$、$m$ が偶数のとき格子点での $\varphi_{11}$ の最大は $\cos^2(\pi h/2) = 0.9915$ なので、積は $2.83\times10^{-3}$ で表と一致します。(d) $u = r^{3/2}$ は <Ref to="thm-five-point-convergence" /> の仮定 $u \in C^4$ を満たさず、観測次数は $1.0$〜$1.2$ にとどまります。

次に、熱方程式を初期値 $u_0 = \sin\pi x + \frac12\sin3\pi x$（厳密解は $e^{-\pi^2t}\sin\pi x + \frac12e^{-9\pi^2t}\sin3\pi x$）で $T = 0.1$ まで解き、$\ell^2$ 誤差を表にします。`FE` は陽的 Euler、`BE` は陰的 Euler、`CN` は Crank–Nicolson です。

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

def theta_step_matrix(m, dt, theta):
    h = 1.0 / (m + 1)
    A = (2 * np.eye(m) - np.eye(m, k=1) - np.eye(m, k=-1)) / h**2
    B = np.eye(m) + theta * dt * A
    E = np.eye(m) - (1 - theta) * dt * A
    return np.linalg.solve(B, E)                # C = B^{-1} E

def heat_error(m, dt, theta, T=0.1):
    h = 1.0 / (m + 1)
    x = h * np.arange(1, m + 1)
    n = round(T / dt)
    dt = T / n                                  # T をちょうど n 等分する
    u = np.sin(np.pi * x) + 0.5 * np.sin(3 * np.pi * x)
    C = theta_step_matrix(m, dt, theta)
    for _ in range(n):
        u = C @ u
    ex = (np.exp(-np.pi**2 * T) * np.sin(np.pi * x)
          + 0.5 * np.exp(-9 * np.pi**2 * T) * np.sin(3 * np.pi * x))
    return np.sqrt(h * np.sum((u - ex)**2))

def table(label, theta, dt_of_h, ms=(9, 19, 39, 79, 159)):   # h = 1/10, ..., 1/160
    e = [heat_error(m, dt_of_h(1.0 / (m + 1)), theta) for m in ms]
    o = [np.log2(e[i] / e[i + 1]) for i in range(len(ms) - 1)]
    print(f"{label:11s}" + " ".join(f"{v:9.2e}" for v in e)
          + "  |" + " ".join(f"{v:7.2f}" for v in o))

table("FE mu=0.4", 0.0, lambda h: 0.4 * h**2)      # FE: 陽的 Euler
table("BE mu=0.4", 1.0, lambda h: 0.4 * h**2)      # BE: 陰的 Euler
table("CN mu=0.4", 0.5, lambda h: 0.4 * h**2)
table("FE mu=1/6", 0.0, lambda h: h**2 / 6)
table("BE dt=h/4", 1.0, lambda h: h / 4)
table("CN dt=h/4", 0.5, lambda h: h / 4)
table("FE mu=0.55", 0.0, lambda h: 0.55 * h**2, ms=(9, 19, 39, 79))
```

```text
FE mu=0.4   3.04e-03  7.51e-04  1.87e-04  4.68e-05  1.17e-05  |   2.01    2.00    2.00    2.00
BE mu=0.4   7.15e-03  1.81e-03  4.54e-04  1.14e-04  2.84e-05  |   1.98    2.00    2.00    2.00
CN mu=0.4   2.11e-03  5.33e-04  1.34e-04  3.34e-05  8.36e-06  |   1.98    2.00    2.00    2.00
FE mu=1/6   4.78e-06  2.97e-07  1.85e-08  1.16e-09  7.25e-11  |   4.01    4.00    4.00    4.00
BE dt=h/4   3.13e-02  1.58e-02  7.95e-03  3.99e-03  2.00e-03  |   0.99    0.99    0.99    1.00
CN dt=h/4   8.36e-04  2.08e-04  5.20e-05  1.30e-05  3.25e-06  |   2.01    2.00    2.00    2.00
FE mu=0.55  5.10e-03  1.23e-03  1.48e+06  1.41e+74  |   2.05  -30.17 -225.82
```

**検算.** (a) 参照解は固有関数展開の閉じた形です。(b) $\mu$ を固定すると $\Delta t = \mu h^2$ なので、$O(\Delta t + h^2)$ も $O(\Delta t^2 + h^2)$ も $O(h^2)$ で、最初の 3 行の次数 $2$ はこれに合います。$\Delta t = h/4$ では陰的 Euler が $1$ 次、Crank–Nicolson が $2$ 次で、<Ref to="thm-theta-convergence" /> の 1 の $|\theta - \frac12|\Delta t$ の項の有無に対応します。$\mu = 1/6$ の陽的 Euler は次数 $4$ で、同じ定理の 3 と合います。(c) $\mu = 0.55 > \frac12$ の陽的 Euler は、$h = 1/40$ 以下で発散します。最も細かい成分の増幅因子は $g_0(\sigma_m) = 1 - 2.2\cos^2(\pi h/2)$ で、$h = 1/80$ では $|g| = 1.1992$、ステップ数は $n = 1164$ なので $|g|^n \approx 10^{92}$ です。初期値のこの成分は丸め誤差の大きさ（$10^{-17}$ 程度）なので、$10^{75}$ 程度まで育つ計算になり、表の $1.4\times10^{74}$ と桁が合います。$h = 1/10$、$1/20$ で発散が見えないのは、ステップ数が少なく（$n = 18$、$73$）、丸め誤差がまだ育ちきっていないからで、安定なわけではありません。

<Example id="ex-cn-oscillation" title="階段状の初期値での Crank–Nicolson の振動">
初期値を $(\frac14, \frac34)$ で $1$、それ以外で $0$ とし、$h = 1/200$、$\Delta t = h$（$\mu = 200$）で解きます。厳密解は Fourier 正弦級数 $\sum_kb_ke^{-k^2\pi^2t}\sin k\pi x$、$b_k = \frac{2}{k\pi}\bigl(\cos\frac{k\pi}{4} - \cos\frac{3k\pi}{4}\bigr)$ です。$t > 0$ では係数が指数的に小さくなるので級数は一様に絶対収束します。4000 項で打ち切って参照解とします。解のぎざぎざの大きさを、2 階差分の絶対値の最大 $\max_j|u_{j+1} - 2u_j + u_{j-1}|$ で測ります。

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

m = 199
h = 1.0 / (m + 1)
x = h * np.arange(1, m + 1)
A = (2 * np.eye(m) - np.eye(m, k=1) - np.eye(m, k=-1)) / h**2
I = np.eye(m)
u0 = np.where((x > 0.25) & (x < 0.75), 1.0, 0.0)        # 階段状の初期値

def exact(t, K=4000):                                    # Fourier 正弦級数
    k = np.arange(1, K + 1)[:, None]
    b = 2 / (k * np.pi) * (np.cos(k * np.pi / 4) - np.cos(3 * k * np.pi / 4))
    return np.sum(b * np.exp(-(k * np.pi)**2 * t) * np.sin(k * np.pi * x), axis=0)

def run(dt, n, n_be):
    BE = np.linalg.solve(I + dt * A, I)
    CN = np.linalg.solve(I + 0.5 * dt * A, I - 0.5 * dt * A)
    u = u0.copy()
    for k in range(n):
        u = (BE if k < n_be else CN) @ u
    return u

wiggle = lambda u: np.max(np.abs(u[2:] - 2 * u[1:-1] + u[:-2]))   # 2 階差分の最大
dt = h                                                   # mu = 1/h = 200
print(f"mu = {dt / h**2:.0f}")
for n in (4, 16, 64):
    ex = exact(n * dt)
    for name, n_be in [("CN", 0), ("BE x2 + CN", 2), ("BE", n)]:
        u = run(dt, n, n_be)
        print(f"n={n:3d} {name:11s} max err={np.max(np.abs(u - ex)):8.2e}"
              f"  wiggle={wiggle(u):8.2e}  min={u.min():9.2e}")
    print(f"      exact       wiggle={wiggle(ex):8.2e}")
```

```text
mu = 200
n=  4 CN          max err=4.25e-01  wiggle=9.64e-01  min= 9.06e-03
n=  4 BE x2 + CN  max err=1.17e-02  wiggle=3.09e-04  min= 8.57e-03
n=  4 BE          max err=2.49e-02  wiggle=2.87e-04  min= 7.82e-03
      exact       wiggle=2.87e-04
n= 16 CN          max err=3.47e-01  wiggle=8.74e-01  min=-5.80e-02
n= 16 BE x2 + CN  max err=2.29e-03  wiggle=1.01e-04  min= 6.37e-03
n= 16 BE          max err=5.12e-03  wiggle=1.04e-04  min= 6.45e-03
      exact       wiggle=1.01e-04
n= 64 CN          max err=2.12e-01  wiggle=6.16e-01  min=-1.85e-01
n= 64 BE x2 + CN  max err=2.33e-04  wiggle=9.38e-06  min= 5.97e-04
n= 64 BE          max err=2.68e-03  wiggle=1.01e-05  min= 6.43e-04
      exact       wiggle=9.44e-06
```

Crank–Nicolson だけの解は、64 ステップ後（$t = 0.32$。厳密解は最大でも $0.04$ 程度まで減衰している時刻）でも 2 階差分が $0.6$ もあり、負の値（$-0.185$）をとります。<Ref to="prop-cn-damping" /> の 1 の下界は $1 - \frac{2\cdot64}{201} = 0.36$ で、最も細かい成分が 64 ステップ後にも 3 割以上残ることを保証しています。$\mu = 200 > 1$ なので <Ref to="thm-theta-maximum" /> の仮定も満たされず、負の値が出ることと矛盾しません。最初の 2 ステップを陰的 Euler にすると、2 階差分は厳密解と同じ大きさになり、誤差は陰的 Euler だけの場合の約 $1/10$（$n = 64$）です。陰的 Euler だけの解は常に滑らかで非負ですが、1 次精度なので誤差が大きくなります。
</Example>

<Example id="ex-explicit-averaging" title="陽的 Euler の μ = 1/2 は隣の平均">
$\theta = 0$、$\mu = \frac12$ のとき、スキームは $u^{n+1}_j = \frac12(u^n_{j-1} + u^n_{j+1})$ です。新しい値が両隣の平均なので最大値原理は一目で分かります。初期値が 1 点だけ $1$ の格子関数なら、$n$ ステップ後の値は二項係数 $\binom{n}{k}2^{-n}$ で、偶数番目と奇数番目の格子点が交互に $0$ になる市松模様が残ります。これは $g_0(\sigma_m) = 1 - 4\mu\cos^2(\pi h/2) \approx -1$ で、最も細かい成分が減衰しないことの現れです。$\mu$ を $\frac12$ より少しでも大きくすると、この成分が $|g| > 1$ で増幅されます（上の表の最後の行）。実用上は $\mu \le \frac12$ の内側、たとえば $\mu = \frac16$ や $\frac14$ を使います。
</Example>

## 7. 演習

<Exercise id="exr-sharp-constant" difficulty="易">
単位区間で $m$ が奇数のとき、<Ref to="thm-five-point-stability" /> の定数 $\frac18$ は改良できないことを示せ。
<Solution>
$u(x) = \frac12x(1-x)$ を格子点に制限した格子関数を考えます。2 次関数では 2 階中心差分が正確なので $\Delta_hu = u'' = -1$、また $u(0) = u(1) = 0$ です。$m$ が奇数なら $x = \frac12 = \frac{m+1}{2}h$ は格子点で、$u(\frac12) = \frac18$ です。したがって $\max|u| = \frac18 = \frac18\max|\Delta_hu|$ で、等号が成り立ちます。
</Solution>
</Exercise>

<Exercise id="exr-nonnegative-rhs" difficulty="標準">
単位正方形の格子領域で、$-\Delta_hu_h = f$、$\Gamma_h$ で $u_h = 0$ とする。$f \ge 0$ ならば $u_h \ge 0$ であること、さらに $f \le f'$ ならば対応する解が $u_h \le u_h'$ を満たすことを示せ。これは $(-\Delta_h)^{-1}$ の成分について何を意味するか。
<Solution>
$-\Delta_hu_h = f \ge 0$ なので、<Ref to="thm-discrete-maximum-principle" /> の最小についての主張から $\min_{\bar\Omega_h}u_h = \min_{\Gamma_h}u_h = 0$、すなわち $u_h \ge 0$ です。$f \le f'$ のときは、差 $u_h' - u_h$ が右辺 $f' - f \ge 0$、境界値 $0$ の解なので、同じ理由で非負です。

$f$ として 1 点 $Q$ だけで $1$、他で $0$ の格子関数をとると、$u_h$ は $(-\Delta_h)^{-1}$ の $Q$ に対応する列です。これが非負なので、$(-\Delta_h)^{-1}$ の成分はすべて $0$ 以上です。$-\Delta_h$ は非対角成分が $0$ 以下の正定値対称行列（Stieltjes 行列。<Ref to="mathematics/numerical-linear-algebra/gmres-and-preconditioning#def-stieltjes" />）で、これはその一般的な性質の特別な場合です。
</Solution>
</Exercise>

<Exercise id="exr-max-norm-be" difficulty="標準">
陰的 Euler 法（$\theta = 1$）について、<Ref to="thm-theta-convergence" /> の仮定のもとで、$\mu$ に何の条件も付けずに $\|w^n - u^n\|_\infty \le t_n\bigl(\frac{\Delta t}{2}M_2 + \frac{\Delta t^2}{6}M_3 + \frac{h^2}{12}M_2\bigr)$ が成り立つことを確かめよ。Crank–Nicolson で同じことを最大値ノルムで主張するには、この章の結果では何が要るか。
<Solution>
$\theta = 1$ では $\mu(1-\theta) = 0 \le \frac12$ が常に成り立つので、<Ref to="thm-theta-convergence" /> の 2 がそのまま使えます。$|\theta - \frac12| = \frac12$ を代入すれば主張の式です。

Crank–Nicolson では $\mu(1-\theta) \le \frac12$ は $\mu \le 1$、すなわち $\Delta t \le h^2$ を意味します。この章の結果で最大値ノルムの評価を得るにはこの制限が要り、$\Delta t$ を $h$ に比例させる使い方には適用できません。その場合に証明されているのは $\ell^2$ ノルムでの評価（同じ定理の 1）です。<Ref to="ex-cn-positivity" /> の表は、$\mu > 1$ でも行和のノルムが $1$ 以下にとどまる範囲があることを示していますが、それはこの章では証明していません。
</Solution>
</Exercise>

<Exercise id="exr-energy-decay" difficulty="難">
<Ref to="prop-energy-estimate" /> の設定で、離散 Poincaré 不等式 $\|v\|_h^2 \le \frac12h\sum_{j=0}^m\bigl(\frac{v_{j+1} - v_j}{h}\bigr)^2$（境界で $0$ の $v$）を示し、陰的 Euler 法の解が $\|v^n\|_h \le (1 + 4a_0\Delta t)^{-n/2}\|v^0\|_h$ を満たすことを導け。
<Solution>
$d_i = v_{i+1} - v_i$ とおくと、$v_0 = 0$ から $v_j = \sum_{i=0}^{j-1}d_i$ です。Cauchy–Schwarz の不等式から $v_j^2 \le j\sum_{i=0}^{m}d_i^2$ で、
$$
\|v\|_h^2 = h\sum_{j=1}^{m}v_j^2 \le h\cdot\frac{m(m+1)}{2}\sum_{i=0}^{m}d_i^2 \le \frac{1}{2h}\sum_{i=0}^md_i^2 = \frac12h\sum_{i=0}^{m}\Bigl(\frac{d_i}{h}\Bigr)^2
$$
です（$m(m+1) \le (m+1)^2 = h^{-2}$）。<Ref to="prop-energy-estimate" /> の 1 と合わせて $\langle A_av, v\rangle_h \ge a_0h\sum_i(d_i/h)^2 \ge 2a_0\|v\|_h^2$ です。

陰的 Euler（$\theta = 1$）では $w = v^{n+1}$ で、同じ命題の 3 の等式から
$$
\|v^{n+1}\|_h^2 - \|v^n\|_h^2 = -2\Delta t\langle A_av^{n+1}, v^{n+1}\rangle_h - \Delta t^2\|A_av^{n+1}\|_h^2 \le -4a_0\Delta t\|v^{n+1}\|_h^2
$$
です。したがって $(1 + 4a_0\Delta t)\|v^{n+1}\|_h^2 \le \|v^n\|_h^2$ で、くり返して結論を得ます。連続の問題での減衰 $\|u(t)\| \le e^{-\pi^2a_0t}\|u_0\|$ に対応する評価で、定数が $\pi^2$ でなく $2$ なのは、上の Poincaré 不等式の定数が最良でないためです。
</Solution>
</Exercise>

## 参考文献

- K. W. Morton, D. F. Mayers『Numerical Solution of Partial Differential Equations』2nd ed., Cambridge University Press, 2005 — 放物型方程式の θ 法、最大値原理による誤差解析、楕円型方程式の差分法。
- J. C. Strikwerda『Finite Difference Schemes and Partial Differential Equations』2nd ed., SIAM, 2004 — 放物型方程式のスキームの安定性、Crank–Nicolson の高波数の振る舞い、ADI 法。
- W. Hackbusch『Elliptic Differential Equations: Theory and Numerical Treatment』Springer, 1992 — Poisson 方程式の差分法、離散最大値原理、一般の領域、多重格子法への導入。
- R. J. LeVeque『Finite Difference Methods for Ordinary and Partial Differential Equations』SIAM, 2007 — 定常問題の安定性の定式化、拡散方程式、硬い問題と L 安定性。
- J. Crank, P. Nicolson, "A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type", *Proceedings of the Cambridge Philosophical Society* 43 (1947), 50–67 — Crank–Nicolson 法の原論文。
- J. F. B. M. Kraaijevanger, "Maximum norm contractivity of discretization schemes for the heat equation", *Applied Numerical Mathematics* 9 (1992) — 熱方程式の離散化スキームが最大値ノルムで縮小的になる条件。