数値積分 II：Euler–Maclaurin 展開・周期関数の台形則・Richardson 外挿・Romberg 法

0. この記事の要点

曲線の下の面積を、細い台形を並べて近似する方法があります。台形を細かくすれば誤差は減りますが、どのくらい減るのか、減り方の形が分かれば、粗い計算を 2 つ組み合わせて誤差を打ち消せるのではないか、というのがこの記事の出発点です。読み終えると、少ない計算で積分の精度を大きく上げる方法（外挿）と、得られた値の誤差を見積もる方法が分かり、あわせて刻みを小さくしすぎると逆効果になる理由も説明できるようになります。

• 複合台形則の誤差は、$f$ が滑らかなら $T(h) - I = c_1h^2 + c_2h^4 + \cdots + c_ph^{2p} + O(h^{2p+2})$ と刻みの偶数冪で展開できます。係数は Bernoulli 数（§2 で使う特別な数列）と端点での奇数階導関数の差で、証明はEuler–Maclaurin の公式です。
• 周期関数を 1 周期にわたって積分するときは端点の項が打ち消し合い、台形則の誤差は $O(h^{2p})$（$C^{2p}$ 級）、帯状領域で正則なら $4\pi M/(e^{aN} - 1)$ 以下と指数的に小さくなります。証明は Fourier 係数の折り返し（高い周波数の波が、少ない点では低い周波数の波と同じに見えてしまうこと）と積分路の平行移動です。
• Richardson 外挿は、誤差の展開の主要項を 2 つの刻みの値の組み合わせで消します。これを繰り返す Romberg 法では、表の第 $j$ 列の誤差が $O(h^{2j+2})$ になり、第 1 列は複合 Simpson 則に一致します。
• 2 つの刻みの値の差から誤差そのものを推定でき、適応型求積の停止判定の根拠になります。仮定（滑らかさ）が崩れると、加速も推定も効かなくなります（$\int_0^1\sqrt x\,dx$）。
• 数値微分では刻みを小さくすると丸め誤差が $u/h$ で増えるので、中心差分の最適な刻みは $h \approx u^{1/3}$ の程度で、得られる精度は $u^{2/3}$ 程度です。

この記事で使う事前知識

• 複合台形則・求積則: 区間を等分して台形の面積を足し、積分の近似値とする方法。積分を関数値の重み付き和で近似する公式を求積則と呼びます。→ 求積則
• 漸近展開: $h \to 0$ での振る舞いを $h$ の冪の和で表す展開。和が収束するとは限らず、$h$ が小さいときに最初の数項が精度を与えます。→ 漸近解析
• 正則関数: 複素平面の領域で複素微分できる関数。→ 正則関数
• Fourier 係数: 周期関数を $e^{in\theta}$ の重ね合わせで表したときの、各 $e^{in\theta}$ の係数 $\hat f(n)$。
• 単位丸め $u$: 浮動小数の 1 回の四則演算で生じる相対誤差の上限（倍精度で $2^{-53}$）。

1. 動機

Newton–Cotes 則と Gauss 求積の章で、複合台形則の誤差は $-\frac{(b-a)h^2}{12}f''(\xi)$ であることを示しました。この式は誤差の大きさを教えてくれますが、$\xi$ が未知なので誤差を消すことはできません。ところが、誤差の式をもっと精密に、

$$T(h) = I + c_1h^2 + c_2h^4 + \cdots$$

と $h$ のべき級数の形で書けるなら話は変わります。$c_1$ が未知でも、$T(h)$ と $T(h/2)$ の 2 つの値から $c_1h^2$ の項を消去できるからです。この「誤差の展開の形だけを使って、係数を知らずに誤差を消す」技法が Richardson 外挿で、数値積分だけでなく、数値微分、常微分方程式の数値解、格子の上の計算の連続極限など、刻みのある計算すべてに使えます。

本章の主役は Euler–Maclaurin の公式です。これが台形則の誤差の展開を与え、外挿が効く理由と、効かない場合（端点で滑らかでない関数）の理由の両方を説明します。さらに同じ公式から、周期関数では台形則が驚くほど正確になるという、一見すると逆説的な事実も出てきます。最後に、刻みを小さくすると逆に精度が落ちる数値微分を例に、打ち切り誤差と丸め誤差の釣り合いを定量的に扱います。

2. 複合台形則の漸近展開

漸近解析の章の Bernoulli 多項式 $B_m(t)$ と Bernoulli 数 $B_m = B_m(0)$（補題 7.1［漸近解析］）、その周期 1 の延長 $\tilde B_m(t) = B_m(t - \lfloor t\rfloor)$ を使います。$B_2 = 1/6$、$B_4 = -1/30$、$B_6 = 1/42$ です。$B_m$ は $[0,1]$ 上の連続関数なので最大値 $\beta_m = \max_{0 \le t \le 1}|B_m(t)|$ をもち（定理 2.3［平均値の定理とテイラーの定理］）、すべての実数 $t$ で $|\tilde B_m(t)| \le \beta_m$ です。

以下 $a \lt b$、$N \ge 1$、$h = (b-a)/N$ とし、複合台形則を $T(h) = h\bigl(\frac12f(a) + \sum_{k=1}^{N-1}f(a + kh) + \frac12f(b)\bigr)$、$I = \int_a^b f(x)\,dx$ と書きます。

次の定理のねらいは、誤差 $T(h) - I$ が $h^2, h^4, \ldots$ の項の和として書け、その係数が両端での奇数階導関数の差だけで決まることを示すことです。Euler–Maclaurin の公式そのものは証明の最初の式で、和 $\sum g(k)$ を、積分、両端の平均、Bernoulli 数を含む端点の項、剰余の積分で書き表したものです。証明では、これを $g(s) = f(a + sh)$ に当てはめて整理します。

定理 2.1（複合台形則の漸近展開）✓

$p \ge 0$ とし、$f \in C^{2p+2}[a,b]$ とする。

$$c_j = \frac{B_{2j}}{(2j)!}\bigl(f^{(2j-1)}(b) - f^{(2j-1)}(a)\bigr) \qquad (j \ge 1)$$

とおくと、すべての $N \ge 1$ について

$$T(h) - I = \sum_{j=1}^{p} c_j h^{2j} + R_p(h), \qquad |R_p(h)| \le K_p h^{2p+2}$$

が成り立つ。ここで $K_p = |c_{p+1}| + \frac{(b-a)\beta_{2p+2}}{(2p+2)!}\max_{[a,b]}|f^{(2p+2)}|$ は $N$ によらない定数である。

証明（定理 2.1）

$g(s) = f(a + sh)$（$s \in [0, N]$）とおくと $g \in C^{2p+2}[0,N]$ で、$g^{(m)}(s) = h^m f^{(m)}(a + sh)$ です。Euler–Maclaurin の公式（定理 7.2［漸近解析］。整数 $0 \lt N$、項の数 $p+1$）を $g$ に使うと

$$\sum_{k=0}^{N}g(k) = \int_0^N g(s)\,ds + \frac{g(0) + g(N)}{2} + \sum_{j=1}^{p+1}\frac{B_{2j}}{(2j)!}\bigl(g^{(2j-1)}(N) - g^{(2j-1)}(0)\bigr) - \int_0^N\frac{\tilde B_{2p+2}(s)}{(2p+2)!}g^{(2p+2)}(s)\,ds$$

