空间分布阶时间分数阶扩散方程的有限体积法

Finite Volume Method for a Space Distributed-Order Time-Fractional Diffusion Equation

摘要: 本文利用有限体积法研究了空间分布阶时间分数阶扩散方程。首先，用中点求积法将空间分布阶项转化为多项空间分数阶项，利用有限体积法对多项空间分数阶项进行离散。而对于时间分数阶导数，我们采用有限差分法。其次，我们证明了迭代格式的无条件稳定性和收敛性。最后通过一个数值例子来证明算法的有效性。

Abstract: In this paper, the space distributed-order time-fractional diffusion equation is considered. We propose a finite volume method to solve the considered equation. Firstly, we use the mid-point quadrature rule to transform the space distributed-order term into a multi-term fractional term, and the multi-term fractional equation is discretized by the finite volume. For the time-fractional derivative, we apply the finite difference method. We prove that the iterative scheme is uncondi-tionally stable and convergent. A numerical example is presented to verify the effectiveness of the proposed method.

1. 引言

随着科学技术的发展，分数阶微分方程因其在工程、物理、经济学、生物学、医学等诸多领域的应用而受到广泛关注。分数阶微分方程可用来模拟复杂扩散现象，包括反常扩散过程和大范围空间相互作用的扩散过程，而这些过程不能通过传统的二阶扩散方程精确模拟。许多研究表明，与单项、多项时间、空间或时空分数阶微分方程相比，分布阶微分方程更适合用来描述在整个时域内缺乏幂律尺度的加速慢扩散和减速超扩散。在Caputo [1] 首次提出用于非弹性介质应力–应变关系的分布阶微分方程之后，越来越多专家关注如何求解分布阶微分方程。

胡和刘 [2] 研究了一种隐式差分方法去求解一个新的时间分布阶和双侧空间分数阶对流–色散方程。卜等人 [3] 使用加权和移位Grünwald差分方法来研究时间分布阶扩散方程。Yamamoto [4] 研究了一类具有Caputo分数阶导数的时间分布阶扩散方程初边值问题，给出了解析解，并讨论了它在反问题中的应用。刘等人 [5] 提出了一种基于径向基函数的无网格方法来求解时间分布阶对流–扩散方程。也有一些关于研究空间分布阶微分方程的文章。如李和刘 [6] 提出了一种有限体积法求解空间分布阶对流–扩散方程。Jia和Wang [7] 发展了一种快速的有限差分法来求解一般凸域上的空间分布阶偏微分方程。

相比较而言，时间分数阶导数被认为是模拟慢扩散过程的有力工具 [8] ，而空间分布阶方程能灵活地表示介质的作用及其与流体的空间相互作用。因此，本文提出了一种有限体积法求解空间分布阶时间分数阶扩散方程：

$\frac{{\partial }^{\beta }u\left(x,t\right)}{\partial t^{\beta }}={\displaystyle {\int }_1^{2}P\left(\alpha \right)}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}\text{d}\alpha +f\left(x,t\right),\left(x,t\right)\in \left(0,L\right)\times \left(0,T\right)$ , (1)

边界值： $u\left(0,t\right)=0,u\left(L,t\right)=0,t\in \left(0,T\right)$ , (2)

初始值： $u\left(x,0\right)=\psi \left(x\right),x\in \left(0,L\right)$ , (3)

$f\left(x,t\right)$ 是源项， $P\left(\alpha \right)$ 是满足下面条件的非负权函数：

$P\left(\alpha \right)\ge 0$ , $P\left(\alpha \right)\cancel{\equiv }0$ , $\alpha \in \left(1,2\right)$ , $0<{\displaystyle {\int }_1^{2}P\left(\alpha \right)\text{d}\alpha }<\infty$ .

在通常情况下 [9] ， $P\left(\alpha \right)$ 可以表示为 $P\left(\alpha \right)=l^{\alpha -2}K\left[A_1\delta \left(\alpha -{\alpha }_1\right)+A_2\delta \left(\alpha -{\alpha }_2\right)\right]$ 。

l和K是正的常数， $A_1>0$ ， $A_2>0$ ， $0<{\alpha }_1<{\alpha }_2\le 2$ 。 $\frac{{\partial }^{\beta }u\left(x,t\right)}{\partial t^{\beta }}$ 定义为Caputo分数阶导数：

