シンプソンの公式で定積分するExcel VBA

夜中に突然、定積分の近似値が必要になることがありませんか？
え？ない？
僕も最近はありません。

でも、約10年前には必要とすることがあったようです。
たまたま、昔のExcelファイルからシンプソンの公式で積分するExcelのマクロが発掘され、2008年の日付が付いていました。

コードを読んでみると、かなり稚拙だったので、全面的に書き直しました。

もし、突然に定積分の値が必要になったときにでも、お使いいただければと思います。

シンプソンの公式マクロの使い方
1. １）マクロをコピペする
2. ２）「Simpson」の初回実行
３）所望の関数に書き換える
1. ４）「Simpson」の再実行
2. ５）グラフを書く

シンプソンの公式マクロの使い方

１）マクロをコピペする

Visual Basicのエディタを開き、マクロをコピペして下さい。
マクロの使い方が分からない方は、こちらをご覧ください。

'-----------------------------------------------------------------------
' Summary : シンプソン法により、定積分を計算する
' Comments: 初めは単純に本関数を呼び出す
'         : 次からは、関数、下限値、上限値、分割数を変えて呼び出す
'         : 他のブックからの呼び出し方法:
'         : Application.Run "Excelのファイル名!Simpson"
'         : 例：個人用マクロブックにこのマクロを入れた場合は下行になる
'         : Application.Run "PERSONAL.XLSB!Simpson"
'-----------------------------------------------------------------------
Public Sub Simpson()
    Const cstOffset     As Integer = 2  ' シート記入上のオフセット値
    Const cstSheetName  As String = "シンプソン"    ' 作業用シート名
    Dim dbla            As Double       ' 下限、
    Dim dblb            As Double       ' 上限
    Dim dblVal          As Double       ' 積分結果
    Dim i               As Integer      ' ループ制御変数
    Dim intN            As Integer      ' 分割数（偶数）
    Dim sht             As Worksheet    ' シート制御変数

'   作業用シートの確認
    For Each sht In Worksheets
        If StrComp(sht.Name, cstSheetName, vbTextCompare) = 0 Then
            Exit For
        End If
    Next sht
    
'   もし、作業用シートが無いならシートを作る
    If sht Is Nothing Then
'   シートを追加して、シート名を設定
        Set sht = Worksheets.Add
        sht.Name = cstSheetName
        With sht
            ' タイトルを付ける
            .Range("A1").Value = "関数 f(x)"
            .Range("A2").Value = "変数 x"
            .Range("A3").Value = "下限 a"
            .Range("A4").Value = "上限 b"
            .Range("A5").Value = "分割数 N"
            .Range("A6").Value = "積分結果"
            
            ' B列のセルに名前を付ける
            .Names.Add Name:="fx", RefersTo:="='" & .Name & "'!$B$1"
            .Names.Add Name:="x", RefersTo:="='" & .Name & "'!$B$2"
            .Names.Add Name:="a", RefersTo:="='" & .Name & "'!$B$3"
            .Names.Add Name:="b", RefersTo:="='" & .Name & "'!$B$4"
            .Names.Add Name:="N", RefersTo:="='" & .Name & "'!$B$5"
            .Names.Add Name:="Result", RefersTo:="='" & .Name & "'!$B$6"
                
            ' 関数や上下限を入れる
            .Range("fx").Value = "=exp(-(x^2)/2)/sqrt(2*pi())"
            .Range("x").Value = "0"
            .Range("a").Value = "-3"
            .Range("b").Value = "3"
            .Range("N").Value = "100"
        End With
    End If
    
    ' シンプソン法による積分
    ' 作業用シートにて作業
    With Worksheets(cstSheetName)
        ' 作業用のD、E列をクリア
        .Range("D:E").Clear
        
        ' 各変数に値をセット
        dbla = .Range("a").Value
        dblb = .Range("b").Value
        intN = .Range("N").Value
        ' 分割数を偶数に丸める
        intN = ((intN + 1) \ 2) * 2
        Range("N").Value = intN
        
        ' D列、E列の１行目に、変数と関数のタイトルを示す
        .Range("D1") = .Range("A2")
        .Range("E1") = .Range("A1")
        
        ' 変数xに下限値をセット
        .Range("x").Value = dbla
        ' 関数値を計算
        Calculate
        ' 積分の初期値を内部変数dblValに保持
        dblVal = Range("fx").Value
        ' xの下限値とその時の関数値をD2とE2にそれぞれ記入
        .Range("D" & cstOffset) = .Range("x")
        .Range("E" & cstOffset) = .Range("fx")
        
        ' ループしながらシンプソン法を実行
        For i = 1 To intN - 1
            ' xの次の値をセット
            .Range("x").Value = dbla + (dblb - dbla) * i / intN
            ' 関数値を計算
            Calculate
            ' 制御変数が奇数か偶数かで計算方法が異なる
            If (i Mod 2 = 1) Then
                ' 奇数の場合は関数値を4倍して積分値に加える
                dblVal = dblVal + 4 * .Range("fx").Value
            Else
                ' 偶数の場合は関数値を2倍して積分値に加える
                dblVal = dblVal + 2 * .Range("fx").Value
            End If
            ' 変数xと関数f(x)の値をそれぞれD、E列に出力する
            .Range("D" & (i + cstOffset)) = .Range("x")
            .Range("E" & (i + cstOffset)) = .Range("fx")
        Next i
        ' 変数xに上限値をセット
        .Range("x").Value = dblb
        ' 関数値を計算
        Calculate
        ' 最終的な積分値を計算する
        dblVal = (dblVal + .Range("fx").Value) * (dblb - dbla) / (3 * intN)
        ' 変数xと関数f(x)の最終値をそれぞれD、E列に出力する
        .Range("D" & (intN + cstOffset)) = .Range("x")
        .Range("E" & (intN + cstOffset)) = .Range("fx")
    
        ' 積分の結果を表示する
        .Range("Result").Value = dblVal
    End With
End Sub

２）「Simpson」の初回実行

「Simpson」というマクロを実行して下さい。

「シンプソン」というシートができて、そこに勝手に値が書きこまれていきます。３秒以内で終わると思います。

Simpsonの初回実行結果

このシートの内容は、正規化した正規分布の確率密度関数の-3σ～3σを積分した結果、つまり次式の結果を示しています。

$\int_a^b f(x) dx=\int_{-3}^{3} \dfrac{1}{\sqrt{2\pi}} e^{\frac{-x^{2}}{2}} dx$

ここで、

B1はf(x)の定義式（=EXP(-(x^2)/2)/SQRT(2*PI())）です。
B2は確率変数xの値で、マクロ終了後には積分範囲の上限値が入っています。
B3は積分範囲の下限a（-3）です。
B4は積分範囲の上限b（3）です。
B5は分割数N（100）です。
B6は積分の結果（厳密には0.997300192454353）です。

ここで、Excelが正解と考えている「=1-(1-NORMSDIST(3))*2」の結果を計算してもらうと「0.99730020393674」なので、差は「-1.14824E-08」です。どちらが正解に近いのか分かりませんが、なかなかの高精度で合っているので、マクロのコードに致命的なミスはなさそうです。