一种用于扩散-反应问题的初始-边界校正分裂方法

《Computers & Fluids》：An initial–boundary corrected splitting method for diffusion-reaction problems

　　蒂·塔姆·当 | 卢卡斯·艾因克默 | 亚历山大·奥斯特曼
芬兰赫尔辛基大学数学与统计系

**摘要**
斯特朗分裂（Strang splitting）是一种广泛用于求解扩散-反应问题的二阶方法。然而，当边界条件为狄利克雷（Dirichlet）时，其收敛阶数通常降为1；而对于诺伊曼（Neumann）和罗宾（Robin）边界条件，收敛阶数降为1.5，这导致精度降低和计算效率下降。在本文中，我们提出了一种新的分裂方法，称为初始-边界校正分裂（initial–boundary corrected splitting），该方法避免了收敛阶数的降低，同时提高了更广泛应用场景下的计算效率。与文献中提出的校正方法不同，它不需要计算依赖于边界条件和边界数据的校正项。通过严格的分析收敛性和数值实验，我们证明了所提方法在精度和性能上的改进效果。

**1. 引言**
我们的目标是开发一种有效的数值方法来求解扩散-反应问题。这一主题受到了相当多的关注，尤其是在时间积分方案的研究领域。特别是分裂方法（splitting methods）引起了极大的兴趣（参见参考文献[1]、[2]、[3]、[4]、[5]、[6]）。这是因为它们使我们能够分别处理扩散和反应问题。线性扩散问题可以使用快速泊松求解器（如多重网格或势方法）高效求解（参见[7]）。如果子求解器具有该特性，分裂方法还具有保持正值的优势，无论时间步长如何（参见[8]）。

斯特朗分裂方法对于具有简单边界条件（如周期性边界条件）的扩散-反应问题可以实现二阶收敛。然而，一般来说，二阶收敛需要扩散边界条件与反应流之间的兼容性。对于非均匀的狄利克雷或诺伊曼边界条件，这通常是不成立的，从而导致收敛阶数的降低（参见参考文献[10]、[11]、[12]、[13]中的实验）。收敛阶数降低的问题引起了研究人员的广泛关注，因此出现了多种克服这一问题的方法。在一系列论文中，艾因克默和奥斯特曼（Einkemmer & Ostermann）[11]、[12]、[13]、[14]引入了一类针对扩散-反应方程的改进分裂方法。他们的方法通过添加或减去校正项来修改两个子过程，而无需改变边界条件。贝尔托利和维尔马尔（Bertoli & Vilmart）[10]提出了一种不同的方法，其中校正项是根据流出的结果构建的，而不是直接根据非线性项构建的。最后，阿尔onso-Mallo、Cano和Reguera（Alonso-Mallo, Cano & Reguera）[15]通过适当调整边界值来避免收敛阶数的降低。

现有的改进分裂方法通常需要预先计算校正项，如上所述，这增加了实现的复杂性并提高了计算复杂度。在本文中，我们提出了一种初始-边界校正（IBC）分裂方案，用于扩散-反应问题，它可以避免收敛阶数的降低并简化实现。通过从数值解中减去一个平滑的校正项z(t)，可以实现零初始数据和均匀边界条件，从而无需为不同的边界条件预先计算校正项。这使得该方法比现有方法具有更广泛的应用性和更高的效率。

除了介绍我们的新型IBC分裂方法外，我们还对所提出的斯特朗分裂方法进行了严格的理论分析。具体来说，定理4.1在适当的平滑性假设下证明了该方法在分析半群框架内的二阶收敛性。为了展示我们方法的实际效果，我们提供了数值实验，证实了其收敛性，并说明了与现有方法相比的优势。

**2. 问题设置**
本节介绍了我们进行初始-边界校正斯特朗分裂收敛性分析的模型问题。设Ω?Rd是一个有界开集，其边界?Ω是光滑的。我们考虑Ω中的以下扩散-反应问题，并具有斜向时间不变边界条件：
$$
\begin{align*}
\partial_t u(t) &= Du(t) + f(t, u(t)), \\
bu | \partial_{\Omega} &= b, \\
u(0) &= u_0.
\end{align*}
$$
二阶线性椭圆微分算子D定义为
$$
D = \sum_{i,j=1}^d \partial_i(a_{ij}(x) \partial_j u) + \sum_{i=1}^d c_i(x) \partial_j u + d(x) u,
$$
其中矩阵值函数$a_{ij}(x) \in \mathbb{R}^d \times \mathbb{R}^d$被假设为对称的，系数$a_{ij}$, $c_i$, $d$足够光滑。当$a_{ij}(x) = \delta_{ij}$, $b_i = 0$, $c = d = 0$时，算子D是拉普拉斯算子。

