数学建模算法的matlab代码

/一，matlab绘图函数汇总基本绘图和图形\o"box"box坐标轴边界\o"errorbar"errorbar沿曲线绘制误差条\o"hold"hold在图形窗口中保留当前图形\o"line"line创建线条对象\o"LineSpec(LineSpecification)"LineSpec(LineSpecification)线条规格字符串语法\o"loglog"loglog对数-对数刻度图\o"plot"plot二维线条图\o"plot3"plot3三维线条图\o"plotyy"plotyyy轴分居左右两侧的线条图\o"polar"polar极坐标图\o"semilogx,semilogy"semilogx,semilogy半对数坐标图\o"subplot"subplot在窗口的平铺位置创建坐标轴绘图工具\o"figurepalette"figurepalette显示或隐藏图形窗口的调色板\o"pan"pan交互式移动图像以多方向浏览\o"plotbrowser"plotbrowser显示或隐藏窗口的图形浏览器\o"plotedit"plotedit交互式编辑和标注图形\o"plottools"plottools显示或隐藏图形工具\o"propertyeditor"propertyeditor显示或隐藏属性编辑器\o"rotate3d"rotate3d使用鼠标旋转三维视图\o"showplottool"showplottool显示或隐藏窗口的图形工具\o"zoom"zoom图形的放大、缩小或按比例缩放标注图形\o"annotation"annotation创建注释对象\o"clabel"clabel等高线高程标签\o"datacursormode"datacursormode使能或禁止交互式数据光标模式\o"datetick"datetick日期格式的刻度标签\o"gtext"gtext在二维视图中利用鼠标放置文本\o"legend"legend线条和补片对象的图例\o"rectangle"rectangle创建二维矩形对象\o"texlabel"texlabel产生Tex格式的字符串\o"title"title为当前坐标轴添加标题\o"xlabel,ylabel,zlabel"xlabel,ylabel,zlabelX，Y和Z轴的标签专业绘图(Area、条形图、圆饼图)\o"area"area填充区域的二维图形\o"bar,barh"bar,barh绘制条形图(垂直和水平)\o"bar3,bar3h"bar3,bar3h绘制三维条形图\o"pareto"pareto帕累托图表\o"pie"pie饼图\o"pie3"pie3三维饼图专业绘图(等高线图)\o"contour"contour等高线图矩阵\o"contour3"contour3三维等高线图\o"contourc"contourc低层次的等高线图计算\o"contourf"contourf填充二维等高线图\o"ezcontour"ezcontour易于使用的轮廓绘图仪\o"ezcontourf"ezcontourf易于使用的填充轮廓绘图仪专业绘图(方向和速度图)\o"comet"comet二维彗星图\o"comet3"comet3三维彗星图\o"compass"compass绘制从原点射出的箭头\o"feather"feather绘制速度矢量图\o"quiver"quiver抖动或速度图\o"quiver3"quiver3三维抖动或速度图专业绘图(离散数据图)\o"stairs"stairs阶梯图\o"stem"stem情节离散序列数据\o"stem3"stem3情节3-D离散序列数据专业绘图(函数绘图)\o"ezcontour"ezcontour易用的等高线绘图仪\o"ezcontourf"ezcontourf易用的填充的等高线绘图仪\o"ezmesh"ezmesh易用的三维网格绘图仪\o"ezmeshc"ezmeshc易用的组合式的网格/等高线绘图仪\o"ezplot"ezplot易用的函数绘图仪\o"ezplot3"ezplot3易用的三维参数化曲线绘图仪\o"ezpolar"ezpolar易用的极坐标绘图仪\o"ezsurf"ezsurf易用的三维彩色绘图仪\o"ezsurfc"ezsurfc易用的曲面/等高线绘图仪\o"fplot"fplot在指定的坐标范围内绘制图形专业绘图(直方图)\o"hist"hist直方图阴谋\o"histc"histc直方图计数\o"rose"rose角直方图阴谋专业绘图(多边形和曲面)\o"cylinder"cylinder生成缸\o"delaunay"delaunayDelaunay三角网\o"delaunay3"delaunay3三维Delaunay镶嵌\o"delaunayn"delaunaynN维Delaunay镶嵌\o"dsearch"dsearch搜索Delaunay三角网的最近点\o"ellipsoid"ellipsoid生成椭球\o"fill"fill填充二维多边形\o"fill3"fill3填充的3-D多边形\o"inpolygon"inpolygon多边形区域内点\o"pcolor"pcolor伪(棋盘)阴谋\o"polyarea"polyarea多边形面积\o"rectint"rectint矩形的交集区\o"ribbon"ribbon丝带阴谋\o"slice"slice容积片情节\o"sphere"sphere生成领域\o"waterfall"waterfall瀑布图专业绘图(散点/气泡图)\o"plotmatrix"plotmatrix散点图矩阵\o"scatter"scatter散点图\o"scatter3"scatter3三维散点图专业绘图(动画)\o"frame2im"frame2im返回图像数据及电影相关的框架\o"getframe"getframe电影帧捕获\o"im2frame"im2frame帧图像转换到电影\o"movie"movie播放录制的电影帧\o"noanimate"noanimate所有对象的变化EraseMode正常位图图像\o"frame2im"frame2im返回图像数据及电影相关的框架\o"im2frame"im2frame帧图像转换到电影\o"im2java"im2java图像转换到Java形象\o"image"image显示图像对象\o"imagesc"imagesc规模的数据和显示图像对象\o"imfinfo"imfinfo信息图形文件\o"imformats"imformats管理图像文件格式的注册表\o"imread"imread阅读图像从图形文件\o"imwrite"imwrite写入图像图形文件\o"ind2rgb"ind2rgb索引图像转换到RGB图像打印\o"hgexport"hgexport出口数字\o"orient"orient印刷纸张方向\o"print,printopt"print,printopt打印数字或保存到文件和打印机的默认配置\o"printdlg"printdlg打印对话框\o"printpreview"printpreview打印预览图\o"saveas"saveas保存数字或Simulink框图使用指定的格式句柄图形(查找和辨识图形对象)\o"allchild"allchild查找指定对象的所有儿童\o"ancestor"ancestor祖先的图形对象\o"copyobj"copyobj复制的图形对象和他们的后代\o"delete"delete删除文件或图形对象\o"findall"findall查找所有的图形对象\o"findfigs"findfigs查找离屏数字可见\o"findobj"findobj找到具有特殊性能的图形对象\o"gca"gca当前轴手柄\o"gcbf"gcbf处理数字包含对象，其回调是执行\o"gcbo"gcbo处理的对象，其回调是执行\o"gco"gco处理当前对象\o"get"get查询处理图形对象的属性\o"ishandle"ishandle确定是否输入是有效的句柄图形处理\o"propedit"propedit打开属性编辑器\o"set"set处理图形对象的属性设置句柄图形(对象创建函数)\o"axes"axes创建轴图形对象\o"figure"figure创建数字图形对象\o"hggroup"hggroup创建hggroup对象\o"hgtransform"hgtransform创建hgtransform图形对象\o"image"image显示图像对象\o"light"light根据创建对象\o"line"line创建线对象\o"patch"patch创建一个或多个填充多边形\o"rectangle"rectangle创建二维矩形对象\o"rootobject"rootobject根\o"surface"surface创建面对象\o"text"text创建文本对象在当前轴\o"uicontextmenu"uicontextmenu创建上下文菜单句柄图形(绘图对象)\o"AnnotationArrowProperties"AnnotationArrowProperties注释定义箭头性能\o"AnnotationDoublearrowProperties"AnnotationDoublearrowProperties定义注释doublearrow性能\o"AnnotationEllipseProperties"AnnotationEllipseProperties椭圆定义批注的属性\o"AnnotationLineProperties"AnnotationLineProperties注释行属性定义\o"AnnotationRectangleProperties"AnnotationRectangleProperties注释定义矩形属性\o"AnnotationTextarrowProperties"AnnotationTextarrowProperties定义注释textarrow性能\o"AnnotationTextboxProperties"AnnotationTextboxProperties定义注释文本框属性\o"AreaseriesProperties"AreaseriesProperties属性定义areaseries\o"BarseriesProperties"BarseriesProperties属性定义barseries\o"ContourgroupProperties"ContourgroupProperties属性定义contourgroup\o"ErrorbarseriesProperties"ErrorbarseriesProperties属性定义errorbarseries\o"ImageProperties"ImageProperties定义图像属性\o"LineseriesProperties"LineseriesProperties属性定义lineseries\o"QuivergroupProperties"QuivergroupProperties属性定义quivergroup\o"ScattergroupProperties"ScattergroupProperties属性定义scattergroup\o"StairseriesProperties"StairseriesProperties属性定义stairseries\o"StemseriesProperties"StemseriesProperties属性定义stemseries\o"SurfaceplotProperties"SurfaceplotProperties属性定义surfaceplot句柄图形(图形窗口)\o"clf"clf清除目前的数字窗口\o"close"close删除指定的数字\o"closereq"closereq默认关闭请求函数图\o"drawnow"drawnow刷新事件队列和更新的数字窗口\o"gcf"gcf目前的数字处理\o"hgload"hgload负载处理图形对象的层次结构从文件\o"hgsave"hgsave保存对象的层次结构来处理图形文件\o"newplot"newplot确定要绘制图形对象\o"opengl"opengl控制OpenGL渲染\o"refresh"refresh目前的数字重绘\o"saveas"saveas保存数字或Simulink框图使用指定的格式句柄图形(坐标轴操作)\o"axis"axis轴缩放和外观\o"box"box轴边界\o"cla"cla清除当前轴\o"gca"gca当前轴手柄\o"grid"grid电网线2-D和3-D图\o"ishold"ishold当前保持状态\o"makehgtform"makehgtform创建4×4变换矩阵句柄图形(对象属性操作)\o"get"get查询处理图形对象的属性\o"linkaxes"linkaxes同步限制指定的2-D轴\o"linkprop"linkprop保持相应的属性值相同\o"refreshdata"refreshdata刷新数据图表数据源时指定\o"set"set处理图形对象的属性设置二，hamiton回路算法提供一种求解最优哈密尔顿的算法三边交换调整法，要求在运行jiaohuan3（三交换法）之前，给定邻接矩阵C和节点个数N，结果路径存放于R中。

