MATLAB教程2012a习题解答1-7章完整版-张志涌-北航[题型借鉴]

《MATLAB教程2012a习题解答1-7章完整版-张志涌-北航[题型借鉴]》由会员分享，可在线阅读，更多相关《MATLAB教程2012a习题解答1-7章完整版-张志涌-北航[题型借鉴]（88页珍藏版）》请在装配图网上搜索。

1、l MATLAB R2012a l 课后习题答案全解 l 第一章 基础准备及入门习题1及解答n 1.数字1.5e2，1.5e3 中的哪个与1500相同吗？解答1.5e3n 2.请指出如下5个变量名中，哪些是合法的？abcd-2 xyz_3 3chan a变量 ABCDefgh解答2、5是合法的。n 3.在MATLAB环境中，比1大的最小数是多少？解答1+epsn 4.设 a = -8 , 运行以下三条指令，问运行结果相同吗？为什么？w1=a(2/3)w2=(a2)(1/3)w3=(a(1/3)2解答（1）不同。具体如下w1=a(2/3)%仅求出主根w2=(a2)(1/3)%求出(-8)2的主根

2、w3=(a(1/3)2%求出(-8)主根后再平方 w1 = -2.0000 + 3.4641iw2 = 4.0000w3 = -2.0000 + 3.4641i （2）复数的多方根的，下面是求取全部方根的两种方法：（A）根据复数方根定义a=-8;n=2;m=3;ma=abs(a);aa=angle(a);for k=1:m%m决定循环次数sa(k)=(aa+2*pi*(k-1)*n/m;%计算各根的相角endresult=(ma(2/3).*exp(j*sa)%计算各根 result = -2.0000 + 3.4641i 4.0000 - 0.0000i -2.0000 - 3.4641i

3、（B）利用多项式求根p=1,0,0,-a2;r=roots(p) r = -2.0000 + 3.4641i -2.0000 - 3.4641i 4.0000 n 5.指令clear, clf, clc各有什么用处？解答clear 清除工作空间中所有的变量。clf 清除当前图形。clc 清除命令窗口中所有显示。n 6.以下两种说法对吗？（1）“MATLAB进行数值的表达精度与其指令窗中的数据显示精度相同。”（2）MATLAB指令窗中显示的数值有效位数不超过7位。”解答（1）否；（2）否。n 7.想要在MATLAB中产生二维数组，下面哪些指令能实现目的？（1） S=1,2,3;4,5,6;7,8

4、;9（2） S=1 2 3;4 5 6;7 8 9（3） S=1，2，3；4，5，6；7，8，9 %整个指令在中文状态下输入解答前两种输入方法可以，后一种方法不行。n 8.试为例1.3-5编写一个解题用的M脚本文件？解答直接点击新文件图标，出现M文件编辑器窗口；在该M文件编辑器中，输入例1.3-5中的全部指令；并另存为p109.m，便得到所需的脚本文件。第2章 符号运算习题2及解答n /1说出以下四条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”符号对象？ 3/7+0.1; sym(3/7+0.1); sym(3/7+0.1); vpa(sym(3/7+0.1)目的l 不能从

5、显示形式判断数据类型，而必须依靠class指令。解答c1=3/7+0.1c2=sym(3/7+0.1)c3=sym(3/7+0.1)c4=vpa(sym(3/7+0.1)Cs1=class(c1)Cs2=class(c2)Cs3=class(c3)Cs4=class(c4) c1 = 0.5286c2 =37/70c3 =0.52857142857142857142857142857143c4 =0.52857142857142857142857142857143Cs1 =doubleCs2 =symCs3 =symCs4 =sym n /2在不加专门指定的情况下，以下符号表达式中的哪一个变量

6、被认为是自由符号变量.sym(sin(w*t),sym(a*exp(-X),sym(z*exp(j*th)目的l 理解自由符号变量的确认规则。解答symvar(sym(sin(w*t),1) ans =w symvar(sym(a*exp(-X),1) ans =a symvar(sym(z*exp(j*th),1) ans =z n /3求以下两个方程的解（1）试写出求三阶方程正实根的程序。注意：只要正实根，不要出现其他根。（2）试求二阶方程在时的根。目的l 体验变量限定假设的影响解答（1）求三阶方程正实根reset(symengine)%确保下面操作不受前面指令运作的影响syms x po

7、sitivesolve(x3-44.5) ans =(2(2/3)*89(1/3)/2 （2）求五阶方程的实根syms a positive%注意：关于x的假设没有去除solve(x2-a*x+a2) Warning: Explicit solution could not be found. In solve at 83ans = empty sym syms x clearsyms a positivesolve(x2-a*x+a2) ans = a/2 + (3(1/2)*a*i)/2 a/2 - (3(1/2)*a*i)/2 n /4观察一个数（在此用记述）在以下四条不同指令作用下的异

8、同。a =, b = sym( ), c = sym( ,d ), d = sym( )在此， 分别代表具体数值 7/3 , pi/3 , pi*3(1/3) ；而异同通过vpa(abs(a-d) , vpa(abs(b-d) , vpa(abs(c-d)等来观察。目的l 理解准确符号数值的创建法。l 高精度误差的观察。解答（1）x=7/3x=7/3;a=x,b=sym(x),c=sym(x,d),d=sym(7/3), a = 2.3333b =7/3c =2.3333333333333334813630699500209d =7/3 v1=vpa(abs(a-d),v2=vpa(abs(b

9、-d),v3=vpa(abs(c-d) v1 =0.0v2 =0.0v3 =0.00000000000000014802973661668756666666667788716 （2）x=pi/3x=pi/3;a=x,b=sym(x),c=sym(x,d),d=sym(pi/3), a = 1.0472b =pi/3c =1.047197551196597631317786181171d =pi/3 v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d) v1 =0.0v2 =0.0v3 =0.0000000000000001148364282799221

10、6762806615818554 （3）x=pi*3(1/3)x=pi*3(1/3);a=x,b=sym(x),c=sym(x,d),d=sym(pi*3(1/3) a = 4.5310b =1275352044764433/281474976710656c =4.5309606547207899041040946030989d =pi*3(1/3) v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d) v1 =0.00000000000000026601114166290944374842393221638v2 =0.0000000000000002

11、6601114166290944374842393221638v3 =0.0000000000000002660111416629094726767991785515 n /5求符号矩阵的行列式值和逆，所得结果应采用“子表达式置换”简洁化。目的l 理解subexpr指令。解答A=sym(a11 a12 a13;a21 a22 a23;a31 a32 a33)DA=det(A)IA=inv(A);IAs,d=subexpr(IA,d) A = a11, a12, a13 a21, a22, a23 a31, a32, a33DA =a11*a22*a33 - a11*a23*a32 - a12*

12、a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31IAs = d*(a22*a33 - a23*a32), -d*(a12*a33 - a13*a32), d*(a12*a23 - a13*a22) -d*(a21*a33 - a23*a31), d*(a11*a33 - a13*a31), -d*(a11*a23 - a13*a21) d*(a21*a32 - a22*a31), -d*(a11*a32 - a12*a31), d*(a11*a22 - a12*a21)d =1/(a11*a22*a33 - a11*a23*a32 - a12*

13、a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31) n /6求的符号解，并进而用该符号解求，的准确值。目的l symsum, subs的应用。l 从实例中，感受指令所给出的关于符号解的含义。解答syms x kf=x(k);Z1=symsum(f,k,0,inf)Z1 =piecewise(1 = x, Inf, abs(x) 1, -1/(x - 1) %piecewise 分段函数subs(Z1,x,sym(-1/3),sym(1/pi),sym(3) ans = 3/4, -1/(1/pi - 1), Inf n /7对于，求。（提示：理

14、论结果为）目的l 符号变量的限定性定义的作用。解答syms k;x=sym(x,positive);f_k=2/(2*k+1)*(x-1)/(x+1)(2*k+1);s=simple(symsum(f_k,k,0,inf) %结果与理论值lnx相符！ s =piecewise(abs(x - 1) x + 1, log(x) %abs 绝对值注意l 解答中，条件abs(x - 1) x + 1意味着：n 约束一：x-10 此式总成立，说明“无约束”。n 情况二：-(x-1)0此为“约束”，满足题意。n /8（1）通过符号计算求的导数。（2）然后根据此结果，求和。目的l diff, limit指

15、令的应用。l 如何理解运行结果。解答syms ty=abs(sin(t)d=diff(y) %求dy/dtd0_=limit(d,t,0,left) %求dy/dt|t=0-dpi_2=limit(d,t,pi/2) %求dy/dt|t=pi/2 y =abs(sin(t)d =sign(sin(t)*cos(t)d0_ =-1dpi_2 =0 n /9求出的具有64位有效数字的积分值。目的l 符号积分的解析解和符号数值解。l 符号计算和数值计算的相互校验。解答（1）符号积分syms x clearsyms xy=exp(-abs(x)*abs(sin(x)si=vpa(int(y,-10*p

16、i,1.7*pi),64) %vpa 指定精确位数y =abs(sin(x)/exp(abs(x)si =1.087849499412904913166671875948174520895458535212845987519414166 （2）数值计算复验xx=-10*pi:pi/100:1.7*pi;sn=trapz(exp(-abs(xx).*abs(sin(xx)*pi/100 %trapz梯形法求积分sn = 1.0877 n /10计算二重积分。目的l 变上限二重积分的符号计算法。解答syms x yf=x2+y2;r=int(int(f,y,1,x2),x,1,2) r =1006

17、/105 n /11在区间，画出曲线，并计算。目的l 在符号计算中，经常遇到计算结果是特殊经典函数的情况。l 如何应用subs获得超过16位有效数字的符号数值结果。l 初步尝试ezplot指令的简便。解答（1）符号计算syms t x;f=sin(t)/t;y=int(f,t,0,x)% 将得到一个特殊经典函数y5=subs(y,x,sym(4.5)ezplot(y,0,2*pi) y =sinint(x)y5 =1.6541404143792439835039224868515 （2）数值计算复验tt=0:0.001:4.5;tt(1)=eps;yn=trapz(sin(tt)./tt)*0

18、.001 yn = 1.6541 n /12在的限制下，求的一般积分表达式，并计算的32位有效数字表达。目的l 一般符号解与高精度符号数值解。解答syms xsyms n positivef=sin(x)n;yn=int(f,x,0,pi/2) y3s=vpa(subs(yn,n,sym(1/3)y3d=vpa(subs(yn,n,1/3) yn =beta(1/2, n/2 + 1/2)/2y3s =1.2935547796148952674767575125656y3d =1.2935547796148951782413405453553 n 13.有序列，（在此，），求这两个序列的卷积。

19、目的l 符号离散卷积直接法和变换法。解答（1）直接法syms a b k nx=ak;h=bk;w=symsum(subs(h,k,n)*subs(x,k,k-n),n,0,k)%据定义y1=simple(w)w =piecewise(a = b, bk + bk*k, a b, (a*ak - b*bk)/(a - b)y1 =piecewise(a = b, bk + bk*k, a b, (a*ak - b*bk)/(a - b) （2）变换法（复验）syms zX=ztrans(ak,k,z);H=ztrans(bk,k,z);y2=iztrans(H*X,z,k)%通过Z变换及反变换

20、 y2 =piecewise(b 0, (a*ak)/(a - b) - (b*bk)/(a - b) 说明l 符号计算不同途径产生的结果在形式上有可能不同，而且往往无法依靠符号计算本身的指令是它们一致。此时，必须通过手工解决。n 14.设系统的冲激响应为，求该系统在输入，作用下的输出。目的l 符号连续函数卷积的直接法和变换法。l 符号变量限定性定义的作用。l laplace, ilaplace指令的应用。解答（1）直接法syms th=exp(-3*t);u=cos(t);syms tao;h_tao=subs(h,t,tao);u_t_tao=subs(u,t,t-tao);hu_tao=

21、h_tao*u_t_tao;hut=simple(int(hu_tao,tao,0,t)%直接卷积hut =(3*cos(t)/10 - 3/(10*exp(3*t) + sin(t)/10 （2）变换法（复验）syms s;HU=laplace(h,t,s)*laplace(u,t,s);huL=simple(ilaplace(HU,s,t) %拉氏变换及反变换 huL =(3*cos(t)/10 - 3/(10*exp(3*t) + sin(t)/10 n 15.求的Fourier变换。目的l 符号变量限定性定义的作用。l fourier指令的应用。解答syms A t wa=sym(a,

22、positive);f=A*exp(-a*abs(t);y=fourier(f,t,w)F=simple(y) y =(2*A*a)/(a2 + w2)F =(2*A*a)/(a2 + w2) n 16.求的Fourier变换，并画出时的幅频谱。目的l 单位阶跃符号函数heaviside的应用。l subs实现多变量置换。l ezplot的使用。解答syms t A w;tao=sym(tao,positive);f=A*(1+t/tao)*(heaviside(t+tao)-heaviside(t)+(1-t/tao)*(heaviside(t)-heaviside(t-tao);Fw=fo

23、urier(f,t,w);Fws=simple(Fw)Fw2=subs(Fws,A,tao,2,2)ezplot(abs(Fw2)grid Fws =-(4*A*(cos(tao*w)/2)2 - 1)/(tao*w2)Fw2 =-(8*cos(w)2 - 8)/(2*w2) n 17.求的Laplace反变换。解答syms s t F=(s+3)/(s3+3*s2+6*s+4);f=simple(ilaplace(F,s,t) f =(3(1/2)*sin(3(1/2)*t) - 2*cos(3(1/2)*t) + 2)/(3*exp(t) n 18.利用符号运算证明Laplace变换的时域

24、求导性质：。目的l 符号计算用于定理证明。解答syms t s;y=sym(f(t);df=diff(y,t);Ldy=laplace(df,t,s) Ldy =s*laplace(f(t), t, s) - f(0) n 19.求的Z变换表达式。目的l 注意：变换中，被变换变量的约定。解答syms lambda k T z;f_k=k*exp(-lambda*k*T);F_z=simple(ztrans(f_k,k,z) F_z =(z*exp(T*lambda)/(z*exp(T*lambda) - 1)2 n 20.求方程的解。目的l solve指令中，被解方程的正确书写，输出量的正确次

25、序。解答eq1=x2+y2=1;eq2=x*y=2;x,y=solve(eq1,eq2,x,y) x = (1/2 + (15(1/2)*i)/2)(1/2)/2 - (1/2 + (15(1/2)*i)/2)(3/2)/2 - (1/2 + (15(1/2)*i)/2)(1/2)/2 + (1/2 + (15(1/2)*i)/2)(3/2)/2 (1/2 - (15(1/2)*i)/2)(1/2)/2 - (1/2 - (15(1/2)*i)/2)(3/2)/2 - (1/2 - (15(1/2)*i)/2)(1/2)/2 + (1/2 - (15(1/2)*i)/2)(3/2)/2y =

26、(1/2 + (15(1/2)*i)/2)(1/2) -(1/2 + (15(1/2)*i)/2)(1/2) (1/2 - (15(1/2)*i)/2)(1/2) -(1/2 - (15(1/2)*i)/2)(1/2) n 21.求图p2-1所示信号流图的系统传递函数，并对照胡寿松主编“自动控制原理”中的例2-21结果，进行局部性验证。图p2-1目的l 理解和掌握信号流图传递函数的“代数状态方程解法”。l 并设法用胡寿松主编的“自动控制原理”的例2-21进行局部性验证。解答（1）求传递函数syms G1 G2 G3 G4 G5 G6 G7 H1 H2 H3 H4 H5A=0 0 0 0 -H3

27、 -H4; G1 0 -H1 0 0 0;0 G2 0 0 -H2 G6; 0 0 G3 0 0 G7; 0 0 0 G4 0 0; 0 G5 0 0 0 -H5;b= 1; 0; 0; 0; 0; 0;c= 0 0 0 0 1 0;Y2U=c*(eye(size(A)-A)b); %求传递函数NN,DD=numden(Y2U);%分离出分子、分母多项式DD=sort(DD);%分母多项式排序disp(blanks(5),传递函数 Y2U 为)pretty(NN/DD) 传递函数 Y2U 为 (G1 G4 (G2 G3 + G5 G7 + G3 G5 G6 + G2 G3 H5) / (H5

28、+ G2 H1 + G3 G4 H2 + G1 G5 H4 + G5 G6 H1 + G2 H1 H5 + G3 G4 H2 H5 + G1 G2 G3 G4 H3 + G1 G4 G5 G7 H3 - G4 G5 G7 H1 H2 + G1 G3 G4 G5 G6 H3 + G1 G2 G3 G4 H3 H5 + G1 G3 G4 G5 H2 H4 + 1) （2）局部性验证syms a b c d e f gy2u=subs(Y2U,G1,G2,G3,G4,G5,G6,G7,H1,H2,H3,H4,H5,a,e,f,1,b,c,0,g,0,0,0,d);nn,dd=numden(y2u);

29、dd=sort(dd);disp(blanks(5),局部性验证用的传递函数y2u)pretty(nn/dd) 局部性验证用的传递函数y2u a (e f + b c f + d e f) - d + e g + b c g + d e g + 1 此结果与胡寿松主编的“自动控制原理”例2-21一致。n 22.采用代数状态方程法求图p2-2所示结构框图的传递函数和。图p2-2目的l 运用“代数状态方程解法”求输入和扰动同时存在的结构框图的传递函数。解答（1）理论演绎对于结构框图写出状态方程（p2-1）此式第一个方程关于x的解可写为 （p2-2）把此式代入式（p2-1）的第二个方程，加以整理后可

30、得据此可写出传递函数 （p2-3） （p2-4）（2）列出“元素级”状态方程值得提醒：在编写M码之前，最好先在草稿纸上，仔细“元素级”状态方程是避免出错的冲要措施。对此，不要掉以轻心。本例的“元素级”状态方程如下 （p2-5）（3）编写相应的M码syms G1 G2 G3 H1 H2A=0 0 0 -G1 -G1; G2 0 -G2 0 0;0 0 0 0 0; 0 H1 0 0 0; 0 H2 0 0 0;b= G1; 0; 0; 0; 0;f= 0; 0; G3; 0; -H2;c= 0 1 0 0 0;d=0;g=-1;R=c/(eye(size(A)-A); %中间变量Y2U=R*b+

31、d;%计算传递函数 Y/U Y2W=R*f+g;%计算传递函数 Y/WNU,DU=numden(Y2U);%分离出分子、分母多项式DU=sort(DU);%分母多项式排序disp(blanks(5),传递函数 Y2U 为)pretty(NU/DU) NW,DW=numden(Y2W);NW=sort(NW);DW=sort(DW);disp(blanks(5),传递函数 Y2W 为)pretty(NW/DW) 传递函数 Y2U 为 G1 G2 - G1 G2 H1 + G1 G2 H2 + 1 传递函数 Y2W 为 G2 G3 + G1 G2 H1 + 1 - - G1 G2 H1 + G1

32、G2 H2 + 1 n 23.求微分方程的通解，并绘制任意常数为1时解的图形。目的l 理解指令dsolve的正确使用。l 对dsolve输出结果的正确理解。l ezplot指令绘图时，如何进行线色控制。l 如何覆盖那些不能反映图形窗内容的图名。解答（1）求通解reset(symengine)clearsyms y xy=dsolve(0.2*y*Dy+0.25*x=0,x) y = 2(1/2)*(C3 - (5*x2)/8)(1/2) -2(1/2)*(C3 - (5*x2)/8)(1/2) （2）根据所得通解中不定常数的符号写出“对其进行数值替代的指令”yy=subs(y,C3,1) %将

33、通解中的C3用1代替 yy = 2(1/2)*(1 - (5*x2)/8)(1/2) -2(1/2)*(1 - (5*x2)/8)(1/2) （3）观察通解中两个分解的平方是否相同yy(1)2=yy(2)2 ans = 1 （4）于是可考虑函数的平方关系syms Yfxy=Y2-yy(1)2 fxy =Y2 + (5*x2)/4 - 2 （5）根据平方关系式画完整曲线clfezplot(fxy,-2,2,-2,2)axis squaregrid on （6）假如直接用“分解”画曲线，那么将是不完整的 ezplot(yy(1),hold oncc=get(gca,Children);set(cc

34、,Color,r)ezplot(yy(2),axis(-2 2 -2 2)legend(y(1),y(2),hold off;title( )%覆盖不完全的图名gridaxis square n 24.求一阶微分方程的解。目的l 初值微分方程的符号解。l pretty指令的使用。解答x=dsolve(Dx=a*t2+b*t,x(0)=2,t)pretty(x)%比较易读的表达形式 x =(t2*(3*b + 2*a*t)/6 + 2 2 t (3 b + 2 a t) - + 2 6 n 25.求边值问题的解。（注意：相应的数值解法比较复杂）。目的l 边值微分方程的符号解。解答f,g=dsol

35、ve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1) f =sin(4*t)*exp(3*t)g =cos(4*t)*exp(3*t) 第3章 数值数组及其运算习题3及解答n 5.要求在闭区间上产生具有10个等距采样点的一维数组。试用两种不同的指令实现。目的l 数值计算中产生自变量采样点的两个常用指令的异同。解答%方法一 t1=linspace(0,2*pi,10)%方法二t2=0:2*pi/9:2*pi %要注意采样间距的选择，如这里的2*pi/9. t1 = Columns 1 through 7 0 0.6981 1.3963 2.0944 2.7925 3.

36、4907 4.1888 Columns 8 through 10 4.8869 5.5851 6.2832t2 = Columns 1 through 7 0 0.6981 1.3963 2.0944 2.7925 3.4907 4.1888 Columns 8 through 10 4.8869 5.5851 6.2832 n 6.由指令rng(default),A=rand(3,5)生成二维数组A，试求该数组中所有大于0.5的元素的位置，分别求出它们的“全下标”和“单下标”。目的l 数组下标的不同描述：全下标和单下标。l sub2ind, int2str, disp的使用。l 随机发生器的

37、状态控制：保证随机数的可复现性。解答rng(default)A=rand(3,5)ri,cj=find(A0.5);id=sub2ind(size(A),ri,cj);ri=ri;cj=cj;disp( )disp(大于0.5的元素的全下标)disp(行号 ,int2str(ri)disp(列号 ,int2str(cj)disp( )disp(大于0.5的元素的单下标)disp(id) A = 0.8147 0.9134 0.2785 0.9649 0.9572 0.9058 0.6324 0.5469 0.1576 0.4854 0.1270 0.0975 0.9575 0.9706 0.8

38、003 大于0.5的元素的全下标行号 1 2 1 2 2 3 1 3 1 3列号 1 1 2 2 3 3 4 4 5 5 大于0.5的元素的单下标 1 2 4 5 8 9 10 12 13 15 n 7.采用默认全局随机流，写出产生长度为1000的“等概率双位（即取-1，+1）取值的随机码”程序指令，并给出 -1码的数目。目的l 两种基本随机发生器的使用。l 关系运算产生逻辑数组可用于数组的元素的标识和寻访。l 逻辑数组的应用。l 如何判断两个整数数组是否相等。解答（1）运用均匀随机数解题法解法1rng default%为以下结果重现而设；产生默认随机流。详见第4.3.2节A=rand(1,1

39、000);a=2*(A0.5)-1;Na=sum(a=-1) Na = 512 （2）运用正态随机数解题法解法2 randn(state,123)B=randn(1,1000);b=2*(B0)-1;Nb=sum(b=-1) Nb = 462 （3）直接发生法解法3c=randsrc(1,1000,-1,1);Nc=sum(c=-1) Nc = 482 n 2.已知矩阵，运行指令B1=A.(0.5), B2=A(0.5), 可以观察到不同运算方法所得结果不同。（1）请分别写出根据B1, B2恢复原矩阵A的程序。（2）用指令检验所得的两个恢复矩阵是否相等。目的l 数组运算和矩阵运算的不同。l 如

40、何判断两个双精度数组是否相等。l norm指令的应用。解答A=1,2;3,4;B1=A.0.5B2=A0.5A1=B1.*B1;A2=B2*B2;norm(A1-A2,fro)% 求误差矩阵的F-范数，当接近eps量级时，就认为实际相等B1 = 1.0000 1.4142 1.7321 2.0000B2 = 0.5537 + 0.4644i 0.8070 - 0.2124i 1.2104 - 0.3186i 1.7641 + 0.1458ians = 8.4961e-016 n 4.在时间区间 0,10中，绘制曲线。要求分别采取“标量循环运算法”和“数组运算法”编写两段程序绘图。 目的l 加强

41、理解数组运算的机理和应用。l 初步使用subplot, plot, xlabel, ylabel等指令绘图。解答%标量循环运算法t=linspace(0,10,200);N=length(t);y1=zeros(size(t);for k=1:Ny1(k)=1-exp(-0.5*t(k)*cos(2*t(k);endsubplot(1,2,1),plot(t,y1),xlabel(t),ylabel(y1),grid on%数组运算法y2=1-exp(-0.5*t).*cos(2*t);subplot(1,2,2),plot(t,y2),xlabel(t),ylabel(y2),grid on

42、 n 8.先运行clear,format long,rng(default),A=rand(3,3)，然后根据A写出两个矩阵：一个对角阵B，其相应元素由A的对角元素构成；另一个矩阵C，其对角元素全为0，而其余元素与对应的A阵元素相同。目的l 常用指令diag的使用场合。解答clear,format longrng(default)A=rand(3,3)B=diag(diag(A)C=A-B A = 0.814723686393179 0.913375856139019 0.278498218867048 0.905791937075619 0.632359246225410 0.5468815

43、19204984 0.126986816293506 0.097540404999410 0.957506835434298B = 0.814723686393179 0 0 0 0.632359246225410 0 0 0 0.957506835434298C = 0 0.913375856139019 0.278498218867048 0.905791937075619 0 0.546881519204984 0.126986816293506 0.097540404999410 0 n 9.先运行指令x=-3*pi:pi/15:3*pi; y=x; X,Y=meshgrid(x,y)

44、; warning off; Z=sin(X).*sin(Y)./X./Y; 产生矩阵Z。（1）请问矩阵Z中有多少个“非数”数据？（2）用指令surf(X,Y,Z); shading interp观察所绘的图形。（3）请写出绘制相应的“无裂缝”图形的全部指令。目的l 初步感受三维曲面的绘制方法。l 非数NaN的产生，非数的检测，和对图形的影响。l sum的应用。l eps如何克服“被零除”的尴尬。解答x=-3*pi:pi/15:3*pi;y=x;X,Y=meshgrid(x,y);warning offZ=sin(X).*sin(Y)./X./Y;NumOfNaN=sum(sum(isnan(

45、Z)%计算“非数”数目subplot(1,2,1),surf(X,Y,Z),shading interp,title(有缝图)%产生无缝图XX=X+(X=0)*eps;YY=Y+(Y=0)*eps;ZZ=sin(XX).*sin(YY)./XX./YY;subplot(1,2,2),surf(XX,YY,ZZ),shading interp,title(无缝图) NumOfNaN = 181 n 10.下面有一段程序，企图用来解决如下计算任务：有矩阵，当依次取10, 9, 8, 7, 6, 5, 4, 3, 2, 1时，计算矩阵“各列元素的和”，并把此求和结果存放为矩阵Sa的第k行。例如时，A

46、阵为，此时它各列元素 的和是一个行数组，并把它保存为Sa的第3行。问题：该段程序的计算结果对吗？假如计算结果不正确，请指出错误发生的根源，并改正之。目的l 正确理解sum的工作机理。l reshape的应用。解答（1）企图用以下程序完成题目要求。for k=10:-1:1A=reshape(1:10*k,k,10);Sa(k,:)=sum(A);endSa Sa = 55 55 55 55 55 55 55 55 55 55 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 122 13

47、8 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 855 955 （2）正确性分析除k=1外，计算所得Sa所有行的结果都正确。但k=1时，Sa的第1行应该与相同。上述程序的错误是对sum理解不正确。

48、sum对二维数组，求和按列施行；而对一维数组，不管行数组或列数组，总是求那数组所有元素的和。正确的程序应该写成for k=10:-1:1A=reshape(1:10*k,k,10);Sa(k,:)=sum(A);if k=1Sa(k,:)=A;endendSa Sa = 1 2 3 4 5 6 7 8 9 10 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 122 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165

49、 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 855 955 第4章 数值运算习题 4 及解答n 1.根据题给的模拟实际测量数据的一组和 试用数值差分diff或数值梯度gradient指令计算，然后把和曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从prob_data401.mat获得）目的l 强调：要

50、非常慎用数值导数计算。l 练习mat数据文件中数据的获取。l 实验数据求导的后果l 把两条曲线绘制在同一图上的一种方法。解答（1）从数据文件获得数据的指令假如prob_data401.mat文件在当前目录或搜索路径上clearload prob_data401.mat （2）用diff求导的指令dt=t(2)-t(1);yc=diff(y)/dt;%注意yc的长度将比y短1plot(t,y,b,t(2:end),yc,r)grid on （3）用gradent求导的指令（图形与上相似）dt=t(2)-t(1);yc=gradient(y)/dt;plot(t,y,b,t,yc,r)grid o

51、n 说明l 不到万不得已，不要进行数值求导。l 假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以上。l 求导会使数据中原有的噪声放大。n 2.采用数值计算方法，画出在区间曲线，并计算。提示l 指定区间内的积分函数可用cumtrapz指令给出。l 在计算要求不太高的地方可用find指令算得。目的l 指定区间内的积分函数的数值计算法和cumtrapz指令。l find指令的应用。解答dt=1e-4;t=0:dt:10;t=t+(t=0)*eps;f=sin(t)./t;s=cumtrapz(f)*dt;plot(t,s,LineWidth,3)ii=find(t=4.5);s45=s(ii) s45 = 1.6541 n 3.求函数的数值积分，并请采用符号计算尝试复算。提示l 数值积分均可尝试。l 符号积分的局限性。目的l 符号积分的局限性。解答dx=pi/2000;x=0:dx:pi;s=trapz(exp(sin(x)