第四十四篇高斯勒让德求解多重积分

多重积分

在工程分析中，经常需要在一个面积或体积上对函数进行积分。多重积分的解析方法在有限的情况下是可能的，但在这一篇中使用数值积分去求解。一维的函数积分详见重复牛顿-科特斯积分，重复高斯勒让德积分
考虑两个变量的函数在二维区域R上的积分，如下图所示。函数f(x, y)可以被认为是在R区域上以直角从平面上伸出的第三维空间。
在这里插入图片描述
需要的积分为

通过前面篇章中描述的一维空间函数积分得到，涉及两个变量的二重积分将产生下面形式的积分法则

样本点位于R的平面上，坐标为(xi, yi)， i = 1，…N是函数的值。每个函数的计算都用wi加权，加权之和必须等于R的面积。
显然，在明确地定义不规则区域的积分计算起来有些困难。在实践中，将形状不规则的区域细分为若干更小的简单子区域(如矩形)就足够了，在这些子区域上可以很容易地进行数值积分。然后将每个子区域的解相加，得到整个区域R上的最终近似结果。这种方法类似于前面篇章中提到的“重复”规则。
首先，考虑在平行于笛卡尔坐标方向的矩形区域上的积分，因为这极大地简化了极限的定义。随后，这些概念被扩展到xy平面上一般的四边形和三角形区域的积分。本节中描述的所有方法都可以很容易地推广到三重积分。

矩形区域积分

在下图中，考虑函数f(x, y)在所示的矩形区域的积分。由于矩形的边界与笛卡尔坐标方向平行，变量之间没有关系，可以直接应用前面描述的积分方法。
在这里插入图片描述

牛顿柯特

例如，在每个方向上应用梯形规则将导致4个样本点(n = 4)，如上图所示，矩形的每个角上都有一个样本点，从而产生该规则
在这里插入图片描述
在每个方向上应用Simpson规则，得到如下图6所示的9个样本点(n = 9)，从而得到该法则为

其实可以猜到，上面得到的的加权系数来自于在每一个方向分别考虑的加权系数的简单乘积

计算实例

在每个方向使用梯形法则去计算下面积分
在这里插入图片描述
在一维中，梯形规则有2个样本点，所以在二重积分中，我们将使用n = 2*2 = 4个样本点。由上图知，h = 1, k = 2，有梯形法则的二维公式得到

精确解为15.3333

计算实例

在每个方向使用辛普森法则去计算下面的积分
在这里插入图片描述
在一维中，辛普森法则有3个样本点，所以在二重积分中，使用n = 3**2 = 9个采样点。由上图得，h = 0.5, k = 1，根据三点辛普森法则得到，二维积分得

在这个例子中，在两个方向上都使用辛普森法是低效的。在x方向上用辛普森法则和在y方向上用梯形法则也可以得到精确的解。可以使用n = 3 × 2=6个样本点，得到

高斯-勒让德

高斯-勒让德法则也可以应用于这种类型的多重积分，但必须注意找到样本点的正确位置。一种方法是进行坐标变换，使每个方向上的积分极限变为±1。就可以使用这篇中提到得权值和采样点值。详细的表格对应数值可见高斯勒让德求积分
矩形区域的另一种方法是将样本点在每个坐标方向上的范围中点左右进行加权和比例系数相乘。
考虑下图所示的矩形区域上的两点高斯-勒让德积分。
在这里插入图片描述
在x方向上，样本点坐标为

对应的加权系数w1 = w2 = h。同样在y方向上，样本点为

对应加权系数w1 = w2 = k，则法则为

计算实例

在两个方向使用高斯-勒让德计算下面的积分
在这里插入图片描述
在这个二重积分中，我们使用n = 2*22 = 4个样本点。由上图可知，h = 1, k = 0.5，根据上面得式子，样本点位于

因此得到

精确解为0.8
在这个例子中，两点法则能够在x方向上精确积分，但不能在y方向上。
通过类似的推理，由图下图所示的三点高斯勒让德法则可以得到x中得样本点和加权系数
x方向为
在这里插入图片描述
在y方向为

