MATLAB 实用教程：第5章MATLAB程序设计.doc

第5章 MATLAB程序设计5.1脚本文件和函数文件M文件有两种形式：M脚本文件和M函数文件。5.1.1 M文本编辑器MATLAB的M文件是通过M文件编辑调试器窗口(EditorDebugger)来创建的。图5.1 M文件编辑/调试器窗口单击MATLAB桌面上的图标，或者单击菜单“File”“New”“M-file”，可打开空白的M文件编辑器，也可以通过打开已有的M文件来打开M文件编辑器。如图5.1所示为打开已创建的M文件。5.1.2 M文件的基本格式下面介绍绘制二阶系统时域曲线的M文件，欠阻尼系统的时域输出y与x的关系为，【例5.1】为M脚本文件，【例5.2】为M函数文件。【例5.1】用M脚本文件绘制二阶系统时域曲线。%EX0501 二阶系统时域曲线%画阻尼系数为0.3的曲线x=0:0.1:20;y1=1-1/sqrt(1-0.32)*exp(-0.3*x).*sin(sqrt(1-0.32)*x+acos(0.3)plot(x,y1,r) 【例5.2】创建一个画二阶系统时域曲线的函数，阻尼系数zeta为函数的输入参数。function y=Ex0502(zeta)% EX0502 Step response of quadratic system.% 二阶系统时域响应曲线% zeta 阻尼系数% y 时域响应% copyright 2003-08-01x=0:0.1:20;y=1-1/sqrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta)plot(x,y)M函数文件的基本格式：函数声明行H1行(用%开头的注释行)在线帮助文本(用%开头)编写和修改记录(用%开头)函数体 例如，在命令窗口输入help和lookfor命令查看帮助信息： help Ex0502 EX0502 Step response of quadratic system. 二阶系统时域响应曲线 zeta 阻尼系数 y 时域响应 lookfor 二阶系统时域响应Ex0502.m: %二阶系统时域响应5.1.3 M脚本文件脚本文件的特点：(1) 脚本文件中的命令格式和前后位置，与在命令窗口中输入的没有任何区别。(2) MATLAB在运行脚本文件时，只是简单地按顺序从文件中读取一条条命令，送到MATLAB命令窗口中去执行。(3) 与在命令窗口中直接运行命令一样，脚本文件运行产生的变量都是驻留在MATLAB的工作空间(workspace)中，可以很方便地查看变量，除非用clear命令清除；脚本文件的命令也可以访问工作空间的所有数据，因此要注意避免变量的覆盖而造成程序出错。【例5.1续】在M文件编辑调试器窗口中编写M脚本文件绘制二阶系统的多条时域曲线。(1) 单击MATLAB桌面上的图标打开M文件编辑器。(2) 将命令全部写入M文件编辑器中，为了能标志该文件的名称，在第一行写入包含文件名的注释。保存文件为Ex0501.m。%EX0501 二阶系统时域曲线x=0:0.1:20;y1=1-1/sqrt(1-0.32)*exp(-0.3*x).*sin(sqrt(1-0.32)*x+acos(0.3)plot(x,y1,r) %画阻尼系数为0.3的曲线hold ony2=1-1/sqrt(1-0.7072)*exp(-0.707*x).*sin(sqrt(1-0.7072)*x+acos(0.707)plot(x,y2,g) %画阻尼系数为0.707的曲线y3=1-exp(-x).*(1+x)plot(x,y3,b) %画阻尼系数为1的曲线图5.2 运行界面(3) 选择M文件编辑器菜单“Debug”“Run”，就可以在图形窗中看到如图5.2所示的曲线。查看工作空间的变量： whos Name Size Bytes Class x 1x201 1608 double array y1 1x201 1608 double array y2 1x201 1608 double array y3 1x201 1608 double arrayGrand total is 804 elements using 6432 bytes 5.1.4 M函数文件函数文件的特点：(1) 第一行总是以“function”引导的函数声明行；函数声明行的格式： function 输出变量列表 = 函数名(输入变量列表)(2) 函数文件在运行过程中产生的变量都存放在函数本身的工作空间；(3) 当文件执行完最后一条命令或遇到“return”命令时，就结束函数文件的运行，同时函数工作空间的变量就被清除；(4) 函数的工作空间随具体的M函数文件调用而产生，随调用结束而删除，是独立的、临时的，在MATLAB运行过程中可以产生任意多个临时的函数空间。【例5.2续】在M文件编辑调试器窗口编写计算二阶系统时域响应的M函数文件，并在MATLAB命令窗口中调用该文件。创建M函数文件并调用的步骤如下：(1) 编写函数代码function y=Ex0502(zeta)%EX0502画二阶系统时域曲线x=0:0.1:20;y=1-1/sqrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta)plot(x,y)(2) 将函数文件保存为“Ex0502.m”。(3) 在MATLAB命令窗口输入以下命令，则会出现f的计算值和绘制的曲线： f=Ex0502(0.3)程序分析： 第一行指定该文件是函数文件，文件名为“Ex0502”，输入参数为阻尼系数zeta，输出参数为时域响应y。 当函数文件调用结束，查看x、y： x? Undefined function or variable x. y? Undefined function or variable y.注意：M脚本文件和M函数文件的文件名及函数名的命名规则与MATLAB变量的命名规则相同。5.2程序流程控制5.2.1 for . end循环结构语法： for 循环变量=array循环体end 说明：循环体被循环执行，执行的次数就是array的列数，array可以是向量也可以是矩阵，循环变量依次取array的各列，每取一次循环体执行一次。【例5.3】使用for . end循环的array向量编程求出 1+3+5.+100 的值。% EX0503 使用向量for循环sum=0;for n=1:2:100 sum=sum+n;endsum sum = 2500 计算的结果为：sum =2500。程序说明：循环变量为n，n对应为向量1:2:100，循环次数为向量的列数，每次循环n取一个元素。【例5.4】使用for . end循环的array矩阵编程将单位阵转换为列向量。% EX0504 使用矩阵for循环sum=zeros(6,1);for n=eye(6,6) sum=sum+n;endsum sum = 1 1 1 1 1 1 程序分析：循环变量n对应为矩阵eye(6,6)的每一列，即第一次n为1;0;0;0;0;0，第一次n为0;1;0;0;0;0；循环次数为矩阵的列数6。5.2.2 while . end循环结构语法： while 表达式循环体end 说明：只要表达式为逻辑真，就执行循环体；一旦表达式为假，就结束循环。表达式可以是向量也可以是矩阵，如果表达式为矩阵则当所有的元素都为真才执行循环体，如果表达式为nan，MATLAB认为是假，不执行循环体。【例5.5】与【例5.3】相同，计算1+3+5.+100 的值。% EX0505 使用while循环sum=0;n=1;while n=100 sum=sum+n; n=n+2 ;endsumn sum = 2500n = 101 程序分析：可以看出while . end循环的循环次数由表达式来决定，当n=101就停止循环。5.2.3 Ifelseend条件转移结构语法：if 条件式1语句段1elseif 条件式2语句段2. else语句段n+1end说明：当有多个条件时，条件式1为假再判断elseif的条件式2，如果所有的条件式都不满足，则执行else的语句段n+1，当条件式为真则执行相应的语句段；Ifelseend结构也可以是没有elseif和else的简单结构。【例5.6】用If结构执行二阶系统时域响应，根据阻尼系数0zeta0)&(zeta *Inner matrix dimensions must agree. 程序分析：试探出矩阵的大小不匹配时，矩阵无法相乘，则再执行catch后面的语句段，将a的子矩阵取出与b矩阵相乘。可以通过这种结构灵活地实现矩阵的乘法运算。5.2.6流程控制语句1. break命令break命令可以使包含break的最内层的for或while语句强制终止，立即跳出该结构，执行end后面的命令，break命令一般和If结构结合使用。【例5.9】将【例5.5】增加条件用If与break命令结合，停止while循环。计算1+3+5.+100 的值，当和大于1000时终止计算。% EX0509 用break终止while循环sum=0;n=1;while n=100 if sum ex0514Too many input arguments.也可以在工作空间查看函数体定义的输入参数个数： nargin(Ex0514)ans = 2【例5.14续】添加以下程序，查看用nargout变量获取输出参数个数。if nargout=0%当输出参数个数为0时，运算结果为0 sum=0;end在命令窗口调用Ex0514函数，当输出参数格式不同时，结果如下： Ex0514(2,3) %当输出参数个数为0时 ans = 0 y=Ex0514(2,3) %当输出参数个数为1时 y = 5 y,n,x=Ex0514%当输出参数个数为太多时 ? Error using = ex0514Too many output arguments. 程序分析：当输出参数个数为0时，即使有两个输入参数，运算结果也为0，结果送给ans变量；当输出的参数个数太多，也会出错。(2) varargin和varargout变量varargin和varargout可以获得输入输出变量的各元素内容。【例5.15】计算所有输入变量的和。function y,n=Ex0515(varargin)% EX0515 使用可变参数vararginif nargin=0%当没有输入变量时输出0 disp(No Input variables.) y=0;elseif nargin=1 %当一个输入变量时，输出该数 y=varargin1;else n=nargin; y=0; for m=1:n y=vararginm+y;%当有多个输入变量时，取输入变量循环相加 endendn=nargin; 在MATLAB的命令窗口中输入不同个数的变量调用函数Ex0515，结果如下： y,n=Ex0515(1,2,3,4)%输入4个参数 y = 10n = 4 y,n=Ex0515(1) %输入1个参数 y = 1n = 1 y,n=Ex0515%无输入参数 No Input variables.y = 0n = 0 程序分析：n为输入参数的个数；y为求和运算的结果。5.3.4程序举例【例5.16】根据阻尼系数绘制不同二阶系统的时域响应，当欠阻尼时，当临界阻尼时，当过阻尼时。M文件的程序代码如下：function y=Ex0516(z1)% EX0516 主函数调用子函数，根据阻尼系数绘制二阶系统时域曲线t=0:0.1:20;if (z1=0)&(z11) y=plotxy1(z1,t);elseif z1=1 y=plotxy2(z1,t);else y=plotxy3(z1,t);endfunction y1=plotxy1(zeta,x)%画欠阻尼二阶系统时域曲线y1=1-1/sqrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta);plot(x,y1)function y2=plotxy2(zeta,x)%画临界阻尼二阶系统时域曲线y2=1-exp(-x).*(1+x);plot(x,y2)function y3=plotxy3(zeta,x)%画过阻尼二阶系统时域曲线y3=1-1/(2*sqrt(zeta2-1)*(exp(-(zeta-sqrt(zeta2-1)*x)./(zeta-sqrt(zeta2-1). -exp(-(zeta+sqrt(zeta2-1)*x)./(zeta+sqrt(zeta2-1);plot(x,y3)程序分析：主函数名为Ex0516，三个子函数名分别为plotxy1、plotxy2、plotxy3，文件保存为Ex0516.m。当在命令窗口中输入以下命令： y=Ex0516(0.3); hold on y=Ex0516(0.707); y=Ex0516(1); y=Ex0516(2); 则产生如图5.4所示时域响应曲线。图5.4 不同阻尼系数的时域响应曲线【例5.17】编写M函数文件，通过流程控制语句，建立如下的矩阵。 function y=Ex0517(m)% EX0517 用循环流程控制语句创建矩阵y=0;m=m-1;for n=1:m y=0,y;%创建全0行endfor n=1:m a=1:1:n; b=a; for k=m:-1:n b=0,b; end y=b;y; n=n+1; end程序分析：将矩阵的行列数用输入参数m来确定，输出参数为矩阵y。使用双重循环来创建矩阵，将文件保存为Ex0517.m。在命令窗口中调用Ex0517函数： y=Ex0517(5) y = 0 1 2 3 4 0 0 1 2 3 0 0 0 1 2 0 0 0 0 1 0 0 0 0 0 5.4 M文件性能的优化和加速5.4.1 P码文件1. P码文件的生成P码文件使用pcode命令生成，生成的P码文件与原M文件名相同，其扩展名为“.p”。语法：pcode Filename.m %在当前目录生成Filename.ppcode Filename.m -inplace %在Filename.m所在目录生成Filename.p pcode Ex0517.m 则在当前目录就生成了P码文件Ex0517.p。2. P码文件的特点(1) P码文件的运行速度比原M文件速度快(2) 存在同名的M文件和P码文件时则P码文件被调用(3) P码文件保密性好用字处理软件打开Ex0517.p文件，看到的是乱码。5.4.2 M文件性能优化1. 使用循环时提高速度的措施循环语句及循环体是MATLAB编程的瓶颈问题，改进这种状况有三种方法：(1) 尽量用向量的运算来代替循环操作。(2) 在必须使用多重循环的情况下，如果两个循环执行的次数不同，则建议在循环的外环执行循环次数少的，内环执行循环次数多的，也可以显著提高速度。(3) 应用Mex技术如果耗时的循环不可避免，就应该考虑用其他语言，如C或Fortran语言，按照Mex技术要求的格式编写相应部分的程序，然后通过编译联接，形成在MATLAB可以直接调用的动态链接库(DLL)文件，这样就可以显著地加快运算速度(在8.1.1小节介绍)。 2. 大型矩阵的预先定维给大型矩阵动态地定维是个很费时间的事。【例5.18】将【例5.17】中的双重循环改为单循环，并先用zeros函数定维来提高运行速度，创建矩阵。function y=Ex0518(m)% EX0518 先定维再创建矩阵m=m-1;y=zeros(m);for n=1:m-1 a=1:m-n; y(n,n+1:m)=a;endy在命令窗口中分别调用Ex0517和Ex0518函数，比较其运行速度的不同： y=Ex0517(200) y=Ex0518(200)3. 优先考虑内在函数矩阵运算应该尽量采用MATLAB的内在函数，因为内在函数是由更底层的C语言构造的，其执行速度显然很快。 4. 采用高效的算法在实际应用中，解决同样的数学问题经常有各种各样的算法。因此，应寻求更高效的算法。5. 尽量使用M函数文件代替M脚本文件由于M脚本文件每次运行时，都必须把程序装入内存，然后逐句解释执行，十分费时。5.4.3 JIT和加速器1. JIT和加速器的加速范围JIT和加速器的应用范围如下： 只对维数不超过3的“非稀疏”数组起作用； 只对“双精度”、“整数”、“字符串”和“逻辑”等四种数据类型起作用； 只对内部函数的调用起作用，对用户的M函数或MEX文件不起作用； 只对控制语句for、while、if、elseif、switch中标量运算的条件表达式起作用； 当一行程序中含有“不可加速的”命令或变量时，整行都不被加速； 当变量改变数据类型或维数，则该语句不被加速； 如果i和j不以虚数单位形式使用，则该语句不被加速； 在Intel系列CPU硬件上，Windows和Linux系统加速能力最强。2. 使用程序性能剖析窗口 (1) 打开程序性能剖析窗口指示灯命令输入栏时间指示器开始剖析图5.5 程序性能剖析窗口选择菜单“View”“Profiler”；或使用在命令窗口输入“profile viewer”命令都可以打开程序性能剖析窗口，如图5.5所示。(2) 对MATLAB的命令进行剖析对MATLAB的命令进行剖析有两种方式： 在上图中的“命令输入栏”中输入需要剖析的命令，然后单击“Start Profiling”按钮，则指示灯变为绿色，“Start Profiling”按钮变灰色，时间指示器的计时在累加；当命令运行结束时，状态恢复，并出现剖析报告。 将“命令输入栏”清空，然后单击“Start Profiling”按钮开始剖析，则“Start Profiling”按钮变为“Stop Profiling”，指示灯变为绿色表示启动剖析，时间指示器的计时在累加；然后在MATLAB命令窗口中输入需要剖析的命令，进行剖析；当命令运行结束，单击“Stop Profiling”按钮停止剖析，则状态恢复，并出现