实用MATLAB教程第五章

第5章 MATLAB程序设计

MATLAB除了可以在命令窗口中编写行命令外，作为一种高级应用软件还可以生成自己的程序文件，要充分发挥MATLAB的功能，必须掌握MATLAB的程序设计。

5.1脚本文件和函数文件

MATLAB程序代码所编写的文件通常以“.m”为扩展名，因此这些文件称为M文件。M文件是一个ASCⅡ码文件，可以用任何字处理软件来编写，MATLAB的程序在第一次运行时由于逐句解释运行程序，故速度比编译型的慢，但M文件一经运行就将编译代码放在内存中，再次运行的速度就大大加快了。

M文件有两种形式：M脚本文件和M函数文件，M 函数文件是MATLAB程序设计的主流。MATLAB本身的一系列工具箱的各种内部函数就是M函数文件，用户可以为某种目的专门编写一组MATLAB函数文件组成工具箱。

5.1.1 M文本编辑器

MATLAB的M文件是通过M文件编辑／调试器窗口(Editor／Debugger)来创建的，M文件编辑／调试器窗口是集编辑与调试功能于一体的环境。

单击MATLAB桌面上的

图标，或者单击菜单“File”——“New”——“M-file”，

可打开空白的M文件编辑器，也可以通过打开已有的M文件来打开M文件编辑器。如图5.1所示为打开已创建的M文件。

图5.1 M文件编辑／调试器窗口

M文件编辑器窗口会以不同的颜色显示注释、关键词、字符串和一般程序代码；可以方便地打开和保存M文件并进行编辑和调试；编辑功能有大多数编辑器都有的复制、粘贴、查找等，还有设书签、定位、清除工作空间和命令窗口、加注释、缩进等功能。

5.1.2 M文件的基本格式

M脚本文件和M函数文件的格式不同，下面介绍绘制二阶系统时域曲线的M文件，欠阻尼系统的时域输出y与x的关系为y?1?11??2【例5.1】为Me??xsin(1??2x?acos?)，

脚本文件，【例5.2】为M函数文件。

【例5.1】用M脚本文件绘制二阶系统时域曲线。

%EX0501 二阶系统时域曲线 %画阻尼系数为0.3的曲线 x=0:0.1:20;

y1=1-1/sqrt(1-0.3^2)*exp(-0.3*x).*sin(sqrt(1-0.3^2)*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-01

x=0:0.1:20;

y=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt(1-zeta^2)*x+acos(zeta)) plot(x,y)

M函数文件的基本格式： 函数声明行

H1行(用%开头的注释行) 在线帮助文本(用%开头) 编写和修改记录(用%开头) 函数体 说明：

(1) 函数声明行是M函数文件必须有的，M脚本文件没有；函数名和文件名一致，当不一致时，MATLAB以文件名为准，Ex0502函数保存为Ex0502.m文件；

(2) H1行通常包含大写的函数文件名，可以提供给help和lookfor关键词查询使用； (3) 在线帮助文本通常包含函数输入输出变量的含义、格式说明； (4) 编写和修改记录一般在空一行后，记录作者、日期和版本记录，用于软件档案管理。 例如，在命令窗口输入help和lookfor命令查看帮助信息：

>> help Ex0502

EX0502 Step response of quadratic system. 二阶系统时域响应曲线 zeta 阻尼系数 y 时域响应

>> lookfor '二阶系统时域响应' Ex0502.m: %二阶系统时域响应

help命令显示M文件的第一个连续注释块，lookfor命令显示第一行注释，因此应在第一行的注释中多包含一些该函数特征的内容。lookfor命令查找的必须是MATLAB搜索路径上的文件。

在后面的例题中，为了节约篇幅，只保留H1注释行。

(5) 函数体由实现M函数文件功能的MATLAB命令组成。

5.1.3 M脚本文件

MATLAB的脚本文件(Script File)比较简单，当用户需要在命令窗口运行大量的命令时，直接从命令窗口输入比较繁琐，就可以将这一组命令存放在脚本文件中，运行时只要输入脚本文件名，MATLAB就会自动执行该文件的命令。

脚本文件的特点：

(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.3^2)*exp(-0.3*x).*sin(sqrt(1-0.3^2)*x+acos(0.3)) plot(x,y1,'r') %画阻尼系数为0.3的曲线 hold on

y2=1-1/sqrt(1-0.707^2)*exp(-0.707*x).*sin(sqrt(1-0.707^2)*x+acos(0.707)) plot(x,y2,'g') %画阻尼系数为0.707的曲线 y3=1-exp(-x).*(1+x) plot(x,y3,'b') %画阻尼系数为1的曲线

