matlab教程(张志涌)课后习题答案

第二章

%分段函数的写法
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]
%122_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)
%122_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)
%122_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,x

0,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 220;
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 供第4道习题使用的程序
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
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）
% ki 指定拍摄照片的瞬间，取 1 到 500 间的任意整数。
% f 存储拍摄的照片数据，可用image(f.cdata)观察照片。
% 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 % <26>
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
%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 functio

n: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 "n" ,while n<3 or not natural number.

% 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') %检查在当前目录下能否看到smoke.m
b_d=b;
b_d(end-4:end)=[]; %在b字符串中去除最后的四个字符，即\work。
str=[b_d,'\toolbox\matlab\elmat\private'];
cd(str) %把smoke.m 所在目录设置为当前目录。
SM=@smoke; %创建smoke.m的函

数句柄
cd(b) %使work恢复为当前目录
disp(' ')
% A=SM(3,0,'double') 此指令适用于2010版本，，只能用三个参数，%产生一个3阶“伪特征根”矩阵
A=feval(SM,3,0) 此指令适用于6.5版本，只能用两个参数
disp(' ')
pwd %显示当前所在目录，以证实符合题意。
which('smoke') %再次检查在当前目录下能否看到smoke.m。disp(' ')