试验三 数值积分

佛山科学技术学院实验报告课程名称数值分析实验项目常微分方程问题初值问题数值解法专业班级姓名学号指导教师成绩日期试验三数值积分一、实验目的1、理解如何在计算机上使用数值方法计算定积分badxxf)(的近似值;2、学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。3、学会利用Matlab提供的积分函数求二重积分。函数名含义函数名含义trapz梯形求积法dblquad二重积分quadSimpson求积法triplequad三重积分二、实验要求（1）按照题目要求完成实验内容；（2）写出相应的Matlab程序；（3）给出实验结果(可以用表格展示实验结果)；（4）分析和讨论实验结果并提出可能的优化实验。（5）写出实验报告。三、实验步骤1、用不同数值方法计算积分10sin0.9460831xdxx（1）取不同的步长h，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两公式的精度，是否存在一个最小的h，使得精度不能被改善？（2）用龙贝格求积计算完成问题（1）。2、计算二重积分xyDedxdy（1）若积分区域{01,01}Dxy；（2）若积分区域22{1,0,0}Dxyxy求积公式参考代码：1、复化梯形求积公式functionI=T_quad(x,y)%复化梯形求积公式，其中%x---向量，被积函数自变量的等距节点%y---向量，被积函数在节点处的函数值n=length(x);m=length(y);ifn~=merror('ThelengthsofXandYmustbeequal');return;endh=(x(n)-x(1))/(n-1);a=[12*ones(1,n-2)1];I=h/2*sum(a.*y);2、复化Simpson求积公式functionI=S_quad(x,y)%复化Simpson求积公式，其中%x---向量，被积函数自变量的等距节点%y---向量，被积函数在节点处的函数值n=length(x);m=length(y);ifn~=merror('ThelengthsofXandYmustbeequal');return;endifrem(n-1,2)~=0%如果n-1不能被2整除，则调用复化梯形公式warning('给出的点数不适用于Simpson公式，调用梯形公式求积分');I=T_quad(x,y);return;endN=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n);%如果等分成N段，且从1开始计算，应该是公式N=(n-1)/2;fork=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1;%如果下标k从1开始计算，则每小段区间的下标应该为1、3、5、7、、、，中间节点%的下标分别是2、4、6、、，对应书本P107公式（3.4），它们的系数分别为1，4，1；endI=h/6*sum(a.*y);3、Romberg求积公式functionI=R_quad_iter(fun,a,b,ep)%Romberg求积公式，其中%fun---被积函数，在调用之前需要用inline()函数定义其为内置函数。%a,b---积分区间的端点，要求ab%ep---精度要求，省缺为1e-5%注意函数feval()及函数inline()的使用。%nargin指设定参数的个数ifnargin4ep=1e-5;end;m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while1N=2^(m-1);h=h/2;I=I/2;fori=1:NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;whileM1;T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);M=M/2;k=k+1;endifabs(T(k,k)-T(k-1,k-1))epbreak;endm=m+1;endI=T(k,k);实验三参考答案：1.1用复合梯形公式计算积分10sinxdxx（用不同的步长h=0.05和0.001,与精确值0.9460831比较）(1)步长h=0.05x=0.0001:0.1255:1;y=sin(x)./x;I=T_quad(x,y)I=0.841326530364063%也可以用eps代替0.0001，表示一个很小的正数（10-15）x=eps:0.125:1;y=sin(x)./x;I=T_quad(x,y)I=0.9456908635827011.2用复合辛普森求积公式计算积分10sinxdxx。x=0.001:0.0001:1;y=sin(x)./x;I=S_quad(x,y)I=0.945083070422740x=0.0001:0.0001:1;y=sin(x)./x;I=S_quad(x,y)Warning:给出的点数不适用于Simpson公式，调用梯形公式求积分InS_quadat13I=0.9459830701162941.3用龙贝格算法计算积分10sinxdxxfun=inline('sin(x)./x')I=R_quad_iter(fun,0.0000001,1)（因x=0时，函数没意义，但求极限为0）T=0.920735400330398000.9397931867634230.94614578224109800.9445134221532120.9460868339498080.9460829040637220.94569076370456400I=0.9460829703872232.1利用dblquad()函数计算矩形区域上的二重积分xyDedxdy，其中{01,01}Dxy。fun=inline('exp(-x.*y)','x','y');I=dblquad(fun,0,1,0,1，1e-6)%dblquad()函数适用于计算矩形区域上的二重积分。I=0.7965995997159532.2利用dblquad()函数计算非矩形区域上的二重积分xyDedxdy，其中22{1,0,0}Dxyxy。fun=inline('exp(-x.*y).*(x.^2+y.^2=1)','x','y');I=dblquad(fun,0,1,0,1,1e-6)I=0.675161691385222四、实验结果MATLABdesktopkeyboardshortcuts,suchasCtrl+S,arenowcustomizable.Inaddition,manykeyboardshortcutshavechangedforimprovedconsistencyacrossthedesktop.Tocustomizekeyboardshortcuts,usePreferences.Fromthere,youcanalsorestorepreviousdefaultsettingsbyselectingR2009aWindowsDefaultSetfromtheactivesettingsdrop-downlist.Formoreinformation,seeHelp.Clickhereifyoudonotwanttoseethismessageagain.h=0.05;x=0.0001:0.1255:1;y=sin(x)./x;I=T_quad(x,y)I=0.841326530364063x=eps:0.125:1;y=sin(x)./x;I=T_quad(x,y)I=0.945690863582701x=0.001:0.0001:1;y=sin(x)./x;I=S_quad(x,y)I=0.945083070422740x=0.0001:0.0001:1;y=sin(x)./x;I=S_quad(x,y)Warning:给出的点数不适用于Simpson公式，调用梯形公式求积分InS_quadat13I=0.945983070116294fun=inline('sin(x)./x')fun=Inlinefunction:fun(x)=sin(x)./xI=R_quad_iter(fun,0.0000001,1)T=0.9207354003303980.939793186763423T=0.92073540033039800.9397931867634230.9461457822410980.9445134221532120T=0.920735400330398000.9397931867634230.94614578224109800.9445134221532120.9460868339498080.9460829040637220.94569076370456400I=0.946082970387223fun=inline('exp(-x.*y)','x','y');I=dblquad(fun,0,1,0,1,1e-6)I=0.796599599715953fun=inline('exp(-x.*y).*(x.^2+y.^2=1)','x','y');I=dblquad(fun,0,1,0,1,1e-6)I=0.675161691385222五、讨论分析求积分10sin0.9460831xdxx是x=0是函数没有意义，实验中采用0.0001逼近0来近似计算积，这样做是存在很大误差的。六、改进实验建议可以采用书本的逐步逼近法取不同的rn=2-n然后使用Romberg算法计算。就可以得到比较理想的计算结果。