[整理]定积分的近似计算

[整理]定积分的近似计算
实验⼆定积分的近似计算
⼀、问题背景与实验⽬的
利⽤⽜顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适⽤于被积函数的原函数能⽤初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的⽅法．在定积分的很多应⽤问题中，被积函数甚⾄没有解析表达式，可能只是⼀条实验记录曲线，或者是⼀组离散的采样值，这时只能应⽤近似⽅法去计算相应的定积分．
本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，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）
例：计算
0sin()d
x 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．矩形法
根据定积分的定义，每⼀个积分和都可以看作是定积分的⼀个近似值，即
1()d ()n
b
i i a i f x x f x ?==?∑?
在⼏何意义上，这是⽤⼀系列⼩矩形⾯积近似⼩曲边梯形的结果，所以把这个近似计算⽅法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才有⼀定的精确度．
针对不同i ?的取法，计算结果会有不同，我们以 12
0d 1x x +?为例（取100=n ），（1）左点法：对等分区间
b x i n
a b a x x a x n i =<<-+=<<<= 10，在区间],[1i i x x -上取左端点，即取1-=i i x ?，
1
2 01d ()1n i i i x f x x ?==?≈+∑?0.78789399673078，理论值 1
2 0d 14x x π=+?，此时计算的相对误差 0.7878939967307840.0031784
ππ-=≈（2）右点法：同（1）中划分区间，在区间],[1i i x x -上取右端点，即取i i x =?，
1
2 01d ()1n i i i x f x x ?==?≈+∑
0.78289399673078，理论值 1
2 0d 14x x π=+?，此时计算的相对误差 0.7828939967307840.0031884ππ-=≈
（3）中点法：同（1）中划分区间，在区间1[,]i i x x -上取中点，即取12
i i i x x ?-+=
， 12 01d ()1n i i i x f x x ?==?≈+∑
0.78540024673078，理论值 1
2 0d 14x x π=+?，此时计算的相对误差 60.785400246730784 2.653104
ππ--=
≈? 如果在分割的每个⼩区间上采⽤⼀次或⼆次多项式来近似代替被积函数，那么可以期望得到⽐矩形法效果好得多的近似计算公式．下⾯介绍的梯形法和抛物线法就是这⼀指导思想的产物．
2．梯形法
等分区间
b x i n a b a x x a x n i =<<-+=<<<= 10，n
a b x -=? 相应函数值为
n y y y ,,,10 （n i x f y i i ,,1,0),( ==）．
曲线)(x f y =上相应的点为
n P P P ,,,10 （n i y x P i i i ,,1,0),,( ==）
将曲线的每⼀段弧i i P P 1-⽤过点1-i P ，i P 的弦i i P P 1-（线性函数）来代替，这使得每个],[1i i x x -上的曲边梯形成为真正的梯形，其⾯积为
x y y i i ??+-2
1，n i ,,2,1 =．于是各个⼩梯形⾯积之和就是曲边梯形⾯积的近似值，
11 11
()d ()22n n b i i i i a i i y y x f x x x y y --==+?≈??=+∑∑?
，即 011 ()d ()22
b
n n a y y b a f x x y y n --≈++++?，称此式为梯形公式．仍⽤ 12
0d 1x x +?的近似计算为例，取100=n ， 10112 0d ()122
n n y y x b a y y x n --≈++++=+?0.78539399673078，理论值 12 0d 14
x x π=+?，此时计算的相对误差 60.785393996730784 5.305104
ππ--=≈? 很显然，这个误差要⽐简单的矩形左点法和右点法的计算误差⼩得多．
3．抛物线法
由梯形法求近似值，当)(x f y =为凹曲线时，它就偏⼩；当)(x f y =为凸曲线时，它就偏⼤．若每段改⽤与它凸性相接近的抛物线来近似时，就可减少上述缺点，这就是抛物线法．
将积分区间],[b a 作n 2等分，分点依次为
b x i n a b a x x a x n i =<<-+=<<<=2102 ，n
a b x 2-=?，对应函数值为
n y y y 210,,, （n i x f y i i 2,,1,0),( ==），
曲线上相应点为
n P P P 210,,, （n i y x P i i i 2,,1,0),,( ==）．
现把区间],[20x x 上的曲线段)(x f y =⽤通过三点),(000y x P ，),(111y x P ，),(222y x P 的抛物线
)(12x p x x y =++=γβα
来近似代替，然后求函数)(1x p 从0x 到2x 的定积分：
2
0 1 ()d x x p x x =?20 2 ()d x x x x x αβγ++=?)()(2)(3022
0223032x x x x x x -+-+-γβ
α
]4)(2)()()[(6
2022022202002γβαγβαγβα++++++++++-=x x x x x x x x x x 由于2
201x x x +=，代⼊上式整理后得 20 1 ()d x x p x x ?)](4)()[(6
12122202002γβαγβαγβα++++++++-=x x x x x x x x )4(621002y y y x x ++-=)4(6210y y y n
a b ++-= 同样也有 42 2 ()d x x p x x ?)4(6432y y y n
a b ++-=
……
222 ()d n n x n x p x x -?)4(621222n n n y y y n
a b ++-=-- 将这n 个积分相加即得原来所要计算的定积分的近似值： 222 22212 11()d ()d (4)6i i n n b x i i i i a x i i b a f x x p x x y y y n ---==-≈=++∑∑?
，即
021******* ()d [4()2()]6b n n n a b a f x x y y y y y y y y n
---≈
++++++++? 这就是抛物线法公式，也称为⾟⼘⽣（Simpson ）公式．仍⽤ 12
0d 1x x +?的近似计算为例，取100=n ， 102132124222 0d [4()2()]16n n n x b a y y y y y y y y x n ---≈+++++++++?
=0.78539816339745，理论值 12 0d 14x x π=+?，此时计算的相对误差 160.785398163397454 2.827104
ππ--=≈?
4. 直接应⽤Matlab 命令计算结果
（1）数值计算 12
0d .1x x +? ⽅法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）数值计算 2 1
2 0 1d d x x y 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
2 0d 1x x +?，取258=n ，并⽐较三种⽅法的精确程度．
2．分别⽤梯形法与抛物线法，计算 2 1d x x
，取120=n ．并尝试直接使⽤函数trapz()、quad()进⾏计算求解，⽐较结果的差异．
3．试计算定积分 0sin d x x x
+∞?．（注意：可以运⽤trapz()、quad()或附录程序求解吗？为什么？）
4．将 12
0d 1x x +?的近似计算结果与Matlab 中各命令的计算结果相⽐较，试猜测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',...。