MATLAB上机答案.doc

2

一 熟悉Matlab工作环境
1、熟悉Matlab的5个基本窗口
思考题：
（1）变量如何声明，变量名须遵守什么规则、是否区分大小写。
答：变量一般不需事先对变量的数据类型进行声明，系统会依据变量被赋值的类型自动进行类型识别，也就是说变量可以直接赋值而不用提前声明。变量名要遵守以下几条规则：
Ø 变量名必须以字母开头，只能由字母、数字或下划线组成。
Ø 变量名区分大小写。
Ø 变量名不能超过63个字符。
Ø 关键字不能作为变量名。
Ø 最好不要用特殊常量作为变量名。
（2）试说明分号、逗号、冒号的用法。
分号：分隔不想显示计算结果的各语句；矩阵行与行的分隔符。
逗号：分隔欲显示计算结果的各语句；变量分隔符；矩阵一行中各元素间的分隔符。
冒号：用于生成一维数值数组；表示一维数组的全部元素或多维数组某一维的全部元素。
（3）linspace（）称为“线性等分”函数，说明它的用法。
LINSPACE Linearly spaced vector. 线性等分函数
LINSPACE(X1, X2) generates a row vector of 100 linearly
equally spaced points between X1 and X2.
以X1为首元素，X2为末元素平均生成100个元素的行向量。
LINSPACE(X1, X2, N) generates N points between X1 and X2.
For N < 2, LINSPACE returns X2.
以X1为首元素，X2为末元素平均生成n个元素的行向量。如果n＜2，返回X2。
Class support for inputs X1,X2:
float: double, single
数据类型：单精度、双精度浮点型。
（4）说明函数ones（）、zeros（）、eye（）的用法。
ones（）生成全1矩阵。
zeros（）生成全0矩阵。
eye（）生成单位矩阵。
2、Matlab的数值显示格式
思考题：
（1）3次执行exist（’pi’）的结果一样吗？如果不一样，试解释为什么？

>> pi
ans =
3.1416
>> sin(pi);
>> exist('pi')
ans =
5
>> pi=0;
>> exist('pi')
ans =
1
>> pi
pi =
0
>> clear
>> exist('pi')
ans =
5
>> pi
ans =
3.1416

19

答：3次执行的结果不一样。exist（）函数是返回变量搜索顺序的一个函数。在第一次执行时返回5代表变量pi是由Matlab构建的变量。在第二次执行时已经通过赋值语句定义了变量pi，返回1代表pi是工作空间变量。第三次执行前清除了工作空间，此时pi为系统默认常量，和第一次执行时性质一样，所以又返回5。
（2）圆周率pi是系统默认常量，为什么会被改变为0。
pi=0 为赋值语句，此时pi不再是系统默认常量，而是定义的变量了。
二 MATLAB语言基础
1、向量的生成和运算
练习：使用logspace（）创建1~4π的有10个元素的行向量。
>> A=logspace(0,1.0992,10)
A =
1.0000 1.3247 1.7550 2.3249 3.0799 4.0801 5.4051 7.1603 9.4856 12.5661

2、矩阵的创建、引用和运算
（1）矩阵的创建和引用
练习：创建以下矩阵：A为3×4的全1矩阵、B为3×3的0矩阵、C为3×3的单位矩阵、D为3×3的魔方阵、E由C和D纵向拼接而成、F抽取E的2~5行元素生成、G由F经变形为3×4的矩阵而得、以G为子矩阵用复制函数生成6×8的大矩阵H。
>> A=ones(3,4),B=zeros(3,3),C=eye(3,3),D=magic(3)

A =
1 1 1 1
1 1 1 1
1 1 1 1
B =
0 0 0
0 0 0
0 0 0
C =
1 0 0
0 1 0
0 0 1
D =
8 1 6
3 5 7
4 9 2
>> E=[C;D], F=E(2:5,:), G=reshape(F,3,4)
E =
1 0 0
0 1 0
0 0 1
8 1 6
3 5 7
4 9 2
F =
0 1 0
0 0 1
8 1 6
3 5 7
G =
0 3 1 1
0 1 5 6
8 0 0 7

