《数值计算方法之数值积分.ppt》由会员分享，可在线阅读，更多相关《数值计算方法之数值积分.ppt（37页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第第12章：数值积分与数值微分章：数值积分与数值微分设设f(x)是是a,b上连续可积的实函数，求上连续可积的实函数，求f(x)在在a,b上的数值积分也就是求上的数值积分也就是求f(x)在在a,b上的定积分的数值上的定积分的数值解。解。即使我们能找到即使我们能找到f(x)的一个原函数的一个原函数F(x)的解析形式的解析形式,并并利用牛顿利用牛顿-来布尼兹公式进行计算，在许多情况下这来布尼兹公式进行计算，在许多情况下这也是非常麻烦的。也是非常麻烦的。问题的关键在于，初等函数的原函数不一定还是初等问题的关键在于，初等函数的原函数不一定还是初等函数，在这种情况下并不能利用牛顿函数，在这种情况下并不能利
2、用牛顿-来布尼兹公式来布尼兹公式进行计算。进行计算。本章的主要任务是寻找求数值积分的最有效的计算方本章的主要任务是寻找求数值积分的最有效的计算方法，最后建议的方法是变步长复化辛卜生加速算法，法，最后建议的方法是变步长复化辛卜生加速算法，为此，需要经过一个曲折的分析过程。为此，需要经过一个曲折的分析过程。12.1 12.1 求积公式与代数精度的概念求积公式与代数精度的概念由定积分的定义可知由定积分的定义可知，连续函数连续函数f(x)在区间在区间a,b上上的定积分近似值可以表示为的定积分近似值可以表示为a,b内的一些点内的一些点x0,x1,xn处的函数值处的函数值f(x0),f(x1),f(xn)
3、的加权和或的加权和或线性组合，即线性组合，即 （1）其中其中w0,w1,wn仅与仅与x0,x1,xn有关而与被积函数有关而与被积函数f(x)无关。我们把这样的公式称为求积公式，也称为无关。我们把这样的公式称为求积公式，也称为机械求积公式。机械求积公式。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
4、为为a,b的一个等份分划那么求积公的一个等份分划那么求积公式式(1)中的中的w0,w1,wn的选取仅仅只与的选取仅仅只与n有关，从而有关，从而可以简化对求积公式的研究。可以简化对求积公式的研究。结论，只要给出了一个如何确定（结论，只要给出了一个如何确定（1）式中的诸）式中的诸w0,w1,wn的机制，我们就可以得到相应的对任何被的机制，我们就可以得到相应的对任何被积函数都有效的计算定积分方法。积函数都有效的计算定积分方法。2.2.求积公式的性质求积公式的性质微积分学中我们曾研究过，定积分保持函数的线性关微积分学中我们曾研究过，定积分保持函数的线性关系不变，它的含义是，若系不变，它的含义是，若f(
5、x),g(x)都是都是a,b上的可上的可积函数，则对任意实数积函数，则对任意实数u,v,我们有我们有uf(x)+vg(x)也是也是a,b上的可积函数，而且上的可积函数，而且不难验证，求积公式也保持函数的线性关系不变，即不难验证，求积公式也保持函数的线性关系不变，即3.3.几种常见的求积公式几种常见的求积公式在后面的讨论中，我们将经常用到下面一些非常简单的求积在后面的讨论中，我们将经常用到下面一些非常简单的求积在后面的讨论中，我们将经常用到下面一些非常简单的求积在后面的讨论中，我们将经常用到下面一些非常简单的求积公式，他们是中点公式、梯形公式和辛卜生公式。公式，他们是中点公式、梯形公式和辛卜生公
6、式。公式，他们是中点公式、梯形公式和辛卜生公式。公式，他们是中点公式、梯形公式和辛卜生公式。中点公式我是我们课程中强调的一个名词，与求积公式中点公式我是我们课程中强调的一个名词，与求积公式中点公式我是我们课程中强调的一个名词，与求积公式中点公式我是我们课程中强调的一个名词，与求积公式（1 1）对比分析，可以认为它是这样一种机制：把积分区间）对比分析，可以认为它是这样一种机制：把积分区间）对比分析，可以认为它是这样一种机制：把积分区间）对比分析，可以认为它是这样一种机制：把积分区间分为分为分为分为2 2等分，取等分，取等分，取等分，取ww0 0=w=w2 2=0=0，ww1 1=1=1所形成的求
7、积公式。从几何所形成的求积公式。从几何所形成的求积公式。从几何所形成的求积公式。从几何上看，它实际上是取区间中点的函数值与区间长的积作为定上看，它实际上是取区间中点的函数值与区间长的积作为定上看，它实际上是取区间中点的函数值与区间长的积作为定上看，它实际上是取区间中点的函数值与区间长的积作为定积分值，类似于用中位线乘以高来计算梯形的面积。积分值，类似于用中位线乘以高来计算梯形的面积。积分值，类似于用中位线乘以高来计算梯形的面积。积分值，类似于用中位线乘以高来计算梯形的面积。4 4截断误差截断误差在求积公式中，我们使用的是近似等号，这是因为，对于一在求积公式中，我们使用的是近似等号，这是因为，对
8、于一在求积公式中，我们使用的是近似等号，这是因为，对于一在求积公式中，我们使用的是近似等号，这是因为，对于一般的被积函数来说，利用这些公式计算所得的结果除了舍入般的被积函数来说，利用这些公式计算所得的结果除了舍入般的被积函数来说，利用这些公式计算所得的结果除了舍入般的被积函数来说，利用这些公式计算所得的结果除了舍入误差外，还有截断误差，因为定积分是用极限来定义的。误差外，还有截断误差，因为定积分是用极限来定义的。误差外，还有截断误差，因为定积分是用极限来定义的。误差外，还有截断误差，因为定积分是用极限来定义的。有时为了进行误差分析，我们可以把上面的（有时为了进行误差分析，我们可以把上面的（有时
9、为了进行误差分析，我们可以把上面的（有时为了进行误差分析，我们可以把上面的（1 1）式写成）式写成）式写成）式写成 （11）其中其中其中其中RfRf表示的就是截断误差。表示的就是截断误差。表示的就是截断误差。表示的就是截断误差。考察前面给出的三个求积公式，如果被积函数是线性函数，考察前面给出的三个求积公式，如果被积函数是线性函数，考察前面给出的三个求积公式，如果被积函数是线性函数，考察前面给出的三个求积公式，如果被积函数是线性函数，那么利用中点公式或梯形公式所得到的结果就是准确值，否那么利用中点公式或梯形公式所得到的结果就是准确值，否那么利用中点公式或梯形公式所得到的结果就是准确值，否那么利用
10、中点公式或梯形公式所得到的结果就是准确值，否则一般不是。对于一般的非线性函数，感觉上辛卜生公式更则一般不是。对于一般的非线性函数，感觉上辛卜生公式更则一般不是。对于一般的非线性函数，感觉上辛卜生公式更则一般不是。对于一般的非线性函数，感觉上辛卜生公式更好一些。好一些。好一些。好一些。为了刻划求积公式对一般的被积函数的精确度，我们引进代为了刻划求积公式对一般的被积函数的精确度，我们引进代为了刻划求积公式对一般的被积函数的精确度，我们引进代为了刻划求积公式对一般的被积函数的精确度，我们引进代数精度的概念。数精度的概念。数精度的概念。数精度的概念。5.5.代数精度的概念代数精度的概念定义：一个求积公
11、式定义：一个求积公式定义：一个求积公式定义：一个求积公式 如果对所有的次数不超过如果对所有的次数不超过如果对所有的次数不超过如果对所有的次数不超过mm的多项式严格相等，而对某些的多项式严格相等，而对某些的多项式严格相等，而对某些的多项式严格相等，而对某些m+1m+1次多项式不相等，则称该公式具有代数精度次多项式不相等，则称该公式具有代数精度次多项式不相等，则称该公式具有代数精度次多项式不相等，则称该公式具有代数精度mm，或该公或该公或该公或该公式的代数精度为式的代数精度为式的代数精度为式的代数精度为mm。利用求积公式的线性性，我们不难证明下面的结论。利用求积公式的线性性，我们不难证明下面的结论
12、。利用求积公式的线性性，我们不难证明下面的结论。利用求积公式的线性性，我们不难证明下面的结论。定理：如果求积公式对定理：如果求积公式对定理：如果求积公式对定理：如果求积公式对1 1，x,xx,xmm严格相等，而对严格相等，而对严格相等，而对严格相等，而对x xm+1m+1不相不相不相不相等，则该公式的代数精度为等，则该公式的代数精度为等，则该公式的代数精度为等，则该公式的代数精度为mm。作为课外练习，鼓励大家给出完整证明。作为课外练习，鼓励大家给出完整证明。作为课外练习，鼓励大家给出完整证明。作为课外练习，鼓励大家给出完整证明。6.6.基本结论基本结论我们可以利用上面的定理所给出的方法证明辛卜
13、生公我们可以利用上面的定理所给出的方法证明辛卜生公式的代数精度是式的代数精度是3，而中点公式和梯形公式的代数精，而中点公式和梯形公式的代数精度是度是1。现在我们可以对这三个公式作一个简单的评价：现在我们可以对这三个公式作一个简单的评价：中点公式和梯形公式的代数精度虽然都是中点公式和梯形公式的代数精度虽然都是中点公式和梯形公式的代数精度虽然都是中点公式和梯形公式的代数精度虽然都是1 1，但中点公式，但中点公式，但中点公式，但中点公式只计算一个点的函数值，而梯形公式却要计算两个点处的只计算一个点的函数值，而梯形公式却要计算两个点处的只计算一个点的函数值，而梯形公式却要计算两个点处的只计算一个点的函
14、数值，而梯形公式却要计算两个点处的函数值，所以中点公式优于梯形公式。函数值，所以中点公式优于梯形公式。函数值，所以中点公式优于梯形公式。函数值，所以中点公式优于梯形公式。与梯形公式相比，辛卜生公式只多计算一个点的函数值，与梯形公式相比，辛卜生公式只多计算一个点的函数值，与梯形公式相比，辛卜生公式只多计算一个点的函数值，与梯形公式相比，辛卜生公式只多计算一个点的函数值，但代数精度却增加到但代数精度却增加到但代数精度却增加到但代数精度却增加到3 3，显然辛卜生公式更为优越。，显然辛卜生公式更为优越。，显然辛卜生公式更为优越。，显然辛卜生公式更为优越。10102 2 牛顿牛顿-柯特斯求积公式柯特斯求
15、积公式牛顿牛顿-柯特斯求积公式就是利用柯特斯求积公式就是利用Lagrange插值多项式插值多项式导出的求积公式。导出的求积公式。把一般的函数的积分转化为相应的插值多项式函数的把一般的函数的积分转化为相应的插值多项式函数的积分也是我们学习插值法的基本目的之一。积分也是我们学习插值法的基本目的之一。1.1.利用插值多项式近似替代被积函数利用插值多项式近似替代被积函数设设f(x)为被积函数，为被积函数，a,b为积分区间，为积分区间，x0,x1,xn为为a,b内的内的n+1个互异的点，记个互异的点，记Ln(x)为相应的拉格朗为相应的拉格朗日插值多项式，那么我们有日插值多项式，那么我们有 f(x)=Ln
16、(x)+Rn(x)2.2.利用插值多项式导出求积公式利用插值多项式导出求积公式提示：在上面给出的公式中，由于诸提示：在上面给出的公式中，由于诸提示：在上面给出的公式中，由于诸提示：在上面给出的公式中，由于诸l lk k(x)(x)都是多项式函数，都是多项式函数，都是多项式函数，都是多项式函数，所以诸所以诸所以诸所以诸wwk k都可以精确地计算出来。从而可以得到一般性的求都可以精确地计算出来。从而可以得到一般性的求都可以精确地计算出来。从而可以得到一般性的求都可以精确地计算出来。从而可以得到一般性的求积公式。积公式。积公式。积公式。2.2.利用插值多项式导出求积公式（注释）利用插值多项式导出求积
17、公式（注释）回顾上一章关于多项式插值的结论，由于任意次数不回顾上一章关于多项式插值的结论，由于任意次数不超过超过n的多项式与它的任意的多项式与它的任意n+1个基点的插值多项式个基点的插值多项式恒等，再由求积公式的代数精度的定义，我们立即得恒等，再由求积公式的代数精度的定义，我们立即得到：到：由由由由n+1n+1个基点的拉格朗日插值多项式所形成的求积公式的个基点的拉格朗日插值多项式所形成的求积公式的个基点的拉格朗日插值多项式所形成的求积公式的个基点的拉格朗日插值多项式所形成的求积公式的代数精度至少式代数精度至少式代数精度至少式代数精度至少式n n，为此，我们上面的为此，我们上面的为此，我们上面的
18、为此，我们上面的wwk k改写为改写为改写为改写为w(n,k),k=0,1,nw(n,k),k=0,1,n。3.3.牛顿牛顿-柯特斯求积公式柯特斯求积公式牛顿牛顿-柯特斯求积公式就是利用等距基点的拉格朗日柯特斯求积公式就是利用等距基点的拉格朗日插值多项式导出的求积公式。插值多项式导出的求积公式。将积分区间将积分区间a,b划分为划分为n等分，记等分，记h=(b-a)/n，取取 x0=a,xk=a+kh.k=0,1,n 我们可以得到我们可以得到 3.3.牛顿牛顿-柯特斯求积公式（注释）柯特斯求积公式（注释）在牛顿在牛顿-柯特斯公式中，我们称柯特斯公式中，我们称c(n,k)为牛顿为牛顿-柯特斯柯特斯
19、系数，一般可通过查表得到。崔国华教材系数，一般可通过查表得到。崔国华教材p60列出了列出了直到直到n=8的所有牛顿的所有牛顿-柯特斯系数，应该说，实用意柯特斯系数，应该说，实用意义不大。义不大。当当n8时牛顿时牛顿-柯特斯公式并没有实际意义。实际上，柯特斯公式并没有实际意义。实际上，我们通常只用到我们通常只用到n=1,2,4的情形，相应的公式分别称的情形，相应的公式分别称为梯形公式，辛卜生公式，和柯特斯公式。对于现代为梯形公式，辛卜生公式，和柯特斯公式。对于现代的计算工具来说，有梯形公式和辛卜生公式也就够用的计算工具来说，有梯形公式和辛卜生公式也就够用了。了。利用牛顿利用牛顿-柯特斯系数，我们
20、可以方便地写出牛顿柯特斯系数，我们可以方便地写出牛顿-柯柯特斯求积公式：特斯求积公式：4.4.梯形公式梯形公式在牛顿在牛顿-柯特斯求积公式中，如果我们取柯特斯求积公式中，如果我们取n=1,那么那么k可可以取以取0和和1。由此所形成的求积公式就是梯形公式。由此所形成的求积公式就是梯形公式。由于得到的结果是梯形面积公式，所以称梯形公式。由于得到的结果是梯形面积公式，所以称梯形公式。5.5.辛卜生公式辛卜生公式在牛顿在牛顿-柯特斯求积公式中，如果我们取柯特斯求积公式中，如果我们取n=2,那么那么k可可以取以取0，1，2，由此所形成的求积公式就是辛卜生公，由此所形成的求积公式就是辛卜生公式。式。6.6
21、.柯特斯公式柯特斯公式作为课外作业，大家可以取作为课外作业，大家可以取作为课外作业，大家可以取作为课外作业，大家可以取n=4,n=4,相应地相应地相应地相应地k k可以取可以取可以取可以取0 0，1 1，2 2，3 3和和和和4 4，仿照上面的方式，可以得到：，仿照上面的方式，可以得到：，仿照上面的方式，可以得到：，仿照上面的方式，可以得到：从而可进一步写出相应的求积公式，这就是柯特斯公式。从而可进一步写出相应的求积公式，这就是柯特斯公式。从而可进一步写出相应的求积公式，这就是柯特斯公式。从而可进一步写出相应的求积公式，这就是柯特斯公式。在后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，在
22、后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，在后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，在后面将要介绍的龙贝格求积算法中，我们将产生梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基辛卜生序列，柯特斯序列和龙贝格序列，前三个序列都是基于牛顿于牛顿于牛顿于牛顿-柯特斯公式产生的序列，而龙贝格序列则不是。柯特斯公式产生的序列，而龙贝格序列则不是。柯特斯公式产生的序列，而龙贝格序列则不是。柯特斯公式产生的序列，而龙贝格序列则不是。10.3 10.3 变步长复化梯形公式
23、法变步长复化梯形公式法在上一节中我们介绍了牛顿在上一节中我们介绍了牛顿-柯特斯公式以及它的特柯特斯公式以及它的特款，并且得到了牛顿款，并且得到了牛顿-柯特斯公式的代数精度为柯特斯公式的代数精度为n.，接下来就是如何利用这个公式。接下来就是如何利用这个公式。一般说来，在一个较大的积分区间上利用较高阶的牛一般说来，在一个较大的积分区间上利用较高阶的牛顿顿-柯特斯公式虽然可以得到较高的代数精度，但实柯特斯公式虽然可以得到较高的代数精度，但实际效果并不好，道理也不难理解。际效果并不好，道理也不难理解。我们可以把积分区间划分为若干等分，在每个子区间我们可以把积分区间划分为若干等分，在每个子区间上利用较低
24、阶的求积公式，上利用较低阶的求积公式，(当然元可以是牛顿当然元可以是牛顿-柯特柯特斯公式，也可以不是，斯公式，也可以不是，)并把这些积分值相加。按照并把这些积分值相加。按照这样的思路所得到的求积公式统称为复化型求积公式。这样的思路所得到的求积公式统称为复化型求积公式。1.1.复化中点公式复化中点公式复化中点公式也许最不为人们所注意，以至在一般的复化中点公式也许最不为人们所注意，以至在一般的教科书中还没有这个名称，我们在后面将会看到，对教科书中还没有这个名称，我们在后面将会看到，对于求数值积分来说，它实际上是最有用的公式。于求数值积分来说，它实际上是最有用的公式。把积分区间把积分区间a,b划分为
25、划分为n等分，记等分，记x0,x1,xn为等分为等分点，记点，记xj-1,xj为第为第j个子区间，个子区间，zj为区间的中点，为区间的中点，j=1,2,n，记记h=(b-a)/n,记记Mn为所有子区间上利用为所有子区间上利用中点公式所求得的积分值的和，那么我们有中点公式所求得的积分值的和，那么我们有 （1）1.1.复化中点公式复化中点公式(等价形式等价形式)现在我们把积分区间现在我们把积分区间现在我们把积分区间现在我们把积分区间 a,ba,b划分为划分为划分为划分为2 2n n等分等分等分等分,记记记记y y0 0,y,y1 1,y,y2n-2n-1 1,y,y2n2n为等分点，此时我们可以把
26、下标为奇数的点为等分点，此时我们可以把下标为奇数的点为等分点，此时我们可以把下标为奇数的点为等分点，此时我们可以把下标为奇数的点y y2k-12k-1,看成是看成是看成是看成是区间区间区间区间 y y2(k-1)2(k-1),y,y2k2k 的中点，的中点，的中点，的中点，k=1,2,nk=1,2,n。(如图如图如图如图)y y0 0 y y1 1 y y2 2 y y3 3 y y4 4 y y2(n-1)2(n-1)y y2n-1 2n-1 y y2n2nx0 x1 x2 xx0 x1 x2 xn-1n-1 x xn n这样，我们也可以把复化中点公式写成这样，我们也可以把复化中点公式写成这
27、样，我们也可以把复化中点公式写成这样，我们也可以把复化中点公式写成 （2 2）提示：这两种表示方法我们在后面都会用到，所以一定要记提示：这两种表示方法我们在后面都会用到，所以一定要记提示：这两种表示方法我们在后面都会用到，所以一定要记提示：这两种表示方法我们在后面都会用到，所以一定要记住。不管公式的形式如何，编写求住。不管公式的形式如何，编写求住。不管公式的形式如何，编写求住。不管公式的形式如何，编写求MMn n的程序还是相同的。的程序还是相同的。的程序还是相同的。的程序还是相同的。1.1.复化中点公式复化中点公式(c c语言代码语言代码)在实际应用中，在实际应用中，在实际应用中，在实际应用中
28、，n n通常取通常取通常取通常取2 2的某个整数幂的某个整数幂的某个整数幂的某个整数幂2 2k k，我们可以把计算我们可以把计算我们可以把计算我们可以把计算结果作为数组结果作为数组结果作为数组结果作为数组MM的第的第的第的第k k个元素，程序段命名为个元素，程序段命名为个元素，程序段命名为个元素，程序段命名为,C,C语言源程序语言源程序语言源程序语言源程序文件可以按如下的方式编写。文件可以按如下的方式编写。文件可以按如下的方式编写。文件可以按如下的方式编写。#include include include include#define F(X)4.0/(1.0+X*X)#define F(X)
29、4.0/(1.0+X*X)#define F(X)4.0/(1.0+X*X)#define F(X)4.0/(1.0+X*X)static double a=0.0,b=1.0;static double a=0.0,b=1.0;static double a=0.0,b=1.0;static double a=0.0,b=1.0;double GetM(int k)double GetM(int k)double GetM(int k)double GetM(int k)double x,y,step;double x,y,step;double x,y,step;double x,y,st
30、ep;int n,j;int n,j;int n,j;int n,j;n=1;n=1;n=1;n=1;for(j=0;jk;j+)n*=2;for(j=0;jk;j+)n*=2;for(j=0;jk;j+)n*=2;for(j=0;jk;j+)n*=2;step=(b-a)/n;step=(b-a)/n;step=(b-a)/n;step=(b-a)/n;x=a+step/2;x=a+step/2;x=a+step/2;x=a+step/2;for(j=0;jn;j+)y+=F(x);x+=step;for(j=0;jn;j+)y+=F(x);x+=step;for(j=0;jn;j+)y+=
31、F(x);x+=step;for(j=0;jn;j+)y+=F(x);x+=step;return(step*y);return(step*y);return(step*y);return(step*y);2.2.复化梯形公式复化梯形公式把积分区间把积分区间a,b划分为划分为n等分，记等分，记x0,x1,xn为等分为等分点，记点，记h=(b-a)/n,则有则有记记Ij,j=1,2,n为为f(x)在第在第j个子区间个子区间xj-1,xj应用梯形应用梯形公式所求得的积分值，则有公式所求得的积分值，则有2.2.复化梯形公式复化梯形公式(续续)这就是我们的复化梯形公式。这就是我们的复化梯形公式。这就是
32、我们的复化梯形公式。这就是我们的复化梯形公式。提示：虽然我们也可以直接编写求提示：虽然我们也可以直接编写求提示：虽然我们也可以直接编写求提示：虽然我们也可以直接编写求T Tn n的计算机程序，但是没的计算机程序，但是没的计算机程序，但是没的计算机程序，但是没有这个必要。有这个必要。有这个必要。有这个必要。3.3.变步长复化梯形公式变步长复化梯形公式假设对某个假设对某个n,我们利用复化梯形公式，也就是上面的我们利用复化梯形公式，也就是上面的（3）式，得到了）式，得到了Tn，如果它不满足我们的精度要求，如果它不满足我们的精度要求，那么我们可以把每个子区间再对分一次，这相当于把那么我们可以把每个子区
33、间再对分一次，这相当于把积分区间划分为积分区间划分为2n等分。等分。记记y0,y1,y2,y2(n-1),y2n-1,y2n为等分点，记为等分点，记t=(b-a)/(2n),则有则有再利用（再利用（3）式即得）式即得3.3.变步长复化梯形公式变步长复化梯形公式(续续)再利用上面的（再利用上面的（再利用上面的（再利用上面的（2 2）式和（）式和（）式和（）式和（3 3）式，并注意到（）式，并注意到（）式，并注意到（）式，并注意到（3 3）式中的诸）式中的诸）式中的诸）式中的诸x xk k即为这里的即为这里的即为这里的即为这里的y y2k2k,k=1,2,n-1,k=1,2,n-1，所以我们有所以
34、我们有所以我们有所以我们有 （4 4）这就是变步长复化梯形公式。我们可以利用它形成一个自动这就是变步长复化梯形公式。我们可以利用它形成一个自动这就是变步长复化梯形公式。我们可以利用它形成一个自动这就是变步长复化梯形公式。我们可以利用它形成一个自动调整精度的算法。调整精度的算法。调整精度的算法。调整精度的算法。4.4.变步长复化梯形公式方法变步长复化梯形公式方法为了方便地利用（为了方便地利用（为了方便地利用（为了方便地利用（4 4）式形成一个算法，我们总是取）式形成一个算法，我们总是取）式形成一个算法，我们总是取）式形成一个算法，我们总是取n n等于等于等于等于2 2的某个指数幂的某个指数幂的某
35、个指数幂的某个指数幂2 2k k,并把并把并把并把T Tn n记为记为记为记为Tk,Tk,由由由由n=2n=2k k可得可得可得可得2 2n=2n=2k+1k+1,从而有从而有从而有从而有 Tk+1=(TK+MK)/2Tk+1=(TK+MK)/2算法说明如下：算法说明如下：算法说明如下：算法说明如下：T0=(b-a)*(f(a)+f(b)/2.0;T0=(b-a)*(f(a)+f(b)/2.0;T0=(b-a)*(f(a)+f(b)/2.0;T0=(b-a)*(f(a)+f(b)/2.0;For(K=0;k20;k+)For(K=0;k20;k+)For(K=0;k20;k+)For(K=0;
36、k20;k+)Mk=GetM(k);Mk=GetM(k);Mk=GetM(k);Mk=GetM(k);Tk+1=(Mk+Tk)/2.0;Tk+1=(Mk+Tk)/2.0;Tk+1=(Mk+Tk)/2.0;Tk+1=(Mk+Tk)/2.0;if(fabs(Mk+1-Mk)EPS)break;if(fabs(Mk+1-Mk)EPS)break;if(fabs(Mk+1-Mk)EPS)break;if(fabs(Mk+1-Mk)EPS)break;提示：如果取提示：如果取提示：如果取提示：如果取k=20k=20k=20k=20，那么那么那么那么GetM(k)GetM(k)GetM(k)GetM(k)
37、中计算函数值的次数超过中计算函数值的次数超过中计算函数值的次数超过中计算函数值的次数超过1000000100000010000001000000次。所以我们一般不让次。所以我们一般不让次。所以我们一般不让次。所以我们一般不让k k k k值超过值超过值超过值超过2020202012.4 变步长复化辛卜生公式法变步长复化辛卜生公式法变步长复化梯形公式法的核心思想是：假设对于变步长复化梯形公式法的核心思想是：假设对于n=2K我们计算出来了我们计算出来了TK,如果如果Tk不满足精度要求，不满足精度要求，则进一步把区间划分为则进一步把区间划分为2K+1等份，也就是在原来的每等份，也就是在原来的每个子区
38、间补上中点的函数值，并利用个子区间补上中点的函数值，并利用Tk+1=(Tk+Mk)/2来计算来计算Tk+1现在我们也可以这样考虑，利用原来每个子区间的端现在我们也可以这样考虑，利用原来每个子区间的端点函数值点函数值f(y2(k-1)和和f(x2k)以及后来补算的中点的函数以及后来补算的中点的函数值值f(y2k-1)采用辛卜生公式计算第采用辛卜生公式计算第k个子区间的积分值，个子区间的积分值，再把他们相加结果计为再把他们相加结果计为SK+1。无论是理论分析还无论是理论分析还是直观考虑，这样做的效果应该更好一些是直观考虑，这样做的效果应该更好一些,问题是能问题是能否也利用否也利用MK,TK来简化我
39、们的计算。来简化我们的计算。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 应用应用辛卜生公式所求得的积分值，则有辛卜生公式所求得的积分值，则有这就是我们的复化辛卜生公式这就是我们
40、的复化辛卜生公式,虽然我们也可以直接虽然我们也可以直接编写求编写求S2n的计算机程序，但是没有这个必要。的计算机程序，但是没有这个必要。2.2.关于复化辛卜生公式的结论关于复化辛卜生公式的结论假设假设x0,x1,x2,x2(k-1),x2k-1,x2k,x2(n-10,x2n-1,x2n把积把积分区间分区间a,b划分为划分为2n等分，那么我们可以得到下面等分，那么我们可以得到下面3个不同的积分值：个不同的积分值：利用利用利用利用x x x x2k2k2k2k,k=0,1,n,k=0,1,n,k=0,1,n,k=0,1,n这这这这n+1n+1n+1n+1个点处的函数值和复化梯形公个点处的函数值和
41、复化梯形公个点处的函数值和复化梯形公个点处的函数值和复化梯形公式计算出式计算出式计算出式计算出T T T Tn n n n;利用利用利用利用x x x xj j j j,j=0,1,2n,j=0,1,2n,j=0,1,2n,j=0,1,2n这这这这2 2 2 2n+1n+1n+1n+1个点处的函数值和复化梯形公个点处的函数值和复化梯形公个点处的函数值和复化梯形公个点处的函数值和复化梯形公式计算出式计算出式计算出式计算出T T T T2n2n2n2n;利用复化辛卜生公式计算出利用复化辛卜生公式计算出利用复化辛卜生公式计算出利用复化辛卜生公式计算出S S S S2n2n2n2n;重要结论：对于上面
42、的约定，我们有重要结论：对于上面的约定，我们有 S2n=(4T2n-Tn)/33 3详细推导详细推导 4.4.变步长复化辛卜生公式加速算法变步长复化辛卜生公式加速算法我们仍然取我们仍然取n为为2的整数幂的形式，即的整数幂的形式，即n=2k,并记并记Sn为为Sk,那么我们有那么我们有Tk+1=(Mk+Tk)/2 Tk+1=(Mk+Tk)/2 Sk+1=(4TK+1-TK)/3.0Sk+1=(4TK+1-TK)/3.0我们可以在变步长复化梯形公式法基础上稍加改进，我们可以在变步长复化梯形公式法基础上稍加改进，即可得到一个相应的加速算法也就是简单地利用即可得到一个相应的加速算法也就是简单地利用TK+
43、1和和TK的线性组合直接得到的线性组合直接得到 Sk+1，额外的，额外的代码和计算量可忽略不计。代码和计算量可忽略不计。5.算法说明算法说明变步长复化辛卜生公式加速算法是一个实用的方法，变步长复化辛卜生公式加速算法是一个实用的方法，（而不是像前面介绍的方法只是过渡性的方法）所以（而不是像前面介绍的方法只是过渡性的方法）所以要求大家都会编程。要求大家都会编程。12.5 12.5 龙贝格求积算法龙贝格求积算法前面我们介绍过变步长复化梯形公式法，利用它可以前面我们介绍过变步长复化梯形公式法，利用它可以形成一个自动调整精度的算法。形成一个自动调整精度的算法。接下来又介绍了如何利用复化辛卜生公式实现变步
44、长接下来又介绍了如何利用复化辛卜生公式实现变步长复化梯形公式法的加速，基本思想是利用复化梯形公式法的加速，基本思想是利用Tk+1和和Tk的线性组合来产生的线性组合来产生Sk+1,我们称为进行一次外我们称为进行一次外推，显然，它比起直接计算推，显然，它比起直接计算Sk+1优越的多。优越的多。类似地，我们还可以利用类似地，我们还可以利用Sk+1和和Sk的线性组合来的线性组合来产生产生Ck+1，称为柯特斯序列；利用称为柯特斯序列；利用Ck+1和和Ck的线性组合来产生的线性组合来产生Rk+1,称为龙贝格序列。显然这称为龙贝格序列。显然这个过程还可以一直继续进行下去。这就是龙贝格算法个过程还可以一直继续
45、进行下去。这就是龙贝格算法的基本思想。的基本思想。1.1.统一的记号与公式统一的记号与公式对于求对于求f(x)在区间在区间a,b上的定积分值，在龙贝格求积上的定积分值，在龙贝格求积算法中，我们采用统一的记号算法中，我们采用统一的记号它们之间的递推关系为它们之间的递推关系为2.2.计算顺序计算顺序对于上面给出的外推序列，我们可按下图所示的顺序进行计对于上面给出的外推序列，我们可按下图所示的顺序进行计对于上面给出的外推序列，我们可按下图所示的顺序进行计对于上面给出的外推序列，我们可按下图所示的顺序进行计算算算算.如果有某个如果有某个如果有某个如果有某个 ,即可停即可停即可停即可停止进一步的计算。止
46、进一步的计算。止进一步的计算。止进一步的计算。提示：为了提示：为了便于编程，便于编程，我们把每行我们把每行的的k值处理得值处理得相同，所以相同，所以约定约定k=m,m+1,，3.3.龙贝格算法的实现方法龙贝格算法的实现方法在龙贝格算法中，我们实际的做法是把梯形序列外推在龙贝格算法中，我们实际的做法是把梯形序列外推3次，即次，即m值取值取0，1，2，3，并把相应的序列称为梯，并把相应的序列称为梯形序列，辛卜生序列，柯特斯序列和龙贝格序列，并形序列，辛卜生序列，柯特斯序列和龙贝格序列，并相应地记为相应地记为Tk,Sk,Ck,Rk。为了使我们的算法为了使我们的算法更为简明，我们再增加一个中点序列，记
47、为更为简明，我们再增加一个中点序列，记为Mk,并并且用一个专门的函数且用一个专门的函数GetM(int k)来计算（返回）来计算（返回）Mk.事实上，我们前面的变步长复化型辛卜生公式加速算事实上，我们前面的变步长复化型辛卜生公式加速算法已经得到了中点序列，梯形序列和辛卜生序列。所法已经得到了中点序列，梯形序列和辛卜生序列。所以稍加改动即可得到龙贝格算法。以稍加改动即可得到龙贝格算法。4.4.实现龙贝格算法的代码实现龙贝格算法的代码具体计算方法是：首先按公式组具体计算方法是：首先按公式组(1)(1)设置设置M0T0,S0,C0,R0,M0T0,S0,C0,R0,接下来再对接下来再对k=0,1,k=0,1,按公按公式组式组(2)(2)计算计算Mk,Tk+1,Sk+1,Ck+1,Rk+1.Mk,Tk+1,Sk+1,Ck+1,Rk+1.