时间分数阶Fokker-Planck方程有限体积法

A Finite Volume Method for Time Fractional Fokker-Planck Equations

摘要: 我们设计一种有限体积法求解变外力场下的时间分数阶Fokker-Planck方程，其中外力场与时间空间相关。实验表明该方法在时间上和空间上分别具有一阶和二阶收敛性。

Abstract: We present a finite volume method for solving the time fractional Fokker-Planck equations with space-and-time-dependent forcing. Numerical test shows that the convergence rates for time and for space are order 1 and order 2 respectively.

1. 研究的问题

考虑如下的时间分数阶Fokker-Planck方程(FFPE)：

$\frac{\partial u}{\partial t}=\left(k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\right)D_t^{1-\alpha }u,a<x<b,0<t\le T,$ (1.1)

初始条件和边值条件为

$u\left(x,0\right)=\phi \left(x\right),a<x<b,u\left(a,t\right)=g_1\left(t\right),u\left(b,t\right)=g_2\left(t\right),0<t\le T,$ (1.2)

其中 $\alpha \in \left(0,1\right)$ ， $k_{\alpha }$ 为正常数， $f\left(x,t\right),\phi \left(x\right),g_1\left(t\right),g_2\left(t\right)$ 为给定函数，方程(1.1)中的分数阶导数为Riemann-Liouville分数阶导数 $D_t^{1-\alpha }u\left(x,t\right)=\frac{1}{\Gamma \left(\alpha \right)}\frac{\partial }{\partial t}{\displaystyle {\int }_0^t\frac{u\left(x,s\right)}{{\left(t-s\right)}^{1-\alpha }}\text{d}s}$ ， $\Gamma \left(x\right)$ 是Gamma函数。方程(1.1)可用于模拟外力场作用下的低扩散现象(参见 [1] )。

对于方程(1.1)的求解，已有数值解法多数是针对f是常数或是x的函数情况(参见 [2] - [10] )：如对情形，Chen [5] 设计有限差分方法，Jiang [3] 针对Chen [5] 中的数值格式给出了稳定性和收敛性的证明；Vong [10] 设计了一种高阶差分格式并证明它的稳定性和收敛性。据我们所知，对于 $f=f\left(x,t\right)$ 情况，仅有Le [11] 研究了两种方法，第一种是时间上连续方法(在空间中使用分段线性Galerkin有限元方法)，第二种是空间连续方法。

有限体积法已开始应用于求解分数阶方程，如 [12] 对空间分数阶方程采用了有限体积方法； [13] 利用有限体积方法数值求解时间–空间分数阶方程； [14] 对二维时间分数阶偏微分方程进行有限体积方法研究(其中 $f=0$ )。

本文设计求解(1.1)的有限体积法，其中空间导数使用中心差分格式进行离散，时间分数阶导数使用一阶线性多步法进行离散。数值实验结果表明该方法在时间上和空间上分别具有一阶和二阶收敛性。文中假定解u充分光滑，C网格大小无关的正常数。

2. 离散

在 $\left[0,T\right]$ 上取分点为 $t_k=k\Delta t$ ， $k=0,1,\cdots ,L$ ，其中时间步长 $\Delta t=T/L$ ，L为正整数。在区间 $\left[a,b\right]$ 上取分点 $x_i=a+ih,i=0,1,\cdots ,N+1$ ，其中空间步长 $h=\left(b-a\right)/\left(N+1\right)$ ，N为正整数。N个空间有限体为 $\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right],i=1,\cdots ,N$ ，其中 $x_{i+\frac{1}{2}}=\frac{x_i+x_{i+1}}{2},i=0,1,\cdots ,N$ 。

在方程(1.1)中取 $t=t_n\left(n=0,1,\cdots ,L\right)$ 得到

${\frac{\partial u}{\partial t}|}_{t_n}=\left[k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t_n\right)\right]D_t^{1-\alpha }u\left(x,t_n\right).$ (2.1)

利用向后差分，(2.1)式左侧可以写成

