Fokker-Planck方程的QUICK离散格式研究

QUICK Discrete Scheme for Fokker-Planck Equation

摘要: 我们研究了一种求解时间分数阶Fokker-Planck方程的有限体积法，时间上用L1近似，在空间对流项利用QUICK格式离散，扩散项用中心差分离散。数值实验结果表明该方法在空间上为二阶收敛。

Abstract: We design a finite volume method for solving time fractional Fokker-Planck equation. The time is dispersed by L1-approximate, the space convection term is discretized by QUICK scheme, and the diffusion term is discretized by central difference. The numerical results show that the method is second-order convergent in space.

1. 研究的问题

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

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

初始条件和边值条件为

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

其中 $\alpha \in \left(0,1\right)$， $k_{\alpha }$ 为正常数， $f,\phi \left(x\right),g_1\left(t\right),g_2\left(t\right)$ 为给定函数，方程(1.1)中的分数阶导数为Riemann-Liouville分数阶导数： $D_t^{1-\alpha }\omega \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] )， $k_{\alpha }$ 表示广义扩散系数，f表示外力场。方程(1.1)可改写为其等价形式：

$\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}=\left(k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\right)\omega \left(x,t\right),a\le x\le b,0\le t\le T.$ (1.3)

其中 ${\partial }^{\alpha }\omega /\partial t^{\alpha }$ 表示 $\alpha \left(0<\alpha <1\right)$ 阶Caputo分数阶导数： $\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}\frac{\partial }{\partial t}{\displaystyle {\int }_0^t\frac{\partial \omega \left(x,\eta \right)}{\partial t}\frac{\text{d}\eta }{{\left(t-\eta \right)}^{\alpha }}}$。

对于分数阶Fokker-Planck方程的求解，已有数值方法，大多是针对f是常数的函数情况(参见 [2] - [10] )。其中对 $f=0$ 情形，Jiang [3] 对Chen等 [5] 中的数值格式给出了稳定性和收敛性的证明；Vong和Wang [10] 开发了一种高阶差分格式求解时间分数阶Fokker-Planck方程，并获得了它的稳定性和收敛性。对于 $f=f\left(x\right)$ (变外力场)情况，已有的数值方法非常少，我们发现Le等 [11] 研究了两种方法，第一种是时间上连续(在空间中使用分段线性Galerkin有限元方法离散化)，第二种是空间连续(采用类似于经典隐式欧拉方法的时间步长方法)，并证明它的稳定性和收敛性。

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

我们研究了一种求解(1.3)的FV方法，在空间上利用对流项的二次向上差分格式离散，时间上利用 $L_1$ 近似离散。数值实验结果表明该方法在空间上为二阶收敛。

本文中，假定解 $\omega$ 充分光滑， $f\left(x\right)$ 满足如下Lipschitz条件，其中C表示正常数，与网格大小无关。

$\left|f\left(x\right)-f\left(x^{\prime }\right)\right|\le Cx-x^{\prime },x,x^{\prime }\in \left[a,b\right].$ (1.4)

2. 离散

设N，L为正整数，取空间步长 $h=\left(b-a\right)/\left(N+1\right)$，时间步长 $\Delta t=T/L$。区间 $\left[a,b\right]$ 上 $N+1$ 等分，分点为 $x_i=a+ih,i=0,1,\cdots ,N+1$，得到N个小区间 $\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right],i=0,1,\cdots ,N$ ；在区间 $\left[0,t\right]$ 上L等分，分点为 $t_k=k\Delta t,k=0,1,\cdots ,L$。为了方便，仅取f在点 $x_{i+\frac{1}{2}},i=0,1,\cdots ,N$ 处的值，记 $f_{x+\frac{1}{2}}=f\left(x_{i+\frac{1}{2}}\right)$。

在方程(1.3)中取 $t=t_n\left(n=0,1,\cdots ,L\right)$，在第i个控制体积上即区间 $\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ 上对方程两边求积得到

${\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}|}_{t_n}\text{d}x}=k_{\alpha }\left({\left(\frac{\partial \omega }{\partial x}\right)}_{i+\frac{1}{2},n}-{\left(\frac{\partial \omega }{\partial x}\right)}_{i-\frac{1}{2},n}\right)-\left({\left(f\omega \right)}_{i+\frac{1}{2},n}-{\left(f\omega \right)}_{i-\frac{1}{2},n}\right),i=1,2,\cdots ,N$ (2.1)

