4 数値積分

4.1 台形公式

定積分、

$\displaystyle S=\int_a^bf(x)dx$ (12)

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

\begin{equation*}\begin{aligned}T&=h\frac{f(a)+f(a+h)}{2}+ h\frac{f(a+h)+f(a+2h)...
...{2}\sum_{j=0}^{N-1}\left[f(a+jh)+f(a+(j+1)h)\right] \end{aligned}\end{equation*}

となる。これが数値積分の台形公式である。なんのことはない、積分を台形の面積に置き 換えているだけである。
図 2: 積分と台形の面積の比較
\includegraphics[keepaspectratio, scale=1.0]{figure/trapezoidal.eps}

4.2 シンプソンの公式

台形公式の考え方は簡単であるが、精度はあまりよくない。そこで、よく似た考え方で精 度が良いシンプソンの公式を説明する。台形公式は、分割点の値を一次関数(直線)で近似 を行い積分を行った。要するに折れ線近似である。ここで、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次関数で近似する
\includegraphics[keepaspectratio, scale=1.0]{figure/simpson.eps}

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

\begin{equation*}
% latex2html id marker 1716
\begin{aligned}\text{式(\ref{eq:si...
...i^2}{2}+2h^2\xi\right]_0^{2h} &=\frac{h}{3}f(x_j) \end{aligned}\end{equation*}

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

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

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

$\displaystyle \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\}$ (17)

となる。

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

\begin{equation*}\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}\end{equation*}

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

この式から、分割数$ N$は偶数でなくてはならないことがわかる。

4.3 モンテカルロ法

積分の境界が複雑な場合、乱数を使うモンテカルロ積分が適している。例えば、関数 $ f(x_1,x_2,x_3,\cdot,x_M)$の体積積分を考える。この体積積分は

$\displaystyle \int_\Omega f(x_1,x_2,x_3,\cdots,x_M)dx_1dx_2dx_3\cdots dx_M= V_{\Omega}\langle f \rangle$ (19)

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

積分を計算するために、まずは体積である。これは体積の計算が容易な単純な形状の内部 に、領域$ \Omega$を包み込み、その内部にランダムに配置されたサンプル点の数を数えれ ば良いのである。単純な形状内部に配置されたランダムな点の数を$ N$とする。そして、 その内部にある積分領域$ \Omega$に含まれる点の数を $ N_{\Omega}$とする。さらに、単純 な形状の体積を$ V_r$、領域$ \Omega$のそれを $ V_{\Omega}$とすると、

$\displaystyle V_{\Omega}\simeq\frac{N_{\Omega}}{N}V_r$ (20)

の関係がある。右辺はコンピューターにより容易に計算できる。ランダムな点の数$ N$が 多くなればなるほど、近似の精度は良くなる。

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

$\displaystyle \langle f \rangle\simeq\frac{1}{N_{\Omega}}\sum_i f(\boldsymbol{r}_i)$   領域$ \Omega$の内部のみ (21)

となる。ここで、 $ \boldsymbol{r}_i$$ i$番目のサンプル点の $ (x_1,x_2,x_3,\cdot,x_M)$座標で ある。また、 $ N_{\Omega}$は領域内部のサンプル数である。この計算も簡単で、コンピューター は難なく、平均値に近似値を求めることができる。

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


ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成17年3月1日


no counter