シンプソン公式

Simpson公式（Simpson's rule，シンプソン法）

直感的な考え方

区分求積法，台形公式からの類推で「より次数の高い関数で近似できた方が誤差が減少する可能性があるのではないか？」と考えられる．つまり，一定（区分求積法），そして直線（台形公式）の次は，２次関数である．以下の図に示すように２次関数を表現する場合には３点を利用する必要があり，２つの区分を利用する必要がありそうである．

ここで，３点を $x_{0}, x_{1}, x_{2}$ ，その関数値を $f_{0}, f_{1}, f_{2}$ として考えよう．まず２次関数が $x_{0}$ で $f_{0}$ を通ることから，２次関数の係数を $a, b, c (= f_{0})$ とするとまずは以下のような式を立てることができる．

\begin{matrix} (1) & f (x) = a (x - x_{0})^{2} + b (x - x_{0}) + f_{0} \end{matrix}

また，同時に他の２点を通るということから，以下の二つの式が得られる．

\begin{matrix} (2) & f_{1} = a (x_{1} - x_{0})^{2} + b (x_{1} - x_{0}) + f_{0} \end{matrix}

\begin{matrix} (3) & f_{2} = a (x_{2} - x_{0})^{2} + b (x_{2} - x_{0}) + f_{0} \end{matrix}

ここで， $x_{1} - x_{0} = h$ , $x_{2} - x_{0} = 2 h$ を考慮して整理すると以下のような $a$ と $b$ に関する連立方程式となり，通常の方法で解くと． $a, b$ は以下のように求められる．

\begin{matrix} (4) & a = \frac{f_{0} - 2 f_{1} + f_{2}}{2 h^{2}}, b = \frac{- 3 f_{0} + 4 f_{1} - f_{2}}{2 h} \end{matrix}

さて，ここで式( $1$ )で示した関数を積分してみる．

\begin{array}{rcl} I & = & \int_{x_{0}}^{x_{2}} (a (x - x_{0})^{2} + b (x - x_{0}) + f_{0}) d x \\ = & {[\frac{1}{3} a {(x - x_{0})}^{3} + \frac{1}{2} b {(x - x_{0})}^{2} + f_{0} x]}_{x_{0}}^{x_{2}} \\ = & \frac{1}{3} a {(x_{2} - x_{0})}^{3} + \frac{1}{2} b {(x_{2} - x_{0})}^{2} + f_{0} (x_{2} - x_{0}) \\ = & \frac{1}{3} \frac{f_{0} - 2 f_{1} + f_{2}}{2 h^{2}} 8 h^{3} + \frac{1}{2} \frac{- 3 f_{0} + 4 f_{1} - f_{2}}{2 h} 4 h^{2} + 2 f_{0} h \\ (5) & = & \frac{h}{3} (f_{0} + 4 f_{1} + f_{2}) \end{array}

Newton-Cotes法からの誘導

ニュートン・コーツ(Newton-Cotes)法の定義は以下のとおりである．

\begin{matrix} (6) & I = \int_{a}^{b} f (x) d x \sim \int_{a}^{b} \sum_{i = 0}^{n} f_{i} \cdot l_{i} d x = \sum_{i = 0}^{n} f_{i} \int_{a}^{b} l_{i} d x = \sum_{i = 0}^{n} α_{i} f_{i} \end{matrix}

ここで $α_{i} = \int_{a}^{b} l_{i} d x$ である．また一般的に $l_{i}$ は以下のように定義されている．

\begin{matrix} (7) & l_{i} (x) = \prod_{j = 0 \sim n, j \neq i} \frac{x - x_{j}}{x_{i} - x_{j}} \end{matrix}

ここで３点，つまりシンプソン則は $n = 2$ の場合に対応するので，３つの係数は以下のようになる．

\begin{matrix} (8) & l_{0} = \frac{(x - x_{1}) (x - x_{2})}{(x_{0} - x_{1}) (x_{0} - x_{2})}, l_{1} = \frac{(x - x_{0}) (x - x_{2})}{(x_{1} - x_{1}) (x_{1} - x_{2})}, l_{2} = \frac{(x - x_{0}) (x - x_{1})}{(x_{2} - x_{0}) (x_{2} - x_{1})} \end{matrix}

例えば $α_{0}$ は以下のような積分になる．

\begin{matrix} (9) & α_{0} = \int_{x_{0}}^{x_{2}} \frac{(x - x_{1}) (x - x_{2})}{(x_{0} - x_{1}) (x_{0} - x_{2})} d x = \frac{h}{3} \end{matrix}

\begin{matrix} (10) & α_{1} = \frac{4 h}{3}, α_{2} = \frac{h}{3} \end{matrix}

ここで $x_{1} - x_{0} = h, x_{2} - x_{0} = 2 h, x_{2} - x_{1} = h$ であることを利用している（上記の $α_{0}, α_{1}, α_{2}$ を実際に積分して求めることを課題にする）．

最終的に，区間内の積分は，２次関数の近似で求めた積分と同じになる．

\begin{matrix} (11) & I = \frac{h}{3} (f_{0} + 4 f_{1} + f_{2}) \end{matrix}

精度の評価

区分求積法や台形公式と同様に，シンプソン公式にもテイラー展開による評価を行おう．

まず例によって，元の関数をテイラー展開後に積分した積分値

I

を考える．

