Ch数值计算方法之数值积分实用.pptx

《Ch数值计算方法之数值积分实用.pptx》由会员分享，可在线阅读，更多相关《Ch数值计算方法之数值积分实用.pptx（36页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1 1术语和记号为了计算f(x)在区间a,b上的定积分近似值，我们通常的做法是，把积分区间a,b划分为n等分，记h=(b-a)/n,x0=a,xk=a+kh,k=0,1,2,n,称x0,x1,xn为a,b的一个等份分划。假如x0,x1,xn为a,b的一个等份分划那么求积公式(1)中的w0,w1,wn的选取仅仅只与n有关，从而可以简化对求积公式的研究。结论，只要给出了一个如何确定（1）式中的诸w0,w1,wn的机制，我们就可以得到相应的对任何被积函数都有效的计算定积分方法。第1页/共36页2.2.求积公式的性质微积分学中我们曾研究过，定积分保持函数的线性关系不变，它的含义是，若f(x),g(x)

2、都是a,b上的可积函数，则对任意实数u,v,我们有uf(x)+vg(x)也是a,b上的可积函数，而且不难验证，求积公式也保持函数的线性关系不变，即第2页/共36页3.3.几种常见的求积公式在后面的讨论中，我们将经常用到下面一些非常简单的求积公式，他们是中点公式、梯形公式和辛卜生公式。中点公式我是我们课程中强调的一个名词，与求积公式（1）对比分析，可以认为它是这样一种机制：把积分区间分为2等分，取w0=w2=0，w1=1所形成的求积公式。从几何上看，它实际上是取区间中点的函数值与区间长的积作为定积分值，类似于用中位线乘以高来计算梯形的面积。第3页/共36页4 4截断误差在求积公式中，我们使用的是

3、近似等号，这是因为，对于一般的被积函数来说，利用这些公式计算所得的结果除了舍入误差外，还有截断误差，因为定积分是用极限来定义的。有时为了进行误差分析，我们可以把上面的（1）式写成 （1）其中Rf表示的就是截断误差。考察前面给出的三个求积公式，如果被积函数是线性函数，那么利用中点公式或梯形公式所得到的结果就是准确值，否则一般不是。对于一般的非线性函数，感觉上辛卜生公式更好一些。为了刻划求积公式对一般的被积函数的精确度，我们引进代数精度的概念。第4页/共36页5.5.代数精度的概念定义：一个求积公式 如果对所有的次数不超过m的多项式严格相等，而对某些m+1次多项式不相等，则称该公式具有代数精度m，

4、或该公式的代数精度为m。利用求积公式的线性性，我们不难证明下面的结论。定理：如果求积公式对1，x,xm严格相等，而对xm+1不相等，则该公式的代数精度为m。作为课外练习，鼓励大家给出完整证明。第5页/共36页6.6.基本结论我们可以利用上面的定理所给出的方法证明辛卜生公式的代数精度是3，而中点公式和梯形公式的代数精度是1。现在我们可以对这三个公式作一个简单的评价：中点公式和梯形公式的代数精度虽然都是1，但中点公式只计算一个点的函数值，而梯形公式却要计算两个点处的函数值，所以中点公式优于梯形公式。与梯形公式相比，辛卜生公式只多计算一个点的函数值，但代数精度却增加到3，显然辛卜生公式更为优越。第6

5、页/共36页10102 2 牛顿-柯特斯求积公式牛顿-柯特斯求积公式就是利用Lagrange插值多项式导出的求积公式。把一般的函数的积分转化为相应的插值多项式函数的积分也是我们学习插值法的基本目的之一。第7页/共36页1.1.利用插值多项式近似替代被积函数设f(x)为被积函数，a,b为积分区间，x0,x1,xn为a,b内的n+1个互异的点，记Ln(x)为相应的拉格朗日插值多项式，那么我们有 f(x)=Ln(x)+Rn(x)第8页/共36页2.2.利用插值多项式导出求积公式提示：在上面给出的公式中，由于诸lk(x)都是多项式函数，所以诸wk都可以精确地计算出来。从而可以得到一般性的求积公式。第9

6、页/共36页2.2.利用插值多项式导出求积公式（注释）回顾上一章关于多项式插值的结论，由于任意次数不超过n的多项式与它的任意n+1个基点的插值多项式恒等，再由求积公式的代数精度的定义，我们立即得到：由n+1个基点的拉格朗日插值多项式所形成的求积公式的代数精度至少式n，为此，我们上面的wk改写为w(n,k),k=0,1,n。第10页/共36页3.3.牛顿-柯特斯求积公式牛顿-柯特斯求积公式就是利用等距基点的拉格朗日插值多项式导出的求积公式。将积分区间a,b划分为n等分，记h=(b-a)/n，取 x0=a,xk=a+kh.k=0,1,n 我们可以得到 第11页/共36页3.3.牛顿-柯特斯求积公式

