《数值积分算法与MATLAB实现毕业论文设计.doc》由会员分享，可在线阅读，更多相关《数值积分算法与MATLAB实现毕业论文设计.doc（63页珍藏版）》请在三一办公上搜索。
1、毕业设计（论文）设计（论文）题目：数值积分算法与MATLAB实现摘 要在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。数值积分就是解决此类问题的一种行之有效的方法。积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。除了研究这些数值积分算法的理
2、论外，本文还将这些数值积分算法在计算机上通过MATLAB软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。【关键词】数值积分 牛顿-科特斯求积公式 高精度求积公式 MATLAB软件ABSTRACTWhen the solution of the definite integral of some function values，because the original function is very complex and difficult to find the elementary function expression, the integral is
3、 difficult to accurately calculate, only managed to find the approximate value, and the case is small that allows to direct interface with the Newton - Leibniz formula to calculate the definite integral. Numerical integration is an effective method to solve such problems. The numerical integration i
4、s an important branch of numerical analysis; therefore, exploring the approximate calculation of the numerical integration method has obvious practical significance. This article departure from the numerical integration problem, described in detail some important numerical integration methods.This p
5、aper has introduced detail the Newton - Coates quadrature formula, and in order to improve the calculation accuracy of numerical integration formulas, More precise formulas have Romberg quadrature formulas and the Gauss - Legendre quadrature formula. In addition to the study of these numerical integ
6、ration algorithm theory, the article also involve what these numerical integration algorithm be programmed by matlab software on the computer, and an example is calculated with a variety of quadrature formulas, finally analysis and comparison to various quadrature formulas calculation error.【Key wor
7、ds】 Numerical integration Newton-Cotes quadrature formula High-precision quadrature formula Matlab software目 录前 言1第一章 牛顿-科特斯求积公式2第一节 数值求积公式的构造2第二节 复化求积公式9第三节 本章小结12第二章 高精度数值积分算法13第一节 梯形法的递推13第二节 龙贝格求积公式14第三节 高斯求积公式17第四节 高斯-勒让德求积公式19第五节 复化两点高斯-勒让德求积公式22第六节 本章小结23第三章 各种求积公式的MATLAB编程实现与应用24第一节 几个低次牛顿-科
8、特斯求积公式的MATLAB实现24第二节 复化求积公式的MATLAB实现28第三节 龙贝格求积公式的MATLAB实现33第三节 高斯-勒让德求积公式的MATLAB实现34第五节 各种求积算法的分析比较36第六节 本章小结38结 论39致 谢40参考文献41附 录43一、英文原文43二、英文翻译52前 言对于定积分，在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式可以计算定积分的值，但在很多情况下的原函数不易求出或非常复杂。被积函数的原函数很难用初等函数表达出来，例如等；有的函数的原函数存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式。因此能够借助牛顿-莱布尼兹公式计算定
9、积分的情形是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值。因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算。而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值。微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节。数值积分是数学上重要的课题之一，是数值分析中重要的内容之一，也是应用数学研究的重点。随着计算机的出现，近几十年来，对于数值积分问题的研究已经
10、成为一个很活跃的研究领域。现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义。国内外众多学者在数值积分应用领域也提出了许多新方法。在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等。通过这个课题的研究，我们将会更好地掌握运用数值积分算法求特殊积分函数的定积分的一些基本方法、理论基础；并且通过matlab软件编程的实现，应用于实际生活中。第一章 牛顿-科特斯求积公式第一节 数值求积公式的构造大多数实际问题的积分是需要用数值积分方法求出近似结果的。数值积
11、分原则上可以用于计算各种被积函数的定积分，无论被积函数是解析解形式还是数表形式，其基本原理都是用多项式函数近似代替被积函数，用多项式的积分结果近似代替对被积函数的积分。由于所选多项式形式的不同，可以有许多种数值积分方法。而利用插值多项式来构造数值求积公式是最常用的一种方法。对于积分，用一个容易积分的函数去代替被积函数，这样的自然以多项式为最佳，因为多项式能很好的逼近任何连续函数，而且容易求出其原函数。一、 求积公式的推导在积分区间上取有限个点，作的次插值多项式，其中，为次插值基函数。用近似代替被积函数，则得 若记 则得数值求积公式 其中称为求积系数，称为求积节点。则称该求积公式为插值型求积公式
12、。知道了插值型求积公式以及其构造方法。为了便于计算与应用，常将积分区间的等分点作为求积节点，这样构造出来的插值型求积公式就称为牛顿-科特斯（Newton-Cotes）求积公式。在积分区间上取个等距节点，其中，做次拉格朗日插值多项式，因为，所以 记 截去第二项得 显然与无关，只与节点有关。令，则当时，于是 而 从而得记 则 故求积公式可写成 这就是牛顿-科特斯求积公式，其中称为科特斯系数。部分科特斯系数取值如下表1.1科特斯系数具有以下特点1（1） （2） （3）当 8 时，出现负数，稳定性得不到保证。而且当 较大时，由于Runge现象，收敛性也无法保证。故一般不采用高阶的牛顿-科特斯求积公式。
13、 （4）当 7 时，牛顿-科特斯公式是稳定的。表1.1 部分科特斯系数表知道了什么是牛顿-科特斯求积公式，下面我们来看它的误差估计，首先来看看牛顿-科特斯求积公式的截断误差。我们知道牛顿-科特斯求积公式是一个插值型数值求积公式，当用插值多项式代替进行积分时，其截断误差即积分真值和近似值之差，推导如下，由插值多项式的误差估计可知，用次拉格朗日多项式逼近函数时产生的误差为 其中。对上式两边从到作定积分，便可得出它的截断误差 二、几个低次牛顿-科特斯求积公式从上面的讨论可知，用多项式近似代替被积函数进行数值积分时，虽然最高次数可是8，但是8次多项式的计算式非常繁杂的。常用的是下面介绍的几种低次多项式
14、。1、 矩形求积公式定义1.1 在牛顿-科特斯求积公式中，如果取，用零次多项式（即常数）代替被积函数，即用矩形面积代替曲边梯形的面积，则有 称式为矩形求积公式根据牛顿-科特斯求积公式的误差理论式，矩形求积公式的误差估计为2、梯形求积公式定义1.21 在牛顿-科特斯求积公式中，如果取，用一次多项式代替被积函数，即用梯形面积代替曲边梯形的面积，则有其中，,查表可得代入上式得出 称式为梯形求积公式由于用一次多项式近似代替被积函数，所以它的精度是1。也就是说，只有当被积函数是一次多项式时，梯形求积公式才是准确的。根据牛顿-科特斯求积公式的误差理论式，梯形求积公式的误差估计为是被积函数二阶导数在点的取值
15、，3、辛浦生求积公式定义1.32 在牛顿-科特斯求积公式中，如果取，用二次多项式代替被积函数，即曲边用抛物线代替，则有其中，,查表可得，代入上式得出 称式为辛浦生求积公式，也称抛物线求积公式。它的几何意义是：用过3个点，的抛物线和，构成的曲边梯形面积，近似地代替了被积函数形成的曲边和，构成的曲线梯形面积。下面对辛浦生求积公式的误差进行估计。由于辛浦生求积公式是用二次多项式逼近被积函数推得的，原则上它的代数精度为2.但因多项式次数是偶数，根据定理1.1可知，它的代数精度为3过，和3个点，构造一个的三次Lagrange插值多项式,且使。根据Lagrange插值余项定理得 对上式两边从到进行积分，即
16、可得到 根据定积分中值定理可知，在上总有一点满足下述关系：通过变量代换,，很容易求得把这个结果代入式，便得出辛浦生求积公式的误差估计式 4、科特斯求积公式由于和时具有相同的迭代精度，但是时计算量小，故的Newton-Cotes积分公式用的很少。定义 1.4 在牛顿-科特斯求积公式中，如果取时，牛顿科特斯公式为 称式为科特斯求积公式。同理根据式可求得其误差估计式 三、求积公式的代数精度如果被积函数为任意一个次数不高于次的多项式时，数值求积公式一般形式的截断误差；而当它是次多项式时，则说明数值求积公式具有次代数精度。一个数值求积公式的代数精度越高，表示用它求数值积分时所需逼近被积函数的多项式次数越
17、高。定理 1.13 牛顿-科特斯求积公式的代数精度等于，当为偶数时，牛顿-科特斯求积公式的代数精度等于证明 如果被积函数是一个不大于次的多项式，则，即；而当时任意一个次多项式时，故 所以，按照代数精度的定义可知，一般情况下，牛顿-科特斯求积公式的代数精度等于。当为次多项式时， 牛顿-科特斯求积公式的代数精度至少等于。若设是一个次多项式，这时为一常数，而因此，只要证明在是偶数时，定理就可得证。为此，设，令，于是由于为偶数，不妨设，为正整数，则，于是再引进变换，则，代入上式右侧，得出 最后的积分中被积函数是奇函数，所以积分结果等于零，定理1.1得证。第二节 复化求积公式前面导出的误差估计式表明，用
18、牛顿-科特斯公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点，亦即提高插值多项式的次数。龙格现象表明，这样做并不一定能提高精度。理论上已经证明，当时，牛顿-科特斯公式所求得的近似值不一定收敛于积分的准确值，而且随着的增大，牛顿-科特斯公式是不稳定的。因此，实际中不常用高阶牛顿-科特斯公式。为了提高计算精度，可考虑对被积函数用分段低次多项式插值，由此导出复化求积公式。一、复化梯形求积公式在实际应用中，若将积分区间分成若干个小区间，在各个小区间上采用低次的求积分式（梯形公式或辛浦生公式），然后再利用积分的区间可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复化积分公式
19、的基本思想。以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。这样求得近似值显然比用梯形公式计算精度高。定积分存在定理表明，只要被积函数连续，当小区间的长度趋于零时，小梯形面积之和即就趋于曲边梯形面积的准确值，即定积分的准确值。定义 1.54 将积分区间进行等分，记为,在每个小区间 上用梯形公式求和，得若将所得的近似值记为，整理得 称式为复化梯形公式。记为 当时， 即收敛于如果，则在小区间上，梯形公式的截断误差为 因此因为区间上连续，由介值定理知存在，使得从而有 这就是复化梯形公式的截断误差。下面讨论复化梯形公式的数值稳定性。设计算函数值时产生的误差为，则用式计算结果的误差为因此
20、，无论为多大，复化梯形公式都是稳定的。二、复化辛浦生求积公式如果用分段二次插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就可导出复化Simpson公式。定义1.65 将积分区间分成等分，分点为， 在每个小区间上。用Simpson公式求积分，则有 求和得 整理后得到 式就称为复化辛浦生求积公式。记为 如果， 则由Simpson插值余项公式可得复化公式的截断误差为 因为为连续，故存在, 使得代入上式得 式表明，步长越小，截断误差越小。与复化梯形公式的分析相类似，可以证明，当时，用复化Simpson公式所求得的近似值收敛于积分值，而且算法具有数值稳定性。三、复化科特斯求积公式
21、定义1.7 将积分区间等分为个子区间，每个子区间的中点,，子区间长度， 在每个子区间上用科特斯公式求和，得 式就称为复化科特斯求积公式，记为 ，式中， ， 类似地可以推出复化科特斯公式的截断误差为 第三节 本章小结本章节开篇介绍了数值求积公式的构造，主要是用运用插值多项式。接着介绍了几个低次的牛顿-科特斯求积公式，即梯形公式、辛浦生公式、科特斯公式，以及牛顿-科特斯求积公式的改进复化求积公式，并对各个求积公式进行了相应的误差分析。第二章 高精度数值积分算法复化求积公式是提高精确度的一种行之有效的方法，但是在使用复化型求积公式之前必须先给出步长。步长太大精度难以保证，步长太小则又会导致计算量的增
22、加，而事先给出一个合适的步长往往是困难的。在实际计算中常常采用变步长的计算方法，即在步长逐次减半的过程中，反复利用复化求积公式进行计算，并同时查看相继两次计算结果的误差是否达到要求，直到所求得的积分值满足要求为止。下面以梯形公式为例第一节 梯形法的递推在变步长的过程中探讨梯形法的计算规律如下：设将积分区间分为等分，则一共有个等分点，这里用表示复化梯形法求得的积分值，其下标表示等分数。由余项公式可知，积分值为 再将各子区间分半，使得区间成等分。此时所得积分近似值记为，则再由余项公式可知，积分值为 假定在上变化不大，即有，于是得，左式也可以写成为 这说明用作为积分的近似值时，其误差近似为。计算过程
23、中常用是否满足作为控制计算精度的条件。如果满足，则取作为的近似值；如果不满足，则再将区间分半，直到满足要求为止。实际计算中的递推公式为 在给定控制参数后，当满足时，则以作为积分的近似值。通过类似的推导，还可以得到下面的结论对于辛浦生公式，假定在上变化不大，则有 对于科特斯公式，假定在上变化不大，则有 第二节 龙贝格求积公式梯形法的算法简单，单精度低，收敛的速度缓慢。如何提高收敛速度以节省计算量，这是人们极为关心的课题。由此引出了龙贝格公式。由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果作为积分的近似值，其误差近似值为。可以设想，如果用这个误差作为的一种补偿，即将 作为积分的
24、近似值，可望提高其精确度。直接根据复化求积公式，不难验证 这说明，将区间对分前后两次复化梯形公式的值，按式作线性组合恰好等于复化辛浦生公式的值，它比更接近于近似值。同样，根据式用于作线性组合会得到比更精确的值，且通过直接验证可得 再由式用与作线性组合，又可得到比更精确的值，通常记为,即 式就称为龙贝格求积公式2。上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推方法。我们将序列,和分别称为梯形序列、辛浦生序列、科特斯序列和龙贝格序列。由龙贝格序列当然还可以继续进行外推，得到新的求积序列。但由于在新的求积序列中，其线性组合的系数分别为与。因此，新的求积序列与前一个序列结果相差不大。
25、故通常外推到龙贝格序列为止。利用龙贝格序列求积的算法称为龙贝格算法。这种算法具有占用内存少、精确度高的优点。因此，成为实际中常用的求积方法。下面给出龙贝格求积算法的计算步骤：第1步：算出和的值，根据公式计算；第2步：将分半，算出后，根据公式 计算,再根据公式计算；第3步：再将区间分半，算出及的值，并根据公式和计算出及，再由公式计算出；第4步：将区间再次分半，计算，并由公式 计算；第5步：将区间再次分半，类似上述过程计算，。重复上述过程可计算得到，一直算到龙贝格序列中前后两项之差的绝对值不超过给定的误差限为止。定义2.16 上述用若干个积分近似值算出更精确的积分近似值的方法，称之为外推法。上述计
26、算步骤也可用表 2.1表示 表 2.1 龙贝格计算步骤表K梯形序列辛浦生序列科特斯序列龙贝格序列012345第三节 高斯求积公式前面介绍的个节点的 Newton -Cotes求积公式，其特征是节点是等距的。这种特点使得求积公式便于构造，复化求积公式易于形成。但同时也限制了公式的精度。是偶数时，代数精度为，是奇数时，代数精度为；我们知道个节点的插值型求积公式的代数精确度不低于。设想：能不能在区间上适当选择个节点 使插值求积公式的代数精度高于？答案是肯定的，适当选择节点，可使公式的精度最高达到，这就是本节所要介绍的高斯型求积公式。例：其中，固定在，可以适当选取，只有两个自由度，得到的是梯形公式，其
27、代数精度只有1。如果对求积节点也进行适当的选取，将有四个自由度，得到如下公式：这个积分公式的代数精确度为3，这就是高斯型求积公式，上面的求积节点称为高斯点。一、高斯型求积公式和高斯点 对于含有个参数，的求积公式：，适当选取这个参数，可以使得数值积分公式的代数精确度达到，我们称这一类求积公式为高斯型求积公式，称这类求积公式的积分节点为高斯点，系数称为高斯系数。定义 2.27 如果个求积节点的求积公式的代数精度为，则这个求积节点称为高斯点。因为高斯求积公式也是插值型求积公式，故有结论 ：个节点的插值型求积公式的代数精度，满足二、高斯点的特征定理 2.17 设是个相异点，以这个点为零点的次多项式为则
28、是高斯点的充要条件是对于任意不超过次的多项式，成立证明 1）必要性。设为高斯点，则对任意不超过次的多项式，均有，则对于任意次数不超过次多项式，是次数不超过的多项式，且注意到，故 2）充分性。设对于一切次数不超过次的多项式，成立，又设是次数不超过次的多项式，用去除，商，余，即，可知和均是不超过次的多项式，从而，又因求积公式是插值多项式的构造导出的，由的选取，其代数精确度可以达到，而是次数不超过次的多项式，因此成立，因，所以，故而也即由于是次数不超过次的多项式，因此该积分公式的代数精确度至少为，因而由定义2.2知节点是高斯点。定义2.3 对于关系，我们称之为正交性，即与任意次多项式正交，而这样的多
29、项式类称为正交多项式。三、高斯型求积公式的余项其中，。积分，可以计算得到 第四节 高斯-勒让德求积公式常用的高斯型求积公式有高斯-勒让德公式、高斯-切比雪夫公式、高斯-拉盖尔公式、高斯-埃尔米特公式，下面着重介绍高斯-勒让德公式一、Legendre多项式， 称式为勒让德（Legendre）多项式，其具有前面提到的正交性性质，即对于任意次数不超过的多项式，成立因此，多项式的零点就是相应的高斯-勒让德求积公式的高斯点。勒让德多项式的前几项如下：，勒让德多项式的性质：性质1 勒让德多项式的首项系数为 性质 2 当为奇数时，为奇函数；当为偶数，为偶函数性质 3 对一切次数不高于次的多项式，有二、高斯-
30、勒让德求积公式1）当时，其零点为，易得，的高斯-勒让德求积公式是，其代数精确度为12）当时，其零点为，此时，的高斯-勒让德求积公式是：，其代数精确度为3 3）当时，其零点为和0，设高斯-勒让德求积公式为：令，得到方程 令，得到方程 令，得到方程 联立方程、可解得，所以，其代数精确度为5 4）当时，由，得，故四个零点为，即，相应的，可以解出， ，高斯-勒让德求积公式为其代数精确度为7综上，高斯-勒让德求积公式的积分节点和积分系数表如 表2.2表2.2 积分节点和积分系数表三、高斯-勒让德求积公式的余项 对阶求积公式（共有个求积节点） 其中，为多项式的零点。积分，可以计算得到： 四、一般区间上的高
31、斯-勒让德求积公式前面讲解的高斯-勒让德求积公式是计算标准区间上的定积分，对于一般的有界区间上的定积分，可以通过变量代换转化为区间的定积分。即在积分中，令，当时，；当时，；且。这个变换为线性变换，其反变换为于是成立然后，我们可以用区间上的高斯-勒让德求积公式。第五节 复化两点高斯-勒让德求积公式本节对两点公式 进行推广和复化，得到了新的数值积分公式 同时分析了它的积分误差和收敛阶定义 2.48 设，是一种复化求积公式，若当时成立 则称求积公式是阶收敛的。例如，复化的Trapezoid公式和复化的Simpson公式分别具有二阶和四阶收敛性。定理 2.28 设，则复化两点高斯-勒让德求积公式为 的
32、余项表达式为 该方法具有四阶收敛性。考虑积分，准确解为 0.74684204用复化两点高斯-勒让德公式可求得=0.74680332 ，其绝对误差为=0.00001960，与复化梯形公式，复化Simpson公式以及原两点Gauss-Legendre 公式相比，这里构造的复化两点高斯-勒让德求积公式方法在运算量和精度方面都有显著的改善。第六节 本章小结本章主要介绍了精确度比较高的两种数值求积公式，即龙贝格求积公式和高斯型求积公式，对其进行了相应的误差分析。其中高斯型求积公式主要介绍了高斯-勒让德求积公式，并对两点高斯-勒让德求积公式进行了复化。第三章 各种求积公式的MATLAB编程实现与应用MAT
33、LAB是由MathWorks公式开发的一种主要用于数值计算及可视化图形处理的工程语言，是当今最优秀的科技应用软件之一。它将数值计算、矩阵运算、图形图像处理、信号处理和仿真等诸多强大的功能集成在较易使用的交互计算机环境中，为科学研究、工程应用提供了一种功能强、效率高的编程工具。下面我们将各种求积算法通过MATLAB软件编程实现，以下程序均用MATLAB7.0编写，运行坏境：1、硬件环境CPU（intel Core i3-2310M,2.1GHz），内存（2GB昱联）,2、软件环境windows7（32位）操作系统。以下总共编写了六个算法程序，部分代码参考文献10-13,为了体现程序的正确性，以下
34、程序都以为例进行运算。原积分的精确值为第一节 几个低次牛顿-科特斯求积公式的MATLAB实现先用M文件定义一个名为f1.m的函数：% i是要调用第几个被积函数g(i)，x是自变量function f=f1(i,x) g(1)=sqrt(x);if x=0 g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x2);f=g(i);程序一：function C,g=NCotes(a,b,n,m)% a，b分别为积分的上下限；% n是子区间的个数；% m是调用上面第几个被积函数；% 当n=1时计算梯形公式；当n=2时计算辛浦生公式，以此类推； i=n; h=(b-a)/i;
35、z=0;for j=0:i x(j+1)=a+j*h; s=1; if j=0 s=s; elsefor k=1:j s=s*k;endendr=1;if i-j=0 r=r;elsefor k=1:(i-j) r=r*k;endendif mod(i-j),2)=1 q=-(i*s*r);else q=i*s*r;endy=1;for k=0:i if k=j y=y*(sym(t)-k); endendl=int(y,0,i);C(j+1)=l/q; z=z+C(j+1)*f1(m,x(j+1);endg=(b-a)*z1）当输入，时，即在MATLAB命令窗口输入 NCotes(0,1,1
36、,2)即可得用梯形公式的积分值和相应科特斯系数 如图3.12）当输入，时，即在MATLAB命令窗口输入 NCotes(0,1,2,2)即可得用辛浦生公式的积分值和相应科特斯系数如图3.23）当输入，时，即在MATLAB命令窗口输入 NCotes(0,1,4,2)即可得用科特斯公式的积分值和相应科特斯系数如图3.3图 3.1 图 3.2 图3.3第二节 复化求积公式的MATLAB实现一、复化梯形求积公式的MATLAB实现通过的个等步长节点逼近积分其中，。程序二：function s=trapr1(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总
37、面积；h=(b-a)/n;s=0;for k=1:(n-1) x=a+h*k; s=s+feval(f,x);endformat long s=h*(feval(f,a)+feval(f,b)/2+h*s;先用M文件定义一个名为f.m的函数：function y=f(x)if x=0 y=1;else y=sin(x)/x;end在MATLAB命令窗口中输入 trapr1(f,0,1,4) 回车得到 如图3.4 图 3.4若取子区间的个数在MATLAB命令窗口中输入 trapr1(f,0,1,8) 回车得到 如图3.5 图3.5二、 复化辛浦生求积公式的MATLAB实现程序三：function
38、 s=simpr1(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积,即所求积分数值；h=(b-a)/(2*n);s1=0;s2=0;for k=1:n x=a+h*(2*k-1); s1=s1+feval(f,x);endfor k=1:(n-1) x=a+h*2*k; s2=s2+feval(f,x);ends=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;先用M文件定义一个名为f.m的函数：function y=f(x)if x=0 y=1;else y=sin(x)/x;end在MATLAB命令窗口中
39、输入 simpr1(f,0,1,4)回车得到 如图3.6图3.6若取子区间个数时在MATLAB命令窗口中输入 simpr1(f,0,1,8)回车得到 如图3.7图3.7三、 复化科特斯求积公式的MATLAB实现程序四：function s=cotespr1(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积,即所求积分数值；h=(b-a)/n;s1=0;s2=0;s3=0;s4=0;for k=1:n x=a+(4*k-3)*h/4; s1=s1+feval(f,x);endfor k=1:n x=a+(4*k-2)*h/4; s2=s2+
40、feval(f,x);endfor k=1:n x=a+(4*k-1)*h/4; s3=s3+feval(f,x);endfor k=1:(n-1) x=a+4*k*h/4; s4=s4+feval(f,x);ends=h*(7*feval(f,a)+7*feval(f,b)+32*s1+12*s2+32*s3+14*s4)/90;同样先用M文件定义一个名为f.m的函数：function y=f(x)if x=0 y=1;else y=sin(x)/x;end在MATLAB命令窗口中输入 cotespr1(f,0,1,4)回车得到 如图3.8 图3.8第三节 龙贝格求积公式的MATLAB实现构造数表来逼近积分其中。表示数表的最后一行，最后一列的值。程序五：function R,quad,err,h=romber(f,a,b,n,delta)% f是被积函数% a,b分别是积分的上下限% n+1是T数表的列数% delta是允许误差% R是T数表% quad是所求积分值M=1;h=b-a;err=1J=0;R=zeros(4,4);R(1,1)=h
