9. 数値積分

はじめに

有限区間あるいは無限区間における定積分の近似計算法を数値積分法 (numerical intergation) という. または 求積法 (quadrature) ということもある.

数値積分プログラムは大きく 2 種類がある. 1 つは関数入力タイプで, 積分を計算する関数 (被積分関数) の値が要求に応じて区間内の任意の点で計算できるプログラムとして与えられる場合でである. もう 1 つはデータ入力タイプで, 被積分関数の値がとびとびの点 (等間隔の場合と不等間隔の場合がある) におけるデータとして与えられる場合である.

ここでは基本的には前者を想定する. 後者については「6. 補間法」で説明予定のスプライン補間が参考になる.

基本的な積分公式

次の積分を数値積分により求める.
\[
I = \int_{0.2}^1 e^x cos(x) dx = 1.1582
\] 最も簡単な方法は, 積分範囲の下端と上端を結ぶ直線で囲まれる台形の領域で近似する方法である. 計算式は次のとおりである.
\[
I = \frac{b – a}{2}(f(a) + f(b))
\] これは台形則とよばれる.

次に, もっと精度をよくするために \(a\) と \(b\) に加え中点 \(x_1 = (a + b)/2\) の 3 点を使って, 積分範囲の下端と上端を結び中点を通る放物線で囲まれる領域で近似する方法が考えられる. 計算式は次のとおりである.
\[
I = \frac{b – a}{6}(f(a) + 4f(x_1) + f(b))
\] これはシンプソン則とよばれる.

シンプソン則よりもさらに精度の良い方法としてガウス公式 (ガウス・ルジャンドル公式ともいう) がある. これは a と b の間のルジャンドル多項式の N 個のゼロ点を標本点として補間して近似する方法で, 計算式は次のとおりである.
\[
I = \frac{b – a}{2} \sum_{i=1}^N w_i f(\frac{a + b}{2} + \frac{b – a}{2} x_i)
\] これは (N 点) ガウス公式とよばれる. \(w_i\) は各分点における重みである.

以上, 3 つの方法を図で示すと次のようになる. ただし, ガウス公式は N = 2 の場合を示した.

これらは, 積分区間 [a, b] を分割する分点列 \(a = x_0 < x_1 < \dots < x_n = b\) または \(a < x_1 < x_2 < \dots < x_N < b\) において関数値 \(f(x_i)\) をとる補間多項式を定め, その積分値を求める方法である. これを補間型積分公式という.

次に, 積分区間を複数の小区間に分割し, その小区間ごとに上の公式を適用するものを複合公式という. 例えば, 積分区間を下図のように 4 分割しそれぞれに台形則を適用すれば単体の台形則より精度が高いことが期待できる.

これを複合台形則というが, 台形則の場合には実用上複合公式を使うのが普通であるためこれを単に台形則ということが多い. シンプソン則も同様に複合シンプソン則を指すことがある.

m 分割したときの複合台形則の計算式は次のようになる.
\[
I_m = h(\frac{f(x_0)}{2} + f(x_1) + f(x_2) + \dots + \frac{f(x_m)}{2}), \space h = \frac{b – a}{m}
\] h は小区間の幅である. ここでは h を一定としているが小区間ごとに可変にする方法も考えられる.

これらの公式を使って実際に最初にあげた例題を計算してみると次のようになった.

公式 計算値 S 誤差 計算量の目安
台形則 1.0663 -7.9% 2
シンプソン則 1.1575 -0.06% 3
ガウス公式 (N = 2) 1.15862 0.04% 2
複合台形則 (m = 4) 1.1523 -0.5% 5

計算量が 3 のシンプソン公式や計算量が 5 の複合台形則よりも, 計算量が 2 のガウス公式のほうが精度がよい.

数値積分を効率よく行うために, これらの公式と種々の技法を組み合わせて使うことが多い. 例えば, 台形則と加速法を組み合わせるロンバーグ積分, 台形則と変数変換を組み合わせる変数変換型積分公式, 高次の公式と複合公式 (区間分割) を組み合わせる適応自動積分などがある.

目次

9.1 補間型積分公式, ロンバーグ積分
9.2 変数変換型積分公式, 適応自動積分
9.3 実用プログラムのベンチマーク


本章の参考文献

[1] Forsythe他「計算機のための数値計算法」科学技術出版社 (1978)
[2] 森正武「FORTRAN77数値計算プログラミング(増補版)」岩波書店 (1987)
[3] 森口繁一「数値計算工学」岩波書店 (1989)
[4] 杉原正顕、室田一雄「数値計算法の数理」岩波書店 (1994)
[5] 森正武「数値解析 (第2版)」(2002) 共立出版
[6] Piessens, de Doncker-Kapenga, Uberhuber, Kahaner, “QUADPACK: A subroutine package for automatic integration”, Springer-Verlag. (1983)
[7] Ooura and Mori, “The double exponential formula for oscillatory functions over the half infinite interval”, J. Comput. Appl. Math., 38 (1991), 353-360.
[8] D. P. Laurie, “Calculation of Gauss-Kronrod quadrature rules”, Mathematics of Computation, Vol.66 No.219 (1997).
[9] Ooura and Mori, “A robust double exponential formula for Fourier-type integrals”, J. Comput. Appl. Math., 112 (1999), 229-241.