です。$T(h) = h\bigl(\sum_{k=0}^N g(k) - \frac{g(0) + g(N)}{2}\bigr)$、置換積分（定理 6.1［積分の基本定理と定積分］）から $h\int_0^N g(s)\,ds = I$ なので、両辺に $h$ を掛けて

$$T(h) - I = \sum_{j=1}^{p+1}c_jh^{2j} - h^{2p+3}\int_0^N\frac{\tilde B_{2p+2}(s)}{(2p+2)!}f^{(2p+2)}(a + sh)\,ds$$

を得ます。$j = p+1$ の項と最後の積分を $R_p(h)$ とします。積分の単調性（命題 4.1［積分の基本定理と定積分］ の 2）と $|\tilde B_{2p+2}| \le \beta_{2p+2}$ から、最後の項の絶対値は $h^{2p+3}N\frac{\beta_{2p+2}}{(2p+2)!}\max|f^{(2p+2)}|$ 以下で、$Nh = b - a$ なので主張を得ます。

∎

$p = 0$ の場合は $|T(h) - I| \le K_0h^2$ で、前の章の複合台形則の評価（系 4.3［数値積分 I］）と同じ次数です。展開の価値は、誤差が偶数冪だけで、しかも係数 $c_j$ が $h$ によらないことにあります。注意すべきは、これは収束する級数ではなく漸近展開だということです。$p$ を増やすと $K_p$ は一般に急激に大きくなり、$p \to \infty$ で右辺の級数が収束するとは限りません。外挿では、固定した $p$ について $h \to 0$ の振る舞いだけを使います。

最初の係数は $c_1 = \frac{1}{12}\bigl(f'(b) - f'(a)\bigr)$ です。$f'(b) = f'(a)$ なら $h^2$ の項が消え、さらに奇数階導関数がすべて両端で一致すれば、どの $c_j$ も消えます。これが次の節の出発点です。

3. 周期関数の台形則

$f$ を周期 $2\pi$ の関数とし、$h = 2\pi/N$、

$$T_N(f) = \frac{2\pi}{N}\sum_{k=0}^{N-1}f\Bigl(\frac{2\pi k}{N}\Bigr), \qquad I(f) = \int_0^{2\pi}f(\theta)\,d\theta$$

とおきます。$f(0) = f(2\pi)$ なので、$T_N$ は $[0, 2\pi]$ の複合台形則そのものです。

次の補題のねらいは、Fourier 係数が $n$ とともに指数的に小さくなることを示すことです。方針は、実軸上の積分を、虚軸方向に $b$ だけずらした直線上の積分に取り替えることです。そこでは $e^{-in\theta}$ の大きさが $e^{-nb}$ に縮むので、係数も同じ割合で小さいと分かります。ずらしてよい理由は、周期性と Cauchy の積分定理（Goursat の定理）で、長方形の周に沿った積分が $0$ になることです。

補題 3.1（帯状領域で正則な周期関数の Fourier 係数）✓

$a \gt 0$ とし、$f$ は帯状領域 $S_a = \{z \in \mathbb{C} : |\mathrm{Im}\,z| \lt a\}$ で正則、$f(z + 2\pi) = f(z)$、$S_a$ 上で $|f| \le M$ とする。このとき $0 \lt b \lt a$ とすべての整数 $n$ について、実軸上の $f$ の Fourier 係数 $\hat f(n) = \frac1{2\pi}\int_{-\pi}^{\pi}f(\theta)e^{-in\theta}\,d\theta$ は

$$|\hat f(n)| \le M e^{-|n|b}$$

を満たす。

証明（補題 3.1）

$n \ge 0$ とします（$n \lt 0$ は下の $-b$ を $b$ に替えて同じです）。$g(z) = f(z)e^{-inz}$ は $S_a$ で正則で、$n$ が整数なので $g(z + 2\pi) = g(z)$ です。頂点 $-\pi - ib$、$\pi - ib$、$\pi$、$-\pi$ の長方形 $R$（周と内部、$S_a$ に含まれる）の周を正の向きに 1 周する積分を考えます。$R$ を対角線で 2 つの閉三角形に分けると、Goursat の定理（定理 4.5［コーシーの積分定理と積分公式］）から各三角形の周の積分は $0$ で、共通の対角線は逆向きに 2 回通るので打ち消し合います（命題 2.3［コーシーの積分定理と積分公式］ の 2・3）。よって $\oint_{\partial R}g\,dz = 0$ です。

周を 4 辺に分けます。右辺 $\pi - ib \to \pi$ の積分は $i\int_0^b g(\pi - ib + it)\,dt = i\int_0^b g(\pi - is)\,ds$（$s = b - t$）、左辺 $-\pi \to -\pi - ib$ の積分は $-i\int_0^b g(-\pi - is)\,ds = -i\int_0^b g(\pi - is)\,ds$（周期性）なので、両者は打ち消し合います。残る上辺（実軸上を $\pi$ から $-\pi$ へ）と下辺（$-\pi - ib$ から $\pi - ib$ へ）から

$$\int_{-\pi}^{\pi}g(\theta)\,d\theta = \int_{-\pi}^{\pi}g(\theta - ib)\,d\theta$$

です。左辺は $2\pi\hat f(n)$ です。右辺の被積分関数は $|f(\theta - ib)||e^{-in(\theta - ib)}| \le Me^{-nb}$ を満たすので、ML 不等式（補題 2.4［コーシーの積分定理と積分公式］。長さ $2\pi$）から $2\pi|\hat f(n)| \le 2\pi Me^{-nb}$ です。

∎

次の定理の 1 は、前節の展開で端点の項が周期性により消えることから従います。2 の方針は、$f$ を Fourier 級数で書いて台形則の和に代入し、$k$ についての和が、$n$ が $N$ の倍数のときだけ $N$、それ以外では $0$ になることを使う、というものです。残るのは周波数 $mN$ の係数だけで、補題がそれらの小ささを与えます。

定理 3.2（周期関数の台形則）✓

$f : \mathbb{R} \to \mathbb{C}$ を周期 $2\pi$ の関数とする。

1. $p \ge 1$ とし、$f$ が $C^{2p}$ 級ならば、すべての $N \ge 1$ について $$|T_N(f) - I(f)| \le \frac{2\pi\beta_{2p}}{(2p)!}\max_{[0,2\pi]}|f^{(2p)}|\Bigl(\frac{2\pi}{N}\Bigr)^{2p}.$$
2. $f$ が 補題 3.1 の仮定（帯状領域 $S_a$ で正則、周期 $2\pi$、$|f| \le M$）を満たすならば、すべての $N \ge 1$ について $$T_N(f) - I(f) = 2\pi\sum_{m \ne 0}\hat f(mN), \qquad |T_N(f) - I(f)| \le \frac{4\pi M}{e^{aN} - 1}.$$

証明（定理 3.2）

1. 定理 2.1 の証明と同じく $g(s) = f(sh)$（$h = 2\pi/N$）に Euler–Maclaurin の公式（定理 7.2［漸近解析］。項の数 $p$）を使います。周期性から $g^{(m)}(N) = h^mf^{(m)}(2\pi) = h^mf^{(m)}(0) = g^{(m)}(0)$ なので、Bernoulli 数の項はすべて消え、$g(N) = g(0)$ から $T_N(f) = h\bigl(\sum_{k=0}^N g(k) - \frac{g(0)+g(N)}{2}\bigr)$ です。よって

