11. 数値積分

11.1 概要

\(f(x)\) は区間 \([a, b]\) 上で定義された関数とし, その定積分
\[
I = \int_a^b f(x) dx
\] の近似計算法を数値積分法 (numerical intergation) または求積法 (quadrature) といいます.

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

11.2 基本的な積分公式

次の積分 I を数値積分により求める例を使って説明します.
\[
\begin{align}
& f(x) = e^x cos(x) \\
& I = \int_{0.2}^1 f(x) dx = 1.15816896857196 \\
\end{align}
\]

最も簡単な方法は, 積分範囲の下端と上端を結ぶ直線で囲まれる台形の領域で近似する方法です. 計算式は次のようになります.
\[
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_N = \frac{b – a}{2} \sum_{i=1}^N w_i f(\frac{a + b}{2} + \frac{b – a}{2} x_i)
\] この公式の \(x_i\) と \(w_i\) は教科書や公式集の数表に載っています. n = 2 と 3 の場合には次のようになります. それぞれ, 2 点ガウス公式, 3 点ガウス公式とよびます.
\[
\begin{array}{ccccccc}
n & x_1 & x_2 & x_3 & w_1 & w_2 & w_3 \\ \hline
2 & -\sqrt{1/3} & \sqrt{1/3} & – & 1 & 1 & – \\
3 & -\sqrt{3/5} & 0 & \sqrt{3/5} & 5/9 & 8/9 & 5/9 \\
\hline
\end{array}
\] \(x_i\) には積分区間の両端 \(a\) と \(b\) が含まれません.

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

次に, 積分区間を複数の小区間に分割し, その小区間ごとに公式を適用すれば単体の公式より高精度が期待できそうです. これを, 複合公式といいます. 例えば, 下図のように積分区間を 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 を一定としていますが小区間ごとに可変にする方法もあります.

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

公式	計算値	誤差	計算量の目安
台形則	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 のガウス公式のほうが精度がよいことがわかります.

11.3 適応求積プログラム

得られた積分値の誤差が推定できたとすると, 必要な精度を満たすように部分区間の数を増減させるなどして自動的に計算量を最適化 (要求精度を満たすが計算量をできるだけ少なくすること) できます. このようなプログラムを適応求積プログラムといいます.

誤差を推定するためには, 例えば部分区間の数を変えて計算し結果を比較する方法などが考えられます. 台形公式やシンプソン公式で等間隔の部分区間に分割することした場合, 分割数が倍 (部分区間の幅が半分) になったときに半分の点では \(f(x)\) の計算済の結果を再利用できるため計算量を減らすことができます.

これに対して, ガウス公式の場合には \(n\) を変えると全部の点で \(f(x)\) の計算をし直す必要があり, 計算効率が良くありません. そこで, 1 つの公式で \(n\) が異なる 2 つの値を得ることができるガウス・クロンロッド公式が提案されました. これは次のようなものです.
\[
I = \frac{b – a}{2}(\sum_{i=1}^n a_i f(\frac{a + b}{2} + \frac{b – a}{2}x_i) + \sum_{j=1}^{n+1} b_j f(\frac{a + b}{2} + \frac{b – a}{2}y_j))
\] ここで, \(a_i\) と \(b_j\) はガウス・クロンロッド公式の係数です. 始点・終点を表す \(a, b\) とは別のものです.

\(x_i\) はガウス公式と同じものです. したがって, ガウス公式の \(w_i\) と組み合わせると \(n\) 点ガウス公式を使った積分値 (\(G_n\) と表す) を求めることができます. なお, \(a_i\) と \(w_i\) の値は異なります.

全体では \(2n + 1\) 点ガウス・クロンロッド公式を使った積分値 (\(K_{2n+1}\) と表す) を求めることができます.