${\frac{\partial u}{\partial t}|}_{t_n}=\frac{u\left(x,t_n\right)-u\left(x,t_{n-1}\right)}{\Delta t}+r_n^{\left(1\right)},$ (2.2)

其中 $r_n^{\left(1\right)}={\frac{\partial u}{\partial t}|}_{t_n}-\frac{u\left(x,t_n\right)-u\left(x,t_{n-1}\right)}{\Delta t}$ 满足 $\left|r_n^{\left(1\right)}\right|\le C\Delta t$ 。

用Lubich一阶线性多步法逼近Riemann-Liouville分数阶导数(参考文献 [15] )

$D_t^{1-\alpha }u\left(x,t_n\right)=\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }u\left(x,t_j\right)}+r_n^{\left(2\right)},$ (2.3)

其中 ${\omega }_k^{\alpha }=\frac{{\left(-1\right)}^k\Gamma \left(\alpha +1\right)}{\Gamma \left(k+1\right)\Gamma \left(\alpha -k+1\right)}\left(k=0,1,2,\cdots \right)$ ， $r_n^{\left(2\right)}=D_t^{1-\alpha }u\left(x,t_n\right)-\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }u\left(x,t_j\right)}$ 满足 $r_n^{\left(2\right)}=C\left(x\right)\Delta t$ ， $C\left(x\right)$ 是与解u有关的光滑函数。

将(2.2)式，(2.3)式带入(2.1)式中，得

$\frac{u\left(x,t_n\right)-u\left(x,t_{n-1}\right)}{\Delta t}+r_n^{\left(1\right)}=\left[k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t_n\right)\right]\left[\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }u\left(x,t_j\right)}+r_n^{\left(2\right)}\right].$ (2.4)

用 $u_n\left(x\right)$ 表示 $u\left(x,t_n\right)$ 的近似，由(2.4)我们得到原问题的时间半离散格式

$\frac{u_n\left(x\right)-u_{n-1}\left(x\right)}{\Delta t}=\left[k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t_n\right)\right]\left[\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }{u}_j\left(x\right)}\right],$ (2.5)

记

$v_n\left(x\right):=\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }{u}_j\left(x\right)}.$ (2.6)

易知

$u_n\left(x\right):=\Delta t^{1-\alpha }{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{\alpha -1}{v}_j\left(x\right)}.$ (2.7)

将(2.6)式，(2.7)式带入(2.5)式中，得到

$\Delta t^{-\alpha }\left[{\displaystyle {\sum }_{j=0}^{n-1}\left({\omega }_{n-j}^{\alpha -1}-{\omega }_{n-j-1}^{\alpha -1}\right)v_j\left(x\right)}+{\omega }_0^{\alpha -1}{v}_n\left(x\right)\right]=\left[k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t_n\right)\right]v_n\left(x\right).$ (2.8)

在有限体 $\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ 上对(2.8)式积分

${\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\Delta t^{-\alpha }\left[{\displaystyle {\sum }_{j=0}^{n-1}\left({\omega }_{n-j}^{\alpha -1}-{\omega }_{n-j-1}^{\alpha -1}\right)v_j\left(x\right)}+{\omega }_0^{\alpha -1}{v}_n\left(x\right)\right]\text{d}x}={\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left[k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t_n\right)\right]v_n\left(x\right)\text{d}x}.$ (2.9)

利用中矩形积分公式，用 $v_i^n$ 表示 $v_n\left(x_i\right)$ ，(2.9)式左侧可以写成

$\begin{array}{l}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\Delta t^{-\alpha }\left[{\displaystyle {\sum }_{j=0}^{n-1}\left({\omega }_{n-j}^{\alpha -1}-{\omega }_{n-j-1}^{\alpha -1}\right)v_j\left(x\right)}+{\omega }_0^{\alpha -1}{v}_n\left(x\right)\right]\text{d}x}\\ =h\Delta t^{-\alpha }\left[{\displaystyle {\sum }_{j=0}^{n-1}\left({\omega }_{n-j}^{\alpha -1}-{\omega }_{n-j-1}^{\alpha -1}\right)v_i^j}+{\omega }_0^{\alpha -1}{v}_i^n\right]-h{r}_{i,n}^{\left(1\right)},\end{array}$ (2.10)

