# 検算と検証：誤差の次数・残差・保存量による検査と、数値実験と証明の区別

> 数値計算の結果を確かめる方法を、観測次数の収束、残差と誤差の関係、保存量のずれ、検算の独立性、反例探索で言える範囲の定理として整理し、数値的な証拠と証明を区別する。誤りを仕込んだ実装が各検査でどう見えるかを実行例で示す。
> https://rikai.mugen-giken.com/computer-science/research-computing/verification-of-computations

## 0. この記事の要点

プログラムが出した数値が正しいかどうかは、答えを知らない問題ほど確かめにくいものです。この記事では、答えを知らないまま計算結果を疑うための手がかりを順に見ます。刻みを半分にしたときの誤差の減り方、方程式に代入し直したときのずれ、変わらないはずの量の変化などです。読み終えると、どの検査がどんな誤りを見つけ、どんな誤りを見逃すかを説明でき、「数値で確かめた」と「証明した」を書き分けられるようになります。

- 誤差が $e(h) = Ch^p + O(h^q)$（$q \gt p$）なら、**観測次数** $\log_2\bigl(e(h)/e(h/2)\bigr)$ は $p$ に収束します。真値が分からなくても、3 つの刻みの値の差の比から同じ $p$ が得られます。実装の誤りは、多くの場合、観測次数の低下として現れます。
- **残差が小さいことは、誤差が小さいことを意味しません。** 相対誤差は条件数 × 相対残差まで大きくなりえて、等号が成り立つ例があります。
- 保存量のずれは、$p$ 次の方法で $T$ までに $O(h^p)$ 以下です。対称性や恒等式による検査は、その対称性を保つ誤りを原理的に検出できません。
- 2 つの検算が同じ誤りを見逃す確率は、独立なら積ですが、共通の原因があると積よりずっと大きくなります。「別手法」は、導出からして別でなければなりません。
- 無作為な $N$ 個の入力で反例が出なかったとき言えるのは「反例の割合は、信頼水準 $95\%$ で $3/N$ 以下」までで、命題が真であることではありません。
- 浮動小数の計算は証拠であって証明ではありません。倍精度で正に見えた量が、実は $0$ のことがあります。数値を証明に使うときは、包含の保証がある計算（区間演算・球算術）に置き換えます。

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