如果对于任何$x \in \Omega$, $\sigma \in \mathbb{R}^d$，存在一个正常数$\eta$，使得
$$
\sum_{i,j=1}^d |a_{ij}(x) \sigma_i \sigma_j| \geq \eta \sum_{i=1}^d |\sigma_i|^2,
$$
则称算子D为强椭圆的。对于足够光滑的$\beta_i(x)$和$\alpha(x)$，我们定义边界算子B为
$$
bu = \sum_{i=1}^d \beta_i(x) \partial_j u + \alpha(x) u.
$$
我们假设B满足[13]中描述的均匀非切线条件。系数$\beta_i$和$\alpha$的选择决定了应用的边界条件类型。例如，设置$\alpha = 1$和$\beta_i(x) = \sum_{j=1}^d a_{ij}(x) n_j(x)$（对于$1 \leq i \leq d$，其中$n_j(x)$表示$\Omega$在点$x \in \partial_{\Omega}$处的外法线），则得到罗宾边界条件。

**3. 提出的分裂方法描述**
在本节中，我们介绍了用于求解(2.1)的IBC分裂方法的构建过程。主要思想是将问题转化为初始数据为零且边界数据均匀的形式。这是通过构造一个分段光滑函数$z(t)$来实现的，该函数满足边界条件同时保留初始数据。我们首先描述如何在时间区间$[t_n, t_{n+1]$内进行分裂步骤，步长为$\tau$。设$u_n$是时间$t = t_n$时$u(t)$的数值近似值，$z_n(t)$是所需的校正项。我们简单地设置$z_n(t) = u_n$，因为数值解$u_n$满足该步骤的边界和初始条件，即$z_n(t_n) = u_n$。变换后的变量$u^\hat{n}(t) = u(t) - z_n(t) = u(t) - u_n$是时间区间$[t_n, t_{n+1}$上以下问题的解：
$$
\begin{align*}
\partial_t u^\hat{n} &= Du^\hat{n} + h(t, u^\hat{n}) + g(t), \\
bu^\hat{n} | \partial_{\Omega} &= 0,
\end{align*}
$$
其中$h$是形式为
$$
h(t, u^\hat{n}) = f(t, u^\hat{n} + u_n) - f(t, u),
$$
的修正非线性项，$g(n)$是以下源项：
$$
g(n) = D u_n + f(t, u_n).
$$
现在我们将(3.1)分为两个子问题：
$$
\begin{align*}
\partial_t v^\hat{n} &= D v^\hat{n} + g(n), \\
bu^\hat{n} | \partial_{\Omega} &= 0,
\end{align*}
$$
并对这些问题应用标准的斯特朗分裂，初始值为$v^\hat{n}(t_n) = 0$。这种方法被称为初始-边界校正分裂方案，适用于时间不变边界条件。相关的数值方案如下（算法1）。

**4. 收敛性分析**
本节的目的是对应用于(2.1)的初始-边界校正斯特朗分裂方法的收敛性质进行全面的分析。我们考虑抽象的演化方程
$$
\partial_t u^\hat(t) + A u^\hat(t) = h(t, u^\hat(t)) + g(t),
$$
其中$h$由(3.2)定义，$-A$是巴拿赫空间$X$上的分析半群的下无限生成元，范数为$|| \cdot ||$。A的作用定义为$Aw = -Dw$，适用于所有函数$w \in D(A)$，即A的定义域。这个子空间包括与问题相关的均匀边界条件。例如，对于$X = L^2(\Omega)$和二阶强椭圆算子D，我们有$D(A) = H^2(\Omega) \cap H^{01}(\Omega)$（对于狄利克雷边界条件）。

在有界时间区间$[t, T]$上，分析半群是有界的（4.2）
$$
\| e^{-tA} \| \leq C,
$$
并且具有抛物线平滑性质（4.3）
$$
\| A^\gamma e^{-tA} \| \leq C^{t-\gamma}
$$
对于所有$\gamma > 0$。

应用常数变分公式到(4.1)，我们可以写出(2.1)的精确解为
$$
u(t_{n+1}) = z(n_{n+1}) + e^{-\tau} A u^\hat(t_n) + \int_{0}^{\tau} e^{-(\tau - s)} A g(t_n + s) ds + \int_{0}^{\tau} e^{-(\tau - s)} A h(t_n + s, u^\hat(t_n + s)) ds.
$$
现在我们将(4.1)分为扩散方程(4.5)
$$
\partial_t v^\hat(n) = D v^\hat{n} + g(n),
$$
和非线性反应方程(4.6)
$$
\partial_t w^\hat(n) = h(t, w^\hat(n)).
$$
收敛性分析的总体思路遵循[13]中概述的方法。主要区别在于对数值解施加了零初始数据，这使得可以消除许多项，从而得到简单简洁的证明。

**4.1. 错误分析**
首先，我们研究了将初始-边界校正斯特朗分裂方法应用于(2.1)时的局部误差。为此，我们从时间$t_n$开始的数值解的一个步骤出发。对初始数据应用校正后得到变换后的初始数据$v^\hat(t_n) = u^\hat(n) - u_n$。接下来，我们以$\tau/2$的步长积分(4.5)，得到
$$
V = v^\hat(t_n + \tau^2) = e^{-\tau^2} A(u(t_n) - u_n) + \int_{0}^{\tau^2} e^{-\tau^2 - s} A g(t_n + s) ds.
$$
随后可以将(4.6)的解表示为泰勒级数展开：
$$
w^\hat(t_n + \tau) = V + \tau h(t_n, V) + \tau^2 \left( \partial_1 h(t_n, V) + \partial_2 h(t_n, V) \right) h(t_n, V) + O(\tau^3).
$$
这里，$\partial_1 h$和$\partial_2 h$分别表示$h$对其第一个和第二个参数的偏导数。我们用$O(\tau^3)$表示余项，它是有界的，并且关于$\tau$是三次方的。在整个小节中，这种表示方法是一致的。最后，我们以$\tau/2$的步长积分(4.5)，使用初始值$v^\hat{n + \tau^2} = w^\hat(t_n + \tau)$得到
$$
u^\hat{n+1} = e^{-\tau^2} A w^\hat(t_n + \tau) + \int_{0}^{\tau^2} e^{-\tau^2 - s} A g(t_n + \tau^2 + s) ds.
$$
因此，时间$t_n + 1$时具有初始值$u(t_n)$的(2.1)的数值解为
$$
S_{\tau} u(t_n) = z(n_{n+1}) + u^\hat{n+1} = z(n_{n+1}) + e^{-\tau} A(u(t_n) - u_n) + \tau e^{-\tau^2} A(h(t_n, V) + \int_{0}^{\tau} e^{-(\tau - s)} A g(t_n + s) ds + \tau^2 \left( \partial_1 h(t_n, V) + \partial_2 h(t_n, V) \right) h(t_n, V) + O(\tau^3).
$$
其中我们用$S_{\tau}$表示IBC斯特朗分裂算子。

设局部误差为$\delta_{n+1} = S_{\tau} u(n) - u(t_n + 1)$。从(4.10)中得到的数值解中减去精确解(4.4)，得到
$$
\delta_{n+1} = \tau e^{-\tau^2} A(h(t_n, V) + \tau^2 \left( \partial_1 h(t_n, V) + \partial_2 h(t_n, V) \right) h(t_n, V) - \int_{0}^{\tau} e^{-(\tau - s)} A h(t_n + s, u^\hat(n + s)} ds + O(\tau^3).
$$
对于我们的分析，我们对(2.1)的数据做以下假设：
**假设4.1**
设$\Omega$是具有光滑边界的域。设线性微分算子D是强椭圆的，系数足够光滑，并且受到具有光滑数据的斜向边界条件的约束。进一步假设非线性项$f$和解$u$也足够光滑。

**引理4.1**
引理4.1是对局部误差(4.11)的精细估计，我们在定理4.1的证明中用它来展示全局误差的收敛性。

根据假设4.1，将初始-边界校正斯特朗分裂应用于(2.1)满足局部误差界限(4.12)
$$
\delta_{n+1} = O(\tau^2), \quad \delta_{n+1} = A \delta^\hat{n+1} + O(\tau^3), \quad \delta_{n+1} = O(\tau^3).
$$

**证明**
我们设置
$$
\phi_n(s) = e^{-(\tau - s) A \phi^\hat{n}(s), \quad \phi^\hat{n}(s) = h(t_n + s, u^\hat{n}(t_n + s)).
$$
使用中点法则，得到以下等式：
$$
\int_{0}^{\tau} \phi_n(s) ds = \tau \phi_n(\tau^2) + \int_{0}^{\tau} \int_{\tau^2}^s \phi_n'(\xi) ds = \tau \phi_n(\tau^2) + \frac{1}{2} \int_{0}^{\tau} M(s, \tau) \phi_n''(s) ds,
$$
其中核函数$M(s, \tau) = s^2$如果$s < \tau/2$，否则$M(s, \tau) = (\tau - s)^2$。使用(4.14)，我们可以将局部误差(4.11)重写为：
$$
\delta_{n+1} = \delta_{n+1}^1 + \delta_{n+1}^2,
$$
其中
$$
\delta_{n+1}^1 = \tau e^{-\tau^2} A(h(t_n, V) - \tau e^{-\tau^2} A(h(t_n + \tau^2, u^\hat{n}(t_n + \tau^2)) + \tau^2 \left( \partial_1 h(t_n, V) + \partial_2 h(t_n, V) \right) h(t_n, V),
$$
以及
$$
\delta_{n+1}^2 = -\frac{1}{2} \int_{0}^{\tau} M(s, \tau) \phi_n''(s) ds + O(\tau^3).
$$
首先我们限制$\| \delta_{n+1}^1\|$。使用(4.4)和(4.7)，得到
$$
u^\hat{n(tn + \tau^2) - V = \int_{0}^{\tau^2} e^{-\tau^2 - s} A(h(t_n + s, u^\hat(n + s}) ds,
$$
其中被积函数$\phi_{n'}(s) = e^{-\tau^2 - s} A(h(t_n + s, u^\hat(n + s))$满足等式：
$$
\int_{0}^{\tau^2} \phi_{n'}(s) ds = \tau^2 \phi_{n'}(0) + \int_{0}^{\tau^2} \int_{0}^s \phi_{n'}(\xi) ds = \tau^2 e^{-\tau^2} A(h(t_n, u^\hat(n + s)) + \int_{0}^{\tau^2} (\tau^2 - \xi) \phi_{n'}(\xi) ds.
$$
使用(4.18)和(4.19)，我们得到(4.20)
$$
h(t_n + \tau^2, u^\hat(n + \tau^2)) - h(t_n, V) = \tau^2 \partial_1 h(t_n, V) + \tau^2 \partial_2 h(t_n, V) e^{-\tau^2} A(h(t_n, u^\hat(n + \tau^2)) + \partial_2 h(t_n, V) \int_{0}^{\tau^2} (\tau^2 - \xi) \phi_{n'}(\xi) ds + O(\tau^2).
$$
将(4.20)代入(4.16)并使用$\phi_{n'}(\xi) =为了方便起见，我们将左边界条件和右边界条件分别称为b1和b2。拉普拉斯算子使用标准的二阶梯中有限差分方案进行离散化，该方案包含500个网格点。为了获得参考解，我们采用了Matlab的ODE45求解器，并将绝对容忍度和相对容忍度设置为10^-9。在我们的模拟中，第3节中介绍的初始-边界校正的Strang分裂方法被称为IBC Strang分裂。我们评估了经典Strang分裂方法和IBC Strang分裂方法在不同边界条件（包括Dirichlet边界条件、Neumann边界条件、Robin边界条件和混合边界条件）下对公式(2.1)的应用效果。

示例5.1 Dirichlet边界条件，即α=1, β1=0
在这个例子中，我们将问题(2.1)的左边界设为b1=2，右边界设为b2=3。初始条件定义为u0=2+sinπ2x。相应的数值结果如图1(a)所示。

示例5.2 Neumann边界条件，即α=0, β1=1
这个例子的设定是b1=1和b2=2。Neumann边界条件通过标准的二阶梯中有限差分方法离散化。初始条件被指定为u0=x+x^2/2+cos(πx)，该条件满足给定的Neumann边界条件。数值结果如图1(b)所示。

示例5.3 Robin边界条件，即α=β1=1
我们考虑具有Robin边界条件的问题(2.1)，将b1设为1/2，b2设为1/2+1/π。初始条件选为u0(x)=1/2+1/2π?1/2πcosπx。图1(c)显示了数值结果。

示例5.4 混合边界条件
这个例子考虑了具有Neumann边界条件b1=0和Dirichlet边界条件b2=2的问题。我们选择u0(x)=2?2cos(1/2πx)作为初始条件。在图1(d)中，我们展示了数值结果。
我们观察到，在Dirichlet边界条件和混合边界条件下，经典Strang方案的阶数降低到大约一阶（见(a)、(d)），而在Neumann边界条件和Robin边界条件下，阶数降低到1.5阶（见(b)、(c)）。相比之下，IBC Strang分裂方法无论选择哪种边界条件都能保持二阶收敛，并且即使对于较大的时间步长也能显著提高精度。

下载：下载高分辨率图像（107KB）
下载：下载完整尺寸图像(a)。Dirichlet边界条件。

下载：下载高分辨率图像（115KB）
下载：下载完整尺寸图像(b)。Neumann边界条件。

下载：下载高分辨率图像（110KB）
下载：下载完整尺寸图像(c)。Robin边界条件。

下载：下载高分辨率图像（111KB）
下载：下载完整尺寸图像(d)。混合边界条件。

图1. 我们解决了具有不同类型边界条件的一维扩散-反应方程。通过在t=0.5时将数值解与参考解进行比较，计算离散无穷范数下的绝对误差。

6. 二维空间中的数值结果
在本节中，我们考虑具有非线性反应项f(t,u(t,x,y))=u(t,x,y)^2的扩散-反应问题(2.1)，其中Ω=(0,1)^2，D是拉普拉斯算子的标准二阶梯中有限差分近似，使用50×50的网格点。通过在t=0.1时将数值解与参考解进行比较，评估离散无穷范数下的误差。这个参考解是使用Matlab ODE45求解器获得的，其绝对容忍度和相对容忍度均设置为10^-9。我们展示了在示例6.1和示例6.2中，分别应用经典Strang分裂方法和IBC Strang分裂方法对公式(2.1)的处理效果，考虑了Dirichlet边界条件和混合边界条件。数值结果分别展示在图2(a)和图2(b)中。

示例6.1
在这个例子中，我们规定问题(2.1)满足u|?Ω=1。初始数据选择为u0=1+sin(πx)sin(πy)，以确保符合边界条件。

示例6.2
在这个例子中，我们研究了由公式(2.1)描述的问题，在顶部和底部应用Dirichlet边界条件，在左侧和右侧应用Neumann边界条件。边界数据b的选择确保了初始条件u0(x,y)与指定的混合边界条件一致。初始值指定如下：u0(x,y)=3+e^(-10y)+e^(-12/2)*cos(2π(x+y))。
与一维情况（见示例5.1、示例5.4）一致，结果表明经典Strang分裂方法降低为一阶方法，而IBC Strang分裂方法明显显示出二阶收敛。

下载：下载高分辨率图像（116KB）
下载：下载完整尺寸图像(a)。Dirichlet边界条件。

下载：下载高分辨率图像（114KB）
下载：下载完整尺寸图像(b)。混合边界条件。

图2. 我们解决了具有不同类型边界条件的二维扩散-反应方程。通过在t=0.1时将数值解与参考解进行比较，计算离散无穷范数下的绝对误差。

7. 结论
我们提出了一种初始-边界校正的分裂方法，该方法有效地避免了在应用分裂方法处理具有非平凡边界条件的扩散-反应问题时常见的阶数降低现象。在许多实际应用中，例如燃烧问题，Neumann边界条件特别值得关注。尽管文献中存在能够恢复二阶收敛性的改进分裂方案，但这些方案通常需要相对复杂的修改[13]。相比之下，我们的新分裂方法可以广泛应用于不同的边界条件，无需重新计算修改参数，因此适用范围更广。

CRediT作者贡献声明
Thi Tam Dang：撰写——审阅与编辑，撰写——原始草稿，可视化，验证，软件开发，方法学研究。
Lukas Einkemmer：撰写——审阅与编辑，监督，方法学研究，形式分析，概念化。
Alexander Ostermann：撰写——审阅与编辑，撰写——原始草稿，项目管理，方法学研究，资金筹集，形式分析，概念化。