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次関数で近似する

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

$$% 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}$$

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

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

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

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

$$\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{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}$$

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

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

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

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

3.2.1 計算アルゴリズム

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

$$T =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)

$$S =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$とす る．すると，それらには，

$$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$となる．これを利用すると，右辺は

$$\frac{4}{3}T_{n+1}-\frac{1}{3}T_n =\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]$$

$$\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]$$

$$=\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: 台形公式とシンプソンの公式の比較．関数の曲線の上が台形公式とシンプソ ンの公式で積分を計算するとき重み．下:台形公式　上:シンプソンの公式．

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$を半分にする．

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

      
      

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

      $$t=f(a+h)+f(a+3h)+f(a+5h)+f(a+7h)+\cdots+f(a+(N-1)h)$$ (21)

      $$T^{(new)}=\frac{T^{(old)}}{2}+ht$$ (22)

      
      

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

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

      
      

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

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

      
      

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

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

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

ホームページ: Yamamoto's laboratory
著者: 山本昌志