Pythonで学ぶ数値計算のアルゴリズムと実装 第4回：高度な数値積分：二重指数型（DE）公式と広義積分

こんにちは、JS2IIUです。
前回の第3回では、数値積分の基本である中点則、台形則、そしてシンプソン則について学びました。これらは「関数が十分に滑らかで、有限の範囲で定義されている」場合には非常に強力な手法です。しかし、実世界のデータ解析や物理シミュレーション、あるいは統計モデルの計算においては、標準的な手法では太刀打ちできない「扱いにくい積分」が頻繁に登場します。

例えば、積分範囲の端で値が無限大に発散してしまう関数や、無限遠（\(-\infty\) から \(\infty\)）まで続く積分、あるいは非常に激しく振動する周期関数の積分などです。これらは「広義積分」と呼ばれ、単純な分割による近似では計算が終わらなかったり、著しく精度が落ちたりします。

今回は、日本で生まれ、今や世界中で「最強の数値積分アルゴリズム」の一つとして知られる「二重指数型（DE）公式（Double Exponential formula）」を中心に、高度な数値積分の世界を探求していきましょう。今回もよろしくお願いします。

1. 標準的な積分則が苦手なもの：広義積分の壁

例えば、次のような積分を考えてみてください。
$$\int_{0}^{1} \frac{1}{\sqrt{x}} dx$$
この関数 \(f(x) = 1/\sqrt{x}\) は、\(x=0\) において値が無限大に発散します。数学的にはこの積分の値は $2$ と定まりますが、台形則やシンプソン則で計算しようとすると、\(x=0\) での評価ができずにエラーになるか、あるいは非常に大きな誤差が生じてしまいます。

また、正規分布の期待値を求める際などに現れる無限区間の積分 \(\int_{-\infty}^{\infty} f(x) dx\) も、どこで計算を打ち切れば良いのかという問題が常に付きまといます。

これらの「扱いにくい問題」を解決するために、座標系そのものを変換して「扱いやすい形」に作り変えてしまおう、というのが高度な数値積分の基本的な考え方です。

2. 二重指数型（DE）公式

二重指数型（DE）公式は、1974年に高橋秀俊教授と森正武教授によって提案された手法です。この手法の核心は、「変数の置き換え」にあります。

変数変換のアイデア

有限区間 \([-1, 1]\) の積分を、変数変換 \(x = \phi(t)\) によって無限区間 \((-\infty, \infty)\) の積分に変換します。このとき、変換関数 \(\phi(t)\) として以下のような「二重指数関数的」に変化するものを選びます。
$$x = \tanh\left( \frac{\pi}{2} \sinh t \right)$$

この変換を行うと、被積分関数は \(t \to \pm \infty\) で驚異的な速さ（\(\exp(-\exp(t))\) のオーダー）で \(0\) に収束するようになります。関数の値がこれほど速く \(0\) に近づくのであれば、無限に続く \(t\) の範囲を適当なところ（例えば \(-4\) から \(4\) 程度）で打ち切っても、コンピュータの限界（浮動小数点誤差）以下の無視できる誤差しか生じません。

結果として、非常に難しいはずの積分が、「等間隔な台形則」を適用するだけで極めて高精度に計算できてしまうのです。

3. PythonによるDE公式の実装

それでは、PyTorchを用いてDE公式を実装してみましょう。このアルゴリズムは、サンプリング点 \(x_i\)と重み \(w_i\) をあらかじめ計算しておき、それらを関数値に掛け合わせる「重み付き和」の形で実行します。

import torch
import numpy as np

