第2讲MATLAB程序设计及其调试

M文件的功能和特点MATLAB是解释性语言语法简单调试容易人机交互性强开放性、可扩展性强第2讲MATLAB程序设计及其调试M文件输入窗口！！！例1：通过M脚本文件，画出下列分段函数所表示的曲面.脚本文件first.m%first.m Thisismyfirstexample.a=2;b=2; clf;x=-a:0.2:a;y=-b:0.2:b;fori=1:length(y)forj=1:length(x)ifx(j)+y(i)>1z(i,j)=0.5457*exp(-0.75*y(i)^2-3.75*x(j)^2-1.5*x(j));elseifx(j)+y(i)<=-1z(i,j)=0.5457*exp(-0.75*y(i)^2-3.75*x(j)^2+1.5*x(j));elsez(i,j)=0.7575*exp(-y(i)^2-6.*x(j)^2);endendendaxis([-a,a,-b,b,min(min(z)),max(max(z))]);colormap(flipud(winter));surf(x,y,z);程序调试1程序调试概述一般来说，应用程序的错误有两类，一类是语法错误，另一类是运行时的错误。语法错误包括词法或文法的错误，例如函数名的拼写错、表达式书写错等。程序运行时的错误是指程序的运行结果有错误，这类错误也称为程序逻辑错误。1．控制单步运行2．断点操作

3调试命令除了采用调试器调试程序外，MATLAB还提供了一些命令用于程序调试。命令的功能和调试器菜单命令类似，具体使用方法请读者查询MATLAB帮助文档。例2：某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度(℃)。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];x=0:2.5:10;h=[0:30:60]';T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)例3数据的输入输出处理方法。Matlab环境下的文件与其它系统一样，也有二类文件组成，一是文件，又称M文件，另一类是数据文件.系统除提供了文件的一般管理功能外，还提供了对数据文件进行操作的特殊功能函数.3.1文件的打开与关闭Matlab提供了对数据文件建立、打开、读、写以及关闭等一系列函数，数据文件一般存放在磁盘等介质上，用文件名标识，系统对文件名没有特殊要求。文件数据格式有二种形式，一是二进制格式文件，二是文本文件，系统对这两类文件提供了不同的读写功能函数.文件的打开

fopen函数的调用格式为：

fid=fopen(文件名，打开方式)

其中文件名用字符串形式，表示待打开的数据文件.

常见的打开方式有：

‘r’

表示对打开的文件读数据，

‘w’

表示对打开的文件写数据，

‘a’

表示在打开的文件末尾添加数据.

fid用于存储文件句柄值，句柄值用来标识该数据文件，其他函数可以利用它对该数据文件进行操作.例如：打开一个名为std.dat的数据文件并进行读操作，其命令格式为：Fid=fopen(‘std.dat’,’r’)上述打开格式均为二进制格式，如果想用ASCII文本格式，则必须在格式字符串中加上字符t，例如用’rt’表示以ASCII格式打开供读操作的数据文件.

其他打开方式：‘r’

打开文件，读数据，文件必须存在.‘w’

打开文件，写数据，若文件不存在，系统会自动建立.‘a’

打开文件，在文件末尾添加数据.‘r+’

打开文件，可以读和写数据，文件必须存在.‘w+’

打开文件，供读与写数据用.‘a+’

打开文件，供读与添加数据用.2.文件的关闭fclose关闭文件文件在进行完读、写等操作后，应及时关闭，以保证文件的安全可靠.关闭文件命令格式为：Sta=fclose(Fid)

关闭Fid所表示的文件Sta表示关闭文件操作的返回代码，若关闭成功，返回0，否则返回–1.二进制文件的读写操作

(1)读二进制文件

fread函数可以读取二进制文件的数据，并将数据存入矩阵。其调用格式为：

[A,Count]=fread(fid,size,precision)

其中A用于存放读取的数据，Count返回所读取的数据元素个数，fid为文件句柄，size为可选项，若不选用则读取整个文件内容，若选用则它的值可以是下列值：3.2文件的读写操作

(1)N表示读取N个元素到一个列向量.

(2)Inf表示读取整个文件.

(3)[M,N]表示读数据到M×N的矩阵中，数据按列存放.

precision代表读写数据的类型

如:'int8’:Integer;8bits'int16’:Integer;16bits(2)写二进制文件

fwrite

函数按照指定的数据类型将矩阵中的元素写入到文件中.其调用格式为：

Count=fwrite(fid,A,precision)

其中Count返回所写的数据元素个数，fid为文件句柄，A用来存放写入文件的数据，precision用于控制所写数据的类型，其形式与fread函数相同.

