数値積分 I：Newton–Cotes 則と Gauss 求積の誤差を Peano 核で証明する

0. この記事の要点

関数の値が何か所かでしか分からないとき、グラフの下の面積（積分）をどうやって見積もればよいでしょうか。点を結んだ台形の面積で近似するのが一番素朴な方法で、この記事ではその近似が実際にどれくらい外れるのか、点を増やせば真の値に近づくのか、少ない点で精度よく求めるにはどう点を選べばよいのか、を扱います。読み終えると、台形則や Simpson 則の誤差が刻み幅の何乗で小さくなるかを式で説明でき、Gauss 求積が少ない点数で高い精度を出せる理由を追えるようになります。

• 求積則 $Q(f) = \sum_i w_i f(x_i)$ の良さの第一の目安は精度次数（何次までの多項式を正確に積分するか）です。$n+1$ 点の則が $n$ 次まで正確であることと、補間多項式を積分した則（補間型）であることは同値です。
• Peano 核定理（誤差を、核と呼ぶ重み関数と導関数の積の積分で書く定理）は、精度次数 $k$ の則の誤差を $\int K(t) f^{(k+1)}(t)\,dt$ と表します。核 $K$ が符号を変えなければ誤差は $f^{(k+1)}(\xi)/(k+1)! \cdot E(x^{k+1})$ の形になり、台形則 $-\frac{(b-a)^3}{12}f''(\xi)$、Simpson 則 $-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$ が出ます。複合則では誤差が $O(h^2)$、$O(h^4)$ です。
• 求積則の列がすべての連続関数で収束するための必要十分条件は、多項式で収束し、重みの絶対値の和が有界であることです（Pólya–Steklov）。必要性の証明は一様有界性原理（関数解析の定理で、どの関数でも値が有界な連続線形汎関数の族は、ノルムも有界になるという主張）です。高次の Newton–Cotes 則は $n = 8$ で負の重みが現れ、重みの絶対値の和が増えていきます。
• Gauss 求積（直交多項式の零点を節点とする $n$ 点の則）は $2n-1$ 次まで正確で、誤差は $\frac{f^{(2n)}(\xi)}{(2n)!}\int P_n^2\,d\mu$ です。証明は Hermite 補間（節点で関数の値と 1 階導関数の値の両方を合わせる補間）です。重みが正なので、Pólya–Steklov の定理からすべての連続関数で収束します。

この記事で使う事前知識

• 補間多項式と Lagrange の基底多項式: 与えられた点で関数の値と一致する多項式と、それを作る部品（その点で 1、他の節点で 0 になる多項式）。→ 多項式補間
• Lebesgue 定数: 補間が関数値の誤差を何倍に増やすかを表す数。→ 多項式補間
• 定積分と置換積分・部分積分: 区間上の面積としての積分と、その計算規則。→ 積分の基本定理と定積分
• 一様ノルム: 連続関数のグラフが絶対値でどこまで大きくなるか（最大値）を測る量。→ 一様ノルム
• 凸関数: グラフの上の 2 点を結んだ弦が、つねにグラフより上にある関数。→ 凸関数

1. 動機

$\int_0^1 e^{-x^2}dx$ のように原始関数が初等関数で書けない積分、測定データや数値解として点の値しかわからない関数の積分は、有限個の値の重み付き和で近似するしかありません。研究で使われる数値計算でも、求積（scipy.integrate.quad、numpy.polynomial.legendre.leggauss など）は固有値計算に次いで多く現れる道具です。

求積則を選ぶとき、知りたいことは次の 3 つです。

1. 誤差の大きさ。刻み $h$ を半分にすると誤差は何分の 1 になるか。誤差の定数は何で決まるか。
2. 収束するか。点を増やしていけば、どんな連続関数でも真の値に近づくか。
3. 安定か。関数値の誤差（丸めや測定誤差）が何倍に増幅されるか。

1 には Peano 核定理が、2 と 3 には重みの絶対値の和 $\sum_i|w_i|$ が答えます。後者は多項式補間の章の Lebesgue 定数（補間が関数値の誤差を増幅する倍率）と同じ役割を果たします。重みがすべて正なら $\sum_i|w_i| = \sum_i w_i$ は区間の長さに等しく、安定性の心配はありません。Gauss 求積が「最も精度が高い」だけでなく「最も安心して使える」のはこのためです。

2. 求積則と精度次数

以下 $a \lt b$ とし、$I(f) = \int_a^b f(x)\,dx$ と書きます。$\mathcal{P}_n$ は $n$ 次以下の実係数多項式全体です。

定義 2.1（求積則）↗

$[a,b]$ の相異なる点 $x_0, \ldots, x_n$（節点）と実数 $w_0, \ldots, w_n$（重み）の組を求積則といい、$[a,b]$ 上の関数 $f$ について $Q(f) = \sum_{i=0}^{n} w_i f(x_i)$ と書く。$E(f) = I(f) - Q(f)$ を誤差という。

定義 2.2（精度次数）↗

すべての $p \in \mathcal{P}_k$ で $E(p) = 0$ となる最大の整数 $k \ge 0$ を求積則の精度次数という（$E(1) \ne 0$ のときは定義しない）。

$E$ は線形なので、精度次数が $k$ であることは $E(1) = E(x) = \cdots = E(x^k) = 0 \ne E(x^{k+1})$ と同値です。さらに $E(x^j)$（$j \le k$）がすべて $0$ なら、任意の $c$ について $E((x - c)^{k+1}) = E(x^{k+1})$ です（展開した低次の項の誤差は $0$）。

命題 2.3（補間型求積）✓

$[a,b]$ の相異なる $n+1$ 点 $x_0, \ldots, x_n$ を節点とする求積則について、次は同値である。

1. 精度次数が $n$ 以上である。
2. 重みが $w_i = \int_a^b \ell_i(x)\,dx$ である。ここで $\ell_i$ は Lagrange の基底多項式 $\ell_i(x) = \prod_{j \ne i}(x - x_j)/(x_i - x_j)$。

2 を満たす則を補間型という。このとき $Q(f)$ は $f$ の補間多項式 $p_n$ の積分 $I(p_n)$ に等しい。

証明（命題 2.3）

2 ⇒ 1: $p \in \mathcal{P}_n$ の補間多項式は $p$ 自身なので（定理 2.2［多項式補間］）$p = \sum_i p(x_i)\ell_i$ で、$I(p) = \sum_i p(x_i)I(\ell_i) = Q(p)$ です。一般の $f$ でも $Q(f) = \sum_i f(x_i)I(\ell_i) = I\bigl(\sum_i f(x_i)\ell_i\bigr) = I(p_n)$ です。

1 ⇒ 2: $\ell_i \in \mathcal{P}_n$ で $\ell_i(x_j) = \delta_{ij}$ なので、$I(\ell_i) = Q(\ell_i) = w_i$ です。

∎

