EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 10, No. 3, 2017, 586-601 ISSN 1307-5543 – www.ejpam.com Published by New York Business Global Analytical Approximate Solution of Fractional Wave Equation by the Optimal Homotopy Analysis Method A. Elsaid1,∗, S. Shamseldeen1, S. Madkour1 1 Mathematics and Engineering Physics Department, Faculty of Engineering Mansoura University, Egypt Abstract. In this article, we study the space-fractional wave equation with Riesz fractional deriva- tive. The continuation of the solution of this space-fractional equation to the solution of the cor- responding integer order equation is proved. The series solution is obtained based on properties of Riesz fractional derivative operator and utilizing the optimal homotopy analysis method (OHAM). Numerical simulations are presented to validate the method and to show the effect of changing the fractional derivative parameter on the solution behavior. 2010 Mathematics Subject Classifications: 35L05, 26A33, 35C10 Key Words and Phrases: Space-fractional wave equation, Riesz, Optimal homotopy analysis method 1. Introduction Fractional derivatives, as generalizations of classical integer order derivatives, are in- creasingly used to model numerous problems in different fields of applied science. In recent years, the fractional derivative models are developed to describe the dissipative attenua- tion in complex materials, such as anomalous diffusion [12] and [15], viscoelastic damping [1] and [11], and wave propagation [4] and [5]. The operators of fractional differentiation and integration are also used for extensions of the diffusion and wave equations [13] and [14]. Studies have been devoted for a type of anomalous diffusion modeled by the frac- tional diffusion equation with spatial Riesz and Riesz-Feller fractional derivatives [6] and [8]. Yet, few articles dealt with applying iterative techniques to Riesz fractional partial differential equations (FPDEs). This is due to the difficulty in repeated application of Riesz fractional derivative to solution components. This work is based on properties that show repetitive behavior for complex exponential function, hence sine and cosine functions, when subjected to the application of Riesz fractional derivative [6] and [7]. In this work, the motivation is to establish the continuation of the solution of the space- fractional wave equation with spatial derivative in Riesz sense to the exact solution of the ∗Corresponding author. http://www.ejpam.com 586 c© 2017 EJPAM All rights reserved. A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 587 corresponding integer-order equation as the order of the fractional derivative approaches its integer limit. This objective is carried out theoretically then via approximate series solution obtained iteratively by applying the optimal homotopy analysis method (OHAM). We consider the space-fractional wave equation of the form ∂2 ∂t2 u(x, t) = Rαxu(x, t) + P (u), −∞ < x <∞ , t > 0, (1) subject to the initial conditions  u(x, 0) = f1(x), ∂ ∂tu(x, 0) = f2(x). (2) where Rαx denotes the Riesz fractional derivative (in space) of order α. The parameter α is restricted to the conditions 0 < α < 2 and α 6= 1. The function P is a continuous function in u, and the two functions f1 and f2 are functions in the space of integrable functions L1(−∞,∞). This paper is organized as follows. In Section two, basic definitions of fractional deriva- tive operators involved are presented. Proof of continuation of solution is presented in Section three. The OHAM is illustrated in Section four. In Section five, the results of nu- merical experiments are presented, considering the space fractional sine-Gordan equation. Section six contains the conclusion of this work. 2. Fractional derivatives and integrals Definition 1. A real function f(x), x > 0, is said to be in the space Cµ, µ ∈ R, if there exists a real number p > µ, such that f(x) = xpf1(x), where f1(x) ∈ C(0,∞), and it is said to be in the space Cmµ if fm ∈ Cµ, m ∈ N. Definition 2. The Riemann-Liouville fractional integral operator of order α ≥ 0 of a function f(x) ∈ Cµ, µ ≥ −1 is defined as Jαf(x) = 1 Γ(α) x∫ 0 (x− τ)α−1f(τ)dτ, α > 0, x > 0, J0f(x) = f(x). (3) Definition 3. The fractional derivative in Riemann-Liouville sense of f(x),m ∈ N, x > 0 is defined as Dβ xf(t) = dm dxm Jm−βf(x), m− 1 < β < m. (4) Definition 4. The fractional derivative in Caputo sense of f(x) ∈ Cm−1, m ∈ N, x > 0 is defined as CDβ xf(x) = { Jm−β dm dxm f(x), m− 1 < β < m, dm dxm f(x), β = m. (5) A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 588 Definition 5. The Riesz partial fractional derivative Rαx is defined as [8] Rαxu(x) = − 1 2 cos(απ/2) [Dα +u(x) +Dα −u(x)], 0 < α < 2, α 6= 1 (6) where Dα ±u(x) are the Weyl fractional derivatives Dα ±u(x) = { ± d dxW 1−α ± u(x), 0 < α < 1 d2 dx2 W 2−α ± u(x), 1 < α < 2 , (7) and W β ± denote the Weyl fractional integrals of order β > 0, given by W β +u(x) = 1 Γ(β) x∫ −∞ (x− z)β−1u(z)dz, W β −u(x) = 1 Γ(β) ∞∫ x (z − x)β−1u(z)dz. (8) When α = 0 the Weyl fractional derivative degenerates into the identity operator D0 ±u(x) = u(x). (9) For continuity we have D1 ±u(x) = ± d dx u(x), D2 ±u(x) = d2 dx2 u(x). (10) Evidently, in case α = 2, we define Rαxu(x) = d2 dx2 u(x). (11) For the case α = 1 we have R1 xu(x) = d dx Hu(x) (12) = d dx 1 π ∞∫ −∞ u(z) z − x dz, (13) where H is the Hilbert transform and the integral is understood in the Cauchy principal value sense. 3. Continuation of the solution In this section, we prove the continuation of the solution to fractional-order wave equation with Riesz spatial derivative to the solution of the corresponding integer-order equation. A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 589 Theorem 1. If f1(x) and f2(x) are functions in the space of integrable functions L1(−∞,∞), then the exact solution uα(x, t) of the space fractional wave equation ∂2 ∂t2 u(x, t) = Rαxu(x, t), −∞ < x <∞ t > 0, (14) with the initial conditions  u(x, 0) = f1(x), ∂ ∂tu(x, 0) = f2(x). (15) is given by uα(x, t) = 1 π ∞∫ −∞ ∞∫ 0 ( E2,1(−ωα t2)f1(v) + t E2,2(−ωα t2)f2(v) ) cos(ω(x− v))dωdv (16) where Eη,γ(z) is the Mittage Leffler function defined by [16] Eη,γ(z) = ∞∑ n=0 zn Γ(ηn+ γ) , (17) where E2,1(−ωα t2) = cos(ωα/2 t), (18) E2,2(−ωα t2) = sin(ωα/2 t) ωα/2t . (19) Theorem 2. Let α ∈ (1, 2), f1(x) and f2(x) are functions in the space of integrable functions L1(−∞,∞), and uα displayed in (16) be the solution of the space-fractional problem (14-15) , then lim α→2 uα(x, t) = u(x, t), where u(x, t) is the exact solution of the integer-order wave equation{ utt(x, t) = uxx(x, t), −∞ < x <∞, t > 0, u(x, 0) = f1(x), ut(x, 0) = f2(x). (20) Proof. Consider the set of functions ϕn(ω) and ψn(ω) for ω ∈ (0,∞), n ∈ N+ by ϕn(ω) = 1 π E2,1(−ω2− 1 n+1 t2) ∞∫ −∞ f1(v) cos(ω(x− v))dv, (21) A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 590 ψn(ω) = 1 π tE2,2(−ω2− 1 n+1 t2) ∞∫ −∞ f2(v) cos(ω(x− v))dv. (22) These two set of functions satisfy Lebesgue dominated convergence theorem as |ϕn(ω)| ≤ 1 π ∣∣∣E2,1(−ω2− 1 n+1 t2) ∣∣∣ ∞∫ −∞ |f1(v)| |cos(ω(x− v))| dv, ≤ 1 π ∣∣∣E2,1(−ω2− 1 n+1 t2) ∣∣∣ ∞∫ −∞ |f1(v)| dv, and since f1 ∈ L1(−∞,∞), there exists M > 0 such that ∞∫ −∞ |f1(v)| dv < M . Hence |ϕn(ω)| ≤ M π ∣∣∣E2,1(−ω2− 1 n+1 t2) ∣∣∣ . (23) From [16] Theorem (1.6), there exits K1 > 0 such that |Eη,γ(−z)| ≤ K1 1 + |z| , (24) then |ϕn(ω)| ≤ MK1 π 1 1 + ∣∣∣ω2− 1 n+1 t2 ∣∣∣ , ω ∈ (0,∞), n = 1, 2, ... (25) For bounded time interval 0 < t < T <∞, there exists K2(ρ) > 0 such that |ϕn(ω)| ≤ g1(ω) = K2(ρ) 1 + ω1+ρ , ρ ∈ (0, 0.5), and g1(ω) ∈ L1(0,∞) since ∞∫ 0 |g1(ω)| dω = K2(ρ)Γ( ρ 1 + ρ )Γ(1 + 1 1 + ρ ). (26) Thus the set of functions ϕn(ω) satisfy Lebesgue dominated convergence theorem. Following the same steps, one can prove that the set of functions ψn(ω) satisfy Lebesgue dominated convergence theorem as well. Now, as lim n→∞ ϕn(ω) = 1 π E2,1(−ω2t2) ∞∫ −∞ f1(v) cos(ω(x− v))dv, (27) A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 591 lim n→∞ ψn(ω) = 1 π tE2,2(−ω2t2) ∞∫ −∞ f2(v) cos(ω(x− v))dv. (28) then setting α = 2− 1 n+1 u2(x, t) = lim α→2 uα(x, t) = lim n→∞ ∞∫ 0 [ϕn(ω) + ψn(ω)]dω (29) = ∞∫ 0 lim n→∞ [ϕn(ω) + ψn(ω)]dω, (30) which yields u2(x, t) = 1 π ∞∫ −∞ ∞∫ 0 ( cos(ω t)q1(v) + sin(ω t) ω q2(v) ) cos(ω(x− v))dωdv, which is the exact solution of the integer-order wave equation (20). 4. Optimal homotopy analysis method (OHAM) We begin by illustrating the classical homotopy analysis method (HAM). Consider the following nonlinear equation N [u(x, t)] = 0, (31) where N is a nonlinear operator, u(x, t) is the unknown function and x and t denote spatial and temporal independent variables, respectively. By generalizing the traditional homotopy method, Liao [9] constructs the so-called zero-order deformation equation (1− p)L[φ(x, t; p)− u0(x, t)] = p~H(x, t)N [φ(x, t; p)], (32) where p ∈ [0, 1] is an embedding parameter, ~ is a nonzero auxiliary parameter, H(x, t) is an auxiliary function, L is an auxiliary linear operator, u0(x, t) is an initial guess of u(x, t) and φ(x, t; p) is an unknown function. Obviously, when p = 0 and p = 1, we have φ(x, t; 0) = u0(x, t), φ(x, t; 1) = u(x, t), respectively. Thus, as p increases from 0 to 1, the solution φ(x, t; p) varies from the initial guess u0(x, t) to the solution u(x, t). By expanding φ(x, t; p) in Taylor series with respect to p, we have φ(x, t; p) = u0(x, t) + ∞∑ m=1 um(x, t)pm, (33) A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 592 where um(x, t) = 1 m! ∂mφ(x, t; p) ∂pm |p=0 . (34) If the auxiliary linear operator, the initial guess and the auxiliary parameter ~ and the auxiliary function are so properly chosen, then, as proved by Liao [9], series (33) converges at p = 1 and one has u(x, t) = u0(x, t) + ∞∑ m=1 um(x, t) (35) which must be one of solutions of the original nonlinear equation, as proved by Liao [9]. Using definition (34), the governing equation of the HAM can be deduced from the zero-order deformation equation (32) as follows. Define the vector −→u n = {u0(x, t), u1(x, t), u2(x, t), ..., un(x, t)} (36) From equation (32), the so-called m th-order deformation equation is given by L[um(x, t)− χmum−1(x, t)] = ~H(x, t) 1. (39) Applying the inverse operator L−1 to both sides of (37), um(x, t) can be easily solved for by symbolic computations software. The HAM has been successfully applied to solve various classes of equations and applied problems [3]-[2]. In the classical HAM, choosing the value of parameter ~ depends on inspecting the graph of the quantity of interest; the solution or one of its derivatives. Yet, when H(x, t) is fixed, it is obvious that um(x, t) contains only one control parameter ~. Thus, by con- structing a formula for the residual error, the OHAM solution is obtained by choosing the value for parameter ~ that minimizes the error. Here, the averaged residual error defined for ordinary differential equations in [10] is generalized to the case of two variable partial differential equations in the following form Em(~) = 1 MK M∑ i=0 K∑ j=0 [ N m∑ n=0 un ( i M , j K )]2 , (40) which is a nonlinear algebraic equation of one unknown; the convergence-control parameter ~. Thus the optimal value of ~ is determined by the minimum of the averaged residual error Em to ensure the fast convergence of the homotopy series. To apply the OHAM recursive technique to the problem, a repeated evaluation of Riesz fractional derivative to solution components is needed. This obstacle is overcome by using property of Riesz fractional derivative in the following lemma. A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 593 Lemma 3. Let α ∈ (0, 2), α 6= 1. Then Rαx(eiωx) = −ωαei(ωx), (41) or in a trigonometric form Rαx sin(ωx) = −ωα sin(ωx), (42) Rαx cos(ωx) = −ωα cos(ωx). (43) Proof. See [6] and [7]. 5. Numerical simulation In this section, we consider linear and nonlinear problems to illustrate the efficiency of the method of solution to this type of problems and to illustrate the continuation of the solution we proved in Section 3. In each problem, a table is presented to show the estimated values the optimal convergence control parameter ~ and the corresponding residual error Em at different values of the fractional derivative α. These estimated values are calculated via minimizing of the averaged residual error Em displayed in (40) in the space domain 0 ≤ x ≤ 2.0 and the time interval 0 ≤ t ≤ 2.0. Example 1. Consider problem (1-2) with p(u) = u, f1(x) = sin(πx/a) and f2(x) = − sin(πx/a) { utt(x, t) = Rαxu(x, t) + u, −∞ < x <∞, t > 0, u(x, 0) = sin(πx/a), ut(x, 0) = sin(πx/a), (44) where a is a real constant. The auxiliary linear operator is chosen as L[φ] = ∂2 ∂t2 (φ), (45) and the nonlinear operator N is chosen as N [φ] = φtt −Rαx(φ)− φ. (46) The m th-order deformation equation, with H(x, t) = 1, for this linear problem is given by ∂2 ∂t2 [um(x, t)− χmum−1(x, t)] = ~ ( ∂2 ∂t2 (um−1)−Rαx(um−1)− um−1 ) , (47) A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 594 with u0(x, t) = f1(x) + t f2(x) = (1− t) sin(πx/a) (48) The inverse integral operator is applied to both sides of equation (47) to obtain the series solution terms. The first three terms are given by u0 = (1− t) sin( πx a ), u1 = −ht 2 6 ( −1 + (π a )α) (−3 + t) sin (πx a ) , u2 = − ht 2 120 ( −1 + (π a )α) ( 20(−3 + t) + h ( −60 + 20t+ ( 5− 5 (π a )α) t2 + ( −1 + (π a )α) t3 )) sin (πx a ) . Table 1 shows the estimated values of the optimal convergence control parameter ~ and the corresponding residual error Em for the linear problem displayed in (44) at different values of the fractional derivative α in the space domain 0 ≤ x ≤ 2.0 and the time interval 0 ≤ t ≤ 2.0. Table 1: The estimated optimal convergence parameter ~ and the corresponding residual error Em for 0 ≤ x ≤ 2.0 and 0 ≤ t ≤ 2.0 at different fractional derivative α for Example (1). α ~ Em Optimal parameter Residual Error 1.7 −0.940496 1.12317E − 5 1.8 −0.938046 8.19812E − 5 1.9 −0.933025 1.75198E − 4 2.0 −0.928713 3.51187E − 4 Α = 1.7 Α = 1.9 Α = 1.8 Α = 2.0 0.5 1.0 1.5 2.0 x 0.02 0.04 0.06 0.08 uΑHx,0.7L Figure 1: The solution of (44) at t = 0.5, 0 ≤ x ≤ 2 and different values of the fractional order α = 1.7, 1.8, 1.9 and 2.0. The series solution is obtained by u = u0 +u1 +u2 +u3 + ..... Figures (1) and (2) show the effect of the fractional order derivative α on the behavior of the solution at fixed time A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 595 Α = 1.7 Α = 1.8 Α = 2.0 Α = 1.9 0.5 1.0 1.5 2.0 x -0.4 -0.3 -0.2 -0.1 uΑHx,1.0L Figure 2: The solution of (44) at t = 1.0, 0 ≤ x ≤ 2 and different values of the fractional order α = 1.7, 1.8, 1.9 and 2.0. t = 0.5 t = 1.0 t = 0.0 t = 1.5 0.5 1.0 1.5 2.0 x -1.0 -0.5 0.5 1.0 u1.9Hx,tL Figure 3: The solution of (44) at different times t = 0.0, 0.5, 1.0, and 1.5, 0 ≤ x ≤ 2 and the fractional order α = 1.9. t = 0.5 and t = 1.0, respectively, while Figure (3) illustrates the temporal behavior of the solution at a fixed fractional order, α = 1.9. The plots represent the sum of the first four terms (u0 to u3) in the OHAM series when a = 2.0 Example 2. Consider problem (1-2) with P (u) = u + c u3, f1(x) = sin(πx/a) and f2(x) = − sin(πx/a) { utt(x, t) = Rαxu(x, t) + u + c u3, −∞ < x <∞, t > 0, u(x, 0) = sin(πx/a), ut(x, 0) = sin(πx/a), (49) where a is a constant. The auxiliary linear operator is chosen as L[φ] = ∂2 ∂t2 (φ), (50) and the nonlinear operator N is chosen as N [φ] = φtt −Rαx(φ)− φ− c φ3. (51) Then, m th-order deformation equation for this problem is given by ∂2 ∂t2 [um(x, t)− χmum−1(x, t)] = ~H(x, t) 0, u(x, 0) = π + ε cos(µx), ut(x, 0) = 0, (55) where ε and µ are real constants. Here the auxiliary linear operator is L[φ] = ∂2 ∂t2 (φ), (56) A. Elsaid, S. Shamseldeen, S. Madkour / Eur. J. Pure Appl. Math, 10 (3) (2017), 586-601 598 and the nonlinear operator N is chosen as N [φ] = φtt −Rαx(φ) + sin(φ). (57) Then, m th-order deformation equation for this problem is given by ∂2 ∂t2 [um(x, t)− χmum−1(x, t)] = ~H(x, t)