Space-time fractional Zener wave equation.

Downloaded from http://rspa.royalsocietypublishing.org/ on March 30, 2015

rspa.royalsocietypublishing.org

Space–time fractional Zener wave equation T. M. Atanackovic1 , M. Janev2 , Lj. Oparnica3 ,

Research Cite this article: Atanackovic TM, Janev M, Oparnica Lj, Pilipovic S, Zorica D. 2015 Space–time fractional Zener wave equation. Proc. R. Soc. A 471: 20140614. http://dx.doi.org/10.1098/rspa.2014.0614 Received: 13 August 2014 Accepted: 1 December 2014

Subject Areas: mechanics, wave motion, mathematical physics Keywords: fractional Zener model, fractional strain measure, Laplace and Fourier transforms, Cauchy problem, generalized solution Author for correspondence: D. Zorica e-mail: [email protected]

S. Pilipovic4 and D. Zorica2 1 Department of Mechanics, Faculty of Technical Sciences, University

of Novi Sad, Trg D. Obradovica 6, 21000 Novi Sad, Serbia 2 Mathematical Institute, Serbian Academy of Arts and Sciences, Kneza Mihaila 36, 11000 Belgrade, Serbia 3 Faculty of Education in Sombor, University of Novi Sad, Podgorička 4, 25000 Sombor, Serbia 4 Department of Mathematics, Faculty of Natural Sciences and Mathematics, University of Novi Sad, Trg D. Obradovica 4, 21000 Novi Sad, Serbia DZ, 0000-0002-9117-8589 The space–time fractional Zener wave equation, describing viscoelastic materials obeying the timefractional Zener model and the space-fractional strain measure, is derived and analysed. This model includes waves with finite speed, as well as nonpropagating disturbances. The existence and the uniqueness of the solution to the generalized Cauchy problem are proved. Special cases are investigated and numerical examples are presented.

1. Introduction Here, we study a class of generalized wave equations. The wave equation can be generalized within the theory of fractional calculus by replacing the second-order derivative (space and/or time) with the fractional ones, as done in [1–7]. The space–time fractional Zener wave equation represents a generalization of the classical wave equation obtained as a system consisting of the equation of motion of the deformable (one-dimensional) body, the time-fractional Zener constitutive equation and the space-fractional strain measure. Our generalization is done by the fractionalization in both space and time variables, based on the physically acceptable concepts. More details on the formulation and the mechanical background will be given in this section, which finishes with the remark related to the analysis of our generalization of the wave equation. 2014 The Author(s) Published by the Royal Society. All rights reserved.


Recall that the classical wave equation describes the waves that occur in an elastic medium. It is obtained from the equations of the deformable body [8]. The wave equation can be written in the form of a system which consists of three equations: equation of motion, constitutive equation and strain measure. Unknown functions depending on time, t > 0, and space, x ∈ R, variables are: displacement u, stress σ , measured in m, Pa, respectively, and (dimensionless) strain ε. We consider an infinite viscoelastic rod (one-dimensional body), positioned along the x-axis, that is not under the influence of body forces. Then, the equation of motion reads ∂x σ (x, t) = ρ ∂t2 u(x, t),

x ∈ R, t > 0,

(1.1)

where ρ > 0 denotes the (constant) density of the rod. The constitutive equation gives the relation between stress and strain, and in the case of elastic media it is the Hooke law. As we consider waves occurring in viscoelastic media, we choose the constitutive equation to be the time-fractional Zener model. The time-fractional Zener model was introduced in [9] in the form α α RL α σ (x, t) + τrα RL 0 Dt σ (x, t) = E(ε(x, t) + τc 0 Dt ε(x, t)),

x ∈ R, t > 0,

(1.2)

where τr and τc are relaxation and creep times (measured in s), respectively, E is the generalized α Young modulus (measured in Pa) and RL 0 Dt denotes the left Riemann–Liouville operator of fractional differentiation of order α ∈ [0, 1). It has been shown in [10] that the Riemann–Liouville fractional derivative in the above constitutive model can be replaced by the Caputo derivative yielding the same model. The fractional Zener model containing the Caputo fractional derivative was independently introduced in [11], where one also finds a discussion on the Mittag–Leffler function. We refer to [12,13] for the background of using the fractional derivatives in the theory of viscoelastic materials. For α = 0, the constitutive equation (1.2) reduces to the Hooke law σ = Er ε,

with Er = E

1 + τε . 1 + τσ

For α = 1, the constitutive equation (1.2) reduces to the classical Zener model. For more details on fractional derivatives see [14,15]. We refer to [16] for a review on the fractional models in viscoelasticity and to [17] for a systematic analysis of the thermodynamical restrictions on parameters in such models. The strain measure gives the connection between strain and displacement. The classical strain εcl (x, t) = ∂x u(x, t),

x ∈ R, t > 0,

...................................................

(a) Model

2

rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20140614

In §2, we show the existence and uniqueness of a distributional and a classical solution to the distributional (2.1) and classical space–time fractional Zener wave equation (1.17), (1.19). For this purpose, we use the Fourier and Laplace transforms in the spaces of distributions. In §2b–d, we examine quite different estimates of double integrals proving their absolute convergence. The analysis presented in §2b implies that the distributional solution of the distributional space–time fractional Zener wave equation (2.1), with u0 , v0 ∈ L1 (R), is a distribution of the second order with respect to t; see theorem 2.1. In §2c, with the assumption u0 = 0 and v0 regular enough, we obtain the classical solution to the space–time fractional Zener wave equation (1.17), (1.19); see theorem 2.2. On the other hand, in §2d, by the regularization of distribution, we show that the solution to (2.1) is given by a distributional limit of a net of approximated solutions, which are continuous with respect to x ∈ R, t > 0, bounded with respect to x ∈ R and exponentially bounded with respect to t > 0. The results in §2 are justified in §3 by discussing the influence of parameters α and β (orders of the time- and space-fractional derivatives) on the solution to (2.1) and in §4 by the numerical examples. Mathematical background is given in appendix A.


1 1 β |x|−β ∗x εcl (x, t), Ex u(x, t) = 1−β 1−β 2 Γ (1 − β)

x ∈ R, t > 0,

(1.3)

β