\(K_{2n+1}\) は \(f(x)\) の計算回数が同じ \(G_{2n+1}\) よりは精度が悪くなります. つまり, ガウス・クロンロッド公式は, ガウス公式よりはやや精度が落ちる代わりに誤差推定機能を追加したものといえます.

\(n = 2\) の場合 (5 点ガウス・クロンロッド公式) の係数は次のようになります.
\[
\begin{array}{cccc}
\hline
x_i & -0.5773502691896258 & 0.5773502691896258 \\
a_i & 0.4909090909090909 & 0.4909090909090909 \\
y_j & -0.9258200997725515 & 0 & 0.9258200997725515 \\
b_j & 0.1979797979797980 & 0.6222222222222222 & 0.1979797979797980 \\
\hline
\end{array}
\] 求められた \(G_n\) と \(K_{2n+1}\) から誤差を推定することができます.

現在では数値積分には精度をある程度保証できるように適応求積プログラムを使うのが一般的になっています. 特に, 各部分区間にガウス・クロンロッド公式を適用して誤差を推定し, 各部分区間の幅を最適化して精度を保ったまま計算量を減らすようにしたプログラムがよく使われています.

11.4 データ入力タイプの数値積分プログラム

被積分関数の値がとびとびの点におけるデータとしてしか与えられていない場合には, 上記の積分公式が使えない場合がほとんどです. ただし, 等間隔に十分なデータが与えられれば複合台形則やシンプソン則を使えそうです.

一般的には, 入力されたデータ点の補間を行い, その補間式の積分を計算することにより積分値を推定します.

多項式補間するものやスプライン補間するものなど種々のプログラムがあります. ただし, 入力データにもよりますが, 一般的には高い精度が期待できるものではありません.

11.5 XLPack を使った数値積分の求め方

VBA サブルーチン Qk15 は分点固定 (15点) のガウス・クロンロッド公式による数値積分プログラムです. 15 回の関数呼び出しで積分値とその推定誤差を計算します. VBA サブルーチン Qag はガウス・クロンロッド公式による適応求積プログラムです. 被積分関数が滑らかな場合には Qk15 で十分な精度が得られます. 被積分関数が滑らかでない場合や, 関数呼び出しが増えてでも所定の目的精度を得たい場合には, Qag を使うとよいでしょう.

積分範囲の片方あるいは両方が有限でない場合には VBA サブルーチン Qagi を使うことができます. これもガウス・クロンロッド公式による適応求積分プログラムです. 半無限積分を有限区間 [0, 1] の積分に変換して計算します.

Qag および Qagi は XLPack ソルバーから使うこともできます.

3 次スプライン補間係数を求める VBA サブルーチン Pchse あるいはワークシート関数 WPchse (「5. 補間法」参照) に VBA サブルーチン Pchia あるいはワークシート関数 WPchia を組み合わせることによりデータ入力タイプの数値積分プログラムとして使用することができます.

以下, 次の積分を数値積分により求める例を使って説明します.
\[
\begin{align}
I = \int_0^4 \frac{1}{x + 1} dx = 1.6094379124341
\\
\end{align}
\]

11.3.1 VBA プログラムを使用した解き方 (関数入力) (1)

Qag および Qk15 を使ったプログラム例を示します.

Function F(X As Double) As Double
    F = 1 / (1 + X)
End Function

Sub Start()
    Dim A As Double, B As Double, S As Double, Info As Long
    Dim AbsErr As Double, Neval As Long
    '--- Input data
    A = Cells(4, 2): B = Cells(5, 2)
    '--- Compute integration by Qag
    Call Qag(AddressOf F, A, B, S, Info, AbsErr, Neval)
    Cells(11, 2) = Info
    If Info = 0 Then
        Cells(8, 2) = S
        Cells(9, 2) = AbsErr
        Cells(10, 2) = Neval
    End If
    '--- Compute integration by Qk15
    Call Qk15(AddressOf F, A, B, S, AbsErr)
    Cells(8, 3) = S
    Cells(9, 3) = AbsErr
