数值分析实验报告4.docVIP

实验四 数值积分 实验目的 了解并掌握matlab软件的基本编程、操作方法； 初步了解matlab中的部分函数，熟悉循环语句的使用； 通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积方法计算积分的原理。 实验题目： 1.用不同数值方法计算积分 ADVANCE ln xdx=-. 取不同的步长h.分别用复合梯形及辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？ 用龙贝格求积计算完成问题（1）。 1.实验原理 复合梯形公式：Tn=2； 复合辛普森公式：Sn=[f(a)+f(b)+2+4]； 以上两种算法都是将a-b之间分成多个小区间（n）,则h=(b-a)/n,xk=a+kh, xk+1/2=a+(k+1/2)h,利用梯形求积根据两公式便可。 龙贝格算法：在指定区间内将步长依次二分的过程中运用如下公式 Sn=T2n-Tn Cn=S2n-Sn Rn=C2n-Cn 2.程序设计 （1）、复合梯形法： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h); t=h*(0.5*(fa+fb)+sum(f)); （2）、复合辛普森法： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0.001*h); f2=feval(fname,a+h/2:h:b-h+0.001*h); t=h/6*(fa+fb+2*sum(f1)+4*sum(f2)); （3）龙贝格法： function [I,step]=Roberg(f,a,b,eps) 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=2*M; 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)); end I=T(k+1,k+1); step=k; 3.实验结果； （1）复合梯形法 format long;natrapz(inline(sqrt(x).*log(x)),eps,1,10),format short; ans = -0.41706283177947 format long;natrapz(inline(sqrt(x).*log(x)),eps,1,100),format short; ans = -0.44311790800816 format long;natrapz(inline(sqrt(x).*log(x)),eps,1,1000),format short; ans = -0.44438753899716 (2)、复合辛普森法结果： format long;natrapz(inline(sqrt(x).*log(x)),eps,1,10),format short; ans = -0.43529789007469 format long;natrapz(inline(sqrt(x).*log(x)),eps,1,100),format short; ans = -0.44416117841567 format long;natrapz(inline(sqrt(x).*log(x)),eps,1,1000),format short; ans = -0.44443411761418 (3)龙贝格法结果 [q,s]=Roberg(sqrt(x)*log(x),0.0000001,1) q =