一、向量、矩阵的表示和使用
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 **f)
10. {
11. int *e = (int *)malloc(sizeof(int)*vertexNum*vertexNum); // 残存网络
12. int *priorMatrix = (int *)malloc(sizeof(int)*vertexNum); // 增广路径的前驱子图
13.
14. // initialize
15. for (int i = 0; i < vertexNum;i++)
16. {
17. for (int j = 0; j < vertexNum; j++)
18. {
19. *(f + i*vertexNum + j) = 0;
20. }
21. }
22.
23. while (1)
24. {
25. calculateENet(c, vertexNum, (int **)f, (int **)e); // 计算残存网络
26. int min;
27. if ((min = findRoute((int **)e, vertexNum, priorMatrix, s, t)) == 0) // 寻找增广路径及其最小流值
28. {
29. break;
30. }
31. int pre = priorMatrix[t - 1];
32. int next = t - 1;
33. while (pre != -1) // 按增广路径更新流网络
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 testFord()
2. {
3. int c[6][6] = { 0, 16, 13, 0, 0, 0,
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 f[6][6];
10. Ford_Fulkerson((int **)c, 6, 1, 6, (int **)f);
11. for (int i = 0; i < 6; i++)
12. {
13. for (int j = 0; j < 6; j++)
14. {
15. int flow = f[i][j];
16. if (flow != 0)
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));