A spectral tau algorithm based on Jacobi operational matrix for ...
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.1 (1-15)
Journal of Computational Physics ••• (••••) •••–•••
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A spectral tau algorithm based on Jacobi operational matrix for numerical solution of time fractional diffusion-wave equations A.H. Bhrawy a,b , E.H. Doha c , D. Baleanu d,e,f,∗ , S.S. Ezz-Eldien g a
Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt d Department of Chemical and Materials Engineering, Faculty of Engineering, King Abdulaziz University, Jeddah, Saudi Arabia e Department of Mathematics and Computer Sciences, Cankaya University, Ankara, Turkey f Institute of Space Sciences, Magurele-Bucharest, R 76900, Romania g Department of Basic Science, Institute of Information Technology, Modern Academy, Cairo, Egypt b c
a r t i c l e
i n f o
Article history: Received 11 January 2014 Received in revised form 7 March 2014 Accepted 23 March 2014 Available online xxxx Keywords: Fractional diffusion-wave equations Tau method Shifted Jacobi polynomials Operational matrix Caputo derivative
a b s t r a c t In this paper, an efficient and accurate spectral numerical method is presented for solving second-, fourth-order fractional diffusion-wave equations and fractional wave equations with damping. The proposed method is based on Jacobi tau spectral procedure together with the Jacobi operational matrix for fractional integrals, described in the Riemann– Liouville sense. The main characteristic behind this approach is to reduce such problems to those of solving systems of algebraic equations in the unknown expansion coefficients of the sought-for spectral approximations. The validity and effectiveness of the method are demonstrated by solving five numerical examples. Numerical examples are presented in the form of tables and graphs to make comparisons with the results obtained by other methods and with the exact solutions more easier. © 2014 Elsevier Inc. All rights reserved.
1. Introduction Due to the accuracy of ordinary and partial fractional differential equations in modeling various engineering and physical phenomena such as bioengineering [25], solid mechanics [36], anomalous transport [28], continuum and statistical mechanics [26], fluid-dynamic [16], economics [1], nonlinear oscillation of earthquakes [15], colored noise [27] and many other problems [33], fractional calculus has become the focus of many researchers in recent years. The time-fractional diffusion-wave equation, which is a mathematical model of a wide class of important physical phenomena, is a linear integro-partial fractional differential equation that obtained from the classical diffusion-wave equation by replacing the second-order time derivative term by a fractional derivative of order ν , 1 < ν 2. Many of the universal mechanical, acoustic and electromagnetic responses may be described accurately by the time-fractional diffusion-wave equation, see [31,32]. The fourth-order space derivative is arising in the wave propagation in beams and modeling formation of grooves on a flat surface, thus considerable attention has been devoted to fourth-order fractional diffusion-wave equation and its applications, see [33].
*
Corresponding author.
http://dx.doi.org/10.1016/j.jcp.2014.03.039 0021-9991/© 2014 Elsevier Inc. All rights reserved.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.2 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
2
Several numerical methods have been proposed in the last few years for solving such equations. In [6,7], Cui presented high-order compact finite difference methods to solve the fractional diffusion and time fractional diffusion equation. A compact difference scheme [17,12] is proposed for solving second- and fourth-order fractional diffusion-wave equations. Chen et al. [5] proposed and developed the method of separation of variables with constructing the implicit difference scheme for fractional diffusion-wave equation with damping. Recently, Hu and Zhang [18] proposed the finite difference scheme for approximating the solution of the fourth-order fractional diffusion-wave equation. An explicit analytical solution is proposed in [14] for space–time fractional wave equation. In [8], Darzi et al. solved fractional diffusion-wave equations, by using the Sumudu transform method. Some other numerical methods are introduced for solving fractional differential equations with sub-diffusion and super-diffusion, see [13,21,22,4,20,35,40]. Spectral methods are a class of important tools for obtaining the numerical solutions of fractional differential equations. They have excellent error properties and they offer exponential rates of convergence for smooth problems [3,38,39]. Considerable attention has been devoted to the use of orthogonal polynomials aiming to deal with various fractional differential equations, by reducing them to systems of algebraic equations that greatly simplify obtaining their solutions. Moreover, the operational matrices of fractional derivatives are derived and used together with spectral methods based on orthogonal polynomials for developing numerical approximations to fractional differential equations, see [9,10]. Recently, and with the help of operational matrices of fractional derivatives for orthogonal polynomials, the tau spectral method is also utilized to solve spectrally space-fractional diffusion equations, see [11,34,37]. In the current work, we investigate and develop some efficient spectral algorithms to solve spectrally second- and fourthorder fractional diffusion-wave equations and fractional diffusion-wave equation with damping. These algorithms are based principally on shifted Jacobi tau-spectral method combined with the generalized shifted Jacobi operational matrix of integrals. The current article is organized as follows. In Section 2, the shifted Jacobi polynomials with their properties and their operational matrix to fractional integration are presented. In Sections 3, 4 and 5, the shifted Jacobi operational matrix of fractional integration is applied together with the shifted Jacobi tau spectral method to give numerical solutions for secondand fourth-order fractional diffusion-wave equations and fractional diffusion-wave equation with damping, respectively. In Section 6, several numerical examples and comparisons between our obtained numerical results and those of other methods are considered. Also a conclusion is given in Section 7. 2. Operational matrices of shifted Jacobi polynomials Riemann–Liouville and Caputo fractional definitions are the two most used from all the other definitions of fractional calculus which have been introduced recently. Definition 2.1. The integral of order
I ν f (x) =
x
1
Γ (ν )
ν 0 (fractional) according to Riemann–Liouville is given by
(x − t )ν −1 f (t ) dt ,
ν > 0, x > 0,
0
I 0 f (x) = f (x),
(2.1)
where
∞ Γ (ν ) =
xν −1 e −x dx
0
is gamma function. The operator J ν satisfies the following properties
I ν I μ f (x) = I ν +μ f (x), I ν I μ f (x) = I μ I ν f (x), Γ (β + 1) β+ν I ν xβ = x . Γ (β + 1 + ν )
(2.2)
Definition 2.2. The Caputo fractional derivative of order
D ν f (x) =
1
Γ (m − ν )
x
(x − t )m−ν −1
dm dt m
f (t ) dt ,
0
where m is the ceiling function of ν . The operator D ν satisfies the following properties
ν is defined by m − 1 < ν m , x > 0,
(2.3)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.3 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
3
D ν I ν f (x) = f (x), I ν D ν f (x) = f (x) −
m −1
f ( i ) 0+
xi i!
i =0
D ν xβ =
,
Γ (β + 1) β−ν x . Γ (β + 1 − ν )
(2.4) (α ,β)
The Jacobi polynomials, denoted by P j
(z); α > −1, β > −1, constitute an orthogonal system with respect to the = [−1, 1], i.e.,
weight function w (α ,β) ( z) = (1 − z)α (1 + z)β over I
1
(α ,β)
Pj
(α ,β)
( z) P k
(α ,β)
( z) w (α ,β) ( z) dz = δ jk γk
,
−1
where δ jk is the Kronecker function and
γk(α ,β) =
2α +β+1 Γ (k + α + 1)Γ (k + β + 1)
(2k + α + β + 1)k!Γ (k + α + β + 1)
.
In order to use Jacobi polynomials on the interval x ∈ [0, L ] we define the so-called shifted Jacobi polynomials by in2x (α ,β) 2x (α ,β) − 1. Let the shifted Jacobi polynomials P j ( − 1) be denoted by P L , j (x). L L (α ,β) (x) = Then the shifted Jacobi polynomials are constituting an orthogonal system with respect to the weight function w L xβ ( L − x)α over I = [0, L ] with the orthogonality property
troducing the change of variable z =
L
(α ,β)
P L, j
(α ,β)
(α ,β)
(x) P L ,k (x) w L
(α ,β)
(x) dx = h L ,k ,
(2.5)
0
where (α ,β)
h L ,k
=
α +β+1 L
(α ,β)
δ jk γ j
2
=
L α +β+1 Γ (k + α + 1)Γ (k + β + 1)
(2k + α + β + 1)k!Γ (k + α + β + 1) (α ,β)
The explicit analytic form of the shifted Jacobi polynomials P L ,i (α ,β)
P L ,i
(x) =
i
(−1)i −k
k =0
δ jk .
(x) of degree i is given by
Γ (i + β + 1)Γ (i + k + α + β + 1) xk , Γ (k + β + 1)Γ (i + α + β + 1)(i − k)!k! Lk
(2.6)
and this in turn, enables one to get
Γ (i + β + 1)(i + α + β + 1)q Γ ( i + β + 1) (α ,β) D q P L ,i (0) = (−1)i −q q , , Γ (β + 1)i ! L Γ (i − q + 1)Γ (q + β + 1) Γ (i + α + 1)(i + α + β + 1)q Γ ( i + α + 1) (α ,β) (α ,β) P L ,i ( L ) = D q P L ,i ( L ) = q , , q i, Γ (α + 1)i ! L Γ (i − q + 1)Γ (q + α + 1) (α ,β)
P L ,i
(0) = (−1)i
q i,
where (a)q = Γ (a + q)/Γ (a), and q is an integer, these relations will be of important use later.
(α ,β)
Assume u (x) is a square integrable function with respect to the Jacobi weight function w L expressed in terms of shifted Jacobi polynomials as
u (x) =
∞
(α ,β)
c j P L, j
in (0, L ), then it can be
(x),
j =0
from which the coefficients c j are given by
1 c j = (α ,β) h L, j
L
(α ,β)
wL
(α ,β)
(x)u (x) P L , j
dx,
j = 0, 1 , . . . .
0
If we approximate u (x) by the first ( N + 1)-terms, then we can write
(2.7)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.4 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
4
N
u N (x) =
(α ,β)
c j P L, j
(x),
(2.8)
j =0
which alternatively may be written in the matrix form:
u N (x) = C T Φ L , N (x),
C T ≡ [c 0 , c 1 , . . . , c N ],
(2.9)
with
(α ,β) T (α ,β) (α ,β) Φ L , N (x) ≡ P L ,0 (x), P L ,1 (x), . . . , P L , N (x) .
(2.10)
Similarly, let u (x, t ) be an infinitely differentiable function defined on 0 < x L and 0 < t τ . Then it is possible to express as
u M , N (x, t ) =
N M
(α ,β)
a i j P τ ,i
(α ,β)
(t ) P L , j (x) = ΦτT, M (t )AΦ L , N (x),
(2.11)
i =0 j =0
with
⎛
a00 ⎜ a10
a01 a11
A=⎜ ⎝ ..
.
.. .
a M0
a M1
⎞
··· ···
a0N a1N ⎟
.. ⎟ ⎠, . ··· · · · aM N
and
1 ai j = (α ,β) (α ,β) h τ ,i h L , j
τL
(α ,β)
u (x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
(x) dx dt ,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N .
(2.12)
The first derivative of the vector Φ L , N (x) can be expressed by
dΦ L , N (x) dx
= DΦ L , N (x),
(2.13)
where D is the ( N + 1) × ( N + 1) operational matrix of derivative given by
D = (di j ) =
ξ(i , j ), i > j , 0
otherwise,
and
ξ(i , j ) =
L α +β (i + α + β + 1)(i + α + β + 2) j ( j + α + 2)i − j −1 Γ ( j + α + β + 1)
× 3 F2
(i − j − 1)!Γ (2 j + α + β + 1) − i + 1 + j , i + j + α + β + 2, j + α + 1 j + α + 2,
2j +α +β +2
;1 .
(For the general definition of a generalized hypergeometric series and special 3 F 2 , see [24], p. 41 and pp. 103–104, respectively.) Repeated use of (2.13) enables one to write
dq Φ L , N (x) dxq
= Dq Φ L , N (x),
(2.14)
where q is a natural number and Dq means matrix powers. The fractional integration of order ν of φτ , M (t ) can be expressed as
I ν Φτ , M (t ) = Pν Φτ , M (t ), where Pν is the ( M + 1) × ( M + 1) operational matrix of fractional integration of order and is defined by
(2.15)
ν in the Riemann–Liouville sense
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.5 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
⎛ Ω (0, 0, α , β) Ω (0, 1, α , β) ν ν ⎜ Ων (1, 0, α , β) Ων (1, 1, α , β) ⎜ .. .. ⎜ . . ⎜ Pν = ⎜ Ων (i , 1, α , β) ⎜ Ων (i , 0, α , β) ⎜ .. .. ⎝ . . Ων ( M , 0, α , β) Ων ( M , 1, α , β)
··· ··· ··· ··· ··· ···
Ων (0, M , α , β) Ων (1, M , α , β) .. .
5
⎞
⎟ ⎟ ⎟ ⎟ ⎟, Ων (i , M , α , β) ⎟ ⎟ .. ⎠ . Ων ( M , M , α , β)
where
Ων (i , j , α , β) =
i (−1)i −k Γ (i + β + 1)Γ (i + k + α + β + 1)Γ (α + 1)τ ν
Γ (k + β + 1)Γ (i + α + β + 1)(i − k)!Γ (k + ν + 1)
k =0
×
j (−1) j − f Γ ( j + f + α + β + 1)Γ ( f + k + β + ν + 1)(2 j + α + β + 1) j ! Γ ( j + α + 1)Γ ( f + β + 1)( j − f )! f !Γ ( f + k + α + β + ν + 2) f =0
(see [2] for proof). 3. Second-order diffusion-wave equation In this section, the spectral tau method together with the shifted Jacobi operational matrix of integrations are applied to solve second-order fractional diffusion-wave equations. Consider the fractional diffusion-wave equation of the form:
∂ 2 u (x, t ) ∂ ν u (x, t ) = c ( x ) , ∂tν ∂ x2
0 < x L , 0 < t τ , 1 < ν 2,
(3.1)
with initial conditions
u (x, 0) = f 0 (x),
ut (x, 0) = f 1 (x),
(3.2)
u ( L , t ) = g 1 (t ).
(3.3)
and boundary conditions
u (0, t ) = g 0 (t ),
If we apply the Riemann–Liouville integral of order Eq. (3.1), namely
ν on Eq. (3.1) and making use of (2.4), we get the integrated form of
∂ 2 u (x, t ) , ∂ x2 u ( L , t ) = g 1 (t ),
u (x, t ) − q(x, t ) = I tν c (x) u (0, t ) = g 0 (t ),
(3.4)
where q(x, t ) = t f 1 (x) + f 0 (x) and I tν is the fractional integration of order
ν with respect to t. In order to use the pure spectral tau method with the shifted Jacobi operational matrix for solving the fully integrated problem (3.4), we assume that c (x) and q(x, t ) are smooth and continuous functions. Therefore, u (x, t ), c (x) and q(x, t ) can be expanded by the shifted Jacobi polynomials as u M , N (x, t ) = ΦτT, M (t )AΦ L , N (x), c N (x) = C T Φ L , N (x), q M , N (x, t ) = ΦτT, M (t )QΦ L , N (x),
(3.5)
where A is the unknown coefficient ( M + 1) × ( N + 1) matrix, C T and Q are known matrices that can be written as
C T = [c 0 , c 1 , . . . , c N ], ⎛q q01 · · · 00 ⎜ q10 q11 · · · Q=⎜ . ⎝ .
..
..
q M0
q M1
⎞
q0N q1N ⎟
.. ⎟ ⎠. . ··· · · · qM N
The coefficients c j ; j = 0, 1, . . . , N can be evaluated exactly by
(3.6)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.6 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
6
1 c j = (α ,β) h L, j
L
(α ,β)
wL
(α ,β)
j = 0, 1, . . . , N ,
(x)c (x) P L , j (x) dx,
(3.7)
0
for a general function c (x) one is unable to compute exactly the previous integral. So, the coefficients c j can be approximated by using the discrete inner product at the Jacobi–Gauss nodes to obtain N 1 (α ,β) (α ,β) (α ,β) (α ,β) c j = (α ,β) N ,i c x N ,i P L , j x N ,i , h L , j i =0
α ,β
where x N ,i , 0 i N are the zeros of Jacobi–Gauss quadrature in the interval (0, L ), with ing Christoffel numbers. Similarly, the coefficients q i j are given by
1 qi j = (α ,β) (α ,β) h τ ,i h L , j
τL
(α ,β)
q(x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
Nα,,βi , 0 i N are correspond-
(x) dx dt ,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N .
(3.8)
It is worthy to mention here that, since the pure spectral tau method would not use any interpolation operator, we approximated the solution and known functions in terms of shifted Jacobi polynomials in order to use orthogonality of such polynomials for computing the spectral tau coefficients. Now, making use of (2.14), (2.15) and (3.5), enables one to write
I tν c (x)
∂ 2 u (x, t ) CT Φ L , N (x) I tν ΦτT, M (t ) A D 2 Φ L , N (x) 2 ∂x = CT Φ L , N (x) ΦτT, M (t )PνT AD(2) Φ L , N (x) = ΦτT, M (t )PνT AD(2) Φ L , N (x)Φ LT, N (x)C.
(3.9)
Let
Φ L , N (x)Φ LT, N (x)C = HT Φ L , N (x),
(3.10)
where H is ( N + 1) × ( N + 1) matrix. Then (3.10) can be put in the form N
(α ,β)
ck P L ,k
(α ,β)
(x) P L , j (x) =
k =0
(α ,β)
k =0
(α ,β)
Hkj P L ,k
(x),
j = 0, 1 , . . . , N .
k =0
Multiplying both sides by P L ,i N
N
L ck
(α ,β)
P L ,i
(α ,β)
(α ,β)
(x) w L
(α ,β)
(x), i = 0, 1, . . . , N and integrating from 0 to L, we obtain (α ,β)
(x) P L ,k (x) P L , j (x) w L
(x) dx =
N k =0
0
L Hi j
(α ,β)
P L ,k
(α ,β)
(x) P L ,i
(α ,β)
(x) w L
(x) dx.
(3.11)
0
Making use of (2.5), yields
1
L N
Hi j = (α ,β) ck h L , i k =0
(α ,β) (α ,β) (α ,β) (α ,β) P L ,i (x) P L ,k (x) P L , j (x) w L (x) dx
,
i , j = 0, 1 , . . . , N .
(3.12)
0
Now, (3.9) and (3.10) give
I tν c (x)
∂ 2 u (x, t ) = ΦτT, M (t )PνT AD(2) HT Φ L , N (x). ∂ x2
(3.13)
Eq. (3.5) with (3.13) enable us to write the residual R N , M (x, t ) for Eq. (3.4) in the form
R M , N (x, t ) = ΦτT, M (t ) A − Q − PνT AD(2) H T Φ L , N (x)
= ΦτT, M (t )EΦ L , N (x),
(3.14)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.7 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
7
where
E = A − Q − PνT AD(2) H T . As in a typical tau method, see [3], we generate ( M + 1) × ( N − 1) linear algebraic equations in the unknown expansion coefficients, ai j , i = 0, 1, . . . , M; j = 0, 1, . . . , N − 2, namely
τ L
(α ,β)
R N , M (x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
(x) dx dt = 0,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N − 2 ,
(3.15)
and the rest of linear algebraic equations are obtained from the boundary conditions (3.3), as
ΦτT, M (t )AΦ L , N (0) = g 0 (t ), ΦτT, M (t )AΦ L , N ( L ) = g 1 (t ),
(3.16) α ,β
respectively. Eqs. (3.16) are collocated at ( M + 1) points. For suitable collocation points we use the shifted Jacobi roots t N ,i , (α ,β) i = 0, 1, . . . , M of P τ , M +1 (t ). The number of the unknown coefficients ai j is equal to ( N + 1) × ( M + 1) and can be obtained from Eqs. (3.15)–(3.16). Consequently u N , M (x, t ) given in Eq. (3.5) can be calculated. Remark 3.1. The cases α = β = 0 and α = β = −1/2 lead to the Legendre and Chebyshev tau methods combined with operational matrices respectively, which are special cases from Jacobi tau method and used most frequently in practice. It is actually more inclusive to study the whole family of Jacobi polynomials with the general parameters α , β > −1 rather than the Legendre or Chebyshev polynomials. 4. Fractional wave equation with damping In this section, we consider the time-fractional diffusion-wave equation with damping, namely
∂ u (x, t ) ∂ 2 u (x, t ) ∂ ν u (x, t ) +λ = + μs(x, t ), ν ∂t ∂t ∂ x2
0 < x L , 0 < t τ , 1 < ν 2,
(4.1)
with initial conditions (3.2) and boundary conditions (3.3). Reason for considering the fractional diffusion-wave equation with damping is that it has the ability of governing wave motion in a string and it arising from fractional diffusion-wave equation by assuming the damping effect due to the term λ ∂ u∂(xt ,t ) . In case of ν = 2, this equation can also be characterized as a telegraph equation which governs electrical transmission in a telegraph cable. After applying the Riemann–Liouville integral of order ν on Eq. (4.1) and making use of (2.4), we get the integrated form, namely ν −1
u (x, t ) − v (x, t ) + λ I t u (0, t ) = g 0 (t ),
∂ 2 u (x, t ) u (x, t ) = I tν ∂ x2
+ μ I tν s(x, t ),
u ( L , t ) = g 1 (t ),
(4.2)
ν −1
where v (x, t ) = q(x, t ) − λ I t f 0 (x). In order to use the tau method with the shifted Jacobi operational matrix for solving the fully integrated problem (4.2), we approximate u (x, t ) and q(x, t ) by the shifted Jacobi polynomials as in Eqs. (3.5) and (3.6) but s(x, t ) will be approximated in the form
s M , N (x, t ) = ΦτT, M (t )SΦ L , N (x),
(4.3)
where S is a known ( M + 1) × ( N + 1) matrix that can be written as
⎛
s00 ⎜ s10
S=⎜ ⎝ ..
s01 s11
.
.. .
s M0
s M1
··· ···
⎞
s0N s1N ⎟
.. ⎟ ⎠, . ··· · · · sM N
where si j are obtained by
(4.4)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.8 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
8
1 si j = (α ,β) (α ,β) h τ ,i h L , j
τL
(α ,β)
s(x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
(x) dx dt ,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N .
(4.5)
Using Eqs. (2.14), (2.15), (3.5) and (4.2), then it is easy to write
I tν −1 u (x, t ) = φτT , M (t )PνT −1 Aφ L , N (x),
I tν
∂ 2 u (x, t ) = φτT, M (t )PνT AD(2) φ L , N (x), ∂ x2
I tν s(x, t ) = φτT , M (t )PνT Sφ L , N (x).
(4.6)
Eqs. (3.5), (4.3) with (4.6), enable us to write the residual R N , M (x, t ) for Eq. (4.2) in the form
R N , M (x, t ) = φτT , M (t ) A − 2V + λPνT −1 A − PνT AD(2) − μPνT S φ L , N (x) = φτT , M (t )Eφ L , N (x),
(4.7)
where
E = A − 2V + λPνT −1 A − PνT AD(2) − μPνT S. We generate ( M + 1) × ( N − 1) linear algebraic equations in the unknown expansion coefficients, ai j , i = 0, 1, . . . , M; j = 0, 1, . . . , N − 2, see (3.15). As in Section 3, the rest of this linear algebraic equations are obtained from the boundary conditions (3.3), see (3.16). Then the approximate solution u N , M (x, t ) can be calculated by solving the generated system of ( M + 1) × ( N + 1) linear algebraic equations. 5. Fourth-order diffusion-wave equation In this section, we consider the fourth-order fractional diffusion-wave equation of the form:
∂ ν u (x, t ) ∂ 4 u (x, t ) + = s(x, t ), ∂tν ∂ x4
0 < x L , 0 < t τ , 1 < ν 2,
(5.1)
with initial conditions (3.2) and boundary conditions
u (0, t ) = g 0 (t ),
u ( L , t ) = g 1 (t ),
u xx (0, t ) = g 2 (t ),
u xx ( L , t ) = g 3 (t ).
(5.2)
We apply the same technique described in Section 3, the residual R N , M (x, t ) for Eq. (5.1) can be written as
R N , M (x, t ) = φτT , M (t ) A − Q + PνT AD(4) − PνT S φ L , N (x) = φτT , M (t )Eφ L , N (x),
(5.3)
where
E = A − Q + PνT AD(4) − PνT S. It is clear that, ( M + 1)( N − 3) linear algebraic equations in the unknown expansion coefficients, ai j , i = 0, 1, . . . , M; j = 0, 1, . . . , N − 4, are generated, and the rest of linear algebraic equations are obtained from the boundary conditions (5.2), as
φτT, M (t )Aφ L , N (0) = g 0 (t i ), φτT, M (t )Aφ L , N ( L ) = g 1 (t i ), φτT, M (t )AD(2) φ L , N (0) = g 2 (t i ), φτT, M (t )AD(2) φ L , N ( L ) = g 3 (t i ).
(5.4)
After collocating Eqs. (5.4) at ( M + 1) points with shifted Jacobi roots, ( M + 1) × ( N + 1) linear algebraic equations will be generated using Eqs. (3.15) and (5.4). Consequently u N , M (x, t ) given in Eq. (3.5) can be calculated.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.9 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
Table 1 MAEs at t = 0.8 with N = M = 8 and various choices of
9
ν , α and β for Example 1.
ν
α=β =0
α=β =1
α = β = 1/2
α = β = −1/2
9/8 10/8 11/8 12/8 13/8 14/8 15/8 16/8
3.55 · 10−7 5.29 · 10−7 5.44 · 10−7 4.45 · 10−7 3.08 · 10−7 1.96 · 10−7 1.06 · 10−7 3.23 · 10−12
2.86 · 10−7 4.19 · 10−7 4.17 · 10−7 3.39 · 10−7 2.80 · 10−7 2.15 · 10−7 7.19 · 10−8 5.36 · 10−12
3.39 · 10−7 4.95 · 10−7 4.94 · 10−7 3.94 · 10−7 2.96 · 10−7 2.13 · 10−7 1.20 · 10−7 2.28 · 10−12
3.24 · 10−7 4.23 · 10−7 1.07 · 10−7 4.61 · 10−7 3.17 · 10−7 1.83 · 10−7 1.07 · 10−7 1.14 · 10−11
Table 2 Absolute errors at x = 1 with N = M = 15,
α = β = 1 and various choices of ν for Example 1.
t
ν = 5/4
ν = 6/4
ν = 7/4
ν = 8/4
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
3.36 · 10−8 3.65 · 10−9 5.31 · 10−9 5.67 · 10−9 4.85 · 10−9 3.94 · 10−9 2.50 · 10−9 1.09 · 10−9 5.37 · 10−9
1.69 · 10−8 1.86 · 10−9 2.73 · 10−9 2.93 · 10−9 2.51 · 10−9 2.05 · 10−9 1.30 · 10−9 5.68 · 10−10 2.80 · 10−9
4.27 · 10−9 4.79 · 10−10 7.07 · 10−10 7.61 · 10−10 6.55 · 10−10 5.35 · 10−10 3.40 · 10−10 1.48 · 10−10 7.34 · 10−10
1.48 · 10−15 1.78 · 10−15 1.60 · 10−15 1.43 · 10−15 1.67 · 10−15 1.67 · 10−15 1.52 · 10−15 1.46 · 10−15 9.61 · 10−16
6. Numerical results To demonstrate the effectiveness of the proposed method in the present paper, five test examples are carried out in this section. Example 1. We consider the following problem [19]
∂ ν u (x, t ) 1 2 ∂ 2 u (x, t ) = x , ∂tν 2 ∂ x2
0 < x 1, 0 < t 1, 1 < ν 2,
(6.1)
with the initial conditions u (x, 0) = x, ut (x, 0) = x2 and boundary conditions u (0, t ) = 0, u (1, t ) = sinh(t ), while the exact solution of this initial value problem for ν = 2, is u (x, t ) = x + x2 sinh(t ), while the analytical solution using the VIM in [29] is given by
u (x, t ) = x + x2
∞ k =0
t k ν +1
Γ (kν + 2)
.
The maximum absolute errors (MAEs) between the exact solution u (x, 0.8) and the approximate solution u 8,8 (x, 0.8) with various choices of Jacobi parameters α , β and four choices of the fractional derivative ν are given in Table 1. In addition, Table 2 also shows the absolute error function |u (1, t ) − u 15,15 (1, t )| with α = β = 1 and various choices of ν , while Fig. 1 plots the error function u (x, t ) − u 10,10 (x, t ) with ν = 2 and α = β = 0. Tables 1 and 2 show the good agreement with the exact solutions obtained by the VIM [19] for different choices of x, t and ν . Example 2. Consider the one-dimensional time-fractional wave equation [30]
∂ 3/2 u (x, t ) ∂ 2 u (x, t ) = , ∂ t 3/2 ∂ x2
0 < x 1, 0 < t 1,
(6.2)
with initial conditions
u (x, 0) = sin(x),
∂ u (x, 0) = − sin(x), ∂t
and boundary conditions
u (0, t ) = 0,
u (1, t ) = sin(1) E 3/2,1 −t 3/2 − t E 3/2,2 −t 3/2 .
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.10 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
10
Fig. 1. The error at
Table 3 Absolute errors with N = M = 10,
ν = 2 with N = M = 10 and α = β = 0 for Example 1.
α = β = 0 and various choices of x and t for Example 2.
t
x = 0.1
x = 0.3
x = 0.5
x = 0.7
x = 0.9
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.19 · 10−6 5.88 · 10−7 1.14 · 10−6 8.97 · 10−7 6.51 · 10−7 2.24 · 10−7 8.43 · 10−7 9.53 · 10−7 1.02 · 10−6
3.48 · 10−6 1.59 · 10−6 3.42 · 10−6 2.82 · 10−6 1.77 · 10−6 7.04 · 10−7 2.29 · 10−6 2.66 · 10−6 3.39 · 10−6
6.05 · 10−6 2.38 · 10−6 5.51 · 10−6 5.04 · 10−6 2.55 · 10−6 1.57 · 10−6 3.64 · 10−6 3.57 · 10−6 4.25 · 10−6
8.56 · 10−6 3.22 · 10−6 7.07 · 10−6 7.62 · 10−6 3.90 · 10−6 2.56 · 10−6 6.39 · 10−6 6.14 · 10−6 5.36 · 10−6
6.17 · 10−6 1.15 · 10−6 5.00 · 10−6 8.32 · 10−6 7.61 · 10−6 2.21 · 10−6 4.48 · 10−6 9.57 · 10−6 1.37 · 10−5
Fig. 2. The error function at N = M = 15 and
α = β = −0.5 for Example 2.
Momani et al. [30] applied the generalized two-dimensional differential transform to both sides of problem (6.2) and showed that its exact solution is E 3/2,1 (−t 3/2 ) sin(x) − t E 3/2,2 (−t 3/2 ) sin(x), where
E ζ,η ( z) =
∞ k =0
zk
Γ (ζ k + η)
.
We solved this problem for different values of N, M and ν and compared our results with the exact solution in Table 3. Fig. 2 shows the error function u (x, t ) − u 15,15 (x, t ) with α = β = −0.5, while Fig. 3 shows the error functions u (x, 0.6) − u 10,10 (x, 0.6) and u (0.1, t ) − u 10,10 (0.1, t ) with α = β = 0.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.11 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
Fig. 3. Errors at t = 0.6 (left) and x = 0.1 (right) with N = M = 10, Table 4 MAEs at t = 1, L = 2,
11
α = β = 0 for Example 2.
τ = 1 with α = β = 0 and various choices of N, M and ν for Example 3.
N=M
ν = 1.1
ν = 1.5
ν = 1.7
ν = 1.9
3 5 7 9 11 13 15 17 19 21 23
1.74 · 10−4 3.76 · 10−5 1.30 · 10−5 5.73 · 10−6 2.90 · 10−6 1.63 · 10−6 9.94 · 10−7 6.26 · 10−7 4.30 · 10−7 3.00 · 10−7 2.16 · 10−7
6.17 · 10−4 1.44 · 10−4 5.84 · 10−5 2.76 · 10−5 1.51 · 10−5 9.03 · 10−6 5.77 · 10−6 3.88 · 10−6 2.72 · 10−6 1.97 · 10−6 1.46 · 10−6
2.72 · 10−4 6.98 · 10−5 3.15 · 10−5 1.55 · 10−5 9.15 · 10−6 5.70 · 10−6 3.80 · 10−6 2.64 · 10−6 1.91 · 10−6 1.41 · 10−6 9.15 · 10−7
6.44 · 10−5 8.57 · 10−6 1.92 · 10−6 2.96 · 10−7 2.29 · 10−7 2.73 · 10−7 2.66 · 10−7 2.26 · 10−7 1.91 · 10−7 1.59 · 10−7 1.32 · 10−7
Example 3. We consider the following fractional wave equation with damping
∂ ν u (x, t ) ∂ u (x, t ) ∂ 2 u (x, t ) + = + q(x, t ), ∂tν ∂t ∂ x2
0 < x L , 0 < t τ , 1 < ν 2,
(6.3)
with initial conditions
u (x, 0) = 0,
∂ u (x, 0) = 0, ∂t
and boundary equations
u (0, t ) = 0,
u ( L , t ) = 0,
and
q(x, t ) =
2x( L − x) 2−ν t + 2tx( L − x) + 2t 2 . Γ (3 − ν )
The exact solution of this problem is u (x, t ) = t 2 x( L − x). In [5], Chen et al. applied the method of separation of variables with constructing the implicit difference approximation to introduce an approximate solution of this problem at L = 2, T = 1 and ν = 1.7 using the Caputo form of the fractional derivative. Regarding problem (6.3), in [5], the best result is achieved at L = 2, T = 1 and ν = 1.7 with 640 steps and the maximum absolute error is 5.31 × 10−5 . In Table 4, we list MAEs at t = 1, L = 2, τ = 1 with α = β = 0 and various choices of N, M and ν . Our method is more accurate than the method of separation of variables combined with the implicit difference approximation [5]. We see in this table that the results are accurate for even small choices of N and M. In addition, to demonstrate the convergence of the proposed method, in Fig. 4, we plot the logarithmic graphs of MAEs (log10 Error) at τ = 1, ν = 1.11, 1.55, 1.77, 1.88 and α = β = 1 with various values of N (N = M); by using the presented algorithm. Clearly, the numerical errors decay rapidly as N and M increase. In order to confirm the accuracy and applicability of the proposed method for all the Jacobi parameters α , β , we introduce in Table 5, the MAEs at t = 1, L = 1, τ = 1 with N = M = 8 and various choices of α , β , and four different choices of ν . From the numerical results of this table, it is observed that neither Chebyshev polynomial (α = β = − 12 ) nor Legendre polynomial (α = β = 0) is the best special case of Jacobi polynomials. In general, it can be concluded that, only a limited number of shifted Jacobi polynomials is needed to obtain a satisfactory result. This numerical experiment demonstrates the utility of the method.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.12 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
12
Fig. 4. Convergence of problem (6.3) at
Table 5 MAEs at t = 1, L = 1,
α 0 1 1/2 −1/2 −1/2
τ = 1, ν = 1.11, 1.55, 1.77, 1.88.
τ = 1 with N = M = 8 and various choices of α , β and ν for Example 3.
β
ν = 1.2
ν = 1.4
ν = 1.6
ν = 1.8
0 1 1/2 −1/2 1/2
5.35 · 10−6
1.01 · 10−5
8.77 · 10−6
2.27 · 10−6 9.75 · 10−6 3.88 · 10−6 7.11 · 10−6 7.36 · 10−6
5.34 · 10−6 5.59 · 10−6 4.00 · 10−6 3.93 · 10−7
8.67 · 10−6 9.87 · 10−6 7.98 · 10−6 1.40 · 10−7
2.01 · 10−6 5.67 · 10−6 9.15 · 10−6 3.00 · 10−6
Example 4. Consider the following two-term wave-diffusion equation with damping with 1 < ν < 2:
∂ ν u (x, t ) ∂ u (x, t ) ∂ 2 u (x, t ) + = + q(x, t ), ∂tν ∂t ∂ x2
0 < x 1, 0 < t τ , 1 < ν < 2,
(6.4)
with initial conditions
u (x, 0) = 0,
∂ u (x, 0) = 0, ∂t
and boundary equations
u (0, t ) = t 3 ,
u (1, t ) = et 3 ,
and
q(x, t ) =
6t 3−ν
Γ (4 − ν )
e x + 3t 2 e x − t 3 e x .
The exact solution of this problem is u (x, t ) = e x t 3 . Recently, Liu et al. [23] proposed the fractional predictor–corrector method (FPCM [23]) to solve this problem with several choices of h and τ = h2 , where h and τ are space and time step sizes, respectively. We make a comparison of the presented method at three choices of Jacobi parameters (α , β) with the numerical method proposed in FPCM [23], where the implicit numerical schemes are implemented for space and time discretizations. In case of ν = 1.85, numerical maximum absolute errors with three different M are listed in Table 6. Meanwhile, the results of FPCM [23] are presented in the sixth column of this table. Obviously, our method for all choices of Jacobi parameters is more accurate than FPCM [23]. Moreover, we can see that for this problem the special choice α = β = 0 is the best one. Accordingly, we choose the special case α = β = 0 to plot in Fig. 5, the logarithmic graphs of MAEs (log10 Error) at τ = 1, ν = 1.09, 1.5, 1.85, 1.95, 1.999 with various values of N (N = M); by using the presented algorithm. From this figure, we conclude that the numerical solution is more accurate when ν tends to its integer values and the numerical errors for all chosen ν decay rapidly as N and M increase. We observe also that the suggested algorithm provides accurate and stable numerical results. Example 5. Consider the following fourth-order fractional diffusion-wave equation
∂ ν u (x, t ) ∂ 4 u (x, t ) x Γ (ν + 4) 3 ν +3 + t , + = e + t t ∂tν 6 ∂ x4
0 < x 1, 0 < t 1, 1 < ν 2,
(6.5)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.13 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
Table 6 Comparing of our scheme with FPCM [23] at
13
ν = 1.85 for Example 4.
N=M
α=β =0
α = β = 0.5
α=β =1
h
FPCM [23]
4 8 16
5.46 · 10−5 5.50 · 10−6 2.42 · 10−7
1.16 · 10−4 1.15 · 10−5 9.75 · 10−7
2.48 · 10−4 2.64 · 10−5 1.25 · 10−6
1/4 1/8 1/16
1.09 · 10−1 2.76 · 10−2 6.72 · 10−3
Fig. 5. Convergence of problem (6.4) at Table 7 Absolute errors at t = 1 with N = M = 9,
τ = 1, ν = 1.09, 1.5, 1.85, 1.95, 1.999.
α = β = 0 and various choices of ν for Example 5.
x
ν = 1.2
ν = 1.3
ν = 1.5
ν = 1.7
ν = 1.8
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.53 · 10−7 2.26 · 10−7 2.88 · 10−7 3.36 · 10−7 3.65 · 10−7 3.75 · 10−7 3.62 · 10−7 3.26 · 10−7 2.69 · 10−7
1.88 · 10−7 2.77 · 10−7 3.54 · 10−7 4.14 · 10−7 4.55 · 10−7 4.72 · 10−7 4.61 · 10−7 4.18 · 10−7 3.42 · 10−7
1.49 · 10−7 2.01 · 10−7 2.50 · 10−7 2.99 · 10−7 3.50 · 10−7 3.95 · 10−7 4.20 · 10−7 4.08 · 10−7 3.47 · 10−7
7.20 · 10−8 7.67 · 10−8 8.47 · 10−8 1.09 · 10−7 1.57 · 10−7 2.20 · 10−7 2.78 · 10−7 2.99 · 10−7 2.61 · 10−7
4.29 · 10−8 3.76 · 10−8 3.44 · 10−8 4.77 · 10−8 8.56 · 10−8 1.43 · 10−7 1.99 · 10−7 2.26 · 10−7 1.95 · 10−7
with initial conditions
u (x, 0) = 0,
∂ u (x, 0) = ex, ∂t
and boundary conditions
u (0, t ) = t ν +3 + t , u xx (0, t ) = t ν +3 + t ,
u (1, t ) = e t ν +3 + t ,
u xx (1, t ) = e t ν +3 + t .
The exact solution of this problem is e x (t ν +3 + t ). Table 7 exhibits the absolute error function |u (x, 1) − u 9,9 (x, 1)|, by applying the shifted Jacobi tau method based on the operational matrices at α = β = 0 and five choices of ν . We present comparisons of our method at α = β = 0 with the numerical method proposed in [18], where the celebrated Crank–Nicolson finite difference (CCNFDM [18]) method is implemented to solve problem (6.5) with different h and τ , where h and τ are space and time step sizes, respectively. In Tables 8 and 9, we introduce maximum absolute errors using the proposed method for different choices of N and those achieved by CCNFDM [18] for various choices of the fractional derivative ν . Again and from numerical results presented in these tables, we clearly observe that the proposed global scheme is accurate for even small values of N and our scheme is more accurate than CCNFDM [18]. 7. Conclusion In this article, the fractional derivatives are described in the Caputo sense to allow traditional initial and boundary conditions to be included in the formulation of the problem. In addition, the Jacobi polynomials are used to approximate the second- and fourth-order fractional diffusion-wave equation and fractional diffusion wave equation with damping and
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.14 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
14
Table 8 Comparing of our scheme with CCNFDM [18] at
ν = 1.3, 1.5, 1.8 for Example 5.
Our scheme (α = β = 0)
CCNFDM [18] with h = 1/500
N=M
ν = 1.3
ν = 1.5
ν = 1.8
τ
ν = 1.3
ν = 1.5
ν = 1.8
4 6 8 10 12
1.80 · 10−3 1.92 · 10−5 1.34 · 10−6 1.59 · 10−7 3.33 · 10−8
3.55 · 10−3 2.41 · 10−5 1.33 · 10−6 1.57 · 10−7 3.36 · 10−8
7.41 · 10−3 1.36 · 10−5 6.57 · 10−7 9.05 · 10−8 1.60 · 10−8
1/5 1/10 1/20 1/40 1/80
7.21 · 10−3 2.25 · 10−3 6.96 · 10−4 2.15 · 10−4 6.63 · 10−5
1.72 · 10−2 6.21 · 10−3 2.21 · 10−3 7.88 · 10−4 2.79 · 10−4
5.73 · 10−2 2.57 · 10−2 1.13 · 10−2 4.99 · 10−3 2.18 · 10−3
Table 9 Comparing of our scheme with CCNFDM [18] at
ν = 1.2, 1.5, 1.7 for Example 5.
Our scheme (α = β = 0) N=M
ν = 1.2
ν = 1.5
ν = 1.7
h
τ = 1/20 000 ν = 1.2
ν = 1.5
ν = 1.7
4 6 8 10 12
1.15 · 10−3 1.48 · 10−5 1.00 · 10−6 1.54 · 10−7 3.18 · 10−8
3.55 · 10−3 2.41 · 10−5 1.33 · 10−6 1.57 · 10−7 3.36 · 10−8
5.97 · 10−3 1.87 · 10−5 8.85 · 10−7 1.09 · 10−7 3.36 · 10−8
1/5 1/10 1/20 1/40 1/80
1.20 · 10−3 3.03 · 10−4 7.66 · 10−5 1.91 · 10−5 4.79 · 10−6
1.18 · 10−3 2.98 · 10−4 7.55 · 10−5 1.89 · 10−5 4.78 · 10−6
1.16 · 10−4 2.94 · 10−4 7.50 · 10−5 1.94 · 10−6 4.90 · 10−6
CCNFDM [18] with
reduce them to systems of linear algebraic equations in the unknown expansion coefficients of the sought-for spectral solutions. The reduction is based on the operational matrix formulation of shifted Jacobi tau method. This proposed method is quite different essentially from those methods presented in literature as it discretized the underline problem spectrally in space and time discretizations. The numerical experiments demonstrated the good accuracy of the described method. Moreover, only a small number of shifted Jacobi polynomials is needed to obtain a satisfactory result. The algorithms presented in this paper can be easily extended to a large amount of space–time fractional differential equations with different boundary conditions. Acknowledgements The authors are very grateful to the referees for carefully reading of the paper and for their comments and suggestions which have improved the paper. References [1] R.T. Baillie, Long memory processes and fractional integration in econometrics, J. Econom. 73 (1996) 5–59. [2] A.H. Bhrawy, M.M. Tharwat, M.A. Alghamdi, A new operational matrix of fractional integration for shifted Jacobi polynomials, Bull. Malays. Math. Soc. (2014), in press. [3] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, New York, 1989. [4] R.C. Charo, C.R. Luz, R. Ikehata, Sharp decay rates for wave equations with a fractional damping via new method in the Fourier space, J. Math. Anal. Appl. 408 (2013) 247–255. [5] J. Chen, F. Liu, V. Anh, S. Shen, Q. Liu, C. Liao, The analytical solution and numerical solution of the fractional diffusion-wave equation with damping, Appl. Math. Comput. 219 (2012) 1737–1748. [6] M.R. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804. [7] M.R. Cui, Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation, Numer. Algorithms 62 (2013) 383–409. [8] R. Darzi, B. Mohammadzade, S. Mousavi, R. Beheshti, Sumudu transform method for solving fractional differential equations and fractional diffusionwave equation, J. Math. Comput. Sci. 6 (2013) 79–84. [9] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order, Comput. Math. Appl. 62 (2011) 2364–2373. [10] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A new Jacobi operational matrix: an application for solving fractional differential equations, Appl. Math. Model. 36 (2012) 4931–4943. [11] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, Numerical approximations for fractional diffusion equations via a Chebyshev spectral-tau method, Cent. Eur. J. Phys. 11 (2013) 1494–1503. [12] R. Du, W.R. Cao, Z.Z. Sun, A compact difference scheme for the fractional diffusion-wave equation, Appl. Math. Model. 34 (2010) 2998–3007. [13] G.H. Gao, Z.Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [14] C.F.L. Godinho, J. Weberszpil, J.A. Helayel-Neto, Extending the D’Alembert solution to space–time modified Riemann–Liouville fractional wave equations, Chaos Solitons Fractals 45 (2012) 765–771. [15] J.H. He, Nonlinear oscillation with fractional derivative and its applications, in: International Conference on Vibrating Engineering, Dalian, China, 1998, pp. 288–291. [16] J.H. He, Some applications of nonlinear fractional differential equations and their applications, Bull. Sci. Technol. 15 (2) (1999) 86–90. [17] X. Hu, L. Zhang, A compact finite difference scheme for the fourth-order fractional diffusion-wave system, Comput. Phys. Commun. 182 (2011) 1645–1650. [18] X. Hu, L. Zhang, On finite difference methods for fourth-order fractional diffusion-wave and subdiffusion systems, Appl. Math. Comput. 218 (2012) 5019–5034. [19] H. Jafari, S. Momani, Solving fractional diffusion and wave equations by modified homotopy perturbation method, Phys. Lett. A 370 (2007) 388–396.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.15 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
15
[20] H. Jiang, F. Liu, I. Turner, K. Burrage, Analytical solutions for the multi-term time-fractional diffusion-wave/diffusion equations in a finite domain, Comput. Math. Appl. 64 (2012) 3377–3388. [21] C.P. Li, A. Chen, J.J. Ye, Numerical approaches to fractional calculus and fractional ordinary differential equation, J. Comput. Phys. 230 (2011) 3352–3368. [22] C.P. Li, Z.G. Zhao, Y.Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (2011) 855–875. [23] F. Liu, M.M. Meerschaert, R.J. McGough, P. Zhuang, Q. Liu, Numerical methods for solving the multi-term time-fractional wave-diffusion equation, Fract. Calc. Appl. Anal. 16 (2013) 9–25. [24] Y. Luke, The Special Functions and Their Approximations, vol. 2, Academic Press, New York, 1969. [25] R.L. Magin, Fractional calculus in bioengineering, Crit. Rev. Biomed. Eng. 32 (1) (2004) 1–104. [26] F. Mainardi, Fractional calculus: ‘Some basic problems in continuum and statistical mechanics’, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer-Verlag, New York, 1997, pp. 291–348. [27] B. Mandelbrot, Some noises with 1/f spectrum, a bridge between direct current and white noise, IEEE Trans. Inf. Theory 13 (2) (1967) 289–298. [28] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37 (2004) 161–208. [29] S. Momani, Z. Odibat, Comparison between the homotopy perturbation method and the variational iteration method for linear fractional partial differential equations, Comput. Math. Appl. 54 (2007) 910–919. [30] S. Momani, Z. Odibat, V.S. Erturk, Generalized differential transform method for solving a space- and time-fractional diffusion-wave equation, Phys. Lett. A 370 (2007) 379–387. [31] R.R. Nigmatullin, To the theoretical explanation of the universal response, Phys. Status Solidi, B Basic Res. 123 (2) (1984) 739–745. [32] R.R. Nigmatullin, Realization of the generalized transfer equation in a medium with fractal geometry, Phys. Status Solidi, B Basic Res. 133 (1) (1986) 425–430. [33] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974. [34] R. Ren, H. Li, W. Jiang, M. Song, An efficient Chebyshev-tau method for solving the space fractional diffusion equations, Appl. Math. Comput. 224 (2013) 259–267. [35] J. Ren, Z.-Z. Sun, X. Zhao, Compact difference scheme for the fractional sub-diffusion equation with Neumann boundary conditions, J. Comput. Phys. 232 (2013) 456–467. [36] Y.A. Rossikhin, M.V. Shitikova, Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids, Appl. Mech. Rev. 50 (1997) 15–67. [37] A. Saadatmandi, M. Dehghan, A tau approach for solution of the space fractional diffusion equation, Comput. Math. Appl. 62 (2011) 1135–1142. [38] Q. Xu, J.S. Hesthaven, Stable multi-domain spectral penalty methods for fractional partial differential equations, J. Comput. Phys. 257 (2014) 241–258. [39] M. Zayernouri, G.E. Karniadakis, Exponentially accurate spectral and spectral element methods for fractional ODEs, J. Comput. Phys. 257 (2014) 460–480. [40] F. Zeng, C. Li, F. Liu, I. Turner, The use of finite difference/element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput. 35 (6) (2013) A2976–A3000.
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.1 (1-15)
Journal of Computational Physics ••• (••••) •••–•••
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A spectral tau algorithm based on Jacobi operational matrix for numerical solution of time fractional diffusion-wave equations A.H. Bhrawy a,b , E.H. Doha c , D. Baleanu d,e,f,∗ , S.S. Ezz-Eldien g a
Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt d Department of Chemical and Materials Engineering, Faculty of Engineering, King Abdulaziz University, Jeddah, Saudi Arabia e Department of Mathematics and Computer Sciences, Cankaya University, Ankara, Turkey f Institute of Space Sciences, Magurele-Bucharest, R 76900, Romania g Department of Basic Science, Institute of Information Technology, Modern Academy, Cairo, Egypt b c
a r t i c l e
i n f o
Article history: Received 11 January 2014 Received in revised form 7 March 2014 Accepted 23 March 2014 Available online xxxx Keywords: Fractional diffusion-wave equations Tau method Shifted Jacobi polynomials Operational matrix Caputo derivative
a b s t r a c t In this paper, an efficient and accurate spectral numerical method is presented for solving second-, fourth-order fractional diffusion-wave equations and fractional wave equations with damping. The proposed method is based on Jacobi tau spectral procedure together with the Jacobi operational matrix for fractional integrals, described in the Riemann– Liouville sense. The main characteristic behind this approach is to reduce such problems to those of solving systems of algebraic equations in the unknown expansion coefficients of the sought-for spectral approximations. The validity and effectiveness of the method are demonstrated by solving five numerical examples. Numerical examples are presented in the form of tables and graphs to make comparisons with the results obtained by other methods and with the exact solutions more easier. © 2014 Elsevier Inc. All rights reserved.
1. Introduction Due to the accuracy of ordinary and partial fractional differential equations in modeling various engineering and physical phenomena such as bioengineering [25], solid mechanics [36], anomalous transport [28], continuum and statistical mechanics [26], fluid-dynamic [16], economics [1], nonlinear oscillation of earthquakes [15], colored noise [27] and many other problems [33], fractional calculus has become the focus of many researchers in recent years. The time-fractional diffusion-wave equation, which is a mathematical model of a wide class of important physical phenomena, is a linear integro-partial fractional differential equation that obtained from the classical diffusion-wave equation by replacing the second-order time derivative term by a fractional derivative of order ν , 1 < ν 2. Many of the universal mechanical, acoustic and electromagnetic responses may be described accurately by the time-fractional diffusion-wave equation, see [31,32]. The fourth-order space derivative is arising in the wave propagation in beams and modeling formation of grooves on a flat surface, thus considerable attention has been devoted to fourth-order fractional diffusion-wave equation and its applications, see [33].
*
Corresponding author.
http://dx.doi.org/10.1016/j.jcp.2014.03.039 0021-9991/© 2014 Elsevier Inc. All rights reserved.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.2 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
2
Several numerical methods have been proposed in the last few years for solving such equations. In [6,7], Cui presented high-order compact finite difference methods to solve the fractional diffusion and time fractional diffusion equation. A compact difference scheme [17,12] is proposed for solving second- and fourth-order fractional diffusion-wave equations. Chen et al. [5] proposed and developed the method of separation of variables with constructing the implicit difference scheme for fractional diffusion-wave equation with damping. Recently, Hu and Zhang [18] proposed the finite difference scheme for approximating the solution of the fourth-order fractional diffusion-wave equation. An explicit analytical solution is proposed in [14] for space–time fractional wave equation. In [8], Darzi et al. solved fractional diffusion-wave equations, by using the Sumudu transform method. Some other numerical methods are introduced for solving fractional differential equations with sub-diffusion and super-diffusion, see [13,21,22,4,20,35,40]. Spectral methods are a class of important tools for obtaining the numerical solutions of fractional differential equations. They have excellent error properties and they offer exponential rates of convergence for smooth problems [3,38,39]. Considerable attention has been devoted to the use of orthogonal polynomials aiming to deal with various fractional differential equations, by reducing them to systems of algebraic equations that greatly simplify obtaining their solutions. Moreover, the operational matrices of fractional derivatives are derived and used together with spectral methods based on orthogonal polynomials for developing numerical approximations to fractional differential equations, see [9,10]. Recently, and with the help of operational matrices of fractional derivatives for orthogonal polynomials, the tau spectral method is also utilized to solve spectrally space-fractional diffusion equations, see [11,34,37]. In the current work, we investigate and develop some efficient spectral algorithms to solve spectrally second- and fourthorder fractional diffusion-wave equations and fractional diffusion-wave equation with damping. These algorithms are based principally on shifted Jacobi tau-spectral method combined with the generalized shifted Jacobi operational matrix of integrals. The current article is organized as follows. In Section 2, the shifted Jacobi polynomials with their properties and their operational matrix to fractional integration are presented. In Sections 3, 4 and 5, the shifted Jacobi operational matrix of fractional integration is applied together with the shifted Jacobi tau spectral method to give numerical solutions for secondand fourth-order fractional diffusion-wave equations and fractional diffusion-wave equation with damping, respectively. In Section 6, several numerical examples and comparisons between our obtained numerical results and those of other methods are considered. Also a conclusion is given in Section 7. 2. Operational matrices of shifted Jacobi polynomials Riemann–Liouville and Caputo fractional definitions are the two most used from all the other definitions of fractional calculus which have been introduced recently. Definition 2.1. The integral of order
I ν f (x) =
x
1
Γ (ν )
ν 0 (fractional) according to Riemann–Liouville is given by
(x − t )ν −1 f (t ) dt ,
ν > 0, x > 0,
0
I 0 f (x) = f (x),
(2.1)
where
∞ Γ (ν ) =
xν −1 e −x dx
0
is gamma function. The operator J ν satisfies the following properties
I ν I μ f (x) = I ν +μ f (x), I ν I μ f (x) = I μ I ν f (x), Γ (β + 1) β+ν I ν xβ = x . Γ (β + 1 + ν )
(2.2)
Definition 2.2. The Caputo fractional derivative of order
D ν f (x) =
1
Γ (m − ν )
x
(x − t )m−ν −1
dm dt m
f (t ) dt ,
0
where m is the ceiling function of ν . The operator D ν satisfies the following properties
ν is defined by m − 1 < ν m , x > 0,
(2.3)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.3 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
3
D ν I ν f (x) = f (x), I ν D ν f (x) = f (x) −
m −1
f ( i ) 0+
xi i!
i =0
D ν xβ =
,
Γ (β + 1) β−ν x . Γ (β + 1 − ν )
(2.4) (α ,β)
The Jacobi polynomials, denoted by P j
(z); α > −1, β > −1, constitute an orthogonal system with respect to the = [−1, 1], i.e.,
weight function w (α ,β) ( z) = (1 − z)α (1 + z)β over I
1
(α ,β)
Pj
(α ,β)
( z) P k
(α ,β)
( z) w (α ,β) ( z) dz = δ jk γk
,
−1
where δ jk is the Kronecker function and
γk(α ,β) =
2α +β+1 Γ (k + α + 1)Γ (k + β + 1)
(2k + α + β + 1)k!Γ (k + α + β + 1)
.
In order to use Jacobi polynomials on the interval x ∈ [0, L ] we define the so-called shifted Jacobi polynomials by in2x (α ,β) 2x (α ,β) − 1. Let the shifted Jacobi polynomials P j ( − 1) be denoted by P L , j (x). L L (α ,β) (x) = Then the shifted Jacobi polynomials are constituting an orthogonal system with respect to the weight function w L xβ ( L − x)α over I = [0, L ] with the orthogonality property
troducing the change of variable z =
L
(α ,β)
P L, j
(α ,β)
(α ,β)
(x) P L ,k (x) w L
(α ,β)
(x) dx = h L ,k ,
(2.5)
0
where (α ,β)
h L ,k
=
α +β+1 L
(α ,β)
δ jk γ j
2
=
L α +β+1 Γ (k + α + 1)Γ (k + β + 1)
(2k + α + β + 1)k!Γ (k + α + β + 1) (α ,β)
The explicit analytic form of the shifted Jacobi polynomials P L ,i (α ,β)
P L ,i
(x) =
i
(−1)i −k
k =0
δ jk .
(x) of degree i is given by
Γ (i + β + 1)Γ (i + k + α + β + 1) xk , Γ (k + β + 1)Γ (i + α + β + 1)(i − k)!k! Lk
(2.6)
and this in turn, enables one to get
Γ (i + β + 1)(i + α + β + 1)q Γ ( i + β + 1) (α ,β) D q P L ,i (0) = (−1)i −q q , , Γ (β + 1)i ! L Γ (i − q + 1)Γ (q + β + 1) Γ (i + α + 1)(i + α + β + 1)q Γ ( i + α + 1) (α ,β) (α ,β) P L ,i ( L ) = D q P L ,i ( L ) = q , , q i, Γ (α + 1)i ! L Γ (i − q + 1)Γ (q + α + 1) (α ,β)
P L ,i
(0) = (−1)i
q i,
where (a)q = Γ (a + q)/Γ (a), and q is an integer, these relations will be of important use later.
(α ,β)
Assume u (x) is a square integrable function with respect to the Jacobi weight function w L expressed in terms of shifted Jacobi polynomials as
u (x) =
∞
(α ,β)
c j P L, j
in (0, L ), then it can be
(x),
j =0
from which the coefficients c j are given by
1 c j = (α ,β) h L, j
L
(α ,β)
wL
(α ,β)
(x)u (x) P L , j
dx,
j = 0, 1 , . . . .
0
If we approximate u (x) by the first ( N + 1)-terms, then we can write
(2.7)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.4 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
4
N
u N (x) =
(α ,β)
c j P L, j
(x),
(2.8)
j =0
which alternatively may be written in the matrix form:
u N (x) = C T Φ L , N (x),
C T ≡ [c 0 , c 1 , . . . , c N ],
(2.9)
with
(α ,β) T (α ,β) (α ,β) Φ L , N (x) ≡ P L ,0 (x), P L ,1 (x), . . . , P L , N (x) .
(2.10)
Similarly, let u (x, t ) be an infinitely differentiable function defined on 0 < x L and 0 < t τ . Then it is possible to express as
u M , N (x, t ) =
N M
(α ,β)
a i j P τ ,i
(α ,β)
(t ) P L , j (x) = ΦτT, M (t )AΦ L , N (x),
(2.11)
i =0 j =0
with
⎛
a00 ⎜ a10
a01 a11
A=⎜ ⎝ ..
.
.. .
a M0
a M1
⎞
··· ···
a0N a1N ⎟
.. ⎟ ⎠, . ··· · · · aM N
and
1 ai j = (α ,β) (α ,β) h τ ,i h L , j
τL
(α ,β)
u (x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
(x) dx dt ,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N .
(2.12)
The first derivative of the vector Φ L , N (x) can be expressed by
dΦ L , N (x) dx
= DΦ L , N (x),
(2.13)
where D is the ( N + 1) × ( N + 1) operational matrix of derivative given by
D = (di j ) =
ξ(i , j ), i > j , 0
otherwise,
and
ξ(i , j ) =
L α +β (i + α + β + 1)(i + α + β + 2) j ( j + α + 2)i − j −1 Γ ( j + α + β + 1)
× 3 F2
(i − j − 1)!Γ (2 j + α + β + 1) − i + 1 + j , i + j + α + β + 2, j + α + 1 j + α + 2,
2j +α +β +2
;1 .
(For the general definition of a generalized hypergeometric series and special 3 F 2 , see [24], p. 41 and pp. 103–104, respectively.) Repeated use of (2.13) enables one to write
dq Φ L , N (x) dxq
= Dq Φ L , N (x),
(2.14)
where q is a natural number and Dq means matrix powers. The fractional integration of order ν of φτ , M (t ) can be expressed as
I ν Φτ , M (t ) = Pν Φτ , M (t ), where Pν is the ( M + 1) × ( M + 1) operational matrix of fractional integration of order and is defined by
(2.15)
ν in the Riemann–Liouville sense
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.5 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
⎛ Ω (0, 0, α , β) Ω (0, 1, α , β) ν ν ⎜ Ων (1, 0, α , β) Ων (1, 1, α , β) ⎜ .. .. ⎜ . . ⎜ Pν = ⎜ Ων (i , 1, α , β) ⎜ Ων (i , 0, α , β) ⎜ .. .. ⎝ . . Ων ( M , 0, α , β) Ων ( M , 1, α , β)
··· ··· ··· ··· ··· ···
Ων (0, M , α , β) Ων (1, M , α , β) .. .
5
⎞
⎟ ⎟ ⎟ ⎟ ⎟, Ων (i , M , α , β) ⎟ ⎟ .. ⎠ . Ων ( M , M , α , β)
where
Ων (i , j , α , β) =
i (−1)i −k Γ (i + β + 1)Γ (i + k + α + β + 1)Γ (α + 1)τ ν
Γ (k + β + 1)Γ (i + α + β + 1)(i − k)!Γ (k + ν + 1)
k =0
×
j (−1) j − f Γ ( j + f + α + β + 1)Γ ( f + k + β + ν + 1)(2 j + α + β + 1) j ! Γ ( j + α + 1)Γ ( f + β + 1)( j − f )! f !Γ ( f + k + α + β + ν + 2) f =0
(see [2] for proof). 3. Second-order diffusion-wave equation In this section, the spectral tau method together with the shifted Jacobi operational matrix of integrations are applied to solve second-order fractional diffusion-wave equations. Consider the fractional diffusion-wave equation of the form:
∂ 2 u (x, t ) ∂ ν u (x, t ) = c ( x ) , ∂tν ∂ x2
0 < x L , 0 < t τ , 1 < ν 2,
(3.1)
with initial conditions
u (x, 0) = f 0 (x),
ut (x, 0) = f 1 (x),
(3.2)
u ( L , t ) = g 1 (t ).
(3.3)
and boundary conditions
u (0, t ) = g 0 (t ),
If we apply the Riemann–Liouville integral of order Eq. (3.1), namely
ν on Eq. (3.1) and making use of (2.4), we get the integrated form of
∂ 2 u (x, t ) , ∂ x2 u ( L , t ) = g 1 (t ),
u (x, t ) − q(x, t ) = I tν c (x) u (0, t ) = g 0 (t ),
(3.4)
where q(x, t ) = t f 1 (x) + f 0 (x) and I tν is the fractional integration of order
ν with respect to t. In order to use the pure spectral tau method with the shifted Jacobi operational matrix for solving the fully integrated problem (3.4), we assume that c (x) and q(x, t ) are smooth and continuous functions. Therefore, u (x, t ), c (x) and q(x, t ) can be expanded by the shifted Jacobi polynomials as u M , N (x, t ) = ΦτT, M (t )AΦ L , N (x), c N (x) = C T Φ L , N (x), q M , N (x, t ) = ΦτT, M (t )QΦ L , N (x),
(3.5)
where A is the unknown coefficient ( M + 1) × ( N + 1) matrix, C T and Q are known matrices that can be written as
C T = [c 0 , c 1 , . . . , c N ], ⎛q q01 · · · 00 ⎜ q10 q11 · · · Q=⎜ . ⎝ .
..
..
q M0
q M1
⎞
q0N q1N ⎟
.. ⎟ ⎠. . ··· · · · qM N
The coefficients c j ; j = 0, 1, . . . , N can be evaluated exactly by
(3.6)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.6 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
6
1 c j = (α ,β) h L, j
L
(α ,β)
wL
(α ,β)
j = 0, 1, . . . , N ,
(x)c (x) P L , j (x) dx,
(3.7)
0
for a general function c (x) one is unable to compute exactly the previous integral. So, the coefficients c j can be approximated by using the discrete inner product at the Jacobi–Gauss nodes to obtain N 1 (α ,β) (α ,β) (α ,β) (α ,β) c j = (α ,β) N ,i c x N ,i P L , j x N ,i , h L , j i =0
α ,β
where x N ,i , 0 i N are the zeros of Jacobi–Gauss quadrature in the interval (0, L ), with ing Christoffel numbers. Similarly, the coefficients q i j are given by
1 qi j = (α ,β) (α ,β) h τ ,i h L , j
τL
(α ,β)
q(x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
Nα,,βi , 0 i N are correspond-
(x) dx dt ,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N .
(3.8)
It is worthy to mention here that, since the pure spectral tau method would not use any interpolation operator, we approximated the solution and known functions in terms of shifted Jacobi polynomials in order to use orthogonality of such polynomials for computing the spectral tau coefficients. Now, making use of (2.14), (2.15) and (3.5), enables one to write
I tν c (x)
∂ 2 u (x, t ) CT Φ L , N (x) I tν ΦτT, M (t ) A D 2 Φ L , N (x) 2 ∂x = CT Φ L , N (x) ΦτT, M (t )PνT AD(2) Φ L , N (x) = ΦτT, M (t )PνT AD(2) Φ L , N (x)Φ LT, N (x)C.
(3.9)
Let
Φ L , N (x)Φ LT, N (x)C = HT Φ L , N (x),
(3.10)
where H is ( N + 1) × ( N + 1) matrix. Then (3.10) can be put in the form N
(α ,β)
ck P L ,k
(α ,β)
(x) P L , j (x) =
k =0
(α ,β)
k =0
(α ,β)
Hkj P L ,k
(x),
j = 0, 1 , . . . , N .
k =0
Multiplying both sides by P L ,i N
N
L ck
(α ,β)
P L ,i
(α ,β)
(α ,β)
(x) w L
(α ,β)
(x), i = 0, 1, . . . , N and integrating from 0 to L, we obtain (α ,β)
(x) P L ,k (x) P L , j (x) w L
(x) dx =
N k =0
0
L Hi j
(α ,β)
P L ,k
(α ,β)
(x) P L ,i
(α ,β)
(x) w L
(x) dx.
(3.11)
0
Making use of (2.5), yields
1
L N
Hi j = (α ,β) ck h L , i k =0
(α ,β) (α ,β) (α ,β) (α ,β) P L ,i (x) P L ,k (x) P L , j (x) w L (x) dx
,
i , j = 0, 1 , . . . , N .
(3.12)
0
Now, (3.9) and (3.10) give
I tν c (x)
∂ 2 u (x, t ) = ΦτT, M (t )PνT AD(2) HT Φ L , N (x). ∂ x2
(3.13)
Eq. (3.5) with (3.13) enable us to write the residual R N , M (x, t ) for Eq. (3.4) in the form
R M , N (x, t ) = ΦτT, M (t ) A − Q − PνT AD(2) H T Φ L , N (x)
= ΦτT, M (t )EΦ L , N (x),
(3.14)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.7 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
7
where
E = A − Q − PνT AD(2) H T . As in a typical tau method, see [3], we generate ( M + 1) × ( N − 1) linear algebraic equations in the unknown expansion coefficients, ai j , i = 0, 1, . . . , M; j = 0, 1, . . . , N − 2, namely
τ L
(α ,β)
R N , M (x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
(x) dx dt = 0,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N − 2 ,
(3.15)
and the rest of linear algebraic equations are obtained from the boundary conditions (3.3), as
ΦτT, M (t )AΦ L , N (0) = g 0 (t ), ΦτT, M (t )AΦ L , N ( L ) = g 1 (t ),
(3.16) α ,β
respectively. Eqs. (3.16) are collocated at ( M + 1) points. For suitable collocation points we use the shifted Jacobi roots t N ,i , (α ,β) i = 0, 1, . . . , M of P τ , M +1 (t ). The number of the unknown coefficients ai j is equal to ( N + 1) × ( M + 1) and can be obtained from Eqs. (3.15)–(3.16). Consequently u N , M (x, t ) given in Eq. (3.5) can be calculated. Remark 3.1. The cases α = β = 0 and α = β = −1/2 lead to the Legendre and Chebyshev tau methods combined with operational matrices respectively, which are special cases from Jacobi tau method and used most frequently in practice. It is actually more inclusive to study the whole family of Jacobi polynomials with the general parameters α , β > −1 rather than the Legendre or Chebyshev polynomials. 4. Fractional wave equation with damping In this section, we consider the time-fractional diffusion-wave equation with damping, namely
∂ u (x, t ) ∂ 2 u (x, t ) ∂ ν u (x, t ) +λ = + μs(x, t ), ν ∂t ∂t ∂ x2
0 < x L , 0 < t τ , 1 < ν 2,
(4.1)
with initial conditions (3.2) and boundary conditions (3.3). Reason for considering the fractional diffusion-wave equation with damping is that it has the ability of governing wave motion in a string and it arising from fractional diffusion-wave equation by assuming the damping effect due to the term λ ∂ u∂(xt ,t ) . In case of ν = 2, this equation can also be characterized as a telegraph equation which governs electrical transmission in a telegraph cable. After applying the Riemann–Liouville integral of order ν on Eq. (4.1) and making use of (2.4), we get the integrated form, namely ν −1
u (x, t ) − v (x, t ) + λ I t u (0, t ) = g 0 (t ),
∂ 2 u (x, t ) u (x, t ) = I tν ∂ x2
+ μ I tν s(x, t ),
u ( L , t ) = g 1 (t ),
(4.2)
ν −1
where v (x, t ) = q(x, t ) − λ I t f 0 (x). In order to use the tau method with the shifted Jacobi operational matrix for solving the fully integrated problem (4.2), we approximate u (x, t ) and q(x, t ) by the shifted Jacobi polynomials as in Eqs. (3.5) and (3.6) but s(x, t ) will be approximated in the form
s M , N (x, t ) = ΦτT, M (t )SΦ L , N (x),
(4.3)
where S is a known ( M + 1) × ( N + 1) matrix that can be written as
⎛
s00 ⎜ s10
S=⎜ ⎝ ..
s01 s11
.
.. .
s M0
s M1
··· ···
⎞
s0N s1N ⎟
.. ⎟ ⎠, . ··· · · · sM N
where si j are obtained by
(4.4)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.8 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
8
1 si j = (α ,β) (α ,β) h τ ,i h L , j
τL
(α ,β)
s(x, t ) P τ ,i
(α ,β)
(α ,β)
(t ) P L , j (x) w τ
(α ,β)
(t ) w L
(x) dx dt ,
0 0
i = 0, 1 , . . . , M , j = 0, 1 , . . . , N .
(4.5)
Using Eqs. (2.14), (2.15), (3.5) and (4.2), then it is easy to write
I tν −1 u (x, t ) = φτT , M (t )PνT −1 Aφ L , N (x),
I tν
∂ 2 u (x, t ) = φτT, M (t )PνT AD(2) φ L , N (x), ∂ x2
I tν s(x, t ) = φτT , M (t )PνT Sφ L , N (x).
(4.6)
Eqs. (3.5), (4.3) with (4.6), enable us to write the residual R N , M (x, t ) for Eq. (4.2) in the form
R N , M (x, t ) = φτT , M (t ) A − 2V + λPνT −1 A − PνT AD(2) − μPνT S φ L , N (x) = φτT , M (t )Eφ L , N (x),
(4.7)
where
E = A − 2V + λPνT −1 A − PνT AD(2) − μPνT S. We generate ( M + 1) × ( N − 1) linear algebraic equations in the unknown expansion coefficients, ai j , i = 0, 1, . . . , M; j = 0, 1, . . . , N − 2, see (3.15). As in Section 3, the rest of this linear algebraic equations are obtained from the boundary conditions (3.3), see (3.16). Then the approximate solution u N , M (x, t ) can be calculated by solving the generated system of ( M + 1) × ( N + 1) linear algebraic equations. 5. Fourth-order diffusion-wave equation In this section, we consider the fourth-order fractional diffusion-wave equation of the form:
∂ ν u (x, t ) ∂ 4 u (x, t ) + = s(x, t ), ∂tν ∂ x4
0 < x L , 0 < t τ , 1 < ν 2,
(5.1)
with initial conditions (3.2) and boundary conditions
u (0, t ) = g 0 (t ),
u ( L , t ) = g 1 (t ),
u xx (0, t ) = g 2 (t ),
u xx ( L , t ) = g 3 (t ).
(5.2)
We apply the same technique described in Section 3, the residual R N , M (x, t ) for Eq. (5.1) can be written as
R N , M (x, t ) = φτT , M (t ) A − Q + PνT AD(4) − PνT S φ L , N (x) = φτT , M (t )Eφ L , N (x),
(5.3)
where
E = A − Q + PνT AD(4) − PνT S. It is clear that, ( M + 1)( N − 3) linear algebraic equations in the unknown expansion coefficients, ai j , i = 0, 1, . . . , M; j = 0, 1, . . . , N − 4, are generated, and the rest of linear algebraic equations are obtained from the boundary conditions (5.2), as
φτT, M (t )Aφ L , N (0) = g 0 (t i ), φτT, M (t )Aφ L , N ( L ) = g 1 (t i ), φτT, M (t )AD(2) φ L , N (0) = g 2 (t i ), φτT, M (t )AD(2) φ L , N ( L ) = g 3 (t i ).
(5.4)
After collocating Eqs. (5.4) at ( M + 1) points with shifted Jacobi roots, ( M + 1) × ( N + 1) linear algebraic equations will be generated using Eqs. (3.15) and (5.4). Consequently u N , M (x, t ) given in Eq. (3.5) can be calculated.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.9 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
Table 1 MAEs at t = 0.8 with N = M = 8 and various choices of
9
ν , α and β for Example 1.
ν
α=β =0
α=β =1
α = β = 1/2
α = β = −1/2
9/8 10/8 11/8 12/8 13/8 14/8 15/8 16/8
3.55 · 10−7 5.29 · 10−7 5.44 · 10−7 4.45 · 10−7 3.08 · 10−7 1.96 · 10−7 1.06 · 10−7 3.23 · 10−12
2.86 · 10−7 4.19 · 10−7 4.17 · 10−7 3.39 · 10−7 2.80 · 10−7 2.15 · 10−7 7.19 · 10−8 5.36 · 10−12
3.39 · 10−7 4.95 · 10−7 4.94 · 10−7 3.94 · 10−7 2.96 · 10−7 2.13 · 10−7 1.20 · 10−7 2.28 · 10−12
3.24 · 10−7 4.23 · 10−7 1.07 · 10−7 4.61 · 10−7 3.17 · 10−7 1.83 · 10−7 1.07 · 10−7 1.14 · 10−11
Table 2 Absolute errors at x = 1 with N = M = 15,
α = β = 1 and various choices of ν for Example 1.
t
ν = 5/4
ν = 6/4
ν = 7/4
ν = 8/4
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
3.36 · 10−8 3.65 · 10−9 5.31 · 10−9 5.67 · 10−9 4.85 · 10−9 3.94 · 10−9 2.50 · 10−9 1.09 · 10−9 5.37 · 10−9
1.69 · 10−8 1.86 · 10−9 2.73 · 10−9 2.93 · 10−9 2.51 · 10−9 2.05 · 10−9 1.30 · 10−9 5.68 · 10−10 2.80 · 10−9
4.27 · 10−9 4.79 · 10−10 7.07 · 10−10 7.61 · 10−10 6.55 · 10−10 5.35 · 10−10 3.40 · 10−10 1.48 · 10−10 7.34 · 10−10
1.48 · 10−15 1.78 · 10−15 1.60 · 10−15 1.43 · 10−15 1.67 · 10−15 1.67 · 10−15 1.52 · 10−15 1.46 · 10−15 9.61 · 10−16
6. Numerical results To demonstrate the effectiveness of the proposed method in the present paper, five test examples are carried out in this section. Example 1. We consider the following problem [19]
∂ ν u (x, t ) 1 2 ∂ 2 u (x, t ) = x , ∂tν 2 ∂ x2
0 < x 1, 0 < t 1, 1 < ν 2,
(6.1)
with the initial conditions u (x, 0) = x, ut (x, 0) = x2 and boundary conditions u (0, t ) = 0, u (1, t ) = sinh(t ), while the exact solution of this initial value problem for ν = 2, is u (x, t ) = x + x2 sinh(t ), while the analytical solution using the VIM in [29] is given by
u (x, t ) = x + x2
∞ k =0
t k ν +1
Γ (kν + 2)
.
The maximum absolute errors (MAEs) between the exact solution u (x, 0.8) and the approximate solution u 8,8 (x, 0.8) with various choices of Jacobi parameters α , β and four choices of the fractional derivative ν are given in Table 1. In addition, Table 2 also shows the absolute error function |u (1, t ) − u 15,15 (1, t )| with α = β = 1 and various choices of ν , while Fig. 1 plots the error function u (x, t ) − u 10,10 (x, t ) with ν = 2 and α = β = 0. Tables 1 and 2 show the good agreement with the exact solutions obtained by the VIM [19] for different choices of x, t and ν . Example 2. Consider the one-dimensional time-fractional wave equation [30]
∂ 3/2 u (x, t ) ∂ 2 u (x, t ) = , ∂ t 3/2 ∂ x2
0 < x 1, 0 < t 1,
(6.2)
with initial conditions
u (x, 0) = sin(x),
∂ u (x, 0) = − sin(x), ∂t
and boundary conditions
u (0, t ) = 0,
u (1, t ) = sin(1) E 3/2,1 −t 3/2 − t E 3/2,2 −t 3/2 .
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.10 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
10
Fig. 1. The error at
Table 3 Absolute errors with N = M = 10,
ν = 2 with N = M = 10 and α = β = 0 for Example 1.
α = β = 0 and various choices of x and t for Example 2.
t
x = 0.1
x = 0.3
x = 0.5
x = 0.7
x = 0.9
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.19 · 10−6 5.88 · 10−7 1.14 · 10−6 8.97 · 10−7 6.51 · 10−7 2.24 · 10−7 8.43 · 10−7 9.53 · 10−7 1.02 · 10−6
3.48 · 10−6 1.59 · 10−6 3.42 · 10−6 2.82 · 10−6 1.77 · 10−6 7.04 · 10−7 2.29 · 10−6 2.66 · 10−6 3.39 · 10−6
6.05 · 10−6 2.38 · 10−6 5.51 · 10−6 5.04 · 10−6 2.55 · 10−6 1.57 · 10−6 3.64 · 10−6 3.57 · 10−6 4.25 · 10−6
8.56 · 10−6 3.22 · 10−6 7.07 · 10−6 7.62 · 10−6 3.90 · 10−6 2.56 · 10−6 6.39 · 10−6 6.14 · 10−6 5.36 · 10−6
6.17 · 10−6 1.15 · 10−6 5.00 · 10−6 8.32 · 10−6 7.61 · 10−6 2.21 · 10−6 4.48 · 10−6 9.57 · 10−6 1.37 · 10−5
Fig. 2. The error function at N = M = 15 and
α = β = −0.5 for Example 2.
Momani et al. [30] applied the generalized two-dimensional differential transform to both sides of problem (6.2) and showed that its exact solution is E 3/2,1 (−t 3/2 ) sin(x) − t E 3/2,2 (−t 3/2 ) sin(x), where
E ζ,η ( z) =
∞ k =0
zk
Γ (ζ k + η)
.
We solved this problem for different values of N, M and ν and compared our results with the exact solution in Table 3. Fig. 2 shows the error function u (x, t ) − u 15,15 (x, t ) with α = β = −0.5, while Fig. 3 shows the error functions u (x, 0.6) − u 10,10 (x, 0.6) and u (0.1, t ) − u 10,10 (0.1, t ) with α = β = 0.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.11 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
Fig. 3. Errors at t = 0.6 (left) and x = 0.1 (right) with N = M = 10, Table 4 MAEs at t = 1, L = 2,
11
α = β = 0 for Example 2.
τ = 1 with α = β = 0 and various choices of N, M and ν for Example 3.
N=M
ν = 1.1
ν = 1.5
ν = 1.7
ν = 1.9
3 5 7 9 11 13 15 17 19 21 23
1.74 · 10−4 3.76 · 10−5 1.30 · 10−5 5.73 · 10−6 2.90 · 10−6 1.63 · 10−6 9.94 · 10−7 6.26 · 10−7 4.30 · 10−7 3.00 · 10−7 2.16 · 10−7
6.17 · 10−4 1.44 · 10−4 5.84 · 10−5 2.76 · 10−5 1.51 · 10−5 9.03 · 10−6 5.77 · 10−6 3.88 · 10−6 2.72 · 10−6 1.97 · 10−6 1.46 · 10−6
2.72 · 10−4 6.98 · 10−5 3.15 · 10−5 1.55 · 10−5 9.15 · 10−6 5.70 · 10−6 3.80 · 10−6 2.64 · 10−6 1.91 · 10−6 1.41 · 10−6 9.15 · 10−7
6.44 · 10−5 8.57 · 10−6 1.92 · 10−6 2.96 · 10−7 2.29 · 10−7 2.73 · 10−7 2.66 · 10−7 2.26 · 10−7 1.91 · 10−7 1.59 · 10−7 1.32 · 10−7
Example 3. We consider the following fractional wave equation with damping
∂ ν u (x, t ) ∂ u (x, t ) ∂ 2 u (x, t ) + = + q(x, t ), ∂tν ∂t ∂ x2
0 < x L , 0 < t τ , 1 < ν 2,
(6.3)
with initial conditions
u (x, 0) = 0,
∂ u (x, 0) = 0, ∂t
and boundary equations
u (0, t ) = 0,
u ( L , t ) = 0,
and
q(x, t ) =
2x( L − x) 2−ν t + 2tx( L − x) + 2t 2 . Γ (3 − ν )
The exact solution of this problem is u (x, t ) = t 2 x( L − x). In [5], Chen et al. applied the method of separation of variables with constructing the implicit difference approximation to introduce an approximate solution of this problem at L = 2, T = 1 and ν = 1.7 using the Caputo form of the fractional derivative. Regarding problem (6.3), in [5], the best result is achieved at L = 2, T = 1 and ν = 1.7 with 640 steps and the maximum absolute error is 5.31 × 10−5 . In Table 4, we list MAEs at t = 1, L = 2, τ = 1 with α = β = 0 and various choices of N, M and ν . Our method is more accurate than the method of separation of variables combined with the implicit difference approximation [5]. We see in this table that the results are accurate for even small choices of N and M. In addition, to demonstrate the convergence of the proposed method, in Fig. 4, we plot the logarithmic graphs of MAEs (log10 Error) at τ = 1, ν = 1.11, 1.55, 1.77, 1.88 and α = β = 1 with various values of N (N = M); by using the presented algorithm. Clearly, the numerical errors decay rapidly as N and M increase. In order to confirm the accuracy and applicability of the proposed method for all the Jacobi parameters α , β , we introduce in Table 5, the MAEs at t = 1, L = 1, τ = 1 with N = M = 8 and various choices of α , β , and four different choices of ν . From the numerical results of this table, it is observed that neither Chebyshev polynomial (α = β = − 12 ) nor Legendre polynomial (α = β = 0) is the best special case of Jacobi polynomials. In general, it can be concluded that, only a limited number of shifted Jacobi polynomials is needed to obtain a satisfactory result. This numerical experiment demonstrates the utility of the method.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.12 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
12
Fig. 4. Convergence of problem (6.3) at
Table 5 MAEs at t = 1, L = 1,
α 0 1 1/2 −1/2 −1/2
τ = 1, ν = 1.11, 1.55, 1.77, 1.88.
τ = 1 with N = M = 8 and various choices of α , β and ν for Example 3.
β
ν = 1.2
ν = 1.4
ν = 1.6
ν = 1.8
0 1 1/2 −1/2 1/2
5.35 · 10−6
1.01 · 10−5
8.77 · 10−6
2.27 · 10−6 9.75 · 10−6 3.88 · 10−6 7.11 · 10−6 7.36 · 10−6
5.34 · 10−6 5.59 · 10−6 4.00 · 10−6 3.93 · 10−7
8.67 · 10−6 9.87 · 10−6 7.98 · 10−6 1.40 · 10−7
2.01 · 10−6 5.67 · 10−6 9.15 · 10−6 3.00 · 10−6
Example 4. Consider the following two-term wave-diffusion equation with damping with 1 < ν < 2:
∂ ν u (x, t ) ∂ u (x, t ) ∂ 2 u (x, t ) + = + q(x, t ), ∂tν ∂t ∂ x2
0 < x 1, 0 < t τ , 1 < ν < 2,
(6.4)
with initial conditions
u (x, 0) = 0,
∂ u (x, 0) = 0, ∂t
and boundary equations
u (0, t ) = t 3 ,
u (1, t ) = et 3 ,
and
q(x, t ) =
6t 3−ν
Γ (4 − ν )
e x + 3t 2 e x − t 3 e x .
The exact solution of this problem is u (x, t ) = e x t 3 . Recently, Liu et al. [23] proposed the fractional predictor–corrector method (FPCM [23]) to solve this problem with several choices of h and τ = h2 , where h and τ are space and time step sizes, respectively. We make a comparison of the presented method at three choices of Jacobi parameters (α , β) with the numerical method proposed in FPCM [23], where the implicit numerical schemes are implemented for space and time discretizations. In case of ν = 1.85, numerical maximum absolute errors with three different M are listed in Table 6. Meanwhile, the results of FPCM [23] are presented in the sixth column of this table. Obviously, our method for all choices of Jacobi parameters is more accurate than FPCM [23]. Moreover, we can see that for this problem the special choice α = β = 0 is the best one. Accordingly, we choose the special case α = β = 0 to plot in Fig. 5, the logarithmic graphs of MAEs (log10 Error) at τ = 1, ν = 1.09, 1.5, 1.85, 1.95, 1.999 with various values of N (N = M); by using the presented algorithm. From this figure, we conclude that the numerical solution is more accurate when ν tends to its integer values and the numerical errors for all chosen ν decay rapidly as N and M increase. We observe also that the suggested algorithm provides accurate and stable numerical results. Example 5. Consider the following fourth-order fractional diffusion-wave equation
∂ ν u (x, t ) ∂ 4 u (x, t ) x Γ (ν + 4) 3 ν +3 + t , + = e + t t ∂tν 6 ∂ x4
0 < x 1, 0 < t 1, 1 < ν 2,
(6.5)
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.13 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
Table 6 Comparing of our scheme with FPCM [23] at
13
ν = 1.85 for Example 4.
N=M
α=β =0
α = β = 0.5
α=β =1
h
FPCM [23]
4 8 16
5.46 · 10−5 5.50 · 10−6 2.42 · 10−7
1.16 · 10−4 1.15 · 10−5 9.75 · 10−7
2.48 · 10−4 2.64 · 10−5 1.25 · 10−6
1/4 1/8 1/16
1.09 · 10−1 2.76 · 10−2 6.72 · 10−3
Fig. 5. Convergence of problem (6.4) at Table 7 Absolute errors at t = 1 with N = M = 9,
τ = 1, ν = 1.09, 1.5, 1.85, 1.95, 1.999.
α = β = 0 and various choices of ν for Example 5.
x
ν = 1.2
ν = 1.3
ν = 1.5
ν = 1.7
ν = 1.8
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.53 · 10−7 2.26 · 10−7 2.88 · 10−7 3.36 · 10−7 3.65 · 10−7 3.75 · 10−7 3.62 · 10−7 3.26 · 10−7 2.69 · 10−7
1.88 · 10−7 2.77 · 10−7 3.54 · 10−7 4.14 · 10−7 4.55 · 10−7 4.72 · 10−7 4.61 · 10−7 4.18 · 10−7 3.42 · 10−7
1.49 · 10−7 2.01 · 10−7 2.50 · 10−7 2.99 · 10−7 3.50 · 10−7 3.95 · 10−7 4.20 · 10−7 4.08 · 10−7 3.47 · 10−7
7.20 · 10−8 7.67 · 10−8 8.47 · 10−8 1.09 · 10−7 1.57 · 10−7 2.20 · 10−7 2.78 · 10−7 2.99 · 10−7 2.61 · 10−7
4.29 · 10−8 3.76 · 10−8 3.44 · 10−8 4.77 · 10−8 8.56 · 10−8 1.43 · 10−7 1.99 · 10−7 2.26 · 10−7 1.95 · 10−7
with initial conditions
u (x, 0) = 0,
∂ u (x, 0) = ex, ∂t
and boundary conditions
u (0, t ) = t ν +3 + t , u xx (0, t ) = t ν +3 + t ,
u (1, t ) = e t ν +3 + t ,
u xx (1, t ) = e t ν +3 + t .
The exact solution of this problem is e x (t ν +3 + t ). Table 7 exhibits the absolute error function |u (x, 1) − u 9,9 (x, 1)|, by applying the shifted Jacobi tau method based on the operational matrices at α = β = 0 and five choices of ν . We present comparisons of our method at α = β = 0 with the numerical method proposed in [18], where the celebrated Crank–Nicolson finite difference (CCNFDM [18]) method is implemented to solve problem (6.5) with different h and τ , where h and τ are space and time step sizes, respectively. In Tables 8 and 9, we introduce maximum absolute errors using the proposed method for different choices of N and those achieved by CCNFDM [18] for various choices of the fractional derivative ν . Again and from numerical results presented in these tables, we clearly observe that the proposed global scheme is accurate for even small values of N and our scheme is more accurate than CCNFDM [18]. 7. Conclusion In this article, the fractional derivatives are described in the Caputo sense to allow traditional initial and boundary conditions to be included in the formulation of the problem. In addition, the Jacobi polynomials are used to approximate the second- and fourth-order fractional diffusion-wave equation and fractional diffusion wave equation with damping and
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.14 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
14
Table 8 Comparing of our scheme with CCNFDM [18] at
ν = 1.3, 1.5, 1.8 for Example 5.
Our scheme (α = β = 0)
CCNFDM [18] with h = 1/500
N=M
ν = 1.3
ν = 1.5
ν = 1.8
τ
ν = 1.3
ν = 1.5
ν = 1.8
4 6 8 10 12
1.80 · 10−3 1.92 · 10−5 1.34 · 10−6 1.59 · 10−7 3.33 · 10−8
3.55 · 10−3 2.41 · 10−5 1.33 · 10−6 1.57 · 10−7 3.36 · 10−8
7.41 · 10−3 1.36 · 10−5 6.57 · 10−7 9.05 · 10−8 1.60 · 10−8
1/5 1/10 1/20 1/40 1/80
7.21 · 10−3 2.25 · 10−3 6.96 · 10−4 2.15 · 10−4 6.63 · 10−5
1.72 · 10−2 6.21 · 10−3 2.21 · 10−3 7.88 · 10−4 2.79 · 10−4
5.73 · 10−2 2.57 · 10−2 1.13 · 10−2 4.99 · 10−3 2.18 · 10−3
Table 9 Comparing of our scheme with CCNFDM [18] at
ν = 1.2, 1.5, 1.7 for Example 5.
Our scheme (α = β = 0) N=M
ν = 1.2
ν = 1.5
ν = 1.7
h
τ = 1/20 000 ν = 1.2
ν = 1.5
ν = 1.7
4 6 8 10 12
1.15 · 10−3 1.48 · 10−5 1.00 · 10−6 1.54 · 10−7 3.18 · 10−8
3.55 · 10−3 2.41 · 10−5 1.33 · 10−6 1.57 · 10−7 3.36 · 10−8
5.97 · 10−3 1.87 · 10−5 8.85 · 10−7 1.09 · 10−7 3.36 · 10−8
1/5 1/10 1/20 1/40 1/80
1.20 · 10−3 3.03 · 10−4 7.66 · 10−5 1.91 · 10−5 4.79 · 10−6
1.18 · 10−3 2.98 · 10−4 7.55 · 10−5 1.89 · 10−5 4.78 · 10−6
1.16 · 10−4 2.94 · 10−4 7.50 · 10−5 1.94 · 10−6 4.90 · 10−6
CCNFDM [18] with
reduce them to systems of linear algebraic equations in the unknown expansion coefficients of the sought-for spectral solutions. The reduction is based on the operational matrix formulation of shifted Jacobi tau method. This proposed method is quite different essentially from those methods presented in literature as it discretized the underline problem spectrally in space and time discretizations. The numerical experiments demonstrated the good accuracy of the described method. Moreover, only a small number of shifted Jacobi polynomials is needed to obtain a satisfactory result. The algorithms presented in this paper can be easily extended to a large amount of space–time fractional differential equations with different boundary conditions. Acknowledgements The authors are very grateful to the referees for carefully reading of the paper and for their comments and suggestions which have improved the paper. References [1] R.T. Baillie, Long memory processes and fractional integration in econometrics, J. Econom. 73 (1996) 5–59. [2] A.H. Bhrawy, M.M. Tharwat, M.A. Alghamdi, A new operational matrix of fractional integration for shifted Jacobi polynomials, Bull. Malays. Math. Soc. (2014), in press. [3] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, New York, 1989. [4] R.C. Charo, C.R. Luz, R. Ikehata, Sharp decay rates for wave equations with a fractional damping via new method in the Fourier space, J. Math. Anal. Appl. 408 (2013) 247–255. [5] J. Chen, F. Liu, V. Anh, S. Shen, Q. Liu, C. Liao, The analytical solution and numerical solution of the fractional diffusion-wave equation with damping, Appl. Math. Comput. 219 (2012) 1737–1748. [6] M.R. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804. [7] M.R. Cui, Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation, Numer. Algorithms 62 (2013) 383–409. [8] R. Darzi, B. Mohammadzade, S. Mousavi, R. Beheshti, Sumudu transform method for solving fractional differential equations and fractional diffusionwave equation, J. Math. Comput. Sci. 6 (2013) 79–84. [9] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order, Comput. Math. Appl. 62 (2011) 2364–2373. [10] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, A new Jacobi operational matrix: an application for solving fractional differential equations, Appl. Math. Model. 36 (2012) 4931–4943. [11] E.H. Doha, A.H. Bhrawy, S.S. Ezz-Eldien, Numerical approximations for fractional diffusion equations via a Chebyshev spectral-tau method, Cent. Eur. J. Phys. 11 (2013) 1494–1503. [12] R. Du, W.R. Cao, Z.Z. Sun, A compact difference scheme for the fractional diffusion-wave equation, Appl. Math. Model. 34 (2010) 2998–3007. [13] G.H. Gao, Z.Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [14] C.F.L. Godinho, J. Weberszpil, J.A. Helayel-Neto, Extending the D’Alembert solution to space–time modified Riemann–Liouville fractional wave equations, Chaos Solitons Fractals 45 (2012) 765–771. [15] J.H. He, Nonlinear oscillation with fractional derivative and its applications, in: International Conference on Vibrating Engineering, Dalian, China, 1998, pp. 288–291. [16] J.H. He, Some applications of nonlinear fractional differential equations and their applications, Bull. Sci. Technol. 15 (2) (1999) 86–90. [17] X. Hu, L. Zhang, A compact finite difference scheme for the fourth-order fractional diffusion-wave system, Comput. Phys. Commun. 182 (2011) 1645–1650. [18] X. Hu, L. Zhang, On finite difference methods for fourth-order fractional diffusion-wave and subdiffusion systems, Appl. Math. Comput. 218 (2012) 5019–5034. [19] H. Jafari, S. Momani, Solving fractional diffusion and wave equations by modified homotopy perturbation method, Phys. Lett. A 370 (2007) 388–396.
JID:YJCPH AID:5161 /FLA
[m3G; v 1.132; Prn:16/04/2014; 7:40] P.15 (1-15)
A.H. Bhrawy et al. / Journal of Computational Physics ••• (••••) •••–•••
15
[20] H. Jiang, F. Liu, I. Turner, K. Burrage, Analytical solutions for the multi-term time-fractional diffusion-wave/diffusion equations in a finite domain, Comput. Math. Appl. 64 (2012) 3377–3388. [21] C.P. Li, A. Chen, J.J. Ye, Numerical approaches to fractional calculus and fractional ordinary differential equation, J. Comput. Phys. 230 (2011) 3352–3368. [22] C.P. Li, Z.G. Zhao, Y.Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (2011) 855–875. [23] F. Liu, M.M. Meerschaert, R.J. McGough, P. Zhuang, Q. Liu, Numerical methods for solving the multi-term time-fractional wave-diffusion equation, Fract. Calc. Appl. Anal. 16 (2013) 9–25. [24] Y. Luke, The Special Functions and Their Approximations, vol. 2, Academic Press, New York, 1969. [25] R.L. Magin, Fractional calculus in bioengineering, Crit. Rev. Biomed. Eng. 32 (1) (2004) 1–104. [26] F. Mainardi, Fractional calculus: ‘Some basic problems in continuum and statistical mechanics’, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer-Verlag, New York, 1997, pp. 291–348. [27] B. Mandelbrot, Some noises with 1/f spectrum, a bridge between direct current and white noise, IEEE Trans. Inf. Theory 13 (2) (1967) 289–298. [28] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37 (2004) 161–208. [29] S. Momani, Z. Odibat, Comparison between the homotopy perturbation method and the variational iteration method for linear fractional partial differential equations, Comput. Math. Appl. 54 (2007) 910–919. [30] S. Momani, Z. Odibat, V.S. Erturk, Generalized differential transform method for solving a space- and time-fractional diffusion-wave equation, Phys. Lett. A 370 (2007) 379–387. [31] R.R. Nigmatullin, To the theoretical explanation of the universal response, Phys. Status Solidi, B Basic Res. 123 (2) (1984) 739–745. [32] R.R. Nigmatullin, Realization of the generalized transfer equation in a medium with fractal geometry, Phys. Status Solidi, B Basic Res. 133 (1) (1986) 425–430. [33] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974. [34] R. Ren, H. Li, W. Jiang, M. Song, An efficient Chebyshev-tau method for solving the space fractional diffusion equations, Appl. Math. Comput. 224 (2013) 259–267. [35] J. Ren, Z.-Z. Sun, X. Zhao, Compact difference scheme for the fractional sub-diffusion equation with Neumann boundary conditions, J. Comput. Phys. 232 (2013) 456–467. [36] Y.A. Rossikhin, M.V. Shitikova, Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids, Appl. Mech. Rev. 50 (1997) 15–67. [37] A. Saadatmandi, M. Dehghan, A tau approach for solution of the space fractional diffusion equation, Comput. Math. Appl. 62 (2011) 1135–1142. [38] Q. Xu, J.S. Hesthaven, Stable multi-domain spectral penalty methods for fractional partial differential equations, J. Comput. Phys. 257 (2014) 241–258. [39] M. Zayernouri, G.E. Karniadakis, Exponentially accurate spectral and spectral element methods for fractional ODEs, J. Comput. Phys. 257 (2014) 460–480. [40] F. Zeng, C. Li, F. Liu, I. Turner, The use of finite difference/element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput. 35 (6) (2013) A2976–A3000.