等間隔の節点 $x_i = a + i(b-a)/n$ をもつ補間型の則を（閉じた）Newton–Cotes 則といいます。$n = 1$ が台形則 $\frac{b-a}{2}\bigl(f(a) + f(b)\bigr)$、$n = 2$ が Simpson 則 $\frac{b-a}{6}\bigl(f(a) + 4f(m) + f(b)\bigr)$（$m = (a+b)/2$）です。1 点 $m$ だけの補間型の則が中点則 $(b-a)f(m)$ です。中点則と Simpson 則は節点が中点について対称なので、奇数次の $(x - m)^{2j+1}$ も正確に積分し、精度次数はそれぞれ $1$、$3$ と、点の数から期待されるより 1 つ高くなります。

求積則は標準の区間 $[-1,1]$ で作って表にしておき、使うときに 1 次の変数変換で $[a,b]$ に移すのが普通です。そのとき誤差の定数がどう変わるかを確かめておきます。

命題 2.4（区間の変換と誤差）✓

$[-1,1]$ 上の求積則 $\hat Q(g) = \sum_i v_i g(t_i)$ について、整数 $k \ge 0$ と定数 $c$ があり、すべての $g \in C^{k+1}[-1,1]$ について $\int_{-1}^1 g - \hat Q(g) = c\,g^{(k+1)}(\eta)$ となる $\eta \in [-1,1]$ が存在するとする。$\varphi(t) = \frac{a+b}{2} + \frac{b-a}{2}t$ とし、$[a,b]$ 上の則を $x_i = \varphi(t_i)$、$w_i = \frac{b-a}{2}v_i$ で定める。このとき $f \in C^{k+1}[a,b]$ について、$\xi \in [a,b]$ があって

$$I(f) - Q(f) = c\Bigl(\frac{b-a}{2}\Bigr)^{k+2} f^{(k+1)}(\xi)$$

が成り立つ。

証明（命題 2.4）

$g = f \circ \varphi$ とおくと、$\varphi' = (b-a)/2$ は定数なので $g \in C^{k+1}[-1,1]$、$g^{(k+1)}(t) = \bigl(\frac{b-a}{2}\bigr)^{k+1}f^{(k+1)}(\varphi(t))$ です。置換積分（定理 6.1［積分の基本定理と定積分］）から $I(f) = \frac{b-a}{2}\int_{-1}^1 g$、定義から $Q(f) = \frac{b-a}{2}\hat Q(g)$ なので、$I(f) - Q(f) = \frac{b-a}{2}\,c\,g^{(k+1)}(\eta) = c\bigl(\frac{b-a}{2}\bigr)^{k+2}f^{(k+1)}(\varphi(\eta))$ です。$\xi = \varphi(\eta)$ とすればよいです。

∎

たとえば $[-1,1]$ の Simpson 則の誤差が $-\frac{1}{90}g^{(4)}(\eta)$ であることがわかれば（後の 定理 4.1 の $b - a = 2$ の場合）、$[a,b]$ では $-\frac{1}{90}\bigl(\frac{b-a}{2}\bigr)^5 f^{(4)}(\xi) = -\frac{(b-a)^5}{2880}f^{(4)}(\xi)$ です。誤差の次数 $(b-a)^{k+2}$ は、精度次数 $k$ から「区間の長さの $k+2$ 乗」として読み取れます。複合則の誤差 $O(h^{k+1})$（小区間 $N \propto 1/h$ 個の和）はこの帰結です。

3. Peano 核定理

誤差 $E(f)$ を $f$ の高階導関数で表すために、Taylor の定理の剰余を積分の形で使います。$s_+ = \max(s, 0)$、$s_+^0$ は $s \gt 0$ で $1$、$s \le 0$ で $0$ と書きます。

ねらいは次のとおりです。$f$ を Taylor の定理で「$k$ 次の多項式部分」と「積分の形の剰余」に分けます。精度次数が $k$ 以上なら、多項式部分の誤差は $0$ になり、誤差は剰余の分だけが残ります。剰余が積分で書けるので、誤差も $f^{(k+1)}$ に重み $K(t)$ をかけた積分で書けます。この $K$ が Peano 核で、式に現れる $(x_i - t)_+^k$ は、節点 $x_i$ が $t$ より右にあるときだけ値をもつ項（左にあれば $0$）です。

補題 3.1（Taylor の定理の積分剰余形）✓

$k \ge 0$、$f \in C^{k+1}[a,b]$ とする。$x \in [a,b]$ について

$$f(x) = \sum_{j=0}^{k}\frac{f^{(j)}(a)}{j!}(x - a)^j + \frac{1}{k!}\int_a^x (x - t)^k f^{(k+1)}(t)\,dt$$

が成り立つ。最後の項は $\frac{1}{k!}\int_a^b (x - t)_+^k f^{(k+1)}(t)\,dt$ とも書ける。

証明（補題 3.1）

$x = a$ では両辺が $f(a)$ なので、$x \gt a$ とします。$k$ についての帰納法です。$k = 0$ は微分積分学の基本定理（定理 5.4［積分の基本定理と定積分］）$f(x) - f(a) = \int_a^x f'(t)\,dt$ です。$k - 1$ で成り立つとし、$f \in C^{k+1}$ とします。$[a,x]$ 上の $C^1$ 級関数 $u(t) = f^{(k)}(t)$、$v(t) = -(x - t)^k/k!$（$v'(t) = (x-t)^{k-1}/(k-1)!$）に部分積分（定理 6.2［積分の基本定理と定積分］）を使うと

$$\frac{1}{(k-1)!}\int_a^x (x-t)^{k-1}f^{(k)}(t)\,dt = \Bigl[-\frac{(x-t)^k}{k!}f^{(k)}(t)\Bigr]_a^x + \frac1{k!}\int_a^x (x-t)^k f^{(k+1)}(t)\,dt$$

で、境界の項は $f^{(k)}(a)(x-a)^k/k!$ です。これを $k-1$ の式に代入すれば $k$ の式を得ます。最後の書き換えは、$t \gt x$ で $(x - t)_+^k = 0$、$t \le x$ で $(x-t)_+^k = (x-t)^k$（$k = 0$ は $t \lt x$ で $1$。1 点の違いは積分に影響しません）であることによります。

∎

定理 3.2（Peano 核定理）✓

$k \ge 1$ とし、$[a,b]$ 上の求積則 $Q$ の精度次数が $k$ 以上であるとする。$f \in C^{k+1}[a,b]$ ならば

$$E(f) = \int_a^b K(t)\, f^{(k+1)}(t)\,dt, \qquad K(t) = \frac{1}{k!}\Bigl(\frac{(b - t)^{k+1}}{k+1} - \sum_{i=0}^{n} w_i (x_i - t)_+^k\Bigr)$$

が成り立つ。$K$ は $[a,b]$ 上連続で、$f$ によらない。$K$ を $Q$ の（$k$ 次の）Peano 核という。

証明（定理 3.2）

