高斯积分

积分法则

考虑实数区间 \([-1,1]\) 上的积分 $$ \int_{-1}^1 f(\xi)\, d\xi. \tag{1} $$ 一个数值积分法则 \( (\xi_i,\ w_i),i=1,\dots,n\), 即积分点和积分权重, 具有如下形式 $$ \int_{-1}^1 f(\xi) \,d\xi \simeq \sum_{i=1}^{n} \omega_i f \left( \xi_i\right). \tag{2} $$ 我们希望选取的积分法则具有尽可能高的数值精度, 即数值积分 (2) 对尽可能高阶的多项式准确成立. 一般来说, \(2n\) 个参数可以唯一确定一个 \(2n-1\) 次多项式. 因此我们希望可以找到积分法则使得下面的等式成立 $$ \int_{-1}^{1} p(\xi) \,d\xi = \sum_{i=1}^{n} \omega_i p \left( \xi_i\right), \quad \forall\, p\in \mathbb{P}_{2n-1}. \tag{3} $$

一个简单的例子

问题: 寻找积分法则 \((\xi_i,\ w_i),i=1,2\) 使如下等式成立 \[ \int_{-1}^{1} p(\xi) \,d\xi = \omega_1 p(\xi_1)+\omega_2 p(\xi_2), \quad \forall\, p\in \mathbb{P}_{3}([-1,1]). \tag{4} \] 求解: 在等式 (4) 中分别另 \(p(x)=1,\ \xi,\ \xi^2,\ \xi^3\), 得到 4 个方程: $$ \omega_1+\omega_2=\int_{-1}^{1} 1 \,d\xi =2,\quad \omega_1 \xi_1+\omega_2\xi_2=\int_{-1}^{1} \xi \,d\xi =0,\\ \omega_1 \xi_1^2+\omega_2\xi_2^2=\int_{-1}^{1} \xi^2 \,d\xi =\frac{3}{2}, \quad \omega_1 \xi_1^3+\omega_2\xi_2^3=\int_{-1}^{1} \xi^3 \,d\xi =0. $$ 求解上面的非线性方程组 (可以利用 Python 符号求解) 可以得到积分法则: \[ (\xi_1,\ \omega_1)=(-\frac{\sqrt{3}}{3},\ 1)=(-0.577350269189626,\ 1)\\ (\xi_2,\ \omega_2)=(\frac{\sqrt{3}}{3},\ 1)=(0.577350269189626,\ 1). \] 可以验证上面得到的积分法则对所有的 3-阶多项式精确成立.
# 在线运行如下代码求解非线性方程组 https://live.sympy.org/
from sympy import Symbol, nonlinsolve

w1, w2 = Symbol('w1'), Symbol('w2')
x1, x2 = Symbol('x1'), Symbol('x2')
eq1 = w1 + w2 - 2
eq2 = w1 * x1 + w2 * x2
eq3 = w1 * x1**2 + w2 * x2**2 - 2/3
eq4 = w1 * x1**3 + w2 * x2**3

nonlinsolve([eq1, eq2, eq3, eq4], [w1, w2, x1, x2])
# print(nonlinsolve([eq1, eq2, eq3, eq4], [w1, w2, x1, x2]))

Gauss 积分

更高阶的积分法则可以通过上面的方法计算得到, 但是最常用的 Gauss 积分可以很容易地通过正交多项式 (比如 Legendre 多项式) 来计算.

Legendre 多项式

在 \([-1,1]\) 上满足如下性质的多项式被称为 Legendre 多项式:

  1. 对任一 \(n\in \mathbb{N}\), \(P_n(x)\) 是 \(n\)-阶首一多项式 (最高阶项系数为 1);
  2. 对所有阶数小于 \(n\) 的多项式 \(P(x)\), \(\int_{-1}^{1}P(x)P_n(x)\ dx=0\).

可以证明每一个 Legendre 多项式的根都是单根且对称分布在 \((-1, 1)\) 内. 重要地是, 如果我们以 \(n\)-阶 Legendre 多项式 \(P_n(x)\) 的根作为积分节点, 再加上合适的权重, 就可以得到 Gauss-Legendre 积分法则 (具有 \(2n-1\) 阶精度).

一个重要定理

