📂数値解析

数値的に広義積分を計算するための変数置換のコツ

定理 ¹

$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("")