11. 数値積分

数値積分のXLPackによる計算例を示します。

例題

次の定積分を求める。

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

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

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

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(3, 2): B = Cells(4, 2)
    '--- Compute integration by Qag
    Call Qag(AddressOf F, A, B, S, Info, AbsErr, Neval)
    Cells(10, 2) = Info
    If Info = 0 Then
        Cells(7, 2) = S
        Cells(8, 2) = AbsErr
        Cells(9, 2) = Neval
    End If
    '--- Compute integration by Qk15
    Call Qk15(AddressOf F, A, B, S, AbsErr)
    Cells(7, 3) = S
    Cells(8, 3) = AbsErr
End Sub

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

このプログラムを実行すると、次の結果が得られました。

Qagは、被積分関数呼び出し105回で15桁いっぱいで正しい結果が得られました。Qk15は被積分関数呼び出し15回で0.00148の推定誤差と返していますが、実際には9桁まで正しい値が求められています。

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

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

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

リバースコミュニケーション版のVBAサブルーチンQag_rおよびQk15_rを使ったプログラム例を示します。

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
    Dim XX As Double, YY As Double, IRev As Long
    '--- Input data
    A = Cells(3, 2): B = Cells(4, 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 = F(XX)
    Loop While IRev <> 0
    Cells(10, 2) = Info
    If Info = 0 Then
        Cells(7, 2) = S
        Cells(8, 2) = AbsErr
        Cells(9, 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 = F(XX)
    Loop While IRev <> 0
    Cells(7, 3) = S
    Cells(8, 3) = AbsErr
End Sub

目的関数を外部関数として与えるのではなく、IRev = 1のときにXXの値を使って関数値を計算しYYに入れて再度Qag_rあるいはQk15_rを呼び出します。リバースコミュニケーション版の詳細についてはこちらを参照してください。

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