基于Gauss积分法计算数值积分2

根据前文的介绍，对于一重积分，我们运用Gauss-Legendre积分公式可以得到：

而对于如下所示的二重积分：

式中，积分区域Ω为：

我们可以先对第一个坐标η积分，然后再对第二个坐标ξ积分，这样可以得到针对二重积分的Gauss-Legendre积分公式：

上式中，在两个坐标方向上高斯点数量不一定需要相等，即n不一定等于m，不过一般为了便于计算，都会采用相同的高斯点。例如，对于两点Gauss-Legendre积分公式，由前文可知其高斯点位置和权重分别为：

如果我们在两个方向上均采用2点Gauss-Legendre积分公式，则在积分区域内的高斯点布置如图1所示。

图1 2点Gauss-Legendre公式的积分点布置

因此，Gauss-Legendre积分公式中的双重求和实际上是遍历积分区域中的每一个高斯点进行求和。

下面采用Gauss-Legendre积分公式计算如下二重积分：

式中，积分区域为：

（1）解析计算

为了便于比较Gauss-Legendre积分公式在计算二重积分时的计算精度，下面首先计算上述积分的精确解。

由于积分区域为矩形区域，因此很容易给出被积函数在积分区域的分布，如图2所示。根据二重积分的定义，该积分实际上代表了在积分区域内被积函数与XOY平面所围成的区域的体积。

图2 被积函数在积分区域的分布

（2）数值计算

由于被积函数中x和y的最高次均为2，已知Gauss型积分公式的代数精度为2n-1，因此预计采用2点Gauss-Legendre积分公式进行计算应满足精度要求，即需要在积分区域内布置4个高斯点。注意到积分区域为Ω={(x,y)|-1≤x≤1,-1≤y≤1}，因此高斯点的坐标为：

并且每个高斯点的权系数A_i (i=1,2)均为1。

由2点Gauss-Legendre公式可得二重积分的计算式为：

因此有：

可见由直接积分和Gauss型积分给出计算结果完全吻合。

需要指出的是，前面示例的二重积分的积分区域均为η=±1和ξ=±1构成的矩形区域，这是因为Legendre多项式是在区间[-1,1]带权ρ(x)=1的正交多项式。对于积分区间不在[-1,1]的一重积分，我们可以通过变换区间的方式使其适用于Gauss-Legendre积分公式；那么对于任意积分区域的二重积分，是否同样可以通过变换区间的方式使其适用于Gauss-Legendre积分公式。为解决该问题，我们首先对所需的一些预备知识做简要介绍。

下面我们首先给出一个简单的示例：

式中，积分区域Ω可表示为：

根据高中知识可以很容易判断积分区域为一个长轴为2a、短轴为2b的椭圆，积分I实际上为椭圆的面积。尽管很容易通过椭圆的面积计算公式给出积分值，但我们还是希望通过二重积分的方式进行求解。

利用二重积分求解椭圆的面积有很多方法，这里我们给出一种非常巧妙的解法。由于椭圆实际上可以看作是一个被压扁的圆，即在圆的X轴和Y轴上对圆的半径进行不同比例的缩放，因此我们定义变换：

则原来的积分区域Ω变换为一个圆形的积分区域：

注意到

因此原积分可表示为：

由于变换后积分区域为一个半径为1的圆，因此积分值为πab。

可以看到，在本例中，我们通过对积分区域进行变化使得该积分更易求解，但需要注意的是在变换前后，两个积分区域的面积并不相同，因为微元的面积满足：

这里的ab可以看作是两个微元面积的比例缩放因子，虽然在本例中这是一个常量，但大多时候可能是一个变量。

为了更清楚地了解变量变换前后微元面积之间的关系，我们取XOY平面中左下角坐标为(x, y)的一个微元进行分析，由于微元非常小，因此假定微元为规则的矩形，如图3所示，则微元其余三个角点的坐标均可通过左下角点的坐标表示。

