9.1 数値積分の基本解法

9.1.1 基本的な積分公式

積分区間 [a, b] を n 分割する分点列を \(a = x_0 < x_1 < \dots < x_n = b\), 各分点における重みを \(w_i\), 分点における関数値を \(f(x_i)\) としたときに, 積分公式の一般系は次のように表される. \[ I = \int_a^b f(x) dx \space \simeq \space I_n = \sum_{i=0}^n w_i f(x_i) \] 分点が積分区間の端点に一致しない公式もあり, N 個の分点 a < x_1 < x_2 < \dots < x_N < b により次のように表される. \[ I = \int_a^b f(x) dx \space \simeq I_N = \sum_{i=1}^N w_i f(x_i) \] 積分公式は重みと分点のとり方により次のように種々の公式が導かれる.

以上からわかるように, 積分公式は大きく 3 つに分けられる.

(1) 分点が等間隔のもの: ニュートン・コーツ公式 (台形則, シンプソン則を含む)
(2) 重みが一定のもの: チェビシェフ公式
(3) 分点間隔も重みも一定でないもの: ガウス型公式, クレンショウ・カーチス公式

いずれも, 分点 \(x_i\) において関数値 \(f(x_i)\) をとる補間多項式を定め, その積分値を求めるものである. これらは補間型積分公式とよばれる.

9.1.1.1 ニュートン・コーツ公式

区間 [a, b] を n 等分して分点を等間隔にとる方法である. n = 1 の場合は台形則, n = 2 の場合はシンプソン則とよばれる.

分点は \(a = x_0, x_1, \dots, x_n = b\) で間隔 \(h = (b – a)/n\) である. すなわち, \(x_k = a + kh\) となる.

被積分関数 \(f(x)\) を近似するラグランジュ補間多項式^(*) \(p_{n+1}(x)\) は次のように与えられる.
\[
p{n+1}(x) = \sum_{k=0}^n L_k^{(n)}(x) f(x_k)
\] (*) ラグランジュ補間多項式については「6. 補間法」を参照せよ.

これを積分すると,
\[
\begin{align}
I_n & = \int_a^b p_{n+1}(x) dx \\
& = \sum_{k=0}^n (\int_a^b L_k^{(n)}(x) dx) f(x_k) \\
& = \frac{b – a}{n} \sum_{k=0}^n w_k f(x_k) \\
& ただし, w_k = \frac{n}{b – a} \int_a^b L_k^{(n)}(x) dx \\
\end{align}
\] と積分公式の形になる. これを n 次のニュートン・コーツ公式という. この公式は n 次までの多項式に対して正確な値を与える.

重み \(w_k\) は多項式の計算により有理数 (分数) で求めることができる.

n = 1 ～ 10 のときの重みの値を「9.1.4 積分公式の係数表」に掲載した. 表のように n = 8 および n = 10 (実際には n \(\ge\) 10) では負の値が現れる. これは桁落ちの原因になるので n = 8 および n \(\ge\) 10 の公式は使わないほうがよいとされている.

一般的には, n を大きくして精度を上げるのではなく, 積分区間を小区間に分割してそれぞれに公式を適用するニュートン・コーツ複合公式が使用される. 例えば, n = 2 (シンプソン則) の場合, [a, b] を m 個の小区間に分割し, 各小区間でシンプソン則を適用すると次の計算式が得られる.
\[
\begin{align}
& I_{2m} = \frac{h}{3}(f(x_0) + 2 \sum_{j=1}^{m-1} f(x_{2j}) + 4 \sum_{j=1}^m f(x_{2j-1}) + f(x_{2m})) \\
& ただし, \space h = \frac{b – a}{2m}, \space x_j = a + jh \space (j = 0, 1, \dots, 2m) \\
\end{align}
\] これは複合シンプソン公式とよばれる.

9.1.1.2 チェビシェフ公式

ニュートン・コーツ公式とは逆に重みが一定になるように分点を選ぶ方法である.

チェビシェフ公式は次のように表される.
\[
I_N = \frac{b – a}{N} \sum_{k=1}^N f(\frac{a + b}{2} + \frac{(b – a)}{2}x_k)
\] \(x_k\) は分点で, 多項式 \(p_N(x) = \prod_{k=1}^N (x – x_k) = \sum_{i=0}^N a_i x^{N-i}\) をゼロとおいた方程式の N 個の解として求められる. 展開した多項式の係数 \(a_i\) は次のように与えられる (導出の詳細は文献[3]を参照).
\[
a_0 = 1 \\
\begin{equation}
a_i =
\begin{cases}
-\frac{N}{i} \sum_{j=2,4,\dots}^{i} \frac{a_{i-j}}{j + 1} & (i = 2, 4, \dots, N2) \\
0 & (i = 1, 3, \dots, N1) \\
\end{cases}
\end{equation}
\] ただし, N1 = N – 1, N2 = N (N が偶数のとき), N1 = N, N2 = N – 1 (N が奇数のとき).

