渗流模拟异常：程序在特定递归深度下无提示终止

我编写了一个基于渗流（Percolation）的相变模拟小程序，原理可参考经典渗流数学讲解视频，建议观看前5分钟理解代码逻辑。

代码逻辑说明 #

• 二维数组main（命名待优化）初始化为全0，每个索引对应网格中的一个节点
• group函数负责为每个节点分配所属“群组”（对应相态/颜色）：检查节点间的连接是否激活，若节点未加入群组则将其纳入，并对该节点重复此遍历操作
• 群组由main中的唯一数值标识，节点间的水平、垂直连接分别由connections_horizontal和connections_vertical两个随机数组表示（连接激活条件为数组值大于参数p）

异常现象 #

当p>0.5时代码运行完全正常，但p<0.5时出现异常：

• 修改sys.setrecursionlimit前，低p值会触发**栈溢出（Stack Overflow）**错误直接崩溃
• 修改递归深度上限为size²后，低p值下程序无任何Python层面报错，直接终止运行

推测是Python在特定递归深度下被系统强制终止，现寻求异常原因分析。

完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import sys

np.random.seed(0)

size=200
x=size**2
sys.setrecursionlimit(x)
main=np.zeros((size,size)) #where everything happens, every index in main is a point on the plane
connections_horizontal=np.random.rand(size-1,size)  #horizontal connections between every points in main
connections_vertical=np.random.rand(size,size-1)  #vertical connections between every points in main
p=0.5 #parameter of percolation

def group(array,index,grp_nbr,p): #index of point for which we want to find the group, the group number (to distinguish the groups), parameter p (percolation)

    array[index]=grp_nbr
    # 1: find if the neighbouring connection is active
    # 2. if the connection is active, put the neighbouring point in the same group
    # 3. for every new point, repeat the procedure until either no connections, or every neighbour is in the same group

    if index[0]<size-1: #checking connections to the right
        if connections_horizontal[index]>p: #if connection is active
            if array[index[0]+1,index[1]]!=grp_nbr: #and neighbour is not in the same group
                array[index[0]+1,index[1]]=grp_nbr #adding the point to the group
                group(array,(index[0]+1,index[1]),grp_nbr,p) #repeat procedure for new point (recursion)
    if index[1]<size-1: #checking the connections up
        if connections_vertical[index]>p:
            if array[index[0],index[1]+1]!=grp_nbr:
                array[index[0],index[1]+1]=grp_nbr
                group(array,(index[0],index[1]+1),grp_nbr,p)
    if index[0]>0: #checking the connections to the left
        if connections_horizontal[index[0]-1,index[1]]>p:
            if array[index[0]-1,index[1]]!=grp_nbr:
                array[index[0]-1,index[1]]=grp_nbr
                group(array,(index[0]-1,index[1]),grp_nbr,p)
    if index[1]>0: #checking the connections down
        if connections_vertical[index[0],index[1]-1]>p:
            if array[index[0],index[1]-1]!=grp_nbr:
                array[index[0],index[1]-1]=grp_nbr
                group(array,(index[0],index[1]-1),grp_nbr,p)


for i in range(size): #finding a group for every point on the plane
    for j in range(size):
        if main[i,j]==0:
            group(main,(i,j),np.random.rand(1),p)

plt.imshow(main,cmap='rainbow',interpolation='none')
plt.savefig('test.png')
plt.show()

异常原因分析 #

1. Python递归限制的本质：sys.setrecursionlimit()仅设置了Python解释器的逻辑栈深度上限，但操作系统会为每个进程分配固定大小的物理栈内存（通常为几MB）。每个递归调用会在物理栈上分配一个栈帧（包含局部变量、返回地址等），当递归深度超过物理栈的承载能力时，会触发操作系统级别的段错误（Segmentation Fault），此时Python解释器会直接崩溃，不会抛出任何Python层面的异常，表现为“无报错直接停止”。
2. 低p值的场景特性：当p<0.5时，渗流系统中会出现超大尺寸的连通团簇（甚至贯穿整个200×200的网格），递归遍历这类团簇的深度会达到数万级（size²=40000），远超操作系统物理栈的承载上限。
3. 高p值正常的原因：p>0.5时，连通团簇的尺寸极小，递归深度较浅，不会触及物理栈的限制，因此运行正常。

修复方案建议 #

• 替换递归实现：将递归式的深度优先搜索改为迭代式BFS/DFS，用手动维护的栈或队列管理遍历节点，完全避开递归栈的限制。
• 使用并查集（Union-Find）：这是渗流问题的标准高效解法，时间复杂度更低，且不存在递归栈的问题。

内容的提问来源于stack exchange，提问作者Ilya Iakoub