📂수치해석

수치적으로 이상적분을 계산하기 위한 변수 치환 트릭

정리 ¹

$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*}$$

■

구현

다음은 예제들을 수치적으로 푼 파이썬 코드다.

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