数学建模常用Matlab/Lingo/c代码总结系列——floyd最短路径
分类： Matlab 数学建模2011-11-16 23:37 281人阅读 评论(0) 收藏 举报
例 9 某公司在六个城市 c1,c2, …c6 中有分公司，从 i
ci 到 cj的直接航程票价记在下述矩阵的 (I,j) 位置上。（∞表示无直接航路），请帮助该公司设计一张城市c1 到其它城市间的票价最便宜的路线图。
[plain] view plaincopy
1. clc,clear 2.
3. a=zeros(6);
4.
5. a(1,2)=50;a(1,4)=40;a(1,5)=25;a(1,6)=10; 6.
7. a(2,3)=15;a(2,4)=20;a(2,6)=25; 8.
9. a(3,4)=10;a(3,5)=20; 10.
11. a(4,5)=10;a(4,6)=25;
12.
13. a(5,6)=55; 14.
15. a=a+a'; 16.
17. a(find(a==0))=inf; 18.
19. pb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a)); 20.
21. d(1:length(a))=inf;d(1)=0;temp=1; 22.
23. while sum(pb) 27. d(tb)=min(d(tb),d(temp)+a(temp,tb)); 28. 29. tmpb=find(d(tb)==min(d(tb))); 30. 31. temp=tb(tmpb(1)); 32. 33. pb(temp)=1; 34. 35. index1=[index1,temp]; 36. 37. temp2=find(d(index1)==d(temp)-a(temp,index1)); 38. 39. index2(temp)=index1(temp2(1)); 40. 41. end 42. 43. d, index1, index2 编写 LINGO 程序如下： [plain] view plaincopy 1. model: 2. 3. sets: 4. 5. cities/A,B1,B2,C1,C2,C3,D/; 6. 7. roads(cities,cities)/A B1,A B2,B1 C1,B1C2,B1 C3,B2 C1, 8. 9. B2 C2,B2 C3,C1 D,C2 D,C3 D/:w,x; 10. 11. endsets 12. 13. data: 14. 15. w=2 4 3 3 1 2 3 1 1 3 4; 16. 17. enddata 18. 19. n=@size(cities); !城市的个数; 20. 21. min=@sum(roads:w*x); 22. 23. @for(cities(i)|i #ne#1 #and# i #ne#n: 24. 25. @sum(roads(i,j):x(i,j))=@sum(roads(j,i):x(j,i))); 26. 27. @sum(roads(i,j)|i #eq#1:x(i,j))=1; 28. 29. @sum(roads(i,j)|j #eq#n:x(i,j))=1; 30. 31. end [plain] view plaincopy 1. model: 2. 3. sets: 4. 5. cities/1..11/; 6. 7. roads(cities,cities):w,x; 8. 9. endsets 10. 11. data: 12. 13. w=0; 14. 15. enddata 16. 17. calc: 18. 19. w(1,2)=2;w(1,3)=8;w(1,4)=1; 20. 21. w(2,3)=6;w(2,5)=1; 22. 23. w(3,4)=7;w(3,5)=5;w(3,6)=1;w(3,7)=2; 24. 25. w(4,7)=9; 26. 27. w(5,6)=3;w(5,8)=2;w(5,9)=9; 28. 29. w(6,7)=4;w(6,9)=6; 30. 31. w(7,9)=3;w(7,10)=1; 32. 33. w(8,9)=7;w(8,11)=9; 34. 35. w(9,10)=1;w(9,11)=2;w(10,11)=4; 36. 37. @for(roads(i,j):w(i,j)=w(i,j)+w(j,i)); 38. 39. @for(roads(i,j):w(i,j)=@if(w(i,j) #eq# 0,1000,w(i,j))); 40. 41. endcalc 42. 43. n=@size(cities); !城市的个数; 44. 45. min=@sum(roads:w*x); 46. 47. @for(cities(i)|i #ne#1 #and# i #ne# 48. 49. n:@sum(cities(j):x(i,j))=@sum(cities(j):x(j,i))); 50. 51. @sum(cities(j):x(1,j))=1; 52. 53. @sum(cities(j):x(j,1))=0; !不能回到顶点1; 54. 55. @sum(cities(j):x(j,n))=1; 56. 57. @for(roads:@bin(x)); 58. 59. end 60. 61. 例12 用Floyd算法求解例9。 矩阵path用来存放每对顶点之间最短路径上所经过的顶点的序号。Floyd算法的 Matlab程序如下： [plain] view plaincopy 1. clear;clc; 2. 3. n=6; a=zeros(n); 4. 5. a(1,2)=50;a(1,4)=40;a(1,5)=25;a(1,6)=10; 6. 7. a(2,3)=15;a(2,4)=20;a(2,6)=25;a(3,4)=10;a(3,5)=20; 8. 9. a(4,5)=10;a(4,6)=25; a(5,6)=55; 10. 11. a=a+a'; M=max(max(a))*n^2; %M为充分大的正实数 12. 13. a=a+((a==0)-eye(n))*M; 14. 15. path=zeros(n); 16. 17. for k=1:n 18. 19. for i=1:n 20. 21. for j=1:n 22. 23. if a(i,j)>a(i,k)+a(k,j) 24. 25. a(i,j)=a(i,k)+a(k,j); 26. 27. path(i,j)=k; 28. 29. end 30. 31. end 32. 33. end 34. 35. end 36. 37. a, path 我们使用LINGO9.0编写的FLOYD算法如下： [plain] view plaincopy 1. model: 2. 3. sets: nodes/c1..c6/; 4. 5. link(nodes,nodes):w,path; !path标志最短路径上走过的顶点; 6. 7. endsets 8. 9. data: 10. 11. path=0; 12. 13. w=0; 14. 15. @text(mydata1.txt)=@writefor(nodes(i):@writefor(nodes(j): 16. 17. @format(w(i,j),' 10.0f')),@newline(1)); 18. 19. @text(mydata1.txt)=@write(@newline(1)); 20. 21. @text(mydata1.txt)=@writefor(nodes(i):@writefor(nodes(j): 22. 23. @format(path(i,j),' 10.0f')),@newline(1)); 24. 25. enddata 26. 27. calc: 28. 29. w(1,2)=50;w(1,4)=40;w(1,5)=25;w(1,6)=10; 30. 31. w(2,3)=15;w(2,4)=20;w(2,6)=25; 32. 33. w(3,4)=10;w(3,5)=20; 34. 35. w(4,5)=10;w(4,6)=25;w(5,6)=55; 36. 37. @for(link(i,j):w(i,j)=w(i,j)+w(j,i)); 38. 39. @for(link(i,j) |i#ne#j:w(i,j)=@if(w(i,j)#eq#0,10000,w(i,j))); 40. 41. @for(nodes(k):@for(nodes(i):@for(nodes(j): 42. 43. tm=@smin(w(i,j),w(i,k)+w(k,j)); 44. 45. path(i,j)=@if(w(i,j)#gt#tm,k,path(i,j));w(i,j)=tm))); 46. 47. endcalc 48. 49. end 数学建模常用Matlab/Lingo/c代码总结系列——最小费用最大流问题 分类： 数学建模 Matlab2011-11-17 12:56 433人阅读 评论(0) 收藏 举报 例 19（最小费用最大流问题）（续例18）由于输油管道的长短不一或地质等原因， 使每条管道上运输费用也不相同，因此，除考虑输油管道的最大流外，还需要考虑输油 管道输送最大流的最小费用。图 8 所示是带有运费的网络，其中第 1 个数字是网络的容 量，第 2 个数字是网络的单位运费。 图8 最小费用最大流问题 解 按照最小费用流的数学规划写出相应的 LINGO 程序如下： [plain] view plaincopy 1. model: 2. 3. sets: 4. 5. nodes/s,1,2,3,4,t/:d; 6. 7. arcs(nodes,nodes)/s 1,s 3,1 2,1 3,2 3,2 t,34,4 2,4 t/:c,u,f; 8. 9. endsets 10. 11. data: 12. 13. d=14 0 0 0 0 -14; !最大流为14; 14. 15. c=2 8 2 5 1 6 3 4 7; 16. 17. u=8 7 9 5 2 5 9 6 10; 18. 19. enddata 20. 21. min=@sum(arcs:c*f); 22. 23. @for(nodes(i):@sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i))=d(i)); 24. 25. @for(arcs:@bnd(0,f,u)); 26. 27. end 求得最大流的最小费用是 205，而原最大流的费用为 210 单位，原方案并不是最优 的。 类似地，可以利用赋权邻接矩阵编程求得最小费用最大流。LINGO 程序如下： [plain] view plaincopy 1. model: 2. 3. sets: 4. 5. nodes/s,1,2,3,4,t/:d; 6. 7. arcs(nodes,nodes):c,u,f; 8. 9. endsets 10. 11. data: 12. 13. d=14 0 0 0 0 -14; 14. 15. c=0; u=0; 16. 17. enddata 18. 19. calc: 20. 21. c(1,2)=2;c(1,4)=8; 22. 23. c(2,3)=2;c(2,4)=5; 24. 25. c(3,4)=1;c(3,6)=6; 26. 27. c(4,5)=3;c(5,3)=4;c(5,6)=7; 28. 29. u(1,2)=8;u(1,4)=7; 30. 31. u(2,3)=9;u(2,4)=5; 32. 33. u(3,4)=2;u(3,6)=5; 34. 35. u(4,5)=9;u(5,3)=6;u(5,6)=10; 36. 37. endcalc 38. 39. min=@sum(arcs:c*f); 40. 41. @for(nodes(i):@sum(nodes(j):f(i,j))-@sum(nodes(j):f(j,i))=d(i)); 42. 43. @for(arcs:@bnd(0,f,u)); 44. 45. end 求最小费用流的一种方法—迭代法 下面我们编写了最小费用最大流函数 mincostmaxflow，其中调用了利用 Floyd 算法 求最短路的函数 floydpath。 求解例 19 具体程序如下（下面的全部程序放在一个文件中）： [plain] view plaincopy 1. function mainexample19 2. 3. clear;clc; 4. 5. global M num 6. 7. c=zeros(6);u=zeros(6); 8. 9. c(1,2)=2;c(1,4)=8;c(2,3)=2;c(2,4)=5; 10. 11. c(3,4)=1;c(3,6)=6;c(4,5)=3;c(5,3)=4;c(5,6)=7; 12. 13. u(1,2)=8;u(1,4)=7;u(2,3)=9;u(2,4)=5; 14. 15. u(3,4)=2;u(3,6)=5;u(4,5)=9;u(5,3)=6;u(5,6)=10; 16. 17. num=size(u,1);M=sum(sum(u))*num^2; 18. 19. [f,val]=mincostmaxflow(u,c) 20. 21. 22. 23. %求最短路径函数 24. 25. function path=floydpath(w); 26. 27. global M num 28. 29. w=w+((w==0)-eye(num))*M; 30. 31. p=zeros(num); 32. 33. for k=1:num 34. 35. for i=1:num 36. 37. for j=1:num 38. 39. if w(i,j)>w(i,k)+w(k,j) 40. 41. w(i,j)=w(i,k)+w(k,j); 42. 43. p(i,j)=k; 44. 45. end 46. 47. end 48. 49. end 50. 51. end 52. 53. if w(1,num) ==M 54. 55. path=[]; 56. 57. else 58. 59. path=zeros(num); 60. 61. s=1;t=num;m=p(s,t); 62. 63. while ~isempty(m) 64. 65. if m(1) 66. 67. s=[s,m(1)];t=[t,t(1)];t(1)=m(1); 68. 69. m(1)=[];m=[p(s(1),t(1)),m,p(s(end),t(end))]; 70. 71. else 72. 73. path(s(1),t(1))=1;s(1)=[];m(1)=[];t(1)=[]; 74. 75. end 76. 77. end 78. 79. end 80. 81. [plain] view plaincopy 1. %最小费用最大流函数 2. 3. function[flow,val]=mincostmaxflow(rongliang,cost,flowvalue); 4. 5. %第一个参数：容量矩阵；第二个参数：费用矩阵； 6. 7. %前两个参数必须在不通路处置零 8. 9. %第三个参数：指定容量值（可以不写，表示求最小费用最大流） 10. 11. %返回值 flow 为可行流矩阵,val 为最小费用值 12. 13. global M 14. 15. flow=zeros(size(rongliang));allflow=sum(flow(1,:)); 16. 17. if nargin<3 18. 19. flowvalue=M; 20. 21. end 22. 23. while allflow 27. path=floydpath(w);%调用 floydpath 函数 28. 29. if isempty(path) 30. 31. val=sum(sum(flow.*cost)); 32. 33. return; 34. 35. end 36. 37. theta=min(min(path.*(rongliang-flow)+(path.*(rongliang-flow)==0).*M)); 38. 39. theta=min([min(path'.*flow+(path'.*flow==0).*M),theta]); 40. 41. flow=flow+(rongliang>0).*(path-path').*theta; 42. 43. allflow=sum(flow(1,:)); 44. 45. end 46. 47. val=sum(sum(flow.*cost)); 数学建模常用Matlab/Lingo/c代码总结系列——整数规划问题 分类： 数学建模2011-11-17 12:59 109人阅读 评论(0) 收藏 举报 [plain] view plaincopy 1. 2. 3. model: 4. sets: 5. row/1..4/:b; 6. col/1..5/:c1,c2,x; 7. link(row,col):a; 8. endsets 9. data: 10. c1=1,1,3,4,2; 11. c2=-8,-2,-3,-1,-2; 12. a=1 1 1 1 1 13. 1 2 2 1 6 14. 2 1 6 0 0 15. 0 0 1 1 5; 16. b=400,800,200,200; 17. enddata 18. max=@sum(col:c1*x^2+c2*x); 19. @for(row(i):@sum(col(j):a(i,j)*x(j)) 数学建模常用Matlab/Lingo/c代码总结系列——旅行商TSP问题 分类： 数学建模 Matlab2011-11-17 13:06 182人阅读 评论(0) 收藏 举报 [plain] view plaincopy 1. Lingo代码： 2. MODEL: 3. 4. SETS: 5. CITY / 1.. 6/: U; ! U( I) = sequence no. of city; 6. LINK( CITY, CITY): 7. DIST, ! The distance matrix; 8. X; ! X( I, J) = 1 if we use link I, J; 9. ENDSETS 10. DATA: !Distance matrix, it need not be symmetric; 11. DIST =0 56 35 21 51 60 12. 56 0 21 57 78 70 13. 35 21 0 36 68 68 14. 21 57 36 0 51 61 15. 51 78 68 51 0 13 16. 60 70 68 61 13 0; 17. ENDDATA 18. !The model:Ref. Desrochers & Laporte, OR Letters, 19. Feb. 91; 20. N = @SIZE( CITY); 21. MIN = @SUM( LINK: DIST * X); 22. @FOR( CITY( K): 23. ! It must be entered; 24. @SUM( CITY( I)| I #NE# K: X( I, K)) = 1; 25. ! It must be departed; 26. @SUM( CITY( J)| J #NE# K: X( K, J)) = 1; 27. ! Weak form of the subtour breaking constraints; 28. ! These are not very powerful for large problems; 29. @FOR( CITY( J)| J #GT# 1 #AND# J #NE# K: 30. U( J) >= U( K) + X ( K, J) - 31. ( N - 2) * ( 1 - X( K, J)) + 32. ( N - 3) * X( J, K))); 33. ! Make the X's 0/1; 34. @FOR( LINK: @BIN( X)); 35. ! For the first and last stop we know...; 36. @FOR( CITY( K)| K #GT# 1: 37. U( K) <= N - 1 - ( N - 2) * X( 1, K); 38. U( K) >= 1 + ( N - 2) * X( K, 1)); 39. END 40. [plain] view plaincopy 1. matlab代码： 2. 3. function main 4. clc,clear 5. global a 6. % a=zeros(6); 7. % a(1,2)=56;a(1,3)=35;a(1,4)=21;a(1,5)=51;a(1,6)=60; 8. % a(2,3)=21;a(2,4)=57;a(2,5)=78;a(2,6)=70; 9. % a(3,4)=36;a(3,5)=68;a(3,6)=68; a(4,5)=51;a(4,6)=61; 10. % a(5,6)=13; a=a+a'; 11. load cost 12. a=Muti_Cost;%边权矩阵 13. L=size(a,1); 14. c1=1:53; %初始圈 15. [circle,long]=modifycircle(c1,L); 16. c2=[1 53 2:52];%改变初始圈，该算法的最后一个顶点不动 17. [circle2,long2]=modifycircle(c2,L); 18. if long2 23. %******************************************* 24. %修改圈的子函数 25. %******************************************* 26. function [circle,long]=modifycircle(c1,L); 27. global a 28. flag=1; 29. while flag>0 30. flag=0; 31. for m=1:L-3 32. for n=m+2:L-1 33. if a(c1(m),c1(n))+a(c1(m+1),c1(n+1))<... 34. a(c1(m),c1(m+1))+a(c1(n),c1(n+1)) 35. flag=1; 36. c1(m+1:n)=c1(n:-1:m+1); 37. end 38. end 39. end 40. end 41. long=a(c1(1),c1(L)); 42. for i=1:L-1 43. long=long+a(c1(i),c1(i+1)); 44. end 45. circle=c1; 数学建模常用Matlab/Lingo/c代码总结系列——参数估计 分类： Matlab 数学建模2011-11-18 12:34 126人阅读 评论(0) 收藏 举报 Matlab中用fminsearch实现参数估计 文章的主要思想来源于Matlab|Simulink仿真世界的一篇类似的文章。我这里把这个思想引入到我们的体系来，并以一个新的例子讲解这一用法。 fminsearch用来求解多维无约束的非线性优化问题，它的基本形式是： [X,FVAL,EXITFLAG,OUTPUT] = FMINSEARCH(FUN,X0,OPTIONS). 大段的Matlab帮助文档我就不翻译解释了，有兴趣的朋友可以参见Matlab联机帮助，我这里只介绍他在参数估计中的作用。 在 参数估计中经常用到正态分布的参数估计。在matlab系统中有一个函数叫做normfit就直接可以完成这样的参数估计，返回均值mu和均方差 sigma的估计，但是这里有一个要求，就是它的输入信息必须是随机的数字序列。如得到1000个服从正态分布的随机数向量R，用命令[phat pci]=normfit(R)，就可以得到参数估计了。然而如果我我们得知某些已经处于pdf函数曲线上的点时，这时需要对函数进行拟合运算。 估计参数的原理是从已知的一序列数据中，对于给定的任何一组参数，计算用其估计数据得到的方差，然后利用fminsearch函数求当方差满足最小的时候的参数，这就是需要估计的参数。 来看一下下面的列子： [plain] view plaincopy 1. smu=10,ssig=25; 2. %假设原来均值方差分别为：10，25 3. R=randn(1000,1)*ssig+smu; 4. 5. %生成满足要求的1000个随机数 6. [y x]=myhist(R); 7. %生成统计信息，x，y分别表示分组中值序列和落入该组的统计数目 8. bar(x,y) 9. %绘制直方图 10. hold on 11. plot(x,y,'ro') 12. %绘制对应点 13. [pms mse]=normpdffit(x,y,8,20); 14. %根据得到的统计信息x,y对其进行参数估计，8，20分别代表均值和方差的初值 15. t=min(x):(max(x)-min(x))/200:max(x); 16. %定义绘图区间 17. ny=normpdf(t,smu,ssig); 18. %真实分布曲线数据 19. nyf=normpdf(t,pms(1),pms(2)); 20. %拟合分布曲线数据 21. plot(t,ny,'r-') 22. plot(t,nyf,'b-.') 23. legend('hist','hist value','ture pdf','fit pdf') 24. %绘制两条曲线作对比 上面例子中所用的几个函数定义如下： [plain] view plaincopy 1. function [h xout]=myhist(data,nbins) 2. %用于统计信息，生成和pdf函数值相同的hist统计方式。 3. if nargin==1 4. nbins=uint32(1+log(length(data))/log(2)); 5. end 6. nbins=double(nbins); 7. data=data(:); 8. [h xout]=hist(data,nbins); 9. ew=xout(2)-xout(1); 10. h=h./(ew*length(data)); 11. [plain] view plaincopy 1. function [pab mse]= normpdffit(x,y,a0,b0) 2. %正态分布pdf参数估计 3. p=[a0 b0]; 4. opt=optimset('fminsearch'); 5. opt.TolX=0.001; 6. opt.Display='off'; 7. [pab mse]=fminsearch(@normpdfse,p,opt,x,y); 8. %这里需要注意，opt参数已经传递给fminsearch，但是对于原计算方差的函数来说，还需要两个参数x,y,这两个参数就写在opt参数的后面，这样可以完成其他参数的传递。 9. %这里说下以前探索的时候的失败经验：用global把参数公有化，然后函数只传递变化的参数（需 要估计的参数），但是失败了。所以了解这种参数传递方法是非常有必要的。 10. function se= normpdfse(pab,X,Y) 11. %计算对于任何一组参数pab(1)，pab(2)，给出当前数据下的方差来。 12. se=var(Y-normpdf(X,pab(1),pab(2))); 运行结果如图所示： 从图中可以看出，随机数在这里变成了统计信息，统计信息反映到了绘制的点信息上，图中圆圈所示。真实的pdf为红色曲线，估计的曲线为蓝色虚线。从图中可知，估计的效果非常满意。 如果在函数中加上： [plain] view plaincopy 1. disp 'the result of normfit function:' 2. 3. [mu sg]=normfit(R) 4. 5. disp 'the result of fminsearch:' 6. 7. [pms mse]=normpdffit(x,y,8,20) 得到结果： [plain] view plaincopy 1. the result of normfit function: 2. 3. mu = 4. 5. 10.44306258428258 6. 7. sg= 8. 9. 25.61945417031251 10. 11. the result of fminsearch: 12. 13. pms = 14. 15. 10.30663244862284 25.32479396733891 16. 17. mse = 18. 19. 7.093014695522283e-008 与真实值相比，我们这里的拟合结果将比直接用normfit的结果更接近真实值。 可 以这么解释：normfit函数是内部通用的拟合函数，适合范围广，而没有任何先验信息加入，而对于我们的fminsearch函数来说，它需要一个先验信息，即参数的初值。我们在调用的时候用了初值8，20.这个先验信息对更进一步的拟合最后的结果有着相当重要的作用！因此，对于参数估计，先验信息还是相当重要的。 由于有这种技巧存在，了解了fminsearch函数之后，可以有很多的扩展应用。 数学建模常用Matlab/Lingo/c代码总结系列——层次分析法 分类： Matlab 数学建模2011-11-18 12:35 133人阅读 评论(0) 收藏 举报 [plain] view plaincopy 1. disp('请输入判断矩阵A(n阶)'); 2. A=input('A='); 3. [n,n]=size(A); 4. x=ones(n,100); 5. y=ones(n,100); 6. m=zeros(1,100); 7. m(1)=max(x(:,1)); 8. y(:,1)=x(:,1); 9. x(:,2)=A*y(:,1); 10. m(2)=max(x(:,2)); 11. y(:,2)=x(:,2)/m(2); 12. p=0.0001;i=2;k=abs(m(2)-m(1)); 13. while k>p 14. i=i+1; 15. x(:,i)=A*y(:,i-1); 16. m(i)=max(x(:,i)); 17. y(:,i)=x(:,i)/m(i); 18. k=abs(m(i)-m(i-1)); 19. end 20. a=sum(y(:,i)); 21. w=y(:,i)/a; 22. t=m(i); 23. disp(w);disp(t); 24. %以下是一致性检验 25. CI=(t-n)/(n-1);RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59]; 26. CR=CI/RI(n); 27. if CR<0.10 28. disp('此矩阵的一致性可以接受!'); 29. disp('CI=');disp(CI); 30. disp('CR=');disp(CR); 31. end 数学建模常用Matlab/Lingo/c代码总结系列——插值拟合 分类： Matlab 数学建模2011-11-18 12:39 221人阅读 评论(0) 收藏 举报 相关知识 在生产和科学实验中，自变量 与因变量 间的函数关系 有时不能写出解析表达式，而只能得到函数在若干点的函数值或导数值，或者表达式过于复杂而需要较大的计算量。当要求知道其它点的函数值时，需要估计函数值在该点的值。 为了完成这样的任务，需要构造一个比较简单的函数 ，使函数在观测点的值等于已知的值，或使函数在该点的导数值等于已知的值，寻找这样的函数 有很多方法。根据测量数据的类型有以下两类处理观测数据的方法。 （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] view plaincopy 1. year=1900:10:2010; 2. 3. 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] 4. 5. p1995=interp1(year,product,1995) 6. 7. x=1900:2010; 8. 9. y=interp1(year,product,x,'cubic'); 10. 11. 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年的劳动报酬表如下： 表：某企业工作人员的月平均工资（元） 20 30 服务年限 10 年份 1950 150.697 169.592 187.652 1960 179.323 195.072 250.287 1970 203.212 239.092 322.767 1980 226.505 273.706 426.730 1990 249.633 370.281 598.243 试计算1975年时，15年工龄的工作人员平均工资。 解：程序如下： [plain] view plaincopy 1. years=1950:10:1990; 2. 3. service=10:10:30; 4. 5. wage=[150.697 169.592 187.652 6. 7. 179.323 195.072 250.287 8. 9. 203.212 239.092 322.767 10. 11. 226.505 273.706 426.730 12. 13. 249.633 370.281 598.243]; 14. 15. mesh(service,years,wage) %绘原始数据图 16. 17. 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,15 16,22,28,35,27,20 18,21,26,32,28,25 20,25,30,33,32,20 求通过这些点的插值曲面。 解：程序为： [plain] view plaincopy 1. x=1:6; 2. 3. y=1:4; 4. 5. t=[12,10,11,11,13,15 6. 7. 16,22,28,35,27,20 8. 9. 18,21,26,32,28,25; 10. 11. 20,25,30,33,32,20] 12. 13. subplot(1,2,1) 14. 15. mesh(x,y,t) 16. 17. x1=1:0.1:6; 18. 19. y1=1:0.1:4; 20. 21. [x2,y2]=meshgrid(x1,y1); 22. 23. t1=interp2(x,y,t,x2,y2,'cubic'); 24. 25. subplot(1,2,2) 26. 27. mesh(x1,y1,t1); 结果如右图。 作业：已知某处山区地形选点测量坐标数据为： x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 55.5 6 海拔高度数据为： z=89 90 87 85 92 91 96 93 90 8782 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 9492 87 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得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则Y DELTA将至少包含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] view plaincopy 1. x=[0.5,1.0,1.5,2.0,2.5,3.0]; 2. 3. y=[1.75,2.45,3.81,4.80,7.00,8.60]; 4. 5. p=polyfit(x,y,2) 6. 7. x1=0.5:0.05:3.0; 8. 9. y1=polyval(p,x1); 10. 11. plot(x,y,'*r',x1,y1,'-b') 计算结果为： p =0.5614 0.8287 1.1560 此结果表示拟合函数为： ，用此函数拟合数据的效果如图所示。 例2：由离散数据 x y 0 .3 .1 .5 .2 1 .3 1.4 .4 1.6 .5 1.9 .6 .6 .7 .4 .8 .8 .9 1.5 1 2 拟合出多项式。 程序： x=0:.1:1; y=[.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.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.7459 10.9802 -0.0035 多项式为：16.7832x3-25.7459x2+10.9802x-0.0035 曲线拟合图形： 也可由函数给出数据。 例3：x=1:20,y=x+3*sin(x) 程序： [plain] view plaincopy 1. x=1:20; 2. 3. y=x+3*sin(x); 4. 5. p=polyfit(x,y,6) 6. 7. xi=linspace(1,20,100); 8. 9. z=polyval(p,xi); %¶àÏîʽÇóÖµº¯Êý 10. 11. plot(x,y,'o',xi,z,'k:',x,y,'b') 结果： p = 0.0000 -0.0021 0.0505 -0.5971 3.6472 -9.729511.3304 再用10阶多项式拟合 程序： [plain] view plaincopy 1. x=1:20; 2. 3. y=x+3*sin(x); 4. 5. p=polyfit(x,y,10) 6. 7. xi=linspace(1,20,100); 8. 9. z=polyval(p,xi); 10. 11. plot(x,y,'o',xi,z,'k:',x,y,'b') 结果：p = Columns 1 through 7 0.0000 -0.0000 0.0004 -0.01140.1814 -1.8065 11.2360 Columns 8 through 11 -42.0861 88.5907 -92.8155 40.267 可用不同阶的多项式来拟合数据，但也不是阶数越高拟合的越好。 数学建模常用Matlab/Lingo/c代码总结系列——非线性拟合 分类： Matlab 数学建模2011-11-18 12:41 130人阅读 评论(0) 收藏 举报 [plain] view plaincopy 1. function f=example1(c,tdata) 2. f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata)); [plain] view plaincopy 1. [plain] view plaincopy 1. x=1:1:12; 2. y=[0 3. 0 4. 0 5. 1 6. 0 7. 3 8. 10 9. 12 10. 8 11. 2 12. 1 13. 2 14. ]'; 15. c0=[2 8]; 16. for i=1:1000 17. c=lsqcurvefit(@zhengtai,c0,x,y); 18. c0=c; 19. end 20. y1=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c(2)^2)); 21. plot(x,y,'r-',x,y1); 22. legend('实验数据','拟合曲线') [plain] view plaincopy 1. x=[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]'; 2. y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4]'; 3. f=@(c,x)c(1)*(exp(-c(2)*x)-exp(-c(3)*x)); 4. c0=[114 0.1 2]'; 5. for i=1:50 6. opt=optimset('TolFun',1e-3); 7. [c R]=nlinfit(x,y,f,c0,opt) 8. c0=c; 9. hold on 10. plot(x,c(1)*(exp(-c(2)*x)-exp(-c(3)*x)),'g') 11. end t=[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];y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4];c0=[1 1 1];for i=1:50 c=lsqcurvefit(@example1,c0,t,y); c0=c;endy1=c(1)*(exp(-c(2)*t)-exp(-c(3)*t));plot(t,y,'+',t,y1);legend('实验数据','拟合曲线') 数学建模常用Matlab/Lingo/c代码总结系列——灰色预测 分类： Matlab 数学建模2011-11-18 12:47 205人阅读 评论(0) 收藏 举报 [plain] view plaincopy 1. 3. X=[136 143 165 152 165 181 204 272 319 491 571 605 665 640 628]; 4. x1(1)=X(1); 5. X1=[]; 6. for i=1:1:14 7. x1(i+1)=x1(i)+X(i+1); 8. X1=[X1,x1(i)]; 9. end 10. X1=[X1,X1(14)+X(15)] 11. for k=3:1:15 12. p(k)=X(k)/X1(k-1); 13. p1(k)=X1(k)/X1(k-1); 14. end 15. p,p1 16. clear k 17. Z=[]; 18. for k=2:1:15 19. z(k)=0.5*X1(k)+0.5*X1(k-1); 20. Z=[Z,z(k)]; 21. end 22. Z 23. B=[-Z',ones(14,1)] 24. 25. Y=[]; 26. clear i 27. for i=2:1:15 28. Y=[Y;X(i)]; 29. end 30. Y 31. A=inv(B'*B)*B'*Y 32. clear k 33. y1=[]; 34. for k=1:1:15 35. y(k)=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1); 36. y1=[y1;y(k)]; 37. end 38. y1 39. clear k 40. X2=[]; 41. for k=2:1:15 42. x2(k)=y1(k)-y1(k-1); 43. X2=[X2;x2(k)]; 44. end 45. X2=[y1(1);X2] 46. e=X'-X2 47. m=abs(e)./X' 48. s=e'*e 49. n=sum(m)/13 50. clear k 51. syms k 52. y=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1) 53. Y1=[]; 54. for j=16:1:21 55. y11=subs(y,k,j)-subs(y,k,j-1); 56. Y1=[Y1;y11]; 57. end 58. Y1 %程序中的变量定义:alpha是包含α、μ值的矩阵;%ago是预测后累加值矩阵;var是预测值矩阵;%error是残差矩阵; c是后验差比值function basicgrey(x,m) %定义函数basicgray(x)if nargin==1 %m为想预测数据的个数，默认为1 m=1;endclc; %清屏,以使计算结果独立显示if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换 x=x';endn=length(x); %取输入数据的样本量x1(:,1)=cumsum(x); %计算累加值,并将值赋与矩阵be for i=2:n %对原始数列平行移位 Y(i-1,:)=x(i,:);endfor i=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; %计算参数α、μ矩阵for i=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,:);for i=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) %评价函数 显示这个模型是否合格 [plain] view plaincopy 1. 3. if rela>0.6 4. '根据经验 关联度检验结果为满意（关联度只是参考 主要看后验差的结果）' 5. else 6. '根据经验 关联度检验结果为不满意（关联度只是参考 主要看后验差的结果）' 7. end 8. if P>0.95&c<0.5 9. '后验差结果显示 这个模型评价为“优”' 10. else if P>0.8&c<0.5 11. '后验差结果显示 这个模型评价为“合格”' 12. else if P>0.7&c<0.65 13. '后验差结果显示 这个模型评价为“勉强合格”' 14. else 15. '后验差结果显示 这个模型评价为“不合格”' 16. end 17. end 18. end [plain] view plaincopy 1. function [P,c,error]=lcheck(x,var) 2. %进行后验差检验 3. n=length(x); 4. for i=1:n 5. error(i,:)=abs(var(i,:)-x(i,:)); %计算绝对残差 6. end 7. c=std(abs(error))/std(x); %调用统计工具箱的标准差函数计算后验差的比 值c 8. s0=0.6745*std(x); 9. ek=abs(error-mean(error)); 10. pk=0; 11. for i=1:n 12. if ek(i,:) [plain] view plaincopy 1. %附带的质料里有一部分讲了关联度 2. function [rela]=relations(x) 3. %以x(1,:)的参考序列求关联度 4. [m,n]=size(x); 5. for i=1:m 6. for j=n:-1:2 7. x(i,j)=x(i,j)/x(i,1); 8. end 9. end 10. for i=2:m 11. x(i,:)=abs(x(i,:)-x(1,:)); %求序列差 12. end 13. c=x(2:m,:); 14. Max=max(max(c)); %求两极差 15. Min=min(min(c)); 16. p=0.5; %p称为分辨率,0 19. r(i,j)=(Min+p*Max)/(c(i,j)+p*Max); %计算关联系数 20. end 21. end 22. for i=1:m-1 23. rela(i)=sum(r(i,:))/n; %求关联度 24. end 25. 数学建模常用Matlab/Lingo/c代码总结系列——Matlab图形绘制函数汇总 分类： 数学建模 Matlab2011-11-20 00:16 438人阅读 评论(0) 收藏 举报 基本绘图和图形 boxerrorbarholdline坐标轴边界 沿曲线绘制误差条 在图形窗口中保留当前图形 创建线条对象 LineSpec (Line Specification)loglogplotplot3plotyypolarsemilogx, semilogysubplot线条规格字符串语法 对数-对数刻度图 二维线条图 三维线条图 y轴分居左右两侧的线条图 极坐标图 半对数坐标图 在窗口的平铺位置创建坐标轴 绘图工具 figurepalettepanplotbrowserploteditplottoolspropertyeditorrotate3dshowplottoolzoom显示或隐藏图形窗口的调色板 交互式移动图像以多方向浏览 显示或隐藏窗口的图形浏览器 交互式编辑和标注图形 显示或隐藏图形工具 显示或隐藏属性编辑器 使用鼠标旋转三维视图 显示或隐藏窗口的图形工具 图形的放大、缩小或按比例缩放 标注图形 annotationclabeldatacursormodedatetickgtextlegendrectangletexlabeltitlexlabel, ylabel, zlabel创建注释对象 等高线高程标签 使能或禁止交互式数据光标模式 日期格式的刻度标签 在二维视图中利用鼠标放置文本 线条和补片对象的图例 创建二维矩形对象 产生Tex格式的字符串 为当前坐标轴添加标题 X，Y和Z轴的标签 专业绘图(Area、条形图、圆饼图) areabar, barh填充区域的二维图形 绘制条形图(垂直和水平) bar3, bar3hparetopiepie3绘制三维条形图 帕累托图表 饼图 三维饼图 专业绘图(等高线图) contourcontour3contourccontourfezcontourezcontourf等高线图矩阵 三维等高线图 低层次的等高线图计算 填充二维等高线图 易于使用的轮廓绘图仪 易于使用的填充轮廓绘图仪 专业绘图(方向和速度图) cometcomet3compassfeatherquiverquiver3二维彗星图 三维彗星图 绘制从原点射出的箭头 绘制速度矢量图 抖动或速度图 三维抖动或速度图 专业绘图(离散数据图) stairsstemstem3阶梯图 情节离散序列数据 情节 3 - D离散序列数据 专业绘图(函数绘图) ezcontourezcontourfezmeshezmeshcezplotezplot3ezpolarezsurf易用的等高线绘图仪 易用的填充的等高线绘图仪 易用的三维网格绘图仪 易用的组合式的网格 /等高线绘图仪 易用的函数绘图仪 易用的三维参数化曲线绘图仪 易用的极坐标绘图仪 易用的三维彩色绘图仪 ezsurfcfplot易用的曲面/等高线绘图仪 在指定的坐标范围内绘制图形 专业绘图(直方图) histhistcrose直方图阴谋 直方图计数 角直方图阴谋 专业绘图(多边形和曲面) cylinderdelaunaydelaunay3delaunayndsearchellipsoidfillfill3inpolygonpcolorpolyarearectintribbonslicespherewaterfall生成缸 Delaunay三角网 三维 Delaunay镶嵌 N维 Delaunay镶嵌 搜索Delaunay三角网的最近点 生成椭球 填充二维多边形 填充的3 - D多边形 多边形区域内点 伪(棋盘)阴谋 多边形面积 矩形的交集区 丝带阴谋 容积片情节 生成领域 瀑布图 专业绘图(散点/气泡图) plotmatrixscatterscatter3散点图矩阵 散点图 三维散点图 专业绘图(动画) frame2imgetframeim2frame返回图像数据与电影相关的框架 电影帧捕获 帧图像转换到电影 movienoanimate播放录制的电影帧 所有对象的变化 EraseMode正常 位图图像 frame2imim2frameim2javaimageimagescimfinfoimformatsimreadimwriteind2rgb返回图像数据与电影相关的框架 帧图像转换到电影 图像转换到Java形象 显示图像对象 规模的数据和显示图像对象 信息图形文件 管理图像文件格式的注册表 阅读图像从图形文件 写入图像图形文件 索引图像转换到RGB图像 打印 hgexportorientprint, printoptprintdlgprintpreviewsaveas出口数字 印刷纸张方向 打印数字或保存到文件和打印机的默认配置 打印对话框 打印预览图 保存数字或Simulink框图使用指定的格式 句柄图形(查找和辨识图形对象) allchildancestorcopyobjdeletefindallfindfigsfindobjgcagcbfgcbo查找指定对象的所有儿童 祖先的图形对象 复制的图形对象和他们的后代 删除文件或图形对象 查找所有的图形对象 查找离屏数字可见 找到具有特殊性能的图形对象 当前轴手柄 处理数字包含对象，其回调是执行 处理的对象，其回调是执行 gcogetishandlepropeditset处理当前对象 查询处理图形对象的属性 确定是否输入是有效的句柄图形处理 打开属性编辑器 处理图形对象的属性设置 句柄图形(对象创建函数) axesfigurehggrouphgtransformimagelightlinepatchrectangleroot objectsurfacetextuicontextmenu创建轴图形对象 创建数字图形对象 创建 hggroup对象 创建 hgtransform图形对象 显示图像对象 根据创建对象 创建线对象 创建一个或多个填充多边形 创建二维矩形对象 根 创建面对象 创建文本对象在当前轴 创建上下文菜单 句柄图形(绘图对象) Annotation Annotation Annotation Annotation Annotation Annotation Annotation Arrow PropertiesDoublearrow PropertiesEllipse PropertiesLine PropertiesRectangle PropertiesTextarrow PropertiesTextbox Properties注释定义箭头性能 定义注释 doublearrow性能 椭圆定义批注的属性 注释行属性定义 注释定义矩形属性 定义注释 textarrow性能 定义注释文本框属性 属性定义 areaseries 属性定义 barseries 属性定义 contourgroup 属性定义 errorbarseries 定义图像属性 Areaseries PropertiesBarseries PropertiesContourgroup PropertiesErrorbarseries PropertiesImage PropertiesLineseries PropertiesQuivergroup PropertiesScattergroup PropertiesStairseries PropertiesStemseries PropertiesSurfaceplot Properties属性定义 lineseries 属性定义 quivergroup 属性定义 scattergroup 属性定义 stairseries 属性定义 stemseries 属性定义 surfaceplot 句柄图形(图形窗口) clfcloseclosereqdrawnowgcfhgloadhgsavenewplotopenglrefreshsaveas清除目前的数字窗口 删除指定的数字 默认关闭请求函数图 刷新事件队列和更新的数字窗口 目前的数字处理 负载处理图形对象的层次结构从文件 保存对象的层次结构来处理图形文件 确定要绘制图形对象 控制OpenGL渲染 目前的数字重绘 保存数字或Simulink框图使用指定的格式 句柄图形(坐标轴操作) axisboxclagcagridisholdmakehgtform轴缩放和外观 轴边界 清除当前轴 当前轴手柄 电网线 2 - D和3 - D图 当前保持状态 创建 4 × 4变换矩阵 句柄图形(对象属性操作) getlinkaxeslinkproprefreshdataset查询处理图形对象的属性 同步限制指定的2 - D轴 保持相应的属性值相同 刷新数据图表数据源时指定 处理图形对象的属性设置 数学建模常用Matlab/Lingo/c代码总结系列——hamilton回路 分类： 数学建模 Matlab2011-11-20 00:20 233人阅读 评论(0) 收藏 举报 提供一种求解最优哈密尔顿的算法---三边交换调整法，要求在运行jiaohuan3（三交换法）之前，给定邻接矩阵C和节点个数N，结果路径存放于R中。 bianquan.m文件给出了一个参数实例，可在命令窗口中输入bianquan,得到邻接矩阵C和节点个数N以及一个任意给出的路径R,，回车后再输入jiaohuan3,得到了最优解。 由于没有经过大量的实验，又是近似算法，对于网络比较复杂的情况，可以尝试多运行几次jiaohuan3,看是否能到进一步的优化结果。 [plain] view plaincopy 1. %%%%%%bianquan.m%%%%%%% 2. N=13; 3. for i=1:N 4. for j=1:N 5. C(i,j)=inf; 6. end 7. end 8. for i=1:N 9. C(i,i)=0; 10. end 11. C(1,2)=6.0;C(1,13)=12.9; 12. C(2,3)=5.9;C(2,4)=10.3; 13. C(3,4)=12.2;C(3,5)=17.6; 14. C(4,13)=8.8;C(4,7)=7.4; 15. C(4,5)=11.5; 16. C(5,2)=17.6;C(5,6)=8.2; 17. C(6,9)=14.9;C(6,7)=20.3; 18. C(7,9)=19.0;C(7,8)=7.3; 19. C(8,9)=8.1;C(8,13)=9.2; 20. C(9,10)=10.3; 21. C(10,11)=7.7; 22. C(11,12)=7.2; 23. C(12,13)=7.9; 24. for i=1:N 25. for j=1:N 26. if C(i,j) < inf 27. C(j,i)=C(i,j); 28. end 29. end 30. end 31. for i=1:N 32. C(i,i)=0; 33. end 34. 35. R=[4 7 6 5 3 2 1 13 12 11 10 9 8]; [plain] view plaincopy 1. 3. for I=1:(N-2) 4. for J=(I+1):(N-1) 5. for K=(J+1):N 6. n=n+1; 7. Z(n,:)=[I J K]; 8. end 9. end 10. end 11. R=1:N 12. for m=1:(N*(N-1)*(N-2)/6) 13. I=Z(m,1);J=Z(m,2);K=Z(m,3); 14. r=R; 15. if J-I~=1&K-J~=1&K-I~=N-1 16. for q=1:(J-I) 17. r(I+q)=R(J+1-q); 18. end 19. for q=1:(K-J) 20. r(J+q)=R(K+1-q); 21. end 22. end 23. if J-I==1&K-J==1 24. r(K)=R(J);r(J)=R(K); 25. end 26. if J-I==1&K-J~=1&K-I~=N-1 27. for q=1:(K-J) 28. r(I+q)=R(I+1+q); 29. end 30. r(K)=R(J); 31. end 32. if K-J==1&J-I~=1&K~=N 33. for q=1:(J-I) 34. r(I+1+q)=R(I+q); 35. end 36. r(I+1)=R(K); 37. end 38. if I==1&J==2&K==N 39. for q=1:(N-2) 40. r(1+q)=R(2+q); 41. end 42. r(N)=R(2); 43. end 44. if I==1&J==(N-1)&K==N 45. for q=1:(N-2) 46. r(q)=R(1+q); 47. end 48. r(N-1)=R(1); 49. end 50. if J-I~=1&K-I==N-1 51. for q=1:(J-1) 52. r(q)=R(1+q); 53. end 54. r(J)=R(1); 55. end 56. if J==(N-1)&K==N&J-I~=1 57. r(J+1)=R(N); 58. for q=1:(N-J-1) 59. r(J+1+q)=R(J+q); 60. end 61. end 62. if cost_sum(r,C,N) 66. fprintf('总长为%f\\n',cost_sum(R,C,N)) %%%%%%cost_sum.m%%%%%%%%function y=cost_sum(x,C,N)y=0;for i=1:(N-1) y=y+C(x(i),x(i+1));endy=y+C(x(N),x(1)); 大白话解析模拟退火算法 分类： 数学建模2011-11-20 00:30 126人阅读 评论(0) 收藏 举报 一. 爬山算法 ( Hill Climbing ) 介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。 爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。 图1 二. 模拟退火(SA,Simulated Annealing)思想 爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局 的最优解。以图1为例，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。 模拟退火算法描述： 若J( Y(i+1) )>= J( Y(i) ) (即移动后得到更优解)，则总是接受该移动 若J( Y(i+1) )< J( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定） 这里的―一定的概率‖的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。 根据热力学的原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为： P(dE) = exp( dE/(kT) ) 其中k是一个常数，exp表示自然指数，且dE<0。这条公式说白了就是：温度越高，出现一次能量差为dE的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于dE总是小于0（否则就不叫退火了），因此dE/kT < 0 ，所以P(dE)的函数取值范围是(0,1) 。 随着温度T的降低，P(dE)会逐渐降低。 我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。 关于爬山算法与模拟退火，有一个有趣的比喻： 爬山算法：兔子朝着比现在高的地方跳去。它找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法，它不能保证局部最优值就是全局最优值。 模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向高处，也可能踏入 平地。但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火。 下面给出模拟退火的伪代码表示。 三. 模拟退火算法伪代码 代码 [plain] view plaincopy 1. /* 2. * J(y)：在状态y时的评价函数值 3. * Y(i)：表示当前状态 4. * Y(i+1)：表示新的状态 5. * r： 用于控制降温的快慢 6. * T： 系统的温度，系统初始应该要处于一个高温的状态 7. * T_min ：温度的下限，若温度T达到T_min，则停止搜索 8. */ 9. while( T > T_min ) 10. { 11. dE = J( Y(i+1) ) - J( Y(i) ) ; 12. 13. if ( dE >=0 ) //表达移动后得到更优解，则总是接受移动 14. Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动 15. else 16. { 17. // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也 18. if ( exp( dE/T ) > random( 0 , 1 ) ) 19. Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动 20. } 21. T = r * T ; //降温退火 ，0 23. * 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。若r过小，则搜索 的过程会很快，但最终可能会达到一个局部最优值 24. */ 25. i ++ ; 26. } 四. 使用模拟退火算法解决旅行商问题 旅行商问题 ( TSP , Traveling Salesman Problem ) ：有N个城市，要求从其中某个问题出发，唯一遍历所有城市，再回到出发的城市，求最短的路线。 旅行商问题属于所谓的NP完全问题，精确的解决TSP只能通过穷举所有的路径组合，其时间复杂度是O(N!) 。 使用模拟退火算法可以比较快的求出TSP的一条近似最优路径。（使用遗传算法也是可以的，我将在下一篇文章中介绍）模拟退火解决TSP的思路： 1. 产生一条新的遍历路径P(i+1)，计算路径P(i+1)的长度L( P(i+1) ) 2. 若L(P(i+1)) < L(P(i))，则接受P(i+1)为新的路径，否则以模拟退火的那个概率接受P(i+1) ，然后降温 3. 重复步骤1，2直到满足退出条件 产生新的遍历路径的方法有很多，下面列举其中3种： 1. 随机选择2个节点，交换路径中的这2个节点的顺序。 2. 随机选择2个节点，将路径中这2个节点间的节点顺序逆转。 3. 随机选择3个节点m，n，k，然后将节点m与n间的节点移位到节点k后面。 五. 算法评价 模拟退火算法是一种随机算法，并不一定能找到全局的最优解，可以比较快的找到问题的近似最优解。 如果参数设置得当，模拟退火算法搜索效率比穷举法要高。 Matlab用xcorr求自相关函数出现的一点问题 分类： Matlab2011-11-30 23:51 1251人阅读 评论(0) 收藏 举报 用Matlab做随机信号分析实验仿真的时候出现一个问题：使用xcorr计算信号的自相关函数并求Rx(0)=max(Rx)，然后通过E(x.^2)求均方值 σx^2，理论二者都是表示平均功率，应该相等，但是实验过程中得出高斯白噪声的均方值0.1495,而Rx(0)可达30以上。如下图： 经过阅读matlab的help文档，发现原因是xcorr缺省计算原始的没有标准化的相关性（求和而非求均值），需要加上‘unbiased‘选项计算自相关函数的无偏估计。 ========================================================================================================================== 以下是整理的xcorr函数原始help文档和部分翻译，供参考： xcorr Cross-correlation Syntax c=xcorr(x,y) c=xcorr(x) c=xcorr(x,y,'option') c=xcorr(x,'option') c=xcorr(x,y,maxlags) c=xcorr(x,maxlags) c=xcorr(x,y,maxlags,'option') c=xcorr(x,maxlags,'option') [c,lags]=xcorr(...) Description xcorr estimates the cross-correlation sequence of a random process.Autocorrelation is handled as a special case. Xcorr估计一个随机过程的互相关序列，自相关是当做一个特例处理的。 The true cross-correlation sequence is 真实的互相关序列是 where xn and yn are jointly stationaryrandom processes, −∞ c=xcorr(x,y) returns the cross-correlation sequence in alength 2*N-1 vector, where x and y are lengthN vectors (N>1).If x and y are not the same length, the shorter vector is zero-padded to thelength of the longer vector. By default, xcorr computes raw correlations with no normalization. 缺省情况下xcorr计算原始的没有标准化的相关性（这里就是开始错误的原因所在，请注意以下表达式，貌似下面计算的是确定性信号的自相关函数，不同于随机信号，求指正） The output vector c has elements given by c(m) = Rxy(m−N),m=1, ..., 2N−1. 输出向量c的元素通过c(m) = Rxy(m−N),m=1, ..., 2N−1.给出 In general, the correlation function requires normalizationto produce an accurate estimate (see below). 一般情况下，相关性函数需要标准化以得出一个精确的估计（见下面） c = xcorr(x) is the autocorrelationsequence for the vector x. If x is anN-by-P matrix, c is amatrix with 2N-1 rows whose Pcolumns contain thecross-correlation sequences for all combinations of the columns of x. For moreinformation on matrix processing with xcorr, seeMultiple Channels. c = xcorr(x)得出向量x的自相关序列，如果x是N*P的矩阵，c是一个2N-1行的矩阵，它的P^2列包含x的所有列组合的互相关序列。xcorr矩阵处理的更多信息，请参阅多Multiple Channels. xcorr produces correlations identically equal to 1.0 at zero lag only whenyou perform an autocorrelation and only when you set the 'coeff' option. Forexample, 只有当你求一个自相关和只有当你设置了'coeff―选项时，xcorr在零滞后产生恒等于1.0的相关性。例如， [plain] view plaincopy 2 1. x=0:0.01:10; 2. 3. X = sin(x); 4. 5. [r,lags]=xcorr(X,'coeff'); 6. 7. max(r) c = xcorr(x,y,'option') specifies anormalization option for the cross-correlation, where'option' is c = xcorr(x,y,'option')为互相关指定一个标准化的选项，当选项是 'biased': Biased estimate of the cross-correlation function 'biased':互相关函数的有偏估计（这里涉及有偏和无偏估计，有兴趣可以参考相关文献，此处不仔细讲解） 'unbiased': Unbiased estimate of the cross-correlation function 'unbiased':互相关函数的无偏估计 'coeff': Normalizes the sequence so the autocorrelations at zero lag are identically 1.0. 'coeff':标准化序列使自相关在零滞后等于1.0 'none', to use the raw, unscaled cross-correlations (default) 'none'，使用原始的未缩放的互相关性（缺省）。 See [1] for moreinformation on the properties of biased and unbiased correlation estimates. c = xcorr(x,'option') specifies one of theabove normalization options for the autocorrelation. c = xcorr(x,y,maxlags) returns thecross-correlation sequence over the lag range [-maxlags:maxlags]. Output c haslength 2*maxlags+1. c = xcorr(x,maxlags) returns theautocorrelation sequence over the lag range [-maxlags:maxlags]. Output c haslength 2*maxlags+1. If x is anN-by-P matrix, c is a matrix with2*maxlags+1 rows whose P2 columns contain the autocorrelationsequences for all combinations of the columns of x. c = xcorr(x,y,maxlags,'option')specifies both a maximum number of lags and a scaling option for thecross-correlation. c = xcorr(x,maxlags,'option') specifiesboth a maximum number of lags and a scaling option for the autocorrelation. [c,lags] = xcorr(...) returns a vector of thelag indices at which c was estimated, with the range [-maxlags:maxlags]. Whenmaxlags is not specified, the range of lags is [-N+1:N-1]. In all cases, the cross-correlation or autocorrelation computed by xcorrhas the zeroth lag in the middle of the sequence, at element or row maxlags+1(element or row N if maxlags is not specified). Examples The second output, lags, is useful for plotting the cross-correlation orautocorrelation. For example, the estimated autocorrelation of zero-meanGaussian white noisecww(m) can be displayed for -10 ≤m ≤ 10 using: [plain] view plaincopy 1. ww = randn(1000,1); 2. 3. [c_ww,lags] = xcorr(ww,10,'coeff'); 4. 5. stem(lags,c_ww) Swapping the x and y input arguments reverses (and conjugates) the outputcorrelation sequence. For row vectors, the resulting sequences are reversedleft to right; for column vectors, up and down. The following example illustratesthis property (mat2str is used for a compact display of complex numbers): [plain] view plaincopy 1. x = [1,2i,3]; y = [4,5,6]; 2. 3. [c1,lags] = xcorr(x,y); 4. 5. c1 = mat2str(c1,2), lags 6. 7. c2 = conj(fliplr(xcorr(y,x))); 8. 9. c2 = mat2str(c2,2) For the case where input argument x is a matrix, the output columns arearranged so that extracting a row and rearranging it into a square arrayproduces the cross-correlation matrix corresponding to the lag of the chosenrow. For example, the cross-correlation at zero lag can be retrieved by: [plain] view plaincopy 1. randn('state',0) 2. 3. X = randn(2,2); 4. 5. [M,P] = size(X); 6. 7. c = xcorr(X); 8. 9. c0 = zeros(P); c0(:) = c(M,:) % Extract zero-lag row You can calculate the matrix of correlation coefficients that the MATLABfunctioncorrcoef generates by substituting: [plain] view plaincopy 1. c = xcov(X,'coef') in the last example. The function xcov subtractsthe mean and then calls xcorr. Use fftshift to move the second half of the sequence starting at the zeroth lag to thefront of the sequence. fftshift swaps the first and second halves of a sequence. Algorithm For more information on estimating covariance and correlation functions,see[1]. References [1] Orfanidis, S.J., Optimum Signal Processing. AnIntroduction. 2nd Edition, Prentice-Hall, Englewood Cliffs, NJ, 1996. 关于数模中编程的一点愚见 分类： 数学建模2011-12-10 20:15 2622人阅读 评论(21) 收藏 举报 2011数模国赛已经过去整整三个月，作为负责编程的队员，在此发表一点愚见，也作为年末一点总结，请各位选择性吸收，欢迎拍砖。 编程语言，这个估计是大家最关心的。数模中编程语言首选Matlab，世界公认加默认。当然c语言和其他高级语言也可以使用，毕竟过于通用，在使用上效率不足Matlab。这个效率不是运行的效率，而是解题的效率。但是如果你是c语言或其他语言的顶级大牛，就当我前面说的都是废话，甚至直接忽略我这篇文章，因为你完全有能力搞定所有这些问题，我所说的不过是给普通大众的一点捷径而已，算得上―旁门左道‖。可惜，90%以上的童鞋不是。 Matlab由于其强大的科学计算功能，以及封装的各种toolbox，成为建模编程的得力助手。c在这点上略逊一筹，很多都得自己动手写，尽管c++提供了各种模板库之类的东西，运用起来也不是像Matlab这样的轻松。至于编程语言方面的比较，我也不想说太多，因为语言各有利弊，Matlab卖的那么贵也是有它的理由的。 Matlab建模优势很大。我仅仅抛砖引玉说几点。首先是各种toolbox和function，它们绝大部分是Mathworks公司的顶级工程师的研究成果，当然有些也是成熟算法的Matlab实现，你也可以自己去Google上寻找业余Matlab爱好者写的toolbox，里面不乏非常出名的算法大牛。Matlab可以轻易的进行矩阵计算，二维三维图形绘制，概率统计，信号（图像也是信号）处理（学通信的童鞋都知道，这个在大三有随机信号分析和数字信号处理等课程要用到），Simulink系统仿真......甚至嵌入式方面也有涉足（记得最初是从奔哥那得知的），还能轻易地其他语言混合编程（本人曾经试过用Matlab里的deploytool将自己写的function编译成托管dll供C#调用，实现简单的C#图像处理应用。具体操作方法不在此介绍，有兴趣者可以参考Matlab高级编程或者Google，Mathworks官网也有简单的视频教程）。 特别要提的是第三方Matlab toolbox。之所以要提得原因就是，数模的时候，时间紧迫，要是碰到编程很复杂的题目，可能三天都无法出来结果，这个是最可怕的。我记得集训期间，做相机标定的那个题，机器视觉这种研究热门，自然少不了Matlab的身影，于是找到一个Camera Calibration Toolbox for Matlab（http://www.vision.caltech.edu/bouguetj/calib_doc/），最终由于出题思路与普通方法有差距而没用到这个toolbox。但是今年国赛的A题，寻找污染源的问题，需要用到Kriging插值，由于普通教程讲解和公式推导晦涩难懂，编程难以下手，当时打算使用ArcGIS软件直接进行插值后图形绘制，然而事情终不是如你所愿，一个软件3G左右，实验室流量有限，网速有限，要是下下来整个实验室的队伍全完了，而且也要等到猴年马月，估计国赛都完了。因此，遇到这种情况得保持镇定，其实只要稍微花上5分钟Google一下，你就会找到一个名叫DACE-A Matlab Kriging Toolbox（ http://www2.imm.dtu.dk/~hbn/dace/）的东西，还有详细的使用例子，当然，前提是你能读懂简单的英文（这个后面再说）。也正是这个toolbox，省掉了一大堆代码和时间，解题效率必然呈指数式增长。 其实平时训练Matlab的过程中，我不建议太过于依赖toolbox，那个是纯粹的比赛技巧。数模国赛组委会和Mathworks公司近两年联合推出了一个Matlab创新奖，我想获得这个奖肯定不会是靠投机取巧所得。 要是你不喜欢Matlab，感觉它没什么意思，那你错了。在Matlab的命令窗口输入demo，在打开的Help browser里你会发现如 Minesweeper （扫雷） World Traveler 3-D Globe （3D地球仪） 之类的GUI小游戏和美妙的2D、3D图形，全部是用Matlab实现的，绝对会让你眼前一亮。当然网络上也有很多发烧友的小程序，也非常有趣（比如萝卜驿站：http://luobo.ycool.com/）。关于如何学习Matlab以及一些Matlab的操作技巧，我也不多说，网上教程很多，这方面的课本和资料也十分丰富。只要培养起兴趣，平时多动手谢谢代码，都会成为高手的。 说了那么多Matlab，提谈谈刚才提到的英语。数学建模中英语也是必过的一关，美赛尤其重要（这个我也没参加，先不说）。因为编程的队员负责的不能只是编程这一部分，也必须懂得资料和数据的查找，海量数据处理，建模相关数学知识，以及论文格式排版和内容修改，会使用绘图等工具软件and so on。你不可能提前得知题目的类型，尽管类型也只有那么几十种，但是资料在哪里找却不能预测。英文的功底决定你可以上国外的网站（也得懂点翻墙什么的技巧）找相关数据，可以看懂相关英文文献，可以使用相关英文软件（毕竟不是所有软件都是国产和汉化后的，国产的一般不给力，但是有能力也可以自己提前汉化一下通用的英文软件，不过这个很费时间，考虑到个人因素，大家看着办，想当初自己费劲心血汉化完Lingo后，发现网上已经有了汉化版，瀑布汗当时。但是后来发现那个汉化版很多还是不太准确的，用自己的实在。这里有侵权之嫌，友情提示最好不要作为商业用途，后果自负）。另外就是Matlab的Help之类的了，有些几乎找不到找到中文解释的，必须得咬咬牙看英文了。 另外Lingo也是一种很好的建模语言，在求解优化问题等方面具有得天独厚的优势，由于本人基本上没有实践过，也毫无经验可言，不过推荐大家学习学习。 前面说的是编程语言，其实数模过程中语言都是次要的。重要的是结果，结果重于过程，有悖于常规观念。但是这是事实。除非碰到类似于公交车调度一样的题目，写出来的代码评委一般不太重视的，代码最终也只是当做附录处理，毕竟数模不是ACM，不是纯粹的程序设计。数模看得是解决的方法和你得出结果，至于你是怎么得出正确结果的，那都无关紧要。编程过程中也不需要写出速率最优化的代码，但是前提是保证能短时间内出结果。这些对于一个专业程序员来说也许是很不好的习惯，但是记住，这是数模。你可以写一千行代码得出结果，也可以写10行得出一样的结果，可以运行1分钟出结果，也可以运行1s出结果，只要结果正确了，其他都不重要了。但是记住，这是比赛，不是训练，训练要有成果，必须练好基本功。 说到结果，数模中可以―不择手段‖。这也算是编程队员的一种必备的能力。比如拟合一个特别复杂的非线性方程，Matlab处理的时候有个―缺点‖是必须选择合适的初始值，否则选择不当，迭代N次后拟合效果可能千差万别。而正确选取初始值也是十分麻烦的事，记得以前做本底趋势线拟合的时候找到一篇专门研究如何估计初始值的论文，不过后来放弃了。不仅是Matlab，专业拟合软件如Origin等也要手工选取初始值。但是使用过1stOpt的人，却可以不费吹灰之力得到很好的结果。1stOpt虽然是国产软件，在拟合和优化方面却力压群雄，提供神经网络、遗传算法、蚁群算法、模拟退火等专业算法的选择，甚至包括一种自创的高效算法，在无需手动设置初始值的情况下短时间内拟合出方程并绘图和预测。题目中拟合只占了很小一部分，若是花太多时间在代码上，无疑会在比赛中处于下风。 这样，学会各种专业软件，具备短时间内学会小众软件的能力也是编程队员的优势。我曾经做过一个小的整理，并对各种主流和很有用的小众软件进行了一个归类，如下（不完整，可以按需要去收集，数学中国上面有）： 科学计算：Matlab、Maple（符号计算） 、Mathematica、Excel（绘图、统计） 概率统计：SPSS、SAS、Eviews、Origion（拟合） 系统动力学仿真：Vensim 优化：Lingo/Lindo、1stOpt 、WinQSB（规划） 绘图 ：SigmaPlot（专业的科学绘图软件） 、ScienceWord、SmartDraw（很强大的绘图工具）、几何画板 、autoCAD 其他小众软件：图论、AHP（层次分析法）… 其他功能软件：tortoiseSVN（版本控制）、dexpot（虚拟桌面）、酷盘（局域网共享协作） 其他：ansys(有限元分析)、comsol、FLAC2D/FLAC3D... 有了这些杀手锏，基本上效率可以提高一倍了。其中要提一下的是Excel，很多人觉得它功能不足，其实Excel在数据处理和绘图方面毫不逊色与专业软件，特别是最新的Excel2010，看过有人用它来画动漫人物的视频，至今膜拜不止，不过数模中不推荐用它绘图，太花哨了。 至于注重自己实现代码，还是使用工具走捷径，很像哲学中的唯物主义和唯心主义，要是你只是功利性地想获奖，那就做个唯物主义者吧。 刚才归类软件的时候，也将版本控制等功能软件放进去了。这个是有目的的。版本控制很有用，这个对很多程序员们都再熟悉不过了。要是没听过，就不厌其烦听我说我吧。编程的队员可能不会一次就完成代码，经常不断修改更新代码，但是要是修改过之后把以前版本删除了，后来又要用到怎么办，手动备份是见吃力不讨好的事。这时候版本控制软件就起作用了，配置好服务器和版本库，接下来的工作就只是简单的更新和提交了，系统自动打上时间戳，以便于后续的版本修改对比。数模中论文丢失莫过于最恐怖的事，这个同样交给版本控制，不能说万无一失（硬盘损坏也可能发生，做好备份很关键），也能减轻一点工作吧。 另外虚拟桌面有时候也有必要，这个得看个人喜好。要是你用的是Linux系统，那就完全可以略过这里。绝大部分童鞋还是用的Windows系统，毕竟Word写论文还是挺方便的（专业排版还是去用Latex吧，不过一般word足够，尽管word经常会出现很多奇怪的问题，Latex听说要写代码的，编程大师弄的东西就是不一样）。虚拟桌面也就是虚拟出几个桌面出来（很废话），以便管理混乱的窗口。数模的时候开个word，开个excel，开个Matlab，开个浏览器、记事本、图片查看器、资源管理器......这都是家常便饭，弄个虚拟桌面，占点资源换来方便也有价值。要是你有多台显示器，接在一台电脑上也不错，就不用什么虚拟桌面了。 当然还有那个局域网共享。数模是三个人的战斗，但是往往要使用多台电脑协作，资源共享是件麻烦事。比如写代码的时候，开两台配置好环境的计算机，以备不测，必能大大减轻单机的负荷。U盘太慢，也容易损坏文件，移动硬盘，估计不一定有，有也不方便。搭好局域网，做好共享很必要。QQ局域网传文件也可快，可惜无法多人同时协作和共享。我们国赛期间用的酷盘，其他的没试过，可能有更好的东西，大家可以去摸索。 编程的人需要做的不仅仅是这些，比赛前做好万全的准备才能有条不紊。等到比赛时再去到处找代码写代码，那就晚了。好的习惯就是平时将各种代码搜集整理好，我这里有个列表仅供参考，如有疏漏，恳求指出： l 规划&优化（lingo）：0-1规划、线性规划、整数规划、非线性规划、动态规划、单目标、多目标、 l 图论：最短路径、hamilton圈、旅行商TSP问题 、最小生成树、网络最大流、最小费用流、 l 插值拟合 ：插值、线性拟合 、非线性拟合、最小二乘拟合 l 概率论&数理统计：概率模型、方差分析、回归分析（二次曲线回归，线性回归）、 假设检验、分布拟合检验、参数估计 l 微分方程：常微分方程、微分方程组、稳定状态、灵敏度分析 l 差分方程： l 时间序列： l 马氏链： l 聚类分析 ： l 智能算法 ：神经网络、遗传算法（gatool）、模拟退火 l 排队论： l 判别分析： l 生存数据分析： l 综合评价：层次分析、综合评分法、综合指数法、Topsis法、秩和比法 l 预测：灰色预测 l 系统仿真：蒙特卡洛 l 模糊数学：模糊聚类、模糊综合评价 l 图像处理：灰度化、二值化、滤波、边缘提取、三维重建 l 数据处理：主成分分析、因子分析 l 解方程： ...... 最好是Lingo和Matlab代码都准备好，尽量收集到并提前做好测试，相关软件前面已经提过，不再赘述。 以及各种数据处理方法： a) 回归分析法(数理统计方法)-用于对函数f（x）的一组观测值（xi,fi）i=1,2,…,n，确定函数的表达式。 b) 时序分析法--处理的是动态的时间序列相关数据，又称为过程统计方法。 c) 多元统计分析（聚类分析、判别分析、因子分析、主成分分析、生存数据分析）。 还有图形绘制：条形图、折线图、散点图、饼图、频率分布直方图...... 当然，资料和数据收集地址也得提前找好，当初比赛前我提前下载了1985-2009年的国家统计年鉴，费了不少流量。现在将有关资料放出来： 资料检索： n CNKI入口（西电图书馆的CNKI账号已经过期，只有出此下策了） http://202.119.208.220:8002/kns50/index.aspx http://61.155.19.94:8011/kns50/ http://cnki1.sztsg.net/kns50 n OA开放图书馆 http://www.oalib.com/search.asp?q=&sa=OA%E5%86%85%E5%AE%B9%E6%90%9C%E7%B4%A2 n 万方数据库（期刊，论文） n 超星图书馆（图书馆入口） n ，豆丁（下载器） l 数据检索 n 中国教育统计网 http://www.stats.edu.cn/ n 中国统计年鉴 http://www.stats.gov.cn/tjsj/ndsj/ n 国家数据统计库 http://219.235.129.58/welcome.do n 中国宏观数据挖掘分析系统 http://number.cnki.net/cyfd/ n 中国基础教育网 http://www.cbe21.com/subject/physics/ n 中国证券网 http://www.cnstock.com/ n 中国科学气象数据网 http://cdc.cma.gov.cn/ n 中国科学文献服务系统 http://sdb.csdl.ac.cn/ n 中国引文数据库 http://ref.cnki.net/knsref/index.aspx n 中科院科学数据库 http://www.csdb.cn/ n 中国动物主题数据库 http://www.zoology.csdb.cn/ n 中国统计年鉴下载： http://lib.njue.edu.cn/libtool/data.htm 最后，说明一点，组队的时候，最好能有两个会变编程的队友，这样遇到像公交车调度的问题，以免出现孤军奋战的局面。两个人可以进行思想的碰撞，必然会提高效率。 希望以上愚见对诸位希望参加数学建模，畏惧数模编程的童鞋起到一丁点的帮助。 祝参加数学建模的童鞋都能取得好成绩~~~~ （初稿。后续可能加入更多经验，尽请期待） Matlab .M文件编译成可执行文件.exe 分类： Matlab2011-08-24 15:42 363人阅读 评论(0) 收藏 举报 如何将MATLAB程序编译成独立可执行的程序？如何将编译好的独立可执行程序发布在没有安装MATLAB的电脑上？下面将一步步实现： 一、生成独立可执行的程序（exe文件）步骤 1、安装编译器。可有多种选择，matlab自带了一个LCC，推荐使用VC++6.0，我基于VS 2003实现。 2、设置编译器。在matlab命令行输入mbuild –setup以及mex –setup，选择安装的c编译器。 3、调用编译器。此处使用MATLAB下的一个GUI平台deploytool下完全实现。在命令窗口输入deploytool即可看到。具体使用方法请Help。 当然，也可以输入mcc -m filaname， filaname为要转成exe的m文件； 注：在以前的版本中，用编译命令mcc -B sglcpp filaname；自2006的版本后，替换为mcc -mfilaname； 4、安装 二、脱离matlab运行可执行程序 MCR是由matlab的运行环境，占用不到300M的对于用不同matlab版本生成的exe文件，MCR版本也会有不同，因此，在程序打包时，最好将相应版本的MCR一起打包。MCR环境的设置文件存放目录如下： 文件名为MCRInstaller.exe。可将其拷贝到自己的文件夹中，(7.0以前的版本是mglinstaller.exe)。 在MATLAB里运行可执行程序的办法是在前面加一个！，比如：!picshow，后缀名可有可无。 在其它没有安装matlab的机器上运行exe文件前： 首先安装matlab的运行环境。在同一机器上可以并存不同版本的matlab环境(换句话说不同版本不兼容)。 其次是要将―MCRinstaller.exe安装目录\\runtime\\win32‖这个路径添加到该计算机的环境变量中，通常是自动加载。 如果没有，也可手动安装，添加的方法是： 右击―我的电脑‖―属性‖―高级‖―环境变量‖―添加‖指定一个变量名，然后将上述路径复制到里面就可以了。 注：在安装过程中会弹出让安装Microsoft.NETFramework可以不用安装。 最后就是将编译生成的相相关文件拷贝到同一目录下，双击即可运行。 问题：目前此方法可完全运行在没有安装MATLAB以及C/C++的电脑上，但是如果是在AMD的CPU可以运行，但是不会出现任何MATLAB编译的界面。 美中不足就是，运行的时候dos的那个黑色地窗口一直存在。下面将实现去除黑屏的办法： 消除运行MATLAB生成的exe程序的dos黑屏的办法 基于MATLAB生成exe文件后，每次运行都存在dos黑屏的问题，现在可以通过以下方法解决： 方法一： 在命令窗口输入 cd(prefdir) edit compopts.bat 在打开的文件最后添加以下语句： A.VC环境下： set LINKFLAGS=%LINKFLAGS%/SUBSYSTEM:WINDOWS /ENTRY:mainCRTStartup B.LCC环境下: set LINKFLAGS=%LINKFLAGS% -subsystemwindows C. Borland: set LINKFLAGS=%LINKFLAGS% -aa 保存以后，再重新编译m文件，生成的exe文件运行起来就没有dos窗口了 方法二：使用suppress工具： 下载附件中的suppress压缩包后解压，(当然您可以自己去Google然后再下载)会看到一个suppress.ini文件，用记事本打开，然后将其中―Name=test.exe‖中text.exe的改为你生成的exe文件名。将suppress.exe(有个关盘和显示器的图标)，改后的 suppress.ini放到你生成的exe的同目录下。执行suppress.exe或者您自己生成的exe可以了。当然您可以自己修改 suppress.exe的名字，改为您自己想要的名字。 其中的方法一在使用后生成的exe再到没有任何安装MATLAB的机子上运行也不会有黑屏了。 方法二的缺点就是要同时存在您生成的exe以及supress.exe，必须在同一目录下。 为了您的安全，请只打开来源可靠的网址 打开网站 取消 来自: http://hi.baidu.com/lzhais/blog/item/18baf20a8beb8237b1351d1c.html 如何使用MATLAB进行USB2.0摄像头的编程 分类： Matlab2011-08-23 22:05 125人阅读 评论(0) 收藏 举报 Matlab中的图像获取工具箱给我们提供了必要的函数，我们直接调用就可以了。在在这帖中我们主要就是简单的介绍如何使用该工具箱进行对USB2.0摄像头的编程 废话不多说，我们开始言归正传了。但是一定记住你必须安装了PC摄像头才可以进行下面的东西，如果说首次安装摄像头最好重启下PC，否则可能出现没法识别摄像头。 整个过程我们需要做如下几件事情： 1、查询USB2.0Camera 的具体参数(imaqhwinfo) 2、创建视频输入对象(videoinput) 3、图像预览和显示(preview、stoppreview、closepreview和image) 4、获取视频图像(getsnapshot) 5、图像获取设备的获取和设置(get和set) 6、关闭视频对象(delete) 在正式讲解之前，我想说明下几个个在图像获取工具箱中的术语： 图像获取设备：比如摄像头、扫描仪 图像获取适配器：主要的目的是通过驱动在Matlab和图像获取设备之间传递信息 ROI：region-of-interest 感兴趣区域 在说说几个常用的函数，我们这里只是说明它的作用，具体如何使用参考帮助系统 getselectedsource imaqfind isvalid peekdata getdata imaqmontage 给我们一个摄像头我们必须知道他的相关参数，才可能进行我们的编程下。当然我们可以查询 商家手册，但是那个累不累人呀。 Matlab的图像获取工具箱为我提供了imaqhwinfo()，来获取PC上以安装的图像获取硬件信息 没有输入参数时，返回一个结构体， 它包含了系统中存在的适配器和Matlab相关的版本信息 (第一次我们一般使用这个) 代码: >> info=imaqhwinfo info = InstalledAdaptors: {'coreco' 'winvideo'}%这里可以看到我的PC上安装了两个适配器 MATLABVersion: '7.6 (R2008a)' ToolboxName: 'Image Acquisition Toolbox' ToolboxVersion: '3.1 (R2008a)' 有输入参数的时候，返回一个结构体，包含了指定的适配器的数据信息 代码: >> win_info=imaqhwinfo('winvideo')%我们看看第二适配器的具体参数 win_info = AdaptorDllName: [1x81 char]%适配器dll文件绝对路径 AdaptorDllVersion: '3.1 (R2008a)'%适配器dll文件版本 AdaptorName: 'winvideo'%s适配器名称 DeviceIDs: {[1]}%设备ID号，这个我们经常需要用到 DeviceInfo: [1x1 struct]%设备信息，这里主要是图像获取设备的一些参数，比较重要 %====================下面我们了解下，这个图像获取设备到底有哪些的详细信息吧 ==================== >> win_info.DeviceIDs ans = [1] >> dev_win_info=win_info.DeviceInfo dev_win_info = DefaultFormat: 'RGB24_320x240'%获取图片的默认格式 DeviceFileSupported: 0 DeviceName: 'USB PC CAMERA P227'%设备名称 DeviceID: 1%设备号 ObjectConstructor: 'videoinput('winvideo', 1)'%对象构建方式，这个绝大部分都是一样的 SupportedFormats: {1x12 cell}%获取的图像支持格式，一般都有好多种，上面的 DefaultFormat只是默认格式而已 %==================================看看图像获取设备支持的图像格式 ================================== >> dev_win_info.SupportedFormats%可以看到我的PC上的摄像头支持下面12中图片格式 ans = Columns 1 through 5 'I420_160x120' 'I420_176x144' 'I420_320x240' 'I420_352x288' 'I420_640x480' Columns 6 through 9 'RGB24_1280x960' 'RGB24_160x120' 'RGB24_176x144' 'RGB24_320x240' Columns 10 through 12 'RGB24_352x288' 'RGB24_640x480' 'RGB24_800x600' 视频预览、采集和保存 (1)创建视频输入对象 obj = videoinput(adaptorname,deviceID,format) adaptorname：适配器名称，首次可以使用不带参数的imaqhwinfo函数获取 deviceID：设备ID号，首次可以通过imaqhwinfo函数获取 format：视频采集格式，可以通过DeviceInfo的SupportedFormats获取，不填写则使用默认 格式 代码: >> obj = videoinput('winvideo',1,'RGB24_320x240')%这里我们使用默认的视频采集格式 Summary of Video Input Object Using 'USB PC CAMERA P227'. Acquisition Source(s): input1 is available. Acquisition Parameters: 'input1' is the current selected source. 10 frames per trigger using the selected source. 'RGB24_320x240' video data to be logged upon START. Grabbing first of every 1 frame(s). Log data to 'memory' on trigger. Trigger Parameters: 1 'immediate' trigger(s) on START. Status: Waiting for START. 0 frames acquired since starting. 0 frames available for GETDATA. (2)打开视频预览窗口 himage=preview(obj,himage) obj：视频采集对象 himage：视频预览窗口对应的句柄，也就是说在指定的句柄对象中预览视频，该参数可以空 缺 至于预览窗口的关闭和停止可以使用colsepreview和stoppreview函数 代码: vidRes = get(obj, 'VideoResolution'); nBands = get(obj, 'NumberOfBands'); figure()%指定预览窗体显示的figure axes()%指定预览窗口显示的坐标系 hImage = image( zeros(vidRes(2), vidRes(1), nBands) ); preview(obj, hImage); (3)图像捕捉、显示和保存 代码: %frame是H×W×B的矩阵。H图像高度，由ROIPosition指定；w图像宽度，由ROIPosition 指定；B索线个数，由NumberOfBands指定 frame = getsnapshot(obj); imshow(frame); imwrite(fame,'snap.jpg','jpg'); 实验步骤 1、 查询USB2.0Camera 的具体参数 输入 imaqInfo = imaqhwinfo 返回信息 InstalledAdaptors: {'winvideo'} MATLABVersion: '7.1 (R14SP3)' ToolboxName: 'Image Acquisition Toolbox' ToolboxVersion: '1.9 (R14SP3)' 输入imaqInfo.InstalledAdaptors 返回信息 ans = 'winvideo' 输入winvideoinfo = imaqhwinfo('winvideo') 返回信息 winvideoinfo = AdaptorDllName: [1x76 char] AdaptorDllVersion: '1.9 (R14SP3)' AdaptorName: 'winvideo' DeviceIDs: {[1]} DeviceInfo: [1x1 struct] 输入 winvideoinfo.DeviceInfo 返回信息 ans = DefaultFormat: 'YUY2_160x120' DeviceFileSupported: 0 DeviceName: 'USB 视频设备' DeviceID: 1 ObjectConstructor: 'videoinput('winvideo', 1)' SupportedFormats: {1x5 cell} 输入device1 = winvideoinfo.DeviceInfo(1) 返回信息 device1 = DefaultFormat: 'YUY2_160x120' DeviceFileSupported: 0 DeviceName: 'USB 视频设备' DeviceID: 1 ObjectConstructor: 'videoinput('winvideo', 1)' SupportedFormats: {1x5 cell} 输入device1.DeviceName 返回信息 ans = USB 视频设备 输入device1.DeviceID 返回信息 ans = 1 输入device1.DefaultFormat 返回信息 ans = YUY2_160x120 输入device1.SupportedFormats 返回信息 Columns 1 through 4 'YUY2_160x120' 'YUY2_176x144' 'YUY2_320x240' 'YUY2_352x288' Column 5 'YUY2_640x480' 2、 最简单采集实验,输入如下代码，可以得到预览下的默认格式的摄像头捕捉窗口obj=videoinput('winvideo',1); preview(obj); 3、 输入如下代码 vidobj = videoinput('winvideo',1,'YUY2_640x480'); sources = vidobj.Source; whos sources set(vidobj,'SelectedSourceName','input1'); sources selectedsrc = getselectedsource(vidobj); get(selectedsrc); delete(vidobj); clear vidobj; 返回信息 Name Size Bytes Class sources 1x1 726 videosource object Grand total is 30 elements using 726 bytes Display Summary for Video Source Object: Index: SourceName: Selected: 1 'input1' 'on' General Settings: Parent = [1x1 videoinput] Selected = on SourceName = input1Tag = Type = videosource Device Specific Properties: BacklightCompensation = on Brightness = -16 Contrast = 120 FrameRate = 30.0000 Gamma = 60 Hue = 0 Saturation = 40 Sharpness = 3 4、 输入如下代码 clc; clf; clear all; imaqmem(30000000); %申请内存空间 %ADAPTOR:MATLAB与视频设备之间的接口，主要的目的是传递信息 vid = videoinput('winvideo', 1, 'YUY2_640x480'); preview(vid); start(vid); h=figure('NumberTitle','off','Name','视频',... 'MenuBar','none','color','c',... 'Position', [0, 0, 1, 1], 'Visible', 'on'); %新建窗口 set(h,'doublebuffer','on','outerposition',get(0,'screensize')); h1=axes('Position', [0.02, 0.1, 0.4, 0.8],'Parent',h); %新建显示窗口 hold on; axis off; while ishandle(h) %判断是否有效的图像对象句柄 a=getsnapshot (vid); % 捕获图像 flushdata(vid); %清除数据获取引擎的所有数据、置属性SamplesAvailable为0 imshow(a); %显示图像 drawnow; % 实时更新图像 end; delete(vid); [原创]Matlab模拟时针和数字时针 分类： Matlab2011-07-22 00:27 90人阅读 评论(0) 收藏 举报 [记录数模集训期间Matlab的学习过程] ①模拟时针 [plain] view plaincopy 1. clear,close all,clc 2. degree=[0:0.01:pi*2]; 3. plot(cos(degree),sin(degree),'r-'); 4. set(gca,'Xtick',[],'Ytick',[]);box on; 5. axis([-1.3 1.3 -1.3 1.3]) 6. axis square 7. title('模拟时钟','Fontsize',22,'Fontname','华文行楷') 8. %时针刻度 9. for m=5*pi/2:-pi/6:pi/2 10. a=line([9*cos(m)/10,cos(m)],[9*sin(m)/10,sin(m)],'color','r'); 11. if(m~=5*pi/2) 12. text(17*cos(m)/20,17*sin(m)/20,int2str(12-(6*m-3*pi)/pi),'Fontsize',10,'Fontname','Times New Roman','Color','B','HorizontalAlignment','center'); 13. end 14. pause(0.005) 15. end 16. text(0,-0.2,'CopyRight By D.C.','Fontsize',8,'Fontname','Times New Roman','Color','B','HorizontalAlignment','center'); 17. %分针刻度 18. for n=0:pi/30:2*pi 19. if mod(n,pi/6)~=0 20. a=line([19*cos(n)/20,cos(n)],[19*sin(n)/20,sin(n)]); 21. end 22. pause(0.005) 23. end 24. clear m,clear n; 25. %秒针 26. while (1) 27. timeNow=clock; 28. timeN=[2*pi-timeNow(4:6)/30*pi+pi/2]; 29. %hou=2*pi-timeNow(5)/30*pi+pi/2; 30. %mini=2*pi-timeNow(5)/30*pi+pi/2; 31. %sec=2*pi-timeNow(6)/30*pi+pi/2; 32. h=timeN(1); 33. i=timeN(2); 34. j=timeN(3); 35. %画时针 36. %ah=line([0,18*cos(h)/20],[0 18*sin(h)/20]); 37. bh=line([1.5*cos(h+pi/2)/20 -3*cos(h)/20],[1.5*sin(h+pi/2)/20 -3*sin(h)/20]); 38. ch=line([1.5*cos(h-pi/2)/20 -3*cos(h)/20],[1.5*sin(h-pi/2)/20 -3*sin(h)/20]); 39. dh=line([1.5*cos(h+pi/2)/20 14*cos(h)/20],[1.5*sin(h+pi/2)/20 14*sin(h)/20]) ; 40. eh=line([1.5*cos(h-pi/2)/20 14*cos(h)/20],[1.5*sin(h-pi/2)/20 14*sin(h)/20]); 41. %画分针 42. a=line([-4*cos(i)/20,17*cos(i)/20],[-4*sin(i)/20 17*sin(i)/20]); 43. b=line([1*cos(i+pi/2)/20 -4*cos(i)/20],[1*sin(i+pi/2)/20 -4*sin(i)/20]); 44. c=line([1*cos(i-pi/2)/20 -4*cos(i)/20],[1*sin(i-pi/2)/20 -4*sin(i)/20]); 45. d=line([1*cos(i+pi/2)/20 17*cos(i)/20],[1*sin(i+pi/2)/20 17*sin(i)/20]); 46. e=line([1*cos(i-pi/2)/20 17*cos(i)/20],[1*sin(i-pi/2)/20 17*sin(i)/20]); 47. %画秒针 48. f=line([-3*cos(j)/20,18*cos(j)/20],[-3*sin(j)/20 18*sin(j)/20]); 49. %定时检测 50. pause(1) 51. delete(bh),delete(ch),delete(dh),delete(eh) 52. delete(a),delete(b),delete(c),delete(d),delete(e),delete(f) 53. end ----------------------------------------------------------------------可恶的分隔线---------------------------------------------------------------------------------- ②数字时针 [plain] view plaincopy 1. while 1 2. c=clock; 3. axes('Position',[0.06,0.4,0.86,0.4]); 4. title('数字时钟','Fontsize',22,'Fontname','华文行楷') 5. set(gca,'Xtick',[],'Ytick',[]);box on; 6. str=[int2strNew(c(4)),':',int2strNew(c(5)),':',int2strNew(c(6))]; 7. Ts=text(0.5,0.5,str,'Fontsize',80,... 8. 'Fontname','Times New Roman','Color','B',... 9. 'HorizontalAlignment','center'); 10. pause(1) 11. clf 12. end 13. 14. %--------------------------------- 15. 16. %int2strNew.m 17. 18. %---------------------------------- 19. 20. function y=int2strNew(a) 21. if size(int2str(a),2)==1 22. y=['0',int2str(a)]; 23. else 24. y=int2str(a); 25. end %----------------------------------------------------------可恶的分隔线-------------------------------------------------------------- CopyRight By D.C. 2011.7.22 0:21 转载请注明出处 数模国赛前一天总结 分类： 数学建模2011-11-16 23:09 271人阅读 评论(0) 收藏 举报 说明：以下部分链接不到的情况原因是本人超链接到相应文件夹下面，谢谢支持 ============================================================================================================ Ø 方法&工具总结 参考资料:算法大全 l 规划&优化（lingo） 1. 0-1规划 2. 线性规划（lingo.exe） 3. 整数规划（lingo代码） 4. 非线性规划 5. 动态规划 6. 单目标 7. 多目标 l 图论 1. 最短路径（Floyd）（图论.exe） 2. hamilton圈 3. 旅行商TSP问题 （matlab）（lingo） 4. 最小生成树（Kruskal）（prim） 5. 网络最大流（matlab）（lingo） 6. 最小费用流 l 插值拟合 （1stOpt） 1. 插值（Matlab） 2. 线性拟合 （1stOpt） 3. 非线性拟合 4. 最小二乘拟合 l 概率论&数理统计（SPSS,eviews） 1. 概率模型 2. 方差分析 3. 回归分析：二次曲线回归，线性回归 4. 假设检验 5. 分布拟合检验 6. 参数估计 l 微分方程（matlab） 1. 常微分方程 2. 微分方程组 3. 稳定状态 4. 灵敏度分析 l 差分方程 l 时间序列（eviews） l 马氏链 l 聚类分析 （SPSS,matlab） l 智能算法 （matlab优化工具箱） 1. 神经网络 2. 遗传算法（gatool） 3. 模拟退火 l 排队论 l 判别分析 l 生存数据分析 l 综合评价 1. 层次分析（AHP.exe） 2. 综合评分法 3. 综合指数法 4. Topsis法 5. 秩和比法 l 预测 1. 灰色预测 l 系统仿真 1. 蒙特卡洛 l 模糊数学 1. 模糊聚类 2. 模糊综合评价 l 图像处理 1. 灰度化 2. 二值化 3. 滤波 4. 边缘提取 5. 三维重建 l 数据处理 1. 主成分分析 2. 因子分析 l 解方程 ============================================================================================================ Ø 处理方法 u 数据分析法 对大量的观测数据进行统计分析，寻求规律建立数学模型，采用的分析方法一般有： a) 回归分析法(数理统计方法)-用于对函数f（x）的一组观测值（xi,fi）i=1,2,…,n，确定函数的表达式。 b) 时序分析法--处理的是动态的时间序列相关数据，又称为过程统计方法。 c) 多元统计分析（聚类分析、判别分析、因子分析、主成分分析、生存数据分析）。 u 图形绘制 ² 类型 1. 条形图 2. 折线图 3. 散点图 4. 饼图 5. 频率分布直方图 1. ² 线性（虚线、实线） 2. ² 颜色 3. ² 点标记 4. ² 坐标轴控制 5. ² legend 6. ² 标题 7. ² grid 8. ² box 9. ² ticks 10. ² 字体 ============================================================================================================ Ø 数学软件 编程计算 概率统计 Matlab SPSS Lingo/Lindo SAS Maple Eviews Mathematica Origion Excel Vensim 优化 绘图 1stOpt SigmaPlot WinQSB ScienceWord SmartDraw 几何画板 autoCAD 其他小众软件：图论、AHP… 其他功能软件：tortoiseSVN（版本控制）、dexpot（虚拟桌面）、酷盘（局域网共享协作） 其他：ansys(有限元分析)、comsol、FLAC2D/FLAC3D、 ============================================================================================================ Ø 数据&资源检索l 资料检索 n CNKI入口 1. http://202.119.208.220:8002/kns50/index.aspx 2. http://61.155.19.94:8011/kns50/ 3. http://cnki1.sztsg.net/kns50 n OA开放图书馆 http://www.oalib.com/search.asp?q=&sa=OA%E5%86%85%E5%AE%B9%E6%90%9C%E7%B4%A2 n 万方数据库（期刊，论文） n 超星图书馆（图书馆入口） n ，豆丁（下载器） l 数据检索 n 中国教育统计网 http://www.stats.edu.cn/ n 中国统计年鉴 http://www.stats.gov.cn/tjsj/ndsj/ n 国家数据统计库 http://219.235.129.58/welcome.do n 中国宏观数据挖掘分析系统 http://number.cnki.net/cyfd/ n 中国基础教育网 http://www.cbe21.com/subject/physics/ n 中国证券网 http://www.cnstock.com/ n 中国科学气象数据网 http://cdc.cma.gov.cn/ n 中国科学文献服务系统 http://sdb.csdl.ac.cn/ n 中国引文数据库 http://ref.cnki.net/knsref/index.aspx n 中科院科学数据库 http://www.csdb.cn/ n 中国动物主题数据库 http://www.zoology.csdb.cn/ n 中国统计年鉴下载： http://lib.njue.edu.cn/libtool/data.htm 中国统计年鉴2009 中国统计年鉴2008 中国统计年鉴2007 中国统计年鉴2006 中国统计年鉴2005 中国统计年鉴2004 中国统计年鉴2003 中国统计年鉴2002 中国统计年鉴2001 中国统计年鉴2000 中国统计年鉴1999 中国统计年鉴1998 中国统计年鉴1997 中国统计年鉴1996 中国统计年鉴1995 中国统计年鉴1994 中国统计年鉴1993 中国统计年鉴1992 中国统计年鉴1991 中国统计年鉴1990 中国统计年鉴1989 中国统计年鉴1988 中国统计年鉴1987 中国统计年鉴1986 中国统计年鉴1985 统计数据使用说明 ============================================================================================================ Ø 注意事项： 1. n 资料备份 2. n 硬件 3. n 软件 4. n 场地 5. n 食物 =====================================以下安排摘自互联网，仅供参考============================================== Ø 分工&流程 第一天 一、上午 1. 各自独立思考1个小时，主要分析题目的问题背景，已知条件，建模目的等问题。至少每人必须提出10到15个问题，并回答自己的问题。 2. 重点用语言的形式表述清楚问题的结构，即用语言描述自己的初步模型。（要自己提出的模型，可能就会产生一些假设。） 3. 再和队友讨论。讨论1个小时。形成自己团队的初步模型，同样是以语言形式描述的。 4. 接下来查找一些文献，讨论修改团队的模型，形成一个最终较完整的模型。并根据讨论最后形成对问题的统一认识，形成问题重述部分的内容。 注：1）如果问题有好几问，可以重点讨论第一个问题，但是也要考虑其他问题与第一问的关系！（一般建模中的几问都是有一定联系得）；也可以同时考虑，同时建模。 2）注意参考文献的处理，参考别人的方法一定要在文中注明！这也是要求一直留意查找文献的目的。 二、下午 将自己团队的模型数学化，用数学符号和数学语言公式的形式，表述自己的模型。此时会继续需要查文献，产生一些假设条件，并产生自己论文中的符号说明。 第二天 一、上午 一个人开始写文章，语言重在逻辑清晰，叙述简洁明了！图、表准确。文章格式正确、内容完整。（问题重述，问题分析，模型假设，符号说明，模型形式，以及参考文献都已经在第一天的讨论中有了一定的共识。） 其余两个人（在不清楚时3人讨论），开始考虑第一个问题的模型的求解，即研究模型的解法。查找文献或者自己提出对模型的求解方法。此时可能需要继续对第一天建立的模型进行修改，简化等处理。（讨论后，及时告诉写文章的队友）。 二、下午 写文章的继续。 编程的开始编程计算模型。此时，可能需要根据所采取的算法对模型的表述重新修改。 另一人帮忙编程，并开始考虑第二个、第三个问题的模型及求解方法。并一起讨论，形成共识，写进文章中。（此时，同样可能需要查文献，符号表示，产生假设） 第三天 一、上午 应该给出所有问题的计算结果了（最迟下午6点前）。 产生论文初稿。 二、下午 进行模型的分析。主要是分析编程计算出的解的现实意义等，通过图、表等形式说明自己的结果。并一定进行误差分析（因为模型是对实际问题的近似，同时在建模中也进行的假设，所以必须进行误差分析。） 注：如果模型的计算在下午才出来，需要加紧进度。晚上不要休息了！ 三、晚上 对模型进行总结推广（３人讨论１个小时，切忌不要在这个问题上过多的讨论，只需写一段。只讨论模型本身的问题，假设的合理性去处和条件的放松，模型的求解方法等。文章的此部分必须有，但是一定不能太多。） 重点精力放在对模型的摘要的书写上，一定3人认真讨论2小时左右。摘要A4纸2/3，主要是模型的目标，方法，结果。用清晰简洁的语言叙述，突出创新的内容。 注： 1）整建模过程中要注意自己数据，文章的电子文件的保存，随时保存副本！ 2）队内交流，不可队外交流，不要和其他队和人交流！以免雷同，抄袭的发生！保护自己的劳动成果。 3）假设要认真考虑，切合实际，又合理，同时，可以使处理的问题简单化。一定不要为了假设而假设（即为了论文中有模型假设这一内容，而做出一些无意的假设!）。 重要算法代码搜集（提前测试） 因篇幅问题不能全部显示，请点此查看更多更全内容function f=zhengtai(c,x) 2. f=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c(2)^2));
clear 2. clc
function judge(P,c,rela) 2. %评价指标 并显示比较结果
%%%%%%%%jiaohuan3.m%%%%%%%%%% 2. n=0;