(3) 选择M文件编辑器菜单“Debug”——“Run”，就可以在图形窗中看到如图5.2所示的曲线。

查看工作空间的变量： >> whos Name Size Bytes Class 图5.2 运行界面

x 1x201 1608 double array y1 1x201 1608 double array y2 1x201 1608 double array y3 1x201 1608 double array Grand total is 804 elements using 6432 bytes

可以看出M脚本文件的变量都存放在MATLAB工作空间中。

5.1.4 M函数文件

MATLAB函数文件(Function File)的第一行必须包含“function”，每一个函数文件都要定义一个函数。函数文件可以接受输入变量，并将运算结果送到输出变量，从外面看函数文件的功能就是，将数据送到函数文件处理后再将结果送出来，易于维护和修改。因此，函数文件可适用于大型程序代码的模块化。

函数文件的特点：

(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-zeta^2)*exp(-zeta*x).*sin(sqrt(1-zeta^2)*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'.

可以看出x和y都未定义，说明都已被清除。

注意：M脚本文件和M函数文件的文件名及函数名的命名规则与MATLAB变量的命名规则相同。

5.2程序流程控制

与大多数计算机语言一样，MATLAB支持各种流程控制结构，如顺序结构、循环结构和条件结构(又称为分支结构)，另外MATLAB还支持一种新的结构——试探结构。

5.2.1 for ... end循环结构

MATLAB的循环结构有两种：for ... end结构和while ... 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 end

计算的结果为：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 end

计算的结果为： 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循环结构

for ... end循环的循环次数确定，而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 end

计算的结果为： sum =

2500 n = 101

程序分析：可以看出while ... end循环的循环次数由表达式来决定，当n=101就停止循环。

5.2.3 If…else…end条件转移结构

If…else…end结构是最常见的条件转移结构。 语法：

if 条件式1 语句段1 elseif 条件式2 语句段2 ... else 语句段n+1 end

说明：当有多个条件时，条件式1为假再判断elseif的条件式2，如果所有的条件式都不满足，则执行else的语句段n+1，当条件式为真则执行相应的语句段；If…else…end结构也可以是没有elseif和else的简单结构。

【例5.6】用If结构执行二阶系统时域响应，根据阻尼系数0

function y=Ex0506(zeta)

% EX0506 使用if结构的二阶系统时域响应 x=0:0.1:20;

if (zeta>0)&(zeta<1)

y=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt(1-zeta^2)*x+acos(zeta)); elseif zeta==1

y=1-exp(-x).*(1+x); end

plot(x,y)

5.2.4 switch…case开关结构

switch…case结构是有多个分支结构的条件转移结构。 语法：

switch 开关表达式 case 表达式1 语句段1 case表达式2 语句段2 ...

otherwise 语句段n end 说明：

(1) 将开关表达式依次与case后面的表达式进行比较，如果表达式1不满足，则与下一个表达式2比较，如果都不满足则执行otherwise后面的语句段n；一旦开关表达式与某个表达式相等，则执行其后面的语句段。

(2) 开关表达式只能是标量或字符串。

(3) case后面的表达式可以是标量、字符串或元胞数组，如果是元胞数组则将开关表达式与元胞数组的所有元素进行比较，只要某个元素与开关表达式相等，就执行其后的语句段。

【例5.7】用switch…case开关结构得出各月份的季节。

% EX0507 使用switch结构 for month=1:12; switch month case{3,4,5}

season='spring' case{6,7,8}

season='summer' case{9,10,11}

season='autumn' otherwise

season='winter' end end

程序分析：开关表达式为向量1:12，case后面的表达式为元胞数组，当元胞数组的某个元素与开关表达式相等，就执行其后的语句段。

5.2.5 try... catch... end试探结构

MATLAB 从 5.3 版本开始提供了一种新的试探结构try... catch... end，这种语句结构是 其他语言如C 语言等中所没有的。

语法： try 语句段1 catch 语句段2 end

说明：首先试探性地执行语句段1，如果在此段语句执行过程中出现错误，则将错误信息赋给保留的lasterr变量，并放弃这段语句，转而执行语句段2中的语句，当执行语句段2又出现错误，则终止该结构。

当试探结构运行结束后，可以调用lasterr函数查询出错原因，空字符串表示成功执行。 【例5.8】用try... catch... end结构来进行矩阵相乘运算。

% EX0508 try结构 n=4;

a=magic(n); m=3; b=eye(3); try

c=a*b catch

c=a(1:m,1:m)*b end lasterr

计算结果是： c =

16 2 3 5 11 10 9 7 6

用lasterr函数查看出错原因显示为：

ans =

Error using ==> *

Inner matrix dimensions must agree.