$$T_N(f) - I(f) = -h^{2p+1}\int_0^N\frac{\tilde B_{2p}(s)}{(2p)!}f^{(2p)}(sh)\,ds$$

で、絶対値は $h^{2p+1}N\beta_{2p}\max|f^{(2p)}|/(2p)! = 2\pi h^{2p}\beta_{2p}\max|f^{(2p)}|/(2p)!$ 以下です（積分の単調性、命題 4.1［積分の基本定理と定積分］ の 2。$f$ が複素数値のときは、$|\alpha| = 1$ を $\alpha\int\varphi = |\int\varphi|$ となるように選べば $|\int\varphi| = \int\mathrm{Re}(\alpha\varphi) \le \int|\varphi|$ となり、同じ評価が成り立ちます）。

2. 正則関数は無限回微分可能なので（系 3.3［正則関数の強力な性質］）、実軸に制限した $f$ は $C^1$ 級の周期関数です。命題 5.4［フーリエ級数］ から $\sum_n|\hat f(n)| \lt \infty$ で、Fourier 級数は $f$ に一様に（とくに各点で）収束します。有限和と収束級数の和は交換できるので

$$T_N(f) = \frac{2\pi}{N}\sum_{k=0}^{N-1}\sum_{n\in\mathbb{Z}}\hat f(n)e^{2\pi ink/N} = 2\pi\sum_{n\in\mathbb{Z}}\hat f(n)\cdot\frac1N\sum_{k=0}^{N-1}e^{2\pi ink/N}$$

です。内側の和は $n$ が $N$ の倍数なら $N$、そうでなければ公比 $e^{2\pi in/N} \ne 1$ の等比級数の和 $\frac{1 - e^{2\pi in}}{1 - e^{2\pi in/N}} = 0$ です。よって $T_N(f) = 2\pi\sum_m\hat f(mN)$ です。一方、周期性から $I(f) = \int_{-\pi}^{\pi}f(\theta)\,d\theta = 2\pi\hat f(0)$ なので、等式を得ます。

$0 \lt b \lt a$ について 補題 3.1 を使うと

$$|T_N(f) - I(f)| \le 2\pi\cdot 2\sum_{m=1}^{\infty}Me^{-mNb} = \frac{4\pi M}{e^{Nb} - 1}$$

です。右辺は $b$ について連続なので、$b \to a$ とすれば主張を得ます。

∎

1 は「滑らかさの次数だけ速く収束する」こと、2 は「解析的なら指数的に速く収束する」ことを述べています。1 では $p$ はいくらでも大きく取れるので、$C^\infty$ 級の周期関数ではどんな多項式の速さよりも速く誤差が減ります。等式 $T_N - I = 2\pi\sum_{m\ne0}\hat f(mN)$ は、$N$ 点の等間隔標本では周波数 $mN$ の成分が周波数 $0$（定数）と区別できない、という**折り返し（エイリアシング）**の表現です。積分の誤差は、折り返されて定数成分に混ざった高周波成分の大きさにほかなりません。

例 3.3（e^{cos θ} の台形則の指数収束）

$\int_0^{2\pi}e^{\cos\theta}d\theta = 2\pi I_0(1)$（$I_0$ は第 1 種変形 Bessel 関数）です。$f(z) = e^{\cos z}$ は整関数で、$z = x + iy$ について $|e^{\cos z}| = e^{\cos x\cosh y} \le e^{\cosh a}$ なので、任意の $a \gt 0$ で 定理 3.2 の 2 が $M = e^{\cosh a}$ として使えます。$f$ の Fourier 係数は $\hat f(n) = I_n(1)$ なので、誤差の主要項は $m = \pm1$ の折り返し $4\pi I_N(1)$ です。

# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0
import math
import numpy as np
from scipy.special import iv
from mpmath import mp, besseli

mp.dps = 30
I_sp = 2 * math.pi * iv(0, 1.0)              # 2π I_0(1)（SciPy）
I_mp = 2 * mp.pi * besseli(0, 1)             # 同じ値を mpmath で
f = lambda t: np.exp(np.cos(t))

def bound(N, a):                             # 4πM/(e^{aN}-1)
    M = math.exp(math.cosh(a))               # |e^{cos z}| ≤ e^{cosh a}
    return 4 * math.pi * M / math.expm1(a * N)

print("  N   誤差       4πI_N(1)   評価(a=1)  評価(a 最適)")
for N in [2, 4, 6, 8, 10, 12]:
    th = 2 * math.pi * np.arange(N) / N
    err = abs(2 * math.pi / N * f(th).sum() - I_sp)
    alias = float(4 * mp.pi * besseli(N, 1))  # 折り返しの主要項
    best = min(bound(N, a) for a in np.linspace(0.1, 6, 600))
    print(f"{N:3d}  {err:.2e}   {alias:.2e}   {bound(N, 1):.2e}"
          f"   {best:.2e}")
print(f"SciPy と mpmath の参照値の差: {abs(I_sp - float(I_mp)):.1e}")

  N   誤差       4πI_N(1)   評価(a=1)  評価(a 最適)
  2  1.74e+00   1.71e+00   9.20e+00   6.92e+00
  4  3.44e-02   3.44e-02   1.10e+00   1.78e-01
  6  2.83e-04   2.83e-04   1.46e-01   1.77e-03
  8  1.25e-06   1.25e-06   1.97e-02   9.00e-06
 10  3.46e-09   3.46e-09   2.67e-03   2.77e-08
 12  6.53e-12   6.53e-12   3.61e-04   5.72e-11
SciPy と mpmath の参照値の差: 0.0e+00

検算. (a) 参照値 $2\pi I_0(1)$ を SciPy の iv と mpmath の besseli（30 桁）の 2 通りで計算し、倍精度で一致しています。(b) 誤差は、証明の等式の主要項 $4\pi I_N(1)$（mpmath で計算）と $N \ge 4$ で表示の桁まで一致します。$N = 2$ でずれるのは、$m = \pm2$ の項 $4\pi I_4(1) \approx 0.03$ が無視できないためです。(c) 誤差はすべて評価 $4\pi M/(e^{aN}-1)$ の内側にあり、$a$ を $N$ ごとに選び直すと評価は真の誤差の 10 倍程度まで近づきます。誤差は $N$ を 2 増やすごとに 2〜3 桁減り、12 点で $10^{-11}$ に達します。整関数なので $a$ をいくらでも大きく取れ、収束は指数関数よりもさらに速くなっています（$I_N(1) \approx (1/2)^N/N!$）。

4. Richardson 外挿

$h$ を $0$ に近づけたときの極限 $A$ を求めたいが、計算できるのは有限の $h$ での近似値 $A(h)$ だけ、という状況を考えます。以下、$A(h)$ は $h \in \{h_0, h_0/2, h_0/4, \ldots\}$ で定義されているとします。

定義 4.1（誤差の漸近展開）↗

$0 \lt p_1 \lt p_2 \lt \cdots \lt p_m \lt q$ とし、$h_k = h_0/2^k$ とする。定数 $c_1, \ldots, c_m$ と $C \ge 0$ があって、すべての $k \ge 0$ で

$$A(h_k) = A + \sum_{i=1}^{m}c_ih_k^{p_i} + r_k, \qquad |r_k| \le Ch_k^q$$

