2 数値積分

2.1 台形公式

2.1.1 公式

定積分,

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

の近似値を数値計算で求めることを考える.積分の計算は面積の計算であるから,図 1のように台形の面積の和で近似ができるであろう.積分の範 囲$ [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*} (2)

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

2.1.2 プログラム例

以下の定積分を台形公式を使い計算するプログラム例をリスト1にしめす.
   1 #include <stdio.h>
   2 #include <math.h>
   3 #define N 1000                      // 分割数
   4 
   5 double f(double x);                 // プロトタイプ宣言
   6 
   7 //======= メイン関数 ==================
   8 int main(void)
   9 {
  10   double a, b, h;
  11   double t=0;
  12   int j;
  13 
  14   a=0;                               // 積分の左端
  15   b=2*M_PI;                          // 積分の右端
  16 
  17   h=(b-a)/N;
  18 
  19   for(j=0; j<N; j++){
  20     t += f(a+j*h)+f(a+(j+1)*h);
  21   }
  22   t *= h/2;
  23 
  24   printf("Integral = %e\n", t);
  25   
  26   return 0;
  27 }
  28 
  29 
  30 //======= 関数 ========================
  31 double f(double x)
  32 {
  33   double y;
  34 
  35   y=sin(x)*sin(x);
  36 
  37   return y;
  38 }

2.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})$ (3)

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

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

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

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

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

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

となる.

これは,ある区間 $ [x_j,\,x_{j+2}]$の積分で,その巾は$ 2h$である.区間$ [a,\,b]$にわ たっての積分$ S$は,式(7)を足し合わせればよい.ただし, $ 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*} (8)

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

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

2.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$ (9)

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

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

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

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

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

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

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

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


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


no counter