龙贝格公式和辛普森公式和复合梯形公式

佛山科学技术学院实验报告课程名称_______________ 数值分析________________________实验项目_______________ 数值积分____________________专业班级机械工程姓名余红杰学号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个区间。

数值积分：梯形规则数值积分：梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式1.问题描述微积分方法求积有很大的局限性，当碰到被积函数很复杂时，找不到相应的原函数。

积分值在几何上可解释为由 x=a,x=b,y=0和y=f(x) 所围成的曲边梯形的面积。

积分计算之所以有困难，就是因为这个曲边梯形有一条边y=f(x)是曲线。

2.理论与方法依据积分中值定理，底为b-a,而高为f(e)的矩形面积恰等于所求曲边梯形的面积I.f(e)称作区间[a,b]上的平均高度。

这样，只要对平均高度f(e)提供一种算法，便相应地获得一种数值求积的算法。

1.梯形规则(Trapezoidal rule)简单选取区间[a ,b]的中点高度作为平均高度。

取h=b-aa0=⌠(a-b)(x-b)/(a-b)dx=(b-a)/2a1=⌠(a-b)(x-a)/(b-a)dx=(b-a)/2得到：2.辛普森规则(Simpson rule)可视作用a , b与c=(a+b)/2三点高度的加权平均值作为平均高度。

3.复合梯形规则(Composite numerical)设将求积区间[a,b]划分为n等份，步长h=(b-a)/2 ,等分点为xi=a+bi , i=0,1,...,n 所谓复化求积法，就是先用低阶求积公式求得每个子段[xi,xi+1]上的积分值,然后再将它们累加求和，用各段积分之和Ii,i=0,1,n-1作为所求积分的近似值。

复化梯形公式：4.复合辛普森规则(Composite Simpson)记子段[xi,xi+1]的中点为则复化公式为复化Simpson公式：5.龙贝格求积公式(Romberg)龙贝格求积公式也称为逐次分半加速法。

它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。

作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。

