A new difference scheme for fractional cable ... - Semantic Scholar

1 downloads 0 Views 401KB Size Report
Ibrahim Karatay 1, Nurdane Kale 2∗. 1Department of Mathematics, Fatih University, Istanbul, Turkey. 2Vocational School, Fatih University, Istanbul, Turkey.
International Journal of Applied Mathematical Research, 4 (1) (2015) 52-57 www.sciencepubco.com/index.php/IJAMR c

Science Publishing Corporation doi: 10.14419/ijamr.v4i1.3875 Research Paper

A new difference scheme for fractional cable equation and stability analysis Ibrahim Karatay 1 , Nurdane Kale 1

2∗

Department of Mathematics, Fatih University, Istanbul, Turkey 2 Vocational School, Fatih University, Istanbul, Turkey *Corresponding author E-mail: [email protected]

c Copyright 2015 Ibrahim Karatay and Nurdane Kale. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract We consider the fractional cable equation. For solution of fractional Cable equation involving Caputo fractional derivative, a new difference scheme is constructed based on Crank Nicholson difference scheme. We prove that the proposed method is unconditionally stable by using spectral stability technique. Keywords: Cable equation; Caputo fractional derivative; Difference scheme; Stability.

1.

Introduction In this study, we consider the following time fractional cable equation;

 

∂ α u(x,t) ∂tα

2

u(x,t) = ∂ ∂x − µ2 u (x, t) + f (x, t) , (0 < x < 1, 0 < t < 1), 2 u(x, 0) = r(x), 0 < x < 1,  u(0, t) = 0, u(1, t) = 0, 0 ≤ t ≤ 1.

Here, the term

∂ α u(t,x) ∂tα

∂ α u(x, t) 1 = ∂tα Γ(1 − α)

Zt

(1)

denotes α-order Caputo derivative with the formula: ut (x, τ ) dτ, where 0 < α < 1, (t − τ )α

(2)

0

where Γ(.) is the Gamma function.

2.

Discretization of problem

We introduce the basic ideas for the numerical solution of the Time Fractional Cable equation by CrankNicholson difference scheme. For some positive integers M and N , the grid sizes in space and time for the finite difference algorithm are defined by h = 1/M and τ = 1/N, respectively. The grid points in the space interval [0, 1] are the numbers

53

International Journal of Applied Mathematical Research

xj = jh, j = 0, 1, 2, ..., M , and the grid points in the time interval [0, 1] are labeled tk = kτ, k = 0, 1, 2, ..., N . The values of the functions u and f at the grid points are denoted ukj = u(xj , tk ) and fjk = f (xj , tk ), respectively. Let u(x, t), ut (x, t) and utt (x, t) are continuous on [0, 1]. As in the classical Crank-Nicholson difference scheme, a discrete approximation to the fractional derivative ∂ α u(x,t) at (xj , tk+ 12 ) can be obtained by the following approximation[12]: ∂tα ∂ α u(xj , tk+ 12 ) ∂tα

"

k−1 X

(uk+1 − ukj ) j = w1 u + (wk−m+1 − wk−m ) u − wk u + σ 21−α m=1 k

m

#

0

(3)

+ O(τ 2−α ).  1 1 1−α Where σ = Γ(2−α) − (j − 1/2)1−α In addition for k = 0 there is no these terms τ α and wj = σ (j + 1/2) w1 uk and wk u0 . On the other hand, we have ∂ 2 u(xj , tk+ 21 ) ∂x2

3.

" k+1 # k+1 k+1 ukj+1 − 2ukj + ukj−1 1 uj+1 − 2uj + uj−1 = + + O(h2 ). 2 h2 h2

(4)

The proposed difference scheme

Using these approximations (3) and (4) into (1), we obtain the following difference scheme for (1) which is accurate of order O(τ 2−α + h2 );

k−1 P (uk+1 −uk ) w1 uk + (wk−m+1 − wk−m ) um − wk u0 + σ j21−α j =  km=1k+1  uj +uj −µ2 + f (xj , tk + τ2 ) 2

1 2



k+1 +uk+1 uk+1 j−1 j+1 −2uj h2

   k−1 P −uk (uk+1  j) j k m 0  w1 uj + (wk−m+1 − wk−m ) uj − wk uj + σ 21−α      k+1 m=1   k k+1   k+1 k k  uj+1 −2uj +uk+1 uj +uj uk  j−1 j+1 −2uj +uj−1 2  + +µ = f (xj , tk + τ2 ),  − 2h2 2h2 2 0 ≤ k ≤ N − 1, 1 ≤ j ≤ M − 1,         u0 = r(xj ), 1 ≤ j ≤ M,    kj u0 = 0, ukM = 0, 0 ≤ k ≤ N.

  h  k+1  σ  k+1 i µ2 k+1 1 1 1  u + u + − uj−1 −  2 1−α + h2 + 2 2 j+1 j 2h 2 2h  h     k  k i  2  µ 1 σ 1 1 k  u + − + + u + − uj−1 + −  2 1−α 2 2 j+1 j 2h 2 h 2    2h   k−1 P 0 + w1 ukj + (wk−m+1 − wk−m ) um j − wk uj  m=1    = f (xj , tk + τ2 ), 0 ≤ k ≤ N − 1, 1 ≤ j ≤ M − 1,    0  u = r(xj ), 1 ≤ j ≤ M,    kj u0 = 0, ukM = 0, 0 ≤ k ≤ N.