$\frac{{\partial }^{\beta }u\left(x,t\right)}{\partial t^{\beta }}=\frac{1}{\Gamma \left(1-\beta \right)}{\displaystyle {\int }_0^t\frac{\partial u\left(x,\xi \right)}{\partial \xi }}\frac{\text{d}\xi }{\left(t-\xi \right)}，0<\beta <1$ .

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}$ 是Riesz分数阶导数，定义为：

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}=-\frac{1}{2\cos \left(\alpha \text{π}/2\right)}(\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial x^{\alpha }}+\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left(-x\right)}^{\alpha }}),1<\alpha \le 2$ .

其中 $$ ， $\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left(-x\right)}^{\alpha }}=\frac{1}{\Gamma \left(2-\alpha \right)}\frac{{\partial }^2}{\partial x^2}{\displaystyle {\int }_x^{\infty }\frac{u\left(s,t\right)}{{\left(s-x\right)}^{\alpha -1}}\text{d}s}$ 。

本文提出了一种有限体积法求解空间分布阶时间分数阶扩散方程(1)~(3)。本文的结构如下：第二节，我们应用有限差分近似时间分数阶导数，然后离散空间分布阶方程为一个多项分数阶方程，进而我们使用有限体积法去求解多项分数阶方程，推导出Crank-Nicolson迭代格式。第三节，我们证明了迭代格式的稳定性和收敛性。第四节给出一个数值例子来说明本文数值方法的有效性。

2. 有限体积法求解空间分布阶时间分数阶扩散方程的Crank-Nicolson迭代格式

令 $t_n=n\tau$ ， $t=0,\text{}1,2,\cdots ,N$ ， $x_i=ih$ ， $i=0,1,2,\cdots ,M$ 离散区间 $\left[0,L\right]\times \left[0,T\right]$ ，其中 $h=L/M$ ， $\tau =T/N$ 。

采用文献 [10] 的离散方法，时间分数阶导数项可以近似为下面格式：

$\begin{array}{c}\frac{{\partial }^{\beta }u\left(x,t_n\right)}{\partial t^{\beta }}=\frac{1}{\Gamma \left(1-\beta \right)}{\displaystyle \underset{j=0}{\overset{n-1}{\sum }}{\displaystyle {\int }_{j\tau }^{\left(j+1\right)\tau }\frac{\partial u\left(x,\xi \right)}{\partial \xi }}}\frac{\text{d}\xi }{{\left(t_n-\xi \right)}^{\beta }}\\ =\frac{1}{\Gamma \left(1-\beta \right)}{\displaystyle \underset{j=0}{\overset{n-1}{\sum }}\frac{u\left(x,t_{j+1}\right)-u\left(x,t_j\right)}{\tau }{\displaystyle {\int }_{j\tau }^{\left(j+1\right)\tau }\frac{\text{d}\xi }{{\left(t_n-\xi \right)}^{\beta }}}}+o\left({\tau }^{2-\beta }\right)\\ =\frac{{\tau }^{1-\beta }}{\Gamma \left(2-\beta \right)}{\displaystyle \underset{j=0}{\overset{n-1}{\sum }}\frac{u\left(x,t_{n-j}\right)-u\left(x,t_{n-j-1}\right)}{\tau }\left[{\left(j+1\right)}^{1-\beta }-j^{1-\beta }\right]}+o\left({\tau }^{2-\beta }\right)\\ =a\left[u\left(x,t_n\right)-u\left(x,t_{n-1}\right)\right]+a{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}\left[u\left(x,t_{n-j}\right)-u\left(x,t_{n-j-1}\right)\right]b_j}+o\left({\tau }^{2-\beta }\right)\end{array}$ (4)

其中 $a=\frac{{\tau }^{-\beta }}{\Gamma \left(2-\beta \right)}$ ， $b_j={\left(j+1\right)}^{1-\beta }-j^{1-\beta }$ ， $j=0,1,2,\cdots ,n$ 。在(1)中，如果 $P\left(1\right)=0$ ，我们可以得到

${\displaystyle {\int }_1^{2}P\left(\alpha \right)}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}\text{d}\alpha ={\displaystyle {\int }_{1^{+}}^{2}P\left(\alpha \right)}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}\text{d}\alpha$ . (5)

