数值微分与数值积分.ppt

关于数值微分与数值积分现在学习的是第1页，共47页（1 1）向前差商）向前差商hxfhxfxf)()()(000由泰勒（由泰勒（TaylorTaylor）展式：）展式：hxxfhxf hxfhxf 002000)(!2)()()(向前差商的截断误差为：向前差商的截断误差为：)()(2)()()()(0000hOfhhxfhxfxfxR 2现在学习的是第2页，共47页hhxfhxfxf2)()()(000（2 2）向后差商）向后差商hhxfxfxf)()()(000向后差商的截断误差阶为：向后差商的截断误差阶为：)()()()()(0000hOhhxfxfxfxR（3 3）中心差商）中心差商3由泰勒（由泰勒（TaylorTaylor）展式：）展式：现在学习的是第3页，共47页中心差商的截断误差阶为：中心差商的截断误差阶为：hxhxhOfhffhhhxfhxfxfxR 00222120000)()(6)()(122)()()()()(!3)(!2)()()(1302000fhxfhxf hxfhxf )(!3)(!2)()()(2302000fhxfhxf hxfhxf 4现在学习的是第4页，共47页2.1.2 2.1.2 插值型数值微分插值型数值微分)(,(iixfx以以为插值点的插值多项式为为插值点的插值多项式为)(xLnniiinxfxlxLxf0)()()()(niiinxfxlxLxf0)()()()(jxx 当当时时njxfxlxLxfniijijnj,.,1,0)()()()(0误差项：误差项：niinxxnfdxdxR0)1()()!1()()(5现在学习的是第5页，共47页)!1()()()()1(0nfxxxRnnjiiijj例例2.2 2.2 给定三点给定三点)(,(),(,(),(,(221100 xfxxfxxfx并且并且hxxxx0112求求)(),(),(210 xfxfxf解解 过三点的插值多项式为过三点的插值多项式为)(2)()()()(2)()(2210122002212xfhxxxxxfhxxxxxfhxxxxxL)2(2)()2()()2(2)()()(1022202121202xxxhxfxxxhxfxxxhxfxLxf6现在学习的是第6页，共47页)()(4)(3(21)(2100 xfxfxfhxf)()(21)(201xfxfhxf)(3)(4)(21)(2102xfxfxfhxf得三点公式：得三点公式：7现在学习的是第7页，共47页对定义在区间对定义在区间a,b上的定积分上的定积分)()()(aFbFdxxffIba以上公式多称为牛顿以上公式多称为牛顿-莱布尼兹公式，莱布尼兹公式，F(x)为为f(x)的原函数的原函数.但有时原但有时原函数不能用初等函数表示，有时原函数又十分复杂，难于求出或计算函数不能用初等函数表示，有时原函数又十分复杂，难于求出或计算.如被如被积函数为积函数为:xxxsin,e2等函数的积分都无法解决，当被积函数为一组数据时，更是无能等函数的积分都无法解决，当被积函数为一组数据时，更是无能为力为力.为解决定积分的近似计算为解决定积分的近似计算,从定积分的定义：从定积分的定义：2.2.1 2.2.1 求积公式求积公式8 2.2 2.2 数值积分数值积分现在学习的是第8页，共47页9有nkkknbaxfdxxf1)(lim)(这样就避开了求原函数的运算这样就避开了求原函数的运算.上式就叫做求积公式，上式就叫做求积公式，Ak(k=0,1,n)与函数与函数f(x)无关，叫做求积系数无关，叫做求积系数,显然要确定一个求显然要确定一个求积公式，要确定求积结点积公式，要确定求积结点xk和求积系数和求积系数Ak，或者说不同的求积结点，或者说不同的求积结点和求积系数将确定不同的求积公式和求积系数将确定不同的求积公式.2.2.2 求积公式的余项和代数精度求积公式的余项和代数精度nkkknbaxfAfIdxxffI0)()()()(一般情况下，上式两端并不相等一般情况下，上式两端并不相等.我们称我们称:nkkkbaxfAdxxffR0)()(为求积公式为求积公式 的的余项余项，或，或截断误差截断误差.现在学习的是第9页，共47页10 为考查一个求积公式的误差，通常用为考查一个求积公式的误差，通常用代数精度代数精度来表示，如果一个来表示，如果一个求积公式对于不超过求积公式对于不超过m次的多项式都能够精确成立次的多项式都能够精确成立(Rf0)，而对，而对m+1次以上的多项式不能精确成立，则称该求积公式的代数精度为次以上的多项式不能精确成立，则称该求积公式的代数精度为m.例如求积公式：例如求积公式：验证当验证当f(x)=xm，m=0,1,2,3,4时，是否有时，是否有Rxm=0)1()0()1()(31343111fRfffdxxf1112,1)(313431右边左边xf10)1(0,)(3134312)1(122右边左边xxf3313343314)1(1310)1(0,)(44右边左边xxf32431434431525)1(1410)1(,)(55右边左边xxf231234231323)1(1210)1(,)(33右边左边xxf现在学习的是第10页，共47页11所以以上求积公式的代数精度为所以以上求积公式的代数精度为3.3.任何一个求积公式的代数精度至少为零任何一个求积公式的代数精度至少为零,即取即取f(x)=1)=1时公式应精确时公式应精确成立成立，这是求积系数应满足的起码条件，可以用它检验一个求积公式的系，这是求积系数应满足的起码条件，可以用它检验一个求积公式的系数的正确性数的正确性.求代数精度的另一种方法：利用余项公式求代数精度的另一种方法：利用余项公式 0 fRn求代数精度求代数精度m.m.如：如：)(90)4(52fhfR若取：若取：,132xxxf 则则0)()4(f0)()4(f而而,4xf 有有0)()4(f0)()4(f所以求积公式的代数精度为所以求积公式的代数精度为3.3.现在学习的是第11页，共47页122.2.3 插值型数值积分插值型数值积分)()()(xRxPxfnn 由插值可知，对任一函数由插值可知，对任一函数f(x)(包括表格形式的函数包括表格形式的函数)可用一可用一n次多项次多项式对其插值，即式对其插值，即babannbadxxRdxxPdxxf)()()(因此当当Pn(x)为拉格朗日插值多项式时，即为拉格朗日插值多项式时，即nkkknxfxlxP0)()()()()()()()!1()()()()(00)1(0fRxfAfRxfdxxldxxnfdxxfxldxxfnnkkknnkkbakbabannkkkba 则现在学习的是第12页，共47页13其中其中:)1.2()()()(0nkkkbanxfAfIdxxf通常将公式通常将公式(2.1)(2.1)叫做叫做插值型求积公式插值型求积公式.（2.22.2）称为求积公式的）称为求积公式的余项余项或误差。或误差。bankkkkkknkkbakkdxxxxxxxxxxxxxxxxxdxxlA)()()()()()()(110110)2.2()(!)1()()()1(dxxnffRbann现在学习的是第13页，共47页例：例：建立建立0,20,2上以上以x x0 0=0,x=0,x1 1=0.5,x=0.5,x2 2=2=2为节点的为节点的20)(dxxf数值积分公式。数值积分公式。解解 bankkkkkknkkbakkdxxxxxxxxxxxxxxxxxdxxlA)()()()()()()(11011095)5.02)(02()5.0)(0()(916)25.0)(05.0()2)(0()(31)2()5.00()2)(5.0()(2020222020112002000dxxxdxxlAdxxxdxxlAdxxxxdxxlA)2(5)5.0(16)0(3(91)()(220ffffIdxxf现在学习的是第14页，共47页 为便于上机计算，通常在求积公式中取等距节点，即将积分区间为便于上机计算，通常在求积公式中取等距节点，即将积分区间a,bn等分，即令等分，即令h=(b-a)/n，且记，且记x0=a,xn=b，则节点为，则节点为xk=x0+kh(k=0,1,n)，作变换，作变换:t=(x-x0)/h，代入求积系数公式：，代入求积系数公式：15nknnnknnbankkkkkknkkbakkdtntktktttkknhhdtkknhntktkttthdxxxxxxxxxxxxxxxxxdxxlA00110110)3.2()()1)(1()1(!)!()1(!)!()1()()1)(1()1()()()()()()()(这种等距节点的求积公式通常叫做这种等距节点的求积公式通常叫做牛顿牛顿-柯特斯公式柯特斯公式,下面介绍几个常下面介绍几个常用的公式用的公式:2.2.4 Newton-Cotes积分积分现在学习的是第15页，共47页取取a=x0,b=x1,(即即n=1)，代入，代入(2.3)式得式得.)(2)1(211)d-(!1)1()1(01021010CababhttthA1 梯形公式梯形公式所以梯形公式为所以梯形公式为)4.2()()(2)()(bfafb-afTdxxfba16.)(221d!0)1()1(11021001CababhttthA从从图图2.1看到，这是用梯形面积近似代替曲边梯形的面积，梯形公式的误看到，这是用梯形面积近似代替曲边梯形的面积，梯形公式的误差估计为：差估计为：现在学习的是第16页，共47页ya0 xb图图2.1梯形求积公式的余项：梯形求积公式的余项：)(12)()(!2)(R31a,bfabdxbxaxfba )()(!2)()()(1a,bbxaxfxLxf 梯形求积公式的代数精度：梯形求积公式的代数精度：取取 2,1)(xxxf0)(,)(;0)(,1)(xfxxfxfxf 012)(R31 fab所以，代数精度为所以，代数精度为1 1。现在学习的是第17页，共47页18例例 利用梯形公式计算利用梯形公式计算dxxI10214解解:.3)24(2111401420122I2 2 抛物形（辛卜生）公式抛物形（辛卜生）公式取取a=x0,(a+b)/2=x1,b=x2,(,(即即n=2)=2)，代入，代入(2.3)(2.3)式得式得.33222)-1)(-(!2)1(2020hhdttthA.34342)-(!1)1(201hhdttthA.3322!22hhdttthA201)-(现在学习的是第18页，共47页19所以抛物形（所以抛物形（SimpsonSimpson）公式为）公式为(2.5)(24)(3)(bfbafafhdxxfba其中其中h=(b-a)/2,上式也可写成上式也可写成:(2.6)(24)(6)(bfbafafabdxxfbaSimpsonSimpson公式的余项：公式的余项：)(2880)()(R)4(52a,bfabf现在学习的是第19页，共47页20或或)(180R)4(42a,bfhabf例例2 2 利用抛物形公式计算利用抛物形公式计算10214dxxI解解:1333.3)24(61)1()5.0(4)0(3564fffhISimpsonSimpson公式的代数精度：公式的代数精度：取取 432,1)(xxxxxf0)()(,0)()()4(4)4(3xfxxfxfxxf 02880)(R)4(52fab所以，代数精度为所以，代数精度为3 3。现在学习的是第20页，共47页(2.3)式给出式给出213 牛顿牛顿-柯特斯公式柯特斯公式)(00)()()1(!)!()1)()()1(!)!()1(nknknnknkCabdtktntttkknnabdtktntttkknhA其中其中:nknknkdtktntttkknnAabC0)()()1(!)!()1(1可以看出，可以看出，C(n)k不依赖函数不依赖函数f(x)和积分区间和积分区间a,b，可以事先计算出来，可以事先计算出来，通常叫做通常叫做牛顿牛顿-柯特斯系数柯特斯系数，下面给出，下面给出n从从16的牛顿的牛顿-柯特斯系数表柯特斯系数表：现在学习的是第21页，共47页22n11/21/221/64/61/631/83/83/81/847/9032/9012/9032/907/90519/28875/28850/28850/28875/28819/288641/840216/84027/840272/84027/840216/84041/840)(nkCCotesCotes系数表系数表现在学习的是第22页，共47页当当n为奇数，为奇数，f(x)C n+1a,b时，求积公式的余项为时，求积公式的余项为:bandxxxxxxxnff)()()!1(R101)(nn23 代数精度为代数精度为n.4 牛顿牛顿-柯特斯公式的余项与代数精度柯特斯公式的余项与代数精度当当n为偶数，为偶数，f(x)C n+2a,b时，求积公式的余项为时，求积公式的余项为:bandxxxxxxxxnff)()()!2(R10)2(nn代数精度为代数精度为n+1.现在学习的是第23页，共47页2.3 2.3 复化求积公式复化求积公式24 高阶的高阶的Newton-CotesNewton-Cotes公式不能保证数值积分的系列的收敛性；同时公式不能保证数值积分的系列的收敛性；同时，高阶，高阶Newton-CotesNewton-Cotes积分的计算是不稳定的。积分的计算是不稳定的。因而，在实际应用中往往采用因而，在实际应用中往往采用 将积分区间划分成若干个小区间，在各小区间上采用低次的求积公将积分区间划分成若干个小区间，在各小区间上采用低次的求积公式式(梯形公式或抛物形公式梯形公式或抛物形公式)，然后再利用积分的可加性，把各区间上，然后再利用积分的可加性，把各区间上的积分加起来，便得到新的求积公式的积分加起来，便得到新的求积公式.这就是复化求积公式的基本思想。这就是复化求积公式的基本思想。现在学习的是第24页，共47页2.3.1.复化梯形公式复化梯形公式用用n+1个分点将区间个分点将区间a，bn等分。每个区间长等分。每个区间长n,kkhax,nnabhk,210,1个分点 101)()(nkxxbakkdxxfdxxf则 在在xk，xk+1上用梯形公式，则上用梯形公式，则10111010k11)(2)()(2)()(2)(nkknkknkknkkkkbafRxfbfafhfRxfxfxxdxxf2511)(2)()(2nknkhafbfafhT记Tn叫做叫做复化梯形求积公式复化梯形求积公式，下标，下标n表示将积分区间等分的份数表示将积分区间等分的份数.现在学习的是第25页，共47页 公式的特点：节点公式的特点：节点xk(k=1,2,n-1)作为小区间的端点参与前、后两个作为小区间的端点参与前、后两个小区间的计算，因而系数为小区间的计算，因而系数为2，端点，端点a与与b只参与一次计算，系数为只参与一次计算，系数为1.26 如果在如果在Tn的基础上，将各小区间对分，这时节点数为的基础上，将各小区间对分，这时节点数为2n+1，分，分段数为段数为2n.记新的分点的函数值的和为记新的分点的函数值的和为n，则，则T2n应为原内节点与新增节应为原内节点与新增节点函数值的和的两倍，加上两端点点函数值的和的两倍，加上两端点a,b的函数值之和再乘上新区间长度的函数值之和再乘上新区间长度的一半，即的一半，即)7.2()(21)(2)(2)()(221)(2)(2)()(21111221111222nnniniiininiiinhTxfxfbfafhxfxfbfafhT)8.2(212nnnhTT或写为为对分新区间的长度.,2hh 这里现在学习的是第26页，共47页定理定理 设设 f(x)C2a,b，复化梯形公式的截断误差，复化梯形公式的截断误差 12R2a,bfhabTn 复化梯形积分的余项复化梯形积分的余项其中：其中：h=(b-a)/n27从这一公式可以看出，将区间对分后，原复化梯形公式的值从这一公式可以看出，将区间对分后，原复化梯形公式的值Tn作为作为一个整体保留一个整体保留.只需计算出新分点的函数值，便可得出对分后的积只需计算出新分点的函数值，便可得出对分后的积分值，不需重复计算原节点的函数值，从而减少了计算量分值，不需重复计算原节点的函数值，从而减少了计算量.利用复化求积公式的余项，我们可以估计出在满足精度的要求下，利用复化求积公式的余项，我们可以估计出在满足精度的要求下，应将积分区间等分多少份，即应将积分区间等分多少份，即n取多少取多少.这种误差估计方法称为事前误这种误差估计方法称为事前误差估计差估计.现在学习的是第27页，共47页例例 利用复化梯形公式计算利用复化梯形公式计算使其误差限为使其误差限为1010-4-4，应将区间，应将区间0,10,1几等分几等分?dxeIx10解解:因为被积函数因为被积函数xexf)(28xkexf)()(10)()(xexfk.10121)(1201422 ehfhTRn6.47n取取n=48=48可满足要求可满足要求.现在学习的是第28页，共47页29将积分区间将积分区间 a,b22m等分，等分，n=2=2m，节点为，节点为xk k=a+kh(k=0,1,2,=0,1,2,2,2m)，h=(b-a)/2m.在每两个小区间在每两个小区间x2 2k k,x2 2k k+2+2(k=0,1,2,=0,1,2,m-1)-1)上用抛物形公式，则有上用抛物形公式，则有:2.3.2 2.3.2 复化抛物形（复化抛物形（SimpsonSimpson）公式）公式 10101211210221221010)(4)(2)()(3)()(4)(3)()(222mkkmkkmkkmkkkkkmkmkxxbafRxfxfbfafhfRxfxfxfhdxxfdxxfkk)8.2()(4)(2)()(310121122mkkmkkmxfxfbfafhSS2 2m m叫做叫做复化抛物形求积公式复化抛物形求积公式，下标，下标2 2m表示积分区间等分的份数，表示积分区间等分的份数，2 2m强强调为偶数份调为偶数份.现在学习的是第29页，共47页30定理定理 设函数设函数f(x)C4a,b，则，则Simpson的余项为的余项为公式的特点为节点公式的特点为节点x2k,(k=1,2,m-1)作为小区间作为小区间x2k,x2k+2的端点，的端点，参与前后两次的辛普生公式的计算，因而系数为参与前后两次的辛普生公式的计算，因而系数为2，而奇数节点，而奇数节点x2k+1,(k=0,1,m-1)因辛普生公式中间点的求积系数为因辛普生公式中间点的求积系数为4而保留而保留4，前，前面的面的h/3为辛普生公式的公共求积系数为辛普生公式的公共求积系数.fhabS)4(42m180R复化抛物形（复化抛物形（SimpsonSimpson）公式的余项和代数精度）公式的余项和代数精度Simpson的代数精度为的代数精度为3现在学习的是第30页，共47页31解解:利用例利用例3的结果的结果.)()(exfk.101801)(1800144)4(42ehfhSRm.5.38526.210012,108526.21hmnh因此只需将区间因此只需将区间0,1四等分，即取四等分，即取m=2(n=4).例例4 利用复化抛物形公式计算利用复化抛物形公式计算使其误差限为使其误差限为10-4，应将区间，应将区间0,1几等分几等分?dxeIx10现在学习的是第31页，共47页2.3.3 复化积分的自动控制误差算法复化积分的自动控制误差算法在计算中常用在计算中常用 来估计误差。来估计误差。)()(2fTfTnn（1）的计算公式的计算公式)(2fTn对于对于 badxxf)(取分点取分点n=1,h1=(b-a),得：得：)()(21)()(2)(11bfafhbfafabfT对分对分(取分点取分点n=2),h2=(b-a)/2 得：得：22221222222),()(21)(2)()(221)()(2)()(2)(haxxfhfTxfbfafabbfxfhxfafhfT32现在学习的是第32页，共47页再对分再对分(取分点取分点n=4),h4=(b-a)/4 得：得：4341314231243,)()()(21)()(4)(21)(haxhaxxfxfhfTxfxfabfTfT33mimmmmhiafhfTfT1222)12()(21)(，再对分，再对分(取分点取分点n=2m),h2m=(b-a)/(2m),得：得：)()(21)(2fHfTfTmmm或：或：其中其中 mabhxfhfHmmiimm)(,)()(1021现在学习的是第33页，共47页34（2）计算公式的误差控制计算公式的误差控制)(2fTn|)()(|)()(|22fTfIfTfTmmm与比较比较 )(212)()()(12)()(222fhabfTfIfhabfTfImm )()(31)()(22fTfTfTfImmm若要若要|)()(|2fTfIm只需只需 3|)()(|2fTfTmm后者比前者便于控制后者比前者便于控制 现在学习的是第34页，共47页复化梯形积分的算法复化梯形积分的算法 1.输入误差控制精度输入误差控制精度tol，初始分点数，初始分点数n，端点，端点a,b,h=(b-a)/n 2.建立函数文件建立函数文件f(x)3.T2=(f(a)+f(b)*h/2;T1=T2+100;6.S=Sum(S1);T2=T1/2+h*S;Clear S1 end 5.For i=1:n/2 计算计算 S1(i)=f(a+(2*i-1)*h)end 4.While|T2-T1|tol T1=T2 h=h/2;n=2n;%区间分半区间分半357.输出输出T2 现在学习的是第35页，共47页362.3.4 龙贝格龙贝格(Romberg)积分公式积分公式我们已知的我们已知的T2n与与Tn的关系的关系1.Romberg积分公式积分公式hTTnn212于是可以逐次对分形成一个序列于是可以逐次对分形成一个序列T1,T2,T4,T8,此序列收敛于积分此序列收敛于积分真值真值I.当当|T2n-Tn|时时,取取T2n为为I的近似值的近似值.以上算法称为以上算法称为复化梯形公式复化梯形公式的逐次分半公式的逐次分半公式.但由于此序列收敛太慢但由于此序列收敛太慢,因此并不实用因此并不实用.现我们试图现我们试图将它改造成为收敛快的序列将它改造成为收敛快的序列.fhabTIn 212 fhabfhabTIn 412212222现在学习的是第36页，共47页10212)(22121nkknnnxfhThTT)()(2)(211bfxfafhTnkknhkaxnabhk)21(21为二分前的步长其中，现在学习的是第37页，共47页38如认为如认为 ff 则有则有42nnTITI于是有于是有:1442nnTTI记记)3134(14422nnnnnTTITTS这样我们从收敛较慢的这样我们从收敛较慢的Tn序列推出了收敛较快的序列推出了收敛较快的Sn序列序列.截断误差由截断误差由O(h2)提高到提高到O(h4).如认为如认为 )4()4(ff则有则有162nnSISI于是有于是有:116162nnSSI记记)1511516(1161622nnnnnSSISSC这样我们从这样我们从Sn序列又推出了收敛更快的序列又推出了收敛更快的Cn柯特斯序列柯特斯序列.截断误差提高到截断误差提高到O(h6),代数精度为代数精度为5.现在学习的是第38页，共47页39如认为如认为 )6()6(ff则有则有642nnCICI于是有于是有:164642nnCCI记记)6316364(1646422nnnnnCCCCR这样我们从这样我们从Cn序列又推出了收敛更快的序列又推出了收敛更快的Rn序列序列.Rn序列也称为龙贝序列也称为龙贝格序列格序列.截断误差提高到截断误差提高到O(h8),代数精度为代数精度为7.,.3,2,14411,11,1,kTTTjjkjkjjk Romberg积分公式：积分公式：对每一个对每一个k，j从从2做到做到k，一直做到，一直做到Tk,k-Tk,k-1tol时停止计算。时停止计算。现在学习的是第39页，共47页40T1T2S1T4S2C1T8S4C2R1T16S8C4R2表表7-27-2该表四个序列都是收敛的该表四个序列都是收敛的.现在学习的是第40页，共47页41例例5 利用龙贝格方法计算利用龙贝格方法计算解解:计算结果列如下表计算结果列如下表:dxxI10214i2iT序列S序列C序列R序列013.00000123.100003.13333243.131183.141573.14212383.138993.141593.141593.141594163.140943.141593.141593.14159这一结果与这一结果与I=相比较已有较好的精度相比较已有较好的精度.现在学习的是第41页，共47页算法公式：算法公式：abhbfafhT111,1),()(2)5.0(21111,11,2hafhTT.)5.1()5.0(212221,21,3hafhafhTT)21(21221111,11,ikiiiihkafhTTi=2,3,其中：其中：11221iiiabhh现在学习的是第42页，共47页jmjTTTjjmjmjjm,.,3,2,14411,11,1,算法：算法：dxxfIba)(求：)()(2,:11,1bfafhTabhstep53,.,3,2:2stepmistep做对)5.0(21:32211,11,ikiihkafhTTstep现在学习的是第43页，共47页144,.,3,2:411,11,1,jjijijjiTTTijstep对停机。若,|:51,jijiTTstep2:6hhstep否则停机。，输出T:7step例例5 利用龙贝格方法计算利用龙贝格方法计算dxxI5.10214程序：先建立函数文件，文件名程序：先建立函数文件，文件名,ff.mfunction y=ff(x)y=4/(1+x2)现在学习的是第44页，共47页主程序：主程序：Romberg1.ma=0,b=1.5,m=5;h=b-a;T=zeros(m);T(1,1)=h*(ff(a)+ff(b)/2;for i=2:m for k=1:2(i-2)ss(k)=ff(a+(k-0.5)*h);end T(i,1)=(T(i-1,1)+h*sum(ss)/2;for j=2:i T(i,j)=(4(j-1)*T(i,j-1)-T(i-1,j-1)/(4(j-1)-1);end if abs(T(i,j)-T(i,j-1)tol break end h=h/2;endT现在学习的是第45页，共47页462.4 重积分计算重积分计算现在学习的是第46页，共47页感谢大家观看现在学习的是第47页，共47页