当前位置：首页 > news >正文

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

news 来源：原创 2024/9/20 9:43:15

原文链接：https://www.cnblogs.com/aksoam/p/18343674
更多精彩，关注博客园主页，不断学习！不断进步！

我的主页

csdn很少看私信，有事请b站私信

博客园主页-发文字笔记-常用

有限元鹰的主页内容：

ABAQUS数值模拟相关
Python科学计算
开源框架，编程学习笔记

哔哩哔哩主页-发视频-常用

FE-有限元鹰的个人空间内容：

模拟案例
网格划分
游戏视频，及其他搬运视频

文章目录

我的主页
- 博客园主页-发文字笔记-常用
- 哔哩哔哩主页-发视频-常用
二维高斯积分
三维高斯积分
验证

本文只给出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}")