数值积分课件(《计算方法》)剖析

数值计算方法第7章数值积分l1 插值型求积公式插值型求积公式 l2 复化求积公式复化求积公式l3 龙贝格龙贝格(Romberg)求积方法求积方法6/23/20241数值计算方法1 插值型求积公式插值型求积公式 l在一元函数的积分学中,我们已经熟知,若函数f(x)在区间a,b上连续且其原函数为F(x),则可用牛顿莱布尼兹公式(71)来求定积分。6/23/20242数值计算方法公式(71)虽然在理论上或在解决实际问题中都起了很大的作用,但它并不能完全解决定积分的计算问题。因为定积分的计算常常会碰到以下三种情况:l(1)被积函数f(x)的原函数F(x)不易找到。许多很简单的函数,例如其原函数都不能用初等函数表示成有限形式。6/23/20243数值计算方法l(2)被积函数f(x)没有具体的解析表达式。其函数关系由表格或图形表示,无法求出原函数。其被积函数的原函数就比较复杂,从数值计算角度来看,计算量太大。l(3)尽管f(x)的原函数能表示成有限形式但其表达式相当复杂。例如定积分6/23/20244数值计算方法图7.1l如图7.1,若用左矩形近似地代替曲边梯形,则得到左矩形公式(72)6/23/20245数值计算方法l同样可得到右矩形公式：(73)6/23/20246数值计算方法图7.2l如图7.2,若用梯形的面积近似地代替曲边梯形的面积,则得到计算定积分的梯形公式(74)6/23/20247数值计算方法l如图7.3,若用抛物线代替曲线f(x),则可得到抛物线公式(或辛普生公式)(75)图7.36/23/20248数值计算方法此外此外,众所周知的众所周知的梯形公式梯形公式:l I(f)(b-a)f(a)+f(b)/2和和 Simpson公式公式:l I(f)(b-a)f(a)+4f(a+b)/2)+f(b)/6则分分别可以看作用可以看作用 a,b,c=(a+b)/2,三点三点高度的加高度的加权平均平均值 f(a)+f(b)/2 和和 f(a)+4f(c)+f(b)/6作作为平均高度平均高度f()的近似的近似值.6/23/20249数值计算方法 更一般地更一般地,取区间取区间a,b内内n+1个点个点 xi,(i=0,1,2,n)处的高度处的高度f(xi)(i=0,1,n)通过通过加权平均加权平均的方的方法近似地得出平均高度法近似地得出平均高度f(),这类求求积方法称方法称为机械机械求求积:6/23/202410数值计算方法 或写成或写成:数值积分公式数值积分公式求积系数求积系数 求积节点求积节点(1)6/23/202411数值计算方法记记称称(2)为数值求积公式为数值求积公式,(3)为求积公式余项为求积公式余项(误误差差).构造或确定一个求积公式，要讨论解决的问构造或确定一个求积公式，要讨论解决的问题有题有(i)确定求积系数确定求积系数Ak和求积节点和求积节点xk；(ii)求积公式的误差估计和收敛性求积公式的误差估计和收敛性 为了构造形如式为了构造形如式(2)的求积公式的求积公式,需要提供需要提供一种一种判定求积方法精度高低准则判定求积方法精度高低准则6/23/202412数值计算方法求积公式的代数精度求积公式的代数精度定义定义1称求积公式称求积公式(2)具有具有m次次代数精度代数精度,如果它满如果它满足如下两个条件足如下两个条件:(i)对所有次数对所有次数 m次的多项式次的多项式 ,有有 (ii)存在存在m+1次多项式次多项式 ,使得使得定义定义1中的条件中的条件(i),(ii)等价于等价于:6/23/202413数值计算方法插值型求积公式插值型求积公式在积分区间在积分区间a,b 上上取取n+1个节点个节点xi，i=0,1,2,n,作作f(x)的的n次代数插值多项式次代数插值多项式（拉格朗（拉格朗日插值公式）日插值公式）:则有则有为插值余项为插值余项于是有于是有6/23/202414数值计算方法取取称称(4)式为插值型求积公式式为插值型求积公式,其中其中求积系数求积系数Ak由由(5)式确定式确定.(4)(5)Ak由由 节点节点 决定，决定，与与 f(x)无关。无关。6/23/202415数值计算方法6/23/202416数值计算方法推论推论1 求积系数满足求积系数满足:误误 差差定理定理1 形如形如 的求积公式至少有的求积公式至少有 n 次代数精次代数精度度 该该公式为公式为插值型插值型（即：（即：）6/23/202417数值计算方法l现用第六章介绍的插值多项式Pn(x)来代替被积函数f(x),即有 1.1 1.1 牛顿牛顿柯特斯公式柯特斯公式(NewtonCotesNewtonCotes)取节点为等距,即a=x0 x1xn=bl建立数值积分公式最基本的思想是选取一个既简单又有足够精度的函数(x),用(x)代替被积函数f(x),于是有6/23/202418数值计算方法l利用拉格朗日插值多项式(76)其中(77)6/23/202419数值计算方法l这里yi=f(xi),对式(76)两边积分得6/23/202420数值计算方法为牛顿牛顿柯特斯柯特斯(Newton-Cotes)(Newton-Cotes)求积公式求积公式,Rn(f)为牛顿柯特斯求积公式的余项。我们称6/23/202421数值计算方法令x=x0+sh,0sndx=hds=(b-a)/nds(711)6/23/202422数值计算方法Newton-Cotes公式的误差为公式的误差为:与与x有关有关注意注意:由由(7-11)式确定的式确定的Cotes系数只与系数只与i和和n有关有关,与与f(x)和积分区间和积分区间a,b无关无关,且且满足满足:(7-9)6/23/202423数值计算方法l称Ci(n)为柯特斯求积系数。很显然,当n=1时,可算得此时式(710)为(712)这是梯形公式梯形公式。6/23/202424数值计算方法l当n=2时,可得于是(713)这是抛物线抛物线(Simpson)公式公式。6/23/202425数值计算方法l当n=3时,代入(710)式得到求积公式6/23/202426数值计算方法类似地可分别求出n=4,5,时的柯特斯系数,从而建立相应的求积公式。具体结果见表71。l从表中可以看出,当n7时,柯特斯系数为正;从n8开始,柯特斯系数有正有负。因此,当n8时,误差有可能传播扩大,牛顿柯特斯求积公式不宜采用。l柯特斯系数Ci(n)仅与n和i有关,与被积函数f(x)无关,且满足(715)事实上,式(710)对f(x)=1是准确成立的。6/23/202427数值计算方法表716/23/202428数值计算方法l定理定理 当阶数当阶数n为偶数时为偶数时,Newton-Cotes公式公式(8)至少具有至少具有n+1次代数精度次代数精度.证明证明只需验证当只需验证当n为偶数时为偶数时,Newton-Cotes公式公式对对f(x)=xn+1的余项为零的余项为零.由于由于f(x)=xn+1,所以所以f(n+1)(x)=(n+1)!.由式由式(7-9)得得引进变换引进变换t=u+n/2,因为因为n为偶数为偶数,故故n/2为整数为整数,于是有于是有据此可断定据此可断定R(f)=0,因为上述被积函数是个奇函数因为上述被积函数是个奇函数.6/23/202429数值计算方法Newton-Cotes公式的数值稳定性公式的数值稳定性现在讨论现在讨论舍入误差舍入误差对计算结果产生的影响对计算结果产生的影响.设设用公式用公式 近似计算积分近似计算积分 时时,其中计算函数值其中计算函数值f(xj)有误差有误差j(j=0,1,2,n).设计算算Cj(n)没有没有误差差,中中间计算算过程中的舍入程中的舍入误差也差也不考不考虑,则在式在式(10)的的计算中算中,由由j引起的引起的误差差为(10)6/23/202430数值计算方法如果如果Cj(n)都是正数都是正数,并设并设故故en是有界的是有界的,即由即由j引起的误差受到控制引起的误差受到控制,不超不超过过的的(b-a)倍倍,保证了保证了数值计算的稳定性数值计算的稳定性.而当而当n7时时,Cj(n)将出现将出现负数负数,保证数值稳定性保证数值稳定性.因此高阶公式不宜采用因此高阶公式不宜采用,有实有实用价值的仅仅是几种用价值的仅仅是几种低阶的求积公式低阶的求积公式.将随将随n增大增大,因而因而不能不能则有则有6/23/202431数值计算方法解解利用梯形公式利用抛物线公式例1试分别用梯形公式和抛物线公式计算积分:原积分的准确值6/23/202432数值计算方法l现对牛顿柯特斯求积公式所产生的误差作一个分析。由式(79),牛顿柯特斯求积公式的余项为1.2 1.2 误差估计误差估计l易知,牛顿柯特斯求积公式(710)对任何不高于n次的多项式多项式是准确成立的。这是因为f(n+1)()0故Rn(f)06/23/202433数值计算方法l牛牛顿顿柯柯特特斯斯求求积积公公式式的的代代数数精精确确度度至至少少为为n n。通常在基点个数相等的情况下,代数精确度愈高,求积公式就愈精确。l一般说来,若某个求积公式对于次数不高于m的多项式都准确成立(即Rn(f)0),而对于某一次数为m+1的多项式并不准确成立,则称这一求积公式的代数精确度代数精确度为m。l定定理理1(梯形公式的误差)设f(x)在区间a,b上具有连续的二阶导数,则梯形求积公式的误差为6/23/202434数值计算方法由于1(x)=(x-a)(x-b)证由式(79)知,梯形公式的余项为l1(x)在区间(a,b)内不变号,f()是x的函数且在a,b上连续,故根据积分第二中值定理知,存在某一(a,b)使6/23/202435数值计算方法l定定理理2 2(抛物线公式的误差)设f(x)在a,b上有连续的四阶导数,则抛物线公式的误差为(717)6/23/202436数值计算方法n=1:Trapezoidal Rule/*令令 x=a+th,h=b a,用用中值定理中值定理*/代数精度代数精度=1n=2:Simpsons Rule代数精度代数精度=3n=4:Cotes Rule,代数精度代数精度=5,6/23/202437数值计算方法复合求积公式复合求积公式高次插值有高次插值有Runge 现象现象,高阶高阶Newton-Cotes公式会出现公式会出现数值不稳定数值不稳定,低阶低阶Newton-Cotes公式有时又不能满足精度公式有时又不能满足精度要求要求.解决这个矛盾的办法是将积分区间解决这个矛盾的办法是将积分区间 a,b 分成若干小区间分成若干小区间,在每个小区间上在每个小区间上用低阶求积公式计算用低阶求积公式计算,然后将它们加起来然后将它们加起来,这就是这就是复合求积方法复合求积方法.6/23/202438数值计算方法2复合求积公式l2.1复合梯形公式对于定积分(71),将积分区间a,b分成n个相等的子区间xi,xi+1,这里步长在每一个子区间xi,xi+1上使用梯形公式,则6/23/202439数值计算方法复化梯形公式积分法复化梯形公式积分法6/23/202440数值计算方法相加后得(718)(719)若f(x)在a,b上连续,由连续函数的介值定理,存在某一(a,b)使得6/23/202441数值计算方法因而于是得到复合梯形公式(721)其余项为6/23/202442数值计算方法例2若用复合梯形公式计算积分则当0 x1时,有因为又解由余项(721)式问积分区间要等分多少才能保证有五位有效数字?6/23/202443数值计算方法l由于原积分的准确值具有一位整数,因此要使近似积分值有五位有效数字,只需取n满足两边取对数得整理后得到取n=68.6/23/202444数值计算方法l类似复合梯形公式的做法,把区间a,b分成n个相等的子区间x2i,x2i+2(i=0,1,n-1),设每个子区间上的中点为x2i+1(i=0,1,n-1),且(722)2.2复合抛物线公式l在每一个子区间x2i,x2i+2上利用抛物线公式得6/23/202445数值计算方法复化复化Simpson公式积分法公式积分法6/23/202446数值计算方法l相加后得(723)6/23/202447数值计算方法l若f(4)(x)在a,b上连续,则从而得到复合抛物线公式(724)l其余项为(725)6/23/202448数值计算方法图7.4复合抛物线公式框图6/23/202449数值计算方法的数据表,分别用复合梯形公式和复合抛物线公式计算l例3已知函数 x f(x)011/80.99739782/80.98961583/80.97672674/80.95885105/80.93615566/80.90885167/80.877192510.84147096/23/202450数值计算方法l解用复合梯形公式,这里6/23/202451数值计算方法l用复合抛物线公式可得 比较上面两个结果比较上面两个结果T8和和S4,它们都需要提供它们都需要提供9个个点上的函数值工作量基本相同点上的函数值工作量基本相同,然而精度却差别很大然而精度却差别很大.同积分的同积分的准确值准确值I(f)=0.9460831比较比较,复化梯形法复化梯形法的结果的结果T8=0.9456909只有只有两位有效数字两位有效数字,而复化而复化Simpson法的结果法的结果S4=0.9460832却有却有六位有效数字六位有效数字.6/23/202452数值计算方法 复化梯形公式：复化梯形公式：在每个在每个 上用梯形公式：上用梯形公式：=Tn/*中值定理中值定理*/6/23/202453数值计算方法 复化复化 Simpson 公式：公式：6/23/202454数值计算方法2.3变步长公式l前面介绍的复合梯形公式和复合抛物线公式的步长都是预先确定的。它的主要缺点是事先很难估计出n的大小(或步长h的大小),使结果达到预先给定的精度。在实际计算中,我们常常借助于计算机来完成积分步长h的自动选择,即采用变步长求积公式。具体地讲,就是将步长逐次折半,反复利用复合求积公式,直到满足精度要求为止。6/23/202455数值计算方法l下面介绍变步长复合抛物线公式(变步长复合梯形公式留给读者作为练习)。(726)其中再把每个子区间分成两半,用逐次将区间a,b分成2,4,2m等分,并按复合抛物线公式逐次计算积分得到S1,S2,Sm,而6/23/202456数值计算方法作步长,按复合抛物线公式计算出积分的近似值S2m。对于相邻两次的积分近似值Sm、S2m,考察当S2m1当S2m1(727)设预先给定的精度为,若d则以S2m作为所要求的积分近似值,否则继续将区间分半,利用复合抛物线公式求积分,直到满足预给的精度为止。6/23/202457数值计算方法图7.5变步长复合抛物线公式6/23/202458数值计算方法图7.5变步长复合抛物线公式6/23/202459数值计算方法3龙贝格(Romberg)积分方法l我 们 已 经 知 道,当 被 积 函 数 f(x)在 区 间a,b上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分l可将分点(即基点)加密,也就是将区间a,b细分,然后利用复合梯形公式或复合抛物线公式求积。6/23/202460数值计算方法l若用Tm表示把a,b作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h=h/2,则T2m与Tm之间有如下关系:(528)其中6/23/202461数值计算方法l另外,若用Sm表示把a,b分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h就有从Tm的定义可得到关系式(529)6/23/202462数值计算方法l我们再举一个计算上半单位圆面积的例子(它的准确面积为/2)。现用内接正多边形的逼近方法来计算。l如图5.6,图(a)、图(b)是用同样的内接正多边形计算上半单位圆的面积。图(a)是用梯形方法计算其面积,图(b)是用三角形方法计算其面积。6/23/202463数值计算方法图5.66/23/202464数值计算方法l设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为同理6/23/202465数值计算方法l如果组合一下,就会得到更精确的结果,即同理6/23/202466数值计算方法l再以类似方法组合得这样继续下去,其值越来越接近上半单位圆面积/2。这种方法可以用到计算定积分6/23/202467数值计算方法l为了推广公式(529)和上述计算上半单位圆面积的组合方法,我们引进龙贝格求积算法。l龙贝格求积算法本来是利用所谓外推法构造出的一种计算积分的方法。为了避免从外推引入而带来理论上的麻烦,我们将直接从构造一个T数表开始。l首 先 将 a,b 依 次 作20,21,22,等分,记6/23/202468数值计算方法l按复合梯形公式(520)算得的值相应地记为T(k)0(k=0,1,2,);把按式(529)算得的S2m依次记为T(k)1(k=0,1,2,崐),而这每一个S2m又理解为由T2m与Tm的线性组合得到的改进值,即我们可按照类似的方法继续进行改进,也即由S2m与Sm的线性组合得到改进值,依次记为T(k)2(k=0,1,2,),即6/23/202469数值计算方法l这样就可构造出一个数表(5-30)6/23/202470数值计算方法l其中除第0列(即最左一列)的T(k)0是按复合梯形公式计算外,其余各列都按下述规则(对m)(531)递推地计算出来。箭头表示计算流程。其计算步骤为:(1)将区间a,b等分为20,用梯形公式计算T(0)0,即6/23/202471数值计算方法l(2)将区间a,b等分为21,用梯形公式算出T(1)0,即再由T(0)0,T(1)0根据公式(531)算出T(0)1,即若T(0)1-T(0)0,(为预给的精度)则停止计算;否则继续往下计算;6/23/202472数值计算方法l(3)依次分别算出T(2)0,T(1)1,T(0)2,这一行地往下推算,每一行算完,就得验证T(0)m(m=1,2,)是否满足预给的精度,即若则停止计算;否则继续进行下一行。为了便于在计算机上实现,可运用下列公式编制程序:6/23/202473数值计算方法6/23/202474数值计算方法图5.76/23/202475数值计算方法图5.76/23/202476数值计算方法l例4计算积分精确到10-4。解6/23/202477数值计算方法6/23/202478数值计算方法6/23/202479数值计算方法6/23/202480数值计算方法于是由于实际上6/23/202481