《第5章 程序设计精选文档.ppt》由会员分享，可在线阅读，更多相关《第5章 程序设计精选文档.ppt（84页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第5章 程序设计本讲稿第一页，共八十四页5.1 脚本文件和函数文件5.2 程序流程控制5.3 函数调用和参数传递5.4 M文件性能的优化和加速5.5 内联函数5.6 利用函数句柄执行函数5.7 利用泛函命令进行数值分析本讲稿第二页，共八十四页5.1 5.1 脚本文件和函数文件脚本文件和函数文件5.1.1 M文本编辑器MATLAB的M文件是通过M文件编辑调试器窗口(EditorDebugger)来创建的。5.1.2 M文件的基本格式下面介绍绘制二阶系统时域曲线的M文件，欠阻尼系统的时域输出y与x的关系为 M文件有两种形式：M脚本文件和M函数文件。本讲稿第三页，共八十四页【例【例5.1】用M脚本文
2、件脚本文件绘制二阶系统时域曲线。%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.2】创建一个画二阶系统时域曲线的M函数文件函数文件，阻尼系数zeta为函数的输入参数。functiony=Ex0502(zeta)%EX0502 Step response of quadratic system.%二阶系统时域响应曲线%zeta 阻尼系数%y 时域响应x=0:0.1:20;y=1-1/s
3、qrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta)plot(x,y)本讲稿第五页，共八十四页M函数文件的基本格式：函数声明行H1行(用%开头的注释行)在线帮助文本(用%开头)编写和修改记录(用%开头)函数体 函数文件必须的，函数名与文件名一致，不一致时以文件名为准可供help和lookfor关键词查询使用通常包含函数输入输出变量的含义、格式说明一般在空一行后，记录作者、日期、版本记录，用于软件档案管理本讲稿第六页，共八十四页 lookfor 二阶系统时域响应Ex0502.m:%二阶系统时域响应例如，在命令窗口输入help和look
4、for命令查看帮助信息：help Ex0502 EX0502 Step response of quadratic system.二阶系统时域响应曲线 zeta 阻尼系数 y 时域响应本讲稿第七页，共八十四页5.1.3 M脚本文件脚本文件的特点：(1)脚本文件中的命令格式和前后位置，与在命令窗口中输入的没有任何区别。(2)MATLAB在运行脚本文件时，只是简单地按顺序从文件中读取一条条命令，送到MATLAB命令窗口中去执行。(3)与在命令窗口中直接运行命令一样，脚本文件运行产生的变量都是驻留在MATLAB的工作空间(workspace)中，可以很方便地查看变量，除非用clear命令清除；脚本文
5、件的命令也可以访问工作空间的所有数据，因此要注意避免变量的覆盖而造成程序出错。本讲稿第八页，共八十四页【例例5.1续续】在M文件编辑调试器窗口中编写M脚本文件绘制二阶系统的多条时域曲线。%EX0501 二阶系统时域曲线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的曲线holdony2=1-1/sqrt(1-0.7072)*exp(-0.707*x).*sin(sqrt(1-0.7072)*x+acos(0.707)plot(x,y2,g)%画阻尼系数为
6、0.707的曲线y3=1-exp(-x).*(1+x)plot(x,y3,b)%画阻尼系数为1的曲线本讲稿第九页，共八十四页5.1.4 M函数文件(3)当文件执行完最后一条命令或遇到“return”命令时，就结束函数文件的运行，同时函数工作空间的变量就被清除；(4)函数的工作空间随具体的M函数文件调用而产生，随调用结束而删除，是独立的、临时的，在MATLAB运行过程中可以产生任意多个临时的函数空间。函数文件的特点：(1)第一行总是以“function”引导的函数声明行；函数声明行的格式：函数声明行的格式：function 输出变量列表输出变量列表=函数名函数名(输入变量列表输入变量列表)(2)
7、函数文件在运行过程中产生的变量都存放在函数本身的工作空间；本讲稿第十页，共八十四页(2)将函数文件保存为“Ex0502.m”。(3)在MATLAB命令窗口输入以下命令，则会出现f的计算值和绘制的曲线：f=Ex0502(0.3)【例例5.2续续】在M文件编辑调试器窗口编写计算二阶系统时域响应的M函数文件，并在MATLAB命令窗口中调用该文件。functiony=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)创建M函数文件
8、并调用的步骤如下：(1)编写函数代码注意：注意：M脚本文件和M函数文件的文件名及函数名的命名规则与MATLAB变量的命名规则相同。本讲稿第十一页，共八十四页5.2 5.2 程序流程控制程序流程控制 支持各种流程控制结构：顺序结构 循环结构：forend whileend 条件结构(分支结构)：ifelseend switchcase 试探结构:trycatchend本讲稿第十二页，共八十四页语法：for 循环变量循环变量=array循环体循环体end 说明：循环体被循环执行，执行的次数就是array的列数，array可以是向量也可以是矩阵，循环变量依次取array的各列，每取一次循环体执行一次
9、。5.2.1 for.end循环结构本讲稿第十三页，共八十四页 【例例5.3】使用for.end循环的array向量编程求出 100内的奇数和。%EX0503 使用向量for循环sum=0;forn=1:2:100sum=sum+n;end本讲稿第十四页，共八十四页 【例例5.4】使用for.end循环的array矩阵编程将单位阵转换为列向量。%EX0504 sum=zeros(6,1);for n=eye(6,6)sum=sum+n;endsum=zeros(6,1)sum=000000n=eye(6,6)n=100000010000001000000100000010000001本讲稿第十
10、五页，共八十四页5.2.2 while.end循环结构说明：只要表达式为逻辑真，就执行循环体；一旦表达式为假，就结束循环。表达式可以是向量也可以是矩阵，如果表达式为矩阵则当所有的元素都为真才执行循环体，如果表达式为nan，MATLAB认为是假，不执行循环体。语法：语法：while 表达式表达式 循环体循环体end 本讲稿第十六页，共八十四页【例例5.5】计算100内的奇数和。%EX0505 使用while循环sum=0;n=1;whilen=100sum=sum+n;n=n+2;end本讲稿第十七页，共八十四页5.2.3 Ifelseend条件转移结构说明：当有多个条件时，条件式1为假再判断e
11、lseif的条件式2，如果所有的条件式都不满足，则执行else的语句段n+1，当条件式为真则执行相应的语句段；ifelseend结构也可以是没有elseif和else的简单结构。语法：if 条件式条件式1 语句段语句段1elseif 条件式条件式2 语句段语句段2 .else 语句段语句段n+1end本讲稿第十八页，共八十四页【例例5.6】用If结构执行二阶系统时域响应，根据阻尼系数0zeta0)&(zeta1)y=1-1/sqrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta);elseifzeta=1y=1-exp(-x).*(1
12、+x);endplot(x,y)本讲稿第十九页，共八十四页5.2.4 switchcase开关结构语法：switch 开关表达式开关表达式case 表达式表达式1 语句段语句段1case表达式表达式2 语句段语句段2.otherwise 语句段语句段nend 本讲稿第二十页，共八十四页说明：(1)将开关表达式依次与case后面的表达式进行比较，如果表达式1不满足，则与下一个表达式2比较，如果都不满足则执行otherwise后面的语句段n；一旦开关表达式与某个表达式相等，则执行其后面的语句段。(2)开关表达式只能是标量或字符串。(3)case后面的表达式可以是标量、字符串或元胞数组，如果是元胞数
13、组则将开关表达式与元胞数组的所有元素进行比较，只要某个元素与开关表达式相等，就执行其后的语句段。本讲稿第二十一页，共八十四页【例例5.7】用switchcase开关结构得出各月份的季节。%EX0507使用使用switch结构结构formonth=1:12;switchmonthcase3,4,5season=springcase6,7,8season=summercase9,10,11season=autumnotherwiseseason=winterendend本讲稿第二十二页，共八十四页5.2.5 try.catch.end试探结构语法：try 语句段语句段1catch 语句段语句段2e
14、nd 说明：首先试探性地执行语句段1，如果在此段语句执行过程中出现错误，则将错误信息赋给保留的lasterr变量，并放弃这段语句，转而执行语句段2中的语句，当执行语句段2又出现错误，则终止该结构。本讲稿第二十三页，共八十四页【例例5.8】用try.catch.end结构来进行矩阵相乘运算。%EX0508 try结构n=4;a=magic(n);m=3;b=eye(3);tryc=a*bcatchc=a(1:m,1:m)*bendlasterr 本讲稿第二十四页，共八十四页5.2.6 流程控制语句1.break命令命令break命令可以使包含break的最内层的for或while语句强制终止，立
15、即跳出该结构，执行end后面的命令，break命令一般和if结构结合使用。break,continue,return,pause,keyboard,input本讲稿第二十五页，共八十四页【例5.9】将【例5.5】增加条件用If与break命令结合，停止while循环。计算1+3+5.+100 的值，当和大于1000时终止计算。%EX0509用用break终止终止while循环循环sum=0;n=1;whilen=100ifsum a=input(inputanumber:)%输入数值给a b=input(inputanumber:,s)%输入字符串给b本讲稿第三十页，共八十四页5.3 5.3
16、函数调用和参数传递函数调用和参数传递使用函数可以把一个比较大的任务分解为多个较小的任务，使程序模块化，每个函数完成特定的功能，通过函数的调用完成整个程序的功能。本讲稿第三十一页，共八十四页5.3.1 子函数和私有函数1.子函数子函数 在一个M函数文件中，可以包含一个以上的函数，其中只有一个是主函数，其它则为子函数。(1)在一个M文件中，主函数必须出现在最上方，其后是子 函数，子函数的次序无任何限制；(2)子函数不能被其它文件的函数调用，只能被同一文件中的函数(可以是主函数或子函数)调用；(3)同一文件的主函数和子函数变量的工作空间相互独立；(4)用help和lookfor命令不能提供子函数的帮
17、助信息。本讲稿第三十二页，共八十四页【例例5.11】将画二阶系统时域曲线的函数作为子函数，编写画多条曲线的程序。functionEx0511()%EX0511 使用函数调用绘制二阶系统时域响应z1=0.3;Ex0502(z1);%调用Ex0502holdonz1=0.5Ex0502(z1)%调用Ex0502z1=0.707;Ex0502(z1)%调用Ex0502function y=Ex0502(zeta)%子函数，画二阶系统时域曲线x=0:0.1:20;y=1-1/sqrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta)plot(x
18、,y)本讲稿第三十三页，共八十四页(2)私有函数父目录的M脚本文件也不可调用私有函数；(3)在函数调用搜索时，私有函数优先于其它MATLAB路径上的函数。(1)在private目录下的私有函数，只能被其父目录的M函数文件所调用，而不能被其它目录的函数调用，对其它目录的文件私有函数是不可见的，私有函数可以和其它目录下的函数重名；2.私有函数私有函数 私有函数是指存放在private子目录中的M函数文件，具有以下性质：本讲稿第三十四页，共八十四页在MATLAB中调用一个函数，搜索的顺序如下：（1）查找是否子函数；（2）查找是否私有函数；（3）从当前路径中搜索此函数；（4）从搜索路径中搜索此函数。3
19、.3.调用函数的搜索顺序调用函数的搜索顺序本讲稿第三十五页，共八十四页5.3.2 局部变量和全局变量1.局部变量局部变量局部变量(Local Variables)是在函数体内部使用的变量，其影响范围只能在本函数内，只在函数执行期间存在。2.全局变量全局变量全局变量(Global Variables)是可以在不同的函数工作空间和MATALB工作空间中共享使用的变量。全局变量在使用前必须用global进行定义，且每个要共享全局变量的函数和工作空间，都必须逐个用global对变量加以定义。本讲稿第三十六页，共八十四页【例【例5.12】修改【例】修改【例5.11】在主函数和子函数中使用全局变量。】在主
20、函数和子函数中使用全局变量。functionEx0512()%EX0512使用全局变量绘制二阶系统时域响应使用全局变量绘制二阶系统时域响应globalXX=0:0.1:20;z1=0.3;Ex0502(z1);holdonz1=0.5;Ex0502(z1);z1=0.707;Ex0502(z1);functionEx0502(zeta)%子函数，画二阶系统时域曲线子函数，画二阶系统时域曲线globalXy=1-1/sqrt(1-zeta2)*exp(-zeta*X).*sin(sqrt(1-zeta2)*X+acos(zeta);plot(X,y);如果在工作空间中定义X为全局变量也可以使用。
21、注意：注意：由于全局变量在任何定义过的函数中都可以修改，因此不提倡使用全局变量；使用时应十分小心，建议把全局变量的定义放在函数体的开始，全局变量用大写字符命名。本讲稿第三十七页，共八十四页5.3.3 函数的参数函数的调用过程即是参数的传递过程。函数的调用格式：输出参数输出参数1 1，输出参数，输出参数2 2，函数名（输入参数函数名（输入参数1 1，输入参数，输入参数2 2，）1.参数传递规则参数传递规则 函数具有自己的工作空间，函数内变量与外界的唯一联系是通过输入、输出参数。输入参数在函数内的变化不会传递出去。本讲稿第三十八页，共八十四页functionx,y=Ex0502(zeta)%子函数
22、，画二阶系统时域曲线子函数，画二阶系统时域曲线x=0:0.1:20;y=1-1/sqrt(1-zeta2)*exp(-zeta*x).*sin(sqrt(1-zeta2)*x+acos(zeta);functionEx0513()z1=0.3;x1,y1=Ex0502(z1);plot(x1,y1)holdonz1=0.5;x2,y2=Ex0502(z1);plot(x2,y2)z1=0.707;x3,y3=Ex0502(z1);plot(x3,y3);【例例5.13】将【例5.11】画二阶系统时域的函数修改，使用输入输出参数来实现参数传递 本讲稿第三十九页，共八十四页2.函数参数的个数函数参
23、数的个数Matlab函数的输入输出参数的个数都可以变化，用户可根据参数的个数来编程。(1)nargin和nargout变量(the number of argument)函数的输入输出参数的个数可以通过变量nargin和nargout获得，nargin用于获得输入参数的个数，nargout用于获得输出参数的个数。语法：语法：nargin%在函数体内获取实际输入变量的个数nargout%在函数体内获取实际输出变量的个数nargin(fun)%在函数体外获取定义的输入参数个数nargout(fun)%在函数体外获取定义的输出参数个数本讲稿第四十页，共八十四页【例例5.14】计算两个数的和，根据输入
24、的参数个数不同使用不同的运算表达式。function sum,n=Ex0514(x,y)%EX0514 参数个数可变，计算x和y的和if nargin=1 sum=x+0;%输入一个参数就计算与0的和elseif nargin=0 sum=0;%无输入参数就输出0else sum=x+y;%输入的是两个数则就计算和end本讲稿第四十一页，共八十四页(2)varargin和varargout变量(variable)varargin和varargout可以获得输入输出变量的各元素内容,它们都是元胞数组。functiony,n=Ex0515(varargin)%EX0515 使用可变参数vararg
25、inifnargin=0%当没有输入变量时输出当没有输入变量时输出0disp(NoInputvariables.)y=0;elseifnargin=1%当一个输入变量时，输出该数当一个输入变量时，输出该数y=varargin1;elsen=nargin;y=0;form=1:ny=vararginm+y;%当有多个输入变量时，取输入变量循环相加当有多个输入变量时，取输入变量循环相加endendn=nargin;本讲稿第四十二页，共八十四页5.3.4 程序举例【例例5.16】根据阻尼系数绘制不同二阶系统的时域响应，当欠阻尼时，当临界阻尼时，当过阻尼时，本讲稿第四十三页，共八十四页M文件的程序代码
26、如下：文件的程序代码如下：functiony=Ex0516(z1)%EX0516主函数调用子函数，%根据阻尼系数绘制二阶系统时域曲线t=0:0.1:20;if(z1=0)&(z1f=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 char(f)ans=sin(x)*exp(-z*x)class(f)ans=inline argnames(f)ans=x z 本讲稿第六十三页，共八十四页3.使内联函数适用于数组运算内联
27、函数的输入变量不能是数组，输入变量不能是数组，但可以使用vectorize命令将内联函数适用于数组运算。语法：语法：vectorize(inline_fun)%使内联函数适用于数组运算 4.执行内联函数 内联函数还可以直接使用feval命令执行。语法：语法：y1，y2,=feval(inline_fun,arg1,arg2)本讲稿第六十四页，共八十四页 x=0:0.1:20;z=0:0.05:10;y=feval(ff,x,z)%执行内联函数 ff=vectorize(f)%使内联函数使内联函数f转换为适合数组运算转换为适合数组运算 ff=Inline function:ff(x,z)=sin
28、(x).*exp(-z.*x)x=0:0.1:20;y=ff(x,0.3);本讲稿第六十五页，共八十四页5.6 5.6 利用函数句柄执行函数利用函数句柄执行函数5.6.1 函数句柄（Function_Handle）的创建 函数句柄包含了函数的路径、函数名、类型以及可能存在的重载方法。它是Matlab6.0后启用的新数据类型。1.函数句柄的创建语法：语法：h_fun=fun%创建函数句柄h_fun=str2func(fun)%创建函数句柄h_array=str2func(fun1,fun2,)%创建函数句柄数组说明：fun是函数名，h_fun是函数句柄，h_array是函数句柄数组。本讲稿第六十
29、六页，共八十四页【例例5.20】创建MATLAB内部函数的句柄。h_sin=sin;%创建函数句柄 h_cos=str2func(cos);%创建函数句柄数组 h_array=str2func(sin,cos,tan)h_array=sin cos tan 本讲稿第六十七页，共八十四页(1)在更大范围调用函数函数句柄包含了函数文件的路径和函数类型，即函数是否为内部函数、M或P文件、子函数、私有函数等，因此无论函数所在的文件是否在搜索路径上，是否是当前路径，是否是子函数或私有函数，只要函数句柄存在，函数只要函数句柄存在，函数就能执行就能执行。2.使用函数句柄的优点(2)提高函数调用的速度不使用函
30、数句柄时，对函数的每次调用都要为该函数进行全面的路径搜索，直接影响了速度。(3)使函数调用象使用变量一样方便、简单。(4)可迅速获得同名重载函数的位置、类型信息。本讲稿第六十八页，共八十四页5.6.2 用feval命令执行函数函数也可以使用feval命令直接执行，feval命令可以使用函数句柄或函数名。语法：语法：y1,y2,=feval(h_fun,arg1,arg2)y1,y2,=feval(funname,arg1,arg2)说明：h_fun是函数句柄，funname是函数名，arg1、arg2是输入参数，y1、y2是输出参数。本讲稿第六十九页，共八十四页function y=Ex052
31、1(z1)%EX0521 利用函数句柄执行函数，二阶系统时域响应t=0:0.1:20;h_plotxy1=str2func(plotxy1)%创建函数句柄h_plotxy2=str2func(plotxy2)%创建函数句柄h_plotxy3=str2func(plotxy3)%创建函数句柄if(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 编写的绘制二阶系统时域响应曲线中的调用各子函数改为利用函数句柄实现
32、。本讲稿第七十页，共八十四页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
33、-sqrt(zeta2-1)-exp(-(zeta+sqrt(zeta2-1)*x)./(zeta+sqrt(zeta2-1);plot(x,y3)本讲稿第七十一页，共八十四页(1)用feval命令利用函数句柄执行 h_Ex0521=str2func(Ex0521)y=feval(h_Ex0521,1);(2)用feval命令利用函数名执行 y=feval(Ex0521,1);(3)直接调用函数 y=Ex0521(1);在MATLAB的命令窗口调用该Ex0521函数有三种格式：本讲稿第七十二页，共八十四页5.7 5.7 利用泛函命令进行数值分析利用泛函命令进行数值分析在MATLAB中，所有以函
34、数为输入变量的命令，都称为泛函命令泛函命令。常见语法：输出变量列表输出变量列表=函数名函数名(h_fun,输入变量列表输入变量列表)输出变量列表输出变量列表=函数名函数名(funname,输入变量列表输入变量列表)说明：h_fun是要被执行的M函数文件的句柄，或者是内联函数和字符串；funname是M函数文件名。下面介绍一些用于数值分析的泛函命令。数值分析数值分析是指对于微分、积分或解析上难以确定的一些特殊值，利用计算机在数值上近似这些结果。本讲稿第七十三页，共八十四页5.7.1 求极小值1.fminbnd函数函数fminbnd函数用来计算单变量非线性单变量非线性函数的最小值。语法：语法：x,
35、y=fminbnd(h_fun,x1,x2,options)x,y=fminbnd(funname,x1,x2,options)说明：h_fun是函数句柄，funname是函数名，必须是单值非线性函数；options是用来控制算法的参数向量，默认值为0可省略；x是fun函数在区间x1xfn=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)y=1.0000 1.0000 【
36、例【例5.23】求著名的Banana测试函数f(x,y)=100(y-x2)2+(1-x)2的最小值，它的理论最小值是x=1，y=1。该测试函数有一片浅谷，很多算法都难以逾越。本讲稿第七十七页，共八十四页5.7.2 求过零点fzero函数可以寻找一维函数的零点，即求f(x)=0的根。语法：语法：x=fzero(h_fun,x0,tol,trace)x=fzero(funname,x0,tol,trace)说明：h_fun是待求零点的函数句柄；x0有两个作用：预定待搜索零点的大致位置和搜索起始点；tol用来控制结果的相对精度，默认值为eps；trace指定迭代信息是否在运算中显示，默认为0，表示
37、不显示迭代信息。tol和trace都可以省略。本讲稿第七十八页，共八十四页5.7.3 数值积分说明：x1和x2分别是积分的上、下限；tol用来控制积分精度，省略时默认为0.001；trace取1用图形展现积分过程，取0则无图形，省略时默认不画图；p1,p2是向函数传递的参数，可省略。函数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,tr
38、ace,p1,p2,)s=quad8(funname,x1,x2,tol,trace,p1,p2,)本讲稿第七十九页，共八十四页5.7.4 微分方程的数值解 当常微分方程式能够解析求解时，符号工具箱中的solve函数可得到精确解；当解析求解困难时，MATLAB提供ode23、ode45和ode113等多个函数求解微分方程的数值解。本讲稿第八十页，共八十四页2.高维方法解一阶常微分方程组高维方法解一阶常微分方程组（最常用）（最常用）语法：语法：t,y=ode45(h_fun,tspan,y0,options,p1,p2)t,y=ode45(funname,tspan,y0,options,p1,
39、p2)3.变维方法解一阶常微分方程组变维方法解一阶常微分方程组语法：语法：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是函数的输入参数，可省略。1.低维方法解一阶常微分方程组低维方法解一阶常微分方程组语法：语法：t,y=ode23(h_fun,tspan,y0,options,p1
40、,p2)t,y=ode23(funname,tspan,y0,options,p1,p2)本讲稿第八十一页，共八十四页【例【例5.26】解经典的范德波尔(Van der Pol)微分方程：(1)必须把高阶微分方程式变换成一阶微分方程组。令y1=x，y2=dx/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;%起始值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)%画相平面图本讲稿第八十四页，共八十四页