N = 2 ～ 7, 9 のときの分点の値を「9.1.4 積分公式の係数表」に掲載した. N = 8 と N \(\ge\) 10 は方程式の解に複素数が含まれ, 計算が複雑になるうえ精度の点でも不利なため通常は使用されない.

9.1.1.3 ガウス型公式

できるだけ精度が高くなるように分点と重みの両方を選ぶ方法である. 被積分関数を直交多項式補間^(*)することにより求められ, 2N – 1 次までの多項式に対して正確な値を与える.

(*) 直交多項式および直交多項式補間については「6. 補間法」を参照せよ.

区間 [a, b] における f(x)w(x) の積分を考える. w(x) は密度関数である.
\[
I = \int_a^b f(x) w(x) dx
\] f(x) を直交多項式補間関数 \(f_N(x)\) で置き換えると次のようになる.
\[
\begin{align}
I_N & = \int_a^b f_N(x)w(x) dx \\
& = \sum_{j=0}^{N-1} \frac{1}{\lambda_j} \sum_{k=1}^N w_k p_j(x_k) f(x_k) \int_a^b p_j(x) w(x) dx \\
\end{align}
\] \(p_j(x)\) の直交性より上式の積分の部分は次のように簡単になる.
\[
\begin{equation}
\int_a^b p_j(x) w(x) dx =
\begin{cases}
\frac{\lambda_0}{\mu_0} & (j = 0) \\
0 & (j \ne 0) \\
\end{cases}
\end{equation}
\] すなわち, \(I_N\) の式において外側の \(\sigma\) は \(j = 0\) 以外の項は \(0\) になり, また \(p_0(x) = \mu_0\) であるから, 次のガウス型積分公式が得られる.
\[
\begin{align}
I_N & = \frac{1}{\lambda_0} \frac{\lambda_0}{\mu_0} \sum_{k=1}^N w_k p_0(x_k) f(x_k) \\
& = \sum_{k=1}^N w_k f(x_k) \\
\end{align}
\] ここで, 分点 \(x_k\) は直交多項式のゼロ点である. 重み \(w_k\) は直交多項式補間の \(w_k\) に等しく, 次式で与えられる.
\[
w_k = \frac{\mu_N \lambda_{N-1}}{\mu_{N-1} p_{N-1}(x_k) p’_N(x_k)}
\] 直交多項式の種類に応じていくつかの積分公式があり, 代表的なものは次のとおりである.

(1) ガウス・ルジャンドル公式 (ガウス公式)

ルジャンドル多項式を使用して, \(p_n(x) = P_n(x), w(x) = 1, [a, b] = [-1, 1]\) とするとガウス・ルジャンドル公式が得られる. ガウス・ルジャンドル公式は単にガウス公式とよばれることが多い.
\[
I = \int_{-1}^{1} f(x) dx \simeq I_N = \sum_{k=1}^N w_k f(x_k)
\] ただし, \(x_k\) は \(P_n(x) = 0\) の N 個のゼロ点である. \(w_k\) は, \(\lambda_N = 2/(2N + 1)\), \(\mu_N = \prod_{k=1}^N (2k – 1)/k\) から, 次のように求めることができる.
\[
\begin{align}
w_k & = \frac{2}{nP_{N-1}(x_k)P’_N(x_k)} \\
& = \frac{2(1 – x_k^2)}{(nP_{N-1}(x_k))^2} \\
\end{align}
\] N = 2 ～ 10 のときの分点および重みの値を「9.1.4 積分公式の係数表」に掲載した.

任意の区間 [a, b] の積分を求めるときは次のように変数変換して計算するとよい.
\[
I = \int_a^b f(x) dx \simeq I_N = \frac{b – a}{2} \sum_{k=1}^N w_k f(\frac{a + b}{2} + \frac{b – a}{2}x_k)
\]

(2) ガウス・ラゲール公式

ラゲール多項式を使用して, \(p_n(x) = L_n(x), w(x) = e^{-x}, [a, b] = [0, +\infty]\) とするとガウス・ラゲール公式が得られる.
\[
I = \int_0^{\infty} e^{-x} f(x) dx \simeq I_N = \sum_{k=1}^N w_k f(x_k)
\] ただし, \(x_k\) は \(L_n(x) = 0\) の N 個のゼロ点である. \(w_k\) は, \(\lambda_N = 1, \mu_N = (-1)^n/n!\) より, 次のように求めることができる.
\[
w_k = -\frac{1}{NL_{N-1}(x_k)L’_N(x_k)} = \frac{x_k}{n^2 L_{N-1}(x_k)^2}
\] N = 2 ～ 10 のときの分点および重みの値を「9.1.4 積分公式の係数表」に掲載した.

(3) ガウス・エルミート公式

エルミート多項式を使用して, \(p_n(x) = H_n(x), w(x) = e^{-x^2}, [a, b] = [-\infty, +\infty]\) とするとガウス・エルミート公式が得られる.
\[
I = \int_{-\infty}^{+\infty} f(x) e^{-x^2} dx \simeq I_N = \sum_{k=1}^N w_k f(x_k)
\] ただし, \(x_k\) は \(H_n(x) = 0\) の N 個のゼロ点である. \(w_k\) は, \(\lambda_N = \pi^{1/2} 2^n n!, \mu_N = 2^n\) より, 次のように求めることができる.
\[
w_k = \frac{\pi^{1/2} 2^n (n – 1)!}{H_{N-1}(x_k) H’_N(x_k)} = \frac{\pi^{1/2} 2^{n-1} (n – 1)!} {nH_{N-1}(x_k)^2}
\] N = 2 ～ 10 のときの分点および重みの値を「9.1.4 積分公式の係数表」に掲載した.

(4) クレンショウ・カーチス公式

(i) 滑らかな関数の積分

被積分関数 f(x) をチェビシェフ級数展開し, 項別に積分する方法である. チェビシェフ級数の収束が速い関数であれば少ない関数計算回数で積分を求めることができる.

区間 [-1, 1] で滑らかな関数 f(x) をチェビシェフ級数展開する.
\[
\begin{align}
& f(x) = \frac{1}{2}a_0 + \sum_{k=1}^{\infty} a_k T_k(x) \\
& T_k(x) = cos(k\theta), x = cos(\theta) \\
& a_k = \frac{2}{\pi} \int_0^{\pi} f(cos(\theta))cos(k\theta) d\theta \\
\end{align}
\] 係数 \(a_k\) の積分を台形則で計算すると次のようになる.
\[
a_k^{(n)} = \frac{1}{n}f(1) + \frac{2}{n}\sum_{j=1}^{n-1} f(cos(\frac{j\pi}{n})) cos(\frac{kj\pi}{n}) + (-1)^k\frac{1}{n}f(-1)
\] 右上添字の (n) は, n 小区間の複合台形則を使ったことを示す. この係数は FFT (コサイン変換) プログラムを使って効率的に計算することができる (XLPack であれば Cost1f が使える).

f(x) のチェビシェフ級数展開を第 n 項までで打ち切って積分を次のように近似する.
\[
I = \int_a^b f(x) dx \simeq I_n = \frac{1}{2}a_0^{(n)} \int_{-1}^1 T_0(x) dx + \sum_{k=1}^{n-1} a_k^{(n)} \int_{-1}^1 T_k(x) dx + \frac{1}{2}a_n^{(n)} \int_{-1}^{1} T_n(x) dx
\] このように項別の積分の形になる. チェビシェフ多項式の項の積分は次の関係を使って計算することができる (積分定数は省略して表示した).
\[
\begin{equation}
\int T_k dx =
\begin{cases}
T_1(x) & (k = 0) \\
\frac{1}{4}(T_2(x) + 1) & (k = 1) \\
\frac{1}{2}(\frac{T_{k+1}(x)}{k + 1} – \frac{T_{k-1}(x)}{k – 1}) & (k \ge 2) \\
\end{cases}
\end{equation}
\] 以上の関係式を使って f(x) の [-1, 1] における積分を計算することができる (この場合, k が奇数の \(\int_{-1}^1 T_k dx\) は 0 になる).

任意の区間 [a, b] の積分を求めるには次のように変数変換して計算するとよい.
\[
x = \frac{b – a}{2} cos(\theta) + \frac{a + b}{2}
\]

(ii) 特異性が強い関数の積分

クレンショウ・カーチス公式で \(f(x)w(x)\) の形の関数を積分することを考える. \(w(x)\) は特異性が強いような関数を想定して, \(f(x)\) だけをチェビシェフ級数展開し \(w(x)\) を含む項別積分は別途行うことにする.
\[
\begin{align}
I & = \int_{-1}^{1} f(x)w(x)dx \\
& \simeq I_n = \frac{1}{2}a_0^{(n)}\int_{-1}^{1}T_0(x)w(x)dx + \sum_{k=1}^{n-1}a_k^{(n)}\int_{-1}^{1}T_k(x)w(x)dx + \frac{1}{2}a_n^{(n)}\int_{-1}^1 T_n(x)w(x)dx \\
& = \frac{1}{2}a_0^{(n)}\mu_0 + \sum_{k=1}^{n-1}a_k^{(n)}\mu_k + \frac{1}{2}a_n^{(n)}\mu_n \\
& ただし, \space \mu_k = \int_{-1}^{1} T_k(x)w(x) dx \\
\end{align}
\] \(\mu_k\) はモーメントとよばれる. モーメントが安定的に計算できる (例えば, 解析的に) ような \(w(x)\) を含む積分については特異性が強くてもクレンショウ・カーチス公式で計算できることになる.

