第5章　数値積分

5.1 数値積分の公式

5.1.1 数値積分の公式

関数 $f(x)$ の定積分

\[I = \int_a^b f(x) \ dx \tag{5.1}\]

の値は、 $f(x)$ が簡単な場合を除いて、解析的に求めることは困難である。 定積分の近似値を数値的に求める公式を導くことが、この章の目的である。

数値積分の公式は、 $f(x)$ を容易に積分できる関数で近似し、 $(5.1)$ の $f(x)$ の代わりに積分し、

\[I \cong I_n \equiv h \sum_{i=0}^n w_i f(x_i) \tag{5.2}\]

の形に求める。 ここに、 $x_i$ は 分点 (point)、 $w_i$ は 重み (weight) という。 $n + 1$ 個の分点は $x_0 < x_1 < x_2 < x_3 < \cdots < x_{n-1} < x_n$ のように、相異なるものとする。 また、容易に積分できる関数としては、 $x$ の多項式を選ぶのが普通である。

端点が積分の上限と下限に一致するとき、すなわち $x_0 = a$ 、 $x_n = b$ のとき、 $I_n$ は 閉型公式 という。 一致しないとき、すなわち $a < x_0$ 、 $x_n < b$ のとき、 開型公式 という。

重み $w_i$ は $x_i$ にはよるが、 $f(x)$ にはよらない定数とする。 分点と重みは、誤差

\[E_n = I - I_n \tag{5.3}\]

が出来るだけ小さくなるように選ぶ。 誤差 $E_n$ は打切り誤差を以て表す。 $(5.2)$ は $n + 1$ 分点の公式である。 通常、 $n$ が大きいほど $E_n$ の小さな公式が導き出される。

なお、 数値微分 については、問題 5-12 に簡単にふれておくに止める。

5.1.2 台形公式

最も簡単な数値積分の閉型公式は 台形公式 (trapezoidal formula) である。 まず、全閉区間 $[\ a,\ b\ ]$ を $N$ 等分し、分点を等間隔に置く。 すなわち

\[x_i = a + ih \qquad (0 \leqq i \leqq N ) \qquad \text{ただし} \quad h = \frac{b - a}{N} \tag{5.4}\]

もちろん、 $x_0 = a$ 、 $x_N = b$ である。 $f(x)$ の近似関数は1次式とする。 2つの分点 $x_i$ と $x_{i+1}$ で $f(x)$ と一致する1次式は

\[p_1 (x) = \frac{x - x_{i+1}}{x_i - x_{i+1}} f(x_i) + \frac{x - x_i}{x_{i+1} - x_i} f(x_{i+1}) \tag{5.5}\]

であることは容易に確かめられる。 $p_1(x)$ を一つの区間 $x_i$ から $x_{i+1}$ まで積分すると、一つの台形の面積が次のように得られる：

\[\begin{aligned} \int_{x_i}^{x_{i+1}} p_1(x) \ dx & = \int_0^h p_1(y + x_i) \ dy \cr & = \int_0^h \left[ \frac{y - h}{-h} f(x_i) + \frac{y}{h} f(x_{i+1}) \right] \ dy \cr & = \frac{h}{2} [f(x_i) + f(x_{i+1})] \end{aligned}\]

これより、

\[\boxed{ \begin{aligned} \quad \cr & \text{台形公式} \cr & \quad I_1 = \frac{h}{2} [ f(x_i) + f(x_{i+1} ) ] \cr & \qquad \text{ただし} \quad h = x_{i+1} - x_i \cr && \quad \cr \end{aligned} } \tag{5.6}\]

台形公式 $(5.6)$ の場合、分点は $x_i$ と $x_{i+1}$ 、重み $w_i$ は

\[w_0 = w_1 = \frac{1}{2} \tag{5.7}\]

ととったことになる。

この面積を $x_0 = a$ から $x_N = b$ の $N$ 区間について加えると、

\[I = \int_a^b f(x) \ dx = \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} f(x) \ dx\]

であるから、 $(5.6)$ より、

\[I \cong \frac{h}{2} \sum_{i=0}^{N-1} [\ f(x_i) + f(x_{i+1})\ ]\]

隣合った各区間の上限と下限は一致するから、

\[\boxed{ \begin{aligned} \quad \cr & \text{台形公式 (複合公式)} \cr & \quad I_1 = h \left[ \frac{f(a) + f(b)}{2} + \sum_{i=1}^{N-1} f(x_i) \right] & \text{ただし} \quad h = \frac{b - a}{N} \quad \end{aligned} } \tag{5.8}\]

が得られる。 この公式が台形の面積の和であるという図上の意味は、図 5.1 より明かである。

5.1.3 シンプソンの公式

台形公式に次いで簡単な閉型公式は、シンプソン (Simpson) の公式である。 先ず、積分範囲を $N$ 個の中区間に分割し、各中区間を2等分し、小区間の幅を

\[h = \frac{b - a}{2N} \tag{5.9}\]

と置く (図 5.2 参照)。 分点は

\[x_i = a + ih \qquad (0 \leqq i \leqq 2N ) \tag{5.10}\]

とする。 各中区間の被積分関数 $f(x)$ を2次式 $p_2(x)$ で近似する。 2次式を完全に (唯一に) 定めるには3分点が必要である。 3分点 $x = x_{2i},\ x_{2i+1},\ x_{2i+2}$ で $p_2(x) = f(x)$ である 2 次式は

\[\begin{aligned} p_2(x) = &\ \frac{(x - x_{2i+1})(x - x_{2i+2})}{(x_{2i} - x_{2i+1})(x_{2i} - x_{2i+2})} f(x_{2i}) \cr & \quad + \frac{(x - x_{2i})(x - x_{2i+2})}{(x_{2i+1} - x_{2i})(x_{2i+1} - x_{2i+2})} f(x_{2i+1}) \cr & \qquad + \frac{(x - x_{2i})(x - x_{2i+1})}{(x_{2i+2} - x_{2i})(x_{2i+2} - x_{2i+1})} f(x_{2i+2}) \cr \end{aligned} \tag{5.11}\]

と書ける。 ( 実際、 $p_2(x_{2i}) = f(x_{2i})$ 、 $p_2(x_{2i+1}) = f(x_{2i+1})$ 、 $p_2(x_{2i+2}) = f(x_{2i+2})$ であることを確かめよ。) 分点は $(5.10)$ を満たす等間隔であることを考慮し、 $x = y + x_{2i+1}$ と置くと

\[p_2(x) = \frac{y(y - h)}{2h^2} f(x_{2i}) + \frac{(y + h)(y - h)}{-h^2} f(x_{2i+1}) + \frac{(y + h) y}{2h^2} f(x_{2i+2}) \tag{5.12}\]

$(5.12)$ を $x_{2i}$ から $x_{2i+2}$ まで積分すると、

\[\begin{aligned} \int_{x_{2i}}^{x_{2i+2}} p_2 (x) \ dx & = \int_{-h}^{h} p_2 (y + x_{2i+1}) \ dy \cr & = \frac{h}{3} [\ f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i+2})\ ] \cr \end{aligned}\]

故に

\[\boxed{ \begin{aligned} \quad & \cr & \text{シンプソンの公式} \cr & \quad I_2 = \frac{h}{3} [\ f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i+2})\ ] \cr & \quad \qquad \text{ただし} \quad h = x_{2i+2} - x_{2i+1} = x_{2i+1} - x_{2i} \cr & & \quad \end{aligned} } \tag{5.13}\]

すなわち、分点は $x_{2i}$ 、 $x_{2i+1}$ 、 $x_{2i+2}$ 、重みは

\[w_0 = w_2 = \frac{1}{3}, w_1 = \frac{4}{3} \tag{5.14}\]

である。

$a$ から $b$ までの $2N$ 小区間の全積分の近似値は、

\[\begin{aligned} \int_{a}^{b} p_2 (x) dx & = \sum_{i=0}^{N-1} \int_{x_{2i}}^{x_{2i+2}} p_2 (y + x_{2i+1}) \ dy \cr & = \frac{h}{3} \sum_{i=0}^{N-1} [ f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i+2}) ] \cr \end{aligned}\]

より、各小区間の上限は隣の小区間の下限と一致するから

\[\boxed{ \begin{aligned} \quad \cr & \text{シンプソンの公式 (複合公式)} \cr & \quad I_2 = \frac{h}{3} \left[ f(a) + f(b) + 4\sum_{i=0}^{N-1} f(x_{2i+1}) + 2 \sum_{i=1}^{N-1} f(x_{2i}) \right] \cr & \quad \qquad \text{ただし} \quad h = \frac{b - a}{2N}, \quad x_i = a + ih \ (0 \leqq i \leqq 2N) \cr & & \quad \end{aligned} } \tag{5.15}\]

となる。 すなわち、全区間の上限と下限の重みは $1$ 、奇数番目の分点の重みは $4$ 、(上限と下限をのぞく) 偶数番目の分点の重みは $2$ にして、加え合わせればよい。

---

例題 5.1

下限を $a$ 、上限を $b$ 、積分範囲を $H = b - a$ とする。 $a$ から $b$ までの全区間を $N = 2^k$ 個の小区間に分割して、台形公式 (複合公式) を用いたときの積分の近似値を $T^k$ とする。 $k = 0$ からはじめて、 $k = 1, 2, 3, \cdots$ と半分割を進めていき、

\[\vert T^{k+1} - T^k \vert < \varepsilon \vert T^{k+1} \vert \qquad \varepsilon \text{ は許容相対誤差}\]

となったとき、 $T^{k+1} を積分の近似値とする。 この手順の PAD を書け。

［解］

$k = 0$ のときは台形公式 $(5.6)$ より

\[T^0 = H \frac{f(a) + f(b)}{2}\]

$k = 1$ のときは

\[\begin{alignedat}{2} T^1 & = \frac{H}{2} \left[ \frac{f(a) + f(b)}{2} + f \right. && \left. \left( a + \frac{H}{2} \right) \right] \cr & = \frac{1}{2} T^0 + h f(a + h), && \quad h = \frac{H}{2} \end{alignedat}\]

この第1項は既知であるから、第2項のみ新たに計算すればよい。 一般に

\[T^k = \frac{H}{2^k} \left[ \frac{f(a) + f(b)}{2} + \sum_{i=1}^{2^k-1} f \left( a + i \frac{H}{2^k} \right) \right] \tag{5.16}\]

である。 これと、 $k$ を $k - 1$ で置き換えた式

\[T^{k-1} = \frac{H}{2^{k-1}} \left[ \frac{f(a) + f(b)}{2} + \sum_{i=1}^{2^{k-1} - 1} f \left( a + i \frac{H}{2^{k-1}} \right) \right] \tag{5.17}\]

とを比べると、 $(5.16)$ の総和の $i$ が偶数の項は $(5.17)$ の総和の項に等しいから

\[T^k = \frac{1}{2} T^{k-1} + h \sum_{i=\text{奇数}}^{2^k -1} f(a + ih) \qquad h = \frac{H}{2^k} \tag{5.18}\]

の関係がある。 各 $k$ では第2項のみ新たに計算し、前回の値の $1/2$ に加えれば、新しい値が求められる。 この手順を PAD で表すと、図 5.3 となる。

このように、台形公式をはじめとする閉型公式では、 $h$ の半分割を繰り返しながらよい精度の積分近似値を求めるとき、新しい分点での $f(x)$ のみ計算するだけで、旧分点での値を計算し直す必要がない。

---

例題 5.2

前問の台形公式の代わりに、シンプソンの公式 (複合公式) $(5.15)$ を用いたときの PAD を書け。

［解］

シンプソンの公式 (複合公式) $(5.15)$ を検討すると、新規に $f(x)$ の値を計算しなければならない分点の重みは $4$ 、すでに $f(x)$ の値がわかっている分点の重みは $2$ である。 このことを使って PAD を書くと図 5.4 が得られる。

---

5.1.4 中点公式

開型公式のもっとも簡単な公式は 中点公式 である。 全区間 $[\ a,\ b\ ]$ を $N$ 個の中区間に分割し、各中区間を2分割して2つの小区間をつく。 各小区間の幅を

\[h = \frac{b - a}{2N} \tag{5.19}\]

と置く。 そして各小区間の中点を分点とする (図 5.5 参照)。 すなわち、第 $i$ 中区間の分点は

\[x_{2i+1} = a + (2i + 1) h \qquad (0 \leqq i \leqq N - 1) \tag{5.20}\]

第 $i$ 中区間 $[\ x_{2i},\ x_{2i+2}\ ]$ では $f(x)$ は区間の中点での値 $f(x_{2i+1})$ であると近似すると、 $x$ 軸と $f(x)$ とによってつくられるこの区間内の矩形の面積 (積分近似値) は、

\[\boxed{ \begin{aligned} \quad \cr & \text{中点公式} \cr & \quad I_0 = 2h f(x_{2i+1}) \quad\text{ただし} \quad h = \frac{b - a}{2N} \cr & & \quad \cr \end{aligned} } \tag{5.21}\]

と得られる。 したがって、全区間 $[\ a,\ b\ ]$ については

\[\boxed{ \begin{align*} \quad \cr & \text{中点公式 (複合公式)} \cr & \quad I_0 = 2h \sum_{i=0}^{N-1} f(x_{2i+1}) \quad \text{ただし} \quad h = \frac{b - a}{2N} \quad \tag{5.22} \cr \end{align*} }\]

を積分の近似値とする。

開型公式 (複合公式) の場合は、小区間同士の分点の重なりはない。

5.2 補間多項式

5.2.1 ラグランジュの補間多項式

前節 5.1 節 で導いた閉型公式 (台形公式、シンプソンの公式) は、被積分関数 $f(x)$ を $x$ の1次式 $p_1(x)$ 、2次式 $p_2(x)$ で近似して、 $f(x)$ の代わりにこれら多項式を積分して導いた。 この節では、任意の関数 $f(x)$ を多項式で近似する方法および近似による誤差をやや一般的に考察しよう。

閉区間 $[\ x_0,\ x_n\ ]$ 内に相異なる分点 $x_0 < x_1 < x_2 < \cdots < x_{n-1} < x_n$ をとる。 これらの $n + 1$ 分点で $f(x)$ と一致する $n$ 次多項式を

\[p_n (x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n \tag{5.23}\]

はただ一通りしかない。 何故なら、分点上で $f(x)$ と $p_n (x)$ は一致するから、

\[a_0 + a_1 x_i + a_2 x_i^2 + \cdots + a_n x_i^n = f( x_i ) \qquad (0 \leqq i \leqq n) \tag{5.24}\]

である。 これを $n + 1$ 個の未知数 $a_0,\ a_1,\ \cdots ,\ a_n$ についての連立1次方程式とみなすと、係数の行列式はファンデルモンド (Vandermonde) の行列式 (問題 5-1 参照)

\[\det \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \cr 1 & x_1 & x_1^2 & \cdots & x_1^n \cr 1 & x_2 & x_2^2 & \cdots & x_2^n \cr \vdots & \vdots & \vdots & \ddots & \vdots \cr 1 & x_n & x_n^2 & \cdots & x_n^n \cr \end{pmatrix} = \prod_{i>j} (x_i - x_j) \tag{5.25}\]

であり、分点はすべて相異なるからこの行列式は $0$ でない。 故に、連立1次方程式の解 $a_0, a_1, \cdots , a_n$ は1組しかなく、多項式 $p_n(x)$ は1つしかない。 解 $a_0, a_1, \cdots , a_n$ を連立1次方程式を解いて求めて $(5.23)$ に代入して $p_n (x)$ を求めてもよいが、 $p_n (x)$ はただ1つしかないから次のようにしてよりスマートに求める。 まず、 $(n + 1)$ 次の多項式

\[\phi_{n+1}(x) = (x - x_0)(x - x_1) \cdots (x - x_n) = \prod_{i=0}^{n} (x - x_i) \tag{5.26}\]

を導入する。 明らかに

\[\phi_{n+1}(x_i) = 0 \qquad (0 \leqq i \leqq n) \tag{5.27}\]

である。 その導関数 $\phi_{n+1}^\prime (x)$ は

\[\phi_{n+1}^\prime (x) =\sum_{i=0}^{n} \prod_{\substack{k = 0 \cr k \neq i}}^{n} (x - x_k) \tag{5.28}\]

であり、したがって

\[\phi_{n+1}^\prime (x_i) = \prod_{\substack{k = 0 \cr k \neq i}}^{n} ( x_i - x_k ) \tag{5.29}\]

である。 いま

\[l_i (x) = \frac{\phi_{n+1}(x)}{(x - x_i) \phi_{n+1}^\prime (x_i)} = \prod_{\substack{k = 0 \cr k \neq i}}^{n} \frac{x - x_k}{x_i - x_k} \tag{5.30}\]

と置くと、 $l_i (x)$ は $x = x_i$ で $1$ 、 $x = x_k \neq x_i$ では $0$ である。 したがって

\[\boxed{ \begin{alignat*}{1.5} \quad \cr & \text{ラグランジュの補間多項式} \cr & \quad p_n(x) = \sum_{i=0}^{n} l_i (x) f(x_i) & \quad \tag{5.31} \cr \end{alignat*} }\]

と置くと、 $p_n (x)$ は $(5.24)$ を満たす $n$ 次多項式である。 ちなみに、 $n = 1$ のときは $(5.5)$ と一致し、また $n = 2$ のときは $(5.11)$ と一致する。 一般に、 $n + 1$ 個の $x$ で一致する $n$ 次多項式は (表し方はいろいろあるかもしれないが) 1つしか存在しないから、 $(5.31)$ の $p_n (x)$ は $(5.24)$ を満たすその $n$ 次多項式である。 $(5.31)$ の $p_n (x)$ を $n$ 次の ラグランジュ (Lagrange) の補間多項式 という。

5.2.2 差分商

ラグランジュの補間多項式は $(5.31)$ にしたがって求められる。 そのほかに $p_0 (x) = f(x_0)$ より始まって、 $p_1 (x)$ 、 $p_2 (x)$ 、 $\cdots$ と次つぎに導く 漸化式 ^$1)$ (recurrence formula) がある。 この漸化式を導く準備として、次のような 差分商 (divided difference) の列を考えよう。

$^{1)}$ ゼンカシキと読む。

\[\begin{align*} & f[x_0] = f(x_0) \tag{5.32} \cr & f[x_1, x_0] = \frac{f[x_1] - f[x_0]}{x_1 - x_0} = \frac{f(x_0)}{x_0 - x_1} + \frac{f(x_1)}{x_1 - x_0} = f[x_0, x_1] \qquad \qquad \qquad \tag{5.33} \cr & f[x_2, x_1, x_0] = \frac{f[x_2, x_1] - f[x_1, x_0]}{x_2 - x_0} \cr \end{align*}\] \[= \frac{f(x_0)}{(x_0 - x_1)(x_0 - x_2)} + \frac{f(x_1)}{(x_1 - x_0)(x_1 - x_2)} + \frac{f(x_2)}{(x_2 - x_0)(x_2 - x_1)}\] \[= f[x_2, x_0, x_1] = f[x_1, x_2, x_0] = f[x_1, x_0, x_2] = \cdots = f[x_0, x_1, x_2] \tag{5.34}\]

一般に

\[\begin{align*} & f[x_n, x_{n-1}, \cdots , x_1, x_0] \cr & \quad = \frac{f[x_n, x_{n-1}, \cdots , x_1] - f[x_{n-1}, \cdots , x_1, x_0]}{x_n - x_0} \cr & \quad = \sum_{i=0}^{n} \frac{f(x_i)}{(x_i - x_0) \cdots (x_i - x_{i-1})(x_i - x_{i+1}) \cdots (x_i - x_n)} \cr & \quad = \sum_{i=0}^{n} \frac{f(x_i)}{\phi_{n+1}^\prime (x_i)} \tag{5.35} \cr \end{align*}\]

差分商 $f[x_n, x_{n-1}, \cdots x_1, x_0]$ は変数 $x_0, x_1, x_2, x_3, \cdots , x_n$ の順番を変えても変わらない。 すなわち、差分商はその変数について対称式である。 このことは $(5.32)$ 、 $(5.33)$ 、 $(5.34)$ 、 $(5.35)$ から明かである。

［例 5.1］

差分商の性質を理解するために、簡単な $f(x)$ について差分商を見てみよう。

\[\begin{aligned} f(x) = &\ x \text{のとき} \cr \cr & f[x_0] = f(x_0) = x_0 \cr & f[x_1, x_0] = \frac{f[x_1] - f[x_0]}{x_1 - x_0} = \frac{x_1 - x_0}{x_1 - x_0} = 1 \cr \cr f(x) = &\ x_2 \text{のとき} \cr \cr & f[x_0] = f(x_0) = x_0^2 \cr & f[x_1, x_0] = \frac{f[x_1] - f[x_0]}{x_1 - x_0} = \frac{x_{1}^{2} - x_{0}^{2}}{x_1 - x_0} = x_1 + x_0 \cr & f[x_2, x_1, x_0] = \frac{f[x_2, x_1] - f[x_1, x_0]}{x_2 - x_0} = \frac{x_2 - x_0}{x_2 - x_0} = 1 \cr \cr f(x) = &\ x_3 \text{のとき} \cr \cr & f[x_0] = x_{0}^{3} \cr & f[x_1, x_0] = \frac{f[x_1] - f[x_0]}{x_1 - x_0} = \frac{x_{1}^{3} - x_{0}^{3}}{x_1 - x_0} = x_{1}^{2} + x_1 x_0 + x_{0}^{2} \cr & f[x_2, x_1, x_0] = \frac{f[x_2, x_1] - f[x_1, x_0]}{x_2 - x} = x_2 + x_1 + x_0 \cr & f[x_3, x_2, x_1, x_0] = 1 \cr \end{aligned}\]

