D_i_j_k_s_t_r_a最短路算法MATLAB程序_

最短路dijkstra算法Matlab程序

function [c0,c,path0,path]=dijkstra(s,t,C,flag)% Use the Dijkstra's algorithm to find the shortest path from% s to t and can also find the shortest path between s and all% the other points.% Reference: Graph Theory with Applications by J. A. Bondy and% U. S. R. Murty.% Input -- s is the starting point and also is the point s.% -- t is the given terminal point and is the point t.% -- C \in R^{n \times n}is the cost matrix, where% C(i,j)>=0 is the cost from point i to point j.% If there is no direct connection between point i and% j, C(i,j)=inf.% -- flag: if flag=1, the function just reports the% shortest path between s and t; if flag~=1, the% function reports the shortest path between s and t,% and the shortest paths between s and other points.% Output -- c0 is the minimal cost from s to t.% -- path0 denotes the shortest path form s to t.% -- c \in R{1\times n} in which the element i is the% minimal cost from s to point i.% -- path \in R^{n \times n} in which the row i denotes% the shortest path from s to point i.% Copyright by MingHua Xu(徐明华), Changhzou University, 27 Jan. 2014. s=floor(s);t=floor(t);n=size(C,1);if s<1 || t < 1 || s > n || t > nerror(' The starting point and the terminal point exceeds the valid range');endif t==sdisp('The starting point and the terminal point are the same points');endlabel=ones(1,n)*inf;label(s)=0;S=[s];Sbar=[1:s-1,s+1:n];c0=0;path=zeros(n,n);path(:,1)=s;c=ones(1,n)*inf;parent=zeros(1,n);i=1; % number of points in point set S.while i<n% for each point in Sbar, replace label(Sbar(j)) by% min(label(Sbar(j)),label(S(k))+C(S(k),Sbar(j)))for j=1:n-ifor k=1:iif label(Sbar(j)) > label(S(k))+C(S(k),Sbar(j))label(Sbar(j))=label(S(k))+C(S(k),Sbar(j));parent(Sbar(j))=S(k);endendend% Find the minmal label(j), j \in Sbar.temp=label(Sbar(1));son=1;for j=2:n-iif label(Sbar(j))< temptemp=label(Sbar(j));son=j;endend% update the point set S and SbarS=[S,Sbar(son)];Sbar=[Sbar(1:son-1),Sbar(son+1:n-i)];i=i+1;% if flag==1, just output the shortest path between s and t.if flag==1 && S(i)==tson=t;temp_path=[son];if son~=swhile parent(son)~=sson=parent(son);temp_path=[temp_path,son];endtemp_path=[temp_path,s];endtemp_path=fliplr(temp_path);m=size(temp_path,2);path0(1:m)=temp_path;c_temp=0;for j=1:m-1c_temp=c_temp+C(temp_path(j),temp_path(j+1));endc0=c_temp;path(t,1:m)=path0;c(t)=c0;returnendend% Form the output resultsfor i=1:nson=i;temp_path=[son];if son~=swhile parent(son)~=sson=parent(son);temp_path=[temp_path,son];endtemp_path=[temp_path,s];endtemp_path=fliplr(temp_path);m=size(temp_path,2);path(i,1:m)=temp_path;c_temp=0;for j=1:m-1c_temp=c_temp+C(temp_path(j),temp_path(j+1));endc(i)=c_temp;c0=c(t);path0=path(t,:);endreturn。

Dijkstra、Floyd算法Matlab_Lingo实现

Dijkstra算法Matlab实现。

