Pythonで学ぶ数値計算のアルゴリズムと実装 第3回：積分の基礎：中点則、台形則からシンプソン則まで

Python
2026/01/01
この記事は約19分で読めます。

こんにちは、JS2IIUです。
連載の第1回、第2回を通して、私たちは多項式を用いた関数の近似や、離散的な点から未知の値を推定する補間技術について学んできました。第3回となる今回は、これらの知識を応用して、数値計算の大きな柱の一つである「数値積分（Numerical Integration）」に挑戦します。

数学の授業で習った積分は、原始関数を求めて値を代入する「解析的な解法」が中心でした。しかし、実世界のデータ解析やAIエンジニアリングにおいて扱う関数の多くは、原始関数が簡単に求まるとは限りません。例えば、統計学で頻出する正規分布の累積分布関数や、複雑な物理シミュレーションの結果として得られるデータなどです。

コンピュータにこうした積分の値を計算させるには、積分範囲を小さく分割し、それぞれの領域の面積を近似的に求めて合計する「数値積分」の手法が必要不可欠です。今回は、その基本となる3つの手法を、PythonとPyTorchを用いて深く理解していきましょう。今回もよろしくお願いします。

1. はじめに：数値積分とは何か？

積分 \(\int_{a}^{b} f(x) dx\) を計算することは、グラフ \(y=f(x)\) と \(x\) 軸、および直線 \(x=a, x=b\) で囲まれた領域の「面積」を求めることに他なりません。

コンピュータは連続的な曲線をそのまま扱うことはできないため、積分範囲 \([a, b]\) を \(n\) 個の小さな区間に分割し、それぞれの小さな区間において「計算しやすい図形（長方形、台形、あるいは放物線など）」で代用して面積を計算します。これを「数値求積（Numerical Quadrature）」とも呼びます。

2. 補間型積分則の考え方

第2回で学んだ「多項式補間」を思い出してください。数値積分の基本的な戦略は以下の通りです。

  1. 積分範囲内のいくつかのサンプリング点において関数の値を取得する。
  2. それらの点を通る「多項式」で関数を近似（補間）する。
  3. その多項式を積分する。

この「補間多項式の次数」によって、中点則（0次）、台形則（1次）、シンプソン則（2次）といった手法の名称が決まります。

3. 中点則（Midpoint Rule）

中点則は、各区間の中央における関数の値 \(f(x_{mid})\) を高さとし、区間の幅 \(\Delta x\) を底辺とする長方形の面積を足し合わせる手法です。

中点則のアルゴリズム

積分区間 \([a, b]\) を \(n\) 分割したとき、各区間の幅 \(h = (b-a)/n\) とすると、各区間の面積は \(h \cdot f(a + (i + 0.5)h)\) となります。

PyTorchによる実装

PyTorchを用いることで、多数の区間の計算を一括（ベクトル演算）して高速に行うことができます。

Python
import torch
import matplotlib.pyplot as plt

def midpoint_rule(f, a, b, n):
    """
    中点則を用いて数値積分を実行する

    Args:
        f (callable): 積分対象の関数
        a (float): 積分範囲の下端
        b (float): 積分範囲の上端
        n (int): 分割数

    Returns:
        torch.Tensor: 積分値の近似
    """
    # 区間の幅を計算
    h = (b - a) / n
    # 各区間の中点を生成
    # 0.5h, 1.5h, ..., (n-0.5)h の位置
    midpoints = torch.linspace(a + h/2, b - h/2, n)

    # 各中点での関数値を計算し、合計してhを掛ける
    integral = h * torch.sum(f(midpoints))
    return integral

# テスト用の関数: f(x) = exp(x) (積分範囲 [0, 1])
# 真の値は exp(1) - exp(0) = 1.71828...
def test_func(x):
    return torch.exp(x)

a, b = 0.0, 1.0
n = 100
result = midpoint_rule(test_func, a, b, n)
print(f"Midpoint Rule result (n={n}): {result.item():.8f}")
Plaintext
Midpoint Rule result (n=100): 1.71827459

中点則は非常にシンプルですが、関数が直線（1次式）であれば誤差なしで計算できるという、意外と強力な性質を持っています。

4. 台形則（Trapezoidal Rule）

台形則は、各区間の両端の点 \((x_i, f(x_i))\) と \((x_{i+1}, f(x_{i+1}))\) を直線で結び、台形の面積として近似する手法です。第2回の言葉で言えば「1次多項式補間」による積分です。

台形則のアルゴリズム

