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

Contents

はじめに
各手法の理論と実装
実装と比較
精度の比較
計算時間の比較
まとめ

はじめに

数値積分は、関数を解析的に積分することが難しい場合や、データ点からの数値的な積分が必要な場合に役立つ手法です。関数 $f(x)$ の定積分 $\int_a^b f(x)dx$ を近似的に計算するために使用されます。

この記事を書いたひと

デジタルリアクタ合同会社代表
望月優輝 / Yuki Mochizuki

機械学習・統計、数値計算などの領域を軸としたソフトウェアエンジニアリングを専門としています。スタートアップからグローバル企業まで、さまざまなスケールの企業にて、事業価値に直結する計算システムを設計・構築してきました。主に機械学習の応用分野において、出版・発表、特許取得等の実績があり、また、IT戦略やデータベース技術等の情報処理に関する専門領域の国家資格を複数保有しています。九州大学理学部卒業。日本ITストラテジスト協会正会員。

この記事では、以下の代表的な数値積分手法について、理論とPythonによる実装を比較します。

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

対象読者:

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

記事のポイント:

各数値積分手法の理論的な背景を、数式を交えて解説します。
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$ は各区間の左端点

を表します。この手法は、各区間において関数を定数関数で近似することに対応します。

Pythonによる実装：

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回ずつ使用されることに注意してください。

Pythonによる実装：

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の比率）。
分割数 $n$ は偶数である必要があります。

Pythonによる実装：

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$ 以下の多項式に対して厳密な積分値を与えます。
分点の位置が等間隔ではなく、最適な位置に配置されます。
少ない評価点数で高い精度が得られます。

Pythonによる実装：

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))

実装と比較

これらの手法をPythonで実装し、比較します。テスト関数として $f(x) = \sin(x)$ を使用し、区間 $[0, \pi]$ での積分を計算します（真の値は2）。

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

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

各手法の理論的な収束率は以下の通りです。

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

ここで、 $n$ は分割数を表します。これらの収束率は、各手法の近似の次数に直接関係しています。長方形法は定数関数で近似、台形法は線形関数で近似、シンプソン法は2次関数で近似するため、それぞれ異なる収束率を示します。特にガウス求積法は、積分区間の分点を最適に選ぶことで、指数関数的な収束率を実現しています。

精度の比較

各手法の精度を比較するため、真の積分値との誤差を計算し、可視化しました。

グラフから、以下のことがわかります。

ガウス求積法が最も高精度であり、少ない分割数でも良好な結果が得られます。
シンプソン法は台形法や長方形法よりも高い精度を示します。
長方形法と台形法は同程度の精度ですが、台形法の方が若干優れています。

計算時間の比較

各手法の計算時間を比較したグラフです。

計算時間に関して、以下の特徴が見られます。

分割数が増加するにつれて、一般的に計算時間も増加します。
長方形法と台形法は最も計算が速く、ほぼ同じ計算時間です。
シンプソン法は若干計算時間が長くなります。
ガウス求積法は分点と重みの計算が必要なため、他の手法よりも計算時間がかかります。

まとめ

各手法の特徴をまとめると以下のようになります。

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

実務では、要求される精度と計算時間のバランスを考慮して、適切な手法を選択することが重要です。以下は、手法選択の一般的な推奨事項です。

高速な計算が必要で、相対誤差 $10^{-3}$ 程度で十分な場合は、台形法を選択します。
相対誤差 $10^{-6}$ 程度の精度が必要な場合は、シンプソン法を選択します。
相対誤差 $10^{-10}$ 以上の高精度が必要で、積分区間が固定の場合は、ガウス求積法を選択します。

近年では、機械学習の分野で、より複雑な関数の積分や、高次元の積分を効率的に行うための研究が進んでいます。例えば、変分推論における周辺尤度の計算や、強化学習における期待収益の計算などに、様々な数値積分手法が応用されています。

お気軽にご相談ください！

弊社のデジタル技術を活用し、貴社のDXをサポートします。

基本的な設計や実装支援はもちろん、機械学習・データ分析・3D計算などの高度な技術が求められる場面でも、最適なソリューションをご提案します。

初回相談は無料で受け付けております。困っていること、是非お聞かせください。

サービスの詳細はこちら

お問い合わせフォームへ

はじめに

各手法の理論と実装

長方形法（矩形法）

台形法

シンプソン法

ガウス求積法

実装と比較

精度の比較

計算時間の比較

まとめ

お気軽にご相談ください！

関連投稿: