数値積分の実装と比較：台形法からガウス求積法まで

はじめに

解析的に積分できない関数や、離散的なデータ点しか手元にない関数を積分したい場面は多くあります。そうしたとき、定積分 $\int_a^b f(x)dx$ は数値的に近似するしかありません。

近似のやり方には、区間を等間隔に刻んで足すだけの素朴なものから、評価点の位置まで最適化する巧妙なものまで幅があります。ここでは代表的な4つの手法を取り上げ、同じ積分を計算させて精度と計算時間を比べます。

長方形法（矩形法）
台形法
シンプソン法
ガウス求積法

対象読者:

学部レベルの理工系の学生
実務で数値計算を扱うエンジニア

記事のポイント:

4つの手法の考え方を、近似に使う多項式の次数と対応づけて説明します。
各手法をPythonで実装し、同じ積分を計算して精度と計算時間を実測します。
実測結果をもとに、要求精度に応じた使い分けの目安を示します。

各手法の理論と実装

長方形法（矩形法）

長方形法は、最も基本的な数値積分法です。積分区間を等間隔に分割し、各小区間の関数値をその区間の左端点（または右端点、中央点）での値で代表させ、長方形の面積で近似します。

\int_a^b f(x)dx \approx \sum_{i=0}^{n-1} f(x_i)\Delta x

ここで、

$\Delta x = \frac{b-a}{n}$ は区間の幅
$x_i = a + i\Delta x$ は各区間の左端点

を表します。各区間で関数を定数で置き換えていることになります。

NumPyで書くと数行で済みます。

from typing import Callable
import numpy as np

def rectangle_method(f: Callable[[float], float], a: float, b: float, n: int) -> float:
    """長方形法による数値積分

    Parameters
    ----------
    f : Callable[[float], float]
        積分する関数
    a : float
        積分区間の下限
    b : float
        積分区間の上限
    n : int
        分割数

    Returns
    -------
    float
        積分値の近似値
    """
    dx = (b - a) / n  # 区間幅
    x = np.linspace(a, b - dx, n)  # 各区間の左端点
    return dx * np.sum(f(x))  # 長方形の面積の総和

台形法

台形法では、各小区間で関数を直線で近似し、台形の面積として積分値を計算します。長方形法よりも精度が高く、実装も比較的容易です。

\int_a^b f(x)dx \approx \frac{\Delta x}{2}\left[f(a) + 2\sum_{i=1}^{n-1}f(x_i) + f(b)\right]

この公式は、各区間で関数を線形関数（一次関数）で補間することから導かれます。積分区間の端点 $a$ 、 $b$ での関数値は1回ずつ、内部の点 $x_i$ での関数値は2回ずつ使われます。隣り合う台形が内部の点を共有するためです。

実装では、この重みの違いをそのまま式にします。

def trapezoidal_method(f: Callable[[float], float], a: float, b: float, n: int) -> float:
    """台形法による数値積分

    Parameters
    ----------
    f : Callable[[float], float]
        積分する関数
    a : float
        積分区間の下限
    b : float
        積分区間の上限
    n : int
        分割数

    Returns
    -------
    float
        積分値の近似値
    """
    dx = (b - a) / n  # 区間幅
    x = np.linspace(a, b, n + 1)  # 分点（端点を含む）
    # 端点の寄与は1回、内部の点の寄与は2回
    return dx / 2 * (f(x[0]) + 2 * np.sum(f(x[1:-1])) + f(x[-1]))

シンプソン法

シンプソン法は、各区間で関数を2次関数で近似します。これにより、長方形法や台形法よりも高次の精度が得られます。具体的には、各小区間をさらに2等分し、3点を通る2次関数で補間します。

\int_a^b f(x)dx \approx \frac{\Delta x}{3}\left[f(a) + 4\sum_{i=1,3,5}^{n-1}f(x_i) + 2\sum_{i=2,4,6}^{n-2}f(x_i) + f(b)\right]

この公式は、ラグランジュ補間多項式を用いて導出されます。3点を通る2次関数で関数を近似し、その積分を計算することで得られます。2次以下の多項式に対しては厳密な積分値を返し、奇数番目の点と偶数番目の点で異なる重み（4:2:1の比率）を使います。隣り合う2区間（3点）で1つの2次関数を作る構造上、分割数 $n$ は偶数でなければなりません。

コードでは、奇数番目と偶数番目の点をスライスで分けて足します。

def simpson_method(f: Callable[[float], float], a: float, b: float, n: int) -> float:
    """シンプソン法による数値積分

    Parameters
    ----------
    f : Callable[[float], float]
        積分する関数
    a : float
        積分区間の下限
    b : float
        積分区間の上限
    n : int
        分割数（偶数に調整されます）

    Returns
    -------
    float
        積分値の近似値

    Notes
    -----
    シンプソン法は2次以下の多項式に対して厳密な積分値を与えます。
    分割数は偶数である必要があるため、奇数が指定された場合は
    自動的に偶数に調整されます。
    """
    if n % 2 != 0:
        n += 1  # nを偶数に調整
    dx = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)  # 関数値を一度に計算
    # 奇数点(係数4)と偶数点(係数2)の和を別々に計算
    return dx / 3 * (y[0] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-1:2]) + y[-1])