7、（注释）在牛顿-柯特斯公式中，我们称c(n,k)为牛顿-柯特斯系数，一般可通过查表得到。崔国华教材p60列出了直到n=8的所有牛顿-柯特斯系数，应该说，实用意义不大。当n8时牛顿-柯特斯公式并没有实际意义。实际上，我们通常只用到n=1,2,4的情形，相应的公式分别称为梯形公式，辛卜生公式，和柯特斯公式。对于现代的计算工具来说，有梯形公式和辛卜生公式也就够用了。利用牛顿-柯特斯系数，我们可以方便地写出牛顿-柯特斯求积公式：第12页/共36页4.4.梯形公式在牛顿-柯特斯求积公式中，如果我们取n=1,那么k可以取0和1。由此所形成的求积公式就是梯形公式。由于得到的结果是梯形面积公式，所以称梯形公式

8、。第13页/共36页5.5.辛卜生公式在牛顿-柯特斯求积公式中，如果我们取n=2,那么k可以取0，1，2，由此所形成的求积公式就是辛卜生公式。第14页/共36页6.6.柯特斯公式作为课外作业，大家可以取n=4,相应地k可以取0，1，2，3和4，仿照上面的方式，可以得到：从而可进一步写出相应的求积公式，这就是柯特斯公式。在后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基于牛顿-柯特斯公式产生的序列，而龙贝格序列则不是。第15页/共36页10.3 10.3 变步长复化梯形公式法在上一节中我们介绍了牛顿-柯特斯公式以及它的特款，并且得到了牛顿-

9、柯特斯公式的代数精度为n.，接下来就是如何利用这个公式。一般说来，在一个较大的积分区间上利用较高阶的牛顿-柯特斯公式虽然可以得到较高的代数精度，但实际效果并不好，道理也不难理解。我们可以把积分区间划分为若干等分，在每个子区间上利用较低阶的求积公式，(当然元可以是牛顿-柯特斯公式，也可以不是，)并把这些积分值相加。按照这样的思路所得到的求积公式统称为复化型求积公式。第16页/共36页1.1.复化中点公式复化中点公式也许最不为人们所注意，以至在一般的教科书中还没有这个名称，我们在后面将会看到，对于求数值积分来说，它实际上是最有用的公式。把积分区间a,b划分为n等分，记x0,x1,xn为等分点，记x

10、j-1,xj为第j个子区间，zj为区间的中点，j=1,2,n，记h=(b-a)/n,记Mn为所有子区间上利用中点公式所求得的积分值的和，那么我们有 （1）第17页/共36页1.1.复化中点公式(等价形式)现在我们把积分区间a,b划分为2n等分,记y0,y1,y2n-1,y2n为等分点，此时我们可以把下标为奇数的点y2k-1,看成是区间y2(k-1),y2k的中点，k=1,2,n。(如图)y0 y1 y2 y3 y4 y2(n-1)y2n-1 y2nx0 x1 x2 xn-1 xn这样，我们也可以把复化中点公式写成 （2）提示：这两种表示方法我们在后面都会用到，所以一定要记住。不管公式的形式如何

11、，编写求Mn的程序还是相同的。第18页/共36页1.1.复化中点公式(c c语言代码)在实际应用中，n通常取2的某个整数幂2k，我们可以把计算结果作为数组M的第k个元素，程序段命名为,C语言源程序文件可以按如下的方式编写。#include#define F(X)4.0/(1.0+X*X)static double a=0.0,b=1.0;double GetM(int k)double x,y,step;int n,j;n=1;for(j=0;jk;j+)n*=2;step=(b-a)/n;x=a+step/2;for(j=0;jn;j+)y+=F(x);x+=step;return(step