例.建立一数据文件magic5.dat，用于存放5阶魔方阵.

程序如下：

fid=fopen(‘magic5.dat’,‘w’);

cnt=fwrite(fid,magic(5),‘int32’);

%int32:Integer;32bitsfclose(fid);

例.建立一数据文件test.dat，用于存放矩阵A的数据.已知A=[-0.6515-0.2727-0.4354-0.3190-0.9047-0.7534-0.4567-0.3212-0.4132-0.3583-0.9264-0.8173-0.7823-0.3265-0.0631-0.1735-0.7373-0.0972-0.3267-0.6298-0.4768-0.6773-0.6574-0.1923-0.4389]Fid=fopen(‘test.dt’,‘w’)cnt=fwrite(Fid,A,‘float’)fclose(Fid)

程序段将矩阵A的数据以二进制浮点数格式写入文件test.dat中.Fid=fopen(‘test.dt’,‘r’)[B,cnt]=fread(Fid,[5,5],‘float’)fclose(Fid)

读取文件test.dat的内容.[5,5]改成[5,4]，有何不同？

2文本文件fscanf

读ASCII文本文件

[A,COUNT]=fscanf(Fid,format,size)其中A为数据矩阵，用以存放读取的数据，COUNT返回所读取的数据元素个数。format用以控制读取的数据格式，由%加上格式符组成，格式符为：

d,i,o,u,x,e,f,g,s,c与[...]

例如：s=fscanf(fid,‘%s’)读取一个字符串a=fscanf(fid,‘%5d’)读取5位数的整数b=fscanf(fid,‘%6.2d’)读取浮点数fprintf

写ASCII数据文件,其格式为：

COUNT=fprintf(Fid,format,A,…)其中A为要写入文件的数据矩阵，先按format格式化数据矩阵A，后写入到Fid所指定的文件.例如：x=0:0.1:1;y=[x;exp(x)];Fid=fopen('exp.dt','w');fprintf(Fid,'%6.2f%12.8f\n',y);fclose(Fid);1.求极限函数limit用于求符号函数f的极限.系统可以根据用户要求，计算变量从不同方向趋近于指定值的极限值.该函数的格式及功能：

limit(f,x,a)：求符号函数f（x）的极限值.即计算当变量x趋近于常数a时，f（x）函数的极限值.3.3微积分函数

limit(f,a):求符号函数f（x）的极限值.由于没有指定符号函数f（x）的自变量，既变量x趋近于a.

limit(f):求符号函数f（x）的极限值.没有指定变量的目标值时，系统默认变量趋近于0，即a=0的情况.

limit(f,x,a,'right'):求符号函数f的极限值.'right'表示变量x从右边趋近于a.

limit(f,x,a,'left'):求符号函数f的极限值.'left'表示变量x从左边趋近于a.2.符号变量的定义参与符号运算的对象可以是符号变量、符号表达式或符号矩阵.

符号变量要先定义，后引用.可以用sym函数、syms函数将运算量定义为符号型数据.

引用符号运算函数时，用户可以指定函数执行过程中的变量参数；若用户没有指定变量参数，则使用findsym函数默认的变量作为函数的变量参数.2.1sym函数

sym函数的主要功能是创建符号变量，以便进行符号运算，也可以用于创建符号表达式或符号矩阵.用sym函数创建符号变量的一般格式为：

x=sym(‘x’)

其目的是将’x’创建为符号变量，以x作为输出变量名.每次调用该函数,可以定义一个符号变量.例作符号计算：

a,b,x,y均为符号运算量.在符号运算前，应先将a,b,x,y定义为符号运算量.a=sym(‘a’);

%定义‘a’为符号运算量，输出变量名为ab=sym(‘b’);x=sym(‘x’);y=sym(‘y’);[x,y]=solve(a*x-b*y-1,a*x+b*y-5,x,y)%以a,b为符号常数，x,y为符号变量即可得到方程组的解：

x=3/ay=2/b2.2syms函数

syms函数的功能与sym函数类似.syms函数可以在一个语句中同时定义多个符号变量，其一般格式为：

symsarg1arg2…argN

例求极限symsx；

%定义符号变量f=(x*(exp(sin(x))+1)-2*(exp(tan(x))-1))/sin(x)^3；

%确定符号表达式w=limit(f)

%求函数的极限w=-1/23.微分函数diff函数用于对符号表达式s求微分.该函数的一般引用格式为：

diff(s,’v’,n)

