가우스 구적법 (Gauss Quadrature

수치 적분법

 

수치 적분(Numerical integration)은 컴퓨터 알고리듬을 이용하여 해석적으로 적분하기 어려운 피적분 함수를 적분하는 방법입니다. 예를 들어,

$$\int_{a}^{b} x^2 + 2x +3 dx$$

와 같이 다항함수를 구간 $[a,b]$에서 적분하기 위해서는 미적분학 기본정리를 이용하면 됩니다. 피적분 함수인 $f(x) = x^2 + 2x +3$의 부정적분 $F(x) = \frac{1}{3}x^3 + x^2 + 3x + C$를 구한다음 $F(b) - F(a)$를 계산하면 이 값은 위 적분의 값이 됩니다. 

 

하지만, 

$$\int_{a}^{b} \sin(x^2) \log\Big|\frac{1}{x^2-\cos x}\Big| dx$$

와 같이 피적분함수가 매우 복잡한 경우에는 부정적분을 구할 수 없고, 해석적인 방법으로는 적분이 불가능 합니다. 이 경우에 수치 적분법을 사용하여 위 적분의 값을 계산하게 됩니다. 

 

수치 적분의 가장 간단한 방법은 고등학교 시절 처음으로 적분의 의미를 배울 때 나왔던 구분구적법입니다. 

위 그림은 함수 $f(x) = x^2$을 $[0,1]$ 구간에서 구분구적법으로 구하는 과정의 개념적인 설명입니다. 적분 하려는 구간을 매우 많은 사각형으로 나누고, 이 사각형의 넓이를 모두 더하면 원래 구하고자 했던 영역의 넓이(적분값)이 됩니다. 이 과정에서 무한대의 개념이 들어오게 되는데, 실제 컴퓨터 계산에서는 무한대를 계산하지 못 하기 때문에 무한대 대신에 매우 큰 값, 예를들면 $10^6$ 정도의 값을 사용할 수 있습니다. 수식으로 나타내면,

$$S = \int_{a}^{b} f(x) dx \approx \sum_{i=1}^{N} \frac{b-a}{N} f(x_i), \text{     here, } x_i = \frac{b-a}{N}i + a$$

입니다. 이 과정에서 적분 영역 $[a,b]$를 $N$ 조각으로 균등하게 자르는 것을 가정하였습니다. 

 

위 방법에서는 적분 영역의 조각 $[x_i, x_{i+1}]$ 구간에서 함수값이 $f(\frac{x_i + x_{i+1}}{2})$로 변하지 않는 값이라고 가정하였습니다. 예를들어 위 그림에서와 같이 $[0.55, 0.65]$ 구간에서 함수값은 실제로 변하지만, 구분구적법에서는 함수값이 $f(0.6)$ 값으로 고정되어 있다고 가정한 것 입니다. 따라서 이 구간에서 영역은 직사각형이 되고 넓이를 쉽게 구할 수 있었던 것 입니다. 매우 간단하게 근사를 한 것인데, 좀 더 정확하게 적분값을 구하기 위해서는 해당 구간에서 함수가 다항함수로 변한다고 가정을 하고 적분을 할 수 있습니다. 이 경우에는 수치적분의 식은 좀 더 복잡해지겠지만, 결과는 더 참값에 유사할 것 입니다. trapezoidal방법, Simpson'r rule이 대표적인 예시 입니다. 

 

식을 일반화 시켜 보자

 

위 구분구적법의 식을 일반화 시켜서 표현하면 다음과 같습니다. 

$$S = \int_{a}^{b} f(x) dx \approx \sum_{i=1}^{N} \omega_i f(x_i), x_i \in [a,b], w_i = \text{weight at } x_i$$

 

 

위 구분구적법에서는 $\omega_i = \frac{b-a}{N}, x_i = \frac{b-a}{N}i +a$ 였습니다. 즉 $x_i$는 균등하게 분포되어 있으며, 각 $x_i$에서의 비중(weight)는 역시 고정된 값, $\frac{b-a}{N}$, 으로 균등했습니다. 

 

여기서 한가지 질문을 할 수 있는데, "$x_i$의 분포는 꼭 균일해야 할까요? 또한 $\omega_i$는 고정된 값이어야할까요?" 혹은 좀 더 일반적으로 표현하면 "만일 우리가 영역 $[a,b]$ 내의 임의의 $n$ 개의 지점 $\{ x_1, x_2, ..., x_n\}$에서의 함수값 $\{ f(x_1), f(x_2), ..., f(x_n) \}$ 을 알 수 있다고 할 때, 이 $\{ x_1, x_2, ..., x_n\}$와 이때의 비중 $\{ \omega_1, \omega_2, ..., \omega_n\}$" 를 어떻게 설정한다면, 실제 적분값과 가장 유사한 적분값을 얻을 수 있을까요?"

 

위 마지막 질문을 또 다르게 해석하면 다음과 같습니다. 

 

(1) 적분 구간 $[a,b]$는 일반성을 잃지 않고 $[-1, 1]$이라고 정할 수 있다. ($a=-1, b=1$이 아닌 경우에는 간단한 치환적분을 통해서 적분 구간을 $[-1,1]$로 바꿀 수 있습니다)

(2) 피적분 함수 $f(x)$에 대한 일반적인 정보는 없다.

(3) 단, 구간 $[-1,1]$ 에 속한 임의의 $n$개의 점 $\{ x_1, x_2, ..., x_n\}$에서 함수값 $\{ f(x_1), f(x_2), ..., f(x_n) \}$은 그 값은 알 수 있다. 

이 경우, $\int_{-1}^1 f(x) dx$를 가장 정확하게 얻기 위해서는 $\{ x_1, x_2, ..., x_n\}$와 $\{ \omega_1, \omega_2, ..., \omega_n\}$를 어떻게 디자인 해야하는가?

 

입니다. (다시 한 번 더 설명하면) 위에서 설명한 구분구적법에서는 $n$개의 점은 $x_i = \frac{b-a}{N}i +a$와 같이 균등한 분포였으며, 비중은 $\omega_i = \frac{b-a}{N}$으로 일정한 값이었습니다. 똑똑한 방법으로 $x_i, \omega_i$를 잘 설계하면, 좀 더 정확한 적분값을 얻을 수 있을까요? 

 

피적분 함수 $f(x)$가 2차 함수라는 사실을 알고 있다면?

 

위에서, 피적분 함수 $f(x)$에 대한 일반적인 정보는 없다고 가정하였습니다. 이는 너무 어려운 가정이니, 일단 조금 쉬운 문제를 풀기 위해서, 피적분 함수 $f(x)$가 2차 함수라고 가정을 하겠습니다. 즉, $f(x) = ax^2 + bx +c$ 의 형태 입니다. 물론, $a,b,c$는 모두 실수이며, $a$는 $0$이 아닙니다. 실제로 적분을 하면,

$$S = \int_{-1}^1 ax^2 + bx +c = \Big[\frac{1}{3}ax^2 + \frac{1}{2}x^2 + cx \Big]^1_{-1} = \frac{2}{3}a + 2c$$

입니다. 최종 결과를 $ \frac{2}{3}a + 2c = \Big(a(\frac{1}{\sqrt{3}})^2 + c\Big) + \Big(a(\frac{1}{\sqrt{3}})^2 + c\Big)$ 와 같이 나눠 쓸 수 있는데, 이는

$$S = f\Big(-\frac{1}{\sqrt{3}}\Big) + f\Big(\frac{1}{\sqrt{3}}\Big)$$

와 같다는 것을 알 수 있습니다! 매우 놀라운 결과인데, 임의의 2차 함수의 $[-1,1]$ 구간에서 적분값은 $x = \pm \frac{1}{\sqrt{3}}$에서의 함수값이 합과 같습니다. 단지 적분 구간에서 두 점의 함수값 만으로 적분값을 정확히 계산할 수 있는 것 입니다. 일반적인 피적분 함수에 대해서는 아직 어떠한 $x_i, \omega_i$ 디자인이 가장 효율적인 디자인이 되는지는 모르지만, 피적분 함수를 2차 함수로 국한 시킨다면 가장 효율적인 디자인은 $n=2$ 이며 $x = \{-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}} \}$ 이고 이때의 비중은 $\omega_i = 1$ 이 됩니다. 

 