が成り立つとき、$A(h)$ は指数 $p_1, \ldots, p_m$、剰余の次数 $q$ の誤差の漸近展開をもつという。係数 $c_i$ は $k$（すなわち $h$）によらない。

定理 2.1 は、$f \in C^{2p+2}[a,b]$ なら複合台形則が指数 $2, 4, \ldots, 2p$、剰余の次数 $2p+2$ の誤差の漸近展開をもつ、と言い換えられます。大切なのは指数が既知で、係数が未知でもよいことです。

定理 4.2（Richardson 外挿）✓

$0 \lt p \lt q$ とし、定数 $c$、$C \ge 0$ があって、ある $h$ と $h/2$ について

$$A(h) = A + ch^p + r(h), \qquad A(h/2) = A + c(h/2)^p + r(h/2), \qquad |r(h)| \le Ch^q,\ |r(h/2)| \le C(h/2)^q$$

が成り立つとする。このとき

$$A_1(h) = \frac{2^pA(h/2) - A(h)}{2^p - 1}$$

は $|A_1(h) - A| \le \dfrac{1 + 2^{p-q}}{2^p - 1}\,Ch^q$ を満たす。

証明（定理 4.2）

定義に代入すると、$c$ の項は $2^pc(h/2)^p - ch^p = 0$ で消え、

$$A_1(h) - A = \frac{2^pr(h/2) - r(h)}{2^p - 1}$$

です。$|2^pr(h/2)| \le 2^pC2^{-q}h^q$、$|r(h)| \le Ch^q$ から主張を得ます。

∎

外挿の要点は、$c$ を知らなくてよいことです。必要なのは誤差の次数 $p$ だけです。次数を誤って仮定すると、主要項は消えずに残ります。外挿を繰り返すと、展開の項を順に消していけます。

次の系は、これを繰り返したときに何が消え、何が残るかを記録するものです。$j$ 回外挿した値 $T_{k,j}$ では、$h_k^{p_1}, \ldots, h_k^{p_j}$ の項が消えて $h_k^{p_{j+1}}$ 以降が残ります。式の中の $\gamma_{i,j}$ と積は、残った項の係数と剰余の大きさが、外挿のたびにどう変わるかを書いただけの記号です。

系 4.3（外挿表の誤差）✓

$A(h)$ が指数 $p_1 \lt \cdots \lt p_m$、剰余の次数 $q$ の誤差の漸近展開（定義 4.1。記号もそこでのもの）をもつとする。外挿表を $T_{k,0} = A(h_k)$、$1 \le j \le \min(k, m)$ について

$$T_{k,j} = \frac{2^{p_j}T_{k,j-1} - T_{k-1,j-1}}{2^{p_j} - 1}$$

で定めると、

$$T_{k,j} = A + \sum_{i=j+1}^{m}c_i\gamma_{i,j}\,h_k^{p_i} + r_{k,j}, \qquad \gamma_{i,j} = \prod_{l=1}^{j}\frac{2^{p_l} - 2^{p_i}}{2^{p_l} - 1}, \qquad |r_{k,j}| \le C\prod_{l=1}^{j}\frac{2^{p_l} + 2^q}{2^{p_l} - 1}\,h_k^q$$

が成り立つ。とくに $T_{k,j} - A = O(h_k^{p_{j+1}})$（$j \lt m$）、$T_{k,m} - A = O(h_k^q)$ である（$k \to \infty$。定数は $k$ によらない）。

証明（系 4.3）

$j$ についての帰納法です。$j = 0$ は仮定です（空の積は $1$）。$j - 1$ で成り立つとし、$k \ge j$ とします。$h_{k-1} = 2h_k$ なので、$T_{k-1,j-1}$ の展開の $h_{k-1}^{p_i}$ は $2^{p_i}h_k^{p_i}$ です。漸化式に代入すると、$h_k^{p_i}$ の係数は

$$c_i\gamma_{i,j-1}\cdot\frac{2^{p_j} - 2^{p_i}}{2^{p_j} - 1} = c_i\gamma_{i,j}$$

で、$i = j$ では $0$ になります。剰余は $r_{k,j} = (2^{p_j}r_{k,j-1} - r_{k-1,j-1})/(2^{p_j} - 1)$ で、帰納法の仮定の定数を $C_{j-1}$ とすると $|r_{k-1,j-1}| \le C_{j-1}h_{k-1}^q = C_{j-1}2^qh_k^q$ なので、$|r_{k,j}| \le C_{j-1}\frac{2^{p_j} + 2^q}{2^{p_j} - 1}h_k^q$ です。最後の主張は、$h_k \le h_0$ から $h_k^{p_i} \le h_0^{p_i - p_{j+1}}h_k^{p_{j+1}}$（$i \ge j+1$）、$h_k^q \le h_0^{q - p_{j+1}}h_k^{p_{j+1}}$ であることによります。

∎

5. Romberg 法

複合台形則に 系 4.3 を $p_j = 2j$ で使うのが Romberg 法です。

定義 5.1（Romberg 表）↗

$N_0 \ge 1$、$H = (b-a)/N_0$、$h_k = H/2^k$ とし、$T_{k,0} = T(h_k)$（複合台形則）、$1 \le j \le k$ について

$$T_{k,j} = \frac{4^jT_{k,j-1} - T_{k-1,j-1}}{4^j - 1}$$

と定める。三角形に並んだ $T_{k,j}$ を Romberg 表、$j$ を固定した並びを第 $j$ 列、$T_{k,k}$ の並びを対角線という。

$T(h_k)$ は、$T(h_{k-1})/2$ に新しい中点での値の和 $h_k\sum f$ を足せば計算でき、それまでの関数値を再利用できます。したがって $T_{K,K}$ までの表全体の計算に要る関数値は $N_02^K + 1$ 個で、最も細かい台形則と同じです。外挿の手間は表の大きさ $K^2/2$ 回の四則だけです。

定理 5.2（Romberg 表の誤差）✓

$m \ge 0$、$f \in C^{2m+2}[a,b]$ とする。$0 \le j \le m$ について、$k$ によらない定数 $D_j$ があって、$k \ge j$ で

$$|T_{k,j} - I| \le D_j\,h_k^{2j+2}$$

が成り立つ。すなわち Romberg 表の第 $j$ 列（$j \le m$）の誤差は $O(h^{2j+2})$ である。

証明（定理 5.2）

定理 2.1（$p = m$）から、すべての $k$ で $T(h_k) = I + \sum_{i=1}^{m}c_ih_k^{2i} + R_m(h_k)$、$|R_m(h_k)| \le K_mh_k^{2m+2}$ です。これは 系 4.3 の仮定（$p_i = 2i$、$q = 2m+2$、$A = I$）そのものなので、$T_{k,j} - I = O(h_k^{2j+2})$（$j \lt m$）、$T_{k,m} - I = O(h_k^{2m+2})$ が、$k$ によらない定数で成り立ちます。

∎

命題 5.3（Romberg 表の第 1 列は複合 Simpson 則）✓

$k \ge 1$ について、$T_{k,1} = \frac{4T(h_k) - T(h_{k-1})}{3}$ は刻み $h_k$ の複合 Simpson 則 $S_{h_k}$ に等しい。

証明（命題 5.3）

