사다리꼴 공식을 이용한 C/C++ 수치 적분

여기서는 사다리꼴 공식 (trapezoidal rule)을 적용하여 주어진 함수의 정적분을 수치적으로 구하는 법에 대해서 알아봅시다. 사다리꼴 공식으로 얻을 수 있는 것은 적분의 근사값입니다만, 원하는 정밀도에 도달할때까지 격자의 갯수를 늘려나가는 방식을 C언어나 C++로 구현하는 것에 대해서도 짚어보겠습니다.

 

반응형

 

적분의 수학적인 정의 및 정적분과 부정적분의 차이에 대해서는 다음 포스팅에 더 자세한 내용이 소개되어 있습니다. 적분의 개념이 생소하게 느껴지는 분들에게 도움이 되리라 생각합니다.

 

수학 상식 : 미분과 적분 이해하기

 

수학적 원리

사다리꼴 공식의 기본 취지를 간단히 말하자면, 적분 대상이 되는 함수의 형태에 따라 여러개의 사다리꼴을 나열한 뒤 그 넓이의 총 합을 통해 정적분의 근사값을 구하는 것입니다. 정적분의 기하학적인 의미가 함수와 가로축 사이의 넓이라는 것과 일맥상통한다고 볼 수 있죠.

 

 

전체 적분구간을 N등분한 뒤, 각 구간의 경계가 되는 두 지점에서의 함수 값을 윗변과 아랫변으로 해서 사다리꼴을 만들 수 있습니다. 일반적인 비선형 함수의 그래프는 곡선의 형태를 띠는 반면에 사다리꼴의 빗변은 직선이기 때문에, 사다리꼴 공식으로 얻을 수 있는것은 정적분의 근사값이라는 점에 유념할 필요가 있습니다. 그럼에도 불구하고 적분 구간을 세밀하게 나눌수록 사다리꼴 공식에 입각한 근사값은 실제 정적분의 값에 근접하게 됩니다. 다시 말해서 나눠진 구간의 갯수 N이 클수록 정밀도가 올라가게 된다는 것입니다.

 

 

복잡한 함수의 정적분을 컴퓨터로 계산하는데 있어서 사다리꼴 공식을 사용할 수 있습니다. 전체 적분구간을 나눈 갯수를 2배로 늘려가면서, 원하는 정밀도에 도달하는지의 여부를 보는 것인데요. 적분구간을 더 세밀하게 나눠서 사다리꼴 공식을 적용해도 근사값이 크게 달라지지 않는다면, 정적분의 실제 값에 근접했다고 볼 수 있는 것입니다.

 

C/C++ 프로그램

사다리꼴 공식을 통해 정적분을 계산하는 C/C++ 코드를 소개합니다.

 

전역 변수와 함수의 프로토타입을 포함하고 있는 헤더 파일입니다.

 

IntegralTrapezoidal.h [다운로드]

 

#ifndef _INTEGRAL_TRAPEZOIDAL_H_
#define _INTEGRAL_TRAPEZOIDAL_H_

#undef __WRAP_CXX_INI
#undef __WRAP_CXX_FIN
#ifdef __cplusplus
#define __WRAP_CXX_INI extern "C" {
#define __WRAP_CXX_FIN }
#else
#define __WRAP_CXX_INI /* empty */
#define __WRAP_CXX_FIN /* empty */
#endif

__WRAP_CXX_INI

extern FILE *log_integral_trapezoidal_;

double get_integral_trapezoidal(double xmin, double xmax,
                                double eps_precision,
                                double (*ptr_func)(double));

__WRAP_CXX_FIN

#endif

 

정적분의 근사값을 리턴하는 함수 get_integral_trapezoidal의 프로토타입이 있는데, 이 함수는 적분구간의 하한 (lower limit) xmin과 상한 (upper limit) xmax 및 정밀도 eps_precision를 매개변수로 받고 있습니다. 적분 대상이 되는 함수 (integrand)의 함수 포인터 ptr_func 역시 매개변수로 들어가고 있으며, 이렇게 하면 함수를 다른 함수의 매개변수로 넘겨주는 것이 가능하죠.

 