台形の面積は「(上底 + 下底) * 高さ / 2」で求められます。積分範囲全体では、両端の点は1回ずつ、それ以外の内部の点は2回ずつ足されることになります。

PyTorchによる実装

PyTorchには torch.trapezoid という組み込み関数も存在しますが、ここでは理解のために基礎的なテンソル操作で実装してみましょう。

Plaintext
def trapezoidal_rule(f, a, b, n):
    """
    台形則を用いて数値積分を実行する
    """
    # 刻み幅の計算
    h = (b - a) / n
    # 評価点の生成（両端を含む n+1 個の点）
    x = torch.linspace(a, b, n + 1)
    y = f(x)

    # 両端の重みは1, 内部の重みは2
    # 積分 = (h/2) * (y[0] + 2*y[1] + ... + 2*y[n-1] + y[n])
    integral = (h / 2.0) * (y[0] + 2 * torch.sum(y[1:-1]) + y[-1])
    return integral

result_trap = trapezoidal_rule(test_func, a, b, n)
print(f"Trapezoidal Rule result (n={n}): {result_trap.item():.8f}")
Plaintext
Trapezoidal Rule result (n=100): 1.71829617

台形則は直感的に理解しやすく、特に周期関数を全周期にわたって積分する場合などに、驚異的な精度を発揮することが知られています。

5. シンプソン則（Simpson’s Rule）

シンプソン則は、各区間（あるいは2つの隣接する区間）を「2次多項式（放物線）」で近似して積分する手法です。中点則や台形則よりも一段高い精度（4次精度）を誇ります。

シンプソン則のアルゴリズム

3つの点を通る放物線の面積を計算すると、重みが「1 : 4 : 1」になるという面白い性質が現れます。区間数が偶数である必要があります。

PyTorchによる実装

スライシングを活用して、重み 4 と 2 の項を効率的に計算します。

Python
def simpsons_rule(f, a, b, n):
    """
    シンプソン則を用いて数値積分を実行する
    注: nは偶数である必要がある
    """
    if n % 2 != 0:
        n += 1 # 奇数の場合は偶数に調整

    h = (b - a) / n
    x = torch.linspace(a, b, n + 1)
    y = f(x)

    # シンプソンの公式: (h/3) * [y_0 + 4*y_1 + 2*y_2 + 4*y_3 + ... + y_n]
    # 奇数番目のインデックス（y_1, y_3, ...）には重み4
    # 偶数番目の内部インデックス（y_2, y_4, ...）には重み2
    sum_odd = torch.sum(y[1:-1:2])
    sum_even = torch.sum(y[2:-1:2])

    integral = (h / 3.0) * (y[0] + 4 * sum_odd + 2 * sum_even + y[-1])
    return integral

result_simp = simpsons_rule(test_func, a, b, n)
print(f"Simpson's Rule result (n={n}): {result_simp.item():.8f}")
Plaintext
Simpson's Rule result (n=10): 6.38911247

シンプソン則は、関数が3次多項式以下であれば、理論上誤差なしで計算が可能です。

6. 精度の比較検証：誤差収束レートの可視化

数値計算において最も重要なのは「計算コスト（分割数 \(n\)）」に対して「どれだけ速く誤差が減るか」という性質です。これを「収束レート」と呼びます。

一般に、分割幅を \(h\) とすると誤差は以下のように評価されます。

  • 中点則・台形則：\(O(h^2)\)（\(h\) を半分にすると誤差は 1/4 になる）
  • シンプソン則：\(O(h^4)\)（\(h\) を半分にすると誤差は 1/16 になる）

実際にグラフを描いて確かめてみましょう。

Python
import numpy as np
import matplotlib.pyplot as plt

# 積分区間
a, b = 0, np.pi
# 真値
true_value = 2.0

def f(x):
    return np.sin(x)

def midpoint_rule(f, a, b, n):
    h = (b - a) / n
    return h * np.sum(f(a + h * (np.arange(n) + 0.5)))

def trapezoidal_rule(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])

def simpson_rule(f, a, b, n):
    if n % 2 == 1:
        n += 1  # nは偶数でなければならない
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return h / 3 * (y[0] + 2 * np.sum(y[2:n:2]) + 4 * np.sum(y[1:n:2]) + y[n])

n_values = [2**i for i in range(1, 11)]
errors_mid = []
errors_trap = []
errors_simp = []

for n in n_values:
    errors_mid.append(abs(midpoint_rule(f, a, b, n) - true_value))
    errors_trap.append(abs(trapezoidal_rule(f, a, b, n) - true_value))
    errors_simp.append(abs(simpson_rule(f, a, b, n) - true_value))

