Multi-dimensional fractional wave equation and some properties of its fundamental solution
MMulti-dimensional fractional wave equationand some properties of its fundamental solution
Yuri Luchko
Department of Mathematics, Physics, and ChemistryBeuth Technical University of Applied Sciences BerlinLuxemburger Str. 10, 13353 Berlin, Germanye-mail: [email protected]
Abstract
In this paper, a multi-dimensional fractional wave equation that de-scribes propagation of the damped waves is introduced and ana-lyzed. In contrast to the fractional diffusion-wave equation, the frac-tional wave equation contains fractional derivatives of the same order α, ≤ α ≤ MSC 2010 : 26A33, 35C05, 35E05, 35L05, 45K05, 60E99
Key Words : Caputo fractional derivative, Riesz fractional derivative, frac-tional wave equation, fundamental solution, Mellin transform, Mellin-Barnesintegral, phase velocity, damped waves
Fractional order differential equations have been successfully employed formodeling of many different processes and systems in physics, chemistry,engineering, medicine, biology, etc. (see e.g. [2], [3], [9]-[15], [19], [22]-[24],[27], [32]) to mention only few of many recent publications).One of the most interesting and important topics in this field thathas been considered by a number of researches since 1980’s are models1 a r X i v : . [ m a t h - ph ] N ov or anomalous transport processes in form of time- and/or space-fractionaladvection-diffusion-wave equations. Anomalous transport represents a verynatural way for description of time behavior of many complex systems thatare characterized by a large diversity of elementary particles that participatein the transport processes, by a strong interactions between them, and by ananomalous evolution of the whole system in time. In this paper, the focuswill be on anomalous diffusion and wave propagation that are described bythe fractional wave equation that is also referred to as the neutral-fractionaldiffusion equation; we do not consider here the effects of the advection orthe degradation components of the general transport processes.It is well known that whereas the diffusion equation describes a pro-cess, where a disturbance of the initial conditions spreads infinitely fast, thepropagation velocity of the disturbance is constant for the wave equation.In this sense, the time-fractional diffusion-wave equation that is obtainedfrom the diffusion equation by substituting the first derivative in time bythe fractional derivative of order α, < α <
2, interpolates between thesetwo different behaviors. Namely, whereas a response of its fundamental solu-tion to a localized disturbance spreads infinitely fast, its maximum locationdisperses with a finite velocity v ( t, α ) that is determined by the formula ([4],[18], [20]) v ( t, α ) = C α t α − . (1)For α = 1 (diffusion), the propagation velocity is equal to zero because of C = 0, for α = 2 (wave propagation) it remains constant, whereas for allintermediate values of α the propagation velocity of the maximum pointdepends on time t and is a decreasing function that varies from + ∞ at time t = 0+ to zero as t → + ∞ . This fact makes it difficult to interpret thesolutions to the time-fractional diffusion-wave equation as waves.The fractional wave equation we deal with in this paper contains frac-tional derivatives of the same order α, ≤ α ≤ α .It is important to note that the anomalous transport models are usuallyfirst formulated in stochastic form in terms of the continuous time random2alk processes. The time- and/or space-fractional differential equations arethen derived from the stochastic models for a special choice of the jumpprobability density functions with the infinite first or/and second moments(see e.g. [11, 15, 27]). The fractional wave equation can be obtained fromthe continuous time random walk processes, too. In [14], the case of thewaiting time probability density function and the jump length probabilitydensity function with the same power low asymptotic behavior has beenconsidered. Under some standard assumptions, the fractional wave equationcab be asymptocially derived from the continuous time random walk modelmentioned above (see [14] for details).From the mathematical viewpoint, the fractional wave equation was con-sidered for the first time in [6], where an explicit formula for the fundamen-tal solution of the one-dimensional fractional wave equation was derivedand where this equation was referred to as the neutral-fractional diffusionequation. In [25], a one-dimensional space-time fractional diffusion-waveequation with the Riesz-Feller derivative of order α ∈ (0 ,
2] and skewness θ and with the Caputo fractional derivative of order β ∈ (0 ,
2] was investigatedin detail. A particular case of this equation called the neutral-fractional dif-fusion equation that for θ = 0 corresponds to our one-dimensional fractionalwave equation has been shortly mentioned in [25]. In [28], a fundamen-tal solution to the one-dimensional neutral-fractional diffusion equation wasdeduced and analyzed in terms of the Fox H-function. Finally, in [12], theone-dimensional fractional wave equation was investigated in detail. Its fun-damental solution G α, was derived in terms of elementary functions for allvalues of α, ≤ α <
2. Moreover, G α, was interpreted as a spatial prob-ability density function evolving in time all whose moments of order lessthan α are finite. For the fundamental solution G α, , both its maximumlocation and its maximum value were determined in closed form. Finally,it was shown in [12] that both the maximum and the gravity and ”mass”centers of the fundamental solution G α, propagate with constant velocitieslike in the case of the wave equation, but in contrast to the wave equa-tion ( α = 2) these velocities are different each to other for a fixed value of α, < α <
2. Moreover, the first, the second, and the Smith centroveloci-ties of the damped waves described by the one-dimensional fractional waveequation were shown to be constants that depend just on the equation order α . In this paper, we deal with the multi-dimensional fractional wave equa-tion with a special focus given to the most important three-dimensionalcase. The fundamental solution to this equation is a spherically symmet-ric function that possesses nice integral representations and can be even3ritten down in explicit form in terms of elementary functions in the one-and three-dimensional cases. In contrast to the one-dimensional case, thefundamental solution for the multi-dimensional case cannot be interpretedas a probability density function and thus these equations cannot be em-ployed for modeling of any diffusion processes. Instead, we show that thefundamental solution can be interpreted as a damped wave with the con-stant phase velocity that just depends on the order α of the fractional waveequation.The rest of the paper is organized as follows. In the 2nd section, ba-sic definitions, problem formulation, and some analytical results for theinitial-value problems for the multi-dimensional fractional wave equationare presented. In particular, the fundamental solutions G α, for the one-dimensional and G α, for the three-dimensional wave equations are derivedin terms of elementary functions for all values of α, ≤ α <
2. Whereas G α, can be interpreted as a spatial probability density function evolving intime whose moments up to order α are finite, it is not true for the higherdimensions. But the interpretation of the fundamental solutions as dampedwaves is still valid for the two- and three-dimensional cases. In particular,we show that the maximum location (phase) velocity of the fundamentalsolution G α, of the three-dimensional wave equation is a constant that de-pends just on the equation order α . To illustrate analytical findings, resultsof numerical calculations, numerous plots, their physical interpretation anddiscussion are presented in the last section along with some conclusions andopen problems for further research. In this paper, we deal with the initial-value problem u ( x,
0) = ϕ ( x ) , ∂u∂t ( x,
0) = 0 , x ∈ IR n (2)for the multi-dimensional fractional wave equation D αt u ( x, t ) = − ( − ∆) α u ( x, t ) , x ∈ IR n , t ∈ IR + , ≤ α ≤ . (3)In (3), u = u ( x, t ) , x ∈ IR n , t ∈ IR + is a real field variable, − ( − ∆) α is theRiesz space-fractional derivative of order α that is defined below, and D αt isthe Caputo time-fractional derivative of order α :( D αt u )( t ) := (cid:18) I n − α ∂ n u∂t n (cid:19) ( t ) , n − < α ≤ n, n ∈ IN (4)4 α , α ≥ I α u )( t ) := (cid:40) α ) (cid:82) t ( t − τ ) α − u ( x, τ ) dτ, α > ,u ( x, t ) , α = 0and Γ the Euler Gamma function. For α = n, n ∈ IN , the Caputo fractionalderivative coincides with the standard derivative of order n .All quantities in (3) are supposed to be dimensionless, so that the coef-ficient by the Riesz space-fractional derivative can be taken to be equal toone without loss of generality.In this paper, we are mostly interested in behavior and properties of thefirst fundamental solution (Green function) G α,n of the equation (3), i.e. itssolution with the initial condition ϕ ( x ) = (cid:81) ni =1 δ ( x i ) , x ∈ IR n , δ being theDirac delta function.For a sufficiently well-behaved function f : IR n → IR , the Riesz space-fractional derivative of order α, < α ≤ −| κ | α ([30], [31]):( F − ( − ∆) α f )( κ ) = −| κ | α ( F f )( κ ) , (5)( F f )( κ ) being the Fourier transform of a function f at the point κ ∈ IR n that is defined by ( F f )( κ ) = ˆ f ( κ ) = (cid:90) IR n e iκ · x f ( x ) dx. (6)For 0 < α < α (cid:54) = 1, the Riesz space-fractional derivative can be repre-sented as a hypersingular integral − ( − ∆) α f ( x ) = − d n,l ( α ) (cid:90) IR n (∆ lh f )( x ) | h | n + α dh with the finite differences operator (∆ lh f )( x ) = (cid:80) lk =0 ( − k (cid:16) lk (cid:17) f ( x − kh )and a suitable normalization constant d n,l ( α ) (see [31] for more details). Inparticular, in the one-dimensional case we get the representation − ( − ∆) α f ( x ) = − − α ) cos( απ ) (cid:90) ∞ f ( x + ξ ) − f ( ξ ) + f ( x − ξ ) ξ α +1 dξ (7)that is valid for 0 < α < α (cid:54) = 1. For α = 1, the Riesz space-fractionalderivative can be interpreted in terms of the Hilbert transform (see e.g. [5]) − ( − ∆) f ( x ) = − π ddx (cid:90) + ∞−∞ f ( ξ ) x − ξ dξ, α = 1 can be rewritten in the form ∂u∂t ( x, t ) = − π ddx (cid:90) + ∞−∞ u ( ξ, t ) x − ξ dξ (8)that is of course not the standard advection equation.For α = 2, equation (3) is reduced to the multi-dimensional wave equa-tion. In what follows, we focus on the case 1 ≤ α < α = 2 (wave equation) is well studied in the literature. We start our analysis by applying the multi-dimensional Fourier integraltransform (6) with respect to the variable x ∈ IR n to the fractional waveequation (3) with 1 < α < ϕ ( x ) = (cid:81) ni =1 δ ( x i ). Using the definition of the Riesz fractional derivative, for theFourier transform ˆ G α,n ( κ, t ) of the fundamental solution G α,n ( x, t ) we getthe initial-value problem (cid:40) ˆ G α,n ( κ,
0) = 1 , ∂ ˆ G α,n ∂t ( κ,
0) = 0 (9)for the fractional differential equation( D αt ˆ G α,n ( κ, t ))( t ) + | κ | α ˆ G α,n ( κ, t ) = 0 . (10)The unique solution of the problem (9)-(10) is given by the expression (seee.g. [16]) ˆ G α,n ( κ, t ) = E α ( −| κ | α t α ) (11)in terms of the Mittag-Leffler function E α ( z ) = ∞ (cid:88) k =0 z k Γ(1 + αk ) , α > . (12)It follows from (11) and the well known asymptotic formula E α ( − x ) = − m (cid:88) k =1 ( − x ) − k Γ(1 − αk ) + O ( | x | − − m ) , m ∈ IN, x → + ∞ , < α < G α,n belongs to the functional space L ( IR n ) with respect to κ for 1 < α <
2. Therefore we can apply the inverseFourier transform and get the representation G α,n ( x, t ) = 1(2 π ) n (cid:90) IR n e − iκ · x E α ( −| κ | α t α ) dκ, x ∈ IR n , t > G α,n . But E α ( −| κ | α t α ) is a radial function(spherically symmetric function) in κ and thus the known formula (see e.g.[31]) 1(2 π ) n (cid:90) IR n e − iκ · x φ ( | κ | ) dκ = | x | − n (2 π ) n (cid:90) ∞ φ ( τ ) τ n J n − ( τ | x | ) dτ (15)for the inverse Fourier transform of the radial functions can be employedunder the condition that the integral at the right hand side of the formula(15) is conditionally or absolutely convergent. In this formula, J ν stays forthe Bessel function with the index ν .The representation (14) can be thus rewritten in the form G α, n ( x, t ) = | x | − n (2 π ) n (cid:90) ∞ E α ( − τ α t α ) τ n J n − ( τ | x | ) dτ (16)under the condition that the integral at the right hand side of the formulais conditionally or absolutely convergent.Before we start with a convergence investigation of the integral in (16),let us separately consider the case x = (0 , , . . . ,
0) = that corresponds tothe case | x | = 0. It follows from (14) that G α,n ( , t ) = 1(2 π ) n (cid:90) IR n E α ( −| κ | α t α ) dκ, t > . Switching to the spherical coordinates, we can represent the last integral inthe form (cid:90) IR n E α ( −| κ | α t α ) dκ = (cid:90) ∞ (cid:90) π . . . (cid:90) π (cid:90) π E α ( − r α t α ) r n − sin n − ( φ ) · · · sin( φ n − ) dr dφ . . . dφ n − . Because of the known formula (see e.g. [29]) (cid:90) π sin α − ( φ ) dφ = √ π Γ (cid:0) α (cid:1) Γ (cid:0) α +12 (cid:1) ,
7e first get the expression G α,n ( , t ) = 12 n − π n Γ (cid:0) n (cid:1) (cid:90) ∞ E α ( − r α t α ) r n − dr (17)that can be interpreted as the Mellin integral transform of the Mittag-Lefflerfunction (see e.g. [17]) and thus evaluated: G α,n ( , t ) = 1 αt n n − π n Γ (cid:0) n (cid:1) Γ (cid:0) nα (cid:1) Γ (cid:0) − nα (cid:1) Γ (1 − n ) . Let us emphasize the convergence condition 0 < n < α for the right-handside integral in (17). For 1 < α <
2, this condition means that the funda-mental solution G α,n ( , t ) is finite just in the one-dimensional case n = 1.Because the Gamma function Γ( z ) has a pole at the point z = 0, in theone-dimensional case we get the formula G α, (0 , t ) = 0that is in accordance with the results presented in [12]. In all other cases,the fundamental solution of the fractional wave equation (3) is not finite atthe point x = (0 , , . . . ,
0) = for t > x (cid:54) = , we employ now the asymptotics (13) of the Mittag-Lefflerfunction and the known asymptotics of the Bessel function J ν ( r ) = O ( r ν ) , as r → , (18) J ν ( r ) = (cid:114) πr cos( r − π (1 + 2 ν ) /
4) + O ( r − ) , as r → + ∞ (19)to derive the condition n < α + 1 for the conditional and the condition n < α − < α <
2, these conditions means thatthis integral is conditionally convergent at least for n = 1 , , n = 1.Thus for 1 < α < n = 1 , , J ν ( z ) = (cid:16) z (cid:17) ν W ,ν +1 (cid:18) − z (cid:19) J ν and the Wright function W α,β that is definedas the convergent series W α,β ( z ) = ∞ (cid:88) m =0 z m m ! Γ( αm + β ) , α > − , β ∈ IR, we get the representation G α, n ( x, t ) = 2(4 π ) n/ (cid:90) ∞ τ n − E α ( − τ α t α ) W ,n/ − (cid:18) − τ | x | (cid:19) dτ of the fundamental solution G α, n in terms of two most known and usedspecial functions of fractional calculus - the Mittag-Leffler function E α andthe Wright function W α,β .Now we derive another important representation of the fundamental so-lution G α,n in terms of a Mellin-Barnes integral.For x (cid:54) = 0, the right-hand side of the representation G α, n ( x, t ) = (2 π ) − n/ | x | − n/ (cid:90) ∞ τ n/ E α ( − τ α t α ) J n/ − ( τ | x | ) dτ can be interpreted as the Mellin convolution of the functions E α ( − t α τ α ) and τ − n/ − J n/ − (1 /τ ) at the point y = 1 /x .Using the Mellin integral transform technique (Mellin transforms of theMittag-Leffler and the Bessel functions, some elementary properties of theMellin transform, and the Mellin convolution theorem) we get the represen-tation ( x (cid:54) = 0) G α,n ( x, t ) = 1 απ n/ | x | n πi (cid:90) L Γ (cid:0) sα (cid:1) Γ (cid:0) − sα (cid:1) Γ (cid:0) n − s (cid:1) Γ(1 − s )2 s Γ (cid:0) s (cid:1) (cid:18) t | x | (cid:19) − s ds (20)of the fundamental solution G α,n in terms of the Mellin-Barnes integral (FoxH-function). For the details regarding evaluation of the improper integralsof the Mellin convolution type we refer the interested readers to the recentpaper [17] or to the book [26].Because of the factor Γ (cid:0) n − s (cid:1) , the Mellin-Barnes representation of G α,n has a different structure for the even and the odd dimensions.In particular, for n = 3 the known formula Γ(1 + z ) = z Γ( z ) for theGamma function leads to the relationΓ (cid:0) sα (cid:1) Γ (cid:0) − sα (cid:1) Γ (cid:0) − s (cid:1) Γ(1 − s )2 s Γ (cid:0) s (cid:1) = (cid:18) − s (cid:19) Γ (cid:0) sα (cid:1) Γ (cid:0) − sα (cid:1) Γ (cid:0) − s (cid:1) Γ(1 − s )2 s Γ (cid:0) s (cid:1) G α, withthe kernel of G α, . This formula along with some elementary properties ofthe Mellin integral transform (see [17] or [26]) induces the relation G α, ( x, t ) = 12 π | x | (cid:18) G α, ( x, t ) + t ∂∂t G α, ( x, t ) (cid:19) , x (cid:54) = 0 (21)between the three-dimensional fundamental solution G α, and the one-dimensional fundamental solution G α, . In the next subsection, we deduceanother relation between these two fundamental solutions via the spatialdifferentiation. In the one-dimensional case ( n = 1), the representation (16) can be rewrittenin the form G α, ( x, t ) = 1 π (cid:90) ∞ E α ( − τ α t α ) cos( τ | x | ) dτ, x ∈ IR, t > J − / ( x ) = (cid:114) πz cos( z ) . The formula (22) is well known; it was used in [6] and then in more detailedform in [12] to get an explicit form of the fundamental solution in terms ofelementary functions: G α, ( x, t ) = 1 π | x | α − t α sin( πα/ t α + 2 | x | α t α cos( πα/
2) + | x | α , t > , x ∈ IR, ≤ α < . (23)For α = 1 (modified advection equation (8)), the fundamental solution G α, is the well known Cauchy kernel G , ( x, t ) = 1 π tt + x (24)that is a spatial probability density function evolving in time.For α = 2 (wave equation), the Green function G , is known to be givenby the formula G , ( x, t ) = 12 ( δ ( x − t ) + δ ( x + t )) . (25)10or 1 < α <
2, the Green function (23) is a spatial probability densityfunction evolving in time, too, as has been shown for the first time in [12].This pdf possesses all finite moments up to the order α . In particular, themean value of G α, (its first moment) exists for all α > β, | β | < α of G α, can be evaluated via the Mellin integral transform (see[12] for details): (cid:90) ∞ G α, ( r, t ) r β dr = t β α sin( πβ/ πβ/α ) . (26)We mention here that in [12] the locations of extrema points, gravity andmass centers of G α, , and location of its energy center were calculated inexplicit form and plotted for different values of the equation order α .Another interesting and important feature of the one-dimensional fun-damental solution G α, is that along with its probabilistic interpretation itcan be interpreted as a damped wave, too. In [12], different velocities ofthis wave were calculated and plotted depending on α . In particular, thepropagation velocity of its maximum that can be interpreted as the phasevelocity, the propagation velocity of its gravity center, the velocity of its”mass”-center, its pulse velocity, and three different kinds of its centrove-locities have been introduced and calculated. It turned out that all thesevelocities are constant in time and depend just on the order α of the frac-tional wave equation.The plots of the propagation velocity v p of the maximum location ofthe fundamental solution G α, (phase velocity), the velocity v g of its gravitycenter, its pulse velocity v m and its centrovelocity v c are presented in Fig.1. As expected, v p = v c = 0, v m = π ≈ .
64 for α = 1 (modified convectionequation), and all velocities smoothly approach the value 1 as α → < α < v p , v m , and v c monotonously increase whereas v g monotonously decreases. It is interesting to note that all velocities of G α are nearly the same as those of the fundamental solution of the waveequation in a small neighborhood of the point α = 2. The velocity v g ofthe gravity center of G α tends to + ∞ for α → t > α, < α < v p , v g , v m , v c are different each to other and fulfill the inequalities v c ( α ) < v p ( α ) < v m ( α ) < v g ( α ). For α = 2, all velocities are equal to 1.In the two-dimensional case ( n = 2), the representation (16) has the11 .0 1.2 1.4 1.6 1.8 2.0 Alpha0.00.51.01.52.02.53.0 V V_cV_pV_mV_g
Figure 1: Plots of the gravity center velocity v g ( α ), the pulse velocity v m ( α ),the phase velocity v p ( α ), and the centrovelocity v c ( α ) for 1 ≤ α ≤ G α, ( x, t ) = 12 π (cid:90) ∞ τ E α ( − τ α t α ) J ( τ | x | ) dτ, x (cid:54) = 0 , t > , (27)where the Bessel function with the index 0 can be represented as e.g. J ( z ) = (cid:90) ∞ cos( z sin( φ )) dφ. The representation (27) can be used for e.g. numerical evaluation of thefundamental solution G α, . The Mittag-Leffler function E α is evaluated byemploying the algorithms suggested in [7] and the MATLAB-programs [33]that implement these algorithms. Some results of the numerical evaluationsof G α, are presented in the next section. Say, the plot presented in Fig. 3clearly shows that the fundamental solution G α, is negative for some valuesof the variables x and t and therefore cannot be interpreted as a probabilitydensity function.Now we focus on the three-dimensional case. For x (cid:54) = 0 , t >
0, therepresentation (16) takes the form G α, ( x, t ) = (2 π ) − / | x | − / (cid:90) ∞ τ / E α ( − τ α t α ) J / ( τ | x | ) dτ (28)that can be represented as G α, ( x, t ) = 12 π | x | (cid:90) ∞ E α ( − τ α t α ) τ sin( τ x ) dτ (29)12ecause of the formula J / ( x ) = (cid:114) πz sin( z ) . Comparing the representations (22) and (29), we easily get the relation( x (cid:54) = 0) G α, ( x, t ) = − π | x | ∂∂ | x | G α, ( x, t ) (30)between the fundamental solutions of the one-dimensional and three-dimensional fractional wave equations. It is worth to mention that the samerelation has been deduced in [8] for the time-space-fractional diffusion equa-tion with the time derivative of order α, < α < β, < β < β (cid:54) = 1 by using a different method.The representation (23) of the fundamental solution G α, along with therelation (21) or the relation (30) between the fundamental solutions G α, and G α, leads to a closed form formula for fundamental solution of thethree-dimensional fractional wave equation ( x (cid:54) = 0 , t > G α, ( x, t ) = sin( πα/ π (cid:0) − ( α − t α + 2 | x | α t α cos( πα/
2) + (1 + α ) | x | α (cid:1) | x | − α t − α ( t α + 2 | x | α t α cos( πα/
2) + | x | α ) . (31)It follows from this formula that G α, is not a pdf because it is not everywherenonnegative. More precisely, the following inequalities are valid with | x | = r > G α, ( r, t ) < r < r ∗ ( α, t ) ,G α, ( r, t ) = 0 for r = r ∗ ( α, t ) ,G α, ( r, t ) > r > r ∗ ( α, t ) , where r ∗ ( α, t ) = z α t, z α = (cid:32) − cos( πα/
2) + (cid:112) α − sin ( πα/ α + 1 (cid:33) α . (32)Because of the relation (30) between the fundamental solutions G α, and G α, , the point r ∗ ( α, t ) = z α t is the only maximum location of the funda-mental solution G α, ( r, t ) with r > , t >
0. In [12], explicit formulas forthe location of the maximum point of G α, , its maximum values and themonotonicity intervals were derived. The inequalities for the fundamentalsolution G α, follow from these results.13rom the closed form formula (31), the asymptotics of G α, ( r, t ) , r = | x | > t, t > α, < α < G α, ( r, t ) = O ( r α − ) , r → ,G α, ( r, t ) = O ( r − α − ) , r → + ∞ . Thus the integral (”moment” of G α, ( r, t ) of the order β ) I α,β ( t ) = (cid:90) ∞ r β G α, ( r, t ) dr exists for 2 − α < β < α, i.e., at least for 1 ≤ β ≤ < α < . Employing the the Mellin-Barnes representation of G α, ( x (cid:54) = 0) in form G α, ( x, t ) = 1 απ / | x | πi (cid:90) L Γ (cid:0) sα (cid:1) Γ (cid:0) − sα (cid:1) Γ (cid:0) − s (cid:1) Γ(1 − s )2 s Γ (cid:0) s (cid:1) (cid:18) t | x | (cid:19) − s ds that is a particular case of the formula (20) and that can be interpreted asthe inverse Mellin integral transform of its kernel, we get an explicit formulafor the Mellin integral transform of G α, and thus its ”moments” of the order β, ≤ β ≤ I α,β ( t ) = (2 t ) β − απ / Γ (cid:16) − βα (cid:17) Γ (cid:16) − − βα (cid:17) Γ (cid:16) β +12 (cid:17) Γ( β − (cid:16) − β (cid:17) . Using the known properties of the Gamma function, a simpler representationfor the ”moments” in form I α,β ( t ) = t β − ( β − απ sin (cid:16) π β (cid:17) sin (cid:16) π (2 − β ) α (cid:17) can be derived. In particular, we get the following important particularcases:1) β = 1 (”‘mean value”’): I α, ( t ) ≡ , for all 1 < α < , t > , β = 2 (2nd ”moment”): I α, ( t ) ≡ π , for all 1 < α < , t > , β = 3 (3rd ”moment”): I α, ( t ) ≡ tαπ sin( π/α ) , for all 1 < α < , t > . It follows from the Mellin-Barnes representation (20) that G α,n ( x, t ) can berepresented via an auxiliary function that depends on a single argument: G α,n ( x, t ) = G α,n ( | x | , t ) = G α,n ( r, t ) = r − n L α,n (cid:16) rt (cid:17) , r > L α,n ( r ) = 1 απ n/ πi (cid:90) L Γ (cid:0) sα (cid:1) Γ (cid:0) − sα (cid:1) Γ (cid:0) n − s (cid:1) Γ(1 − s )2 s Γ (cid:0) s (cid:1) ( r ) s ds. (34)Let the function r − n L α,n ( r ) attain its local maximum or minimum at thepoint r ∗ = c ( α, n ). Then the fundamental solution G α,n can be representedas G α,n ( r, t ) = t − n (cid:16) rt (cid:17) − n L α,n (cid:16) rt (cid:17) and therefore for a fixed t > G α,n ( r, t ) attains its local maximum or mini-mum at the point r ∗ ( t, α, n ) = r ∗ t = c ( α, n ) t. (35)For n = 1, the function r − n L α,n ( r ) is connected with the stable unimodaldistributions and is unimodal, i.e., it has only one maximum point thatmeans that the fundamental solution G α, of the one-dimensional fractionalwave equation has only one maximum point (see [12] for details). In thetwo-dimensional case, the fundamental solution G α, has many (probablyinfinitely many) local minima and maxima as can be seen on the plot of Fig.3. Finally, the fundamental solution G α, has only one maximum point ascan be seen on the plots of Fig. 4 and is proved by direct calculations basedon the closed form formula (31).The value G ∗ α,n ( t ) of the fundamental solution G α,n ( r, t ) at a local min-imum or maximum point given by (35) is equal to G ∗ α,n ( t ) = ( c ( α, n ) t ) − n L α,n ( c ( α, n )) = 1 t n ( c ( α, n )) − n L α,n ( c ( α, n )) , G . t=0.1t=0.2 t=0.3 Figure 2: Plots of G α, ( x, t ) with α = 1 . t where the function L α,n is defined as in (33) and (34) and thus the product( r ∗ ( t, α, n )) n G ∗ α,n ( t ) = L α,n ( c ( α, n ))is time independent. In particular, for n = 1 we get the relation r ∗ ( t, α, G ∗ α, ( t ) = L α, ( c ( α, , i.e., the product of the maximum location and the maximum value of thefundamental solution G α, is time independent and depends just on α . Forthe numerical calculations and plots of the maximum location, maximumvalue and their product we refer the interested reader to [12]. In this section, some results of numerical evaluations of the fundamentalsolutions G α,n , n = 1 , ,
3, their plots, discussions, and open problems arepresented.
We start with some plots of the fundamental solution G α, that can beeasily obtained using the closed form formula (23) and are shown in Fig.2. For more results regarding G α, , its probabilistic and wave propagationinterpretations we refer to [12].On Fig. 3, a plot of the fundamental solution G α, of the two-dimensionalfractional wave equation is presented. To produce the plot, we employed theintegral representation (27) and the MATLAB-programs [33] for numericalevaluation of the Mittag-Leffler function E α that implement the algorithms16 .2 0.4 0.6 0.8 1.0 x (cid:45) (cid:45) (cid:45) Figure 3: Plots of G α, ( r, t ) , r = | x | > α = 1 . t = 1suggested in [7]. An important information that can be read from the ploton Fig. 3 is that the fundamental solution G α, is negative for some valuesof the variables x and t and therefore cannot be interpreted as a probabilitydensity function. Another important point is that G α, is not unimodal andhas many (probably infinitely many) local minimum and maximum points.On the next Fig. 4, some plots of the fundamental solution G α, of thethree-dimensional fractional wave equation are presented.As we can see on the plots, for a fixed value of t , the fundamental solution G α, ( r, t ) has only one maximum point that depends on the time instant t (and of course on the value of the parameter α ).The velocity of the maximum location of G α, ( r, t ) (its phase velocity)can be calculated based on the formula (35): v p ( α ) := dr ∗ ( t, α, dt = d ( c ( α, t ) dt = c ( α, , where c ( α,
3) is the location of the maximum point of the fundamental so-lution G α, ( r, t ) at the time instant t = 1. As we see, the phase velocityis time independent. It can be numerically calculated based on the closedform formula (31) for the fundamental solution G α, ( r, t ). The results ofnumerical calculations are presented in Fig. 5.It is interesting to note that the phase velocity v p ( α ) is not monotone in α and attains a maximum at the point α ≈ . α , α , α (cid:54) = α that fulfill the equality v p ( α ) = v p ( α ), i.e., the corresponding damped waves propagate with thesame phase velocity. 17 = 0.2t=0.3 t=0.4 x (cid:45) (cid:45) Figure 4: Plots of G α, ( r, t ) with α = 1 . t =0 . , . , , a=1.575, max =1.12 Figure 5: Plot of the phase velocity v p ( α ) of G α, ( r, t ) , r = | x | = 0.3r = 0.5 r = 0.7 t (cid:45) Figure 6: Plots of the fundamental solution G α, ( r, t ) , r = | x | > α = 1 . r = 0 . , . , , G α, with thefixed α = 1 . | x | = r are presented.As we see on the plots, the profiles of the fundamental solution G α, look like those of the damped waves that justifies our interpretation of theequation (3) as the fractional wave equation. In this paper, we discussed different representations and properties of thefundamental solution G α,n of the fractional wave equation and interpreted itas a damped wave. In particular, the maximum location of the fundamentalsolution and its velocity were determined. The maximum location velocitythat can be interpreted as the phase velocity was shown to be a constantthat depends just on the order α of the equation.Among open questions for further research, investigation of other veloci-ties of the damped waves that are described by the fractional wave equationlike the velocity of the mass-center, the pulse velocity, the group velocity,the centrovelocities, etc. should be first mentioned. Of course, the fractionalwave equations with fractional derivatives of other types and with the non-constant coefficients should be investigated, too. In particular, the fractionalwave equation with the Riesz-Feller derivative of order α ∈ (0 ,
2] and skew-ness θ would be for sure an interesting research object. We note here that in[25] an one-dimensional space-time fractional diffusion-wave equation withthe Caputo derivative of order β ∈ (0 ,
2] in time and the Riesz-Feller deriva-19ive of order α ∈ (0 ,
2] and skewness θ in space has been investigated indetail. A particular case of this equation called neutral-fractional diffu-sion equation that for θ = 0 corresponds to our one-dimensional fractionalwave equation was introduced in [25]. It would be interesting to considera multi-dimensional analog of the neutral-fractional diffusion equation thatcorresponds to the one-dimensional case considered in [25].Another potentially interesting research direction would be to employ thefractional wave equation introduced in this paper in the theory of the frac-tional Schr¨odinger equation (see e.g. [1] or [13] for more details). Until now,different kinds of the space-, time, and space-time-fractional Schr¨odingerequations were introduced and analyzed, but not a fractional Schr¨odingerequation that corresponds to our fractional wave equation.Finally, we mention a very recent research field of fractional calculusthat deals with the inverse problems for the fractional differential equations(see e.g. [21] and references there). It would be very interesting to considersome inverse problems for the multi-dimensional fractional wave equationintroduced in this paper and their applications. References [1] B. Al-Saqabi, L. Boyadjiev, and Yu. Luchko, Comments on employingthe Riesz-Feller derivative in the Schr¨odinger equation.
The EuropeanPhysical Journal Special Topics (2013), 1779-1794.[2] J.L.A. Dubbeldam, A. Milchev, V.G. Rostiashvili, and T.A. Vilgis,Polymer translocation through a nanopore: A showcase of anomalousdiffusion.
Physical Review E (2007), 010801 (R).[3] A. Freed, K. Diethelm, and Yu. Luchko, Fractional-order viscoelasticity(FOV): Constitutive development using the fractional calculus . NASA’sGlenn Research Center, Ohio (2002).[4] Y. Fujita, Integrodifferential equation which interpolates the heat equa-tion and the wave equation, I, II.
Osaka J. Math. (1990), 309–321,797–804.[5] R. Gorenflo and F. Mainardi, Random walk models for space-fractionaldiffusion processes. Fract. Calc. Appl. Anal. (1998), 167-191.[6] R. Gorenflo, A. Iskenderov, and Yu. Luchko, Mapping between solutionsof fractional diffusion-wave equations. Fract. Calc. Appl. Anal. (2000),75-86. 207] R. Gorenflo, J. Loutchko, and Yu. Luchko, Computation of the Mittag-Leffler function and its derivatives. Fract. Calc. Appl. Anal. (2002),491-518.[8] A. Hanyga, Multi-dimensional solutions of space-time-fractional diffu-sion equations. Proc. R. Soc. London. A (2002), 429-450.[9] R. Herrmann,
Fractional Calculus: An Introduction for Physicists .World Scientific, Singapore (2011).[10] R. Hilfer (Ed.),
Applications of Fractional Calculus in Physics . WorldScientific, Singapore (2000).[11] R. Klages, G. Radons, and I.M. Sokolov (Eds.),
Anomalous Transport:Foundations and Applications , Wiley-VCH, Weinheim (2008).[12] Yu. Luchko, Fractional wave equation and damped waves.
J. Math.Phys. (2013), 031505.[13] Yu. Luchko, Fractional Schr¨odinger equation for a particle moving in apotential well. J. Math. Phys. (2013), 012111.[14] Yu. Luchko, Models of the neutral-fractional anomalous diffusionand their analysis. In: AIP Conf. Proc. 1493 (2012), 626-632, doi:http://dx.doi.org/10.1063/1.4765552.[15] Yu. Luchko, Anomalous diffusion models and their analysis. Forum derBerliner mathematischen Gesellschaft (2011), 53-85.[16] Yu. Luchko, Operational method in fractional calculus. Fract. Calc.Appl. Anal. (1999), 463-489.[17] Yu. Luchko and V. Kiryakova, The Mellin integral transform in frac-tional calculus. Fract. Calc. Appl. Anal. (2013), 405-430.[18] Yu. Luchko and F. Mainardi, Some properties of the fundamental solu-tion to the signalling problem for the fractional diffusion-wave equation. Central European Journal of Physics (2013).[19] Yu. Luchko and A. Punzi, Modeling anomalous heat transport ingeothermal reservoirs via fractional diffusion equations. InternationalJournal on Geomathematics (2011), 257-276.2120] Yu. Luchko, F. Mainardi, and Yu. Povstenko, Propagation speed of themaximum of the fundamental solution to the fractional diffusion-waveequation. Computers and Mathematics with Applications (2013),774-784.[21] Yu. Luchko, W. Rundell, M. Yamamoto and L. Zuo, Uniqueness andreconstruction of an unknown semilinear term in a time-fractionalreaction-diffusion equation. Inverse Problems (2013), 065019.[22] R.L. Magin, Fractional Calculus in Bioengineering: Part1, Part 2 andPart 3. Critical Reviews in Biomedical Engineering (2004), 1-104,105-193, 195-377.[23] F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena. Chaos, Solitons and Fractals (1996), 1461-1477.[24] F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity .Imperial College Press, London (2010).[25] F. Mainardi, Yu. Luchko, and G. Pagnini, The fundamental solutionof the space-time fractional diffusion equation.
Fract. Calc. Appl. Anal. (2001), 153-192. [E-print http://arxiv.org/abs/cond-mat/0702419].[26] O.I. Marichev, Handbook of Integral Transforms of Higher Transcen-dental Functions, Theory and Algorithmic Tables . Chichester, Ellis Hor-wood (1983).[27] R. Metzler and J. Klafter, The restaurant at the end of the randomwalk: Recent developments in the description of anomalous transportby fractional dynamics.
J. Phys. A. Math. Gen. (2004), R161-R208.[28] R. Metzler and T. F. Nonnenmacher, Space- and time-fractional diffu-sion and wave equations, fractional Fokker-Planck equations, and phys-ical motivation. Chem. Phys. (2002), 67-90.[29] A.P. Prudnikov, Yu.A. Brychkov, and O.I. Marichev,
Integrals and Se-ries, Vol. 1: Elementary Functions . Gordon and Breach, New York(1986).[30] A. Saichev, G. Zaslavsky, Fractional kinetic equations: solutions andapplications.