cms.integration.integ_scipy_quad のソースコード

"""
SciPyの積分機能の使用例を示すスクリプト。

このスクリプトは、`scipy.integrate` モジュールを用いて、
様々な数値積分手法をデモンストレーションします。
具体的には、解析的積分（QUADPACK）、ロンバーグ積分、広義積分、
そして離散データの累積積分（`cumulative_trapezoid`）の利用方法を示します。

関連リンク: :doc:`integ_scipy_quad_usage`
"""
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt



def computePi(x):
    """
    円周率計算のための被積分関数を定義します。

    この関数は、`integrate.quad` を用いて円周率を計算する際の被積分関数 `4/(1+x^2)` を返します。
    0から1まで積分すると円周率（π）になります。

    :param x: float: 積分変数の値。
    :returns: float: 被積分関数の評価値。
    """
    return 4/(1+x**2)




def main():
    """
    スクリプトの主要な処理を実行します。

    この関数は、`scipy.integrate` モジュールの様々な積分手法をデモンストレーションします。
    具体的には、以下の処理を実行します。

    1.  `integrate.quad` を用いた解析的な積分（QUADPACK）で円周率を計算します。
    2.  `integrate.romberg` を用いたロンバーグ積分で円周率を計算します。
    3.  `integrate.quad` を用いた広義積分（0から無限大まで）を実行します。
    4.  `integrate.cumulative_trapezoid` を用いて離散データの累積積分を行い、
        その結果を解析解と比較してグラフ表示します。

    :returns: None
    """
    # 1. 解析的な積分 (QUADPACK)
    res = integrate.quad(computePi, 0.0, 1.0)
    print("(res, error)=", res)

    # 2. ロンバーグ積分
    res = integrate.romberg(computePi, 0, 1)
    print("(res, error)=", res)

    # 3. 広義積分 (0 to infinity)
    res = integrate.quad(lambda x: x**5*np.exp(-x)*np.sin(x), 0, np.inf)
    print("(res, error)=", res)

    # ---------------------------------------------------------
    # 離散データの累積積分
    # ---------------------------------------------------------
    x = np.linspace(0, 2, num=2**4+1)
    y = x**4

    # cumtrapz ではなく cumulative_trapezoid を使用
    # initial=0 とすることで、入力 x と同じ長さの配列が返ります
    y_int = integrate.cumulative_trapezoid(y, x, initial=0)

    # グラフ表示
    plt.plot(x, y_int, 'ro', label='cumulative_trapezoid')
    plt.plot(x, y[0] + 0.2 * x**5, 'b-', label='analytical (x^5 / 5)')
    plt.title("Comparison of Cumulative Integration")
    plt.legend()
    plt.show()


if __name__ == "__main__":
    main()