高斯求积简介


#1

高斯求积

高斯求积是常用的数值求积方法,这里主要介绍一下高斯-勒让得求积法则。我们都知道函数 f(x) 的积分可以用 n 个矩形的面积来逼近(如下图),一般会在函数 f(x) 上等间距采 n 个点的值 f(x_1), f(x_2), ... , f(x_n) ,然后分别乘以间隔(权重)得到面积,再累加求和得到积分。
image

而数值求积的思路是找到一些点 x_i 以及合适的权重 w_i ,用 \sum w_if(x_i) 来逼近 f(x)[-1,1] 的定积分。倘若我们随机选 4 个点 x_1, x_2, x_3, x_4 ,那么要使👇对所有 f(x) 都成立:

\int_{-1}^{1} f(x)dx \approx w_1f(x_1) + w_2f(x_2) + w_3f(x_3) + w_4f(x_4)

现在,我们的目标是用 n 个点完美积分 n-1 阶多项式。可以分别令 f(x) 为张成多项式空间的一组基 (1, x, x^2, x^3)

2 = \int_{-1}^{1} 1dx \approx w_1 + w_2 + w_3 + w_4 \\ 0 = \int_{-1}^{1} xdx \approx w_1x_1 + w_2x_2 + w_3x_3 + w_4x_4 \\ \frac{2}{3} = \int_{-1}^{1} x^2dx \approx w_1x_1^2 + w_2x_2^2 + w_3x_3^2 + w_4x_4^2 \\ 0 = \int_{-1}^{1} x^3dx \approx w_1x_1^3 + w_2x_2^3 + w_3x_3^3 + w_4x_4^3 \\

这样就得到了一个线性方程组,4 个方程 4个权重 w_1, w_2, w_3, w_4 未知数,正好可以解出。如果我们选择 x_i[-1,1] 上均匀分布,那么这就是牛顿-柯特斯积分。这里有一个 Julia 的例子:

julia> V(x) = [x[j]^(i-1) for i in eachindex(x), j in eachindex(x)]
V (generic function with 1 method)

julia> x = range(-1, stop=1, length=4)
-1.0:0.6666666666666666:1.0

julia> w = V(x)\[2, 0, 2/3, 0]
4-element Array{Float64,1}:
 0.24999999999999978
 0.7500000000000002
 0.7500000000000002
 0.24999999999999983

julia> f(x) = 1+15x+2x^2+12x^3
f (generic function with 1 method)

julia> w'f.(x)
3.333333333333334

F(x)|_{-1}^{1} = (x + \frac{15}{2} x^2 + \frac{2}{3} x^3 + 3x^4)|_{-1}^{1} = \frac{10}{3}

可以看到对于多项式,可以完美积分。下面是一个非多项式的例子,也可以得到不错的精度。

julia> w'cos.(x)  # cos(x) 在 [-1,1] 的积分是 2sin(1)
1.6875865724061767

julia> 2sin(1)
1.682941969615793

julia> w'sin.(x)  # sin(x) 在 [-1,1] 的积分是 0
2.7755575615628914e-17

高斯-勒让得求积

上面选取的点,权重以及基 (1, x, x^2, x^3) 都比较随意,合理的选择就可以做到用更少的点来获得更精确的估计。首先 (1, x, x^2, x^3) 不是正交基(我们现在定义多项式内积的定义是: \langle p, q \rangle = \int_{-1}^{1} p(x) q(x) dx ),我们可以对其作格拉姆-施密特正交化(Gram-Schmidt)正交化得到另一组基 (L_1, L_2, ..., L_n)

任何一个 2n-1 阶的多项式都可以被分解成 p_{2n-1}(x) = q_{n-1}(x)L_n(x) + r_{n-1}(x) ,其中下标表示多项式的阶数。由于 L_n(x) 是由一组正交基 (L_1, L_2, ..., L_n) 的线性组合, q_{n-1}(x) 又可以表示为 (L_1, L_2, ..., L_{n-1}) 的线性组合,因此 L_n(x)q_{n-1}(x) 正交,那么我们有积分

