数值分析课程设计--Romberg积分法.docVIP

数值分析上机实验报告 (十) 学院：数学与信息科学 专业：信息与计算科学 课程设计题目 10. 设计题目：Romberg积分法 设计内容：（1）用几个低阶的Newton-Cotes公式计算积分的近似值，自己设置不同精度要求，对结果进行比较分析。 （2）用Romberg积分法计算积分的近似值，自己设置不同精度要求，对结果进行比较分析；与（1）的结果进行比较分析，谈谈你的体会。 ①Cotes公式积分计算MatLab代码如下： function s=cotes(f,a,b,n) %f:òa?óoˉêy %a,b:???? %n:μè·?êy %μ÷ó???ê?￡o format long %éè?????è syms t h=(b-a)/n; %????2?3¤ C=zeros(1,n+1); %?óNewton-Cotes?μêy j=1; for k=0:n d=1; J=((-1)^(n-k))/(n*prod(1:k)*prod(1:n-k)); for i=0:n if i~=k d=d*(t-i); end end C(j)=J*int(d,0,n); j=j+1; end %?ó?¨?y·? s=0; j=1; for k=0:n s=s+C(j)*feval(f,a+k*h); j=j+1; end s=s*(b-a); %ê?3??á1? %f=inline(exp(-x.^2)); %a=0;b=1;n=4; %cotes(f,a,b,n) 运行结果如下图窗口截图： 流程图如下： ②辛普森（Simpson）积分公式求解MatLab代码如下： function s=simpson(f,a,b,n) h=(b-a)/(2*n); s1=0; s2=0; for k=1:n x=a+h*(2*k-1); s1=s1+feval(f,x); end for k=1:(n-1) x=a+h*2*k; s2=s2+feval(f,x); end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3; syms x f=exp(-x.^2); %f=inline(exp(-x.^2)) %a=0;b=1; %simpson(f,0,1,4) 运行结果如下图窗口截图： 流程图如下： ③梯形公式求解MatLab代码如下： function s=trap(f,a,b,n) h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; syms x f=exp(-x.^2); %f=inline(exp(-x.^2)) %a=0;b=1; %trap(f,a,b,4) 运行结果如下图窗口截图： 龙贝格（Romberg）计算方法MatLab代码如下： %Romberg?y·? ??exp(-x.^2)?ˉ fprintf(??????ó?Romberg??·¨?ó?3??oˉêy′óaμ?b?y·??μμ??ü???a?￡\n) %ê?è?2?êy syms x f=input(ê?è?±??yoˉêyμ?±í′?ê?:) a=input(ê?è??y·????T:) b=input(ê?è??y·?é??T:) epsilon=input(ê?è????èeps:) %±?á?3?ê??ˉ?￡ h=b-a k=1 T(1,1)=(h/2)*(subs(f,a)+subs(f,b)) %μü′ú1y3ì. while(1) T(k+1,1)=0.5*(T(k,1)+h*sum(subs(f,a+(2*[1:(2^(k-1))]-1)*h/2))) for m=1:k T(k-m+1,m+1)=(4^m*T(k-m+2,m)-T(k-m+1,m))/(4^m-1) end if abs(T(1,k+1)-T(1,k))epsilon break; else h=h/2 k=k+1 end end fprintf(??è·?μ?a%9.8f?￡\n,T(1,k+1)) 运行结果如下图窗口截图： 流程图如下： 体会和感想： 在实验中体会到编写程序需要仔细，熟悉matlab的使用以及matlab程序设计，熟悉并掌握复化辛普森公式的原理以及所设计的程序，以及龙贝格方法的原理以及所设计的程序。 调用函数 输入a、b、n 步长h=（b-a）/n X=a+0.5*h S=4*f（x） K=1 K=k+1 X=x+0.5*h s=s+4*f（x） X=x+0.5*h s=s+2*f（x） Kn 得到结果