如何在Reaction Diffusion Simulation中实现2D Laplacian算子

二维拉普拉斯是作用在二维标量场上的二阶微分算子，核心衡量的是某个点的数值和其周围邻域加权平均值的差异。在反应扩散模拟场景下，它用来描述化学物质的扩散效应：如果某点浓度高于周围的加权平均水平，拉普拉斯计算结果为负，物质会向外扩散；反之拉普拉斯结果为正，物质会向该点聚集。

你推导的计算式完全正确，对应的就是Karl Sims反应扩散方案中用到的3x3离散拉普拉斯卷积核，权重分布如下：

[
 [0.05, 0.2, 0.05],
 [0.2,  -1,  0.2],
 [0.05, 0.2, 0.05]
]

该卷积核所有权重相加和为0，符合拉普拉斯算子的基本特性：如果某个区域所有点的数值完全相等，卷积输出结果为0，代表没有浓度差、不存在扩散趋势。

实现注意事项 #

• 边界处理：网格边缘的点没有完整的3x3邻域，通常有三种常用处理方式：
  1. 边界固定：固定边缘所有点的数值为常量（比如0）
  2. 循环边界：网格上下、左右连通形成环面，是反应扩散模拟最常用的方案，能避免边界截断导致的图案断裂
  3. 边缘复制：用边缘像素的值填充邻域外的空缺位置
• 计算一致性：必须基于原始网格的数值计算所有点的拉普拉斯结果，全程不能修改原网格的数值，全部计算完成后再统一更新，否则会覆盖后续点计算需要的邻域原始值，导致结果错误。

拉普拉斯算子伪代码 #

假设输入是二维浓度网格grid，尺寸为width * height，采用循环边界处理：

def laplacian(grid, width, height):
    # 初始化结果数组
    result = [[0.0 for _ in range(height)] for _ in range(width)]
    # 遍历网格所有点
    for x in range(width):
        for y in range(height):
            # 计算循环边界下的邻域坐标
            left = (x - 1 + width) % width
            right = (x + 1) % width
            up = (y - 1 + height) % height
            down = (y + 1) % height
            
            # 提取3x3邻域的9个值
            a = grid[left][up]    # 左上
            b = grid[x][up]       # 正上
            c = grid[right][up]   # 右上
            d = grid[left][y]     # 正左
            e = grid[x][y]        # 中心
            f = grid[right][y]    # 正右
            g = grid[left][down]  # 左下
            h = grid[x][down]     # 正下
            i = grid[right][down] # 右下
            
            # 代入公式计算
            result[x][y] = 0.2 * (b + f + h + d) + 0.05 * (a + c + i + g) - e
    return result

反应扩散迭代调用示例 #

以常用的Gray-Scott反应扩散模型为例，单次迭代的调用逻辑如下：

def reaction_diffusion_step(gridA, gridB, diffA, diffB, feed_rate, kill_rate, time_step, width, height):
    # 分别计算两种物质的拉普拉斯结果
    lapA = laplacian(gridA, width, height)
    lapB = laplacian(gridB, width, height)
    # 临时数组存储更新后的浓度
    newA = [[0.0 for _ in range(height)] for _ in range(width)]
    newB = [[0.0 for _ in range(height)] for _ in range(width)]
    
    for x in range(width):
        for y in range(height):
            a = gridA[x][y]
            b = gridB[x][y]
            # 代入Gray-Scott方程更新浓度
            newA[x][y] = a + (diffA * lapA[x][y] - a * b * b + feed_rate * (1 - a)) * time_step
            newB[x][y] = b + (diffB * lapB[x][y] + a * b * b - (feed_rate + kill_rate) * b) * time_step
    
    return newA, newB

内容的提问来源于stack exchange，提问作者Phineas Greene