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

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$ のことがあります。数値を証明に使うときは、包含の保証がある計算（区間演算・球算術）に置き換えます。

この記事で使う事前知識

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

1. 動機

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

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

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

2. 何を確かめるのか

定義 2.1（検証と妥当性確認）↗

数値計算の正しさには 2 つの別の問いがある。

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

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

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

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

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

定義 2.2（数値的な証拠・反例・証明）↗

命題 $P$ についての数値計算の結果を、次の 3 つに区別する。

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

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

3.1. 観測次数

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

定義 3.1（観測次数）↗

刻み $h$ と $h/2$ での誤差 $e(h), e(h/2)$ がともに $0$ でなく同符号のとき、

$$\hat p(h) = \log_2\frac{e(h)}{e(h/2)}$$

を観測次数という。

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

補題 3.2（対数の評価）✓

$|t| \le 1/2$ ならば $|\ln(1 + t)| \le 2|t|$。

証明（補題 3.2）

$t = 0$ なら両辺は $0$ です。$t \ne 0$ のとき、関数 $s \mapsto \ln(1 + s)$ は $[-1/2, 1/2]$ で微分可能で導関数は $1/(1 + s)$ なので、$0$ と $t$ を端点とする閉区間に平均値の定理（定理 3.3［平均値の定理とテイラーの定理］）を使うと、$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|$ です。

∎

定理 3.3（観測次数の収束）✓

$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$ である。

証明（定理 3.3）

$\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}$$

です。補題 3.2 から $|\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$ に収束します。

∎

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

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

例 3.4（正しい Simpson 則、重みを 1 つ誤った Simpson 則、滑らかでない被積分関数）

複合 Simpson 則の誤差は、被積分関数が $C^4$ 級なら $O(h^4)$ です（系 4.3［数値積分 I］ の 2）。正しい実装と、重みを設定するスライスの終わりを 1 つ取り違えて終点の重みを $2$ にしてしまった実装を、$\int_0^1 e^x\,dx = e - 1$ で比べます。あわせて、正しい実装で $\int_0^1\sqrt x\,dx = 2/3$ を計算します。

# 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))

正しい実装, 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$ に収束しています（定理 3.3）。(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 つあります。

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$ に近づきます。長い条件式は、高次の項が先頭の項に比べて十分小さい、という仮定を書き下したものです。

定理 3.5（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$$ の精度で推定できる。

証明（定理 3.5）

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$ で、定理 3.3 の $\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 外挿の定理（定理 4.2［数値積分 II］）の仮定（定数 $c = C$、剰余の定数 $K$）を満たしています。2 つの刻みによる誤差の推定（命題 6.1［数値積分 II］）で $h$ を $h/2$ に置き換えると、主張の不等式そのものです。

∎

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

例 3.6（3 つの刻みからの次数と誤差の推定、丸めによる次数の崩れ）

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

# 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}")

台形則
    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$ に収束し、理論と合います（定理 3.5 の 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}$ を底にして、それより細かい刻みでは増えています。次の小節で説明します。

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

系 3.7（丸めの雑音があるときの観測次数）✓

定理 3.3 の $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$ を満たす。

証明（系 3.7）

$r(h') + \delta(h')$ を改めて $r(h')$ とみなして 定理 3.3 の前半を使います。$|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$ 以下です。

∎

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

中心差分はその典型です。打ち切り誤差は $\frac{M_3}{6}h^2$、丸めによる誤差は $\varepsilon M_0/h$ の程度で（命題 7.1［数値積分 II］）、$\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}$ を上回ります。表で次数が崩れ始める位置と合っています。次数を観測する実験は、丸めが見え始める手前で止めるか、作業精度を上げて行います。

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

$f=e^x$、$x=1$ の中心差分。$h$ を $2^{-26}$ まで小さくすると、誤差は打ち切り誤差の線に沿って下がったあと、丸めの雑音で下げ止まり、観測次数は $2$ から崩れる。

3.4. 製造解の方法

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

定義 3.8（製造解の方法）↗

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

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

例 3.9（製造解で 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 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)))

正しい差分
   最大誤差   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}$ の手計算を誤ると、正しいコードが不合格になるからです。

4. 残差は誤差ではない

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

命題 4.1（残差と誤差）✓

$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$ である。

証明（命題 4.1）

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}$ です。誘導ノルムの定義（定義 4.1［行列のノルムと線形方程式の条件数］）から、任意の行列 $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$ は相対誤差に一致します。

∎

