積分

山本昌志^{∗}

2005 年 1 月 25 日

1 はじめに

数学の授業で学習したように、どんな複雑な関数でも微分は可能である。一方、積分となるとそうはいかない。積分の学習では、どのようにして積分を行うかといういろいろなテクニックを学んだはずである。微分に比べれば、圧倒的に計算が難しいことも経験済みであろう。

例えば、ガウス分布を表す以下の関数を考える。

f(x) = e^{−x^{2}}  (1)

この関数の微分すること (導関数) は、簡単で

f′(x) = −2xe^{−x^{2}}  (2)

となることは説明の必要がないであろう。今まで学習してきたように、初等関数で表すことができる関数の微分は、初等関数で表現できるのである。要するに微分 (導関数) の値を求めたいときには、人間が実際に微分をして、初等関数を計算すればよいのである。

一方、積分となるとそうはいかない。先の式 (1) の不定積分を初等関数で表すことができない。初等関数からできた関数であろうとも、不定積分は初等関数の範囲を超えることがある。だからといって、定積分の値 (数値) が不定というわけではない。

いろいろ計算をしていると、不定積分はできないが、定積分の値が必要な場合がしばしば訪れる。そのときに、ここで学習する定積分を数値計算で求めるテクニックが使われる。

定積分は、図 1 に示すように面積の計算になる。したがって、直感的にわかりやすく、アルキメデスの時代からあった。一方、微分法はニュートンの時代とすると、およそ 2000 年の開きがあるのである。

∗独立行政法人 秋田工業高等専門学校 電気工学科

図 1: 定積分と面積

2 台形公式

2.1 台形公式の求め方

定積分、

S = ∫_{a}^{b} f(x)dx  (3)

の近似値を数値計算で求めることを考える。積分の計算は、先に示したように面積の計算であるから、図 2 のように台形の面積の和で近似ができるであろう。積分の範囲 [a, b] を N 等分した台形で近似した面積 T は、

\begin{aligned}T & = h\frac{f(a) + f(a + h)}{2} + h\frac{f(a + h) + f(a + 2h)}{2} + h\frac{f(a + 2h) + f(a + 3h)}{2} + ⋯ \\ & + h\frac{f(a + (N − 1)h) + f(a + Nh)}{2} \\ & = \frac{h}{2} ∑_{j=0}^{N−1} [f(a + jh) + f(a + (j + 1)h)] \end{aligned} (4)

となる。これが数値積分の台形公式である。なんのことはない、積分を台形の面積に置き換えているだけである。

2.2 台形公式の誤差について

台形公式による数値積分では、分割数 N を大きくするとその誤差は小さくなることは直感で分かる。それでは、分割数を増やしていくとどのように精度が良くなるのか考えてみよう。

図 2: 積分と台形の面積の比較

まずは、式 4 のある一つの台形の面積と実際の積分の値を比較する。台形の面積 t は、台形公式より、

t = \frac{h}{2}[f(x_{1}) + f(x_{1} + h)]  (5)

となる。これを実際の積分

s = ∫_{x_{1}}^{x_{1}+h} f(x)dx  (6)

と比較することにする。これら 2 つの式の形がぜんぜん違うので比較できないと考えるかもしれないが、このような場合の常套手段がある。このようなときには、テーラー展開をすれば良いのである。式 (5) を x_{1} の周りで、テイラー展開すると

\begin{aligned}t & = \frac{h}{2}[f(x_{1}) + f(x_{1} + h)] \\ & = \frac{h}{2} [f(x_{1}) + f(x_{1}) + f′(x_{1})h + \frac{f′′(x_{1})}{2!}h^{2} + \frac{f′′′(x_{1})}{3!}h^{3} + ⋯] \\ & = f(x_{1})h + \frac{f′(x_{1})}{2}h^{2} + \frac{f′′(x_{1})}{2 × 2!}h^{3} + \frac{f′′′(x_{1})}{2 × 3!}h^{4} + ⋯ \end{aligned}  (7)

となる。これが台形の面積のテイラー展開である。一方、積分の式 (6) もテイラー展開する。これは、

