应用数学进展  >> Vol. 8 No. 9 (September 2019)

Fokker-Planck方程的QUICK离散格式研究
QUICK Discrete Scheme for Fokker-Planck Equation

DOI: 10.12677/AAM.2019.89177, PDF, 下载: 259  浏览: 1,356  国家自然科学基金支持

作者: 尹海发:长沙理工大学数学与统计学院,湖南 长沙

关键词: 时间分数阶Fokker-Planck方程有限体积法QUICK离散格式Time Fractional Fokker-Planck Equations Finite Volume Method QUICK Scheme

摘要: 我们研究了一种求解时间分数阶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.

文章引用: 尹海发. Fokker-Planck方程的QUICK离散格式研究[J]. 应用数学进展, 2019, 8(9): 1515-1521. https://doi.org/10.12677/AAM.2019.89177

1. 研究的问题

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

ω t = ( k α 2 x 2 x f ) D t 1 α ω , a x b , 0 t T , (1.1)

初始条件和边值条件为

ω ( x , 0 ) = φ ( x ) , a x b , ω ( a , t ) = g 1 ( t ) , ω ( b , t ) = g 2 ( t ) , 0 t T , (1.2)

其中 α ( 0 , 1 ) k α 为正常数, f , φ ( x ) , g 1 ( t ) , g 2 ( t ) 为给定函数,方程(1.1)中的分数阶导数为Riemann-Liouville分数阶导数: D t 1 α ω ( x , t ) = 1 Γ ( α ) t 0 t u ( x , s ) ( t s ) 1 α d s ,其中 Γ ( x ) 是Gamma函数。

方程(1.1)常用来模拟受外力场作用下的反常扩散(如 [1] ), k α 表示广义扩散系数,f表示外力场。方程(1.1)可改写为其等价形式:

α ω t α = ( k α 2 x 2 x f ) ω ( x , t ) , a x b , 0 t T . (1.3)

其中 α ω / t α 表示 α ( 0 < α < 1 ) 阶Caputo分数阶导数: α ω t α = 1 Γ ( 1 α ) t 0 t ω ( x , η ) t d η ( t η ) α

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

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

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

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

| f ( x ) f ( x ) | C x x , x , x [ a , b ] . (1.4)

2. 离散

设N,L为正整数,取空间步长 h = ( b a ) / ( N + 1 ) ,时间步长 Δ t = T / L 。区间 [ a , b ] N + 1 等分,分点为 x i = a + i h , i = 0 , 1 , , N + 1 ,得到N个小区间 [ x i 1 2 , x i + 1 2 ] , i = 0 , 1 , , N ;在区间 [ 0 , t ] 上L等分,分点为 t k = k Δ t , k = 0 , 1 , , L 。为了方便,仅取f在点 x i + 1 2 , i = 0 , 1 , , N 处的值,记 f x + 1 2 = f ( x i + 1 2 )

在方程(1.3)中取 t = t n ( n = 0 , 1 , , L ) ,在第i个控制体积上即区间 [ x i 1 2 , x i + 1 2 ] 上对方程两边求积得到

x i 1 2 x i + 1 2 α ω t α | t n d x = k α ( ( ω x ) i + 1 2 , n ( ω x ) i 1 2 , n ) ( ( f ω ) i + 1 2 , n ( f ω ) i 1 2 , n ) , i = 1 , 2 , , N (2.1)

ω ( x i , t n ) 记做 ω i n ( i = 0 , 1 , , N + 1 ; n = 0 , 1 , , L ) ,(2.1)式左边可以改写为:

x i 1 2 x i + 1 2 α ω t α | t n d x = ( α ω t α ) ( x i , t n ) h h γ i , n ( 1 ) = h Δ t α Γ ( 2 α ) ( ω i , n k = 1 n 1 ( a n k 1 a n k ) ω i k a n 1 ω i 0 ) h γ i , n ( 1 ) h γ i , n ( 2 ) , (2.2)

