教学PPTMATLAB程序设计

5 . 1 5 . 1 程序流程控制程序流程控制 5 . 25 . 2 M文件文件 5 . 3 5 . 3 函数调用和参数传递函数调用和参数传递 5 . 4 5 . 4 利用函数句柄执行函数利用函数句柄执行函数 5 . 5 5 . 5 利用泛函命令进行数值分析利用泛函命令进行数值分析 5 . 6 5 . 6 内联函数内联函数 5 . 7 M5 . 7 M文件性能的优化和加速文件性能的优化和加速 5.1.1 for5.1.1 forendend循环结构循环结构 MATLAB的循环结构有2种：forend结构和whileend结构。 语法： for循环变量=array 循环体 end 【例5.1】 使用forend循环编程求出1+3+5+99的值。 sum=0; for n=1:2:100 sum=sum+n end 【例5.2】 使用forend循环将单位阵转换为列向量。 sum=zeros(6,1); for n=eye(6,6) sum=sum+n end 计算结果如下。 sum= 1 1 1 1 1 1 forend循环的循环次数确定，而whileend循环的循环次数不确定。 语法： while表达式 循环体 end 【例5.3】 根据 ，求y3时的最大n值和y值。 y=0; n=1; while y3 y=y+1/(2*n-1); n=n+1; z(n)=y; end mn=n-2 %y=8 y(n)=1; elseif x(n)=6 y(n)=5-x(n)/2; elseif x(n)=4 y(n)=2; else y(n)=sqrt(x(n); end end plot(x,y) axis(0 10 0 2.5); switchcase结构是有多个分支结构的条件转移结构。 语法： switch开关表达式 case表达式1 语句段1 case表达式2 语句段2 otherwise 语句段n end 【例5.5】 用switchcase开关结构得出各月份的季节。 for month=1:12; switch month case3,4,5 season=spring case6,7,8 season=summer case9,10,11 season=autumn otherwise season=winter end end MATLAB提供一种试探结构trycatchend，这种语句结构是其他很多语言所没有的。 语法： try 语句段1 catch 语句段2 end 【例5.6】 用trycatchend结构进行矩阵相乘运算。 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= Errorusing=* Inner matrix dimensions must agree. 程序分析：试探出矩阵的大小不匹配时，矩阵无法相乘，则再执行catch后面的语句段，将A A的子矩阵取出与B B矩阵相乘。 1break命令命令 break命令可以使包含break的最内层的for或while语句强制终止，立即跳出该结构，执行end后面的命令，break命令一般和if结构结合使用。 【例5.7】 将【例5.3】增加条件，用if与break命令结合，停止while循环。计算 值，当y=3时终止计算。 y=0; n=1; while n=100 if y3 y=y+1/(2*n-1); n=n+1; z(n)=y; else break end end mn=n-2 %ya=input(input a number:) %输入数值给a input a number:45 a= 45 b=input(input a number:,s) %输入字符串给b input a number:45 b= 45 input(input a number:) %将输入值进行运算 input a number:2+3 ans= 5 1 1以电影方式产生动画以电影方式产生动画 以电影方式产生动画，有2个步骤。 （1）使用getframe命令抓取图形作为画面，每个画面都是以1个列向量的方式，置于存放整个电影的矩阵M M中。 （2）使用movie(M,k)命令播放电影，并可指定矩阵M M播放的重复次数k。 【例5.9】 使用电影方式制作动画并显示，最后一幅画面如图5.2所示。 t=-pi:0.1:pi; n=length(t) for ii=1:n %循环次数由横坐标长度决定 x=-pi:0.1:t(ii); y=sin(2*x); z=cos(x); plot3(x,y,z) %绘制曲线 axis(-4,4,-1,1,-1,1); %确定坐标轴范围 M(ii)=getframe; %抓取画面 end movie( M) %播放数组M的画面 2 2以对象方式产生动画以对象方式产生动画 （1）擦除属性EraseMode。以对象方式产生动画需要设置EraseMode属性。EraseMode为1个字符串，代表对象的擦除方式，即对于旧对象的处理方式。EraseMode属性有以下几种。 normal：计算整个画面的数据，重画整个图形。 xor：将旧对象的点以xor的方式还原，即只画与屏幕色不一致的新对象点，擦除不一致的原对象点，这种方式不会擦除被擦对象下面的其他图像。 background：将旧对象的点变成背景颜色，实现擦除，这种方式会擦除被擦对象下的其他图像。 none：保留旧对象的点，不进行任何擦除。 （2）对象的位置属性。通常在动画过程中，会改变对象的位置、尺寸或颜色等外观属性。位置属性有如下2种。 xdata：为1个向量，代表对象的x坐标值。 ydata：为1个向量，代表对象的y坐标值。 （3）屏幕刷新。当新对象的属性设置后，应刷新屏幕，使新对象显示出来，刷新屏幕用drawnow命令实现。drawnow命令使MATLAB暂停当前的任务序列而去刷新屏幕，如果没有drawnow命令，MATLAB则会等当前的任务序列执行完才去刷新屏幕。 （4）产生动画。产生动画的具体步骤是：首先产生1个对象，其EraseMode属性为xor、background或none；然后在循环中产生动画，每次循环改变此对象的xdata或ydata（或两者）；最后使用drawnow命令刷新屏幕。 【例5.10】 使用对象方式产生1个红色的小球沿着曲线运行的界面，如图5.3所示。 x=0:0.1:20; y=11/sqrt(10.32)*exp(0.3*x).*sin(sqrt(10.32)*x+acos(0.3); plot(x,y) h=line(0,0,color,red,marker,.,markersize,40,erasemode,xor);%定义红色的小球 for i=1:length(x) set(h,xdata,x(i),ydata,y(i); %设置小球的新位置 pause(0.005) %暂停0.005 s drawnow %刷新屏幕 end 运行界面 5.2.1 M5.2.1 M文件编辑器文件编辑器 单击MATLAB桌面上的 图标，或者选择菜单“File”“New”“Script”命令，或选择“Function”可以分别打开空白的M脚本文件编辑器和M函数文件编辑器。如 图5.4（a）所示为空白的M脚本文件编辑器，如图5.4（b）所示为M函数文件编辑器。 脚本文件具有如下特点。 （1）脚本文件中的命令格式和前后位置，与在命令窗口中输入时没有任何区别。 （2）MATLAB在运行脚本文件时，只是简单地按顺序从文件中读取一条条命令，送到MATLAB命令窗口中去执行。 （3）与在命令窗口中直接运行命令一样，脚本文件运行产生的变量都是驻留在MATLAB的工作空间（workspace）中，可以很方便地查看变量，除非用clear命令清除；脚本文件的命令也可以访问工作空间的所有数据，为此要注意避免变量的覆盖而造成程序出错。 【例5.11】 在M文件编辑/调试器窗口中编写M脚本文件绘制二阶系统的多条时域曲线。欠阻尼系统的时域输出y与x的关系为 。 （1）单击MATLAB桌面上的 图标打开M文件编辑器。 （2）将命令全部写入M文件编辑器中，为了能够标志该文件的名称，在第1行写入包含文件名的注释。保存文件为Ex0501.m。 % EX0511 绘制二阶系统的时域曲线 x=0:0.1:20; y1=11/sqrt(10.32)*exp(0.3*x).*sin(sqrt(10.32)*x+acos(0.3) plot(x,y1,r) %画阻尼系数为0.3的曲线 hold on y2=11/sqrt(10.7072)*exp(0.707*x).*sin(sqrt(10.7072)*x+acos(0.707) plot(x,y2,g) %画阻尼系数为0.707的曲线 y3=1exp(x).*(1+x) plot(x,y3,b) %画阻尼系数为1的曲线 2211esin( 1arccos )1xyx （3）选择M文件编辑器菜单“Debug”“Run”命令，运行界面如图5.5所示。 查看工作空间的变量如下： whos Name Size Byte Class x 1x201 1608 doublearray y1 1x201 1608 doublearray y2 1x201 1608 doublearray y3 1x201 1608 doublearray Grand total is 804 elements using 6432 Byte M函数文件的基本格式如下： 函数声明行 H1行（用%开头的注释行） 在线帮助文本（用%开头） 编写和修改记录（用%开头） 函数体 说明： （1）函数声明行以“function”引导，是M函数文件必须有的，M脚本文件没有；函数名和文件名一致，当不一致时，MATLAB以文件名为准，Ex0502函数保存在Ex0502.m文件中。 函数声明行的格式： function 输出变量列表=函数名（输入变量列表） （2）H1行通常包含大写的函数文件名，可以提供help和look for关键词用于查询。 （3）在线帮助文本通常包含函数输入、输出变量的含义和格式说明。 （4）编写和修改记录一般在空1行后，记录作者、日期和版本记录，用于软件档案管理。 （5）函数体由MATLAB的命令或者通过流程控制结构组织的命令组成。通过函数体实现函数的功能。 在MATLAB界面中选择菜单“File”“New”“Function”命令，则会创建一个新的函数文件并自动生成函数语句，函数文件内容如下： function output_args = Untitled1( input_args ) %UNTITLED1 Summary of this function goes here % Detailed explanation goes here end 【例5.12】 在M文件编辑/调试器窗口编写计算二阶系统时域响应的M函数文件，并在MATLAB命令窗口中调用该文件。 创建M函数文件并调用的步骤如下。 （1）修改上面的程序内容，修改output_args，Untitled1和input_args的名称。 function y=Ex0512(zeta) %EX0512 Step response of quadratic system. %Second order output curve %copyright 20110501 x=0:0.1:20; y=11/sqrt(1zeta2)*exp(zeta*x).*sin(sqrt(1zeta2)*x+acos(zeta) plot(x,y) （2）将函数文件保存为“Ex0512.m”。 在命令窗口输入help和lookfor命令查看帮助信息： help Ex0512 EX0512 Step response of quadratic system. Second order output curve lookfor Second order output curve Ex0512.m: %二阶系统时域响应 （3）在MATLAB命令窗口输入以下命令，则会出现f的计算值和绘制的曲线： f=Ex0512(0.3) 第1行指定该文件是函数文件，文件名为“Ex0512”，输入参数为阻尼系数zeta，输出参数为时域响应y。 当函数文件调用结束，查看x、y。 x ?Undefined function or variablex. y ?Undefined function or variabley. 5.3.1 5.3.1 子函数和私有函数子函数和私有函数 1 1子函数子函数 在1个M函数文件中，可以包含1个以上的函数，其中只有一个是主函数，其他均为子函数。 （1）在1个M文件中，主函数必须出现在最上方，其后是子函数，子函数的次序无任何限制。 （2）子函数不能被其他文件的函数调用，只能被同一文件中的函数（可以是主函数或子函数）调用。 （3）同一文件的主函数和子函数变量的工作空间相互独立。 （4）用help和lookfor命令不能提供子函数的帮助信息。 【例5.13】 将画二阶系统时域曲线的函数作为子函数，编写画多条曲线的程序。 function Ex0513() %EX0513 使用函数调用绘制二阶系统时域响应 z1=0.3; Ex0512(z1); %调用Ex0502 hold on z1=0.5 Ex0512(z1) %调用Ex0502 z1=0.707; Ex0512(z1) %调用Ex0502 functiony=Ex0512(zeta) %子函数，画二阶系统时域曲线 x=0:0.1:20; y=11/sqrt(1zeta2)*exp(zeta*x).*sin(sqrt(1zeta2)*x+acos(zeta) plot(x,y) 2 2私有函数私有函数 私有函数是指存放在private子目录中的M函数文件，具有以下性质。 （1）在private目录下的私有函数，只能被其父目录的M函数文件所调用，而不能被其他目录的函数调用。私有函数对其他目录的文件是不可见的，私有函数可以和其他目录下的函数重名。 （2）私有函数父目录的M脚本文件也不可调用私有函数。 （3）在函数调用搜索时，私有函数优先于其他MATLAB路径上的函数。 3 3调用函数的搜索顺序调用函数的搜索顺序 在MATLAB中调用1个函数，搜索的顺序如下。 （1）查找是否为子函数。 （2）查找是否为私有函数。 （3）从当前路径中搜索此函数。 （4）从搜索路径中搜索此函数。 1 1局部变量局部变量 局部变量（Local Variables）是在函数体内部使用的变量，其影响范围只能在本函数内；每个函数在运行时，都占用独立的函数工作空间，此工作空间和MATLAB的工作空间是相互独立的。局部变量仅存在于函数的工作空间内，只在函数执行期间存在，当函数执行完变量就消失。 2 2全局变量全局变量 全局变量（Global Variables）是可以在不同的函数工作空间和MATLAB工作空间中共享使用的变量。 【例5.14】 修改【例5.13】在主函数和子函数中使用全局变量。 function Ex0514() %EX0514 使用全局变量绘制二阶系统时域响应 global X X=0:0.1:20; z1=0.3; Ex0512(z1); hold on z1=0.5; Ex0512(z1); z1=0.707; Ex0512(z1); function Ex0512(zeta) %子函数，画二阶系统时域曲线 global X y=11/sqrt(1zeta2)*exp(zeta*X).*sin(sqrt(1zeta2)*X+acos(zeta); plot(X,y); 程序分析：X变量为全局变量，在需要使用的主函数和子函数中都需要用global定义；同样，如果在工作空间中定义X为全局变量后也可以使用。 global X who Your variables are: X 函数调用的格式如下： 输出参数1，输出参数2，=函数名（输入参数1，输入参数2，） 1 1参数传递规则参数传递规则 在MATLAB中，函数具有自己的工作空间，函数内变量与外界（包括其他函数和工作空间）的唯一联系就是通过函数的输入/输出参数。输入参数在函数中的任何变化，都仅在函数内进行，不会传递回去。 【例5.15】 将画二阶系统时域的函数修改，使用输入/输出参数来实现参数传递，如图5.6所示。 2 2函数参数的个数函数参数的个数 MATLAB函数的调用有一个与其他语言不同的特点是：函数的输入/输出参数的数目都可以变化，用户可以根据参数的个数编程。 （1）nargin变量和nargout变量。在MATLAB中有2个特殊变量：nargin和nargout。函数的输入/输出参数的个数可以通过变量nargin和nargout获得，nargin用于获得输入参数的个数，nargout用于获得输出参数的个数。 语法： nargin %在函数体内获取实际输入变量的个数 nargout %在函数体内获取实际输出变量的个数 nargin(fun) %在函数体外获取定义的输入参数个数 nargout(fun) %在函数体外获取定义的输出参数个数 【例5.16】 计算2个数的和，根据输入参数个数的不同使用不同的运算表达式。 functionsum,n=Ex0516(x,y) %EX0516 参数个数可变，计算x和y的和 if nargin=1 sum=x+0; %若输入1个参数就计算与0的和 else if nargin=0 sum=0; %若无输入参数就输出0 else sum=x+y; %若输入的是2个数则计算和 end 在命令窗口调用Ex0516函数，分别使用2个、1个和无输入参数，结果如下所示。 y,n=Ex0516(2,3) y= 5 n= 2 y,n=Ex0516(2) y= 2 n= 1 y,n=Ex0516 y= 0 n= 0 如果输入的参数多于输入参数个数，则会出错。 y,n=Ex0516(1,2,3) ?Errorusing=ex0516 Toomanyinputarguments. 也可以在工作空间查看函数体定义的输入参数个数。 nargin(Ex0516) ans= 2 【例5.16续】 添加以下程序，查看用nargout变量获取输出参数个数。 if nargout=0 %当输出参数个数为0时，运算结果为0 sum=0; end 在命令窗口调用Ex0516函数，当输出参数格式不同时，结果如下。 Ex0516(2,3) %当输出参数个数为0时 ans= 0 y=Ex0516(2,3) %当输出参数个数为1时 y= 5 y,n,x=Ex0516 %当输出参数个数太多时 ?Errorusing=ex0516 Too many out put arguments （2）varargin变量和varargout变量。MATLAB还有2个特殊变量：varargin和varargout，可以获得输入/输出变量的各元素内容，varargin和varargout都是元胞数组。 【例5.16续】 计算所有输入变量的和。 functiony,n=Ex0516(varargin) %EX0516 使用可变参数varargin if nargin=0 %当没有输入变量时输出0 disp(NoInputvariables.) y=0; elseif nargin=1 %当一个输入变量时，输出该数 y=varargin1; else n=nargin; y=0; for m=1:n y=vararginm+y; %当有多个输入变量时，取输入变量循环相加 end end n=nargin; 在MATLAB的命令窗口中输入不同个数的变量调用函数Ex0516，结果如下。 y,n=Ex0516(1,2,3,4) %输入4个参数 y= 10 n= 4 y,n=Ex0516(1) %输入1个参数 y= 1 n= 1 y,n=Ex0516 %无输入参数 NoInputvariables. y= 0 n= 0 【例5.17】 根据输入的参数绘制不同的曲线，每个子函数绘制一种曲线。 M文件的程序代码。 程序分析：主函数名为Ex0517，4个子函数名分别为po1、po2、po3和po4，文件保存为Ex0517.m。若在命令窗口中输入以下命令： y=Ex0517(1);y=Ex0517(1); 则产生如图5.7所示的调用第一个子函数的极坐标曲线。 【例5.18】 编写M函数文件，通过流程控制语句，建立如下的矩阵。 0123.0012.10001.2.0000.0nnnYfunctiony=Ex0518(m) %EX0518 用循环流程控制语句创建矩阵 y=0; m=m1; 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 在命令窗口中调用Ex0518函数。 y=Ex0518(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.19】 在本书第4章例4.30中，在界面中绘制正弦波形，在单选按钮中选择振幅，在列表框中选择频率，然后绘制正弦曲线。 1 1设计界面设计界面 界面设计如图5.8所示。 2 2回调函数回调函数 回调函数需要在单击“Draw”按钮（pushbutton1）时绘制正弦曲线，单击“Close”按钮（pushbutton1）时关闭窗口。 编写回调函数的步骤：选择“Draw”按钮（pushbutton1）单击鼠标右键，然后选择菜单“View Callbacks”“Callback”命令，进入pushbutton1_Callback函数，在函数中添加函数内容如下： function pushbutton1_Callback(hObject, eventdata, handles) am_1=get(handles.radiobutton1,value);%获取单选按钮radiobutton1的值 am_2=get(handles.radiobutton2,value); if am_1=1 %如果单选按钮radiobutton1被选中 am=5; else am=10; end fre1=get(handles.listbox1,string); %获取下拉列表listbox1的string属性 n=get(handles.listbox1,value); %获取下拉列表listbox1被选中的列表索引 fre=eval(fre1n); %得出下拉列表listbox1被选中项的内容 x=0:0.1:10; plot(am*sin(fre*x); %绘制正弦曲线，振幅和频率根据选择得出 handles变量非常重要，它是1个结构数组，包括以下两部分内容。 （1）存储所有在图形界面中的控件、菜单、坐标轴对象的句柄，每个对象的句柄名称与对象的“Tag”名相同，如下拉列表listbox1的句柄是“handles.listbox1”。 （2）通过get函数可以获得不同对象的属性值，例如，get(handles.listbox1,string)得出91的元胞数组，记录了所有下拉项的值，fre1n用来获取对应下拉列表项的值；get(handles.listbox1,value)得出用户所选的下拉项的索引，如果选择第一项则为1。 运行界面如图5.9所示。 5.4.1 5.4.1 函数句柄的创建函数句柄的创建 1 1函数句柄的创建函数句柄的创建 与图形对象的句柄不同，函数句柄不是在函数文件创建时自动创建的，而是必须通过专门的定义。创建函数句柄使用“”符号或str2func命令实现。 语法： h_fun=fun %创建函数句柄 h_fun=str2func(fun) %创建函数句柄 h_array=str2func(fun1, fun2, ) %创建函数句柄数组 【例5.20】 创建MATLAB内部函数的句柄。 h_sin=sin; %创建函数句柄 h_cos=str2func(cos); %创建函数句柄数组 h_array=str2func(sin,cos,tan) h_array= sin cos tan 2 2使用函数句柄的优点使用函数句柄的优点 利用函数句柄执行函数的优点有以下几点。 （1）在更大范围内调用函数。函数句柄包含了函数文件的路径和函数类型，即函数是否为内部函数、M或P文件、子函数、私有函数等。 （2）提高函数调用的速度。不使用函数句柄时，对函数的每次调用都要为该函数进行全面的路径搜索，直接影响了速度。 （3）使函数调用像使用变量一样方便、简单。 （4）可迅速获得同名重载函数的位置、类型信息。 函数也可以使用feval命令直接执行，feval命令可以使用函数句柄或函数名。 语法： y1,y2, =feval(h_fun,arg1,arg2) y1,y2, =feval(funname,arg1,arg2) 【例例5.21】 根据阻尼系数绘制不同二阶系统的时域响应，当 时， 01在MATLAB的命令窗口调用该Ex0520函数有3种格式。 （1）用feval命令利用函数句柄执行。 h_Ex0520=str2func(Ex0520) h_Ex0520= Ex0520 y=feval(h_Ex0520,1); （2）用feval命令利用函数名执行。 y=feval(Ex0520,1); （3）直接调用函数。 y=Ex0520(1); 5.5.1 5.5.1 求极小值求极小值 1 1fminbndfminbnd函数函数 fminbnd函数用来计算单变量非线性函数的极小值。 语法： x,y=fminbnd(h_fun,x1,x2,options) x,y=fminbnd(funname,x1,x2,options) 【例5.22】 用fminbnd求解humps函数的极小值。 x,y=fminbnd(humps,0.5,0.8) %求在0.50.8之间的极小值 x= 0.6370 y= 11.2528 程序分析：humps函数是MATLAB提供的M文件，保存为humps.m文件；humps表示humps函数的句柄，humps的函数最小值曲线如图5.10所示，最小值为图中的圆点(0.6370，11.2528)，误差小于104。 2 2fminsearchfminsearch函数函数 fminsearch函数用于求多变量无束缚非线性最小值，采用Nelder-Mead单纯形算法求解多变量函数的最小值。 语法： x=fminsearch(h_fun,x0) x=fminsearch(funname,x0) 【例5.23】 求著名的Banana测试函数f(x,y)=100(yx2)2+(1x)2的最小值，它的理论最小值是x=1，y=1。该测试函数有一片浅谷，很多算法都难以逾越。 fn=inline(100*(x(2) x(1)2)2+(1x(1)2,x) %用inline产生内联函数，x和y用二元数组表示 fn= Inline function: fn(x)=100*(x(2) x(1)2)2+(1x(1)2 y=fminsearch(fn,0.5, 1) %从(0.5,1)为初始值开始搜索求最小值 y= 1.0000 1.0000 fzero函数可以寻找一维函数的零点，即求f(x)=0的根。 语法： x=fzero(h_fun,x0,tol,trace) x=fzero(funname,x0,tol,trace) 【例5.24】 求解humps函数的过零点，humps函数的过零点用圆点表示，如图5.11所示。 xzero=fzero(humps,1) %求在1附近的零点 xzero= 1.2995 xzero=fzero(humps,0.5,1.5) %求在0.51.5范围内的零点 xzero= 1.2995 xzero=fzero(humps,0.5,1) %求在0.51范围内的零点 ?Errorusing=fzero The function valuesat the interval end points must differin sign. humps函数的过零点 语法： 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, ) 【例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计算积分 area2= 29.8583 MATLAB提供了ode23、ode45和ode113等多个函数求解微分方程的数值解。以下介绍低维方法、高维方法和变维方法解一阶常微方程组。 语法： t,y=ode45(h_fun,tspan,y0,options,p1,p2) t,y=ode45(funname,tspan,y0,options,p1,p2) 【例5.26】 解经典的范德波尔（Van derPol）微分方程： 222dd(1)0ddxxxxtt（1）必须把高阶微分方程式变换成为一阶微分方程组。 令y1=x，y2=dx/dt，则将二阶微分方程变为一阶微分方程组： （2）编写1个函数vdpol.m文件，设定=2，该函数返回上述导数值。输出结果由1个列向量yprime给出。y1和y2合并写成列向量y。 函数M文件vdpol.m： %范德波尔方程 functiony prime=vdpol(t,y) yprime=y(2);2*(1y(1)2)*y(2) y(1) 1221212dd(1)ddyytyyyyt（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) %画相平面图 设y1为横坐标，y2为纵坐标，则画出相平面图如图5.12所示。 1 1内联函数的创建内联函数的创建 创建内联函数可以使用inline命令实现。 语法： inline(string,arg1,arg2, ) %创建内联函数 【例5.26】 创建内联函数实现 。 f=inline(sin(x)*exp(z*x),x,z) %创建内联函数 f= Inline function: f(x,z)=sin(x)*exp(z*x) y=f(5,0.3) %调用函数f y= 0.2140 2 2查看内联函数查看内联函数 MATLAB可以用char、class和argnames命令方便地查看内联函数的信息。 语法： char(inline_fun) %查看内联函数的内容 class(inline_fun) %查看内联函数的类型 argnames(inline_fun) %查看内联函数的变量 【例5.26续】 查看内联函数的信息。 char(f) ans= sin(x)*exp(z*x) class(f) ans= inline argnames(f) ans= x z 3 3使内联函数适用于数组运算使内联函数适用于数组运算 内联函数的输入变量不能是数组，但可以使用vectorize命令将内联函数适用于数组运算。 语法： vectorize(inline_fun) %使内联函数适用于数组运算 【例5.26续】 使内联函数适用于数组运算。 ff=vectorize(f) %使内联函数f转换为适合于数组运算 ff= Inline function: ff(x,z)=sin(x).*exp(z.*x) x=0:0.1:20; y=ff(x,0.3); 4 4执行内联函数执行内联函数 内联函数还可以直接使用feval命令执行。 语法： y1，y2, =feval(inline_fun,arg1,arg2) 【例5.26续】 执行内联函数。 x=0:0.1:20; z=0:0.05:10; y=feval(ff,x,z) 5.7.1 M5.7.1 M文件性能优化文件性能优化 1 1使用循环时提高速度的措施使用循环时提高速度的措施 循环语句及循环体是MATLAB编程的瓶颈问题。MATLAB与其他编程语言不同，它的基本数据是向量和矩阵，编程时应尽量对向量和矩阵编程，而不要对矩阵的元素编程。 2 2大型矩阵的预先定维大型矩阵的预先定维 由于MATLAB变量在使用之前不需要定义和指定维数，当变量新赋值的元素下标超出数组的维数时，MATLAB就为该数组扩维1次，大大地降低了运行的效率。 【例5.27】 将【例5.18】中的双重循环改为单循环，并用zeros函数定维来提高运行速度，创建矩阵 。 0123.0012.10001.2.0000.0nnnYfunctiony=Ex0518(m) %EX0518先定维再创建矩阵 m=m1; y=zeros(m); for n=1:m1 a=1:mn; y(n,n+1:m)=a; end y 3优先考虑内在函数优先考虑内在函数 矩阵运算应该尽量采用MATLAB的内在函数。内在函数是由更底层的C语言构造的，其执行速度显然很快。 4采用高效的算法采用高效的算法 在实际应用中，解决同样的数学问题经常有各种各样的算法。例如，求解定积分的数值解法在MATLAB中就提供了2个函数：quad和quad8，其中后者在精度、速度上都明显高于前者。可见，应寻求更高效的算法。 5尽量使用尽量使用M函数文件代替函数文件代替M脚本文件脚本文件 由于M脚本文件每次运行时，都必须把程序装入内存，然后逐句解释执行，十分费时。因此，要尽量使用M函数文件代替M脚本文件。 1 1P P码文件的生成码文件的生成 P码文件使用pcode命令生成，生成的P码文件与原M文件名相同，其扩展名为“.p”。 语法： pcode Filename.m %在当前目录生成Filename.p pcode Filename.m-inplace %在Filename.m所在目录生成Filename.p 例如，将【例5.27】生成P码文件，在命令窗口中输入代码将Ex0527.m文件转换为P码文件。 pcode Ex0527.m 则在当前目录就生成了P码文件Ex0518.p。 2 2P P码文件的特点码文件的特点 （1）P码文件的运行速度比原M文件速度快。 （2）若存在同名的M文件和P码文件时则P码文件被调用。 为了测试，可以将Ex0527.m文件和Ex0527.p文件放在同一目录，将Ex0527.m文件修改并保存，然后在命令窗口运行该文件，结果仍然与原来一样。 y=Ex0527(5) （3）P码文件保密性好。用字处理软件打开Ex0527.p文件，看到的是乱码。 （1）打开程序性能剖析窗口。选择MATLAB菜单“Desktop”“Profiler”命令，或在命令窗口输入“profile viewer”命令都可以打开程序性能剖析窗口，如图5.13所示。或者在M文件编辑器窗口中选择菜单“Tools”“Open Profiler”命令进行剖析。 将本书的【例5.27】和【例5.18】进行剖析比较，在图5.13中输入“y=Ex0527(100)”，然后单击“Start Profiling”按钮，查看剖析报告。图5.14（a）显示的是单循环【例5.27】的剖析报告，图5.14（b）显示的是双循环【例5.18】的剖析报告。 （3）查看剖析报告。在程序性能分析窗口中以表格显示“剖析分析汇总表”（Profile Summary），从上到下按占用时间的多少排列。 在图5.14中可以看出创建100100的矩阵，单循环只需要0.020s，而双重循环需要0.032s。 如果需要得出分析的详细内容，则可以单击表格中的超链接，打开相应的内容页面。在图5.14（a）图中单击Ex0527则打开图5.15，可以看到Ex0527的详细运行细节。