差分商は、その変数について対称式である。 たとえば、 $f(x) = x_3$ のとき、

\[\begin{aligned} & f[x_0, x_1] = f[x_1, x_0] = x_{1}^{2} + x_1 x_0 + x_{0}^{2} \cr & f[x_0, x_1, x_2] = f[x_0, x_2, x_1] = \cdots = f[x_2, x_1, x_0] = x_2 + x_1 + x_0 \cr & f[x_0, x_1, x_2, x_3] = f[x_0, x_1, x_3, x_2] = \cdots = f[x_3, x_2, x_1, x_0] = 1 \cr \end{aligned}\]

であり、 $x_0$ 、 $x_1$ 、 $x_2$ 、 $x_3$ の交換に対して不変である。

5.2.3 前進差分

特に $x_0, x_1, x_2, \cdots, x_n$ が等間隔であるときは、

\[x_i = x_0 + ih \qquad (0 \leqq i \leqq n)\]

である。 このとき、

\[\varDelta f(x) \equiv f(x + h) - f(x) \tag{5.36}\]

を 前進差分 (forward difference) という。 前進差分の前進差分は

\[\begin{aligned} \varDelta (\varDelta f (x)) & = \varDelta (f(x + h) - f(x)) \cr & = \varDelta f(x + h) - \varDelta f(x) \cr \end{aligned}\]

である。 これを $\varDelta^2 f(x)$ と書き 第2階前進差分 と呼ぶ。 くわしく書くと、

\[\begin{aligned} \varDelta^2 f(x) & = (f(x + 2h) - f(x + h)) - (f(x + h) - f(x)) \cr & = f(x + 2h) - 2f(x + h) + f(x) \cr \end{aligned}\]

である。 これらの前進差分は差分商によって

\[\begin{alignedat}{2} & \varDelta f(x_0) & & = f(x_1) - f(x_0) = h f[x_1, x_0] \cr & \varDelta^2 f(x_0) & & = \varDelta f(x_1) - \varDelta f(x_0) \cr & & & = h f[x_2, x_1] - h f[x_1, x_0] = 2 h^2 f[x_2, x_1, x_0] \cr \end{alignedat}\]

とも書ける。 一般に第 $m$ 階前進差分を

\[\begin{align*} \varDelta^m f(x) & = \varDelta (\varDelta^{m-1} f(x)) \cr & = \varDelta^{m-1} f(x + h) - \varDelta^{m-1} f(x) \tag{5.37} \cr \end{align*}\]

と定義すれば、 \varDelta mf(x_0) = m ! hmf[ x_0, x_1, x_2, \cdots , xm ] \tag{5.38} であることがわかる。

なお、前進差分の他に 後退差分 (backward difference) や 中心差分 (centered dif- ference) もある。

\[\begin{alignat*}{3} & \nabla f(x) & & \equiv f(x) - f(x - h) & \qquad & \text{後退差分} \tag{5.39} \cr & \delta f(x) & & \equiv f \left( x + \frac{1}{2} h \right) - f \left( x - \frac{1}{2} h \right) & \qquad & \text{中心差分} \tag{5.40} \cr \end{alignat*}\]

後退差分は、第6章常微分方程式の数値解法のとき、とくに重要になる。 このとき、 $(5.38)$ に対応する関係式は

\[\begin{align*} \nabla^m f(x) & = \nabla (\nabla^{m-1} f(x)) \cr & = \nabla^{m-1} f(x) - \nabla^{m-1} f(x - h) \tag{5.41} \cr \end{align*}\]

と $m$ 階後退差分を定義すれば、

\[\nabla^m f(x_m) = m!\ h^m f[x_0, x_1, x_2, \cdots , x_m] \tag{5.42}\]

である。

5.2.4 差分商の拡張

差分商は、 $(5.32)$ から $(5.35)$ のように、変数 $x_0, x_1, \cdots, x_n$ に等しいものがあれば、分子も分母も $0$ となり不定形になるように見える。 しかし、[例 5.1] の場合に見られるように、差分商は対称式であり、変数を等しくしてみると、 $f[x, x] = f^\prime (x)$ 、 $f[x, x, x] = f^{\prime\prime}(x) / 2!$ が成り立っている。 以下、

\[f[\underbrace{x, x, \cdots , x, x, x}_ {k+1 \text{ 個}}] = \frac{f^{(k)} (x)}{k!} \tag{5.43}\]

が任意の $k$ に対して成立することを証明する。 いま、ある $k$ にたいして成立しているとする。 $k + 1$ に対して $(5.43)$ は

\[\begin{align*} f[\underbrace{x, x, \cdots , x, x, x}_ {k+2 \text{ 個}} ] & = \lim_{h \to 0} \frac{1}{k + 1} \sum_{i=1}^{k+1} f[\underbrace{x + h, \cdots , x + h}_ {k+2-i \text{ 個}} , \underbrace{x, \cdots , x}_ {i \text{ 個}} ] \cr & = \frac{1}{k + 1} \lim_{h \to 0} \sum_{i=1}^{k+1} \frac{1}{h} \left[ f[ \underbrace{x + h, \cdots, x + h}_ {k+2-i \text{ 個}}, \underbrace{x, \cdots, x}_ {i-1 \text{ 個}} ] - f[ \underbrace{x + h, \cdots, x + h}_ {k+1-i \text{ 個}}, \underbrace{x, \cdots , x}_ {i \text{ 個}} ] \right] \cr & = \frac{1}{k + 1} \lim_{h \to 0} \frac{1}{h} [ f[ \underbrace{x + h, \cdots, x + h}_ {k+1 \text{ 個}} ] - f[ \underbrace{x, \cdots, x}_ {k+1 \text{ 個}} ] ] \cr & = \frac{1}{(k + 1)!} \lim_{h \to 0} \frac{f^{(k)} (x + h) - f^{(k)}(x)}{h} \cr & = \frac{f^{(k+1)}(x)}{(k + 1)!} \tag{5.44} \cr \end{align*}\]

である。 故に、数学的帰納法より $(5.43)$ は任意の $k$ に対して成立する。

一般に、 $D = d / dx$ と置くと、

\[f[x_0, x_1, x_2, \cdots , x_n, \underbrace{x, x, \cdots , x}_ {k+1 \text{ 個}} ] = \frac{1}{k!} D^k f[x_0, x_1, x_2, \cdots , x_n, x] \tag{5.45}\]

が任意の $x$ の (微分可能な) 関数 $f(x)$ に対して成立することが、 $(5.43)$ と同様に証明できる。 その証明のためには、 $(5.43)$ の $f(x)$ を $f[x_0, x_1, x_2, \cdots , x_n, x]$ で置き換えるだけでよい。 $(5.45)$ はさらに $x_0, x_1, \cdots x_n$ のそれぞれが複数個存在する場合に拡張できることは明らかであろう (エルミートの補間公式 $(5.60)$ 参照)。

また、差分商をいくつかの変数が等しい場合でも成立するように

\[\begin{align*} f & [x_n, x_{n-1}, \cdots , x_1, x_0] \cr & = \int_0^1 \int_0^{t_1} \cdots \int_0^{t_{n-1}} f^{(n)} \left( x_0 + \sum_{k=1}^{n} t_k (x_k - x_{k-1}) \right) dt_n \cdots dt_2 dt_1 \tag{5.46} \cr \end{align*}\]

と定義することができる (問題 5-3 参照)。

5.2.5 ニュートンの補間公式

以上の差分商を用いれば、 $p_n (x)$ は次のように表される。

\[\boxed{ \begin{align*} \quad \cr & \text{ニュートンの補間公式} \cr & \quad \begin{alignat*}{5} p_n (x) & = & & f[x_0] + (x - x_0) f[x_0, x_1] + \cdots \cr & & & + (x - x_0)(x - x_1) \cdots (x - x_{n-1})f[x_0, x_1, \cdots , x_n] \cr & = & & p_{n-1} (x) + \phi_n (x) f[x_0, x_1, \cdots , x_{n-1}, x_n] \tag{5.47} \cr \end{alignat*} \cr \end{align*} \quad }\]

これを ニュートンの補間公式 という。 この公式は、ある $x$ に対する $p_n (x)$ を $f(x_i)$ $(0 \leqq i \leqq n)$ から数値的に求めるのによく使われる。 特に等間隔分点のときは、差分商を $(5.38)$ から前進差分で表して

\[\begin{align*} p_n (x) = &\ f(x_0) + \frac{x - x_0}{h} \varDelta f(x_0) + \cdots \cr & + \frac{(x - x_0)(x - x_1) \cdots (x - x_{n-1})}{n!\ h^n} \varDelta ^n f(x_0) \tag{5.48} \cr \end{align*}\]

である。 次の 例題 5.3 で公式 $(5.47)$ の証明をしよう。

---

例題 5.3

ニュートンの補間公式 $(5.47)$ を証明せよ。

［解］

この第1式を証明するために、 $n$ 次多項式 $p_n (x)$ を

\[\begin{align*} p_n (x) = c_0 & + c_1(x - x_0) + c_2(x - x_0)(x - x_1) + \cdots \cr & + c_n(x - x_0)(x - x_1) \cdots (x - x_{n-1}) \tag{5.49} \cr \end{align*}\]

と表し、 $p_n (x_i) = f(x_i) \ (0 \leqq i \leqq n)$ となるように $c_i$ を定めよう。 まず、 $x = x_0$ と置けば、 $c_0 = f(x_0) = f[x_0]$ が得られる。 これを $(5.49)$ に代入し 変形して

\[\begin{align*} \frac{p_n (x) - f[x_0]}{x - x_0} = c_1 & + c_2(x - x_1) + \cdots \cr & + c_n (x - x_1) \cdots (x - x_{n-1}) \tag{5.50} \cr \end{align*}\]

が得られる。 この式で $x = x_1$ と置くと、 $c_1 = f[x_1, x_0]$ が得られる。 $c_1$ を $(5.50)$ に代入して、

\[\begin{align*} & \frac{p_n (x) - f[x_0] - (x - x_0) f[x_1, x_0]}{(x - x_0)(x - x_1)} \cr & \qquad = c_2 + c_3 (x - x_2) + \cdots + c_n (x - x_2) \cdots (x - x_{n-1}) \tag{5.51} \cr \end{align*}\]

である。 この式で $x = x_2$ と置くと、右辺は $c_2$ 、左辺は $f[x_2, x_1, x_0]$ であるから、 $c_2 = f[x_2, x_1, x_0]$ が得られる。 同様の繰り返しによって、 $(5.47)$ の第1式の証明が完結する。 第1式の最後の項を除いたものが $p_{n-1} (x)$ であるから、第2式が得られる。

なお、 $p_0 (x) = f[x_0] = f(x_0)$ は $x = x_0$ で $f(x)$ と一致する0次式であり、 $p_1 (x) = f(x_0) + (x - x_0) f[x_0, x_1]$ は $x = x_0$ と $x = x_1$ で $f(x)$ と一致する1次式、 $\cdots$ と順に確かめて行くことによっても、 $(5.47)$ は確かめられる。 ( もう一つのラグランジュの補間式の形 ネビユの補間公式 については 問題 5-4 を参照。)

---

［例 5.2］

等間隔4分点における $f(x)$ の値が次のようにわかっているとき、これらの点を通る3次多項式を求めよう。

\[\begin{array}{rrrrr} x & 0 & 1 & 2 & 3 \cr f(x) & -5 & 0 & 6 & 20 \cr \end{array}\]

次のような階差表をつくり、 $\varDelta^m f(0)$ を求め、 $(5.48)$ に代入する。

\[\begin{array}{c|cccc} x & f(x) & \varDelta f(x) & \varDelta^2 f(x) & \varDelta^3 f(x) \cr \hline 0 & \underline{-5} & \underline{5} & \underline{1} & \underline{7} \cr 1 & 0 & 6 & 8 \cr 2 & 6 & 14 \cr 3 & 20 \cr \end{array}\]

この表より、 $f(0) = -5$ 、 $\varDelta f(0) = 5$ 、 $\varDelta^2 f(0) = 1$ 、 $\varDelta^3 f(0) = 7$ 。 そして、 $h = 1$ で等間隔だから、前進差分 $(5.38)$ とニュートンの補間公式 $(5.47)$ より

\[\begin{aligned} p_0(x) & = -5, \frac{}{} \cr p_1(x) & = p_0(x) + (x - 0) \varDelta f (0) = 5x - 5, \frac{}{} \cr p_2(x) & = p_1(x) + \frac{(x - 0)(x - 1)}{2} \varDelta^2 f(0) \cr & = \frac{1}{2} x^2 + \frac{9}{2} x - 5 \cr p_3(x) & = p_2(x) + \frac{(x - 0)(x - 1)(x - 2)}{6} \varDelta^3 f(0) \cr & = \frac{7}{6} x^3 - 3x^2 + \frac{41}{6} x - 5 \cr \end{aligned}\]

なお、ニュートンの補間公式を用いた数値補間の例を第8章図 8.1 に示す。 図 8.1 の例は、補間領域が広すぎて、ルンゲの現象が起こる例である。

5.2.6 補間多項式の誤差

$f(x)$ を $n$ 次多項式 $p_n (x)$ で近似したときの誤差 $e_n (x)$ を

\[e_n (x) \equiv f(x) - p_n (x) \tag{5.52}\]

と定義しよう。 $n = 0$ のときは、

\[\begin{aligned} e_0 (x) & = f(x) - p_0 (x) \cr & = f(x) - f(x_0) \cr & = f[x] - f[x_0] = (x - x_0) f[x_0, x] \cr \end{aligned}\]

と書ける。 $n = 1$ のときは

\[\begin{aligned} e_1 (x) & = f(x) - p_1 (x) \cr & = f[x] - f[x_0] - (x - x_0) f[x_0, x_1] \cr & = (x - x_0)(f[x, x_0] - f[x_0, x_1]) \cr & = (x - x_0)(x - x_1) f[x_0, x_1, x] \cr \end{aligned}\]

同様にして、一般に $n$ 次補間多項式で任意の関数 $f(x)$ 近似したときの誤差は

\[\boxed{ \begin{aligned} \quad \cr & n \text{ 次補間多項式による誤差 } (1) \cr & \quad \begin{align*} e_n (x) & \equiv f(x) - p_n (x) \cr & = \phi_{n+1} (x)f[x_0, x_1, \cdots , x_n, x] \tag{5.53} \cr \end{align*} \quad \end{aligned} }\]

と表される。 ただし、 $\phi_{n+1} (x)$ は $(5.26)$ で定義した $n + 1$ 次多項式である。

$(5.53)$ の分かりやすい表現を得るために、やや天下りであるが、

\[\begin{align*} q_n (t) = &\ f(t) - p_n (t) \cr & - (t - x_0) \cdots (t - x_n) f[x_0, x_1, \cdots , x_n, x] \tag{5.54} \cr \end{align*}\]

なる関数を考える。 $(5.53)$ より明らかに $q_n (x) = 0 \ (x \neq x_i)$ 、かつ、 $q_n (x_i) = 0 \ (0 \leqq i \leqq n)$ である。 以後 $n + 2$ 個の点 $x_0, x_1, \cdots , x_n, x_{n+1}$ を小さい方から順に番号をつけ変え

\[x_0 < x_1 < \cdots < x_n < x_{n+1}\]

とする。 微分積分学のロール (Rolle) の定理により開区間 $(x_i, x_{i+1})$ の中に

\[q_n^\prime (\xi_i^{(1)} ) = 0 \qquad (0 \leqq i \leqq n)\]

であるような点 $\xi_i^{(1)}$ が存在する。 したがって、開区間 $(\xi_i^{(1)}, \xi_{i+1}^{(1)})$ 内に、

\[q_n^{\prime\prime} (\xi_i^{(2)}) = 0 \qquad (0 \leqq i \leqq n - 1)\]

であるような点 $\xi_i^{(2)}$ が存在する。 同様な繰り返しにより、開区間 $(\xi_0^{(n)}, \xi_1^{(n)})$ の中に

\[q_n^{(n+1)} (\xi) = 0 \qquad (x_0 < \xi < x_{n+1}) \tag{5.55}\]

であるような点 $\xi$ が存在する。 $n$ 次多項式の $(n + 1)$ 階微分は $p_n^{(n+1)} (t) = 0$ だか ら、 $(5.54)$ と $(5.55)$ より、差分商は

\[f[x_0, x_1, \cdots, x_n, x] = \frac{f^{(n+1)}(\xi)}{(n + 1)!} \tag{5.56}\]

とも表される。 ここの $\xi = \xi (x)$ は $x$ の関数であり、 $\xi (x)$ は $x_0 < \xi < x_{n+1}$ の範 囲にある。 ( もちろん、この式が成り立つのは、 $f(x)$ の $n + 1$ 階の導関数が存在する場合に限る。) これを $(5.47)$ に代入して考えればニュートンの補間公式は $f(x)$ のテーラー展開に似た形をしていることがわかる。 ( 特に、 $x_0 = x_1 = \cdots = x_n$ なら、 $f(x)$ のニュートンの補間公式はテーラー展開そのものである。) また、 $(5.53)$ の $e_n (x)$ は

\[\boxed{ \begin{align*} \quad \cr & n \text{ 次補間多項式による誤差 } (2) & \quad \cr & \quad \begin{align*} e_n (x) & \equiv f(x) - p_n (x) \cr & = \phi_{n+1} (x) \frac{f^{(n+1)} (\xi)}{(n + 1)!} \tag{5.57} \cr \end{align*} \end{align*} }\]

とも表される。 $\xi$ は $x$ の関数だから、 $f^{(n+1)} (\xi)$ は $x$ の関数であるが、 $f(x)$ が $n$ 次以下の多項式なら $f^{(n+1)} (\xi) = 0$ だから、 $e_n (x) = 0$ である (誤差はない)。

5.2.7 エルミート の補間公式

$n$ 次のラグランジュの補間公式 $(5.31)$ は、 $n + 1$ 点で関数 $f(x)$ と一致する $n$ 次の多項式である。 さらにその上に、この $n + 1$ 点で1階導関数が $f^\prime (x)$ と一致するような多項式は、より高次の多項式でなければならないだろう。 関数値だけのときは、 $(5.24)$ の $n + 1$ 個の条件を満たせばよいが、導関数の値までが一致するためには、さらに $n + 1$ 個の導関数値についての条件が加えられるから、合計 $2n + 2$ 個の条件があり、この多項式は $2n + 1$ 次式になるであろう。

このような多項式の表現の仕方はいろいろあるにしても、多項式としてはやはり1つしかなく、表現の1つは $(5.30)$ の $l (x)$ を用いて

\[\begin{align*} p_{2n+1}^\ast (x) = & \sum_{i=0}^{n} l_i (x)^2 {1 - 2l_i^\prime (x_i)(x - x_i)}f(x_i) \cr & + \sum_{i=0}^{n} l_i (x)^2 (x - x_i)f^\prime (x_i) \tag{5.58} \cr \end{align*}\]

なる $2n + 1$ 次式がある。 この式は確かに $x = x_i$ において

\[p_{2n+1}^{\ast}(x_i) = f(x_i) \quad \text{および} \quad p_{2n+1}^{\ast\prime}(x_i) = f^\prime (x_i) \tag{5.59}\]

を満足する。 ( 読者自ら確かめよ。)

$(5.58)$ は1階の導関数を含むから、差分商で表せば、変数が2つずつ等しい差分商によって表されることが予想される。 実際に、次のように表される。

\[\boxed{ \begin{aligned} \quad \cr & \text{エルミートの補間公式 } (n + 1 \text{ 点、 } 2n + 1 \text{ 次多項式}) \cr & \quad \begin{align*} & p_{2n+1}^{\ast}(x) \cr & = f[x_0] + (x - x_0) f[x_0, x_0] + (x - x_0)^2 f[x_0, x_0, x_1] \cr & \quad + (x - x_0)^2 (x - x_1) f[x_0, x_0, x_1, x_1] \cr & \qquad \cdots \qquad \cdots \cr & \quad + (x - x_0)^2 \cdots (x - x_{n-1})^2 (x - x_n) f[x_0, x_0, \cdots , x_n, x_n] \tag{5.60} \cr \end{align*} \quad \end{aligned} }\]

このエルミート (Hermite) の補間公式の証明は次の例題で行う。

---

例題 5.4

エルミートの補間公式 $(5.60)$ の $p_{2n+1}^\ast (x)$ は、 $(5.59)$ を満たす $2n + 1$ 次多項式であることを示せ。

［解］

一般性を失うことなく、 $x_0 < x_1 < \cdots < x_n$ とする。 $n + 1$ 点を間につけ加えて $x_0 < \xi_0 < x_1 < \xi_1 < \cdots < x_n < \xi_n$ とする。 $(5.47)$ より、この $2n + 2$ 点を通る $2n + 1$ 次の多項式