>> H=repmat(G,2)
H =
0 3 1 1 0 3 1 1
0 1 5 6 0 1 5 6
8 0 0 7 8 0 0 7
0 3 1 1 0 3 1 1
0 1 5 6 0 1 5 6
8 0 0 7 8 0 0 7
2）矩阵运算

练习：1）用矩阵除法求下列方程组的解x=

>> A=[6 3 4;-2 5 7;8 -1 -3],B=[3;-4;-7]

A =
6 3 4
-2 5 7
8 -1 -3
B =
3
-4
-7

>> x=A\B
x =
1.0200
-14.0000
9.7200
2）求矩阵的秩；
>> r=rank(A)
r =
3
3）求矩阵的特征值与特征向量
>> [X,Lamda]=eig(A)

X =
0.8013 -0.1094 -0.1606
0.3638 -0.6564 0.8669
0.4749 0.7464 -0.4719
Lamda =
9.7326 0 0
0 -3.2928 0
0 0 1.5602

4）矩阵的乘幂（平方）与开方
>> A^2
ans =
62 29 33
34 12 6
26 22 34
>> A1=sqrtm(A)
A1 =
2.2447 + 0.2706i 0.6974 - 0.1400i 0.9422 - 0.3494i
-0.5815 + 1.6244i 2.1005 - 0.8405i 1.7620 - 2.0970i
1.9719 - 1.8471i -0.3017 + 0.9557i 0.0236 + 2.3845i
5）矩阵的指数与对数（以e为底）
>> Ae=expm(A)
Ae =
1.0e+004 *
1.0653 0.5415 0.6323
0.4830 0.2465 0.2876
0.6316 0.3206 0.3745

>> Ael=logm(A)
Ael =
1.7129 + 0.4686i 0.5305 - 0.2425i 0.5429 - 0.6049i
1.1938 + 2.8123i 0.3658 - 1.4552i -0.5514 - 3.6305i
-0.0748 - 3.1978i 0.7419 + 1.6546i 1.8333 + 4.1282i
6）矩阵的提取（取右上三角）与翻转（逆时针转90度）

>> a=triu(A)
a =
6 3 4
0 5 7
0 0 -3
>> a1=rot90(A)
a1 =
4 7 -3
3 5 -1
6 -2 8

3、多维数组的创建及运算

练习：创建三维数组A，第一页为，第二页为，第三页为。然后用reshape函数重排为数组B，B为3行、2列、2页。
>> a=[1 3;4 2],b=[1 2;2 1],c=[3 5;7 1]
>> A=cat(3,a,b,c)

A(:,:,1) =
1 3
4 2
A(:,:,2) =
1 2
2 1
A(:,:,3) =
3 5
7 1

>> B=reshape(A,3,2,2)

B(:,:,1) =
1 2
4 1
3 2
B(:,:,2) =
2 7
1 5
3 1

三 Matlab数值运算
1、多项式运算
练习：求的商及余多项式。
>> p1=conv([1 0 1],conv([1 3],[1 1]))
p1 =
1 4 4 4 3
>> [q r]=deconv(p1,[1 0 2 1])
q =
1 4
r =
0 0 2 -5 -1
2、多形式插值和拟合
有一组实验数据如附表1-1所示。请分别用拟合（二阶至三阶）和插值（线性和三次样条）的方法来估测X=9.5时Y的值
X
1
2
3
4
5
6
7
8
9
10
Y
16
32
70
142
260
436
682
1010
1432
1960
>> x=1:10;y=[16 32 70 142 260 436 682 1010 1432 1960];

>> p1=polyfit(x,y,1)
p1 =
204.8000 -522.4000
>> y1=polyval(p1,9.5)
y1 =
1.4232e+003
>> p2=polyfit(x,y,2),y2=polyval(p2,9.5)
p2 =
32.0000 -147.2000 181.6000
y2 =
1.6712e+003

>> p3=polyfit(x,y,3),y3=polyval(p3,9.5)
p3 =
2.0000 -1.0000 5.0000 10.0000
y3 =
1.6820e+003
>> y4=interp1(x,y,9.5)
y4 =
1696
>> y5=spline(x,y,9.5)
y5 =
1682

