非均匀网格上分数拉普拉斯模型快速配置方法的数值分析及高效实现

《Communications in Nonlinear Science and Numerical Simulation》：Numerical analysis and efficient implementation of fast collocation methods for fractional Laplacian model on nonuniform grids

　　孔美洁|傅鸿飞 中国海洋大学数学科学学院，山东青岛266100 摘要：针对分数拉普拉斯问题，我们提出了一种基于Krylov子空间迭代求解器的快速配置方法，该方法适用于一般非均匀网格，且分数算子以奇异积分形式表示。只要指数和近似足够精确，该算法在α∈(0,1)时的一般非均匀网格上被证明具有唯一解。此外，我们还开发了一种改进方案，并证明其在α∈(0,2)时的均匀网格上同样具有唯一解。文中详细探讨了基于快速矩阵-向量乘法的快速配置方案的实现方式，包括计算复杂度和内存需求。为进一步提升计算效率，还在Krylov子空间迭代求解器中加入了带状预处理器。针对α∈(0,1)的情况，在特定的分级网格上进行了严格的最大范数误差分析，结果表明收敛阶与分级参数有关。数值实验验证了预测的收敛性，并展示了这些快速配置方案的高效性。

引言：分数拉普拉斯算子是用于模拟非局部扩散现象的典型算子。这类包含长程相互作用的非局部扩散方程近年来引起了广泛关注[1][2][3]。例如，Dipierro和Valdinoci[4]描述了一个模型，其中非局部算子与局部算子相结合，用以刻画那些在远距离跳跃和局部搜索之间交替行为的觅食动物。此外，分数拉普拉斯算子还被广泛应用于各种数学模型中，如分数Allen–Cahn方程、Cahn–Hilliard方程以及多孔介质方程[5][6]。它还是模拟复杂物理现象的强大工具，比如异常扩散[7]、准地转流[8]以及药物传输[9]。与此同时，在数值求解相关非局部模型方面也取得了不少进展，包括有限差分法[10][11][12][13]、有限元法[14][15]、配置法[16][17][18][19]、无网格方法[20][21]等等。