We can arrange the system above to obtain

+

k k uk j+1 −2uj +uj−1 h2



54

International Journal of Applied Mathematical Research

  i  h   µ2 σ 1  uk+1 + − 2h1 2 uk+1 − 2h1 2 uk+1  j j+1 + 21−α + h2 + 2 j−1   h    k  k i   µ2 1 σ 1 1 k  u + − + u + − + + uj−1 −  2 1−α 2 2 j j+1 2h 2 h 2 2h     P k  + wm uk−m+1 − uk−m j j  m=1   τ  = f (xj , tk + 2 ), 0 ≤ k ≤ N − 1, 1 ≤ j ≤ M − 1,    0  u = r(xj ), 1 ≤ j ≤ M,    kj u0 = 0, ukM = 0, 0 ≤ k ≤ N. Then we rewrite the equation following type      − 2h1 2 ukj+1 + − 2h1 2 uk+1  j+1        k   P  µ2 µ2 k+1 k−m+1 k−m σ 1 σ 1 k  + − + + u + + + u + w u − u  1−α 2 1−α 2 m j j j j  2 h 2 2 h 2   m=1  k  k+1  1 1 + − 2h2 uj−1 + − 2h2 uj−1    = f (xj , tk + τ2 ), 0 ≤ k ≤ N − 1, 1 ≤ j ≤ M − 1,     u0 = r(xj ), 1 ≤ j ≤ M,   kj u0 = 0, ukM = 0, 0 ≤ k ≤ N.

3.1.

(5)

Spectral Stability of the Difference Method

The difference scheme above (5) can be written in matrix form:  T 0 DUj+1 + EUj + DUj−1 = ϕj where ϕj = ϕ0j , ϕ1j , ϕ2j , ..., ϕN , ϕj = r(xj ), ϕkj = f (xj , tk+ 12 ), 1 ≤ j ≤ M, 1 ≤ j   0 1 2 T k ≤ N, and Uj = UJ , UJ , UJ , ..., UJN . Here D(N +1)x(N +1) and E(N +1)x(N +1) are the matrices of the form 



0  1     1 D= − 2  2h    

     E=   

1 1

1 .. .

..

.

..

.

        

..

. 1

1



1 b −w1 −w2 .. .

a b + w1 w2 − w1

−wN −1

wN −1 − wN −2 2

a b + w1 .. . ···

a .. . w2 − w1

       

..

. b + w1

a

2

σ σ where a = 21−α + h12 + µ2 , b = − 21−α + h12 + µ2 Using the idea on the modified Gauss-Elimination method, we can convert into the following form: Uj = ψj+1 Uj+1 + µj+1 , j = M, ..., 2, 1, 0.

Then, we write D + Eψj+1 + Dψj ψj+1 = 0, Eµj+1 + Dψj µj+1 + Dµj = ϕj , where 1 ≤ j ≤ M. So, we obtain the following pair of formulas: −1 −1 ψj+1 = − (E + Dψj ) D, µj+1 = (E + Dψj ) (ϕj − Dµj ) , where 1 ≤ j ≤ M.

55

International Journal of Applied Mathematical Research

We will prove that ρ (ψj ) < 1, 1 ≤ j ≤ M , by induction. Since ψ1 is a zero matrix ρ (ψ1 ) = 0 < 1.Moreover,   2 ψ2 = −E −1 D, ρ (ψ2 ) = ρ −E −1 D = σ −11 µ2 . − 2h1 2 =  σ 1/h 1 µ2  ,since ψ2 is of the form + + + 2 2 1−α 2+ 1−α 2 h

2

h

2



2



0   ∗    ∗ ψ2 =       ∗

1/h  2

σ 21−α

2 2

+ h12 + µ2

 1/h2



 2

σ 21−α

2

+ h12 + µ2



.. ∗



.



1/h2 2



σ 21−α

2

+ h12 + µ2



           (M +1)x(M +1)

1 1 σ = Γ(2−α) τ α > 0, therefore, ρ (ψ2 ) < 1. Now, assume ρ (ψj ) < 1. We find that −1 ψj+1 = − (E + Dψj ) D

0  ∗    ∗ 1  = 2 2h    ∗ 

 1 E2,2 −(1/2h2 )ψj2,2

1 E3,3 −(1/2h2 )ψj3,3



.. ∗

and we already know that Ej,j = ρ (ψj+1 ) =

. ∗

∗ σ 21−α

+

1 h2

+

σ 21−α

+

1 h2

+

