数学实验报告实验序号：4 日期：2012 年12 月13 日实验名称定积分的近似计算问题背景描述：利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．实验目的：本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法。
对于定积分的近似数值计算，Matlab有专门函数可用。
实验原理与数学模型：1．矩形法根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度．针对不同的取法，计算结果会有不同。
（1）左点法：对等分区间，在区间上取左端点，即取。
（2）右点法：同（1）中划分区间，在区间上取右端点，即取。
（3）中点法：同（1）中划分区间，在区间上取中点，即取。
2．梯形法等分区间，相应函数值为（）．曲线上相应的点为（）将曲线的每一段弧用过点，的弦（线性函数）来代替，这使得每个上的曲边梯形成为真正的梯形，其面积为，．于是各个小梯形面积之和就是曲边梯形面积的近似值，，即，称此式为梯形公式。
3．抛物线法将积分区间作等分，分点依次为，，对应函数值为（），曲线上相应点为（）．现把区间上的曲线段用通过三点，，的抛物线来近似代替，然后求函数从到的定积分：由于，代入上式整理后得同样也有……将这个积分相加即得原来所要计算的定积分的近似值：，即这就是抛物线法公式，也称为辛卜生（Simpson）公式．实验所用软件及版本：Matlab 7.0主要内容（要点）：1．分别用梯形法与抛物线法，计算，取．并尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异．2．试计算定积分．（注意：可以运用trapz()、quad()或附录程序求解吗？为什么？）3．学习fulu2sum.m的程序设计方法，尝试用函数sum 改写附录1和附录3的程序，避免for 循环。
实验过程记录（含基本步骤、主要程序清单及异常情况记录等）：1：○1梯形法format longn=120;a=1;b=2;syms x fxfx=1/x;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); %所有右点值f=(fxi+fxj)/2*(b-a)/n; %梯形面积inum=sum(f) %加和梯形面积求解integrate=int(fx,1,2);integrate=double(integrate)fprintf('The relative error between inum and real-value is about:%g/n/n',...abs((inum-integrate)/integrate))【调试结果】inum =0.69315152080005integrate =0.69314718055995The relative error between inum and real-value is about:6.26164e-006/n/n○2抛物线法：%抛物线法format longn=120;a=1;b=2;inum=0;syms x fxfx=1/x;for i=1:nxj=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);inum=inum+(fxj+4*fxk+fxi)*(b-a)/(6*n);endinumintegrate=int(fx,1,2);integrate=double(integrate);fprintf('The relative error between inum and real-value is about:%g/n/n',...abs((inum-integrate)/integrate))【调试结果】inum =0.69314718056936The relative error between inum and real-value is about:1.35886e-011/n/n>>○3使用函数trapz()x=1:1/120:2;y=1./x;trapz(x,y)【调试结果】ans =0.69315152080005○4使用函数quad()quad('1./x',1,2)【调试结果】ans =0.693147199862972:○1使用函数trapz()x=1:1/120:inf;y=sin(x)./x;trapz(x,y)【调试结果】??? Error using ==> colonMaximum variable size allowed by the program is exceeded.○2使用函数quad()quad('sin(x)./x',0,inf)【调试结果】ans =NaN○3程序法%矩阵法format longn=inf;a=0;b=inf;syms x fxfx=sin(x)./x;i=1:n;xj=a+(i-1)*(b-a)/n; %左点xi=a+i*(b-a)/n; %右点xij=(xi+xj)/2;fxj=subs(fx,'x',xj); %左点值fxi=subs(fx,'x',xi); %右点值fxij=subs(fx,'x',xij); %中点值f1=fxj*(b-a)/n;f2=fxi*(b-a)/n;f3=fxij*(b-a)/n;inum1=sum(f1)inum2=sum(f2)inum3=sum(f3)integrate=int(fx,0,inf);integrate=double(integrate);fprintf('the relative error between inum1 and real-value is about: %g\n\n',...abs((inum1-integrate)/integrate))fprintf('the relative error between inum2 and real-value is about: %g\n\n',...abs((inum2-integrate)/integrate))fprintf('the relative error between inum3 and real-value is about: %g\n\n',...abs((inum3-integrate)/integrate))【调试结果】??? Maximum variable size allowed by the program is exceeded.○4使用matlab命令syms x;f=sin(x)/x;I=int(f,0,inf)【调试结果】I =1/2*pi3：○1矩形法：利用求和函数%矩阵法format longn=100;a=0;b=1;syms x fxfx=1/(1+x^2);i=1:n;xj=a+(i-1)*(b-a)/n; %左点xi=a+i*(b-a)/n; %右点xij=(xi+xj)/2;fxj=subs(fx,'x',xj); %左点值fxi=subs(fx,'x',xi); %右点值fxij=subs(fx,'x',xij); %中点值f1=fxj*(b-a)/n;f2=fxi*(b-a)/n;f3=fxij*(b-a)/n;inum1=sum(f1)inum2=sum(f2)inum3=sum(f3)integrate=int(fx,0,1);integrate=double(integrate);fprintf('the relative error between inum1 and real-value is about: %g\n\n',...abs((inum1-integrate)/integrate))fprintf('the relative error between inum2 and real-value is about: %g\n\n',...abs((inum2-integrate)/integrate))fprintf('the relative error between inum3 and real-value is about: %g\n\n',...abs((inum3-integrate)/integrate))【调试结果】inum1 =0.78789399673078inum2 =0.78289399673078inum3 =0.78540024673078the relative error between inum1 and real-value is about: 0.00317779the relative error between inum2 and real-value is about: 0.0031884the relative error between inum3 and real-value is about: 2.65258e-006○2抛物线法：使用求和函数%抛物线format longn=100;a=0;b=1;syms x fxfx=1/(1+x^2);i=1:n;xj=a+(i-1)*(b-a)/n; %左点xi=a+i*(b-a)/n; %右点xij=(xi+xj)/2;fxj=subs(fx,'x',xj); %左点值fxi=subs(fx,'x',xi); %右点值fxij=subs(fx,'x',xij); %中点值f=(fxj+4*fxij+fxi)*(b-a)/(6*n);inum=sum(f)integrate=int(fx,0,1);integrate=double(integrate);fprintf('the relative error between inum and real-value is about: %g\n\n',...abs((inum-integrate)/integrate))【调试结果】inum =0.78539816339745the relative error between inum and real-value is about: 2.82716e-016【情况记录】1、梯形法和抛物线法程序设计较为顺利。