# 高斯求积简介

#1

## 高斯求积

\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

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")

## 延伸

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

#2

#3