End Sub

目的関数を外部関数として用意し, 区間と必要に応じて他のパラメータを指定して呼び出します.

このプログラムを実行すると, 次の結果が得られます.

Qag は, 被積分関数呼び出し 75 回 (内部で Qk15 を 5 回呼び出し) で推定誤差 \(7.38 \times 10^{-13}\) の結果が得られました. Qk15 は被積分関数呼び出し 15 回で推定誤差 \(2.07 \times 10^{-5}\) となりました. 推定誤差は安全方向に算出されており, 実際にはそれぞれ \(2.22 \times 10^{-16}, 1.20 \times 10^{-11}\) でした.

11.3.2 VBA プログラムを使用した解き方 (関数入力) (2)

リバースコミュニケーション版 (RCI) の Qag_r および Qk15_r を使ったプログラム例を示します.

Sub Start()
    Dim A As Double, B As Double, S As Double, Info As Long
    Dim AbsErr As Double, Neval As Long
    Dim XX As Double, YY As Double, IRev As Long
    '--- Input data
    A = Cells(4, 2): B = Cells(5, 2)
    '--- Compute integration by Qag_r
    IRev = 0
    Do
        Call Qag_r(A, B, S, Info, XX, YY, IRev, AbsErr, Neval)
        If IRev = 1 Then YY = 1 / (1 + XX ^ 2)
    Loop While IRev <> 0
    Cells(11, 2) = Info
    If Info = 0 Then
        Cells(8, 2) = S
        Cells(9, 2) = AbsErr
        Cells(10, 2) = Neval
    End If
    '--- Compute integration by Qk15_r
    IRev = 0
    Do
        Call Qk15_r(A, B, S, XX, YY, IRev, AbsErr)
        If IRev = 1 Then YY = 1 / (1 + XX ^ 2)
    Loop While IRev <> 0
    Cells(8, 3) = S
    Cells(9, 3) = AbsErr
End Sub

目的関数を外部関数として与えるのではなく, IRev = 1 のときに XX の値を使って関数値を計算し YY に入れて再度 Qag_r あるいは Qk15_r を呼び出します. RCI の詳細についてはこちらを参照してください.

このプログラムを実行すると上と同じ結果が得られます.

11.3.3 ワークシート関数を使用した解き方 (データ入力)

データ入力の数値積分ルーチンとして, スプライン補間による積分の求め方を説明します.

上と同じ例において, 関数形は不明で, \(x = 0, 0.5, 1, \dots, 4\) における関数値のみが与えられているものとします.
\[
\begin{array}{cccc}
\hline
x & f(x) \\ \hline
0 & 1 \\
0.5 & 0.666667 \\
1 & 0.5 \\
1.5 & 0.4 \\
2 & 0.333333 \\
2.5 & 0.285714 \\
3 & 0.25 \\
3.5 & 0.222222 \\
4 & 0.2 \\
\hline
\end{array}
\] まずは, WPchse 関数を使ってスプライン関数の係数を求めます(「5. 補間」を参照).

いったん係数が求められれば, WPchia 関数を使って積分値を計算することができます. WPchia の必要なパラメータは A, B, N, X, Y, D です. A と B は積分範囲(この例では 0 と 4). N はデータ数(この例では 9), X と Y は補間データのセル範囲です. D はスプライン係数のセル範囲です.

Enter を押すと, 1.61082… で, 9点のデータから3桁程度の精度で積分値を求めることができました.

注 – 例題ワークシートには VBA を使ったスプライン補間による積分例も入っています.

11.3.4 ソルバーを使用した解き方

XLPack ソルバーアドインの「数値積分 (有限区間)」を使って解くこともできます. B4セルに被積分関数式 (=1/(1+B3^2)) が入力されています.

ソルバーについてはこちらも参照ください.

公開日: 2025/03/31 (更新日: 2026/07/19)