下载文档原格式
佛山科学技术学院实验报告课程名称_______________ 数值分析________________________实验项目_______________ 数值积分____________________专业班级机械工程姓名余红杰学号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积分，建立文件。
近似求积公式在数学中，求解一个函数的积分是一项重要的任务。
但是，对于某些函数，我们可能无法找到其精确的积分，或者即使找到了，也不方便使用。
这时，我们可以使用近似求积公式来估算函数的积分值。
近似求积公式是一种数值积分方法，其基本思想是将要积分的函数在一定区间上近似为某个简单函数的和或积，然后计算这个简单函数的积分值，从而得到原函数的近似积分值。
这种方法常常用于数值计算和科学工程中。
常见的近似求积公式有梯形公式、辛普森公式和龙贝格公式等。
下面我们将分别介绍这些公式的原理和应用。
一、梯形公式梯形公式是最简单的近似求积公式之一，其基本思想是将要积分的函数在积分区间上近似为一个梯形，然后计算梯形面积得到近似积分值。
具体来说，梯形公式的计算公式为：∫a~bf(x)dx ≈ (b-a)×[f(a)+f(b)]/2其中，f(x)是要积分的函数，a和b分别是积分区间的下限和上限。
梯形公式的误差随着积分区间的缩小而减小，但随着步长的增加而增大。
因此，在实际应用中，我们需要根据所需的精度和计算资源，选择合适的步长来进行计算。
二、辛普森公式辛普森公式是一种二次近似求积公式，其基本思想是将要积分的函数在积分区间上近似为一个二次函数，然后计算二次函数的积分值得到近似积分值。
具体来说，辛普森公式的计算公式为：∫a~bf(x)dx ≈ (b-a)×[f(a)+4f((a+b)/2)+f(b)]/6其中，f(x)是要积分的函数，a和b分别是积分区间的下限和上限。
辛普森公式的误差随着积分区间的缩小而减小，但随着步长的增加而增大。
因此，在实际应用中，我们需要根据所需的精度和计算资源，选择合适的步长来进行计算。
三、龙贝格公式龙贝格公式是一种多项式近似求积公式，其基本思想是将要积分的函数在积分区间上近似为一个多项式函数，然后计算多项式函数的积分值得到近似积分值。
具体来说，龙贝格公式的计算公式为：B(m,n) = (4^nB(m,n-1)-B(m-1,n-1))/(4^n-1)其中，m和n分别表示递归计算的次数和精度级别，B(m,n)表示第m次递归、精度为n的近似积分值。
数值计算方法上机题目3一、计算定积分的近似值：221x e xe dx =⎰ 要求：（1）若用复化梯形公式和复化Simpson 公式计算，要求误差限71021-⨯=ε，分别利用他们的余项估计对每种算法做出步长的事前估计；（2）分别利用复化梯形公式和复化Simpson 公式计算定积分；（3）将计算结果与精确解比较，并比较两种算法的计算量。
1.复化梯形公式程序：程序1（求f （x ）的n 阶导数：syms xf=x*exp(x) %定义函数f （x ）n=input('输入所求导数阶数:')f2=diff(f,x,n) %求f(x)的n 阶导数结果1输入n=2f2 =2*exp(x) + x*exp(x)程序2：clcclearsyms x%定义自变量xf=inline('x*exp(x)','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('(2*exp(x) + x*exp(x))','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。
f3='-(2*exp(x) + x*exp(x))'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=5*10^(-8) %精度要求值a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求连加和xk=a+k*hTn1=Tn1+f(xk)endTn=h/2*((f(a)+2*Tn1+f(b)))z=exp(2)R=Tn-z %求已知值与计算值的差fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn)fprintf('等分数 n=')disp(n) %输出等分数fprintf('已知值与计算值的误差 R=')disp(R)输出结果显示：用复化梯形算法计算的结果Tn= 7.3891等分数n=7019已知值与计算值的误差R= 2.8300e-0082. Simpson公式程序：程序1：（求f（x）的n阶导数）：syms xf=x*exp(x) %定义函数f（x）n=input('输入所求导数阶数:')f2=diff(f,x,n) %求f(x)的n阶导数结果1输入n=4f2 =4*exp(x) + x*exp(x)程序2：clcclearsyms x%定义自变量xf=inline('x*exp(x)','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('(4*exp(x) + x*exp(x))','x') %定义f(x)的四阶导数，输入程序1里求出的f2即可f3='-(4*exp(x) + x*exp(x))'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10^(-8) %精度要求值a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endendh=(b-a)/n %求hSn1=0Sn2=0for k=0:n-1 %求两组连加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b)) %因Sn2多加了k=0时的值，故减去f(a)z=exp(2)R=Sn-z %求已知值与计算值的差fprintf('用Simpson公式计算的结果 Sn=')disp(Sn)fprintf('等分数 n=')disp(n)fprintf('已知值与计算值的误差 R=')disp(R)输出结果显示：用Simpson公式计算的结果Sn= 7.3891等分数n=24已知值与计算值的误差R= 2.7284e-008用复化梯形公式计算的结果为：7.3891，与精确解的误差为：2.8300e-008。
湖南农业大学综合设计报告综合设计五多方法求解数值积分学生姓名：学号：年级专业：指导老师：学院：评阅成绩：评阅意见：成绩评定教师签名：时间：湖南·长沙提交日期：2014年6月多方法求解数值积分具体题目要求：用不同数值方法计算积分49xdx=-⎰(1) 取不同的步长h，分别用复合梯形及复合辛普森公式计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？(2) 用龙贝格求积计算完成问题(1);(3) 用自适应辛普森积分，使其精度达到410-。
1设计目的、要求由积分学基本理论，定积分可由Newton Leibniz-公式计算，但是对于一些无法找到原函数的函数（如2xe-等）不能通过牛顿—莱布尼兹公式计算，就必须得另寻它法。
因此需要我们能够熟练地应用常用的数值积分计算方法（如机械求积、Newton Cotes-公式等）并掌握结合数值计算软件（Matlab、Lingo 等）及计算机高级语言()c java、进行对应算法实现的技能。
熟练数学软件求解数学问题，掌握各种数学问题的求解方法。
本设计主要是通过多种复合求积公式求解积分，主要包括复化梯度法、复化辛普森法、龙贝格以及自适应辛普森法等求解方法，利用Matlab软件编写相对应的算法进行求解，大大地提高了解题的速度。
2设计原理由积分中值定理我们可以知道在积分区间[],a b内存在一点ξ，使得式子()()()baf x dx b a fξ=-⎰成立。
这个式子在于对于点ς的具体位置一般是不知道的，因此难以准确算出()fξ的值。
也就是不同算法求得平均高度()fξ，对应的就是一种不同的数值求积方法。
更一般地，我们可以在区间[],a b上适当选取某些节点k x，然后用的加权平均得到平均高度()fξ的近似值，这样构造出的求积公式具有下列形式：()()nbk k ak f x dx A f x =≈∑⎰称为机械求积公式。
龙贝格算法matlab程序一、龙贝格算法简介龙贝格算法是一种数值积分的方法，它可以用来计算函数在给定区间上的定积分。
该算法基于复合梯形公式和复合辛普森公式，通过逐步逼近真实值来得到数值解。
它是一种自适应方法，即在每个子区间上都采用不同的步长以获得更高的精度。
二、龙贝格算法的原理1. 复合梯形公式复合梯形公式是将一个区间分成若干个小区间，在每个小区间上应用梯形公式求出积分值，最后将所有小区间的积分值相加得到整个区间上的积分值。
具体公式如下：$$\int_{a}^{b}f(x)dx \approx \frac{h}{2}(f(a)+2\sum_{i=1}^{n-1}f(a+ih)+f(b))$$其中，h为步长，n为子区间数。
2. 复合辛普森公式复合辛普森公式是将一个区间分成若干个小区间，在每个小区间上应用辛普森公式求出积分值，最后将所有小区间的积分值相加得到整个区间上的积分值。
具体公式如下：$$\int_{a}^{b}f(x)dx \approx \frac{h}{6}(f(a)+4\sum_{i=1}^{n/2-1}f(a+(2i)h)+2\sum_{i=1}^{n/2}f(a+(2i-1)h)+f(b))$$其中，h为步长，n为子区间数。
3. 龙贝格算法龙贝格算法是通过不断加密网格来逼近真实值的方法。
首先，将整个区间分成若干个小区间，并在每个小区间上应用复合梯形公式求出初始积分值T(0,0)，然后将这些积分值相加得到整个区间上的积分值T(0,1)。
接着，在每个小区间上应用复合辛普森公式求出更精确的积分值T(1,0)，并将所有小区间的积分值相加得到整个区间上的积分值T(1,1)。
然后，计算两次结果之差E(1,0)=T(1,1)-T(0,1)，如果E(1,0)小于给定误差限，则直接输出T(1,1)作为积分结果；否则，在每个子区间中再次应用复合辛普森公式求出更精确的积分值，并计算两次结果之差，直到满足误差限为止。
实验八数值积分信息与计算科学金融崔振威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通过利用龙贝格公式积分求得结果，两个定积分得出结果的积分速度明显不同。
