수치적분_Numerical Integral

$\sin(x^2), 1/\ln x, \sqrt{1+x^4}$과 같은 함수는 역도함수(Antiderivative)를 간단하게 나타낼 수 없다. 적분해야 하는 함수 $f$가 쓸만한 역도함수(Antiderivative)를 찾을 수 없을 때, 주어진 부분구간(suinterval)에서 함수 $f$와 비슷한 다항함수(polynomial)로 대체한 다음 이 다항함수를 적분한 값을 더하여 함수 $f$의 적분값을 근사한다. 이런 과정을 수치적분으로 부른다. 여기서는 사다리꼴 규칙(Trapezoidal rule)과 심슨 규칙(Simpson's rule)을 공부하자. 구간 $[a,b]$에서 함수 $f$는 양이고 연속이라고 가정하자.

먼저 아래와 같은 영역의 넓이를 구하는 것을 생각해 보자.

정적분의 정의와 연관지어 생각한다면 아래와 같은 직사각형의 넓이를 더한 것으로 근삿값을 구할 수 있다. 분할(partition) $P$에서 $||P||\rightarrow 0$일 때 극한값으로 넓이를 정의한다.

$$A=2\cdot 7+2\cdot 6 +2\cdot 5+2\cdot 6 +2\cdot 7 +2\cdot 5=2(7+6+5+6+7+5)=72$$

이것보다 더 빠르게 근사시키는 방법으로 아래와 같은 사다리꼴을 쓸 수 있다.

$$A=\frac{1}{2}\cdot 2\{(5+7)+(7+6)+(6+5)+(5+6)+(6+7)+(7+5)\}=72$$

분할점이 많아지면 점점 가까운 근삿값을 구할 수 있지만 속도가 문제다. 두 가지 방법으로 구한 값이 우연히 같아졌지만 직관적으로 사다리꼴을 쓰는 것이 빠르다.

이제 일반화해 보자.

사다리꼴로 구한 근삿값(Trapezoidal Approximations)

구간 $[a,b]$를 등분한 분할을 생각하자. $a=x_0 <x_1 <x_2 <\cdots <x_{n-1}<x_n =b$

부분구간의 길이(step size 또는 mesh size)는 아래와 같다. 

$$\Delta x=\frac{b-a}{n}$$

$y_i=f(x_i)$일 때, $i$번 째 구간에 있는 사다리꼴 넓이는 아래와 같다.

$$\Delta x \big(\frac{y_{i-1}+y_i}{2}\big)=\frac{\Delta x}{2}(y_{i-1}+y_i)$$

$\Delta x$를 높이(altitude)로 삼는다면 수직으로 놓인 두 변의 평균과 곱한 값이 넓이다. 곡선 아래 넓이는 아래와 같이 모든 사다리꼴 넓이를 더하여 근삿값을 구한다.

$$T=\frac{1}{2}(y_0+y_1)\Delta x +\frac{1}{2}(y_1+y_2)\Delta x+\cdots+\frac{1}{2}(y_{n-2}+y_{n-1})\Delta x+\frac{1}{2}(y_{n-1}+y_n)\Delta x$$

$$=\frac{\Delta x}{2}(y_0 +2y_1 +2y_2 +2y_3+\cdots+2y_{n-1}+y_n)$$

The Trapezoidal Rule

$$\int_{a}^{b} f(x)dx \approx \frac{\Delta x}{2}(y_0 +2y_1 +2y_2 +2y_3+\cdots+2y_{n-1}+y_n)$$

포물선으로 구하는 근삿값(Simpson's Rule : Approximation Using Parabolar)

아래 그림과 같이 세점을 지나는 포물선을 구하여 적분한 값을 더하여 근사값을 구할 수 있다.

먼저 세 점 $(x_{i-1},y_{i-1}),(x_i,y_i),(x_{i+1},y_{i+1})$을 지나는 포물선을 구해 보자. 계산을 간단하게 하기 위하여 $x_0 =-h, x_1 =0 ,x_2=h$이고 $h=\Delta x=(b-a)/n$이라고 하자.

포물선의 방정식을

$$y=Ax^2 +Bx+C$$라고 한다면 아래 면적은

$$A_p =\int_{-h}^{h} (Ax^2 +Bx+C)dx=2\int_{0}^{h} (Ax^2 +C)dx=2\Big[\frac{1}{3}Ax^3 +Cx\Big]_{0}^{h}= \frac{h}{3}(2Ah^2 +6C)$$

세 점 $(-h,y_0),(0,y_1),(h,y_2)$가 포물선 위에 있으므로

$$y_0=Ah^2 -Bh +C,\;\; y_1=C,\;\;y_2 =Ah^2 +Bh+C$$이다.

여기서 $C=y_1, 2Ah^2 =y_0+y_2 -2y_1$임을 쉽게 얻을 수 있다. 따라서 위에서 구한 넓이는

$$A_p =\frac{h}{3}(2Ah^2 +6C)=\frac{h}{3}((y_0+y_2 -2y_1)+6y_1)=\frac{h}{3}(y_0+ 4y_1 +y_2)$$이다.

정리하면 $ (x_0,y_0),(x_1,y_1),(x_2,y_2)$를 지나는 포물선 아래 넓이는

$$\frac{h}{3}(y_0+ 4y_1 +y_2)$$

이다. 마찬가지로 $ (x_2,y_2),(x_3,y_3),(x_4,y_4)$를 지나는 포물선 아래 넓이는

$$\frac{h}{3}(y_2+ 4y_3 +y_4)$$

이다. 그러므로 아래와 같이 근삿값을 구할 수 있다.

$$\int_{a}^{b} f(x)dx \approx \frac{h}{3}(y_0+ 4y_1 +y_2)+\frac{h}{3}(y_2+ 4y_3 +y_4)\cdots +\frac{h}{3}(y_{n-2}+ 4y_{n-1} +y_{n})$$

$$=\frac{h}{3}(y_0+ 4y_1 +2y_2+4y_3 +2y_4+\cdots+2y_{n-2}+4y_{n-1}+y_n)$$

이 방법을 심프슨 규칙이라고 한다. Thomas Simpson(1720~1761)

Simpson's Rule

$$\int_{a}^{b} f(x)dx \approx \frac{\Delta x}{3}(y_0+ 4y_1 +2y_2+4y_3 +2y_4+\cdots+2y_{n-2}+4y_{n-1}+y_n)$$