MATLAB基础
一、向量、矩阵的表示和使用
format long %小数很多
format short %默认4位小数
format rat %显示最近的分数
format short e %指数格式的数 尾数多少exp(1) %自然对数的底
exp(10)
log(x) %x的自然对数 底是e
log10(x) %以10为底asin(pi)
atan(pi/3) %反三角函数向量（vector）一维数值数组。MATLAB 允许你创建列向量和行向量，列向量通过在方
括号内把数值用分号（;）隔开来创建，对元素的个数没有限制。linspace(a, b)创建了 a、b 之间含有 100 个等差元素的向量，
linspace(a, b, n)创建了 a、b 之间含有 n 个等差元素的向量。
logspace(a, b, n)创建 n 个对数值相隔相同的行向量 这创建了 10^a 和 10^b 之间 n 个对数值等差的向量。(对应元素相乘得到的向量)：A=[1 2 3] A.*A
求向量A的模： sqrt(sum(A.*A))abs 命令返回向量的绝对值 即在原位置返回向量中每个元素的绝对
值max(A):A为矩阵，返回列向量最大值
cumsum(A):A为矩阵，列向量累和运算向量点乘：逐个数相乘相加得到一个数：
dot(a,b) / sum(a.*b)
对于带有复数元素的向量，dot 操作也能正确计算：dot 可以用来求模，只要再 sqrt 一下 复数也可以
cross 叉乘：得到的依旧是向量
A(4:6)选出第 4 到第 6 个元素组成一个新的、含有 3 个元素的向量
数组乘法【相应元素相乘】：A.*B A和B的shape完全相同
eye(n)：n维单位矩阵
引用矩阵或者向量的元素要用【圆括号】
对于矩阵A
要引用第 i 列的所有元素，我们输入 A(:,i)
要选出从第 i 列到第 j 列之间的所有元素，我们输入 A(:,i:j)。通过引用矩阵中的行或列来创建新的矩阵:
我们复制 A 矩阵的第一行四次来创建一个新矩阵：
>> E = A([1,1,1,1],:)
复制 A 矩阵的第3,2,2,2,5行形成新矩阵：
>> E = A([3,2,2,2,5],:)求行列式：det(A) A方阵
行列式可以用来找出一个线性系统方程是否有解求线性方程组的解：
5x + 2y - 9z = -18
-9x - 2y + 2z = -7
6x + 7y + 3z = 29A = [5,2,-9;-9,2,2;6,7,3];
b = [-18 -7 29]';
A\b %不是点乘，因为点乘是逐个元素运算的意思
相当于inv(A)*b
如果方程欠定但解存在，则A\b返回：通过把其中一个变量（本例中是 z）设为零产生一组解
解存在：rank(A) = rank([A b])通过伪逆矩阵解欠定方程组：
x = pinv(A) * brank 求方阵的秩
求逆矩阵：inv(A) A必须为方阵
pinv(A):伪逆矩阵rref(A) :A是系数矩阵，产生A用guass消元法的梯形矩阵
magic(n):产生n阶魔方矩阵
魔方矩阵（幻方）是一个 n×n矩阵。矩阵的元素从 1 到 n^2 之间，并且行元素的和等于列元素的和.对矩阵进行 LU、QR 或
SVD 分解也是可行的
LU分解：
[L,U] = lu(A) ->求解：x = U\(L\b)
二、绘图
当一个函数是由二个或更多个函数相乘构成，别忘记在相乘时加上“.”以便告诉 MATLAB 我们是对两个矩阵进行相乘。
基本绘图：
plot(x,y),title(''),xlabel(),ylabel();
fplot(@(x)function(x),[start,end])
grid on 添加网格
在绘图命令行中加进 axis square，这会使得 MATLAB 产生正方形图象。
如果我们输入 axis equal，那么 MATLAB会产生一个两坐标轴比例和间距都相同的图象。
axis auto：自动选择
plot的调整：
'r--' 线条红色且是虚线
xlabel : 设置x坐标
ylabel : 设置y坐标
title : 设置标题
legend('A','B') :为每一条曲线设置图例，可自由拖动
axis([-5,5,0,1]) ：设置坐标轴比例，x是-5到5，y是-1到1
grid on 添加网格
hold on 保持绘图窗口，后续继续向同一图内画线
实线 '-'
虚线 '--'
虚点线 '-.'
点线 ':'
白色 w
黑色 k
蓝色 b
红色 r
青色 c
绿色 g
洋红 m
蓝色 y创建子图
subplot(1,2,1),plot,hold on：1行2列，当前绘图窗口在1
subplot(1,2,2),plot
极坐标绘图：polar(x,y)
log坐标绘图：loglog(x,y)
柱状图：bar(x,y)
针状图：stem(x,y,'--dg','fill'):'--dg':依次是线的形状，针端点的类型，和颜色；'fill'表示填充
包括方块(s)、菱形(d)、五角星(p)、圆圈(o)、叉号(x)、星号(*)和点号(.)用来产生 x 数集的命令，即 linspace 命令
x = linspace(a,b) a到b之间100个点
x = linspace(a,b,n) n个meshgrid 是一个可以为我们建立独立变量的一个易用的函数，
它所做的工作是为我们产生矩阵元素，元素x和y按照我们分别所指定的范围和增量来产生。
如[x,y] = meshgrid(-5:0.1:5,2:0.1:3)，产生矩阵，x-5到5y2到3的组合，间隔0.1等高线图：
先用meshgrid产生x，y数据，他是x-o-y平面上的每一对点
[x,y] = meshgrid(1:0.1:2,1:0.1:3)
接下来定义高程函数z；比如z = sin(x).*cos(x)
可以使用[C,h] = contour(x,y,z)，xlabel,ylabel 产生等高线图，但这是二维的
进一步操作标记等高线(不关闭图形)，上面C是数据，h是该图形的对象
使用：set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2) 标记等高线
'ShowText'标记，'on'等高线内部，get命令是每条线一标记，上面是get(h,'LevelStep')*2是每隔一条线一标记三维等高线图：contour3(x,y,z)
美化三维等高线图:设置EdgeColor、FaceColor以及底面网格
surface(x,y,z,'EdgeColor',[.8 .8 .8],'FaceColor','none')
grid off
view(-15,10)散点图：scatter
cylinder 函数：绘制三维圆柱图
[x,y,z]=cylinder 函数返回一半径和高度都为1的圆柱体x，y，z轴的坐标值，圆柱体沿其周长有20个等距分布的点
[x,y,z]=cylinder(r) 函数一个半径为r、高度为1的圆柱体的x，y，z轴的坐标值，圆柱体沿其周长有20个等距分布的点
[x,y,z]=cylinder(r,n) 函数一个半径为r、高度为1的圆柱体的x，y，z轴的坐标值，圆柱体沿其周长有n个等距分布的点绘制三维图像：
上面surface可以美化三维图
绘制三维图像：mesh 或者 surf 系列的命令
mesh(x,y,z) :表面没有颜色
surf(x,y,z) : 表面渐变的颜色
surfc(x,y,z): 会在图象中留下映像
surfl(x,y,z): l告诉我们这是一个光照表面（lighted surface））是另一个好的选择，它给了我们显示三维光照物体的表面。shading 命令,surf 后跟上
设置图中的阴影。shading interp/shading flat/shading faceted
flat 是用同一颜色为每个网格进行着色并隐藏网格线，而 facted 则显示网格，使用 interp 是
告诉 MATLAB 使用颜色插值的办法进行着色，因此显得非常平滑colormap(gray) ：设置为灰度图
axis square：XY坐标面正方形三维图像总结：
plot3 三维曲线图；
mesh 三维网格图；
meshc 除了生成网格图外，还在xy平面生成曲面的等高线；
meshz 除了生成网格图外，还在曲线下面加上个矩形垂帘；
surf 三维着色曲面图；
surfc 同时画出三维着色曲面图与等高线；
surfl 带光照的三维着色曲面图。
三、统计和编程基础
柱状图命令bar
bar(x,y,'grouped') %默认，组合柱状图
bar(x,y,'stacked') %堆积柱状图
barh %横向
bar3 %三维
bar3h %横向三维多组数据：如在一个图中统计三个班的分数
bar(x,y)
其中 x依旧是横轴，但y是矩阵，矩阵的每一列是一个班级的数据统计：
mean：平均
std：标准差
median：中位数计算加权平均，如给定的数据70分有3人，80分2人等 需要自己写个程序
x = [70,80];
n = [3,2];
me = sum(x.*n)/sum(n)编程基础：
for循环
for index = start:increment:end
statements
end
if-else-end:
if condition
body
elseif condition
body
else
body
end
while循环：
while condition
body
end
switch-case-end
switch switch_expression
case case_expression
body
case case_expression2
body
otherwise
body
end
Other:
输入：x = input('请输入x：');
输出：disp(x)
disp('要显示的信息')
四、代数方程、方程组的求解和简化
solve 函数
匿名函数：solve(@(x) x+2==0);%注意是双等号
f = 'x+2=0'; %其中可以带常量a，可以是任意可求解的方程
x = solve(f);绘图方程的绘制：
f ='...';%注意，绘制f仅仅是这个方程，不带=0
ezplot(f);
比如：
f = 'sin(x)'
ezplot(f)
ezplot(f,[start,end]); %指定绘图区间
%但 f = 'sin(x) =0';ezplot(f)不可行使用符号变量
syms x;
f = x^2 +x -1;
solv(x); %这种求解方法更加推荐有时求得解很麻烦，直接double强制转换一下
方程组求解
solve(f,g); %fg都是方程，如：
syms x,y;
f = x^2+y - 2
g = x+y-1
s = solve(x,y)
s.x
s.y %使用s作为返回变量，s.x返回x的求解值
注意不可直接solve不返回x直接double强制转换syms x
方程展开：expand(f)
方程合并：collect(f)
泰勒展开：taylor(f)
五、微积分相关
符号微分方程求解：dsolve 函数
dsolve('eqn', 'cond1', 'cond2', ...)。
统统用 '' 括起来
'' 中的方程是微分方程，导数用D指示，二阶导数用D2指示E.g.
dsolve('Df = -2*f + cos(t)')
dsolve('D2f + f = 4*sin(t)','f(0) = 0','Df(0) = 1')dsolve默认用t作为独立变量。可以指定独立变量：最后一项参数'x'
g = dsolve('D2f - sin(x) = x^3','f(0) = 2','x')dsolve 求解方程组：只需要依次将各个方程作为参数
常微分方程:
dsolve('Dy = a*y') %返回带常数C的通解计算符号导数导数 diff
syms x
f = x^2;
g = sin(10*t);
diff(f) %求得2*x
diff(g)
高阶导数：
求n阶导数：diff(f,n)令表达式好看一些：
pretty(f)两变量值是否相等：isequal
求极限
limit(f,a) x趋向于a处的极限
E.g.limit(x+5,3)
左/右极限：limit(f,3,'left/right')
limit命令还可以得到渐近线->用到再查求符号方程在某点c处的值： subs(f,c)
syms x
f = x^2+2*x+1
y = subs(f,3) %求f在x=3处的值
%另外：
方程有多个未知变量时，可以用subs指定带入的变量
f2 = subs(f,'x',2)
E.g.
syms x y
f = x^2+y
f2 = subs(f,'x',3)在一个图中曲线上标记处特定值点：
plot(x1,subs(f,x1),'o',x,f(x))
在一个图中添加文本text说明：
text(0.8,3.2,'文本内容') %0.8 3.3是文本起始点坐标位置常微分的数值解用 ode23 ode45 命令求解，在此略
积分命令：int
不定积分：int('f') %其中f是用''括起来的
指定积分变量(f多变量时),int('f',x)
多重积分：int(int('f'))
定积分： int('f',start,end)
E.g.
int('3*y^2+sec(x)',x)
梯形公式积分：trapz(x,y) %x-y是对应的数据点对
正交积分：
quad('f',start,end) %采用辛普森积分
quadl('f'.start,end) %采用洛巴托积分拟合命令：
ployfit(x,y,n)
记得计算r^2 r^2->1则拟合效果好MATLAB集成了许多特殊函数，如伽马函数，gamma()，勒让德函数等
六、图论
最短路:
Dijkstra算法：
C实现：（使用说明：假设n个点m条边，输入m条边的信息，即入点、出点、权值，调用计算函数，得出start点到每一个点的最短距离，如果是规划路径，要再加一些操作）
#include<iostream>
#include<stdio.h>
#include<queue>
#include<algorithm>
#define INF 100000
#define MAX 300
using namespace std;
int N,M;
int d[MAX]; //用来存放最短距离
typedef pair<int,int> P ; //用于优先队列中距离与顶点的对应，其中first为距离 second为编号
int path[MAX]; //记录路径
struct edge //边的结构体
{
int from; //从这个点
int to; //到这个点 的一条边
int cost; //权值
edge(int t, int c):to(t), cost(c) {} //构造函数
};
void Dijkstra(int s,vector<edge> g[MAX]){
priority_queue<P, vector<P>, greater<P> > que; //优先队列
fill(d, d+MAX, INF);
d[s] = 0;
que.push(P(0, s));
while(!que.empty()){
P p = que.top();que.pop();
int v = p.second;
if(d[v] < p.first) continue;
for(int i = 0; i < g[v].size(); i++){
edge e = g[v][i];
if(d[e.to] > d[v] + e.cost){
d[e.to] = d[v] + e.cost;
que.push(P(d[e.to], e.to));
}
}
}
}
int main(){
while(scanf("%d%d",&N,&M)!=EOF){ // n个顶点，m条边
fill(d, d+MAX, INF); //初始化
vector<edge> g[MAX]; //邻接表法
for(int i = 0;i<M;i++){
int a,b,d;
scanf("%d%d%d",&a,&b,&d); //a到b有一条边，边的权值是d
g[a].push_back(edge(b,d)); //点a后接b
g[b].push_back(edge(a,d)); //点b后接a；无向图
}
int s,t;
cin>>s>>t;
Dijkstra(s,g); //s是起点，g是图的邻接表；用d数组保存s到所有点的最小距离
if(d[t] >= INF) d[t] = -1;
cout<<d[t]<<endl;
}
}Matlab实现：（使用：输入邻接矩阵，起点终点，输出mydistance是最短距离，）
function [mydistance,mypath]=mydijkstra(a,sb,db);
% 输入：a—邻接矩阵(aij)是指i到j之间的距离，可以是有向的
% sb—起点的标号, db—终点的标号
% 输出：mydistance—最短路的距离, mypath—最短路的路径
a = zeros(6);
a(1,1) = 2; …//aij是权值
n=size(a,1); visited(1:n) = 0;
distance(1:n) = inf; % 保存起点到各顶点的最短距离
sb = 1;
db = 5; //起点终点
distance(sb) = 0; parent(1:n) = 0;
for i = 1: n-1
temp=distance;
id1=find(visited==1); %查找已经标号的点
temp(id1)=inf; %已标号点的距离换成无穷
[t, u] = min(temp); %找标号值最小的顶点
visited(u) = 1; %标记已经标号的顶点
id2=find(visited==0); %查找未标号的顶点
for v = id2
if a(u, v) + distance(u) < distance(v)
distance(v) = distance(u) + a(u, v); %修改标号值
parent(v) = u;
end
end
end
mypath = [];
if parent(db) ~= 0 %如果存在路!
t = db; mypath = [db];
while t ~= sb
p = parent(t);
mypath = [p mypath];
t = p;
end
end
mydistance = distance(db);
return
Flod算法，单源最短路，权可正可负，但不允许负环。迪杰斯特拉只允许正权
C实现：
int d[10010];
int pre[10010]; //记录路径
struct edge{
int from,to,cost;
}edges[10010];
int V,E,original; //V顶点数，E边数 起点
bool Bellmax_ford(int s){
//初始化
for(int i = 0;i<V;i++){
d[i] = INF;
}
d[s] = 0;
//进行循环，循环下标为从1到n－1（n等于图中点的个数）。在循环内部，遍历所有的边，进行松弛计算。
for(int j = 1;j <= V;j++)
for(int i = 1;i <= E;i++){
edge e = edges[i];
if(d[e.to] > d[e.from] + e.cost && d[e.to] != INF){
d[e.to] = d[e.from] + e.cost;
pre[e.to] = e.from;
}
}
//判断有无负环
/*遍历途中所有的边（edge（u，v）），判断是否存在这样情况：
d（v） > d (u) + w(u,v)
则返回false，表示途中存在从源点可达的权为负的回路。*/
bool flag = 0;
for(int i = 0;i<E;i++){
edge e = edges[i];
if(d[e.to] > d[e.from] + e.cost){
flag = 1;
break;
}
}
return flag;
}
void print_path(int root) //打印最短路的路径（反向）
{
while(root != pre[root]) //前驱
{
printf("%d-->", root);
root = pre[root];
}
if(root == pre[root])
printf("%d\n", root);
}
int main()
{
scanf("%d%d%d", &V, &E, &original);
pre[original] = original;
for(int i = 1; i <= E; ++i)
{
scanf("%d%d%d", &edge[i].from, &edge[i].to, &edge[i].cost);
}
if(Bellman_Ford())
for(int i = 1; i <= V; ++i) //每个点最短路
{
printf("%d\n", d[i]);
printf("Path:");
print_path(i);
}
else
printf("have negative circle\n");
return 0;
}
Matlab实现：
function [dist,mypath]=myfloyd(a,sb,db);
% 输入：a—邻接矩阵(aij)是指 i 到 j 之间的距离，可以是有向的
% sb—起点的标号；db—终点的标号
% 输出：dist—最短路的距离；% mypath—最短路的路径
n=size(a,1); path=zeros(n);
for i=1:n
for j=1:n
if a(i,j)~=inf
path(i,j)=j; %j 是 i 的后续点
end
end
end
for k=1:n
for i=1:n
for j=1:n
if a(i,j)>a(i,k)+a(k,j)
a(i,j)=a(i,k)+a(k,j);
path(i,j)=path(i,k);
end
end
end
end
dist=a(sb,db);
mypath=sb; t=sb;
while t~=db
temp=path(t,db);
mypath=[mypath,temp];
t=temp;
mypath//打印路径
end
return
次短路：
#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<vector>
#include<queue>
#include<algorithm>
#define MAX_V 5010
#define MAx_E 100010
using namespace std;
int V,E;
const int INF = 99999999;
int first_short[MAX_V],second_short[MAX_V];//代表最短距离和次短距离
typedef pair<int,int> P;//first代表到该点距离,second代表该点编号,这里也可用结构体代替
struct Edge{int to,cost;};
vector<Edge>G[MAX_V];
void Dijkstra(int a){
priority_queue<P,vector<P>,greater<P> >q;
fill(first_short+1,first_short+1+V,INF);
fill(first_short+1,second_short+1+V,INF);
first_short[1] = 0;
q.push(P(0,1));
while(!q.empty()){
P p = q.top();q.pop();
int v = p.second;
int d = p.first;
if(second_short[v] < d) continue;//这一句话也可以去掉，因为放入堆中的最长也是次短路的距离
for(int i=0;i<G[v].size();i++){
Edge e = G[v][i];
int d2 = d + e.cost;//d2代表走这条边的情况
if(first_short[e.to] > d2){//如果走这条边后到e.to这个点的距离比最短路还小，就交换他们
swap(first_short[e.to],d2);
q.push(P(first_short[e.to],e.to));
}
if(second_short[e.to] > d2 && d2 > first_short[e.to]){//如果d2比最短小，比次短路长，就更新次短路的值
second_short[e.to] = d2;
q.push(P(second_short[e.to],e.to));//这里次短路也要放入堆中，思考思路中的前几句话。
}
}
}
}
//这里求的是从点1 到最后一点n的次短路
int main(void)
{
while(~scanf("%d %d",&V,&E)){
Edge e;
int a,b,c;
for(int i=1;i<=V;i++)
G[i].clear();
for(int i=1;i<=E;i++){//由于是无向图，故这里要存两次
scanf("%d %d %d",&a,&b,&c);
e.to = b,e.cost = c;
G[a].push_back(e);
e.to = a,e.cost = c;
G[b].push_back(e);
}
Dijkstra(1);
printf("%d\n",second_short[V]);
}
return 0;
}
FLOYD算法
可以求出任意两点间的最短路径。
V顶点数 d[i][j] 表示i到j的最短路径长度
void floyd(){
for(int k = 1; k <= V; k++)
for(int i = 1; i <= V; i++)
for(int j = 1; j <= V; j++)
d[i][j] = min(d[i][j], d[i][k] + d[k][j]);
}
树：
树的应用：欲修筑连接 n 个城市的铁路，已知 i 城与 j 城之间的铁路造价为
Cij，设计一个线路图，使总造价最低。
这就是构造最小生成树：（赋权图具有最小的权值）
构造最小生成树：
Matlab实现：
%result的第一、二、三行分别表示生成树边的起点、终点、权集合！
clc;clear;
a=zeros(7);
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;
a=a+a';a(find(a==0))=inf;
result=[];p=1;tb=2:length(a);
while length(result)~=length(a)-1
temp=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))=[];
end
result
C实现:
#include <iostream>
#include <cstdio>
#include <string.h>
#define INF 100000
#define MAXN 1010
using namespace std;
int m, n;
int edge[MAXN][MAXN];
int lowcost[MAXN]; //存权值，存未找到的集合中的各个点到已找到集合(看做一个点)的权值，要更新
//找最小，对应的点，用过是-1
int nearvex[MAXN]; //nearvex[j]记录未找到集合中的点j与 已找到集合的所有点 最邻近的点(存的值) ，-1表示已用
void Prim(int u0){
int i, j;
int sumweight = 0;
for(i = 1; i <= n; i++){
lowcost[i] = edge[u0][i];
nearvex[i] = u0;
}
nearvex[u0] = -1;
for (i = 1; i < n; ++i){ //n个结点 n-1次循环解决问题
int minn = INF;
int v = -1;
//在lowcost数组（的nearvex[]值不为-1的元素中）找最小值
for (j = 1; j <= n; ++j){
if(nearvex[j] != -1 && lowcost[j] < minn){
v = j; //用v保存最小权值的编号
minn = lowcost[j];
}
}
//找到了，更新
if(v != -1){
printf("%d %d %d\n", nearvex[v], v, lowcost[v]); //打印路径的
nearvex[v] = -1;
sumweight += lowcost[v];
for(j = 1; j <= n; j++){
if(nearvex[j]!=-1 && edge[v][j] < lowcost[j]){
lowcost[j] = edge[v][j];
nearvex[j] = v;
}
}
}
}
printf("weight of MST is %d\n", sumweight);
}
int main(int argc, char const *argv[])
{
int i, j;
int u, v, w;
scanf("%d%d", &n, &m);
memset(edge, 0, sizeof(edge));
//用 邻接矩阵 来存放 edge[i][j] 表示 i 点到 j 点的权值 自己到自己0 不可直接达为INF
for (int i = 1; i <= m; ++i)
{
scanf("%d%d%d", &u, &v, &w);
edge[u][v] = edge[v][u] = w;
}
for (int i = 1; i <= n; ++i)
{
for (int j = 1; j <= n; ++j)
{
if(i == j) edge[i][j] = 0;
else if(edge[i][j] == 0)edge[i][j] = INF;
}
}
Prim(1);
return 0;
}流网络：一张图，对于每一个顶点，都有顶点权值（权函数）；每一条边，都有最小、最大两个权函数。
流：流网络中，一个流是弧集A到R的一个函数，即対每一条弧赋一个实数fij，称作这条弧的流量。
如果：fij之和-fji之和 = di，则称f为一个可行流；至少存在一个可行流的流网络成为可行网络。即对于弧集，对于任意一个顶点i，它到任意顶点j的流之和减去j到它的流之和，等于顶点的权值（权函数值）；则就是可行流。理解为对于该点流量守恒。因此，对于可行网络，所有顶点的流量和（权）等于0.
最大流问题：单源单汇点，可行流网络中，找到流值最大的可行流，即最大流。（特殊的线性规划）
最大流的Ford-Fullkerson算法（计算最大流）：
clc,clear;%用u矩阵存放最大限制，（最小默认0）；用f存当前流
u(1,2)=1;u(1,3)=1;u(1,4)=2;u(2,3)=1;u(2,5)=2;
u(3,5)=1;u(4,3)=3;u(4,5)=3;
f(1,2)=1;f(1,3)=0;f(1,4)=1;f(2,3)=0;f(2,5)=1;
f(3,5)=1;f(4,3)=1;f(4,5)=0;
n=length(u);list=[];maxf(n)=1;
while maxf(n)>0
maxf=zeros(1,n);pred=zeros(1,n);
list=1;record=list;maxf(1)=inf;
%list是未检查邻接点的标号点，record是已标号点
while (~isempty(list))&(maxf(n)==0)
flag=list(1);list(1)=[];
label1= find(u(flag,:)-f(flag,:));
label1=setdiff(label1,record);
list=union(list,label1);
pred(label1)=flag;
maxf(label1)=min(maxf(flag),u(flag,label1)...
-f(flag,label1));
record=union(record,label1);
label2=find(f(:,flag));
label2=label2';
label2=setdiff(label2,record);
list=union(list,label2);
pred(label2)=-flag;
maxf(label2)=min(maxf(flag),f(label2,flag));
record=union(record,label2);
end
if maxf(n)>0
v2=n; v1=pred(v2);
while v2~=1
-95-
if v1>0
f(v1,v2)=f(v1,v2)+maxf(n);
else
v1=abs(v1);
f(v2,v1)=f(v2,v1)-maxf(n);
end
v2=v1; v1=pred(v2);
end
end
end
f %打印最后结果
c实现：
1. /**
2. * Ford-Fulkerson方法的一种实现
3. * @param c 二维矩阵，记录每条边的容量
4. * @param vertexNum 顶点个数，包括起点和终点
5. * @param s 起点编号，编号从1开始
6. * @param t 终点编号
7. * @param f 输出流网络矩阵，二维矩阵，记录每条边的流量
8. */
9. void Ford_Fulkerson(int **c, int vertexNum, int s, int t, int
10. {
11. int *e = (int *)malloc(sizeof(int)*vertexNum*vertexNum); // 残存网络
12. int *priorMatrix = (int *)malloc(sizeof(int)*vertexNum); // 增广路径的前驱子图
13.
14. // initialize
15. for (int
16. {
17. for (int
18. {
19. *(f + i*vertexNum + j) = 0;
20. }
21. }
22.
23. while
24. {
25. calculateENet(c, vertexNum, (int **)f, (int
26. int
27. if ((min = findRoute((int
28. {
29. break;
30. }
31. int
32. int
33. while
34. {
35. if (*((int*)c + pre * vertexNum + next) != 0)
36. {
37. *((int*)f + pre * vertexNum + next) += min;
38. }
39. else
40. {
41. *((int*)f + next * vertexNum + pre) -= min;
42. }
43. next = pre;
44. pre = priorMatrix[pre];
45. }
46. }
47. }
测试：
1. void
2. {
3. int
4. 0, 0, 0, 12, 0, 0,
5. 0, 4, 0, 0, 14, 0,
6. 0, 0, 9, 0, 0, 20,
7. 0, 0, 0, 7, 0, 4,
8. 0, 0, 0, 0, 0, 0 };
9. int
10. Ford_Fulkerson((int **)c, 6, 1, 6, (int
11. for (int
12. {
13. for (int
14. {
15. int
16. if
17. {
18. printf("%d -> %d : %d\n", i + 1, j + 1, flow);
19. }
20. }
21. }
22. }
最小费用流：
在运输问题中，人们总是希望在完成运输任务的同时，寻求一个使总的运输费用最小的运输方案。这个就是最小费用流。
最小费用最大流问题：
function mainexample19
clear;clc;
global M num
c=zeros(6);u=zeros(6);
%c是费用，u是最大容量
c(1,2)=2;c(1,4)=8;c(2,3)=2;c(2,4)=5;
c(3,4)=1;c(3,6)=6;c(4,5)=3;c(5,3)=4;c(5,6)=7;
u(1,2)=8;u(1,4)=7;u(2,3)=9;u(2,4)=5;
u(3,4)=2;u(3,6)=5;u(4,5)=9;u(5,3)=6;u(5,6)=10;
num=size(u,1);M=sum(sum(u))*num^2;
[f,val]=mincostmaxflow(u,c) %调用，返回f是可行流矩阵，val是最小费用的值；
%u是容量矩阵，c是费用矩阵。可以再加一个指定流量的参数。
%求最短路径函数
function path=floydpath(w);
global M num
w=w+((w==0)-eye(num))*M;
p=zeros(num);
for k=1:num
for i=1:num
for j=1:num
if w(i,j)>w(i,k)+w(k,j)
w(i,j)=w(i,k)+w(k,j);
p(i,j)=k;
end
end
end
end
if w(1,num) ==M
path=[];
else
path=zeros(num);
s=1;t=num;m=p(s,t);
while ~isempty(m)
if m(1)
s=[s,m(1)];t=[t,t(1)];t(1)=m(1);
m(1)=[];m=[p(s(1),t(1)),m,p(s(end),t(end))];
else
path(s(1),t(1))=1;s(1)=[];m(1)=[];t(1)=[];
end
end
end%最小费用最大流函数函数
function [flow,val]=mincostmaxflow(rongliang,cost,flowvalue);
%第一个参数：容量矩阵；第二个参数：费用矩阵；
%前两个参数必须在不通路处置零
%第三个参数：指定容量值（可以不写，表示求最小费用最大流）
%返回值 flow 为可行流矩阵,val 为最小费用值
global M
flow=zeros(size(rongliang));allflow=sum(flow(1,:));
if nargin<3
flowvalue=M;
end
while allflow<flowvalue
w=(flow<rongliang).*cost-((flow>0).*cost)';
path=floydpath(w);%调用 floydpath 函数
if isempty(path)
val=sum(sum(flow.*cost));
return;
end
theta=min(min(path.*(rongliang-flow)+(path.*(rongliang-flow)==0).*M));
theta=min([min(path'.*flow+(path'.*flow==0).*M),theta]);
flow=flow+(rongliang>0).*(path-path').*theta;
allflow=sum(flow(1,:));
end
val=sum(sum(flow.*cost));以上，请批评指正。