%求一个点到其他各点的最短路径function [min,path]=dijkstra(w,start,terminal)%W是邻接矩阵%start是起始点Array %terminal是终止点%min是最短路径长度%path是最短路径n=size(w,1);label(start)=0;f(start)=start;for i=1:nif i~=startlabel(i)=inf;endends(1)=start;u=start;while length(s)<nfor i=1:nins=0;forif i==s(j)ins=1;endendif ins==0v=i;if label(v)>(label(u)+w(u,v))label(v)=(label(u)+w(u,v));f(v)=u;endendendv1=0;k=inf;for i=1:nins=0;for j=1:length(s)if i==s(j)ins=1;endend-if ins==0v=i;if k>label(v)k=label(v);v1=v;endendends(length(s)+1)=v1;u=v1;endmin=label(terminal);path(1)=terminal;i=1;while path(i)~=startpath(i+1)=f(path(i));i=i+1 ;endpath(i)=start;L=length(path);path=path(L:-1:1);Floyd算法：matlab程序：%floyd算法，function [D,path,min1,path1]=floyd(a,start,terminal)%a是邻接矩阵%start是起始点%terminal是终止点%D是最小权值表D=a;n=size(D,1);path=zeros(n,n);for i=1:nfor j=1:nif D(i,j)~=infpath(i,j)=j;endendendfor 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);path(i,j)=path(i,k);endendendendif nargin==3min1=D(start,terminal);m(1)=start;i=1;path1=[ ];while path(m(i),terminal)~=terminalk=i+1;m(k)=path(m(i),terminal);i=i+1;endm(i+1)=terminal;path1=m;end1 6 5 5 5 66 2 3 4 4 65 2 3 4 5 45 2 3 4 5 61 4 3 4 5 11 2 4 4 1 6Floyd算法：Lingo程序：!用LINGO11.0编写的FLOYD算法如下;model:sets:nodes/c1..c6/;link(nodes,nodes):w,path; !path标志最短路径上走过的顶点;endsetsdata:path=0;w=0;@text(mydata1.txt)=@writefor(nodes(i):@writefor(nodes(j):-@format(w(i,j),' 10.0f')),@newline(1));@text(mydata1.txt)=@write(@newline(1));@text(mydata1.txt)=@writefor(nodes(i):@writefor(nodes(j):@format(path(i,j),' 10.0f')),@newline(1));enddatacalc:w(1,2)=50;w(1,4)=40;w(1,5)=25;w(1,6)=10;w(2,3)=15;w(2,4)=20;w(2,6)=25;w(3,4)=10;w(3,5)=20;w(4,5)=10;w(4,6)=25;w(5,6)=55;@for(link(i,j):w(i,j)=w(i,j)+w(j,i));@for(link(i,j) |i#ne#j:w(i,j)=@if(w(i,j)#eq#0,10000,w(i,j)));@for(nodes(k):@for(nodes(i):@for(nodes(j):tm=@smin(w(i,j),w(i,k)+w(k,j));path(i,j)=@if(w(i,j)#gt# tm,k,path(i,j));w(i,j)=tm)));endcalcend无向图的最短路问题Lingomodel:sets:cities/1..5/;roads(cities,cities):w,x;endsetsdata:w=0;enddatacalc:w(1,2)=41;w(1,3)=59;w(1,4)=189;w(1,5)=81;w(2,3)=27;w(2,4)=238;w(2,5)=94;w(3,4)=212;w(3,5)=89;w(4,5)=171;@for(roads(i,j):w(i,j)=w(i,j)+w(j,i));@for(roads(i,j):w(i,j)=@if(w(i,j) #eq# 0, 1000,w(i,j)));endcalcn=@size(cities); !城市的个数;min=@sum(roads:w*x);@for(cities(i)|i #ne#1 #and# i #ne#n:@sum(cities(j):x(i,j))=@sum(cities(j):x(j,i)));@sum(cities(j):x(1,j))=1;-@sum(cities(j):x(j,1))=0; !不能回到顶点1;@sum(cities(j):x(j,n))=1;@for(roads:@bin(x));endLingo编的sets:dian/a b1 b2 c1 c2 c3 d/:;link(dian,dian)/a,b1 a,b2 b1,c1 b1,c2 b1,c3 b2,c1 b2,c2 b2,c3 c1,d c2,d c3,d/:x,w;endsetsdata:w=2 4 3 3 1 2 3 1 1 3 4;enddatamin=@sum(link:w*x);@for(link:@bin(x));n=@size(dian);@sum(link(i,j)|i#eq#1:x(i,j))=1;@sum(link(j,i)|i#eq#n:x(j,i))=1;@for(dian(k)|k#ne#1#and#k#ne#n:@sum(link(i,k):x(i,k))=@sum(link(k,i):x(k,i)));- sets:dian/1..5/:level; !level(i)表示点i的水平，用来防止生产圈;link(dian,dian):d,x;endsetsdata:d=0 41 59 189 8141 0 27 238 9459 27 0 212 89189 238 212 0 17181 94 89 171 0;enddatan=@size(dian);min=@sum(link(i,j)|i#ne#j:d(i,j)*x(i,j));@sum(dian(j)|j#gt#1:x(1,j))>1;@for(dian(i)|i#gt#1:@sum(dian(j)|j#ne#i:x(j,i))=1);@for(dian(i)|i#gt#1:@for(dian(j)|j#ne#i#and#j#gt#1:level(j)>level(i)+x(i,j)-(n-2)*(1-x(i,j))+(n-3)*x(j, i)));@for(dian(i)|i#gt#1:level(i)<n-1-(n-2)*x(1,i));@for(dian(i)|i#gt#1:@bnd(1,level(i),100000));@for(link:@bin(x));。

电力系统三相短路计算的MATLAB代码

. . . .电力系统三相短路计算a main.mclear tim%打开文件[dfile,pathname]=uigetfile('*.m','Select Data File');if pathname == 0error(' you must select a valid data file')elselfile =length(dfile);eval(dfile(1:lfile-2));end%定义输出文件output_file=fopen('output.dat','w');%开始计时tic;%求解节点导纳矩阵，其中Ymatrix1是考虑了变比，且支路未近似的导纳矩阵；Ymatrix2是近似变比为1，但是支路未近似计算的节点导纳矩阵；Ymatrix3是近似变比为1，采取近似支路参数1的导纳矩阵；Ymatrix4是近似变比为1，采取近似支路参数2的导纳矩阵。

Y = Ymatrix2(bus,line);%对故障点进行导纳修正fixY = FixY(Y,bus,fault);%求注入电流Iinj = Inode(bus,calcSettings);U = fixY\Iinj;%得到故障支路与其他支路电流Bcurrent = Ibranch( line,U,fault,Y );%如果发生支路三相短路，那么对应该支路的电流修正为-999999-j999999Ib = ReviseBcurrent( fault,Bcurrent );%结束计时tim=toc;fprintf('\n程序运行结果');fprintf('\n计算完成，共用时%4.4fs,相关结果已保存在output.dat\n',tim);%输出结果fprintf_result(output_file, Ib);fprintf_result1(Ib);b FixY.mfunction fixY = FixY( Y,bus,fault )%对形成的导纳矩阵进行故障点的修正[nb,mb]=size(bus);[nf,mf]= size(fault);fixY = Y;%对发电机节点导纳修正for k=1:nbbusType=bus(k,7);if (busType==1)fixY(bus(k,1),bus(k,1)) = fixY(bus(k,1),bus(k,1)) + 1/1i/bus(k,8);endend%对节点短路和支路短路的导纳矩阵进行修正for k=1:nfnodeI=fault(k,1);nodeJ=fault(k,2);dis=fault(k,3);if (nodeI==0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue;endif (nodeJ==0)fixY(nodeI,nodeI) = 999999+1i*999999;continue;endif (dis==0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 999999+1i*999999;continue;endif (dis==1)&&(nodeI*nodeJ~=0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue;endif (dis~=1)&&(dis~=0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = fixY(nodeI,nodeI) - fixY(nodeI,nodeJ)/dis;fixY(nodeJ,nodeJ) = fixY(nodeJ,nodeJ) - fixY(nodeI,nodeJ)/(1-dis);fixY(nodeI,nodeJ)=0;fixY(nodeJ,nodeI)=0;endendendc fprintf_result.mfunction [ output_args ] = fprintf_result( output_file, Ib )%将得到的短路电流输入到输出文件中[n,m]=size(Ib);fprintf( output_file, ' No. No. vector of I value of I\n');for k=1:nI=Ib(k,1);J=Ib(k,2);I01=real(Ib(k,3));I02=imag(Ib(k,3));I1=Ib(k,4);if(I02>=0)fprintf( output_file, '%3d %3d %10.6f+j%10.6f %10.6f',I,J,I01,I02,I1);endif(I02<0)I02=abs(I02);fprintf( output_file, '%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf( output_file, '\n');endendd fprintf_result1.mfunction [ output_args ] = fprintf_result1( Ib )%UNTITLED ÇëÔÚ´Ë´¦ÊäÈëº¯Êý¸ÅÒª[n,m]=size(Ib);fprintf(' No. No. vector of I value of I\n');for k=1:nI=Ib(k,1);J=Ib(k,2);I01=real(Ib(k,3));I02=imag(Ib(k,3));I1=Ib(k,4);if(I02>=0)fprintf('%3d %3d %10.6f+j%10.6f %10.6f',I,J,I01,I02,I1);endif(I02<0)I02=abs(I02);fprintf('%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf('\n');endende Ibranch.mfunction Bcurrent = Ibranch( line,U,fault,Y )%计算短路电流%记录短路故障参数，如短路节点，如为支路短路，记录距离节点的距离%此段计算采用的支路参数未近似，如果计算近似的时候需要修改[nl,ml]=size(line);Bcurrent=zeros(nl+1,4);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);faultNode = 0;if(faultI==0)faultNode = faultJ;endif(faultJ==0)faultNode = faultI;endif(dis==1)&&(faultI*faultJ~=0)faultNode = faultJ;endif(dis==0)&&(faultI*faultJ~=0)faultNode = faultI;endif(faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Iij = 0;Iij1=0;end%计算非故障支路的短路电流for k=1:nli=line(k,1);j=line(k,2);Ui=U(i);if j~=0Uj=U(j);elseUj=0;endif line(k,2)==0Ym=line(k,5)+1i*line(k,6);Iij=Ui*Ym;Iij1=abs(Iij);endif line(k,2)~=0Zt=line(k,3)+1i*line(k,4);Yt=1/Zt;Ym=line(k,5)+1i*line(k,6);Iij=(Ui-Uj)*Yt+Ui*Ym;Iij1=abs(Iij);endBcurrent(k,1)=i;Bcurrent(k,2)=j;Bcurrent(k,3)=Iij;Bcurrent(k,4)=Iij1;end%如果为节点短路，修正短路点的电流大小if(faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Ifault = 0;branchCurrent=0;for k=1:nlI=line(k,1);J=line(k,2);if(I*J==0)continue;endbranchCurrent = (U(I)-U(J))/(line(k,3)+1i*line(k,4));if (I==faultNode)Ifault = Ifault - branchCurrent ;elseif (J==faultNode)Ifault = Ifault + branchCurrent ;endendendBcurrent(nl+1,3) = Ifault;Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));end%如果为支路短路，修正短路支路的短路电流大小if(dis~=0)&&(dis~=1)&&(faultI*faultJ~=0)Bcurrent(nl+1,1) = faultI;Bcurrent(nl+1,2) = faultJ;Bcurrent(nl+1,3) = U(faultI)*Y(faultI,faultJ)/dis + U(faultJ)*Y(faultI,faultJ)/(1-dis);Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));endendf Inode.mfunction Iinj = Inode( bus,calcSettings )%计算节点注入电流[nb,mb]=size(bus);Iinj = zeros(nb,1);for k=1:nbbusType=bus(k,7);if(calcSettings(1)==1)v = 1;elsev = bus(k,2);end%对发电机节点电流进行修正if (busType==1)Iinj(bus(k,1),1) = Iinj(bus(k,1),1) + v/1i/bus(k,8);endendendg ReviseBcurrent.mfunction Ib = ReviseBcurrent( fault,Bcurrent )%如果发生支路短路，对原来的计算电流进行修正，使该支路短路电流输出为-999999-j999999 clear faultIfaultJdis[nt,mt]=size(Bcurrent);Ib=zeros(nt,mt);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);for k=1:nt-1i=Bcurrent(k,1);j=Bcurrent(k,2);Ib(k,:)=Bcurrent(k,:);if (faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultI)&&(j==faultJ)Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endif (faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultJ)&&(j==faultI)Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endIb(nt,:)=Bcurrent(nt,:);endh Ymatrix1.mfunction Y = Ymatrix1( bus,line )%考虑变压器，并且支路参数不近似的节点导纳矩阵[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=1/Zt;Ym=line(k,5)+1i*line(k,6);K=line(k,7);if (K==0)&&(J~=0)Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif (K==0)&&(J==0)Y(I,I)=Y(I,I)+Ym;endif K>0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/(K*K);Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endif K<0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)+K*Yt;Y(J,I)=Y(I,J);endendendi Ymatrix2.mfunction Y = Ymatrix2( bus,line )%考虑变压器变比近似为1，支路参数不等效[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=1/Zt;Ym=line(k,5)+1i*line(k,6);if J~=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif J==0Y(I,I)=Y(I,I)+Ym;endendendj Ymatrix3.mfunction Y = Ymatrix3( bus,line )%考虑变压器变比为1，采用支路参数近似1[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=imag(1/Zt);Ym=imag(line(k,5)+1i*line(k,6));if J~=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif J==0Y(I,I)=Y(I,I)+Ym;endendendk Ymatrix4.mfunction Y = Ymatrix4( bus,line )%变压器变比近似为1，采用支路等效参数2 [nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=1i*line(k,4);Yt=1/Zt;Ym=1i*line(k,6);if J~=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif J==0Y(I,I)=Y(I,I)+Ym;endendend。

matlab dijkstra算法求解最短路径例题

matlab dijkstra算法求解最短路径例题摘要：一、Dijkstra 算法简介1.Dijkstra 算法背景2.Dijkstra 算法原理二、MATLAB 实现Dijkstra 算法求解最短路径1.创建图对象2.计算最短路径3.可视化结果三、Dijkstra 算法应用示例1.例题描述2.解题步骤3.结果分析正文：一、Dijkstra 算法简介Dijkstra 算法是一种经典的图论算法，用于计算图中两个节点之间的最短路径。

它是由荷兰计算机科学家Edsger W.Dijkstra 于1956 年提出的，其基本思想是以起始点为中心向外层层扩展，直到扩展到终点为止。

Dijkstra 算法能得出最短路径的最优解，但由于它遍历计算的节点很多，所以效率低。

可以用堆优化来提高效率。

二、MATLAB 实现Dijkstra 算法求解最短路径1.创建图对象首先，我们需要使用MATLAB 的graph 函数创建一个图对象，指定节点和边的信息。

例如，我们创建一个简单的图，包含4 个节点和3 条边：```matlabG = graph(4, 3);```其中，4 表示图中有4 个节点，3 表示图中有3 条边。

2.计算最短路径接下来，我们可以使用MATLAB 的shortestpath 函数计算两个节点之间的最短路径。

例如，我们计算节点1 到节点3 的最短路径：```matlabSP = shortestpath(G, 1, 3);```3.可视化结果最后，我们可以使用MATLAB 的plot 函数将最短路径可视化。

例如，我们绘制节点和边以及最短路径：```matlabplot(G, SP);```三、Dijkstra 算法应用示例以下是一个使用Dijkstra 算法求解最短路径的例题：在一个图中，有4 个节点和3 条边，如下所示：```1 --2 -- 3| /| /| /| /|/4```请问，节点1 到节点4 的最短路径是多少？。

最短路dijkstra算法Matlab程序

function [c0,c,path0,path]=dijkstra(s,t,C,flag)% Use the Dijkstra's algorithm to find the shortest path from% s to t and can also find the shortest path between s and all% the other points.% Reference: Graph Theory with Applications by J. A. Bondy and% U. S. R. Murty.% Input -- s is the starting point and also is the point s.% -- t is the given terminal point and is the point t.% -- C \in R^{n \times n}is the cost matrix, where% C(i,j)>=0 is the cost from point i to point j.% If there is no direct connection between point i and% j, C(i,j)=inf.% -- flag: if flag=1, the function just reports the% shortest path between s and t; if flag~=1, the% function reports the shortest path between s and t,% and the shortest paths between s and other points.% Output -- c0 is the minimal cost from s to t.% -- path0 denotes the shortest path form s to t.% -- c \in R{1\times n} in which the element i is the% minimal cost from s to point i.% -- path \in R^{n \times n} in which the row i denotes% the shortest path from s to point i.% Copyright by MingHua Xu(徐明华), Changhzou University, 27 Jan. 2014. s=floor(s);t=floor(t);n=size(C,1);if s<1 || t < 1 || s > n || t > nerror(' The starting point and the terminal point exceeds the valid range');endif t==sdisp('The starting point and the terminal point are the same points');endlabel=ones(1,n)*inf;label(s)=0;S=[s];Sbar=[1:s-1,s+1:n];c0=0;path=zeros(n,n);path(:,1)=s;c=ones(1,n)*inf;parent=zeros(1,n);i=1; % number of points in point set S.while i<n% for each point in Sbar, replace label(Sbar(j)) by% min(label(Sbar(j)),label(S(k))+C(S(k),Sbar(j)))for j=1:n-ifor k=1:iif label(Sbar(j)) > label(S(k))+C(S(k),Sbar(j))label(Sbar(j))=label(S(k))+C(S(k),Sbar(j));parent(Sbar(j))=S(k);endendend% Find the minmal label(j), j \in Sbar.temp=label(Sbar(1));son=1;for j=2:n-iif label(Sbar(j))< temptemp=label(Sbar(j));son=j;endend% update the point set S and SbarS=[S,Sbar(son)];Sbar=[Sbar(1:son-1),Sbar(son+1:n-i)];i=i+1;% if flag==1, just output the shortest path between s and t.if flag==1 && S(i)==tson=t;temp_path=[son];if son~=swhile parent(son)~=sson=parent(son);temp_path=[temp_path,son];endtemp_path=[temp_path,s];endtemp_path=fliplr(temp_path);m=size(temp_path,2);path0(1:m)=temp_path;c_temp=0;for j=1:m-1c_temp=c_temp+C(temp_path(j),temp_path(j+1));endc0=c_temp;path(t,1:m)=path0;c(t)=c0;returnendend% Form the output resultsfor i=1:nson=i;temp_path=[son];if son~=swhile parent(son)~=sson=parent(son);temp_path=[temp_path,son];endtemp_path=[temp_path,s];endtemp_path=fliplr(temp_path);m=size(temp_path,2);path(i,1:m)=temp_path;c_temp=0;for j=1:m-1c_temp=c_temp+C(temp_path(j),temp_path(j+1));endc(i)=c_temp;c0=c(t);path0=path(t,:);endreturn。

MATLAb最短路问题

称为 G 的由 E 1 导出的子图 , 记为 G [ E 1 ] .
[ { v 1 , v 4 , v 5 } ] [ { e 1 , e 2 , e 3 } ]
MATLAb最短路问题
返回
关联矩阵
对无向图Ｇ，其关联矩阵Ｍ＝ ( m i) j ，其中：
G 的图解如图.
MATLAb最短路问题
定义在图 G 中，与 V 中的有序偶 ( v i ， v j ) 对应的边 e ，称为图的有向
边（或弧），而与 V 中顶点的无序偶 v i v j相对应的边 e ，称为图的无向边 . 每一条边都是无向边的图，叫无向图；每一条边都是有向边的图，称为有向图；既有无向边又有有向边的图称为混合图 .
则称w(P) w(e)为路径P的权． eE(P)
(2) 在赋权图G中，从顶点u到顶点v的具有最小权的路
P*(u,v)，称为u到v的最短路．
MATLAb最短路问题
返回
固定起点的最短路
最短路是一条路径，且最短路的任一段也是最短路．假设在u0-v0的最短路中只取一条，则从u0到其余顶点的最短路将构成一棵以u0为根的树．
称为相邻的边．（４）边和它的端点称为互相关联的．（５）既没有环也没有平行边的图，称为简单图．（６）任意两顶点都相邻的简单图，称为完备图，记为Kn，其中n
为顶点的数目．
( 7)若V=X Y，X Y= ，X 中任两顶点不相邻，Y 中任两顶
点不相邻，称G为二元图；若X 中每一顶点皆与Y 中一切顶点相邻，称为完备二元图，记为Km,n，其中m,n 分别为X 与Y 的顶点数目．

Floyd最短路算法的MATLAB程序

Floyd最短路算法的MATLAB程序Floyd最短路算法的MATLAB程序%floyd.m%采用floyd算法计算图a中每对顶点最短路%d是矩离矩阵%r是路由矩阵function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendrfor 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)endendendkdrend数学算法（2）实例2：三角函数曲线（2）function shili02h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例02');x=-pi:0.05:pi;y=sin(x)+cos(x);plot(x,y,'-*r','linewidth',1); grid onxlabel('自变量X');ylabel('函数值Y');title('三角函数');实例3：图形的叠加function shili03h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例03');x=-pi:0.05:pi;y1=sin(x);y2=cos(x);plot(x,y1,...'-*r',...x,y2,...'--og');grid onxlabel('自变量X');ylabel('函数值Y');title('三角函数');实例4：双y轴图形的绘制function shili04h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例04');x=0:900;a=1000;b=0.005;y1=2*x;y2=cos(b*x);[haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot'); axes(haxes(1))ylabel('semilog plot');axes(haxes(2))ylabel('linear plot');实例6：图形标注function shili06h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例06');t=0:pi/10:2*pi;h=plot(t,sin(t));xlabel('t=0到2\pi','fontsize',16);ylabel('sin(t)','fontsize',16);title('\it{从 0to2\pi 的正弦曲线}','fontsize',16)x=get(h,'xdata');y=get(h,'ydata');imin=find(min(y)==y);imax=find(max(y)==y);text(x(imin),y(imin),...['\leftarrow最小值=',num2str(y(imin))],...'fontsize',16)text(x(imax),y(imax),...['\leftarrow最大值=',num2str(y(imax))],...'fontsize',16)实例7：条形图形function shili07h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例07');tiao1=[562 548 224 545 41 445 745 512];tiao2=[47 48 57 58 54 52 65 48];t=0:7;bar(t,tiao1)xlabel('X轴');ylabel('TIAO1值');h1=gca;h2=axes('position',get(h1,'position'));plot(t,tiao2,'linewidth',3)set(h2,'yaxislocation','right','color','none','xticklabel',[]) 实例8：区域图形function shili08h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例08');x=91:95;profits1=[88 75 84 93 77];profits2=[51 64 54 56 68];profits3=[42 54 34 25 24];profits4=[26 38 18 15 4];area(x,profits1,'facecolor',[0.5 0.9 0.6],...'edgecolor','b',...'linewidth',3)hold onarea(x,profits2,'facecolor',[0.9 0.85 0.7],... 'edgecolor','y',...'linewidth',3)hold onarea(x,profits3,'facecolor',[0.3 0.6 0.7],... 'edgecolor','r',...'linewidth',3)hold onarea(x,profits4,'facecolor',[0.6 0.5 0.9],... 'edgecolor','m',...'linewidth',3)hold offset(gca,'xtick',[91:95])set(gca,'layer','top')gtext('\leftarrow第一季度销量')gtext('\leftarrow第二季度销量')gtext('\leftarrow第三季度销量')gtext('\leftarrow第四季度销量')xlabel('年','fontsize',16);ylabel('销售量','fontsize',16);实例9：饼图的绘制function shili09h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例09');t=[54 21 35;68 54 35;45 25 12;48 68 45;68 54 69];x=sum(t);h=pie(x);textobjs=findobj(h,'type','text');str1=get(textobjs,{'string'});val1=get(textobjs,{'extent'});oldext=cat(1,val1{:});names={'商品一：';'商品二：';'商品三：'};str2=strcat(names,str1);set(textobjs,{'string'},str2)val2=get(textobjs,{'extent'});newext=cat(1,val2{:});offset=sign(oldext(:,1)).*(newext(:,3)-oldext(:,3))/2; pos=get(textobjs,{'position'});textpos=cat(1,pos{:});textpos(:,1)=textpos(:,1)+offset;set(textobjs,{'position'},num2cell(textpos,[3,2])实例10：阶梯图function shili10h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例10');a=0.01;b=0.5;t=0:10;f=exp(-a*t).*sin(b*t);stairs(t,f)hold onplot(t,f,':*')hold offglabel='函数e^{-(\alpha*t)}sin\beta*t的阶梯图';gtext(glabel,'fontsize',16)xlabel('t=0:10','fontsize',16)axis([0 10 -1.2 1.2])实例11：枝干图function shili11h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例11');x=0:pi/20:2*pi;y1=sin(x);y2=cos(x);h1=stem(x,y1+y2);hold onh2=plot(x,y1,'^r',x,y2,'*g');hold offh3=[h1(1);h2];legend(h3,'y1+y2','y1=sin(x)','y2=cos(x)') xlabel('自变量X');ylabel('函数值Y');title('正弦函数与余弦函数的线性组合');实例12：罗盘图function shili12h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例12');winddirection=[54 24 65 84256 12 235 62125 324 34 254];windpower=[2 5 5 36 8 12 76 14 10 8];rdirection=winddirection*pi/180;[x,y]=pol2cart(rdirection,windpower); compass(x,y);desc={'风向和风力','北京气象台','10月1日0：00到','10月1日12：00'};gtext(desc)实例13：轮廓图function shili13h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例13');[th,r]=meshgrid((0:10:360)*pi/180,0:0.05:1); [x,y]=pol2cart(th,r);z=x+i*y;f=(z.^4-1).^(0.25);contour(x,y,abs(f),20)axis equalxlabel('实部','fontsize',16);ylabel('虚部','fontsize',16);h=polar([0 2*pi],[0 1]);delete(h)hold oncontour(x,y,abs(f),20)实例14：交互式图形function shili14h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例14');axis([0 10 0 10]);hold onx=[];y=[];n=0;disp('单击鼠标左键点取需要的点'); disp('单击鼠标右键点取最后一个点'); but=1;while but==1[xi,yi,but]=ginput(1);plot(xi,yi,'bo')n=n+1;disp('单击鼠标左键点取下一个点'); x(n,1)=xi;y(n,1)=yi;endt=1:n;ts=1:0.1:n;xs=spline(t,x,ts);ys=spline(t,y,ts);plot(xs,ys,'r-');hold off实例15：变换的傅立叶函数曲线function shili15h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例15');axis equalm=moviein(20,gcf);set(gca,'nextplot','replacechildren') h=uicontrol('style','slider','position',... [100 10 500 20],'min',1,'max',20)for j=1:20plot(fft(eye(j+16)))set(h,'value',j)m(:,j)=getframe(gcf);endclf;axes('position',[0 0 1 1]);movie(m,30)实例16：劳伦兹非线形方程的无序活动function shili15h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例15');axis equalm=moviein(20,gcf);set(gca,'nextplot','replacechildren') h=uicontrol('style','slider','position',... [100 10 500 20],'min',1,'max',20)for j=1:20plot(fft(eye(j+16)))set(h,'value',j)m(:,j)=getframe(gcf);endclf;axes('position',[0 0 1 1]);movie(m,30)实例17：填充图function shili17h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例17');t=(1:2:15)*pi/8;x=sin(t);y=cos(t);fill(x,y,'r')axis square offtext(0,0,'STOP',...'color',[1 1 1],...'fontsize',50,...'horizontalalignment','center') 实例18：条形图和阶梯形图function shili18h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例18');subplot(2,2,1)x=-3:0.2:3;y=exp(-x.*x);bar(x,y)title('2-D Bar Chart') subplot(2,2,2)x=-3:0.2:3;y=exp(-x.*x);bar3(x,y,'r')title('3-D Bar Chart') subplot(2,2,3)x=-3:0.2:3;y=exp(-x.*x);stairs(x,y)title('Stair Chart')subplot(2,2,4)x=-3:0.2:3;y=exp(-x.*x);barh(x,y)title('Horizontal Bar Chart') 实例19：三维曲线图function shili19h0=figure('toolbar','none',... 'position',[200 150 450 400],... 'name','实例19');subplot(2,1,1)x=linspace(0,2*pi);y1=sin(x);y2=cos(x);y3=sin(x)+cos(x);z1=zeros(size(x));z2=0.5*z1;z3=z1;plot3(x,y1,z1,x,y2,z2,x,y3,z3) grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure1:3-D Plot') subplot(2,1,2)x=linspace(0,2*pi);y1=sin(x);y3=sin(x)+cos(x);z1=zeros(size(x));z2=0.5*z1;z3=z1;plot3(x,z1,y1,x,z2,y2,x,z3,y3)grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure2:3-D Plot')实例21：PEAKS函数曲线function shili21h0=figure('toolbar','none',...'position',[200 100 450 450],...'name','实例21');[x,y,z]=peaks(30);subplot(2,1,1)x=x(1,:);y=y(:,1);i=find(y>0.8&y<1.2);j=find(x>-0.6&x<0.5);z(i,j)=nan*z(i,j);surfc(x,y,z)xlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure1:surfc函数形成的曲面') subplot(2,1,2)x=x(1,:);i=find(y>0.8&y<1.2);j=find(x>-0.6&x<0.5);z(i,j)=nan*z(i,j);surfl(x,y,z)xlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure2:surfl函数形成的曲面')实例22：片状图function shili22h0=figure('toolbar','none',...'position',[200 150 550 350],...'name','实例22');subplot(1,2,1)x=rand(1,20);y=rand(1,20);z=peaks(x,y*pi);t=delaunay(x,y);trimesh(t,x,y,z)hidden offtitle('Figure1:Triangular Surface Plot'); subplot(1,2,2)x=rand(1,20);y=rand(1,20);z=peaks(x,y*pi);t=delaunay(x,y);trisurf(t,x,y,z)title('Figure1:Triangular Surface Plot'); 实例23：视角的调整function shili23h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例23');x=-5:0.5:5;[x,y]=meshgrid(x);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;subplot(2,2,1)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure1')view(-37.5,30)subplot(2,2,2)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure2')view(-37.5+90,30)subplot(2,2,3)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure3')view(-37.5,60)subplot(2,2,4)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure4')view(180,0)实例24：向量场的绘制function shili24h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例24');subplot(2,2,1)z=peaks;ribbon(z)title('Figure1')subplot(2,2,2)[x,y,z]=peaks(15);[dx,dy]=gradient(z,0.5,0.5); contour(x,y,z,10)hold onquiver(x,y,dx,dy)hold offtitle('Figure2')subplot(2,2,3)[x,y,z]=peaks(15);[nx,ny,nz]=surfnorm(x,y,z); surf(x,y,z)hold onquiver3(x,y,z,nx,ny,nz)hold offtitle('Figure3')subplot(2,2,4)x=rand(3,5);y=rand(3,5);z=rand(3,5);c=rand(3,5);fill3(x,y,z,c)grid ontitle('Figure4')实例26：柱状图function shili26h0=figure('toolbar','none',... 'position',[200 50 450 450],... 'name','实例26');subplot(2,1,1)x=[5 2 18 7 39 8 65 5 54 3 2];bar(x)xlabel('X轴');ylabel('Y轴');title('第一子图');subplot(2,1,2)y=[5 2 18 7 39 8 65 5 54 3 2];xlabel('X轴');ylabel('Y轴');title('第二子图');实例28：羽状图function shili28h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例28');subplot(2,1,1)alpha=90:-10:0;r=ones(size(alpha));m=alpha*pi/180;n=r*10;[u,v]=pol2cart(m,n);feather(u,v)title('羽状图')axis([0 20 0 10])subplot(2,1,2)t=0:0.5:10;x=0.05+i;y=exp(-x*t);feather(y)title('复数矩阵的羽状图')实例29：立体透视（1）function shili29h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例29');[x,y,z]=meshgrid(-2:0.1:2,...-2:0.1:2);v=x.*exp(-x.^2-y.^2-z.^2); grid onfor i=-2:0.5:2;h1=surf(linspace(-2,2,20),... linspace(-2,2,20),...zeros(20)+i);rotate(h1,[1 -1 1],30)dx=get(h1,'xdata');dy=get(h1,'ydata');dz=get(h1,'zdata');delete(h1)slice(x,y,z,v,[-2 2],2,-2)hold onslice(x,y,z,v,dx,dy,dz)hold offaxis tightview(-5,10)drawnowend实例30：立体透视（2）function shili30h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例30');[x,y,z]=meshgrid(-2:0.1:2,...-2:0.1:2,...-2:0.1:2);v=x.*exp(-x.^2-y.^2-z.^2);[dx,dy,dz]=cylinder;slice(x,y,z,v,[-2 2],2,-2)for i=-2:0.2:2h=surface(dx+i,dy,dz);rotate(h,[1 0 0],90)xp=get(h,'xdata');yp=get(h,'ydata');zp=get(h,'zdata');delete(h)hold onhs=slice(x,y,z,v,xp,yp,zp);axis tightxlim([-3 3])view(-10,35)drawnowdelete(hs)hold offend实例31：表面图形function shili31h0=figure('toolbar','none',...'position',[200 150 550 250],... 'name','实例31');subplot(1,2,1)x=rand(100,1)*16-8;y=rand(100,1)*16-8;r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;xlin=linspace(min(x),max(x),33); ylin=linspace(min(y),max(y),33);[X,Y]=meshgrid(xlin,ylin);Z=griddata(x,y,z,X,Y,'cubic'); mesh(X,Y,Z)axis tighthold onplot3(x,y,z,'.','Markersize',20) subplot(1,2,2)k=5;n=2^k-1;theta=pi*(-n:2:n)/n;phi=(pi/2)*(-n:2:n)'/n;X=cos(phi)*cos(theta);Y=cos(phi)*sin(theta);Z=sin(phi)*ones(size(theta)); colormap([0 0 0;1 1 1])C=hadamard(2^k);surf(X,Y,Z,C)axis square实例33：曲线转换按钮h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例33');x=0:0.5:2*pi;y=sin(x);h=plot(x,y);grid onhuidiao=[...'if i==1,',...'i=0;,',...'y=cos(x);,',...'set(hm,''string'',''正弦函数''),',...'h=plot(x,y);,',...'grid on,',...'else if i==0,',...'i=1;,',...'y=sin(x);,',...'set(hm,''string'',''余弦函数''),',...'delete(h),',...'h=plot(x,y);,',...'grid on,',...'end,',...'end'];hm=uicontrol(gcf,'style','pushbutton',... 'string','余弦函数',...'callback',huidiao);i=1;set(hm,'position',[250 20 60 20]);set(gca,'position',[0.2 0.2 0.6 0.6])title('按钮的使用')hold on实例34：栅格控制按钮h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例34');x=0:0.5:2*pi;y=sin(x);plot(x,y)huidiao1=[...'set(h_toggle2,''value'',0),',...];huidiao2=[...'set(h_toggle1,''value'',0),',...'grid off,',...];h_toggle1=uicontrol(gcf,'style','togglebutton',... 'string','grid on',...'value',0,...'position',[20 45 50 20],...'callback',huidiao1);h_toggle2=uicontrol(gcf,'style','togglebutton',... 'string','grid off',...'value',0,...'position',[20 20 50 20],...'callback',huidiao2);set(gca,'position',[0.2 0.2 0.6 0.6])title('开关按钮的使用')实例35：编辑框的使用h0=figure('toolbar','none',...'position',[200 150 350 250],...'name','实例35');f='Please input the letter';huidiao1=[...'g=upper(f);,',...'set(h2_edit,''string'',g),',...];huidiao2=[...'g=lower(f);,',...'set(h2_edit,''string'',g),',...];h1_edit=uicontrol(gcf,'style','edit',...'position',[100 200 100 50],...'HorizontalAlignment','left',...'string','Please input the letter',...'callback','f=get(h1_edit,''string'');',...'background','w',...'max',5,...'min',1);h2_edit=uicontrol(gcf,'style','edit',...'HorizontalAlignment','left',...'position',[100 100 100 50],...'background','w',...'max',5,...'min',1);h1_button=uicontrol(gcf,'style','pushbutton',... 'string','小写变大写',...'position',[100 45 100 20],...'callback',huidiao1);h2_button=uicontrol(gcf,'style','pushbutton',... 'string','大写变小写',...'position',[100 20 100 20],...'callback',huidiao2);实例36：弹出式菜单h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例36');x=0:0.5:2*pi;y=sin(x);h=plot(x,y);grid onhm=uicontrol(gcf,'style','popupmenu',... 'string',...'sin(x)|cos(x)|sin(x)+cos(x)|exp(-sin(x))',... 'position',[250 20 50 20]);set(hm,'value',1)huidiao=[...'v=get(hm,''value'');,',...'switch v,',...'case 1,',...'delete(h),',...'y=sin(x);,',...'h=plot(x,y);,',...'grid on,',...'case 2,',...'delete(h),',...'y=cos(x);,',...'h=plot(x,y);,',...'grid on,',...'case 3,',...'delete(h),',...'y=sin(x)+cos(x);,',...'h=plot(x,y);,',...'grid on,',...'case 4,',...'delete(h),',...'y=exp(-sin(x));,',...'h=plot(x,y);,',...'grid on,',...'end'];set(hm,'callback',huidiao)set(gca,'position',[0.2 0.2 0.6 0.6]) title('弹出式菜单的使用')hold on实例37：滑标的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例37');[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;h0=mesh(x,y,z);h1=axes('position',...[0.2 0.2 0.5 0.5],...'visible','off');htext=uicontrol(gcf,...'units','points',...'position',[20 30 45 15],...'string','brightness',...'style','text');hslider=uicontrol(gcf,...'units','points',...'position',[10 10 300 15],...'min',-1,...'max',1,...'style','slider',...'callback',...'brighten(get(hslider,''value''))'); 实例38：多选菜单h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例38');[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;h0=mesh(x,y,z);hlist=uicontrol(gcf,'style','listbox',...'string','default|spring|summer|autumn|winter',... 'max',5,...'min',1,...'position',[20 20 80 100],...'callback',[...'k=get(hlist,''value'');,',...'switch k,',...'case 1,',...'colormap default,',...'case 2,',...'colormap spring,',...'case 3,',...'colormap summer,',...'case 4,',...'colormap autumn,',...'case 5,',...'colormap winter,',...'end']);实例39：菜单控制的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例39');x=0:0.5:2*pi;y=cos(x);h=plot(x,y);grid onset(gcf,'toolbar','none')hm=uimenu('label','example');huidiao1=[...'set(hm_gridon,''checked'',''on''),',...'set(hm_gridoff,''checked'',''off''),',...'grid on'];huidiao2=[...'set(hm_gridoff,''checked'',''on''),',...'set(hm_gridon,''checked'',''off''),',...'grid off'];hm_gridon=uimenu(hm,'label','grid on',... 'checked','on',...'callback',huidiao1);hm_gridoff=uimenu(hm,'label','grid off',... 'checked','off',...'callback',huidiao2);实例40：UIMENU菜单的应用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例40');h1=uimenu(gcf,'label','函数');h11=uimenu(h1,'label','轮廓图',...'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'[x,y,z]=peaks;,',...'contour3(x,y,z,30)']);h12=uimenu(h1,'label','高斯分布',... 'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'mesh(peaks);,',...'axis tight']);h13=uimenu(h1,'label','Sinc函数',... 'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'[x,y]=meshgrid(-8:0.5:8);,',...'r=sqrt(x.^2+y.^2)+eps;,',...'z=sin(r)./r;,',...'mesh(x,y,z)']);h2=uimenu(gcf,'label','色彩');hl2(1)=uimenu(h2,'label','Default',... 'checked','on',...'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(1),''checked'',''on''),',...'colormap(''default'')']);hl2(2)=uimenu(h2,'label','spring',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(2),''checked'',''on''),',...'colormap(spring)']);hl2(3)=uimenu(h2,'label','Summer',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(3),''checked'',''on''),',...'colormap(summer)']);hl2(4)=uimenu(h2,'label','Autumn',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(4),''checked'',''on''),',...'colormap(autumn)']);hl2(5)=uimenu(h2,'label','Winter',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(5),''checked'',''on''),',...'colormap(winter)']);h3=uimenu(gcf,'label','坐标选项');h31=uimenu(h3,'label','Axis on',... 'callback',...[...'axis on,',...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off'')']);h32=uimenu(h3,'label','Axis off',... 'callback',...[...'axis off,',...'set(h32,''checked'',''on''),',...'set(h31,''checked'',''off'')']);实例41：除法计算器h=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例41');h1=uicontrol(gcf,'style','edit',...'position',[80 200 100 20],...'HorizontalAlignment','right',...'callback',['m=get(h1,''string'');,',...'a=str2num(m);']);h2=uicontrol(gcf,'style','edit',...'HorizontalAlignment','right',...'position',[80 150 100 20],...'callback',['n=get(h2,''string'');,',...'b=str2num(n);']);h3=uicontrol(gcf,'style','text',...'string','被除数',...'position',[80 230 100 20]);h4=uicontrol(gcf,'style','edit',...'position',[80 50 100 20]);h5=uicontrol(gcf,'style','pushbutton',...'position',[80 100 100 20],...'string','=',...'callback',[...'if b==0,',...'h7=errordlg(''除数不能为0！'',''error'',''on'');,',... 'else,',...'k=a/b;,',...'c=num2str(k);,',...'set(h4,''string'',c),',...'end']);h8=uicontrol(gcf,'style','text',...'string','除数',...'position',[80 175 100 20]);h9=uicontrol(gcf,'style','text',...'string','商',...'position',[80 75 100 20]);实例42：单选框的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例42');x=0:0.5:2*pi;y=sin(x);plot(x,y)grid onset(gcf,'toolbar','none')g=set(gca,'position',[0.2 0.2 0.6 0.6]); huidiao1=[...'grid on,',...'set(box_on,''value'',1),',...'set(box_off,''value'',0),'];huidiao2=[...'grid off,',...'set(box_off,''value'',1),',...'set(box_on,''value'',0),'];box_on=uicontrol(gcf,'style','radio',... 'position',[5 50 50 20],...'string','grid on',...'value',1,...'callback',huidiao1);box_off=uicontrol(gcf,'style','radio',... 'position',[5 20 50 20],...。

电力系统短路故障分析的MATLAB辅助程序设计短路计算程序

电力系统短路故障分析的MATLAB辅助程序设计电力系统短路故障可分为三相对称短路故障(three-phase balanced faults)和不对称短路故障(unbalanced faults )。

不对称短路故障又分为单相接地短路故障(single line-to-ground fault)、两相短路故障 (line-to-line fault)以及两相接地短路故障(double line-to-ground fault)。

根据故障分析结果可以对继电保护装置、自动装置进行整定计算，我们可以建立算法来形成节点阻抗矩阵,利用节点阻抗矩阵来计算短路故障情况下的节点电压和线路电流。

一、三相对称短路故障进行三相短路计算需要两个程序：zbuild /zbuildpi和symfault程序，zbuild、zbuildpi用来在MATLAB中形成节点阻抗矩阵，symfault用来计算三相对称故障。

Zbus=zbuild(zdata)这里的参数zdata是一个(e×4)阶矩阵，e是拓扑图的总支路数目。

第一列和第二列为元素两端的节点编号，第三列和第四列分别是线路的电阻、电抗的标幺值。

连接在0节点和发电机节点之间的发电机阻抗可能是次暂态电抗、暂态电抗或同步电抗，而且这个矩阵中还包含并联电抗器和负荷阻抗。

Zbus=zbuildpi(linedata,gendata,yload)这个函数与潮流计算程序是相容的，第一个参数linedat a与潮流计算程序中的文件是一致的。

第一列和第二列为节点编号；第三列到第五列分别是线路的电阻、电抗以及1/2线路电纳值，这三项都为在统一基准容量下的标幺值；最后一列是变压器分接头位置，对线路来说，必须输入1；线路无输入顺序。

发电机参数不包含在Linedata参数中，而是包含在第二个参数gendata中，gendata是一个g×4阶矩阵,g是发电机总数。

第一列和第二列为0节点、发电机节点编号，第三列和第四列为发电机的暂态电阻和暂态电抗。

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程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。

k短路算法代码matlab

k短路算法代码matlab"k短路"算法可能指的是Kruskal算法，这是一种用于找出最小生成树的算法。

以下是使用MATLAB实现的Kruskal 算法：matlab复制代码:function [tree, weight] = kruskal(graph)% graph是一个邻接矩阵，表示图的结构% 图的顶点被编号为1到n，其中n是顶点的数量% 图的边由(u, v, w)三元组表示，其中u和v是顶点，w是边的权重% 输入:% graph - 邻接矩阵% 输出:% tree - 最小生成树的边% weight - 最小生成树的权重n = size(graph, 1); % 获取顶点数量edges = find(~eye(n) & graph); % 获取所有边edges = sortrows(edges); % 对边进行排序tree = []; % 初始化最小生成树weight = 0; % 初始化最小生成树的权重while numel(tree) < n - 1u = edges(1, 1);v = edges(1, 2);w = edges(1, 3);edges(1, :) = []; % 从边列表中移除已经使用的边if ~ismember(v, [tree(:,1) tree(:,2)], 'binary') && ~ismember(u, [tree(:,1) tree(:,2)], 'binary')tree = [tree; u v w]; % 将边添加到最小生成树中weight = weight + w; % 更新最小生成树的权重endendend注意：这个函数假设图是连通的，也就是说，从任意一个顶点都可以到达其他所有顶点。

如果图不是连通的，那么这个函数将返回一个子图的最小生成树，而不是整个图的最小生成树。

Floyd最短路算法的MATLAB程序

Floyd最短路算法的MATLAB程序2006-08-17 20:14%floyd.m%采用floyd算法计算图a中每对顶点最短路 %d是矩离矩阵%r是路由矩阵function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendrfor 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)endendendkdrendvoid Dijkstral(int v0){int i;bool s[MAX_VEX];for(i=0;i<dim;i++){d[v0][i]=map[v0][i];s[i]=false;if((i!=0)&&(d[v0][i]<INF))p[v0][i]=v0;elsep[v0][i]=-1;}s[v0]=true;d[v0][v0]=0;for(i=0;i<dim;i++){double min=INF;int u=v0;for(int j=0;j<dim;j++)if(!s[j]&&d[v0][j]<min){u=j;min=d[v0][j];}s[u]=true;for(int w=0;w<dim;w++){if((!s[w])&&(d[v0][w]>d[v0][u]+map[u][w])){d[v0][w]=d[v0][u]+map[u][w];p[v0][w]=u;}}}}Justin Hou介绍寻找最有价值路径（c语言）描述：从上（入口）往下行走，直到最下节点（出口）结束，将所经节点上的数值相加，要求找到一条最有价值路径（既是路径总数值最大）并输出总数值。

