高斯-勒让德积分求解函数积分
前言
梯度和辛普森是经典的几何求积分方法，简单易懂，那如果要更加高档（复杂难懂）的求积分方法找哪家了？高斯-勒让德积分当仁不让。举例来说，下面这个公式看着很高档，但真的要用C来实现还真的有些令人头痛。
A
=
1
2
∫
v
∞
e
−
t
t
d
t
\begin{aligned} A=\frac{1}{2}\int_{v}^{\infty}\frac{e^{-t}}{t} {\rm d}t \end{aligned}
A=21∫v∞te−tdt
高斯-勒让德积分
关于勒让德复杂的表征方法暂时不想多做分析学习，暂且记住勒让德多项式
P
l
(
x
)
，
l
∈
(
0
,
1
,
.
.
.
)
P_l(x)，l\in(0,1,...)
Pl(x)，l∈(0,1,...)为本征函数族，正交和完备性可以令其作为广义傅里叶级数的基，即假设函数
f
(
x
)
f(x)
f(x)定义在区间
(
−
1
,
1
)
(-1,1)
(−1,1)上，则：
f
(
x
)
=
∑
l
=
0
∞
f
l
P
l
(
x
)
f(x)=\sum_{l=0}^{\infty}f_lP_l(x)
f(x)=l=0∑∞flPl(x)系数
f
l
f_l
fl为
f
l
=
2
l
+
1
2
∫
−
1
1
f
(
x
)
P
l
(
x
)
d
x
f_l=\frac{2l+1}{2}\int_{-1}^{1}f(x)P_l(x)dx
fl=22l+1∫−11f(x)Pl(x)dx
而勒让德多项式的微分表征：
P
n
(
x
)
=
1
2
n
n
!
d
n
d
x
n
(
x
2
−
1
)
n
P_n(x)=\frac{1}{2^nn!}\frac{d^n}{dx^n}(x^2-1)^n
Pn(x)=2nn!1dxndn(x2−1)n
看着有点晕，如果展开不同阶数
n
n
n的多项式，可能会直观一些【5】
另外【5】还将5阶的多项式画在同一个图上，虽然看不懂，但也挺震撼。
高斯积分（Gaussion Quadrature）公式如下：
∫
−
1
1
f
(
x
)
d
x
≈
∑
i
=
1
n
w
i
f
(
x
i
)
\int_{-1}^{1}f(x)dx \approx \sum_{i=1}^n w_if(x_i)
∫−11f(x)dx≈i=1∑nwif(xi)
即选择一定的插值点和插值权重去逼近积分值。高斯勒让德给出的权重值为【7】
w
i
=
2
(
1
−
x
2
)
[
P
n
′
(
x
i
)
]
2
w_i=\frac{2}{(1-x^2)[P_n^\prime(x_i)]^2}
wi=(1−x2)[Pn′(xi)]22
而对应的插值点和低阶权重值已经被算好，可以直接拿来用。
一般积分区间的归一化
翻译一下wiki的公式：
∫
a
b
f
(
x
)
d
x
=
∫
−
1
1
f
(
b
−
a
2
ξ
+
a
+
b
2
)
d
x
d
ξ
d
ξ
≈
b
−
a
2
∑
i
=
1
n
w
i
f
(
b
−
a
2
ξ
+
a
+
b
2
)
,
d
x
d
ξ
=
b
−
a
2
\int_a^bf(x)dx = \int_{-1}^{1}f(\frac{b-a}{2}\xi+\frac{a+b}{2})\frac{dx}{d\xi}d\xi\\ \approx\frac{b-a}{2}\sum_{i=1}^nw_if(\frac{b-a}{2}\xi+\frac{a+b}{2}),\ \frac{dx}{d\xi}=\frac{b-a}{2}
∫abf(x)dx=∫−11f(2b−aξ+2a+b)dξdxdξ≈2b−ai=1∑nwif(2b−aξ+2a+b), dξdx=2b−a
Exponential Integral
文章开头的指数积分，也有比较成熟的计算方法：
E
1
(
z
)
=
∫
z
∞
e
−
t
t
d
t
=
−
γ
−
l
n
z
−
∑
k
=
1
∞
(
−
z
)
k
k
k
!
E_1(z)=\int_z^\infty\frac{e^{-t}}{t}dt=-\gamma-ln\ z-\sum_{k=1}^\infty\frac{(-z)^k}{kk!}
E1(z)=∫z∞te−tdt=−γ−ln z−k=1∑∞kk!(−z)k
更通用的表达，定义
E
i
n
(
z
)
=
∫
0
z
(
1
−
e
−
t
)
d
t
t
=
∑
k
=
1
∞
(
−
1
)
k
+
1
z
k
k
k
!
Ein(z)=\int_0^z(1-e^{-t})\frac{dt}{t}=\sum_{k=1}^{\infty}\frac{(-1)^{k+1}z^k}{kk!}
Ein(z)=∫0z(1−e−t)tdt=k=1∑∞kk!(−1)k+1zk
那么
E
1
(
z
)
=
−
γ
−
l
n
z
+
E
i
n
(
z
)
E_1(z)=-\gamma-ln\ z+Ein(z)
E1(z)=−γ−ln z+Ein(z)
第一个参数叫做欧拉常数，从wiki上抄写推导：
γ
=
lim
n
→
∞
(
−
l
o
g
n
+
∑
k
=
1
n
1
k
)
≈
0.577215664901532860606512
\gamma=\lim_{n \to \infty}(-log\ n +\sum_{k=1}^n\frac{1}{k})\approx0.577215664901532860606512
γ=n→∞lim(−log n+k=1∑nk1)≈0.577215664901532860606512
实验
最后用【10】里提供的算法，计算 − E 1 ( z ) -E_1(z) −E1(z)的几个积分值：
x=0.05 Ei(x)=-2.4678985
x=0.25 Ei(x)=-1.0442826
x=0.45 Ei(x)=-0.6253313
x=0.65 Ei(x)=-0.4115170
x=0.85 Ei(x)=-0.2840193
x=1.05 Ei(x)=-0.2018728
x=1.25 Ei(x)=-0.1464134
x=1.45 Ei(x)=-0.1077774
x=1.65 Ei(x)=-0.0802476
x=1.85 Ei(x)=-0.0602950
利用python的scipy.specials 库exp1计算的 E 1 ( z ) E_1(z) E1(z)对应值：
>>> exp1(0.05) 2.467898488509974
>>> exp1(0.25) 1.0442826344437384
>>> exp1(0.45) 0.6253313163232692
>>> exp1(0.65) 0.41151697594947967
>>> exp1(0.85) 0.28401926850246173
>>> exp1(1.05) 0.2018728132201966
>>> exp1(1.25) 0.14641337252591016
>>> exp1(1.45) 0.10777743986897659
>>> exp1(1.65) 0.08024762667334323
>>> exp1(1.85) 0.06029496682837346
参考
1.勒让德多项式及性质
2.高斯-勒让德数值积分公式
3.高斯-勒让德求积公式及Matlab实现
4.高斯积分
5.Legendre polynomials from wikipedia
6.6.2(i) Exponential and Logarithmic Integrals
7.Legendre-Gauss Quadrature
8.Exponential Integral from mathworld
9.Exponential integral from wikepedia
10.常用算法程序集(C语言描述)第三版
9933
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
