近年来，随着分数阶偏微分方程(FDEs)的研究受到越来越多的关注，证实了带分数阶导数的模型比传统的整数阶导数模型能更精确地描述科学与工程领域中的系统现象. 目前 FDEs已应用于许多方面，如电磁学^[1]、经济学^[2]、资产期权^[3]、物理学^[4]等. 求解FDEs有多种方法，对于具有简单边界条件的偏微分方程，解析解可以通过分离变量法^[5]或Laplace^[6]变换得到. 但因问题的复杂性，多数解析解难以求出,因此研究FDEs数值解显得非常必要.

对于时间分数阶微分方程的数值计算，目前已有不少研究成果. 其中L1式是时间Caputo分数阶导数的一个常用近似，具体可以参见文献[7-12]了解更多关于L1近似的应用信息. 文献[8]中利用时间有限差分方法和空间勒让德光谱法，建立了有限分差格式. Cui^[9]给出了方程的一个具有 $O\left( {\tau + {h^4}} \right)$ 精度的紧格式. 在文献[10-11]研究了一个具有非线性源项的修正异常时间分数阶慢扩散方程，提出了一阶时间精度和四阶空间精度的有限差分格式. 基于Grünwald-Letnikov近似，Yuste和Acdido^[13]构造了分数阶微分方程的显式差分格式，并利用冯诺依曼方法研究了稳定性. 而Zhuang等^[14]用能量法和隐式数值方法研究异常子渗流方程的稳定性及收敛性. 近来，Alikhanov^[15]提出了一个新的数值推导公式—— $L2 - {1_\sigma }$ 公式，在一个特殊的点上近似Caputo分数导数. 文献[15]和[16]构造了一种求解二维时空分数阶微分方程的高阶格式.

然而，上述工作均涉及常系数分数阶方程，通常为复杂模型的简化模型，而很多实际问题常用变系数方程来描述，如文献[17-20]. 因此研究变系数问题显得非常重要. 基于 $L2 - {1_\sigma }$ 离散公式，本文考虑带变系数的次扩散方程(见式(1))，并严格分析该差分格式的稳定性及收敛性。

$\begin{split}& {}_o^CD_t^\gamma u(x,t) = K(x)u_{xx}^{} + P(x)u_{xxx}^{} + Q(x)u_{xxxx}^{} + \\ &{\rm{ }}f(x,t),{\rm{ }}x \in {\bf{R}},{\rm{ }}t \in (0,T],{\rm{ }}0 < \gamma < 1. \end{split} $

(1)

方程(1)为满足如下周期边界条件的解：