将 $\omega \left(x_i,t_n\right)$ 记做 ${\omega }_i^n\left(i=0,1,\cdots ,N+1;n=0,1,\cdots ,L\right)$，(2.1)式左边可以改写为：

$\begin{array}{l}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}|}_{t_n}\text{d}x}={\left(\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}\right)}_{\left(x_i,t_n\right)}h-h{\gamma }_{i,n}^{\left(1\right)}\\ =\frac{h\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left({\omega }_{i,n}-{\displaystyle {\sum }_{k=1}^{n-1}\left(a_{n-k-1}-a_{n-k}\right){\omega }_i^k}-a_{n-1}{\omega }_i^0\right)-h{\gamma }_{i,n}^{\left(1\right)}-h{\gamma }_{i,n}^{\left(2\right)},\end{array}$ (2.2)

(2.3)式中第一个等式应用中矩形公式，第二个等式应用了 $L_1$ 近似，其中 $a_k={\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }$，(见 [15] [16] )，其中误差项

${\gamma }_{i,n}^{\left(1\right)}:=\left({\left(\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}\right)}_{\left(x_i,t_n\right)}h-{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}|}_{t_n}}\right)/h$ (2.3)

${\gamma }_{i,n}^{\left(2\right)}:=\left(\frac{h\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left({\omega }_{i,n}-{\displaystyle {\sum }_{k=1}^{n-1}\left(a_{n-k-1}-a_{n-k}\right){\omega }_i^k}-a_{n-1}{\omega }_i^0\right)-h{\left(\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}\right)}_{\left(x_i,t_n\right)}\right)/h$ (2.4)

根据式(2.3)和(2.4)，我们可以得到截断误差

$h\left|{\gamma }_{i,n}^{\left(1\right)}\right|\le C{h}^3,h\left|{\gamma }_{i,n}^{\left(2\right)}\right|\le Ch\cdot t^{2-\alpha },i=1,\cdots ,N;n=1,\cdots ,L$

(2.1)式右侧第一项可以根据中点差分公式写作

$k_{\alpha }\left({\left(\frac{\partial \omega }{\partial x}\right)}_{i+\frac{1}{2},n}-{\left(\frac{\partial \omega }{\partial x}\right)}_{i-\frac{1}{2},n}\right)=k_{\alpha }\left(\widehat{{\left(\frac{\partial \omega }{\partial x}\right)}_{i+\frac{1}{2},n}}-\widehat{{\left(\frac{\partial \omega }{\partial x}\right)}_{i-\frac{1}{2},n}}\right)+h{\gamma }_{i,n}^{\left(3\right)}$ (2.5)

其中

$\widehat{{\left(\frac{\partial \omega }{\partial x}\right)}_{i+\frac{1}{2},n}}:=\frac{{\omega }_{i+1}^n-{\omega }_i^n}{h},i=0,1,\cdots ,N$ (2.6)

根据泰勒公式，容易得到

$h\left|{\gamma }_{i,n}^{\left(3\right)}\right|\le C{h}^3,i=1,2,\cdots ,N$ (2.7)

式(2.3)右侧第二项与对流速度 $f\left(x\right)$ 有关，我们利用对流项的二次向上差分即QUICK格式

$f_{i+\frac{1}{2}}{\omega }_{i+\frac{1}{2}}^n=\widetilde{f_{i+\frac{1}{2}}{\omega }_{i+\frac{1}{2}}^n}+r_{i+\frac{1}{2}},$ (2.8)

其中 $i=0,1,\cdots ,N$，

$\widetilde{f_{i+\frac{1}{2}}{\omega }_{i+\frac{1}{2}}^n}:=f_{i+\frac{1}{2}}\frac{3{\omega }_{i+1}^n+6{\omega }_i^n-{\omega }_{i-1}^n}{8},$ (2.9)

$r_{i+\frac{1}{2}}:=f_{i+\frac{1}{2}}{\omega }_{i+\frac{1}{2}}^n-\widetilde{f_{i+\frac{1}{2}}{\omega }_{i+\frac{1}{2}}^n}=-\frac{1}{2}{f}_{i+\frac{1}{2}}\left(\frac{{\partial }^3{\omega }_{i+\frac{1}{2}}^n}{\partial x^3}\right)h^3+o\left(h^4\right),$ (2.10)

式子(2.10)中是利用泰勒公式得到的，据此可以将(2.1)式右侧改写为

