MATLAB程序设计与应用（第二版）
实验参考答案
%实验一 MATLAB运算基%第四题 %(1):
A=100:999; B=rem(A,21);
础
%第一题 %（1）
z1=2*sin(85*pi/180)/(1+exp(2)) %（2）
x=[2,1+2i;-0.45,5];
z2=0.5*log(x+sqrt(1+x.^2)) %(3)
a=-3.0:0.1:3.0;
z3=(exp(0.3*a)-exp(-0.3*a))/2.*sin(a+0.3)+log((0.3+a)/2) %(4)
t=0:0.5:2.5;
z4=t.^2.*(t>=0&t<1)+(t.^2-1).*(t>=1&t<2)+(t.^2-2*t+1).*(t>=2&t<3)
%第二题
A=[12 34 -4;34 7 87;3 65 7]; B=[1 3 -1;2 0 3;3 -2 7]; A+6*B
A-B+eye(size(A)) A*B A.*B A^3 A.^3 A/B B\\A [A,B]
[A([1,3],:);B^2]
%第三题
A=[1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20;21 22 23 24 25]
B=[3 0 16;17 -6 9;0 23 -4;9 7 0;4 13 11] C=A*B
F=size(C)
D=C(F(1)-2:F(1),F(2)-1:F(2)) whos
C=length(find(B==0)) %(2):
A='lsdhKSDLKklsdkl'; k=find(A>='A'&A<='Z'); A(k)=[]
%实验二 MATLAB矩阵分析与处理
%第一题
E=eye(3); R=rand(3,2); O=zeros(2,3); S=diag([2,3]); A=[E,R;O,S]; A^2
B=[E,(R+R*S);O,S^2]
%第二题
H=hilb(5) P=pascal(5) Hh=det(H) Hp=det(P) Th=cond(H) Tp=cond(P)
%第三题:
A=fix(10*rand(5)) H=det(A) Trace=trace(A) Rank=rank(A) Norm=norm(A)
%第四题:
A=[-29,6,18;20,5,12;-8,8,5] [V,D]=eig(A) %数学意义略
%第五题方法一:
%(1):
A=[1/2,1/3,1/4;1/3,1/4,1/5;1/4,1/5,1/6];
b=[0.95,0.67,0.52]'; x=inv(A)*b
1
%(2):
B=[0.95,0.67,0.53]'; x=inv(A)*B %(3): cond(A)
0];
y=[]; %建立存放所有y值的矩阵 for x0=x
if x0<0&x0~=-3
y=[y,x0*x0+x0-6]; elseif
x0>=0&x0<5&x0~=2&x0~=3
y=[y,x0*x0-5*x0+6]; else
y=[y,x0*x0-x0-1]; end end
x %输出所有x
y %输出所有y
%第一题程序二
x=[-5,-3,1,2,2.5,3,5]; y=[];
for a=1:7
if x(a)<0&x(a)~=-3
y=[y,(x(a))^2+x(a)-6]; elseif
x(a)>=0&x(a)<5&x(a)~=2&x(a)~=3 y=[y,(x(a))^2-5*x(a)+6]; else
y=[y,x(a)*x(a)-x(a)-1]; end end
%第二题程序一:
score=input('请输入一个分数:'); if score<0&score>100 disp('输入的分数不合理'); else
if score>=90&score<=100 disp('等级为:A');
elseif score>=80&score<=89 disp('等级为:B');
elseif score>=70&score<=79 disp('等级为:C');
elseif score>=60&score<=69 disp('等级为:D'); else
%第五题方法二:
A=hilb(4) A(:,1)=[] A(4,:)=[]
B=[0.95,0.67,0.52]'; X=inv(A)*B
B1=[0.95,0.67,0.53]'; X1=inv(A)*B1 N=cond(B) N1=cond(B1)
Na=cond(A) %矩阵A为病态矩阵
%第六题
A=[1,4,9;16,25,36;49,64,81] B=sqrtm(A)
C=sqrt(A) %sqrtm函数是以矩阵为单位进行计算，sqrt函数是以矩阵中的元素进行计算
%实验三 选择程序结构设计
%第一题程序一
x=[-5.0,-3.0,1.0,2.0,2.5,3.0,5.0]; y=[]; for x0=x
if x0<0&x0~=-3
y=[y,x0*x0+x0-6];
elseif x0>=0&x0<5&x0~=2&x0~=3 y=[y,x0*x0-5*x0+6]; else
y=[y,x0*x0-x0-1]; end end x
y
x=[-5.0,-3.0,1.0,2.0,2.5,3.0,5.
2
disp('等级为:E'); end end
x=input('请输入一个百分制成绩：'); if x>100|x<0
disp('您输入的成绩不是百分制成绩，s=input('请输入一个成绩（0分到100分之间）：'); %s用于存放成绩 while
1 %判断输入成绩的合理性 if s<0|s>100
disp('输入的成绩需在0到100之间，请重新输入：')
s=input('请输入一个成绩（0请重新输入。'); else
if x<=100&x>=90 disp('A');
elseif x<=89&x>=80 disp('B');
elseif x<=79&x>=70 disp('C');
elseif x<=69&x>60 disp('D'); else
disp('E'); end end
%第二题程序二：
score=input('请输入一个分数:'); if score<0&score>100 disp('输入的分数不合理'); else
switch fix(score/10) case {9,10} disp('A'); case {8}
disp('B'); case {7}
disp('C'); case {6}
disp('D'); otherwise
disp('E'); end end
分到100分之间）：'); else
break; end end switch
fix(s/10) %对成绩做出等级判断 case {9,10} disp('A') case 8
disp('B') case 7
disp('C') case 6
disp('D') otherwise disp('E') end
%第三题
a=input('工号:'); b=input('工时:'); if b>120
c=120*84+(b-120)*84*(1+0.15); elseif b<60 c=84*b-700; else
c=84*b; end
disp(['工号为',num2str(a),'的员工的工资为',num2str(c)]);
3
n=input('请输入员工工号：'); h=input('该员工工作时数是：'); if h>120
x=(h-120)*84*(1+0.15)+120*84; elseif h<60 x=h*84-700; else
x=h*84; end
disp([num2str(n),'号员工','的应发工资为',num2str(x)]);
%第四题（还可以用switch语句实现） a=input('a='); b=input('b=');
c=input('输入一个四则运算符号：','s'); if c=='+' d=a+b;
elseif c=='-' d=a-b; elseif c=='*' d=a*b;
elseif c=='/' d=a/b; end
disp(['结果为:',num2str(d)]);
a=fix(10+(99-10)*rand(1,2)) %产生两个随机整数 x=a(1); y=a(2);
t=input('请输入运算符号：','s'); if t=='+' z=x+y;
elseif t=='-' z=x-y;
elseif t=='*' z=x*y;
elseif t=='/' z=x/y; end
disp([num2str(x),t,num2str(y),'=',num2str(z)]) %输出运算结果 %第五题
a=rand(5,6);
n=input('输出矩阵的第几行:'); if n>5
disp(['超出了矩阵的行数,输出矩阵的最后一行为:',num2str(a(5,:))]); else
disp(['矩阵的第',num2str(n),'行为:',num2str(a(n,:))]); end
a=rand(5,6) %产生5x6的随机矩阵
n=input('请输入您要输出矩阵的第几行：'); if n>5
disp('超出了矩阵的行数，矩阵的最后一行为：') a(5,:) else
disp(['矩阵的第',num2str(n),'行为：']) a(n,:) end
%实验四 循环结构程序设计
%第一题程序一 s=0;
n=input('n=?'); for i=1:n
s=s+1/i/i; end
PI=sqrt(6*s) pi
%第一题程序二
n=input('n=?'); a=1:n;
4
b=1./a.^2;
PI=sqrt(6*sum(b)) pi
%第二题 y=0; n=1;
while(y<3)
y=y+1/(2*n-1); n=n+1; end
y=y-1/(2*(n-1)-1) n=n-2
%第三题
a=input('a=?'); b=input('b=?'); Xn=1;
Xn1=a/(b+Xn); n=0;
while abs(Xn1-Xn)>1e-5 Xn=Xn1;
Xn1=a/(b+Xn); n=n+1; if n==500 break; end end n Xn1
r1=(-b+sqrt(b*b+4*a))/2 r2=(-b-sqrt(b*b+4*a))/2 %第四题
for i=1:100 if i==1 f(i)=1; elseif i==2 f(i)=0; elseif i==3 f(i)=1; else
f(i)=f(i-1)-2*f(i-2)+f(i-3);
end end max(f) min(f) sum(f)
length(find(f>0)) length(find(f==0)) length(find(f<0)) %第五题： s=0;n=0; for i=2:49
b=i*(i+1)-1; m=fix(sqrt(b)); for j=2:m
if rem(b,j)==0 break end end
if j==m n=n+1; s=s+b; end end n s
%实验五 函数文件
%第一题
function y=mat1(x) %建立函数文件mat1.m
y=[exp(x),log(x),sin(x),cos(x)];
%在命令窗口调用上述函数文件： y=mat1(1+i)
%第二题程序一 function
[a,b,N,M]=shiyanwu2(m,n,t)
A=[m*cos(t*pi/180),-m,-sin(t*pi/180),0;m*sin(t*pi/180),0,cos(t*pi/180),0;0,n,-sin(t*pi/180),0;0,0,-cos(t*pi/180),1]; B=[0,9.8*m,0,9.8*n]; C=inv(A)*B'; a=C(1);
5
b=C(2); N=C(3); M=C(4);
%在命令窗口调用该函数文件: m1=input('m1='); m2=input('m2=');
theta=input('theta=');
[a1,a2,N1,N2]=shiyanwu2(m1,m2,t%第五题 %(1)
function f1=mat5(n) f1=n+10*log(n*n+5);
%在命令窗口中调用该函数文件：
y=mat5(40)/(mat5(30)+mat5(20)) %(2)方法一
function f2=mat6(n) heta)
%第二题程序二
function X=mat2(m1,m2,t) g=9.8;
A=[m1*cos(t*pi/180),-m1,-sin(t*pi/180),0;m1*sin(t*pi/180),0,cos(t*pi/180),0;0,m2,-sin(t*pi/180),0;0,0,-cos(t*pi/180),1]; B=[0;m1*g;0;m2*g]; X=inv(A)*B;
%在命令窗口调用该函数文件： X=mat2(1,1,60)
%第三题
function flag=mat3(x) flag=1;
for i=2:sqrt(x) if rem(x,i)==0 flag=0; break; end end
%在命令窗口调用该函数文件： for i=10:99
j=10*rem(i,10)+fix(i/10); if mat3(i)&mat3(j) disp(i) end end
%第四题
function y=fx(x)
y=1./((x-2).^2+0.1)+1./((x-3).^4+0.01);
%在命令窗口调用该函数文件： y=fx(2)
a=[1,2;3,4]; y=fx(a)
f2=0;
for i=1:n
f2=f2+i*(i+1); end
%在命令窗口中调用该函数文件如：
y=mat6(40)/(mat6(30)+mat6(20)) %(2)方法二
function f2=mat7(n) i=1:n;
m=i.*(i+1); f2=sum(m); end
%在命令窗口中调用该函数文件如：
y=mat7(40)/(mat7(30)+mat7(20))
%实验六 高层绘图操作
%第一题：
x=linspace(0,2*pi,101);
y=(0.5+3*sin(x)./(1+x.^2)).*cos(x);
plot(x,y)
%第二题： %（1）
x=linspace(-2*pi,2*pi,100); y1=x.^2;
y2=cos(2*x); y3=y1.*y2;
plot(x,y1,'b-',x,y2,'r:',x,y3,'y--');
text(4,16,'\\leftarrow y1=x^2'); text(6*pi/4,-1,'\\downarrow y2=cos(2*x)');
text(-1.5*pi,-2.25*pi*pi,'\%uparrow y3=y1*y2');
%（2）
x=linspace(-2*pi,2*pi,100);
6
y1=x.^2;
y2=cos(2*x); y3=y1.*y2;
subplot(1,3,1);%分区 plot(x,y1);
title('y1=x^2');%设置标题 subplot(1,3,2); plot(x,y2);
title('y2=cos(2*x)'); subplot(1,3,3); plot(x,y3);
title('y3=x^2*cos(2*x)'); %（3）
x=linspace(-2*pi,2*pi,20); y1=x.^2;
subplot(2,2,1);%分区 bar(x,y1);
title('y1=x^2的条形图');%设置标题 subplot(2,2,2); stairs(x,y1);
title('y1=x^2的阶梯图'); subplot(2,2,3); stem(x,y1);
title('y1=x^2的杆图'); subplot(2,2,4);
fill(x,y1,'r');%如果少了'r'则会出错
title('y1=x^2的填充图'); %其他的函数照样做。 %第三题
x=-5:0.01:5;
y=[];%起始设y为空向量 for x0=x
if x0<=0 %不能写成x0=<0
y=[y,(x0+sqrt(pi))/exp(2)]; %将x对应的函数值放到y中 else
y=[y,0.5*log(x0+sqrt(1+x0^2))]; end end
plot(x,y) %第四题：
a=input('a=');
b=input('b='); n=input('n=');
t=-2*pi:0.01:2*pi; r=a*sin(b+n*t); polar(t,r)
%第五题
x=linspace(-5,5,21); y=linspace(0,10,31);
[x,y]=meshgrid(x,y);%在[-5,5]*[0,10]的范围内生成网格坐标 z=cos(x).*cos(y).*exp(-sqrt(x.^2+y.^2)/4);
subplot(2,1,1); surf(x,y,z); subplot(2,1,2);
contour3(x,y,z,50);%其中50为高度的等级数，越大越密 %第六题
ezsurf('cos(s)*cos(t)','cos(s)*sin(t)','sin(s)',[0,0.5*pi,0,1.5*pi]); %利用ezsurf隐函数
shading interp %进行插值着色处理
%实验七 低层绘图操作
%第一题
h=figure('MenuBar','figure','color','r','WindowButtonDownFcn','disp(''Left Button Pressed'')') %第二题
x=-2:0.01:2;
y=x.^2.*exp(2*x); h=line(x,y);
set(h,'color','r','linestyle',':','linewidth',2)
text(1,exp(2),'y=x^2*exp(2*x)') %第三题
t=0:0.00001:0.001; [t,x]=meshgrid(t);
v=10*exp(-0.01*x).*sin(2000*pi*t-0.2*x+pi);
axes('view',[-37.5,30]); h=surface(t,x,v);
title('v=10*exp(-0.01*x).*sin(2000*pi*t-0.2*x+pi)');
7
xlabel(Ct'),ylabel('x'),zlabel('v')
%第四题
x=0:0.01:2*pi; y1=sin(x); y2=cos(x); y3=tan(x); y4=cot(x);
subplot(2,2,1); plot(x,y1);
subplot(2,2,2); plot(x,y2);
subplot(2,2,3); plot(x,y3);
subplot(2,2,4); plot(x,y4);
%第五题
cylinder(5);
light('Position',[0,1,1]); material shiny
%实验八 数据处理与多项式运算
%第一题 %(1)
A=rand(1,30000); b=mean(A) std(A,0,2) %(2) max(A) min(A) %(3) n=0;
for i=1:30000 if A(i)>0.5 n=n+1; end end
p=n/30000
%第二题 %(1)
A=45+51*rand(100,5); [Y,U]=max(A)
[a,b]=min(A) %(2)
m=mean(A) s=std(A) %(3)
sum(A,2)
[Y,U]=max(ans) [a,b]=min(ans) %(4)
[zcj,xsxh]=sort(ans) %第三题 h=6:2:18;
x=6.5:2:17.5;
t1=[18,20,22,25,30,28,24]; t2=[15,19,24,28,34,32,30]; T1=spline(h,t1,x) T2=spline(h,t2,x) %第四题
x=1:0.1:101; y1=log10(x);
p=polyfit(x,y1,5) y2=polyval(p,x);
plot(x,y1,':',x,y2,'-') %第五题 %(1)
p1=[1,2,4,0,5]; p2=[1,2]; p3=[1,2,3];
p=p1+[0,conv(p2,p3)] %为使两向量大小相同，所以补0 %(2)
A=roots(p) %(3)
A=[-1,1.2,-1.4;0.75,2,3.5;0,5,2.5];
polyval(p,A) %(4)
polyvalm(p,A)
实验十
程序：
x=sym('6'); y=sym('5');
z=(x+1)/(sqrt(3+x)-sqrt(y))
8
1、 分解因式 （1） 程序： syms x y; A=x^4-y^4; factor(A) （2） 程序：
factor(sym('5135')) 3、化简表达式 （1） 程序：
syms beta1 beta2
y=sin(beta1)*cos(beta2)-cos(beta1)*sin(beta2) simple(y) （2） 程序：
syms x
y=(4*x^2+8*x+3)/(2*x+1) simple(y)
5、用符号方法求下列极限或导数 （1） 程序：
syms x
f=(x*(exp(sin(x))+1)-2*(exp(tan(x))-1))/(sin(x)) limit(f) （2） 程序：
syms x
y=(sqrt(pi)-sqrt(acos(x)))/(sqrt(x+1));
limit(f,x,-1,'right') （3） 程序：
syms x
y=(1-cos(2*x))/x; y1=diff(y) y2=diff(y,x,2)
6、用符号方法求下列积分 （1） 程序： syms x
f=1/(1+x^4+x^8) int(f) （2） 程序：
syms x
f=1/(((asin(x))^2)*sqrt(1-x^2)) int(f) （3） 程序：
syms x
f=(x^2+1)/(x^4+1) int(f,x,0,inf) （4） 程序：
syms x
f=exp(x)*(1+exp(x))^2 y=int(f,x,0,log(2)) double(y)
实验十一 级数与方程符号求解 1. 级数符号求和。 (1) 计算 。
(2) 求级数 的和函数，并求 之和。 解：
M文件： clear all;clc;
n=sym('n');x=sym('x');
S1=symsum(1/(2*n-1),n,1,10) S2=symsum(n^2*x^(n-1),n,1,inf)
S3=symsum(n^2/5^n,n,1,inf) %vpa(S3)可以转化成小数 运行结果： S1 =
31037876/14549535 S2 =
piecewise([abs(x) < 1, -(x^2 + x)/(x*(x - 1)^3)]) S3 = 15/32
2. 将lnx在x=1处按5次多项式展开为泰勒级数。 解：
M文件： clear all;clc; x=sym('x');
9
f=log(x); clear all;clc; taylor(f,x,6,1) dsolve('D2y+4*Dy+29*y','y(0)=0','Dy(0)=15','
x') 运行结果：
ans = x - (x - 1)^2/2 + (x - 1)^3/3 - (x - 1)^4/4 + (x - 运行结果： 1)^5/5 - 1 ans =
(3*sin(5*x))/exp(2*x) 3. 求下列方程的符号解。
解： 5. 求微分方程组的通解。 M文件： 解： clear all;clc; M文件： x1=solve('log(x+1)-5/(1+sin(x))=2') clear all;clc; x2=solve('x^2+9*sqrt(x+1)-1') [x y z]=dsolve('Dx=2*x-3*y+3*z',... x3=solve('3*x*exp(x)+5*sin(x)-78.5') 'Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t') [x4 运行结果： y4]=solve('sqrt(x^2+y^2)-100','3*x+5*y-8') x =
C1/exp(t) + C2*exp(2*t) 运行结果：
x1 = y = 521.67926389905839979437366649258 C1/exp(t) + C2*exp(2*t) + C3/exp(2*t) x2 = z =
C2*exp(2*t) + C3/exp(2*t) -1 实验九 数值微积分与方程数值求解 (3^(1/2)*i*(4/(9*(6465^(1/2)/2 + 1. 求函数在指定点的数值导数。 2171/54)^(1/3)) - (1/2*6465^(1/2) + 解：M文件： 2171/54)^(1/3)))/2 - (6465^(1/2)/2 + clc;clear; 2171/54)^(1/3)/2 - 2/(9*(6465^(1/2)/2 + x=1; 2171/54)^(1/3)) + 1/3 i=1; 1/3 - (6465^(1/2)/2 + 2171/54)^(1/3)/2 - f=inline('det([x x^2 x^3;1 2*x 3*x^2;0 2 (3^(1/2)*i*(4/(9*(6465^(1/2)/2 + 6*x])'); 2171/54)^(1/3)) - (1/2*6465^(1/2) + while x<=3.01 2171/54)^(1/3)))/2 - 2/(9*(6465^(1/2)/2 + g(i)=f(x); 2171/54)^(1/3)) i=i+1; x3 = x=x+0.01; %以0.01的步2.3599419584772910151699327715486 长增加，可再缩小步长提高精度 x4 = end 12/17 - (10*21246^(1/2))/17 g; (10*21246^(1/2))/17 + 12/17 t=1:0.01:3.01; y4 = dx=diff(g)/0.01; %差分法近似求导 (6*21246^(1/2))/17 + 20/17 f1=dx(1) %x=1的数值倒数 20/17 - (6*21246^(1/2))/17 f2=dx(101) %x=2的数值倒数 4. 求微分方程初值问题的符号解，并与数值f3=dx(length(g)-1) %x=3的数值倒数 解进行比较。 运行结果： f1 =
6.0602 解：
f2 = M文件：
10
24.1202 f3 =
54.1802
2. 用数值方法求定积分。 (1) 的近似值。 (2)
解：M文件： clc;clear;
f=inline('sqrt(cos(t.^2)+4*sin(2*t).^2+1)'); I1=quad(f,0,2*pi)
g=inline('log(1+x)./(1+x.^2)'); I2=quad(g,0,2*pi) 运行结果：
3. 分别用3种不同的数值方法解线性方程组。
解：M文件： clc;clear;
A=[6 5 -2 5;9 -1 4 -1;3 4 2 -2;3 -9 0 2]; b=[-4 13 1 11]'; x=A\\b
y=inv(A)*b [L,U]=lu(A); z=U\\(L\\b) 运行结果：
4. 求非齐次线性方程组的通解。 解：M文件
function [x,y]=line_solution(A,b) [m,n]=size(A); y=[ ];
if norm(b)>0 %非齐次方程组 if rank(A)==rank([A,b]) if rank(A)==n
disp('有唯一解x'); x=A\\b; else
disp('有无穷个解，特解x，基础解系y');
x=A\\b;
y=null(A,'r'); end else
disp('无解'); x=[ ];
11
end
else %齐次方程组 disp('有零解x'); x=zeros(n,1); if rank(A) disp('有无穷个解，基础解系y'); y=null(A,'r'); end end clc;clear; format rat A=[2 7 3 1;3 5 2 2;9 4 1 7]; b=[6 4 2]'; [x,y]=line_solution(A,b) 运行结果： 有无穷个解，特解x，基础解系y Warning: Rank deficient, rank = 2, tol = 8.6112e-015. > In line_solution at 11 x = -2/11 10/11 0 0 y = 1/11 -9/11 -5/11 1/11 1 0 0 1 所以原方程组的通解是： ，其中 为任意常数。 5. 求代数方程的数值解。 (1) 3x+sinx-ex=0在x0=1.5附近的根。 (2) 在给定的初值x0=1，y0=1，z0=1下，求方程组的数值解。 解：M文件： function g=f(x) g=3*x+sin(x)-exp(x); clc;clear; fzero('f',1.5) 结果是： ans = 1289/682 (2). M文件： function F=fun(X) x=X(1); y=X(2); z=X(3); F(1)=sin(x)+y^2+log(z)-7; F(2)=3*x+2-z^3+1; F(3)=x+y+z-5; X=fsolve('myfun',[1,1,1]',optimset('Display','off')) 运行结果： 6. 求函数在指定区间的极值。 (1) 在(0,1)内的最小值。 (2) 在[0,0]附近的最小值点和最小值。 解：M文件： function f=g(u) x=u(1); y=u(2); f=2*x.^3+4*x.*y^3-10*x.*y+y.^2; clc;clear; format long f=inline('(x^3+cos(x)+x*log(x))/exp(x)'); [x,fmin1]=fminbnd(f,0,1) [U,fmin2]=fminsearch('g',[0,0]) 运行结果 7. 求微分方程的数值解。 解：M文件： function xdot= sys( x,y) xdot=[y(2);(5*y(2)-y(1))/x]; clc;clear; x0=1.0e-9;xf=20; [x,y]=ode45('sys',[x0,xf],[0 0]); [x,y] 运行结果： 8. 求微分方程组的数值解，并绘制解的曲线。 解： 令y1=x,y2=y,y3=z; 这样方程变为: ,自变量是t M文件： function xdot=sys(x,y) xdot=[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]; clc;clear; t0=0;tf=8; 12 [x,y]=ode23('sys',[t0,tf],[0,1,1]) plot(x,y)