一种针对复杂域中具有Neumann边界条件的反应扩散系统的高效高阶时间积分虚拟域谱方法

《Computers & Mathematics with Applications》：An efficient fictitious-domain spectral method with high-order time integration for reaction-diffusion systems with Neumann boundary conditions in complex domains

　　

**摘要**

我们针对定义在具有曲边界的单连通平面域上、满足齐次Neumann边界条件的反应扩散系统，发展了一种虚拟域谱方法。通过将物理域嵌入一个辅助圆盘中，并利用极坐标变量替换及角向的Fourier展开，将每个时间步简化为一族一维径向修正Helmholtz问题。我们将每个径向问题分解为齐次部分和特解部分。特解部分采用依赖于模式的Legendre-Galerkin基进行离散化，该基在极点处满足极点条件，并通过截断奇异值分解稳定化的边界最小二乘公式，从物理Neumann条件中恢复未知的虚拟边界数据。对于具有交叉扩散的反应扩散系统，我们通过实Schur分解将扩散矩阵化为上三角形式，从而得到顺序解耦的标量问题，使得在无需修改径向离散化的情况下可以复用同一空间求解器。对标量和耦合测试问题的数值实验验证了广义BDF方案在理论上预测的时间精度，而对Schnakenberg模型和CIMA模型在圆盘及花瓣形域上的模拟则表明该方法能够捕捉几何诱导的模式选择以及交叉扩散驱动的模式转换。

**引言**

反应扩散系统为化学、生物学和材料科学中模式形成、输运及非线性时空动力学的建模提供了一个经典框架。自Turing关于扩散驱动不稳定性开创性工作[1]以来，此类系统的精确数值模拟一直是科学计算中的核心挑战，尤其是当模式在有曲边或非平凡边界的有界域上演化时[2]。当需要高阶空间精度、稳健的时间步进以及齐次Neumann边界条件的忠实施加时，这一挑战更为突出——这些挑战已在数值偏微分方程领域引起持续关注[3]、[4]。

两个典型模型驱动了本研究。第一个是Schnakenberg系统[5]，它作为最简自催化模型被提出，广泛用于研究扩散驱动的模式形成。第二个是CIMA（氯酸盐-碘-丙二酸）模型，在该模型中De Kepper等人[6]首次实验观测到了Turing模式。其简化的二组分形式——Lengyel-Epstein模型[7]、[8]——此后已成为模式选择和模态相互作用的标准基准。在两种情形中，计算域的几何形状强烈影响着观测到的形态、起始时间及斑点或条纹的最终排列方式。这些几何效应在非圆形域上尤为显著，此时传统的张量积谱离散化方法不再直接适用。

对于规则的笛卡尔域或张量积域，标准的有限差分法、有限元法和谱方法均已成熟[3]。然而，在不规则域上，真正的高阶谱方法要难得多。这一困难源于边界条件必须在弯曲的物理边界上施加，而模态表示最自然地建立在规则的嵌入域上。有限元方法在处理复杂几何时具有灵活性，但其精度通常受限于代数阶收敛[9]。尽管p型细化或自适应网格细化可以提高精度，但这些方法仍无法与谱方法在光滑解上的效率相匹配[10]。高精度与复杂域适应性之间固有的取舍是计算数学中长期存在的挑战。

为解决这一张力，已发展出两种主要方法来应对该挑战。第一种是域分解策略，它将复杂域划分为多个规则子域，并在每个子域内应用谱方法，同时施加适当的界面耦合[11]、[12]。虽然有效，但此类方法可能在耦合条件的实现上引入额外复杂性，特别是在处理不可分离的椭圆算子或各向异性扩散过程时[13]、[14]。第二种方法聚焦于直接为复杂域构造谱方法。映射技术通过光滑的显式映射将物理域变换为标准计算域[15]、[16]。然而，为任意复杂几何构造这样的映射往往很困难，且对于高阶算子，所得的高阶算子系数可能变得繁琐。现代计算框架促进了这些复杂技术的实现，使得谱方法在实用应用中更易获取[17]。

虚拟域谱方法[18]、[19]、[20]将不规则域嵌入一个更大的、简单形状的辅助域中；然后在该辅助域上应用快速的标准谱算法，同时通过附加约束来施加物理边界条件。对于涉及复杂微分算子的问题，具有协调性质的专门有限元离散化方法也显示出前景，但其扩展到谱方法领域仍是活跃的研究领域[21]、[22]。Gu和Shen为椭圆方程[18]、[19]及外部Helmholtz问题[20]发展了一种系统的、计算高效的且数学上适定的虚拟域谱方法，为此类方法提供了坚实的理论基础。最近，Shen和Wen[23]将该圆形嵌入策略扩展到了四阶问题——包括广义Stokes问题和Cahn-Hilliard方程——而Li和Sun[24]则利用Huang和Shen[25]的广义BDF/IMEX框架将其扩展到了含时Sobolev型方程。Li和An[26]为极坐标几何发展的降维思想进一步为反应扩散系统的高效径向求解器提供了自然途径。

