MATLAB上机答案

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 pians = 3.1416 sin(pi); exist(pi)ans = 5 pi=0; exist(pi)ans = 1 pipi = 0 clear exist(pi)ans = 5 pians = 3.141619答：3次执行的结果不一样。exist（）函数是返回变量搜索顺序的一个函数。在第一次执行时返回5代表变量pi是由Matlab构建的变量。在第二次执行时已经通过赋值语句定义了变量pi，返回1代表pi是工作空间变量。第三次执行前清除了工作空间，此时pi为系统默认常量，和第一次执行时性质一样，所以又返回5。（2）圆周率pi是系统默认常量，为什么会被改变为0。pi=0 为赋值语句，此时pi不再是系统默认常量，而是定义的变量了。二 MATLAB语言基础1、向量的生成和运算练习：使用logspace（）创建14的有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.56612、矩阵的创建、引用和运算（1）矩阵的创建和引用练习：创建以下矩阵：A为34的全1矩阵、B为33的0矩阵、C为33的单位矩阵、D为33的魔方阵、E由C和D纵向拼接而成、F抽取E的25行元素生成、G由F经变形为34的矩阵而得、以G为子矩阵用复制函数生成68的大矩阵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 1B = 0 0 0 0 0 0 0 0 0C = 1 0 0 0 1 0 0 0 1D = 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 2F = 0 1 0 0 0 1 8 1 6 3 5 7G = 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 72）矩阵运算练习：1）用矩阵除法求下列方程组的解x= A=6 3 4;-2 5 7;8 -1 -3,B=3;-4;-7A = 6 3 4 -2 5 7 8 -1 -3B = 3 -4 -7 x=ABx = 1.0200 -14.0000 9.72002）求矩阵的秩； r=rank(A)r = 33）求矩阵的特征值与特征向量 X,Lamda=eig(A)X = 0.8013 -0.1094 -0.1606 0.3638 -0.6564 0.8669 0.4749 0.7464 -0.4719Lamda = 9.7326 0 0 0 -3.2928 0 0 0 1.56024）矩阵的乘幂（平方）与开方 A2ans = 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.3845i5）矩阵的指数与对数（以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.1282i6）矩阵的提取（取右上三角）与翻转（逆时针转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 83、多维数组的创建及运算练习：创建三维数组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 2A(:,:,2) = 1 2 2 1A(:,:,3) = 3 5 7 1 B=reshape(A,3,2,2)B(:,:,1) = 1 2 4 1 3 2B(:,:,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 4r = 0 0 2 -5 -12、多形式插值和拟合有一组实验数据如附表1-1所示。请分别用拟合（二阶至三阶）和插值（线性和三次样条）的方法来估测X=9.5时Y的值X12345678910Y163270142260436682101014321960 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.6000y2 = 1.6712e+003 p3=polyfit(x,y,3),y3=polyval(p3,9.5)p3 = 2.0000 -1.0000 5.0000 10.0000y3 = 1.6820e+003 y4=interp1(x,y,9.5)y4 = 1696 y5=spline(x,y,9.5)y5 = 16823、习题（1）用函数roots求方程的根 roots(1 -1 -1)ans = -0.6180 1.6180（2），在n个节点（n不要太大，如取511）上用分段线性和三次样条插值方法，计算m个插值点（m可取50100）的函数值。通过数值和图形输出，将两种插值结果与精度进行比较。适当增加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.2832y = 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变化的理论公式为，为验证这一公式，测得某地大气压强随高度变化的一组数据如表所示。试用插值法和拟合法进行计算并绘图，看那种方法较为合理，且总误差最小。高度/m0300600100015002000压强/Pa0.96890.93220.89690.85190.79890.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-004a = 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比较。X710.51317.53440.544.548566168.576.580.591Y14445475050383030343634414546Y24459707293100110110110117118116118118X96101104106.5111.5118123.5136.5142146150157158Y143373328326555545250666668Y2121124121121121116122838182868568 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-t1t1 = 3.3819e+004t2 = 7.6328e+004t = 4.2510e+004 expt=t-41228expt = 1.2819e+0032、习题（4）利用梯形法和辛普森法求定积分的值，并对结果进行比较。如果积分区间改为-55结果有何不同？梯形积分中改变自变量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*x2+b*x+c-3)3-a*(c*x2+4*b*x-1)f =(a*x2+b*x+c-3)3-a*(c*x2+4*b*x-1) f1=collect(f),f2=collect(f,a)f1 =a3*x6+3*b*a2*x5+(c-3)*a2+2*b2*a+a*(2*(c-3)*a+b2)*x4+(4*(c-3)*b*a+b*(2*(c-3)*a+b2)*x3+(c-3)*(2*(c-3)*a+b2)+2*b2*(c-3)+a*(c-3)2-a*c)*x2+(3*(c-3)2*b-4*b*a)*x+(c-3)3+af2 =a3*x6+3*(b*x+c-3)*x4*a2+(3*(b*x+c-3)2*x2-c*x2-4*b*x+1)*a+(b*x+c-3)38、符号方程的求解练习：（1）求 f=sym(x2-1)/(x2-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.xans =1 g.yans =2 g.zans =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+xg_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=10e-6 k=k+1; if rem(k,2)=0 jspi=jspi+1/i; else jspi=jspi-1/i; end i=i+2;endp=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;end3、习题（3）编写一个函数project1.m，其功能是判断某一年是否为闰年。function ryear(year)s=0;if rem(year,4)=0 s=s+1;endif rem(year,100)=0 s=s-1;endif rem(year,400)=0 s=s+1;endif s=1 fprintf(%4d 是闰年.n,year)else fprintf(%4d 不是闰年.n,year)end（4）编制一个函数，使得该函数能对输入的两个数值进行比较并返回其中的最小值。function c = bijiao(a,b)if nargin=2 if a while var var=3; while var100var=var2;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)