def de_integration(f, a, b, n=100, h=0.1):
    """
    二重指数型(DE)公式を用いた数値積分
    Args:
        f (callable): 積分対象の関数
        a (float): 積分範囲の下端
        b (float): 積分範囲の上端
        n (int): サンプリング点の数 (片側)
        h (float): ステップサイズ
    Returns:
        torch.Tensor: 積分値
    """
    # tの範囲を生成 (-n*h から n*h まで)
    k = torch.arange(-n, n + 1)
    t = k * h
    
    # 変数変換 phi(t) = tanh(pi/2 * sinh(t))
    # これは [-1, 1] の範囲を取る
    sinh_t = torch.sinh(t)
    u = 0.5 * torch.pi * sinh_t
    phi = torch.tanh(u)
    
    # phi'(t) = d/dt [tanh(pi/2 * sinh(t))]
    #         = (pi/2) * cosh(t) * sech^2(pi/2 * sinh(t))
    cosh_t = torch.cosh(t)
    cosh_u = torch.cosh(u)
    phi_prime = 0.5 * torch.pi * cosh_t / (cosh_u * cosh_u)
    
    # 区間 [-1, 1] から [a, b] への線形変換
    # x = (b-a)/2 * phi + (b+a)/2
    # dx/dphi = (b-a)/2
    x_mapped = 0.5 * (b - a) * phi + 0.5 * (b + a)
    
    # ヤコビアン: dx/dt = dx/dphi * dphi/dt = (b-a)/2 * phi'(t)
    jacobian = 0.5 * (b - a) * phi_prime
    
    # 関数値の計算（特異点の安全な処理）
    # x=0での評価を避けるため、極端に小さい値をクリップ
    x_safe = torch.clamp(x_mapped, min=1e-300)
    
    try:
        f_values = f(x_safe)
    except:
        # エラーが発生した場合の処理
        f_values = torch.zeros_like(x_safe)
        for i, x_val in enumerate(x_safe):
            try:
                f_values[i] = f(x_val)
            except:
                f_values[i] = 0.0
    
    # inf/nanのチェックと除外
    valid_mask = torch.isfinite(f_values) & torch.isfinite(jacobian)
    
    # 台形公式: I ≈ h * Σ f(x(t)) * (dx/dt)
    integrand = f_values * jacobian
    integrand = torch.where(valid_mask, integrand, torch.zeros_like(integrand))
    
    integral = h * torch.sum(integrand)
    
    return integral


# テスト1：端点に特異点を持つ関数 f(x) = 1/sqrt(x) の [0, 1] での積分
# 理論値 = 2.0
def singular_func(x):
    return 1.0 / torch.sqrt(x)

print("=== Test 1: ∫[0,1] 1/√x dx (理論値: 2.0) ===")
a, b = 0.0, 1.0
result = de_integration(singular_func, a, b, n=150, h=0.05)
print(f"DE Formula result: {result.item():.12f}")
print(f"Theoretical value: 2.0")
print(f"Error: {abs(result.item() - 2.0):.12e}")

# テスト2：通常の関数での検証
def normal_func(x):
    return torch.exp(-x)

print("\n=== Test 2: ∫[0,1] e^(-x) dx (理論値: 0.632120558...) ===")
theoretical = 1 - np.exp(-1)
result2 = de_integration(normal_func, 0.0, 1.0, n=100, h=0.1)
print(f"DE Formula result: {result2.item():.12f}")
print(f"Theoretical value: {theoretical:.12f}")
print(f"Error: {abs(result2.item() - theoretical):.12e}")

# テスト3：両端特異点
def double_singular(x):
    return 1.0 / torch.sqrt(x * (1 - x))

print("\n=== Test 3: ∫[0,1] 1/√(x(1-x)) dx (理論値: π) ===")
result3 = de_integration(double_singular, 0.0, 1.0, n=200, h=0.05)
print(f"DE Formula result: {result3.item():.12f}")
print(f"Theoretical value: {np.pi:.12f}")
print(f"Error: {abs(result3.item() - np.pi):.12e}")

import torch
import numpy as np

def de_integration(f, a, b, n=100, h=0.1):
    """
    二重指数型(DE)公式を用いた数値積分
    Args:
        f (callable): 積分対象の関数
        a (float): 積分範囲の下端
        b (float): 積分範囲の上端
        n (int): サンプリング点の数 (片側)
        h (float): ステップサイズ
    Returns:
        torch.Tensor: 積分値
    """
    # tの範囲を生成 (-n*h から n*h まで)
    k = torch.arange(-n, n + 1)
    t = k * h
    
    # 変数変換 phi(t) = tanh(pi/2 * sinh(t))
    # これは [-1, 1] の範囲を取る
    sinh_t = torch.sinh(t)
    u = 0.5 * torch.pi * sinh_t
    phi = torch.tanh(u)
    
    # phi'(t) = d/dt [tanh(pi/2 * sinh(t))]
    #         = (pi/2) * cosh(t) * sech^2(pi/2 * sinh(t))
    cosh_t = torch.cosh(t)
    cosh_u = torch.cosh(u)
    phi_prime = 0.5 * torch.pi * cosh_t / (cosh_u * cosh_u)
    
    # 区間 [-1, 1] から [a, b] への線形変換
    # x = (b-a)/2 * phi + (b+a)/2
    # dx/dphi = (b-a)/2
    x_mapped = 0.5 * (b - a) * phi + 0.5 * (b + a)
    
    # ヤコビアン: dx/dt = dx/dphi * dphi/dt = (b-a)/2 * phi'(t)
    jacobian = 0.5 * (b - a) * phi_prime
    
    # 関数値の計算（特異点の安全な処理）
    # x=0での評価を避けるため、極端に小さい値をクリップ
    x_safe = torch.clamp(x_mapped, min=1e-300)
    
    try:
        f_values = f(x_safe)
    except:
        # エラーが発生した場合の処理
        f_values = torch.zeros_like(x_safe)
        for i, x_val in enumerate(x_safe):
            try:
                f_values[i] = f(x_val)
            except:
                f_values[i] = 0.0
    
    # inf/nanのチェックと除外
    valid_mask = torch.isfinite(f_values) & torch.isfinite(jacobian)
    
    # 台形公式: I ≈ h * Σ f(x(t)) * (dx/dt)
    integrand = f_values * jacobian
    integrand = torch.where(valid_mask, integrand, torch.zeros_like(integrand))
    
    integral = h * torch.sum(integrand)
    
    return integral


