佛山科学技术学院实验报告课程名称数值分析实验项目数值积分专业班级机械工程姓名余红杰学号2111505010指导教师陈剑成绩日期月日一、实验目的1、理解如何在计算机上使用数值方法计算定积分badxxf)(的近似值;2、学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。3、探索二重积分Ddxdyyxf),(在矩形区域},|),{(dycbxayxD的数值积分方法。二、实验要求(1)按照题目要求完成实验内容;(2)写出相应的Matlab程序;(3)给出实验结果(可以用表格展示实验结果);(4)分析和讨论实验结果并提出可能的优化实验。(5)写出实验报告。三、实验步骤1、用不同数值方法计算积分104ln9xxdx(1)取不同的步长h,分别用复合梯形及复合辛普森求积计算积分,给出误差中关于h的函数,并与积分精确值比较两公式的精度。(2)用龙贝格求积计算完成问题(1)。2、给出一种求矩形区域上二重积分的复化求积方法,然后计算二重积分xyDedxdy,其中积分区域{01,01}Dxy。1.%Int_t.m复化梯形:functionF=Int_t(x1,x2,n)%复化梯形求积公式%x1,x2为积分起点和中点%分为n个区间,没选用步长可以防止区间数为非整数。%样点矩阵及其函数值:x=linspace(x1,x2,n+1);y=f(x);m=length(x);%本题中用Matlab计算端点位置函数值为NaN,故化为零:y(1)=0;y(m)=0;%算出区间长度,步长h:h=(x2-x1)/n;a=[12*ones(1,m-2)1];%计算估计的积分值:F=h/2*sum(a.*y);%f.mfunctiony=f(x)y=sqrt(x).*log(x);%run11.mclc,clear;%分为10个区间,步长0.1的积分值:F=Int_t(0,1,10);F10=F%分为100个区间F=Int_t(0,1,100);F100=F%误差计算W10=abs((-4/9)-F10);W100=abs((-4/9)-F100);W=[W10W100]%复化辛普森:%Int_s.mfunctionF=Int_s(x1,x2,n)%复化梯形求积公式%x1,x2区间,分为n个区间。%样点矩阵及其函数值:x=linspace(x1,x2,n+1);y=f(x);m=length(x);h=(x2-x1)/n;y(1)=0;y(m)=0;%本题中用Matlab计算端点位置函数值为NaN,故化为零:F1=sum(y);xo=x+h/2;xo(m)=[];y=f(xo);F2=sum(y);F=(h/6)*(2*F1+4*F2);%run112.mclc,clear;%分为10个区间,步长0.1的积分值:F=Int_s(0,1,10);S10=F%分为100个区间F=Int_s(0,1,100);S100=F%误差计算W10=abs((-4/9)-S10);W100=abs((-4/9)-S100);W=[W10W100]%run113.m拟合误差和步长之间的三次曲线关系。clc,clear;%建立梯形误差、辛普森误差、步长矩阵:T=zeros(1,10);S=zeros(1,10);h=zeros(1,10);fori=1:10F=Int_t(0,1,10*i);T(i)=-4/9-F;F=Int_s(0,1,10*i);S(i)=-4/9-F;h(i)=1/(10*i);endTP=polyfit(h,T,3)SP=polyfit(h,S,3)%龙贝格:%Romberg.m:functionF=Romberg(x1,x2,n)可以明显看出其精度高于复化梯形也可以明显看出辛普森误差曲线各项系数都较小的%建立龙贝格推算矩阵、求最初步长:R=zeros(4);h=(x2-x1)/n;x=linspace(x1,x2,n+1);%计算矩阵第一列:复化梯形结果:fori=1:4F=Int_t(x1,x2,n*i);R(i,1)=F;end%计算第二列:辛普森fori=1:3R(i,2)=(4/3)*R(i+1,1)-(1/3)*R(i,1);end%计算第三列:复化斯科特fori=1:2R(i,3)=(16/15)*R(i+1,2)-(1/15)*R(i,2);endR(1,4)=(64/63)*R(2,3)-(1/63)*R(1,3);F=R(1,4);%run12.mclc,clear;F=Romberg(0,1,10);F10=F;F=Romberg(0,1,100);F100=F;[F10,F100;-4/9-F10,-4/9-F100]%右图为初始划分为10个区间和100个区间进行运算的结果,可以看出初始10次划分的精度比辛普森和梯形结果高出不少2,我选用的是类似于复化梯形的求积方法,由于是正方形的积分区域,把它划分为n*n块,在其中1.(x(k),y(k)),2.(x(k+1),y(k)),3.(x(k),y(k+1)),4.(x(k+1),y(k+1))的小区间上ds=dxdy是已知的,而且存在一e个f(e)ds=该函数在此小区域的数值,此处就选用(f(1),f(2),f(3),f(4))/4的值来近似,可知随着n的增大,是会越来越逼近真实值的。程序如下:%qes2.m被积分函数:functionf=qes2(x,y)f=exp(-x*y);%Int_xy.m积分代码:functionF=Int_xy(x1,x2,y1,y2,n)%仅适合此题目的正方形区域二重积分%x1,x2,y1,y2分别为横坐标和纵坐标起点与中点;%n为划分为n*n的小正方形网格%建立所有网格交点处的横坐标和纵坐标:x=linspace(x1,x2,n+1);y=linspace(y1,y2,n+1);%计算每个小方格的面积,dxdy:ds=(x2-x1)*(y2-y1)/(n^2);%建立一个矩阵存放各个小方块的积分值:Z=zeros(n);%初始化积分值:F=0;%把小方块的积分值存入矩阵Z中,第一行第一列存放最左下角方格数据,以此类推:%并依次把小方格区域的积分值累加,得到最终积分值。fori=1:nforj=1:nf=qes2(x(i),y(j));f1=f;f=qes2(x(i+1),y(j));f2=f;f=qes2(x(i),y(j+1));f3=f;f=qes2(x(i+1),y(j+1));f4=f;Z(i,j)=(ds*(f1+f2+f3+f4))/4;F=F+Z(i,j);endend%run2.m运行代码:clc,clear;F=Int_xy(0,1,0,1,2);%划分为4个小格F2=F;F=Int_xy(0,1,0,1,10);%划分为100个小格F10=F;F=Int_xy(0,1,0,1,100);%划分为10000个小格:F100=F;[F2;F10;F100]%得到的积分值如右图,可见是逐渐精确的: