河南城建学院MATLAB上机实验答案(新)

一熟悉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 =

>> pi

pi =

>> clear

>> exist('pi')

ans =

>> pi

ans =

3.1416

答：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

8 -1 -3 B =

3 -

4 -7

>> x=A\B

x =

1.0200

-14.0000

9.7200

2）求矩阵的秩；

>> r=rank(A)

r =

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: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)

-0.1

-0.08-0.06-0.04-0.0200.020.040.060.080.1

>> plot(xi,e2)

01234567

-0.015

-0.01

-0.005

0.005

0.01

0.015

（3）大气压强p 随高度x 变化的理论公式为

，为验证这一公式，

测得某地大气压强随高度变化的一组数据如表所示。试用插值法和拟合法进行计算并绘图，看那种方法较为合理，且总误差最小。

>> 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,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,6 6,68];

y2=[44,59,70,72,93,100,110,110,110,117,118,116,118,118,121,124,121,121,121,116,1 22,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

（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('èy???ú??')

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)