河南城建学院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

答：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 >> 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.65 0.8669 0.4749 0.74 -0.4719 4）矩阵的乘幂（平方）与开方 >> A^2 ans =

62 29 33 34 12 6

B = 3 -4 -7

Lamda =

9.7326 0 0 0 -3.2928 0 0 0 1.5602 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.15 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.29 - 0.6049i 1.1938 + 2.8123i 0.3658 - 1.4552i -0.5514 - 3.6305i -0.0748 - 3.1978i 0.7419 + 1.66i 1.8333 + 4.1282i 6）矩阵的提取（取右上三角）与翻转（逆时针转90度） >> a=triu(A) >> a1=rot90(A) a =

a1 =

6 3 4 4 7 -3 0 5 7 3 5 -1 0 0 -3 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

>> B=reshape(A,3,2,2) B(:,:,1) = 1 2 4 1 3 2

A(:,:,2) = 1 2 2 1

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

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 Y 1 16 2 32 3 70 4 142 5 260 6 436 7 682 8 1010 9 1432 10 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

3、习题

（1）用函数roots求方程>> roots([1 -1 -1]) ans = -0.6180 1.6180 （2）

的根

y3 =

1.6820e+003

>> y4=interp1(x,y,9.5) y4 = 1696

>> y5=spline(x,y,9.5) y5 = 1682

，在n个节点（n不要太大，如取5~11）上用分段线性和三次

样条插值方法，计算m个插值点（m可取50~100）的函数值。通过数值和图形输出，将两种插值结果与精度进行比较。适当增加n，再作比较。

>> x=linspace(0,2*pi,8),y=sin(x) x =

0 0.76 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)

10.80.60.40.20-0.2-0.4-0.6-0.8-101234567

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

0.10.080.060.040.020-0.02-0.04-0.06-0.08-0.101234567

>> plot(xi,e2)

0.0150.010.0050-0.005-0.01-0.01501234567

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

，为验证这一公式，

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

高度/m 压强/Pa 0 0.96 300 0.9322 600 0.69 1000 0.8519 1500 0.79 2000 0.7491 插值法： >> x=[0 300 600 1000 1500 2000]; p=[0.96 0.9322 0.69 0.8519 0.79 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.96 0.9322 0.69 0.8519 0.79 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.96

>> 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，试由测量数据计算瑞士国土近似面积，与其精确值41228km比较。 X Y1 Y2 X Y1 Y2 >>

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,,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 = 96 43 101 104 106.5 111.5 118 123.5 136.5 142 146 150 157 158 37 33 28 32 121 65 55 83 52 81 50 82 66 86 66 85 68 68 7 44 44 10.5 13 45 59 47 70 17.5 34 50 72 50 93 40.5 44.5 48 38 30 30 56 34 61 36 68.5 76.5 80.5 91 34 41 45 46 2

100 110 110 110 117 118 116 118 118 121 124 121 121 116 122 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)

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

>> 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 =

，

q = 3.9374 q = 3.9268

0.1126 梯形法：

>> x=linspace(0,1);y=x./(x.^2+4); 辛普森法：

>> q=quad('x./(x.^2+4)',0,1) q = 0.1116

>> t=trapz(x,y) t = 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）已知达式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)

,按照自变量x和自变量a，对表

>> 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）用

为止，试编写其M脚本文件。

公式求pi的近似值，直到最后一项的绝对值小于10

-6

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)