首页/文章/ 详情

有限元中的高斯-勒让德积分

1年前浏览732
高斯-勒让德积分是有限元中最常见的数值积分方法之一,在有限元中有着广泛的应用。实际上,关于该积分方法的书籍和公 众 号文章也已经较多,本文主要是基于现有的理论,基本上把该方法的具体理论和使用重复了一遍,另外基于常用的数值计算库PETSc描述下在PETSc中如何使用高斯勒让德积分(本文后续都将高斯勒让德积分简称为高斯积分)。
高斯积分的具体公式如下:

上式即是n点高斯积分的计算公式,其中xi一般称作高斯点的位置,wi称作权重。
当积分是三维空间时,其表达式如下:

从该公式可知,其本质上是将上下界为-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库也直接给出了高斯积分的值和权重。
以下是高斯积分点积分多项式的一个例子:

上式中,x的最高次数是3,因此我们采用2点高斯积分进行积分(2点高斯积分可以准确积分2*2-1阶多项式),积分点位置和权重分别为(+-0.577350269189626)和(1.0,1.0) 。

而准确解:

显然,高斯积分给出了准确解。当然,如果采用多于2点的高斯积分,也能给出准确解。
高斯积分点的位置和权重可以采用多项式待定系数求出,还是以两点高斯积分为例,具体过程如下:

由于对任意的a,b,c,e均成立,因此:

求解上述方程组即可得到积分点位置和权重值。
在PETSc中,高斯积分函数如下:
    #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;}
        输出:

        以上,就是本文的主要内容,感谢您的阅读!

        来源:有限元术
        理论
        著作权归作者所有,欢迎分享,未经许可,不得转载
        首次发布时间:2023-05-25
        最近编辑:1年前
        寒江雪_123
        硕士 | cae工程师 签名征集中
        获赞 48粉丝 98文章 46课程 9
        点赞
        收藏
        未登录
        还没有评论
        课程
        培训
        服务
        行家
        VIP会员 学习 福利任务 兑换礼品
        下载APP
        联系我们
        帮助与反馈