Advances in Applied Mathematics
Vol. 08  No. 02 ( 2019 ), Article ID: 28840 , 7 pages
10.12677/AAM.2019.82023

A Finite Volume Method for Time Fractional Fokker-Planck Equations

Lan Huang

School of Mathematics and Statistics, Changsha University of Science and Technology, Changsha Hunan

Received: Jan. 22nd, 2019; accepted: Feb. 6th, 2019; published: Feb. 13th, 2019

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.

Keywords:Time Fractional Fokker-Planck Equations, Finite Volume Method

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

黄兰

长沙理工大学数学与统计学院,湖南 长沙

收稿日期:2019年1月22日;录用日期:2019年2月6日;发布日期:2019年2月13日

摘 要

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

关键词 :时间分数阶Fokker-Planck方程,有限体积法

Copyright © 2019 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 研究的问题

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

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

初始条件和边值条件为

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

其中 α ( 0 , 1 ) k α 为正常数, f ( x , t ) , φ ( x ) , g 1 ( t ) , g 2 ( t ) 为给定函数,方程(1.1)中的分数阶导数为Riemann-Liouville分数阶导数 D t 1 α u ( x , t ) = 1 Γ ( α ) t 0 t u ( x , s ) ( t s ) 1 α d s Γ ( x ) 是Gamma函数。方程(1.1)可用于模拟外力场作用下的低扩散现象(参见 [1] )。

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

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

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

2. 离散

[ 0 , T ] 上取分点为 t k = k Δ t k = 0 , 1 , , L ,其中时间步长 Δ t = T / L ,L为正整数。在区间 [ a , b ] 上取分点 x i = a + i h , i = 0 , 1 , , N + 1 ,其中空间步长 h = ( b a ) / ( N + 1 ) ,N为正整数。N个空间有限体为 [ x i 1 2 , x i + 1 2 ] , i = 1 , , N ,其中 x i + 1 2 = x i + x i + 1 2 , i = 0 , 1 , , N

在方程(1.1)中取 t = t n ( n = 0 , 1 , , L ) 得到

u t | t n = [ k α 2 x 2 x f ( x , t n ) ] D t 1 α u ( x , t n ) . (2.1)

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

u t | t n = u ( x , t n ) u ( x , t n 1 ) Δ t + r n ( 1 ) , (2.2)

其中 r n ( 1 ) = u t | t n u ( x , t n ) u ( x , t n 1 ) Δ t 满足 | r n ( 1 ) | C Δ t

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

D t 1 α u ( x , t n ) = Δ t α 1 j = 0 n ω n j 1 α u ( x , t j ) + r n ( 2 ) , (2.3)

其中 ω k α = ( 1 ) k Γ ( α + 1 ) Γ ( k + 1 ) Γ ( α k + 1 ) ( k = 0 , 1 , 2 , ) r n ( 2 ) = D t 1 α u ( x , t n ) Δ t α 1 j = 0 n ω n j 1 α u ( x , t j ) 满足 r n ( 2 ) = C ( x ) Δ t C ( x ) 是与解u有关的光滑函数。

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

u ( x , t n ) u ( x , t n 1 ) Δ t + r n ( 1 ) = [ k α 2 x 2 x f ( x , t n ) ] [ Δ t α 1 j = 0 n ω n j 1 α u ( x , t j ) + r n ( 2 ) ] . (2.4)

u n ( x ) 表示 u ( x , t n ) 的近似,由(2.4)我们得到原问题的时间半离散格式

u n ( x ) u n 1 ( x ) Δ t = [ k α 2 x 2 x f ( x , t n ) ] [ Δ t α 1 j = 0 n ω n j 1 α u j ( x ) ] , (2.5)

v n ( x ) : = Δ t α 1 j = 0 n ω n j 1 α u j ( x ) . (2.6)

易知

u n ( x ) : = Δ t 1 α j = 0 n ω n j α 1 v j ( x ) . (2.7)

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

Δ t α [ j = 0 n 1 ( ω n j α 1 ω n j 1 α 1 ) v j ( x ) + ω 0 α 1 v n ( x ) ] = [ k α 2 x 2 x f ( x , t n ) ] v n ( x ) . (2.8)

在有限体 [ x i 1 2 , x i + 1 2 ] 上对(2.8)式积分

x i 1 2 x i + 1 2 Δ t α [ j = 0 n 1 ( ω n j α 1 ω n j 1 α 1 ) v j ( x ) + ω 0 α 1 v n ( x ) ] d x = x i 1 2 x i + 1 2 [ k α 2 x 2 x f ( x , t n ) ] v n ( x ) d x . (2.9)

利用中矩形积分公式,用 v i n 表示 v n ( x i ) ,(2.9)式左侧可以写成