補題 3.1 で $f = T + R$、$T$ は $k$ 次の Taylor 多項式、$R(x) = \frac1{k!}\int_a^b (x-t)_+^k f^{(k+1)}(t)\,dt$ と書きます。精度次数の仮定から $E(T) = 0$ なので $E(f) = E(R) = I(R) - Q(R)$ です。

$k \ge 1$ なので $(x, t) \mapsto (x - t)_+^k f^{(k+1)}(t)$ は正方形 $[a,b]^2$ 上で連続です。連続関数に対する Fubini の定理（系 4.3［重積分と累次積分］）から積分の順序を交換でき、

$$I(R) = \frac1{k!}\int_a^b f^{(k+1)}(t)\Bigl(\int_a^b (x - t)_+^k\,dx\Bigr)dt = \frac1{k!}\int_a^b f^{(k+1)}(t)\,\frac{(b-t)^{k+1}}{k+1}\,dt$$

です（内側の積分は $\int_t^b (x-t)^k dx$ で、微分積分学の基本定理から計算できます）。一方 $Q(R) = \frac1{k!}\int_a^b f^{(k+1)}(t)\sum_i w_i(x_i - t)_+^k\,dt$ は有限和と積分の交換（線形性）です。差をとれば主張を得ます。$K$ の連続性は $(x_i - t)_+^k$（$k \ge 1$）の連続性から従います。

∎

系 3.3（核が符号を変えない場合）✓

定理 3.2 の状況で、$K$ が $[a,b]$ 上で符号を変えない（つねに $K \ge 0$ またはつねに $K \le 0$）ならば、$\xi \in [a,b]$ があって

$$E(f) = \frac{f^{(k+1)}(\xi)}{(k+1)!}\,E(x^{k+1})$$

が成り立つ。

証明（系 3.3）

$K \ge 0$ としてよいです（$K \le 0$ なら $-K$ で同じ議論をします）。$g = f^{(k+1)}$ は $[a,b]$ 上連続なので最小値 $m$ と最大値 $M$ をとり（定理 2.3［平均値の定理とテイラーの定理］）、$mK \le Kg \le MK$ です。積分の単調性（命題 4.1［積分の基本定理と定積分］ の 2）から $m\int K \le \int Kg \le M\int K$ です。$\int K = 0$ なら $\int Kg = 0$ で、どの $\xi$ でも $E(f) = g(\xi)\int K$ が成り立ちます。$\int K \gt 0$ なら $\int Kg/\int K \in [m, M]$ なので、中間値の定理（定理 5.7［極限と連続性］）から $g(\xi) = \int Kg/\int K$ となる $\xi \in [a,b]$ があります。いずれの場合も $E(f) = g(\xi)\int_a^b K$ です。最後に $f = x^{k+1}$ に 定理 3.2 を使うと $f^{(k+1)} = (k+1)!$ なので $E(x^{k+1}) = (k+1)!\int_a^b K$ で、代入すれば主張を得ます。

∎

この系の価値は、誤差の定数 $E(x^{k+1})/(k+1)!$ を 1 回の計算で求められることです。核の符号さえ確かめれば、あとは $x^{k+1}$ を積分して求積則と比べるだけです。

次の図は、区間 $[-1,1]$ での台形則と Simpson 則の核です（Simpson 則の核の式は次節で求めます）。どちらも区間全体で $0$ 以下で、符号が変わりません。

台形則（左）と Simpson 則（右）の Peano 核。区間は [-1,1]。縦軸は各図で最小値が同じ深さになるようにとってあり、台形則の核の最小値は -1/2、Simpson 則の核の最小値は -1/72。

4. 台形則・Simpson 則と複合則

定理 4.1（中点則・台形則・Simpson 則の誤差）✓

$m = (a+b)/2$ とする。

1. $f \in C^2[a,b]$ なら、$\xi_1, \xi_2 \in [a,b]$ があって $$I(f) - (b-a)f(m) = \frac{(b-a)^3}{24}f''(\xi_1), \qquad I(f) - \frac{b-a}{2}\bigl(f(a) + f(b)\bigr) = -\frac{(b-a)^3}{12}f''(\xi_2).$$
2. $f \in C^4[a,b]$ なら、$\xi_3 \in [a,b]$ があって $$I(f) - \frac{b-a}{6}\bigl(f(a) + 4f(m) + f(b)\bigr) = -\frac{(b-a)^5}{2880}f^{(4)}(\xi_3).$$

証明（定理 4.1）

3 つの則の精度次数はそれぞれ $1, 1, 3$ 以上です（命題 2.3 と、中点について対称な則が $(x-m)^{2j+1}$ を正確に積分すること）。$L = b - a$ とおき、定理 3.2 の核を求めて 系 3.3 を使います。$t \in [a,b]$ では $(a - t)_+^k = 0$（$k \ge 1$）です。

中点則（$k = 1$）: $K(t) = \frac{(b-t)^2}{2} - L(m - t)_+$ です。$t \ge m$ では $K(t) = (b-t)^2/2 \ge 0$、$t \lt m$ では $s = t - a$ として $K(t) = \frac{(L - s)^2}{2} - L\bigl(\frac L2 - s\bigr) = \frac{s^2}{2} \ge 0$ です。$E(x^2) = E((x - m)^2) = \int_a^b (x-m)^2dx - 0 = L^3/12$ なので、誤差は $\frac{f''(\xi_1)}{2}\cdot\frac{L^3}{12}$ です。

台形則（$k = 1$）: $K(t) = \frac{(b-t)^2}{2} - \frac L2(b - t) = -\frac{(b-t)(t-a)}{2} \le 0$ です。$E((x-m)^2) = \frac{L^3}{12} - \frac L2\cdot 2\cdot\frac{L^2}{4} = -\frac{L^3}{6}$ なので、誤差は $\frac{f''(\xi_2)}{2}\cdot\bigl(-\frac{L^3}{6}\bigr)$ です。

Simpson 則（$k = 3$）: $h = L/2$ とします。$t \in [m, b]$ で $v = (t - m)/h \in [0,1]$ とおくと、$b - t = h(1-v)$、$(m - t)_+ = 0$ なので

$$K(t) = \frac16\Bigl(\frac{h^4(1-v)^4}{4} - \frac{h}{3}h^3(1-v)^3\Bigr) = \frac{h^4}{6}\,g(v), \qquad g(v) = \frac{(1-v)^4}{4} - \frac{(1-v)^3}{3}$$

です。$t \in [a, m]$ で $u = (m - t)/h \in [0,1]$ とおくと、$b - t = h(1+u)$、$m - t = hu$ なので $K(t) = \frac{h^4}{6}\bigl(\frac{(1+u)^4}{4} - \frac{4u^3 + (1+u)^3}{3}\bigr)$ です。両方を展開すると、どちらも $\frac{h^4}{6}\cdot\frac{3s^4 - 8s^3 + 6s^2 - 1}{12}$（$s = v$ または $u$）に等しく、$3s^4 - 8s^3 + 6s^2 - 1 = (s-1)^3(3s+1)$ なので

