反应扩散方程二维有限差分格式

Question

我想可视化以下偏微分方程的解：

u_t=u_xx+u_yy+f(u)

在不同时间，例如

u(t=1)

和

u(t=3)

。我使用有限差分方案和以下 Python 代码：

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()

我的目标是拥有这样的东西：

但是我的代码没有给我这样的结果。相反，我得到以下信息：

我不确定我的代码在哪里失败。我不是一个在数值方法方面有丰富经验的人，所以我想知道我可以做些什么来纠正这个结果？

Answer 1

主要问题是数值不稳定。您可以通过增加时间分辨率或降低空间分辨率来解决此问题。我发现，对于常数的组合，将时间分辨率提高 5 倍可以解决发散问题。另外，我认为

axis=

和

u_xx

的

u_yy

论点是错误的。

我已经包含了更准确的二阶中心差异（

np.gradient()

）供您参考，尽管它们没有被使用。我还包括了用于调试的绘图，以及您想要的曲面图。