数值分析课程设计

《数值分析》课程设计第六组学号10801010106姓名侯杰学号10801010121姓名王豪学号10801010206姓名付渝翔学号10801010221姓名孙珊珊2011年6月数学与统计学院专业课程设计成绩评定书设计题目：题目6信息与计算科学专业108010101108010102班0621学号指导教师刘瑞华牛普许安见指导教师评语成绩：指导教师时间：答辩小组意见设计成绩：答辩组长：审定系主任：摘要:用复化梯形公式、复化Simpson公式和复化Cotes公式以及用Romberg积分法计算积分，精度计算。关键字梯形公式SimpsonCotes题目（六）6、用熟悉的计算机语言编程上机完成（1）用复化梯形公式、复化Simpson公式和复化Cotes公式计算积分10sindxxx，自己设置不同精度要求，对结果进行比较分析。（2）用Romberg积分法计算积分10sindxxx，自己设置不同精度要求，对结果进行比较分析。（3）记xxxfsin)(，在上面的计算中)(xf只取4位有效数字或7位有效数字，计算结果又有什么不同。（4）上面计算精度可达8-20位有效数字吗？若可以，请说明实现过程，并举例。一、复化求积公式的基本思想1.柯特斯系数与辛普森公式设将积分区间ba,划分为n等份，步长hnab，选取等距节点khaxk构造出的插值型求积公式In=)()(0)(knknkxfcab(1.1)称为牛顿-柯特斯公式，式中)(nkc称为柯特斯系数，引进变换，则有khax)(nkc)!(!100knnkdtjkjtabhnkjjnkn)1(nnkjjdtjt00)((1.2)由于是多项式的积分，柯特斯系数的计算不会遇到实质性的困难，则当1n时，21)1(1)1(0cc,这时的求积公式就是我们所熟悉的梯形公式。当2n时，按（1.2）式，这时的柯特斯系数为)2(0c=20)2)(1(41dttt=61,)2(1c=-2064)2(21dttt,)2(2c=61)1(4120dttt,相应的求积公式时下列的辛普森公式：s)]()2(4)([6bfbafafab，（1.3）而n=4时的牛顿-柯特斯公式则特别称为柯特斯公式，而形式是c)(7)(32)(12)(32)(79043210xfxfxfxfxfab(1.4)这里h4ab,khax(4,3,2,1,0k)2.复合梯形公式将区间ba,划分为n等份，分点khaxk，nabh,nk,2,1,0,在每个子区间1,kkxx(1,2,1,0nk)上采用梯形公式（1.1），则得badxxfI)(=10nk1)(xkxkdxxf=102nkh)()()(1fRxfxfnkk(2.1)记102nknhT)()(1kkxfxf=2h)]()(2)([11bfxfafknk(2.2)称为复合梯形公式，其余项可由（1.10）式得10)(nknnTIfR)](12[''3kfh,),(1kkkxx,由于],[)(2baCxf，且)(10minkfnk101nkn)(10max)(kkfnkf，所以),(ba使)(1)(10knkfnf于是复合梯形公式的余项为)(12)(2fhabfRn（2.3）3.复合辛普森求积公式将区间ba,分为n等份，在每个子区间1kkxx上采用辛普森公式，若hxxkk212/1,则得I=badxxf)(=101)(nkxkxkdxxf=106nkh)()(2)(4)(12/1fRxfxfxfnkkk（3.4）记ns=106nkh112/110)()(2)(4)(nkkknkbfxfxfaf=6h104)([nkaf)(2/1kxf+210nk)()(bfxfk]（3.5）称为复合辛普森求积公式，其余项由（2.5）得).,(),()2(180)(110)4(4kkkknknnxxfhhSIfR于是当baCxf,)(4时，与复合梯形公式相似有).,(),()2(180)(44bafhabSIfRnn（3.6）由（3.6）式看出，误差阶为h4，收敛性是显然的，实际上，只要baCxf,)(则可得到收敛性，即nlimbandxxfS)(此外，由于Sn中的求积系数均为正数，故知复合辛普森公式计算稳定二、算法实现1、对积分10sinxIdxx,使其精度达到610解：先建立三种复化公式的函数文件,它们分别为复化梯形公式trap.m、复化Simpson公式为simpson.m、Cotes公式为cotes.m,三个函数的源程序如下:(1)复化梯形公式trap.mfunctionT=trap(f,a,b,n)%trap.m%复化梯形公式求积分值%f为积分函数%[a,b]为积分区间%n是等分区间份数h=(b-a)/n;%步长T=0;fork=1:(n-1)x0=a+h*k;T=T+limit(f,x0);endT=h*(limit(f,a)+limit(f,b))/2+h*T;T=double(T);(2)复化Simpson公式simpson.m:functionS=simpson(f,a,b,n)%simpson.m%Simpson公式求积分值%f为积分函数%[a,b]为积分区间%n是等分区间份数h=(b-a)/(2*n);%步长s1=0;s2=0;fork=1:nx0=a+h*(2*k-1);s1=s1+limit(f,x0);endfork=1:(n-1)x0=a+h*2*k;s2=s2+limit(f,x0);endS=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3;S=double(S);(3)复化Cotes公式cotes.m:functionC=cote(f,a,b,n)%cote.m%复化cotes公式求积分值%f为积分函数%[a,b]为积分区间%n是等分区间份数h=(b-a)/n;%步长C=0;fori=1:(n-1)x0=a+i*h;C=C+14*limit(f,x0);endfork=0:(n-1)x0=a+h*k;s=32*limit(f,x0+h*1/4)+12*limit(f,x0+h*1/2)+32*limit(f,x0+h*3/4);C=C+s;endC=C+7*(limit(f,a)+limit(f,b));C=C*h/90;C=double(C);输入命令其结果显示如下：symsx;f=sym('sin(x)/x');simpson(f,0,1,20)ans=0.94608307075153simpson(f,0,1,1)ans=0.94614588227359simpson(f,0,1,10)ans=0.94608307651773trap(f,0,1,20)ans=0.94602032535503trap(f,0,1,1)ans=0.92073549240395trap(f,0,1,10)ans=0.94583207186691cote(f,0,1,20)ans=0.94608307036718cote(f,0,1,1)ans=0.94608300406367cote(f,0,1,10)ans=0.946083070367122、用Romberg积分法计算10sinxIdxx,精度610.解：编写Romberg积分法的函数M文件,源程序如下(romberg.m):function[I,T]=romberg(f,a,b,n,Eps)%Romberg积分计算%f为积分函数%[a,b]为积分区间%n+1是T数表的列数目%Eps为迭代精度%返回值中I为积分结果,T是积分表ifnargin5Eps=1E-6;endm=1;h=(b-a);err=1;j=0;T=zeros(4,4);T(1,1)=h*(1+feval(f,b))/2;while((errEps)&(jn))|(j4)j=j+1;h=h/2;s=0;forp=1:mx0=a+h*(2*p-1);s=s+feval(f,x0);endT(j+1,1)=T(j,1)/2+h*s;m=2*m;fork=1:jT(j+1,k+1)=T(j+1,k)+(T(j+1,k)-T(j,k))/(4^k-1);enderr=abs(T(j,j)-T(j+1,k+1));endI=T(j+1,j+1);ifnargout==1T=[];end将上述源程序另存为romberg.m后,进入计算:symsx;%创建符号变量f=sym('sin(x)/x')%符号函数f=sin(x)/x[I,T]=romberg(f,0,1,3,1E-6)%积分计算I=0.9461T=0.920700000.93980.94610000.94450.94610.9461000.94570.94610.94610.946100.94600.94610.94610.94610.9461其中T为Romberg积分表,由输出结果可知计算结果为10sin0.9461xIdxx.记xxxfsin)(，在上面的计算中)(xf只取4位有效数字或7位有效数字，计算结果又有什么不同。symsx;f=sym('sin(x)/x');trap(f,0,1,5)ans=0.9451取四位有效数字得：vpa(trap(f,0,1,5),4)ans=0.9451取7为有效数字得：vpa(trap(f,0,1,6),7)ans=0.9453857上面计算精度可达8-20位有效数字吗？若可以，请说明实现过程，并举例。可以实现举例如下：formatlongansans=.9453857trap(f,0,1,6)ans=0.94538573076686vpa(ans,2)ans=0.95实验总结：通过本次数值分析课程设计分析实验，教我们学会了如何使用MATLAB实现复化梯形、辛普森、复化cotes以及龙贝格公式计算积分。参考文献【1】《数值计算方法》合肥工业大学数学与信息科学系编合肥工业大学【2】黄友谦，程诗杰，陈浙鹏，《数值试验》，北京：高等教育出版社，1989【3】蔡大用，《数值分析与实验学习指导》，北京：清华大学出版社与施普林格出版社，2001【4】肖筱南，《值计算方法与上机实习指导》，北京：北京大学出版社，2004【5】刘萍.数值计算方法北京：人民邮电出版社，2002.2【6】大用.数值分析与实验学习指导北京：清华大学出版社，2001