\begin{aligned}s & = ∫_{x_{1}}^{x_{1}+h} f(x)dx \\ & = ∫_{0}^{h} f(x_{1} + ξ)dξ \\ & = ∫_{0}^{h} [f(x_{1}) + f′(x_{1})ξ + \frac{f′′(x_{1})}{2!}ξ^{2} + \frac{f′′′(x_{1})}{3!}ξ^{3} + ⋯] dξ \\ & = [f(x_{1})ξ + \frac{f′(x_{1})}{2}ξ^{2} + \frac{f′′(x_{1})}{3 × 2!}ξ^{3} + \frac{f′′′(x_{1})}{4 × 3!}ξ^{4} + ⋯]_{0}^{h} \\ & = f(x_{1})h + \frac{f′(x_{1})}{2}h^{2} + \frac{f′′(x_{1})}{3 × 2!}h^{3} + \frac{f′′′(x_{1})}{4 × 3!}h^{4} + ⋯ \end{aligned}  (8)

となる。この 2 つの式 (7) と ( 8) が台形での近似とまっとうに積分を行ったときのテイラー展開を表す。これらの式を比べると、刻み巾 h の 2 次まで一致している。異なるのは 3 次以降で、積分の誤差は、

\begin{aligned}|s − t| & = \frac{f′′(x_{1})}{2 × 2!}h^{3} − \frac{f′′(x_{1})}{3 × 2!}h^{3} + O(h^{4}) \\ & = \frac{f′′(x_{1})}{12}h^{3} + O(h^{4}) \end{aligned}  (9)

と表せる。即ち、積分を台形で近似したひとつの区間の誤差は、刻み幅の h^{3} で効いてくるのである。従って、積分のトータルの誤差は、それを区間の個数 N を乗じた

\begin{aligned}|S − T | & = N|s − t| \\ & = N \frac{f′′(x_{1})}{12}h^{3} \\ & = N \frac{f′′(x_{1})}{12} (\frac{b − a}{N})^{3} \\ & = \frac{f′′(x_{1})}{12} \frac{(b − a)^{3}}{N^{2}} \end{aligned}  (10)

となる。要するに積分の誤差は、分割数 N^{2} に反比例する。分割数を 10 倍にすれば、積分の誤差は 1/100 になるわけである。

3 シンプソンの公式

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

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

P (x) = \frac{(x − x_{j+1})(x − x_{j+2})}{(x_{j} − x_{j+1})(x_{j} − x_{j+2})}f(x_{j}) + \frac{(x − x_{j})(x − x_{j+2})}{(x_{j+1} − x_{j})(x_{j+1} − x_{j+2})}f(x_{j+1})

+ \frac{(x − x_{j})(x − x_{j+1})}{(x_{j+2} − x_{j})(x_{j+2} − x_{j+1})}f(x_{j+2}) (11)

となる。図 3 に示すとおりである。

図 3: 元の関数を区間 [x_{j}, x_{j+2}] を 2 次関数で近似する

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

\begin{matrix}式 (11) 右辺第 1 項の積分 & = ∫_{x_{j}}^{x_{j+2}} \frac{(x − x_{j+1})(x − x_{j+2})}{(x_{j} − x_{j+1})(x_{j} − x_{j+2})}f(x_{j})dx \\ = ∫_{0}^{2h} \frac{(x_{j} + ξ − x_{j+1})(x_{j} + ξ − x_{j+2})}{(x_{j} − x_{j+1})(x_{j} − x_{j+2})}f(x_{j})dξ \\ = ∫_{0}^{2h} \frac{(ξ − h)(ξ − 2h)}{(−h)(−2h)}f(x_{j})dξ \\ = \frac{f(x_{j})}{2h^{2}} ∫_{0}^{2h} ξ^{2} − 3hξ + 2h^{2}dξ \\ = \frac{f(x_{j})}{2h^{2}} [\frac{ξ^{3}}{3} − 3h\frac{ξ^{2}}{2} + 2h^{2}ξ]_{0}^{2h} \\ = \frac{h}{3}f(x_{j}) \end{matrix}  (12)

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

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

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