\[\begin{aligned} & p_{2n+1} (x) = f[x_0] + (x - x_0) f[x_0, \xi_0] + \cdots \cr & \qquad + (x - x_0)(x - \xi_0) \cdots (x - \xi_{n-1})(x - x_n)f[x_0, \xi_0, \cdots x_n, \xi_n] \tag{5.61} \cr \end{aligned}\]

を使って、 $(5.53)$ より

\[f(x) = p_{2n+1} (x) + \prod_{i=0}^{n} (x - x_i)(x - \xi_i) \cdot f[x_0, \xi_0, x_1, \xi_1, \cdots , x_n, \xi_n, x] \tag{5.62}\]

である。 ここで、 $\xi_i \to x_i \ (0 \leqq i \leqq n)$ とすると、 $p_{2n+1} (x) \to p_{2n+1}^\ast (x)$ が得られる。

また $(5.62)$ は

\[f(x) = p_{2n+1}^\ast (x) + \phi_{n+1} (x)^2 f[x_0, x_0, x_1, x_1, \cdots x_n, x_n, x] \tag{5.63}\]

となる。 ここに $\phi_{n+1} (x)$ は $(5.26)$ で導入した $n + 1$ 次式である。

$(5.63)$ の左辺は $x_i \ (0 \leqq i \leqq n)$ によらないし、右辺第2項は $x_i \ (0 \leqq i \leqq n)$ の交換をしても変わらない。 したがって $p_{2n+1}^\ast (x)$ も $x_i \ (0 \leqq i \leqq n)$ の入れ換えに対して変わらない。 いま、 $x_0$ と $x_i$ とを入れ換えれば、 $(5.60)$ より

\[p_{2n+1}^\ast (x) = f[x_i] + (x - x_i)f[x_i, x_i] + (x - x_i)^2 g(x)\]

と書ける。 $g(x)$ は $x$ の $2n - 1$ 次式である。 これより

\[\begin{aligned} & p_{2n+1}^\ast (x_i) = f[x_i] = f(x_i) \cr & p_{2n+1}^{\ast\prime}(x_i) = f[x_i, x_i] = f^\prime (x_i) \cr \end{aligned} \tag{5.64}\]

である。

---

$(5.63)$ はエルミートの補間による誤差が

\[\boxed{ \begin{aligned} \quad \cr & \text{エルミートの補間公式による誤差} \cr & \quad \begin{aligned} e_{2n+1} (x) & \equiv f(x) - p_{2n+1}^\ast (x) \cr & = \phi_{n+1} (x)^2 f[x_0, x_0, x_1, x_1, \cdots x_n, x_n, x] \cr & = \frac{\phi_{n+1} (x)^2}{(2n + 2)!} f^{(2n+2)} (\xi) x_0 < \xi < x_n \cr \end{aligned} \end{aligned} \quad } \tag{5.65}\]

であることを示している。

5.3 ニュートン・コーツの公式

5.3.1 ニュートン・コーツの閉型公式

5-1 節で導いた等間隔分点の閉型公式 (台形公式、シンプソンの公式) を一般化しよう。 全区間 $[ a, b ]$ を $N$ 等分し、おのおのを中区間とよぶことにする。 各中区間を $n$ 個の小区間に $n$ 等分する。 小区間の幅は、

\[h = \frac{b - a}{nN} \tag{5.66}\]

である。 第 $j$ 中区間 $(0 \leqq j \leqq N - 1)$ 内の分点は

\[x_{ij} = a + (jn + i) h \qquad (0 \leqq i \leqq n) \tag{5.67}\]

の $n + 1$ 個の分点を含んでいる。 第 $j$ 中区間の $n + 1$ 個の分点 $x_{ij}$ において、被積分関数 $f(x)$ と一致する $n$ 次補間多項式はラグランジュの補間多項式 $p_n (x)$ $(5.31)$ である。

\[p_n (x_{ij}) = f(x_{ij}) \qquad (0 \leqq i \leqq n) \tag{5.68}\]

1つの中区間 $[ x_{0j}, x_{nj} ]$ の積分の近似値は、 $f(x)$ の代わりに $p_n (x)$ をこの中区間にわたって積分することによって $(5.2)$ の形に得られる。 すなわち、 $(5.28)$ から $(5.30)$ までの式において、まず、変数変換 $x = x_{0j} + th$ によって、

\[\begin{aligned} & \phi_{n+1} (x) = h^{n+1} \prod_{k=0}^n (t - k) \cr & \phi_{n+1}^\prime (x_{ij}) = h^n \prod_{\substack{k=0 \cr k \neq i}} (i - k) = h^n (-1)^{n-i} (n - i)!\ i! \cr \end{aligned}\]

とし、これらを $(5.31)$ に代入すると、 $(5.2)$ より

\[\boxed{ \begin{aligned} \quad \cr & \text{ニュートン・コーツの公式 (閉型公式)} \cr & \quad I_n = h \sum_{i=0}^n w_i f(x_{ij}) \cr & \quad \qquad w_i = \frac{(-1)^{n-i}}{(n - i)!\ i!} \int_0^n \prod_{\substack{k=0 \cr k \neq i}} (t - k) \ dt \cr & \quad \qquad h = \frac{b - a}{nN}, \qquad x_{ij} \text{ は } (5.67) & \quad \end{aligned} } \tag{5.69}\]

が得られる。 これを ニュートン・コーツ (Newton - Cotes) の閉型公式 という。 $(5.69)$ は、第 $j$ 中区間 $[ x_{0j}, x_{nj} ]$ についての積分公式であるから、全区間 $[ a, b ]$ については、中区間について加え合わせて、次の複合閉型公式が得られる。

\[\boxed{ \begin{aligned} \quad \cr & \text{ニュートン・コーツの公式 (複合閉型公式)} \cr & \quad I_n = h \sum_{i=0}^{n} w_i \sum_{j=0}^{N-1} f(x_{ij}) \cr & \quad \qquad w_i,\ h,\ x_{ij} \text{ は } (5.69) \text{ と同じ} & \quad \cr \end{aligned} } \tag{5.70}\]

［例 5.3］

$(5.69)$ の重み $w_i$ は、 $n = 1$ のときは台形公式、 $n = 2$ のときはシンプソンの公式と一致することを確かめる。

\[\begin{aligned} & n = 1 & & w_0 = \frac{-1}{0!\ 1!} \int_{0}^{1} (t - 1) \ dt = \frac{1}{2}, \quad w_1 = \frac{1}{1!\ 0!} \int_{0}^{1} t \ dt = \frac{1}{2} = w_0 \cr & n = 2 & & w_0 = \frac{1}{2!\ 0!} \int_{0}^{2} (t - 1)(t - 2) \ dt = \frac{1}{3}, \quad w_1 = -\frac{1}{1!\ 1!} \int_{0}^{2} t (t - 2) \ dt = \frac{4}{3}, \cr & & & w_2 = \frac{1}{0!\ 2!} \int_{0}^{2} t(t - 1) \ dt = \frac{1}{3} = w_0 \end{aligned}\]

［例 5.4］

$n = 3$ のときのニュートン・コーツの公式を シンプソンの 3/8 公式 という。 この公式を求めてみると、重みは

\[w_0 = w_3 = \frac{3}{8}, \qquad w_1 = w_2 = \frac{9}{8} \tag{5.71}\]

である。

ニュートン・コーツの公式 $(5.69)$ は、中区間の $f(x)$ を $n$ 次の多項式で近似して得られた公式であるから、もし $f(x)$ が $n$ 次以下の関数であれば、(丸めの誤差をのぞいて) 正確なはずである。 すなわち、 $(5.69)$ の打切り誤差は $h^{n+1}$ の程度かもっと小さいかのはずである。 (このことから、積分しないで $w_i$ を求めることが出来る。問題 5-5 参照)

打切り誤差は $(5.3)$ より

\[E_n = I - I_n = \int \{f(x) - p_n (x)\} \ dx = \int e_n (x) \ dx \tag{5.72}\]

である。 $(5.53)$ より、1つの中区間 $[ x_0, x_n ]$ では

\[E_n = \int_{x_0}^{x_n} \phi_{n+1} (x) \ f [x_0, x_1, \cdots x_n, x] \ dx \tag{5.73}\]

この積分に、積分の平均値の第2定理を適用する。 この定理が成立するためには、被積分関数には、積分領域内では定符号の部分が存在する必要がある。 定符号の部分は、部分積分によって作り出すが、 $n$ が偶数のときと奇数のときとでは異なる。 結果は、 $n$ が偶数のときは、

\[E_n = \frac{f^{(n+2)} (\xi)}{(n + 2)!} \int_{x_0}^{x_n} x \phi_{n+1} (x) \ dx \qquad n \text{ は偶数} \tag{5.74}\]

$n$ が奇数のときは、

\[E_n = \frac{f^{(n+1)} (\xi)}{(n + 1)!} \int_{x_0}^{x_n} \phi_{n+1} (x) \ dx \qquad n \text{ は奇数} \tag{5.75}\]

が得られる (問題 5-7 参照)。 ここに $\xi$ は、 $x_0 < \xi < x_n$ なるある定数である。

$x_i (0 \leqq i \leqq n)$ は等間隔であるから、変数変換 $x = x_0 + ht$ および $(5.26)$ の $\phi_{n+1} (x)$ の定義により、次の結果を得る。

\[\boxed{ \begin{align*} \quad \cr & \text{ニュートン・コーツの公式の打切り誤差 (閉型公式)} & \quad \cr & \quad \begin{align*} E_n = \frac{f^{(n+2)} (\xi)}{(n + 2)!} h^{n+3} \int_{0}^{n} t^2 \prod_{i=1}^{n}(t - i) \ dt & \qquad & n : \text{偶数} \tag{5.76} \cr = \frac{f^{(n+1)} (\xi)}{(n + 1)!} h^{n+2} \int_{0}^{n} t \prod_{i=1}^{n}(t - i) \ dt & \qquad & n : \text{奇数} \tag{5.77} \cr \end{align*} \end{align*} }\]

これらは、中区間にわたる打切り誤差である。 $E_n$ は、 $n$ が偶数の時も奇数の時も、 $h$ の奇数べき乗であることに注意せよ (問題 5-9 も参照)。 全区間にわたる複合公式の打切り誤差はこの式を加えることによって求められる。

$n$ の小さい場合の重み $w_i (0 \leqq i \leqq n \leqq 10)$ と、打切り誤差 $E_n$ の係数を公式 $(5.69)$ および $(5.76)$ 、 $(5.77)$ によって求めた結果を表 5-1 に示す。 ( また、 問題 5-5 も参照せよ。) ただし、重みの対称性 $w_{n-i} = w_i$ より、 $w_i (i \geqq n/2)$ は一部省略している。 重み $w_i$ は $w_i = A W_i$ の形に共通因数 $A$ をくくりだし、 $W_i$ は整数になるように表した。 また、「 $E_n$ の係数」とは $E_n = \gamma_m h^{m+1} f^{(m)}(\xi)$ ( $n$ が奇数の時 $m = n + 1$ 、または $n$ が偶数の時 $m = n + 2$ ) と表したときの $\gamma_m$ である。 この表からわかるように、重みの符号が負になるものが現れる $n \geqq 8$ の公式は、桁落ちによる誤差の生じる可能性があるので、滅多に用いられない。

\[\begin{array}{c} \begin{array}{c||l||r|r|r|r|r|r|c} \hline n & A & W_0 & W_1 & W_2 & W_3 & W_4 & W_5 & E_n \text{ の係数} \\ \hline \hline 1 & 1/2 & 1 & 1 & & & & & -1/12 \cr 2 & 1/3 & 1 & 4 & 1 & & & & -1/90 \cr 3 & 3/8 & 1 & 3 & 3 & 1 & & & -3/80 \cr 4 & 2/45 & 7 & 32 & 12 & 32 & 7 & & -8/945 \cr 5 & 5/288 & 19 & 75 & 50 & 50 & 75 & 19 & -275/12096 \cr 6 & 1/140 & 41 & 216 & 27 & 272 & 27 & 216 & -9/1400 \cr 7 & 7/17280 & 751 & 3577 & 1323 & 2989 & 2989 & 1323 & -8183/518400 \cr 8 & 4/14175 & 989 & 5888 & -928 & 10496 & -4540 & 10496 & -2368/467775 \cr 9 & 9/89600 & 2857 & 15741 & 1080 & 19344 & 5778 & 5778 & -4671/394240 \cr \end{array} \cr \text{表 } 5-1 \text{ ニュートン・コーツの閉型公式の} \cr \text{重み } w_i = A W_i \text{ と打切り誤差 } E_n \text{ の係数} \cr \end{array}\]

5.3.2 ニュートン・コーツの開型公式

閉型公式についての以上の考察は、あまり変更なしに、開型公式に当てはまる。 開型公式の場合は、各小区間の端点は分点ではないことが大きな違いである。 いま、全区間 $[ a, b ]$ を $N$ 中区間に分ける。 $n$ 次式で近似するために各中区間内に端点ではない $n + 1$ 個の分点を置く。 すなわち各中区間を $n + 2$ 等分する。 こうして得た小区間の幅は、

\[h = \frac{b - a}{(n + 2) N} \tag{5.78}\]

である。 第 $j$ 中区間の分点は、この中区間の端点を除く点で、

\[x_{ij} = a + [ j(n + 2) + i + 1 ] h \qquad (0 \leqq j \leqq N - 1,\ 0 \leqq i \leqq n) \tag{5.79}\]

である。 いま、第 $j$ 中区間 $[ x_{-1j}, x_{n+1\ j} ]$ 内の $n + 1$ 個の分点 $x_{0j}, x_{1j} , \cdots , x_{nj}$ で $f(x)$ と一致する $n$ 次式 $p_n (x)$ $(5.31)$ をつくり、 $f(x)$ の代わりに用いる。 故に

\[w_i = \frac{1}{h} \int_{x_{-1j}}^{x_{n+1\ j}} \phi_{n+1} (x) (x - x_{ij}) \phi_{n+1}^\prime (x_{ij}) \ dx \tag{5.80}\]

$x = x_{ij} + ht$ と変数変換をすると、一つの中区間について

\[\boxed{ \begin{align*} \quad \cr & \text{ニュートン・コーツの公式 (開型公式)} \cr & \quad I_n = h \sum_{i=0}^{n} w_i f(x_{ij}) \tag{5.81} \cr & \quad \qquad w_i = \frac{(-1)^{n-i}}{(n - i)!\ i!} \int_{-1}^{n+1} \prod_{\substack{k=0 \cr k \neq i}}^n (t - k) \ dt \cr & \quad \qquad h = \frac{b - a}{(n + 2) N} , \qquad x_{ij} \text{ は } (5.79) \cr \end{align*} }\]

［例 5.5］

$n = 0$ のときの中点公式では

\[w_0 = \frac{(-1)^0}{0!\ 0!} \int_{-1}^{1} 1 \cdot dt = 2\]

が得られる。

複合公式は、 $(5.81)$ を加え合わせて

\[\boxed{ \begin{align*} \quad \cr & \text{ニュートン・コーツの公式 (複合開型公式)} & \quad \cr & \quad I_n = h \sum_{i=0}^{n} w_i \sum_{j=0}^{N-1} f(x_{ij} ) \tag{5.82} \cr & \quad \qquad w_i,\ h,\ x_{ij} \text{ は } (5.81) \text{ と同じ} \cr \end{align*} }\]

である。 また、打切り誤差は閉型公式とほぼ同じ (積分範囲に注意) であり、結果は

\[\boxed{ \begin{align*} & \text{ニュートン・コーツの公式の打切り誤差 (開型公式)} \cr & \quad \begin{align*} E_n & = \frac{f^{(n+2)} (\xi)}{(n + 2)!} h^{n+3} \int_{-1}^{n+1} t^2 \prod_{i=1}^{n} (t - i) \ dt & n : \text{偶数} \tag{5.83} \cr & = \frac{f^{(n+1)} (\xi)}{(n + 1)!} h^{n+2} \int_{-1}^{n+1} t \prod_{i=1}^{n} (t - i) \ dt & n : \text{奇数} \tag{5.84} \cr \end{align*} \end{align*} }\]

である。 閉型公式と同じく、 $E_n$ は $h$ の奇数べき乗であることに注意。

開型公式の重み $(5.81)$ と打切り誤差 $(5.83)$ 、 $(5.84)$ を $n$ の小さな場合を表 5-2 に示す。 重みは $n \geqq 2$ で負の値が現れるから、開型公式は小さい $n$ のときにすでに桁落ちによる誤差の懸念が生じる。 したがって、中点公式 $(n = 0)$ 以外の開型公式は閉型公式と比べてとくにメリットはない。 なお、 $E_n$ の係数はすべて正であることが、閉型公式と対照的である。

\[\begin{array}{c} \begin{array}{c||l||r|r|r|r|r|c} \hline n & A & W_0 & W_1 & W_2 & W_3 & W_4 & E_n \text{ の係数} \cr 0 & 2 1 & 1/3 & \cr 1 & 3/2 & 1 & 1 & 3/4 \cr 2 & 4/3 & 2 & -1 & 2 & 14/45 \cr 3 & 5/24 & 11 & 1 & 1 & 11 & 95/144 \cr 4 & 3/10 & 11 & -14 & 26 & -14 & 11 & 41/140 \cr 5 & 7/1440 & 611 & -453 & 562 & 562 & -453 & 5257/8640 \cr 6 & 8/945 & 460 & -954 & 2196 & -2459 & 2196 & 3956/14175 \cr \end{array} \cr \text{表 } 5-2 \text{ ニユートン・コーツの開型公式の} \cr \text{重み } w_i = A W_i \text{ と打切り誤差 } E_n \text{ の係数} \cr \end{array}\]

5.4 ロンバーグの積分法

5.4.1 ロンバーグの積分法

閉型公式の打切り誤差 $E_n$ は、被積分関数を近似する多項式 $p_n (x)$ の次数を $n$ とすれば、 $h^{n+2}$ ( $n$ は奇数 $(5.77)$ ) または $h^{n+3}$ ( $n$ は偶数 $(5.76)$ ) に比例する。 したがって、 $h$ が小さいほど打切り誤差は小さい。 一方、丸めの誤差は $h$ が小さいほど(演算回数が $h$ に反比例して増えるから) 大きい。 実際、 $h$ を $h = H = b - a$ から出発して半分割を繰り返すと、打切り誤差は小さくなり、丸めの誤差は大きくなる。

その大体の様子を 図 5.6 に示す。 横軸は半分割回数、縦軸は誤差を最初の打切り誤差で割った値の対数である。 2種類の誤差の和は、半分割を進めていくと次第に小さくなるが、ある半分割回数を越すと、丸めの誤差が大きくなるためにかえって上昇に転ずる。 この半分割回数を 最適半分割回数 とよぶことにする。 最適半分割回数をこえて半分割を進めることは意味がない。 むしろ高次の公式 (近似多項式の次数 $n$ の大きな公式) を使うべきである。 この節では、半分割を繰り返しながら、同時に高次の公式を自動的に使うようにした数値積分法を考察しよう。

台形公式は、近似多項式の次数 $n = 1$ の公式であり最低次の公式であるが、 $h = H = b - a$ からはじめて、 $k$ 回の半分割を行ったときの近似積分値を $T_0^k$ と表わそう (ここでの $T_0^k$ は 5-2 節の 例題 5.1 の $(5.16)$ の $T^k$ と全く同じものである)。 すなわち、 $h = H = H / 2^0$ のときは

\[T_0^0 = \frac{H}{2} (f(a) + f(b)) \tag{5.85}\]

であり、以下 $k$ 回の半分割を $k = 1, 2, 3, \cdots$ と繰り返したときの近似積分値 $T_0^k$ は、 $h = H / 2^k$ として、 $(5.18)$ と同じく

\[T_0^k = \frac{1}{2} T_0^{k-1} + h \sum_{j=\text{奇数}}^{2^k - 1} f(a + jh) \tag{5.86}\]

である。

さて、 $T_0^{k+1}$ と $T_0^k$ とを使って、次の $T_1^k$ を求める：

\[T_1^k = T_0^{k+1} + \frac{T_0^{k+1} - T_0^k}{4^1 - 1}\]

同様に、 $m = 1, 2, 3, \cdots$ の順に、 $T_{m-1}^k$ と $T_{m-1}^{k+1}$ とを使って $T_m^k$ を求める：

\[T_m^k = T_{m-1}^{k+1} + \frac{T_{m-1}^{k+1} - T_{m-1}^k}{4^m - 1} \tag{5.87}\]

この $m$ を順次大きくする手続きを、 補外 (extrapolation) と呼ぶ。 $k$ が大きくなり半分割が進めば、 $(T_{m-1}^{k+1} - T_{m-1}^k)$ は小さくなり、 $m$ が大きくなれば $(4^m - 1)$ が大きくなるから、やがて $T_m^k$ と $T_{m-1}^{k+1}$ が一致するようになり、そのとき計算を停止する。

