第4章 MATLAB的数值计算 41数值微积分 42矩阵和代数方程 43概率分布和统计分析 44多项式运算和卷积 2021/1/26 第1页
2021/1/26 第1页 第4章 MATLAB 的数值计算 4.1 数值微积分 4.2 矩阵和代数方程 4.3 概率分布和统计分析 4.4 多项式运算和卷积
41数值微积分 411近似数值极限及导数 412数值求和与近似数值积分 413计算精度可控的数值积分 414函数极值的数值求解 415常微分方程的数值解 2021/1/26 第2页
2021/1/26 第2页 4.1数值微积分 4.1.1 近似数值极限及导数 4.1.2 数值求和与近似数值积分 4.1.3 计算精度可控的数值积分 4.1.4 函数极值的数值求解 4.1.5 常微分方程的数值解
411近似数值极限及导数 ◆在 MATLAB数值计算中,既没有专门的求极限指令, 也没有专门的求导指令。但 MATLAB提供了与“求导” 概念有关的“求差分”指令。 ◆dx=difx)%计算向量X的前向差分 ◆FX= gradient(F)%求一元(函数梯度 ◆[FX,FY]= gradien(F)%求二元(函数)梯度 ◆对df言,当X是向量时,dx=X(2:n-X(1:n-1);当X 是矩阵时,dx=X(2:n,;)-X(1:n1,;)。dx的长度比x 的长度少1个元素 2021/1/26 第3页
2021/1/26 第3页 ◆在MATLAB数值计算中，既没有专门的求极限指令， 也没有专门的求导指令。但MATLAB提供了与“求导” 概念有关的“求差分”指令。 ◆dx=diff(X) %计算向量X的前向差分 ◆FX=gradient(F) %求一元(函数)梯度 ◆[FX, FY] =gradient(F) %求二元（函数）梯度 ◆对diff而言，当X是向量时，dx= X(2:n)-X(1:n-1) ;当X 是矩阵时，dx= X(2:n, :)-X(1:n-1, :) 。 dx的长度比x 的长度少1个元素。 4.1.1 近似数值极限及导数
411近似数值极限及导数 ◆dx=dif(x)%计算向量X的前向差分 ◆FX= gradient(F)%求一元(函数)梯度 ◆[FX,FY]= gradient(F)%求二元(函数)梯度 ◆对 gradien而言,当F是向量时,FX(1)=F(2)F(1 FX(2: end-1)=(F(3: end)-F(I: end-2))/2, FX(end)=F(end)F(end-1);FX长度与F的长度相同 ◆当F是矩阵时,FX,FY是与F同样大小的矩阵。FX的 每行对应F相应行元素间的梯度;FY的每列对应F相应 列元素间的梯度 2021/1/26 第4页
2021/1/26 第4页 ◆dx=diff(X) %计算向量X的前向差分 ◆FX=gradient(F) %求一元(函数)梯度 ◆[FX, FY] =gradient(F) %求二元（函数）梯度 ◆对gradient而言，当F是向量时，FX(1) = F(2)-F(1), FX(2:end-1) = (F(3:end)-F(1:end-2))/2 , FX(end) = F(end)-F(end-1) ; FX长度与F的长度相同 ◆当F是矩阵时， FX, FY是与F同样大小的矩阵。 FX的 每行对应F相应行元素间的梯度 ; FY的每列对应F相应 列元素间的梯度 ; 4.1.1 近似数值极限及导数
数值极限和导数的应用应十分谨慎 1-cos 2x sii元 【例4-1】设f1(x)= f2(x) ,试用机器零阈值eps替代理论0计算极 xsin x 限L1(0)=1imf1(x),L2(0)=limf2(x) x→0 x→0 xeps syms t LI=(1-coS(2*x))/(x*sin(x)) cOs(2 t/t*sin(t)); L2=sin(x)x, f2=sin(t)/t L1= LsI=limit(fl, t. 52=limit(f2, t,0) L2 LS Ls2 X=pi/1000;,%可得到正确结果 2021/1/26 第5页
2021/1/26 第5页 数值极限和导数的应用应十分谨慎 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 x=pi/1000; %可得到正确结果
数值极限和导数的应用应十分谨慎 【例41-2】已知x=si(t),求该函数在区间[02z] 中的近似导函数。 %(1)自变量的增量取得过小(eps数量级) d=pi/100 x(t) t=0:d:2*pi; X=sIn(t, dt=eps; x eps=sin(t+dt) dxdt eps=(x_eps-x)/dt plot(t, x, Line Width, 5) hold on plot(t, dxdt eps hold off 数值导数受计算中有限精度困扰,当增量d过小时 legend(x()ydxt)frt)与f的数值+分接近,高位有效数字完全相 同,d=f(t+l)-f(0)造成高位有效数字消失,精度 xlabel(t 急剧变差。 2021/1/26
2021/1/26 第6页 数值极限和导数的应用应十分谨慎 %（1）自变量的增量取得过小（eps数量级） 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-2】已知 x = sin( t) [0, 2 ] , 求该函数在区间 中的近似导函数。 0 1 2 3 4 5 6 7 -1.5 -1 -0.5 0 0.5 1 t x(t) dx/dt 数值导数受计算中有限精度困扰，当增量dt过小时， f(t+dt)与f(t)的数值十分接近，高位有效数字完全相 同， df =f(t+dt)-f(t) 造成高位有效数字消失，精度 急剧变差
数值极限和导数的应用应十分谨慎 【例41-2】已知x=si(t),求该函数在区间[02z] 中的近似导函数。 %(2)自变量的增量取得适当d=p100 x d=sin(t+d dxdt d=x_ d-x)/d plot(t, x, Line Width, 5) = wat hold on plot(t, dxdt d) hold off 0.6 legend (x(t),dx/, xlabel(,) 2021/1/26 第7页
2021/1/26 第7页 数值极限和导数的应用应十分谨慎 %（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】已知 x = sin( t) [0, 2 ] , 求该函数在区间 中的近似导函数。 0 1 2 3 4 5 6 7 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 t x(t) dx/dt d=pi/100;
end-10, end 【例41-3】已知x=sin d 函数在区间[0,2z]中的 098 0.97 d=pi/100;t=0:d:2pi; x=sin(t); 096 dxdt diff=diff(x)/d; 0950 ○ddt dxdt grad=gradient(x)/d; subplot(1, 2, 1); plot(t, x, b); hold on plot(t, dxdt grad, 'm, Line Width, 8) 宏观上, difffigradient结果大致相同 plot(t(1: end-1), dxdt diff, k, Marker Size, 8) axis([0, 2*pi, 1.1, 1.1 D; title( 0, 2\pi) legend(x(t),dxdt grad , dxdt diff, Location, North) xlabel('t), box off; hold off subplot(1, 2, 2) 细节上, diffFHgradient数值有差异, kk(length(t)-10):length(t); df没有给出最后一点导数 hold on; plot(t(kk), dxdt grad(kk),om, MarkerSize, 8) plot(t(kk-1), dxdt diff(kk-1), k, Marker Size, 8) titleflend-10, end) legend(dxdt grad ,dxdt diff, Location, South East) xlabel('t), box off; hold off 页
2021/1/26 第8页 d=pi/100; t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; dxdt_grad=gradient(x)/d; 【例4.1-3】已知 采用diff和gradient计算该 函数在区间 中的近似导函数。 x = sin( t) [0, 2 ] 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); 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 宏观上, diff和gradient结果大致相同 细节上, diff和gradient数值有差异， diff没有给出最后一点导数
412数值求和与近似数值积分 S RESum( S= dt*sum(X),S=Sum(t2X)%近似矩形法求积 %sum按列向求和得Xn)数组Sx(k)=∑Xm(,k) Scs=cumsum(X) Scs=cumsum(t, x)=dt* cumsum(X) %沿X列向求累计和,仍是数组,第(〔,k)个元素是X数组第k列前f个 元素的和。最后一行等于Sx 12 21 22 (x1+x2)/2 (x2+x2)/2 Xr X (x1+x2+x21+x)2(x2+x21+x2+x2)/2 St= trapz(t2X)或St=dt*rapz(X)%梯形法求積分 Sctcumtrapz(t, X)E Sctdt*cumtrapz(X) %梯形法沿列向求X关于x的累计积分最后一个值等于St 20∠I
2021/1/26 第9页 Sx=sum(X) % sum按列向求和得(1×n)数组 Scs=cumsum(X) %沿X列向求累计和, 仍是数组, 第(i, k)个元素是X数组第k列前i个 元素的和。最后一行等于Sx 4.1.2 数值求和与近似数值积分 1 ( ) ( , ) m m n i Sx k X i k = = 11 12 21 22 31 32 x x X x x x x = St=trapz(t,X) 或 St=dt*trapz(X) %梯形法求積分 Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X关于x的累计积分,最后一个值等于St 11 12 11 21 12 22 11 21 31 12 22 32 x x Scs x x x x x x x x x x = + + + + + + 11 21 12 22 11 21 21 31 12 22 22 32 0 0 ( ) / 2 ( ) / 2 ( ) / 2 ( ) / 2 Sct x x x x x x x x x x x x = + + + + + + + + S=dt*sum(X), S=sum(t,X) %近似矩形法求积分 Scs=cumsum(t,X) =dt*cumsum(X)
丌/2 【例414】求积分s(x)=y(t)dl,其中y=02+S( clear; d=pi/8; t=0: d: pi/2; y=0. 2+sin(t); sum(y);ssa=d*s;%ssa=sum(t,y),近似矩形法积分 sta= trapa,y);%梯形法积分 disp('Smm求得积分,bank3), trapz求得积分) disp(ls sa, s taD) t2=t, t(end)+d]; y2=y, nan]; stairs(t2, y2,: k); hold on plot(t, y,'r, ' Line width, 3) h=stem(t, y, Line width, 2) set(h(i), MarkerSize, 10) axis(|0,pi/2+d,0,1S|) hold off; shg sum求得积分 trapz求得积分 15762 1.3013 第10页
2021/1/26 第10页 clear; d=pi/8; t=0:d:pi/2; y=0.2+sin(t); s=sum(y); s_sa=d*s; % s_sa=sum(t, y), 近似矩形法积分 s_ta=trapz(t,y); %梯形法积分 = / 2 0 ( ) ( ) 【例4.1 s x y t dt y = 0.2 + sin( t) -4】求积分 , 其中 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 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.5 1 1.5