数学建模竞赛中常用软件的操作

数学建模竞赛中常用软件的操作本节主要介绍数学建模竞赛中常用软件MATLAB和Lingo的一些基本操作。一、Desktop简介在桌面双击MATLABb图标，或双击安装目录C:Program FilesMATLAB R2012abin下的MATLAB文件。启动后默认界面如下图。图1 Desktop操作桌面的外貌1. Command Window该窗口是进行MATLAB各种操作的主要窗口。在该窗内可以输入各类指令、函数、表达式；显示除了图形外所有的运算结果，错误时，给出相关出错提示。指令输入完后只有按回车键【Enter】才能执行；如果输入的指令不含赋值号，计算结果被赋于默认的变量ans。变量名和函数名对大小写敏感，变量第一个字符必须是英文字母，最多包含63个字符（英文、数字和下划线），不能包括空格、标点、运算符；不能使MATLAB的关键词和自用的变量名（eps,pi等）函数名（sin,exp等）、文件夹名（rwt,toolbox等）。在Matlab中有一些固定变量，例如(1) ans：在没有定义变量名时，系统默认变量名为ans；(2) eps：容许误差，非常小的数；(3) pi：即圆周率；(4) i, j：虚数单位；(5) inf：表示正无穷大，由1/0运算产生；(6) NaN（Not A Number）：表示不定值，由inf/inf或0/0运算产生；(7) nargin：函数的输入变量数目；(8) nargout：函数的输出变量数目。在MATLAB中，控制流关键字if, for, end等用蓝色字体表示；输入指令中的非控制指令、数字显示为黑色字体；字符串显示为紫色字体；注释为绿色字体；警告信息为红色字体。2 工作空间浏览器工作空间（Workspace）窗口用于浏览MATLAB中的变量。在工作空间窗口内，用户可以方便地查看、编辑存储的数据变量。表1 工作空间浏览器主要功能及其操作方法功 能操 作 方 法新变量创建单击该图标，在工作空间产生unnamed新变量；双击该变量，引出Variable Editor变量编辑器；可输入数据；可重新命名变量内容显示选中变量，点击该图标，则变量内容显示在Variable Editor变量编辑器中向内存装载文件数据选择MAT数据文件，单击该图标，引出Import Wizard界面，选择需要装载的数据把变量保存进文件选择一个或多个内存变量，单击该图标或单击右键，选择Save as，把这些变量保存到MAT数据文件图形绘制启动键绘制选定类型的图形绘图类型菜单引出键点击引出绘图类型菜单以供选择工作空间常用的管理指令有：（1）who及whos：查询指令（2）clear：清除工作空间中的所有变量 clear var1 var2：清除工作空间中的变量var1和var2 （3）save FileName：把全部内存变量保存为Filename.mat文件 save FileName v1 v2：把变量v1和v2保存为Filename.mat文件 save FileName v1 v2 append：把变量v1和v2添加到Filename.mat文件（4）load FileName：把Filename.mat文件中全部变量调入内存 load FileName v1 v2：把Filename.mat文件中v1 , v2调入内存3. Current Folder用户保存文件时，如果不专门指定目录名，则所存文件将保存在当前文件夹下。注意尽量不要把MATLAB所在的根目录或其任何子目录作为当前目录，以免破坏MATLAB原有文件的完整性。4. 向量与矩阵的产生MATLAB的操作对象为矩阵。标量被看作是11的矩阵，维行向量或列向量可看作是或维矩阵。行向量生成基本格式：start: step: end（step缺省时为1）例如：在命令窗口中提示符 后输入，x=1:5，按回车【Enter】后，显示输出结果如下：x = 1 2 3 4 5Matlab的每条命令后，若为逗号“，”或无标点符号，则显示命令的结果；若命令后为分号“；”，则禁止显示结果，用“%”来实现注释功能。利用下面两个函数也可以产生向量。linspace：线性等分向量，例如，a=linspace(0,2*pi,100)。logspace：对数等分向量，例如，a=logspace(0,2*pi,100)。100为插入点数（99等分），默认值为100。在矩阵的输入中，约定：(1) 元素之间用空格或逗号隔开；(2) 用中括号方 把所有元素括起来；(3) 用分号（；）说明行结束。在命令窗口，输入矩阵时，按Shift+Enter键表示开始输入新的一行，且要求所有的行具有相同的列。例如：a=1 2 3; 4,5 6; 7 8,9该矩阵一直保存在工作空间，直至被修改。二、M文件与编程1. M文件M文件有两种类型：脚本类M文件和函数类脚本文件。脚本类M文件是一串按照用户意图排列而成的Matlab指令集合。脚本文件运行后，产生的所有变量都驻留在Matlab的基本工作空间（base workspace）中。只有用clear清除或关闭指令窗清除，否则一直保留。基本工作空间随Matlab的启动而产生，只有当Matlab关闭时，基本空间才被删除。函数类M文件的第一行总是以function引导的函数申明行，还包括函数与外界交换数据的全部标称输入输出量（可有可无）。Matlab允许使用比标称数目少的输入输出量，实现对函数的调用。运行函数文件时，Matlab开辟一个临时工作空间，称之为函数工作空间（function workspace），所有中间变量都存放在工作空间中，该函数文件运行结束后，该函数临时工作空间及其中间变量立即删除。在函数文件中，对某脚本文件的调用而产生的所有变量也存放于函数工作空间，而不存放在基本工作空间。函数工作空间随M函数文件的被调用而产生，调用结束后删除，可产生任意多个，相对于基本工作空间而言是独立的、临时的。2. 编程MATLAB程序的编写主要有以下几种控制结构。（1）if-else-end条件控制表2 if-else-end分支结构的使用方法单分支双分支多分支 if expr (commands) Endif expr (commands1)else (commands2)endif expr1 (commands1)elseif expr2 (commands2) .else (commandsk)end当expr为逻辑1时, 指令组(commands)被执行。当expr为逻辑1时,指令组1执行;否则,指令组2执行。expr1,expr2,中,首先给出逻辑1的分支被执行;否则,分支k被执行。常被switch-case替代。expr为控制分支的条件表达式，通常为关系、逻辑表达式，其运算结果为“标量逻辑值1或0”。也可以为一般代数表达式，则任何非0值均等同于“逻辑1”。expr也可以进行数组间的关系、逻辑运算，其运算结果为逻辑数组，只有该数组不包含任何0元素时，expr控制的分支才被执行。expr可为空数组，Matlab认为条件为“假”，分支不被执行。（2）switch-case控制结构表3 switch-case控制结构的使用方法指 令 格 式含 义switch expr case value_1 (commands1) case value_2 (command2) . . . case value_k (commandsk) otherwise (commands)endexpr为根据此前给定的变量进行计算的表达式；value_1是给定的数值、字符串标量（或胞元数组）若expr结果与value_1（或其中的胞元元素）相等，则执行commands1该项是以上各项的“并”的“补”；即以上所有的case均不发生，则执行改组指令。【说明】switch后的表达式expr的值只能是标量数值或标量字符串。对于标量数值，比较：表达式=检测值i；对于字符串，利用strcmp（表达式，检测值i）比较。case后的检测值不仅可以是一个标量数值或字符串，还可以是一个胞元数组。此时将expr和胞元数组中每一个元素比较，如与某一胞元数组元素相等，认为比较结果为真，从而执行与该检测值相应的一组命令。（3） for循环和while循环表4 循环结构的使用for循环while循环for ix=array (commands) endwhile expression (commands)endix为循环变量，commands为循环体ix依次取array中元素，每取一个元素，执行循环体依次，直到ix大于array的最后一个元素跳出循环for循环的循环次数确定检测expression的值，若为逻辑真（非0）则执行循环体，否则结束循环。while循环的循环次数不确定三、读写Excel文件现在竞赛的附件大多是以Excel电子表格的形式给出，所以读写Excel电子表格文件的方法就显得十分重要。MATLAB提供了读取excel文件（后缀为.xls）的函数xlsread，其调用格式有：num = xlsread(filename)：从工作文件夹中，从A1单元格开始读取filename文件，把数据返回给num。num = xlsread(filename, sheet)：选择所在的sheet，比如说sheet1。num = xlsread(filename ,range)：所在的单元格范围，比如说range=A1:A8。num = xlsread(filename, sheet, range)：sheet和单元格范围同时限制。num,txt= xlsread(filename, .) ：把返回的数据与文本分开保存。num,txt,raw = xlsread(filename, .)：分开保存的同时，有把num和txt保存到raw里，形成一个单一变量。输出结果num中存储的是数据，txt存储的是文本，raw为所有的数据和文件。例如2011年全国大学生数学建模竞赛的A题城市表层土壤重金属污染分析中的附件，可以按如下的方法读取附件数据。 读取cumcm2011A附件_数据.xls中附件1中有关取样点位置及其所属功能区数据：A=xlsread(cumcm2011A附件_数据.xls,附件1,B4:E322)读取附件2中8种主要重金属元素的浓度数据：format short gB=xlsread(cumcm2011A附件_数据.xls,附件2,B4:I322)读取附件3中8种主要重金属元素的背景值（平均值，标准偏差和范围）：C=xlsread(cumcm2011A附件_数据.xls,附件3,B4:C11)MATLAB提供了读取excel写入的函数xlswrite，其调用格式有：status,message = xlswrite(filename, array, sheet, range)其中array为要保存的数据矩阵，sheet为工作表，range制定存储的区域，例如C1:C2。status为返回完成状态值，如果写入成功，则status为1；反之写入失败，则status为0。message包括了错误、警告信息和标识信息。例如：将数据写入testdata.xls文件的Sheet1。 xlswrite(testdata.xls, 12.7, 5.02, -98, 63.9, 0, -.2, 56)将文本和数值混合数据从E1单元格开始写入testdata.xls文件的Sheet2： d = Time,Temperature; 12,98; 13,99; 14,97;xlswrite(testdata.xls, d, 2, E1)在testdata.xls中新建表格Temperature： status,msginfo = xlswrite(testdata.xls, d,Temperatures, E1)在程序中，经常要对数据进行排序，常有以下函数。（1）sort函数sort(A)：若A是列（或行）向量，默认对A进行升序排列，而sort(A,descend)是降序排序。若A是矩阵，默认对A的各列进行升序排列。sort(A,dim)：dim=1时等效sort(A)，dim=2时表示对A中的各行元素升序排列。（2）sortrow函数B=sortrows(A,n)：A为要排序的矩阵，n表示按第几列进行排序，正整数代表升序，负整数代表降序。B=sortrows(a,2)表示将矩阵A按照第2列的升序排列，其他列不作升序排列，而是将跟随着第2列的变化，保持一一对应。B=sortrows(a,-3)表示将矩阵A按照第3列的降序排列，其他列将跟随着第3列的变化，保持一一对应。四、插值与拟合、回归分析在生产和实验中，由于函数的表达式不便于计算或这没有表达式而只有在给定点的函数值（或其导数值），为此，我们希望建立一个简单而且便于计算的近似函数，来逼近函数，这就用到差值和拟合方法。1 插值插值就是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。MATLAB中的插值函数主要有以下几个。（1） interp1函数Matlab中用于一维数据插值的函数是interp1，其调用格式为：yi = interp1(x, y, xi, method , extrap )该命令用于找出由参量x决定的一元函数y=y(x)在点xi处的值yi。其中x,y为插值节点的横坐标和纵坐标，yi为在被插值点xi处的插值结果；x,y为向量， method 表示采用的插值方法，Matlab提供的插值方法有几种：nearest：最近邻点插值算法；linear：线性插值spline：三次样条函数插值；pchip：分段三次Hermite插值；cubic：与pchip操作相同；缺省时表示线性插值。 extrap表示对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。注意：所有的插值方法都要求x是单调的。例如，在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为12，9，9，1，0，18 ，24，28，27，25，20，18，15，13推测中午13时的温度。我们采用三次样条插值，在命令窗口输入： x=0:2:24;a=13;y=12 9 9 10 18 24 28 27 25 20 18 15 13;y1=interp1(x,y,a,spline)输出结果为：y1 = 27.8725 （2） interp2函数Matlab中用于二维数据插值的函数是interp2，其调用格式为：zi = interp2(x, y, z, xi, yi, method)该命令用于找出由参量x,y决定的二元函数z=z(x,y)在点（xi,yi）处的值zi。其中返回矩阵为zi，其元素为对应于参量xi与yi（可以是向量、或同型矩阵）的元素，若xi与yi 中有在x与y范围之外的点，则相应地返回NaN（Not a Number），method和interp1一样，常用的是cubic （双三次插值），缺省为linear（双线性插值算法）。例如，在命令窗口输入X,Y = meshgrid(-3:.25:3);Z = peaks(X,Y); %具有两个变量的采样函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点。XI,YI = meshgrid(-3:.125:3);ZZ = interp2(X,Y,Z,XI,YI);surf(X,Y,Z);hold on;surf(XI,YI,ZZ+15) %为作比较，将插值曲面向上平移15单位。hold off（3） griddata函数 griddata也是一种常用的二维插值方法，其调用格式为：zi =griddata (x, y, z, xi, yi, method)该命令用于找出由参量x,y决定的二元函数z=z(x,y)在点（xi,yi）处的值zi。它和interp2的区别在于，interp2的插值数据必须是矩形域，即已知数据点（x,y）组成规则的矩阵，可使用meshgid生成。而griddata函数的已知数据点（x,y）不要求规则排列，特别是对试验中随机没有规律采取的数据进行插值具有很好的效果。method包括：nearest（最近邻点插值）；linear（线性插值，为缺省算法）； cubic（ 基于三角形的三次插值）和v4（Matlab4中的griddata 算法）。例如，有一组散乱数据点矩阵如下：A=1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1;4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1;1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2;4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2;2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3;2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4;3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5;3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6;3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7;x=A(:,1);y=A(:,2);z=A(:,3);X,Y,Z=griddata(x,y,z,linspace(min(x),max(x),20),linspace(min(y),max(y),20),v4);mesh(X,Y,Z); hold onplot3(x,y,z,o); hold off;2011年的A题城市表层土壤重金属污染分析，可用插值拟合的方法获得各重金属污染物浓度的空间分布。由于空间数据是不规则的，较好的方法是用散乱数据插值，例如Kriging插值、Shepard插值等。也可以用其他方法插值拟合，但应明确所使用的方法，并作出分析，不能只简单套用软件。4spline函数该函数是利用三次样条对数据进行插值，其调用格式为：yy = spline(x,y,xx)该命令用三次样条插值计算出由向量x与y 确定的一元函数y=f(x)在点xx 处的值。若参量y 是一矩阵，则以y 的每一列和x 配对，再分别计算由它们确定的函数在点xx 处的值。则yy 是一个阶数为length(xx)*size(y,2)的矩阵。例如，对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值。在命令窗口输入：x = 0 2 4 5 8 12 12.8 17.2 19.9 20;y = exp(x).*sin(x);xx = 0:.25:20;yy = spline(x,y,xx);plot(x,y,o,xx,yy) 2拟合拟合就是用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。如果已知某函数的若干离散函数值f1,f2,fn，通过调整该函数中若干待定系数f(1, 2,n)，使得该函数与已知点集的差别（最小二乘意义）最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫做非线性拟合或者非线性回归。MATLAB中提供了线性最小二乘拟合和非线性最小二乘拟合函数。（1） polyfit函数ployfit函数是多项式拟合函数，其调用格式为：p = polyfit(x,y,n)其中x,y为长度相同的向量，n为拟合多项式的次数，返回值p是拟合多项式的系数向量，幂次由高到低，拟合多项式在x处的值可以通过y=polyval(p,x)来计算。例如，2004年全国大学生数学建模竞赛C题，饮酒驾车问题，时间和血液中酒精浓度的函数关系，可以利用polyfit进行拟合。可输入： time=0.25, 0.5 0.75, 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16;vol=20 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4;plot(time,vol,o)p=polyfit(time,vol,7) % 7次多项式拟合f=polyval(p,time); %求多项式的值hold on plot(time,f)hold off输出7多项次的系数向量为：p = 0.0002 -0.0113 0.2817 -3.6307 25.5345 -94.6710 153.8757 -1.2592如果要求血液中酒精浓度的最大值，可用以下指令：roots(polyder(p)重庆文理学院的一个参赛队于2004年关于“饮酒驾车”问题获得一等奖后，研制了“酒后安全驾车时刻表”和“人体内酒精浓度反推软件”两个产品，并在重庆永川交警大队得到使用。Matlab提供了两个非线性最小二乘拟合函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的。（2）lsqcurvefit函数该函数用来进行非线性拟合，其调用格式为：x = lsqcurvefit (fun,x0,xdata,ydata,options);其中，fun为事先建立的拟合函数F(x,xdata)，其中自变量x表示拟合函数中的待定参数，xdata为已知拟合节点的x坐标，x0为待定参数x的迭代初始值，xdata,ydata为已知数据点的x和y坐标，options是一些控制参数。lsqcurvefit函数用来求含参数x（向量）的向量值函数F(x,xdata)=f(x,data1), f(x,data2), f(x,datan)中的参数x（向量），使得最小。例如，根据表5.22中的数据，利用lsqcurvefit函数拟合。表5已知数据点x1002003004005006007008009001000y10-34.544.995.355.655.906.106.266.396.506.59首先建立拟合函数的M文件myfit1.m，其内容如下：function f = myfit1(x,xdata) f = x(1)+x(2)*exp(-0.02*x(3)*xdata); 其中x(1),x(2),x(3)分别表示拟合曲线中的参数a,b,k。然后在命令窗口输入： xdata = 100:100:1000;ydata = 1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;x0=0.2,0.05,0.05;x = lsqcurvefit (myfit1,x0,xdata,ydata)f= myfit1(x,xdata) （3） lsqnonlin函数该函数用来进行非线性拟合，其调用格式为：x = lsqnonlin (fun,x0, options);其中fun为事先建立的拟合函数f(x)，其中自变量x表示拟合函数中的待定参数，x0为待定参数x的迭代初始值，options是一些控制参数。由于lsqnonlin中定义的拟合函数的自变量是x，所以已知参数xdata,ydata应写在该函数中。lsqnonlin函数用来求含参量x（向量）的向量值函数f(x)=f1(x), f2(x),，fn(x)中的参量x，使得最小。其中fi(x)=f(x,xdatai,ydaytai)=F(x, xdatai)- ydatai。例如，根据表5中的数据，利用lsqnonlin函数拟合。首先建立拟合函数的M文件myfit2.m，其内容如下：function f=myfit2(x)xdata=100:100:1000;ydata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=x(1)+x(2)*exp(-0.02*x(3)*xdata)- ydata;然后在命令窗口输入： x0=0.2,0.05,0.05;x = lsqnonlin(myfit2,x0)f = myfit2(x)该结果与lsqcurvefit函数拟合的结果相同。3. 回归分析MATLAB中提供了一些线性和非线性回归分析函数。（1） regress函数一般地，称为高斯马尔柯夫线性模型（k元线性回归模型）。，Matlab提供了多元线性回归函数regress，采用的是最小二乘估计，其调用格式有：b = regress (y,x)返回值为线性模型y = x*b的回归系数向量。其中x为n(k+1)矩阵，行对应于观测值，列对应于预测变量，y为n1 向量，为因变量，一元线性回归可取k=1。b,bint,r,rint,stats = regress(y,x,alpha)其中bint是回归系数的区间估计，r是残差，rint是置信区间，stats是用于检验回归模型的统计量，有四个数值：相关系数r2，F值，与F对应的概率P和误差方差估计，alpha是显著性水平（缺省的时候为0.05）。相关系数r2越接近1，说明回归方程越显著；与F对应的概率P y=8.8818 8.9487 9.0541 9.1545 9.2693 9.4289 9.6160 9.8150 9.9825 10.1558 10.3193;k= 7.8381 7.9167 8.0048 8.1026 8.2556 8.5822 8.8287 9.0756 9.2175 9.4148 9.6198;l= 8.3871 8.3872 8.3935 8.3971 8.4025 8.4048 8.4079 8.4141 8.4261 8.4377 8.4444;t= 9.9551 9.9057 10.0972 9.9537 9.9370 9.9449 9.9636 10.1291 10.1573 10.2944 10.2093;x= ones(size(k) k l t; b,bint,r,rint,stats=regress(y,x)输出结果为b = -55.4988 0.5644 7.1254 0.0222stats = 1.0e+003 * 0.0010 2.1022 0.0000 0.0000该结果说明y=-55.4988+0.5644k+7.1254l+0.0222t，stats中的数据说明r2=1，F=2102.2，p=0，由于p x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;y=100 75 80 70 50 65 90 100 110 60;x=x1 x2;rstool(x,y,purequadratic)其输出如图2。图2多元二项式回归在图2中x1上面的方框中输入1000，x2上面的方框中输入6，在图形框左侧的“Predicted Y1”下方的数据变为88.47981，即预测出平均收入为1000，价格为6时的商品需求量为88.4791。单击图形框左边Export，则出现图3对话框，可以将回归参数beta、剩余标准差rmse和残差residuals传送到Matlab的工作区中。图3输出对话框在命令窗口输入： betarmse可知故回归模型为：，剩余标准差为4.5362, 说明此回归模型的显著性较好。还可将该模型转化为多元线性回归模型，利用regress函数进求解。可输入： x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;X=ones(10,1) x1 x2 (x1.2) (x2.2);b,bint,r,rint,stats=regress(y,X);b,stats输出的回归参数和rstool函数的结果相同。（3） nlinfit函数nlinfit函数用来确定非线性回归系数，调用格式为：beta, r, J = nlinfit (x, y, modelfun, beta0)其中输入数据x,y分别为np维矩阵和n维列向量，对于一元非线性回归，取p=1即可；modelfun为事先定义的非线性回归函数的M文件，是回归系数beta和x的函数；beta0是回归系数的初值，输出参数beta是估计出的回归系数，r为残差，J为Jacobain矩阵。例如，根据表5中的数据，利用nlinfit函数进行非线性回归，回归函数为。首先建立回归函数的M文件myfit3.m，内容如下：function f = myfit3(beta,xdata) f = beta(1)+beta(2)*exp(-0.02*beta(3)*xdata); 在命令窗口输入： xdata = 100:100:1000;ydata = 1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;beta0 = 0.1 0.1 0.2;beta,r,J = nlinfit(xdata,ydata,myfit3,beta0) ;beta输出结果为：beta = 0.0070 -0.0030 0.1012这与例9结果（ 0.0069 -0.0029 0.0809）很接近。将其和前例的拟合曲线对比后可见，nlinfit函数的回归效果明显优于lsqcurvefit函数，对比结果见图4（红色曲线为nlinfit函数得到的回归曲线，蓝颜色的为lsqcurvefit函数得到的拟合曲线）。 100200300400500600700800900100044.555.566.57x 10-3图4 对比图五、胞元数组与构架数组1. 胞元数组胞元数组的基本组分是胞元（Cell）。每一胞元在数组中平等，以下标区分。同一胞元数组中不同胞元可以存放不同类型和不同大小的数据，如任意维数值数组，字符串数组、符号对象等。胞元标识（Cell Indexing）：例A(2,3)表示A胞元数组第2行第3列的胞元；胞元内容编址（Content Adressing）：例A2,3表示A胞元数组第2行第3列的胞元存放的内容。例如clearC_str=这是胞元数组创建例子;R=reshape(1:9,3,3);%按列排序Cn=1+2i;S_sym=sym(sin(-3*t)*exp(-t);%创建胞元数组方法之一B1,1=C_str;B1,2=R;B2,1=Cn;B2,2=S_sym;%胞元的援引a=B(1,2)%class(a)a = 3x3 doubleans =cell%胞元内容的援引b=B1,2% class(b)b = 1 4 7 2 5 8 3 6 9ans =double2. 构架数组组织数据能力比胞元数组更强，富于变化。架构数组的基本组分是构架（structure）。数组中每个架构平等，以下标区分。构架必须在划分“域”后才能使用。数据不能直接存放于构架，而只能存放在域中。构架的域可以存放任何类型、任何大小的数组（如任意维数值数组，字符串数组、符号对象等），不同构架的同名域中存放的内容可以不同。构架和域名间有小黑点“.”。构架数组演示（1）直接对域赋值产生“单构架”clearG.name=一号房;G.volume=2000;G.temperature=31.2,30.4,31.6,28.7; G.humidity=62.1,59.5,57.7,61.5;63,60,58.1,62.3;（2）向第二个构架赋值，形成构架数组G(2).name=二号房;G(2).volume=2400;（3）构架数组显示GG = 1x2 struct array with fields: name volume temperature humidity（4）构架元素显示G(1)ans = name: 一号房 volume: 2000 temperature: 31.2000 30.4000 31.6000 28.7000 humidity: 2x4 double（5）构架域的显示G.humidityans = 62.1000 59.5000 57.7000 61.5000 63.0000 60.0000 58.1000 62.3000ans = 对于单构架，输入构架名可输出该构架的域中的内容，非单构架通常只能得到该构架的结构信息；指令窗输入构架域名，显示每一构架的该域内容。一般代数方程组的符号求解函数S=solve(eq1,eq2,eqn,v1,v2,vn)的结果S就是构架数组形式，例如：S = solve(x+y=5,x-y=1,x,y)S = x: 1x1 sym y: 1x1 sym S.xans =3 S.yans =2六、图像处理在全国大学生数学建模竞赛中，与图像处理有关的赛题有三道，分别是2001年A题三维血管重建，2008年A题数码相机定位，2013年B题碎纸片的拼接复原。下面简单介绍图像处理中的几个简单函数。MATLAB支持的图像格式有*.bmp,*.jpg,*.ico,*.pcx,*.png,*.tif等。1. imread函数imread函数用于读入各种图像文件，其调用格式为：A = imread(filename,fmt)其中filename为指定工作目录下的图像文件名，fmt为图像格式，可选值为bmp,jpg,ico,pcx,png,tif等，如果不指定，Matlab自动根据文件头确定。A为读取的图像矩阵。2. imwrite函数imwrite函数用于输出图像，其调用格式为：imwrite(A, filename,fmt)将矩阵A写入当前工作目录下的filename.fmt文件。3. imshow函数该函数是最常用的显示各种图像的函数，其调用格式有：imshow(I,n)：用于显示灰度图像，n为灰度级，默认为256；imshow(I,low,high)：用于显示灰度图像，low,high为图像的数域，如不清楚，可使用空向量。例如，输入 I = imread(cameraman.tif); imshow(I,)4. imnoise函数该函数可以对图像添加一些典型的噪声。其调用格式为：imnoise(I, type,parameters)其中type为噪声类型，包括：gaussian：高斯白噪声，localvar ：0均值高斯白噪声poisson：泊松噪声salt & pepper：椒盐噪声speckle：乘性噪声在命令窗口输入： I1 = imread(cameraman.tif); subplot(2,2,1);imshow(I1);title(原始图像);I2 = imnoise(I1, gaussian); subplot(2,2,2);imshow(I2);title(高斯白噪声);I3 = imnoise(I1,salt & pepper);subplot(2,2,3); imshow(I3);title(椒盐噪声);I4 = imnoise(I1,speckle);subplot(2,2,4);imshow(I4);title(乘性噪声);5. fspecial函数fspecial函数用于建立预定义的滤波算子，其语法格式为：h = fspecial(type)h = fspecial(type,para)其中type指定算子的类型，para指定相应的参数；type的类型有：(1) average：为均值滤波，参数为hsize代表模板尺寸，默认值为3 3。(2) disk：为圆形区域均值滤波，参数为radius代表区域半径，默认值为5。(3) gaussian：为高斯低通滤波，有两个参数，hsize表示模板尺寸，默认值为3 3，sigma为滤波器的标准值，单位为像素，默认值为0.5。(4) laplacian ：为拉普拉斯算子，参数alpha用于控制算子形状，取值范围为0,1，默认值为0.2。(5) log ：为拉普拉斯高斯算子，有两个参数，hsize表示模板尺寸，默认值为3 3，sigma为滤波器的标准差，单位为像素，默认值为0.5。(6) motion：为运动模糊算子，有两个参数，表示摄像物体逆时针方向以theta角度运动了len个像素，len的默认值为9，theta的默认值为0。(7) prewitt：用于边缘增强，大小为3 3，无参数。(8) sobel：用于边缘提取，无参数。(9) unsharp：为对比度增强滤波器。参数alpha用于控制滤波器的形状，范围为0,1，默认值为0.2。6. filter2函数该函数是对图像进行二维线性滤波，经常和fspecial函数配合使用，其调用格式为：A = filter2(type,I,shape)利用指定滤波器type对图像I进行滤波，其大小由参数 shape 确定，其取值如下：full：返回二维相关的全部结果，size(A)size(I)；same：返回二维互相关结果的中间部分，A与I大小相同（默认）； valid：返回在二维互相关过程中，未使用边缘补 0进行计算的结果， size(A) I = imread(cameraman.tif); imshow(I)In = imnoise(I, gauss);figure; imshow(In);f1 = fspecial(average); I1 = filter2(f1,In); figure; imshow(I1,) f2 = fspecial(disk); I2 = filter2(f2,In); figure; imshow(I2,)f3 = fspecial(gaussian); I3 = filter2(f3,In); figure; imshow(I3,)f4 = fspecial(motion); I4 = filter2(f4,In); figure; imshow(I4,)7. medfilt2函数该函数是对图像进行中值滤波，其调用格式为：A = medfilt2(I,m,n)用指定大小为mn的窗口对图像I进行中值滤波，默认为33。8. wiener2函数该函数实现图像自适应滤波。wiener2函数根据图像的局部方差来调整滤波器的输出。其调用格式为：A = wiener2(I,m,n)用指定大小为mn的窗口对图像I进行自适应滤波，默认为33。例如对上例中添加高斯白噪声的图像进行中值滤波和自适应滤波，在命令窗口输入 A1 = medfilt2(In); figure; imshow(A1,)A2 = wiener2 (In); figure; imshow(A2,)9. edge函数边缘检测技术对于图像处理非常重要，Matlab提供的edge函数可以实现图像的边缘检测，其调用格式为：IE = edge(I, type)其中type为边缘检测算子，包括：sobel,prewitt,roberts,log,zerocross,canny。例如，对于cameraman图像，不加噪声时进行边缘检测，在命令窗口输入： I = imread(cameraman.tif);BW_sobel = edge(I,sobel); figure; imshow(BW_sobel,)BW_prewitt = edge(I,prewitt); figure; imshow(BW_prewitt,)BW_roberts = edge(I,roberts); figure; imshow(BW_roberts,)BW_laplace = edge(I,log); figure; imshow(BW_laplace,)BW_zerocross=edge(I, zerocross ); figure; imshow(BW_zerocross,)BW_canny = edge(I,canny); figure; imshow(BW_canny,)在2013高教社杯全国大学生数学建模竞赛B题碎纸片的拼接复原中，每一附件中都含有大量图像，由于文件名一般具有一定的规律性，故可采用胞元数组的方式一次读入。例如，可用以下语句：Files = dir(*.bmp); %获取当前文件夹下后缀为bmp的文件for i=1:length(Files) %得到图像数量 Ai =im2bw(imread(Files(i).name); %将图像存入胞元数组A中 end胞元数组A就存储了该路径下所有的后缀为bmp的图像文件。七、优化问题MATLAB软件可以线性规划、0-1规划、非线性规划、二次规划、多目标规划等优化问题。（1） 线性规划（Linear Programming，简写为LP）问题线性规划问题就是求多变量线性函数在线性约束条件下的最优值。满足约束条件的解称为可行解，所有可行解构成的集合称为可行域，满足目标式的可行解称为最优解。MATLAB解决的线性规划问题的标准形式为：min 其中为列向量，为矩阵。其它形式的线性规划问题都可经过适当变换化为此标准形