가우스 구적법 (Gauss Quadrature

수치 적분법

 

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

$${\int }_a^b{x}^2+2x+3dx$$

와 같이 다항함수를 구간 $[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 |\frac{1}{x^2-\cos x}|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)=a{x}^2+bx+c$ 의 형태 입니다. 물론, $a,b,c$는 모두 실수이며, $a$는 $0$이 아닙니다. 실제로 적분을 하면,

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

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

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

와 같다는 것을 알 수 있습니다! 매우 놀라운 결과인데, 임의의 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)-2x{P}_n^{{}^{\prime }}(x)+n(n+1)P_n(x)=0$$

의 해 $P_n(x)$ 입니다. 여기서 ${}^{{}^{\prime }}$는 미분연산자이고 $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<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 (\sum _{k=0}^n{a}_k{P}_k(x))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)({\int }_{-1}^1{l}_i(x)dx)$$

가 됩니다. 

 

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

$${\int }_{-1}^{1}f(x)dx=\sum _{i=0}^{n}f(x_i)A_i^{n+1}$$

가 됩니다! 임의의 $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_i^{n+1}$와 같다는 것을 증명하였습니다. 위에서 2차 함수의 적분값은 단지 $x=\pm \frac{1}{\sqrt{3}}$에서의 함수값만 알면 구할 수 있다는 것의 일반화된 증명입니다. 위에서 $A_i^{n+1}$가 우리가 찾고자 한 ${\omega }_i$ 였고, $n+1$차 르장드르 다항식의 해가 우리가 찾고자 한 $x_i$였습니다. 

 

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

 

예시

 

$${\int }_1^{1.5}{x}^2\log xdx=?$$

가우스 구적법 $n=1$을 써서 계산해 보도록 하겠습니다. 우선 $P_2(x)$의 해를 알아야 하는데, $P_2(x)=\frac{1}{2}(3x^2-1)$이므로 $x_i=\pm \frac{1}{\sqrt{3}}$ 이 됩니다. 또한 $A_i^2$를 구해야하는데, 이 값은 (직접 계산해 보는 것은 숙제로 남깁니다) 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 xdx=\frac{1}{4}{\int }_{-1}^{1}f(\frac{x+5}{4})dx=\frac{1}{4}{\int }_{-1}^1(\frac{x+5}{4}{)}^2\log (\frac{x+5}{4})dx$$

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

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

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