Subsections

4 数値積分

4.1 台形公式

定積分,

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

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

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

となる.これが数値積分の台形公式である.なんのことはない,積分を台形の面積に置き 換えているだけである.
図 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)$ はラグランジュ補間に他ならないので,

$\displaystyle P(x)$ $\displaystyle =\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})$    
  $\displaystyle \hspace{60mm} +\frac{(x-x_j)(x-x_{j+1})}{(x_{j+2}-x_j)(x_{j+2}-x_{j+1})}f(x_{j+2})$ (16)

となる.図3に示すとおりである.
図 3: 元の関数を区間 $ [x_j,x_{j+2}]$ を2次関数で近似する
\includegraphics[keepaspectratio, scale=1.0]{figure/simpson.eps}

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

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

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

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

となる.以上より,近似した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\}$ (20)

となる.

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

\begin{align*}\begin{aligned}S&=\frac{h}{3}\left\{f(x_0)+4f(x_1)+f(x_2)\right\} ...
...mm}\left.\cdots+2f(x_{N-2})+4f(x_{N-1})+f(x_N)\right\} \end{aligned}\end{align*} (21)

これが,シンプソンの公式と呼ばれるもので,先ほどの台形公式よりも精度が良い.精度 は,$ 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$ (22)

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

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

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

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

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

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

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

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


ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
2006-02-08


no counter