$h = h_k$、$2M = (b-a)/h$、$f_i = f(a + ih)$ とします。$T(h) = h\bigl(\frac{f_0}{2} + f_1 + \cdots + f_{2M-1} + \frac{f_{2M}}{2}\bigr)$、$T(2h) = 2h\bigl(\frac{f_0}{2} + f_2 + f_4 + \cdots + f_{2M-2} + \frac{f_{2M}}{2}\bigr)$ なので

$$4T(h) - T(2h) = h\Bigl(f_0 + 4\sum_{i\ \text{奇数}}f_i + 2\sum_{0 \lt i \lt 2M,\ i\ \text{偶数}}f_i + f_{2M}\Bigr)$$

で、$3$ で割ると複合 Simpson 則です。

∎

この命題と 定理 2.1 から、複合 Simpson 則の誤差の主要項は $-4c_2h^4 = \frac{h^4}{180}\bigl(f'''(b) - f'''(a)\bigr)$（$S_h - I$）です。これは 前の章 の評価 $S_h - I = \frac{(b-a)h^4}{180}f^{(4)}(\xi)$（系 4.3［数値積分 I］）と、$f'''(b) - f'''(a) = \int_a^b f^{(4)}$ を通じて整合しています。第 2 列は Boole 則（5 点の Newton–Cotes 則）の複合則になります。

例 5.4（∫₀¹ eˣ dx の Romberg 表）

10 行余りの関数で Romberg 表を作り、各列の誤差の減り方を見ます。$N_0 = 1$ から 6 回刻みを半分にします。

# Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0
import math
import numpy as np
from mpmath import mp, quad, exp as mexp

def romberg(f, a, b, K):
    """T[k][j]: 刻み (b-a)/2^k の台形則から j 回外挿した値"""
    T = [[(b - a) / 2 * (f(a) + f(b))]]
    for k in range(1, K + 1):
        h = (b - a) / 2**k
        mid = f(a + h * np.arange(1, 2**k, 2)).sum()   # 新しい点だけ
        row = [T[-1][0] / 2 + h * mid]
        for j in range(1, k + 1):
            row.append((4**j * row[j - 1] - T[-1][j - 1]) / (4**j - 1))
        T.append(row)
    return T

mp.dps = 30
I = math.e - 1
fmt = lambda x: f"{x:.2e}" if x > 1e-13 else "<1e-13  "
T = romberg(np.exp, 0.0, 1.0, 6)
print(" k   第0列     第1列     第2列     対角")
for k in range(7):
    e = [abs(T[k][j] - I) for j in range(min(k, 2) + 1)]
    s = "  ".join(fmt(x) for x in e).ljust(30)
    print(f"{k:2d}  {s}  {fmt(abs(T[k][k] - I))}")
d = float(abs(I - quad(mexp, [0, 1])))
print(f"e-1 と mpmath.quad の差: {d:.1e}")

 k   第0列     第1列     第2列     対角
 0  1.41e-01                        1.41e-01
 1  3.56e-02  5.79e-04              5.79e-04
 2  8.94e-03  3.70e-05  8.59e-07    8.59e-07
 3  2.24e-03  2.33e-06  1.38e-08    3.35e-10
 4  5.59e-04  1.46e-07  2.16e-10    <1e-13
 5  1.40e-04  9.10e-09  3.39e-12    <1e-13
 6  3.50e-05  5.69e-10  <1e-13      <1e-13
e-1 と mpmath.quad の差: 1.4e-16

検算. (a) 真の値は $e - 1$ で、倍精度の $e - 1$ と mpmath の quad（30 桁）の差は $1.4 \times 10^{-16}$（倍精度の丸めの水準）です。(b) 第 1 列は 命題 5.3 のとおり複合 Simpson 則で、前の章の例の Simpson 則の誤差（$N = 4, 8, \ldots, 64$ で $3.70 \times 10^{-5}, \ldots, 5.69 \times 10^{-10}$）と表示の桁まで一致します。(c) 第 0 列・第 1 列・第 2 列の誤差は $k$ が 1 増えるごとにほぼ $4$、$16$、$64$ 分の 1（たとえば第 2 列で $1.38 \times 10^{-8}/2.16 \times 10^{-10} \approx 64$）になり、定理 5.2 の次数 $h^2$、$h^4$、$h^6$ と合います。対角線は 17 点（$k = 4$）で丸めの水準に達します。

例 5.5（端点で滑らかでない関数では加速しない）

$\int_0^1\sqrt x\,dx = 2/3$ では、$f'(x) = 1/(2\sqrt x)$ が $x = 0$ で発散し、定理 2.1 の仮定 $f \in C^2[0,1]$ すら満たされません。Romberg 表の対角線と、二重指数型（tanh–sinh）求積（注意 6.3）を使う mpmath の quad を比べます。

# Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0
import numpy as np
from mpmath import mp, mpf, quad, sqrt as msqrt

def romberg(f, a, b, K):              # 前の例と同じ Romberg 表
    T = [[(b - a) / 2 * (f(a) + f(b))]]
    for k in range(1, K + 1):
        h = (b - a) / 2**k
        mid = f(a + h * np.arange(1, 2**k, 2)).sum()
        row = [T[-1][0] / 2 + h * mid]
        for j in range(1, k + 1):
            row.append((4**j * row[j - 1] - T[-1][j - 1]) / (4**j - 1))
        T.append(row)
    return T

T = romberg(np.sqrt, 0.0, 1.0, 10)
print("  k   点の数  台形の誤差  対角の誤差  比（台形）")
prev = None
for k in range(0, 11, 2):
    et, ed = abs(T[k][0] - 2 / 3), abs(T[k][k] - 2 / 3)
    r = f"{prev / et:.2f}" if prev else "  - "
    print(f"{k:3d}   {2**k + 1:5d}   {et:.2e}    {ed:.2e}    {r}")
    prev = et
mp.dps = 30
v = quad(msqrt, [0, 1])               # tanh-sinh（二重指数型）
print(f"mpmath.quad（30 桁）の誤差: {float(abs(v - mpf(2) / 3)):.1e}")
print(f"理論の比 4^1.5 = {4**1.5:.2f}")

  k   点の数  台形の誤差  対角の誤差  比（台形）
  0       2   1.67e-01    1.67e-01      -
  2       5   2.34e-02    8.91e-03    7.13
  4      17   3.09e-03    1.07e-03    7.58
  6      65   3.96e-04    1.34e-04    7.79
  8     257   5.01e-05    1.67e-05    7.90
 10    1025   6.30e-06    2.09e-06    7.95
mpmath.quad（30 桁）の誤差: 0.0e+00
理論の比 4^1.5 = 8.00

検算. (a) 真の値は $2/3$ です。(b) mpmath の quad は端点の特異性を変数変換で吸収し、30 桁の計算で $2/3$ に一致します（差は 30 桁の作業精度で $0$）。(c) 台形則の誤差は刻みを $1/4$ にするごとに約 $8 = 4^{1.5}$ 分の 1 で、$O(h^{3/2})$ の速さです。Euler–Maclaurin 展開が $h^2, h^4, \ldots$ を仮定して消そうとしても、実際の主要項は $h^{3/2}$ なので消えず、対角線も同じ $O(h^{3/2})$ の速さでしか減りません（定数が 3 分の 1 程度に縮むだけです）。端点に $x^{\alpha}$ 型の特異性がある関数では、誤差の展開に $h^{1+\alpha}, h^{2+\alpha}, \ldots$ の項が加わることが知られており（I. Navot, J. Math. and Phys. 40 (1961) の Euler–Maclaurin 公式の拡張）、その指数を使って外挿すれば加速を取り戻せます。