次の積分を考える.
\[
I = \int_{-1}^{1} \frac{f(x)}{x – c} dx
\] これはコーシーの主値積分である.

\(w(x) = 1/(x – c)\) についてモーメントは次のように計算できる (文献[6]参照).
\[
\mu_0 = ln|\frac{1 – c}{1 + c}| \\
\mu_1 = 2 + c\mu_0 \\
\begin{equation}
\mu_{k+1} – 2c\mu_k + \mu_{k-1} =
\begin{cases}
0 & (k が奇数) \\
\frac{4}{1 – k^2} & (k が偶数) \\
\end{cases}
\end{equation}
\] \(c\) の値が \([-1-\epsilon, 1+\epsilon]\) (\(\epsilon\) はおおむね 0.1) の中にあれば上式は安定に計算できる. そうでなければ, \(c\) は積分範囲から十分に遠くにあるとみなして通常のガウス公式などで計算できる.

文献[6]では他に以下の関数についての計算法も示されており, これらの計算法は QUADPACK のサブルーチンで採用されている.

• \(sin(\omega x)\) および \(cos(\omega x)\)
• \((x – a)^{\alpha} (b – x)^{\beta} log^{\mu}(x – a) log^{\nu}(b – x)\) (\(\mu\), \(\nu\) = 0 または 1)

9.1.2 ロンバーグ積分

数列 \(a_k\) が定数 \(a\) に収束するものとして, 収束率 \(\lambda\)を
\[
\lambda = \lim_{k \to \infty} \frac{a_{k+1} – a}{a_k – a}
\] と定義する. \(\lambda\) の値が事前にわかっているとして, これを \(a\) について解いて求めた値
\[
b_k = \frac{a_{k+1} – \lambda a_k}{1 – \lambda}
\] による数列 {\(b_k\)} は元の数列 {\(a_k\)} よりも速く収束すると予想される. これをリチャードソンの加速法という.

複合台形則において分点数を倍々にした系列 \(I_1, I_2, \dots\) を作れば, h は半々になってゆき, 誤差は 1 段進むごとに 1/4 になる (すなわち収束率が 1/4: 台形則の誤差の項を参照せよ). これにリチャードソンの加速法を適用した補外を行い,
\[
I_k^{(1)} = \frac{I_{k+1} – \lambda I_k}{1 – \lambda} = \frac{4 I_{k+1} – I_k}{4 – 1}
\] とした系列 { \(I_1^{(1)}, I_2^{(1)}, \dots\) } を得ることができる. これはシンプソン則に一致していることがわかる. 収束率 \(\lambda^{(1)}\) は \(1/4^2\) になる. 更に収束を加速するために新しい系列に第 2 回補外を行うと
\[
I_k^{(2)} = \frac{I_{k+1}^{(1)} – \lambda^{(1)} I_k^{(1)}}{1 – \lambda^{(1)}} = \frac{4^2 I_{k+1}^{(1)} – I_k^{(1)}}{4^2 – 1}
\] とした系列 { \(I_1^{(2)}, I_2^{(2)}, \dots\) } を得る. このように補外を繰り返して収束を加速する方法をロンバーグ積分法という.

複合台形則において分点数を倍々にしていく際には, 毎回半分は倍にする前の分点と一致するので半分だけ計算すればよいことになる.

9.1.3 数値実験

ここまで説明した積分公式についていくつかの例題関数を使って計算精度を比較する. なお, 計算は VBA の Double 型を使って 64 ビット倍精度 (10 進 15 ～ 16 桁) で行った.

9.1.3.1 有限区間の積分

数値実験 (1) 性質の良い関数

次の積分を考える.
\[
I = \int_0^1 e^x cos(x) dx = 1.378024613547364
\] これは性質の良い (解析的な) 関数で, グラフに表すと次のような形をしている.

これを種々の積分公式で分点数 N (= 関数評価回数) (横軸) を変えて計算し, その時の相対誤差 (縦軸) を対数グラフにすると次のようになった.

台形則は小区間数を N – 1 とした複合公式の結果を示した. シンプソン則は小区間数を (N – 1)/2 とした複合公式の結果を示した. ニュートン・コーツ公式およびクレンショウ・カーチス公式では N = n + 1 としてプロットした.

ニュートン・コーツ公式およびチェビシェフ公式では分点数が偶数よりも奇数のほうが成績がよい.

ガウス公式は他に比べて精度がよく, 可能であればこの公式を使うのがよい. N = 7 以上で頭打ちなのは倍精度計算の丸め誤差のためである.