12、*y);第19页/共36页2.2.复化梯形公式把积分区间a,b划分为n等分，记x0,x1,xn为等分点，记h=(b-a)/n,则有记Ij,j=1,2,n为f(x)在第j个子区间xj-1,xj应用梯形公式所求得的积分值，则有第20页/共36页2.2.复化梯形公式(续)这就是我们的复化梯形公式。提示：虽然我们也可以直接编写求Tn的计算机程序，但是没有这个必要。第21页/共36页3.3.变步长复化梯形公式假设对某个n,我们利用复化梯形公式，也就是上面的（3）式，得到了Tn，如果它不满足我们的精度要求，那么我们可以把每个子区间再对分一次，这相当于把积分区间划分为2n等分。记y0,y1,y2,y2(n-

13、1),y2n-1,y2n为等分点，记t=(b-a)/(2n),则有再利用（3）式即得第22页/共36页3.3.变步长复化梯形公式(续)再利用上面的（2）式和（3）式，并注意到（3）式中的诸xk即为这里的y2k,k=1,2,n-1，所以我们有 （4）这就是变步长复化梯形公式。我们可以利用它形成一个自动调整精度的算法。第23页/共36页4.4.变步长复化梯形公式方法为了方便地利用（4）式形成一个算法，我们总是取n等于2的某个指数幂2k,并把Tn记为Tk,由n=2k可得2n=2k+1,从而有 Tk+1=(TK+MK)/2算法说明如下：T0=(b-a)*(f(a)+f(b)/2.0;For(K=0;k

14、20;k+)Mk=GetM(k);Tk+1=(Mk+Tk)/2.0;if(fabs(Mk+1-Mk)EPS)break;提示：如果取k=20，那么GetM(k)中计算函数值的次数超过1000000次。所以我们一般不让k值超过20第24页/共36页12.4 变步长复化辛卜生公式法变步长复化梯形公式法的核心思想是：假设对于n=2K我们计算出来了TK,如果Tk不满足精度要求，则进一步把区间划分为2K+1等份，也就是在原来的每个子区间补上中点的函数值，并利用Tk+1=(Tk+Mk)/2来计算Tk+1现在我们也可以这样考虑，利用原来每个子区间的端点函数值f(y2(k-1)和f(x2k)以及后来补算的中点

15、的函数值f(y2k-1)采用辛卜生公式计算第k个子区间的积分值，再把他们相加结果计为SK+1。无论是理论分析还是直观考虑，这样做的效果应该更好一些,问题是能否也利用MK,TK来简化我们的计算。第25页/共36页1.1.复化辛卜生公式假设x0,x1,x2,x2(k-1),x2k-1,x2k,x2(n-1),x2n-1,x2n把积分区间a,b划分为2n等分，我们可以认为其中的x0,x2,x2k,x2n把区间划a,b划分为n等份，并且x2k-1就是第k个子区间x2(k-1),x2k的中点。记Ik,k=1,2,n为f(x)在第k个子区间x2(k-1),x2k 应用辛卜生公式所求得的积分值，则有这就是我

16、们的复化辛卜生公式,虽然我们也可以直接编写求S2n的计算机程序，但是没有这个必要。第26页/共36页2.2.关于复化辛卜生公式的结论假设x0,x1,x2,x2(k-1),x2k-1,x2k,x2(n-10,x2n-1,x2n把积分区间a,b划分为2n等分，那么我们可以得到下面3个不同的积分值：利用x2k,k=0,1,n这n+1个点处的函数值和复化梯形公式计算出Tn;利用xj,j=0,1,2n这2n+1个点处的函数值和复化梯形公式计算出T2n;利用复化辛卜生公式计算出S2n;重要结论：对于上面的约定，我们有 S2n=(4T2n-Tn)/3第27页/共36页3 3详细推导 第28页/共36页4.4

17、.变步长复化辛卜生公式加速算法我们仍然取n为2的整数幂的形式，即n=2k,并记Sn为Sk,那么我们有Tk+1=(Mk+Tk)/2 Sk+1=(4TK+1-TK)/3.0我们可以在变步长复化梯形公式法基础上稍加改进，即可得到一个相应的加速算法也就是简单地利用TK+1和TK的线性组合直接得到 Sk+1，额外的代码和计算量可忽略不计。第29页/共36页5.算法说明变步长复化辛卜生公式加速算法是一个实用的方法，（而不是像前面介绍的方法只是过渡性的方法）所以要求大家都会编程。第30页/共36页12.5 12.5 龙贝格求积算法前面我们介绍过变步长复化梯形公式法，利用它可以形成一个自动调整精度的算法。接下

18、来又介绍了如何利用复化辛卜生公式实现变步长复化梯形公式法的加速，基本思想是利用Tk+1和Tk的线性组合来产生Sk+1,我们称为进行一次外推，显然，它比起直接计算Sk+1优越的多。类似地，我们还可以利用Sk+1和Sk的线性组合来产生Ck+1，称为柯特斯序列；利用Ck+1和Ck的线性组合来产生Rk+1,称为龙贝格序列。显然这个过程还可以一直继续进行下去。这就是龙贝格算法的基本思想。第31页/共36页1.1.统一的记号与公式对于求f(x)在区间a,b上的定积分值，在龙贝格求积算法中，我们采用统一的记号它们之间的递推关系为第32页/共36页2.2.计算顺序对于上面给出的外推序列，我们可按下图所示的顺序

19、进行计算.如果有某个 ,即可停止进一步的计算。提示：为了提示：为了便于编程，便于编程，我们把每行我们把每行的的k值处理得值处理得相同，所以相同，所以约定约定k=m,m+1,，第33页/共36页3.3.龙贝格算法的实现方法在龙贝格算法中，我们实际的做法是把梯形序列外推3次，即m值取0，1，2，3，并把相应的序列称为梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，并相应地记为Tk,Sk,Ck,Rk。为了使我们的算法更为简明，我们再增加一个中点序列，记为Mk,并且用一个专门的函数GetM(int k)来计算（返回）Mk.事实上，我们前面的变步长复化型辛卜生公式加速算法已经得到了中点序列，梯形序列和辛卜生序列。所以稍加改动即可得到龙贝格算法。第34页/共36页4.4.实现龙贝格算法的代码具体计算方法是：首先按公式组(1)设置M0T0,S0,C0,R0,接下来再对k=0,1,按公式组(2)计算Mk,Tk+1,Sk+1,Ck+1,Rk+1.第35页/共36页感谢您的欣赏！第36页/共36页