11. 数値積分

注 – 例題ワークシートをダウンロードすることができます. 1. 例題ワークシートの入手と使用方法をご覧ください.

数値積分の XLPack による計算例を説明します.

例題

次の定積分を求める.

  
  ∫ 1/(1 + x²) dx [0, 4]

関数の形は次のとおりで, 積分の真値は arctan(4) になります.

11-01

数値積分の汎用ルーチンは, 積分区間内の任意の点において関数値が計算できる関数入力タイプと, 積分区間を含む有限個の点の関数値が与えられるデータ入力タイプがあります.

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

VBA サブルーチン Qag および Qk15 を使ったプログラム例を示します.

Qag は適応求積ルーチンとよばれ, 要求精度を与えるとそれを満たすように積分範囲の部分区間を調節して積分則を適用し積分値を求めます. Qk15 は固定分点数(15点)の求積ルーチンで, 15回の関数評価を行い積分値を求めるとともに推定誤差を返します.

Function F(X As Double) As Double
    F = 1 / (1 + X ^ 2)
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 は, 被積分関数呼び出し105回で15桁いっぱいで正しい結果が得られました. Qk15 は被積分関数呼び出し15回で 0.00148 の推定誤差と返していますが, 実際には9桁まで正しい値が求められています.

要求精度を満たす結果を得たい場合には適応求積ルーチン, 計算速度を優先しある程度の精度でよい場合には固定分点求積ルーチンが適しています.

積分区間の片方, あるいは, 両方が(正あるいは負の)無限大の場合には, 適応求積ルーチン Qagi を使うことができます.

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

リバースコミュニケーション版の VBA サブルーチン 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 を呼び出します. リバースコミュニケーション版の詳細についてはこちらを参照してください.

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

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

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

上と同じ例において, 関数形は不明で, x = 0, 0.5, 1, …, 4 における関数値のみが与えられているものとします.

x	f(x)
0	1
0.5	0.8
1	0.5
1.5	0.3076923
2	0.2
2.5	0.137931
3	0.1
3.5	0.0754717
4	0.0588235

まずは, WPchse 関数を使ってスプライン関数の係数を求めます(「5. 補間」を参照).

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