# テスト1：端点に特異点を持つ関数 f(x) = 1/sqrt(x) の [0, 1] での積分
# 理論値 = 2.0
def singular_func(x):
    return 1.0 / torch.sqrt(x)

print("=== Test 1: ∫[0,1] 1/√x dx (理論値: 2.0) ===")
a, b = 0.0, 1.0
result = de_integration(singular_func, a, b, n=150, h=0.05)
print(f"DE Formula result: {result.item():.12f}")
print(f"Theoretical value: 2.0")
print(f"Error: {abs(result.item() - 2.0):.12e}")

# テスト2：通常の関数での検証
def normal_func(x):
    return torch.exp(-x)

print("\n=== Test 2: ∫[0,1] e^(-x) dx (理論値: 0.632120558...) ===")
theoretical = 1 - np.exp(-1)
result2 = de_integration(normal_func, 0.0, 1.0, n=100, h=0.1)
print(f"DE Formula result: {result2.item():.12f}")
print(f"Theoretical value: {theoretical:.12f}")
print(f"Error: {abs(result2.item() - theoretical):.12e}")

# テスト3：両端特異点
def double_singular(x):
    return 1.0 / torch.sqrt(x * (1 - x))

print("\n=== Test 3: ∫[0,1] 1/√(x(1-x)) dx (理論値: π) ===")
result3 = de_integration(double_singular, 0.0, 1.0, n=200, h=0.05)
print(f"DE Formula result: {result3.item():.12f}")
print(f"Theoretical value: {np.pi:.12f}")
print(f"Error: {abs(result3.item() - np.pi):.12e}")

=== Test 1: ∫[0,1] 1/√x dx (理論値: 2.0) ===
DE Formula result: 1.999716758728
Theoretical value: 2.0
Error: 2.832412719727e-04

=== Test 2: ∫[0,1] e^(-x) dx (理論値: 0.632120558...) ===
DE Formula result: 0.632120549679
Theoretical value: 0.632120558829
Error: 9.149755175741e-09

=== Test 3: ∫[0,1] 1/√(x(1-x)) dx (理論値: π) ===
DE Formula result: 3.140845775604
Theoretical value: 3.141592653590
Error: 7.468779855451e-04

=== Test 1: ∫[0,1] 1/√x dx (理論値: 2.0) ===
DE Formula result: 1.999716758728
Theoretical value: 2.0
Error: 2.832412719727e-04

=== Test 2: ∫[0,1] e^(-x) dx (理論値: 0.632120558...) ===
DE Formula result: 0.632120549679
Theoretical value: 0.632120558829
Error: 9.149755175741e-09

=== Test 3: ∫[0,1] 1/√(x(1-x)) dx (理論値: π) ===
DE Formula result: 3.140845775604
Theoretical value: 3.141592653590
Error: 7.468779855451e-04

DE公式のコード解説

1. 元の問題

このコードが近似しているのは，区間 \([a,b]\) 上の定積分
$$
I \;=\; \int_a^b f(x)\,dx
$$

です。

2. 区間変換 \([a,b] \to [-1,1]\)

まず，線形変換を考えます：
$$
x = \frac{b-a}{2}\,\xi + \frac{b+a}{2},
\quad \xi \in [-1,1]
$$
このとき
$$
dx = \frac{b-a}{2}\, d\xi
$$
したがって，
$$
I
= \frac{b-a}{2}
\int_{-1}^{1}
f!\left(
\frac{b-a}{2}\,\xi + \frac{b+a}{2}
\right)
d\xi
$$

3. DE（tanh–sinh）変換 \([-1,1] \to (-\infty,\infty)\)

DE公式では，次の変数変換を用います：
$$
\xi = \phi(t)
= \tanh!\left(\frac{\pi}{2}\sinh t\right),
\quad t \in (-\infty,\infty)
$$

導関数（ヤコビアン）
$$
\phi’(t)
= \frac{d}{dt}
\tanh!\left(\frac{\pi}{2}\sinh t\right)
= \frac{\pi}{2}
\cosh t \,
\operatorname{sech}^2!\left(\frac{\pi}{2}\sinh t\right)
$$

4. 積分の完全な変換

