📂수치해석

수치적 적분

정의 ¹

$f : [a,b] \to \mathbb{R}$ 가 $[a,b]$ 에서 적분가능하고 $[a,b]$ 를 $a = x_{0} < \cdots < x_{n} = b$ 와 같은 노드 포인트들로 쪼갰다고 하자.

1. 적분 오퍼레이터 $I$ 를 $\displaystyle I(f) := \int_{a}^{b} f(x) dx$ 와 같이 정의한다.
2. 적분 오퍼레이터 $I_{n}$ 을 $\displaystyle I_{n} (f) := \sum_{k=1}^{n} \int_{x_{k-1}}^{x_{k}} f(x) dx$ 와 같이 정의한다.
3. 에러 $E_{n}$ 을 $E_{n} (f) := I (f) - I_{n} ( f )$ 와 같이 정의한다.
4. $\displaystyle \lim_{n \to \infty} {{\tilde{E}_{n} (f) } \over { E_{n} (f) }} = 1$ 을 만족하는 $\tilde{E}_{n}$ 을 $E_{n}$ 에 대한 어심토틱 에러라고 한다.

설명

정적분을 계산하는 것이 수치해석의 목표 중 하나라는 점에는 그 누구도 이의를 제기하지 않을 것이다. 그런데 사실 구분구적법의 아이디어를 떠올려보면 공간을 가능한 많이 쪼개면 쪼갤수록 적분값을 잘 근사할 수 있을 것이라는 짐작은 든다. 그러니 수치해석에서 관심을 갖는 수치적 적분이라하면 구체적으로 그 오차가 얼마나 되는지, 폐구간이 아닐 땐 적분을 어떻게 해야하는지 등이라고 볼 수 있겠다.

예시

사다리꼴 룰

수치적 적분의 예로써 $f \in C^2 [a,b]$ 에 대한 사다리꼴 룰을 살펴보도록 하자. 사다리꼴 룰^{trapezoidal rule} 이란 $I(f)$ 에 대해 다음과 같이 근사값을 계산하는 방법을 말한다.

$$ I_{1}^{1} (f) := \left( {{ b - a } \over { 2 }} \right) [ f(a) + f(b) ] $$ 이는 구간 $[a,b]$ 에서 $f$ 를 리니어 인터폴레이션 해서 그 함수의 적분값으로 $I(f)$ 를 근사한 것으로 볼 수 있다. 그렇다면 실제 $I(f)$ 와 $I_{1}^{1} (f)$ 의 오차 $E_{1}^{1} (f)$ 는 어떤 $\xi \in [a,b]$ 에 대해 다음과 같이 계산된다.

폴리노미얼 인터폴레이션:

• [4] 실제 함수와의 오차: $(n+1)$번 미분가능한 $f : \mathbb{R} \to \mathbb{R}$ 와 어떤 $\xi \in \mathscr{H} \left\{ x_{0} , \cdots , x_{n} \right\}$ 에 대해 $f$ 의 폴리노미얼 인터폴레이션 $p_{n}$ 은 어떤 $t \in \mathbb{R}$ 에 대해 $\displaystyle f(t) - p_{n} (t) = {{ (t - x_{0}) \cdots (t - x_{n}) } \over { (n+1)! }} f^{(n+1)} ( \xi )$ 을 만족한다.

$$ \begin{align*} \displaystyle E_{1}^{1} (f) :=& I(f) - I_{1}^{1} (f) \\ =& \int_{a}^{b} \left[ f(x) - {{ f(b) ( x - a ) - f(a) (x - b) } \over { b - a }} \right] dx \\ =& \int_{a}^{b} \left[ f(x) - p_{1} (x) \right] dx \\ =& {{1} \over {2}} f '' ( \xi ) \int_{a}^{b} (x-a) (x-b) dx \\ =& \left[ {{1} \over {2}} f '' ( \xi ) \right] \left[ - {{1} \over {6}} (b-a)^{3} \right] \\ =& - {{1} \over {12}} (b-a)^{3} f '' ( \xi ) \end{align*} $$

컴포짓 사다리꼴 룰