$$K(t) = \frac{h^4}{72}(s - 1)^3(3s + 1) \le 0 \qquad \Bigl(s = \frac{|t - m|}{h} \in [0,1]\Bigr)$$

です。$E(x^4) = E((x-m)^4) = \frac{2h^5}{5} - \frac{2h}{6}\cdot 2h^4 = -\frac{4h^5}{15}$ なので、誤差は $\frac{f^{(4)}(\xi_3)}{24}\cdot\bigl(-\frac{4h^5}{15}\bigr) = -\frac{h^5}{90}f^{(4)}(\xi_3) = -\frac{L^5}{2880}f^{(4)}(\xi_3)$ です。

∎

Simpson 則の証明で使った因数分解は、右辺を展開して確かめられます。$(s-1)^3 = s^3 - 3s^2 + 3s - 1$ に $3s + 1$ をかけると $3s^4 + s^3 - 9s^3 - 3s^2 + 9s^2 + 3s - 3s - 1 = 3s^4 - 8s^3 + 6s^2 - 1$ です。

例 4.2（中点則と台形則の誤差は符号が逆）

$f'' \gt 0$（凸関数）なら、中点則の誤差 $I - M = \frac{L^3}{24}f''(\xi_1) \gt 0$ で中点則は真の値を下回り、台形則の誤差 $I - T = -\frac{L^3}{12}f''(\xi_2) \lt 0$ で台形則は上回ります。図形的には、中点での接線の下の台形（面積は中点則と同じ）がグラフの下にあり、弦の下の台形がグラフの上にあるからです。$f''$ がほぼ一定なら誤差の比は $-1:2$ なので、$\frac{2M + T}{3}$ では主要項が打ち消し合います。実際 $\frac{2M + T}{3} = \frac L6(f(a) + 4f(m) + f(b))$ はちょうど Simpson 則です。

区間を $N$ 等分して、各小区間に同じ則を使うのが複合則です。$h = (b-a)/N$、$f_i = f(a + ih)$ とおくと、複合台形則と（$N$ が偶数のときの）複合 Simpson 則は

$$T_h(f) = h\Bigl(\frac{f_0}{2} + f_1 + \cdots + f_{N-1} + \frac{f_N}{2}\Bigr), \qquad S_h(f) = \frac h3\bigl(f_0 + 4f_1 + 2f_2 + 4f_3 + \cdots + 4f_{N-1} + f_N\bigr)$$

です。複合中点則は $M_h(f) = h\sum_{i=1}^{N} f\bigl(a + (i - \tfrac12)h\bigr)$ です。

系 4.3（複合則の誤差）✓

$h = (b-a)/N$ とする。

1. $f \in C^2[a,b]$ なら、$\xi, \xi' \in [a,b]$ があって $I(f) - T_h(f) = -\frac{(b-a)h^2}{12}f''(\xi)$、$I(f) - M_h(f) = \frac{(b-a)h^2}{24}f''(\xi')$。
2. $N$ が偶数で $f \in C^4[a,b]$ なら、$\xi'' \in [a,b]$ があって $I(f) - S_h(f) = -\frac{(b-a)h^4}{180}f^{(4)}(\xi'')$。

証明（系 4.3）

複合台形則の誤差は、各小区間 $[a + (i-1)h, a + ih]$ の台形則の誤差の和です（積分の区間加法性、命題 4.1［積分の基本定理と定積分］ の 3）。定理 4.1 から小区間ごとに $-\frac{h^3}{12}f''(\xi_i)$ と書けるので、和は $-\frac{h^3 N}{12}\cdot\frac1N\sum_i f''(\xi_i)$ です。平均 $\frac1N\sum_i f''(\xi_i)$ は $f''$ の $[a,b]$ 上の最小値と最大値（定理 2.3［平均値の定理とテイラーの定理］）の間にあるので、中間値の定理（定理 5.7［極限と連続性］）から $f''(\xi)$ と書けます。$Nh = b - a$ なので主張を得ます。中点則も同じです。Simpson 則は幅 $2h$ の小区間 $N/2$ 個の和で、小区間ごとの誤差は $-\frac{(2h)^5}{2880}f^{(4)}(\xi_i) = -\frac{h^5}{90}f^{(4)}(\xi_i)$ なので、同じ議論で $-\frac{N}{2}\cdot\frac{h^5}{90}f^{(4)}(\xi'') = -\frac{(b-a)h^4}{180}f^{(4)}(\xi'')$ です。

∎

例 4.4（刻みを半分にすると誤差は 1/4 と 1/16）

$\int_0^1 e^x\,dx = e - 1$ を複合台形則と複合 Simpson 則で計算し、刻みを半分にしたときの誤差の比を見ます。

# Python 3.12.10 / NumPy 2.5.2
import math
import numpy as np

def trap(f, a, b, N):                 # 複合台形則（N 等分）
    x = np.linspace(a, b, N + 1)
    y = f(x)
    return (b - a) / N * (y.sum() - (y[0] + y[-1]) / 2)

def simpson(f, a, b, N):              # 複合 Simpson 則（N は偶数）
    x = np.linspace(a, b, N + 1)
    y = f(x)
    return (b - a) / N / 3 * (y[0] + y[-1] + 4 * y[1:-1:2].sum()
                              + 2 * y[2:-1:2].sum())

I = math.e - 1                        # ∫_0^1 e^x dx
print("   N   台形の誤差  比     Simpson の誤差  比")
pt = ps = None
for N in [4, 8, 16, 32, 64]:
    et, es = trap(np.exp, 0, 1, N) - I, simpson(np.exp, 0, 1, N) - I
    rt = f"{pt / et:5.2f}" if pt else "  -  "
    rs = f"{ps / es:6.2f}" if ps else "   -  "
    print(f"{N:4d}   {et:.2e}  {rt}   {es:.2e}       {rs}")
    pt, ps = et, es
h = 1 / 64                            # 誤差の主要項の予測
print(f"予測 h^2(e-1)/12 = {h**2 * I / 12:.2e},"
      f" h^4(e-1)/180 = {h**4 * I / 180:.2e}")

   N   台形の誤差  比     Simpson の誤差  比
   4   8.94e-03    -     3.70e-05          -
   8   2.24e-03   4.00   2.33e-06        15.91
  16   5.59e-04   4.00   1.46e-07        15.98
  32   1.40e-04   4.00   9.10e-09        15.99
  64   3.50e-05   4.00   5.69e-10        16.00
予測 h^2(e-1)/12 = 3.50e-05, h^4(e-1)/180 = 5.69e-10