bianquan.m文件给出了一个参数实例，可在命令窗口中输入bianquan,得到邻接矩阵C和节点个数N以及一个任意给出的路径R,，回车后再输入jiaohuan3,得到了最优解。由于没有经过大量的实验，又是近似算法，对于网络比较复杂的情况，可以尝试多运行几次jiaohuan3,看是否能到进一步的优化结果。%%%%%%bianquan.m%%%%%%%N=13;fori=1:Nforj=1:NC(i,j)=inf;endendfori=1:NC(i,i)=0;endC(1,2)=6.0;C(1,13)=12.9;C(2,3)=5.9;C(2,4)=10.3;C(3,4)=12.2;C(3,5)=17.6;C(4,13)=8.8;C(4,7)=7.4;C(4,5)=11.5;C(5,2)=17.6;C(5,6)=8.2;C(6,9)=14.9;C(6,7)=20.3;C(7,9)=19.0;C(7,8)=7.3;C(8,9)=8.1;C(8,13)=9.2;C(9,10)=10.3;C(10,11)=7.7;C(11,12)=7.2;C(12,13)=7.9;fori=1:Nforj=1:NifC(i,j)<infC(j,i)=C(i,j);endendendfori=1:NC(i,i)=0;endR=[47653211312111098];<prename="code"class="plain">%%%%%%%%jiaohuan3.m%%%%%%%%%%n=0;forI=1:(N-2)forJ=(I+1):(N-1)forK=(J+1):Nn=n+1;Z(n,:)=[IJK];endendendR=1:Nform=1:(N*(N-1)*(N-2)/6)I=Z(m,1);J=Z(m,2);K=Z(m,3);r=R;ifJ-I~=1&K-J~=1&K-I~=N-1forq=1:(J-I)r(I+q)=R(J+1-q);endforq=1:(K-J)r(J+q)=R(K+1-q);endendifJ-I==1&K-J==1r(K)=R(J);r(J)=R(K);endifJ-I==1&K-J~=1&K-I~=N-1forq=1:(K-J)r(I+q)=R(I+1+q);endr(K)=R(J);endifK-J==1&J-I~=1&K~=Nforq=1:(J-I)r(I+1+q)=R(I+q);endr(I+1)=R(K);endifI==1&J==2&K==Nforq=1:(N-2)r(1+q)=R(2+q);endr(N)=R(2);endifI==1&J==(N-1)&K==Nforq=1:(N-2)r(q)=R(1+q);endr(N-1)=R(1);endifJ-I~=1&K-I==N-1forq=1:(J-1)r(q)=R(1+q);endr(J)=R(1);endifJ==(N-1)&K==N&J-I~=1r(J+1)=R(N);forq=1:(N-J-1)r(J+1+q)=R(J+q);endendifcost_sum(r,C,N)<cost_sum(R,C,N)R=rendendfprintf('总长为%f\n',cost_sum(R,C,N))%%%%%%cost_sum.m%%%%%%%%functiony=cost_sum(x,C,N)y=0;fori=1:(N-1)y=y+C(x(i),x(i+1));endy=y+C(x(N),x(1));三，灰色预测代码<prename="code"class="plain">clearclcX=[136143165152165181204272319491571605665640628];x1(1)=X(1);X1=[];fori=1:1:14x1(i+1)=x1(i)+X(i+1);X1=[X1,x1(i)];endX1=[X1,X1(14)+X(15)]fork=3:1:15p(k)=X(k)/X1(k-1);p1(k)=X1(k)/X1(k-1);endp,p1clearkZ=[];fork=2:1:15z(k)=0.5*X1(k)+0.5*X1(k-1);Z=[Z,z(k)];endZB=[-Z',ones(14,1)]Y=[];clearifori=2:1:15Y=[Y;X(i)];endYA=inv(B'*B)*B'*Yclearky1=[];fork=1:1:15y(k)=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1);y1=[y1;y(k)];endy1clearkX2=[];fork=2:1:15x2(k)=y1(k)-y1(k-1);X2=[X2;x2(k)];endX2=[y1(1);X2]e=X'-X2m=abs(e)./X's=e'*en=sum(m)/13clearksymsky=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1)Y1=[];forj=16:1:21y11=subs(y,k,j)-subs(y,k,j-1);Y1=[Y1;y11];endY1%程序中的变量定义:alpha是包含α、μ值的矩阵;%ago是预测后累加值矩阵;var是预测值矩阵;%error是残差矩阵;c是后验差比值functionbasicgrey(x,m)%定义函数basicgray(x)ifnargin==1%m为想预测数据的个数，默认为1m=1;endclc;%清屏,以使计算结果独立显示iflength(x(:,1))==1%对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x';endn=length(x);%取输入数据的样本量x1(:,1)=cumsum(x);%计算累加值,并将值赋及矩阵befori=2:n%对原始数列平行移位Y(i-1,:)=x(i,:);endfori=2:n%计算数据矩阵B的第一列数据z(i,1)=0.5*x1(i-1,:)+0.5*x1(i,:);endB=ones(n-1,2);%构造数据矩阵BB(:,1)=-z(2:n,1);alpha=inv(B'*B)*B'*Y;%计算参数α、μ矩阵fori=1:n+m%计算数据估计值的累加数列,如改n+1为n+m可预测后m个值ago(i,:)=(x1(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:);endvar(1,:)=ago(1,:);fori=1:n+m-1%可预测后m个值var(i+1,:)=ago(i+1,:)-ago(i,:);%估计值的累加数列的还原,并计算出下m个预测值end[P,c,error]=lcheck(x,var);%进行后验差检验[rela]=relations([x';var(1:n)']);%关联度检验ago%显示输出预测值的累加数列alpha%显示输出参数α、μ数列var%显示输出预测值error%显示输出误差P%显示计算小残差概率c%显示后验差的比值crela%显示关联度judge(P,c,rela)%评价函数显示这个模型是否合格<prename="code"class="plain">functionjudge(P,c,rela)%评价指标并显示比较结果ifrela>0.6'根据经验关联度检验结果为满意（关联度只是参考主要看后验差的结果）'else'根据经验关联度检验结果为不满意（关联度只是参考主要看后验差的结果）'endifP>0.95&c<0.5'后验差结果显示这个模型评价为“优”'elseifP>0.8&c<0.5'后验差结果显示这个模型评价为“合格”'elseifP>0.7&c<0.65'后验差结果显示这个模型评价为“勉强合格”'else'后验差结果显示这个模型评价为“不合格”'endendendfunction[P,c,error]=lcheck(x,var)%进行后验差检验n=length(x);fori=1:nerror(i,:)=abs(var(i,:)-x(i,:));%计算绝对残差endc=std(abs(error))/std(x);%调用统计工具箱的标准差函数计算后验差的比值cs0=0.6745*std(x);ek=abs(error-mean(error));pk=0;fori=1:nifek(i,:)<s0pk=pk+1;endendP=pk/n;%计算小残差概率%附带的质料里有一部分讲了关联度function[rela]=relations(x)%以x(1,:)的参考序列求关联度[m,n]=size(x);fori=1:mforj=n:-1:2x(i,j)=x(i,j)/x(i,1);endendfori=2:mx(i,:)=abs(x(i,:)-x(1,:));%求序列差endc=x(2:m,:);Max=max(max(c));%求两极差Min=min(min(c));p=0.5;%p称为分辨率,0<p<1,一般取p=0.5fori=1:m-1forj=1:nr(i,j)=(Min+p*Max)/(c(i,j)+p*Max);%计算关联系数endendfori=1:m-1rela(i)=sum(r(i,:))/n;%求关联度end四，非线性拟合functionf=example1(c,tdata)f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata));<prename="code"class="plain">functionf=zhengtai(c,x)f=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c(2)^2));x=1:1:12;y=[00010310128212]';c0=[28];fori=1:1000c=lsqcurvefit(@zhengtai,c0,x,y);c0=c;endy1=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c(2)^2));plot(x,y,'r-',x,y1);legend('实验数据','拟合曲线')x=[0.250.50.7511.522.533.544.55678910111213141516]';y=[3068758282776868585150413835282518151210774]';f=@(c,x)c(1)*(exp(-c(2)*x)-exp(-c(3)*x));c0=[1140.12]';fori=1:50opt=optimset('TolFun',1e-3);[cR]=nlinfit(x,y,f,c0,opt)c0=c;holdonplot(x,c(1)*(exp(-c(2)*x)-exp(-c(3)*x)),'g')endt=[0.250.50.7511.522.533.544.55678910111213141516];y=[3068758282776868585150413835282518151210774];c0=[111];fori=1:50c=lsqcurvefit(@example1,c0,t,y);c0=c;endy1=c(1)*(exp(-c(2)*t)-exp(-c(3)*t));plot(t,y,'+',t,y1);legend('实验数据','拟合曲线')五，插值拟合相关知识在生产和科学实验中，自变量及因变量间的函数关系有时不能写出解析表达式，而只能得到函数在若干点的函数值或导数值，或者表达式过于复杂而需要较大的计算量。当要求知道其它点的函数值时，需要估计函数值在该点的值。为了完成这样的任务，需要构造一个比较简单的函数，使函数在观测点的值等于已知的值，或使函数在该点的导数值等于已知的值，寻找这样的函数有很多方法。根据测量数据的类型有以下两类处理观测数据的方法。（1）测量值是准确的，没有误差，一般用插值。（2）测量值及真实值有误差，一般用曲线拟合。在MATLAB中，无论是插值还是拟合，都有相应的函数来处理。一、插值1、一维插值：已知离散点上的数据集，即已知在点集X=上的函数值Y=，构造一个解析函数（其图形为一曲线）通过这些点，并能够求出这些点之间的值，这一过程称为一维插值。MATLAB命令：yi=interp1(X,Y,xi,method)该命令用指定的算法找出一个一元函数，然后以给出处的值。xi可以是一个标量，也可以是一个向量，是向量时，必须单调，method可以下列方法之一：‘nearest’：最近邻点插值，直接完成计算；‘spline’：三次样条函数插值；‘linear’：线性插值（缺省方式），直接完成计算；‘cubic’：三次函数插值；对于[min{xi},max{xi}]外的值，MATLAB使用外推的方法计算数值。例1：已知某产品从1900年到2010年每隔10年的产量为：75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,249.633,256.344,267.893，计算出1995年的产量，用三次样条插值的方法，画出每隔一年的插值曲线图形，同时将原始的数据画在同一图上。解：程序如下[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?year=1900:10:2010;product=[75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,249.633,256.344,267.893]p1995=interp1(year,product,1995)x=1900:2010;y=interp1(year,product,x,'cubic');plot(year,product,'o',x,y);计算结果为：p1995=252.9885。2、二维插值已知离散点上的数据集，即已知在点集上的函数值，构造一个解析函数（其图形为一曲面）通过这些点，并能够求出这些已知点以外的点的函数值，这一过程称为二维插值。MATLAB函数：Zi=interp2(X,Y,Z,Xi,Yi,method)该命令用指定的算法找出一个二元函数，然后以给出处的值。返回数据矩阵，Xi，Yi是向量，且必须单调，和meshgrid(Xi,Yi)是同类型的。method可以下列方法之一：‘nearest’：最近邻点插值，直接完成计算；‘spline’：三次样条函数插值；‘linear’：线性插值（缺省方式），直接完成计算；‘cubic’：三次函数插值；例2：已知1950年到1990年间每隔10年，服务年限从10年到30年每隔10年的劳动报酬表如下：表：某企业工作人员的月平均工资（元）服务年限年份1020301950150.697169.592187.6521960179.323195.072250.2871970203.212239.092322.7671980226.505273.706426.7301990249.633370.281598.243试计算1975年时，15年工龄的工作人员平均工资。解：程序如下：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?years=1950:10:1990;service=10:10:30;wage=[150.697169.592187.652179.323195.072250.287203.212239.092322.767226.505273.706426.730249.633370.281598.243];mesh(service,years,wage)%绘原始数据图w=interp2(service,years,wage,15,1975);%求点(15,1975)处的值计算结果为：235.6288例3：设有数据x=1,2,3,4,5,6，y=1,2,3,4，在由x,y构成的网格上，数据为：12,10,11,11,13,1516,22,28,35,27,2018,21,26,32,28,2520,25,30,33,32,20求通过这些点的插值曲面。解：程序为：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?x=1:6;y=1:4;t=[12,10,11,11,13,1516,22,28,35,27,2018,21,26,32,28,25;20,25,30,33,32,20]subplot(1,2,1)mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;[x2,y2]=meshgrid(x1,y1);t1=interp2(x,y,t,x2,y2,'cubic');subplot(1,2,2)mesh(x1,y1,t1);结果如右图。作业：已知某处山区地形选点测量坐标数据为：x=00.511.522.533.544.55y=00.511.522.533.544.555.56海拔高度数据为：z=8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287