추가로 log_integral_trapezoidal_이라는 파일 포인터가 전역 변수로 등장하는데요. 적분 구간을 나눈 갯수를 늘려가는 과정에서, 각 단계별로 정적분의 근사값이 어떻게 달라지는지를 출력하기 위한 목적이 있습니다. 만약에 이 파일 포인터의 값이 NULL인 경우, 아무것도 출력하지 않습니다. 만약에 터미널 콘솔 창에 출력을 하고 싶다면, 표준 출력을 나타내는 stdout 등을 사용하는 것이 가능하죠.

 

get_integral_trapezoidal 함수는 C언어로 정의되지만, extern "C" 코드블록과 __cplusplus 식별자가 사용된 덕분에 C++ 프로그램에서도 이 함수를 사용할 수 있습니다.

 

get_integral_trapezoidal 함수가 정의된 소스 파일입니다.

 

IntegralTrapezoidal.c [다운로드]

 

#include<stdio.h>
#include<math.h>
#include"IntegralTrapezoidal.h"

FILE *log_integral_trapezoidal_ = NULL;

double get_integral_trapezoidal(double xmin, double xmax,
                                double eps_precision,
                                double (*ptr_func)(double)) {

    double ret_now = 0.;
    double ret_prev;
    int converging = 0;

    unsigned long int nbin_x = 1;
    double delta_x =
        (xmax - xmin) / (double)nbin_x;

    int istep = 1;
    while (converging == 0) {
        ret_prev = ret_now;

        unsigned int ix;
        if (istep == 1) {
            ret_now = 0.;
            for (ix = 0; ix < nbin_x; ix++) {
                double x0 = xmin + delta_x * (double)ix;
                double x1 = x0 + delta_x;
                ret_now +=
                    0.5 * delta_x * ((*ptr_func)(x0) +
                                     (*ptr_func)(x1));
            }
        } else {
            ret_now = 0.5 * ret_prev;
            for (ix = 0; ix <= nbin_x; ix++) {
                if (ix % 2 == 0) {
                    continue;
                }

                double x_now = xmin + delta_x * (double)ix;
                ret_now +=
                    delta_x * (*ptr_func)(x_now);
            }
        }

        if (log_integral_trapezoidal_ != NULL) {
            fprintf(log_integral_trapezoidal_,
                "      step %d : integral = %.8f\n",
                istep, ret_now);
    	}
            
        if (fabs(ret_now - ret_prev) <
            0.5 * eps_precision *
            fabs(ret_now + ret_prev)) {
            converging = 1;
        }

        istep += 1;
        nbin_x = 2 * nbin_x;
        delta_x = 0.5 * delta_x;
    }

    return ret_now;
}

 

앞에서 언급한대로 적분구간을 더 세밀하게 나눠가면서, 사다리꼴 공식에 따라 근사값을 계산합니다. 적분구간을 나눈 갯수 nbin_x의 초기값은 1이지만, 다음 단계로 넘어갈때마다 그 값이 2배로 증가합니다. 사다리꼴의 갯수가 2배로 늘어나는거죠. 이전 단계와 현재 단계의 근사값을 ret_prev, ret_now라는 변수들에 각각 저장하고 이들의 값을 비교해서 수렴 여부를 결정합니다. 적분 구간을 더 세밀하게 나눠도 적분 값이 별로 달라지지 않는다면, 수렴한다고 가정하고 제어변수 converging의 값을 1로 바꿔서 반복문을 빠져나오게 됩니다.

 

간단한 예시로서 가우스 함수의 정적분에 대해서 살펴봅시다. 가우스 함수는 정규분포 (normal distribution)의 확률밀도함수로서, 평균값과 표준편차의 값에 따라 그 형태가 특정됩니다. 가우스 함수의 정의와 특징에 대해서는 다음 포스팅에 더 자세한 내용이 소개되어 있습니다.

 

수학 상식 : 정규 분포와 가우스 함수

 

사다리꼴 공식을 적용해서 계산한 정적분의 근사값을 누적분포함수로부터 얻어진 것과 비교해 보겠습니다.

 

main 함수가 정의된 C++ 소스 파일입니다.

 

test_integral_gaussian.cpp [다운로드]

 

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include"IntegralTrapezoidal.h"

double mean_gauss = 0.;
double std_gauss = 1.;

double eps_precision = 1e-8;

