IBN AL- HAITHAM J. FO R PURE & APPL. SC I. VO L.24 (1 ) 2011 Explicit Finite Difference Approximation for the Two- Dimensional Fractional Dispersion Equation I. I. Gorial De partment of Mathematics, Ibn Al–Haitham College Education, Baghdad University Recei ved in Feb,7,2010 Accepted in May,13,2010 Abstract In this p aper, we introduce and discuss an algorithm for the numerical solution of two- dimensional fractional disp ersion equation. The algorithm for the numerical solution of this equation is based on exp licit finite difference app roximation. Consist ency , conditional st ability , and convergence of this numerical method are described. Finally, numerical examp le is p resented to show the disp ersion behavior according to the order of the fractional derivative and we demonst rate that our exp licit finite difference app roximation is a comp utationally efficient method for solving two-dimensional fractional disp ersion equation. Key words: Fractional derivative, Two-dimensional probem, Explicit Eul er method, fractional dispersi on e quati on, S tability, Convergence Introduction The sp ace fractional disp ersion equation is obtained from the classical disp ersion equation by replacing the second sp ace derivative by a fractional derivative. Numerical methods associated with integer-order differential equations have been treated extensively in the literature. On other hand, st udies of the numerical methods and error estimates of fractional order differential equations are quite limited to date [1, 2, 3]. M any works by researchers from various fields of science and engineering deal with dy namical sy st ems described by fractional p artial differential equations, which have been used to represent many natural processes in phy sics [4], finance [5,6], and hydrology [7,8]. In this p aper, we find the numerical solution of the two-dimensional fractional disp ersion equation of the form: ),,( ),,( ),( ),,( ),( ),,( tyxq y tyxu yxb x tyxu yxa t tyxu              …………(1) subject to the initial condition u (x,y ,0) = f(x,y ), for x0  x  xR and y 0  y y R …………………(2) and the boundary conditions u (x0,y ,t) = 0, for y 0  y y R and 0 t  Τ u (x,y 0,t) = 0, for x0  x  xR and 0 t  Τ …………(3) u (xR,y ,t) = g(y ,t), for y 0  y y R and 0 t  Τ u (x,y R,t) = k(x,t), for x0  x  xR and 0 t  Τ IBN AL- HAITHAM J. FO R PURE & APPL. SC I. VO L.24 (1 ) 2011 where a, b and f are known functions of x and y, g is a known function of y and t, k is aknwon function of x and t.  and  are given fractional number. q is a knwon function of x, y and t. We use a variation on the classical exp licit Euler method. We p rove this method by using a novel shifted version of the usual Grunwaled finite difference an app roximation for the non-local fractional derivative op erator. Explicit Finite Difference Approximation for Solving the Two-Dimensional Fractional Dispersion Equation In this section, we p rop ose exp licit finite difference app roximation for solving the initial and boundary value problem two-dimensional fractional disp ersion equation (1)-(3). The finite difference method st arts by dividing the x-interval [x0, xR] into n subintervals to get the grid p oints xi = x0 + ix, where   nxxx R 0 and i=0,1,…,n. Also we divide the y- interval [y0, yR] into m subintervals to get the grid p oints yj = y0 + jy, where   myyy R 0 and j=0,1,…,m. Also, t he t-interval [0,T] is divided into M subintervals to get t he grid p oints ts = st, s = 0, …,M, where MTt  . Now, we evaluate eq.(1) at (x i, y j, t s) and we use the exp licit finite difference app roximation to get ),,( ),,( ),( ),,( ),( ),,(),,( 1 sji sji ji sji ji sjisji tyxq y tyxu yxb x tyxu yxa t tyxutyxu              …..(4) Then use the shifted Grunwald est imate to the  , - the fractional derivative, [9]: )(),,)1(( )( 1),,( 0 , xOtyxkxug xx tyxu M k k          …………….(5) )(),)1(,( )( 1),,( 0 , yOtykyxug yy tyxu M k k          to reduce eq.(4) as in the following form ),,( 1 ),( 1 ),( 1 0 1,,,1 1 0 , , 1 , sji j k s kjikji s jki i k kji s ji s ji tyxqug y yxbug x yxa t uu                          s ji s kji j k k ji i k s jkik ji s ji s ji qug y b ug x a t uu ,1, 1 0 , , 1 0 ,1, ,, 1 ,                 , 1,...,1  ni , Msmj ,...,0,1,...,1  ……………. (6) where ),,(, sji s ji tyxuu  , ),(, jiji yxaa  , ),(, jiji yxbb  , ),,(, sji s ji tyxqq  , ! )1()1( )1(, k k g k k      , k=0,1,2,… and ! )1()1( )1(, k k g k k      , k=0,1,2,… The resulting equation can be exp licitly solved for 1 , s jiu to give s ji s ji s kji j k kji i k s jkikji s ji utqug y t bug x t au ,,1, 1 0 ,, 1 0 ,1,, 1 ,                 ...............(7) Also form the initial condition and boundary conditions one can get IBN AL- HAITHAM J. FO R PURE & APPL. SC I. VO L.24 (1 ) 2011 jiji fu , 0 ,  , i=0,…, n 0,0  s ju , j=0,…, m and s=1,…,M 00,  s iu , i=0,…, n and s=1,…,M sj s jR gu , , j=0,…, m and s=1,…,M si s Ri ku , , i=0,…, n and s=1,…,M where ),,(, sjiji tyxff  , ),( sj s j tygg  and ),( si s i txk  Stability Analysis of the Explicit Finite Difference Approximation Define the following fractional p artial difference op erators:       1 0 ,1, , ,, i k s jkik jis jix ug x a u  and s kji j k k jis jiy ug y b u 1, 1 0 , , ,,        which is an )( xO  app roximation to the  th fractional derivative and )( yO  app roximation to the  th fractional derivative term. Then eq.(7) may be writt en in the operator form s ji s jiyx s ji tquttu ,,,, 1 , )1(     ................(8) eq.(8) may be writt en in form s ji s jiyx s ji tquttu ,,,, 1 , )1)(1(     ..................(9) where Ts mn s m s m s n sss n sss uuuuuuuuuU ],,,,,,,,,,,,[ 1,11,21,12,12,22,11,11,21,1   To solve the p roblem for each fixed yj to obtain an intermediate solution  jiu , from 1,,,, )1(   s ji s jijix utqut  ......................(10) Then solve for each fixed xi  ji s jiy uut ,,, )1(  .......................(11) Now, we must p rove each one-dimensional exp licit sy st em defined by the linear difference eqs. (10) and (11) is conditionally st able for all 1 < < 2, 1 <  < 2. The orem: The exp licit sy st em defined by the linear difference eqs.(10) and (11) with 1< < 2, 1 <  < 2 is conditionally st able if max 1 ax t     and max 1 by t     Proof: At each grid p oint y k, for 1,,1  mk  , the sy st em of equation defined by eq.(10) can be written in the exp licit matrix form skkk s k QtUCU  1 where IBN AL- HAITHAM J. FO R PURE & APPL. SC I. VO L.24 (1 ) 2011 ,],,,[ 1,1 1 ,2 1 ,1 1 Ts kn s k s k s k uuuU      ,],,,[ ,1,2,1 T knkkk uuuU      Ts kn s k s k s k qqqQ ],,,[ ,1,2,1   kC is the matrix of coefficients, and is the sum of a lower triangular matrix and a sup er diagonal matrix at t he grid p oint y k, where the matrix entries along the ith row are defined from eq.(10). For examp le, for i = 1 the equation becomes s kkkkkkk s k tqugugugu ,1,20,,1,11,,1,02,,1 1 ,1 )1(     for i = 2 we have s kkkkkkkkk s k tqugugugugu ,2,30,,2,21,,2,12,,2,03,,2 1 ,2 )1(     and for i = n - 1 we get          knknknknknkn s kn ugugugu ,11,,1,22,,1,0,,1 1 ,1 )1(    s knknkn tqug ,1,0,,1     Where the coefficients   x t a jiki    ,, , Therefore the resulting matrix entries jiC , for 1,,1  ni  and 1,,1  nj  by are defined by          1,, 0,, 2,, 1,, , 1 jki ki ki ki ji g g g g C         for for for for 1 1 1     ij ij ij ij According to the Greshgorin theorem [9], the eigenvalues of the matrix C lie in the union of the circles centered at iic , with radius     n il l lii cr 0 , . Here we have   kikiii gc ,1,,, 11  and          n il l i il l kilikilii gcr 0 1 0 ,1,,,   and therefore 1,  iii rc . We also have  kikikiiii rc ,,,, 211             x t a ki,21         x t a max 21 Therefore, for the sp ectral radius of t he matrix C to be at most one, it suffices t o have 121 max           x t a          1max  x t a      1 max   x t a max 1 ax t     Same method above, resulting the sy st em of equation defined by eq.(11) is t hen defined by , k s kk UUS where ,],,,[ 1,2,1, Ts mk s k s k s k uuuU   ,],,,[ 1,2,1, T mkkkk uuuU      IBN AL- HAITHAM J. FO R PURE & APPL. SC I. VO L.24 (1 ) 2011 kS is the matrix of coefficients, and is the sum of a lower triangular matrix and a sup er diagonal matrix at the grid p oint xk for 1,,1  nk  . Therefore the resulting matrix entries kS for 1,,2,1  mi  and 1,,1  mj  by are defined by          1,, 0,, 2,, 1,, , 1 jik ik ik ik ji g g g g S         for for for for 1 1 1     ij ij ij ij where the coefficients   y t b jiki    ,, So, and in the same way , according to the Greshgorin theorem [9], to get max 1 by t     Consiste ncy and Convergent Analysis of the Explicit Finite Difference Approximation We note that the three difference op erators used in eq.(6) each have a local truncation error with ),( tO  ),( xO  and )( yO  resp ectively. The )( tO  , for the time derivative term, is obtained from the classical Taylor's exp ansion. The )( xO  and )( yO  for the local truncation error of the fractional derivative terms was p roved in [10]. Therefore, the exp licit finite difference app roximation is consistency . Theorem above shows that s jiyx u ,,, )(   converges to the mixed fractional derivative linearly, as )()( yOxO  . Therefore, the local truncation error of the exp licit Euler method eq.(8) is )()()( yOxOtO  . This consistency of the exp licit finite difference app roximation together with the above result on conditional stability imp lies that t he exp licit finite difference app roximation is convergent and this convergence is )( tyxO  . Numerical Example In this section, numerical examp le is p resented which confirm our theoretical results. Example : Consider the two-dimensional fractional disp ersion equation: tt exyeyx y tyxuy x tyxu x t tyxu 22225.0 6.1 6.16.1 5.1 5.1 ),,( 2 )4.1(),,( )5.0( ),,(          subject to the initial condition u (x,y ,0) = xy 2 , 0  x  0.5, 0 < y < 0.5 and the boundary conditions u (0,y ,t) = 0, 0 < y < 0.5, 0  t 0.025 IB N AL- HAITHAM J. FOR PURE & APPL. SC I. VO L.24 (1) 2011 u (x,0,t) = 0, 0 < x < 0.5, 0  t 0.025 u (,0.5,y ,t) = 0.5 e 2t , 0 < y < 0.5, 0  t 0.025 u (x,, 0.5,t) = 0.25e 2t x , 0 < x < 0.5, 0  t 0.025 This fractional disp ersion equation together with the above initial and boundary condition is constructed such that t he exact solution is u(x,y ,t) = e 2t xy 2 . Table (1) and (2) show the numerical solution using the exp licit finite difference app roximation. From table (1) and (2), it can be seen that thereisa good agreement between the numerical solution and exact solution. 6. Conclusi ons In this p aper, a numerical method for solving the two-dimensional fractional disp ersion equation has been described and demonst rated. The exp licit difference app roximation is p roved to be conditionally st able and converges. Furt hermore numerical examp le is p resented to show that good agreement between the numerical solution and exact solution has been noted. Re ferences 1. Huang, F. and Liu, F. (2005), Anziam Journal, 46: 1-14. 2. Shen, S. and Lin,F. (2005),ANZIAM J., 46(E): 871-887. 3. M eerschaert, M .M .; Scheffler, H.P. and Tadjeran,C. (2006), J. Comp ut. Phy s., 211:249–261. 4. M etz ler, R. and Klafter, J. (2004), J. Phy sics, A 37: R161-R208. 5. Sabatelli, L.; Keating, S.; Dudley,J. and Richmond, P. (2002), Eur. Phy s. J., B 27: 273–275. 6. Scalas,E.; Gorenflo, R. and M ainardi, F. (2000), Phy s., A 284: 376–384. 7. Benson, D.; Wheatcraft, S. and M eerschaert, M . (2000), Water Resour. Res., 36:1403–1412. 8. Schumer,R.; Benson, D.A.; M eerschaert, M .M . and Baeumer, B. (2003), Water Resour. Res., 39: 1022–1032. 9. Isaacson, E. and Keller, H.B. (1966), Wiley, New York. 10.M eerschaert, M .M . and Tadjeran,C. (2004), J. Comput. Ap p l. M ath., 172: 65–77. Table (1) The numerical sol uti on of e xample by using the explicit fi ni te difference approximati on for 1.0,1.0  yx and 0125.0t Numerical S oluti on Exact S oluti on Error 9.730E-4 1.02532 E -3 6.72814 E -2 7.876E-3 8.20252 E -3 3.26521 E -4 0.027 2.76835 E -2 6.83508 E -4 0.064 6.56202 E -2 1.62017 E -3 9.795E-4 1.05127 E -3 7.17711 E -5 7.759E-3 8.41017 E -3 6.51169 E -4 0.027 2.83843 E -2 1.38432 E 3 0.036 6.72814 E -2 3.12814 E -2 Table (2) The numerical sol uti on of e xample by using the explicit fi ni te difference approximati on for 125.0,125.0  yx and 0125.0t Numerical S oluti on Exact S oluti on Error 1.908E-3 2.00257 E -3 9.45686 E -5 0.015 1.60205 E -2 1.02055 E -3 0.052 5.40693 E -2 2.06935 E -3 1.929E-3 2.05326 E -3 1.24264 E -4 0.015 1.64261 E -2 1.42611 E -3 0.036 5.54381 E -2 1.94381 E -2 2011) 1( 24مجلة ابن الھیثم للعلوم الصرفة والتطبیقیة المجلد تتشتال دلةلمعا الفروق المنتهیة الصریحةّ تقریب ذات البعدین الكسریة ایمان ایشو كولایر جامعة بغداد ،ابن الهیثم -كلیة التربیة ،قسم الریاضیات 2010 شباط 7استلم البحث في 2010 ایار 13قبل البحث في الخالصة وارزمیــة الحـل العــددي وان خ. البعـدین ذي الكســریة تتشـتال ةلمعادلــفـي هـذا البحــث قـدمنا وناقشــنا خوارزمیـة للحــل العـددي شـرطي، والتقـارب للطریقـة السـتقرار االو ،تسـاق اال اكمـا ناقشـن. الصـریحة الفـروق المنتهیـة تقریـب قائمة على اساس دلةالتلك المع .العددیة وبینـا أن تقریـب الفـروق المنتهیـة الصـریحة ةحسـب مرتبـة االشـتقاق الكسـری التشـتتسلوك ظهرلی اعددی اخیرا قدمنا مثاال .البعدین الكسریة ذي التشتت هو طریقة فعالة حسابیا لمعادلة .التقارب االستقرار، التشتت الكسریةمعادلة ، طریقة اوبلر الصریحة، مسألة ذي بعدین، كسریة مشتقة -:الكلمات المفتاحیة