3、习题
（1）用函数roots求方程的根
>> roots([1 -1 -1])
ans =
-0.6180
1.6180
（2），在n个节点（n不要太大，如取5~11）上用分段线性和三次样条插值方法，计算m个插值点（m可取50~100）的函数值。通过数值和图形输出，将两种插值结果与精度进行比较。适当增加n，再作比较。
>> x=linspace(0,2*pi,8),y=sin(x)
x =
0 0.8976 1.7952 2.6928 3.5904 4.4880 5.3856 6.2832
y =
0 0.7818 0.9749 0.4339 -0.4339 -0.9749 -0.7818 -0.0000
>> xi=linspace(0,2*pi,100);y0=sin(xi);y1=interp1(x,y,xi);y2=interp1(x,y,xi,'spline');
>> plot(xi,y0,'*',xi,y1,'-.',xi,y2)

>> e1=y1-y0;e2=y2-y0;
>> plot(xi,e1)

>> plot(xi,e2)

（3）大气压强p随高度x变化的理论公式为，为验证这一公式，测得某地大气压强随高度变化的一组数据如表所示。试用插值法和拟合法进行计算并绘图，看那种方法较为合理，且总误差最小。
高度/m
0
300
600
1000
1500
2000
压强/Pa
0.9689
0.9322
0.8969
0.8519
0.7989
0.7491
插值法：
>> x=[0 300 600 1000 1500 2000]; p=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491];
>> xi=linspace(0,2000);p0=1.0332*exp(-(xi+500)/7756);
>> p1=interp1(x,p,xi,'spline');
>> plot(xi,p0,'*',xi,p1)
>> e1=p1-p0;
>> e=sum(e1.^2)
e =
1.8652e-005
拟合法：
>> x=[0 300 600 1000 1500 2000]; p=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491];
>> P=log10(p)
P =
-0.0137 -0.0305 -0.0473 -0.0696 -0.0975 -0.1255
>> p1=polyfit(x,P,1)
p1 =
-0.0001 -0.0137
>> b=p1(1)/0.4343,a=10.^p1(2)
b =
-1.2863e-004
a =
0.9689
>> xi=linspace(0,2000);p0=1.0332*exp(-(xi+500)/7756);
>> p2=polyval(p1,xi);P2=10.^p2;
>> e2=P2-p0;e=sum(e2.^2)
e =
1.8116e-005

四 Matlab数值运算
1、数值微积分
练习：瑞士地图如图所示，为了算出其国土面积，首先对地图作如下测量：以由西向东方向为X轴，由南到北方向为Y轴，选择方便的原点，并将从最西边界点到最东边界点在X轴上的区间适当划分为若干段，在每个分点的Y方向测出南边界点和北边界点的Y坐标Y1和Y2，根据地图比例尺知道18mm相当于40km，试由测量数据计算瑞士国土近似面积，与其精确值41228km2比较。
X
7
10.5
13
17.5
34
40.5
44.5
48
56
61
68.5
76.5
80.5
91
Y1
44
45
47
50
50
38
30
30
34
36
34
41
45
46
Y2
44
59
70
72
93
100
110
110
110
117
118
116
118
118

X
96
101
104
106.5
111.5
118
123.5
136.5
142
146
150
157
158
Y1
43
37
33
28
32
65
55
54
52
50
66
66
68
Y2
121
124
121
121
121
116
122
83
81
82
86
85
68

>> x=[7,10.5,13,17.5,34,40.5,44.5,48,56,61,68.5,76.5,80.5,91,96,101,104,106.5,111.5,118,123.5,136.5,142,146,150,157,158];
>> y1=[44,45,47,50,50,38,30,30,34,36,34,41,45,46,43,37,33,28,32,65,55,54,52,50,66,66,68];
>> y2=[44,59,70,72,93,100,110,110,110,117,118,116,118,118,121,124,121,121,121,116,122,83,81,82,86,85,68];
>> X=x./18*40;Y1=y1./18*40;Y2=y2./18*40;
>> t1=trapz(X,Y1),t2=trapz(X,Y2), t=t2-t1
t1 =
3.3819e+004
t2 =
7.6328e+004
t =
4.2510e+004
>> expt=t-41228
expt =
1.2819e+003
2、习题
（4）利用梯形法和辛普森法求定积分的值，并对结果进行比较。如果积分区间改为-5~5结果有何不同？梯形积分中改变自变量x的维数，结果有何不同？

