数值积分与数值微分21599学习教案

1、数值数值(shz)积分与数值积分与数值(shz)微分微分21599第一页，共51页。第1页/共51页第二页，共51页。第2页/共51页第三页，共51页。1 引言引言(ynyn)一、数值一、数值(shz)求积的基本思想求积的基本思想1、牛顿(ni dn)-莱布尼兹公式但是求函数f(x)的原函数F(x)不一定比计算积分容易，例如函数 找不到用初等函数表示的原函数。另外若给出的函数f(x)是数据表，也不好求函数的积分。计算定积分的方法：计算定积分的方法：第3页/共51页第四页，共51页。2、积分中值定理( )() ( )baf x dxba f但是点 的具体位置一般不知道，故难以准确算出 的值。有两

2、种近似(jn s)方法：得到(d do)梯形(txng)公式另一种是矩形法矩形法，左矩形公式一种是梯形法梯形法，用 代替中矩形公式右矩形公式简称矩形公式简称矩形公式第4页/共51页第五页，共51页。定义定义(dngy)：公式：公式叫做叫做(jiozu)数值求积公式（机械求积）数值求积公式（机械求积）其中其中xk称为求积节点，称为求积节点， Ak称为求积系数称为求积系数(xsh),亦称伴随节点亦称伴随节点xk的权。的权。第5页/共51页第六页，共51页。二、二、 代数代数(dish)精度的概精度的概念念定义定义若某个求积公式对于次数不超过若某个求积公式对于次数不超过m的多项式均能准确成的多项式均

3、能准确成立，但立，但对于对于m+1m+1次的多项式次的多项式就不能准确成立，则称此求积公式的就不能准确成立，则称此求积公式的代数精度代数精度为为 m 。一般要使求积公式0( )()nbkkakfx dxA fx具有m次次精度，只要令它使2( )1, ,.,mf xx xx都能准确成立，即：第6页/共51页第七页，共51页。解：逐次检查解：逐次检查(jinch)公式是否精确成立公式是否精确成立当当 f(x)=1：当当 f(x)=x ：当当 f(x)= x2 ： 故：代数故：代数(dish)精精度度 = 1例例1 1 考察考察1112( ) ( 1)2 (0)(1)f x dxfff有几次代数精度

4、。有几次代数精度。第7页/共51页第八页，共51页。三、三、 插值型的求积公式插值型的求积公式(gngsh)设给定一组节点:012naxxxxb且已知函数在这些节点处的函数值f(xi)(i=0,1,n),由第二章知可以作插值函数Ln(x),由于Ln(x)为多项式，所以其积分可以很容易求得：( )bkkaAlx dx记记 则上式则上式= =0()()nbkkakAfxfx dx称公式：称公式：其中其中 为求积系数。为求积系数。定义定义0( )()nbkkaknf x dxA fI fIxf( )bkkaAlx dx(k=0,1,.,n)为插值型的求积公式。为插值型的求积公式。其余项为：其余项为：

5、(1)( )( )(1)!nbnafR fIIx dxn第8页/共51页第九页，共51页。证明：证明：显然对于次数不超过显然对于次数不超过n次的多项式，次的多项式，Rf等于等于0，所以说明，插值型求，所以说明，插值型求积至少是积至少是n次代数次代数(dish)精度的精度的至少具有n次代数(dish)精度，所以用插值基函数lk(x)当作(dn zu)f(x)代入，上式精确成立，即：所以为插值型的求积定理定理 含有n+1个节点的插值型求积公式 至少具有n次代数精度的充分必要条件是它是插值型的。0( )()nnkkkIfA f x假设第9页/共51页第十页，共51页。例例2 对于对于(duy)a,