应用diff（s）没有指定微分变量和微分阶数，则系统按findsym函数指示的默认变量对符号表达式s求一阶微分.

应用diff（s，‘v’）或diff（s，sym（‘v’））格式，表示以v为自变量，对符号表达式s求一阶微分.

应用diff（s，n）格式，表示对符号表达式s求n阶微分，n为正整数.

应用diff（s，‘v’，n）格式，表示以v为自变量，对符号表达式s求n阶微分.例求导数：x=sym('x');%定义符号变量diff(sin(x^2))%求导运算ans=2*cos(x^2)*x4．积分函数积分函数int（s，v，a，b)可以对被积函数或符号表达式s求积分。其引用格式为：

int（s，v，a，b)说明：应用int（s）格式，表示没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求一阶积分.应用int（s，v）格式，表示以v为自变量，对被积函数或符号表达式s求一阶不定积分.应用积分函数时，如果给定a、b两项，表示是进行定积分运算.a、b分别表示定积分的下限和上限。不指定积分的下限和上限表示求不定积分.例求下述积分.求积分：symsxint(1/(1+x^2))ans=atan(x)解方程函数的格式为：

solve(expr1,expr2,...,exprN,var1,var2,...varN)或

solve(expr1,expr2,...,exprN)其功能为：求解代数方程组expr1,expr2,...,exprN的根，未知数为var1,var2,...varN.说明：若不指明符号表达式expr1,expr2,...,exprN的值，系统默认为0.例如给出一个表达式x^2-3*x-8,则系统将按x^2-3*x-8=0进行运算.5.解方程例解代数方程：a*x2-b*x-6=0symsabxsolve(a*x^2-b*x-6)ans=[1/2/a*(b+(b^2+24*a)^(1/2))][1/2/a*(b-(b^2+24*a)^(1/2))]即该方程有两个根:x1=1/2/a*(b+(b^2+24*a)^(1/2))；

x2=1/2/a*(b-(b^2+24*a)^(1/2))6.级数(级数求和)级数求和运算是数学中常见的一种运算。例如：

f(x)=a0+a1x+a2x2+a3x3+…+anxn函数symsum可以用于此类对符号函数f的求和运算.该函数的引用时，应确定级数的通项式s，变量的变化范围a和b.该函数的引用格式为：

symsum(s,a,b)例求级数的和:

键入：1/12+1/22+1/32+1/42+……symsksymsum(1/k^2,1,Inf)%k值为1到无穷大ans=1/6*pi^2其结果为：1/12+1/22+1/32+1/42+……=π2/6MATLAB求解优化问题的主要函数优化函数的输入变量

优化函数的输出变量

非线性优化有约束的一元函数的最小值单变量函数求最小值的标准形式为

sub.to函数fminbnd格式x=fminbnd(fun,x1,x2)%返回自变量x在区间上函数fun取最小值时x值，fun为目标函数的表达式字符串或MATLAB自定义函数的函数柄.x=fminbnd(fun,x1,x2,options)[x,fval]=fminbnd(…)[x,fval,exitflag]=fminbnd(…)[x,fval,exitflag,output]=fminbnd(…)例4计算下面函数在区间(0，1)内的最小值.

解：>>[x,fval,exitflag,output]=fminbnd('(x^3+cos(x)+x*log(x))/exp(x)',0,1)x=0.5223fval=0.3974exitflag=1output=iterations:9funcCount:9algorithm:'goldensectionsearch,parabolicinterpolation'例5在[0，5]上求下面函数的最小值解：先自定义函数：在MATLAB编辑器中建立M文件为：functionf=myfun(x)f=(x-3).^2-1;保存为myfun.m，然后在命令窗口键入命令：>>x=fminbnd(@myfun,0,5)则结果显示为：x=3无约束多元函数最小值多元函数最小值的标准形式为其中：x为向量.命令利用函数fminsearch求无约束多元函数最小值.函数fminsearch格式x=fminsearch(fun,x0)%x0为初始点，fun为目标函数的表达式字符串或MATLAB自定义函数的函数柄.x=fminsearch(fun,x0,options)%options查optimset.[x,fval]=fminsearch(…)%最优点的函数值.[x,fval,exitflag]=fminsearch(…)%exitflag与单变量情形一致.[x,fval,exitflag,output]=fminsearch(…)%output与单变量情形一致.例6求的最小值点.解：>>X=fminsearch('2*x(1)^3+4*x(1)*x(2)^3-10*x(1)*x(2)+x(2)^2',[0,0])结果为