注意 5.6（Romberg 法の対角線の収束）

Romberg 表の対角線 $T_{k,k}$ は、$f$ が $[a,b]$ で Riemann 可積分でありさえすれば $I$ に収束することが知られています（F. L. Bauer, H. Rutishauser, E. Stiefel, “New aspects in numerical quadrature”, Proc. Symp. Appl. Math. 15, AMS, 1963）。証明は、$T_{k,k}$ が正の重みをもつ求積則であることと、その重みの性質についての別の議論を要します。参考書としては P. J. Davis, P. Rabinowitz『Methods of Numerical Integration』（2nd ed., 1984）があります。本章ではこの事実を定理として使いません。上の例のように、収束しても速さは $f$ の滑らかさで決まります。

6. 2 つの刻みによる誤差の推定

外挿の式は、誤差を消すだけでなく、誤差の大きさを推定するのにも使えます。研究で数値を報告するとき「何桁正しいか」を示す根拠として最もよく使われる方法です。

命題 6.1（2 つの刻みによる誤差の推定）✓

定理 4.2 の仮定のもとで、

$$\Bigl|\bigl(A - A(h/2)\bigr) - \frac{A(h/2) - A(h)}{2^p - 1}\Bigr| \le \Bigl(2^{-q} + \frac{1 + 2^{-q}}{2^p - 1}\Bigr)Ch^q$$

が成り立つ。とくに $c \ne 0$ で、仮定がすべての十分小さい $h = h_0/2^k$ について同じ $c$、$C$ で成り立つならば、$k \to \infty$ で

$$\frac{\bigl(A(h/2) - A(h)\bigr)/(2^p - 1)}{A - A(h/2)} \to 1$$

である。

証明（命題 6.1）

仮定から $A - A(h/2) = -c(h/2)^p - r(h/2)$、$A(h/2) - A(h) = c h^p(2^{-p} - 1) + r(h/2) - r(h)$ なので

$$\frac{A(h/2) - A(h)}{2^p - 1} = -c\Bigl(\frac h2\Bigr)^p + \frac{r(h/2) - r(h)}{2^p - 1}$$

です。差は $-r(h/2) - \frac{r(h/2) - r(h)}{2^p - 1}$ で、$|r(h/2)| \le C2^{-q}h^q$、$|r(h)| \le Ch^q$ から最初の評価を得ます。後半は、分子と分母がともに $-c(h/2)^p + O(h^q)$ で、$q \gt p$、$c \ne 0$ から比が $1$ に近づくことによります。

∎

台形則（$p = 2$）なら、$\bigl(T(h/2) - T(h)\bigr)/3$ が $T(h/2)$ の誤差の推定値です。これは誤差の上界ではなく近似値で、$h$ が十分小さく、展開の仮定が成り立っているときにだけ信頼できます。推定値を使う前に、刻みを変えた 3 つ以上の値で比 $\bigl(A(h) - A(h/2)\bigr)/\bigl(A(h/2) - A(h/4)\bigr)$ が $2^p$ に近いことを確かめるのがよいと思います。例 5.5 の $\sqrt x$ では、この比が $4$ ではなく $2^{1.5} \approx 2.8$ に近づくので、$p = 2$ を仮定した推定が誤っていることを見抜けます。

注意 6.2（適応型求積と Gauss–Kronrod 則）

実用の自動求積（scipy.integrate.quad が呼び出す QUADPACK など）は、区間ごとに 2 通りの近似値を計算してその差を誤差の推定値とし、推定値の大きい区間だけを細かく分ける適応型の算法です。2 通りの近似値として、$n$ 点の Gauss 則に $n+1$ 点を足して関数値を再利用する Gauss–Kronrod 則（たとえば 7 点と 15 点、10 点と 21 点）がよく使われます。算法と実装は R. Piessens ほか『QUADPACK: A Subroutine Package for Automatic Integration』（Springer, 1983）、Kronrod 則の節点の計算は D. P. Laurie, “Calculation of Gauss–Kronrod quadrature rules”, Math. Comp. 66 (1997) にあります。誤差の推定値は上界ではないので、被積分関数が振動する、特異性がある、などの場合には推定値が小さくても誤差が大きいことがあります。

注意 6.3（二重指数型求積）

$x = \tanh\bigl(\frac\pi2\sinh t\bigr)$ と変数変換すると、$\int_{-1}^1 f(x)\,dx$ は $t \in \mathbb{R}$ 上の積分になり、被積分関数は $|t| \to \infty$ で二重指数的に減衰します。これに刻み $h$ の台形則を使うのが二重指数型（DE、tanh–sinh）求積で、端点に $x^{\alpha}$ 型の特異性があっても、$f$ が区間の内部の近くで解析的なら誤差がほぼ $\exp(-cN/\log N)$ の速さで減ります。仕組みは本章の周期関数の台形則と同じく「無限区間の上で解析的な関数の台形則は指数的に正確」という事実です。高橋秀俊・森正武, “Double exponential formulas for numerical integration”, Publ. RIMS Kyoto Univ. 9 (1974), 721–741 が原論文で、mpmath の quad の既定の方法です。

7. 数値微分：打ち切り誤差と丸め誤差

刻みを小さくするほど打ち切り誤差は減りますが、浮動小数で計算すると、近い 2 つの値の差に含まれる丸め誤差が刻みで割られて拡大されます。中心差分 $D(h) = \frac{f(x+h) - f(x-h)}{2h}$ でこの釣り合いを定量的に調べます。$u$ は単位丸め（倍精度で $2^{-53}$）で、浮動小数の四則は標準モデル（系 3.1［浮動小数点数と丸め誤差］）に従うとします。標準モデルとは、浮動小数の 1 回の四則演算の結果が、真の値の $(1 + \delta)$ 倍（$|\delta| \le u$）になる、という取り決めです。

次の図は、後の例の数値を使って、刻みを小さくしたときの全体の誤差が V 字になる様子を描いたものです。

中心差分の誤差（後の例 x³ の表の値）。横軸は h = 2^-k、縦軸は誤差で、どちらも対数目盛。実線が全体の誤差、破線が打ち切り誤差 h²、点線が丸めの評価。

命題 7.1（中心差分の打ち切り誤差と丸め誤差）✓

$f$ は $[x - h_0, x + h_0]$ を含む開区間で $C^3$ 級とし、$M_j = \max_{[x-h_0, x+h_0]}|f^{(j)}|$（$j = 0, 1, 3$）とおく。$0 \lt h \le h_0$ で $x \pm h$ と $2h$ が浮動小数で正確に表され、計算された関数値が $\hat f_{\pm} = f(x \pm h)(1 + \varepsilon_{\pm})$、$|\varepsilon_{\pm}| \le \varepsilon$ を満たすとする。$\hat D = \mathrm{fl}\bigl(\mathrm{fl}(\hat f_+ - \hat f_-)/(2h)\bigr)$ の計算でオーバーフローとアンダーフローが起きなければ

$$|\hat D - f'(x)| \le \frac{M_3}{6}h^2 + \frac{(1+u)^2\varepsilon M_0}{h} + (2u + u^2)M_1$$

