一种针对二维各向异性空间分数扩散方程的残差控制方法

《Mathematics and Computers in Simulation》：A residual-controlled method for two-dimensional anisotropic space-fractional diffusion equations

　　

• 针对带有偏斜双侧通量的各向异性空间分数阶扩散问题，提出了一种JSPG方法。
• 雅可比分数恒等式导出了对角刚度矩阵和Kronecker结构算子。
• 残差控制的有限维稳定性与误差分析给出了非精确CN格式的显式空间收敛阶。

引言
二维各向异性空间分数阶扩散模型描述了具有方向依赖性的非局部输运。其双侧通量、端点行为以及非精确时间步进系统促使我们寻求与算子结构和加权分析相兼容的离散化方案。分数阶扩散方程用于建模复杂与非均匀介质中的反常扩散和长程输运[1]、[2]。其Riemann–Liouville形式在文献[3]中有描述。相应的积分与微分构成了左、右双侧扩散通量的基础[4]、[5]，其中偏置参数对两侧进行加权。

本文考虑的模型是在矩形域Ω=(x_L, x_R)×(y_L, y_R)上的二维各向异性空间分数阶扩散方程：
???
?_t u ? ?_x(a_x D^{s_x, ?_x}_{x} u) ? ?_y(a_y D^{s_y, ?_y}_{y} u) = f(x, y, t), (x, y) ∈ Ω, 0 < t ≤ T,
u(x, y, 0) = u_0(x, y), (x, y) ∈ Ω,
u(x, y, t) = 0, (x, y) ∈ ?Ω, 0 < t ≤ T.

其中a_x, a_y > 0，s_x, s_y ∈ (0,1)，?_x, ?_y ∈ (0,1)。算子D^{s_x, ?_x}_{x}和D^{s_y, ?_y}_{y}是下文定义的偏斜双侧分数阶通量算子。令δ_x = 2 ? s_x，δ_y = 2 ? s_y，其散度为阶数在(1,2)内的左、右Riemann–Liouville微分的偏斜组合。相关的有限差分格式在文献[6]中有研究，变分格式在文献[7]、[8]中有研究。

有限差分方法在规则网格上逼近Riemann–Liouville型或Riesz型微分[9]，其二维稳定性与收敛性在文献[10]中有研究。其结构化系统促使了Toeplitz表示、ADI迭代、带状预条件子以及快速高阶求解器的提出[11]、[12]、[13]、[14]、[15]、[16]。其他近期的高阶发展包括用于多阶次和分布阶次时间分数阶扩散的基于核的伪谱逼近、用于多维非线性分数阶积分-微分方程的分级网格与紧格式离散化，以及用于三维空间分数阶中性延迟模型的紧格式差分方案[17]、[18]、[19]。

谱方法与算子结构相匹配的逼近空间[20]。分数阶谱配点法和Petrov–Galerkin方法使用与算子适应的基函数[21]、[22]，广义雅可比函数将分数阶微积分与加权端点逼近联系起来[23]、[24]。相关的Fourier–Galerkin和Jacobi–Galerkin离散化出现在文献[25]、[26]中，并启发了下文所使用的试探函数与检验函数。

三个问题相互耦合。左、右Riemann–Liouville分量必须由双侧参考算子表示。不同的试探空间与检验空间需要兼容的质量、刚度和投影算子。每个非精确CN求解的代数残差也必须进入全离散分析。相关的快速求解器和带状预条件子在文献[27]、[28]、[29]中有研究。物理区间被映射到参考区间，其中双侧雅可比恒等式[30]给出对角的各向异性刚度矩阵和Kronecker结构的全局算子。CN格式的分析使用由检验空间Gram矩阵导出的对偶残差范数，采用抛物型问题的经典Galerkin框架[31]。

实验评估了与算子对齐和非对齐端点配置下的空间误差、时间收敛性和CPU时间。在所列离散化下，与文献[10]、[13]的比较表明，本文方法报告了更小的误差和CPU时间。

本文包含两个主要部分。首先，算子适应的一维雅可比恒等式通过张量积推广到二维各向异性通量模型，产生对角的各向异性刚度因子和Kronecker结构算子。其次，非精确CN格式与可直接计算的对偶残差停止准则耦合。离散下界和Petrov映射常数由有限维Rayleigh商定义。在预定的有限空间自由度数集上，其有限最大值给出了稳定性与误差估计，量化了显式的图范数空间收敛阶、时间缺陷以及累积对偶残差（见表1）。

第2节介绍分数阶算子和雅可比恒等式。第3节空间离散化和第4节时间离散化与残差控制求解器给出离散化方案。第5节稳定性和第6节误差分析给出稳定性与误差分析。第7节呈现实验结果，随后是第8节。

