matlab定积分的近似计算 实验报告二

实验序号： 实验二 日期： 2015 年 05 月 09 日 班级 132132002 姓名 高馨 学号 1321320041 实验名称 定积分的近似计算 问题背景描述 定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分． 实验目的 本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算 实验原理与数学模型 MATLAB7.11.0 主要内容（要点） 1． 实现实验内容中的例子，即分别采用矩形法、梯形法、抛 1dx物线法计算，取n258，并比较三种方法的精确程 01x2度． 2dx2． 分别用梯形法与抛物线法，计算，取n120．并尝试 1x直接使用函数trapz()、quad()进行计算求解，比较结果的差异． 1dx3． 将的近似计算结果与Matlab中各命令的计算结果 01x2相比较，试猜测Matlab中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例子支持你的观点． 4． 学习fulu2sum.m的程序设计方法，尝试用函数 sum 改写附录1和附录3的程序，避免for 循环。 实验过程记录（含基本步骤、主要程序清单及异常情况记录等） 一、 实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算dx，取n258，并比较三 01x2 1种方法的精确程度． （一） 矩形法（中点法）程序： clear format long; a=0;b=1;n=258;h=(b-a)/n; C=0; for i=1:n; xi=(i-1)*h; xj=i*h; xz=(xi+xj)/2; C=C+h*(1/(1+xz*xz)); end disp(C) E=abs((C-pi/4)/(pi/4)); fprintf('The relative error between C and real-value is about: %d\\n',E) 答案： 0.785398476379441 The relative error between C and real-value is about: 3.985010e-007 （二）梯形法程序： clear format long; a=0;b=1;n=258;h=(b-a)/n; B=0; x1=0; y1=1/(1+x1*x1); B=h*y1/2; for i=1:n-1 xi=i*h; fxi=1/(1+xi*xi); B=B+h*fxi; end xn=1; yn=1/(1+xn*xn); B=B+h*yn/2; disp(B) E=abs((B-pi/4)/(pi/4)); fprintf('The relative error between B and real-value is about: %d\\n',E) 答案： 0.785397537433464 The relative error between B and real-value is about: 7.970021e-007 （三）抛物线法程序： clear format long a=0;b=1;n=258;h=(b-a)/n; A=0; for i=1:n xj=(i-1)*h; xi=i*h; xk=(xi+xj)/2; fxi=1/(1+xi*xi); fxj=1/(1+xj*xj); fxk=1/(1+xk*xk); A=A+(h/6)*(fxj+4*fxk+fxi); end disp(A) E=abs((B-pi/4)/(pi/4)); fprintf('The relative error between A and real-value is about: %d\\n',E) 答案： 0.785398163397449 The relative error between A and real-value is about: 2.827160e-016 从他们的相对误差值，我们可以看出，抛物线法精确程度最高，其次是矩形法，最后是梯形法。 dx，取n120．并 1x尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异． 二、分别用梯形法与抛物线法，计算 2（一） 梯形法： （1） 程序 clear format long; a=1;b=2;n=120;h=(b-a)/n; B=0; x1=1; y1=1/x1; B=h*y1/2; for i=1:n-1 xi=1+i*h; fxi=1/xi; B=B+h*fxi; end xn=2; yn=1/xn; B=B+h*yn/2; disp(B) 答案 0.684887057990131 （2）trapz x=1:1/120:2; y=1./x; trapz(x, y) 答案：ans =0.693151520800048 （二）抛物线法： （1） 程序 clear format long a=1;b=2;n=120;h=(b-a)/n; A=0; for i=1:n xj=1+(i-1)*h; xi=1+i*h; xk=(xi+xj)/2; fxi=1/xi; fxj=1/xj; fxk=1/xk; A=A+(h/6)*(fxj+4*fxk+fxi); end disp(A) 答案： 0.684848377754341 （2） quad f=inline('1./x'); I=quad(f, 1,2) 答案： I =0.693147199862970 从他们的最终结果可以看出，用梯形法与抛物线法求解的结果比直接使用函数trapz()、quad()进行计算求解的结果都要来得小。 使用函数trapz()与梯形法命令的结果接近，说明函数trapz()后台计算时用的是梯形法的算法。。。。。。同理） dx的近似计算结果与Matlab中各命令的计算结 01x2果相比较，试猜测Matlab中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例子支持你的观点． 三、将 1 Matlab中的数值积分命令最可能采用了中点法近似计算。 如下的梯形法中的左点法、中点法与右点法中中点法更接近Matlab中各命令的计算结果。 （1）clear format long; a=0;b=1;n=258;h=(b-a)/n; A=0;B=0;C=0; for i=1:n; xi=(i-1)*h; xj=i*h; xz=(xi+xj)/2; A=A+h*(1/(1+xi*xi)); B=B+h*(1/(1+xj*xj)); C=C+h*(1/(1+xz*xz)); end disp(A) disp(B) disp(C) 答案： 0.786366529681526 0.784428545185402 0.785398476379441 （2）x=0:1/100:1; y=1./(1+x.^2); trapz(x, y) 答案： ans = 0.785393996730782 可以看出中点法与答案更接近。（算出误差） 4、学习fulu2sum.m的程序设计方法，尝试用函数 sum 改写附录1和附录3的程序，避免for 循环． 附录1 format long n=100;a=0;b=1; inum1=0;inum2=0;inum3=0; syms x fx fx=1/(1+x^2); i=1:n；（不然会把结果直接输出来） xj=a+(i-1)*(b-a)/n; %左点 xi=a+i*(b-a)/n; %右点 fxj=subs(fx,'x',xj); %左点值 fxi=subs(fx,'x',xi); %右点值 fxij=subs(fx,'x',(xi+xj)/2); %中点值 e=fxj*(b-a)/n; f=fxi*(b-a)/n; g=fxij*(b-a)/n; inum1=sum(e); inum2=sum(f); inum3=sum(g); integrate=int(fx,0,1) integrate=double(integrate) fprintf('The relative error between inum1 and real-value is about: %d\\n\\n',... abs((inum1-integrate)/integrate)) fprintf('The relative error between inum2 and real-value is about: %d\\n\\n',... abs((inum2-integrate)/integrate)) fprintf('The relative error between inum3 and real-value is about: %d\\n\\n',... abs((inum3-integrate)/integrate)) 答案： integrate =pi/4 integrate = 0.785398163397448 The relative error between inum1 and real-value is about: 3.177794e-003 The relative error between inum2 and real-value is about: 3.188404e-003 The relative error between inum3 and real-value is about: 2.652582e-006 附录3 clear format long n=100;a=0;b=1;inum=0;ie=0 syms x fx fx=1/(1+x^2); for i=1:n xj=a+(i-1)*(b-a)/n; %左点 xi=a+i*(b-a)/n; %右点 xk=(xi+xj)/2; %中点 fxj=subs(fx,'x',xj); fxi=subs(fx,'x',xi); fxk=subs(fx,'x',xk); ie=ie+(fxj+4*fxk+fxi)*(b-a)/(6*n); end inum=sum(ie); integrate=int(fx,0,1) integrate=double(integrate) fprintf('The relative error between inum and real-value is about: %d\\n\\n',... abs((inum-integrate)/integrate)) 答案： ie = 0 integrate =pi/4 integrate = 0.785398163397448 The relative error between inum and real-value is about: 2.827160e-016 实验结果报告与实验总结 其实最重要的就是弄清楚第一题，第一题可以让我们对梯形法、矩形法与抛物线法的原理与程序有明确的概念，这样对后面的题目就有了一个基础。第二题要注意的就是x1与xn的值,应该令x1=1.xn=2，不然无法得出结论，还有就是xi=1+i*h中的加1，之前因为没加1调试了很久都无法成功，后来还是参考了同学的答案才恍然大悟。第四题比较难找不同，不过发现他们的规律（主要是需要用sum求得和面积）也就不太难调试了。 总体来说这些程序的调试就是需要仔细认真地对待每一个语句，一旦有一步出错，就会变得很头疼。 教师评语 实验成绩