1、1第四章 数值计算4.1 引言本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程（初值和边值问题）求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源 MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和 MATLAB 计算指令之间的内在联系及区别。对于那些熟悉其他高级语言（如 FORTRAN，Pascal，C+）的读者来说，通过本章，MATLAB 卓越的数组处理能力、浩瀚而灵活的 M 函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开
2、朗，感受到摆脱烦琐编程后的眉眼舒展。对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB 精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。对于那些熟悉 MATLAB 基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。由于 MATLAB 的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier 分
3、析，和一般常微分方程初值、边值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。MATLAB 从 5.3 版升级到 6.x 版后，本章内容的变化如下： MATLAB 从 6.0 版起，其矩阵和特征值计算指令不再以 LINPACK 和 EISPACK 库为基础，而建筑在计算速度更快、运行更可靠的 LAPACK 和 ARPACK 程序库的新基础上。因此，虽然各种矩阵计算指令没有变化，
4、但计算结果却可能有某些不同。这尤其突出地表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。对此，用户不必诧异，因为构成空间的基向量时不唯一的，且新版的更可信。本书新版全部算例结果是在 6.x 版上给出的。 在 5.3 版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。这种调用方式在 6.x 版中已被宣布为“过渡期内允许使用但即将被淘汰的调用方式”；而新的调用方式是借助“函数句柄”进行的。因此，关于述泛函指令，本章新版着重讲述如何使用“函数句柄”，同时兼顾“函数名字符串”调用法。 MATLAB 从 6.0 版起，提供了一组专门求微分方程 “边值问题”数值解的指令。适应这种变化，本章
5、新增第 4.14.5 节，用 2 个算例阐述求解细节。 5.3 版中的积分指令 quad8 已经废止；6.x 版启用新积分指令 quadl ；6.5 版新增三重积分指令 triplequad。本章新版对此作了相应的改变。 4.2 LU 分解和恰定方程组的解4.2.1 LU 分解、行列式和逆24.2.2 恰定方程组的解【例 4.2.2-1】“求逆”法和“ 左除”法解恰定方程的性能对比（1）randn(state,0);A=gallery(randsvd,100,2e13,2);x=ones(100,1);b=A*x;cond(A) ans =1.9990e+013 （2）ticxi=inv(A)
6、*b;ti=toceri=norm(x-xi)rei=norm(A*xi-b)/norm(b) ti =0.7700eri =0.0469rei =0.0047 （3）tic;xd=Ab;td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b) td =0erd =0.0078red =2.6829e-015 4.2.3 范数、条件数和方程解的精度【例 4.2.3-1】Hilbert 矩阵是著名的病态矩阵。MATLAB 中有专门的 Hilbert 矩阵及其准确逆矩阵的生成函数。本例将对方程 近似解和准确解进行比较。bHxN=6 8 10 12 14;for
7、k=1:length(N)n=N(k);H=hilb(n);Hi=invhilb(n);b=ones(n,1);x_approx=Hb; x_exact=Hi*b;ndb=norm(H*x_approx-b);nb=norm(b);ndx=norm(x_approx - x_exact);nx=norm(x_approx);er_actual(k)=ndx/nx;K=cond(H);er_approx(k)=K*eps;er_max(k)=K*ndb/nb;end disp(Hilbert 矩阵阶数),disp(N)format short e3disp(实际误差 er_actual),dis
8、p(er_actual),disp()disp(近似的最大可能误差 er_approx),disp(er_approx),disp()disp(最大可能误差 er_max),disp(er_max),disp() Hilbert 矩阵阶数6 8 10 12 14实际误差 er_actual1.5410e-010 1.7310e-007 1.9489e-004 9.1251e-002 2.1257e+000近似的最大可能误差 er_approx3.3198e-009 3.3879e-006 3.5583e-003 3.9846e+000 9.0475e+001最大可能误差 er_max7.949
9、8e-007 3.8709e-002 1.2703e+003 4.7791e+007 4.0622e+010 4.3 矩阵特征值和矩阵函数4.3.1 特征值和特征向量的求取【例 4.3.1-1】简单实阵的特征值问题。A=1,-3;2,2/3;V,D=eig(A) V =0.7746 0.7746 0.0430 - 0.6310i 0.0430 + 0.6310iD =0.8333 + 2.4438i 0 0 0.8333 - 2.4438i 【例 4.3.1-2】本例演示：如矩阵中有元素与截断误差相当时的特性值问题。A=3 -2 -0.9 2*eps-2 4 -1 -eps-eps/4 eps
10、/2 -1 0-0.5 -0.5 0.1 1 ;V1,D1=eig(A);ER1=A*V1-V1*D1V2,D2=eig(A,nobalance);ER2=A*V2-V2*D2 ER1 =0.0000 0.0000 0.0000 0.00000 -0.0000 -0.0000 -0.00000.0000 -0.0000 -0.0000 0.00000.0000 0.0000 0.0000 -0.5216ER2 =1.0e-014 *-0.2665 0.0111 -0.0559 -0.10550.4441 0.1221 0.0343 0.08330.0022 0.0002 0.0007 00.0
11、194 -0.0222 0.0222 0.0333 【例 4.3.1-3】指令 eig 与 eigs 的比较。rand(state,1),A=rand(100,100)-0.5;t0=clock;V,D=eig(A);T_full=etime(clock,t0)options.tol=1e-8;options.disp=0;t0=clock;v,d=eigs(A,1,lr,options);T_part=etime(clock,t0)Dmr,k=max(real(diag(D);d,D(1,1) T_full =40.2200T_part =3.1300d =3.0140 + 0.2555ia
12、ns =3.0140 + 0.2555i vk1=V(:,k);vk1=vk1/norm(vk1);v=v/norm(v);V_err=acos(norm(vk1*v)*180/piD_err=abs(D(k,k)-d)/abs(d) V_err =1.2074e-006D_err =4.2324e-010 4.3.2 特征值问题的条件数【例 4.3.2-1】矩阵的代数方程条件数和特征值条件数。B=eye(4,4);B(3,4)=1;Bformat short e,c_equ=cond(B),c_eig=condeig(B) B =1 0 0 00 1 0 00 0 1 10 0 0 1c_e
13、qu =2.6180e+000Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 1.110223e-016. In D:MATLAB6P1toolboxmatlabmatfuncondeig.m at line 30c_eig =1.0000e+0001.0000e+0004.5036e+0154.5036e+015 【例 4.3.2-2】对亏损矩阵进行 Jordan 分解。A=gallery(5)VJ,DJ=jordan(A);V,D,c_eig=condeig(A)
14、;c_equ=cond(A);DJ,D,c_eig,c_equ A =-9 11 -21 63 -25270 -69 141 -421 1684-575 575 -1149 3451 -138013891 -3891 7782 -23345 933651024 -1024 2048 -6144 24572DJ =0 1 0 0 00 0 1 0 00 0 0 1 00 0 0 0 10 0 0 0 0D =Columns 1 through 4 5-0.0408 0 0 0 0 -0.0119 + 0.0386i 0 0 0 0 -0.0119 - 0.0386i 0 0 0 0 0.0323
15、 + 0.0230i0 0 0 0 Column 5 0 0 0 0 0.0323 - 0.0230ic_eig =1.0e+010 *2.12932.07962.07962.00202.0020c_equ =2.0253e+018 4.3.3 复数特征值对角阵与实数块特征值对角阵的转化【例 4.3.3-1】把例 4.3.1-1 中的复数特征值对角阵 D 转换成实数块对角阵，使VR*DR/VR=A。A=1,-3;2,2/3;V,D=eig(A);VR,DR=cdf2rdf(V,D) VR =0.7746 00.0430 -0.6310DR =0.8333 2.4438-2.4438 0.833
16、3 4.3.4 矩阵的谱分解和矩阵函数【例 4.3.4-1】数组乘方与矩阵乘方的比较。clear,A=1 2 3;4 5 6;7 8 9;A_Ap=A.0.3A_Mp=A0.3 A_Ap =1.0000 1.2311 1.39041.5157 1.6207 1.71181.7928 1.8661 1.9332A_Mp =0.6962 + 0.6032i 0.4358 + 0.1636i 0.1755 - 0.2759i0.6325 + 0.0666i 0.7309 + 0.0181i 0.8292 - 0.0305i0.5688 - 0.4700i 1.0259 - 0.1275i 1.483
17、0 + 0.2150i 【例 4.3.4- 2】标量的数组乘方和矩阵乘方的比较。（ A 取自例 4.3.4-1）pA_A=(0.3).ApA_M=(0.3)A pA_A =0.3000 0.0900 0.02700.0081 0.0024 0.000760.0002 0.0001 0.0000pA_M =2.9342 0.4175 -1.0993-0.0278 0.7495 -0.4731-1.9898 -0.9184 1.1531 【例 4.3.4-3】sin 的数组运算和矩阵运算比较。（A 取自例 4.3.4-1）A_sinA=sin(A)A_sinM=funm(A,sin) A_sinA
18、 =0.8415 0.9093 0.1411-0.7568 -0.9589 -0.27940.6570 0.9894 0.4121A_sinM =-0.6928 -0.2306 0.2316-0.1724 -0.1434 -0.11430.3479 -0.0561 -0.4602 4.4 奇异值分解4.4.1 奇异值分解和矩阵结构4.4.1.1 奇异值分解定义4.4.1.2 矩阵结构的奇异值分解描述4.4.2 线性二乘问题的解4.4.2.1 矩阵除运算的广义化4.4.2.2 线性模型的最小二乘解【例 4.4.2.2-1】对于超定方程 ，进行三种解法比较。其中 取 MATLAB 库中的特AxyA
19、殊函数生成。（1）A=gallery(5);A(:,1)=;y=1.7 7.5 6.3 0.83 -0.082;x=inv(A*A)*A*y,xx=pinv(A)*y,xxx=Ay Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 5.405078e-018.x =3.4811e+0005.1595e+0009.5340e-001-4.6569e-002xx =3.4759e+0005.1948e+0007.1207e-001-1.1007e-001Warning: R
20、ank deficient, rank = 3 tol = 1.0829e-010.xxx =3.4605e+00075.2987e+0000-2.9742e-001 （2）nx=norm(x),nxx=norm(xx),nxxx=norm(xxx) nx =6.2968e+000nxx =6.2918e+000nxxx =6.3356e+000 （3）e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx) e =6.9863e-001ee =4.7424e-002eee =4.7424e-002 4.5 函数的数值导数和切平面4.5.1 法线【例 4.
21、5.1-1】曲面法线演示。y=-1:0.1:1;x=2*cos(asin(y);X,Y,Z=cylinder(x,20);surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21);view(120,18) 05 0 1 02 4 6 8 12 图 4.5.1-1 4.5.2 偏导数和梯度4.5.2.1 理论定义4.5.2.2 数值计算指令【例 4.5.2.2-1】用一个简单矩阵表现 diff 和 gradient 指令计算方式。F=1,2,3;4,5,6;7,8,98Dx=diff(F)Dx_2=diff(F,1,2)FX,FY=gradient(F)FX_2,FY_
22、2=gradient(F,0.5) F =1 2 34 5 67 8 9Dx =3 3 33 3 3Dx_2 =1 11 11 1FX =1 1 11 1 11 1 1FY =3 3 33 3 33 3 3FX_2 =2 2 22 2 22 2 2FY_2 =6 6 66 6 66 6 6 【例 4.5.2.2-2】研究偶极子 (Dipole)的电势（Electric potential）和电场强度（Electric field density）。设在 处有电荷 ，在 处有电荷 。那么在电荷所在平面上任),(baq),(baq何一点的电势和场强分别为 ， 。其中14),(0ryxVVE。 。又
23、设电荷2222 )(,)()( yryxr 9014， ， 。610q5.a.1bclear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x;X,Y=meshgrid(x,y);rp=sqrt(X-a).2+(Y-b).2);rm=sqrt(X+a).2+(Y+b).2);V=q*k*(1./rp-1./rm);Ex,Ey=gradient(-V);AE=sqrt(Ex.2+Ey.2);Ex=Ex./AE;Ey=Ey./AE;cv=linspace(min(min(V),max(max(V),49);contourf(X,Y,V,cv,k-)%axis
24、(square)title(fontname隶书fontsize22偶极子的场),hold onquiver(X,Y,Ex,Ey,0.7)plot(a,b,wo,a,b,w+)plot(-a,-b,wo,-a,-b,w-)xlabel(x);ylabel(y),hold off 9 （2）a=0.1;b=0.5;t=-10:0.01:10;y_char=vectorize(y); % Y=feval(y_char,t,a,b);clf,plot(t,Y,r);hold on,plot(t,zeros(size(t),k);xlabel(t);ylabel(y(t),hold off 5 图 4
25、.6-1 （3）由于 Notebook 中无法实现 zoom、ginput 指令涉及的图形和鼠标交互操作，因此下面指令必须在 MATLAB 指令窗中运行，并得到如图 4.6-2 所示的局部放大图及鼠标操作线。zoom on10tt,yy=ginput(5);zoom off图 4.6-2 tt tt =-2.0032-0.5415-0.00720.58761.6561 （4）t4,y4,exitflag=fzero(y,tt(4),a,b) % t4 =0.5993y4 =0exitflag =1 （5）t3,y3,exitflag=fzero(y,tt(3),a,b) t3 =-0.5198
26、y3 =-5.5511e-017exitflag =1 （6）op=optimset(fzero) op = ActiveConstrTol: Display: notifyTolX: 2.2204e-016TypicalX: op=optimset(tolx,0.01);op.TolX ans =0.0100 （7）t4n,y4n,exitflag=fzero(y,tt(4),op,a,b) t4n =110.6042y4n =0.0017exitflag =1 4.6.3 多元函数的零点【例 4.6.3-1】求解二元函数方程组 的零点。0)cos(),(in21yxyf（1）x=-2:0.
27、5:2;y=x;X,Y=meshgrid(x,y);F1=sin(X-Y);F2=cos(X+Y);v=-0.2, 0, 0.2;contour(X,Y,F1,v)hold on,contour(X,Y,F2,v),hold off 0 5 1 5 2 05 15 2图 4.6-3 （2）x0,y0=ginput(2); disp(x0,y0) -0.7926 -0.78430.7926 0.7843 （3）fun=sin(x(1)-x(2),cos(x(1)+x(2); %xy,f,exit=fsolve(fun,x0(2),y0(2) % Optimization terminated s
28、uccessfully:First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detectedxy =0.7854 0.7854f =1.0e-006 *-0.0984 0.2011exit =1 说明fun.mfunction ff=fun(x)ff(1)=sin(x(1)-x(2);12ff(2)=cos(x(1)+x(2);4.7 函数极值点4.7.1 一元函数的极小值点4.7.2 多元函数的极小值点【例 4.7.2-1】求 的极小值点。它即是著名的22)1()(10),( xyxf
29、 Rosenbrocks “Banana“ 测试函数。该测试函数有一片浅谷，许多算法难以越过此谷。（演示本例搜索过程的文件名为 exm04072_1_1.m 。）（1）ff=inline(100*(x(2)-x(1)2)2+(1-x(1)2,x); （2）x0=-1.2,1;sx,sfval,sexit,soutput=fminsearch(ff,x0) sx =1.0000 1.0000sfval =8.1777e-010sexit =1soutput = iterations: 85funcCount: 159algorithm: Nelder-Mead simplex direct se
30、arch （3）ux,sfval,uexit,uoutput,grid,hess=fminunc(ff,x0) Warning: Gradient must be provided for trust-region method; using line-search method instead. In D:MATLAB6P1toolboxoptimfminunc.m at line 211Optimization terminated successfully:Current search direction is a descent direction, and magnitude of
31、directional derivative in search direction less than 2*options.TolFunux =1.0000 1.0000sfval =1.9116e-011uexit =1uoutput = iterations: 26funcCount: 162stepsize: 1.2992firstorderopt: 5.0020e-004algorithm: medium-scale: Quasi-Newton line searchgrid =1.0e-003 *-0.5002-0.1888hess =820.4028 -409.5496-409.
32、5496 204.7720 134.8 数值积分4.8.1 一元函数的数值积分4.8.1.1 闭型数值积分【例 4.8.1.1-1】求 ，其精确值为 。dxeI102 746820.（1）syms x;IS=int(exp(-x*x),x,0,1)vpa(IS) IS =1/2*erf(1)*pi(1/2)ans =.74682413281242702539946743613185 （2）fun=inline(exp(-x.*x),x);Isim=quad(fun,0,1),IL=quadl(fun,0,1) Isim =0.7468IL =0.7468 （3）Ig=gauss10(fun,0
33、,1) Ig =0.7463 （4）xx=0:0.1:1.5;ff=exp(-xx.2);pp=spline(xx,ff);int_pp=fnint(pp);Ssp=ppval(int_pp,0,1)*-1;1 Ssp =0.7468 （5）图 4.8-1 4.8.1.2 开型数值积分gauss10.mfunction g = gauss10(fun,a,b)%GAUSS10(fun,a,b)% fun%=14x = 0.1488743390;0.4333953941;0.6974095683;.0.8650633667;0.9739065285;w = 0.2955242247;0.2692
34、667193;0.2190863625;.0.1494513492;0.0666713443;t = .5*(b+a)+.5*(b-a)*-flipud(x);x;W = flipud(w);w;g = sum(W.*feval(fun,t)*(b-a)/2;【例 4.8.1.2-1】当 时，比较解析积分和近似积分。)cos()(xf（1）syms x;F=int(cos(x),x,-1,1),vpa(F) F =2*sin(1)ans =1.6829419696157930133050046432606 （2）aF=cos(1/sqrt(3)+cos(-1/sqrt(3) aF =1.675
35、8 【例 4.8.1.2-2】求 ，准确结果是 。10lndxI .8629（1）syms x;IS=vpa(int(sqrt(log(1/x),x,0,1) Warning: Explicit integral could not be found. In D:MATLAB6P5toolboxsymbolicsymint.m at line 58In D:MATLAB6P5toolboxsymboliccharint.m at line 9IS =.88622692545275801364908374167057 （2）用 quad 指令求积ff=inline(sqrt(log(1./x),
36、x);Isim=quad(ff,0,1) Warning: Divide by zero. In D:MATLAB6P5toolboxmatlabfunfuninlineeval.m at line 13In D:MATLAB6P5toolboxmatlabfunfuninlinefeval.m at line 34In D:MATLAB6P5toolboxmatlabfunfunquad.m at line 62Isim =0.8862 （3）Ig=gauss10(ff,0,1) Ig =0.8861 4.8.2 多重数值积分4.8.2.1 积分限为常数的二重积分指令【例 4.8.2.1-1
37、】计算 和 。2100dyxSyx 1022dyxSyx（1）syms x yssx01=vpa(int(int(xy,x,0,1),y,1,2)15ssx12=vpa(int(int(xy,y,0,1),x,1,2) Warning: Explicit integral could not be found. In D:MATLAB6P5toolboxsymbolicsymint.m at line 58ssx01 =.40546510810816438197801311546435ssx12 =1.2292741343616127837489278679215 （2）zz=inline(x
38、.y,x,y);nsx01=dblquad(zz,0,1,1,2)nsx12=dblquad(zz,1,2,0,1) nsx01 =0.4055nsx12 =1.2293 4.8.2.2 内积分限为函数的二重积分double_int.mfunction SS=double_int(fun,innlow,innhi,outlow,outhi)%double_int%fun%innlow,innhi%outlow,outhiy1=outlow;y2=outhi;x1=innlow;x2=innhi;f_p=fun;SS=quad(G_yi,y1,y2,x1,x2,f_p); G_yi.mfunct
39、ion f=G_yi(y,x1,x2,f_p)%G_yi%y%x1,x2%f_py=y(:);n=length(y);if isnumeric(x1)=1;xx1=x1*ones(size(y);else xx1=feval(x1,y);endif isnumeric(x2)=1;xx2=x2*ones(size(y);else xx2=feval(x2,y);endfor i=1:nf(i)=quad(f_p,xx1(i),xx2(i),y(i);endf=f(:);【例 4.8.2.2-1】计算 。dyxIy4122)(（1）x_low.mfunction f=x_low(y)f=sqrt
40、(y);（2）16（3）ff=inline(x.2+y.2,x,y);SS=double_int(ff,x_low,2,1,4) Warning: Minimum step size reached; singularity possible. In D:MATLAB6P5toolboxmatlabfunfunquad.m at line 88In D:MATLAB6p5workG_yi.m at line 11In D:MATLAB6P5toolboxmatlabfunfuninlinefeval.m at line 20In D:MATLAB6P5toolboxmatlabfunfunqu
41、ad.m at line 62In D:MATLAB6p5workdouble_int.m at line 8SS =9.5810 （4）Ssym=vpa(int(int(x2+y2,x,sqrt(y),2),y,1,4) Ssym =9.5809523809523809523809523809524 4.8.3 卷积4.8.3.1 “完整”离散序列的数值卷积4.8.3.2 “截尾”离散序列的数值卷积4.8.3.3 多项式乘法与离散卷积的算法同构4.8.3.4 连续时间函数的数值卷积4.8.3.5 卷积的 MATLAB 实现【例 4.8.3.5-1】有序列 和 。elsnA12,4301)(
42、elsnB9,3201)(（A）%a=ones(1,10);n1=3;n2=12;b=ones(1,8);n3=2;n4=9;c=conv(a,b);nc1=n1+n3;nc2=n2+n4;kc=nc1:nc2;%aa=a(1:6);nn1=3;nn2=8;cc=conv(aa,b);ncc1=nn1+n3;nx=nn2+n4;ncc2=min(nn1+n4,nn2+n3);kx=(ncc2+1):nx;kcc=ncc1:ncc2;N=length(kcc);stem(kcc,cc(1:N),r,filled)axis(nc1-2,nc2+2,0,10),grid,hold onstem(kc
43、,c,b),stem(kx,cc(N+1:end),g),hold off 17 6 图 4.8-2 【例 4.8.3.5-2】求函数 和 的卷积。本例展)()(tUetu)(2/tUteh示：（A）符号 Laplace 变换求卷积的理论表示；（B）SIMULINK 卷积法的执行过程和它的快速精确性。（C）从理论符号解产生相应的理论数值序列。（1）syms tao;t=sym(t,positive);US1=laplace(exp(-t);HS1=laplace(t*exp(-t/2)yt1=simple(ilaplace(US1*HS1) HS1 =1/(1/2+s)2yt1 =4*exp(
44、-t)+(2*t-4)*exp(-1/2*t) （2）图 4.8-3 （3）t=yt2(:,1);yyt1=eval(vectorize(char(yt1);dy,kd=max(abs(yyt1-yt2(:,2);dy12=dy/abs(yyt1(kd) dy12 =2.8978e-006 【例 4.8.3.5-3】用“零阶”近似法求 和 的卷积。本例演)()(tUetu)(2/tteh示：（A）连续函数的有限长度采样。（ B）卷积数值计算三个误差（“截尾”误差、“零阶”近似误差、计算机字长误差）的影响。（C）卷积“无截尾误差”区间、“非平凡”区间端点的确定。（D）离散卷积和连续卷积之间的关系
45、。（E）指令 conv 的使用。（F）绘图分格线的运用。（1）（2）%t2=3;t4=11;T=0.01;tu=0:T:t2;N2=length(tu);th=0:T:t4;N4=length(th);18u=exp(-tu);h=th.*exp(-th/2);tx=0:T:(t2+t4);Nx=length(tx);yt3=T*conv(u,h);%t=tx;yyt1=eval(vectorize(char(yt1);dy,kd=max(abs(yyt1(1:N2)-yt3(1:N2);dy13(1)=dy/abs(yyt1(kd);dy,kd=max(abs(yyt1(N2+1:N4)-y
46、t3(N2+1:N4);dy13(2)=dy/abs(yyt1(N2+kd);dy,kd=max(abs(yyt1(N4+1:Nx)-yt3(N4+1:Nx);dy13(3)=dy/abs(yyt1(N4+kd); （3）disp(与理论结果的相对误差)disp(blanks(4),0,3段 3,11段 11,14段),disp(dy13)plot(t,yyt1,:b,tx,yt3,r)set(gca,Xtick,0,3,11,14),grid 与理论结果的相对误差0,3段 3,11 段 11,14 段0.0068 0.0810 0.69740 3 01 2 3 4 5 图 4.8.3.5-3-1 4.9 随机数据的统计描述4.9.1 统计分布的数字特征【