$\epsilon = 10^{-10}$ なら、相対残差が $10^{-10}$ でも解は $100\%$ ずれています。残差が教えてくれるのは「$\boldsymbol{y}$ は、少しだけ違う問題の厳密解である」ということ、つまり後退誤差（問題のほうをどれだけ動かせば、この $\boldsymbol{y}$ が厳密解になるか）が小さいことです。これに対し前進誤差は、答えそのもののずれ $\boldsymbol{y} - \boldsymbol{x}$ のことです（残差と後退誤差の関係は Rigal–Gaches の公式 定理 7.2［行列のノルムと線形方程式の条件数］）。それを前進誤差に直すには条件数が要ります（系 7.3［行列のノルムと線形方程式の条件数］、非線形の問題では 定理 4.3［条件数と安定性］）。残差による検算は、条件数の見積もりと組にして初めて誤差の主張になります。

残差がそのまま誤差の上界になる場合もあります。実対称行列の固有値では、単位ベクトル $\boldsymbol{x}$ と Rayleigh 商の残差 $\boldsymbol{r}$ について、$\|\boldsymbol{r}\|$ 以内に必ず固有値があります（系 7.2［固有値の計算 I］）。条件数が現れないのは、対称行列の固有値が摂動に対してつねに良条件だからです。対称でない行列ではこの評価は成り立ちません。

例 4.2（固有値計算の検算を重ねる）

乱数で作った $8$ 次の実対称行列の固有値を計算し、5 通りの検算をします。

# 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}")

トレース      |Σλ - 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 つずつ見ますが、系 7.2［固有値の計算 I］ が保証するのは「近くに固有値がある」ことで、8 個の固有値を漏れなく求めたことではありません。同じ固有値を 2 回数えていないことは直交性が補います。別の算法との一致は最も強い検査に見えますが、2 つの算法は同じライブラリの同じ入力 A を使っているので、A を作る段階の誤りには無力です（§6）。

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

5.1. 保存量のずれ

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

命題 5.1（保存量のずれの次数）✓

$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}$$

である。

証明（命題 5.1）

$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$ です。

∎

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

例 5.2（調和振動子のエネルギーのずれ）

$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 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}")

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$ に近づいています。命題 5.1 の $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 つを取り違えないようにします。

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

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

命題 5.3（不変性の検査が検出できる誤り）✓

$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)$ で、正しい値とのどんな比較も合格する。

証明（命題 5.3）

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)$ です。

∎

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

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

例 5.4（転置を取り違えた行列ベクトル積はどの検査で見つかるか）

行列ベクトル積 $A\boldsymbol{x}$ を、添字を取り違えて $A^{\mathsf{T}}\boldsymbol{x}$ と計算している実装を、5 通りの検査にかけます。

# 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 '見逃し'}")

線形性 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 つの恒等式を保つ誤りです（命題 5.3 の 1 に当たる状況）。対称行列での比較が合格するのは同じ命題の 2（$\sigma$ が転置、$S^{\mathsf{T}} = S$）です。検出できたのは、対称でない入力での別手法との比較と、手で答えの分かる最小の入力だけでした。乱数で作った行列は確率 1 で対称でないので、検査の入力を乱数で作ることには、入力の対称性による見逃しを避ける効果があります。

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

6.1. 2 つの検算は独立か

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

命題 6.1（共通の原因があると見逃しの確率は積にならない）✓

ある誤りを検算 $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$$ である。

証明（命題 6.1）

1. 事象の独立性の定義（定義 6.1［確率空間とコルモゴロフの公理］）そのものです。

2. 全確率の公式（定理 5.3［確率空間とコルモゴロフの公理］ の 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$ です。

∎

数値を入れてみます。それぞれ単独では $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. 反例が見つからなかったとき

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

定理 6.2（無作為な反例探索で言える範囲）✓

入力の集合の上の確率分布を 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$ 以上で反例が出ないまま検査を通過する。

証明（定理 6.2）

1. 反例の集合を $B$ とすると、反例が出ないことは $X_1 \notin B, \ldots, X_N \notin B$ で、独立性（定義 2.3［確率変数と期待値］）から確率は $\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)$ で広義単調増加で（系 3.4［平均値の定理とテイラーの定理］ の 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$ です。

∎

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$ であっても、命題が真であるとは限りません。

例 6.3（無作為な探索が見逃す反例）

測度 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$ は直接の計算が届かないほど大きいところにあります。計算できる範囲での一致が、命題の証拠としてどれほど弱いことがあるかを示す例です。

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

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