定义 $1^{+}=1+\epsilon$ ， $\epsilon \to 0$ ，用网格 $1+\epsilon ={\xi }_0<{\xi }_1<\cdots <{\xi }_S=2$ 离散区间 $\left[1^{+},2\right]$ 且 $\Delta {\xi }_k={\xi }_k-{\xi }_{k-1}=\frac{1}{S}=\sigma$ ， $k=1,2,\cdots ,S$ ， ${\alpha }_k=\frac{{\xi }_k+{\xi }_{k+1}}{2}=1+\frac{2k-1}{2S}$ 。用中点求积法，得到

$\begin{array}{l}{\displaystyle {\int }_{1^{+}}^{2}P\left(\alpha \right)}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}\text{d}\alpha \\ ={\displaystyle \underset{k=1}{\overset{S}{\sum }}P\left({\alpha }_k\right)\frac{{\partial }^{{\alpha }_k}u\left(x,t\right)}{\partial {\left|x\right|}^{{\alpha }_k}}}\Delta {\xi }_k+o\left({\sigma }^2\right)\end{array}$ . (6)

结合(1)，(4)，(6)我们可以得到

$\begin{array}{l}a\left[u\left(x,t_n\right)-u\left(x,t_{n-1}\right)\right]+a{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{b}_j}\left[u\left(x,t_{n-j}\right)-u\left(x,t_{n-j-1}\right)\right]\\ =\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}P\left({\alpha }_k\right)}\frac{{\partial }^{{\alpha }_k}u\left(x,t_n\right)}{\partial {\left|x\right|}^{{\alpha }_k}}+f\left(x,t_n\right)+o\left({\tau }^{2-\beta }+{\sigma }^2\right)\\ =\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}{a}_{k}P\left({\alpha }_k\right)}\frac{\partial }{\partial x}[\frac{{\partial }^{{\gamma }_k}u\left(x,t_n\right)}{\partial x^{{\gamma }_k}}-\frac{{\partial }^{{\gamma }_k}u\left(x,t_n\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}]+f\left(x,t_n\right)+o\left({\tau }^{2-\beta }+{\sigma }^2\right)\end{array}$ (7)

其中 $1<{\alpha }_k<2$ ， $a_k=-\frac{1}{2\cos \left(\text{π}{\alpha }_k/2\right)}>0$ ， $0<{\gamma }_k={\alpha }_k-1<1$ ， $a_k=-\frac{1}{2\cos \left(\text{π}\left(1+{\gamma }_k\right)/2\right)}>0$ 。

定义 $x_{i-1/2}=\frac{x_i+x_{i-1}}{2}$ ， $$ 为区间 $\left[x_{i-1},x_i\right]$ 的中点。方程(7)在 $\left[x_{i-1/2},x_{i+1/2}\right]$ ， $i=1,2,\cdots ,M-1$ 下积分，得到

$\begin{array}{l}a{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}u\left(x,t_n\right)\text{d}x}-\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}{a}_{k}P\left({\alpha }_k\right)}{[\frac{{\partial }^{{\gamma }_k}u\left(x,t_n\right)}{\partial x^{{\gamma }_k}}-\frac{{\partial }^{{\gamma }_k}u\left(x,t_n\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}]}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\\ =a{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}u\left(x,t_{n-1}\right)\text{d}x}-a{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{b}_j}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left[u\left(x,t_{n-j}\right)-u\left(x,t_{n-j-1}\right)\right]\text{d}x}+{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f\left(x,t_n\right)\text{d}x}+o\left(h\left({\tau }^{2-\beta }+\sigma \right)\right)\end{array}$ (8)

现在，我们将空间 $V_h$ 定义为网格 ${\left\{\left[x_{i-1},x_i\right]\right\}}_{i=1,2,\cdots ,M}$ 上的分段线性多项式的集合。然后，近似解 $u_h\left(x,t_n\right)\in P\left(0,1\right)$ 可以用分段多项式表示为 $u_h\left(x,t_n\right)\approx {\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{u}_m^n}{\varphi }_m\left(x\right)$ 。

