【数学建模】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 = 29

A = [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) * b

rank 求方阵的秩
求逆矩阵：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));