記号の世界ゟ

このブログでは, 数学書などの書評を書きます。また、受験などの勉強法をまとめます。

シンプソンの公式と誤差評価

今回は数値積分の一つの手法であるシンプソンの公式を紹介します. この公式が単なる近似公式ではないことを見ていきます.

(数学が苦手でない人はシンプソンの公式から読むと良いかもしれません)

数値積分とは


関数を定積分することを考えましょう。原始関数を求めれば計算ができるわけですが、実は、普段使う関数の原始関数は普段使う関数では書けないことが多々あります。例えば、 \displaystyle \frac{\sin x}{x} や \displaystyle e^{-x^2/2} などは、重要であるにも関わらず、原始関数が"簡単に"書けないことが知られています。

今回は定積分, つまり, 実数を求めたいわけなので, その近似値を求めるという方向で考えていきます. 特に, 数値計算で積分値を計算するということを考えます. 定積分の数値計算において, ポイントは二つあると考えています.
・有限個のデータを使う
・任意の(有理数での)関数値が分かる
コンピュータを使うことから前者は自明でしょう. 一方で, 関数 f(x) の積分をするときに, 好きな有理数 x での値が分かる, つまり, "使う"データは有限個でも"使える"データは無限個というのは重要です. というのも, 使えるデータが有限個で決まっているとなると, 使える手法に制限がかかります. また, 実際に決まった有限個のデータから推定値を出すという方向で考えることもあるので, それとは区別しましょうという意味です. よって, 問題は以下のようになります.


与えられるもの: 関数 f(x), x \in \mathbb{Q} と積分区間 [a, b] , a, b \in \mathbb{Q}


求められるもの: 有限個の点 c_1< \dots < c_n \in \mathbb{Q} での値 f(c_1), \dots, f(c_n) から, \int_a^b f(x) dx の近似値を求めるアルゴリズム.
特に, 使う点の個数 n が大きくなるほど精度が良くなるもの.


さらに, 問題を絞りましょう.
使う点の間隔は一定, つまり, a_2 - a_1 = a_3-a_2 = \dots =a_n - a_{n-1} = h とします.　すると, k + 1個のデータが並ぶ区間 [ a_1, a_k] における定積分アルゴリズムを与えると、それを他の区間にも繰り返し適用することで全体の積分値が得られます。以下では、台形公式とシンプソンの公式と呼ばれるものを紹介しますが, それぞれ, k = 2, k=3 の場合の一つの手法であり, それを全区間に繰り返すことは省略します.

台形公式


まず, 2 点でのデータが分かったときに, その間を積分区間とする定積分を計算することにしましょう. 関数 f(x) を積分するかわりに, 二点 (a, f ( a ) ), (b, f(b) ) を通る直線
\displaystyle \quad \hat{f} (x) = \frac{f(b) - f(a)}{b - a} (x - a) + f(a)
を積分します. すると, 一次関数の積分は簡単にできるので, 近似値が分かるわけです.

実際に計算してみると,
\displaystyle \begin{align*} \quad \int_a^b \hat{f} (x ) dx &= \frac{f(b) - f(a) }{ b - a } \frac{(b-a)^2}{2} + f(a) (b-a) \\ &= \frac{b-a}{2} (f(a) + f(b) ) \end{align*}
となるので以下の近似公式が得られます.

(台形公式)
関数 f(x) の aから bまでの定積分を
\displaystyle \quad \int_a^b f(x) dx \approx \frac{b-a}{2} (f(a) + f(b) )
の右辺で求めることを台形公式という.

図を書いてみるとすぐに分かるように, 求める図形を台形で近似したものになっています.


