Gauss Quadrature数值求积基础与R语言实践
Gauss Quadrature的基本概念
回顾一下Riemann积分的定义
∫
a
b
f
(
x
)
d
x
=
lim
P
→
0
∑
i
=
1
N
f
(
ξ
i
)
Δ
x
i
\int_a^b f(x)dx = \lim_{P \to 0} \sum_{i=1}^N f(\xi_i)\Delta x_i
∫abf(x)dx=P→0limi=1∑Nf(ξi)Δxi
其中
[
a
,
b
]
=
[
x
0
,
x
1
]
∪
[
x
1
,
x
2
]
∪
⋯
∪
[
x
N
−
1
,
x
N
]
[a,b]=[x_0,x_1] \cup [x_1,x_2] \cup \cdots \cup [x_{N-1},x_N]
[a,b]=[x0,x1]∪[x1,x2]∪⋯∪[xN−1,xN]并且
Δ
x
i
=
x
i
−
x
i
−
1
,
i
=
1
,
⋯
,
N
\Delta x_i=x_i-x_{i-1},i=1,\cdots,N
Δxi=xi−xi−1,i=1,⋯,N，
P
=
sup
i
Δ
x
i
P=\sup_{i}\Delta x_i
P=supiΔxi，
∀
ξ
i
∈
[
x
i
−
1
,
x
i
]
\forall \xi_i \in [x_{i-1},x_i]
∀ξi∈[xi−1,xi]；不考虑极限，只是把
N
N
N视为一个比较大的整数，那么上述定积分的近似公式为
∫
a
b
f
(
x
)
d
x
≈
∑
i
=
1
N
f
(
ξ
i
)
Δ
x
i
\int_a^b f(x)dx \approx \sum_{i=1}^N f(\xi_i)\Delta x_i
∫abf(x)dx≈i=1∑Nf(ξi)Δxi
从一般性角度理解，这个近似公式本质上就是一系列函数值
{
f
(
ξ
i
)
}
i
=
1
N
\{f(\xi_i)\}_{i=1}^N
{f(ξi)}i=1N的线性组合。Gauss积分秉承这个思路，它的基本公式为
∫
−
1
1
f
(
x
)
d
x
=
∑
i
=
1
n
w
i
f
(
x
i
)
\int_{-1}^1 f(x)dx=\sum_{i=1}^n w_i f(x_i)
∫−11f(x)dx=i=1∑nwif(xi)
其中
w
i
w_i
wi为权重。如果积分域为
[
a
,
b
]
[a,b]
[a,b]，则可以通过如下变换解决：
∫
a
b
f
(
x
)
d
x
=
∫
−
1
1
b
−
a
2
f
(
b
−
a
2
y
+
a
+
b
2
)
d
y
\int_a^b f(x)dx = \int_{-1}^1 \frac{b-a}{2}f \left( \frac{b-a}{2}y+\frac{a+b}{2} \right)dy
∫abf(x)dx=∫−112b−af(2b−ay+2a+b)dy
R语言实践
在R语言中，Gauss Quadrature求积可以用gaussquad包实现，没有下载的同学可以用下面的代码下载
install.packages("gaussquad")
下面用三个例子介绍这个包的用法。
例1：Gauss-Legendre Quadrature
形如 ∫ − 1 1 f ( x ) d x \int_{-1}^1 f(x)dx ∫−11f(x)dx
的积分可以用Gauss-Legendre Quadrature进行数值求解，用
P
n
P_n
Pn表示Legendre多项式，Gauss-Legendre Quadrature的公式为
∫
−
1
1
f
(
x
)
d
x
≈
∑
i
=
1
n
w
i
f
(
x
i
)
=
∑
i
=
1
n
2
w
i
f
(
x
i
)
(
1
−
x
i
2
)
[
P
n
′
(
x
i
)
]
2
\int_{-1}^1 f(x)dx\approx \sum_{i=1}^n w_if(x_i) = \sum_{i=1}^n \frac{2w_if(x_i)}{(1-x_i^2)[P_n'(x_i)]^2}
∫−11f(x)dx≈i=1∑nwif(xi)=i=1∑n(1−xi2)[Pn′(xi)]22wif(xi)
其中 x i x_i xi代表 P n P_n Pn的第 i i i个零点。
考虑积分 ∫ 0 2 π d θ 2 + cos 2 θ = 2 π 6 ≈ 2.56509966 \int_0^{2\pi} \frac{d \theta}{2+\cos^2 \theta}=\frac{2\pi}{\sqrt 6}\approx 2.56509966 ∫02π2+cos2θdθ=62π≈2.56509966
我们用Gauss-Legendre Quadrature验证一下这个结果。
对其积分域做变换，
∫
0
2
π
d
θ
2
+
cos
2
θ
=
π
∫
−
1
1
1
2
+
cos
2
(
π
x
+
π
)
⏟
被
积
函
数
d
x
\int_0^{2\pi} \frac{d \theta}{2+\cos^2 \theta}=\pi \int_{-1}^{1} \underbrace{\frac{1}{2+\cos^2(\pi x+\pi)}}_{被积函数}dx
∫02π2+cos2θdθ=π∫−11被积函数
2+cos2(πx+π)1dx
第一步 调用gaussquad包
library(gaussquad)
第二步 创建被积函数
myfunc <- function(x){
A = cos(pi*x+pi);Value = 1/(2+A^2)
return(Value)}
第三步 创建 n n n阶Legendre求积规则的Data frame（包含我们求积分需要的权重与零点）
n <- 20
rules <- legendre.quadrature.rules(n)
因为我选择的是用前20阶的Legendre多项式进行近似，所以rules是一个包含20个元素的list，第
i
i
i个元素是一个data frame，包含零点和权重，比如查看第2个data frame
> rules[[2]]
x w
1 0.5773503 1
2 -0.5773503 1
返回的信息中x列代表2阶Legendre多项式的两个零点，w列表示在Gauss-Legendre Quadrature公式需要的权重。
第四步 计算数值积分
legendre.inte = legendre.quadrature(myfunc,rules[[20]])
因为我们要用20阶的Legendre多项式展开被积函数进行近似，所以输入参数的时候选择第20个data frame。验证一下数值积分的结果
> pi*legendre.inte
[1] 2.565099
发现和理论值确实一致。另外，因为已经有零点和权重的信息了，我们也可以自己写一行代码计算积分
> pi*sum(myfunc(rules[[40]]$x)*(rules[[40]]$w))
[1] 2.5651
也就是把零点处被积函数的取值与权重相乘再加总，可以发现结果与调用函数的结果一致。
例2：Gauss-Laguerre Quadrature
形如 ∫ 0 + ∞ f ( x ) e − x d x \int_{0}^{+\infty} f(x)e^{-x}dx ∫0+∞f(x)e−xdx
的积分可以用Gauss-Laguerre Quadrature进行数值求解，用
L
n
L_n
Ln表示Laguerre多项式，Gauss-Laguerre Quadrature的公式为
∫
0
+
∞
f
(
x
)
d
x
≈
∑
i
=
1
n
w
i
f
(
x
i
)
=
∑
i
=
1
n
x
i
(
n
+
1
)
2
[
L
n
+
1
(
x
i
)
]
2
\int_{0}^{+\infty} f(x)dx\approx \sum_{i=1}^n w_if(x_i) = \sum_{i=1}^n \frac{x_i}{(n+1)^2[L_{n+1}(x_i)]^2}
∫0+∞f(x)dx≈i=1∑nwif(xi)=i=1∑n(n+1)2[Ln+1(xi)]2xi
其中 x i x_i xi代表 L n L_n Ln的第 i i i个零点。
考虑积分
∫
0
+
∞
1
x
(
1
+
x
2
)
e
−
1
2
x
2
d
x
\int_0^{+\infty} \frac{1}{x(1+x^2)}e^{-\frac{1}{2x^2}}dx
∫0+∞x(1+x2)1e−2x21dx
做变量代换 y = 1 2 x 2 y=\frac{1}{2x^2} y=2x21，则 x = 1 2 y d y = − 1 x 3 d x ⇒ d x = − x 3 d y = − 2 − 3 / 2 y − 3 / 2 d y x =\frac{1}{\sqrt{2y}} \\ dy=-\frac{1}{x^3}dx \Rightarrow dx=-x^3dy=-2^{-3/2}y^{-3/2}dy x=2y1dy=−x31dx⇒dx=−x3dy=−2−3/2y−3/2dy
所以
∫
0
+
∞
1
x
(
1
+
x
2
)
e
−
1
2
x
2
d
x
=
∫
0
+
∞
1
2
y
+
1
e
−
y
d
y
=
e
E
1
(
0.5
)
2
≈
0.46145532
\int_0^{+\infty} \frac{1}{x(1+x^2)}e^{-\frac{1}{2x^2}}dx = \int_0^{+\infty} \frac{1}{2y+1}e^{-y}dy=\frac{\sqrt e E_1(0.5)}{2} \approx 0.46145532
∫0+∞x(1+x2)1e−2x21dx=∫0+∞2y+11e−ydy=2eE1(0.5)≈0.46145532
其中 E 1 E_1 E1是exponential integral，一个特殊函数。
myfunc <- function(x){
A = 1/(2*x+1)
return(A)
}
n <- 20
rules <- glaguerre.quadrature.rules(n,0)
laguerre.inte = glaguerre.quadrature(myfunc,rules[[20]])
Gauss-Laguerre Quadrature的glaguerre.quadrature.rules函数比Gauss-Legendre Quadrature的legendre.quadrature.rules函数多一个参数，也就是上面代码中n后面那个参数，因为glaguerre.quadrature.rules函数默认使用的是广义Laguerre多项式，所以多出来这个参数其实是广义Laguerre多项式的参数，比如1阶Laguerre多项式为
−
x
+
1
-x+1
−x+1，一阶广义Laguerre多项式为
−
x
+
α
+
1
-x+\alpha+1
−x+α+1。通常我们不需要用广义Laguerre多项式，所以把那个参数设为0即可。
> laguerre.inte
[1] 0.4614418
可以发现，数值积分的结果与理论结果一致。
例3：Gauss-Hermite Quadrature
形如 ∫ − ∞ + ∞ f ( x ) e − x 2 d x \int_{-\infty}^{+\infty} f(x)e^{-x^2}dx ∫−∞+∞f(x)e−x2dx
的积分可以用Gauss-Hermite Quadrature进行数值求解，用
H
n
H_n
Hn表示Hermite多项式，Gauss-Hermite Quadrature的公式为
∫
0
+
∞
f
(
x
)
d
x
≈
∑
i
=
1
n
w
i
f
(
x
i
)
=
∑
i
=
1
n
2
n
n
!
π
f
(
x
i
)
n
2
[
H
n
−
1
(
x
i
)
]
2
\int_{0}^{+\infty} f(x)dx\approx \sum_{i=1}^n w_if(x_i) = \sum_{i=1}^n \frac{2^n n!\sqrt{\pi}f(x_i)}{n^2[H_{n-1}(x_i)]^2}
∫0+∞f(x)dx≈i=1∑nwif(xi)=i=1∑nn2[Hn−1(xi)]22nn!πf(xi)
其中 x i x_i xi代表 H n H_n Hn的第 i i i个零点。
考虑积分
∫
−
∞
+
∞
1
1
+
x
2
e
−
x
2
2
d
x
=
2
∫
−
∞
+
∞
1
1
+
2
y
2
e
−
y
2
d
y
\int_{-\infty}^{+\infty} \frac{1}{1+x^2}e^{-\frac{x^2}{2}}dx=\sqrt 2 \int_{-\infty}^{+\infty} \frac{1}{1+2y^2}e^{-y^2}dy
∫−∞+∞1+x21e−2x2dx=2∫−∞+∞1+2y21e−y2dy
myfunc <- function(x){
A = 1/(2*x^2+1)
return(A)
}
n <- 20
rules <- hermite.h.quadrature.rules(n)
hermite.inte = hermite.h.quadrature(myfunc,rules[[20]])
结果为
> hermite.inte
[1] 1.161323