이제 $[a,b]$ 를 $n$ 개의 구간으로 자르고 $\displaystyle h:= {{b-a} \over {2}} $, $x_{k} := x_{0} + h k$ 에 대해 각각의 $[x_{k-1} , x_{k} ]$ 에서 사다리꼴 룰로 적분을 근사한다고 생각해보자. 다음과 같이 근사값을 계산하는 방법을 (컴포짓) 사다리꼴 룰^{(Composite) Trapezoidal rule}이라 한다. $$ I_{n}^{1} (f) := \sum_{k=1}^{n} {{ h } \over { 2 }} \left( f (x_{k-1}) + f( x_{k} ) \right) $$ 그렇다면 실제 $I(f)$ 와 $I_{n}^{1} (f)$ 의 오차는 어떤 $\xi_{k} \in [x_{k-1}, x_{k} ]$ 들에 대해 다음과 같이 계산된다. $$ \begin{align*} \displaystyle E_{n}^{1} (f) =& I (f) - I_{n}^{1} (f) \\ =& \sum_{k=1}^{n} \left( - {{ h^3 } \over { 12 }} f '' ( \xi_{k} ) \right) \\ =& - {{ h^3 n } \over { 12 }} \left( {{ 1 } \over { n }} \sum_{k=1}^{n} f '' ( \xi_{k} ) \right) \end{align*} $$ 여기서 $\displaystyle \left( {{ 1 } \over { n }} \sum_{k=1}^{n} f '' ( \xi_{k} ) \right)$ 은 $f’’ ( \xi_{k} )$ 들의 평균이므로, 중간값 정리에 의해 $\displaystyle \min_{x \in [a,b]} f ''(x) \le f ''( \xi ) \le \max_{x \in [a,b]} f ''(x)$ 를 만족하는 어떤 $\xi \in [a,b]$ 가 존재해야한다. 따라서 $$ E_{n}^{1} (f) = - {{ (b-a) h^2 } \over {12}} f '' (\xi ) $$ 이렇게 에러를 직접 식으로 보니까 적분 구간 $[a,b]$ 가 넓어질수록 오차가 커지며 $h$ 가 작을수록 정확해지는 것을 알 수 있다. 그러나 $f’’ ( \xi ) $ 때문에 구체적으로 얼마나 틀렸는지는 알 수가 없다. 이를 극복하기 위해 다음과 같이 어심토틱 에러를 구한다. $$ \begin{align*} \lim_{n \to \infty} {{ E_{n}^{1} (f) } \over { h^2 }} =& \lim_{n \to \infty} {{1} \over {h^2}} \sum_{k=1}^{n} \left( - {{ h^3 } \over { 12 }} f '' ( \xi_{k} ) \right) \\ =& - {{ 1 } \over { 12 }} \lim_{ n \to \infty} \sum_{k=1}^{n} h f '' ( \xi_{k} ) \\ =& - {{ 1 } \over { 12 }} \int_{a}^{b} f ''(x) dx \\ =& - {{ 1 } \over { 12 }} [ f '(b) - f '(a) ] \end{align*} $$ 에러 $E_{n}^{1} (f)$ 만 알 땐 정리에 의해서 $f’’ ( \xi )$ 가 존재하는 건 알지만 구체적으로 수치를 댈 수는 없었는데, 이제 $$ \tilde{E}_{n}^{1} (f) := \lim_{n \to \infty} E_{n}^{1} (f) $$ 이라고 하면 충분히 큰 $n$ 에 대해 $$ E_{n}^{1} (f) \approx - {{ h^2 } \over { 12 }} [ f '(b) - f '(a) ] = \tilde{E}_{n}^{1} (f) $$ 임을 알 수 있고 $f ' (a) , f '(b)$ 를 계산하는 것은 비교적 간단한 일이므로 상당히 쉽고 정확하게 오차를 파악할 수 있다. 게다가 여기까지 오면 오차가 그정도 되는구나 하는 정도로 끝날 이유도 없다. 계산을 하기도 전부터 이렇게 오차가 날 것이라고 예측할 수 있다면 아예 처음부터 오차를 보정할 수도 있다. $$ I(f) - I_{n}^{1} (f) \approx \tilde{E}_{n}^{1} (f) \implies I(f) \approx I_{n}^{1} (f) + \tilde{E}_{n}^{1} (f) $$

커렉티드 사다리꼴 룰

이렇게 $I(f)$ 를 계산하는 오퍼레이터로써 $CT_{n}$ 을 다음과 같이 정의하자. $$ CT_{n} (f) := h \left[ {{1} \over {2}} f (x_{0} ) + f (x_{1} ) + \cdots + f (x_{n-1} ) + {{1} \over {2}} f (x_{n} ) \right] - {{ h^2 } \over { 12 }} [ f '(b) - f '(a) ] $$

이 방법을 (커렉티드) 사다리꼴 룰(^{corrected) Trapezoidal rule}이라 한다. [ NOTE: 위에서 사다리꼴 룰이 개선되는 과정을 봤다면 알겠지만 굳이 이들을 구분할 필요는 없다. 그냥 사다리꼴 룰이라고 해도 의미는 다 통한다. ]

구현

실제로 $\displaystyle \int_{0}^{\pi} e^{x} \cos (x) dx = -{{e^{\pi} + 1} \over {2}} \approx - 12.0703$ 을 $I_{n}^{1}$ 과 $CT_{n}$ 으로 계산해본 결과는 다음과 같다.

첫번째 열은 $n$, 두번째 열은 $\left| I(f) - I_{n}^{1} (f) \right|$, 세번째 열은 $\left| I(f) - CT (f) \right|$ 를 나타낸다. 컴포짓 사다리꼴 룰은 노드를 두배로 늘릴 때마다 오차가 $\displaystyle {{1} \over {4}}$ 배로 줄어들고, 커렉티드 사다리꼴 룰은 $\displaystyle {{1} \over {16}}$ 배로 줄어들고 있음을 확인할 수 있다. 단순히 수치만 보아도 커렉티드 사다리꼴 룰이 훨씬 좋은 걸 알 수 있지만, 특히 좋은 점은 그 속도가 엄청나게 빠르다는 것이다.

다음은 파이썬으로 작성된 예제 코드다.

import numpy as np
def f(x) :
    return np.exp(x)*np.cos(x)
 
def d(f,x,tol=10**(-8)) :
    h=1
    d1=1
    d2=0
 
    while(abs(d1-d2)>tol) :
        d1=d2
        d2=(f(x+h)-f(x))/h
        h=h/2
 
    return(d2)
 
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 CT(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)-(h**2)*(d(f,b)-d(f,a))/12)
 
 
TV=-(np.exp(np.pi)+1)/2
print("True Value : %2.4f" % TV)
 
print("|   n |  I(f) Err | CT(f) Err |")
print("="*30)
for n in range(1,10) :
    print("| %3d |  %8.2e |  %8.2e |" % (2**n,abs(TV-I(f,0,np.pi,2**n)),abs(TV-CT(f,0,np.pi,2**n))))