# 波動方程式と Schrödinger 方程式の時間発展：leapfrog・Crank–Nicolson・分割演算子法

> 保存量をもつ発展方程式の時間離散化として、leapfrog 法の安定条件・離散エネルギーの保存・数値分散、Crank–Nicolson 法のユニタリ性と誤差、Lie–Trotter・Strang 分割の交換子による誤差評価、箱による切断の誤差を証明する。
> https://rikai.mugen-giken.com/computer-science/pde-discretization/wave-and-schrodinger-time-stepping

## 0. この記事の要点

波の動きや量子力学の粒子の動きをコンピュータで追うときは、時間を小さな刻みに区切って少しずつ進めます。熱の拡散と違って、これらの現象は時間がたっても衰えず、全体の大きさ（エネルギーや確率の総和）が一定に保たれます。ところが刻み方を誤ると、計算の上では大きさが勝手に膨らんだり、波の進み方がずれたりします。この記事を読み終えると、どんな刻み幅なら計算が壊れないか、大きさを正確に保つ刻み方は何か、そして誤差がどれだけになるかを説明できるようになります。

- 波動方程式を空間について離散化した系 $u'' = -Au$（$A$ は[正定値](/mathematics/linear-algebra/spectral-theorem#def-positive-definite)対称行列）に対する leapfrog 法は、$\Delta t^2\lambda_{\max}(A) \lt 4$ のとき、かつそのときに限り解が有界にとどまります。1 次元の[中心差分](/computer-science/pde-discretization/finite-difference-convergence#def-difference-quotients)では CFL 条件（波が 1 ステップで格子 1 マス分より遠くへ進まないための刻み幅の上限）$c\Delta t/h \le 1$ です。
- leapfrog 法は離散エネルギー $E^n = \frac12\|(u^{n+1}-u^n)/\Delta t\|^2 + \frac12\langle u^{n+1}, Au^n\rangle$ を厳密に保存します。安定条件のもとで $E^n$ は正定値で、これから誤差評価 $O(\Delta t^2 + h^2)$ が出ます。振幅は保たれますが、位相速度（波の山が進む速さ）は $-(1-\nu^2)(kh)^2/24$ の割合でずれ、長時間では位相の誤差が蓄積します。
- Schrödinger 方程式を離散化した系 $i\psi' = H\psi$ に対する Crank–Nicolson 法（Cayley 変換）は[ユニタリ](/mathematics/linear-algebra/spectral-theorem#def-hermitian-unitary)（ベクトルの長さを変えない変換）で、ノルムとエネルギーを厳密に保存し、誤差は $\frac{t\,\Delta t^2}{12}\|H^3\psi_0\|$ 以下です。陽的 Euler 法はどんな刻み幅でもノルムを増やします。
- $H = A + B$ を分けて $e^{-i\tau A}$ と $e^{-i\tau B}$ を交互に掛ける分割法は、各因子がユニタリなので積もユニタリです。誤差は交換子（$A$ と $B$ を掛ける順序を入れ替えたときのずれ $AB - BA$）で決まり、Lie–Trotter 分割で $\frac{t\tau}{2}\|[A,B]\|$、Strang 分割で $t\tau^2\bigl(\frac1{12}\|[A,[A,B]]\| + \frac1{24}\|[B,[B,A]]\|\bigr)$ 以下です。
- 全空間の問題を有限の箱で計算したときの誤差は、箱の外に出た波動関数の大きさを時間について積分したもので抑えられます。

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

- **ユニタリ行列**: ベクトルの長さと内積を変えない行列（$U^{*}U = I$ を満たす）。→ [ユニタリ行列](/mathematics/linear-algebra/spectral-theorem#def-hermitian-unitary)
- **エルミート行列**: 転置して複素共役をとっても変わらない行列。実対称行列の複素版で、固有値は実数。→ [エルミート行列](/mathematics/linear-algebra/spectral-theorem#def-hermitian-unitary)
- **正定値**: 零でないどのベクトル $x$ についても $\langle x, Ax\rangle \gt 0$ となる対称行列。→ [正定値](/mathematics/linear-algebra/spectral-theorem#def-positive-definite)
- **行列の指数関数**: $e^{tA}$ を冪級数で定めたもので、$u' = Au$ の解を与える。→ [行列指数関数](/mathematics/ode/linear-systems#def-matrix-exp)
- **中心差分**: 点 $x$ の前後の格子点の値の差で微分を近似する方法。→ [中心差分](/computer-science/pde-discretization/finite-difference-convergence#def-difference-quotients)

</Aside>

## 1. 動機

熱方程式の解は時間とともに滑らかになり、高い波数の成分はすぐに減衰します。そのため、数値解法が高い波数の成分を少し間違えても、誤差は自然に消えていきます。波動方程式と Schrödinger 方程式はそうではありません。波動方程式はエネルギー $\frac12\int(u_t^2 + c^2u_x^2)dx$ を保存し（<Ref to="mathematics/pde/wave-equation#thm-energy-conservation" />）、Schrödinger 方程式は確率の総和 $\int|\psi|^2dx$ を保存します。どの波数の成分も減衰せず、位相を変えながら永久に残ります。

このことは数値解法に 2 つの要求を課します。1 つは、保存量を壊さないことです。Schrödinger 方程式に陽的 Euler 法を使うと、$\|\psi\|$ は毎ステップ必ず増え、確率の総和が 1 を超えていきます（§6）。もう 1 つは、位相を正しく進めることです。振幅が正確でも、位相速度が波数によって少しずつずれると、波束は時間とともに形を崩します（§5）。

この章では、空間の離散化（差分法・スペクトル法）は済んでいるものとして、行列の常微分方程式 $u'' = -Au$ と $i\psi' = H\psi$ の時間積分を扱います。$A$ と $H$ は対称・エルミートなので、固有分解（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）によって、すべての主張はスカラーの漸化式と不等式に帰着します。分割法の誤差評価は、虚時間の場合を含む一般の形を[経路積分 Monte Carlo 法](/computer-science/stochastic-computation/path-integral-monte-carlo)の章で証明し、ここではユニタリな場合の系として使います。

## 2. 設定

<Definition id="def-semidiscrete-systems" title="空間を離散化した波動方程式と Schrödinger 方程式">
$A$ を $d$ 次の実対称な正定値行列、$H$ を $d$ 次のエルミート行列とする。

$$
u''(t) = -Au(t),\quad u(0) = g,\ u'(0) = g_1\qquad\text{と}\qquad i\psi'(t) = H\psi(t),\quad \psi(0) = \psi_0
$$

を、それぞれ**半離散の波動方程式**、**半離散の Schrödinger 方程式**という。$A$ の固有値を $0 \lt \lambda_1 \le \cdots \le \lambda_d$ と書く。$\langle x, y\rangle = \sum_j\bar x_jy_j$、$\|x\| = \langle x, x\rangle^{1/2}$ は $\mathbb{C}^d$ の標準内積とノルム、行列のノルム $\|\cdot\|$ はそれから誘導される作用素ノルム（2 ノルム）とする。交換子を $[X, Y] = XY - YX$ と書く。
</Definition>

1 次元の波動方程式 $u_{tt} = c^2u_{xx}$（$0 \lt x \lt 1$、両端で 0）を $h = 1/M$ の格子で中心差分により離散化すると、$A$ は $(M-1)$ 次の三重対角行列 $\frac{c^2}{h^2}\mathrm{tridiag}(-1, 2, -1)$ で、固有値は

$$
\lambda_k = \frac{4c^2}{h^2}\sin^2\frac{k\pi h}{2}\qquad(k = 1, \ldots, M-1)
$$

です（<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#lem-tridiagonal-toeplitz-eigen" /> で $\alpha = 2c^2/h^2$、$\beta = -c^2/h^2$）。Schrödinger 方程式 $i\psi_t = -\frac12\psi_{xx} + V\psi$ では $H = \frac{1}{2h^2}\mathrm{tridiag}(-1,2,-1) + \mathrm{diag}(V(x_j))$ です。格子関数については、$h$ の重みをつけたノルム $\|v\|_h = \sqrt h\|v\|$ と内積 $\langle v, w\rangle_h = h\langle v, w\rangle$ も使います。

<Lemma id="lem-unitary-group" title="エルミート行列の生成するユニタリ群">
$H$ をエルミート行列、$H = U\Lambda U^{*}$（$U$ はユニタリ、$\Lambda = \mathrm{diag}(\mu_1, \ldots, \mu_d)$ は実対角）をその固有分解とする。実数 $t$ について $e^{-itH} = U\,\mathrm{diag}(e^{-it\mu_k})\,U^{*}$ で、これはユニタリ行列であり、$H$ と可換である。半離散の Schrödinger 方程式の解は $\psi(t) = e^{-itH}\psi_0$ で、$\|\psi(t)\| = \|\psi_0\|$、$\langle\psi(t), H\psi(t)\rangle = \langle\psi_0, H\psi_0\rangle$ を満たす。また、$\|X\| \le 1$、$\|Y\| \le 1$ である行列について $\|X^n - Y^n\| \le n\|X - Y\|$（$n \ge 1$）である。
</Lemma>

<Proof of="lem-unitary-group">
固有分解は <Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" /> です。<Ref to="mathematics/ode/linear-systems#prop-exp-properties" /> の 4 から $U^{*}e^{-itH}U = e^{-itU^{*}HU} = e^{-it\Lambda}$ で、対角行列の指数関数は定義の級数（<Ref to="mathematics/ode/linear-systems#def-matrix-exp" />）から $\mathrm{diag}(e^{-it\mu_k})$ です。$|e^{-it\mu_k}| = 1$ なので $e^{-itH}$ はユニタリです（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-normal" />）。$H = U\Lambda U^{*}$ と $e^{-itH}$ は同じ $U$ で対角化されるので可換です。同じ命題の 3 から $\frac{d}{dt}e^{-itH}\psi_0 = -iHe^{-itH}\psi_0$ で、解の一意性（<Ref to="mathematics/ode/linear-systems#thm-variation-of-constants" />）から $\psi(t) = e^{-itH}\psi_0$ です。ユニタリ行列はノルムと内積を保つので（<Ref to="mathematics/linear-algebra/spectral-theorem#prop-unitary" />）、$\|\psi(t)\| = \|\psi_0\|$、$\langle\psi(t), H\psi(t)\rangle = \langle e^{-itH}\psi_0, e^{-itH}H\psi_0\rangle = \langle\psi_0, H\psi_0\rangle$ です。

最後の主張は $X^n - Y^n = \sum_{j=0}^{n-1}X^{n-1-j}(X - Y)Y^j$ と、各項のノルムが $\|X\|^{n-1-j}\|X - Y\|\|Y\|^j \le \|X - Y\|$ であることから出ます。
</Proof>

## 3. leapfrog 法と安定条件

<Definition id="def-leapfrog" title="leapfrog 法">
半離散の波動方程式に対し、$t_n = n\Delta t$ での近似 $u^n$ を

$$
u^{n+1} - 2u^n + u^{n-1} = -\Delta t^2Au^n\qquad(n \ge 1)
$$

で定める方法を **leapfrog 法**（Störmer–Verlet 法）という。出発値 $u^0$、$u^1$ は別に与える。すべての出発値について $\sup_n\|u^n\| \lt \infty$ であるとき、（その $A$ と $\Delta t$ について）**安定**であるという。
</Definition>

左辺を $\Delta t^2$ で割ったものは $u''(t_n)$ の 2 階中心差分です。3 つの時刻の値を使う 2 段のスキームなので、[差分法の整合性・安定性・収束](/computer-science/pde-discretization/finite-difference-convergence)の 1 段スキームの枠組みはそのままでは使えませんが、固有分解でスカラーの 3 項漸化式に帰着させれば直接調べられます。

次の補題は、各固有成分が従う 3 項漸化式の解を、$\mu$ の大きさで場合分けして書き下すものです。$\mu \lt 4$ なら解は振動して有界にとどまり、$\mu \ge 4$ なら増大する、という分かれ目がそのまま安定条件になります。証明は解の形を与えて漸化式と出発値を満たすことを確かめるだけですが、$\cos\varphi = 1 - \mu/2$ という形は、$z^n$ の形の解を探したときの特性方程式 $z^2 - (2-\mu)z + 1 = 0$ の根が、$\mu \lt 4$ では絶対値 1 の複素数 $e^{\pm i\varphi}$ になることから来ています。

<Lemma id="lem-three-term-recurrence" title="3 項漸化式の解">
$\mu \gt 0$ とし、複素数列 $(y^n)_{n\ge0}$ が $y^{n+1} - (2 - \mu)y^n + y^{n-1} = 0$（$n \ge 1$）を満たすとする。

1. $0 \lt \mu \lt 4$ のとき、$\cos\varphi = 1 - \mu/2$ で定まる $\varphi \in (0, \pi)$ について
   $$
   y^n = y^0\cos n\varphi + \frac{y^1 - y^0\cos\varphi}{\sin\varphi}\sin n\varphi,\qquad |y^n| \le |y^0| + \frac{|y^1 - y^0\cos\varphi|}{\sin\varphi}.
   $$
2. $\mu = 4$ のとき $y^n = (-1)^n\bigl(y^0 - n(y^0 + y^1)\bigr)$。$y^1 \ne -y^0$ なら $|y^n| \to \infty$（1 次の増大）。
3. $\mu \gt 4$ のとき、$z = \frac12\bigl(2 - \mu - \sqrt{\mu^2 - 4\mu}\bigr)$ は $z \lt -1$ を満たす実数で、$y^n = z^n$ は解である（幾何級数的な増大）。
</Lemma>

<Proof of="lem-three-term-recurrence">
漸化式は $y^0$、$y^1$ から $y^n$ を一意に決めるので、各場合に、右辺の式が漸化式を満たし、$n = 0, 1$ で $y^0$、$y^1$ に一致することを確かめれば十分です。

1. 加法定理から $\cos(n+1)\varphi + \cos(n-1)\varphi = 2\cos\varphi\cos n\varphi$、$\sin(n+1)\varphi + \sin(n-1)\varphi = 2\cos\varphi\sin n\varphi$ で、$2\cos\varphi = 2 - \mu$ なので、$\cos n\varphi$ と $\sin n\varphi$ はどちらも漸化式を満たします。$n = 0$ で右辺は $y^0$、$n = 1$ で $y^0\cos\varphi + (y^1 - y^0\cos\varphi) = y^1$ です。不等式は $|\cos n\varphi| \le 1$、$|\sin n\varphi| \le 1$ から出ます。

2. $y^n = (-1)^n(a + bn)$ を $y^{n+1} + 2y^n + y^{n-1}$ に入れると $(-1)^{n+1}\bigl[(a + b(n+1)) - 2(a + bn) + (a + b(n-1))\bigr] = 0$ です。$a = y^0$、$b = -(y^0 + y^1)$ とすると $n = 1$ で $-(y^0 - y^0 - y^1) = y^1$ です。

3. $z$ は $z^2 - (2-\mu)z + 1 = 0$ の解なので $z^n$ は漸化式を満たします。$\mu \gt 4$ なら $\mu^2 - 4\mu \gt 0$ で $z$ は実数、$2 - \mu \lt -2$ から $z \lt \frac12(2 - \mu) \lt -1$ です。
</Proof>

<Theorem id="thm-leapfrog-stability" title="leapfrog 法の安定条件">
<Ref to="def-leapfrog" /> の leapfrog 法が安定であるためには、$\Delta t^2\lambda_d \lt 4$ が必要十分である。$\Delta t^2\lambda_d = 4$ のときは $\|u^n\|$ が $n$ に比例して増える解が、$\Delta t^2\lambda_d \gt 4$ のときは幾何級数的に増える解が存在する。
</Theorem>

<Proof of="thm-leapfrog-stability">
$A$ の固有ベクトルからなる正規直交基底 $q_1, \ldots, q_d$（$Aq_k = \lambda_kq_k$）をとり（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）、$y_k^n = \langle q_k, u^n\rangle$ とおきます。leapfrog 法の式と $q_k$ の内積をとると、$\langle q_k, Au^n\rangle = \langle Aq_k, u^n\rangle = \lambda_ky_k^n$ から

$$
y_k^{n+1} - (2 - \mu_k)\,y_k^n + y_k^{n-1} = 0,\qquad \mu_k = \Delta t^2\lambda_k \gt 0
$$

です。$\|u^n\|^2 = \sum_k|y_k^n|^2$ です。

$\Delta t^2\lambda_d \lt 4$ なら、すべての $k$ で $0 \lt \mu_k \lt 4$ なので、<Ref to="lem-three-term-recurrence" /> の 1 から各 $|y_k^n|$ は $n$ について有界で、$\|u^n\|$ も有界です。$\Delta t^2\lambda_d = 4$ なら、$u^0 = q_d$、$u^1 = 0$ とすると $y_k^n = 0$（$k \ne d$）、同じ補題の 2 から $y_d^n = (-1)^n(1 - n)$ で、$\|u^n\| = |n - 1|$ です。$\Delta t^2\lambda_d \gt 4$ なら、$u^0 = q_d$、$u^1 = zq_d$（$z$ は同じ補題の 3 のもの）とすると $u^n = z^nq_d$、$\|u^n\| = |z|^n$ です。
</Proof>

<Corollary id="cor-wave-cfl" title="波動方程式の中心差分と CFL 条件">
1 次元の波動方程式 $u_{tt} = c^2u_{xx}$ の中心差分と leapfrog 法について、$\nu = c\Delta t/h$ とおく。

1. 区間 $(0,1)$ で両端 0（$h = 1/M$）のとき、安定であるための必要十分条件は $\nu\sin\bigl((M-1)\pi/(2M)\bigr) \lt 1$ である。したがって $\nu \le 1$ ならすべての $M$ で安定であり、$\nu \gt 1$ なら十分大きいすべての $M$ で不安定である。
2. 周期境界で格子点の数 $N$ が偶数のとき、$A$ は固有値 0（定数の格子関数）をもつ半正定値行列である。定数の成分を除いた部分（成分の和が 0 である出発値 $u^0$、$u^1$）について、$\sup_n\|u^n\| \lt \infty$ がつねに成り立つための必要十分条件は $\nu \lt 1$ である。$\nu = 1$ のとき $u_j^n = (-1)^{j+n}n$ は解で、1 次で増大する。
</Corollary>

<Proof of="cor-wave-cfl">
1. §2 の固有値の式から $\Delta t^2\lambda_{M-1} = 4\nu^2\sin^2((M-1)\pi/(2M))$ で、<Ref to="thm-leapfrog-stability" /> の条件は $\nu\sin((M-1)\pi/(2M)) \lt 1$ です。$\sin \lt 1$ なので $\nu \le 1$ なら成り立ちます。$M \to \infty$ で $\sin((M-1)\pi/(2M)) \to 1$ なので、$\nu \gt 1$ なら十分大きい $M$ で破れます。

2. 周期境界の行列 $A$ は半正定値で、$v = ((-1)^j)_j$（$N$ が偶数なので周期的）について $(Av)_j = -\frac{c^2}{h^2}(v_{j+1} - 2v_j + v_{j-1}) = \frac{4c^2}{h^2}v_j$ です。一般に $v_j = e^{2\pi ijm/N}$ は固有値 $\frac{4c^2}{h^2}\sin^2(\pi m/N)$ の固有ベクトルで、$m = 0, \ldots, N-1$ の $N$ 個は互いに直交して $\mathbb{C}^N$ の基底をなすので（$m \ne m'$ のとき $\sum_{j=0}^{N-1}e^{2\pi ij(m-m')/N} = 0$。等比数列の和の公式による）、最大固有値は $m = N/2$ の $4c^2/h^2$ です。$m = 0$ の固有値 0 の成分は $y^{n+1} - 2y^n + y^{n-1} = 0$ を満たし、出発値が $y^1 = y^0$ でなければ 1 次で増えますが、これは方程式 $u_{tt} = c^2u_{xx}$ の解 $u = t$ に対応する正しい振る舞いなので、安定性は $m \ne 0$ の成分について考えます。その条件は <Ref to="lem-three-term-recurrence" /> から $4\nu^2 \lt 4$ です。$\nu = 1$ のとき、$u_j^n = (-1)^{j+n}n$ を代入すると、左辺は $(-1)^{j+n}\bigl[-(n+1) - 2n - (n-1)\bigr] = -4n(-1)^{j+n}$、右辺は $\nu^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n) = (-1)^{j+n}n(-1 - 2 - 1) = -4n(-1)^{j+n}$ で一致します。
</Proof>

$\nu \le 1$ は、数値的な依存領域が波動方程式の依存領域（<Ref to="mathematics/pde/wave-equation#thm-finite-speed" />）を含むという CFL 条件（<Ref to="computer-science/pde-discretization/finite-difference-convergence#thm-cfl" />）と同じ条件です。1 の $\nu = 1$ の場合、各 $M$ では安定ですが、<Ref to="lem-three-term-recurrence" /> の 1 の評価の分母 $\sin\varphi$ は最高波数で $\sin(\pi h) \approx \pi h$ まで小さくなり、解の大きさの評価は $h \to 0$ で一様ではありません。$h$ について一様な評価を得るには、次節のエネルギーを使います。

## 4. 離散エネルギーと収束

前節の安定条件は固有値ごとの議論から出ましたが、$h \to 0$ で行列の次元が増える問題では、次元によらない形の評価が欲しくなります。そこで、連続のエネルギー保存に対応する量を離散化して、leapfrog 法が厳密に保存する量を探します。位置エネルギーを同じ時刻の $\langle u^n, Au^n\rangle$ ではなく、隣り合う時刻の積 $\langle u^{n+1}, Au^n\rangle$ にするのが、次の定理で保存を成り立たせる工夫です。

<Theorem id="thm-leapfrog-energy" title="leapfrog 法の離散エネルギー">
$A$ を実対称な正定値行列、$\theta = \Delta t^2\lambda_d/4$ とする。実ベクトルの列 $u^n$ が、与えられた $r^n$ について

$$
u^{n+1} - 2u^n + u^{n-1} = -\Delta t^2Au^n + \Delta t^2r^n\qquad(n \ge 1)
$$

を満たすとし、$\delta^n = (u^{n+1} - u^n)/\Delta t$、$\bar u^n = (u^{n+1} + u^n)/2$、

$$
E^n = \frac12\|\delta^n\|^2 + \frac12\langle u^{n+1}, Au^n\rangle
$$

とおく。

1. $E^n - E^{n-1} = \frac{\Delta t}{2}\langle r^n, \delta^n + \delta^{n-1}\rangle$。とくに $r^n = 0$ なら $E^n = E^0$（すべての $n$）である。
2. $E^n = \frac12\|\delta^n\|^2 - \frac{\Delta t^2}{8}\langle A\delta^n, \delta^n\rangle + \frac12\langle A\bar u^n, \bar u^n\rangle$ で、
   $$
   \frac{1-\theta}{2}\|\delta^n\|^2 + \frac12\langle A\bar u^n, \bar u^n\rangle \le E^n \le \frac12\|\delta^n\|^2 + \frac12\langle A\bar u^n, \bar u^n\rangle.
   $$
3. $\theta \lt 1$ のとき $\sqrt{E^n} \le \sqrt{E^0} + \dfrac{\Delta t}{\sqrt{2(1-\theta)}}\displaystyle\sum_{m=1}^{n}\|r^m\|$。$r^n = 0$ なら、すべての $n$ で
   $$
   \|\delta^n\|^2 \le \frac{2E^0}{1-\theta},\qquad \|u^n\| \le \sqrt{\frac{2E^0}{\lambda_1}} + \frac{\Delta t}{2}\sqrt{\frac{2E^0}{1-\theta}}.
   $$
</Theorem>

<Proof of="thm-leapfrog-energy">
1. $a = u^{n+1} - u^n$、$b = u^n - u^{n-1}$ とおくと、$a - b = -\Delta t^2Au^n + \Delta t^2r^n$、$a + b = u^{n+1} - u^{n-1}$ です。$\|a\|^2 - \|b\|^2 = \langle a - b, a + b\rangle$ なので

$$
\frac{1}{2\Delta t^2}\bigl(\|a\|^2 - \|b\|^2\bigr) = -\frac12\langle Au^n, u^{n+1} - u^{n-1}\rangle + \frac12\langle r^n, u^{n+1} - u^{n-1}\rangle
$$

です。一方 $A$ の対称性から $\langle u^{n+1}, Au^n\rangle - \langle u^n, Au^{n-1}\rangle = \langle Au^n, u^{n+1} - u^{n-1}\rangle$ です。2 つを足すと $E^n - E^{n-1} = \frac12\langle r^n, u^{n+1} - u^{n-1}\rangle$ で、$u^{n+1} - u^{n-1} = \Delta t(\delta^n + \delta^{n-1})$ です。

2. $u^{n+1} = \bar u^n + \frac{\Delta t}{2}\delta^n$、$u^n = \bar u^n - \frac{\Delta t}{2}\delta^n$ を入れると、$A$ の対称性から交差項が消えて $\langle u^{n+1}, Au^n\rangle = \langle A\bar u^n, \bar u^n\rangle - \frac{\Delta t^2}{4}\langle A\delta^n, \delta^n\rangle$ です。$0 \le \langle A\delta^n, \delta^n\rangle \le \lambda_d\|\delta^n\|^2$ から $0 \le \frac{\Delta t^2}{8}\langle A\delta^n,\delta^n\rangle \le \frac\theta2\|\delta^n\|^2$ で、不等式が出ます。

3. $\theta \lt 1$ なら 2 から $E^n \ge 0$、$\|\delta^n\| \le \sqrt{2E^n/(1-\theta)}$ です。1 と Cauchy–Schwarz の不等式から

$$
\bigl(\sqrt{E^n} - \sqrt{E^{n-1}}\bigr)\bigl(\sqrt{E^n} + \sqrt{E^{n-1}}\bigr) = E^n - E^{n-1} \le \frac{\Delta t}{2}\|r^n\|\sqrt{\frac{2}{1-\theta}}\bigl(\sqrt{E^n} + \sqrt{E^{n-1}}\bigr)
$$

です。$\sqrt{E^n} + \sqrt{E^{n-1}} \gt 0$ なら割って $\sqrt{E^n} - \sqrt{E^{n-1}} \le \Delta t\|r^n\|/\sqrt{2(1-\theta)}$、0 なら両辺 0 でやはり成り立ちます。$n$ について足すと最初の式です。$r^n = 0$ なら $E^n = E^0$ で、2 から $\|\delta^n\|^2 \le 2E^0/(1-\theta)$、$\lambda_1\|\bar u^n\|^2 \le \langle A\bar u^n, \bar u^n\rangle \le 2E^0$ です。$u^n = \bar u^n - \frac{\Delta t}{2}\delta^n$ に三角不等式を使うと最後の式が出ます。
</Proof>

$E^n$ は連続のエネルギー $\frac12\|u'\|^2 + \frac12\langle u, Au\rangle$ の離散版です。位置エネルギーの項が $\langle u^n, Au^n\rangle$ ではなく、隣り合う時刻の積 $\langle u^{n+1}, Au^n\rangle$ になっている点が要で、このために 2 の第 2 項（負の項）が現れ、正定値であるための条件 $\theta \lt 1$ が <Ref to="thm-leapfrog-stability" /> の安定条件と一致します。3 の定数は $A$ の次元によらず $\theta$ と $\lambda_1$ だけで決まるので、$\nu \le \nu_0 \lt 1$ と固定すれば、$h \to 0$ で一様な安定性の評価になります。

<Theorem id="thm-leapfrog-convergence" title="波動方程式の leapfrog 法の収束">
$u \in C^4([0,1]\times[0,T])$ が $u_{tt} = c^2u_{xx}$、$u(0,t) = u(1,t) = 0$、$u(x,0) = g(x)$、$u_t(x,0) = g_1(x)$ を満たすとする。$h = 1/M$、$\nu = c\Delta t/h \le \nu_0 \lt 1$ とし、$U^0_j = g(x_j)$、

$$
U^1_j = g(x_j) + \Delta t\,g_1(x_j) + \frac{c^2\Delta t^2}{2}\cdot\frac{g(x_{j+1}) - 2g(x_j) + g(x_{j-1})}{h^2}
$$

を出発値とする leapfrog・中心差分の解を $U^n$ とする。$B_t = \max|\partial_t^4u|$、$B_x = \max|\partial_x^4u|$、$B_3 = \max|\partial_t^3u|$ とおくと、$t_n \le T$ である $n$ について

$$
\|u(\cdot, t_n) - U^n\|_h \le \frac{t_n}{\sqrt{1-\nu_0^2}}\Bigl(\frac{c^2\Delta t\,h^2}{24}B_x + \frac{\Delta t^2}{6}B_3\Bigr) + \frac{t_n^2}{1-\nu_0^2}\Bigl(\frac{\Delta t^2}{12}B_t + \frac{c^2h^2}{12}B_x\Bigr)
$$

が成り立つ（左辺の $u(\cdot,t_n)$ は内部の格子点での値のベクトル）。右辺は $O(\Delta t^2 + h^2)$ である。

また、格子 $h\mathbb{Z}$ の上で $c\Delta t = h$（$\nu = 1$）のとき、$u \in C^2(\mathbb{R}\times[0,\infty))$ が $u_{tt} = c^2u_{xx}$ を満たし、$U^0_j = u(x_j, 0)$、$U^1_j = u(x_j, \Delta t)$ ならば、すべての $n \ge 0$ と $j$ で $U^n_j = u(x_j, t_n)$ である。
</Theorem>

<Proof of="thm-leapfrog-convergence">
**誤差の方程式。** $e^n_j = u(x_j, t_n) - U^n_j$ とおきます。2 階中心差分の打ち切り誤差（<Ref to="computer-science/pde-discretization/finite-difference-convergence#lem-difference-truncation" /> の 3）を $t$ と $x$ のそれぞれについて使うと、内部の格子点と $n \ge 1$ で

$$
\frac{u(x_j,t_{n+1}) - 2u(x_j,t_n) + u(x_j,t_{n-1})}{\Delta t^2} = c^2\frac{u(x_{j+1},t_n) - 2u(x_j,t_n) + u(x_{j-1},t_n)}{h^2} + r^n_j,\qquad |r^n_j| \le \rho := \frac{\Delta t^2}{12}B_t + \frac{c^2h^2}{12}B_x
$$

です（$u_{tt} = c^2u_{xx}$ を使った）。$U^n$ の式を引くと、$e^n$ は <Ref to="thm-leapfrog-energy" /> の式を $A = \frac{c^2}{h^2}\mathrm{tridiag}(-1,2,-1)$ と $r^n$ について満たします。内部の格子点は $M - 1$ 個で $h(M-1) \lt 1$ なので $\|r^n\|_h \le \rho$ です。

**出発値の誤差。** $e^0 = 0$ です。$t$ についての Taylor の定理（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-taylor" />）と $u_{tt}(x,0) = c^2g''(x)$ から $u(x_j,\Delta t) = g(x_j) + \Delta t\,g_1(x_j) + \frac{c^2\Delta t^2}{2}g''(x_j) + \frac{\Delta t^3}{6}\partial_t^3u(x_j,\xi)$ なので、再び打ち切り誤差の補題から

$$
|e^1_j| \le \frac{c^2\Delta t^2}{2}\cdot\frac{h^2}{12}B_x + \frac{\Delta t^3}{6}B_3,\qquad \|\delta^0\|_h = \frac{\|e^1\|_h}{\Delta t} \le d_0 := \frac{c^2\Delta t\,h^2}{24}B_x + \frac{\Delta t^2}{6}B_3
$$

です。

**エネルギーによる評価。** <Ref to="thm-leapfrog-energy" /> を、内積を $\langle\cdot,\cdot\rangle_h$ に替えて使います（すべての式が $h$ 倍されるだけです）。$e^0 = 0$ なので $E^0 = \frac12\|\delta^0\|_h^2 \le d_0^2/2$。$\theta = \Delta t^2\lambda_{M-1}/4 \le \nu^2 \le \nu_0^2$ なので、同じ定理の 3 と 2 から

$$
\|\delta^n\|_h \le \sqrt{\frac{2E^n}{1-\theta}} \le \frac{1}{\sqrt{1-\nu_0^2}}\Bigl(d_0 + \frac{t_n\rho}{\sqrt{1-\nu_0^2}}\Bigr)
$$

です。$e^n = e^0 + \Delta t\sum_{m=0}^{n-1}\delta^m$ なので $\|e^n\|_h \le t_n\max_{m \lt n}\|\delta^m\|_h$ で、主張の式が出ます。

**$\nu = 1$ の場合。** $s = ct$ を時間の変数にとると $u_{ss} = u_{xx}$ なので、d'Alembert の公式（<Ref to="mathematics/pde/wave-equation#thm-dalembert" />）から、$C^2$ 級の関数 $F, G$ があって $u(x,t) = F(x + ct) + G(x - ct)$ と書けます（$g = u(\cdot,0)$、$g_1 = u_t(\cdot,0)$ として $F(y) = \frac12g(y) + \frac1{2c}\int_0^yg_1$、$G(y) = \frac12g(y) - \frac1{2c}\int_0^yg_1$）。$c\Delta t = h$ なので $x_j + ct_{n\pm1} = x_{j\pm1} + ct_n$、$x_j - ct_{n\pm1} = x_{j\mp1} - ct_n$ で、

$$
u(x_j,t_{n+1}) + u(x_j,t_{n-1}) = F(x_{j+1} + ct_n) + G(x_{j-1} - ct_n) + F(x_{j-1} + ct_n) + G(x_{j+1} - ct_n) = u(x_{j+1},t_n) + u(x_{j-1},t_n)
$$

です。$\nu = 1$ のスキームは $U^{n+1}_j = U^n_{j+1} + U^n_{j-1} - U^{n-1}_j$ なので、厳密解の格子値は同じ漸化式を満たし、$n = 0, 1$ で一致するので、すべての $n$ で一致します。
</Proof>

<Example id="ex-leapfrog-wave" title="CFL 数と離散エネルギー">
$c = 1$、初期値 $g(x) = \sin\pi x + \frac12\sin3\pi x$、初速 0 の定在波（厳密解は $\sin\pi x\cos\pi t + \frac12\sin3\pi x\cos3\pi t$）を、$\nu = 0.5,\ 1.0,\ 1.01$ で解きます。

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

def leapfrog(M, nu, T):
    """u_tt = u_xx（0 < x < 1、両端 0）を leapfrog・中心差分で解く。h = 1/M、dt = nu h。"""
    h, dt = 1.0 / M, nu / M
    x = np.linspace(0, 1, M + 1)
    g = np.sin(np.pi * x) + 0.5 * np.sin(3 * np.pi * x)
    lap = lambda v: np.concatenate(([0.0], v[2:] - 2 * v[1:-1] + v[:-2], [0.0]))
    u_old, u = g, g + 0.5 * nu**2 * lap(g)            # 初速 0、2 次精度の出発値
    energy = lambda a, b: 0.5 * h * np.sum(((b - a) / dt) ** 2) - 0.5 * np.sum(b * lap(a)) / h
    E0, n, maxima = energy(u_old, u), 1, {}
    steps = int(round(T / dt))
    while n < steps:
        u_old, u = u, 2 * u - u_old + nu**2 * lap(u)
        n += 1
        if n % (steps // 4) == 0:
            maxima[round(n * dt, 6)] = np.max(np.abs(u))
    exact = np.sin(np.pi * x) * np.cos(np.pi * n * dt) + 0.5 * np.sin(3 * np.pi * x) * np.cos(3 * np.pi * n * dt)
    return maxima, abs(energy(u_old, u) - E0) / E0, np.sqrt(h * np.sum((u - exact) ** 2))

print("M = 200, T = 4:   max|u| at t = 1, 2, 3, 4            |E^n - E^0| / E^0")
for nu in [0.5, 1.0, 1.01]:
    maxima, dE, _ = leapfrog(200, nu, 4.0)
    print(f"nu = {nu:4.2f}:  " + "  ".join(f"{v:.3e}" for v in maxima.values()) + f"      {dE:.1e}")

print("nu = 0.5, T = 0.7:   M     l2 error     order")
prev = None
for M in [25, 50, 100, 200, 400]:
    err = leapfrog(M, 0.5, 0.7)[2]
    print(f"                {M:4d}    {err:.3e}    " + ("" if prev is None else f"{np.log2(prev / err):.2f}"))
    prev = err
print("nu = 1.0, T = 0.7, M = 100: l2 error =", f"{leapfrog(100, 1.0, 0.7)[2]:.1e}")
```

```
M = 200, T = 4:   max|u| at t = 1, 2, 3, 4            |E^n - E^0| / E^0
nu = 0.50:  1.076e+00  1.076e+00  1.076e+00  1.076e+00      5.3e-15
nu = 1.00:  1.076e+00  1.076e+00  1.076e+00  1.076e+00      3.3e-15
nu = 1.01:  4.271e+07  6.712e+31  1.135e+56  1.962e+80      1.0e+00
nu = 0.5, T = 0.7:   M     l2 error     order
                  25    3.125e-03    
                  50    8.068e-04    1.95
                 100    2.033e-04    1.99
                 200    5.092e-05    2.00
                 400    1.274e-05    2.00
nu = 1.0, T = 0.7, M = 100: l2 error = 3.5e-16
```

$\nu = 0.5$ と $1.0$ では最大値は変わらず、離散エネルギーは丸めの大きさで保存されています。$\nu = 1.01$ では、丸め誤差から入った最高波数の成分が増幅されます。<Ref to="lem-three-term-recurrence" /> の 3 で $\mu = 4\nu^2\sin^2(199\pi/400) \approx 4.0801$ とすると $|z| \approx 1.326$ で、時間 1（198 ステップ）あたり $|z|^{198} \approx 1.8\times10^{24}$ 倍です。表の隣り合う値の比（$1.6\times10^{24}$、$1.7\times10^{24}$、$1.7\times10^{24}$）と合います。このとき $E^n$ は厳密な計算では保存されますが、正定値ではないので解の大きさを抑えません（増大する成分は $\frac12\|\delta\|^2$ と $\frac12\langle u^{n+1}, Au^n\rangle$ が打ち消し合ってエネルギーが 0 です）。表の $1.0$ は、$10^{160}$ の大きさの 2 つの項の差が浮動小数点数で計算できなくなった結果です。$\nu = 0.5$ の誤差の次数は 2、$\nu = 1$ では <Ref to="thm-leapfrog-convergence" /> の後半のとおり丸めの大きさです（区間の問題は、奇関数として周期 2 で延長すれば全直線の場合に帰着します。初速が 0 のとき、出発値 $U^1_j = (g(x_{j+1}) + g(x_{j-1}))/2$ も d'Alembert の公式の値に一致します）。
</Example>

## 5. 数値分散とシンプレクティック性

leapfrog 法は各固有成分の振幅を保ちます（<Ref to="lem-three-term-recurrence" /> の 1 で、解は $\cos n\varphi$ と $\sin n\varphi$ の 1 次結合）。誤差はすべて位相に現れます。

<Proposition id="prop-numerical-dispersion" title="leapfrog・中心差分の数値分散関係">
格子 $h\mathbb{Z}$ 上のスキーム $u^{n+1}_j - 2u^n_j + u^{n-1}_j = \nu^2(u^n_{j+1} - 2u^n_j + u^n_{j-1})$（$\nu = c\Delta t/h$）について次が成り立つ。

1. 実数 $k$、$\omega$ について、平面波 $u^n_j = e^{i(kx_j - \omega t_n)}$ が解であることは
   $$
   \sin^2\frac{\omega\Delta t}{2} = \nu^2\sin^2\frac{kh}{2}
   $$
   と同値である。
2. $0 \lt \nu \le 1$ を固定し、$|kh| \le \pi$ について $\omega_h(k) = \frac{2}{\Delta t}\arcsin\bigl(\nu\sin\frac{kh}{2}\bigr)$ とおくと、$kh \to 0$ で
   $$
   \frac{\omega_h(k)}{ck} = 1 - \frac{1-\nu^2}{24}(kh)^2 + O\bigl((kh)^4\bigr)
   $$
   である。$\nu = 1$ では $\omega_h(k) = ck$ が厳密に成り立つ。
</Proposition>

<Proof of="prop-numerical-dispersion">
1. 平面波を代入して $e^{i(kx_j - \omega t_n)}$ で割ると、左辺は $e^{-i\omega\Delta t} - 2 + e^{i\omega\Delta t} = 2\cos\omega\Delta t - 2 = -4\sin^2(\omega\Delta t/2)$、右辺は $\nu^2(2\cos kh - 2) = -4\nu^2\sin^2(kh/2)$ です。

2. $G(\vartheta) = \arcsin(\nu\sin\vartheta)$ とおきます。$\nu = 1$ なら $|\vartheta| \le \pi/2$ で $G(\vartheta) = \vartheta$ なので $\omega_h(k) = \frac{2}{\Delta t}\cdot\frac{kh}{2} = ck$（$h/\Delta t = c$）です。$\nu \lt 1$ のとき、$G$ は $\mathbb{R}$ 上で $C^\infty$ 級の奇関数で、$G'(\vartheta) = \nu\cos\vartheta\,(1 - \nu^2\sin^2\vartheta)^{-1/2}$ です。$\cos\vartheta = 1 - \vartheta^2/2 + O(\vartheta^4)$、$(1 - \nu^2\sin^2\vartheta)^{-1/2} = 1 + \nu^2\vartheta^2/2 + O(\vartheta^4)$ から $G'(\vartheta) = \nu - \frac{\nu(1-\nu^2)}{2}\vartheta^2 + O(\vartheta^4)$ で、$G'(0) = \nu$、$G'''(0) = -\nu(1-\nu^2)$ です。$G$ は奇関数なので $G(0) = G''(0) = G''''(0) = 0$ で、Taylor の定理（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-taylor" />）から $G(\vartheta) = \nu\vartheta - \frac{\nu(1-\nu^2)}{6}\vartheta^3 + O(\vartheta^5)$ です。$\vartheta = kh/2$ とし、$2\nu/\Delta t\cdot kh/2 = ck$ を使うと

$$
\omega_h(k) = \frac{2}{\Delta t}G\Bigl(\frac{kh}{2}\Bigr) = ck\Bigl(1 - \frac{1-\nu^2}{24}(kh)^2 + O\bigl((kh)^4\bigr)\Bigr)
$$

です。
</Proof>

厳密な分散関係は $\omega = ck$ で、位相速度 $\omega/k = c$ は波数によりません。数値解では位相速度が $c\bigl(1 - \frac{1-\nu^2}{24}(kh)^2\bigr)$ で、波長が短いほど遅く進みます（**数値分散**）。時刻 $t$ までの位相の誤差は $(\omega_h - ck)t \approx -\frac{1-\nu^2}{24}c\,k^3h^2t$ で、$t$ に比例して増えます。波長あたりの格子点の数を $n_\lambda = 2\pi/(kh)$ とすると、位相の誤差を $\varepsilon$ 以下に保つには、波が進む波長の数 $ckt/(2\pi)$ が $\frac{3n_\lambda^2\varepsilon}{\pi^3(1-\nu^2)}$ 以下でなければなりません（$n_\lambda = 20$、$\nu = 0.5$、$\varepsilon = 0.1$ なら約 5 波長）。長い距離を伝わる波の計算で、2 次の差分より高次の差分やスペクトル法が使われる理由です。

leapfrog 法が振幅を保つことは、対称行列の線形系に限らない構造から来ています。その構造とは、位置と運動量の組 $(q, p)$ を相空間の 1 点とみたとき、1 ステップの写像が体積を保つだけでなく、正準 2 形式という量までまるごと保つこと（シンプレクティック性。$d = 1$ では面積の保存と同じことです）です。次の命題は、leapfrog 法がこの性質をもつことと、調和振動子で数値解がどのような振動になるかを示します。

<Proposition id="prop-leapfrog-symplectic" title="leapfrog 法のシンプレクティック性">
$A$ を $d$ 次の実対称行列とし、$\mathbb{R}^{2d}$ 上の写像 $\Phi_{\Delta t} : (q, p) \mapsto (q', p')$ を

$$
p^{1/2} = p - \tfrac{\Delta t}{2}Aq,\qquad q' = q + \Delta t\,p^{1/2},\qquad p' = p^{1/2} - \tfrac{\Delta t}{2}Aq'
$$

で定める（速度 Verlet 法）。$(q^{n+1}, p^{n+1}) = \Phi_{\Delta t}(q^n, p^n)$ とする。

1. $q^n$ は leapfrog 法の式 $q^{n+1} - 2q^n + q^{n-1} = -\Delta t^2Aq^n$（$n \ge 1$）を満たす。
2. $\Phi_{\Delta t}$ は線形で、その行列は $J = \begin{pmatrix}0 & I\\ -I & 0\end{pmatrix}$ について $\Phi_{\Delta t}^{\mathsf T}J\Phi_{\Delta t} = J$ を満たし（シンプレクティック）、$\det\Phi_{\Delta t} = 1$、$\Phi_{-\Delta t} = \Phi_{\Delta t}^{-1}$（時間反転対称）である。
3. $d = 1$、$A = \omega^2$、$0 \lt \omega\Delta t \lt 2$ のとき、$q^n$ は振動数 $\tilde\omega = \frac{2}{\Delta t}\arcsin\frac{\omega\Delta t}{2}$ の調和振動 $q^n = a\cos\tilde\omega t_n + b\sin\tilde\omega t_n$ で、$p^2 + \omega^2\bigl(1 - \frac{\omega^2\Delta t^2}{4}\bigr)q^2$ は $\Phi_{\Delta t}$ で不変である。
</Proposition>

<Proof of="prop-leapfrog-symplectic">
1. 第 $n$ ステップの中間の値を $p^{n+1/2}$ と書くと、$q^{n+1} - q^n = \Delta t\,p^{n+1/2}$、$q^n - q^{n-1} = \Delta t\,p^{n-1/2}$ です。$p^n = p^{n-1/2} - \frac{\Delta t}{2}Aq^n$、$p^{n+1/2} = p^n - \frac{\Delta t}{2}Aq^n$ から $p^{n+1/2} - p^{n-1/2} = -\Delta tAq^n$ で、差をとると主張の式です。

2. 3 つの段階はそれぞれ線形で、行列は

$$
L = \begin{pmatrix}I & 0\\ -\frac{\Delta t}{2}A & I\end{pmatrix},\qquad W = \begin{pmatrix}I & \Delta t\,I\\ 0 & I\end{pmatrix},\qquad \Phi_{\Delta t} = L\,W\,L
$$

です。$C$ を対称行列として $L_C = \begin{pmatrix}I & 0\\ C & I\end{pmatrix}$ とすると

$$
L_C^{\mathsf T}JL_C = \begin{pmatrix}I & C\\ 0 & I\end{pmatrix}\begin{pmatrix}0 & I\\ -I & 0\end{pmatrix}\begin{pmatrix}I & 0\\ C & I\end{pmatrix} = \begin{pmatrix}-C & I\\ -I & 0\end{pmatrix}\begin{pmatrix}I & 0\\ C & I\end{pmatrix} = \begin{pmatrix}0 & I\\ -I & 0\end{pmatrix}
$$

で、同様に $W_B = \begin{pmatrix}I & B\\ 0 & I\end{pmatrix}$（$B$ は対称）について $W_B^{\mathsf T}JW_B = \begin{pmatrix}0 & I\\ -I & B\end{pmatrix}\begin{pmatrix}I & B\\ 0 & I\end{pmatrix} = J$ です。$X^{\mathsf T}JX = J$、$Y^{\mathsf T}JY = J$ なら $(XY)^{\mathsf T}J(XY) = Y^{\mathsf T}JY = J$ なので、積 $LWL$ もこの式を満たします。$L$ と $W$ は対角成分が 1 のブロック三角行列なので行列式は 1 です。$L_C^{-1} = L_{-C}$、$W_B^{-1} = W_{-B}$ なので、$\Delta t$ を $-\Delta t$ に替えた $\Phi_{-\Delta t} = L^{-1}W^{-1}L^{-1} = (LWL)^{-1}$ です。

3. 1 と <Ref to="lem-three-term-recurrence" /> の 1（$\mu = \omega^2\Delta t^2 \lt 4$）から、$q^n$ は $\cos n\varphi$ と $\sin n\varphi$ の 1 次結合で、$\cos\varphi = 1 - \omega^2\Delta t^2/2$ です。$1 - \cos\varphi = 2\sin^2(\varphi/2)$ から $\sin(\varphi/2) = \omega\Delta t/2$、$\varphi = \tilde\omega\Delta t$ で、$n\varphi = \tilde\omega t_n$ です。不変量: $s = \omega^2\Delta t^2/2$、$\kappa = \omega^2(1 - s/2)$ とおくと、定義から $q' = (1 - s)q + \Delta t\,p$、$p' = -\Delta t\,\kappa\,q + (1 - s)p$ で、

$$
p'^2 + \kappa q'^2 = \kappa\bigl(\kappa\Delta t^2 + (1-s)^2\bigr)q^2 + \bigl((1-s)^2 + \kappa\Delta t^2\bigr)p^2
$$

です（$pq$ の項は $p'^2$ から $-2\Delta t\,\kappa(1-s)pq$、$\kappa q'^2$ から $+2\Delta t\,\kappa(1-s)pq$ が出て消える）。$\kappa\Delta t^2 = 2s(1 - s/2) = 2s - s^2$ なので $(1-s)^2 + \kappa\Delta t^2 = 1$ で、右辺は $\kappa q^2 + p^2$ です。
</Proof>

2 の $\Phi^{\mathsf T}J\Phi = J$ は、相空間の正準 2 形式 $\sum_jdq_j\wedge dp_j$ を保つことを意味します。連続の流れ（Hamilton 系の解）がもつ性質で、行列式が 1 であること（相空間の体積の保存。連続の場合は <Ref to="mathematics/ode/linear-systems#thm-liouville" />）より強い条件です。3 は、leapfrog 法の数値解が「振動数を $\tilde\omega = \omega(1 + \omega^2\Delta t^2/24 + \cdots)$ に取り替えた調和振動子」の厳密解であることを言っています。誤差を「少し違う問題を厳密に解いている」と読む見方（後退誤差解析）の最も簡単な場合で、保存量 $p^2 + \omega^2(1 - \omega^2\Delta t^2/4)q^2$ は修正されたエネルギーです。もとのエネルギー $p^2 + \omega^2q^2$ は保存されませんが、この不変量があるので $O(\Delta t^2)$ の幅で振動するだけで、時間とともにずれていくことはありません。

<Remark id="rem-backward-error" title="非線形の系での後退誤差解析（文献）">
非線形のポテンシャル $q'' = -\nabla V(q)$ でも、速度 Verlet 法はシンプレクティックで時間反転対称です。このとき数値解は、修正されたハミルトニアン $\tilde H = H + \Delta t^2H_2 + \cdots$ の流れに指数的に小さい誤差で従い、エネルギーの誤差は指数的に長い時間 $O(\Delta t^2)$ にとどまります（Benettin–Giorgilli、Hairer–Lubich）。証明は Hairer・Lubich・Wanner の教科書（参考文献 [1]）の第 IX 章にあります。この章では線形の場合に限り、定理環境にはしません。
</Remark>

## 6. Schrödinger 方程式：Cayley 変換と陽的 Euler 法

<Definition id="def-cayley" title="Crank–Nicolson 法（Cayley 変換）">
半離散の Schrödinger 方程式に対し

$$
\psi^{n+1} = C_{\Delta t}\psi^n,\qquad C_{\Delta t} = \Bigl(I + \frac{i\Delta t}{2}H\Bigr)^{-1}\Bigl(I - \frac{i\Delta t}{2}H\Bigr)
$$

で定まる方法を **Crank–Nicolson 法**といい、$C_{\Delta t}$ を $H$ の **Cayley 変換**という。$\psi^{n+1} = (I - i\Delta tH)\psi^n$ を**陽的 Euler 法**、$\psi^{n+1} = (I + i\Delta tH)^{-1}\psi^n$ を**陰的 Euler 法**という。
</Definition>

Crank–Nicolson 法は、$i(\psi^{n+1} - \psi^n)/\Delta t = H(\psi^{n+1} + \psi^n)/2$ を $\psi^{n+1}$ について解いたものです。熱方程式に対する $\theta = 1/2$ の $\theta$ 法（<Ref to="computer-science/pde-discretization/elliptic-parabolic-difference#thm-theta-stability" />）と同じ形で、係数が純虚数になっています。$H$ が三重対角なら、1 ステップは三重対角の連立方程式 1 回です。

<Lemma id="lem-cayley-scalar" title="Cayley 変換のスカラーの評価">
実数 $x$ について $r(x) = \dfrac{1 - ix/2}{1 + ix/2}$ とおくと、$|r(x)| = 1$、$r(x) = e^{-2i\arctan(x/2)}$ で、

$$
\bigl|r(x) - e^{-ix}\bigr| \le \frac{|x|^3}{12}
$$

である。
</Lemma>

<Proof of="lem-cayley-scalar">
$1 - ix/2 = \sqrt{1 + x^2/4}\,e^{-i\arctan(x/2)}$ で、$1 + ix/2$ はその複素共役なので、$r(x) = e^{-2i\arctan(x/2)}$、$|r(x)| = 1$ です。実数 $a, b$ について $|e^{-ia} - e^{-ib}| = \bigl|\int_b^ae^{-is}ds\bigr| \le |a - b|$ なので、$|r(x) - e^{-ix}| \le |x - 2\arctan(x/2)|$ です。$f(x) = x - 2\arctan(x/2)$ は $f(0) = 0$、

$$
f'(x) = 1 - \frac{1}{1 + x^2/4} = \frac{x^2/4}{1 + x^2/4} \in \Bigl[0, \frac{x^2}{4}\Bigr]
$$

なので $|f(x)| \le \bigl|\int_0^xs^2/4\,ds\bigr| = |x|^3/12$ です。
</Proof>

<Theorem id="thm-cayley-unitary" title="Crank–Nicolson 法のユニタリ性と誤差">
$H$ をエルミート行列、$\Delta t \gt 0$ とする。

1. $C_{\Delta t}$ はユニタリで $H$ と可換である。したがって $\|\psi^n\| = \|\psi_0\|$、$\langle\psi^n, H\psi^n\rangle = \langle\psi_0, H\psi_0\rangle$（すべての $n$）。
2. $t_n = n\Delta t$ について
   $$
   \bigl\|\psi^n - e^{-it_nH}\psi_0\bigr\| \le \frac{t_n\,\Delta t^2}{12}\,\|H^3\psi_0\|.
   $$
</Theorem>

<Proof of="thm-cayley-unitary">
$H = U\Lambda U^{*}$、$\Lambda = \mathrm{diag}(\mu_k)$ とします（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）。$I \pm \frac{i\Delta t}{2}H = U(I \pm \frac{i\Delta t}{2}\Lambda)U^{*}$ で、$1 + i\Delta t\mu_k/2 \ne 0$ なので逆行列が存在し、$C_{\Delta t} = U\,\mathrm{diag}\bigl(r(\Delta t\mu_k)\bigr)U^{*}$ です。

1. <Ref to="lem-cayley-scalar" /> から $|r(\Delta t\mu_k)| = 1$ なので $C_{\Delta t}$ はユニタリです（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-normal" />）。$H$ と同じ $U$ で対角化されるので可換で、$\langle C_{\Delta t}\psi, HC_{\Delta t}\psi\rangle = \langle C_{\Delta t}\psi, C_{\Delta t}H\psi\rangle = \langle\psi, H\psi\rangle$（<Ref to="mathematics/linear-algebra/spectral-theorem#prop-unitary" />）です。

2. $\psi_0 = \sum_kc_kq_k$（$q_k$ は $U$ の列）と展開すると、<Ref to="lem-unitary-group" /> から

$$
\psi^n - e^{-it_nH}\psi_0 = \sum_kc_k\bigl(r(\Delta t\mu_k)^n - e^{-in\Delta t\mu_k}\bigr)q_k
$$

です。絶対値 1 の複素数 $a, b$ について $|a^n - b^n| \le n|a - b|$（<Ref to="lem-unitary-group" /> の最後の主張の $1\times1$ の場合）なので、<Ref to="lem-cayley-scalar" /> から各係数の絶対値は $|c_k|\,n|\Delta t\mu_k|^3/12$ 以下です。よって

$$
\bigl\|\psi^n - e^{-it_nH}\psi_0\bigr\|^2 \le \Bigl(\frac{n\Delta t^3}{12}\Bigr)^2\sum_k\mu_k^6|c_k|^2 = \Bigl(\frac{t_n\Delta t^2}{12}\Bigr)^2\|H^3\psi_0\|^2
$$

です。
</Proof>

誤差の評価には $\|H\|$ ではなく $\|H^3\psi_0\|$ が現れます。$H$ が微分作用素の離散化なら $\|H\| = O(h^{-2})$ ですが、$\psi_0$ が滑らかな関数の格子値なら $\|H^3\psi_0\|$ は $h$ によらず有界で（$\psi_0$ の 6 階までの導関数で抑えられる）、$h$ と $\Delta t$ を独立に小さくできます。安定性に刻み幅の条件がないこと（無条件安定）と合わせて、これが陰的な方法の利点です。一方、高い波数の成分については $\Delta t\mu_k$ が大きく、$r(\Delta t\mu_k) \to -1$ となって位相がまったく合いません。Crank–Nicolson 法は高い波数の成分を減衰させず、間違った位相で保ち続けます（熱方程式の場合の <Ref to="computer-science/pde-discretization/elliptic-parabolic-difference#prop-cn-damping" /> と同じ現象です）。

<Proposition id="prop-explicit-euler-unstable" title="陽的 Euler 法と陰的 Euler 法">
$H$ をエルミート行列、$\varrho = \max_k|\mu_k|$ をそのスペクトル半径とする。

1. $\|(I - i\Delta tH)^n\| = (1 + \Delta t^2\varrho^2)^{n/2}$。$H \ne 0$ なら、どんな $\Delta t \gt 0$ でも $n \to \infty$ で発散する。$n\Delta t = t$ を固定すると $(1 + \Delta t^2\varrho^2)^{n/2} \le e^{t\Delta t\varrho^2/2}$ である。
2. $\|(I + i\Delta tH)^{-n}\psi\| \le \|\psi\|$ で、$H\psi = \mu\psi$（$\mu \ne 0$）なら $\|(I + i\Delta tH)^{-n}\psi\| = (1 + \Delta t^2\mu^2)^{-n/2}\|\psi\|$ である。
</Proposition>

<Proof of="prop-explicit-euler-unstable">
$H = U\Lambda U^{*}$ から $(I - i\Delta tH)^n = U\,\mathrm{diag}\bigl((1 - i\Delta t\mu_k)^n\bigr)U^{*}$ で、ユニタリ行列を掛けても 2 ノルムは変わらないので、ノルムは $\max_k|1 - i\Delta t\mu_k|^n = (1 + \Delta t^2\varrho^2)^{n/2}$ です。$1 + y \le e^y$ から $(1 + \Delta t^2\varrho^2)^{n/2} \le e^{n\Delta t^2\varrho^2/2} = e^{t\Delta t\varrho^2/2}$ です。陰的 Euler 法では対角成分が $(1 + i\Delta t\mu_k)^{-n}$ で、絶対値は $(1 + \Delta t^2\mu_k^2)^{-n/2} \le 1$ です。
</Proof>

陽的 Euler 法の増幅因子 $1 - i\Delta t\mu$ は、単位円周上の点 $e^{-i\Delta t\mu}$ を接線の方向に近似したもので、必ず円の外に出ます。$\varrho = O(h^{-2})$ なので、時刻 $t$ までノルムの増大を抑えるには $\Delta t = O(h^4)$ が要り、実用になりません。陰的 Euler 法は逆に円の内側に入り、ノルムが減ります。どちらも確率の保存を壊します。

<Figure caption="1 ステップの増幅因子（刻み幅と固有値の積が 1.5 の場合）。厳密な解は単位円周上に載り、Crank–Nicolson 法も円周上（位相だけがずれる）、陽的 Euler 法は円の外、陰的 Euler 法は円の内側に載る。">
<svg viewBox="0 0 520 300" width="100%" role="img" aria-label="複素平面の単位円と、増幅因子の位置。厳密な解と Crank–Nicolson 法は円周上、陽的 Euler 法は円の外、陰的 Euler 法は円の内側にある">
  <line x1="60" y1="120" x2="360" y2="120" stroke="currentColor" stroke-width="1" opacity="0.5" />
  <line x1="200" y1="10" x2="200" y2="285" stroke="currentColor" stroke-width="1" opacity="0.5" />
  <circle cx="200" cy="120" r="100" fill="none" stroke="currentColor" stroke-width="1.5" />
  <text x="283" y="44" font-size="13" fill="currentColor">単位円（絶対値 1）</text>
  <text x="364" y="116" font-size="12" fill="currentColor">実軸</text>
  <text x="206" y="20" font-size="12" fill="currentColor">虚軸</text>
  <circle cx="207" cy="220" r="5" fill="currentColor" />
  <text x="196" y="240" text-anchor="end" font-size="13" fill="currentColor">厳密</text>
  <circle cx="228" cy="216" r="5" fill="var(--sl-color-accent)" />
  <text x="240" y="214" font-size="13" fill="currentColor">Crank–Nicolson</text>
  <circle cx="231" cy="166" r="5" fill="none" stroke="currentColor" stroke-width="2" />
  <text x="243" y="166" font-size="13" fill="currentColor">陰的 Euler</text>
  <circle cx="300" cy="270" r="5" fill="none" stroke="var(--sl-color-accent)" stroke-width="2" />
  <text x="312" y="274" font-size="13" fill="currentColor">陽的 Euler</text>
</svg>
</Figure>

<Example id="ex-rabi" title="2 準位系の Rabi 振動">
$H = \frac\Omega2\begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}$、$\psi_0 = (1, 0)$ とすると、固有値は $\pm\Omega/2$ で、第 2 成分の確率は $|\psi_2(t)|^2 = \sin^2(\Omega t/2)$（Rabi 振動）です。$\Omega = 1$、$T = 10$ で 3 つの方法のノルムを比べます。

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

Omega, T = 1.0, 10.0
H = 0.5 * Omega * np.array([[0, 1], [1, 0]], dtype=complex)     # 固有値 ±Omega/2
I = np.eye(2)
print(" steps    explicit Euler ||psi||   (1 + dt^2 Omega^2/4)^(n/2)   implicit Euler ||psi||   Crank-Nicolson ||psi||   P_2 error (CN)")
for steps in [100, 1000, 10000]:
    dt = T / steps
    pe = pi = pc = np.array([1.0, 0.0], dtype=complex)
    Mi = np.linalg.inv(I + 1j * dt * H)
    Mc = np.linalg.solve(I + 0.5j * dt * H, I - 0.5j * dt * H)
    for _ in range(steps):
        pe, pi, pc = pe - 1j * dt * (H @ pe), Mi @ pi, Mc @ pc
    print(f"{steps:6d}    {np.linalg.norm(pe):.10f}             {(1 + dt**2 * Omega**2 / 4) ** (steps / 2):.10f}"
          f"                 {np.linalg.norm(pi):.10f}             {np.linalg.norm(pc):.10f}"
          f"             {abs(abs(pc[1])**2 - np.sin(Omega * T / 2)**2):.2e}")
```

```
 steps    explicit Euler ||psi||   (1 + dt^2 Omega^2/4)^(n/2)   implicit Euler ||psi||   Crank-Nicolson ||psi||   P_2 error (CN)
   100    1.1329717069             1.1329717069                 0.8826345741             1.0000000000             5.66e-04
  1000    1.0125782933             1.0125782933                 0.9875779548             1.0000000000             5.67e-06
 10000    1.0012507814             1.0012507814                 0.9987507811             1.0000000000             5.67e-08
```

この例では 2 つの固有値の絶対値が等しいので、陽的 Euler 法のノルムは <Ref to="prop-explicit-euler-unstable" /> の式 $(1 + \Delta t^2\Omega^2/4)^{n/2}$ に一致します。ノルムのずれは $\Delta t$ に比例して減るだけで、100 ステップでは確率の総和が $1.13^2 \approx 1.28$ になります。Crank–Nicolson 法はノルムを保ち、確率 $|\psi_2|^2$ の誤差はステップ数を 10 倍にするごとに $1/100$ です。
</Example>

## 7. 分割法

$H = A + B$ で、$e^{-i\tau A}$ と $e^{-i\tau B}$ はそれぞれ安く計算できるが $e^{-i\tau(A+B)}$ はそうでない、という状況がよくあります。典型は $A$ が運動エネルギー（Fourier 変換で対角化される）、$B$ がポテンシャル（位置の表示で対角）の場合です。

<Definition id="def-splitting" title="Lie–Trotter 分割と Strang 分割">
エルミート行列 $A, B$ と $\tau \gt 0$ について

$$
S_1(\tau) = e^{-i\tau A}e^{-i\tau B},\qquad S_2(\tau) = e^{-i\tau B/2}\,e^{-i\tau A}\,e^{-i\tau B/2}
$$

を **Lie–Trotter 分割**、**Strang 分割**（対称分割）という。$t = n\tau$ のとき、$S_1(\tau)^n$、$S_2(\tau)^n$ を $e^{-it(A+B)}$ の近似とする。
</Definition>

$S_1(\tau)$ と $S_2(\tau)$ は、ユニタリ行列（<Ref to="lem-unitary-group" />）の積なのでユニタリです。したがって、分割法はどんな $\tau$ でもノルムを厳密に保存します。$S_2(\tau)^n$ では隣り合う $e^{-i\tau B/2}$ がまとまって $e^{-i\tau B}$ になるので、最初と最後の半ステップを除けば、1 ステップの手間は $S_1$ と同じです。

<Lemma id="lem-conjugation-commutator" title="共役による変化と交換子">
$A$ をエルミート行列、$X$ を同じ次数の行列、$s$ を実数とすると

$$
\bigl\|e^{isA}Xe^{-isA} - X\bigr\| \le |s|\,\|[A, X]\|
$$

である。
</Lemma>

<Proof of="lem-conjugation-commutator">
$Y(s) = e^{isA}Xe^{-isA}$ とおくと、<Ref to="mathematics/ode/linear-systems#prop-exp-properties" /> の 3 から $Y'(s) = e^{isA}(iAX - iXA)e^{-isA} = i\,e^{isA}[A, X]e^{-isA}$ です。$e^{\pm isA}$ はユニタリなので $\|Y'(s)\| = \|[A, X]\|$ で、$Y(s) - X = \int_0^sY'(\sigma)\,d\sigma$ から主張が出ます。
</Proof>

次の定理は、2 つの因子を交互に掛けたときの誤差を $A$ と $B$ の交換子で評価します。$A$ と $B$ が可換なら分割は厳密で、可換からのずれの大きさが誤差になる、というのが直感です。1 は誤差を積分で直接評価し、2 は[経路積分 Monte Carlo 法](/computer-science/stochastic-computation/path-integral-monte-carlo)の章の定理に $iA$、$iB$ を代入して得ます。

<Theorem id="thm-splitting-unitary-error" title="ユニタリな分割の誤差">
$A, B$ を $d$ 次のエルミート行列、$t \ge 0$、$n \ge 1$、$\tau = t/n$ とする。

1. $\bigl\|S_1(\tau)^n - e^{-it(A+B)}\bigr\| \le \dfrac{t\,\tau}{2}\,\|[A, B]\|$。
2. $\bigl\|S_2(\tau)^n - e^{-it(A+B)}\bigr\| \le t\,\tau^2\Bigl(\dfrac{\|[A,[A,B]]\|}{12} + \dfrac{\|[B,[B,A]]\|}{24}\Bigr)$。
3. $[A, B] = 0$ なら $S_1(\tau)^n = S_2(\tau)^n = e^{-it(A+B)}$ である。
</Theorem>

<Proof of="thm-splitting-unitary-error">
1. $F(s) = e^{-isA}e^{-isB}$ とおきます。<Ref to="mathematics/ode/linear-systems#prop-exp-properties" /> の 3 から

$$
F'(s) = -iAF(s) - ie^{-isA}Be^{-isB} = -i(A + B)F(s) + R(s),\qquad R(s) = i\bigl(Be^{-isA} - e^{-isA}B\bigr)e^{-isB}
$$

です。$F(0) = I$ なので、定数変化法（<Ref to="mathematics/ode/linear-systems#thm-variation-of-constants" /> を各列に使う。基本解行列は $e^{-is(A+B)}$）から

$$
F(\tau) - e^{-i\tau(A+B)} = \int_0^\tau e^{-i(\tau - s)(A+B)}R(s)\,ds
$$

です。ユニタリ行列を掛けてもノルムは変わらないので、<Ref to="lem-conjugation-commutator" />（$X = B$）から $\|R(s)\| = \|Be^{-isA} - e^{-isA}B\| = \|e^{isA}Be^{-isA} - B\| \le s\|[A, B]\|$ で、

$$
\bigl\|S_1(\tau) - e^{-i\tau(A+B)}\bigr\| \le \int_0^\tau s\,\|[A, B]\|\,ds = \frac{\tau^2}{2}\|[A, B]\|
$$

です。$S_1(\tau)$ と $e^{-i\tau(A+B)}$ はユニタリなので、<Ref to="lem-unitary-group" /> の最後の主張から $\|S_1(\tau)^n - e^{-it(A+B)}\| \le n\tau^2\|[A,B]\|/2 = t\tau\|[A,B]\|/2$ です。

2. <Ref to="computer-science/stochastic-computation/path-integral-monte-carlo#thm-symmetric-trotter-error" /> の 1 を、その定理の $A$、$B$、$\beta$、$M$ として $iB$、$iA$、$t$、$n$ をとって使います（参照先の定理では外側に半ステップで置く行列が $A$ なので、こちらで外側に置く $B$ を向こうの $A$ に対応させます）。$iA$、$iB$、$i(A+B)$ は反エルミート行列で、$e^{-s(iA)} = e^{-isA}$ などはユニタリなのでノルムは 1 です。したがって、その定理の仮定である 3 つの半群の縮小性 $\|e^{-s(iB)}\| \le 1$、$\|e^{-s(iA)}\| \le 1$、$\|e^{-s(iA + iB)}\| \le 1$（$s \ge 0$）が満たされます（<Ref to="computer-science/stochastic-computation/path-integral-monte-carlo#lem-contraction-cases" /> の 2）。参照先の対称分割 $e^{-\tau(iB)/2}e^{-\tau(iA)}e^{-\tau(iB)/2}$ は $S_2(\tau)$ で、その $n$ 乗が $S_2(\tau)^n$ です。参照先の評価では外側の行列の二重交換子に $1/24$、内側の行列の二重交換子に $1/12$ が付きます。交換子は $[iB, [iB, iA]] = i^3[B,[B,A]]$、$[iA,[iA,iB]] = i^3[A,[A,B]]$ なので、右辺は

$$
\frac{t^3}{n^2}\Bigl(\frac{\|[B,[B,A]]\|}{24} + \frac{\|[A,[A,B]]\|}{12}\Bigr) = t\,\tau^2\Bigl(\frac{\|[A,[A,B]]\|}{12} + \frac{\|[B,[B,A]]\|}{24}\Bigr)
$$

です。

3. $A$ と $B$ が可換なら $e^{-i\tau(A+B)} = e^{-i\tau A}e^{-i\tau B}$（<Ref to="mathematics/ode/linear-systems#prop-exp-properties" /> の 1）で、$e^{-i\tau B/2}$ は $e^{-i\tau A}$ と可換（$A$ と可換な $B$ は $A$ の冪級数 $e^{-i\tau A}$ とも可換）なので $S_2(\tau) = e^{-i\tau A}e^{-i\tau B}$ です。
</Proof>

1 は参照先の定理の 2（対称でない分割。$\frac{\beta^2}{2M}\|[A,B]\|$）に $iA$、$iB$、$\beta = t$、$M = n$ を入れたものと同じ評価で、上の証明はそれをユニタリな場合に直接書いたものです。2 では、半ステップで外側に置く $B$（分割演算子法ではポテンシャル）の二重交換子 $[B,[B,A]]$ に $1/24$、内側の $A$（運動項）の二重交換子 $[A,[A,B]]$ に $1/12$ が付きます。誤差の定数が $\|A\|$、$\|B\|$ ではなく交換子で書けることが重要です。$A$ と $B$ が「ほとんど可換」なら、それぞれのノルムが大きくても誤差は小さくなります。対称分割で誤差の次数が 1 つ上がるのは、$S_2(-\tau) = S_2(\tau)^{-1}$（時間反転対称。<Ref to="prop-leapfrog-symplectic" /> の $\Phi_{\Delta t} = LWL$ と同じ構造）のため、1 ステップの誤差の $\tau^2$ の項が消えるからです。実際、速度 Verlet 法は、$q' = p$ と $p' = -Aq$ のそれぞれの厳密な流れを Strang の順に合成したものです。

<Example id="ex-splitting-bounds" title="分割の誤差と交換子による上界">
乱数で作った 6 次のエルミート行列 $A, B$ と、可換な組（$B = A^2 - 2A$）について、$t = 1$ での誤差を上界と比べます。参照は `scipy.linalg.expm` です。

```python
# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.linalg import expm

rng = np.random.default_rng(1)
def hermitian(n):
    X = rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))
    return (X + X.conj().T) / 2

comm = lambda X, Y: X @ Y - Y @ X
norm2 = lambda X: np.linalg.norm(X, 2)

def report(A, B, t=1.0):
    exact = expm(-1j * t * (A + B))
    c1 = norm2(comm(A, B))
    c2 = norm2(comm(A, comm(A, B))) / 12 + norm2(comm(B, comm(B, A))) / 24
    print("     n    Lie-Trotter   bound t tau c1/2    Strang       bound t tau^2 c2")
    for n in [10, 20, 40, 80, 160]:
        tau = t / n
        lie = np.linalg.matrix_power(expm(-1j * tau * A) @ expm(-1j * tau * B), n)
        half = expm(-0.5j * tau * B)
        strang = np.linalg.matrix_power(half @ expm(-1j * tau * A) @ half, n)
        print(f"{n:6d}    {norm2(lie - exact):.3e}     {t * tau * c1 / 2:.3e}         "
              f"{norm2(strang - exact):.3e}    {t * tau**2 * c2:.3e}")

A, B = hermitian(6), hermitian(6)
print("random Hermitian A, B (6 x 6)")
report(A, B)
print("commuting case: B = A^2 - 2A")
report(A, A @ A - 2 * A)
```

```
random Hermitian A, B (6 x 6)
     n    Lie-Trotter   bound t tau c1/2    Strang       bound t tau^2 c2
    10    2.095e-01     4.804e-01         1.345e-02    4.802e-02
    20    1.051e-01     2.402e-01         3.355e-03    1.200e-02
    40    5.268e-02     1.201e-01         8.382e-04    3.001e-03
    80    2.637e-02     6.005e-02         2.095e-04    7.503e-04
   160    1.319e-02     3.002e-02         5.238e-05    1.876e-04
commuting case: B = A^2 - 2A
     n    Lie-Trotter   bound t tau c1/2    Strang       bound t tau^2 c2
    10    2.214e-15     3.199e-16         2.919e-15    5.675e-17
    20    3.928e-15     1.599e-16         3.692e-15    1.419e-17
    40    4.076e-15     7.996e-17         6.807e-15    3.547e-18
    80    5.594e-15     3.998e-17         1.138e-14    8.866e-19
   160    2.686e-14     1.999e-17         3.370e-14    2.217e-19
```

誤差は上界の $1/2$〜$1/4$ 程度で、$n$ を 2 倍にすると Lie–Trotter は $1/2$、Strang は $1/4$ になります。可換な場合は誤差が丸めの大きさです（上界の列も、交換子を浮動小数点数で計算した丸め誤差だけが残っています。誤差が上界より大きく見えるのは、行列指数関数と冪の計算の丸め誤差です）。
</Example>

<Definition id="def-split-operator-fourier" title="分割演算子 Fourier 法">
区間 $[-L, L)$ を周期的とみなし、$N$ 点の格子 $x_j = -L + j\Delta x$（$\Delta x = 2L/N$、$j = 0, \ldots, N-1$）と波数 $\kappa_m = \pi m/L$（$m$ は $N$ 個の連続する整数。$N$ が奇数なら $|m| \le (N-1)/2$）をとる。$F$ を $F_{mj} = N^{-1/2}e^{-i\kappa_mx_j}$ で定まる行列（離散 Fourier 変換）とし、

$$
T_N = F^{*}\,\mathrm{diag}\bigl(\tfrac12\kappa_m^2\bigr)\,F,\qquad V_N = \mathrm{diag}\bigl(V(x_j)\bigr)
$$

とおく。$i\psi_t = -\frac12\psi_{xx} + V(x)\psi$ の近似として $i\psi' = (T_N + V_N)\psi$ を Strang 分割

$$
\psi^{n+1} = e^{-i\tau V_N/2}\,F^{*}\,\mathrm{diag}\bigl(e^{-i\tau\kappa_m^2/2}\bigr)\,F\,e^{-i\tau V_N/2}\,\psi^n
$$

で進める方法を**分割演算子 Fourier 法**（split-step Fourier 法）という。
</Definition>

$T_N$ は、[スペクトル法](/computer-science/pde-discretization/spectral-methods)の 2 階のスペクトル微分行列（<Ref to="computer-science/pde-discretization/spectral-methods#prop-spectral-differentiation-matrix" /> の 4。区間の長さを $2\pi$ から $2L$ に変えたもの）の $-\frac12$ 倍です。空間の離散化の誤差は、その章の補間誤差の評価に従います。

<Proposition id="prop-split-operator-fourier" title="分割演算子 Fourier 法の性質">
<Ref to="def-split-operator-fourier" /> について次が成り立つ。

1. $F$ はユニタリで、$T_N$ はエルミート、$e^{-i\tau T_N} = F^{*}\mathrm{diag}(e^{-i\tau\kappa_m^2/2})F$ である。1 ステップの行列は $S_2(\tau)$（$A = T_N$、$B = V_N$）で、ユニタリである。したがって $\|\psi^n\| = \|\psi^0\|$。
2. $t = n\tau$ について
   $$
   \bigl\|\psi^n - e^{-it(T_N + V_N)}\psi^0\bigr\| \le t\,\tau^2\Bigl(\frac{\|[T_N,[T_N,V_N]]\|}{12} + \frac{\|[V_N,[V_N,T_N]]\|}{24}\Bigr)\|\psi^0\|.
   $$
3. 連続の作用素 $T = -\frac12\frac{d^2}{dx^2}$ と $C^4$ 級の関数 $V$ について、$C^4$ 級の関数 $\psi$ に対し
   $$
   [T, V]\psi = -\tfrac12V''\psi - V'\psi',\qquad [V,[V,T]]\psi = -(V')^2\psi
   $$
   である。
</Proposition>

<Proof of="prop-split-operator-fourier">
1. $(FF^{*})_{mm'} = \frac1N\sum_je^{-i(\kappa_m - \kappa_{m'})x_j} = \frac1Ne^{i(\kappa_m - \kappa_{m'})L}\sum_{j=0}^{N-1}e^{-2\pi i(m - m')j/N}$ で、$m, m'$ は $N$ 個の連続する整数の中にあるので $m \ne m'$ なら $m - m'$ は $N$ の倍数でなく、和は 0（等比数列の和の公式）、$m = m'$ なら 1 です。よって $F$ はユニタリです。$T_N$ は実対角行列のユニタリ共役なのでエルミートで、<Ref to="lem-unitary-group" /> の証明と同じく $e^{-i\tau T_N} = F^{*}\mathrm{diag}(e^{-i\tau\kappa_m^2/2})F$ です。$V_N$ は実対角行列で $e^{-i\tau V_N/2} = \mathrm{diag}(e^{-i\tau V(x_j)/2})$ なので、1 ステップの行列は $e^{-i\tau V_N/2}e^{-i\tau T_N}e^{-i\tau V_N/2}$ です。

2. <Ref to="thm-splitting-unitary-error" /> の 2 です。

3. $T(V\psi) = -\frac12(V''\psi + 2V'\psi' + V\psi'')$、$V(T\psi) = -\frac12V\psi''$ の差が最初の式です。$[V, T] = -[T, V]$ なので $[V, T]\psi = \frac12V''\psi + V'\psi'$ で、

$$
[V,[V,T]]\psi = V\bigl(\tfrac12V''\psi + V'\psi'\bigr) - \bigl(\tfrac12V''V\psi + V'(V\psi)'\bigr) = VV'\psi' - V'(V'\psi + V\psi') = -(V')^2\psi
$$

です。
</Proof>

1 ステップの演算量は、離散 Fourier 変換 2 回（FFT で各 $O(N\log N)$）と、長さ $N$ のベクトルの成分ごとの積 3 回です。行列を作る必要も連立方程式を解く必要もありません。

2 の右辺の交換子のノルムは、$N$ を大きくすると増えます（$\|T_N\| = O(N^2)$ なので）。しかし 3 が示すように、交換子は「$V$ の導関数 × 低階の微分作用素」で、$[V,[V,T]]$ は微分を含まない掛け算です。そのため、滑らかな $V$ と滑らかな波動関数については、誤差は $N$ によらず $O(\tau^2)$ になります。これを作用素ノルムではなく $\psi$ の滑らかさ（Sobolev ノルム）で評価する定理は、非有界作用素の分割についての Jahnke–Lubich の結果（参考文献 [4]）で、証明には作用素の定義域の議論が要るので、ここでは定理環境にしません。次数のつかない収束だけなら、自己共役作用素の Trotter 積公式（<Ref to="mathematics/functional-analysis/one-parameter-groups#thm-trotter" />）と、行列の場合の Lie–Trotter の積公式（<Ref to="mathematics/lie-theory/matrix-lie-groups#prop-trotter" />）から従います。

次の例のコヒーレント状態は、形を変えずに古典的な粒子のように揺れる波束のことで、厳密解が分かっているので誤差を測るのに向いています。

<Example id="ex-coherent-state" title="調和振動子のコヒーレント状態">
$H = -\frac12\frac{d^2}{dx^2} + \frac12x^2$ の初期値 $\psi_0(x) = \pi^{-1/4}e^{-(x - x_0)^2/2}$ に対する解は

$$
\psi(x, t) = \pi^{-1/4}\exp\Bigl(-\tfrac12(x - x_0\cos t)^2 - ix_0x\sin t - \tfrac{i}{2}t + \tfrac{i}{4}x_0^2\sin2t\Bigr)
$$

で、形を変えずに古典的な軌道 $x_0\cos t$ の上を往復します（代入して確かめられます）。$x_0 = 1$、$T = 2$ とし、(a) 3 点差分の $H$ に対する Crank–Nicolson 法を、同じ行列 $H$ の厳密な時間発展と比べ、(b) 分割演算子 Fourier 法を、上の厳密解と比べます。

```python
# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.linalg import eigh_tridiagonal, solve_banded

x0 = 1.0
def coherent(x, t):
    """H = -(1/2) d^2/dx^2 + x^2/2 のコヒーレント状態（t = 0 で中心 x0 の Gauss 関数）。"""
    return np.pi ** -0.25 * np.exp(-0.5 * (x - x0 * np.cos(t)) ** 2 - 1j * x0 * x * np.sin(t)
                                   - 0.5j * t + 0.25j * x0**2 * np.sin(2 * t))

T = 2.0
# (a) Crank–Nicolson（3 点差分）。参照は同じ行列 H の厳密な時間発展（固有分解）。
L, M = 10.0, 400
h = 2 * L / M
x = -L + h * np.arange(1, M)
d, e = 1.0 / h**2 + 0.5 * x**2, -0.5 / h**2 * np.ones(M - 2)
lam, U = eigh_tridiagonal(d, e)
psi0 = coherent(x, 0.0).astype(complex)
ref = U @ (np.exp(-1j * lam * T) * (U.T @ psi0))
H = lambda v: d * v + np.concatenate(([0], e * v[:-1])) + np.concatenate((e * v[1:], [0]))
H3 = np.sqrt(h) * np.linalg.norm(H(H(H(psi0))))
norm = lambda v: np.sqrt(h) * np.linalg.norm(v)
print("Crank-Nicolson:  steps   | ||psi|| - 1 |   |<H> - <H>_0|    error      order   bound T dt^2 ||H^3 psi_0|| / 12")
prev = None
for steps in [100, 200, 400, 800]:
    dt = T / steps
    ab = np.zeros((3, M - 1), dtype=complex)           # I + i dt H / 2 の帯行列
    ab[0, 1:], ab[1], ab[2, :-1] = 0.5j * dt * e, 1 + 0.5j * dt * d, 0.5j * dt * e
    psi = psi0.copy()
    for _ in range(steps):
        psi = solve_banded((1, 1), ab, psi - 0.5j * dt * H(psi))
    err = norm(psi - ref)
    dE = abs(h * np.vdot(psi, H(psi)) - h * np.vdot(psi0, H(psi0)))
    print(f"                 {steps:5d}     {abs(norm(psi) - norm(psi0)):.1e}        {dE:.1e}       {err:.3e}   "
          + ("    " if prev is None else f"{np.log2(prev / err):.2f}") + f"    {T * dt**2 * H3 / 12:.3e}")
    prev = err

# (b) 分割演算子 Fourier 法（Strang）。参照は連続の厳密解。
N = 128
h = 2 * L / N
x = -L + h * np.arange(N)
k = 2 * np.pi * np.fft.fftfreq(N, d=h)
V = 0.5 * x**2
print("split-operator:  steps   | ||psi|| - 1 |   error      order")
prev = None
for steps in [25, 50, 100, 200, 400]:
    dt = T / steps
    half, kin = np.exp(-0.5j * dt * V), np.exp(-0.5j * dt * k**2)
    psi = coherent(x, 0.0).astype(complex)
    for _ in range(steps):
        psi = half * np.fft.ifft(kin * np.fft.fft(half * psi))
    err = np.sqrt(h) * np.linalg.norm(psi - coherent(x, T))
    print(f"                 {steps:5d}     {abs(np.sqrt(h) * np.linalg.norm(psi) - 1):.1e}        {err:.3e}   "
          + ("" if prev is None else f"{np.log2(prev / err):.2f}"))
    prev = err
```

```
Crank-Nicolson:  steps   | ||psi|| - 1 |   |<H> - <H>_0|    error      order   bound T dt^2 ||H^3 psi_0|| / 12
                   100     5.8e-15        1.4e-14       5.325e-04           5.329e-04
                   200     1.3e-15        1.6e-15       1.332e-04   2.00    1.332e-04
                   400     4.3e-15        4.1e-15       3.331e-05   2.00    3.331e-05
                   800     6.9e-15        1.6e-14       8.327e-06   2.00    8.327e-06
split-operator:  steps   | ||psi|| - 1 |   error      order
                    25     1.3e-15        1.141e-03   
                    50     1.6e-15        2.853e-04   2.00
                   100     3.8e-15        7.131e-05   2.00
                   200     6.9e-15        1.783e-05   2.00
                   400     1.6e-14        4.457e-06   2.00
```

どちらもノルムは丸めの大きさで保存され、誤差の次数は 2 です。Crank–Nicolson 法の誤差は <Ref to="thm-cayley-unitary" /> の上界とほぼ一致しています。<Ref to="lem-cayley-scalar" /> の不等式は $x \to 0$ で等号に近づき（$x - 2\arctan(x/2) = x^3/12 + O(x^5)$）、この例では $\psi_0$ の主な成分が $\Delta t\mu_k \ll 1$ の範囲にあるからです。(a) は時間の誤差だけを見ています。連続の解との差には、これに 3 点差分の空間の誤差 $O(h^2)$ が加わります。(b) は $N = 128$ 点で空間の誤差が $10^{-6}$ より十分小さいので、連続の解との差がそのまま時間の誤差になっています。同じステップ数で比べると、分割演算子法の誤差は Crank–Nicolson 法の約 $1/7$ です（100 ステップで $7.1\times10^{-5}$ と $5.3\times10^{-4}$）。
</Example>

## 8. 箱による切断

全空間の問題を計算するには、有限の区間で打ち切ります。打ち切りの誤差は、波動関数が箱の端に届くまでは小さく、届いた分だけ生じます。

次の命題は、箱の中だけで計算した解と真の解の差が、箱の内側と外側をつなぐ項の時間積分で抑えられることを示します。$P$ は箱の中の成分だけを残す[直交射影](/mathematics/functional-analysis/hilbert-spaces#def-orthogonal-projection)、$Q = I - P$ は箱の外の成分だけを残す射影です。

<Proposition id="prop-box-truncation-dynamics" title="箱による切断の誤差">
$H$ を $d$ 次のエルミート行列、$P$ を直交射影（$P = P^{*} = P^2$）、$Q = I - P$ とし、$H_P = PHP$ とおく。$\psi(t) = e^{-itH}\psi_0$、$\varphi(t) = e^{-itH_P}P\psi_0$ とすると、$t \ge 0$ について

$$
\|P\psi(t) - \varphi(t)\| \le \int_0^t\|PHQ\,\psi(s)\|\,ds,\qquad \|\psi(t) - \varphi(t)\| \le \|Q\psi(t)\| + \int_0^t\|PHQ\,\psi(s)\|\,ds
$$

が成り立つ。$H = \frac{1}{2h^2}\mathrm{tridiag}(-1,2,-1) + \mathrm{diag}(V_j)$（格子点 $j = -J', \ldots, J'$）で、$P$ が $|j| \le J$（$J \lt J'$）の成分への射影のとき、

$$
\|PHQ\,\psi(s)\| = \frac{1}{2h^2}\bigl(|\psi_{J+1}(s)|^2 + |\psi_{-J-1}(s)|^2\bigr)^{1/2}
$$

である。
</Proposition>

<Proof of="prop-box-truncation-dynamics">
$H_P$ はエルミートで、$PH_P = H_PP = H_P$ なので $e^{-itH_P}$ は $P$ と可換で、$\varphi(t) = P\varphi(t)$ です。$w(t) = P\psi(t) - \varphi(t)$ とおくと、$\psi' = -iH\psi$、$\varphi' = -iH_P\varphi$、$I = P + Q$ から

$$
w' = -iPH(P + Q)\psi + iH_P\varphi = -iH_P(P\psi - \varphi) - iPHQ\psi = -iH_Pw - iPHQ\psi,\qquad w(0) = 0
$$

です（$PHP\psi = H_PP\psi$ を使った）。定数変化法（<Ref to="mathematics/ode/linear-systems#thm-variation-of-constants" />）から $w(t) = -i\int_0^te^{-i(t-s)H_P}PHQ\psi(s)\,ds$ で、$e^{-i(t-s)H_P}$ はユニタリ（<Ref to="lem-unitary-group" />）なので最初の不等式が出ます。2 番目は $\psi - \varphi = Q\psi + w$ と三角不等式です。

三重対角の場合、$(HQ\psi)_j = \sum_lH_{jl}(Q\psi)_l$ で、$|j| \le J$ かつ $|l| \gt J$ で $H_{jl} \ne 0$ となるのは $(j, l) = (J, J+1), (-J, -J-1)$ だけで、その値は $-1/(2h^2)$ です。よって $PHQ\psi$ の 0 でない成分は $j = \pm J$ の $-\psi_{\pm(J+1)}/(2h^2)$ だけです。
</Proof>

$\varphi$ は、箱の中だけで（箱の外を 0 とする Dirichlet 境界条件で）計算した解です。命題は、その誤差が、箱のすぐ外の格子点での真の解の大きさ $|\psi_{\pm(J+1)}(s)|$ の時間積分で抑えられることを言っています。波動関数の裾が時刻 $t$ まで箱の端で小さければ、切断の誤差は小さいままです。逆に、波束が端に届くと、Dirichlet 境界でも周期境界でも反射や回り込みが起こり、誤差は $O(1)$ になります。散乱の計算で端に吸収ポテンシャルを置くのはこのためです。

<Example id="ex-gaussian-box" title="自由な Gauss 波束と箱の大きさ">
自由粒子（$V = 0$）の初期値 $\psi_0(x) = (2\pi\sigma^2)^{-1/4}e^{-x^2/(4\sigma^2) + ik_0x}$ に対する解は、確率密度が平均 $k_0t$、分散 $\sigma_t^2 = \sigma^2 + t^2/(4\sigma^2)$ の正規分布になります（Fourier 変換で $\hat\psi(\xi, t) = e^{-it\xi^2/2}\hat\psi_0(\xi)$ とし、Gauss 関数の Fourier 変換 <Ref to="mathematics/fourier-analysis/fourier-transform#prop-gaussian-ft" /> を使う。ノルムの保存は <Ref to="mathematics/fourier-analysis/fourier-transform#thm-plancherel" />）。$L \gt |k_0|t$ のとき、箱 $[-L, L]$ の外にある確率は、正規分布の裾の評価から $e^{-(L - |k_0|t)^2/(2\sigma_t^2)}$ 以下です。$\sigma = 1$、$k_0 = 2$、$t = 3$ で、周期境界の箱で（運動項だけなので時間の誤差なしに）計算した解と厳密解の差を、裾の確率の上界の平方根と比べます。

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

sigma, k0, T, steps = 1.0, 2.0, 3.0, 600
def packet(x, t):
    """自由粒子（H = -(1/2) d^2/dx^2）の Gauss 波束。t = 0 で幅 sigma、平均運動量 k0。"""
    s = sigma**2 + 0.5j * t
    return (2 * np.pi * sigma**2) ** -0.25 * np.sqrt(sigma**2 / s) * np.exp(
        -(x - k0 * t) ** 2 / (4 * s) + 1j * k0 * x - 0.5j * k0**2 * t)

print("    L     N     error at t = 3    sqrt of the tail bound exp(-(L - k0 t)^2 / (2 sigma_t^2))")
for L in [8.0, 10.0, 12.0, 16.0, 20.0]:
    N = int(16 * L)
    h = 2 * L / N
    x = -L + h * np.arange(N)
    k = 2 * np.pi * np.fft.fftfreq(N, d=h)
    psi = packet(x, 0.0)
    kin = np.exp(-0.5j * (T / steps) * k**2)
    for _ in range(steps):
        psi = np.fft.ifft(kin * np.fft.fft(psi))        # V = 0 なので運動項だけ（時間の誤差はない）
    err = np.sqrt(h) * np.linalg.norm(psi - packet(x, T))
    st2 = sigma**2 + T**2 / (4 * sigma**2)
    print(f"  {L:4.1f}  {N:4d}     {err:.3e}         {np.exp(-(L - k0 * T) ** 2 / (4 * st2)):.3e}")
```

```
    L     N     error at t = 3    sqrt of the tail bound exp(-(L - k0 t)^2 / (2 sigma_t^2))
   8.0   128     3.758e-01         7.351e-01
  10.0   160     1.203e-01         2.921e-01
  12.0   192     2.222e-02         6.271e-02
  16.0   256     1.326e-04         4.563e-04
  20.0   320     7.256e-08         2.833e-07
```

時刻 3 で波束の中心は $x = 6$、幅は $\sigma_3 \approx 1.8$ です。$L = 8$ では裾が端を越えて反対側から回り込み、誤差は $0.38$ です。$L = 20$ では $10^{-7}$ です。誤差は裾の確率の平方根（振幅の大きさ）と同じ減り方をしています。必要な箱の大きさは、初期の幅ではなく、最終時刻での位置 $|k_0|t$ と広がった幅 $\sigma_t$ で決まります。
</Example>

<Remark id="rem-time-stepping-further" title="この章で扱わなかった方法">
行列指数関数の作用 $e^{-itH}\psi$ を Krylov 部分空間で近似する方法（Hochbruck–Lubich、参考文献 [5]）、Chebyshev 多項式で展開する方法（Tal-Ezer–Kosloff）、4 次以上の分割（Suzuki、Yoshida の合成法）、吸収境界と複素吸収ポテンシャルは、Lubich のモノグラフ（参考文献 [2]）と Hairer・Lubich・Wanner の教科書（参考文献 [1]）にあります。分割の誤差を $\tau$ の冪級数として形式的に求めるには Baker–Campbell–Hausdorff の公式（<Ref to="mathematics/lie-theory/matrix-lie-groups#thm-bch" />）が使えますが、<Ref to="thm-splitting-unitary-error" /> のような $\tau$ の大きさによらない上界は、この章のように積分表示から導きます。
</Remark>

## 9. 演習

<Exercise id="exr-leapfrog-oscillator" difficulty="易">
調和振動子 $q'' = -\omega^2q$ に leapfrog 法を $\omega\Delta t = 1$ で使う。数値解の振動数 $\tilde\omega$ と、厳密な 1 周期 $2\pi/\omega$ のあいだに生じる位相の誤差を求めよ。
<Solution>
<Ref to="prop-leapfrog-symplectic" /> の 3 から $\tilde\omega = \frac{2}{\Delta t}\arcsin\frac12 = \frac{2}{\Delta t}\cdot\frac\pi6 = \frac{\pi}{3}\omega \approx 1.0472\,\omega$ です。1 周期 $t = 2\pi/\omega$ のあいだの位相の誤差は $(\tilde\omega - \omega)t = 2\pi(\pi/3 - 1) \approx 0.297$ ラジアン（約 17 度）で、数値解は進みすぎます。$\omega\Delta t = 1$ は安定条件 $\omega\Delta t \lt 2$ を十分に満たしていますが、1 周期あたり $2\pi \approx 6.3$ ステップしかなく、位相の精度は低いです。先頭項の近似 $\omega^2\Delta t^2/24 = 1/24 \approx 0.0417$ は、正確な相対誤差 $0.0472$ より少し小さい値を与えます。
</Solution>
</Exercise>

<Exercise id="exr-cn-energy" difficulty="標準">
Crank–Nicolson 法について、$\psi^{n+1} - \psi^n$ と $\psi^{n+1} + \psi^n$ の内積を使って、固有分解を使わずに $\|\psi^{n+1}\| = \|\psi^n\|$ を示せ。同じ方法で、陰的 Euler 法では $\|\psi^{n+1}\|^2 = \|\psi^n\|^2 - \|\psi^{n+1} - \psi^n\|^2$ となることを示せ。
<Solution>
Crank–Nicolson 法の式は $\psi^{n+1} - \psi^n = -\frac{i\Delta t}{2}H(\psi^{n+1} + \psi^n)$ です。$\phi = \psi^{n+1} + \psi^n$ との内積をとると $\langle\phi, \psi^{n+1} - \psi^n\rangle = -\frac{i\Delta t}{2}\langle\phi, H\phi\rangle$ で、$H$ がエルミートなので $\langle\phi, H\phi\rangle$ は実数、右辺は純虚数です。左辺の実部は $\mathrm{Re}\langle\psi^{n+1} + \psi^n, \psi^{n+1} - \psi^n\rangle = \|\psi^{n+1}\|^2 - \|\psi^n\|^2$（交差項 $\langle\psi^n, \psi^{n+1}\rangle - \langle\psi^{n+1}, \psi^n\rangle$ は純虚数）なので、これが 0 です。

陰的 Euler 法は $\psi^{n+1} - \psi^n = -i\Delta tH\psi^{n+1}$ です。$\psi^{n+1}$ との内積の実部をとると $\mathrm{Re}\langle\psi^{n+1}, \psi^{n+1} - \psi^n\rangle = 0$ です。恒等式 $2\mathrm{Re}\langle a, a - b\rangle = \|a\|^2 - \|b\|^2 + \|a - b\|^2$ から $\|\psi^{n+1}\|^2 - \|\psi^n\|^2 + \|\psi^{n+1} - \psi^n\|^2 = 0$ です。このエネルギー法の議論は $H$ が変数係数でも、時間に依存しても通ります（<Ref to="computer-science/pde-discretization/elliptic-parabolic-difference#prop-energy-estimate" /> と同じ考え方です）。
</Solution>
</Exercise>

<Exercise id="exr-strang-oscillator" difficulty="標準">
$T = -\frac12\frac{d^2}{dx^2}$、$V = \frac12x^2$ について、滑らかな関数に作用する作用素として $[T,[T,V]] = -2T$、$[V,[V,T]] = -2V$ を示せ。これから、調和振動子の分割演算子法の誤差の定数が、作用素のノルムではなく状態のエネルギーの大きさで決まることを説明せよ。
<Solution>
<Ref to="prop-split-operator-fourier" /> の 3 で $V' = x$、$V'' = 1$ とすると $[T, V]\psi = -\frac12\psi - x\psi'$、$[V,[V,T]]\psi = -x^2\psi = -2V\psi$ です。次に $[T, x\frac{d}{dx}]$ を計算します。$\frac{d^2}{dx^2}(x\psi') = 2\psi'' + x\psi'''$、$x\frac{d}{dx}\psi'' = x\psi'''$ なので $[\frac{d^2}{dx^2}, x\frac{d}{dx}] = 2\frac{d^2}{dx^2}$、$[T, x\frac{d}{dx}] = -\frac{d^2}{dx^2} = 2T$ です。定数 $-\frac12$ との交換子は 0 なので $[T,[T,V]] = -[T, x\frac{d}{dx}] = -2T$ です。

したがって <Ref to="thm-splitting-unitary-error" /> の 2 の形の評価で交換子を状態 $\psi$ に作用させたものは $\|T\psi\|$ と $\|V\psi\|$ の大きさで、どちらも $\psi$ の運動エネルギーと位置エネルギーの程度です。<Ref to="ex-coherent-state" /> のコヒーレント状態ではこれらは 1 の程度なので、誤差は $t\tau^2 \times O(1)$ で、格子を細かくしても増えません。行列 $T_N$ のノルムは $O(N^2)$ ですが、それは波動関数がほとんど成分をもたない高い波数から来ています。
</Solution>
</Exercise>

<Exercise id="exr-wave-two-step" difficulty="難">
<Ref to="thm-leapfrog-energy" /> の離散エネルギーを、<Ref to="prop-leapfrog-symplectic" /> の $(q^n, p^n)$ を使って $E^n = \frac12\|p^{n+1/2}\|^2 + \frac12\langle q^{n+1}, Aq^n\rangle$ と書く。$A$ が $\Delta t^2\lambda_d \lt 4$ を満たすとき、$\tilde E(q, p) = \frac12\|p\|^2 + \frac12\langle q, A(I - \frac{\Delta t^2}{4}A)q\rangle$ が $\Phi_{\Delta t}$ の不変量であることを、$A$ の固有分解で <Ref to="prop-leapfrog-symplectic" /> の 3 に帰着させて示せ。また $\tilde E(q^n, p^n) = E^n$ を確かめよ。
<Solution>
$A$ の正規直交固有基底 $q_k$（固有値 $\lambda_k = \omega_k^2$）で $q = \sum_ka_kq_k$、$p = \sum_kb_kq_k$ と展開すると、$\Phi_{\Delta t}$ は $A$ とスカラー倍だけからなるので、各 $k$ の成分 $(a_k, b_k)$ に、$d = 1$、$A = \omega_k^2$ の場合の写像として作用します。$\tilde E = \sum_k\frac12\bigl(b_k^2 + \omega_k^2(1 - \omega_k^2\Delta t^2/4)a_k^2\bigr)$ で、各項は <Ref to="prop-leapfrog-symplectic" /> の 3 の不変量の $1/2$ 倍です（その証明の計算は $\omega\Delta t \lt 2$ を使っていないので、$\kappa$ の符号によらず成り立ちます）。

$E^n$ との一致: $p^{n+1/2} = p^n - \frac{\Delta t}{2}Aq^n$、$q^{n+1} = q^n + \Delta t\,p^{n+1/2} = (I - \frac{\Delta t^2}{2}A)q^n + \Delta t\,p^n$ を代入すると

$$
\tfrac12\|p^{n+1/2}\|^2 = \tfrac12\|p^n\|^2 - \tfrac{\Delta t}{2}\langle p^n, Aq^n\rangle + \tfrac{\Delta t^2}{8}\|Aq^n\|^2,\qquad \tfrac12\langle q^{n+1}, Aq^n\rangle = \tfrac12\langle q^n, Aq^n\rangle - \tfrac{\Delta t^2}{4}\|Aq^n\|^2 + \tfrac{\Delta t}{2}\langle p^n, Aq^n\rangle
$$

で、和は $\frac12\|p^n\|^2 + \frac12\langle q^n, Aq^n\rangle - \frac{\Delta t^2}{8}\|Aq^n\|^2 = \tilde E(q^n, p^n)$ です。$\Delta t^2\lambda_d \lt 4$ は $A(I - \frac{\Delta t^2}{4}A)$ が正定値である条件で、<Ref to="thm-leapfrog-energy" /> の 2 の正定値性の条件と同じです。離散エネルギーは、時刻 $t_n$ の位置と運動量で書くと、ポテンシャルを $O(\Delta t^2)$ だけ修正したエネルギーになっています。
</Solution>
</Exercise>

## 参考文献

1. E. Hairer, C. Lubich, G. Wanner, *Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations*, 2nd ed., Springer, 2006 — Störmer–Verlet 法、シンプレクティック性、後退誤差解析（第 IX 章）。
2. C. Lubich, *From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis*, European Mathematical Society, 2008 — Schrödinger 方程式の時間積分（分割法・Krylov 法・Chebyshev 法）。
3. R. J. LeVeque, *Finite Difference Methods for Ordinary and Partial Differential Equations*, SIAM, 2007 — 波動方程式の差分法、数値分散。
4. T. Jahnke, C. Lubich, "Error bounds for exponential operator splittings", *BIT Numerical Mathematics* 40 (2000), 735–744.
5. M. Hochbruck, C. Lubich, "On Krylov subspace approximations to the matrix exponential operator", *SIAM Journal on Numerical Analysis* 34 (1997), 1911–1925.
6. G. Strang, "On the construction and comparison of difference schemes", *SIAM Journal on Numerical Analysis* 5 (1968), 506–517.