数値実験で見つけた事実を、論文や報告の中で「成り立つ」と書くまでには、いくつかの段があります。

(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$ の根拠（定理 3.5 の誤差の推定、条件数と残差、誤差予算）を示します。推定値と保証された上界は区別します。Richardson の推定は漸近的に正しい推定値で、上界ではありません。

(iii) 探索の範囲・件数・種・失敗例を記録する。 定理 6.2 のとおり、反例が出なかったという報告の価値は、分布と件数で決まります。途中で止まった試行、収束しなかった試行も消さずに件数を残します（長時間ジョブの管理の章の打ち切りの偏り）。

(iv) 証明に使う数値は、包含の保証がある計算に置き換える。 浮動小数の計算結果は、真の値の近くにあると期待できるだけで、どちら側にどれだけずれているかを教えてくれません。証明の中で「この量は正である」「この区間に零点がある」という事実を使うなら、結果が真の値を含む区間であることが保証された計算を使います。区間演算では、算術式の値が計算された区間に必ず含まれます（定理 4.4［浮動小数と区間演算］）。球算術でも同じです（定理 4.9［浮動小数と区間演算］）。区間が $0$ を含まなければ符号が確定し、区間 Newton 法の条件が満たされれば零点の存在と一意性が確定します（定理 5.2［浮動小数と区間演算］）。

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

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

例 7.1（倍精度で正に見える量が、実は 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 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)}")

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$ の程度）は上の説明と合います。

この例の球算術の保証について、正確に述べておきます。Rikai で証明してあるのは、球の加算と乗算が真の値を含むこと（定理 4.9［浮動小数と区間演算］）です。平方根・指数関数・円周率の球が真の値を含むことは、ライブラリ（Arb）の実装が保証する性質で、ここではそれを信頼しています。計算機を使った証明では、「何を信頼したか」（ライブラリ、その版、コンパイラ、ハードウェアの演算）を明記するのが作法です。この先の、数値を証明に変える技法は、浮動小数点数と区間演算の章から始まる精度保証つき計算のテーマで扱います。

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

8. 演習

演習 8.1易

ある求積のコードの誤差が、刻み $h, h/2, h/4$ でそれぞれ $3.2 \times 10^{-3}$、$8.1 \times 10^{-4}$、$2.0 \times 10^{-4}$ でした。観測次数を求めてください。このコードが複合 Simpson 則のつもりで書かれていたとすると、何を疑いますか。

解答

$\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 つは仮定で、被積分関数の高階の導関数が積分区間で有界でない場合です（例 3.4 の $\sqrt x$ は次数 $1.5$ でした）。滑らかな関数（$e^x$ など）で同じコードを試し、それでも $2$ なら実装の誤りです。

演習 8.2標準

真値の分からない量について、刻み $h, h/2, h/4$ での計算値が $A(h) = 1.2000$、$A(h/2) = 1.0500$、$A(h/4) = 1.0125$ でした。定理 3.5 の方法で次数を推定し、$A(h/4)$ の誤差の推定値と、それを使って補正した値を求めてください。

解答

$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$ に近いことを確かめるのが安全です。

演習 8.3標準

ある不等式を、種を固定した $10^6$ 個の一様乱数の入力で調べ、反例は見つかりませんでした。(i) 信頼水準 $95\%$ で、反例の割合について何が言えますか。(ii) 信頼水準 $99\%$ で「反例の割合は $10^{-9}$ 以下」と言うには、何個の標本が必要ですか。(iii) (ii) の個数を調べて反例がなかったとして、不等式が証明されたと言えない理由を 2 つ挙げてください。

解答

(i) 定理 6.2 の 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$ の反例、標本の分布の外の反例。例 6.3）。第二に、各標本での判定は浮動小数の計算によるので、不等式の両辺の差が丸め誤差の程度の入力では、判定そのものが信頼できません（例 7.1）。等号に近い場合がある不等式ほど、この問題は深刻です。証明にするには、入力の集合を有限個の区間に分け、各区間の上で区間演算により不等式を確かめる、といった方法をとります。

演習 8.4標準

例 5.2 の調和振動子に 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$ のときの値を求め、例の表と比べてください。

解答

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 次です（命題 5.1 の $p = 1$）。Euler 法のエネルギーは増える一方で、時間について指数関数的に増えます。命題の上界が $T$ について線形なのは、数値解が有界な集合 $U$ にとどまる範囲での評価だからです。

参考文献

• 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
• W. Tucker, Validated Numerics: A Short Introduction to Rigorous Computations, Princeton University Press, 2011.