3 シンプソンの公式

3.1 基本的な考え方

台形公式の考え方は簡単であるが,精度はあまりよくない.そこで,よく似た考え方で精 度が良いシンプソンの公式を説明する.台形公式は,分割点の値を一次関数(直線)で近似 を行い積分を行った.要するに折れ線近似である.ここで,1次関数ではなく,高次の関 数で近似を行えばより精度が上がることは,直感的に分かる.

2次関数で近似を行うことを考える.2次関数で近似するためには,3点必要である.3つの 分点をそれぞれ, $ (x_j,\,x_{j+1},\,x_{j+2})$とする.そして,この2次関数を$ P(x)$と する.$ P(x)$はラグランジュ補間に他ならないので,


となる.図3に示すとおりである.
図 3: 元の関数を区間 $ [x_j,x_{j+2}]$を2次関数で近似する
\includegraphics[keepaspectratio, scale=1.0]{figure/simpson.eps}

これを,区間 $ [x_j,\,x_{j+2}]$で積分する.紙面の都合上,式 (11)の右辺を各項毎に積分を行う.まず,右辺第1項で あるが,それは以下のようになる.

\begin{equation*}
% latex2html id marker 1155
\begin{aligned}\text{式(\ref{eq:si...
...i^2}{2}+2h^2\xi\right]_0^{2h}\\ &=\frac{h}{3}f(x_j) \end{aligned}\end{equation*}

同様に,第2,3項を計算すると

式(11)右辺第2項の積分 $\displaystyle =\frac{4h}{3}f(x_{j+1})$ (12)
式(11)右辺第3項の積分 $\displaystyle =\frac{h}{3}f(x_{j+2})$ (13)

となる.以上より,近似した2次関数$ P(x)$の範囲 $ [x_j,\,x_{j+2}]$の積分は,

$\displaystyle \int_{x_j}^{x_{j+2}}P(x)dx =\frac{h}{3}\left\{f(x_j)+4f(x_{j+1})+f(x_{j+2})\right\}$ (14)

となる.

これは,ある区間 $ [x_j,\,x_{j+2}]$の積分で,その巾は$ 2h$である.区間$ [a,\,b]$にわ たっての積分$ S$は,式(15)を足し合わせればよい.ただし, $ j=0,2,4,6$と足し合わせる.

\begin{equation*}\begin{aligned}S&=\frac{h}{3}\left[f(x_0)+4f(x_1)+f(x_2)\right]...
...}\left.\cdots+2f(x_{N-2})+4f(x_{N-1})+f(x_N)\right] \end{aligned}\end{equation*}

これが,シンプソンの公式と呼ばれるもので,先ほどの台形公式よりも精度が良い.誤差 は,$ N^{-4}$に比例する.分割数$ N$を10倍にすると,誤差は$ 1/10000$程度にになること が期待できる.

注意シンプソンの公式を使う場合,分割数$ N$は偶数でなくてはならない.これに 注意して,プログラムを作成しよう.

3.2 シンプソンの公式の実際の数値積分

実際にシンプソンの公式を使って数値積分を行う場合,台形公式を上手に使う.台形公式からど うやってシンプソンの公式を計算するのか?--という驚きの疑問が湧くだろう.それをこ れから示す.ここでは,基礎となる計算のアルゴリズムとフローチャートを示す.ここの お話は,主に文献 [1]を参考にした2

3.2.1 計算アルゴリズム

区間$ [a,b]$$ N$等分して関数$ f(x)$を数値積分する.この場合の台形公式とシンプソンの公式 を比較する.比較しやすくするために,台形公式(4)とシンプソンの公式 (16)を少し書き直す.

$\displaystyle T$ $\displaystyle =h\left[ \frac{1}{2}f(a)+f(a+h)+f(a+2h)+f(a+3h)+f(a+4h)+\cdots+\frac{1}{2}f(a)f(a+Nh) \right]$ (16)
$\displaystyle S$ $\displaystyle =h\left[ \frac{1}{3}f(a)+\frac{4}{3}f(a+h)+\frac{2}{3}f(a+2h)+\frac{4}{3}f(a+3h)+\frac{2}{3}f(a+4h)+\cdots+\frac{1}{3}f(a+Nh) \right]$ (17)

これから,面白いことに気がつくだろう.どちらも,関数の値を計算するポイントは$ N+1$個 である.しかし,計算された値の加算の方法が異なる.そのようすを図 4にしめす.加算の方法が異なるだけで,台形公式の精度は 2次,シンプソンの公式は4次になる.

