Numerial integration by Gauss-Legendre method
Download script from .\integ_gauss_legendre.py
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()
