扩散过程是自然界中一种常见的现象. 早在1785年, Jan Ingenhousz就描述过煤尘在酒精表面的扩散现象, 后来罗伯特布朗称这种现象为布朗运动, 傅里叶在《热的解析理论》一文中实现了对热量扩散过程的数学刻画. 之后, 经典扩散过程的数学模型成为了分析许多自然现象的有利工具. 例如, 菲克利用二阶扩散方程用来描述营养物质在细胞膜内的传输过程; 爱因斯坦从第一原理出发间接证明了分子和原子的存在; 巴舍里结合随机分析将二阶扩散方程成功运用到了股票和期权市场. 这些模型都有一个共同的假设: 随机过程都满足经典的高斯分布.

然而, 自然界中有许多复杂的扩散现象并不满足上述假设. 大量的实验和研究表明, 许多复杂系统的扩散过程不再具有高斯统计性. 例如, 生物学家观察到在某些海岛如西西里岛中每只海鸟的运动轨迹就无法用布朗运动模型来描述[1], 这样的扩散现象通常被称为反常扩散[2]. 反常扩散广泛存在于多孔介质[3], 图像分析[4], 固体表面扩散[5], 非菲克扩散湍流[6, 7, 8]等领域中. 对于这些反常扩散过程, 通常可以用分数阶扩散方程来进行数学刻画. 因此, 研究分数阶扩散方程具有重要的物理意义和现实意义: 一方面, 可以更加深刻地理解反常扩散中Lévy飞行[1, 9]等物理现象; 另一方面, 还可以为实际工程应用提供决策依据. 其中分数阶拉普拉斯算子$ -(-\Delta)^{\alpha/2}, \alpha\in(0, 2) $作为模拟反常扩散过程的原型算子[10, 11, 12], 近年来逐渐成为理论和工程领域中研究的热点.

在实际应用上, 微分方程一般都很难求出解析解, 分数阶微分方程也是一样. 即使是线性分数阶微分方程, 其解也大多也都含有特殊函数, 但计算这些函数往往非常困难; 而对大部分的非线性分数阶微分方程, 其解析解完全不可给出. 正因如此, 如何有效地对分数阶微积分方程进行数值模拟逐渐成为相关领域研究的前沿问题. 目前, 针对$ n $维分数阶拉普拉斯算子, 文献[13]讨论了Dirichlet齐次分数阶拉普拉斯方程基于积分形式的标准线性有限元方法, 得到了拟均匀网格和渐变网格下该方法的最优收敛阶, 最后给出了一些数值算例来证明其理论结果; 文献[14]在文献[13]的基础上进一步给出了二维情形的一种简洁的Matalb代码实现; 文献[15]则主要研究分数阶拉普拉斯方程的Caputo型发展方程, 结合延拓技巧将非局部的问题局部化后利用有限元方法进行数值逼近, 其时间离散采用Diethelm方法, 并且证明了时间分数阶导数与其有限部分积分的等价性, 还研究了在长时间积分下分数阶导数的某些性质; 基于$ -(-\Delta)^{\alpha/2} $的奇异积分形式, 文献[16]首先将分数阶拉普拉斯算子分为奇异部分和正则部分，然后用数值积分公式分别对其进行离散，最后构造出一种有限差分公式，类似地，文献[17]通过将分数阶拉普拉斯算子表示为弱奇异函数的加权积分, 然后利用加权的梯形公式对其进行近似, 从而提出了一种新颖的有限差分方法. 本文主要针对一维分数阶拉普拉斯算子构造一种新型的数值离散格式, 并在此基础上提出分数阶拉普拉斯方程的一种新型有限差分格式. 本文的算法思想类似于文献[16], 但分数阶拉普拉斯算子的离散方式略有不同, 本文深入分析了数值积分公式的误差展开式, 从而构造了更高阶的数值格式. 相比于文献[16]的$ O(h^{3-\alpha}) $和文献[17]的$ O(h^2) $，我们所提出的有限差分格式的精度可达到$ O(h^{4-\alpha}) $.

本文的组织如下: 第二节提出分数阶拉普拉斯算子的数值积分公式，讨论其误差展开式并给出了最优的误差估计; 第三节基于该数值积分公式构造分数阶拉普拉斯方程的一种新型有限差分格式，并讨论其误差分析; 第四节通过一些数值实验以验证有限差分格式的有效性和理论分析的正确性; 最后给出对全文的总结和对未来工作的展望。

在介绍分数阶拉普拉斯算子之前, 我们先引入一些定义和记号.

定义2.1   称$ {{\boldsymbol{\alpha}}} = (\alpha_1, \alpha_2, {\cdots}, \alpha_n) $称为多重指数, 其中$ \alpha_i $为非负整数. 称

$ D^{{\boldsymbol{\alpha}}} = \frac{\partial^{\alpha_1}}{\partial x_1^{\alpha_1}}\frac{\partial^{\alpha_2}}{\partial x_2^{\alpha_2}}{\cdots}\frac{\partial^{\alpha_n}}{\partial x_n^{\alpha_n}} $

为多重微分算子, 其次数为$ |{{\boldsymbol{\alpha}}}| = \alpha_1+\alpha_2+{\cdots}+\alpha_n $.

定义2.2   $ \mathbb R^n $中的速降函数空间$ {\mathcal{S}} $是满足以下条件的函数的集合

$ {\mathcal{S}}({\mathbb{R}}^n) = \left\{f\in C^\infty({\mathbb{R}}^n) \; |\; \|f\|_{{{\boldsymbol{\alpha}}}, {{\boldsymbol{\beta}}}} < \infty, \; \forall{{\boldsymbol{\alpha}}}, {{\boldsymbol{\beta}}} \right\}, $