4や式(17)や(18)を見 ると台形公式とシンプソンの公式はよく似ており,なんらかの関係があることは容易に推 測がつく.それらの関係を考えることにする.区間$ [a,b]$$ N=2^n$で分割して,数値積 分を行うとする.この時の台形公式の結果を$ T_n$,シンプソンの公式の結果を$ S_n$とす る.すると,それらには,

$\displaystyle S_{n+1}=\frac{4}{3}T_{n+1}-\frac{1}{3}T_n$ (18)

の関係がある.これは,次のように右辺を計算すれば簡単に証明できる.これを計算する ために,$ S_{n+1}$$ T_{n+1}$のときのステップ幅を $ h=(b-a)/2^{n+1}$とする.こうする と,$ T_{n}$のときのステップ幅は$ 2h$となる.これを利用すると,右辺は

$\displaystyle \frac{4}{3}T_{n+1}-\frac{1}{3}T_n$ $\displaystyle =\frac{4h}{3}\left[ \frac{1}{2}f(a)+f(a+h)+f(a+2h)+f(a+3h)+\cdots+f(b-h)+\frac{1}{2}f(a)f(b) \right]$    
  $\displaystyle \qquad-\frac{2h}{3}\left[ \frac{1}{2}f(a)+f(a+2h)+f(a+4h)+f(a+6h)+\cdots+f(b-2h)+\frac{1}{2}f(a)f(b) \right]$    
  $\displaystyle =\frac{h}{3}\left[ f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+\cdots+2f(b-2h)+4f(b-h)+f(b) \right]$ (19)

となる.

これまでの結果から,$ N=2^{n}$$ N=2^{n+1}$に分割したときの台形公式の積分の値が分 かると,$ N=2^{n+1}$のときのシンプソンの公式での積分値が計算できると言える.式 (19)のとおりである.しかし,式(18)を 直に計算する方法とどれほどのメリットがあるのだろうか?--という疑問も湧くだろう. これだけだと,少し計算が早いかもしれないが,それほど大きなメリットはない.ただ, 台形公式と比べると,大変大きなメリットがある.台形公式とほとんど同じ計算で,シン プソンの公式の精度が得られるのである.ロンバーグ積分まで考えると,さらに大きなメ リットが生じる.

図 4: 台形公式とシンプソンの公式の比較.関数の曲線の上が台形公式とシンプソ ンの公式で積分を計算するとき重み.下:台形公式 上:シンプソンの公式.
\includegraphics[keepaspectratio,scale=1.0]{figure/trapezoid_simpson.eps}

3.2.2 計算アルゴリズム

シンプソンの公式で数値積分を行う場合,式(18)を計算しても良いが, 台形公式を使う漸化式(19)を使うことを勧める.台形公式 からシンプソンの公式の値を求めるには,次のような手順で計算する.ただし,計算精度 は $ \varepsilon$とする.
  1. 区間の区切りを$ N=2$として台形公式で積分$ T_1$を計算する.計算結果は, $ T^{old}$とする.
  2. 区間の区切りを$ N=2$としてシンプソンの公式で積分$ S_1$を計算する.計算結果 は$ S^{old}$とする.
  3. 目的の精度に達するまで,以下を繰り返す.
    1. $ N$の値を2倍にする.それにともない,ステップ幅$ h$を半分にする.

        $\displaystyle N\leftarrow 2N$   $\displaystyle h\leftarrow h/2$   (20)

    2. 新しい台形公式の積分値を計算する.

        $\displaystyle t=f(a+h)+f(a+3h)+f(a+5h)+f(a+7h)+\cdots+f(a+(N-1)h)$ (21)
        $\displaystyle T^{(new)}=\frac{T^{(old)}}{2}+ht$ (22)

    3. 新しいシンプソンの公式の積分値を計算する.

      $\displaystyle S^{(new)}=\frac{4T^{(new)}-T^{(old)}}{3}$ (23)

    4. 以下の条件を満たすならば,目的の計算精度に達したと判断して,繰り返 し文から脱出する.

      $\displaystyle \frac{\vert S^{(new)}-S^{(old)}\vert}{S^{(new)}}<\varepsilon$ (24)

    5. 目的の精度に達していなければ,

        $\displaystyle S^{(old)}\leftarrow S^{(new)}$   $\displaystyle T^{(old)}\leftarrow T^{(new)}$   (25)

      として,最初から繰り返す.
  4. シンプソンの公式による数値積分の結果$ S^{(new)}$を表示する.

ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成19年1月19日


no counter