- **条件数**: 入力のわずかなずれが、答えのずれをどれだけ拡大するかを表す数。→ [精度の選び方](/computer-science/research-computing/precision-and-error-budget)
- **単位丸め**: 浮動小数で 1 回の演算が生む相対誤差の上限（倍精度で約 $10^{-16}$）。→ [精度の選び方](/computer-science/research-computing/precision-and-error-budget)
- **作業精度**: 計算の途中で使う桁数や精度のこと。→ [作業精度](/computer-science/research-computing/precision-and-error-budget#def-precision-terms)
- **誤差予算**: 目標の精度のために、各段階の誤差にどれだけ割り当てるかの見積もり。→ [誤差予算](/computer-science/research-computing/precision-and-error-budget#def-error-budget)
- **Richardson 外挿**: 刻みを変えた 2 つの結果から、誤差の主要な項を消して精度を上げる方法。→ [数値積分 II](/mathematics/numerical-analysis/extrapolation-and-romberg)

</Aside>

## 1. 動機

数値計算のプログラムは、誤っていても数を返します。構文の誤りは実行前に、型の誤りは実行時に見つかりますが、重みを 1 つ取り違えた求積公式も、添字を入れ替えた行列積も、何事もなかったように「それらしい数」を出力します。出力がそれらしいかどうかを眺めて判断するのは、答えを知っている問題でしか通用しません。そして、答えを知っている問題は、わざわざ計算する必要のない問題です。

そこで、答えを知らない問題の計算結果を、**答えを知らないまま**確かめる方法が要ります。原理は単純で、正しい計算ならば必ず成り立つ関係を探し、それが成り立っているかを調べます。刻みを半分にすれば誤差は $2^p$ 分の 1 になるはずだ、方程式に代入すれば残差は小さいはずだ、エネルギーは保存されるはずだ、別の方法で計算しても同じ値になるはずだ、という具合です。

ただし、これらの検査はどれも必要条件の検査です。合格しても正しいとは限りません。この章では、それぞれの検査について「何が保証され、何が保証されないか」を命題として述べます。どの検査がどの種類の誤りを捕まえ、どの種類を見逃すかを知っていれば、検査を組み合わせて穴をふさぐことができます。最後に、こうして集めた数値的な証拠と、数学の証明との間にある段差をはっきりさせ、その段差を越える方法（精度保証つき計算）への入口を示します。このテーマの最後の章として、[再現性](/computer-science/research-computing/reproducible-experiments)と[精度の選び方](/computer-science/research-computing/precision-and-error-budget)の章の道具もここで使います。

## 2. 何を確かめるのか

<Definition id="def-verification-validation" title="検証と妥当性確認">
数値計算の正しさには 2 つの別の問いがある。

- **検証**（verification）: プログラムは、意図した数式（方程式、算法、公式）を正しく計算しているか。
- **妥当性確認**（validation）: その数式は、調べたい対象（物理現象、データ）を正しく表しているか。

検証は数学とプログラムの間の問いで、対象についての知識がなくても行える。妥当性確認は数式と対象の間の問いで、実験や観測との比較を要する。
</Definition>

この章が扱うのは検証です。検証が済んでいない計算で実験との一致・不一致を論じても、モデルの誤りと実装の誤りを区別できません。検証のための検査（**検算**）を、使う情報の種類で分けると次の表になります。

| 検算の種類 | 調べること | 捕まえやすい誤り | 見逃しやすい誤り |
|---|---|---|---|
| 閉形式・既知値 | 答えの分かる問題での一致 | 公式・定数の誤り | その問題で効かない項の誤り |
| 別手法 | 独立な算法との一致 | 算法ごとの実装の誤り | 共通の入力・共通の公式の誤り |
| 別精度 | 作業精度を上げた結果との一致 | 丸め誤差・桁落ち | 打ち切り誤差、公式の誤り |
| 誤差の次数 | 刻みを変えたときの収束の速さ | 重み・添字・境界の誤り | 次数を保つ誤り |
| 残差 | 解を方程式に戻したときのずれ | 解く手続きの誤り | 方程式そのものの誤り、悪条件 |
| 極限・対称性・保存量・恒等式 | 成り立つはずの関係 | その関係を破る誤り | その関係を保つ誤り |

右端の列がこの章の主題です。どの検査にも見逃す誤りがあり、それは検査の原理から決まっています。

<Definition id="def-evidence-levels" title="数値的な証拠・反例・証明">
命題 $P$ についての数値計算の結果を、次の 3 つに区別する。

- **数値的な証拠**: 有限個の入力について、浮動小数の計算で $P$ と矛盾しない結果を得たこと。$P$ が真であることを含意しない。
- **反例**: $P$ が偽であることを示す具体的な入力。ただし浮動小数の計算で見つけた「反例」は、丸め誤差の影響を除く検査（厳密な計算または包含の保証がある計算）を経て初めて反例になる。
- **証明**: すべての入力について $P$ が成り立つことの論証。計算機を使う場合は、計算の結果が真の値を含むことが保証された計算（**精度保証つき計算**。区間演算・球算術・厳密な有理数演算）と、有限個の場合への帰着の論証からなる。
</Definition>

## 3. 誤差の次数を観測する

### 3.1. 観測次数

刻み $h$ をもつ近似（数値積分、差分、常微分方程式の数値解）の誤差は、多くの場合 $h$ のべきで減ります。そのべきを数値実験から読み取ります。

<Definition id="def-observed-order" title="観測次数">
刻み $h$ と $h/2$ での誤差 $e(h), e(h/2)$ がともに $0$ でなく同符号のとき、
$$
\hat p(h) = \log_2\frac{e(h)}{e(h/2)}
$$
を**観測次数**という。
</Definition>

証明で使う対数の評価を先に示します。

<Lemma id="lem-log-bound" title="対数の評価">
$|t| \le 1/2$ ならば $|\ln(1 + t)| \le 2|t|$。
</Lemma>

<Proof of="lem-log-bound">
$t = 0$ なら両辺は $0$ です。$t \ne 0$ のとき、関数 $s \mapsto \ln(1 + s)$ は $[-1/2, 1/2]$ で微分可能で導関数は $1/(1 + s)$ なので、$0$ と $t$ を端点とする閉区間に平均値の定理（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-mvt" />）を使うと、$0$ と $t$ の間の $\xi$ があって $\ln(1 + t) - \ln 1 = \frac{t}{1 + \xi}$ です。$|\xi| \lt 1/2$ なので $1 + \xi \gt 1/2$ で、$|\ln(1 + t)| = |t|/(1 + \xi) \le 2|t|$ です。
</Proof>

<Theorem id="thm-observed-order" title="観測次数の収束">
$p \gt 0$、$C \ne 0$ とし、ある $h \gt 0$ について
$$
e(h) = Ch^p + r(h), \qquad e(h/2) = C\Bigl(\frac h2\Bigr)^p + r(h/2)
$$
と書けるとする。$\varepsilon = \dfrac{\max\bigl(|r(h)|,\ 2^p|r(h/2)|\bigr)}{|C|h^p}$ とおく。$\varepsilon \le 1/2$ ならば、$e(h)$ と $e(h/2)$ は $0$ でなく同符号で、
$$
\bigl|\hat p(h) - p\bigr| \le \frac{4}{\ln 2}\,\varepsilon
$$
が成り立つ。とくに、$q \gt p$ と $K \ge 0$ があって、すべての $0 \lt h \le h_0$ で $|r(h)| \le Kh^q$ ならば、$Kh^{q-p} \le |C|/2$ を満たす $h \le h_0$ について
$$
\bigl|\hat p(h) - p\bigr| \le \frac{4K}{|C|\ln 2}\,h^{q-p}
$$
で、$h \to 0$ のとき $\hat p(h) \to p$ である。
</Theorem>

<Proof of="thm-observed-order">
$\alpha = r(h)/(Ch^p)$、$\beta = 2^pr(h/2)/(Ch^p)$ とおくと、$|\alpha| \le \varepsilon$、$|\beta| \le \varepsilon$ で、
$$
e(h) = Ch^p(1 + \alpha), \qquad e(h/2) = \frac{Ch^p}{2^p}(1 + \beta)
$$
です。$\varepsilon \le 1/2$ なら $1 + \alpha \ge 1/2 \gt 0$、$1 + \beta \ge 1/2 \gt 0$ なので、$e(h)$ と $e(h/2)$ はどちらも $0$ でなく、$C$ と同じ符号です。比をとると $e(h)/e(h/2) = 2^p(1 + \alpha)/(1 + \beta) \gt 0$ で、
$$
\hat p(h) = p + \frac{\ln(1 + \alpha) - \ln(1 + \beta)}{\ln 2}
$$
です。<Ref to="lem-log-bound" /> から $|\ln(1 + \alpha) - \ln(1 + \beta)| \le 2|\alpha| + 2|\beta| \le 4\varepsilon$ で、最初の評価を得ます。

後半は、$|r(h)| \le Kh^q$ と $2^p|r(h/2)| \le 2^pK(h/2)^q = 2^{p-q}Kh^q \le Kh^q$（$q \gt p$ から $2^{p-q} \lt 1$）から $\varepsilon \le Kh^{q-p}/|C|$ で、仮定 $Kh^{q-p} \le |C|/2$ のもとで $\varepsilon \le 1/2$ です。前半に代入して評価を得ます。$q \gt p$ なので右辺は $h \to 0$ で $0$ に収束します。
</Proof>

定理の前半は、$r$ の中身を問わない形で書いてあります。$r$ に打ち切り誤差の高次の項だけでなく、丸め誤差も含められるようにするためです（§3.3）。

観測次数による検査の手順は次のとおりです。理論から期待される次数 $p$ を先に書き出す。刻みを半分ずつにして 5〜6 段階で誤差を測る。観測次数が $p$ に近づけば合格、別の値に落ち着くか、どこにも落ち着かなければ、実装か仮定のどちらかが誤っています。

<Example id="ex-simpson-orders" title="正しい Simpson 則、重みを 1 つ誤った Simpson 則、滑らかでない被積分関数">
複合 Simpson 則の誤差は、被積分関数が $C^4$ 級なら $O(h^4)$ です（<Ref to="mathematics/numerical-analysis/newton-cotes-and-gauss-quadrature#cor-composite-error" /> の 2）。正しい実装と、重みを設定するスライスの終わりを 1 つ取り違えて終点の重みを $2$ にしてしまった実装を、$\int_0^1 e^x\,dx = e - 1$ で比べます。あわせて、正しい実装で $\int_0^1\sqrt x\,dx = 2/3$ を計算します。

```python
# Python 3.12.14 / NumPy 2.5.2
import math
import numpy as np

def simpson(f, a, b, N, buggy=False):
    """複合 Simpson 則（N は偶数）。重みは 1, 4, 2, 4, ..., 2, 4, 1"""
    x = np.linspace(a, b, N + 1)
    w = np.ones(N + 1)
    w[1:-1:2] = 4
    if buggy:
        w[2::2] = 2          # 誤り: 範囲の終わりを 1 つ取り違え、終点の重みまで 2 にする
    else:
        w[2:-1:2] = 2
    return (b - a) / N / 3 * np.dot(w, f(x))

cases = [("正しい実装, e^x", np.exp, math.e - 1, False),
         ("誤った実装, e^x", np.exp, math.e - 1, True),
         ("正しい実装, √x", np.sqrt, 2 / 3, False)]
Ns = [8, 16, 32, 64, 128, 256]
for name, f, exact, buggy in cases:
    errs = [abs(simpson(f, 0.0, 1.0, N, buggy) - exact) for N in Ns]
    orders = [math.log2(errs[i] / errs[i + 1]) for i in range(len(Ns) - 1)]
    print(name)
    print("   N      " + "".join(f"{N:>10d}" for N in Ns))
    print("   誤差   " + "".join(f"{e:10.2e}" for e in errs))
    print("   観測次数         " + "".join(f"{o:10.3f}" for o in orders))
```

```text
正しい実装, e^x
   N               8        16        32        64       128       256
   誤差     2.33e-06  1.46e-07  9.10e-09  5.69e-10  3.56e-11  2.22e-12
   観測次数              3.998     3.999     4.000     4.000     4.000
誤った実装, e^x
   N               8        16        32        64       128       256
   誤差     1.13e-01  5.66e-02  2.83e-02  1.42e-02  7.08e-03  3.54e-03
   観測次数              1.000     1.000     1.000     1.000     1.000
正しい実装, √x
   N               8        16        32        64       128       256
   誤差     3.59e-03  1.27e-03  4.48e-04  1.59e-04  5.61e-05  1.98e-05
   観測次数              1.500     1.500     1.500     1.500     1.500
```

**検算.** (a) 正しい実装の誤差は、評価 $\frac{(b-a)h^4}{180}\max|f^{(4)}|$ の内側にあります。$N = 8$ で $(1/8)^4 \times e/180 \approx 3.7 \times 10^{-6}$、実際は $2.33 \times 10^{-6}$ です。観測次数は $4$ に収束しています（<Ref to="thm-observed-order" />）。(b) 誤った実装は、終点の重みが $1$ 多いので、正しい値に $\frac h3f(1)$ が加わります。$N = 8$ で $e/24 \approx 0.113$ で、表の値と一致します。誤差が $h$ に比例するので観測次数は $1$ です。**誤った実装も $N$ を増やせば正しい値に収束する**ことに注意してください。「刻みを細かくしたら値が落ち着いた」ことは、実装が正しいことの証拠になりません。次数まで見て初めて誤りが分かります。(c) $\sqrt x$ の観測次数 $1.5$ は実装の誤りではありません。$\sqrt x$ は $x = 0$ で微分可能でなく、$C^4$ 級という仮定が破れています。次数が期待と違うとき、疑う先は実装と仮定の 2 つあります。
</Example>

### 3.2. 真値が分からないとき

実際の問題では真値 $A$ が分からないので、誤差 $e(h) = A(h) - A$ そのものは測れません。それでも、3 つの刻みの値があれば次数は観測できます。

次の定理のねらいは、真値 $A$ が出てこない形で次数 $p$ を取り出すことです。隣り合う刻みの値の差 $A(h) - A(h/2)$ では $A$ が打ち消し合うので、差の比 $\bigl(A(h) - A(h/2)\bigr)/\bigl(A(h/2) - A(h/4)\bigr)$ が $2^p$ に近づきます。長い条件式は、高次の項が先頭の項に比べて十分小さい、という仮定を書き下したものです。

<Theorem id="thm-three-level-order" title="3 つの刻みによる次数の推定と誤差の推定">
$0 \lt p \lt q$、$C \ne 0$、$K \ge 0$ とし、$A(h)$ が $h' \in \{h, h/2, h/4\}$ について
$$
A(h') = A + Ch'^{\,p} + r(h'), \qquad |r(h')| \le Kh'^{\,q}
$$
を満たすとする。$D(h) = A(h) - A(h/2)$ とおく。

1. $\dfrac{2Kh^{q-p}}{|C|(1 - 2^{-p})} \le \dfrac12$ ならば、$D(h)$ と $D(h/2)$ は $0$ でなく同符号で、
   $$
   \Bigl|\log_2\frac{A(h) - A(h/2)}{A(h/2) - A(h/4)} - p\Bigr| \le \frac{8K}{|C|(1 - 2^{-p})\ln 2}\,h^{q-p}.
   $$
2. 最も細かい刻みの値の誤差は、差から
   $$
   \Bigl|\bigl(A - A(h/4)\bigr) - \frac{A(h/4) - A(h/2)}{2^p - 1}\Bigr| \le \Bigl(2^{-q} + \frac{1 + 2^{-q}}{2^p - 1}\Bigr)K\Bigl(\frac h2\Bigr)^q
   $$
   の精度で推定できる。
</Theorem>

<Proof of="thm-three-level-order">
1. $h' \in \{h, h/2\}$ について
   $$
   D(h') = Ch'^{\,p} - C\Bigl(\frac{h'}{2}\Bigr)^p + r(h') - r(h'/2) = C'h'^{\,p} + \rho(h'), \qquad C' = C(1 - 2^{-p}),\quad \rho(h') = r(h') - r(h'/2)
   $$
   です。$p \gt 0$ なので $C' \ne 0$ で、$|\rho(h')| \le Kh'^{\,q} + K(h'/2)^q \le 2Kh'^{\,q}$ です。したがって $|\rho(h)| \le 2Kh^q$、$2^p|\rho(h/2)| \le 2^p \cdot 2K(h/2)^q = 2^{p-q}\cdot2Kh^q \le 2Kh^q$ で、<Ref to="thm-observed-order" /> の $\varepsilon$（$e$ を $D$、$C$ を $C'$、$r$ を $\rho$ としたもの）は
   $$
   \varepsilon \le \frac{2Kh^q}{|C'|h^p} = \frac{2Kh^{q-p}}{|C|(1 - 2^{-p})}
   $$
   を満たします。仮定からこれは $1/2$ 以下なので、同じ定理の前半から、$D(h)$ と $D(h/2)$ は $0$ でなく同符号で、$|\log_2(D(h)/D(h/2)) - p| \le \frac{4}{\ln 2}\varepsilon$ です。$\varepsilon$ の評価を代入して主張を得ます。

2. 仮定は、刻み $h/2$ と $h/4$ について Richardson 外挿の定理（<Ref to="mathematics/numerical-analysis/extrapolation-and-romberg#thm-richardson" />）の仮定（定数 $c = C$、剰余の定数 $K$）を満たしています。2 つの刻みによる誤差の推定（<Ref to="mathematics/numerical-analysis/extrapolation-and-romberg#prop-error-estimate" />）で $h$ を $h/2$ に置き換えると、主張の不等式そのものです。
</Proof>

1 で求めた次数の推定値を 2 の式の $p$ に使えば、真値も理論上の次数も知らずに、誤差の推定値まで得られます。推定した次数が理論と合っていれば実装の検証になり、誤差の推定値は結果に付ける誤差の見積もりになります。

<Example id="ex-three-level" title="3 つの刻みからの次数と誤差の推定、丸めによる次数の崩れ">
$\int_0^1 e^{-x^2}dx$ を台形則と Simpson 則で計算し、3 つの刻みの値だけから次数と誤差を推定します。真値 $\frac{\sqrt\pi}{2}\mathrm{erf}(1)$ は推定には使わず、推定の良し悪しを確かめるためだけに使います。後半では、中心差分で同じ推定を行い、刻みを小さくしすぎると何が起こるかを見ます。

```python
# Python 3.12.14 / NumPy 2.5.2
import math
import numpy as np

f = lambda x: np.exp(-x * x)
exact = math.sqrt(math.pi) / 2 * math.erf(1.0)      # 比較用の真値（推定には使わない）

def trapezoid(N):
    x = np.linspace(0.0, 1.0, N + 1)
    y = f(x)
    return (y.sum() - (y[0] + y[-1]) / 2) / N

def simpson(N):
    x = np.linspace(0.0, 1.0, N + 1)
    w = np.ones(N + 1); w[1:-1:2] = 4; w[2:-1:2] = 2
    return np.dot(w, f(x)) / (3 * N)

for name, rule in (("台形則", trapezoid), ("Simpson 則", simpson)):
    print(name)
    print("    N(h/4)  推定次数   誤差の推定    真の誤差     比")
    for N in (8, 16, 32, 64, 128):
        A1, A2, A4 = rule(N), rule(2 * N), rule(4 * N)
        p = math.log2((A1 - A2) / (A2 - A4))
        est = (A4 - A2) / (2 ** p - 1)               # A - A(h/4) の推定
        true = exact - A4
        print(f"  {4 * N:6d}   {p:8.4f}   {est:11.3e}  {true:11.3e}  {est / true:6.3f}")

print("中心差分 (f = exp, x = 1)")
print("   h=2^-k   D(h) の誤差    3 点からの推定次数")
D = lambda h: (math.exp(1 + h) - math.exp(1 - h)) / (2 * h)
for k in range(2, 26, 2):
    h = 2.0 ** -k
    num, den = D(h) - D(h / 2), D(h / 2) - D(h / 4)
    p = f"{math.log2(abs(num / den)):8.3f}" if num != 0 and den != 0 else "  （差が 0）"
    print(f"   {k:4d}    {abs(D(h) - math.e):10.2e}   {p}")
```

```text
台形則
    N(h/4)  推定次数   誤差の推定    真の誤差     比
      32     2.0007     5.985e-05    5.988e-05   0.999
      64     2.0002     1.497e-05    1.497e-05   1.000
     128     2.0000     3.742e-06    3.742e-06   1.000
     256     2.0000     9.356e-07    9.356e-07   1.000
     512     2.0000     2.339e-07    2.339e-07   1.000
Simpson 則
    N(h/4)  推定次数   誤差の推定    真の誤差     比
      32     3.9952    -7.816e-09   -7.795e-09   1.003
      64     3.9989    -4.875e-10   -4.872e-10   1.001
     128     3.9997    -3.046e-11   -3.045e-11   1.000
     256     3.9999    -1.904e-12   -1.903e-12   1.000
     512     4.0001    -1.190e-13   -1.191e-13   0.999
中心差分 (f = exp, x = 1)
   h=2^-k   D(h) の誤差    3 点からの推定次数
      2      2.84e-02      2.004
      4      1.77e-03      2.000
      6      1.11e-04      2.000
      8      6.91e-06      2.000
     10      4.32e-07      2.000
     12      2.70e-08      2.000
     14      1.69e-09      2.021
     16      1.06e-10      1.585
     18      9.94e-12     （差が 0）
     20      1.26e-10     （差が 0）
     22      3.39e-10     （差が 0）
     24      5.92e-10     （差が 0）
```

**検算.** (a) 推定次数は台形則で $2$、Simpson 則で $4$ に収束し、理論と合います（<Ref to="thm-three-level-order" /> の 1）。(b) 誤差の推定値は、符号も含めて真の誤差と 3 桁一致します（同じ定理の 2）。真値を使わずに「この値の誤差は $-1.2 \times 10^{-13}$ 程度」と言えたことになります。(c) 中心差分では、$h \ge 2^{-12}$ で推定次数が $2$ ですが、$h = 2^{-16}$ で崩れ、$h \le 2^{-18}$ では隣り合う刻みの値が完全に一致して比が定義できません。誤差も $h = 2^{-18}$ の $10^{-11}$ を底にして、それより細かい刻みでは増えています。次の小節で説明します。
</Example>

### 3.3. 丸めが次数を壊す領域

<Corollary id="cor-order-noise-floor" title="丸めの雑音があるときの観測次数">
<Ref to="thm-observed-order" /> の $e(h)$ の代わりに、計算で得た値 $\hat e(h) = Ch^p + r(h) + \delta(h)$ を使うとする。$|r(h')| \le Kh'^{\,q}$（$q \gt p$）、$|\delta(h')| \le \eta$ が $h' \in \{h, h/2\}$ で成り立ち、
$$
\varepsilon_\eta = \frac{Kh^{q-p}}{|C|} + \frac{2^p\eta}{|C|h^p} \le \frac12
$$
ならば、$\hat e$ から計算した観測次数 $\hat p$ は $|\hat p - p| \le \frac{4}{\ln 2}\varepsilon_\eta$ を満たす。
</Corollary>

<Proof of="cor-order-noise-floor">
$r(h') + \delta(h')$ を改めて $r(h')$ とみなして <Ref to="thm-observed-order" /> の前半を使います。$|r(h) + \delta(h)| \le Kh^q + \eta$、$2^p|r(h/2) + \delta(h/2)| \le 2^{p-q}Kh^q + 2^p\eta \le Kh^q + 2^p\eta$ なので、同じ定理の $\varepsilon$ は $(Kh^q + 2^p\eta)/(|C|h^p) = \varepsilon_\eta$ 以下です。
</Proof>

$\varepsilon_\eta$ の第 1 項は $h$ を小さくすると減り、第 2 項は増えます。観測次数が信頼できるのは、両方が小さい**中間の領域**だけです。$h$ が大きすぎると高次の項が邪魔をし（漸近領域、つまり先頭の項だけで誤差が決まる範囲に入っていない）、小さすぎると丸めの雑音 $\eta$ が打ち切り誤差 $|C|h^p$ を上回ります。雑音の大きさは、単位丸め $u$ に計算の条件数を掛けた程度で、これが誤差の下限になります。

中心差分はその典型です。打ち切り誤差は $\frac{M_3}{6}h^2$、丸めによる誤差は $\varepsilon M_0/h$ の程度で（<Ref to="mathematics/numerical-analysis/extrapolation-and-romberg#prop-difference-quotient-error" />）、$\eta$ 自体が $1/h$ で増えます。上の例では $\varepsilon \approx u \approx 10^{-16}$、$M_0 \approx e$ なので、$h = 2^{-18} \approx 4 \times 10^{-6}$ で丸めの項が $10^{-10}$ の程度になり、打ち切り誤差 $\frac{e}{6}h^2 \approx 7 \times 10^{-12}$ を上回ります。表で次数が崩れ始める位置と合っています。次数を観測する実験は、丸めが見え始める手前で止めるか、[作業精度を上げて](/computer-science/research-computing/precision-and-error-budget)行います。

中心差分で $h$ を小さくしていくと、誤差と観測次数がどう動くかを見てください。

<Figure caption="$f = e^x$、$x = 1$ の中心差分。$h$ を $2^{-26}$ まで小さくすると、誤差は打ち切り誤差の線に沿って下がったあと、丸めの雑音で下げ止まり、観測次数は $2$ から崩れる。">
<img src="/figures/central-diff-order-collapse.gif" alt="刻み h を小さくしていく動画。上の図では中心差分の誤差が打ち切り誤差の直線に沿って減り、h が 10^-5 付近を過ぎると下げ止まって増える。下の図では観測次数が 2 のまま進み、同じ付近で崩れ、さらに細かい刻みでは定義できなくなる。" loading="lazy" />
</Figure>

### 3.4. 製造解の方法

微分方程式を解くコードでは、厳密解の分かる問題が限られます。そこで順序を逆にします。

<Definition id="def-manufactured-solution" title="製造解の方法">
方程式 $\mathcal{L}u = f$ を解くコードを検証するために、解 $u^{\dagger}$ を先に自由に選び、$f^{\dagger} = \mathcal{L}u^{\dagger}$ を手計算（または数式処理）で求めて右辺とし、境界条件も $u^{\dagger}$ から作る。コードに $f^{\dagger}$ を与えて得た数値解と $u^{\dagger}$ の差を、刻みを変えて測る。この $u^{\dagger}$ を**製造解**という。
</Definition>

製造解は物理的に意味のある解である必要はありません。むしろ、方程式のすべての項が効き、どの導関数も $0$ にならない、対称性のない関数を選びます。多項式を選ぶと差分が厳密になって誤差が $0$ になり、次数が観測できないことがあります。

<Example id="ex-manufactured" title="製造解で Poisson 方程式の差分コードを検証する">
$-u'' = f$（$0 \lt x \lt 1$）、$u(0) = u(1) = 0$ を 3 点の中心差分で解くコードを検証します。製造解を $u^{\dagger}(x) = e^x\sin\pi x$ とすると、境界条件を満たし、$u^{\dagger\prime\prime} = e^x\bigl((1 - \pi^2)\sin\pi x + 2\pi\cos\pi x\bigr)$ から $f^{\dagger} = e^x\bigl((\pi^2 - 1)\sin\pi x - 2\pi\cos\pi x\bigr)$ です。期待される次数は $2$ です。比較のため、行列の対角成分を $2$ ではなく $2 + h^2$ とした実装も試します。

```python
# Python 3.12.14 / NumPy 2.5.2
import math
import numpy as np

# 製造解: u(x) = e^x sin(πx) を先に決め、-u'' = f となる f を手で作る
u = lambda x: np.exp(x) * np.sin(np.pi * x)
f = lambda x: np.exp(x) * ((np.pi**2 - 1) * np.sin(np.pi * x) - 2 * np.pi * np.cos(np.pi * x))

def solve_poisson(N, coef=2.0):
    """-u'' = f, u(0) = u(1) = 0 を 3 点差分で解く。coef は対角成分（正しくは 2）"""
    h = 1.0 / N
    x = np.linspace(0.0, 1.0, N + 1)[1:-1]
    A = (np.diag(np.full(N - 1, coef)) - np.diag(np.ones(N - 2), 1) - np.diag(np.ones(N - 2), -1)) / h**2
    return x, np.linalg.solve(A, f(x))

for name, coef in (("正しい差分", 2.0), ("対角成分を 2 + h² と誤った差分", None)):
    errs = []
    for N in (16, 32, 64, 128, 256):
        x, v = solve_poisson(N, 2.0 if coef else 2.0 + (1.0 / N) ** 2)
        errs.append(float(np.max(np.abs(v - u(x)))))
    print(name)
    print("   最大誤差 " + "".join(f"{e:10.2e}" for e in errs))
    print("   観測次数           " + "".join(f"{math.log2(errs[i] / errs[i + 1]):10.3f}" for i in range(4)))
```

```text
正しい差分
   最大誤差   3.85e-03  9.63e-04  2.41e-04  6.02e-05  1.50e-05
   観測次数                2.000     2.000     2.000     2.000
対角成分を 2 + h² と誤った差分
   最大誤差   1.52e-01  1.54e-01  1.54e-01  1.54e-01  1.55e-01
   観測次数               -0.022    -0.004    -0.001    -0.000
```

**検算.** (a) 正しい差分は観測次数 $2$ です。(b) 誤った差分の誤差は刻みを細かくしても減りません（観測次数 $0$）。対角成分に $h^2$ を足して $h^2$ で割ると、方程式に $+u$ の項が加わるので、このコードは別の方程式 $-u'' + u = f$ を 2 次の精度で解いています。数値解は滑らかで、境界条件も満たし、見た目には何の異常もありません。製造解との比較が、方程式そのものの取り違えを検出しています。(c) 製造解の導関数は、$u^{\dagger}$ を数値微分した値との比較でも確かめられます。$f^{\dagger}$ の手計算を誤ると、正しいコードが不合格になるからです。
</Example>

## 4. 残差は誤差ではない

方程式を解いた結果は、解を方程式に代入し直せば確かめられます。代入して残る量が残差です。残差は真の解を知らなくても計算できる便利な量ですが、誤差そのものではありません。

<Proposition id="prop-residual-vs-error" title="残差と誤差">
$A$ を $n$ 次の正則な実行列、$\boldsymbol{b} \ne \boldsymbol{0}$、$A\boldsymbol{x} = \boldsymbol{b}$ とし、$\boldsymbol{y} \in \mathbb{R}^n$ の残差を $\boldsymbol{r} = \boldsymbol{b} - A\boldsymbol{y}$ とする。ベクトルのノルムを 1 つ固定し、行列にはその誘導ノルムを使い、$\kappa(A) = \|A\|\,\|A^{-1}\|$ とする。

1. $\dfrac{\|\boldsymbol{y} - \boldsymbol{x}\|}{\|\boldsymbol{x}\|} \le \kappa(A)\,\dfrac{\|\boldsymbol{r}\|}{\|\boldsymbol{b}\|}$。
2. この不等式で等号が成り立つ例がある。2 ノルムで、$0 \lt \epsilon \le 1$、$A = \begin{pmatrix}1 & 0\\ 0 & \epsilon\end{pmatrix}$、$\boldsymbol{b} = \begin{pmatrix}1\\ 0\end{pmatrix}$、$\boldsymbol{y} = \begin{pmatrix}1\\ 1\end{pmatrix}$ とすると、相対残差は $\epsilon$、相対誤差は $1$、$\kappa_2(A) = 1/\epsilon$ である。
</Proposition>

<Proof of="prop-residual-vs-error">
1. $\boldsymbol{r} = \boldsymbol{b} - A\boldsymbol{y} = A\boldsymbol{x} - A\boldsymbol{y} = A(\boldsymbol{x} - \boldsymbol{y})$ なので $\boldsymbol{x} - \boldsymbol{y} = A^{-1}\boldsymbol{r}$ です。誘導ノルムの定義（<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#def-induced-norm" />）から、任意の行列 $M$ とベクトル $\boldsymbol{v}$ について $\|M\boldsymbol{v}\| \le \|M\|\,\|\boldsymbol{v}\|$ なので、$\|\boldsymbol{y} - \boldsymbol{x}\| \le \|A^{-1}\|\,\|\boldsymbol{r}\|$ です。同じく $\|\boldsymbol{b}\| = \|A\boldsymbol{x}\| \le \|A\|\,\|\boldsymbol{x}\|$ で、$\boldsymbol{b} \ne \boldsymbol{0}$ から $\boldsymbol{x} \ne \boldsymbol{0}$ なので $1/\|\boldsymbol{x}\| \le \|A\|/\|\boldsymbol{b}\|$ です。2 つを掛け合わせると主張の不等式です。

2. $\boldsymbol{x} = (1, 0)^{\mathsf{T}}$ で、$\boldsymbol{r} = \boldsymbol{b} - A\boldsymbol{y} = (1, 0)^{\mathsf{T}} - (1, \epsilon)^{\mathsf{T}} = (0, -\epsilon)^{\mathsf{T}}$ です。$\|\boldsymbol{r}\|_2/\|\boldsymbol{b}\|_2 = \epsilon$、$\|\boldsymbol{y} - \boldsymbol{x}\|_2/\|\boldsymbol{x}\|_2 = \|(0, 1)^{\mathsf{T}}\|_2 = 1$ です。対角行列 $\mathrm{diag}(d_1, d_2)$ の 2 ノルムは $\max_i|d_i|$ です（$\|D\boldsymbol{v}\|_2^2 = d_1^2v_1^2 + d_2^2v_2^2 \le \max_i d_i^2\,\|\boldsymbol{v}\|_2^2$ で、$|d_i|$ が最大の成分に対応する単位ベクトルで等号）。よって $\|A\|_2 = 1$、$\|A^{-1}\|_2 = 1/\epsilon$、$\kappa_2(A) = 1/\epsilon$ で、$\kappa_2(A) \times \epsilon = 1$ は相対誤差に一致します。
</Proof>

$\epsilon = 10^{-10}$ なら、相対残差が $10^{-10}$ でも解は $100\%$ ずれています。残差が教えてくれるのは「$\boldsymbol{y}$ は、少しだけ違う問題の厳密解である」ということ、つまり後退誤差（問題のほうをどれだけ動かせば、この $\boldsymbol{y}$ が厳密解になるか）が小さいことです。これに対し前進誤差は、答えそのもののずれ $\boldsymbol{y} - \boldsymbol{x}$ のことです（残差と後退誤差の関係は Rigal–Gaches の公式 <Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#thm-rigal-gaches" />）。それを前進誤差に直すには条件数が要ります（<Ref to="mathematics/numerical-linear-algebra/norms-and-conditioning#cor-forward-error-bound" />、非線形の問題では <Ref to="mathematics/numerical-analysis/conditioning-and-stability#thm-forward-backward" />）。残差による検算は、条件数の見積もりと組にして初めて誤差の主張になります。

残差がそのまま誤差の上界になる場合もあります。実対称行列の固有値では、単位ベクトル $\boldsymbol{x}$ と Rayleigh 商の残差 $\boldsymbol{r}$ について、$\|\boldsymbol{r}\|$ 以内に必ず固有値があります（<Ref to="mathematics/numerical-linear-algebra/power-and-inverse-iteration#cor-residual-eigen" />）。条件数が現れないのは、対称行列の固有値が摂動に対してつねに良条件だからです。対称でない行列ではこの評価は成り立ちません。

<Example id="ex-eigen-checks" title="固有値計算の検算を重ねる">
乱数で作った $8$ 次の実対称行列の固有値を計算し、5 通りの検算をします。

```python
# Python 3.12.14 / NumPy 2.5.2
import numpy as np

rng = np.random.default_rng(20261008)
n = 8
B = rng.standard_normal((n, n))
A = (B + B.T) / 2                                  # 実対称行列
lam, V = np.linalg.eigh(A)                         # 対称行列用の算法
lam2 = np.sort(np.linalg.eigvals(A).real)          # 別の算法（一般行列用）
scale = np.abs(lam).max()
print(f"トレース      |Σλ - tr A| / max|λ|        = {abs(lam.sum() - np.trace(A)) / scale:.1e}")
print(f"行列式        |Πλ - det A| / |det A|      = {abs(lam.prod() - np.linalg.det(A)) / abs(np.linalg.det(A)):.1e}")
print(f"残差          max ||A v - λ v|| / max|λ|  = {np.linalg.norm(A @ V - V * lam, axis=0).max() / scale:.1e}")
print(f"直交性        ||V^T V - I||_F             = {np.linalg.norm(V.T @ V - np.eye(n)):.1e}")
print(f"別の算法      max |λ - λ'| / max|λ|       = {np.abs(lam - lam2).max() / scale:.1e}")
```

```text
トレース      |Σλ - tr A| / max|λ|        = 3.8e-16
行列式        |Πλ - det A| / |det A|      = 1.1e-15
残差          max ||A v - λ v|| / max|λ|  = 6.6e-16
直交性        ||V^T V - I||_F             = 2.2e-15
別の算法      max |λ - λ'| / max|λ|       = 1.9e-15
```

5 つの検算は、単位丸め $u\approx1.1\times10^{-16}$ の数倍から 20 倍程度の範囲で合格しています。それぞれが見ているものは違います。トレースは固有値の和、行列式は積だけを見るので、2 つの固有値が逆向きにずれる誤りを見逃します。残差は固有値と固有ベクトルの組を 1 つずつ見ますが、<Ref to="mathematics/numerical-linear-algebra/power-and-inverse-iteration#cor-residual-eigen" /> が保証するのは「近くに固有値がある」ことで、8 個の固有値を漏れなく求めたことではありません。同じ固有値を 2 回数えていないことは直交性が補います。別の算法との一致は最も強い検査に見えますが、2 つの算法は同じライブラリの同じ入力 `A` を使っているので、`A` を作る段階の誤りには無力です（§6）。
</Example>

## 5. 保存量・対称性・恒等式

### 5.1. 保存量のずれ

微分方程式の厳密解が保つ量（エネルギー、運動量、確率の総和）は、数値解がどれだけ保っているかを調べることで検算に使えます。数値解は一般に保存量を厳密には保ちませんが、ずれの大きさには方法の次数で決まる上界があります。

<Proposition id="prop-conservation-drift" title="保存量のずれの次数">
$U \subset \mathbb{R}^d$ とし、各 $h \in (0, h_0]$ について写像 $\varphi_h, \Phi_h : U \to \mathbb{R}^d$ が与えられているとする（$\varphi_h$ は微分方程式の厳密解を時間 $h$ だけ進める写像、$\Phi_h$ は数値解法の 1 ステップ）。$E$ は $\mathbb{R}^d$ の部分集合 $V$ の上の実数値関数で、次を仮定する。

- （保存）$y \in U$ ならば $\varphi_h(y) \in V$ で $E(\varphi_h(y)) = E(y)$。
- （局所誤差が $p + 1$ 次）定数 $c \ge 0$、$p \gt 0$ があって、$y \in U$、$0 \lt h \le h_0$ について $\Phi_h(y) \in V$ で $\|\Phi_h(y) - \varphi_h(y)\| \le ch^{p+1}$。
- （Lipschitz 連続）定数 $L_E \ge 0$ があって、$z, w \in V$ について $|E(z) - E(w)| \le L_E\|z - w\|$。

$y_{k+1} = \Phi_h(y_k)$ で定まる数値解が $y_0, \ldots, y_{n-1} \in U$ を満たし、$T = nh$ とすると
$$
|E(y_{k+1}) - E(y_k)| \le L_Ec\,h^{p+1}\quad(0 \le k \lt n), \qquad |E(y_n) - E(y_0)| \le L_Ec\,T\,h^{p}
$$
である。
</Proposition>

<Proof of="prop-conservation-drift">
$0 \le k \lt n$ とします。保存の仮定から $E(y_k) = E(\varphi_h(y_k))$ なので、
$$
|E(y_{k+1}) - E(y_k)| = |E(\Phi_h(y_k)) - E(\varphi_h(y_k))| \le L_E\|\Phi_h(y_k) - \varphi_h(y_k)\| \le L_Ec\,h^{p+1}
$$
です（$\Phi_h(y_k), \varphi_h(y_k) \in V$ なので Lipschitz 連続性が使えます）。$E(y_n) - E(y_0) = \sum_{k=0}^{n-1}\bigl(E(y_{k+1}) - E(y_k)\bigr)$ に三角不等式を使うと $|E(y_n) - E(y_0)| \le nL_Ec\,h^{p+1} = L_Ec\,(nh)\,h^p$ です。
</Proof>

証明は、数値解の誤差そのもの（$y_n$ と厳密解の差）の評価を経由していません。1 ステップのずれを足しただけです。そのため、この上界は長い時間 $T$ について線形にしか増えません。検算としての使い方は、ずれの大きさと次数を測り、**上界より大きいか、次数が $p$ より低ければ**実装の誤りを疑う、というものです。上界なので、次数が $p$ より高く出ることはあります。

<Example id="ex-energy-drift" title="調和振動子のエネルギーのずれ">
$q' = p$、$p' = -q$（保存量 $E = (q^2 + p^2)/2$）を、初期値 $(1, 0)$ から $T = 8$ まで、$n$ ステップで解きます。Euler 法（局所誤差 $O(h^2)$）、4 段 4 次の Runge–Kutta 法（局所誤差 $O(h^5)$）、厳密解の 1 ステップ（角 $h$ の回転）の 3 通りです。

```python
# Python 3.12.14 / NumPy 2.5.2
import math
import numpy as np

# 調和振動子 q' = p, p' = -q。保存量 E = (q^2 + p^2)/2、初期値 (1, 0)、T = 8
J = np.array([[0.0, 1.0], [-1.0, 0.0]])
f = lambda y: J @ y
E = lambda y: 0.5 * float(y @ y)

def euler(y, h):
    return y + h * f(y)

def rk4(y, h):
    k1 = f(y); k2 = f(y + h / 2 * k1); k3 = f(y + h / 2 * k2); k4 = f(y + h * k3)
    return y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

def rotation(y, h):                    # 厳密解の 1 ステップ（回転）
    c, s = math.cos(h), math.sin(h)
    return np.array([c * y[0] + s * y[1], -s * y[0] + c * y[1]])

def drift(step, n, T=8.0):
    y = np.array([1.0, 0.0])
    for _ in range(n):
        y = step(y, T / n)
    return abs(E(y) - 0.5)

ns = [64, 128, 256, 512, 1024]
for name, step in (("Euler 法", euler), ("RK4", rk4), ("厳密な回転", rotation)):
    d = [drift(step, n) for n in ns]
    orders = ["     -"] + [f"{math.log2(d[i] / d[i + 1]):6.2f}" if d[i] > 0 and d[i + 1] > 0 else "     -" for i in range(len(ns) - 1)]
    print(name)
    print("   n       " + "".join(f"{n:>10d}" for n in ns))
    print("   |ΔE|    " + "".join(f"{v:10.2e}" for v in d))
    print("   観測次数" + "".join(f"{o:>10s}" for o in orders))
print(f"n·u (n = 1024) = {1024 * 2.0**-53:.2e}")
```

```text
Euler 法
   n               64       128       256       512      1024
   |ΔE|      8.49e-01  3.24e-01  1.42e-01  6.66e-02  3.22e-02
   観測次数         -      1.39      1.19      1.09      1.05
RK4
   n               64       128       256       512      1024
   |ΔE|      1.69e-06  5.30e-08  1.66e-09  5.17e-11  1.62e-12
   観測次数         -      5.00      5.00      5.00      5.00
厳密な回転
   n               64       128       256       512      1024
   |ΔE|      2.66e-15  5.00e-15  9.55e-15  1.63e-14  5.60e-14
   観測次数         -     -0.91     -0.93     -0.77     -1.78
n·u (n = 1024) = 1.14e-13
```

**検算.** (a) Euler 法のずれの観測次数は $1$ に近づいています。<Ref to="prop-conservation-drift" /> の $p = 1$ と合います。この問題では 1 ステップで $q^2 + p^2$ がちょうど $1 + h^2$ 倍になるので、ずれは $\frac12\bigl((1 + h^2)^n - 1\bigr)$ と閉じた形で書け、$n = 1024$ で $0.0322$ です（演習）。(b) RK4 の観測次数は $5$ で、命題の $p = 4$ より高く出ています。命題は上界なので矛盾しません。この線形の問題では、RK4 の 1 ステップは複素数 $R(ih)$（$R(z) = 1 + z + z^2/2 + z^3/6 + z^4/24$）を掛けることに当たり、$|R(ih)|^2 = (1 - h^2/2 + h^4/24)^2 + (h - h^3/6)^2 = 1 - h^6/72 + h^8/576$ なので、1 ステップのずれが $O(h^5)$ ではなく $O(h^6)$ になるためです。$n = 1024$ での値は $\frac12 \cdot n \cdot h^6/72 = \frac{T h^5}{144} \approx 1.6 \times 10^{-12}$ で、表と一致します。(c) 厳密な回転のずれは刻みを細かくすると**増え**、大きさは $n \cdot u$ の程度以下です。これは打ち切り誤差ではなく、1 ステップごとの丸め誤差が積もったものです。保存量を厳密に保つ方法のずれは、このように $h$ のべきではなく、ステップ数に比例する丸めの大きさで現れます。2 つを取り違えないようにします。
</Example>

### 5.2. 対称性と恒等式による検査の限界

行列式の乗法性、Parseval の等式、確率の和が 1 であること、回転しても変わらないはずの量、といった恒等式や対称性は、答えを知らなくても検査できる便利な性質です。その限界は、次の形で特徴づけられます。

<Proposition id="prop-invariant-test" title="不変性の検査が検出できる誤り">
$X$ を集合、$\mathcal{G}$ を $X$ から $X$ への写像の集合、$W$ をベクトル空間とし、$f : X \to W$ がすべての $g \in \mathcal{G}$、$x \in X$ について $f(g(x)) = f(x)$ を満たすとする（$f$ は $\mathcal{G}$ について不変）。$f$ の実装を $\tilde f = f + \delta$（$\delta : X \to W$ は実装の誤り）とする。

1. 検査「$\tilde f(g(x)) = \tilde f(x)$」が入力 $(g, x)$ で合格することは、$\delta(g(x)) = \delta(x)$ と同値である。したがって、$\mathcal{G}$ について不変な誤り $\delta$ は、検査の入力をどう選んでも検出されない。
2. $\sigma : X \to X$ を写像とし、実装が $\tilde f = f \circ \sigma$（入力を $\sigma$ で変換してから正しい計算をする、という誤り）であるとする。このとき、$\sigma(x) = x$ を満たす入力 $x$ の上では $\tilde f(x) = f(x)$ で、正しい値とのどんな比較も合格する。
</Proposition>

<Proof of="prop-invariant-test">
1. $\tilde f(g(x)) - \tilde f(x) = \bigl(f(g(x)) - f(x)\bigr) + \bigl(\delta(g(x)) - \delta(x)\bigr) = \delta(g(x)) - \delta(x)$ です（$f$ の不変性）。左辺が $0$ であることと $\delta(g(x)) = \delta(x)$ は同値です。$\delta$ が不変なら、これはすべての $(g, x)$ で成り立ちます。

2. $\sigma(x) = x$ なら $\tilde f(x) = f(\sigma(x)) = f(x)$ です。
</Proof>

1 は、不変性の検査が見ているのは「誤りのうち、対称性を破る成分」だけだと言っています。全体に定数倍が掛かる誤り、規格化の誤りは、たいていの対称性を保つので素通りします。恒等式でも同じです。行列式の乗法性 $\det(AB) = \det A\det B$ の検査は、$\det$ の代わりに $(\det)^2$ や $|\det|$ や定数 $1$ を返す実装も合格させます。どれも乗法的だからです。

2 は、**検査に使う入力が対称すぎる**と誤りが隠れることを言っています。転置の取り違えは対称行列では見えず、添字の順序の誤りは正方行列や立方体の格子では見えず、符号の誤りは偶関数では見えません。検査の入力は、対称性のない一般的なものを選びます。

<Example id="ex-transpose-bug" title="転置を取り違えた行列ベクトル積はどの検査で見つかるか">
行列ベクトル積 $A\boldsymbol{x}$ を、添字を取り違えて $A^{\mathsf{T}}\boldsymbol{x}$ と計算している実装を、5 通りの検査にかけます。

```python
# Python 3.12.14 / NumPy 2.5.2
import numpy as np

def matvec_ref(A, x):                 # 定義どおりの二重ループ（別手法）
    n, m = A.shape
    return np.array([sum(A[i, j] * x[j] for j in range(m)) for i in range(n)])

def matvec_bug(A, x):                 # 誤り: 添字を取り違えて A^T x を計算している
    return np.einsum("ji,j->i", A, x)

def rel(a, b):                        # 相対的な食い違い
    return float(np.linalg.norm(a - b) / max(np.linalg.norm(a), np.linalg.norm(b)))

rng = np.random.default_rng(20261008)
n = 6
A = rng.standard_normal((n, n))
S = (A + A.T) / 2                     # 対称行列
x, y = rng.standard_normal(n), rng.standard_normal(n)
f = matvec_bug

tests = {
    "線形性 f(A,2x-3y) = 2f(A,x) - 3f(A,y)": rel(f(A, 2 * x - 3 * y), 2 * f(A, x) - 3 * f(A, y)),
    "随伴の恒等式 <f(A,x),y> = <x,f(A^T,y)>": abs(f(A, x) @ y - x @ f(A.T, y)) / abs(f(A, x) @ y),
    "別手法との比較（対称行列 S）": rel(f(S, x), matvec_ref(S, x)),
    "別手法との比較（一般の行列 A）": rel(f(A, x), matvec_ref(A, x)),
    "既知の値 (E_12 e_2 = e_1)": rel(f(np.outer(np.eye(n)[0], np.eye(n)[1]), np.eye(n)[1]), np.eye(n)[0]),
}
for name, d in tests.items():
    print(f"{name:42s} 食い違い = {d:.2e}  {'検出' if d > 1e-10 else '見逃し'}")
```

```text
線形性 f(A,2x-3y) = 2f(A,x) - 3f(A,y)         食い違い = 2.84e-16  見逃し
随伴の恒等式 <f(A,x),y> = <x,f(A^T,y)>           食い違い = 2.43e-16  見逃し
別手法との比較（対称行列 S）                            食い違い = 3.46e-17  見逃し
別手法との比較（一般の行列 A）                           食い違い = 1.27e+00  検出
既知の値 (E_12 e_2 = e_1)                      食い違い = 1.00e+00  検出
```

線形性と随伴の恒等式は、誤った実装でも**すべての入力で**厳密に成り立ちます。$A^{\mathsf{T}}\boldsymbol{x}$ も $\boldsymbol{x}$ について線形で、$\langle A^{\mathsf{T}}\boldsymbol{x}, \boldsymbol{y}\rangle = \langle\boldsymbol{x}, A\boldsymbol{y}\rangle = \langle\boldsymbol{x}, (A^{\mathsf{T}})^{\mathsf{T}}\boldsymbol{y}\rangle$ だからです。転置の取り違えは、この 2 つの恒等式を保つ誤りです（<Ref to="prop-invariant-test" /> の 1 に当たる状況）。対称行列での比較が合格するのは同じ命題の 2（$\sigma$ が転置、$S^{\mathsf{T}} = S$）です。検出できたのは、対称でない入力での別手法との比較と、手で答えの分かる最小の入力だけでした。乱数で作った行列は確率 1 で対称でないので、**検査の入力を乱数で作る**ことには、入力の対称性による見逃しを避ける効果があります。
</Example>

## 6. 検算の独立性と、反例探索で言えること

### 6.1. 2 つの検算は独立か

検算を 2 つ重ねれば、両方が見逃す確率は小さくなるはずです。どれだけ小さくなるかは、2 つの検算の「見逃し方」が独立かどうかで決まります。

<Proposition id="prop-independent-checks" title="共通の原因があると見逃しの確率は積にならない">
ある誤りを検算 $i$（$i = 1, 2$）が見逃すという事象を $M_i$ とし、$m_i = P(M_i)$ とおく。

1. $M_1$ と $M_2$ が独立ならば $P(M_1 \cap M_2) = m_1m_2$。
2. 事象 $Z$（共通の原因。確率 $c$、$0 \lt c \lt 1$）が起こると両方の検算が必ず見逃し、$Z$ が起こらないときは、2 つの検算が条件付き確率 $P(\,\cdot \mid Z^{c})$ のもとで独立に確率 $a_1, a_2$ で見逃すとする。このとき
   $$
   P(M_1 \cap M_2) = c + (1 - c)a_1a_2 \ge c, \qquad P(M_1 \cap M_2) - m_1m_2 = c(1 - c)(1 - a_1)(1 - a_2) \ge 0
   $$
   である。
</Proposition>

<Proof of="prop-independent-checks">
1. 事象の独立性の定義（<Ref to="mathematics/probability/probability-spaces#def-independence" />）そのものです。

2. 全確率の公式（<Ref to="mathematics/probability/probability-spaces#thm-bayes" /> の 2）を分割 $\{Z, Z^{c}\}$ に使います。仮定から $P(M_1 \cap M_2 \mid Z) = 1$、$P(M_1 \cap M_2 \mid Z^{c}) = a_1a_2$ なので $P(M_1 \cap M_2) = c + (1 - c)a_1a_2$ で、第 2 項は非負なので $c$ 以上です。同じく $m_i = P(M_i) = c + (1 - c)a_i$ です。積を展開すると
   $$
   m_1m_2 = c^2 + c(1 - c)(a_1 + a_2) + (1 - c)^2a_1a_2
   $$
   で、
   $$
   P(M_1 \cap M_2) - m_1m_2 = c(1 - c) - c(1 - c)(a_1 + a_2) + (1 - c)a_1a_2\bigl(1 - (1 - c)\bigr) = c(1 - c)\bigl(1 - a_1 - a_2 + a_1a_2\bigr)
   $$
   です。括弧の中は $(1 - a_1)(1 - a_2) \ge 0$ です。
</Proof>

数値を入れてみます。それぞれ単独では $1\%$ しか見逃さない検算（$a_1 = a_2 = 0.01$）が 2 つあるとき、独立なら両方が見逃す確率は $10^{-4}$ です。しかし共通の原因が $5\%$ の確率で起こる（$c = 0.05$）なら、$0.05 + 0.95 \times 10^{-4} \approx 0.050$ で、500 倍です。2 つ目の検算は、共通の原因に対しては何も足していません。

数値計算での共通の原因は、ありふれています。同じ（誤って導いた）公式を 2 通りに実装する。2 つの算法が同じライブラリの同じ関数を呼ぶ。同じ入力ファイル、同じ単位の取り違え、同じ仕様の誤解を共有する。同じ人が、同じ思い込みのまま 2 回書く。「別手法による検算」が意味をもつのは、**公式の導出まで遡って別**であるときです。差分法の結果をスペクトル法で、数値積分を級数展開で、数値解を保存則で確かめる、というように、誤りの入り方が違う方法を選びます。

### 6.2. 反例が見つからなかったとき

命題を多数の入力で試して反例が出なかった、という報告は数値実験の定番です。そこから言えることを正確に述べます。

<Theorem id="thm-random-search-bound" title="無作為な反例探索で言える範囲">
入力の集合の上の確率分布を 1 つ固定し、命題の反例である入力の集合の確率を $\theta \in [0, 1]$ とする。この分布に従う独立な $N \ge 1$ 個の標本 $X_1, \ldots, X_N$ を調べる。

1. 反例が 1 つも出ない確率は $(1 - \theta)^N$ で、$(1 - \theta)^N \le e^{-\theta N}$。
2. $0 \lt \delta \lt 1$ とする。$\theta \gt \ln(1/\delta)/N$ ならば、反例が 1 つも出ない確率は $\delta$ より小さい。
3. どの $N$ についても $(1 - \theta)^N \ge 1 - N\theta$。とくに $\theta \le 1/(100N)$ である偽の命題（$\theta \gt 0$）は、確率 $0.99$ 以上で反例が出ないまま検査を通過する。
</Theorem>

<Proof of="thm-random-search-bound">
1. 反例の集合を $B$ とすると、反例が出ないことは $X_1 \notin B, \ldots, X_N \notin B$ で、独立性（<Ref to="mathematics/probability/random-variables#def-independence" />）から確率は $\prod_{i=1}^N P(X_i \notin B) = (1 - \theta)^N$ です。$\phi(t) = e^{-t} - 1 + t$ は $\phi(0) = 0$、$\phi'(t) = 1 - e^{-t} \ge 0$（$t \ge 0$）なので、$[0, \infty)$ で広義単調増加で（<Ref to="mathematics/calculus/mean-value-and-taylor#cor-monotonicity" /> の 2）$\phi(\theta) \ge 0$、すなわち $0 \le 1 - \theta \le e^{-\theta}$ です。両辺を $N$ 乗して主張を得ます。

2. 1 から、確率は $e^{-\theta N}$ 以下で、$\theta N \gt \ln(1/\delta)$ なら $e^{-\theta N} \lt e^{-\ln(1/\delta)} = \delta$ です。

3. $N$ についての帰納法です。$N = 1$ は等号です。$N$ で成り立つとすると、両辺に $1 - \theta \ge 0$ を掛けて $(1 - \theta)^{N+1} \ge (1 - N\theta)(1 - \theta) = 1 - (N + 1)\theta + N\theta^2 \ge 1 - (N + 1)\theta$ です。$\theta \le 1/(100N)$ なら $1 - N\theta \ge 0.99$ です。
</Proof>

2 は次のように読みます。$N$ 個の標本で反例が出なかったとき、「$\theta \le \ln(1/\delta)/N$」という主張は、信頼水準 $1 - \delta$ で成り立ちます（$\theta$ がそれより大きければ、反例が出ないという観測は確率 $\delta$ 未満でしか起こらないからです）。$\delta = 0.05$ なら $\ln 20 \approx 3.0$ で、「反例の割合は $3/N$ 以下」です（**3/N 則**）。$10^6$ 個で反例なしなら、信頼水準 $95\%$ で $\theta \le 3 \times 10^{-6}$ です。

言えるのはここまでです。3 のとおり、反例の割合が $1/N$ よりずっと小さい偽の命題は、ほぼ確実に検査を通過します。そして $\theta = 0$ であっても、命題が真であるとは限りません。

<Example id="ex-missed-counterexamples" title="無作為な探索が見逃す反例">
**測度 0 の反例.** 「すべての $n$ 次実正方行列は正則である」は偽ですが、成分を連続分布から独立に引いた行列が特異になる確率は $0$ です。特異行列の集合は多項式 $\det A = 0$ の零点集合で、体積が $0$ だからです。$\theta = 0$ で、何個試しても反例は出ません。重根をもつ多項式、退化した固有値、分母が $0$ になるパラメータなど、数値計算で問題を起こす入力の多くはこの型です。

**標本の分布に入らない反例.** 「binary64 の数 $x$ について $\mathrm{fl}(x + 1) \gt x$」は、$[0, 1)$ の一様乱数で試す限り反例が出ません。反例は $x \ge 2^{53}$ にあります（$2^{53} + 1$ は表現できず、$2^{53}$ に丸められます）。探索の範囲の外にある反例については、探索は何も言いません。

**はるか遠くにある反例.** 素数の個数 $\pi(x)$ と対数積分 $\mathrm{li}(x)$ について、$\pi(x) \lt \mathrm{li}(x)$ は計算された範囲のすべての $x$ で成り立ちますが、Littlewood は 1914 年に、この不等号が無限回入れ替わることを証明しました。最初に入れ替わる $x$ は直接の計算が届かないほど大きいところにあります。計算できる範囲での一致が、命題の証拠としてどれほど弱いことがあるかを示す例です。

したがって、探索の報告には、標本の分布・範囲・個数・乱数種を必ず添えます。それらが「何を調べていないか」を読者に伝える情報だからです。
</Example>

## 7. 数値実験を命題に昇格させる

[数値実験](/computer-science/research-computing/reproducible-experiments#def-experiment)で見つけた事実を、論文や報告の中で「成り立つ」と書くまでには、いくつかの段があります。

**(i) 主張を量化子つきで書く。** 「誤差は小さい」ではなく、「すべての $N \in \{8, 16, \ldots, 256\}$ について誤差は $10^{-5}$ 以下」「ある定数 $C$ があって、すべての $h \le h_0$ で誤差は $Ch^4$ 以下」のように、どの範囲のすべてについての主張か、何の存在の主張かを書きます。有限個の場合についての主張は計算で確かめられます。無限個の場合についての主張は、計算だけでは確かめられません。

**(ii) 数値の主張に誤差の上界を付ける。** 「値は $0.7468$」ではなく「値は $0.746824 \pm 10^{-6}$」と書き、$\pm$ の根拠（<Ref to="thm-three-level-order" /> の誤差の推定、条件数と残差、[誤差予算](/computer-science/research-computing/precision-and-error-budget)）を示します。推定値と保証された上界は区別します。Richardson の推定は漸近的に正しい推定値で、上界ではありません。

**(iii) 探索の範囲・件数・種・失敗例を記録する。** <Ref to="thm-random-search-bound" /> のとおり、反例が出なかったという報告の価値は、分布と件数で決まります。途中で止まった試行、収束しなかった試行も消さずに件数を残します（[長時間ジョブの管理](/computer-science/research-computing/long-running-jobs)の章の打ち切りの偏り）。

**(iv) 証明に使う数値は、包含の保証がある計算に置き換える。** 浮動小数の計算結果は、真の値の近くにあると期待できるだけで、どちら側にどれだけずれているかを教えてくれません。証明の中で「この量は正である」「この区間に零点がある」という事実を使うなら、結果が真の値を含む区間であることが保証された計算を使います。区間演算では、算術式の値が計算された区間に必ず含まれます（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-inclusion" />）。球算術でも同じです（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-ball-inclusion" />）。区間が $0$ を含まなければ符号が確定し、区間 Newton 法の条件が満たされれば零点の存在と一意性が確定します（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-interval-newton" />）。

**(v) 「数値的に確かめた」と「証明した」を書き分ける。** 「$N \le 10^6$ のすべてについて計算で確かめた」「$10^6$ 個の乱数（種 20261008、$[0,1)$ の一様分布）で反例は見つからなかった」「球算術により証明した」は、別の強さの主張です（<Ref to="def-evidence-levels" />）。読者が区別できるように書きます。

とくに危ういのは、**ほぼ $0$ の量の符号や、等しいかどうかを浮動小数で判定すること**です。

<Example id="ex-sign-of-near-zero" title="倍精度で正に見える量が、実は 0 である例">
2 つの量を考えます。$q_1 = \sqrt2\sqrt3 - \sqrt6$ は厳密に $0$ です（$\sqrt2\sqrt3$ も $\sqrt6$ も正で、2 乗するとどちらも $6$）。$q_2 = e^{\pi\sqrt{163}} - 262537412640768744$ は、$0$ にきわめて近い、$0$ でない数です。倍精度、mpmath の 50 桁、python-flint の球（`arb`）で計算します。

```python
# Python 3.12.14 / mpmath 1.4.1 / python-flint 0.9.0
import math
from mpmath import mp, mpf, sqrt as msqrt, exp as mexp, pi as mpi, nstr
from flint import arb, ctx

N163 = 262537412640768744

# 倍精度
q1 = math.sqrt(2) * math.sqrt(3) - math.sqrt(6)
q2 = math.exp(math.pi * math.sqrt(163)) - N163
print(f"float64 : q1 = {q1:+.3e}   q2 = {q2:+.3e}")

# mpmath（10 進 50 桁の浮動小数。包含の保証はない）
mp.dps = 50
m1 = msqrt(2) * msqrt(3) - msqrt(6)
m2 = mexp(mpi * msqrt(163)) - N163
print(f"mpmath50: q1 = {nstr(m1, 4)}   q2 = {nstr(m2, 8)}")

# python-flint の球（真の値を含む区間）
for prec in (53, 128, 256, 1024):
    ctx.prec = prec
    b1 = arb(2).sqrt() * arb(3).sqrt() - arb(6).sqrt()
    b2 = (arb.pi() * arb(163).sqrt()).exp() - arb(N163)
    print(f"arb {prec:4d} ビット: q1 = {b1.str(5)}  0 を含む: {b1.contains(0)}")
    print(f"               q2 = {b2.str(8)}  0 を含む: {b2.contains(0)}  負と確定: {bool(b2 < 0)}")
```

```text
float64 : q1 = +4.441e-16   q2 = -4.800e+02
mpmath50: q1 = -5.346e-51   q2 = -7.499274e-13
arb   53 ビット: q1 = [+/- 2.04e-15]  0 を含む: True
               q2 = [+/- 6.49e+3]  0 を含む: True  負と確定: False
arb  128 ビット: q1 = [+/- 4.21e-38]  0 を含む: True
               q2 = [-7.49927e-13 +/- 5.70e-19]  0 を含む: False  負と確定: True
arb  256 ビット: q1 = [+/- 1.24e-76]  0 を含む: True
               q2 = [-7.4992740e-13 +/- 2.81e-21]  0 を含む: False  負と確定: True
arb 1024 ビット: q1 = [+/- 1.02e-307]  0 を含む: True
               q2 = [-7.4992740e-13 +/- 2.81e-21]  0 を含む: False  負と確定: True
```

**$q_1$ について.** 倍精度は $+4.4 \times 10^{-16}$（正）、mpmath の 50 桁は $-5.3 \times 10^{-51}$（負）を返します。真の値は $0$ で、どちらの符号も丸めの産物です。「2 つの精度で計算して、どちらも $0$ でなかった」ことは、$0$ でないことの証拠になりません。球算術は、どの精度でも「$0$ を含む」と答えます。これは「分からない」という正直な答えです。精度を上げると半径は縮みますが、$0$ を含むことは変わらず、**有限の精度の計算で等号を証明することはできません**。$q_1 = 0$ は、冒頭に書いた代数的な議論で証明します。

**$q_2$ について.** 倍精度の $-480$ は、桁が合っていません。$e^{\pi\sqrt{163}} \approx 2.6 \times 10^{17}$ の付近では隣り合う binary64 の数の間隔が $32$ で、さらに指数関数の引数 $\pi\sqrt{163} \approx 40$ の丸め誤差が約 $40$ 倍に拡大されるからです。mpmath の 50 桁は $-7.499274 \times 10^{-13}$ を返します。これはもっともらしい値ですが、包含の保証はありません。球算術は、53 ビットでは「分からない」、128 ビット以上では「$0$ を含まない、負である」と答えます。これは、ライブラリの球の演算が真の値を含むという保証のもとでの、$q_2 \lt 0$ の証明です。$0$ でない量の符号は、十分な精度の球算術で確定できます。

**検算.** $q_2$ は、3 通りの計算（mpmath の 50 桁、`arb` の 128 ビットと 256 ビット）で $-7.4992740 \times 10^{-13}$ の表示の桁まで一致します。倍精度の値だけが外れており、その外れ方（$32$ の倍数で、$10^3$ の程度）は上の説明と合います。
</Example>

この例の球算術の保証について、正確に述べておきます。Rikai で証明してあるのは、球の加算と乗算が真の値を含むこと（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-ball-inclusion" />）です。平方根・指数関数・円周率の球が真の値を含むことは、ライブラリ（Arb）の実装が保証する性質で、ここではそれを信頼しています。計算機を使った証明では、「何を信頼したか」（ライブラリ、その版、コンパイラ、ハードウェアの演算）を明記するのが作法です。この先の、数値を証明に変える技法は、[浮動小数点数と区間演算](/mathematics/computer-assisted-proofs/floating-point-interval-arithmetic)の章から始まる精度保証つき計算のテーマで扱います。

なお、この章で扱ったのは数値の検算です。整列の結果が整列されているか、データ構造の不変条件が保たれているか、といった離散的な性質を無作為な入力で調べる方法（入力の生成と、失敗した入力の縮約）は、ソフトウェアのテストの理論に属する別の話題です。

## 8. 演習

<Exercise id="exr-observed-order-read" difficulty="易">
ある求積のコードの誤差が、刻み $h, h/2, h/4$ でそれぞれ $3.2 \times 10^{-3}$、$8.1 \times 10^{-4}$、$2.0 \times 10^{-4}$ でした。観測次数を求めてください。このコードが複合 Simpson 則のつもりで書かれていたとすると、何を疑いますか。
<Solution>
$\log_2(3.2/0.81) \approx 1.98$、$\log_2(8.1/2.0) \approx 2.02$ で、観測次数は $2$ です。

複合 Simpson 則なら、被積分関数が $C^4$ 級のとき次数は $4$ のはずです。疑う先は 2 つあります。1 つは実装で、たとえば重み $4$ と $2$ の位置が入れ替わっている、台形則の重みになっている、などです。もう 1 つは仮定で、被積分関数の高階の導関数が積分区間で有界でない場合です（<Ref to="ex-simpson-orders" /> の $\sqrt x$ は次数 $1.5$ でした）。滑らかな関数（$e^x$ など）で同じコードを試し、それでも $2$ なら実装の誤りです。
</Solution>
</Exercise>

<Exercise id="exr-three-level-hand" difficulty="標準">
真値の分からない量について、刻み $h, h/2, h/4$ での計算値が $A(h) = 1.2000$、$A(h/2) = 1.0500$、$A(h/4) = 1.0125$ でした。<Ref to="thm-three-level-order" /> の方法で次数を推定し、$A(h/4)$ の誤差の推定値と、それを使って補正した値を求めてください。
<Solution>
$A(h) - A(h/2) = 0.1500$、$A(h/2) - A(h/4) = 0.0375$ で、比は $4$、推定次数は $\log_2 4 = 2$ です。

誤差の推定は $A - A(h/4) \approx \bigl(A(h/4) - A(h/2)\bigr)/(2^2 - 1) = -0.0375/3 = -0.0125$ です。補正した値は $1.0125 - 0.0125 = 1.0000$ です。これは Richardson 外挿の値 $(4A(h/4) - A(h/2))/3$ と同じものです。

ただし、3 つの値から次数を 1 つ推定しただけなので、この次数が漸近領域のものかどうかは確かめられていません。もう 1 段細かい $A(h/8)$ を計算し、$A(h/2), A(h/4), A(h/8)$ からの推定次数も $2$ に近いことを確かめるのが安全です。
</Solution>
</Exercise>

<Exercise id="exr-rule-of-three" difficulty="標準">
ある不等式を、種を固定した $10^6$ 個の一様乱数の入力で調べ、反例は見つかりませんでした。(i) 信頼水準 $95\%$ で、反例の割合について何が言えますか。(ii) 信頼水準 $99\%$ で「反例の割合は $10^{-9}$ 以下」と言うには、何個の標本が必要ですか。(iii) (ii) の個数を調べて反例がなかったとして、不等式が証明されたと言えない理由を 2 つ挙げてください。
<Solution>
(i) <Ref to="thm-random-search-bound" /> の 2 で $\delta = 0.05$ とすると、$\theta \le \ln 20/10^6 \approx 3.0 \times 10^{-6}$ です。

(ii) $\ln(1/0.01)/N \le 10^{-9}$ から $N \ge \ln 100 \times 10^9 \approx 4.6 \times 10^9$ 個です。

(iii) 第一に、$\theta \le 10^{-9}$ は $\theta = 0$ ではなく、$\theta = 0$ であっても反例の集合が空とは限りません（測度 $0$ の反例、標本の分布の外の反例。<Ref to="ex-missed-counterexamples" />）。第二に、各標本での判定は浮動小数の計算によるので、不等式の両辺の差が丸め誤差の程度の入力では、判定そのものが信頼できません（<Ref to="ex-sign-of-near-zero" />）。等号に近い場合がある不等式ほど、この問題は深刻です。証明にするには、入力の集合を有限個の区間に分け、各区間の上で区間演算により不等式を確かめる、といった方法をとります。
</Solution>
</Exercise>

<Exercise id="exr-euler-energy" difficulty="標準">
<Ref to="ex-energy-drift" /> の調和振動子に Euler 法を使うと、1 ステップで $q^2 + p^2$ がちょうど $1 + h^2$ 倍になることを示し、$E(y_0) = 1/2$、$T = nh$ のときのエネルギーのずれが $\frac12\bigl((1 + h^2)^n - 1\bigr) \le \frac12(e^{Th} - 1)$ であることを導いてください。$T = 8$、$n = 1024$ のときの値を求め、例の表と比べてください。
<Solution>
Euler 法の 1 ステップは $(q, p) \mapsto (q + hp,\ p - hq)$ です。$(q + hp)^2 + (p - hq)^2 = q^2 + 2hqp + h^2p^2 + p^2 - 2hpq + h^2q^2 = (1 + h^2)(q^2 + p^2)$ です。したがって $E(y_n) = (1 + h^2)^nE(y_0)$ で、ずれは $\frac12\bigl((1 + h^2)^n - 1\bigr)$ です。$\psi(t) = e^t - 1 - t$ は $\psi(0) = 0$、$\psi'(t) = e^t - 1 \ge 0$（$t \ge 0$）を満たすので $t \ge 0$ で $1 + t \le e^t$ で、$t = h^2$ として $(1 + h^2)^n \le e^{nh^2} = e^{Th}$ です。

$T = 8$、$n = 1024$ では $h = 1/128$、$nh^2 = Th = 0.0625$ で、$(1 + h^2)^n = \exp\bigl(n\ln(1 + h^2)\bigr) \approx \exp(0.062498) \approx 1.06449$、ずれは $\approx 0.0322$ です。表の $3.22 \times 10^{-2}$ と一致します。$Th$ が小さいとき $\frac12(e^{Th} - 1) \approx Th/2$ で、ずれは $h$ の 1 次です（<Ref to="prop-conservation-drift" /> の $p = 1$）。Euler 法のエネルギーは増える一方で、時間について指数関数的に増えます。命題の上界が $T$ について線形なのは、数値解が有界な集合 $U$ にとどまる範囲での評価だからです。
</Solution>
</Exercise>

## 参考文献

- P. J. Roache, *Verification and Validation in Computational Science and Engineering*, Hermosa Publishers, 1998 — 検証と妥当性確認の区別、格子収束による次数の検査。
- W. L. Oberkampf, C. J. Roy, *Verification and Validation in Scientific Computing*, Cambridge University Press, 2010 — 製造解の方法、観測次数。
- K. Salari, P. Knupp, "Code Verification by the Method of Manufactured Solutions", Sandia National Laboratories, SAND2000-1444, 2000.
- N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — Chapter 1（残差・後退誤差・前進誤差）と Chapter 7（線形方程式の摂動論）。
- E. Hairer, C. Lubich, G. Wanner, *Geometric Numerical Integration*, 2nd ed., Springer, 2006 — 保存量と数値解法。
- F. Johansson, "Arb: efficient arbitrary-precision midpoint-radius interval arithmetic", *IEEE Transactions on Computers* 66 (2017), 1281–1292. [doi:10.1109/TC.2017.2690633](https://doi.org/10.1109/TC.2017.2690633)
- W. Tucker, *Validated Numerics: A Short Introduction to Rigorous Computations*, Princeton University Press, 2011.