本文致力于高效数值求解以下一维非局部模型问题：{(?Δ)α2u=f(x)，在(a,b)区域内；u=0，在(?∞,a]∪[b,+∞)上。这里的分数拉普拉斯算子(?Δ)α2（α∈(0,2)）可以在整个空间R上以多种等价方式定义[22]，此处我们采用其超奇异积分形式：(?Δ)α2u(x)=CαP.V.∫Ru(x)?u(y)|x?y|1?αdy，其中P.V.表示柯西主值，Cα=α2α?1Γ(1+α2)π1/?Γ(1?α2)为归一化常数。需要注意的是，公式(1.2)所给出的定义体现了分数拉普拉斯算子的非局部特性，因此要计算位置x处函数u对应的算子值，就需要准确计算在整个空间R上的奇异积分。然而众所周知，任何对分数拉普拉斯算子的数值离散化都会导致稠密的刚度矩阵，从而在大型建模与仿真中带来巨大的计算成本。已有不少研究致力于针对在均匀空间网格上离散化的此类非局部模型(1.1)开发快速求解方法。例如，Duo和Zhang[10]将分数算子重新表述为中心差商的加权积分，并提出了针对二维和三维模型问题(1.1)的分数中心差分方案。他们的离散化方法能得到对称的多层托普利茨刚度矩阵。Minden和Ying[23]也提出了一种简化积分算子离散化的方法，同样可以得到托普利茨矩阵。结合基于有限差分拉普拉斯算子的预处理器，可通过预处理过的Krylov方法高效求解对应的线性系统。不过，这类收敛性分析依赖于对精确解的强正则性假设，而这一假设往往并不成立。Antil等人[18]没有直接对积分核进行离散化，而是采用了傅里叶-辛克近似来表示分数拉普拉斯算子。所有这些方法都使得基于快速傅里叶变换(FFT)的快速算法得以发展，这类算法可将每次迭代的计算成本降至O(NlogN)，同时内存需求也降至O(N)，其中N为未知数的总数。

在设计非局部问题的数值方法时，另一个重要挑战在于解的正则性较低。尽管现有大多数数值方案的收敛性分析都基于精确解的平滑性假设，但实际上，解往往存在弱边界奇异性。这类奇异行为常见于具有奇异积分核的问题中，比如时间分数方程[24]、空间分数方程[25][26]，以及非局部/区域分数拉普拉斯算子相关的问题[27][28][29][30]。根据[31]中的推论1：若f∈Cβ(Ω)，其中β=2?α+γ，且α<γ<1，那么存在常数c0=c(ω,α,β,f)，使得对于模型(1.1)，有|?? ?x?u(x)|≤c0[(x?a)(b?x)]α2??，?=>< /><α<2时，同样的正则性估计(1.3)也成立。因此，为了解决类似(1.3)中的奇异性问题并提高数值精度，人们通常会使用分级网格[29][30][33]。不过，对于非均匀空间网格，相应数值方法得到的系数矩阵会失去托普利茨结构，这样一来，前述基于fft的快速算法就不再适用。而传统的直接求解方法，如高斯消元法，无论用于何种数值离散化，都需要o(n2)的内存和o(n3)的计算复杂度。即便采用krylov子空间迭代方法，每次迭代的计算复杂度依然为o(n2)。我们最近提出了一些针对一般非均匀网格上的空间分数扩散方程的快速有限差分方法[34][35]，这些方法将每次迭代的计算成本和内存需求大幅降低至最多o(nlog?n)。不过，目前还缺乏基于正则性假设(1.3)的严格误差分析。因此，为任意非均匀网格上的非局部模型(1.1)构建并分析高效算法仍然是一项挑战。

鉴于在分数拉普拉斯方程(1.1)建模过程中存在的这些难题，本文致力于开发适用于一般非均匀网格的快速配置方法，并在低正则性假设下给出严格的误差分析。这些方法是基于所谓的指数和近似方法构建的，详见下文引理2.1，该方法能够高效地在非均匀网格上表示非局部核函数。通过充分利用指数和结构，可针对得到的稠密系数矩阵实现快速的矩阵-向量乘法运算。这样一来，无需进行任何矩阵组装操作，就能通过每次o(nlog?n)次运算得到数值解。我们严格证明了，只要选择足够精确的指数和近似误差，所提出的快速配置方案在α∈(0,1)时的非均匀网格上具有唯一解。不过，对于α∈(1,2)的情况，即便在均匀网格上，也无法从理论上保证原有快速方案的唯一解性质。为解决这一问题，我们通过改进局部积分近似，进一步提出了改进后的快速配置方案，并证明了其在均匀网格上对于更宽的α∈(0,2)范围也具有唯一解性质。此外，还引入了带状预处理器，以减少迭代次数并加快迭代速度。最重要的是，我们对原始方案和改进方案在α∈(0,1)的对称分级网格上的最大范数误差进行了严格分析，结果表明收敛阶与分级参数有关。文中还通过数值实验测试了这两种快速方案的性能，结果显示：对于α∈(0,1)，改进后的方案在均匀网格上产生的误差更小，而原始方案在非均匀网格上表现更好——尤其是当α值较小时；对于α∈[1,2)，改进方案能够达到预期的最优阶收敛性，其结果比原始方案的近似结果更为精确和稳定。

本文的其余结构如下：第2节中，我们提出了针对在一般非均匀网格上离散化的模型(1.1)的快速配置方案，并分析了其在α∈(0,1)时的唯一解性质。随后我们开发了改进后的快速方案，并研究了其在均匀网格上的唯一解性质。最后，我们详细探讨了基于快速矩阵-向量乘法的这些快速配置方案的高效实现方式，包括计算复杂度和内存需求方面的内容。第3节中，我们对所提快速方法的截断误差进行了详细分析，由此得到了在较低正则性假设下的主要误差估计。第4节则通过一系列数值实验来验证理论结果。最后一部分为结论部分。

章节要点：
快速配置方法与高效实现：本节首先针对模型问题(1.1)，提出了基于指数和近似的、与分数阶α相关的快速配置方案，然后分析了该方案在α∈(0,1)时的一般非均匀网格上的唯一解性质。随后，通过开发新的局部近似，提出了改进后的快速配置方案，并分析了其在α∈[1,2)时的均匀网格上的唯一解性质。最后，还研究了这些配置方法的高效实现方式。
误差分析：本节将误差分析范围限定在α∈(0,1)的情况。不过α∈[1,2)的情况更为复杂，未来将会进一步研究。为简化问题，我们取(a,b)=(0,2l)，给定一个偶数n，考虑以下对称分级空间网格：xi={l(2in)κ，i=0,1,?,n/2}，以及{2l?l(2?2in)κ=2l?xn?i，i=n/2+1,?,n}，其中常数κ≥1为用户自定义的分级参数。特别地，当κ=1时，即对应于均匀空间划分。
数值实验：本节通过数值实验来展示新开发的快速配置方案(2.14)和改进方案(2.29)的性能。模拟将在满足以下条件之一时终止：相对残差误差降至10??以下，或者迭代次数超过允许的最大值。指数和近似中的容差设定为?=10??。所有数值实验均使用matlab r2019b在某台笔记本电脑上完成。
结论：本文提出了一种用于数值模拟积分分数拉普拉斯问题的快速配置方案。该方案非常适合一般非均匀网格。通过仔细分析系数矩阵的数学结构，我们证明了在α∈(0,1)时，该方案在一般非均匀网格上具有唯一解。此外，我们还通过改进局部部分近似提出了改进方案，并分析了其在均匀网格上的唯一解性质。

credit作者贡献说明：
孔美洁：撰写——初稿、软件应用、方法设计、形式化分析。
傅鸿飞：撰写——审阅与编辑、指导监督、方法设计、资金获取、概念构思。

利益冲突声明：作者声明不存在任何利益冲突。 鉴于在分数拉普拉斯方程(1.1)建模过程中存在的这些难题，本文致力于开发适用于一般非均匀网格的快速配置方法，并在低正则性假设下给出严格的误差分析。这些方法是基于所谓的指数和近似方法构建的，详见下文引理2.1，该方法能够高效地在非均匀网格上表示非局部核函数。通过充分利用指数和结构，可针对得到的稠密系数矩阵实现快速的矩阵-向量乘法运算。这样一来，无需进行任何矩阵组装操作，就能通过每次o(nlog?n)次运算得到数值解。我们严格证明了，只要选择足够精确的指数和近似误差，所提出的快速配置方案在α∈(0,1)时的非均匀网格上具有唯一解。不过，对于α∈(1,2)的情况，即便在均匀网格上，也无法从理论上保证原有快速方案的唯一解性质。为解决这一问题，我们通过改进局部积分近似，进一步提出了改进后的快速配置方案，并证明了其在均匀网格上对于更宽的α∈(0,2)范围也具有唯一解性质。此外，还引入了带状预处理器，以减少迭代次数并加快迭代速度。最重要的是，我们对原始方案和改进方案在α∈(0,1)的对称分级网格上的最大范数误差进行了严格分析，结果表明收敛阶与分级参数有关。文中还通过数值实验测试了这两种快速方案的性能，结果显示：对于α∈(0,1)，改进后的方案在均匀网格上产生的误差更小，而原始方案在非均匀网格上表现更好——尤其是当α值较小时；对于α∈[1,2)，改进方案能够达到预期的最优阶收敛性，其结果比原始方案的近似结果更为精确和稳定。 本文的其余结构如下：第2节中，我们提出了针对在一般非均匀网格上离散化的模型(1.1)的快速配置方案，并分析了其在α∈(0,1)时的唯一解性质。随后我们开发了改进后的快速方案，并研究了其在均匀网格上的唯一解性质。最后，我们详细探讨了基于快速矩阵-向量乘法的这些快速配置方案的高效实现方式，包括计算复杂度和内存需求方面的内容。第3节中，我们对所提快速方法的截断误差进行了详细分析，由此得到了在较低正则性假设下的主要误差估计。第4节则通过一系列数值实验来验证理论结果。最后一部分为结论部分。 章节要点： 快速配置方法与高效实现：本节首先针对模型问题(1.1)，提出了基于指数和近似的、与分数阶α相关的快速配置方案，然后分析了该方案在α∈(0,1)时的一般非均匀网格上的唯一解性质。随后，通过开发新的局部近似，提出了改进后的快速配置方案，并分析了其在α∈[1,2)时的均匀网格上的唯一解性质。最后，还研究了这些配置方法的高效实现方式。 误差分析：本节将误差分析范围限定在α∈(0,1)的情况。不过α∈[1,2)的情况更为复杂，未来将会进一步研究。为简化问题，我们取(a,b)=(0,2L)，给定一个偶数N，考虑以下对称分级空间网格：xi={L(2iN)κ，i=0,1,?,N/2}，以及{2L?L(2?2iN)κ=2L?xN?i，i=N/2+1,?,N}，其中常数κ≥1为用户自定义的分级参数。特别地，当κ=1时，即对应于均匀空间划分。 数值实验：本节通过数值实验来展示新开发的快速配置方案(2.14)和改进方案(2.29)的性能。模拟将在满足以下条件之一时终止：相对残差误差降至10??以下，或者迭代次数超过允许的最大值。指数和近似中的容差设定为?=10??。所有数值实验均使用Matlab r2019b在某台笔记本电脑上完成。 结论：本文提出了一种用于数值模拟积分分数拉普拉斯问题的快速配置方案。该方案非常适合一般非均匀网格。通过仔细分析系数矩阵的数学结构，我们证明了在α∈(0,1)时，该方案在一般非均匀网格上具有唯一解。此外，我们还通过改进局部部分近似提出了改进方案，并分析了其在均匀网格上的唯一解性质。 credit作者贡献说明： 孔美洁：撰写——初稿、软件应用、方法设计、形式化分析。 傅鸿飞：撰写——审阅与编辑、指导监督、方法设计、资金获取、概念构思。>
鉴于在分数拉普拉斯方程(1.1)建模过程中存在的这些难题，本文致力于开发适用于一般非均匀网格的快速配置方法，并在低正则性假设下给出严格的误差分析。这些方法是基于所谓的指数和近似方法构建的，详见下文引理2.1，该方法能够高效地在非均匀网格上表示非局部核函数。通过充分利用指数和结构，可针对得到的稠密系数矩阵实现快速的矩阵-向量乘法运算。这样一来，无需进行任何矩阵组装操作，就能通过每次o(nlog?n)次运算得到数值解。我们严格证明了，只要选择足够精确的指数和近似误差，所提出的快速配置方案在α∈(0,1)时的非均匀网格上具有唯一解。不过，对于α∈(1,2)的情况，即便在均匀网格上，也无法从理论上保证原有快速方案的唯一解性质。为解决这一问题，我们通过改进局部积分近似，进一步提出了改进后的快速配置方案，并证明了其在均匀网格上对于更宽的α∈(0,2)范围也具有唯一解性质。此外，还引入了带状预处理器，以减少迭代次数并加快迭代速度。最重要的是，我们对原始方案和改进方案在α∈(0,1)的对称分级网格上的最大范数误差进行了严格分析，结果表明收敛阶与分级参数有关。文中还通过数值实验测试了这两种快速方案的性能，结果显示：对于α∈(0,1)，改进后的方案在均匀网格上产生的误差更小，而原始方案在非均匀网格上表现更好——尤其是当α值较小时；对于α∈[1,2)，改进方案能够达到预期的最优阶收敛性，其结果比原始方案的近似结果更为精确和稳定。

本文的其余结构如下：第2节中，我们提出了针对在一般非均匀网格上离散化的模型(1.1)的快速配置方案，并分析了其在α∈(0,1)时的唯一解性质。随后我们开发了改进后的快速方案，并研究了其在均匀网格上的唯一解性质。最后，我们详细探讨了基于快速矩阵-向量乘法的这些快速配置方案的高效实现方式，包括计算复杂度和内存需求方面的内容。第3节中，我们对所提快速方法的截断误差进行了详细分析，由此得到了在较低正则性假设下的主要误差估计。第4节则通过一系列数值实验来验证理论结果。最后一部分为结论部分。

章节要点：
快速配置方法与高效实现：本节首先针对模型问题(1.1)，提出了基于指数和近似的、与分数阶α相关的快速配置方案，然后分析了该方案在α∈(0,1)时的一般非均匀网格上的唯一解性质。随后，通过开发新的局部近似，提出了改进后的快速配置方案，并分析了其在α∈[1,2)时的均匀网格上的唯一解性质。最后，还研究了这些配置方法的高效实现方式。
误差分析：本节将误差分析范围限定在α∈(0,1)的情况。不过α∈[1,2)的情况更为复杂，未来将会进一步研究。为简化问题，我们取(a,b)=(0,2l)，给定一个偶数n，考虑以下对称分级空间网格：xi={l(2in)κ，i=0,1,?,n/2}，以及{2l?l(2?2in)κ=2l?xn?i，i=n/2+1,?,n}，其中常数κ≥1为用户自定义的分级参数。特别地，当κ=1时，即对应于均匀空间划分。
数值实验：本节通过数值实验来展示新开发的快速配置方案(2.14)和改进方案(2.29)的性能。模拟将在满足以下条件之一时终止：相对残差误差降至10??以下，或者迭代次数超过允许的最大值。指数和近似中的容差设定为?=10??。所有数值实验均使用matlab r2019b在某台笔记本电脑上完成。
结论：本文提出了一种用于数值模拟积分分数拉普拉斯问题的快速配置方案。该方案非常适合一般非均匀网格。通过仔细分析系数矩阵的数学结构，我们证明了在α∈(0,1)时，该方案在一般非均匀网格上具有唯一解。此外，我们还通过改进局部部分近似提出了改进方案，并分析了其在均匀网格上的唯一解性质。

credit作者贡献说明：
孔美洁：撰写——初稿、软件应用、方法设计、形式化分析。
傅鸿飞：撰写——审阅与编辑、指导监督、方法设计、资金获取、概念构思。

利益冲突声明：作者声明不存在任何利益冲突。>