Stochastic Solution of Fractional Fokker-Planck Equations with Space-Time-Dependent Coefficients
aa r X i v : . [ m a t h . P R ] O c t STOCHASTIC SOLUTION OF FRACTIONAL FOKKER-PLANCKEQUATIONS WITH SPACE-TIME-DEPENDENT COEFFICIENTS
ERKAN NANE AND YINAN NI
Abstract.
This paper develops solutions of fractional Fokker-Planck equations describing sub-diffusion of probability densities of stochastic dynamical systems driven by non-Gaussian L´evyprocesses, with space-time-dependent drift, diffusion and jump coefficients, thus significantlyextends Magdziarz and Zorawik’s result in [14]. Fractional Fokker-Planck equation describingsubdiffusion is solved by our result in full generality from perspective of stochastic representa-tion. Introduction
Fractional Fokker-Planck equation has shown its application in diverse scientific areas, in-cluding biology [9], physics [6], [28], finance [7]. For example, in physics, it has been broadlyused to describe phenomena related to anomalous diffusion [27], [29]. Many different typesof fractional Fokker-Planck equation have been solved in terms of probability density function(PDF) of corresponding process [20], [21], related researches are also growing rapidly in differentbranches, e.g., [15], [16], [17], [18], [22]. Recently, Magdziarz and Zorawik [14] provide solutionto an extended type of Factional Fokker-Planck equation.To state the result, let’s introduce subordinator T Ψ ( t ) with Laplace transform(1.1) E ( e − uT Ψ ( t ) ) = e − t Ψ( u ) , where Laplace exponent(1.2) Ψ( u ) = Z ∞ (1 − e − ux ) ν ( dx ) ,ν is L´evy measure satisfying(1.3) Z R −{ } min {| x | , } ν ( dx ) < ∞ and ν ((0 , ∞ )) = ∞ . Then, the first passenge time process defined as(1.4) S Ψ ( t ) = inf { γ > T Ψ ( γ ) > t } , t ≥ , is called the inverse subordinator.Define the integro-differential operator Φ t as(1.5) Φ t f ( t ) = ∂∂t Z t M ( t − y ) f ( y ) dy, where the function f is smooth enough and kernel M ( t ) is defined by its Laplace transform as(1.6) ˜ M ( u ) = Z ∞ e − ut M ( t ) dt = 1Ψ( u ) . A stochastic process X = { X ( t ) , t ≥ } is a L´evy process if (a) X (0) = 0 , a.s. , (b) X hasindependent and stationary increments, (c) X is stochastically continuous in time. X − ( t ) isused to denote left limit, X − ( t ) = lim s → t − X ( s ). Key words and phrases.
Fokker-Plank equation, space-time dependent coefficients, L´evy process, subdiffusion. y Theorem 1.2.14 and Proposition 1.3.1 of [1], a L´evy process X has characteristics ( b, A, ν ),that’s,(1.7) E ( e iuX ( t ) ) = exp " t ibu − Au + Z R −{ } ( e iuy − − iuy B ( y )) ν ( dy ) ! , where ν is L´evy measure. In the remaining part of this paper, for convenience, we will decomposea L´evy process X ( t ) into 3 parts: drift, Brownian motion B ( t ) and pure jump L´evy process L ( t ).Improving on the methods in the papers [5] [10] [12], Magdziarz and Zorawik [14] proved thatthe PDF of process X ( t ) = Y − ( S Ψ ( t )), where(1.8) dY ( t ) = F ( Y − ( t ) , T − Ψ ( t )) dt + σ ( Y − ( t ) , T − Ψ ( t )) dB ( t ) + E ( T − Ψ ( t )) dL ( t ) , t ≥ ,Y (0) = 0 , T Ψ ( t ) = 0 , with F ( x, t ) , σ ( x, t ) , E ∈ C ( R ) satisfying Lipschitz condition, solves fractional Fokker-Planckequation(1.9) ∂q ( x, t ) ∂t = (cid:20) − ∂∂x F ( x, t ) + 12 ∂ ∂x σ ( x, t ) (cid:21) Φ t q ( x, t )+ Z R −{ } (cid:20) Φ t q ( r, t ) | r = x − E ( t ) y − Φ t q ( x, t ) + E ( t ) y ∂∂x Φ t q ( x, t ) B ( y ) (cid:21) ν ( dy ) , with q ( x,
0) = δ ( x ), where is an indicator function, B = { y, | y | < } .This result extends the following celebrated fractional Fokker-Planck equation introduced byMetzler and Klafter [21] in 2000,(1.10) ∂q ( x, t ) ∂t = D − αt (cid:20) − ∂∂x F ( x ) + σ ∂ ∂x (cid:21) q ( x, t ) , with σ > q ( x,
0) = δ ( x ), which describes anomalous diffusion in the presence of an space-dependent force F ( x ). Note that the operator D − αt , α ∈ (0 ,
1) is fractional derivative ofRiemann-Liouville type [25],(1.11) D − αt f ( t ) = 1Γ( α ) ddt Z t ( t − s ) α − f ( s ) ds, for f ∈ C ([0 , ∞ )). Magdziarz et al. [13] showed that the PDF of the subordinated process X ( t ) = Y ( S α ( t )) solves equation (1.10), where(1.12) dY ( t ) = F ( Y ( t )) dt + σdB ( t ) , Y (0) = 0 . Here, let D α be an α -stable subordinator with Laplace transform E [ e − uD α ( γ ) ] = e − tu α , itsinverse S α ( t ) is defined as(1.13) S α ( t ) = inf { γ > D α ( γ ) > t } . Magdziarz and Zorawik’s result [14] also extends the following fractional Fokker-Planck equa-tion introduced by Sokolov and Klafter [28] in 2009,(1.14) ∂q ( x, t ) ∂t = (cid:20) − F ( t ) ∂∂x + σ ∂ ∂x (cid:21) D − αt q ( x, t ) , with σ > q ( x,
0) = δ ( x ), where external force F ( t ) is time-dependent. Its solution isobtained by PDF of subordinated process X ( t ) = Y ( S α ( t )), where(1.15) dY ( t ) = F ( T α ( t )) dt + σdB ( t ) , Y (0) = 0 . ne of the main results of this paper is that the PDF of X ( t ) = Y − ( S Ψ ( t )), where(1.16) dY ( t ) = F ( Y − ( t ) , T − Ψ ( t )) dt + σ ( Y − ( t ) , T − Ψ ( t )) dB ( t ) + h ( Y − ( t ) , T − Ψ ( t )) dL ( t ) , t ≥ ,Y (0) = 0 , T Ψ ( t ) = 0 , solves the following fractional Fokker-Planck equation, which involves fractional Laplacian op-erator,(1.17) ∂q ( x, t ) ∂t = (cid:20) − ∂∂x F ( x, t ) + 12 ∂ ∂x σ ( x, t ) − ( − ∆) α (sgn( h ( x, t )) | h ( x, t ) | α ) (cid:21) Φ t q ( x, t ) , with q ( x,
0) = δ ( x ), note that α ∈ (0 ,
2) and L ( t ) is the α -stable L´evy process.Furthermore, in section 3, we extend Magdziarz and Zorawik’s result [14] by solving(1.18) ∂∂t w ( x, t ) = (cid:20) − ∂∂x ( F ( x, t ) + ∂ ∂x σ ( x, t ) + T l + (cid:21) Φ t w ( x, t ) , with q ( x,
0) = δ ( x ), where(1.19) T l + f ( x, t ) = Z R −{ } " ∞ X k =1 ( − r ) k k ! ∂ k ∂y k ( h ( x, t ) k f ( x, t )) + B ( r, t ) r ∂∂y ( h ( x, t ) f ( x, t )) ν ( dr ) , for any f ( x, t ) ∈ C ∞ ( R ).Note that coefficient E ( t ) of pure jump L´evy process in (1.9) is time-dependent, while co-efficient h ( x, t ) of pure jump L´evy process in (1.18) is space-time-dependent. Thus, (1.9) is aspecial case of (1.18) when h ( x, t ) only depends on time t . For more details on this, see remark3.9.In the remaining of this paper, necessary concepts will be given in the Preliminaries section; inthe Main Results section, we will solve 3 different fractional Fokker-Planck equations involvingoperators of α -stable, symmetric and general L´evy processes, respectively, one by one.2. Preliminaries
Let X = { X ( t ) , t ≥ } be a L´evy process with characteristics ( b, A, ν ), by Theorem 6.7.4 of[1], it has infinitesimal generator(2.1) Af ( x ) = b ∂∂x f ( x ) + 12 ∂ ∂x f ( x ) + Z R −{ } (cid:20) f ( x + y ) − f ( x ) − y ∂∂x f ( x ) B ( y ) (cid:21) ν ( dy ) , for each f ∈ C ( R ) , x ∈ R .The following L´evy processes and their generators will be used in this paper later.An α -stable L´evy process X ( t ) has characteristics (0 , , ν ) and L´evy symbol η ( u ) = −| u | α , α ∈ (0 , Af ( x ) = Z R −{ } [ f ( x + y ) − f ( x )] ν ( dy ) = − ( − ∆) α/ f ( x ) , where ν ( dy ) = C α dy | y | α , C α = α Γ(1 − α ) .A symmetric L´evy process X ( t ) has characteristics (0 , , ν ) where ν is symmetric L´evy mea-sure, that’s, ν ( B ) = ν ( − B ), for each B ⊂ R , and infinitesimal generator(2.3) Af ( x ) = Z R −{ } [ f ( x + y ) − f ( x )] ν ( dy ) . A general L´evy process X ( t ) with characteristics (0 , , ν ) has infinitesimal generator(2.4) Af ( x ) = Z R −{ } (cid:20) f ( x + y ) − f ( x ) − y ∂∂x f ( x ) B ( y ) (cid:21) ν ( dy ) , n fact, such X ( t ) is a pure jump L´evy process since drift and Brownian motion parts are goneas b = A = 0.Let G ( w ) = ν (( w, ∞ )), define the following operator [14](2.5) Θ w g ( w ) = Z w G ( w − z ) g ( z ) dz, with Laplace transform of its kernel(2.6) ˜ G ( u ) = Ψ( u ) u . Main Results
In this section, we will solve fractional Fokker-Planck equations with infinitesimal generatorof L´evy processes with space-time-dependent coefficients for drift, diffusion and jump parts. Weanalyze the cases when the jump process is α -stable and symmetric L´evy process, respectively,before analyzing general L´evy process, since the results for the cases of α -stable and symmetricL´evy process are more explicit.3.1. Fractional Fokker-Planck equation with α -stable L´evy generator.Definition 3.1. Let L ( t ) be α -stable L´evy process satisfying the following SDE (3.1) dL ( t ) = Z B x ˜ N ( dt, dx ) + Z B c xN ( dt, dx ) , where L´evy symbol η ( u ) = R R −{ } [ e iux − − iux B ( x )] ν ( dx ) , with ν ( dx ) = C α dx | x | α , B = { x, | x | < } . Theorem 3.2.
Suppose that the standard Brownian motion B ( t ) , the α -stable L´evy process L ( t ) as defined in definition 3.1 and subordinator T Ψ ( t ) are independent, where T Ψ ( t ) has Laplaceexponent Ψ( u ) and its inverse S Ψ ( t ) . Let Y ( t ) be the solution of the stochastic equation (3.2) dY ( t ) = F ( Y − ( t ) , T − Ψ ( t )) dt + σ ( Y − ( t ) , T − Ψ ( t )) dB ( t ) + h ( Y − ( t ) , T − Ψ ( t )) dL ( t ) , t ≥ ,Y (0) = 0 , T Ψ (0) = 0 where the function F ( x, t ) , σ ( x, t ) , h ( x, t ) ∈ C ( R ) satisfy the Lipschitz condition. As-sume that the PDF of the process ( Y ( t ) , T Ψ ( t )) , p t ( y, z ) exists. Furthermore, we assume that ∂∂t p t ( y, z ) , ∂∂y p t ( y, z ) , ∂ ∂y p t ( y, z ) exist, (3.3) Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂t p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdt < ∞ for each t > , (3.4) Z x x Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂x p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdx < ∞ , Z x x Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ ∂x p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdx < ∞ , for each x , x ∈ R and (3.5) Z ∞ Z R −{ } (cid:12)(cid:12)(cid:12)(cid:12) p s ( x + y, t ) sgn ( h ( x + y, t )) | h ( x + y, t ) | α | y | α − p s ( x, t ) sgn ( h ( x, t )) | h ( x, t ) | α | y | α (cid:12)(cid:12)(cid:12)(cid:12) C α dyds < ∞ . for each x ∈ R .Then, the PDF of the process X ( t ) = Y − ( S Ψ ( t )) is the weak solution of fractionalFokker-Planck equation (1.17) . roof. This proof uses methods in [5] and [14] with crucial changes.Equation (3.2) can be represented as the following stochastic equations(3.6) dY ( t ) = F ( Y − ( t ) , Z − ( t )) dt + σ ( Y − ( t ) , Z − ( t )) dB ( t )+ h ( Y − ( t ) , Z − ( t )) Z B x ˜ N ( dt, dx ) + h ( Y − ( t ) , Z − ( t )) Z B c xN ( dt, dx ) dZ ( t ) = dT Ψ ( t ) . By Theorem 6.7.4 of [1], the infinitesimal generator of the process ( Y ( t ) , Z ( t )) that operateson functions f ∈ C ( R ) is(3.7) Γ f ( y, z ) = F ( y, z ) ∂∂y f ( y, z ) + 12 σ ( y, z ) ∂ ∂y f ( y, z )+ Z R −{ } (cid:20) f ( y + xh ( y, z ) , z ) − f ( y, z ) − B ( x ) xh ( y, z ) ∂∂y f ( y, z ) (cid:21) ν ( dx )+ Z ∞ [ f ( y, z + u ) − f ( y, z )] µ ( du ) , where µ is L´evy measure of T Ψ ( t ).Decompose Γ = A + T α , where(3.8) Af ( y, z ) = F ( y, z ) ∂∂y f ( y, z ) + 12 σ ( y, z ) ∂ ∂y f ( y, z ) , + Z ∞ [ f ( y, z + u ) − f ( y, z )] µ ( du ) T α f ( y, z ) = Z R −{ } (cid:20) f ( y + xh ( y, z ) , z ) − f ( y, z ) − B ( x ) xh ( y, z ) ∂∂y f ( y, z ) (cid:21) ν ( dx ) . By setting E ( t ) = 0 in (2.17) of [14],(3.9) A + p t ( y, z ) = − ∂∂y ( F ( y, z ) p t ( y, z )) + ∂ ∂y ( 12 σ ( y, z ) p t ( y, z )) − ∂∂z Θ z p t ( y, z ) , where A + is the Hermitian adjoint of A .Next, T + α is derived below(3.10) T α f ( y, z ) = Z R −{ } [ f ( y + xh ( y, z ) , z ) − f ( y, z )] C α dx | x | α =sgn( h ( y, z )) | h ( y, z ) | α Z R −{ } [ f ( y + x, z ) − f ( y, z )] C α dx | x | α =sgn( h ( y, z )) | h ( y, z ) | α [ − ( − ∆) α/ f ( y, z )] . Note that the second equation is the result of change of variable.Since the infinitesimal generator of α -stable L´evy process is self-adjoint,(3.11) Z R f ( y, z ) T + α p t ( y, z ) dydz = Z R p t ( y, z ) T α f ( y, z ) dydz = Z R p t ( y, z )sgn( h ( y, z )) | h ( y, z ) | α [ − ( − ∆) α/ f ( y, z )] dydz = Z R f ( y, z )[ − ( − ∆) α/ ( p t ( y, z )sgn( h ( y, z )) | h ( y, z ) | α )] dydz. hat’s,(3.12) T + α p t ( y, z ) = − ( − ∆) α/ ( p t ( y, z )sgn( h ( y, z )) | h ( y, z ) | α )Since Γ + = A + + T + α and ∂∂t p t ( y, z ) = Γ + p t ( y, z ), we have(3.13) ∂∂t p t ( y, z ) =Γ + p t ( y, z ) = A + p t ( y, z ) + T + α p t ( y, z )= − ∂∂y ( F ( y, z ) p t ( y, z )) + ∂ ∂y ( 12 σ ( y, z ) p t ( y, z )) − ∂∂z Θ z p t ( y, z ) − ( − ∆) α/ ( p t ( y, z )sgn( h ( y, z )) | h ( y, z ) | α ) . Next, we establish the relationship between q ( x, t ) and p t ( y, z ), the probability density func-tions of X ( t ) and ( Y ( t ) , Z ( t )), respectively. Let w denote a random path of stochastic process,for each fixed interval I , define indicator function(3.14) I ( x ) = ( x ∈ I ,0 otherwise.and the auxiliary function [5](3.15) H t ( s, w, u ) = ( I ( Y − ( s, w )) if Z − ( s, w ) ≤ t ≤ Z − ( s, w ) + u ,0 otherwise.Then, we get(3.16) Z I q ( x, t ) dx = E [ I ( X ( t, w ))] . Let ∆ Z ( t, w ) = Z ( t, w ) − Z − ( t, w ) and s = S Ψ ( t, w ) = inf { τ > T Ψ ( τ ) > t } , then H t ( S Ψ ( t ) , w, ∆ Z ( t, w )) = I ( X ( t, w )). We can prove this as follows, by definition of H t ( s, w, u ),(3.17) H t ( S Ψ ( t ) , w, ∆ Z ( t, w )) = ( I ( Y − ( S Ψ ( t, w ) , w )) if Z − ( S Ψ ( t, w ) , w ) ≤ t ≤ Z ( S Ψ ( t, w ) , w ),0 otherwise.Since T Ψ ( t, w ) is a subordinator, T − Ψ ( S Ψ ( t, w ) , w ) ≤ t ≤ T Ψ ( S Ψ ( t, w ) , w ) is always true, thus(3.18) H t ( S Ψ ( t, w ) , w, ∆ Z ( t, w )) = I ( Y − ( S Ψ ( t, w ) , w )) . To avoid Z ( t, w ) being a compound Poisson process, we set ν ([0 , ∞ ) = ∞ , see Remark 27.3and 27.4 of [26], thus, jumping times of Z ( t, w ) are dense in [0 , ∞ ] almost surely, see Theorem21.3 of [26]. Then, we can derive that H t ( s, w, ∆ Z ( t, w )) = 0 , if s = S Ψ ( t, w )Since if s < S Ψ ( t, w ), as Z ( t ) = T Ψ ( t ) is a subordinator, Z ( s ) < Z − ( S Ψ ( t, w )) ≤ t. Similarly, if s > S Ψ ( t, w ), then Z ( s ) > Z ( S Ψ ( t, w )) ≥ t . In both cases, H t ( s, w, ∆ Z ( t, w )) = 0.Hence,(3.19) I ( X ( t, w )) = X s> H t ( s, w, ∆ Z ( t, w )) . y Compensation Formula in Ch. XII, Proposition (1.10) of [24],(3.20) E "X s> H t ( s, w, ∆ Z ( s, w )) = E (cid:20)Z ∞ Z ∞ H t ( s, w, u ) ν ( du ) ds (cid:21) = E (cid:20)Z ∞ Z ∞ I ( Y ( s, w )) [ Z ( s − ,w ) ,Z ( s − w )+ u ] ( t ) ν ( du ) ds (cid:21) = E (cid:20)Z ∞ I ( Y ( s, w )) Z ∞ [ t − z, ∞ ] ( u ) [0 ,t ] ( z ) ν ( du ) ds (cid:21) = E (cid:20)Z ∞ I ( Y ( s, w )) [0 ,t ] ( Z ( s, w )) ν ( t − Z ( s, w ) , ∞ ) ds (cid:21) = E (cid:20)Z ∞ I ( Y ( s, w )) [0 ,t ] ( Z ( s, w )) G ( t − Z ( s, w ) ds (cid:21) = Z I Z ∞ Z t G ( t − z ) p s ( y, z ) dzdsdy = Z I Z ∞ Θ t p s ( y, z ) dsdy. By (3.16), (3.19) and (3.20), we have(3.21) Z I q ( x, t ) dx = Z I Z ∞ Θ t p s ( y, t ) dsdy. By the arbitrariness of interval I ,(3.22) q ( x, t ) = Z ∞ Θ t p s ( y, t ) ds. Next we claim that(3.23) ∂∂t Θ t p s ( x, t ) = Θ t ∂∂t p s ( x, t ) . To see this, let P s ( x, u ) and G ( u ) be the Laplace transform ( t → u ) of p s ( x, t ) and g ( t ),respectively. Then Laplace transform of ∂∂t Θ t p s ( x, t ) is given as(3.24) L (cid:20) ∂∂t Θ t p s ( x, t ) (cid:21) = u L [Θ t p s ( x, t )] − Θ p s ( x, u L (cid:20)Z t G ( t − z ) p s ( x, z ) dz (cid:21) = uG ( u ) P s ( x, u ) . On the other hand, since T Ψ (0) = 0 a.s. , Laplace transform of Θ t ∂∂t p s ( x, t ) is as below,(3.25) L (cid:20) Θ t ∂∂t p s ( x, t ) (cid:21) = L (cid:20)Z t G ( t − z ) ∂∂z p s ( x, z ) dz (cid:21) = G ( u ) L (cid:20) ∂∂t p s ( x, t ) (cid:21) = G ( u ) { uP s ( x, u ) − p s ( x, } = uG ( u ) P s ( x, u ) . Notice that(3.26) Z t G ( u ) du = Z t Z ( u, ∞ ) ν ( dw ) du = Z (0 , ∞ ) min( w, t ) ν ( dw ) = K < ∞ , y assumption (3.3) and (3.23), we derive that(3.27) Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂t Θ t p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdt = Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) Θ t ∂∂t p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdt = Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12)Z t G ( t − u ) ∂∂z p s ( x, z ) (cid:12)(cid:12)(cid:12)(cid:12) dsdt ≤ Z t Z ∞ Z t G ( t − u ) (cid:12)(cid:12)(cid:12)(cid:12) ∂∂z p s ( x, z ) (cid:12)(cid:12)(cid:12)(cid:12) dudsdt = Z t Z ∞ Z tu G ( t − u ) (cid:12)(cid:12)(cid:12)(cid:12) ∂∂z p s ( x, z ) (cid:12)(cid:12)(cid:12)(cid:12) dtdsdu = Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂z p s ( x, z ) (cid:12)(cid:12)(cid:12)(cid:12) Z tu G ( t − u ) dtdsdu ≤ K Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂z p s ( x, z ) (cid:12)(cid:12)(cid:12)(cid:12) dsdu< ∞ Thus, we can put differentiation on both side of (3.22) and move it inside the integral on therighthand side as below,(3.28) ∂∂t q ( x, t ) = Z ∞ ∂∂t Θ t p s ( x, t ) ds. Next, we claim that(3.29) Z ∞ p s ( x, t ) ds = Φ t q ( x, t ) . By Fubini theorem,(3.30) q ( x, t ) = Z ∞ Θ t p s ( x, t ) ds = Z ∞ Z t G ( t − z ) p s ( x, z ) dzds = Z t G ( t − z ) Z ∞ p s ( x, z ) dsdz = Θ t Z ∞ p s ( x, t ) ds, thus,(3.31) Z ∞ p s ( x, t ) ds = Θ − t q ( x, t ) . To prove Θ − t = Φ t , let ˜ M ( u ), ˜ G ( u ) and ˜ Q ( x, u ) be the Laplace transform of M ( t ), G ( t ) and q ( x, t ), respectively. Since, by (1.2)(3.32) Z ∞ e − ut G ( t ) dt = Z ∞ Z ( t, ∞ ) e − ut ν ( ds ) dt = Z ∞ − e us u ν ( ds ) = Ψ( u ) u , and(3.33) L [ q ( x, z )] = L [ G ( t )] L [Θ − t q ( x, t )] , we have(3.34) L [Θ − t q ( x, t )] = L [ q ( x, z )] L [ G ( t )] = ˜ Q ( x, u )˜ G ( u ) = u ˜ Q ( x, u )Ψ( u ) . lso(3.35) L [Φ t q ( x, t )] = L (cid:20) ddt Z t M ( t − y ) q ( x, y ) dy (cid:21) = u L (cid:20)Z t M ( t − y ) q ( x, y ) dy (cid:21) − u ˜ M ( u ) ˜ Q ( x, u ) = u ˜ Q ( x, u )Ψ( u ) , this shows Φ t q ( x, t ) = Θ − t q ( x, t ), hence R ∞ p s ( x, t ) ds = Φ t q ( x, t ).Since lim s →∞ p s ( x, t ) = 0 and p ( x, t ) = (0 , ( x, t ), by (3.13), (3.28) and (3.29),(3.36) ∂∂t q ( x, t ) = Z ∞ (cid:20) ∂ ∂x ( 12 σ ( x, t ) p s ( x, t )) − ∂∂x ( F ( x, t ) p s ( x, t )) − ( − ∆) α/ ( p s ( x, t )(sgn( h ( x, t )) | h ( x, t ) | α )) − ∂∂s p s ( x, t ) (cid:21) ds = Z ∞ (cid:20) ∂ ∂x ( 12 σ ( x, t ) p s ( x, t )) − ∂∂x ( F ( x, t ) p s ( x, t )) (cid:21) ds + Z ∞ Z R −{ } (cid:18) p s ( x + y, t ) sgn( h ( x + y, t )) | h ( x + y, t ) | α | y | α − p s ( x, t ) sgn( h ( x, t )) | h ( x, t ) | α | y | α (cid:19) C α dyds = Z ∞ (cid:20) ∂ ∂x ( 12 σ ( x, t ) p s ( x, t )) − ∂∂x ( F ( x, t ) p s ( x, t )) (cid:21) ds + Z R −{ } (cid:20)Z ∞ p s ( x + y, t ) ds sgn( h ( x + y, t )) | h ( x + y, t ) | α | y | α − Z ∞ p s ( x, t ) ds sgn( h ( x, t )) | h ( x, t ) | α | y | α (cid:21) C α dy = ∂ ∂x (cid:18) σ ( x, t ) Z ∞ p s ( x, t ) ds (cid:19) − ∂∂x (cid:18) F ( x, t ) Z ∞ p s ( x, t ) ds (cid:19) − ( − ∆) α/ (cid:18) (sgn( h ( x, t )) | h ( x, t ) | α ) Z ∞ p s ( x, t ) ds (cid:19) = ∂ ∂x ( 12 σ ( x, t )Φ t q ( x, t )) − ∂∂x ( F ( x, t )Φ t q ( x, t )) − ( − ∆) α/ ((sgn( h ( x, t )) | h ( x, t ) | α )Φ t q ( x, t )) . In summary,(3.37) ∂∂t q ( x, t ) = (cid:20) ∂∂x σ ( x, t ) − ∂∂x F ( x, t ) − ( − ∆) α/ (sgn( h ( x, t )) | h ( x, t ) | α ) (cid:21) Φ t q ( x, t ) (cid:3) Fractional Fokker-Planck equation with symmetric L´evy generator.
As is known, α -stable L´evy process is a special case of symmetric L´evy process, this subsection solves frac-tional Fokker-Planck equation associated with symmetric L´evy process with space-time-dependentcoefficient, which extends result of previous subsection. Definition 3.3.
Let L ( t ) be a symmetric L ´ evy process satisfying the following SDE (3.38) dL ( t ) = Z B x ˜ N ( dt, dx ) + Z B c xN ( dt, dx ) here L ´ evy symbol η ( u ) = R R [ cos ( u, x ) − ν ( dx ) and ν is a symmetric L ´ evy measure. Theorem 3.4.
Suppose that the standard Brownian motion B ( t ) , the symmetric L´evy process L ( t ) as defined in definition 3.3 and subordinator T Ψ ( t ) are independent, where T Ψ ( t ) has Laplaceexponent Ψ( u ) and its inverse S Ψ ( t ) . Let Y ( t ) be the solution of the stochastic equation (3.39) dY ( t ) = F ( Y − ( t ) , T − Ψ ( t )) dt + σ ( Y − ( t ) , T − Ψ ( t )) dB ( t ) + h ( Y − ( t ) , T − Ψ ( t )) dL ( t ) , t ≥ ,Y (0) = 0 , T Ψ (0) = 0 where the function F ( x, t ) , σ ( x, t ) , h ( x, t ) ∈ C ( R ) satisfy the Lipschitz condition. As-sume that the PDF of the process ( Y ( t ) , T Ψ ( t )) , p t ( y, z ) exists. Furthermore, we assume that ∂∂t p t ( y, z ) , ∂∂y p t ( y, z ) , ∂ ∂y p t ( y, z ) exist, (3.40) Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂t p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdt < ∞ for each t > , (3.41) Z x x Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂x p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdx < ∞ , Z x x Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ ∂x p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdx < ∞ , for each x , x ∈ R and (3.42) Z ∞ Z R −{ } | p s ( x + r, t ) − p s ( x, t ) | ν ′ ( dr ) ds < ∞ . for each x ∈ R , where ν ′ ( B ) = ν ( { x ; xh ( y, z ) ∈ B } ) .Then, the PDF of the process X ( T ) = Y − ( S Ψ ( t )) is the weak solution of fractional Fokker-Planck equation (3.43) ∂q ( x, t ) ∂t = (cid:20) − ∂∂x F ( x, t ) + 12 ∂ ∂x σ ( x, t ) + T s (cid:21) Φ t q ( x, t ) , with q ( x,
0) = δ ( x ) , where (3.44) T s f ( x, t ) = Z R −{ } [ f ( x + r, t ) − f ( x, t )] ν ′ ( dr ) , for any f ( x, t ) ∈ C ( R ) .Proof. We follow the same steps as in the proof of Theorem 3.2, and modify the part relatedto T α . Let L ( t ) be a symmetric L´evy process as defined in definition 3.3, it has self-adjointinfinitesimal generator(3.45) ( T f )( y ) = Z R −{ } [ f ( y + x ) − f ( y )] ν ( dx )for each f ∈ C ( R ), and L´evy symbol(3.46) η l ( u ) = Z R −{ } [ cos ( u, x ) − ν ( dx ) . Define L h ( t ) = L ( t ) h ( Y ( t ) , Z ( t )), since theorem 3.4 involves h ( Y ( t ) , T Ψ ( t )) dL ( t ) instead ofjust dL ( t ), by Proposition 11.10 of [26] and Corollary 3.4.11 of [1], L h ( t ) has L´evy symbol(3.47) η l h ( u ) = Z R −{ } [ cos ( u, x ) − ν ′ ( dx ) , also ν ′ is symmetric, hence L h has self-adjoint generator(3.48) T s f ( y, z ) = Z R −{ } [ f ( y + x, z ) − f ( y, z )] ν ′ ( dx ) . t follows that(3.49) Z R p t ( y, z ) T s f ( y, z ) dydz = Z R f ( y, z ) T s p t ( y, z ) dydz = Z R f ( y, z ) Z R −{ } [ p t ( y + x, z ) − p t ( y, z )] ν ′ ( dx ) dydz, so,(3.50) T + s p t ( y, z ) = Z R −{ } [ p t ( y + x, z ) − p t ( y, z )] ν ′ ( dx )Since Γ + = A + + T s + and ∂∂t p t ( y, z ) = Γ + p t ( y, z ), we have(3.51) ∂∂t p t ( y, z ) = Γ + p t ( y, z )= A + p t ( y, z ) + T s + p t ( y, z )= − ∂∂y ( F ( y, z ) p t ( y, z )) + ∂ ∂y ( 12 σ ( y, z ) p t ( y, z )) − ∂∂z Θ z p t ( y, z )+ Z R −{ } [ p t ( y + x, z ) − p t ( y, z )] ν ′ ( dx )Similar as proof of Theorem 3.2, we have (3.28) and (3.29) , plus lim s →∞ p s ( x, t ) = 0 and p ( x, t ) = (0 , ( x, t ),(3.52) ∂∂t q ( x, t ) = Z ∞ (cid:20) − ∂∂x ( F ( x, t ) p s ( x, t )) + ∂ ∂x ( 12 σ ( x, t ) p s ( x, t )) − ∂∂s p s ( x, t ) (cid:21) ds + Z ∞ Z R −{ } [ p s ( x + r, t ) − p s ( x, t )] ν ′ ( dr ) ds = − ∂∂x (cid:20) F ( x, t ) Z ∞ p s ( x, t ) ds (cid:21) + ∂ ∂x (cid:20) σ ( x, t ) Z ∞ p s ( x, t ) ds (cid:21) + Z R −{ } (cid:20)Z ∞ p s ( x + r, t ) ds − Z ∞ p s ( x, t ) ds (cid:21) ν ′ ( dr )= − ∂∂x (cid:20) F ( x, t ) Z ∞ p s ( x, t )) ds (cid:21) + ∂ ∂x (cid:20) σ ( x, t ) Z ∞ p s ( x, t ) ds (cid:21) + T s Z ∞ p s ( x, t ) ds = − ∂∂x [ F ( x, t )Φ t w ( x, t )] + ∂ ∂x (cid:20) σ ( x, t )Φ t q ( x, t ) (cid:21) + T s Φ t q ( x, t )= (cid:20) − ∂∂x F ( x, t ) + 12 ∂ ∂x σ ( x, t ) + T s (cid:21) Φ t q ( x, t ) (cid:3) Remark 3.5. (1.17) is a special case of (3.43) when L´evy measure in Theorem 3.4 is definedas in Definition 3.1.
Fractional Fokker-Planck equation with a general L´evy generator.
Now we solvefractional Fokker-Planck equation associated with a general L´evy process.
Definition 3.6.
Let L ( t ) be a L ´ evy process satisfying the following SDE (3.53) dL ( t ) = Z B x ˜ N ( dt, dx ) + Z B c xN ( dt, dx ) ith L´evy symbol η ( u ) = R R [ e iux − − iux B ( x )] ν ( dx ) , where ν is a L ´ evy measure. Theorem 3.7.
Suppose that the standard Brownian motion B ( t ) , the L´evy process L ( t ) as de-fined in definition 3.6 and subordinator T Ψ ( t ) are independent, where T Ψ ( t ) has Laplace exponent Ψ( u ) and its inverse S Ψ ( t ) . Let Y ( t ) be the solution of the stochastic equation (3.54) dY ( t ) = F ( Y − ( t ) , T − Ψ ( t )) dt + σ ( Y − ( t ) , T − Ψ ( t )) dB ( t ) + h ( Y − ( t ) , T − Ψ ( t )) dL ( t ) , t ≥ ,Y (0) = 0 , T Ψ (0) = 0 where the function F ( x, t ) , σ ( x, t ) ∈ C ( R ) , h ( x, t ) ∈ C ∞ ( R ) satisfy the Lipschitz condition.Assume that the PDF of the process ( Y ( t ) , T Ψ ( T )) , p t ( y, z ) exists. Furthermore, we assumethat ∂∂t p t ( y, z ) , ∂ k ∂y k p t ( y, z ) , k ∈ N + exist, (3.55) Z t Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂t p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdt < ∞ for each t > , (3.56) Z x x Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂∂x p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdx < ∞ , Z x x Z ∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ ∂x p s ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) dsdx < ∞ , for each x , x ∈ R and (3.57) Z ∞ Z R −{ } ∞ X k =1 (cid:12)(cid:12)(cid:12)(cid:12) ( − r ) k k ! ∂ k ∂x k ( p s ( x, t ) h ( x, t ) k ) (cid:12)(cid:12)(cid:12)(cid:12) ν ( dr ) ds < ∞ . for each x ∈ R .Then, the PDF of the process X ( T ) = Y − ( S Ψ ( t )) is the weak solution of fractional Fokker-Planck equation (3.58) ∂∂t q ( x, t ) = (cid:20) − ∂∂x F ( x, t ) + ∂ ∂x σ ( x, t ) + T l + (cid:21) Φ t q ( x, t ) , with q ( x,
0) = δ ( x ) , where (3.59) T l + f ( x, t ) = Z R −{ } " ∞ X k =1 ( − r ) k k ! ∂ k ∂x k ( h ( x, t ) k f ( x, t )) + B ( r, t ) r ∂∂x ( h ( x, t ) f ( x, t )) ν ( dr ) , for any f ( x, t ) ∈ C ∞ ( R ) . Remark 3.8.
The operator T l + appears naturally in Sun and Duan [30] . Theorem 2-14 in [3] guarantees the existence of density p t ( y, z ) by requiring some regularity on coefficients asfollows, F ( x, t ) , σ ( x, t ) , h ( x, t ) ∈ C b ( R ) have bounded partial derivatives from order 0 to 3, sup x | D nx ( xh ( x, t )) | ∈ L p ( µ, ν ) for all p ≥ .Proof of Theorem 3.7. We follow the same steps as in the proof of Theorem 3.2, and modifythe part related to T α .Let L = ( L ( t ) , t ≥
0) be a L´evy process as defined in definition 3.6, it has infinitesimalgenerator(3.60) T l f ( y, z ) = Z R −{ } (cid:20) f ( y + xh ( y, z ) , z ) − f ( y, z ) + B ( x, z ) xh ( y, z ) ∂∂y f ( y, z ) (cid:21) ν ( dx )By equation (32) in [30],(3.61) T l + p t ( y, z ) = Z R −{ } " ∞ X k =1 ( − x ) k k ! ∂ k ∂y k ( h ( y, z ) k p t ( y, z )) + B ( x, z ) x ∂∂y ( h ( y, z ) p t ( y, z )) ν ( dx ) ince Γ + = A + + T l + and ∂∂t p t ( y, z ) = Γ + p t ( y, z ), we have(3.62) ∂∂t p t ( y, z ) = − ∂∂y ( F ( y, z ) p t ( y, z )) + ∂ ∂y (cid:20) σ ( y, z ) p t ( y, z ) (cid:21) − ∂∂z Θ z p t ( y, z )+ Z R −{ } " ∞ X k =1 ( − x ) k k ! ∂ k ∂y k ( h ( y, z ) k p t ( y, z )) + B ( x, z ) x ∂∂y ( h ( y, z ) p t ( y, z )) ν ( dx )Similar as proof of Theorem 3.2, we can derive (3.28) and (3.29), combined with lim s →∞ p s ( x, t ) =0 and p ( x, t ) = (0 , ( x, t ),(3.63) ∂∂t w ( x, t ) = Z ∞ (cid:20) ( − ∂∂x ( F ( x, t ) p s ( x, t )) + ∂ ∂x (cid:18) σ ( x, t ) p s ( x, t ) (cid:19) − ∂∂s p s ( x, t )+ Z R −{ } ∞ X k =1 ( − r ) k k ! ∂ k ∂x k ( h ( x, t ) k p s ( x, t )) + B ( r, t ) r ∂∂x ( h ( x, t ) p s ( x, t )) ! ν ( dr ) ds = − ∂∂x (cid:20) F ( x, t ) Z ∞ p s ( x, t ) ds (cid:21) + ∂ ∂x (cid:20) σ ( x, t ) Z ∞ p s ( x, t ) ds (cid:21) + Z R −{ } " ∞ X k =1 ( − r ) k k ! ∂ k ∂x k (cid:18) h ( x, t ) k Z ∞ p s ( x, t ) ds (cid:19) + B ( r, t ) r ∂∂x (cid:18) h ( x, t ) Z ∞ p s ( x, t ) ds (cid:19) ν ( dr )= (cid:20) − ∂∂x F ( x, t ) + ∂ ∂x σ ( x, t ) (cid:21) Φ t q ( x, t )+ Z R −{ } " ∞ X k =1 ( − r ) k k ! ∂ k ∂x k ( h ( x, t ) k Φ t q ( x, t )) + B ( r, t ) r ∂∂x ( h ( x, t )Φ t q ( x, t )) ν ( dr )= (cid:20) − ∂∂x ( F ( x, t ) + ∂ ∂x σ ( x, t ) + T l (cid:21) Φ t w ( x, t ) (cid:3) Remark 3.9.
In Theorem 3.7, when the space-time-dependent coefficient h ( x, t ) of pure jumpL´evy process only depends on time t , say h ( t ) , then (3.64) T l + Φ t q ( x, t ) = Z R −{ } " ∞ X k =1 ( − r ) k k ! ∂ k ∂x k ( h ( t ) k Φ t q ( x, t )) + B ( r, t ) r ∂∂x ( h ( t )Φ t q ( x, t )) ν ( dr )= Z R −{ } " ∞ X k =1 ( − rh ( t )) k k ! ∂ k ∂x k (Φ t q ( x, t )) + B ( r, t ) rh ( t ) ∂∂x Φ t q ( x, t ) ν ( dr )= Z R −{ } (cid:20) Φ t q ( x − rh ( t ) , t ) − Φ t q ( x, t ) + B ( r, t ) rh ( t ) ∂∂x Φ t q ( x, t ) (cid:21) ν ( dr ) , thus, (3.58) becomes (3.65) ∂∂t q ( x, t ) = (cid:20) − ∂∂x ( F ( x, t ) + ∂ ∂x σ ( x, t ) (cid:21) Φ t q ( x, t ) , + Z R −{ } (cid:20) Φ t q ( x − rh ( t ) , t ) − Φ t q ( x, t ) + B ( r, t ) rh ( t ) ∂∂x Φ t q ( x, t ) (cid:21) ν ( dr ) , which corresponds (1.9) . emark 3.10. The methods that Magdziarz and Zorawik used to calculate the adjoint of in-finitesimal generator of the L´evy process when the coefficient of the L´evy noise depends only ontime is substitution and integration by parts, (2.11) in [14] . Such a method does not work whencoefficient of L´evy noise depends on both time and space; however, using the self-adjointness ofthe infinitesimall generator of symmetric L´evy process and the method in [30] for general L´evyprocess, we can figure out the adjoint operator for L´evy noise with space-time-dependent coef-ficients. Thus, Theorem 3.7 extends [14] and provides stochastic solution of fractional Fokker-Plank equation (3.58) describing subdiffusion in full generality.
Remark 3.11.
Simulations of paths of stochastic processes play an important role in applica-tions. Results in this paper provide a useful way for obtaining approximate solutions of fractionalFokker-Planck equations mentioned above. Using Monte Carlo methods based on realization of X ( t ) , our results can be used to approximate solutions of fractional Fokker-Planck equations (1.17) , (3.43) , and (3.58) , see [8] , [10] , [11] , [23] . Also our results can be used to obtain solutionof equations with particle tracking methods, see [2] , [4] , [19] . References [1] D. Applebaum, L´evy Process and Stochastic Calculus, Cambridge University Press, Cambridge, 2004.[2] D.A. Benson, M.M. Meerschaert, Simulation of chemical reaction via particle tracking: diffusion-limitedversus thermodynamic rate-limited regimes, Water Resources Research, Vol. 44 (2008), p. W12201,doi:10.1029/2008WR007111.[3] K. Bichteler and J. B. Gravereaux, Malliavin calculus for processes with jumps, Stochastic monographs v.2,Gordon and Breach Science Publishers, Glasgow, 1987.[4] P. Chakraborty, M. M. Meerschaert, and C. Y. Lim, Parameter Estimation for Fractional Transport: Aparticle tracking approach, Water Resources Research, Vol. 45 (2009), W10415, doi:10.1029/2008WR007577.[5] B.I. Henry, T.A.M. Langlands and P. Straka, Fractional Fokker-Planck equations for subdiffusion with space-and time-dependent forces, Phys. Rev. Letter., 105(2010), 170602.[6] R. Hilfer, Applications of Fractional Calculus in Physics (World Scientific, 2000).[7] J. Janczura and A. Wylomanska, Subdynamics of financial data from fractional Fokker-Planck equation,Physica Polonica B 5.40(2009): pp. 1341-1351.[8] A. Janicki and A. Weron, Simulation and Chaotic Behavior of α -Stable Stochastic Processes, Marcel Dekker,New York, 1994.[9] Y. Liu and B. Xin, Numerical solutions of a fractional predatorprey system, Advances in Difference Equa-tions, vol. 2011, Article ID 190475, 11 pages, 2011.[10] M. Magdziarz, Stochastic representation of subdiffusion processes with time-dependent drift, Stoch. Proc.Appl., 119 (2009), 3238-3252.[11] M. Magdziarz, Langevin picture of subdiffusion with infinitely divisible waiting times, J. Stat. Phys., 135(2009), 763-772.[12] M. Magdziarz, J. Gajda and T. Zorawik, Comment on Fractional FokkerPlanck Equation with Space andTime Dependent Drift and Diffusion, J. Stat. Phys. 154 (2014), 1241-1250.[13] M. Magdziarz, A. Weron and K. Weron, Fractional Fokker-Planck dynamics: Stochastic representation andcomputer simulation, Phys. Rev. E, 75 (2007), 016708.[14] M. Magdziarz and T. Zorawik, Stochastic representation of fractional subdiffusion equation. The case ofinfinitely divisible waiting times, L´evy noise and space-time-dependent coefficients. Proc. Amer. Math. Soc.,Accepted (2015).[15] M. M. Meerschaert and Alla Sikorskii, Stochastic models for fractional calculus,De Gruyter, 2012.[16] M.M. Meerschaert, E. Nane and P. Vellaisamy, Fractional Cauchy problems on bounded domains, Ann.Probab., 37 (2009), 979-1007.[17] M.M. Meerschaert, E. Nane and P. Vellaisamy, Distributed-order fractional Cauchy problems on boundeddomains, J. Math. Anal. Appl., 379 (2011), 216-228.[18] M.M. Meerschaert, E. Nane and P. Vellaisamy, Transient Anomalous Sub-diffusion on bounded domains,Proc. Amer. Math. Soc., 141 (2013), 699-710.[19] M. M. Meerschaert, Y. Zhang, and B. Baeumer Particle tracking for fractional diffusion with two timescales, Computers and Mathematics with Applications, Special Issue on Advances in Fractional DifferentialEquations, Vol. 59 (2010), No. 3, pp. 1078-1086, doi:10.1016/j.camwa.2009.05.009.[20] R. Metzler, E. Barkai and J. Klafter., Anomalous diffusion and relaxation close to thermal equilibrium: afractional Fokker-Planck equation approach. Phys. Rev. Lett.
21] R. Metzler and J.Klafter, The Random Walk’s Guide to Anomalous Diffusion: A Fractional DynamicsApproach, Phys. Rep. 339.1 (2000).[22] J. Mijena and E. Nane, Strong analytic solutions of fractional Cauchy problems, Proc. Amer. Math. Soc.,Proc. Amer. Math. Soc. 142 (2014), 1717-1731.[23] P. Protter, D. Talay, The Euler scheme for L´ e vy driven stochastic differential equations, Ann. Prob. 25,393?423 (1997).[24] D.Revuz and M.Yor, Continuous Martingales and Brownian Motion, Springer-Verlag, 1991.[25] S.G. Samko, A.A. Kilbas and D.I. Maritchev, Integral and Derivatives of the Fractional Order and Some ofTheir applications, Gordon and Breach, Amsterdam, 1993.[26] K.I. Sato, L´evy Processes and Infinitely Divisible Distributions, Cambridge University Press, Cambridge,1999.[27] H. Scher and M. Lax, Stochastic transport in a disordered solid. I. Theory, Phys, Rev. B 7, 4491-4502.[28] I. M. Sokolov and J. Klafter, Field-Induced Dispersion in Subdiffusion, Phys. Rev. Lett., 97 (2006), 140602.[29] T.H. Solomon, E.R. Weeks and H.L. Swinney, Observation of anomalous diffusion and Lvy flights in atwo-dimensional rotating flow. Phys. Rev. Lett. 71, 3975-3978 (1993).[30] X. Sun and J. Duan, Fokker-Planck equation for nonlinear dynamical systems driven by non-Gaussian L´evyprocesses, Journal of Mathematical Physics, Volume 53, Issue 7, pp. 072701-072701-10 (2012). Department of Mathematics and Statistics, Auburn University, Auburn, AL 36849 USA
E-mail address : [email protected] Department of Mathematics and Statistics, Auburn University, Auburn, AL 36849 USA
E-mail address : [email protected]@auburn.edu