plt.figure(figsize=(8, 6))
plt.loglog(n_values, errors_mid, 'o-', label='Midpoint Rule')
plt.loglog(n_values, errors_trap, 's-', label='Trapezoidal Rule')
plt.loglog(n_values, errors_simp, '^-', label="Simpson's Rule")
plt.xlabel('Number of intervals (n)')
plt.ylabel('Absolute Error')
plt.title('Convergence of Numerical Integration Methods')
plt.legend()
plt.grid(True, which="both", ls="--")
plt.show()

対数グラフにおいて、直線の傾きが急であるほど精度が高い（収束が速い）ことを意味します。シンプソン則がいかに圧倒的な速さで真の値に近づいているかが視覚的に確認できたはずです。

7. 実践的なアドバイス：SciPyでの積分

実務において、自作の実装を使用するのは小規模な計算や学習目的に限られます。大規模、あるいは高い信頼性が求められるプロジェクトでは、SciPy ライブラリの integrate モジュールを使用するのが一般的です。

Python
from scipy import integrate

# SciPyのquad関数（適応型ガウス・クロンロッド則という非常に高度な手法を使用）
val, err = integrate.quad(lambda x: torch.exp(torch.tensor(x)).item(), 0, 1)
print(f"SciPy quad result: {val:.12f} (estimated error: {err:.2e})")
Plaintext
SciPy quad result: 1.718281827670 (estimated error: 4.46e-09)

integrate.quad は、関数の形状に応じて自動的にサンプリング点の位置や密度を調整する「適応型」のアルゴリズムを採用しており、非常に堅牢です。

まとめ

第3回となる今回は、数値積分の基礎となる3つのアルゴリズムを学びました。

  • 中点則は、シンプルながらも実装が容易で、0次近似に基づいています。
  • 台形則は、1次多項式補間に基づき、直感的な面積計算を行います。
  • シンプソン則は、2次多項式補間を利用し、少ない分割数で非常に高い精度を得ることができます。
  • 誤差の収束を理解することは、計算リソースを最適化する上で極めて重要です。

積分は、AI分野においても「期待値の計算」や「事後分布の正規化」など、至る所で登場します。今回学んだ「連続なものを離散的なサンプルの集まりとして捉える」という視点は、今後の学習の大きな助けになるでしょう。

次回は、さらに高度な積分の世界へ進みます。通常の積分則では太刀打ちできない「無限遠までの積分」や「端点で値が発散する関数」を鮮やかに解く「二重指数型（DE）公式」について解説します。お楽しみに。

最後まで読んでいただきありがとうございます。

連載記事、各回へのリンク

Pythonで学ぶ数値計算のアルゴリズムと実装 第1回：多項式計算とテイラー展開 | アマチュア無線局JS2IIU
https://js2iiu.com/2025/12/31/python_numerical_01/

Pythonで学ぶ数値計算のアルゴリズムと実装 第2回：多項式補間とチェビシェフ補間の戦略 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_02/

Pythonで学ぶ数値計算のアルゴリズムと実装 第3回：積分の基礎：中点則、台形則からシンプソン則まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_03/

Pythonで学ぶ数値計算のアルゴリズムと実装 第4回：高度な数値積分：二重指数型（DE）公式と広義積分 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_04/

Pythonで学ぶ数値計算のアルゴリズムと実装 第5回：計算の「正しさ」を評価する：ノルムと条件数 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numeric_05/

Pythonで学ぶ数値計算のアルゴリズムと実装 第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_06/

Pythonで学ぶ数値計算のアルゴリズムと実装 第8回：スパース性を活かす：帯行列の高速計算アルゴリズム | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerical_08/

Pythonで学ぶ数値計算のアルゴリズムと実装 第9回：最小二乗法とQR分解：ハウスホルダー変換による安定化 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerial_09/

Pythonで学ぶ数値計算のアルゴリズムと実装 第10回：非線形方程式の探索：2分法からニュートン法まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_10/

Pythonで学ぶ数値計算のアルゴリズムと実装 第11回：行列の本質を捉える：固有値問題の数値解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_11/

Pythonで学ぶ数値計算のアルゴリズムと実装 第12回：現象の変化を追う：ルンゲ・クッタ法による常微分方程式の解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_12/

Pythonで学ぶ数値計算のアルゴリズムと実装 第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_13/

コメント

タイトルとURLをコピーしました