其中

$h\left|r_{i,n}^{\left(1\right)}\right|\le C{h}^3.$ (2.11)

用 $f_{i+1/2,n}$ 表示 $f\left(x_{i+1/2},t_n\right)$ ，(2.9)式右侧可以写成

${\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left[k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t_n\right)\right]v_n\left(x\right)\text{d}x}=k_{\alpha }\left({\frac{\partial v_n}{\partial x}|}_{x_{i+\frac{1}{2}}}-{\frac{\partial v_n}{\partial x}|}_{x_{i-\frac{1}{2}}}\right)-\left(f_{i+\frac{1}{2},n}{v}_{i+\frac{1}{2}}^n-f_{i-\frac{1}{2},n}{v}_{i-\frac{1}{2}}^n\right).$ (2.12)

用中心差分格式，(2.12)式右侧第一项可以写成

$k_{\alpha }\left({\frac{\partial v_n}{\partial x}|}_{x_{i+\frac{1}{2}}}-{\frac{\partial v_n}{\partial x}|}_{x_{i-\frac{1}{2}}}\right)=k_{\alpha }\left(\frac{v_{i+1}^n-v_i^n}{h}-\frac{v_i^n-v_{i-1}^n}{h}\right)+h{r}_{i,n}^{\left(2\right)},$ (2.13)

其中

$h\left|r_{i,n}^{\left(2\right)}\right|\le C{h}^3;$ (2.14)

由于(2.12)式右侧第二项可以写成

$\begin{array}{l}{f}_{i+\frac{1}{2},n}{v}_{i+\frac{1}{2}}^n-f_{i-\frac{1}{2},n}{v}_{i-\frac{1}{2}}^n\\ =f_{i+\frac{1}{2},n}\frac{v_i^n+v_{i+1}^n}{2}-f_{i-\frac{1}{2},n}\frac{v_i^n+v_{i-1}^n}{2}-h{r}_{i,n}^{\left(3\right)},\end{array}$ (2.15)

其中

$h\left|r_{i,n}^{\left(3\right)}\right|\le C{h}^3.$ (2.16)

将(2.10)，(2.12)，(2.13)，(2.15)式带入(2.9)式中，得 $i=1,2,\cdots ,N;n=1,2,\cdots ,L$ ，

$\begin{array}{l}h\Delta t^{-\alpha }\left[{\displaystyle {\sum }_{j=0}^{n-1}\left({\omega }_{n-j}^{\alpha -1}-{\omega }_{n-j-1}^{\alpha -1}\right)v_i^j}+{\omega }_0^{\alpha -1}{v}_i^n\right]-h{r}_{i,n}^{\left(1\right)}\\ =k_{\alpha }\left(\frac{v_{i+1}^n-v_i^n}{h}-\frac{v_i^n-v_{i-1}^n}{h}\right)+h{r}_{i,n}^{\left(2\right)}-\left(f_{i+\frac{1}{2},n}\frac{v_i^n+v_{i+1}^n}{2}-f_{i-\frac{1}{2},n}\frac{v_i^n+v_{i-1}^n}{2}\right)+h{r}_{i,n}^{\left(3\right)},\end{array}$ (2.17)

用 $V_i^n$ 近似 $v_i^n$ ，由(2.17)式我们可以得到如下的有限体积法(FV)： $i=1,2,\cdots ,N;n=1,2,\cdots ,L$ ，

