- Email: [email protected]
A compact finite difference scheme for the fourth-order fractional diffusion-wave system
Computer Physics Communications 182 (2011) 1645–1650
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
A compact finite difference scheme for the fourth-order fractional diffusion-wave system Xiuling Hu a,b , Luming Zhang a,∗ a b
College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China School of Mathematical Sciences, Xuzhou Normal University, Xuzhou 221116, China
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 1 August 2010 Received in revised form 13 March 2011 Accepted 19 April 2011 Available online 28 April 2011 Keywords: Diffusion-wave system Compact difference scheme Extrapolation Solvability Stability Convergence
A high-order compact finite difference scheme combined with the temporal extrapolation technique is investigated for the fourth-order fractional diffusion-wave system in this paper. The solvability, stability and convergence of the scheme are analyzed simultaneously by the energy method. Numerical experiments show that the proposed compact scheme is more accurate and efficient than the Crank– Nicolson scheme. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Fractional diffusion-wave (FDW) equations have important applications to mathematical physics and have been investigated widely in recent years. The FDW equation is a linear integrodifferential equation obtained from the classical diffusion or wave equation by replacing the first- or second-order time derivative by a fractional derivative of order α > 0 [1,2]. These equations arised from basic random walk models or from generalized master equations [3] and many electromagnetic, acoustic, mechanical responses can be described accurately by the fractional diffusionwave equations [4,5]. Nigmatullin [5] modeled diffusion in media with fractal geometry by the fractional diffusion equation. Ginoa et al. [6] explored a fractional diffusion equation describing relaxation phenomena in complex viscoelastic materials. Huang and Liu [7] investigated the time fractional diffusion equation in n-dimensional whole and half-space. Mbodje and Montseny [8] studied the wave equation with fractional derivative feedback and proved the existence, uniqueness, as well as the asymptotic decay of the solution towards zero and pointed out that the method developed could easily be applied to a wide class of problems involving fractional derivative or integral operators of the time variable. Mainardi [9] constructed the fundamental solution for a fractional order diffusion-wave equation using the method of the Laplace transform. Wess [10] and Schneider and Wess [11] pro-
*
Corresponding author. Tel.: +86 516 83500387. E-mail addresses: [email protected] (X. Hu), [email protected] (L. Zhang).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2011.04.013
©
2011 Elsevier B.V. All rights reserved.
vided solutions of Fox’s H -functions for FDW equations. Rossikhin and Shitikova [12] offered a review application of fractional derivatives and other fractional operators to problems connected with vibrations and waves in solids having hereditarily elastic properties and demonstrated the merits and demerits of the simplest fractional calculus viscoelastic models with numerous examples and compared the results with integer derivatives for the similar viscoelastic models. This review article included 174 references. Increasing interest has been paid to the FDW equations in recent years [13–18]. Many numerical researches have been done considerably on the FDW equations in the meanwhile. Tadjeran et al. [19] presented a Crank–Nicolson method for the one-dimensional fractional diffusion equation and improved the spacial accuracy by using a spacial extrapolation. For the two-dimensional fractional diffusion equation, Tadjeran and Meerschaert [20] presented a numerical method which combined the alternating directions implicit (ADI) method with a Crank–Nicolson discretization and also improved the spacial accuracy by using a Richardson extrapolation. Chen et al. [21,22] proposed a Fourier method for the fractional diffusion equation describing sub-diffusion and the Galilei invariant fractional advection–diffusion equation. Liu et al. [23] established a new implicit difference scheme for a modified anomalous subdiffusion equation with a nonlinear source term and analyzed the stability and convergence by using a new energy method. Sun and Wu [24] investigated a fully discrete scheme by the method of order reduction for a FDW system and proved the solvability, stability and convergence by the energy method. Also many
1646
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
scholars concentrated on the improving of the accuracy and convergence order of the numerical methods. Cui [25] constructed a compact finite difference scheme with spatial accuracy of fourthorder for the one-dimensional fractional diffusion equation. Du et al. [26] provided a compact difference scheme for the FDW system on the basis of [24]. Zhuang et al. [3] considered an anomalous subdiffusion equation and proposed a new implicit numerical method (INM) and two solution techniques for improving the order of convergence of the INM. Liu et al. [27] derived two new implicit numerical methods with convergence order O (τ + h2 ) and O (τ 2 + h2 ) and analyzed the convergence and stability using the energy method and pointed out that these techniques could also be adapted to other anomalous subdiffusions. Chen et al. [28] investigated a variable-order anomalous subdiffusion equation and proposed a numerical scheme with order O (τ + h4 ) and discussed the convergence, stability and solvability by the Fourier analysis. Then they improved the convergence order to O (τ 2 + h4 ). In some applications, a fourth-order space derivative term must be indispensable. For example, fourth-order space derivative terms are required during wave propagation in their formulations in beams and modeling formation of grooves on a flat surface because of grain [29,30]. Jafari et al. [31] solved a fourth-order FDW equation in a bounded domain by decomposition method. Agarwal [32,33] presented a general solution to FDW equations containing fourth-order space derivative in unbounded and bounded domains. Golbabai and Sayevand [34] investigated the generalized fourth-order fractional diffusion-wave equations in the Caputo sense by the homotopy perturbation method and proposed a reliable scheme for calculating nonlinear operators. So far, the finite difference research for the fourth-order FDW system is few to be seen. Because compact finite difference schemes are developed to approximate the second-order derivatives by the fourth-order accuracy and keep the desirable tridiagonal nature, the compact methods have been received great attention [25,26,28]. Here, we explore the compact method for the fourth-order FDW system. In this paper, we consider the mixed initial-boundary value problem of the fourth-order fractional diffusion-wave equation [33]
∂αu ∂ 4u + b2 4 = f (x, t ), x ∈ [0, L ], t ∈ [0, T ], α ∂t ∂x u (0, t ) = u ( L , t ) = 0, ∂ 2 u (0, t ) ∂ 2 u ( L , t ) = = 0, t ∈ [0, T ], ∂ x2 ∂ x2 ∂ u (x, 0) = ψ(x), x ∈ [0, L ], u (x, 0) = φ(x), ∂t
(1)
where h = L / M is the space step and τ = T / N is the time step. Suppose U = {U nj | 0 j M , 0 n N } is a grid function. We introduce the following notations:
unj = u (x j , tn ), n− 12
n
Uj
t 0
∂ 2 u (x, s) ds , 2 ∂s (t − s)α −1
α ∈ (1, 2).
Let
x j = jh,
j = 0, 1 , . . . , M ;
t n = nτ ,
n = 0, 1, . . . , N ,
2 U nj+1 − U nj
=
h
,
W nj ≈ w (x j , tn ),
U nj t¯ =
n
U j x¯ =
1
τ
U nj − U nj −1 ,
U nj − U nj−1 h
,
1 U nj xx¯ = 2 U nj+1 − 2U nj + U nj−1 , h n 1 U j xxx¯ x¯ = 4 U nj+2 − 4U nj+1 + 6U nj − 4U nj−1 + U nj−2 . h Let U h = {U | U = (U 0 , U 1 , . . . , U M ), U 0 = U M = 0}. For V , U ∈ U h , we introduce the following discrete l2 -inner product (·,·), discrete l2 -norm · and discrete l∞ -norm · ∞ :
U n, V n = h
n U
M −1
U nj V nj ,
j =1
∞
=
max
1 j M −1
n 2 n n U = U , U ,
n U . j
In the paper, C denotes a general positive constant, which may have different values in different occurrences. Now we consider the compact scheme for the problem (1)–(3). 6, 3 Suppose (1)–(3) has solution u (x, t ) ∈ C x,t ([0, L ] × [0, T ]). Let
w (x, t ) =
t
1
(2 − α )
0
∂ 2 u (x, s) ds , ∂ s2 (t − s)α −1
(2) (3)
(5)
and
z(x, t ) =
∂ 4 u (x, t ) . ∂ x4
(6)
Then (1) becomes
b2 z(x, t ) = − w (x, t ) + f (x, t ).
(7)
From Taylor expansion and (6), we have
unj
xxx¯ x¯
(4)
When α = 1 and α = 2, Eq. (1) represents a diffusion and a wave equation, respectively. This paper is organized as follows. In Section 2, the spacial high-order compact difference scheme for the problem (1)–(3) is derived. Section 3 is devoted to prove the solvability, stability and convergence of the scheme. In Section 4, the temporal extrapolation technique is introduced. In Section 5, numerical experiments are provided to verify the accuracy and efficiency of the scheme. In the end, a concise conclusion is made. 2. The compact finite difference scheme
x
U nj + U nj −1 ,
where b denotes a constant coefficient, f (x, t ) is a known function. Here we consider the Caputo fractional derivative [35]:
∂ α u (x, t ) 1 = α ∂t (2 − α )
1
=
Uj
U nj ≈ u (x j , tn ),
∂ 4 u (x j , t n ) 1 2 ∂ 6 u (x j , t n ) + h + O h4 4 6 6 ∂x ∂x 2 ∂ z ( x , t ) 1 j n = znj + h2 + O h4 2 6 ∂x 1 2 n h2 ∂ 4 z(η jk , tn ) n + O h4 = z j + h z j xx¯ − 6 12 ∂ x4 1 = znj + h2 znj xx¯ + O h4 =
= where
6
6
znj+1
+ 4znj + znj−1 + O h4 ,
(8)
η jk ∈ (x j−1 , x j+1 ). Hence, taking into (7) and (8), we obtain
n− 12
uj
1
xxx¯ x¯
1 n− 12 n− 1 n− 1 z j +1 + 4z j 2 + z j −12 + O h4 6 1 n− 1 n− 1 n− 1 = − 2 w j +12 + 4w j 2 + w j −12 6b 1 n− 1 n− 1 n− 1 + 2 f j +12 + 4 f j 2 + f j −12 + O h4 . 6b
=
According to [24], we have
(9)
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
n− 1 wj 2
1
=
−
Lemma 3.2. For U n ∈ U h , and (U 0n )xx¯ = (U nM )xx¯ = 0, we have
a0 unj t¯
(2 − α )τ n −1
(an−k−1 − an−k ) ukj t¯ − an−1 ψ j + O τ 3−α , (10)
where
τ 2 −α ( j + 1)2−α − j 2−α , aj = 2−α
j 0.
n U n , U xx x¯ x¯ = h
(11)
=h
Substituting (10) into (9), omitting the small terms and noticing the initial-boundary conditions (2)–(3), we can get the following high-order compact scheme:
6
n− 1
+ W j −12
1
= −b
2
n− 1 U j 2 xxx¯ x¯
U 0n = U nM = 0,
=
n− 1 +4fj 2
+
n− 1 f j −12 ,
1 n N,
0 j M,
(2 − α )τ
=h
1
M −1
M −1
j =1
=h
n −1
(an−k−1 − an−k ) U kj t¯ − an−1 ψ j .
j =2
(15)
k =1
+
From the process of the establish of the compact scheme (12)– (14), we can easily see the local truncation error of the scheme (12)–(14) is R nj = O (τ 3−α + h4 ). Therefore, the scheme proposed here is consistent with the problem (1)–(3).
M −2 j =0
=h
M −1
j =1
3. Solvability, stability and convergence of the compact scheme
⎧ n if j = 0, ⎪ ⎨ U0, 1 2 n n (LU ) j = U j + 6 h (U j )xx¯ , if 1 j M − 1, ⎪ ⎩ n UM if j = M .
=h
1 j M − 1, U 0n = U nM = 0, U 0j = φ j ,
1 n N,
n
1 n N,
0 j M.
(17) (18)
Lemma 3.1. (See [24].) For any G = {G 1 , G 2 , G 3 , . . .} and q, we have
N
n a0 G n − (an−k−1 − an−k )G k − an−1 q G n
n =1
k =1
t 1N−α 2
τ
N n =1
G n2 −
t 2N−α 2(2 − α )
where a j is defined in (11).
q2 ,
U nj−1 xx¯ U nj
U nj xx¯ U nj+1
M −1 j =1
1 h2
U nj xx¯ U nj
M −1 j =1
1 h2
U nj xx¯ U nj
h2
−
1
Un Un h 1 xx¯ 0
M −1
2
U nj xx¯
2
Lemma 3.3. (See [36].) For U n ∈ U h , the following inequality is tenable:
(16)
U 0 xx¯ = U nM xx¯ = 0,
U nj−1 − 2U nj + U nj+1
1
j =1
n− 1 n− 1 = −b2 U j 2 xxx¯ x¯ + L f j 2 ,
U nj xx¯ − U nj−1 xx¯ U nj
U nj xx¯
2 = U xnx¯ .
So the compact scheme (12)–(14) can be rewritten as follows: n− 12
U nj xx¯ x¯ U nj
U nj xx¯ U nj−1 − 2
1 h2
Un Un h M −1 xx¯ M 1 1 + U nM xx¯ U nM −1 + U 0n xx¯ U 1n h h
−
Sometimes, the induction of an operator can make the proof easier [26]. To prove the solvability, stability and convergence of the compact scheme, we also induct a compact operator L:
LW j
h2
h2
−
U nj+1 xx¯ U nj − 2h
1
M 1
x
U nj+1 xx¯ − U nj xx¯ U nj
1
1
M −1
h2
h2
U nj
U nj xx¯
h2
j =1
a0 U nj t¯
xxx¯ x¯
h
j =1
(13)
1
M −1
+h
U nj
j =1 M −1
−h
(14)
1
−
=h
(12)
U 0n xx¯ = U nM xx¯ = 0,
where n− 1 Wj 2
n− 1 f j +12
j =1
6 1 n N,
1 j M − 1, U 0j = φ j ,
+
M −1
j =1
n− 1
n− 12
W j +12 + 4W j
2
U n , (U n )xxx¯ x¯ = U xnx¯ .
Proof.
k =1
1
1647
N = 1, 2, 3, . . . ,
2 sin π2h h
n n U U . xx¯ x¯
Lemma 3.4. (See [24].) For U n ∈ U h , the following inequality holds:
n U
√
∞
L
2
U n . x¯
Theorem 3.1. Suppose φ(x) ∈ C 3 [0, L ], the solution of the difference scheme (16)–(18) depends continuously on the initial data and the inhomogeneous term in the l∞ -norm. Proof. Taking inner product of (16) with hτ L(U nj )t¯ and summing up for n from 1 to m, we can obtain
1648
hτ
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
m M −1
n− 12
LW j
LU nj
m M −1
n− 12
t¯
j =1 n =1
= −b2 hτ
Uj
xxx¯ x¯
n
L Uj
t¯
hτ
n− 1 L f j 2 L U nj t¯ .
j =1 n =1
m M −1
n− 12
LU nj
LW j
j =1 n =1
=
M −1
h
(2 − α )
j =1
(19)
h
−
a0 L U j t¯
(2 − α )
=
2 −α tm
2(2 − α )(2 − α ) 1 −α tm τ
2(2 − α )
−
n =1
m M −1
n− 12
Uj
n =1 j =1
= −b2 hτ
m M −1
Uj
n =1 j =1
= −b2 hτ
n− 12
Uj
n =1 j =1 m M −1
= −b2 hτ
n =1 j =1
− b2 hτ
=−
xx¯
(20)
−
1 U nj t¯ + h2 U nj t¯xx¯ 6
n− 12
1 2
2
1 xxx¯ x¯ 6
xx¯
m 2 2
12
U n 2 − U n−1 2
n =1
xxx¯
[(U nj )xx¯
U nj xx¯ + U nj −1 xx¯ xx¯
U n 2 − U n −1 2 xx¯
U nj xx¯ + U nj −1 xx¯
τ
b h
Lψ2
m n− 12 2 L f . j
n =1
n 1.
(23)
Combining with Lemmas 3.3 and 3.4, for sufficiently small h, we get ∞
2 2 C φx0x¯ + C φxxx¯ 2 + C φ2 + C f n ,
n 1.
(24)
This completes the proof.
2
Proof. Since (16)–(18) is a system of linear algebraic equations at each time level, we consider the following corresponding homogeneous equations:
h2 U nj t¯xx¯
m 2 n =1
6
2
2(2 − α )(2 − α )
n 2 U C φ 0 2 + C φxxx¯ 2 + C φ2 + C f n 2 , xx¯ xx¯
n U
U nj xx¯ t¯
[(U nj )xx¯ − (U nj −1 )xx¯ ]
b
12
2 −α tm
For m is an arbitrary integer, from the above inequality, we have
1
× h2
U 0 2 + b h U 0 2 + xx¯ xxx¯ 2
xxx¯ x¯
Uj
m M −1 1 n =1 j =1
2
xxx¯
2 2
1
b
12
2
α −1 τ + (2 − α )tm
L U nj t¯
n =1 j =1 m M −1
xx¯
Theorem 3.2. The difference scheme (16)–(18) is uniquely solvable.
m M −1
− b2 hτ
2 2 U m 2 + b h U m 2
2
Lψ .
n− 12
(22)
j
n =1
b2
(Lψ j )2
xxx¯ x¯
m n− 12 2 L f .
Substituting (20)–(22) into (19), we obtain
m n 2 L U j t¯
2
2(2 − α )(2 − α )
j t¯
n =1
2
j t¯
2 −α tm
n 2 L U j t¯
m n 2 L U
1
m n 2 L U n =1
2(2 − α )
α −1 τ + (2 − α )tm
According to the definition of the operator L and Lemma 3.2, we have
−b2 hτ
2(2 − α )
n
2
j =1
1 −α tm τ
=
τ
1 −α tm
2
t¯
1 −α tm
m M −1
hτ
m n =1
L U nj t¯
1 1 α −1 L f n− 2 2 + (2 − α )tm j
k =1
(21)
n− 12
Lfj
j =1 n =1
n −1 k n − (an−k−1 − an−k )L U j t¯ − an−1 Lψ j L U j t¯ M −1
12
m M −1
Applying Lemma 3.1 and (15), we have
hτ
2
j =1 n =1
+ hτ
2 2 U m 2 − U 0 2 − b h U m 2 − U 0 2 . xx¯ xxx¯ xx¯ xxx¯
b2
From Young’s inequality, we have
n =1 j =1 m M −1
=−
xxx¯
n− 1 = −b2 U j 2 xxx¯ x¯ , 1 j M − 1, 1 n N , n U 0 xx¯ = U nM xx¯ = 0, 1 n N , U 0n = U nM = 0, n− 12
LW j − (U nj −1 )xx¯ ]
τ
U 0j
= 0,
0 j M,
(25) (26) (27)
n− 1
where W j 2 satisfies (15). From the process of the proof for Theorem 3.1, we can easily know that the solutions of (25)–(27) satisfy
U ∞ = 0, that is
U nj = 0,
0 j M, 1 n N.
This completes the proof.
2
For arbitrary h and τ , the compact scheme (16)–(18) is consistent with the problem (1)–(3) since it has the local truncation error
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
O (τ 3−α + h4 ). Therefor, according to the Lax’s Equivalence Theorem, the compact scheme (16)–(18) converges at the same rate, too. So we have Theorem 3.3. Let u (x, t ) be the solution of the problem (1)–(3), {U nj | 0 j M , 1 n N } be the solution of the scheme (16)–(18) and enj = u (x j , tn ) − U nj . Then for nτ T , we have
n e
∞
4
O τ 3 −α + h .
Table 1 Comparison of some numerical and exact solutions with different M at t = 1. M 32 64 128 256
4. Improving the convergence order in time direction by extrapolation From Theorem 3.3, we know the compact scheme is convergence with the order O (τ 3−α + h4 ). Now we use the extrapolation technique to improve the temporal order of the scheme. The derivation of the scheme shows that the truncation error of the scheme is C τ 3−α + C τ 4−α + O (h4 ). Then if the compact scheme is applied on a coarse grid size h = τ and a fine grid size τ /2, the Richardson extrapolation solution can be obtained (τ /2) − U nj (τ )]/(23−α − 1) with convergence order from [23−α U 2n j
O (τ 4−α + h4 ). Here U 2n (τ /2), U nj (τ ) signify the numerical soluj tions on the fine and coarse grid size, respectively. 5. Numerical experiments
In this section, numerical results are provided to demonstrate the accuracy and efficiency of the compact scheme (12)–(14). In the following experiments, we take b = 1, L = 2, T = 1. The exact solution of the system (1)–(3) with f (x) = 0, α = 1.5 is ∞
u (x, t ) =
2 L
L
E α −b2 a4n4 t 4 sin(anx)
n =1
0
∞ n =0
zn
(αn + 1)
α > 0, z ∈ C.
,
The Crank–Nicolson (C–N) scheme of the problem (1)–(3) is
1
n −1 a0 U j t¯ − (an−k−1 − an−k ) U nj t¯ − an−1 ψ j
(2 − α )τ = −b
2
n
n− 1 U j 2 xxx¯ x¯
U 0j = φ j ,
n− 12
, 1 j M − 1, 1 n N , (28) U 0 xx¯ = U nM xx¯ = 0, 1 n N , (29)
0 j M,
(30)
where a j satisfies (11). Let {U nj | 1 j M − 1, 0 n N } be the solution of the scheme (12)–(14) and denote
e ∞ (τ , h) =
max
1 j M −1
Suppose that
e ∞ (τ , h) = O
( 18 , 1)
( 38 , 1)
( 48 , 1)
Compact scheme C–N scheme Compact scheme C–N scheme Compact scheme C–N scheme Compact scheme C–N scheme Exact solution
−0.435179 −0.434155 −0.435743 −0.435285 −0.435946 −0.435770 −0.436026 −0.435953 −0.436104
−1.239565 −1.235815 −1.240891 −1.239869 −1.24150 −1.240970 −1.241728 −1.241519 −1.241919
−1.577353 −1.574060 −1.579427 −1.577753 −1.580160 −1.579526 −1.580450 −1.580185 −1.580663
Table 2 Comparison of e ∞ by the compact scheme with the C–N scheme at t = 1. N
Compact scheme
C–N scheme
32 64 128 256
4.296e−003 1.843e−003 6.941e−004 2.197e−004
2.863e−002 7.231e−003 1.691e−003 3.915e−004
Table 3a Numerical convergence order of the compact scheme. 8
8
τ, h
e ∞ (τ , h)
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
log2
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
τ = h = 1/3 h = 1/6, τ = 1/19 h = 1/12, τ = 1/121
8.639e−002 9.895e−003 6.991e−004
– 8.731 14.154
– 3.126 3.823
τ = h = 1/4 h = 1/8, τ = 1/25 h = 1/16, τ = 1/159
7.507e−002 6.784e−003 4.638e−004
– 11.066 15.489
– 3.468 3.953
τ = h = 1/8 h = 1/16, τ = 1/51 h = 1/32, τ = 1/324
4.500e−002 2.648e−003 1.636e−004
– 16.993 16.182
– 4.0869 4.0163
Table 3b Numerical convergence order of the C–N scheme. 4
4
τ, h
e ∞ (τ , h)
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
log2
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
τ = h = 1/3 h = 1/6, τ = 1/8 h = 1/12, τ = 1/20
7.499e−002 1.934e−003 7.820e−003
– 3.877 2.473
– 1.955 1.306
τ = h = 1/4 h = 1/8, τ = 1/10 h = 1/16, τ = 1/28
9.792e−002 2.344e−002 5.098e−003
– 4.177 4.606
– 2.063 2.204
τ = h = 1/8 h = 1/16, τ = 1/20 h = 1/32, τ = 1/50
3.666e−002 9.635e−003 2.478e−003
– 3.805 3.888
– 1.928 1.959
k =1
+ fj n
U 0n = U nM = 0,
Numerical scheme
φ(x) sin(anx) dx,
where n is a wave number and E α is the Mittag–Leffler function. The Mittag–Leffler function E α with α > 0 is defined by the following series representation, and valid in the whole complex plane [37]:
E α ( z) =
1649
N u − U N . j j
τ q + hp .
If O (τ q ) ≈ O (h p ), that is τ ≈ h p /q , then e ∞ (τ , h) ≈ C 1 h p or e ∞ (τ , h) ≈ C 2 τ q . So we have e ∞ (2 p /q τ , 2h)/e ∞ (τ , h) ≈ C 1 (2h) p /C 1 h p = 2 p and log2 (e ∞ (2 p /q τ , 2h)/e ∞ (τ , h)) ≈ p.
Example 1. The initial condition φ(x) is taken as [33]
φ(x) =
12x − 4x3 ,
for 0 x 1, −8 + 36x − 24x2 + 4x3 , for 1 x 2.
This function represents the deflection of a simply supported beam subjected to a point load at the center of the beam. Thus the initial data is symmetric. Table 1 gives the comparison of the numerical and exact solutions at some nodes with different M at t = 1. Table 2 shows the e ∞ computed by the compact scheme and the C–N scheme for different mesh sizes at t = 1. Tables 3a and 3b provide the convergence order of the compact scheme and the C–N scheme. From these tables, we can see the compact difference scheme is more accurate and efficient than the C–N scheme. The convergence order of the compact scheme is higher than the C–N scheme, which meets our anticipations. Table 4 gives the results
1650
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
Fig. 1. Numerical solutions (left) of the compact scheme and the exact solutions (right) at t = 1. Table 4 Errors and temporal orders of the compact scheme and the extrapolation.
τ =h
e ∞ (τ , h) (compact)
Order
e ∞ (τ , h) (extrapolation)
Order
1/10 1/20 1/40 1/80
3.124e−002 3.192e−002 4.282e−003 1.506e−003
– 1.390 1.477 1.508
5.649e−003 1.711e−003 4.577e−004 9.219e−005
– 1.723 1.902 2.312
of the compact scheme and extrapolated solutions. We can see the temporal accuracy and order after extrapolation are higher than non-extrapolation. Example 2. Take the initial condition φ(x) as [33]
φ(x) =
21x/4 − 3x3 ,
for 0 x 0.5,
−0.5 + 33x/4 − 6x + x , for 0.5 x 2. 2
3
This function represents the deflection of a simply supported beam subjected to a point load at L /4 from the left end of the beam. Thus the initial conditions are not symmetric. Take h = τ = 1/32. Fig. 1 shows the numerical solutions (left) of the compact scheme and the exact solutions (right) for the fourthorder fractional diffusion system. Fig. 2 plots the curves of the numerical solutions of the compact scheme and the exact solutions at t = 1. From these figures, it is easy to see that the compact scheme constructed in this paper imitates the exact solutions effectively. References [1] F. Mainardi, in: J.I. Wegner, F.R. Norwood (Eds.), Nonlinear Waves in Solids, ASME, Fairfield, NJ, 1995, p. 93. [2] F. Mainardi, F. Mainardi, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wien, 1997, p. 291. [3] P. Zhuang, F. Liu, V. Anh, I. Turnner, SIAM J. Numer. Anal. 46 (2008) 1079. [4] R.R. Nigmatullin, Physica B 123 (1984) 739. [5] R.R. Nigmatullin, Physica B 133 (1986) 425. [6] M. Ginoa, S. Cerbelli, H.E. Roman, Physica A 191 (1992) 449. [7] F. Huang, F. Liu, ANZIAM J. 46 (2005) 1. [8] B. Mbodje, G. Montseny, IEEE Trans. Automat. Control 40 (1995) 378. [9] F. Mainardi, Appl. Math. Lett. 9 (1996) 23. [10] W. Wess, J. Math. Phys. 27 (1996) 2782. [11] W.R. Schneider, W. Wess, J. Math. Phys. 30 (1989) 134. [12] Y.A. Rossikhin, M.V. Shitikova, Appl. Mech. Rev. 50 (1997) 15. [13] Ralf Metzler, Theo F. Nonnenmacher, Chem. Phys. 284 (2002) 67. [14] J. Hossein, D.-G. Varsha, Appl. Math. Comptut. 180 (2006) 488. [15] S. Momani, Z. Odibat, V.S. Erturk, Phys. Lett. A 370 (2007) 379. [16] H. Jafari, S. Seifi, Commun. Nonlinear Sci. Numer. Simulat. 14 (2009) 2006.
Fig. 2. Curves of the numerical solutions of the compact scheme and the exact solutions at t = 1. [17] Y.Z. Povstenko, Physica A 389 (2010) 4696. [18] F. Zhang, X. Jiang, Nonlinear Analysis: Real World Applications 12 (2011) 1841. [19] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, J. Comput. Phys. 213 (2006) 205. [20] C. Tadjeran, M.M. Meerschaert, J. Comput. Phys. 220 (2007) 813. [21] C.M. Chen, F. Liu, I. Turner, V. Anh, J. Comput. Phys. 227 (2007) 886. [22] C.M. Chen, F. Liu, I. Turner, V. Anh, ANZIAM J. 48 (2007) C605. [23] F. Liu, C. Yang, K. Burrage, J. Comput. Appl. Math. 231 (2009) 160. [24] Z.Z. Sun, X. Wu, Appl. Numer. Math. 56 (2006) 193. [25] M. Cui, J. Comput. Phys. 228 (2009) 7792. [26] R. Du, W.R. Cao, Z.Z. Sun, Appl. Math. Modelling 34 (2010) 2998. [27] F. Liu, Q. Yang, I. Turner, in: Proceedings of 2009 Design Engineering Technical Conferences and Computers and Information in Engineering Conference, San Diego, CA, USA, 2009. [28] C. Chen, F. Liu, I. Turner, SIAM J. Scientific Comput. 32 (2010) 1740. [29] K.B. Oldhan, J. Spainer, The Fractional Calculus, Academic Press, New York, 1974. [30] I.N. Sneddon, Fourier Transforms, McGraw–Hill, New York, 1951. [31] H. Jafari, M. Dehghan, K. Sayevand, Numer. Methods Partial Differen. Equat. 24 (2008) 1115. [32] O.P. Agrawal, Fract. Calculus Appl. Anal. 3 (2000) 1. [33] O.P. Agrawal, Comput. Struct. 79 (2001) 1497. [34] A. Golbabai, K. Sayevand, Comput. Math. Appl. 61 (2011) 2227. [35] R. Gorenflo, F. Mainardi, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, New York, 1997, p. 223. [36] R.P. Agarwal, Difference Equations and Inequalities, Marcel Dekker Inc., New York, 1992. [37] F. Mainardi, R. Gorenflo, J. Comput. Appl. Math. 118 (2000) 283.
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
A compact finite difference scheme for the fourth-order fractional diffusion-wave system Xiuling Hu a,b , Luming Zhang a,∗ a b
College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China School of Mathematical Sciences, Xuzhou Normal University, Xuzhou 221116, China
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 1 August 2010 Received in revised form 13 March 2011 Accepted 19 April 2011 Available online 28 April 2011 Keywords: Diffusion-wave system Compact difference scheme Extrapolation Solvability Stability Convergence
A high-order compact finite difference scheme combined with the temporal extrapolation technique is investigated for the fourth-order fractional diffusion-wave system in this paper. The solvability, stability and convergence of the scheme are analyzed simultaneously by the energy method. Numerical experiments show that the proposed compact scheme is more accurate and efficient than the Crank– Nicolson scheme. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Fractional diffusion-wave (FDW) equations have important applications to mathematical physics and have been investigated widely in recent years. The FDW equation is a linear integrodifferential equation obtained from the classical diffusion or wave equation by replacing the first- or second-order time derivative by a fractional derivative of order α > 0 [1,2]. These equations arised from basic random walk models or from generalized master equations [3] and many electromagnetic, acoustic, mechanical responses can be described accurately by the fractional diffusionwave equations [4,5]. Nigmatullin [5] modeled diffusion in media with fractal geometry by the fractional diffusion equation. Ginoa et al. [6] explored a fractional diffusion equation describing relaxation phenomena in complex viscoelastic materials. Huang and Liu [7] investigated the time fractional diffusion equation in n-dimensional whole and half-space. Mbodje and Montseny [8] studied the wave equation with fractional derivative feedback and proved the existence, uniqueness, as well as the asymptotic decay of the solution towards zero and pointed out that the method developed could easily be applied to a wide class of problems involving fractional derivative or integral operators of the time variable. Mainardi [9] constructed the fundamental solution for a fractional order diffusion-wave equation using the method of the Laplace transform. Wess [10] and Schneider and Wess [11] pro-
*
Corresponding author. Tel.: +86 516 83500387. E-mail addresses: [email protected] (X. Hu), [email protected] (L. Zhang).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2011.04.013
©
2011 Elsevier B.V. All rights reserved.
vided solutions of Fox’s H -functions for FDW equations. Rossikhin and Shitikova [12] offered a review application of fractional derivatives and other fractional operators to problems connected with vibrations and waves in solids having hereditarily elastic properties and demonstrated the merits and demerits of the simplest fractional calculus viscoelastic models with numerous examples and compared the results with integer derivatives for the similar viscoelastic models. This review article included 174 references. Increasing interest has been paid to the FDW equations in recent years [13–18]. Many numerical researches have been done considerably on the FDW equations in the meanwhile. Tadjeran et al. [19] presented a Crank–Nicolson method for the one-dimensional fractional diffusion equation and improved the spacial accuracy by using a spacial extrapolation. For the two-dimensional fractional diffusion equation, Tadjeran and Meerschaert [20] presented a numerical method which combined the alternating directions implicit (ADI) method with a Crank–Nicolson discretization and also improved the spacial accuracy by using a Richardson extrapolation. Chen et al. [21,22] proposed a Fourier method for the fractional diffusion equation describing sub-diffusion and the Galilei invariant fractional advection–diffusion equation. Liu et al. [23] established a new implicit difference scheme for a modified anomalous subdiffusion equation with a nonlinear source term and analyzed the stability and convergence by using a new energy method. Sun and Wu [24] investigated a fully discrete scheme by the method of order reduction for a FDW system and proved the solvability, stability and convergence by the energy method. Also many
1646
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
scholars concentrated on the improving of the accuracy and convergence order of the numerical methods. Cui [25] constructed a compact finite difference scheme with spatial accuracy of fourthorder for the one-dimensional fractional diffusion equation. Du et al. [26] provided a compact difference scheme for the FDW system on the basis of [24]. Zhuang et al. [3] considered an anomalous subdiffusion equation and proposed a new implicit numerical method (INM) and two solution techniques for improving the order of convergence of the INM. Liu et al. [27] derived two new implicit numerical methods with convergence order O (τ + h2 ) and O (τ 2 + h2 ) and analyzed the convergence and stability using the energy method and pointed out that these techniques could also be adapted to other anomalous subdiffusions. Chen et al. [28] investigated a variable-order anomalous subdiffusion equation and proposed a numerical scheme with order O (τ + h4 ) and discussed the convergence, stability and solvability by the Fourier analysis. Then they improved the convergence order to O (τ 2 + h4 ). In some applications, a fourth-order space derivative term must be indispensable. For example, fourth-order space derivative terms are required during wave propagation in their formulations in beams and modeling formation of grooves on a flat surface because of grain [29,30]. Jafari et al. [31] solved a fourth-order FDW equation in a bounded domain by decomposition method. Agarwal [32,33] presented a general solution to FDW equations containing fourth-order space derivative in unbounded and bounded domains. Golbabai and Sayevand [34] investigated the generalized fourth-order fractional diffusion-wave equations in the Caputo sense by the homotopy perturbation method and proposed a reliable scheme for calculating nonlinear operators. So far, the finite difference research for the fourth-order FDW system is few to be seen. Because compact finite difference schemes are developed to approximate the second-order derivatives by the fourth-order accuracy and keep the desirable tridiagonal nature, the compact methods have been received great attention [25,26,28]. Here, we explore the compact method for the fourth-order FDW system. In this paper, we consider the mixed initial-boundary value problem of the fourth-order fractional diffusion-wave equation [33]
∂αu ∂ 4u + b2 4 = f (x, t ), x ∈ [0, L ], t ∈ [0, T ], α ∂t ∂x u (0, t ) = u ( L , t ) = 0, ∂ 2 u (0, t ) ∂ 2 u ( L , t ) = = 0, t ∈ [0, T ], ∂ x2 ∂ x2 ∂ u (x, 0) = ψ(x), x ∈ [0, L ], u (x, 0) = φ(x), ∂t
(1)
where h = L / M is the space step and τ = T / N is the time step. Suppose U = {U nj | 0 j M , 0 n N } is a grid function. We introduce the following notations:
unj = u (x j , tn ), n− 12
n
Uj
t 0
∂ 2 u (x, s) ds , 2 ∂s (t − s)α −1
α ∈ (1, 2).
Let
x j = jh,
j = 0, 1 , . . . , M ;
t n = nτ ,
n = 0, 1, . . . , N ,
2 U nj+1 − U nj
=
h
,
W nj ≈ w (x j , tn ),
U nj t¯ =
n
U j x¯ =
1
τ
U nj − U nj −1 ,
U nj − U nj−1 h
,
1 U nj xx¯ = 2 U nj+1 − 2U nj + U nj−1 , h n 1 U j xxx¯ x¯ = 4 U nj+2 − 4U nj+1 + 6U nj − 4U nj−1 + U nj−2 . h Let U h = {U | U = (U 0 , U 1 , . . . , U M ), U 0 = U M = 0}. For V , U ∈ U h , we introduce the following discrete l2 -inner product (·,·), discrete l2 -norm · and discrete l∞ -norm · ∞ :
U n, V n = h
n U
M −1
U nj V nj ,
j =1
∞
=
max
1 j M −1
n 2 n n U = U , U ,
n U . j
In the paper, C denotes a general positive constant, which may have different values in different occurrences. Now we consider the compact scheme for the problem (1)–(3). 6, 3 Suppose (1)–(3) has solution u (x, t ) ∈ C x,t ([0, L ] × [0, T ]). Let
w (x, t ) =
t
1
(2 − α )
0
∂ 2 u (x, s) ds , ∂ s2 (t − s)α −1
(2) (3)
(5)
and
z(x, t ) =
∂ 4 u (x, t ) . ∂ x4
(6)
Then (1) becomes
b2 z(x, t ) = − w (x, t ) + f (x, t ).
(7)
From Taylor expansion and (6), we have
unj
xxx¯ x¯
(4)
When α = 1 and α = 2, Eq. (1) represents a diffusion and a wave equation, respectively. This paper is organized as follows. In Section 2, the spacial high-order compact difference scheme for the problem (1)–(3) is derived. Section 3 is devoted to prove the solvability, stability and convergence of the scheme. In Section 4, the temporal extrapolation technique is introduced. In Section 5, numerical experiments are provided to verify the accuracy and efficiency of the scheme. In the end, a concise conclusion is made. 2. The compact finite difference scheme
x
U nj + U nj −1 ,
where b denotes a constant coefficient, f (x, t ) is a known function. Here we consider the Caputo fractional derivative [35]:
∂ α u (x, t ) 1 = α ∂t (2 − α )
1
=
Uj
U nj ≈ u (x j , tn ),
∂ 4 u (x j , t n ) 1 2 ∂ 6 u (x j , t n ) + h + O h4 4 6 6 ∂x ∂x 2 ∂ z ( x , t ) 1 j n = znj + h2 + O h4 2 6 ∂x 1 2 n h2 ∂ 4 z(η jk , tn ) n + O h4 = z j + h z j xx¯ − 6 12 ∂ x4 1 = znj + h2 znj xx¯ + O h4 =
= where
6
6
znj+1
+ 4znj + znj−1 + O h4 ,
(8)
η jk ∈ (x j−1 , x j+1 ). Hence, taking into (7) and (8), we obtain
n− 12
uj
1
xxx¯ x¯
1 n− 12 n− 1 n− 1 z j +1 + 4z j 2 + z j −12 + O h4 6 1 n− 1 n− 1 n− 1 = − 2 w j +12 + 4w j 2 + w j −12 6b 1 n− 1 n− 1 n− 1 + 2 f j +12 + 4 f j 2 + f j −12 + O h4 . 6b
=
According to [24], we have
(9)
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
n− 1 wj 2
1
=
−
Lemma 3.2. For U n ∈ U h , and (U 0n )xx¯ = (U nM )xx¯ = 0, we have
a0 unj t¯
(2 − α )τ n −1
(an−k−1 − an−k ) ukj t¯ − an−1 ψ j + O τ 3−α , (10)
where
τ 2 −α ( j + 1)2−α − j 2−α , aj = 2−α
j 0.
n U n , U xx x¯ x¯ = h
(11)
=h
Substituting (10) into (9), omitting the small terms and noticing the initial-boundary conditions (2)–(3), we can get the following high-order compact scheme:
6
n− 1
+ W j −12
1
= −b
2
n− 1 U j 2 xxx¯ x¯
U 0n = U nM = 0,
=
n− 1 +4fj 2
+
n− 1 f j −12 ,
1 n N,
0 j M,
(2 − α )τ
=h
1
M −1
M −1
j =1
=h
n −1
(an−k−1 − an−k ) U kj t¯ − an−1 ψ j .
j =2
(15)
k =1
+
From the process of the establish of the compact scheme (12)– (14), we can easily see the local truncation error of the scheme (12)–(14) is R nj = O (τ 3−α + h4 ). Therefore, the scheme proposed here is consistent with the problem (1)–(3).
M −2 j =0
=h
M −1
j =1
3. Solvability, stability and convergence of the compact scheme
⎧ n if j = 0, ⎪ ⎨ U0, 1 2 n n (LU ) j = U j + 6 h (U j )xx¯ , if 1 j M − 1, ⎪ ⎩ n UM if j = M .
=h
1 j M − 1, U 0n = U nM = 0, U 0j = φ j ,
1 n N,
n
1 n N,
0 j M.
(17) (18)
Lemma 3.1. (See [24].) For any G = {G 1 , G 2 , G 3 , . . .} and q, we have
N
n a0 G n − (an−k−1 − an−k )G k − an−1 q G n
n =1
k =1
t 1N−α 2
τ
N n =1
G n2 −
t 2N−α 2(2 − α )
where a j is defined in (11).
q2 ,
U nj−1 xx¯ U nj
U nj xx¯ U nj+1
M −1 j =1
1 h2
U nj xx¯ U nj
M −1 j =1
1 h2
U nj xx¯ U nj
h2
−
1
Un Un h 1 xx¯ 0
M −1
2
U nj xx¯
2
Lemma 3.3. (See [36].) For U n ∈ U h , the following inequality is tenable:
(16)
U 0 xx¯ = U nM xx¯ = 0,
U nj−1 − 2U nj + U nj+1
1
j =1
n− 1 n− 1 = −b2 U j 2 xxx¯ x¯ + L f j 2 ,
U nj xx¯ − U nj−1 xx¯ U nj
U nj xx¯
2 = U xnx¯ .
So the compact scheme (12)–(14) can be rewritten as follows: n− 12
U nj xx¯ x¯ U nj
U nj xx¯ U nj−1 − 2
1 h2
Un Un h M −1 xx¯ M 1 1 + U nM xx¯ U nM −1 + U 0n xx¯ U 1n h h
−
Sometimes, the induction of an operator can make the proof easier [26]. To prove the solvability, stability and convergence of the compact scheme, we also induct a compact operator L:
LW j
h2
h2
−
U nj+1 xx¯ U nj − 2h
1
M 1
x
U nj+1 xx¯ − U nj xx¯ U nj
1
1
M −1
h2
h2
U nj
U nj xx¯
h2
j =1
a0 U nj t¯
xxx¯ x¯
h
j =1
(13)
1
M −1
+h
U nj
j =1 M −1
−h
(14)
1
−
=h
(12)
U 0n xx¯ = U nM xx¯ = 0,
where n− 1 Wj 2
n− 1 f j +12
j =1
6 1 n N,
1 j M − 1, U 0j = φ j ,
+
M −1
j =1
n− 1
n− 12
W j +12 + 4W j
2
U n , (U n )xxx¯ x¯ = U xnx¯ .
Proof.
k =1
1
1647
N = 1, 2, 3, . . . ,
2 sin π2h h
n n U U . xx¯ x¯
Lemma 3.4. (See [24].) For U n ∈ U h , the following inequality holds:
n U
√
∞
L
2
U n . x¯
Theorem 3.1. Suppose φ(x) ∈ C 3 [0, L ], the solution of the difference scheme (16)–(18) depends continuously on the initial data and the inhomogeneous term in the l∞ -norm. Proof. Taking inner product of (16) with hτ L(U nj )t¯ and summing up for n from 1 to m, we can obtain
1648
hτ
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
m M −1
n− 12
LW j
LU nj
m M −1
n− 12
t¯
j =1 n =1
= −b2 hτ
Uj
xxx¯ x¯
n
L Uj
t¯
hτ
n− 1 L f j 2 L U nj t¯ .
j =1 n =1
m M −1
n− 12
LU nj
LW j
j =1 n =1
=
M −1
h
(2 − α )
j =1
(19)
h
−
a0 L U j t¯
(2 − α )
=
2 −α tm
2(2 − α )(2 − α ) 1 −α tm τ
2(2 − α )
−
n =1
m M −1
n− 12
Uj
n =1 j =1
= −b2 hτ
m M −1
Uj
n =1 j =1
= −b2 hτ
n− 12
Uj
n =1 j =1 m M −1
= −b2 hτ
n =1 j =1
− b2 hτ
=−
xx¯
(20)
−
1 U nj t¯ + h2 U nj t¯xx¯ 6
n− 12
1 2
2
1 xxx¯ x¯ 6
xx¯
m 2 2
12
U n 2 − U n−1 2
n =1
xxx¯
[(U nj )xx¯
U nj xx¯ + U nj −1 xx¯ xx¯
U n 2 − U n −1 2 xx¯
U nj xx¯ + U nj −1 xx¯
τ
b h
Lψ2
m n− 12 2 L f . j
n =1
n 1.
(23)
Combining with Lemmas 3.3 and 3.4, for sufficiently small h, we get ∞
2 2 C φx0x¯ + C φxxx¯ 2 + C φ2 + C f n ,
n 1.
(24)
This completes the proof.
2
Proof. Since (16)–(18) is a system of linear algebraic equations at each time level, we consider the following corresponding homogeneous equations:
h2 U nj t¯xx¯
m 2 n =1
6
2
2(2 − α )(2 − α )
n 2 U C φ 0 2 + C φxxx¯ 2 + C φ2 + C f n 2 , xx¯ xx¯
n U
U nj xx¯ t¯
[(U nj )xx¯ − (U nj −1 )xx¯ ]
b
12
2 −α tm
For m is an arbitrary integer, from the above inequality, we have
1
× h2
U 0 2 + b h U 0 2 + xx¯ xxx¯ 2
xxx¯ x¯
Uj
m M −1 1 n =1 j =1
2
xxx¯
2 2
1
b
12
2
α −1 τ + (2 − α )tm
L U nj t¯
n =1 j =1 m M −1
xx¯
Theorem 3.2. The difference scheme (16)–(18) is uniquely solvable.
m M −1
− b2 hτ
2 2 U m 2 + b h U m 2
2
Lψ .
n− 12
(22)
j
n =1
b2
(Lψ j )2
xxx¯ x¯
m n− 12 2 L f .
Substituting (20)–(22) into (19), we obtain
m n 2 L U j t¯
2
2(2 − α )(2 − α )
j t¯
n =1
2
j t¯
2 −α tm
n 2 L U j t¯
m n 2 L U
1
m n 2 L U n =1
2(2 − α )
α −1 τ + (2 − α )tm
According to the definition of the operator L and Lemma 3.2, we have
−b2 hτ
2(2 − α )
n
2
j =1
1 −α tm τ
=
τ
1 −α tm
2
t¯
1 −α tm
m M −1
hτ
m n =1
L U nj t¯
1 1 α −1 L f n− 2 2 + (2 − α )tm j
k =1
(21)
n− 12
Lfj
j =1 n =1
n −1 k n − (an−k−1 − an−k )L U j t¯ − an−1 Lψ j L U j t¯ M −1
12
m M −1
Applying Lemma 3.1 and (15), we have
hτ
2
j =1 n =1
+ hτ
2 2 U m 2 − U 0 2 − b h U m 2 − U 0 2 . xx¯ xxx¯ xx¯ xxx¯
b2
From Young’s inequality, we have
n =1 j =1 m M −1
=−
xxx¯
n− 1 = −b2 U j 2 xxx¯ x¯ , 1 j M − 1, 1 n N , n U 0 xx¯ = U nM xx¯ = 0, 1 n N , U 0n = U nM = 0, n− 12
LW j − (U nj −1 )xx¯ ]
τ
U 0j
= 0,
0 j M,
(25) (26) (27)
n− 1
where W j 2 satisfies (15). From the process of the proof for Theorem 3.1, we can easily know that the solutions of (25)–(27) satisfy
U ∞ = 0, that is
U nj = 0,
0 j M, 1 n N.
This completes the proof.
2
For arbitrary h and τ , the compact scheme (16)–(18) is consistent with the problem (1)–(3) since it has the local truncation error
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
O (τ 3−α + h4 ). Therefor, according to the Lax’s Equivalence Theorem, the compact scheme (16)–(18) converges at the same rate, too. So we have Theorem 3.3. Let u (x, t ) be the solution of the problem (1)–(3), {U nj | 0 j M , 1 n N } be the solution of the scheme (16)–(18) and enj = u (x j , tn ) − U nj . Then for nτ T , we have
n e
∞
4
O τ 3 −α + h .
Table 1 Comparison of some numerical and exact solutions with different M at t = 1. M 32 64 128 256
4. Improving the convergence order in time direction by extrapolation From Theorem 3.3, we know the compact scheme is convergence with the order O (τ 3−α + h4 ). Now we use the extrapolation technique to improve the temporal order of the scheme. The derivation of the scheme shows that the truncation error of the scheme is C τ 3−α + C τ 4−α + O (h4 ). Then if the compact scheme is applied on a coarse grid size h = τ and a fine grid size τ /2, the Richardson extrapolation solution can be obtained (τ /2) − U nj (τ )]/(23−α − 1) with convergence order from [23−α U 2n j
O (τ 4−α + h4 ). Here U 2n (τ /2), U nj (τ ) signify the numerical soluj tions on the fine and coarse grid size, respectively. 5. Numerical experiments
In this section, numerical results are provided to demonstrate the accuracy and efficiency of the compact scheme (12)–(14). In the following experiments, we take b = 1, L = 2, T = 1. The exact solution of the system (1)–(3) with f (x) = 0, α = 1.5 is ∞
u (x, t ) =
2 L
L
E α −b2 a4n4 t 4 sin(anx)
n =1
0
∞ n =0
zn
(αn + 1)
α > 0, z ∈ C.
,
The Crank–Nicolson (C–N) scheme of the problem (1)–(3) is
1
n −1 a0 U j t¯ − (an−k−1 − an−k ) U nj t¯ − an−1 ψ j
(2 − α )τ = −b
2
n
n− 1 U j 2 xxx¯ x¯
U 0j = φ j ,
n− 12
, 1 j M − 1, 1 n N , (28) U 0 xx¯ = U nM xx¯ = 0, 1 n N , (29)
0 j M,
(30)
where a j satisfies (11). Let {U nj | 1 j M − 1, 0 n N } be the solution of the scheme (12)–(14) and denote
e ∞ (τ , h) =
max
1 j M −1
Suppose that
e ∞ (τ , h) = O
( 18 , 1)
( 38 , 1)
( 48 , 1)
Compact scheme C–N scheme Compact scheme C–N scheme Compact scheme C–N scheme Compact scheme C–N scheme Exact solution
−0.435179 −0.434155 −0.435743 −0.435285 −0.435946 −0.435770 −0.436026 −0.435953 −0.436104
−1.239565 −1.235815 −1.240891 −1.239869 −1.24150 −1.240970 −1.241728 −1.241519 −1.241919
−1.577353 −1.574060 −1.579427 −1.577753 −1.580160 −1.579526 −1.580450 −1.580185 −1.580663
Table 2 Comparison of e ∞ by the compact scheme with the C–N scheme at t = 1. N
Compact scheme
C–N scheme
32 64 128 256
4.296e−003 1.843e−003 6.941e−004 2.197e−004
2.863e−002 7.231e−003 1.691e−003 3.915e−004
Table 3a Numerical convergence order of the compact scheme. 8
8
τ, h
e ∞ (τ , h)
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
log2
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
τ = h = 1/3 h = 1/6, τ = 1/19 h = 1/12, τ = 1/121
8.639e−002 9.895e−003 6.991e−004
– 8.731 14.154
– 3.126 3.823
τ = h = 1/4 h = 1/8, τ = 1/25 h = 1/16, τ = 1/159
7.507e−002 6.784e−003 4.638e−004
– 11.066 15.489
– 3.468 3.953
τ = h = 1/8 h = 1/16, τ = 1/51 h = 1/32, τ = 1/324
4.500e−002 2.648e−003 1.636e−004
– 16.993 16.182
– 4.0869 4.0163
Table 3b Numerical convergence order of the C–N scheme. 4
4
τ, h
e ∞ (τ , h)
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
log2
e ∞ (2 3 τ ,2h) e ∞ (τ ,h)
τ = h = 1/3 h = 1/6, τ = 1/8 h = 1/12, τ = 1/20
7.499e−002 1.934e−003 7.820e−003
– 3.877 2.473
– 1.955 1.306
τ = h = 1/4 h = 1/8, τ = 1/10 h = 1/16, τ = 1/28
9.792e−002 2.344e−002 5.098e−003
– 4.177 4.606
– 2.063 2.204
τ = h = 1/8 h = 1/16, τ = 1/20 h = 1/32, τ = 1/50
3.666e−002 9.635e−003 2.478e−003
– 3.805 3.888
– 1.928 1.959
k =1
+ fj n
U 0n = U nM = 0,
Numerical scheme
φ(x) sin(anx) dx,
where n is a wave number and E α is the Mittag–Leffler function. The Mittag–Leffler function E α with α > 0 is defined by the following series representation, and valid in the whole complex plane [37]:
E α ( z) =
1649
N u − U N . j j
τ q + hp .
If O (τ q ) ≈ O (h p ), that is τ ≈ h p /q , then e ∞ (τ , h) ≈ C 1 h p or e ∞ (τ , h) ≈ C 2 τ q . So we have e ∞ (2 p /q τ , 2h)/e ∞ (τ , h) ≈ C 1 (2h) p /C 1 h p = 2 p and log2 (e ∞ (2 p /q τ , 2h)/e ∞ (τ , h)) ≈ p.
Example 1. The initial condition φ(x) is taken as [33]
φ(x) =
12x − 4x3 ,
for 0 x 1, −8 + 36x − 24x2 + 4x3 , for 1 x 2.
This function represents the deflection of a simply supported beam subjected to a point load at the center of the beam. Thus the initial data is symmetric. Table 1 gives the comparison of the numerical and exact solutions at some nodes with different M at t = 1. Table 2 shows the e ∞ computed by the compact scheme and the C–N scheme for different mesh sizes at t = 1. Tables 3a and 3b provide the convergence order of the compact scheme and the C–N scheme. From these tables, we can see the compact difference scheme is more accurate and efficient than the C–N scheme. The convergence order of the compact scheme is higher than the C–N scheme, which meets our anticipations. Table 4 gives the results
1650
X. Hu, L. Zhang / Computer Physics Communications 182 (2011) 1645–1650
Fig. 1. Numerical solutions (left) of the compact scheme and the exact solutions (right) at t = 1. Table 4 Errors and temporal orders of the compact scheme and the extrapolation.
τ =h
e ∞ (τ , h) (compact)
Order
e ∞ (τ , h) (extrapolation)
Order
1/10 1/20 1/40 1/80
3.124e−002 3.192e−002 4.282e−003 1.506e−003
– 1.390 1.477 1.508
5.649e−003 1.711e−003 4.577e−004 9.219e−005
– 1.723 1.902 2.312
of the compact scheme and extrapolated solutions. We can see the temporal accuracy and order after extrapolation are higher than non-extrapolation. Example 2. Take the initial condition φ(x) as [33]
φ(x) =
21x/4 − 3x3 ,
for 0 x 0.5,
−0.5 + 33x/4 − 6x + x , for 0.5 x 2. 2
3
This function represents the deflection of a simply supported beam subjected to a point load at L /4 from the left end of the beam. Thus the initial conditions are not symmetric. Take h = τ = 1/32. Fig. 1 shows the numerical solutions (left) of the compact scheme and the exact solutions (right) for the fourthorder fractional diffusion system. Fig. 2 plots the curves of the numerical solutions of the compact scheme and the exact solutions at t = 1. From these figures, it is easy to see that the compact scheme constructed in this paper imitates the exact solutions effectively. References [1] F. Mainardi, in: J.I. Wegner, F.R. Norwood (Eds.), Nonlinear Waves in Solids, ASME, Fairfield, NJ, 1995, p. 93. [2] F. Mainardi, F. Mainardi, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wien, 1997, p. 291. [3] P. Zhuang, F. Liu, V. Anh, I. Turnner, SIAM J. Numer. Anal. 46 (2008) 1079. [4] R.R. Nigmatullin, Physica B 123 (1984) 739. [5] R.R. Nigmatullin, Physica B 133 (1986) 425. [6] M. Ginoa, S. Cerbelli, H.E. Roman, Physica A 191 (1992) 449. [7] F. Huang, F. Liu, ANZIAM J. 46 (2005) 1. [8] B. Mbodje, G. Montseny, IEEE Trans. Automat. Control 40 (1995) 378. [9] F. Mainardi, Appl. Math. Lett. 9 (1996) 23. [10] W. Wess, J. Math. Phys. 27 (1996) 2782. [11] W.R. Schneider, W. Wess, J. Math. Phys. 30 (1989) 134. [12] Y.A. Rossikhin, M.V. Shitikova, Appl. Mech. Rev. 50 (1997) 15. [13] Ralf Metzler, Theo F. Nonnenmacher, Chem. Phys. 284 (2002) 67. [14] J. Hossein, D.-G. Varsha, Appl. Math. Comptut. 180 (2006) 488. [15] S. Momani, Z. Odibat, V.S. Erturk, Phys. Lett. A 370 (2007) 379. [16] H. Jafari, S. Seifi, Commun. Nonlinear Sci. Numer. Simulat. 14 (2009) 2006.
Fig. 2. Curves of the numerical solutions of the compact scheme and the exact solutions at t = 1. [17] Y.Z. Povstenko, Physica A 389 (2010) 4696. [18] F. Zhang, X. Jiang, Nonlinear Analysis: Real World Applications 12 (2011) 1841. [19] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, J. Comput. Phys. 213 (2006) 205. [20] C. Tadjeran, M.M. Meerschaert, J. Comput. Phys. 220 (2007) 813. [21] C.M. Chen, F. Liu, I. Turner, V. Anh, J. Comput. Phys. 227 (2007) 886. [22] C.M. Chen, F. Liu, I. Turner, V. Anh, ANZIAM J. 48 (2007) C605. [23] F. Liu, C. Yang, K. Burrage, J. Comput. Appl. Math. 231 (2009) 160. [24] Z.Z. Sun, X. Wu, Appl. Numer. Math. 56 (2006) 193. [25] M. Cui, J. Comput. Phys. 228 (2009) 7792. [26] R. Du, W.R. Cao, Z.Z. Sun, Appl. Math. Modelling 34 (2010) 2998. [27] F. Liu, Q. Yang, I. Turner, in: Proceedings of 2009 Design Engineering Technical Conferences and Computers and Information in Engineering Conference, San Diego, CA, USA, 2009. [28] C. Chen, F. Liu, I. Turner, SIAM J. Scientific Comput. 32 (2010) 1740. [29] K.B. Oldhan, J. Spainer, The Fractional Calculus, Academic Press, New York, 1974. [30] I.N. Sneddon, Fourier Transforms, McGraw–Hill, New York, 1951. [31] H. Jafari, M. Dehghan, K. Sayevand, Numer. Methods Partial Differen. Equat. 24 (2008) 1115. [32] O.P. Agrawal, Fract. Calculus Appl. Anal. 3 (2000) 1. [33] O.P. Agrawal, Comput. Struct. 79 (2001) 1497. [34] A. Golbabai, K. Sayevand, Comput. Math. Appl. 61 (2011) 2227. [35] R. Gorenflo, F. Mainardi, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, New York, 1997, p. 223. [36] R.P. Agarwal, Difference Equations and Inequalities, Marcel Dekker Inc., New York, 1992. [37] F. Mainardi, R. Gorenflo, J. Comput. Appl. Math. 118 (2000) 283.