基于有限差分法在Python中实现二维扩散方程Neumann边界条件

给定扩散方程：
$$\frac{\partial S}{\partial t} = c \Big(\frac{\partial^2 S}{\partial x^2} + \frac{\partial^2 S}{\partial y^2} \Big)$$
及齐次Neumann边界条件：
$$\frac{\partial S}{\partial \nu} = 0,$$
其中$\nu$为边界外法向。使用有限差分法求解时出现结果偏向右侧边界的问题，核心原因是边界条件的实现方式繁琐，存在潜在的索引逻辑漏洞，同时初始化的中心位置未严格对齐坐标原点，加剧了不对称性。以下是具体分析和修正方案：

原代码核心问题 #

1. 边界条件实现冗余且易出错：原代码通过大量分支判断处理边界和角点，逻辑复杂，容易出现索引或差分格式的隐性错误，导致边界处的扩散行为不符合Neumann条件。
2. 初始化位置未严格对齐原点：直接使用nx//2和ny//2作为初始点索引，由于arange生成的数组并非严格对称，初始点偏离了坐标中心，进一步放大了扩散的不对称性。

修正后的代码 #

采用虚拟边界法实现Neumann条件，逻辑更简洁清晰，避免分支判断的漏洞，同时精准定位初始点：

import numpy as np
import matplotlib.pyplot as plt

dt = 0.1
D = 0.6
dx = np.sqrt(D * dt / 0.1)
dy = dx
t = np.arange(0, 100 + dt, dt)
nt = len(t)
x = np.arange(-15, 15 + dx, dx)
nx = len(x)
y = np.arange(-20, 20 + dy, dy)
ny = len(y)

X, Y = np.meshgrid(x, y)

S = np.zeros((nt, nx, ny))
# 精准定位坐标原点对应的索引
x0_idx = np.argmin(np.abs(x - 0))
y0_idx = np.argmin(np.abs(y - 0))
S[0, x0_idx, y0_idx] = 100

for k in range(nt - 1):
    # 创建带虚拟边界的数组，方便统一处理边界条件
    S_pad = np.zeros((nx + 2, ny + 2))
    S_pad[1:-1, 1:-1] = S[k]
    
    # 应用齐次Neumann边界条件：虚拟点值等于对称内部点值
    S_pad[0, 1:-1] = S_pad[2, 1:-1]    # 左边界
    S_pad[-1, 1:-1] = S_pad[-3, 1:-1]  # 右边界
    S_pad[1:-1, 0] = S_pad[1:-1, 2]    # 下边界
    S_pad[1:-1, -1] = S_pad[1:-1, -3]  # 上边界
    # 处理四个角点
    S_pad[0, 0] = S_pad[2, 2]
    S_pad[0, -1] = S_pad[2, -3]
    S_pad[-1, 0] = S_pad[-3, 2]
    S_pad[-1, -1] = S_pad[-3, -3]
    
    # 统一计算所有点的二阶导数
    S_xx = (S_pad[2:, 1:-1] - 2 * S_pad[1:-1, 1:-1] + S_pad[:-2, 1:-1]) / dx**2
    S_yy = (S_pad[1:-1, 2:] - 2 * S_pad[1:-1, 1:-1] + S_pad[1:-1, :-2]) / dy**2
    
    # 更新下一时间步
    S[k+1] = S[k] + D * dt * (S_xx + S_yy)

plt.contourf(X, Y, S[-1].transpose())
plt.colorbar()
plt.show()

修正说明 #

1. 虚拟边界法：通过在原数组外围添加一层虚拟点，利用Neumann条件的对称性（虚拟点值等于对应内部点值），将所有点的差分计算统一为内部点格式，彻底避免了分支判断的逻辑漏洞。
2. 精准初始点：使用np.argmin(np.abs(x-0))找到最接近坐标原点的索引，确保初始点严格位于区域中心，保证扩散的对称性。

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