无斜率选择分子束外延模型隐显龙格 - 库塔方法的稳定性与收敛性

《Journal of Mathematical Analysis and Applications》：Stability and convergence of implicit explicit Runge-Kutta methods for molecular beam epitaxy model without slope selection

　　

摘要

本文建立了隐式-显式 Runge-Kutta（IERK）方法应用于无斜率选择分子束外延（MBENSS）模型时原始能量稳定性和收敛性的统一证明。利用离散正交核性质和逐阶段 Gr?nwall 不等式，我们推导出了 H1 范数误差估计，且对空间步长与时间步长之间的比例关系没有限制。数值实验证明了 IERK 方法的精度以及原始能量耗散行为。此外，对小扩散系数下粗化动力学的模拟揭示了能量耗散和表面粗糙度及丘宽度幂律增长的标度律，与现有的理论标度律一致。

引言

本文研究无斜率选择的分子束外延模型（MBENSS），也称外延薄膜生长模型，其形式为

?Φ/?t = -ε2Δ2Φ - ?·f(?Φ), x = (x, y) ∈ Ω, 0 < t ≤ T,

初始条件为 Φ(x, y, 0) = Φ?(x, y), (x, y) ∈ Ω ∈ ?2，以及周期边界条件，其中非线性力矢量 f(v) = v/(1 + |v|2) 描述了 Ehrlich-Schwoebel 效应 [11], [24]，ε > 0 是过渡层宽度的常数参数，Φ(·, t) 是随动坐标系中薄膜的标度高度函数，四阶项模拟了表面扩散。众所周知，模型 (1.1) 可视为 Ehrlich-Schwoebel 能量泛函的 L2 梯度流

E[Φ] := ?_Ω (ε2/2)|ΔΦ|2 + G(Φ) dx dy,

其中 G(Φ) := -(1/2)ln(1 + |?Φ|2)。对数项上方界为零，但不存在相对极小值，这意味着 ?Φ 没有能量上偏好的值，从能量角度来看，无法为高度函数 Φ 设定优先斜率。从物理角度来看，这意味着外延生长动力学中缺乏斜率选择机制。更多细节可参阅 [24], [25] 及其参考文献。众所周知，模型 (1.1) 是质量守恒的，即 (Φ, 1) = (Φ?, 1) 对 t > 0 成立，并满足能量耗散律

dE[Φ]/dt = ?_Ω (δE[Φ]/δΦ)(?Φ/?t) dx dy = -?_Ω (?Φ/?t)2 dx dy ≤ 0。

为了控制由非线性项显式离散化引入的可能非线性不稳定性，我们遵循 [10], [14] 引入稳定化参数 κ ≥ 0 和线性稳定化项 κΔΦ，使模型 (1.1) 变为稳定化版本

?Φ/?t = -L_κΦ + g_κ(Φ), (x, y) ∈ Ω,

其中 L_κΦ := ε2Δ2Φ - κΔΦ 和 g_κ(Φ) := -κΔΦ - ?·f(?Φ)。

MBENSS 模型已成为材料科学中的一项有力技术，特别适用于高精度、高纯度化合物半导体的制造。它也被广泛用于研究单晶生长的薄膜沉积 [8], [17], [39]。由于该模型以粗糙-光滑-粗糙转变特征化的多尺度行为，特别是在粗糙表面上外延生长早期阶段，许多工作致力于对模型 (1.1) 的数学和数值分析。数学上，Li 和 Liu [24] 利用摄动分析建立了模型 (1.1) 初边值问题的适定性。在后续工作 [25] 中，他们进一步研究了粗化动力学，推导了关键的物理标度律，包括能量衰减为 O(-ln t)、表面粗糙度增长为 O(√t)、丘宽度标度为 O(t^(1/4))，这些在小 ε 的空间扩展系统中在长时间尺度上被观察到。

大量研究致力于开发模型 (1.1) 的数值方法，包括稳定化半隐式格式 [44]、Crank-Nicolson 型格式 [34]、凸分割方法 [1], [4], [36], [41]、不变能量正交化（IEQ）方法 [16], [45]、标量辅助变量（SAV）方法 [6], [7]、后向差分公式（BDF）方法 [22], [28], [29]、指数时间差分（ETD）格式 [13], [15], [20], [37] 及其参考文献。这些方法的发展显示出显著的进步。Wang 等人 [41] 在 Eyre 的基础工作 [12] 之上建立了一阶凸分割格式。Shen 等人 [36] 随后扩展构建了二阶格式。Chen 等人 [4] 进一步通过引入 Douglas-Dupont 型正则化项增强了该方法，使得非线性系统可以分解为通过基于 FFT 的实现高效求解的线性子问题。这些数值技术在各种相场模型中展示了广泛的应用性。例如，凸分割方法已成功应用于相场晶体模型 [43] 和 Cahn-Hilliard 模型 [27], [46]。同时，基于能量重构的替代方法也已出现，特别是 IEQ 方法 [16], [45] 及其改进变体 SAV 方法 [6], [7]，尽管这些方法保持的是修正形式的离散能量。