1、画出原始数据图;2、画出加密后的地貌图，并在图中标出原始数据二、拟合曲线拟合已知离散点上的数据集，即已知在点集上的函数值，构造一个解析函数（其图形为一曲线）使在原离散点上尽可能接近给定的值，这一过程称为曲线拟合。最常用的曲线拟合方法是最小二乘法，该方法是寻找函数使得最小。MATLAB函数：p=polyfit(x,y,n)[p,s]=polyfit(x,y,n)说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s用于生成预测值的误差估计。(见下一函数polyval)多项式曲线求值函数：polyval()调用格式：y=polyval(p,x)[y,DELTA]=polyval(p,x,s)说明：y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。[y,DELTA]=polyval(p,x,s)使用polyfit函数的选项输出s得出误差估计YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则YDELTA将至少包含50%的预测值。例5：求如下给定数据的拟合曲线，x=[0.5,1.0,1.5,2.0,2.5,3.0]，y=[1.75,2.45,3.81,4.80,7.00,8.60]。解：MATLAB程序如下：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?x=[0.5,1.0,1.5,2.0,2.5,3.0];y=[1.75,2.45,3.81,4.80,7.00,8.60];p=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'-b')计算结果为：p=0.56140.82871.1560此结果表示拟合函数为：，用此函数拟合数据的效果如图所示。例2：由离散数据x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52拟合出多项式。程序：x=0:.1:1;y=[.3.511.41.61.9.6.4.81.52]n=3;p=polyfit(x,y,n)xi=linspace(0,1,100);z=polyval(p,xi);%多项式求值plot(x,y,’o’,xi,z,’k:’,x,y,’b’)legend(‘原始数据’,’3阶曲线’)结果：p=16.7832-25.745910.9802-0.0035多项式为：16.7832x3-25.7459x2+10.9802x-0.0035曲线拟合图形：也可由函数给出数据。例3：x=1:20,y=x+3*sin(x)程序：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?x=1:20;y=x+3*sin(x);p=polyfit(x,y,6)xi=linspace(1,20,100);z=polyval(p,xi);%¶àÏîʽÇóÖµº¯Êýplot(x,y,'o',xi,z,'k:',x,y,'b')结果：p=0.0000-0.00210.0505-0.59713.6472-9.729511.3304

