MATLAB数值分析实验二（复合梯形、辛普森和龙贝格求积，以及二重

佛山科学技术学院

实 验 报 告

课程名称 数值分析 实验项目 数值积分

专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日

一、实验目的

1、理解如何在计算机上使用数值方法计算定积分?baf(x)dx的近似值;

2、学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。 3、探索二重积分??f(x,y)dxdy在矩形区域D?{(x,y)|a?x?b,c?y?d}的数值

D积分方法。

二、实验要求

（1） 按照题目要求完成实验内容； （2） 写出相应的Matlab 程序；

（3） 给出实验结果(可以用表格展示实验结果)； （4） 分析和讨论实验结果并提出可能的优化实验。 （5） 写出实验报告。

三、实验步骤

11、用不同数值方法计算积分

?04xlnxdx??

9（1）取不同的步长h，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两公式的精度。 （2）用龙贝格求积计算完成问题（1）。

2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分其中积分区域D?{0?x?1,0?y?1}。

?xye??dxdy， D

1.

%Int_t.m 复化梯形：

function F = 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 = [1 2*ones(1,m-2) 1]; %计算估计的积分值： F = h/2*sum(a.*y); %f.m

function y = f(x) y = sqrt(x).*log(x);

%run11.m clc,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 = [W10 W100]

%复化辛普森：

%Int_s.m

function F = 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.m clc,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 = [W10 W100]

%run113.m 拟合误差和步长之间的三次曲线关系。 clc,clear;

%建立梯形误差、辛普森误差、步长矩阵： T=zeros(1,10); S=zeros(1,10); h=zeros(1,10); for i=1:10

F = 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); end

TP = polyfit(h,T,3) SP = polyfit(h,S,3)

%龙贝格：

%Romberg.m：

function F=Romberg(x1,x2,n)

可以明显看出其精度高于复化梯形 也可以明显看出辛普森误差曲线各项系数都较小的 %建立龙贝格推算矩阵、求最初步长： R = zeros(4); h = (x2-x1)/n;

x = linspace(x1,x2,n+1);

%计算矩阵第一列：复化梯形结果： for i=1:4

F = Int_t(x1,x2,n*i); R(i,1) = F; end

%计算第二列：辛普森 for i =1:3

R(i,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]

%右图为初始划分为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 被积分函数： function f=qes2(x,y) f = exp(-x*y);

%Int_xy.m 积分代码：

function F = 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中，第一行第一列存放最左下角方格数据，以此类推： %并依次把小方格区域的积分值累加，得到最终积分值。 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 = 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]

%得到的积分值如右图，可见是逐渐精确的：

本文来源：https://www.bwwdw.com/article/xada.html