cms.integ_imt のソースコード

import numpy as np
from math import exp, sqrt, pi
from integ_gauss_legendre import integ_GL

"""
  Numerial integration by Imoto-Moriguchi-Takazawa (I.M.T.) formula
"""


# define function to be integrated


def func(x):
    return sqrt(1.0 - x*x)


# define analytical integration of f(x)
#def F(x):
#    return exp(x)

Q = 0.00702985840661; # int(exp(-1/t - 1/(1-t) by Gauss-Legendre method


def IMTfunc(x):
    return exp(-1.0 / x - 1.0 / (1.0-x)) / Q;


#===================
# parameters
#===================
# Gauss-Legendre integration parameters to convert x to u
n_gl_order = 5
ngl = 10

# Integration parameters
xmin, xmax = -1.0, 1.0
nblock_max = 10

#===================
# main routine
#===================
print("Numerical integration using by IMT formula")

fx = func(xmin)
#Fx = F(xmax) - F(xmin)
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
            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
        S += yi * phi1

    S *= (xmax - xmin) / nblock
    print("{}\t{}\t{}".format(nblock, S, S - Fx))