${\left(f\omega \right)}_{i+\frac{1}{2},n}-{\left(f\omega \right)}_{i-\frac{1}{2},n}=\widetilde{{\left(f\omega \right)}_{i+\frac{1}{2},n}}-\widetilde{{\left(f\omega \right)}_{i-\frac{1}{2},n}}-h{\gamma }_{i,n}^{\left(4\right)}$ (2.11)

其中

$h{\gamma }_{i,n}^{\left(4\right)}=\left(r_{i-\frac{1}{2}}-r_{i+\frac{1}{2}}\right)$

我们很容易得到

$\left|r_{i-\frac{1}{2}}-r_{i+\frac{1}{2}}\right|\le C{h}^3$ .

将(2.2)、(2.5)和(2.11)代入(2.1)式可以得到

$\begin{array}{l}\frac{h\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left({\omega }_i^n-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right){\omega }_i^k}-a_{n-1}{\omega }_i^0\right)\\ =k_{\alpha }\left(\widehat{{\left(\frac{\partial \omega }{\partial x}\right)}_{i+\frac{1}{2},n}}-\widehat{{\left(\frac{\partial \omega }{\partial x}\right)}_{i-\frac{1}{2},n}}\right)+\widetilde{{\left(f\omega \right)}_{i+\frac{1}{2},n}}-\widetilde{{\left(f\omega \right)}_{i-\frac{1}{2},n}}+h{\gamma }_i^n\end{array}$ (2.12)

$i=1,2,\cdots ,N;n=1,2,\cdots ,L$，其中 $r_i^n={\displaystyle {\sum }_{j=1}^4{\gamma }_{i,n}^{\left(j\right)}}$。我们用 $W_i^n$ 近似 ${\omega }_i^n$，由(2.12)式我们可以得到如下的有限体积法(FV)： $i=1,2,\cdots ,N;n=1,2,\cdots ,L$

