ch4_数值计算2010a.docVIP

数值计算 近似数值极限及导数 【例4.1-1】 x=eps; L1=(1-cos(2*x))/(x*sin(x)), % L2=sin(x)/x, % L1 = 0 L2 = 1 syms t f1=(1-cos(2*t))/(t*sin(t)); f2=sin(t)/t; Ls1=limit(f1,t,0) Ls2=limit(f2,t,0) Ls1 = 2 Ls2 = 1 【例4.1-2】 （1） d=pi/100; t=0:d:2*pi; x=sin(t); dt=5*eps; % x_eps=sin(t+dt); dxdt_eps=(x_eps-x)/dt; % plot(t,x,LineWidth,5) hold on plot(t,dxdt_eps) hold off legend(x(t),dx/dt) xlabel(t) 图 4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线 （2） x_d=sin(t+d); dxdt_d=(x_d-x)/d; % plot(t,x,LineWidth,5) hold on plot(t,dxdt_d) hold off legend(x(t),dx/dt) xlabel(t) 图 4.1-2 增量适当所得导函数比较光滑 【例4.1-3】 clf d=pi/100; % t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; % dxdt_grad=gradient(x)/d; % subplot(1,2,1) plot(t,x,b) hold on plot(t,dxdt_grad,m,LineWidth,8) plot(t(1:end-1),dxdt_diff,.k,MarkerSize,8) axis([0,2*pi,-1.1,1.1]) title([0, 2\pi]) legend(x(t),dxdt_{grad},dxdt_{diff},Location,North) xlabel(t),box off hold off subplot(1,2,2) kk=(length(t)-10):length(t);%t hold on plot(t(kk),dxdt_grad(kk),om,MarkerSize,8) plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSize,8) title([end-10, end]) legend(dxdt_{grad},dxdt_{diff},Location,SouthEast) xlabel(t),box off hold off 图 4.1-3 diff和gradient求数值近似导数的异同比较 数值求和与近似数值积分 【例 4.1-4】 clear d=pi/8; % t=0:d:pi/2; % y=0.2+sin(t); % s=sum(y); % s_sa=d*s; % 6 s_ta=d*trapz(y); % 7 disp([sum求得积分,blanks(3),trapz求得积分]) disp([s_sa, s_ta]) t2=[t,t(end)+d]; % y2=[y,nan]; % stairs(t2,y2,:k) % hold on plot(t,y,r,LineWidth,3) % h=stem(t,y,LineWidth,2); % set(h(1),MarkerSize,10) axis([0,pi/2+d,0,1.5]) % hold off shg sum求得积分 trapz求得积分 1.5762 1.3013 图 4.1-4 sum 和trapz求积模式示意 计算精度可控的数值积分 【例 4.1-5】 （1） syms x Isym=vpa(int(exp(-x^2),x,0,1)) Isym = 0.74682413281242702539946743613185 （2） format long % d=0.001;x=0:d:1; Itrapz=d*trapz(exp(-x.*x)) Itrapz = 0.746824071499185 （3） fx=exp(-x.^2); % Ic=quad(fx,0,1,1e-8) % 7 Ic = 0.746824132854452 【例 4.1-6】 （1） syms x y s=vpa(int(int(