数值积分常用算法设计与实现

. 毕 业论文 题 目 数值积分常用算法设计与实现 学 生 * 指导教师 * 教授 年 级 2005级 专 业 计算机科学与技术 系 别 计算机系 学 院 计算机科学技术与信息工程学院*2019年06月 论文提要本文应用插值积分法和逼近论的思想，简单重述了推导复化梯形公式、复化Simpson公式和Newton-Cotes公式的过程，以及这三个公式的系数、精度等问题。并以这两种数值积分的求解方法为基础，应用quad、guass函数编写具体Matlab程序，通过计算机软件计算出所给题目的近似数值积分。对三者所得的结果进行比较，从而研究了用梯形公式、Simpson公式和Newton-Cotes公式求积分的方法和三者的精确度问题。得知：梯形公式的代数精度虽然比较低，但它对函数的光滑性要求也比较低，特别是当被积函数是周期函数时，梯形公式是较理想的方法。Simpson公式具有较高的代数精度，且程序设计也比较简单，因此是使用比较广泛的一种方法。但应当指出，不能指望通过无限增加区间的等分数n来提高数值积分的精确度。而且，当等分数n不恰当的增大时，梯形公式和Simpson公式的计算结果都很差。这两种求积公式所得的结果在精度上的确存在差异，结合理论部分更加充分地说明了，n相同时比Newton-Cotes公式具有更高的代数精度，但当代数精度相同时，三者之间计算的结果仍存在细微的差异。 数值积分常用算法设计与实现 * 摘 要：在数值分析中，给定函数的定积分的计算不总是可行的。许多定积分不能得到精确值的实际问题可以利用数值积分求解。本文应用插值积分法和逼近论的思想，构造一个简单函数p(x)近似表示f(x),然后对 p(x)求积分得到 f(x)的积分的近似值。基于插值原理，推导出数值积分的基本公式。利用Matlab程序设计实现，进而对三大公式运算结果进行分析以及进一步认识。 关键字：插值积分、代数精度、插值型求积公式、复化梯形公式、复化辛普森公式、Newton-Cotes公式。1、数值积分简介： 在数值分析中，数值积分是计算定积分数值的方法和理论。在数学分析中，给定函数的定积分的计算不总是可行的。许多定积分不能用已知的积分公式得到精确值。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解。由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题。数值积分是利用黎曼积分等数学定义，用数值逼近的方法近似计算给定的定积分值。借助于电子计算设备，数值积分可以快速而有效地计算复杂的积分。对微积分学作出杰出贡献的数学大师，如牛顿、.欧拉、高斯等人也在数值积分这个领域作出了各自的贡献，并奠定了它的理论基础。 1、1数值积分公式一般是形如 的近似公式，又称求积公式,xj和Aj(i=0,1,m)分别称为求积结点和求积系数，通常 a,b,式(2)右端称为求积和；两端之差 称为求积余项或求积误差；区间，b可以是有限的或无限的。构造求积公式的问题就是确定xj和Aj使得E()在某种意义下尽可能地小。1、2代数精度若公式对(x)xk(k=0,1,d)精确成立，亦即E()=0，而当(x)=x时该公式不再是确等式，则说该公式的代数精度是d。根据K.外尔斯特拉斯的多项式逼近定理,就一般的连续函数而言,d越大E()越小，因此可以用代数精度的高低说明求积公式的优劣。1、3插值型求积公式通过插值途径构成的求积公式。用(x)的以x0,x1,xm为结点的插值多项式 公式 公式 公式近似替代(x)后,经过积分可以得到形如(2)的插值型求积公式，其中求积系数 公式特别，若所有的xj都属于a,b,则称它为内插型求积公式。这是一类最基本的求积公式。由于m+1个结点的插值型求积公式的代数精度至少是m,所以具有一定代数精度的求积公式总是存在的。2三大常用数值积分公式推导2.1复化梯形公式以及复化辛普森公式的推导： 在区间 a,b不大时,用梯形公式、辛卜生公式计算定积分是简单实用的,但当区间a,b较大时,用梯形公式、Simpson公式计算定积分达不到精确度要求. 为了提高计算的精确度，我们将 a,b 区间n等分，在每个小区间上应用梯形公式、Simpson公式计算定积分, 然后将其结果相加，这样就得到了复化梯形公式和复化Simpson公式。2.11. 复化梯形公式将积分区间 a,b 做n等分,设 , 则节点为 对每个小区间上应用梯形公式, 然后将其结果相加，则得以上公式则为复化梯形公式当 在a,b上有连续的二阶导数时，则复化梯形公式的余项推导如下：因为 所以在区间a,b上公式(3.14)的误差为又因为 在区间a,b上连续，由连续函数的性质知，在区间a,b上存在一点 ，于是 称上式为复化梯形公式的余项。2.12复化Simpson公式与复化梯形公式类似,可以推导出复化Simpson公式.不同之处在于将区间 分为 等分,对每个小区间 应用Simpson公式,然后相加, 则得称为复化Simpson公式当 在a,b上有连续的四阶导数时，则复化辛卜生公式(的余项推导如下：在区间 上辛卜生公式的误差已知为所以，在区间 上公式的误差为又因为 在a,b上连续, 由连续函数的性质知，在区间a,b上存在一点 ，于是 称上式为复化Simpson 公式的余项2.21NewtonCotes公式的推导当1.1插值求积公式的插值节点为等距节点时，就得到NewtonCotes公式。将区间a,bn等分，n+1个节点为xk=a+kh (k=0,1,n)在节点上对f(x)的Lagrange插值多项式是：用Pn(x)代替f(x)构造求积公式：记，(k=0,1,n)作代换x=a+th带入上式，变为：其中： (k=0,1,n) （1-1）这个积分是有理多项式积分，它与被积函数f(x)和区间a,b无关。只要确定n就能计算出系数。于是得到称为NewtonCotes公式的求积公式： （1-2）其中称为NewtonCotes系数。如表1所示。表1 NewtonCotes系数n11/21/221/64/61/631/83/83/81/847/9032/9012/9032/907/90519/28825/9625/14425/14425/9019/288641/8409/359/28034/1059/2809/3541/8402.22NewtonCotes公式误差和稳定性在积分公式中用插值多项式Pn(x)代替f(x)的插值误差是因此，NewtonCotes公式的截断误差是 （1-3）讨论舍入误差对计算结果产生的影响，设（1-2）式近似计算其中计算函数值f(xn)有误差值（k=0,1,2, ,n）。在（1-2）式中令设计算无误差，舍入误差也忽略，则，由（1-2）式计算时引式的误差为如果皆为正，并设，则，故有界，即引起的误差受控制，不超过倍。保证了数值计算的稳定性。但当n8时，将出现负数，这时，数值计算的稳定性不能保证，所以节点超过8时NewtonCotes公式不能用。当n为偶数时，NewtonCotes积分公式具有n+1次代数精度。2.23经典NewtonCotes公式当n=4,5点公式称为经典NewtonCotes公式其中 (k=0,1,4)，它具有5次代数精度。3、实例运算应用与分析3.1 数值分析举例3.2 程序设计流程图 分别利用复化梯形公式、复化Simpson公式、牛顿-科特斯公式对题目求解 ，并做出流程设计图如下 开始0.0 LowLimit1 0.0 UpLimit1输出f(x)=e(-x)请输入积分上限输入UpLimit1输出请输入积分下限输入LowLimit1 tixing1.SetValue(LowLimit1,UpLimit1);/梯形赋值 xinpusheng1.SetValue(LowLimit1,UpLimit1);/辛浦生赋值 输出积分结果：梯形tixing1.Result()输出积分结果: 辛浦生xinpusheng1.Result() 结束 图1 复化梯形和复化辛普森公式流程图 开始输出Input the xi(16): i=0 h=0.2 f=0I=15 输入ai 是i+1 iNC(a,h,n,&r,f) ntf 否ntf= =0 是 否输出Wrong!Return code=ntf输出R=r 结束 图 2 牛顿-科特斯公式流程设计图用Matlab 编程计算，和真解比较误差，结果看下三表，（程序附后。） 表 1 复化梯形公式，y(t)- y n (t)表 2 复化Simpson 公式 ，y(t)- y n (t)表 3 NewtonCotes公式 ，y(t)- y n (t) 比较表1 和表2 和表 3可以看出，用复化Simpson 公式离散比复化梯形公式离散得到结果的精度明显高出很多，这和复化Simpson 公式的代数精度比复化梯形公式高是一致的。 对于 NewtonCotes公式来说随着n 值的增大 其结果精确度也就越高，但不是越大越精确，一般来说n=4最为适合。4、结论 一般情况下可以采用复化梯形公式，复化Simpson 公式和 NewtonCotes公式离散积分方程求得近似解，采用复化梯形公式和复化Simpson公式的结果可以进一步外推提高精度。在不要求很高精度（1.0e-10）的时候，采用复化梯形公式相对较为简单，尤其是当被积函数是周期函数时，由于梯形公式对函数的光滑性要求比较低所以梯形公式是较理想的方法。复化Simpson 公式，具有较高的代数精度，且程序设计也比较简单，因此是使用比较广泛的一种方法，但取点相对较多。NewtonCotes公式具有较高的代数精度，但由于NewtonCotes积分是通过拉格朗日多项式插值变化而来的，我们都知道高次多项式插值会出现Runge振荡现象，因此会导致高阶的NewtonCotes公式不稳定。由于计算机舍入误差的存在，NewtonCotes公式的节点不能无限制增加，从而限制了精度的提高，只能采用复化梯形公式或复化Simpson 公式或它们的外推，通过增加节点个数提高精度使达到所需要的精度。 5、结束语 通过本次设计，我更进一步地了解了复化梯形公式、复化Simpson公式、牛顿科特斯公式在求解定值函数积分及一些实际问题上的应用，尤其是对实例进行程序设计方法上有了更深的理解。在书本上学习有关进程的同步的问题只是理论上的说法，并没有一个实际的东西，至少应是看得见的，对其含义也只是处于表面上的理解。这次实验我才体会到进程同步的真正涵义。 在本论文完成之际，我要向所有帮助过我的老师、同学表示衷心的感谢！我要特别感谢我的指导老师*老师的热情关怀和悉心指导。在我撰写论文的过程中，*老师他给予了殷切的指导，提出了许多宝贵的意见。无论是在论文的选题、构思和资料的收集方面，还是在论文的研究方法以及成文定稿方面，我都得到了*老师悉心细致的教诲和无私的帮助，特别是他广博的学识、严谨的治学精神和一丝不苟的工作作风使我受益匪浅，在此表示真诚地感谢和深深的谢意。 参考文献：1 王世儒、王金金、冯有前：计算方法（第二版），西安电子科技大学出版社。2 王正林、刘 明：精通Matlab7.， 北京:电子工业出版社, 2006.7。 3 李庆扬：数值分析（第四版），华中科技大学出版社。附录，程序1.1 复化梯形公式和复化Simpson公式MATLAB实现#includeiostream.h#includemath.hclass tixing /梯形类public: tixing();/初始化 void SetValue(double LowLimit1,double UpLimit1);/设置积分上下限值 double f(double x);/函数f(x) virtual double Result(void);/积分结果,虚函数 virtual void wucha(void);/误差范围,虚函数protected: double LowLimit;/上限 double UpLimit;/下限;tixing:tixing() LowLimit=0.0;/初始值 UpLimit=0.0;void tixing:SetValue(double LowLimit1,double UpLimit1)/设置函数 LowLimit=LowLimit1; UpLimit=UpLimit1;double tixing:f(double x) return exp(-x);/本例f(x)为e(-x)double tixing:Result(void) return (UpLimit-LowLimit)/2*(f(LowLimit)+f(UpLimit);/梯形公式void tixing:wucha(void) double x; x=-pow(UpLimit-LowLimit),3)/12*(-f(LowLimit);/梯形误差公式 cout梯形公式误差范围=xendl;class xinpusheng:public tixing/辛浦生类,从梯形继承public: virtual double Result(void);/对结果重新定义(因为公式不同) virtual void wucha(void);/同上,对误差重新定义;double xinpusheng:Result(void)/结果 return (UpLimit-LowLimit)*(f(LowLimit)+4*f(LowLimit+UpLimit)/2.0)+f(UpLimit)/6.0;void xinpusheng:wucha(void)/误差 double x; x=-pow(UpLimit-LowLimit),5)/2880*(-f(LowLimit); coutwucha:; if (x0) x=-x; cout辛浦生公式误差范围=:xendl;int main(void) double LowLimit1=0.0,UpLimit1=0.0; tixing tixing1;/梯形对象 xinpusheng xinpusheng1;/辛浦生对象 coutf(x)=e(-x)endl; coutUpLimit1; coutLowLimit1; tixing1.SetValue(LowLimit1,UpLimit1);/梯形赋值 xinpusheng1.SetValue(LowLimit1,UpLimit1);/辛浦生赋值 cout积分结果: 梯形tixing1.Result()endl; tixing1.wucha(); cout积分结果: 辛浦生xinpusheng1.Result()endl; xinpusheng1.wucha(); return 0;/计算方法,数值分析中的梯形积分公式,辛浦生积分公式,简称梯形公式,辛浦生公式求函数f(x)=e(-x)在积分上下限/为10的范围内的积分,并算出了各自算法的误差./下面是对上面的代码简化#includeiostream.h#includemath.hdouble f(double x) return exp(-x);/本例f(x)为e(-x)double tixingResult(double UpLimit,double LowLimit) return (UpLimit-LowLimit)/2*(f(LowLimit)+f(UpLimit);/梯形公式void tixingwucha(double UpLimit,double LowLimit) double x; x=-pow(UpLimit-LowLimit),3)/12*(-f(LowLimit);/梯形误差公式 cout梯形公式误差范围=xendl;double xinpushengResult(double UpLimit,double LowLimit)/结果 return (UpLimit-LowLimit)*(f(LowLimit)+4*f(LowLimit+UpLimit)/2.0)+f(UpLimit)/6.0;void xinpushengwucha(double UpLimit,double LowLimit)/误差 double x; x=-pow(UpLimit-LowLimit),5)/2880*(-f(LowLimit); coutwucha:; if (x0) x=-x; cout辛浦生公式误差范围=:xendl;int main(void) double LowLimit,UpLimit; coutf(x)=e(-x)endl; coutUpLimit; coutLowLimit; cout积分结果: 梯形tixingResult(UpLimit,LowLimit)endl; tixingwucha(UpLimit,LowLimit); cout积分结