(2.3)式中第一个等式应用中矩形公式,第二个等式应用了 L 1 近似,其中 a k = ( k + 1 ) 1 α k 1 α ,(见 [15] [16] ),其中误差项

γ i , n ( 1 ) : = ( ( α ω t α ) ( x i , t n ) h x i 1 2 x i + 1 2 α ω t α | t n ) / h (2.3)

γ i , n ( 2 ) : = ( h Δ t α Γ ( 2 α ) ( ω i , n k = 1 n 1 ( a n k 1 a n k ) ω i k a n 1 ω i 0 ) h ( α ω t α ) ( x i , t n ) ) / h (2.4)

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

h | γ i , n ( 1 ) | C h 3 , h | γ i , n ( 2 ) | C h t 2 α , i = 1 , , N ; n = 1 , , L

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

k α ( ( ω x ) i + 1 2 , n ( ω x ) i 1 2 , n ) = k α ( ( ω x ) i + 1 2 , n ^ ( ω x ) i 1 2 , n ^ ) + h γ i , n ( 3 ) (2.5)

其中

( ω x ) i + 1 2 , n ^ : = ω i + 1 n ω i n h , i = 0 , 1 , , N (2.6)

根据泰勒公式,容易得到

h | γ i , n ( 3 ) | C h 3 , i = 1 , 2 , , N (2.7)

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

f i + 1 2 ω i + 1 2 n = f i + 1 2 ω i + 1 2 n ˜ + r i + 1 2 , (2.8)

其中 i = 0 , 1 , , N

f i + 1 2 ω i + 1 2 n ˜ : = f i + 1 2 3 ω i + 1 n + 6 ω i n ω i 1 n 8 , (2.9)

r i + 1 2 : = f i + 1 2 ω i + 1 2 n f i + 1 2 ω i + 1 2 n ˜ = 1 2 f i + 1 2 ( 3 ω i + 1 2 n x 3 ) h 3 + ο ( h 4 ) , (2.10)

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

( f ω ) i + 1 2 , n ( f ω ) i 1 2 , n = ( f ω ) i + 1 2 , n ˜ ( f ω ) i 1 2 , n ˜ h γ i , n ( 4 ) (2.11)

其中

h γ i , n ( 4 ) = ( r i 1 2 r i + 1 2 )

我们很容易得到

| r i 1 2 r i + 1 2 | C h 3 .

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

h Δ t α Γ ( 2 α ) ( ω i n k = 1 n 1 ( a n k 1 a n k ) ω i k a n 1 ω i 0 ) = k α ( ( ω x ) i + 1 2 , n ^ ( ω x ) i 1 2 , n ^ ) + ( f ω ) i + 1 2 , n ˜ ( f ω ) i 1 2 , n ˜ + h γ i n (2.12)

i = 1 , 2 , , N ; n = 1 , 2 , , L ,其中 r i n = j = 1 4 γ i , n ( j ) 。我们用 W i n 近似 ω i n ,由(2.12)式我们可以得到如下的有限体积法(FV): i = 1 , 2 , , N ; n = 1 , 2 , , L

h Δ t α Γ ( 2 α ) ( W i n k = 1 n 1 ( a n k 1 a n k ) W i k a n 1 W i 0 ) = k α ( ( W x ) i + 1 2 , n ^ ( W x ) i 1 2 , n ^ ) + ( f W ) i + 1 2 , n ˜ ( f W ) i 1 2 , n ˜ + h γ i n (2.13)

边界条件和初始条件为

W 0 n = g 1 ( t n ) , W N + 1 n = g 2 ( t n ) , W i 0 = φ ( x i ) (2.14)

其中 ( W x ) i + 1 2 , n ^ 是直接在(2.4)式中用W替换 ω ( f W ) i + 1 2 , n ˜ 是直接在(2.8)、(2.9)和(2.11)中用W替换 ω

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