$m + k$ を縦に、 $m$ を横に並べた $T_m^k$ の表を ロンバーグの表 (Romberg table) と呼ぶ。 表 5-3 にロンバーグの表を示す。 この表の各要素と $(5.86)$ 、 $(5.87)$ を見比べて、表の作成法を考えよう。 まず台形公式から $T_0^0$ と $T_0^1$ を求め、この $T_0^0$ と $T_0^1$ から $T_1^0$ を求める。 次に $T_0^2$ を台形公式から求め、 $T_0^1$ と $T_0^2$ から $T_1^1$ を、 $T_1^0$ と $T_1^1$ から $T_2^0$ を求める。 以下同様に、最左の列は半分割して台形公式から求めて、2列目以降は補外によって左隣と左上隣の要素から $(5.87)$ によって求める。 左隣の要素との値の差が十分小さくなったら、その値を積分値とするというわけである。

\[\begin{array}{c} \begin{array}{|c||c|c|c|c|c|c|} \hline m + k & m = 0 & m = 1 & m = 2 & m = 3 & \cdots & m = n \cr \hline\hline 0 & T_0^0 \cr 1 & T_0^1 & T_1^0 \cr 2 & T_0^2 & T_1^1 & T_2^0 \cr 3 & T_0^3 & T_1^2 & T_2^1 & T_3^0 \cr \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \cr n & T_0^n & T_1^{n-1} & T_2^{n-2} & T_3^{n-3} & \cdots & T_n^0 \cr \hline \end{array} \cr \text{表 } 5-3 \text{ ロンバーグの表} \end{array}\]

さて、補外によってつくられた $T_m^k$ $(m > 0)$ は何を表すのであろうか。 $T_0^k$ $(m = 0)$ は台形公式で $h = H / 2^k$ として得られた積分の近似値であった。 その打切り誤差は

\[\varDelta T_0^k = c_1 h^2 + c_2 h^4 + c_3 h^6 + \cdots \tag{5.88}\]

のように $h^2$ の (漸近) 級数であることが証明される (問題 5-8 参照)。 ここに $c_1, c_2, c_3, \cdots$ は $h$ によらない $0$ でない定数である。 このような場合には、半分割をさらに1回行ったときの $T_0^{k+1}$ の打切り誤差は、同じ定数 $c_1, c_2, c3, \cdots$ を用いて ( $h$ を $h / 2$ で置き換えて)

\[\varDelta T_0^{k+1} = c_1 \left( \frac{h}{2} \right)^2 + c_2 \left( \frac{h}{2} \right)^4 + c_3 \left( \frac{h}{2} \right)^6 + \cdots \tag{5.89}\]

と書ける。 したがって $T_1^k$ の打切り誤差は

\[\begin{align*} \varDelta T_1^k & = \varDelta T_0^{k+1} + \frac{\varDelta T_0^{k+1} - \varDelta T_0^k}{4^1 - 1} \cr & = c_2^\prime h^4 + c_3^\prime h^6 + \cdots \tag{5.90} \cr \end{align*}\]

となり、 $h_2$ の項は消えてしまう。 ここに $c_2^\prime, c_3^\prime, \cdots$ も $h$ によらない定数である ( $c_{l+1}^\prime = (1 - 4^l) / (3 \cdot 4^l) c_{l+1}$ )。 同様にして、$T_m^k$ の打切り誤差は $(5.87)$ より $h^{2(m+1)}$ の項から始まることが示せる ( $c_{m+l}^\prime = (1 - 4^l) / (4^l (4^m - 1)) c_{m+l}$ )。 こうして、 $m$ が大きいほど高次の近似多項式による積分公式を使ったことになる。 実際、 $T_1^k$ は2次のシンプソンの公式、 $T_2^k$ は4次のニュートン・コーツの公式に一致する ( $k > 0$ なら複合公式)。 $m > 2$ では、ニュートン・コーツの公式に対応するものはなく、例えば公式 $T_3^0$ は、 $n = 2(m + 1) =$ 8 次公式であり、係数は

\[A = \frac{4}{2835}, \quad W_0 = 217, \quad W_1 = W_3 = 1024, \quad W_2 = 352, \quad W_4 = 436, \quad W_i = W_{8-i}\]

そして打ち切り誤差は $E_8 = -1/338688 \cdot h^9 f^{(8)}(\xi)$ であることが確かめられる (表 5-1 と対照せよ)。 いずれにせよ $T_m^k$ の打ち切り誤差は $h^{2m+3}$ のオーダーである。

補外の計算は $(5.87)$ に見られるように簡単であり、わずかな計算で高次の積分値が得られるのが、ロンバーグの積分法の特徴である。

$T_0^0$ と $T_0^1$ を台形公式で求めておき、1回だけ補外を行い $T_1^1$ を積分近似値とする積分法が古くから用いられており、 リチャードソン (Richardson) の補外法 と呼ばれている。

5.4.2 ロンバーグ積分法の手順

表 5-3 のロンバーグ表と $(5.87)$ を参照して手順をつくる。 $T_{m-1}^k$ と $T_{m-1}^{k+1}$ から $T_m^k$ を求めたら、以後 $T_{m-1}^k$ は参照されることはない。 したがって、 $T_m^k$ は $T_{m-1}^k$ と同じ記憶場所に記憶して置けばよい。 すなわち、すべての $T_m^k$ は $T_0^k$ に記憶させる。

PAD ではすべての $T_m^k$ は T(k) と記してある。

丸めの誤差が大きくなると、補外は意味がなくなる。 何故なら丸めの誤差は $h^2$ の級数では表せないから、丸めの誤差の大きい $T_m^k$ に対して補外をおこなっても誤差は小さくならないからである。 そこである $m$ の時、 $T_{m-1}^{k+1}$ と $T_m^k$ の差が小さくなったら、以後、この $m$ 以上に補外は行わないことにし、 $k$ を増やす ( $h$ を半分割する) ことだけを行う。

PAD ではすべての $T_m^k$ は T(k) に記憶されているから、 T(k) と T(k+1) の差が小さくなったら、この k+1 を変数 kmin に覚えておき、以後補外は k $\geqq$ kmin までとし、$h$ の半分割ごとに kmin を $1$ づつ増やす (ロンバーグの表の下方に進む) ことにしてある。

ロンバーグ 表の大きさは十分大きくとっておかないと、収束する前に足らなくな ることがある (図 5.7 の PAD では MXHLF )。 $H$ が $1$ の程度の時は、 $f(x)$ が $\sin x$ のような初等超越関数でも $10$ くらいであるから、 $20$ もあれば十分である (MXHLF=20)。

$H$ は $1$ の程度であることが望ましい。 積分範囲が大きいときには、これをさらに小区間に区分しておいて、各々の小区間に対して、ロンバーグ積分を実行して、加え合わせる必要がある。

［例 5.6］

円周率 $\pi$ を求める積分

\[\int_0^1 \frac{4}{1 + x^2} \ dx = \pi\]

を倍精度計算 (10 進 15 桁) の精度で行った結果、次のようなロンバーグ表を得た。

\[\begin{array}{c} \begin{array}{c||c|c|c} m + k & m = 0,3,6 & m = 1,4 & m = 2,5 \cr \hline \hline 0 & 3.000000000000000 \cr 1 & 3.100000000000000 & 3.133333333333334 \cr 2 & 3.131176470588236 & 3.141568627450980 & 3.142117647058824 \cr 3 & 3.138988494491090 & 3.141592502458707 & 3.141594094125889 \cr & 3.141585783761874 \cr 4 & 3.140941612041389 & 3.141592651224823 & 3.141592661142564 \cr & 3.141592638396796 & 3.141592665277718 \cr 5 & 3.141429893174975 & 3.141592653552837 & 3.141592653708037 \cr & 3.141592653590030 & 3.141592653649611 & 3.141592653638244 \cr 6 & 3.141551963485657 & 3.141592653589217 & 3.141592653591642 \cr & \underline{3.141592653589793} & \underline{3.141592653589793} & 3.141592653589735 \cr & 3.141592653589723 \cr 7 & 3.141582481063753 & 3.141592653589785 & 3.141592653589823 \cr & \underline{3.141592653589793} & \underline{3.141592653589793} & \ast \cr \end{array} \cr \text{積分値} = 3.141592653589793 \end{array}\]

$T_3^3$ と $T_4^2$ ( $m + k = 6$ の2つのアンダーライン) が一致したので、以下 $m > 4$ の補外 ( $\ast$ 印以下) は行わないことにして、半分割を進めたところ、 $T_3^4$ と $T_4^3$ ( $m + k = 7$ の2つのアンダーライン) が一致した。 そこで、 $T_4^3$ を積分近似値とした。 この値は倍精度計算で達し得る最高の精度を持っている。 補外を行わないときの値は $T_0^7$ であるが、5桁の精度しか持っていない。 したがって補外によって精度は10桁上がったことになる。

5.5 ガウス型積分公式

5.5.1 ガウス型積分公式

積分公式は $n$ 分点の場合、 $n$ 個の分点と $n$ 個の重みがある。 ニュートン・コーツの公式のときは $n$ 個の分点は等間隔に決められていた。 したがって、精度を高くする (打切り誤差を小さくする) ためには $n$ 個の重みだけが選べ、近似多項式の次数は $n - 1$ 次多項式であった。 もし等間隔という制限をなくして、重みと分点の $2n$ 個を近似多項式の次数をあげ、精度を高くするように選べば、 $2n - 1$ 次まで可能なはずである。 この節では、分点の座標と重みを、近似多項式の次数が $2n - 1$ 次であるように選ぶ公式を導こう。

求める積分は、重み関数 $w(x)$ をふくめて

\[I = \int_a^b w(x) f(x) dx \tag{5.91}\]

とする。 ( ニュートン・コーツの公式のときは $w(x) = 1$ である。) この積分範囲の中に $n$ 個の分点を置く。 その分点の座標 $x_1, x_2, \cdots , x_n$ を

\[a < x_1 < x_2 < \cdots < x_n < b\]

とする。 ( この分点座標は以下で定める。) ニュートン・コーツの公式のときと同じように、 $f(x)$ とこの $n$ 個の分点で一致する $n - 1$ 次多項式 $p_{n-1} (x)$ は

\[p_{n-1} (x) = \sum_{i=1}^n \frac{\phi_n(x)}{(x - x_i) \phi_n^\prime (x_i)} f(x_i) \tag{5.92}\]

である。 ただし、

\[\phi_n(x) = \prod_{i=1}^n (x - x_i) \tag{5.93}\]

積分公式は

\[I_n = \int_a^b w(x) p_{n-1} (x) \ dx = \sum_{i=1}^n w_i f(x_i) \tag{5.94}\]

ここに、

\[w_i = \int_a^b \frac{w(x) \phi_n(x)}{(x - x_i) \phi_n^\prime (x_i)} \ dx \tag{5.95}\]

また、打切り誤差は、 $(5.53)$ より

\[\begin{align*} E_n & = \int_a^b { p_{n-1} (x) - f(x) } w(x) \ dx \cr & = \int_a^b \phi_n (x) f[x_1, x_2, \cdots , x_n, x] w(x) \ dx \tag{5.96} \end{align*}\]

したがって、分点が不等間隔であり、 $(5.95)$ の $w$ に $(1/h)$ の係数がつかないこと除けば、ニュートン・コーツの公式の場合と同じである。

以下、分点を決めていくことになる。 ところで分点は $\phi_n(x)$ の零点であるから、分点を決めることは $\phi_n(x)$ を決めることになる。 $\phi_n(x)$ として次に述べる直交多項式をとったときの積分公式を ガウス型積分公式 という。

5.5.2 直交多項式

$n$ 次の ルジャンドル (Legendre) の多項式 $P_n (x)$ は、

\[P_n (x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} (x^2 - 1)^n \qquad (-1 \leqq x \leqq 1) \tag{5.97}\]

で定義される $n$ 次多項式である。 例えば、

\[P_0(x) = 1, \qquad P_1(x) = x, \qquad P_2(x) = \frac{1}{2} (3x^2 - 1), \qquad \cdots\]

などである。 $(5.97)$ から、低い次数のルジャンドル多項式から高い次数のルジャンドル多項式を求める漸化式

\[(n + 1) P_{n+1} (x) = (2n + 1) x P_n (x) - n P_{n-1} (x) \tag{5.98}\]

を導くことができる。 すなわち、 $P_0 (x) = 1$ 、 $P_1 (x) = x$ を与えると、 $P_n (x)\ (n > 1)$ を次々と順に求められる。 多項式の数値計算法であるホーナー法によるより $(5.98)$ によって $P_n (x)$ を求めるほうが、演算量も少なく精度もよい。 また、この多項式の導関数は、

\[(x^2 - 1)P_n^\prime (x) = n (x P_n (x) - P_{n-1} (x)) \tag{5.99}\]

によって求められる。 この関係式は同じく $(5.97)$ から導き出せる。

ルジャンドルの多項式は、次の直交関係を満たすので、直交多項式の1つである。

\[\int_{-1}^1 P_n (x) P_m (x) \ dx = \frac{2}{2n + 1} \delta_{nm} \tag{5.100}\]

ここに $\delta_{nm}$ はクロネッカーのデルタである。 この直交関係は、数値積分にとって重要な関係であることが間もなくわかる。

このような直交多項式としては、他に次のようなものがある。

ラゲル (Laguerre) の多項式

