import numpy as np
from scipy. integrate import quad
from sympy import *
init_printing( )
import matplotlib. pyplot as plt
数值积分
数值积分的本质都是采用了微积分的思想:化整为零(将积分区间离散化),以直代曲(基于积分中值定理,用函数值的加权平均数近似表示区间的平均高度)
闭型牛顿——科斯特公式:梯形公式、辛普森(simpson)公式和布尔(Boole)公式,都是选择等距节点 ;
高斯积分公式:高斯——勒让德、高斯——切比雪夫、高斯——拉盖尔和高斯——埃尔米特。高斯——勒让德(Gauss-Legendre)公式选择某些勒让德多项式的根作为不等距节点 。类似的还有高斯——切比雪夫积分、高斯——拉盖尔积分和高斯——埃尔米特积分等等。
积分中值定理——平均高度的各种近似法
共同点:
(a) 都依赖于步长(梯形宽度),权重一般呈中间高,两边低的对称分布。 (b) 用基于各个节点
x
0
,
x
1
,
.
.
.
,
x
M
x_0, x_1,...,x_M
x 0 , x 1 , . . . , x M 的
f
(
x
)
f(x)
f ( x ) 的拉格朗日逼近多项式
P
M
(
x
)
P_M(x)
P M ( x ) 来代替被积函数
f
(
x
)
f(x)
f ( x ) ,本质上就是取各个节点的"加权平均数",这里的权重之和不一定为1。
不同点
(a) 闭型牛顿——科斯特公式选择等距节点,这限制了求积公式的代数精度;而高斯——勒让德(Gauss-Legendre)公式则取消了这个限制条件,因此能够极大地提高了代数精度。从数值分析可以知道, 高斯方法的确比梯形法, 辛普森, 龙格库塔法有优越性: 对被积函数的拟合精度高, 收敛快。
1.闭型牛顿——科斯特公式
假设积分区间固定为
[
a
,
b
]
[a,b]
[ a , b ] ,那么不同公式要使用不一样的步长
h
h
h 。
1.梯形公式:2个点,将积分区间分为
1
1
1 等分,精度为
n
=
1
n=1
n = 1 ,步长
h
=
b
−
a
h=b-a
h = b − a
2.辛普森公式:3个点,将积分区间分为
2
2
2 等分,精度为
n
=
3
n=3
n = 3 ,步长
h
=
b
−
a
2
h=\frac{b-a}{2}
h = 2 b − a
3.辛普森
3
8
\frac{3}{8}
8 3 公式:4个点,将积分区间分为
3
3
3 等分,精度为
n
=
3
n=3
n = 3 ,步长
h
=
b
−
a
3
h=\frac{b-a}{3}
h = 3 b − a
4.布尔公式:5个点,将积分区间分为
4
4
4 等分,精度为
n
=
5
n=5
n = 5 , 步长
h
=
b
−
a
4
h=\frac{b-a}{4}
h = 4 b − a
5.组合梯形:将积分区间分为
M
M
M 等分,
M
M
M 为正整数,共
(
M
+
1
)
(M+1)
( M + 1 ) 个点,步长为
h
=
b
−
a
M
h=\frac{b-a}{M}
h = M b − a ,误差的阶为
O
(
h
2
)
O(h^{2})
O ( h 2 )
6.组合辛普森公式:将积分区间分为
2
M
2M
2 M 等分,共
(
2
M
+
1
)
(2M+1)
( 2 M + 1 ) 个点,步长为
h
=
b
−
a
2
M
h=\frac{b-a}{2M}
h = 2 M b − a ,误差的阶为
O
(
h
4
)
O(h^{4})
O ( h 4 )
各公式的具体表达式如下:
Python 3.6 代码
def trapezoid ( h , f) :
Integral_appro = ( h / 2 ) * ( f[ 0 ] + f[ 1 ] )
return Integral_appro
def simpson ( h, f) :
Integral_appro = ( h / 3 ) * ( f[ 0 ] + 4 * f[ 1 ] + f[ 2 ] )
return Integral_appro
def simpson_3_8 ( h, f) :
Integral_appro = ( 3 / 8 * h) * ( f[ 0 ] + 3 * f[ 1 ] + 3 * f[ 2 ] + f[ 3 ] )
return Integral_appro
def boole ( h, f) :
Integral_appro = ( 2 / 45 * h) * ( 7 * f[ 0 ] + 32 * f[ 1 ] + 12 * f[ 2 ] + 32 * f[ 3 ] + 7 * f[ 4 ] )
return Integral_appro
闭型牛顿——科斯特公式假设积分区间为[0,1],那么各公式的步长 h 依次为1,1/2, 1/3, 1/4, 也就是分别将区间等分为1、2、3、4份,取各个节点的函数值
f
(
x
i
)
f(x_i)
f ( x i ) 的"加权平均数"
w
(
x
i
)
f
(
x
i
)
w(x_i)f(x_i)
w ( x i ) f ( x i ) 作为整个区间平均高度
f
(
e
)
f(e)
f ( e ) 的近似
def func ( x) :
f = 1 + np. exp( - x) * np. sin( 4 * x)
return f
I_precise = quad( func, 0 , 1 )
I_precise
(
1.3082506046426685
,
1.4524499432540732
e
−
14
)
\left ( 1.3082506046426685, \quad 1.4524499432540732e-14\right )
( 1 . 3 0 8 2 5 0 6 0 4 6 4 2 6 6 8 5 , 1 . 4 5 2 4 4 9 9 4 3 2 5 4 0 7 3 2 e − 1 4 )
x_s = Symbol( 'x_s' , real= True )
f_s = 1 + exp( - x_s) * sin( 4 * x_s)
f_s
1
+
e
−
x
s
sin
(
4
x
s
)
1 + e^{- x_{s}} \sin{\left (4 x_{s} \right )}
1 + e − x s sin ( 4 x s )
Inte_f_symbol = Integral( f_s, ( x_s, 0 , 1 ) ) # 积分表达式(只是符号,没有求出原函数)
Inte_f_s = integrate( f_s, x_s) # 积分表达式(求出原函数)
Inte_f_s_n1 = integrate( f_s, ( x_s, 0 , 1 ) ) # 定积分
Inte_f_s_n2 = integrate( f_s, ( x_s, 0 , 1 ) ) . evalf( ) # 定积分,并将结果化成小数
Inte_f_symbol, Inte_f_s, Inte_f_s_n1, Inte_f_s_n2
(
∫
0
1
(
1
+
e
−
x
s
sin
(
4
x
s
)
)
 
d
x
s
,
x
s
−
e
−
x
s
17
sin
(
4
x
s
)
−
4
17
e
−
x
s
cos
(
4
x
s
)
,
−
sin
(
4
)
17
e
−
4
17
e
cos
(
4
)
+
21
17
,
1.30825060464267
)
\left ( \int_{0}^{1} \left(1 + e^{- x_{s}} \sin{\left (4 x_{s} \right )}\right)\, dx_{s}, \quad x_{s} - \frac{e^{- x_{s}}}{17} \sin{\left (4 x_{s} \right )} - \frac{4}{17} e^{- x_{s}} \cos{\left (4 x_{s} \right )}, \quad - \frac{\sin{\left (4 \right )}}{17 e} - \frac{4}{17 e} \cos{\left (4 \right )} + \frac{21}{17}, \quad 1.30825060464267\right )
( ∫ 0 1 ( 1 + e − x s sin ( 4 x s ) ) d x s , x s − 1 7 e − x s sin ( 4 x s ) − 1 7 4 e − x s cos ( 4 x s ) , − 1 7 e sin ( 4 ) − 1 7 e 4 cos ( 4 ) + 1 7 2 1 , 1 . 3 0 8 2 5 0 6 0 4 6 4 2 6 7 )
# 假设积分区间为[0,1],那么各公式的步长 h 依次为1,1/2, 1/3, 1/4,分别将区间分为1、2、3、4份,
# 取各个节点的函数值的加权平均数作为整个区间平均高度的近似
I_trapezoid = trapezoid( 1 , [ func( 0 ) , func( 1 ) ] )
I_simpson = simpson( 1 / 2 , [ func( 0 ) , func( 1 / 2 ) , func( 1 ) ] )
I_simpson_3_8 = simpson_3_8( 1 / 3 , [ func( 0 ) , func( 1 / 3 ) , func( 2 / 3 ) , func( 1 ) ] )
I_boole = boole( 1 / 4 , [ func( 0 ) , func( 1 / 4 ) , func( 2 / 4 ) , func( 3 / 4 ) , func( 1 ) ] )
I_trapezoid, I_simpson, I_simpson_3_8, I_boole
(
0.8607939604744832
,
1.3212758322698814
,
1.3143968149336276
,
1.3085919215646966
)
\left ( 0.8607939604744832, \quad 1.3212758322698814, \quad 1.3143968149336276, \quad 1.3085919215646966\right )
( 0 . 8 6 0 7 9 3 9 6 0 4 7 4 4 8 3 2 , 1 . 3 2 1 2 7 5 8 3 2 2 6 9 8 8 1 4 , 1 . 3 1 4 3 9 6 8 1 4 9 3 3 6 2 7 6 , 1 . 3 0 8 5 9 1 9 2 1 5 6 4 6 9 6 6 )
def func ( x) :
f = 2 + np. sin( 2 * np. sqrt( x) )
return f
# M 等分
def composite_trapezoid ( f, a, b, M) :
if M < 1 :
print ( 'M must be larger than or equal to 1' )
return
h = ( b - a) / M
s1 = h / 2 * ( f( a) + f( b) )
s2 = 0
for k in range ( 1 , M) :
x = a + h * k
s2 += f( x)
s = s1 + s2 * h
return s
# 2M 等分
def composite_simpson ( f, a, b, M) :
if M < 1 :
print ( 'M must be larger than or equal to 1' )
return
h = ( b - a) / ( 2 * M)
s1 = h / 3 * ( f( a) + f( b) )
s2 = 0
# k = 1, 2,...,M
for k in range ( 1 , M+ 1 ) :
x = a + h * ( 2 * k - 1 )
s2 += f( x)
s2 = 4 * h / 3 * s2
s3 = 0
# k = 1,2,...,M-1
for k in range ( 1 , M) :
x = a + h * k * 2
s3 += f( x)
s3 = 2 * h / 3 * s3
s = s1 + s2 + s3
return s
I_true = quad( func, 1 , 6 )
I_tra_M_10 = composite_trapezoid( func, 1 , 6 , 10 )
I_tra_M_20 = composite_trapezoid( func, 1 , 6 , 20 )
I_tra_M_40 = composite_trapezoid( func, 1 , 6 , 40 )
I_tra_M_60 = composite_trapezoid( func, 1 , 6 , 60 )
I_tra_M_80 = composite_trapezoid( func, 1 , 6 , 80 )
I_tra_M_100 = composite_trapezoid( func, 1 , 6 , 100 )
I_true, I_tra_M_10, I_tra_M_20, I_tra_M_40, I_tra_M_60, I_tra_M_80, I_tra_M_100
(
(
8.183479207662728
,
1.79863015057564
e
−
10
)
,
8.193854565172531
,
8.186049263770313
,
8.184120191790313
,
8.183763962635512
,
8.183639357318619
,
8.183581696027387
)
\left ( \left ( 8.183479207662728, \quad 1.79863015057564e-10\right ), \quad 8.193854565172531, \quad 8.186049263770313, \quad 8.184120191790313, \quad 8.183763962635512, \quad 8.183639357318619, \quad 8.183581696027387\right )
( ( 8 . 1 8 3 4 7 9 2 0 7 6 6 2 7 2 8 , 1 . 7 9 8 6 3 0 1 5 0 5 7 5 6 4 e − 1 0 ) , 8 . 1 9 3 8 5 4 5 6 5 1 7 2 5 3 1 , 8 . 1 8 6 0 4 9 2 6 3 7 7 0 3 1 3 , 8 . 1 8 4 1 2 0 1 9 1 7 9 0 3 1 3 , 8 . 1 8 3 7 6 3 9 6 2 6 3 5 5 1 2 , 8 . 1 8 3 6 3 9 3 5 7 3 1 8 6 1 9 , 8 . 1 8 3 5 8 1 6 9 6 0 2 7 3 8 7 )
I_sim_M_10 = composite_simpson( func, 1 , 6 , 10 )
I_sim_M_20 = composite_simpson( func, 1 , 6 , 20 )
I_sim_M_10, I_sim_M_20
(
8.183447496636239
,
8.183477167796982
)
\left ( 8.183447496636239, \quad 8.183477167796982\right )
( 8 . 1 8 3 4 4 7 4 9 6 6 3 6 2 3 9 , 8 . 1 8 3 4 7 7 1 6 7 7 9 6 9 8 2 )
2. 高斯——勒让德公式
2.1 高斯积分
高斯积分法是精度最高的插值型数值积分,具有
2
n
−
1
2n-1
2 n − 1 阶精度,并且高斯积分总是稳定。而高斯求积系数,可以由Lagrange多项式插值系数进行积分得到。
高斯积分有很多种,高斯——勒让德只是其中常用的一种。如果对具体的推导过程及其它类型的Gauss积分有兴趣,请查阅参考资料 [4] 高斯——勒让德公式–中南大学 和 [5] 各种类型的高斯积分 ,讲得很详细!
如果你想更进一步,请看这里!
目前最高效的数值积分方法:Gauss–Kronrod quadrature formula及其变种,这是一种自适应的数值积分方法,在高斯——勒让德的
n
n
n 个节点中,根据Stieltjes多项式增加了
n
+
1
n+1
n + 1 个节点,因此一共有
2
n
+
1
2n+1
2 n + 1 个节点,特别适合于高度振荡的函数。 详情请看[3] wiki:Gauss–Kronrod quadrature formula ,目前最高效的实现为Laurie(1997)的论文[8] 及后来Calvetti 等人(2000)提出的改进版本[9](在百度学术上都能找到可免费下载的资源)
Virginia Tech 的 John Burkardt 老师在2016年给出了 Gauss–Kronrod 积分的 Python 3实现:http://people.sc.fsu.edu/~jburkardt/py_src/kronrod/kronrod.html
下面摘录Calvetti 等人(2000)的文章中一些关键的定义:
好啦,来聊些简单的东西。 首先看看高斯积分的定义:
注意:这张PPT是我从参考资料中随手截取的,其实它取了
n
+
1
n+1
n + 1 个点,因此精度为
2
n
+
1
2n+1
2 n + 1 ,而下面的讨论都是取
n
n
n 个点,别蒙圈了兄弟
注:关于正交多项式的定义,在一个国外的课件找到另外一种定义,暂时还没搞清楚哪一种更准确,先记录一下:
感叹一下:百度文库进步了,连文献和国外的PPT都有!
https://wenku.baidu.com/view/1d88b420bcd126fff7050b0c.html?rec_flag=default
不同的权函数
ρ
(
x
)
\rho(x)
ρ ( x ) 对应不同的正交多项式
P
(
x
)
P(x)
P ( x ) 。
高斯——勒让德积分的权函数
ρ
(
x
)
=
1
\rho(x)=1
ρ ( x ) = 1 ,对应的正交多项式
P
(
x
)
P(x)
P ( x ) 为勒让德多项式
P
n
(
x
)
P_n(x)
P n ( x ) ,积分区间为
[
−
1
,
1
]
[-1,1]
[ − 1 , 1 ] 。
高斯——勒让德积分的本质就是被积函数
f
(
x
)
f(x)
f ( x ) 在区间
(
−
1
,
1
)
(-1,1)
( − 1 , 1 ) 内的积分近似于区间
(
−
1
,
1
)
(-1,1)
( − 1 , 1 ) 内的
n
n
n 个勒让德——多项式零点的函数值
f
(
x
k
)
f(x_k)
f ( x k ) 的加权平均数
s
u
m
(
w
k
f
(
x
k
)
)
sum(w_k f(x_k))
s u m ( w k f ( x k ) ) ,权重之和不一定为1,权重依然呈现中间高,两边低的趋势。
高斯——切比雪夫积分的权函数
ρ
(
x
)
=
1
1
−
x
2
\rho(x)=\frac{1}{1-x^2}
ρ ( x ) = 1 − x 2 1 ,对应的正交多项式
P
(
x
)
P(x)
P ( x ) 为切比雪夫多项式
T
n
(
x
)
T_n(x)
T n ( x ) ,积分区间为
[
−
1
,
1
]
[-1,1]
[ − 1 , 1 ] 。**
高斯——拉盖尔积分的权函数
ρ
(
x
)
=
e
−
x
\rho(x)=e^{-x}
ρ ( x ) = e − x ,对应的正交多项式
P
(
x
)
P(x)
P ( x ) 为拉盖尔多项式
L
n
(
x
)
L_n(x)
L n ( x ) ,积分区间为
[
0
,
+
∞
]
[0,+\infty]
[ 0 , + ∞ ] 。
2.1.4 高斯——埃尔米特积分 高斯——埃尔米特积分的权函数
ρ
(
x
)
=
e
−
x
2
\rho(x)=e^{-x^{2}}
ρ ( x ) = e − x 2 ,对应的正交多项式
P
(
x
)
P(x)
P ( x ) 为埃尔米特多项式
H
n
(
x
)
H_n(x)
H n ( x ) ,积分区间为
[
−
∞
,
+
∞
]
[-\infty,+\infty]
[ − ∞ , + ∞ ] 。
numpy有对应的模块:
Gauss-Legendre: 节点和权值:numpy.polynomial.legendre.leggauss 权函数
ρ
(
x
)
=
1
\rho(x)=1
ρ ( x ) = 1 :numpy.polynomial.legendre.legweight
Gauss-Chebyshev: 节点和权值:numpy.polynomial.chebyshev.chebgauss 权函数
ρ
(
x
)
=
1
1
−
x
2
\rho(x)=\frac{1}{1-x^2}
ρ ( x ) = 1 − x 2 1 :numpy.polynomial.chebyshev.chebweight
Gauss-Laguerre: 节点和权值:numpy.polynomial.laguerre.laggauss 权函数
ρ
(
x
)
=
e
−
x
\rho(x)=e^{-x}
ρ ( x ) = e − x :numpy.polynomial.laguerre.lagweight
Gauss-Hermite: 节点和权值:numpy.polynomial.hermite.hermgauss 权函数
ρ
(
x
)
=
e
−
x
2
\rho(x)=e^{-x^{2}}
ρ ( x ) = e − x 2 :numpy.polynomial.hermite.hermweight
对于任意有限区间
[
a
,
b
]
,
a
,
b
≠
±
∞
[a,b], a,b≠±\infty
[ a , b ] , a , b ̸ = ± ∞ 内的积分,可以借助线性变换将积分区间变换到
[
−
1
,
1
]
[-1,1]
[ − 1 , 1 ] ,从而利用高斯——勒让德积分和高斯——切比雪夫积分进行计算。
2.2 高斯-勒让德积分
高斯-勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定,让函数
f
(
x
)
f(x)
f ( x ) 依次取
i
=
0
,
1
,
2....
n
i=0,1,2....n
i = 0 , 1 , 2 . . . . n 次多项式,使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是
2
n
−
1
2n−1
2 n − 1 ,而且是最高的。通常运用的是
(
−
1
,
1
)
(−1,1)
( − 1 , 1 ) 的积分节点和积分系数,其他积分域是通过一定的变换变换到-1到1之间积分。
勒让德多项式 首先复习一下勒让德多项式:
一般
N
N
N 点的高斯——勒让德公式对于次数小于等于
2
N
−
1
2N-1
2 N − 1 次的多项式是精确的,即
G
N
(
f
)
G_N(f)
G N ( f ) 精度为
n
=
2
N
−
1
n=2N-1
n = 2 N − 1 ,误差项
E
[
f
]
E[f]
E [ f ] 包含
2
N
2N
2 N 阶导数。
Q:为什么
N
N
N 点的高斯——勒让德公式对于次数小于等于
2
N
−
1
2N-1
2 N − 1 次的多项式是精确的呢?
A:因为
N
N
N 点的高斯——勒让德公式包含
N
N
N 个未知节点和相应的
N
N
N 个未知权值,一个有
2
N
2N
2 N 个未知数。我们可以从0阶多项式开始构造方程,则
2
N
2N
2 N 个未知数需要
2
N
2N
2 N 个等式方程,需要多项式
x
c
,
c
=
0
,
1
,
2
,
.
.
.
2
N
−
1
x^c, c=0,1,2,...2N-1
x c , c = 0 , 1 , 2 , . . . 2 N − 1 ,所以对次数小于等于
2
N
−
1
2N-1
2 N − 1 次的多项式是可以取等号的,因此是精确的,而对
2
N
−
1
2N-1
2 N − 1 次及以上的多项式则只能取约等于号,存在误差。
因此
2
2
2 点的高斯——勒让德多项式对于
3
3
3 次的多项式是精确的,
3
3
3 点的高斯——勒让德多项式对于
5
5
5 次的多项式是精确的。
因此,
2
2
2 点的高斯——勒让德多项式的误差项
E
[
f
]
E[f]
E [ f ] 包含
4
4
4 阶导数,
3
3
3 点的高斯——勒让德多项式的误差项
E
[
f
]
E[f]
E [ f ] 包含
6
6
6 阶导数。
N
N
N 点的高斯——勒让德公式的节点就是
n
n
n 阶勒让德多项式
P
n
(
x
)
P_n(x)
P n ( x ) 的根。
想要在区间
t
∈
[
a
,
b
]
t∈[a,b]
t ∈ [ a , b ] 上使用高斯——勒让德公式,只需作变量替换
t
=
(
a
+
b
)
/
2
+
(
b
−
a
)
/
2
∗
x
t=(a+b)/2 + (b-a)/2 * x
t = ( a + b ) / 2 + ( b − a ) / 2 ∗ x ,节点和权值必须从高斯——勒让德节点与权值表中获得。
以下是Python3.6的代码实现。
2.1 用Sympy和numpy求勒让德多项式的节点
n阶勒让德多项式
x = Symbol( 'x' , real= True )
n = Symbol( 'n' , real= True )
legendre( n, x)
P
n
(
x
)
P_{n}\left(x\right)
P n ( x )
0——4阶勒让德多项式及多项式的根
legendre( 0 , x) , solve( legendre( 0 , x) , x) # 0阶
(
1
,
[
]
)
\left ( 1, \quad \left [ \right ]\right )
( 1 , [ ] )
legendre( 1 , x) , solve( legendre( 1 , x) , x) # 1阶
(
x
,
[
0
]
)
\left ( x, \quad \left [ 0\right ]\right )
( x , [ 0 ] )
legendre( 2 , x) , solve( legendre( 2 , x) , x) # 2阶
(
3
x
2
2
−
1
2
,
[
−
3
3
,
3
3
]
)
\left ( \frac{3 x^{2}}{2} - \frac{1}{2}, \quad \left [ - \frac{\sqrt{3}}{3}, \quad \frac{\sqrt{3}}{3}\right ]\right )
( 2 3 x 2 − 2 1 , [ − 3 3
, 3 3
] )
legendre( 3 , x) , solve( legendre( 3 , x) , x) # 3阶
(
5
x
3
2
−
3
x
2
,
[
0
,
−
15
5
,
15
5
]
)
\left ( \frac{5 x^{3}}{2} - \frac{3 x}{2}, \quad \left [ 0, \quad - \frac{\sqrt{15}}{5}, \quad \frac{\sqrt{15}}{5}\right ]\right )
( 2 5 x 3 − 2 3 x , [ 0 , − 5 1 5
, 5 1 5
] )
legendre( 4 , x) , solve( legendre( 4 , x) , x) # 4阶
(
35
x
4
8
−
15
x
2
4
+
3
8
,
[
−
−
2
30
35
+
3
7
,
−
2
30
35
+
3
7
,
−
2
30
35
+
3
7
,
2
30
35
+
3
7
]
)
\left ( \frac{35 x^{4}}{8} - \frac{15 x^{2}}{4} + \frac{3}{8}, \quad \left [ - \sqrt{- \frac{2 \sqrt{30}}{35} + \frac{3}{7}}, \quad \sqrt{- \frac{2 \sqrt{30}}{35} + \frac{3}{7}}, \quad - \sqrt{\frac{2 \sqrt{30}}{35} + \frac{3}{7}}, \quad \sqrt{\frac{2 \sqrt{30}}{35} + \frac{3}{7}}\right ]\right )
⎝ ⎛ 8 3 5 x 4 − 4 1 5 x 2 + 8 3 , ⎣ ⎡ − − 3 5 2 3 0
+ 7 3
, − 3 5 2 3 0
+ 7 3
, − 3 5 2 3 0
+ 7 3
, 3 5 2 3 0
+ 7 3
⎦ ⎤ ⎠ ⎞
求
f
=
1
/
x
f = 1 / x
f = 1 / x 在[1,5]上的积分
def func ( x) :
f = 1 / x
return f
# nodes是标准的勒让德多项式的根
# T是替换变量
# 对于不同阶数的勒让德多项式,nodes和weights是不相同的
def gauss_lengendre ( fun, a, b, nodes, weights) :
nodes = np. array( nodes)
weights = np. array( weights)
T = ( a+ b) / 2 + ( b- a) / 2 * nodes # 变量替换
quad = ( b- a) / 2 * np. sum ( weights * fun( T) ) # element-wise multiplication
return quad
精确值
I_true = quad( func, 1 , 5 )
I_true
(
1.6094379124341014
,
3.6599536780638603
e
−
09
)
\left ( 1.6094379124341014, \quad 3.6599536780638603e-09\right )
( 1 . 6 0 9 4 3 7 9 1 2 4 3 4 1 0 1 4 , 3 . 6 5 9 9 5 3 6 7 8 0 6 3 8 6 0 3 e − 0 9 )
用3点勒让德多项式近似
le_3 = legendre( 3 , x)
nodes_3 = solve( le_3, x)
nodes_3
[
0
,
−
15
5
,
15
5
]
\left [ 0, \quad - \frac{\sqrt{15}}{5}, \quad \frac{\sqrt{15}}{5}\right ]
[ 0 , − 5 1 5
, 5 1 5
]
weights = np. array( [ 5 / 9 , 8 / 9 , 5 / 9 ] )
nodes = np. array( [ - np. sqrt( 15 ) / 5 , 0 , np. sqrt( 15 ) / 5 ] )
a = 1
b = 5
I_n_3 = gauss_lengendre( func, a, b, nodes, weights)
error_n_3 = abs ( I_true[ 0 ] - I_n_3)
I_n_3, error_n_3
(
1.6026936026936027
,
0.006744309740498666
)
\left ( 1.6026936026936027, \quad 0.006744309740498666\right )
( 1 . 6 0 2 6 9 3 6 0 2 6 9 3 6 0 2 7 , 0 . 0 0 6 7 4 4 3 0 9 7 4 0 4 9 8 6 6 6 )
用8点勒让德多项式近似
le_8 = legendre( 8 , x)
nodes_8 = solve( le_8, x) # 求出根节点
nodes_8_f = [ i. evalf( ) for i in nodes_8] # list nodes_8_f 的每个元素都是 sympy.core.numbers.Float
nodes_8_arr = np. array( nodes_8_f, dtype= float ) # 转换成ndarray,将数据类型改为float
nodes_8_arr
array([-0.18343464, 0.18343464, -0.52553241, 0.52553241, -0.79666648,
0.79666648, -0.96028986, 0.96028986])
weights_8_arr = np. array( [ 0.3626837834 , 0.3626837834 , 0.3137066459 , 0.3137066459 ,
0.2223810345 , 0.2223810345 , 0.1012285363 , 0.1012285363 ] , dtype= np. float64)
weights_8_arr
array([0.36268378, 0.36268378, 0.31370665, 0.31370665, 0.22238103,
0.22238103, 0.10122854, 0.10122854])
a = 1
b = 5
I_n_8 = gauss_lengendre( func, a, b, nodes_8_arr, weights_8_arr)
error_n_8 = abs ( I_true[ 0 ] - I_n_8)
I_n_8, error_n_8
(
1.609437439052705
,
4.733813963042621
e
−
07
)
\left ( 1.609437439052705, \quad 4.733813963042621e-07\right )
( 1 . 6 0 9 4 3 7 4 3 9 0 5 2 7 0 5 , 4 . 7 3 3 8 1 3 9 6 3 0 4 2 6 2 1 e − 0 7 )
2.2 用numpy.polynomial.legendre.leggauss(n)一键求出n阶高斯——勒让德节点和权值
numpy.polynomial.legendre.leggauss(n) 函数可以一键求出
n
n
n 阶高斯——勒让德节点和权值,在
n
<
=
100
n<=100
n < = 1 0 0 的范围内,这个函数的结果是准确的。
n_100, w_100 = np. polynomial. legendre. leggauss( 100 )
n_100, w_100
(array([-0.99971373, -0.99849195, -0.99629513, -0.99312494, -0.9889844 ,
-0.98387754, -0.97780936, -0.97078578, -0.96281365, -0.95390078,
-0.94405587, -0.93328854, -0.9216093 , -0.90902957, -0.89556164,
-0.88121868, -0.86601469, -0.84996453, -0.83308388, -0.81538924,
-0.79689789, -0.77762791, -0.75759812, -0.73682809, -0.71533812,
-0.6931492 , -0.67028302, -0.64676191, -0.62260886, -0.59784747,
-0.57250193, -0.54659701, -0.52015802, -0.49321079, -0.46578165,
-0.4378974 , -0.40958529, -0.38087298, -0.35178853, -0.32236034,
-0.29261719, -0.26258812, -0.23230248, -0.20178986, -0.17108008,
-0.14020314, -0.1091892 , -0.07806858, -0.04687168, -0.01562898,
0.01562898, 0.04687168, 0.07806858, 0.1091892 , 0.14020314,
0.17108008, 0.20178986, 0.23230248, 0.26258812, 0.29261719,
0.32236034, 0.35178853, 0.38087298, 0.40958529, 0.4378974 ,
0.46578165, 0.49321079, 0.52015802, 0.54659701, 0.57250193,
0.59784747, 0.62260886, 0.64676191, 0.67028302, 0.6931492 ,
0.71533812, 0.73682809, 0.75759812, 0.77762791, 0.79689789,
0.81538924, 0.83308388, 0.84996453, 0.86601469, 0.88121868,
0.89556164, 0.90902957, 0.9216093 , 0.93328854, 0.94405587,
0.95390078, 0.96281365, 0.97078578, 0.97780936, 0.98387754,
0.9889844 , 0.99312494, 0.99629513, 0.99849195, 0.99971373]),
array([0.00073463, 0.00170939, 0.00268393, 0.00365596, 0.00462445,
0.00558843, 0.00654695, 0.00749907, 0.00844387, 0.00938042,
0.0103078 , 0.01122511, 0.01213146, 0.01302595, 0.01390771,
0.01477588, 0.01562962, 0.01646809, 0.01729046, 0.01809594,
0.01888374, 0.01965309, 0.02040323, 0.02113344, 0.021843 ,
0.02253122, 0.02319742, 0.02384096, 0.0244612 , 0.02505754,
0.0256294 , 0.02617622, 0.02669746, 0.02719261, 0.0276612 ,
0.02810276, 0.02851685, 0.02890309, 0.02926108, 0.02959049,
0.02989098, 0.03016227, 0.03040408, 0.03061619, 0.03079838,
0.03095048, 0.03107234, 0.03116384, 0.03122488, 0.03125542,
0.03125542, 0.03122488, 0.03116384, 0.03107234, 0.03095048,
0.03079838, 0.03061619, 0.03040408, 0.03016227, 0.02989098,
0.02959049, 0.02926108, 0.02890309, 0.02851685, 0.02810276,
0.0276612 , 0.02719261, 0.02669746, 0.02617622, 0.0256294 ,
0.02505754, 0.0244612 , 0.02384096, 0.02319742, 0.02253122,
0.021843 , 0.02113344, 0.02040323, 0.01965309, 0.01888374,
0.01809594, 0.01729046, 0.01646809, 0.01562962, 0.01477588,
0.01390771, 0.01302595, 0.01213146, 0.01122511, 0.0103078 ,
0.00938042, 0.00844387, 0.00749907, 0.00654695, 0.00558843,
0.00462445, 0.00365596, 0.00268393, 0.00170939, 0.00073463]))
I_n_100 = gauss_lengendre( func, a, b, n_100, w_100)
error_n_100 = abs ( I_true[ 0 ] - I_n_100)
I_n_100, error_n_100
(
1.6094379124340965
,
4.884981308350689
e
−
15
)
\left ( 1.6094379124340965, \quad 4.884981308350689e-15\right )
( 1 . 6 0 9 4 3 7 9 1 2 4 3 4 0 9 6 5 , 4 . 8 8 4 9 8 1 3 0 8 3 5 0 6 8 9 e − 1 5 )
error_n_3, error_n_8, error_n_100
(
0.006744309740498666
,
4.733813963042621
e
−
07
,
4.884981308350689
e
−
15
)
\left ( 0.006744309740498666, \quad 4.733813963042621e-07, \quad 4.884981308350689e-15\right )
( 0 . 0 0 6 7 4 4 3 0 9 7 4 0 4 9 8 6 6 6 , 4 . 7 3 3 8 1 3 9 6 3 0 4 2 6 2 1 e − 0 7 , 4 . 8 8 4 9 8 1 3 0 8 3 5 0 6 8 9 e − 1 5 )
比较不同阶数的计算误差,可见高阶的一般能取得较高精度。
[1]高斯——勒让德积分中不同阶数下最大高斯节点间距的关系 节点数N与最大间距呈反相关;帖子给出了2-8阶勒让德多项式的权值 [2] numpy.polynomial.legendre.leggauss [3] wiki:Gauss–Kronrod quadrature formula [4] 高斯——勒让德公式–中南大学 [5] 各种类型的高斯积分 [6] 高斯型求积公式 [7] Calculation of Gauss Quadrature Rules [8] Laurie D P. Calculation of Gauss-Kronrod quadrature rules[J]. Mathematics of Computation, 1997, 66(219):1133-1145 [9] Calvetti D, Golub G H, Gragg W B, et al. Computation of Gauss-Kronrod Quadrature Rules[J]. Mathematics of Computation, 2000, 69(231):1035-1052