数值计算方法-2
- 格式:pdf
- 大小:715.54 KB
- 文档页数:66
毕业论文文献综述信息与计算科学定积分的数值计算方法一、 前言部分在科学与工程计算中,经常要计算定积分()()().baI f f x dx a b =-∞≤≤≤∞⎰ (1.1)这个积分的计算似乎很简单,只要求出f 的原函数F 就可以得出积分(1.1)的值,即()()().I f F b F a =- (1.2)如果原函数F 非常简单又便于使用,那么式(1.2)就提供了计算起来最快的积分法.但是,积分过程往往将导出新的超越函数,例如,简单积分1dx x ⎰可引出对数函数,它已不是代数函数了;而积分2x edx -⎰,将引出一个无法用有限个代数运算、对数运算或指数运算组合表示的函数.有些积分虽然容易求解,并且原函数仍然是一个初等函数,但可能过于复杂,以致于人们采用(1.2)来计算之前还得三思而行[1].例如411dx C x =++⎰, (1.3) 采用式(1.3)这种“精确”表达式时,所需运算次数是个根本问题.由式(1.3)看出,需计算对数和反正切,因此只能计算到一定的近似程度.因此可以看出,这类表面上是“精确”的方法,实际上也是近似的.因此,我们常常需要探讨一些近似计算定积分的数值方法[2].通过人们的研究和发现,得出了很多数值计算的方法,比如利用牛顿-科茨求积公式,复合求积公式,龙贝格积分法,高斯求积公式,切比雪夫求积法等来解决定积分的数值计算问题.构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数,由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-柯茨公式,例如梯形公式与抛物线公式就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中,对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法,它具有公式简练、计算结果准确、使用方便、稳定性好等优点,因此在等距情形宜采用龙贝格求积公式.当用不等距节点进行计算时,常用高斯型求积公式计算,它在节点数目相同情况下,准确程度较高,稳定性好,而且还可以计算无穷积分[3].二、 主题部分2.1 牛顿-科茨求积公式[4]2.1.1 公式的一般形式[4]将积分(1.1)中的积分区间[],a b 分成n 等分,其节点k x 为1,()k x a kh h b a n=+=- (0,1,,)k n =L . 对于给定的函数f ,在节点k x (0,1,,)k n =L 上的值()k f x 为已知.那么f 在n+1个节点01,,,n x x x L 上的n 次代数插值多项式为00()().n nj n kk j k j j k x x p x f x x x ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏ 如果记x a th =+,则上式可以写为00()().n nn kk j j k t j p x f x k j ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏ (2.1) 在积分(1.1)中的被积函数f 用其n+1个节点的代数插值多项式()n p x 来代替,可 得 ()()()()bbn n aaI f f x dx I f p x dx =≈=⎰⎰.多项式的积分是容易求出的,因此把上式写为()()()nn n k k I f I f A f x =≈=∑, (2.2)其中 ()00(),n n n k k j j kb a t j A dt b ac n k j=≠--==--∏⎰ (2.3) ()00(1)().!()!n kn n n kj j kct j dt k n k n -=≠-=--∏⎰ (2.4) 公式(2.2)称为牛顿-科茨求积公式或称为等距节点求积公式,k A 称为求积公式系数,()n k c 称为科茨求积系数.牛顿-科茨求积公式的误差估计()n E f ()()n I f I f =-,由下面定理给出 定理2.1 (1) 如果n 为偶数,(2)n f +在[],a b 上连续,则有[]3(2)()(),,n n n n E f c hf a b ηη++=∈, (2.5)其中 201(1)(2)()(2)!n n c t t t t n dt n =---+⎰L . (2) 如果n 为奇数,(1)n f+在[],a b 上连续,则有[]2(1)()(),,n n n n E f c h f a b ηη++=∈, (2.6)其中 01(1)(2)()(1)!n n c t t t t n dt n =---+⎰L . 定义2.1 如果求积公式()()nbk k ak f x dx A f x =≈∑⎰对所有次数不高于n 的代数多项式等式精确成立,但存在n+1次的代数多项式使等式不成立,则称上式求积公式具有n 次代数精度.由定理2.1可知,牛顿-科茨求积公式(2.2)的代数精度至少是n 次,而当n 是偶数时,(2.2)的代数精度可达n+1次.2.1.2 梯形公式[5]在牛顿-科茨公式(2.2)中,取n=1时(1)(1)011,2c c ==所以有 []1()()()().2b aI f I f f a f b -≈=+ (2.7) 公式(2.7)称为梯形公式,如果用连接(),()a f a 和(),()b f b 的直线来逼近f ,并对这线性函数进行积分可得到1()I f .再用1()I f 来逼近()I f . 定理 2.2 若[]2,f Ca b ∈,则梯形公式(2.7)的误差为[]3111()()()()''(),,.12E f I f I f b a f a b ηη=-=--∈ 2.1.3 辛普森公式[6]在牛顿-科茨公式(2.2)中,取n=2,则有220011(1)(2),46c t t dt =--=⎰221014(2),26c t t dt =--=⎰ 222011(1),46c t t dt =-=⎰有此得到2()()()4()().32h a b I f I f f a f f b +⎡⎤≈=++⎢⎥⎣⎦(2.8) 其中1()2h b a =-.式(2.8)称为辛普森公式. 定理2.3 若[]4,f Ca b ∈,则辛普森公式(2.8)的误差为[]5(4)221()()()(),,.90E f I f I f h f a b ηη=-=-∈2.2 复化求积公式[7]上面已经给出了计算积分()()baI f f x dx =⎰的3个基本的求积公式:梯形公式,辛普森公式,牛顿-科茨公式,并给出了它们误差的表达式.由这些表达式可知其截断误差依赖于求积区间的长度.若积分区间的长度是小量的话,则这些求积公式的截断误差是该长度的高阶小量.但若积分区间的长度比较大,直接使用这些公式,则精度难以保证.为了提高计算积分的精度,可把积分区间分为若干个小区间,()I f 等于这些小区间上的积分和,然后对每个小区间上的积分应用上述求积公式,并把每个小区间上的结果累加,所得到的求积公式称为复化求积公式.将积分区间[],a b 作n 等分,并记,,0,1,,k b ah x a kh k n n-==+=L ,于是 11()()k kn x x k I f f x dx +-==∑⎰.2.2.1 复化梯形求积公式[8]如果需要求出一个已知函数()f x 在一个很大区间[],a b 上的积分,那么我们可以把区间分成n 个长度为x h ∆=的小区间,对每一个小区间用梯形法则,然后再把这些小区间上的积分值相加.于是就得到了计算定积分的复化梯形公式:1101210()()(222)22n bi i n n ai h hf x dx f f f f f f f -+-=≈+=+++++∑⎰L (2.9)整体积分误差等于n 个小区间上的积分误差之和:整体误差= []312''()''()''()12n h f f f ξξξ-+++L ,其中i ξ是第i 个小区间上的某一点.如果''()f x 在区间[],a b 上连续,那么由连续函数的性质可知,在区间[],a b 上存在点ξ使得''()i f ξ的平均值等于()f ξ.于是由于nh b a =-,有整体误差= 322''()''()()1212nh b a f h f O h ξξ--=-=, 局部误差是3()O h ,整体误差是2()O h .2.2.2 复化辛普森求积公式[9]对于积分()baf x dx ⎰,将[],a b 等分,每个小区间长度b ah n-=,节点记为 (0,1,2,,)k x a kh k n =+=L ,第k 个小区间记为[]1,(1,2,,)k k x x k n -=L .记[]1,k k x x -的中点为1121()2k k k xx x --=+,则复化辛普森公式为 1112()()()4()()6n bk k ak k h f x dx S h f x f x f x --=⎡⎤≈=++⎢⎥⎣⎦∑⎰.2.3 龙贝格积分[10]现在要介绍用龙贝格(Romberg )命名的一个算法,龙贝格首先给出了这种算法的递推形式,假设需要积分()baI f x dx =⎰ (2.10)的近似值.在讨论过程中函数()f x 和区间[],a b 将保持不变.2.3.1 递推梯形法则[10]设()T n 表示在长度是()/h b a n =-的n 个子区间上积分I 的梯形法则.根据()''()nbai f x dx h f a ih =≈+∑⎰,我们有 00()()''()''()nn n i i b a b a T h f a ih f a i n n ==--=+=+∑∑, (2.11) 这里求和符号中的两撇表示和式中第一项和最后一项减半. 2.3.2 龙贝格算法[10]在龙贝格算法中使用上述公式.设(,0)R n 表示具有2n个子区间的梯形估计,我们有[]1211(0,0)()()()21(,0)(1,0)((21))2n n n i R b a f a f b R n R n hf a i h -=⎧=-+⎪⎪⎨⎪=-++-⎪⎩∑ , (2.12) 对于一个适度的M 值,计算(0,0),(1,0),(2,0),,(,0)R R R R M L ,并且其中没有重复的函数值的计算.在龙贝格算法的其余部分中,还要计算附加值(,)R n m .所有这些都可以被理解为积分I 的估计.计算出(,0)R M 后,不再需要被积函数f 值的计算.根据公式[]1(,)(,1)(,1)(1,1)41m R n m R n m R n m R n m =-+-----, (2.13)对于1n ≥和1m ≥构造R 阵列的各列.定理 2.4(龙贝格算法收敛性定理)[10]若[],f C a b ∈,则龙贝格阵列中每一列都收敛于f 的积分.因此,对每个m ,lim (,)()baR n m f x dx =⎰.2.4高斯求积[11]前面研究的求积公式都是事先确定了n 个节点,然后按使求积公式阶数达到最大的原 则选取最佳权.由于自由参数为n 个,所以阶数一般为n-1,但如果节点的位置也自由选择,则自由参数的个数将变为2n ,因此求积公式的阶数可达到2n-1.高斯求积公式就是通过选择最佳的节点和权,使求积公式的阶数最大化.一般地,对每个n ,n 点高斯公式都是唯一的,而且阶数为2n-1.因而,对一定的节点个数,高斯求积公式的精度是最高的.但它的求得比牛顿—柯特斯公式要困难得多.虽然它的节点和权也可由待定系数法确定,但得到的方程是非线性的.2.4.1 高斯求积公式[11]为说明高斯求积公式,推导区间[]1,1-上的两点公式1112221()()()()()I f f x dx w f x w f x G f -=≈+=⎰,其中的节点1x 、2x 及权1w 、2w 按使求积公式阶数最大化的原则选取.令公式对前四个单项式精确成立,得力矩方程组112111122112221122113331122112,0,2,30.w w dx w x w x xdx w x w x x dx w x w x x dx ----⎧+==⎪⎪+==⎪⎪⎨⎪+==⎪⎪+==⎪⎩⎰⎰⎰⎰这个非线性方程组的一个解为12121,1,x x w w =-===另一个解可通过改变1x ,2x 的符号而得到.这样,两点高斯求积公式为2()(G f f f =-+,阶数为3.另外,高斯求积公式的节点也可以由正交多项式得到.若p 是n 次多项式,且满足()0,0,,1,bk ap x x dx k n ==-⎰L 则p 与[],a b 区间上所有次数小于n 的多项式正交,容易证明:1. p 的n 个零点都是实的、单的,且位于开区间(,)a b .2. 区间[],a b 上以p 的零点为节点的n 点插值型求积公式的阶数为2n-1,是唯一的n 点高斯公式.定义2.2[12] 如果1n +个节点的求积公式()()()nbk k ak x f x dx A f x ρ=≈∑⎰(2.14)的代数精度达到21n +,则称式(2.14)为高斯型求积公式,此时称节点k x 为高斯点,系数k A 称为高斯系数.定理2.5[12] 以01,,,n x x x L 为高斯点的插值型求积公式具有21n +次代数精确度的充要条件是以这些节点为零点的多项式101()()()()n n x x x x x x x ω+=---L与任意次数不超过n 的多项式()p x 带权()x ρ均在区间[],a b 上正交,即1()()()0bn ax p x x dx ρω+=⎰. (2.14)定理2.6 高斯公式()()nbi i ai f x dx A f x =≈∑⎰(2.15)的求积系数k A 全为正,且 2()(),0,1,,bbk k k aaA l x dx l x dx k n ===⎰⎰L . (2.16)定理2.7 对于高斯公式(2.14),其余项为 (22)211()()()()(22)!b n n a R f f x x dx n ξρω++=+⎰ , (2.17) 其中[]101,,()()()().n n a b x x x x x x x ξω+∈=---L2.4.2 高斯—勒让德(Gauss-Legendre )公式[13] 对于任意求积区间[],a b ,通过变换22a b b ax t +-=+,可化为区间[]1,1-,这时11()()222bab a a b b af x dx f t dt --+-=+⎰⎰. 因此,不失一般性,可取1,1a b =-=,考查区间[]1,1-上的高斯公式 11()()ni i i f x dx A f x -==∑⎰. (4.5)我们知道,勒让德多项式1211111()(1)2(1)!n n n n n d L x x n dx+++++⎡⎤=-⎣⎦+, (4.6) 是区间[]1,1-上的正交多项式,因此,1()n L x +的n+1个零点就是高斯公式(4.5)的n+1个节点.特别地,称1()n L x +的零点为高斯点,形如(4.5)的高斯公式称为高斯—勒让德公式.以上这些公式中的节点和求积系数可查表得到. 2.4.3 高斯—哈米特求积公式(Gauss-Hermite )[14] Gauss-Hermite 求积公式2()0()()nx n k k k ef x dx f x ω∞--∞=≈∑⎰, (4.7)其余项为(22)1(().2(22)!n n n n R f f n ξ+++=+ (4.8)2.4.4 高斯—切比雪夫(Gauss-Chebyshev )求积公式[15] 区间为[]1,1-,权函数()x ρ=Gauss 型求积公式,其节点k x 是Chebyshev多项式1()n T x +的零点,即21cos (0,1,,)2(1)k k x k n n π⎡⎤+==⎢⎥+⎣⎦L ,而(0,1,,)1k A k n n π==+L于是得到1021cos 12(1)nk k f n n ππ-=⎡⎤+≈⎢⎥++⎣⎦∑⎰(4.9) 称为Gauss-Chebyshev 求积公式,公式的余项为 (22)2(1)2()(),(1,1)2(22)!n n n R f f n πηη++=∈-+ , (4.10) 这种求积公式可用于计算奇异积分.2.5 递推型高斯求积[10]高斯求积公式不具有递推性:当节点个数一定时,如果自由选择所有的节点和权以达到最高的阶数,则节点个数不同的公式一般没有公共节点,这意味着与一组节点对应的积分值,在用另外一组节点计算积分值时不能被利用.Kronrod 求积公式避免了这种工作量的增加,这类公式是对称的,n 点高斯公式n G 与2n+1个点Kronrod 公式21n K +对应.21n K +节点的约束条件为:以n G 的节点作为21n K +的节点,按求积公式达到最高阶数的要求确定21n K +中剩下的n+1个节点及2n+1个权(其中包括n G 的节点的权).这样,求积公式的阶数可达到3n+1,而真正2n+1个点高斯公式应该是4n+1阶的,所以精度和效率是一对矛盾.使用两个节点个数不同的求积公式的主要原因是可以用它们的差估计积分近似值的误差.使用Gauss-Kronrod 公式对时,若以21n K +的值作为积分的近似值,则一半基于理论,一半基于经验,可以得到关于误差的保守估计: 1.521(200)n n G K +-.Gauss-Kronrod 公式不仅有效地提供了较高的精度,还给出了可靠误差估计,所以它被认为是最有效的求积公式之一,并且构成了主要软件库中求积程序的基础,特别地,公式715(,)G K 已被普遍使用.三、 总结部分因为一些定积分的求解比较复杂,所以数值积分的理论与方法一直是计算数学研究的基本课题.各种定积分的数值计算方法的出现和发展,加快和简化了求解定积分的效率和步骤.以上主要介绍了各种数值积分的方法——牛顿-科茨求积公式,复合求积公式,龙贝格积分法,高斯求积公式等.每种方法都有各自的优缺点,针对不同的积分函数采用不同的方法,所以在实际计算时,要做适当的采取.相信随着理论分析和研究的日益深入,求定积分的数值计算方法将更加简单和完善,为我们的计算带来前所未有的方便,在数学领域也将会更上一层楼.四、参考文献[1] 孙志忠,吴宏伟,袁慰平,闻震初.计算方法与实习(第4版)[M].南京:东南大学出版社,2009,(2): 128~129.[2]Micheal T .Heath . 张威,贺华,冷爱萍译.科学计算导论(第2版)[M].北京:清华大学出版社,2005,(10): 396~297.[3]李桂成.计算方法[M].北京:电子工业出版社,2005,(10):186.[4] 现代应用数学手册编委会. 现代应用数学手册——计算与数值分析卷[M]. 北京:清华大学出版社,2005,(1): 163~168.[5] 林成森. 数值计算方法(上)[M]. 北京:科学出版社,2004,(5): 220~221.[6]冯康.数值计算方法[M].北京:国防工业出版社,1978,(12): 45~47.[7]孙志忠,袁慰平,闻震初.数值分析(第2版)[M].南京:东南大学出版社,2002,(1): 191~194.[8] (美)柯蒂斯F .杰拉尔德 帕特里克O .惠特莱. 应用数值分析(第7版)[M].北京:机械工业出版社,2006,(8): 222~225.[9]夏爱生,胡宝安,孙利民,夏凌辉.复化Simpson 数值求积公式的外推算法[J].军事交通学院学报.2006,第8卷(第1期): 66~68.[10](美)David Kincaid, Ward Cheney .王国荣,俞耀明,徐兆亮译.数值分析(原书第三版)[M].北京:机械工业出版社,2005,(9): 400~403.[11]M.T.Heath. Scientific Computing:An Introductory Survey, Sscond Edition[M].清华大学出版社.英文影印版. 2001,(10): 351~355.[12]封建湖,车刚明,聂玉.数值分析原理[M].北京:科学出版社,2001,(9): 111~114.[13]杨大地,涂光裕.数值分析[M].重庆:重庆大学出版社,,2006,(9): 139~142.[14] 黄明游,刘播,徐涛.数值计算方法[M].北京:科学出版社,2005,(8):137~138.[15]Jeffery J.Leader. Numerical Analysis and Scientific Computation[M].英文影印本.北京:清华大学出版社,2005,(8): 342~349。
(4) 北京理工大学函大2004-2005学年第1学期计算机科学与技术专业专升本数值计算方法思考题和习题教科书:《科学与工程计算》廖晓钟赖汝编国防工业出版社 2003年版第1 章思考题p26 1,2,3,4,5第1 章习题pp26-27 1,3,4,5,6,11第2 章思考题p66 1,3,6,7,8,9,12.13第2 章习题pp67-68 2,3,4,5,7,11,12,13,14,17,18第3 章思考题p119 1,3,4,5,6,10,18,19第3 章习题pp119-121 1,2,3,4,5,12,13第4 章思考题p144 1,2,3,4,5,7,8第4 章习题pp144-146 1,2,3,4,5,6,7,10,11,12,13第5 章思考题p207 1,2,3,4,5,6,7,9,10,11,12.13第5 章习题pp208-209 1,2,3,4,5,6,7,8,9,10,11,12,13,15第6 章思考题p257 1,2,3,4,5,6,7,8,10,11,12.14第6 章习题pp257-259 1,2,3,4,5,6,7,8,11,12,13,15,16,17,18第7 章思考题p292 1,2,3,4,5,6,8,9第7 章习题pp293-295 1,2,3,4,5,6,7,8,11,12,20作业题第1 章习题pp26-27 1(1),(2),3(3),5,6第2 章习题pp67-68 2,4,5,11,13,17第3 章习题pp119-121 1(1),2(1),5(2),12第4 章习题pp144-146 1(1),2,10,11,12,13第5 章习题pp208-209 1,3,4,7,10,13,,15第6 章习题pp257-259 1(2),3,6(1),12,16第7 章习题pp293-295 1,3,6,11,20数值计算方法复习题第1 章绪论1.说明数值算法的意义,计算机解题步骤和数值算法的特点。
数值分析(p11页)4 试证:对任给初值x 0,0)a >的牛顿迭代公式112(),0,1,2,......k ak k x x x k +=+= 恒成立下列关系式:2112(1)(,0,1,2,....(2)1,2,......kk k x k x x k x k +-=-=≥=证明:(1)(21122k k k k k kx a x x x x +-⎫⎛-=+==⎪ ⎝⎭(2) 取初值00>x ,显然有0>k x ,对任意0≥k ,a a x a x x a x x k k k k k ≥+⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎭⎫ ⎝⎛+=+2121216 证明:若k x 有n 位有效数字,则n k x -⨯≤-110218, 而()k k k k k x x x x x 288821821-=-⎪⎪⎭⎫⎝⎛+=-+ nnk k x x 2122110215.22104185.28--+⨯=⨯⨯<-∴>≥ 1k x +∴必有2n 位有效数字。
8 解:此题的相对误差限通常有两种解法. ①根据本章中所给出的定理:(设x 的近似数*x 可表示为m n a a a x 10......021*⨯±=,如果*x 具有l 位有效数字,则其相对误差限为()11**1021--⨯≤-l a x x x ,其中1a 为*x 中第一个非零数)则7.21=x ,有两位有效数字,相对误差限为025.010221111=⨯⨯≤--x x e 71.22=x ,有两位有效数字,相对误差限为025.010221122=⨯⨯≤--x x e 3 2.718x =,有两位有效数字,其相对误差限为:00025.010221333=⨯⨯≤--x e x ②第二种方法直接根据相对误差限的定义式求解 对于7.21=x ,0183.01<-e x∴其相对误差限为00678.07.20183.011≈<-x e x 同理对于71.22=x ,有003063.071.20083.022≈<-x e x 对于718.23=x ,有00012.0718.20003.033≈<-x e x备注:(1)两种方法均可得出相对误差限,但第一种是对于所有具有n 位有效数字的近似数都成立的正确结论,故他对误差限的估计偏大,但计算略简单些;而第二种方法给出较好的误差限估计,但计算稍复杂。
《计算方法》期中复习试题一、填空题:1、已知3.1)3(,2.1)2(,0.1)1(===f f f ,则用辛普生(辛卜生)公式计算求得⎰≈31_________)(dx x f ,用三点式求得≈')1(f 。
答案:,2、1)3(,2)2(,1)1(==-=f f f ,则过这三点的二次插值多项式中2x 的系数为 ,拉格朗日插值多项式为 。
答案:-1,)2)(1(21)3)(1(2)3)(2(21)(2--------=x x x x x x x L3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字;4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( );答案)(1)(1n n n n n x f x f x x x '---=+5、对1)(3++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 );6、计算方法主要研究( 截断 )误差和( 舍入 )误差;7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时,二分n 次后的误差限为( 12+-n a b );8、已知f (1)=2,f (2)=3,f (4)=,则二次Newton 插值多项式中x 2系数为( ); 11、 两点式高斯型求积公式⎰1d )(xx f ≈(⎰++-≈1)]3213()3213([21d )(f f x x f ),代数精度为( 5 );12、 为了使计算32)1(6)1(41310---+-+=x x x y 的乘除法次数尽量地少,应将该表达式改写为11,))64(3(10-=-++=x t t t t y ,为了减少舍入误差,应将表达式19992001-改写为199920012+ 。
13、 用二分法求方程01)(3=-+=x x x f 在区间[0,1]内的根,进行一步后根的所在区间为 ,1 ,进行两步后根的所在区间为 , 。
C语言编程习题第二章习题2-25.用二分法编程求6x4 -40x2+9=0 的所有实根。
#include <stdio.h>#include <math.h>#define N 10000double A,B,C;double f(double x){return (A*x*x*x*x+B*x*x+C);}void BM(double a,double b,double eps1,double eps2){int k;double x,xe;double valuea = f(a);double valueb = f(b);if (valuea > 0 && valueb > 0 || valuea <0 && valueb < 0) return;printf("Finding root in the range: [%.3lf, %.3lf]\n", a, b);for(k=1;k<=N;k++) {x=(a+b)/2;xe=(b-a)/2;if(fabs(xe)<eps2 || fabs(f(x))<eps1) {printf("The x value is:%g\n",x);printf("f(x)=%g\n\n",f(x));return;}if(f(a)*f(x)<0) b=x;else a=x;}printf("No convergence!\n");}int main(){double a,b,eps1,eps2,step,start;printf("Please input A,B,C:\n");scanf("%lf %lf %lf",&A,&B,&C);printf("Please input a,b, step, eps1,eps2:\n");scanf("%lf %lf %lf %lf %lf",&a,&b,&step,&eps1,&eps2);for (start=a; (start+step) <= b; start += step) { double left = start;double right = start + step;BM(left, right, eps1, eps2);}return 0;}运行:Please input A,B,C:6 -40 9Please input a,b, step, eps1,eps2:-10 10 1 1e-5 1e-5Finding root in the range: [-3.000, -2.000]The x value is:-2.53643f(x)=-0.00124902Finding root in the range: [-1.000, 0.000]The x value is:-0.482857f(x)=0.00012967Finding root in the range: [0.000, 1.000]The x value is:0.482857f(x)=0.00012967Finding root in the range: [2.000, 3.000]The x value is:2.53643f(x)=-0.00124902有时若把判别语句if(fabs(xe)<eps2 || fabs(f(x))<eps1)改为if(fabs(xe)<eps2 && fabs(f(x))<eps1)会提高精度,对同一题运行结果:Finding root in the range: [-3.000, -2.000]The x value is:-2.53644f(x)=-4.26496e-007Finding root in the range: [-1.000, 0.000]The x value is:-0.482861f(x)=-7.3797e-006Finding root in the range: [0.000, 1.000]The x value is:0.482861f(x)=-7.3797e-006Finding root in the range: [2.000, 3.000]The x value is:2.53644f(x)=-4.26496e-007习题2-35. 请用埃特金方法编程求出x=tgx在4.5(弧度)附近的根。
《计算方法》期中复习试题一、填空题:1、已知3.1)3(,2.1)2(,0.1)1(===f f f ,则用辛普生(辛卜生)公式计算求得⎰≈31_________)(dx x f ,用三点式求得≈')1(f 。
答案:2.367,0.252、1)3(,2)2(,1)1(==-=f f f ,则过这三点的二次插值多项式中2x 的系数为 ,拉格朗日插值多项式为 。
答案:-1,)2)(1(21)3)(1(2)3)(2(21)(2--------=x x x x x x x L3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字;4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( );答案)(1)(1n n n n n x f x f x x x '---=+5、对1)(3++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 );6、计算方法主要研究( 截断 )误差和( 舍入 )误差;7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时,二分n 次后的误差限为( 12+-n a b );8、已知f (1)=2,f (2)=3,f (4)=5.9,则二次Newton 插值多项式中x 2系数为( 0.15 ); 11、 两点式高斯型求积公式⎰1d )(xx f ≈(⎰++-≈1)]3213()3213([21d )(f f x x f ),代数精度为( 5 );12、 为了使计算32)1(6)1(41310---+-+=x x x y 的乘除法次数尽量地少,应将该表达式改写为11,))64(3(10-=-++=x t t t t y ,为了减少舍入误差,应将表达式19992001-改写为199920012+ 。
13、 用二分法求方程01)(3=-+=x x x f 在区间[0,1]内的根,进行一步后根的所在区间为 0.5,1 ,进行两步后根的所在区间为 0.5,0.75 。
数值计算方法第章数值积分第2章-数值积分电子科技大学自动化工程学院彭晓明1要主要内容机械求积z1.机械求积z2.牛顿-柯特斯公式z3.龙贝格算法z4.高斯求积2机械求积z1.机械求积z2.牛顿-柯特斯公式z3.龙贝格算法z4.高斯求积3z 对于定积分问题()bI f x dx =¾如果很容易求得不定积分a ∫则I=F(b)F(a)()F f x dx=∫I=F(b)-F(a)¾但是,现实中这往往很困难,举例sin xi 4x 2sin xz 根据积分中值定理bx ξ=−即:存在一个底为(b a)()()()a f dx b a f ∫即:存在个底为(b-a)而高为f(ξ)的矩形,其面积等于所求积分。
z 现在的问题是,ξ是未知的!z 数值积分的基本思想是找到一种近的方法似f(ξ)的方法。
5z 举例¾梯形公式()[]()()b b a f x dx f a f b −=+∫2a ⎛¾中矩形公式()()2b aa b f x dx b a f +⎞=−⎜⎟∫Si ⎝⎠¾Simpson 公式bb a a b d −+⎡⎤⎛6()()4()62a f x dx f a f f b ⎞=++⎜⎟⎢⎥⎝⎠⎣⎦∫一般形式z 般形式取[a, b]内若干个节点x 处的函数值[,]k f(x k ),通过加权平均的方法生成)这类求积公式称f(ξ) ,这类求积公式称机械求积公式,可以表示为n b()()0k k a k f x dx A f x =≈∑∫求积节点7求积系数z特点¾求积系数A仅与积分节点x的选取k k有关,而不依赖于被积函数f(x)的具体形式¾利用积分节点上的函数值计算积分值,从而把积分问题转化为函数值的计算8z数值求积方法是近似方法,自然希望对于“尽可能多”的函数是准确的。
z如果求积公式对于次数≤m的多项k(k=0,1,2,…,m)式x(k0, 1, 2,…,m)均能准确成立,但对x m+1不准确,则称机械求积公式具有m次代数精度。
式有z问题:梯形公式具有几次代数精度问题梯形公式具有几次代数精度?9n ⎧z 如果()()0bk k a k f x dx A f x ==⎪⎪∑∫()()n b k k a g x dx A g x ⎨⎪=∑∫0k =⎪⎩n ⎡¾f()()()()0()()b k k k a k f x g x dx A f x g x αβαβ=+=+⎤⎡⎤⎣⎦⎣⎦∑∫如果求积公式对于f(x)和g(x)准确成立,则求积公式对于f(x)和g(x)的线性组合也准确成立。
10代数精度z 以代数精度作为标准构造求积公式¾令求积公式对f(x)=x k (k=0, 1,均能准确成立构造求积2,…,m)均能准确成立,即构造求积系数满足系数,满足求积公式具有m 次代数精度⎧0122n A A A b a ""+++=−⎪=−0011()/2n n A x A x A x b a ""+++⎪⎨⎪11110011()/(1)m m m m m n nA x A x A x b a m "++⎪+++=−+⎩z (k=0,1,设已给出f(x)在插值节点x k (k 0, 1, 2,…,n)处的函数值,作插值多项式()()()nn k k p x f x l x =∑0k =nj x x l x −=¾()0k j k jj kx x =≠−∏p n (x)是次数≤n 的多项式)=f(x ¾p n (x k )f(x k )12(x)z 利用p n (x) 取代f(x)做积分bbx ≈¾由于()()n aaf dx p x dx∫∫()()()()()nnbb bn kkkkp x dx f x l x dx f x l x dx==∑∑∫∫∫对比aaa k k ==()nkkA f x ∑0k =bd可知这时13()k k aA l x dx =∫z 所谓“插值型的求积公式”有两重所谓插值型的求积公式有两重含义¾插值函数p n (x) 取代f(x)做积分bbx dx p x dx≈()()n aaf ∫∫¾积分系数A k 可以表示为bd14()k k a A l x dx =∫z 由于¾任意次数≤n 的多项式f(x)的插值多()项式p n (x)就是其本身,所以b bd d 至少具有n 次代数()()n a a f x dx p x dx ≈∫∫精度¾如果至少具有n 次()()bbnaaf x dx p x dx ≈∫∫代数精度,则n15()()b k j k j kaj l x dx A l x A ===∑∫•以下回顾拉格朗日插值多项式16z 一般情形般情形¾l (x)的构造k ()则()0,k ≠⎧如果满足,则p n (x)()1,k j j l x j k=⎨=⎩满足条件,意味着()()nk j l x c x x =−∏17标量0j j k=≠z 一般情形般情形¾l (x)的构造k ()根据()1n和l k (x k )=1,()0()k j j l x c x x ==−∏我们有n−j k≠()0j j k x x l x x x ==−∏18jj k ≠z 所以求积公式具有n 次代数精度的充要条件是它是插值型的。
z 求解A k 的方法⎧0122n A A A b a ""+++=−⎪=−方法10011()/2n n A x A x A x b a ""+++⎪⎨⎪110011()/(1)m m m m m n nA x A x A x b a m "++⎪+++=−+⎩19()bk k a A l x dx=∫方法2z 方法2举例构造下列形式的插值型求积公式1113x dx A f A f A f ⎛⎞⎛⎞⎛⎞⎟⎟⎟⎜⎜⎜≈++⎟⎟⎟⎜⎜⎜0120()424f ⎟⎟⎟⎜⎜⎜⎝⎠⎝⎠⎝⎠∫2bd ()003a A l x dx ==∫()1113b a A l x dx ==−∫20()2223b a A l x dx ==∫机械求积z1.机械求积z2.牛顿-柯特斯公式z3.龙贝格算法z4.高斯求积21形式z [a,b]h=(b-将区间[a, b]分为n 等份,步长h (b a)/n ,积分节点选取为等分点构出x k =a+kh(k=0, 1, 2,…,n),构造出下面的插值型求积公式nx 称为(N t ()()0n k k k I b a C f ==−∑柯特斯系数n 阶“牛顿-柯特斯(Newton-”公式。
nCotes)公式。
22()nk kk I A f x ==∑z 令x=a+th, 则1n −0b jk aj kjx xC b a xx ==−−∏∫j k≠=−=∏¾”The Calculus of ()(),0,1,2,...,!!t j dt k nn k n k ⋅−∫具体推导过程见” The Calculus of Observations A Treatise on152156p.152-156z -1阶牛顿柯特斯公式对应梯形公式(代数精度为1)z 2阶牛顿-柯特斯公式对应Simpson 公式(代数精度为3)z 4阶牛顿-柯特斯公式称为柯特斯公式(代数精度为5)73212()327b a C x f x f x f x f x −=24[]01234()()()()90f ++++积分中值定理(Weighted Meanz(Weighted Mean Value Theorem for Integrals) [Bradie, Appendix A]25z 梯形公式的余项¾根据牛顿插值公式,我们有1()()[,,]()()f x p x f a b x x a x b =+−−1()()[,,]()()b b bf x dx p x dx f a b x x a x b dx=+−−由在间aaa∫∫∫梯形公式由于在区间[a, b]上(x-a)(x-b)≤0,应用积分中值定理我们有应用积分中值定理,我们有26•以下回忆牛顿插值公式27牛顿插值公式顿z 牛顿插值公式[][]0010(),()f x f x f x x x x =+−+[]01201,,()()...f x x x x x x x −−++−−−[][]010110101,,...,()()...(),,...,,()()...()n n n n f x x x x x x x x x f x x x x x x x x x x −+−−−()()n p x R x =+28z 梯形公式的余项(,,)()()(,,)()()b b f a b x x a x b dx f a b x a x b dxξ−−=−−)6''(=−[](,,...,,f f x x x x ξ=()12f ξz 一般地,如果积分公式的余项对应阶导数则该积分被积函数的第m 阶导数,则该积分公式的代数精度为m-129z Simpson 公式的余项可以表示为(详见Bradie, p.462)()5(4)ˆˆ)b a −¾因此Si (),2880fa bξξ−<<因此,Simpson 公式的代数精度为3z 类似地分析可知,柯特斯公式代数精度为530z 不能通过提高阶数来提高()ba f x dx ∫的求解精度。
种有效的方法是(it )z 一种有效的方法是复化(composite)-[a,b]牛顿柯特斯求积:将[a, b]分为n 等份,步长h=(b-a)/n ,积分节点为x k =a+kh (0≤k ≤n),先用低阶的求x 积公式求得每个子段[x k , x k+1]上的积分值I ,然后用作为所求1n I −k 积分的近似值。
31k k =∑z 复化梯形公式¾回顾一下,采用梯形公式近似积分时我们有b []3()()()()''()212b ab a a f x dx f a f b f ξ−−=+−∫¾按照前面的“复化牛顿-柯特斯求积”思路,我们有−3211()()k kn b x ax k f x dx f x dx+==∑∫∫z 复化梯形公式¾而11k n b x+−余项)()3k ξ12''=−()()2()()k k f a f x f b f ξ⎢⎥++⎢⎥∑∑其中x k <ξk <x k+133复化梯形公式z 复化梯形公式¾对上面余项进一步分析,利用极值(B di A di A)定理(Bradie, Appendix A)可知,区间[a, b]内存在两个数c 1和对于所有的c 2,对于所有的k 有下式成立''''''34()()()21k f c f f c ξ≤≤z 复化梯形公式−¾也即(B di A di A)()()()12101''''''n k k f c f f c n ξ=≤≤∑根据介值定理(Bradie, Appendix A) 可知,区间[a, b]内存在数ξ,使得()11''()''n k f f ξξ−=∑35k n =z 复化梯形公式¾于是,余项可以表示为(((3321()''''''n h nh b a h f f f ξξξ−−−=−=−项我们称复化)))0121212k k =∑¾注意上式中包含h 2项,我们称复化梯形公式的收敛速度(rate of (convergence)为O(h 2)36z复化Simpson公式¾由于在运用Simpson公式时已经把区间[a, b]分为了两段,所以在应用[b]分为了两段所以在应用复化Simpson公式时需要把区间[a,[a, b]分为偶数段¾令n=2m,步长h=(b-a)/n=(b-a)/(2m),xi=a+ih (0≤i 2m),则复)/(2)i+ih(0i≤2)则复(Bradie,化Simpson公式可以写为(Bradie, p.470)37z 复化Simpson 公式b 复化Simpson 公式42()()a f x f x f b ⎢⎥2124()()3j j f −+++⎢⎥∑∑()b a h −¾收敛速度为O(h 4)()38z 举例计算,分别采用复化梯形公式对区间0sin xdx π∫公式和复化Simpson 公式,对区间[0, π都分成8段。