图3 位于XOY平面的微元

因此，微元的面积为：

下面，我们引入线性变换：

如果将该变换应用于上面的微元，则变换之后在UOV平面的微元如图4所示，注意到由于采用的是线性变换，因此变换之后微元的四条边仍为直线，只是相临边的角度可能会发生改变。

图4 变换后位于UOV平面的微元

图4中，变换之后微元四个角点的坐标可分别表示为：

因此有：

由多元微分学的知识可知：

因此可以将向量重新表示为：

因此微元的面积可以表示为：

如果设

则可以得到两个微元面积之间的关系为：

这里的J₂的绝对值相当于就是两个微元面积之间的缩放因子，J₂即为2阶Jacobian行列式。利用上述关系，如果某个二重积分在当前积分域难以积分时，我们可以通过构造合适的变换将其变换到另一个较容易积分的积分域中：

同理，如果定义逆变换：

则同样可以得到：

注意此时的2阶Jacobian行列式为：

2.2.1 二重积分算例2

下面以一个简单的例子来说明上述变换的正确性，已知一个二重积分为：

很容易知道该二重积分的值代表一个半径为1的圆的面积。我们引入变换将其变换到极坐标系中：

其2阶Jacobian行列式应表示为：

因此变换前后微元面积之间的关系为：

因此原积分可表示为：

对于更为一般的情况，如果我们针对积分变量y₁,…,y_n定义如下变换：

若它们对每个自变量x_j都存在偏导数∂f_i/∂x_j（i=1,2,…,n; j=1,2,Λ,n），则n阶Jacobian行列式为：

因此，对于如下所示的三重积分：

如果我们定义变换：

则可将上述积分变换为：

式中，3阶Jacobian行列式为：

当然，上述变换也可以参照2.2节中的方法，通过计算变换前后微元体积之间的关系来获得，即证明：

计算过程微元的体积可以通过混合积计算得到：

式中：a、b和c分别为微元三条相临边构成的矢量。这里不再给出具体的证明过程。

2.3.1 三重积分算例1

下面以一个简单的例子来说明上述变换的正确性，已知一个三重积分为：

很容易知道该二重积分的值代表一个半径为1的球的体积。我们引入变换将其变换到球坐标系中：

则3阶Jacobian行列式J₃为：

即

因此可以得到：

从节2中的分析可以看到，对于一个难以进行数值积分的任意积分区域，如果我们能够将其变换到Gauss-Legendre积分公式的积分区域（Ω={(x, y)|-1≤x≤1, -1≤y≤1}），并且利用Jacobian行列式建立变换前后的微元面积之间的关系，就可以使用Gauss-Legendre积分公式来数值计算该积分。因此，下面首先讨论如何建立该变换关系。

假定变换之后的积分区域位于自然坐标系η-ξ中，按照二重积分的Gauss-Legendre积分公式对积分区域的要求，将坐标原点取在积分区域的中心，并且将积分区域的四个角点和边界限制在+1或-1上，如图5所示。

图5 位于自然坐标系中的积分区域

假定变换之前位于原来的积分区域的全局坐标x和y与变换之后的位于Gauss-Legendre积分公式的积分区域的自然坐标η和ξ之间满足如下变换关系：

图6给出了经过上述变换后映射到全局坐标系x-y上的四边形积分区域，四边形的大小和形状可以通过四个角点处的坐标来确定。

图6 映射到全局坐标系的积分区域

根据(x_i, u_i)和(η_i, ξ_i)（i=1,..,4）之间的关系，可得：

解得：

当然，另一种表示形式为：

将其表示为矩阵形式为：

式中：

熟悉有限元方法的朋友应该可以看出，实际上N_i (i=1,…,4)为4节点线性平面四边形单元的形函数。

因此，其二阶Jacobian行列式为：

即有：

