# 高斯求积简介

## 高斯求积

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

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

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


## 高斯-勒让得求积

\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.

julia> using LinearAlgebra

β = @. .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)

([-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
β = @. .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
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")


## 延伸

\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).

## 感谢

8赞

julia> using LinearAlgebra

β = @. .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
gauss_quad (generic function with 1 method)

julia> f(x) = ℯ^(-x^2)
f (generic function with 1 method)

julia> b = 10.; a = 0.;

julia> (b-a)/2 * w' * @. f((b-a)/2*x + (a+b)/2)
0.8862269254527585


\displaystyle \int _{-1}^{1}\tau^m\,d\tau = \sum_{i=1}^{n}w_{i} x_{i}^m, \quad m = 0, 1, ..., n-1

\displaystyle \ell_i(t) = \prod^{n}_{k=1,k\ne i} \frac{t-x_k}{x_i - x_k}, \quad i=1,2,...,n,

\displaystyle \sum^n_{i=1} \ell_i(t)p(x_i) = p(t)

\displaystyle \sum_{i=1}^n w_i x_i^m = \sum_{i=1}^n \int_{-1}^{1} \ell_i(\tau)\, d\tau \cdot x_i^m = \int_{-1}^1 \left(\sum^n_{i=1} \ell_i(\tau) x_i^m \right)\, d\tau = \int_{-1}^{1} \tau^m \, d\tau, \quad m=0,1,...,n-1.

\displaystyle w_i = \int_{-1}^{1} \prod^{m}_{k=1,k\ne i} \frac{t-x_k}{x_i - x_k} \, dt, \quad i=1,2,...,n.

2赞