[수치해석] Numerical integration (2) - Simpson's rule, Romberg integration, Adaptive quadrature

이번에는 수치해석 적분 중 하나인 Simpson's rule에 대해서 공부하겠습니다.

(mid point rule과 trapezoidal rule에 대해서는 이전 글을 참고해주세요.)

 

 

1. Simpson's rule

 

 

$\text{Rectangular (mid point) rule error : }R(f)=\frac{1}{24}h_{i}^{3}f^{''}(y_{i})+\theta (h_{i}^{5})$

$\text{Trapezoidal rule error : }T(f)=-\frac{1}{12}h_{i}^{3}f^{''}(y_{i})+\theta (h_{i}^{5})$

$y_{i}=\frac{x_{i}+x_{i+1}}{2}\text{ : mid point}$

$h_{i}=x_{i+1}-x_{i}$

 

 

위의 두 rule의 error를 알고 있는 상태에서, 이 rule을 조합해서 아래의 simpson's rule을 구합니다.

 

$\text{Simpson's rule : }S(f)=\frac{2}{3}R(f)+\frac{1}{3}T(f)=\theta (h_{i}^{5})$

 

$\int_{x_{i}}^{x_{i+1}}f(x)dx=\frac{2}{3}(h_{i}f(y_{i}))+\frac{1}{3}\left ( \frac{h_{i}}{2}(f(x_{i})+f(x_{i+1})) \right )=\frac{1}{6}\left \{ f(x_{i})+4f(y_{i})+f(x_{i+1}) \right \}+\theta(h_{i}^{5})$

 

 

이 식으로 $\theta(h^{5})$의 order of accuracy를 달성할 수 있지만 문제가 있습니다.

위 식에서 $f(y_{i})$가 사실은 모르는 값이라는 것입니다.

mid point rule을 사용하면 발생하는 것이고, 결국 mid point를 만들기 위해 원래의 $h_{i}$간격을 포기하고 $2h_{i}$를 택하면서 $x_{i},x_{i+1},x_{i+2}$가 있을 때 $x_{i+1}$을 mid point로 잡게 됩니다.

 

두 점 사이의 간격이 2배로 증가했기 때문에 5차 order라고 해도 오차가 증가할 수 밖에 없고, 실제로 simpson's rule을 적용했을 때 오차가 더 큰 경우도 있다고 합니다.

 

 

 

2. Romberg intergration (only works for uniform grid)

 

integration scheme에 Richardson extrapolation을 적용한 형태

 

일단 처음에는 Trapezoidal integration으로 시작하겠습니다. trapezoidal은 locally third order error를 가지고 있지만, globally second order입니다. (mean value theorem으로 증명 가능)

 

$I=\int_{a}^{b} f(x)dx = \frac{h}{2}\left [ f(a)+2\sum_{i=1}^{n-1}f(x_{i})+f(b)\right ] +c_{1}h^{2}+c_{2}h^{4}+\cdots$

$c_{1}h^{2}=-\frac{(b-a)}{12}h^{2}f(\overline{x})\rightarrow $ 이전 포스팅 참고

 

우변의 수치해석적으로 구한 integral 값($\frac{h}{2}\left [ f(a)+2\sum_{i=1}^{n-1}f(x_{i})+f(b)\right ]$)를 $\widetilde{I_{1}}$라고 할 때, 식을 다음과 같이 바꿀 수 있습니다.

 

$\widetilde{I_{1}}=I-c_{1}h^{2}-c_{2}h^{4}-c_{3}h^{6}+\cdots (1)$

이 때, $h \rightarrow \frac{h}{2}$로 바꾸면

 

$\widetilde{I_{2}}=I-c_{1}\left(\frac{h}{2}\right)^{2}-c_{2}\left(\frac{h}{2}\right)^{4}-c_{3}\left(\frac{h}{2}\right)^{6}+\cdots (2)$

 

우리의 목표는 order of accuracy를 증가시키는 것이므로 (1),(2) 식에서 $h^{2}$를 제거하는 방향으로 식을 변형합니다.

 

$4 \times (2) -(1)=3I+\frac{3}{4}c_{2}h^{4}+\cdots$

$I=\frac{4 \widetilde{I_{2}}-\widetilde{I_{1}}}{3}-\frac{1}{4}c_{2}h^{4}-\cdots (3)$

 

$\text{If }h\rightarrow \frac{h}{4}$

$\widetilde{I_{3}}=I-c_{1}\left(\frac{h}{4}\right)^{2}-c_{2}\left(\frac{h}{4}\right)^{4}-c_{3}\left(\frac{h}{4}\right)^{6}+\cdots (4)$

 

 

 

각 점 사이의 간격을 4분의 1로 줄인다면 나오는 결과는 (3)과 같고, 위에서 했던 것과 마찬가지로 

$4\times (4)-(2)=3I+\frac{3}{64}h^{4}+\cdots (5)$

 

다음과 같이 구할 수 있습니다.

여기서 구한 (3)식과 (5)식은 공통적으로 $\theta(h^{4})$입니다. 따라서 (3)과 (5)식을 조합하여 $h^{4}$\항을 제거하고 $h^{6}$항만 남겨서 order of accuracy를 증가할 수 있습니다.

 

이렇게 하는 것의 장점은, 에러를 줄이기 위해 $\widetilde{I}$에서 $h$를 계속 작게 만드는 것보다($h\rightarrow \frac{h}{128}$로 줄이는 것처럼) 적분 결과를 조합해서 위 그림의 $(iii)$처럼 order를 낮추는 것이 더 효과적일 수 있다는 점입니다.

 

 

이렇게 romberg integration을 적용한 예시 :

 

 

3. Adaptive quadrature

 

 

이 방법은, 더 그리드를 촘촘하게 해야할 부분은 그리드를 더 만들고, 덜 필요한 부분은 중단하는 방법을 통해 효율적으로 계산하는 방법을 말합니다.

 

위의 Romberg 방식은 항상 uniform grid로 해야한다는 전제가 있는데 이 방법은 그 전제에서 벗어날 수 있습니다.

점 사이의 간격 $h$는 gradient integral에 따라서 결정됩니다.

 

 이 방법의 장점은 error tolerance를 가질 수 있다는 점이고, 단점은 error analysis가 필요하다는 점입니다.

 

$\text{error }e = \left | I-\int_{b}^{a} f(x) dx \right | \leq \epsilon \text{(tolerance)}$

 

위 식에서 $I$는 numerical한 방법으로 얻은 적분값이고 $\int_{b}^{a} f(x) dx$은 exact한 값입니다.

어느정도의 에러를 허용할 것인지를 결정하는 것이 $\epsilon$, tolerance입니다.

 

 

그렇다면 이 방법이 구체적으로 어떤 식으로 진행되는지 보여드리겠습니다.

 

$\text{Simpson's rule for }[x_{i},x_{i+1}]\text{ and }h_{i}=x_{i+1}-x_{i}$에서 시작합니다. 즉, midpoint를 가지고 scheme을 만드는 것부터 시작입니다.

 

simpson's rule 식은 아래와 같은데 $S_{i}^{(1)}$은 첫 번째 mid point인 $x_{i+\frac{1}{2}}$을 기준으로 세운 식이고, 그 아래의 식은 $S_{i}^{(1)}$식에서 $x_{i+\frac{1}{2}}$을 중심으로 양 쪽에 또 mid point를 만들어서 $S_{i}^{(2)}$식을 만든 것입니다. (위 그림 참고)

 

$S_{i}^{(1)}=\frac{h_{i}}{6}[f(x_{i})+4f(x_{i+\frac{1}{2}})+f(x_{i+1})]$

$S_{i}^{(2)}=\frac{h_{i}}{12}[f(x_{i})+2f(x_{i+\frac{1}{2}})+4f(x_{i+\frac{1}{4}})+4f(x_{i+\frac{3}{4}})+f(x_{i+1})]$

 

$S_{i}^{(2)}$식을 유도하는 방법 : $h_{i}\rightarrow \frac{h_{i}}{2}$로 바꾸고 $[x_{i},x_{i+\frac{1}{2}}],[x_{i+\frac{1}{2}},x_{i+1}]$ 두 부분으로 나눠서 각각에 대해 simpson's rule을 적용하면 됩니다.

$\frac{h_{i}}{2}\times \frac{1}{6} \left\{ f(x_{i})+4f(x_{i+\frac{1}{4}})+f(x_{i+\frac{1}{2}})+f(x_{i+\frac{1}{2}})+4f(x_{i+\frac{3}{4}})+f(x_{i+1})\right\}$

이 식을 정리하면 위와 같이 나옵니다.

 

$I_{i}\text{ : exact integral in }[x_{i},x_{i+1}]$

$I_{i}-S_{i}^{(1)}=cf^{(4)}(x_{i+\frac{1}{2}})h_{i}^{5}$

$I_{i}-S_{i}^{(2)}=cf^{(4)}(x_{i+\frac{1}{4}})\left( \frac{h_{i}}{2}\right)^{5}+cf^{(4)}(x_{i+\frac{3}{4}})\left(\frac{h_{i}}{2}\right)^{5}$

 

여기서 $f^{(4)}$는 f를 4번 미분한 것을 의미하고 $I_{i}$는 exact한 적분값입니다.

 

higher order를 위해 $I_{i}-S_{i}^{(1)}$식과 $I_{i}-S_{i}^{(2)}$식을 cancel하고자하는데, $f^{(4)}(x_{i+\frac{1}{2}})$과 $f^{(4)}(x_{i+\frac{1}{4}})$,$f^{(4)}(x_{i+\frac{3}{4}})$이 저마다 다르기 때문에 하나로 통일해줘야 합니다.

 

$(+)\left\{\begin{matrix}
f^{(4)}(x_{i+\frac{1}{4}})=f^{(4)}(x_{i+\frac{1}{2}})-\frac{h_{i}}{4}f^{(5)}(x_{i+\frac{1}{2}})+\left ( \frac{h_{i}}{4} \right )^{2}f^{(6)}(x_{i+\frac{1}{2}})\\ 
f^{(4)}(x_{i+\frac{3}{4}})=f^{(4)}(x_{i+\frac{1}{2}})+\frac{h_{i}}{4}f^{(5)}(x_{i+\frac{1}{2}})+\left ( \frac{h_{i}}{4} \right )^{2}f^{(6)}(x_{i+\frac{1}{2}})
\end{matrix}\right.$

이 두 식을 더해주면 $f^{(4)}(x_{i+\frac{1}{4}})+f^{(4)}(x_{i+\frac{3}{4}})=2f^{(4)}(x_{i+\frac{1}{2}})+\frac{h_{i}^{2}}{8}f^{(6)}(x_{i+\frac{1}{2}})+\cdots$이고

 

이제 $f^{(4)}(x_{i+\frac{1}{4}})+f^{(4)}(x_{i+\frac{3}{4}})$식을 $f^{(4)}(x_{i+\frac{1}{2}})$에 대해서 표현할 수 있으므로

$I_{i}-S_{i}^{(2)}=c\left( \frac{h_{i}}{2}\right)^{5}\left[ 2f^{(4)}(x_{i+\frac{1}{2}})+\theta(h_{i}^{2})\right]$

 

$I_{i}-S_{i}^{(1)}-(I_{i}-S_{i}^{(2)})=S_{i}^{(2)}-S_{i}^{(1)}=\frac{15}{16}\left( h_{i}^{5}f^{(4)}(x_{i+\frac{1}{2}}) \right)+\cdots$

 

이제 이 식을 $I_{i}-S_{i}^{(2)}=c\frac{h_{i}^{5}}{16} f^{(4)}(x_{i+\frac{1}{2}}) + \theta$를 변형하는 데에 쓸 수 있습니다.

 

$I_{i}-S_{i}^{(2)}=\frac{1}{15}(S_{i}^{(2)}-S_{i}^{(1)})$

 

두 점의 간격을 $\frac{1}{4}$배로 해서 구한 $S_{i}^{(2)}$의 error가 $\frac{1}{15}(S_{i}^{(2)}-S_{i}^{(1)})$라면 이 error가 특정 값보다 작을 때는 그 만큼 에러가 작은 것이므로 멈추고, 특정 값보다 클 때는 $\frac{h_{i}}{4}$에서 $\frac{h_{i}}{8}$로 더 쪼갤 수 있습니다.

 

왜 $I_{i}-S_{i}^{(2)}$를 구해서 error를 알아내는 것이 아닐까요? $I_{i}$는 exact한 적분값이기 때문에 실제로는 우리가 모르는 값이기 때문입니다.

 

 

기준 :

$\left\{\begin{matrix}
\left | \frac{1}{15}(S_{i}^{(2)}-S_{i}^{(1)}) \right |\leq \frac{h_{i}}{b-a}\epsilon & \rightarrow\text{Stop}\\ 
\left | \frac{1}{15}(S_{i}^{(2)}-S_{i}^{(1)}) \right |> \frac{h_{i}}{b-a}\epsilon & \rightarrow \text{Subdivide}
\end{matrix}\right.$

적분 전체에 대한 tolerance가 $\epsilon$일 때 각 integral에 대한 tolerance는 $\frac{h_{i}}{b-a}\epsilon$입니다.