2022年第章MATLAB程序设计 .pdf

第 5章 MATLAB 程序设计5.1 脚本文件和函数文件M 文件有两种形式：M 脚本文件和M 函数文件。5.1.1 M 文本编辑器MATLAB的 M 文件是通过M 文件编辑调试器窗口(EditorDebugger)来创建的。单击MATLAB桌面上的图标，或者单击菜单“File”“New”“M-file”，可打开空白的M 文件编辑器，也可以通过打开已有的M 文件来打开M 文件编辑器。如图5.1 所示为打开已创建的M 文件。5.1.2 M 文件的基本格式下面介绍绘制二阶系统时域曲线的M 文件，欠阻尼系统的时域输出y 与 x 的关系为)cosax1sin(e111y2x2，【例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.1 M 文件编辑 /调试器窗口名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 31 页 - - - - - - - - - 【例 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-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 二阶系统时域曲线名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 31 页 - - - - - - - - - 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 on y2=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 的曲线(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 array Grand total is 804 elements using 6432 bytes 5.1.4 M 函数文件函数文件的特点：(1) 第一行总是以“function ”引导的函数声明行；函数声明行的格式：function 输出变量列表 = 函数名 (输入变量列表) (2) 函数文件在运行过程中产生的变量都存放在函数本身的工作空间；图 5.2 运行界面名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 31 页 - - - - - - - - - (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; end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 31 页 - - - - - - - - - sum 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; end sum 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 ; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 31 页 - - - - - - - - - end sum n sum = 2500 n = 101 程序分析：可以看出while . end 循环的循环次数由表达式来决定，当n=101 就停止循环。5.2.3 Ifelse end条件转移结构语法：if 条件式 1 语句段 1 elseif 条件式 2 语句段 2 . else 语句段 n+1 end 说明：当有多个条件时，条件式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 ex0514 Too many input arguments. 也可以在工作空间查看函数体定义的输入参数个数： nargin(Ex0514) ans = 2 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 14 页，共 31 页 - - - - - - - - - 【例 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 = ex0514 Too many output arguments. 程序分析：当输出参数个数为0 时，即使有两个输入参数，运算结果也为0，结果送给 ans变量；当输出的参数个数太多，也会出错。(2) 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=varargin1; else n=nargin; y=0; for m=1:n y=vararginm+y; %当有多个输入变量时，取输入变量循环相加 end end n=nargin; 在 MATLAB的命令窗口中输入不同个数的变量调用函数Ex0515，结果如下： y,n=Ex0515(1,2,3,4) %输入 4 个参数y = 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 15 页，共 31 页 - - - - - - - - - 10 n = 4 y,n=Ex0515(1) % 输入 1 个参数y = 1 n = 1 y,n=Ex0515 %无输入参数No Input variables. y = 0 n = 0 程序分析： n 为输入参数的个数；y 为求和运算的结果。5.3.4程序举例【 例5.16 】 根 据 阻 尼 系 数 绘 制 不 同 二 阶 系 统 的 时 域 响 应 ， 当 欠 阻 尼 时)cosax1sin(e111y2x2， 当 临 界 阻 尼 时xe)x1(1y， 当 过 阻 尼 时2x)1(2x)1(21e1e1211y22。M 文件的程序代码如下：function y=Ex0516(z1) % EX0516 主函数调用子函数，根据阻尼系数绘制二阶系统时域曲线t=0:0.1:20; if (z1=0)&(z1=0)&(z11) y=feval(h_plotxy1,z1,t); % 执行函数elseif z1=1 y=feval(h_plotxy2,z1,t); % 执行函数else y=feval(h_plotxy3,z1,t); % 执行函数end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 24 页，共 31 页 - - - - - - - - - function 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) 程序分析：在主函数中使用函数句柄执行各子函数，运行的结果与【例5.17】相同。在 MATLAB的命令窗口调用该Ex0521 函数有三种格式：(1) 用 feval 命令利用函数句柄执行 h_Ex0521=str2func(Ex0521) h_Ex0521 = Ex0521 y=feval(h_Ex0521,1); h_plotxy1 = plotxy1 h_plotxy2 = plotxy2 h_plotxy3 = plotxy3 (2) 用 feval 命令利用函数名执行 y=feval(Ex0521,1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 25 页，共 31 页 - - - - - - - - - 0246810121416182000.20.40.60.811.21.4(3) 直接调用函数 y=Ex0521(1); 5.7 利用泛函命令进行数值分析在 MATLAB中，所有以函数为输入变量的命令，都称为泛函命令。常见语法：输出变量列表=函数名 (h_fun, 输入变量列表) 输出变量列表=函数名 ( funname ,输入变量列表) 说明： h_fun 是要被执行的M 函数文件的句柄，或者是内联函数和字符串； funname是 M 函数文件名。5.7.1求极小值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 函数在区间x1x fzero The function values at the interval endpoints must differ in sign. 程序分析：当在0.51的范围找不到零点，提示出错信息。图 5.8 求 humps函数的 过零点名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 28 页，共 31 页 - - - - - - - - - 5.7.3数值积分函数 quad 和 quad8是基于数学上的正方形概念来计算函数的面积，quad8 比 quad更精确速度更快。语法：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) % 用 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 35 area2 = 29.8583 程序分析：用trapz 函数梯形近似可能低估了实际面积，如果当梯形的宽度变窄时，就能够得到更精确的结果。quad和 quad8 函数返回非常相近的估计面积。5.7.4微分方程的数值解MATLAB提供 ode23、ode45 和 ode113等多个函数求解微分方程的数值解：低维方法解一阶常微分方程组名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 29 页，共 31 页 - - - - - - - - - 语法：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, 是函数的输入参数，可省略。【例 5.26】解经典的范德波尔(Van der Pol)微分方程：0 xdtdx)x(1dtxd222(1) 必须把高阶微分方程式变换成一阶微分方程组。令 y1=x，y2=dx/dt ，则将二阶微分方程变为一阶微分方程组：1221221y)yy(1ydtdydtdy(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 和终止值 30 y0=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) %画相平面图名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 30 页，共 31 页 - - - - - - - - - 则微分方程y1 和 y2 在时间域的曲线如图5.9 所示。将 y(1)为横坐标， y(2)为纵坐标则为相平面图，如图5.10 所示。图 5.10 相平面图(t,y(2) 曲 线(t,y(1) 曲线图 5.9 微分方程解名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 31 页，共 31 页 - - - - - - - - -