${\varphi }_i,0\le i\le M$ 是 $V_h$ 的节点基函数。 ${\varphi }_0,{\varphi }_1,\cdots ,{\varphi }_M$ 表示为下面形式：

${\varphi }_i\left(x\right)=\{\begin{array}{l}\frac{x-x_{i-1}}{h},x\in \left[x_{i-1},x_i\right]\\ \frac{x_{i+1}-x}{h},x\in \left[x_i,x_{i+1}\right],i=1,2,\cdots ,M-1\\ 0,其它\end{array}$ 和 ${\varphi }_0\left(x\right)=\{\begin{array}{l}\frac{x_1-x}{h},x\in \left[x_0,x_1\right]\\ 0,其它\end{array}$ ， ${\varphi }_M\left(x\right)=\{\begin{array}{l}\frac{x-x_{M-1}}{h},x\in \left[x_{M-1},x_M\right]\\ 0,其它\end{array}$

因此我们可以得到Crank-Nicolson格式：

$\begin{array}{l}a{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{u}_m^n}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\varphi }_m}\left(x\right)\text{d}x-\sigma {\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{u}_m^n}{\displaystyle \underset{k=1}{\overset{S}{\sum }}{a}_{k}P\left({\alpha }_k\right)}{[\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x\right)}{\partial x^{{\gamma }_k}}-\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}]}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\\ =a{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{u}_m^{n-1}}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\varphi }_m\left(x\right)\text{d}x}-a{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{b}_j}{u}_m^{n-j}}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\varphi }_m\left(x\right)\text{d}x}+a{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{b}_j}{u}_m^{n-j-1}}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\varphi }_m\left(x\right)\text{d}x}+{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f\left(x,t_n\right)\text{d}x}\end{array}$ (9)

通过直接计算，当 $1\le i$ ， $m\le M-1$ ，得到

${\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\varphi }_m\left(x\right)\text{d}x}=\{\begin{array}{l}h/8,\left|i-m\right|=1\\ 3h/4,i=m\\ 0,其它\end{array}$ 和

$\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i+1/2}\right)}{\partial x^{{\gamma }_k}}=\frac{1}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\{\begin{array}{l}0,m>i+1\\ 2^{{\gamma }_k-1},m=i+1\\ {\left(3/2\right)}^{1-{\gamma }_k}-2^{{\gamma }_k},m=i\\ c_{i-m+1}^k,m<i\end{array}$

$\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i-1/2}\right)}{\partial x^{{\gamma }_k}}=\frac{1}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\{\begin{array}{l}0,m>i\\ 2^{{\gamma }_k-1},m=i\\ {\left(3/2\right)}^{1-{\gamma }_k}-2^{{\gamma }_k},m=i-1\\ c_{i-m}^k,m<i-1\end{array}$

$\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i+1/2}\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}=\frac{1}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\{\begin{array}{l}{c}_{m-i}^k,m>i+1\\ {\left(3/2\right)}^{1-{\gamma }_k}-2^{{\gamma }_k},m=i-1\\ 2^{{\gamma }_k-1},m=i\\ 0,m<i\end{array}$

$\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i-1/2}\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}=\frac{1}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\{\begin{array}{l}{c}_{m-i+1}^k,\text{}m>i\\ {\left(3/2\right)}^{1-{\gamma }_k}-2^{{\gamma }_k},m=i\\ 2^{{\gamma }_k-1},m=i-1\\ 0,m<i-1\end{array}$

其中， $c_i^k={\left(i-3/2\right)}^{1-{\gamma }_k}-2{\left(i-1/2\right)}^{1-{\gamma }_k}+{\left(i+1/2\right)}^{1-{\gamma }_k},i=2,3,\cdots$ 。

因此，方程(9)可以表示为 $$ 令 $G_{im}=G_{1,im}-G_{2,im}$ ，其中

$G_{1,im}=\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}{a}_{k}P\left({\alpha }_k\right)}[\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i+1/2}\right)}{\partial x^{{\gamma }_k}}-\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i-1/2}\right)}{\partial x^{{\gamma }_k}}]$ , $G_{2,im}=\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}{a}_{k}P\left({\alpha }_k\right)}[\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i+1/2}\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}-\frac{{\partial }^{{\gamma }_k}{\varphi }_m\left(x_{i-1/2}\right)}{\partial {\left(-x\right)}^{{\gamma }_k}}]$ .

