VBAで数値積分　シンプソンの公式(シンプソン法)

16115 ワード
VBA 数値積分 数学 積分 VBA テキストリンク

1．はじめに

本記事では，数値積分の手法の1つであるシンプソンの公式による方法をVBAで実装してみました．
数値積分については他にもいろんな方が記事を書いているのでここでは説明を省略したいと思います．

2．シンプソンの公式

関数f(x)のaからbまでの積分はシンプソンの公式(または，シンプソン法)を使うと次式のように表されます．

\int _{a} ^{b} f(x)dx=\dfrac{b-a}{6} \left\{ f(a) + 4 f\left( \dfrac{a+b}{2}\right) + f(b) \right\}

高校数学の美しい物語：シンプソンの公式の証明と例題
https://mathtrain.jp/simpson

例えば，pの関数f(p)をp1からp7まで積分する場合を考えます．

\int _{p1} ^{p7} f(x)dx=\dfrac{p7-p1}{6} \left\{ f(p1) + 4 f\left( \dfrac{p1+p7}{2}\right) + f(p7) \right\}

これでもいいのですが，ざっくりしすぎな気もしますね．
数値積分などをするときは積分の計算結果をできるだけ正しく得たいため，積分区間を細かく分割することが多いです．
例えば，p1からp7の積分を6つに分けて考えてみましょう．
ちなみに図中の記号(x1，x2，D)は無視しておいてください．

\begin{align}
\int _{p1} ^{p7} f(x)dx&=\int _{p1} ^{p2} f(x)dx+\int _{p2} ^{p3} f(x)dx+\int _{p3} ^{p4} f(x)dx\\
&+\int _{p4} ^{p5} f(x)dx+\int _{p5} ^{p6} f(x)dx+\int _{p6} ^{p7} f(x)dx\\
\end{align}

これにシンプソン公式を当てはめると次式になります．

\begin{align}
\int _{p1} ^{p7} f(x)dx&=\dfrac{p2-p1}{6} \left\{ f(p1) + 4 f\left( \dfrac{p1+p2}{2}\right) + f(p2) \right\}\\
&+\dfrac{p3-p2}{6} \left\{ f(p2) + 4 f\left( \dfrac{p2+p3}{2}\right) + f(p3) \right\}\\
&+\dfrac{p4-p3}{6} \left\{ f(p3) + 4 f\left( \dfrac{p3+p4}{2}\right) + f(p4) \right\}\\
&+\dfrac{p5-p4}{6} \left\{ f(p4) + 4 f\left( \dfrac{p4+p5}{2}\right) + f(p5) \right\}\\
&+\dfrac{p6-p5}{6} \left\{ f(p5) + 4 f\left( \dfrac{p5+p6}{2}\right) + f(p6) \right\}\\
&+\dfrac{p7-p6}{6} \left\{ f(p6) + 4 f\left( \dfrac{p6+p7}{2}\right) + f(p7) \right\}\\
\end{align}

いま分割幅は一定とすると，

\dfrac{p2-p1}{6}=\dfrac{p3-p2}{6}=\dfrac{p4-p3}{6}=\dfrac{p5-p4}{6} =\dfrac{p6-p5}{6}=\dfrac{p7-p6}{6}

となる．ここで，

d=\dfrac{p2-p1}{6}

とおくと，

\begin{align}
\dfrac{1}{d} \int _{p1} ^{p7} f(x)dx&=\left( 4 f\left( \dfrac{p1+p2}{2}\right) + 2f(p2) \right)\\
&+ \left( 4 f\left( \dfrac{p2+p3}{2}\right) + 2f(p3) \right)\\
&+ \left( 4 f\left( \dfrac{p3+p4}{2}\right) + 2f(p4) \right)\\
&+ \left( 4 f\left( \dfrac{p4+p5}{2}\right) + 2f(p5) \right)\\
&+ \left( 4 f\left( \dfrac{p5+p6}{2}\right) + 2f(p6) \right)\\
&+ 4 f\left( \dfrac{p6+p7}{2}\right)\\
&+ f(p1)+f(p7) \\
\end{align}

このようなおかしな表記をしている理由は，プログラムで実際に計算するときの考え方に合わせているためです．
ここでは，この式を「式S」と呼ぶことにします．

では，実際にプログラムを見てみましょう．

3．プログラムを読み解く

プログラムは以下の通りです．
後ほど説明するので今は読み飛ばして構いません．

Option Explicit
Function FNE(ByRef θ As Double)
'引数(角度θ(rad))が与えられるとsinθの値を返す関数
FNE = Sin(θ)
End Function
Sub NumericalIntegration()
'最終更新日：2021/2/15
'シンプソンの公式を使ってsinθの積分を数値解析的に求めるプログラム
'積分区間[A1,B1]をND個に分割した小区間にそれぞれシンプソン公式を当てはめる

