数学建模——详解弗洛伊德(Folyd)算法【分别用 C/C++ 和 matlab 实现】

大家好，欢迎来到IT知识分享网。

今天势必拿下 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$ ( distance 的首字母)，如下图所示：

  说明：
    ①第 $i$ 行第 $j$ 列的数字代表结点 $i$ 到结点 $j$ 的直达距离。即 $D[i][j]$ 。

    ② $D[i][j]=\infty$ 表示的是结点 $i$ 没有直达结点 $j$ 的路径。

    ③出现 $D[i][j]=0$ 的情况都是在对角线上，即 $i=j$ 。这里的意思是 结点自身到自身的距离我们认为是 0。(不考虑环)

第二步：“解封”死结点 🔒→再动态更新距离矩阵

  怎么理解 “解封死结点” 呢？emmm…为了方便大家理解。大家先看看下面这张图。

  我们总共有4个结点，蓝色的结点我通通都给它们换了一个新名字。我把它们称之为 “死结点” 。 而 1 号结点就是我最先 解封 的结点，一解封后，它就成了 “活结点”。

  “活”，代表“灵活、活动”。“死”，代表“死板，僵硬”。之所有把 1号结点 称之为 “活结点”，就是因为，当解封它后，原先 2号结点 不能直达 4号结点的，但现在可以通过 1号“活”结点 中转 一下就到了。

  所以 $D[2][3]$ 从原先的 $\infty$ 变为了红框里面的 $13$ 。同理 $D[2][4]$ 也从原先的 $\infty$ 变为了红框里面的 $5$

  而 “死结点” 就不能像 “活结点” 那样实现中转这个功能。比如说，此时 3号结点 就不能通过现有路径达到 1号结点。很显然是因为 2号结点 现在还是 “死结点”，它现在还不能实现中转功能。

🍔 🍔 🍔

  类似地，我们再依次解封，就能慢慢实现我们的目标了。

  但是在解封的过程中，我们要遵循一个原则。

  那就是在每次解封一个结点后，我们要判断一下：除该解封结点外，其他任意两个结点之间的距离是否有缩短。⭐️ ⭐️

  比如说，我们再次解封 2号结点，如下图所示：

  同理，在右边的距离矩阵中，$D[3][1]$ 从原先的 $\infty$ 变为了红框里面的 $4$ 。

  又因为，在没有解封2号结点前，3号结点能 直达 4号结点，故 $D[3][4]=7$ 。但现在解封了2号结点后，3号结点能在 2号结点 中转 一下到 1号结点，再 中转 一下就可以抵达4号结点。

  虽然感觉绕了好大一圈。但是我们一计算路径，发现 “ $1+3+2=6<7$ ”，3号结点→4号结点 的路程变得更短了！！！所以红框里面的 $D[3][4]$ 从原先的 $7$ 变为了大红框里面的 $6$。

🔑 🔑 🔑

  到这里，我们即可得到下面这个状态转移方程。$$D[i][j]=min(D[i][j]，D[i][k]+D[k][j])$$  这里面的 $k$ 即代表 新解封 的结点标号。

  举个栗子，如上面所说， $D[3][1]$ 从原先的 $\infty$ 变为了红框里面的 $4$， 就是通过 “ $D[3][1]=min(D[3][1]，D[3][2]+D[2][1])$ ”，即 “ $D[3][1]=min(\infty ，1+3)$ ” 得到。

  之所以叫 状态转移方程 ，是因为 每次新解封一个死结点后，整个矩阵 $D$ 上面的所有项都要更新一遍 “状态”。看在新解封结点后， $D[i][j]$ 小，还是 $D[i][k]+D[k][j]$ 更小，哪个更小，就转移到那个状态。【注：如果是 $D[i][j]$ 转移到 $D[i][j]$，即代表状态未变】

  我们接着解封剩下的两个结点，过程如下：【红框部分是变化了的元素值】

第三步：输出最短距离

  当我们解封完所有结点后，我们即可得到最终的答案矩阵。该矩阵上的记录的便是两两结点之间的最短距离。用两层 $for$ 循环输出即可。

第四步：输出最短路径（补充内容）

  通过第三步，我们只知道 3 号结点到 4 号结点的最短距离为 6 。 但并不知道 3 号结点是怎么走到 4 号结点的。为此，我们将在
