1、河南城建学院第 1 页一 熟悉 Matlab 工作环境1、熟悉 Matlab 的 5 个基本窗口思考题：（1）变量如何声明，变量名须遵守什么规则、是否区分大小写。答：变量一般不需事先对变量的数据类型进行声明，系统会依据变量被赋值的类型自动进行类型识别，也就是说变量可以直接赋值而不用提前声明。变量名要遵守以下几条规则：变量名必须以字母开头，只能由字母、数字或下划线组成。变量名区分大小写。变量名不能超过 63 个字符。关键字不能作为变量名。最好不要用特殊常量作为变量名。（2）试说明分号、逗号、冒号的用法。分号：分隔不想显示计算结果的各语句；矩阵行与行的分隔符。逗号：分隔欲显示计算结果的各语句；变量
2、分隔符；矩阵一行中各元素间的分隔符。冒号：用于生成一维数值数组；表示一维数组的全部元素或多维数组某一维的全部元素。（3）linspace（）称为“线性等分”函数，说明它的用法。LINSPACE Linearly spaced vector. 线性等分函数LINSPACE(X1, X2) generates a row vector of 100 linearlyequally spaced points between X1 and X2.以 X1 为首元素，X2 为末元素平均生成 100 个元素的行向量。LINSPACE(X1, X2, N) generates N points betwe
3、en 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.1416答：3 次执行的结果不一样。exist（）函数是返回变量搜索顺序的一个函数。在第一次执行时返回 5 代表变量 pi 是由 Matlab 构建的变量。在第二次执行时已经通过赋值语句定义了变量 pi，返回 1 代表 pi 是工作空间变量。第三次执行前清除了工作空间，此时 pi 为系统默认常量，和第一次执行时性质一样，所以又返回 5。（2）圆周率 pi 是系
4、统默认常量，为什么会被改变为 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 由
5、 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河南城建学院第 3 页1 1 1 11 1 1 1B =0 0 00 0 00 0 0C =1 0 00 1 00 0 1D =8 1 63 5 74 9 2 E=C;D, F=E(2:5,:), G=reshape(F,3,4)E =1 0 00 1 00 0 18 1 63 5 74 9 2F =0 1 00 0 18 1
6、 63 5 7G =0 3 1 10 1 5 68 0 0 7 H=repmat(G,2)H =0 3 1 1 0 3 1 10 1 5 6 0 1 5 68 0 0 7 8 0 0 70 3 1 1 0 3 1 1河南城建学院第 4 页0 1 5 6 0 1 5 68 0 0 7 8 0 0 72）矩阵运算练习：1）用矩阵除法求下列方程组的解 x=1;2;3;61+32+43=321+52+73=481233=7 A=6 3 4;-2 5 7;8 -1 -3,B=3;-4;-7A =6 3 4-2 5 78 -1 -3B =3-4-7 x=ABx =1.0200-14.00009.72002
7、）求矩阵的秩； r=rank(A)r =33）求矩阵的特征值与特征向量 X,Lamda=eig(A)X =0.8013 -0.1094 -0.16060.3638 -0.6564 0.86690.4749 0.7464 -0.4719Lamda =9.7326 0 00 -3.2928 00 0 1.56024）矩阵的乘幂（平方）与开方 A2ans =62 29 33河南城建学院第 5 页34 12 626 22 34 A1=sqrtm(A)A1 =2.2447 + 0.2706i 0.6974 - 0.1400i 0.9422 - 0.3494i-0.5815 + 1.6244i 2.100
8、5 - 0.8405i 1.7620 - 2.0970i1.9719 - 1.8471i -0.3017 + 0.9557i 0.0236 + 2.3845i5）矩阵的指数与对数（以 e 为底） Ae=expm(A)Ae =1.0e+004 *1.0653 0.5415 0.63230.4830 0.2465 0.28760.6316 0.3206 0.3745 Ael=logm(A)Ael =1.7129 + 0.4686i 0.5305 - 0.2425i 0.5429 - 0.6049i1.1938 + 2.8123i 0.3658 - 1.4552i -0.5514 - 3.6305i
9、-0.0748 - 3.1978i 0.7419 + 1.6546i 1.8333 + 4.1282i6）矩阵的提取（取右上三角）与翻转（逆时针转 90 度） a=triu(A)a =6 3 40 5 70 0 -3 a1=rot90(A)a1 =4 7 -33 5 -16 -2 83、多维数组的创建及运算练习：创建三维数组 A，第一页为 ，第二页为 ，第三页为 。然后用1 34 2 1 22 1 3 57 1reshape 函数重排为数组 B，B 为 3 行、2 列、2 页。 a=1 3;4 2,b=1 2;2 1,c=3 5;7 1河南城建学院第 6 页 A=cat(3,a,b,c)A(:
10、,:,1) =1 34 2A(:,:,2) =1 22 1A(:,:,3) =3 57 1 B=reshape(A,3,2,2)B(:,:,1) =1 24 13 2B(:,:,2) =2 71 53 1三 Matlab 数值运算1、多项式运算练习：求 的商及余多项式。(2+1)(+3)(+1)3+2+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时
11、Y的值X 1 2 3 4 5 6 7 8 9 10Y 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)河南城建学院第 7 页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
12、),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 求方程 的根21=0 roots(1 -1 -1)ans =-0.61801.6180（2） ，在 n 个节点（n 不要太大，如取 511）上用分段线性和三=， 02次样条插值方法，计算 m 个插值点（m 可取 50100）的函数值。通过数值和图形输出，将两种插值结果与精度进行比较。适当增加 n，再作比较。 x=lin
13、space(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)河南城建学院第 8 页0 1 2 3 4 5 6 7-1-0.8-0.6-0.4-0.200.20.40.60.
14、81 e1=y1-y0;e2=y2-y0; plot(xi,e1)0 1 2 3 4 5 6 7-0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1 plot(xi,e2)河南城建学院第 9 页0 1 2 3 4 5 6 7-0.015-0.01-0.00500.0050.010.015（3）大气压强 p 随高度 x 变化的理论公式为 ，为验证这一公=1.0332(+500)7756式，测得某地大气压强随高度变化的一组数据如表所示。试用插值法和拟合法进行计算并绘图，看那种方法较为合理，且总误差最小。高度/m 0 300 600 1000 1500 2000压
15、强/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
16、.9322 0.8969 0.8519 0.7989 河南城建学院第 10 页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.8
17、116e-005河南城建学院第 11 页四 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 91Y1 44 45 47 50 50 38 30 30 3
18、4 36 34 41 45 46Y2 44 59 70 72 93 100 110 110 110 117 118 116 118 118X 96 101 104 106.5 111.5 118 123.5 136.5 142 146 150 157 158Y1 43 37 33 28 32 65 55 54 52 50 66 66 68Y2 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,1
19、18,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
20、 =3.3819e+004t2 =7.6328e+004t =河南城建学院第 12 页4.2510e+004 expt=t-41228expt =1.2819e+0032、习题（4）利用梯形法和辛普森法求定积分的值，并对结果进行比较。如果积分区间改为-55 结果有何不同？梯形积分中改变自变量 x 的维数，结果有何不同？123322 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
21、/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 位有效数字。102+4矩形法： x=linspace(0,1);y=x./(x.2+4); t=cumsum(y)*1/99;T=t(100)T =0.1126梯形法：河南城建学院第 13 页 x=linspac
22、e(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 创建符号表达式： ， 。1=+ sin22=2 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 个矩阵：=
23、123456789 = 5 92 6 103 7 114 8 126、符号表达式的变量替换练习：（1）已知 ,按照自变量x和自变量a，对=(2+3)3(2+41)表达式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)河南城建学院第 14 页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
24、-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）求lim22123+2 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*s
25、in(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）计算定积分 60（ +2） f=sym(sin(x)+2); int(f,x,0,pi/6)ans =-1/2*3(1/2)+1/3*pi+1（4）求下列线性方程组的解+=103+2+=142+3=1 f1=sym(x+y+z=10); f2=sym(3*x+2*y+z=14);河南城建学院第 15 页 f3=sym(2*x+3*y-z=1); g=solve(f1,f2
26、,f3,x,y,z)g = x: 1x1 symy: 1x1 symz: 1x1 sym g.xans =1 g.yans =2 g.zans =7（5）求解当y（0）=2，z（0）=7时，微分方程组的解=+=1+ 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=e
27、xp3(x)n=1;s=0;while n=10e-6k=k+1;if rem(k,2)=0jspi=jspi+1/i;elsejspi=jspi-1/i;endi=i+2;endp=4*jspi ,k 2、子函数和参数传递练习：编写求矩形面积函数rect，当没有输入参数时，显示提示信息；当只输入一个参数时，则以该参数作为正方形的边长计算其面积；当有两个参数时，则以这两个参数为长和宽计算其面积。function s=mianji(a,b)switch nargincase 0error(没有输入参数)case 1s=a*a;case 2s=a*b;end3、习题（3）编写一个函数project
28、1.m，其功能是判断某一年是否为闰年。function ryear(year)s=0;if rem(year,4)=0s=s+1;end河南城建学院第 17 页if rem(year,100)=0s=s-1;endif rem(year,400)=0s=s+1;endif s=1fprintf(%4d 是闰年.n,year)elsefprintf(%4d 不是闰年.n,year)end（4）编制一个函数，使得该函数能对输入的两个数值进行比较并返回其中的最小值。function c = bijiao(a,b)if nargin=2if a while var var=3; while var10
29、0var=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=0line(xdate,ydate,Col
30、or,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()河南城建学院第 19 页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） 。请=2(2+2)把 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)