章节片段

Riemann–Liouville算子与通量形式算子
假设0 < s_x, s_y < 1，并定义δ_x := 2 ? s_x ∈ (1,2)，δ_y := 2 ? s_y ∈ (1,2)。
引入互补阶次，因为阶数为1 ? s_ξ的通量微分的散度总阶数为2 ? s_ξ = δ_ξ。通过(2.10)，δ_ξ固定α_ξ + β_ξ = δ_ξ，并与?_ξ一起确定广义雅可比试探函数以及刚度构造与分析中使用的关联双侧分数阶雅可比算子。为连接物理扩散算子与参考域上的对应算子，我们使用

空间离散化
JSPG空间离散化在L^2(Ω; w)中建立，检验在(?, ?)_w下进行。引理2.2给出对角的各向异性算子，张量积Gauss–Jacobi求积给出矩阵条目。

时间离散化与残差控制求解器
(3.3)的CN离散化在每个时间层给出一个代数系统。其被接受的BiCGSTAB迭代定义了全离散解，代数残差保留在稳定性与误差分析中。

令Δt := T/N_t，t_n = nΔt，t_{n+1/2} = t_n + Δt/2。对任意时间依赖量□(t)，记□_n := □(t_n)，□_{n+1/2} := □(t_{n+1/2})。

设u_{n+1/2}^N := (u_{n+1}^N + u_n^N)/2，f_{n+1/2}(x, y) := f(x, y, t_{n+1/2})。

选取u_0^N ∈ U_N作为u_0的加权Petrov–Galerkin投影，即(u_0^N, v_N)_w = (u_0, v_N)_w，?v_N ∈ V_N。

设c_0为系数

稳定性
对非精确CN解(4.2)在L^2(Ω; w)中测量稳定性，使用第2.3节的内积与范数以及第4.2节的对偶残差范数。

定义试探Gram矩阵Q = (Q_{mn})_{m,n=1}^{N_tot}，其中Q_{mn} = (B_m, B_n)_w，1 ≤ m, n ≤ N_tot。

对d_N = Σ_{m=1}^{N_tot} d_m B_m，系数向量为d，(Qd)_m = Σ_{n=1}^{N_tot} Q_{mn} d_n，‖d_N‖_w^2 = d^? Q d。

回忆Petrov质量矩阵M（见(3.4)），其中M_{im} = (B_m, B*_i)_w。则对d_N = Σ_{m=1}^{N_tot} d_m B_m ∈ U_N和v_N = Σ_{i=1}^{N_tot} v_i B*_i ∈ V_N，(d_N, v_N)_w = v^? M d。

引理5.1 设M由(3.4)定义，试探/检验基为

误差分析
设u(?, t)为(1.1)的解。分析首先建立单个空间对(N_x, N_y)的估计，然后在预定的空间自由度数集上取有限最大值。带下标N的常数可能依赖于所选空间，且与Δt无关。

设P_N: L^2(Ω; w) → U_N为Petrov投影，定义为(P_N v ? v, v_N)_w = 0，?v_N ∈ V_N。

由引理5.1，此投影是良定义的。对任意序列{a_n}，定义ā_{n+1/2} := (a_{n+1} + a_n)/2。

特别地，u_{n+1/2}^N = ū_{n+1/2}^N

数值实验
计算在MATLAB R2025a中进行，运行于配备Intel(R) Core(TM) Ultra 9 285H处理器（主频2.90 GHz）和32 GB内存的64位Windows计算机上。除特别说明外，a_x = a_y = 1，CN系统通过BiCGSTAB求解，最大迭代次数为400。欧几里得残差容差由BiCGSTAB内部使用，而接受准则通过直接计算的对偶残差√(r_{n+1}^? H^{-1} r_{n+1}) / √(b_n^? H^{-1} b_n) ≤ η_res检查。检验空间Gram矩阵在张量积雅可比基下为对角阵。

结论
针对带有偏斜双侧分数阶通量的二维各向异性空间分数阶扩散方程，开发了一种JSPG方法。广义雅可比函数用作试探函数，雅可比多项式用作检验函数。双侧雅可比分数恒等式给出了一维对角的各向异性刚度矩阵。张量积构造随后导出了二维质量矩阵与刚度矩阵的Kronecker表示。空间离散化被

CRediT作者贡献声明
Chaoyue Guan：原创草稿撰写，概念构思。Yuli Sun：审阅与编辑，监督指导。Jing Niu：审阅与编辑，监督指导。

资助
本研究由HSDBSCX2025-20和中国黑龙江省自然科学基金（项目编号LH2024A018）资助。

利益冲突声明
作者声明不存在利益冲突。

Chaoyue Guan | Yuli Sun | Jing Niu