Gauss-Legendre Quadrature - Python实现

2022-04-08 02:05:30

算法特征:
①. 插值型数值积分; ②. 求积节点取Legendre多项式之零点; ③. $n + 1$个求积节点对应$2n + 1$的代数精度

算法推导:
积分区间$[a, b]$上带权函数的插值型数值积分公式如下:
\begin{equation}
\int_a^b\rho (x)f(x)\mathrm{d}x \approx \sum_{i=0}^n A_i f(x_i)
\label{eq_1}
\end{equation}
根据Gauss Quadrature之定义, 当存在$n+1$个求积节点时, 数值积分式$\eqref{eq_1}$具有最高代数精度$2n + 1$. 由此可建立如下代数方程组(此处取权函数$\rho (x) = 1$):
\begin{equation}
\int_a^b x^k dx = \sum_{i=0}^n A_i x_i^k, \quad \text{where }k = 0, 1, 2, \cdots, 2n + 1
\label{eq_2}
\end{equation}
求解上述方程组, 即可获得$n+1$个求积节点$x_i$及求积系数$A_i$. 显然, 硬解此非线性方程组并不是一个好主意, 尤其是当$n$取值较大时.
不过, 在积分区间$[-1, 1]$上, 此$n+1$个求积节点$x_i$可取$n+1$阶Legendre多项式$P_{n+1}(x)$之零点, 相应数值积分也称Gauss-Legendre Quadrature. $n$阶Legendre多项式表示如下:
\begin{equation}
P_n(x) = \frac{1}{2^n n!}\frac{\mathrm{d}^n}{\mathrm{d}x^n}[(x^2 - 1)^n]
\label{eq_3}
\end{equation}
根据上述形式计算并获取$n+1$个求积节点$x_i$后, 代入式$\eqref{eq_2}$, 任选其中$n+1$个方程, 组成新的线性方程组, 求解此线性方程组, 即可获取$n+1$个求积系数$A_i$.
实际使用Gauss-Legendre Quadrature时, 需将任意积分区间$[a, b]$换算至$[-1, 1]$, 即:
\begin{equation}
\int_a^b f(x) \mathrm{d}x = \frac{b - a}{2}\int_{-1}^1 f\left (\frac{b - a}{2}t + \frac{a + b}{2} \right ) \mathrm{d}t
\label{eq_4}
\end{equation}
以下给出5组Legendre多项式之零点及相关求积系数:

阶数$n$	零点$x_i$	求积系数$A_i$
1	$0$	$2$
2	$\pm \sqrt{1 / 3}$	$1$
3	$\pm \sqrt{3 / 5}$ $0$	$5 / 9$ $8 / 9$
6	$\pm 0.238619186081526$ $\pm 0.661209386472165$ $\pm 0.932469514199394$	$0.467913934574257$ $0.360761573028128$ $0.171324492415988$
10	$\pm 0.973906528517240$ $\pm 0.433395394129334$ $\pm 0.865063366688893$ $\pm 0.148874338981367$ $\pm 0.679409568299053$	$0.066671344307672$ $0.269266719309847$ $0.149451349151147$ $0.295524224714896$ $0.219086362515885$

代码实现:
Part Ⅰ:
现以如下定积分为例进行算法实施:
\begin{equation*}
\int_{0}^1 x^2\mathrm{e}^x \mathrm{d}x = \int_{-1}^1\frac{1}{8}(t + 1)^2\mathrm{e}^{(t + 1) / 2}\mathrm{d}t = \mathrm{e} - 2
\end{equation*}
Part Ⅱ:
Gauss-Legendre Quadrature实现如下:

 1 # Gauss-Legendre Quadrature之实现
 2 
 3 import numpy
 4 
 5 
 6 def getGaussParams(num=6):
 7     if num == 1:
 8         X = numpy.array([0])
 9         A = numpy.array([2])
10     elif num == 2:
11         X = numpy.array([numpy.math.sqrt(1/3), -numpy.math.sqrt(1/3)])
12         A = numpy.array([1, 1])
13     elif num == 3:
14         X = numpy.array([numpy.math.sqrt(3/5), -numpy.math.sqrt(3/5), 0])
15         A = numpy.array([5/9, 5/9, 8/9])
16     elif num == 6:
17         X = numpy.array([0.238619186081526, -0.238619186081526, 0.661209386472165, -0.661209386472165, 0.932469514199394, -0.932469514199394])
18         A = numpy.array([0.467913934574257, 0.467913934574257, 0.360761573028128, 0.360761573028128, 0.171324492415988, 0.171324492415988])
19     elif num == 10:
20         X = numpy.array([0.973906528517240, -0.973906528517240, 0.433395394129334, -0.433395394129334, 0.865063366688893, -0.865063366688893, \
21             0.148874338981367, -0.148874338981367, 0.679409568299053, -0.679409568299053])
22         A = numpy.array([0.066671344307672, 0.066671344307672, 0.269266719309847, 0.269266719309847, 0.149451349151147, 0.149451349151147, \
23             0.295524224714896, 0.295524224714896, 0.219086362515885, 0.219086362515885])
24     else:
25         raise Exception(">>> Unsupported num = {} <<<".format(num))
26     return X, A
27     
28 
29 def getGaussQuadrature(func, a, b, X, A):
30     '''
31     func: 原始被积函数
32     a: 积分区间下边界
33     b: 积分区间上边界
34     X: 求积节点数组
35     A: 求积系数数组
36     '''
37     term1 = (b - a) / 2
38     term2 = (a + b) / 2
39     term3 = func(term1 * X + term2)
40     term4 = numpy.sum(A * term3)
41     val = term1 * term4
42     return val
43 
44 
45 def testFunc(x):
46     val = x ** 2 * numpy.exp(x)
47     return val    
48     
49     
50     
51 if __name__ == "__main__":
52     X, A = getGaussParams(10)
53     a, b = 0, 1
54     numerSol = getGaussQuadrature(testFunc, 0, 1, X, A)
55     exactSol = numpy.e - 2
56     relatErr = (exactSol - numerSol) / exactSol
57     print("numerical: {}; exact: {}; relative error: {}".format(numerSol, exactSol, relatErr))

View Code

结果展示:

阶数$n$	numerical integral	exact integral	relative error
1	$0.412180317675032$	$0.718281828459045$	$0.4261579489497921$
2	$0.711941774242270$	$0.718281828459045$	$0.008826694433265773$
3	$0.718251779040964$	$0.718281828459045$	$4.1835136140953615e-05$
6	$0.718281828489842$	$0.718281828459045$	$-4.2875662903868903e-11$
10	$0.718281828458231$	$0.718281828459045$	$1.1338997850082094e-12$

可以看到, 随着Legendre多项式阶数$n$的提升, 即增加插值节点数, Gauss-Legendre Quadrature之相对误差逐渐缩小, 并最终获得较好的数值积分效果.

使用建议:
①. 被积函数在求积节点上的函数值是可得的.
参考文档:
吕书龙, 刘文丽, 梁飞豹. Legendre多项式零点的一种求解方法及应用. 福州大学学报(自然科学版), 2011, 39(3): 334-338.

码农公寓

相关文章