试验定积分的近似计算

实验二定积分的近似计算一、问题背景与实验目的利用牛顿一莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数达出来的情形.如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法.在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分.本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法.对于定积分的近似数值计算，Matlab有专门函数可用.二、相关函数(命令)及简介sum(a)：求数组a的和.formatlong：长格式，即屏幕显示15位有效数字.(注：由于本实验要比较近似解法和精确求解间的误差，需要更高的精度).double():若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值.quad()：抛物线法求数值积分.格式：quad(fun,a,b)，注意此处的fun是函数，并且为数值形式的，所以使用*、人A等运算时要在其前加上小数点，即.*、./、a等.例：Q=quad('1./(x.A3-2*x-5)',0,2);trapz()：梯形法求数值积分.格式：trapz(x,y)其中x为带有步长的积分区间；y为数值形式的运算(相当于上面介绍的函数fun)例：计算°sin(x)dxx=0:pi/100:pi;y=sin(x);trapz(x,y)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：functionz=integrnd(x,y)z=y*sin(x);fprintf(文件地址，格式，写入的变量)：把数据写入指定文件.例：x=0:.1:1;y=[x;exp(x)];fid=fopen('exp.txt','w');%打开文件fprintf(fid,'%6.2f%12.8f\n',y);%写入fclose(fid)%关闭文件syms变量1变量2…：定义变量为符号.sym('表达式'):将表达式定义为符号.解释：Matlab中的符号运算事实上是借用了Maple的软件包，所以当在Matlab中要对符号进行运算时，必须先把要用到的变量定义为符号.int(f,v,a,b):求f关于v积分，积分区间由a到b.subs(f,'x',a)：将a的值赋给符号表达式f中的x，并计算出值.若简单地使用subs(f)，则将f的所有符号变量用可能的数值代入，并计算出值.三、实验矩形法根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即nf(i)逐f(x)dx八iT在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法.不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度.1Hv针对不同i的取法，计算结果会有不同，我们以半为例(取n=100)，‘°1+x(1)左点法：对等分区间::Xn=b，b—a.x。—a■■■X1叮Xj=ai二n在区间[Xi」，Xi]上取左端点，即取i，1dx01■X20.78789399673078n工為f(i)Xi=1理论值.:丹盲，此时计算的相对误差0.7878939967307〜珥4吋40.003178(2)右点法：同(1)中划分区间，在区间[Xy,Xi]上取右端点，即取i=Xj，1dx01x2八f(J.：Xj：0.78289399673078i4(3)中点法：同0.78289399673079jt/4叫4:0.003188(1)中划分区间，在区间［Xy,X］上取中点，即取i二Xi」Xi21dx1X2n八f(J.Xj0.78540024673078i4理论值匸寻冷，此时计算的相对误差0.78540024673078-二4u/4:2.65310》理论值匸無冷，此时计算的相对误差如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那么可以期望得到比矩形法效果好得多的近似计算公式.下面介绍的梯形法和抛物线法就是这一指导思想的产物.梯形法等分区间b—a•….血b—ax°二a::.x〔：：•：:•Xj二ai：：•::.Xn—b,—x=nn相应函数值为y。，％,,yn(%=f(xj,i=0,1,,n).曲线y二f(x)上相应的点为Po,R,…,Pn(Pi=(Xi,y」,i=0,1，…,n)将曲线的每一段弧P^Pi用过点Pi4,P的弦PjR(线性函数)来代替，这使得每个［X2,Xi］上的曲边梯形成为真正的梯形，其面积为i=1,2,,n.于是各个小梯形面积之和就是曲边梯形面积的近似值,n(YijYi)f(x)dx：'i4即bb_ayoynTOC\o"1-5"\h\zaf(x)dx(片yi•|l(-Ynj-),an22称此式为梯形公式.idx仍用.的近似计算为例'取“00,01x2：上^(西％yn」北)=0.78539399673078n22理论值.01溼=4，此时计算的相对误差0.78539399673078-叫46=st5.305x10二4很显然，这个误差要比简单的矩形左点法和右点法的计算误差小得多.抛物线法由梯形法求近似值，当y=f(x)为凹曲线时，它就偏小；当y=f(x)为凸曲线时，它就偏大•若每段改用与它凸性相接近的抛物线来近似时，就可减少上述缺点，这就是抛物线法.将积分区间［a,b］作2n等分，分点依次为b—a•….血b—ax°二a：：•x〔：：•：:.Xj二ai：：•::.x?n二b，=x=■2n2n对应函数值为y°,yi,，Y2n(Yi=f(xj,i=0,1,…,2n)，曲线上相应点为P°,Pi,…,P2n(P=(Xi,yi),i=0,1，…,2n).现把区间［x°,X2］上的曲线段y二f(x)用通过三点P)(X0,y0)，Pi(Xi,yi)，B(X2,y2)的抛物线y=x2x=pi(x)来近似代替，然后求函数Pi(x)从X0到X2的定积分：X2勺2门0(33^22■Pi(x)dx二Cx2】x)dx(X2-X0)(X2-X0)(X2-X。)■X0■X032X2x°[(:xf亠.x0)(：x；亠〉x2)亠：(x0x2)22：(x0x2)4]X0产，代入上式整理后得X2j1(X)dxX2Xo[(：x2：x0)(_：ix；】；X2：：；'")4(：xj】jx—)]6X2-X。b-a丁(yo4y1y»6n(y04y1y2)同样也有x4b-ax2nb-aJPn(x)dx=^^(y22+4y2n_L+y2n)x2nz6n将这n个积分相加即得原来所要计算的定积分的近似值:jbf(x)dx肚£「Pi(x)dx=£宁(y^+Ayz+yzi)，ay勺隹i16n即baf（x）dx[y°y2n-4（力y3川y2n」）2皿丫4川丫22）]a6n这就是抛物线法公式，也称为辛卜生（Simpson）公式.仍用1-dxT的近似计算为例，取n=100,L01+x21dx1x2b-as?b-a:右[y0y2n4(八山"山2(y2"d]=0.78539816339745,蔦，此时计算的相对误差0.7853981633974&叼42.82710'6（符号求积分）（抛物线法求数值积分）4.直接应用Matlab命令计算结果(1)数值计算1-^.'01+x2方法1:int('1/(1+xA2)','x',0,1)方法2：quad('1./(1+x.A2)',0,1)方法3:x=0:0.001:1;y=1./(1+x.A2);trapz(x,y)(梯形法求数值积分)(2)数值计算0dx^xy2dy方法1:int(int('x+yA2','y',-1,1),'x',0,2)(符号求积分)方法2：dblquad(inline('x+yA2'),0,2,-1,1)(抛物线法二重数值积分)四、自己动手实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算1%,01+x2取n=258，并比较三种方法的精确程度.2dx分别用梯形法与抛物线法，计算坐，取n=120.并尝试直接使用函数'1xtrapz()、quad()进行计算求解，比较结果的差异.试计算定积分sin-xdx.(注意：可以运用trapz()、quad()或附录程序求解、0x吗？为什么？)将1%的近似计算结果与Matlab中各命令的计算结果相比较，试猜测」01+x2Matlab中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例子支持你的观点.通过整个实验内容及练习，你能否作出一些理论上的小结，即针对什么类型的函数(具有某种单调特性或凹凸特性)，用某种近似计算方法所得结果更接近于实际值？学习fulu2sum.m的程序方法，尝试用函数sum改写附录1和附录3的程序，避免for循环.五、附录附录1:矩形法(左点法、右点法、中点法)(fulu1.m)formatlongn=100;a=0;b=1;inum1=0;inum2=0;inum3=0;symsxfxfx=1/(1+xA2);%左点%右点%左点值fori=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('Therelativeerrorbetweeninum1andreal-valueisabout:%d\n\n',...abs((inum1-integrate)/integrate))fprintf('Therelativeerrorbetweeninum2andreal-valueisabout:%d\n\n',...abs((inum2-integrate)/integrate))fprintf('Therelativeerrorbetweeninum3andreal-valueisabout:%d\n\n',...abs((inum3-integrate)/integrate))附录2：梯形法(fulu2.m)formatlongn=100;a=0;b=1;inum=0;symsxfxfx=1/(1+xA2);fori=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('Therelativeerrorbetweeninumandreal-valueisabout:%d\n\n',...abs((inum-integrate)/integrate))附录2sum：梯形法(fulu2sum.m),利用求和函数，避免for循环formatlongn=100;a=0;b=1;symsxfxfx=1/(1+xA2);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('Therelativeerrorbetweeninumandreal-valueisabout:%d\n\nabs((inum-integrate)/integrate))附录3：抛物线法(fulu3.m)formatlongn=100;a=0;b=1;inum=0;symsxfxfx=1/(1+xA2);fori=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('Therelativeerrorbetweeninumandreal-valueisabout:%d\n\n',...abs((inum-integrate)/integrate))