因此法则将变成

计算实例

在每个方向使用3点高斯-勒让德去计算下面的积分
在这里插入图片描述
在这个二重积分中，我们将使用n = 3**2 = 9个采样点。由上图可知，h = 1, k = 0.5，然后根据样本点坐标值得方程求得样本点为

因此根据积分方程，得到I = 0。8，这是精确解。
在这种情况下，四次项在y上的积分的精确解可以通过三点法则得到，在x方向上使用两点法则，总共有6个样本点。尽管在不同的方向上实现不同的数量积分法则是很容易的，但是大多数数值积分软件包在所有方向上都使用相同的数量法则。

程序如下

其中有一个主程序，和两个子程序，分别为高斯勒让德的样本点值和权值的子程序gauss_legendre;和局部坐标的形状函数的子程序fun_der。

#高斯拉盖尔法则
import numpy as np
import B
import math
ndim=2;nsp=9
if ndim==1:nod=2
elif ndim==2:nod=4
elif ndim==3:nod=8
coord=np.zeros((nod,ndim));der=np.zeros((ndim,nod));fun=np.zeros(nod)
coord[:,0]=(0.0,1.0,4.0,6.0);coord[:,1]=(0.0,2.0,3.0,1.0)
samp=np.zeros((nsp,ndim));wt=np.zeros((nsp))
B.gauss_legendre(samp,wt)
res=0
def f65(point):x=point[0];y=point[1]#;z=point[2]f65=x**2*y**2return f65
for i in range(1,nsp+1):B.fun_der(fun,der,samp,i)res=res+np.linalg.det(np.dot(der,coord))*wt[i-1]*f65(np.dot(np.transpose(coord),fun))
print('高斯-勒让德法则的多重积分')
print('维数',ndim)
for i in range(1,nod+1):print('坐标 (x,y[,z])',coord[i-1,:])
print('样本点的数量',nsp)
print('计算结果   ','{:13.4e}'.format(res))

gauss-legendre

