【数学建模】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. } 

clip_image001[6]

最小费用流:

在运输问题中,人们总是希望在完成运输任务的同时,寻求一个使总的运输费用最小的运输方案。这个就是最小费用流。

最小费用最大流问题:

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));

posted @ 2018-08-07 11:40  pigcv  阅读(1275)  评论(1编辑  收藏  举报