検算. (a) 真の値は閉じた形 $e - 1$ です。(b) 表の誤差は「数値 − 真の値」なので、系 4.3 の符号（$I - T_h \lt 0$）と合って正です。$f'' = f^{(4)} = e^x$ の平均値は $\int_0^1 e^x dx = e - 1$ なので、$N = 64$ の誤差は主要項の予測 $h^2(e-1)/12$、$h^4(e-1)/180$ と表示の桁まで一致します。(c) 比は $4 = 2^2$ と $16 = 2^4$ に近づき、誤差が $O(h^2)$、$O(h^4)$ であることを示します。

5. 求積列の収束：Pólya–Steklov の定理

点の数を増やした求積則の列 $Q_n(f) = \sum_{i=0}^{m_n} w_{n,i} f(x_{n,i})$（各 $n$ で節点は相異なる）が、すべての連続関数について真の値に収束するかを考えます。$C[a,b]$ には一様ノルム $\|f\| = \max_{[a,b]}|f|$ を入れます。$C[a,b]$ はこのノルムで完備です（補題 2.1［Arzelà–Ascoli と Stone–Weierstrass］）。

この節の言葉を先に整理します。線形汎関数とは、関数に実数を対応させる線形な写像のことで、$I$ や $Q$ がその例です。その作用素ノルムは、$\|f\| = 1$ の関数 $f$ について $|Q(f)|$ がどこまで大きくなるか（上限）です。ねらいは、収束するかどうかを重みの絶対値の和 $\sum_i|w_i|$ という 1 つの数だけで判定できるようにすることです。「多項式で収束すれば十分」は、連続関数を多項式でいくらでも近似できること（Weierstrass の多項式近似定理）から出ます。「重みの絶対値の和が有界」が必要になる理由は、一様有界性原理を使う証明の中で示します。

補題 5.1（点評価の線形結合のノルム）✓

$[a,b]$ の相異なる点 $x_0, \ldots, x_m$ と実数 $w_0, \ldots, w_m$ について、$C[a,b]$ 上の線形汎関数 $Q(f) = \sum_i w_i f(x_i)$ の作用素ノルムは $\sum_i |w_i|$ である。

証明（補題 5.1）

$|Q(f)| \le \sum_i|w_i||f(x_i)| \le \bigl(\sum_i|w_i|\bigr)\|f\|$ です。逆向きのために、$s_i$ を $w_i$ の符号（$w_i = 0$ なら $1$）とし、点を小さい順に並べて、$f$ を各 $x_i$ で値 $s_i$ をとる折れ線（両端の点の外では一定）とします。$f$ は連続で $\|f\| = 1$、$Q(f) = \sum_i w_is_i = \sum_i|w_i|$ です。

∎

定理 5.2（Pólya–Steklov の定理）✓

$L$ を $C[a,b]$ 上の線形汎関数で、ある定数 $C$ について $|L(f)| \le C\|f\|$ を満たすものとする（たとえば $L(f) = \int_a^b f(x)\,dx$）。求積則の列 $Q_n$ について、次は同値である。

1. すべての $f \in C[a,b]$ について $Q_n(f) \to L(f)$。
2. (i) すべての多項式 $p$ について $Q_n(p) \to L(p)$、かつ (ii) $\sup_n \sum_i |w_{n,i}| \lt \infty$。

証明（定理 5.2）

2 ⇒ 1. $W = \sup_n\sum_i|w_{n,i}|$ とします。$f \in C[a,b]$ と $\varepsilon \gt 0$ について、Weierstrass の多項式近似定理（系 4.4［Arzelà–Ascoli と Stone–Weierstrass］）から $\|f - p\| \lt \varepsilon$ となる多項式 $p$ があります。

$$|L(f) - Q_n(f)| \le |L(f - p)| + |L(p) - Q_n(p)| + |Q_n(p - f)| \le C\varepsilon + |L(p) - Q_n(p)| + W\varepsilon$$

で、(i) から $\limsup_n |L(f) - Q_n(f)| \le (C + W)\varepsilon$ です。$\varepsilon$ は任意なので $Q_n(f) \to L(f)$ です。

1 ⇒ 2. (i) は 1 の特別な場合です。(ii) を示します。$Q_n$ は Banach 空間 $C[a,b]$ 上の有界線形汎関数で、補題 5.1 からノルムは $\sum_i|w_{n,i}|$ です。各 $f$ について $Q_n(f)$ は収束列なので有界、すなわち $\sup_n|Q_n(f)| \lt \infty$ です。一様有界性原理（定理 5.1［Banach 空間の基本定理］。$Y = \mathbb{R}$）から $\sup_n\|Q_n\| \lt \infty$ で、これが (ii) です。

∎

$\sum_i |w_{n,i}|$ は求積則の安定性の指標でもあります。関数値に大きさ $\varepsilon$ 以下の誤差 $\delta_i$ が乗ると、$Q_n$ の値は $|\sum_i w_{n,i}\delta_i| \le \varepsilon\sum_i|w_{n,i}|$ だけずれ、補題 5.1 からこの上界は達成されます。補間型の則で $Q_n(1) = b - a$ なので、重みがすべて正なら $\sum_i|w_{n,i}| = b - a$ で、(ii) は自動的に成り立ちます。複合則（重みは正）はこの場合です。問題は重みが負になるときです。

例 5.3（Newton–Cotes 則の負の重み）

閉じた Newton–Cotes 則の重みを、fractions による正確な有理数計算で求めます。区間は $[0,1]$ に正規化し（重みの和は $1$）、点 $0, 1, \ldots, n$ の Lagrange 基底を係数のまま積分します。

# Python 3.12.10（標準ライブラリのみ）
from fractions import Fraction as F

def nc_weights(n):
    """閉じた Newton–Cotes 則（点 0,1,…,n）の重みを [0,1] 用に正規化"""
    ws = []
    for i in range(n + 1):
        c = [F(1)]                      # ℓ_i の係数（低次から）
        for k in range(n + 1):
            if k == i:
                continue
            d = F(i - k)
            c = [F(0)] + c              # t を掛ける
            for j in range(len(c) - 1):
                c[j] -= k * c[j + 1]    # -k を掛けて足す
            c = [x / d for x in c]
        ws.append(sum(cj * F(n) ** (j + 1) / (j + 1)
                      for j, cj in enumerate(c)) / n)
    return ws

print(" n   最小の重み  重みの和  |重み|の和")
for n in [1, 2, 4, 6, 7, 8, 9, 10, 12, 16, 20]:
    w = nc_weights(n)
    print(f"{n:2d}   {float(min(w)):9.5f}   {str(sum(w)):6s}"
          f"  {float(sum(abs(x) for x in w)):.3f}")
print("n = 2:", [str(x) for x in nc_weights(2)])
print("n = 8:", [str(x) for x in nc_weights(8)[:5]], "…")

 n   最小の重み  重みの和  |重み|の和
 1     0.50000   1       1.000
 2     0.16667   1       1.000
 4     0.07778   1       1.000
 6     0.03214   1       1.000
 7     0.04346   1       1.000
 8    -0.16014   1       1.451
 9     0.01205   1       1.000