定理: 令 \(\xi_i,i=1,2,\dots,n\) 为 \(P_n(\xi)\) 的 \(n\) 个根, 选取权重 $$ \omega_i = \int_{-1}^{1} L_i(\xi)\ d\xi, $$ 这里 \(L_i(\xi) = \sum_{j=1,j\neq i}^n \frac{\xi-\xi_j} {\xi_i-\xi_j}\) 为 Lagrange 插值基函数. 那么, 对所有小于 \(2n-1\) 阶的多项式 \(P(\xi)\), 我们有 $$ \int_{-1}^{1} P(\xi) \ d\xi = \sum_{i=1}^{n} \omega_i P(\xi_i). $$ 证明: 由多项式分解理论可知, 存在 \(Q(\xi),R(\xi)\in \mathbb{P}{n-1}\) 使得 \(P(\xi)\) 具有如下分解 $$ P(\xi) = Q(\xi)\, P_n(\xi)+R(\xi). $$ 因为 \(P_n(\xi_i)=0\), 所以 \(P(\xi_i)=R(\xi_i)\). 由 Legendre 多项式的正交性质, 有 $$ \int_{-1}^{1} P(\xi)\ d\xi = \int_{-1}^{1} Q(\xi)\, P_n(\xi)+R(\xi) \ dx = \int_{-1}^{1} R(\xi) \ d\xi. \tag{2.1} $$ 因为 \(\mathbb{P}_{n-1}\ni R(\xi)=\sum_{i=1}^n R(\xi_i)L_i(\xi) =\sum_{i=1}^n P(\xi_i)L_i(\xi)\), 可知 $$ \int_{-1}^{1} R(\xi) \ d\xi = \int_{-1}^{1} \left( \sum_{i=1}^n P(\xi_i)L_i(\xi) \right) \ d\xi =\sum_{i=1}^n \left(\int_{-1}^{1} L_i(\xi) \ d\xi \right) P(\xi_i) = \sum_{i=1}^n \omega_i P(\xi_i). $$ 代入 (2.1) 完成定理证明.

一般区间上的数值积分

通常的做法是将一般区间 \([a,b], a,b\in \mathbb{R}\) 上的积分线性变换到参考单元 \([-1,1]\) 上. 具体来说, 积分 (1) 有如下数值积分公式 $$ \int_{a}^b f(x) dx=\frac{b-a}{2} \int_{-1}^{1}f \left( \frac{b-a}{2}\xi+\frac{a+b}{2} \right) \, d\xi \simeq\frac{b-a}{2} \sum_{i}\omega_i f \left( \frac{b-a}{2}\xi_i+\frac{a+b}{2} \right) $$ 因此, 区间 \([a,b]\) 上的积分法则 \((x_i, w_i)\) 为 $$ x_i = \frac{b-a}{2}\xi_i+\frac{a+b}{2},\quad w_i = \frac{b-a}{2} \omega_i, $$ 即 $$ \int_{a}^b f(x) dx \simeq\sum_{i}w_i f (x_i) . $$

附录 A: Gauss-Legendre

n积分点 \(\pm\xi_i\)权重 \(\omega_i\)
21.01.000000000000000
301.333333333333333
1.00.333333333333333
4 0.4472135954999580.833333333333333
1.00.166666666666667
50.00.56888888888888889
0.53846931010568300.47862867049936647
0.906179845938663990.23692688505618909
60.23861918608319690.46791393457269105
0.66120938646626450.36076157304813860
0.93246951420315200.17132449237917035
70.00.41795918367346939
0.40584515137739720.38183005050511894
0.74153118559939440.27970539148927667
0.94910791234275850.12948496616886969
80.18343464249564980.36268378337836198
0.52553240991632900.31370664587788729
0.79666647741362670.22238103445337447
0.96028985649753620.10122853629037626

附录 B: Gauss-Lobatto

n积分点 \(\pm\xi_i\)权重 \(\omega_i\)
21.01.0
30.01.333333333333333
1.00.333333333333333
40.4472135954999580.833333333333333
1.00.166666666666667
50.00.711111111111111
0.6546536707079770.544444444444444
1.00.100000000000000
60.2852315164806450.554858377035486
0.7650553239294650.378474956297847
1.00.066666666666667
70.00.487619047619048
0.4688487934707140.431745381209863
0.8302238962785670.276826047361566
1.00.047619047619048
80.2092992179024790.412458794658704
0.5917001814331420.341122692483504
0.8717401485096070.210704227143506
1.00.035714285714286

参考文献

发表评论

电子邮件地址不会被公开。 必填项已用*标注