因此，当 $n=1$ 时，

$a\frac{h}{8}\left(u_{i-1}^1+6u_i^1+u_{i+1}^1\right)-{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{u}_m^1{G}_{im}}=a\frac{h}{8}\left(u_{i-1}^0+6u_i^0+u_{i+1}^0\right)+{\left(F^1\right)}_i$ , (10)

当 $n>1$ 时，

(11)

其中 $i=1,2,\cdots ,M-1$ ，令 ${\left(F^n\right)}_i={\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f\left(x,t_n\right)\text{d}x}$ ， $A_{im}={\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\varphi }_m\left(x\right)\text{d}x}$ ， $r=\frac{1}{a}={\tau }^{\beta }\Gamma \left(2-\beta \right)$ ， $d_j=b_{j-1}-b_j$ ， $j=1,2,\cdots ,n-1$ ， $U^n={\left(u_1^n,u_2^n,\cdots ,u_{M-1}^n\right)}^{\text{T}}$ 。因此方程(10)，(11)可以表示为

$\{\begin{array}{l}\left(A-rG\right)U^1=A{U}^0+rF{}^1,n=1\\ \left(A-rG\right)U^n=A\left(d_1{U}^{n-1}+d_2{U}^{n-2}+\cdots +d_{n-1}{U}^1+b_{n-1}{U}^0\right)+rF{}^n,n>1\end{array}$ (12)

边界值和初始值分别离散为 ${\psi }_i^0=\psi \left(ih\right)$ ，其中 $i=0,1,2,\cdots ,M$ ， $U^0={\left(u_1^0,u_2^0,\cdots ,u_{M-1}^0\right)}^{\text{T}}$ 。

3. 迭代方法的稳定性和收敛性

定理1. 当 $0<{\gamma }_k<1$ 时，系数 $G_{im}$ 满足：

$\left|G_{ii}\right|>{\displaystyle \underset{m=1,m\ne i}{\overset{M-1}{\sum }}\left|G_{im}\right|},i=1,2,\cdots ,M-1$ ，则G是严格对角占优的。

定理2. 当 $0<{\gamma }_k<1$ 时， $B=A-rG$ ，B是严格对角占优的，并且 $B^{-1}$ 的谱半径满足： $\rho \left(B^{-1}\right)<\frac{2}{h}$ 。

证明：

$B_{im}=\{\begin{array}{l}r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left(c_{m-i}^k-c_{m-i+1}^k\right),m>i+1}\\ \frac{h}{8}-r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left[3{\left(\frac{1}{2}\right)}^{1-{\gamma }_k}-{\left(\frac{3}{2}\right)}^{1-{\gamma }_k}+c_2^k\right]，m=i+1}\\ \frac{3h}{4}-r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}2\left[{\left(\frac{3}{2}\right)}^{1-{\gamma }_k}-3{\left(\frac{1}{2}\right)}^{1-{\gamma }_k}\right],m=i}\\ \frac{h}{8}-r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left[3{\left(\frac{1}{2}\right)}^{1-{\gamma }_k}-{\left(\frac{3}{2}\right)}^{1-{\gamma }_k}+c_2^k\right],m=i-1}\\ r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left(c_{i-m}^k-c_{i-m+1}^k\right),m<i-1}\end{array}$

类似与 [6] 定理2的证明，可得该定理的证明。

定理 3. 令 $D=\left(\lambda -1\right)A-\lambda rG$ ，当时，如果 $\lambda >1$ 或 $\lambda <-1$ ，则D是严格对角占优的，且

$\left|D_{ii}\right|>{\displaystyle \underset{m=1,m\ne i}{\overset{M-1}{\sum }}\left|D_{im}\right|},i=1,2,\cdots ,M-1$ .

证明：

