高斯积分及python实现

Author : zbzhen,        Modified : Tue Feb 21 10:54:44 2023

1. 高斯积分及python实现

1.1. 返回首页

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

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

1.2. 高斯积分理论

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

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

1.2.1. 定理及结论

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

• 高斯权为正数，即
  $$A_k=\int^{1}_{-1} l_k^2(x) {\rm d}x>0, \quad k=0,1,\cdots,n$$

  其中

  $$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)$上的求积，只需做变量替换

$$x(t)=\frac{1}{2}(a+b) + \frac{1}{2}(b-a)t$$

有

$$\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$上的积分
$$\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) .$$

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

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

给出矩阵$\mathbb G$定义如下
$$\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$$

可以证明，矩阵$\mathbb G$所有特征值构成的集合$\bm \lambda = \{\lambda_0,\lambda_1,\cdots,\lambda_n\}$恰是所有的高斯点集合，因此只要将$\bm\lambda$按从小到大的顺序排序就可得到排好序的高斯点，设$\bm v$是矩阵$\mathbb G$最小的特征值的特征向量，将特征向量$\bm v$的所有分量元素平方后乘以2构成的集合$W$就是高斯系数集合，只要将$W$按照$\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. 问题

用高斯积分计算
$$\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. 问题

用高斯积分计算
$$\int_{-1}^1\int_{0}^{y}y \cos(xy){\rm d} x {\rm d}y$$

1.4.2. 二维数值积分分析

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

$$\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}$$

因此有

$$\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) = \frac{1}{2}y^2 \cos\left[\frac{1}{2}y^2(t+1)\right]$，$\{ x_i,w_i\}_{i=0}^n$及$\{ x_j,w_j\}_{j=0}^n$为前面讨论的一维高斯积分点和权。

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

当然，不难计算该积分可以转化为一维数值积分，因此可以很方便对比数值结果
$$\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)^2$上的积分
• 无论是一维还是二维数值积分，都可以归为矩阵向量的运算

1.6. 返回首页

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

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