• {{article.mytitle}}

    • 高斯积分及python实现

    1. 高斯积分及python实现

    1.1. 返回首页

    首页地址https://kz16.top/na
    本文来源https://kz16.top/na/Guass

    在一维情况下，由于等距插值会有龙格现象存在，因此等距插值型求积公式不太稳定。 这里我们将讨论非等距插值型求积公式。

    1.2. 高斯积分理论

    不失一般性，考虑区间(−1,1)(-1,1)上的插值求积公式
    ∫−11f(x)dx≈∑k=0nAkf(xk)(1)\tag{1} \int^{1}_{-1} f(x) {\rm d}x \approx \sum\limits^{n}_{k=0}A_kf(x_k)

    其中{xk,Ak}k=0,1,⋯ ,n\{x_k,A_k\}_{k=0,1,\cdots,n}为2n+22n+2个待定参数，如果这2n+22n+2个参数，使得(1)式对任何次数不超过2n+12n+1的多项式f(x)f(x)准确成立，称该公式具有2n+12n+1次代数精度，也称这类求积公式为高斯型积分公式，同时称其节点{xk}k=0,1,⋯ ,n\{x_k\}_{k=0,1,\cdots,n}为高斯点，相应的权系数{Ak}k=0,1,⋯ ,n\{A_k\}_{k=0,1,\cdots,n}为高斯权。

    1.2.1. 定理及结论

    • 区间(−1,1)(-1,1)上的n+1n+1次正交多项式（勒让德多项式）的零点恰好是高斯点{xk}k=0,1,⋯ ,n\{x_k\}_{k=0,1,\cdots,n}.

    • 高斯权为正数，即
      Ak=∫−11lk2(x)dx>0,k=0,1,⋯ ,nA_k=\int^{1}_{-1} l_k^2(x) {\rm d}x>0, \quad k=0,1,\cdots,n

      其中

      lk(x)=(x−x0)(x−x1)⋯(x−xk−1)(x−xk+1)⋯(x−xn)(xk−x0)(xk−x1)⋯(xk−xk−1)(xk−xk+1)⋯(xk−xn)l_k(x)=\frac{(x - x_0)(x - x_1)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_{n})}{(x_k - x_0)(x_k - x_1)\cdots (x_k -x_{k-1})(x_k -x_{k+1})\cdots (x_k -x_{n})}

    1.2.2. 一维一般有限区间上的高斯积分

    对于一般区间(a,b)(a,b)上的求积，只需做变量替换

    x(t)=12(a+b)+12(b−a)tx(t)=\frac{1}{2}(a+b) + \frac{1}{2}(b-a)t

    有

    ∫abf(x)dx=b−a2∫−11f[12(a+b)+12(b−a)t]dt\int^{b}_{a}f(x){\rm d}x=\frac{b-a}{2}\int^{1}_{-1}f \left[\frac{1}{2}(a+b) + \frac{1}{2}(b-a)t \right] {\rm d}t

    对于上式右端可以用前面的求积公式。

    1.2.3. 二维情形

    不失一般性，区域(−1,1)2(-1,1)^2上的积分
    ∫−11∫−11f(x,y)dxdy≈∑i=0n∑j=0nωiωjf(xi,xj).(2)\tag{2} \int_{-1}^1 \int_{-1}^1 f(x,y) {\rm d}x {\rm d}y \approx \sum_{i=0}^n \sum_{j=0}^n \omega_i \omega_j f(x_i, x_j) .

    这里{xi,wi}i=0n\{ x_i,w_i\}_{i=0}^n及{xj,wj}j=0n\{ x_j,w_j\}_{j=0}^n为前面讨论的一维高斯积分点和权。

    1.2.4. 特征值法求高斯点和权

    给出矩阵G\mathbb G定义如下
    G=[0u10⋯00u10u2⋯000u20⋯00⋮⋮⋮⋮⋮000⋯0un000⋯un0],uk=k4k2−1,k=1,2,⋯ ,n\mathbb G = \begin{bmatrix} 0 & u_1 & 0 & \cdots & 0 & 0\\ u_1 & 0 & u_2 & \cdots & 0 & 0\\ 0 & u_2 & 0 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & u_n\\ 0 & 0 & 0 & \cdots & u_n & 0\\ \end{bmatrix} ,u_k = \frac{k}{\sqrt{4k^2 - 1}},k=1,2,\cdots,n

    可以证明，矩阵G\mathbb G所有特征值构成的集合λ={λ0,λ1,⋯ ,λn}\bm \lambda = \{\lambda_0,\lambda_1,\cdots,\lambda_n\}恰是所有的高斯点集合，因此只要将λ\bm\lambda按从小到大的顺序排序就可得到排好序的高斯点，设v\bm v是矩阵G\mathbb G最小的特征值的特征向量，将特征向量v\bm v的所有分量元素平方后乘以2构成的集合WW就是高斯系数集合，只要将WW按照λ\bm \lambda 的排序角标排序就能得到对应的高斯系数。

    python代码如下

    import numpy as np
