辛普森求积公式.doc

摘要在工程实验及研究中，实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.可以说，曲线拟合模型与我们的生活生产密切相关.本课题着重介绍曲线拟合模型及其应用，其中包括它的基本思想、模型的建立、以及具体应用.为了更好的了解曲线拟合模型，可以将它分为线性与非线性模型，在模型建立的基础上我们可以用最小二乘法来解决一些我们日常所应用的问题.关键词 曲线拟合；线性与非线性模型；最小二乘发 目录引言1第一章 曲线拟合21.1 基本思想及基本概念21.1.1 方法思想21.1.2几个基本概念21.2辛普森算法基本定义及其应用41.2.1辛普森求积公式的定义41.2.2辛普森求积公式的几何意义51.2.3辛普森求积公式的代数精度及其余项51.2.4辛普森公式的应用6第二章 辛普森求积公式的拓展及其应用72.1 复化辛普森求积公式72.1.1问题的提出72.1.2复化辛普森公式及其分析72.1.3复化辛普森公式计算流程图82.1.4复化辛普森公式的应用92.2 变步长辛普森求积公式102.2.1变步长辛普森求积公式的导出过程102.2.2变步长辛普森求积公式的加速过程122.2.3变步长辛普森求积公式的算法流程图132.2.4变步长辛普森公式算法程序代码142.2.5变步长辛普森求积公式的应用142.2.6小结142.2.7数值求积公式在实际工程中的应用14参考文献16附录A17附录B18附录C21引言辛普森是英国数学家.1710年8月20日生于波士沃希；1761年5月14日卒于波士沃希.在定积分近似计算中，以他的姓来命名的“辛普森公式”，虽早在他之前牛顿的学生柯特斯（Cotes）和斯特林就已经得出了（包括一些更高阶的近似公式），但真正广泛地为人所知并加以应用，则是1743年辛普森重新发现之后的事了.辛普森的工作使牛顿的微积分学说得到了进一步完善.在我们的日常生活中计算积分与我们的生活生产密切相关.所以掌握数值积分方法是学生储备知识能量的武器.数值积分的一个基本的计算策略，用易于积分的简单函数来逼近曲线.简单曲线下面的面积近似等于下面的面积.如果涉及初等函数的积分找不到其他由初等函数构成的解析表达式，或者只在一些离散的点上知道函数的值，在多数情况下，被积函数的原函数很难用初等函数表达出来， 因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.那么就必须对定积分进行数值逼近.数值积分实现是将整个闭区间划分为个小段，在每个小段上对进行低阶分段多项式逼近.对每个小段上的逼近多项式积分时，就得到基本公式.基本公式只涉及足够多的对来定义分段多项式的某一段，将此公式应用到个小段并把结果相加得到复合公式，或称为扩展公式. 在一个小段中节点的位置和数目决定了基本公式的很多重要特性.当节点均匀分布时，所有的积分公式称为牛顿柯特斯公式.例如，梯形、辛普森、柯特斯求积公式等.经典辛普森求积公式来源于Lagrange插值多项式的应用，它的代数精度高达3阶，其形变后的代数精度高达4阶，且二者都具有良好的稳定性与收敛性，从而提高了计算效率及准确度，是定积分近似计算常使用的方法，一直是理工科大学生必修的内容. 下面将给出具体辛普森求积公式的具体思想以及其算法程序设计并给出将其拓展后在实际工程问题中的应用.第一章 辛普森求积公式的理论实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系.依据人们所熟知的微积分基本定理，对积分只要找到被积函数的原函数,便有下列牛顿-莱布尼茨公式：，但实际计算往往遇到一些困难，如: 1)的原函数不能用初等函数表示，故不能用牛顿-莱布尼茨公式计算.2) 虽然找到了 的原函数, 但因表达式过于复杂而不便应用牛顿-莱布尼茨公式.3) 在许多实际问题中是以列表函数的形式给出, 即仅仅知道其在一些节点处的函数值， 牛顿-莱布尼茨公式也不能直接运用，因此有必要研究积分的数值计算问题，数值积分是解决上述困难的一种有效方法.1.1基本思想及基本概念1.1.1 方法思想由定积分中值定理：可知: 积分可以通过被积函数在处的值得到. 由于积分中值定理仅仅告诉我们在一定条件下是存在的, 但并没有给出确定的方法. 一个很自然的想法就是利用被积函数在节点处函数值的加权平均来替代(近似), 按此思想有 （1-1）这就是数值求积的思想（有效地解决了本章开始提出的问题），权因子和节点 的不同确定方法就对应不同的数值求积公式.1.1.2 几个基本概念定义1.1 称形如(1-1)式的求积公式为机械求积公式,其中仅节点的选择与无关，称为求积节点，（）称为求积系数.定义1.2 如果某个求积公式对于次数不超过的多项式均能准确地成立,而对于次多项式就不准确成立, 则称该求积公式具有次代数精度（或代数精确度）.注1.1 a) 越大近似程度越高，标志着使函数准确成立的“个数”越多，但代数精度不是唯一衡量标准. b) 若机械求积公式的代数精度，则有.c) 若机械求积公式的代数精度为，即当时，由（1.1）式可得，对任意次数不超过的次多项式有.d) 代精度的高低, 从侧面反映求积公式的精度高低.定义1.3 称求积公式为插值型求积公式，式中求积系数通过插值基函数积分求得，即 （1-2）定理1.1 插值型求积公式的代数精度至少为次.定义1.4 若节点将被积区间等分成等分, 即则相应的插值求积公式称为Newton-Cotes (牛顿-柯特斯)求积公式. 即等距节点情形下的插值求积公式称为牛顿-柯特斯公式, 相应的求积系数称为Cotes系数.常见的几个简单求积公式( Newton-Cotes公式)，如表1-1所示：表1-1 几种简单N-C求积公式总结表名称形式梯形求积公式辛普森求积公式柯特斯求积公式其中注1.2 a）时，N-C公式出现数值不稳定. b）为偶数时，N-C公式的代数精度至少为次，为奇数时，N-C公式的代数精度至少为次.定义1.5 截断误差: 由 （1-3）当时可得梯形求积公式的截断误差类似的，可得当，时的截断误差注1.3 从截断误差公式可知，当区间长度较大时，求积公式误差较大.1.2辛普森算法基本定义及其应用1.2.1 辛普森求积公式的定义设计积分区间划分为等份，步长，选取等距节点构造出的插值型求积公式为牛顿柯特斯（Newton-Cotes）公式，式中称为柯特斯系数.根据插值型求积公式系数（1-2），引进变换，则有当时，由上式有 则相应的求积公式是辛普森求积公式： （1-4）1.2.2辛普森求积公式的几何意义辛普森公式的几何意义就是用通过A,B,C三点的抛物线代替所得曲边梯形面积，如图1.1所示.Y=f（x）图1.1 辛普森求积公式的几何意义图 yx O 0a bABC1.2.3辛普森求积公式的代数精度及其余项由N-C公式的特点知，当为偶数时N-C公式的代数精度至少为次，由于Simpson求积公式为时的N-C公式，故它的代数精度至少为3次，即将代入Simpson公式（1-4）左边右边左边由此可知使得Simpson求积公式不准确成立，所以即Simpson公式代数精度为3次由N-C公式的余项公式（1-3）知，当时可得辛普森求积公式的截断误差 （1-5）1.2.4辛普森公式的应用例1.1 用辛普森求积公式计算积分.由积分形式可知 用辛普森公式计算有下式其中.计算流程图图1.2 例1.1流程图开始定义函数f（x）输入n，a，b的值计算h=（b-a）/n调用函数f（x），计算s的值输出s的值结束C语言程序代码及其运算结果详见附录A分析附录A可知第二章 辛普森求积公式的拓展及其应用为了提高精度，通常在实际应用中往往采用将积分区间划分成若干个小区间，在各小区间上采用低次的求积公式，如：梯形公式或辛普森公式，然后再利用积分的可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复化求积公式，本章重点介绍复化辛普森求积公式.2.1 复化辛普森求积公式2.1.1问题的提出由截断误差可知，当区间长度较大时，Newton-Cotes求积公式的误差较大. 为构造更高精度的数值积分公式，可以采用分段低次多项式替代整体高次多项式，为此，利用积分关于区间具有可加性，将区间上的积分，分成若干小区间上的积分，以此来减少积分区间长度引起的误差.这就引用了复化求积公式. 其基本思想是：先把积分区间分成一些长度较小的子区间，在每个子区间上使用低阶的牛顿-柯特斯公式，即利用并把小区间上的积分用前面的方法近似求得，由此即可得到相应的复化求积公式. 最常用的是复化梯形公式和复化辛普森公式，下面学习辛普森求积公式.2.1.2复化辛普森公式及其分析定义2.1 将小区间上的积分分别用辛普森公式计算，即可得到复化辛普森公式其中.另一种定义形式为：用分段二次插值函数代替，记在第段的两个小区间上，用三个结点作二次插值函数，然后积分，求段之和可得整个区间上的近似积分称该求积公式为复化辛普森求积公式（抛物线公式）.定理2.1 若则复化辛普森公式的截断误差为且.注2.1 从误差公式可以看出当时，比的精度一般要高，但他们的计算量几乎一样.注2.2 属于机械型求积公式，但不属于插值型、也不属于N-C求积公式.的代数精度为4次，具有稳定性和收敛性即（或）.2.1.3复化辛普森公式计算流程图为了减少计算工作量，优化程序设计，将复化辛普森公式改写为则于此相对应的辛普森流程图为：图2.1 复化辛普森算法流程图开始输入A,B,N H=（B-A）/（2*N）S=0.5*（F（A）- F（B），调用函数FS=S+2*FA+（2*I-1*H）+（F（A+2*I*H），调用函数FS=（B-A）/（3*N）S输出S结束I=1,N定义函数F2.1.4 复化辛普森公式的应用例2.1 用复化辛普森公式计算正弦积分的近似值.分析该积分可知 ，则为步长C语言程序代码及其运算结果详见附录B由此可知例2.2 用复化辛普森公式计算定积分.分析该积分可知，则为步长C语言程序代码及其运算结果详见附录B.由此可知在利用插值型求积公式求积分时，为了提高精度有两种途径.一是提高积分区间上的插值多项式的阶数，从而也就提高了求积的阶数.但是，由于插值多项式的阶数越高，其逼近性质未必好（即精度未必能提高），因此，牛顿-柯特斯公式的阶数越高，其积分精度也未必越高，工程上一般只作到六阶牛顿-柯特斯公式（即龙贝格公式）为止.二是采用复化公式，尽量减小每个求积小区间的长度.在实际应用时，往往将两种方法混合使用，以便提高求积的精度.2.2 变步长辛普森求积公式在数值积分中，精度是一个很重要的问题，如果误差太大，就没有实际意义.为了提高精度，通过需要在复化求积公式中尽量减少各细分小区间的长度，即减少步长.显然，如果步长取得太大，则精度就难以保证.但是，如果步长取得太小，则计算工作量就随之增大，并且，由于项数增加，其误差积累也就增大.因此，在采用复化公式求积时，关键的问题是合理地选择步长（即合理选择对整个积分区间的细分数），以便既能满足精度要求，又不至于引起过多的误差积累和过大的计算工作量.在实际计算过程中，通常采用变步长的求积法.2.2.1变步长辛普森求积公式的导出过程变步长辛普森求积公式是建立在变步长梯形公式的基础上，同时它又是龙贝格算法导出的中间过程，我们知道, 若被积函数具有一定的光滑性, 则增加节点可以降低复化求积公式的截断误差.这里需要解决的问题是增加节点后的复化求积方法能否充分利用已有的计算工作量. 譬如: 若将作为的近似精度不够, 需减少步长(增加节点数)计算相应的来近似, 当然我们想要充分利用已经求得的.为此, 设区间等分后, 利用复化梯形公式已经求得这一结果, 为了得到精度更高的数值结果, 我们将原有的步长折半, 即把区间分为等分, 然后应用复化梯形公式求得.下面将会看到这样既提高了精度, 又能充分利用已经求得的.事实上, 我们可以建立与的下述递推关系. 设则其中即，注2.3 由上述公式可知在的基础上计算只需调用次函数即可，最大限度地节省了的计算量.加速公式的导出：由前面的误差分析，我们可以得到复化梯形公式的截断误差为，即类似根据复化梯形公式的截断误差为，有两式相比可得， 其中即 （2-1）注2.4 公式（2-1）说明的误差可以近似地由与表现， 这样就给出了复化梯形公式估计误差的事后估计法. 由公式（2-1）还可以得到校正公式（加速公式）数值实验结果表明，在一定条件下，上式计算出来的值比原来的好得多，上述公式称为梯形公式的加速公式.梯形求积公式的实质：假设已知，则即上式表明与通过上面公式处理后，可得精度更高的.即复化辛普森公式，这也是加速的实质.2.2.2变步长辛普森求积公式的加速过程类似梯形加速公式的推导，由的截断误差公式（1-5）可得即注2.5 上述两个公式分别称为复化辛普森公式估计误差的事后估计公式及复化辛普森公式的加速公式. 类似地可以证明：在求得，的基础上，可以进一步加速得：龙贝格公式2.2.3变步长辛普森求积公式的算法流程图开始N=1，H=B-AIP=F(A)+F(B)FIC=0,X=A-H/2K=1,NX=X+HIC=IC+F(x)FI2=(4*IC+IP)*H/6N=1|I2-I|ESPI1=I2,IP=IP+2*ICN=N+NH=0.5*HYNNI=2I输出结束图2.2 变步长辛普森算法流程图2.2.4变步长辛普森公式算法程序代码详见附录C2.2.5变步长辛普森求积公式的应用例2.3 用变步长辛普森求积公式计算定积分取.C语言程序代码及其运算结果详见附录C.分析结果可知 2.2.6小结通过分析例1.1、2.2、2.3有下表2-1表2-1 三种算法比较算法名称代数精度积分形式计算结果余项辛普森求积30.111765复化辛普森求积40.11157变步长辛普森求积0.111572由表2-1可以得出用变步长辛普森求积公式求得的结果偏离准确值的程度最小，即其计算结果最接近准确值，其次是复化辛普森求积方法，辛普森求积方法较前述两种方法误差较大.但三种算法均具有良好的稳定性与收敛性，从而提高了计算效率及准确度在工程技术中有较为广泛的应用.2.2.7 数值求积公式在实际工程中的应用例 2.4人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球表面439Km，远地点距地球表面2384Km，地球半径为6371Km.求：该问题的轨道长度.模型 1），分别是长半轴和短半轴；2）焦距为；3）地球半径为；近地点和远地点与地球表面的距离分别是和.图2.3 卫星轨道的示意图 椭圆的参数方程为，弧长的公式椭圆长度利用数值积分公式计算卫星轨道长度问题，编写Mathematica语句如下：结果为:48707.4参考文献1李庆扬 王能超 易大义.数值分析基础M.北京:清华大学出版,2008年. 2徐士良.数值方法与计算机实现M.北京:清华大学出版社，2006年.3封建湖 聂玉峰 王振海. 数值分析M.西安:西北工业大学出版社，2003年.4朝伦巴根 贾德彬. 数值计算方法M.北京:中国水利水电出版社，2007年.5王立秋 魏焕彩 周学圣.工程数值分析题解M.山东：山东大学出版社，2004年.6张韶华 奚梅成 陈效群 .数值计算方法与算法.北京：科学出版社，2006年.附录A本附录介绍例1.1 用辛普森求积公式计算过程中的变量说明、C语言程序代码、以及运行结果.表 变量说明表变量名变量说明a积分上限b积分下限n区间等分后的个数h步长s结果该积分用Simpson公式计算的C语言程序代码为：#include#includefloat f (float x) float t; t=x/(4+pow(x,2); return t;main()float n,a,b;float h,s;printf(putn=);scanf(%f,&n);printf(puta=);scanf(%f,&a);printf(putb=);scanf(%f,&b);printf(a=%fb=%f,n=%fn,a,b,n);h=(b-a)/n;printf(h=%fn,h);s=(h/3)*(f(a)+4*f(a+h)+f(b);printf(s=%fn,s);其结果为：图 例1.1计算结果附录B本附录介绍例2.1和例2. 2 用复化辛普森求积公式计算过程中的变量说明、C语言程序代码、以及运行结果.表 变量说明变量名变量说明a积分上限b积分下限n区间等分后的个数h步长s结果s1中间变量s2中间变量m中间变量i中间变量j中间变量l中间变量d中间变量例2.1程序代码：#include#includefloat f (float x) float t; t=sin(x)/x; return t;main()int n=8,i,j;float m,s,s1=0,s2=0,h,a=0,b=1,c,l,d;h=(b-a)/4;c=h/6;printf(a=%f,a);printf(b=%f,b);printf(h=%f,h);printf(f(0.000000)=1.000000n);l=1;/f(a)=1.000000;l=f(a);d=f(b);for(m=0.125;m=1;m=m+0.125)printf(f(%f)=%fn,m,f(m);for(i=1;i8;i=i+2)s1=s1+4*f(a+i*h/2);for(j=2;j8;j=j+2)s2=s2+2*f(a+j*h/2);s=c*(l+d+s1+s2);printf(s1=%fn,s1);printf(s2=%fn,s2);printf(s=%6.5fn,s);例2.1程序运行结果为：图 例2.1计算结果例2.2变量说明同例2.1例2.2程序代码：#include#includefloat f (float x) float t; t=t=x/(4+pow(x,2); return t;main()int n=8,i,j;float m,s,s1=0,s2=0,h,a=0,b=1,c,l,d;h=(b-a)/4;c=h/6;printf(a=%f,a);printf(b=%f,b);printf(h=%f,h);printf(f(0.000000)=1.000000n);l=f(a);d=f(b);for(m=0.125;m=1;m=m+0.125)printf(f(%f)=%fn,m,f(m);for(i=1;i8;i=i+2)s1=s1+4*f(a+i*h/2);for(j=2;jadouble eps; /积分精度要求double (*f)(); /指向计算积分函数值 double simp(); /函数返回一个双精度实型积分值/计算被积函数f(x)值的函数形式为double f(x);double x; double y; y=被积函数f(x)的表达式； return(y); #include#includedouble simp(a,b,esp,f)double a,b,esp,(*f)(); int n,k; double h,t1,t2,s1,s2,ep,p,x; n=1; h=