程序分析：试探出矩阵的大小不匹配时，矩阵无法相乘，则再执行catch后面的语句段，将a的子矩阵取出与b矩阵相乘。可以通过这种结构灵活地实现矩阵的乘法运算。

上述每种循环结构、分支结构和试探结构都可以相互嵌套。

5.2.6流程控制语句

在程序执行中，有一些可以控制程序流程的命令，下面主要介绍break、continue、return、pause、keyboard和imput命令。

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

sum=sum+n n=n+2 else

break end end

计算结果为： sum =

1024 n= 65

程序分析：while…end循环结构嵌套If…else…end分支结构，当sum为1024时跳出while循环结构，终止循环。

2. continue命令

continue命令用于结束本次for或while循环，与break命令不同的是continue只结束本次循环而继续进行下次循环。

【例5.10】将If命令与continue命令结合，计算的1～100中所有素数的和，判断是否为素数是将100以内的每个数都被2～n整除，不能被整除的就是素数。

% EX0510 用continue终止while循环 sum=2;ss=0; for n=3:100

for m=2:fix(sqrt(n)) if mod(n,m)==0 ss=1; %能被整除就用ss为1表示 break; %能被整除就跳出内循环 else

ss=0; %不能被整除就用ss为0表示 end end if ss==1

continue; %能被整除就跳出本次循环 end

sum=sum+n; end sum

计算结果是：

sum =

1060

程序分析：fix(sqrt(n))是将n取整；本程序为双重循环，两个for循环嵌套还嵌套一个if结构；当mod(n,m)==0时就用break跳出判断是否为素数的内循环，并继续用continue跳出求素数和的外循环而继续下次外循环。

3. return命令

return命令是终止当前命令的执行，并且立即返回到上一级调用函数或等待键盘输入命令，可以用来提前结束程序的运行。

return命令也可以用来终止键盘方式。

注意：当程序进入死循环，则按Ctrl+break键来终止程序的运行。

4. pause命令

pause命令用来使程序运行暂停，等待用户按任意键继续。用于程序调试或查看中间结果，也可以用来控制动画执行的速度。

语法：

pause %暂停 pause(n) %暂停n秒

5. keyboard命令

keyboard命令用来使程序暂停运行，等待键盘命令，执行完自己的工作后，输入return语句，程序就继续运行。keyboard命令可以用来在程序调试或程序执行时修改变量。

6. input命令

input命令用来提示用户应该从键盘输入数值、字符串和表达式，并接受该输入。 Input命令常用的几种格式如下：

>> a=input('input a number:') input a number:45 a = 45

>> b=input('input a number:','s') input a number:45 b = 45

>> input('input a number:') input a number:2+3 ans = 5

%输入数值给a

%输入字符串给b

%将输入值进行运算

5.3函数调用和参数传递

在MATLAB中，使用函数可以把一个比较大的任务分解为多个较小的任务，使程序模块化，每个函数完成特定的功能，通过函数的调用完成整个程序。

5.3.1子函数和私有函数

1. 子函数

在一个M函数文件中，可以包含一个以上的函数，其中只有一个是主函数，其它则为子函数。

(1) 在一个M文件中，主函数必须出现在最上方，其后是子函数，子函数的次序无任何限制；

(2) 子函数不能被其它文件的函数调用，只能被同一文件中的函数(可以是主函数或子函数)调用；

(3) 同一文件的主函数和子函数变量的工作空间相互独立； (4) 用help和lookfor命令不能提供子函数的帮助信息。

【例5.11】将【例5.2】画二阶系统时域曲线的函数作为子函数，编写画多条曲线的程序。

function Ex0511()

% EX0511 使用函数调用绘制二阶系统时域响应 z1=0.3;

Ex0502(z1); %调用Ex0502 hold on z1=0.5

Ex0502(z1) %调用Ex0502 z1=0.707;

Ex0502(z1) %调用Ex0502

function y=Ex0502(zeta)

%子函数，画二阶系统时域曲线 x=0:0.1:20;

y=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt(1-zeta^2)*x+acos(zeta)) plot(x,y)

程序分析：主函数是Ex0511，子函数是Ex0502，在主函数中三次调用子函数。程序保存为Ex0511.m文件。

2. 私有函数

私有函数是指存放在private子目录中的M函数文件，具有以下性质：

(1) 在private目录下的私有函数，只能被其父目录的M函数文件所调用，而不能被其它目录的函数调用，对其它目录的文件私有函数是不可见的，私有函数可以和其它目录下的函数重名；

(2) 私有函数父目录的M脚本文件也不可调用私有函数；

(3) 在函数调用搜索时，私有函数优先于其它MATLAB路径上的函数。

根据私有函数的特点，可在自己的工作目录下建立一个private子目录，但不要添加到MATLAB的搜索路径中。

