from scipy import integrate
def Gaussion_Quadrature1D(f,a,b):
x_y_w=np.array([[-0.9491079123427585,0.1294849661688697],
[-0.7415311855993945,0.2797053914892766],
[-0.4058451513773972,0.3818300505051189],
[ 0.0000000000000000,0.4179591836734694],
[ 0.4058451513773972,0.3818300505051189],
[ 0.7415311855993945,0.2797053914892766],
[ 0.9491079123427585,0.1294849661688697]])
S=0
for i in range(0,np.size(x_y_w,0)):
xq = x_y_w[i,0]
wq = x_y_w[i,1]
S=S+wq*f(((b-a)/2)*xq+(b+a)/2)
return S*(b-a)/2
def Gaussion_Quadrature2D1(f,a,b,c,d):
x_y_w=np.array([[-0.9491079123427585,0.1294849661688697],
[-0.7415311855993945,0.2797053914892766],
[-0.4058451513773972,0.3818300505051189],
[ 0.0000000000000000,0.4179591836734694],
[ 0.4058451513773972,0.3818300505051189],
[ 0.7415311855993945,0.2797053914892766],
[ 0.9491079123427585,0.1294849661688697]])
S=0
for i in range(0,np.size(x_y_w,0)):
xq = x_y_w[i,0]
xw = x_y_w[i,1]
for j in range(0,np.size(x_y_w,0)):
yq = x_y_w[j,0]
yw = x_y_w[j,1]
S=S+xw*yw*f(((b-a)/2)*xq+(b+a)/2,((d-c)/2)*yq+(d+c)/2)
return S*((b-a)/2)*((d-c)/2)
def Gaussion_Quadrature2D2(f,a,b,c,d):
x_y_w=np.array([[-0.9491079123427585,0.1294849661688697],
[-0.7415311855993945,0.2797053914892766],
[-0.4058451513773972,0.3818300505051189],
[ 0.0000000000000000,0.4179591836734694],
[ 0.4058451513773972,0.3818300505051189],
[ 0.7415311855993945,0.2797053914892766],
[ 0.9491079123427585,0.1294849661688697]])
S=0
for i in range(0,np.size(x_y_w,0)):
xq = x_y_w[i,0]
xw = x_y_w[i,1]
yq = x_y_w[:,0]
yw = x_y_w[:,1]
S=S+xw*yw@f(((b-a)/2)*xq+(b+a)/2,((d-c)/2)*yq+(d+c)/2)
return S*((b-a)/2)*((d-c)/2)
def Gaussion_Quadrature2Dv(f,h):
x_y_w=np.array([[-0.9491079123427585,0.1294849661688697],
[-0.7415311855993945,0.2797053914892766],
[-0.4058451513773972,0.3818300505051189],
[ 0.0000000000000000,0.4179591836734694],
[ 0.4058451513773972,0.3818300505051189],
[ 0.7415311855993945,0.2797053914892766],
[ 0.9491079123427585,0.1294849661688697]])
S=0
a=0
b=h
d=h
for i in range(0,np.size(x_y_w,0)):
xq = x_y_w[i,0]
xw = x_y_w[i,1]
kasi=((b-a)/2)*xq+(b+a)/2
c=(h**2-kasi**2)**(0.5)
for j in range(0,np.size(x_y_w,0)):
yq = x_y_w[j,0]
yw = x_y_w[j,1]
eta=((d-c)/2)*yq+(d+c)/2
S=S+xw*yw*f(kasi,eta)*(d-c)/2
return S*(b-a)/2
def Gaussion_Quadrature3D(F,a,b,c,d,e,f):
x_y_w=np.array([[-0.9491079123427585,0.1294849661688697],
[-0.7415311855993945,0.2797053914892766],
[-0.4058451513773972,0.3818300505051189],
[ 0.0000000000000000,0.4179591836734694],
[ 0.4058451513773972,0.3818300505051189],
[ 0.7415311855993945,0.2797053914892766],
[ 0.9491079123427585,0.1294849661688697]])
S=0
N=np.size(x_y_w,0)
for k in range(0,N):
zq = x_y_w[k,0]
zw = x_y_w[k,1]
for i in range(0,N):
xq = x_y_w[i,0]
xw = x_y_w[i,1]
yq = x_y_w[:,0]
yw = x_y_w[:,1]
S=S+zw*xw*yw@F(((b-a)/2)*xq+(b+a)/2,((d-c)/2)*yq+(d+c)/2,((f-e)/2)*zq+(f+e)/2)
return S*((b-a)/2)*((d-c)/2)*((f-e)/2)
if(1):
f=lambda x:x**2
print(Gaussion_Quadrature1D(f,-1,1))
print(integrate.quad(f,-1,1))
f=lambda x,y:x**2*y**2
print(Gaussion_Quadrature2D2(f,-1,1,-1,1))
print(integrate.dblquad(f,-1,1,-1,1))
if(1):
f=lambda x,y,z:x**2*y**2*z**2
print(Gaussion_Quadrature3D(f,-1,1,-1,1,-1,1))
print(integrate.tplquad(f,-1,1,-1,1,-1,1))
0.6666666666666665
(0.6666666666666666, 7.401486830834376e-15)
0.4444444444444443
(0.44444444444444436, 7.337339523254554e-15)
0.29629629629629617
(0.2962962962962963, 7.273748168439867e-15)