Python：使用scipy和自定义高斯积分计算多重积分

最新推荐文章于 2024-05-26 16:43:38 发布

自我o救赎

最新推荐文章于 2024-05-26 16:43:38 发布

阅读量369

点赞数

分类专栏： Numerical Method For PDE 文章标签： python scipy

本文链接：https://blog.csdn.net/y15520833229/article/details/130475676

版权

Numerical Method For PDE 专栏收录该内容

12 篇文章 0 订阅

订阅专栏

from scipy import integrate

def Gaussion_Quadrature1D(f,a,b):#定义一元函数的7点高斯积分区间在-1 到 1上
    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)
    #print(S*(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)
    #print(S*(b-a)/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]
         #for j in range(0,np.size(x_y_w,0)):
         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)
    #print(S*(b-a)/2)
    return S*((b-a)/2)*((d-c)/2)
#定义一个变限2D积分 -高斯积分  第二个积分号  的下极限为可变函数
#这个数值积分不可以随便  调用
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
         #print(h**2,kasi**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]
             #for j in range(0,np.size(x_y_w,0)):
             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)
    #print(S*(b-a)/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)