$\begin{array}{l}\frac{h\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left(W_i^n-\underset{k=1}{\overset{n-1}{{\displaystyle \sum }}}\left(a_{n-k-1}-a_{n-k}\right)W_i^k-a_{n-1}{W}_i^0\right)\\ =k_{\alpha }\left(\widehat{{\left(\frac{\partial W}{\partial x}\right)}_{i+\frac{1}{2},n}}-\widehat{{\left(\frac{\partial W}{\partial x}\right)}_{i-\frac{1}{2},n}}\right)+\widetilde{{\left(fW\right)}_{i+\frac{1}{2},n}}-\widetilde{{\left(fW\right)}_{i-\frac{1}{2},n}}+h{\gamma }_i^n\end{array}$ (2.13)

边界条件和初始条件为

$W_0^n=g_1\left(t_n\right),W_{N+1}^n=g_2\left(t_n\right),W_i^0=\phi \left(x_i\right)$ (2.14)

其中 $\widehat{{\left(\frac{\partial W}{\partial x}\right)}_{i+\frac{1}{2},n}}$ 是直接在(2.4)式中用W替换 $\omega$， $\widetilde{{\left(fW\right)}_{i+\frac{1}{2},n}}$ 是直接在(2.8)、(2.9)和(2.11)中用W替换 $\omega$。

有限体积法(FV)的矩阵形式如下

$\left(\frac{h\Delta t^{-\alpha }}{\text{Γ}\left(2-\alpha \right)}I+A+B\right)W^n=\frac{h\Delta t^{-\alpha }}{\text{Γ}\left(2-\alpha \right)}\left({\displaystyle {\sum }_{k=1}^n\left(a_{n-k-1}-a_{n-k}\right)W^k}+a_{n-1}{W}^0\right)+d^n$ (2.15)

$n=1,2,\cdots ,L$， $W^k={\left(W_1^k,W_2^k,\cdots ,W_N^k\right)}^{\text{T}}$， $d^n=\left(d_1^n,d_2^n,\cdots ,d_N^n\right)\in R^N$， $I\in R^{N\times N}$ 是单位矩阵。 $A=\left(a_{ij}\right)\in R^{N\times N}$ 是(2.13)式右侧第一项系数矩阵， $B=\left(b_{ij}\right)\in R^{N\times N}$ 是(2.13)式右侧第二项系数矩阵，矩阵 $A,B,d^n$ 如下：

矩阵A的第1列

$a_{11}=\frac{2k_{\alpha }}{h},a_{21}=-\frac{k_{\alpha }}{h},a_{i1}=0,i\ge 3；$ (2.16)

矩阵A的第N列

$a_{NN}=\frac{2k_{\alpha }}{h},a_{\left(N-1\right)N}=-\frac{k_{\alpha }}{h},a_{iN}=0,i\le N-2；$ (2.17)

矩阵A的第i列 $\left(i=2,\cdots ,N-1\right)$

$a_{ii}=\frac{2k_{\alpha }}{h},a_{\left(i+1\right)i}=a_{\left(i-1\right)i}=-\frac{k_{\alpha }}{h},a_{ij}=0,i\ge j+2,i\le j-2；$ (2.18)

矩阵B的第1列

$b_{11}=\frac{6}{8}{f}_{\frac{3}{2}}-\frac{1}{2}{f}_{\frac{1}{2}},b_{21}=-\frac{1}{8}{f}_{\frac{5}{2}}-\frac{6}{8}{f}_{\frac{3}{2}},b_{31}=\frac{1}{8}{f}_{\frac{5}{2}},b_{i1}=0,i\ge 4$ (2.19)

矩阵B的第 $N-1$ 列

$b_{\left(N-1\right)\left(N-1\right)}=\frac{6}{8}{f}_{N-\frac{1}{2}}-\frac{3}{8}{f}_{N-\frac{3}{2}},b_{\left(N-2\right)\left(N-1\right)}=\frac{3}{8}{f}_{N-\frac{3}{2}}$ (2.20)

$b_{N\left(N-1\right)}=-\frac{1}{8}{f}_{N+\frac{1}{2}}-f_{N-\frac{1}{2}},b_{i1}=0,i\le N-3$ (2.21)

矩阵B的第N列

$b_{NN}=\frac{6}{8}{f}_{N+\frac{1}{2}}-\frac{3}{8}{f}_{N-\frac{1}{2}},b_{\left(N-1\right)N}=\frac{3}{8}{f}_{N-\frac{1}{2}}$ (2.22)

矩阵B的第j列 $2\le j\le N-2$ 列

$b_{jj}=\frac{6}{8}{f}_{j+\frac{1}{2}}-\frac{3}{8}{f}_{j-\frac{1}{2}},b_{\left(j-1\right)j}=\frac{3}{8}{f}_{j-\frac{1}{2}}$ (2.23)

$b_{\left(j+1\right)j}=-\frac{1}{8}{f}_{j+\frac{3}{2}}-\frac{6}{8}{f}_{j+\frac{1}{2}},b_{\left(j+2\right)j}=\frac{1}{8}{f}_{j+\frac{3}{2}}$ (2.24)

矩阵 $d^n$

$d_1^n=\left(\frac{1}{8}{f}_{\frac{3}{2}}+\frac{1}{2}{f}_{\frac{1}{2}}+\frac{k_{\alpha }}{h}\right)g_1\left(t_n\right),d_2^n=\left(-\frac{1}{8}{f}_{\frac{3}{2}}\right)g_1(tn)$

$d_i^n=0,i=3,\cdots ,N-1，$ (2.25)

$d_N^n=\left[-\frac{3}{8}{f}_{N+\frac{1}{2}}+\frac{k_{\alpha }}{h}\right]g_2\left(t_n\right)。$ (2.26)

3. 数值实验及结论

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

$\frac{{\partial }^{\alpha }\omega }{\partial t^{\alpha }}=\left(k_{\alpha }\frac{{\partial }^2}{\partial x^2}-\frac{\partial }{\partial x}f\left(x\right)\right)+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$，其中

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

并且 $k_{\alpha }=1$， $f\left(x\right)=x-x^2+1500$。此方程的精确解为 $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|$

空间收敛阶的数值结果列于表1~表3中。数值实验表明该方法空间上可以达到二阶收敛。

Table 1. Convergence rate for space with a = 0.2, L = 5000

表1. 空间收敛阶a = 0.2，L = 5000

Table 2. Convergence rate for space with a = 0.5, L = 5000

表2. 空间收敛阶a = 0.5，L = 5000

Table 3. Convergence rate for space with a = 0.8, L = 5000

表3. 空间收敛阶a = 0.8，L = 5000

致谢

感谢我的导师姜英军教授和师兄周帅虎、师姐黄兰对我论文完成中的悉心指导和帮助。

基金项目

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