シンプソンの公式で積分近似計算を行います

 

シンプソンの公式 Simpson's rule

 ある曲線上の 3 点
 
\[(\alpha,\:f_0),\quad \left(\frac{\alpha+\beta}{2},\:f_1\right),\quad (\beta,\:f_2)\]
を通る放物線と $x=a,\:x=b$ および $x$ 軸とで囲まれる面積を求めてみます。

 VBA曲線を放物線で近似

放物線を $y=Ax^2+Bx+C$ とおいて面積を計算すると
 
\[s=\int_{\alpha}^{\beta}(Ax^2+Bx+C)dx=\frac{A}{3}(\beta^3-\alpha^3)+\frac{B}{2}(\beta^2-\alpha^2)+C(\beta-\alpha)\]
 $\beta^3-\alpha^3$ を因数分解すると
 
\[\beta^3-\alpha^3=(\beta-\alpha)(\alpha^2+\alpha\beta+\beta^2)\]
となるので、式を整理すると
 
\[s=\frac{\beta-\alpha}{3}\left\{ A(\alpha^2+\alpha\ \beta+\beta^2)+\frac{3B}{2}(\alpha+\beta)+3C \right\}\]
となります。放物線は
 
\[(\alpha,\:f_0),\quad \left(\frac{\alpha+\beta}{2},\:f_1\right),\quad (\beta,\:f_2)\]
を通るので、
 
\[\begin{align*}f_0&=A\alpha^2+B\alpha+C\\[6pt]
f_1&=A\left( \frac{\alpha+\beta}{2} \right) ^2+B \left( \frac{\alpha+\beta}{2} \right) +C\\[6pt]
f_2&=A\beta^2+B\beta+C\end{align*}\]
が成り立ちます。そこで $f_0+4f_1+f_2$ を計算してみると
 
\[f_0+4f_1+f_2=2\left\{ A(\alpha^2+\alpha\ \beta+\beta^2)+\frac{3B}{2}(\alpha+\beta)+3C\right\}\]
となっていることがわかります。よって面積は
 
\[s=\frac{\beta-\alpha}{6}(f_0+4f_1+f_2)\]
と計算できます。今度は下図のように曲線の区間 $[a,\:b]$ の面積を $2m$ 個の帯に分割して曲線上の点 $P_0,\:P_1,\:P_2,\:\cdots,\:P_{2m-2},\:P_{2m-1},\:P_{2m}$ とします。

 VBAシンプソン積分公式

 1 区分の幅は
 
\[h=\frac{b-a}{2m}\]
で与えられます。そして帯を 2 個ずつまとめて
 
\[P_0P_1P_2,\:P_2P_3P_4,\:P_4P_5P_6,\:\cdots\]
という弧をもつ $m$ 個の帯をつくります。区間 $[a,\:b]$ において、曲線 $y=f(x)$ と $x=a,\:x=b$ および $x$ 軸に囲まれる面積を
 
\[S=\int_{a}^{b}S(x)dx\]
とおいて、それぞれの帯の面積に先ほどの放物線近似を適用すると

\[S=\frac{h}{3}\sum_{k=0}^{m-1}\left[f(a+2kh)+4f(a+(2k+1)h)+f(a+(2k+2)h)\right]\]

となります。これが シンプソンの公式 とよばれる積分近似式です。
 

シンプソンの公式を使った積分

 シンプソンの公式を使って
 
\[\int_{a}^{b}\exp (-x^2)dx\]
を計算する VBA マクロを作ってみます。

 Sub シンプソン公式()

 Dim a As Double, b As Double
 Dim h As Double, s As Double
 Dim m As Integer, k As Integer

 区間[a,b]の入力を促します
 a = Application.InputBox("左端を入力してください", Type:=1)
 b = Application.InputBox("右端を入力してください", Type:=1)

 If a > b Then
  MsgBox "数値が正しくありません"
 Exit Sub

 Else

  mを指定します
  m = 10

  1 区分の長さです
  h = (b - a) / (2 * m)

  帯を足し合わせます
  For k = 0 To m - 1
   s = s + exp(-(a + 2 * k * h) ^ 2) _
   + 4 * exp(-(a + (2 * k + 1) * h) ^ 2) _
   + exp(-(a + (2 * k + 2) * h) ^ 2)
  Next k

  s = s * h / 3

  MsgBox s

 End If

 End Sub

 積分を実行すると端点 $a$ と $b$ を入力するように促されますので、
 たとえば $a=0,\:b=1$ としてみると 0.74682 という値が返ってきます。

 ≫ VBA 数値計算

スポンサーリンク
スポンサーリンク
末尾広告
末尾広告

コメントをどうぞ

メールアドレスが公開されることはありません。

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください