7.3 Romberg积分
Romberg求积方法是在积分区间逐次分半的过程中利用外推法产生的一种数值积分方法,当被积函数的光滑性条件满足时,
可以得到较精确的积分近似法,
7.3.1 Richardson外推法外推法是一种精确度较低的近似公式组合成精确度较高的近似公式的方法,
设 h≠0是任意数,F(h) 是关于步长 h逼近
F*的近似公式,它们的误差估计式为
(1)
* 2 31 2 3( ),..F F h k h k h k h
″
这里,k1,k2,k3… 是一组常数,
按 (1)式,称 F(h)逼近 F*的误差为 O(h),把 h
的幂次称为误差的阶,例如,称为二阶误差,
等等,
我们希望找到一种简便的方法,用近似公式的组合,得到误差阶较高的近似公式,使
(2)
此时,逼近 F*的误差为类似地,用 组合产生逼近 F*的误差的近似公式等,下面我们给出一种具体的组合方法,
*F
2()Oh
()Fh ~ ()Fh
~* ' 2 ' 3
23( ),,,F F h k h k h
~ ()Fh
2()Oh
~ ()Fh
3()Oh
把 (1)式改写为
(3)
用 h/2代替 (3)式中的 h,得
(4)
用 2乘 (4)式再减去 (3)式,消去含 h的项,得
(5)
令,且记
* 2 31 2 3( ),..F F h k h k h k h
23
*
1 2 3( ),.,2 2 4 8
h h h hF F k k k
23
* 2 3
23[ ( ) ( ( ) ( ) ) ] ( ) ( ),,,2 2 2 2
h h h hF F F F h k h k h
1 ( ) ( )F h F h?
2 1 1 1( ) ( ) [ ( ) ( ) ]22
hhF h F F F h
那么 (5)式可写为
(6)
这里,逼近 F*的误差为再用 h/2 代替 h,使 (6)式变为
(7)
用 4乘 (7)式减去 (6)式,消去含 的项,得
(8)
同样记
* 2 3
2 2 3
13( ),,,
24
F F h k h k h
2 ()Fh 2()Oh
* 2 3
2 2 3
13( ),,,
8 3 2
F F h k h k h
2h
*3 22
23
( / 2 ) ( ) 1[ ( ) ],,,
2 3 8
F h F hhF F k h
22
32
( / 2 ) ( )( ) ( )
23
F h F hhF h F
(8)式可以写为
(9)
这里 逼近 F*的误差为还是用 h/2代替 h代入 (9)式后,类似上述过程,可以得到误差为 的一般地,对 k=2,3,…,n,有逼近 F* 的误差为的递推公式
(10)
也称为关于步长 h的外推公式,
表 7-1列出了 k=2,3,4时,按 (10)式产生的计算次序,表中各列左边黑体数字表示序号,
*3
33
1( ),.,
8
F F h k h
3 ()Fh
3()Oh
4()Oh
4 ()Fh
()kOh
11
1 1
( / 2 ) ( )( ) ( )
2 2 1
kk
kk k
F h F hhF h F
()kFh
表 7-1
例 1 设 带余项的差分公式为
()Oh 2()Oh 3()Oh 4()Oh
11,( ) ( )F h F h?
122,( ) ( ) 3,( )22
hhF F F h?
1 2 34,( ) ( ) 5,( ) 6,( )4 4 2
h h hF F F F h?
1 2 3 47,( ) ( ) 8,( ) 9,( ) 1 0,( )8 8 4 2
h h h hF F F F F h?
'
0()fx
(11)
导出具有误差为 的外推公式,
解 令用 h/2代替 h,得
(12)
为消去含 的项,用 4乘 (12)式减去 (11)式,得
2
' '''
0 0 0 0
1( ) [ ( ) ( ) ] ( )
26
hf x f x h f x h f x
h
4
( 5 )
0( ),,,120
h fx
2()jOh
1 0 0
1( ) ( ) [ ( ) ( ) ]
2F h F h f x h f x hh
24
' ''' ( 5 )
0 1 0 0( ) ( ) ( ) ( ),..2 2 4 1 9 2 0
h h hf x F f x f x
2h
2h
从而有
(13)
这里这时,逼近 的误差为,
重复用 h/2代替 h并消去含 的项
,得到逼近 的误差为 的外推公式为
( 4 )
' ( 5 )
0 1 1 03 ( ) 4 ( ) ( ) ( ),..2 16 0
hhf x F F h f x
( 4 )
' ( 5 )
0 2 0( ) ( ) ( ),..480
hf x F h f x
11
21
( / 2 ) ( )( ) ( )
23
F h F hhF h F
2 ()Fh ' 0()fx ( 4 )()Oh
2ih ( 2,3,.,,,1 )ij
46(,,..,)hh如 ' 0()fx 2()jOh
注意 (14)式中第二项的分母为 而不是 (10)式中的,这是由于 (11)式中的余项为关于 的幂次而不是关于 h的幂次,
7.3.2 Romberg求积方法
Romberg求积方法是以复化梯形公式为基础,应用 Richardson外推法导出的数值求积方法,
回忆 7.2.1节的复化梯形公式,分别把积分区
11
1 1
( / 2 ) ( )( ) ( ) 2,3,...,
2 4 1
jj
jj j
F h F hhF h F j k
141j
121j
2h
间 [a,b]分为 1,2,4等分的结果列入表 7-2.
表 7-2
k
1 1
2 2
3 4
12 kkm
k
k
bah
m
1h b a
2
1 ()
2h b a
3
1 ()
4h b a
1 [ ( ) ( ) ]
2
h f a f b?
2
2[ ( ) ( ) 2 ( ) ]2
h f a f b f a h
3
3{ ( ) ( ) 2 [ ( )2
h f a f b f a h
33( 2 ) ( 3 ) ] }f a h f a h
我们还可以进一步推导出它们的递推关系,由可以化为类似地,有一般地,把区间 [a,b]逐次分半 k-1次,
区间长度 (步长 )为,其中,为
2
22[ ( ) ( ) 2 ( ) ]2
hT f a f b f a h
1
2 1 2
1 { [ ( ) ( ) ] ( ) ] }
22
hT f a f b h f a h
1 1 2
1 [ ( ) ]
2 T h f a h
3 2 2 3 3
1 [ ( ( ) ( 3 ) ) ]
2T T h f a h f a h
( 1,2,.,,,)kn?
k
k
bah
m
12 kkm
叙述方便起见,记,那么,
(15)
或
(16)
从而有
(17)
其中,
按外推法的思想,可以把 (15)看成是关于
(1)kkTT?
1
(1)
1
[ ( ) ( ) 2 ( ( )) ]
2
km
k
kk
j
hT f a f b f a j h?
( 1 ) 2 ''( ) ( )
12
b
k k ka
baf x d x T h f
[,]k ab
kh
/2
( 1 ) ( 1 )
11
1
1 [ ( ( 2 1 ) )]
2
km
k k k k
j
T T h f a j h
误差为 的一个近似公式,因此,复化梯形公式的误差公式为
(18)
为消去 项,再取 代替 (18)式中的,
得
(19)
2()kOh
( 1 ) 2 4
12( ),.,
b
k k ka f x d x T K h K h
2 2 2
1
12
ii
i k k i k
ii
K h K h K h
2
kh 1
1
2kkhh kh
( 1 ) 2 2
11 2
1
1()
2
b ii
k i k i kia
i
f x dx T K h K h
22
1
2
11
44
i
k i ki
i
K h K h
用 4乘 (19)式再减去 (18)式,得
(20)
记
(21)
这是误差为 的外推公式,
重复上述过程,将区间逐次分半 k-1次后,可
( 1 ) ( 1 )
( 1 ) 1
1( ) ( )41
b kk
ka
TTf x d x T?
2
2
2
1 ()
3 4 1
i
ik
iki
i
hKh?
1
2
1
2
1 2 4()
3 4 1
i
i
iki
i
Kh
( 1 ) ( 1 )
( 2 ) ( 1 ) 1,2,3,.,,,.
3
kk
kk
TTT T k n
4()Oh
以得到误差为 的外推公式
(22)
当 j=2时
(23)
当 K=2时,有这是 n=2的复化 Simpson公式的,不难验证,对一般的 k,,这里,是 的复化
2()jOh
2,3,.,,,,2,3,.,,,k n j n
( 1 ) ( 1 )
( 2 ) ( 1 ) ( 1 ) ( 1 )1
1
1 ( 4 )
33
kk
k k k k
TTT T T T?
( 2 ) 2
2 2 2
11[ 4 ( ( ) ( ) 2 ( ) ) 2 ( ( ) ( ) ) ]
3 2 2
hT f a f b f a h h f a f b
22
1 [ ( ) 4 ( ) ( ) ]
3 h f a f a h f b
2S
( 2 ) 12kkTS 1
2kS?
12 kn
( 1 ) ( 1 )
( ) ( 1 ) 1
141
jj
jj kk
kk j
TTTT
Simpson公式,
类似地,当 j=3时,
(24)
在实际计算中,经常直接应用 (23)式和形式与 (24)
式相类似的公式进行计算,
所谓 Romberg求积方法,就是由上述两部分组成,第一部分,对积分区间逐次分半 k-1次,用复化梯形求积公式 (16)计算,第二部分,用外推公式 (22)计算
( 2 ) ( 2 )
( 3 ) ( 2 ) ( 2 ) ( 2 )1
1
1 ( 1 6 )
1 5 1 5
kk
k k k k
TTT T T T?
( 1 ) ( 1,2,..,)kTk?
() ( 2,3,...,2,3,...,)jkT k j k
用 Romberg求积方法计算 的计算值的过程如下,首先,令 k=1,区间长度,用梯形求积公式计算 (表 7-3中第一行 );区间分半,令
k=2,区间长度,先按 (16)式计算,再按外推公式 (22)式计算 (表 7-3中第二行 );再区间分半,令 k=3,区间长度,先按 (16)
式计算,再按 (22)式计算 (表 7-3中第三行 )等等,逐次分半区间 k次后的计算结果如表
7-3所示 (见下页 ).
()ba f x dx?
1h b a
(1)1T
21
1
2hh? (1)2T
(2)2T
3 2 12
11
22h h h
(1)3T
( 2 ) ( 3 )33,TT
表 7-3
:
:
:
,…….
表 7-3中 的计算按行 (k的序号 )进行,每行第
1个元素 用复化梯形公式 (16)计算,其他元素 均按 (22)式用 与 的组合得到,在实际应用中,往往根据实际问题对计
(1)1T
(1)2T
(2)2T
(1)3T (2)
3T (3)3T
1h
2h
3h
kh (1)kT
(2)kT (3)kT ()kkT
(1)kT
()jkT
( 1)jkT? ( 1 )1 ( 2,3,...,)jkT j k()j
kT
算精确度的要求来确定区间逐次分半的次数,
常用不等式
(25)
作为达到精确度要求的判断准则,这里,?是给定的一个小的正数,
例 2 用 Romberg求积方法计算
(26)
的近似值,给定?=0.001
解 首先令区间长度 h1=1,用梯形求积公式计算
1
0
1 l n 2 0,6 9 3 1 4 7 1,,,
1 dxx
( 1 ) 1
1 [ ( 0 ) ( 1 ) ] 0,7 5 0 0 0 0 02
hT f f
( ) ( 1 )[]jj
kkTT?
区间 [0,1]分半,令区间长度,按 16式计算再按 (23)式计算这时未达到精确度要求,
为此,再将区间分半,令区间长度按 (16)式计算
21
11
22hh
( 1 ) ( 1 )
2 1 2
1 ( 0,5 ) 0,7 0 8 3 3 3 3,
2T T h f
( 2 )
2
1 ( 4 0,7 0 8 3 3 3 3 0,7 5 0 0 0 0 0 ) 0,6 9 4 4 4 4 4
3T
( 1 ) ( 2 )22[ ] 0,01 38 88 9TT
32
11
24hh
按 j=2和 j=3的外推公式 (23)和 (24),分别用和 的组合得到 以及用 和 的组合得到,即以及这时,
已满足不等式 (25)的要求,作为积分 (26)式的近似,其误差为,
( 1 ) ( 1 )
3 2 3
1 [ (0,2 5 ) (0,7 5 ) ] 0,6 9 7 0 2 3 7
2T T h f f
(1)2T
(1)3T (2)
3T
(2)3T(2)2T
(3)3T
( 2 ) ( 1 ) ( 1 )
3 3 2
1 ( 4 ) 0,6 9 3 2 5 3 8
3T T T
( 3 ) ( 2 ) ( 2 )
3 3 2
1 ( 1 6 )
15T T T
( 1 ) ( 2 )
33[ ] 0,00 06 36 4TT
63()Oh
(3)3T
下面给出用 Romberg求积方法计算近似值的计算步骤,用二维数组 T的元素 存放表 7-3中的,
1 输入,积分区间端点
2 令,计算
3 令
4 令,计算
5 for j=2,3,…,k
5.1 计算
end for (j)
()ba f x dx?
jkT
()jkT
,;ab?
1h b a 1
11 [ ( ) ( ) ]2
hT f a f b
1
12,
2kk h h
12 kkm
/2
( 1 ) ( 1 )
11
1
1 [ ( ( 2 1 ) )]
2
km
k k k k
j
T T h f a j h
( 1 ) ( 1 )
( ) ( 1 ) 1
141
jj
jj kk
kk j
TTTT
6 if,then
goto 8
end if
7 令 k=k+1,,goto 4
8 输出
9 end
( ) ( 1 )[]jjkkTT
1
1
2kkhh
jkT
7.5 Gauss求积公式
7.5.1 引言求积公式
(1)
当求积系数,求积节点 都可以自由选取时,其代数精确度最高可以达到多少次?
下面的引理可以回答上述问题,
1
( ) ( ) ( ) ( ) (,)
nb
kka
k
I f x f x dx A f x R f
1{}nkkA? 1{}nkkx?
引理 1 当求积系数 和求积节点都可以自由选取时,n点的求积公式 (1)的代数精确度最高可以达到 2n-1次,
证 假设求积公式 (1)具有 m次代数精确度,
即对任意的 m次代数多项式求积公式 (1)的精确成立,于是成立等式即若记 (2)
0
()
m
i
m i m
i
p x a x P
1
( ) ( ) ( )
mb
m k m ka
k
x p x dx A p x?
1
1 1 0
01
( ) (,.,)
mn b
i m m
i k m k m k ka
ik
a x x dx A a x a x a x a
( ),0,1,2,..,,b i i
a
x x d x i m
1{}nkkA? 1{}nkkx?
则 (2)式成为
(3)
由于系数 的任意性,故使 (3)式成为恒等式的充要条件是
(4)
(4)式的待定系数有 2n个,所以确定待定系数的
1 1 1 1 0 0...m m m ma a a a
10
1 1 1
...
n n n
m
m k k k k k
k k k
a A x a A x a A
1 1 0,,.,,,,mma a a a?
1 2 0
1 1 2 2 1
1 1 2 2
n
nn
m m m
n n m
A A A
A x A x A x
A x A x A x
独立条件至多给出 2n个,从而可知 m至多为
2n-1.
定义 1 n点的求积公式 (1)具有 2n-1次代数精确度 (或称为具有最高的代数精确度 )时,称为 Gauss型求积公式,
Gauss型求积公式的求积节点,称为
Gauss点,它们可以通过求区间 [a,b]上带权?(x)
的 n次正交多项式 的 n个根获得,所以先介绍正交多项式及其性质,然后讨论 Gauss型求积公式的构造,等等,
1{}nkkx?
()ngx
7.5.2 正交多项式及其性质定义 2 若
(1),则称函数 f(x)和 g(x)在区间
[a,b]上正交,
(2),则称函数 f(x)和 g(x)在区间 [a,b]上带权?(x)正交,
(3)代数多项式序列 (下标 k为多项式的次数,表示 k次多项式 ),在区间 [a,b]上满足当 m?n
当 m=n
则称多项式序列 为区间 [a,b]上带权
( ) ( ) 0ba f x g x d x
( ) ( ) ( ) 0ba x f x g x d x
0{ ( ) }kkgx
()kgx
2
0,
( ) ( ) ( )
( ) ( ( ) ) 0
b
bnm
a
na
x g x g x d x
x g x d x
0{ ( ) }kkgx
(x)的正交多项式序列,
定义 3 若 n次记多项式 中含 项的系数为,则称 为 的 首次系数 ; 时,称为 首次系数为 1的 n次多项式,
正交多项式有如下性质,
性质 1 若 是区间 [a,b]上带权?(x)的正交多项式序列,则它们线行无关,
证 对任意的 x?[a,b],若,在式子两边同乘?(x)g1(x)(L=0,1,..n),并从 a到 b积分,由的正交性定义 1中的 (3)可知必有
()ngx nx
nd nd
()ngx 0
nd?
* ()() n
n
n
gxgx
d?
0{ ( ) } nkkgx?
0
( ) 0
n
kk
k
c g x
0{ ( )} nkkgx? 0
lc?
L=0,1,..,n.故正交多项式序列 线性无关,
由性质 1可知,若 为 [a,b]上带权?(x)的正交多项式序列,则序列 可以作为空间的一组基函数,即 中的任一元素可由它们线性表出,
其中 为组合系数,
性质 2 若 为 [a,b]上带权的正交?(x)多项式序列,且,则
0{ ( )} nkkgx?
0{ ( )} nkkgx?
[,]nP a b [,]nP a b ()
npx
0
( ) ( )
n
n k k
k
p x a g x
0{}nkka?
0{ ( )} nkkgx?
0{ ( ) }kkgx
( ) [,]nq x P a b?
(1)
(2)
事实上,由性质 1,.由 的正交性定义容易证得 (1).证 (2)也是类似的,
为方便起见,记下面,不加证明地给出正交多项式如下的性质,
性质 3 [a,b]上带权函数?(x)的正交多项式序列 相邻三项的递推关系为其中,,
0
( ) ( )
n
kk
k
q x a g x
{ ( )}ngx
(,) ( ) ( ) ( )baf g x f x g x d x
0{ ( )}kkgx
1 1 1( ) ( ( ) ) ( ) ( ),1,2,..,,n n n n n ng x a x g x g x n
1n
n
n
d
a
d
1 (,)
(,)
n n n
n
n n n
d g x g
d g g?
( ) ( ) ( ) 0,1,2,..,b ka x q x g x d x k n n
( ) ( ) 0,0,1,2,...,1b ina x g x x dx i n
为 的首项系数,
即为性质 4 [a,b]上带权函数 的正交多项式序列 中任意相邻两个正交多项式 和的根相间,
若记,的根分别为,
则所谓 与 的根相同,即是指这两个正交多项式的根有如下的关系,
11
1 2
11
(,),
(,)
n n n n
n
n n n
d d g g
d g g?
11,,n n nd d d 11( ),( ),( )n n ng x g x g x
( ) (,1 )kg x k n n
0{ ( )}kkgx ()ngx
1 ()ngx?
()ngx 1 ()ngx? () 1{}nniix? ( 1) 11{}nnjjx
1 ()ngx?()ngx
()x?
性质 5 (1) 区间 [a,b]上带权函数?(x)的正交多项式序列 与 对应元素之间只相差一个比例常数,
(2)区间 [a,b]上带权函数?(x)首项系数为 1的正交多项式序列 唯一,
常见的正交多项式有 Legendre(勒让德 )多项式,Hermite多项式,Chebyshev多项式以及 Jacobi多项式,Chebyshev多项式在 5.2节已
( 1 ) ( ) ( 1 ) ( ) ( 1 )
1 1 2,1,2,...,1
n n n n n
i i i i ix x x x x i n
0{ ( ) }nnfx 0{ ( ) }nngx
* 0{ ( )}nngx
详细讨论,这里主要介绍 Legendre多项式一,Legendre多项式隐式表达式显式表达式
()nPx
0
2
( ) 1,
1 ( 1 )
( ),1,2,.,,
2!
nn
n nn
Px
dx
P x n
n d x
0
2
0
( ) 1,
1 ( 2 2 )
( ) ( 1 ),
2 ! ( ) ! ( 2 ) !
N
j n j
n n
j
Px
nj
P x x
j n j n j
1,2,...n?
当 n为偶数时,
当 n为奇数时,
Legendre多项式的主要性质有
(1)n次 Legendre多项式 的首项系数当 x=1,
当 x=-1.
(3)正交性为,为区间 [-1,1]上带权函数?(x)?1的正交多项式序列,且有
/ 2,
( 1 ) / 2,
n
N
n
其 中
()nPx
2
( 2 ) !( ) ;
2 ( ! )n n
nx
n?d
1,
( 2 ) ( 1 )
( 1 ),n n
P
0{ ( ) }nnPx
当 m?n
当 m=n
(4) Legendre多项式相邻三项的递推关系为二,Legendre多项式将隐式表达式
1
1
0,
(,) ( ) ( ) 2
,
21
m n m nP P P x P x d x
n
0
1
11
( ) 1,
( ),
21
( ) ( ) ( ),1,2,.,,
11
n n n
Px
P x x
nn
P x x P x P x n
nn
()Ln x
()() n n xx
n n
d x eL x e
dx
将隐式表达式中 n阶导数用乘积导数的
Leibniz公式可得显式表达式
Legendre多项式的主要性质有
(1)n次 Legendre多项式 的首项系数
(2)正交性为,为区间 [0,+?)上带权函数 的正交多项式序列,且有
0
!( ) ( 1 )
( ) !
n
n k k n k
nn
k
nL x C x
nk
()nLx ( 1 ) nnd
0{ ( ) }nnLx
() xxe
20
0,
(,) ( ) ( )
( ! ),
x
n m n mL L e L x L x d x n
当 mn
当 m=n
权函数 的正交多项式序列,且有
(3) Hermite多项式相邻三项的递推关系为四,Jacobi多项式
Jacobi多项式是在区间 [-1,1]上带权函数的正交多项式,其中?>-1,?>-1
2() xxe
2 0,,
(,) ( ) ( )
2 !,
x
m n m n n
mn
H H e H x H x d x
n?
当当 m=n;
0
1
11
( ) 1,
( ) 2,
( ) 2 ( ) 2 ( ),1,2,...n n n
Hx
H x x
H x x H x nH x n
( ) ( 1 ) ( 1 )x x x
(3) Legendre多项式相邻三项的递推关系为三,Hermite多项式表达式
Hermite多项式的主要性质有
(1)n次 Hermite多项式 的首项系数
(2)正交性为,为区间 上带
()nHx
()nHx
0
1
2
11
( ) 1,
( ) 1,
( ) ( 1 2 ) ( ) ( ),1,2,...n n n
Lx
L x x
L x n x L x n L x n
2
2 ()( ) ( 1 )
nx
nx
n n
deH x e
dx
( 2 ) nnd
0{ ( )}nnHx (,)
有的书籍文献把 Jacobi多项式记为即 n次 Jacobi多项式表示为其中 或,两种系数推出两种 Jacobi多项式,详细的情形请参阅文献 [27].
(,) ()nJx
(,) ( ) ( 1 ) ( 1 ) [ ( 1 ) ( 1 ) ]
n
nn
nn n
dJ x K x x x x
dx
1 2 1
( 2 ) ! 2n
nK
n
( 1 )
2!
n
n nK n
7.5.3 Gauss型求积公式由 7.5.1中的引理 1和定义 1可知 n点的求积公式 (1)若具有最高的代数精确度,或具有 2n-1
次的代数精确度成为 Gauss型求积公式,到底求积公式 (1)的求积节点 和求积系数如何选取,才能使之成为 Gauss型求积公式?
定理 7.5.1 求积公式 (1)中的 n个求积节点
,取在区间 [a,b]上带权函数?(x)的 n次正交多项式 的 n个根成为 Gauss型求积公式,
证 设,[a,b]上带权函数?(x)的
1{}nkkx? 1{}nkkA?
1{}nkkx?
()ngx
21( ) [,]nf x P a b
n次正交多项式 的 n个根记为,记的首项系数为,由定义 2有因此,
(5)
其中,在 (5)式两边同乘?(x),并从
a到 b积分,由正交多项式的性质可知,含项的积分为零,所以
(6)
注意到当 作为插值节点时建立的 n点插值
()ngx 1{}nkkx?
1{}nkkx?
nd
*
1
()( ) ( )n n
nk
k n
gxg x x x
d?
*( ) ( ) ( ) ( ),nf x q x g x r x
1( ),( ) [,]nq x r x P a b
* ( ) ( )ng x q x
( ) ( ) ( ) ( )bbaax f x d x x r x d x
求积公式至少具有 n-1次代数精确度,而,
所以
(7)
又由 (5)式可知,
即 (8)
综合 (6),(7),(8)式可知,当 时,求积公式 (1)
成立,
1
( ) ( )
n
n k k
k
I f A f x
1( ) [,]nr x P a b
1
( ) ( ) ( )
n
k k n
k
I r A r x I r
( ) ( ) ( 1,2,..,,)kkf x r x k n
( ) ( )nnI f I r?
21( ) [,]nf x P a b
1
( ) ( ) ( )
nb
kka
k
x f x d x A r x?
用 n点 Gauss求积公式
(9)
之值近似积分值 有下面的误差估计,
定理 7.5.2 若,则 Gauss型求积公式 (1)的误差估计 R(?,f)为其中证明略,
在稍后讨论 Gauss积分值数列的收敛性
1
( ) ( )
n
n k k
k
I f A f x
()nIf
2( ) [,]nf x C a b?
2
*2()(,) ( ) [ ( ) ]
( 2 ) !
n b
na
fR f x g x dx
n
*
1
( ) ( )
n
nk
k
g x x x
等问题时,需要用到 Gauss型求积公式的求积系数 大于零的结论,这里用下面的定理给出,
定理 7.5.3 Gauss型求积公式的求积系数大于零,
证 令,这里 为区间 [a,b]
上带权函数?(x)的 n次正交多项式 的 n个根,显然
1{}nkkA?
1{}nkkx?
()ngx
1{}nkkA?
2
1()
n
i
i
k
xx
fx
xx
1,
0,
() ()
n
j
ki
i
ik
fx xx
当 jk
,当 j=k
由于,所以对 f(x)求积公式 (1)精确成立,即因为所以在 7.2节,我们讨论了复化梯形求积公式和复化 Simpson求积公式的收敛性,那么 Gauss
1
( ) ( ) ( ) ( ) ( )
nb
n j ja
j
I f x f x dx I f A f x?
2
1,
()
n
k k i
i
ik
A x x
,
2
1,
( ) 0,( ) ( ) 0
n b
ki a
i
ik
x x x f x d x?
0,1,2,3,.,,,kA k n
22() nf x P
型求积公式被积函数 f(x)应当满足什么条件才收敛呢? Gauss型求积公式的收敛性问题由下面的定理给出,
定理 7.5.4 若 f(x)?[a,b],则 Gauss型求积公式所求积分值序列 收敛于积分值 I(f),即证 因为 f(x)?[a,b],由 Weierstrass定理对任意的?1>0,存在,使得
1
{ ( ) ( ) }
n
n k k
k
I f A f x
1
l im ( ) ( ) ( )
n b
kk an
k
A f x x f x d x?
( ) [,]mmp x P a b?
(10)
对任意的 x?[a,b]成立,
由于公式 (1)为 Gauss型求积公式时具有
2n-1次代数精确度,取 N>(m+1)/2,故当 n>N时,
即 m<2N-1<2n-1时,有
(11)
成立,于是由 (11)式可知,而由 (10)式,有
1( ) ( )mf x p x
( ) ( )n m mI p I p?
( ) ( ) ( ) ( ) ( ) ( )n m m n mI f I f I f I p I p I p
( ) ( )n m nI p I f
( ) ( ) 0m n mI p I p
和,从而因为,记
( ) ( ) ( ) [ ( ) ( ) ]bmm
a
I f I p x f x p x d x
1
1
n
k
k
A?
1( ) ( )k m kf x p x
1
( ) ( ) [ ( ) ( ) ]
n
n m n k m k k
k
I p I f A p x f x
1
1
n
k
k
A?
0,1,2,..,,kA k n
故即
1 1 1
11
( ) ( ) 0 2
nn
n k k
kk
I f I f A A C
l i m ( ) ( ),nn I f I f
1
()
nb
ka
k
C x dx A?
7.5.4 Gauss型求积公式的构造与应用定理 7.5.1实际上给出了构造 Gauss型求积公式的一种方法,这种方法,当给定了积分区间
[a,b]和权函数?(x)以后构造 n个点的 Gauss型求积公式,先求出区间 [a,b]上带权函数?(x)的 n次正交多项式,然后用多项式求根的方法求出的 n个根,从而获得了求积节点为了求得求积系数,将 n个求积节点代入方程组 (4)中的前 n个方程并加以求解,
即解线性代数方程组
()ngx
()ngx
1{}nkkx? 1{}nkkx?
1{}nkkA?
1{}nkkx?
求得求积系数,完成 Gauss型求积公式的构造,
表 7-5为 Gauss-Legendre求积公式的求积系数 和求积节点 的一个表,而表 7-6和表
7-7则分别是 Gauss-Legendre和 Gauss-Hermit
-te求积公式的求积系数 和求积节点 的一
1{}nkkA?
1 2 0
1 1 2 2 1
1 1 2 2
n
nn
m m m
n n m
A A A
A x A x A x
A x A x A x
kA
kx
kxkA
个表,有些数值分析的书籍还给出了 Gauss-
Chebyshev求积公式的求积节点与求积系数,
表见 P366.
需要指出的是,不同求积公式的求积系数与求积节点,积分区间和权函数是不同的,Gauss-
Lagendre求积公式的积分区间为 [-1,1],权函数
(x)=1.而 Gauss-Legendre求积公式的积分区间为 [0,+?),权函数,Gauss-Hermitte
求积公式的积分区间为 (-?,+?),权函数
Gauss-Chebyshev求积公式的积分区间为 [-1,1]
() xxe
2() xxe
权函数这里需要指出的另一点是 Gauss-
Chebyshev求积公式的求积系数是相同的,例如,n点的 Gauss-Chebyshev求积公式,它的 n个求积系数 都是,即,而 n个求积节点则为正是因为 Gauss-Chebyshev求积公式的求积系数相同,所以在实际计算时,乘法的次数只需一次,节省了 n-1次的乘法运算,
2
1()
1x x
kA
n
12 nA A A n
21c o s,1,2,.,,,
2k
kx k n
n?
例 1 求 使求积公式具有三次代数精确度,
问题是构造区间 [0,1]上带权函数的两点 Gauss型求积公式,
解 方法 1 容易计算出当时 的积分值分别为,所求公式具有 3次代数精确度,故可得为未知数的方程组为
1 2 1 2,,,,A A x x
1
1 1 2 20 ( ) ( ) ( ) (,)x f x d x A f x A f x R x f
()xx
23( ) 1,,,f x x x x?
1
0 ()x f x dx?
2222,,,,
3 5 7 9
1 2 1 2,,,,A A x x
(1)
(2)
(3)
(4)
又因为 为 Gauss型求积公式的求积节点,
所以它们是区间 [0,1]上带权函数 且首项系数为 1的二次正交多项式 的两个根,
不妨记,为此又因为 必须满足方程 (2)(3)(4),所以
12
1 1 2 2
22
1 1 2 2
33
1 1 2 2
2 / 3
2 / 5
2 / 7
2 / 9
AA
A x A x
A x A x
A x A x
12,xx
()xx
*2 ()gx
*22 ()g x x px q
**
2 1 2 2( ) ( ) 0g x g x
12,xx
由可得关于 p,q的方程组为解此方程组得
**
1 2 1 2 2 2
**
1 1 2 1 2 2 2 2
( 1 ) ( 2 ) ( 3 ) ( ) ( ) 0
( 2 ) ( 3 ) ( 4 ) ( ) ( ) 0
q p A g x A g x
q p A x g x A x g x
2 2 2
3 5 7
2 2 2
5 7 9
qp
qp
10
9
5
21
p
q
将所求 p,q代入,求得其根为再将所求 代入方程 (1)(2),联立解得为此,公式为所求具有 3次代数精确度的求积公式,
*22 ()g x x px q
1
2
0,2 8 9 9 4 9 1 9 8
0,8 2 1 1 6 1 9 1 2
x
x
12,xx
1
0 ( ) 0,2 7 7 5 5 5 9 9 7 ( 0,2 8 9 9 4 9 1 9 8 )x f x d x f
1
2
0,2 7 7 5 5 5 9 9 7
0,3 8 9 1 1 0 6 6 9
A
A
0.389110669 ( 0.821161912) (,)f R x f
方法 2 求出区间 [0,1]上带权函数的二次正交多项式,并求出其根,获得求积节点,再求方程组 (12)前两个方程组成的方程组获得求积系数,
由于 和 都是的线性无关组,所以考虑由 求得所需的正交多项式,这种方法可以称作将它们正交化,
记
()xx
*2()gx
12,xx
12,AA
***0 1 2{ ( ),( ),( ) }g x g x g x
2{1,,}xx
2 [0,1]P
1 1 / 2
20(,) ( ) ( ),(,)f g x f x g x d x f f f
令所以又令 这样即如此作出的 与 正交,也与 正交,
由于所以这样,
* * * *0 0 0 0
2( ) 1,(,) 2 / 3,2 / 3g x g g g
*1 0 0( ) (,),g x x x e e
*
0
0 *
02
() 3,
2()
gxe
gx
*1 0 0 0 0 0 0( ( ),) ( (,),) (,) (,) 0g x e x x e e e x e x e
*1 ()gx
0e
0 ()gx
1
0 0
3 2 3(,)
2 5 2x e x x d x
**
11 2
38( ),( )
5 175
g x x g x
同理可令并且容易验证,与 正交,容易求得所以
*
1
1 *
1 2
() 17 5 3 17 5
8 5 8()
gx
ex
gx
* 2 2 22 1 1 0 0( ) (,) (,),g x x x e e x e e
*2 ()gx 10,ee
2
1
1 7 5 1 6(,)
8 3 1 5xe 2 0 32(,),
27
xe
*2
2
1 0 5()
9 2 1g x x x
第二种方法与第一种方法求出的 是一样的,后面的求解过程相同,这里略去,
方法 2是一种将一组线性无关函数组正交化而得到正交多项式 的方法,高于 2次的正交多项式用这样的方法同样可以得到,这样求正交多项式在实际应用时是方便的,
这里需要附带说明的是,Gauss-Legendre
,Gauss-Chebyshev求积公式作数值求积精度不够时,可以采取将积分区间 [a,b]若干等分后,
将每一个子区间映射到区间 [-1,1]上再用相同
*2 ()gx
*2()gx
2{1,,}xx
节点数的求积公式进行数值求积的计算,通常会得到精度较好的计算结果,
Gauss型求积公式具有数值结果精度高,收敛得以保证、计算简便、易于在计算机上实现等优点,并且在积分区间 [a,b]有限时便于推广到高维数值积分,不足之处是公式的构造比较困难,另一个是由于相邻次数的正交多项式的根,从而造成增加求积节点以提高计算结果的精度时,原先所有求积节点上的函数值全部无用,所以在具体应用 Gauss求积公式计算数值积分时,n取得都较小,
Romberg求积方法是在积分区间逐次分半的过程中利用外推法产生的一种数值积分方法,当被积函数的光滑性条件满足时,
可以得到较精确的积分近似法,
7.3.1 Richardson外推法外推法是一种精确度较低的近似公式组合成精确度较高的近似公式的方法,
设 h≠0是任意数,F(h) 是关于步长 h逼近
F*的近似公式,它们的误差估计式为
(1)
* 2 31 2 3( ),..F F h k h k h k h
″
这里,k1,k2,k3… 是一组常数,
按 (1)式,称 F(h)逼近 F*的误差为 O(h),把 h
的幂次称为误差的阶,例如,称为二阶误差,
等等,
我们希望找到一种简便的方法,用近似公式的组合,得到误差阶较高的近似公式,使
(2)
此时,逼近 F*的误差为类似地,用 组合产生逼近 F*的误差的近似公式等,下面我们给出一种具体的组合方法,
*F
2()Oh
()Fh ~ ()Fh
~* ' 2 ' 3
23( ),,,F F h k h k h
~ ()Fh
2()Oh
~ ()Fh
3()Oh
把 (1)式改写为
(3)
用 h/2代替 (3)式中的 h,得
(4)
用 2乘 (4)式再减去 (3)式,消去含 h的项,得
(5)
令,且记
* 2 31 2 3( ),..F F h k h k h k h
23
*
1 2 3( ),.,2 2 4 8
h h h hF F k k k
23
* 2 3
23[ ( ) ( ( ) ( ) ) ] ( ) ( ),,,2 2 2 2
h h h hF F F F h k h k h
1 ( ) ( )F h F h?
2 1 1 1( ) ( ) [ ( ) ( ) ]22
hhF h F F F h
那么 (5)式可写为
(6)
这里,逼近 F*的误差为再用 h/2 代替 h,使 (6)式变为
(7)
用 4乘 (7)式减去 (6)式,消去含 的项,得
(8)
同样记
* 2 3
2 2 3
13( ),,,
24
F F h k h k h
2 ()Fh 2()Oh
* 2 3
2 2 3
13( ),,,
8 3 2
F F h k h k h
2h
*3 22
23
( / 2 ) ( ) 1[ ( ) ],,,
2 3 8
F h F hhF F k h
22
32
( / 2 ) ( )( ) ( )
23
F h F hhF h F
(8)式可以写为
(9)
这里 逼近 F*的误差为还是用 h/2代替 h代入 (9)式后,类似上述过程,可以得到误差为 的一般地,对 k=2,3,…,n,有逼近 F* 的误差为的递推公式
(10)
也称为关于步长 h的外推公式,
表 7-1列出了 k=2,3,4时,按 (10)式产生的计算次序,表中各列左边黑体数字表示序号,
*3
33
1( ),.,
8
F F h k h
3 ()Fh
3()Oh
4()Oh
4 ()Fh
()kOh
11
1 1
( / 2 ) ( )( ) ( )
2 2 1
kk
kk k
F h F hhF h F
()kFh
表 7-1
例 1 设 带余项的差分公式为
()Oh 2()Oh 3()Oh 4()Oh
11,( ) ( )F h F h?
122,( ) ( ) 3,( )22
hhF F F h?
1 2 34,( ) ( ) 5,( ) 6,( )4 4 2
h h hF F F F h?
1 2 3 47,( ) ( ) 8,( ) 9,( ) 1 0,( )8 8 4 2
h h h hF F F F F h?
'
0()fx
(11)
导出具有误差为 的外推公式,
解 令用 h/2代替 h,得
(12)
为消去含 的项,用 4乘 (12)式减去 (11)式,得
2
' '''
0 0 0 0
1( ) [ ( ) ( ) ] ( )
26
hf x f x h f x h f x
h
4
( 5 )
0( ),,,120
h fx
2()jOh
1 0 0
1( ) ( ) [ ( ) ( ) ]
2F h F h f x h f x hh
24
' ''' ( 5 )
0 1 0 0( ) ( ) ( ) ( ),..2 2 4 1 9 2 0
h h hf x F f x f x
2h
2h
从而有
(13)
这里这时,逼近 的误差为,
重复用 h/2代替 h并消去含 的项
,得到逼近 的误差为 的外推公式为
( 4 )
' ( 5 )
0 1 1 03 ( ) 4 ( ) ( ) ( ),..2 16 0
hhf x F F h f x
( 4 )
' ( 5 )
0 2 0( ) ( ) ( ),..480
hf x F h f x
11
21
( / 2 ) ( )( ) ( )
23
F h F hhF h F
2 ()Fh ' 0()fx ( 4 )()Oh
2ih ( 2,3,.,,,1 )ij
46(,,..,)hh如 ' 0()fx 2()jOh
注意 (14)式中第二项的分母为 而不是 (10)式中的,这是由于 (11)式中的余项为关于 的幂次而不是关于 h的幂次,
7.3.2 Romberg求积方法
Romberg求积方法是以复化梯形公式为基础,应用 Richardson外推法导出的数值求积方法,
回忆 7.2.1节的复化梯形公式,分别把积分区
11
1 1
( / 2 ) ( )( ) ( ) 2,3,...,
2 4 1
jj
jj j
F h F hhF h F j k
141j
121j
2h
间 [a,b]分为 1,2,4等分的结果列入表 7-2.
表 7-2
k
1 1
2 2
3 4
12 kkm
k
k
bah
m
1h b a
2
1 ()
2h b a
3
1 ()
4h b a
1 [ ( ) ( ) ]
2
h f a f b?
2
2[ ( ) ( ) 2 ( ) ]2
h f a f b f a h
3
3{ ( ) ( ) 2 [ ( )2
h f a f b f a h
33( 2 ) ( 3 ) ] }f a h f a h
我们还可以进一步推导出它们的递推关系,由可以化为类似地,有一般地,把区间 [a,b]逐次分半 k-1次,
区间长度 (步长 )为,其中,为
2
22[ ( ) ( ) 2 ( ) ]2
hT f a f b f a h
1
2 1 2
1 { [ ( ) ( ) ] ( ) ] }
22
hT f a f b h f a h
1 1 2
1 [ ( ) ]
2 T h f a h
3 2 2 3 3
1 [ ( ( ) ( 3 ) ) ]
2T T h f a h f a h
( 1,2,.,,,)kn?
k
k
bah
m
12 kkm
叙述方便起见,记,那么,
(15)
或
(16)
从而有
(17)
其中,
按外推法的思想,可以把 (15)看成是关于
(1)kkTT?
1
(1)
1
[ ( ) ( ) 2 ( ( )) ]
2
km
k
kk
j
hT f a f b f a j h?
( 1 ) 2 ''( ) ( )
12
b
k k ka
baf x d x T h f
[,]k ab
kh
/2
( 1 ) ( 1 )
11
1
1 [ ( ( 2 1 ) )]
2
km
k k k k
j
T T h f a j h
误差为 的一个近似公式,因此,复化梯形公式的误差公式为
(18)
为消去 项,再取 代替 (18)式中的,
得
(19)
2()kOh
( 1 ) 2 4
12( ),.,
b
k k ka f x d x T K h K h
2 2 2
1
12
ii
i k k i k
ii
K h K h K h
2
kh 1
1
2kkhh kh
( 1 ) 2 2
11 2
1
1()
2
b ii
k i k i kia
i
f x dx T K h K h
22
1
2
11
44
i
k i ki
i
K h K h
用 4乘 (19)式再减去 (18)式,得
(20)
记
(21)
这是误差为 的外推公式,
重复上述过程,将区间逐次分半 k-1次后,可
( 1 ) ( 1 )
( 1 ) 1
1( ) ( )41
b kk
ka
TTf x d x T?
2
2
2
1 ()
3 4 1
i
ik
iki
i
hKh?
1
2
1
2
1 2 4()
3 4 1
i
i
iki
i
Kh
( 1 ) ( 1 )
( 2 ) ( 1 ) 1,2,3,.,,,.
3
kk
kk
TTT T k n
4()Oh
以得到误差为 的外推公式
(22)
当 j=2时
(23)
当 K=2时,有这是 n=2的复化 Simpson公式的,不难验证,对一般的 k,,这里,是 的复化
2()jOh
2,3,.,,,,2,3,.,,,k n j n
( 1 ) ( 1 )
( 2 ) ( 1 ) ( 1 ) ( 1 )1
1
1 ( 4 )
33
kk
k k k k
TTT T T T?
( 2 ) 2
2 2 2
11[ 4 ( ( ) ( ) 2 ( ) ) 2 ( ( ) ( ) ) ]
3 2 2
hT f a f b f a h h f a f b
22
1 [ ( ) 4 ( ) ( ) ]
3 h f a f a h f b
2S
( 2 ) 12kkTS 1
2kS?
12 kn
( 1 ) ( 1 )
( ) ( 1 ) 1
141
jj
jj kk
kk j
TTTT
Simpson公式,
类似地,当 j=3时,
(24)
在实际计算中,经常直接应用 (23)式和形式与 (24)
式相类似的公式进行计算,
所谓 Romberg求积方法,就是由上述两部分组成,第一部分,对积分区间逐次分半 k-1次,用复化梯形求积公式 (16)计算,第二部分,用外推公式 (22)计算
( 2 ) ( 2 )
( 3 ) ( 2 ) ( 2 ) ( 2 )1
1
1 ( 1 6 )
1 5 1 5
kk
k k k k
TTT T T T?
( 1 ) ( 1,2,..,)kTk?
() ( 2,3,...,2,3,...,)jkT k j k
用 Romberg求积方法计算 的计算值的过程如下,首先,令 k=1,区间长度,用梯形求积公式计算 (表 7-3中第一行 );区间分半,令
k=2,区间长度,先按 (16)式计算,再按外推公式 (22)式计算 (表 7-3中第二行 );再区间分半,令 k=3,区间长度,先按 (16)
式计算,再按 (22)式计算 (表 7-3中第三行 )等等,逐次分半区间 k次后的计算结果如表
7-3所示 (见下页 ).
()ba f x dx?
1h b a
(1)1T
21
1
2hh? (1)2T
(2)2T
3 2 12
11
22h h h
(1)3T
( 2 ) ( 3 )33,TT
表 7-3
:
:
:
,…….
表 7-3中 的计算按行 (k的序号 )进行,每行第
1个元素 用复化梯形公式 (16)计算,其他元素 均按 (22)式用 与 的组合得到,在实际应用中,往往根据实际问题对计
(1)1T
(1)2T
(2)2T
(1)3T (2)
3T (3)3T
1h
2h
3h
kh (1)kT
(2)kT (3)kT ()kkT
(1)kT
()jkT
( 1)jkT? ( 1 )1 ( 2,3,...,)jkT j k()j
kT
算精确度的要求来确定区间逐次分半的次数,
常用不等式
(25)
作为达到精确度要求的判断准则,这里,?是给定的一个小的正数,
例 2 用 Romberg求积方法计算
(26)
的近似值,给定?=0.001
解 首先令区间长度 h1=1,用梯形求积公式计算
1
0
1 l n 2 0,6 9 3 1 4 7 1,,,
1 dxx
( 1 ) 1
1 [ ( 0 ) ( 1 ) ] 0,7 5 0 0 0 0 02
hT f f
( ) ( 1 )[]jj
kkTT?
区间 [0,1]分半,令区间长度,按 16式计算再按 (23)式计算这时未达到精确度要求,
为此,再将区间分半,令区间长度按 (16)式计算
21
11
22hh
( 1 ) ( 1 )
2 1 2
1 ( 0,5 ) 0,7 0 8 3 3 3 3,
2T T h f
( 2 )
2
1 ( 4 0,7 0 8 3 3 3 3 0,7 5 0 0 0 0 0 ) 0,6 9 4 4 4 4 4
3T
( 1 ) ( 2 )22[ ] 0,01 38 88 9TT
32
11
24hh
按 j=2和 j=3的外推公式 (23)和 (24),分别用和 的组合得到 以及用 和 的组合得到,即以及这时,
已满足不等式 (25)的要求,作为积分 (26)式的近似,其误差为,
( 1 ) ( 1 )
3 2 3
1 [ (0,2 5 ) (0,7 5 ) ] 0,6 9 7 0 2 3 7
2T T h f f
(1)2T
(1)3T (2)
3T
(2)3T(2)2T
(3)3T
( 2 ) ( 1 ) ( 1 )
3 3 2
1 ( 4 ) 0,6 9 3 2 5 3 8
3T T T
( 3 ) ( 2 ) ( 2 )
3 3 2
1 ( 1 6 )
15T T T
( 1 ) ( 2 )
33[ ] 0,00 06 36 4TT
63()Oh
(3)3T
下面给出用 Romberg求积方法计算近似值的计算步骤,用二维数组 T的元素 存放表 7-3中的,
1 输入,积分区间端点
2 令,计算
3 令
4 令,计算
5 for j=2,3,…,k
5.1 计算
end for (j)
()ba f x dx?
jkT
()jkT
,;ab?
1h b a 1
11 [ ( ) ( ) ]2
hT f a f b
1
12,
2kk h h
12 kkm
/2
( 1 ) ( 1 )
11
1
1 [ ( ( 2 1 ) )]
2
km
k k k k
j
T T h f a j h
( 1 ) ( 1 )
( ) ( 1 ) 1
141
jj
jj kk
kk j
TTTT
6 if,then
goto 8
end if
7 令 k=k+1,,goto 4
8 输出
9 end
( ) ( 1 )[]jjkkTT
1
1
2kkhh
jkT
7.5 Gauss求积公式
7.5.1 引言求积公式
(1)
当求积系数,求积节点 都可以自由选取时,其代数精确度最高可以达到多少次?
下面的引理可以回答上述问题,
1
( ) ( ) ( ) ( ) (,)
nb
kka
k
I f x f x dx A f x R f
1{}nkkA? 1{}nkkx?
引理 1 当求积系数 和求积节点都可以自由选取时,n点的求积公式 (1)的代数精确度最高可以达到 2n-1次,
证 假设求积公式 (1)具有 m次代数精确度,
即对任意的 m次代数多项式求积公式 (1)的精确成立,于是成立等式即若记 (2)
0
()
m
i
m i m
i
p x a x P
1
( ) ( ) ( )
mb
m k m ka
k
x p x dx A p x?
1
1 1 0
01
( ) (,.,)
mn b
i m m
i k m k m k ka
ik
a x x dx A a x a x a x a
( ),0,1,2,..,,b i i
a
x x d x i m
1{}nkkA? 1{}nkkx?
则 (2)式成为
(3)
由于系数 的任意性,故使 (3)式成为恒等式的充要条件是
(4)
(4)式的待定系数有 2n个,所以确定待定系数的
1 1 1 1 0 0...m m m ma a a a
10
1 1 1
...
n n n
m
m k k k k k
k k k
a A x a A x a A
1 1 0,,.,,,,mma a a a?
1 2 0
1 1 2 2 1
1 1 2 2
n
nn
m m m
n n m
A A A
A x A x A x
A x A x A x
独立条件至多给出 2n个,从而可知 m至多为
2n-1.
定义 1 n点的求积公式 (1)具有 2n-1次代数精确度 (或称为具有最高的代数精确度 )时,称为 Gauss型求积公式,
Gauss型求积公式的求积节点,称为
Gauss点,它们可以通过求区间 [a,b]上带权?(x)
的 n次正交多项式 的 n个根获得,所以先介绍正交多项式及其性质,然后讨论 Gauss型求积公式的构造,等等,
1{}nkkx?
()ngx
7.5.2 正交多项式及其性质定义 2 若
(1),则称函数 f(x)和 g(x)在区间
[a,b]上正交,
(2),则称函数 f(x)和 g(x)在区间 [a,b]上带权?(x)正交,
(3)代数多项式序列 (下标 k为多项式的次数,表示 k次多项式 ),在区间 [a,b]上满足当 m?n
当 m=n
则称多项式序列 为区间 [a,b]上带权
( ) ( ) 0ba f x g x d x
( ) ( ) ( ) 0ba x f x g x d x
0{ ( ) }kkgx
()kgx
2
0,
( ) ( ) ( )
( ) ( ( ) ) 0
b
bnm
a
na
x g x g x d x
x g x d x
0{ ( ) }kkgx
(x)的正交多项式序列,
定义 3 若 n次记多项式 中含 项的系数为,则称 为 的 首次系数 ; 时,称为 首次系数为 1的 n次多项式,
正交多项式有如下性质,
性质 1 若 是区间 [a,b]上带权?(x)的正交多项式序列,则它们线行无关,
证 对任意的 x?[a,b],若,在式子两边同乘?(x)g1(x)(L=0,1,..n),并从 a到 b积分,由的正交性定义 1中的 (3)可知必有
()ngx nx
nd nd
()ngx 0
nd?
* ()() n
n
n
gxgx
d?
0{ ( ) } nkkgx?
0
( ) 0
n
kk
k
c g x
0{ ( )} nkkgx? 0
lc?
L=0,1,..,n.故正交多项式序列 线性无关,
由性质 1可知,若 为 [a,b]上带权?(x)的正交多项式序列,则序列 可以作为空间的一组基函数,即 中的任一元素可由它们线性表出,
其中 为组合系数,
性质 2 若 为 [a,b]上带权的正交?(x)多项式序列,且,则
0{ ( )} nkkgx?
0{ ( )} nkkgx?
[,]nP a b [,]nP a b ()
npx
0
( ) ( )
n
n k k
k
p x a g x
0{}nkka?
0{ ( )} nkkgx?
0{ ( ) }kkgx
( ) [,]nq x P a b?
(1)
(2)
事实上,由性质 1,.由 的正交性定义容易证得 (1).证 (2)也是类似的,
为方便起见,记下面,不加证明地给出正交多项式如下的性质,
性质 3 [a,b]上带权函数?(x)的正交多项式序列 相邻三项的递推关系为其中,,
0
( ) ( )
n
kk
k
q x a g x
{ ( )}ngx
(,) ( ) ( ) ( )baf g x f x g x d x
0{ ( )}kkgx
1 1 1( ) ( ( ) ) ( ) ( ),1,2,..,,n n n n n ng x a x g x g x n
1n
n
n
d
a
d
1 (,)
(,)
n n n
n
n n n
d g x g
d g g?
( ) ( ) ( ) 0,1,2,..,b ka x q x g x d x k n n
( ) ( ) 0,0,1,2,...,1b ina x g x x dx i n
为 的首项系数,
即为性质 4 [a,b]上带权函数 的正交多项式序列 中任意相邻两个正交多项式 和的根相间,
若记,的根分别为,
则所谓 与 的根相同,即是指这两个正交多项式的根有如下的关系,
11
1 2
11
(,),
(,)
n n n n
n
n n n
d d g g
d g g?
11,,n n nd d d 11( ),( ),( )n n ng x g x g x
( ) (,1 )kg x k n n
0{ ( )}kkgx ()ngx
1 ()ngx?
()ngx 1 ()ngx? () 1{}nniix? ( 1) 11{}nnjjx
1 ()ngx?()ngx
()x?
性质 5 (1) 区间 [a,b]上带权函数?(x)的正交多项式序列 与 对应元素之间只相差一个比例常数,
(2)区间 [a,b]上带权函数?(x)首项系数为 1的正交多项式序列 唯一,
常见的正交多项式有 Legendre(勒让德 )多项式,Hermite多项式,Chebyshev多项式以及 Jacobi多项式,Chebyshev多项式在 5.2节已
( 1 ) ( ) ( 1 ) ( ) ( 1 )
1 1 2,1,2,...,1
n n n n n
i i i i ix x x x x i n
0{ ( ) }nnfx 0{ ( ) }nngx
* 0{ ( )}nngx
详细讨论,这里主要介绍 Legendre多项式一,Legendre多项式隐式表达式显式表达式
()nPx
0
2
( ) 1,
1 ( 1 )
( ),1,2,.,,
2!
nn
n nn
Px
dx
P x n
n d x
0
2
0
( ) 1,
1 ( 2 2 )
( ) ( 1 ),
2 ! ( ) ! ( 2 ) !
N
j n j
n n
j
Px
nj
P x x
j n j n j
1,2,...n?
当 n为偶数时,
当 n为奇数时,
Legendre多项式的主要性质有
(1)n次 Legendre多项式 的首项系数当 x=1,
当 x=-1.
(3)正交性为,为区间 [-1,1]上带权函数?(x)?1的正交多项式序列,且有
/ 2,
( 1 ) / 2,
n
N
n
其 中
()nPx
2
( 2 ) !( ) ;
2 ( ! )n n
nx
n?d
1,
( 2 ) ( 1 )
( 1 ),n n
P
0{ ( ) }nnPx
当 m?n
当 m=n
(4) Legendre多项式相邻三项的递推关系为二,Legendre多项式将隐式表达式
1
1
0,
(,) ( ) ( ) 2
,
21
m n m nP P P x P x d x
n
0
1
11
( ) 1,
( ),
21
( ) ( ) ( ),1,2,.,,
11
n n n
Px
P x x
nn
P x x P x P x n
nn
()Ln x
()() n n xx
n n
d x eL x e
dx
将隐式表达式中 n阶导数用乘积导数的
Leibniz公式可得显式表达式
Legendre多项式的主要性质有
(1)n次 Legendre多项式 的首项系数
(2)正交性为,为区间 [0,+?)上带权函数 的正交多项式序列,且有
0
!( ) ( 1 )
( ) !
n
n k k n k
nn
k
nL x C x
nk
()nLx ( 1 ) nnd
0{ ( ) }nnLx
() xxe
20
0,
(,) ( ) ( )
( ! ),
x
n m n mL L e L x L x d x n
当 mn
当 m=n
权函数 的正交多项式序列,且有
(3) Hermite多项式相邻三项的递推关系为四,Jacobi多项式
Jacobi多项式是在区间 [-1,1]上带权函数的正交多项式,其中?>-1,?>-1
2() xxe
2 0,,
(,) ( ) ( )
2 !,
x
m n m n n
mn
H H e H x H x d x
n?
当当 m=n;
0
1
11
( ) 1,
( ) 2,
( ) 2 ( ) 2 ( ),1,2,...n n n
Hx
H x x
H x x H x nH x n
( ) ( 1 ) ( 1 )x x x
(3) Legendre多项式相邻三项的递推关系为三,Hermite多项式表达式
Hermite多项式的主要性质有
(1)n次 Hermite多项式 的首项系数
(2)正交性为,为区间 上带
()nHx
()nHx
0
1
2
11
( ) 1,
( ) 1,
( ) ( 1 2 ) ( ) ( ),1,2,...n n n
Lx
L x x
L x n x L x n L x n
2
2 ()( ) ( 1 )
nx
nx
n n
deH x e
dx
( 2 ) nnd
0{ ( )}nnHx (,)
有的书籍文献把 Jacobi多项式记为即 n次 Jacobi多项式表示为其中 或,两种系数推出两种 Jacobi多项式,详细的情形请参阅文献 [27].
(,) ()nJx
(,) ( ) ( 1 ) ( 1 ) [ ( 1 ) ( 1 ) ]
n
nn
nn n
dJ x K x x x x
dx
1 2 1
( 2 ) ! 2n
nK
n
( 1 )
2!
n
n nK n
7.5.3 Gauss型求积公式由 7.5.1中的引理 1和定义 1可知 n点的求积公式 (1)若具有最高的代数精确度,或具有 2n-1
次的代数精确度成为 Gauss型求积公式,到底求积公式 (1)的求积节点 和求积系数如何选取,才能使之成为 Gauss型求积公式?
定理 7.5.1 求积公式 (1)中的 n个求积节点
,取在区间 [a,b]上带权函数?(x)的 n次正交多项式 的 n个根成为 Gauss型求积公式,
证 设,[a,b]上带权函数?(x)的
1{}nkkx? 1{}nkkA?
1{}nkkx?
()ngx
21( ) [,]nf x P a b
n次正交多项式 的 n个根记为,记的首项系数为,由定义 2有因此,
(5)
其中,在 (5)式两边同乘?(x),并从
a到 b积分,由正交多项式的性质可知,含项的积分为零,所以
(6)
注意到当 作为插值节点时建立的 n点插值
()ngx 1{}nkkx?
1{}nkkx?
nd
*
1
()( ) ( )n n
nk
k n
gxg x x x
d?
*( ) ( ) ( ) ( ),nf x q x g x r x
1( ),( ) [,]nq x r x P a b
* ( ) ( )ng x q x
( ) ( ) ( ) ( )bbaax f x d x x r x d x
求积公式至少具有 n-1次代数精确度,而,
所以
(7)
又由 (5)式可知,
即 (8)
综合 (6),(7),(8)式可知,当 时,求积公式 (1)
成立,
1
( ) ( )
n
n k k
k
I f A f x
1( ) [,]nr x P a b
1
( ) ( ) ( )
n
k k n
k
I r A r x I r
( ) ( ) ( 1,2,..,,)kkf x r x k n
( ) ( )nnI f I r?
21( ) [,]nf x P a b
1
( ) ( ) ( )
nb
kka
k
x f x d x A r x?
用 n点 Gauss求积公式
(9)
之值近似积分值 有下面的误差估计,
定理 7.5.2 若,则 Gauss型求积公式 (1)的误差估计 R(?,f)为其中证明略,
在稍后讨论 Gauss积分值数列的收敛性
1
( ) ( )
n
n k k
k
I f A f x
()nIf
2( ) [,]nf x C a b?
2
*2()(,) ( ) [ ( ) ]
( 2 ) !
n b
na
fR f x g x dx
n
*
1
( ) ( )
n
nk
k
g x x x
等问题时,需要用到 Gauss型求积公式的求积系数 大于零的结论,这里用下面的定理给出,
定理 7.5.3 Gauss型求积公式的求积系数大于零,
证 令,这里 为区间 [a,b]
上带权函数?(x)的 n次正交多项式 的 n个根,显然
1{}nkkA?
1{}nkkx?
()ngx
1{}nkkA?
2
1()
n
i
i
k
xx
fx
xx
1,
0,
() ()
n
j
ki
i
ik
fx xx
当 jk
,当 j=k
由于,所以对 f(x)求积公式 (1)精确成立,即因为所以在 7.2节,我们讨论了复化梯形求积公式和复化 Simpson求积公式的收敛性,那么 Gauss
1
( ) ( ) ( ) ( ) ( )
nb
n j ja
j
I f x f x dx I f A f x?
2
1,
()
n
k k i
i
ik
A x x
,
2
1,
( ) 0,( ) ( ) 0
n b
ki a
i
ik
x x x f x d x?
0,1,2,3,.,,,kA k n
22() nf x P
型求积公式被积函数 f(x)应当满足什么条件才收敛呢? Gauss型求积公式的收敛性问题由下面的定理给出,
定理 7.5.4 若 f(x)?[a,b],则 Gauss型求积公式所求积分值序列 收敛于积分值 I(f),即证 因为 f(x)?[a,b],由 Weierstrass定理对任意的?1>0,存在,使得
1
{ ( ) ( ) }
n
n k k
k
I f A f x
1
l im ( ) ( ) ( )
n b
kk an
k
A f x x f x d x?
( ) [,]mmp x P a b?
(10)
对任意的 x?[a,b]成立,
由于公式 (1)为 Gauss型求积公式时具有
2n-1次代数精确度,取 N>(m+1)/2,故当 n>N时,
即 m<2N-1<2n-1时,有
(11)
成立,于是由 (11)式可知,而由 (10)式,有
1( ) ( )mf x p x
( ) ( )n m mI p I p?
( ) ( ) ( ) ( ) ( ) ( )n m m n mI f I f I f I p I p I p
( ) ( )n m nI p I f
( ) ( ) 0m n mI p I p
和,从而因为,记
( ) ( ) ( ) [ ( ) ( ) ]bmm
a
I f I p x f x p x d x
1
1
n
k
k
A?
1( ) ( )k m kf x p x
1
( ) ( ) [ ( ) ( ) ]
n
n m n k m k k
k
I p I f A p x f x
1
1
n
k
k
A?
0,1,2,..,,kA k n
故即
1 1 1
11
( ) ( ) 0 2
nn
n k k
kk
I f I f A A C
l i m ( ) ( ),nn I f I f
1
()
nb
ka
k
C x dx A?
7.5.4 Gauss型求积公式的构造与应用定理 7.5.1实际上给出了构造 Gauss型求积公式的一种方法,这种方法,当给定了积分区间
[a,b]和权函数?(x)以后构造 n个点的 Gauss型求积公式,先求出区间 [a,b]上带权函数?(x)的 n次正交多项式,然后用多项式求根的方法求出的 n个根,从而获得了求积节点为了求得求积系数,将 n个求积节点代入方程组 (4)中的前 n个方程并加以求解,
即解线性代数方程组
()ngx
()ngx
1{}nkkx? 1{}nkkx?
1{}nkkA?
1{}nkkx?
求得求积系数,完成 Gauss型求积公式的构造,
表 7-5为 Gauss-Legendre求积公式的求积系数 和求积节点 的一个表,而表 7-6和表
7-7则分别是 Gauss-Legendre和 Gauss-Hermit
-te求积公式的求积系数 和求积节点 的一
1{}nkkA?
1 2 0
1 1 2 2 1
1 1 2 2
n
nn
m m m
n n m
A A A
A x A x A x
A x A x A x
kA
kx
kxkA
个表,有些数值分析的书籍还给出了 Gauss-
Chebyshev求积公式的求积节点与求积系数,
表见 P366.
需要指出的是,不同求积公式的求积系数与求积节点,积分区间和权函数是不同的,Gauss-
Lagendre求积公式的积分区间为 [-1,1],权函数
(x)=1.而 Gauss-Legendre求积公式的积分区间为 [0,+?),权函数,Gauss-Hermitte
求积公式的积分区间为 (-?,+?),权函数
Gauss-Chebyshev求积公式的积分区间为 [-1,1]
() xxe
2() xxe
权函数这里需要指出的另一点是 Gauss-
Chebyshev求积公式的求积系数是相同的,例如,n点的 Gauss-Chebyshev求积公式,它的 n个求积系数 都是,即,而 n个求积节点则为正是因为 Gauss-Chebyshev求积公式的求积系数相同,所以在实际计算时,乘法的次数只需一次,节省了 n-1次的乘法运算,
2
1()
1x x
kA
n
12 nA A A n
21c o s,1,2,.,,,
2k
kx k n
n?
例 1 求 使求积公式具有三次代数精确度,
问题是构造区间 [0,1]上带权函数的两点 Gauss型求积公式,
解 方法 1 容易计算出当时 的积分值分别为,所求公式具有 3次代数精确度,故可得为未知数的方程组为
1 2 1 2,,,,A A x x
1
1 1 2 20 ( ) ( ) ( ) (,)x f x d x A f x A f x R x f
()xx
23( ) 1,,,f x x x x?
1
0 ()x f x dx?
2222,,,,
3 5 7 9
1 2 1 2,,,,A A x x
(1)
(2)
(3)
(4)
又因为 为 Gauss型求积公式的求积节点,
所以它们是区间 [0,1]上带权函数 且首项系数为 1的二次正交多项式 的两个根,
不妨记,为此又因为 必须满足方程 (2)(3)(4),所以
12
1 1 2 2
22
1 1 2 2
33
1 1 2 2
2 / 3
2 / 5
2 / 7
2 / 9
AA
A x A x
A x A x
A x A x
12,xx
()xx
*2 ()gx
*22 ()g x x px q
**
2 1 2 2( ) ( ) 0g x g x
12,xx
由可得关于 p,q的方程组为解此方程组得
**
1 2 1 2 2 2
**
1 1 2 1 2 2 2 2
( 1 ) ( 2 ) ( 3 ) ( ) ( ) 0
( 2 ) ( 3 ) ( 4 ) ( ) ( ) 0
q p A g x A g x
q p A x g x A x g x
2 2 2
3 5 7
2 2 2
5 7 9
qp
qp
10
9
5
21
p
q
将所求 p,q代入,求得其根为再将所求 代入方程 (1)(2),联立解得为此,公式为所求具有 3次代数精确度的求积公式,
*22 ()g x x px q
1
2
0,2 8 9 9 4 9 1 9 8
0,8 2 1 1 6 1 9 1 2
x
x
12,xx
1
0 ( ) 0,2 7 7 5 5 5 9 9 7 ( 0,2 8 9 9 4 9 1 9 8 )x f x d x f
1
2
0,2 7 7 5 5 5 9 9 7
0,3 8 9 1 1 0 6 6 9
A
A
0.389110669 ( 0.821161912) (,)f R x f
方法 2 求出区间 [0,1]上带权函数的二次正交多项式,并求出其根,获得求积节点,再求方程组 (12)前两个方程组成的方程组获得求积系数,
由于 和 都是的线性无关组,所以考虑由 求得所需的正交多项式,这种方法可以称作将它们正交化,
记
()xx
*2()gx
12,xx
12,AA
***0 1 2{ ( ),( ),( ) }g x g x g x
2{1,,}xx
2 [0,1]P
1 1 / 2
20(,) ( ) ( ),(,)f g x f x g x d x f f f
令所以又令 这样即如此作出的 与 正交,也与 正交,
由于所以这样,
* * * *0 0 0 0
2( ) 1,(,) 2 / 3,2 / 3g x g g g
*1 0 0( ) (,),g x x x e e
*
0
0 *
02
() 3,
2()
gxe
gx
*1 0 0 0 0 0 0( ( ),) ( (,),) (,) (,) 0g x e x x e e e x e x e
*1 ()gx
0e
0 ()gx
1
0 0
3 2 3(,)
2 5 2x e x x d x
**
11 2
38( ),( )
5 175
g x x g x
同理可令并且容易验证,与 正交,容易求得所以
*
1
1 *
1 2
() 17 5 3 17 5
8 5 8()
gx
ex
gx
* 2 2 22 1 1 0 0( ) (,) (,),g x x x e e x e e
*2 ()gx 10,ee
2
1
1 7 5 1 6(,)
8 3 1 5xe 2 0 32(,),
27
xe
*2
2
1 0 5()
9 2 1g x x x
第二种方法与第一种方法求出的 是一样的,后面的求解过程相同,这里略去,
方法 2是一种将一组线性无关函数组正交化而得到正交多项式 的方法,高于 2次的正交多项式用这样的方法同样可以得到,这样求正交多项式在实际应用时是方便的,
这里需要附带说明的是,Gauss-Legendre
,Gauss-Chebyshev求积公式作数值求积精度不够时,可以采取将积分区间 [a,b]若干等分后,
将每一个子区间映射到区间 [-1,1]上再用相同
*2 ()gx
*2()gx
2{1,,}xx
节点数的求积公式进行数值求积的计算,通常会得到精度较好的计算结果,
Gauss型求积公式具有数值结果精度高,收敛得以保证、计算简便、易于在计算机上实现等优点,并且在积分区间 [a,b]有限时便于推广到高维数值积分,不足之处是公式的构造比较困难,另一个是由于相邻次数的正交多项式的根,从而造成增加求积节点以提高计算结果的精度时,原先所有求积节点上的函数值全部无用,所以在具体应用 Gauss求积公式计算数值积分时,n取得都较小,
