MATLAB程序设计与应用(第二版)
实验参考答案
%实验一MATLAB运算基础
%第一题
%(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%上面依次获得(3:5,2:3)%第四题
%(1):
A=100:999;
B=rem(A,21);%求余
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]);%对角矩阵[2,0;0,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))%产生5*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 %A的逆*b
%(2):
B=[0.95,0.67,0.53];
x=inv(A)*B
%(3):
cond(A)
%第五题方法二:
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=[]; %建立存放所有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
%第二题程序一:
x=input('请输入一个百分制成绩:'); if x>100|x<0
disp('您输入的成绩不是百分制成绩,请重新输入。');
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
%第二题程序二:
s=input('请输入一个成绩(0分到100分之间):'); %s用于存放成绩
while
1 %判断输入成绩的合理性
if s<0|s>100
disp('输入的成绩需在0到100之间,请重新输入:')
s=input('请输入一个成绩(0分到100分之间):');
else
break;
%第二题程序二:
s=input('请输入一个成绩(0分到100分之间):'); %s用于存放成绩
while
1 %判断输入成绩的合理性
if s<0|s>100
disp('输入的成绩需在0到100之间,请重新输入:')
s=input('请输入一个成绩(0分到100分之间):');
else
break;
%第三题
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=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) %产生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;
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);
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 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,co s(t*pi/180),0;0,m2,-sin(t*pi/18 0),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)%x为数
a=[1,2;3,4];%x为矩阵
y=fx(a)
%第五题
%(1)
function f1=mat5(n)
f1=n+10*log(n*n+5);
%在命令窗口中调用该函数文件:
y=mat5(40)/(mat5(30)+mat5(20)) %(2)方法一
function f2=mat6(n)
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,'\upar row y3=y1*y2');
%(2)
x=linspace(-2*pi,2*pi,100);
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','co lor','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(2 000*pi*t-0.2*x+pi)');
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))
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(beta 1)*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');
f=log(x);
taylor(f,x,6,1)
运行结果:
ans =
x - (x - 1)^2/2 + (x - 1)^3/3 - (x - 1)^4/4 + (x - 1)^5/5 - 1
3. 求下列方程的符号解。
解:
M文件:
clear all;clc;
x1=solve('log(x+1)-5/(1+sin(x))=2')
x2=solve('x^2+9*sqrt(x+1)-1')
x3=solve('3*x*exp(x)+5*sin(x)-78.5')
[x4
y4]=solve('sqrt(x^2+y^2)-100','3*x+5*y-8')
运行结果:
x1 =
521.67926389905839979437366649258
x2 =
-1
(3^(1/2)*i*(4/(9*(6465^(1/2)/2 + 2171/54)^(1/3)) - (1/2*6465^(1/2) + 2171/54)^(1/3)))/2 - (6465^(1/2)/2 + 2171/54)^(1/3)/2 - 2/(9*(6465^(1/2)/2 + 2171/54)^(1/3)) + 1/3
1/3 - (6465^(1/2)/2 + 2171/54)^(1/3)/2 - (3^(1/2)*i*(4/(9*(6465^(1/2)/2 + 2171/54)^(1/3)) - (1/2*6465^(1/2) + 2171/54)^(1/3)))/2 - 2/(9*(6465^(1/2)/2 + 2171/54)^(1/3))
x3 =
2.3599419584772910151699327715486
x4 =
12/17 - (10*21246^(1/2))/17
(10*21246^(1/2))/17 + 12/17
y4 =
(6*21246^(1/2))/17 + 20/17
20/17 - (6*21246^(1/2))/17
4. 求微分方程初值问题的符号解,并与数值解进行比较。
解:
M文件:
clear all;clc;
dsolve('D2y+4*Dy+29*y','y(0)=0','Dy(0)=15',' x')
运行结果:
ans =
(3*sin(5*x))/exp(2*x)
5. 求微分方程组的通解。
解:
M文件:
clear all;clc;
[x y z]=dsolve('Dx=2*x-3*y+3*z',...
'Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t') 运行结果:
x =
C1/exp(t) + C2*exp(2*t)
y =
C1/exp(t) + C2*exp(2*t) + C3/exp(2*t)
z = C2*exp(2*t) + C3/exp(2*t)
实验九数值微积分与方程数值求解
1. 求函数在指定点的数值导数。
解:M文件:
clc;clear;
x=1;
i=1;
f=inline('det([x x^2 x^3;1 2*x 3*x^2;0 2 6*x])');
while x<=3.01
g(i)=f(x);
i=i+1;
x=x+0.01; %以0.01的步长增加,可再缩小步长提高精度
end
g;
t=1:0.01:3.01;
dx=diff(g)/0.01; %差分法近似求导
f1=dx(1) %x=1的数值倒数
f2=dx(101) %x=2的数值倒数
f3=dx(length(g)-1) %x=3的数值倒数
运行结果:
f1 =
6.0602
f2 =
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=[ ];
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 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','o ff')) 运行结果: 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; [x,y]=ode23('sys',[t0,tf],[0,1,1]) plot(x,y)