高斯勒讓德積分是有限元中最常見的數(shù)值積分方法之一，在有限元中有著廣泛的應(yīng)用。實(shí)際上，關(guān)于該積分方法的書籍和公眾號(hào)文章也已經(jīng)較多，本文主要是基于現(xiàn)有的教程，基本上把該方法的具體理論和使用重復(fù)了一遍，另外基于常用的數(shù)值計(jì)算庫(kù)PETSc描述下在PETSc中如何使用高斯勒讓德積分（本文后續(xù)都將高斯勒讓德積分簡(jiǎn)稱為高斯積分）。

高斯積分的具體公式如下：

上式即是n點(diǎn)高斯積分的計(jì)算公式，其中xi一般稱作高斯點(diǎn)的位置，wi稱作權(quán)重。

當(dāng)積分是三維空間時(shí)，其表達(dá)式如下：

從該公式可知，其本質(zhì)上是將上下界為-1和1的積分計(jì)算轉(zhuǎn)換為多項(xiàng)求和計(jì)算。當(dāng)函數(shù)表達(dá)式比較復(fù)雜時(shí)，f(x)的原函數(shù)可能難以求出，而采用高斯積分，其省去了求f(x)原函數(shù)，只需要將數(shù)值代入f(x)的表達(dá)式即可求解。

到目前為止，高斯積分的公式已經(jīng)介紹完成，那么有兩個(gè)最直接最現(xiàn)實(shí)的問題出現(xiàn)了：（1）f(x)的表達(dá)式是什么形式時(shí)適合采用高斯積分，精度怎么樣；（2）xi和wi的取值是多少。

關(guān)于（1），實(shí)踐表明，當(dāng)f（x）的表達(dá)式為多項(xiàng)式時(shí)，高斯積分是合適的，并且，n點(diǎn)高斯積分可以準(zhǔn)確積分2n-1次多項(xiàng)式。

關(guān)于（2），xi和wi的取值一般較多的有限元教科書中會(huì)給出數(shù)值，如果沒有給出數(shù)值，也可以用多項(xiàng)式手動(dòng)算出具體值，另外，scipy庫(kù)，PETSc庫(kù)也直接給出了高斯積分的值和權(quán)重。

以下是高斯積分點(diǎn)積分多項(xiàng)式的一個(gè)例子：

上式中，x的最高次數(shù)是3，因此我們采用2點(diǎn)高斯積分進(jìn)行積分（2點(diǎn)高斯積分可以準(zhǔn)確積分2*2-1階多項(xiàng)式），積分點(diǎn)位置和權(quán)重分別為（+-0.577350269189626）和(1.0,1.0) 。

而準(zhǔn)確解：

顯然，高斯積分給出了準(zhǔn)確解。當(dāng)然，如果采用多于2點(diǎn)的高斯積分，也能給出準(zhǔn)確解。

高斯積分點(diǎn)的位置和權(quán)重可以采用多項(xiàng)式待定系數(shù)求出，還是以兩點(diǎn)高斯積分為例，具體過程如下：

令

則

由于對(duì)任意的a,b,c,e均成立，因此：

求解上述方程組即可得到積分點(diǎn)位置和權(quán)重值。

在PETSc中，高斯積分函數(shù)如下：

#include "petscdt.h"
PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)

為了方便使用，可以對(duì)其進(jìn)行封裝，封裝之后的類如下：

#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;
}

輸出：

以上，就是本文的主要內(nèi)容，感謝您的閱讀！

歡迎關(guān)注公眾號(hào)  有限元術(shù)