( h Δ t α Γ ( 2 α ) I + A + B ) W n = h Δ t α Γ ( 2 α ) ( k = 1 n ( a n k 1 a n k ) W k + a n 1 W 0 ) + d n (2.15)

n = 1 , 2 , , L W k = ( W 1 k , W 2 k , , W N k ) T d n = ( d 1 n , d 2 n , , d N n ) R N I R N × N 是单位矩阵。 A = ( a i j ) R N × N 是(2.13)式右侧第一项系数矩阵, B = ( b i j ) R N × N 是(2.13)式右侧第二项系数矩阵,矩阵 A , B , d n 如下:

矩阵A的第1列

a 11 = 2 k α h , a 21 = k α h , a i 1 = 0 , i 3 (2.16)

矩阵A的第N列

a N N = 2 k α h , a ( N 1 ) N = k α h , a i N = 0 , i N 2 (2.17)

矩阵A的第i列 ( i = 2 , , N 1 )

a i i = 2 k α h , a ( i + 1 ) i = a ( i 1 ) i = k α h , a i j = 0 , i j + 2 , i j 2 (2.18)

矩阵B的第1列

b 11 = 6 8 f 3 2 1 2 f 1 2 , b 21 = 1 8 f 5 2 6 8 f 3 2 , b 31 = 1 8 f 5 2 , b i 1 = 0 , i 4 (2.19)

矩阵B的第 N 1

b ( N 1 ) ( N 1 ) = 6 8 f N 1 2 3 8 f N 3 2 , b ( N 2 ) ( N 1 ) = 3 8 f N 3 2 (2.20)

b N ( N 1 ) = 1 8 f N + 1 2 f N 1 2 , b i 1 = 0 , i N 3 (2.21)

矩阵B的第N列

b N N = 6 8 f N + 1 2 3 8 f N 1 2 , b ( N 1 ) N = 3 8 f N 1 2 (2.22)

矩阵B的第j列 2 j N 2

b j j = 6 8 f j + 1 2 3 8 f j 1 2 , b ( j 1 ) j = 3 8 f j 1 2 (2.23)

b ( j + 1 ) j = 1 8 f j + 3 2 6 8 f j + 1 2 , b ( j + 2 ) j = 1 8 f j + 3 2 (2.24)

矩阵 d n

d 1 n = ( 1 8 f 3 2 + 1 2 f 1 2 + k α h ) g 1 ( t n ) , d 2 n = ( 1 8 f 3 2 ) g 1 (tn)

d i n = 0 , i = 3 , , N 1 (2.25)

d N n = [ 3 8 f N + 1 2 + k α h ] g 2 ( t n ) (2.26)

3. 数值实验及结论

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

α ω t α = ( k α 2 x 2 x f ( x ) ) + g ( x , t ) , 0 x 1 , 0 t 1 (3.1)

初始条件和边值条件为 u ( x , 0 ) = 0 g 1 ( t ) = t 2 g 2 ( t ) = t 2 ,其中

g ( x ) = Γ ( 3 ) Γ ( 3 α ) t 2 α cos ( π x ) + k α π 2 t 2 cos π x + t 2 [ ( 1 2 x ) cos ( π x ) f ( x ) π sin ( π x ) ] , (3.2)

并且 k α = 1 f ( x ) = x x 2 + 1500 。此方程的精确解为 u ( x , t ) = t 2 cos ( π x )

我们定义空间收敛阶如下:

= | ln ( 1 / 1 ) ln ( N + 1 / N + 1 ) |

空间收敛阶的数值结果列于表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)。

参考文献