∫_{x_{j}}^{x_{j+2}} P (x)dx = \frac{h}{3} \{f(x_{j}) + 4f(x_{j+1}) + f(x_{j+2})\}  (15)

となる。

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

S = \frac{h}{3} \{f(x_{0}) + 4f(x_{1}) + f(x_{2})\} + \frac{h}{3} \{f(x_{2}) + 4f(x_{3}) + f(x_{4})\} + \frac{h}{3} \{f(x_{4}) + 4f(x_{5}) + f(x_{6})\} + ⋯

⋯ + \frac{h}{3} \{f(x_{N−2}) + 4f(x_{N−1}) + f(x_{N})\}

= \frac{h}{3} \{f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + 4f(x_{3}) + 2f(x_{4}) + 4f(x_{5}) + 2f(x_{6}) + ⋯

⋯ + 2f(x_{N−2}) + 4f(x_{N−1}) + f(x_{N})}

(16)

これが、シンプソンの公式と呼ばれるもので、先ほどの台形公式よりも精度が良い。精度は、 N^{4} に反比例する。

この式から、分割数 N は偶数でなくてはならないことがわかる。これに注意して、プログラムを作成しよう。

4 ロンバーグ積分法

興味のある者は、各自勉強せよ。

5 モンテカルロ積分法

5.1 単純な場合 (２次元の面積)

ここでは、モンテカルロ法 (Monte Carlo method) 呼ばれる乱数を使った一風変わった数値積分の方法を示す。モンテカルロとは有名なモナコのカジノで、博打をするところである。この方法を少し学習すると、名前の由来が分かったような気になるはずである。。

まずは、図 4 の様な円の面積を求めることにする。もちろん、この円の面積は π と分かっているが、これを乱数を用いて数値計算しようと言うのである。円を囲む正方形内部には、ランダムな点がある。正方形内部にある点の個数を N_{s}、円の内部にある個数を N_{c} とする。そして、正方形の面積を S_{s}、円を S_{c} とする。すると、これらには

\frac{S_{c}}{S_{s}} ≃ \frac{N_{c}}{N_{s}}  (17)

の関係がある。ランダムな点の数増加させると、これらは、ほとんど等式で結ばれるであろう。形状が単純な正方形の面積 S_{s} は分かっているので、円の面積 S_{c} は

S_{c} ≃ \frac{N_{c}}{N_{s}} S_{s}  (18)

と計算できる。点の数を増やすと、それがランダムに分布する限り、右辺は円の面積に近づく。正方形の面積と点の数を計算するだけで円の面積が分かるのである。このようにランダムな点 (乱数) を用いて、積分する方法をモンテカルロ積分と言う。

円の面積を求めても大したことはないが、例えば図 5 のような複雑な形状の面積が計算できるとなるとぐっと御利益はある。先ほどの円の面積と同じようにして、計算できるのである。計算したい部分内の点の数 N_{in} は、

N_{in} = (大円Ω_{0}内部の点数) − (大円Ω_{0}内部かつ小円Ω_{1}内部の点数)
− (大円Ω_{0}内部かつ楕円Ω_{2}内部の点数) − (大円Ω_{0}内部かつ正方形Ω_{3}内部の点数) (19)

こうやって、内部の点の数を調べて、後は、式 (18) のようにすれば面積を求めることができる。

このように複雑な領域の積分でも適用でき、応用範囲の広い方法である。積分の精度を上げようとすると、サンプル数 N を増やすことになる。しかし、その精度は、たかだかサンプル数の平方根に比例するのみである。誤差を 1/10 にしたければ、サンプル数を 100 倍にする必要がある。これが、この方法の最大の弱点である。

精度がサンプル数の平方根に比例するのは、完全な乱数を仮定しているからである。それを計算機で発生させることは不可能なので、その生成も問題となる。できるだけ良い乱数を使わなくてはならないが、この問題は奥が深い。ここではこれ以上立ち入らないが、興味のある者は調べてみると良い。