其中$ C^\infty({\mathbb{R}}^n) $是光滑函数空间, $ \|f\|_{{{\boldsymbol{\alpha}}}, {{\boldsymbol{\beta}}}} = \sup_{{\boldsymbol{x}}\in{\mathbb{R}}^n} |{\boldsymbol{x}}^{{\boldsymbol{\alpha}}} D^{{\boldsymbol{\beta}}} f({\boldsymbol{x}})|, $且$ {\boldsymbol{x}}^{{\boldsymbol{\alpha}}} = x_1^{\alpha_1}x_2^{\alpha_2}{\cdots} x_n^{\alpha_n}. $

定义2.3   设函数$ u: \mathbb R^n \to \mathbb R $为光滑速降函数, 即$ u\in{\mathcal{S}}({\mathbb{R}}^n) $, 则分数阶拉普拉斯算子$ (-\Delta)^{\alpha/2} $定义为

$ \begin{equation} (-\Delta)^{\alpha/2}u({\boldsymbol{x}}) = c_{n, \alpha}\int_{\mathbb{R}^n}\frac{u({\boldsymbol{x}})-u({\boldsymbol{y}})}{{\left\|{\boldsymbol{x}}-{\boldsymbol{y}}\right\|}^{n+\alpha}}\, d{\boldsymbol{y}}, \quad \alpha \in (0, 2), \end{equation} $

(2.1)

其中系数$ c_{n, \alpha} $由下式给出

$ c_{n, \alpha} = \frac{\alpha2^{\alpha-1}\Gamma(\frac{\alpha+n}{2})} {\pi^\frac{n}{2}\Gamma(\frac{2-\alpha}{2})}, $

这里$ \Gamma(x) $是Gamma函数, 即

$ \Gamma(x) = \int_0^\infty t^{x-1} e^{-t}\, dt, \quad x>0. $

本文主要关注一维分数阶拉普拉斯算子, 即

$ \begin{equation} (-\Delta)^{\alpha/2}u(x) = c_{1, \alpha}\int_{\mathbb{R}}\frac{u(x)-u(y)}{{\left\lvert x-y\right\rvert}^{1+\alpha}}dy. \end{equation} $

(2.2)

的数值离散格式. 引入变量替换$ y = x+z $, 式(2.2) 可改写为

$ \begin{equation} \begin{aligned} (-\Delta)^{\alpha/2}u(x) & = -c_{1, \alpha}\int_{0}^{\infty}\frac{u(x-z)-2u(x)+u(x+z)}{z^{1+\alpha}}dz. \end{aligned} \end{equation} $

(2.3)

定义

$ \begin{equation} v(x, z) = u(x-z)-2u(x)+u(x+z), \end{equation} $

(2.4)

则式(2.3) 可简化为

$ \begin{equation} (-\Delta)^{\alpha/2}u(x) = -c_{1, \alpha}\int_{0}^{\infty}\frac{v(x, z)}{z^{1+\alpha}}dz. \end{equation} $

(2.5)

以下将基于式(2.5) 构造分数阶拉普拉斯算子(2.2) 的一种新型数值求积公式. 为此, 我们对$ {\mathbb{R}}^+ $进行等距的网格剖分, 得到等距节点$ x_j = jh, \; j = 0, 1, \cdots, $其中$ h $为步长. 另外, 我们还将式(2.5) 分成两部分进行处理, 即

$ \begin{equation} \mathcal{Q}^S u(x) = -c_{1, \alpha}\int_{0}^{h}\frac{v(x, z)}{z^{1+\alpha}}dz \end{equation} $

(2.6)

和

$ \begin{equation} \mathcal{Q}^Ru(x) = -c_{1, \alpha}\int_{h}^{\infty}\frac{v(x, z)}{z^{1+\alpha}}dz. \end{equation} $

(2.7)

由于积分$ \mathcal{Q}^Su(x) $具有奇异性, 故称之为奇异部分(Singular Part); 而积分$ \mathcal{Q}^Ru(x) $不包含奇异性, 故称之为正则部分(Regular Part).

设$ u\in C^4({\mathbb{R}}) $, 利用泰勒展开式可得

$ u(x+z) = u(x)+u'(x)z+\frac{u^{\prime\prime}(x)}{2!}z^2+\frac{u^{\prime\prime\prime}(x)}{3!}z^3+\frac{u^{(4)}(\xi_1)}{4!}z^4, \quad \xi_1\in(0, h), $

和

$ u(x-z) = u(x)-u'(x)z+\frac{u^{\prime\prime}(x)}{2!}z^2-\frac{u^{\prime\prime\prime}(x)}{3!}z^3+\frac{u^{(4)}(\xi_2)}{4!}z^4, \quad \xi_2\in(-h, 0). $

将以上两式代入式(2.4) 得

$ \begin{equation} v(x, z) = u^{\prime\prime}(x)z^2+\frac{u^{(4)}(\xi_1)+u^{(4)}(\xi_2)}{24}z^4, \quad \xi\in (0, h). \end{equation} $

(2.8)

再将上式代入(2.6) 可得

$ \begin{equation} \begin{aligned} \mathcal{Q}^Su(x) & = -c_{1, \alpha}\left(u^{\prime\prime}(x)\int_{0}^{h}z^{1-\alpha}dz+\frac{1}{24} \int_{0}^{h}\left(u^{(4)}(\xi_1)+u^{(4)}(\xi_2)\right)z^{3-\alpha}dz\right)\\ & = -c_{1, \alpha}\sigma_0u^{\prime\prime}(x)+O(h^{4-\alpha}), \end{aligned} \end{equation} $

(2.9)

其中

$ \begin{equation} \sigma_0 = \frac{h^{2-\alpha}}{2-\alpha}. \end{equation} $