6、b上一次上一次插值，有插值，有即即 。考察。考察(koch)其代数其代数精度。精度。f(x)abf(a)f(b)解：逐次检查解：逐次检查(jinch)公式是否精公式是否精确成立确成立代入代入 P0 = 1：=代入代入 P1 = x ：=代入代入 P2 = x2 ： 故：代数精度故：代数精度 = 1则积分公式则积分公式第10页/共51页第十一页，共51页。四、四、 求积公式求积公式(gngsh)有收敛性与稳定性有收敛性与稳定性在求积公式0( )()nbkkakf x dxA f x中，若1100lim()( ),(max()nbkkiianinkhA f xf x dxhxx 其中则称求积公式是

7、收敛的定义定义定义定义(dngy)()kkkf xfkfk在求积公式0( )()nbkkakf x dxA f x中，由于计算 f(xk)可能产生误差，实际得到即：00()()()nnnkknkkkkIfA fxIfA f，如果对任给的小正数0，只要误差充分小就有0()()()nnnkkkkIfIfAfxf则称求积公式是稳定的记第11页/共51页第十二页，共51页。求积公式0( )()nbkkakfx dxA fx中的系数k（k=0,1,n)则此求积公式是稳定的定理定理证明：若取ba,则当()kkfxf有000()()()()()nnnkkkknnkkkkkkIfIfAfxfAfxfAba所以

8、求积公式是稳定的第12页/共51页第十三页，共51页。一、柯特斯系数一、柯特斯系数(xsh)注：注：Cotes 系数系数(xsh)仅取决于仅取决于 n 和和 i，可查表得到。与，可查表得到。与 f (x) 及区间及区间a, b均无关。均无关。 2 Newton-Cotes 公式公式(gngsh)为为n阶阶Newton-cotes求积公式。其中：求积公式。其中：定义：等距节点下的插值型求积公定义：等距节点下的插值型求积公式式为为cotes系数。系数。第13页/共51页第十四页，共51页。21,21)1(1)1(0 CCn = 1:)()(2)(bfafabdxxfba 梯形梯形(txng)公式公

9、式代数代数(dish)精度精度 = 1n = 2:61,32,61)2(2)2(1)2(0 CCC)()(4)(6)(2bffafabdxxfbaba 辛普森公式辛普森公式(gngsh)代数精度代数精度 = 3余项余项余项余项第14页/共51页第十五页，共51页。由书中表知，当时柯特斯系数出了负值(f zh)，所以故时时Newton-Cotes 公式公式(gngsh)不适用。不适用。n = 4:(4)(4)(4)(4)(4)012347162167,9045154590CCCCC01234( )7 () 32 ( ) 12 () 32 ( ) 7 ()90bab af x dxf xf xf

10、xf xf x柯特斯公式柯特斯公式(gngsh)代数精度代数精度=5余项余项第15页/共51页第十六页，共51页。二、偶数阶求积公式的代数二、偶数阶求积公式的代数(dish)(dish)精度精度证明： 当n 为偶数时，由于有(1)1 (1)( )()(1)!nnnfxxn(1)0( ) ( )( )( )()(1)!nnbbnjaajfR fI fIfx dxxx dxn引进变量 x=a+th,并且显然有x j=a+jh2220002/ 222/ 22()()2()0nnnnnnnjjnnnnjnnRfhtj dthuj duhtj dt n 为为偶数阶偶数阶的的Newton-Cotes 公式

11、至少有公式至少有 n+1 次代数精度。次代数精度。定理定理t=u+n/2第16页/共51页第十七页，共51页。几种几种(j zhn(j zhn) )低阶求积公式的余项低阶求积公式的余项梯形梯形(txng)公式：公式：辛普森公式辛普森公式(gngsh)：构造三次插值函数H(x)满足条件：则该函数满足三次代数精度，故辛普森公式柯特斯公式：柯特斯公式：第17页/共51页第十八页，共51页。1 1、 复化梯形复化梯形(txng)(txng)公式公式将区间a,b划分为n等分，分点xk=a+kh,h=(b-a)/n,k=0,1,n，在每个小区间xk,xk+1上采用梯形公式(gngsh)计算.称为(chn

12、wi)复化梯形公式其误差为 3 复化求积公式复化求积公式第18页/共51页第十九页，共51页。由于(yuy)f(x)属于C2a,b,且所以(suy)存在即有复化梯形求积公式是稳定(wndng)的和收敛的。=O(h=O(h2 2) )第19页/共51页第二十页，共51页。2 2 复化辛普森公式复化辛普森公式(gngsh)(gngsh)若记则得：第20页/共51页第二十一页，共51页。其余(qy)项与复化梯形公式相似有：显然(xinrn)有：而由其系数为正数，也为稳定(wndng)的，收敛阶数为4阶，梯形公式的收敛阶数为2阶的。=O(h4)4(4)( )( )( )1802nba hbaRffhn

13、 2( )( )12nbabaRfh fhn 31( )( )121baR fh fh 4(4) ( )1802babaR fh fh 所以复化辛普生求积公式是收敛的。第21页/共51页第二十二页，共51页。例例1 对于对于(duy)函数函数，给出n=8的函数表，试用(shyng)复化梯形公式和复化辛普生公式(gngsh)求积分：xif(xi)011/80.99739782/80.98961583/80.97672674/80.95885105/80.93615566/80.90885167/80.877192510.8414709Excel求解求解第22页/共51页第二十三页，共51页。例2

14、 试用复化辛普森公式求积分：1204xIdxxfunction S=FSimpson(f,a,b,N) h=(b-a)/N; fa=feval(f,a); fb=feval(f,b); S=fb+fa; x=a; for i=1:N x=x+h/2; fx=feval(f,x); S=S+4*fx; x=x+h/2; fx=feval(f,x); S=S+2*fx; end S=h*S/6;function f=f1(x)f=x/(4+x2);f=f1;a=0;b=1;N=64;S=FSimpson(f,a,b,N)第23页/共51页第二十四页，共51页。一、梯形一、梯形(txng)(txng

15、)法的递推化法的递推化为了提高求积精度(jn d)，可在复化求积的基础上将积分小区间xk , xk+1二分一次，增加了一个分点x k+1/2= (x k+x k+1) /2,记h=(b-a)/n，则：4 龙贝格求积公式龙贝格求积公式(gngsh)称之为梯形法的递推化称之为梯形法的递推化第24页/共51页第二十五页，共51页。二、龙贝格算法二、龙贝格算法(sun (sun f)f)由复化梯形(txng)公式的余项知：有即：所以说明(shumng)复化梯形公式二分前后两次计算值的线性组合就为复化辛普森求积公式。第25页/共51页第二十六页，共51页。同理对辛普生公式进行二分处理(chl)，前后两次

16、计算值的误差进行比较，复化柯特斯公式复化柯特斯公式(gngsh)有即：第26页/共51页第二十七页，共51页。例3、用加速公式加工例2得到(d do)的梯形值，计算结果如：kT2kS2k-1C2k-2R2k-301230.92073550.93979330.94451350.94569090.94614590.94608690.94608830.94608300.9460831 0.9460831重复(chngf)同样的操作，可进一步得到龙贝格（Romberg)公式：第27页/共51页第二十八页，共51页。三、理查森外推加速三、理查森外推加速(ji s)(ji s)法法设 ，则有24212()