\int_{-1}^{1} p_{2n-1}(x) dx = \int_{-1}^{1} q_{n-1}(x)L_n(x)dx + \int_{-1}^{1}r_{n-1}(x) dx \\ = 0 + \int_{-1}^{1} r_{n-1}(x) dx.

在这个基下, p_{2n-1}r_{n-1}(x) 的积分相等,我们只需要对 n-1 阶多项式积分就能得到 2n-1 阶多项式的积分。

现在我们已经选好了基,下一步是选点,按照高斯求积的思路,选取的这些点最终要使 \sum w_if(x_i) 完美计算 p_{2n-1}(x) 的积分,所以 p_{2n-1}(x)r_{n-1}(x) 多项式应该在这些点上相等。

如果我们取 L_n(x) 的根,这个时候我们有 q_{n-1}(x_i)L_n(x_i) = 0 ,和 p_{2n-1}(x_i)=r_{n-1}(x_i) ,其中 x_i 就是积分要取值的点。 那么就应该将点选为满足 L_n(x_i)=0 条件的。

下面是用Golub-Welsch算法计算高斯-勒让得的点和权重。

julia> using LinearAlgebra

julia> function gauss_quad(n)
           β = @. .5/sqrt(1-(2*(1:n-1))^(-2.))
           T = SymTridiagonal(zeros(n), β)
           D, V = eigen(T)
           i = sortperm(D); x = D[i]
           w = 2*view(V, 1, i).^2
           x, w
       end
gauss_quad (generic function with 1 method)

julia> x, w = gauss_quad(4)
([-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855])

julia> w'cos.(x)
1.6829416886959736

julia> 2sin(1)
1.682941969615793

我们可以看到虽然仅仅用了四个点,但是计算得十分准确。

数值实验

using Plots; gr()
using SpecialFunctions
using LinearAlgebra
function newton_cotes(m)
    V(x) = [x[j]^(i-1) for i in eachindex(x), j in eachindex(x)] # transposed Vandermonde matrix
    rhs(m) = [((-1)^n + 1)/(n+1) for n in 0:m-1]
    x = range(-1, stop=1, length=m)
    w = V(x)\rhs(m)
    return x, w
end
function gauss_quad(n)
    β = @. .5/sqrt(1-(2*(1:n-1))^(-2.))
    T = SymTridiagonal(zeros(n), β)
    D, V = eigen(T)
    i = sortperm(D); x = D[i]
    w = 2*view(V, 1, i).^2
    return x, w
end
f(x) = ℯ^(-x^2)
sol = √pi*erf(1) # integral of f from -1 to 1
ns = 2:50
newton_errs = map(ns) do n
    x, w = newton_cotes(n)
    return err = abs(sol - w'f.(x))
end
gauss_errs = map(ns) do n
    x, w = gauss_quad(n)
    return err = abs(sol - w'f.(x)) + eps() # machine epsilon to aviod exact 0
end
plot(ns, newton_errs, lab="Newton-Cotes")
plot!(ns, gauss_errs, lab="Gauss-Legendre", dpi=300)
yaxis!("Error", :log10)
xaxis!("Number of nodes")
savefig("$(homedir())/quad_plt.png")

值得注意的是牛顿-柯特斯(Newton-Cotes)收敛的速度不但比高斯-勒让得(Gauss-Legendre)收敛的速度慢,而且因为龙格现象Runge’s phenomenon)导致牛顿-柯特斯在浮点运算下不会收敛。

延伸

注意高斯求积要求积分区间为 [-1,1] ,实际应用需要用变区间法则将积分区间转换到这个范围。

\int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}x_{i}+{\frac {a+b}{2}}\right).

感谢

本文大部分由 @Gnimuc 整理和撰写, @scheme 修改。


#2

#3

有没有人对上面写的Golub-Welsch算法或者计算Newton-Cotes的方法有疑问?