\[\begin{align*} & L_n (x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (e^{-x} x^n) \tag{5.101} \cr & (n + 1) L_{n+1} (x) = -(x - 2n - 1) L_n (x) - n L_{n-1} (x) \tag{5.102} \cr & x L_n^\prime (x) = n(L_n (x) - L_{n-1} (x)) \tag{5.103} \cr & \int_0^\infty e^{-x} L_n (x) L_m (x) \ dx = \delta_{nm} \tag{5.104} \cr & L_0 (x) = 1, \quad L_1 (x) = -x + 1, \quad L_2 (x) = \frac{x^2}{2} - 2x + 1, \quad \cdots \cr \end{align*}\]

エルミート (Hermite) の多項式

\[\begin{align*} & H_n (x) = (-1)^n e^{x^2} \frac{d^n}{dx^n} e^{-x^2} \tag{5.105} \cr & H_{n+1} (x) = 2x H_n (x) - 2n H_{n-1} (x) \tag{5.106} \cr & H_n^\prime (x) = 2n H_{n-1} (x) \tag{5.107} \cr & \int_{-\infty}^\infty e^{-x^2} H_n (x) H_m (x) \ dx = 2^n n! \sqrt{\pi} \delta_{nm} \tag{5.108} \cr & H_0 (x) = 1, \quad H_1 (x) = 2x, \quad H_2 (x) = 4x^2 - 2, \cdots \end{align*}\]

直交関係 $(5.104)$ と $(5.108)$ は、積分範囲が無限大であることと、 $e^{-x}$ とか $e^{-x^2}$ という $x$ が無限大ですばやく $0$ に収束する重み関数 $w(x)$ を持つことが特徴である。

5.5.3 ガウスの積分公式

いま、近似多項式 $p_{n-1} (x)$ $(5.92)$ の $\phi_n (x)$ $(5.93)$ としてルジャンドルの多項式をとる。 ただし、ルジャンドルの多項式は変数の範囲が $[ -1, 1 ]$ で定義されているから、変数変換

\[x = \frac{b - a}{2} t + \frac{a + b}{2} \tag{5.109}\]

によって、 $a \leqq x \leqq b$ を $-1 \leqq t \leqq 1$ に変換する。 さらに、 $\phi_n (x)$ の $x_n$ の係数は $1$ だから、正確には

\[\phi_n (x) = \frac{2^n (n!)^2}{(2n)!} \left[ \frac{b - a}{2} \right]^n P_n (t) \tag{5.110}\]

ととる。 このようにとったことは、分点は $P_n (t)$ の零点 ( $P_n (t) = 0$ となる $t$ ) に決めたことになる。 すなわち、分点 $x_i (1 \leqq i \leqq n)$ は

\[x_i = \frac{b - a}{2} t_i + \frac{a + b}{2} \qquad \text{ただし} \quad P_n (t_i) = 0 \tag{5.111}\]

で、 $t_i$ は零点である。 零点は実数で $n$ 個あることが知られている。 こうして、次の公式が得られる。 これを ガウス・ルジャンドルの積分公式 、または簡単に ガウスの積分公式 という。 $(5.112)$ における $w_i$ は、 $P_n (t)$ の性質を使って解析的に積分した結果である (後述の $(5.117)$ を参照)。 また、 $E_n$ はこの公式の打切り誤差である (すぐ後の $(5.115)$ を参照)。

\[\boxed{ \begin{align*} \quad \cr & \text{ガウスの積分公式} & \quad \cr & \quad \int_a^b f(x) \ dx = \frac{b - a}{2} \sum_{i=1}^n w_i f(x_i) + E_n \tag{5.112} \cr & \quad \qquad x_i = \frac{b - a}{2} t_i + \frac{a + b}{2} \qquad (P_n (t_i) = 0) \cr & \quad \qquad w_i = \frac{2(1 - t_i^2)}{n^2 P_{n-1}^2 (t_i)} \qquad (0 \leqq i \leqq n) \cr \end{align*} }\]

5.5.4 ガウスの積分公式の打切り誤差

打切り誤差は $(5.96)$ 、 $(5.110)$ から得られる。 まず $(5.100)$ の意味するところは、 $m < n$ なら

\[\int_{-1}^1 P_n (x) \cdot [\ m \text{ 次多項式 } ] \ dx = 0 \qquad m < n \tag{5.113}\]

ということである。 なぜなら、 $x$ のべきは

\[\begin{alignedat}{2} & x^0 & = &\ 1 = P_0 (x) \cr & x & = &\ P_1 (x), \cr & x^2 & = &\ \frac{2}{3} P_2 (x) + \frac{1}{3} P_0 (x), \cr & x^3 & = &\ \frac{2}{5} P_3 (x) + \frac{3}{5} P_1 (x), \cr & & &\ \cdots \cr \end{alignedat}\]

と繰り返せば、一般に $x^m$ は $P_m (x), P_{m-1} (x), \cdots, P_0 (x)$ の1次結合で表せるから、 $m$ 次多項式も $m$ 次および $m$ 次以下のルジャンドル多項式の1次結合で表せる。 故に、 $(5.100)$ から $(5.113)$ が得られる。

打切り誤差 $E_n$ は $(5.96)$ で $w(x) = 1$ として与えられる。 $(5.110)$ より $g(x) \equiv f[x_1, x_2, \cdots , x_n, x]$ の $x$ について $n - 1$ 次以下の項は、 $(5.113)$ より $(5.96)$ に寄与しない。 寄与するのは $n$ 次以上の項である。 $g(x)$ の $x$ について $n$ 次以上の項とは、 $g(x)$ を $n - 1$ 次式 $p_{n-1} (x)$ で近似したときの「誤差」 $e_{n-1}$ と見なせる。 $(5.53)$ の $f(x)$ を $g(x)$ で置き換えて考えれば、この「誤差」 $e_{n-1}$

\[\begin{aligned} e_{n-1} & = f[x_1, x_2, \cdots , x_n, x] - p_{n-1} (x) \cr & = g(x) - p_{n-1} (x) \cr & = \phi_n (x)g[x_1, x_2, \cdots , x_n, x] \cr & = \phi_n (x)f[x_1, x_1, x_2, x_2, \cdots , x_n, x_n, x] \cr \end{aligned}\]

である。 したがって、

\[f[x_1, x_2, \cdots , x_n, x] = p_{n-1} (x) + \phi_n (x)f[x_1, x_1, x_2, x_2, \cdots , x_n, x_n, x]\]

である。 これを $(5.96)$ に代入すると、 $p_{n-1} (x)$ は寄与しないから、積分の平均値の定理を適用すれば

\[E_n = \frac{f^{(2n)} (\xi)}{(2n)!} \int_a^b \phi_n^2 (x) \ dx \tag{5.114}\]

が得られる。 $\phi_n (x)$ として $(5.110)$ をとって

\[\boxed{ \begin{align*} \quad \cr & \text{ガウスの積分公式の打切り誤差} & \quad \cr & \quad E_n =\frac{2^{2n+1} (n!)^4}{(2n + 1){(2n)!}^3} \left[ \frac{b - a}{2} \right]^{2n+1} f^{(2n)} (\xi) \tag{5.115} \cr & \qquad \qquad \qquad \qquad \text{ただし} \quad a < \xi < b \cr \end{align*} }\]

$f(x)$ が $2n - 1$ 次以下の多項式なら、ガウスの積分公式は (丸めの誤差を除いて) 正確であることがわかる。

5.5.5 ガウスの積分公式の分点と重み

ガウスの積分公式の欠点は分点と重みが簡単ではないことである。 分点や重みを数値的に求めるには次のようにすればよかろう。 分点 $t_i$ はルジャンドルの多項式 $P_n (x)$ の零点であるが、すべて実数で単根である。 そこで分点は

\[t_{i0} = \cos \left( \frac{\pi}{2} \frac{4i - 1}{2n + 1} \right) \qquad (1 \leqq i \leqq n) \tag{5.116}\]

を初期値とする実数のニュートン法により求める。 このとき、零点 $t_i$ は絶対値の等しい正負の値 $t_i = -t_{n-i}$ であり、 $n$ が奇数のときは、 $t_{(n+1)/2} = 0$ である。 したがって、いつでも、正の零点だけを求めればよい。 このとき必要なルジャンドルの多項式 $P_n (x)$ は、 $P_0 (x) = 1$ 、 $P_1 (x) = x$ として、漸化式 $(5.98)$ を用いる。 さらに導関数 $P_n^\prime (x)$ が必要だが、 $P_n (x)$ を求める際に必要である $P_{n-1} (x)$ を使って $(5.99)$ より計算できる。 重みは $(5.112)$ の $w_i$ である。 以下に4倍精度で求めた分点数 $n$ が $2 \leqq n \leqq 5$ の分点と重みの結果の表を掲げる (小数点以下21桁目を4捨5入)。

\[\begin{array}{|c||l|l|} \hline n & \text{分点 } (t_i) & \text{重み } (w_i) \cr \hline \hline 2 & \pm\ 0.57735\ 02691\ 89625\ 76451 & 1.0 \cr \cr 3 & \pm\ 0.77459\ 66692\ 41483\ 37704 & 0.55555\ 55555\ 55555\ 55556 \cr & \quad 0.0 & 0.88888\ 88888\ 88888\ 88889 \cr \cr 4 & \pm\ 0.86113\ 63115\ 94052\ 57522\ & 0.34785\ 48451\ 37453\ 85737 \cr & \pm\ 0.33998\ 10435\ 84856\ 26480\ & 0.65214\ 51548\ 62546\ 14263 \cr \cr 5 & \pm\ 0.90617\ 98459\ 38663\ 99280\ & 0.23692\ 68850\ 56189\ 08751 \cr & \pm\ 0.53846\ 93101\ 05683\ 09104\ & 0.47862\ 86704\ 99366\ 46804 \cr & \quad 0.0 & 0.56888\ 88888\ 88888\ 88889 \cr \hline \end{array}\]

ただし 実用的には、4倍精度は必要でなく倍精度でよく、許容絶対誤差を $\varepsilon = 10^{-15}$ とすれば、数回 ( $n \leqq 100$ でたかだか4回) の反復で収束する。 分点および重みの計算は正の分点のみについて行えばよく、負の分点は正の分点と絶対値が等しく、重みは同じである。 なお、 $\varepsilon = 10^{-16}$ とすると収束しなかったり、収束しても小数点以下15桁目が正しくない場合がある。 これは、使用計算機の丸めの誤差によって異なる。 4倍精度で求めた上の表を参照して精度を調べておいた方がよかろう。 いづれにせよ、 $\varepsilon = 10^{-15}$ で実用的には満足できる。

ガウスの積分公式を用いて、相対誤差が $\varepsilon ($ eps $)$ 以下の積分値を得る手順を図 5.8 の PAD に示す。

5.5.6 ガウス型積分公式の重み

ガウス型積分公式の重みは、 $(5.95)$ で与えられている。 直交多項式の性質を利用すれば、この重みの積分はガウス型積分公式の重み

\[w_i = - \frac{A_{n+1} \gamma_n}{A_n Q_{n+1} (x_i) Q_n^\prime (x_i)} \qquad (1 \leqq i \leqq n) \tag{5.117}\]

で与えられる (問題 5-10 参照)。 ここで、 $Q_n (x)$ はルジャンドルの多項式のような $n$ 次の直交多項式で、直交関係

\[\int_a^b w(x) Q_n (x) Q_m (x) = \delta_{nm} \gamma_n \tag{5.118}\]

を満たすものである。 $A_n$ は $Q_n (x)$ の $x_n$ の係数で、 $Q_n (x) = A_n \phi_n (x)$ と書ける。

［例 5.7］

ガウスの積分公式の場合には、 $Q_n (x) = P_n (x)$ 、 $w(x) = 1$ である。 $(5.98)$ の両辺の $x_{n+1}$ の係数をくらべて

\[A_{n+1} = \frac{2n + 1}{n + 1} A_n \tag{5.119}\]

$(5.98)$ より

\[Q_{n+1} (x_i) = P_{n+1} (x_i) = - \frac{n}{n + 1} P_{n-1} (x_i) \tag{5.120}\]

$(5.99)$ より

\[Q_n^\prime (x_i) = P_n^\prime (x_i) = \frac{n}{1 - x_i^2} P_{n-1} (x_i) \tag{5.121}\]

$(5.100)$ より

\[\gamma_n = \frac{2}{2n + 1} \tag{5.122}\]

これらを $(5.117)$ に代入すると $(5.112)$ の wi が得られる。

なお、 $(5.112)$ の wi は変数変換 $(5.109)$ 後の値であることに注意すること。 そうでなくては、全体にかかる係数 $(b - a) / 2$ は現れない。

5.5.7 ガウス型積分公式の打切り誤差

ガウス型積分公式の打切り誤差は $(5.114)$ に対応して、次のように与えられる。

\[\boxed{ \begin{align*} \quad && \quad \cr & \text{ガウス型積分公式の打切り誤差} \cr & \quad E_n = \frac{\gamma_n f^{(2n)} (\xi)}{(2n)! A_n^2} \tag{5.123} \end{align*} }\]

ここに $\xi$ は積分範囲中のある座標値である。

---

例題 5.5

直交多項式としてラゲルの多項式 $(5.101)$ を用い、次のガウス・ラゲルの積分公式を導け。

\[\boxed{ \begin{align*} \quad && \quad \cr & \text{ガウス・ラゲルの積分公式} \cr & \quad \int_0^\infty e^{-x} f(x) \ dx = \sum_{i=1}^n w_i f(x_i) + E_n \tag{5.124} \cr & \quad \qquad w_i = \frac{x_i}{[n L_{n-1} (x_i)]^2} \qquad L_n (x_i) = 0\ (\text{零点}) \cr & \quad \qquad E_n = \frac{(n!)^2}{(2n)!} f^{(2n)} (\xi) \qquad x_0 < \xi < x_1 \cr \end{align*} }\]

［解］

$x_i$ は $L_n (x)$ の零点 (すなわち $L_n (x_i) = 0$ ) とする。

\[\begin{aligned} & (n + 1) A_{n+1} = -A_n, \cr & \gamma_n = 1, \cr & L_{n+1} (x_i) = - \frac{n}{n + 1} L_{n-1} (x_i), \cr & L_n^\prime (x_i) = - \frac{n}{x_i} L_{n-1} (x_i) \cr \end{aligned}\]

---

例題 5.6

直交多項式としてエルミートの多項式 $(5.105)$ を用い、次のガウス・エルミートの積分公式を導け。

\[\boxed{ \begin{align*} \quad & & \quad \cr & \text{ガウス・エルミートの積分公式} \cr & \quad \int_{-\infty}^\infty e^{-x^2} f(x) \ dx = \sum_{i=1}^n w_i f(x_i) + E_n \tag{5.125} \cr & \quad \qquad w_i = \frac{2^{n-1} n! \sqrt{\pi}}{[nH_{n-1} (x_i)]^2} \qquad H_n(x_i) = 0 \quad (x_i \text{ は零点}) \cr & \quad \qquad E_n = \frac{n! \sqrt{n}}{2^n (2n)!} f^{(2n)} (\xi) \qquad x_0 < \xi < x_n \cr \end{align*} }\]

［解］

$x_i$ は $H_n (x)$ の零点 (すなわち $H_n (x_i) = 0$ ) とする。

\[\begin{aligned} & A_{n+1} = 2 A_n, \cr & \gamma_n = 2^n n! \sqrt{\pi} \cr & H_{n+1} (x_i) = -2n H_{n-1} (x_i), \cr & H_n^\prime (x_i) = 2n H_{n-1} (x_i) \cr \end{aligned}\]

---

例題 5.7

ガウス・ラゲルの積分公式およびガウス・エルミートの積分公式について、分点数 $n \leqq 6$ について、分点を小数点以下20桁、重みを有効数字21桁求めよ。

［解］

\[\begin{array}{c} \text{ガウス・ラゲルの積分公式の分点と重み} \cr \begin{array}{|c||r|c|} \hline n & \text{分点 } (t_i) & \text{重み } (w_i) \cr \hline 2 & 0.58578 ~ 64376 ~ 26904 ~ 95120 & 8.53553 ~ 39059 ~ 32737 ~ 62200 ~ 10^{-1} \cr & 3.41421 ~ 35623 ~ 73095 ~ 04880 & 1.46446 ~ 60940 ~ 67262 ~ 37800 ~ 10^{-1} \cr \cr 3 & 0.41577 ~ 45567 ~ 83479 ~ 08331 & 7.11093 ~ 00992 ~ 91730 ~ 15450 ~ 10^{-1} \cr & 2.29428 ~ 03602 ~ 79041 ~ 71982 & 2.78517 ~ 73356 ~ 92408 ~ 48801 ~ 10^{-1} \cr & 6.28994 ~ 50829 ~ 37479 ~ 19687 & 1.03892 ~ 56501 ~ 58613 ~ 57490 ~ 10^{-2} \cr \cr 4 & 0.32254 ~ 76896 ~ 19392 ~ 31180 & 6.03154 ~ 10434 ~ 16336 ~ 01636 ~ 10^{-1} \cr & 1.74576 ~ 11011 ~ 58346 ~ 57569 & 3.57418 ~ 69243 ~ 77996 ~ 86641 ~ 10^{-1} \cr & 4.53662 ~ 02969 ~ 21127 ~ 98328 & 3.88879 ~ 08515 ~ 00538 ~ 42724 ~ 10^{-2} \cr & 9.39507 ~ 09123 ~ 01133 ~ 12923 & 5.39294 ~ 70556 ~ 13274 ~ 50104 ~ 10^{-4} \cr \cr 5 & 0.26356 ~ 03197 ~ 18140 ~ 91020 & 5.21755 ~ 61058 ~ 28086 ~ 52476 ~ 10^{-1} \cr & 1.41340 ~ 30591 ~ 06516 ~ 79222 & 3.98666 ~ 81108 ~ 31759 ~ 27454 ~ 10^{-1} \cr & 3.59642 ~ 57710 ~ 40722 ~ 08122 & 7.59424 ~ 49681 ~ 70759 ~ 53877 ~ 10^{-2} \cr & 7.08581 ~ 00058 ~ 58837 ~ 55692 & 3.61175 ~ 86799 ~ 22048 ~ 45446 ~ 10^{-3} \cr & 12.64080 ~ 08442 ~ 75782 ~ 65943 & 2.33699 ~ 72385 ~ 77622 ~ 78911 ~ 10^{-5} \cr \cr 6 & 0.22284 ~ 66041 ~ 79260 ~ 68946 & 4.58964 ~ 67394 ~ 99635 ~ 93568 ~ 10^{-1} \cr & 1.18893 ~ 21016 ~ 72623 ~ 03074 & 4.17000 ~ 83077 ~ 21209 ~ 94113 ~ 10^{-1} \cr & 2.99273 ~ 63260 ~ 59314 ~ 07769 & 1.13373 ~ 38207 ~ 40449 ~ 75739 ~ 10^{-1} \cr & 5.77514 ~ 35691 ~ 04510 ~ 50184 & 1.03991 ~ 97453 ~ 14907 ~ 48989 ~ 10^{-2} \cr & 9.83746 ~ 74183 ~ 82589 ~ 91772 & 2.61017 ~ 20281 ~ 49320 ~ 59479 ~ 10^{-4} \cr & 15.98287 ~ 39806 ~ 01701 ~ 78255 & 8.98547 ~ 90642 ~ 96212 ~ 38825 ~ 10^{-7} \cr \hline \end{array} \end{array}\] \[\begin{array}{c} \text{ガウス・エルミート の積分公式の分点と重み} \cr \begin{array}{|c||l|c|} \hline n & \text{分点 } (t_i) & \text{重み } (w_i) \cr \hline 2 & \pm 0.70710 ~ 67811 ~ 86547 ~ 52440 & 8.86226 ~ 92545 ~ 27580 ~ 13649 ~ 10^{-1} \cr \cr 3 & \pm 1.22474 ~ 48713 ~ 91589 ~ 04910 & 2.95408 ~ 97515 ~ 09193 ~ 37883 ~ 10^{-1} \cr & \quad 0.0 & 1.18163 ~ 59006 ~ 03677 ~ 35153 \cr \cr 4 & \pm 1.65068 ~ 01238 ~ 85784 ~ 55588 & 8.13128 ~ 35447 ~ 24517 ~ 71430 ~ 10^{-2} \cr & \pm 0.52464 ~ 76232 ~ 75290 ~ 31788 & 8.04914 ~ 09000 ~ 55128 ~ 36506 ~ 10^{-1} \cr \cr 5 & \pm 2.02018 ~ 28704 ~ 56085 ~ 63293 & 1.99532 ~ 42059 ~ 04591 ~ 32077 ~ 10^{-2} \cr & \pm 0.95857 ~ 24646 ~ 13818 ~ 50711 & 3.93619 ~ 32315 ~ 22411 ~ 59828 ~ 10^{-1} \cr & \quad 0.0 & 9.45308 ~ 72048 ~ 29418 ~ 81226 ~ 10^{-1} \cr \cr 6 & \pm 2.35060 ~ 49736 ~ 74492 ~ 22283 & 4.53000 ~ 99055 ~ 08845 ~ 64086 ~ 10^{-3} \cr & \pm 1.33584 ~ 90740 ~ 13696 ~ 94971 & 1.57067 ~ 32032 ~ 28566 ~ 43916 ~ 10^{-1} \cr & \pm 0.43607 ~ 74119 ~ 27616 ~ 50868 & 7.24629 ~ 59522 ~ 43925 ~ 24092 ~ 10^{-1} \cr \hline \end{array} \end{array}\]

---

例題 5.8

積分の範囲を $N$ 区間に分割し、各区間にガウスの積分公式を適用する。

各区間の分点数 $n$ は $n = 1$ からはじめて順次ふやす。 $n$ のときの積分値と $n - 1$ のときの積分値とを比べて、必要な精度 (許容相対誤差 $\varepsilon$ ) に達したときその区間の積分は終了し、次の区間の積分に移る。 ある区間の積分が $n = 100$ になっても必要精度に達しなかったら、メッセージを出して次の区間に移る。 全積分値は各区間の積分値の和とする。 この手順の PAD を書け。

［解］

図 5.8 の PAD を見よ。 この中で、分点 axis(*) と重み weight(*) は分点数 $n$ を増やしたときに1回だけ計算する。 なお、計算するのは分点 $t_i \geqq 0$ だけでよい。 各区間での被積分関数はほぼ同じ性質をもつとする。 全積分範囲で同じ性質の被積分関数のときは、 $N = 1$ として分点数 $n$ だけを大きくしていく方が収束は早い。

［例 5.8］

[例 5.6] の積分を 図 5.8 の PAD の手順で行う。 まず積分区間を $N = 8$ 等分しそれぞれを積分した時、収束したときの分点数 $n$ と、そのときまでの分点数の合計 $n(n + 1)/2$ と積分値 dS 、および積分値の和 S は次のようになった。 ただし $\varepsilon = 10^{-15}$ とした。

\[\begin{array}{rcl} n & n(n + 1)/2 & \qquad \texttt{dS} \cr 6 & 21 & \qquad 0.49741997818704614 \cr 8 & 36 & \qquad 0.48249467432040998 \cr 8 & 36 & \qquad 0.45516802857483140 \cr 6 & 21 & \qquad 0.41950775492093594 \cr 6 & 21 & \qquad 0.37980682537102555 \cr 6 & 21 & \qquad 0.33960717379888805 \cr 6 & 21 & \qquad 0.30131556331336067 \cr 6 & 21 & \qquad 0.26627265510329539 \cr \text{合計} & 198 & \texttt{S} = 3.141592653589793 \cr \end{array}\]

一方、 $N = 1$ として全区間を一挙に積分すると、次のように分点数は約半分ですむ。

\[\begin{matrix} n & n(n + 1)/2 & \texttt{S} \cr 14 & 105 & 3.141592653589795 \cr \end{matrix}\]

例題 5.9

例題 5.8 のように、ガウス・ラゲルの積分公式を用いて、分点数を必要 精度が得られるまで増加させる手順を考えて、PAD を書け。

［解］

ラゲルの多項式の零点が必要になるが、零点を求めるときのニュートン法の初期値はルジャンドルの多項式のように簡単でよい初期値は見あたらない。 そこで、DKA 法により許容相対誤差 $\varepsilon = 10^{-3}$ として複素近似値を求め、その実数部をニュートン法の初期値とする。 また、ラゲルの多項式の値はホーナー法によらず漸化式 $(5.102)$ を用いるから、DKA 法の解の存在範囲は、多項式の係数から求めるわけには行かない。 そこで $n$ 次のラゲルの多項式の零点は、

\[0 \leqq x_i \leqq 4n - 3 \qquad (1 \leqq i \leqq n) \tag{5.126}\]

の範囲にあることを利用する。 積分値の許容相対誤差は $\varepsilon = 10^{-15}$ とする。 このように考えて書いた PAD は 図 5.9 のようになる。

［例 5.9］

図 5.9 の PAD にしたがって、

\[\int_0^\infty e^{-x} \frac{x^m}{m!} = 1 \tag{5.127}\]

を求める。 収束値に達したときの分点数 $n$ (実際の計算は収束確認のためもう1分点余計に行った) とそのときの積分値は以下のようになった。

\[\begin{array}{ccccrcc} m & n & \text{収束値} & & m & n & \text{収束値} \cr 1 & 1 & 1.0000000000000000 & & 6 & 4 & 0.9999999999999998 \cr 2 & 2 & 1.0000000000000000 & & 7 & 4 & 0.9999999999999998 \cr 3 & 2 & 0.9999999999999998 & & 8 & 5 & 1.0000000000000000 \cr 4 & 3 & 1.0000000000000000 & & 9 & 5 & 1.0000000000000000 \cr 5 & 3 & 0.9999999999999998 & & 10 & 6 & 0.9999999999999995 \cr \end{array}\]

例題 5.10

例題 5.9 のように、ガウス・エルミートの積分公式を用いて、分点数を必要精度が得られるまで増加させる手順を考えて、PAD を書け。

［解］

エルミートの多項式の零点の存在範囲は、

\[0 \leqq x_i^2 \leqq 4n + 3 \qquad (1 \leqq i \leqq n) \tag{5.128}\]

である。 そのときの PAD を 図 5.10 に示す。 なお、 $n$ が奇数のときは $0$ も零点であるから、収束判定には許容絶対誤差も必要である。

［例 5.10］

図 5.10 の PAD にしたがって、

\[\int_{-\infty}^\infty e^{-x^2} \frac{2^m x^{2m}}{(2m - 1)!!} \ dx = \sqrt{\pi} = 1.772453850905516027 \cdots \tag{5.129}\]

を求める。 ただし、 $(2m-1)!! = (2m-1)(2m-3) \cdots 3 \cdot 1$ 。 結果は次の通りであった。

\[\begin{matrix} m & n & \text{収束値} \cr 1 & 2 & 1.77245385090551 \underline{54} \cr 2 & 3 & 1.77245385090551 \underline{57} \cr 3 & 4 & 1.772453850905516 \underline{2} \cr 4 & 5 & 1.77245385090551 \underline{59} \cr 5 & 6 & 1.77245385090551 \underline{71} \cr \end{matrix}\]

---

5.6 2重指数関数型公式

5.6.1 DE 公式

$(5.8)$ の台形公式 $I_1$ はもっとも簡単な数値積分公式であるが、無限区間の積分

\[I = \int_{-\infty}^\infty e^{-x^2} f(x) \ dx\]

に適用すると、精度の非常によい結果が得られる。 例えば、 $f(x) = 1$ の場合 $I = \sqrt{\pi}$ であるが、台形公式 $I_1$ を適用すると、

\[I_1 = h \sum_{j=-\infty}^\infty e^{-(jh)^2}\]

である。 この公式による絶対誤差は $\vert I - I_1 \vert \cong 2 \sqrt{\pi} exp(-\pi / h^2)$ であり、 $h$ が小さいときは精度はきわめてよい。

高精度数値積分公式として、最近注目されている数値積分公式に 2重指数関数型公式 (double exponential formula)、あるいは簡単に $\text{DE}$ 公式 と呼ばれる公式がある。 この公式は、有限区間の積分を無限区間の積分に変換し、かつ被積分関数を2重指数関数的に小さくなる関数に変換し、台形公式を適用することによって得られる公式である。

$(5.1)$ の積分を、ガウスの積分公式と同じく、

\[x = \frac{1}{2} (b - a) s + \frac{1}{2} (a + b) \tag{5.130}\]

と変数変換して、積分範囲を $[a, b]$ から $[-1, 1]$ に変え、

\[I = \int_a^b f(x) ~dx = \frac{1}{2} (b - a) \int_{-1}^1 f(x(s)) ~ds \tag{5.131}\]

と置く。 さらに積分範囲を $[-1, 1]$ から $(-\infty, \infty)$ に変えるように、

\[s = \tanh( \sinh \xi ) \tag{5.132}\]

と変数変換すると^$2)$

$^{2)}$ 通常の教科書では、 $s = \tanh( (\pi / 2) \sinh \xi )$ としてあるが、係数 $\pi / 2$ は無意味であるので省いた。

\[I = \frac{1}{2} (b - a) \int_{-\infty}^{\infty} f(x(s(\xi))) w(\xi) ~d\xi \tag{5.133}\]

が得られる。 ただし、

\[w(\xi) = \frac{\cosh \xi}{\cosh^2 [ \sinh \xi ]} \tag{5.134}\]

である。 ここに用いた双曲線関数 (hyperbolic function) は指数関数によって