3. 调用函数的搜索顺序

在MATLAB中调用一个函数，搜索的顺序如下： ? 查找是否子函数； ? 查找是否私有函数；

? 从当前路径中搜索此函数； ? 从搜索路径中搜索此函数。

5.3.2局部变量和全局变量

根据变量的作用域不同，可以将MATLAB程序中的变量分为局部变量和全局变量。

1. 局部变量

局部变量(Local Variables)是在函数体内部使用的变量，其影响范围只能在本函数内；每个函数在运行时，都占用独立的函数工作空间，此工作空间和MATLAB的工作空间是相互独立的，局部变量仅存在于函数的工作空间内。局部变量只在函数执行期间存在，当函数执行完变量就消失。

2. 全局变量

全局变量(Global Variables)是可以在不同的函数工作空间和MATALB工作空间中共享使用的变量。

全局变量在使用前必须用global定义，而且每个要共享全局变量的函数和工作空间，都必须逐个用global对变量加以定义。

【例5.12】修改【例5.11】在主函数和子函数中使用全局变量。

function Ex0512()

% EX0512 使用全局变量绘制二阶系统时域响应 global X X=0:0.1:20; z1=0.3;

Ex0502(z1); hold on z1=0.5;

Ex0502(z1); z1=0.707; Ex0502(z1);

function Ex0502(zeta)

%子函数，画二阶系统时域曲线 global X

y=1-1/sqrt(1-zeta^2)*exp(-zeta*X).*sin(sqrt(1-zeta^2)*X+acos(zeta)); plot(X,y);

程序分析：X变量为全局变量，在需要使用的主函数和子函数中都需要用global定义；同样如果在工作空间中定义X为全局变量后也可以使用：

>> global X who

Your variables are: X

注意：由于全局变量在任何定义过的函数中都可以修改，因此不提倡使用全局变量；必须使用时应十分小心，建议把全局变量的定义放在函数体的开始，全局变量用大写字符命名。

5.3.3函数的参数

MATLAB的函数调用过程实际上就是参数传递的过程。 函数调用的格式：

[输出参数1，输出参数2，…]=函数名(输入参数1，输入参数2，…)

1. 参数传递规则

在MATLAB中，函数具有自己的工作空间，函数内变量与外界(包括其它函数和工作空间)的唯一联系就是通过函数的输入输出参数。输入参数在函数中的任何变化，都仅在函数内进行，不会传递回去。

【例5.13】将【例5.11】画二阶系统时域的函数修改，使用输入输出参数来实现参数传递，如图5.3所示。

function Ex0513() % EX0513 参数传递绘制二阶系统时域响应 z1=0.3; [x1,y1]=Ex0502(z1); plot(x1,y1) hold on z1=0.5; [x2,y2]=Ex0502(z1); plot(x2,y2) z1=0.707; [x3,y3]=Ex0502(z1); function [x,y]=Ex0502(zeta) %子函数，画二阶系统时域曲线 x=0:0.1:20; y=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin... (sqrt(1-zeta^2)*x+acos(zeta)); 图5.3 参数传递 程序分析：主函数Ex0513调用子函数Ex0502，子函数中的zeta为输入参数，函数调用时将z1传递给子函数zeta，子函数计算后将输出参数x和y传回给主函数的x1、y1；主函数调用子函数三次，后面两次参数的传递也是同样。

2. 函数参数的个数

MATLAB函数的调用有一个与其他语言不同的特点是：函数的输入输出参数的数目都可以变化，用户可以根据参数的个数来编程。

(1) nargin和nargout变量

在MATLAB中有两个特殊变量(2.1.3小节介绍过)nargin和nargout，函数的输入输出参数的个数可以通过变量nargin和nargout获得，nargin用于获得输入参数的个数，nargout用于获得输出参数的个数。

语法：

nargin %在函数体内获取实际输入变量的个数 nargout %在函数体内获取实际输出变量的个数 nargin(’fun’) %在函数体外获取定义的输入参数个数 nargout(‘fun’) %在函数体外获取定义的输出参数个数

【例5.14】计算两个数的和，根据输入的参数个数不同使用不同的运算表达式。

function [sum,n]=Ex0514(x,y)

% EX0514 参数个数可变，计算x和y的和 if nargin==1 sum=x+0; %输入一个参数就计算与0的和 elseif nargin==0 sum=0; %无输入参数就输出0 else

sum=x+y; %输入的是两个数则就计算和 end

在命令窗口调用Ex0514函数，分别使用2个、1个和无输入参数结果如下：

>> [y,n]=Ex0514(2,3) y = 5 n = 2

>> [y,n]=Ex0514(2) y = 2 n = 1

>> [y,n]=Ex0514 y = 0 n = 0

