资源描述
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的详细运行细节。
展开阅读全文