1、北京科技大学数理学院卫宏儒 ,科学与工程计算,数值积分与数值微分,微积分学的创始人:,英国数学家 Newton,德国数学家 Leibniz,引言,在实际问题中，往往会遇到被积函数f(x)的原函数无法用初等函数来表示，或函数只能用表格表示，或有的虽然能用初等函数表示，但太复杂，所以这些情形都需要去建立定积分的近似计算公式。在数值积分方面，最容易得到的是用f(x)的代数插值函数p(x)来代替它，即： 将积分 区间细分，在每小区间内用简单函数代替复杂函数，这是数值积分的基本思想。 对替代函数的要求：1：精度要高。2：计算量要小。,求积系数,求积节点,上述公式称为数值求积公式。其中 仅与,一、代数精度
2、的定义及确定:,定义：若求积公式 ， 对一切不高于m次的多项式都准确成立，而对于m+1次多项式等号不成立，则称此公式的代数精度为m。代数精度越高，公式越精确。代数精度的求法：从 依次验证求积公式是否成立，若第一个不成立的等式是 ，则该求积公式的代数精度就是 。,二、牛顿柯特斯求积公式（等距结点）,将积分区间a,b n等分，其节点xk为xk=a+kh, k=0,1,2,n式中h=(ba)/n。在n+1个节点上建立插值于f(x)的n次代数多项式（拉格朗日插值公式）Pn(x)， 并引进变换则有,于是插值型积分公式为:,(1.1),这里:,几个常用的Newton-Cotes公式,柯特斯系数,下面分别考
3、虑几种特殊请况。,（一）梯形公式,若积分区间x0,x1两端点处的函数值f0,f1为已知，可应用线性插 值公式P1(x)在区间x0,x1上的积分来近似,这就是(1.1)式中n=1的 情况。当n =1时，C0(1)=1/2，于是有(1.2)(1.2)式称为梯形公式。积分的这种近似计算方法称为梯形法则。 它的几何意义是用四边梯形x0 ABx1的面积(x1x0)(f0+f1)/2代替曲边梯形的面积 。,当n=1时，为梯形公式:,（二）辛普生(Simpson)公式,如果已知步长的三个等距节点x0x1x2处的函数值f0、f1和f2， 则可应用二次插值公式(抛物线插) P2(x)值在区间x0,x2上进行积分
4、。这就是牛顿柯特斯求积公式(1.1)中n=2的情况。 这里，C0(2) =1/6 , C1(2) =2/3 , C2(2) =1/6 代入(1.2)式，可得(1.3)式中h=(x2 -x0)/2, 它通常称式为辛普生公式或抛物线公式。它的几何意义是用抛物线y=P2(x)围成的曲边梯形面积。,代替由y=f(x)围成的曲边梯形面积图(1.2)。这种用抛物线弧来逼近的数值积分方法称为辛普生法则。类似的，如在积分区间x0,x3内取n=3的等间距插值多项式P3(x)进行积分，则将得到(1.4)此公式也称为Simpson（辛普生3/8公式）。n=4时所得到的求积公式称为柯特斯（Cotes）公式。,当n=2
5、时，为抛物线公式,等距结点求积公式的余项问题,误差分析,1：梯形公式的误差取决于插值多项式P1(x)的误差。 根据 ,R1(x)为插值余项,可得于是得到梯形公式的截断误差为,积分广义中值定理,上式中，(xx0)(xx1)在区间x0,x1内不会改变符号， 如果f“(x)在区间内连续，则由积分的广义中值定理，在x0,x1内至少存在一个点，使得所以，带误差项的梯形公式是,因此(1.8) 上式表明，可以用三次多项式G3(x)来估计误差。由于G3(x)是满足插值条件(1.9) 的任意三次多项式，所以A要由另一个附加条件来确定。,由于我们只讨论误差，故不需要知道G3(x)的具体函数形式。 因此做辅助函数显
6、然上式符合,于是，由式(1.8)得到因子(xx0)(xx1)2(xx2)在区间x0,x2内不会变号，故可以应用广义中值定理，即在x0,x2内存在，使,3、Newton-Cotes公式的代数精度 定理： 由(n+1)个相异节点x0 、x1 、x n构造的求积公式的代数精度至少为n。 证明：记Ln(x)为x0,x1,x2.xn的Lagrange 插值多项式,即因为 则有,由于Newton-Cotes公式是特殊情形(等距节点),他的代数精度至少是n,还可以证明当n为偶数时Newton-Cotes公式的代数精度至少是n+1. 例如 当n=1时，R1=（ba)/12f(n), anb。它的代数精度为1。
7、 当n=2时，R2=(ba )5/2880f(4)(n), anb 。它的代数精度为3初步看来似乎n值越大，精确度越大。是不是n越大越好呢？答案是否定的。,N-C求积公式的数值稳定性,例2：用Newton-Cotes公式计算解：当n取不同值时，计算结果如下所示。I准=0.9460831,0.9270354 0.9461359 0.9461109 0.9460830 0.9460830,近似结果,n,1 2 3 4 5,复合求积公式,从余项的讨论看到，积分区间越小，也可使求积公式的截断误差变小。因此，我们经常把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，
8、构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式，这就是复合求积公式的基本思想。 复合求积公式克服了高次Newton-Cotes公式计算不稳定的问题，其运算简单且易于在计算机上实现。 常用的复合求积公式是复合梯形公式和复合抛物线公式。,复合梯形公式、复合抛物线公式,复合梯形公式把区间a,b n等分，取节点xk=a+kh,k=0,1,.n, h=(b-a)/n,对每个小区间xk,xk+1用梯形求积公式，再累加起来得：（1） 复合抛物线Simpson公式令n=2m,m为正整数，在每个小区间x2k-2,x2k上用抛物线求积公式，有：（2）,内点,奇点,偶点,把一个大的积分区间a,b分
9、成若干小区间，在每个小区间上用次数不高的插值公式，如梯形公式和辛普生公式。构造出相应的积分公式，然后再把它们加起来得到整个区间上的求积公式。为了方便计算，通常把大区间等分成若干小区间。,1、复合积分的有关思想,2、有关定义*步长将积分区间a,b n等分，每一小区间的长度为 h=(b-a)/n ,h即为步长。*复合梯形求积公式将积分区间a,b n等分，h=(b-a)/n,则,为复合梯形公式。,3、有关定理 定理1：复合梯形求积公式的误差：,其中h=(b-a)/n,证明：将积分区间a,bn等分，并在每个小区间上直接用梯形求积公式的误差估计，则,定理2： 复合辛普生求积公式的误差：,证明：将积分区间
10、a,bn等分，并在每个小区间上直接用辛普生求积公式的误差估计。,这样就得证。,4、有关方法：自动选取积分步长我们看到，截断误差随着n的增大而减小，如何确定n的值，使误差在允许范围之内。利用误差公式来估计虽可以，但实际运用上很困难。自动选取积分步长，可在求积过程中，根据精度要求，自动确定n ，并算出近似值。具体方法如下。（假设使用复合积分公式）,用复化梯形公式计算令h=1/8=0.125，n=8用复化抛物线计算令h=1/8=0.125,m=4,n=8,三、外推法与Romberg求积公式由低代数精度求积公式构造高阶代数精度求积公式的方法。在区间a,b上，利用复合梯形公式，求T1=(b-a)f(a
11、)/2+f(b)/2 ， T2=T1/2+ (b-a)fa+(b-a)/2令S1=pT1+qT2，p、q为待定系数，由于T1、T2对一次多项式精确成立，我们可以确定p、q使S1对f(x)=x,x2精确成立，为简单记，在0，1上考虑。,Richardson外推法,假设有一个量F*，用一个步长为h的函数F1(h)去逼近，F*与h无关，既有h 0 时， F1(h) F*，其与F1(h)的截断误差有估计式：R(F*)=F*-F1(h)=a1h p1+a2hp2+a k p k+ (1)p kp k-1p2p10,a i(i=1,2, )都是与h无关的常数，也就是说，F1(h)逼近F*的阶是hp1 ，现
12、在提出的问题是能否通过构造出一个新的序列，它逼近F*的阶要比hp1更高，如为hp2 。,将(1)中的h用qh来代替， q0， 则有 F*-F1(qh)=a1(qh)p1+a2(qh)p2+a k(qh)pk+ 现在用hp1乘(1)的两边后和上式相减,整理得(1-qp1)F*-(F1(qh)-qp1F1(h) =a2 (qp2 - qp1)hp2+a k(q p k- qp1 ) h p k+ 因为(1-qp1)不等于零，用(1-qp1)除等式两边有,F* -,F1(qh)-qp1F1(h),1-qp1,=a2(2)hp2 +a k(2)h pk+ （2）,a2(2)= , . . . ,a k
13、(2) = , . . . 都是与h无关的常数，令F2(h)= (3)那么F2(h)逼近F*的误差由(2)知道为hp2. 依次做下去,计算公式为F m+1(H)= ，m=1,2, . . . (4),F1(qh)-qp1F1(h),1-qp1,F m(qh)-q p m F m(h),1-q p m,其中：,用归纳法容易证明，由(4)得到的Fm(h)逼近F*的误差为 上面的这种方法，称为Richardson外推法。,Romberg求积公式,注： 外推法不只用在积分的计算，也可用于计算其他量，只要能把该量看成f(h)当 h 0时的极限为f(0),而f(0)不易直接计算，要通过hi序列对应的函数值
14、序列f(hi)来外推。如果f(h)有类似的展式，即f(h)= 则一切的算法和讨论可类似进行，更一般地设为f(h)=其中 0r1r2rm，ri不必是整数，I是不依赖于h的常数，0= limf(h)，可以类似地推出外推计算公式。,例：单位圆内接正n边形的面积为 S= sin 设h=2/n,当n 有h 0可以类似地推出外推计算公式。 解:将S与h的关系写成 S(h)= 。 当h 0时,S(h)的极限是，即单位圆的面积。将S(h)作Taylor展开： ,即S(h)有如上展式的形式，注意是无理数，所以可按Romberg算法计算积分的过程来计算到任意精度。 取n=6,12,24,即h=1/3,1/6,1/
15、12,这时S(h)值很容易计算。经过计算S(h)及两次外推法，可计算出S(0)= 的近似值。 取更多的h值外推过程还可继续下去。我国古代数学家刘徽在公元263年就是用这种“割圆术”加上特殊的外推技巧，计算得3.1416。,例1：取e=0.00001,用龙贝格方法计算积分,例2:分别用不同方法计算如下积分,并做比较,令I=各种做法比较如下：（1）、Newton-Cotes公式当n=1时：即用梯形公式，I=0.9270354当n=2时：即用Simpson公式，I=0.9461359当n=3时：I=0.9461090当n=4时：I=0.9460830当n=5时：I=0.9460831,（2）用复化梯
16、形公式：令h=1/8=0.125（3）用复化抛物线令h=1/8=0.125,（4） Romberg公式（P218）K T1 T2 T3 T40 0.9207355 1 0.9397933 0.94614592 0.9445135 0.9460869 0.94008303 0.9456906 0.9460833 0.9460831 0.9460831,比较,此例题的精确值为0.9460831.由例题的各种算法可知： 对Newton-cotes公式，当n=1时只有1位有效数字，当n=2时有3位有效数字，当n=5时有7位有效数字。 对复化梯形公式有2位有效数字，对复化抛物线公式有6位有效数字。 用复
17、合梯形公式，对积分区间0，1二分了11次用2049个函数值，才可得到7位准确数字。 用Romberg公式对区间二分3次，用了9个函数值，得到同样的结果。,外推法在数值微分中的应用,总结,1：梯形求积公式和抛物线求积公式是低精度的方法，但对于光滑性较差的函数有时比用高精度方法能得到更好的效果。复化梯形公式和抛物线求积公式，精度较高，计算较简，使用非常广泛。2：Romberg求积方法，算法简单，当节点加密提高积分近似程度时，前面的计算结果可以为后面的计算使用，因此，对减少计算量很有好处。并有比较简单的误差估计方法。3:外推法不仅用于计算积分的近似值，而且可以作为一般方法用于类似问题的解决，比如用来
18、进行微分的数值计算等。,作业9：,牛顿(1642 1727),伟大的英国数学家 , 物理学家, 天文,学家和自然科学家.,他在数学上的卓越,贡献是创立了微积分.,1665年他提出正,流数 (微分) 术 ,次年又提出反流数(积分)术,并于1671,年完成流数术与无穷级数一书 (1736年出版).,他,还著有自然哲学的数学原理和广义算术等 .,莱布尼兹(1646 1716),德国数学家, 哲学家.,他和牛顿同为,微积分的创始人 ,他在学艺杂志,上发表的几篇有关微积分学的论文中,有的早于牛顿,所用微积分符号也远远优于牛顿 .,他还设计了作乘法的计算机 ,系统地阐述二进制计,数法 ,并把它与中国的八卦联系起来 .,