因此可以将二重积分的原积分区域Ω变换至Gauss-Legendre积分公式所要求的积分区域，这样即可应用Gauss-Legendre积分公式来求解二重积分：

3.1.1 二重积分算例3

下面以一个简单的二重积分来介绍如何运用上述方法来完成积分的数值计算：

这里我们将积分区域Ω设为：

（1）解析计算

这里我们还是首先采用解析方法计算上面的二重积分。注意到积分区域为一个中心在原点的正方形，并且被积函数关于变量x和y均为偶函数，因此可以充分利用对称性将其化简为：

如果沿X轴对其积分，则有：

（2）数值计算

已知原积分区域四个角点的坐标分别为：

因此如果引入变换

可以解得：

因此有

不难验证变换之后的积分区域为：

因此有：

实际上对照1.1.1节中的二重积分算例1已经可以看出，本算例的积分值为算例1中积分值的1/4，即积分值I=2/3，因此这里不再给出运用Gauss-Legendre积分公式求解上述积分的完整过程。可见解析积分和数值积分的计算结果完全吻合。

在实际运用中，积分区域很可能并不是规则的，如积分区域可能为三角形、多边形或任意形状，此时无法直接运用上述方法进行积分区域的变换。因此，本文参照有限元法的思想，首先积分区域离散为若干个四边形子积分区域Ω_i，然后针对每个子积分区域进行积分区域的变换，并通过Gauss-Legendre积分公式进行数值积分，最后将每个子积分区域的积分值进行累加，获得总的积分值：

3.2.1 二重积分算例4

下面以一个简单的示例来验证上述数值计算方法的有效性，我们构造一个二重积分，并使得积分区域为圆形：

（1）解析计算

令

则有：

（2）数值计算

为构造子积分区域，我们首先在有限元前处理软件HyperMesh中利用四边形单元划分积分区域，并将每个单元作为一个子积分区域。例如，图7给出了当子积分区域数量为384时的积分区域划分情况。然后，将该网格模型输出为Abaqus求解文件（.inp）以获得单元和节点信息。本文基于MATLAB编制了二重积分的数值计算程序（参见附录A），通过读取单元信息和节点信息，以获取每个子积分区域的角点坐标，然后基于3.1节中的方法完成子积分的数值积分，最终实现任意积分区域的二重积分的数值计算。

图7 积分区域划分示意图（子积分区域数量为384）

为研究子积分区域数量与计算精度之间的关系，本文分别选取子积分区域数量为4、32、384和1536，积分点数量为4、9、16和25，计算得到子积分域数量、积分点数量和计算误差之间的关系如图8所示。

图8 子积分域数量、高斯点与积分值之间的关系

从图8中可以看到，随着子积分域数量的增加，数值积分的精度有了明显提升（虽然采用了对数坐标，但提升趋势仍然非常明显），而高斯点数量对积分精度的提升影响不大。因此，正如在前文中所提到的，在计算二重积分时同样应尽量采用较少数量的积分点以避免截断误差，转而以增加子积分区域数量的方法来提高计算精度。

从本质上来说，本文提到的划分子积分区域的方法属于复合积分公式的计算方法，虽然采用梯形公式或其他积分公式同样可以实现任意积分区域的二重积分的数值积分，但由于Gauss型积分的代数精度高于梯形公式等积分公式。因此，在子积分区域数量相同的情况下，采用Gauss型积分公式计算二重积分可以获得更高的计算精度。

上面讨论的情况都是被积函数的形式已知的情况，而在某些特殊情况下，有可能被积函数的具体形式未知，而仅有子积分区域角点处的函数值（这种情况多发生于对有限元模型进行后处理时）。此时，如果我们需要计算子积分区域高斯点位置处的函数值，则只能通过插值的方式来获取。由于子积分区域角点处的函数值是已知的，因此我们可以直接使用定义子积分区域形状的形函数来进行插值，即有：