分点を倍々に増やしていったときの台形則の計算結果に着目すると次のようになっている.

収束率はほぼ 1/4 になっており, 「9.1.2 ロンバーグ積分」で示した収束率が成り立っていることがわかる.

ロンバーグ積分は簡単な仕組みで, 台形則やシンプソン則と比べると手間はあまり変わらないが精度はよいので, これを使う価値はある.

数値実験 (2) 積分区間外に特異点を持つ関数

積分区間外に特異点を持つ次の関数の積分を考える.
\[
I = \int_0^4 \frac{1}{1 + x} dx = ln(1 + 4) = ln(5)
\] 次のように積分区間 [0, 4] では滑らかな関数に見えるが, 近くの x = -1 に特異点を持つ.

これを上と同じように種々の積分公式で計算し誤差を対数グラフにすると次のようになった.

図のようにどの公式でも精度が落ちたのがわかる. ガウス公式でも N = 18 以上でないと最高精度は得られない. 積分区間外の特異点が計算精度に明らかに影響しているのがわかる.

数値実験 (3) 複素平面に特異点を持つ関数

次に以下の積分を考える.
\[
I = \int_0^4 \frac{1}{1 + x^2} dx = arctan(4)
\] これは実数軸上には特異点がないが, 複素平面では近く(x = i) に特異点がある.

これを上と同じように種々の積分公式で計算し誤差を対数グラフにすると次のようになった.

これも一つ前の例と同じようにどの公式でも精度が落ちた. 複素平面の特異点も影響するのがわかる.

数値実験 (4) 端点に特異点を持つ関数

次に端点に特異点があるとどうなるか見てみる.
\[
I = \int_0^1 \sqrt{1 – x^2} dx = \frac{\pi}{4}
\] x = 1 が特異点である.

これを上と同じように種々の積分公式で計算し誤差を対数グラフにすると次のようになった.

このように端点とはいえ区間内に特異点がある場合にはどの公式でも精度が得られなかった.

この例における台形則の収束率は次のようになった.

収束率は 1/4 ではなく 1/2.8 にしかなっておらず, 正常な動きをしていないことを示唆している. これに伴ってロンバーグ積分の加速の効果が小さくなっていることがわかる.

9.1.3.2 無限区間, 半無限区間の積分

数値実験 (5) 半無限区間の積分

次の積分を考える.
\[
I = \int_a^{\infty} \frac{e^{-x}}{x} = E1(a)
\] E1(a) はよく使われる特殊関数のひとつの指数積分である. グラフに表すと次のような形をしている.

ガウス・ラゲール公式による a = 1 のとき (E1(1)) の計算結果は次のとおりである.

数値実験 (6) 無限区間の積分

次の積分を考える.
\[
I = \int_{-\infty}^{+\infty} \frac{e^{-x^2}}{1 + x^2} = 1.3432934216
\] 関数のグラフは次のとおりである.

ガウス・エルミート公式による計算結果は次のとおりである.

9.1.4 積分公式の係数表

(1) ニュートン・コーツ公式

\[
\begin{align}
& \int_a^b f(x) dx = h \sum_{k=0}^n w_k f(x_k) \\
& ただし, h = \frac{b – a}{n} \\
\end{align}
\] ニュートン・コーツ公式における係数 \(w_i\) を以下に示す. 本文で説明した式を数式処理ソフトを使って有理数のまま計算した.

  
         1  1
nc(1) = [-, -]
         2  2

         1  4  1
nc(2) = [-, -, -]
         3  3  3

         3  9  9  3
nc(3) = [-, -, -, -]
         8  8  8  8

         14  64  8   64  14
nc(4) = [--, --, --, --, --]
         45  45  15  45  45

         95   125  125  125  125  95
nc(5) = [---, ---, ---, ---, ---, ---]
         288  96   144  144  96   288

         41   54  27   68  27   54  41
nc(6) = [---, --, ---, --, ---, --, ---]
         140  35  140  35  140  35  140

         5257   25039  343  20923  20923  343  25039  5257
nc(7) = [-----, -----, ---, -----, -----, ---, -----, -----]
         17280  17280  640  17280  17280  640  17280  17280

         3956   23552    3712   41984    3632  41984    3712   23552  3956
nc(8) = [-----, -----, - -----, -----, - ----, -----, - -----, -----, -----]
         14175  14175    14175  14175    2835  14175    14175  14175  14175

         25713  141669  243   10881  26001  26001  10881  243   141669  25713
nc(9) = [-----, ------, ----, -----, -----, -----, -----, ----, ------, -----]
         89600  89600   2240  5600   44800  44800  5600   2240  89600   89600

          80335   132875    80875  28375    24125  89035    24125  28375    80875  132875  80335
