1、 基金项目：国家自然科学基金项目（51174236） ；国家重点基础研究发展计划（2011CB606306）; 中南大学研究生自由探索项目（2013zzts146）通信作者：郑洲顺，教授，博士，中南大学数学与统计学院，E-mail: Riesz 空间分数阶扩散方程的分数阶中心差分加权离散格式邓娟，郑洲顺 *中南大学数学与统计学院，湖南长沙，410083摘要：在有限区域内考虑带齐次 Dirichlet 边界条件的 Riesz 空间分数阶扩散方程的初边值问题，利用分数阶中心差分对空间方向进行离散，在时间方向上用隐式和显式 Euler 格式的加权平均进行离散，构造了空间 2 阶、时间 阶( )的全离
2、散加权差分格式。利用函1,2数的单调性证明了当加权因子 时差分离散格式是无条件稳定的，当 时差分120 12离散格式是条件稳定的，并给出了稳定的条件。证明了相应差分离散格式的收敛性。用实际数值算例验证了差分离散格式的有效性。关键词：Riesz 导数；分数阶扩散方程；分数阶中心差分；稳定性分析；收敛性分析中图分类号：O 241.82 文献标识码：A分数阶微分方程越来越多地出现在不同的科学领域的应用中，例如分数阶扩散方程、波动方程、薛定谔方程、电报方程、弛豫振荡方程等等 1-5，应用领域涵盖了流体力学 3,6、地下水模拟 7、 、湍流 8、生物数学、统计力学、分形介质中的扩散问题等。由此，分数阶微
3、分方程求解也受到越来越多的关注，尤其是数值方法。以往的差分方法大多数是基于Grunwald-Letnikov 导数得到的逼近格式 9-13，或基于数值积分的思想 13-15，或利用Richardson 外推技巧 10,16等方法来得到的。作为一个新的方法，Ortigueira 17给出了分数阶中心差分（Fractional Centred Difference）的定义，证明了解析函数的 Riesz 导数可以由分数阶中心差分来表示。C.Celik 18等用分数阶中心差分的定义构造了一种二阶精度的 Crank-Nicholson 差分算法。本文将利用分数阶中心差分构造 Riesz 空间分数阶扩散方
4、程的加权平均差分格式，分析相应格式的稳定性和收敛性，用数值算例验证算法的可行性。这里考虑如下的带齐次 Dirichlet 边界条件的 Riesz 空间分数阶扩散方程的初边值问题（1）(,)(,), in,(0,uxtuxtDptabT初边值条件为： ， 。(,),0atbt0(,)(,xux其中， 为扩散系数， 为源项，分数阶导数为 Riesz 分数阶导数，即0D(,)pxt。21()1(),22cos()(f xfdx 1 差分格式的构造Manuel D. Ortigueira 在文献17中给出了分数阶中心差分的如下定义， （2）(1)() ()12khkfx fxkh其中， ，并给出了分数
5、阶中心差分与分数阶导数的逼近关系如下：12。00() ()lim()lim()12khhhkdfxfx fxkh 为了方便起见，记 ，则上述逼近关系即为(1)2kg。 （3）0()lim()khdfxgfxh在有限区域 内，对于给定的均匀网格剖分的空间步长 和时间步长 ，有 ，,ab baMh从而有（4）(,)(,)jMj hkjxuthguxktR其中， 为截断误差。hR在文献18中，Cem elik 对 Riesz 分数阶导数的中心差分逼近格式（4）的精度进行了分析，证明了其逼近精度为二阶，为叙述方便将其结论作为引理如下。引理 1 令 ，且 ，令 为分数阶中心差分，5()fCR()1,5n
6、fLR()()hkfxgfxh则当 时，有0h2()()hfxfOh其中， 为 Riesz 导数，且 。()fx12一般地，对扩散方程（1）的时间导数差分离散的隐式和显式 Euler 格式分别为， ，1 1nMjjj nkjjjuDhgup1nMjjj nkjjjuDhgup根据上述隐式和显式 Euler 格式，用加权平均的方法给出其如下差分格式：。 （5）1 1 1() ()nMj Mjjj nnnnkj kjjjj juhguhgup 其中， 为加权因子。0对带齐次 Dirichlet 边界条件的 Riesz 空间分数阶扩散方程的初边值问题，用分数阶中心差分离散 Riesz 空间分数阶导数
7、，再将时间导数用隐式和显式 Euler 格式的加权平均的进行离散，从而得到其全离散格式为：（6）1 10() ,1,1,2,1nMj Mjjj nnkj kjjj jjjnMuDhguDhgupMnxn 其中， 。()nj jjpp2 误差分析由 的表达式，容易得出下面的结论 18。kg引理 2 关于系数 有如下性质：k（i） ；02(1)0g（ii） ；, kfor（iii ） ；kg（iv） 。00k在文献18中给出了数值逼近 的系数与生成函数的关0lim()()hkhkfxgfxh系，即 ，其中 。2sin()ikzkzgeR当 时，上式变成 ，根据性质（i）和（ii）中 的正负情况，则
8、有00k kg。于是（iv ）得证。 000kkgg定理 1 Riesz 空间分数阶扩散方程的初边值问题的全离散格式（6）的截断误差阶为，其中，当 时， ，当 时， 。2()Oh1212证：令 为方程（1）的解析解，则差分格式的局部截断误差 为：(,uxt njR，1 11,)(,)()(,)(,) ,jnjnnj kjnkjnjnjntuxtRDhguxtDhguxtptpt 由于 ，且由引理 1，有21 2(,)(,)()jnjnjjuxttuOtt 121()(,)(,)() ()nnnkjnkjjkj juuDhgxtDhguxtDDOhxx ，则由方程（1）有， ，上式可变换为：(,
9、)nnjnjjupxttx1(1)(,)(,)kjnkjnDhgtDhguxt 12() ()nj juuOxx。1 21(1)(,)(,)()n njnjnj jptpxtht t由此，可以得出截断误差为 12 1221()(,)(,) (,)1,)nn n nnj jnjnj j jjjnjnuuuRpxtpxtt t ttpxpxOh 然后将上述微分在 处展开，有,jt2 222(1) ()n nn nnnjj j jj juuuR Oht t tt t221()njOht因此，当 时， ，当 时， 。 2()njR122()njROh根据定理 1，得到所构造的全离散格式的逼近精度在空间
10、方向是二阶的，当 时，12在时间方向是二阶的，当 时，在时间方向是一阶的。123 稳定性及收敛性分析为了方便进行差分格式稳定性和收敛性分析，先将离散格式用矩阵形式表示。由齐次Dirichlet 边界条件，可令 ，将方程（1）中的源项 的离散向量记12(,)nnTMUu (,)pxt为 ，即nP，11111()() nnnnnnnMpPP考虑在节点 处的分数阶算子的离散形式中的和式部分，即(,)jnxt，101kj jjjgugugu 根据区间的剖分，记。012320(1)MMggB 令 ，记矩阵 ，从而得全离散格式（6）的矩阵形式为：DhA。 （7）1() ,12,nnIUIAP令 、 分别为
11、矩阵 A、B 的特征值，则由圆盘定理（Gerschgorins Circle Theorem）AB有 00Bikgrg由引理 2 中关于系数 的性质（iv ） 有kg0k00Bikgrg即 ，则根据 有， 。02BgA002A如果要 Riesz 空间分数阶扩散方程初边值问题的差分格式（7）是稳定的，则需要证明矩阵 的特征值 满足条件1()()PIAI 。1()A定理 2 针对 Riesz 空间分数阶扩散方程初边值问题的差分格式（7） ，当 时，差分102格式（7）是无条件稳定的；当 时，差分格式（7）是条件稳定的，其稳定的条件12是 。01Dgh证：记 ， ，则有()()Af01 2 2()(
12、)0(1)AAAf即 在区间 上是单调递减函数，故有 。()f0,1 (0)ff又因 ，于是有Af11(0) (1)()AAAf f因此，当 时，要 ，0()1()Af只要 ，即 ，亦即 。1()()AfAA2(1)A由于 ，故当 时， 恒成立，于是当 时，恒有002Ag2102(1)02成立。从而得到当 时，差分格式（ 7）是无条件稳定的。1()A,而当 时，由 在区间 上是单调递减函数，故有，12()f1,2，()(1)()AAfff因此，要 ，只要 ，即 。1()()Af1Af 2A又因为 ，所以只要 ，即 ，亦即网格比 时，则有002Ag02Ag0g01Dgh恒成立。从而得到当 时，差
13、分格式（7）是条件稳定的，其稳1()A 1,定的条件是 。 01Dgh综上，当 时，差分格式（7）是无条件稳定的，当 时，差分格式2 12（7）是条件稳定的，其稳定的条件是 。01Dgh下面分析 Riesz 空间分数阶扩散方程初边值问题的差分格式（7）的收敛性。令 为微分方程在网格点 出的精确解，数值解为 ，则误差为 ，njv(,)jnxt njunnjjjevu，则 ，代入差分格式（7） ，由初始误差为 0，即 ，有121(,)nMEe njjjuve 0E， （8）1(I)(I)nnAER其中， 。2()nRCh定理 3 Riesz 空间分数阶扩散方程初边值问题的差分格式（7）的收敛性分析
14、如下：当时，差分格式（7）是无条件收敛的，且 ，其中，当102 2(),1,nECThn时， ，当 时， ；当 ，且 时，差分格式是条件收102120Dg敛的，且 。(),nECThn证：由（8）式，有。11 21(I)(I)()(I)In nAEChA 不妨记 ， ，则1(I)P 12QA。1nnP所以有， ，10E，2PBQ，3 2PQ依次类推，有 。12nnEI根据定理 1，可以得到，当 时，或 且 时，矩阵 的特征值10120DghP满足 ，此时，有 成立。1P因此，。11 1nn nEPIQPPIQ 根据（8）式，可以得到 ，且 ，则2()Ch n，11 22()()nnIChCTh
15、于是得到，当 且 时，有 。 0 h10nE因此，当 时，差分格式（7）是无条件收敛的，且12，其中，当 时， ，当 时， ；当(),nECThn 12102，且 时，差分格式是条件收敛的，且 。 1201Dg (),nECThn4 数值算例在方程（1）中，取 ， ，令函数 ，,0,1abD20()1)ux，12 22222(,)() ()(6()(5) 67Sectfxttxxt 则相应的初边值问题如下（9）2(,)(,)(,)01,0,1, ()(),uxttfxttTttTx容易验证该初边值问题的精确解为 。2(,)1()utx分别取不同的加权系数 、 分数阶导数的阶数 ， 在不同的步长
16、下，用所构造的差分方法对方程进行数值求解。当 、 时，时间步长取 ，空间步长取12.5150，在时间点 时方程的解析解和数值解对比如图 1 所示。150ht图 1 时的数值结果对比（ ）t1.5,0,15hFig.1 the comparison between numerical solution and analytical solution at ( )t.,150,h对不同的加权系数 、分数阶导数的阶数 ，在不同的网格剖分下，用所构造的差分方法对方程（9）进行数值求解，空间方向的误差分析结果如表 13，时间方向的误差分析结果如 45。表 1 差分格式的数值结果的最大绝对误差（ ， ，
17、）1.70N2Tab.1 the Maximum absolute error of the numerical solution ( ， ， )1.7012空间步长 150hh0h最大绝对误差 4.352.396.5收敛阶 1063412表 2 差分格式的数值结果的最大绝对误差（ ， ， ）.5NTab.2 the Maximum absolute error of the numerical solution ( ， ， ).50空间步长 150h1h120h最大绝对误差 4.752.97465.38收敛阶 34表 3 差分格式的数值结果的最大绝对误差（ ， ， ）1.50NTab.3 t
18、he Maximum absolute error of the numerical solution ( ， ， )1.503空间步长 150hh120h最大绝对误差 4.625.97025.896收敛阶 1467表 1-3 的数值结果表明，在时间步长固定，且满足稳定性条件和收敛性条件的情况下，空间方向的收敛阶约等于 2，与上述的理论分析结果是相符合的。表 4 差分格式的数值结果的最大绝对误差（ ， ， ）1.50M12Tab.4 the Maximum absolute error of the numerical solution ( ， ， ).5012时间步长 12486最大绝对误差
19、 59.0353.7106.193.8940收敛阶 1.50962.35601.64表 5 差分格式的数值结果的最大绝对误差（ ， ， ）.1M4Tab.5 the Maximum absolute error of the numerical solution ( ， ， ).5104时间步长 124816最大绝对误差 46.033.8104.3759.80收敛阶 92092表 4-5 的结果表明，在空间步长固定，且满足稳定性条件和收敛性条件的情况下，当时，时间方向的收敛阶接近于 1，当 ，时间方向的收敛阶也比较接近 2，与理论12 2分析结果基本相符。5 本文的主要结论（1）本文针对带齐次
20、 Dirichlet 边界条件的 Riesz 空间分数阶扩散方程的初边值问题，在空间方向上采用分数阶中心差分格式进行离散，将隐式和显示 Euler 格式进行加权平均，对时间进行离散，构造了 Riesz 空间分数阶扩散方程的全离散差分格式，其截断误差在空间方向是二阶的，在时间方向是 阶，其中，当 时， ，当 时， 。1212（2）根据圆盘定理，利用函数的单调性，证明了当 时所构造的差分离散格式0是无条件稳定的，而当 时，所构造的差分离散格式是条件稳定的，其稳定的条件12是 。01Dgh（3）针对给定的剖分，将方程的全离散差分格式用矩阵形式表示，给出了误差传递公式，分析了所构造的差分离散格式的收敛
21、性。（4）针对实际的算例，用所构造的差分方法进行了数值求解，通过数值解和解析解的比较，结果表明差分方法是可行的。参考文献1W. Deng. Finite element method for the space and time fractional Fokker-Planck equationJ. SIAM Journal on Numerical Analysis, 2008,47(1):204-226.2胡劲松，胡朝浪，郑茂渡. 广义对称正则长波方程的一个拟紧致守恒差分算法 J. 四川大学学报：自然科学版，2010,47(2):221.3S. Yao, Y. Zhang. Numeric
22、al solution of generlized Maxwell fluid fractional differential equationJ. Computer Engineering and Applications, 2013,49(12):41-46.4杨晨航,刘发旺. 分数阶Relaxation-Oscillation 方程的一种分数阶预估-校正方法J.厦门大学学报（自然科学版）,2005,44(6):761-765.5章红梅,刘发旺.时间分数阶电报方程的一种解技巧J. 厦门大学学报（自然科学版） ,2007,46(1):10-13.6L. W. Gelhar, C. Wehy,
23、 Kenneth R. Rehfeldth. A critical review of data on field-scale dispersion in aquifersJ. Water Resources Research, 1992,28(7):1955-1974.7 M. Xu，W. Tan. Theory analysis of general second order anomalous diffusion velocity field, stress field and eddy field in liquidJ. Science in China(Series A), 2001
24、,31:626-637.8D.del-Castillo-Negrete, B. A. Carreras, V. E. Lynch. Fractional diffusion in plasma turbulenceJ. Physics of Plasmas, 2004,11(8):3854-3864. 9M. M. Meerschaert, Charles Tadjeran. Finite difference approximations for two-sided space- fractional partial differential equationsJ. Applied Nume
25、rical Mathematics, 2006,56:80-90.10C. Tadjeran, M. M. Meerschaert, H. P. Sehefller. A second-order accurate numerical approximation for the fractional diffusion equationJJournal of Computational Physics, 2006, 213:205-213.11S. Shen, F. Liu, V. Anh, Numerical approximations and solution techniques fo
26、r the space-time Riesz-Caputo fractional advection-diffusion equationJ, Numerical Algorithms 56(2011) 383-40312林然,刘发旺.分数阶常微分方程初值问题的高阶近似J.厦门大学学报(自然科学版),2004,43(1):21-25. 13Q. Yang, F. Liu, I Turner, Numerical method for fractional partial differential equations with Riesz space fractional derivatives
27、J, Applied Mathematical Modeling34(1)(2010)200-21814Q. Yang, I turner, F. Liu, Analytical and numerical solutions for the time and space- symmetric fractional diffusion equationJ, Australian and New Zealand Industrial and Applied Mathematics Journal(ANIZIAM) ,50(2009):800-814.15于强,刘发旺.时间分数阶反应- 扩散方程的
28、隐式差分近似J.厦门大学学报（自然科学版）,2006,45(3):315-319.16蔡新,刘发旺.空间分数阶扩散方程的隐式高精度方法J.厦门大学学报（自然科学版）,2007,46(3):317-321.17M. D. Ortigueira, Riesz potential operators and inverses via fractional centred derivativesJ. International Journal of Mathematics and Mathematical Sciences(2006)1-12.18C.Celik, M. Duman. Crank-Ni
29、colson method for the fractional diffusion equation with the Riesz fractional derivativeJ. Journal of Computational Physics 231(2012)1743-1750.A Weighted Fractional Centred Difference Method for the Space Fractional Diffusion Equation based on the Riesz DerivativesJuan DENG, Zhoushun ZHENG *School o
30、f Mathematics and Statistics, Central South University, Changsha, 410083, ChinaAbstract: In this paper, the space fractional diffusion equation with homogeneous Dirichlet boundary condition based on Riesz derivative on a finite domain is considered. A Weighted difference method with the numerical er
31、ror is derived, when the fractional 2()12)Ohorcentred difference is applied to obtain the finite difference approach in space and the weighted average method based on the implicit and explicit Euler method is proposed in time direction. The stability analysis of this scheme is carried out by means o
32、f the monotonicity of function. It is proved that the method is unconditionally stable for , and conditionally stable for 120, here is the weighting parameter subjected to . The convergence of the 12corresponding difference method is proved. A numerical example with known exact solution is examined to verify the effectiveness of the method.Keywords: Riesz fractional derivative; Fractional diffusion equation; Fractional centred difference; Stability analysis; convergence analysis