图：入口↓③/\⑤④/ \ / \①②⑤\ / \ /⑨⑧\ /③出口↓输入文件：（abc.in)第一行只有一个数n(1<=n<=199),且n为奇数，说明节点的层。

(完整版)短路计算matlab程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。

matlab图论程序算法大全

精心整理图论算法matlab实现求最小费用最大流算法的 MATLAB 程序代码如下：n=5;C=[0 15 16 0 00 0 0 13 14forwhileforfor(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 的最短路, 算法终止. 注意在求最小费用最大流时构造有whileifelseifififpd=0;值t=n;ifelseifif(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 %显示最小费用, 程序结束__Kruskal 避圈法：k=1; %forifkk=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;whileforforifif(pd)pd=0;forifelseT %用2 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); %k %D %R %pd=0;ifif(pd)end %0 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)whileforfor弧时ifelseifififif(pd)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 %显示标号, 由此可得最小割, 程序结束ifforendelseifforW(a(1),a(2))=1;W(a(2),a(1))=1;endelsefprint('Please imput the right value of f');endW;程序二：可达矩阵算法function P=dgraf(A)精心整理n=size(A,1);P=A;for i=2:nP=P+A^i;endP(P~=0)=1;P;程序三：有向图关联矩阵和邻接矩阵互换算法function W=mattransf(F,f)if f==0forendelseifforendelsefprint('Please imput the right value of f');endW;第二讲：最短路问题程序一：Dijkstra算法（计算两点间的最短路）function [l,z]=Dijkstra(W)l(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;endendendendendend程序三：n2short.m 计算指定两点间的最短距离function [P u]=n2short(W,k1,k2)n=length(W);U=W;m=1;for 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);k=1;kk=k2;whileforendendk=1;forendP;)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;elseendP;d;ifforendendelseb=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;endc;k=1:l;e=1;whileformin=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 算法构造最小生成树）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;endendresult第四讲：Euler图和Hamilton图程序一：Fleury算法（在一个Euler图中找出Euler环游）注：包括三个文件；fleuf1.m, edf.m, flecvexf.mfunction [T c]=fleuf1(d)%注：必须保证是Euler环游，否则输出T=0,c=0b(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:nifendendif m~=0endif m==0forfor g=1:edsc=c+d(T(1,g),T(2,g));endflagg=0;break;endendendendfunction[flag ed]=edf(matr,eds,vexs,ed,tem)[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;endendf=0;dvex=0;ded=[];ifelseforendif 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;endend精心整理if 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);elseendend路）%d%Cendendendd2=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):jd1;endfor 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;endifendendJ;ifendb=a;ifendfor 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');whileendwhiletp=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))');for%求初始匹配Mif(M(i,j))break;end;end %获得仅含一条边的初始匹配M while(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 表示0for取Xif标记for;e nd %if(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 M-增广路forM去掉Mif*,M %程序3 Kuhn-Munkres算法（即赋权完备二元图的最佳匹配程序）function kk=kexingdian(A)%求赋权完备二元图最佳匹配A=[4 5 5 1;2 2 4 6;4 2 3 3;5 0 2 1]; %A为邻接矩阵n=length(A);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; %初始可行点标记L M(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;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;end if(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;endfor(j=1:jsn)pd=1; %取y∈NL(S)\Tfor(k=1:jst)if(T(k)==NlS(j))pd=0;break;end;end if(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 %%C%f1%f%wfifelsef=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);endifendendt=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);whileforendforendforendforfor (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;end精心整理endif (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));endendifendendif (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);endend。

MATLAB编程：最短路问题

（ 2 ）更新 l ( v ) 、 z ( v ) ： v S V \ S ,若 l ( v ) > l ( u ) W ( u , v )
z 则令 l(v ) = l(u ) W (u , v ) , (v ) = u
（ 3）设 v 是使 l(v ) 取最小值的 S
定义３（１）设 P (u ,v)是赋权图 G 中从 u 到 v 的路径，则称 w(P)
e E ( P )
w (e) 为路径
P 的权．
(2 )
在赋权图 G 中，从顶点 u 到顶点 v 的具有最小权的路
P (u , v ) ，称为 u 到 v 的最短路．
u2
u 6
6
u5
图 G 的边为边集的图 G 的子图，称为 G 的由 V 1 导出的子图，记为 G[V 1 ]. (3)设 E 1 E ,且 E 1 ,以 E 1 为边集 ,E 1 的端点集为顶点集的图 G 的子图 , 称为 G 的由 E 1 导出的子图 ,记为 G[E 1 ].
返回
邻接矩阵
对无向图Ｇ，其邻接矩阵 A ( a ij ) ，其中：
a ij
1 0
若 v i 与 v j 相邻若 v i 与 v j 不相邻
v1 A= 0 1 0 1 v2 1 0 1 1 0 1 0 1
注：假设图为简单图
返回
顶点的次数
定义（１）在无向图中，与顶点 v 关联的边的数目（环算两次）称为 v 的次数，记为 d (v)．（２）在有向图中，从顶点 v 引出的边的数目称为 v 的出度，记为 d + ( v)，从顶点 v 引入的边的数目称为的入度，记为 d - (v)， d ( v)= d + ( v)+ d - ( v) 称为 v 的次数．