$D_{im}=\{\begin{array}{l}\lambda r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left(c_{m-i}^k-c_{m-i+1}^k\right),m>i+1}\\ \left(\lambda -1\right)\frac{h}{8}-\lambda r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left[3{\left(\frac{1}{2}\right)}^{1-{\gamma }_k}-{\left(\frac{3}{2}\right)}^{1-{\gamma }_k}+c_2^k\right],m=i+1}\\ \left(\lambda -1\right)\frac{3h}{4}-\lambda r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}2\left[{\left(\frac{3}{2}\right)}^{1-{\gamma }_k}-3{\left(\frac{1}{2}\right)}^{1-{\gamma }_k}\right],m=i}\\ \left(\lambda -1\right)\frac{h}{8}-\lambda r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left[3{\left(\frac{1}{2}\right)}^{1-{\gamma }_k}-{\left(\frac{3}{2}\right)}^{1-{\gamma }_k}+c_2^k\right],m=i-1}\\ \lambda r\sigma {\displaystyle \underset{k=1}{\overset{S}{\sum }}\frac{a_{k}P\left({\alpha }_k\right)}{\Gamma \left(2-{\gamma }_k\right)h^{{\gamma }_k}}\left(c_{i-m}^k-c_{i-m+1}^k\right),m<i-1}\end{array}$

类似与 [11] 定理3的证明，可得该定理的证明。

定理4. 矩阵 $B^{-1}A$ 的谱半径满足 $\rho \left(B^{-1}A\right)<1$ 。

证明：类似与 [11] 定理4的证明，可得该定理的证明。

定理5. 假定 ${\tilde{u}}_i^n$ 和 $u_i^n$ $\left(i=1,2,\cdots ,M-1;n=0,1,\cdots ,N\right)$ 分别为差分格式(12)的近似解和数值解且 ${\epsilon }_i^n={\tilde{u}}_i^n-u_i^n$ ， $Y^n={\left({\epsilon }_1^n,{\epsilon }_2^n,\cdots ,{\epsilon }_{M-1}^n\right)}^{\text{T}}$ ，则 $\Vert Y^n\Vert \le \Vert Y^0\Vert$ ，因此差分格式(12)是无条件稳定的。

证明：现在我们用数学归纳法来证明稳定性。

当 $n=1$ 时，将 ${\epsilon }_i^n={\tilde{u}}_i^n-u_i^n$ 代入(12)，得到 $Y^1={\left(A-rG\right)}^{-1}A{Y}^0$ 。令 $Q={\left(A-rG\right)}^{-1}A$ ，根据定理4，我们知道 $\rho \left(Q\right)<1$ ，则存在向量范数和诱导矩阵范数 $\Vert .\Vert$ ，使 $\Vert Q\Vert <1$ 。使用以上范数，得到：

.

假设当 $n\le k-1$ 时，成立 $\Vert Y^{k-1}\Vert \le \Vert Y^0\Vert$ 。

则当 $n=k$ 时，由文献 [12] 可知 ${\displaystyle \underset{j=1}{\overset{k}{\sum }}{d}_j}=1-b_k,j=1,\cdots ,n$ 。因此

$\begin{array}{c}\Vert Y^k\Vert =\Vert {\left(A-rG\right)}^{-1}A\left(d_1{Y}^{k-1}+d_2{Y}^{k-2}+\cdots +d_{k-1}{Y}^1+b_{k-1}{Y}^0\right)\Vert \\ \le \Vert Q\Vert \left(d_1+d_2+\cdots +d_{k-1}+b_{k-1}\right)\Vert Y^0\Vert \\ \le \left(d_1+d_2+\cdots +d_{k-1}+b_{k-1}\right)\Vert Y^0\Vert \\ =\left(1-b_{k-1}+b_{k-1}\right)\Vert Y^0\Vert \\ =\Vert Y^0\Vert \end{array}$

得到 $\Vert Y^n\Vert \le \Vert Y^0\Vert$ ，定理证毕。

定理6. 假定 $u^n$ 为(1)~(3)的精确解，则当 $\sigma$ ， $\tau$ 和h趋近于零时，数值解 $U^n$ 无条件收敛到精确解 $u^n$ ，并且，

$\Vert u^n-U^n\Vert \le C\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)$ .

其中 $u^n={\left(u\left(x_1,t_n\right),u\left(x_2,t_n\right),\cdots ,u\left(x_{M-1},t_n\right)\right)}^{\text{T}}$ ， $U^n={\left(u_1^n,u_2^n,\cdots ,u_{M-1}^n\right)}^{\text{T}}$ .

