수치 적분

구적법의 아이디어: 함숫값의 가중합

모든 구적법은 같은 꼴입니다. 노드 $x_i$에서 함숫값을 계산하고 가중치 $w_i$를 곱해 더합니다. 노드와 가중치를 고르는 가장 자연스러운 방법은 노드를 지나는 보간 다항식 $p$(5장)를 만들고 그것을 정확히 적분하는 것입니다. 라그랑주 기저 $\ell_i$로 $p = \sum f(x_i)\ell_i$이므로 $\int p = \sum f(x_i)\int\ell_i$, 곧 가중치는 다음과 같습니다.

정확도 차수는 공식을 검사하는 가장 쉬운 방법이기도 합니다. $1, x, x^2, \dots$를 차례로 넣어 보아 맞지 않는 첫 지수를 찾으면 됩니다. 가중치의 합이 $b - a$여야 한다($f = 1$)는 것이 첫 번째 확인입니다.

뉴턴–코츠 공식: 사다리꼴, 심프슨, 3/8

노드를 $[a, b]$에 등간격으로 놓고 양 끝을 포함한 공식이 닫힌 뉴턴–코츠 공식(Newton–Cotes formula)입니다. 점 두 개(직선)는 사다리꼴, 세 개(포물선)는 심프슨, 네 개(3차식)는 심프슨 3/8 공식입니다. 한가운데 점 하나만 쓰면(열린 공식) 중점 공식입니다.

심프슨이 점 세 개로 3차식까지 맞추는 것은 짝수 개의 구간을 대칭으로 쓰는 공식의 공통된 덤입니다. 그래서 점 네 개를 쓰는 3/8 공식은 차수를 올려 주지 못하고 상수만 조금 바꿉니다. 실무에서 심프슨이 3/8보다 흔한 이유입니다.

합성 공식과 수렴 차수

소구간마다 사다리꼴을 쓰고 더하면, 안쪽 점은 이웃한 두 소구간에 한 번씩 들어가므로 가중치가 $h$, 양 끝은 $h/2$입니다. 심프슨은 소구간 두 개씩 묶어 가중치 $\frac h3(1, 4, 2, 4, \dots, 2, 4, 1)$이 됩니다($n$ 짝수).

위 표는 $\int_0^1 e^x\,dx$의 합성 공식 오차(근삿값 − 참값)입니다. $n$을 두 배로 할 때 비율이 4($=2^2$)와 16($=2^4$)으로 다가갑니다. 이 비율의 $\log_2$가 실험으로 잰 수렴 차수입니다. 심프슨 $n = 32$(함숫값 33개)에서 오차 $10^{-8}$인데, 사다리꼴로 같은 오차를 얻으려면 $n$이 약 4000이어야 합니다. 같은 비용에서는 차수가 높은 공식이 압도적입니다. 단, $f$가 매끄러울 때만입니다.

오일러–매클로린 공식과 롬베르크 적분

이 전개는 $\int_0^1 4/(1 + x^2)dx = \pi$에서 확인할 수 있습니다. $f'(1) - f'(0) = -2$라 첫 항이 $-h^2/6$이고, $n = 4$이면 $-0.0104167$, 실제 오차 $T_{1/4} - \pi = -0.0104162$입니다. 계수가 미리 정해져 있으니 리처드슨 외삽을 반복할 수 있습니다. $R_{i,0} = T_{h/2^i}$에서 출발해

로 표를 채우는 것이 롬베르크 적분(Romberg integration)입니다. 둘째 열 $R_{i,1} = (4T_{h/2} - T_h)/3$은 정확히 합성 심프슨 공식이고, 셋째 열은 불(Boole) 공식입니다. $h$를 반으로 줄일 때 새로 계산할 함숫값은 새 중점들뿐이라 $T_{h/2} = \frac12 T_h + \frac h2\sum f(\text{새 중점})$으로 앞의 계산을 그대로 재사용합니다.

