Gauss-Legendre Quadrature - Python实现

• 算法特征:
  ①. 插值型数值积分; ②. 求积节点取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.

标签：Legendre,Python,equation,求积,Quadrature,array,numpy,pm
来源： https://www.cnblogs.com/xxhbdk/p/14203877.html