1、第 8 章 MATLAB 程序设计8.1 脚本文件和函数文件M 文件有两种形式：M 脚本文件和 M 函数文件。8.1.1 M 文本编辑器MATLAB 的 M 文件是通过 M 文件编辑调试器窗口(EditorDebugger)来创建的。单击 MATLAB 桌面上的 图标，或者单击菜单 “File”“New”“M-file”，可打开空白的 M 文件编辑器，也可以通过打开已有的 M 文件来打开 M 文件编辑器。如图8.1 所示为打开已创建的 M 文件。8.1.2 M 文件的基本格式下面介绍绘制二阶系统时域曲线的 M 文件，欠阻尼系统的时域输出 y 与 x 的关系为，【例 8.1】为 M 脚本文件，【
2、例 8.2】为 M 函数文)cosax1sin(e1y2x2件。【例 8.1】用 M 脚本文件绘制二阶系统时域曲线。%EX0801 二阶系统时域曲线%画阻尼系数为 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) 图 8.1 M 文件编辑/调试器窗口【例 8.2】创建一个画二阶系统时域曲线的函数，阻尼系数 zeta 为函数的输入参数。function y=Ex0802(zeta)% EX0802 Step response of quadratic system
3、.% 二阶系统时域响应曲线% 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 Ex0802EX0802 Step response of quadratic system.二阶系统时域响应曲线zet
4、a 阻尼系数y 时域响应lookfor 二阶系统时域响应Ex0802.m: %二阶系统时域响应8.1.3 M 脚本文件脚本文件的特点：(1) 脚本文件中的命令格式和前后位置，与在命令窗口中输入的没有任何区别。(2) MATLAB 在运行脚本文件时，只是简单地按顺序从文件中读取一条条命令，送到MATLAB 命令窗口中去执行。(3) 与在命令窗口中直接运行命令一样，脚本文件运行产生的变量都是驻留在MATLAB 的工作空间(workspace)中，可以很方便地查看变量，除非用 clear 命令清除；脚本文件的命令也可以访问工作空间的所有数据，因此要注意避免变量的覆盖而造成程序出错。【例 8.1 续】
5、在 M 文件编辑调试器窗口中编写 M 脚本文件绘制二阶系统的多条时域曲线。(1) 单击 MATLAB 桌面上的 图标打开 M 文件编辑器。(2) 将命令全部写入 M 文件编辑器中，为了能标志该文件的名称，在第一行写入包含文件名的注释。保存文件为 Ex0801.m。%EX0801 二阶系统时域曲线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
6、(1-0.7072)*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”，就可以在图形窗中看到如图 8.2所示的曲线。查看工作空间的变量：whos Name Size Bytes Classx 1x201 1608 double arrayy1 1x201 1608 double arrayy2 1x201 1608 double arrayy3 1x201 1608 double arrayGrand total
7、 is 804 elements using 6432 bytes 8.1.4 M 函数文件函数文件的特点：(1) 第一行总是以“function ”引导的函数声明行；函数声明行的格式： function 输出变量列表 = 函数名(输入变量列表) (2) 函数文件在运行过程中产生的变量都存放在函数本身的工作空间；(3) 当文件执行完最后一条命令或遇到“return”命令时，就结束函数文件的运行，同时函数工作空间的变量就被清除；图 8.2 运行界面(4) 函数的工作空间随具体的 M 函数文件调用而产生，随调用结束而删除，是独立的、临时的，在 MATLAB 运行过程中可以产生任意多个临时的函数空间
8、。【例 8.2 续】在 M 文件编辑调试器窗口编写计算二阶系统时域响应的 M 函数文件，并在 MATLAB 命令窗口中调用该文件。创建 M 函数文件并调用的步骤如下：(1) 编写函数代码function y=Ex0802(zeta)%EX0802 画二阶系统时域曲线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) 将函数文件保存为“Ex0802.m”。(3) 在 MATLAB 命令窗口输入以下命令，则会出现 f 的计算值和绘制的曲线：f=Ex0802(0.3)程序分析：
9、 第一行指定该文件是函数文件，文件名为“Ex0802” ，输入参数为阻尼系数 zeta，输出参数为时域响应 y。 当函数文件调用结束，查看 x、y：x? Undefined function or variable x.y? Undefined function or variable y.注意：M 脚本文件和 M 函数文件的文件名及函数名的命名规则与 MATLAB 变量的命名规则相同。8.2 程序流程控制8.2.1 for . end 循环结构语法： for 循环变量=array循环体end 说明：循环体被循环执行，执行的次数就是 array 的列数，array 可以是向量也可以是矩阵，循环
10、变量依次取 array 的各列，每取一次循环体执行一次。【例 8.3】使用 for . end 循环的 array 向量编程求出 1+3+5.+100 的值。% EX0803 使用向量 for 循环sum=0;for n=1:2:100sum=sum+n;endsum sum =2500 计算的结果为：sum =2500。程序说明：循环变量为 n，n 对应为向量 1:2:100，循环次数为向量的列数，每次循环n 取一个元素。【例 8.4】使用 for . end 循环的 array 矩阵编程将单位阵转换为列向量。% EX0804 使用矩阵 for 循环sum=zeros(6,1);for n=
11、eye(6,6)sum=sum+n;endsum sum =111111 程序分析：循环变量 n 对应为矩阵 eye(6,6)的每一列，即第一次 n 为1;0;0;0;0;0，第一次 n 为0;1;0;0;0;0 ；循环次数为矩阵的列数 6。8.2.2 while . end 循环结构语法： while 表达式循环体end 说明：只要表达式为逻辑真，就执行循环体；一旦表达式为假，就结束循环。表达式可以是向量也可以是矩阵，如果表达式为矩阵则当所有的元素都为真才执行循环体，如果表达式为 nan， MATLAB 认为是假，不执行循环体。【例 8.5】与【例 8.3】相同，计算 1+3+5.+100
12、的值。% EX0805 使用 while 循环sum=0;n=1;while n0)n=1;while n Ex0814Too many input arguments.也可以在工作空间查看函数体定义的输入参数个数：nargin(Ex0814)ans =2【例 8.14 续】添加以下程序，查看用 nargout 变量获取输出参数个数。if nargout=0 %当输出参数个数为 0 时，运算结果为 0sum=0;end在命令窗口调用 Ex0814 函数，当输出参数格式不同时，结果如下：Ex0814(2,3) %当输出参数个数为 0 时 ans =0 y=Ex0814(2,3) %当输出参数个数
13、为 1 时 y =5 y,n,x=Ex0814 %当输出参数个数为太多时 ? Error using = Ex0814Too many output arguments. 程序分析：当输出参数个数为 0 时，即使有两个输入参数，运算结果也为 0，结果送给 ans 变量；当输出的参数个数太多，也会出错。(2) varargin 和 varargout 变量varargin 和 varargout 可以获得输入输出变量的各元素内容。【例 8.15】计算所有输入变量的和。function y,n=Ex0815(varargin)% EX0815 使用可变参数 vararginif nargin=0
14、%当没有输入变量时输出 0disp(No Input variables.)y=0;elseif nargin=1 %当一个输入变量时，输出该数y=varargin1;else n=nargin;y=0;for m=1:ny=vararginm+y; %当有多个输入变量时，取输入变量循环相加endendn=nargin; 在 MATLAB 的命令窗口中输入不同个数的变量调用函数 Ex0815，结果如下：y,n=Ex0815(1,2,3,4) %输入 4 个参数 y =10n =4 y,n=Ex0815(1) %输入 1 个参数 y =1n =1 y,n=Ex0815 %无输入参数 No Inp
15、ut variables.y =0n =0 程序分析：n 为输入参数的个数；y 为求和运算的结果。8.3.4 程序举例【例 8.16】根据阻尼系数绘制不同二阶系统的时域响应，当欠阻尼时，当临界阻尼时 ，当过阻尼时)cosax1sin(e1y2x2 xe)1(y。2x)1(2x)1(2eM 文件的程序代码如下：function y=Ex0816(z1)% EX0816 主函数调用子函数，根据阻尼系数绘制二阶系统时域曲线t=0:0.1:20;if (z1=0)y=humps(x);area=trapz(x,y) %用梯形计算积分 area =29.8571 area1=quad(humps,0,1
16、) %用 quad 计算积分 area1 =29.8583 area2=quad8(humps,0,1) %用 quad8 计算积分 Warning: QUAD8 is obsolete. QUADL is its recommended replacement.(Type “warning off MATLAB:quad8:ObsoleteFunction“ to suppress this warning.) In D:MATLABtoolboxmatlabfunfunquad8.m at line 35area2 =29.8583 程序分析：用 trapz 函数梯形近似可能低估了实际面积
17、，如果当梯形的宽度变窄时，就能够得到更精确的结果。quad 和 quad8 函数返回非常相近的估计面积。8.7.4 微分方程的数值解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
18、,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是函数的输入参数，可省略。【例 8.26】解经典的范德波尔(Van der Pol) 微分方程：0xdt)(1dtx22(1) 必须把高阶微分方程式变换成一阶微分方程组。令 y
19、1=x，y2=dx/dt，则将二阶微分方程变为一阶微分方程组： 1221y)(dty(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; %起始值 0 和终止值 30y0=1;0; %初始值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) %画相平面图
