matlab教程2011a（张志涌）-课后答案

第二章

%分段函数的写法 syms y x z=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, y^2/4], [y <= 1/2, y]) %习题2-1

class(vpa(sym(3/7+0.1),4)) %习题2-2

findsym(sym('sin(w*t)'),1) findsym(sym('a*exp(-X)'),1) findsym(sym('z*exp(j*th)'),1) %习题2-3 syms x a positive x_1=solve('x^3-44.5') syms x unreal

x_2=solve('x^2-a*x+a^2') syms a unreal

x_2=solve('x^2-a*x+a^2') x_3=solve('x^2-a*x-a^2') syms x real

evalin(symengine,'anames(Properties)') evalin(symengine,'getprop(x)')

%习题2-4

7/3, pi/3 , pi*3^(1/3) a = pi*3^(1/3) , b = sym( pi*3^(1/3) ), c = sym( pi*3^(1/3),'d'), d = sym( 'pi*3^(1/3) ' )

vpa(abs(a-d)) , vpa(abs(b-d)) , vpa(abs(c-d)) %习题2-5

syms a11 a12 a13 a21 a22 a23 a31 a32 a33 A = [a11 a12 a13;a21 a22 a23;a31 a32 a33] a=det(A) B=inv(A) C=subexpr(B) [RS,w]=subexpr(B,'w') %习题2-6 syms k syms a positive

% fk =a^k * heaviside(k) fk =a^k

s=symsum(fk,k,0,inf) %习题2-7 clear all syms k syms x positive

fk =2/(2*k+1)*(((x-1)/(x+1))^(2*k+1)) s=symsum(fk,k,0,inf) s1=simple(s) %习题2-8 clear all, syms t

y=abs(sin(t)) df=diff(y),class(df) df1=limit(df,t,0,'left') df2=subs(df,'t',sym(pi/2)) %习题2-9 clear all, syms x;

f=exp(-abs(x)).*abs(sin(x));

fint=int(f,x,-5*pi,1.7*pi), digits(64), vpa(fint) class(fint) %习题2-10

clear all,syms x y,f=x.^2+y.^2,

fint=(int(int(f,y,1,x.^2),x,1,2)), double(fint) %习题2-11

clear all, syms t x; f=sin(t)/t, yx=int(f,t,0,x), ezplot(yx,[0 2*pi]) fint=subs(yx,x,4.5),%或yxd=int(f,t,0,4.5),fint=double(yxd) hold on, plot(4.5,fint,'*r') %习题2-12

% clear all, syms x n; f=(sin(x))^n; yn=int(f,'x',0,pi/2), class(yn)

clear all, syms x n; syms n positive ; f=(sin(x))^n; yn=int(f,'x',0,pi/2), class(yn) % y(1/3)=?

yn1=subs(yn,'n',sym(1/3)),vpa(yn1) %或yn=limit(yn,n,1/3),vpa(yn)

%或yy=int(sin(x).^(1/3),x,0,pi/2) ,vpa(yy) %习题2-23 clear, syms x y S

S = dsolve('Dy*y/5+x/4=0','x')

ezplot(subs(S(1),'C3',1),[-2,2 -2,2],1), hold on ezplot(subs(S(2),'C3',1),[-2,2 -2,2],1)

%解为 S =

% 2^(1/2)*(C3 - (5*x^2)/8)^(1/2) % -2^(1/2)*(C3 - (5*x^2)/8)^(1/2) ezplot(subs(S(1),'C3',1),[-2,2 -2,2],1), hold on

ezplot(subs(S(2),'C3',1),[-2,2 -2,2],1) % 用此两条指令绘圆，在 y=0处有间隙 ezplot(subs(y^2-(S(1))^2, 'C3', 1),[-2,2 -2,2],2) %用椭圆方程绘图不产生间隙 colormap([0 0 1]) %用ezplot(fun)绘图时，如果fun中只有一个参数，绘图的颜色是蓝色；如果fun中有两个参数，绘图的颜色是绿色，此指令设置图形颜色为蓝。 grid on

S1=subs(S(1),'C3',1) subs(S1,{x},{1.6^(1/2)}) y=double(solve(S1)) t=linspace(y(2),y(1),100) S2=subs(S(2),'C3',1) figure