話は変わるが、ひとつ注意を与えておく。実際のプログラムでは、乱数を使って (x, y) を発生させて、それが内部にあるか否かをその場で判断する。配列に座標データを入れて、最後に検査するとメモリーがいくらあっても足りなくなるからである。

図 4: モンテカルロ積分。正方形内にランダムな 点をばらまく。ランダムな点の数を増やせば、正 方形内と円内の点の数の比が、正方形と円の面 積の比になる。 {図 5: この面積を解析的に計算するのは、大変面 \atop 倒である。}

5.2 複雑な積分

もう少し複雑な積分が必要な場合を考えてみよう。ある関数 f(x_{1}, x_{2}, x_{3}, ·, x_{M}) の体積積分を考える。この体積積分は

∫_{Ω} f(x_{1}, x_{2}, x_{3}, ⋯, x_{N})dx_{1}dx_{2}dx_{3} ⋯ dx_{M} = V ⟨f⟩  (20)

となる。ここで、 V は M 次元体の領域 Ω の体積、 ⟨f⟩ はその内部の関数の平均値である。これは 3 次元体の質量を考えると簡単である。その質量は、密度の体積積分となる。これは体積に平均密度を乗じた値に等しい。当たり前の式である。このことから、式 (20) の積分の値が欲しければ、体積と平均値が分かればよいことになる。

積分を計算するために、まずは体積である。先ほどは、二次元問題であるが、はかなり複雑な形状の面積をモンテカルロ積分で計算する方法を示した。 2 次元と同じ考えで、 3 次元、 4 次元、 ⋯ 、 M 次元の体積も比較的、容易に求めることができる。分かり切った体積に、領域 Ω を包み込み、その内部にランダムに配置されたサンプル点の数を数えれば良いのである。体積の計算は簡単である。

残りは、体積内部の平均 ⟨f⟩ である。これも簡単で、領域 Ω 内部にあるサンプル点の平均より求めることができる。即ち、

⟨f⟩ ≃ \frac{1}{N_{in}} ∑_{i} f(r_{i})  領域 Ω の内部のみ  (21)

となる。ここで、 r_{i} は i 番目のサンプル点の (x_{1}, x_{2}, x_{3}, ·, x_{N}) 座標である。また、 N_{in} は領域内部のサンプル数である。この計算も簡単で、コンピューターは難なく、平均値に近似値を求めることができる。

以上より、モンテカルロ法を用いると、体積 V と平均値 ⟨f⟩ の近似値が計算できることが分かる。従って、式 (20) の近似値を求めることができる。

数値計算で有名な本「NUMERICAL RECIPES in C」 [1] には、図 6 の様な形状の重心を求める問題が出ている。このように複雑な積分が必要なときにモンテカルロ法は威力を発揮する。

図 6: 「NUMERICAL RECIPES in C」の例題。一部を切り取られたトーラスの重心を求める。

6 練習問題とレポート

6.1 問題

6.1.1 基本

本日学習したの 3 通りの方法で半径 1 の円の面積を求めよ。実際の計算は、図 8 のように、 1/4 を計算して、それを 4 倍する。

• 台形公式を使った積分

• シンプソンの公式を使った積分

• モンテカルロ積分

図 7: 円の面積

6.1.2 ちょっと応用

モンテカルロ積分を使って、図 8 の体積を計算せよ。これは、半径 1 の円柱が直交して重なった部分の体積である。

6.2 レポート

基本の 3 問とちょっと応用の 1 問の計 4 問のうち、少なくともひとつの問題のソースプログラムを提出すること。。

期限  2 月 4 日 (金)PM5:00 まで
用紙  A4
提出場所山本研究室の入口のポスト
表紙表紙を 1 枚つけて、以下の項目を分かりやすく記述すること。

授業科目名「計算機応用」課題名「積分」 5E 学籍番号氏名提出日

図 8: 直交する円柱。 3 角法で示している。

参考文献

[1] Willam H. Press et al. NUMERICAL RECIPES in C [日本語版]. 技術評論社, 1996.