数値的に広義積分を計算するための変数置換のコツ
定理 1
$0 < a < b < \infty$ としよう。
- [1]: $ 0 < p < 1$ ならば $$\int_{0}^{b} {{ f(x) } \over {x^{p} }} dx = \int_{0}^{{{ 1 } \over { 1-p }} b^{1-p} } f \left( \left[ ( 1- p ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm$$
- [2]: $ 1 < p$ ならば $$\int_{a}^{ \infty } {{ f(x) } \over {x^{p} }} dx = \int_{0}^{{{ 1 } \over { p-1 }} a^{1-p}}f \left( \left[ ( p-1 ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm$$
説明
不適切積分はいつも頭痛の種だけど、解くのを諦めるわけにはいかない。積分が通常難しい理由は、どんなトリックを使っても、それらのトリックが結局新しい問題を作り出すからだ。幸い数値積分では、関数を正確に持っていれば、詳細なプロセスがなくてもとにかく計算を進めることができる。代表的な方法は代入で、関数が特異であっても、$x^{p}$ のような項を分母から引くことができれば、簡単に解決できる。
例
回避積分
例として、次の不適切積分を考えてみよう。 $$ \int_{0}^{1} {{e^{-x}} \over {\sqrt{x}}} dx \approx 1.493648265625 $$ $\displaystyle m(x) := \int_{0}^{x} {{1} \over {\sqrt{t} }} dt = 2 \sqrt{x}$ とすると $\displaystyle \sqrt{x} dm = dx$ となるので $$ \begin{align*} \int_{0}^{1} {{e^{-x}} \over {\sqrt{x}}} dx =& \int_{m=0}^{2} {{e^{-x}} \over {\sqrt{x}}} \sqrt{x} dm \\ =& \int_{0}^{2} e^{-x} dm \\ =& \int_{0}^{2} e^{- m^2 / 4} dm \end{align*} $$ このようなトリックを使って閉区間に対する積分を作れば、積分技術を適用できる。実際に台形則を使って数値を計算してみると、次のような結果が得られる。
速度の向上
一方で、より興味深い例として、次の場合を見てみよう。 $$ \int_{0}^{1} e^{-x} \sqrt{x} dx \approx 0.378944691641 $$ 最初に思い浮かぶのは「特に代入が必要ない」という点だろう。$e^{-x} \sqrt{x}$ は$[0,1]$ で特異性を持たないので、そのまま 台形則 を適用できる。しかし、$x=0$ で微分不可能であるため、エラーのバウンドがなくなり、収束速度が遅くなる。同様に $\displaystyle m(x) := \int_{0}^{x} {{1} \over {\sqrt{t} }} dt = 2 \sqrt{x}$ とすると $\displaystyle \sqrt{x} dm = dx$ となるので $$ \begin{align*} \int_{0}^{1} e^{-x} \sqrt{x} dx =& \int_{m=0}^{2} e^{-x} \sqrt{x} \sqrt{x} dm \\ =& \int_{0}^{2} e^{-x} x dm \\ =& \int_{0}^{2} e^{- m^2 / 4} {{m^2} \over {4}} dm \end{align*} $$ 見ての通り、代入積分を行うと、$m=0$ の特異性がなくなる。同じように台形則 を使って数値を計算してみると、次のような結果が得られる。左がそのまま計算したもので、右が代入を行って計算したものだ。
代入なしでも収束はするが、代入を通じた計算結果が怖いほど真値に収束するのに対して、代入なしで計算した値は収束速度が遅いことが確認できる。数値解析の特徴として速度を重視するため、このようなロスはできるだけ減らすのが良い。いつか選ばなくてはならないときが来たら、コードの簡潔さ、美しさよりも、数学者の細かくて丁寧なトリックを選んだ方が良いだろう。
証明
戦略:実際の計算はコンピュータに任せるため、閉区間での積分に変換するだけで良い。方法は本質的に同じだ。
[1]
$$ m := \int_{0}^{x} {{ 1 } \over { t^{p} }} dt = {{ 1 } \over { 1-p }} x^{1-p} $$ のようにすると $\displaystyle {{ dm } \over { dx }} = x^{-p} \iff x^{p} dm = dx$ となるので $$ \begin{align*} \int_{0}^{b} {{ f(x) } \over {x^{p} }} dx =& \int_{m=0}^{{{ 1 } \over { 1-p }} b^{1-p}} {{ f(x) } \over {x^{p} }} x^{p} dm \\ =& \int_{0}^{{{ 1 } \over { 1-p }} b^{1-p}}f \left( \left[ ( 1- p ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm \end{align*} $$
■
[2]
$$ m := \int_{x}^{\infty} {{ 1 } \over { t^{p} }} dt = - {{ 1 } \over { 1-p }} x^{1-p} $$ のようにすると $\displaystyle {{ dm } \over { dx }} = - x^{-p} \iff - x^{p} dm = dx$ となるので $$ \begin{align*} \int_{a}^{ \infty } {{ f(x) } \over {x^{p} }} dx =& \int_{m= - {{ 1 } \over { 1-p }} a^{1-p}}^{0} {{ f(x) } \over {x^{p} }} ( - x^{p} ) dm \\ =& \int_{0}^{{{ 1 } \over { p-1 }} a^{1-p}}f \left( \left[ ( p-1 ) m \right]^{{{ 1 } \over { 1-p }}} \right) dm \end{align*} $$
■
実装
次は、例を数値的に解くためのPythonのコードだ。
import numpy as np
def f1(x) :
return np.exp(-x)
def f2(x) :
return np.exp(-x)*np.sqrt(x)
def f3(x) :
return np.exp(-x)*x
def I(f,a,b,n) :
h=(b-a)/n
x=np.linspace(a,b,n+1)
x=f(x)
x[0]=x[0]/2
x[-1]=x[-1]/2
return(h*sum(x))
def CVI(f,a,b,n) :
b=2*np.sqrt(b)
h=(b-a)/n
x=np.linspace(a,b,n+1)
x=f((x**2)/4)
x[0]=x[0]/2
x[-1]=x[-1]/2
return(h*sum(x))
TV1=1.49364826562485405079
TV2=0.37894469164098470380
print("")
print("True Value 1 : %2.12f" % TV1)
print("| n | In | I-In |")
print("="*32)
for n in range(3,10) :
print("| %3d | %8.6f | %9.2e |" % (
2**n,
CVI(f1,0,1,2**n),
TV1-CVI(f1,0,1,2**n)))
print("")
print("True Value 2 : %2.12f" % TV2)
print("| n | In | I-In |")
print("="*32)
for n in range(3,10) :
print("| %3d | %8.6f | %9.2e |" % (
2**n,
I(f2,0,1,2**n),
TV2-I(f2,0,1,2**n)))
print("")
print("True Value 2 : %2.12f" % TV2)
print("| n | In | I-In |")
print("="*32)
for n in range(3,10) :
print("| %3d | %8.6f | %9.2e |" % (
2**n,
CVI(f3,0,1,2**n),
TV2-CVI(f3,0,1,2**n)))
print("")
Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p305~306. ↩︎