Gray-Scott模拟生成Turing patterns时无法对称扩展的问题求助

问题概况 #

编写的Gray-Scott模拟程序生成Turing模式时存在以下异常：

• 模式无法对称扩展，始终偏向右下角
• 仅能在kill rate和feed rate均为0.192时运行
• 模式存在无法定位的随机性

用户代码如下：

def Laplace(loc0,loc1,loc2,loc3,loc4,loc5,loc6,loc7,loc8):
    lap = ((loc0*0.05) + (loc1*0.2) + (loc2*0.05) + (loc3*0.2) + (loc4*-1) + (loc5*0.2) + (loc6*0.05) + (loc7*0.2) + (loc8*0.05))
    return lap

for reps in range(epochs):
    for i in range(height-1):
        for j in range(width-1):
            
            aChem = aConc[i, j] + (aDiffRate * (Laplace(aConc[i-1,j-1],aConc[i-1,j],aConc[i-1,j+1],aConc[i,j-1],aConc[i,j],aConc[i,j+1],aConc[i+1,j-1],aConc[i+1,j],aConc[i+1,j+1])) - aConc[i, j]* (bConc[i,j]**2)) + feedRate * (1 - aConc[i, j]) * dT
            bChem = bConc[i,j] + (bDiffRate * (Laplace(bConc[i-1,j-1],bConc[i-1,j],bConc[i-1,j+1],bConc[i,j-1],bConc[i,j],bConc[i,j+1],bConc[i+1,j-1],bConc[i+1,j],bConc[i+1,j+1])) + aConc[i, j]* (bConc[i,j]*2)) - (killRate + feedRate)* bConc[i,j] * dT

            aConc[i, j] = max(0, min(1, aChem))
            bConc[i, j] = max(0, min(1, bChem))

核心问题与修复方案 #

1. Gray-Scott方程反应项错误 #

标准Gray-Scott模型中，b物质的反应项为a*b²，但代码中写成了aConc[i,j]*(bConc[i,j]*2)（线性项），这直接扭曲了反应动力学，导致仅特定参数能勉强运行，模式行为异常。

修复：将bChem中的+ aConc[i, j]* (bConc[i,j]*2)改为+ aConc[i, j] * (bConc[i,j] ** 2)

2. 原地更新导致顺序依赖破坏对称性 #

直接在原数组aConc和bConc上更新浓度，计算后续细胞时会用到前面已更新的值，引入左上到右下的方向偏倚，导致模式往右下偏移。

修复：创建临时数组存储新浓度，全部计算完成后再替换原数组：

for reps in range(epochs):
    # 创建临时数组复制当前浓度
    new_a = aConc.copy()
    new_b = bConc.copy()
    
    for i in range(height):
        for j in range(width):
            # 周期性边界处理，避免越界并保证对称
            i_prev = (i - 1) % height
            i_next = (i + 1) % height
            j_prev = (j - 1) % width
            j_next = (j + 1) % width
            
            # 计算Laplace值
            lap_a = Laplace(
                aConc[i_prev, j_prev], aConc[i_prev, j], aConc[i_prev, j_next],
                aConc[i, j_prev], aConc[i, j], aConc[i, j_next],
                aConc[i_next, j_prev], aConc[i_next, j], aConc[i_next, j_next]
            )
            lap_b = Laplace(
                bConc[i_prev, j_prev], bConc[i_prev, j], bConc[i_prev, j_next],
                bConc[i, j_prev], bConc[i, j], bConc[i, j_next],
                bConc[i_next, j_prev], bConc[i_next, j], bConc[i_next, j_next]
            )
            
            # 正确的Gray-Scott方程
            new_a[i,j] = aConc[i,j] + (aDiffRate * lap_a - aConc[i,j] * bConc[i,j]**2 + feedRate * (1 - aConc[i,j])) * dT
            new_b[i,j] = bConc[i,j] + (bDiffRate * lap_b + aConc[i,j] * bConc[i,j]**2 - (killRate + feedRate)*bConc[i,j]) * dT
            
            # 限制浓度在合理范围
            new_a[i,j] = max(0, min(1, new_a[i,j]))
            new_b[i,j] = max(0, min(1, new_b[i,j]))
    
    # 更新原数组
    aConc = new_a
    bConc = new_b

3. 循环范围未覆盖全部细胞 #

原循环range(height-1)和range(width-1)导致最后一行、最后一列的细胞完全未更新，这些区域保持初始值，吸引模式往右下偏移。

修复：循环改为range(height)和range(width)，配合周期性边界处理避免越界。

4. 边界条件缺失 #

边缘细胞的Laplace计算需要对称边界条件，否则会出现非对称边缘效应。使用周期性边界（边缘邻居为对面边缘细胞）可保证全局对称性，这也是参考模拟的标准设置。

其他优化建议 #

• Laplace函数可简化为直接在循环中计算，或使用numpy卷积操作提升效率，减少参数传递错误
• 初始条件建议使用中心对称扰动（如中心小区域b浓度较高），更易生成对称模式；随机初始扰动需配合对称边界条件

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