10    -0.43516   1       3.065
12    -1.39221   1       7.532
16   -10.47598   1       58.457
20   -90.00537   1       544.177
n = 2: ['1/6', '2/3', '1/6']
n = 8: ['989/28350', '2944/14175', '-464/14175', '5248/14175', '-454/2835'] …

検算. (a) $n = 2$ の重みは Simpson 則の $1/6, 2/3, 1/6$ に一致し、重みの和はどの $n$ でも正確に $1$（定数の積分が正確）です。$n = 8$ の重み（9 点の則。対称なので前半だけ表示）は分母をそろえると $989, 5888, -928, 10496, -4540$（$/28350$）で、数表（Abramowitz–Stegun『Handbook of Mathematical Functions』§25.4 の閉じた Newton–Cotes 公式）と一致します。(b) 計算はすべて有理数で行っており、丸め誤差はありません。(c) 負の重みは $n = 8$ で初めて現れます（同じ関数で $n = 1, \ldots, 24$ を調べると、負の重みがないのは $n \le 7$ と $n = 9$ だけでした）。重みの絶対値の和は、表の $n = 12, 16, 20$ で $n$ を 4 増やすごとに 7〜10 倍に増えています。

注意 5.4（Newton–Cotes 則の重みの絶対値の和の発散）

閉じた Newton–Cotes 則の重みの絶対値の和は $n \to \infty$ で発散し、その増え方は指数的であることが R. O. Kuzmin（1931）により示されています。定理 5.2 により、このとき Newton–Cotes 則の列が収束しない連続関数が存在します。漸近評価の証明は長いので、P. J. Davis, P. Rabinowitz『Methods of Numerical Integration』（2nd ed., Academic Press, 1984）、H. Brass, K. Petras『Quadrature Theory』（AMS, 2011）を参照してください。本章では負の重みが現れることを上の例で正確に計算するにとどめ、発散は定理として使いません。実用で高い精度が要るときは、Newton–Cotes 則の次数を上げるのではなく、低次の複合則に外挿を組み合わせるか（Euler–Maclaurin 展開・Richardson 外挿・Romberg 法）、Gauss 求積を使います。

6. Gauss 求積

節点も自由に選べるなら、$n$ 点で $2n$ 個の自由度（節点と重み）があり、$2n - 1$ 次までの多項式を正確に積分できると期待できます。これを実現するのが Gauss 求積で、その基本的な性質は直交多項式の章で証明されています。この節では重み付きの積分も扱うため、$\mathbb{R}$ 上の正の測度 $\mu$（重みつき積分の「重み」を一般にしたもの）で直交多項式の章の条件（定義 2.1［直交多項式］。すべてのモーメントが有限で、有限個の点に集中しない）を満たし、さらに $\mu(\mathbb{R} \setminus [a,b]) = 0$ となるものを考えます。$d\mu = dx$（$[-1,1]$ 上、Gauss–Legendre）、$d\mu = dx/\sqrt{1-x^2}$（Gauss–Chebyshev）などです。

$P_n$ を $\mu$ のモニック（最高次の係数が $1$）な直交多項式（定理 2.2［直交多項式］）とし、$h_n = \int P_n^2\,d\mu \gt 0$ と書きます。$P_n$ の零点 $x_1 \lt \cdots \lt x_n$ はすべて実数で相異なり、$(a,b)$ に入ります（定理 5.1［直交多項式］）。これを節点とし、重みを $w_i = \int \ell_i\,d\mu$ とした則

$$G_n(f) = \sum_{i=1}^{n} w_i f(x_i)$$

を $n$ 点 Gauss 求積といいます。定理 6.1［直交多項式］ により、$G_n$ は $\mathcal{P}_{2n-1}$ の元を正確に積分し、重みはすべて正で、$n$ 点のどんな則も $\mathcal{P}_{2n}$ のすべてを正確に積分することはできません。節点と重みは Jacobi 行列（三項漸化式の係数から作る三重対角行列）の固有値問題として計算できます（Golub–Welsch、定理 6.3［直交多項式］）。以下で足すのは、誤差の公式と収束です。

誤差公式の方針は次のとおりです。節点で値と 1 階導関数を $f$ に合わせた $2n-1$ 次以下の多項式 $H$ をとると、$G_n$ は $H$ を正確に積分します。すると誤差は $f - H$ の積分になり、Hermite 補間の誤差公式から $f - H$ は $P_n^2$ に $f^{(2n)}/(2n)!$ のある値をかけたものです。$P_n^2$ は 0 以上なので、積分は $f^{(2n)}/(2n)!$ の最小値と最大値の間の定数に $\int P_n^2\,d\mu$ をかけたものとして書けます。

定理 6.1（Gauss 求積の誤差公式）✓

$n \ge 1$ とし、$f$ は $[a,b]$ を含む開区間で $2n$ 回微分可能で、$f^{(2n)}$ は $[a,b]$ 上連続とする。このとき $\xi \in [a,b]$ があって

$$\int f\,d\mu - G_n(f) = \frac{f^{(2n)}(\xi)}{(2n)!}\int P_n^2\,d\mu$$

が成り立つ。

証明（定理 6.1）

$H \in \mathcal{P}_{2n-1}$ を、節点 $x_1, \ldots, x_n$ で $f$ と値・1 階導関数が一致する Hermite 補間多項式とします（定理 4.2［多項式補間］。$n$ 点なので次数は $2n-1$ 以下）。$\omega(x) = \prod_i(x - x_i)$ と $P_n$ はともにモニックな $n$ 次式で同じ $n$ 個の零点をもつので、差は $\mathcal{P}_{n-1}$ の元で $n$ 個の零点をもち、系 2.3［多項式補間］ から $\omega = P_n$ です。

$G_n$ は $H$ を正確に積分し、$H(x_i) = f(x_i)$ なので、$G_n(f) = G_n(H) = \int H\,d\mu$ です。よって誤差は $\int(f - H)\,d\mu$ です。$m$、$M$ を $f^{(2n)}/(2n)!$ の $[a,b]$ 上の最小値・最大値とします（定理 2.3［平均値の定理とテイラーの定理］）。$x \in [a,b]$ について、Hermite 補間の誤差公式（定理 4.2［多項式補間］ の 2）の $\xi$ は $[a,b]$ に入るので

$$m\,P_n(x)^2 \le f(x) - H(x) \le M\,P_n(x)^2 \qquad (x \in [a,b])$$

です。$\mu$ は $[a,b]$ の外に重さをもたず、$[a,b]$ 上の連続関数は有界で $\mu$ は有限なので、これらの関数は $\mu$ について可積分です。非負の関数 $(f - H) - mP_n^2$ と $MP_n^2 - (f - H)$ の積分は定義から $0$ 以上（定義 5.1［ルベーグ積分の定義と収束定理］）で、積分の線型性（命題 7.2［ルベーグ積分の定義と収束定理］）と合わせて

$$m\,h_n \le \int(f - H)\,d\mu \le M\,h_n$$

