如何修复Python生成二维隐式对流扩散矩阵时出现的索引错误

代码错误原因与修正方案 #

核心错误点 #

1. 遍历范围缺失
   你当前设置range(Nx-1)和range(Ny-1)，遍历的i最大值为Nx-2、j最大值为Ny-2，自然不会处理最后一列和最后一行的网格点，导致矩阵对应行全零。
2. 网格索引计算错误
   二维网格按行展开时，第j行第i列的全局索引应为 k = j * Nx + i，你写成了k = j * Ny + i，当Nx≠Ny时直接触发索引越界，即使Nx=Ny，后续边界访问也会出错。
3. 边界处非法索引访问
   你在处理边界节点时，直接访问了边界外的索引：比如j=0（第一行节点）时访问k-Nx，得到的是负数索引；j=Ny-1（最后一行节点）时访问k+Nx，索引超过矩阵规模；这些操作是你直接改range(Nx)、range(Ny)后报错的直接原因。
4. 步长计算错误
   x方向步长dx应为width / Nx，y方向步长dy应为height / Ny，你写反了计算逻辑，得到的步长数值完全错误。
5. 角落节点边界条件逻辑错误
   你使用elif判断边界，会导致同时属于两个边界的角落节点（比如i=0且j=0的左上角节点）只触发第一个匹配的边界分支，无法正确叠加两个方向的边界条件。

修正步骤 #

1. 调整遍历范围为全网格：

for i in range(Nx):
    for j in range(Ny):

2. 修正全局索引计算：

k = j * Nx + i

3. 修正步长计算：

dx = width / Nx
dy = height / Ny

4. 重构边界判断逻辑，先初始化四个方向的系数，再根据边界状态调整，避免非法索引访问：

# 先初始化内部节点的系数
A_e = D*dx - U_x/2
A_w = D*dx + U_x/2
A_n = D*dy + U_y/2
A_s = D*dy - U_y/2

# 处理x方向边界
if i == 0: # 左边界
    A_w = 0
    if left_bc_type == 0:
        B[k] = left_bc_value * (D*dx + U_x/2) # 狄利克雷边界修正源项
elif i == Nx-1: # 右边界
    A_e = 0
    if right_bc_type == 0:
        B[k] = right_bc_value * (D*dx - U_x/2)

# 处理y方向边界
if j == 0: # 上边界
    A_n = 0
    if top_bc_type == 0:
        B[k] = top_bc_value * (D*dy + U_y/2)
elif j == Ny-1: # 下边界
    A_s = 0
    if bottom_bc_type == 0:
        B[k] = bottom_bc_value * (D*dy - U_y/2)

# 计算主对角线系数
A_p = A_e + A_w + A_n + A_s
Amat[k, k] = A_p

# 仅当系数不为0时赋值相邻索引，避免越界
if A_e != 0:
    Amat[k, k+1] = -A_e # 可根据你的离散格式调整符号
if A_w != 0:
    Amat[k, k-1] = -A_w
if A_n != 0:
    Amat[k, k-Nx] = -A_n
if A_s != 0:
    Amat[k, k+Nx] = -A_s

5. 注意矩阵系数的符号，大部分有限体积离散的对流扩散矩阵，非主对角线系数为负值，你可根据自己的离散格式调整。

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