第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
matlab定积分的近似计算.doc
实验二 定积分的近似计算
一、问题背景与实验目的
利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被
积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办
到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚
至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这
时只能应用近似方法去计算相应的定积分．
本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线
法．对于定积分的近似数值计算，Matlab 有专门函数可用．
二、相关函数（命令）及简介
1．sum(a)：求数组 a 的和．
2．format long：长格式，即屏幕显示 15 位有效数字．
（注：由于本实验要比较近似解法和精确求解间的误差，需要更高的精度）．
3．double()：若输入的是字符则转化为相应的 ASCII 码；若输入的是整型数值则
转化为相应的实型数值．
4．quad()：抛物线法求数值积分．
格式： quad(fun,a,b) ，注意此处的 fun 是函数，并且为数值形式的，所以使
用*、/、^等运算时要在其前加上小数点，即 .*、./、.^等．
例：Q = quad('1./(x.^3-2*x-5)',0,2);
5．trapz()：梯形法求数值积分．
格式：trapz(x,y)
其中 x 为带有步长的积分区间；y 为数值形式的运算（相当于上面介绍的函
数 fun）
例：计算
0
sin( )dx x
x=0:pi/100:pi;y=sin(x);
trapz(x,y)
6．dblquad()：抛物线法求二重数值积分．
格式：dblquad(fun,xmin,xmax,ymin,ymax)，fun 可以用 inline 定义，也可以通
过某个函数文件的句柄传递．
例 1：Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)
顺便计算下面的 Q2，通过计算，比较 Q1 与 Q2 结果（或加上手工验算），
找出积分变量 x、y 的上下限的函数代入方法．
Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi)
例 2：Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)
这时必须存在一个函数文件 integrnd.m：
５３
function z = integrnd(x, y)
z = y*sin(x);
7．fprintf（文件地址，格式，写入的变量）：把数据写入指定文件．
例：x = 0:.1:1;
y = [x; exp(x)];
fid = fopen('exp.txt','w'); %打开文件
fprintf(fid,'%6.2f %12.8f\n',y);
fclose(fid)
%写入
%关闭文件
8．syms 变量 1 变量 2 …：定义变量为符号．
9．sym('表达式')：将表达式定义为符号．
解释：Matlab 中的符号运算事实上是借用了 Maple 的软件包，所以当在 Matlab
中要对符号进行运算时，必须先把要用到的变量定义为符号．
10．int(f,v,a,b)：求 f 关于 v 积分，积分区间由 a 到 b．
11．subs(f，'x'，a)：将 a 的值赋给符号表达式 f 中的 x，并计算出值．若简单
地使用 subs(f)，则将 f 的所有符号变量用可能的数值代入，并计算出值．
三、实验内容
1． 矩形法
根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即
b
a
( )d
f x x
n
i
1
f
(
x
i
)
i
在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这
个近似计算方法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才
有一定的精确度．
针对不同 i 的取法，计算结果会有不同，我们以
1
0
d
x
x
2
1
为例（取 100
n
），
（1）左点法：对等分区间
x
1
a
x
0
x
i
a
在区间
[
x
1
i
,
x
i
]
上取左端点，即取
x
i
1 i
i
x
n
b
，
ab
n
，
理论值
1
0
1
0
d
x
x
2
1
n
i
1
f
(
i
)
x
i
0.78789399673078，
d
x
x
2
1
4
，此时计算的相对误差
0.78789399673078
4
4
0.003178
（2）右点法：同（1）中划分区间，在区间
[
x
1
i
,
x
i
]
上取右端点，即取
x
i
i
，
５４
理论值
1
0
1
0
d
x
x
2
1
n
i
1
f
(
i
)
x
i
0.78289399673078，
d
x
x
2
1
4
，此时计算的相对误差
0.78289399673078
4
4
0.003188
（3）中点法：同（1）中划分区间，在区间 1
x
i
[
,
x
i
]
上取中点，即取
i
x
i
x
i
1
2
，
理论值
1
0
1
0
d
x
x
2
1
n
i
1
f
(
i
)
x
i
0.78540024673078，
d
x
x
2
1
4
，此时计算的相对误差
0.78540024673078
4
4
2.653 10
6
如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那
么可以期望得到比矩形法效果好得多的近似计算公式．下面介绍的梯形法和抛物
线法就是这一指导思想的产物．
2． 梯形法
等分区间
相应函数值为
x
0
a
x
1
x
i
a
ab
n
i
x
n
b
，
x
ab
n
y
,
0 （
ny
y
1
,
,
曲线
y
)(xf
上相应的点为
,
PP
0 （
1
nP
,
,
y
i
(
xf
i
),
i
,1,0
,
n
）．
P
i
(
,
yx
i
i
),
i
,1,0
,
n
）
将曲线的每一段弧
i PP 1 用过点 1iP ， iP 的弦
i
i PP 1 （线性函数）来代替，这使得
i
每个
[
x
1
i
,
x
i
]
上的曲边梯形成为真正的梯形，其面积为
y
i
y
i
1
2
x
，
i
,2,1
,
n
．
于是各个小梯形面积之和就是曲边梯形面积的近似值，
５５
b
a
( )d
f x x
n
i
1
y
i
y
i
1
2
x
x
2
n
i
1
(
y
y
)
i
i
1
，
即
b
称此式为梯形公式．
a
( )d
f x x
b a
n
(
y
0
2
y
1
y
n
1
y
n
2
)
，
仍用
1
0
d
x
2
x
d
x
x
2
1
1
0
d
1
x
x
2
4
理论值
1
0
1
n
的近似计算为例，取 100
y
n
2
b a
n
y
0
2
y
1
y
(
1
n
，
)
0.78539399673078，
，此时计算的相对误差
0.78539399673078
4
4
5.305 10
6
很显然，这个误差要比简单的矩形左点法和右点法的计算误差小得多．
3． 抛物线法
由梯形法求近似值，当
y
)(xf
为凹曲线时，它就偏小；当
y
)(xf
为凸曲
线时，它就偏大．若每段改用与它凸性相接近的抛物线来近似时，就可减少上述
缺点，这就是抛物线法．
将积分区间 ],[ ba 作 n2 等分，分点依次为
x
0
a
x
1
x
i
a
ab
2
n
i
x
2
n
b
，
x
ab
2
n
，
对应函数值为
曲线上相应点为
y
0
, （
y
1
ny
2
,
,
y
i
(
xf
i
),
i
,1,0
2,
n
），
PP
0
1
, （
nP
2
,
,
P
i
(
,
yx
i
i
),
i
,1,0
2,
n
）．
现 把 区 间
[
,
0 xx
2
]
上 的 曲 线 段
y
)(xf
用 通 过 三 点
(
xP
0
0
,
y
0
)
，
,
yxP
1
1
(
1
)
，
(
xP
2
2
,
y
2
)
的抛物线
y
2
x
x
)(1
xp
来近似代替，然后求函数
x
2
x
0
p x x
1
( )d
x
2
x
0
2
(
x
x
)(1 xp 从 0x 到 2x 的定积分：
2
(
3
)d
x
x
x
(
)
3
0
3
2
x
2
2
2
x
0
)
(
x
2
x
0
)
５６
x
2
x
0
6
[(
(
x
x
x
x
(
)
)
2
0
2
2
0
2
x
0
x
2
2
)
(2
x
0
x
2
]4)
由于
x
1
x
0
x
2
2
，代入上式整理后得
x
2
x
0
( )d
p x x
1
x
2
x
0
6
[(
)]
(4)
x
x
x
x
(
)
x
1
2
0
2
0
2
2
2
x
1
x
2
x
0
6
(
y
0
4
y
1
y
2
)
ab
6
n
(
y
0
4
y
1
y
2
)
同样也有
x
4
x
2
( )d
p x x
2
x
2
n
x
2
n
2
( )d
p x x
n
ab
6
n
……
ab
6
n
y
(
(
y
2
4
y
3
y
4
)
4
y
2
n
1
y
2
n
)
2
n
2
将这 n 个积分相加即得原来所要计算的定积分的近似值：
b
a
( )d
f x x
n
i
1
即
x
2
i
x
2 2
i
( )d
p x x
i
n
i
1
b a
6
n
(
y
4
y
y
)
2
i
2 1
i
，
2
i
2
b
( )d
f x x
b a
6
n
a
2
这就是抛物线法公式，也称为辛卜生（Simpson）公式．
1
0
2
3
2
n
n
[
y
y
4(
y
1
y
y
) 2(
y
y
4
y
2
n
2
)]
) 2(
y
2
y
4
y
2
n
2
)]
仍用
1
0
1
0
1
d
x
x
x
x
2
1
d
2
n
的近似计算为例，取 100
b a y
[
6
n
4(
y
3
y
1
y
0
y
，
2
n
=0.78539816339745，
1
2
n
理论值
1
0
d
x
x
2
1
4
，此时计算的相对误差
0.78539816339745
4
4
2.827 10
16
4. 直接应用 Matlab 命令计算结果
（1） 数值计算
1
0
d
x
x
1
.
2
方法 1：int('1/(1+x^2)','x',0,1)
方法 2：quad('1./(1+x.^2)',0,1) （抛物线法求数值积分）
(符号求积分)
５７
方法 3：x=0:0.001:1;
y=1./(1+x.^2);
trapz(x,y)
2
d
（2）数值计算
x
0
1
1
（梯形法求数值积分）
x
y
2
d
y
方法 1：int(int('x+y^2','y',-1,1),'x',0,2) （符号求积分）
方法 2：dblquad(inline('x+y^2'),0,2,-1,1) （抛物线法二重数值积分）
四、自己动手
1． 实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算
1
0
d
x
x
2
1
，
取 258
n
，并比较三种方法的精确程度．
2． 分别用梯形法与抛物线法，计算
2
1
dx
x
trapz()、quad()进行计算求解，比较结果的差异．
，取 120
n
．并尝试直接使用函数
．（注意：可以运用 trapz()、quad()或附录程序求解
3． 试计算定积分
0
吗？为什么？）
sin dx x
x
4． 将
1
0
d
x
x
1
的近似计算结果与 Matlab 中各命令的计算结果相比较，试猜测
2
Matlab 中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例
子支持你的观点．
5． 通过整个实验内容及练习，你能否作出一些理论上的小结，即针对什么类型
的函数（具有某种单调特性或凹凸特性），用某种近似计算方法所得结果更接
近于实际值？
6． 学习 fulu2sum.m 的程序设计方法，尝试用函数 sum 改写附录 1 和附录 3 的
程序，避免 for 循环．
五、附录
附录 1：矩形法（左点法、右点法、中点法）（fulu1.m）
format long
n=100;a=0;b=1;
inum1=0;inum2=0;inum3=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;
fxj=subs(fx,'x',xj);
%左点
%右点
%左点值
５８
fxi=subs(fx,'x',xi);
fxij=subs(fx,'x',(xi+xj)/2);
inum1=inum1+fxj*(b-a)/n;
inum2=inum2+fxi*(b-a)/n;
inum3=inum3+fxij*(b-a)/n;
%右点值
%中点值
end
inum1
inum2
inum3
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))
附录 2：梯形法（fulu2.m）
format long
n=100;a=0;b=1;inum=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;
fxj=subs(fx,'x',xj);
fxi=subs(fx,'x',xi);
inum=inum+(fxj+fxi)*(b-a)/(2*n);
end
inum
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))
５９
附录 2sum：梯形法（fulu2sum.m），利用求和函数，避免 for 循环
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;
fxj=subs(fx,'x',xj);
fxi=subs(fx,'x',xi);
f=(fxi+fxj)/2*(b-a)/n;
inum=sum(f)
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))
附录 3：抛物线法（fulu3.m）
format long
n=100;a=0;b=1;inum=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);
inum=inum+(fxj+4*fxk+fxi)*(b-a)/(6*n);
end
inum
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))
６０
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
