みなさんこんばんにちは、ユウです。
今回はプログラミングで数値積分する方法の一つ、シンプソン公式について紹介していきたいと思うよ！
このシンプソン公式は、前回紹介したニュートン・コーツ公式を応用して使っていくのでぜひそちらの記事で、ニュートン・コーツ公式を理解した上でこちらの記事を見ることをお勧めします。
ニュートン・コーツ公式（おさらい）
冒頭でも述べた通り、今回はニュートン・コーツ公式について簡単に紹介していきたいと思います。
詳しく知りたい方は、以前の記事を見てください。
ニュートン・コーツ公式は以下のように表されます。
\[\int_{a}^{b} P_{n} (x) \, dx = \sum_{k=0}^{n} f_{k} \int_{a}^{b} l_{k} (x) \, dx = \sum_{k=0}^{n} a_{k} f_{k} \]
ただし、
\[a = \int_{a}^{b} l_{k} (x) \, dx \]
ニュートン・コーツ公式はラグランジュの補間多項式を利用して数値積分を行っています。
ラグランジュの補間多項式は元の曲線を多項式で近似するので、その近似した関数を用いて面積を求めていきます。
ニュートン・コーツ公式の特徴として、関数f(x)のx軸を区切る分点を以下のように取ります。
\[a=x_{0}<x_{1}<x_{2}<\dots<x_{n}=b\]
この際、分点の間隔はaとbの差をある大きい値nで割ったものとなります。
このnの大きさは大きければ大きいほど、近似した関数との誤差が少ないので精度がよくなります。（ただし、処理速度を考慮すると大きければ良いというわけではありません）
このaとbの差をnで割ったものをhとすると、数式では以下のようになります。
\[h=\frac{b-a}{n}\]
これから、hを使ってxの値等を表現していくと次の通りです。
\[x_{0}=a, x_{1}=a+h, x_{2}=a+2h, \dots ,x_{n}=a+nh=b \]
今回紹介する、シンプソン公式はこれらを応用して導出していくので、ここをきちんと理解していてください。
シンプソン公式
シンプソン公式を理解するには、ラグランジュの補間公式を理解しておく必要があります。
ですので、まずはラグランジュの補間公式について少し紹介していきたいと思います。
ラグランジュの補間公式
ラグランジュの補間公式は以下のように表すことができます。
ラグランジュの補間公式：
ｘ座標が相異なるn+1点(x1,y1),(x2,y2),・・・,(xn+1,yn+1)を通るn次以下の多項式P(ｘ)が１つ定まり、以下の式で表される：
\[P(x) = \sum_{i=1}^{n+1} y_{i} \frac{ f_{i} (x)}{f_{i} (x_{i})} \]
ただし、
\[f_{i} (x) = \prod_{k \neq i} (x – x_{k}) \]
この公式により、関数を簡略化でき計算しやすくできるようにしていく。
ただ、いきなり公式を出されても少しわかりにくいと感じるかと思うので、少し計算をしてみようと思います。
ここでは、後でわかりやすくするためにn=2とした場合のラグランジュ補間をやっていきたいと思います。
\[P(x) = y_{0} \frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})}\]
\[+ y_{1} \frac{(x-x_{0})(x-x_{2})}{(x_{1}-x_{0})(x_{1}-x_{2})}\]
\[+ y_{2} \frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})}\]
この場合、右辺はn=2次多項式になっていることがわかります。
まず、この式がどのように上の式を使っているのかわからない方はまずは、上の式を頭の中で計算するか、紙に起こして計算してみてください。
ちなみに下の記号は、「相乗」の意味なのでかけ算を行っています。
\[\prod \]
シンプソン公式導出
それでは、シンプソン公式について紹介していきたいと思います。
シンプソン公式は先ほど紹介したラグランジュの補間公式を使用して考えていくので、ぜひ覚えるようにしてください。
最初に紹介したニュートン・コーツ公式において、n=2とするとhは以下のようになります。
\[h= \frac{b-a}{2} \]
さらに、x0などの値もhから以下のようになります。
\[x_{0} = a, x_{1} = a+h, x_{2} = a+2h = b \]
これからから、xをx=a+shとおくと、ラグランジュの補間公式から関数を求め、その関数を定積分すると、以下のようになる。
\[\int_{a}^{b} \frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})} + \frac{(x-x_{0})(x-x_{1})}{(x_{1}-x_{0})(x_{1}-x_{2})} + \frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})} \, dx\]
これを、分割して考えていくと以下のようになる。
\[\alpha_{0} = \int_{a}^{b} \frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})} ,\ dx\]
\[\alpha_{1} = \int_{a}^{b} \frac{(x-x_{0})(x-x_{1})}{(x_{1}-x_{0})(x_{1}-x_{2})} ,\ dx\]
\[\alpha_{2} = \int_{a}^{b} \frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})} \, dx\]
下の条件から、式変形を行っていく。
\[s = \frac{x-a}{h}より、\]
\[dx = hds\]
\[x:[a,b] → s:[0,2] \]
式変形：
\[\alpha_{0} = \frac{1}{2h^{2}} \int_{0}^{2} (s-1)(s-2)h^{3} \, ds = \frac{h}{3} \]
\[\alpha_{1} = -\frac{1}{h^{2}} \int_{0}^{2} s(s-2)h^{3} \, ds = \frac{4}{3}h \]
\[\alpha_{2} = \frac{1}{2h^{2}} \int_{0}^{2} s(s-1)h^{3} \, ds = \frac{3}{h} \]
以上のようになる。
これらより、以下のことが求まる。
\[\int_{a}^{b} f(x) \, dx \approx \frac{h}{3} (f_{0} + 4f_{1} + f_{2}) \]
これにより、シンプソン公式を求めることができました。
シンプソン公式での実際の計算
１．区間[a,b]を2n等分します。
２．等分したものを端から2つで1つの組にする。
３．２で等分したものからは、３点の点が得られる。
4．この3点から元の曲線を放物戦で近似する。(シンプソン公式の適用)
これを行うことで、以下のような計算式が求まる。
\[\int_{a}^{b} f(x)dx \approx \sum_{k=0}^{n-1} \frac{h}{3} f_{2k}+4f_{2k+1}+f_{2k+2}) \]
これをアルゴリズムで表現していきます。
シンプソン公式のアルゴリズム
Input a,b,n
h ← (b-a)/2n
S ← f(a)+f(b)
For i=1,2,…,n-1
S ← S+4f(a+(2i-1)h)+2f(a+2ih)
end for
S ← S+4f(a+(2n-1)h)
S ← hS/3
Output S
シンプソン公式を用いたソースコード
最後に
今回は、シンプソン公式について紹介してきました。
シンプソン公式は、ラグランジュの補間公式やニュートン・コーツ公式を理解しておかないといけません。
ですので、初見で見た場合はかなり難しく感じるかもしれません。
ですが、きちんとノートに書いてみたりするとわかりやすいので、わかりづらいかもと感じる方は一度ノートでまとめてみるとわかりやすいかと思います。
今回はこれで以上となります。
最後まで読んでいただきありがとうございました！
