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

"""
倍精度指数関数型数値積分モジュール。

二重指数関数型変換 (Double Exponential Transformation, DE変換) を用いて、
特定区間の定積分を数値的に評価します。
本モジュールは、特に被積分関数 `sqrt(1 - x^2)` の `[-1, 1]` 区間における
数値積分を実行し、解析解と比較します。

:doc:`integ_doubleexp_usage`
"""
import numpy as np
from math import exp, sqrt, pi, tanh, sinh, cosh
# 既存の自作モジュール（必要に応じてパスを通してください）
# from integ_gauss_legendre import integ_GL

# define function to be integrated


def func(x):
    """
    積分対象の関数 f(x) = sqrt(1 - x^2) を定義します。

    :param x: float型の独立変数。
    :returns: float型の `sqrt(1 - x^2)` の計算結果。
    """
    return sqrt(1.0 - x*x)


#===================
# parameters (Default values)
#===================
# Integration parameters
xmin_def, xmax_def = -1.0, 1.0
umin_def, umax_def = -2.0, 2.0
nblock_max_def = 32



def main():
    """
    倍精度指数関数型数値積分を実行し、結果を出力するメイン関数です。

    `[-1, 1]` の区間で `sqrt(1 - x^2)` を積分します。
    各ブロック数 (nblock) における数値積分結果と解析解との誤差を表示します。
    """
    # パラメータのローカル変数化
    xmin, xmax = xmin_def, xmax_def
    umin, umax = umin_def, umax_def
    nblock_max = nblock_max_def

    print("Numerical integration using by double exponential integra for [-1, 1]")

    fx = func(xmin)
    # Analytical integration of sqrt(1-x^2) from -1 to 1 is pi/2
    Fx = pi / 2.0
    
    print("")
    print("Analytical values:")
    print("  f({})={}".format(xmin, fx))
    print("  integ(f)[{}, {}]={}".format(xmin, xmax, Fx))
    print("")

    print("nBlock\tS\tError")
    for nblock in range(2, nblock_max + 1):
        h = (umax - umin) / nblock
        S = 0.0
        for i in range(0, nblock + 1):
            # ui is in the converted range [-u, u]
            ui = umin + i * h
            # x' = phi(ui) is in [-1, 1]
            phiui = tanh(pi / 2.0 * sinh(ui))

            # obtain the original x from x' = phi(ui)
            xi = xmin + (phiui + 1.0) / 2.0 * (xmax - xmin)
            yi = func(xi)
            
            # Derivative of phi(ui) for the weight
            coshsinh = cosh(pi / 2.0 * sinh(ui))
            coshnh = cosh(ui)
            phi1 = coshnh / (coshsinh * coshsinh)
            S += yi * phi1

        # Final scaling
        S *= pi / 2.0 * h * (xmax - xmin) / 2.0
        print("{}\t{}\t{}".format(nblock, S, S - Fx))


if __name__ == "__main__":
    main()