cms.integ_doubleexp のソースコード

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

"""
  Numerial integration by Double Exponential Integration
"""


# 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)

#===================
# parameters
#===================
# Integration parameters
xmin, xmax = -1.0, 1.0
umin, umax = -2.0, 2.0
nblock_max = 32

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

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 = (umax - umin) / nblock;
    S = 0.0;
    for i in range(0, nblock+1):
# ui is in the converted range [-1, 1]
        ui = umin + i * h; 
# x' = phi(ui) is in [0,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);
        coshsinh = cosh(pi / 2.0 * sinh(ui))
        coshnh = cosh(ui)
        phi1 = coshnh / (coshsinh*coshsinh)
        S += yi * phi1

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