$ 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}} $
$ {\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\}, $
$ \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} = \frac{\alpha2^{\alpha-1}\Gamma(\frac{\alpha+n}{2})} {\pi^\frac{n}{2}\Gamma(\frac{2-\alpha}{2})}, $
$ \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)
$ \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)
2.1 奇异部分
$ \mathcal{Q}^S $
的离散及误差估计
设$ 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) 直接推出.
2.2 正则部分$ \mathcal{Q}^Ru(x) $的求积公式及误差估计
在每个子区间
$ [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.3 分数阶拉普拉斯算子的数值求积公式
联立(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)
3 分数阶拉普拉斯方程的有限差分格式及其误差分析
本节将基于分数阶拉普拉斯算子(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}. $
定理得证.
4 数值实验
注意到分数阶拉普拉斯算子的积分区域为整个实数域, 在实际计算中需要对其进行截断, 为此我们将其改写为如下形式
$ \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) $为超几何函数.
4.1 分数阶拉普拉斯算子的数值求积公式
本节将通过数值实验来分析分数阶拉普拉斯算子
$ (-\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中的结论是吻合的.
4.2 分数阶拉普拉斯方程的有限差分格式
本节考虑一类具有齐次边界条件($ 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中的结论是吻合的.
Focardi S, Cecere J, Hays G.
The Lévy flight foraging hypothesis in a pelagic seabird[J]. Journal of Animal Ecology, 2014, 83(2): 353–364.
DOI:10.1111/1365-2656.12147
Benson D, Wheatcraft S, Meerschaert M.
Application of a fractional advection-dispersion equation[J]. Water Resources Research, 2000, 36(6): 1403–1412.
DOI:10.1029/2000WR900031
Cushman J, Ginne T.
Nonlocal dispersion in media with continuously evolving scales of heterogeneity[J]. Transport in Porous Media, 1993, 13(1): 123–138.
DOI:10.1007/BF00613273
Buades A, Coll B, Morel J.
Image denoising methods. A new nonlocal principle[J]. SIAM Review, 2010, 52(1): 113–147.
DOI:10.1137/090773908
Butzer P and Westphal U. An introduction to fractional calculus[C]//Hilfer R. Applications of Fractional Calculus in Physics. World Scientific, 2000: 1-85.
Carreras B, Lynch V, Zaslavsky G.
Anomalous diffusion and exit time distribution of particle tracers in plasma turbulence models[J]. Physics of Plasmas, 2001, 8(12): 5096–5103.
DOI:10.1063/1.1416180
Del-Castillo-Negrete D, Carreras B, Lynch V.
Fractional diffusion in plasma turbulence[J]. Physics of Plasmas, 2004, 11(8): 3854–3864.
DOI:10.1063/1.1767097
Applebaum D.
Lévy processes and stochastic calculus[M]. Cambridge: Cambridge University Press, 2009.
Saichev A, Zaslavsky G.
Fractional kinetic equations: solutions and applications[J]. Chaos, 1997, 7(4): 753–764.
DOI:10.1063/1.166272
Nezza E, Palatucci G, Valdinoci E.
Hitchhiker's guide to the fractional Sobolev spaces[J]. Bulletin des Sciences Mathématiques, 2012, 136(5): 521–573.
DOI:10.1016/j.bulsci.2011.12.004
Acosta G, Borthagaray J.
A fractional Laplace equation: regularity of solutions and finite element approximations[J]. SIAM Journal on Numerical Analysis, 2015, 55(2): 472–495.
Huang Y, Oberman A.
Numerical methods for the fractional Laplacian: A finite differencequadrature approach[J]. SIAM Journal on Numerical Analysis, 2014, 52(6): 3056–3084.
DOI:10.1137/140954040
Duo S, Wyk H, Zhang Y.
A novel and accurate finite difference method for the fractional Laplacian and the fractional Poisson problem[J]. Journal of Computational Physics, 2017, 355(15): 233–252.
Biler P, Funaki T, Woyczynski W.
Fractal burgers equations[J]. Journal of Differential Equations, 1998, 148(1): 9–46.
DOI:10.1006/jdeq.1998.3458
Pablo A, Quirós F, Rodríguez A, Vázquez J.
A fractional porous medium equation[J]. Advances in Mathematics, 2011, 226(2): 1378–1409.
DOI:10.1016/j.aim.2010.07.017