实验八数值积分信息与计算科学金融崔振威201002034031一、实验目的：1、掌握数据积分算法设计及程序实现二、实验内容：1、p290-1、p301-2三、实验要求：主程序：复合梯形公式：function [I,step,h2] = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;h2=(b-a)/n;function [I,step,h] = IntSimpson(f,a,b,type,eps)%type = 1 辛普森公式%type = 2 辛普森3/8公式%type = 3 复合辛普森公式if(type==3 && nargin==4)eps=1.0e-4; %精度为0.0001endI=0;switch typecase 1,I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+...4*subs(sym(f),findsym(sym(f)),(a+b)/2)+...subs(sym(f),findsym(sym(f)),b));step=1;case 2,I=((b-a)/8)*(subs(sym(f),findsym(sym(f)),a)+...3*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+ ...3*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+subs(sym(f),findsym(sym(f)),b));step=1;case 3,n=2;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+...4*subs(sym(f),findsym(sym(f)),(x+x1)/2)+...subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;endfunction [q,step]=Roberg(f,a,b,eps)%f是被积函数%a积分上限%b积分下限%eps是精度%q输出结果%step循环次数if(nargin==3)eps=1.0e-4;end;M=1;tol=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b)); while tol>epsk=k+1;h=h/2;Q=0;for i=1:Mx=a+h*(2*i-1);Q=Q+subs(sym(f),findsym(sym(f)),x);end;T(k+1,1)=T(k,1)/2+h*Q;M=M*2;for j=1:k;T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);end;tol=abs(T(k+1,j+1)-T(k,j));endq=T(k+1,k+1);step=k;p290-11、（i ）用组合梯形公式和M=10求下列每个积分(a)、dx x 1112)1(--⎰+解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('(1+x^2)^(-1)',-1,1)得出结果：q =1.56996299445358s =20h =0.10000000000000 所以dx x 1112)1(--⎰+值约为1.569962994，步长为0.1，M 为20(b)、dx x ⎰+10))2sin(2(解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('(2+sin(2*x^(1/2)))^(-1)',0,1)得出结果：q =0.35117779429220s =19h =0.05263157894737 所以dx x ⎰+10))2sin(2(值约为 0.351177794，步长为0.052631579，M 为19(c)、⎰425.0/x dx解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('1/(x^(1/2))',0.25,4)得出结果：q =3.00216646875717s =46h =0.08152173913043 所以⎰425.0/x dx 值约为 3.002166469，步长为0.081521739，M 为46(d)、dx e x x -⎰402解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('x^2*exp(-1)',0,4)得出结果：q =7.85012162896417s =44h =0.09090909090909所以dx e x x -⎰402值约为7.850121629，步长为0.090909090，M 为44(e)、⎰20)(2dx x xcox解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('2*x*cos(x)',0,2)得出结果：q =0.80323191187607s =36h =0.05555555555556所以⎰20)(2dx x xcox 值约为 0.803231912，步长为0.055555556，M 为36(f)、dx e x x ⎰-π0)2sin(解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('sin(2*x)*exp(-x)',0,pi)得出结果：q =-6.738211157945265e-018s =1h =3.14159265358979所以dx e x x ⎰-π0)2sin(值约为 -6.7382111578，步长为3.141592653，M 为12、（ii ）用组合辛普森公式和M=10求下列每个积分(a)、dx x 1112)1(--⎰+解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('(1+x^2)^(-1)',-1,1,3)得出结果：I =1.57079538809119step =5s =0.400000000000000所以dx x 1112)1(--⎰+值约为1.570795388，步长为0.4，M 为5(b)、dx x ⎰+10))2sin(2(解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('(2+sin(2*x^(1/2)))^(-1)',0,1,3)得出结果：I =0.35047636367125step =10s =0.10000000000000 所以dx x ⎰+10))2sin(2(值约为0.350476364，步长为0.1，M 为10(c)、⎰425.0/x dx解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('1/(x^(1/2))',0.25,4,3)得出结果：I =3.00030161515673step =14s =0.26785714285714 所以⎰425.0/x dx 值约为3.000301615，步长为0.267857143，M 为14(d)、dx e x x -⎰402解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('x^2*exp(-1)',0,4,3)得出结果：I =7.84809474499077step =4s =1所以dx e x x -⎰402值约为7.848094745，步长为1，M 为4(e)、⎰20)(2dx x xcox解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('2*x*cos(x)',0,2,3)得出结果：I =0.80494830253761step =6s =0.33333333333333所以⎰20)(2dx x xcox 值约为0.804948303，步长为0.333333333，M 为6(f)、dx e x x ⎰-π0)2sin(解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('sin(2*x)*exp(-x)',0,pi,3)得出结果：I =-6.738211157945265e-018step =2s =1.57079632679490所以dx e x x ⎰-π0)2sin(值约为 -6.738211158，步长为1.570796327，M 为2P301-12、（a ）、dx x x ⎰-2024解：在matlab 窗口中输入>> [q,s]=Roberg('sqrt(4*x-x^2)',0,2)得出结果：>> [q,s]=Roberg('sqrt(4*x-x^2)',0,2)q =3.14155917503659s =9所以，该积分值约为3.1415591750（b ）、dx x ⎰+1021/4解：在matlab 窗口中输入>> [q,s]=Roberg('4/(1+x^2)',0,1)得出结果：q =3.14159266527772s =4所以，该积分值约为3.1415926653通过利用龙贝格公式积分求得结果，两个定积分得出结果的积分速度明显不同。

（一） 实验目的熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化simpson 公式以及romberg 积分。

（二） 问题描述问题三数值积分椭圆周长的计算。

考虑椭圆22221x y a b+=,为计算其周长,只要计算其第一象限的长度即可.用参数方程可以表示为cos (0/2)sin x a t t y b t π=⎧≤≤⎨=⎩,计算公式为/0π⎰为计算方便,我们可以令1a =,即计算下面的积分/0π⎰/0π=⎰(/0π⎰/0a π=⎰可以归结为上面的形式)采用复化梯形公式，复化Simpson 公式以及Romberg 积分的方法计算积分/0()I b π=⎰给出通用程序，该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分。

程序输出为计算出的数值积分值以及计算函数值的次数。

（三） 算法介绍首先利用给出的各迭代公式，设计程序。

在matlab 对话框中输入要计算的函数，给出区间和精度。

复化梯形的迭代公式为：;复化simpson迭代公式为：;Romberg迭代公式为：。

（四）程序对于复化梯形公式和复化simpson公式，我们放在中。

（%标记后的程序可用来把b看为变量时的算法实现）%复化梯形公式function y=jifenn(f,n,a,b) （说明：f表示任一函数，n精度，a，b为区间）fi=f(a)+f(b);h=(b-a)/n;d=1;%function f=jifen(n,a,b,c)%syms t%y=sqrt(1+(c^2-1)*cos(t)^2);%ya=subs(y,t,a);%yb=subs(y,t,b);%fi=ya+yb;for i=1:n-1x=a+i*h;fi=fi+2*f(x);d=d+1;%yx=subs(y,t,x);%fi=fi+2*yx;endf4=h/2*fi,d%复化simposon公式f1=0;f2=0;dd=1;for i=1:n-1dd=dd+1;if rem(i,2)~=0;x1=a+i*h;f1=f1+f(x1);else rem(i,2)==0;x2=a+i*h;f2=f2+f(x2) ;endendf3=(h/3)*(f(a)+4*f1+2*f2+f(b)),dd对于romberg积分，建立文件。