ガウス求積法

ガウス求積法は、積分区間の分点 $x_i$ と重み $w_i$ を特別に選ぶことで、同じ評価点数でも近似の精度を大きく引き上げます。理論の土台は直交多項式で、なかでもガウス・ルジャンドル求積法はルジャンドル多項式を用います。

\int_a^b f(x)dx \approx \frac{b-a}{2}\sum_{i=1}^n w_i f\left(\frac{b-a}{2}x_i + \frac{b+a}{2}\right)

ここで、 $x_i$ と $w_i$ はガウス・ルジャンドル求積の分点と重みです。 $n$ 点のガウス求積法は、次数 $2n-1$ 以下の多項式に対して厳密な積分値を与えます。分点は等間隔ではなく、精度が最大になる位置（ルジャンドル多項式の零点）に置かれるため、少ない評価点数で高い精度が出ます。

分点と重みはNumPyの leggauss で得られます。

def gauss_quadrature(f: Callable[[float], float], a: float, b: float, n: int) -> float:
    """ガウス求積法による数値積分

    Parameters
    ----------
    f : Callable[[float], float]
        積分する関数
    a : float
        積分区間の下限
    b : float
        積分区間の上限
    n : int
        分点の数（最大50に制限されます）

    Returns
    -------
    float
        積分値の近似値

    Notes
    -----
    この実装では、ガウス・ルジャンドル求積を使用しています。
    n点のガウス求積法は、2n-1次以下の多項式に対して
    厳密な積分値を与えます。
    """
    n = min(n, 50)  # 精度と計算時間のバランスを考慮
    # ガウス・ルジャンドル求積の分点と重みを取得
    x, w = np.polynomial.legendre.leggauss(n)
    # [-1,1]から[a,b]への変換
    t = 0.5 * (b - a) * x + 0.5 * (b + a)
    return 0.5 * (b - a) * np.sum(w * f(t))

同じ積分を4手法で計算する

4つの実装に同じ積分を計算させます。テスト関数は $f(x) = \sin(x)$ 、区間は $[0, \pi]$ で、真の値は2です。

分割数 $n = 100$ での各手法の結果：

長方形法  : 1.9998355039
台形法   : 1.9998355039
シンプソン法: 2.0000000108
ガウス求積法: 2.0000000000
真の値   : 2.0000000000

長方形法と台形法が同じ値になっているのは偶然ではありません。この積分では両端で $\sin(0) = \sin(\pi) = 0$ となり、2つの手法の差である端点の扱いが結果に効かないためです。

誤差が分割数 $n$ とともにどう減るかは、理論的に次のように知られています。

長方形法： $O(1/n)$
台形法： $O(1/n^2)$
シンプソン法： $O(1/n^4)$
ガウス求積法：指数的収束 $O(e^{-cn})$

収束率の差は、近似に使う多項式の次数から来ます。定数で近似する長方形法、直線で近似する台形法、2次関数で近似するシンプソン法の順に速くなり、分点の位置まで選べるガウス求積法は指数的に収束します。

誤差は分割数とともにどう減るか

真の積分値との誤差を、分割数を変えながら計測しました。

ガウス求積法の誤差は他の手法と桁違いに小さく、少ない分割数でも十分な精度が出ています。シンプソン法がそれに続き、台形法と長方形法は同程度ですが、台形法の方が若干小さい誤差を保ちます。理論の収束率どおりの序列です。

計算時間はどこで差がつくか

同じ条件で計算時間も計測しました。

分割数が増えれば計算時間も増えます。長方形法と台形法が最も速く、両者はほぼ同じです。シンプソン法はそれより若干長くなります。ガウス求積法だけは、関数評価の前に分点と重みの計算が入るため、他の手法より時間がかかります。精度あたりのコストで見れば依然として有利ですが、分点数を毎回計算し直す使い方では、この前処理が無視できません。

まとめ

4つの手法の特徴を表に整理します。

手法	近似	精度	計算コスト	特徴
長方形法	定数関数近似	$O(1/n)$	低	実装が容易だが、精度は低い。
台形法	線形近似	$O(1/n^2)$	低	実装が簡単で、長方形法より高精度。
シンプソン法	2次関数近似	$O(1/n^4)$	中	高い精度が得られるが、計算コストは若干増加。分割数は偶数である必要がある。
ガウス求積法	最適な分点選択	$O(e^{-cn})$	高	指数関数的な収束を示し、少ない評価回数で高精度な結果が得られる。ただし、分点と重みの計算が必要。

使い分けは要求精度で決まります。相対誤差 $10^{-3}$ 程度で足りるなら台形法で十分で、実装も検証も簡単です。 $10^{-6}$ 程度まで必要ならシンプソン法、 $10^{-10}$ 以上の高精度が必要で積分区間が固定ならガウス求積法を選びます。数値積分は機械学習の内部でも随所に現れ、変分推論の周辺尤度や強化学習の期待収益の計算がその例です。ただし次元が高くなると、ここで扱った求積法よりモンテカルロ法系の手法が使われることが多くなります。