[1] So, F. and Liu, K.L. (2004) A Study of the Subdiffusive Fractional Fokker-Planck Equation of Bistable Systems. Physica A: Statistical Mechanics and Its Applications, 331, 378-390.
https://doi.org/10.1016/j.physa.2003.09.026
[2] Jiang, Y. and Xu, X. (2018) A Monotone Finite Volume Meth-od for Time Fractional Fokker-Planck Equations. Science China Mathematics, 62, 783-794.
https://doi.org/10.1007/s11425-017-9179-x
[3] Jiang, Y. (2015) A New Analysis of Stability and Convergence for Finite Difference Schemes Solving the Time Fractional Fokker-Planck Equation. Applied Mathematical Modelling, 39, 1163-1171.
https://doi.org/10.1016/j.apm.2014.07.029
[4] Jiang, Y. and Ma, J. (2011) High-Order Finite Element Methods for Time-Fractional Partial Differential Equations. Journal of Computational and Applied Mathematics, 235, 3285-3290.
https://doi.org/10.1016/j.cam.2011.01.011
[5] Chen, S., Liu, F., Zhuang, P. and Anh, V. (2009) Finite Difference Approximations for the Fractional Fokker-Planck Equation. Applied Mathematical Modelling, 33, 256-273.
https://doi.org/10.1016/j.apm.2007.11.005
[6] Deng, W. (2007) Numerical Algorithm for the Time Fractional Fokker-Planck Equation. Journal of Computational Physics, 227, 1510-1522.
https://doi.org/10.1016/j.jcp.2007.09.015
[7] Cao, X., Fu, J. and Huang, H. (2012) Numerical Method for the Time Fractional Fokker Planck Equation. Advances in Applied Mathematics and Mechanics, 4, 848-863.
[8] Saadatmandi, A., Dehghan, M. and Azizi, M. (2012) The Sinc-Legendre Collocation Method for a Class of Fractional Convection-Diffusion Equations with Variable Coefficients. Communications in Nonlinear Science and Numerical Simulation, 17, 4125-4136.
https://doi.org/10.1016/j.cnsns.2012.03.003
[9] Fairweather, G., Zhang, H., Yang, X. and Xu, D. (2014) A Backward Euler Orthogonal Spline Collocation Method for the Time-Fractional, Fokker-Planck Equation. Numerical Methods for Partial Differential Equations, 31, 1534-1550.
https://doi.org/10.1002/num.21958
[10] Vong, S. and Wang, Z. (2015) A High Order Compact Finite Difference Scheme for Time Fractional Fokker-Planck Equation. Applied Mathematics Letters, 43, 38-43.
https://doi.org/10.1016/j.aml.2014.11.007
[11] Le, K.N., Mclean, W. and Mustapha, K. (2016) Numerical So-lution of the Time-Fractional Fokker-Planck Equation with General Forcing. SIAM Journal on Numerical Analysis, 54, 1763-1784.
https://doi.org/10.1137/15M1031734
[12] Feng, L.B., Zhuang, P., Liu, F. and Turne, I. (2015) Stabil-ity and Convergence of a New Finite Volume Method for a Two-Sided Space Fractional Diffusion Equation. Applied Mathematics and Computation, 257, 52-65.
https://doi.org/10.1016/j.amc.2014.12.060
[13] Hejazi, H., Moroney, T. and Liu, F. (2013) A Finite Volume Method for Solving the Two-Sided Time-Space Fractional Advection-Dispersion Equation. Central European Journal of Physics, 11, 1275-1283.
https://doi.org/10.2478/s11534-013-0317-y
[14] Karaa, S., Mustapha, K. and Pani, A.K. (2016) Finite Volume Element Method for Two-Dimensional Fractional Sub-Diffusion Problems. IMA Journal of Numerical Analysis, 37, 945-964.
https://doi.org/10.1093/imanum/drw010
[15] Langlands, T. and Henry, B. (2005) The Accuracy and Stability of an Implicit Solution Method for the Fractional Diffffusion Equation. Journal of Computational Physics, 205, 719-736.
https://doi.org/10.1016/j.jcp.2004.11.025
[16] Oldham, K.B. and Spanier, J. (1974) The Fractional Calculus. Academic Press, New York.