FRACTIONAL REACTION-DIFFUSION EQUATIONS 7 0 R.K. SAXENA 0 Department of Mathematics and Statistics, Jai Narain Vyas University 2 Jodhpur, 342005, India n a J A.M. MATHAI 3 Department of Mathematics and Statistics, McGill University ] Montreal, Canada H3A 2K6 A C H.J. HAUBOLD . h Office for Outer Space Affairs, United Nations t a P.O.Box 500, A1400 Vienna, Austria m [ 2 Abstract. In a series of papers, Saxena, Mathai, and Haubold (2002, 2004a, v 2004b) derived solutions of a number of fractional kinetic equations in terms 3 7 of generalized Mittag-Leffler functions which provide the extension of the 4 work of Haubold and Mathai (1995, 2000). The subject of the present paper 4 is to investigate the solution of a fractional reaction-diffusion equation. The 0 6 results derived are of general nature and include the results reported earlier 0 by many authors, notably by Jespersen, Metzler, and Fogedby (1999) for / h anomalousdiffusion anddel-Castillo-Negrete, Carreras, andLynch (2003)for t a reaction-diffusionsystems withLévyflights. Thesolutionhasbeendeveloped m in terms of the H-function in a compact form with the help of Laplace and : v Fourier transforms. Most of the results obtained are in a form suitable for i X numerical computation. r a 1 Introduction Reaction-diffusion models have found numerous applications in pattern formation in biology, chemistry, and physics, see Murray (2003), Kuramoto (2003),WilhelmssonandLazzaro(2001),andHundsdorferandVerwer(2003). These systems indicate that diffusion can produce spontaneous formation of spatio-temporalpatterns. Fordetails, onecanrefertotheworkofNicolis and Prigogine (1977) and Haken (2004). A general model for reaction-diffusion systems is investigated by Henry and Wearne (2000, 2002) and Henry, Lang- lands, and Wearne (2005). The simplest reaction-diffusion models are of the form ∂N ∂2N = d +F(N),N = N(x,t), (1) ∂t ∂x2 where d is the diffusion coefficient and F(N) is a nonlinear function rep- resenting reaction kinetics. It is interesting to observe that for F(N) = γN(1 − N2), eq. (1) reduces to the Fisher-Kolmogorov equation and if we set F(N) = N(1 − N2), it gives rise to the real Ginsburg-Landau equation. Del-Castillo-Negrete, Carreras, and Lynch (2002) studied the front propagation and segregation in a system of reaction-diffusion equations with cross-diffusion. Recently, del-Castillo-Negrete, Carreras, and Lynch (2003) discussed the dynamics in reaction-diffusion systems with non-Gaussian diffusion caused by asymmetric Lévy flights and solved the following model ∂N = ηDαN +F(N),N = N(x,t), (2) ∂t x with F = 0. In this paper we present a solution of a more general model of reaction- diffusion systems (2) in which ∂N has been replaced by ∂βN,β > 0. This ∂t ∂tβ new model extends the work of Jespersen, Metzler, and Fogedby (1999) and del-Castillo-Negrete, Carreras, and Lynch (2003). Most of the results are obtained in a compact form suitable for numerical computation. AgeneralizationoftheMittag-Lefflerfunction(Mittag-Leffler,1903,1905) ∞ zn E (z) := ,α ∈ C,Re(α) > 0, (3) α Γ(nα+1) n=0 X was introduced by Wiman (1905) in the generalized form ∞ zn E (z) := ,α,β ∈ C,Re(α) > 0. (4) α,β Γ(nα+β) n=0 X The main results of these functions are available in the handbook of Erdélyi, Magnus, Oberhettinger, and Tricomi (1955, Section 18.1) and monographs by Dzherbashyan (1966, 1993). The H-function is defined by means of a Mellin-Barnes type integral in the following manner (Mathai and Saxena, 1978) Hm,n(z) = Hm,n z (ap,Ap) (5) p,q p,q (bq,Bq) h (cid:12) i 1 = Hm,n z(cid:12)(cid:12)(a1,A1),...,(ap,Ap) = Θ(ξ)z−ξdξ, p,q h (cid:12)(b1,B1),...,(Bq,Bq)i 2πi ZΩ (cid:12) where i = (−1)1/2, (cid:12) Πm Γ(b +B ξ) Πn Γ(1−a −A ξ) j=1 j j j=1 j j Θ(ξ) = , (6) Πhq Γ(1−b −Bihξ) Πp Γ(a +Aiξ) j=m+1 j j j=n+1 j j h ih i and an empty product is always interpreted as unity; m,n,p,q ∈ N with 0 0 ≤ n ≤ p,1 ≤ m ≤ q,A ,B ∈ R ,a ,b ∈ R or C(i = 1,...,p;j = 1,...,q) i j + i j such that A (b +k) 6= B (a −l−1),k,l ∈ N ;i = 1,...,n;j = 1,...,m, (7) i j j i 0 where we employ the usual notations:N = (0,1,2...),R = (−∞,∞),R = 0 + (0,∞) and C being the complex number field. The contour Ω is either L ,L , or L . These contours are defined in the monographs by Prud- −∞ +∞ iγ∞ nikov, Brychkov, and Marichev (1989), Mathai (1993), and Kilbas and Saigo (2004). However, the explicit definitions of these contours are given below to complete the definition of the H-function. (i) Ω = L is a left loop situated in a horizontal strip starting at the −∞ point −∞+iϕ−1 and terminating at the point −∞+iϕ with −∞ < ϕ < 2 1 ϕ < +∞; 2 (ii) Ω = L is a right loop situated in a horizontal strip starting at the +∞ point +∞+iϕ and terminating at the point +∞+iϕ with −∞ < ϕ < 1 2 1 ϕ < +∞. 2 (iii) Ω = L is a contour starting at the point γ −i∞ and terminating iγ∞ at the point γ +i∞, where γ ∈ R = (−∞,+∞). A detailed and comprehensive account of the H-function is available from the monograph by Mathai and Saxena (1978), Prudnikov, Brychkov, and Marichev (1989), and Kilbas and Saigo (2004). The relation connecting Ψ (z) and the H-function is given for the first time in the monograph by p q Mathai and Saxena (1978, p.11, Eq.1.7.8) as Ψ (a1,A1),...,(ap,Ap)|z = H1,p −z (1−a1,A1),...,(1−ap,Ap) , (8) p q (b1,B1),...,(bq,Bq) p,q+1 (0,1),(1−b1,B1),...,(1−bq,Bq) h i h (cid:12) i (cid:12) where pΨq(z) is Wright’s generalized hyperg(cid:12)eometric function (Wright, 1940); also see (Erdélyi, Magnus, Oberhettinger, and Tricomi, 1953, Section 4.1), defined by means of the series representation in the form ∞ [Πp Γ(a +A r)]zr Ψ (z) = Ψ (ap,Ap)|z = j=1 j j , (9) p q p q (bq,Bq) [Πq Γ(b +B r)(r)!] h i rX=0 j=1 j j where z ∈ C,a ,b ∈ C,A ,B ∈ R = (−∞,∞),A ,B 6= 0 (i = 1,...,p;j = i j i j i j 1,...,q), q B − p A > −1;C being the set of complex numbers and j=1 j j=1 j Γ(z) is Euler’s gamma function. This function includes many special func- P P tions. It is interesting to observe that for A = B = 1,∀i and j, eq. (9) i j reduces to a generalized hypergeometric function F (z) as p q Πp Γ(a ) Ψ (ap,1) z = j=1 j F (a ,...,a ;b ,...,b ;z), (10) p q (bq,1) Πq Γ(b ) p q 1 p 1 q h (cid:12) i j=1 j (cid:12) (cid:12) where a 6= −ν(j = 1,...,p and ν = 0,1,2,....);p < q or p = q,|z| < 1. j Prior to (9), Wright (1933) introduced a special case of (9) in the form ∞ 1 zr Φ(a,b;z) = Ψ − z = , 0 1 (b,a) Γ(ar +b)(r)! h (cid:12) i rX=0 (cid:12) (cid:12) which widely occurs in problems of fractional diffusion. It has been shown by Saxena, Mathai, and Haubold (2004b) that (1,1) E (z) = Ψ z (11) α,β 1 1 (β,α) = H1,1h−z (cid:12)(cid:12)(0,i1) . (12) 1,2 (cid:12)(0,1),(1−β,α) h (cid:12) i (cid:12) If we further take β = 1 in (11) and (12), w(cid:12)e find that (1,1) E (z) = E (z) = Ψ z (13) α,1 α 1 1 (1,α) = H1,1 −z (0,1) h (cid:12)(cid:12), i (14) 1,2 (0,1),(0,α) (cid:12) h (cid:12) i (cid:12) (cid:12) where Re(α) > 0,α ∈ C. FromMathai andSaxena (1978)and Prudnikov, Brychkov, and Marichev (1989, p.355, Eq.2.25.3), it follows that the Laplace transform of the H- function is given by L tρ−1Hm,n ztσ (ap,Ap) = s−ρHm,n+1 zs−σ (1−ρ,σ),(ap,Ap) , (15) p,q (bq,Bq) p+1,q (bq,Bq) n h (cid:12) io h (cid:12) i (cid:12) (cid:12) where σ > 0,Re(s) >(cid:12) 0,Re[ρ + σmin (bj )] > 0,(cid:12)|argz| < [π/2]θ,θ > 0; 1≤j≤m Bj θ = n A − p A + m B − q B . j=1 j j=n+1 j j=1 j j=m+1 j By virtue of the cancellation law for the H-function (Mathai and Saxena, P P P P 1978), it can be readily seen that L−1 s−ρHm,n zsσ (ap,Ap) = tρ−1Hm,n zt−σ (ap,Ap),(ρ,σ) , (16) p,q (bq,Bq) p+1,q (bq,Bq) n h (cid:12) io h (cid:12) i (cid:12) (cid:12) where σ > 0,Re(s) > 0,Re[(cid:12)ρ+σmax (1−aj)] > 0,|argz(cid:12)| < 1πθ ,θ > 0; 1≤j≤n Aj 2 1 1 θ = θ−a. Two interesting special cases of (16) are worth mentioning. If we employ the identity (Mathai and Saxena, 1978) H1,0 x = xαexp(−x), (17) 0,1 (α,1) h (cid:12) i (cid:12) we obtain (cid:12) L−1[s−ρexp(−zsσ)] = tρ−1H1,0 zt−σ (ρ,σ) , (18) 1,1 (0,1) h (cid:12) i where Re(s) > 0,σ > 0. (cid:12) (cid:12) Further if we use the identity (Mathai and Saxena, 1978) H2,0 x = 2K (2x1/2), (19) 0,2 (ν,1)(−ν,1) ν 2 2 h (cid:12) i (cid:12) eq.(16) yields (cid:12) z2t−2σ 2L−1[s−ρK (zsσ)] = tρ−1H2,0 (ρ,2σ) , (20) ν 1,2 " 4 (ν2,1),(−ν2,1)# (cid:12) (cid:12) (cid:12) where Re(ρ) > 0,Re(z2) > 0,Re(s) > 0, and K (.) is the Bessel function of ν the third kind . In view of the result of Saxena, Mathai, and Haubold (2004a, p.49), also see Prudnikov, Brychkov, and Marichev (1989, p.355, eq.(2.25.3.2)), the cosine transform of the H-function is given by ∞ tρ−1cos(kt)Hm,n atµ (ap,Ap) dt (21) Z0 p,q h (cid:12)(bq,Bq)i = π Hn+1,m kµ (1−bq,Bq),((cid:12)(cid:12)1+2ρ,µ2) , kρ q+1,p+2" a (cid:12)(ρ,µ),(1−ap,Ap),(1+2ρ,µ2)# (cid:12) (cid:12) where Re[ρ+µmin ] > 1,|arg a| < 1(cid:12)πθ;θ > 0,θ is defined with the result 1≤j≤m 2 eq. (15). The Riemann-Liouville fractional integral of order ν is defined by (Miller and Ross, 1993, p.45) 1 t D−νf(t) = (t−u)ν−1f(u)du, (22) 0 t Γ(ν) Z0 where Re(ν) > 0. Following Samko, Kilbas, and Marichev (1990, p.37), we define the fractional derivative for α > 0 in the form 1 dn t f(u)du Dαf(t) = , n = [α]+1, (23) 0 t Γ(n−α)dtn (t−u)α−n+1 Z0 where [α] means the integral part of the number α. In particular, if 0 < α < 1, d 1 t f(u)du Dαf(t) = , (24) 0 t dtΓ(1−α) (t−u)α Z0 and if α = n ∈ N = {1,2,...,}, then Dnf(t) = Dnf(t)(D = d/dt), (25) 0 t is the usual derivative of order n. From Erdélyi, Magnus, Oberhettinger, and Tricomi (1954b, p. 182), we have L D−νf(t) = s−νF(s), (26) 0 t n o where ∞ F(s) = L{f(t);s} = f∗(s) = exp(−st)f(t)dt,Re(s) > 0. (27) Z0 The Laplace transform of the fractional derivative is given by Oldham and Spanier (1974, p.134, eq.(8.1.3)) n L{ Dαf(t)} = sαF(s)− sr−1 Dα−rf(t)| . (28) 0 t 0 t t=0 r=1 X In certain boundary-value problems, the following fractional derivative of order α > 0 is introduced by Caputo (1969) in the form 1 t f(m)(τ)dτ Dαf(t) = , (29) t Γ(m−α) (t−τ)α+1−m Z0 m−1 < α ≤ m,Re(α) > 0,m ∈ N. dmf = ,if α = m. (30) dtm Caputo (1969) has given the Laplace transform of this derivative as m−1 L{Dαf(t);s} = sαF(s)− sα−r−1fr(0+),m−1 < a ≤ m. (31) t r=0 X The above formula is very useful in deriving the solution of differ-integral equations of fractional order governing certain physical problems of reaction and diffusion. We also need the Weyl fractional operator defined by 1 dn t f(u)du Dµf(t) = , (32) −∞ x Γ(n−µ)dtn (t−u)µ−n+1 Z−∞ where n = [µ] is an integral part of µ > 0. Its Fourier transform is (Metzler and Klafter, 2000, p.59, A.11, 2004) F { Dµf(x)} = (ik)µf˜(k), (33) −∞ x where we define the Fourier transform as ∞ ˜ h(q) = h(x)exp(iqx)dx. (34) Z−∞ Following the convention initiated by Compte (1996), we suppress the imag- inary unit in Fourier space by adopting the slightly modified form of above result in our investigations (Metzler and Klafter, 2000, p.59, A.12, 2004) F { Dµf(x)} = −|k|µf˜(k) (35) −∞ x instead of (33). Finally we also need the following property of the H-function (Mathai and Saxena, 1978) 1 Hm,n xδ (ap,Ap) = Hm,n x (ap,Ap/δ) , (36) p.q (bq,Bq) δ p,q (bq,Bq/δ) h (cid:12) i h (cid:12) i (cid:12) (cid:12) where δ > 0. (cid:12) (cid:12) 2 The Fractional Reaction-Diffusion Equation In this section we will investigate the solution of the reaction-diffusion equation (37). The result is given in the form of the following theorem. Theorem. Consider the following fractional reaction-diffusion model ∂βN(x,t) = η DαN(x,t)+ϕ(x,t);η,t > 0,x ∈ R,0 < β ≤ 2, (37) ∂tβ −∞ x with the initial condition N(x,0) = f(x),N (x,0) = g(x) for x ∈ R, (38) t where N (x,0) means the first partial derivative of N(x,t) with respect to t ϕ evaluated at t = 0,η is a diffusion constant and ϕ(x,t) is a nonlinear function belonging to the area of reaction-diffusion. Then for the solution of (37), subject to the initial conditions (38), there holds the formula 1 ∞ N(x,t) = f˜(k)E (−η|k|αtβ)exp(−ikx)dk (39) β,1 2π Z−∞ 1 ∞ + tg˜(k)E (η|k|αtβ)exp(−ikx)dk β,2 2π Z−∞ 1 t ∞ + ξβ−1 ϕ˜(k,t−ξ)E (−η|k|αξβ)exp(−ikx)dkdξ. β,β 2π Z0 Z−∞ Proof. If we apply the Laplace transform with respect to the time variable t and Fourier transform with respect to space variable x and use the initial conditions and (32), then the given equation transforms into the form f˜(k)sβ−1 g˜(k)sβ−2 ϕ˜∗(k) N˜∗(k,s) = + + . (40) sβ +η|k|α sβ +η|k|α sβ +η|k|α On taking the inverse Laplace transform of (40) and applying the result sa−1 L−1 = tβ−αE (−atβ), (41) (a+sβ) β,β−α+1 where Re(s) > 0,Re(β −α+1) > 0, it is seen that N˜(k,t) = f˜(k)E (−η|k|αtβ)+g˜(k)tE (−η|k|αtβ) (42) β,1 β,2 t + ϕ˜(k,t−ξ)ξβ−1E (−η|k|αξβ)dξ. β,β Z0 The required solution (39) is now obtained by taking the inverse Fourier transform of (38). Thus, we have 1 ∞ N(x,t) = f˜(k)E (−η|k|αtβ)exp(−ikx)dk β,1 2π Z−∞ 1 ∞ + tg˜(k)E (−η|k|αtβ)exp(−ikx)dk β,2 2π Z∞ 1 t ∞ + ξβ−1 ϕ∗(k,t−ξ)E (−η|k|αξβ)exp(−ikx)dkdξ. β,β 2π Z0 Z−∞ This completes the proof of the theorem. If we set α = 2, we obtain the result given by Debnath (2003). Note. By virtue of the identity (12), the solution (39) can be expressed in terms of the H-function as can be seen from the solutions given in the special cases of the theorem in the next section. 3 Special Cases When g(x) = 0, then applying the convolution theorem of the Fourier transform to the solution (35), the theorem yields Corollary 1.1. The solution of fractional reaction-diffusion equation ∂β N(x,t)−η DαN(x,t) = ϕ(x,t), x ∈ R,t > 0,η > 0, (43) ∂tβ −∞ x subject to the initial conditions N(x,0) = f(x),N (x,0) = 0 for x ∈ R,1 < β ≤ 2, (44) t where η is a diffusion constant and ϕ(x,t) is a nonlinear function belonging to the area of reaction-diffusion, is given by ∞ N(x,t) = G (x−τ,t)f(τ)dτ (45) 1 Z−∞ t x + (t−ξ)β−1 G (x−τ,t−ξ)ϕ(τ,ξ)dτdξ, 2 Z0 Z0 where 1 ∞ G (x,t) = exp(−ikx)E (−η|k|αtβ)dk (46) 1 β,1 2π Z−∞ 1 ∞ = cos(kx)H1,1 kη1/αtβ/α (0,1/α) dk πα 1,2 (0,1/α),(0,β/α) Z0 h (cid:12) i 1 |x| (cid:12) = H2,1 (1,1/α),(1,β/α)(cid:12),(1,1/2) ,α > 0, α|x| 3,3 "η1/αtβ/α (1,1),(1,1/α),(1,1/2) # (cid:12) (cid:12) (cid:12) 1 ∞ G (x,t) = exp(−ikx)E (−η|k|αtβ)dk (47) 2 β,β 2π Z−∞ 1 ∞ = cos(kx)H1,1 kη1/αtβ/α (0,1/α) dk πα 1,2 (0,1/α),(1−β,β/α) Z0 h (cid:12) i 1 |x| (cid:12) = H2,1 (1,1/α),(β,β/α(cid:12)),(1,1/2) ,α > 0. α|x| 3,3 "η1/αtβ/α (1,1),(1,1/α),(1,1/2) # (cid:12) (cid:12) (cid:12) If we set f(x) = δ(x),ϕ ≡ 0,g(x) = 0, where δ(x) is the Dirac-delta function, then we arrive at the following Corollary 1.2. Consider the following reaction-diffusion model ∂βN(x,t) = η DαN(x,t),η > 0,x ∈ R,0 < β ≤ 1, (48) ∂tβ −∞ x with the initial condition N(x,t = 0) = δ(x), where η is a diffusion constant and δ(x) is the Dirac-delta function. Then the solution of (44) is given by 1 |x| N(x,t) = H2,1 (1,1/α),(1,β/α),(1,1/2) . (49) α|x| 3,3 "(ηtβ)1/α (1,1),(1,1/α),(1,1/2) # (cid:12) (cid:12) (cid:12) In the case β = 1, then in view of the cancellation law for the H-function (Mathai and Saxena, 1978), (49) gives rise to the following result given by

Fractional reaction-diffusion equations

Similar Fractional reaction-diffusion equations