\[\begin{aligned} \sinh x & = \frac{e^x - e^{-x}}{2} \cr \cosh x & = \frac{e^x + e^{-x}}{2} \cr \tanh x & = \frac{\sinh x}{\cosh x} = \frac{e^x - e^{-x}}{e^x + e^{-x}} \cr \end{aligned} \tag{5.135}\]

と定義された関数で、次の性質を持つ。

\[\begin{aligned} & \cosh^2 x - \sinh^2 x = 1 \cr & \sinh(x + y) = \sinh x \cosh y + \cosh x \sinh y \cr & \cosh(x + y) = \cosh x \cosh y + \sinh x \sinh y \cr \end{aligned} \tag{5.136}\] \[\begin{aligned} \frac{d}{dx} \sinh x & = \cosh x \cr \frac{d}{dx} \cosh x & = \sinh x \cr \frac{d}{dx} \tanh x & = \frac{1}{\cosh^2 x} \cr \end{aligned} \tag{5.137}\]

積分 $(5.133)$ の上限と下限は無限大であるが、被積分関数の重み $w(\xi)$ は $\xi \to \pm \infty$ で

\[w(\xi) = \frac{\cosh \xi}{\cosh^2 [ \sinh \xi ]} \cong 2 \exp[ - \exp( \vert \xi \vert ) ] \to 0\]

すなわち、2重指数関数的に $0$ に近づく。

$w(\xi)$ は $\vert \xi \vert$ が大きいときアンダーフローするが、このとき $w(\xi) = 0$ とおいてもあまり問題はない。 しかし、 $w(\xi)$ の算出のときに $w(\xi)$ の分母は $\vert \xi \vert$ が大きくなると2重指数関数的にオーバーフローする。 オーバーフローすれば計算機は通常実行を停止してしまう。 オーバーフローを避けるために、 $F$ をオーバーフローしない程度の大きな数とすると、 $\vert \xi \vert$ の最大値 $H$ は

\[\cosh^2 [\sinh H] < F \tag{5.138}\]

を満たさなくてはならない。 すなわち、 $H$ は

\[H < \sinh^{-1} [ \cosh^{-1} \sqrt{F} ] \tag{5.139}\]

の右辺より小さい必要がある。 $F$ の値は、指数部7ビット仮数部56ビットの16進計算機 (IBM 方式) では $F = 16^{63} (1 - 2^{-56}) = 7.237 \times 10^{75}$ 、指数部11ビット仮数部52ビットの2進計算機 (IEEE 方式) では $F = 2^{1023} (2 - 2^{-52}) = 1.798 \times 10^{308}$ であるから、

\[H < 5.171 \text{ (IBM 方式),} \quad 6.567 \text{ (IEEE 方式)}\]

である。

分母のオーバーフローを避けるもっとも簡単な方法は、 $(5.134)$ の $w(\xi)$ の分母と分子に2重指数関数的に減少する関数を掛けてやればよい。 そのために分子がアンダーフローするかも知れないが、そのときは $w(\xi) = 0$ と置けばよい (これは計算機が自動的にそうしてくれる)。 ここでは簡単に

\[w(\xi) = \frac{\cosh \xi \exp[ -2 \sinh \xi ]}{\{\cosh[ \sinh \xi ] \exp[ - \sinh \xi ]\}^2} \tag{5.140}\]

とする。 $\xi$ が大きいときもっとも早くオーバーフローが起こる可能性があるのは、分母の $\cosh[ \sinh \xi ]$ の演算のときであるが、これがオーバーフローしないための条件は、

\[H < \sinh^{-1} [ \cosh^{-1} F ] \tag{5.141}\]

この値は

\[H < 5.860 \text{ (IBM 方式),} \quad 7.259 \text{ (IEEE 方式)}\]

となり、 $H$ は若干大きくとれる。

$(5.132)$ より $s(\xi)$ は $\xi \to \infty$ のとき $s(\xi) \to 1$ 、また $\xi \to -\infty$ のとき $s(\xi) \to -1$ となるから、 $(5.130)$ を書変えた

\[x = \frac{1}{2} a (1 - s) + \frac{1}{2} b (1 + s)\]

より、 $1 \pm s(\xi)$ は $\vert \xi \vert \gg 1$ のとき桁落ちが起こりやすい。 被積分関数 $f(x)$ が上限や下限で特異性を持つときには、この桁落ちは致命的な誤差を産み出す。 桁落ちを避けるためには、 $(5.130)$ に $(5.132)$ の $s(\xi)$ を代入して得られる

\[x = \frac{a \exp(- \sinh \xi) + b \exp(\sinh \xi)}{2 \cosh(\sinh \xi)} \tag{5.142}\]

によって $x$ を計算することが必要である。 なお、特異関数の積分については、後の 5.6.3 節に詳述する。

$(5.133)$ に台形公式 $(5.8)$ を適用して得られる公式は

\[\boxed{ \begin{align*} \quad & & \quad \cr & \text{DE 公式 (2重指数関数型公式)} \cr & \quad I_1 = \frac{1}{2} (b - a) h \sum_{j=-N}^N w_j f(x_j) \tag{5.143} \cr & \qquad x_j = \frac{a \exp ( - \sinh jh ) + b \exp ( \sinh jh )}{2 \cosh ( \sinh jh )} \tag{5.144} \cr & \qquad w_j = \frac{\cosh jh \exp(-2 \sinh jh)}{[ \cosh(\sinh jh) \exp(- \sinh jh) ]^2} = w_{-j} \tag{5.145} \cr \end{align*} }\]

である。 $N$ は理論計算上は無限大であるが、上に述べた数値計算上の理由から、 $N h < H$ によって決まる整数である。

なお、 $w_j$ は $\vert jh \vert$ が大きくなると2重指数関数的に小さくなるから、 $(5.143)$ の和は $\vert j \vert$ の大きな項から加え合わせねばならない。

---

例題 5.11

DE 公式を利用して、例題 5.1 と同じように、必要精度の積分値を得るまで、台形公式のきざみを半分割して行く手順を求めて PAD を書け。

［解］

最大半分割回数を MXHLF とする。 このときのの区間数は $2 \cdots 2^\texttt{MXHLF}$ である。

\[I_1 = \frac{1}{2} (b - a) S \tag{5.146}\]

とおき、分点が ${-H, 0, H}$ の3点のときの $S$ を $S^0$ とすると、 $h = H$ として

\[S^0 = H \left[ f \left( \frac{1}{2} (a + b) \right) w_0 + [ f(x_{-1}) + f(x_1) ] w_1 \right] \tag{5.147}\]

ただし、 $w_0 = 1$ である。 一般に $k$ 回目の半分割のときの区間幅は $h = H / 2^k$ であり、このときの $S$ を $S_k$ とすると、例題 5.1 と同じように、

\[S^k = \frac{1}{2} S^{k-1} + h \sum_{j=\text{奇数}}^{2^k - 1} [ f(x_j) + f(x_{-j}) ]w_j \tag{5.148}\]

$H = 5.860$ とした手順の PAD を図 5.11 に示す。

［例 5.11］

次の積分を図 5.11 の PAD にしたがって計算したときの結果を示す。 許容相対誤差は $\varepsilon = 10^{-16}$ とした。

\[I = \int_{-1}^1 \frac{2}{1 + x^2} ~dx = \pi\] \[\begin{matrix} \text{半分割回数} & \cr 0 & \underline{11.72000000000000} \cr 1 & \underline{5.860001708167349} \cr 2 & 3.\underline{374160156023132} \cr 3 & 3.14\underline{6962440347332} \cr 4 & 3.14159\underline{4991730010} \cr 5 & 3.1415926535\underline{90228} \cr 6 & 3.141592653589793 \cr 7 & 3.141592653589793 \cr \end{matrix}\]

このように2次収束していることがわかる。 有効数字16桁目まで正しい。 半分割回数7回目は収束確認のための計算である。

参考のため述べると、同じ問題を $(5.8)$ のままの台形公式を使ったとき (例題 5.1 参照) は、半分割回数8回では6桁しか合わない。 誤差は $h^2$ に比例するから、各回ごとに約1/4の割合で小さくなって行く。 この割合で行くと15桁正しくなるためには、丸めの誤差はないとしても、半分割回数25回位必要である。

一方、ロンバーグ積分法 (補外法) によると、半分割回数7回で15桁正しく収束する ([例 5.6 ] 参照)。

5.6.2 DE 公式の収束性

$h$ が小さいとき、誤差 $I_1 - I$ はほぼ $\exp(-c / h)$ に比例することがわかっている。 したがって、半分割すれば $\exp(-2c / h)$ に比例することになり、誤差は2乗倍される。 すなわち2次収束である。

この公式は、変数変換 $(5.132)$ によって有限区間を無限区間に引き延ばし、被積分関数 $f(x)$ の複雑な振る舞いをゆっくりした振る舞いに変換して、台形公式を適用したと見ることが出来る。 引き延ばしの倍率 $d \xi / dx$ は、 $(5.130)$ および $(5.132)$ より

\[\frac{dx}{d\xi} = \frac{1}{2} (b - a) w(\xi) \tag{5.149}\]

に反比例するから、上下限に近い点が大きい。 したがって、上下限に特異性がある関数に対しては、変換なしの単純な台形公式にくらべて、精度は圧倒的によい。

逆に、単純な $f(x)$ に対しては変換なしの台形公式の方がよい。 たとえば、 $f(x)$ が1次式のときは、単純な台形公式は正確な解を区間半分割なしに与える。 一方、 DE 公式はかなりの分割数が必要である。 これは極端な例であるが、この公式は複雑な振る舞い、特に次に述べるように、特異点を持つ被積分関数のときに有効な公式である。

5.6.3 特異関数の数値積分

特異関数の積分は、先ず特異性を除くように変数変換をしてから、数値積分の公式を適用すべきである。 例えば、次の積分は変数変換 $x = \cos \theta$ によって、特異性の除去ができる。 (ここでの $f(x)$ は特異性を持たないとする。)

\[I = \int_{-1}^{1} \frac{f(x)}{\sqrt{1 - x^2}} ~dx = \int_{0}^{\pi} f(\cos \theta)d\theta\]

除去が難しい場合には、特異点を上下限に移動させる。 DE 公式の重み $w_j$ は上下限近傍では小さく、かつ、積分変数 $x$ は上下限に近づくが上下限に達することはないから精度は非常によい。 もし、上下限にない特異点があるときは、上下限のみに特異点があるいくつかの積分に分けて積分する。 以下、特異点は上下限にのみあるとする。

いま、 $(5.1)$ において、下限 $a$ のみが $f(x)$ の特異点とする。 このとき、特異点の近傍では $x - a$ の計算で桁落ちが起こるから、変数変換 $y = x - a$ を行って

\[I = \int_{a}^{b} f(x) ~dx = \int_{0}^{b-a} f(a + y) ~dy = \int_{0}^{b-a} \bar{f}(y) ~dy \tag{5.150}\]

とする。 ただし $\bar{f} (y) = f(a + y)$ である。 $\bar{f}(y)$ の特異点は $y = 0$ である。

上限 $b$ のみに特異点があるときには、 $(5.1)$ に変数変換 $y = b - x$ を行って

\[I = \int_{a}^{b} f(x) ~dx = \int_{0}^{b-a} f(b - y) ~dy = \int_{0}^{b-a} \bar{f} (y) ~dy \tag{5.151}\]

ここに $\bar{f} (y) = f(b - y)$ である。 $\bar{f} (y)$ の特異点は $y = 0$ である。

下限 $a$ と上限 $b$ の両方に特異点が存在し、その他には存在しないときには、

\[\begin{aligned} I & = \int_{a}^{b} f(x) ~dx = \int_{a}^{(a+b)/2} f(x) ~dx + \int_{(a+b)/2}^{b} f(x) ~dx \cr & = \int_{0}^{(b-a)/2} [ f(a + y) + f(b - y) ] ~dy = \int_{0}^{(b-a)/2} \bar{f} (y) ~dy \cr \end{aligned} \tag{5.152}\]

ただし $\bar{f} (y) = f(a + y) + f(b - y)$ である。 $\bar{f} (y)$ の特異点は $y = 0$ である。

こうして、いずれの場合も特異点は $y = 0$ のみにある被積分関数 $\bar{f} (y)$ とすることができる。 これに DE 公式を適用すれば

\[I_1 = \frac{1}{2} (b - a) h \sum_{j=-N}^N w_j \bar{f} (y_j) \tag{5.153}\]

ここで大切なことは、特異点 $y = 0$ の近傍で $\vert y \vert < \varepsilon_M \ll 1$ のとき、 $1 \pm y = 1$ となり $y \neq 0$ は丸めの誤差のため無視されてしまい、そのため $\bar{f} (y)$ は無限大 $\bar{f}(0)$ になってしまうことである。 次の例でこのことを注意しよう。

［例 5.12］

次の積分

\[I = \int_{-1}^{1} \frac{dx}{\sqrt{1 - x^2}} = \pi\]

を考える。 特異点は上限 $+1$ と下限 $-1$ の両方にある。 したがって、 $a = -1$ 、 $b = 1$ として $(5.152)$ を適用する。 被積分関数は

\[\bar{f} (y) = f(a + y) + f(b - y) = \frac{2}{\sqrt{1 - (1 - y)^2}}\]

しかし、特異点のごく近傍 $y < \varepsilon_M \ll 1$ では $1 - y = 1$ となり、次に $1 -(1 - y)^2 = 0$ となり、 $\bar{f} (y)$ は無限大 $\bar{f} (0)$ となる。 ここでは $1 - (1 - y)^2 = y(2 - y)$ としなければならない。 こうすれば、 $\vert j \vert$ は大きくても有限であるから $y_j = 0$ となることはなく、したがって、 $y_j (2 - y_j) = 0$ となることはない。 この配慮をした後に DE 公式を適用すると

\[I = \int_{-1}^{1} \frac{dx}{\sqrt{1 - x^2}} = \int_{0}^{1} \frac{2~dy}{\sqrt{y(2 - y)}} \cong h \sum_{j=-N}^N \frac{w_j}{\sqrt{y_j (2 - y_j)}} \tag{5.154}\]

が得られる。 $y_j$ は $(5.144)$ において $a = 0$ 、 $b = 1$ として

\[y_j = \frac{\exp(\sinh j_h)}{2 \cosh(\sinh j_h)}\]

である。

図 5.11 の PAD による結果を次に示す。

\[\begin{matrix} \text{半分割回数} \cr 0 & \underline{6.766545154902415} \cr 1 & 3.\underline{390129003450192} \cr 2 & 3.1\underline{03156970830038} \cr 3 & 3.1415\underline{71845776584} \cr 4 & 3.14159265\underline{2854071} \cr 5 & 3.141592653589793 \cr 6 & 3.141592653589793 \cr \end{matrix}\]

すなわち、[例 5.11] と同じように2次収束する。

上に述べた桁落ち防止処理をしないで DE 公式を使うとき、特異点近傍での無限大をさけるには、 $(5.141)$ の $H$ を小さくとって $H = 1$ とすることも考えられる。 このときには特異点における無限大は避けられるが、しかし、 $H$ を小さくとったことにより積分値の精度は悪くなり、半分割を進めても振動して収束しない。

なお、 DE 公式以外の公式では被積分関数が特異点で無限大となり計算不可能である。

5.7 第5章の問題

問題 5-1

ファンデルモンドの行列式

\[\det A = \det \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \cr 1 & x_1 & x_1^2 & \cdots & x_1^n \cr 1 & x_2 & x_2^2 & \cdots & x_2^n \cr \vdots & \vdots & \vdots & \ddots & \vdots \cr 1 & x_n & x_n^2 & \cdots & x_n^n \cr \end{pmatrix} = \prod_{i>j} (x_i - x_j) \tag{5.25}\]

を証明せよ。

［解］

この行列式は $x_0, x_1, \cdots x_n$ の多項式である。 $x_i = x_j$ とすると、行列 $A$ の第 $i$ 行と第 $j$ 行が等しくなり行列式は $0$ となるから、行列式は因子 $x_i - x_j$ を含む (剰余定理)。 したがって、

\[\begin{aligned} \det A = \alpha & (x_1 - x_0) \cr \times & (x_2 - x_0)(x_2 - x_1) \cr \times & (x_3 - x_0)(x_3 - x_1)(x_3 - x_2) \cr & \cdots \quad \cdots \qquad \cdots \cr \times & (x_n - x_0)(x_n - x_1)(x_n - x_2) \cdots (x_n - x_n-1) \cr \end{aligned}\]

と書ける。 この両辺の次数は $n(n + 1) / 2$ であるから、 $\alpha$ は $x_i$ を含まない定数である。 右辺の各因子の第1項の積は $\alpha x_1 x_2^2 x_3^3 \cdots x_n^n$ であり、これは左辺の行列の対角要素の積であるから、係数は $+1$ である。 したがって $\alpha = 1$ である。

別解として、差分商を使う方法を示す。 まず、 $f_k (x) = x^k$ として $n - k + 1$ 次行列 $A_k$ を差分商を使って

\[A_k = \begin{pmatrix} 1 & f_{k+1} [x_k, \quad x_{k-1}, \cdots, x_0] & \cdots & f_n [x_k, \quad x_{k-1}, \cdots, x_0] \cr 1 & f_{k+1}[x_{k+1}, x_{k-1}, \cdots, x_0] & \cdots & f_n[x_{k+1}, x_{k-1}, \cdots, x_0] \cr \vdots & \vdots & \ddots & \vdots \cr 1 & f_{k+1} [x_n, \quad x_{k-1}, \cdots, x_0] & \cdots & f_n [x_n, \quad x_{k-1}, \cdots, x_0] \cr \end{pmatrix}\]

と定義する。 このとき、 $A_0 = A$ 、 $A_n = (1)$ である。 差分商の関係

\[\begin{aligned} f_l [x_j, x_{k-1}, \cdots, x_0] & - f_l [x_k, x_{k-1}, \cdots, x_0] \cr & = (x_j - x_k) fl[x_j, x_k, x_{k-1}, \cdots, x_0] \cr \end{aligned}\]

を考慮して $\det A_k$ の第1行を第2行以下から引くと

\[\begin{aligned} \det A_k =~ & (x_{k+1} - x_k)(x_{k+2} - x_k) \cdots (x_n - x_k) \cr & \times \det \begin{pmatrix} 1 & f_{k+1} [x_k, x_{k-1}, \cdots, x_0] & \cdots & f_n [x_k, x_{k-1}, \cdots, x_0] \cr 0 & f_{k+1} [x_{k+1}, x_k, \cdots, x_0] & \cdots & f_n [x_{k+1}, x_k, \cdots, x_0] \cr \vdots & \vdots & \ddots & \vdots \cr 0 & f_{k+1} [x_n, x_k, \cdots, x_0] & \cdots & f_n [x_n, x_k, \cdots, x_0] \cr \end{pmatrix} \end{aligned}\]

行列式を第1列で展開し、 $f_{k+1} [x_j, x_k, \cdots, x_0] = 1$ であることから

\[\begin{aligned} =~ & (x_{k+1} - x_k)(x_{k+2} - x_k) \cdots (x_n - x_k) \cr & \times \det \begin{pmatrix} 1 & f_{k+2} [x_{k+1}, x_k, \cdots, x_0] & \cdots & f_n [x_{k+1}, x_k, \cdots, x_0] \cr 1 & f_{k+2} [x_{k+2}, x_k, \cdots, x_0] & \cdots & f_n [x_{k+2}, x_k, \cdots, x_0] \cr \vdots & \vdots & \ddots & \vdots \cr 1 & f_{k+2} [x_n, x_k, \cdots, x_0] & \cdots & f_n [x_n, x_k, \cdots, x_0] \cr \end{pmatrix} \cr =~ & (x_{k+1} - x_k)(x_{k+2} - x_k) \cdots (x_n - x_k) \cdot \det A_{k+1} \cr =~ & \prod_{j=k+1}^{n} (x_j - x_k) \cdot \det A_{k+1} \end{aligned}\]

この関係式を $k = 0, 1, \cdots n$ に適用すれば、

\[\begin{aligned} \det A = & \prod_{j=1}^n (x_j - x_0) \cdot \det A_1 = \prod_{j=1}^n (x_j - x_0) \cdot \prod_{j=2}^n (x_j - x_1) \cdot \det A_2 \cr & \cdots \qquad \cdots \cr = & \prod_{k=0}^{n-1} \prod_{j=k+1}^n (x_j - x_k) \cdot \det A_n = \prod_{j>k} (x_j - x_k) \cr \end{aligned}\]

---

問題 5-2

いま $\phi_n (x) = (x - x_1)(x - x_2) \cdots (x - x_n)$ とする。 このとき $f(x) = {\phi_n (x)}^2$ の $n$ 階差分商は

\[f[x, x_1, x_2, \cdots x_n] = \phi_n (x) \tag{5.155}\]

［解］

差分商の性質 $(5.35)$ を利用する。 $\phi_{n+1} (y) = (y-x) \phi_n (y)$ と置くと、 $\phi_{n+1}^\prime (y) = \phi_n (y) + (y - x) \phi_n^\prime (y)$ 。 したがって