where (measured in m) denotes the length-scale parameter and Ex is the symmetrized Caputo fractional derivative of order β ∈ [0, 1); see appendix A. We assume that the rod is of finite crosssection √ and that the cross-sectional area A characterizes the length-scale parameter , such that = A. Note that both strain measure ε and classical strain εcl are dimensionless quantities. For the derivation of the multi-dimensional fractional strain measure using the standard continuum mechanics approach and its geometrical interpretation see [18] and references therein. The use of strain measures other than classical strain is acceptable in continuum field theory [19, p. 268]. In [20], it was shown that (1.3) can be used as a strain measure, as the displacement field is only a function of time (corresponding to rigid body motion), if and only if ε from (1.3) equals zero. The ˆ , t), Fourier transform of the strain measure (1.3) reads ε(ξ ˆ , t) = (1/ 1−β )i(ξ/|ξ |1−β ) sin(βπ/2)u(ξ ˆ , t), which is the Fourier transform of the see appendix A, which for β = 1 becomes ε(ξ ˆ , t) = iξ u(ξ classical strain εˆcl . Regarding the fractionalization of the strain measure, we follow the approach presented in [20], where the symmetrized fractional derivative is introduced in order to describe the non-local effects of the material. Note that, in [21], the same type of fractional derivative is used in the framework of the heat conduction problem of the space–time fractional Cattaneotype equation. One may also treat the non-locality in viscoelastic media by a different approach. Namely, contrary to (1.3), one may retain the classical strain measure and introduce the non-locality in the constitutive equation. In the classical setting, this was done by Eringen [22]. In the framework of the fractional calculus, this approach is followed in [23–27]. The wave equation, obtained from a system consisting of the equation of motion, the fractional Eringen-type constitutive equation and the classical strain measure, is studied in [28,29]. The approach of a construction of spatially fractional viscoelastic wave equations in three dimensions based on fractional constitutive equations in the spirit of Kunin’s stress–strain relations for crystalline solids, or Edelen’s and Eringen’s approach to non-locality, can be found in [30]. Regarding the applications of the anisotropic multi-dimensional operator considered in [30] in diffusion in biological tissues, we refer to [31,32]. We also refer to [33–36] for some recent mechanical models of fractional-order viscoelasticity. Following the approach where the classical strain measure is retained and non-locality is introduced in the constitutive equation, the system of equations (1.1)–(1.3) may be reformulated. Using the procedure described in §1b, the constitutive equation (1.2) solved with respect to σ reads α α τε τε δ(t) + − 1 eα (t) ∗t ε(x, t), x ∈ R, t > 0, (1.4) σ (x, t) = E τσ τσ where δ is the Dirac distribution, eα = (d/dt)eα , with eα being the Mittag–Leffler function eα (t) = Eα (−(t/τσ )α ); see also (1.22). When the constitutive equation, written in the form (1.4), is combined with the strain measure (1.3), a single constitutive relation between stress σ and classical strain εcl , non-local in space and time, is obtained as (x ∈ R, t > 0) α α τε τε E |x|−β δ(t) + − 1 |x|−β eα (t) ∗x,t εcl (x, t). (1.5) σ (x, t) = 1−β τσ τσ 2 Γ (1 − β) Constitutive equation (1.5) corresponds to constitutive modelling of materials that show both history-dependent properties [37] and length-scale or non-local properties [22]. Model parameters, corresponding to materials non-local in space and time described by (1.5), can be experimentally evaluated by fitting data obtained by standard measurement techniques to constitutive relation (1.5). We refer to, for example, [38] for an estimation of the model parameters appearing in history-dependent models of materials used in dentistry.

...................................................

ε(x, t) =

3


describes local deformations. As we consider non-local effects in a material, we use, at all points of the body, the fractional model of strain measure, assumed in the form of spatially averaged (weighted) classical strain as


The system of equations, equivalent to (1.1)–(1.3), consisting of the equation of motion, the space–time non-local constitutive equation and the classical strain εcl reads (x ∈ R, t > 0)

and

E 2 1−β Γ (1 − β)

τε τσ

α

|x|−β δ(t) +

τε τσ

α

(1.6)

− 1 |x|−β eα (t) ∗x,t εcl (x, t)

(1.7)

εcl (x, t) = ∂x u(x, t).

(1.8)

The initial conditions corresponding to the system of equations (1.1)–(1.3) are ∂t u(x, 0) = v0 (x),

u(x, 0) = u0 (x),

σ (x, 0) = 0

and

ε(x, 0) = 0, x ∈ R,

(1.9)

where u0 and v0 are the initial displacement and velocity, while the boundary conditions are lim u(x, t) = 0

x→±∞

and

lim σ (x, t) = 0, t > 0.

(1.10)

x→±∞

Alternatively, if the system of equations (1.6)–(1.8) is analysed, then it is subject to initial conditions (1.9)1,2 and boundary condition (1.10)1 . Note that boundary conditions (1.10) are the natural choice for the case of the unbounded domain, while in the case of the bounded domain there can be a large variety of different boundary conditions depending on the type of problem one is faced with. In the case of the local, time-fractional wave equation on a bounded domain we refer to [16,17,39] and references therein.

(b) System of equations Introducing the dimensionless quantities x x¯ = , where =

t¯ = √

t , τε

u¯ =

u ,

σ¯ =

σ , E

τ=

τσ τε

α

u¯0 =

,

u0 ,

v¯0 = v0

τε

and c˜ =

τε

ρ , E

A, in either (1.1)–(1.3) or (1.6)–(1.8), and, omitting the bar, we obtain ∂x σ (x, t) = c˜2 ∂t2 u(x, t),

x ∈ R, t > 0,

α C α σ (x, t) + τ C 0 Dt σ (x, t) = ε(x, t) + 0 Dt ε(x, t), β ε(x, t) = Ex u(x, t),

and

(1.11)

x ∈ R, t > 0,

x ∈ R, t > 0,

(1.12) (1.13)

or (x ∈ R, t > 0) ∂x σ (x, t) = c˜2 ∂t2 u(x, t), σ (x, t) = and

1 2Γ (1 − β)

1 −β |x| δ(t) + τ

1 −β − 1 |x| eα (t) ∗x,t εcl (x, t) τ

εcl (x, t) = ∂x u(x, t).

(1.14) (1.15) (1.16)

In the sequel, we take c˜ = 1. The system of equations (1.11)–(1.13), or alternatively (1.14)–(1.16), can be reduced to the space–time fractional Zener wave equation β

∂t2 u(x, t) = Lαt ∂x Ex u(x, t),

x ∈ R, t > 0,

(1.17)

where Lαt is a linear operator (of convolution type) given by Lαt = L−1

1 1 + sα 1 δ(t) + − 1 e ∗ = (t) ∗t , t α 1 + τ sα τ τ

t > 0,

(1.18)

...................................................

σ (x, t) =


∂x σ (x, t) = ρ∂t2 u(x, t),

4


∂t u(x, 0) = v0 (x),

σ (x, 0) = 0

and

ε(x, 0) = 0, x ∈ R,

(1.19)

and lim u(x, t) = 0

x→±∞

and

lim σ (x, t) = 0, t > 0.

(1.20)

x→±∞

The procedure of obtaining (1.17) is as follows. Applying the Laplace transform to (1.12) with respect to time variable t, one obtains ˜ s), (1 + τ sα )σ˜ (x, s) = (1 + sα )ε(x,

x ∈ R, Re s > 0.

[40], The inverse Laplace transform, as L−1 [(1 + sα )/(1 + τ sα )] is a well-defined element in S+ gives 1 + sα ∗t ε. (1.21) σ = L−1 1 + τ sα

Setting Lαt = L−1 [(1 + sα )/(1 + τ sα )]∗t , inserting ε, given by (1.13), into (1.21) and then inserting the obtained σ into (1.11), we obtain (1.17). Note that Lαt = L−1 [(1 + sα )/(1 + τ sα )]∗t can be explicitly expressed via the Mittag–Leffler function. Recall that for the Mittag–Leffler function eα , defined by α t , t > 0, α ∈ (0, 1), (1.22) eα (t) = Eα − τ k ∞ where Eα (z) = ∞ k=0 (z /Γ (αk + 1)), z ∈ C, we have that eα ∈ C ((0, ∞)) ∩ C([0, ∞)) and eα ∈ 1 ∞ C ((0, ∞)) ∩ Lloc ([0, ∞)), where eα (t) = (d/dt)eα (t), t > 0, and L[eα (t)](s) =

sα−1 , sα + 1/τ

Re s > 0;

see [41]. Therefore, 1 1 + sα (1 − τ )sα 1 −1 −1 (t) = δ(t) + (t) = L 1+ L − 1 eα (t), 1 + τ sα τ (sα + 1/τ ) τ τ

t > 0,

and thus we obtain Lαt as given by (1.18). For α = 0 and β = 1, i.e. when the Hooke law and the classical strain measure are used, equation (1.17) is the classical wave equation 2 . ∂t2 u = c2 ∂x2 u, with c = 1+τ Therefore, the system of equations (1.11)–(1.13), or equivalently (1.17), generalize the classical wave equation. We collect other special cases of (1.17) in the following remark. Remark 1.1. Generalizations of the classical wave equation, given by the system of equations (1.11)–(1.13), or (1.17), are distinguished and classified according to parameter β as follows. (i) Case β = 0. We obtain the non-propagating disturbance if v0 = 0. Namely, for β = 0, we obtain ε = 0, owing to (1.13) and the property of the symmetrized fractional derivative that Ex0 u = 0 (see appendix A). This and (1.12) imply σ = 0, so that from (1.11), (1.19) and (1.20) one obtains u(x, t) = u0 (x) + v0 (x)t, x ∈ R, t ≥ 0. (1.23) Note, for v0 = 0, we have u(x, t) = u0 (x), x ∈ R, t ≥ 0. (ii) Case β ∈ (0, 1). For α = 0, we obtain the space-fractional wave equation 2 β , x ∈ R, t > 0, ∂t2 u(x, t) = c2 ∂x Ex u(x, t) and c = 1+τ

(1.24)

...................................................

u(x, 0) = u0 (x),

5


and L−1 denotes the inverse Laplace transform (see appendix A). The dimensionless quantities give that initial and boundary conditions, (1.10) and (1.9), for the space–time fractional Zener wave equation (1.17) again become


x ∈ R, t > 0.

(1.25)

For all α ∈ [0, 1], when β tends to zero, the solution to the system of equations (1.11)–(1.13), (1.19) and (1.20) tends to (1.23); see §3. This suggests that the parameter β measures the resistance of the material to the propagation of initial disturbance. (iii) In the case when β = 1, α ∈ (0, 1), equation (1.17) reduces to the time-fractional Zener wave equation (1.26) ∂t2 u(x, t) = Lαt ∂x2 u(x, t), x ∈ R, t > 0, studied in [42,43]. For α = 0, as already mentioned above, we obtain the classical wave equation and for α = 1 the Zener wave equation ∂t2 u(x, t) = L1t ∂x2 u(x, t),

x ∈ R, t > 0.

2. Cauchy problems (2.1) and (1.17), (1.19) (a) Framework The framework for our analysis is the spaces of the distributions: S (R) (or shortly S ) and K (R) (or K ), the duals of the Schwartz space S(R) (or S) and of the space K(R) (or K); K is the space of smooth functions ϕ with the property that there exists m ∈ N0 , such that supx∈R,α≤m |ϕ (α) (x)| em|x| < ∞. The elements of S , respectively of K , are of the form f = r (α) k0 α=0 Φα , where Φα are continuous functions on R and |Φα (t)| ≤ C(1 + |t|) , respectively k |t| 0 |Φα (t)| ≤ C e , α ≤ r, t ∈ R, for some C > 0, r ∈ N0 and k0 ∈ N0 . The space S+ (K+ ) is a subspace , respectively of K , are of S (K ) consisting of elements supported by [0, ∞). The elements of S+ + k (p) kt (p) of the form f (t) = (Φ(t)(1 + |t|) ) , respectively f (t) = (Φ(t) e ) , t ∈ R, where Φ is a continuous are subspaces of K and K , bounded function such that Φ(t) = 0, t ≤ 0. Note that S and S+ + respectively. The elements of K+ have the Laplace transform, which are analytic functions in the domain Re s > s0 > 0. We also recall that for the Lebesgue spaces of integrable and bounded functions L1 (R) and L∞ (R), f ∗ g ∈ L∞ (R) if f ∈ L1 (R) and g ∈ L∞ (R). We shall apply the Fourier transform with respect to x and the Laplace transform with respect , which is the subspace to t. Actually, we shall consider the distributions within the space S ⊗ K+ 2 of K (R ), consisting of distributions having support in R × [0, ∞). For the background of the tensor product, we refer to [44]. We shall obtain the solution u as an element of C(R) ∩ L∞ (R) for for fixed ψ ∈ K, i.e. u(x, t), ψ(x) ∈ fixed ϕ ∈ K, i.e. u(x, t), ϕ(t) ∈ C(R) ∩ L∞ (R) and elements of K+ K+ .

(b) Existence and uniqueness of a generalized solution We consider the existence and uniqueness of the solution to the Cauchy problem (1.17), (1.19) and (1.20). If u0 ∈ C1 (R) and v0 ∈ C(R), then the classical solution to the Cauchy problem (1.17), (1.19) and (1.20) is a function u(x, t) of class C2 for t > 0, of class C1 for t ≥ 0, which satisfies equation (1.17) for t > 0 and initial conditions (1.19) when t = 0, as well as the boundary conditions (1.20). If the function ucl is continued by zero for t < 0, then putting u(x, t) = ucl (x, t)H(t),

x, t ∈ R,

we obtain the distributional space–time fractional Zener wave equation β

∂t2 u(x, t) = Lαt ∂x Ex u(x, t) + u0 (x)δ (t) + v0 (x)δ(t),

in K (R2 ).

We shall refer to solution u of (2.1) as a generalized, or distributional, solution.

(2.1)

...................................................

β

∂t2 u(x, t) = L1t ∂x Ex u(x, t),

6


studied in [20]. Case α ∈ (0, 1), to the best of the authors’ knowledge, has not been studied in the literature, so it is the subject of the analysis presented in this work. For α = 1, (1.17) becomes the space-fractional Zener wave equation


Proof. Formally applying the Laplace transform to (2.1) with respect to t, we obtain β

˜ s) − s2 ∂x Ex u(x,

1 + τ sα 1 + τ sα ˜ (su0 (x) + v0 (x)), u(x, s) = − 1 + sα 1 + sα

x ∈ R, Re s > s0 ,

(2.3)

for arbitrary, but fixed, s0 > 0, where u˜ is an analytic function with respect to s, Re s > s0 . Equation (2.3) is of the type β

∂x Ex u(x) − ωu(x) = −νu0 (x) − μv0 (x),

x ∈ R,

(2.4)

where ω = ω(s) = s2

1 + τ sα , 1 + sα

ν = ν(s) = s

1 + τ sα 1 + sα

and μ = μ(s) =

1 + τ sα , 1 + sα

Re s > s0 .

(2.5)

We have shown in [42, theorem 4.2] that ω(s) ∈ C \ (−∞, 0] for Re s > 0. For fixed s, Re s > s0 , the unique solution u ∈ C(R) ∩ L∞ (R) to (2.4), given by ∞ 1 1 cos(ρx) dρ, x ∈ R, (2.6) u(x) = (νu0 (x) + μv0 (x)) ∗x 1+β π sin(βπ/2) + ω 0 ρ is obtained as the inverse Fourier transform of ˆ )= u(ξ

ν uˆ0 (ξ ) + μvˆ0 (ξ ) , |ξ |1+β sin(βπ/2) + ω

ξ ∈ R.

The previous expression, with ω, ν and μ given by (2.5), takes the form ˆ˜ , s) = u(ξ

s2

suˆ0 (ξ ) + vˆ0 (ξ ) α + ((1 + s )/(1 + τ sα ))|ξ |1+β

sin(βπ/2)

,

ξ ∈ R, Re s > s0 .

(2.7)

The Fourier and Laplace transform of u (2.7), with u0 = v0 = 0, shows the uniqueness ˆ s) ∈ L∞ (R) ∩ C(R), as x → of the distributional solution to (2.1). For fixed s, Re s > s0 , u(·, ∞ 1+β sin(βπ/2) + ω)) cos(ρx) dρ is a continuous bounded function and u0 , v0 ∈ L1 (R). 0 (1/(ρ Therefore, by (2.6), we have the solution to (2.3) in the form ˜ s) = u(x,

=

1 1 + τ sα (su0 (x) + v0 (x)) π 1 + sα ∞ 1 cos(ρx) dρ ∗x 1+β sin(βπ/2) + s2 ((1 + τ sα )/(1 + sα )) ρ 0 1 (su0 (x) + v0 (x)) π ∞ 1 cos(ρx) dρ, ∗x 2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2) s 0

x ∈ R, Re s > s0 .

(2.8)

...................................................

Functions I, J1+ , J1− , J2+ and J2− , given by (2.12), (2.15), (2.19), (2.16) and (2.20), respectively, are bounded and continuous on R, for fixed t ≥ 0, and continuous, exponentially bounded on [0, ∞), for fixed x ∈ R.

7


Theorem 2.1. Let α ∈ [0, 1), β ∈ [0, 1), τ ∈ (0, 1) and let u0 , v0 ∈ L1 (R). Then there exists a unique generalized solution u ∈ K (R2 ), supp u ⊂ R × [0, ∞), to the distributional space–time fractional Zener wave equation (2.1). More precisely, u is of the form ⎫ 1 ⎪ ⎪ (δ (t)u (x) + δ(t)v (x)) ∗ P(x, t), x ∈ R, t > 0, u(x, t) = ⎪ x,t 0 0 ⎪ ⎪ 2π 2 ⎪ ⎪ ⎬ 2 ∂ ∂ s0 t (2.2) P(x, t) = I(x, t) − J1 (x, t) + 2 J2 (x, t) e , x ∈ R, t > 0,⎪ ⎪ ∂t ∂t ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ and J1 = i(J1+ − J1− ), J2 = J2+ + J2− .


Thus, the justification of the previously presented procedure is based on the analysis of the inverse Laplace transform. Formally, when applied to (2.8), the inverse Laplace transform gives

∞

P(x, t) = 2π L−1 = −i

0

cos(ρx) dρ (x, t) s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)

