《[理学]matlab教程2011a张志涌-课后答案.doc》由会员分享,可在线阅读,更多相关《[理学]matlab教程2011a张志涌-课后答案.doc(38页珍藏版)》请在三一文库上搜索。
1、第二章 %分段函数的写法syms y xz=int(x*y,x,0,1);g = evalin(symengine, piecewise(y 1/2, char(z) , y = 1/2, 1)G=int(g,y)g = piecewise(1/2 y, y/2, y = 1/2, 1) G = piecewise(1/2 y, y2/4, y 0.5)subindex_A=sub2ind(size(A),I1,J1)subindex_A=find(A0.5)I,J=ind2sub(size(A),subindex_A)index_A=I J%122_2rand(twister,0),A=ra
2、nd(3,5),B=A0.5,C=find(B)ci,cj=ind2sub(size(A),C) %此法太繁Ci,Cj=find(B)%122_5t=linspace(1,10,100) (1) y=1-exp(-0.5.*t).*cos(2.*t),plot(t,y)(2) L=length(t)for k=1:L, yy(k)=1-exp(-0.5*t(k)*cos(2*t(k), end, plot(t,yy)%122_6clear,format long, rand(twister,1),A=rand(3,3),B=diag(diag(A),C=A.*(B)%或C=A-B%或C=tri
3、u(A,1)+tril(A,-1)%习题3-3% s=sign(randint(1,1000,123)-.5);% n=sum(s=-1)rand(state,123),s=sign(rand(1,1000)-.5),n=sum(s=-1)%习题3-4clear allA=1 2;3 4B1=A.(0.5), B2=A(0.5),A1=B1.2A2=B22A1-Anorm(A1-A)A1-A2norm(A1-A2)clear allA=sym(1 2;3 4)B1=A.(0.5), B2=A(0.5),A1=simple(B1.2)A2=simple(B22)A1-Avpa(A1-A)A1-A
4、2vpa(A1-A2)%习题3-5t=linspace(1,10,100) %(1)L=length(t)for k=1:L yy(k)=1-exp(-0.5*t(k)*cos(2*t(k) endfigure(1),plot(t,yy)%(2)y=1-exp(-0.5.*t).*cos(2.*t),figure(2),plot(t,y) figure(3),subplot(1,2,1),plot(t,y)subplot(1,2,2),plot(t,yy)%习题3-6clear,format long, rand(twister,1),A=rand(3,3),B=diag(diag(A),C=
5、A.*(B)%或C=A-B%或C=triu(A,1)+tril(A,-1)%习题3-7clear,x=-3*pi:pi/15:3*pi; y=x; X,Y=meshgrid(x,y); % X=ones(size(y)*x;Y=y*ones(size(x);warning off; Z=sin(X).*sin(Y)./X./Y;% (1)“非数”数目m=sum(sum(isnan(Z); %k=Z(isnan(Z);m=length(k)% (2)绘图surf(X,Y,Z); shading interp% (3)“无裂缝”图形的全部指令:x=-3*pi:pi/15:3*pi;y=x;X,Y=
6、meshgrid(x,y); % X=ones(size(y)*x;Y=y*ones(size(x);X=X+(X=0)*eps;Y=Y+(Y=0)*eps;Z=sin(X).*sin(Y)./X./Y;surf(X,Y,Z);shading interp%习题3-8clearfor k=10:-1:1; A=reshape(1:10*k,k,10); Sa(k,:)=sum(A,1);endSa clearfor k=10:-1:1; A=reshape(1:10*k,k,10); if k=1 Sa(k,:)=A; else Sa(k,:)=sum(A); endendSa 第四章%习题4
7、-1clear allload prob_401;diff_y=diff(y)./diff(t);gradient_y=gradient(y)./gradient(t);% plot(t(1:end-1),diff_y,t,gradient_y)figure(1)subplot(1,2,1),plot(t,y,t(1:end-1),diff_y)subplot(1,2,2),plot(t,y,t,gradient_y)%上面结果不好因自变量 采样间距太小,将间距增大后结果较好N=20diff_y1=(diff(y(1:N:end)./diff(t(1:N:end);gradient_y1=(g
8、radient(y(1:N:end)./gradient(t(1:N:end);% plot(t(1:end-1),diff_y,t,gradient_y)t1=t(1:N:end);length(t1)figure(2)subplot(1,2,1),plot(t,y,t1(1:end-1),diff_y1)subplot(1,2,2),plot(t,y,t1,gradient_y1)%习题4-2d=0.5; tt=0:d:10; t=tt+(tt=0)*eps; y=sin(t)./t; s=d*trapz(y) % 计算出积分值ss=d*(cumtrapz(y) %计算梯形法累计积分并绘积
9、 分曲线plot(t,y,t,ss,r),hold on%用find指令计算y(4.5),并绘出该点y4_5=ss(find(t=4.5)%插值法计算y(4.5),并绘出该点yi=interp1(t,ss,4.5), plot(4.5,yi,r+)% yi=interp1(t,ss,4.2 4.3 4.5), plot(4.2 4.3 4.5,yi,r+)% clf%用精度可控指令quad 即Simpson法,quadl 即Lobatto法计算y(4.5)yy=quad(sin(t)./t,0,4.5)yy=quadl(sin(t)./t,0,4.5)warning off % 取消警告性提示
10、时用%此法可用,但有警告性提示,取消提示加 warning offcleartt=0:0.1:10; warning offfor i=1:101 q(i)=quad(sin(t)./t,0,tt(i);endplot(tt, q)y=quad(sin(t)./t,0,4.5)%符号解法,匿名函数求y(4.5)f=(x)(int(sin(t)/t,0,x),vpa(f(4.5)%符号解法syms x t y1 y2 y1i, y1=sin(t)./t, y1i=int(y1,t,0,x), y2=subs(y1i,x,4.5)hold on ,plot(4.5,y2,*m)%习题4-3clea
11、r alld=pi/20; x=0:d:pi; fx=exp(sin(x).3); s=d*trapz(fx) % 梯形法计算积分值%用精度可控指令quad 即Simpson法,quadl 即Lobatto法计算积分s1=quad(exp(sin(x).3),0,pi)s2=quadl(exp(sin(x).3),0,pi)%符号计算解法% s3=vpa(int(exp(sin(x)3),x,0,pi)s3=vpa(int(exp(sin(x)3),0,pi)s4=vpa(int(sym(exp(sin(x)3),0,pi)%习题4-4clear all%用精度可控指令quad 即Simpso
12、n法,quadl 即Lobatto法计算积分s1=quad(exp(-abs(x).*abs(sin(x),-5*pi,1.7*pi,1e-10)s2=quadl(exp(-abs(x).*abs(sin(x),-5*pi,1.7*pi)%符号计算解法在Matlab6.5版本上下式计算结果多加了值1% s3=vpa(int(exp(-abs(x)*abs(sin(x),-5*pi,1.7*pi)syms x;s3=vpa(int(exp(-abs(x)*abs(sin(x),-5*pi,1.7*pi)% s3=vpa(int(sin(x),-pi/2,pi/2)% 梯形法计算积分值d=pi/1
13、000; x=-5*pi:d:1.7*pi; fx=exp(-abs(x).*abs(sin(x); s=d*trapz(fx)%习题4-5clear allx1=-5;x2=5;%采用内联函数或字符串函数求极值yx=inline(sin(5*t).2.*exp(0.06*t.2)-1.5.*t.*cos(2*t)+1.8.*abs(t+0.5)xn0,fval=fminbnd(yx,x1,x2)%绘函数图并标出最小值点t=x1:0.1:x2; plot(t,yx(t),hold on ,plot(xn0,fval,r*)%字符串函数求极值xn0,fval=fminbnd(sin(5.*x).
14、2.*exp(0.06.*x.2)-1.5.*x.*cos(2.*x)+1.8.*abs(x+0.5),-5,5)% yx=(sin(5.*x).2.*exp(0.06.*x.2)-1.5.*x.*cos(2.*x)+1.8.*abs(x+0.5)% xn0,fval=fminbnd(yx,x1,x2)下一条指令即字符串函数在无论在2010a版本还是Matlab6.5版本上不适用,因为对字符串函数只是别符号x,不识别t% xn0,fval=fminbnd(sin(5.*t).2.*exp(0.06.*t.2)-1.5.*t.*cos(2.*t)+1.8.*abs(t+0.5),-5,5)下一条
15、指令即匿名函数在Matlab6.5版本上不适用% yx=(t)(sin(5*t).2.*exp(0.06*t.2)-1.5.*t.*cos(2*t)+1.8.*abs(t+0.5) %本条指令即匿名函数在Matlab6.5版本上不适用% yx=(t)(sin(5.*x).2.*exp(0.06.*x.2)-1.5.*x.*cos(2.*x)+1.8.*abs(x+0.5) %本条指令即匿名函数在Matlab6.5版本上不适用% xn0,fval=fminbnd(yx,x1,x2)%习题4-6tspan=0,0.5;%y0=1;0;%tt,yy=ode45(DyDt_6,tspan,y0);y0
16、_5=yy(end,1)% figure(1)% plot(tt,yy(:,1)% xlabel(t),title(x(t) % figure(2)% plot(yy(:,1),yy(:,2)% xlabel(位移),ylabel(速度)函数DyDt_6另存为一个文件function ydot=DyDt_6(t,y)mu=3;ydot=y(2);mu*y(2)-2*y(1)+1;%符号计算解法S = dsolve(D2y-3*Dy+2*y = 1,y(0) = 1,Dy(0) = 0)ys0_5=subs(S,0.5)%习题4-7clear allA=magic(8)B=orth(A)rref
17、(A)rref(B)A = 64 2 3 61 60 6 7 57 9 55 54 12 13 51 50 16 17 47 46 20 21 43 42 24 40 26 27 37 36 30 31 33 32 34 35 29 28 38 39 25 41 23 22 44 45 19 18 48 49 15 14 52 53 11 10 56 8 58 59 5 4 62 63 1B = -0.35355339059327 0.54006172486732 0.35355339059327 -0.35355339059327 -0.38575837490523 -0.353553390
18、59327 -0.35355339059327 -0.23145502494314 -0.35355339059327 -0.35355339059327 0.07715167498105 0.35355339059327 -0.35355339059327 -0.07715167498105 0.35355339059327 -0.35355339059327 0.23145502494314 -0.35355339059327 -0.35355339059327 0.38575837490523 -0.35355339059327 -0.35355339059327 -0.54006172
19、486732 0.35355339059327ans = 1 0 0 1 1 0 0 1 0 1 0 3 4 -3 -4 7 0 0 1 -3 -4 4 5 -7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0%习题4-8A = sym(gallery(5)v, lambda = eig(A)cond(A)clear all, A=gallery(5)V,D,s=condeig(A)V,D=eig(A)cond(A)jordan(A)ans = 0 1 0 0 0 0 0 1 0 0
20、 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0%习题4-9clear all,A=magic(3)b=ones(3,1)x=Abx = 0.06666666666667 0.06666666666667 0.06666666666667 x=inv(A)*b rref(A,b)ans = 1.00000000000000 0 0 0.06666666666667 0 1.00000000000000 0 0.06666666666667 0 0 1.00000000000000 0.06666666666667%习题4-10解不唯一 clear all,A=magic(4)b=o
21、nes(4,1)rref(A,b)ans = 1.00000000000000 0 0 1.00000000000000 0.05882352941176 0 1.00000000000000 0 3.00000000000000 0.11764705882353 0 0 1.00000000000000 -3.00000000000000 -0.05882352941176 0 0 0 0 0x=Abx = 0.05882352941176 0.11764705882353 -0.05882352941176 0 xg=null(A) xg = 0.22360679774998 0.6708
22、2039324994 -0.67082039324994 -0.22360679774998x+xg 也是方程的解ans = 0.28243032716174 0.78846745207347 -0.72964392266170 -0.22360679774998下面逆阵法的解不正确x=inv(A)*b x = 0.03125000000000 0.12500000000000 -0.06250000000000 -0.01562500000000因为A*xans = 0.35937500000000 0.78125000000000 0.59375000000000 0.9218750000
23、0000%习题4-11clear all,A=magic(4)b=(1:4)rref(A,b)ans = 1 0 0 1 0 0 1 0 3 0 0 0 1 -3 0 0 0 0 0 1x=Abx = 1.0e+015 * -0.56294995342131 -1.68884986026394 1.68884986026394 0.56294995342131 A*xans = 0 6.50000000000000 4.00000000000000 7.81250000000000x=inv(A)*b x = 1.0e+015 * -0.56294995342131 -1.6888498602
24、6394 1.68884986026394 0.56294995342131没有准确解,以下是伪逆阵求得的最小二乘解 x=pinv(A)*bx = 0.02352941176471 0.12352941176471 0.12352941176471 0.02352941176471A*xans = 1.30000000000000 2.90000000000000 2.10000000000000 3.70000000000000 x=lsqnonneg(A,b) x = 0.04705882352941 0.19411764705882 0.05294117647059 0 A*xans =
25、 1.30000000000000 2.90000000000000 2.10000000000000 3.70000000000000%习题4-12clear all,y_C=inline(-0.5+t-10.*exp(-0.2.*t).*abs(sin(sin(t),t); t=-10:0.01:10;Y=y_C(t); plot(t,Y,r); hold onplot(t,zeros(size(t),k);xlabel(t);ylabel(y(t) zoom ontt,yy=ginput(1),zoom offt1,y1=fzero(y_C,tt) t1 = 2.734117322080
26、59y1 = -4.440892098500626e-016t2,y2=fsolve(y_C,tt) t2 = 2.73411732208069y2 = 6.279421427279885e-013%习题4-13%符号计算解法S=solve(sin(x-y)=0,cos(x+y)=0,x,y)S.x, S.y ans = 1/4*pi -1/4*pi ans = 1/4*pi -1/4*pi %数值计算解法x=-2:0.05: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) h
27、old on,contour(X,Y,F2,v),hold off x0,y0=ginput(2); disp(x0,y0) % -0.77880184331797 -0.77777777777778% 0.79723502304147 0.77777777777778 fun=sin(x(1)-x(2),cos(x(1)+x(2); xy,f,exit=fsolve(fun,x0(2),y0(2) % xy =% % 0.78539816339745 0.78539816339745% % % f =% % 1.0e-016 *% % 0 0.61232339957368% % % exit
28、 =% % 1或者采用下面的程序x0 = -1/4*pi; -1/4*pi; % Make a starting guess at the solutionoptions=optimset(Display,iter); % Option to display outputx,fval = fsolve(Myfsove,x0,options) % Call solverx = -5.49778714378216 -5.49778714378216fval = 1.0e-013 * 0 0.43980294605305x0 = -pi/10; -pi/2; % Make a starting gu
29、ess at the solutionoptions=optimset(Display,iter); % Option to display outputx,fval = fsolve(Myfsove,x0,options) % Call solver函数Myfsove另存为一个文件function F=Myfsove(x)F = sin(x(1)-x(2); cos(x(1)+x(2); %习题4-14y = binopdf(0:100,100,0.157);stem(1:length(y),y),axis(0 length(y) 0 .12 ) %习题4-15a=normrnd(4,2,1
30、0000,1);hist(a)histfit(a) hist(a,sqrt(10000)figure(2)subplot(3,1,1),hist(a)subplot(3,1,2),histfit(a) subplot(3,1,3),hist(a,sqrt(10000) %习题4-16 load prob_416; Mx=max(max(R), Me=mean(mean(R), St=std(R(:), %习题4-17format rat, NX=conv(3,0,1,0,1,0,0,0.5), DX=conv(1,2,-2,5,2,0,1)q,r=deconv(NX,DX), cq=商多项式为
31、 ;cr=余多项式为 ;disp(cq,poly2str(q,s),disp(cr,poly2str(r,s)%验算qp2=conv(q,DX), pp1=qp2+r, pp1=NX%或验算norm( pp1-NX,fro)%习题4-18load prob_418, who,xP=polyfit(x,y,5), Pt=poly2str(P,t)xx=-1:0.01:4, yy=polyval(P,xx), plot(xx,yy,x,y,*r)legend(拟合曲线,原始曲线,Location,SouthEast)%解法二(拟合的基本思想)x0=x;y0=y;m=length(x);n=5;X=zeros(m,n+1);for k=1:n X(:,n-k+1)=(x0.k);endX(:,n+1)=ones(m,1);aT=(Xy0),norm(P-aT)%习题4-18h=0.05,0.24,0.40,0.24,0.15,-0.1,0.1 randn(state,1);u=2*(randn(1,100)0.5)-1;y=conv(u,h);figure(2)subplot(2,1,1),stem(u,filled)axis(
链接地址:https://www.31doc.com/p-1986551.html