'*********************変数の定義*********************
Dim i As Integer    'カウンタ
Dim A1 As Double    '積分区間[A1,B1]の左端
Dim B1 As Double    '積分区間[A1,B1]の右端
Dim ND As Integer   '積分区間の分割数
Dim x1 As Double    'シンプソン公式第二項の独立変数
Dim x2 As Double    'シンプソン公式第三項の独立変数
Dim S As Double     '数値積分の計算結果
Dim D As Double     '分割幅の1/2

'*********************シンプソン公式による数値積分*********************
'積分区間を定義
A1 = 0      '単位：°
B1 = 180    '単位：°

'積分区間の単位を°からradに変換
A1 = 2 * 4 * Atn(1) / 360 * A1
B1 = 2 * 4 * Atn(1) / 360 * B1

'分割数をテキトーに設定
ND = 500

D = (B1 - A1) / (2 * ND)    '分割幅の1/2
x1 = A1 + D                 '当該小区間の中央に相当する独立変数の値
x2 = x1 + D                 '当該小区間の右端に相当する独立変数の値
S = 0                       '初期化
For i = 1 To ND - 1
    S = S + 4 * FNE(x1) + 2 * FNE(x2)
    x1 = x2 + D
    x2 = x1 + D
Next i
S = S + 4 * FNE(x1)
S = S + FNE(A1) + FNE(B1)
S = S * (D / 6 * 2)
Debug.Print "シンプソン公式による積分値："; S

End Sub

プログラムの核になるのは以下の部分です．
まずはこの部分と変数の定義のところだけを見てください．

D = (B1 - A1) / (2 * ND)    '分割幅の1/2
x1 = A1 + D                 '当該小区間の中央に相当する独立変数の値
x2 = x1 + D                 '当該小区間の右端に相当する独立変数の値
S = 0                       '初期化
For i = 1 To ND - 1
    S = S + 4 * FNE(x1) + 2 * FNE(x2)
    x1 = x2 + D
    x2 = x1 + D
Next i
S = S + 4 * FNE(x1)
S = S + FNE(A1) + FNE(B1)
S = S * (D / 6 * 2)

はじめに，小区間p1～p2について考えます．
分割幅の1/2を表す変数D，当該小区間の中央に相当する独立変数の値x1，当該小区間の右端に相当する独立変数の値x2を図に表すと以下の図のようになります．
なお，積分区間[A1,B1]=[p1,p7]です．分割数NDは6となります．

つまり，

x1=\dfrac{p2-p1}{2}\\
x2=p2\\

for分では，下図のようにx1とx2を1小区間ずつずらしながら式Sの第1項～第5項に相当する計算をします．
図はfor文のカウンタiが2のときを表しています．

そしてfor分の次の文は式Sの第6項を，さらに次の文は式Sの第7項，第8項を表しています．

最後に次に示す文が，式Sの両辺にdをかけることを表しておりプログラム中の変数Sが積分値になります．

S = S * (D / 6 * 2)

4．プログラムの実行結果についての考察

さて，実は今回紹介したプログラムはsin θをθ =0°～180°まで積分するプログラムになっています．
分割数はテキトーに500としました．
ちなみに，解析解は2になります．

\int _{0} ^{\pi} sin \theta \cdot d\theta=2

プログラムの計算結果は2.00000000000111ぐらいになると思います．
ちなみにシンプソン公式をそのまま適用すると(小区間に分割しない)，計算結果は2π/3になります．

\begin{align}
\int _{0} ^{\pi} sin \theta d \theta &=\dfrac{\pi-0}{6} \left\{ f(0) + 4 f\left( \dfrac{0+\pi}{2}\right) + f(\pi) \right\}\\
&=\dfrac{\pi}{6} \left\{ sin(0) + 4 sin\left( \dfrac{\pi}{2}\right) + sin(\pi) \right\}\\
&=\dfrac{\pi}{6} \left( 0+ 4 \cdot  1 + 0 \right)\\
&=\dfrac{2}{3}\pi\\
&\fallingdotseq 2.094

\end{align}

分割することでより正確になることが確認できますね．
ただし，分割数を上げれば必ずしもいい結果が得られるとは限らないので注意してください(分割数を上げすぎると計算過程で数値を丸めるときの誤差によりむしろ精度が悪くなることがある)．

5．おわりに

最後までお付き合いいただきありがとうございました．

今回はいつもと趣向が違う記事でしたが，今後も農業農村工学(水文学，かんがい排水，土壌物理学，水理学)を中心に記事を執筆していきたいと思います．
リクエスト等も受け付けておりますので，ご遠慮なく連絡ください．
Twitter Acount：https://twitter.com/6LxAi9GCOmRigUI
メール：[email protected]
ツイッターで絡んでくれると喜びます(*´･ω･｡)σｶﾏﾁｮ

引用・参考文献

(1) 河村三郎：土木工学プログラム集1　水文・水理1　実用プログラムとグラフィック表示，森北出版，1984，pp.6-11．
(2) Keiji Harada：「数値積分法」
https://www-np.acs.i.kyoto-u.ac.jp/~harada/education/open_courses/integral/index.html
(3) 重田出：数値積分と数値微分（基礎）
https://www.sci.kagoshima-u.ac.jp/shigeta/AC-NumericalDiffAndInt-Basic.pdf