数学建模培训Matlab.ppt

Matlab软件的使用,主要内容,一、Matlab初步Matlab环境Matlab语言Matlab矩阵Matlab程序设计二、Matlab绘图二维作图三维作图其它作图三、数据分析数据的导入导出多项式的处理插值与拟合四、方程求解与优化问题求解方程求解无约束问题求解（极值与最值）约束问题求解,一、Matlab初步,Matlab的历史Matlab是Matrix和Laboratory的组合八十年代初：作为免费软件，C语言编写现已发展成为适合众多学科、多种工作平台的功能强劲的大型软件。Matlab的构成Matlab的语言、Matlab的工作环境、Matlab的工具箱、Matlab的APIMatlab的特点数值运算全面、表示方法简单、丰富的工具箱Matlab的目录结构开放的编程系统大部功能是通过后缀为.m的文本文件实现的,1、MATLAB环境,工作界面的显示属性调整,“File”|”Preferences”选项,命令窗口(Command Window),输入及运行分号的作用命令窗口常用快捷键,工作界面,历史命令窗口（Command History）记录用户在Matlab命令窗口中输入的所有的命令包括每次启动Matlab的时间和运行的所有命令行单击右键，对历史命令进行编辑（剪切/复制/运行/创建m文件/快捷方式/profile code等）工作空间浏览器（Workspace Browser）运行Matlab的程序或命令所生成的所有变量和Matlab提供的常量构成的空间查询和编辑已定义变量通过右键菜单进行编辑或绘图等相关操作who(whos)、clear、save命令查看、清除、保存工作空间的所有变量程序编辑窗口（Editor）一个内置的具有编辑和调试功能的程序编辑器,工作界面,当前目录窗口（Current Directory）MATLAB搜索路径机制和搜索顺序修改Matlab的搜索路径查看任一路径下的所有文件添加自己的搜索路径MATLAB当前目录Matlab的帮助系统输入帮助命令（help、lookfor、demo）通过“Help”菜单,工作界面,2、Matlab语言,变量:字母开头的数字、字母、下划线序列长度不超过63个字符，区分大小写特殊变量与常量列表,语言初步,数据类型,MATLAB的数据类型数值型（整型、单精度型、双精度型、复数类型）逻辑型字符型（单引号括起来的一个或多个字符）函数句柄MATLAB的数据结构矩阵、数组 表现形式相同，用方括号作为界定符 运算不同，数组运算针对每个元素结构体 每个域可保存不同维数、不同类型的数据 使用点号“.”运算元胞（单元）数组 各单元可保存不同维数、不同类型的数据 使用大括号,语言初步,基本函数,语言初步,运算符,(1)算术运算基本运算:+-*/(右除)(左除)点运算(数组运算):+-.*./.(2)关系运算运算符：、=、=、=两个标量比较，直接比较，结果为1或0;两个同维矩阵比较，按对应元素比较，结果是0、1矩阵；标量和矩阵比较，标量与矩阵每一元素比较，结果是0、1矩阵。(3)逻辑运算运算符：&、|、捷径运算符：&、|,数据输出格式,format 命令指定 数值型数据的输出格式,3、矩阵_生成,Matlab运算的基本单元是矩阵(1)行向量的生成直接生成 矩阵元素列入方括号中,元素之间用逗号或空格分隔利用冒号生成行向量 格式1：初值:终值 产生从初值开始到终值结束增量为1的行向量 格式2：初值:步长:终值 产生从初值开始到终值结束增量为步长的行向量利用linspace函数生成向量 格式：linspace(a,b,n)生成n个元素的行向量，元素在a,b之间平均分布，n的默认值为100.,向量和矩阵,矩阵_生成,(2)矩阵的生成直接生成:矩阵元素列入方括号中,元素之间用逗号或空格分隔,行与行之间用分号分开通过函数zeros、ones、rand、randn分别产生元素全为零、全为1、随机数、正态分布随机数的矩阵通过函数Compan、hilb、magic、pascal分别产生伴随阵、Hilbert阵、魔方阵、Pascal阵通过函数eye产生单位矩阵,矩阵_处理,(1)矩阵元素的引用格式1：A(u,v)其中u,v的取值：向量、标量A(i,j)、A(i,:)、A(:,j)、A(m1:m2,n1:n2)格式2：A(n)矩阵一个按列优先排列的列向量,A(n)表示序号为n的元素格式3：A(逻辑数组)A(A10)(2)矩阵元素的修改A(i,j)=新值A(m1:m2,n1:n2)=B注意B的行数与列数与A中相应的块匹配(3)矩阵维数的缩维（删除行与列）A(:,n)=A(n,:)=(4)矩阵维数的扩维A(i,j)=值i(j)的值大于行数(列数),向量和矩阵,矩阵_操作,矩阵的大小size(A)、length(A)分别返回行列数、行列中大的数改变矩阵的形状reshape(A,m,n)将矩阵调整为m行n列构造对角阵和三角矩阵diag(A),diag(v)提取对角线无素组成列向量、以向量作对角线元素产生对角矩阵tril(A),triu(A)：提取下三角矩阵、提取上三角矩阵矩阵的转置与旋转 rot90(A)：逆时针旋转90度flipud(A),fliplr(A)：分别将矩阵上下翻转、左右翻转行列式、逆、秩、迹det、inv、rank、trace特征值和特征向量Eeig(A)求A的特征值V,D=eig(A)求A的特征值和特征向量,向量和矩阵,数集_操作,数集在MATLAB中表现为元素互斥的向量数集的交、并、余intersect、union、setdiff判断元素是否属于数集ismember去掉数集中重复的元素unique,向量和矩阵,4、Matlab程序设计,用Matlab语言编写的可在Matlab中运行的程序,称为M文件。扩展名为.m，可用任何文本编辑器编辑种类：命令文件(脚本文件)、函数文件,基本语句input(string)pausepause(n)breakreturndisp,条件语句if 语句swhich语句 switch e case c1 语句1 case c2 语句2 otherwise 语句 endc1,c2的值若有多个用“”隔开或用大括号括起来。,循环语句forend循环for 变量=矩阵 循环体end矩阵为向量,将元素依次赋给变量,否则将矩阵每一列赋给变量whileend循环while 表达式e 语句块end,命令文件,命令文件是包含一系列MATLAB语句的简单文件它不接受输入参数，输出结果返回在命令窗口变量保存在工作空间中运行方式：命令窗口中输入命令文件的文件名编辑窗口中选中部分内容，通过菜单可运行例：随机产生20个两位整数，将小于平均值的偶数输出。%ls.mt=floor(10+10*rand(1,20);m=mean(t);t(tm&rem(t,2),例：%ffibno.mf=1,1;i=1;while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1);i=i+1;endf,i例：当n=1000时，求1/4+1/16+1/(4n)的值。方法s=0;for i=1:1000 s=s+1/4i;enddisp(s);方法 t=1:1000;sum(1/4).t)方法2属于向量法，较方法1优,函数文件,函数文件是按照一定格式编写的M文件，也称为永久性函数。使用自已的局部变量，接受输入参数，也能返回输出参数定义function y1,y2,=函数名(n1,n2,)注释说明部分：%开头函数体：包括进行运算和赋值操作的所有Matlab代码保存：一般以函数名为文件名保存为.m文件调用：输出参数表=函数名(输入参数表）stat(1 2 3 4)m,s=stat(1 2 3 4),function f=ffibno(n)%FFIBNO%f=ffibno(n)%2003/5/20f=1,1;i=1;while f(i)+f(i+1)n f(i+2)=f(i)+f(i+1);i=i+1;end,function mean,stdev=stat(x)n=length(x);mean=sum(x)/n;stdev=sqrt(sum(x-mean).2)/n);,临时函数临时性函数又包括内联函数和匿名函数。内联函数由inline函数建立，其格式为：f=inline(expr,arg1,arg2,.)例：f1=inline(sin(x)*cos(y),x,y);y=f1(pi/8,pi/9)匿名函数由符建立，其格式为：f=(arg1,arg2,.)expr例：f2=(x,y)sin(x)*cos(y);y=f2(pi/8,pi/9),2006年B题,%z=load(z1.txt);s=zeros(1,5);n=zeros(1,5);for i=1:length(z)if(z(i,2)=0|z(i,2)=1)s(1)=s(1)+z(i,3);n(1)=n(1)+1;elseif(z(i,2)=2 disp(t),二、Matlab绘图,1、二维图形plot(y)、plot(x,y)、plot(x,y,选项)x与y为向量 x作为横坐标,y作为纵坐标,绘制连线图y为有一维与向量x同维的矩阵 绘制多条色彩不同的连线图 x=linspace(0,2*pi,100);y=sin(x);1+sin(x);2+sin(x);plot(x,y)x,y均为矩阵 各取列向量,绘制多条色彩不同的连线图 x1=linspace(0,2*pi,100);x2=linspace(0,3*pi,100);x3=linspace(0,4*pi,100);x=x1;x2;x3;y=sin(x1);1+sin(x2);2+sin(x3);plot(x,y),二维图形,plot(x1,y1,x2,y2,)plot(x1,y1,选项1,x2,y2,选项2,)每对x、y按选项要求绘制曲线选项为单引号作为定界符的字符序列,用来控制线型、颜色、数据点和标记符等x=linspace(0,2*pi,100);y1=sin(x);y2=1+sin(x);y3=2+sin(x);plot(x,y1,r,x,y2,g,x,y3,b),二维图形,设置曲线的样式,二维图形,例:x=linspace(0,2*pi,1000);y1=0.2*exp(-0.5*x).*cos(4*pi*x);y2=2*exp(-0.5*x).*cos(pi*x);k=find(abs(y1-y2)1e-2);x1=x(k);y3=0.2*exp(-0.5*x1).*cos(4*pi*x1);plot(x,y1,x,y2,k:,x1,y3,bp);,图形保持与图形窗口,图形保持(同一绘图窗口陆续添加图形)hold on|off例 将上述三个图形依次添加到同一窗口中plot(x,y1);hold on;plot(x,y2,k:);plot(x1,y3,bp);图形窗口的创建(不同的图形在不同窗口显示)figure 新建绘图窗口figure(n)将第n个窗口作为当前窗口例 将上述三个图形在不同的窗口中显示plot(x,y1);figure;plot(x,y2,k:);figure;plot(x1,y3,bp);,二维图形,子图形的创建和控制,subpolt(m,n,p)把图形窗口分成m*n个小区域，并指定第p个为当前的绘制区域例t=0:pi/20:2*pi;subplot(2,2,1);plot(cos(t),sin(t);subplot(2,2,2);plot(t,sin(t);subplot(2,2,3);z=sin(t).*cos(t);plot(t,z);subplot(2,2,4);z=sin(t).3+cos(t).3;plot(t,z);subplot(111)将绘图窗口还成原成一个,二维图形,图形标注和坐标控制,图标标注title(字符串)xlabel(字符串)、ylabel(字符串)text(x,y,字符串)纵横比的调整axis squareaxis equalaxis normalt=0:0.1:2*pi;plot(sin(t),2*cos(t)坐标系的调整axis on|off格式：axis(xmin,xmax,ymin,ymax)例：x=0:1/3000:1;y=cos(tan(pi*x);figure(1);subplot(2,1,1);plot(x,y)title(复杂函数)subplot(2,1,2);plot(x,y)axis(0.4 0.6-1 1);title(复杂函数的局部透视),二维图形,Text,二维绘图的图例标注说明,二维图形,图形的可视化编辑,Matlab在图形窗口中提供了可视化的图形编辑工具利用这些编辑工具可完成各种图形对象的编辑处理Insert菜单的使用 添加坐标轴说明、图形标题、图例、箭头、线条、文本等属性编辑器的使用 进入图形编辑状态（“Tools-Edit Plot”）双击图形对象,打开属性编辑器,在相应的栏目中设置,二维图形,2、绘制三维图形,三维曲线图（plot）格式：plot3(x,y,z)plot3(x,y,z,选项)x,y,z为同维向量或矩阵：向量：它们相应的元素构成三维曲线的数据点坐标；矩阵：它们相应的列构成三维曲线数据点的坐标。选项同二维，控制线型、色彩、数据点标号类型。例：t=0:pi/50:10*pi;plot3(sin(t),cos(t),t)生成区域的坐标meshgrid(m:l:n)将区域m,n*m,n以步长l 划分为矩形网格meshgrid(1:3,4:7)例x,y=meshgrid(1:5);z=x+y;plot3(x,y,z)x,y=meshgrid(-2:0.1:2);z=x.*exp(-x.4-y.4);plot3(x,y,z),三维图形,三维网线图(mesh),原理：在x-y平面内一个长方形区域,采用与坐标轴平行的直线将其分隔;计算矩形网格上的函数值(即z的坐标值),得到三维空间中的数据点;将点分别用平行于xz平面和平行于yz平面内的曲线连接,形成网线图.格式：mesh(z)以z矩阵元素以及下标为数据点,绘制网线图mesh(x,y,z)x,y是向量，则x的长度=矩阵z的列数，y的长度=矩阵z的行数。x,y,z是同维矩阵，数据点分别取自三个矩阵meshc(x,y,z):生成具有基本等高线的网格图例：x,y=meshgrid(-8:0.5:8);z=sqrt(x.2+y.2);z1=sin(z)./z;mesh(x,y,z)和 mesh(x,y,z1)z0=zeros(size(z);plot3(x,y,z0,r*),三维图形,着色表面图(surf),原理：表面图是指把网线图表面的网格围成的小区域（或者叫补片）用不同的颜色填充形成的彩色表面。函数surf的格式同meshsurfc、surfl分别绘制带等高线的着色表面图和可以控制光照效应的着色表面图附加命令shading flat:去掉连接线条,平滑当前图形的颜色shading interp:去掉连接线条,在各片之间使用颜色插值,使得片与片之间以及片的内部的颜色过渡很平滑。shading faceted：带有连接线条的曲线。例：x,y=meshgrid(-8:0.5:8);z=sqrt(x.2+y.2);z1=sin(z)./z;surf(x,y,z1)hold on;shading flat;shading interp;axis off,三维图形,3、特殊图形,条形图bar(y)、bar(x,y)、bar(x,y,width)、bar3(y)x为横坐标向量,y为向量或矩阵;y为向量时,每一元素对应的竖条y是m行n列矩阵时,将画出m组竖条,每组包括n个条例：y=5 2 1;9 5 6;8 7 3;5 1 5;4 3 2;subplot(1,2,1);bar(y);subplot(1,2,2);bar3(y);饼图pie(x)、pie(x,explode)x为向量,绘制出各个元素在向量所有元素之和中所占的比例x 为矩阵,绘制出各个元素在矩阵所有元素之和中所占的比例参数explode指定饼图中某些片是否和整个饼图脱开,非零时脱开。例：a=100,150,400,250;pie(a);pie(a,0,0,0,1);colormap(hot),阶梯图、杆图和填充图例：分别用条形图、阶梯图、杆图和填充图形式绘制曲线y=2sin(x).解：绘图程序如下：x=0:pi/10:2*pi;y=2*sin(x);subplot(2,2,1);bar(x,y,g);title(bar(x,y);axis(0,7,-2,2);subplot(2,2,2);stairs(x,y,b);title(stairs(x,y);axis(0,7,-2,2);subplot(2,2,3);stem(x,y,k);title(stem(x,y);axis(0,7,-2,2);subplot(2,2,4);fill(x,y,y);title(fill(x,y);axis(0,7,-2,2);,柱面(cylinder函数)x,y,z=cylinder(r,n)得到柱面的三维坐标cylinder(r,n)绘出柱面r:母线,向量n:旋转圆周上分割线条数例：cylinder(3)cylinder(10,1)t=0:pi/100:4*pi;R=sin(t);cylinder(R,30),球面(sphere函数)x,y,z=sphere(n)得到球面的三维坐标sphere(n)绘出球面n:圆周上分割线条数例：sphere(30);x,y,z=sphere(30);surf(x,y,z),三维图形,多峰函数(peaks函数)格式：Z=peaks(n)用来生成绘图矩阵，矩阵元素由特殊函数在矩形区域-3,3*-3,3的等分网格上的函数值确定 n表示等分的数peaks(n)绘制出多峰函数曲面图 例：peaks(50);x,y,z=peaks(30);surf(x,y,z);,4、其它作图函数,fplot fplot(fname,lims,tol,选项)fplot(sin(x),0,2*pi,*)fplot(sin(x),cos(x),0,2*pi,-1.5,1,5,1e-3,r.)ezplot ezplot(f,x1,x2,y1,y2)在区间x1xx2和y1yy2绘制f(x,y)=0的图形。例：subplot(2,2,1);ezplot(x2+y2-9)subplot(2,2,2);ezplot(cos(tan(pi*x),0,1);subplot(2,2,3);ezplot(cos(tan(pi*x),0,1,0,1);subplot(2,2,4);ezplot(8*cos(t),4*sqrt(2)*sin(t),0,2*pi);,01年数学建模A题,str1=sprintf(W%02d.txt,0:99);str2=reshape(str1,7,100);filename=str2;x=;y=;z=;for i=1:100 ls=load(filename(i,:);ls(1,3)=0;ls(:,3)=2*i;x=x,ls(:,1);y=y,ls(:,2);z=z,ls(:,3);endplot3(x,y,z,.),x=;y=;z=;for i=0:99ff=sprintf(W%02d.txt,i);ls=load(ff);ls(1,3)=0;ls(:,3)=i;x=x,ls(:,1);y=y,ls(:,2);z=z,ls(:,3);endplot3(x,y,z,.),三、数据分析,1、多项式计算2、数据的导入与导出3、数据的统计分析4、插值和拟合,1、多项式计算,多项式可看做是符号表达式,利用符号运算可进行处理多项式的向量表示 设有多项式 表示为:(1)多项式的加法对应系数的加减运算(2)多项式的乘法：conv(P1,P2)(3)多项式除法：Q,r=deconv(P1,P2)其中Q返回相除的商，r返回余式。Q和r仍是多项式的系数向量(4)代数多项式求值：Y=polyval(P,x)若x为一数值，则求多项式P在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求多项式P的值(5)矩阵多项式求值：Y=polyvalm(P,X)X为向方阵，结果为矩阵X的乘与和(6)多项式的求根：x=roots(P)在复数范围内求多项式P的全部根,2、数据的导入与导出,Matlab支持的数据文件Matlab特制的数据文件：.mat文件通用数据文件:文本文件、EXCEL文件、数据库文件 等特制数据文件的打开与保存保存右击需保存的变量，通过菜单操作save命令保存变量到数据文件打开通过菜单“File|open”load 命令通用数据文件通过“File|Import Data”打开导入向导,3、数据的统计分析,求最大值和最小值（max和min）Y,U=max(A)A为向量,将A的最大值存入U,最大值序号存入UA为矩阵,Y表示每列的最大值,U记录每列最大值的行号 U=max(A,B)由A和B中对应元素的较大者组成矩阵U。求和、积、均值、中值、累加和、累乘sum和prod、mean和median、cumsum和cumprod 例：x=1,ones(1,10)*2;y=cumprod(x);sum(y)求得1+2+22+210排序Y,I=sort(A,dim,mode)Y是排序后的矩阵,I记录Y中的元素在A中的位置.dim取1(对列排序);取2(对行排序);默认列排序.mode的取值ascend(升序);descend(降序);默认升序.,向量和矩阵,4、插值与拟合,引例1:在一天24小时内,从零点开始每间隔2小时测得的环境温度数据分别为(度)12,9,9,10,18,24,28,27,25,20,18,15,13.请推测中午1点（即13点）时的温度。引例2:在化学反应中,为研究某化合物的浓度随时间的变化规律,测得一组数据如下：请给出变化规律,由此推测t为1、1.5、2、2.5、10、10.5、11分钟时的值.利用已给出的值,计算相关的值设有一组实验观测数据(xi,yi),i=0,1,n,如何揭示自变量x与因变量y之间的关系？寻找近似的函数关系表达式y=f(x)常用的方法:插值与拟合,插值,设测得的n个点的数据为(xi,yi)(i=1,2,n),构造一个函数y=g(x),使得在xi(i=1,2,n)有g(xi)=yi。g(x)称为插值函数。插值函数自变量的个数：一维插值、二维插值、多维插值等构造插值函数方法:线性插值、多项式插值、样条插值等一维数据的插值Y1=interpl(X,Y,X1,mothod)根据X和Y的值插值，并求插值函数在X1处的值Y是与X等长的向量,Y1是插值函数在X1处的值X1是向量或标量,X1中的元素应在X的范围内method的取值 nearest：最近插值、linear：线性插值 spline：三次样条插值、cubic：三次插值 二维数据插值 Z1=interp2(X,Y,Z,X1,Y1,method)根据X、Y和Z的值插值，并求插值函数在X1、Y1处的值X是行向量、Y是列向量，Z为函数值(size(Z)=length(Y)length(X),引例1的求解,在一天24小时内,从零点开始每间隔2小时测得的环境温度分别为(度)12,9,9,10,18,24,28,27,25,20,18,15,13.请推测中午1点（即13点）时的温度。解：(1)构造数据h=0:2:24;T=12,9,9,10,18,24,28,27,25,20,18,15,13;(2)了解数据的分布绘出数据的图形：plot(h,T,*)(3)选取适当插值方法，计算结果interp1(h,T,13,spline)interp1(h,T,13,cubic)插值函数的评价X=0:24;%选取更多的数据Y1=interp1(h,T,X,spline);Y2=interp1(h,T,X,cubic);subplot(2,1,1);plot(h,T,:,X,Y1);title(spline)subplot(2,1,2);plot(h,T,:,X,Y2,r);title(cubic),例：某实验对一根长为10米的钢轨进行热源的温度传播测试。如下表，其中x表示测量点，h表示测量时间，T表示测得的温度。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度。解：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;x1=0:10;h1=0:20:60;T1=interp2(x,h,T,x1,h1),拟合,根据一组数据(xi,yi)(i=1,2,n)，要求确定一个函数y=f(x)，使这些数据点与曲线y=f(x)总体来说尽量接近，称为曲线拟合。拟合的原理：最小二乘法 最小若拟合函数f(x)是一个多项式,称为多项式拟合.Matlab中多项式拟合命令：P=polyfit(X,Y,m)求数据X与Y的m阶拟合多项式P为拟合多项式的系数,长为m+1的向量X与Y是等长的向量通过polyval求拟合函数在自变量处的值,引例2的求解,引例2:在化学反应中,为研究某化合物的浓度随时间的变化规律,测得一组数据如下：请给出变化规律,推测t为1、1.5、2、2.5、10、10.5、11分钟时的值.解：(1)构造数据t=1:10;y=4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2;(2)选取拟合多项式的阶plot(t,y,o)(3)构造拟合多项式 p=polyfit(t,y,2)p=-0.1109 1.7922 2.9467(4)计算相关的值 t1=1:0.5:11;y1=polyval(p,t1);,拟合函数的评价选取不同的阶做拟合根据散点图直观判断计算最小二乘指标例t=1:10;y=4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2;p=polyfit(t,y,2);y2=polyval(p,t);plot(t,y,:o,t,y2,-*)sum(y-y2).2)类似做3次、4次多项式拟合，从中进行比较,例：在彩色显影中，由经验得知，形成染料的光学密度与析出银的光学密度由公式确定，试验测得如下一批数据：求y关于x的拟合函数。解：由给定的经验公式来求拟合函数，无法直接用polyfit函数，我们换一种思路，将经验公式两边取对数，得 令Y=lny，X=1/x只需求Y关于X的线性拟合x=0.05,0.06,0.07,0.10,0.14,0.20,0.25,0.31,0.38,0.43,0.47;y=0.10,0.14,0.23,0.37,0.59,0.79,1.00,1.12,1.19,1.25,1.29;X=1./x;Y=log(y)P=polyfit(X,Y,1)P=-0.1459 0.5476则Y=-0.1459X+0.5476y=exp(Y)=exp(0.5476)e(-0.1459X)=exp(0.5476)e(-0.1459/x),整理即可,四、方程求解与优化问题求解,方程求解多项式方程求解线性方程组求解 利用左除运算符 利用矩阵的分解 非线性方程求解 单个自变量的非线性方程求解 非线性方程组求解 微分方程求解,线性方程组求解,若有线性方程组利用左除运算符A=2,1,-1,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;x=Ab 利用矩阵的分解LU分解：将矩阵表示为一个下三角矩阵与一个上三角矩阵的乘积。L,U=lu(X):产生L和U,使得X=LU。A=2,1,-1,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;L,U=lu(A);x=U(Lb)QR分解:是将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积.：Q,R=qr(X):产生Q和R,使得X=QR。,非线性方程求解(近似解),例：考察函数作图观察其根的分布尝试用solve求解单个自变量的非线性方程求解 z=fzero(fun,x0)z=fzero(fname,x0)在x0附近寻找函数fun的近似根(fname是待求根的函数名)求上述函数在区间-3,1内的根（1）画图确定根的范围 f=abs(x*sin(x)-exp(x)-1;fplot(f,-3,1)观察图形得出在x=0、-1.2、-2.6附近有根（2）调用fzero函数分别求出其根 fzero(f,0)fzero(f,-1.2)fzero(f,-2.6),上述问题也可这样求解(1)定义函数文件 funx.m,function y=funx(x)y=abs(x.*sin(x)-exp(x)-1;(2)画图确定根的位置 fplot(funx,-3,1)%funx为函数funx的句炳观察图形得到出在x=0、-1.2、-2.6附近有根(3)调用fzero函数分别求出其根 fzero(funx,0)fzero(funx,-1.2)fzero(funx,-2.6)上述问题还可这样求解(1)定义内联（匿名）函数ff=(x)abs(x.*sin(x)-exp(x)-1;(2)画图确定根的位置 fplot(ff,-3,1)(3)调用fzero函数分别求出其根 fzero(ff,-2.5)fzero(ff,-1.5)fzero(ff,0),非线性方程组求解,求解方程组 初始值为-5,5MATLAB的优化工具箱（Optimization Toolbox）对于非线性方程组F(x)=0,求解命令为：X=fsolve(fname,X0,options)解(1)定义函数文件myfun.m function f=myFun(x)f(1)=2*x(1)-x(2)-exp(-x(1);f(2)=-x(1)+2*x(2)-exp(-x(2);(2)提交给求解方程的函数-5,-5%初始值 op=optimset(display,off);%修改参数,中间结果不显示 x=fsolve(myfun,x0,op)x=0.5671 0.5671,微分方程求解,求解微分方程：MATLAB提供了多个求常微分方程数值解的函数,格式为:t,y=solver(fname,tspan,y0)solver为求常微分方程数值解的函数t和y分别给出时间向量和相应的状态向量。fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan为求解区间形式为t0,tf,y0是初始状态列向量。解：建立函数文件funt.mfunction yp=funt(t,y)yp=(y2-t-2)/4/(t+1);求解微分方程t,y=ode23(funt,0,10,2);%求数值解,优化问题求解,1、无约束问题求解（求极值、求最值）模型:minxf(x)求一组x(x=x1,x2,xnT)使得目标函数f(x)为最小2、有约束问题求解模型：求一组x(x=x1,x2,xnT)使得目标函数f(x)为最小,且x满足约束条件G(x)0.约束条件可表示为：线性不等式约束：Ax b线性等式约束：Aeqx=beq非线性不等式约束：Cx b非线性等式约束：Ceqx=0X的上界和下界：Lbnd x Ubng 3、线性规划问题求解线性约束条件下线性目标函数的极值问题,1、无约束问题求解（求极值、最值）,若函数可导,换转为求导数为零的点。例:求y=x*sin(x)-ex在区间-3,2内的极值点。解:fplot(x*sin(x)-exp(x),-3,2)%图形显示在x=-2附近有一极大值点 syms x f=x*sin(x)-exp(x);g=diff(f)%求导函数 sin(x)+x*cos(x)-exp(x)fzero(sin(x)+x*cos(x)-exp(x),-2)结果：-2.0745故在x-2.0745处函数取得极大值。上述问题也求解如下：syms x f=x*sin(x)-exp(x);fplot(char(f),-3,2);%将符号表达式转换为字符串 x0=fzero(char(diff(f),-2);y0=subs(f,x,x0);%将符号表达式f中的x替换为x0,求得极大值.x0,y0=-2.0745,1.6912,Matlab专门提供了求极小值的函数，格式为：(1)x,fval=fminbnd(fname,x1,x2,option)求一元函数在区间(x1,x2)中的极小值点x和极小值。(2)x,fval=fminsearch(fname,x0,option)用单纯形法求多元函数在x0附近的极小值点x和极小值.(3)x,fval=fminunc(fname,x0,option)用拟牛顿法求多元函数在x0附近的极小值点x和极小值.例:求y=x*sin(x)-ex在区间-3,2内的极值点。(1)建立函数文件myfin.m,命令如下：function y=myfin(x)y=-(x.*sin(x)-exp(x);(2)调用fminbnd函数求极大值,命令如下：x,fval=fminbnd(myfin,-3,-1.5)结果：x=-2.0745，fval=-1.6912y=x*sin(x)-ex在x=-2.0745处取得极大值1.6912。等价于f=(x)-(x*sin(x)-exp(x);x,fval=fminbnd(f,-3,-1.5),2、有约束问题求解,模型：求一组x(x=x1,x2,xnT)使得目标函数f(x)为最小,且x满足约束条件G(x)0.约束条件可表示为：线性不等式约束：Ax b线性等式约束：Aeqx=beq非线性不等式约束：Cx b,非线性等式约束：Ceqx=0X的上界和下界：Lbnd x Ubng Matlab优化工具箱提供了解各种约束下的最优化问题。格式：x,fval=fmincon(fname,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,options)求目标函数在x0附近满足约束条件的极小值点和极小值A,b：构成线性不等式约束Aeq,beq：线性等式约束Lbnd,Ubnd：x的上下界常用：x,fval=fmincon(fname,x0,A,b)x,fval=fmincon(fname,x0,A,b,Lbnd,Ubnd),例:设求解有约束最优化问题：解:(1)编写目标函数文件fop.mat,命令如下：function f=fop(x)f=0.4*x(2)+x(1)2+x(2)2-x(1)*x(2)+1/30*x(1)3;(2)设定约束条件，并调用fmincon求解。x0=0.5;0.5;%列向量A=-1,-0.5;-0.5,-1;b=-0.4;-0.5;lb=0;0;%列向量options=optimset(Display,off);x,f=fmincon(fop,x0,A,b,lb,options),3、线性规划问题求解,线性规划问题是研究线性约束条件下线性目标函数极值问题的数学理论和方法。线性规划问题的标准形式为：在Matlab中求解线性规划问题使用函数linprog，格式为 x,fval=linprog(f,A,b,Aeq,beq,lb,ub)x为最优解，fval为目标函数的最优值，f为目标函数的系数向量。例：求解线性规划问题：解：f=2;1;A=-3,-1;-4,-3;-1,-2;b=-3;-6;-2;lb=0;0;x,fval=linprog(f,A,b,lb)x=0.6000 1.2000fval=2.4000,07年B题,原始数据的处理每一条线路均有上、下行，可看作两条线路每行是一条线路：线路、计算方式、各站名Matlab中以表xl表示过某站的所有线路%aa=load(xl.xls)function r,y=allxl(zhan)global aa;bb=strcmp(aa,zhan);r,c=find(bb);r=sort(r);y=aa(r,1);,直达,%aa=load(xl.xls)function y=z0(zhan1,zhan2)global aa;r1=allxl(zhan1);r2=allxl(zhan2);