近年来，指数时间差分（ETD）方法，特别是多步（ETDMS）[2], [3], [5], [20] 和 ETDRK 或显式指数 Runge-Kutta（EERK）变体 [9], [13], [15], [18], [19], [30], [37]，在数值分析中引起了广泛关注。Ju 等人 [20] 为模型 (1.1) 引入了一类 ETDMS 格式，将空间离散化的 Fourier 谱配点法与使用显式多步近似的 ETD 时间积分相结合。Fu 和 Yang [15] 开发了结构保持的 ETDRK 格式，最高可达二阶，在保持原始能量耗散的同时适用于包括模型 (1.1) 在内的广泛相场模型。事实上，保持梯度流模型原始能量耗散的线性时间步长方法在理论和实践方面都具有很大的吸引力。通过计算相邻时间层之间的能量差，Fu 等人 [13] 推导出了确保 EERK 方法保持原始能量递减的充分条件。Liao 和 Wang [30] 从微分形式的角度，利用一般 EERK 方法的相关微分矩阵推导出了相同的充分条件。[30] 中提出的理论框架不仅适用于 EERK 方法，如果相关微分矩阵是正定的且非线性体 f 是全局 Lipschitz 连续的，也适用于积分因子 Runge-Kutta（CIFRK）和 IERK 方法的修正。

作为与 EERK 方法密切相关的方案，积分因子 RK（IFRK）或所谓的 Lawson 方法 [23] 近年来也被应用于 Allen-Cahn 型梯度流，因为它们具有保持最大界定性原理的明显优势 [21], [26], [48]。然而，[32], [40] 中的数值结果表明 IFRK 方法不能保持梯度流问题的稳态解。为此，[32], [40] 中随后开发了一些 CIFRK 方法。Liao 等人 [32] 提出了两类修正策略，包括望远镜修正（TIFRK）和非线性项平移修正（NIFRK），CIFRK 的修正系数通过强制在所有阶段保持稳态解来确定。Wang 等人 [40] 通过将指数函数和线性函数结合来近似 IFRK 每个阶段中的指数函数，开发了一种新颖的 CIFRK 方法，该方法不仅消除了 IFRK 方法中的指数衰减/增长现象，而且保持了原始稳态。事实上，[40] 中的 CIFRK 方法和 [32] 中的 CIFRK 方法仅在于其构造机制上有所不同，但它们具有相同的数学形式。

EERK 和 CIFRK 方法的计算效率严重依赖于矩阵指数的计算。目前，由于与 EERK 和 CIFRK 方法相比计算复杂度较低，IERK 方法在实际数值模拟中具有自身优势。最近，Teng 等人 [38] 和 Zhang 等人 [47] 建立了梯度流模型中 Radau 型 IERK 方法的误差估计，通过将误差方程的每个阶段重构为两个子阶段，并在 Fourier 伪谱空间中对每个子阶段进行精细的特征值估计。因此，逐阶段分析方法与 IERK 方法的阶段数有关。建立 IERK 方法收敛性的统一理论仍然是一项具有挑战性的工作。幸运的是，Liao 等人 [31] 建立了 Lobatto 型 IERK 方法原始能量耗散的统一框架，如果非线性项是全局 Lipschitz 连续的且相应的微分矩阵是正定的。

本研究的动机来自两个方面。首先，(1.1) 中的非线性体 f 本质上满足 Lipschitz 连续性，使模型 (1.1) 特别适合 [30] 中的理论框架。其次，统一微分形式 (2.4) 和 DOC 形式 (2.5) 与模型 (1.1) 表现出密切的数学相似性。这使我们能够采用更类似于梯度流问题本身的标准化离散能量方法来分析 Runge-Kutta 方法的能量稳定性和收敛性。与 [38], [47] 中的理论分析方法相比，本文提出的方法独立于 Runge-Kutta（RK）阶段的数量，因此更加简洁。

