\(\sqrt{x}\)を\([0\sim b]\)の範囲で数値的に積分したいとします。
この時に問題となるのは\(\sqrt{x}\)が\(x=0\)で特異点(この場合、一階微分が発散する)を持ってしまうことです。
解決するには変数変換を行い、
\(
\displaystyle \int_0^{b}\sqrt{x} dx = \int_0^{b^{1/2}} 2x^2 dx
\)
として右辺を計算するのが良いでしょう。

端点特異点を持つ場合、数値積分は工夫しなければなりません。
ここでは\(\sqrt{x}\)を\(0\)から\(b\)まで積分することを考えます。

問題点

台形則、シンプソン則などのニュートン・コーツ型の数値積分は点と点をラグランジュ補間(多項式による補間)を行って積分します。
そのため、多項式で表せない関数を数値積分するときは著しく積分精度が落ちてしまいます。
ガウス・ルジャンドル求積法は次数が高いだけでニュートンコーツ型と同じラグランジュ補間を行うため精度が良くありません。
こういった時には台形則を組み合わせて精度評価が出来るロンバーグ積分法による方法が推奨されますが、結局は台形則なので、多項式の表現から抜け出ることはできず、被積分関数の評価回数がかなり多くなってしまいます。

ガウス求積法の一部にラゲール多項式を用いるガウス・ラゲール求積法があります。
これは
\(x=0\)で\(x^{\alpha},(0\lt \alpha)\)の振る舞い
かつ
\(x=\infty\)で被積分関数が\(e^{-x}\)で減衰する関数
に対して非常によい方法ですが、これを使う際の積分範囲は\([0\sim \infty]\)であり、今回の場合うまく用いることが出来ません。

最善の方法は二重指数関数型数値積分公式(DE公式)という積分法でしょう。
補間するのではなく、変数変換によって積分する空間を変えて積分する方法で、端点特異点がある場合に強い方法です二重指数関数型数値積分(理論とプログラム)。
端点でDE公式も使えないほどの非常に強い特異性を持つ場合、超函数法(hyperfunction method)と呼ばれる数値積分法が良いそうです。

ある程度の精度で簡単に計算するのならば、積分区間を分割します。
被積分関数が特異性を発揮するのは\(x=0\)付近ですから、積分
\(
\displaystyle \int_0^{b}\sqrt{x} dx = \int_0^{\Delta}\sqrt{x} dx + \int_{\Delta}^b\sqrt{x} dx
\)
にし、右辺第1項を低次の方法(例えば台形則)で分点数を多くとり計算します。第2項には特異性を含まないので荒い積分方法で良いでしょう。

実際に\(\sqrt{x}\)の数値積分が困難であることを数値計算で示します。
\(
\displaystyle \int_0^1 \sqrt{x} dx = \frac{2}{3} = 0.66666666666\cdots
\)
を各々の方法で計算しますと、

という結果を得ます。
いかに特異点が存在する際の数値積分が難しいか分かるでしょう。

解決法

変数変換を行います。
これにより被積分関数を多項式に置き換えることが出来ます。
変数変換により、
\(
\displaystyle \int_0^{b}\sqrt{x} dx = \int_0^{b^{1/2}} 2x^2 dx
\)
と、簡単になります。
実際、右辺を数値的に計算してみますと、

となります。短冊近似、台形則で精度が悪いのは\(x^2\)の多項式はこの二つの方法では元々出来ないものなのでこれは問題ありません。数値計算法の想定と一致する結果が導かれています。

Tweet