Operations Research and Fuzziolgy 运筹与模糊学, 2012, 2, 1-7 http://dx.doi.org/10.12677/orf.2012.21001 Published Online February 2012 (http://www.hanspub.org/journal/orf) Finite Difference Methods for Space-Time Fractional Two-Sided Space Partial Differential Equations* Yang Zhang, Ningping Li, Lu Chen School of Mathematical Science, Nankai University, Tianjin Email: zhangyang@nankai.edu.cn Received: Dec. 4th, 2011; revised: Jan. 6th, 2012; accepted: Jan. 8th, 20 12 Abstract: Fractional order differential equations are generalizations of classical differential equations. Now, they are widely used in the fields of physics, information; finance and others. In this paper, an explicit finite difference method for space-time fractional two-sided space partial differential equations is established. Be- sides, the stability and co nvergence order are analyzed. Keywords: Fractional Derivative; Explicit Methods; Stability; Convergence; Error Estimates 分数阶导数双边空间微分方程的显式差分解法* 张 阳,李宁平,陈 璐 南开大学数学科学学院,天津 Email: zhangyang@nankai.edu.cn 收稿日期:2011 年12 月4日;修回日期:2012 年1月6日;录用日期:2012 年1月8日 摘 要:分数阶微分方程作为广义的微分方程,被广泛地应用于物理,信息处理,金融等领域。本文 给出了数值求解时间空间分数阶导数的双边空间微分方程的一种显式差分格式,并对其稳定性和收敛 性进行了理论分析。 关键词:分数阶导数;显格式;稳定性;收敛性;误差估计 1. 引言 分数阶导数微分方程,是传统的整数阶微分方程的推广。相比整数阶导数微分方程,分数阶导数微分方程 能够更好拟合某些自然物理过程和动态系统过程,近年来被广泛应用于反常扩散,粘弹性本构建模,信号处理, 控制,流体力学,图像处理,软物质研究,金融等领域[1-8]。近年来,分数阶导数微积分的研究与兴起,伴随着 数字计算机计算技术的提高而迅速发展。双边空间分数阶对流–扩散方程是一种反常扩散,在地下水溶质运移 方面,它可描述含水层中溶质运移过程中的反常扩散现象。关于这类问题的数值解法 Meerschaert 等人分别对单 边对流–扩散方程和双边扩散方程,进行差分求解,文[9]针对变系数反应扩散方程 ,,, d, crtcrt crt vrrf rt trr *资助信息:该工作由中国国家自然科学基金资助(11171168,11071132)和中国高校博士点基金资助(20100031110002)。 Copyright © 2012 Hanspub 1 张阳 等 分数阶导数双边空间微分方程的显式差分解法 建立显式和隐式两种差分解法,文[10,11]建立了双边空间分数阶导数微分方程 ,,, ,, uxtuxtuxt cxtcsxt txx 数值解法并进行了稳定性分析,文[12]针对二维分数阶扩散方程 ,,,, ,, ,, uxytuxyt uxyt dxy exyqxyt txx ,, 建立了交替方向隐式差分解法;夏源、吴吉春在[13]中对一维时间分数阶对流–弥散方程 2 2 ,,cxt cxtcxt vD , x tx 给出了差分算法;刘发旺等人用另一种方法对这个问题进行离散计算[14]。周璐莹等在[15]中针对二维分数阶对 流–弥散方程 22 22 ,,,, ,,,,,,cxytcxyt cxytcxytcxyt vD xytx y 建立了数值解法。[16]则将时间空间分数阶导数微分方程 ,,, ,, uxtuxt uxt bxaxcxuxt qxt txx 的隐式差分解法进行了研究。目前,关于带时间分数阶导数的双边空间分数阶导数微分方程的数值解法尚未见 到,本文正是针对该类问题建立了相应的显式数值解法并进行理论分析,涵盖了[10]和[13]的类型,且采用两种 理论方法分析,较已有的方法更为完善。采用显式差分解法可以减少计算量,但理论分析较为复杂,从而本文 得到的结论丰富了这一领域的研究成果。 本文结构如下:第二节针对研究的对象进行了描述,并给出了两个常见的引理;第三节建立了所描述方程 的显式差分解法并进行了方法的稳定性分析;第四节给出了方法的收敛性分析和收敛阶的估计;第五节总结了 结论与展望。 2. 问题描述 本文针对带分数阶时间空间导数的双边空间微分方程建立显式差分解法并进行理论分析,研究如下方程: ,,,, ,, uxtuxtuxt uxt kxtcxtcs xt txxx ,, (1) 这里, 。其中 LxR 0tT 01 ,01 ,12 , ,kxt, ,cxt , 均为非负有界函数, 并且假定初始条件: ,以及边界条件 ,cxt ,0t f xuxLxR ,,Rt 0ux Ltux 。 方程(1)右端的左边(+)以及右边(–)的分数阶导数为 Riemann-Liouville 分数阶导数。按如下方式定义[10]: 1 d1d d dd nx n L fx f n xx x n (2) 和 1 d1 dd. dd nnR n x fx f nxx x n (3) 这里 是正整数,并且满足n1nn 。空间分数阶导数 ,uxt x 是 次Riemman-Liouville分数阶导数, 2 Copyright © 2012 Hanspub 张阳 等 分数阶导数双边空间微分方程的显式差分解法 定义为: 0 , (,) 1d, 01. 1 xut uxt x xx (4) 方程(1)左端的 ,uxt t 由Caputo 分数阶导数给出定义[17]: 0 ,, 1d, 1 t uxt ux t t (5) 这里是 Gamma 函数。 () 根据 Grünwald 估计,为了得到稳定的数值方法[9],定义左右两边分数阶导数分别按如下形式近似(Shifted Grünwald): 0 d1 lim 1 d M k Mk fx g fx kh xh (6) 和 0 d()1 lim 1 d M k Mk fx g fx kh xh (7) 这里 M , M 均为正整数, x L h M ,Rx h M , 11 1, ! k k k gk k 0,1,2 。 引理 1. 上式定义的 k g 满足如下条件: 1) ,2) 01g10g ,3) 11 10 ! k k k, >2 g k k 。 对于引理 1,结论 1)、2)是显然的,下面只说明结论 3)。 因为12 ,所以 0 ,10 ,0k , 。从而有 2k 2 22 11 121 111 !! 12 1 10 ! kkk k k kk gkk k k , 0 引理成立。对于 ,根据二项式定理, 有 0 1, m m zz m 1z 1z。令 ,有。由引理 1知 0 0 j j g 10, 1 , 1,2,, i j jj g gi m (8) 3. 显格式及其稳定性分析 先建立方程(1)在前文所述初边值条件下的显格式并通过圆盘定理证明该差分格式是无条件稳定的,并且还 将采用另一种方法来证明格式的稳定性及收敛性,并给出收敛阶的估计。 引入步长 h,时间步长 ,则有 n , >0; ,0 j xLjhh tn 。其中 ,0 TRL h K m 。对于 , a a uxt x , , a a uxt x 用shifted Grünwald形式,有如下表达式[9]: 1 0 ,11, , i nki n k uxt g ux khth xh (9) Copyright © 2012 Hanspub 3 张阳 等 分数阶导数双边空间微分方程的显式差分解法 1 0 ,11, , mi nki n k uxt g ux khth xh (10) 12 0 ,1, i ki n k uxt , g ux khth xh (11) 其中 2 (2) 0 11 1, 1,1,2, ! k k j gg k k 易验证 2 k g 具有类似 k g 的性质: 且 22 2 01 1,0,0 2 k ggg k 2 1 k j jg 2,3, ,km。 对于时间分数阶导数有如下估计[17]: 1 11 0 1 1 01 11 11 0 ,, 1d 1 ,, 1d 1 ,, 1 2 n t in i n nj ij ij j jn ninj inj j uxt ux t t uxt uxt t uxtuxt jj (12) 故可对方程(1)建立如下的显格式: 1 111 12 1,1 ,1 000 1 1 2 nj njn niim nnnnn ii i 0 i n j ikki ikki iki jkkk uu k jjgugcu gcu hh s (13) 易知,该方法是相容的,其舍入误差为 h 。记 11 1 j cj j ,则有 引理 2. 上面所定义的 11 1 j cj j 满足以下条件 1) , 2) , 3) 0, 0,1, j cj 1,0,1, jj cc j 1 1+1 1 22+1 k jjk j cc c ,即 1 1 1 + 1 k jk j dc 1 22 对于 1),3) 是显然的,2) 利用单调性也可以得到。 令 ,, 22 2 ,, nn n ii i iii cc b rhhh , 则当 时,格式(13)为 0n 11 100000 02111 20111 33 1+() imi n iiiiiiiiiiiikiikkiik kk uggurg gurggugugus i (14) 当 时,格式(13)为 0 n 11 02111 20 111 0 111 331 22 nn n iiii iiiiiii imin nnnjn kiikkiikjinii kk j uggurg gurggu gugudu cus 1 n i , (15) 因此格式(14)、(15)有如下等价形式 100 0 1110 2 0 , ,0 , nnn n nn UHU S UHUdU dUcUSn Uf (16) 4 Copyright © 2012 Hanspub 张阳 等 分数阶导数双边空间微分方程的显式差分解法 其中矩阵 ij H , 有如下表达式 1,1;1, ,1;imj m 111 20 02 1 1 22 1 1 1 iii iii iji i iji ji i rg gji rg gji g gji g ji g ji , 当 ,当 , 当 , 当 , 当 H 1 其中 ,,; 00 1H00 j H1, 2,,jm1 mm H ,0 mj H ,0,1, ,1jm 。 定理 1. 当01 ,12 时,对于方程(1)所建立的显格式(16)在满足如下条件: 1,, ,, 22 22maxmaxmax 22 maxmaxmax1 0 ii ii cc b h h cc b h h (17) 时是稳定的。 证明. 根据 Gerschgorin 定理,以及引理 1,知矩阵 H的特征值所在中心为 11 1 22 22 iiii iii i rgr H , 半径 11 11 0,0, 10, 1 () ()() mimi iikikikiiiiii kkikkkk Hggrrggr. i 根据(8)式,有 11 22 2222 1, ii i r 由(17)式知,,即 H的特征值不超过1。 1 22 221, ii i r 对于 时的情形,由于0n0 H 的特征值 满足 122 1 ii i r ,r0 ,又因为 故 ii 2i r 1 122 221 i i 2 i H 的特征值亦不超过 1,从而稳定性得证。 下面用另一种方法给出关于最大模的稳定性分析。假设 1, 2,,1;1, 2,, j i uimjK ,是 (15)所得的近似 值,误差(扰动)为 1, 2,,1;1, 2,, jjj iii uuim jK ,并且满足: 当 时, 0n 11 101000 021112011 33 22 ; imi iiiiiiii iiiikiikkiik kk ggrgg rgggg 0 1 当 时, 0n 11 0211120 111 0 111 331 22 . nn n iiii iiiiiii imin nnnj ki ikki ikjini kk j ggrgg rgg ggdc 1 n i 令,。 12 1 ,,,T kkk km E 0,1,2, ,kK 定理 2. 在条件(16)下,有 0,0,1,2, k EEk n;这里11 11 max nn il im E1n 。 证明. 用数学归纳法。 1) 当时,令0n1 11 max l im 1 i ,由(16)知 1 22 0 ll l r ,故有 1 22 0 ll l r 。再 Copyright © 2012 Hanspub 5 张阳 等 分数阶导数双边空间微分方程的显式差分解法 由引理 1,则有 10 0 021 1120 11 11 00 0 11 33 00 1 1. llllll llllll lml lml lk lklk lklklk kk kk ggrgg rgg gg ggE 0 1 由引理 1-3)知 ,,故 1 0 0 l k k g 1 0 0 ml k k g 10 lE ,即 0k EE 。 2) 假设 0,0,1,2, n EEn k成立,考虑 1nk 的情形。 令11 11 max kk il im E1k ,类似于 1),有 11 02111 201 111 111 01 1111 331 001 11 22 22 (2 2) kkk k llllllllllll lmlklmlk kkkj ls lslslsjlkllklkjk ss jkkj jk ggrgg rgg ggdc ggdc dc 0 E 100 1 k j EE 定理得证。 由上述定理知格式(16)条件稳定。 4. 显格式的收敛性分析与收敛阶估计 下面进行收敛性分析。令 是方程(1)在相应初边值条件下于点 ,1,2,,1;1,2, ik uxt imkK , ik x t的 精确解。定义 ,1,2,,1;1,2, kk iiki euxt uimkK并且 12 1 ,,, T kkk k m Yeee ,则有如下表达式 11 02111 20 11 0 111 331 0 22( ) ,1, 0, kk k iiiiiii iiii i imin kkkj siissiisjiki ssj i eggerg gergge gegedeck e 1 k K 1, 2,,1;0,1, 2,imk。记 为局部截断误差,则存在常数,使得 1k i R0C 1, 1,2,,1;0,1,2,. k i RC himk K 定理 3. 当满足条件(15)时, 11kk YCc h0,1,2, 1kK , 。这 里11 max kk e i im Y ,C为正常数。 证明. 用数学归纳法。 1) 当 时,令0k1 11 max l im e 1 i e,根据引理 1,则有 10 1 0ll eRCc h ,即 11 0 YCc h。 2) 假设 11,1,2, nn YCc hn k;成立,由(17)知 1 22 0 ll l r ,考 虑的情形。 令 1nk 11 11 max kk iL im Ye 1k e ,利用引理 1,引理2,可以得到 11 02111 20 11 1 0 11 33 1 111 11 1 001 22 1 22 kk k lllllll l llll lml k kkjk ls lsjlkll ss j lmlk lkl kjkkk kkj eggerg gergge k geged eceR lsl s ggdcCchCc ,h 1 1 k 6 Copyright © 2012 Hanspub 张阳 等 分数阶导数双边空间微分方程的显式差分解法 Copyright © 2012 Hanspub 7 故 11 0 k YCc h,定理得证。 根据定理 3,并注意到 11 11 1 1 lim limlim1 11 11 k kk k ckk kkk k ,因此存在常数 C,使得 11kk YCc h ,因 为kT ,进而 1k YCT h ,从而有下述定理 定理 4. 假设 是方程(1)在相应初边值条件下在点 ,1,2,,1;1,2, ik uxt imkK , ik x t处的精确解。 是格式(16)所得到的数值解,则在满足条件(17)时,存在正常数,使得 k i u * C * ,,1,2,, 1; k ik i uxtuCh imkK 1,2,. 5. 结论与展望 本文将[10]中讨论的双边空间分数阶导数微分方程的时间导数推广到分数阶,建立了一种显式差分求解格式 并进行了相应的理论分析,得到了丰满的一阶时间和空间的收敛阶逼近。所得结论亦可推广到方程右端含项及 更为复杂的情形,所建立的方法也不难推广到高维情形,至于高维问题如何进行算子分裂尚有待进一步探讨, 对于所讨论方程隐格式的建立与理论分析作者将另文讨论。最后,作者对审稿人对稿件提出的修改意见表示感 谢。 参考文献 (References) [1] R. Metzler, J. Klafter. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. Journal of Physics, 2004, A37: 161-208. [2] Z. Deng, V. P. Singh and L. Bengtsson. Numerical solution of fractional advection-dispersion equation. Journal of Hydraulic Engineering, 2004, 130(5): 422-431. [3] R. Gorenflo, F. Mainardi, E. Scalas and M. Raberto. Fractional calculus and continuous-time finance III: The diffusion limit. In: Kolhmann, S. Tang, Eds., Mathematical Finance, Basel: Birkhauser Verlag, 2001: 171-180. [4] 苏丽娟, 王文洽. 双边分数阶对流–扩散方程的一种有限差分解法[J]. 山东大学学报(理学版), 2009, 10: 29-32. [5] F. Liu, P. Zhuang, V. Anh, I. Turner and K. Burrage. Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation. Applied Mathematics and Computation, 2007, 191(1): 12-21. [6] 周激流, 蒲亦非, 廖科. 分数阶微积分原理及其在现代信号分析与处理中的应用[M]. 北京: 科学出版社, 2010. [7] R. Hilfer. Application of fractional calculus in physics. Singapore, New Jersey, London and Hong Kong: World Scientific Publication Company, 2000. [8] A. A. Kibas, H. M. Srivastava and J. J. Trujillo. Theory and application of fractional differential equations. Amsterdam: Elser vier, 2006. [9] M. M. Meerschaert, C. Tadjeran. Finite difference approximations for fractional advection-dispersion flow equations. Journal of Computational & Applied Mathematics, 2004, 172(1): 65-77. [10] M. M. Meerschaert, C. Tadjeran. Finite difference approximations for two-sided space-fractional partial differential equations. Applied Nu- merical Mathematics, 2006, 56(1): 80-90. [11] V. K. Tuan, R. Gorenflo. Extrapolation to the limit for numerical fractional differentiation. Zeitschrift für Angewandte Mathematik und Me- chanik, 1995, 75: 646-648. [12] M. M. Meerschaert, H. P. Scheffler and C. Tadjeran. Finite difference method for two dimensional fractional dispersion equations. Journal of Computational Physics, 2006, 211: 249-261. [13] 夏源, 吴吉春. 分数阶对流—弥散方程的数值求解[J]. 南京大学学报(自然科学版), 2007, 43(4): 44-446. [14] F. Liu, V. Ahn and I. Turner. Numerical solution of the space fractional Fokker-Planck equation. Journal of Computational & Applied Mathematics, 2004, 166(1): 209-219. [15] 周璐莹, 吴吉春, 夏源. 二维分数阶对流–弥散方程的数值解[J]. 高校地质学报, 2009, 15(4): 569-575. [16] Y. Zhang. A finite difference method for fractional partial differential equation. Applied Mathematics and Computation, 2009, 215(2): 524- 529. [17] M. M. Meerschaert, D. A. Benson, H. P. Scheffler and B. Baeumer. Stochastic solution of space-time fractional diffusion equations. Physical Review, 2002, E65: 1103-1106. |