2 数値積分

2.1 台形公式

2.1.1 公式

定積分，

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

(1)

の近似値を数値計算で求めることを考える．積分の計算は面積の計算であるから，図 1のように台形の面積の和で近似ができるであろう．積分の範囲

を

等分した台形で近似した面積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次関数をとする．はラグランジュ補間に他ならないので，

$\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次関数

の範囲 $[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}]$ の積分で，その巾はである．区間 $[a,\,b]$ にわたっての積分は，式(7)を足し合わせればよい．ただし，と足し合わせる．

$\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)

これが，シンプソンの公式と呼ばれるもので，先ほどの台形公式よりも精度が良い．精度は，

に反比例する．

この式から，分割数は偶数でなくてはならないことがわかる．

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)

となる．ここで，

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

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