证明：设为点 $\left(x_i,t_n\right)$ 处的误差。将 $e_i^n=u\left(x_i,t_n\right)-u_i^n$ 代入(10)和(11)，并结合(7)和文献 [13] 的引理3，推论1，得到：

$\{\begin{array}{l}\frac{h}{8}\left(e_{i-1}^1+6e_i^1+e_{i+1}^1\right)-r{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{e}_m^1{G}_{im}}=\frac{h}{8}\left(e_{i-1}^0+6e_i^0+e_{i+1}^0\right)+ro\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right),n=1\\ \frac{h}{8}\left(e_{i-1}^n+6e_i^n+e_{i+1}^n\right)-r{\displaystyle \underset{m=1}{\overset{M-1}{\sum }}{e}_m^n{G}_{im}}=\left(1-b_1\right)\frac{h}{8}\left(e_{i-1}^{n-1}+6e_i^{n-1}+e_{i+1}^{n-1}\right)\\ +{\displaystyle \underset{j=1}{\overset{n-2}{\sum }}\frac{h}{8}\left(e_{i-1}^{n-j-1}+6e_i^{n-j-1}+e_{i+1}^{n-j-1}\right)}\left(b_j-b_{j+1}\right)+\frac{h}{8}\left(e_{i-1}^0+6e_i^0+e_{i+1}^0\right)b_{n-1}+ro\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right),n>1\end{array}$

其中 $e_0^n=e_M^n=0$ ， $e_i^0=0,i=1,2,\cdots ,M-1$ ，且矩阵形式为：