르장드르 다항식

 

수치 적분과는 직접적으로 관련이 없어 보이는 이야기를 잠시 하도록 하겠습니다. 르장드르 다항식에 대한 이야기 입니다. 르장드르 다항식은 수학 곳곳에서 접할 수 있는데, 물리학이 전공인 필자는 전자기학 시간에 구면좌표계에서 라플라스 방정식을 푸는 과정에서 처음 접하였습니다. (자세한 설명을 생략하고....) 르장드르 다항식은 미분방정식

$$(1-x^2) P_n^{''}(x) - 2xP_n^{'}(x) + n(n+1) P_n(x) = 0$$

의 해 $P_n(x)$ 입니다. 여기서 $^{'}$는 미분연산자이고 $n$은 0과 양의 정수 입니다. 함수의 정의역은 $[-1,1]$입니다. 르장드르 다항식의 특징을 언급하면,

 

(1) 구체적인 예를 들면 $P_0 = 1, P_1(x) = x, P_2(x) = \frac{1}{2}(3x^2-1)$ 으로, $P_n(x)$는 n차 다항식이다.

(2)

$$\int_{-1}^{1} P_n(x) P_m(x) = \frac{2}{2n+1} \delta_{nm}$$

이고, 이 성질을 이용하면, 일반적인 m차 다항식은 $n \le m$인 르장드르 다항식 $P_n(x)$의 합으로 쓸 수 있다. 즉, 임의의 $n$차 다항식 $f(x)$는

$$f(x) = \sum_{n=0}^{m} a_n P_n(x)$$ 

가 되는 $a_n$가 존재한다.

(3) (2)의 성질을 이용하면, $m \lt n$ 일때,

$$\int_{-1}^{1} x^m P_n(x) dx = 0$$

이다.

 

입니다. 이 밖에도 매우 재미난 성질이 여럿 있지만, 이번 포스팅과는 크게 관련이 없는 사항이라서 생략하도록 하겠습니다. 르장드르 다항식에 대해서 익숙하지 않으신 분이라면, 이후 글을 읽는 것 보다는 르장드르 다항식에 대해서 좀 더 공부를 하는 편이 좋습니다. 

 

라그랑쥬 다항식

 

한 번 더 수치 적분과는 직접적으로 관련이 없어 보이는 이야기를 하겠습니다. 르장드르가 나왔으니 라그랑주가 나와야겠죠. 르장드르만 자신의 다항식을 갖고 있는 것은 아닙니다. 라그랑주 역시 자신의 이름을 딴 다항식을 갖고 있습니다. 좌표 평면 위의 점들 $\{(x_0, y_0), (x_1, y_1), (x_2, y_2), ..., (x_n, y_n) \}$이 주어졌을 때, 이 점을 지나는 $n$ 차 다항식 $y(x) = L(x)$ 는 단지 하나 존재합니다. $n$차 방정식의 각 항의 계수가 $n+1$개 라는것을 생각하면 자연스러운 결론입니다. 실제로 $\{(x_0, y_0), (x_1, y_1), (x_2, y_2), ..., (x_n, y_n) \}$이 주어졌을 때, 이를 만족하는 $L(x)$는 어떻게 구할 수 있을까요? 이 $L(x)$가 바로 라그랑쥬 다항식으로

$$L(x) = \sum_{j=0}^n y_j l_j(x)$$

$$l_j(x) = \prod_{0 \le m \le n, m \ne j} \frac{x-x_m}{x_j-x_m}$$

와 같습니다. 

 

다항식의 몫과 나머지 표현

 

정수 $x$를 $p$로 나눈 몫을 $q$, 나머지를 $r$이라고 할 때, $x = p \cdot q +r$로 표현할 수 있습니다. 다항식도 비슷한 방법으로 표현할 수 있는데, $m$차 방정식 $f(x)$를 $n$차 방정식 $p(x)$로 나눈 몫과 나머지를 각각 $q(x), r(x)$라 한다면, 

$$f(x) = p(x) \cdot q(x) + r(x)$$

와 같이 표현할 수 있습니다. 몫 $q(x)$는 $m-n$차 방정식이며, 나머지 $r(x)$의 차수는 $n$보다 작습니다. 예를들어서, $f(x)$가 $2n+1$차 방정식이고, $p(x)$가 $n+1$차 방정식이라면, 몫 $q(x)$와 나머지 $r(x)$는 모두 $n$차 이하의 방정식이 됩니다. $f(x)$를 구간 $[-1,1]$에서 적분하면, 그 값은

$$\int_{-1}^{1}f(x) dx = \int_{-1}^{1} p(x) \cdot q(x) dx + \int_{-1}^{1} r(x)dx$$

가 됩니다. 

 

앞에서 예시로 들었던 것과 같이 $2n+1$차 방정식 $f(x)$에 대해서, $p(x)$가 $n+1$차 르장드르 다항식이라면,

$$\int_{-1}^{1}f(x) dx = \int_{-1}^{1} P_{n+1}(x) \cdot q(x) dx + \int_{-1}^{1} r(x)dx$$

가 됩니다. 여기서 $n$차 다항식인 $q(x)$는 위 르장드르 다항식의 성질 (2)로 부터 

$$q(x) = \sum_{k=0}^{n} a_k P_k(x)$$

와 같이 쓸 수 있고, 따라서 $f(x)$의 적분은

$$\int_{-1}^{1} P_{n+1}(x) \cdot q(x) dx + \int_{-1}^{1} r(x)dx = \int_{-1}^{1} P_{n+1}(x) \cdot \Big( \sum_{k=0}^{n} a_k P_k(x)\Big)dx + \int_{-1}^{1} r(x)dx$$

$$ = \sum_{k=0}^{n} a_k\int_{-1}^{1} P_{n+1}(x) P_k(x) dx +  \int_{-1}^{1} r(x)dx $$

가 됩니다. 르장드르 다항식의 성질 (3)으로 부터 마지막 식의 첫번째 항의 적분 값 $\int_{-1}^{1} P_{n+1}(x) P_k(x) dx$은 $0$이 되는 것을 알 수 있고,

$$\int_{-1}^{1}f(x) dx = \int_{-1}^{1} r(x)dx$$

우리가 구하려는 적분 값은 나머지 $r(x)$의 적분값과 일치하게 됩니다! $2n+1$차 방정식을 적분한 값과, 이 방정식을 $n+1$차 르장드르 다항식 $P_{n+1}(x)$로 나눈 나머지 값을 적분한 값이 같다는 것이 (수식을 통해서 보면 쉽게 이해가 되지만) 매우 신기합니다. 

 

가우스 구적법 정리

 

$n$차 다항식인 $r(x)$의 서로 다른 $x$에서의 함수값 $n+1$개가 있다면, $r(x)$를 라그랑쥬 방정식을 통해 표현할 수 있고, $r(x) = \sum_{i=0}^{n}r(x_i) l_i(x)$ 가 됩니다. 라그랑쥬 방정식을 얻는데 이용할 $n+1$ 개의 $x$값 $x_i$를 $P_{n+1}(x)$ 의 $n+1$개의 해로 선정한다면, 즉, $P_{n+1}(x_i) = 0$ 이라고 한다면,

$$f(x_i) = P_{n+1}(x_i) q(x_i) + r(x_i)$$

에서 $P_{n+1}(x_i) =0$ 이므로 $f(x_i) = r(x_i)$가 되고, 

$$\int_{-1}^{1}f(x) dx = \int_{-1}^{1} r(x)dx = \int_{-1}^{1} \sum_{i=0}^{n}r(x_i) l_i(x)dx$$

$$=\int_{-1}^{1} \sum_{i=0}^{n}f(x_i) l_i(x)dx = \sum_{i=0}^{n}f(x_i) \Big(\int_{-1}^{1}  l_i(x)dx\Big)$$

가 됩니다. 

 

여기서 $l_i(x)$는 $P_{n+1}(x)$ 의 $n+1$개의 해로 부터 정해지는 방정식으로 적분값 $\int_{-1}^{1}  l_i(x)dx$ 는 $n+1$이 정해진다면 정해지는 상수값($f(x)$와 상관 없이 주어지는 값) 입니다. 이 값을 $A^{n+1}_i$ 이라고 놓는다면, 최종적인 적분값을 정리하면

$$\int_{-1}^{1}f(x) dx = \sum_{i=0}^{n}f(x_i) A^{n+1}_i$$

가 됩니다! 임의의 $2n+1$차 다항식에 대해서, $n+1$차 르장드르 다항식 $P_{n+1}(x)$의 $n+1$개의 해 $\{x_i | i=0,1, 2, ..., n\}$에서의 합수값 $\{ f(x_i) | i=0, 1, 2, ..., n\}$을 알고 있다면 이 값만을 이용하여 적분값은 $\sum_{i=0}^{n}f(x_i) A^{n+1}_i$와 같다는 것을 증명하였습니다. 위에서 2차 함수의 적분값은 단지 $x = \pm \frac{1}{\sqrt{3}}$에서의 함수값만 알면 구할 수 있다는 것의 일반화된 증명입니다. 위에서 $A^{n+1}_i$가 우리가 찾고자 한 $\omega_i$ 였고, $n+1$차 르장드르 다항식의 해가 우리가 찾고자 한 $x_i$였습니다. 

 

이 방법을 이용하여 수치 적분을 하는 것을 가우스 구적법 이라고 합니다. 르장드르 다항식, 라그랑쥬 다항식 정도만 알고 있다면 충분히 증명할 수 있는 간단한 정리인데, 이와 같이 적분을 하는 방법을 처음 알았을 때 매우 충격적이었습니다. 제가 생각하는 적분은 적분 구간의 모든 점에 대해서 피적분함수의 함수값을 알아야 하는 것 이었지만, 매우 적은 갯수의 $x$값에서의 피적분함수의 함수값만 알고 있다면 매우 정확하게 적분값을 구할 수 있으니까요. 한편, $n$차 다항식은 $n+1$개의 실수만 있으면 그 정보가 완전히 정해진다는 것을 생각하여 조금만(?) 확장을 시키면 얻을 수 있는 당연한(?) 결과 인것 같기도 합니다. 

 

예시

 

$$\int_{1}^{1.5} x^2 \log x dx = ?$$

가우스 구적법 $n=1$을 써서 계산해 보도록 하겠습니다. 우선 $P_2(x)$의 해를 알아야 하는데, $P_2(x) = \frac{1}{2} (3x^2-1)$이므로 $x_i = \pm \frac{1}{\sqrt{3}}$ 이 됩니다. 또한 $A^2_i$를 구해야하는데, 이 값은 (직접 계산해 보는 것은 숙제로 남깁니다) 1이 된다는 것을 알 수 있습니다. 

 

위 적분의 적분 영역은 $[1, 1.5]$로 지금까지 우리가 해 오던 적분 영역 $[-1,1]$와는 다릅니다. 따라서 변수 변환을 통해서 적분 영역을 원하는 영역으로 변경해 주어야 하는데, $t = (0.5\cdot t + 2.5)/2$ 의 변수 변환을 통해서, $f(x) = x^2 \log x$ 일때,

$$\int_{1}^{1.5} x^2 \log x dx = \frac{1}{4} \int_{-1}^{1} f\Big(\frac{x+5}{4}\Big) dx= \frac{1}{4} \int_{-1}^1 \Big(\frac{x+5}{4} \Big)^2 \log \Big(\frac{x+5}{4} \Big) dx$$

가 됩니다. 위 식의 가장 마지막항은 우리가 원하는 형태의 적분 형태이고, 따라서 적분값은

$$\frac{1}{4} \Big[ \Big( \frac{\frac{-1}{\sqrt{3}}+5}{4} \Big)^2 \log \Big( \frac{\frac{-1}{\sqrt{3}}+5}{4} \Big) +\Big( \frac{\frac{1}{\sqrt{3}}+5}{4} \Big)^2 \log \Big( \frac{\frac{1}{\sqrt{3}}+5}{4} \Big)  \Big] = 0.192269$$

로 근사할 수 있게 됩니다. 실제 더 정확한 방식으로 얻은 적분값은 0.192259로 참값과 우리가 구한 근사값과의 오차는 0.00001 밖에 되지 않습니다. 단지 두 점에서의 피적분 함수만을 가지고 이렇게 정확한 적분을 할 수 있다는 것이 다시금 놀라워집니다!