【数值计算方法】2&3维高斯积分的python实现

合集 - 数值计算(18)

1.【数学公式】mathtype和word2016集成2023-04-192.【MathType】word2016数学公式编号2023-04-193.【Matlab】基于KDtree的最近邻搜索和范围搜索2023-07-304.【Matlab】判断点和多面体位置关系的两种方法实现2023-07-295.《空间三角面片对相交判断算法》的matlab实现_ 0.2微秒2023-11-116.Mpmath库-学习笔记2024-07-027.有限元方法[Matlab]-笔记2024-04-138.结构动力学教材-学习笔记2024-06-199.复合材料力学基础及有限元分析2024-05-3010.【数值计算方法】数值积分&微分-python实现2024-07-30

11.【数值计算方法】2&3维高斯积分的python实现2024-08-05

12.【数值计算方法】线性方程组的迭代解法2024-08-0613.【数值计算方法】线性方程组迭代算法的Python实现2024-08-0714.【数值计算方法】线性方程组的迭代解法-数值实验2024-08-0715.【数值计算方法】非线性方程求根2024-08-0816.【数值计算方法】非线性方程求根-数值实验2024-08-1217.【数值计算方法】常微分方程初值问题的数值解2024-08-1418.Note_Fem边界条件的处理和numpy实现的四种方法03-08

收起

目录

• 二维高斯积分
• 三维高斯积分
• 验证

本文只给出pythont实现和例题,数学推导见【数值计算方法】数值积分&微分-python实现 - FE-有限元鹰 - 博客园

二维高斯积分

python实现二维高斯积分:

 
def Integ2dGuassLegendre(f,lowLimit:List[float]=[-1,-1], 
                            upLimit:List[float]=[1,1], 
                            n:int=3)->float:
    """给定积分区域[lowLimit,upLimit]和高斯积分点个数n(>=1),计算二维高斯-勒让德积分公式"""
    a,b,c,d=lowLimit[0],upLimit[0],lowLimit[1],upLimit[1]
    if n<=0:
        raise ValueError("高斯-勒让德积分时,n必须大于0")
    if n==1:
        return 4*f(0,0)
    if a==-1 and b==1 and c==-1 and d==1:
        # 标准型积分
        #计算权重和积分点位置
        x_is,w_is=legendre_gauss_points_and_weights(n)
        y_js,w_js=legendre_gauss_points_and_weights(n)
 
        return np.sum([w_is[ind_x]*w_js[ind_y]*f(xi,yj) for ind_x,xi in enumerate(x_is) for ind_y,yj in enumerate(y_js) ])
    else:
        # 非标准型积分,积分区域:[a,b]x[c,d]
        xt1=lambda t1: 0.5*(b-a)*t1+0.5*(b+a)
        yt2=lambda t2: 0.5*(d-c)*t2+0.5*(d+c)
        
        t1_is,w_is=legendre_gauss_points_and_weights(n)
        t2_js,w_js=legendre_gauss_points_and_weights(n)
        
        return np.sum([w_is[ind_x]*w_js[ind_y]*f(xt1(t1i),yt2(t2j))*(b-a)*(d-c)*0.25 for ind_x,t1i in enumerate(t1_is) for ind_y,t2j in enumerate(t2_js) ])
 

三维高斯积分

python实现:

def Integ3dGuassLegendre(f,lowLimit:List[float]=[-1,-1,-1], 
                            upLimit:List[float]=[1,1,1], 
                            n:int=3)->float:
    """给定积分区域[lowLimit,upLimit]和高斯积分点个数n(>=1),计算二维高斯-勒让德积分公式"""
    a,b=lowLimit[0],upLimit[0]
    c,d=lowLimit[1],upLimit[1]
    g,h=lowLimit[2],upLimit[2]
    if n<=0:
        raise ValueError("高斯-勒让德积分时,n必须大于0")
    if n==1:
        return 8*f(0,0)
    if a==-1 and b==1 and c==-1 and d==1:
        # 标准型积分
        #计算权重和积分点位置
        x_is,w_is=legendre_gauss_points_and_weights(n)
        y_js,w_js=legendre_gauss_points_and_weights(n)
        z_js,w_ks=legendre_gauss_points_and_weights(n)
 
        return np.sum([w_is[ind_x]*w_js[ind_y]*w_ks[ind_z]*f(xi,yj,zk) for ind_x,xi in enumerate(x_is) for ind_y,yj in enumerate(y_js) for ind_z,zk in enumerate(z_js)])
    else:
        # 非标准型积分,积分区域:[a,b]x[c,d]x[g,h]
        xt1=lambda t1: 0.5*(b-a)*t1+0.5*(b+a)
        yt2=lambda t2: 0.5*(d-c)*t2+0.5*(d+c)
        zt3=lambda t3: 0.5*(h-g)*t3+0.5*(h+g)
        
        t1_is,w_is=legendre_gauss_points_and_weights(n)
        t2_js,w_js=legendre_gauss_points_and_weights(n)
        t3_ks,w_ks=legendre_gauss_points_and_weights(n)
        
        return np.sum([w_is[ind_x]*w_js[ind_y]*w_ks[ind_z]*f(xt1(t1i),yt2(t2j),zt3(t3k))*(b-a)*(d-c)*(h-g)*0.125 for ind_x,t1i in enumerate(t1_is) for ind_y,t2j in enumerate(t2_js) for ind_z,t3k in enumerate(t3_ks)])

验证

from formu_lib import *
import numpy as np
from scipy.integrate import dblquad
# 定义被积函数
def integrand(x, y):
    return np.exp(x*x+y*y)
 
# 计算二重积分
result, error = dblquad(integrand, -1, 1, lambda x: -1, lambda x: 1)
print("numpy 二重积分结果:", result)
ans1=Integ2dGuassLegendre(integrand,[-1, -1],[1, 1],n=5)
print(f"本地二重积分结果:{ans1}")
 
from scipy.integrate import tplquad
# 定义被积函数
def integrand3(x, y, z):
    return np.exp(x*x+y*y+z*z)
 
# 计算三重积分
result3, error = tplquad(integrand3, -1, 1, lambda x: -1, lambda x: 1, lambda x, y: -1, lambda x, y: 1)
ans2=Integ3dGuassLegendre(integrand3,[-1,-1,-1],[1,1,1],n=5)
 
print("numpy三重积分结果:", result3)
print(f"本地三重积分结果:{ans2}")

输出:

• numpy 二重积分结果: 8.557400519221307
• 本地二重积分结果:8.557173227239266
• numpy三重积分结果: 25.03299361973213
• 本地三重积分结果:25.03199627931168