$\begin{array}{l}h\Delta t^{-\alpha }\left[{\displaystyle {\sum }_{j=0}^{n-1}\left({\omega }_{n-j}^{\alpha -1}-{\omega }_{n-j-1}^{\alpha -1}\right)V_i^j}+{\omega }_0^{\alpha -1}{V}_i^n\right]\\ =k_{\alpha }\left(\frac{V_{i+1}^n-V_i^n}{h}-\frac{V_i^n-V_{i-1}^n}{h}\right)-\left(f_{i+\frac{1}{2},n}\frac{V_i^n+V_{i+1}^n}{2}-f_{i-\frac{1}{2},n}\frac{V_i^n+V_{i-1}^n}{2}\right),\end{array}$ (2.18)

初始条件和边值条件为 $\left(i=1,2,\cdots ,N;n=1,2,\cdots ,L\right)$

$\begin{array}{l}{V}_i^0=\phi \left(x_i\right),V_0^n=G_1\left(t_n\right)=\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }{g}_1\left(t_n\right)},\\ V_{N+1}^n=G_2\left(t_n\right)=\Delta t^{\alpha -1}{\displaystyle {\sum }_{j=0}^n{\omega }_{n-j}^{1-\alpha }{g}_2\left(t_n\right)}.\end{array}$ (2.19)

3. 数值实验

本节利用我们设计的有限体积法解决{(1.1)(1.2)}问题。考虑下列具有精确解的方程

$\frac{\partial u}{\partial t}=\left(\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x,t\right)\right)D_t^{1-\alpha }u+g\left(x,t\right)，0\le x\le 1,0\le t\le 1,$ (3.1)

初始条件和边值条件为 $u\left(x,0\right)=0,g_1\left(t\right)=t^2,g_2\left(t\right)=-t^2$ ，其中

$\begin{array}{c}g\left(x,t\right)=2t\cos \left(\text{π}x\right)+\frac{2{\text{π}}^2}{\Gamma \left(2+\alpha \right)}{t}^{1+\alpha }\cos \left(\text{π}x\right)\\ +\frac{2}{\Gamma \left(2+\alpha \right)}{t}^{1+\alpha }\left[\left(1-2x+t\right)\cos \left(\text{π}x\right)-f\left(x,t\right)\text{π}\sin \left(\text{π}x\right)\right],\end{array}$ (3.2)

Table 1. Convergence rate for space with α = 0.2, L = 50,000

表1. 空间收敛阶α = 0.2，L = 50,000

Table 2. Convergence rate for space with α = 0.5, L = 50,000

表2. 空间收敛阶α = 0.5，L = 50,000

Table 3. Convergence rate for space with α = 0.8, L = 50,000

表3. 空间收敛阶α = 0.8，L = 50,000

Table 4. Convergence rate for time with α = 0.2, N = 5000

表4. 空间收敛阶α = 0.2，N = 5000

Table 5. Convergence rate for time α = 0.5, N = 5000

表5. 空间收敛阶α = 0.5，N = 5000

Table 6. Convergence rate for time with α = 0.8, N = 5000

表6. 空间收敛阶α = 0.8，N = 5000

$f\left(x,t\right)=x-x^2+t+xt$ 。此方程的精确解为 $u\left(x,t\right)=t^2\cos \left(\text{π}x\right)$ 。

我们定义空间和时间收敛阶如下：

$空间收敛阶=\left|\frac{\ln \left({\Vert 细网格误差\Vert }_1/{\Vert 粗网格误差\Vert }_1\right)}{\ln \left(细网格划分数N+1/粗网格划分数N+1\right)}\right|,$

$时间收敛阶=\left|\frac{\ln \left({\Vert 细网格误差\Vert }_1/{\Vert 粗网格误差\Vert }_1\right)}{\ln \left(细网格划分数L/粗网格划分数L\right)}\right|.$

空间收敛阶的数值结果列于表1~表3中，时间收敛阶的数值结果列于表4~表6中。数值实验表明该方法在时间和空间上皆为一阶收敛，当空间网格足够细时，空间上可以达到二阶收敛。

致谢

感谢我的导师姜英军教授和我的师兄周帅虎对我论文写作的悉心指导和热情帮助。

基金项目

国家自然科学基金(11571053)。