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

实验八 数值积分信息与计算科学金融 崔振威 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)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=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 type case 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)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=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); end end I=I2; step=n;end龙贝格公式function 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 toleps k=k+1; h=h/2; Q=0; for i=1:M x=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)/(4j-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) 、解：在matlab窗口中输入 q,s,h=CombineTraprl(1+x2)(-1),-1,1) 得出结果：q = 1.56996299445358s = 20h = 0.10000000000000所以值约为1.569962994，步长为0.1，M为20(b) 、解：在matlab窗口中输入 q,s,h=CombineTraprl(2+sin(2*x(1/2)(-1),0,1) 得出结果：q = 0.35117779429220s = 19h = 0.05263157894737所以值约为 0.351177794，步长为0.052631579，M为19(c) 、解：在matlab窗口中输入 q,s,h=CombineTraprl(1/(x(1/2),0.25,4) 得出结果：q = 3.00216646875717s = 46h = 0.08152173913043所以值约为 3.002166469，步长为0.081521739，M为46(d) 、解：在matlab窗口中输入 q,s,h=CombineTraprl(x2*exp(-1),0,4) 得出结果：q = 7.85012162896417s = 44h = 0.09090909090909所以值约为7.850121629，步长为0.090909090，M为44(e) 、解：在matlab窗口中输入 q,s,h=CombineTraprl(2*x*cos(x),0,2) 得出结果：q = 0.80323191187607s = 36h = 0.05555555555556所以值约为 0.803231912，步长为0.055555556，M为36(f) 、解：在matlab窗口中输入 q,s,h=CombineTraprl(sin(2*x)*exp(-x),0,pi) 得出结果：q = -6.738211157945265e-018s = 1h = 3.14159265358979所以值约为 -6.7382111578，步长为3.141592653，M为12、 （ii）用组合辛普森公式和M=10求下列每个积分(a) 、解：在matlab窗口中输入 I,step,s = IntSimpson(1+x2)(-1),-1,1,3) 得出结果：I = 1.57079538809119step = 5s = 0.400000000000000所以值约为1.570795388，步长为0.4，M为5(b)、解：在matlab窗口中输入 I,step,s = IntSimpson(2+sin(2*x(1/2)(-1),0,1,3) 得出结果：I = 0.35047636367125step = 10s = 0.10000000000000所以值约为0.350476364，步长为0.1，M为10(c)、解：在matlab窗口中输入 I,step,s = IntSimpson(1/(x(1/2),0.25,4,3) 得出结果：I = 3.00030161515673step = 14s = 0.26785714285714所以值约为3.000301615，步长为0.267857143，M为14(d)、解：在matlab窗口中输入 I,step,s = IntSimpson(x2*exp(-1),0,4,3) 得出结果：I = 7.84809474499077step = 4s = 1所以值约为7.848094745，步长为1，M为4(e)、解：在matlab窗口中输入 I,step,s = IntSimpson(2*x*cos(x),0,2,3) 得出结果：I = 0.80494830253761step = 6s = 0.33333333333333所以值约为0.804948303，步长为0.333333333，M为6(f)、解：在matlab窗口中输入 I,step,s = IntSimpson(sin(2*x)*exp(-x),0,pi,3) 得出结果：I = -6.738211157945265e-018step = 2s = 1.57079632679490所以值约为 -6.738211158，步长为1.570796327，M为2P301-12、（a） 、解：在matlab窗口中输入 q,s=Roberg(sqrt(4*x-x2),0,2)得出结果： q,s=Roberg(sqrt(4*x-x2),0,2)q = 3.14155917503659s = 9所以，该积分值约为3.1415591