logo

シンプソンの公式 📂数値解析

シンプソンの公式

定義

20190611_151154.png

f:[a,b]Rf : [a,b] \to \mathbb{R}[a,b][a,b]で積分可能で、[a,b][a,b]h:=ban\displaystyle h:= {{b-a} \over {n}}の間隔で等しく分割したa=x0<<xn=ba = x_{0} < \cdots < x_{n} = bみたいなノードポイントで割ったとする。以下で定義された数値積分オペレータIn2I_{n}^{2}シンプソンのルールという。 In2(f):=k=1n/2h3[f(x2k2)+4f(x2k1)+f(x2k)] I_{n}^{2} (f) := \sum_{k=1}^{n/2} {{h} \over {3}} \left[ f(x_{2k-2}) + 4 f( x_{2k-1} ) + f(x_{2k} ) \right]

定理

fC4[a,b]f \in C^4 [a,b]とする。シンプソンのルールのエラーE12E_{1}^{2}と漸近エラーE~n2\tilde{E}_{n}^{2}は以下の通り。

  • [1]: E12(f)=h590f(4)(ξ)E_{1}^{2} (f) = - {{h^5} \over {90}} f^{(4)} ( \xi )
  • [2]: E~n2(f)=h4180[f(3)(b)f(3)(a)]\tilde{E}_{n}^{2} (f) = - {{ h^4 } \over {180}} [ f^{(3)} (b) - f^{(3)} (a) ]

説明

In2(f)I_{n}^{2} (f)を展開すると次のようになる。 In2(f)=h3[f(x0)+4f(x1)+2f(x2)+4f(x3)+2f(x4)++2f(xn2)+4f(xn1)+f(xn)] \begin{align*} I_{n}^{2} (f) =& {{h} \over {3}} [ f(x_{0}) + 4 f ( x_{1} ) + 2 f( x_{2} ) + 4 f ( x_{3} ) + 2 f ( x_{4} ) + \cdots \\ & + 2 f (x_{n-2} ) + 4 f ( x_{n-1} ) + f(x_{n} ) ] \end{align*} 確定積分I(f)=abf(x)dx\displaystyle I (f) = \int_{a}^{b} f(x) dx数値積分を求めるために、線形補間を使った台形則とは異なり、二次補間をした方法だ。

特に注目すべき点は、エラーがE12(f)=h590f(4)(ξ)\displaystyle E_{1}^{2} (f) = - {{h^5} \over {90}} f^{(4)} ( \xi )として計算されることだが、これは積分のために22次元で補間したにも関わらずff33次未満の多項式なら、エラーが正確に00になるという意味だ。

証明 1

[1]

戦略: 二次関数が与えられた関数の二次補間なので、多項式補間の性質を使える。


便宜上c:=(a+b2)\displaystyle c:= \left( {{a+b} \over {2}} \right)と置こう。 I12(f):=(ba6)[f(a)+4f(c)+f(b)] I_{1}^{2} (f) := \left( {{ b - a } \over { 6 }} \right) \left[ f(a) + 4 f \left( c \right) + f(b) \right] これは区間[a,b][a,b]ff二次補間してその関数の積分値としてI(f)I(f)を近似したものだ。

ニュートンの差分公式: pn(x)=i=0nf[x0,,xi]j=0i1(xxj)p_{n} (x) =\sum_{i=0}^{n} f [ x_{0} , \cdots , x_{i} ] \prod_{j=0}^{i-1} (x - x_{j} )

3点a,c,b\displaystyle a , c , bとすべてのx[a,b]x \in [a,b]に対してf(x)=p2+1(x)f(x) = p_{2+1} (x)とする。 E12(f):=I(f)I12(f)=ab[f(x)p2(x)]dx=ab[p2+1(x)p2(x)]dx=ab(xa)(xc)(xb)f[a,c,b,x]dx \begin{align*} \displaystyle E_{1}^{2} (f) :=& I(f) - I_{1}^{2} (f) \\ =& \int_{a}^{b} \left[ f(x) - p_{2} (x) \right] dx \\ =& \int_{a}^{b} \left[ p_{2+1} (x) - p_{2} (x) \right] dx \\ =& \int_{a}^{b} (x-a)(x-c)(x-b) f [a,c,b,x] dx \end{align*} ここでwwを次のように定義する。 w(x):=ab(ta)(tc)(tb)dt w(x) := \int_{a}^{b} (t-a)(t-c)(t-b) dt (ta)(tc)(tb)(t-a)(t-c)(t-b)t=ct = cを中心とし奇関数であるため、w(b)=0w(b)=0であり、上端と下端が等しい場合、積分値が00になるので、w(a)=0w(a) = 0だ。すると部分積分法と差分に関する微分式に従って E12(f)=abw(x)f[a,c,b,c]dx=[w(x)f[a,c,b,x]]ababw(x)ddxf[a,c,b,x]dx=abw(x)f[a,c,b,x,x]dx \begin{align*} \displaystyle E_{1}^{2} (f) =& \int_{a}^{b} w’(x) f [a,c,b,c] dx \\ =& \left[ w(x) f [ a,c,b, x] \right]_{a}^{b} - \int_{a}^{b} w(x) {{d} \over {dx}} f [a,c,b,x] dx \\ =& - \int_{a}^{b} w(x) f [a,c,b,x,x] dx \end{align*} 微積分学の基本定理によりw(x)0w(x) \ge 0であるので、積分の平均値定理を使える。

