内江师范学院数值分析实验三

1、内江师范学院数 值 分 析实 验 报 告 册编制 张莉 审定 牟廉明专业： 班级： 级 班学号： 姓名： 数学与信息科学学院2013年9月说 明一、学生在做实验之前必须要准备实验，主要包括预习与本次实验相关的理论知识，熟练与本次实验相关的软件操作，收集整理相关的实验参考资料，要求学生在做实验时能带上充足的参考资料；若准备不充分，则学生不得参加本次实验，不得书写实验报告；二、要求学生要认真做实验，主要是指不得迟到、早退和旷课，在做实验过程中要严格遵守实验室规章制度，认真完成实验内容，极积主动地向实验教师提问等； 三、各个实验按照学生水平分别设置了A、B、C、D四个等级，其中对应的难度系数为1、0

2、.8、0.7、0.6，也可根据实际完成情况制定相应地的难度系数，但总体保证难度排序为A级难度最大，B级次之，C级较易，D级最简单。四、学生可以根据自己对各个实验涉及到的知识点掌握的程度自由选取A、B、C、D等级的实验题目。五、学生要认真工整地书写实验报告，实验报告的内容要紧扣实验的要求和目的，不得抄袭他人的实验报告；四、根据实验准备、实验态度、实验报告的书写、实验报告的内容进行综合评定，并给出实验成绩评定分。实验名称： 实验三 数值积分 指导教师： 吴开腾 张莉 实验时数： 2 实验设备：安装了Matlab、C+、VF软件的计算机 实验日期：2013年 10 月 29 日 实验地点： 第五教学

3、楼北902 实验目的：1. 掌握数值积分的基本思想和基本步骤。2. 理解各类数值积分方法的优缺点，并能自行编程求解。3. 比较各类数值积分的代数精度，体会复化、变步长及其加速的思想和实现步骤。实验准备：1. 在开始本实验之前，请回顾教科书的相关内容；2. 需要一台准备安装Windows XP Professional操作系统和装有数学软件的计算机。实验内容及要求A题 计算积分。B题 用Romberg法求函数积分，精度为。C题 用复化梯形公式计算积分并分析剖分区间数对误差的影响，取，积分精确值。D题 计算。说明：实验过程应包括对问题的简要分析、求解方法、求解步骤、程序及其必要的图表等内容。实验过

4、程：A题的实验过程1、 实验中问题的重述计算积分2、 实验分析由于积分为是从负无穷到正无穷积分，用常规的求积分的int函数求解时，会显示Warning:Explicit integral could not be found。说明函数y=exp(sin(x)-x.2/100)不是初等原函数，可用数值积分进行分析。为方便计算，将积分区间分为两个区间，以0为分界点，分为和，则有：所以在两个区间分别计算，然后将结果相加。3、 实验求解过程首先利用辛普生法，从结果得到：在积分区间积分上限增加，积分值趋于12.3599；在积分区间积分下限减小，积分值趋于10.0804；然后利用牛顿柯特斯法，观察结果得到

5、，的积分区间，积分上限增加，积分值趋于12.3599；在的积分区间，积分下限减小，积分值趋于10.0804；4、 实验求解结果从上述两种方法的实验结果表明，辛普生法和牛顿柯特斯法在的积分区间，积分上限增加，积分值均趋于12.3599；在的积分区间，积分下限减小，积分值均趋于10.0804，故在两个积分区间都是可积的，所以故I=10.0804+12.3599=22.4403B题的实验过程1、 实验中问题的重述用Romberg法求函数积分，精度为。2、 实验求解程序M文件function quad,R = Romberg(f,a,b,eps)%f表示被积函数%a，b表示积分区间的端点%eps表示精

