基于Gauss积分法计算数值积分
在有限元分析求解单元的刚度矩阵时，通常需要计算在整个单元域的积分值，这些积分往往没有解析解，只能通过数值方法进行求解。对于等参单元，由于单元的刚度矩阵采用等参坐标进行描述，因此采用Gauss积分法计算刚度矩阵时具有非常明显的优势。本文首先介绍了一些常见的数值积分方法，然后重点介绍了如何基于任意积分点数量的Gauss-Legendre公式完成数值积分，并给出了相应的计算程序。
当某些积分无法用解析方法进行求解时，就需要借助数值方法进行离散化后计算积分的近似值，我们可以采用简单可积的函数来近似替代被积函数，如多项式函数。以一元函数为例，我们可以利用多项式的积分来逼近所求的积分：
式中：Pn(x)为n阶多项式，可表示为：
为使得多项式Pn(x)能尽可能替代被积函数f(x)，根据多项式插值的原理，应使得在[a,b]的积分区间内多项式在n+1个点（也被称为采样点或积分点）处的取值与被积函数的取值相同。在求得多项式的系数后，即可将被积函数的积分转化为对多项式的积分，由于多项式简单可积，因此很容易即可给出函数的积分值。
-1-
Newton-Cotes公式
Newton-Cotes公式的特点在于在积分区间内均匀地取点来构造多项式函数，根据取点个数的不同，可分为以下几类：
该方法通过选取一个积分点来构造多项式P0(x)，通常这个点选为积分区间的中点，如图1所示，图中阴影区域的面积即为多项式的积分值，因此被积函数的积分近似值可表示为：
显然这种方法计算得到的积分值并不精确，只有当被积函数f(x)同为零阶多项式（即一条平行于x轴的直线）时积分才能取得精确值。
图1 矩形公式积分示意图
梯形公式（n=2）
该方法通过选取两个点来构造一阶多项式P1(x)，并将两个积分点分别选为积分区间的下限和上限值，如图2所示。
图2 梯形公式积分示意图
图中阴影区域的面积即为一阶多项式P1(x)的精确积分值，因此根据梯形面积的计算公式有：
同样地，只有当被积函数同为一阶多项式时，上述积分才能取得精确值。
Simpson公式（n=3）
用直线代替被积函数的精度可能并不好，因此继续增加积分点的数量，即采用抛物线来替代被积函数。由于用于近似被积函数的近似函数为抛物线，因此可设为：
因此原积分可近似表示为：
与前面类似，将积分区间[a,b]等分为两份，从而得到三个积分点用于构造二次插值多项式，因此在三个积分点处有：
根据矩形公式和梯形公式，可以推测采用抛物线作为近似函数的积分公式应具备如下所示的形式：
联立两式可以得到：
可以解得：
因此积分公式可表示为：
以上积分公式也被称为Simpson公式，Simpson公式对于不高于三次的多项式积分都是精确成立的。为了比较这些积分公式的精度，一般会给出代数精度定义，即如果积分公式对任意不高于n阶的多项式积分均精确成立，而对高于n+1阶的多项式不能精确积分，则称该积分公式具有n阶代数精度。确定积分公式的代数精度需要对其进行误差分析，这里不再详细讨论。
Newton-Cotes公式
基于前面的方法，下面采用n次多项式Pn(x)来近似代替被积函数，即将积分区间均匀分割为n个子区间（共n+1个插值节点），并使得在每个采样点处多项式的函数值与被积函数的函数值相等，根据Lagrange插值法可以得到穿过这些点的n次Lagrange插值多项式为：
因此原积分可表示为：
式中：
式中：Cn(k)被称为Cotes系数，它与被积函数和积分的上下限均无关，只与插值点的数量有关。
这种通过在积分区间上等间距布置插值节点，采用n次插值多项式替代被积函数来近似求解积分的方法即为Newton-Cotes公式。当多项式阶次n为偶数时，求积公式有n+1次代数精度；而当n为奇数时，求积公式只有n次代数精度。通过提高多项式的阶次可以提高求积公式的代数精度，但过高阶次的积分公式会变得不稳定，因此高次积分最多用到7次。
如果需要提高积分精度，且不能随意更改插值多项式的阶次，则一种方法为采用复合积分公式，即将积分区域均匀分割为若干个子区间，在每个子区间上使用低次积分公式，这不在本文的讨论范围；而另一种方法为去掉插值节点等间距的限制，尝试通过移动插值节点的位置来提高积分公式的代数精度，即Gauss型积分。
-2-
Gauss型积分
基于上述内容可以看到，对于插值型积分，可以统一表示为：
如果插值节点是等距的，即为Newton-Cotes公式，并且当插值节点个数n为奇数时，代数精度为是n次。在Newton-Cotes公式中，只要插值节点个数n确定下来，则插值节点xi的位置以及对应的权系数Ai都是确定。而Gauss型积分可以通过适当改变插值节点的位置，在有限的插值节点下进一步提高积分公式的代数精度。
Gauss型积分定义
Gauss型积分的定义为：
gn(x)是[a,b]上带权ρ(x)的n次正交多项式，xk(k=1,2,…,n)是gn(x)的零点，则积分公式：
有(2n-1)次代数精度，其中xk(k=1,…,n)称为高斯点，并且权系数Ak可表示为：
式中：lk(x)为Lagrange插值基函数。
可以看到，Gauss型积分在使用时并没有太多限制，只需要给定积分区间[a,b]，权函数ρ(x)，以及代数精度（假定为2n-1次），即可构造出相应的Gauss型积分。构造Gauss型积分可分为两步，第一步为构造n次的正交多项式，并求解其零点作为高斯点；第二步为确定对应的权系数。下面以在区间[-1,1]带权ρ(x)=1的正交多项式g0(x)、g1(x)以及g2(x)为例，介绍如何构造其Gauss型积分公式。
构造正交多项式
下面首先给出正交多项式的定义：
如果多项式
满足内积
则称gn(x)为[a,b]上带权ρ(x)的n次正交多项式。
因此，对于在区间[-1,1]带权ρ(x)=1的正交多项式g0(x)、g1(x)以及g2(x)，应满足：
根据正交多项式的定义有：
可以解得：
因此高斯点为：
对于更高次的正交多项式，可采用上述方法继续求解。
求解权系数
权系数Ak可通过Lagrange插值基函数进行求解，但此积分计算较为繁琐。需要注意到Gauss型积分具有2n-1次代数精度，故求积公式对任意不高于2n-1次的代数多项式积分应精确成立，因此当f(x)=1,x时，由于积分精确成立，故可以得到：
可以解得权系数为：
对于更高次的正交多项式，同样可以以此类推。
因此，在[-1,1]的积分区间内，上述Gauss型积分公式为：
上式即为最常见的两点Gauss积分公式。
Legendre多项式的应用
当积分计算时要求的代数精度较高时，用上述方法确定正交多项式需要解析计算积分，计算过程非常繁琐。下面介绍一个特殊的正交多项式，Legendre多项式，利用该正交多项式可以简化Gauss型积分的构造过程。
Legendre多项式是[-1,1]上带权ρ(x)=1的正交多项式，可表示为：
Legendre多项式具有奇偶性质，即当n为偶数时Pn(x)为偶函数，n为奇数时Pn(x)为奇函数。并且满足如下递推公式：
由于Legendre多项式为正交多项式，因此可以用来构造Gauss型积分，并且根据Gauss型积分的定义，Pn(x)=0的根xk（k=1,2,…,n）即为积分点，并且
具有（2n-1）次代数精度，其中
以上构造得到的Gauss型积分也被称为Gauss-Legendre积分公式。
需要注意的是，对于积分限在[a,b]并且a，b为有限数的定积分，通过将其区间变换到[-1,1]，通过可以采用Gauss-Legendre公式进行计算，即需要构造变换使得：
令x=α+βt，根据上式有：
因此有：
代入到原积分有：
式中：tk为Legendre多项式Pn(t)的根，并且
算例
下面采用Gauss-Legendre积分公式计算如下积分：
（1）解析计算
为了便于比较Gauss型积分的计算精度，下面首先计算上述积分的精确解。
令x=sin θ，则有：
因此有：
（2）数值计算
这里采用3点Gauss-Legendre积分公式计算上述积分，该公式具有5阶代数精度。高次的Legendre多项式不易直接计算，可以先计算第零次和第一次的Legendre多项式，然后利用递推公式计算高次的Legendre多项式。根据Legendre多项式的定义有：
由递推公式可以得到：
由P3(x)=0可得对应的高斯点为：
而由
代入xk（k=1,2,3）权系数为：
代入3点Gauss-Legendre积分公式
可得积分的近似值为：
而理论值与近似值的误差为：
可以看到理论值与近似值之间的误差较大，观察被积函数可以发现，在积分上下限处的函数值为无穷大，这说明3点Gauss-Legendre积分公式的代数精度可能并不满足要求。针对此种情况有两种解决方法，第一种方法为在保持当前代数精度不变的情况下，将积分区间分割为若干个间隔相等的小区间，在小区间上分别应用3点Gauss-Legendre积分公式，然后将小区间上的积分值进行累积获得被积函数的积分值；第二种方法为继续增加积分点，以提高积分公式的代数精度。第一种方法属于复合积分公式的内容，这里主要讨论第二种方法的实现。
多点Gauss-Legendre积分公式的实现
如果继续提高积分点数量，则需要构造更高次的Legendre多项式，并求解其零根。当积分点数量超过3个后，Legendre多项式非常复杂，往往难以通过解析方法确定其零根，只能通过数值方法进行求解，或者查询已经编制好的高斯积分表，但显然查表不具有通用性。因此，如何确定多点Gauss-Legendre积分公式的积分点位置以及权系数，是本节将要讨论的内容。
n阶Legendre多项式及其导数的计算
设高斯点数量为n，为了计算n阶Legendre多项式，可以首先根据Legendre多项式的定义得到1阶和2阶Legendre多项式，然后再利用递推公式得到更高次的Legendre多项式。需要注意的是，在实际数值计算中，按照Legendre多项式的定义来计算会涉及到大量的乘除法运算，同时由于数据字节长度的限制，对于基数较大阶乘的运算会导致数据的溢出，因此使用递推公式进行计算更为合适。根据定义有：
而递推公式为：
为了便于编程实现，令n=k+1，得到第n阶Legendre多项式Pn(x)的递推公式为：
除此之外，在计算权系数时，还需要用到第n阶Legendre多项式的导数P'n(x)。Legendre多项式的导数同样存在如下所示的递推关系：
可以看到第四个递推公式只需用到n-1阶和n阶Legendre多项式即可推出n阶Legendre多项式的导数，计算过程最为简单，因此本文采用该递推公式来计算第n阶Legendre多项式的导数：
n阶Legendre多项式零根的计算
在获得n阶Legendre多项式及其导数后，就可以开始计算Legendre多项式的零根，即高斯点的位置。为了提高计算效率，需要注意到，当n为偶数时，Legendre多项式为偶函数；当n为奇数时，Legendre多项式为奇函数，因此只需计算区间[0,1]内的零根即可，并且当n为奇数时必有一个零根为x=0。例如，图3给出了前5阶Legendre多项式的曲线。
图3 前5阶Legendre多项式曲线
此外，Legendre多项式也不存在重根及虚根的情况。这些限制条件极大地方便了Legendre多项式零根的计算。对于非线性方程求根的问题，常用的数值方法有二分法（Bisection Method）、不等点迭代法以及Newton-Raphson法等。
本文采用易于实现的二分法来求根，首先将积分区间[-1,1]划分为若干个等间距的子区间，为了保证可以找到Legendre多项式在积分区间内的所有零根，应使得每个子区间应尽可能的小，这里将子区间数量Ns取为高斯点数量n的50倍。然后遍历每一个子区间，假定第i个子区间的范围为[Bi, Bi+1]，如果该子区间满足：
则表明该区间必然存在零根，然后对该子区间使用二分法，即可确定零根。需要注意的是，二分法在存在零根的子区间内是必然收敛的，但收敛速度不快，如果对计算效率有要求，可以使用收敛速度更快的Newton-Raphson法。
在求得零根后，即可根据零根以及Legendre多项式在零根处的导数值，得到相应的权系数，然后根据Gauss-Legendre积分公式即可计算被积函数的积分值。
本文基于上述方法编制了函数数值积分的Fortran程序（参见附录A），支持输入任意数量的高斯点，计算相应的高斯点位置，权系数以及被积函数的积分值。例如，下面给出了由程序计算得到的5点Gauss-Legendre积分的高斯点位置以及权系数：
高斯点序号: 1
高斯点位置: -0.90617984593866341
高斯点权值: 0.23692688505619044
高斯点序号: 2
高斯点位置: -0.53846931010568233
高斯点权值: 0.47862867049936686
高斯点序号: 3
高斯点位置: 0.0000000000000000
高斯点权值: 0.56888888888888889
高斯点序号: 4
高斯点位置: 0.53846931010568233
高斯点权值: 0.47862867049936686
高斯点序号: 5
高斯点位置: 0.90617984593866341
高斯点权值: 0.23692688505619044算例
针对前面算例中的被积函数，本文再次采用更多积分点的Gauss-Legendre积分公式求解其数值积分，并研究了积分点数量与积分值之间的关系，结果如图4所示。
图4 示例函数积分值随高斯点数量的变化
从图4中可以看到，随着高斯点数量的增加，由数值方法确定的积分近似值逐渐趋近于理论值，当高斯点数量达到200个后，近似值与理论值基本吻合。但也可以看到，Gauss-Legendre公式并非高斯点越多越精确，高斯点变多后会有很多的截断误差。因此，最好的方法是使用复合积分公式，即将积分区间划分为若干个子区间，然后在每个子区间中使用少量的高斯点进行积分运算。
参考文献
[1] 金一庆,陈越,王冬梅.数值方法.机械工业出版社.
[2] 余海洋,方世跃.关于勒让德多项式递推公式的研究.四川理工学院学报(自然科学版),2008,21(2).
END
附录A
!!程序使用Gauss-Legendre积分公式求解积分问题MODULE Module_Gauss_Legendre_Weight!!模块用于求解n点Gauss_Legendre积分公式的权重系数IMPLICIT NONE::DP=8 !设置浮点型精度::LONG_INT=4 !设置整型精度::Num_Iter_Max=1000 !最大迭代次数=DP)::EPS=1.0D-15 !设置求零根时的迭代容差=DP),PARAMETER::ZR=0D0=DP),PARAMETER::ONE=1.0D0=DP),PARAMETER::TWO=2.0D0=DP),PARAMETER::THR=3.0D0=DP),PARAMETER::HLF=0.5D0CONTAINS=DP) FUNCTION Poly_Legendre(n,x)!!函数用于构造n阶Legendre多项式IMPLICIT NONE!!函数变量声明=LONG_INT)::n !Legendre多项式阶次(高斯点数)=LONG_INT)::I !循环索引变量=DP)::x !用于迭代计算的采样点=DP)::P(n) !n阶Legendre多项式在采样点x处的取值==1)THEN !阶次为1的特殊情况=x==2)THEN !阶次为2的特殊情况=HLF*(THR*x*x-ONE)ELSE !阶次大于2的一般情况!!构造1阶和2阶Legendre多项式=x !1阶Legendre多项式=HLF*(THR*x*x-ONE) !2阶Legendre多项式!!利用递推公式构造n阶Legendre多项式DO I=3,n=((2*I-1)*x*P(I-1)-(I-1)*P(I-2))/IEND DOEND IF!!生成指定阶次的Legendre多项式Poly_Legendre=P(n)END FUNCTION Poly_Legendre=DP) FUNCTION DPoly_Legendre(n,x)!!函数用于构造n阶Legendre多项式的导数IMPLICIT NONE!!函数变量声明=LONG_INT)::n !Legendre多项式阶次(高斯点数)=LONG_INT)::I !循环索引变量=DP)::x !用于迭代计算的采样点=DP)::P(n) !n阶Legendre多项式在采样点x处的取值==1)THEN !阶次为1的特殊情况DPoly_Legendre=ONERETURN==2)THEN !阶次为2的特殊情况DPoly_Legendre=THR*xRETURNELSE !阶次大于2的一般情况!!构造1阶和2阶Legendre公式=x !1阶Legendre公式=HLF*(THR*x*x-ONE) !2阶Legendre多项式!!利用递推公式构造n阶Legendre公式DO I=3,n=((2*I-1)*x*P(I-1)-(I-1)*P(I-2))/IEND DO!!利用递推公式构造n阶Legendre公式的导数DPoly_Legendre=n*(x*P(n)-P(n-1))/(x*x-ONE)END IFEND FUNCTION DPoly_Legendre=DP) FUNCTION Legendre_Root_Bisection(n,a,b)!!函数用于二分法求解Legendre多项式在有根区间[a,b]内的零根IMPLICIT NONE=LONG_INT),INTENT(IN)::n !高斯点数量=DP)::a,b,cDOc=(a+b)*HLFIF(Poly_Legendre(n,a)*Poly_Legendre(n,c)<0)THEN!!零根位于区间[a,c]b=cELSE!!零根位于区间[c,b]a=cEND IFIF((b-a)<EPS)THENEXITEND IFEND DOLegendre_Root_Bisection=cEND FUNCTION Legendre_Root_BisectionSUBROUTINE Gauss_Legendre_Weight(n,X,A)!!基于二分法求解Gauss-Legendre积分公式的高斯点位置并计算权重IMPLICIT NONE!!函数变量声明=LONG_INT),INTENT(IN)::n !高斯点数量=LONG_INT)::Ns !子区间数量(高斯点数量十倍)=DP)::X(n) !高斯点位置=DP)::A(n) !高斯点权重=DP)::B(50*n+1) !子区间位置=DP)::Delta_B !子区间范围=LONG_INT)::I !循环索引变量=LONG_INT)::J !高斯点计数索引!!设置子区间数量(高斯点数量的10倍)Ns=50*n!!直接计算只有一个高斯点的特殊情况==1)THEN=ZR=TWO/((ONE-X(n)*X(n))*&(DPoly_Legendre(n,X(n)))*&(DPoly_Legendre(n,X(n))))RETURNEND IF!!考虑高斯点大于等于2的一般情况!!根据高斯点数量计算子区间Delta_B=TWO/REAL(Ns,DP)=-ONE=ONEDO I=2,Ns=B(I-1)+Delta_BEND DO!!初始化高斯点位置及权值:n)=ZR:n)=ZR!!采用二分法迭代计算有零根子区间内的零根J=0DO I=1,NsIF(Poly_Legendre(n,B(I))*Poly_Legendre(n,B(I+1))<0)THEN!!当前区间存在零根 采用二分法迭代求解J=J+1=Legendre_Root_Bisection(n,B(I),B(I+1))END IFEND DO!!根据Legendre多项式的对称性整理所得零根==0)THEN !高斯点为偶数DO I=n/2+1,n=HLF*(X(I)-X(n-I+1)) !根据对称点将零根平均化=-X(I)END DOELSE !高斯点为奇数=ZR !中间的零根必为零DO I=(n+1)/2+1,n=HLF*(X(I)-X(n-I+1))=-X(I)END DOEND IF!!计算权值DO I=1,n=TWO/((1-X(I)*X(I))*(DPoly_Legendre(n,X(I)))*&(DPoly_Legendre(n,X(I))))END DOEND SUBROUTINE Gauss_Legendre_WeightEND MODULE Module_Gauss_Legendre_WeightPROGRAM Gauss_Legendre_IntegrationUSE Module_Gauss_Legendre_WeightIMPLICIT NONE!!变量声明=LONG_INT),PARAMETER::Num_Gauss=5 !高斯点数量(≥1)=DP)::Gauss_X(Num_Gauss) !高斯点位置=DP)::Gauss_Weight(Num_Gauss) !高斯点权重=LONG_INT)::I !循环索引变量=DP)::Sum !积分结果!!计算高斯点位置及权值CALL Gauss_Legendre_Weight(Num_Gauss,Gauss_X,Gauss_Weight)!!输出高斯点位置及权值DO I=1,Num_Gauss"高斯点序号:",I"高斯点位置:",Gauss_X(I)"高斯点权值:",Gauss_Weight(I)END DO!!计算积分Sum=ZRDO I=1,Num_GaussSum=Sum+y(Gauss_X(I))*Gauss_Weight(I)END DO"积分结果为:",Sum" Press any key to exist "READ(*,*)CONTAINS!定义被积函数=DP) FUNCTION y(x)IMPLICIT NONE=DP)::Xy=x*x/SQRT(1-x*x)END FUNCTION yEND PROGRAM Gauss_Legendre_Integration
