图论算法及Matlab程序代码
图论常用算法matlab程序
运筹学算法matlab程序西北工业大学数学系2009级1.顺向Dijkstra 算法M=[ 0 5 9 Inf Inf Inf InfInf 0 Inf Inf 12 Inf InfInf 3 0 15 Inf 23 InfInf 6 Inf 0 Inf 8 7Inf 12 Inf 5 0 Inf 14Inf Inf Inf Inf Inf 0 10Inf Inf Inf Inf Inf Inf 0];first=1;last=7;[m,n]=size(M);L=zeros(1,m);symbol=zeros(1,m);direction=zeros(1,m);for i=1:mif(i~=first)L(i)=inf;enddirection(i)=first;endjudge=1;while judgefor i=1:mif(symbol(i)==0)min=L(i);temporary=i;breakendendfor i=1:mif(symbol(i)==0)if(L(i)<min)min=L(i);temporary=i;endendendk=temporary;for j=1:mif(symbol(1,j)==0)if(M(k,j)==inf)continue;elseif(L(k)+M(k,j)<L(j))L(j)=L(k)+M(k,j);direction(j)=k;endendendendsymbol(k)=1;num=0;for i=1:mif(symbol(i)==1)num=num+1;endendif(num==m)judge=0;endendp=last;arrow=zeros(1,m);arrow(1)=last;i=2;while p~=firstarrow(1,i)=direction(p);i=i+1;p=direction(p);enddistance=L(last);M=[ 0 5 9 Inf Inf Inf Inf Inf 0 Inf Inf 12 Inf InfInf 3 0 15 Inf 23 Inf Inf 6 Inf 0 Inf 8 7 Inf 12 Inf 5 0 Inf 14 Inf Inf Inf Inf Inf 0 10Inf Inf Inf Inf Inf Inf 0]; [m,n]=size(M);first=1;last=7;L=zeros(1,m);direction=zeros(1,m);symbol=zeros(1,m);for i=1:mdirection(i)=last;if(i~=last)L(i)=inf;endendjudge=1;while judgefor i=1:mif(symbol(i)==0)min=L(i);temporary=i;breakendendfor i=1:mif(symbol(i)==0)if(L(i)<min)min=L(i);temporary=i;endendendk=temporary;for i=1:mif(M(i,k)==inf)continueelseif(M(i,k)+L(k)<L(i))L(i)=L(k)+M(i,k);direction(i)=k;endendendsymbol(k)=1;sum=0;for i=1:mif(symbol(i)==1)sum=sum+1;endendif(sum==m)judge=0;endendp=first;i=2;arrow=zeros(1,m);arrow(1)=first;while p~=lastarrow(i)=direction(p);i=i+1;p=direction(p);endd=[0 7 5 12 inf infinf 0 inf 3 inf infinf inf 0 6 inf 1512 inf 6 0 inf 86 inf 13 inf 0 infinf 4 15 inf 9 0];[m,n]=size(d);p=zeros(m,n);for i=1:np(:,i)=i;endfor k=1:nfor i=1:mfor j=1:nif(d(i,k)+d(k,j)<d(i,j))d(i,j)=d(i,k)+d(k,j);p(i,j)=p(i,k);endendendend4.仿floyd 算法d=[inf 6 0 4 0 0 00 inf 0 0 5 0 04 7 inf 0 05 00 0 4 inf 0 3 00 0 2 0 inf 0 00 0 0 0 4 inf 50 0 0 0 6 0 inf];[m,n]=size(d);first=1;last=7;direction=zeros(m,m);for i=1:mdirection(:,i)=i;endfor i=1:mfor j=1:mfor k=1:msmall=min(d(i,k),d(k,j));if d(i,j)<smalld(i,j)=small;direction(i,j)=direction(i,k);endendendendarrow=zeros(1,m);arrow(1)=first;i=2;p=first;while p~=lastp=direction(p,last);arrow(i)=p;i=i+1;end—dijkstra算法d=[0 inf 3 5 inf10 0 14 inf 8inf inf 0 7 -6inf inf inf 0 infinf inf inf -1 0];[m,n]=size(d);first=2;last=4;L=zeros(1,n);z=zeros(m,n);symbol=zeros(1,n);direction=zeros(1,n);for i=1:nfor j=1:mif d(i,j)~=0if d(i,j)~=infz(i,j)=1;endendenddirection(i)=first;if i~=firstL(i)=inf;endendjudge=1;while judgemini=10;for j=1:nif symbol(j)==0sum=0;for i=1:mp=z(i,j)*(1-symbol(i));sum=sum+p;endif(sum==0)mini=j;breakendendendfor j=1:nif symbol(j)==0&&z(mini,j)==1if L(mini)+d(mini,j)<L(j)L(j)=L(mini)+d(mini,j);direction(j)=mini;endendendsymbol(mini)=1;num=0;for i=1:nif symbol(i)==1num=num+1;endendif num==m;judge=0;endendarrow=zeros(1,m);p=last;arrow(1)=last;i=2;while p~=firstp=direction(p);arrow(i)=p;i=i+1;end—dijkstra算法d=[0 inf 3 5 inf10 0 14 inf 8inf inf 0 7 -6inf inf inf 0 infinf inf inf -1 0];[m,n]=size(d);first=2;last=4;L=zeros(1,n);z=zeros(m,n);symbol=zeros(1,n);direction=zeros(1,n);for i=1:nfor j=1:mif d(i,j)~=0if d(i,j)~=infz(i,j)=1;endendenddirection(i)=last;if i~=lastL(i)=inf;endendjudge=1;while judgemini=10;for i=1:nif symbol(i)==0sum=0;for j=1:mp=z(i,j)*(1-symbol(j));sum=sum+p;endif(sum==0)mini=i;breakendendendfor i=1:nif symbol(i)==0&&z(i,mini)==1if L(mini)+d(i,mini)<L(i)L(i)=L(mini)+d(i,mini);direction(i)=mini;endendendsymbol(mini)=1;num=0;for i=1:nif symbol(i)==1num=num+1;endendif num==m;judge=0;endendarrow=zeros(1,m);p=first;arrow(1)=first;i=2;while p~=lastp=direction(p);arrow(i)=p;i=i+1;endM=[ 0 17 11 inf inf inf17 0 13 12 28 1511 13 0 inf 19 infinf 12 inf 0 inf 16inf 28 19 inf 0 10inf 15 inf 16 10 0];[m,n]=size(M);X=zeros(m,n);Y=zeros(m);Z=zeros(m);Y(1)=1;for i=2:mZ(i)=i;endjudge=1;while judgefor i=1:mif(Y(i)~=0)for j=1:mif(Z(j)~=0)min=M(i,j);a=i;b=j;endendendendfor i=1:mif(Y(i)~=0)for j=1:mif(Z(j)~=0)if(M(i,j)<min)min=M(i,j);a=i;b=j;endendendendendY(b)=b;Z(b)=0;X(a,b)=1;X(b,a)=1;c=0;for i=1:mif(Y(i)~=0)c=c+1;endendif(c==m)judge=0;endend网络最大流Ford—Fulkersen算法d=[inf 12 17 0 0 00 inf 0 8 0 00 6 inf 0 12 00 0 5 inf 0 150 0 0 4 inf 90 0 0 0 0 inf];[m,n]=size(d);X=zeros(m,n);first=1;last=6;recognize=1;while recognizeL=zeros(1,m);L(first)=inf;direction=ones(1,m);symbol=zeros(1,m);judge=1;while judgefor i=1:mif symbol(i)==0big=L(i);k=i;break;endendfor i=1:mif symbol(i)==0if L(i)>bigbig=L(i);k=i;endendendif k==nif L(n)==0breakendelsefor j=1:mif d(k,j)>0u=min(L(k),d(k,j)-X(k,j));if u>L(j)L(j)=u;direction(j)=k;endelseif d(j,k)>0u=min(L(k),X(j,k));if u>L(j)L(j)=u;direction(j)=k;endendendendendsymbol(k)=1;num=0;for i=1:mif symbol(i)==1num=num+1;endendif num==mjudge=0;endendafter=last;before=after;while before~=firstbefore=direction(after);if d(before,after)>0X(before,after)=X(before,after)+L(n); elseX(before,after)=X(before,after)-L(n); endafter=before;endif L(m)==0recognize=0;end end。
超全图论matlab程序-可解决图论方面的绝大多数问题
程序三:有向图关联矩阵和邻接矩阵互换算法
function W=mattransf(F,f) if f==0 m=sum(sum(F)); n=size(F,1); W=zeros(n,m); k=1; for i=1:n for j=i:n if F(i,j)~=0 W(i,k)=1; W(j,k)=-1; k=k+1; end end end elseif f==1 m=size(F,2); n=size(F,1); W=zeros(n,n); for i=1:m a=find(F(:,i)~=0); if F(a(1),i)==1 W(a(1),a(2))=1; else W(a(2),a(1))=1; end end else fprint('Please imput the right value of f'); end W;
第二讲:最短路问题
程序一:Dijkstra算法(计算两点间的最短路)
图论matlab有关程序
超全的图论程序程序一:可达矩阵算法function P=dgraf(A)n=size(A,1);P=A;for i=2:nP=P+A^i;endP(P~=0)=1;P;程序二:关联矩阵和邻接矩阵互换算法function W=incandadf(F,f)if f==0m=sum(sum(F))/2;n=size(F,1);W=zeros(n,m);k=1;for i=1:nfor j=i:nif F(i,j)~=0W(i,k)=1;W(j,k)=1;k=k+1;endendendelseif f==1m=size(F,2);n=size(F,1);W=zeros(n,n);for i=1:ma=find(F(:,i)~=0);W(a(1),a(2))=1;W(a(2),a(1))=1;endelsefprint('Please imput the right value of f');endW;程序三:有向图关联矩阵和邻接矩阵互换算法function W=mattransf(F,f)if f==0m=sum(sum(F));n=size(F,1);W=zeros(n,m);k=1;for i=1:nfor j=i:nif F(i,j)~=0W(i,k)=1;W(j,k)=-1;k=k+1;endendendelseif f==1m=size(F,2);n=size(F,1);W=zeros(n,n);for i=1:ma=find(F(:,i)~=0);if F(a(1),i)==1W(a(1),a(2))=1;elseW(a(2),a(1))=1;endendelsefprint('Please imput the right value of f'); endW;第二讲:最短路问题程序一:Dijkstra算法(计算两点间的最短路)function [l,z]=Dijkstra(W)n = size (W,1);for i = 1 :nl(i)=W(1,i);z(i)=0;endi=1;while i<=nfor j =1 :nif l(i)>l(j)+W(j,i)l(i)=l(j)+W(j,i);z(i)=j-1;if j<ii=j-1;endendendi=i+1;end程序二:floyd算法(计算任意两点间的最短距离)function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendr;for k=1:nfor i=1:nfor j=1:nif d(i,k)+d(k,j)<d(i,j)d(i,j)=d(i,k)+d(k,j); r(i,j)=r(i,k);endendendend程序三:n2short.m 计算指定两点间的最短距离function [P u]=n2short(W,k1,k2)n=length(W);U=W;m=1;while m<=nfor i=1:nfor j=1:nif U(i,j)>U(i,m)+U(m,j)U(i,j)=U(i,m)+U(m,j);endendendm=m+1;endu=U(k1,k2);P1=zeros(1,n);k=1;P1(k)=k2;V=ones(1,n)*inf;kk=k2;while kk~=k1for i=1:nV(1,i)=U(k1,kk)-W(i,kk);if V(1,i)==U(k1,i)P1(k+1)=i;kk=i;k=k+1;endendendk=1;wrow=find(P1~=0);for j=length(wrow):-1:1P(k)=P1(wrow(j));k=k+1;endP;程序四、n1short.m(计算某点到其它所有点的最短距离) function[Pm D]=n1short(W,k)n=size(W,1);D=zeros(1,n);for i=1:n[P d]=n2short(W,k,i);Pm{i}=P;D(i)=d;end程序五:pass2short.m(计算经过某两点的最短距离) function [P d]=pass2short(W,k1,k2,t1,t2)[p1 d1]=n2short(W,k1,t1);[p2 d2]=n2short(W,t1,t2);[p3 d3]=n2short(W,t2,k2);dt1=d1+d2+d3;[p4 d4]=n2short(W,k1,t2);[p5 d5]=n2short(W,t2,t1);[p6 d6]=n2short(W,t1,k2);dt2=d4+d5+d6;if dt1<dt2d=dt1;P=[p1 p2(2:length(p2)) p3(2:length(p3))]; elsed=dt1;p=[p4 p5(2:length(p5)) p6(2:length(p6))]; endP;d;第三讲:最小生成树程序一:最小生成树的Kruskal算法function [T c]=krusf(d,flag)if nargin==1n=size(d,2);m=sum(sum(d~=0))/2;b=zeros(3,m);k=1;for i=1:nfor j=(i+1):nif d(i,j)~=0b(1,k)=i;b(2,k)=j;b(3,k)=d(i,j);k=k+1;endendendelseb=d;endn=max(max(b(1:2,:)));m=size(b,2);[B,i]=sortrows(b',3);B=B';c=0;T=[];k=1;t=1:n;for i=1:mif t(B(1,i))~=t(B(2,i))T(1:2,k)=B(1:2,i);c=c+B(3,i);k=k+1;tmin=min(t(B(1,i)),t(B(2,i)));tmax=max(t(B(1,i)),t(B(2,i)));for j=1:nif t(j)==tmaxt(j)=tmin;endendendif k==nbreak;endendT;c;程序二:最小生成树的Prim算法function [T c]=Primf(a)l=length(a);a(a==0)=inf;k=1:l;listV(k)=0;listV(1)=1;e=1;while (e<l)min=inf;for i=1:lif listV(i)==1for j=1:lif listV(j)==0 & min>a(i,j)min=a(i,j);b=a(i,j);s=i;d=j;endendendendlistV(d)=1;distance(e)=b;source(e)=s;destination(e)=d;e=e+1;endT=[source;destination];for g=1:e-1c(g)=a(T(1,g),T(2,g));endc;另外两种程序最小生成树程序1(prim 算法构造最小生成树)a=[inf 50 60 inf inf inf inf;50 inf inf 65 40 inf inf;60 inf inf 52 inf inf 45;...inf 65 52 inf 50 30 42;inf 40 inf 50 inf 70 inf;inf inf inf 30 70 inf inf;...inf inf 45 42 inf inf inf];result=[];p=1;tb=2:length(a);while length(result)~=length(a)-1temp=a(p,tb);temp=temp(:);d=min(temp);[jb,kb]=find(a(p,tb)==d);j=p(jb(1));k=tb(kb(1));result=[result,[j;k;d]];p=[p,k];tb(find(tb==k))=[];endresult最小生成树程序2(Kruskal 算法构造最小生成树)clc;clear;a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40;a(3,4)=52;a(3,7)=45; a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70;[i,j,b]=find(a);data=[i';j';b'];index=data(1:2,:);loop=max(size(a))-1;result=[];while length(result)<looptemp=min(data(3,:));flag=find(data(3,:)==temp);flag=flag(1);v1=data(1,flag);v2=data(2,flag);if index(1,flag)~=index(2,flag)result=[result,data(:,flag)];endindex(find(index==v2))=v1;data(:,flag)=[];index(:,flag)=[];endresult第四讲:Euler图和Hamilton图程序一:Fleury算法(在一个Euler图中找出Euler环游)注:包括三个文件;fleuf1.m, edf.m, flecvexf.mfunction [T c]=fleuf1(d)%注:必须保证是Euler环游,否则输出T=0,c=0n=length(d);b=d;b(b==inf)=0;b(b~=0)=1;m=0;a=sum(b);eds=sum(a)/2;ed=zeros(2,eds);vexs=zeros(1,eds+1);matr=b;for i=1:nif mod(a(i),2)==1m=m+1;endendif m~=0fprintf('there is not exit Euler path.\n') T=0;c=0;endif m==0vet=1;flag=0;t1=find(matr(vet,:)==1);for ii=1:length(t1)ed(:,1)=[vet,t1(ii)];vexs(1,1)=vet;vexs(1,2)=t1(ii);matr(vexs(1,2),vexs(1,1))=0;flagg=1;tem=1;while flagg[flagg ed]=edf(matr,eds,vexs,ed,tem); tem=tem+1;if ed(1,eds)~=0 & ed(2,eds)~=0T=ed;T(2,eds)=1;c=0;for g=1:edsc=c+d(T(1,g),T(2,g));endflagg=0;break;endendendendfunction[flag ed]=edf(matr,eds,vexs,ed,tem)flag=1;for i=2:eds[dvex f]=flecvexf(matr,i,vexs,eds,ed,tem);if f==1flag=0;break;endif dvex~=0ed(:,i)=[vexs(1,i) dvex];vexs(1,i+1)=dvex;matr(vexs(1,i+1),vexs(1,i))=0;elsebreak;endendfunction [dvex f]=flecvexf(matr,i,vexs,eds,ed,temp) f=0;edd=find(matr(vexs(1,i),:)==1);dvex=0;dvex1=[];ded=[];if length(edd)==1dvex=edd;elsedd=1;dd1=0;kkk=0;for kk=1:length(edd)m1=find(vexs==edd(kk));if sum(m1)==0dvex1(dd)=edd(kk);dd=dd+1;dd1=1;elsekkk=kkk+1;endendif kkk==length(edd)tem=vexs(1,i)*ones(1,kkk);edd1=[tem;edd];for l1=1:kkklt=0;ddd=1;for l2=1:edsif edd1(1:2,l1)==ed(1:2,l2)lt=lt+1;endendif lt==0ded(ddd)=edd(l1);ddd=ddd+1;endendendif temp<=length(dvex1)dvex=dvex1(temp);elseif temp>length(dvex1) & temp<=length(ded)dvex=ded(temp);elsef=1;endend程序二:Hamilton改良圈算法(找出比较好的Hamilton路,但并不是最优的)function [C d1]= hamiltonglf(v)%d表示权值矩阵%C表示算法最终找到的Hamilton圈。
图论算法及matlab程序的三个案例
图论实验三个案例单源最短路径问题 1.1 Dijkstra 算法Dijkstra 算法是解单源最短路径问题的一个贪心算法。
其基本思想是,设置 一个顶点集合S 并不断地作贪心选择来扩充这个集合。
一个顶点属于集合S 当且 仅当从源到该顶点的最短路径长度已知。
设 v 是图中的一个顶点,记l(v)为顶点 v 到源点V 1的最短距离,V i,V jV ,若(V i,V j)E ,记“到百的权w 。
Dijkstra 算法:① S {V J I(V J 0 ; V V {可 1(V ) i i S V {V J ;J7JJJ7②S,停止,否则转③;l(v) min{ l(v) , d(V j ,v)}V j S④ 存在Vi 1,使l (V i l) min{l(V)},V S ;⑤SSU{v i 1}S S {v i 1}i i 1实际上,Dijkstra 算法也是最优化原理的应用:如果V 1V 2LV n1Vn是从V1到Vn的最短路径,贝UV 1V 2L Vn1也必然是从V1到Vn 1的最优路径。
在下面的MATLA 实现代码中,我们用到了距离矩阵,矩阵第 i 行第j 行元 素表示顶点Vi到Vj的权Wj,若v 到V j无边,则W ijrealmax,其中realmax 是 MATLA 常量,表示最大的实数(1.7977e+308)function re=Dijkstra(ma)%用Dijkstra 算法求单源最短路径%俞入参量ma是距离矩阵%输出参量是一个三行n 列矩阵,每列表示顶点号及顶点到源的最短距离和前顶点n=size(ma,1);% 得到距离矩阵的维数s=ones(1,n);s(1)=0;% 标记集合S和S 的补r=zeros(3,n);r(1,:)=1:n;r(2,2:end)=realmax;% 初始化for i=2:n;% 控制循环次数mm=realmax;for j=find(s==0);% 集合S中的顶点for k=find(s==1);% 集合S补中的顶点if(r(2,j)+ma(j,k)<r(2,k))r(2,k)=r(2,j)+ma(j,k);r(3,k)=j;endif(mm>r(2,k))mm=r(2,k);t=k;endendends(1,t)=0;%找到最小的顶点加入集合Send re=r;1.2动态规划求解最短路径动态规划是美国数学家 Richard Bellman 在1951年提出来的分析一类多阶 段决策过程的最优化方法,在工程技术、工业生产、经济管理、军事及现代化控 制工程等方面均有着广泛的应用。
图论算法
Dijkstra 算法:用矩阵n n a ⨯(n 为顶点个数)存放各边权的邻接矩阵,行向量pb 、1index 、2index、d 分别用来存放P 标号信息、标号顶点顺序、标号顶点索引、最短通路的值。
其中分量⎩⎨⎧=顶点未标号当第顶点已标号当第i i i pb 01)(;)(2i index存放始点到第i 点最短通路中第i 顶点前一顶点的序号;)(i d 存放由始点到第i 点最短通路的值。
求第一个城市到其它城市的最短路径的Matlab 程序如下: clear; clc; M=10000;a(1,:)=[0,50,M,40,25,10]; a(2,:)=[zeros(1,2),15,20,M,25]; a(3,:)=[zeros(1,3),10,20,M]; a(4,:)=[zeros(1,4),10,25]; a(5,:)=[zeros(1,5),55]; a(6,:)=zeros(1,6); a=a+a';pb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a)); d(1:length(a))=M;d(1)=0;temp=1; while sum(pb)<length(a) tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb)); tmpb=find(d(tb)==min(d(tb))); temp=tb(tmpb(1)); pb(temp)=1;index1=[index1,temp];index=index1(find(d(index1)==d(temp)-a(temp,index1))); if length(index)>=2 index=index(1); endindex2(temp)=index; endd, index1, index2%dijkstra 最短路算法通用程序,用于求从起始点s 到其它各点的最短路%D 为赋权邻接矩阵,d 为s 到其它各点最短路径的长度,DD 记载了最短路径生成树 function [d,DD]=dijkstra_aiwa(D,s) [m,n]=size(D); d=inf.*ones(1,m); d(1,s)=0;dd=zeros(1,m);dd(1,s)=1;y=s;DD=zeros(m,m);DD(y,y)=1;counter=1;while length(find(dd==1))<mfor i=1:mif dd(i)==0d(i)=min(d(i),d(y)+D(y,i)); endendddd=inf;for i=1:mif dd(i)==0&&d(i)<dddddd=d(i);endendyy=find(d==ddd);counter=counter+1;DD(y,yy(1,1))=counter;DD(yy(1,1),y)=counter;y=yy(1,1);dd(1,y)=1;endFloyd算法:Matlab程序如下:clear;clc;M=10000;a(1,:)=[0,50,M,40,25,10];a(2,:)=[zeros(1,2),15,20,M,25];a(3,:)=[zeros(1,3),10,20,M];a(4,:)=[zeros(1,4),10,25];a(5,:)=[zeros(1,5),55];a(6,:)=zeros(1,6);b=a+a';path=zeros(length(b));for k=1:6for i=1:6for j=1:6if b(i,j)>b(i,k)+b(k,j)b(i,j)=b(i,k)+b(k,j);path(i,j)=k;end end end end b, pathprim 算法构造最小生成树:prim 算法如下:(i )}{1v P =,Φ=Q ; (ii )while V P =~},,min(P V v P p w pv pv -∈∈= }{v P P += }{pv Q Q += end用prim 算法求右图的最小生成树。
图论算法及Matlab程序代码
图论算法及其MATLAB 程序代码求赋权图G = (V , E , F )中任意两点间的最短路的Warshall-Floyd 算法:设A = (a ij )n ×n 为赋权图G = (V , E , F )的矩阵, 当v i v j ∈E 时a ij = F (v i v j ), 否则取a ii =0, a ij = +∞(i ≠j ), d ij 表示从v i 到v j 点的距离, r ij 表示从v i 到v j 点的最短路中一个点的编号.① 赋初值. 对所有i , j , d ij = a ij , r ij = j . k = 1. 转向②② 更新d ij , r ij . 对所有i , j , 若d ik + d k j <d ij , 则令d ij = d ik + d k j , r ij = k , 转向③.③ 终止判断. 若d ii <0, 则存在一条含有顶点v i 的负回路, 终止; 或者k = n 终止; 否则令k = k + 1, 转向②.最短路线可由r ij 得到.例1 求图6-4中任意两点间的最短路.解:用Warshall-Floyd 算法, MATLAB 程序代码如下:n=8;A=[0 2 8 1 Inf Inf Inf Inf2 0 6 Inf 1 Inf Inf Inf8 6 0 7 5 1 2 Inf1 Inf 7 0 Inf Inf 9 InfInf 1 5 Inf 0 3 Inf 8Inf Inf 1 Inf 3 0 4 6Inf Inf 2 9 Inf 4 0 3Inf Inf Inf Inf 8 6 3 0]; % MATLAB 中, Inf 表示∞D=A; %赋初值for (i=1:n)for (j=1:n)R(i,j)=j;end ;end %赋路径初值for (k=1:n)for (i=1:n)for (j=1:n)if (D(i,k)+D(k,j)<D(i,j))D(i,j)=D(i,k)+D(k,j); %更新dijR(i,j)=k;end ;end ;end %更新rijk %显示迭代步数D %显示每步迭代后的路长R %显示每步迭代后的路径pd=0;for i=1:n %含有负权时if (D(i,i)<0)pd=1;break ;end ;end %存在一条含有顶点vi 的负回路if (pd)break ;end %存在一条负回路, 终止程序end %程序结束图6-4Kruskal避圈法:将图G中的边按权数从小到大逐条考察, 按不构成圈的原则加入到T 中(若有选择时, 不同的选择可能会导致最后生成树的权数不同), 直到q (T ) = p (G ) − 1为止, 即T的边数= G的顶点数− 1为止.Kruskal避圈法的MATLAB程序代码如下:n=8;A=[0 2 8 1 0 0 0 02 0 6 0 1 0 0 08 6 0 7 5 1 2 01 0 7 0 0 0 9 00 1 5 0 0 3 0 80 0 1 0 3 0 4 60 0 2 9 0 4 0 30 0 0 0 8 6 3 0];k=1; %记录A中不同正数的个数for(i=1:n-1)for(j=i+1:n) %此循环是查找A中所有不同的正数if(A(i,j)>0)x(k)=A(i,j); %数组x记录A中不同的正数kk=1; %临时变量for(s=1:k-1)if(x(k)==x(s))kk=0;break;end;end%排除相同的正数k=k+kk;end;end;endk=k-1 %显示A中所有不同正数的个数for(i=1:k-1)for(j=i+1:k) %将x中不同的正数从小到大排序if(x(j)<x(i))xx=x(j);x(j)=x(i);x(i)=xx;end;end;endT(n,n)=0; %将矩阵T中所有的元素赋值为0q=0; %记录加入到树T中的边数for(s=1:k)if(q==n)break;end%获得最小生成树T, 算法终止for(i=1:n-1)for(j=i+1:n)if (A(i,j)==x(s))T(i,j)=x(s);T(j,i)=x(s); %加入边到树T中TT=T; %临时记录Twhile(1)pd=1;%砍掉TT中所有的树枝for(y=1:n)kk=0;for(z=1:n)if(TT(y,z)>0)kk=kk+1;zz=z;end;end%寻找TT中的树枝if(kk==1)TT(y,zz)=0;TT(zz,y)=0;pd=0;end;end%砍掉TT中的树枝if(pd)break;end;end%已砍掉了TT中所有的树枝pd=0;%判断TT中是否有圈for(y=1:n-1)for(z=y+1:n)if(TT(y,z)>0)pd=1;break;end;end;endif(pd)T(i,j)=0;T(j,i)=0;%假如TT中有圈else q=q+1;end;end;end;end;endT %显示近似最小生成树T, 程序结束求二部图G的最大匹配的算法(匈牙利算法), 其基本思想是:从G的任意匹配M开始, 对X中所有M的非饱和点, 寻找M−增广路. 若不存在M−增广路, 则M为最大匹配; 若存在M−增广路P, 则将P中M与非M的边互换得到比M多一边的匹配M1 , 再对M1重复上述过程.设G = ( X, Y, E )为二部图, 其中X = {x1, x2, … , x n }, Y = { y1, y2, … , y n}. 任取G的一初始匹配M (如任取e∈E, 则M = {e}是一个匹配).①令S = φ , T = φ , 转向②.②若M饱和X \S的所有点, 则M是二部图G的最大匹配. 否则, 任取M的非饱和点u∈X \ S , 令S = S ∪{ u }, 转向③.③记N (S ) = {v | u∈S, uv∈E}. 若N (S ) = T, 转向②. 否则取y∈N (S ) \T. 若y是M 的饱和点, 转向④, 否则转向⑤.④设x y∈M, 则令S = S ∪{ x }, T = T ∪{ y }, 转向③.⑤u −y路是M−增广路, 设为P, 并令M = M⊕P, 转向①. 这里M⊕P = M∪P \M∩P, 是对称差.由于计算M−增广路P比较麻烦, 因此将迭代步骤改为:①将X中M的所有非饱和点(不是M中某条边的端点)都给以标号0和标记*, 转向②.②若X中所有有标号的点都已去掉了标记*, 则M是G的最大匹配. 否则任取X中一个既有标号又有标记*的点x i , 去掉x i的标记*, 转向③.③找出在G中所有与x i邻接的点y j (即x i y j∈E ), 若所有这样的y j都已有标号, 则转向②, 否则转向④.④对与x i邻接且尚未给标号的y j都给定标号i. 若所有的y j都是M的饱和点, 则转向⑤, 否则逆向返回. 即由其中M的任一个非饱和点y j的标号i找到x i, 再由x i的标号k找到y k , … , 最后由y t的标号s找到标号为0的x s时结束, 获得M−增广路x s y t…x i y j, 记P = {x s y t, …, x i y j }, 重新记M为M⊕P, 转向①.⑤将y j在M中与之邻接的点x k (即x k y j∈M), 给以标号j和标记*, 转向②.例1求图6-9中所示的二部图G的最大匹配.图6-9匈牙利算法的MATLAB程序代码如下:m=5;n=5;A=[0 1 1 0 01 1 0 1 10 1 1 0 00 1 1 0 00 0 0 1 1];M(m,n)=0;for(i=1:m)for(j=1:n)if(A(i,j))M(i,j)=1;break;end;end%求初始匹配Mif(M(i,j))break;end;end%获得仅含一条边的初始匹配Mwhile(1)for(i=1:m)x(i)=0;end%将记录X中点的标号和标记*for(i=1:n)y(i)=0;end%将记录Y中点的标号和标记*for(i=1:m)pd=1;%寻找X中M的所有非饱和点for(j=1:n)if(M(i,j))pd=0;end;endif(pd)x(i)=-n-1;end;end%将X中M的所有非饱和点都给以标号0和标记*, 程序中用n+1表示0标号, 标号为负数时表示标记*pd=0;while(1)xi=0;for(i=1:m)if(x(i)<0)xi=i;break;end;end%假如X中存在一个既有标号又有标记*的点, 则任取X中一个既有标号又有标记*的点xiif(xi==0)pd=1;break;end%假如X中所有有标号的点都已去掉了标记*, 算法终止x(xi)=x(xi)*(-1); %去掉xi的标记*k=1;for(j=1:n)if(A(xi,j)&y(j)==0)y(j)=xi;yy(k)=j;k=k+1;end;end%对与xi邻接且尚未给标号的yj都给以标号iif(k>1)k=k-1;for(j=1:k)pdd=1;for(i=1:m)if(M(i,yy(j)))x(i)=-yy(j);pdd=0;break;end;end%将yj在M中与之邻接的点xk (即xkyj∈M), 给以标号j和标记*if(pdd)break;end;endif(pdd)k=1;j=yy(j); %yj不是M的饱和点while(1)P(k,2)=j;P(k,1)=y(j);j=abs(x(y(j))); %任取M的一个非饱和点yj, 逆向返回if(j==n+1)break;end%找到X中标号为0的点时结束, 获得M-增广路Pk=k+1;endfor(i=1:k)if(M(P(i,1),P(i,2)))M(P(i,1),P(i,2))=0; %将匹配M在增广路P中出现的边去掉else M(P(i,1),P(i,2))=1;end;end%将增广路P中没有在匹配M中出现的边加入到匹配M中break;end;end;endif(pd)break;end;end%假如X中所有有标号的点都已去掉了标记*, 算法终止M %显示最大匹配M, 程序结束利用可行点标记求最佳匹配的算法步骤如下:设G = ( X , Y , E , F )为完备的二部赋权图, L 是其一个初始可行点标记, 通常取.,,0)(},|)(max{)(Y y X x y L Y y xy F x L ∈∈ =∈= M 是G L 的一个匹配. ① 若X 的每个点都是M 的饱和点, 则M 是最佳匹配. 否则取M 的非饱和点u ∈X , 令S = {u }, T = φ , 转向②.② 记N L (S ) = {v | u ∈S , uv ∈E L }. 若N L ( S ) = T , 则G L 没有完美匹配, 转向③. 否则转向④.③ 调整可行点标记, 计算a L = min { L ( x ) + L ( y ) − F (x y ) | x ∈S , y ∈Y \T }.由此得新的可行顶点标记H (v ) =,,),(,)(,)(T v S v v L a v L a v L L L ∈∈+−令L = H , G L = G H , 重新给出G L 的一个匹配M , 转向①.④ 取y ∈N L ( S ) \T , 若y 是M 的饱和点, 转向⑤. 否则, 转向⑥.⑤ 设x y ∈M , 则令S = S ∪{ x }, T = T ∪{ y }, 转向②.⑥ 在G L 中的u − y 路是M −增广路, 记为P , 并令 M = M ⊕P , 转向①.利用可行点标记求最佳匹配算法的MATLAB 程序代码如下:n=4;A=[4 5 5 12 2 4 64 2 3 35 0 2 1];for (i=1:n)L(i,1)=0;L(i,2)=0;endfor (i=1:n)for (j=1:n)if (L(i,1)<A(i,j))L(i,1)=A(i,j);end ; %初始可行点标记LM(i,j)=0;end ;endfor (i=1:n)for (j=1:n) %生成子图Glif (L(i,1)+L(j,2)==A(i,j))Gl(i,j)=1;else Gl(i,j)=0;end ;end ;endii=0;jj=0;for (i=1:n)for (j=1:n)if (Gl(i,j))ii=i;jj=j;break ;end ;endif (ii)break ;end ;end %获得仅含Gl 的一条边的初始匹配MM(ii,jj)=1;for (i=1:n)S(i)=0;T(i)=0;NlS(i)=0;endwhile (1)for (i=1:n)k=1;否则.for(j=1:n)if(M(i,j))k=0;break;end;endif(k)break;end;endif(k==0)break;end%获得最佳匹配M, 算法终止S(1)=i;jss=1;jst=0;%S={xi}, T=φwhile(1)jsn=0;for(i=1:jss)for(j=1:n)if(Gl(S(i),j))jsn=jsn+1;NlS(jsn)=j;%NL(S)={v|u∈S,uv∈EL}for(k=1:jsn-1)if(NlS(k)==j)jsn=jsn-1;end;end;end;end;endif(jsn==jst)pd=1; %判断NL(S)=T?for(j=1:jsn)if(NlS(j)~=T(j))pd=0;break;end;end;endif(jsn==jst&pd)al=Inf; %如果NL(S)=T, 计算al, Inf为∞for(i=1:jss)for(j=1:n)pd=1;for(k=1:jst)if(T(k)==j)pd=0;break;end;endif(pd&al>L(S(i),1)+L(j,2)-A(S(i),j))al=L(S(i),1)+L(j,2)-A(S(i),j);end;end;end for(i=1:jss)L(S(i),1)=L(S(i),1)-al;end%调整可行点标记for(j=1:jst)L(T(j),2)=L(T(j),2)+al;end%调整可行点标记for(i=1:n)for(j=1:n) %生成子图GLif(L(i,1)+L(j,2)==A(i,j))Gl(i,j)=1;else Gl(i,j)=0;endM(i,j)=0;k=0;end;endii=0;jj=0;for(i=1:n)for(j=1:n)if(Gl(i,j))ii=i;jj=j;break;end;endif(ii)break;end;end%获得仅含Gl的一条边的初始匹配MM(ii,jj)=1;breakelse%NL(S)≠Tfor(j=1:jsn)pd=1;%取y∈NL(S)\Tfor(k=1:jst)if(T(k)==NlS(j))pd=0;break;end;endif(pd)jj=j;break;end;endpd=0;%判断y是否为M的饱和点for(i=1:n)if(M(i,NlS(jj)))pd=1;ii=i;break;end;endif(pd)jss=jss+1;S(jss)=ii;jst=jst+1;T(jst)=NlS(jj); %S=S∪{x}, T=T∪{y}else%获得Gl的一条M-增广路, 调整匹配Mfor(k=1:jst)M(S(k),T(k))=1;M(S(k+1),T(k))=0;endif(jst==0)k=0;endM(S(k+1),NlS(jj))=1;break;end;end;end;endMaxZjpp=0;for(i=1:n)for(j=1:n)if(M(i,j))MaxZjpp=MaxZjpp+A(i,j);end;end;endM %显示最佳匹配MMaxZjpp %显示最佳匹配M的权, 程序结束从一个可行流f 开始, 求最大流的Ford--Fulkerson 标号算法的基本步骤:⑴ 标号过程① 给发点v s 以标号(+, +∞) , δ s = +∞.② 选择一个已标号的点x , 对于x 的所有未给标号的邻接点y , 按下列规则处理:当yx ∈E , 且f yx >0时, 令δ y = min { f yx , δ x }, 并给y 以标号 ( x − , δ y ).当xy ∈E , 且f xy <C xy 时, 令δ y = min {C xy − f xy , δ x }, 并给y 以标号 ( x + , δ y ). ③ 重复②直到收点v t 被标号或不再有点可标号时为止. 若v t 得到标号, 说明存在一条可增广链, 转⑵调整过程; 若v t 未得到标号, 标号过程已无法进行时, 说明f 已经是最大流.⑵ 调整过程④ 决定调整量δ =δ vt , 令u = v t .⑤ 若u 点标号为( v +, δ u ), 则以f vu + δ 代替f vu ; 若u 点标号为( v −, δ u ), 则以 f vu − δ 代替f vu .⑥ 若v = v s , 则去掉所有标号转⑴重新标号; 否则令u = v , 转⑤.算法终止后, 令已有标号的点集为S , 则割集(S , S c )为最小割, 从而W f = C (S , S c ). 例1 求图6-19所示网络的最大流.利用Ford--Fulkerson 标号法求最大流算法的MATLAB 程序代码如下:n=8;C=[0 5 4 3 0 0 0 00 0 0 0 5 3 0 00 0 0 0 0 3 2 00 0 0 0 0 0 2 00 0 0 0 0 0 0 40 0 0 0 0 0 0 30 0 0 0 0 0 0 50 0 0 0 0 0 0 0]; %弧容量for (i=1:n)for (j=1:n)f(i,j)=0;end ;end %取初始可行流f 为零流for (i=1:n)No(i)=0;d(i)=0;end %No,d 记录标号图6-19while(1)No(1)=n+1;d(1)=Inf; %给发点vs标号while(1)pd=1;%标号过程for(i=1:n)if(No(i)) %选择一个已标号的点vifor(j=1:n)if(No(j)==0&f(i,j)<C(i,j)) %对于未给标号的点vj, 当vivj为非饱和弧时No(j)=i;d(j)=C(i,j)-f(i,j);pd=0;if(d(j)>d(i))d(j)=d(i);endelseif(No(j)==0&f(j,i)>0) %对于未给标号的点vj, 当vjvi为非零流弧时No(j)=-i;d(j)=f(j,i);pd=0;if(d(j)>d(i))d(j)=d(i);end;end;end;end;endif(No(n)|pd)break;end;end%若收点vt得到标号或者无法标号, 终止标号过程if(pd)break;end%vt未得到标号, f已是最大流, 算法终止dvt=d(n);t=n; %进入调整过程, dvt表示调整量while(1)if(No(t)>0)f(No(t),t)=f(No(t),t)+dvt; %前向弧调整elseif(No(t)<0)f(No(t),t)=f(No(t),t)-dvt;end%后向弧调整if(No(t)==1)for(i=1:n)No(i)=0;d(i)=0; end;break;end%当t的标号为vs时, 终止调整过程t=No(t);end;end; %继续调整前一段弧上的流fwf=0;for(j=1:n)wf=wf+f(1,j);end%计算最大流量f %显示最大流wf %显示最大流量No %显示标号, 由此可得最小割, 程序结束设网络G = ( V , E , C ), 取初始可行流 f 为零流, 求解最小费用流问题的迭代步骤: ① 构造有向赋权图 G f = ( V , E f , F ), 对于任意的v i v j ∈E , E f , F 的定义如下:当f ij = 0时, v i v j ∈E f , F ( v i v j ) = b ij ;当f ij = C ij 时, v j v i ∈E f , F ( v j v i ) = −b ij ;当0< f ij <C ij 时, v i v j ∈E f , F ( v i v j ) = b ij , v j v i ∈E f , F ( v j v i ) = −b ij .转向②.② 求出有向赋权图G f = (V , E f , F )中发点v s 到收点v t 的最短路µ , 若最短路µ存在转向③; 否则f 是所求的最小费用最大流, 停止.③ 增流. 同求最大流的方法一样, 重述如下:令.,,,−+∈∈ −=µµδj i j i ij ij ij ij v v v v f f C δ = min {δ ij | v i v j ∈µ}, 重新定义流f = { f ij }为 f ij =,,,,−+∈∈ −+µµδδj i j i ijij ij v v v v f f f如果W f 大于或等于预定的流量值, 则适当减少δ 值, 使W f 等于预定的流量值, 那么 f 是所求的最小费用流, 停止; 否则转向①.求解含有负权的有向赋权图G = ( V , E , F )中某一点到其它各点最短路的Ford 算法. 当v i v j ∈E 时记w ij = F (v i v j ), 否则取w ii =0, w ij = +∞(i ≠j ). v 1到v i 的最短路长记为π ( i ), v 1到v i 的最短路中v i 的前一个点记为θ ( i ). Ford 算法的迭代步骤:① 赋初值π (1) = 0, π ( i ) = +∞, θ ( i ) = i , i = 2, 3, … , n .② 更新π ( i ), θ ( i ). 对于i = 2, 3, … , n 和j = 1, 2, … , n , 如果π ( i )<π ( j ) + w ji , 则令π ( i ) = π ( j ) , θ ( i ) = j . ③ 终止判断:若所有的π ( i )都无变化, 停止; 否则转向②. 在算法的每一步中, π ( i )都是从v 1到v i 的最短路长度的上界. 若不存在负长回路, 则从v 1到v i 的最短路长度是π ( i )的下界, 经过n −1次迭代后π ( i )将保持不变. 若在第n 次迭代后π ( i )仍在变化时, 说明存在负长回路.其它.例2 在图6-22所示运输网络上, 求s 到t 的最小费用最大流, 括号内为(C ij , b ij ).求最小费用最大流算法的MATLAB 程序代码如下:n=5;C=[0 15 16 0 00 0 0 13 140 11 0 17 00 0 0 0 80 0 0 0 0]; %弧容量b=[0 4 1 0 00 0 0 6 10 2 0 3 00 0 0 0 20 0 0 0 0]; %弧上单位流量的费用wf=0;wf0=Inf; %wf 表示最大流量, wf0表示预定的流量值for (i=1:n)for (j=1:n)f(i,j)=0;end ;end %取初始可行流f 为零流while (1)for (i=1:n)for (j=1:n)if (j~=i)a(i,j)=Inf;end ;end ;end %构造有向赋权图for (i=1:n)for (j=1:n)if (C(i,j)>0&f(i,j)==0)a(i,j)=b(i,j);elseif (C(i,j)>0&f(i,j)==C(i,j))a(j,i)=-b(i,j);elseif (C(i,j)>0)a(i,j)=b(i,j);a(j,i)=-b(i,j);end ;end ;endfor (i=2:n)p(i)=Inf;s(i)=i;end %用Ford 算法求最短路, 赋初值for (k=1:n)pd=1; %求有向赋权图中vs 到vt 的最短路for (i=2:n)for (j=1:n)if (p(i)>p(j)+a(j,i))p(i)=p(j)+a(j,i);s(i)=j;pd=0;end ;end ;endif (pd)break ;end ;end %求最短路的Ford 算法结束if (p(n)==Inf)break ;end %不存在vs 到vt 的最短路, 算法终止. 注意在求最小费用最大流时构造有向赋权图中不会含负权回路, 所以不会出现k=ndvt=Inf;t=n; %进入调整过程, dvt 表示调整量while (1) %计算调整量if (a(s(t),t)>0)dvtt=C(s(t),t)-f(s(t),t); %前向弧调整量elseif (a(s(t),t)<0)dvtt=f(t,s(t));end %后向弧调整量if (dvt>dvtt)dvt=dvtt;endif (s(t)==1)break ;end %当t 的标号为vs 时, 终止计算调整量t=s(t);end %继续调整前一段弧上的流fpd=0;if (wf+dvt>=wf0)dvt=wf0-wf;pd=1;end %如果最大流量大于或等于预定的流量值t=n;while (1) %调整过程if (a(s(t),t)>0)f(s(t),t)=f(s(t),t)+dvt; %前向弧调整elseif (a(s(t),t)<0)f(t,s(t))=f(t,s(t))-dvt;end %后向弧调整if (s(t)==1)break ;end %当t 的标号为vs 时, 终止调整过程t=s(t);endif (pd)break ;end %如果最大流量达到预定的流量值wf=0; for (j=1:n)wf=wf+f(1,j);end ;end %计算最大流量zwf=0;for (i=1:n)for (j=1:n)zwf=zwf+b(i,j)*f(i,j);end ;end %计算最小费用f %显示最小费用最大流图6-22wf %显示最小费用最大流量zwf %显示最小费用, 程序结束。
超全图论matlab程序
超全的图论程序关注微信公众号“超级数学建模”,教你做有料、有趣的数模人程序一:可达矩阵算法function P=dgraf(A)n=size(A,1);P=A;for i=2:nP=P+A^i;endP(P~=0)=1;P;程序二:关联矩阵和邻接矩阵互换算法function W=incandadf(F,f)if f==0m=sum(sum(F))/2;n=size(F,1);W=zeros(n,m);k=1;for i=1:nfor j=i:nif F(i,j)~=0W(i,k)=1;W(j,k)=1;k=k+1;endendendelseif f==1m=size(F,2);n=size(F,1);W=zeros(n,n);for i=1:ma=find(F(:,i)~=0);W(a(1),a(2))=1;W(a(2),a(1))=1;endelsefprint('Please imput the right value of f');endW;程序三:有向图关联矩阵和邻接矩阵互换算法function W=mattransf(F,f)if f==0m=sum(sum(F));n=size(F,1);W=zeros(n,m);k=1;for i=1:nfor j=i:nif F(i,j)~=0W(i,k)=1;W(j,k)=-1;k=k+1;endendendelseif f==1m=size(F,2);n=size(F,1);W=zeros(n,n);for i=1:ma=find(F(:,i)~=0);if F(a(1),i)==1W(a(1),a(2))=1;elseW(a(2),a(1))=1;endendelsefprint('Please imput the right value of f'); endW;第二讲:最短路问题程序一:Dijkstra算法(计算两点间的最短路)function [l,z]=Dijkstra(W)n = size (W,1);for i = 1 :nl(i)=W(1,i);z(i)=0;endi=1;while i<=nfor j =1 :nif l(i)>l(j)+W(j,i)l(i)=l(j)+W(j,i);z(i)=j-1;if j<ii=j-1;endendendi=i+1;end程序二:floyd算法(计算任意两点间的最短距离)function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendr;for k=1:nfor i=1:nfor j=1:nif d(i,k)+d(k,j)<d(i,j)d(i,j)=d(i,k)+d(k,j); r(i,j)=r(i,k);endendendend程序三:n2short.m 计算指定两点间的最短距离function [P u]=n2short(W,k1,k2)n=length(W);U=W;m=1;while m<=nfor i=1:nfor j=1:nif U(i,j)>U(i,m)+U(m,j)U(i,j)=U(i,m)+U(m,j);endendendm=m+1;endu=U(k1,k2);P1=zeros(1,n);k=1;P1(k)=k2;V=ones(1,n)*inf;kk=k2;while kk~=k1for i=1:nV(1,i)=U(k1,kk)-W(i,kk);if V(1,i)==U(k1,i)P1(k+1)=i;kk=i;k=k+1;endendendk=1;wrow=find(P1~=0);for j=length(wrow):-1:1P(k)=P1(wrow(j));k=k+1;endP;程序四、n1short.m(计算某点到其它所有点的最短距离)function[Pm D]=n1short(W,k)n=size(W,1);D=zeros(1,n);for i=1:n[P d]=n2short(W,k,i);Pm{i}=P;D(i)=d;end程序五:pass2short.m(计算经过某两点的最短距离) function [P d]=pass2short(W,k1,k2,t1,t2)[p1 d1]=n2short(W,k1,t1);[p2 d2]=n2short(W,t1,t2);[p3 d3]=n2short(W,t2,k2);dt1=d1+d2+d3;[p4 d4]=n2short(W,k1,t2);[p5 d5]=n2short(W,t2,t1);[p6 d6]=n2short(W,t1,k2);dt2=d4+d5+d6;if dt1<dt2d=dt1;P=[p1 p2(2:length(p2)) p3(2:length(p3))]; elsed=dt1;p=[p4 p5(2:length(p5)) p6(2:length(p6))]; endP;d;第三讲:最小生成树程序一:最小生成树的Kruskal算法function [T c]=krusf(d,flag)if nargin==1n=size(d,2);m=sum(sum(d~=0))/2;b=zeros(3,m);k=1;for i=1:nfor j=(i+1):nif d(i,j)~=0b(1,k)=i;b(2,k)=j;b(3,k)=d(i,j);k=k+1;endendendelseb=d;endn=max(max(b(1:2,:)));m=size(b,2);[B,i]=sortrows(b',3);B=B';c=0;T=[];k=1;t=1:n;for i=1:mif t(B(1,i))~=t(B(2,i))T(1:2,k)=B(1:2,i);c=c+B(3,i);k=k+1;tmin=min(t(B(1,i)),t(B(2,i)));tmax=max(t(B(1,i)),t(B(2,i)));for j=1:nif t(j)==tmaxt(j)=tmin;endendendif k==nbreak;endendT;c;程序二:最小生成树的Prim算法function [T c]=Primf(a)l=length(a);a(a==0)=inf;k=1:l;listV(k)=0;listV(1)=1;e=1;while (e<l)min=inf;for i=1:lif listV(i)==1for j=1:lif listV(j)==0 & min>a(i,j)min=a(i,j);b=a(i,j);s=i;d=j;endendendendlistV(d)=1;distance(e)=b;source(e)=s;destination(e)=d;e=e+1;endT=[source;destination];for g=1:e-1c(g)=a(T(1,g),T(2,g));endc;另外两种程序最小生成树程序1(prim 算法构造最小生成树)a=[inf 50 60 inf inf inf inf;50 inf inf 65 40 inf inf;60 inf inf 52 inf inf 45;...inf 65 52 inf 50 30 42;inf 40 inf 50 inf 70 inf;inf inf inf 30 70 inf inf;...inf inf 45 42 inf inf inf];result=[];p=1;tb=2:length(a);while length(result)~=length(a)-1temp=a(p,tb);temp=temp(:);d=min(temp);[jb,kb]=find(a(p,tb)==d);j=p(jb(1));k=tb(kb(1));result=[result,[j;k;d]];p=[p,k];tb(find(tb==k))=[];endresult最小生成树程序2(Kruskal 算法构造最小生成树)clc;clear;a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40;a(3,4)=52;a(3,7)=45; a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70;[i,j,b]=find(a);data=[i';j';b'];index=data(1:2,:);loop=max(size(a))-1;result=[];while length(result)<looptemp=min(data(3,:));flag=find(data(3,:)==temp);flag=flag(1);v1=data(1,flag);v2=data(2,flag);if index(1,flag)~=index(2,flag)result=[result,data(:,flag)];endindex(find(index==v2))=v1;data(:,flag)=[];index(:,flag)=[];endresult第四讲:Euler图和Hamilton图程序一:Fleury算法(在一个Euler图中找出Euler环游)注:包括三个文件;fleuf1.m, edf.m, flecvexf.m function [T c]=fleuf1(d)%注:必须保证是Euler环游,否则输出T=0,c=0n=length(d);b=d;b(b==inf)=0;b(b~=0)=1;m=0;a=sum(b);eds=sum(a)/2;ed=zeros(2,eds);vexs=zeros(1,eds+1);matr=b;for i=1:nif mod(a(i),2)==1m=m+1;endendif m~=0fprintf('there is not exit Euler path.\n') T=0;c=0;endif m==0vet=1;flag=0;t1=find(matr(vet,:)==1);for ii=1:length(t1)ed(:,1)=[vet,t1(ii)];vexs(1,1)=vet;vexs(1,2)=t1(ii);matr(vexs(1,2),vexs(1,1))=0;flagg=1;tem=1;while flagg[flagg ed]=edf(matr,eds,vexs,ed,tem); tem=tem+1;if ed(1,eds)~=0 & ed(2,eds)~=0T=ed;T(2,eds)=1;c=0;for g=1:edsc=c+d(T(1,g),T(2,g));endflagg=0;break;endendendendfunction[flag ed]=edf(matr,eds,vexs,ed,tem)flag=1;for i=2:eds[dvex f]=flecvexf(matr,i,vexs,eds,ed,tem);if f==1flag=0;break;endif dvex~=0ed(:,i)=[vexs(1,i) dvex];vexs(1,i+1)=dvex;matr(vexs(1,i+1),vexs(1,i))=0;elsebreak;endendfunction [dvex f]=flecvexf(matr,i,vexs,eds,ed,temp) f=0;edd=find(matr(vexs(1,i),:)==1);dvex=0;dvex1=[];ded=[];if length(edd)==1dvex=edd;elsedd=1;dd1=0;kkk=0;for kk=1:length(edd)m1=find(vexs==edd(kk));if sum(m1)==0dvex1(dd)=edd(kk);dd=dd+1;dd1=1;elsekkk=kkk+1;endendif kkk==length(edd)tem=vexs(1,i)*ones(1,kkk);edd1=[tem;edd];for l1=1:kkklt=0;ddd=1;for l2=1:edsif edd1(1:2,l1)==ed(1:2,l2)lt=lt+1;endendif lt==0ded(ddd)=edd(l1);ddd=ddd+1;endendendif temp<=length(dvex1)dvex=dvex1(temp);elseif temp>length(dvex1) & temp<=length(ded)dvex=ded(temp);elsef=1;endend程序二:Hamilton改良圈算法(找出比较好的Hamilton路)function [C d1]= hamiltonglf(v)%d表示权值矩阵%C表示算法最终找到的Hamilton圈。
应用篇-第14章-图论算法及其MATLAB实现
对每个图G=(V,E),均有
d(v)=2|E|
vV
证明:根据顶点度的定义,在计算点度时每条边对于它所关联的顶点被计 算了两次。因此,图G 中点度的总和恰为边数|V|的2倍。证毕。 推论14.1 在任何图G=(V,E)中,奇点的个数为偶数。
定理14.4对任意有向图D=(V,A)均有
14.8 Dijkstra 算法及其MATLAB实现
14.8.1 问题描述与算法思想
Dijkstra算法是解单源最短路径问题的一个贪心算法。其基本思想是,设置 一个顶点集合S并不断地作贪心选择来扩充这个集合。一个顶点属于集合S当 且仅当从源到该顶点的最短路径长度已知。设v是图中的一个顶点,记L(v) 为顶点v到源点v1的最短距离, vi , v j V 若 (vi , v j ) E ,记vi到vj的权。
① A(G)为对称矩阵;
② 若G 为无环图,则A(G)中第i行(列)的元素之和等于顶点vi的度; ③ 两图G 和H 同构的充分必要条件是存在置换矩阵P 使得A(G)=PTA(H)P。
类似地,有向图D 的邻接矩阵A(D)=(aij)n×n的元素aij定义为:元素aij表示从始点vi到 终点vj的有向边的条,,其中vi和vj为D 的顶点。
14.3.2 关联矩=(V,E),其中顶点集V={v1,v2,…,vn},边集E={e1,e2,…,eε}。 用 mij表示顶点vi与边ej关联的次数,可能取值为0,1,2,称所得矩阵M(G)=(mij)n×ε为 图G 的关联矩阵。 类似地,有向图D 的关联矩阵M(D)=(mij)n×ε的元素mij定义为:
14.2.4 路
在图论理论中,路具有特殊的重要性,古往今来,许多学者均对它进行过深入研究。本 节主要介绍简单图G=(V,E)中有关路和连通性的简单性质。 定理14.1 若图G 中有一条(u,v)途径,则G 中也存在一条(u,v)路。
图论和网络分析算法及Matlab实现(Graph_and_Network_Analysis)
2017/11/5
问题的两个共同特点
(1)目的都是从若干可能的安排或方案中寻求 某种意义下的最优安排或方案,数学问题称 为最优化或优化问题。 (2)它们都可用图形形式直观描述,数学上把这 种与图相关的结构称为网络。图和网络相关 的最优化问题就是网络最优化。 网络优化问题是以网络流为研究的对象,常 常被称为网络流或网络流规划等。
v2
2
v1 3 5 1
v3
2
7 5 3 5 v5
v6 1 7
5
v7
v4
• 2. 方法:Dijkstra算法(Dijkstra,1959)
Dijkstra, E.W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik 1, 269–271.
2017/11/5
5 、旅行商问题 Traveling salesman problem
一名推销员准备前往若干城市推销产 品。如何为他设计一条最短的旅行 路线? (从驻地出发,经过每个城 市恰好一次,最后返回驻地)
2017/11/5
6、运输问题 Transportation problem
某种原材料有 M个产地,现在需要将原材料从产 地运往 N个使用这些原材料的工厂。假定 M个产 地的产量和 N家工厂的需要量已知,单位产品从 任一产地到任一工厂的运费已知,那么如何安排 运输方案可以使总运输成本最低?
wij,i能一步到达j d ij j ,i不能一步到达
2017/11/5
Dijkstra 算法
由图G建立一步可达距离阵D=(dij)n×n
给V1(Vs)括号(l1,Vk)=(0,s)给出已标号集合 I和未标号集合J的元素
(图论)matlab模板程序
第一讲:图论模型程序一:可达矩阵算法%根据邻接矩阵A〔有向图〕求可达矩阵P〔有向图〕function P=dgraf<A>n=size<A,1>;P=A;for i=2:nP=P+A^i;endP<P~=0>=1; %将不为0的元素变为1P;程序二:无向图关联矩阵和邻接矩阵互换算法F表示所给出的图的相应矩阵W表示程序运行结束后的结果f=0表示把邻接矩阵转换为关联矩阵f=1表示把关联矩阵转换为邻接矩阵%无向图的关联矩阵和邻接矩阵的相互转换function W=incandadf<F,f>if f==0 %邻接矩阵转换为关联矩阵m=sum<sum<F>>/2; %计算图的边数n=size<F,1>;W=zeros<n,m>;k=1;for i=1:nfor j=i:nif F<i,j>~=0W<i,k>=1; %给边的始点赋值为1W<j,k>=1; %给边的终点赋值为1k=k+1;endendendelseif f==1 %关联矩阵转换为邻接矩阵m=size<F,2>;n=size<F,1>;W=zeros<n,n>;for i=1:ma=find<F<:,i>~=0>;W<a<1>,a<2>>=1; %存在边,则邻接矩阵的对应值为1 W<a<2>,a<1>>=1;endelsefprint<'Please imput the right value of f'>;W;程序三:有向图关联矩阵和邻接矩阵互换算法%有向图的关联矩阵和邻接矩阵的转换function W=mattransf<F,f>if f==0 %邻接矩阵转换为关联矩阵m=sum<sum<F>>;n=size<F,1>;W=zeros<n,m>;k=1;for i=1:nfor j=i:nif F<i,j>~=0 %由i发出的边,有向边的始点W<i,k>=1; %关联矩阵始点值为1W<j,k>=-1; %关联矩阵终点值为-1k=k+1;endendendelseif f==1 %关联矩阵转换为邻接矩阵m=size<F,2>;n=size<F,1>;W=zeros<n,n>;for i=1:ma=find<F<:,i>~=0>; %有向边的两个顶点if F<a<1>,i>==1W<a<1>,a<2>>=1; %有向边由a<1>指向a<2>elseW<a<2>,a<1>>=1; %有向边由a<2>指向a<1>endendelsefprint<'Please imput the right value of f'>;endW;第二讲:最短路问题程序0:最短距离矩阵W表示图的权值矩阵D表示图的最短距离矩阵%连通图中各项顶点间最短距离的计算function D=shortdf<W>%对于W<i,j>,若两顶点间存在弧,则为弧的权值,否则为inf;当i=j时W<i,j>=0 n=length<W>;m=1;while m<=nfor i=1:nfor j=1:nif D<i,j>>D<i,m>+D<m,j>D<i,j>+D<i,m>+D<m,j>; %距离进行更新 endendendm=m+1;endD;程序一:Dijkstra算法〔计算两点间的最短路〕function [l,z]=Dijkstra<W>n = size <W,1>;for i = 1 :nl<i>=W<1,i>;z<i>=0;endi=1;while i<=nfor j =1 :nif l<i>>l<j>+W<j,i>l<i>=l<j>+W<j,i>;z<i>=j-1;if j<ii=j-1;endendendi=i+1;end程序二:floyd算法〔计算任意两点间的最短距离〕function [d,r]=floyd<a>n=size<a,1>;d=a;for i=1:nfor j=1:nr<i,j>=j;endendr;for k=1:nfor i=1:nfor j=1:nif d<i,k>+d<k,j><d<i,j>d<i,j>=d<i,k>+d<k,j>; r<i,j>=r<i,k>;endendendend程序三:n2short.m 计算指定两点间的最短距离function [P u]=n2short<W,k1,k2>n=length<W>;U=W;m=1;while m<=nfor i=1:nfor j=1:nif U<i,j>>U<i,m>+U<m,j>U<i,j>=U<i,m>+U<m,j>;endendendm=m+1;endu=U<k1,k2>;P1=zeros<1,n>;k=1;P1<k>=k2;V=ones<1,n>*inf;kk=k2;while kk~=k1for i=1:nV<1,i>=U<k1,kk>-W<i,kk>;if V<1,i>==U<k1,i>P1<k+1>=i;kk=i;k=k+1;endendendk=1;wrow=find<P1~=0>;for j=length<wrow>:-1:1P<k>=P1<wrow<j>>;k=k+1;endP;程序四、n1short.m<计算某点到其它所有点的最短距离> function[Pm D]=n1short<W,k>n=size<W,1>;D=zeros<1,n>;for i=1:n[P d]=n2short<W,k,i>;Pm{i}=P;D<i>=d;end程序五:pass2short.m<计算经过某两点的最短距离> function [P d]=pass2short<W,k1,k2,t1,t2>[p1 d1]=n2short<W,k1,t1>;[p2 d2]=n2short<W,t1,t2>;[p3 d3]=n2short<W,t2,k2>;dt1=d1+d2+d3;[p4 d4]=n2short<W,k1,t2>;[p5 d5]=n2short<W,t2,t1>;[p6 d6]=n2short<W,t1,k2>;dt2=d4+d5+d6;if dt1<dt2d=dt1;P=[p1 p2<2:length<p2>> p3<2:length<p3>>]; elsed=dt1;p=[p4 p5<2:length<p5>> p6<2:length<p6>>]; endP;d;第三讲:最小生成树程序一:最小生成树的Kruskal算法function [T c]=krusf<d,flag>if nargin==1n=size<d,2>;m=sum<sum<d~=0>>/2;b=zeros<3,m>;k=1;for i=1:nfor j=<i+1>:nif d<i,j>~=0b<1,k>=i;b<2,k>=j;b<3,k>=d<i,j>;k=k+1;endendendelseb=d;endn=max<max<b<1:2,:>>>;m=size<b,2>;[B,i]=sortrows<b',3>;B=B';c=0;T=[];k=1;t=1:n;for i=1:mif t<B<1,i>>~=t<B<2,i>>T<1:2,k>=B<1:2,i>;c=c+B<3,i>;k=k+1;tmin=min<t<B<1,i>>,t<B<2,i>>>; tmax=max<t<B<1,i>>,t<B<2,i>>>; for j=1:nif t<j>==tmaxt<j>=tmin;endendendif k==nbreak;endendT;c;程序二:最小生成树的Prim算法function [T c]=Primf<a>l=length<a>;a<a==0>=inf;k=1:l;listV<k>=0;listV<1>=1;e=1;while <e<l>min=inf;for i=1:lif listV<i>==1for j=1:lif listV<j>==0 & min>a<i,j>min=a<i,j>;b=a<i,j>;s=i;d=j;endendendendlistV<d>=1;distance<e>=b;source<e>=s;destination<e>=d;e=e+1;endT=[source;destination];for g=1:e-1c<g>=a<T<1,g>,T<2,g>>;endc;第四讲:Euler图和Hamilton图程序一:Fleury算法〔在一个Euler图中找出Euler环游〕注:包括三个文件;fleuf1.m, edf.m, flecvexf.mfunction [T c]=fleuf1<d>%注:必须保证是Euler环游,否则输出T=0,c=0n=length<d>;b=d;b<b==inf>=0;b<b~=0>=1;m=0;a=sum<b>;eds=sum<a>/2;ed=zeros<2,eds>;vexs=zeros<1,eds+1>;matr=b;for i=1:nif mod<a<i>,2>==1m=m+1;endendif m~=0fprintf<'there is not exit Euler path.\n'>T=0;c=0;endif m==0vet=1;flag=0;t1=find<matr<vet,:>==1>;for ii=1:length<t1>ed<:,1>=[vet,t1<ii>];vexs<1,1>=vet;vexs<1,2>=t1<ii>;matr<vexs<1,2>,vexs<1,1>>=0;flagg=1;tem=1;while flagg[flagg ed]=edf<matr,eds,vexs,ed,tem>;tem=tem+1;if ed<1,eds>~=0 & ed<2,eds>~=0T=ed;T<2,eds>=1;c=0;for g=1:edsc=c+d<T<1,g>,T<2,g>>;endflagg=0;break;endendendendfunction[flag ed]=edf<matr,eds,vexs,ed,tem>flag=1;for i=2:eds[dvex f]=flecvexf<matr,i,vexs,eds,ed,tem>;if f==1flag=0;break;endif dvex~=0ed<:,i>=[vexs<1,i> dvex];vexs<1,i+1>=dvex;matr<vexs<1,i+1>,vexs<1,i>>=0;elsebreak;endendfunction [dvex f]=flecvexf<matr,i,vexs,eds,ed,temp> f=0;edd=find<matr<vexs<1,i>,:>==1>;dvex=0;dvex1=[];ded=[];if length<edd>==1dvex=edd;elsedd=1;dd1=0;kkk=0;for kk=1:length<edd>m1=find<vexs==edd<kk>>;if sum<m1>==0dvex1<dd>=edd<kk>;dd=dd+1;dd1=1;elsekkk=kkk+1;endendif kkk==length<edd>tem=vexs<1,i>*ones<1,kkk>;edd1=[tem;edd];for l1=1:kkklt=0;ddd=1;for l2=1:edsif edd1<1:2,l1>==ed<1:2,l2>lt=lt+1;endendif lt==0ded<ddd>=edd<l1>;ddd=ddd+1;endendendif temp<=length<dvex1>dvex=dvex1<temp>;elseif temp>length<dvex1> & temp<=length<ded>dvex=ded<temp>;elsef=1;endend程序二:Hamilton改良圈算法〔找出比较好的Hamilton路〕function [C d1]= hamiltonglf<v>%d表示权值矩阵%C表示算法最终找到的Hamilton圈.%v =[ 51 67;37 84;41 94;2 99;18 54;4 50;24 42;25 38;13 40;7 64;22 60;25 62;18 40;41 26];n=size<v,1>;subplot<1,2,1>hold on;plot <v<:,1>,v<:,2>,'*'>; %描点for i=1:nstr1='V';str2=num2str<i>;dot=[str1,str2];text<v<i,1>-1,v<i,2>-2,dot>; %给点命名endplot <v<:,1>,v<:,2>>;%连线plot<[v<n,1>,v<1,1>],[v<n,2>,v<1,2>]>;for i =1:nfor j=1:nd<i,j>=sqrt<<v<i,1>-v<j,1>>^2+<v<i,2>-v<j,2>>^2>;endendd2=0;for i=1:nif i<nd2=d2+d<i,i+1>;elsed2=d2+d<n,1>;endendtext<10,30,num2str<d2>>;n=size<d,2>;C=[linspace<1,n,n> 1];for nnn=1:20C1=C;if n>3for m=4:n+1for i=1:<m-3>for j=<i+2>:<m-1>if<d<C<i>,C<j>>+d<C<i+1>,C<j+1>><d<C<i>,C<i+1>>+d<C<j>,C<j+1>>>C1<1:i>=C<1:i>;for k=<i+1>:jC1<k>=C<j+i+1-k>;endC1<<j+1>:m>=C<<j+1>:m>;endendendendelseif n<=3if n<=2fprint<'It does not exist Hamilton circle.'>; elsefprint<'Any cirlce is the right answer.'>;endendC=C1;d1=0;for i=1:nd1=d1+d<C<i>,C<i+1>>;endd1;endsubplot<1,2,2>;hold on;plot <v<:,1>,v<:,2>,'*'>; %描点for i=1:nstr1='V';str2=num2str<i>;dot=[str1,str2];text<v<i,1>-1,v<i,2>-2,dot>; %给点命名endv2=[v;v<1,1>,v<1,2>];plot<v<C<:>,1>,v<C<:>,2>,'r'>;text<10,30,num2str<d1>>;第五讲:匹配问题与算法程序一:较大基础匹配算法function J=matgraf<W>n=size<W,1>;J=zeros<n,n>;while sum<sum<W>>~=0a=find<W~=0>;t1=mod<a<1>,n>;if t1==0t1=n;endif a<1>/n>floor<a<1>/n>t2=floor<a<1>/n>+1;elset2=floor<a<1>/n>;endJ<t1,t2>=1,J<t2,t1>=1;W<t1,:>=0;W<t2,:>=0;W<:,t1>=0;W<:,t2>=0;endJ;程序二:匈牙利算法〔完美匹配算法,包括三个文件fc01,fc02,fc03〕function [e,s]=fc01<a,flag>if nargin==1flag=0;endb=a;if flag==0cmax=max<max<b>'>;b=cmax-b;endm=size<b>;for i =1:m<1>b<i,:>=b<i,:>-min<b<i,:>>;endfor j=1:m<2>b<:,j>=b<:,j>-min<b<:,j>>;endd=<b==0>;[e,total]=fc02<d>;while total~=m<1>b=fc03<b,e>;d=<b==0>;[e,total]=fc02<d>;endinx=sub2ind<size<a>,e<:,1>,e<:,2>>;e=[e,a<inx>];s=sum<a<inx>>;function [e,total]=fc02<d>total=0;m=size<d>;e=zeros<m<1>,2>;t=sum<sum<d>'>;nump=sum<d'>;while t~=0[s,inp]=sort<nump>;inq=find<s>;ep=inp<inq<1>>;inp=find<d<ep,:>>;numq=sum<d<:,inp>>;[s,inq]=sort<numq>;eq=inp<inq<1>>;total=total+1;e<total,:>=[ep,eq];inp=find<d<:,eq>>;nump<inp>=nump<inp>-1;nump<ep>=0;t=t-sum<d<ep,:>>-sum<d<:,eq>>+1;d<ep,:>=0*d<ep,:>;d<:,eq>=0*d<:,eq>;endfunction b=fc03<b,e>m=size<b>;t=1;p=ones<m<1>,1>;q=zeros<m<1>,1>;inp=find<e<:,1>~=0>;p<e<inp,1>>=0;while t~=0tp=sum<p+q>;inp=find<p==1>;n=size<inp>;for i=1:n<1>inq=find<b<inp<i>,:>==0>;q<inq>=1;endinp=find<q==1>;n=size<inp>;for i=1:n<1>if all<e<:,2>-inp<i>>==0inq=find<<e<:,2>-inp<i>>==0>;p<e<inq>>=1;endendtq=sum<p+q>;t=tq-tp;endinp=find<p==1>;inq=find<q==0>;cmin=min<min<b<inp,inq>>'>;inq=find<q==1>;b<inp,:>=b<inp,:>-cmin;b<:,inq>=b<:,inq>+cmin;第六讲:最大流最小费用问题程序一:2F算法<Ford-Fulkerson算法>,求最大流%C=[0 5 4 3 0 0 0 0;0 0 0 0 5 3 0 0;0 0 0 0 0 3 2 0;0 0 0 0 0 0 2 0; %0 0 0 0 0 0 0 4;0 0 0 0 0 0 0 3;0 0 0 0 0 0 0 5;0 0 0 0 0 0 0 0 ] function [f wf]=fulkersonf<C,f1>%C表示容量%f1表示当前流量,默认为0%f表示最大流±íʾ×î´óÁ÷%wf表示最大流的流量n=length<C>;if nargin==1;f=zeros<n,n>;elsef=f1;endNo=zeros<1,n>;d=zeros<1,n>;while <1>No<1>=n+1;d<1>=Inf;while <1>pd=1;for <i=1:n>if <No<i>>for <j=1:n>if <No<j>==0 & f<i,j><C<i,j>>No<j>=i;d<j>=C<i,j>-f<i,j>;pd=0;if <d<j>>d<i>>d<j>=d<i>;endelseif <No<j>==0 & f<j,i>>0>No<j>=-i;d<j>=f<j,i>;pd=0;if <d<j>>d<i>>d<j>=d<i>;endendendendendif <No<n>|pd>break;endendif <pd>break;enddvt=d<n>;t=n;while <1>if<No<t>>0>f<No<t>,t>=f<No<t>,t>+dvt;elseif <No<t><0>f<No<t>,t>=f<No<t>,t>-dvt;endif <No<t>==1>for <i=1:n>No<i>=0;d<i>=0;endbreakendt=No<t>;endendwf=0;for <j=1:n>wf=wf+f<1,j>;endf;wf;程序二:Busacker-Gowan算法<求最大流最小费用>%C=[0 15 16 0 0;0 0 0 13 14;0 11 0 17 0;0 0 0 0 8;0 0 0 0 0] %b=[0 4 1 0 0;0 0 0 6 1;0 2 0 3 0;0 0 0 0 2;0 0 0 0 0]%function [f wf zwf]=BGf<C,b>%C表示弧容量矩阵%b表示弧上单位流量的费用%f表示最大流最小费用矩阵%wf最大流量%zwf表示最小费用n=size<C,2>;wf=0;wf0=inf;f=zeros<n,n>;while <1>a=ones<n,n>*inf;for <i=1:n>a<i,i>=0;endfor <i=1:n>for <j=1:n>if<C<i,j>>0 & f<i,j>==0>a<i,j>=b<i,j>;elseif <C<i,j>>0 & f<i,j>==C<i,j>>a<j,i>=-b<i,j>;elseif <C<i,j>>0>a<i,j>=b<i,j>;a<j,i>=-b<i,j>;endendendfor <i=2:n>p<i>=inf;s<i>=i;endfor <k=1:n>pd=1;for <i=2:n>for <j=1:n>if <p<i>>p<j>+a<j,i>>p<i>=p<j>+a<j,i>;s<i>=j;pd=0; endendendif <pd>break;endendif <p<n>==inf>break;enddvt=inf;t=n;while <1>if <a<s<t>,t>>0>dvtt=C<s<t>,t>-f<s<t>,t>;elseif <a<s<t>,t><0>dvtt=f<t,s<t>>;endif <dvt>dvtt>dvt=dvtt;endif <s<t>==1>break;endt=s<t>;endpd=0;if <wf+dvt>=wf0>dvt=wf0-wf;pd=1;endt=n;while <1>if <a<s<t>,t>>0>f<s<t>,t>=f<s<t>,t>+dvt; elseif <a<s<t>,t><0>f<<t>,s<t>>=f<t,s<t>>-dvt; endif <s<t>==1>break;endt=s<t>;endif <pd>break;endwf=0;for <j=1:n>wf=wf+f<1,j>;endendzwf=0;for <i=1:n>for <j=1:n>zwf=zwf+b<i,j>*f<i,j>;endendf;。
图论编程实现连通图实验报告MATLAB
连通图的判断
一、实验目的
了解连通图的判定方法。
二、实验内容
1、设计一个算法判断图是否连通;
2、在matlab中编程来实现此算法。
用下面的实例来调试程序:
三、使用环境
个人计算机,MATLAB软件
四、编程思路
从邻接矩阵中得到顶点数N,对邻接矩阵求N次方,并求其秩,可以通过邻接矩阵N次幂的秩来判断是否为连通,连通图的秩+1=顶点数,则连通。
五、调试过程
1.程序代码:
function liantong(F)
n=size(F,1);%计算矩阵行数
P=zeros(n,n);%生成同维0矩阵
k=1;
for k=1:n
F1=F^k;
P=P+F1;%将原矩阵的n次方传给P
end
S=n-rank(P);%求矩阵的秩
if S==1
disp('连通');
else
disp('不连通');
end
2.运行窗口:
在运行窗口输入:。
matlab、图论及最短路径算法
mesh(x,y,z)
colormap ezmesh(z)
e03.m
函数画图: ezsurf(z)
数学建模课件
主讲人:孙云龙
4、统计分析
概率分布
pdf cdf inv stat rnd
unif bino poiss norm exp t chi2 F
描述统计
的月降水量为data.m中:矩阵B——行为月数据、列为 12个观测站数据。求:
(1)各观测站月平均降雨量、标准差。 (2)月最大、最小降雨量为多少?发生在哪个气象站的 哪一月份?
mean max
std
min
数学建模课件
主讲人:孙云龙
例2
级数 1 2 3 4 5 6 7 8 9
个人所得税:工资、薪金所得适用
1 0 0 1 1 0 0 0 1 1 0 0 0 1 0 0 1 i点为j边端点 rij R 0 1 1 0 0 0 1 0 0 否则 0 0 1 1 0 0 0 1 0 0 0 0 1 1 1 1
数学建模课件
主讲人:孙云龙
2、连通性
x =‘x’
x =sym(‘x’)
syms x y z
定义符号表达式
字符串 f =’……’ 符号型 运算符 syms x ,f = … f=inline(‘…’)
M函数文件 e02.m
数学建模课件
微积分
主讲人:孙云龙
极限 limit (f,x,a,option) 导数 diff (f,x,n)
税率 (%) 5 10 15 20 25 30 35 40 45
年奖金有什么问题?
数学建模课件
图算法的应用以及在Matlab中的实现
图算法的应用以及在Matlab中的实现图算法是图论的基础,广泛应用于各个领域。
图算法可以用来解决很多实际问题,例如社交网络分析、路网优化、数据挖掘等。
在Matlab中,图算法的实现主要依赖于图对象和图函数库。
一、图算法的应用1.社交网络分析:通过图算法可以进行社交网络的分析和挖掘,例如寻找网络中的关键节点、查找社区结构、计算网络的中心性指标等。
2.路网优化:图算法可以用来解决路网中的最短路径问题、最小生成树问题、最大流问题等,以优化交通运输和资源分配。
3.数据挖掘:图算法可以用于发现模式和关联规则,例如通过挖掘网页链接关系发现重要网页、通过分析推荐系统中用户的行为图谱进行个性化推荐等。
二、Matlab中图算法的实现在Matlab中,图算法的实现主要依赖于图对象(Graph Object)和图函数库(Graph and Network Algorithms),Matlab提供了一套完整的图算法库供用户使用。
1. 图对象(Graph Object)Matlab中的图对象是用来存储图数据的数据结构,可以表示有向图或无向图。
可以通过创建图对象,设置节点和边的属性,以及添加节点和边来构建图。
图对象提供了很多有用的方法,用于访问和操作图数据。
2. 图函数库(Graph and Network Algorithms)Matlab提供了一系列的图函数,用于解决常见的图算法问题。
其中包括最短路径算法、最小生成树算法、最大流算法等。
这些图函数可以用来解决各种实际问题,例如计算两节点之间的最短路径、查找网络中的关键节点等。
在Matlab中,通过以下几个步骤可以实现图算法:1.创建图对象:使用图对象的构造函数可以创建一个空的图对象,并指定图的类型(有向图或无向图)。
2.添加节点和边:使用图对象的方法可以添加节点和边,设置节点和边的属性。
3.访问图数据:可以通过图对象的方法访问和操作图数据,例如获取节点数、获取边数、获取节点的邻居等。
matlab 最短路距离
在MATLAB中,可以使用图论算法来求解最短路问题。
其中,Dijkstra算法是一种常用的最短路算法。
假设我们有一个有向图,其中每条边的权重非负,那么可以使用Dijkstra算法来求解单源最短路问题,即求解从一个顶点到其他所有顶点的最短路径。
以下是一个使用Dijkstra算法求解最短路问题的MATLAB代码示例:matlab复制代码function[dist, path] = dijkstra(adjMatrix, startNode)% 输入:% adjMatrix:邻接矩阵,表示有向图的边权值% startNode:起始节点编号% 输出:% dist:距离矩阵,dist(i,j)表示从起始节点到第i个节点的最短距离% path:路径矩阵,path(i,j)表示从起始节点到第i个节点的前一个节点编号n = size(adjMatrix,1); % 获取顶点数zero_row = find(adjMatrix == 0); % 找到所有不与起始节点相连的行dist = inf(1,n); % 初始化距离矩阵为无穷大dist(startNode) = 0; % 起始节点到自己的距离为0path = zeros(1,n); % 初始化路径矩阵为0prev = zeros(1,n); % 记录前一个节点编号prev(startNode) = -1; % 起始节点的前一个节点编号为-1Q = 1:n; % 待处理的节点集合,初始时为所有节点while ~isempty(Q)[~,min_ind] = min(dist(Q)); % 选择距离最短的节点u = Q(min_ind); % 当前处理的节点编号Q(min_ind) = []; % 从集合中删除该节点neighbors = find(adjMatrix(u,:) > 0); % 找到所有与当前节点相连的节点编号for v = neighborsalt = dist(u) + adjMatrix(u,v); % 计算从起始节点经过u到v的距离if alt < dist(v) % 如果更短,则更新距离和路径dist(v) = alt;path(v) = u;prev(v) = u;if ~ismember(v,Q) % 如果该节点还没有处理过,则加入集合中Q = [Q v]; endendendend。
图论算法及其Matlab程序
求单源最短路径的Dijkstra算法的Matlab程序function [d index1 index2]=Dijkf(a)M=max(max(a));pb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a));d(1:length(a))=M;d(1)=0;temp=1;while sum(pb)<length(a)tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb));tmpb=find(d(tb)==min(d(tb)));temp=tb(tmpb(1));pb(temp)=1;index1=[index1,temp];index=index1(find(d(index1)==d(temp)-a(temp,index1)));if length(index)>=2index=index(1);endindex2(temp)=index;endd;index1;index2;求任意两点间最短路的Floyd算法的Matlab程序function [D,R]=floyd(a)n=size(a,1); D=a;for i=1:nfor j=1:nR(i,j)=j;endendfor k=1:nfor i=1:nfor j=1:nif D(i,k)+D(k,j)<D(i,j)D(i,j)=D(i,k)+D(k,j);R(i,j)=R(i,k);endendendendD;R ;求Euler回路的Fleury算法的Matlab程序function [eu,cEu]=arEuler(E)eu=0;cEu=[];ncV=arComp(E);if max(ncV)>1returnendn=max(max(E(:,1:2)));m=size(E,1);for i=1:nb(i)=0;for j=1:mif E(j,1)==i|E(j,2)==ib(i)=b(i)+1;endendendrp=rem(b,2);srp=sum(rp);switch srpcase 0,eu=1;case 2,eu=0.5;otherwise,returnendif srp==0v1=1;elsev1=find(rp);v1=v1(1);endvc=v1;m=size(E,1);E1=[E(:,1:2),[1:m]'];while ~isempty(E1)evc=find((E1(:,1)==vc)|(E1(:,2)==vc));levc=length(evc);if levc==1cEu=[cEu;E1(evc,3)];vcold=vc;vc=sum(E1(evc,1:2))-vc;E1=E1(setdiff([1:size(E1,1)],evc),:);E2=E1(:,1:2);E2gv=E2>vcold;E2(E2gv)=E2(E2gv)-1;E1(:,1:2)=E2;if vc>vcoldvc=vc-1;endif v1>vcoldv1=v1-1;endelsefor k=1:levcE2=E1(setdiff([1:size(E1,1)],evc(k)),:); ncv=arComp(E2);nco=max(ncv);if (max(ncv)==1)cEu=[cEu;E1(evc(k),3)];vc=sum(E1(evc(k),1:2))-vc;E1=E2;break;endendendendreturn求最小生成树的Prim算法的Matlab程序function [T e]=prim(a)T=[];e=0;v=1;n=size(a,1);c=2:n;for j=2:nb(1,j-1)=1;b(2,j-1)=j;b(3,j-1)=a(1,j);endwhile size(T,2)<n-1[m,i]=min(b(3,:));T(:,size(T,2)+1)=b(:,i);e=e+b(3,i);v=b(2,i);t=find(c==b(2,i));c(t)=[];b(:,i)=[];for j=1:length(c)d=a(v,b(2,j));if d<b(3,j)b(1,j)=v;b(3,j)=d;endendendT;e;求最小生成树的Kruskal算法的Matlab程序function [T c]=kruskal(a)n=size(a,1);m=0;for i=1:n-1for j=i+1:nif a(i,j)>0&a(i,j)<infm=m+1;b(1,m)=i;b(2,m)=j;b(3,m)=a(i,j);endendend[B,i]=sortrows(b',3);B=B';k=0;t=1:n;T=[];c=0;for i=1:mif t(B(1,i))~=t(B(2,i))k=k+1;T(:,k)=B(:,i);c=c+B(3,i);tmin=min(t(B(1,i)),t(B(2,i)));tmax=max(t(B(1,i)),t(B(2,i)));for j=1:nif t(j)==tmaxt(j)=tmin;endendendif k==n-1break;endendt,c求Huffman树的Matlab程序function [h,l]=huffman(p)if (length(find(p<0))~=0)error('Not a prob,negative component');endif (abs(sum(p)-1)>10e-10)error('Not a prob.vector,component do not add to 1') endn=length(p);q=p;m=zeros(n-1,n);for i=1:n-1[q,l]=sort(q);m(i,:)=[l(1:n-i+1),zeros(1,i-1)];q=[q(1)+q(2),q(3:n),1];endfor i=1:n-1c(i,:)=blanks(n*n);endc(n-1,n)='0';c(n-1,2*n)='1';for i=2:n-1c(n-i,1:n-1)=c(n-i+1,n*(find(m(n-i+1,:)==1))...-(n-2):n*(find(m(n-i+1,:)==1)));c(n-i,n)='0';c(n-i,n+1:2*n-1)=c(n-i,1:n-1);c(n-i,2*n)='1';for j=1:i-1c(n-i,(j+1)*n+1:(j+2)*n)=c(n-i+1,...n*(find(m(n-i+1,:)==j+1)-1)+1:n*find(m(n-i+1,:)==j+1));endendfor i=1:nh(i,1:n)=c(1,n*(find(m(1,:)==i)-1)+1:find(m(1,:)==i)*n); ll(i)=length(find(abs(h(i,:))~=32));endl=sum(p.*ll);hl最大流算法Matlab程序function [f wf No]=fofuf(C,f1)n=length(C);if nargin==1;f=zeros(n,n);elsef=f1;endNo=zeros(1,n);d=zeros(1,n);while (1)No(1)=n+1;d(1)=Inf;while (1)pd=1;for (i=1:n)if (No(i))for (j=1:n)if (No(j)==0&f(i,j)<C(i,j))No(j)=i;d(j)=C(i,j)-f(i,j);pd=0;if (d(j)>d(i))d(j)=d(i);endelseif (No(j)==0&f(j,i)>0)No(j)=-i;d(j)=f(j,i);pd=0;if (d(j)>d(i))d(j)=d(i);endendendendendif (No(n)|pd)break;endendif (pd)break;enddvt=d(n);t=n;while (1)if(No(t)>0)f(No(t),t)=f(No(t),t)+dvt;elseif (No(t)<0)f(No(t),t)=f(No(t),t)-dvt;endif (No(t)==1)for (i=1:n)No(i)=0;d(i)=0;endbreakendt=No(t);endendwf=0;for (j=1:n)wf=wf+f(1,j);endf;wf;No;。
图论MATLAB算法
第一章:Dijkstra 算法开始?dot i ≤输入,确定邻接矩阵a 确定邻接矩阵a 的节点数dot1=i算每一个节点到U 中每一个节点的最小值 输出第n i 个节点到第一个节点的最小距离i L ,i=1∧dot结束算这dot 个最小值的最小值l ,并确定其节点位置i nl L i n =将第一个节点放入集合U 中将已经确定的第i n 节点到所有节点的权值赋为∞ 将所有节点到第i n 节点的权值加上l 并代替之 ?2>i1+=i i将第i n 节点放入集合U 中 YESNOYESNO求下面赋权图(左图)中顶点u0到其余顶点的最短路。
其邻接矩阵W 为:⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞∞=024782063446046340357630135102273201847210W)(i u l迭 代 次 数0u 1u 2u 3u 4u 5u 6u 7u1 2 3 4 5 6 7 8 0 ∞ ∞ ∞ ∞ ∞ ∞ ∞1 2 ∞ 7 ∞ 4 82 4 7 ∞ 4 83 7 ∞4 86 9 4 86 9 69 6 9 最后标记)(v l )(v z0 1 2 3 6 9 4 6 0u 0u 0u 2u 3u 3u 0u 6u1u 2u 3u 4u 5u 6u 7u 0ufunction dijkstra%注:此程序仅作参考,欢迎批评指正。
clcclear%Dijkstra算法:%%%%给邻接矩阵赋值%%%%%%%%%%%%a=[0,1,2,inf,7,inf,4,8;1,0,2,3,inf,inf,inf,7;0,0,0,1,5,inf,inf,inf;0,0,0,0,3,6,inf,inf;0,0,0,0,0,4,3,inf;0,0,0,0,0,0,6,4;0,0,0,0,0,0,0,2;];for i=2:8for j=1:i-1a(i,j)=a(j,i);endenddot=size(a,1);%节点数fprintf('\t邻接矩阵的标准形式:');afuquantu=a;%在赋权图中用到fprintf('\t其中,inf代表无穷大∞,a(i,j)代表第i个节点到第j个节点的权。
图论在matlab中的实现
用 index2×n 存放各边端点的信息, 当选中某一边之后,就将此边对 应的顶点序号中较大序号改记为 此边的另一序号,同时把后面边 中所有序号为的改记为。此方法 的几何意义是:将序号的这个顶 点收缩到顶点,顶点不复存在。 后面继续寻查时,发现某边的两 个顶点序号相同时,认为已被收 缩掉,失去了被选取的资格。
while length(result)<loop temp=min(data(3,:)); flag=find(data(3,:)==temp); flag=flag(1); v1=data(1,flag);v2=data(2,flag); if index(1,flag)~=index(2,flag) result=[result,data(:,flag)]; end if v1>v2 index(find(index==v1))=v2; else index(find(index==v2))=v1; end data(:,flag)=[]; index(:,flag)=[]; end result
Kruskal算法如下 Kruskal算法如下 : clc;clear; M=1000; a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40; a(3,4)=52;a(3,7)=45; a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70; [i,j]=find((a~=0)&(a~=M)); b=a(find((a~=0)&(a~=M))); data=[i';j';b'];index=data(1:2,:); loop=max(size(a))loop=max(size(a))-1; result=[];
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
图论算法及其MATLAB 程序代码求赋权图G = (V , E , F )中任意两点间的最短路的Warshall-Floyd 算法:设A = (a ij )n ×n 为赋权图G = (V , E , F )的矩阵, 当v i v j ∈E 时a ij = F (v i v j ), 否则取a ii =0, a ij = +∞(i ≠j ), d ij 表示从v i 到v j 点的距离, r ij 表示从v i 到v j 点的最短路中一个点的编号.① 赋初值. 对所有i , j , d ij = a ij , r ij = j . k = 1. 转向②② 更新d ij , r ij . 对所有i , j , 若d ik + d k j <d ij , 则令d ij = d ik + d k j , r ij = k , 转向③.③ 终止判断. 若d ii <0, 则存在一条含有顶点v i 的负回路, 终止; 或者k = n 终止; 否则令k = k + 1, 转向②.最短路线可由r ij 得到.例1 求图6-4中任意两点间的最短路.解:用Warshall-Floyd 算法, MATLAB 程序代码如下:n=8;A=[0 2 8 1 Inf Inf Inf Inf2 0 6 Inf 1 Inf Inf Inf8 6 0 7 5 1 2 Inf1 Inf 7 0 Inf Inf 9 InfInf 1 5 Inf 0 3 Inf 8Inf Inf 1 Inf 3 0 4 6Inf Inf 2 9 Inf 4 0 3Inf Inf Inf Inf 8 6 3 0]; % MATLAB 中, Inf 表示∞D=A; %赋初值for (i=1:n)for (j=1:n)R(i,j)=j;end ;end %赋路径初值for (k=1:n)for (i=1:n)for (j=1:n)if (D(i,k)+D(k,j)<D(i,j))D(i,j)=D(i,k)+D(k,j); %更新dijR(i,j)=k;end ;end ;end %更新rijk %显示迭代步数D %显示每步迭代后的路长R %显示每步迭代后的路径pd=0;for i=1:n %含有负权时if (D(i,i)<0)pd=1;break ;end ;end %存在一条含有顶点vi 的负回路if (pd)break ;end %存在一条负回路, 终止程序end %程序结束图6-4Kruskal避圈法:将图G中的边按权数从小到大逐条考察, 按不构成圈的原则加入到T 中(若有选择时, 不同的选择可能会导致最后生成树的权数不同), 直到q (T ) = p (G ) − 1为止, 即T的边数= G的顶点数− 1为止.Kruskal避圈法的MATLAB程序代码如下:n=8;A=[0 2 8 1 0 0 0 02 0 6 0 1 0 0 08 6 0 7 5 1 2 01 0 7 0 0 0 9 00 1 5 0 0 3 0 80 0 1 0 3 0 4 60 0 2 9 0 4 0 30 0 0 0 8 6 3 0];k=1; %记录A中不同正数的个数for(i=1:n-1)for(j=i+1:n) %此循环是查找A中所有不同的正数if(A(i,j)>0)x(k)=A(i,j); %数组x记录A中不同的正数kk=1; %临时变量for(s=1:k-1)if(x(k)==x(s))kk=0;break;end;end%排除相同的正数k=k+kk;end;end;endk=k-1 %显示A中所有不同正数的个数for(i=1:k-1)for(j=i+1:k) %将x中不同的正数从小到大排序if(x(j)<x(i))xx=x(j);x(j)=x(i);x(i)=xx;end;end;endT(n,n)=0; %将矩阵T中所有的元素赋值为0q=0; %记录加入到树T中的边数for(s=1:k)if(q==n)break;end%获得最小生成树T, 算法终止for(i=1:n-1)for(j=i+1:n)if (A(i,j)==x(s))T(i,j)=x(s);T(j,i)=x(s); %加入边到树T中TT=T; %临时记录Twhile(1)pd=1;%砍掉TT中所有的树枝for(y=1:n)kk=0;for(z=1:n)if(TT(y,z)>0)kk=kk+1;zz=z;end;end%寻找TT中的树枝if(kk==1)TT(y,zz)=0;TT(zz,y)=0;pd=0;end;end%砍掉TT中的树枝if(pd)break;end;end%已砍掉了TT中所有的树枝pd=0;%判断TT中是否有圈for(y=1:n-1)for(z=y+1:n)if(TT(y,z)>0)pd=1;break;end;end;endif(pd)T(i,j)=0;T(j,i)=0;%假如TT中有圈else q=q+1;end;end;end;end;endT %显示近似最小生成树T, 程序结束求二部图G的最大匹配的算法(匈牙利算法), 其基本思想是:从G的任意匹配M开始, 对X中所有M的非饱和点, 寻找M−增广路. 若不存在M−增广路, 则M为最大匹配; 若存在M−增广路P, 则将P中M与非M的边互换得到比M多一边的匹配M1 , 再对M1重复上述过程.设G = ( X, Y, E )为二部图, 其中X = {x1, x2, … , x n }, Y = { y1, y2, … , y n}. 任取G的一初始匹配M (如任取e∈E, 则M = {e}是一个匹配).①令S = φ , T = φ , 转向②.②若M饱和X \S的所有点, 则M是二部图G的最大匹配. 否则, 任取M的非饱和点u∈X \ S , 令S = S ∪{ u }, 转向③.③记N (S ) = {v | u∈S, uv∈E}. 若N (S ) = T, 转向②. 否则取y∈N (S ) \T. 若y是M 的饱和点, 转向④, 否则转向⑤.④设x y∈M, 则令S = S ∪{ x }, T = T ∪{ y }, 转向③.⑤u −y路是M−增广路, 设为P, 并令M = M⊕P, 转向①. 这里M⊕P = M∪P \M∩P, 是对称差.由于计算M−增广路P比较麻烦, 因此将迭代步骤改为:①将X中M的所有非饱和点(不是M中某条边的端点)都给以标号0和标记*, 转向②.②若X中所有有标号的点都已去掉了标记*, 则M是G的最大匹配. 否则任取X中一个既有标号又有标记*的点x i , 去掉x i的标记*, 转向③.③找出在G中所有与x i邻接的点y j (即x i y j∈E ), 若所有这样的y j都已有标号, 则转向②, 否则转向④.④对与x i邻接且尚未给标号的y j都给定标号i. 若所有的y j都是M的饱和点, 则转向⑤, 否则逆向返回. 即由其中M的任一个非饱和点y j的标号i找到x i, 再由x i的标号k找到y k , … , 最后由y t的标号s找到标号为0的x s时结束, 获得M−增广路x s y t…x i y j, 记P = {x s y t, …, x i y j }, 重新记M为M⊕P, 转向①.⑤将y j在M中与之邻接的点x k (即x k y j∈M), 给以标号j和标记*, 转向②.例1求图6-9中所示的二部图G的最大匹配.图6-9匈牙利算法的MATLAB程序代码如下:m=5;n=5;A=[0 1 1 0 01 1 0 1 10 1 1 0 00 1 1 0 00 0 0 1 1];M(m,n)=0;for(i=1:m)for(j=1:n)if(A(i,j))M(i,j)=1;break;end;end%求初始匹配Mif(M(i,j))break;end;end%获得仅含一条边的初始匹配Mwhile(1)for(i=1:m)x(i)=0;end%将记录X中点的标号和标记*for(i=1:n)y(i)=0;end%将记录Y中点的标号和标记*for(i=1:m)pd=1;%寻找X中M的所有非饱和点for(j=1:n)if(M(i,j))pd=0;end;endif(pd)x(i)=-n-1;end;end%将X中M的所有非饱和点都给以标号0和标记*, 程序中用n+1表示0标号, 标号为负数时表示标记*pd=0;while(1)xi=0;for(i=1:m)if(x(i)<0)xi=i;break;end;end%假如X中存在一个既有标号又有标记*的点, 则任取X中一个既有标号又有标记*的点xiif(xi==0)pd=1;break;end%假如X中所有有标号的点都已去掉了标记*, 算法终止x(xi)=x(xi)*(-1); %去掉xi的标记*k=1;for(j=1:n)if(A(xi,j)&y(j)==0)y(j)=xi;yy(k)=j;k=k+1;end;end%对与xi邻接且尚未给标号的yj都给以标号iif(k>1)k=k-1;for(j=1:k)pdd=1;for(i=1:m)if(M(i,yy(j)))x(i)=-yy(j);pdd=0;break;end;end%将yj在M中与之邻接的点xk (即xkyj∈M), 给以标号j和标记*if(pdd)break;end;endif(pdd)k=1;j=yy(j); %yj不是M的饱和点while(1)P(k,2)=j;P(k,1)=y(j);j=abs(x(y(j))); %任取M的一个非饱和点yj, 逆向返回if(j==n+1)break;end%找到X中标号为0的点时结束, 获得M-增广路Pk=k+1;endfor(i=1:k)if(M(P(i,1),P(i,2)))M(P(i,1),P(i,2))=0; %将匹配M在增广路P中出现的边去掉else M(P(i,1),P(i,2))=1;end;end%将增广路P中没有在匹配M中出现的边加入到匹配M中break;end;end;endif(pd)break;end;end%假如X中所有有标号的点都已去掉了标记*, 算法终止M %显示最大匹配M, 程序结束利用可行点标记求最佳匹配的算法步骤如下:设G = ( X , Y , E , F )为完备的二部赋权图, L 是其一个初始可行点标记, 通常取.,,0)(},|)(max{)(Y y X x y L Y y xy F x L ∈∈ =∈= M 是G L 的一个匹配. ① 若X 的每个点都是M 的饱和点, 则M 是最佳匹配. 否则取M 的非饱和点u ∈X , 令S = {u }, T = φ , 转向②.② 记N L (S ) = {v | u ∈S , uv ∈E L }. 若N L ( S ) = T , 则G L 没有完美匹配, 转向③. 否则转向④.③ 调整可行点标记, 计算a L = min { L ( x ) + L ( y ) − F (x y ) | x ∈S , y ∈Y \T }.由此得新的可行顶点标记H (v ) =,,),(,)(,)(T v S v v L a v L a v L L L ∈∈+−令L = H , G L = G H , 重新给出G L 的一个匹配M , 转向①.④ 取y ∈N L ( S ) \T , 若y 是M 的饱和点, 转向⑤. 否则, 转向⑥.⑤ 设x y ∈M , 则令S = S ∪{ x }, T = T ∪{ y }, 转向②.⑥ 在G L 中的u − y 路是M −增广路, 记为P , 并令 M = M ⊕P , 转向①.利用可行点标记求最佳匹配算法的MATLAB 程序代码如下:n=4;A=[4 5 5 12 2 4 64 2 3 35 0 2 1];for (i=1:n)L(i,1)=0;L(i,2)=0;endfor (i=1:n)for (j=1:n)if (L(i,1)<A(i,j))L(i,1)=A(i,j);end ; %初始可行点标记LM(i,j)=0;end ;endfor (i=1:n)for (j=1:n) %生成子图Glif (L(i,1)+L(j,2)==A(i,j))Gl(i,j)=1;else Gl(i,j)=0;end ;end ;endii=0;jj=0;for (i=1:n)for (j=1:n)if (Gl(i,j))ii=i;jj=j;break ;end ;endif (ii)break ;end ;end %获得仅含Gl 的一条边的初始匹配MM(ii,jj)=1;for (i=1:n)S(i)=0;T(i)=0;NlS(i)=0;endwhile (1)for (i=1:n)k=1;否则.for(j=1:n)if(M(i,j))k=0;break;end;endif(k)break;end;endif(k==0)break;end%获得最佳匹配M, 算法终止S(1)=i;jss=1;jst=0;%S={xi}, T=φwhile(1)jsn=0;for(i=1:jss)for(j=1:n)if(Gl(S(i),j))jsn=jsn+1;NlS(jsn)=j;%NL(S)={v|u∈S,uv∈EL}for(k=1:jsn-1)if(NlS(k)==j)jsn=jsn-1;end;end;end;end;endif(jsn==jst)pd=1; %判断NL(S)=T?for(j=1:jsn)if(NlS(j)~=T(j))pd=0;break;end;end;endif(jsn==jst&pd)al=Inf; %如果NL(S)=T, 计算al, Inf为∞for(i=1:jss)for(j=1:n)pd=1;for(k=1:jst)if(T(k)==j)pd=0;break;end;endif(pd&al>L(S(i),1)+L(j,2)-A(S(i),j))al=L(S(i),1)+L(j,2)-A(S(i),j);end;end;end for(i=1:jss)L(S(i),1)=L(S(i),1)-al;end%调整可行点标记for(j=1:jst)L(T(j),2)=L(T(j),2)+al;end%调整可行点标记for(i=1:n)for(j=1:n) %生成子图GLif(L(i,1)+L(j,2)==A(i,j))Gl(i,j)=1;else Gl(i,j)=0;endM(i,j)=0;k=0;end;endii=0;jj=0;for(i=1:n)for(j=1:n)if(Gl(i,j))ii=i;jj=j;break;end;endif(ii)break;end;end%获得仅含Gl的一条边的初始匹配MM(ii,jj)=1;breakelse%NL(S)≠Tfor(j=1:jsn)pd=1;%取y∈NL(S)\Tfor(k=1:jst)if(T(k)==NlS(j))pd=0;break;end;endif(pd)jj=j;break;end;endpd=0;%判断y是否为M的饱和点for(i=1:n)if(M(i,NlS(jj)))pd=1;ii=i;break;end;endif(pd)jss=jss+1;S(jss)=ii;jst=jst+1;T(jst)=NlS(jj); %S=S∪{x}, T=T∪{y}else%获得Gl的一条M-增广路, 调整匹配Mfor(k=1:jst)M(S(k),T(k))=1;M(S(k+1),T(k))=0;endif(jst==0)k=0;endM(S(k+1),NlS(jj))=1;break;end;end;end;endMaxZjpp=0;for(i=1:n)for(j=1:n)if(M(i,j))MaxZjpp=MaxZjpp+A(i,j);end;end;endM %显示最佳匹配MMaxZjpp %显示最佳匹配M的权, 程序结束从一个可行流f 开始, 求最大流的Ford--Fulkerson 标号算法的基本步骤:⑴ 标号过程① 给发点v s 以标号(+, +∞) , δ s = +∞.② 选择一个已标号的点x , 对于x 的所有未给标号的邻接点y , 按下列规则处理:当yx ∈E , 且f yx >0时, 令δ y = min { f yx , δ x }, 并给y 以标号 ( x − , δ y ).当xy ∈E , 且f xy <C xy 时, 令δ y = min {C xy − f xy , δ x }, 并给y 以标号 ( x + , δ y ). ③ 重复②直到收点v t 被标号或不再有点可标号时为止. 若v t 得到标号, 说明存在一条可增广链, 转⑵调整过程; 若v t 未得到标号, 标号过程已无法进行时, 说明f 已经是最大流.⑵ 调整过程④ 决定调整量δ =δ vt , 令u = v t .⑤ 若u 点标号为( v +, δ u ), 则以f vu + δ 代替f vu ; 若u 点标号为( v −, δ u ), 则以 f vu − δ 代替f vu .⑥ 若v = v s , 则去掉所有标号转⑴重新标号; 否则令u = v , 转⑤.算法终止后, 令已有标号的点集为S , 则割集(S , S c )为最小割, 从而W f = C (S , S c ). 例1 求图6-19所示网络的最大流.利用Ford--Fulkerson 标号法求最大流算法的MATLAB 程序代码如下:n=8;C=[0 5 4 3 0 0 0 00 0 0 0 5 3 0 00 0 0 0 0 3 2 00 0 0 0 0 0 2 00 0 0 0 0 0 0 40 0 0 0 0 0 0 30 0 0 0 0 0 0 50 0 0 0 0 0 0 0]; %弧容量for (i=1:n)for (j=1:n)f(i,j)=0;end ;end %取初始可行流f 为零流for (i=1:n)No(i)=0;d(i)=0;end %No,d 记录标号图6-19while(1)No(1)=n+1;d(1)=Inf; %给发点vs标号while(1)pd=1;%标号过程for(i=1:n)if(No(i)) %选择一个已标号的点vifor(j=1:n)if(No(j)==0&f(i,j)<C(i,j)) %对于未给标号的点vj, 当vivj为非饱和弧时No(j)=i;d(j)=C(i,j)-f(i,j);pd=0;if(d(j)>d(i))d(j)=d(i);endelseif(No(j)==0&f(j,i)>0) %对于未给标号的点vj, 当vjvi为非零流弧时No(j)=-i;d(j)=f(j,i);pd=0;if(d(j)>d(i))d(j)=d(i);end;end;end;end;endif(No(n)|pd)break;end;end%若收点vt得到标号或者无法标号, 终止标号过程if(pd)break;end%vt未得到标号, f已是最大流, 算法终止dvt=d(n);t=n; %进入调整过程, dvt表示调整量while(1)if(No(t)>0)f(No(t),t)=f(No(t),t)+dvt; %前向弧调整elseif(No(t)<0)f(No(t),t)=f(No(t),t)-dvt;end%后向弧调整if(No(t)==1)for(i=1:n)No(i)=0;d(i)=0; end;break;end%当t的标号为vs时, 终止调整过程t=No(t);end;end; %继续调整前一段弧上的流fwf=0;for(j=1:n)wf=wf+f(1,j);end%计算最大流量f %显示最大流wf %显示最大流量No %显示标号, 由此可得最小割, 程序结束设网络G = ( V , E , C ), 取初始可行流 f 为零流, 求解最小费用流问题的迭代步骤: ① 构造有向赋权图 G f = ( V , E f , F ), 对于任意的v i v j ∈E , E f , F 的定义如下:当f ij = 0时, v i v j ∈E f , F ( v i v j ) = b ij ;当f ij = C ij 时, v j v i ∈E f , F ( v j v i ) = −b ij ;当0< f ij <C ij 时, v i v j ∈E f , F ( v i v j ) = b ij , v j v i ∈E f , F ( v j v i ) = −b ij .转向②.② 求出有向赋权图G f = (V , E f , F )中发点v s 到收点v t 的最短路µ , 若最短路µ存在转向③; 否则f 是所求的最小费用最大流, 停止.③ 增流. 同求最大流的方法一样, 重述如下:令.,,,−+∈∈ −=µµδj i j i ij ij ij ij v v v v f f C δ = min {δ ij | v i v j ∈µ}, 重新定义流f = { f ij }为 f ij =,,,,−+∈∈ −+µµδδj i j i ijij ij v v v v f f f如果W f 大于或等于预定的流量值, 则适当减少δ 值, 使W f 等于预定的流量值, 那么 f 是所求的最小费用流, 停止; 否则转向①.求解含有负权的有向赋权图G = ( V , E , F )中某一点到其它各点最短路的Ford 算法. 当v i v j ∈E 时记w ij = F (v i v j ), 否则取w ii =0, w ij = +∞(i ≠j ). v 1到v i 的最短路长记为π ( i ), v 1到v i 的最短路中v i 的前一个点记为θ ( i ). Ford 算法的迭代步骤:① 赋初值π (1) = 0, π ( i ) = +∞, θ ( i ) = i , i = 2, 3, … , n .② 更新π ( i ), θ ( i ). 对于i = 2, 3, … , n 和j = 1, 2, … , n , 如果π ( i )<π ( j ) + w ji , 则令π ( i ) = π ( j ) , θ ( i ) = j . ③ 终止判断:若所有的π ( i )都无变化, 停止; 否则转向②. 在算法的每一步中, π ( i )都是从v 1到v i 的最短路长度的上界. 若不存在负长回路, 则从v 1到v i 的最短路长度是π ( i )的下界, 经过n −1次迭代后π ( i )将保持不变. 若在第n 次迭代后π ( i )仍在变化时, 说明存在负长回路.其它.例2 在图6-22所示运输网络上, 求s 到t 的最小费用最大流, 括号内为(C ij , b ij ).求最小费用最大流算法的MATLAB 程序代码如下:n=5;C=[0 15 16 0 00 0 0 13 140 11 0 17 00 0 0 0 80 0 0 0 0]; %弧容量b=[0 4 1 0 00 0 0 6 10 2 0 3 00 0 0 0 20 0 0 0 0]; %弧上单位流量的费用wf=0;wf0=Inf; %wf 表示最大流量, wf0表示预定的流量值for (i=1:n)for (j=1:n)f(i,j)=0;end ;end %取初始可行流f 为零流while (1)for (i=1:n)for (j=1:n)if (j~=i)a(i,j)=Inf;end ;end ;end %构造有向赋权图for (i=1:n)for (j=1:n)if (C(i,j)>0&f(i,j)==0)a(i,j)=b(i,j);elseif (C(i,j)>0&f(i,j)==C(i,j))a(j,i)=-b(i,j);elseif (C(i,j)>0)a(i,j)=b(i,j);a(j,i)=-b(i,j);end ;end ;endfor (i=2:n)p(i)=Inf;s(i)=i;end %用Ford 算法求最短路, 赋初值for (k=1:n)pd=1; %求有向赋权图中vs 到vt 的最短路for (i=2:n)for (j=1:n)if (p(i)>p(j)+a(j,i))p(i)=p(j)+a(j,i);s(i)=j;pd=0;end ;end ;endif (pd)break ;end ;end %求最短路的Ford 算法结束if (p(n)==Inf)break ;end %不存在vs 到vt 的最短路, 算法终止. 注意在求最小费用最大流时构造有向赋权图中不会含负权回路, 所以不会出现k=ndvt=Inf;t=n; %进入调整过程, dvt 表示调整量while (1) %计算调整量if (a(s(t),t)>0)dvtt=C(s(t),t)-f(s(t),t); %前向弧调整量elseif (a(s(t),t)<0)dvtt=f(t,s(t));end %后向弧调整量if (dvt>dvtt)dvt=dvtt;endif (s(t)==1)break ;end %当t 的标号为vs 时, 终止计算调整量t=s(t);end %继续调整前一段弧上的流fpd=0;if (wf+dvt>=wf0)dvt=wf0-wf;pd=1;end %如果最大流量大于或等于预定的流量值t=n;while (1) %调整过程if (a(s(t),t)>0)f(s(t),t)=f(s(t),t)+dvt; %前向弧调整elseif (a(s(t),t)<0)f(t,s(t))=f(t,s(t))-dvt;end %后向弧调整if (s(t)==1)break ;end %当t 的标号为vs 时, 终止调整过程t=s(t);endif (pd)break ;end %如果最大流量达到预定的流量值wf=0; for (j=1:n)wf=wf+f(1,j);end ;end %计算最大流量zwf=0;for (i=1:n)for (j=1:n)zwf=zwf+b(i,j)*f(i,j);end ;end %计算最小费用f %显示最小费用最大流图6-22wf %显示最小费用最大流量zwf %显示最小费用, 程序结束。