6、度%quad是用Romberg加速算法求得的积分值%R为Romberg表%err表示误差的估计h=b-a;R(1,1)=h*(feval(f,a)+feval(f,b)/2;M=1;J=0;err=1;while err>eps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; M=2*M; for k=1:J R(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J)-R(J+1,J+1);e

7、ndquad=R(J+1,J+1)主程序f=(x)(sin(x)./x);quad,R=Romberg(f,eps,1,0.5*10(-6)3、 实验求解结果用Romberg加速算法求得的积分值quad = 0.9461Romberg表C题的实验过程1、 实验中问题的重述用复化梯形公式计算积分并分析剖分区间数对误差的影响，取，积分精确值。2、 实验分析在许多实际问题中，常常需要计算定积分的值。根据微积分学基本定理，若被积函数f(x)在区间a,b上连续，只要能找到f(x)的一个原函数F(x)，便可利用牛顿-莱布尼兹公式求得积分值。但是在实际使用中，往往遇到如下困难，而不能使用牛顿-莱布尼兹公式。

8、a) 找不到用初等函数表示的原函数b) 虽然找到了原函数，但因表达式过于复杂而不便计算c) f(x)是由测量或计算得到的表格函数由于以上种种困难，有必要研究积分的数值计算问题。利用插值多项式则积分转化为，显然易算。称为插值型求积公式。最简单的插值型求积公式是梯形公式和Simpson公式，。当求积结点提供较多，可以分段使用少结点的梯形公式和Simpson公式，并称为复化梯形公式、复化Simpson公式。如步长未知，可以通过误差限的控制用区间逐次分半的策略自动选取步长的方法称自适应算法。梯形递推公式给出了区间分半前后的递推关系。由梯形递推公式求得梯形序列，相邻序列值作线性组合得Simpson序列,

9、 Simpson序列作线性组合得柯特斯序列, 柯特斯序列作线性组合的龙贝格序列。若|R2-R1|<e,则输出R2;否则依此类推。如此加工数据的过程叫龙贝格算法。3、 实验结果利用龙贝格算法，得到当分别取时，产生的误差见下表： 表1 误差结果表248误差-0.013239352830249-0.003315574025695-0.000828929586313163264误差-0.000207232961173-0.000051808249116-0.000012952062417128256512误差-0.000003238015606-0.000000809503902-0.00000

10、0202375974D题的实验过程1、 实验分析为了计算，设将积分区间划分为n等份，则一共有个等分点，则可从0到累加求和，即可导出如下递推算式：式中为二分前的步长，而。2、 实验求解过程已知被积函数，求积区间为。首先，准备初值： ，然后，根据上式推导公式：有最后，利用复化梯形公式求积分。3、 实验结果利用复化梯形公式求出积分为：=1.275067526872575e+002附 录A题的求解程序% 辛普生法clcclearfun=inline('exp(sin(x)-x.2/100)','x'); Isim=quad(fun,0,100)Isim=quad(fun

11、,0,1000) Isim=quad(fun,0,1000000)Isim=quad(fun,-100,0)Isim=quad(fun,-1000,0)Isim=quad(fun,-10000,0)% 牛顿柯特斯法clcclearfun=inline('exp(sin(x)-x.2/100)','x'); IL=quadl(fun,0,100) IL=quadl(fun,0,1000)IL=quadl(fun,0,10000)IL=quadl(fun,-10000,0)IL=quadl(fun,-1000,0)IL=quadl(fun,-100,0)B题的求解程

12、序function quad,R = Romberg(f,a,b,eps)%f表示被积函数%a，b表示积分区间的端点%eps表示精度%quad是用Romberg加速算法求得的积分值%R为Romberg表%err表示误差的估计h=b-a;R(1,1)=h*(feval(f,a)+feval(f,b)/2;M=1;J=0;err=1;while err>eps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; M=2*M; for k=1:J R(J+1,k+1)=R(

13、J+1,k)+(R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J)-R(J+1,J+1);endquad=R(J+1,J+1)主程序：f=(x)(sin(x)./x);quad,R=Romberg(f,eps,1,0.5*10(-6)C题的求解程序function T = Complexification(a,b,Y,n)% 复化梯形公式求解积分% a:下限% b：上限% n：点的个数c = 2*sum(Y(2:n-1);T = (b-a)/(2*n-2)*(Y(1)+c+Y(n);clcclear% 主程序a = 0; b = 1;n = 2 4 8 16 32 64 128 256 512;T = ;for i = 1:length(n) x = a:(b-a)/n(i):b; y = 4./(1+x.*x); tt = Complexification(a,b,y,length(x); T(i) = tt;endfprintf('分别产生的误差为： ')format longErr = (T - pi)./piD题的求解程序function T = Complexification(a,b,Y,n)% 复化梯形公式求解积分% a:下限% b：上限% n：点的个数c = 2*su