である。さらに $M_3 \gt 0$、$E = (1+u)^2\varepsilon M_0 \gt 0$ のとき、右辺の最初の 2 項の和は $h^* = (3E/M_3)^{1/3}$ で最小になり、最小値は $\frac{3^{2/3}}{2}E^{2/3}M_3^{1/3}$ である。

証明（命題 7.1）

打ち切り誤差. Taylor の定理（定理 5.3［平均値の定理とテイラーの定理］。$a = x$、$n = 2$）から $f(x \pm h) = f(x) \pm hf'(x) + \frac{h^2}{2}f''(x) \pm \frac{h^3}{6}f'''(\xi_\pm)$ となる $\xi_\pm$ があり、

$$D(h) - f'(x) = \frac{h^2}{12}\bigl(f'''(\xi_+) + f'''(\xi_-)\bigr), \qquad |D(h) - f'(x)| \le \frac{M_3}{6}h^2$$

です。

丸め誤差. 標準モデル（系 3.1［浮動小数点数と丸め誤差］）から $|\delta_1|, |\delta_2| \le u$ があって $\hat D = (\hat f_+ - \hat f_-)(1 + \delta_1)(1 + \delta_2)/(2h)$ です（$2h$ は正確なので割り算の丸めは 1 回）。$f_\pm = f(x \pm h)$ と書くと

$$\hat D - D(h) = \frac{(f_+ - f_-)\bigl((1+\delta_1)(1+\delta_2) - 1\bigr) + (f_+\varepsilon_+ - f_-\varepsilon_-)(1+\delta_1)(1+\delta_2)}{2h}$$

です。平均値の定理（定理 3.3［平均値の定理とテイラーの定理］）から $|f_+ - f_-| \le 2hM_1$、また $|(1+\delta_1)(1+\delta_2) - 1| \le 2u + u^2$、$|f_+\varepsilon_+ - f_-\varepsilon_-| \le 2\varepsilon M_0$ なので、$|\hat D - D(h)| \le (2u + u^2)M_1 + (1+u)^2\varepsilon M_0/h$ です。三角不等式で 2 つを合わせます。

最適な刻み. $g(h) = \frac{M_3}{6}h^2 + \frac Eh$ は $h \gt 0$ で微分可能で $g'(h) = \frac{M_3}{3}h - \frac{E}{h^2}$ は $h \lt h^*$ で負、$h \gt h^*$ で正なので、$g$ は $(0, h^*]$ で減少、$[h^*, \infty)$ で増加し（系 3.4［平均値の定理とテイラーの定理］）、$h^*$ で最小です。$E = M_3(h^*)^3/3$ を代入すると $g(h^*) = \frac{M_3}{6}(h^*)^2 + \frac{M_3}{3}(h^*)^2 = \frac{M_3}{2}(h^*)^2 = \frac{3^{2/3}}{2}E^{2/3}M_3^{1/3}$ です。

∎

$\varepsilon$ が $u$ の数倍なら、最適な刻みは $h^* \propto u^{1/3} \approx 5 \times 10^{-6}$ の程度で、得られる精度は $u^{2/3} \approx 10^{-11}$ の程度です。倍精度の 16 桁のうち、中心差分で得られる微分係数は 10〜11 桁にとどまります。前進差分 $\bigl(f(x+h) - f(x)\bigr)/h$ では打ち切り誤差が $O(h)$ なので、同じ議論で $h^* \propto u^{1/2}$、精度は $u^{1/2} \approx 10^{-8}$ になります（演習 8.4）。刻みを小さくしすぎるのは、大きすぎるのと同じくらい有害です。

例 7.2（中心差分の誤差の V 字）

$f(x) = x^3$ を $x = 0.7$（に最も近い倍精度数）で中心差分します。$x^3$ は掛け算 2 回で計算するので $\varepsilon \le 2u + u^2$ で、刻みは $2^{-k}$ にして $x \pm h$、$2h$ を正確にします。この関数では打ち切り誤差がちょうど $h^2$ になるので、mpmath（50 桁）で同じ点の差分を厳密に計算し、全体の誤差を打ち切りと丸めに分けます。

# Python 3.12.10 / mpmath 1.3.0
from mpmath import mp, mpf

mp.dps = 50
f = lambda x: x * x * x              # f(x) = x^3
x, u = 0.7, 2.0**-53                 # x は 0.7 に最も近い倍精度数
d = 3 * mpf(x)**2                    # f'(x) の厳密な値
print("  h       全体の誤差  打ち切り   丸めの寄与  丸めの評価")
for k in [4, 8, 12, 14, 16, 18, 20, 24, 28, 32, 36, 40]:
    h = 2.0**-k                      # x±h と 2h は浮動小数で正確
    D = (f(x + h) - f(x - h)) / (2 * h)          # 倍精度の中心差分
    Dm = (mpf(x + h)**3 - mpf(x - h)**3) / (2 * mpf(h))   # 厳密
    trunc = float(Dm - d)                        # 打ち切り（= h^2）
    rnd = abs(float(mpf(D) - Dm))                # 丸めによる差
    M0, M1 = (x + h)**3, 3 * (x + h)**2
    eps = 2 * u + u * u                          # 掛け算 2 回の誤差
    rb = (1 + u)**2 * eps * M0 / h + (2 * u + u * u) * M1
    print(f"2^-{k:<2d}   {float(abs(mpf(D) - d)):.2e}   {trunc:.2e}"
          f"   {rnd:.2e}    {rb:.2e}")
print(f"最適な刻みの予測 (3 eps M0/M3)^(1/3) = "
      f"{(3 * 2 * u * x**3 / 6) ** (1 / 3):.1e}")

  h       全体の誤差  打ち切り   丸めの寄与  丸めの評価
2^-4    3.91e-03   3.91e-03   5.06e-16    1.96e-15
2^-8    1.53e-05   1.53e-05   6.16e-15    2.02e-14
2^-12   5.96e-08   5.96e-08   2.75e-14    3.13e-13
2^-14   3.73e-09   3.73e-09   2.00e-13    1.25e-12
2^-16   2.34e-10   2.33e-10   1.16e-12    4.99e-12
2^-18   8.44e-12   1.46e-11   6.11e-12    2.00e-11
2^-20   1.16e-12   9.09e-13   2.55e-13    7.99e-11
2^-24   2.05e-10   3.55e-15   2.05e-10    1.28e-09
2^-28   1.19e-09   1.39e-17   1.19e-09    2.04e-08
2^-32   2.86e-08   5.42e-20   2.86e-08    3.27e-07
2^-36   1.22e-06   2.12e-22   1.22e-06    5.23e-06
2^-40   1.22e-06   8.27e-25   1.22e-06    8.37e-05
最適な刻みの予測 (3 eps M0/M3)^(1/3) = 3.4e-06

検算. (a) $x^3$ の中心差分は厳密には $3x^2 + h^2$ なので、打ち切りの列は $h^2$（$2^{-8}$ で $1.53 \times 10^{-5}$ など）に一致します。(b) 丸めの寄与は倍精度の値と mpmath の 50 桁の値の差で、すべての行で 命題 7.1 の丸めの評価（$(1+u)^2\varepsilon M_0/h + (2u+u^2)M_1$）以下です。(c) 全体の誤差は $h$ が大きいうちは $h^2$ に比例して減り（$h$ を $1/4$ にするごとに $1/16$）、$h \approx 2^{-18}$〜$2^{-20}$（$4 \times 10^{-6}$〜$10^{-6}$）で最小になって、その後は丸めの寄与が $1/h$ に比例して増えます。最小の位置は予測 $h^* \approx 3.4 \times 10^{-6}$ の近くです。$2^{-20}$ で丸めの寄与が評価より 2 桁以上小さいのは、丸めがたまたま打ち消し合ったためで、評価は最悪の場合の上界です。

8. 演習

演習 8.1易

$f \in C^5$ のとき中心差分 $D(h)$ は $D(h) = f'(x) + \frac{f'''(x)}{6}h^2 + O(h^4)$ を満たします（Taylor の定理を 4 次まで使えば確かめられます）。定理 4.2（$p = 2$）の外挿 $\bigl(4D(h/2) - D(h)\bigr)/3$ を $f$ の値で書き下し、誤差が $O(h^4)$ であることを確かめてください。

