積分法

不定積分

不定積分の表を上記公式集から抜粋する :

不定積分の表

なお、`I_A, I_B` は次の積分を表わす :

`I_A = int (dx)/(x^2+c)`
`= { (1/sqrt(c)arctan{:x/sqrt(c):}, (c gt 0)), (1/(2sqrt(abs(c))) log abs((x-sqrt(abs(c)))/(x+sqrt(abs(c)))) , (c lt 0)) :}`
`I_B = int (dx)/(ax^2+bx+c)`
` = {(1/sqrt(b^2-4ac) log abs((2ax+b-sqrt(b^2-4ac))/(2ax+b+sqrt(b^2-4ac))), (b^2 gt 4ac)), (2/sqrt(4ac-b^2) arctan{: (2ax+b)/sqrt(4ac-b^2):}, (b^2 lt 4ac)) , (-2/(2ax+b), (b^2=4ac)) :}`

数値積分における基礎的な積分則

不定積分が求められない関数でも、定積分では厳密解が求められる場合がある。 たとえば、複素関数論を利用する方法がある。 複素関数論が適用できない場合でも、数値積分により精度の高い数値解を得ることができる。 このための数値解析は研究が進んだ分野である。 分点数 `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)}]`

台形則の計算

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

以前は関数 eval を使っていたが、セキュリティ上の問題などで使用は推奨されていない。 その代わりに、Function を使うことにした。eval から Function へを参照。

また、with 構文も推奨されなくなった。そのため、JavaScript にある関数を使うときは、 上記の Math.exp(x) のように、Math モジュールを明示的に指定してほしい。

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

積分範囲が端点で特異性があるときなどに有効なのが、 二重指数関数型数値積分公式 (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))`

JavaScript にある関数を使うときは、台形公式のときと同様に、 上記の Math.exp(x) のように、Math モジュールを明示的に指定してほしい。

例１．円の面積

`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_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`

また、k = -m の欄は、以前は square(Math.cosh (0.5 * Math.PI * Math.sinh (k * h) ) ) としていた（ square(x) は x の二乗を計算する関数）としていたが、 ここだけでしか使わない二乗を計算する関数に新たな sqaure という名前をつける必要はない。 同じ意味となるように、無名関数を用いた表記 (function(t){return t*t;})(Math.cosh (0.5 * Math.PI * Math.sinh (k * h) ) ) に改めた。

関連

• 複素積分

リンク集

• 長田直樹先生のホームページ＞ お話 : 数値解析（理系への数学）から 数値積分(2008年8月号)

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