です。$h_n \gt 0$ なので $\int(f-H)\,d\mu/h_n \in [m, M]$ で、中間値の定理（定理 5.7［極限と連続性］）から $f^{(2n)}(\xi)/(2n)!$ がこの値に等しい $\xi \in [a,b]$ があります。

∎

$[-1,1]$ 上の Gauss–Legendre（$d\mu = dx$）では、モニックな Legendre 多項式について $h_n = \frac{2^{2n+1}(n!)^4}{(2n+1)((2n)!)^2}$ で、$n = 2$ なら $h_2 = 8/45$、誤差は $f^{(4)}(\xi)/135$ です。同じ 4 次の誤差をもつ Simpson 則（3 点、誤差 $-f^{(4)}(\xi)/90$）より 1 点少なく、定数も小さくなっています。$f$ が解析関数なら、導関数を使わない幾何級数的な誤差の上界（Bernstein 楕円による評価）が 定理 4.3［誤差保証つきの求積と関数近似］ にあります。

系 6.2（Gauss 求積の収束）✓

すべての $f \in C[a,b]$ について、$n \to \infty$ で $G_n(f) \to \int f\,d\mu$ である。

証明（系 6.2）

$L(f) = \int f\,d\mu$ は $C[a,b]$ 上の線形汎関数で、$|L(f)| \le \int|f|\,d\mu \le \mu(\mathbb{R})\|f\|$ です（命題 7.2［ルベーグ積分の定義と収束定理］ の 2 と、定理 6.1 の証明と同じ単調性）。定理 5.2 の 2 を確かめます。(i) 多項式 $p$ について、$2n - 1 \ge \deg p$ となる $n$ では $G_n(p) = L(p)$ です（定理 6.1［直交多項式］ の 1）。(ii) 重みは正なので（同じ定理の 2）、$\sum_i|w_i| = \sum_i w_i = G_n(1) = \mu(\mathbb{R})$ で、$n$ によりません。よって $G_n(f) \to L(f)$ です。

∎

この収束はどんな連続関数でも成り立ちますが、速さは $f$ の滑らかさで決まります。$f \in C^{2n}$ なら 定理 6.1 から誤差は $f^{(2n)}$ の大きさと $h_n/(2n)!$ で決まり、解析関数なら幾何級数的に減ります。$f$ が滑らかでないと、誤差は $n$ の多項式の速さでしか減りません。

例 6.3（Gauss–Chebyshev 求積は角度の中点則）

$d\mu = dx/\sqrt{1 - x^2}$ では、例 6.4［直交多項式］ のとおり、節点は $x_k = \cos\frac{(2k-1)\pi}{2n}$、重みはすべて $\pi/n$ です。$x = \cos\theta$ と置換すると $\int_{-1}^1 \frac{f(x)}{\sqrt{1-x^2}}dx = \int_0^\pi f(\cos\theta)\,d\theta$ で、$\theta_k = \frac{(2k-1)\pi}{2n}$ は $[0, \pi]$ を $n$ 等分した小区間の中点です。つまり Gauss–Chebyshev 求積は、$\theta$ の変数で見ると等間隔の複合中点則にほかなりません。$\theta \mapsto f(\cos\theta)$ は周期 $2\pi$ の偶関数なので、外挿と Romberg 法の章で示す周期関数の台形則（中点則も同じ）の速い収束がここで効きます。

例 6.4（Golub–Welsch と、滑らかな関数・滑らかでない関数での収束）

Legendre 多項式（モニック）の三項漸化式の係数 $\beta_k = k^2/(4k^2 - 1)$ から Jacobi 行列を作り、numpy.linalg.eigh の固有値・固有ベクトルで節点と重みを求めて（定理 6.3［直交多項式］）、NumPy の leggauss と照合します。そのうえで $\int_{-1}^1 e^x dx = e - 1/e$ と $\int_{-1}^1 |x|\,dx = 1$ の誤差を比べます。

# Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0
from math import factorial as fac
import numpy as np
from numpy.polynomial.legendre import leggauss
from mpmath import mp, quad, exp as mexp

def golub_welsch(n):                  # Legendre の Jacobi 行列から
    k = np.arange(1, n)
    b = k / np.sqrt(4.0 * k * k - 1)  # 副対角 √β_k
    lam, U = np.linalg.eigh(np.diag(b, 1) + np.diag(b, -1))
    return lam, 2 * U[0] ** 2         # 節点と重み（μ(R) = 2）

mp.dps = 30
ref_exp = float(quad(mexp, [-1, 1]))           # e - 1/e
ref_abs = float(quad(abs, [-1, 0, 1]))         # 1
print("  n  一致   e^x の誤差  比     |x| の誤差  n^2×誤差")
for n in [2, 3, 4, 6, 8, 16, 32, 64]:
    x, w = golub_welsch(n)
    xl, wl = leggauss(n)
    ok = max(np.abs(x - xl).max(), np.abs(w - wl).max()) < 1e-13
    ee = abs(w @ np.exp(x) - ref_exp)
    hn = (2**n * fac(n)**2 / fac(2 * n))**2 * 2 / (2 * n + 1)
    r = f"{ee / (hn / fac(2 * n)):.3f}" if ee > 1e-12 else "  -  "
    s = f"{ee:.2e}" if ee > 1e-12 else "< 1e-12 "
    ea = abs(w @ np.abs(x) - ref_abs)
    print(f"{n:3d}  {ok!s:5}  {s}    {r}  {ea:.2e}    {n*n*ea:.3f}")
print(f"参照値と閉じた形の差: {abs(ref_exp - (np.e - 1/np.e)):.1e},"
      f" {abs(ref_abs - 1):.1e}")

  n  一致   e^x の誤差  比     |x| の誤差  n^2×誤差
  2  True   7.71e-03    1.040  1.55e-01    0.619
  3  True   6.55e-05    1.031  1.39e-01    1.254
  4  True   2.95e-07    1.025  4.25e-02    0.681
  6  True   1.57e-12    1.017  1.99e-02    0.716
  8  True   < 1e-12       -    1.15e-02    0.738
 16  True   < 1e-12       -    3.03e-03    0.776
 32  True   < 1e-12       -    7.79e-04    0.798
 64  True   < 1e-12       -    1.98e-04    0.810
参照値と閉じた形の差: 0.0e+00, 0.0e+00

検算. (a) 参照値は mpmath の quad（30 桁）で計算し、閉じた形 $e - 1/e$ と $1$ に倍精度で一致します。(b) Golub–Welsch の節点・重みは、別の実装である NumPy の leggauss（Legendre 多項式の随伴行列の固有値を Newton 法で 1 回改良し、重みは別の公式で求める）の値と $10^{-13}$ 以内で一致します。(c) 「比」の列は、誤差を 定理 6.1 の $h_n/(2n)!$ で割った値で、定理からこれは $e^{\xi}$（$\xi \in [-1,1]$）、すなわち $[e^{-1}, e]$ に入るはずです。実際 $1.02$〜$1.04$ で、$\xi$ が $0$ の近くにあることを示しています。$n = 8$ 以降は誤差が丸めの水準なので表示していません。一方 $|x|$ は $x = 0$ で微分できず、誤差は $n^2 \times$ 誤差がほぼ一定、すなわち $O(n^{-2})$ の速さでしか減りません（$n$ が奇数のときは $0$ が節点になり、定数が変わります）。系 6.2 の収束は保証されていますが、速さは $f$ の滑らかさで大きく変わります。