function romberg(f, a, b, kmax, tol)
  h = b - a
  R[0][0] = h * (f(a) + f(b)) / 2
  for i = 1 .. kmax
    h = h / 2
    s = 0
    for j = 1 .. 2^(i-1):  s = s + f(a + (2j - 1) * h)   // 새 중점만
    R[i][0] = R[i-1][0] / 2 + h * s
    for j = 1 .. i
      R[i][j] = (4^j * R[i][j-1] - R[i-1][j-1]) / (4^j - 1)
    if |R[i][i] - R[i-1][i-1]| < tol then return R[i][i]
  return R[kmax][kmax]

가우스 구적법: 노드의 위치까지 최적으로

$[-1, 1]$에서 2점 공식 $w_1f(x_1) + w_2f(x_2)$가 $1, x, x^2, x^3$을 정확히 적분하라는 조건은 $w_1 + w_2 = 2$, $w_1x_1 + w_2x_2 = 0$, $w_1x_1^2 + w_2x_2^2 = \frac23$, $w_1x_1^3 + w_2x_2^3 = 0$입니다. 대칭 $x_2 = -x_1$, $w_1 = w_2 = 1$을 넣으면 $2x_1^2 = \frac23$, 곧 $x_{1,2} = \mp1/\sqrt3 = \mp0.5773503$입니다. 이 점들은 2차 르장드르 다항식 $P_2(x) = \frac12(3x^2 - 1)$의 근입니다. 우연이 아닙니다.

구간 $[a, b]$에서는 $x = \frac{a+b}{2} + \frac{b-a}{2}t$로 옮겨 $\int_a^b f\,dx = \frac{b-a}{2}\sum w_if\bigl(\frac{a+b}{2} + \frac{b-a}{2}t_i\bigr)$로 씁니다. 노드는 $P_n$의 3항 점화식 $kP_k = (2k - 1)xP_{k-1} - (k - 1)P_{k-2}$과 뉴턴 방법(3장)으로 구하거나, 점화식 계수로 만든 대칭 삼중대각 행렬의 고유값으로 구합니다(골럽–웰시 알고리즘, 14장). 노드가 양 끝으로 몰리는 모양은 체비쇼프 노드(5장)와 닮았습니다.

적응 구적법: 어려운 곳만 촘촘히

구간 $[a, b]$에서 심프슨을 한 번($S_1$), 반으로 나눠 두 번($S_2 = S_{\text{왼}} + S_{\text{오}}$) 계산합니다. 심프슨 오차는 $h^5$에 비례하므로 $I - S_2 \approx (I - S_1)/16$. 두 식을 빼면

참값 없이 오차를 추정했습니다. 이 추정이 허용 오차보다 작으면 받아들이고(보정값 $S_2 + (S_2 - S_1)/15$을 쓰면 리처드슨 외삽으로 한 차수를 더 얻습니다), 크면 두 반쪽에 대해 허용 오차를 반씩 나눠 주고 재귀적으로 같은 일을 합니다.

function adapt(f, a, b, fa, fm, fb, S1, tol, depth)
  m = (a + b) / 2;  lm = (a + m) / 2;  rm = (m + b) / 2
  flm = f(lm);  frm = f(rm)                       // 새 함숫값은 두 개뿐
  L = (m - a) / 6 * (fa + 4*flm + fm)
  R = (b - m) / 6 * (fm + 4*frm + fb)
  S2 = L + R
  if depth >= maxdepth or |S2 - S1| <= 15 * tol
    return S2 + (S2 - S1) / 15                   // 받아들이고 외삽
  return adapt(f, a, m, fa, flm, fm, L, tol/2, depth+1)
       + adapt(f, m, b, fm, frm, fb, R, tol/2, depth+1)

$\int_0^1\sqrt x\,dx$에 허용 오차 $10^{-6}$을 주면 함숫값 105개, 소구간 26개로 오차 $2.0\times10^{-8}$을 얻습니다. 가장 작은 소구간은 폭 $4.8\times10^{-7}$으로 $x = 0$에 붙어 있고, 가장 큰 것은 $1/4$입니다. 같은 105개를 균등 분할 심프슨에 쓰면 오차가 $7.7\times10^{-5}$, 4000배 나쁩니다. 반면 어디나 똑같이 매끄러운 함수라면 적응 방법은 오차 추정에 쓰는 계산만큼 손해이고, 추정이 보수적이라 허용 오차보다 훨씬 정확한(그만큼 비싼) 답을 내는 경우가 많습니다.