x i 1 2 x i + 1 2 Δ t α [ j = 0 n 1 ( ω n j α 1 ω n j 1 α 1 ) v j ( x ) + ω 0 α 1 v n ( x ) ] d x = h Δ t α [ j = 0 n 1 ( ω n j α 1 ω n j 1 α 1 ) v i j + ω 0 α 1 v i n ] h r i , n ( 1 ) , (2.10)

其中

h | r i , n ( 1 ) | C h 3 . (2.11)

f i + 1 / 2 , n 表示 f ( x i + 1 / 2 , t n ) ,(2.9)式右侧可以写成

x i 1 2 x i + 1 2 [ k α 2 x 2 x f ( x , t n ) ] v n ( x ) d x = k α ( v n x | x i + 1 2 v n x | x i 1 2 ) ( f i + 1 2 , n v i + 1 2 n f i 1 2 , n v i 1 2 n ) . (2.12)

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

k α ( v n x | x i + 1 2 v n x | x i 1 2 ) = k α ( v i + 1 n v i n h v i n v i 1 n h ) + h r i , n ( 2 ) , (2.13)

其中

h | r i , n ( 2 ) | C h 3 ; (2.14)

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

f i + 1 2 , n v i + 1 2 n f i 1 2 , n v i 1 2 n = f i + 1 2 , n v i n + v i + 1 n 2 f i 1 2 , n v i n + v i 1 n 2 h r i , n ( 3 ) , (2.15)

其中

h | r i , n ( 3 ) | C h 3 . (2.16)

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

h Δ t α [ j = 0 n 1 ( ω n j α 1 ω n j 1 α 1 ) v i j + ω 0 α 1 v i n ] h r i , n ( 1 ) = k α ( v i + 1 n v i n h v i n v i 1 n h ) + h r i , n ( 2 ) ( f i + 1 2 , n v i n + v i + 1 n 2 f i 1 2 , n v i n + v i 1 n 2 ) + h r i , n ( 3 ) , (2.17)

V i n 近似 v i n ,由(2.17)式我们可以得到如下的有限体积法(FV): i = 1 , 2 , , N ; n = 1 , 2 , , L

h Δ t α [ j = 0 n 1 ( ω n j α 1 ω n j 1 α 1 ) V i j + ω 0 α 1 V i n ] = k α ( V i + 1 n V i n h V i n V i 1 n h ) ( f i + 1 2 , n V i n + V i + 1 n 2 f i 1 2 , n V i n + V i 1 n 2 ) , (2.18)

初始条件和边值条件为 ( i = 1 , 2 , , N ; n = 1 , 2 , , L )

V i 0 = φ ( x i ) , V 0 n = G 1 ( t n ) = Δ t α 1 j = 0 n ω n j 1 α g 1 ( t n ) , V N + 1 n = G 2 ( t n ) = Δ t α 1 j = 0 n ω n j 1 α g 2 ( t n ) . (2.19)

3. 数值实验

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

u t = ( 2 x 2 x f ( x , t ) ) D t 1 α u + 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 , t ) = 2 t cos ( π x ) + 2 π 2 Γ ( 2 + α ) t 1 + α cos ( π x ) + 2 Γ ( 2 + α ) t 1 + α [ ( 1 2 x + t ) cos ( π x ) f ( x , t ) π sin ( π x ) ] , (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 ( x , t ) = x x 2 + t + x t 。此方程的精确解为 u ( x , t ) = t 2 cos ( π x )

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

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

= | ln ( 1 / 1 ) ln ( L / L ) | .

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

致谢

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

基金项目

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

文章引用

黄 兰. 时间分数阶Fokker-Planck方程有限体积法
A Finite Volume Method for Time Fractional Fokker-Planck Equations[J]. 应用数学进展, 2019, 08(02): 203-209. https://doi.org/10.12677/AAM.2019.82023

参考文献

  1. 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. 2. Jiang, Y. and Xu, X. (2018) A Monotone Finite Volume Method for Time Fractional Fokker-Planck Equations. Science China Mathematics, 1-12. https://doi.org/10.1007/s11425-017-9179-x

  3. 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. 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. 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. 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. 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. 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. 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. 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. 11. Le, K.N., Mclean, W. and Mustapha, K. (2016) Numerical Solution 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. 12. Feng, L.B., Zhuang, P., Liu, F. and Turne, I. (2015) Stability and Convergence of a New Finite Volume Method for a Two-Sided Space Fractional Diffusion Equation. Applied Mathe-matics and Computation, 257, 52-65. https://doi.org/10.1016/j.amc.2014.12.060

  13. 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.

  14. 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. 15. 郭柏灵, 蒲学科, 黄凤辉. 分数阶偏微分方程及其数值解[M]. 北京: 科学出版社, 2011: 198-219.

期刊菜单