本工作的目的是为曲平面域上的反应扩散及交叉扩散系统发展一套连贯的、可复现的高阶框架。我们的重点主要在于算法：我们旨在构造一个对于长时间模式计算足够精确、从基本原理即可实现且足够透明的公式，并具有支持后续关于稳定性阈值、模式选择和分支形成等解析研究的结构。虽然完整非线性方案的完整收敛性理论超出了本研究的范围，但该数值公式旨在便于在明确的计算框架基础上开展此类分析。

主要贡献总结如下：

（i）通过将物理域嵌入辅助圆盘并将每个时间步简化为径向修正Helmholtz子问题，我们建立了具有齐次Neumann边界条件的单连通复杂域上反应扩散系统的虚拟域谱求解器。

（ii）我们为三种Fourier模式类$m=0$、$|m|=1$和$|m|\geq 2$构造了依赖于模式的Legendre-Galerkin基，确保极点处的极点条件和特解的边界约束被纳入试探空间[26]、[27]、[28]。

（iii）我们通过边界最小二乘公式[19]恢复未知的虚拟边界数据，并利用截断奇异值分解（TSVD）[24]稳定所得的病态系统。

（iv）通过扩散矩阵的实Schur分解，我们将框架扩展到具有正实谱的交叉扩散系统。该方法得到顺序解耦的标量问题，使得在无需修改径向离散化的情况下可以复用空间求解器。

本文其余部分组织如下。第2节介绍模型问题和嵌入域重新表述。第3节发展广义BDF时间离散化、极坐标-Fourier降维、Legendre-Galerkin径向求解器以及基于TSVD的边界重构。第4节呈现时间验证测试以及Schnakenberg和CIMA模型有无交叉扩散的模式模拟，并将解释性评论直接纳入数值讨论中。第5节对论文进行总结。

**章节片段**

**反应扩散系统**

我们考虑以下定义在$\Omega$上的反应扩散系统：

$$\begin{cases} \dfrac{\partial \mathbf{V}(x,y,t)}{\partial t} = \mathbf{D}\Delta \mathbf{V}(x,y,t) + \mathbf{F}(\mathbf{V}), & (x,y)\in\Omega,\; t\in(0,T], \\ \mathbf{V}(x,y,0) = \mathbf{V}_0(x,y), & (x,y)\in\Omega, \\ \dfrac{\partial \mathbf{V}(x,y,t)}{\partial \mathbf{n}} = \mathbf{G}(x,y,t), & (x,y)\in\partial\Omega,\; t\in(0,T], \end{cases}$$

其中$\mathbf{V}=[V_1,\cdots,V_N]^\top$，$\mathbf{F}(\mathbf{V})$表示非线性反应动力学，$\mathbf{D}\in\mathbb{R}^{N\times N}$是扩散矩阵，$\mathbf{n}$是$\Gamma$上的外法线单位向量。设$\Omega\subset\mathbb{R}^2$是具有光滑边界$\partial\Omega$的有界域，令$\tilde{\Omega}$为包含$\Omega$的圆形域，即$\Omega\Subset\tilde{\Omega}$（见图1）。此外，$\mathbf{G}$是给定的Neumann数据；$\mathbf{G}\equiv 0$对应于

**广义BDF时间离散化**

在Huang和Shen[25]的广义BDF公式中引入正的插值参数$\kappa>1$。在均匀网格$t_n=n\tau$上，时间导数和步外解值由下式逼近：

$$\frac{1}{\tau}\sum_{q=0}^{k} a_{k,q}(\kappa) U_{n+1-k+q} = \partial_t U_{n+\kappa} + \mathcal{O}(\tau^k),$$

$$\sum_{q=0}^{k-1} b_{k,q}(\kappa) U_{n+2-k+q} = U_{n+\kappa} + \mathcal{O}(\tau^k),$$

$$\sum_{q=0}^{k-1} c_{k,q}(\kappa) U_{n+1-k+q} = U_{n+\kappa} + \mathcal{O}(\tau^k),$$

其中系数$a_{k,q}(\kappa)$、$b_{k,q}(\kappa)$和$c_{k,q}(\kappa)$由求解相应的Vandermonde线性方程组唯一确定。将广义BDF-$k$方案(6)应用于离散化(3)，

**数值实验**

相应地，在例1–3中，给定的非齐次Neumann数据为

$$\left.\frac{\partial u}{\partial \mathbf{n}}\right|_{\Gamma} = g(\theta,t), \quad g(\theta,t) := \mathbf{J}_u e^{\mathrm{ex}\,\mathrm{ac}\,t} \cdot \mathbf{n}\Big|_{r=\rho(\theta)},$$

边界重构按照注3进行。根据该注，数据仅通过最小二乘系统的右端项$G$进入；矩阵$\mathbf{A}_{\mathrm{LS}}$、其TSVD正则化、依赖于模式的Legendre-Galerkin径向求解器以及广义BDF时间离散化与齐次情形完全相同。因此这些测试

**结论**

我们为具有齐次Neumann边界条件的单连通平面域上的反应扩散和交叉扩散系统发展了一种虚拟域谱方法。通过将物理域嵌入辅助圆盘并结合极坐标-Fourier降维与齐次-特解分解，该方法将原始二维问题简化为一系列一维径向修正Helmholtz方程。依赖于模式的Legendre-Galerkin基

Shan Li | Tianyu Yang

上海理工大学基础科学学部数学系，上海，200093，中国