s0 +i∞ ∞ s0 −i∞

0

x ∈ R, t > 0,

cos(ρx) est s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β

sin(βπ/2)

dρ ds,

(2.9)

(2.10) x ∈ R, t > 0,

(2.11)

is a fundamental solution to (2.1). Consider the divergent integral (2.11). We introduce the parametrization s = s0 + ip, p ∈ R, in (2.11), so that P(x, t) =

∞ ∞ −∞ 0

cos(ρx) es0 t eipt dρ dp [s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s0 +ip

= I(x, t) + I+ (x, t) + I− (x, t),

x ∈ R, t > 0,

where I, I+ and I− are given below. The integral I(x, t) =

p0 ∞ −p0 0

[s2

cos(ρx) es0 t eipt dρ dp, sin(βπ/2)]s=s0 +ip

+ ((1 + sα )/(1 + τ sα ))ρ 1+β

is absolutely convergent, as p0 ∞ |I(x, t)| ≤ es0 t −p0 0

x ∈ R, t > 0.

In (2.13), p0 is chosen so that s2 +

(2.12)

1 Re([s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s0 +ip )

× dρ dp < ∞,

Re

x ∈ R, t > 0,

1 + sα 1+β βπ ρ sin 1 + τ sα 2

= r2 cos(2ϕ) +

(2.13)

s=s0 +ip

1 + (1 + τ )rα 1 + 2τ rα

cos(αϕ) + τ r2α 1+β βπ ρ sin > 0, 2 cos(αϕ) + τ 2 r2α

with r = s20 + p2 and tan ϕ = p/s0 . Note that for p = 0 and ρ = 0 the integrand is well defined, due to s0 > 0. Thus, the integral I exists and belongs to C(R) ∩ L∞ (R) with respect to x, and it is a continuous exponentially bounded function with respect to t ≥ 0. Next, we consider ∞ ∞

cos(ρx) eipt dρ dp 2 α α 1+β sin(βπ/2)] p0 0 [s + ((1 + s )/(1 + τ s ))ρ s=s0 +ip ∂ + ∂2 + = − i J1 (x, t) + 2 J2 (x, t) es0 t , x ∈ R, t > 0, ∂t ∂t

I+ (x, t) = es0 t

(2.14)

where J1+ (x, t) = and J2+ (x, t) =

∞ 1 p0

0

cos(ρx) eipt p[s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β

∞ ∞ p0

1

sin(βπ/2)]s=s0 +ip

dρ dp,

cos(ρx) eipt p2 [s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β


dρ dp,

x ∈ R, t > 0,

(2.15)

x ∈ R, t > 0. (2.16)

...................................................

where

1 (δ (t)u0 (x) + δ(t)v0 (x)) ∗x,t P(x, t), 2π 2


u(x, t) =

8


Function J1+ , given by (2.15), is continuous on R × [0, ∞), bounded with respect to x and exponentially bounded with respect to t, as

0

where we used s2 +

Im


s=s0 +ip

rα sin(αϕ) βπ = r2 sin(2ϕ) + (1 − τ ) ρ 1+β sin >0 α 2 2α 2 1 + 2τ r cos(αϕ) + τ r 1 − τ 1 1+β βπ , as r → ∞, ρ sin(αϕ) sin 2 τ 2 rα βπ 1 − τ 1 1+β απ sin , as p → ∞. ∼ 2s0 p + ρ sin α 2 2 2 τ p

∼ r2 sin(2ϕ) +

(2.17)

In obtaining (2.17), we used: r2 sin(2ϕ) = r2 (2 tan ϕ/(1 + tan2 ϕ)) = 2s0 p, as well as ϕ ∼ π/2 and r ∼ p, as p → ∞. Function J2+ , given by (2.16), is continuous on R × [0, ∞), bounded with respect to x and exponentially bounded with respect to t, as, by (2.17) and by the Fubini theorem, ∞ ∞

1 dρ dp p2 Im([s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s0 +ip ) ∞ ∞ 1 ≤ dρ dp 3 2 2−α ρ 1+β sin(απ/2) sin(βπ/2) p0 1 2s0 p + ((1 − τ )/τ )p ∞ ∞ τ2 1 1 1 ≤ dρ dp < ∞, x ∈ R, t > 0. 1+β 1 − τ sin(απ/2) sin(βπ/2) p0 p2−α 1 ρ

|J2+ (x, t)| ≤

p0

1

Thus, we have that I+ , given by (2.14), belongs to C(R) ∩ L∞ (R) with respect to x, and it is a derivative of a continuous exponentially bounded function with respect to t ≥ 0. Similarly as I+ , we analyse I− (x, t) = es0 t

−p0 ∞ −∞ 0

∞ ∞

cos(ρx) eipt [s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β


dρ dp

cos(ρx) e−iqt dρ dq 2 α α 1+β sin(βπ/2)] p0 0 [s + ((1 + s )/(1 + τ s ))ρ s=s0 −iq ∂ − ∂2 − = i J1 (x, t) − 2 J2 (x, t) es0 t , x ∈ R, t > 0, ∂t ∂t = es0 t

(2.18)

where J1− (x, t) =

∞ 1 p0

0

cos(ρx) e−iqt q[s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β

q2 [s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β

sin(βπ/2)]s=s0 −iq

dρ dq

x ∈ R, t > 0,

(2.19)

and J2− (x, t) =

∞ ∞ p0

1

cos(ρx) e−iqt sin(βπ/2)]s=s0 −iq

dρ dq,

x ∈ R, t > 0. (2.20)

...................................................

p0


∞ 1

1 dρ dp p Im([s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s0 +ip ) 1 ∞ 1 1 ≤ dρ dp < ∞, x ∈ R, t > 0, 2s0 p0 0 p2

|J1+ (x, t)| ≤

9


By (2.19), we have ∞ 1 0

1 dρ dq, q|Im([s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s0 −iq )|

x ∈ R, t > 0. Function J1− is continuous on R × [0, ∞), bounded with respect to x and exponentially bounded with respect to t, as 1 ∞ 1 1 dρ dq < ∞, x ∈ R, t > 0, |J1− (x, t)| ≤ 2s0 p0 0 q2 where r = s20 + q2 , tan ϕ = −q/s0 and where we used Im

1 + sα 1+β βπ s + ρ sin 1 + τ sα 2

2

= r2 sin(2ϕ) + (1 − τ )

s=s0 −iq

rα sin(αϕ) βπ ρ 1+β sin 2 cos(αϕ) + τ 2 r2α

1 + 2τ rα

1 − τ 1 1+β βπ < 0, as r → ∞, ρ sin(αϕ) sin ∼ r2 sin(2ϕ) + α 2 2 τ r βπ 1 − τ 1 1+β απ sin , as q → ∞. ∼ − 2s0 q + ρ sin 2 2 τ 2 qα Consider (2.20): |J2− (x, t)| ≤

∞ ∞

1

dρ dq sin(βπ/2)]s=s0 −iq )| ∞ ∞ 1 ≤ dρ dq 3 2 2−α ρ 1+β sin(απ/2) sin(βπ/2) p0 1 2s0 q + ((1 − τ )/τ )q ∞ ∞ τ2 1 1 1 ≤ dρ dq, x ∈ R, t > 0. 1+β 1 − τ sin(απ/2) sin(βπ/2) p0 q2−α 1 ρ p0

1

q2 |Im([s2

+ ((1 + sα )/(1 + τ sα ))ρ 1+β

Function J2− is continuous on R × [0, ∞), bounded with respect to x and exponentially bounded with respect to t. Thus, we have that I− , given by (2.18), belongs to C(R) ∩ L∞ (R) with respect to x, and it is a derivative of a continuous exponentially bounded function with respect to t ≥ 0. Thus, P has the form ∂2 + ∂ + − − P(x, t) = I(x, t) − i (J1 (x, t) − J1 (x, t)) + 2 (J2 (x, t) + J2 (x, t)) es0 t , x ∈ R, t > 0, ∂t ∂t where I, J1+ , J1− , J2+ and J2− are bounded and continuous functions with respect to x and continuous exponentially bounded functions with respect to t.

(c) Existence and uniqueness of a classical solution If the initial condition u0 is assumed to be zero, one may also obtain the existence and uniqueness result for the classical solution. Theorem 2.2. Assume u0 = 0 and v0 ∈ C5 (R) is compactly supported. Then the space–time fractional Zener wave equation (1.17), subject to initial conditions (1.19)1,2 , admits the (classical) solution of the C1 class with respect to x ∈ R and of the C2 class with respect to t ∈ [0, ∞), given in the form ⎞ ⎞ ⎛⎛ 2 4 1 ∂ (2.21) ucl (x, t) = − 2 (tH(t)) ∗t ⎝⎝ 1 − 2 v0 (x)⎠ ∗x I(x, t)⎠ + v0 (x)t, x ∈ R, t > 0. 2π ∂x

...................................................

p0


|J1− (x, t)| ≤

10


u(x, t) =

1 (δ(t)v0 (x)) ∗x,t P(x, t), 2π 2

x ∈ R, t > 0.

(2.22)

Further assume that P, given by (2.10), or alternatively by (2.11), with P(x, t) = 0, x ∈ R, t > 0, is considered in the sense of distributions. Let Q(x, t) = −i

s0 +i∞ ∞ s0 −i∞

cos(ρx) est dρ ds, (1 + ρ 2 )4 (s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2))

0

so that

∂2 P(x, t) = 1 − 2 ∂x

4 Q(x, t).

(2.23)

In the sense of distributions, we have ∂2 Q(x, t) ∂t2 s0 +i∞ ∞ = −i 1−1+ s0 −i∞

= −i

0

s0 +i∞ ∞ s0 −i∞

0

s2 s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin

βπ 2

cos(ρx) est dρ ds (1 + ρ 2 )4

cos(ρx) est cos(ρx) est ρ 1+β sin(βπ/2) − 2 2 4 α α 1+β (1 + ρ ) s ((1 + τ s )/(1 + s )) + ρ sin(βπ/2) (1 + ρ 2 )4

= 2π g(x)δ(t) − I(x, t),

dρ ds

x, t ∈ R,

(2.24)

where I(x, t) = −i

s0 +i∞ ∞ s0 −i∞

0

cos(ρx) est ρ 1+β sin(βπ/2) dρ ds, s2 ((1 + τ sα )/(1 + sα )) + ρ 1+β sin(βπ/2) (1 + ρ 2 )4

x ∈ R, t > 0, (2.25)

∞ and g(x) = 0 (cos(ρx)/(1 + ρ 2 )4 ) dρ, x ∈ R, is a continuous function. We used that 2π δ(t) = s +i∞ −i s00−i∞ est ds, t > 0, in the sense of distributions. Our formal calculations will be justified by showing that the double integral (2.25) exists as the Lebesgue integral and that it is a continuous function for x ∈ R, t > 0. We introduce the parametrization s = s0 + ip, p ∈ R, in (2.25), so that I(x, t) =

∞ ∞ −∞ 0

ρ 1+β sin(βπ/2) cos(ρx) es0 t eipt dρ dp [s2 ((1 + τ sα )/(1 + sα )) + ρ 1+β sin(βπ/2)]s=s0 +ip (1 + ρ 2 )4

= I1 (x, t) + I2 (x, t) + I3 (x, t),

x ∈ R, t > 0,

where I1 , I2 and I3 will be defined below. We put I1 (x, t) =

∞ ∞ p0

×

0

ρ 1+β sin(βπ/2) sin(βπ/2)]s=s0 +ip

[s2 ((1 + τ sα )/(1 + sα )) + ρ 1+β

cos(ρx) es0 t eipt dρ dp, (1 + ρ 2 )4

x ∈ R, t > 0.

...................................................

Proof. Assume that u0 = 0 in (2.2), so that

11


Note that the assumptions on v0 imply that ((1 − ∂ 2 /∂x2 )4 v0 (x)) ∗x I(x, t) is a C1 function with respect to x. Moreover, we can assume less regular conditions on v0 , which would lead to the use of the pseudodifferential operators.


(1 + (1 + τ )rα cos(αϕ) + τ r2α ) sin(2ϕ) + (τ − 1)rα sin(αϕ) cos(2ϕ) . 1 + 2rα cos(αϕ) + r2α

Since ϕ ∈ (π/4, π/2), we have cos(2ϕ) < 0. Moreover, by assumption, τ < 1. Thus, we have α Im s2 1 + τ s + ρ 1+β sin βπ ≥ cr2+α | cos(αϕ) sin(2ϕ) − sin(αϕ) cos(2ϕ)| 1 + sα 2 r2α √ 2 . ≥ Cr2−α , r ≥ s0 2 Therefore, |I1 (x, t)| ≤ C es0 t

∞ ∞ s0

0

1 s20 + p2

2−α

ρ 1+β dρ dp, (1 + ρ 2 )4

(2.26)

x ∈ R, t > 0,

so that I1 is a continuous, bounded function with respect to x and a continuous, exponentially bounded function with respect to t. We put −s0 ∞ ρ 1+β sin(βπ/2) cos(ρx) es0 t eipt dρ dp I2 (x, t) = 2 α α 1+β sin(βπ/2)] (1 + ρ 2 )4 −∞ 0 [s ((1 + τ s )/(1 + s )) + ρ s=s0 +ip ∞ ∞ ρ 1+β sin(βπ/2) = 2 α α 1+β sin(βπ/2)] s0 0 [s ((1 + τ s )/(1 + s )) + ρ s=s0 −ip ×

cos(ρx) es0 t e−ipt dρ dp, (1 + ρ 2 )4

x ∈ R, t > 0.

A similar estimate holds for (2.26). Thus, ∞ ∞ 1 ρ 1+β dρ dp, |I2 (x, t)| ≤ C es0 t 2−α (1 + ρ 2 )4 s0 0 s20 + p2

x ∈ R, t > 0,

so that I2 is a continuous, bounded function with respect to x and a continuous, exponentially bounded function with respect to t. We put s0 ∞ ρ 1+β sin(βπ/2) cos(ρx) es0 t eipt dρ dp, I3 (x, t) = 2 α α 1+β sin(βπ/2)]s=s0 +ip (1 + ρ 2 )4 −s0 0 [s ((1 + τ s )/(1 + s )) + ρ x ∈ R, t > 0, and write 1 + τ sα βπ 1+β + ρ sin Re s2 1 + sα 2

1 + (1 + τ )rα cos(αϕ) + τ r2α (τ − 1)rα sin(αϕ) = Re (r cos(2ϕ) + ir sin(2ϕ)) + i 1 + 2rα cos(αϕ) + r2α 1 + 2rα cos(αϕ) + r2α 2

= r2

2

(1 + (1 + τ )rα cos(αϕ) + τ r2α ) cos(2ϕ) + (1 − τ )rα sin(αϕ) cos(2ϕ) , 1 + 2rα cos(αϕ) + r2α

...................................................

= r2

12


We put s = r eiϕ , where r = s20 + p2 , tan ϕ = p/s0 and also put p0 = s0 , so that ϕ ∈ (π/4, π/2). This implies 1 + τ sα βπ 1+β + ρ sin Im s2 1 + sα 2 1 + (1 + τ )rα cos(αϕ) + τ r2α (τ − 1)rα sin(αϕ) 2 2 = Im (r cos(2ϕ) + ir sin(2ϕ)) +i 1 + 2rα cos(αϕ) + r2α 1 + 2rα cos(αϕ) + r2α


so that |I3 (x, t)| ≤ C es0 t

s0 ∞

x ∈ R, t > 0.

∂2 Q(x, t) = 2π g(x)δ(t) − I(x, t), ∂t2

x ∈ R, t > 0,

(2.27)

where g is a continuous function and I is a continuous, bounded function with respect to x and a continuous, exponentially bounded function with respect to t. By (2.22) and (2.23), we have ⎞ ⎛ 2 4 ∂ 1 (δ(t)v0 (x)) ∗x,t ⎝ 1 − 2 Q(x, t)⎠ u(x, t) = 2π 2 ∂x ⎞ ⎛ 4 1 ⎝ ∂2 = δ(t) 1 − 2 v0 (x)⎠ ∗x,t Q(x, t), 2π 2 ∂x

x ∈ R, t > 0,

so that, by (2.27),

⎞ ⎛ 4 2 ∂2 1 ⎝ ∂2 ⎠ ∗x,t ∂ Q(x, t) u(x, t) = v (x) δ(t) 1 − 0 ∂t2 2π 2 ∂x2 ∂t2 ⎞ ⎛ 4 1 ⎝ ∂2 = δ(t)v0 (x) − 1 − 2 v0 (x)⎠ ∗x I(x, t), 2π 2 ∂x

as (1 − ∂ 2 /∂x2 )4 g(x) = π δ(x). Using u(x, t) = ucl (x, t)H(t),

∂2 i.e. 2 u(x, t) = ∂t

x ∈ R, t > 0,

∂2 ucl (x, t) H(t) + δ(t)v0 (x) ∂t2

and the previous equation, we have

⎞ ⎛ 4 ∂2 1 ⎝ ∂2 ucl (x, t) = − 2 1 − 2 v0 (x)⎠ ∗x I(x, t), ∂t2 2π ∂x

x ∈ R, t > 0,

yielding (2.21).

(d) Regularization of a generalized solution We give a regularization of the generalized solution u to the space–time fractional Zener wave equation (2.1), which is of particular importance for the numerical analysis of the problem. We start from the Fourier and Laplace transform of the solution given by (2.7) and write it as ˆ˜ , s), ξ ∈ R, Re s > s , ˆ˜ , s) = uˆ0 (ξ ) + 1 vˆ0 (ξ ) K(ξ (2.28) u(ξ 0 s where ˆ˜ , s) = K(ξ ˆ˜ , s) = Q(ξ Note that

s

=

s2

+ ((1 + sα /(1 + τ sα ))|ξ |1+β

s3

((1 + sα )/(1 + τ sα ))|ξ |1+β sin(βπ/2) , + s((1 + sα )/(1 + τ sα ))|ξ |1+β sin(βπ/2)

ˆ˜ , s), ˆ˜ , s) = 1 sP(ξ K(ξ 2π

i.e. K(x, t) =

sin(βπ/2)

1 ∂ P(x, t), 2π ∂t

1 ˆ˜ , s), − Q(ξ s

with

(2.29)

ξ ∈ R, Re s > s0 .

(2.30)

x, ξ ∈ R, t > 0,

...................................................

The integral I3 is a continuous, bounded function with respect to x and a continuous, exponentially bounded function with respect to t. By (2.24), in the sense of distributions, we have


−s0 0

13

ρ 1+β sin(βπ/2) dρ dp, (1 + ρ 2 )4


where P is given by (2.10). We already know from theorem 2.1 that P (and also therefore K) is a distribution. We regularize Kˆ˜ by multiplying it by the Fourier transform of the Gaussian

which is a δ-net, i.e. the Gaussian in a limiting process ε → 0 represents the Dirac delta distribution. Thus, we have that ˆ˜ , s) e−(εξ )2 /4 , Kˆ˜ ε (ξ , s) = K(ξ

where F [δε (x)](ξ ) = e−(εξ )

2

/4

, ξ ∈ R, Re s > s0 , ε ∈ (0, 1],

(2.31)

has the inverse Laplace and Fourier transforms which is a function and, in a distributional limit, gives the solution kernel K as a distribution. We summarize these observations in the following theorem, which is given after we state the lemma. Lemma 2.3. Let α ∈ [0, 1), τ ∈ (0, 1) and θ > 0. Then Ψα (s) = s2 + θ

1 + sα , 1 + τ sα

s ∈ C,

admits exactly two zeros. They are complex-conjugate, located in the left complex half-plane and each of them is of multiplicity one. Theorem 2.4. Let all conditions of theorem 2.1 be satisfied. Let u ∈ K (R2 ), with support in R × [0, ∞), be the distributional solution to the distributional space–time fractional Zener wave equation (2.1). Then u is of the form u(x, t) = (u0 (x)δ(t) + v0 (x)H(t)) ∗x,t K(x, t), where K is a distributional limit in K (R2 ), K(x, t) = lim Kε (x, t) and ε→0

with

Kε (x, t) =

1 π

∞

S(ρ, t) cos(ρx) e−(ερ)

2

/4

(2.32)

dρ,

0

x ∈ R, t > 0,

∞

1 2 + ((1 + qα eiαπ )/(1 + τ qα eiαπ ))ρ 1+β sin(βπ/2) q 0 1 qe−qt dq − 2 q + ((1 + qα e−iαπ )/(1 + τ qα e−iαπ ))ρ 1+β sin(βπ/2) s est + α−1 α 2 1+β 2s + (α(1 − τ )s /(1 + τ s ) )ρ sin(βπ/2) s=sz (ρ) s est + , 2s + (α(1 − τ )sα−1 /(1 + τ sα )2 )ρ 1+β sin(βπ/2) s=s¯z (ρ)

1 S(ρ, t) = 2π i

(2.33)

(2.34)

and sz are zeros of Ψα from lemma 2.3. In particular, for suitable s0 > 0, Kε (x, t) e−s0 t is bounded and continuous with respect to x ∈ R, t > 0, for every ε ∈ (0, 1]. Proof of lemma 2.3. Let s = reiϕ , r > 0, ϕ ∈ (−π , π ). We have Re Ψα (s) = r2 cos(2ϕ) + θ

1 + (1 + τ )rα cos(αϕ) + τ r2α 1 + 2τ rα cos(αϕ) + τ 2 r2α

and Im Ψα (s) = r2 sin(2ϕ) + θ (1 − τ )

rα sin(αϕ) . cos(αϕ) + τ 2 r2α

1 + 2τ rα

From (2.35) and (2.36), one can easily see that Ψα (sz ) = 0 implies Ψα (s¯z ) = 0.

(2.35)

(2.36)

...................................................

x ∈ R, ε ∈ (0, 1],


1 2 2 δε (x) = √ e−x /ε , ε π

14


C3 : s = x eiπ ; x ∈ [r, R],

π ,π , C2 : s = R eiϕ ; ϕ ∈ 2 π C4 : s = r eiϕ ; ϕ ∈ ,π , 2

where r < r0 , R > R0 and r0 , R0 are chosen as follows: r0 is small enough such that for all r < r0 it holds that Re Ψ (s) ∼ θ and Im Ψ (s) ∼ θ(1 − τ )rα sin(απ ) and, therefore, there are no zeros for r < r0 , and R is large enough such that for all R > R0 it holds that Re Ψ (s) ∼ R2 cos(2ϕ) and Im Ψ (s) ∼ R2 sin(2ϕ). On the contour C1 , we have that Im Ψα (s) ≥ 0 (as τ ∈ (0, 1) and α ∈ [0, 1) implies sin(απ/2) ≥ 0), and Im Ψα (s) → 0 for r, x → 0, as well as for R, x → ∞. The real part of Ψα varies from θ (for r, x → 0) to −∞ (for R, x → ∞). Therefore, on C1 we have Ψα (s) = −π . On the contour C2 , for R > R0 , R0 large enough, we have Im Ψ ∼ R2 sin(2ϕ) +

θ (1 − τ ) sin(αϕ) 1 ∼ R2 sin(2ϕ) ≤ 0, Rα τ2

and we have Im Ψα (s) → 0 for both ϕ = π/2 and ϕ = π . The real part of Ψα changes from −∞ (for ϕ = π/2) to ∞ (for ϕ = π ) because Re Ψ (s) ∼ R2 cos(2ϕ) +

θ ∼ R2 cos(2ϕ). τ

So, on C2 the change of the argument is Ψα (s) = −π . On the contours C3 and C4 the argument does not change. On C3 the imaginary part of Ψα is always positive and it tends to zero for both R → ∞ and r → 0, while the real part changes from ∞ (for R → ∞) to θ (for r → 0), and even if it changes the sign it does not change the argument of Ψα . On C4 it holds that Ψα (s) ∼ θ and so there is no change of argument. Taking everything together we have Ψα (s) = −2π as s ∈ C, and by the argument principle there is one zero inside of the contour C. Therefore, there is a unique pair of complex-conjugate numbers in the left complex plane which are zeros of Ψα . Proof of theorem 2.4. Let Kˆ˜ ε (ξ , s) =

s 1 2 2 e−(εξ ) /4 = e−(εξ ) /4 s s2 + ((1 + sα )/(1 + τ sα ))|ξ |1+β sin(βπ/2) ˆ˜ (ξ , s), −Q ε

ξ ∈ R, Re s > s0 ,

by (2.29) and (2.31). We shall prove that ˆ˜ , s) e−(εξ )2 /4 = ˆ˜ (ξ , s) = Q(ξ Q ε ξ ∈ R, Re s > s0 ,

s3

((1 + sα )/(1 + τ sα ))|ξ |1+β sin(βπ/2) 2 e−(εξ ) /4 , + s((1 + sα )/(1 + τ sα ))|ξ |1+β sin(βπ/2)

(2.37)

...................................................

C1 : s = x eiπ/2 ; x ∈ [r, R],

15


Next we show that if Re sz > 0, then such sz cannot be a zero of Ψα and therefore zeros must lie in the left complex half-plane. Suppose Re s > 0, i.e. ϕ ∈ (−π/2, π/2). Since zeros appear in complex-conjugate pairs, we can suppose that ϕ ∈ [0, π/2). For ϕ = 0, we have Ψα (s) > 0. Since α ∈ [0, 1), we have that αϕ, 2ϕ ∈ (0, π ) and therefore sin(2ϕ) > 0 and sin(αϕ) > 0, which together with θ > 0 and τ ∈ (0, 1) implies Im Ψα (s) > 0, so such an s cannot be a zero. It is left to show that there is only one pair of zeros of Ψα . We use the argument principle. Recall that if f is an analytic function inside and on a regular closed curve C, and non-zero on C, then the number of zeros of f (counted as many times as its multiplicity) inside the contour C is equal to the total change in the argument of f (s) as s travels around C. For our purpose, we choose contour C = C1 ∪ C2 ∪ C3 ∪ C4 , parametrized as


1 (Jε (x, t) + Jε+ (x, t) + Jε− (x, t)) es0 t , 2π 2

with (x ∈ R, t > 0, ε ∈ (0, 1]) Jε (x, t) =

p0 ∞ −p0 0

(2.38)

((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2) s3 + s((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)

s=s0 +ip

× cos(ρx) e−(ερ) /4 eipt dρ dp, ∞ ∞ ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2) + Jε (x, t) = 3 α α 1+β sin(βπ/2) p0 0 s + s((1 + s )/(1 + τ s ))ρ 2

(2.39)

s=s0 +ip

−(ερ)2 /4 ipt

Jε− (x, t) =

and

× cos(ρx) e p0 ∞ −∞ 0

e

dρ dp

((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2) s3 + s((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)

(2.40)

s=s0 +ip

−(ερ)2 /4 ipt

× cos(ρx)e

e dρ dp,

(2.41)

where we introduced the parametrization s = s0 + ip, p ∈ (−∞, ∞) in (2.38) and used the fact that ˆ˜ is an even function in ξ . From (2.39), we have (x ∈ R, t > 0, ε ∈ (0, 1]) Q ε |Jε (x, t)| ≤

p0 ∞ −p0 0

|(1 + sα )/(1 + τ sα )|s=s0 +ip |ρ 1+β [s3 + s((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s

e−(ερ) +ip

2

0

/4

dρ dp < ∞.

Let us estimate the integral given by (2.40) as (x ∈ R, t > 0, ε ∈ (0, 1]) |Jε+ (x, t)| ≤ ≤

∞ ∞ p0

e−(ερ)

|(1 + sα )/(1 + τ sα )|s=s0 +ip |ρ 1+β Im([s3 + s((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s

e−(ερ) +ip )

2

0 +ip

0

∞ ∞ p0

|(1 + sα )/(1 + τ sα )|s=s0 +ip |ρ 1+β [s3 + s((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)]s=s

0

0

Re

1 + sα 1 + τ sα

=

1 + (1 + τ )rα cos(αϕ) + τ r2α 1 ∼ , τ 1 + 2τ rα cos(αϕ) + τ 2 r2α

rα sin(αϕ) 1−τ 1 1 + sα = (1 − τ ) ∼ sin(αϕ), 1 + τ sα 1 + 2τ rα cos(αϕ) + τ 2 r2α τ 2 rα and thus for r = s20 + p2 , tan ϕ = s0 /p, Im

1 + sα 1+β βπ s +s ρ sin 1 + τ sα 2

/4

dρ dp.

(2.42)

as p → ∞,

as r → ∞,

Im

dρ dp 2

We have Im(s0 + ip)3 = −p3 + 3s20 p, |(1 + sα )/(1 + τ sα )|s=s0 +ip | ∼ 1/τ and 1 + sα 1+β 1−τ 1 βπ απ 1 , ρ sin Im s ∼ p + s0 2 α sin 1 + τ sα 2 s=s0 +ip τ 2 τ p as

/4

as r → ∞,

3

∼ −p3 + 3s20 p +

s=s0 +ip

1−τ 1 1 απ βπ p + s0 2 α ρ 1+β sin sin , τ 2 2 τ p

as p → ∞.

(2.43)

...................................................

=

16


see (2.30), has the inverse Laplace and Fourier transforms by examining the convergence of the double integral (x ∈ R, t > 0, ε ∈ (0, 1]) s0 +i∞ ∞ 1 ˆ˜ (ξ , s) eiξ x dξ est ds Qε (x, t) = Q ε (2π )2 i s0 −i∞ −∞


We choose p0 so that (2.43) becomes

s=s0 +ip

∼ −p3 ,

as p → ∞.

Thus, for (2.42) we have |Jε+ (x, t)| ≤

∞ ∞ p0

0

ρ 1+β −(ερ)2 /4 e dρ dp < ∞. p3

Using the same arguments as for (2.40), we can prove that Jε− , given by (2.41), is also absolutely integrable. We proved that Qε , given by (2.38), has the inverse Laplace and Fourier transforms and therefore, by (2.37), we have 1 2 ˆ˜ (ξ , s), Kˆ˜ ε (ξ , s) = e−(εξ ) /4 − Q ε s

ξ ∈ R, Re s > s0 , ε ∈ (0, 1],

Kε (x, t) = H(t)δε (x) − Qε (x, t),

i.e.

x ∈ R, t > 0.

Thus, in (2.37) we can first invert the Laplace transform and subsequently the Fourier transform. The Fourier transform of the solution kernel Kˆ ε is obtained by the use of the inversion formula of the Laplace transform 1 Kˆ ε (ρ, t) = 2π i

s0 +i∞ s0 −i∞

Kˆ˜ ε (ρ, s) est ds,

a ≥ 0,

(2.44)

where Kˆ˜ ε is given by (2.37) and the complex integration along the contour Γ = Γ1 ∪ Γ2 ∪ Γr ∪ Γ3 ∪ Γ4 ∪ γ0 . The contour Γ is parametrized by Γ1 : s = R eiϕ , ϕ0 < ϕ < π ;

Γ2 : s = q eiπ , −R < −q < −r;

Γr : r eiϕ , −π < −ϕ < π ; Γ4 : s = R eiϕ , −π < ϕ < ϕ0 ;

Γ3 : s = q e−iπ , r < q < R; γ0 : s = s0 (1 + i tan ϕ), −ϕ0 < ϕ < ϕ0 ,

for arbitrary chosen R > 0 and 0 < r < R, and ϕ0 = arccos(s0 /R). By the Cauchy residues theorem and the results of lemma 2.3 we obtain: 1 (2.45) Kˆ˜ ε (ρ, s) est ds = Res(Kˆ˜ ε (ρ, s) est , sz (ρ)) + Res(Kˆ˜ ε (ρ, s) est , s¯z (ρ)). 2π i Γ Now, one shows (see, for example, [42] for similar calculations) that in (2.45), when R tends to infinity and r tends to zero, integrals along contours Γ1 , Γ4 and Γr tend to zero. The integrals along contours Γ2 and Γ3 in the limiting process (when R tends to infinity and r tends to zero) read (ρ ≥ 0, t > 0)

Kˆ˜ ε (ρ, s) est ds

lim

R→∞,r→0 Γ2

= −e−(εξ )

2

/4

∞ 0

q e−qt dq, q2 + ((1 + qα eiαπ )/(1 + τ qα eiαπ ))ρ 1+β sin(βπ/2)

Kˆ˜ ε (ρ, s) est ds

lim

R→∞,r→0 Γ3

= e−(εξ )

2

/4

∞ 0

q e−qt dq. q2 + ((1 + qα e−iαπ )/(1 + τ qα e−iαπ ))ρ 1+β sin(βπ/2)

...................................................



s3 + s

Im

17


By lemma 2.3, we have that the residues in (2.45) read (ρ ≥ 0, t > 0)

18

Res(Kˆ˜ ε (ρ, s) est , s¯z (ρ)) =

s est 2 e−(εξ ) /4 . 2s + (α(1 − τ )sα−1 /(1 + τ sα )2 )ρ 1+β sin(βπ/2) s=s¯z (ρ)

The integral along the contour γ0 in the limiting process tends to the integral on the right-hand side of (2.44) and therefore, putting all together in (2.45), we obtain 2 Kˆ ε (ρ, t) = S(ρ, t) e−(ερ) /4 ,

x ∈ R, t > 0,

with S given by (2.34). The inverse Fourier transform of such obtained Kˆ ε reads 1 ∞ 2 S(ρ, t) cos(ρx) e−(ερ) /4 dρ, x ∈ R, t > 0, Kε (x, t) = π 0 which in the distributional limit when ε → 0 gives the solution kernel K in the form (2.33).

3. Dependence of a solution on parameters α and β We examine the solutions in the limiting cases of the system of equations (1.11)–(1.13), or equivalently (1.17), subject to (1.19), (1.20) in view of remark 1.1. In all cases, we write the solution u of theorem 2.4 as u(x, t) = (u0 (x)δ(t) + v0 (x)H(t)) ∗x,t Kα,β (x, t),

x ∈ R, t > 0,

where the inverse Fourier transform of Kˆ˜ α,β , (2.29), is given by s 1 ∞ cos(ρx) dρ, K˜ α,β (x, s) = π 0 s2 + ((1 + sα )/(1 + τ sα ))ρ 1+β sin(βπ/2)

x ∈ R, Re s > 0.

(3.1)

Note that the integral in (3.1) is written formally and it denotes the inverse Fourier transform. As will be seen below, it may either converge or diverge, representing a distribution. We are interested in the behaviour of Kα,β for α and β tending to zero and one. We expect that the solution kernel Kα,β tends to solution kernels in specific cases. From the form of Kα,β given by (2.33), this cannot be easily seen. However, numerical examples (see §4) suggest that this holds true. What can be seen analytically is that the Laplace transform of Kα,β tends to the Laplace transforms of solution kernels in specific cases. When β → 0, then, in the sense of distributions, 1 ∞1 cos(ρx) dρ, x ∈ R, Re s > 0, K˜ α,0 (x, s) = π 0 s and the solution kernel Kα,0 is of the form Kα,0 (x, t) = δ(x)H(t),

x ∈ R, t > 0.

(3.2)

This is the case of the non-propagating disturbance, if the initial velocity is zero and the solution is given by (1.23) with v0 = 0. Therefore, regardless of the parameter α, when β tends to zero, the solution kernel tends to (3.2). This supports the idea from remark 1.1 (ii) that our system of equations can be useful in modelling materials which resist the propagation of the initial disturbance.

...................................................

s est 2 = e−(εξ ) /4 , 2s + (α(1 − τ )sα−1 /(1 + τ sα )2 )ρ 1+β sin(βπ/2) s=sz (ρ)


Res(Kˆ˜ ε (ρ, s) est , sz (ρ))


with

f+ (q) =

1 + τ qα eiαπ 1 + qα eiαπ

and f− (q) =

x ∈ R, t > 0,

(3.3)

1 + τ qα e−iαπ , q > 0. 1 + qα e−iαπ

We note that in [42] the solution is given in a slightly different form. When α → 0, then s 1 ∞ cos(ρx) dρ, K˜ 0,β (x, s) = π 0 s2 + (2/(1 + τ ))ρ 1+β sin(βπ/2) Using L−1 [s/(s2 + ω2 )](t) = cos(ωt) one easily comes to 2 1 ∞ βπ ρ 1+β sin K0,β (x, t) = cos t cos(ρx) dρ, π 0 1+τ 2

x ∈ R, Re s > 0.

x ∈ R, t > 0,

in the sense of distributions, which can be transformed to 1 1 ∞ βπ sin cos x + ct ρ K0,β (x, t) = 2π 0 2 ρ 1−β 1 βπ sin ρ dρ, x ∈ R, t > 0, + cos x − ct 2 ρ 1−β

(3.4)

where c = 2/(1 + τ ). This is the case of the space-fractional wave equation studied in [20]. Note that the solution in [20] is given in a different form. In this case, from (3.4), one can recover solution kernels for β = 0 and β = 1. If we put β = 0 in (3.4), we obtain 1 ∞ cos(xρ) dρ = δ(x), x ∈ R, t > 0, K0,0 (x, t) = π 0 i.e. (3.2). For β = 1 in (3.4), we obtain, in the sense of distributions, 1 ∞ (cos((x + ct)ρ) + cos((x − ct)ρ)) dρ, K0,1 (x, t) = 2π 0

x ∈ R, t > 0,

1 = (δ(x + ct) + δ(x − ct)), 2

with c =

2/(1 + τ ). This is the solution kernel for the classical wave equation.

4. Numerical examples We examine the qualitative properties of the solution to the space–time fractional Zener wave equation (1.17). Furthermore, we investigate the influence of the orders α and β of the, respective, time and space fractionalization of the constitutive equation and strain measure. Also, we numerically compare the solution to (1.17) with the solutions to the time-fractional Zener wave

...................................................

and the calculation similar to one presented in [42] leads to 1 ∞ (f− (q) e|x|qf− (q) − f+ (q) e|x|qf+ (q) ) e−qt dq, Kα,1 (x, t) = 4π i 0

19


When β → 1, we obtain the case of the time-fractional Zener wave equation, studied in [42]. In this case, s 1 ∞ cos(ρx) dρ, x ∈ R, Re s > 0, K˜ α,1 (x, s) = π 0 s2 + ((1 + sα )/(1 + τ sα ))ρ 2


u(x, t) 4

20

t = 0.5

t=2 x

0.5

1.0

1.5

2.0

2.5

–2

–4

Figure 1. Snapshots of the solution u(x, t) of (1.17) at t ∈ {0.5, 1, 1.5, 2} as a function of x. equation (1.26), which represents the limiting case β = 1 in (1.17). Both equations are subject to initial conditions u0 = δ, v0 = 0. In this case, the solution to (1.17), given by (2.32), (2.33), becomes u(x, t) = δ(x) ∗x K(x, t) = K(x, t),

x ∈ R, t > 0.

(4.1)

Since K, and therefore u, is a distribution in x, it cannot be plotted. Thus, we use the regularization Kε of the solution kernel so that (4.1) becomes 1 ∞ˆ 2 uε (x, t) = (4.2) K(ρ, t) e−(ερ) /4 cos(ρx) dρ, x ∈ R, t > 0. π 0 In all figures, we present the displacement field only on the half-axis x ≥ 0, as the field is symmetric with respect to the displacement axis. Figure 1 presents the plot of the displacement versus coordinate, obtained according to (4.2) for several time instants, while the other parameters of the model are α = 0.25, β = 0.45, τ = 0.1 and ε = 0.01. From figure 1, we see that, as time increases, the height of the peaks decreases, as the energy introduced by the initial disturbance field is being dissipated. This is the consequence of the viscoelastic properties of the material. In figure 2, we compare the displacements obtained as solutions for non-local (1.17) and local (1.26) wave equations, given by (4.2) and (3.3), respectively. Apart from β = 0.45 in (1.17) and β = 1 in (1.26) the other parameters in both models are as above. The effect of non-locality introduced in the strain measure is observed, as at a fixed time-instant, apart from the primary peak that exists in both models, in the non-local one there are secondary peaks in the displacement field. These secondary peaks reflect the influence of the oscillations of a certain material point on other material points in a medium. Thus, when the initial disturbance propagates (this is reflected by the existence of the primary peak), due to the non-locality, the secondary peaks reflect the residual influence of the disturbance transported by the primary peak. The non-locality changes not only the number of peaks at a certain time-instant but also the shape of the primary peak. The primary peak in the non-local model (compared with the local one) is higher and placed closer to the origin—the point where the initial Dirac-type disturbance field is introduced. The aim of the following figures is to show the influence of changing the non-locality parameter β. All other parameters are as above. Figure 3 presents plots of displacements for various values of β. From figure 3, one sees that, as the non-locality parameter β increases, the effects of non-locality decrease, as the height of the secondary peaks decreases and eventually the secondary peaks cease to exist. Also, as β increases, the height of the primary peak decreases and its position increases, being further from the point of the initial Dirac-type disturbance. In the limiting case β = 1, the displacement curve of non-local model (1.17) overlaps with the displacement curve of the local model (1.26). Figures 4 and 5 present the displacement field for

...................................................

0

t = 1.5


t=1

2


u(x, t) 4

21

t = 0.5

0

0.5

1.0

1.5

2.0

t=2 2.5

3.0

3.5

x

–2

–4

Figure 2. Snapshots of the solution u(x, t) at t ∈ {0.5, 1, 1.5, 2} as a function of x: (1.17), solid line, and (1.26), dashed line.

u(x, t) 2 b = 0.8 1 0

0.5

1.0

1.5

b = 0.95 b =1 x 2.0

b = 0.65 –1 b = 0.5 –2 –3

b = 0.35

Figure 3. Snapshots of the solution u(x, t) at t = 1 as a function of x: (1.17), solid line, and (1.26), dashed line.

u(x, t) 5

0

0.1

0.2

b = 0.2

–5

b = 0.15 –10

b = 0.1

Figure 4. Snapshots of the solution u(x, t) of (1.17) at t = 1 as a function of x.

0.3

x 0.4

...................................................

t = 1.5


t=1

2


u(x, t) 50

30 20

b = 0.01 b = 0.05

10 0

0.05

b = 0.1 0.10

x 0.20

0.15

–10

Figure 5. Snapshots of the solution u(x, t) of (1.17) at t = 1 as a function of x.

smaller values of β. In figure 4, one notices the significant influence of non-local effects. Namely, the secondary peaks are more prominent in height than the primary peak. Finally, in the limiting case when β → 0 one expects to obtain the displacement field in the form (1.23), i.e. u(x, t) = δ(x), x ∈ R, t > 0. This can be seen from figure 5. We might, thus, say that these numerical examples support the claim that β, as the non-locality parameter, measures the resistance of material to the disturbance propagation. Namely, at the same time-instant, as β decreases, the primary peak is placed closer to the point where the Dirac-type disturbance occurred and for β → 0 we obtain the non-propagating disturbance. Also, the shape of the primary peaks changes and, as β increases, they become more like the peak of the local model and for β = 1 the curves overlap. Acknowledgement The authors thank Radovan Obradović for his valuable suggestions concerning the calculations in numerical examples. Funding statement. This research is supported by the Serbian Ministry of Education and Science projects 174005, 174024, III44003 and TR32035, as well as by the Secretariat for Science of Vojvodina project 114-451-1084.

Appendix A. Mathematical background This section serves as a mathematical survey needed in the analysis that we have presented. We single out definitions and properties of fractional derivatives and, as our main tools are integral transforms, we recall the, more or less well known, main definitions and properties used. For a detailed exposition of the theory of fractional calculus see [14,15], and for the spaces and integral transforms we refer to [15,45]. Let 0 ≤ α < 1, −∞ ≤ a < b ≤ ∞. The left and right Caputo derivatives, of order α, of an absolutely continuous function u are defined by C α a Dt u(t) =

1 Γ (1 − α)

t a

u (θ ) dθ (t − θ )α

and

C α t Db u(t) = −

1 Γ (1 − α)

b t

u (θ) dθ, (θ − t)α

(A 1)

0 C 0 where Γ is the Euler gamma function and u = (d/dt)u. Note that C a Dt u(t) = t Db u(t) = u(t), and for 1 continuously differentiable functions and distributions we have that, as α → 1− , C a Dt u(t) → u (t), C D1 u(t) → −u (t). Therefore, the Caputo derivatives generalize integer order derivatives. t b Let 0 ≤ β < 1, −∞ ≤ a < b ≤ ∞. The symmetrized fractional derivative of an absolutely continuous function u is defined as b 1 1 u (θ) 1 C β C β C β dθ. a Eb u(x) = a Dx − x Db u(x) = 2 2 Γ (1 − β) a |x − θ|β

...................................................

b = 0.001


40

22


β

β

For a = −∞ and b = ∞, we write Ex instead of C a Eb and then

23

as fˇ (t) = f (−t), where H is the Heaviside function. Then f ∗ and fˇ ∗ are and {fˇα }α∈R ∈ S− α α α α convolution operators and for α < 0 they are operators of left and right fractional differentiation, α C α ˇ so that for u absolutely continuous we have C 0 Dt u = f1−α ∗ u and t Da u = −f1−α ∗ u . ˆ ϕ = u, ϕ , ˆ ϕ ∈ S(R), where for ϕ ∈ S For u ∈ S , the Fourier transform is defined as u, ∞ ϕ(ξ ˆ ) = F [ϕ(x)](ξ ) = ϕ(x) e−iξ x dx, ξ ∈ R. −∞

The Laplace transform of

u ∈ S

is defined by

˜ = L[u(t)](s) = F [e−ξ t u(t)](η), u(s)

s = ξ + iη.

It is well known that the function u˜ is holomorphic in the half-plane Re s > 0 (e.g. [45]). In particular, for u ∈ L1 (R) such that u(t) = 0, for t < 0, and |u(t)| ≤ A eat (a, A > 0) the Laplace transform is ˜ = 0∞ u(t) e−st dt, Re s > 0. u(s) We recall the main properties of the Fourier and Laplace transforms. Let u, u1 , u2 ∈ S F [u1 ∗ u2 ](ξ ) = F u1 (ξ ) · F u1 (ξ ), L[u1 ∗ u2 ](s) = Lu1 (s) · Lu1 (s),

F [u(n) ](ξ ) = (iξ )n F y(ξ ), n ∈ N, F δ(s) = 1, L[0 Dαt u](s) = sα Lu(s), α ≥ 0, Lδ(s) = 1,

where (·)(n) denotes the nth derivative. For β ∈ [0, 1) it holds that F [|x|−β ](ξ ) = 2Γ (1 − β) sin

βπ 1 2 |ξ |1−β

β

and F [Ex u(x)](ξ ) = i

ξ βπ ˆ ). u(ξ sin 1−β 2 |ξ |

References 1. Atanackovic TM, Pilipovic S, Zorica D. 2009 Time distributed-order diffusion-wave equation. I. Volterra type equation. Proc. R. Soc. A 465, 1869–1891. (doi:10.1098/rspa.2008.0445) 2. Atanackovic TM, Pilipovic S, Zorica D. 2009 Time distributed-order diffusion-wave equation. II. Applications of the Laplace and Fourier transformations. Proc. R. Soc. A 465, 1893–1917. (doi:10.1098/rspa.2008.0446) 3. Hanyga A. 2002 Multi-dimensional solutions of space-time-fractional diffusion equations. Proc. R. Soc. Lond. A 458, 429–450. (doi:10.1098/rspa.2001.0893) 4. Kochubei AN. 2008 Distributed order calculus and equations of ultraslow diffusion. J. Math. Anal. Appl. 340, 252–281. (doi:10.1016/j.jmaa.2007.08.024) 5. Mainardi F, Pagnini G, Gorenflo R. 2007 Some aspects of fractional diffusion equations of single and distributed order. Appl. Math. Comput. 187, 295–305. (doi:10.1016/j.amc.2006. 08.126) 6. Mainardi F, Pagnini G, Mura A, Gorenflo R. 2008 Time-fractional diffusion of distributed order. J. Vib. Control 14, 1267–1290. (doi:10.1177/1077546307087452) 7. Naber M. 2004 Distributed order fractional sub-diffusion. Fractals 12, 23–32. (doi:10.1142/ S0218348X04002410) 8. Atanackovic TM, Guran A. 2000 Theory of elasticity for scientists and engineers. Boston, MA: Birkhäuser.

...................................................

β

Note that Ex0 u(x) = 0 and Ex u(x) → u (x), as β → 1. So, the symmetrized fractional derivative generalizes the first derivative of a function. The zeroth order symmetrized fractional derivative of a function is zero (not a function itself). as For fractional operators in the distributional setting, one introduces a family {fα }α∈R ∈ S+ ⎧ tα−1 ⎪ ⎪ ⎨H(t) , α > 0, Γ (α) fα (t) = N ⎪ ⎪d f ⎩ α+N (t), α ≤ 0, α + N > 0, N ∈ N, dtN


1 1 β |x|−β ∗ u (x). Ex u(x) = 2 Γ (1 − β)


24 ...................................................


9. Meshkov SI. 1967 Description of internal friction in the memory theory of elasticity using kernels with a weak singularity. Zh. Prikl. Mekh. Tekh. Fiz. 8, 147–151. 10. Bagley R. 2007 On the equivalence of the Riemann-Liouville and the Caputo fractional order derivatives in modeling of linear viscoelastic materials. Fractional Calc. Appl. Anal. 10, 123–126. 11. Caputo M, Mainardi F. 1971 A new dissipation model based on memory mechanism. Pure Appl. Geophys. 91, 134–147. (doi:10.1007/BF00879562) 12. Mainardi F. 2010 Fractional calculus and waves in linear viscoelasticity. London, UK: Imperial College Press. 13. Rossikhin YuA. 2010 Reflections on two parallel ways in the progress of fractional calculus in mechanics of solids. Appl. Mech. Rev. 63, 010701. (doi:10.1115/1.4000246) 14. Podlubny I. 1999 Fractional differential equations. San Diego, CA: Academic Press. 15. Samko SG, Kilbas AA, Marichev OI. 1993 Fractional integrals and derivatives. Amsterdam, The Netherlands: Gordon and Breach. 16. Rossikhin YuA, Shitikova MV. 2010 Application of fractional calculus for dynamic problems of solid mechanics: novel trends and recent results. Appl. Mech. Rev. 63, 010801. (doi:10.1115/1.4000563) 17. Atanackovic TM, Konjik S, Oparnica Lj, Zorica D. 2011 Thermodynamical restrictions and wave propagation for a class of fractional order viscoelastic rods. Abstr. Appl. Anal. 2011, 32. (doi:10.1155/2011/975694) 18. Sumelka W, Blaszczyk T. 2014 Fractional continua for linear elasticity. Arch. Mech. 66, 147–172. 19. Truesdell C, Toupin RA. 1960 The classical field theories. In Encyclopedia of physics. Vol. III/1. Principles of classical mechanics and field theory (ed. S. Flügge), pp. 226–859. Berlin, Germany: Springer Verlag. 20. Atanackovic TM, Stankovic B. 2009 Generalized wave equation in nonlocal elasticity. Acta Mech. 208, 1–10. (doi:10.1007/s00707-008-0120-9) 21. Atanackovic TM, Konjik S, Oparnica Lj, Zorica D. 2012 The Cattaneo type space-time fractional heat conduction equation. Contin. Mech. Thermodyn. 24, 293–311. (doi:10.1007/ s00161-011-0199-4) 22. Eringen AC. 2002 Nonlocal continuum field theories. New York, NY: Springer. 23. Carpinteri A, Cornetti P, Sapora A. 2011 A fractional calculus approach to nonlocal elasticity. Eur. Phys. J. 193, 193–204. (doi:10.1140/epjst/e2011-01391-5) 24. Challamel N, Zorica D, Atanacković TM, Spasić DT. 2013 On the fractional generalization of Eringen’s nonlocal elasticity for wave propagation. C. R. Méc. 341, 298–303. (doi:10.1016/ j.crme.2012.11.013) 25. Di Paola M, Failla G, Pirrotta A, Sofi A, Zingales M. 2013 The mechanically based nonlocal elasticity: an overview of main results and future challenges. Phil. Trans. R. Soc. A 371, 20120433–1–16. (doi: 10.1098/rsta.2012.0433) 26. Di Paola M, Failla G, Zingales M. 2009 Physically-based approach to the mechanics of strong non-local linear elasticity theory. J. Elast. 97, 103–130. (doi:10.1007/s10659-009-9211-7) 27. Di Paola M, Zingales M. 2008 Long-range cohesive interactions of non-local continuum faced by fractional calculus. Int. J. Solids Struct. 45, 5642–5659. (doi:10.1016/j.ijsolstr.2008.06.004) 28. Cottone G, Di Paola M, Zingales M. 2009 Elastic waves propagation in 1D fractional non-local continuum. Phys. E 42, 95–103. (doi:10.1016/j.physe.2009.09.006) 29. Sapora A, Cornetti P, Carpinteri A. 2013 Wave propagation in nonlocal elastic continua modelled by a fractional calculus approach. Commun. Nonlinear Sci. Numer. Simul. 18, 63–74. (doi:10.1016/j.cnsns.2012.06.017) 30. Hanyga A, Seredynska ´ M. 2012 Spatially fractional-order viscoelasticity, non-locality, and a new kind of anisotropy. J. Math. Phys. 53, 052902. (doi:10.1063/1.4712300) 31. Hanyga A, Magin RL. 2014 A new anisotropic fractional model of diffusion suitable for applications of diffusion tensor imaging in biological tissues. Proc. R. Soc. A 470, 20140319. (doi:10.1098/rspa.2014.0319) 32. Hanyga A, Seredynska ´ M. 2012 Anisotropy in high-resolution diffusion-weighted MRI and anomalous diffusion. J. Magn. Reson. 220, 85–93. (doi:10.1016/j.jmr.2012.05.001) 33. Lazopulos KA. 2006 Non-local continuum mechanics and fractional calculus. Mech. Res. Commun. 33, 753–757. (doi:10.1016/j.mechrescom.2006.05.001) 34. Di Paola M. 2012 Exact mechanical models of fractional hereditary materials. J. Rheol. 56, 983–1004. (doi:10.1122/1.4717492) 35. Di Paola M, Paolo Pinnola F, Zingales M. 2013 A discrete mechanical model of fractional hereditary materials. Meccanica 48, 1573–1586. (doi:10.1007/s11012-012-9685-4)


25 ...................................................


36. Di Paola M, Paolo Pinnola F, Zingales M. 2013 Fractional differential equations and related exact mechanical models. Comput. Math. Appl. 66, 608–620. (doi:10.1016/j.camwa.2013.03.012) 37. Fabrizio M, Morro A. 1992 Mathematical problems in linear viscoelasticity, vol. 12. SIAM Studies in Applied Mathematics. Philadelphia, PA: SIAM. 38. Petrovic LjM, Zorica DM, Stojanac LjI, Krstonosic VS, Hadnadjev MS, Atanackovic TM. 2013 A model of the viscoelastic behavior of flowable resin composites prior to setting. Dental Mater. 29, 929–934. (doi:10.1016/j.dental.2013.06.005) 39. Atanackovic TM, Pilipovic S, Zorica D. 2012 On a system of equations arising in viscoelasticity theory of fractional type. (http://arxiv.org/abs/1205.5343) 40. Oparnica Lj. 2002 Generalized fractional calculus with applications in mechanics. Mat. Vesnik 53, 151–158. 41. Mainardi F, Gorenflo R. 2000 On Mittag-Leffler-type functions in fractional evolution processes. J. Comput. Appl. Math. 118, 283–299. (doi:10.1016/S0377-0427(00)00294-6) 42. Konjik S, Oparnica Lj, Zorica D. 2010 Waves in fractional Zener type viscoelastic media. J. Math. Anal. Appl. 365, 259–268. (doi:10.1016/j.jmaa.2009.10.043) 43. Näsholm SP, Holm S. 2013 On a fractional Zener elastic wave equation. Fractional Calc. Appl. Anal. 16, 26–50. (doi:10.2478/s13540-013-0003-1) 44. Treves F. 1975 Basic linear partial differential equations. New York, NY: Academic Press. 45. Vladimirov VS. 1984 Equations of mathematical physics. Moscow, Russia: Mir Publishers.

Diffraction-free beams in fractional Schrödinger equation.

Davey-Stewartson equation with fractional coordinate derivatives.

FRACTIONAL WAVE EQUATIONS WITH ATTENUATION.

Modified fractional variational iteration method for solving the generalized time-space fractional Schrödinger equation.

Approximate solution of time-fractional advection-dispersion equation via fractional variational iteration method.

Numerical simulation of fractional Cable equation of spiny neuronal dendrites.

finite element method for fractional diffusion equation.

An efficient technique for higher order fractional differential equation.

Transient aging in fractional Brownian and Langevin-equation motion.

Beam propagation management in a fractional Schrödinger equation.

A domain decomposition method for time fractional reaction-diffusion equation.

On the fractional Fourier and continuous fractional wave packet transforms of almost periodic functions.

Exact traveling wave solutions of modified KdV-Zakharov-Kuznetsov equation and viscous Burgers equation.

Hybrid Instantaneous Wave-Free Ratio-Fractional Flow Reserve versus Fractional Flow Reserve in the Real World.

On the stability of one-dimensional wave equation.

The pulsating orb: solving the wave equation outside a ball.

Special relativity derived from spacetime magma.

An efficient method for systems of variable coefficient coupled Burgers' equation with time-fractional derivative.

A parallel algorithm for the two-dimensional time fractional diffusion equation with implicit difference method.

Lie symmetry analysis and explicit solutions of the time fractional fifth-order KdV equation.

New Insights into the Fractional Order Diffusion Equation Using Entropy and Kurtosis.

A new Sumudu transform iterative method for time-fractional Cauchy reaction-diffusion equation.

Landau-Zener type surface hopping algorithms.

Propagation Dynamics of a Light Beam in a Fractional Schrödinger Equation.