1、第第5章章MATLAB程序设计程序设计 MATLAB有两种常用工作方式：（1)直接交互的指令行操作方式；（2)M文件中的编程工作方式。5.1M文件中的功能和特点文件中的功能和特点(1)MATLAB程序文件中一个ASC码文件，扩展名一律为“.m”，可用任何字处理软件进行编辑和修改。(2)MATLAB是程序性的语言是解释性语言。(3)M文件大大扩展了MATLAB的能力 5.2 M文件的建立与使用文件的建立与使用1、M文件的形式文件的形式M文件中有两种形式(1)命令文件 （Script File）也叫做脚本文件(2)函数文件（Function File）2、脚本、脚本M文件文件（1）建立 用记事本（
2、Notebook）编写程序；在MATLAB系统下选择：“file”“new”“m-file”（2）内容：按照程序的功能，依据MATLAB的程序结构，组织的合法的MATLAB命令。(4)注意注意（1）命令文件中的语句可以访问MATLAB工作区(workspace)中的所有变量。运行过程中所有产生的变量均是全局变量。（2）“%”开始的行为注释行，不予执行。（3）若将文件未存放在系统的搜索目录下时，运用该文件之前，应先进入该目录。（如cd d:mywok）在命令窗口直接键入 文件名.m(5)举例：画出函数图形举例：画出函数图形n脚本m文件的内容：y=1./(1+exp(-x)，保存文件名jb1.mn
3、在工作窗口下定义：x=0:1/20:20 运行脚本文件（jb1）得到函数值yn命令窗口下，用plot(x,y)画出对应的图形 3、函数文件、函数文件(1)建立用记事本（Notebook）编写程序在MATLAB系统下选择：“file”“new”“m-file”(2)函数文件的格式 Function 输出参数=函数名（输入参数表）函数体(3)函数体的内容 函数定义行 Hi行：帮助信息的第一行，用于提示函数的功能，当用命令 Lookfor查询该函数的帮助信息时将显示该行内容。帮助体 函数体 n函数名命名规则同变量名。以字母、下划命名规则同变量名。以字母、下划线和数字线和数字 组成，不识别函数名共组成
4、，不识别函数名共3131个字符。个字符。n输入参数：用中括号：用中括号“”括起来，参数括起来，参数两两之间用逗号隔开。两两之间用逗号隔开。n输出参数：函数：函数m m文件的运算结果传到调用文件的运算结果传到调用处，当参数不只一个时，用逗号隔开。处，当参数不只一个时，用逗号隔开。n函数文件名由函数名再加后缀由函数名再加后缀“m”m”组成。组成。n函数文件中定义变量为函数文件中定义变量为局部变量局部变量，函数，函数文件命令运行结束，该类变量会自动释文件命令运行结束，该类变量会自动释放。如果变量在程序运行前就已存在的放。如果变量在程序运行前就已存在的话，程序运行后它不会受到影响话，程序运行后它不会受
5、到影响(4)注意注意（5）应用举例）应用举例n函数表示函数表示 函数m文件的内容为：function y=foft(x)y=1./(1+exp(-x)文件名为：foft fplot(foft,0 20,1e-04)n数学函数曲线的绘图数学函数曲线的绘图 fplot(fun，lims)其中：fun表示函数名，定义函数的M文件名lims=XMIN XMAX YMIN YMAXn函数的极小值点和零值点函数的极小值点和零值点l 单变量函数的局部极小值点 fmin(F，x1，x2，options)F函数名，x1,x2指定极小值区间，options第一个分量为正时，则显示函数的运行步骤。Options（1
6、）屏幕上是否显示最小值的迭代过程，非0时，显示，为0时不显示，默认为0 Options（2）迭代时自变量的误差控程限，默认为1.0e-04 Options（3）为迭代时函数值的误差控程限，函数fmins使用该参数 Options（4）为迭代时最大的迭代次数，函数fmin默认值为500，fmins为200步。举例：定义函数文件humps.mfunction y=humps(x)y=1./(x-0.3).2+0.01)+1./(x-0.9).2+0.04)-6figure(1)subplot(2,2,1),fplot(humps,-5,5),grid onsubplot(2,2,2),fplot(
7、humps,-5,5 -10 25),grid onsubplot(2,2,3),fplot(5*sin(x),-1,1),grid onsubplot(2,2,3),fplot(5*sin(x),humps,-1,1)grid onx=fmin(humps,0.3,1)极小值点 x函数的极小值为：humps（x）l 函数的零点fzero(F，猜测值)得到零点的区间l 多变量函数的局部极小值点fmins(函数名，初始极值向量，options)得到向量附近的一个极小值举例：定义函数 function b=two_var(v)x=v(1);y=v(2)b=x.2+y.2-x.2*y.2给定初始值：
8、v=-0.7 1.2调用函数：a=fmins(two_var,v)options与单变量时相同 n 数值积分数值积分l quad(F,A,B,TOL,TRACE)其中：F为函数名，A和B指定积分区间，TOL迭代误差，TRACE，采用的递归的辛浦生方法。l quad 8(F,A,B,TOL,TRACE)其中：F为函数名，A和B指定积分区间，TOL迭代误差，TRACE，采用的newton-cotes方法。l dblquad(F,inmin,inmax,outnmin,outinmax)其中：inmin与inmax，分别为内层积分的上下限，outnmin与outinmax，分别为外层积分的上下限举例
9、：定义函数 function y=ex542f(x)y=-x.*x+115 调用函数：s=quad(ex542f,0,10)例 建立一个命令文件将变量a,b的值互换，然后运行该命令文件。首先建立命令文件并以文件名exch.m存盘：clear;a=1:10;b=11,12,13,14;15,16,17,18;c=a;a=b;b=c;a b 然后在MATLAB的命令窗口中输入exch，将会执行该命令文件。例 建立一个函数文件将变量a,b的值互换，然后在命令窗口调用该函数文件。首先建立函数文件fexch.m：function a,b=exch(a,b)c=a;a=b;b=c;然后在MATLAB的命令
10、窗口调用该函数文件：clear;x=1:10;y=11,12,13,14;15,16,17,18;x,y=fexch(x,y)5.2 数据的输入输出1、input函数格式：A=input(提示信息，选项)；其中：提示信息为一个字符串，用于提示用户输入什么样的数据。如果在input函数调用时采用s选项，则允许用户输入一个字符串。例如，想输入一个人的姓名，可采用命令：xm=input(Whats your name?,s)2、disp函数函数格式格式:disp(输出项输出项)其中输出项既可以为字符串，也可以为矩阵。其中输出项既可以为字符串，也可以为矩阵。注意：注意：用用disp函数显示矩阵时函数显
11、示矩阵时将不显示矩阵将不显示矩阵的名字，而且其格式更紧密，且不留任何没的名字，而且其格式更紧密，且不留任何没有意义的空行。有意义的空行。例 1 求一元二次方程ax2+bx+c=0的根。程序如下：a=input(a=?);b=input(b=?);c=input(c=?);d=b*b-4*a*c;x=(-b+sqrt(d)/(2*a),(-b-sqrt(d)/(2*a);disp(x1=,num2str(x(1),x2=,num2str(x(2);5.3 M文件的程序结构文件的程序结构1、顺序结构顺序结构2、分支结构分支结构 （1）if （条件）（命令组）end （2)if（条件）（命令组1）e
12、lse （命令2）end当条件成立时，则执行语句组，执行完之后继续执行if语句的后继语句，若条件不成立，则直接执行if语句的后继语句 例例2 删除数据中的无关数据：一般认为，误差大于倍标准差的数据项，就认为是无关项，因此可以从数据中删除，为此可删除所找列的无关项if exist(c:matlabbksj.txt)load c:matlabbksj.txt i,j=size(sj)mu=mean(sj)sigma=std(sj)outliers=abs(sj-ones(i,1)*mu)3*ones(i,1)*sigma sj(outliers)=else sj=rand(40,3)end例例3：
13、输入一个字符，若为大写字母，则输出其后继字符，若为小写输入一个字符，若为大写字母，则输出其后继字符，若为小写字母，则输出其前导字符，若为数字字符则输出其对应的数值，字母，则输出其前导字符，若为数字字符则输出其对应的数值，若为其他字符则原样输出。若为其他字符则原样输出。程序如下程序如下：c=input(c=input(请输入一个字符请输入一个字符,s);,s);if c=A&c=A&c=a&c=a&c=0&c=0&c1Z=0.54*exp(-0.75*X.2-3.75*Y.2-1.5*Y);elseifT-1andT=1Z=0.7575*exp(-X.2-6*Y.2);elseZ=0.5457*
14、exp(-0.75*X.2-3.75*Y.2+1.5*Y);endsurf(X,Y,Z)%plot3(X,Y,Z)%mesh(X,Y,Z)15457.0117575.0154.0),(5.175.375.065.175.375.0222222yxeyxeyxeyxpyyxyxyyx 4、switch语句格式为：switch 表达式 case 表达式1 语句组1 case 表达式2 语句组2 case 表达式m 语句组m otherwise 语句组m+1 end例例5:5:某商场对顾客所购买的商品实行打折销售，标准如下(商品价格用price来表示)：price200 没有折扣 200price5
15、00 3%折扣 500price1000 5%折扣 1000price2500 8%折扣 2500price0 y=x.2+1 subplot(1,2,1)plot(x,y)hold on else y=x.2-1 subplot(1,2,2)plot(x,y)hold onendpauseend例例4：用循环语句编写用循环语句编写M文件计算文件计算ex的值，其中的值，其中x,n为输入变量，为输入变量，ex的近似表达式为的近似表达式为2312!3!nxxxxexn function y=e(x,n)y=1;s=1;for i=1:n s=s*i;y=y+xi/s;endyy=e(1,100)y
16、=2.7183 for语句更一般的格式：for 循环变量=矩阵表达式 循环体语句 end 执行过程是依次将矩阵的各列元素赋给循环变量，然后执行循环体语句，直至各列元素处理完毕。实际上，“表达式1:表达式2:表达式3”是一个仅为一行的矩阵(行向量)，因而列向量是单个数据。例5：已知5个学生4门功课的成绩，求每名学生的总成绩。程序如下：s=0;a=65,76,56,78;98,83,74,85;76,67,78,79;98,58,42,73;67,89,76,87;for k=a s=s+k;enddisp(s);（2）While 表达式 end 其执行过程为：若条件成立，则执行循环体语句，执行后
17、再判断条件是否成立，如果不成立则跳出循环例例6：根据矩阵指数的幂级数展开式求矩阵指数。根据矩阵指数的幂级数展开式求矩阵指数。程序如下：程序如下：X=input(X=);E=zeros(size(X);F=eye(size(X);n=1;whilenorm(F,1)0E=E+F;F=F*X/n;n=n+1;endEexpm(X)%调用调用MATLAB矩阵指数函数求矩阵矩阵指数函数求矩阵指数指数 循环的嵌套n如果一个循环结构的循环体又包括一个循环结构，就称为循环的嵌套，或称为多重循环结构。可以按照嵌套层数，分别叫做二重循环、三重循环等。处于内部的循环叫作内循环，处于外部的循环叫作外循环。n在设计多
18、重循环时，要特别注意内、外循环之间的关系，以及各语句放置的位置，不要搞错。例7 用筛选法求某自然数范围内的全部素数。程序如下：m=input(m=);p=2:m;for i=2:sqrt(m)n=find(rem(p,i)=0&p=i);p(n)=;end p综合实例综合实例 用用 MATLAB 求解问题时，一般要经历求解问题时，一般要经历建模建模和和编程编程两个过两个过程，只有在建模正确的前提下，方能得出正确的结果。程，只有在建模正确的前提下，方能得出正确的结果。n单自由度系统有阻尼自由振动单自由度系统有阻尼自由振动1.1.建立计算模型建立计算模型 由动力学可知，单自由度有阻尼自由振动的由动
19、力学可知，单自由度有阻尼自由振动的振动方程振动方程为：为：无量刚化无量刚化后有：后有：其中其中,上述上述方程的解为方程的解为：其中其中 x0 表示初始位置，表示初始位置，v0 表示初始速度。表示初始速度。0kxxcxm 0 xxxnn mknkmc2)sin()(dtnAetx220200)()(ddnxxvA0001xvxtgnd21nd参数参数n 10，x0 1，v0 0，计算的终止时间，计算的终止时间 t=2。试求。试求 从从 0.1 到到 1运动方程的解，并画出波形运动方程的解，并画出波形。2.MATLAB 编程编程 编写编写 M 文件文件 ex1.m%首先清空首先清空 MATLAB
20、的工作空间的工作空间clear;%给定初值给定初值 wn=10;tf=2;x0=1;v0=0;%计算不同的计算不同的 值所对应的振型值所对应的振型 for j=1:10;eta(j)=0.1*j;wd(j)=wn*sqrt(1-eta(j)2);%求振幅求振幅 A a=sqrt(wn*x0*eta(j)+v0)2+(x0*wd(j)2)/wd(j);综合实例综合实例%求相位角求相位角phi=atan2(wd(j)*x0,v0+eta(j)*wn*x0);%设定自变量数组设定自变量数组 t t=0:tf/1000:tf;%求过渡过程求过渡过程 x(j,:)=a*exp(-eta(j)*wn*t)
21、.*sin(wd(j)*t+phi);end%在同一个图形窗口中绘制不同的在同一个图形窗口中绘制不同的 值所对应的振型值所对应的振型 plot(t,x(1,:),t,x(2,:),t,x(3,:),t,x(4,:),.t,x(5,:),t,x(6,:),t,x(7,:),t,x(8,:),.t,x(9,:),t,x(10,:)grid on%新建一个图形窗口，绘制三维网格图新建一个图形窗口，绘制三维网格图 figure mesh(x)综合实例综合实例综合实例综合实例如果改变初始条件令如果改变初始条件令x00，v01，其运动曲线实际上就是，其运动曲线实际上就是系统的脉冲过渡函数系统的脉冲过渡函数
22、。综合实例综合实例2、气体分子运动的麦克斯韦分布曲线、气体分子运动的麦克斯韦分布曲线 通过本例说明如何通过本例说明如何用复杂的数学公式绘制曲线用复杂的数学公式绘制曲线。利用气体分子运动的麦克斯韦速度分布律，求氯分子运动利用气体分子运动的麦克斯韦速度分布律，求氯分子运动的的速度分布曲线速度分布曲线，并讨论，并讨论温度温度T及及分子量分子量mu对速度分布曲对速度分布曲线的影响。线的影响。(1)建立计算模型建立计算模型 麦克斯韦速度分布律为麦克斯韦速度分布律为:其中，其中，m-分子质量分子质量,m=mu/NA,mu-分子量分子量,NA-阿伏加德罗数阿伏加德罗数 k-波尔茨曼常数波尔茨曼常数 T-气体
23、的绝对温度气体的绝对温度 v-分子速度分子速度kTmvvkTmf2exp242223综合实例综合实例为研究单个参数的影响，先把麦克斯韦分布律编为一个为研究单个参数的影响，先把麦克斯韦分布律编为一个函函数子程序数子程序，以便重复调用，同时将，以便重复调用，同时将常数项常数项也放在子程序中。也放在子程序中。需要强调的是：需要强调的是：子程序不得与主程序放在同一个子程序不得与主程序放在同一个 M 文件中文件中，只能将子，只能将子程序单独做成程序单独做成 M 文件，并放在与主程序同一个工作路径中。文件，并放在与主程序同一个工作路径中。2.MATLAB 编程编程 首先建立计算麦克斯韦分布律的子程序首先建
24、立计算麦克斯韦分布律的子程序 mxw.mfunction f=mxw(T,mu,v)%The subfunction mxw.m of ex2 利用麦克斯韦速度分布律求分子的速度分布曲线的子程利用麦克斯韦速度分布律求分子的速度分布曲线的子程序序%mu、v、T分别是分子量、分子速度和气体的绝对温度分别是分子量、分子速度和气体的绝对温度 k=1.381*10(-23);%波尔茨曼常数波尔茨曼常数 NA=6.022*1023;%阿伏加德罗数阿伏加德罗数 m=mu/NA%分子质量分子质量 f=4*pi*(m/2*pi*k*T).(3/2).*v.*v.*exp(-m*v.2./(2*k*T);综合实例
25、综合实例编写主程序编写主程序 ex2.m T=300;mu=28e-3;%给出给出T和和mu的值的值 v=0:1500;%调出自变量数组调出自变量数组 y=mxw(T,mu,v);%调用子程序调用子程序 plot(v,y,r)%绘制分布曲线绘制分布曲线 hold on%为了看出不同的为了看出不同的T和和mu对曲线形状的对曲线形状的影响，再次给定影响，再次给定T和和mu，在同一幅图中，在同一幅图中绘制分布律曲线的图形绘制分布律曲线的图形 T=200;mu=28e-3;y=mxw(T,mu,v);plot(v,y,b)hold onT=300;mu=2e-3;y=mxw(T,mu,v);plot(
26、v,y,g)综合实例综合实例3 3、方波的分解、方波的分解 在连续信号系统中，在连续信号系统中，方波方波可以用可以用相应频率相应频率的的基波及其奇次基波及其奇次谐波合成谐波合成，这也是将，这也是将方波方波展开为正弦级数的出发点。本节将演展开为正弦级数的出发点。本节将演示这一现象。示这一现象。（1 1）建立计算模型建立计算模型 一个以一个以原点为奇对称中心原点为奇对称中心的的方波方波 y(t)可以用奇次正弦波的可以用奇次正弦波的叠加来逼近：叠加来逼近：方波的宽度为方波的宽度为，周期为，周期为 2。.)12sin(121.)5sin(51)3sin(31)sin()(tkktttty2.MATLA
27、B编程编程建立建立 M 文件文件 ex3.m%演示基波和奇次谐波合成方波演示基波和奇次谐波合成方波 t=0:0.1:10;%首先设定一个有首先设定一个有101个点的时间数组个点的时间数组%绘制频率绘制频率w1（f=1/2）的正弦基波，并设置暂停）的正弦基波，并设置暂停 y=sin(t);plot(t,y)pause%叠加叠加3次谐波，绘图并设置暂停次谐波，绘图并设置暂停 y=sin(t)+sin(3*t)/3;plot(t,y)pause%叠加叠加1、3、5、7、9次谐波，绘图并设置暂停次谐波，绘图并设置暂停 y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+
28、sin(9*t)/9;plot(t,y)pause综合实例综合实例%为了绘制三维曲面，需要将各次波形数据存储为一个三维数组，因此需要重新定义为了绘制三维曲面，需要将各次波形数据存储为一个三维数组，因此需要重新定义y，重新编程，本例将求至重新编程，本例将求至 19 次谐波次谐波 t=0:0.031:3.14;y=zeros(10,max(size(t);x=zeros(size(t);for k=1:2:19 x=x+sin(k*t)/k;y(k+1)/2,:)=x;end pause%将各个波形叠合绘出，并设置暂停将各个波形叠合绘出，并设置暂停 plot(t,y(1,:),t,y(2,:),t
29、,y(3,:),t,y(4,:),t,y(5,:),.t,y(7,:),t,y(8,:),t,y(9,:)pause%将各个波形绘制成三维网格图将各个波形绘制成三维网格图 mesh(y)pause综合实例综合实例、程序流控制、程序流控制交互包括两个方面，输入数据和暂停。input命令 格式：n=input（显示信息字符串）通过键盘输入的数据可以是数值也可以是字符串。pause命令：系统暂停状态，也进入键盘主控状态 keyboard 进入键盘主控状态，直接修改或输入变量，用return可退出键盘控制状态，系统继续执行后续操作 echo echo on/off 显示/不显示命令文件中的控制 ech
30、o 两种状态下切换 5.4 函数变量及变量的作用域函数变量及变量的作用域n 控制输入变量的函数（1）nargin：控制输入变量的个数，使得在编程过程中，可以实现不定参数输入的操作（2）varargin：实现不定数目输入变量的函数设计。对函数的一切输入变量均存储在varargin命名的细胞数组中。例：当输入一个矩阵时，求行列式值，当输入二个矩阵时，例：当输入一个矩阵时，求行列式值，当输入二个矩阵时，求矩阵和求矩阵和nFunction c=test69(a,b)if(nargin=1)c=det(a)elseif(nargin=2)c=a+b end5.4 函数变量及变量的作用域函数变量及变量的作
31、用域Function mathavg,engavg,chineseavg=test610(varargin)l=length(varargin)mathsum=0 engsum=0 chinesesum=0 for i=1:l mathsum=mathsum+varargini(1)engsum=engsum+varargini(2)chinesesum=chinesesum+varargini(3)endmathavg=mathsum/l engavg=engsum/l chineseavg=chinesesum/l调用：test610（60,70,80,45,80,90）n控制输出变量的函
32、数控制输出变量的函数lNargout：控制输出变量的个数，使得在编程过程中，可以实现不定参数输入的操作。lVarargout：现不定数目输出变量的函数设计。5.4函数变量及变量的作用域函数变量及变量的作用域全局变量与局部变量全局变量与局部变量 函数M文件的内部变量是局部的，其他函数文件无法使用。为解决函数之间的传递数据的问题，就需要在函数文件执行之前，定义变量为全局变量全局变量，其作用域是整个MATLAB的工作区，即全程有效。命令的格式：global x y zglobal x y z举例：z=(x-1)2+(y+1)2定义过程：函数M文件的内容：function z=fun1(x,y)glo
33、bal alpha beta n=length(y)m=length(x)x1=x*ones(1,n)y1=(y*ones(1,m)z=alpha*(x1-1).2+beta*(y1+1).2 调用过程:global alpha betaglobal alpha betax=0:0.2:2;x=0:0.2:2;y=0:0.2:2;y=0:0.2:2;subplot(2,2,1)subplot(2,2,1)alpha=1alpha=1beta=1beta=1z=fun1(x,y)z=fun1(x,y)mesh(z)mesh(z)改变改变alpha,beta值值,在不同的图形区域下在不同的图形区域
34、下,画画出图形出图形,观察图形的变化情况观察图形的变化情况 5.5子函数与局部函数子函数与局部函数nMATLAB可以定义子函数，扩充函数的功能n在函数文件中题头定义的函数为主函数，在函数体内定义的其他函数均被视为子函数。n把放置在private下的函数称为局部函数，这些函数只有在private目录的父目录中的函数才可以调用。n局部函数和子函数的不同：调用不同，即子函数只能被主函数调用，局部函数只能被父目录下的函数调用。Function c=test(a,b)c=test1(a,b)*test2(a,b)Function c=test1(a,b)c=a+bFunction c=test2(a,b
35、)c=a-b5.6 程序设计的优化程序设计的优化、以矩阵为操作主体，应尽量避免循环运算，转化、以矩阵为操作主体，应尽量避免循环运算，转化为向量运算。为向量运算。function y=test2_5_6(x)x=1 for I=1 to 10000 y(i)=sin(x)x=x+0.1*pi end 若编写为无参数的子程序：function test x=1:0.1:pi*10000 y=sin(x)、数据的与定义、数据的与定义 矩阵或向量在使用时，不一定预先定义，系统可根据使用情况不断扩大或改变矩阵的维数，但这样大大降低了速度，应适当估计变量可能的最大维数，进行预先定义。、内存的管理 使用内存
36、管理命令clear quit save load 等命令。5.7 程序调试 一类是语法错误，包括词法或文法的错误，例如函数 名的拼写错、表达式书写错等。另一类是运行时的错误，指程序的运行结果有错误，这类错误也称为程序逻辑错误。调试器1Debug菜单项 该菜单项用于程序调试，需要与Breakpoints菜单项配合使用。2Breakpoints菜单项 该菜单项共有6个菜单命令，前2个是用于在程序中设置和清除断点的，后4个是设置停止条件的，用于临时停止M文件的执行，并给用户一个检查局部变量的机会，相当于在M文件指定的行号前加入了一个keyboard命令。调试命令 除了采用调试器调试程序外，MATLAB还提供了一些命令用于程序调试。命令的功能和调试器菜单命令类似，具体使用方法请读者查询MATLAB帮助文档。
