高斯-勒让德积分是有限元中最常见的数值积分方法之一,在有限元中有着广泛的应用。实际上,关于该积分方法的书籍和公 众 号文章也已经较多,本文主要是基于现有的理论,基本上把该方法的具体理论和使用重复了一遍,另外基于常用的数值计算库PETSc描述下在PETSc中如何使用高斯勒让德积分(本文后续都将高斯勒让德积分简称为高斯积分)。![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462334920_d6j5n1.jpg?imageView2/0)
上式即是n点高斯积分的计算公式,其中xi一般称作高斯点的位置,wi称作权重。![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462334656_x0at3w.jpg?imageView2/0)
从该公式可知,其本质上是将上下界为-1和1的积分计算转换为多项求和计算。当函数表达式比较复杂时,f(x)的原函数可能难以求出,而采用高斯积分,其省去了求f(x)原函数,只需要将数值代入f(x)的表达式即可求解。到目前为止,高斯积分的公式已经介绍完成,那么有两个最直接最现实的问题出现了:(1)f(x)的表达式是什么形式时适合采用高斯积分,精度怎么样;(2)xi和wi的取值是多少。关于(1),实践表明,当f(x)的表达式为多项式时,高斯积分是合适的,并且,n点高斯积分可以准确积分2n-1次多项式。关于(2),xi和wi的取值一般较多的有限元教科书中会给出数值,如果没有给出数值,也可以用多项式手动算出具体值,另外,scipy库,PETSc库也直接给出了高斯积分的值和权重。![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462336714_01y6zp.jpg?imageView2/0)
上式中,x的最高次数是3,因此我们采用2点高斯积分进行积分(2点高斯积分可以准确积分2*2-1阶多项式),积分点位置和权重分别为(+-0.577350269189626)和(1.0,1.0) 。![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462335174_nqq67f.jpg?imageView2/0)
![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462336089_oyrdor.jpg?imageView2/0)
显然,高斯积分给出了准确解。当然,如果采用多于2点的高斯积分,也能给出准确解。高斯积分点的位置和权重可以采用多项式待定系数求出,还是以两点高斯积分为例,具体过程如下:![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462336977_enwq4t.jpg?imageView2/0)
![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462335491_qvomjh.jpg?imageView2/0)
![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462334163_7q3zm2.jpg?imageView2/0)
![](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1685462335805_h3jedn.jpg?imageView2/0)
#include "petscdt.h"
PetscErrorCode PetscDTGaussQuadrature
(PetscInt npoints, PetscReal a, PetscReal b,
PetscReal *x, PetscReal *w)
为了方便使用,可以对其进行封装,封装之后的类如下:#include<iostream>
static char help[] = "Basic gausslegend routines.\n\n";
#include<petscdt.h>
template<int N>
class gausslegend
{
public:
gausslegend<N>()
{
PetscDTGaussQuadrature(npoints, a, b, x, w);
}
void print()
{
std::cout<<npoints<<" points gauss-legend Quadrature:\n";
std::cout<<"position\n";
for(int i=0;i<N;i++)
{
std::cout<<x[i]<<" ";
}
std::cout<<"\nweight\n";
for(int i=0;i<N;i++)
{
std::cout<<w[i]<<" ";
}
std::cout<<std::endl;
}
private:
PetscInt npoints=N;
PetscReal a=-1,b=1;
PetscReal x[N],w[N];
};
int main(int argc, char **argv)
{
PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
gausslegend<2> g1;
g1.print();
gausslegend<4> g2;
g2.print();
PetscCall(PetscFinalize());
return 0;
}