第8章MATLAB软件入门.doc

第8章 Matlab软件入门计算机技术的迅速发展，产生了许多功能强大的数学软件，为应用数学和工程技术的发展提供了威力无比的武器和工具，本章介绍由美国MathWorks 公司开发的集数值计算、符号计算和图形可视化三大基本功能于一体的，功能强大、简单易学的软件MATLAB。它是国际公认的优秀数学应用软件之一。从20世纪80年代产生到现在，MATLAB已经发展成为适合多学科的大型软件，在世界各高校，MATLAB已经成为线性代数、数值分析、数理统计、优化方法、自动控制、数字信号处理、动态系统仿真等高级课程的基本教学工具。特别是最近几年，MATLAB在我国大学生数学建模竞赛中的应用，为参赛者在有限的时间内准确、有效的解决问题提供了有力的保证。本章以MATLAB6.5.1版本为例，主要通过举例介绍MATLAB的初步应用。先简单介绍MATLAB的工作界面，启动MATLAB软件将看到以下界面图8-1命令窗口（Command Window）：其中“”为命令提示符，表示MATLAB已经准备好，可以接受用户在此输入命令,命令执行的结果也显示在这个窗口命令历史窗口（Command History）：已经执行过的命令依次显示在这个窗口中，可以备查工作区窗口（Workspace）：显示当前计算机内存中有哪些变量及相关信息新建m文档按钮：点此按钮会出现如下图窗口，为了解决某个问题我们可以在此窗口中编写程序，命名、保存后让计算机运行，这样的文档叫m文档8.1 基本运算例1：求输入(12+2*(7-4)/32 执行，结果为ans = 2 （系统默认变量名为ans）或输入x=(12+2*(7-4)/32执行，结果为x = 2x为变量名，命名变量可以用一个字母，也可用多个字母，中间可以有数字，但不能以数字开头，如a，dgs，s1，ff2等都是合法变量名；区分大小写，如a，A就是两个不同的变量。例2：计算 形式1： 56+sin(pi)+exp(3) ans = 1.5645e+004形式2： a=56+sin(pi)+exp(3)a = 1.5645e+004如果在表达式后加分号“；”，即 a=56+sin(pi)+exp(3)；运算后不显示运算结果。表8-1 常用运算符关系运算符功 能+加号-减号*乘号/除法乘幂= =等于 =不等于= =等于=大于等于&逻辑与|逻辑或逻辑非 表8-2 常用函数名称含义sin(x)x的正弦cos(x)x的余弦tan(x)x的正切cot(x)x的余切asin(x)x的反正弦acos(x)x的反余弦exp(x)pow2(x)log(x)log2(x)abs(x)x的绝对值sqrt(x)x的平方根round(x)四舍五入到最近整数rem(x,y)x除以y的余数gcd(x,y)x和y的最大公因数lcm(x,y)x和y的最小公倍数例3：执行 y=sin(10)*exp(-0.3*42) 结果 y = -0.0045系统默认4位有效数字，若想提高精度可用命令：vpa(y,10) 结果为 -.4477156132e-2，即执行就可精确到10位有效数字例4： x=sqrt(2)/2 % “”为命令提示符，不用输入x = 0.7071 y=asin(x) % x的反正弦y = 0.7854 % 此处y为弧度，系统默认 y_deg=y*180/pi % 将弧度转化为角度，y_deg为变量名y_deg = 45.0000例5： Inf、-Inf表示正负无穷大，例如以下命令 2/0ans = Inf log(0)Warning: Log of zero.ans = -Inf exp(3000)ans = Inf % 超过计算机能表示的最大的数，故显示为无穷大 NaN表示无意义，Not a Nunber的缩写，例如 0/0ans = NaN8.2 矩阵与数组8.2.1 输入矩阵例1：输入一维矩阵（即数组或向量）a=1,4,6,8,10，运行a = 1 4 6 8 10例2：要得到矩阵，输入A = 1,4,3; 6,5,-2; 7,8,5 执行， A = 1 4 3 6 5 -2 7 8 5例3：MATLAB中有一些函数可用来直接生成特殊矩阵，如 u=zeros(2,3) %生成23全零矩阵u = 0 0 0 0 0 0 v=ones(3) %生成33全1矩阵v = 1 1 1 1 1 1 1 1 1 w=eye(3) %生成33对角线为1的单位阵w = 1 0 0 0 1 0 0 0 1 x=(1:2:121); % 起始值为1,每次增加2,终止值为121的向量，去掉分号可在屏幕上显示结果 x=(1:121); % 起始值为1,系统默认每次增加1,终止值为121的向量 x=linspace(0,1,101); % 利用linspace，生成以0为起始值,1为终止值,元素数目为101的向量8.2.2对矩阵中的元素可以进行访问例4： A = 1,4,3; 6,5,-2; 7,8,5 % 先生成矩阵A s=A(2,3) % s为A第二行第三列的元素s =-2 s1=A(3,:) % s1为A的第三行s1 = 7 8 5 s2=A(:,2) % s2为A的第二列s2 = 4 5 8 s3=A(1:2,2:3) % s3取A的一、二行，二、三列s3 = 4 3 5 9注：A输入后一直保存在内存工作区（工作空间，Workspace）中,同时会显示矩阵的大小（几行几列），工作区内的的变量可以随时调用，除非被清除（用命令clear清除）矩阵中的元素可以修改，如 A(2,3)=9A = 1 4 3 6 5 9 7 8 58.2.3 矩阵的拼接例5： a=1:7，b=11:2:23 % 先生成矩阵a及ba = 1 2 3 4 5 6 7b = 11 13 15 17 19 21 23 c=b a % 利用阵列 a 及b ，组成新矩阵cc =11 13 15 17 19 21 23 1 2 3 4 5 6 7 d=b; a % 利用a及b，组成新矩阵dd = 11 13 15 17 19 21 23 1 2 3 4 5 6 7 f=b,a % b表示b的转置，将行变成列f = 11 1 13 2 15 3 17 4 19 5 21 6 23 7 g=f(:,1),zeros(7,1),ones(7,2)g = 11 0 1 1 13 0 1 1 15 0 1 1 17 0 1 1 19 0 1 1 21 0 1 1 23 0 1 1例6: 输入矩阵x=4,8,12,10,23;6,3,15,13,19;9,1,2,18,14;11,7,5,21,17运行以下命令，观察结果m=max(x) % x每列的最大值,如果x只是一行（或一列），则m就是这一行（或一列）向量的最大值n=min(x) % 每列的最小值（问：如何得到整个矩阵的最小值与最大值？）m,n=size(x) % 给出x的大小，即行数和列数l=length(x) % 给出x的列数，如果x是一行或一列向量，则给出这个向量的元素的个数，即长度y=x % 矩阵x的转置，即行与列互换d=x(1:3,3:5) % x的第1到3行，第3到5列送给d，产生一个新矩阵y(2,3)=y(2,3)/2 % y的第二行第三列的元素除以2，再替换原来的数y(2,:)=y(2,:)/2 %第二行的元素全部除以二y(:,4)=y(:,4)+y(:,2) %第四列加第二列替换原来的第四列矩阵的函数运算是矩阵运算中最实用的部分，常用的主要有以下几个：det(a) 求矩阵a的行列式eig(a) 求矩阵a的特征值inv(a)或a(-1) 求矩阵a的逆矩阵rank(a) 求矩阵a的秩trace(a) 求矩阵a的迹（对角线元素之和）8.3 函数作图下面画函数图像，最好打开m文档，在文档里做，在命令窗口里做不方便，见8.4节程序设计8.3.1二维图形例1：下面画sin(x)的图像，matlab画图的基本原理是描点连线x=0:0.1:10; % 0到10的101个点(每隔0.1画一个点)的x座标y=sin(x); % 对应的y座标，也是101个点plot(x,y); % 以x为自变量，y为函数值画函数图像图8-2注：matlab画图就是描点连线，因此如果点取得不密，画出来就成了折线图，请试验之 继续输入命令 hold on % 用此命令使刚才画的图像保持不动，下面继续在画好的图形上添加新的图形Y=sin(10*x);plot(x,y,r:,x,Y,b) % 同时画两个函数，其中的r:，b表示给图形设置颜色、线型等属性，见表8-3，若要改变颜色，在座标对后面加上相关字串即可： x=0:0.1:10; plot(x,sin(x),r)若要同时改变颜色及图线型态（Line style），也是在坐标对后面加上相关字串即可：plot(x,sin(x),r*)用axis(xmin,xmax,ymin,ymax)函数来调整坐标轴的范围，对上面已经画出的图形继续执行下面的命令，观察结果axis(0,6,-1.5,1.6)再继续，对图形加上各种注解与处理：（见上表）xlabel(x轴); % x轴注解 ylabel(y轴); % y轴注解 title(余弦函数); % 加标题 legend(y = cos(x); % 加图形注解 gtext(y = cos(x); % 加图形注解 ,用鼠标定位注解位置grid on; % 显示格线表8-3命 令含 义plot绘图函数的参数 plot建立向量或矩阵各队队向量的图形字元 颜色字元 图线型态loglogx、y轴都取对数标度建立图形y 黄色. 点semilogxx轴用于对数标度，y轴线性标度绘制图形k 黑色o 圆semilogyy轴用于对数标度，x轴线性标度绘制图形w 白色x Xtitle给图形加标题b 蓝色+ +xlabel给x轴加标记g 绿色* *ylabel给y轴加标记r 红色- 实线text在图形指定的位置上加文本字符串c 亮青色: 点线gtext在鼠标的位置上加文本字符串m 锰紫色-. 点虚线grid打开网格线- 虚线例2：画椭圆a=0:pi/50:2*pi;% 角度从0到2，每隔/50一个点X=cos(a)*3;% 参数方程Y=sin(a)*2;plot(X,Y);xlabel(x), ylabel(y);title(椭圆)图8-3例3：绘制函数 在时的曲线。x=0:0.1:1y=x.*exp(-x) % 注意点运算plot(x,y),xlabel(x),ylabel(y),title(y=x*exp(-x)图略例4：画出衰减振荡曲线与它的包络线及。t 的取值范围是0, 4 。t=0:pi/50:4*pi;y0=exp(-t/3);y=exp(-t/3).*sin(3*t); plot(t,y,-r,t,y0,:b,t,-y0,:b) % -r表示红色实线，:b表示蓝色点线，见表8-3grid图8-4例5：在同一个画面上建立几个坐标系, 用subplot(m,n,p)命令；把一个画面分成mn个图形区域, p代表当前的区域号，在每个区域中分别画一个图,下面在一个画面里同时画四个图像，注意哪个图像在哪个区域x=linspace(0,2*pi,30); y=sin(x); z=cos(x);u=2*sin(x).*cos(x); v=sin(x)./cos(x);subplot(2,2,1),plot(x,y),axis(0 2*pi -1 1),title(sin(x)subplot(2,2,2),plot(x,z),axis(0 2*pi -1 1),title(cos(x)subplot(2,2,3),plot(x,u),axis(0 2*pi -1 1),title(2sin(x)cos(x)subplot(2,2,4),plot(x,v),axis(0 2*pi -20 20),title(sin(x)/cos(x)图8-58.3.2 三维图形例6：三维螺旋线： t=0:pi/50:10*pi; plot3(sin(t),cos(t),t) %参数方程 grid %添加网格线图8-6例7：t=linspace(0,20*pi, 501); plot3(t.*sin(t), t.*cos(t), t); %注意点乘图略也可以同时画出两条曲线，格式与二维情况类似，这里不举例了，可以试着把上面两个图形画在一个坐标系里。例8：画空间中的曲面要用mesh命令，画出由函数形成的立体网状图: a=linspace(-2, 2, 25); % 在x轴上从（-2，2）取25点 b=linspace(-2, 2, 25); % 在y轴上取25点 x,y=meshgrid(a, b); % x和y都是2525的矩阵 z=x.*exp(-x.2-y.2); % 计算函数值，z也是2525的矩阵 mesh(x, y, z); % 画出立体网状图 图8-7例9：surf和mesh的用法类似： a=linspace(-2, 2, 25); % 在x轴上取25点 b=linspace(-2, 2, 25); % 在y轴上取25点 x,y=meshgrid(a, b); % x和y都是2525的矩阵 z= x.2-y.2; % 计算函数值，z也是2525的矩阵，马鞍面 surf(x, y, z); % 画出立体曲面图图8-88.4 程序设计MATLAB作为一种高级语言，它不仅可以如前几节所介绍的那样，以一种人机交互式的命令行的方式工作，还可以像BASIC、FORTRAN、C等其他高级计算机语言一样进行控制流的程序设计，即编制一种以 .m为扩展名的MATLAB程序,简称M文件,所谓M文件就是由MATLAB语言编写的可在MATLAB语言环境下运行程序源代码文件。由于商用的MATLAB软件是用C语言编写而成。因此，M文件的语法与C语言十分相似。对广大参加建模竞赛且学过C语言的同学来说，M文件的编写是相当容易的。8.4.1 M文件可以分为脚本文件(Script)和函数文件(Function)两种。1 脚本文件：（1）编写文档:点击MATLAB指令窗口上面最左端的图标 ，即新建文件，就可打开MATLAB文件编辑器，像word一样。用户即可在空白窗口中编写程序。例如输入下面的程序： x=linspace(0,2*pi,20); y=sin(x); plot(x,y,r+) title(正弦曲线) (2）点击文件编辑器上面工具条中的保存，命名（例如将上面的程序命名为picture），然后保存。这样的文件就是M-文件：picture.m (3)运行：i）在命令窗口中输入文件名(如上面的picture),然后执行。ii）或直接在文件编辑器上面的工具条中找到debug里面的run（即运行），点击即可（或直接按F5）。iii）如果发现错误，在写好的程序里直接改正，然后再保存，再运行，注意一定要先保存，后运行。2 函数文件：我们经常用到的像sin、cos、exp这样的一些函数都是MATLAB软件自身所带的函数，因此直接应用即可，但有时我们为了解决一些问题需要自己编写函数。自己编写函数有两个基本要求i）必须在MATLAB文件编辑器中编写，也是M-文件。ii）函数名和文件名必须相同。 例1: 编写函数 , 计算 （1）打开MATLAB文件编辑器，即点图标，输入 function y= fun1(x) % 表示y是x的函数，x是自变量，fun1是我们自己命名的函数名y=(x3 - 2*x2 + x - 6.3)/(x2 + 0.05*x - 3.14);然后保存。注：在自己编写的函数前都要写上function（关键字），表示这是自己定义的函数。fun1表示函数名，那么最后保存文件也要命名为fun1。 （2）这样在命令窗口中就可以像应用sin、cos那样来使用函数fun1，如：在命令窗口中输入 fun1(1)*fun1(2)+fun1(3)*fun1(3) ans = -12.6023例2：也可以多输入、多输出，如输入向量x，y，要求编写函数输出x、y各自的长度、x与y的数量积、向量积。(1)建立M文件如下function L1,L2,nei,wai=fun2(x,y); % 定义函数名为fun2L1=sqrt(sum(x.2), % x的长度L2=sqrt(sum(y.2), % y的长度nei=dot(x,y), % x与y的数量积wai=cross(x,y), % x与y的向量积（2）保存命名为fun2.m（3）在命令窗口中输入如下语句 x=4 3 1;y=-2 1 6;fun2(x,y)运行结果为：L1 =5.0990L2 =6.4031nei = 1wai =17 -26 108.4.2 程序结构和任何计算机语言一样，matlab有三种程序结构：顺序结构、选择结构（或分支结构）、循环结构1顺序结构，即按照程序语句的顺序逐条运行命令，前面函数作图都是顺序结构，这里不再叙述2选择结构例3：if-end语句，例：cost=10;number=12; % cost和number都是变量名if number8 % 如果number8执行后面的语句，否则跳出sums=number*0.95*cost;end,sums 例4：if-else-end语句，例：cost=10;number=5; % 改变number的初值，看结果有何不同if number8 sums=number*0.95*cost;else sums=number*0.5*cost;end,sums例5：定义函数，计算，建立函数文件fff.mfunction y=fff(x)if x=2 & x8 & x y=fff(0.1),fff(1),fff(9),fff(22),fff(2*pi)结果y =1.1000 2.0000 31.0000 -1.0088 18.84963循环结构for循环语句（这里的for语句与C语言中的for语句不同，要更简单一些）例6：一个简单的for循环示例for i=1:10; % i依次取1,2,10,.x(i)=2*i; % 对每个i值，重复执行该指令end; % 表示循环结束，每一个for要对应一个endx % 要求显示运行后数组x的值。输入后观察结果，体会for语句的作用。注：在MATLAB里（在C语言中也一样），“”的作用表示把等号右边的值送给左边的变量，这和数学中相等的含义不同。while循环语句例7: Fibonacci 数列：1，1，2，3，5，8， 即：，（ 1，2，3）现要求该数列中第一个大于10000 的元素。a(1)=1;a(2)=1;i=2;while a(i)=10000 % 当a(i)=10000a(i),break; % 表示跳出循环end;end, i8.5 插值与拟合作为Matlab的简单应用，本节我们讲插值与拟合。插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。8.5.1 一维插值在实际中遇到的函数，有许多我们并不知道其表达式，而只知道它在某些点上的函数值，如，但其它点上的值就不知道了。现在为了获得这些其它点上的函数值，我们设法构造一个比较简单的函数（例如多项式）来近似代替，并要求经过那些已知点，即与在这些点上的函数值相同。这样我们想要知道在某点的函数值就可以通过计算在该点的值来近似。这种方法就叫插值法，最简单的方法就是把各相邻已知点用线段连起来，叫做线性插值，工程上常用的还有二次插值、三次插值、三次样条插值等等，具体方法可参看任何一本数值方法的教材，使用这些插值方法不需编制函数程序，MATLAB自身提供了内部函数interp1供我们直接调用，格式如下：yi=interp1(x,y,xi,method) 对一组点(x,y) 进行插值，计算插值点xi的函数值。x为节点向量值，y为对应的节点函数值。 method用来指定插值的算法。默认为线性算法。其值常用的可以是如下的字符串： nearest 线性最近项插值。 linear 线性插值。 spline 三次样条插值。 cubic 三次插值。所有的插值方法要求x是单调的（从小到大），但不要求等距。例1：考虑下列问题，12小时内，一小时测量一次室外温度。数据如下：时间：1，2，3， 4， 5， 6， 7， 8， 9，10，11，12温度：5，8，9，15，25，29，31，30，22，25，27，24现在根据以上数据估计3.2，4.7等时刻的温度hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24;t=interp1(hours, temps, 3.2,4.7) % 一阶线性插值，如果只估计一个点的值，则无须加方括号运行结果t =10.2000 22.0000输入如下程序，画出插值曲线hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12; %每隔0.1一个点，共111个点t=interp1(hours, temps, h) ; %给出这111个点的插值的结果plot(hours, temps, + , h, t)图8-9工程上常用所谓三次样条插值，效果比较好，只需在上面插值命令后加上spline即可t=interp1(hours, temps, 3.2，4.7, spline);t=interp1(hours, temps, h，spline) ;图像如下图8-10Matlab也能够完成二维插值运算，函数为interp2，用法与interp1基本相同，只是输入和输出参数均为矩阵，对应于平面上的数值点。8.5.2 曲线拟合在科学实验的统计方法研究中,往往要从一组实验数据中寻找出自变量x 和因变量y之间的函数关系。由于观测数据往往不够准确，因此并不要求经过所有的点 （前述插值法要求插值函数必须经过所有已知点）,而只要求在给定点上误差按照某种标准达到最小，通常采用欧氏范数作为误差量度的标准。这就是所谓的最小二乘法。我们找到的这条曲线叫做所给数据点的拟合曲线，一般先在平面直角坐标系中描出所给数据点的图形，称为散点图，其次观察散点图，根据经验确定一条曲线，尽可能好的反映x与y的变化关系（通常用多项式），设出其一般形式（如含待定系数的多项式），然后根据已知数据，具体求出这条曲线。在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行（多项式拟合）。例2：两组数据如下：x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1;y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;先做出x与y的散点图（用命令plot(x,y,o)）图8-11估计变量x与y之间大概是什么样的函数关系，比如估计是一次函数，则用下面的命令：n=1;p1=polyfit(x,y,