3 数値積分

3.1 台形公式

定積分，

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

の近似値を数値計算で求めることを考える．積分の計算は面積の計算であるから，図 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*}$$ (6)

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

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

3.2 シンプソンの公式

台形公式の考え方は簡単であるが，精度はあまりよくない．そこで，よく似た考え方で精 度が良いシンプソンの公式を説明する．台形公式は，分割点の値を一次関数(直線)で近似 を行い積分を行った．要するに折れ線近似である．ここで，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})$$

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

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

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

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

$$\begin{align*} % latex2html id marker 1148 \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*}$$ (8)

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

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

$$=\frac{h}{3}f(x_{j+2})$$ 式(7)右辺第3項の積分 (10)

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

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

となる．

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

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

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

3.3 モンテカルロ法

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

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

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

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

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

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

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

$$\langle f \rangle\simeq\frac{1}{N_{\Omega}}\sum_i f(\boldsymbol{r}_i) \Omega$$ (15)

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

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

ホームページ: Yamamoto's laboratory
著者: 山本昌志