本章の求積則を、1 区間あたりの点の数と誤差の形でまとめます（$L = b - a$、$\xi \in [a,b]$）。

則	点の数	精度次数	1 区間の誤差 $I - Q$	重み
中点則	1	1	$\frac{L^3}{24}f''(\xi)$	正
台形則	2	1	$-\frac{L^3}{12}f''(\xi)$	正
Simpson 則	3	3	$-\frac{L^5}{2880}f^{(4)}(\xi)$	正
2 点 Gauss–Legendre	2	3	$\frac{L^5}{4320}f^{(4)}(\xi)$	正
$n$ 点 Gauss–Legendre	$n$	$2n-1$	$\frac{f^{(2n)}(\xi)}{(2n)!}\bigl(\frac L2\bigr)^{2n+1}h_n$	正
9 点 Newton–Cotes	9	9	（本章では扱わない）	負を含む

2 点 Gauss–Legendre の行は、$[-1,1]$ での誤差 $f^{(4)}(\xi)/135$（演習 7.3 で確かめます）を 命題 2.4 で移したものです。$n$ 点の行も同じ変換で、$h_n$ は $[-1,1]$ 上のモニックな Legendre 多項式の $\int P_n^2$ です。実務での使い分けについての私の考えは次のとおりです。関数が滑らかで任意の点で評価できるなら Gauss–Legendre を第一の候補にします。等間隔の点の値しかないなら複合 Simpson 則か、次の章の外挿を組み合わせた台形則を使います。被積分関数に特異性や急な変化があるときは、区間を分けるか、特異性を重みに取り込んだ Gauss 求積（Gauss–Chebyshev など）を使います。どの場合も、点を増やしたときの値の変化から誤差を見積もる習慣が大切です。

7. 演習

演習 7.1易

$[-1,1]$ 上の中点則 $2f(0)$ と Simpson 則 $\frac13(f(-1) + 4f(0) + f(1))$ の精度次数を、定義に従って $E(x^j)$ を計算して求めてください。

解答

中点則: $E(1) = 2 - 2 = 0$、$E(x) = 0 - 0 = 0$、$E(x^2) = 2/3 - 0 \ne 0$ なので精度次数は $1$ です。Simpson 則: $E(1) = 2 - 2 = 0$、$E(x) = 0$、$E(x^2) = 2/3 - 2/3 = 0$、$E(x^3) = 0 - 0 = 0$、$E(x^4) = 2/5 - 2/3 = -4/15 \ne 0$ なので精度次数は $3$ です。奇数次の項は対称性で消えます。

演習 7.2標準

$\int_0^1 e^x dx$ を複合台形則と複合 Simpson 則で計算するとき、系 4.3 の評価で誤差 $10^{-8}$ 以下が保証される最小の分割数 $N$ をそれぞれ求めてください（Simpson 則の $N$ は偶数）。

解答

$[0,1]$ 上で $|f''| = |f^{(4)}| \le e$ です。台形則は $h^2 e/12 \le 10^{-8}$、すなわち $h \le \sqrt{12 \times 10^{-8}/e} \approx 2.101 \times 10^{-4}$ で、$N \ge 4759.4$ から $N = 4760$ です。Simpson 則は $h^4 e/180 \le 10^{-8}$、すなわち $h \le (180 \times 10^{-8}/e)^{1/4} \approx 0.02853$ で、$N \ge 35.06$、偶数なので $N = 36$ です。Simpson 則は約 130 分の 1 の関数値の計算で同じ保証を得ます。

演習 7.3標準

$[-1,1]$ 上の 2 点 Gauss–Legendre 則の節点と重みを $P_2(x) = x^2 - 1/3$ から求め、定理 6.1 から誤差が $f^{(4)}(\xi)/135$ であることを示してください。

解答

節点は $P_2$ の零点 $\pm1/\sqrt3$ です。重みは $w_1 = \int_{-1}^1 \frac{x - 1/\sqrt3}{-2/\sqrt3}dx = 1$、同様に $w_2 = 1$ です。$h_2 = \int_{-1}^1(x^2 - 1/3)^2dx = \frac25 - \frac49 + \frac29 = \frac{8}{45}$ なので、誤差は $\frac{f^{(4)}(\xi)}{4!}\cdot\frac{8}{45} = \frac{f^{(4)}(\xi)}{135}$ です。確認として $f = x^4$ では $\int x^4 = 2/5$、$G_2(x^4) = 2/9$ で差は $8/45 = 24/135$ と一致します。

演習 7.4難

複合台形則 $T_h$（$h = (b-a)/N$）が、$N \to \infty$ ですべての $f \in C[a,b]$ について $I(f)$ に収束することを、定理 5.2 を使って示してください。また、$f$ がどの点でも微分できない連続関数でもこの結論が成り立つことと、系 4.3 が使えないことの違いを説明してください。

解答

重み $h/2, h, \ldots, h, h/2$ は正で、和は $Nh = b - a$ なので、定理 5.2 の (ii) は $W = b - a$ で成り立ちます。(i) は、多項式 $p$ が $C^2$ 級なので 系 4.3 から $|I(p) - T_h(p)| \le \frac{(b-a)h^2}{12}\max_{[a,b]}|p''| \to 0$ です。よって 1 が成り立ちます。系 4.3 は $f \in C^2$ を仮定し、収束の速さ（$O(h^2)$）まで与えます。Pólya–Steklov の定理は滑らかさを仮定しない代わりに、収束の速さについては何も言いません。どこでも微分できない連続関数（Weierstrass 型の関数）では、$O(h^2)$ の速さは一般には期待できず、収束することだけが保証されます。

参考文献

• K. E. Atkinson, An Introduction to Numerical Analysis, 2nd ed., Wiley, 1989 — Chapter 5（数値積分、Peano 核、Gauss 求積）。
• P. J. Davis, P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, 1984。
• H. Brass, K. Petras, Quadrature Theory: The Theory of Numerical Integration on a Compact Interval, Mathematical Surveys and Monographs 178, AMS, 2011。
• G. H. Golub, J. H. Welsch, “Calculation of Gauss quadrature rules”, Mathematics of Computation 23 (1969), 221–230.
• M. Abramowitz, I. A. Stegun (eds.), Handbook of Mathematical Functions, National Bureau of Standards, 1964 — §25.4（数値積分の公式）。
• 杉原正顯・室田一雄『数値計算法の数理』岩波書店、1994 — 数値積分の章。