解答

$$\frac{4D(h/2) - D(h)}{3} = \frac13\Bigl(\frac{4\bigl(f(x+\frac h2) - f(x-\frac h2)\bigr)}{h} - \frac{f(x+h) - f(x-h)}{2h}\Bigr) = \frac{8\bigl(f(x+\frac h2) - f(x-\frac h2)\bigr) - \bigl(f(x+h) - f(x-h)\bigr)}{6h}$$

です。$D(h) = f'(x) + ch^2 + r(h)$、$c = f'''(x)/6$、$|r(h)| \le Ch^4$ なので、定理 4.2 から誤差は $\frac{1 + 2^{-2}}{3}Ch^4 = \frac{5}{12}Ch^4$ 以下です。刻みを $h/2$ に読み替えれば、よく使われる 5 点公式 $\frac{-f(x+2k) + 8f(x+k) - 8f(x-k) + f(x-2k)}{12k}$（$k = h/2$）です。

演習 8.2標準

$\int_0^1 x^4\,dx = 1/5$ について、$N_0 = 1$ の Romberg 表の $T_{0,0}, T_{1,0}, T_{2,0}, T_{1,1}, T_{2,1}, T_{2,2}$ を有理数で求め、$T_{2,2}$ が正確に $1/5$ になる理由を 定理 2.1 から説明してください。

解答

$T_{0,0} = \frac12(0 + 1) = \frac12$、$T_{1,0} = \frac12\bigl(0 + \frac1{16} + \frac12\bigr) = \frac{9}{32}$、$T_{2,0} = \frac14\bigl(\frac{1}{256} + \frac{16}{256} + \frac{81}{256} + \frac{128}{256}\bigr) = \frac{113}{512}$ です。$T_{1,1} = (4 \cdot \frac{9}{32} - \frac12)/3 = \frac{5}{24}$（Simpson 則）、$T_{2,1} = (4 \cdot \frac{113}{512} - \frac{9}{32})/3 = \frac{77}{384}$、$T_{2,2} = (16 \cdot \frac{77}{384} - \frac{5}{24})/15 = \frac{1152/384}{15} = \frac15$ です。$f = x^4$ では $f^{(5)} = 0$ なので、定理 2.1 を $p = 2$ で使うと $c_3 = 0$、剰余の積分も $f^{(6)} = 0$ で $0$ となり、$T(h) - I = c_1h^2 + c_2h^4$ が正確に成り立ちます。2 回の外挿で $h^2$ と $h^4$ の項が消えるので（系 4.3 で $r_k = 0$）、$T_{2,2} = I$ です。

演習 8.3標準

$f(\theta) = 1/(2 - \cos\theta)$ について $\int_0^{2\pi}f = 2\pi/\sqrt3$ です。(i) $|\mathrm{Im}\,z| \lt a$ で $|2 - \cos z| \ge 2 - \cosh a$ を示し、$a = 1$ として 定理 3.2 の 2 から $N = 20$ の台形則の誤差の上界を求めてください。(ii) $\hat f(n) = r^{|n|}/\sqrt3$（$r = 2 - \sqrt3$）を使って、実際の誤差を求めてください。

解答

(i) $z = x + iy$ で $|\cos z|^2 = \cos^2x + \sinh^2y \le \cosh^2y$ なので $|2 - \cos z| \ge 2 - |\cos z| \ge 2 - \cosh a$ です。$a = 1$ では $\cosh 1 \approx 1.543 \lt 2$ なので $f$ は帯状領域で正則で、$M = 1/(2 - \cosh 1) \approx 2.19$ です。上界は $4\pi M/(e^{20} - 1) \approx 5.7 \times 10^{-8}$ です。

(ii) 定理 3.2 の等式から、誤差は $2\pi\sum_{m \ne 0}\hat f(20m) = \frac{4\pi}{\sqrt3}\cdot\frac{r^{20}}{1 - r^{20}} \approx 2.6 \times 10^{-11}$ です。(i) の上界は 3 桁大きめです。$f$ の極は $\cos z = 2$、すなわち $z = \pm i\log(2 + \sqrt3) \approx \pm1.317i$ にあるので、$a$ をこの値に近づけるほど上界は真の誤差の減り方 $r^N = e^{-1.317N}$ に近づきます（ただし $M$ は大きくなります）。

演習 8.4難

前進差分 $D_+(h) = \bigl(f(x+h) - f(x)\bigr)/h$ について、命題 7.1 と同じ仮定（$f \in C^2$、$M_2 = \max|f''|$、関数値の相対誤差 $\varepsilon$、$x + h$ と $h$ は正確）のもとで誤差の上界を導き、最適な刻みが $u^{1/2}$ の程度、得られる精度も $u^{1/2}$ の程度であることを示してください。

解答

Taylor の定理（$n = 1$）から $D_+(h) - f'(x) = \frac h2f''(\xi)$ で、打ち切り誤差は $\frac{M_2}{2}h$ 以下です。丸めは中心差分と同じ計算で $(1+u)^2\cdot 2\varepsilon M_0/h + (2u + u^2)M_1$ 以下です（分母が $2h$ ではなく $h$ なので $\varepsilon$ の項が 2 倍）。$g(h) = \frac{M_2}{2}h + \frac{E'}{h}$（$E' = 2(1+u)^2\varepsilon M_0$）は $g'(h) = \frac{M_2}{2} - \frac{E'}{h^2}$ の符号から $h^* = \sqrt{2E'/M_2}$ で最小になり、最小値は $\sqrt{2M_2E'}$ です。$\varepsilon$ が $u$ の数倍なら $h^* \propto u^{1/2} \approx 10^{-8}$、誤差も $u^{1/2} \approx 10^{-8}$ の程度で、中心差分（$u^{2/3}$）より 3 桁ほど悪くなります。

参考文献

• K. E. Atkinson, An Introduction to Numerical Analysis, 2nd ed., Wiley, 1989 — Chapter 5（Euler–Maclaurin 公式、Richardson 外挿、Romberg 法）。
• L. N. Trefethen, J. A. C. Weideman, “The exponentially convergent trapezoidal rule”, SIAM Review 56 (2014), 385–458. doi:10.1137/130932132
• P. J. Davis, P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, 1984。
• R. Piessens, E. de Doncker-Kapenga, C. W. Überhuber, D. K. Kahaner, QUADPACK: A Subroutine Package for Automatic Integration, Springer, 1983。
• H. Takahasi, M. Mori, “Double exponential formulas for numerical integration”, Publ. RIMS Kyoto Univ. 9 (1974), 721–741.