实验三定积分的近似计算一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，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()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)这时必须存在一个函数文件：function z = integrnd(x, y)z = y*sin(x);7．fprintf（文件地址，格式，写入的变量）：把数据写入指定文件．例：x = 0:.1:1;y = [x; exp(x)];fid = fopen('','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 ()nbi i a i f x x f x ς==∆∑⎰在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度．针对不同i ς的取法，计算结果会有不同，我们以 12 0d 1x x +⎰为例（取100=n ）， （1） 左点法：对等分区间 b x i na b a x x a x n i =<<-+=<<<=ΛΛ10， 在区间],[1i i x x -上取左端点，即取1-=i i x ς，12 01d ()1n i i i x f x x ς==∆≈+∑⎰， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 0.7878939967307840.0031784ππ-=≈ （2）右点法：同（1）中划分区间，在区间],[1i i x x -上取右端点，即取i i x =ς，12 01d ()1n i i i x f x x ς==∆≈+∑⎰， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 0.7828939967307840.0031884ππ-=≈（3）中点法：同（1）中划分区间，在区间1[,]i i x x -上取中点，即取12i i i x x ς-+=， 12 01d ()1n i i i x f x x ς==∆≈+∑⎰，理论值 12 0d 14x x π=+⎰，此时计算的相对误差 60.785400246730784 2.653104ππ--=≈⨯ 如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那么可以期望得到比矩形法效果好得多的近似计算公式．下面介绍的梯形法和抛物线法就是这一指导思想的产物．2． 梯形法等分区间b x i n a b a x x a x n i =<<-+=<<<=ΛΛ10，na 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 ∆⨯+-21，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 ()22b n n a y y b a f x x y y n --≈++++⎰L ， 称此式为梯形公式． 仍用 120d 1x x +⎰的近似计算为例，取100=n ，10112 0d ()122n n y y x b a y y x n --≈++++=+⎰L ， 理论值 12 0d 14x 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ΛΛ，na 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 的定积分：20 1 ()d x x p x x =⎰20 2 ()d x x x x x αβγ++=⎰)()(2)(30220223032x x x x x x -+-+-γβα ]4)(2)()()[(62022022202002γβαγβαγβα++++++++++-=x x x x x x x x x x 由于2201x x x +=，代入上式整理后得20 1 ()d x x p x x ⎰)](4)()[(612122202002γβαγβαγβα++++++++-=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 na 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---≈++++++++⎰L L 这就是抛物线法公式，也称为辛卜生（Simpson ）公式． 仍用 120d 1x x +⎰的近似计算为例，取100=n ， 102132124222 0d [4()2()]16n n n x b a y y y y y y y y x n ---≈+++++++++⎰L L =， 理论值 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::1;y=1./(1+x.^2);trapz(x,y) （梯形法求数值积分）（2）数值计算 2 12 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． 实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算 12 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． 学习的程序设计方法，尝试用函数 sum 改写附录1和附录3的程序，避免for 循环．五、附录附录1：矩形法（左点法、右点法、中点法）（）format longn=100;a=0;b=1;inum1=0;inum2=0;inum3=0;syms x fxfx=1/(1+x^2);for i=1:nxj=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;endinum1inum2inum3integrate=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：梯形法（）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=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);endinumintegrate=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：梯形法（），利用求和函数，避免for 循环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; %所有右点的数组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：抛物线法（）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);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,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))。