実際に使ってみましょう. 公式の作り方から一次関数を積分すると厳密な答えが出ることは自明です. そこで x^2 を積分しましょう. x = 1 から x=2 まで積分すると, 近似値は
\displaystyle \quad \frac{2-1}{2} (2^2 + 1^2) = \frac{5}{2} = 2.5
であり, 実際の値は
\displaystyle \quad \int_1^2 x^2 dx = \left[ \frac{x^3}{3} \right]_1^2 = \frac{7}{3} = 2.333\dots
となり, 悪くはないですが, それなりに実際の値からはズレますね.
多項式以外にも適用してみましょう. e^x を x = 1 から x=2 まで積分すると, 近似値は
\displaystyle \quad \frac{2-1}{2} (e^2 + e^1) = \frac{1}{2}(e^2 + e) = 5.053\dots
であり, 実際の値は
\displaystyle \quad \int_1^2 e^x dx = \left[ e^x \right]_1^2 = e^2 - e = 4.670\dots
となります. 式の形はずいぶん違いますが, それなりの近似値が出ている気がします.
上で述べたように, たくさんの点を取り細かく近似することで, いくらでも実際の値に近づくので, 問題があるわけではありません. しかし, もっと良い方法を考えたいものです.

シンプソンの公式

次は, 等間隔に並ぶ 3 点でのデータを使う場合を考えます. 3 点が分かればそれらを通る 2次関数が定まるので, 関数を 2次関数で近似することにしましょう. (a ,\, f(a)) ,\, (b ,\, f(b) ) に加えてそれらの中点 ( (a+b)/2 ,\, f( (a+b)/2) ) の3点を通る 2次関数を \hat{f} (x)と書くことにします. すると, \int_a^b \hat{f} (x) dxを \int_a^b f(x) dx の近似値であると考えることができます.

\int_a^b \hat{f} (x) dx を計算するために, 例えば, ラグランジュ補間なんかを使って \hat{f} (x)を計算することができます. 結構計算は大変です*1. しかし, 積分の結果を知っていれば, その証明は簡単にできるので, 今回はそうします.

補題
g (x) を 2次関数とすると,
\displaystyle \quad \int_a^b g (x) dx = \frac{b-a}{6} \left( g (a) + 4 g \left(\frac{a+b}{2} \right) + g (b) \right)
となる.

（証明） g(x) = \alpha x^2 + \beta x + \gammaとおく. すると,
\displaystyle \begin{align*} \quad \int_a^b g(x) dx &= \frac{\alpha}{3} (b^3 - a^3) + \frac{\beta}{2} (b^2 - a^2) + \gamma (b - a)\\ &= \frac{b-a}{6} \left\{ 2\alpha (b^2 + ba + a^2) + 3\beta (b + a) + 6 \gamma \right\} \end{align*}
であり,
\displaystyle \begin{align*} \quad g(a) + 4 g \left(\frac{a+b}{2} \right) + g (b) =&\, \alpha a^2 + \beta a + \gamma \\ &+\alpha (a+b)^2 + e \beta ( a + b) + 4 \gamma \\ &+ \alpha b^2 + \beta b + \gamma \\ =&\, 2\alpha (b^2 + ba + a^2) + 3\beta (b + a) + 6 \gamma \end{align*}
となることから, 示すべき等式が証明できた. □

今回は, f (x) と \hat{f} (x)が x = a ,\, (a+b)/2 ,\, bでの値が等しいとするので, 補題より
\displaystyle \quad \int_a^b f (x) dx \approx \quad \int_a^b \hat{f} (x) dx = \frac{b-a}{6} \left \{ f(a) + 4 f \left( \frac{a+b}{2} \right) + f(b) \right \}
が分かります. よって, f(x)を 2次関数で近似することによって以下の近似式を得ることができました.

(シンプソンの公式)
関数 f(x) の aから bまでの定積分を
\displaystyle \quad \int_a^b f(x) dx \approx \frac{b-a}{6} \left \{ f(a) + 4 f \left( \frac{a+b}{2} \right) + f(b) \right \}
の右辺で求めることをシンプソンの公式と呼ぶ.

