数学建模——详解弗洛伊德(Folyd)算法【分别用 C/C++ 和 matlab 实现】
关键字:Folyd算法的背景小故事、算法实现流程(配合样例+图例)、用C/C++和 Matlab分别实现。
文章目录
今天势必拿下 Folyd !🍋 🍋
数学建模系列文章——总结篇:《数模美一国一退役选手的经验分享[2021纪念版]》.
一、Folyd算法的背景小故事
学习一门东西,如果想要记得久,那就得了解一些关于它的历史。
感性的人文 + 理性的推导 = 长久的记忆 感性的人文 + 理性的推导 = 长久的记忆 感性的人文+理性的推导=长久的记忆
我就 简单 地讲讲创造该算法的 弗洛伊德 的传奇小故事吧。
(1936-2001)Robert W.Floyd
历届图灵奖得主基本上都有高学历、高学位,而且绝大多数都有博士头衔。这是可以理解的,因为创新型人才需要有很好的文化素养,丰富的知识底蕴,因而必须接受良好的教育。
但上帝总是有遗漏的骰子。1978年图灵奖获得者就是 斯坦福大学计算机科学系教授 罗伯特·弗洛伊德 ,而他却是一位 “自学成才的计算机科学家”。
弗洛伊德1936年6月8日生于纽约。说他 “自学成才” 并不是说他没有接受过高等教育,他是芝加哥大学的毕业生,但学的不是数学或电气工程等与计算机密切相关的专业,而是 文学 ,1953年 获得 文学士学位 。【17岁那时他才!】
20世纪50年代初期美国经济不太景气,找工作比较困难,因学习文学而 没有任何专门技能 的弗洛伊德在就业上遇到很大麻烦。
无奈之中到西屋电气公司当了二年 计算机操作员 ,在IBM650机房值夜班。我们知道,早期的计算机都是以批处理方式工作的,计算机操作员的任务就是把程序员编写好的程序在卡片穿孔机(这是脱机的辅助外部设备)上穿成卡片,然后把卡片叠放在读卡机上输入计算机,以便运行程序。【早期的搞计算机看来 还得有点体力值才行】
因此,操作员的工作比较简单,同打字员类似,不需要懂计算机,也不需要懂程序设计。但弗洛伊德毕竟是一个受过高等教育的人,又是一个有心人,干了一段时间的操作员,很快对计算机产生了兴趣,决心弄懂它,掌握它。
于是他借了有关书籍资料在值班空闲时间刻苦学习钻研,有问题就虚心向程序员请教。白天不值班,他又回母校去听讲有关课程。这样,他不但在 1958年 又获得了 理科学士学位 ,而且逐渐从计算机的门外汉变成计算机的行家里手。
【时隔五年,直接来一波“孔夫子学Java——能文能码”】
1956年他离开西屋电气公司,到芝加哥的装甲研究基金会。开始还是当操作员,后来就当了程序员。1962年他被马萨诸塞州的Computer Associates公司聘为分析员。此时与 Warsall 合作发布 Floyd-Warshall 算法。1965年他应聘成为卡内基—梅隆大学的副教授,3年后转至斯坦福大学。1970年被聘任为教授。
之所以能这样快地步步高升,关键就在于弗洛伊德通过 勤奋学习 和 深入研究 ,在计算机科学的诸多领域:算法,程序设计语言的逻辑和语义,自动程序综合,自动程序验证,编译器的理论和实现等方面都作出创造性的贡献。
前辈尚如此,吾辈当自强。淦淦淦!好好啃下前辈传递下来的人类智慧吧!
二、Folyd算法简介
Folyd算法又称为插点法,是一种利用 动态规划 的思想,寻找给定的加权图中 任意两点之间的最短路径 的算法。
它能可以正确处理 有向图 或 无向图 或 负权(但不可存在负权回路) 的最短路径问题,同时也被用于计算有向图的传递闭包。
三、算法实现流程(配合样例)
假如,我们有下面这张图。现在要求找出从某一点到达任意一点的最短距离。怎么做?
当然,这么简单的图,大家 肯定一👀就能看出来。(专门把数值大的线加粗了)
第一步:初始化距离矩阵
首先,我们要构造一个距离矩阵 D D D ( distance 的首字母),如下图所示:
说明:
①第
i
i
i 行第
j
j
j 列的数字代表结点
i
i
i 到结点
j
j
j 的直达距离。即
D
[
i
]
[
j
]
D[i][j]
D[i][j] 。
② D [ i ] [ j ] = ∞ D[i][j]=∞ D[i][j]=∞ 表示的是结点 i i i 没有直达结点 j j j 的路径。
③出现 D [ i ] [ j ] = 0 D[i][j]=0 D[i][j]=0 的情况都是在对角线上,即 i = j i=j i=j 。这里的意思是 结点自身到自身的距离我们认为是 0。(不考虑环)
第二步:“解封”死结点 🔒→再动态更新距离矩阵
怎么理解 “解封死结点” 呢?emmm…为了方便大家理解。大家先看看下面这张图。
我们总共有4个结点,蓝色的结点我通通都给它们换了一个新名字。我把它们称之为 “死结点” 。 而 1 号结点就是我最先 解封 的结点,一解封后,它就成了 “活结点”。
“活”,代表“灵活、活动”。“死”,代表“死板,僵硬”。之所有把 1号结点 称之为 “活结点”,就是因为,当解封它后,原先 2号结点 不能直达 4号结点的,但现在可以通过 1号“活”结点 中转 一下就到了。
所以 D [ 2 ] [ 3 ] D[2][3] D[2][3] 从原先的 ∞ ∞ ∞ 变为了红框里面的 13 13 13 。同理 D [ 2 ] [ 4 ] D[2][4] D[2][4] 也从原先的 ∞ ∞ ∞ 变为了红框里面的 5 5 5
而 “死结点” 就不能像 “活结点” 那样实现中转这个功能。比如说,此时 3号结点 就不能通过现有路径达到 1号结点。很显然是因为 2号结点 现在还是 “死结点”,它现在还不能实现中转功能。
🍔 🍔 🍔
类似地,我们再依次解封,就能慢慢实现我们的目标了。
但是在解封的过程中,我们要遵循一个原则。
那就是在每次解封一个结点后,我们要判断一下:除该解封结点外,其他任意两个结点之间的距离是否有缩短。⭐️ ⭐️
比如说,我们再次解封 2号结点,如下图所示:
同理,在右边的距离矩阵中, D [ 3 ] [ 1 ] D[3][1] D[3][1] 从原先的 ∞ ∞ ∞ 变为了红框里面的 4 4 4 。
又因为,在没有解封2号结点前,3号结点能 直达 4号结点,故 D [ 3 ] [ 4 ] = 7 D[3][4]=7 D[3][4]=7 。但现在解封了2号结点后,3号结点能在 2号结点 中转 一下到 1号结点,再 中转 一下就可以抵达4号结点。
虽然感觉绕了好大一圈。但是我们一计算路径,发现 “ 1 + 3 + 2 = 6 < 7 1+3+2=6<7 1+3+2=6<7 ”,3号结点→4号结点 的路程变得更短了!!!所以红框里面的 D [ 3 ] [ 4 ] D[3][4] D[3][4] 从原先的 7 7 7 变为了大红框里面的 6 6 6。
🔑 🔑 🔑
到这里,我们即可得到下面这个状态转移方程。 D [ i ] [ j ] = m i n ( D [ i ] [ j ] , D [ i ] [ k ] + D [ k ] [ j ] ) D[i][j]=min(D[i][j],D[i][k]+D[k][j]) D[i][j]=min(D[i][j],D[i][k]+D[k][j]) 这里面的 k k k 即代表 新解封 的结点标号。
举个栗子,如上面所说, D [ 3 ] [ 1 ] D[3][1] D[3][1] 从原先的 ∞ ∞ ∞ 变为了红框里面的 4 4 4, 就是通过 “ D [ 3 ] [ 1 ] = m i n ( D [ 3 ] [ 1 ] , D [ 3 ] [ 2 ] + D [ 2 ] [ 1 ] ) D[3][1]=min(D[3][1],D[3][2]+D[2][1]) D[3][1]=min(D[3][1],D[3][2]+D[2][1]) ”,即 “ D [ 3 ] [ 1 ] = m i n ( ∞ , 1 + 3 ) D[3][1]=min(∞,1+3) D[3][1]=min(∞,1+3) ” 得到。
之所以叫 状态转移方程 ,是因为 每次新解封一个死结点后,整个矩阵 D D D 上面的所有项都要更新一遍 “状态”。看在新解封结点后, D [ i ] [ j ] D[i][j] D[i][j] 小,还是 D [ i ] [ k ] + D [ k ] [ j ] D[i][k]+D[k][j] D[i][k]+D[k][j] 更小,哪个更小,就转移到那个状态。【注:如果是 D [ i ] [ j ] D[i][j] D[i][j] 转移到 D [ i ] [ j ] D[i][j] D[i][j],即代表状态未变】
我们接着解封剩下的两个结点,过程如下:【红框部分是变化了的元素值】
第三步:输出最短距离
当我们解封完所有结点后,我们即可得到最终的答案矩阵。该矩阵上的记录的便是两两结点之间的最短距离。用两层 f o r for for 循环输出即可。
第四步:输出最短路径(补充内容)
通过第三步,我们只知道 3 号结点到 4 号结点的最短距离为 6 。 但并不知道 3 号结点是怎么走到 4 号结点的。为此,我们将在
第二步 “ ‘解封’ 死结点 🔒→再动态更新距离矩阵 ” 之中,再加一个新公式,用它来记录最短路径。
在新解封一个结点 k 后 ; 在新解封一个 结点k 后; 在新解封一个结点k后;
i
f
(
D
[
i
]
[
k
]
+
D
[
k
]
[
j
]
<
D
[
i
]
[
j
]
)
if(D[i][k]+D[k][j]<D[i][j])
if(D[i][k]+D[k][j]<D[i][j])
{
\{
{
D
[
i
]
[
j
]
=
D
[
i
]
[
k
]
+
D
[
k
]
[
j
]
; //更新矩阵
D[i][j]=D[i][k]+D[k][j]\text{; //更新矩阵}
D[i][j]=D[i][k]+D[k][j]; //更新矩阵
p
a
t
h
[
i
]
[
j
]
=
p
a
t
h
[
k
]
[
j
]
; //更新路径
path[i][j]=path[k][j]\text{; //更新路径}
path[i][j]=path[k][j]; //更新路径
}
\}
}
在上式中, p a t h [ i ] [ j ] path[i][j] path[i][j] 表示:在 i i i 号结点到 j j j 号结点最短路径上的 j j j 号结点的 前驱结点 。
当然 p a t h path path 矩阵和 矩阵 D D D 一样,都是需要初始化、动态更新的。见下述步骤:
s t e p 1 step1 step1:初始化路径矩阵
路径矩阵 p a t h path path 和 矩阵 D D D 是在同一个 f o r for for 循环里进行初始化。得到结果如下图所示:
说明:
①
p
a
t
h
[
i
]
[
j
]
=
0
path[i][j]=0
path[i][j]=0 代表结点
i
i
i 到结点
j
j
j 之间目前还没有最短路径(因为还没有解封死结点)。
②因为结点自身到自身的最短距离为0,所以在该路径上,该结点的前驱结点就是它本身。因此 p a t h path path 矩阵的对角线为 1 、 2 、 3 、 4 1、2、3、4 1、2、3、4 。
③ p a t h [ i ] [ j ] path[i][j] path[i][j] 存储的是 i→j 的最短路径上的 j 结点 的前一个结点编号。
s t e p 2 step2 step2:动态更新路径矩阵
然后,我们在更新距离矩阵时,也同时更新路径矩阵。演示图如下:
从图中我们可以看到,当
D
[
i
]
[
j
]
D[i][j]
D[i][j] 更新时,
p
a
t
h
[
i
]
[
j
]
path[i][j]
path[i][j] 也要在 “相应的” 位置更新。更新的规则就是先前所述公式:
在新解封一个结点 k 后 ; 在新解封一个 结点k 后; 在新解封一个结点k后;
i
f
(
D
[
i
]
[
k
]
+
D
[
k
]
[
j
]
<
D
[
i
]
[
j
]
)
if(D[i][k]+D[k][j]<D[i][j])
if(D[i][k]+D[k][j]<D[i][j])
{
\{
{
D
[
i
]
[
j
]
=
D
[
i
]
[
k
]
+
D
[
k
]
[
j
]
; //更新矩阵
D[i][j]=D[i][k]+D[k][j]\text{; //更新矩阵}
D[i][j]=D[i][k]+D[k][j]; //更新矩阵
p
a
t
h
[
i
]
[
j
]
=
p
a
t
h
[
k
]
[
j
]
; //更新路径
path[i][j]=path[k][j]\text{; //更新路径}
path[i][j]=path[k][j]; //更新路径
}
\}
}
到这里,有人可能会问,“为什么等号右边的是 path[k][j] , 不是path[i][k]、也不是 k 等等等? ”
我来举个栗子:🙌 🥜 (其实是花生 ) 便于大家理解。
如上图所示,红色结点都是已经被解封的 “活结点”,且现在有一个即将被解封的 “k号死结点”。
通过观察,易知,k号结点没解封前,活结点 i i i 到活结点 j j j 的最短路径便是直达路径(即是 “ i → j i → j i→j ”)。直达距离为 15 ,此时 p a t h [ i ] [ j ] = i path[i][j]=i path[i][j]=i。
当结点 k k k 被解封后,易知,结点 i i i 到结点 j j j 的最短路径现在就更新为 “ i → k → 2 → 3 → j i → k →2→3→j i→k→2→3→j ”。那么 p a t h [ i ] [ j ] path[i][j] path[i][j] 现在就更新为 3 3 3 。
那这个 3 3 3 哪里来的,是哪里的 3 3 3 ?
不就是在解封结点 2 和结点 3 后, p a t h [ 2 ] [ j ] path[2][j] path[2][j]、 p a t h [ 3 ] [ j ] path[3][j] path[3][j] 更新得到并存储下来的 3 3 3 ,是从这里来的。 ⭐️ ⭐️ ⭐️ 关键点
四、代码实现:
① C/C++语言:
#include<cstdio> // 如果是 C 环境,请注释掉该语句
#include<stdio.h> // 如果是 C++ 环境,请注释掉该语句
const int N =100;
int main()
{
int i, j, k, D[N][N], path[N][N];
int u, v, JieDian_num; // u代表起始结点,v代表终止结点, JieDian_num表示结点个数
int e_num, w; // e_num代表边的条数,w代表边的权值
printf("请输入分别输入结点个数和边的条数:");
scanf("%d%d",&JieDian_num,&e_num);
for( i=1;i<=JieDian_num;i++ )
{
for( j=1;j<=JieDian_num;j++ )
{
D[i][j] = 999999; // 距离矩阵初始化(这里用999999假装代表无穷大)
D[i][i] = 0; // 但自己到自己的距离还是 0
path[i][i] = i; // 路径矩阵初始化
}
}
printf("\n***请输入每条边的起始结点编号、终止结点编号和连接它俩的边的权值***\n");
for( i=1;i<=e_num;i++ )
{
printf("第%d条边:",i);
scanf("%d%d%d",&u,&v,&w);
D[u][v] = w; // 距离矩阵初始化
path[u][v] = u; // 路径矩阵初始化
}
for( k=1;k<=JieDian_num;k++ ) // 每次新“解封”一个结点
{
for( i=1;i<=JieDian_num;i++ )
{
for( j=1;j<=JieDian_num;j++ )
{
if( D[i][k] + D[k][j] < D[i][j] )
{
D[i][j] = D[i][k] + D[k][j]; // 动态更新距离矩阵
path[i][j] = path[k][j]; // 动态更新路径矩阵
}
}
}
}
printf("\n距离矩阵的结果如下:\n");
for( i=1;i<=JieDian_num;i++ ) // 输出
{
for( j=1;j<=JieDian_num;j++ )
{
printf("%d ",D[i][j]);
}
printf("\n");
}
printf("\n路径矩阵的结果如下:\n");
for( i=1;i<=JieDian_num;i++ ) // 输出
{
for( j=1;j<=JieDian_num;j++ )
{
printf("%d ",path[i][j]);
}
printf("\n");
}
return 0;
}
/* 样例:
4 7
1 2 8
1 3 10
1 4 2
2 1 3
3 2 1
3 4 7
4 3 4
*/
运行结果:
可以看到,距离矩阵 D D D 的答案和原先手算的一直。 p a t h path path 矩阵通过手推,也可以得到上述答案。
② matlab:
clc,clear,close all;
v = 4; % 结点个数
e = 7; % 边的条数
D = [ 0 8 10 2 ;
3 0 inf inf;
inf 1 0 7;
inf inf 4 0] % 初始化距离矩阵
path = [1 1 1 1;
2 2 0 0;
0 3 3 3;
0 0 4 4] % 初始化路径矩阵
% Folyd算法:
for k=1:v
for i=1:v
for j=1:v
if D(i,j) > D(i,k) + D(k,j)
D(i,j) = D(i,k) + D(k,j); % 更新最短距离
path(i,j) = path(k,j); % 更新最短路径
end
end
end
end
disp('最终的距离矩阵如下:')
D
disp('最终的路径矩阵如下:')
path
运行结果:
可以看到,答案与 C/C++ 跑出来的一致!😲
五、Folyd算法的特点
①Floyd算法适用于APSP(All Pairs Shortest Paths),是一种动态规划算法,稠密图效果最佳,边权可正可负。
②此算法简单有效,由于三重循环结构紧凑,对于稠密图,效率要高于执行 ∣ V ∣ |V| ∣V∣ 次的Dijkstra算法。
③优点:容易理解,可以算出任意两个节点之间的最短距离,代码编写简单。
④缺点:时间复杂度比较高,不适合计算大量数据。
⑤时间复杂度: O(n3);空间复杂度: O(n2);
六、总结:
首先从 Folyd算法的背景小故事 讲起,弗洛伊德 的传奇经历属实非同凡响。
然后讲解了 Folyd算法的简介 和 具体的实现过程,并搭配了详细的图例,最后用两种语言进行测试。最后简单列了一下算法的特点。
整个Folyd算法里面,最重要的就是要理解那 两个状态转移方程,为此我再列一下:
在新解封一个结点 k 后 ; 在新解封一个 结点k 后; 在新解封一个结点k后;
i
f
(
D
[
i
]
[
k
]
+
D
[
k
]
[
j
]
<
D
[
i
]
[
j
]
)
if(D[i][k]+D[k][j]<D[i][j])
if(D[i][k]+D[k][j]<D[i][j])
{
\{
{
D
[
i
]
[
j
]
=
D
[
i
]
[
k
]
+
D
[
k
]
[
j
]
; //更新矩阵
D[i][j]=D[i][k]+D[k][j]\text{; //更新矩阵}
D[i][j]=D[i][k]+D[k][j]; //更新矩阵
p
a
t
h
[
i
]
[
j
]
=
p
a
t
h
[
k
]
[
j
]
; //更新路径
path[i][j]=path[k][j]\text{; //更新路径}
path[i][j]=path[k][j]; //更新路径
}
\}
}
七、部分参考附录:
[1] 《总结一下最短路径的弗洛伊德算法(Floyd)》
链接: https://blog.csdn.net/riba2534/article/details/54562440.
[2] 《最短路径模板+解析——(FLoyd算法)》
链接: https://blog.csdn.net/ytuyzh/article/details/88617987.
数学建模系列文章——总结篇:《数模美一国一退役选手的经验分享[2021纪念版]》.
🤒 终于写完这一篇了!!!从白天写到黑夜,12000字,哎…夜深人静了
码字码图不易,觉得好的话点个 👍 ,满足一下我的虚荣心~ ❤️。我还会继续创作,给我一点继续写文的动力吧~ 感谢。🍀
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)