注意：如果输入的参数多于输入参数个数，则会出错。

>> [y,n]=Ex0514(1,2,3) ??? Error using ==> ex0514 Too many input arguments.

也可以在工作空间查看函数体定义的输入参数个数：

>> nargin('Ex0514') ans = 2

【例5.14续】添加以下程序，查看用nargout变量获取输出参数个数。

if nargout==0 sum=0; end

%当输出参数个数为0时，运算结果为0

在命令窗口调用Ex0514函数，当输出参数格式不同时，结果如下：

>> Ex0514(2,3) %当输出参数个数为0时 ans = 0

>> y=Ex0514(2,3) %当输出参数个数为1时 y = 5

>> [y,n,x]=Ex0514 %当输出参数个数为太多时 ??? Error using ==> ex0514 Too many output arguments

程序分析：当输出参数个数为0时，即使有两个输入参数，运算结果也为0，结果送给ans变量；当输出的参数个数太多，也会出错。

(2) varargin和varargout变量

MATLAB还有两个特殊变量“varargin”和“varargout”可以获得输入输出变量的各元素内容，varargin和varargout都是元胞数组。

【例5.15】计算所有输入变量的和。

function [y,n]=Ex0515(varargin)

% EX0515 使用可变参数varargin if nargin==0 %当没有输入变量时输出0 disp('No Input variables.') y=0;

elseif nargin==1 %当一个输入变量时，输出该数 y=varargin{1}; else

n=nargin; y=0;

for m=1:n

y=varargin{m}+y; %当有多个输入变量时，取输入变量循环相加 end end

n=nargin;

在MATLAB的命令窗口中输入不同个数的变量调用函数Ex0515，结果如下：

>> [y,n]=Ex0515(1,2,3,4) y = 10 n = 4

>> [y,n]=Ex0515(1) y = 1 n = 1

>> [y,n]=Ex0515 No Input variables. y = 0 n = 0

%输入4个参数

%输入1个参数

%无输入参数

程序分析：n为输入参数的个数；y为求和运算的结果。

5.3.4程序举例

利用M文件的流程控制和函数调用，实现各种复杂的程序设计。

【例5.16】根据阻尼系数绘制不同二阶系统的时域响应，当欠阻尼时

y?1?11??21e??xsin(1??2x?acos?)，当临界阻尼时y?1?(1?x)e?x，当过阻尼时

2??(???2?1)xe?(????1)x?ey?1???22?2??1???1????1??2??。 ???

M文件的程序代码如下：

function y=Ex0516(z1)

% EX0516 主函数调用子函数，根据阻尼系数绘制二阶系统时域曲线 t=0:0.1:20;

if (z1>=0)&(z1<1) y=plotxy1(z1,t); elseif z1==1

y=plotxy2(z1,t); else

y=plotxy3(z1,t); end

function y1=plotxy1(zeta,x) %画欠阻尼二阶系统时域曲线

