数值积分的MATLAB实现.doc

数值分析实验报告实验名称使用matlab编写数值计算程序实验时间*姓名*班级*学号*成绩实验报告内容要求：一、实验目的与内容；二、算法描述(数学原理或设计思路、计算公式、计算步骤)；三、程序代码；四、数值结果；五、计算结果分析（如初值对结果的影响；不同方法的比较；该方法的特点和改进等）；六、实验中出现的问题，解决方法及体会（整个实验过程中（包括程序编写，上机调试等）出现的问题及其处理等广泛的问题）. 实验四 数值积分的Matlab实现一、实验目的与要求1熟练梯形公式、Simpson公式、复化梯形公式、复化Simpson公式、Romberg求积公式；2熟悉符号积分，能数值积分的值与精确值的比较;3培养编程与上机调试能力。二、实验原理设将积分区间a,b划分为n等份，步长h=，选取等距节点x=a+kh构造出的插值型求积公式 I=（b-a）称为牛顿-科特斯公式，式中C称为科特斯系数.当n=1时，C= C=，即为梯形公式T=f(a)+f(b)当n=2时，同理可求Simpson公式：S=.对应的复合公式分别为T，S，运用理查德外推加速方法可得龙贝格求积算法：T，k=1,2,. 三、实验内容与步骤分别用梯形公式、Simpson公式、复化梯形公式、复化Simpson公式、Romberg求积公式求下列定积分,观察所得结果,并与精确值进行比较.(4)(5)1、 依照实验原理编写各种公式的的程序.2、 首先在电脑上安装matlab，然后启动matlab，分别建立不同公式的M文件，实验程序如下； 梯形公式程序代码 程序代码说明function T=chen_trap(fun,a,b)；T=(b-a)*(feval(fun,a)+feval(fun,b)/2；%递推梯形公式主程序%fun:被积函数%a,b分别为积分区间的左右端点 Simpson公式程序代码 程序代码说明function y=chen_simpson(fun,a,b)；c=(a+b)/2;y=(b-a)/6*(feval(fun,a)+4*feval(fun,c)+feval(fun,b);%simpson公式主程序%fun表示被积函数句柄%a,b分别为积分区间的左右端点复合梯形公式程序代码 程序代码说明function y=comtrap(fun,a,b,n)；z1=feval(fun,a)+feval(fun,b);h=(b-a)/n; z2=0; x=a;for k=1:n-1 x1=x+k*h; z2= z2+2*feval(fun,x1); endy=(z1+z2)*h/2;%复化simpson公式主程序%fun表示被积函数句柄%a,b分别为积分区间的左右端点%n为等分积分区间的个数%y为积分近似值Sn 复合Simpson公式程序代码 程序代码说明 function y=comsimpson(fun,a,b,n)； z1=feval(fun,a)+feval(fun,b);m=n/2;h=(b-a)/(2*m); x=a;z2=0; z3=0; x2=0; x3=0;for k=1:m-1 x2=x+2*k*h; z2=z2+2*feval(fun,x2); endfor k=1:m x3=x+(2*k-1)*h; z3= z3+4*feval(fun,x3); endy=(z1+z2+z3)*h/3;%复化simpson公式主程序%fun表示被积函数句柄%a,b分别为积分区间的左右端点%n=2m,m为等分积分区间的个数%y为积分近似值Sm Romberg求积公式程序代码 程序代码说明 function quad,R=Romberg(f,a,b,eps)； h=b-a;R(1,1)=h*(feval(f,a)+feval(f,b)/2;M=1; J=0; err=1;while Jeps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; M=2*M; for k=1:J R(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J)-R(J+1,J+1);endquad=R(J+1,J+1);%f表示被积函数句柄%a,b表示被积区间a,b的端点%eps表示精度%quad是用Romberg加速算法求得的积分值%R为Romberg表%err表示误差的估计3、 具体函数编程如下：（1）； 首先编写被积函数的M文件f1.mfunction y=f1(x) 然后编写如下程序，并保存为t1.m%t1.mclear;a=0;b=1;T=chen_trap(f1,a,b);S=chen_simpson(f1,a,b);FT8=comtrap(f1,a,b,8);FS4=comsimpson(f1,a,b,8);quad=Romberg(f1,a,b,10e-6);syms x f=x2*sqrt(2*x2+3);S1=int(f, x, 0,1); format long, Sj=double(S1);disp(积分准确值为:),Sjdisp(积分近似值分别为:),T,S,FT8,FS4,quaddisp(与精确值比较的误差为:),wuchi=double(abs(Sj-T),double(abs(Sj-S),double(abs(Sj-FT8),double(abs(Sj-FS4),double(abs(Sj-quad);首先编写被积函数的M文件f2.mfunction y=f2(x)y=x*sin(x);然后编写如下程序,并保存为t2.m%t2.mclear;a=0;b=2*pi;T=chen_trap(f2,a,b);S=chen_simpson(f2,a,b);FT8=comtrap(f2,a,b,8);FS4=comsimpson(f2,a,b,8);quad=Romberg(f2,a,b,10e-6);syms x f=x*sin(x);S1=int(f, x, 0, 2*pi); format long, Sj=double(S1);disp(积分准确值为:),Sjdisp(积分近似值分别为:),T,S,FT8,FS4,quaddisp(与精确值比较的误差为:),wuchi=double(abs(Sj-T),double(abs(Sj-S),double(abs(Sj-FT8),double(abs(Sj-FS4),double(abs(Sj-quad) 首先编写被积函数的M文件f3.mfunction y=f3(x)y=x*sqrt(1+x2);然后编写如下程序,并保存为t3.m%t3.mclear;a=0;b=3;T=chen_trap(f3,a,b);S=chen_simpson(f3,a,b);FT8=comtrap(f3,a,b,8);FS4=comsimpson(f3,a,b,8);quad=Romberg(f3,a,b,10e-6);syms x f= x*sqrt(1+x2);S1=int(f, x, 0,3); format long, Sj=double(S1);disp(积分准确值为:),Sjdisp(积分近似值分别为:),T,S,FT8,FS4,quaddisp(与精确值比较的误差为:),wuchi=double(abs(Sj-T),double(abs(Sj-S),double(abs(Sj-FT8),double(abs(Sj-FS4),double(abs(Sj-quad) (4); 首先编写被积函数的M文件f4.mfunction y=f4(x)y=4(1+x2);然后编写如下程序,并保存为t4.m%t4.mclear;a=0;b=1;T=chen_trap(f4,a,b);S=chen_simpson(f4,a,b);FT8=comtrap(f4,a,b,8);FS4=comsimpson(f4,a,b,8);quad=Romberg(f4,a,b,10e-6);syms x f= 4(1+x2);S1=int(f, x, 0,1); format long, Sj=double(S1);disp(积分准确值为:),Sjdisp(积分近似值分别为:),T,S,FT8,FS4,quaddisp(与精确值比较的误差为:),wuchi=double(abs(Sj-T),double(abs(Sj-S),double(abs(Sj-FT8),double(abs(Sj-FS4),double(abs(Sj-quad) (5)首先编写被积函数的M文件f5.mfunction y=f5(x)y=sqrt1+(cosx)2;然后编写如下程序,并保存为t5.m%t5.mclear;a=0;b=48;T=chen_trap(f5,a,b);S=chen_simpson(f5,a,b);FT8=comtrap(f5,a,b,8);FS4=comsimpson(f5,a,b,8);quad=Romberg(f5,a,b,10e-6);syms x f= sqrt1+(cosx)2;S1=int(f, x, 0, 48); format long, Sj=double(S1);disp(积分准确值为:),Sjdisp(积分近似值分别为:),T,S,FT8,FS4,quaddisp(与精确值比较的误差为:),wuchi=double(abs(Sj-T),double(abs(Sj-S),double(abs(Sj-FT8),double(abs(Sj-FS4),double(abs(Sj-quad)五、实验数据及结果：运行t1得t1积分准确值为:Sj = 0.68175966492712积分近似值分别为:ans =1.11803398874989 0.68448277848113 0.68874468339797 0.68177045308821 0.68175956505170与精确值比较的误差为:wuchi =0.43627432382277 0.00272311355400 0.00698501847084 0.00001078816109 0.00000009987542运行t2得t2积分准确值为:Sj =-6.28318530717959积分近似值分别为:ans =-0.00000000000000 0 -5.95683320009178 -6.29751019994081 -6.28318521082678与精确值比较的误差为:wuchi = 6.28318530717958 6.28318530717959 0.32635210708781 0.01432489276123 0.00000009635280运行t3得 t3积分准确值为:Sj =10.20759220056126积分近似值分别为:ans =14.23024947075771 10.15174340344855 10.26636718387094 10.20722396340437 10.20759362810855与精确值比较的误差为:wuchi = 4.02265727019644 0.05584879711271 0.05877498330968 0.00036823715689 0.00000142754728运行t4得t4积分准确值为:Sj =0.33333333333333积分近似值分别为:ans =0.37500000000000 0.33333333333333 0.33398437500000 0.33333333333333 0.33333333333333与精确值比较的误差为:wuchi = 0.04166666666667 0 0.00065104166667 0 0运行t5得t5积分准确值为:Sj = 58.47046915489933积分近似值分别为:ans = 62.43737140065479 55.57229168712463 56.26305465206046 56.20282264402146 56.20406238610796与精确值比较的误差为:wuchi =3.96690224575546 2.89817746777469 2.20741450283887 2.26764651087787 2.26640676879137 从积分近似值和比较误差可以看出Simpson公式要比梯形公式精确，一般情况下复合公式要比单公式的精确度高，有时Simpson公司反而比复合梯形公式好，Romberg公式是所有