数値積分

 作成日：2013-09-15
最終更新日：

基礎的な積分則

数値積分は研究が進んだ分野である。 分点数 `n` があたえられたときに `n` 個の関数値と重みを使い、 積和計算によって積分値を求める方法を、ここでは基礎的な積分則と呼ぶ。

台形則

区間 `[a, b]` を `n` 等分する。 `h = (b-a)/ n` として、それぞれの区間 `[x_0, x_1], [x_1, x_2], ......, [x_(n-1), x_n]` の各台形の面積の和を求める。 これが台形(Trapezoidal)則である。

`int_a^b f(x)dx = h / 2 [sum_(i=1)^n{f(x_(i-1)) + f(x_i)}]`

台形則の計算

実際の計算を下記の台形則で実行することができる。

下限 a：
上限 b：
関数 f：
分点数 n：
積分値I：

二重指数関数型数値積分公式

積分範囲が端点で特異性があるときなどに有効なのが、 二重指数関数型数値積分公式 (DE 公式) である。 英米では tanh-sinh 求積法(quadrature) と呼ばれることもある。 求める定積分の下限を `a`, 上限を `b` とするとき、`x = a ` から`b`までの定積分を `I(a, b)` と記す。

`I(-1, 1) = int_-1^1 f(x) dx ~= sum_(k=-oo)^oo w_k f(x_k)　~= sum_(k=-N)^N w_k f(x_k)`

`N` は適当に十分大きな数である。 分点 `x_k` と 重み `w_k` は任意のきざみ幅 `h` を用いて次の式で与えられる。

`x_k = tanh (1/2 pi sinh kh)`
`w_k = (1/2 h pi cosh kh) / (cosh^2 (1/2 pi sinh kh))`

任意の区間`I(a, b)`のときの分点と重みは次のようになる。

`x_k = (b - a)/2 tanh (1/2 pi sinh kh) + (b + a) / 2`
`w_k = ((b - a)/2) (1/2 h pi cosh kh) / (cosh^2 (1/2 pi sinh kh))`

上下限 N：
刻み幅 h：
関数 f：
積分値I：

例１．円の面積

`I_1 = 2 int_-1^1 sqrt(1-x^2)dx`

`y = sqrt(1-x^2)` は円の方程式の y が正の部分だから、単位円の面積の 1/2 である。 したがって上記積分は単位円の面積であり、厳密解は `pi` である。

例２．一方の端点で特異点を持つ例

`I_2 = 1/2 int _-1^1 dx / (sqrt(x + 1))`

被積分関数は `x -> -1+0 ` で発散する。一方

` 1/2 int _a^b dx / (sqrt(x + 1)) = [sqrt(x+1)]_a^b`

であるから、`I_2 = sqrt(2)` である。

上記に入れて試してみると、`N` が 4以上では積分値が Infinity と表示され、計算できない。 対策は追って考える。

例３．両方の端点で特異点を持つ例

`I_3 = int _-1^1 dx / sqrt(1-x^2)`

`I_3 = pi` である。 やはり、`N=4` 以上では積分値が Infinity となる。

原因と対策

積分値が Infinity と表示される原因と対策を述べる。以下は、 長田直樹先生のホームページの資料に基づく。

Infinity となる原因は、被積分関数が`1+x`や`1-x`を因数にもつと、 `1 +- tanh (pi/2 sinh t)` の計算において桁落ちが起こるためという。

桁落ちを避ける対策としては式を変形することである。具体的には、双曲線関数で成り立つ下記の等式

`1 + tanh x = exp x / cosh x, 1 - tanh^2 x = 1 / (cosh^2 x) , 1 - tanh x = (1 - tanh^2 x) / (1 + tanh x) = (1 / (cosh^2 x)) * (cosh x / exp x) = 1 / (cosh x exp x)`

を利用する。すると、例２の積分値は次のように書ける。

`I_2 = 1/2 sum_(k=-m)^n (1/2 h pi cosh kh) / (cosh^2 (1/2 pi sinh kh)) sqrt (exp(1/2 pi sinh kh) cosh (1/2 pi sinh kh))`

例３の積分値は次の通りとなる。

`I_3 = sum_(k=-m)^n (1/2 h pi cosh kh) / (cosh^2 (1/2 pi sinh kh)) cosh (1/2 pi sinh kh) = sum_(k=-m)^n (1/2 h pi cosh kh) cosh (1/2 pi sinh kh) `

これを計算するフォームを用意した。下記のフォームはデフォルトでは例２を計算するもので、 `h = 0.25` にするとかなりの桁まで合っている。m と n は長田先生のリンクにある下記の方法で、非常に小さな数 `epsilon` を与えて自動的に求めている。

`|w_m f_m(-mh)| + |w_m f_m(-(m+1)h)| < epsilon`
`|w_n f_n(nh)| + |w_n f_n((n+1)h)| < epsilon`

また、 square(x) は x の二乗を計算する関数である。

`epsilon` =
n
`Sigma`
k = -m
刻み幅 h：
下限ｍ：
上限ｎ：
積分値I：

リンク集

まりんきょ学問所 ＞ JavaScript 手習い ＞ 数値積分