특이점과 무한 구간: 변수 변환

먼저 매끄럽지 않은 함수에서 무슨 일이 생기는지 봅시다. $\int_0^1\sqrt x\,dx$를 합성 공식으로 계산하면 $n$을 두 배로 할 때 오차가 사다리꼴과 심프슨 모두 약 $2.83 = 2^{1.5}$배씩만 줄어듭니다. 수렴 차수가 $h^{1.5}$로, 공식의 차수와 상관없이 함수의 매끄러움($\sqrt x$의 도함수가 0에서 $x^{-1/2}$로 발산)이 결정합니다. $1/\sqrt x$를 중점 공식(끝점을 쓰지 않는 열린 공식)으로 계산하면 오차가 $\sqrt2$배씩, $h^{0.5}$로 줄어 $n = 128$에서도 $5\times10^{-2}$입니다.

변환 뒤의 두 적분은 다항식이라 심프슨이 정확합니다. 무한 구간 $\int_0^\infty e^{-x^2}dx = \sqrt\pi/2$에 $x = t/(1 - t)$를 쓰면 새 피적분 함수는 $t \to 1$에서 0으로 매끄럽게 사라지고, 합성 심프슨 오차가 $n = 16$에서 $1.4\times10^{-4}$, $n = 64$에서 $1.5\times10^{-11}$입니다. 반면 $\int_0^\infty\frac{dx}{1 + x^2}$을 $[0, 10]$에서 잘라 버리면 꼬리 $\int_{10}^\infty \approx 0.1$을 통째로 잃습니다. 같은 변환을 쓰면 피적분 함수가 $1/\bigl((1 - t)^2 + t^2\bigr)$, $[0, 1]$의 매끄러운 함수가 됩니다.

다차원 적분, 차원의 저주, 몬테카를로

곱 공식의 오차는 축마다 $O(h^p)$이고 총점수가 $N = n^d = h^{-d}$이므로 $O(N^{-p/d})$입니다. 심프슨($p = 4$)도 $d = 8$이면 $N^{-1/2}$, $d = 20$이면 $N^{-1/5}$로 느려집니다. 차원이 높아질수록 같은 정확도에 필요한 점이 지수적으로 늘어나는 이 현상을 차원의 저주(curse of dimensionality)라 합니다. 또 고차원에서는 직관이 틀립니다. 단위 정육면체 $[-1, 1]^d$ 안에서 단위 공이 차지하는 비율은 $d = 2$에서 0.785, $d = 5$에서 0.164, $d = 10$에서 0.0025, $d = 20$에서 $2.5\times10^{-8}$입니다. 부피의 거의 전부가 "모서리"에 있습니다.

몬테카를로 적분(Monte Carlo integration)의 오차는 확률적이고 느립니다. 한 자리 더 맞추려면 표본이 100배 필요합니다. 그러나 $d \gt 8$ 정도부터는 어떤 곱 공식보다 빠르고, 피적분 함수의 매끄러움도 요구하지 않습니다. 정사각형에 점을 뿌려 사분원 안의 비율로 $\pi$를 구하면 $f$는 0 또는 1이라 $\sigma_f = \sqrt{p(1 - p)}$, $p = \pi/4$이므로 $\pi$ 추정의 표준오차는 $4\sqrt{p(1 - p)}/\sqrt N = 1.64/\sqrt N$입니다. $N = 10^6$이면 약 $0.0016$, 소수 셋째 자리 근처까지입니다.

이 도구가 쓰이는 곳

확률 계산에서 정규분포의 누적 확률은 오차 함수로, 오차 함수는 $\operatorname{erf}(x) = \frac{2}{\sqrt\pi}\int_0^xe^{-t^2}dt$로 정의됩니다. 라이브러리는 이를 구간별 유리 근사로 계산하지만, 그 근사의 검증은 고정밀 구적으로 합니다. 물리량 적분에서는 일 $W = \int F\,dx$, 질량 중심, 관성 모멘트, 곡선의 길이가 모두 이 장의 공식으로 계산됩니다.