plot(t,subs(S1,x,t)), hold on plot(t,subs(S2,x,t)) axis([-2,2 -2,2]) %习题2-24

x=dsolve('Dx=a*t^2+b*t','x(0)=2','t') %习题2-25

[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0','g(0)=1','x') 第三章 %习题3-1 a=0:2*pi/9:2*pi b=linspace(0,2*pi,10) %习题3-2 rand( 'twister',0),

A=rand(3,5) [I1,J1]=find(A>0.5)

subindex_A=sub2ind(size(A),I1,J1) subindex_A=find(A>0.5)

[I,J]=ind2sub(size(A),subindex_A) index_A=[I J] 2_2

rand('twister',0),A=rand(3,5),B=A>0.5,C=find(B) [ci,cj]=ind2sub(size(A),C) %此法太繁 [Ci,Cj]=find(B) 2_5

t=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) 2_6

clear,format long, rand('twister',1),A=rand(3,3),

B=diag(diag(A)),C=A.*(~B)%或C=A-B%或C=triu(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-4 clear all A=[1 2;3 4]

B1=A.^(0.5), B2=A^(0.5), A1=B1.^2

A2=B2^2 A1-A norm(A1-A) A1-A2 norm(A1-A2) clear all

A=sym('[1 2;3 4]') B1=A.^(0.5), B2=A^(0.5), A1=simple(B1.^2) A2=simple(B2^2) A1-A vpa(A1-A) A1-A2 vpa(A1-A2) %习题3-5

t=linspace(1,10,100) %(1) L=length(t) for k=1:L

yy(k)=1-exp(-0.5*t(k))*cos(2*t(k)) end

figure(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-6 clear,format long,

rand('twister',1),A=rand(3,3), B=diag(diag(A)), C=A.*(~B)

%或C=A-B%或C=triu(A,1)+tril(A,-1) %习题3-7 clear,

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]=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-8 clear

for k=10:-1:1;

A=reshape(1:10*k,k,10); Sa(k,:)=sum(A,1); end Sa

clear

for k=10:-1:1;

A=reshape(1:10*k,k,10); if k==1 Sa(k,:)=A; else

Sa(k,:)=sum(A); end end Sa 第四章 %习题4-1 clear all load 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=20

diff_y1=(diff(y(1:N:end)))./diff(t(1:N:end));

gradient_y1=(gradient(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-2

d=0.5; tt=0:d:10; t=tt+(tt==0)*eps; y=sin(t)./t; s=d*trapz(y) % 计算出积分值 ss=d*(cumtrapz(y)) %计算梯形法累计积分并绘积 分曲线 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 % 取消警告性提示时用

%此法可用，但有警告性提示,取消提示加 warning off clear

tt=0:0.1:10; warning off for i=1:101

q(i)=quad('sin(t)./t',0,tt(i)); end plot(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-3 clear all

d=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-4 clear all

%用精度可控指令quad 即Simpson法,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/1000; x=-5*pi:d:1.7*pi; fx=exp(-abs(x)).*abs(sin(x)); s=d*trapz(fx) %习题4-5 clear all x1=-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)).^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)

下一条指令即匿名函数在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-6 tspan=[0,0.5]; y0=[1;0];

%

%

[tt,yy]=ode45(@DyDt_6,tspan,y0); y0_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-7 clear all A=magic(8) B=orth(A) rref(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 1 B =

-0.35355339059327 0.54006172486732 0.35355339059327 -0.35355339059327 -0.38575837490523 -0.35355339059327

-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.54006172486732 0.35355339059327 ans =

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-8 A = 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 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 %习题4-9 clear all, A=magic(3) b=ones(3,1) x=A\\b x =

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=ones(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 0 x=A\\b x =

0.05882352941176 0.11764705882353 -0.05882352941176 0 xg=null(A) xg =

0.22360679774998 0.67082039324994 -0.67082039324994 -0.22360679774998 x+xg 也是方程的解 ans =

0.28243032716174 0.78846745207347 -0.72964392266170

-0.22360679774998 下面逆阵法的解不正确 x=inv(A)*b x =

0.03125000000000 0.12500000000000 -0.06250000000000 -0.01562500000000 因为 A*x ans =

0.35937500000000 0.78125000000000 0.59375000000000 0.92187500000000 %习题4-11 clear 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 1 x=A\\b

x =

1.0e+015 *

-0.56294995342131 -1.68884986026394 1.68884986026394 0.56294995342131 A*x ans =

0 6.50000000000000 4.00000000000000 7.81250000000000 x=inv(A)*b x =

1.0e+015 *

-0.56294995342131 -1.68884986026394 1.68884986026394 0.56294995342131

没有准确解，以下是伪逆阵求得的最小二乘解 x=pinv(A)*b x =

0.02352941176471 0.12352941176471

0.12352941176471 0.02352941176471 A*x ans =

1.30000000000000 2.90000000000000 2.10000000000000 3.70000000000000

x=lsqnonneg(A,b) x =

0.04705882352941 0.19411764705882 0.05294117647059 0 A*x ans =

1.30000000000000 2.90000000000000 2.10000000000000 3.70000000000000 %习题4-12 clear 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 on

plot(t,zeros(size(t)),'k'); xlabel('t');ylabel('y(t)') zoom on

[tt,yy]=ginput(1),zoom off [t1,y1]=fzero(y_C,tt) t1 =

2.73411732208059 y1 =

-4.440892098500626e-016 [t2,y2]=fsolve(y_C,tt) t2 =

2.73411732208069 y2 =

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)

hold 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 = % % 1

或者采用下面的程序

x0 = [-1/4*pi; -1/4*pi]; % Make a starting guess at the solution options=optimset('Display','iter'); % Option to display output [x,fval] = fsolve(@Myfsove,x0,options) % Call solver x =

-5.49778714378216 -5.49778714378216 fval =

1.0e-013 *

0 0.43980294605305

x0 = [-pi/10; -pi/2]; % Make a starting guess at the solution options=optimset('Display','iter'); % Option to display output [x,fval] = fsolve(@Myfsove,x0,options) % Call solver 函数Myfsove另存为一个文件 function F=Myfsove(x) F = [sin(x(1)-x(2)); cos(x(1)+x(2))]; %习题4-14

y = binopdf([0:100],100,0.157); stem(1:length(y),y), axis([0 length(y) 0 .12 ]) %习题4-15

a=normrnd(4,2,10000,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-17

format 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='商多项式为 ';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-18

load prob_418, who,x

P=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); end

X(:,n+1)=ones(m,1); aT=(X\\y0)', norm(P-aT) %习题4-18

h=[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([0 length(y) -1 1 ]) subplot(2,1,2),stem(y,'filled') axis([0 length(y) -1 1 ]) 第五章

第五章 1,2,3,4,5,6,7 %习题5_1

t=2*pi*(0:199)/199; a=4;b=2;

x=a*cos(t);y=b*sin(t); plot(x,y,'r.','MarkerSize',15) axis equal

%习题5_2

t = 0:.01:2*pi; P=1-cos(t); pline=polar(t,P,'r'), set(pline,'LineWidth',5) title('P=1-cos\\theta')

%习题5_3 clear x=1:6; Y=[170 120 180 200 190

120 100 110 180 170 180; 70 50

80

100

95

120]

bar(x,Y','style'); % bar(x,Y','grouped'); % bar(x,Y','stacked'); colormap(cool);

% legend('A','B','C','Location','NorthWest') legend('A','B','C',2)

%习题5_4 % exmp504.m

clc,clf,clear; t=(0:0.05:18)'; N=length(t); zeta=0.2:0.2:1.4; L=length(zeta); y=zeros(N,L); hold on

220;

供第4道习题使用的程序

for k=1:L zk=zeta(k);

beta=sqrt(abs(1-zk^2)); if zk<1

缺陷在此，由于计算机的精度，zeta(5)<1 ,可

改为zk-1<-2*eps

%满足此条件，绘蓝色线

y=1/beta*exp(-zk*t).*sin(beta*t); plot(t,y,'b') if zk<0.4

text(2.2,0.63,'\\zeta = 0.2') end elseif

zk==1

缺陷在此，由于计算机的精度，zk！=1

改为(zk-1)<2*eps%满足此条件，绘黑色线 y=t.*exp(-t);

plot(t,y,'k','LineWidth',2) else

%其余，绘红色线

y=(exp(-(zk-beta)*t)-exp(-(zk+beta)*t))/(2*beta); plot(t,y,'r') if zk>1.2

text(0.3,0.14,'\\zeta = 1.4') end end end

text(10,0.7,'\\Delta\\zeta=0.2') axis([0,18,-0.4,0.8]) hold off box on grid on

,可

%习题5_5

t=4*pi*(0:100)/100; x=sin(t);y=cos(t);z=t;

plot3(x,y,z,'g','LineWidth',3),box on

%习题5_6

x=-3:0.1:3;y=x;[X,Y]=meshgrid(x,y); Z=4.*X.*exp(-X.^2-Y.^2); mesh(X,Y,Z) hidden off % colormap(cool), % shading interp,

syms x y z

% z=4.*x.*exp(-x.^2-y.^2); z=4*x*exp(-x^2-y^2); ezmesh(z,[-3,3]) hidden off

%习题5_7 clear all

x=4*pi*(-50:50)/50;y=x;[X,Y]=meshgrid(x,y); Z=sin(X+Y)./(X+Y+(X+Y==0)*eps); surf(X,Y,Z)

view([21,32]) %图形界面旋转图形，认为合适后记下方位角和俯视角，再写出命令 shading interp

% size(find(isnan(Z))) % sum(sum(isnan(Z)))

%习题5_8

ezplot('y/(1+x^2+y^2)-0.1',[-2*pi,2*pi,-pi/3,3.5*pi]) hold on

ezplot('sin(x+cos(y))',[-2*pi,2*pi,-pi/3,3.5*pi])

可看到6个交点，及方程组有6个实数解 hold off grid on

[xx,yy]=ginput1)

要解最接近x=0,y=0的解，首先将‘myfun8’另存为一个文件 function F=myfun8(x,y)

F=[y/(1+x^2+y^2)-0.1;sin(x+cos(y))] end zoom on xy= ginput(1) f=fsolve(@myfun8,xy)

[X,Y] = ginput(1) f=fsolve(@myfun8,[X,Y]) clear all syms x y

s=solve('y/(1+x^2+y^2)-0.1','sin(x+cos(y))') 在2010版求解即得最接近x=0,y=0的解， %习题5_9

在6.5版本无法运行%习题5_9

function f=prob5_9 clf;

[X,Y,Z]=sphere(40);

surf(X,Y,Z),axis equal off,shading flat set(gcf,'Color','w') material shiny

light('position',[-1,0.5,1]) light('position',[2,1,1]) colormap([jet; flipud(jet)])

disp('按任意键，观察色图变幻。退出按Ctrl+C') pause spinmap(80,9) %习题5_10p

function f=prob5_10(K,ki)

%prob5_10函数产生动态衰减正弦函数，K控制动态曲线动态变化的循环次数，ki控制曲线动态变化的快慢 nargin=0 clf t=0;k=0;

x=(0:100)/100*4*pi; n=length(x); y0=zeros(n);

plot(x,y0,'b','LineWidth',2.5);hold on plot(x(1),y0(1),'d','MarkerFace',[0,0,1]); axis off; if nargin==0 K=1;ki=1; elseif nargin==1 ki=1;

end while 1 t=t+1;

y1=sin(pi*t/24); y(2:n)=y0(1:n-1); y(1)=y1; y0=y;

y=exp(-0.2.*x).*y;

plot(x,y,'LineWidth',2.5);hold on

plot(x(1),y(1),'p','MarkerSize',20,'MarkerFace',[0,0,1]);hold off axis([ 0 14 -1.2 1.2]) axis off;

if round(t/240)==K break end

pause(0.001*ki) end

f=getframe(gcf)

%习题5_10p,用line函数 function f=prob5_10(K,ki)

%prob5_10函数产生动态衰减正弦函数，K控制动态曲线动态变化的循环次数，ki控制曲线动态变化的快慢 nargin=0 clf t=0;k=0;

x=(0:100)/100*4*pi; n=length(x); y0=zeros(n);

h1=line(x(1),y0(1),'Marker', 'd','MarkerSize',20,'MarkerFace',[0,0,1]); h=line(x,y0,'LineWidth',2.5,'Color',[0 0 1]);

axis off;

% set(h,'erasemode','background') %h是需要执行动画图像的句柄，一般都是由line或者plot创建 if nargin==0 K=1;ki=1; elseif nargin==1 ki=1; end while 1 t=t+1;

y1=sin(pi*t/24); y(2:n)=y0(1:n-1); y(1)=y1; y0=y;

y=exp(-0.2.*x).*y; axis([ 0 14 -1.2 1.2])

set(h1,'xdata',x(1),'ydata',y(1),'erasemode','xor')%更新图像的坐标数据 drawnow %刷新屏幕

set(h,'xdata',x,'ydata',y) %更新图像的坐标数据 drawnow %刷新屏幕 if round(t/240)==K break end

pause(0.001*ki) end

f=getframe(gcf) %习题5_11

function f=prob5_11(K,ki)

%prob5_10函数产生驻波动画，K控制驻波动画动态变化的次数，ki控制曲线动态变化的快慢 nargin=0 clf t=0;k=0;

x=(0:100)/100*3*pi; n=length(x); x0=[0 pi 2*pi 3*pi]; y0=zeros(1,4); y=sin(t).*sin(x);

h=line(x,y,'LineWidth',4,'Color',[0 0 1]);

h1=line(x0,y0,'LineWidth',4,'Marker', '.','MarkerSize',50,'Color',[1 0 0]); axis off;

% set(h,'erasemode','background') %h是需要执行动画图像的句柄，一般都是由line或者plot创建 if nargin==0 K=1;ki=1; elseif nargin==1 ki=1; end while 1 t=t+0.5; y=sin(t).*sin(x); axis([ 0 10 -1.1 1.1])

set(h,'xdata',x,'ydata',y) %更新图像的坐标数据 drawnow %刷新屏幕

set(h1,'xdata',x0,'ydata',y0,'erasemode','background')%更新图像的坐标数据 drawnow %刷新屏幕 if round(t/48)==K

break end

pause(0.1*ki) end

f=getframe(gcf) %习题5_12

function f=prob5_12(K,ki) %prob5_12 % K

演示红色小球沿一条封闭三叶线运动的实时动画 红球运动的循环数（不小于1）

指定拍摄照片的瞬间，取 1 到 500 间的任意整数。 存储拍摄的照片数据，可用image(f.cdata)观察照片。

% ki % f

% nargin=0 if nargin==0 K=1;ki=1; elseif nargin==1 ki=1; end

theta1=(0:1000)/1000*2*pi; rho=cos(3*theta1); % h=polar(theta1,rho,'r')

% set(h,'Color',[0.5,0.5,0.5],'LineWidth',2.5); y=-rho.*cos(theta1);x=rho.*sin(theta1); % h=line(x,y,'LineWidth',4,'Color',[0.5,0.5,0.5]);

% h1=line(x(251),y(251),'LineWidth',4,'Marker', '.','MarkerSize',50,'Color',[1 0 0]); yt=y;

y(1:751)=yt(251:end); y(752:end)=yt(1:250); xt=x;

x(1:751)=xt(251:end);

x(752:end)=xt(1:250);

h1=line(x(1),y(1),'LineWidth',4,'Marker', '.','MarkerSize',50,'Color',[1 0 0]); h=line(x,y,'LineWidth',4,'Color',[0.5,0.5,0.5]);

axis off

n=length(x);i=1;j=1;

while 1

set(h1,'xdata',x(i),'ydata',y(i),'erasemode','xor'); drawnow; <22>

pause(0.0005)

% <23> i=i+1;

if nargin==2 & nargout==1

if(i==ki&j==1);f=getframe(gcf);end

end

if i>round(n/2) i=1;j=j+1; if j>K;break;end end end 第六章

第六章1,2,3 %习题6_1 %for循环 S=0.2;SUM=1; tic for i=0:1000000

SUM=SUM+S; S=S*0.2; end t=toc, SUM

% <26>

%

%while循环

S=0.2; SUM=1; i=0; tic while i<=1000000

SUM=SUM+S; S=S*0.2; i=i+1; end t=toc, SUM %sym

syms i; tic, SUM=vpa(symsum(0.2.^i,0,1000000),6), t=toc %向量法(用时最短)

N=1000000; S=zeros(1,N);tic k=0:N; S=0.2.^k; sum(S),t=toc %函数法

function y=serisum(x,n) y=(1-x.^n)./(1-x);

调用：serisum(0.2,1000001)

%习题6_2 function prob6_2(n)

%prob6_2 function:draw unit circle or regular polygon of n sides. % n The number of sides, while n>=3; % draw circle, while no parameter;

% not a appropriate \

% edited by zhongguo liu, 12 Dec. 2010

if nargin~=0 & (n<0 | fix(n)~=n)

disp('Input parameter error! Please input natural number n larger than 2.')

else

switch nargin case 0

m=100;t=(0:m)/m*2*pi; x=sin(t);y=cos(t);

plot(x,y,'r','LineWidth',2.5);title('Circle') % ss=exp(i*t);

% plot(real(ss),imag(ss)) axis equal, axis off otherwise

if n==1 | n==2,disp('please input n>=3') else

t=(0:n)/n*2*pi; x=sin(t);y=cos(t);

plot(x,y,'r','LineWidth',2.5); title(['Poly gon nit',int2str(n),' edges'])

% ss='exp(i*t)'; % se=eval(ss);

% plot(real(se),imag(se)) axis equal, axis off end end end %调用

prob6_2(59)，prob6_2(-2)

%%习题6_2另解 function prob6_2_(N)

if nargin~=0 &(N<3 | fix(N)~=N)

disp('input parameter error!Please input natural number n larger than 2.') elseif nargin==0

t=-pi:pi/100:pi;plot(sin(t),cos(t)), axis equal square title('Circle') else

t=-pi:2*pi/(N):pi;plot(sin(t),cos(t)), axis equal square title(['Poly gon nit ',int2str(n),' edges']); end

%习题6_2（学生）解 function aa(n) if nargin==0

t=0:pi/100:2*pi; x=exp(i*t); str='Circle'; else

if(nargin~=0)&(n<=2) error('输入边数少了') end

if n-round(n)~=0 error('输入了非自然数') end

t=(0:n)/n*2*pi; x=exp(i*t);

str=['Poly gon nit ',int2str(n),' edges'] end

plot(real(x),imag(x),'r','LineWidth',4) title(str)

axis square image off shg

%习题6_3

匿名函数表达式法在6.5版本不能用 %习题6_3匿名函数表达式法，@(x)(exp)

y0=@(x)(-exp(-x)*abs(sin(cos(x)))); fplot(y0,[-6,4]) 或 ezplot(y0,[-6,4]) [x,fval]=fminbnd(y0,-1,1); hold on, plot(x,fval,'*r')

%习题6_3子程序函数句柄法_function function y=prob6_3(x) y=-exp(-x)*abs(sin(cos(x))); end

%main program

hh=@prob6_3; ezplot(hh,[-6,4]) , [x,fval]=fminbnd(hh,-1,1); hold on, plot(x,fval,'*r')

%习题6_3_inline内联函数法

y1=inline('-exp(-x)*abs(sin(cos(x)))'); [x,fval]=fminbnd(y1,-1,1) x0=[-6,4] ; ezplot(y1,x0), hold on, plot(x,fval,'*r')

char字符串法

y2='-exp(-x)*abs(sin(cos(x)))'; [x,fval]=fminbnd(y2,-1,1) x0=[-6,4] ; ezplot(y2,x0), hold on, plot(x,fval,'*r')

%习题6_4_ clear b=pwd

%获取当前目录名字符串

which('smoke') b_d=b;

%检查在当前目录下能否看到smoke.m

b_d(end-4:end)=[]; %在b字符串中去除最后的四个字符，即\\work。

str=[b_d,'\\toolbox\\matlab\\elmat\\private']; cd(str)

%把smoke.m 所在目录设置为当前目录。 %创建smoke.m的函数句柄

%使work恢复为当前目录

SM=@smoke; cd(b) disp(' ')

% A=SM(3,0,'double') 此指令适用于2010版本，，只能用三个参数，%

产生一个3阶“伪特征根”矩阵

A=feval(SM,3,0) 此指令适用于6.5版本，只能用两个参数 disp(' ') pwd

%显示当前所在目录，以证实符

合题意。 which('smoke') smoke.m。disp(' ')

%再次检查在当前目录下能否看到

本文来源：https://www.bwwdw.com/article/my2r.html