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

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