$\{\begin{array}{l}\left(A-rG\right)E^1=A{E}^0+ro\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\chi ,n=1\\ \left(A-rG\right)E^n=A\left(d_1{E}^{n-1}+d_2{E}^{n-2}+\cdots +d_{n-1}{E}^1+b_{n-1}{E}^0\right)+ro\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\chi ,n>1\end{array}$

$\chi ={\left(1,1,\cdots ,1\right)}^{\text{T}},E^n={\left(e_1^n,e_2^n,\cdots ,e_{M-1}^n\right)}^{\text{T}},Q={\left(A-rG\right)}^{-1}A,B=A-rG$ .

下面，我们用数学归纳法来证明这个定理。

当 $n=1$ 时，得到

$E^1=Q{E}^0+r{\left(A-rG\right)}^{-1}o\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\chi =r{B}^{-1}o\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\chi$ .

根据定理2和定理4已知 $\rho \left(B^{-1}\right)<\frac{2}{h}$ ， $\rho \left(Q\right)<1$ ，则存在向量范数和诱导矩阵范数 $\Vert .\Vert$ ，使得 $\Vert Q\Vert <1$ ， $\Vert B^{-1}\Vert <c{h}^{-1}$ ，且 $r={\tau }^{\beta }\Gamma \left(2-\beta \right)$ ， $\tau =\frac{T}{n}$ 。根据以上范数，得到

$\Vert E^1\Vert =\Vert r{B}^{-1}o\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\chi \Vert \le o\left({\tau }^{\beta }\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)$ .

假设当 $n\le k-1$ ，成立

$\Vert E^{k-1}\Vert \le o\left({\tau }^{\beta }\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)$ .

则，当 $n=k$ 时，有

$\begin{array}{c}\Vert E^k\Vert =\Vert Q\left(d_1{E}^{k-1}+d_2{E}^{k-2}+\cdots +d_{k-1}{E}^1+b_{k-1}{E}^0\right)+r{B}^{-1}o\left(h\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\chi \Vert \\ \le \Vert Q\Vert \left(d_1+d_2+\cdots +d_{k-1}+b_{k-1}\right)o\left({\tau }^{\beta }\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)+o\left({\tau }^{\beta }\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\\ \le \left(1-b{}_{k-1}+b_{k-1}\right)o\left({\tau }^{\beta }\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\\ =o\left({\tau }^{\beta }\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)\right)\end{array}$

因此， $\Vert E^n\Vert \le C\left({\tau }^{2-\beta }+{\sigma }^2+h^2\right)$ ，定理证毕。

4. 数值例子

考虑下面的空间分布阶时间分数阶方程：

$\frac{{\partial }^{\beta }u\left(x,t\right)}{\partial t^{\beta }}={\displaystyle {\int }_1^{2}P\left(\alpha \right)}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}\text{d}\alpha +f\left(x,t\right),\left(x,t\right)\in \left(0,1\right)\times \left(0,1\right]$ ,

边界值： $u\left(0,t\right)=0,u\left(1,t\right)=0$ ,

初始值： $u\left(x,0\right)=x^2{\left(1-x\right)}^2$ .

且 $P\left(\alpha \right)=-2\Gamma \left(5-\alpha \right)\cos \left(\frac{\text{π}\alpha }{2}\right)$ ,

$\begin{array}{c}f\left(x,t\right)=\frac{2}{\Gamma \left(3-\beta \right)}{t}^{2-\beta }{x}^2{\left(1-x\right)}^2-{\displaystyle {\int }_1^{2}P\left(\alpha \right)}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial {\left|x\right|}^{\alpha }}\text{d}\alpha \\ =\frac{2}{\Gamma \left(3-\beta \right)}{t}^{2-\beta }{x}^2{\left(1-x\right)}^2-\left(t^2+1\right)\left[R\left(x\right)+R\left(1-x\right)\right]\end{array}$ ,

其中

$\begin{array}{c}R\left(x\right)=\Gamma \left(5\right)\frac{1}{\ln x}\left(x^3-x^2\right)-2\Gamma \left(4\right)\left[\frac{1}{\ln x}\left(3x^2-2x\right)-\frac{1}{{\left(\ln x\right)}^2}\left(x^2-x\right)\right]\\ +\Gamma \left(3\right)\frac{1}{\ln x}\left[6x-2-\frac{5x}{\ln x}+\frac{3}{\ln x}+\frac{2x}{{\left(\ln x\right)}^2}-\frac{2}{{\left(\ln x\right)}^2}\right]\end{array}$

上述方程的精确解为 $u\left(x,t\right)=\left(t^2+1\right)x^2{\left(1-x\right)}^2$ 。

图1显示了这个例子的精确解和数值解的对比，可以看出，数值解与精确解非常吻合。表1显示了我们的数值方法在 $\tau ,h\left(\tau =h=1/1000\right)$ 时关于 $\sigma$ 的误差和收敛性。随着 $\sigma$ 的减少， $\sigma$ 的收敛阶达到二阶。表2表示了在 $\sigma ,\tau \left(\sigma =\tau =1/500\right)$ 时关于h的误差和收敛性，随着h的减少，h的收敛阶也达到二阶。表3表示了在 $\sigma =h=1/400$ 时 $\tau$ 的收敛阶，可以看出当 $\tau$ 减少时， $\tau$ 的收敛阶是 $2-\beta$ 。从下面图表的误差和收敛阶可以看出，我们的数值方法是有效的，稳定的。

Figure 1. Exact solution and numerical solution with $\sigma =\tau =h=1/100$ at $t=1.0$ and $\beta =0.7$

图1. 在 $\sigma =\tau =h=1/100,t=1.0,\beta =0.7$ 时的数值解和精确解

Table 1. The error and the convergence order of σ for τ = h = 1 / 1000 at t = 1 and β = 0.7

表1. 在 $\tau =h=1/1000,\beta =0.7,t=1$ 时 $\sigma$ 的误差和收敛阶

Table 2. The error and the convergence order of h for τ = σ = 1 / 500 at and β = 0.7

表2. 在 $\tau =\sigma =1/500,\beta =0.7,t=1$ 时h的误差和收敛阶

Table 3. The error and the convergence order of τ for h = σ = 1 / 400 at t = 1 and β = 0.7

表3. 在 $h=\sigma =1/400,\beta =0.7,t=1$ 时 $\tau$ 的误差和收敛阶

5. 结论

本文讨论了空间分布阶时间分数阶扩散方程的数值解。首先，基于有限体积法，得到了该方程的离散格式。其次，证明了离散格式的无条件稳定性和收敛性。最后，通过一个算例验证了该方法的有效性。在未来，我们希望能提出一种新的数值方法来求解二维的空间分布阶时间分数阶扩散方程，并提高其收敛速度。

基金项目

本研究由湖南省自然科学基金资助(批准号：No.2018JJ3519)和湖南省教育厅科研项目(批准号：NO.17B003)。

NOTES

^*通讯作者。