nc(10) = [------, ------, - -----, -----, - -----, -----, - -----, -----, - -----, ------, ------]
          299376  74844     99792  6237     5544   12474    5544   6237     99792  74844   299376

 

(2) チェビシェフ公式

\[
\int_a^b f(x) dx = \frac{b – a}{N} \sum_{i=1}^N f(\frac{a + b}{2} + \frac{(b – a)x_i}{2})
\] チェビシェフ公式における分点 \(x_i\) を以下に示す. 拡張倍精度 (80 ビット浮動小数) の計算により本文で説明した多項式の方程式の解を求め, 結果の 16 桁を表示した.

  
N = 2
  xi: -5.773502691896258E-0001  5.773502691896258E-0001
N = 3
  xi: -7.071067811865475E-0001  0.000000000000000E+0000  7.071067811865475E-0001
N = 4
  xi: -7.946544722917661E-0001 -1.875924740850799E-0001  1.875924740850799E-0001
       7.946544722917661E-0001
N = 5
  xi: -8.324974870009819E-0001 -3.745414095535811E-0001  0.000000000000000E+0000
       3.745414095535811E-0001  8.324974870009819E-0001
N = 6
  xi: -8.662468181078206E-0001 -4.225186537611115E-0001 -2.666354015167047E-0001
       2.666354015167047E-0001  4.225186537611115E-0001  8.662468181078206E-0001
N = 7
  xi: -8.838617007580490E-0001 -5.296567752851568E-0001 -3.239118105199076E-0001
       0.000000000000000E+0000  3.239118105199076E-0001  5.296567752851568E-0001
       8.838617007580490E-0001
N = 9
  xi: -9.115893077284345E-0001 -6.010186553802381E-0001 -5.287617830578800E-0001
      -1.679061842148039E-0001  0.000000000000000E+0000  1.679061842148039E-0001
       5.287617830578800E-0001  6.010186553802381E-0001  9.115893077284345E-0001

 

(3) ガウス・ルジャンドル公式 (ガウス公式)