式中：f(η_j, ξ_k)为自然坐标为η_j和ξ_k时的函数值，N_j为子积分区域第j个角点的形函数。由于这里定义自然坐标和全局坐标的变换关系与计算积分点处的函数值采用了同一种插值方式（即使用了同一套形函数），因此常被称为平面四边形单元的等参公式描述。上述公式只需对附录A中的程序稍加改进即可嵌入，因此这里不再给出具体的数值实现算法。

从前面的数值计算方法可以看到，对于具有任意积分区域的三重积分，只要采用六面体单元将积分区域划分为若干个子积分区域，即可采用相应的形函数将其变换到基于Gauss-Legendre积分公式计算三重积分所需的积分域中。对于8节点线性六面体单元，参照前面的方法，可以假定变换之前位于原积分区域的全局坐标x、y和z与变换之后位于Gauss-Legendre积分公式的积分区域的自然坐标η、ξ和γ之间满足如下变换关系：

上式同样可以通过形函数表示为：

式中的形函数可表示为：

其中η_i、ζ_i和γ_i（i=1,…,8）分别取值为+1或者-1。

例如，假定六面体单元的第一个节点的坐标为(x₁,y₁,z₁)，该节点映射到自然坐标系下的坐标为(1,-1,1)，则该节点对应的形函数为：

其他节点处的形函数可以此类推。

根据上述变换，只要计算得到对应的3阶Jacobian行列式，即建立变换前后微元体积之间的关系：

这样即可应用Gauss-Legendre积分公式来求解三重积分：

式中3阶Jacobian行列式的具体形式这里不再给出，有需要的读者可以参照3.1节的推导过程自行推导。

本文基于Gauss-Legendre积分公式，提供了一种针对任意积分区域的二重积分数值积分方法。该方法通过将积分区域划分为若干个四边形子积分区域，然后在每个子积分区域内运用Gauss-Legendre积分公式，最后累加每个子积分区域的积分值以获得总的积分值。

需要主要的是，如果积分区域存在较多曲边，有时候仅采用四边形单元进行划分不一定能取得较好的效果，因此后续将继续发展针对三角形区域的Gauss-Legendre数值积分方法，通过采用四边形和三角形混合单元划分积分区域，以实现对积分区域更好的逼近效果。

[1] 邹泽民, 刘志伟. 雅可比行列式在积分坐标变换中的应用. 广西梧州师范高等专科学校学报. 2005, 21(3).

[2] Daryl L. Logan. A first course in the finite element method, SI edition, fifth edition. Cengae Learning Asia Pet. Ltd.

附录A：二重积分MATLAB计算程序

