Numerical integration/Gauss-Legendre Quadrature 高斯勒让德积分

https://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature#Fortran

==============

Gauss 积分部分

from numpy.polynomial.legendre import leggauss
x, w = leggauss(5)
print('x=\n',x)
print('w=\n',w)

x=
[-0.90617985 -0.53846931 0. 0.53846931 0.90617985]
w=
[0.23692689 0.47862867 0.56888889 0.47862867 0.23692689]

2-nodes  形函数-刚度矩阵

import numpy as np
from numpy.polynomial.legendre import leggauss
import matplotlib.pyplot as plt

N1 = lambda x: -x/2 +1/2
N2 = lambda x:  x/2 +1/2
dN1 = lambda x: -1/2
dN2 = lambda x:  1/2
x=np.linspace(-1,1,200)
"""
Numerical quadrature (Gauss Legendre) 
L 
NB_E = 1
Element matrix
"""
def gauss_legendre_quad(f,n,a,b):
    x,w = leggauss(n)
    sum_ = 0
    for k in range(len(x)):
        sum_ += w[k] * f(0.5*(b-a)*x[k]+0.5*(b+a))
    return 0.5*(b-a)*sum_

# Element matrix
L = 1
nb_e=2
h = L/nb_e       # element_size
nb_dof = nb_e +1 # number of dof

N = [N1, N2]
dN = [dN1, dN2]
K_e=np.zeros((2,2))
M_e=np.zeros((2,2))

for i in range(len(dN)):
    for j in range(len(dN)):
        f = lambda x:dN[i](x)*dN[j](x)
        K_e[i, j] = (2/h)*gauss_legendre_quad(f, 5, -1, 1)
        g = lambda x:N[i](x)*N[j](x)
        M_e[i, j]  = (2/h)*gauss_legendre_quad(g, 5, -1, 1)

# Assembly the element matrix
K = np.zeros((nb_dof,nb_dof))
for i in range(nb_e):
    K[i:i+2,i:i+2] +=K_e