第二步 “ ‘解封’ 死结点 🔒→再动态更新距离矩阵 ” 之中，再加一个新公式，用它来记录最短路径。

  $在新解封一个结点k后;$

  $if(D[i][k]+D[k][j]<D[i][j])$
   { \{ {

    $D[i][j]=D[i][k]+D[k][j]\text{; //更新矩阵}$
    $path[i][j]=path[k][j]\text{; //更新路径}$
  $\}$

  在上式中， $path[i][j]$ 表示：在 $i$ 号结点到 $j$ 号结点最短路径上的 $j$ 号结点的 前驱结点 。

  当然 $path$ 矩阵和 矩阵 $D$ 一样，都是需要初始化、动态更新的。见下述步骤：

$step1$：初始化路径矩阵

  路径矩阵 $path$ 和 矩阵 $D$ 是在同一个 $for$ 循环里进行初始化。得到结果如下图所示：

  说明：
    ①$path[i][j]=0$ 代表结点 $i$ 到结点 $j$ 之间目前还没有最短路径(因为还没有解封死结点)。

    ②因为结点自身到自身的最短距离为0，所以在该路径上，该结点的前驱结点就是它本身。因此 $path$ 矩阵的对角线为 $1、2、3、4$ 。

    ③$path[i][j]$ 存储的是 i→j 的最短路径上的 j 结点 的前一个结点编号。

$step2$：动态更新路径矩阵

  然后，我们在更新距离矩阵时，也同时更新路径矩阵。演示图如下：

  从图中我们可以看到，当 $D[i][j]$ 更新时，$path[i][j]$ 也要在 “相应的” 位置更新。更新的规则就是先前所述公式：

  $在新解封一个结点k后;$

  $if(D[i][k]+D[k][j]<D[i][j])$
   { \{ {

    $D[i][j]=D[i][k]+D[k][j]\text{; //更新矩阵}$
    $path[i][j]=path[k][j]\text{; //更新路径}$
  $\}$

  到这里，有人可能会问，“为什么等号右边的是 path[k][j] ， 不是path[i][k]、也不是 k 等等等? ”

  我来举个栗子：🙌 🥜 (其实是花生 ) 便于大家理解。

  如上图所示，红色结点都是已经被解封的 “活结点”，且现在有一个即将被解封的 “k号死结点”。

  通过观察，易知，k号结点没解封前，活结点 $i$ 到活结点 $j$ 的最短路径便是直达路径(即是 “ $i\to j$ ”)。直达距离为 15 ，此时 $path[i][j]=i$。

  当结点 $k$ 被解封后，易知，结点 $i$ 到结点 $j$ 的最短路径现在就更新为 “ $i\to k\to 2\to 3\to j$ ”。那么 $path[i][j]$ 现在就更新为 $3$ 。

  那这个 $3$ 哪里来的，是哪里的 $3$ ？

  不就是在解封结点 2 和结点 3 后，$path[2][j]$、$path[3][j]$ 更新得到并存储下来的 $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] = ; // 距离矩阵初始化（这里用假装代表无穷大） 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$ 的答案和原先手算的一直。$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|$ 次的Dijkstra算法。

  ③优点：容易理解，可以算出任意两个节点之间的最短距离，代码编写简单。

  ④缺点：时间复杂度比较高，不适合计算大量数据。

  ⑤时间复杂度: O(n³)；空间复杂度: O(n²)；

六、总结：

  首先从 Folyd算法的背景小故事 讲起，弗洛伊德 的传奇经历属实非同凡响。

  然后讲解了 Folyd算法的简介 和 具体的实现过程，并搭配了详细的图例，最后用两种语言进行测试。最后简单列了一下算法的特点。

  整个Folyd算法里面，最重要的就是要理解那 两个状态转移方程，为此我再列一下：

  $在新解封一个结点k后;$

  $if(D[i][k]+D[k][j]<D[i][j])$
   { \{ {

    $D[i][j]=D[i][k]+D[k][j]\text{; //更新矩阵}$
    $path[i][j]=path[k][j]\text{; //更新路径}$
  $\}$

七、部分参考附录：

[1] 《总结一下最短路径的弗洛伊德算法（Floyd）》
链接: 数学建模——详解弗洛伊德(Folyd)算法【分别用 C/C++ 和 matlab 实现】.

[2] 《最短路径模板+解析——（FLoyd算法）》
链接: 数学建模——详解弗洛伊德(Folyd)算法【分别用 C/C++ 和 matlab 实现】.

数学建模系列文章——总结篇：《数模美一国一退役选手的经验分享[2021纪念版]》.