- 算法特征:
①. 插值型数值积分; ②. 求积节点取Legendre多项式之零点; ③. $n + 1$个求积节点对应$2n + 1$的代数精度 - 算法推导:
积分区间$[a, b]$上带权函数的插值型数值积分公式如下:
egin{equation}
int_a^b ho (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$. 由此可建立如下代数方程组(此处取权函数$ ho (x) = 1$):
egin{equation}
int_a^b x^k dx = sum_{i=0}^n A_i x_i^k, quad ext{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多项式表示如下:
egin{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]$, 即:
egin{equation}
int_a^b f(x) mathrm{d}x = frac{b - a}{2}int_{-1}^1 fleft (frac{b - a}{2}t + frac{a + b}{2} ight ) 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 Ⅰ:
现以如下定积分为例进行算法实施:
egin{equation*}
int_{0}^1 x^2mathrm{e}^x mathrm{d}x = int_{-1}^1frac{1}{8}(t + 1)^2mathrm{e}^{(t + 1) / 2}mathrm{d}t = mathrm{e} - 2
end{equation*}
Part Ⅱ:
Gauss-Legendre Quadrature实现如下:
View Code1 # 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))
- 结果展示:
可以看到, 随着Legendre多项式阶数$n$的提升, 即增加插值节点数, Gauss-Legendre Quadrature之相对误差逐渐缩小, 并最终获得较好的数值积分效果.阶数$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多项式零点的一种求解方法及应用. 福州大学学报(自然科学版), 2011, 39(3): 334-338.