def gauss_legendre(samp,wt):nsp=samp.shape[0]ndim=samp.shape[1]if ndim==1:if nsp==1:samp[0,0]=  0.0 wt[0]  =    2.0 elif nsp==2:samp[0,0]= -0.57735026918962576449 samp[1,0]=  0.57735026918962576449 wt[0]=      1.0 wt[1]=      1.0 elif nsp==3:samp[0,0]= -0.77459666924148337704 samp[1,0]=  0.0 samp[2,0]=  0.77459666924148337704 wt[0]=      0.55555555555555555556 wt[1]=      0.88888888888888888889 wt[2]=      0.55555555555555555556 elif nsp==4:samp[0,0]= -0.86113631159405257524 samp[1,0]= -0.33998104358485626481 samp[2,0]=  0.33998104358485626481 samp[3,0]=  0.86113631159405257524 wt[0]=      0.34785484513745385737 wt[1]=      0.65214515486254614271 wt[2]=      0.65214515486254614271 wt[3]=      0.34785484513745385737 elif nsp==5:samp[0,0]= -0.90617984593866399282 samp[1,0]= -0.53846931010568309105 samp[2,0]=  0.0 samp[3,0]=  0.53846931010568309105 samp[4,0]=  0.90617984593866399282 wt[0]=      0.23692688505618908751 wt[1]=      0.47862867049936646804 wt[2]=      0.56888888888888888889 wt[3]=      0.47862867049936646804 wt[4]=      0.23692688505618908751    else:print('积分点数量错误')elif ndim==2:if nsp==1:samp[0,0]=  0.0  samp[0,1]=  0.0 wt[0]=      4.0 elif nsp==4:samp[0,0]= -0.57735026918962576449 samp[1,0]=  0.57735026918962576449 samp[2,0]= -0.57735026918962576449 samp[3,0]=  0.57735026918962576449 samp[0,1]= -0.57735026918962576449 samp[1,1]= -0.57735026918962576449 samp[2,1]=  0.57735026918962576449 samp[3,1]=  0.57735026918962576449 wt[0]=      1.0 wt[1]=      1.0 wt[2]=      1.0 wt[3]=      1.0 elif nsp==9:samp[0,0]= -0.77459666924148337704 samp[1,0]=  0.0 samp[2,0]=  0.77459666924148337704 samp[3,0]= -0.77459666924148337704 samp[4,0]=  0.0 samp[5,0]=  0.77459666924148337704 samp[6,0]= -0.77459666924148337704 samp[7,0]=  0.0 samp[8,0]=  0.77459666924148337704 samp[0,1]= -0.77459666924148337704 samp[1,1]= -0.77459666924148337704 samp[2,1]= -0.77459666924148337704 samp[3,1]=  0.0 samp[4,1]=  0.0 samp[5,1]=  0.0 samp[6,1]=  0.77459666924148337704 samp[7,1]=  0.77459666924148337704 samp[8,1]=  0.77459666924148337704 wt[0] =     0.30864197530864197531 wt[1] =     0.49382716049382716049 wt[2] =     0.30864197530864197531 wt[3] =     0.49382716049382716049 wt[4] =     0.79012345679012345679 wt[5] =     0.49382716049382716049 wt[6] =     0.30864197530864197531 wt[7] =     0.49382716049382716049 wt[8] =     0.30864197530864197531 elif nsp==16:samp[0,0] =-0.86113631159405257524 samp[1,0] =-0.33998104358485626481 samp[2,0] = 0.33998104358485626481 samp[3,0] = 0.86113631159405257524 samp[4,0] =-0.86113631159405257524 samp[5,0] =-0.33998104358485626481 samp[6,0] = 0.33998104358485626481 samp[7,0] = 0.86113631159405257524 samp[8,0] =-0.86113631159405257524 samp[9,0]=-0.33998104358485626481 samp[10,0]= 0.33998104358485626481 samp[11,0]= 0.86113631159405257524 samp[12,0]=-0.86113631159405257524 samp[13,0]=-0.33998104358485626481 samp[14,0]= 0.33998104358485626481 samp[15,0]= 0.86113631159405257524 samp[0,1] =-0.86113631159405257524 samp[1,1] =-0.86113631159405257524 samp[2,1] =-0.86113631159405257524 samp[3,1] =-0.86113631159405257524 samp[4,1] =-0.33998104358485626481 samp[5,1] =-0.33998104358485626481 samp[6,1] =-0.33998104358485626481 samp[7,1] =-0.33998104358485626481 samp[8,1] = 0.33998104358485626481 samp[9,1]= 0.33998104358485626481 samp[10,1]= 0.33998104358485626481 samp[11,1]= 0.33998104358485626481 samp[12,1]= 0.86113631159405257524 samp[13,1]= 0.86113631159405257524 samp[14,1]= 0.86113631159405257524 samp[15,1]= 0.86113631159405257524 wt[0]=      0.12100299328560200551 wt[1]=      0.22685185185185185185 wt[2]=      0.22685185185185185185 wt[3]=      0.12100299328560200551 wt[4]=      0.22685185185185185185 wt[5]=      0.42529330301069429082 wt[6]=      0.42529330301069429082 wt[7]=      0.22685185185185185185 wt[8]=      0.22685185185185185185 wt[9]=     0.42529330301069429082 wt[10]=     0.42529330301069429082 wt[11]=     0.22685185185185185185 wt[12]=     0.12100299328560200551 wt[13]=     0.22685185185185185185 wt[14]=     0.22685185185185185185 wt[15]=     0.12100299328560200551 elif nsp==25:samp[0,0] =-0.90617984593866399282 samp[1,0] =-0.53846931010568309105 samp[2,0] = 0.0 samp[3,0] = 0.53846931010568309105 samp[4,0] = 0.90617984593866399282 samp[5,0] =-0.90617984593866399282 samp[6,0] =-0.53846931010568309105 samp[7,0] = 0.0 samp[8,0] = 0.53846931010568309105 samp[9,0]= 0.90617984593866399282 samp[10,0]=-0.90617984593866399282 samp[11,0]=-0.53846931010568309105 samp[12,0]= 0.0 samp[13,0]= 0.53846931010568309105 samp[14,0]= 0.90617984593866399282 samp[15,0]=-0.90617984593866399282 samp[16,0]=-0.53846931010568309105 samp[17,0]= 0.0 samp[18,0]= 0.53846931010568309105 samp[19,0]= 0.90617984593866399282 samp[20,0]=-0.90617984593866399282 samp[21,0]=-0.53846931010568309105 samp[22,0]= 0.0 samp[23,0]= 0.53846931010568309105 samp[24,0]= 0.90617984593866399282 samp[0,1] =-0.90617984593866399282 samp[1,1] =-0.90617984593866399282 samp[2,1] =-0.90617984593866399282 samp[3,1] =-0.90617984593866399282 samp[4,1] =-0.90617984593866399282 samp[5,1] =-0.53846931010568309105 samp[6,1] =-0.53846931010568309105 samp[7,1] =-0.53846931010568309105 samp[8,1] =-0.53846931010568309105 samp[9,1]=-0.53846931010568309105 samp[10,1]= 0.0 samp[11,1]= 0.0 samp[12,1]= 0.0 samp[13,1]= 0.0 samp[14,1]= 0.0 samp[15,1]= 0.53846931010568309105 samp[16,1]= 0.53846931010568309105 samp[17,1]= 0.53846931010568309105 samp[18,1]= 0.53846931010568309105 samp[19,1]= 0.53846931010568309105 samp[20,1]= 0.90617984593866399282 samp[21,1]= 0.90617984593866399282 samp[22,1]= 0.90617984593866399282 samp[23,1]= 0.90617984593866399282 samp[24,1]= 0.90617984593866399282 wt[0] =     0.05613434886242863595 wt[1] =     0.1134 wt[2] =     0.13478507238752090312 wt[3] =     0.1134 wt[4] =     0.05613434886242863595 wt[5] =     0.1134 wt[6] =     0.22908540422399111713 wt[7] =     0.27228653255075070182 wt[8] =     0.22908540422399111713 wt[9]=     0.1134 wt[10]=     0.13478507238752090305  wt[11]=     0.27228653255075070171  wt[12]=     0.32363456790123456757  wt[13]=     0.27228653255075070171  wt[14]=     0.13478507238752090305  wt[15]=     0.1134 wt[16]=     0.22908540422399111713 wt[17]=     0.27228653255075070182 wt[18]=     0.22908540422399111713 wt[19]=     0.1134 wt[20]=     0.05613434886242863595 wt[21]=     0.1134 wt[22]=     0.13478507238752090312 wt[23]=     0.1134 wt[24]=     0.05613434886242863595 else:print('积分点数量错误')elif ndim==3:if nsp==1:samp[0,0]=  0.0  samp[0,1]=  0.0 samp[0,2]=  0.0 wt[0]=      8.0 elif nsp==8:samp[0,0]=-0.57735026918962576449 samp[1,0]= 0.57735026918962576449 samp[2,0]=-0.57735026918962576449 samp[3,0]= 0.57735026918962576449 samp[4,0]=-0.57735026918962576449 samp[5,0]= 0.57735026918962576449 samp[6,0]=-0.57735026918962576449 samp[7,0]= 0.57735026918962576449 samp[0,1]=-0.57735026918962576449 samp[1,1]=-0.57735026918962576449 samp[2,1]=-0.57735026918962576449 samp[3,1]=-0.57735026918962576449 samp[4,1]= 0.57735026918962576449 samp[5,1]= 0.57735026918962576449 samp[6,1]= 0.57735026918962576449 samp[7,1]= 0.57735026918962576449 samp[0,2]=-0.57735026918962576449 samp[1,2]=-0.57735026918962576449 samp[2,2]= 0.57735026918962576449 samp[3,2]= 0.57735026918962576449 samp[4,2]=-0.57735026918962576449 samp[5,2]=-0.57735026918962576449 samp[6,2]= 0.57735026918962576449 samp[7,2]= 0.57735026918962576449 wt[0]=     1.0                                          wt[1]=     1.0                                          wt[2]=     1.0                                          wt[3]=     1.0                                          wt[4]=     1.0                                          wt[5]=     1.0                                          wt[6]=     1.0                                          wt[7]=     1.0                                          elif nsp==27:samp[0,0]= -0.77459666924148337704 samp[1,0]=  0.0 samp[2,0]=  0.77459666924148337704 samp[3,0]= -0.77459666924148337704 samp[4,0]=  0.0 samp[5,0]=  0.77459666924148337704 samp[6,0]= -0.77459666924148337704 samp[7,0]=  0.0 samp[8,0]=  0.77459666924148337704 samp[9,0]=-0.77459666924148337704 samp[10,0]= 0.0 samp[11,0]= 0.77459666924148337704 samp[12,0]=-0.77459666924148337704 samp[13,0]= 0.0 samp[14,0]= 0.77459666924148337704 samp[15,0]=-0.77459666924148337704 samp[16,0]= 0.0 samp[17,0]= 0.77459666924148337704 samp[18,0]=-0.77459666924148337704 samp[19,0]= 0.0 samp[20,0]= 0.77459666924148337704 samp[21,0]=-0.77459666924148337704 samp[22,0]= 0.0 samp[23,0]= 0.77459666924148337704 samp[24,0]=-0.77459666924148337704 samp[25,0]= 0.0 samp[26,0]= 0.77459666924148337704 samp[0,1]= -0.77459666924148337704 samp[1,1]= -0.77459666924148337704 samp[2,1]= -0.77459666924148337704 samp[3,1]= -0.77459666924148337704 samp[4,1]= -0.77459666924148337704 samp[5,1]= -0.77459666924148337704 samp[6,1]= -0.77459666924148337704 samp[7,1]= -0.77459666924148337704 samp[8,1]= -0.77459666924148337704 samp[9,1]= 0.0 samp[10,1]= 0.0 samp[11,1]= 0.0 samp[12,1]= 0.0 samp[13,1]= 0.0 samp[14,1]= 0.0 samp[15,1]= 0.0 samp[16,1]= 0.0 samp[17,1]= 0.0 samp[18,1]= 0.77459666924148337704 samp[19,1]= 0.77459666924148337704 samp[20,1]= 0.77459666924148337704 samp[21,1]= 0.77459666924148337704 samp[22,1]= 0.77459666924148337704 samp[23,1]= 0.77459666924148337704 samp[24,1]= 0.77459666924148337704 samp[25,1]= 0.77459666924148337704 samp[26,1]= 0.77459666924148337704 samp[0,2]= -0.77459666924148337704 samp[1,2]= -0.77459666924148337704 samp[2,2]= -0.77459666924148337704 samp[3,2]=  0.0 samp[4,2]=  0.0 samp[5,2]=  0.0 samp[6,2]=  0.77459666924148337704 samp[7,2]=  0.77459666924148337704 samp[8,2]=  0.77459666924148337704 samp[9,2]=-0.77459666924148337704 samp[10,2]=-0.77459666924148337704 samp[11,2]=-0.77459666924148337704 samp[12,2]= 0.0 samp[13,2]= 0.0 samp[14,2]= 0.0 samp[15,2]= 0.77459666924148337704 samp[16,2]= 0.77459666924148337704 samp[17,2]= 0.77459666924148337704 samp[18,2]=-0.77459666924148337704 samp[19,2]=-0.77459666924148337704 samp[20,2]=-0.77459666924148337704 samp[21,2]= 0.0 samp[22,2]= 0.0 samp[23,2]= 0.0 samp[24,2]= 0.77459666924148337704 samp[25,2]= 0.77459666924148337704 samp[26,2]= 0.77459666924148337704 wt[0]     = 0.17146776406035665295 wt[1]     = 0.27434842249657064472 wt[2]     = 0.17146776406035665295 wt[3]     = 0.27434842249657064472 wt[4]     = 0.43895747599451303155 wt[5]     = 0.27434842249657064472 wt[6]     = 0.17146776406035665295 wt[7]     = 0.27434842249657064472 wt[8]     = 0.17146776406035665295 wt[9]    = 0.27434842249657064472 wt[10]    = 0.43895747599451303155 wt[11]    = 0.27434842249657064472 wt[12]    = 0.43895747599451303155 wt[13]    = 0.70233196159122085048 wt[14]    = 0.43895747599451303155 wt[15]    = 0.27434842249657064472 wt[16]    = 0.43895747599451303155 wt[17]    = 0.27434842249657064472 wt[18]    = 0.17146776406035665295 wt[19]    = 0.27434842249657064472 wt[20]    = 0.17146776406035665295 wt[21]    = 0.27434842249657064472 wt[22]    = 0.43895747599451303155 wt[23]    = 0.27434842249657064472 wt[24]    = 0.17146776406035665295 wt[25]    = 0.27434842249657064472 wt[26]    = 0.17146776406035665295 else:print('积分点数量错误')else:print('维度不合理')