double func_gaussian(double x) {
    double x_scaled =
        (x - mean_gauss) / std_gauss;
    double ret =
        exp(-0.5 * fabs(x_scaled * x_scaled)) /
        (sqrt(2. * M_PI) * std_gauss);
    return ret;
}

double get_cdf_gaussian(double x) {
    double x_scaled =
        (x - mean_gauss) / std_gauss;
    double ret =
        0.5 * (1. + erf(x_scaled / sqrt(2.)));
    return ret;
}

void make_integral_gaussian(double xmin,
                            double xmax) {
    double (*ptr_integrand)(double);
    ptr_integrand = &func_gaussian;

    fprintf(stdout, "\n");
    fprintf(stdout,
        "    xmin = %e, xmax = %e\n",
        xmin, xmax);

    double integral_num =
        get_integral_trapezoidal(xmin, xmax,
                                 eps_precision,
                                 ptr_integrand);
    fprintf(stdout,
        "    integral (numerical) = %e\n",
        integral_num);

    double integral_cdf =
        get_cdf_gaussian(xmax) -
        get_cdf_gaussian(xmin);
    fprintf(stdout,
        "    integral (using cdf) = %e\n",
        integral_cdf);
}

int main(int argc, char *argv[]) {
    log_integral_trapezoidal_ = NULL;

    fprintf(stdout, "Intgration of gaussian\n");
    fprintf(stdout,
        "  mean = %e, standard deviation = %e\n",
        mean_gauss, std_gauss);

    make_integral_gaussian(-1., 1.);

    make_integral_gaussian(-2., 2.);

    make_integral_gaussian(-3., 3.);

    return 0;
}

 

가우스 함수를 정의한 func_gaussian 함수와 더불어, 가우스 함수의 누적분포함수를 계산하는 get_cdf_gaussian 함수가 있는데요. 누적분포함수를 계산하기 위해 오차함수 (error function) erf를 사용했고, 이는 C언어 수학 함수 라이브러리에도 정의되어 있습니다.

 

코드를 간결하게 하기 위해서 make_integral_gaussian이라는 함수를 추가로 정의했습니다. 이 함수는 사다리꼴 공식과 누적분포함수로부터 구한 적분값을 integral_num과 integral_cdf라는 변수에 각각 저장한 다음, 이들을 적분 범위와 함께 출력하는 기능을 가지고 있죠. main 함수에서 적분구간을 바꿔가면서 make_integral_gaussian 함수를 호출하고 있습니다.

 

GNU C/C++ 컴파일러를 사용하는 경우, 다음과 같이 컴파일 및 빌드가 가능합니다.

 

• 컴파일
  gcc IntegralTrapezoidal.c -c
  g++ test_integral_gaussian.cpp -c
• 링크
  g++ test_integral_gaussian.o IntegralTrapezoidal.o -lm -o [실행파일 이름]

 

프로그램을 실행시켜보면 주어진 범위에 대해서 가우스 함수를 적분한 근사값을 얻을 수 있게 되는데요. 사다리꼴 규칙을 적용한 근사값과 누적분포함수로부터 구한 적분값이 상당히 근접해 있다는 것을 확인할 수 있습니다.

 

 

---

 

C언어 또는 C++ 소스 코드로부터 프로그램을 만드는 과정에 대해서는 다음 포스팅에 더 자세한 내용이 소개되어 있습니다.

 

C/C++ 코드가 프로그램이 되는 과정

 

여기서는 적분 대상이 되는 함수 자체를 매개변수로 넘겨주기 위해서 함수 포인터를 사용하고 있는데요. 이 방법을 통해 함수를 다른 함수의 매개변수로 사용하는 방법에 대한 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

C/C++ 에서 함수를 매개변수로 사용하기

 

C언어로 정의된 함수를 C++ 프로그램에서 사용하기 위해 필요한 extern "C" 코드블록 및 __cplusplus 식별자에 대해서는 다음 포스팅에 더 자세한 내용이 소개되어 있습니다.

 

C언어와 C++를 조합한 프로그래밍

 

여기서는 가우스 함수의 정적분에 대해서 다루었지만, 정적분을 통해 원주율의 근사값을 구하는 것도 가능합니다.

 

C/C++ 원주율 계산 프로그램