\[\phi_{n+1}^\prime (x) = \phi_n (x) \qquad \text{および} \qquad \phi_{n+1}^\prime (x_i) = (x_i - x) \phi_n^\prime (x_i)\]

$f(x_i) = 0 (1 \leqq i \leqq n)$ に注意して

\[\begin{aligned} f[x, x_1, x_2, \cdots , x_n] & = \frac{f(x)}{\phi_{n+1}^\prime (x)} + \sum_{i=1}^n \frac{f(x_i)}{\phi_{n+1}^\prime (x_i)} \cr & = \frac{f(x)}{\phi_n (x)} = \phi_n (x) \end{aligned}\]

---

問題 5-3

$(5.46)$ は $x_0, x_1, x_2, \cdots , x_n$ が相異なるとき $(5.35)$ に一致し、すべて同じとき $(5.43)$ に一致することを示せ。

［解］

$u_k = x_0 + \sum_{l=1}^k t_l (x_l - x_{l-1})$ と置く。 $n = 1$ のとき $(5.46)$ は

\[\begin{aligned} \int_{0}^{1} & f^\prime (x_0 + t_1 (x_1 - x_0)) ~dt_1 \cr =~ & \frac{1}{x_1 - x_0} \int_{x_0}^{x_1} f^\prime (u_1) ~du_1 = \frac{f[x_1] - f[x_0]}{x_1 - x_0} = f[x_1, x_0] \end{aligned}\]

であるから成立している。 $n = k$ のとき成立していると仮定する。 $n = k + 1$ のときを考えると

\[\begin{aligned} \int_{0}^{t_k} f^{(k+1)} (u_{k+1}) dt_{k+1} & = \int_{u_k}^{u_{k-1} + t_k (x_{k+1} - x_{k-1})} \frac{f^{(k+1)} (u_{k+1})}{x_{k+1} - x_k} ~du_{k+1} \cr & = \frac{f^{(k)} (u_{k-1} + t_k (x_{k+1} - x_{k-1})) - f^{(k)} (u_k)}{x_{k+1} - x_k} \cr \end{aligned}\]

これを $n = k + 1$ のときの $(5.46)$ の右辺に代入すると

\[\begin{aligned} \int_{0}^{1} & \int_0^{t_1} \cdots \int_{0}^{t_{k-1}} \int_0^{t_k} f^{(k+1)}(u_{k+1}) ~dt_{k+1} ~dt_k ~\cdots ~dt_1 \cr =~ & \int_{0}^{1} \int_{0}^{t_1} \cdots \int_{0}^{t_{k-1}} \frac{f^{(k)}(u_{k-1} + t_k (x_{k+1} - x_{k-1})) - f^{(k)}(u_k)}{x_{k+1} - x_k} ~dt_k ~dt_{k-1} ~ \cdots ~dt_1 \cr =~ & \frac{f[x_{k+1}, x_{k-1}, \cdots , x_1, x_0] - f[x_k , x_{k-1}, \cdots , x_1, x_0]}{x_{k+1} - x_k} \cr =~ & f[x_{k+1}, x_k, x_{k-1}, \cdots , x_1, x_0] \cr \end{aligned}\]

であるから、 $(5.46)$ は $n = k + 1$ のときにも成立する。 故に、 $(5.46)$ は任意の $n$ に対して成立する。

$x_n = x_{n-1} = \cdots = x_1 = x_0 = x$ のときには、 $(5.46)$ より

\[\begin{aligned} & \int_{0}^{1} \int_{0}^{t_1} \int_{0}^{t_2} \cdots \int_{0}^{t_{n-1}} f^{(n)}(x_0)~dt_n ~dt_{n-1} ~\cdots ~dt_2 ~dt_1 \cr & = f^{(n)}(x_0) \int_{0}^{1} \int_{0}^{t_1} \int_{0}^{t_2} \cdots \int_{0}^{t_{n-1}} ~dt_n ~dt_{n-1} ~\cdots ~dt_2 ~dt_1 \cr & = \frac{1}{n!} f^{(n)}(x_0) \end{aligned}\]

---

問題 5-4

$x = x_k, x_{k+1}, \cdots , x_{k+n-1}, x_{k+n}$ で $f(x)$ と一致する $n$ 次多項式 $p_{n}(x)$ をあらためて $p_{nk}(x)$ と書くと、 $(5.47)$ は、

\[\boxed{ \begin{aligned} \quad & & \quad \cr & \text{ネビユの補間公式} \cr & \quad p_{nk}(x) = \frac{(x_{k+n} - x)p_{n-1~k}(x) + (x - x_k)p_{n-1~k}+1(x)}{x_{k+n} - x_k} \tag{5.156} \end{aligned} }\]

と書けることを示せ。 この漸化式は ネビユ (Neville) の補間公式 または ネビユの算法 と呼ばれる。 この公式による数値補間の手順を考えよ。

［解］

$(5.47)$ の第2式の第1項は $p_{n-1~k}(x)$ 、第2項は

\[\begin{aligned} \phi_n (x) f & [x_{k+n}, x_{k+n-1}, \cdots , x_{k+1}, x_k] \cr =~ & \frac{x - x_k}{x_{k+n} - x_k} \phi_{n-1} (x) {f[x_{k+1}, x_{k+2}, \cdots , x_{k+n}] - f[x_k , x_{k+1}, \cdots , x_{k+n-1}]} \cr =~ & \frac{x - x_k}{x_{k+n} - x_k} \{ p_{n-1~k}+1(x) - p_{n-1~k}(x) \} \end{aligned}\]

と書ける。 これらを加えれば、 $(5.156)$ を得る。

$x_k \leqq x < x_{k+1}$ とする。 $p_{0j} (x) = f(x_j)$ $(j = k - 1, k, k + 1)$ から2つの差分商 $f[x_k , x_{k+1}]$ と $f[x_{k-1}, x_k]$ を求め、 $(5.156)$ により、

\[\begin{array}{ccl} \vert f[x_{k-1}, x_k] \vert < \vert f[x_k, x_{k+1}] \vert & \text{のとき} & p_{1~k-1}(x) \cr \vert f[x_{k-1}, x_k] \vert \geqq \vert f[x_k, x_{k+1}] \vert & \text{のとき} & p_{1k}(x) \cr \end{array}\]

を求める。 すなわち、分点は $p(x)$ が区分的に滑らかに変化する方向に付加する。 同様に、次々と分点をつけ加えつつ補間多項式の精度を上げていく。 もし差分商があまりに大きいときは補間を中止する。

---

問題 5-5

ニュートン・コーツの閉型公式 $(5.69)$ の $w_i$ および $\gamma_m$ (表 5-1 の「 $E_n$ の係数」) の間には、対称性 $w_i = w_{n-i}$ 、および $n + 2$ 個の関係

\[\begin{alignedat}{3} & \sum_{i=0}^n i^k w_i & & = \frac{n^{k+1}}{k + 1} & \qquad & 0 \leqq k \leqq n \cr & \sum_{i=0}^n i^m w_i + m!~ \gamma_m & & = \frac{n^{m+1}}{m + 1} & \qquad & m = \begin{cases} n + 2 & n: \text{偶数} \cr n + 1 & n: \text{奇数} \cr \end{cases} \end{alignedat} \tag{5.157}\]

が成り立つことを示せ。 関係 $(5.157)$ は、 $w_i$ および $\gamma_m$ を積分でなく連立1次方程式の解として求めるのに利用できる^$3)$ 。

$^{3)}$ 記号処理言語 Mathematica を使って有理数で求めるか、 FORTRAN や C のような言語を使って浮動小数点で求める。 浮動小数点で求めるときにはピボット選択なしのガウスの消去法で十分である。

［ヒント］

対称性は $(5.69)$ の積分から直接証明できる。 $(5.157)$ は、 $(5.69)$ 、 $(5.76)$ 、 $(5.77)$ より、 $f(x) = x^k$ $(0 \leqq k \leqq n)$ のとき $E_n = 0$ 、 $f(x) = x^m$ のとき $E_n = m!~\gamma_m / n^{m+1}$ であることに注意する。

---

問題 5-6

ニュートン・コーツの開型公式 $(5.81)$ について、問題 5-5 と同様な関係を導け。

［ヒント］

\[\begin{alignedat}{3} & \sum_{i=0}^n (i + 1)^k w_i && = \frac{(n + 2)^{k+1}}{k + 1} & \qquad & 0 \leqq k \leqq n \cr & \sum_{i=0}^n (i + 1)^m w_i + m!~\gamma_m && = \frac{(n + 2)^{m+1}}{m + 1} & \qquad & m = \begin{cases} n + 2 & n: \text{偶数} \cr n + 1 & n: \text{奇数} \cr \end{cases} \end{alignedat}\]

---

問題 5-7

ニュートン・コーツの公式の打切り誤差 $(5.74)$ と $(5.75)$ を $(5.73)$ から導け。

［解］

$n + 1$ 個の分点を $x_k = x_0 + kh(0 \leqq k \leqq n)$ とする。

\[\phi_{n+1} (x) = (x - x_0)(x - x_1) \cdots (x - x_{n-1})(x - x_n) \tag{5.158}\]

において、中点 $x_m = \frac{1}{2} (x_0 + x_n) = x_0 + \frac{1}{2} nh$ を原点とするような変数変換 $y = x - x_m$ をおこない、 $y_k = (k - \frac{1}{2} n) h$ と置くと、 $y_k = - y_{n-k}$ だから、

\[\begin{aligned} \phi_{n+1} (x) & = (y - y_0)(y - y_1) \cdots (y - y_{n-1})(y - y_n) \cr & = (y - y_0)(y - y_1) \cdots (y + y_1)(y + y_0) \cr & = \begin{cases} (y^2 - y_0^2)(y^2 - y_1^2) \cdots (y^2 - (h/2)^2) & ( n \text{ は奇数} ) \cr (y^2 - y_0^2)(y^2 - y_1^2) \cdots (y^2 - h^2)y & ( n \text{ は偶数} ) \cr \end{cases} \end{aligned}\]

故に、 $\phi_{n+1} (x)$ は中点 $x = x_m (y = 0)$ に関して、 $n$ が奇数ならば偶関数、偶数ならば奇関数である。

また、 $x$ が $x_0$ または $x_n$ から中点 $x_m$ に向かうにしたがって、 $\phi_{n+1} (x)$ は分点で $0$ になりながら振動し、振幅は減少する。 何故ならば、左半分 $( y < 0 )$ の分点でないところ $( y \neq y_k )$ では、 $h - y_k = - y_{k-1}$ なる関係を使って、

\[\begin{aligned} \left\vert \frac{\phi_{n+1} (x + h)}{\phi_{n+1} (x)}\right\vert & = \left\vert \frac{(y + h - y_0)(y - y_0) \cdots (y - y_{n-1})}{(y - y_0)(y - y_1) \cdots (y - y_n)} \right\vert \cr & = \frac{y + h - y_0}{y_n - y} \leqq \frac{\frac{1}{2} nh}{nh - ( \frac{1}{2} nh - h)} \cr & = \frac{n}{n + 2} < 1 \cr \end{aligned}\]

だからである。 右半分については、 $\phi_{n+1} (x)$ が $x_m$ に関して偶関数または奇関数であることから、明らかである。

次に、不定積分

\[\varPhi_{n+1} (x) = \int_{x_0}^{x} \phi_{n+1} (x) ~dx\]

を考える。 まず、 $n$ が奇数の時も偶数の時も、明らかに

\[\varPhi_{n+1} (x_0) = 0 \tag{5.159}\]

$n$ は偶数とする。

$\phi_{n+1} (x)$ は中点に対して奇関数だから、

\[\varPhi_{n+1} (x_n) = 0 \qquad (n \text{ は偶数}) \tag{5.160}\]

また、 $\varPhi_{n+1} (x)$ は中点に関して偶関数である。 何故なら、 $\phi_{n+1} (x)$ は中点に関して奇関数だから、

\[\begin{align*} \varPhi_{n+1} (x_m - x) & = \int_{x_0}^{x_{m-x}} \phi_{n+1} (x) ~dx \cr & = \int_{x_0}^{x_m - x} \phi_{n+1} (x) ~dx + \int_{x_m - x}^{x_m + x} \phi_{n+1} (x) ~dx \cr & = \int_{x_0}^{x_m + x} \phi_{n+1} (x) ~dx = \varPhi_{n+1} (x_m + x) \qquad (n \text{ は偶数}) \tag{5.161} \cr \end{align*}\]

次に、 $n$ が奇数のとき $\varPhi_{n+1} (x) > 0$ であることを示す。 まず $x_0 < x < x_1$ において $\phi_{n+1} (x) > 0$ であるから、 $\varPhi_{n+1} (x) > 0$ である。 $x_1 < x < x_2$ のとき $\phi_{n+1} (x) < 0$ であり、 $\varPhi_{n+1} (x)$ は減少するが、 $\phi_{n+1} (x)$ は中点に向かって減衰振動をしているから、 $x_0 < x < x_1$ のときの $\phi_{n+1} (x) > 0$ からの寄与が大きく、やはり $\varPhi_{n+1} (x) > 0$ である。 これを繰り返して、 $\varPhi_{n+1} (x)$ は増減を繰り返すが中点 $x = x_m$ まで $\varPhi_{n+1} (x) > 0$ を保つことがわかる。 $x > x_m$ のときは、 $\varPhi_{n+1} (x)$ が中点 $x_m$ に関して偶関数であることから、 $\varPhi_{n+1} (x) > 0$ である。 故に、 $x_0$ と $x_n$ を除く全区間について

\[\varPhi_{n+1} (x) > 0 \qquad x_0 < x < x_n \qquad (n \text{ は偶数}) \tag{5.162}\]

である。

以上の結果を使って、 $(5.73)$ を変形する。

$(1)$ $n$ が偶数のとき、 $(5.73)$ を部分積分して

\[\begin{aligned} \int_{x_0}^{x_n} & \phi_{n+1} (x) f [ x_0, x_1, \cdots , x_n, x ] ~dx \cr & = \Bigl[\varPhi_{n+1} (x) f [x_0, x_1, \cdots x_n, x] \Bigr]_{x_0}^{x_n} \cr & \quad - \int_{x_0}^{x_n} \varPhi_{n+1} (x) f [ x_0, x_1, \cdots , x_n, x, x ] ~dx \end{aligned}\]

この第1項は $0$ である。 また第2項において $\varPhi_{n+1} (x) > 0$ (定符号) だから、積分平均値の第2定理より、

\[\begin{aligned} & = - f [x_0, x_1, \cdots , x_n, \xi, \xi] \int_{x_0}^{x_n} \varPhi_{n+1} (x) ~dx \quad (x_0 < \xi < x_n) \cr & = - \frac{f^{(n+2)}(\xi)}{(n + 2)!} \int_{x_0}^{x_n} \varPhi_{n+1} (x) ~dx \quad (x_0 < \xi < x_n) \end{aligned}\]

部分積分により、

\[\begin{aligned} \int_{x_0}^{x_n} \varPhi_{n+1} (x) ~dx & = \Bigl[ x \varPhi_{n+1} (x) \Bigr]_{x_0}^{x_n} - \int_{x_0}^{x_n} x \phi_{n+1} (x) ~dx \cr & = - \int_{x_0}^{x_n} x \varPhi_{n+1} (x) ~dx \cr \end{aligned}\]

これより、 $(5.74)$ が得られた。

$(2)$ $n$ が奇数のときは、 $n - 1$ が偶数であることを利用する。 $(5.73)$ は差分商の定義より、

\[\begin{aligned} \int_{x_0}^{x_n} & \phi_{n+1} (x)f[x_0, x_1, \cdots , x_n, x] ~dx \cr & \quad = \int_{x_0}^{x_n} \phi_{n+1} (x) \frac{f[x_0, \cdots , x_{n-1}, x] - f[x_0, \cdots , x_n]}{x - x_n} ~dx \cr & \quad = \int_{x_0}^{x_n} \phi_n (x) \{ f[x_0, \cdots , x_{n-1}, x] - f[x_0, \cdots , x_n] \} ~dx \cr & \quad = \left( \int_{x_0}^{x_{n-1}} + \int_{x_{n-1}}^{x_n} \right) \phi_n (x) \{f[x_0, \cdots , x_{n-1}, x] - f[x_0, \cdots , x_n] \} ~dx \end{aligned}\]

$n - 1$ は偶数であるから、第1の積分は部分積分により、前項の $n$ を $n - 1$ で置き換えて、

\[\begin{aligned} \int_{x_0}^{x_{n-1}} & \phi_n (x) \{ f[x_0, \cdots x_{n-1}, x] - f[x_0, \cdots x_n]\} ~dx \cr & = - f [x_0, x_1, \cdots , x_{n-1}, \eta, \eta] \int_{x_0}^{x_{n-1}} \varPhi_n (x) ~dx \quad (x_0 < \eta < x_{n-1}) \cr \end{aligned}\]

また、部分積分により

\[\begin{aligned} \int_{x_0}^{x_{n-1}} \varPhi_n (x) ~dx & = - \int_{x_0}^{x_{n-1}} x \phi_n (x) ~dx = - \int_{x_0}^{x_{n-1}} (x - x_n) \phi_n (x) ~dx \cr & = - \int_{x_0}^{x_{n-1}} \phi_{n+1} (x) ~dx = - \varPhi (x_{n-1}) \end{aligned}\]

故に、第1の積分は

\[\begin{aligned} f[x_0, & x_1, \cdots , x_{n-1}, \xi, \xi]\varPhi(x_{n-1}) \cr & = \frac{f^{(n+1)}(\xi)}{(n + 1)!} \varPhi (x_{n-1}) \qquad (x_0 < \xi < x_{n-1}) \cr \end{aligned}\]

第2の積分は、積分範囲 $( x_{n-1} < x < x_n )$ で $\phi_{n+1} (x) < 0$ (定符号) だから、積分平均値の第2定理より、

\[\begin{aligned} \int_{x_{n-1}}^{x_n} & \phi_n (x) \{ f[x_0, \cdots , x_{n-1}, x] f [x_0, \cdots , x_{n-1}, x_n] \} ~dx \cr & = \int_{x_{n-1}}^{x_n} \phi_{n+1} (x)f[x_0, \cdots , x_n, x] ~dx \cr & = f[x_0, \cdots x_n, \xi^\prime] \int_{x_{n-1}}^{x_n} \phi_{n+1} (x) ~dx \quad (x_{n-1} < \xi^\prime < x_n) \cr & = \frac{f^{(n+1)} (\xi)}{(n + 1)!} \{ \varPhi(x_n) - \varPhi(x_{n-1}) \} \qquad (x_0 < \xi < x_n) \end{aligned}\]

こうして、 $(5.75)$ が得られる。

---

問題 5-8

台形公式 (複合公式) の打切り誤差が $(5.88)$ のように、 $h^2$ の級数で表されることを示せ。

［解］

$D$ を微分演算子とする： $D f(x) = (d / dx) f (x)$ 。 そして、 $D f(a) = D f(x) \vert_{x=a}$ と約束する。 $D$ を用いると、 $f(a + h)$ のテーラー展開は

\[\begin{align*} f(a + h) & = f(a) + hDf(a) + \frac{(hD)^2}{2!} f(a) + \cdots + \frac{(hD)^n}{n!} f(a) + \cdots \cr & = [ 1 + hD + \frac{(hD)^2}{2!} + \cdots + \frac{(hD)^n}{n!} + \cdots ] f(a) \cr & = e^{hD} f(a) \tag{5.163} \end{align*}\]

と書ける。 また、 $D^{-1}$ は微分の逆演算である積分とする： $D^{-1} f(x) = \int^x f(x)~dx$ 。 微分演算子 $D$ を用いると、台形公式 $(5.8)$ は

\[\begin{align*} I_1 & = \frac{h}{2} \left[ f(a) + 2 \sum_{n=1}^{N-1} f(a + nh) + f(b) \right] \cr & = \frac{h}{2} \left[ 1 + 2 \sum_{n=1}^{N-1} e^{nhD} + e^{NhD} \right] f(a) \cr & = \frac{h}{2} (e^{NhD} - 1) \frac{e^{hD/2} + e^{-hD/2}}{e^{hD/2} - e^{-hD/2}} f(a) \tag{5.164} \end{align*}\]

と書ける。 ここで

\[\coth x \equiv \frac{e^x + e^{-x}}{e^x - e^{-x}} = \frac{1}{x} + \frac{x}{3} - \frac{x^3}{45} + \frac{2x^5}{945} - \cdots = \sum_{k=0}^\infty \frac{2^{2k} B_{2k}}{(2k)!} x^{2k-1}\]

であることに注意する。 ただし、 $B_{2k}$ は ベルヌーイ数 (Bernoulli number) と呼ばれる定数で

\[\begin{aligned} & B_0 = 1, \quad B_1 = - \frac{1}{2}, \quad B_2 = \frac{1}{6}, \quad B_4 = - \frac{1}{30}, \quad \cdots \cr & B_{2k+1} = 0 \quad (k \geqq 1) \end{aligned}\]

などである。 また、演算子 $e^{NhD} - 1$ は、 $b = a + Nh$ に注意すれば、

\[(e^{NhD} - 1) f(a) = f(b) - f(a)\]

であるから、 $(5.164)$ は