(2.10)

利用二阶中心差商$ \delta_{xx}u(x) = \frac{1}{h^2}[u(x+h)-2u(x)+u(x-h)] $近似(2.9) 中的$ u^{\prime\prime}(x) $, 即得$ \mathcal{Q}^Su(x) $的近似

$ \begin{equation} \mathcal{Q}_{h}^Su(x) = -c_{1, \alpha}\sigma_0\delta_{xx}u(x), \end{equation} $

(2.11)

其误差估计由如下定理给出.

定理2.1    设$ u\in C^4(\mathbb R) $, 则对于奇异部分(2.6) 的近似公式(2.11), 以下误差估计式成立:

$ \begin{equation} \mathcal{Q}^Su(x)-\mathcal{Q}_{h}^Su(x) = O(h^{4-\alpha}). \end{equation} $

(2.12)

证   该结论可由二阶中心差商的误差估计式

$ \begin{equation} u^{\prime\prime}(x) = \delta_{xx}u(x)+O(h^2), \end{equation} $

(2.13)

及(2.9) 和(2.11) 直接推出.

在每个子区间$ [x_j, x_{j+1}], \; j = 1, 2, \cdots $中, 用$ v(x, z) $的分段线性插值

$ \begin{equation} \mathcal L_{h}v(x, z) = v(x, x_j)\frac{x_{j+1}-z}{h}+v(x, x_{j+1})\frac{z-x_j}{h} \end{equation} $

(2.14)

代替正则部分(2.7) 中的$ v(x, z) $可得

$ \begin{equation} \begin{aligned} \widetilde{\mathcal{Q}}_{h}^R u(x) & = -c_{1, \alpha}h^{-1}\sum\limits_{j = 1}^{\infty}\left[v(x, x_j)\int_{x_j}^{x_{j+1}}\frac{x_{j+1}-z}{z^{1+\alpha}}dz+v(x, x_{j+1})\int_{x_j}^{x_{j+1}}\frac{z-x_j}{z^{1+\alpha}} dz\right]\\ &: = -c_{1, \alpha}\sum\limits_{j = 1}^{\infty}\widetilde{w}_{j}\left[u(x-x_j)-2u(x)+u(x+x_j) \right], \end{aligned} \end{equation} $

(2.15)

其中$ \widetilde{w}_{j} $为求积系数, 其具体形式为

