8.1.1MATLAB数值计算功能
MATLAB强大的数值计算功能使其在诸多数学计算软件中傲视群雄,数值计算是MATLAB软件的基础.本节将简要介绍MATLAB的数据类型、矩阵的建立及运算.
8.1.1.1MATLAB数据类型
MATLAB的数据类型主要包括:数字、字符串、矩阵、单元型数据及结构型数据等,限于篇幅我们将重点介绍其中几个常用类型.
1. 变量与常量
变量是任何程序设计语言的基本要素之一,MATLAB语言当然也不例外.与常规的程序设计语言不同的是MATLAB并不要求事先对所使用的变量进行声明,也不需要指定变量类型,MATLAB语言会自动依据所赋予变量的值或对变量所进行的操作来识别变量的类型.在赋值过程中如果赋值变量已存在时,MATLAB语言将使用新值代替旧值,并以新值类型代替旧值类型.
在MATLAB语言中变量的命名应遵循如下规则:
(1) 变量名区分大小写.
(2) 变量名长度不超31位,第31个字符之后的字符将被MATLAB语言所忽略.
(3) 变量名以字母开头,可以是字母、数字、下划线组成,但不能使用标点.
与其他的程序设计语言相同,在MATLAB语言中也存在变量作用域的问题.在未加特殊说明的情况下,MATLAB语言将所识别的一切变量视为局部变量,即仅在其使用的M文件内有效.若要将变量定义为全局变量,则应当对变量进行说明,即在该变量前加关键字global.一般来说全局变量均用大写的英文字符表示.
MATLAB语言本身也具有一些预定义的变量,这些特殊的变量称为常量.表81给出了MATLAB语言中经常使用的一些常量值.
表81
常量表示数值
pi圆周率
eps浮点运算的相对精度
inf正无穷大
NaN表示不定值
realmax最大的浮点数
i,j虚数单位
在MATLAB语言中,定义变量时应避免与常量名重复,以防改变这些常量的值,如果已改变了某外常量的值,可以通过“clear+常量名”命令恢复该常量的初始设定值(当然,也可通过重新启动MATLAB系统来恢复这些常量值).
2. 数字变量的运算及显示格式
MALAB是以矩阵为基本运算单元的,而构成数值矩阵的基本单元是数字.为了更好地学习和掌握矩阵的运算,首先对数字的基本知识作简单的介绍.
对于简单的数字运算,可以直接在命令窗口中以平常惯用的形式输入,如计算2和3的乘积再加1时,可以直接输入:
1+2*3
ans=
7
这里“ans”是指当前的计算结果,若计算时用户没有对表达式设定变量,系统就自动赋当前结果给“ans”变量.用户也可以输入:
a=1+2*3
a=
7
此时系统就把计算结果赋给指定的变量a了.
MATLAB语言中数值有多种显示形式,在缺省情况下,若数据为整数,则就以整数表示;若数据为实数,则以保留小数点后4位的精度近似表示.MATLAB语言提供了10种数据显示格式,常用的有下述几种格式:
short小数点后4位(系统默认值)
long小数点后14位
shorte 5位指数形式
longe 15位指数形式
MATLAB语言还提供了复数的表达和运算功能.在MATLAB语言中,复数的基本单位表示为i或j.在表达数值时虚部的数值与i、j之间可以不使用乘号,但是如果是表达式,则必须使用乘号以识别虚部符号.
3. 字符串
字符和字符串运算是各种高级语言必不可少的部分,MATLAB中的字符串是其进行符号运算表达式的基本构成单元.
在MATLAB中,字符串和字符数组基本上是等价的;所有的字符串都用单引号进行输入或赋值(当然也可以用函数char来生成).字符串的每个字符(包括空格)都是字符数组的一个元素.例如:
s=??matrixlaboratory??;
s=
matrixlaboratory
size(s)%size查看数组的维数
ans=
117
另外,由于MATLAB对字符串的操作与C语言几乎完全相同这里不在赘述.
8.1.1.2矩阵及其运算
矩阵是MATLAB数据存储的基本单元,而矩阵的运算是MATLAB语言的核心,在MATLAB语言系统中几乎一切运算均是以对矩阵的操作为基础的.下面重点介绍矩阵的生成、矩阵的基本运算和矩阵的数组运算.
1. 矩阵的生成
(1) 直接输入法
从键盘上直接输入矩阵是最方便、最常用的创建数值矩阵的方法,尤其适合较小的简单矩阵.在用此方法创建矩阵时,应当注意以下几点:
输入矩阵时要以“[]”为其标识符号,矩阵的所有元素必须都在括号内.
矩阵同行元素之间由空格或逗号分隔,行与行之间用分号或回车键分隔.
矩阵大小不需要预先定义.
矩阵元素可以是运算表达式.
若“[]”中无元素表示空矩阵.
另外,在MATLAB语言中冒号的作用是最为丰富的.首先,可以用冒号来定义行向量.
例如:
a=1:0.5:4
a=
11.522.533.54
其次,通过使用冒号,可以截取指定矩阵中的部分.
例如:
A=[123;456;789]
A=
123
456
789
B=A(1:2,:)
B=
123
456
通过上例可以看到B是由矩阵A的1到2行和相应的所有列的元素构成的一个新的矩阵.在这里,冒号代替了矩阵A的所有列.
(2) 外部文件读入法
MATLAB语言也允许用户调用在MATLAB环境之外定义的矩阵.可以利用任意的文本编辑器编辑所要使用的矩阵,矩阵元素之间以特定分断符分开,并按行列布置.读入矩阵的一种方法可参考3.3节数据交换系统.另外也可以利用load函数,其调用方法为:Load+文件名[参数]
Load函数将会从文件名所指定的文件中读取数据,并将输入的数据赋给以文件名命名的变量,如果不给定文件名,则将自动认为matlab.mat文件为操作对象,如果该文件在MATLAB搜索路径中不存在时,系统将会报错.
例如:事先在记事本中建立文件(并以data1.txt保存):111
123
136
在MATLAB命令窗口中输入:
load data1.txt
data1
data1=
111
123
136
(3) 特殊矩阵的生成
对于一些比较特殊的矩阵(单位阵、矩阵中含1或0较多),由于其具有特殊的结构,MATLAB提供了一些函数用于生成这些矩阵.常用的有下面几个:
zeros(m)生成m阶全0矩阵
eye(m)生成m阶单位矩阵
ones(m)生成m阶全1矩阵
rand(m)生成m阶均匀分布的随机阵
randn(m)生成m阶正态分布的随机矩阵
2. 矩阵的基本数学运算
矩阵的基本数学运算包括矩阵的四则运算、与常数的运算、逆运算、行列式运算、秩运算、特征值运算等基本函数运算,这里进行简单介绍.
(1) 四则运算
矩阵的加、减、乘运算符分别为“+,—,*”,用法与数字运算几乎相同,但计算时要满足其数学要求(如:同型矩阵才可以加、减).
在MATLAB中矩阵的除法有两种形式:左除“\”和右除“/”.在传统的MATLAB算法中,右除是先计算矩阵的逆再相乘,而左除则不需要计算逆矩阵直接进行除运算.通常右除要快一点,但左除可避免被除矩阵的奇异性所带来的麻烦.在MATLAB6中两者的区别不太大.
(2) 与常数的运算
常数与矩阵的运算即是同该矩阵的每一元素进行运算.但需注意进行数除时,常数通常只能做除数.
(3) 基本函数运算
矩阵的函数运算是矩阵运算中最实用的部分,常用的主要有以下几个:
det(a)求矩阵a的行列式
eig(a)求矩阵a的特征值
inv(a)或a^(-1)求矩阵a的逆矩阵
rank(a)求矩阵a的秩
trace(a)求矩阵a的迹(对角线元素之和)
例如:a=[21-3-1;3107;-124-2;10-15];
a1=det(a);
a2=det(inv(a));
a1*a2
ans=
1
注意:命令行后加“;”表示该命令执行但不显示执行结果.
3. 矩阵的数组运算
我们在进行工程计算时常常遇到矩阵对应元素之间的运算.这种运算不同于前面讲的数学运算,我们称之为数组运算.
(1) 基本数学运算
数组的加、减与矩阵的加、减运算完全相同.而乘除法运算有相当大的区别,数组的乘除法是指两同维数组对应元素之间的乘除法,它们的运算符为“.*”和“./”或“.\”.前面讲过常数与矩阵的除法运算中常数只能做除数.在数组运算中有了“对应关系”的规定,数组与常数之间的除法运算没有任何限制.
另外,矩阵的数组运算中还有幂运算(运算符为.^)、指数运算(exp)、对数运算(log)、和开方运算(sqrt)等.有了“对应元素”的规定,数组的运算实质上就是针对数组内部的每个元素进行的.
例如:
a=[21-3-1;3107;-124-2;10-15];
a^3
ans=
32-28-10134
99-12-151239
-149938
51-17-98139
a.^3
ans=
81-27-1
2710343
-1864-8
10-1125
由上例可见矩阵的幂运算与数组的幂运算有很大的区别.
(2) 逻辑关系运算
逻辑运算是MATLAB中数组运算所特有的一种运算形式,也是几乎所有的高级语言普遍适用的一种运算.它们的具体符号、功能及用法见表82.
表82
符号运算符功能函数名
==等于eq
~=不等于ne
大于gt
=大于等于ge
&逻辑与and
|逻辑或or
~逻辑非not
说明:
?? 在关系比较中,若比较的双方为同维数组,则比较的结果也是同维数组.它的元素值由0和1组成.当比较双方对应位置上的元素值满足比较关系时,它的对应值为1,否则为0.
?? 当比较的双方中一方为常数,另一方为一数组,则比较的结果与数组同维.
?? 在算术运算、比较运算和逻辑与、或、非运算中,它们的优先级关系先后为:
比较运算、算术运算、逻辑与或非运算.
例如:
a=[123;456;789];
x=5;
y=ones(3)*5;
xa=xb=[010;101;001];
ab=a&b
ab=
010
101
001
8.1.2MATLAB图形功能
MATLAB有很强的图形功能,可以方便地实现数据的视觉化.强大的计算功能与图形功能相结合为MATLAB在科学技术和教学方面的应用提供了更加广阔的天地.下面着重介绍二维图形的画法,对三维图形只作简单叙述.
8.1.2.1二维图形的绘制
1. 基本形式
二维图形的绘制是MATLAB语言图形处理的基础,MATLAB最常用的画二维图形的命令是plot.看两个简单的例子:
y=[0.580.700.950.830.25];
plot(y)
生成的图形见图81,是以序号1,2,…,6为横坐标、数组y的数值为纵坐标画出的折线.
x=linspace(0,2*pi,30);%生成一组线性等距的数值
y=sin(x);
plot(x,y)
生成的图形见图82,是[0,2π]上30个点连成的光滑的正弦曲线.
图81
图82
2. 多重线
在同一个画面上可以画许多条曲线,只需多给出几个数组,例如
x=0:pi/15:2*pi;
y1=sin(x);
y2=cos(x);
plot(x,y1,x,y2)
则可以画出图83.多重线的另一种画法是利用hold命令.在已经画好的图形上,若设置hold on,MATLA将把新的plot命令产生的图形画在原来的图形上.而命令hold off将结束这个过程.例如:
x=linspace(0,2*pi,30);y=sin(x);plot(x,y)
先画好图82,然后用下述命令增加cos(x)的图形,也可得到图83.
hold on
z=cos(x);plot(x,z)
hold off
图83
3. 线型和颜色
MATLAB对曲线的线型和颜色有许多选择,标注的方法是在每一对数组后加一个字符串参数,说明如下:
线型线方式:-实线:点线-.虚点线--波折线.
线型点方式:.圆点+加号*星号xx形o小圆
颜色:y黄;r红;g绿;b蓝;w白;k黑;m紫;c青.
以下面的例子说明用法:
x=0:pi/15:2*pi;
y1=sin(x);y2=cos(x);
plot(x,y1,??b:+??,x,y2,??g-.*??)
可得图形84.
图84
4. 网格和标记
在一个图形上可以加网格、标题、x轴标记、y轴标记,用下列命令完成这些工作.
x=linspace(0,2*pi,30);y=sin(x);z=cos(x);
plot(x,y,x,z)
grid
xlabel(??IndependentVariableX??)
ylabel(??DependentVariablesYandZ??)
title(??SineandCosineCurves??)
可得到图85:
图85
也可以在图形的任何位置加上一个字符串,如用:
text(2.5,0.7,??sinx??)
表示在坐标x=2.5,y=0.7处加上字符串sinx.更方便的是用鼠标来确定字符串的位置,方法是输入命令:
gtext(??sinx??)
在图形窗口十字线的交点是字符串的位置,用鼠标点一下就可以将字符串放在那里.
5. 坐标系的控制
在缺省情况下MATLAB自动选择图形的横、纵坐标的比例,如果你对这个比例不满意,可以用axis命令控制,常用的有:
axis([xmin xmax ymin ymax])[]中分别给出x轴和y轴的最大值、最小值
axisequal或axis(??equal??)x轴和y轴的单位长度相同
axissquare或axis(??square??)图框呈方形
axisoff或axis(??off??)清除坐标刻度
还有axisauto, axisimage, axisxy, axisij, axisnormal, axison, axis(axis)
用法可参考在线帮助系统.
6. 多幅图形
可以在同一个画面上建立几个坐标系,用subplot(m,n,p)命令;把一个画面分成m×n个图形区域,p代表当前的区域号,在每个区域中分别画一个图,如
x=linspace(0,2*pi,30);y=sin(x);z=cos(x);
u=2*sin(x).*cos(x);v=sin(x)./cos(x);
subplot(2,2,1),plot(x,y),axis([02*pi-11]),title(??sin(x)??)
subplot(2,2,2),plot(x,z),axis([02*pi-11]),title(??cos(x)??)
subplot(2,2,3),plot(x,u),axis([02*pi-11]),title(??2sin(x)cos(x)??)
subplot(2,2,4),plot(x,v),axis([02*pi-2020]),title(??sin(x)/cos(x)??)
共得到4幅图形,见图86.
图86
8.1.2.2三维图形
限于篇幅这里只对几种常用的命令通过例子作简单介绍.
1. 带网格的曲面
作曲面z=f(x,y)的图形
z=sinx2+y2x2+y2,-7.5≤x≤7.5,-7.5≤y≤7.5
用以下程序实现:
x=-7.5:0.5;7.5;
y=x;
[X,Y]=meshgrid(x,y);(3维图形的X,Y数组)
R=sqrt(X.^2+Y.^2)+eps;(加eps是防止出现0/0)
Z=sin(R)./R;
mesh(X,Y,Z)(3维网格表面)
画出的图形如图87.mesh命令也可以改为surf,只是图形效果有所不同,读者可以上机查看结果.
图87
2. 空间曲线
作螺旋线x=sint,y=cost,z=t
用以下程序实现:
t=0:pi/50:10*pi;
plot3(sin(t),cos(t),t)(空间曲线作图函数,用法类似于plot)
画出的图形如图88
图88
3. 等高线
用contour或contour3画曲面的等高线,如对于图87的曲面,在上面的程序后接contour(X,Y,Z,10)即可得到10条等高线.
4. 其他
较有用的是给三维图形指定观察点的命令view(azi,ele),azi是方位角,ele是仰角.缺省时azi=-37.5°,ele=30°.
8.1.2.3图形的输出
在数学建模中,往往需要将产生的图形输出到Word文档中.通常可采用下述方法:
首先,在MATLAB图形窗口中选择【File】菜单中的【Export】选项,将打开图形输出对话框,在该对话框中可以把图形以emf、bmp、jpg、pgm等格式保存.然后,再打开相应的文档,并在该文档中选择【插入】菜单中的【图片】选项插入相应的图片即可.
8.1.3程序设计
MATLAB作为一种高级语言,它不仅可以如前几节所介绍的那样,以一种人机交互式的命令行的方式工作,还可以像BASIC、FORTRAN、C等其他高级计算机语言一样进行控制流的程序设计,即编制一种以.m为扩展名的MATLAB程序(简称M文件).而且,由于MATLAB本身的一些特点,M文件的编制同上述几种高级语言比较起来,有许多无法比拟的优点.
8.1.3.1M文件
所谓M文件就是由MATLAB语言编写的可在MATLAB语言环境下运行程序源代码文件.由于商用的MATLAB软件是用C语言编写而成.因此,M文件的语法与C语言十分相似.对广大参加建模竞赛且学过C语言的同学来说,M文件的编写是相当容易的.M文件可以分为脚本文件(Script)和函数文件(Function)两种.M文件不仅可以在MATLAB的程序编辑器中编写,也可以在其他的文本编辑器中编写,并以“m”为扩展名加以存储.
1. 脚本文件
脚本类似于DOS下的批处理文件,不需要在其中输入参数,也不需要给出输出变量来接受处理结果,脚本仅是若干命令或函数的集合,用于执行特定的功能.脚本的操作对象为MATLAB工作空间内的变量,并且在脚本执行结束后,脚本中对变量的一切操作均会被保留.在MATLAB语言中也可以在脚本内部定义变量,并且该变量将会自动地被加入到当前的MATLAB工作空间中,并可以为其他的脚本或函数引用,直到MATLAB被关闭或采用一定的命令将其删除.
例如:
%命令窗口中定义矩阵a,b
a=pascal(3)
a=
111
123
136
b=magic(3)
b=
816
357
492
%在编辑器中编写下述命令
a=a+b
b=a-b
a=a-b
在编辑器中编辑完上例的脚本文件后,保存至文件scripts—example中,然后在工作窗口中调用该脚本文件,
scripts—example
a
a=
816
357
492
b
b=
111
123
136
其中矩阵a、b均是在工作空间中已定义完毕的,脚本运行时直接使用该变量,并对其进行操作,然后在命令窗口中调用该脚本,可以看到变量a、b已经进行了相互交换.
2. 函数文件
MATLAB语言中,相对于脚本文件而言,函数文件是较为复杂的.函数需要给定输入参数,并能够对输入变量进行若干操作,实现特定的功能,最后给出一定的输出结果或图形等,其操作对象为函数的输入变量和函数内的局部变量等.
MATLAB语言的函数文件包含如下5个部分.
函数题头:指函数的定义行,是函数语句的第一行,在该行中将定义函数名、输入变量列表及输出变量列表等.
HI行:指函数帮助文本的第一行,为该函数文件的帮助主题,当使用lookfor命令时,可以查看到该行信息.
帮助信息:这部分提供了函数的完整的帮助信息,包括HI之后至第一个可执行行或空行为止的所有注释语句,通过MATLAB语言的帮助系统查看函数的帮助信息时,将显示该部分.
函数体;指函数代码段,也是函数的主体部分.
注释部分:指对函数体中各语句的解释和说明文本,注释语句是以%引导的.
例如:
function[output,output2]=function—example(input1,input2)%函数题头
%This is function to exchange two matrices%HI行
%input1,input2 are inputvariables%帮助信息
%output1,output2 are outputvariables%帮助信息
output1=input2;%函数体
output2=input1;%函数体
%Theend of this example function
[a,b]=function—example(a,b)
a=
816
357
492
b=
111
123
136
可以看到通过使用函数可以和上一节中的示例一样对矩阵a、b进行了相互交换.在该函数题头中,“function”为MATLAB语言中函数的标示符,而function??example为函数名,input1、input2为输入变量,而output1、output2为输出变量,实际调用过程中,可以用有意义的变量替代使用.题头的定义是有一定的格式要求的,输出变量是由中括号标识的,而输入变量是由小括号标识的,各变量间用逗号间隔,应该注意到,函数的输入变量引用的只是该变量的值而非其他值,所以函数内部对输入变量的操作不会带回到工作空间中.
函数题头下的第一行注释语句为HI行,可以通过lookfor命令查看;函数的帮助信息可以通过help命令查看.
函数体是函数的主体部分,也是实现编程目的的核心所在,它包括所有可执行的一切MATLAB语言代码.
在函数体中“%”后的部分为注释语句,注释语句主要是对程序代码进行说明解释,使程序易于理解,也有利于程序的维护.MATLAB语言中将一行内百分号后所有文本均视为注释部分,在程序的执行过程中不被解释,并且百分号出现的位置也没有明确的规定,可以是一行的首位,这样,整行文本均为注释语句,也可以是在行中的某个位置,这样其后所有文本将被视为注释语句,这也展示了MATLAB语言在编程中的灵活性.
尽管在上文中介绍了函数文件的5个组成部分,但是并不是所有的函数文件都需要全部的这5个部分,实际上,这5部分中只有函数题头是一个函数文件所必需的,而其他的4个部分均可省略.当然,如果没有函数体则为一空函数,不能产生任何作用.
在MATLAB语言中,存储M文件时文件名应当与文件内主函数名相一致,这是因为在调用M文件时,系统查询的相应的文件而不是函数名,如果两者不一致,则或者打不开目的文件,或者打开的是其他文件.鉴于这种查询文件的方式与以往程序设计语言不同,在其他的语言系统中,函数的调用都是指对函数名本身的,所以,建议在存储M文件时,应将文件名与主函数名统一起来,以便于理解和使用.
8.1.3.2函数变量及变量作用域
在MATLAB语言的函数中,变量主要有输入变量、输出变量及函数内所使用的变量.输入变量相当于函数入口数据,是一个函数操作的主要对象.某种程度上讲,函数的作用就是对输入变量进行加工以实现一定的功能.如前节所述,函数的输入变量为形式参数,即只传递变量的值而不传递变量的地址,函数对输入变量的一切操作和修改如果不依靠输出变量传出的话,将不会影响工作空间中该变量的值.
MATLAB语言提供了函数nargin和函数varargin来控制输入变量的个数,以实现不定个数参数输入的操作.
函数对于函数变量而言,还应当指出的是其作用域的问题.在MATLAB语言中,函数内定义的变量均被视为局部变量,即不加载到工作空间中,如果希望使用全局变量,则应当使用命令global定义,而且在任何使用该全局变量的函数中都应加以定义.在命令窗口中也不例外.
例如:
%这里一个全局变量的示例
function[num1,num2,num3]=text(varargin)
global firstlevel secondlevel%定义全局变量
num1=0;
num2=0;
num3=0;
list=zeros(nargin);
fori=1:nargin
list(i)=sum(varargin{i}(:));
list(i)=list(i)/length(varargin{i});
iflist(i)>firstlevel
num1=num1+1
elseiflist(i)>secondlevel
num2=num2+1;
else
num3=num3+1;
end
end
%在命令窗口中也应定义相应的全局变量
global firstlevel secondlevel
firstlevel=85;
secondlevel=75;(程序运行结果略)
从该例中可以看到,定义全局变量时,与定义输入变量和输出变量不同,变量之间必须用空格分隔,而不能用逗号分隔,否则系统将不能识别逗号后的全局变量.
8.1.3.3子函数与局部函数
在MATLAB语言中,与其他的程序设计语言类似,也可以定义子函数,以扩充函数的功能.在函数文件中题头中所定义的函数为主函数,而在函数体内定义的其他函数均被视为子函数.子函数只能被主函数或同一主函数下其他的子函数所调用.
在MATLAB语言中将放置在目录private下的函数称为局部函数,这些函数只能被private目录的父目录中函数调用,而不能被其他的目录的函数调用.
局部函数与子函数所不同的是局部函数可以被其父目录下的所有函数所调用,而子函数则只能为其所在的M文件的主函数所调用,所以局部函数可应用范围大于子函数;在函数编辑的结构上,局部函数与一般的函数文件的编辑相同,而子函数则只能在主函数文件中编辑.
当在MATLAB的M文件中调用函数时,首先将检测该函数是否为此文件的子函数;如果不是的话,再检测是否为可用的局部函数;当结果仍然为否定时,再检测该函数是否为MATLAB搜索路径上的其他M文件.
8.1.3.4流程控制语句
如其他的程序设计语言一样,MATLAB语言也给出了丰富的流程控制语句,以实现具体的程序设计.在命令窗口中的操作虽然可以实现人面交互,但是所能实现的功能却相对简单,虽然也可以在命令窗口中使用流程控制语句,但是由于命令窗口中交互式的执行方式,使用这样的操作极为不方便;而在M文件中,通过对流程控制语句的组合使用,可以实现多种复杂功能.MATLAB语言的流程控制语句主要有for、while、if??else??end及switch??case等4种语句.
1. for语句
for循环语句是流程控制语句中的基础,使用该循环语句可以以指定的次数重复执行循环体内的语句.
for循环语句的调用形式为:
for循环控制变量=〈循环次数设定〉
循环体
end
例如:
for i=1∶2∶12
s=s+i;
end
在上例中,循环次数由数组1∶2∶12决定,设定循环次数的数组可以是已定义的数组,也可以在for循环语句中定义,此时定义的格式为:
〈初始值〉:〈步长〉:〈终值〉
初始值为循环变量的初始设定值,每执行循环体一次,循环控制变量将增加步长大小,直至循环控制变量的值大于终值时循环结束,这里步长是可以为负的.在for循环语句中,循环体内不能出现对循环控制变量的重新设置,否则将会出错,for循环允许嵌套使用.
2. while语句
while循环语句与for循环语句不同的是,前者是以条件的满足与否来判断循环是否结束的,而后者则是以执行次数是否达到指定值为判断的.
while循环语句的一般形式为:
while〈循环判断的语句〉
循环体
end
其中循环判断语句为某种形式的逻辑判断表达式,当该表达式的值为真时,就执行循环体内的语句;当表达式的逻辑值为假时,就退出当前的循环体.如果循环判断语句为矩阵时,当且仅当所有的矩阵元素非零时,逻辑表达式的值为真.
在while循环语句中,在语句内必须有可以修改循环控制变量的命令,否则该循环语言将陷入死循环中,除非循环语句中有控制退出循环的命令,如break语句.当程序流程运行至该命令时,则不论循环控制变量是否满足循环判断语句均将退出当前循环,执行循环后的其他语句.
与break语句对应,MATLAB还提供了continue命令用于控制循环,当程序流运行至该命令时会忽略其后的循环体操作转而执行下一层次的循环.当循环控制语句为一空矩阵时,将不执行循环体的操作而直接执行其后的其他命令语句,即空矩阵被认为是假.
3. if??else??end语句
条件判断语句也是程序设计语言中流程控制语句之一.使用该语句,可以选择执行指定的命令,MATLAB语言中的条件判断语句是if??else??end语句.
if??else??end语句的一般形式为:
if〈逻辑判断语句〉
逻辑值为“真”时执行的语句
else
逻辑值为“假”时执行的语句
end
当逻辑判断表达式为“真”时,将执行if与else语句间的命令,否则将执行else与end语句间的命令.
例如:
if a=1
a=a+1
else
a=a+2
end
在MATLAB语言的if??else??end语句中的eles子句是可选项,即语句中可以不包括else子句的条件判断.在程序设计中,也经常碰到需要进行多重逻辑选择的问题,这时可以采用if??else??end语句的嵌套形式:
if〈逻辑判断语句1〉
逻辑值1为“真”时的执行语句
elseif〈逻辑判断语句2〉
逻辑值2“真”时的执行语句
elseif〈逻辑判断语句3〉
……
else
当以上所有的逻辑值均为假时的执行语句
end
在以上的各层次的逻辑判断中,若其中任意一层逻辑判断为真,则将执行对应的执行语句,并跳出该条件判断语句,其后的逻辑判断语句均不进行检查.
4. switch??case语句
if??else??end语句所对应的是多重判断选择,而有时也会遇到多分支判断选择的问题.MATLAB语言为解决多分支判断选择提供了switch??case语句.
switch??case语句的一般表达形式为:
switch〈选择判断量〉
case选择判断值1
选择判断语句1
case选择判断值2
选择判断语句2
……
otherwise
判断执行语句
end
与其他的程序设计语言的switch??case语句不同的是,在MATLAB语言中,当其中一个case语句后的条件为真时,switch??case语句不对其后的case语句进行判断,也就是说在MATLAB语言中,即使有多条case判断语句为真,也只执行所遇到的第一条为真的语句.这样就不必像C语言那样,在每条case语句后加上break语句以防止继续执行后面为真的case条件语句.
8.1.4MATLAB的应用
8.1.4.1MATLAB在数值分析中的应用
插值与拟合是来源于实际、又广泛应用于实际的两种重要方法.随着计算机的不断发展及计算水平的不断提高,它们已在国民生产和科学研究等方面扮演着越来越重要的角色.下面对插值中分段线性插值、拟合中的最为重要的最小二乘法拟合加以介绍.
1. 分段线性插值
所谓分段线性插值就是通过插值点用折线段连接起来逼近原曲线,这也是计算机绘制图形的基本原理.实现分段线性插值不需编制函数程序,MATLAB自身提供了内部函数interp1其主要用法如下:
interp1(x,y,xi)一维插值
◆ yi=interp1(x,y,xi)
对一组点(x,y)进行插值,计算插值点xi的函数值.x为节点向量值,y为对应的节点函数值.如果y为矩阵,则插值对y的每一列进行,若y的维数超出x或xi的维数,则返回NaN.
◆ yi=interp1(y,xi)
此格式默认x=1:n,n为向量y的元素个数值,或等于矩阵y的size(y,1).
◆ yi=interp1(x,y,xi,??method??)
method用来指定插值的算法.默认为线性算法.其值常用的可以是如下的字符串.
?? nearest线性最近项插值.
?? linear线性插值.
?? spline三次样条插值.
?? cubic三次插值.
所有的插值方法要求x是单调的.x也可能并非连续等距的.
正弦曲线的插值示例:
x=0:0.1:10;
y=sin(x);
xi=0:0.25:10;
yi=interp1(x,y,xi);
plot(x,y,??0??,xi,yi)
则可以得到相应的插值曲线(读者可自己上机实验).
Matlab也能够完成二维插值的运算,相应的函数为interp2,使用方法与interpl基本相同,只是输入和输出的参数为矩阵,对应于二维平面上的数据点,详细的用法见Matlab联机帮助.
2. 最小二乘法拟合
在科学实验的统计方法研究中,往往要从一组实验数据(xi,yi)中寻找出自变量x和因变量y之间的函数关系y=f(x).由于观测数据往往不够准确,因此并不要求y=f(x)经过所有的点(xi,yi),而只要求在给定点xi上误差δi=f(xi)-yi按照某种标准达到最小,通常采用欧氏范数‖δ‖2作为误差量度的标准.这就是所谓的最小二乘法.在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行.
函数polyfit是指用一个多项式函数来对已知数据进行拟合,我们以下列数据为例介绍这个函数的用法:
x=0:0.1:1;
y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2]
为了使用polyfit,首先必须指定我们希望以多少阶多项式对以上数据进行拟合,如果我们指定一阶多项式,结果为线性近似,通常称为线性回归.我们选择二阶多项式进行拟合.
P=polyfit(x,y,2)
P=
-9.810820.1293-0.0317
函数返回的是一个多项式系数的行向量,写成多项式形式为:
-9.8108x2+20.1293x-0.0317
为了比较拟合结果,我们绘制两者的图形:
xi=linspace(0,1,100);%绘图的X-轴数据.
Z=polyval(p,xi);%得到多项式在数据点处的值.
当然,我们也可以选择更高幂次的多项式进行拟合,如10阶:
p=polyfit(x,y,10);
xi=linspace(0,1,100);
z=ployval(p,xi);
读者可以上机绘图进行比较,曲线在数据点附近更加接近数据点的测量值了,但从整体上来说,曲线波动比较大,并不一定适合实际使用的需要,所以在进行高阶曲线拟合时,“越高越好”的观点不一定是对的.
8.1.4.2符号工具箱及其应用
在数学应用中,常常需要做极限、微分、求导等运算,MATLAB称这些运算为符号运算.MATLAB的符号运算功能是通过调用符号运算工具箱(SymbolicMathToolbox)内的工具实现,其内核是借用Maple数学软件的.MATLAB的符号运算工具箱包含了微积分运算、化简和代换、解方程等几个方面的工具,其详细内容可通过MATLAB系统的联机帮助查阅,本节仅对它的常用功能做简单介绍.
1. 符号变量与符号表达式
MATLAB符号运算工具箱处理的对象主要是符号变量与符号表达式.要实现其符号运算,首先需要将处理对象定义为符号变量或符号表达式,其定义格式如下:
格式1:sym(‘变量名’)或sym(‘表达式’)
功能:定义一个符号变量或符号表达式.
例如:
sym(??x??)%定义变量x为符号变量
sym(??x+1??)%定义表达式x+1为符号表达式
格式2:syms变量名1变量名2……变量名n
功能:定义变量名1、变量2……、变量名n为符号变量.
例如:
symsabxt%定义a,b,x,t均为符号变量
2. 微积分运算
(1) 极限
格式:limit(f,t,a,??left?? or ??right??)
功能:求符号变量t趋近a时,函数f的(左或右)极限.??left??表示求左极限,??right??表示求右极限,省略时表示求一般极限;a省略时变量t趋近0;t省略时默认变量为x,若无x则寻找(字母表上)最接近字母x的变量.
例如:求极限limx→∞1+2tx3x的命令及结果为:
symsxt
limit((1+2*t/x)^(3*x),x,inf)
ans=
exp(6*t)
再如求函数x/|x|,当x→0时的左极限和右极限,命令及结果为:
symsx
limit(x/abs(x),x,0,??left??)ans=-1
limit(x/abs(x),x,0,??right??)ans=1
(2) 导数
格式:diff(f,t,n)
功能:求函数f对变量t的n阶导数.当n省略时,默认n=1;当t省略时,默认变量x,若无x时则查找字母表上最接近字母x的字母.
例如:求函数f=a*x^2+b*x+c对变量x的一阶导数,命令及结果为
symsabcx
f=a*x^2+b*x+c;
diff(f)
ans=
2*a*x+b
求函数f对变量b的一阶导数(可看作求偏导),命令及结果为
diff(f,b)
ans=x
求函数f对变量x的二阶导数,命令及结果为
diff(f,2)
ans=2*a
(3) 积分
格式:int(f,t,a,b)
功能:求函数f对变量t从a到b的定积分.当a和b省略时求不定积分;当t省略时,默认变量为(字母表上)最接近字母x的变量.
例如:求函数f=a*x^2+b*x+c对变量x不定积分,命令及结果为
symsabcx
f=a*x^2+b*x+c;
int(f)
ans=
1/3*a*x^3+1/2*b*x^2+c*x
求函数f对变量b不定积分,命令及结果为
int(f,b)
ans=
a*x^2*b+1/2*b^2*x+c*b
求函数f对变量x从1到5的定积分,命令及结果为
int(f,1,5)
ans=
124/3*a+12*b+4*c
(4) 级数求和
格式:symsum(s,t,a,b)
功能:求表达式s中的符号变量t从第a项到第b项的级数和.
例如:求级数1/1+1/2+1/3+…+1/x的前三项的和,命令及结果为
symsum(1/x,1,3)ans=11/6
3. 化简和代换
MATLAB符号运算工具箱中,包括了较多的代数式化简和代换功能,下面仅举出部分常见运算.
simplify利用各种恒等式化简代数式
expand将乘积展开为和式
factor把多项式转换为乘积形式
collect合并同类项
horner把多项式转换为嵌套表示形式
例如:进行合并同类项执行
symsx
collect(3*x^3-0.5*x^3+3*x^2)
ans=
5/2*x^3+3*x^2)
进行因式分解执行
factor(3*x^3-0.5*x^3+3*x^2)
ans=
1/2*x^2*(5*x+6)
4. 解方程
(1) 代数方程
格式:solve(f,t)
功能:对变量t解方程f=0,t缺省时默认为x或最接近字母x的符号变量.
例如:求解一元二次方程f=a*x^2+b*x+c的实根,
symsabcx
f=a*x^2+b*x+c;
solve(f,x)
ans=
[1/2/a*(-b+(b^2-4*a*c)^(1/2))]
[1/2/a*(-b-(b^2-4*a*c)^(1/2))]
(2) 微分方程
格式:dsolve(??s??,??s1??,??s2??,…,??x??)
其中s为方程;s1,s2,……为初始条件,缺省时给出含任意常数c1,c2,……的通解;x为自变量,缺省时默认为t.
例如:求微分方程y′=1+y2的通解
dsolve(??Dy=1+y^2??)
ans=
tan(t+c1)
8.1.4.3优化工具箱及其应用
在工程设计、经济管理和科学研究等诸多领域中,人们常常会遇到这样的问题:如何从一切可能的方案中选择最好、最优的方案,在数学上把这类问题称为最优化问题.这类问题很多,例如当设计一个机械零件时如何在保证强度的前提下使重量最轻或用量最省(当然偷工减料除外);如何确定参数,使其承载能力最高;在安排生产时,如何在现有的人力、设备的条件下,合理安排生产,使其产品的总产值最高;在确定库存时如何在保证销售量的前提下,使库存成本最小;在物资调配时,如何组织运输使运输费用最少.这些都属于最优化问题所研究的对象.
MATLAB的优化工具箱被放在toolbox目录下的optim子目录中,其中包括有若干个常用的求解函数最优化问题的程序.MATLAB的优化工具箱也在不断地完善.不同版本的MATLAB,其工具箱不完全相同.在MATLAB5.3版本中,对优化工具箱作了全面的改进.每个原有的常用程序都重新编制了一个新的程序.除fzero和fsolve外都重新起了名字.这些新程序使用一套新的控制算法的选项.与原有的程序相比,新程序的功能增强了.在MATLAB5.3和6.0版本中,原有的优化程序(除fzero和fsolve外)仍然保留并且可以使用,但是它们迟早会被撤消的.鉴于上述情况,本书将只介绍那些新的常用的几个优化程序.
1. 线性规划问题
线性规划是最优化理论发展最成熟,应用最广泛的一个分支.在MATLAB的优化工具箱中用于求解下述线性规划的问题
minz=cx
s.t.Ax≤b(线性不等式约束)
A1x=b1(线性等式约束)
LB≤x≤UB(有界约束)
的函数是linprog,其主要格式为:
[x,fval,exitflag,output,lambda]=linprog(c,A,b,A1,b1,LB,UB,x0,options)
其中,linprog为函数名,中括号及小括号中所含的参数都是输入或输出变量,这些参数的主要用法及说明如下:
(1) c,A和b是不可缺省的输入变量;x是不可缺省的输出变量,它是问题的解.
(2) 当x无下界时,在LB处放置[].当无上界时,在UB处放置[].如果x的某个分量xi无下界,则置LB(i)=-inf.如果xi无上界,则置UB(i)=inf.如果无线性不等式约束,则在A和b处都放置[].
(3) x0是解的初始近似值.
(4) options是用来控制算法的选项参数向量.
(5) 输出变量fval是目标函数在解x处的值.
(6) 输出变量exitflag的值描述了程序的运行情况.如果exitflag的值大于0,则程序收敛于解x;如果exitflag的值等于0,则函数的计算达到了最大次数;如果exitflag的值小于0,则问题无可行解,或程序运行失败.
(7) 输出变量output输出程序运行的某些信息.
(8) 输出变量Lambda为在解x处的值Lagrange乘子.
例:求解线性规划问题
minz=-2x1-x2+x3,
s.t.x1+x2+2x3=6,
x1+4x2-x3≤4,
2x1-2x2+x3≤12,
x1≥0,x2≥0,x3≤5.
解:在命令窗口中键入
c=[-2,-1,1];a=[1,4,-1;2,-2,1];b=[4;12];a1=[1,1,2];b1=6;
lb=[0;0;-inf];ub=[inf;inf;5];
[x,z]=linprog(c,a,b,a1,b1,1b,ub)
运行后得到:
x=
4.6667
0.0000
0.6667
z=
-8.6667
2. 非线性约束最优化
在MATLAB的优化工具箱中有一个求解下述非线性规划的问题
minf(x)
s.t.Ax≤b(线性不等式约束)
A1x=b1(线性等式约束)
C(x)≤0(非线性不等式约束)
C1(x)=0(非线性等式约束)
LB≤x≤UB(有界约束)
的函数是fmincon,其主要格式为:
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(??fun??,x0,A,b,A1,b1,LB,
UB,??nonlcon??,options,p1,p2,……)
其中,fmincon为函数名,参数的主要用法与线性规划中的相同,下面介绍几个非线性规划特有的参数:
(1) ??fun??和x0是不可缺省的输入变量.fun是给出目标函数的M文件的名字,x0是极小值点的初始近似值.x是不可缺省的输出变量,它是问题的解.
(2) nonlcon是给出非线性约束函数C(x)和C1(x)的M文件的文件名.
(3) 变量p1,p2…是向目标函数传送的参数的值.
(4) 输出变量grad为目标函数在解x处的梯度.
(5) 输出变量hessian为目标函数在解x处的Hessian矩阵.
例:求解非线性规划问题
minf(x)=ex1(4x21+2x22+4x1x2+2x2+1),
s.t.x1-x2≤1,
x1+x2=0,
1.5+x1x2-x1-x2≤0,
-x1x2-10≤0
解:建立目标函数的M文件
functiony=nline(x)
y=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
建立非线性约束条件的M文件
function[c1,c2]=nyueshu(x)
c1=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];
c2=0;
在命令窗口中键入
x0=[-1,1];a=[1,-1];b=1;a1=[1,1];b1=0;
[x,f]=fmincon(??nline??,x0,a,b,a1,b1,[],[],??nyueshu??)
运行后得到:
x=
-1.22471.2247
f=
1.8951
3. 二次规划问题
二次规划数学模型的一般形式为:
min12xTHx+cx
s.t.Ax≤b
A1x=b1
LB≤x≤UB
其中H为对称矩阵,约束条件与线型规划相同.在MATLAB的优化工具箱中有一个求解上述规划问题的程序:
[x,fval,exitflag,output,lambda]=quadprog(H,c,A,b,A1,b1,LB,UB,x0,options)
其中,quadprog为函数名,参数的主要用法及说明同线性规划,这里不再赘述.
例求解如下二次优化问题.
minf(x)=x21+x22-8x1-10x2
s.t.3x1+2x2≤6
x1,x2≥0
解:将目标函数化为标准形式
f(x)=12(x1x2)20
02
x1
x2+(-8-10)x1
x2θ
在命令窗口中键入
H=[2,0;0,2];c=[-8,-10];a=[3,2];b=6;lb=[0,0];x0=[1,1];
x=quadprog(H,c,a,b,[],[],lb,[],x0)
运行后得到:
x=
0.3077
2.5385
4. foptions函数
对于优化的控制,MATLAB共提供了18个参数,这些参数对优化的进行起者很关键的作用.下面就对参数选择函数foptions作详细介绍.
?? foptions优化函数调用中的参数选择.参数具体意义如下:
options(1)参数显示控制(默认值为0).等于1时显示一些结果.
options(2)优化点x的精度控制(默认值为1e4).
options(3)优化函数F的精度控制(默认值为1e4).
options(4)违反约束的结束标准(默认值为1e6).
options(5)策略选择.不常用.
options(6)优化程序方法的选择.值为0时为BFGS算法,值为1时采用DFP算法.
options(7)线性插值算法选择.值为0时为混合插值算法,值为1时采用立方插值算法.
options(8)函数值显示(目标-达到问题中的Lambda).
options(9)若需要检测用户提供的导数则设为1.
options(10)函数和约束求值的数目.
options(11)函数导数求值的个数
options(12)约束求值的数目.
options(13)等式约束的数目.
options(14)函数求值的最大次数(默认值为100×变量个数).
options(15)用于目标-达到问题中的特殊目标.
options(16)优化过程中变量的最小梯度值.
options(17)优化过程中变量的最大梯度值.
options(18)步长设置(默认值为1或更小).