import numpy as np
from math import exp, sqrt, pi
# 既存の自作モジュール（パスが通っていることを前提としています）
from integ_gauss_legendre import integ_GL
"""
井本・森口・高澤 (I.M.T.) 公式を用いた数値積分モジュール。
概要:
井本・森口・高澤 (I.M.T.) 公式を実装し、指定された関数を数値積分します。
詳細説明:
I.M.T.変換は、積分区間の端点で特異性を持つ関数を、変換によって滑らかな関数にすることで、
高精度な数値積分を可能にする手法です。
このモジュールでは、特に`sqrt(1 - x^2)`のような区間端点で微分不可能となる関数に対して、
ガウス・ルジャンドル積分と組み合わせて使用します。
モジュールの主要な機能は`main`関数にまとめられており、
`nblock`の数を変えながら積分の精度を確認することができます。
関連リンク:
:doc:`integ_imt_usage`
"""
# define function to be integrated
[ドキュメント]
def func(x):
"""
被積分関数 sqrt(1 - x^2) を定義します。
概要:
与えられた`x`に対して`sqrt(1 - x*x)`を計算します。
詳細説明:
この関数は、積分される対象の具体的な式を表現します。
円の方程式の一部であり、`x`が`[-1, 1]`の範囲で実数を返します。
:param x: float 積分変数の値。
:returns: float `sqrt(1 - x*x)` の計算結果。
"""
return sqrt(1.0 - x*x)
# I.M.T. weight constants
Q_VAL = 0.00702985840661 # int(exp(-1/t - 1/(1-t)) by Gauss-Legendre method
[ドキュメント]
def IMTfunc(x):
"""
I.M.T.変換における重み関数を生成します。
概要:
I.M.T.変換で使用される`exp(-1/x - 1/(1-x)) / Q_VAL`の重み関数を計算します。
詳細説明:
この関数は、I.M.T.変換の核となる部分であり、特異点近傍での関数の振る舞いを調整し、
積分区間全体で関数を滑らかにする役割を果たします。
`Q_VAL`は重み関数の正規化定数で、この関数の`[0, 1]`区間での積分値が1になるように調整します。
:param x: float 変換された積分変数の値 (通常は `[0, 1]` の範囲)。
:returns: float I.M.T.重み関数の値。
"""
return exp(-1.0 / x - 1.0 / (1.0 - x)) / Q_VAL
#===================
# parameters (Default values)
#===================
# Gauss-Legendre integration parameters to convert x to u
n_gl_order_def = 5
ngl_def = 10
# Integration parameters
xmin_def, xmax_def = -1.0, 1.0
nblock_max_def = 10
[ドキュメント]
def main():
"""
I.M.T.公式を用いて`sqrt(1 - x^2)`の数値積分を実行し、結果を表示します。
概要:
`nblock`の数を変化させながら、I.M.T.公式による数値積分を実行し、
解析解との誤差を比較して表示します。
詳細説明:
この関数は、以下のステップで数値積分を行います。
1. パラメータをデフォルト値で初期化します。
2. 解析解と被積分関数`func(x)`の特定点での値を表示します。
3. `nblock`を2から`nblock_max_def`まで増加させながら以下の処理を繰り返します。
a. 積分区間`[0, 1]`を`nblock`個の小区間に分割します。
b. 各小区間において、I.M.T.変換`phi(ui)`によって元の積分変数を求めます。
この`phi(ui)`の計算には、`integ_GL`（ガウス・ルジャンドル積分）を用いて`IMTfunc`を積分します。
c. 変換された積分変数`xi`を用いて`func(xi)`を評価し、`IMTfunc`の重みを乗じて積算します。
d. 最後にスケーリングを行い、数値積分結果`S`と解析解`Fx`との誤差を表示します。
:returns: None
"""
# パラメータのローカル変数化
n_gl_order = n_gl_order_def
ngl = ngl_def
xmin, xmax = xmin_def, xmax_def
nblock_max = nblock_max_def
print("Numerical integration using by IMT formula")
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 = (1.0 - 0.0) / nblock
S = 0.0
for i in range(1, nblock):
# ui is in the converted range [0, 1]
ui = i * h
# x' = phi(ui) is in [0,1]
phiui = 0.0
for j in range(ngl):
hgl = (ui - 0.0) / ngl
x0_gl = 0.0 + j * hgl
x1_gl = x0_gl + hgl
# Gauss-Legendre積分を用いて座標変換を行う
phiui += integ_GL(n_gl_order, IMTfunc, x0_gl, x1_gl)
# obtain the original x from x' = phi(ui)
xi = xmin + (xmax - xmin) * phiui
yi = func(xi)
phi1 = exp(-1.0 / ui - 1.0 / (1.0 - ui)) / Q_VAL
S += yi * phi1
# Final scaling
S *= (xmax - xmin) / nblock
print("{}\t{}\t{}".format(nblock, S, S - Fx))
if __name__ == "__main__":
main()