复化积分法(复化梯形求积,复化Simpson公式,变步长求积法)MATLAB编程实验报告

MATLAB复化梯形法与龙贝格法计算定积分复化梯形法和龙贝格法是常用的数值积分方法，用于计算定积分的近似值。

在MATLAB中，可以使用这两种方法来计算定积分。

1.复化梯形法：复化梯形法是基于将积分区间等分成若干子区间，并用梯形面积来近似每个子区间上的积分值。

整个积分的近似值等于所有子区间上的梯形面积之和。

首先，将积分区间[a,b]等分成N个子区间，每个子区间的长度为h=(b-a)/N。

然后，可以用以下公式计算每个子区间上的梯形面积：S=(f(x_i)+f(x_{i+1}))*h/2其中，f(x_i)和f(x_{i+1})是在子区间上取的两个点的函数值，x_i 和x_{i+1}分别是子区间的起始点和终止点。

```matlabh=(b-a)/N;x=a:h:b;result = 0;for i = 1:Nresult = result + (f(x(i)) + f(x(i + 1))) * h / 2;endend```2.龙贝格法：龙贝格法是一种自适应的数值积分方法，它基于逐次加密网格和不同阶数的梯形法来提高近似的精度。

首先，将积分区间[a, b]分割成n个子区间。

然后，可以使用复化梯形法来计算每个子区间上的积分值。

接下来，应用Richardson外推方法，通过逐次加密网格并对不同阶数的梯形法进行迭代，以获得更精确的近似值。

```matlabfunction result = romberg(f, a, b, max_depth, tol)R = zeros(max_depth, max_depth);h=b-a;R(1,1)=(f(a)+f(b))*h/2;for j = 2:max_depthh=h/2;R(j, 1) = R(j - 1, 1) / 2 + sum(f(a + (0:2^(j - 2)) * h) * h);for k = 2:jR(j,k)=R(j,k-1)+(R(j,k-1)-R(j-1,k-1))/(4^(k-1)-1);endif abs(R(j, j) - R(j - 1, j - 1)) < tolresult = R(j, j);return;endendresult = R(max_depth, max_depth);end```在使用上述代码计算定积分时，需要定义一个函数f来表示被积函数，并提供积分区间[a, b]的起始点和终止点。

佛山科学技术学院实验报告课程名称_______________ 数值分析________________________实验项目_______________ 数值积分____________________专业班级机械工程姓名余红杰学号2111505010 指导教师陈剑成绩日期月日一、实验目的b1、理解如何在计算机上使用数值方法计算定积分 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）分析和讨论实验结果并提出可能的优化实验。

（5）写出实验报告。

三、实验步骤1、用不同数值方法计算积xln xdx =-- 0 9（1）取不同的步长h，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两公式的精度。

（2）用龙贝格求积计算完成问题（1 ）。

2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分..e"y dxdy，其中积分区域D二｛0乞x岂1,0岂y乞1｝。

1.%lnt_t.m复化梯形：function F = Int_t(x1,x2,n)%复化梯形求积公式% x1,x2为积分起点和中点%分为n个区间，没选用步长可以防止区间数为非整数。

%样点矩阵及其函数值：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.mfun cti on y = f(x)y = sqrt(x).*log(x);%run 11.mclc,clear;%分为10个区间，步长0.1的积分值：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]%复化辛普森：%l nt_s.mfun cti on F = In t_s(x1,x2 ,n)%复化梯形求积公式% x1,x2区间，分为n个区间。