%程序采用Gauss-Legendre积分公式计算二重积分clear,clc%********************************参数输入**********************************Gauss_num=5;                               %高斯点数量%**************************************************************************
%***************************被积函数及积分域定义***************************%被积函数定义Integrand_func=@(x,y) exp(x^2+y^2);%导入积分域定义信息Elem_data=importdata("Elem_id.txt");       %子积分域定义Node_data=importdata("Node_id.txt");       %子积分域节点定义%积分区域信息统计Elem_num=size(Elem_data,1);                %子积分域总数Node_num=size(Node_data,1);                %节点总数%**************************************************************************
%***************************求解Gauss点位置及权重**************************%确定Gauss点位置syms xRoots=vpasolve(legendreP(Gauss_num,x)==0);Gauss_Pos=double(Roots);clear x%根据递推公式计算Legendre多项式导数Legendre_DPoly=Gauss_num*(Gauss_Pos.*legendreP(Gauss_num,Gauss_Pos)-...  legendreP(Gauss_num-1,Gauss_Pos))./(Gauss_Pos.^2-1);%计算Gauss点权重Gauss_Weight=2./((1-Gauss_Pos.^2).*(Legendre_DPoly.^2));%**************************************************************************
%*********************************数值积分*********************************I_elem=zeros(Elem_num,1);                     %存放每个子积分域的积分值for i=1:Elem_num                            %子积分域的循环  %查找当前积分域的节点坐标  Coord_X=zeros(4,1);  Coord_Y=zeros(4,1);  for j=1:4    index=find(Node_data(:,1)==Elem_data(i,j+1));    Coord_X(j)=Node_data(index,2);            %节点坐标X    Coord_Y(j)=Node_data(index,3);            %节点坐标Y  end  %计算变量变换公式中的系数a和b  Coeff_a=zeros(4,1);  Coeff_b=zeros(4,1);  Coeff_a(1)=0.25*(Coord_X(1)+Coord_X(2)+Coord_X(3)+Coord_X(4));  Coeff_b(1)=0.25*(Coord_Y(1)+Coord_Y(2)+Coord_Y(3)+Coord_Y(4));  Coeff_a(2)=0.25*(Coord_X(1)+Coord_X(2)-Coord_X(3)-Coord_X(4));  Coeff_b(2)=0.25*(Coord_Y(1)+Coord_Y(2)-Coord_Y(3)-Coord_Y(4));  Coeff_a(3)=0.25*(Coord_X(1)-Coord_X(2)-Coord_X(3)+Coord_X(4));  Coeff_b(3)=0.25*(Coord_Y(1)-Coord_Y(2)-Coord_Y(3)+Coord_Y(4));  Coeff_a(4)=0.25*(Coord_X(1)-Coord_X(2)+Coord_X(3)-Coord_X(4));  Coeff_b(4)=0.25*(Coord_Y(1)-Coord_Y(2)+Coord_Y(3)-Coord_Y(4));  %基于Gauss-Legendre积分公式进行数值积分  for j=1:Gauss_num    for k=1:Gauss_num      %计算当前Gauss点的2阶Jacobian行列式      Jacobian_Matrix=...        [Coeff_a(2)+Coeff_a(4)*Gauss_Pos(k) Coeff_b(2)+Coeff_b(4)*Gauss_Pos(k);        Coeff_a(3)+Coeff_a(4)*Gauss_Pos(j) Coeff_b(3)+Coeff_b(4)*Gauss_Pos(j)];      Jacobian_det=det(Jacobian_Matrix);      %计算当前Gauss点对应的全局坐标x和y      Coord_X_Gauss=Coeff_a(1)+Coeff_a(2)*Gauss_Pos(j)+...        Coeff_a(3)*Gauss_Pos(k)+Coeff_a(4)*Gauss_Pos(j)*Gauss_Pos(k);      Coord_Y_Gauss=Coeff_b(1)+Coeff_b(2)*Gauss_Pos(j)+...        Coeff_b(3)*Gauss_Pos(k)+Coeff_b(4)*Gauss_Pos(j)*Gauss_Pos(k);      %调用Gauss-Legendre积分公式      I_elem(i)=I_elem(i)+Gauss_Weight(j)*Gauss_Weight(k)*...        Integrand_func(Coord_X_Gauss,Coord_Y_Gauss)*abs(Jacobian_det);    end  endend%累加所有子积分域的积分值I=sum(I_elem);fprintf("积分值为%f\n",I);%**************************************************************************

附录B：积分域节点数据输入格式(Node_id.txt)

         1,  0.7071067811865,  -0.707106781186         2,  -0.707106781186,  -0.707106781186         3,  -1.22464680E-16,  -1.0                    4,  -1.0           ,  -4.44089210E-16         5,  -0.707106781186,  0.7071067811865         6,  0.7071067811865,  0.7071067811865         7,  1.0            ,  -4.44089210E-16         8,  0.0            ,  1.0                     9,  0.0            ,  -2.22044605E-16

附录C：积分域单元数据输入格式(Elem_id.txt)

         1,         3,         1,         7,         9         2,         3,         9,         4,         2         3,         8,         5,         4,         9         4,         8,         9,         7,         6