2D反应扩散方程有限差分格式可视化结果不符求助

问题描述 #

需要可视化以下偏微分方程的解：

u_t = u_xx + u_yy + f(u)

其中反应项 f(u) = u(1-u)，边界条件为u=0，初始条件是中心高斯分布。预期在t=1、t=3时得到类似目标形态的结果，但当前代码运行结果与预期不符。

现有代码 #

import numpy as np
import matplotlib.pyplot as plt

# Parameters
D = 0.1  # Diffusion coefficient
L = 1.0  # Length of domain
T = 3.0  # Total time
Nx = 100  # Number of grid points in x-direction
Ny = 100  # Number of grid points in y-direction
Nt = 3000  # Number of time steps
dx = L / (Nx - 1)  # Grid spacing in x-direction
dy = L / (Ny - 1)  # Grid spacing in y-direction
dt = T / Nt  # Time step size

# Initial condition
def initial_condition(x, y):
    return np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.1)

# Reaction term
def reaction_term(u):
    return u * (1 - u)

# Set boundary condition (u = 0 at boundary)
def apply_boundary_condition(u):
    u[0, :] = 0  # Bottom boundary
    u[-1, :] = 0  # Top boundary
    u[:, 0] = 0  # Left boundary
    u[:, -1] = 0  # Right boundary

# Create grid
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# Initialize solution array
u = np.zeros((Nx, Ny))

# Set initial condition
u = initial_condition(X, Y)

# Time points to visualize
times_to_visualize = [1, 3]  # Time points at which to visualize solution

# Iterate over time
for n in range(Nt + 1):
    # Compute spatial derivatives using central differences
    u_xx = (np.roll(u, -1, axis=0) - 2*u + np.roll(u, 1, axis=0)) / dx**2
    u_yy = (np.roll(u, -1, axis=1) - 2*u + np.roll(u, 1, axis=1)) / dy**2
    
    # Compute reaction term
    f_u = reaction_term(u)
    
    # Update solution using forward Euler method
    u += dt * (D * (u_xx + u_yy) + f_u)
    
    # Apply boundary condition
    apply_boundary_condition(u)
    
    # Check if current time is in times_to_visualize
    if n * dt in times_to_visualize:
        # Plot solution
        plt.figure()
        plt.contourf(X, Y, u, cmap='coolwarm')
        plt.colorbar(label='u')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title('2D Reaction-Diffusion Equation at t = {:.2f}'.format(n * dt))
        plt.show()

结果对比 #

• 预期结果：
  
• 实际运行结果：
  

修正措施 #

1. 修复空间差分的边界计算错误 #

当前用np.roll计算二阶导数会导致边界点与对面边界循环，完全不符合Dirichlet边界条件（u=0）的物理意义。正确的做法是仅对内部网格点使用中心差分：

# 初始化导数数组
u_xx = np.zeros_like(u)
u_yy = np.zeros_like(u)

# 仅对内部点计算中心差分
u_xx[1:-1, 1:-1] = (u[2:, 1:-1] - 2*u[1:-1, 1:-1] + u[:-2, 1:-1]) / dx**2
u_yy[1:-1, 1:-1] = (u[1:-1, 2:] - 2*u[1:-1, 1:-1] + u[1:-1, :-2]) / dy**2

边界点的导数无需计算，因为后续会强制设置u=0，不参与内部点的更新。

2. 调整扩散系数匹配原方程 #

原方程中没有扩散系数D，但代码里设置了D=0.1，这会大幅降低扩散速度，与预期结果的扩散程度不符。将D改为1，或者直接从更新公式中移除D：

# 更新公式修改为（对应原方程）
u += dt * ((u_xx + u_yy) + f_u)

3. 修复时间点判断的精度问题 #

浮点数运算存在精度误差，直接用n*dt in times_to_visualize可能会错过目标时间点。改用绝对值判断：

current_time = n * dt
for t in times_to_visualize:
    if abs(current_time - t) < 1e-6:
        # 绘图代码
        plt.figure()
        plt.contourf(X, Y, u, cmap='coolwarm')
        plt.colorbar(label='u')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title('2D Reaction-Diffusion Equation at t = {:.2f}'.format(current_time))
        plt.show()
        break

4. 保证数值稳定性（Forward Euler的稳定性条件） #

Forward Euler方法用于2D扩散方程时，必须满足稳定性条件：

dt ≤ (dx² * dy²) / (2*D*(dx² + dy²))

代入当前参数（D=1，dx=dy≈0.0101），计算得最大稳定dt约为2.5e-5。当前dt=0.001远大于该值，会导致数值不稳定。解决方法：

• 增加时间步数Nt，比如设为Nt=120000，使dt≈2.5e-5；
• 改用更稳定的隐式方法（如Crank-Nicolson），但实现复杂度更高。

修正后的示例代码 #

整合以上修改后的代码：

import numpy as np
import matplotlib.pyplot as plt

# Parameters
D = 1.0  # 匹配原方程的扩散系数
L = 1.0  # Length of domain
T = 3.0  # Total time
Nx = 100  # Number of grid points in x-direction
Ny = 100  # Number of grid points in y-direction
Nt = 120000  # 增加时间步数满足稳定性条件
dx = L / (Nx - 1)  # Grid spacing in x-direction
dy = L / (Ny - 1)  # Grid spacing in y-direction
dt = T / Nt  # Time step size

# Initial condition
def initial_condition(x, y):
    return np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.1)

# Reaction term
def reaction_term(u):
    return u * (1 - u)

# Set boundary condition (u = 0 at boundary)
def apply_boundary_condition(u):
    u[0, :] = 0  # Bottom boundary
    u[-1, :] = 0  # Top boundary
    u[:, 0] = 0  # Left boundary
    u[:, -1] = 0  # Right boundary

# Create grid
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# Initialize solution array
u = initial_condition(X, Y)

# Time points to visualize
times_to_visualize = [1, 3]  # Time points at which to visualize solution

# Iterate over time
for n in range(Nt + 1):
    current_time = n * dt
    
    # Compute spatial derivatives using central differences (仅内部点)
    u_xx = np.zeros_like(u)
    u_yy = np.zeros_like(u)
    u_xx[1:-1, 1:-1] = (u[2:, 1:-1] - 2*u[1:-1, 1:-1] + u[:-2, 1:-1]) / dx**2
    u_yy[1:-1, 1:-1] = (u[1:-1, 2:] - 2*u[1:-1, 1:-1] + u[1:-1, :-2]) / dy**2
    
    # Compute reaction term
    f_u = reaction_term(u)
    
    # Update solution using forward Euler method
    u += dt * (D * (u_xx + u_yy) + f_u)
    
    # Apply boundary condition
    apply_boundary_condition(u)
    
    # Check if current time is in times_to_visualize
    for t in times_to_visualize:
        if abs(current_time - t) < 1e-6:
            # Plot solution
            plt.figure()
            plt.contourf(X, Y, u, cmap='coolwarm')
            plt.colorbar(label='u')
            plt.xlabel('x')
            plt.ylabel('y')
            plt.title('2D Reaction-Diffusion Equation at t = {:.2f}'.format(current_time))
            plt.show()
            break

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