如何在Python二维反应扩散系统中实现周期性边界条件（PBC）

你当前代码采用的是反射边界，要改成周期性边界，需从拉普拉斯算子计算和数组更新逻辑两个核心部分修改，具体如下：

1. 修改拉普拉斯算子的计算逻辑 #

原laplacian_numba只处理了数组内部（1到N-2）的点，完全忽略边界。周期性边界要求：

• 顶部边界（i=0）的上邻点是底部边界（i=N-1）
• 底部边界（i=N-1）的下邻点是顶部边界（i=0）
• 左侧边界（j=0）的左邻点是右侧边界（j=N-1）
• 右侧边界（j=N-1）的右邻点是左侧边界（j=0）

修改后的拉普拉斯函数需遍历所有点，用模运算实现周期性索引：

@jit(nopython=True, fastmath=True)
def laplacian_numba(field, dh2_inv, out):
    """
    带周期性边界的5点模板拉普拉斯计算
    """
    N = field.shape[0]
    for i in range(N):
        for j in range(N):
            # 用模运算获取周期性邻点索引
            i_up = (i - 1) % N
            i_down = (i + 1) % N
            j_left = (j - 1) % N
            j_right = (j + 1) % N
            
            out[i, j] = (
                field[i_down, j]
                + field[i_up, j]
                + field[i, j_right]
                + field[i, j_left]
                - 4 * field[i, j]
            ) * dh2_inv
    return out

2. 调整数组更新逻辑，移除反射边界代码 #

原update函数中通过a_copy[0,:] = a_copy[1,:]这类操作实现反射边界，现在需完全删除这部分代码。同时，数组更新不能只取内部区域a[1:-1,1:-1]，要直接对整个数组操作：

修改后的update核心片段：

@jit(nopython=True, fastmath=True)
def update(a, b, ..., n):
    # 计算所有点的拉普拉斯（包括边界）
    laplacian_numba(a, dh2_inv, out=lap_a)
    laplacian_numba(b, dh2_inv, out=lap_b)

    # 直接更新整个数组，不再限制内部区域
    a_copy = a + dt * (ODE项 + lap_a * a)
    # 其他物种的更新逻辑同理...

    # 彻底删除原来的反射边界赋值代码

    # 交换数组完成时间步更新
    a, a_copy = a_copy, a
    b, b_copy = b_copy, b  

    return a, b, ..., n  

额外注意事项 #

• 初始化数组时，原代码np.array((N,N), dtype=np.float64)会创建形状为(2,)的错误数组，需改为np.zeros((N,N), dtype=np.float64)或np.random.rand(N,N).astype(np.float64)来生成正确的二维数组。
• 拉普拉斯计算的输出数组lap_a、lap_b需提前初始化好与输入同形状的数组，避免在jit函数内动态创建数组影响性能。

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