\[\begin{align*} I_1 & = (e^{NhD} - 1) \left[ \frac{1}{D} + \sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} h^{2k} D^{2k-1} \right] f(a) \cr & = \int_{a}^{b} f(x)~dx + \sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} h^{2k} [ f^{(2k-1)} (b) - f^{(2k-1)} (a) ] \tag{5.165} \end{align*}\]

$(5.165)$ は、 オイラー・マクローリン (Euler-Maclaurin) の公式 と呼ばれる。 台形公式 (複合公式) の打切り誤差 $E_1 = I - I_1$ は $h^2$ の級数

\[E_1 = \sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} h^{2k} [ f^{(2k-1)} (a) - f^{(2k-1)} (b) ] \tag{5.166}\]

である。 $(5.166)$ の総和は、一般に $k \to \infty$ で収束しない。 このような級数を漸近級数と呼んでいる。 しかし、 $k$ を有限に限ったとき、 $h \to 0$ で $0$ に近づく。 特に $f(x)$ が周期関数ならば、すなわち、 $f^{(2k-1)} (a) = f^{(2k-1)} (b)$ が成り立つとき、台形公式は極めて精度がよい (打切り誤差は、 $\exp(-c/h)$ で減少する)。 このように、周期関数に対してはロンバーグ積分の補外は必要がないから、半分割 ( $h \to 0$ ) のみを行うがよい。

---

問題 5-9

一般にニュートン・コーツの閉型公式 (複合公式) について、 問題 5-8 と同様な議論から打切り誤差を求めよ。

［解］

準備として、 ベルヌーイの多項式 と呼ばれる多項式 $B_m (x)$ を次のように、 母関数 $\psi(x, t)$ の $t^m / m!$ の係数として定義する。

\[\psi(x, t) = \frac{te^{xt}}{e^t - 1} = \sum_{m=0}^\infty \frac{B_m (x)}{m!} t^m \tag{5.167}\]

ベルヌーイの関数 $B_m (x)$ はこの定義から次の性質を持つことが示される。

$(1)$ $B_m (x)$ は $m$ 次の多項式、 $B_m$ を 問題 5-8 のベルヌーイ数とすると

$\phantom{(1)~\qquad} B_m (x) = \sum_{l=0}^m \begin{pmatrix} m \cr l \end{pmatrix} B_l x^{m-l} \qquad \text{故に} \quad B_m (0) = B_m$

$\phantom{(1)}$ たとえば

$\phantom{(1)~\qquad} \displaystyle{B_0(x) = 1, \quad B_1(x) = x - \frac{1}{2}, \quad B_2(x) = x^2 - x + \frac{1}{6},}$

$\phantom{(1)~\qquad} \displaystyle{B_3(x) = x^3 - \frac{3}{2} x^2 + \frac{1}{2} x, \quad B_4(x) = x^4 - 2x^3 + x^2 - \frac{1}{30}, \quad \cdots}$

$(2)$ $B_m(1 - x) = (-1)^m B_m(x), \quad B_m (1 + x) - B_m (x) = mx^{m-1} \quad (m \geqq 0)$

$\phantom{(2)}$ これから

$\phantom{(2)~\qquad} B_{2m+1} = 0 \quad (m \geqq 1)$

$\phantom{(2)~\qquad} B_1 (1) = -B_1, \quad B_m (1) = B_m \quad (m \geqq 2)$

$(3)$ $B^\prime_m (x) = mB_{m-1} (x)$

$\phantom{(3)}$ これから $\int_0^1 B_m(x)~dx = 0 \quad (m \geqq 1)$

$(1)$ および $(3)$ から、

\[h \sum_{k=0}^n w_k B_m \left( \frac{k}{n} \right) = \int_{0}^{1} B_m (x)~dx = 0 \qquad (n \geqq m \geqq 1)\]

また $(2)$ および $w_{n-k} = w_k$ から

\[h \sum_{k=0}^n w_k B_{2m+1} \left( \frac{k}{n} \right) = 0 \qquad (m \geqq 0)\]

以上の関係を用いれば、ニュートン・コーツの閉型公式 (複合公式) $(5.69)$ より

\[\begin{align*} I_n & = h \sum_{k=0}^n w_k \sum_{j=0}^{N-1} f(a + (jn + k)h) \cr & = h \sum_{k=0}^n w_k [ e_{nNhD} - 1 ] \frac{e^{khD}}{e^{nhD} - 1} f(a) \cr & = [ e^{nNhD} - 1 ] \frac{h}{nhD} \sum_{k=0}^n w_k \psi \left( \frac{k}{n}, nhD \right) f(a) \cr & = [ e^{nNhD} - 1 ] \left\{ \frac{1}{D} + \sum_{l=[n/2]+1}^\infty \frac{n^{2l-1}}{(2l)!} h^{2l} \sum_{k=0}^n w_k B_{2l} \left( \frac{k}{n} \right) D^{2l-1} \right\} f(a) \cr & = \int_{a}^{b} f(x)~dx \cr & \phantom{=}+ \sum_{l=[n/2]+1}^\infty \frac{n^{2l-1}}{(2l)!} h^{2l} \sum_{k=0}^n w_k B_{2l} \left( \frac{k}{n} \right) [ f^{(2l-1)} (b) - f^{(2l-1)} (a) ] \tag{5.168} \end{align*}\]

ただし、 $[n / 2]$ は $n / 2$ を越えない最大の整数 (ガウスの記号) である。 この第2項の総和が打切り誤差 $(-En = In - I)$ である。 $n = 1$ とすれば $(5.166)$ が得られる。

また $b = a + nh$ として平均値の定理を適用すれば $(5.76)$ 、 $(5.77)$ が得られる。

---

問題 5-10

次の漸化式と直交関係を満たす $n$ 次直交多項式 $Q_n (x)~(n \geqq 0)$ がある。

\[\begin{align*} & Q_{n+1} (x) = (a_n x + b_n) Q_n (x) + c_n Q_{n-1} (x) \tag{5.169} \cr & \int_{a}^{b} w(x) Q_n (x) Q_m (x)~dx = \delta_{nm} \gamma_n \tag{5.170} \cr \end{align*}\]

$Q_n (x)$ を用いたガウス型積分公式の重みは、 $Q_n (x)$ の $x_n$ の係数を $A_n$ として、 $(5.117)$ 、すなわち

\[w_i = \int_{a}^{b} \frac{Q_n (x) w(x)}{(x - x_i) Q_n^\prime (x_i)}~dx = - \frac{ A_{n+1} \gamma_n}{A_n Q_{n+1} (x_i) Q_n^\prime (x_i)}\]

で与えられることを示せ (本文 $(5.117)$ )。 ただし、 $x_i$ は $Q_n (x)$ の $i$ 番目の零点 $Q_n (x_i) = 0$ である。

［解］

$Q_n (x)$ の最高次である $x_n$ の係数を $A_n$ とする。 $(5.169)$ の両辺の $x_{n+1}$ の係数を比べることにより

\[A_{n+1} = a_n A_n \qquad \therefore \quad a_n = A_{n+1} / A_n\]

漸化式 $(5.169)$ に $w(x) Q_{n-1} (x)$ を乗じて積分すると、

\[0 = \frac{a_n}{a_{n-1}} \gamma_{n} + c_n \gamma_{n-1} \qquad \therefore \quad c_n = - \frac{a_n \gamma_{n}}{a_{n-1} \gamma_{n-1}}\]

漸化式 $(5.169)$ に $Q_n (y)$ を乗じ、 $a_n \gamma_{n}$ で割れば

\[\frac{Q_{n+1} (x) Q_n (y)}{a_n \gamma_{n}} = \frac{(a_n x + bn) Q_n (x) Q_n (y)}{a_n \gamma_{n}} - \frac{Q_{n-1} (x) Q_n (y)}{a_{n-1} \gamma_{n-1}}\]

この式の $x$ と $y$ とを交換した式をつくり、この式から引くと、

\[\begin{aligned} & \frac{Q_{n+1} (x) Q_n (y) - Q_n (x) Q_{n+1} (y)}{a_n \gamma_{n}} \cr & - \frac{Q_n (x) Q_{n-1} (y) - Q_{n-1} (x) Q_n (y)}{a_{n-1} \gamma_{n-1}} = (x - y) \frac{ Q_n (x) Q_n (y)}{\gamma_{n}} \end{aligned}\]

左辺第1項を $S_n$ と書くと第2項は $S_{n-1}$ と書ける。 両辺の $n = 0$ から $n - 1$ までの和を取ると、 $S_{-1} = 0$ に注意して、

\[\frac{Q_{n+1} (x) Q_n (y) - Q_n (x) Q_{n+1} (y)}{a_n \gamma_{n}} = (x - y) \sum_{k=0}^{n-1} \frac{Q_k (x) Q_k (y)}{\gamma_{k}} \tag{5.171}\]

この式は クリストッフェル・ダルブー (Christoffel-Darbout) の恒等式 として知られている。 ここで、 $y = x_i$ とおいて、 $w(x) Q_0 (x) / (x - x_i)$ を乗じて積分すると、

\[- \frac{Q_{n+1} (x_i)}{a_n \gamma_{n}} \int_{a}^{b} \frac{Q_0 (x) Q_n (x) w(x)}{x - x_i}~dx = \frac{Q_0 (x_i)}{\gamma_{0}} \gamma_{0} = Q_0 (x_i)\]

ここに、 $Q_0 (x) = Q_0 (x_i) = \text{定数}$ であるから、

\[w_i = \int_{a}^{b} \frac{Q_n (x) w(x)}{(x - x_i) Q_n^\prime (x_i)}~dx = -\frac{a_n \gamma_{n}}{Q_{n+1} (x_i) Q_n^\prime (x_i)}\]

すなわち、 $(5.117)$ を得る。

---

問題 5-11

$\cos n \theta$ は $x = \cos \theta$ の多項式である。 $T_n (x) = \cos n \theta$ を $x = \cos \theta$ についてのチェビシェフ (Chebyshev) の多項式という。 たとえば、

\[T_0 (x) = \cos 0 = 1, \quad T_1(x) = \cos \theta = x, \quad T_2 (x) = \cos 2\theta = 2x^2 - 1, \cdots\]

などである。 チェビシェフの多項式の漸化式は、

\[\cos (n + 1) \theta = 2 \cos \theta \cos n \theta - \cos (n - 1) \theta\]

より、

\[T_{n+1}(x) = 2x T_n (x) - T_{n-1} (x) \tag{5.172}\]

である。 チェビシェフの多項式は、直交多項式である。 何故なら

\[\int_{0}^{\pi} \cos n\theta \cos m\theta d\theta = \begin{cases} 0 & n \neq m \cr \pi & n = m = 0 \cr \pi / 2 & n = m \neq 0 \cr \end{cases}\]

より、

\[\int_{-1}^{1} \frac{T_n (x) T_m (x)}{\sqrt{1 - x^2}}~dx = \frac{\pi}{2} (\delta_{nm} + \delta_{n0} \delta_{m0} ) \tag{5.173}\]

であるからである。 チェビシェフの多項式を用いたガウス型積分公式をガウス・チェビシェフの積分公式という。 この公式は次の通りであることを示せ。

\[\begin{array}{lcr} \hphantom{(5.888)} & \boxed{ \begin{aligned} \quad \cr & \text{ガウス・チェビシェフの積分公式} \cr & \quad \int_{-1}^{1} \frac{f(x)}{\sqrt{1 - x^2}}~dx = \frac{\pi}{n} \sum_{i=1}^n f(x_i) + E_n \cr & \quad \qquad x_i = \cos \theta_i, \quad \theta_i = \frac{2i - 1}{2n} \pi \quad (1 \leqq i \leqq n) \cr & \quad \qquad E_n = \frac{\pi f^{(2n)} (\xi)}{2^{2n-1} (2n)!} \qquad (-1 \leqq \xi \leqq 1) & \quad \cr \end{aligned} } & \begin{align*} \vphantom{\quad} \cr \vphantom{\text{ガウス・チェビシェフの積分公式}} \cr \vphantom{\sum_1^1} \tag{5.174} \cr \vphantom{\frac{()}{1}} \tag{5.175} \cr \vphantom{\frac{(^)}{()}} \tag{5.176} \cr \end{align*} \end{array}\]

［解］

チェビシェフの多項式 $T_n(x) = \cos n \theta$ の零点は $(5.175)$ で与えられる。 零点における値は

\[\begin{aligned} & T_n^\prime (x_i) = \frac{-n \sin n \theta_i}{- \sin \theta_i} = \frac{(-1)^{i+1} n}{\sin \theta_i} \cr & T_{n+1} (x_i) = \cos((n + 1) \theta_i) = (-1)^i \sin \theta_i \cr \end{aligned}\]

また $(5.172)$ より、 $A_{n+1} = 2 A_n$ 、 $A_n = 2^{n-1}$ 。 次に $(5.173)$ より $\gamma_{n} = \pi / 2$ 。 これらを、 $(5.117)$ および $(5.123)$ に代入すれば $(5.174)$ および $(5.176)$ が得られる。

この公式は分点である零点が正確にわかっていて、重みが全分点に共通であるから使いやすい。 上下限に特異性を持つが、大きな $n$ を用いない限り分点は上下限にないので問題はない。

この公式の PAD の図 5.12 は、分点と重みの計算がいらない分だけ、これまでのガウス型の積分公式にくらべて簡単になる。

---

問題 5-12

$n + 1$ 個の分点 $x_1, x_2, \cdots , x_{n+1}$ で $f(x)$ と一致する $n$ 次多項式 $p_{n}(x)$ を用いて、 $f(x)$ の $x = x_0$ における導関数 $f^{(l)} (x_0)(0 \leqq l \leqq n)$ の近似値を求める手順をつくれ。

［解］

$p_{n}(x)$ は、ニュートンの補間公式 $(5.47)$ より

\[p_{n}(x) = \sum_{k=0}^n \phi_k (x)f[x_1, x_2, \cdots x_{k+1}] \tag{5.177}\]

とする。 ここに、 $f[x_1, x_2, \cdots , x_{k+1}]$ は $k$ 階差分商で、また

\[\phi_k (x) = \prod_{l=1}^{k} (x - x_l) \tag{5.178}\]

である。 $\phi_k (x)$ を $(x - x_0)$ の多項式に書き直す。 $Z_l = x_0 - x_l$ と置くと、

\[\begin{align*} \phi_k (x) & = \prod_{l=1}^{k} (x - x_l) = \prod_{l=1}^{k} (x - x_0 + Z_l) \cr & = a_{k0} (x - x_0)^k + a_{k1} (x - x_0)^{k-1} + \cdots + a_{kk} \tag{5.179} \end{align*}\]

の係数は、 $0 \leqq k \leqq n$ に対して

\[\begin{aligned} & a_{k0} = 1 \cr & a_{k1} = Z_1 + Z_2 + \cdots + Z_k = a_{k-1~1} + a_{k-1~0} Z_k \cr & a_{k2} = Z_1 Z_2 + Z_1 Z_3 + \cdots + Z_{k-2} Z_k + Z_{k-1} Z_k = a_{k-1~2} + a_{k-1~1} Z_k \cr & \cdots \qquad \quad \qquad \cdots \cr & a_{kl} = Z_1 \cdots Z_l + \cdots + Z_{k-l+1} \cdots Z_k = a_{k-1~l} + a_{k-1~l-1} Z_k \cr & \cdots \qquad \quad \qquad \cdots \cr & a_{kk} = Z_1 Z_2 \cdots Z_k = a_{k-1~k-1} Z_k \cr \end{aligned}\]

であるから、 $a_{kl}(0 \leqq l \leqq k)$ は次の漸化式に従う：

\[\begin{aligned} & a_{k0} = 1 \cr & a_{kl} = a_{k-1~l} + a_{k-1~l-1} Z_k & (1 \leqq l \leqq k - 1) \cr & a_{kk} = a_{k-1~k-1} Z_k \cr \end{aligned} \tag{5.180}\]

$(5.179)$ を $(5.177)$ に代入すると

\[\begin{aligned} p_{n}(x) & = \sum_{k=0}^n \sum_{l=0}^{k} a_{kl}(x - x_0)^{k-l} f[x_1, x_2, \cdots x_{k+1}] \cr & = \sum_{k=0}^n \sum_{l=0}^{k} a_{k~k-l} f[x_1, x_2, \cdots x_{k+1}](x - x_0)^l \cr & = \sum_{l=0}^n \sum_{k=l}^{n} a_{k~k-l} f[x_1, x_2, \cdots x_{k+1}](x - x_0)^l \cr \end{aligned}\]

一方、 $f(x)$ を $n$ 次多項式として $x = x_0$ のまわりにテーラー展開すると

\[f(x) = \sum_{l=0}^n \frac{f^{(l)} (x_0)}{l!} (x - x_0)^l\]

であるから、 $p_{n}(x)$ と $f(x)$ の $(x - x_0)^l$ の係数をくらべて、導関数の近似値として

\[f^{(l)}(x_0) = l! \sum_{k=l}^n a_{k~k-l} f[x_1, x_2, \cdots x_{k+1}] \tag{5.181}\]

を得る。 $(5.180)$ と $(5.181)$ から、手順をつくればよい。

以上を基もとにして、 $x_1, x_2, \cdots x_N$ の $N$ 点のデータ $f(x_k)$ が与えられているとき、区間 $[ x_1, x_N ]$ 内の点 $x_0$ における $f(x_0)$ から $f^{(\text{L} \max)} (x_0)$ までの値を、 $N$ 点の中から $x_0$ の近傍の $m + 1$ 点を選び出して求める手順を求める。 当然 $\text{L} \max < m < N$ とする。 まず、与えられた点 $x_0$ の最近接点 $x_j$ を求める。 すなわち、 $x_j \leqq x_0 < x_{j+1}$ (この点 $x_j$ が $(5.181)$ の $x_1$ である)。 次の点 ( $(5.181)$ の $x_2$ ) は、 $x_{j+1}$ または $x_{j-1}$ である。 どちらの点を次の点とするかは、どちらを使った補間がより滑らかになるかによって決める。 すなわち、1階差分商の大きさにしたがって

$\qquad \vert f[x_{j-1}, x_j ] \vert < \vert f[x_j , x_{j+1}] \vert \qquad \text{ならば} \qquad x_{j-1}$

$\qquad \vert f[x_j , x_{j+1}] \vert \leqq \vert f[x_{j-1}, x_j ] \vert \qquad \text{ならば} \qquad x_{j+1}$

とする。 差分商は対称式であることに注意。 以下同様に、2階、3階、 $\cdots$ の差分商の小さい方の点を選択する。 いま、 $x_j < x_{j+1} < \cdots < x_{j+k-1}$ の $k$ 個の分点が選ばれてきたとき、次の分点は $x_{j-1}$ または $x_{j+k}$ であるが、

$\qquad \vert f[x_{j-1}, x_j , \cdots , x_{j+k-1}] \vert < \vert f[x_j , x_{j+1}, \cdots , x_{j+k}] \vert \text{ ならば}$

$\qquad \qquad \text{分点は } x_{j-1}\text{ 、差分商は } f[x_{j-1}, x_j , \cdots , x_{j+k-1}] \text{ を選択}$

$\qquad \vert f[x_j , x_{j+1}, \cdots , x_{j+k}] \vert \leqq \vert f[x_{j-1}, x_j , \cdots , x_{j+k-1}] \vert \text{ ならば}$

$\qquad \qquad \text{分点は } x_{j+k} \text{ 、差分商は } f[x_j , x_{j+1}, \cdots , x_{j+k}] \text{ を選択}$

する。 それでも、小さい方の差分商がある値より大きくなってしまったら、補間を中止する。 このような分点選択法は、 $k$ 階差分商が $k$ 階導関数に比例するから、変化のよりゆるやかな (より滑らかな) 方向の分点を選択するということである。 また、こうして得られる多項式は、(全領域ではなく) 評価点 $x_0$ を含む小区間で区分的に連続な区分多項式となる。 高次多項式による全域補間より、低次多項式による区分補間の方がかえって精度がよいことが多い。 (第8章図 8.1 はその例である。)

図 5.13 の PAD に手順を示す。 配列 DVDF(k,L) $=f [x_k, x_{k+1}, \cdots , x_{k+L}]$ はあらかじめ求めておく。 このとき、 $x < x_1$ および $x > x_N$ である (存在しない) 分点を探しに行かないように番兵 (Sentinel) を立てる (Sentnl>>1)。 $(5.181)$ の $a_{k~k-l}$ は、 $(5.180)$ の漸化式に従うが、すべての $k$ に共通に配列 b(L) $=a_{k~k-L}$ (一般に $a_{ij} =$ b(i-j) ) と置けば、

\[\begin{array}{lcl} a_{kk} = a_{k-1~k-1} Z_k & \text{は} & \texttt{b(0)* = Zk} \cr a_{k~k-L} = a_{k-1~k-L} + a_{k-1~k-L-1} Z_k & \text{は} & \texttt{b(L)* = Zk + b(L - 1)} \cr \end{array}\]

となる。 b(-1)=0 と置けば、第1の漸化式は第2の漸化式に含めることができる。

以上の手順が終了すれば、配列 df(i,L) に df(i,L) $=f^{(L)}(x_0)$ (ただし、 1 $\leqq$ i $\leqq$ Nd 、 0 $\leqq$ L $\leqq$ Lmax 、 $x_0$ は $i$ 番目の評価点) が求められている。