and ψjr,r = ρ (ψj ) for 2 ≤ r ≤ M + 1 :

 = h − 2h1 2 ρ (ψj ) 2 M 2 1 −

1/2h2 µ2 2

µ2 2

1 EM +1,M +1 −(1/2h2 )ψjM +1,M +1

      

M2 

ρ(ψj ) 2

+

σ 21−α

+

µ2 2

i

Since 0 ≤ ρ (ψj ) < 1, it follows that ρ (ψj+1 ) < 1. So, ρ (ψj ) < 1 for any j, where 1 ≤ j ≤ M.

4.

Numerical example

Consider this problem, ∂ α u(t,x) ∂tα

2

2−α

u(t,x) 2t = ∂ ∂x − u(t, x) + Γ(3−α) (1 − x) sin(x) + 2t2 [cos(x) + (1 − x) sin(x))] , 2 (0 < x < 1, 0 < t < 1),  u(0, x) = 0, 0 ≤ x ≤ 1,   u(t, 0) = 0, u(t, 1) = 0, 0 ≤ t ≤ 1.

   

Exact solution of this problem is u(t, x) = t2 (1 − x) sin(x). The errors for some M and N are given in figure 1. The errors when solving this problem are listed in the table1 for various values of time and space nodes. Table 1: The errors for some values of M, N and α

N 8 16 32

M 32 32 32

α = 0.3 Error(α, τ ) 0.001811212 0.000449950 0.000111687

Err. rate 4.02 4.02

α = 0.5 Error(α, τ ) 0.001688126 0.000409875 0.000099150

Err. rate 4.1 4.1

α = 0.8 Error(α, τ ) 0.001265365 0.000301407 0.000086960

Err. rate 4.1 3.46

56

International Journal of Applied Mathematical Research

Figure 1: The errors when t=1 for some M and N

5.

Conclusion

In this work, O(τ 2−α + h2 ) order approximation for the Caputo fractional derivative based on the CrankNicholson difference scheme was successfully applied to solve the time-fractional cable equation. It is proven that the time-fractional Crank-Nicholson difference scheme is unconditionally stable by spectral stability analysis.

References [1] Xikui Li, Xianhong Han, Xuanping Wang, ”Numerical modeling of viscoelastic flows using equal low-order finite elements”, Comput. Methods Appl. Mech. Engrg., Vol.199, (2010), pp.570-581. [2] M. Raberto, E. Scalas, F. Mainardi, ”Waiting-times returns in high frequency financial data: an empirical study”,Physica A, Vol.314, (2002), pp.749-755. [3] D.A. Benson, S. Wheatcraft, M.M. Meerschaert,”Application of a fractional advection–dispersion equation”, Water Resour. Res., Vol.36, (2000), pp.1403–1412. [4] X. Li, M. Xu, X. Jiang,”Homotopy perturbation method to time-fractional diffusion equation with a moving boundary”,Appl. Math. Comput. , Vol. 208, (2009), pp.434–439. [5] J. A. T. Machado, No.1,(2001),pp.47–66.

”Discrete-time fractional-order controllers”,

Fractional Calculus Applied Analysis,Vol.4,

[6] Z. Deng, V.P. Singh, L. Bengtsson,” Numerical solution of fractional advection-dispersion equation”,J. Hydraulic Eng. Vol.130, (2004), pp. 422–431. [7] V.E. Lynch, B.A. Carreras, D. del-Castillo-Negrete, K.M. Ferreira-Mejias, H.R. Hicks,” Numerical methods for the solution of partial differential equations of fractional order”,J. Comput. Phys. , Vol. 192, (2003), pp. 406–421. [8] I. Podlubny, Fractional Differential Equations, Academic Press, New York, (1999). [9] Chang-Ming Chen , F. Liu , I. Turner b, V. Anh,” A Fourier method for the fractional diffusion equation describing sub-diffusion”, Journal of Computational Physics, Vol. 227, (2007), pp. 886-897. [10] Ibrahim Karatay, S ¸ erife Rabia Bayramo˘ glu,” An Efficient Difference Scheme for Time Fractional Advection Dispersion Equations”,Applied Mathematical Sciences , Vol. 6, No. 98,(2012), pp. 4869 - 4878. [11] Zafer Cakir, ”Stability of Difference Schemes for Fractional Parabolic PDE with the Dirichlet-Neumann Conditions”,Abstract and Applied Analysis, Volume 2012, Article ID 463746, 17 pages

International Journal of Applied Mathematical Research

57

˘ [12] Ibrahim KARATAY, Nurdane KALE, Serife Rabia BAYRAMOGLU,” A new difference scheme for time fractional heat equations based on the Crank-Nicholson method”, Volume 16, Issue 4,(2013), pp 892-910. [13] A. Ashyralyev and Z. Cakir,” On the numerical solution of fractional parabolic partial differential equations with the Dirichlet condition”,In Proceedings of the 2nd International Symposium on Computing in Science and Engineering (ISCSE ’11) , (2011), pp. 529-530.