>> x=linspace(-3,3);y=exp(-x.^2/2);

>> t=(1/2*pi)*trapz(x,y)
t =
3.9267
>> q=(1/2*pi)*quad('exp(-x.^2/2)',-3,3)
q =
3.9268

>> x=linspace(-5,5);y=exp(-x.^2/2);

>> t=(1/2*pi)*trapz(x,y)
t =
3.9374
>> q=(1/2*pi)*quad('exp(-x.^2/2)',-5,5)
q =
3.9374

>> x=linspace(-3,3,150);y=exp(-x.^2/2);
>> t=(1/2*pi)*trapz(x,y)
t =
3.9268
（5）分别用矩形法、梯形法、辛普森法和牛顿-科茨4种方法近似计算定积分，取n=4，保留4位有效数字。

矩形法：
>> x=linspace(0,1);y=x./(x.^2+4);
>> t=cumsum(y)*1/99;T=t(100)
T =
0.1126
梯形法：
>> x=linspace(0,1);y=x./(x.^2+4);
>> t=trapz(x,y)
t =
0.1116

辛普森法：
>> q=quad('x./(x.^2+4)',0,1)
q =
0.1116
牛顿-科茨法：
>> q=quadl('x./(x.^2+4)',0,1)
q =
0.1116

五 Matlab符号运算
1、符号矩阵创建

练习：分别用sym和syms创建符号表达式：，。
>> f1=sym('cos(x)+(-(sin(x)^2))^(1/2)')
f1 =cos(x)+(-(sin(x)^2))^(1/2)

>> syms y e t
>> f2=y/exp(-2*t)
f2 =y/exp(-2*t)
2、习题
（2）试创建以下2个矩阵：

6、符号表达式的变量替换
练习：（1）已知,按照自变量x和自变量a，对表达式f分别进行降幂排列。
>> f=sym('(a*x^2+b*x+c-3)^3-a*(c*x^2+4*b*x-1)')
f =(a*x^2+b*x+c-3)^3-a*(c*x^2+4*b*x-1)
>> f1=collect(f),f2=collect(f,'a')
f1 =
a^3*x^6+3*b*a^2*x^5+((c-3)*a^2+2*b^2*a+a*(2*(c-3)*a+b^2))*x^4+(4*(c-3)*b*a+b*(2*(c-3)*a+b^2))*x^3+((c-3)*(2*(c-3)*a+b^2)+2*b^2*(c-3)+a*(c-3)^2-a*c)*x^2+(3*(c-3)^2*b-4*b*a)*x+(c-3)^3+a
f2 =
a^3*x^6+3*(b*x+c-3)*x^4*a^2+(3*(b*x+c-3)^2*x^2-c*x^2-4*b*x+1)*a+(b*x+c-3)^3
8、符号方程的求解
练习：（1）求
>> f=sym('(x^2-1)/(x^2-3*x+2)');
>> limit(f,'x',2)
ans =NaN
（2）求函数f（x）=cos2x-sin2x的积分；求函数的导数。
>> f=sym('cos(2*x)-sin(2*x)');
>> int(f)
ans =1/2*sin(2*x)+1/2*cos(2*x)

>> g=sym('(exp(x)+x*sin(x))^(1/2)');
>> diff(g)
ans =1/2/(exp(x)+x*sin(x))^(1/2)*(exp(x)+sin(x)+x*cos(x))
（3）计算定积分
>> f=sym('sin(x)+2');
>> int(f,'x',0,pi/6)
ans =-1/2*3^(1/2)+1/3*pi+1
（4）求下列线性方程组的解