上記を組み合わせると，
$$
I
= \frac{b-a}{2}
\int_{-\infty}^{\infty}
f!\left(
\frac{b-a}{2}\phi(t) + \frac{b+a}{2}
\right)
\phi’(t)\,dt
$$
これが コードが連続的に表している積分式です。

5. 数値積分（台形公式）

コードでは \(t\) を等間隔で離散化しています：
$$
t_k = k h, \quad k = -n, \ldots, n
$$
これに対し，台形公式（実質的には単純和）で近似：
$$
I \;\approx\;
h
\sum_{k=-n}^{n}
f!\left(
\frac{b-a}{2}\phi(t_k) + \frac{b+a}{2}
\right)
\left[
\frac{b-a}{2}\phi’(t_k)
\right]

$$

6. コードとの対応表

コードのポイント

• 変数変換の威力: 通常の台形則では $x=0$ を評価しようとしてエラーになりますが、DE公式では \(x\) が \(0\) や \(1\) に「限りなく近づくが、決して到達しない」ため、特異点を安全に、かつ精密に扱うことができます。
• ベクトル化: torch.sinh や torch.tanh をテンソル全体に一括適用することで、Pythonの遅いループを回避しています。

4. 周期関数と急減少関数の積分

DE公式のバリエーションは、無限区間の積分 \(\int_{0}^{\infty} f(x) dx\) や周期関数（\(\sin, \cos\) を含む積分）にも存在します。

特に物理学や信号処理では、急速に減衰しながら激しく振動する関数の積分が求められます。このような場合、通常のサンプリングでは振動の「山」と「谷」を正確に捉えきれません。DE公式の考え方を応用すると、振動の周期に合わせてサンプリング点を自動的に最適化するような戦略が可能になります。

AIの分野においても、ガウス過程のカーネル計算や、ベイズ推論における複雑な正規化定数の計算など、この種の「急減少関数」の積分は精度を左右する重要な要素です。

5. 実践：SciPyの quad との比較

数値計算の専門家を目指すなら、自作アルゴリズムだけでなく、標準ライブラリが内部で何をしているかを知っておくべきです。SciPyの scipy.integrate.quad は、内部で「QUADPACK」という歴史あるFortranライブラリを呼び出しています。

quad は、関数の変化が激しい場所を自動的に細かく分割する「適応型ガウス・クロンロッド則（Adaptive Gauss-Kronrod rule）」を採用しています。

from scipy import integrate

# SciPyによる計算
# quadは特異点があってもある程度自動で対処してくれるが、
# 非常に厳しい特異点の場合はDE公式の方が有利なこともある
val, err = integrate.quad(lambda x: 1.0/np.sqrt(x), 0, 1)

print(f"SciPy quad result: {val:.12f}")
print(f"Estimated error: {err:.12e}")

from scipy import integrate

# SciPyによる計算
# quadは特異点があってもある程度自動で対処してくれるが、
# 非常に厳しい特異点の場合はDE公式の方が有利なこともある
val, err = integrate.quad(lambda x: 1.0/np.sqrt(x), 0, 1)

print(f"SciPy quad result: {val:.12f}")
print(f"Estimated error: {err:.12e}")

どちらを使うべきか？

• SciPy quad: ほとんどのケースで最適な第一選択です。適応型アルゴリズムなので、関数の形状を事前に知らなくても高い精度を出してくれます。
• 自作DE公式: 積分の端点に非常に強い特異点がある場合や、計算を完全にPyTorchの計算グラフ内に載せて、積分値そのもので微分（勾配計算）を行いたい場合（微分可能プログラミング）に威力を発揮します。

6. 計算量とオーダーの視点

DE公式の驚くべき点は、その効率性です。精度（有効桁数）を $D$ 桁出すために必要なサンプリング点数 $N$ は、おおよそ $N \propto D \log D$ 程度であることが知られています。これは、他の手法が $N \propto 10^D$ のような膨大な計算を必要とするのと比較して、圧倒的に高速であることを意味します。

数値計算における「オーダー $O(n)$」の議論は通常、データの数に対する計算量を指しますが、積分の世界では「必要な精度に対する計算量」という視点が欠かせません。

まとめ

第4回では、標準的な積分手法では解決できない問題を突破するための「高度な数値積分」について解説しました。

• 広義積分は、端点の発散や無限区間を含む、数値計算上の難所である。
• 二重指数型（DE）公式は、変数変換によって被積分関数を急激に減衰させ、等間隔なサンプリングでの高精度計算を可能にする。
• 変数変換というアイデアは、数値計算における「座標系の変更による問題の単純化」という普遍的な戦略の代表例である。
• SciPyの quad は非常に強力な汎用ツールだが、特定の条件下やディープラーニングとの統合においては、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/