콘텐츠로 이동

Simpson 공식을 이용한 적분

다음 정적분의 값을 계산하려고 합니다.

$$\int_a ^ b f (x) dx$$

여기서 설명하는 방법은 1743년 Thomas Simpson의 학위 논문 중 하나에 발표되었습니다.

Simpson 공식

$n$을 어떤 자연수라고 합시다. 적분 구간 $[a, b]$를 $2n$개의 같은 부분으로 나눕니다:

$$x_i = a + i h, ~~ i = 0 \ldots 2n,$$
$$h = \frac {b-a} {2n}.$$

이제 각 구간 $[x_ {2i-2}, x_ {2i}]$, $i = 1 \ldots n$에서 적분을 따로 계산한 다음, 그 값을 모두 더합니다.

그러면 다음 구간 $[x_ {2i-2}, x_ {2i}], i = 1 \ldots n$을 살펴본다고 합시다. 이 구간에서 함수 $f(x)$를 그래프 위의 세 점 $(x_{2i-2}, f(x_{2i-2}))$, $(x_{2i-1}, f(x_{2i-1}))$, $(x_{2i}, f(x_{2i}))$를 지나는 포물선 $P(x)$로 바꿉니다. 이러한 포물선은 항상 존재하며 유일하고, 해석적으로 구할 수 있습니다. 예를 들어 Lagrange 다항식 보간법을 사용하여 이 포물선을 구성할 수 있습니다. 이제 남은 일은 이 다항식을 적분하는 것뿐입니다. 일반적인 함수 $f$에 대해 이를 수행하면 놀라울 만큼 간단한 식을 얻습니다:

$$\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1}) + f(x_{2i})\right)\frac {h} {3} $$

이 값들을 모든 구간에 대해 더하면 최종적인 Simpson 공식을 얻습니다:

$$\int_a ^ b f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2n-1}) + f(x_{2n}) \right)\frac {h} {3} $$

오차

구간 $[a,b]$에 걸친 Simpson 패널 하나의 오차는 다음과 같습니다.

$$ -\tfrac{1}{90} \left(\tfrac{b-a}{2}\right)^5 f^{(4)}(\xi)$$

여기서 $\xi$는 $a$와 $b$ 사이의 어떤 수입니다.

따라서 패널 하나의 오차는 패널 폭에 대해 5차입니다. 위의 복합 공식에서는 구간을 폭이 $2h$인 패널들로 나눕니다. 4계 도함수가 유계라고 가정하면, 패널 오차를 모두 더했을 때 전역 오차는 $O((b-a)h^4)$ 차수가 됩니다. Simpson 공식은 피적분 함수의 값을 계산하는 점들이 각 패널 안에서 대칭적으로 분포하기 때문에, 단순한 보간 오차 추정보다 차수가 하나 더 높습니다.

구현

여기서 $f(x)$는 사용자가 정의한 어떤 함수입니다. 단계 수 $N$은 짝수여야 합니다.

const int N = 1000 * 1000; // number of steps (already multiplied by 2)

double simpson_integration(double a, double b){
    double h = (b - a) / N;
    double s = f(a) + f(b); // a = x_0 and b = x_2n
    for (int i = 1; i <= N - 1; ++i) { // Refer to final Simpson's formula
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}

연습 문제