\[
\int_a^b f(x) dx = \frac{(b – a}{2} \sum_{k=1}^N w_kf(\frac{a + b}{2} + \frac{b – a}{2}x_k)
\] ガウス・ルジャンドル公式における分点 \(x_i\) および重み \(w_i\) を以下に示す. 分点としてルジャンドル多項式のゼロ点を求め, それから重みを計算した. 計算は拡張倍精度 (80 ビット浮動小数) により行い結果の 16 桁を表示した.

  
N = 2
  xi: -5.773502691896258E-0001  5.773502691896258E-0001
  wi:  1.000000000000000E+0000  1.000000000000000E+0000
N = 3
  xi: -7.745966692414834E-0001  0.000000000000000E+0000  7.745966692414834E-0001
  wi:  5.555555555555556E-0001  8.888888888888889E-0001  5.555555555555556E-0001
N = 4
  xi: -8.611363115940526E-0001 -3.399810435848563E-0001  3.399810435848563E-0001
       8.611363115940526E-0001
  wi:  3.478548451374539E-0001  6.521451548625461E-0001  6.521451548625461E-0001
       3.478548451374539E-0001
N = 5
  xi: -9.061798459386640E-0001 -5.384693101056831E-0001  0.000000000000000E+0000
       5.384693101056831E-0001  9.061798459386640E-0001
  wi:  2.369268850561891E-0001  4.786286704993665E-0001  5.688888888888889E-0001
       4.786286704993665E-0001  2.369268850561891E-0001
N = 6
  xi: -9.324695142031520E-0001 -6.612093864662645E-0001 -2.386191860831969E-0001
       2.386191860831969E-0001  6.612093864662645E-0001  9.324695142031520E-0001
  wi:  1.713244923791703E-0001  3.607615730481386E-0001  4.679139345726910E-0001
       4.679139345726910E-0001  3.607615730481386E-0001  1.713244923791703E-0001
N = 7
  xi: -9.491079123427585E-0001 -7.415311855993944E-0001 -4.058451513773972E-0001
       0.000000000000000E+0000  4.058451513773972E-0001  7.415311855993944E-0001
       9.491079123427585E-0001
  wi:  1.294849661688697E-0001  2.797053914892767E-0001  3.818300505051189E-0001
       4.179591836734694E-0001  3.818300505051189E-0001  2.797053914892767E-0001
       1.294849661688697E-0001
N = 8
  xi: -9.602898564975362E-0001 -7.966664774136267E-0001 -5.255324099163290E-0001
      -1.834346424956498E-0001  1.834346424956498E-0001  5.255324099163290E-0001
       7.966664774136267E-0001  9.602898564975362E-0001
  wi:  1.012285362903763E-0001  2.223810344533745E-0001  3.137066458778873E-0001
       3.626837833783620E-0001  3.626837833783620E-0001  3.137066458778873E-0001
       2.223810344533745E-0001  1.012285362903763E-0001
N = 9
  xi: -9.681602395076261E-0001 -8.360311073266358E-0001 -6.133714327005904E-0001
      -3.242534234038089E-0001  0.000000000000000E+0000  3.242534234038089E-0001
       6.133714327005904E-0001  8.360311073266358E-0001  9.681602395076261E-0001
  wi:  8.127438836157441E-0002  1.806481606948574E-0001  2.606106964029355E-0001
       3.123470770400028E-0001  3.302393550012598E-0001  3.123470770400028E-0001
       2.606106964029355E-0001  1.806481606948574E-0001  8.127438836157441E-0002
N = 10
  xi: -9.739065285171717E-0001 -8.650633666889845E-0001 -6.794095682990244E-0001
      -4.333953941292472E-0001 -1.488743389816312E-0001  1.488743389816312E-0001
       4.333953941292472E-0001  6.794095682990244E-0001  8.650633666889845E-0001
       9.739065285171717E-0001
  wi:  6.667134430868814E-0002  1.494513491505806E-0001  2.190863625159820E-0001
       2.692667193099964E-0001  2.955242247147529E-0001  2.955242247147529E-0001
       2.692667193099964E-0001  2.190863625159820E-0001  1.494513491505806E-0001
       6.667134430868814E-0002

 

(4) ガウス・ラゲール公式

\[
\int_0^{\infty} f(x)e^{-x} dx = \sum_{k=1}^N w_k f(x_k)
\] ガウス・ラゲール公式における分点 \(x_i\) および重み \(w_i\) を以下に示す. 分点としてラゲール多項式のゼロ点を求め, それから重みを計算した. 計算は拡張倍精度 (80 ビット浮動小数) により行い結果の 16 桁を表示した.

  
N = 2
  xi: 5.857864376269050E-0001 3.414213562373095E+0000
  wi: 8.535533905932738E-0001 1.464466094067262E-0001
N = 3
  xi: 4.157745567834791E-0001 2.294280360279042E+0000 6.289945082937479E+0000
  wi: 7.110930099291730E-0001 2.785177335692408E-0001 1.038925650158614E-0002
N = 4
  xi: 3.225476896193923E-0001 1.745761101158347E+0000 4.536620296921128E+0000
      9.395070912301133E+0000
  wi: 6.031541043416336E-0001 3.574186924377997E-0001 3.888790851500538E-0002
      5.392947055613275E-0004
N = 5
  xi: 2.635603197181409E-0001 1.413403059106517E+0000 3.596425771040722E+0000
      7.085810005858838E+0000 1.264080084427578E+0001
  wi: 5.217556105828087E-0001 3.986668110831759E-0001 7.594244968170760E-0002
      3.611758679922048E-0003 2.336997238577623E-0005
N = 6
  xi: 2.228466041792607E-0001 1.188932101672623E+0000 2.992736326059314E+0000
      5.775143569104511E+0000 9.837467418382590E+0000 1.598287398060170E+0001
  wi: 4.589646739499636E-0001 4.170008307721210E-0001 1.133733820740450E-0001
      1.039919745314907E-0002 2.610172028149321E-0004 8.985479064296212E-0007
N = 7
  xi: 1.930436765603624E-0001 1.026664895339192E+0000 2.567876744950746E+0000
      4.900353084526485E+0000 8.182153444562861E+0000 1.273418029179781E+0001
      1.939572786226254E+0001
  wi: 4.093189517012739E-0001 4.218312778617198E-0001 1.471263486575053E-0001
      2.063351446871694E-0002 1.074010143280746E-0003 1.586546434856420E-0005
      3.170315478995581E-0008
N = 8
  xi: 1.702796323051010E-0001 9.037017767993799E-0001 2.251086629866131E+0000
      4.266700170287659E+0000 7.045905402393466E+0000 1.075851601018100E+0001
      1.574067864127800E+0001 2.286313173688926E+0001
  wi: 3.691885893416375E-0001 4.187867808143430E-0001 1.757949866371718E-0001
      3.334349226121565E-0002 2.794536235225673E-0003 9.076508773358213E-0005
      8.485746716272532E-0007 1.048001174871510E-0009
N = 9
  xi: 1.523222277318082E-0001 8.072200227422558E-0001 2.005135155619347E+0000
      3.783473973331233E+0000 6.204956777876613E+0000 9.372985251687576E+0000
      1.346623691109209E+0001 1.883359778899170E+0001 2.637407189092738E+0001
  wi: 3.361264217979625E-0001 4.112139804239844E-0001 1.992875253708856E-0001
      4.746056276565160E-0002 5.599626610794583E-0003 3.052497670932106E-0004
      6.592123026075352E-0006 4.110769330349548E-0008 3.290874030350708E-0011
N = 10
  xi: 1.377934705404924E-0001 7.294545495031705E-0001 1.808342901740316E+0000
      3.401433697854900E+0000 5.552496140063804E+0000 8.330152746764497E+0000
      1.184378583790007E+0001 1.627925783137810E+0001 2.199658581198076E+0001
      2.992069701227389E+0001
  wi: 3.084411157650201E-0001 4.011199291552736E-0001 2.180682876118094E-0001
      6.208745609867775E-0002 9.501516975181101E-0003 7.530083885875388E-0004
      2.825923349599566E-0005 4.249313984962686E-0007 1.839564823979631E-0009
      9.911827219609009E-0013

 

(5) ガウス・エルミート公式

\[
\int_{-\infty}^{+\infty} f(x) e^{-x^2} dx = \sum_{k=1}^N w_k f(x_k)
\] ガウス・エルミート公式における分点 \(x_i\) および重み \(w_i\) を以下に示す. 分点としてエルミート多項式のゼロ点を求め, それから重みを計算した. 計算は拡張倍精度 (80 ビット浮動小数) により行い結果の 16 桁を表示した.

  
N = 2
  xi: -7.071067811865475E-0001  7.071067811865475E-0001
  wi:  8.862269254527580E-0001  8.862269254527580E-0001
N = 3
  xi: -1.224744871391589E+0000  0.000000000000000E+0000  1.224744871391589E+0000
  wi:  2.954089751509193E-0001  1.181635900603677E+0000  2.954089751509193E-0001
N = 4
  xi: -1.650680123885785E+0000 -5.246476232752903E-0001  5.246476232752903E-0001
       1.650680123885785E+0000
  wi:  8.131283544724518E-0002  8.049140900055128E-0001  8.049140900055128E-0001
       8.131283544724518E-0002
N = 5
  xi: -2.020182870456086E+0000 -9.585724646138185E-0001  0.000000000000000E+0000
       9.585724646138185E-0001  2.020182870456086E+0000
  wi:  1.995324205904591E-0002  3.936193231522412E-0001  9.453087204829419E-0001
       3.936193231522412E-0001  1.995324205904591E-0002
N = 6
  xi: -2.350604973674492E+0000 -1.335849074013697E+0000 -4.360774119276165E-0001
       4.360774119276165E-0001  1.335849074013697E+0000  2.350604973674492E+0000
  wi:  4.530009905508846E-0003  1.570673203228566E-0001  7.246295952243925E-0001
       7.246295952243925E-0001  1.570673203228566E-0001  4.530009905508846E-0003
N = 7
  xi: -2.651961356835233E+0000 -1.673551628767471E+0000 -8.162878828589647E-0001
       0.000000000000000E+0000  8.162878828589647E-0001  1.673551628767471E+0000
       2.651961356835233E+0000
  wi:  9.717812450995192E-0004  5.451558281912703E-0002  4.256072526101278E-0001
       8.102646175568073E-0001  4.256072526101278E-0001  5.451558281912703E-0002
       9.717812450995192E-0004
N = 8
  xi: -2.930637420257244E+0000 -1.981656756695843E+0000 -1.157193712446780E+0000
      -3.811869902073221E-0001  3.811869902073221E-0001  1.157193712446780E+0000
       1.981656756695843E+0000  2.930637420257244E+0000
  wi:  1.996040722113676E-0004  1.707798300741348E-0002  2.078023258148919E-0001
       6.611470125582413E-0001  6.611470125582413E-0001  2.078023258148919E-0001
       1.707798300741348E-0002  1.996040722113676E-0004
N = 9
  xi: -3.190993201781528E+0000 -2.266580584531843E+0000 -1.468553289216668E+0000
      -7.235510187528376E-0001  0.000000000000000E+0000  7.235510187528376E-0001
       1.468553289216668E+0000  2.266580584531843E+0000  3.190993201781528E+0000
  wi:  3.960697726326438E-0005  4.943624275536947E-0003  8.847452739437657E-0002
       4.326515590025558E-0001  7.202352156060510E-0001  4.326515590025558E-0001
       8.847452739437657E-0002  4.943624275536947E-0003  3.960697726326438E-0005
N = 10
  xi: -3.436159118837738E+0000 -2.532731674232790E+0000 -1.756683649299882E+0000
      -1.036610829789514E+0000 -3.429013272237046E-0001  3.429013272237046E-0001
       1.036610829789514E+0000  1.756683649299882E+0000  2.532731674232790E+0000
       3.436159118837738E+0000
  wi:  7.640432855232621E-0006  1.343645746781233E-0003  3.387439445548106E-0002
       2.401386110823147E-0001  6.108626337353258E-0001  6.108626337353258E-0001
       2.401386110823147E-0001  3.387439445548106E-0002  1.343645746781233E-0003
       7.640432855232621E-0006