17、llThIhhh其中系数l与h无关。( ) , f xCa b定理定理4 4记Tn=T(h) ，则T2n=T(h/2) 第28页/共51页第二十九页，共51页。同理由(lyu)可得：如此(rc)下去，每加速一次，误差的量级便提高2阶。第29页/共51页第三十页，共51页。一般(ybn)地：上述处理(chl)方法称为理查森外推加速法。m=1,2,表示加速(ji s)次数第30页/共51页第三十一页，共51页。设 表示(biosh)二分次后求得的梯形值：计算计算(j sun)步骤：步骤：k表示(biosh)二分次数，m表示(biosh)加速次数表示表示序列 的m次加速值。则：第31页/共51页第三

18、十二页，共51页。例如例如(lr):(lr):地球卫星轨道是一个椭圆，椭圆周长的计算公式是地球卫星轨道是一个椭圆，椭圆周长的计算公式是这里这里a a是椭圆的半径轴，是椭圆的半径轴，c c是地球中心与轨道是地球中心与轨道(gudo)(gudo)中心（椭圆中心）中心（椭圆中心）的距离，记的距离，记h h为近地点距离，为近地点距离，H H为远地点距离，为远地点距离，R=6371R=6371（kmkm）为）为地球半径，则地球半径，则 我国第一颗地球我国第一颗地球(dqi)(dqi)卫星近地点距离卫星近地点距离h=439(km)h=439(km)，远，远地点距离地点距离H=2384(kmH=2384(k

19、m）。试求卫星轨道的周长误差不超过）。试求卫星轨道的周长误差不超过10-510-5。第32页/共51页第三十三页，共51页。01.56464011.5646461.56464821.5646461.5646461.564646即人造卫星即人造卫星(rnzowixng)轨道的周长为轨道的周长为48708km从而从而(cng r)有有解：解：()0kT第33页/共51页第三十四页，共51页。5 5 高斯高斯( (o s)o s)求求积公式积公式定义：若一组节点定义：若一组节点 x0 xn a,b,是使插值型求积公式是使插值型求积公式具有具有2n+1次代数精度次代数精度(jn d)。这样的节点称为。

20、这样的节点称为Gauss 点，点，Ak称为称为Gauss系数，求积公式称为系数，求积公式称为Gauss 型求积公式。型求积公式。一、高斯一、高斯( (o s)o s)求积的求积的一般理论一般理论节点节点 x0 xn 以及系数以及系数 A0 An 都作为待定系数。都作为待定系数。要使求积公式具有要使求积公式具有2n+1 次代数精度，令次代数精度，令 f (x) = 1, x, x2, , x2n+1 代代入求积公式精确成立，解出入求积公式精确成立，解出xk和和 Ak .第34页/共51页第三十五页，共51页。例：例：求求 的的 2 点点 Gauss 公式。公式。dxxfx)(10 解：解：设设

21、，应有，应有 3 次代数精度。次代数精度。 101100)()()(xfAxfAdxxfx代入代入 f (x) = 1, x, x2, x3不是线性方程组，不是线性方程组，不易不易(b y)求解。求解。从求积过程从求积过程(guchng)知需求解非线性方程组，可以利用知需求解非线性方程组，可以利用正交多项式的特性来构造求积公式。正交多项式的特性来构造求积公式。第35页/共51页第三十六页，共51页。证明证明(zhngmng)： “” x0 xn 为为 Gauss 点点, 则公式则公式 至少有至少有 2n+1 次代数精度。次代数精度。 bankkkxfAdxxfx0)()()( 对任意次数不大于

22、对任意次数不大于n 的多项式的多项式 P (x)， P (x) wn+1(x)的次数不的次数不大于大于2n+1，则代入公式，则代入公式(gngsh)应精确成立：应精确成立：0= 0 “” 要证明要证明 x0 xn 为为 Gauss 点，即要证公式对任意次点，即要证公式对任意次数数不大于不大于2n+1 的多项式的多项式 f (x) 精确成立，即证明：精确成立，即证明：0( ) ( )()nbkkakx f x dxA f x设设1( )( ) ( )( )nf xwx P xq x0 x0 xn 为为 Gauss 点点 与任意次数与任意次数不大于不大于n 的多项式的多项式 P(x) （带权）正交

23、（带权）正交。即。即10( )()nnkkwxxx定理定理求积公式(gngsh)是插值型第36页/共51页第三十七页，共51页。 Gauss求积公式的求积系数求积公式的求积系数Ak(k=0,1,n)全为正全为正。定理定理证明证明(zhngmng)：是n次多项式，所以(suy) 为2n次多项式，故高斯(o s)求积公式对其能准确成立，即有：证毕。高斯求积公式是稳定的。推论推论定理定理设则高斯求积公式是收敛的，即有0lim()( )( )nbiianiA fxfxx dx 2( ) , f xCa b第37页/共51页第三十八页，共51页。以勒让德多项式以勒让德多项式Pn+1 (x)的根就是的根就

24、是(jish)求积公式的求积公式的Gauss点，点，相应的求积公式称为相应的求积公式称为Gauss-Legendre 公式。公式。若取若取P2(x)=1/2 (3x2-1)的两个零点的两个零点 做节点构造求积公式做节点构造求积公式13101111( )()()33fx dxA fA f二、高斯二、高斯(o s)-勒让德求积公式勒让德求积公式两点高斯两点高斯(o s)-勒让德求积公式勒让德求积公式(n=1)一点高斯一点高斯-勒让德求积公式（勒让德求积公式（n=0)11()2(0)fx dxf若取若取P1(x)=x的零点的零点x0=0做节点构造求积公式做节点构造求积公式令其对令其对f(x)=1 准

25、确成立，则有：准确成立，则有：A0=2。得到得到一点高斯一点高斯-勒让德求积公式勒让德求积公式101( )(0)f x dxA f第38页/共51页第三十九页，共51页。令其对令其对f(x)=1,x都能准确成立，则有都能准确成立，则有01012 ;110 .33AAAA解得解得A0=1,A1=1,从而得到两点高斯从而得到两点高斯-勒让德求积公式。勒让德求积公式。1111()()()33fx dxff同理可以求得同理可以求得三点高斯三点高斯- -勒让德求积公式勒让德求积公式如下：如下：115158515()()(0)()95995fx dxfff高斯高斯- -勒让德求积公式勒让德求积公式(gng

26、sh)(gngsh)节点和系数见书节点和系数见书中表中表4-74-7。第39页/共51页第四十页，共51页。例例6: 用用4点高斯点高斯-勒让德求积公式勒让德求积公式(gngsh)计算：计算：解解: 先将区间先将区间(q jin)0，/2化为化为-1,1，有：，有：再由表再由表4-74-7中中n=3n=3的节点及系数值的节点及系数值(shz)(shz)可以求得：可以求得：对于一般区间对于一般区间a,b应用高斯求积公式时，先用变量置换：应用高斯求积公式时，先用变量置换：将它转变为将它转变为-1,1上的积分上的积分.第40页/共51页第四十一页，共51页。6 6 数值数值(shz)(shz)微分微

27、分一、中点方法与误差一、中点方法与误差(wch)分析分析其中(qzhng)h称为步长。中点方法中点方法考查f(ah)在x=a处的泰勒公式：按导数定义可以用差商近似导数，如下：第41页/共51页第四十二页，共51页。所以(suy)中点公式有：则有 其中(qzhng)1. 从截断误差的角度看，步长越小，计算结果越精确(jngqu)。2. 从舍入误差上来看，由于当h很小时，f(a+h)与f(a-h)很接 近，所以在中点公式中，出现了两个接近的数相减，会造成有效数字的损失，所以步长h又不能太小。例如，用中点公式求在x=2处的一阶导数。计算结果如书中表4-8,显然在步长为0.1时逼近效果最好，因为当f(

28、a+h)及f(a-h)分别有舍入误差时，f (a)的舍入误差上限为说明步长越小则舍入误差越大.第42页/共51页第四十三页，共51页。一般用中点公式计算(j sun)的误差为：所以(suy)步长不宜太小，也不宜太大，最优步长为：第43页/共51页第四十四页，共51页。二、插值型求导公式二、插值型求导公式(gngsh)用插值多项式来构造(guzo)数值微分的基本方法：对给定的f(x)的函数表，构造(guzo)对应的插值多项式Pn(x),再令 去求函数微商。常用(chn yn)的是在节点处带余项的数值微分公式：(1)1(1)(1)11( )( )( )( )( )(1)!( )( )( )( )(

29、 )( )(1)!(1)!nnnnnnnnnfRxfxPxxnxfdfxPxxfnndx由于1()0nkx分析分析:(1)1( )()()()(1)!nknknkffxPxxn故第44页/共51页第四十五页，共51页。1、两点公式设两个(lin )节点x0 , x1处的函数值为：f(x0), f(x1),求f(x0), f(x1)。线性插值公式(gngsh)为记h= x1 -x0 有余项：第45页/共51页第四十六页，共51页。2、三点、三点(sn din)公式公式设已知三点设已知三点(sn din)x0,x1=x0+h,x2=x0+2h上的函数值上的函数值f(x0), f(x1), f(x2),求求f(x0), f(x1), f(x2) 。对t求导可得：令由故令t=0,1,2，便得到(d do)三节点处的导数：余项：第46页/共51页第四十七页，共51页。同理可以得到同理可以得到(d do)多点插值求导公式，也可以多点插值求导公式，也可以用插值函数求高阶导数。用插值函数求高阶导数。例如(lr)：已知函数y=f(x)的数值如下表：x1.01.