X=1.00160.8335或在MATLAB编辑器中建立函数文件.functionf=myfun(x)f=2*x(1)^3+4*x(1)*x(2)^3-10*x(1)*x(2)+x(2)^2;保存为myfun.m，在命令窗口键入>>X=fminsearch('myfun',[0,0])或>>X=fminsearch(@myfun,[0,0])结果为：

X=1.00160.8335有约束的多元函数最小值非线性有约束的多元函数的标准形式为：sub.to其中：x、b、beq、lb、ub是向量，A、Aeq为矩阵，C(x)、Ceq(x)是返回向量的函数，f(x)为目标函数，f(x)、C(x)、Ceq(x)可以是非线性函数.在MATLAB5.x中，它的求解由函数constr实现.函数

fmincon格式x=fmincon(fun,x0,A,b)x=fmincon(fun,x0,A,b,Aeq,beq)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)[x,fval]=fmincon(…)[x,fval,exitflag]=fmincon(…)[x,fval,exitflag,output]=fmincon(…)[x,fval,exitflag,output,lambda]=fmincon(…)[x,fval,exitflag,output,lambda,grad]=fmincon(…)[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(…)参数说明：fun为目标函数，它可用前面的方法定义；nonlcon的作用是通过接受的向量x来计算非线性不等约束和等式约束分别在x处的估计C和Ceq，通过指定函数柄来使用，如：>>x=fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)，先建立非线性约束函数，并保存为mycon.m：function[C,Ceq]=mycon(x)C=…%计算x处的非线性不等约束的函数值.Ceq=…%计算x处的非线性等式约束的函数值.lambda是Lagrange乘子，它体现哪一个约束有效.output输出优化信息；grad表示目标函数在x处的梯度；hessian表示目标函数在x处的Hessian值.例7求下面问题在初始点（0，1）处的最优解minsub.to解：约束条件的标准形式为：sub.to先在MATLAB编辑器中建立非线性约束函数文件：

function[c,ceq]=mycon(x)c=(x(1)-1)^2-x(2);ceq=[];%无等式约束.然后，在命令窗口键入如下命令或建立M文件：>>fun='x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2)';%目标函数.>>x0=[01];>>A=[-23];%线性不等式约束.>>b=6;>>Aeq=[];%无线性等式约束.>>beq=[];>>lb=[];%x没有下、上界.>>ub=[];>>[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,@mycon)则结果为x=34fval=-13exitflag=1%解收敛.

output=iterations:2funcCount:9stepsize:1algorithm:'medium-scale:SQP,Quasi-Newton,line-search'firstorderopt:[]cgiterations:[]lambda=lower:[2x1double]%x下界有效情况，通过lambda.lower可查看.upper:[2x1double]%x上界有效情况，为0表示约束无效.

eqlin:[0x1double]%线性等式约束有效情况，不为0表示约束有效.eqnonlin:[0x1double]%非线性等式约束有效情况.ineqlin:2.5081e-008%线性不等式约束有效情况.neqnonlin:6.1938e-008%非线性不等式约束有效情况.grad=%目标函数在最小值点的梯度.1.0e-006*-0.1776hessian=%目标函数在最小值点的Hessian值.1.0000-0.0000-0.00001.0000二次规划问题二次规划问题（quadraticprogramming）的标准形式为：

sub.to

其中，H、A、Aeq为矩阵，f、b、beq、lb、ub、x为向量其它形式的二次规划问题都可转化为标准形式.MATLAB5.x版中的qp函数已被6.0版中的函数quadprog取代。函数

quadprog格式x=quadprog(H,f,A,b)%其中H,f,A,b为标准形中的参数，x为目标函数的最小值.x=quadprog(H,f,A,b,Aeq,beq)%Aeq,beq满足等约束条件.x=quadprog(H,f,A,b,Aeq,beq,lb,ub)%lb,ub分别为解x的下界与上界.x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)%x0为设置的初值x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)%options为指定的优化参数.[x,fval]=quadprog(…)%fval为目标函数最优值.[x,fval,exitflag]=quadprog(…)%exitflag与线性规划中参数意义相同.[x,fval,exitflag,output]=quadprog(…)%output与线性规划中参数意义相同.[x,fval,exitflag,output,lambda]=quadprog(…)%lambda与线性规划中参数意义相同.例8求二次规划的最优解

maxf(x1,x2)=x1x2+3sub.tox1+x2-2=0解：化成标准形式：

sub.tox1+x2=2在Matlab中实现如下：H=[0,-1;-1,0];f=[0;0];Aeq=[11];b=2;[x,fval,exitflag,output,lambda]=quadprog