C++で二重積分を計算するプログラム｜シンプソン1/3則による数値積分の実装

変数xの下限・上限、変数yの下限・上限、そしてx・yそれぞれの刻み幅（ステップ幅）が与えられたとき、二重積分を数値的に計算し、その結果を表示するのが本記事のテーマです。

入出力の例

入力:
xの刻み幅 = 1.2
yの刻み幅 = 0.54
xの下限 = 1.3
xの上限 = 2.1
yの下限 = 1.0
yの上限 = 2.1

出力:
double integration is : 2.1

計算のアプローチ

本プログラムでは、以下の手順で二重積分を求めます。

• xとyの上限・下限の値に加えて、x・yそれぞれの刻み幅を入力として受け取ります。
• 二重積分の計算にはシンプソン1/3則（Simpson 1/3 rule）を採用します。
• 計算に先立ち、x方向の分割点を行・y方向の分割点を列とする二次元の表（配列）に、各格子点における被積分関数の値を格納します。
• まずy方向（各行）に対してシンプソン1/3則を適用して内側の積分を求め、続いてx方向にもう一度同じ規則を適用することで二重積分を完成させます。
• 最終的な積分値を出力します。

シンプソン1/3則とは

シンプソン1/3則は、積分区間を偶数個の小区間に分割し、各区間を放物線（2次関数）で近似する数値積分の手法です。各格子点の重みは、両端が1、偶数番目の点が2、奇数番目の点が4となり、合計に「刻み幅 ÷ 3」を乗じることで積分値が得られます。二重積分の場合は、この処理を一方の変数方向に適用した後、もう一方の変数方向に再度適用します。

アルゴリズム

開始
ステップ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 とする
    ループ（int i = 0; i < n1; ++i）:
        ループ（int j = 0; j < n2; ++j）:
            arr[i][j] = fun(lower_x + i * step_x, lower_y + j * step_y) とする
    ループ（int i = 0; i < n1; ++i）:
        arr_2[i] = 0 とする
        ループ（int j = 0; j < n2; ++j）:
            j == 0 または j == n2 - 1 の場合:
                arr_2[i] += arr[i][j]
            それ以外で j % 2 == 0 の場合:
                arr_2[i] += 2 * arr[i][j]
            それ以外の場合:
                arr_2[i] += 4 * arr[i][j]
        arr_2[i] *= (step_y / 3) とする
    result = 0 とする
    ループ（int i = 0; i < n1; ++i）:
        i == 0 または i == n1 - 1 の場合:
            result += arr_2[i]
        それ以外で i % 2 == 0 の場合:
            result += 2 * arr_2[i]
        それ以外の場合:
            result += 4 * arr_2[i]
    result *= (step_x / 3) とする
    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) を呼び出す
終了

サンプルコード

#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);
        }
    }
    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;
}

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 << "double integration is : " << doubleIntegral(step_x, step_y, lower_x, upper_x, lower_y, upper_y);
    return 0;
}

出力

double integration is : 2.1

このように、シンプソン1/3則をy方向とx方向の2段階で順に適用することで、二重積分の近似値を効率的に求めることができます。刻み幅を小さくするほど分割数が増え、より正確な結果が得られます。