積分の平均値定理: 閉区間[a,b][a,b]で関数ffが連続でw(x)0w(x) \ge 0が積分可能ならば、abf(x)w(x)dx=f(η)abw(x)dx\displaystyle \int_{a}^{b} f(x) w(x) dx = f( \eta ) \int_{a}^{b} w(x) dxを満たすη\eta[a,b][a,b]に少なくとも一つ存在する。

そして、差分の性質により、 E12(f)=f[a,c,b,η,η]abw(x)dx=f(4)(ξ)24[415h5]=h590f(4)(ξ) \begin{align*} \displaystyle E_{1}^{2} (f) =& - f [a,c,b,\eta,\eta] \int_{a}^{b} w(x) dx \\ =& - {{ f^{(4)} ( \xi ) } \over {24}} \left[ {{4} \over {15}} h^5 \right] \\ =& - {{h^5} \over {90}} f^{(4)} ( \xi ) \end{align*} 何らかのη,ξ[a,b] \eta , \xi \in [a,b]が存在することになる。

証明[2]

戦略: リーマン和を導き出せば、後は微積分学の基本定理によって自然と導かれる。ただし、シンプソンのルールはk=1n/2\displaystyle \sum_{k=1}^{n/2}のように加算を行うので、リーマン和を導く過程で22\displaystyle {{2} \over {2}}を掛けるトリックが使われる。


定理[1]により、実際のI(f)I(f)In2(f)I_{n}^{2} (f)の誤差はいくつかのξk[x2(k1),x2k]\xi_{k} \in [x_{2(k-1)}, x_{2k} ]に対して次のように計算される。 En2(f)=I(f)In2(f)=k=1n/2(h590f(4)(ξk)) \begin{align*} \displaystyle E_{n}^{2} (f) =& I (f) - I_{n}^{2} (f) \\ =& \sum_{k=1}^{n/2} \left( - {{ h^5 } \over { 90 }} f^{(4)} ( \xi_{k} ) \right) \end{align*} これについて、 limnEn2(f)h4=limn1h4k=1n/2(h590f(4)(ξk))=limn1h422k=1n/2(h590f(4)(ξk))=1180limnk=1n/22hf(4)(ξk)=1180limnk=1n/2ban/2f(4)(ξk)=1180abf(4)(x)dx=1180[f(3)(b)f(3)(a)] \begin{align*} \lim_{n \to \infty} {{ E_{n}^{2} (f) } \over { h^4 }} =& \lim_{n \to \infty} {{1} \over {h^4}} \sum_{k=1}^{n/2} \left( - {{ h^5 } \over { 90 }} f^{(4)} ( \xi_{k} ) \right) \\ =& \lim_{n \to \infty} {{1} \over {h^4}} {{2} \over {2}} \sum_{k=1}^{n/2} \left( - {{ h^5 } \over { 90 }} f^{(4)} ( \xi_{k} ) \right) \\ =& - {{ 1 } \over { 180 }} \lim_{ n \to \infty} \sum_{k=1}^{n/2} 2 h f^{(4)} ( \xi_{k} ) \\ =& - {{ 1 } \over { 180 }} \lim_{ n \to \infty} \sum_{k=1}^{n/2} {{b-a} \over {n/2}} f^{(4)} ( \xi_{k} ) \\ =& - {{ 1 } \over { 180 }} \int_{a}^{b} f^{(4)} (x) dx \\ =& - {{ 1 } \over { 180 }} [ f^{(3)} (b) - f^{(3)} (a) ] \end{align*} したがって、 limnE~n(f)En(f)=1 \lim_{n \to \infty} {{\tilde{E}_{n} (f) } \over { E_{n} (f) }} = 1

En2(f)E~n2(f)=h4180[f(3)(b)f(3)(a)] E_{n}^{2} (f) \approx \tilde{E}_{n}^{2} (f) = - {{ h^4 } \over { 180 }} [ f^{(3)}(b) - f^{(3)}(a) ]


  1. Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p257~258. ↩︎