$\left\{ \begin{split} & u(x,0) = \varphi (x),{\rm{ }}x \in {\bf{R}}; \\ & u(x,t) = u(x + L,t),{\rm{ }}x \in {\bf{R}},{\rm{ }}t \in (0,T]. \end{split} \right.$

(2)

其中 $K(x),{\rm{ }}P(x),{\rm{ }}Q(x)$ 足够光滑，并且满足 $K(x) > 0, $ $Q(x) < 0.$ Caputo分数阶导数的定义：

${}_o^CD_t^\gamma u(x,t) = \frac{1}{{\Gamma (1 - \gamma )}}\int\nolimits_0^t {\frac{{{\text{∂}} u(x,s)}}{{{\text{∂}} s}}\frac{{{\rm d}s}}{{{{(t - s)}^\gamma }}}} .$

本文构局如下：第2节，基于文献[15]中所提出的 $L2 - {1_\sigma }$ 离散公式，构建了有限差分格式；第3节，证明该差分格式的唯一可解性，并给出局部截断误；第4节，用Fourier方法研究了该格式的稳定性和收敛性；第5节中给出了一些与理论分析一致的数值结果；最后得出简要结论.

1 差分格式的建立及分析

给定正整数 $N,{\rm{ }}M,$ 记 $\tau = T/N,$ $h = L/M,$ 时间方向的网格剖分为 ${\Omega _\tau } = \{ {\rm{ }}{t_n} = n\tau ,{\rm{ }}0 \leqslant n \leqslant N\} .$ 记 $ {{\rm{\mathcal{V}}}_\tau } = \{ {\rm{ }}v|v =$ $ ({v^{\rm{0}}},{v^{\rm{1}}}, \cdots ,{v^N})\}$ , 为定义在 ${\Omega _\tau }$ 上的网格函数空间. 对任意 $v \in {{\rm{\mathcal{V}}}_\tau }$ ,空间上的网格剖分为 ${\Omega _h} = \{ {x_i}|0 \leqslant i \leqslant M\} .$ 记 ${{\rm{\mathcal{V}}}_h} = \{ u|u = ({u_0},{u_1}, \cdots ,$ ${u_M})\} $ 为定义在 ${\Omega _h}$ 上网格函数空间. 考虑到周期边界条件有 ${\rm{ }}u_i^n = u_{i \pm M}^n.$

本文给出了一个序列 $\left\{ {c_m^{(k)}} \right\}$ ，如文献[15]中的定义，假设 $0 < \gamma < 1,{\rm{ }}\sigma = 1 - \displaystyle\frac{\gamma }{2}$ ，则有，

$\begin{gathered} {a_0} = {\sigma ^{1 - \gamma }},{\rm{ }}{a_l} = {\left( {l + \sigma } \right)^{1 - \gamma }} - {\left( {l - 1 + \sigma } \right)^{1 - \gamma }},{\rm{ }}l \geqslant 1 ; \\ {b_l} = \frac{1}{{2 - \gamma }}\left[ {{{\left( {l + \sigma } \right)}^{2 - \gamma }} - {{\left( {l - 1 + \sigma } \right)}^{2 - \gamma }}} \right]- \\ \frac{1}{2}\left[ {{{\left( {l + \sigma } \right)}^{1 - \gamma }} - {{\left( {l - 1 + \sigma } \right)}^{1 - \gamma }}} \right],{\rm{ }}l \geqslant 1. \end{gathered} $

当 $k = 0$ 时， $c_0^{\left( k \right)} = {a_0}$ ；当 $k \geqslant 1$ 时，

$c_m^{\left( k \right)} = \left\{ \begin{array}{l} {a_0} + {b_1},{\rm{ }}m = 0; \\ {a_m} + {b_{m + 1}} - {b_m},{\rm{ }}1 \leqslant m \leqslant k - 1; \\ {a_k} - {b_k},{\rm{ }}m = k. \\ \end{array} \right.$

它的定义符合

$\left\{ \begin{split} & c_m^{\left( k \right)} < \displaystyle\frac{{1 - \gamma }}{2}{\left( {m + \sigma } \right)^{ - \gamma }}, \\ & c_0^{\left( k \right)} > c_2^{\left( k \right)} > c_2^{\left( k \right)} > \cdots > c_k^{\left( k \right)}. \end{split} \right.$

(3)

当 $0 < \gamma < 1{\rm{ }},{\rm{ }}u(x) \in {C^2}[0,T]$ 时，

${}_0^CD_t^\gamma u({t_{n + \sigma }}) - \Delta _t^\gamma u({t_{n + \sigma }}) = O({\tau ^{3 - \gamma }}),$

其中

$\Delta _t^\gamma u({t_{n + \sigma }}) = \frac{{{\tau ^{ - \gamma }}}}{{\Gamma (2 - \gamma )}}\sum\limits_{k = 0}^n {c_k^{(n + 1)}} [u({t_{n - k}}) - u({t_{n - k + 1}})].$

引理1 　假设 $U(x) \in {C^6}[{x_{i - 2}},{x_{i + 2}}]$ ，定义如下3个算子 ${{\rm{\mathcal{H}}}_{{\rm{ }}1}},{{\rm{\mathcal{H}}}_{{\rm{ }}2}},{{\rm{\mathcal{H}}}_{{\rm{ }}3}},$ 通过Taylor展式运算得：

$\left\{\begin{array}{l} {{\rm{\mathcal{H}}}_1}{U_i} = \displaystyle\frac{1}{{{h^2}}}({U_{i - 1}} - 2{U_i} + {U_{i + 1}}), \\ {{\rm{\mathcal{H}}}_2}{U_i} = \displaystyle\frac{1}{{{h^3}}}(\frac{1}{2}{U_{i - 2}} + {U_{i - 1}} - {U_{i + 1}} + \frac{1}{2}{U_{i + 2}}),{\rm{ }} \\ {{\rm{\mathcal{H}}}_3}{U_i} = \displaystyle\frac{1}{{{h^4}}}({U_{i - 2}} - 4{U_{i - 1}} + 6{U_i} - 4{U_{i + 1}} + {U_{i + 2}}). \end{array}\right. $

(4)

$\left\{\begin{gathered} \left| {{u_{xx}}({x_i}) - {{\rm{\mathcal{H}}}_1}{U_i}} \right| \leqslant \frac{1}{{12}}\mathop {\max }\limits_{{x_{i - 2}} \leqslant x \leqslant {x_{i + 2}}} \left| {u_x^{(4)}(x)} \right|{h^2}, \\ \left| {{u_{xxx}}({x_i}) - {{\rm{\mathcal{H}}}_2}{U_i}} \right| \leqslant \frac{1}{4}\mathop {\max }\limits_{{x_{i - 2}} \leqslant x \leqslant {x_{i + 2}}} \left| {u_x^{(5)}(x)} \right|{h^2}, \\ \left| {{u_{xxxx}}({x_i}) - {{\rm{\mathcal{H}}}_3}{U_i}} \right| \leqslant \frac{1}{6}\mathop {\max }\limits_{{x_{i - 2}} \leqslant x \leqslant {x_{i + 2}}} \left| {u_x^{(6)}(x)} \right|{h^2}. \end{gathered} \quad\quad\right.$

(5)

引理2 　根据文献[17-19]，假设 $0 < \gamma < 1{\rm{ }},\;u(t) \in $ ${C^2}[0,T],$ 则有

$u({t_{n + \sigma }}) = \sigma u({t_{n + 1}}) + (1 - \sigma )u({t_n}) + O({\tau ^2}).$

(6)

定义 ${u^{({\sigma _{k + 1}})}} = \sigma u({t_{k + 1}}) + (1 - \sigma )u({t_k})$ . 根据引理1和引理2，并用数值解 $u_i^n$ 代替精确解 $U_i^n$ ，可以得到如下差分格式：

$\left\{\begin{array}{l} \Delta _t^\gamma u_i^{n + \sigma } = {K_i}{{\rm{\mathcal{H}}}_1}u_i^{({\sigma _{n + 1}})} + {P_i}{{\rm{\mathcal{H}}}_2}u_i^{({\sigma _{n + 1}})} + \\ {\rm{ }}{Q_i}{{\rm{\mathcal{H}}}_3}u_i^{({\sigma _{n + 1}})} + f_i^{n + \sigma },{\rm{ }}0 \leqslant n \leqslant N,{\rm{ }}1 \leqslant i \leqslant M; \\ u_i^0 = \varphi ({x_i}),{\rm{ }}1 \leqslant i \leqslant M; \\ u_i^{n + 1} = u_{i \pm M}^{n + 1},{\rm{ }}1 \leqslant i \leqslant M,{\rm{ }}0 \leqslant n \leqslant N. \end{array}\quad\quad\right. $

(7)

基于式(5)和(6)，可以很容易得到该差分格式的截断误差为 $O({\tau ^2} + {h^2}).$

2 差分格式的唯一可解性

在这一部分，用Fourier方法讨论差分格式(7)是唯一可解的. 首先给出一些Fourier的定义，定义 ${\nu ^{n + 1}}(x),{\rm{ }}{z^{n + \sigma }}(x)$ 为具有周期为 $L$ 的2个周期函数，其中 $[0,L]$ 为它们的周期：

${\nu ^{n + 1}}(x) = \left\{ \begin{array}{l} u_0^{n + 1},x \in [{x_0},{x_{1/2}}]; \\ u_i^{n + 1},x \in ({x_{i - 1/2}},{x_{i + 1/2}}],{\rm{ }}i = 1,2, \cdots ,M - 1; \\ u_M^{n + 1},x \in [{x_{M - 1}},{x_M}]. \end{array} \right.$

${z^{n + \sigma }}(x) = \left\{ \begin{array}{l} f_0^{n + \sigma },x \in [{x_0},{x_{1/2}}]; \\ f_i^{n + \sigma },x \in ({x_{i - 1/2}},{x_{i + 1/2}}],{\rm{ }}i = 1,2, \cdots ,M - 1;\\ f_M^{n + \sigma },x \in [{x_{M - 1}},{x_M}]. \end{array}\right.$

由Fourier 级数，定义 ${\nu ^{n + 1}}(x)$ 和 ${z^{n + \sigma }}(x)$ 如下：

${\nu ^{n + 1}}(x) = \sum\limits_{l = - \infty }^\infty {\tilde v_l^{n + 1}} {{\rm e}^{\beta x{\rm j}}},$

${z^{n + \sigma }}(x) = \sum\limits_{l = - \infty }^\infty {\tilde z_l^{n + \sigma }} {{\rm e}^{\beta x{\rm j}}}.$

其中

${\rm{ }}{{\rm j}^2} = - 1,{\rm{ }}\beta = \frac{{2{\text{π}} l}}{L},$

$\tilde v_l^{n + 1}(x) = \frac{1}{L}\int_0^L {{v^{n + 1}}(x)} {\rm{ }}{{\rm e}^{ - \beta x{\rm j}}}{\rm d}x,$

$\tilde z_l^{n + \sigma }(x) = \frac{1}{L}\int_0^L {{z^{n + \sigma }}(x){\rm{ }}} {{\rm e}^{ - \beta x{\rm j}}}{\rm d}x.$

定义 ${L_2}$ 范数为

${\left\| {{u^k}} \right\|_2} = {\left( {h\sum\limits_{i = 1}^M {{{\left| {u_i^k} \right|}^2}} } \right)^{\frac{1}{2}}}.$

根据 Parseval 等式可以得到如下等式：

${\left\| {{u^k}} \right\|_2} = {\left( {\int_0^L {{{\left| {{v^k}(x)} \right|}^2}{\rm d}x} } \right)^{\frac{1}{2}}} = {\left( {\sum\limits_{l = - \infty }^\infty {{{\left| {\tilde v_l^k(x)} \right|}^2}} } \right)^{\frac{1}{2}}},$

${\left\| {{f^{n + \sigma }}} \right\|_2} = {\left( {\int_0^L {{{\left| {{z^{n + \sigma }}(x)} \right|}^2}{\rm d}x} } \right)^{\frac{1}{2}}} = {\left( {\sum\limits_{l = - \infty }^\infty {{{\left| {\tilde z_l^{n + \sigma }(x)} \right|}^2}} } \right)^{\frac{1}{2}}}.$

现开始证明差分格式(7)的唯一可解性，先不考虑周期边界问题，通过计算，差分格式(7)可以等价为如下形式：

$\begin{split} & \left[ {sc_0^{\left( {n + 1} \right)} - \sigma {\rm{\mathcal{H}}}} \right]u_i^{n + 1} = \\ & \sum\limits_{k = 0}^N {g_k^{n + 1}u_i^k + 2{h^4}f_i^{n + \sigma } + } (1 - \sigma ){\rm{\mathcal{H}}}u_i^n,{\rm{ }} \\ & {\rm{ }}1 \leqslant i \leqslant M,{\rm{ }}0 \leqslant n \leqslant N. \end{split} $

(8)

其中定义

$\begin{split} & sc_0^{\left( {n + 1} \right)} = \frac{{2{h^4}{\tau ^{ - \gamma }}}}{{\Gamma (2 - \gamma )}},{\rm{ }}g_0^{n + 1} = sc_n^{(n + 1)},{\rm{ }} \\ & g_k^{n + 1} = s\left[ {c_{n - k}^{(n + 1)} - c_{n - k + 1}^{(n + 1)}} \right],{\rm{ }}k = 0,1,2, \cdots ,n. \end{split} $

$\begin{split} &{\rm{\mathcal{H}}}u_i^{n + 1} = 2{h^4}({K_i}{{\rm{\mathcal{H}}}_1} + {P_i}{{\rm{\mathcal{H}}}_2} + {Q_i}{{\rm{\mathcal{H}}}_3})u_i^{n + 1} = \\ & \left[ {{\lambda _{i.i - 2}}u_{i - 2}^{n + 1} + {\lambda _{i.i - 1}}u_{i - 1}^{n + 1} + {\lambda _{i.i}}u_i^{n + 1}} \right. + \left. { {\lambda _{i.i + 1}}u_{i + 1}^{n + 1} + {\lambda _{i.i + 2}}u_{i + 2}^{n + 1}} \right]. \end{split} $

(9)

将通过Taylor展式运算得到的式(4)代入式(9)中，通过计算，可以得到:

$\left\{ \begin{split} &{\lambda _{i,i}} = - 4{h^2}{K_i} + 12{Q_i}, \\ &{\lambda _{i,i + 1}} = 2{h^2}{K_i} - 2h{P_i} - 8{Q_i}, \\ & {\lambda _{i,i + 2}} = h{P_i} + 2{Q_i}, \\ & {\lambda _{i,i - 1}} = 2{h^2}{K_i} + 2h{P_i} - 8{Q_i}, \\ & {\lambda _{i,i - 2}} = - h{P_i} + 2{Q_i}. \end{split} \right.$

由Fourier 级数的定义和式(8)，得到如下方程：

$\begin{split} & \sum\limits_{l = - \infty }^\infty {\left\{\left( {sc_0^{(n + 1)} - \sigma {\rm{\mathcal{H}}}} \right) {{{\rm e}^{\beta x{\rm j}}}} \right\}\tilde v_l^{n + 1}} = \\ & \sum\limits_{l = - \infty }^\infty {\left[ {\sum\limits_{k = 0}^n {g_k^{n + 1}\tilde v_l^n{{\rm e}^{\beta x{\rm j}}} + 12{h^4}\tilde z_l^{n + \sigma }{{\rm e}^{\beta x{\rm j}}}} } \right]} + \\ & \sum\limits_{l = - \infty }^\infty {(1 - \sigma ){\rm{\mathcal{H}}}{{\rm e}^{\beta x{\rm j}}}\tilde v_l^n}.\end{split} $

(10)

其中定义

$\begin{split} & {\rm{\mathcal{H}}}{{\rm e}^{\beta x{\rm j}}} = {{\rm e}^{\beta x{\rm j}}}\left[ {{\lambda _{i,i - 2}}{{\rm e}^{ - 2\beta x{\rm j}}} + {\lambda _{i,i - 1}}{{\rm e}^{ - \beta x{\rm j}}} + {\lambda _{i,i}}} \right. + \\ & \left. {{\lambda _{i,i + 1}}{{\rm e}^{\beta x{\rm j}}} + {\lambda _{i,i + 2}}{{\rm e}^{2\beta x{\rm j}}}} \right] = {{\rm e}^{\beta x{\rm j}}}({y_1} + {y_2}{\rm j}). \end{split} $

(11)

通过进一步计算求出

$\begin{split} &{y_1} = 4{h^2}[\cos (\beta h) - 1]{K_i} + {\rm{8}}[{\cos ^2}(\beta h) - 2\cos (\beta h) + 1]{Q_i}. \\ &{y_2} = 2h\left[ { - 2\sin (\beta h) + \sin (2\beta h)} \right]{P_i}{\rm j}. \end{split}$

由式(10)和式(11)可以推算得到

$\begin{split} &\sum\limits_{l = - \infty }^\infty {\left[ {sc_0^{(n + 1)} - \sigma ({y_1} + {y_2}{\rm j})} \right]} {\rm{ }}\tilde v_l^{n + 1}{{\rm e}^{\beta x{\rm j}}} = \\ & \sum\limits_{l = - \infty }^\infty {\left[ {\sum\limits_{k = 0}^n {g_k^{n + 1}\tilde v_l^k + 2{h^4}} \tilde z_l^{n + \sigma }} \right]} {\rm{ }}{{\rm e}^{\beta x{\rm j}}} + \\ & \sum\limits_{l = - \infty }^\infty {(1 - \sigma )({y_1} + {y_2}{\rm j})} {\rm{ }}\tilde v_l^n{{\rm e}^{\beta x{\rm j}}}. \end{split} $

(12)

将式(12)的两边同时乘以 ${{\rm e}^{ - \beta x{\rm j}}}$ ，再0到L积分，可以得到

$\begin{split} & \left[ {sc_0^{(n + 1)} - \sigma ({y_1} + {y_2}{\rm j})} \right]{\rm{ }}\tilde v_l^{n + 1} = \\ & \sum\limits_{k = 0}^n {g_k^{n + 1}\tilde v_l^k + 2{h^4}} \tilde z_l^{n + \sigma } + (1 - \sigma )({y_1} + {y_2}{\rm j})\tilde v_l^n, \\ & l = 0, \pm 1, \pm 2, \cdots , \pm \infty . \end{split} $

(13)

为了证明差分格式(7)是唯一可解的，只需要证明下面方程是唯一可解即可：

$\left\{ \begin{split}& \left( {sc_0^{(n + 1)} - \sigma {\rm{\mathcal{H}}}} \right)u_i^{n + 1} = 0,{\rm{ }}0 \leqslant n \leqslant N; \\ & u_i^0 = \varphi ({x_i}),{\rm{ }}1 \leqslant i \leqslant M; \\ & u_i^{n + 1} = u_{i \pm M}^{n + 1},{\rm{ }}1 \leqslant i \leqslant M,{\rm{ }}0 \leqslant n \leqslant N. \end{split} \right.$

(14)

因为 ${K_i} < 0,{\rm{ }}{Q_i} < 0$ ，且 ${y_1}$ 不随 $\cos (\beta h)$ 递增，所以求出 ${y_1} < 0$ . 因此由式(14)可以得到

$\begin{split} & \left( {sc_0^{(n + 1)} - \sigma ({y_1} + {y_2}{\rm j}} \right)\tilde v_l^{n + 1} = 0,{\rm{ }} \\ &l = 0, \pm 1, \pm 2, \cdots , \pm \infty. \end{split} $

(15)

由上面证明，结合式(15)，可知 $\tilde v_l^{n + 1} = 0,$ 因为 $\left| {sc_0^{(n + 1)} - \sigma ({y_1} + {y_2}{\rm j})} \right| = \sqrt {{{\left( {sc_0^{(n + 1)} - \sigma {y_1}} \right)}^2} + {{(\sigma {y_2})}^2}} $ $ > 0.$ 所以可以得到 ${v^{n + 1}}(x) = 0$ ，进而知 ${\left\| {{u^{n + 1}}} \right\|_2} = 0$ ,此时，由以上证明可知差分格式(7)是唯一可解的. 定理证毕.

3 差分格式的稳定性分析

这一部分分析差分格式 $(7)$ 的稳定性和收敛性.

定理1 　对于方程 $(13)$ 可以得到：

$\left| {\tilde v_l^{n + 1}} \right| \leqslant \left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n} \left| {\tilde z_l^{m + \sigma }} \right|.$

(16)

证明 　用数学归纳法证明方程(16)：当 $n = 0$ , 基于式(13), 可以得到

$\begin{split} &{\rm{ }}\left( {sc_0^{(1)} - \sigma ({y_1} + {y_2}{\rm j}} \right)\tilde v_l^1 = \\ & g_0^1\tilde v_l^0 + 2{h^4}\tilde z_l^\sigma + (1 - \sigma )({y_1} + {y_2}{\rm j})\tilde v_l^0 = \\ & sc_0^{\left( 1 \right)}\tilde v_l^0 + 2{h^4}\tilde z_l^\sigma + (1 - \sigma )({y_1} + {y_2}{\rm j})\tilde v_l^0. \end{split} $

注意到 ${y_1} < 0,{\rm{ }}\sigma > \displaystyle\frac{1}{2},$ 结合式 $(3)$ 的结论(参见文献[15])，可以推测：

$\begin{split} & \left| {\tilde v_l^1} \right| = \left| {\frac{{sc_0^{\left( 1 \right)}\tilde v_l^0 + 2{h^4}\tilde z_l^\sigma + (1 - \sigma )({y_1} + {y_2}{\rm j})\tilde v_l^0}}{{sc_0^{(1)} - \sigma ({y_1} + {y_2}){\rm j}}}} \right| \leqslant\\ & \sqrt {\frac{{{{\left[ {sc_0^{\left( 1 \right)} + (1 - \sigma ){y_1}} \right]}^2} + {{(1 - \sigma )}^2}{y_2}^2}}{{{{\left( {sc_0^{(1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} \left| {\tilde v_l^0} \right| + \\ & {\rm{ }}\frac{{2{h^4}}}{{\sqrt {{{\left( {sc_0^{(1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2} }}\left| {\tilde z_l^\sigma } \right| = \\ & \sqrt {\frac{{{{\left( {sc_0^{(1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}{{{{\left( {sc_0^{(1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} \left| {\tilde v_l^0} \right| + \\ & \sqrt {\frac{{\left[ {2sc_0^{(1)} + (1 - 2\sigma ){y_1}} \right]{y_1} + (1 - 2\sigma ){y_2}^2}}{{{{\left( {sc_0^{(1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} \left| {\tilde v_l^0} \right| + \\ & \frac{{\Gamma (2 - \gamma ){\tau ^\gamma }}}{{c_0^{(1)}}}\left| {\tilde z_l^\sigma } \right|{\rm{ }} = \sqrt {1 + \frac{{2sc_0^{\left( 1 \right)}{y_1} + (1 - 2\sigma )({y_1}^2 + {y_2}^2)}}{{{{\left( {sc_0^{(1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} {\rm{ }}\left| {\tilde v_l^0} \right| +\\ & \frac{{\Gamma (2 - \gamma ){\tau ^\gamma }}}{{c_0^{(1)}}}\left| {\tilde z_l^\sigma } \right|{\rm{ }} {\rm{ }} \leqslant \left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\left| {\tilde z_l^\sigma } \right|. \end{split} $

因此 $(16)$ 对于 $0,1,2, \cdots ,n - 1$ 均成立，现在讨论 $n$ 时

$\begin{split} &\left| {\tilde v_l^{n + 1}} \right| = \left| {\displaystyle\frac{{\sum\limits_{k = 0}^n {g_k^{n + 1}} \tilde v_l^k + 2{h^4}\tilde z_l^{n + \sigma } + (1 - \sigma )({y_1} + {y_2}{\rm j})\tilde v_l^n}}{{sc_0^{(n + 1)} - \sigma ({y_1} + {y_2}){\rm j}}}} \right|\leqslant\\ & \sqrt {\displaystyle\frac{{{{\left[ {\sum\limits_{k = 0}^n {\left| {g_k^{n + 1}} \right|} + (1 - \sigma ){y_1}} \right]}^2} + {{(1 - \sigma )}^2}{y_2}^2}}{{{{\left( {sc_0^{(n + 1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} \times \\ &\left[ {\left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n - 1} \left| {\tilde z_l^{m + \sigma }} \right|} \right] + \\ &\displaystyle\frac{{2{h^4}\left| {\tilde z_l^{n + \sigma }} \right|}}{{\sqrt {{{\left( {sc_0^{(n + 1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2} }}. \end{split}$

因为 $\sum\limits_{k = 0}^n {\left| {g_k^{n + 1}} \right|} = sc_0^{(n + 1)},$ 再结合式(3)可以得到：

$\begin{split} &\left| {\tilde v_l^{n + 1}} \right| \leqslant \sqrt {\displaystyle\frac{{{{\left[ {sc_0^{\left( {n + 1} \right)} + (1 - \sigma ){y_1}} \right]}^2} + {{(1 - \sigma )}^2}{y_2}^2}}{{{{\left( {sc_0^{(n + 1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} \times \\ &{\rm{ }}\left[ {\left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n - 1} \left| {\tilde z_l^{m + \sigma }} \right|} \right] + \\ &{\rm{ }}\displaystyle\frac{{2{h^4}\left| {\tilde z_l^{n + \sigma }} \right|}}{{\sqrt {{{\left( {sc_0^{(n + 1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2} }}=\\ &{\rm{ }} \sqrt {1 + \displaystyle\frac{{2sc_0^{\left( {n + 1} \right)}{y_1} + (1 - 2\sigma )({y_1}^2 + {y_2}^2)}}{{{{\left( {sc_0^{(n + 1)} - \sigma {y_1}} \right)}^2} + {\sigma ^2}{y_2}^2}}} \times \\ &{\rm{ }}\left[ {\left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n - 1} \left| {\tilde z_l^{m + \sigma }} \right|} \right] + \\ &{\rm{ }}\displaystyle\frac{{\Gamma (2 - \gamma ){\tau ^\gamma }}}{{c_0^{(n + 1)}}}\left| {\tilde z_l^{n + \sigma }} \right|\leqslant \\ &{\rm{ }} \left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n} \left| {\tilde z_l^{m + \sigma }} \right|. \end{split}$

定理1证毕.

定理2 　结合方程(16)，可证差分格式收敛，即式(17)的估计成立.

${\left\| {{u^{n + 1}}} \right\|_2} \leqslant {\left\| {{u^0}} \right\|_2} + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n} {\left\| {{f^{m + \sigma }}} \right\|_2}.$

(17)

证明 　根据定理1和 Minkowski 不等式可以推测：

$\begin{split}& {\rm{ }}{\left[ {\sum\limits_{l = - \infty }^\infty {{{\left| {\tilde v_l^{n + 1}} \right|}^2}} } \right]^{\frac{1}{2}}}\leqslant\\ & {\left[ {\sum\limits_{l = - \infty }^\infty {{{\left( {\left| {\tilde v_l^0} \right| + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n} \left| {\tilde z_l^{m + \sigma }} \right|} \right)}^2}} } \right]^{\frac{1}{2}}} \leqslant \\ & {\left( {\sum\limits_{l = - \infty }^\infty {{{\left| {\tilde v_l^0} \right|}^2}} } \right)^{\frac{1}{2}}} + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n} {\left( {\sum\limits_{l = - \infty }^\infty {{{\left| {\tilde z_l^{m + \sigma }} \right|}^2}} } \right)^{\frac{1}{2}}}. \end{split} $

结合Parseval等式，证明 $(17)$ 成立.

定理3 　记 $U_i^{n + 1} = {[U_1^{n + 1},U_2^{n + 1}, \cdots ,U_m^{n + 1}]^T}$ ，其中 $0 \leqslant n \leqslant N,$ 为方程(1)~(2)的解， $u_i^{n + 1} = [u_1^{n + 1},$ $u_2^{n + 1}, \ldots ,u_m^{n + 1}{]^{\rm T}},{\rm{ }}0 \leqslant n \leqslant N,$ 为差分格式(7)的解，记 ${\rm e}_i^{n + 1} = U_i^{n + 1} - u_i^{n + 1},{\rm{ }}0 \leqslant n \leqslant N,{\rm{ }}1 \leqslant i \leqslant M.$ 存在下面的 ${c_1},{\rm{ }}{c_2}$ ，使得

$\left\| {{{\rm e}^{n + 1}}} \right\| \leqslant {c_1}{\tau ^2} + {c_2}{h^2},{\rm{ }}0 \leqslant n \leqslant N.$

证明 　由误差方程可以得到

$\begin{split} &\left[ {sc_0^{\left( {n + 1} \right)} - \sigma {\rm{\mathcal{H}}}} \right]{\rm e}_i^{n + 1} = \\ & \sum\limits_{k = 0}^n {g_k^{n + 1}{\rm e}_i^k} + 2{h^4}R_i^{n + \sigma } + (1 - \sigma ){\rm{\mathcal{H}}}{\rm e}_i^{n - 1}, \\ &1 \leqslant i \leqslant M,{\rm{ }}0 \leqslant n \leqslant N. \\ \end{split} $

其中 $R_i^{n + \sigma } = O({\tau ^2} + {h^2}).$

根据之前的分析和定理3可以得到：

$\begin{split} & {\left\| {{{\rm e}^{n + 1}}} \right\|_2} \leqslant {\left\| {{{\rm e}^0}} \right\|_2} + 2\Gamma (1 - \gamma ){T^\gamma }\mathop {\max }\limits_{0 \leqslant m \leqslant n} {\left\| {{R^{m + \sigma }}} \right\|_2} \leqslant \\ & {c_1}{\tau ^2} + {c_2}{h^2},{\rm{ }}0 \leqslant n \leqslant N. \\ \end{split} $

其中

${c_1} = \frac{{(1 - \sigma )\sigma }}{2}\mathop {\max }\limits_{{t_0} \leqslant t \leqslant {t_n}} \left| {u_t^{(2)}(x,t)} \right|,$

$\begin{gathered} {c_2} = 2\Gamma (1 - \gamma ){T^\gamma }\left( {\frac{1}{{12}}\mathop {\max }\limits_{{x_{i - 2}} \leqslant x \leqslant {x_{i + 2}}} \left| {u_x^{(4)}(x,t)} \right|{h^2}} \right. + \\ \left. {{\rm{ }}\frac{1}{4}\mathop {\max }\limits_{{x_{i - 2}} \leqslant x \leqslant {x_{i + 2}}} \left| {u_x^{(5)}(x,t)} \right|{h^2} + \frac{1}{6}\mathop {\max }\limits_{{x_{i - 2}} \leqslant x \leqslant {x_{i + 2}}} \left| {u_x^{(6)}(x,t)} \right|{h^2}} \right). \\ \end{gathered} $

定理2证毕.

4 数值实验

本节考虑了带空间变系数的时间分数阶偏微分方程的周期函数问题，算例验证了差分格式 $(7)$ 对Dirichelet边界值的周期函数问题的有效性. 记

${E_2}(\tau ,h) = {\left\| {U_i^{n + 1} - u_i^{n + 1}} \right\|_2},$

时间方向和空间方向的收敛阶分别为 $ {\rm Rate}_\tau ,$ ${\rm Rate}{_h}$ ：

$\begin{gathered} {{\rm Rate}_\tau } = {\log _2}\left( {\frac{{{E_2}(h,2\tau )}}{{{E_2}(h,\tau )}}} \right), \\ {\rm Rate}_h = {\log _2}\left( {\frac{{{E_2}(2h,\tau )}}{{{E_2}(h,\tau )}}} \right). \\ \end{gathered} $

当计算时间方向的收敛阶 ${\rm{Rat}}{{\rm{e}}_\tau }$ 时，我们要求 $h$ 足够小使得空间方向的误差不影响整体误差. 反之，当计算空间方向的收敛阶 ${\rm{Rat}}{{\rm{e}}_h}$ 时，我们要求 $\tau $ 足够小.

首先用差分格式(7)计算有已知精确解的问题(1)~(3)以验证理论结果.

例1 　在 $[0,1] \times [0,1]$ 上考虑问题(1)~(3)精确解.

$u(x,t) = ({t^2} + 1)\sin (2{\text{π}} x),{\rm{ }}x \in {\bf{R}}.$

系数

$K(x) = {x^4} + 0.01,{\rm{ }}P(x) = - {x^4},{\rm{ }}Q(x) = - ({x^4} + 0.01),$

源项

$\begin{split} & f(x,t) = - 8{{\text{π}} ^3}{x^4}({t^2} + 1)\cos 2{\text{π}} x + \\ & \left[ {4{{\text{π}} ^2}({t^2} + 1)({x^4} + 0.01)(1 + 4{{\text{π}} ^2}) + \frac{{2{t^{2 - \gamma }}}}{{\Gamma (3 - \gamma )}}} \right]\sin 2{\text{π}} x. \end{split} $

表1中列出了当空间步长 $h = 1/1\;500$ 分别对 $\gamma = 0.2,{\rm{ }}0.5,{\rm{ }}0.8$ 不同时间步长下的误差及收敛阶. 容易看出数值收敛阶验证了差分格式在时间方向的收敛阶为 $2$ .

当固定时间步长 $\tau = 1/5\;000$ ，表2中给出了当 $\gamma = 0.5$ 时不同空间步长下的误差及收敛结果. 从表2可以看出，该算例验证了空间方向的二阶收敛精度.

5 结论

本文主要研究带空间变系数的时间分数阶扩散方程的差分方法. 时间方面应用 $L2 - {1_\sigma }$ 离散公式以达到二阶精度，空间方面采用二阶差分格式. 证明了该差分格式的无条件稳定性和收敛性,并通过数值实验来验证理论结果.