fun_der

def fun_der(fun,der,points,i):
#计算局部坐标形状函数的偏导ndim=der.shape[0]if ndim==1:xi= points[i-1,0]der[0,0]=-0.5; der[0,1]= 0.5fun[:]=(0.5*(1.0-xi),0.5*(1.0+xi))elif ndim==2:      xi= points[i-1,0]; eta=points[i-1,1]etam=0.25*(1.0-eta); etap=0.25*(1.0+eta)xim= 0.25*(1.0-xi); xip= 0.25*(1.0+xi)der[0,0]=-etam; der[0,1]=-etap; der[0,2]=etap; der[0,3]=etamder[1,0]=-xim; der[1,1]=xim; der[1,2]=xip; der[1,3]=-xipfun[:]=(4.0*xim*etam,4.0*xim*etap,4.0*xip*etap,4.0*xip*etam)elif ndim==3:   xi  =points[i-1,0]; eta =points[i-1,1]; zeta=points[i-1,2]etam=1.0-eta; xim=1.0-xi; zetam=1.0-zetaetap=eta+1.0; xip=xi+1.0; zetap=zeta+1.0der[0,0]=-0.125*etam*zetam; der[0,1]=-0.125*etam*zetapder[0,2]= 0.125*etam*zetap; der[0,3]= 0.125*etam*zetamder[0,4]=-0.125*etap*zetam; der[0,5]=-0.125*etap*zetapder[0,6]= 0.125*etap*zetap; der[0,7]= 0.125*etap*zetamder[1,0]=-0.125*xim*zetam; der[1,1]=-0.125*xim*zetapder[1,2]=-0.125*xip*zetap; der[1,3]=-0.125*xip*zetamder[1,4]= 0.125*xim*zetam; der[1,5]= 0.125*xim*zetapder[1,6]= 0.125*xip*zetap; der[1,7]= 0.125*xip*zetamder[2,0]=-0.125*xim*etam; der[2,1]= 0.125*xim*etamder[2,2]= 0.125*xip*etam; der[2,3]=-0.125*xip*etamder[2,4]=-0.125*xim*etap; der[2,5]= 0.125*xim*etapder[2,6]= 0.125*xip*etap; der[2,7]=-0.125*xip*etapfun=(0.125*xim*etam*zetam,0.125*xim*etam*zetap,                       \0.125*xip*etam*zetap,0.125*xip*etam*zetam,                     \0.125*xim*etap*zetam,0.125*xim*etap*zetap,                     \0.125*xip*etap*zetap,0.125*xip*etap*zetam)else:print('维数错误')