def guasspw(n):
    ran = np.arange(1, n+1)
    u = ran/np.sqrt((2*ran)**2 - 1)
    [bp, vc] = np.linalg.eig(np.diag(u, 1) + np.diag(u, -1))
    wf = 2*vc[0] ** 2
    wf[n - n//2:] = wf[n//2::-1]
    return np.sort(bp), wf
print(guasspw(3))

    当然我们也可以直接用numpy中的函数

    import numpy as np
print(np.polynomial.legendre.leggauss(4))

    两者的结果都是

    (array([-0.86113631, -0.33998104,  0.33998104,  0.86113631]), 
 array([0.34785485, 0.65214515, 0.65214515, 0.34785485]))

    1.3. 一维情形数值例子

    1.3.1. 问题

    用高斯积分计算
    ∫0141+x2dx\int_{0}^1 \frac{4}{1+x^2} {\rm d}x

    1.3.2. 一维数值积分分析

    高斯数值积分运算实质是加权求和，也可以认为是两个向量点乘，因此可以直接用python第三方包numpy中的dot函数很高效计算。当然也可以用for循环去实现，但是matlab或python的for循环比起向量点乘效率要低很多。

    1.3.3. 一维数值积分python代码

    import numpy as np
p, w = np.polynomial.legendre.leggauss(100)
f = lambda x: 4.0/(1 + x*x)
ans = 0.5*f(0.5*p + 0.5).dot(w)
print(ans)

    结果为

    3.1415926535897927

    1.4. 二维情形数值例子

    1.4.1. 问题

    用高斯积分计算
    ∫−11∫0yycos⁡(xy)dxdy\int_{-1}^1\int_{0}^{y}y \cos(xy){\rm d} x {\rm d}y

    1.4.2. 二维数值积分分析

    对于该二重积分，我们可以分两步，先考虑内层xx方向的积分，
    利用变换x=12y(t+1)x=\frac{1}{2}y(t+1)，则有

    ∫0yycos⁡(xy)dx=∫−1112y2cos⁡[12y2(t+1)]dt\begin{aligned} \int_{0}^{y}y \cos(xy){\rm d} x = \int_{-1}^1 \frac{1}{2}y^2 \cos\left[\frac{1}{2}y^2(t+1)\right] {\rm d} t \end{aligned}

    因此有

    ∫−11∫0y2ycos⁡(xy)dxdy=∫−11∫−1112y2cos⁡[12y2(t+1)]dtdy≈∑i=0n∑j=0nωiωjf(xi,xj).(3)\tag{3} \begin{aligned} \int_{-1}^1\int_{0}^{y^2}y \cos(xy){\rm d} x {\rm d}y &= \int_{-1}^1 \int_{-1}^1 \frac{1}{2}y^2 \cos\left[\frac{1}{2}y^2(t+1)\right] {\rm d} t {\rm d}y \\ & \approx \sum_{i=0}^n \sum_{j=0}^n \omega_i \omega_j f(x_i, x_j) . \end{aligned}

    这里f(y,t)=12y2cos⁡[12y2(t+1)]f(y, t) = \frac{1}{2}y^2 \cos\left[\frac{1}{2}y^2(t+1)\right]，{xi,wi}i=0n\{ x_i,w_i\}_{i=0}^n及{xj,wj}j=0n\{ x_j,w_j\}_{j=0}^n为前面讨论的一维高斯积分点和权。

    注： 上面(3)式的最后一项其实也是可以矩阵乘法形式，因此依旧可以免去for循环高效计算

    当然，不难计算该积分可以转化为一维数值积分，因此可以很方便对比数值结果
    ∫−11∫0yycos⁡(xy)dxdy=∫−11sin⁡(xy)∣0ydy=∫−11sin⁡(y2)dy.\begin{aligned} \int_{-1}^1\int_{0}^{y}y \cos(xy){\rm d} x {\rm d}y &= \int_{-1}^1 \sin(xy) \big|_{0}^y {\rm d}y \\ & = \int_{-1}^1 \sin(y^2) {\rm d}y . \end{aligned}

    1.4.3. 二维数值积分python代码

    import numpy as np
p, w = np.polynomial.legendre.leggauss(100)
f = lambda y: np.sin(y*y)
ans = f(p).dot(w)
print("integrate 1d: ", ans)

p2d = np.meshgrid(p, p)
w2d = np.outer(w, w)
f2d = lambda x: 0.5 * x[1]**2 * np.cos(0.5 * x[1]**2 * (x[0]+1))
ans = f2d(p2d)*w2d
print("integrate 2d: ", ans.sum())

    结果为

    integrate 1d:  0.6205366034467572
integrate 2d:  0.6205366034467573

    1.5. 温故知新

    通过本次教程，我们知道了

    • 高斯积分是非等距插值型积分公式，并且积分精度最高
    • 高斯点是勒让德多项式的零点，高斯权都是正数
    • 一维一般积分区间，可通过积分区间变换转化为标准积分区间(−1,1)(-1,1)上的积分
    • 二维一般积分区间，可通过积分区域变换转化为标准积分区域(−1,1)2(-1,1)^2上的积分
    • 无论是一维还是二维数值积分，都可以归为矩阵向量的运算

    1.6. 返回首页

    首页地址https://kz16.top/na
    本文来源https://kz16.top/na/Guass

    仅供参考，部分漏洞在所难免
    欢迎转载，转载请指明来源，
    请勿用于商业
    作者：Zhou Bingzhen
    2019.11.21