>> f1=sym('x+y+z=10');
>> f2=sym('3*x+2*y+z=14');
>> f3=sym('2*x+3*y-z=1');
>> g=solve(f1,f2,f3,'x','y','z')
g =
x: [1x1 sym]
y: [1x1 sym]
z: [1x1 sym]
>> g.x
ans =1
>> g.y
ans =2
>> g.z
ans =7
（5）求解当y（0）=2，z（0）=7时，微分方程组的解

>> [g_y,g_z]=dsolve('Dy-z=sin(x)','Dz+y=1+x','y(0)=2','z(0)=7','x')
g_y =cos(x)+6*sin(x)+1/2*sin(x)*x+1+x
g_z =-3/2*sin(x)+6*cos(x)+1+1/2*cos(x)*x
六 Matlab程序设计
1、程序流程控制结构
练习：（1）请把exp2.函数文件用while循环改写。
function s=exp3(x)
n=1;s=0;
while n<=x
s=s+n;
n=n+1;
end
s
（2）用公式求pi的近似值，直到最后一项的绝对值小于10-6为止，试编写其M脚本文件。
k=0;jspi=1;i=3;
while (1/i)>=10e-6
k=k+1;
if rem(k,2)==0
jspi=jspi+1/i;
else
jspi=jspi-1/i;
end
i=i+2;
end
p=4*jspi ,k
2、子函数和参数传递
练习：编写求矩形面积函数rect，当没有输入参数时，显示提示信息；当只输入一个参数时，则以该参数作为正方形的边长计算其面积；当有两个参数时，则以这两个参数为长和宽计算其面积。
function s=mianji(a,b)
switch nargin
case 0
error('没有输入参数')
case 1
s=a*a;
case 2
s=a*b;
end
3、习题
（3）编写一个函数project1.m，其功能是判断某一年是否为闰年。
function ryear(year)
s=0;
if rem(year,4)==0
s=s+1;
end
if rem(year,100)==0
s=s-1;
end
if rem(year,400)==0
s=s+1;
end
if s==1
fprintf('%4d 是闰年.\n',year)
else
fprintf('%4d 不是闰年.\n',year)
end
（4）编制一个函数，使得该函数能对输入的两个数值进行比较并返回其中的最小值。
function c = bijiao(a,b)
if nargin==2
if a < b
c=a;
else
c=b;
end
else
error('输入参数不正确')
end
（6）观察以下循环语句，计算每个循环的循环次数和循环结束之后var的值。
var=1;
while mod(var,10)~=0
var=var+1
end
循环次数10，var=10。
var=2;
>> while var<=100
var=var^2;
end
循环次数4，var=256。
>> var=3;
>> while var>100
var=var^2;
end
循环次数0次，var=3。
七 Matlab数据可视化
1、二维图形绘制
练习：写出图A2的绘制方法。
y1=sin(x);y2=cos(x);
plot(x,y1,'r -',x,y2,'m --')
x=linspace(0,4*pi);
y1=sin(x);y2=cos(x);
plot(x,y1,'r --',x,y2,'m -')
ylabel('·ù¶È'),xlabel('ʱ¼ä')
legend('sinx','cosx')
gtext('\leftarrowsinx')
gtext('\leftarrowcosx')
axis([0 16 -1 1])
xdate=0:0.5:16;ydate=0
line(xdate,ydate,'Color','k','Marker','.')
2、三维曲线和三维曲面绘制
练习：（1）绘制以上空间螺旋线的俯视图、左侧视图和前视图。
z=0:0.1:6*pi;
x=cos(z);y=sin(z);
subplot(2,2,1);plot3(cos(z),sin(z),z);
title('ÈýάÇúÏß')
subplot(2,2,2);plot(cos(z),sin(z));
title('¸©ÊÓͼ')
subplot(2,2,3);plot(cos(z),z);
title('×óÊÓͼ')
subplot(2,2,4);plot(sin(z),z);
title('Ç°ÊÓͼ')
（2）设，求定义域x=[-2，2]，y=[-2，2]内的z值（网格取0.1）。请把z的值用网面图形象的表示出来，如图A3所示。
x=-2:0.1:2;y=x;
[X,Y]=meshgrid(x,y);
Z=X.^2.*exp(-(X.^2+Y.^2));
surf(X,Y,Z)