高斯求积公式-数值分析课程设计2

  • 格式:doc
  • 大小:291.50 KB
  • 文档页数:14

一、 引言介绍高斯型求积公式,并使用其求积分⎰=1sin I xdx 。

要求:数值实验结果要体现出随高斯点的增加误差的变化。

我们知道,求积公式⎰∑=≈bani i ix f Adx x f 0)()( (1.1)含有22+n 个待定常数i x 及),,2,1,0(n i A i =,如果它具有n 次代数精确度,则它应使1+m 个方程mk dx x x A bakni ki i ,,2,1,0,==⎰∑= (1.2)精确成立。

作为插值型求积公式(1.1)它至少具有n 次代数精确度;另一方面,令)())(()(101n n x x x x x x x ---=+ ω,则对22+n 次多项式)()(21x x f n +=ω而言,(7.5.1)右端为零,而左端严格大于零,即(7.5.1)式对22+n 次多项式)(21x n +ω不准确成立。

但要确定方程组(7.5.2)中的22+n 个待定常数i x 与i A ,最多需要给出22+n 个独立条件,所以m最大取12+n 。

因此,插值型求积公式(1.1)的代数精确度最小是n ,最大是12+n .由此可见,高斯公式的代数精度比牛顿-科特斯公式高,求解高斯求积公式的关键就是解出上述2n+2个待定常数。

为解决上述问题,首先要先给出三个定理:定理一:以n x x x ,,,10 为节点的插值型求积公式(7.5.1)具有12+n 次代数精确度的充要条件是以这些节点为零点的多项式)())(()(101n n x x x x x x x ---=+ ω与任意次数不超过n 的多项式)(x P 均在区间],[b a 上正交,即⎰=+ban dx x x P 0)()(1ω (1.3)定理二:高斯公式(1.1)的求积系数k A 全为正,且nk dx x l dx x l A bak bak k ,1,0,)()(2===⎰⎰(1.4)定理三:对于高斯公式(1.1),其余项为dxx fn f R ban n ⎰+++=)()()!22(1)(21)22(ωη (1.5)其中).())(()(],,[101n n x x x x x x x b a ---=∈+ ωη证明 以n x x x ,,,10 为节点构造)(x f 的埃尔米特插值多项式)(x H),()(i i x f x H = ni x f x H i i ,1,0),()(='='因为)(x H 是12+n 次多项式,而它的余项是)()()!22(1)()(21)22(x fn x H x f n n +++=-ωξ所以高斯公式(7.5.1)对)(x H 能准确成立,即∑∑⎰====ni i in i iibax f Ax H A dx x H 0)()()(从而dxx fn dxx H dx x f x f A dx x f f R n ban babani i i ba)()()!22(1)()()()()(21)22(0++=⎰⎰⎰∑⎰+=-=-=ωξ若)()22(x fn +在区间],[b a 上连续,由于)(21x n +ω在],[b a 上不变号,故应用积分中值定理可得],[,)()()!22(1)(21)22(b a dx x fn f R ban n ∈+=⎰++ηωη上述定理说明,与牛顿—科兹公式进行比较,高斯公式不但具有高精度,而且它还是数值稳定的,但是节点和求积系数的计算比较麻烦。

注:由于篇幅有限以及定理三的重要性,故略去定理一、二的证明。

二、 方法描述:2.1、高斯—勒让德(Gauss-Legendre )公式对于任意求积区间],[b a ,通过变换ta b b a x 22-++=可化为区间]1,1[-,这时⎰⎰--++-=11)22(2)(dtt a b b a f a b dx x f ba因此,不失一般性,可取,1,1=-=b a 考查区间]1,1[-上的高斯公式∑⎰=-≈ni i ix f Adx x f 011)()( (2.1)我们知道,勒让德(Legendre)多项式])1[()!1(21)(121111+++++-+=n n n n n x dxdn x L (2.2)是区间]1,1[-上的正交多项式,因此, )(1x L n +的1+n 个零点就是高斯公式(2.1)的1+n 个节点。

特别地,称)(1x L n +的零点为高斯点,形如(2.1)的高斯公式称为高斯-勒让德公式。

利用勒让德多项式的一个性质)]()()[1()()1(112x xL x L n x L x n n n++-+='-可得,高斯-勒让德求积系数i A 为ni x L n x A i n i i ,2,1,0,)]()1[()1(222=+-=(2.3)按(1.5)式,可推得其余项为)(])!22)[(32()]1[(2)()22(3432η+++++=n n fn n n f R (2.4)若取x x L =)(1的零点00=x 为节点,则2)]0([)01(2200=-=L A从而一点高斯-勒让德公式(中矩形公式)为 )0(2)(11f dx x f ≈⎰- (2.5)其余项为)(31)(ηf f R ''=若取)13(21)(22-=x x L 的两个零点31±为节点,则1)]31(2[])31(1[22120=---=L A1)]31(2[])31(1[22121=-=L A从而二点高斯-勒让德公式为)31()31()(11f f dx x f +-≈⎰- (2.6)其余项为)(1351)(24522)()4()4(345ηηfff R =⋅⋅=同理,三点高斯-勒让德公式为)515(95)0(98)515(95)(11f f f dx x f ++-≈⎰- (2.7)其余项为)(157501)()6(ηff R =一般地,高斯-勒让德公式(2.1)的节点可以通过勒让德多项式的零点确定,而求积系数通过(2.3)式确定。

表2-1给出了高斯—勒让德公式在节点数为6,5,4,3,2,1时的节点、求积系数及余项。

表2-1例 2.1 用二点高斯-勒让德公式计算积分⎰=2sin πxdxI解 作变量代换),1(4+=t x π则⎰⎰-+==1124)1(sin4sin dtt xdx I πππ记4)1(sin)(+=t t f π,因为节点5773503.0±=i t 得,32589.0)(0=t f 94541.0)(1=t f所以,由二点高斯公式94541.0)94541.032589.0(4)]()([104=+=+≈ππt f t f I计算结果比用复合梯形公式7个节点计算的结果还要好。

2.2、高斯-切比雪夫(Gauss-Chebyshev)求积公式区间为[-1,1],权函数的Gauss型求积公式,其节点是Chebyshev多项式的零点,即,而,于是得到(2.8)称为Gauss-Chebyshev求积公式,公式的余项为(2.9)这种求积公式可用于计算奇异积分.例2.2用三点和四点Gauss-Chebyshev求积公式计算积分,并估计误差.解这里,由Gauss-Chebyshev求积公式(5.4.9)可得当n=2时,,求得代入上式得估计误差可用余项表达式(2.8),因,故当n=3时,,求得误差小结:Gauss型求积公式是上带权的求积公式,它具有最高代数精确度2n+1,实际上由于求积系数及节点都是待定系数,它共有2n+2个,可使(5.4.1)对任何2n+1次多项式精确成立。

具有2n+1次代数精确度的求积公式节点就是Gauss点,实际上它就是在上带权正交多项式的零点。

得到求积节点以后,同样可利用(5.4.1)对精确成立,得到关于的线性方程组解此方程组得到的求积公式系数,它是稳定的,也是收敛的,具有较高的精度。

通常使用的具体公式是Gauss-Legendre求积公式(简称Gauss求积公式),它是区间为,权函数为的公式,当n=1时可得,比n=2(三点)simpson公式好,当n=2时可得比n=4(五点)的Cotes公式好,而计算量却减少。

另一个Gauss型求积公式时Gauss-Chebyshev求积公式,它除了精度高,还可计算反常积分,如例2.2。

三、数值实验为了观察高斯积分公式随着高斯点的增加积分值变化的规律,本试验利用MA TLAB数学软件,选取了27组数据,逐个算出⎰=10sinI xdx的积分值,且依次算出了各个不同节点处的误差。

数据表格如下:由上表可以得出这样的结论:当高斯点数N为奇数是,按照高斯求积公式算出来的结果和真实值相差较大,并且随着N的增加误差逐渐减小,但是收敛效果不明显;当N为偶数时,按照高斯求积公式算出来的结果和真实值相差不大,并且随着N的增加误差逐渐减小,且收敛性明显,如表所示,当N仅仅为8时,误差为610 数量级了。

四、参考文献[1] 李庆杨王能超易大义.数值分析(第4版).华中科技大学出版社.2006年7月[2] 邓建中主编. 计算方法(第二版)[M]. 交通大学出版社,2001.[3] 刘琼荪编. 数学实验[M]. 北京:高等教育出版社,2004年7月[4] 王建卫曲中水凌滨编著.MATLAB 7.X程序设计.中国水利水电出版社,2007年9月[5] 王沫然编著MATLAB与科学计算(第二版)电子工业大学出版社附录:此课程设计的部分原程序代码:程序部分:function s=guassl(a,b,n)h=(b-a)/n;s=0.0;for m=0:(1*n/2-1)s=s+h*(guassf(a+h*((1-1/sqrt(3))+2*m))+guassf(a+h*((1+1/sqrt(3))+2*m) ));endsI=int('sin(x)',0,1);c=(I-s)/I;d=vpa(c,10)%Îó²î%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%function y=guassf(x)y=sin(x);输出结果:>> guassl(0,1,2)s =0.459587812395265d =.23903042803558660869976310583379e-3ans =0.459587812395265>> guassl(0,1,2)s =0.459587812395265d =.2390303484e-3ans =0.459587812395265>> guassl(0,1,3)s =0.214102809190798d =.5342530277ans =0.214102809190798>> guassl(0,1,4)s =0.459690990276648d =.1458306206e-4ans =0.459690990276648>> guassl(0,1,5)s =0.303291484191192d =.3402370991ans =0.303291484191192>> guassl(0,1,6)s =0.459696375761941d =.2867754215e-5ans =0.459696375761941>> guassl(0,1,7)s =0.345399399138380d =.2486379559ans =0.345399399138380。