实验三数值积分与数值微分.docVIP

《数值分析》课程设计实验报告 实验三 数值积分与数值微分 一、问题提出 选用复合梯形公式，复合Simpson公式，Romberg算法，计算 （1） （2） （3） 二、问题分析 由上可知这个积分找不到用初等函数表示的原函数，直接计算起来很困难，因此我们考虑利用函数在若干点得函数值，近似地计算该函数在一个区间上得定积分。这里采用的方法有三种：复合梯形公式，复合Simpson公式，Romberg算法。 1、复合梯形公式MATLAB程序： function I=T_quad(x,y) n=length(x);m=length(y); if n~=m error(the length of X and Y must be equal!); return; end h=(x(n)-x(1))/(n-1); a=[1 2*ones(1,n-2) 1]; I=h/2 * sum(a.*y); 2、复合Simpson公式MATLAB程序： function I=S_quad(x,y) n=length(x);m=length(y); if n~=m error(the length of X and Y must be equal!); return; end if rem(n-1,2)~=0 %如果n-1不能被2整除，则调用复化梯形公式 I=T_quad(x,y); return; end N=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n); for k=1:N a(2*k-1)=a(2*k-1)+1; a(2*k)=a(2*k)+4; a(2*k+1)=a(2*k+1)+1; end I=h/6*sum(a.*y); 3、Romberg算法MATLAB程序： function I=R_quad_iter(fun,a,b,ep) if nargin 4 ep=1e-5; end m=1;h=b-a; I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I; while 1 N=2^(m-1);h=h/2;I=I/2; for i=1:N; I=I+h*feval(fun,a+(2*i-1)*h); end T(m+1,1)=I;M=2*N;k=1; while M1; T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k)/(4^k-1)); M=M/2;k=k+1; end if abs(T(k,k)-T(k-1,k-1))ep break; end m=m+1; end I=T(k,k); 自适应步长梯形公式： function I=R_quad_iter(fun,a,b,ep) if nargin 4 ep=1e-5;end; N=1;h=b-a; T=h/2*(feval(fun,a)+feval(fun,b)); while 1 h=h/2;I=T/2; for k=1:N; I=I+h*feval(fun,a+(2*k-1)*h); end if abs(I-T)ep break; end N=2*N;T=I; end 对积分求解： %对求解 x=0:1/40:1/4; y=sqrt(4-(sin(x)).^2); format long I=T_quad(x,y) %调用复化梯形公式求得积分值I =0.49870482652029 I=S_quad(x,y) %调用复化辛普森公式求得积分值I =0.49871111844568 fun=inline(sqrt(4-(sin(x)).^2)); I=R_quad_iter(fun,0,1) %调用龙贝格算法可得积分值I=0.498711118445678 x=0:1/400:1/4; y=sqrt(4-(sin(x)).^2); I=T_quad(x,y) %调用复化梯形公式求得积分值I = 0.49871105466684 I=S_quad(x,y) %调用复化辛普森公式求得积分值I = 0.49871111757532 %对求解： x=0.1:0.1:1; y=sin(x)./x %可以得到y值，由题可知f(0)=1。 y=[1 0.99833416646828 0.99334665397531 0.98506735553780 0.97354585577163 0.95885107720841 0.