\begin{array}{rcl} I & = & \int_{0}^{2 h} (f (x_{i}) + \frac{t}{1!} f^{'} (x_{i}) + \frac{t^{2}}{2!} f^{″} (x_{i}) + \frac{t^{3}}{3!} f^{‴} (x_{i}) + \dots) d t \\ (12) & = & 2 f (x_{i}) h + \frac{4 h^{2}}{2} f^{'} (x) + \frac{8 h^{3}}{6} f^{″} (x_{i}) + \frac{16 h^{4}}{24} f^{‴} (x) + \frac{32 h^{5}}{120} f^{⁗} (x) + \dots \end{array}

次にシンプソン公式の括弧内の係数を除いた各項を計算する（テイラー展開する）．

\begin{array}{rcl} (13) & f (x_{i}) & = & f (x_{i}) \\ (14) & f (x_{i} + h) & = & f (x_{i}) + h f^{'} (x) + \frac{h^{2}}{2} f^{″} (x) + \frac{h^{3}}{6} f^{‴} (x) + \frac{h^{4}}{24} f^{⁗} (x) + \dots \\ (15) & f (x_{i} + 2 h) & = & f (x_{i}) + 2 h f^{'} (x) + \frac{4 h^{2}}{2} f^{″} (x) + \frac{8 h^{3}}{6} f^{‴} (x) + \frac{16 h^{4}}{24} f^{⁗} (x) + \dots \end{array}

これらを係数を合わせて総合すると以下のようになる

\begin{matrix} (16) & I^{'} = 2 h f (x_{i}) + \frac{4 h^{2}}{2} f^{'} (x) + \frac{8 h^{3}}{6} f^{″} (x_{i}) + \frac{16 h^{4}}{24} f^{‴} (x) + \frac{20 h^{5}}{72} f^{⁗} (x) + \dots \end{matrix}

これらの引き算して誤差を評価すると以下のようになる．

\begin{matrix} (17) & I - I^{'} = - \frac{h^{5}}{90} f^{⁗} (x) + \dots \end{matrix}

区分の誤差は $h^{5}$ に比例することがわかった．これまでと同様な議論をすれば，全区間の誤差は $h^{4}$ に比例する．つまり，Simpson公式は４次精度のスキームであることがわかる．まずは，これまでと同様に $y = 3 x^{2}$ の積分で誤差を評価してみよう．

プログラムを確認すると，区分求積法や台形公式との違いは，14-17行の４行である(別途f2などの定義は加えているが）．このプログラムも積分区分が両端が重複するので，より効率的にプログラムを作成できる．興味のある学生さんは効率化を試みてみると良い．

    program simpson_rule
  implicit none
  integer i, j, n
  double precision a, b, h, s, s_exact, x, f0, f1, f2
  a = 0.0
  b = 5.0
  s_exact = 125.0
  do i=1,6
    x = 0.0
    s = 0.0
    n = 2**i
    h = (b-a)/n
    do j=0,n/2-1
      f0 = 3.0*x**2
      f1 = 3.0*(x+h)**2
      f2 = 3.0*(x+2*h)**2
      s = s+h/3.d0*(f0+4*f1+f2)
      x = x+2*h
    end do
    print *, n, abs((s-s_exact)/s_exact)
  end do
end program

実行結果を以下に示す．驚くことに全ての分割数で誤差がゼロである．これはシンプソン公式自身が２次関数を構築するようにできており，積分対象が２次関数の場合は数値解と解析解は一致する(さらに驚くことに３次関数でも誤差はゼロとなる．解析解と数値解の誤差か中点 $x_{i + 1}$ で点対象になっている)．

  2  0.0000000000000000     
 4  0.0000000000000000     
 8  0.0000000000000000     
16  0.0000000000000000     
32  0.0000000000000000     
64  0.0000000000000000

これでは面白くないので，4次関数である

y = 5 x^{4}

と

x

軸に囲まれる範囲

x = 0 \sim 5

の面積で試してみる（厳密解は

S = 3125.0

である）．

  2   4.1666666666666713E-002
 4  2.6041666666667151E-003
 8  1.6276041666671518E-004
16  1.0172526041715172E-005
32  6.3578287779819223E-007
64  3.9736429898766802E-008

これまでと同様にGnuplotによって誤差と分割数の関係を見てみよう．新たに $y = 1 / x^{4}$ という４次精度のラインも加えた結果を示す．計算結果（紫線）は４次精度のスキームであることがわかる．これまで学んだ他のスキームよりも圧倒的に高精度であることがわかる．

	program simpson_rule
	implicit none
	integer i, j, n
	double precision a, b, h, s, s_exact, x, f0, f1, f2
	a = 0.0
	b = 5.0
	s_exact = 125.0
	do i=1,6
	x = 0.0
	s = 0.0
	n = 2**i
	h = (b-a)/n
	do j=0,n/2-1
	f0 = 3.0x*2
	f1 = 3.0(x+h)*2
	f2 = 3.0(x+2h)**2
	s = s+h/3.d0(f0+4f1+f2)
	x = x+2*h
	end do
	print *, n, abs((s-s_exact)/s_exact)
	end do
	end program

	2 0.0000000000000000
	4 0.0000000000000000
	8 0.0000000000000000
	16 0.0000000000000000
	32 0.0000000000000000
	64 0.0000000000000000

	2 4.1666666666666713E-002
	4 2.6041666666667151E-003
	8 1.6276041666671518E-004
	16 1.0172526041715172E-005
	32 6.3578287779819223E-007
	64 3.9736429898766802E-008