再用10阶多项式拟合程序：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?x=1:20;y=x+3*sin(x);p=polyfit(x,y,10)xi=linspace(1,20,100);z=polyval(p,xi);plot(x,y,'o',xi,z,'k:',x,y,'b')结果：p=Columns1through70.0000-0.00000.0004-0.01140.1814-1.806511.2360Columns8through11-42.086188.5907-92.815540.267

可用不同阶的多项式来拟合数据，但也不是阶数越高拟合的越好。六，层次分析法disp('请输入判断矩阵A(n阶)');A=input('A=');[n,n]=size(A);x=ones(n,100);y=ones(n,100);m=zeros(1,100);m(1)=max(x(:,1));y(:,1)=x(:,1);x(:,2)=A*y(:,1);m(2)=max(x(:,2));y(:,2)=x(:,2)/m(2);p=0.0001;i=2;k=abs(m(2)-m(1));whilek>pi=i+1;x(:,i)=A*y(:,i-1);m(i)=max(x(:,i));y(:,i)=x(:,i)/m(i);k=abs(m(i)-m(i-1));enda=sum(y(:,i));w=y(:,i)/a;t=m(i);disp(w);disp(t);%以下是一致性检验CI=(t-n)/(n-1);RI=[000.520.891.121.261.361.411.461.491.521.541.561.581.59];CR=CI/RI(n);ifCR<0.10disp('此矩阵的一致性可以接受!');disp('CI=');disp(CI);disp('CR=');disp(CR);end七，参数估计文章的主要思想来源于Matlab|Simulink仿真世界的一篇类似的文章。我这里把这个思想引入到我们的体系来，并以一个新的例子讲解这一用法。fminsearch用来求解多维无约束的非线性优化问题，它的基本形式是：[X,FVAL,EXITFLAG,OUTPUT]=FMINSEARCH(FUN,X0,OPTIONS).大段的Matlab帮助文档我就不翻译解释了，有兴趣的朋友可以参见Matlab联机帮助，我这里只介绍他在参数估计中的作用。在参数估计中经常用到正态分布的参数估计。在matlab系统中有一个函数叫做normfit就直接可以完成这样的参数估计，返回均值mu和均方差sigma的估计，但是这里有一个要求，就是它的输入信息必须是随机的数字序列。如得到1000个服从正态分布的随机数向量R，用命令[phatpci]=normfit(R)，就可以得到参数估计了。然而如果我我们得知某些已经处于pdf函数曲线上的点时，这时需要对函数进行拟合运算。估计参数的原理是从已知的一序列数据中，对于给定的任何一组参数，计算用其估计数据得到的方差，然后利用fminsearch函数求当方差满足最小的时候的参数，这就是需要估计的参数。来看一下下面的列子：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?smu=10,ssig=25;%假设原来均值方差分别为：10，25R=randn(1000,1)*ssig+smu;%生成满足要求的1000个随机数[yx]=myhist(R);%生成统计信息，x，y分别表示分组中值序列和落入该组的统计数目bar(x,y)%绘制直方图holdonplot(x,y,'ro')%绘制对应点[pmsmse]=normpdffit(x,y,8,20);%根据得到的统计信息x,y对其进行参数估计，8，20分别代表均值和方差的初值t=min(x):(max(x)-min(x))/200:max(x);%定义绘图区间ny=normpdf(t,smu,ssig);%真实分布曲线数据nyf=normpdf(t,pms(1),pms(2));%拟合分布曲线数据plot(t,ny,'r-')plot(t,nyf,'b-.')legend('hist','histvalue','turepdf','fitpdf')%绘制两条曲线作对比上面例子中所用的几个函数定义如下：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?function[hxout]=myhist(data,nbins)%用于统计信息，生成和pdf函数值相同的hist统计方式。ifnargin==1nbins=uint32(1+log(length(data))/log(2));endnbins=double(nbins);data=data(:);[hxout]=hist(data,nbins);ew=xout(2)-xout(1);h=h./(ew*length(data));[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?function[pabmse]=normpdffit(x,y,a0,b0)%正态分布pdf参数估计p=[a0b0];opt=optimset('fminsearch');opt.TolX=0.001;opt.Display='off';[pabmse]=fminsearch(@normpdfse,p,opt,x,y);%这里需要注意，opt参数已经传递给fminsearch，但是对于原计算方差的函数来说，还需要两个参数x,y,这两个参数就写在opt参数的后面，这样可以完成其他参数的传递。%这里说下以前探索的时候的失败经验：用global把参数公有化，然后函数只传递变化的参数（需要估计的参数），但是失败了。所以了解这种参数传递方法是非常有必要的。functionse=normpdfse(pab,X,Y)%计算对于任何一组参数pab(1)，pab(2)，给出当前数据下的方差来。se=var(Y-normpdf(X,pab(1),pab(2)));运行结果如图所示：从图中可以看出，随机数在这里变成了统计信息，统计信息反映到了绘制的点信息上，图中圆圈所示。真实的pdf为红色曲线，估计的曲线为蓝色虚线。从图中可知，估计的效果非常满意。如果在函数中加上：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?disp'theresultofnormfitfunction:'[musg]=normfit(R)disp'theresultoffminsearch:'[pmsmse]=normpdffit(x,y,8,20)得到结果：[plain]\o"viewplain"viewplain\o"copy"copy\o"print"print\o"?"?theresultofnormfitfunction:mu=10.443sg=25.61945417031251theresultoffminsearch:pms=10.38425.32479396733891mse=7.228