基于有限差分法求解反应扩散方程生成图灵图案的问题排查

作为学校项目，我使用Python实现有限差分法以数值求解如下反应扩散PDE系统：

目前已完成如下代码实现：

import numpy as np
import matplotlib.pyplot as plt

D_u = 1.0
D_v = .03
k1 = 5.0
k2 = 11.0
dx = .01
dy = .01
dt = .0001 

# Initialize domain
x = np.arange(0, 3, dx)
y = np.arange(0, 3, dy)

plate_length = len(x)
max_iter_time = 10000

# initialize solution matrix
u = np.zeros((max_iter_time, plate_length, plate_length)) #u(t,x,y) 
v = np.zeros((max_iter_time, plate_length, plate_length)) #v(t,x,y)

# Initial condition : (u,v) = (u_s, v_s) + r
random_perturbation_u = np.random.normal(loc = 0.01, size=(plate_length, plate_length))
random_perturbation_v = np.random.normal(loc = 0.01, size=(plate_length, plate_length))

u[0,:,:] = 1 + (k2**2)/25 + random_perturbation_u
v[0,:,:] = k2/5 + random_perturbation_v

def periodic_bc(matrix):
    matrix[:, 0, :] = matrix[:, -2, :] 
    matrix[:, -1, :] = matrix[:, 1, :]  
    matrix[:, :, 0] = matrix[:, :, -2]  
    matrix[:, :, -1] = matrix[:, :, 1]  

    return matrix

#implement finite difference
def calculate_u_and_v(u,v):
    for k in range(0, max_iter_time-1,1):
        for j in range(1, plate_length-1):
            for i in range(1, plate_length-1): 
                u[k + 1, i, j] = u[k][i][j] + dt * ((D_u * ((u[k][i+1][j] - 2*u[k][i][j]+ u[k][i-1][j])/dx**2 + (u[k][i][j+1] - 2*u[k][i][j] + u[k][i][j-1])/dy**2))
                + k1 * (v[k][i][j] - (u[k][i][j] * v[k][i][j])/(1 + v[k][i][j]**2)))

                v[k + 1, i, j] = v[k][i][j] + dt * ((D_v * ((v[k][i+1][j] - 2*v[k][i][j]+ v[k][i-1][j])/dx**2 + (v[k][i][j+1] - 2*v[k][i][j] + v[k][i][j-1])/dy**2))
                + k2 - v[k][i][j] - (4 * u[k][j][i] * v[k][j][i])/(1 + v[k][i][j]**2))

        u = periodic_bc(u)
        v = periodic_bc(v)
        
    return u,v

u,v = calculate_u_and_v(u,v)

# plot the final image
plt.figure()
plt.pcolormesh(x, y, v[max_iter_time-1], cmap=plt.cm.jet, vmin=np.min(v), vmax=np.max(v))
plt.title(f"t = {max_iter_time}")
plt.ylabel(f"$k_1=${k1}")
plt.xlabel(f"$D_v=${D_v}")
plt.colorbar(ticks = [])

plt.show()

运行后得到均匀分布的结果，未生成预期的斑点状图灵图案，恳请帮忙排查实现中的问题。

问题排查与修正建议 #

• 反应项索引错误：在v的更新公式中，u[k][j][i] * v[k][j][i]的索引顺序颠倒，应该是u[k][i][j] * v[k][i][j]。索引错位会破坏空间对称性，导致无法形成图案。修正后代码行：

  v[k + 1, i, j] = v[k][i][j] + dt * ((D_v * ((v[k][i+1][j] - 2*v[k][i][j]+ v[k][i-1][j])/dx**2 + (v[k][i][j+1] - 2*v[k][i][j] + v[k][i][j-1])/dy**2))
  + k2 - v[k][i][j] - (4 * u[k][i][j] * v[k][i][j])/(1 + v[k][i][j]**2))
  

• 初始扰动幅度不足：当前使用loc=0.01的正态分布扰动，均值非零且幅度太小，难以触发图灵不稳定性。建议改为零均值、更大标准差的扰动：

  random_perturbation_u = np.random.normal(loc=0.0, scale=0.05, size=(plate_length, plate_length))
  random_perturbation_v = np.random.normal(loc=0.0, scale=0.05, size=(plate_length, plate_length))
  

• 迭代时间不足：10000次迭代对应实际时间为10000*0.0001=1.0，系统可能还未演化到稳定图案阶段。可将max_iter_time增大至50000或100000，或在保证稳定性的前提下适当增大dt。

• 周期性边界条件实现错误：当前periodic_bc函数修改整个时间维度的矩阵，这会覆盖未计算的未来时间步数据。应只对当前迭代的时间步切片应用边界条件：

  def apply_periodic_bc(matrix_slice):
      matrix_slice[0, :] = matrix_slice[-2, :] 
      matrix_slice[-1, :] = matrix_slice[1, :]  
      matrix_slice[:, 0] = matrix_slice[:, -2]  
      matrix_slice[:, -1] = matrix_slice[:, 1]  
      return matrix_slice
  

  在循环内替换为：

  u[k+1] = apply_periodic_bc(u[k+1])
  v[k+1] = apply_periodic_bc(v[k+1])
  

• 参数范围验证：确认k1=5.0、k2=11.0、D_v=0.03是否处于图灵空间内。可微调参数（比如将k2调整为10.0，D_v调整为0.02），观察是否能触发图案形成。

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