以下では右辺を
\displaystyle \quad S_a^b (f) = \frac{b-a}{6} \left \{ f(a) + 4 f \left( \frac{a+b}{2} \right) + f(b) \right \}
と書くことにしましょう.

3次の多項式に適用→驚きの結果


さて, シンプソンの公式の作り方から, ２次関数までは厳密な結果が出ることが明らかなので, 三次関数 x^3 の積分を計算してみましょう. x = 1 から x=2 まで積分すると, 近似値は
\displaystyle \quad \frac{2-1}{6} \left( 1^3 + 4 \left( \frac{3}{2} \right)^3 + 2^3 \right) = \frac{1}{6} \left( ( 1 + \frac{27}{2} + 8 \right) =\frac{15}{4}
であり, 実際の値は
\displaystyle \quad \int_1^2 x^3 dx = \left[ \frac{x^4}{4} \right]_1^2 = \frac{16 - 1}{4} = \frac{15}{4}
となり, 結果が一致してしまいます. たまたまかな?と思って, 他の３次関数に適用してみても, 必ず一致することが分かります. (ぜひ自分でも計算してみてください.)

シンプソンの公式は, 関数を２次関数で近似したので, ２次関数の積分値が正しくでることは当然です. ３次関数を２次関数で近似すると当然関数は違います. しかし, その積分値は２次関数で近似しても正しくでるということが分かったのです. 今回の目標はこの現象を理解することです.

少し脱線しますが, シンプソンの公式でも e^x の積分を計算しておきましょう. e^x を x = 1 から x=2 まで積分すると, 近似値は
\displaystyle \quad \frac{2-1}{6} (e^1 + 4 e^{\frac{3}{2} } + e^2) = \frac{e + 4 e^{\frac{3}{2} } + e^2 }{6} = 4.6723\dots
であり, 実際の値は
\displaystyle \quad \int_1^2 e^x dx = \left[ e^x \right]_1^2 = e^2 - e = 4.6707\dots
となり, 小数第２位までの結果が一致しています. 非常に精度が良いです.

結果の理由1 (計算+幾何的説明)

３次関数に対しては, 近似値ではなく, 本当の公式になっているということは, 実際に代入して計算してみれば簡単に分かります.

また, y = f(x) - \hat{f} (x) のグラフを書いてみると, \displaystyle x = \frac{a + b}{2} で点対称になっていることが分かります. このことから, 本当の関数と近似した関数の差 f(x) - \hat{f} (x) の積分値が 0 になることが分かるので, シンプソンの公式は厳密な値を返すということが分かります.

このように, 理由はいくらでも与えられるのですが, シンプソンの公式の本質に迫るような理由づけがしたいものです. それが次に説明する誤差評価です.

結果の理由2 (シンプソンの公式の誤差評価)

定理（シンプソンの公式の誤差）
関数 f(x) は 4回微分可能で f^{(4)} (x)は連続とする. （つまり, f(x) は C^4 級とする.）
このとき, 区間 [a ,\, b ] における |f^{(4)} (x)|の最大値を M, h = (b- a)/2とすると
\displaystyle \quad \left| \int_a^b f(x) dx - S_a^b (f) \right| = \frac{M}{90} h^5
が成り立つ.

誤差評価の証明にはテイラー展開が必要です. ここで使うテイラー展開は数IIIの部分積分さえ分かっていれば証明ができます.

まず, 積分の基本定理より
\displaystyle \quad \int_a^x f^{\prime} (t) dt = f(x) - f(a)
なので,
\displaystyle \quad f(x) = f(a) + \int_a^x f^{\prime} (t) dt
となります. さらに, 部分積分の公式より
\displaystyle \begin{align*} \int_a^x f^{\prime} (t) dt &= - \int_a^x f^{\prime} (t) (x - t)^{\prime} dt \\ &= - \left [f^{\prime} (t) (x-t) \right ]_a^x + \int_a^x f^{\prime \prime} (t) (x-t) dt \\ &= f^{\prime} (a) (x-a)+ \int_a^x f^{\prime \prime} (t) (x-t) dt \end{align*}
が成立するので,
\displaystyle \quad f(x) = f(a) + f^{\prime} (a) (x-a) + \int_a^xf^{\prime \prime} (t) (x-t) dt
となります. これを繰り返すわけですが, もう一回分だけ書いておくと,
\displaystyle \begin{align*} \int_a^xf^{\prime \prime} (t) (x-t) dt &= - \int_a^x f^{\prime \prime} (t) \left( \frac{(x-t)^2}{2} \right)^{\prime} dt \\ &= \frac{ f^{\prime \prime} (a)}{2} (x-a)^2 + \frac{1}{2} \int_a^x f^{\prime \prime \prime } (t) \, (x-t)^2dt \end{align*}
となるので,
\displaystyle \quad f(x) = f(a) + f^{\prime} (a) \, (x-a) + \frac{ f^{\prime \prime} (a)}{2} \, (x-a)^2 + \frac{1}{2} \int_a^x f^{\prime \prime \prime } (t) \, (x-t)^2 dt
と書けます. これを繰り返して得られる式を a周りのテイラー展開と言います.

補題（テイラー展開）
関数 f(x)が n+1回微分可能で f^{(n+1)} (x)が連続なら
\displaystyle \quad f(x) = \sum_{k=0} ^n \frac{f^{(k)} (a)}{k !} (x - a) ^k + \frac{1}{n!} \int_a^x f^{(n+1)} (t) \, (x-t)^ndt
と書ける.

この補題を用いて, 誤差評価の証明をします.

（証明）
計算を簡単にするために, m = (a+b)/2,\, h = (b - a)/2 (= b - m = m - a) とおく.
また,
\displaystyle \quad e(h) = \int_{m - h} ^{m+h} f(x) dx - S_{m - h}^{m + h} (f) \, ( = \int_a^b f(x) dx - S_a^b (f))
とおいて, この関数をテイラー展開を用いて評価する. そのために, この関数の微分を計算する.
ここで,
\displaystyle \quad S_{m - h}^{m + h} (f) = \frac{h}{3} \left\{ f(m-h) + 4 f(m) + f(m+h) \right\}
であることに注意する. 以下で, 計算が中途半端に見える部分があるが, 機械的に計算できることと計算量が少なくなることを意図している. まず一階の導関数は
\displaystyle e^{\prime} (h ) = f(m+h) + f(m- h) - \frac{1}{3} \left\{ f(m-h) + 4 f(m) + f(m+h) \right\} - \frac{h}{3} \left\{ - f^{\prime} (m-h) + f^{\prime} (m+h) \right\}
であり, 微分係数は \displaystyle\quad e^{\prime} (0) = 0 となる. 2階の導関数は
\displaystyle \begin{align*} e^{\prime \prime} (h) =& \, f^{\prime} (m + h) - f^{\prime} (m - h) - \frac{1}{3} \left\{ - f^{\prime} (m-h) + f^{\prime} (m+h) \right\} \\ &- \frac{1}{3} \left\{ - f^{\prime} (m-h) + f^{\prime} (m+h) \right\} - \frac{h}{3} \left\{ f^{\prime \prime } (m-h) + f^{\prime \prime} (m+h) \right\} \end{align*}
であり, 微分係数は e^{\prime \prime} (0) = 0 となる. 3階の導関数は
\displaystyle \begin{align*} e^{\prime \prime \prime } (h) =& \, f^{\prime \prime } (m + h) + f^{\prime \prime} (m - h) - \frac{1}{3} \left\{ f^{\prime \prime } (m-h) + f^{\prime \prime} (m+h) \right\} \\ &- \frac{1}{3} \left\{ f^{\prime \prime} (m-h) + f^{\prime \prime } (m+h) \right\} - \frac{1}{3} \left\{ f^{\prime \prime } (m-h) + f^{\prime \prime} (m+h) \right\} \\ &- \frac{h}{3} \left\{ - f^{\prime \prime \prime} (m-h) + f^{\prime \prime \prime } (m+h) \right\} \\ = & - \frac{h}{3} \left\{ - f^{\prime \prime \prime} (m-h) + f^{\prime \prime \prime } (m+h) \right\} \end{align*}
となり, やはり微分係数は e^{\prime \prime \prime} (0) = 0 となる. のちの不等式評価で便利なように fの 4 次の導関数のみが現れるように, e (h)の 4 次の導関数を計算すると
\displaystyle \begin{align*} e^{(4)} (h) & = - \frac{1}{3} \left\{ - f^{\prime \prime \prime} (m-h) + f^{\prime \prime \prime } (m+h) \right\} - \frac{h}{3} \left\{ f^{(4)} (m-h) + f^{(4) } (m+h) \right\} \\ & = - \frac{1}{3} \int_{m - h}^{m + h} f^{(4)} (t) dt - \frac{h}{3} \left\{ f^{(4)} (m-h) + f^{(4) } (m+h) \right\} \end{align*}
と書ける.

以上の結果, テイラー展開で n = 3とした場合にあてはめると,
\displaystyle \quad e(h) = \int_0^h \frac{(h - t)^3}{3!} e^{(4)} (t) dt
となる. 絶対値をとり, 不等式を評価していくと
\displaystyle \quad \left| e(h) \right| \leq \int_0^h \frac{(h - t)^3}{3!} \left|e^{(4)} (t) \right| dt
であり, M = max \left| f^{(4)} (t) \right|とおくと,
\displaystyle \begin{align*} \left| e^{(4)} (h) \right| & \leq \frac{1}{3} \int_{m - h}^{m + h} \left| f^{(4)} (t) \right|dt - \frac{h}{3} \left| f^{(4)} (m-h) + f^{(4) } (m+h) \right|\\ &\leq \frac{1}{3} \int_{m-h}^{m+h} M dt + \frac{h}{3}2M\\ &= \frac{4Mh}{3} \end{align*}
となるので,
\displaystyle \quad | e(h) | = \frac{4M}{3! \times 3} \int_0^h (h - t)^3 t dt
となる. 最後に積分を計算すると,
\ \displaystyle \begin{align*} \int_0^h (h - t)^3 t dt &= \int_0^h \left\{ - \frac{(h - t)^4}{4} \right\}^{\prime} t dt \\ &= \left[ - \frac{(h - t)^4}{4} t \right ]_0^h + \int_0^h \frac{(h - t)^4}{4} dt \\ &= \left[- \frac{(h - t)^5}{5\cdot 4} \right ]_0^h \\ &=\frac{h^5}{4} \end{align*}
となるので, 求めたかった不等式
\displaystyle \| e(h) \| \leq \frac{4M}{3!\cdot 3 \cdot 5\cdot 4} h^5 = \frac{M}{90} h^5
を得ることができた.□


さて, 3 次関数は 4 階微分が 0 です. よって, シンプソンの公式を 3 次関数に適用すると, 誤差は 0 となり, 厳密に面積を求める公式となるのです.

まとめ

数値計算では, いろんな考え方で公式を導く方法があります. シンプソンの公式も他の見方をすることができます. どの考え方を採用するかで, 他の問題への一般化や新しい問題に対する適用範囲が変わってくるので, 実は公式だけでなくその背後の考え方=アルゴリズムの構造も大切だったりします. 今回の記事では, 数値計算という分野の面白さの一端を見せることができていれたなら, 幸いです.
今回の証明では,
吉田耕作『私の微積分法』
を参考にしました. 高校生でも大学の解析学を学べる良書だと思います. 気になった人はぜひ手にとってみてください.

↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓

↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑

*1:高校生に教えるときには, 計算練習もかねて実際に計算してもらいます