数值分析课件数值微积分.pdf

第五章第五章数值微积分数值微积分数值微积分数值微积分数值积分数值积分数值积分数值积分Newton-Leibniz 公式公式()()().bf x dxF bF a=Newton Leibniz 公式公式()()()af实际问题中的局限性：实际问题中的局限性：实际问题中的局限性：实际问题中的局限性：1）F(x)非常复杂非常复杂2）没有初等函数形式的）没有初等函数形式的F(x)数值积分数值积分3）函数由离散数据组成）函数由离散数据组成3）函数由离散数据组成）函数由离散数据组成5.1 数值求积公式的基本思想数值求积公式的基本思想()lim()nbf x dxf xx=由定积分的定义，由定积分的定义，max0()lim(),iiiaxif x dxf xx=axxxbxxx其中其中()()nbf x dxA f x011.,.niiiaxxxbxxx=其中其中，0()().iiaif x dxA f x=()()()nbfdA fR f或或0()()().iiaif x dxA f xR f=+或或机械求积法机械求积法：节点处函数值的线性组合作为积分近似值：节点处函数值的线性组合作为积分近似值两个问题两个问题两个问题两个问题：1）求积系数如何选取？1）求积系数如何选取？2）若求积节点可以自由选取，取哪些节点好？2）若求积节点可以自由选取，取哪些节点好？求积公式的求积公式的代数精度代数精度0()()nbiiaif x dxA f x=定义定义：若求积公式对一切：若求积公式对一切不高于不高于m次的多项式都准确成立，而对于至少一个次的多项式都准确成立，而对于至少一个m+1次多项式等号不成立，则称此公式的代数次多项式等号不成立，则称此公式的代数代数精度的求法代数精度的求法精度为精度为m.代数精度的求法代数精度的求法：依次取依次取(x)=1,x,x2,x3,，验证求积公式是，验证求积公式是依次取依次取(x)1,x,x,x,，验证求积公式是，验证求积公式是否成立。若第一个不成立的等式是否成立。若第一个不成立的等式是xm，则其代，则其代数精度是数精度是m-1.数精度是数精度是m 1.例例1：确定如下求积公式的待定系数：确定如下求积公式的待定系数A0A1A2，1()(1)(0)(1)f x dxA fA fA f+0、1、2，使其代数精度尽可能高。使其代数精度尽可能高。0121()(1)(0)(1).f x dxA fA fA f+解：取解：取 (x)=1，x，x2，有，有()01220AAAAA+=+=01/34/3AA=020202/3AAAA+=+=124/3.1/3AA=11()(1)4(0)(1)/3.fx dxfff=+则则取取 (x)=x3，左，左=右右=0，但，但 (x)=x4时，左时，左=-11x4dx=2/5 右右=2/3,()1,所以求积公式具有所以求积公式具有3次代数精度。次代数精度。例例2：选择常数：选择常数a，使如下求积公式代数精度尽量高：，使如下求积公式代数精度尽量高：0h(x)dx h(0)+(h)/2+ah2(0)-(h).例例2：选择常数：选择常数a，使如下求积公式代数精度尽量高：，使如下求积公式代数精度尽量高：解：取解：取 (x)=1，左，左=h=h（1+1）/2=右；右；取取 ()左左 h2/2 h（0 h）/2 右右取取 (x)=x，左左=h2/2=h（0+h）/2=右右；取；取 (x)=x2，左，左=h3/3，()令令 h3/2 2 h3h3/3得得1/12右右=h（0+h2）/2+ah2（-2h）=h3/2-2ah3，令令 h3/2-2ah3=h3/3，得得 a=1/12.取取 (x)=x3，左，左=h4/4，右，右=h(h3)/2+h2(-3h2)/12=h4/4；取取 ()，左，左，右，右()()；取取 (x)=x4，左，左=h5/5，右，右=h5/2+h2(-4h3)/12=h5/6，左左 右，所以右，所以 a=1/12 时，求积公式有时，求积公式有3次代数精度。次代数精度。5 2 Newton Cotes求积公式求积公式1 1插值型求积公式插值型求积公式5.2 Newton-Cotes求积公式求积公式1 1、插值型求积公式插值型求积公式给定一组节点给定一组节点a=x0 x1xn=b,构造n次Lagrange型插值多项式构造n次Lagrange型插值多项式01n0()()().nnkkkL xf x lx=0k其中为n次Lagrange插值基函数,其中为n次Lagrange插值基函数,0()()nkjjlxxx=0jj k=则则()()()f xL xR x=+则则()()().nf xL xR x=+()()()()bbnnaaI ff x dxIfL x dx=近似函数近似函数()()().nbiil x dx f x=()()nbiil x f x dx=iA()0()()iiaif=0()()iiaif=其中，只与求积节点有关，与被积函数无关.其中，只与求积节点有关，与被积函数无关.iA(1)()nbbf+误差误差()()()()().(1)!bbnnnaafI fIfR x dxx dxn=+可以看出，至少可以看出，至少 n 阶代数精度阶代数精度.2、等距节点情形（、等距节点情形（Newton-Cotes）h=ban,xi=a+ih,i=0,n,取步长取步长2、等距节点情形（、等距节点情形（Newton Cotes）nAi=li(x)bdx=(xx0)(xxi1)(xxi+1)(xxn)bdxii()a(xix0)(xixi1)(xixi+1)(xixn)ax=a+tht(t1)(ti+1)(ti1)(tn)hnd1(1)nin=t(t1)(ti+1)(ti1)(tn)i!(ni)!(1)nih0dt=nh1n(1)i!(ni)!(tk)k=0ki0ndt=(ba)ci(n).ki与区间选取无关与区间选取无关nb()0()()().bniiaif x dxbacf x=Newton-Cotes系数表系数表n11/21/2()nkc11/21/221/62/31/631/83/83/81/847/9016/452/1516/457/90()1.nnkc=注注：0kk=注注梯形公式梯形公式(n1)梯形公式梯形公式(n1)：11(1)(1)11(1)CtdtCtdt=0100(1),22CtdtCtdt=11()()()T fbff b()()().22T fb af af b=+yBy=L1(x)y=f(x)Ay1()f0f1x0 x0 x1Simpson公式公式(n2)：2(2)0011(1)(2),46Cttdt=2(2)1014(2),26Ct tdt=2(2)2011(1),46Ct tdt=46141 ()()()()b aS fbfff b+()()()().6626S fb af aff b=+y yy=L2(x)f()y=f(x)xx0 x2x101例例3：计算：计算211.Idxx=解：由解：由Newton-Leibniz公式得公式得21l 20 69314 18d11ln20.69314718.Idxx=由梯形公式由梯形公式1 110 75I（）由梯形公式由梯形公式，0.75,2 21I+=（）由由公式公式1 111(4)0 6944I由由Simpson公式公式，1 111(4)0.6944,6 11.52I+=11111111(1 33)0.69375,84/35/32I+=由由Newton公式，公式，由由公式公式0.693175.I 由由Cotes公式公式，3N C公式的截断误差和代数精度公式的截断误差和代数精度梯形公式梯形公式3、N-C公式的截断误差和代数精度公式的截断误差和代数精度梯形公式梯形公式()()bRffdT f()()TaRff x dx T f=()bf()()()2!bafx a x b dx=()()()2!bafx a x b dx=(积分中值定理)(积分中值定理)3()().(,)12b afa b=12Simpson公式公式()()()()H af a H bf b=构造三次构造三次Hermite插值多项式插值多项式H(x)，满足，满足3333()(),()(),()(),()().H af a H bf bababababHfHf=+=33()(),()()2222ff3()()()()bbSRfI fS ff x dxHx dx=3()()()()Saaffff2(4)()()()bfabxaxxb dx+=()()4!2axaxxb dx2(4)()()()bfabb d+55(4)(4)1()baba()()()4!2afxaxxb dx=(4)(4)1()()().(,)9022880babaffa b=梯形公式的代数精度梯形公式的代数精度一次代数精度一次代数精度22211()()bbbfddb当当(),f xx=bbaba22211()(),22baaaf x dxxdxxba=()()()(),22bababaT ff af babf x dx=+=+=2()f当当233311()()bbbfddb2(),f xx=当当2333()(),33baaaf x dxx dxxba=bbaba22()()()(),22bababaT ff af babf x dx=+=+Simpson公式的代数精度公式的代数精度三次代数精度三次代数精度()f xx=，当当221()(),2bbaaf x dxxdxba=2()4()()babaS ff aff b+=+1baab+()4()()62S ff aff b=+221(4)().622baababba+=+=2()f xx=，当当2331()(),3bbaaf x dxx dxba=3()4()()babaS ff aff b+=+()4()()62S ff aff b+222baab+222(4()62baabab+=+22331(242)().63baaabbba=+=32441bb3()f xx=，当当2441()(),4bbaaf x dxx dxba=()4()()62babaS ff aff b+=+333(4()62baabab+=[phone](33)baaaa babbb=+31ba(33)62aaa babbb=[phone]()().624baaa babbba=+=定理定理：含有：含有n个求积节点的插值型求个求积节点的插值型求积公式至少具有积公式至少具有1次次代数精度代数精度积公式至少具有积公式至少具有n-1次次代数精度代数精度。定理定理：对于：对于N-C公式，当公式，当n为奇数时为奇数时至少具有至少具有 次代数精度当次代数精度当 为偶数为偶数至少具有至少具有n次代数精度次代数精度，当当n为偶数为偶数时至少有时至少有n+1次代数精度。次代数精度。注实际计算中般只用低阶注实际计算中般只用低阶N C公式公式注注：实际计算中实际计算中一一般只用低阶般只用低阶N-C公式公式（高阶不稳定）（高阶不稳定）3x例例5.用用n=2和和n=3的的N-C公式近似计算公式近似计算321.edx解：解：n=2时，时，12332x3222212(4)0.766575505,6edxeee+=n=3时，时，571332x336622212(33)0.766916279,8xe dxeeee+=321(0.7668010.)xedx15 3 复化求积公式复化求积公式5.3 复化求积公式复化求积公式从余项的形式可以看到，积分区间越小，求积从余项的形式可以看到，积分区间越小，求积公式的截断误差越小。因此，可以把积分区间分成若干个小区间，在每个小区间上采用低次公式的截断误差越小。因此，可以把积分区间分成若干个小区间，在每个小区间上采用低次求积公式（梯形公式或抛物线公式）计算积分，然后相加得到整个区间上的积分近似值，这就求积公式（梯形公式或抛物线公式）计算积分，然后相加得到整个区间上的积分近似值，这就是是复化求积复化求积的基本思想。的基本思想。ba1、复化梯形公式1、复化梯形公式3hh把区间（把区间（a,b）n等分，等分，,0,.ibahxaih inn=+=131()()()(),212iixiiixhhf xf xf xf+=+3110()()()()212nbiiiaihhf x dxf xf xf+=+0212i=311()2()()()nniihhf af xf bf=+()()TfRf10()2()()()212iiiif af xf bf=+()().nnTTfRf=+331nhnhba20()()()().121212nTiihnhbaRfffh f=,0,()0,().TnhRfTfI时收收敛敛,0,()0,().nTnnhRfTfI时敛敛0 lim,phR fCh=22,sin10,24(2)962TRf hnn=80.n4(4)4(4)()/2()ibaRf hh fh4(4)4(4)(),()sin,28802880SRf hh fh=451,10,57602SRf hh57602/20 315 1h0.315.1,0.31hn6.n总结：总结：梯形求积公式和抛物线求积公式是低精度方法，但梯形求积公式和抛物线求积公式是低精度方法，但对于光滑性较差的函数有时比用高精度方法能得到对于光滑性较差的函数有时比用高精度方法能得到对于光滑性较差的函数有时比用高精度方法能得到对于光滑性较差的函数有时比用高精度方法能得到更好的效果。复化梯形公式和更好的效果。复化梯形公式和 Simpson公式，精公式，精度较高，计算简单，使用非常广泛。度较高，计算简单，使用非常广泛。度较高，计算简单，使用非常广泛。度较高，计算简单，使用非常广泛。5.4 变步长求积法变步长求积法实际应用中，往往事先很难估计合适的划分数实际应用中，往往事先很难估计合适的划分数 n，使结果达到预期精度。，使结果达到预期精度。因为分点的加密会改善精度因此可以采用因为分点的加密会改善精度因此可以采用因为分点的加密会改善精度因为分点的加密会改善精度，因此可以采用因此可以采用自动加密分点的方式，并利用事后估计判断自动加密分点的方式，并利用事后估计判断精度是否足够从而停止计算精度是否足够从而停止计算精度是否足够精度是否足够，从而停止计算从而停止计算。n等分区间（等分区间（a,b）,1()()2()=+nhTf af bf x1()()2(),2=+niiTf af bf x1211210()()2()2()4+=+nnniiiihTf af bf xf x11nh104=ii1201().22+=+niihTf x2()ba2()()()(),12nbaI fTfh f=2()bh2()()()(),122nbahI fTff=1()(),ff若有若有()221()()()().3nnnI fTfTfTf事后误差估计事后误差估计()21()()3nnTfTf时可停止计算.时可停止计算.3例：计算10.1dxI=01x+1)分别用梯形公式和Si公式计算并估计误差1)分别用梯形公式和Simpson公式计算并估计误差;2)用n=5的复化梯形公式计算并估计误差;3)用变步长梯形公式计算，使其误差小于510.1111)(0)(1)(1)0.75.222Tff+=+=11111252221111125(0)4()(1)(14)0.69444.6263/2236Sfff=+=+=()33331212(1 0)(1 0)0.166667,12(1)12(10)R T=+()5512 3 412 3 4 12(1)12(10)+()555512 3 412 3 4(1 0)(1 0)0.008333.2880(1)2880(10)R S =+12)0.2,5bahn=令5n0,(0,1,2,3,4,5),ixihi=+=构造节点有：5111111()2()2 1 01 0 21 0 41 0 61 0 81 1hT h=+2 1 01 0.21 0.41 0.61 0.81 10.21111110.695635,20 60 70 80 92+=+()221 021 02|0 20 20 066667R T20.60.70.80.92()22533|0.20.20.066667.12(1)12(10)R T=[phone]),1T T T T TT计算521103nnTT直至为止,12n1/201T().22nniibaTf xn+=+其中,022in=数值微分数值微分当函数当函数 f(x)以离散列表形式给出，或者函数以离散列表形式给出，或者函数f(x)过于复杂时，要求用数值的方法计算节点过于复杂时，要求用数值的方法计算节点f(x)过于复杂时，要求用数值的方法计算节点过于复杂时，要求用数值的方法计算节点处的导数值。处的导数值。()()()()()lilif xhf xf xf xhf+由导数的定义：由导数的定义：00()()()()()limlim()()lihhfffffxhhf xhf xh=+0()()lim,2hffh=个自然且简单的方法取极限的近似值即差商个自然且简单的方法取极限的近似值即差商一一个自然且简单的方法个自然且简单的方法：取极限的近似值取极限的近似值，即差商即差商.xfhxf)()(+向前差商向前差商hxfhxfxf)()()(000+由由Taylor展式展式200000()()()(),2!hf xhf xhfxfxxh+=+2!因此，有误差因此，有误差000()()()()()().2!f xhf xhR xfxfO hh+=2!hhff)()(向后差商向后差商hhxfxfxf)()()(000由由Taylor展式展式200000()()()(),2!hf xhf xhfxfxhx=+y2!因此，有误差因此，有误差000()()()()()().2!f xf xhhR xfxfO hh=2!h中心差商中心差商hxfhxfxf)()()(000+中心差商中心差商hxf2)(0由由Tl展式展式2300001010()()()()(),hhf xhf xhfxfxfxxh+=+由由Taylor展式展式000010102300002020()()()()(),2!3!()()()()(),2!3!fffffhhf xhf xhfxfxfxhx=+00002020()()()()(),2!3!fffff()()fhfh因此，有误差因此，有误差00022()()()()2f xhf xhR xfxhhh+=22212()()()().126hhfffO h=+=注：由误差表达式，注：由误差表达式，h越小，误差越小，但同时越小，误差越小，但同时舍入误差增大所以有个最佳步长舍入误差增大所以有个最佳步长舍入误差增大舍入误差增大，所以所以，有个最佳步长有个最佳步长.事后误差估计：设事后误差估计：设D(h),D(h/2)分别为步长为分别为步长为h,h/2的差商公式。则的差商公式。则h,h/2的差商公式。则的差商公式。则)()(hDhD)2()(DhD时的步长时的步长h/2就是合适的步长就是合适的步长时的步长时的步长h/2就是合适的步长就是合适的步长.例例12：f(x)=exp(x).采用中心差分格式，采用中心差分格式，hf(1.15)R(x)hf(1.1R(x)()()(5)()0.103.1630-0.00480.053.1590-0.00080.093.1622-0.00400.043.1588-0.00060 083 1613-0 00310 033 1583-0 00010.083.1613-0.00310.033.1583-0.00010.073.1607-0.00250.023.1575-0.00070.063.1600-0.00180.013.1550-0.0032插值型微分公式：插值型微分公式：用插值函数用插值函数p(x)的导数近似原函数的导数近似原函数f(x)的导数的导数()()f()()nfxpx()()f一一般只考虑节点上般只考虑节点上误差误差()().inifxpx般只考虑节点上般只考虑节点上导数的近似值导数的近似值误差误差()(1)()()()()(),nddfR xf xpxx+=()()()()(),(1)!nnR xf xpxxdxdxn+(1)()nf+()()().(1)!inifR xxn=+两点公式两点公式(n=1)()01001()()(),()(),2hfxf xf xR xfh=两点公式两点公式()()()0100110121()()(),()().hhfxf xf xR xf=()1101()()(),()()2ffffh三点公式三点公式(n=2)2(3)1()(34)()()hfxyyyR xf=+=三点公式三点公式(n=2)001202(3)()(34),()(),231()()()()fxyyyR xfhhfRf=+=+(3)10212(3)()(),()(),261()()()()fxyyR xfhhff=+=(3)201201()(43),()().23hfxyyyR xfh=+=例例13：设，取：设，取h=0.05，用三点公式，用三点公式计算计算 f(2)的近似值。的近似值。()lnf xx=计算计算 f(2)的近似值。的近似值。解解()1(2)3(2)4(2.05)(2.10)0.49980286,ffff+=解解：()()()()(),2 0.05ffff()1(2)(1.95)(2.05)0.50010421,fff+=()()()(),2 0.05fff()1(2)(1.90)4(1.95)3(2)0.49977938,ffff+=()(2)(1.90)4(1.95)3(2)0.49977938,2 0.05ffff+真值真值 真值真值 f(2)=0.5.