y1=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt(1-zeta^2)*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(zeta^2-1))*(exp(-((zeta-sqrt(zeta^2-1))*x))./(zeta-sqrt(zeta^2-1))... -exp(-((zeta+sqrt(zeta^2-1))*x))./(zeta+sqrt(zeta^2-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函数文件，通过流程控制语句，建立如下的矩阵。 ?0?0? y??0?????0100?0210?0321?0?n??n?1???n?2?

????0??function y=Ex0517(m)

% EX0517 用循环流程控制语句创建矩阵 y=0; m=m-1; for n=1:m y=[0,y]; %创建全0行 end

for 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码文件

一般的M文件都是文本文件，因此MATLAB源代码都能看到，如果希望程序的保密性好而且运行速度快，可以将M脚本文件和M函数文件转换成P码(Psedocode)文件。

1. P码文件的生成

P码文件使用pcode命令生成，生成的P码文件与原M文件名相同，其扩展名为“.p”。 语法：

pcode Filename.m %在当前目录生成Filename.p

pcode Filename.m -inplace %在Filename.m所在目录生成Filename.p

例如，将【例5.17】生成P码文件，在命令窗口中输入代码将Ex0517.m文件转换为P码文件：

>> pcode Ex0517.m

则在当前目录就生成了P码文件Ex0517.p。

2. P码文件的特点

(1) P码文件的运行速度比原M文件速度快

一个M文件第一次被调用时，MATLAB就将其进行语法分析，并生成P码文件(又称伪代码)存放在内存中。以后再次调用该M文件时，MATLAB就直接调用内存中的P码文件，而不需要再对原M文件重新语法分析，因此当再次调用M文件时运行的速度就高于第一次运行。

在命令窗口中两次调用Ex0517.m文件时，比较执行的时间：

>> y=Ex0517(100) >> y=Ex0517(100)

(2) 存在同名的M文件和P码文件时则P码文件被调用

当在MATLAB的同一目录或搜索路径中，同时存在同名的M文件和P码文件，则该文件名被调用时，执行的是P码文件。

为了测试，可以将Ex0517.m文件和Ex0517.p文件放在同一目录，将Ex0517.m文件修改并保存，然后在命令窗口运行该文件，结果仍然是与原来一样：

>> y=Ex0517(5)

(3) P码文件保密性好

用字处理软件打开Ex0517.p文件，看到的是乱码。

5.4.2 M文件性能优化

MATLAB语言是解释性语言，所以有时MATLAB程序的执行速度不是很理想，为了能够提高程序的运行效率，编程时应注意以下几个方面：

1. 使用循环时提高速度的措施

循环语句及循环体是MATLAB编程的瓶颈问题，MATLAB与其它编程语言不同，MATLAB的基本数据是向量和矩阵，所以编程时应尽量对向量和矩阵编程，而不要对矩阵的元素编程。改进这种状况有三种方法：

(1) 尽量用向量的运算来代替循环操作。

(2) 在必须使用多重循环的情况下，如果两个循环执行的次数不同，则建议在循环的外环执行循环次数少的，内环执行循环次数多的，也可以显著提高速度。

(3) 应用Mex技术

如果耗时的循环不可避免，就应该考虑用其他语言，如C或Fortran语言，按照Mex技术要求的格式编写相应部分的程序，然后通过编译联接，形成在MATLAB可以直接调用的动态链接库(DLL)文件，这样就可以显著地加快运算速度(在8.1.1小节介绍)。

2. 大型矩阵的预先定维

给大型矩阵动态地定维是个很费时间的事，由于MATLAB变量的使用之前不需要定义和指定维数，当变量新赋值的元素下标超出数组的维数时，MATLAB就为该数组扩维一次，大大的降低了运行的效率。

建议在定义大矩阵时，首先用MATLAB的内在函数，如zeros() 或 ones() 对其先进行定维，然后再进行赋值处理，这样会显著减少所需的时间的。

【例5.18】将【例5.17】中的双重循环改为单循环，并先用zeros函数定维来提高运行?0?0?速度，创建矩阵y??0?????0100?0210?0321?0?n??n?1???n?2?。

????0??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; end y

在命令窗口中分别调用Ex0517和Ex0518函数，比较其运行速度的不同：

>> y=Ex0517(200) >> y=Ex0518(200)

3. 优先考虑内在函数

矩阵运算应该尽量采用MATLAB的内在函数，因为内在函数是由更底层的C语言构造的，其执行速度显然很快。

4. 采用高效的算法 在实际应用中，解决同样的数学问题经常有各种各样的算法。例如求解定积分的数值解法在MATLAB中就提供了两个函数：quad和quad8，其中后一种算法在精度、速度上都明显高于前一种方法。因此，应寻求更高效的算法。

5. 尽量使用M函数文件代替M脚本文件

由于M脚本文件每次运行时，都必须把程序装入内存，然后逐句解释执行，十分费时。

5.4.3 JIT和加速器

由于MATLAB是以解释方式工作，其最主要的缺点就是运行速度慢，从MATLAB6.5版开始引入了JIT(Just In Time)和加速器(Accelerator)，从而使一些MATLAB程序的运行速度可以提高到C/Fortran程序的速度。

1. JIT和加速器的加速范围

并不是所有的MATLAB命令都可以通过JIT和加速器实现加速，JIT和加速器的应用范围如下：

? 只对维数不超过3的“非稀疏”数组起作用； ? 只对“双精度”、“整数”、“字符串”和“逻辑”等四种数据类型起作用； ? 只对内部函数的调用起作用，对用户的M函数或MEX文件不起作用；

? 只对控制语句for、while、if、elseif、switch中标量运算的条件表达式起作用； ? 当一行程序中含有“不可加速的”命令或变量时，整行都不被加速； ? 当变量改变数据类型或维数，则该语句不被加速；

? 如果i和j不以虚数单位形式使用，则该语句不被加速；

? 在Intel系列CPU硬件上，Windows和Linux系统加速能力最强。

2. 使用程序性能剖析窗口

使用程序性能剖析窗口可以对程序中命令的运行时间进行分析，找出运行的“瓶颈”，并对程序的“瓶颈”进行改写，直到剖析报告表明，绝大多数运行时间花费在MATLAB的内部函数上，才是比较优良的文件性能。

(1) 打开程序性能剖析窗口

选择菜单“View”——“Profiler”；或使用在命令窗口输入“profile viewer”命令都可以打开程序性能剖析窗口，如图5.5所示。

命令输入栏 开始剖析 时间指示器 指示灯 图5.5 程序性能剖析窗口 (2) 对MATLAB的命令进行剖析

对MATLAB的命令进行剖析有两种方式：

[输出变量列表]=函数名(‘funname’,输入变量列表)

说明：h_fun是要被执行的M函数文件的句柄，或者是内联函数和字符串；‘funname’是M函数文件名。

泛函命令有很多，下面介绍常用于数值分析的一些命令。数值分析是指对于积分、微分或解析上难以确定的一些特殊值，就利用计算机在数值上来近似这些结果。

5.7.1求极小值

在许多应用中，确定函数的极值即极大值(峰值)和极小值(谷值)是较常见的。数学上，可通过确定函数导数(斜率)为零的点，解析求出这些极值点。然而，有很多函数很难找到导数为零的点，因此，必须通过在数值上来找函数的极值点。

MATLAB的函数fminbnd和fminsearch可以用来寻找最小值。这两个函数在MATLAB的早期版本中没有，早期版本中对应的分别为fmin和fmins函数。

最大值的求法则可以通过“f(x)的最大值=－f(x)的最小值”来得出。

1. fminbnd函数

fminbnd函数用来计算单变量非线性函数的最小值。 语法：

[x,y]=fminbnd(h_fun,x1,x2,options) [x,y]=fminbnd(‘funname’,x1,x2,options)

说明：h_fun是函数句柄，’funname’是函数名，必须是单值非线性函数；options是用来控制算法的参数向量，默认值为0可省略；x是fun函数在区间x1

【例5.22】用fminbnd求解humps函数的极小值。

>> [x,y]=fminbnd(@humps,0.5,0.8) x =

0.6370 y =

11.2528

%求在0.5—0.8之间极小值

程序分析：humps函数是MATLAB提供的M文件，保存为humps.m文件，@humps 表示humps函数的句柄，humps的曲线如图5.7所示，最小值为图中的圆点(0.6370，11.2528)，误差小于10-4。

图5.7 求humps函数最小值

2. fminsearch函数

fminsearch函数是求多变量无束缚非线性最小值，是采用Nelder-Mead单纯形算法求解多变量函数的最小值。

语法：

x=fminsearch(h_fun,x0) x=fminsearch(‘funname’,x0)

说明：x0是最小值点的初始猜测值。

【例5.23】求著名的Banana测试函数f(x,y)=100(y-x2)2+(1-x)2的最小值，它的理论最小值是x=1，y=1。该测试函数有一片浅谷，很多算法都难以逾越。

>> fn=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2' ,'x') %用inline产生内联函数,x和y用二元数组表示 fn =

Inline function:

fn(x) = 100*(x(2)-x(1)^2)^2+(1-x(1))^2 >> y=fminsearch(fn,[0.5,-1]) %从(0.5,-1)为初始值开始搜索求最小值 y =

1.0000 1.0000

5.7.2求过零点

函数f(x)的过零点或等于某个常数点的求解也非常重要。f(x)可能有一个零点，也可能有多个或无数个零点，也可能没有零点，用一般的解析方法寻找这类点非常困难或不可能。MATLAB提供了该类问题的数值解法。

fzero函数可以寻找一维函数的零点，即求f(x)=0的根。 语法：

x=fzero(h_fun,x0,tol,trace) x=fzero(‘funname’,x0,tol,trace)

说明：h_fun是待求零点的函数句柄；x0有两个作用：预定待搜索零点的大致位置和搜索起始点；tol用来控制结果的相对精度，默认值为eps；trace指定迭代信息是否在运算中显示，默认为0，表示不显示迭代信息。tol和trace都可以省略。

注意：一个函数可能会有多个零点，但fzero函数的结果只给出离x0最近的那个零点，如果fzero无法找出零点，它将停止运行并提供解释。fzero函数不能接受以x为自变量的字符串来描述的函数。fzero函数比求多项式根的roots函数能求出更一般的单变量函数的零点。

【例5.24】求解humps函数的过零点，humps函数如图5.8所示，过零点用圆点表示。

>> xzero=fzero(@humps,1) %求在1附近的零点 xzero = 1.2995

>> xzero=fzero(@humps,[0.5,1.5]) %求在0.5—1.5范围的零点 xzero = 1.2995

>> xzero=fzero(@humps,[0.5,1]) %求在0.5—1范围的零点 ??? Error using ==> fzero

The function values at the interval endpoints must differ in sign.

程序分析：当在0.5～1的范围找不到零点，提示出错信息。

图5.8 求humps函数的过零点

fzero函数不仅能寻找零点，它还可以寻找函数等于任何常数值的点。寻找f(x)=c的点，通过定义函数“g(x)=f(x)-c”，然后使用fzero找出g(x)为零的x值。

5.7.3数值积分

MATLAB提供了在有限区间内，计算出某个函数的积分或面积。

函数quad和quad8是基于数学上的正方形概念来计算函数的面积。这两个函数在所需的区间计算被积函数都应用递归调用的方法，quad8比quad更精确速度更快。与trapz梯形比较这两个函数能进行更高阶的近似。

语法：

s=quad(h_fun,x1,x2,tol,trace,p1,p2,…) s=quad(‘funname’,x1,x2,tol,trace,p1,p2,…) s=quad8(h_fun,x1,x2,tol,trace,p1,p2,…) s=quad8(‘funname’,x1,x2,tol,trace,p1,p2,…)

说明：x1和x2分别是积分的上、下限；tol用来控制积分精度，省略时默认为0.001；trace取1用图形展现积分过程，取0则无图形，省略时默认不画图；p1,p2?是向函数传递的参数，可省略。

【例5.25】计算y=humps(x)曲线下面的面积。

>> x=0:0.01:1; >> y=humps(x);

>> area=trapz(x,y) area =

29.8571

>> area1=quad(@humps,0,1) area1 = 29.8583

>> area2=quad8(@humps,0,1) area2 = 29.8583

%用梯形计算积分

%用quad计算积分

%用quad8计算积分

程序分析：用trapz函数梯形近似可能低估了实际面积，如果当梯形的宽度变窄时，就能够得到更精确的结果。quad和quad8函数返回非常相近的估计面积。

5.7.4微分方程的数值解

微分方程一般用来描述系统内部变量如何受系统内部结构和外部激励(如输入)的影响。当常微分方程式能够解析求解时，MATLAB提供了符号工具箱中的solve函数(3.6小节介绍)来找到精确解；如果求取解析解困难时，MATLAB为解决常微分方程提供了数值求解的方法。

MATLAB提供ode23、ode45和ode113等多个函数求解微分方程的数值解： ? 低维方法解一阶常微分方程组 语法：

[t,y]=ode23(h_fun,tspan,y0,options,p1,p2…) [t,y]=ode23(‘funname’,tspan,y0,options,p1,p2…) ? 高维方法解一阶常微分方程组 语法：

[t,y]=ode45(h_fun,tspan,y0,options,p1,p2…) [t,y]=ode45(‘funname’,tspan,y0,options,p1,p2…) ? 变维方法解一阶常微分方程组 语法：

[t,y]=ode113(h_fun,tspan,y0,options,p1,p2…) [t,y]=ode113(‘funname’,tspan,y0,options,p1,p2…)

说明：h_fun是函数句柄，函数以dx为输出，以t,y为输入量；tspan=[起始值 终止值]，表示积分的起始值和终止值；y0是初始状态列向量；options可以定义函数运行时的参数，可省略；p1,p2?是函数的输入参数，可省略。

一般来说，ode45求解算法最常用。ode23和ode45都运用了基本的龙格－库塔(Runge-Kutta)数值积分法的变形，ode23运用组合2／3阶龙格－库塔法，而ode45运用组合的4／5阶龙格－库塔法，ode45可取较多的时间步。

【例5.26】解经典的范德波尔(Van der Pol)微分方程：

d2xdt2?μ(1?x2)dx?x?0 dt(1) 必须把高阶微分方程式变换成一阶微分方程组。

令y1=x，y2=dx/dt，则将二阶微分方程变为一阶微分方程组： ?dy1??dt??y2??dy???? 2μ(1?y)y?y121??2???dt???(2) 编写一个函数vdpol.m文件，设定μ=2，该函数返回上述导数值。输出结果由一个列向量yprime给出。y1和y2合并写成列向量y。

函数M文件vdpol.m：

%范德波尔方程

function yprime=vdpol(t,y)

yprime=[y(2);2*(1-y(1)^2)*y(2)-y(1)]

(3) 给定当前时间及y1和y2的初始值，解微分方程：

>> tspan=[0,30]; >> y0=[1;0];

%起始值0和终止值30 %初始值

>> [t,y]=ode45(@vdpol,tspan,y0); %解微分方程 >> y1=y(:,1); >> y2=y(:,2); >> figure(1)

>>plot(t,y1,':b',t,y2,'-r') %画微分方程解 >>figure(2)

>>plot(y1,y2) %画相平面图

则微分方程y1和y2在时间域的曲线如图5.9所示。

(t,y(1))曲线 (t,y(2))曲线 图5.9 微分方程解

将y(1)为横坐标，y(2)为纵坐标则为相平面图，如图5.10所示。

图5.10 相平面图

本文来源：https://www.bwwdw.com/article/w4p8.html