cms.integ_gauss_legendre のソースコード

import numpy as np
from math import exp, sqrt, sin, cos, pi

"""
  Numerial integration by Gauss-Legendre method
"""


# define function to be integrated


def func(x):
    return exp(x)

#    return sqrt(1.0 - x*x)

# define analytical integration of f(x)


def F(x):
    return exp(x)

#   return pi / 2.0

#===================
# parameters
#===================
xmin, xmax = -1.0, 1.0

xfrac = [
    [], # 0-point, not dfined
    [], # 1-point, not dfined
    [-1.0 / sqrt(3), + 1.0 / sqrt(3)],           # 2-point order
    [-sqrt(3.0/5.0), 0.0, +sqrt(3.0/5.0)],       # 3-point order
    [-0.861136311594052, -0.339981043584856,      
        +0.339981043584856, +0.861136311594052], # 4-point order
    [-0.906179845938663992797626878299, -0.538469310105683091036314420700, 0.0, 
        +0.538469310105683091036314420700, +0.906179845938663992797626878299] # 5-point order
    ]
weight = [
    [], # 0-point, not dfined
    [], # 1-point, not dfined
    [1.0, 1.0],                  # 2-point order
    [5.0/9.0, 8.0/9.0, 5.0/9.0], # 3-point order
    [0.347854845137453, 0.652145154862546,
        0.652145154862546, 0.347854845137453], # 4-point order
    [0.236926885056189087514264040719, 0.478628670499366468041291514835, 0.568888888888888888888888888888, 
        0.478628670499366468041291514835, 0.236926885056189087514264040719] # 5-point order
	]



def integ_GL(np, func, xmin, xmax):
    px = xfrac[np]
    pw = weight[np]
    S = 0.0
    n = len(px)
    xh = (xmin + xmax) / 2.0
    h = xmax - xmin
    for j in range(n):
        x = xh + px[j] * h / 2.0;
        S += pw[j] * func(x) * h;
    S *= 0.5
    return S


#===================
# main routine
#===================


def main():
    print("Numerical integration using by Gauss-Legendre method")

    fx = func(xmin)
    Fx = F(xmax) - F(xmin)
    print("")
    print("Analytical values:")
    print("  f({})={}".format(xmin, fx))
    print("  integ(f)[{}, {}]={}".format(xmin, xmax, Fx))
    print("")

    for np in [3, 4, 5]:
        print("{}-point formulat".format(np))
        px = xfrac[np]
        pw = weight[np]
        print("  x=", px)
        print("  w=", pw)
        S = integ_GL(np, func, xmin, xmax)
        print("  S={}".format(S))
        print("")


if __name__ == '__main__':
    main()