本文其余部分组织如下。在第 2 节中，我们简要描述了 Fourier 空间离散化并介绍了 IERK 方法的统一微分形式和 DOC 形式。然后在第 3 节中，我们提供了 IERK 方法的能量稳定性。作为主要理论部分，我们在第 4 节中利用引理 3.3 中的离散正交卷积核性质和引理 4.2 中的逐阶段 Gr?nwall 不等式建立了 IERK 方法的 H1 误差估计。第 5 节提供了数值结果以验证我们的理论结果。我们最后在本文最后一节中给出一些结论。

Fourier 伪谱方法

我们回顾 Fourier 伪谱空间离散化如下。考虑 2D 空间 Ω = (0, L_x) × (0, L_y)。设 N_x 和 N_y 分别表示 x 和 y 方向上的网格数的两个偶整数。我们在各维度中定义均匀网格尺寸为 h_x = L_x/N_x 和 h_y = L_y/N_y。然后我们获得离散网格 Ω_h = {(jh_x, kh_y) | 0 ≤ j ≤ N_x - 1, 0 ≤ k ≤ N_y - 1} 和 Ω?_h = {(jh_x, kh_y) | 0 ≤ j ≤ N_x, 0 ≤ k ≤ N_y}。

记 L_x × L_y 周期网格函数空间 V_h := {v | v = [v_h(x_j, y_k)], (x_j, y_k) ∈ Ω?_h}。然后对于任意 u, v ∈ V_h 和

能量稳定性

在本节中，我们考察 IERK 方法 (2.1) 的原始能量稳定性以及数值解的有界性。

引理 3.1 [20] f(v) 在 v = (v?, v?)? 处的 Jacobian 矩阵为

?f(v) = (1/(1 + |v|2)2) * (1 - v?2 + v?2, -2v?v?; -2v?v?, 1 - v?2 + v?2)。

?f(v) 的特征值为 μ?(v) = (1 - |v|2)/(1 + |v|2)2，μ?(v) = 1/(1 + |v|2) 满足 -1/8 ≤ μ?(v) ≤ μ?(v) ≤ 1。此外，对于任意 v, u ∈ ?2 有 |f(v) - f(u)| ≤ |v - u| 成立。

引理 3.2 假设稳定化参数 κ ≥ 1/8，对于任意 u_h, v_h ∈ V_h，我们有

?u_h - v_h, g_{κ,h}(v_h) - (1/2)L_{κ,h}(u_h + v_h)? ≤ E[v_h] - E[...

误差估计

我们使用 Sobolev 空间 H^m(Ω) 中 m ≥ 0 的标准半范数和范数。设 C^∞(Ω) 为定义在 Ω 上的无限可微 L_x × L_y 周期函数的集合，H^m_per(Ω) 为 C^∞ 在 H^m(Ω) 中的闭包，配备半范数 |·|_{H^m_per} 和范数 ‖·‖_{H^m_per}。为简便起见，我们记 |·|_{H^m} = |·|_{H^m_per}，‖·‖_{H^m} = ‖·‖_{H^m_per}，‖·‖_{H^0} = ‖·‖_{L2}。子节 2.1 中定义的 L2 投影算子 P_N 和三角插值算子 I_N 的理论结果如下

引理 4.1 [35]

数值实验

首先，我们考虑一参数二阶 IERK 方法，即 IERK2-1，其 Butcher 表如下，见 [31]

隐式和显式部分的系数矩阵分别为

A_I = (a??, 0; 1 - 2a??, a??/√2)，A_E = (√2/2, 0; 2 - √2/2, √2/2)。

根据 (2.3)，我们有

Γ = (√2, 0; 2√2, -2√2)，

D(z, a??) = (√2, -(1/2 + √2 a??)z; 0, 2(√2 - 1 - z + (2 - √2)a?? z)/√2)

和

Θ(z, a??) =
(1/√2, -(1/2 + √2 a??)z;
0, 8(1 - √2 + z + (√2 - 2)a?? z);
2√2, -(1 + 2√2 a??)z;
2, 1/√2 - (1/2 + √2 a??)z)

如 [31] 所示，IERK2-1 的微分矩阵 D(z, a??) 是正的

结束语

本文提出了 Lobatto 型 IERK 方法应用于模型 (1.1) 的原始能量稳定性和收敛性的统一理论框架。本工作的特征如下关键点：

(1) 如果微分矩阵 D(z) 对 z < 0 是正定的，则模型 (1.1) 的 IERK 方法保持原始能量耗散。

(2) 基于 DOC 技术，IERK 方法的误差方程可以重构为 DOC 形式 (2.5)。

(3) 借助引理 3.3 和 Feng Liao...