目录
本文只给出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