MATLAB数值分析试验二复合梯形、辛普森和龙贝格求积,以及二重积分计算等

1、 佛山科学技术学院 实验报告 课程名称 _ 数值分析 实验工程 _ 数值积分 专业班级机械工程姓名余红杰学号2111505010指导教师陈剑成绩 日期月曰 一、实验目的 b 1、理解如何在计算机上使用数值方法计算定积分 a f X的近似值； 2、 学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。 3、 探索二重积分.11 f (x, y) dxdy 在矩形区域 D = ( x, y) | a _ x _ b, c _ y d的数值 D 积分方法。 实验要求 (1)按照题目要求完成实验内容； (2)写出相应的Matlab程序； (3)给出实验结果(可以用表格展示实验结果)； (4

2、)分析和讨论实验结果并提出可能的优化实验。 (5)写出实验报告。 三、实验步骤 1、 用不同数值方法计算积 xln xdx =- 0 9 (1)取不同的步长h,分别用复合梯形及复合辛普森求积计算积分，给出误差中关于 h的 函数，弁与积分精确值比拟两公式的精度。 (2) 用龙贝格求积计算完成问题(1 )。 2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分 .eydxdy, 其中积分区域D二0乞x岂1,0岂y乞1 1. %lnt_t.m复化梯形： function F = Int_t(x1,x2,n) %复化梯形求积公式 % x1,x2为积分起点和中点 %分为n个区间，没选用步长可

3、以防止区间数为非整数。 %样点矩阵及其函数值： x = lin space(x1,x2 ,n+1); y = f(x); m = len gth(x); %此题中用Matlab计算端点位置函数值为 NaN,故化为零: y(1) = 0; y(m) = 0; %算出区间长度，步长 h： h = (x2 -x1)/n; a = 1 2*o nes(1,m-2) 1; %计算估计的积分值： F = h/2*sum(a.*y); %f.m fun cti on y = f(x) y = sqrt(x).*log(x); %run 11.m clc,clear; %分为10个区间，步长0.1的积分值：

4、F = In t_t(0,1,10); F10 = F %分为100个区间 F = In t_t(0,1,100); F100 = F %误差计算 W10 = abs(-4/9)-F10); W100 = abs(-4/9)-F100); W = W10 W100 F10 - -0. 4171 F10Q = -0.4421 W = 0.0八74 0.0013 %复化辛普森： %l nt_s.m fun cti on F = In t_s(x1,x2 ,n) %复化梯形求积公式 % x1,x2区间，分为n个区间 %样点矩阵及其函数值： x = lin space(x1,x2 ,n+1); y =

5、 f(x); m = len gth(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); %run 112.m clc,clear; %分为10个区间，步长0.1的积分值: F = In t_s(0,1,10); S10 = F %分为100个区间 F = In t_s(0,1,100); S100 = F %误差计算 W10 = abs(-4/9)

6、-S10); W100 = abs(-4/9)-S100); f - 可以 W = W10 W100 - %run 113.m拟合误差和步长之间的三次曲线关系 clc,clear; %建立梯形误差、辛普森误差、步长矩阵： T=zeros(1,10); S=zeros(1,10); h=zeros(1,10); fo门=1:10 - F = In t_t(0,1,10*i); fp - F = In t_s(0,1,10*i); S(i) = -4/9 - F; end 也可以明显看出辛普森误差曲线 SP = polyfit(h,S,3) %龙贝格: %Romberg.m fun cti on

7、F=Romberg(x1,x2 ,n)h(i) = 1/(10*i); L.5SGI -0.4651 -0.0276 O.DC DI SI Ou T(i) = -4/9 - F; r. 1831 -D.1502 d,OGDl TP = polyfit(h,T,3) %建立龙贝格推算矩阵、求最初步长： R = zeros(4); h = (x2-x1)/n; x = lin space(x1,x2 ,n+1); %计算矩阵第一列：复化梯形结果： for i=1:4 F = In t_t(x1,x2, n*i); R(i,1) 二 F； end %计算第二列：辛普森 for i =1:3 R(i,

8、2) = (4/3)*R(i+1,1)-(1/3)*R(i,1); end %计算第三列：复化斯科特 for i =1:2 R(i,3) = (16/15)*R(i+1,2)-(1/15)*R(i,2); end R(1,4) = (64/63)*R(2,3)-(1/63)*R(1,3); F = R(1,4); %run12.m clc,clear; F = Romberg(0,1,10); F10 = F; F = Romberg(0,1,100); F100 = F; F10,F100;-4/9-F10,-4/9-F100 ans - -0- 44C7 -0- 0037 -0. 4443

9、 -0. 0002 %右图为初始划分为10个区间和100个区间进行运算的结果， 可以看出初始 辛普森和梯形结果高出不少 10次划分的精度比 2, 我选用的是类似于复化梯形的求积方法，由于是正方形的积分区域，把它划分为 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) 随着n的增大，是会越来越逼近真实值的。程序如下： %qes2.m被积分函数： fun cti on f=qes2(x,y) f

10、 = exp(-x*y); %I nt_xy.m积分代码： function F = Int_xy(x1,x2,y1,y2,n) %仅适合此题目的正方形区域二重积分 %x1,x2,y1,y2分别为横坐标和纵坐标起点与中点； n*n块， 在其中 是的，而且存在 /4的值来近似，可知 %门 为划分为 n*n 的小正方形网格 %建立所有网格交点处的横坐标和纵坐标 : x = lin space(x1,x2 ,n+1); y = lin space(y1,y2 ,n+1); %计算每个小方格的面积， dxdy : ds = (x2-x1)*(y2-y1)/( n A2); %建立一个矩阵存放各个小方块

11、的积分值 : Z = zeros( n); %初始化积分值： F = 0; %把小方块的积分值存入矩阵 Z中，第一行第一列存放最左下角方格数据，以此类推 %并依次把小方格区域的积分值累加，得到最终积分值。 for i = 1:n for j = 1:n f = 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); end end %run2.m 运行代码： clc,clear; F = In t_xy(0,1,0,1,2); %划