复合梯形公式、复合⾟普森公式matlab 1. ⽤1阶⾄4阶Newton-Cotes公式计算积分程序：function I = NewtonCotes(f,a,b,type)%syms t;t=findsym(sym(f));I=0;switch typecase 1,I=((b-a)/2)*(subs(sym(f),t,a)+subs(sym(f),t,b));case 2,I=((b-a)/6)*(subs(sym(f),t,a)+4*subs(sym(f),t,(a+b)/2)+...subs(sym(f),t,b));case 3,I=((b-a)/8)*(subs(sym(f),t,a)+3*subs(sym(f),t,(2*a+b)/3)+...3*subs(sym(f),t,(a+2*b)/3)+subs(sym(f),t,b));case 4,I=((b-a)/90)*(7*subs(sym(f),t,a)+...32*subs(sym(f),t,(3*a+b)/4)+...12*subs(sym(f),t,(a+b)/2)+...32*subs(sym(f),t,(a+3*b)/4)+7*subs(sym(f),t,b));case 5,I=((b-a)/288)*(19*subs(sym(f),t,a)+...75*subs(sym(f),t,(4*a+b)/5)+...50*subs(sym(f),t,(3*a+2*b)/5)+...50*subs(sym(f),t,(2*a+3*b)/5)+...75*subs(sym(f),t,(a+4*b)/5)+19*subs(sym(f),t,b));case 6,I=((b-a)/840)*(41*subs(sym(f),t,a)+...216*subs(sym(f),t,(5*a+b)/6)+...27*subs(sym(f),t,(2*a+b)/3)+...272*subs(sym(f),t,(a+b)/2)+...27*subs(sym(f),t,(a+2*b)/3)+...216*subs(sym(f),t,(a+5*b)/6)+...41*subs(sym(f),t,b));case 7,I=((b-a)/17280)*(751*subs(sym(f),t,a)+...3577*subs(sym(f),t,(6*a+b)/7)+...1323*subs(sym(f),t,(5*a+2*b)/7)+...2989*subs(sym(f),t,(3*a+4*b)/7)+...1323*subs(sym(f),t,(2*a+5*b)/7)+...3577*subs(sym(f),t,(a+6*b)/7)+751*subs(sym(f),t,b));endsyms xf=exp(-x).*sin(x);a=0;b=2*pi;I = NewtonCotes(f,a,b,1)N=1:I =N=2:I =N=3:I =(pi*((3*3^(1/2)*exp(-(2*pi)/3))/2 - (3*3^(1/2)*exp(-(4*pi)/3))/2))/4N=4:I =(pi*(32*exp(-pi/2) - 32*exp(-(3*pi)/2)))/452. 已知，因此可以通过数值积分计算的近似值。

佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 数值积分专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、理解如何在计算机上使用数值方法计算定积分⎰badxx f )(的近似值;2、学会复合梯形、复合Simpson 和龙贝格求积分公式的编程与应用。

3、探索二重积分⎰⎰Ddxdy y x f ),(在矩形区域},|),{(d y c b x a y x D ≤≤≤≤=的数值积分方法。

二、实验要求（1） 按照题目要求完成实验内容； （2） 写出相应的Matlab 程序；（3） 给出实验结果(可以用表格展示实验结果)； （4） 分析和讨论实验结果并提出可能的优化实验。

（5） 写出实验报告。

三、实验步骤1、用不同数值方法计算积分149xdx =-（1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两公式的精度。

（2）用龙贝格求积计算完成问题（1）。

2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分xyDedxdy -⎰⎰，其中积分区域{01,01}D x y =≤≤≤≤。

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.mfunction y = 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 = [W10 W100]%复化辛普森：%Int_s.mfunction F = Int_s(x1,x2,n)% 复化梯形求积公式% x1,x2 区间，分为n个区间。

陕西科技大学机械教改班用C++的积分其实积分的思想就是,微分—>求和—>取极限，如果是用纯手工法那就是先对一个函数微分，再求出它的面积，在取极限，因为我们的计算速度和计算量有限，现在有了计算机这个速度很快的机器，我们可以把微分后的每个小的面积加起来，为了满足精度，我们可以加大分区，即使实现不了微分出无限小的极限情况，我们也至少可以用有限次去接近他，下面我分析了四种不同的积分方法，和一个综合通用程序。

一．积分的基本思想1、思路：微分—>求和—>取极限。

2、Newton —Leibniz 公式 ⎰-=ba a Fb F dx x f )()()( 其中，)(x F 被积函数)(x f的原函数。

3、用计算机积分的思路在积分区间内“微分—>求和—>控制精度”。

因为计算机求和不可以取极限，也就是不可以无限次的加下去，所以要控制精度。

二．现有的理论1、一阶求积公式---梯形公式⎰=+-=b a T b f a f a b dx x f )]()([2)( 他只能精确计算被积函数为0、1次多项式时的积分。

2、二阶求积分公式——牛顿、科特斯公式 ⎰=+++-=ba Sb f a b f a f a b dx x f )]()2(4)([6)(他只能精确计算被积函数为0、1、2、3次多项式时的积分。

三．四种实现方法1.复化矩形法将积分区间[a,b]等分成n 个子区间:],[],[],[],[],[112322110n n n n x x x x x x x x x x ---、、、 则h=(b-a)/n,区间端点值k x =a+kh)hf(x ))f(x x (x I 11121=-=)()()x (22232x hf x f x I =-=............................)()()(111n ---=-=n n n n x hf x f x x I∑==ni i x hf T 1n )(源程序：#include <iostream.h>#include<math.h>double f(double x) //计算被积函数{double y;y=log(1+x)/(1+x*x); //被积函数return y;}double Tn(double a,double b,int n) //求Tn{double t=0.0;double xk; //区间端点值double t1,t2; //用来判断精度do{double h=(b-a)/n;for(int k=1;k<=n-1;k++) //每一小段的矩形叠加 {t1=t;xk=a+k*h;t+=h*f(xk);t2=t;}n++; //如果精度不够就对区间再次细分，直到达到精度要求 }while(fabs(t1-t2)<=1e-7); //判断计算精度return t;}void main(){double a=0.0; //积分下线double b=2.0; //积分上限int n=1024; //把区间分为1024段cout<<Tn(a,b,n)<<endl; //输出积分结果}执行结果：2.复化梯形法方法和复化矩形法类似，只是把原来的矩形小面积变成了梯形小面积，但是精确度明显提高了，也就是说达到同样的精度需要的时间少了。

东莞理工学院《数值分析》实验报告实验名称：牛顿插值法系别：计算机学院专业：2013级信息与计算科学班级：1班姓名：学号：实验日期：1、实验内容用不同数值方法计算积分104ln 9x xdx =-⎰。

2、算法说明 梯形求积公式算法：将积分区间[,]a b 划分为n 等份，步长b ah n-=分点为 ,1,2,...,k x a kh k n =+=。

积分10()2nk k k h I x x +==+∑。

辛普森求积公式算法：5(4)012()(4)()390xx h h f x dx y y y f c =++-⎰其中h 为步长。

3、Matlab 软件程序清单梯形求积公式TiXing_quad(a,b,h)：function t = TiXing_quad(a,b,h) %a 为积分下界，b 为积分上界，h 为步长。

format long x = a:h:b;y = sqrt(x).*log(x); y(1) = 0;t = 0;for k=1:(b-a)/h,t=t+y(k)+y(k+1);endt=t*h/2;辛普森求积公式Sinpson_quad(a,b,h)：function s=Sinpson_quad(a,b,h)format longx=a:h:b;y=sqrt(x).*log(x);z=sqrt(x+h/2).*log(x+h/2);y(1)=0;s=0;for k=1:(b-a)/h,s=s+y(k)+y(k+1)+4*z(k);ends=s*h./6;4、运行结果真值I=-4/9=-0.444444444444444次数I（梯形公式）I（辛普森公式）50 -0.441090226387332 -0.443793798301150100 -0.443117905322695 -0.444194********* 200 -0.443925359444891 -0.4443490454521025、分析与思考“”“”At the end, Xiao Bian gives you a passage. Minand once said, "people who learn to learn are very happy people.". In every wonderful life, learning is an eternal theme. As a professional clerical and teaching position, I understand the importance of continuous learning, "life is diligent, nothing can be gained", only continuous learning can achieve better self. Only by constantly learning and mastering the latest relevant knowledge, can employees from all walks of life keep up with the pace of enterprise development and innovate to meet the needs of the market. This document is also edited by my studio professionals, there may be errors in the document, if there are errors, please correct, thank you!。

佛山科学技术学院实 验 报 告课程名称 数值分析实验项目 数值积分专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、理解如何在计算机上使用数值方法计算定积分⎰ba dx x f )(的近似值;2、学会复合梯形、复合Simpson 和龙贝格求积分公式的编程与应用。

3、探索二重积分⎰⎰Ddxdy y x f ),(在矩形区域},|),{(d y c b x a y x D ≤≤≤≤=的数值积分方法。

二、实验要求（1） 按照题目要求完成实验内容；（2） 写出相应的Matlab 程序；（3） 给出实验结果(可以用表格展示实验结果)；（4） 分析和讨论实验结果并提出可能的优化实验。

（5） 写出实验报告。

三、实验步骤1、用不同数值方法计算积分1049xdx =- （1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两公式的精度。

（2）用龙贝格求积计算完成问题（1）。

2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分xy De dxdy -⎰⎰， 其中积分区域{01,01}D x y =≤≤≤≤。

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.mfunction y = 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 = [W10 W100]%复化辛普森：%Int_s.mfunction F = Int_s(x1,x2,n)% 复化梯形求积公式% x1,x2 区间，分为n个区间。

变步长梯形求积matlab
变步长梯形求积法是一种数值积分的方法，用于计算一个函数在给定区间上的定积分。

在MATLAB中，可以使用以下代码来实现变步长梯形求积法：
```matlab
function integral = trapezoidal_rule(a, b, n, f)
h = (b - a) / n;
x = a:h:b;
y = f(x);
integral = (h/2) * (y(1) + 2*sum(y(2:end-1)) + y(end)); end
```
其中，`a`和`b`是积分区间的上下限，`n`是划分的梯形数目，`f`是要求积的函数。

函数内部的变量`h`计算出每个梯形的宽度，`x`生成了在区间上等间隔的节点，`y`根据`f`函数计算了这些节点上的函数值。

最后，利用梯形公式，将每个梯形面积加总并乘以梯形宽度的一半，得到了定积分的近似值。

使用时，可以将要求积的函数定义为一个匿名函数，并将其作为参数传递给`trapezoidal_rule`函数。

例如，计算函数`f(x) = x^2`在区间`[0, 1]`上的定积分，可以执行以下代码：
```matlab
f = @(x) x.^2;
a = 0;
b = 1;
n = 100;
integral = trapezoidal_rule(a, b, n, f);
disp(integral);
```
上述代码将会输出`0.3333`，这是函数`f(x) = x^2`在区间`[0, 1]`上的定积分的近似值。

数值积分 1.求积分dx e x ⎰--112，在积分区间中，点与点之间的间隔取为0.1.解：（一）用MATBLE 编写复化梯形求积函数：function I=T_quad(x,y)n=length(x);m=length(y);if n ~=merrorendh=(x(n)-x(1))/(n-1);a=[1 2*ones(1,n-2) 1];I =h/2*sum(a.*y);输入：x=-1:0.1:1;y=exp(-x.^2);I =T_quad(x,y)运行得到：I =1.4924(二)用MATBL 编写复化Simpson 求积函数：function I=S_quad(x,y);n=length(x);m=length(y);if n ~=merrorendif rem(n-1,2) ~=0I=T_quad(x,y);return;endN=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n);for k=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;endI=h/6*sum(a.*y);输入：x= -1:0.1:1;y=exp(-x.^2);I= S_quad(x,y)运行得到：I =1.4936（三）用MATBL 编写复化Cotes 求积函数：function I=C_quad(x,y);n=length(x);m=length(y);if n ~=merrorendif rem(n-1,4) ~=0I=S_quad(x,y);returnendN=(n-1)/4;h=(x(n)-x(1))/N;a=zeros(1,n);for k=1:Na(4*k-3)=a(4*k-3)+7;a(4*k-2)=a(4*k-2)+32;a(4*k-1)=a(4*k-1)+12;a(4*k)=a(4*k)+32;a(4*k+1)=a(4*k+1)+7;endI=h/90*sum(a.*y);输入：x= -1:0.1:1;y=exp(-x.^2);I= C_quad(x,y)运行得：I =1.4936（四）利用trapz()函数，采用复化梯形公式求积分输入：x=-1:0.1:1;y=exp(-x.^2);I=trapz(x,y)输出：I =1.49242. 求积分dx e x ⎰--112，取精度要求510-=ε （一）用MATLAB 编写自适应步长的梯形公式function I=T_quad_iter(fun,a,b,ep)if nargin<4 ep=1e-5;endN=1;h=b-aT=h/2*(feval(fun,a)+feval(fun,b));while 1h=h/2;I=T/2;for k=1:NI=I+h*feval(fun,a+(2*k-1)*h);endif abs(I-T)<ep breakendN=2*N;T=IEnd输入：I=T_quad_iter(@(x)exp(-x.^2),-1,1)输出：I =1.4936（二）用MATLAB编写自适应步长的Simpson公式function I=S_quad_iter(fun,a,b,ep)if nargin<4 ep=1e-5endN=1;h=b-aT1=h/2*(feval(fun,a)+feval(fun,b));S0=T1;while 1h=h/2;T2=T1/2;for k=1:NT2=T2+h*feval(fun,a+(2*k-1)*h)endI=(4*T2-T1)/3;if abs(I-S0)<ep breakendN=2*N;T1=T2;S0=IEnd输入：I=S_quad_iter(@(x) exp(-x.^2),-1,-1)输出：I =1.4936（三）用quad（）函数，采用自适应步长的Simposon求积分输入：I=quad(@(x) exp(-x.^2) ,-1,1)输出：I =1.4936。

matlab编程积分复合辛普森公式编程求解复合辛普森公式是一种常用的数值积分方法，在MATLAB中可以通过编写相应的代码实现。

本文将详细介绍如何使用MATLAB编程计算复合辛普森公式，并给出具体的代码实现。

我们需要了解什么是复合辛普森公式。

复合辛普森公式是一种数值积分方法，用于近似计算函数的定积分。

它是在区间[a, b]上使用多个小区间进行逼近，而不是直接在整个区间上进行逼近。

这种方法的优势在于可以提高计算精度，并且对于复杂的函数也能够得到较好的近似结果。

我们需要将整个区间[a, b]划分为n个小区间。

每个小区间的长度为h=(b-a)/n。

然后，我们可以使用复合辛普森公式来近似计算每个小区间上的定积分。

复合辛普森公式的表达式为：I = (h/6)*(f(a) + 4*f((a+b)/2) + f(b))其中，f(x)是要求解的函数。

根据复合辛普森公式的定义，我们需要对每个小区间应用该公式进行求解，并将结果累加得到最终的积分值。

在MATLAB中，我们可以通过编写以下代码来实现复合辛普森公式的求解：```matlabfunction I = composite_simpson(f, a, b, n)h = (b-a)/n;x = a:h:b;y = f(x);I = 0;for i = 1:nI = I + h/6*(y(i) + 4*y(i+1) + y(i+2));endend```上述代码中，函数composite_simpson接受四个参数：函数f、积分区间的起点a、终点b和划分的小区间数n。

其中，函数f是一个函数句柄，表示要求解的函数。

在代码中，我们首先计算出每个小区间的长度h，并生成对应的x值。

然后，通过调用函数f计算出对应的y值。

接下来，我们使用循环对每个小区间应用复合辛普森公式，并将结果累加到变量I中。

最后，我们将得到的积分值I作为函数的输出。

在使用该函数时，我们需要先定义要求解的函数，并将其作为参数传递给composite_simpson函数。