C++로 이중 적분 계산하기 – 심프슨 1/3 법칙 완벽 가이드

문제 개요

변수 x의 하한과 상한, 변수 y의 하한과 상한, 그리고 각 변수에 대한 구간 간격(step)이 주어졌을 때, 이중 적분을 계산하여 그 결과를 출력하는 프로그램을 작성하는 것이 이번 글의 목표입니다.

접근 방식

이 프로그램에서는 심프슨 1/3 법칙(Simpson's 1/3 Rule)을 사용합니다. 심프슨 1/3 법칙은 구간을 짝수 개의 소구간으로 나눈 뒤, 인접한 세 점을 지나는 2차 다항식으로 함수를 근사하여 정적분 값을 구하는 수치 적분 기법입니다. 이 방법을 x축 방향과 y축 방향으로 각각 적용하면 이중 적분을 계산할 수 있습니다.

x와 y의 상한·하한 값, 그리고 각 변수의 구간 간격을 입력받습니다.
x와 y 방향 모두에 심프슨 1/3 법칙을 적용하여 이중 적분을 계산합니다.
본격적인 계산에 앞서 아래와 같은 함수값 표를 먼저 작성합니다.
각 행에 심프슨 1/3 규칙을 적용해 첫 번째 적분을 구하고, 같은 과정을 한 번 더 반복하여 이중 적분을 완성합니다.
최종 결과를 출력합니다.

C++로 이중 적분 계산하기 – 심프슨 1/3 법칙 완벽 가이드

예시

입력 -:
x의 간격 = 1.2
y의 간격 = 0.54
x의 하한 = 1.3
x의 상한 = 2.1
y의 하한 = 1.0
y의 상한 = 2.1

출력 -:
이중 적분 값 : 2.1

알고리즘

시작
1단계 -> 적분 대상 함수의 값을 계산하는 함수 선언
    float fun(float x, float y)
    return pow(pow(x, 4) + pow(y, 5), 0.5)
2단계 -> 이중 적분 값을 구하는 함수 선언
    float doubleIntegral(float step_x, float step_y, float lower_x, float upper_x, float lower_y, float upper_y)
    int n1, n2 선언
    float arr[50][50], arr_2[50], result 선언
    n1 = (upper_x - lower_x) / step_x + 1
    n2 = (upper_y - lower_y) / step_y + 1
    반복 For int i = 0, i < n1, ++i
        반복 For int j = 0, j < n2, ++j
            arr[i][j] = fun(lower_x + i * step_x, lower_y + j * step_y)
        끝
    끝
    반복 For int i = 0, i < n1, ++i
        arr_2[i] = 0
        반복 For int j = 0, j < n2, ++j
            IF (j == 0 || j == n2 - 1)
                arr_2[i] += arr[i][j]
            ELSE IF (j % 2 == 0)
                arr_2[i] += 2 * arr[i][j]
            ELSE
                arr_2[i] += 4 * arr[i][j]
            끝
        끝
        arr_2[i] *= (step_y / 3)
    끝
    result = 0
    반복 For int i = 0, i < n1, ++i
        IF (i == 0 || i == n1 - 1)
            result += arr_2[i]
        ELSE IF (i % 2 == 0)
            result += 2 * arr_2[i]
        ELSE
            result += 4 * arr_2[i]
        끝
        result *= (step_x / 3)
    끝
    return result
3단계 -> main() 함수에서
    x의 간격: float step_x = 1.2
    y의 간격: float step_y = 0.54
    x의 하한: float lower_x = 1.3
    x의 상한: float upper_x = 2.1
    y의 하한: float lower_y = 1.0
    y의 상한: float upper_y = 2.1
    doubleIntegral(step_x, step_y, lower_x, upper_x, lower_y, upper_y) 호출
종료

C++ 예제 코드

#include <bits/stdc++.h>
using namespace std;

float fun(float x, float y) {
    return pow(pow(x, 4) + pow(y, 5), 0.5);
}

// 이중 적분 값을 구하는 함수
float doubleIntegral(float step_x, float step_y, float lower_x, float upper_x, float lower_y, float upper_y) {
    int n1, n2;
    float arr[50][50], arr_2[50], result;
    n1 = (upper_x - lower_x) / step_x + 1;
    n2 = (upper_y - lower_y) / step_y + 1;

    // 각 격자점에서의 함수값 계산
    for (int i = 0; i < n1; ++i) {
        for (int j = 0; j < n2; ++j) {
            arr[i][j] = fun(lower_x + i * step_x, lower_y + j * step_y);
        }
    }

    // y 방향으로 심프슨 1/3 법칙 적용
    for (int i = 0; i < n1; ++i) {
        arr_2[i] = 0;
        for (int j = 0; j < n2; ++j) {
            if (j == 0 || j == n2 - 1)
                arr_2[i] += arr[i][j];
            else if (j % 2 == 0)
                arr_2[i] += 2 * arr[i][j];
            else
                arr_2[i] += 4 * arr[i][j];
        }
        arr_2[i] *= (step_y / 3);
    }

    // x 방향으로 심프슨 1/3 법칙 적용
    result = 0;
    for (int i = 0; i < n1; ++i) {
        if (i == 0 || i == n1 - 1)
            result += arr_2[i];
        else if (i % 2 == 0)
            result += 2 * arr_2[i];
        else
            result += 4 * arr_2[i];
    }
    result *= (step_x / 3);
    return result;
}

int main() {
    float step_x = 1.2;   // x의 간격
    float step_y = 0.54;  // y의 간격
    float lower_x = 1.3;  // x의 하한
    float upper_x = 2.1;  // x의 상한
    float lower_y = 1.0;  // y의 하한
    float upper_y = 2.1;  // y의 상한

    cout << "이중 적분 값 : " << doubleIntegral(step_x, step_y, lower_x, upper_x, lower_y, upper_y);
    return 0;
}

출력 결과

이중 적분 값 : 2.1

정리

이 프로그램은 먼저 x와 y의 격자점마다 함수값을 2차원 배열에 저장한 뒤, y 방향으로 심프슨 1/3 법칙을 적용해 각 행의 적분 값을 구합니다. 이어서 그 결과들에 다시 x 방향으로 심프슨 1/3 법칙을 적용하여 최종적인 이중 적분 값을 얻습니다. 심프슨 1/3 법칙을 사용하면 단순 사다리꼴 공식보다 훨씬 높은 정확도로 수치 적분을 수행할 수 있다는 점이 큰 장점입니다.