$ \begin{equation} \begin{aligned} \widetilde{w}_{j}& = h^{-1}\left\{ \begin{aligned} &\int_{x_1}^{x_2}\frac{x_2-z}{z^{1+\alpha}}dz, \quad &j = 1, \\ &\int_{x_{j-1}}^{x_j}\frac{z-x_{j-1}}{z^{1+\alpha}}dz+\int_{x_j}^{x_{j+1}}\frac{x_{j+1}-z}{z^{1+\alpha}}dz, \quad &j\geq 2.\\ \end{aligned} \right.\\ \end{aligned} \end{equation} $

(2.16)

定义函数空间$ C_0^k({\mathbb{R}}) = \left\{ u \in C^k({\mathbb{R}}) \; |\; u^{(i)} \in C_0({\mathbb{R}}), i = 0, {\cdots}, k\right\}, $其中$ C_0({\mathbb{R}}) = \{u: {\mathbb{R}} \to {\mathbb{R}}, u\mbox{连续且具有紧支集}\}. $

定理2.2    设$ u\in C^4_0(\mathbb R) $, 则对于正则部分(2.7) 的数值求积公式(2.15), 有如下误差展开式:

$ \begin{equation} \mathcal{Q}^R u(x)-\widetilde{\mathcal{Q}}_{h}^R u(x) = -c_{1, \alpha}\sum\limits_{j = -\infty, j\neq0}^{\infty} \sigma_j u^{\prime\prime}(x+x_j)+O(h^{4-\alpha}), \end{equation} $

(2.17)

其中

$ \begin{equation} \begin{aligned} \sigma_j& = \frac{1}{2}\int_{x_j}^{x_{j+1}}\frac{(z-x_j)(z-x_{j+1})}{z^{1+\alpha}}dz, \quad j \ge 1, \end{aligned} \end{equation} $

(2.18)

证   定义误差函数$ e_j(z) = v(x, z)- \mathcal L_{h}v(x, z), \; \; z\in [x_j, x_{j+1}], \; \; j\geq 1. $因为$ u\in C^{4}(\mathbb R) $, 由$ v(x, z) $的四阶泰勒展开式

$ \begin{equation} \left\{ \begin{aligned} &v(x, z)-v(x, x_j) = \sum\limits_{k = 1}^{3}\frac{v^{(k)}(x, x_j)(z-x_j)^k}{k!}+\frac{v^{(4)}(x, \xi)(z-x_j)^4}{4!}, \quad &\xi \in (x_j, z), \\ &v(x, x_{j+1})-v(x, x_j) = \sum\limits_{k = 1}^{3}\frac{v^{(k)}(x, x_j)h^k}{k!}+\frac{v^{(4)}(x, \eta)h^4}{4!}, \quad &\eta\in(x_j, x_{j+1}), \end{aligned} \right. \end{equation} $

(2.19)

可得

$ \begin{equation} \begin{aligned} e_j(z) = &v(x, z)-v(x, x_j)-\frac{v(x, x_{j+1})-v(x, x_j)}{h}(z-x_j)\\ = &\frac{v^{\prime\prime}(x, x_j)}{2!}(z-x_j)(z-x_{j+1})+\frac{v^{\prime\prime\prime}(x, x_j)}{3！}(z-x_{j-1})(z-x_j)(z-x_{j+1})+\mathcal R_j, \end{aligned} \end{equation} $

(2.20)

其中

$ \left|\mathcal R_j\right| = \left|\frac{v^{(4)}(x, \xi)(z-x_j)^4}{4!}-\frac{v^{(4)}(x, \eta)(z-x_j)h^3}{4!}\right|\le Ch^4. $

于是

$ \begin{equation} \begin{aligned} \mathcal{Q}^R u(x)-\widetilde{\mathcal{Q}}_{h}^Ru(x) & = \sum\limits_{j = 1}^{\infty}\frac{v^{\prime\prime}(x, x_j)}{2！}\int_{x_j}^{x_{j+1}}\frac{(z-x_j)(z-x_{j+1})}{z^{1+\alpha }}dz\\ &\quad +\sum\limits_{j = 1}^{\infty}\frac{v^{\prime\prime\prime}(x, x_j)}{3!}\int_{x_j}^{x_{j+1}}\frac{(z-x_{j-1})(z-x_j)(z-x_{j+1})}{z^{1+\alpha}}dz+\mathcal R, \end{aligned} \end{equation} $

(2.21)

其中

$ \begin{equation} \left|\mathcal R\right| = \left|\sum\limits_{j = 1}^{\infty}\int_{x_j}^{x_{j+1}}\mathcal R_jz^{-1-\alpha}dz\right| \leq Ch^4\sum\limits_{j = 1}^{\infty}\int_{x_j}^{x_{j+1}}z^{-1-\alpha}dz = C h^{4-\alpha}. \end{equation} $

(2.22)

对$ v(x, z) $在节点$ x_j $处关于$ z $求三阶导数, 可得

$ v^{\prime\prime\prime}(x, x_j) = u^{\prime\prime\prime}(x+x_j)-u^{\prime\prime\prime}(x-x_j) = 2u^{(4)}(\xi)x_j, \; \; \; j\ge 1, $

从而有

$ \begin{equation} \begin{aligned} &\left\lvert \sum\limits_{j = 1}^{\infty}\frac{v^{\prime\prime\prime}(x, x_j)}{3!}\int_{x_j}^{x_{j+1}}\frac{(z-x_{j-1})(z-x_j)(z-x_{j+1})}{z^{1+\alpha}}dz\right \rvert\\ \leq& Ch^3\sum\limits_{j = 1}^{\infty}\int_{x_j}^{x_{j+1}}\frac{x_j}{z^{1+\alpha}}dz\\ \leq& Ch^3\sum\limits_{j = 1}^{\infty}\int_{x_j}^{x_{j+1}}\frac{z-x_j}{z^{1+\alpha}}dz+Ch^3\sum\limits_{j = 1}^{\infty} \int_{x_j}^{x_{j+1}}\frac{1}{z^\alpha}dz\leq Ch^{4-\alpha}. \end{aligned} \end{equation} $

(2.23)

联立(2.21)– (2.23)易得(2.19), 定理得证.

更进一步, 用中心差商

$ \delta_{xx}u(x+x_j) = \frac{u(x+x_{j-1})-2u(x+x_j)+u(x+x_{j+1})}{h^2} $

近似$ (2.17) $中的$ u^{\prime\prime}(x+x_j) $, 即得正则部分(2.7)的一种改进的数值积分公式:

$ \begin{equation} \mathcal{Q}_{h}^R u(x) = \widetilde{\mathcal{Q}}_{h}^Ru(x) -c_{1, \alpha} \sum\limits_{j = -\infty, j\neq0}^{\infty} \sigma_j \delta_{xx}u(x+x_j), \end{equation} $

(2.24)

其误差估计可通过如下定理给出.

定理2.3    设$ u\in C^4_0(\mathbb R) $, 则对于正则部分(2.7) 的数值求积公式(2.24), 以下误差估计式成立:

$ \begin{equation} \mathcal{Q}^R u(x)-\mathcal{Q}_{h}^R u(x) = O(h^{4-\alpha}). \end{equation} $

(2.25)

证   由定理2.2和(2.24) 可知

$ \begin{equation} \mathcal{Q}^R u(x)-\mathcal{Q}_{h}^R u(x) = -c_{1, \alpha}\sum\limits_{j = -\infty, \; j\ne 0}^\infty \sigma_j \left[u^{\prime\prime}(x+x_j)-\delta_{xx} u(x+x_j)\right]+O(h^{4-\alpha}). \end{equation} $

(2.26)

又由$ \sigma_j $的定义可知$ \sigma_j = \sigma_{-j} $, 从而有

$ -\sum\limits_{j = -\infty, j\ne 0}^\infty \sigma_j = -2\sum\limits_{j = 1}^\infty \sigma_j \le Ch^2\int_h^\infty z^{-1-\alpha}\, dz = Ch^{2-\alpha}. $

再由中心差商格式的误差估计$ u^{\prime\prime}(x+x_j) = \delta_{xx} u(x+x_j)+O(h^2) $可得

$ \begin{aligned} \left|\mathcal{Q}^R u(x)-\mathcal{Q}_{h}^R u(x)\right|&\le -Ch^2\sum\limits_{j = -\infty, \; j\ne 0}^\infty \sigma_j +O(h^{4-\alpha}) \le Ch^{4-\alpha}. \end{aligned} $

证毕.

联立(2.11) 和(2.24) 可得一维分数阶拉普拉斯算子(2.2) 的数值求积公式

$ \begin{equation} (-\Delta_{h})^{\alpha/2}u(x) = \mathcal{Q}_{h}^S u(x) + \mathcal{Q}_{h}^R u(x) = -c_{1, \alpha}\sum\limits_{j = 1}^\infty\omega_{j} \left[u(x-x_j)-2u(x)+u(x+x_j)\right], \end{equation} $

(2.27)

其中求积系数$ w_{j} $为

$ \begin{equation} \omega_{j} = \widetilde{ \omega}_{j} + h^{-2}\rho_j, \end{equation} $

(2.28)

$ \begin{equation} \rho_j = \sigma_{j-1}-2\sigma_j+\sigma_{j+1}, \end{equation} $

(2.29)

而$ \widetilde{ \omega}_{j} $由(2.16)给出. 其误差估计由如下定理给出, 其证明可结合定理2.1和2.3直接得到.

定理2.4    设$ u\in C^4_0(\mathbb R) $, 则对于分数阶拉普拉斯算子(2.2) 的数值求积公式(2.27), 以下误差估计式成立:

$ \begin{equation} (-\Delta)^{\alpha/2}u(x)-(-\Delta_{h})^{\alpha/2}u(x) = O(h^{4-\alpha}). \end{equation} $

(2.30)

本节将基于分数阶拉普拉斯算子(2.2) 的数值求积公式(2.27), 来构造如下分数阶拉普拉斯方程Dirichlet边值问题的有限差分格式:

$ \begin{equation} \left\{ \begin{array}{rl} (-\Delta)^{\alpha/2}u(x) = f(x), & x\in \Omega = (-1, 1), \\ u(x) = g(x), & x \in \mathbb R\backslash \Omega, \end{array} \right. \end{equation} $

(3.1)

其中$ f $和$ g $已知.

首先我们将求解区域$ (-1, 1) $等分为$ 2J $个子区间, 其步长为$ h = 1/J $, 则网格节点可表示为$ x_k = kh, k \in \mathbb Z. $记$ \Omega_h = \{-J+1, \cdots, 0, \cdots, J-1\}, \quad \Omega_h^c = \mathbb Z \backslash \Omega_h, $在节点$ x_i, \; i \in \Omega_h $处, 用求积公式(2.27) 代替(3.1) 中的分数阶拉普拉斯算子, 即得有限差分格式

$ \begin{equation} \left\{ \begin{array}{rl} (-\Delta_h)^{\alpha/2}u_i = f(x_i), & i \in \Omega_h, \\ u_i = g(x_i), & i\in \Omega_h^c. \end{array} \right. \end{equation} $

(3.2)

有限差分格式(3.2) 可表示为矩阵形式

$ \begin{equation} {\boldsymbol{A}} {\boldsymbol{u}} = {\boldsymbol{f}}, \end{equation} $

(3.3)

其中

$ {\boldsymbol{u}} = \begin{bmatrix} u_0&u_1&\cdots&u_{2J-2} \end{bmatrix}^T, \quad {\boldsymbol{f}} = \begin{bmatrix} f_0&f_1&\cdots&f_{2J-2} \end{bmatrix}^T, $

而$ f_i $由右端项$ f(x_i) $及边界条件$ g(x_i) $确定. 显然, 系数矩阵$ {\boldsymbol{A}} $为对称的Toeplitz矩阵, 即

$ \begin{equation} {\boldsymbol{A}} = \begin{bmatrix} a_0&a_1&\cdots&a_{2J-2}\\ a_1&a_0&\cdots&a_{2J-3}\\ \vdots&\vdots&\ddots&\vdots\\ a_{2J-2}&a_{2J-3}&\cdots&a_{0} \end{bmatrix}, \end{equation} $

(3.4)

其中

$ \begin{equation} a_0 = -2\sum\limits_{j = 1}^\infty a_j, \quad a_j = -c_{1, \alpha}w_j, \quad j \ge 1, \end{equation} $

(3.5)

而$ w_j(j\ge 1) $由式(2.28) 定义.

以下将讨论差分格式(3.2)的误差估计. 在此之前, 我们先介绍几个重要的定义和引理.

定义3.1   称$ {\boldsymbol{A}}\in\mathbb R^{n\times n} $为L矩阵, 若$ a_{ii}>0 $而$ a_{ij}\le 0(i\ne j) $; 称$ {\boldsymbol{A}}\in\mathbb R^{n\times n} $为M矩阵, 若$ {\boldsymbol{A}} $为L矩阵, 且其所有特征值的实部皆为正.

引理3.1   (Gershgorin圆盘定理) 设$ {\boldsymbol{A}} = (a_{ij}) $为$ n $阶方阵, $ G_i({\boldsymbol{A}}) $是复平面上以$ a_{ii} $为圆心, 半径为$ \sum_{j = 1, j\neq i}^{n}|a_{ij}| $的圆盘, 即

$ G_i({\boldsymbol{A}}) = \{ z\in \mathbb{C}\big| |z-a_{ii}|\leq \sum\limits_{j = 1, j\neq i}^{n}|a_{ij}| \}, \; \; i = 1, 2, ....n, $

则$ {\boldsymbol{A}} $的任一特征值$ \lambda $皆满足$ \lambda\in G = \bigcup_{i = 1}^n G_i(A). $

引理3.2   对于定义于式(2.28) 中的$ w_j $, 有如下性质

$ \begin{equation} w_j \ge 0, \quad j \ge 1. \end{equation} $

(3.6)

证   由(2.28) 可知$ w_j = \widetilde{w}_j + h^{-2}\sigma_{j-1} - h^{-2}\sigma_{j} + h^{-2} \left(\sigma_{j+1}-\sigma_{j}\right). $由(2.16) 和(2.18) 可知,

(1) 当$ j = 1 $时,

$ \begin{equation} w_1 = \widetilde{w}_1+ h^{-2}\sigma_{0} = h^{-1}\int_{x_{1}}^{x_2}\frac{x_{2}-z}{z^{1+\alpha}}dz + \frac{h^{-\alpha}}{2-\alpha} > 0; \end{equation} $

(3.7)

(2) 当$ j\ge 2 $时,

$ \begin{equation} \widetilde{w}_j + h^{-2}\sigma_{j-1} = \frac{h^{-2}}2\int_{x_{j-1}}^{x_{j}}\frac{(z-x_{j-1})(z-x_{j-2})}{z^{1+\alpha}}dz + h^{-1}\int_{x_j}^{x_{j+1}}\frac{x_{j+1}-z}{z^{1+\alpha}}dz > 0. \end{equation} $

(3.8)

由$ \sigma_j $的定义(2.18) 可知对任意的$ j\ge 1 $皆有$ \sigma_j \le 0 $. 再由

$ \begin{equation} \begin{aligned} \sigma_{j+1}-\sigma_{j} = \frac{h^{2-\alpha}}{2}\;\int_{0}^{1}\left[\frac{1}{(\tau+j+1)^{1+\alpha}}-\frac{1}{(\tau+j)^{1+\alpha}}\right]\tau(\tau-1)\, d\tau > 0 \end{aligned} \end{equation} $

(3.9)

可知(3.6) 成立.

定理3.1    对于有限差分格式(3.2) 的系数矩阵$ {\boldsymbol{A}} $, 它为对称正定的M矩阵.

证   由引理3.2和(3.5) 可知$ {\boldsymbol{A}} $是$ L $矩阵. 由M矩阵的定义可知, 欲证$ {\boldsymbol{A}} $为M矩阵, 只需说明其特征值均为正. 事实上, 由(3.4) 的定义以及(3.5) 容易推导出$ {\boldsymbol{A}} $为严格对角占优矩阵, 再由引理3.1可知, $ {\boldsymbol{A}} $的任一特征值皆为正.

引理3.3   (极值原理) 若$ (-\Delta_h)^{{\alpha}/{2}}u_i\leq 0, i \in \Omega_h $, 则对任意的$ i\in \Omega_h $, 有$ \max\limits_{i\in \Omega_h} u_i\leq \max\limits_{i\in \Omega_h^\bf{C}}u_i. $类似地, 若$ {(-\Delta_h)}^{\alpha/2}u\geq 0 $, 则对任意的$ i\in \Omega_h $, 有$ \min\limits_{i\in\Omega_h}u_i\geq \min\limits_{i\in \Omega^C_h}u_i. $

证   我们只证明第一种情形, 第二种情形的证明过程类似. 利用反证法, 假设$ u_{i_0} = \max\limits_{i\in \Omega_h}u_i>\max\limits_{i\in\Omega^C_h} u_i $这说明$ u_{i_0} $为$ u $在所有节点上的最大值, 但是由$ w_j(j\ge 1) $的非负性可知

$ (-\Delta_h)^{\alpha/2}u_{i_0} = -c_{1, \alpha}\sum\limits_{j = 1}^\infty w_j (u_{i_0-j}-2u_{i_i}+u_{i_0-j})>0. $

这与已知条件矛盾, 故假设不成立, 也就是说$ u $在$ x_i, i\in \Omega_h $上无法取得最大值, 定理得证.

为证明有限差分方程的收敛性, 我们引入一个相容函数$ v(x) $:

$ \begin{equation} v(x) = \left\{ \begin{array}{cl} 4-x^2, & |x|<1, \\ 0, & \mbox{otherwise}. \end{array} \right. \end{equation} $

(3.10)

引理3.4   对于函数$ v(x) $, 它满足

$ \begin{equation} \left\{ \begin{array}{rl} (-\Delta_h)^{\alpha/2} v_i \geq 1, & i\in \Omega_h, \\ v_i = 0, & i \in \Omega^C_h, \end{array} \right. \end{equation} $

(3.11)

其中$ v_i = v(x_i) $.

证   定义$ \delta_ju_i = u_{i+j} - 2u_i + u_{i-j}. $首先, 证明$ -\delta_jv_i\geq \min\left\{2, 2(jh)^2\right\}, \; \forall i\in \Omega_h, \; \forall j\geq 1. $事实上,

(1) 当$ i+j $和$ i-j $都属于$ \Omega_h $, 有$ -\delta_jv_i = h^2((i+j)^2-2i^2+(i-j)^2) = 2(jh)^2; $

(2) 当$ i+j $和$ i-j $都不属于$ \Omega_h $, 有$ -\delta_jv_i = 2v_i = 8-2(ih)^2\geq 6; $

(3) 当$ i\pm j $中的一个属于$ \Omega_h $, $ -\delta_jv_i = -v((i\pm j)h)+2u(x_i)\geq 2. $

因此,

$ \begin{equation} \begin{aligned} (-\Delta_h)^{\alpha/2}v_i\geq c_{1, \alpha}\sum\limits_{j = 1}^{\infty}w_j\min\left\{2, 2(jh)^2\right\} & = 2c_{1, \alpha}\left(\sum\limits_{j\ge J}w_j+\sum\limits_{j< J}w_jx_j^2\right). \end{aligned} \end{equation} $

(3.12)

接下来, 我们将证明

$ \begin{equation} \sum\limits_{j \ge J}w_j \ge \frac{1}{\alpha}, \; \sum\limits_{j<J}w_jx_j^2 \ge \frac{1}{2-\alpha}. \end{equation} $

(3.13)

一方面, 由$ \widetilde{\omega}_{j} $的定义可知

$ \sum\limits_{j\ge J} \widetilde{\omega}_{j} = h^{-1}\sum\limits_{j\ge J} \left[\int_{x_{j-1}}^{x_j} \frac{z-x_{j-1}}{z^{1+\alpha}}\, dz +\int_{x_j}^{x_{j+1}} \frac{x_{j+1}-z}{z^{1+\alpha}}\, dz\right] = \int_1^\infty z^{-1-\alpha}\, dz + h^{-1}\int_{x_{J-1}}^{x_{J}} \frac{x_{J}-z}{z^{1+\alpha}}\, dz $

从而有

$ \begin{equation*} \begin{aligned} \sum\limits_{j\ge J}w_j & = \sum\limits_{j\ge J} \left[\widetilde{\omega}_{j}+h^{-2}(\sigma_{j-1}-2\sigma_j+\sigma_{j+1})\right]\\ & = \int_{1}^{\infty}z^{-1-\alpha}dz + h^{-1}\int_{x_{J-1}}^{x_{J}} \frac{x_{J}-z}{z^{1+\alpha}}\, dz + h^{-2}(\sigma_{J-1} - \sigma_J) \\ & = \frac{1}{\alpha} +h^{-2}\int_{x_{J-1}}^{x_{J}}\frac{(x_{J}-z)^2}{z^{1+\alpha}}\, dz - h^{-2}(\sigma_{J-1} + \sigma_J) \ge \frac{1}{\alpha}, \end{aligned} \end{equation*} $

这里用到了$ \sigma_j $的非正性. 另一方面, 再由$ \widetilde{\omega}_{j} $的定义可知

$ \begin{equation*} \begin{aligned} \sum\limits_{0\le j < J}w_jx_j^2 = & \sum\limits_{1\le j< J}\left[\widetilde{\omega}_{j}+ \frac{\sigma_{j-1}-2\sigma_j+\sigma_{j+1}}{h^2}\right]x_j^2 \\ = &\sum\limits_{1\le j<J} \int_{x_j}^{x_{j+1}} \frac{jx_j(x_{j+1}-z)+(j+1)x_{j+1}(z-x_j)+(z-x_{j-1})(z-x_j)}{z^{1+\alpha}}dz\\ &+\sigma_0h^{-2}x_1^2 +x_{J-1}^2\left[h^{-2}\sigma_J-h^{-1}\int_{x_{J-1}}^{x_J}\frac{(z-x_{J-1})}{z^{1+\alpha}}dz\right] -x_J^2h^{-2}\sigma_{J-1} \\ = &\int_0^1z^{1-\alpha}dz +\frac{1}{2}\int_{x_{J-1}}^{x_J}\frac{(z-x_{J-1})(x_{J+1}-z)}{z^{1+\alpha}}dz -J^2 \sigma_{J-1} \ge \frac{1}{2-\alpha}, \end{aligned} \end{equation*} $

这里同样用到了$ \sigma_j $的非正性. 最后, 由(3.12) 和(3.13) 可知,

$ (-\Delta_h)^{\alpha/2}v_i\geq 2c_{1, \alpha}(\frac{1}{2-\alpha}+\frac{1}{\alpha}) = \frac{2^{\alpha}\Gamma(\frac{\alpha+1}{2})} {\Gamma(\frac{1}{2})\Gamma(2-\frac{\alpha}{2})}. $

再由$ 2^{\alpha}>1, \Gamma(\frac{\alpha+1}{2})>\pi^\frac{1}{2} = \Gamma(\frac{1}{2}) $, 以及$ \Gamma(2-\frac{\alpha}{2})<\Gamma(2) = 1 $可知

$ (-\Delta_h)^{\alpha/2}v_i\ge1. $

定理得证.

定理3.2    对于有限差分方程(3.2), 有误差估计$ \|{\boldsymbol{u}}-{\boldsymbol{u}}_h\|_\infty \le Ch^{4-\alpha}, $其中$ {\boldsymbol{u}} = \left[u(x_{-J+1}), \cdots, u(x_{J-1})\right]^T $和$ {\boldsymbol{u}}_h = \left[u_{-J+1}, \cdots, u_{J-1}\right]^T $分别表示(3.1)的真解和(3.2)的近似解.

证   由于$ {\boldsymbol{A}} $是M矩阵, 故$ (-\Delta_h)^{\alpha/2} $可逆, 不妨记其逆为$ (-\Delta_h)^{-\alpha/2} $. 将(3.11)中的第一式改写为向量形式$ (-\Delta_h)^{\alpha/2} {\boldsymbol{v}}_h \ge {\boldsymbol{1}}, \quad {\boldsymbol{1}} = \left( 1, 1, \cdots, 1 \right)^T, $由极值原理可知

$ {\boldsymbol{v}}_h \ge (-\Delta_h)^{-\alpha/2}{\boldsymbol{1}}, $

于是,

$ \left\|\left(-\Delta_h\right)^{-\alpha/2}\right\|_{\infty}\leq \|{\boldsymbol{v}}\|_{\infty} = 4. $

由(3.1) 和(3.2) 可知

$ (-\Delta)^{\alpha/2}{\boldsymbol{u}} = (-\Delta_h)^{\alpha/2}{\boldsymbol{u}}_h, $

于是

$ (-\Delta_h)^{\alpha/2}({\boldsymbol{u}} - {\boldsymbol{u}}_h) = (-\Delta_h)^{\alpha/2}{\boldsymbol{u}}-(-\Delta)^{\alpha/2}{\boldsymbol{u}}, $

亦即

$ {\boldsymbol{u}} - {\boldsymbol{u}}_h = (-\Delta_h)^{\alpha/2}\left[(-\Delta_h)^{\alpha/2}{\boldsymbol{u}}-(-\Delta)^{\alpha/2}{\boldsymbol{u}}\right]. $

利用求积公式的误差估计(2.30) 可得

$ \|{\boldsymbol{u}}-{\boldsymbol{u}}_h\|_{\infty} \le \left\|\left(-\Delta_h\right)^{-\alpha/2}\right\|_{\infty}\cdot{\|(-\Delta_h)^{\alpha/2}{\boldsymbol{u}}-(-\Delta_h)^{\alpha/2}{\boldsymbol{u}}_h\|}_{\infty} \leq C h^{4-\alpha}. $

定理得证.

注意到分数阶拉普拉斯算子的积分区域为整个实数域, 在实际计算中需要对其进行截断, 为此我们将其改写为如下形式

$ \begin{aligned} (-\Delta)^{\alpha/2}u(x) = &-c_{1, \alpha}\int_{0}^{\infty}\left(u(x-z)-2u(x)+u(x+z)\right)z^{-1-\alpha}dz\\ = &\underbrace{-c_{1, \alpha}\int_{0}^{L}(u(x-z)-2u(x)+u(x+z))z^{-1-\alpha}dz}_{\mathcal Q_1}+\underbrace{2c_{1, \alpha}\int_{L}^{\infty}u(x)z^{-1-\alpha}dz}_{\mathcal Q_2}\\ &-\underbrace{c_{1, \alpha}\int_{L}^{\infty}(u(x-z)+u(x+z))z^{-1-\alpha}dz}_{\mathcal Q_3}\\ \end{aligned} $

其中$ [-L, L] $为截断区间, 而$ L $通常选取为一个比较大的正数.

以下将讨论$ \mathcal Q_1, \mathcal Q_2, \mathcal Q_3 $对分数阶拉普拉斯算子近似计算的影响.

(1) 可直接用数值积分公式近似$ \mathcal Q_1 $时, 此时其截断误差$ O(h^{4-\alpha}) $, 它与$ L $的大小无关.

(2) $ \mathcal Q_2 $可直接计算出来, 即

$ \mathcal Q_2 = 2c_{1, \alpha}\int_{L}^{\infty} u(x)z^{-1-\alpha}dz = 2c_{1, \alpha}u(x)\int_{L}^{\infty}z^{-1-\alpha}dz = \frac{2c_{1, \alpha}}{\alpha L^{\alpha}}u(x). $

(3) $ \mathcal Q_3 $需要具体问题具体分析. 对于分数阶拉普拉斯方程(3.1), 若边界条件是齐次的, 即$ g = 0 $, 则$ \mathcal Q_3 $可以忽略；而对于有的分数阶拉普拉斯方程, 如分数阶Burgers方程[18], 分数阶分数多孔介质方程[19]等, 它们的解具有代数衰减性, 此时即使$ L $取很大, $ \mathcal Q_3 $也不能忽略不计. 事实上, 若$ u(x)\rightarrow 0 $具有代数衰减性, 即$ u(z)\sim |z|^{-\beta} $, 更具体地,

$ \begin{equation} u(z)\approx u(L)L^\beta|z|^{-\beta}, \; \; z\rightarrow +\infty; \quad u(z)\approx u(-L)L^\beta |z|^{-\beta}, \; \; z\rightarrow -\infty, \end{equation} $

(4.1)

则$ \mathcal Q_3 $可近似为

$ \begin{eqnarray*} \mathcal Q_3 &\approx& \frac{L^\beta}{(\alpha+\beta)L^{\alpha+\beta}}\\&&\left[u(-L)\sideset{_2}{_1}{\mathop{\mathrm{H}}}(\beta, \alpha+\beta, \alpha+\beta+1, \frac{x}{L}) +u(L)\sideset{_2}{_1}{\mathop{\mathrm{H}}}(\beta, \alpha+\beta, \alpha+\beta+1, -\frac{x}{L})\right], \end{eqnarray*} $

其中的$ \sideset{_2}{_1}{\mathop{\mathrm{H}}}(\cdot, \cdot, \cdot, \cdot) $为超几何函数.

本节将通过数值实验来分析分数阶拉普拉斯算子$ (-\Delta)^{\alpha/2}u $的数值求积公式(2.27)的误差. 为此, 我们取$ u(x) = e^{-x^2} $, 它具有指数衰减性. 利用傅里叶变换容易得到

$ (-\Delta)^{\alpha/2}u(0) = \frac{1}{\sqrt{\pi}}\int_{0}^{\infty}x^\alpha e^{-\frac{x^2}{4}} = \frac{2^\alpha\Gamma(\frac{1+\alpha}{2})}{\sqrt{\pi}}, $

相应的数值结果见表 1.由该表的结果可知, 当$ \alpha $取为$ 1.8 $时, 数值求积公式(2.27)的收敛阶基本上达到$ O(h^{2.2}) $, 这与定理2.3中的结论是吻合的.

本节考虑一类具有齐次边界条件($ g\equiv 0 $)的分数阶拉普拉斯方程(3.1). 若取右端项$ f\equiv1 $, 则其精确解为

$ u(x) = \frac{2^{-\alpha}\Gamma(\frac12)}{\Gamma(1+\frac\alpha2)\Gamma(\frac{1+\alpha}2)} (1-x^2)_+^{\alpha/2}. $

我们使用有限差分格式(3.2)求解该问题, 其数值结果见表 2. 由该表的结果可知, 当$ \alpha $取为$ 1.8 $时, 有限差分格式(3.2)的收敛阶接近$ O(h^{2.2}) $, 这与定理3.2中的结论是吻合的.

本文主要讨论了一维分数阶拉普拉斯算子的数值积分公式, 并基于此构造了相应分数阶拉普拉斯方程的一种新型有限差分格式。相应的数值实验表明，相应的误差估计皆为最优。我们的下一步工作将把该方法推广到二维分数阶拉普拉斯算子和方程，即借助双线性插值在矩形网格上构造相应的数值积分公式和有限差分格式。