- Email: [email protected]
A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations
Applied Mathematics Letters 26 (2013) 306–312
Contents lists available at SciVerse ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations✩ Qifeng Zhang, Chengjian Zhang ∗ School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China
article
info
Article history: Received 28 July 2012 Received in revised form 22 September 2012 Accepted 23 September 2012 Keywords: Neutral delay parabolic differential equations Compact difference scheme Local truncation error Stability
abstract In this work, we present an implicit compact difference scheme for solving a class of neutral delay parabolic differential equations (NDPDEs). The unique solvability and unconditional stability of the scheme are proved. The temporal accuracy of the scheme is improved by using different Richardson extrapolation techniques for linear and nonlinear problems, and fourth-order accuracy in both temporal and spatial dimensions is obtained. Finally, numerical experiments are conducted to verify the accuracy and efficiency of the algorithms. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction From the twentieth century, research into the theory of delay differential equations (DDEs) has attracted more and more attention from scholars; see e.g. [1,2]. As everyone knows, most DDEs cannot be solved analytically. Developing efficient numerical methods for solving DDEs, especially delay partial differential equations (DPDEs), is necessary and important. To our knowledge, the compact difference scheme combined with extrapolation techniques has been applied to partial differential equations (PDEs), and compact difference methods for solving DPDEs have been developed recently; see e.g. [3–5]. Nevertheless, few numerical works have touched on the NDPDEs. Up to now, the research into such equations has mainly focused on the oscillation, existence and stability of the theoretical solution, which refer to papers [6–9]. Jin et al. [10] studied the NDPDEs with an implicit difference scheme and achieved second-order accuracy in both temporal and spatial dimensions. To improve the numerical accuracy and efficiency, in this work, a linearized difference scheme combined with Richardson extrapolation techniques is considered for the NDPDEs as follows:
2 ∂ u(x, t ) ∂ 2 u(x, t − s) ∂ u(x, t ) ∂ u(x, t − s) + + f (x, t ), + =α ∂t ∂t ∂ x2 ∂ x2 u(x, t ) = g (x, t ), (x, t ) ∈ [a, b] × (−s, 0], u(a, t ) = ua (t ), u(b, t ) = ub (t ), t ∈ (0, T ],
(x, t ) ∈ (a, b) × (0, T ],
where α, s, T are all positive constants, and g (x, t ), ua (t ) and ub (t ) are given smooth functions.
✩ This work was supported by NSFC (Grant No. 11171125, 91130003) and HSF (Grant No. 2011CDB289).
∗
Corresponding author. E-mail addresses: [email protected] (Q. Zhang), [email protected] (C. Zhang).
0893-9659/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.aml.2012.09.015
(1.1)
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
307
The main work is proving that the difference scheme is uniquely solvable and unconditionally stable. Furthermore, Richardson extrapolation techniques are applied to improve the temporal accuracy, and we obtain fourth-order accuracy in both temporal and spatial dimensions. 2. The compact difference scheme and the local truncation error Referring to the work of Sun [4], we firstly divide the region Ω × (0, T ], where Ω = (a, b). Take a positive number M, −a and let h = bM . Suppose the lag s is an integer multiple of the time step τ , τ = ns . Define xi = a + ih, tk = kτ , tk+ 1 =
( + tk+1 ), fi
1 t 2 k T
2
k+ 12
= f (xi , tk+ 1 ). Define Ωhτ = Ωh × Ωτ , where Ωh = {xi |0 ≤ i ≤ M }, Ωτ = {tk | − n ≤ k ≤ N }, N = 2
[ τ ], Uik = u(xi , tk ), 0 ≤ i ≤ M , −n ≤ k ≤ N. Suppose W = {vik |0 ≤ i ≤ M , −n ≤ k ≤ N } is the grid function space defined on Ωhτ . We define 1 k 1 k+1 1 k+ 21 k+ 12 k+ 21 k+ 21 k+ 12 k+1 k δ t vi δx vi vi+1 − vi vi = (vi + vi ), = (vi − vi ), = , 2 τ h 1 k 1 δx vik+ 1 − δx vik− 1 , Avik = (vi+1 + 10vik + vik−1 ). δx2 vik = h
2
12
2
Lemma 2.1 ([4]). Suppose g (x) ∈ C 6 [xi−1 , xi+1 ]. Then 1 12
[g ′′ (xi−1 ) + 10g ′′ (xi ) + g ′′ (xi+1 )] −
1 h2
[g (xi−1 ) − 2g (xi ) + g (xi+1 )] =
h4 240
g (6) (ωi ),
where ω ∈ (xi−1 , xi+1 ). Considering (1.1) at the point (xi , tk+ 1 ), we have 2
2 ∂ 2u ∂ u ∂ k+ 1 + − s + fi 2 . u xi , tk+ 1 + u xi , tk+ 1 − s = α x , t x , t 1 1 i k+ i k+ 2 2 2 2 2 2 ∂t ∂x ∂x
(2.1)
Taylor expansion yields
∂u τ 2 ∂ 3u k+ 1 (xi , ξik ), xi , tk+ 1 = δt Ui 2 − 2 ∂t 24 ∂ t 3 1 ∂ 2u ∂ 2u τ 2 ∂ 4u ∂ 2u = x , t ( x , t ) + ( x , t ) − (xi , ηik ), 1 i k+ i k i k+1 2 2 2 2 ∂x 2 ∂x ∂x 8 ∂ x2 ∂ t 2 ∂u τ 2 ∂ 3u k+ 1 −n k (xi , ξ i ), xi , tk+ 1 − s = δt Ui 2 − 2 ∂t 24 ∂ t 3 1 ∂ 2u ∂ 2u ∂ 2u τ 2 ∂ 4u = ( x , t − s ) + ( x , t − s ) − (xi , ηki ), x , t 1 − s i k i k + 1 i k + 2 ∂ x2 2 ∂ x2 ∂ x2 8 ∂ x2 ∂ t 2
(2.2) (2.3)
(2.4) (2.5)
k
where ξik , ηik ∈ (tk , tk+1 ), ξ i , ηki ∈ (tk − s, tk+1 − s). Substituting (2.2)–(2.5) into (2.1) and acting with A on both sides of (2.1), we get k+ 12
Aδt Ui
k+ 12 −n
+ Aδt Ui
∂ 2u ∂ 2u = A 2 (xi , tk ) + A 2 (xi , tk+1 ) 2 ∂x ∂x 2 a ∂ u ∂ 2u k+ 1 + A 2 (xi , tk − s) + A 2 (xi , tk+1 − s) + Afi 2 + τ 2 Arik 2 ∂x ∂x a
(2.6)
where rik =
1 ∂ 3u
(xi , ξik ) − 3
α ∂ 4u 1 ∂ 3u α ∂ 4u k k ( x , η ) + ( x , ξ ) − (xi , ηki ). i i i i 8 ∂ x2 ∂ t 2 24 ∂ t 3 8 ∂ x2 ∂ t 2
24 ∂ t According to Lemma 2.1 and Taylor expansion, we have
A
h4 ∂ 6 u k ∂ 2u 2 k ( x , t ) = δ U + (θi , tk ), i k x i ∂ x2 240 ∂ x6
θik ∈ (xi−1 , xi+1 ).
(2.7)
(2.8)
Substituting (2.8) into (2.6), we can get
Aδt
k+ 12
Ui
k+ 12 −n
+ Ui
k+ 1 k + 1 −n k+ 1 = α δx2 Ui 2 + δx2 Ui 2 + Afi 2 + Rki ,
(2.9)
308
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
where Rki = τ 2 Arik +
∂ 6u ∂ 6u k+ 21 k+ 12 −n , t , t θ + θ . 1 1 i k+ 2 k+ 2 −n 480 ∂ x6 ∂ x6 i h4
(2.10)
According to (2.10), there exists a constant c such that
|Rki | ≤ c (τ 2 + h4 ), We omit the small term
Aδt
k+ 21
ui
1 ≤ i ≤ M − 1, 0 ≤ k ≤ N − 1.
Rki ,
k+ 12 −n
replace
+ ui
with
uki
(2.11)
in Eq. (2.9), and construct the compact difference scheme as follows:
k+ 1 k+ 1 −n k+ 1 = α δx2 ui 2 + δx2 ui 2 + Afi 2 ,
uk0 = ua (tk ), ukM = ub (tk ), uki = g (xi , tk ),
Uik
(2.12)
1 ≤ k ≤ N,
0 ≤ i ≤ M , − n ≤ k ≤ 0.
For nonlinear f (u, x, t ), we use the following linearized difference scheme (cf. [4]). Suppose that f (u, x, t ) has a first-order continuous derivative with respect to the first component. Taylor expansion at point (xi , tk ) yields
u xi , tk+ 1
2
τ 2 1 ∂ 2 u τ ∂u 3 1 3 − Uik + Uik−1 = u(xi , tk ) + (xi , tk ) + (xi , tk ) + o(τ 3 ) − Uik 2 2 2 ∂t 2 2! ∂ t 2 2 2 2 1 τ ∂ u ∂ u + Uik − τ (xi , tk ) + (xi , tk ) + o(τ 3 ) 2 ∂t 2 ∂t2 3
=
8
τ2
∂ 2u (xi , tk ) + o(τ 3 ). ∂t2
(2.13)
Similarly, Taylor expansion at point ( 23 Uik − 12 Uik−1 , xi , tk ) yields f
u xi , tk+ 1
2
, xi , tk+ 1
=f
2
+
3 2
Uik −
1 2
Uik−1 , xi , tk+ 1
2
u xi , tk+ 1
2
3
1
2
2
− Uik + Uik−1
∂f k ξi , xi , tk+ 1 , 2 ∂u
(2.14)
where ξik is between 23 Uik − 12 Uik−1 and u(xi , tk+ 1 ). 2
Substituting (2.13) into (2.14), we have f
u xi , tk+ 1 2
f
3 2
Uik
−
1 2
, xi , tk+ 1
2
Uik−1
, xi , tk+ 1
+
3τ 2 ∂ 2 u
2
8 ∂t
(xi , tk ) 2
(2.15)
∂f k ξi , xi , tk+ 1 + o(τ 3 ). 2 ∂u
The local truncation error of scheme (2.12), at present, is also O(τ 2 + h4 ). 3. The solvability and stability Let V = {v|v = (v0 , v1 , . . . , vM ), v0 = vM = 0} be the grid function space defined on Ωh . Define ∥v∥ = v ∈ V.
M −1 h
i =1
(vi )2 ,
Lemma 3.1 ([11]). Let Q be a nonnegative definite matrix. For an arbitrary constant λ > 0, we have
∥(λI − Q )(λI + Q )−1 ∥ = ∥(λI + Q )−1 (λI − Q )∥ ≤ 1. We rewrite the compact scheme (2.12) in vector form: 1
(A + rB)(uk+1 + uk+1−n ) = (A − rB)(uk + uk−n ) + AFk+ 2 + dk , ατ 1 5 where r = 2h 2 > 0, A = tridiag( 12 , 6 , k+ 21
ukN −1 )T , F
k+ 12
= (f1
k+ 12
, f2
k+ 21
1 12
(3.1)
) and B = tridiag(−1, 2, −1) are tridiagonal matrices. uk = (uk1 , uk2 , . . . ,
, . . . , fN −1 )T and dk = (dk1 , dk2 , . . . , dkN −1 )T are column vectors. The components of the vector
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
309
dk are, respectively, dk1 = −
1
(u0k+1 − uk0+1−n + uk0+1−n − uk0−n ) + r (uk0+1 + uk0 + u0k+1−n + uk0−n ) +
12 dki = 0, i = 2, . . . , N − 2, dkN −1 = −
1 12
τ 12
(ukN+1 − ukN+1−n + ukN+1−n − ukN−n ) + r (ukN+1 + ukN + uNk+1−n + uNk−n ) +
k+ 12
f0
τ 12
,
k+ 21
fN
.
Define error vectors ek = uk − vk , −n ≤ k ≤ N, where uk is the exact solution of difference scheme (2.12) and vk is the numerical approximation of uk . Suppose there is no error before time level k = 0, i.e. ek = 0,
−n ≤ k ≤ −1.
(3.2)
Then the error equation corresponding to Eq. (2.12) can be written as
(A + rB)(ek+1 + ek+1−n ) = (A − rB)(ek + ek−n ).
(3.3)
It easily see that the coefficient matrix A + rB of system (3.1) is a positive definite matrix; thus, we have the following results. Theorem 3.2. The compact scheme (2.12) has a unique solution. Proof. Mathematical induction will be used to prove this theorem. When −n ≤ k ≤ 0, uk is determined according to the initial condition of (2.12). Suppose that ul has been determined. Let k = l in (3.1). Then we have get the linear algebraic system with respect to ul+1 . Since the coefficient matrix A + rB of this system is a positive definite matrix, there exists a unique solution ul+1 . By the inductive principle, this completes the proof. Theorem 3.3. The compact scheme (2.12) is unconditionally stable, i.e., there exists a constant C independent of the spatial step h and the temporal step τ such that ∥ek ∥ ≤ C ∥e0 ∥, k = 1, 2, . . . , N. Proof. The error equation (3.3) can be rewritten as
(ek+1 + ek+1−n ) = (A + rB)−1 (A − rB)(ek + ek−n ).
(3.4)
Furthermore,
(A + rB)−1 (A − rB) =
1 r
−1
I + A −1 B
1 r
I − A−1 B .
(3.5)
It is easy to see that A−1 B is a positive semidefinite matrix and 1r > 0. According to Lemma 3.1, we can get ∥(A + rB)−1 (A − rB)∥ ≤ 1. Taking the norm on both sides of (3.5), we have the following recurrence relation:
∥ek+1 + ek+1−n ∥ ≤ ∥ek + ek−n ∥.
(3.6)
For any positive number k, there exist nonnegative integers m and p such that k = mL + p, where 0 ≤ p < L. In the following, we will use mathematical induction to prove that ∥emL+p ∥ ≤ C ∥e0 ∥, where 0 ≤ m ≤ [T /s], 0 ≤ p < L. According to (3.6), we have
∥ek + ek−L ∥ ≤ ∥emL+p + e(m−1)L+p ∥ ≤ · · · ≤ ∥emL + e(m−1)L ∥.
(3.7)
Thus
∥emL+p ∥ = ∥emL+p + e(m−1)L+p − e(m−1)L+p ∥ ≤ ∥emL+p + e(m−1)L+p ∥ + ∥e(m−1)L+p ∥ ≤ ∥emL ∥ + ∥e(m−1)L ∥ + ∥e(m−1)L+p ∥.
(3.8)
When m = 0, for any p = 1, 2, . . . , L, combining (3.2) and (3.8), we have
∥ep ∥ ≤ ∥e0 ∥ + ∥e−L ∥ + ∥e−L+p ∥ ≤ 3∥e0 ∥ ≡ γ0 .
(3.9)
Suppose there exists some constant number γm−1 satisfying the following inequality:
∥e(m−1)L+p ∥ ≤ γm−1 ,
(3.10)
where p = 0, 1, 2, . . . , L. Because of (3.9) and (3.10), for any p = 1, 2, . . . , L, we have
∥emL+p ∥ ≤ ∥emL ∥ + ∥e(m−1)L ∥ + ∥e(m−1)L+p ∥ ≤ 3γm−1 ≡ γm .
(3.11)
310
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
Fig. 4.1. Error surface maps of difference scheme (2.12); Example 1 (left), Example 2 (right).
When p = 0, according to (3.10), we have
∥emL+p ∥ = ∥emL ∥ = ∥e(m−1)L+L ∥ ≤ γm−1 ≤ γm .
(3.12)
Thus ∥emL+p ∥ ≤ γm , where p = 0, 1, 2, . . . , L,
γm = 3γm−1 = 3m γ0 = 3m+1 ∥e0 ∥.
(3.13)
By the inductive principle, for any positive integer m, we have
∥emL+p ∥ ≤ 3m+1 ∥e0 ∥,
0 ≤ p < L.
(3.14) [T /s]+1
As a result, when n = mL + p, 0 ≤ m ≤ [T /s], 0 ≤ p < L, according to (3.14), we have ∥e ∥ ≤ 3 C = 3[T /s]+1 , this completes the proof. n
0
∥e ∥. Taking
4. A numerical experiment We define the discrete L2 -norm and L∞ -norm of the compact scheme (2.12) and Richardson extrapolation algorithms, respectively, by
N M E2 (h, τ ) = |eki |2 , k=0 i=0
N M ET 2 (h, τ ) = |ˆeki |2 ,
E∞ (h, τ ) =
max
0≤i≤M ,0≤k≤N
ET ∞ (h, τ ) =
k=0 i=0
max
|eki |,
0≤i≤M ,0≤k≤N
|ˆeki |,
4 2k u (h, τ2 ) − 31 uki (h, τ ) for linear problems and eˆ ki = u(xi , tk ) − 3 i 1 k u4k (h, τ4 ) − 12 u2k (h, τ2 ) + 21 ui (h, τ ) for nonlinear problems, see e.g. [12–14]. We expect compact scheme (2.12) to 21 i 21 i
where eki = u(xi , tk ) − uki (h, τ ), eˆ ki = u(xi , tk ) −
32
be of second-order accuracy and the extrapolation algorithm to be of fourth-order accuracy.
Example 1. Consider the linear NDPDE with a = 0, b = 1, T = 2, s = 0.1, α = 1/2. The initial condition, boundary condition and f (x, t ) are determined by the theoretical solution t 2 cos π x. Example 2. A nonlinear NDPDE is tested with a = 0, b = 1, T = 1, s = 0.1, α = 2/3, whose theoretical solution of the π system is u(x, t ) = (x2 + x + 1) cos t + . 1 f (u, x, t ) =
π (x2 + x + 1) π (x2 + x + 1) π π sin + sin 2 2 (t + 1) t +1 (t − s + 1)2 t −s+1 4 π π 3 π − cos + cos − (x2 + x + 1)2 cos2 . 3 t +1 t −s+1 2 t +1 3
u2 (x, t ) +
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
311
Fig. 4.2. Error curves of difference scheme (2.12) with different sizes of the time step at the final time step; Example 1 (left), Example 2 (right).
Fig. 4.3. The convergence rates of compact scheme (2.12) and the Richardson extrapolation algorithm with h = τ ; Example 1 (left), Example 2 (right).
As the nonlinear term f (u, x, t ) has been linearized using formula (2.15), the corresponding numerical scheme is linear, we expect our theoretical results to also apply to the nonlinear problem. The error surfaces of numerical solutions are shown in Fig. 4.1 for Example 1 (left) and Example 2 (right), where the step 1 1 , τ = 1600 . Fig. 4.2 denotes error curves of the difference scheme (2.12) with different sizes of the time sizes are h = 40 step at the final time step for Example 1 (left) and Example 2 (right) respectively. The convergence rates of compact scheme (2.12) and the extrapolation algorithm with h = τ are shown in Fig. 4.3. The results suggest that the compact scheme is of second-order accuracy and the extrapolation algorithm is of fourth-order accuracy in both time and space with h = τ . It can be concluded that the numerical results are consistent with the theoretical results. References [1] J.H. Wu, Theory and Applications of Partial Functional Differential Equations, Springer-Verlag, New York, 1996. [2] D.F. Li, C.J. Zhang, Nonlinear stability of discontinuous Galerkin methods for delay differential equations, Appl. Math. Lett. 23 (2010) 457–461. [3] L.Z. Qian, H.B. Gu, High order compact scheme combined with extrapolation technique for solving convection–diffusion equations, J. Shandong Univ. Nat. Sci. 46 (12) (2011) 39–43. [4] Z.Z. Sun, A linearized compact difference scheme for a class of nonlinear delay partial differential equations, Appl. Math. Model. (2012) http://dx.doi.org/10.1016/j.bbr.2011.03.031. [5] Y.R. Chi, A high-precision compact difference scheme for a class of nonlinear delay parabolic partial differential equations, J. Yanbian Univ. Nat. Sci. 36 (4) (2010) 287–290. [6] M.K. Grammatikopoulos, Oscillations of first-order neutral delay differential equations, J. Math. Anal. Appl. 120 (2) (1986) 510–520. [7] X.L. Fu, X.Z. Liu, S. Sivaloganathan, Oscillation criteria for impulsive parabolic differential equations with delay, J. Math. Anal. Appl. 268 (2002) 647–664. [8] E. Hernández, H.R. Henríquez, Existence results for partial neutral functional differential equations with unbounded delay, J. Math. Anal. Appl. 221 (2003) 452–475. [9] M. Adimy, K. Ezzinbi, A class of linear partial neutral functional differential equations with nondense domain, J. Differential Equations 147 (1998) 285–332.
312
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
[10] C.R. Jin, Z.H. Yu, R.N. Qu, An implicit difference scheme for solving the neutral delay parabolic differential equation, J. Shandong Univ. Nat. Sci. 46 (8) (2011) 13–16. [11] G.C. Feng, G.P. Yu, J.F. Gui, Introduction to Matrix Iterative Analysis, Jilin University Press, Jilin, 1991 (in Chinese). [12] J. Zhang, X.Y. Geng, R.X. Dai, Analysis on two approaches for high order accuracy finite difference computation, Appl. Math. Lett. 25 (2012) 2081–2085. [13] H.L. Liao, Z.Z. Sun, Maximum norm error estimates of efficient difference schemes for second-order wave equations, J. Comput. Appl. Math. 235 (2011) 2217–2233. [14] C. Zhang, G. Sun, Nonlinear stability of Runge–Kutta methods applied to infinite-delay-differential equations, Math. Comput. Model. 39 (2004) 495–503.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations✩ Qifeng Zhang, Chengjian Zhang ∗ School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China
article
info
Article history: Received 28 July 2012 Received in revised form 22 September 2012 Accepted 23 September 2012 Keywords: Neutral delay parabolic differential equations Compact difference scheme Local truncation error Stability
abstract In this work, we present an implicit compact difference scheme for solving a class of neutral delay parabolic differential equations (NDPDEs). The unique solvability and unconditional stability of the scheme are proved. The temporal accuracy of the scheme is improved by using different Richardson extrapolation techniques for linear and nonlinear problems, and fourth-order accuracy in both temporal and spatial dimensions is obtained. Finally, numerical experiments are conducted to verify the accuracy and efficiency of the algorithms. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction From the twentieth century, research into the theory of delay differential equations (DDEs) has attracted more and more attention from scholars; see e.g. [1,2]. As everyone knows, most DDEs cannot be solved analytically. Developing efficient numerical methods for solving DDEs, especially delay partial differential equations (DPDEs), is necessary and important. To our knowledge, the compact difference scheme combined with extrapolation techniques has been applied to partial differential equations (PDEs), and compact difference methods for solving DPDEs have been developed recently; see e.g. [3–5]. Nevertheless, few numerical works have touched on the NDPDEs. Up to now, the research into such equations has mainly focused on the oscillation, existence and stability of the theoretical solution, which refer to papers [6–9]. Jin et al. [10] studied the NDPDEs with an implicit difference scheme and achieved second-order accuracy in both temporal and spatial dimensions. To improve the numerical accuracy and efficiency, in this work, a linearized difference scheme combined with Richardson extrapolation techniques is considered for the NDPDEs as follows:
2 ∂ u(x, t ) ∂ 2 u(x, t − s) ∂ u(x, t ) ∂ u(x, t − s) + + f (x, t ), + =α ∂t ∂t ∂ x2 ∂ x2 u(x, t ) = g (x, t ), (x, t ) ∈ [a, b] × (−s, 0], u(a, t ) = ua (t ), u(b, t ) = ub (t ), t ∈ (0, T ],
(x, t ) ∈ (a, b) × (0, T ],
where α, s, T are all positive constants, and g (x, t ), ua (t ) and ub (t ) are given smooth functions.
✩ This work was supported by NSFC (Grant No. 11171125, 91130003) and HSF (Grant No. 2011CDB289).
∗
Corresponding author. E-mail addresses: [email protected] (Q. Zhang), [email protected] (C. Zhang).
0893-9659/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.aml.2012.09.015
(1.1)
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
307
The main work is proving that the difference scheme is uniquely solvable and unconditionally stable. Furthermore, Richardson extrapolation techniques are applied to improve the temporal accuracy, and we obtain fourth-order accuracy in both temporal and spatial dimensions. 2. The compact difference scheme and the local truncation error Referring to the work of Sun [4], we firstly divide the region Ω × (0, T ], where Ω = (a, b). Take a positive number M, −a and let h = bM . Suppose the lag s is an integer multiple of the time step τ , τ = ns . Define xi = a + ih, tk = kτ , tk+ 1 =
( + tk+1 ), fi
1 t 2 k T
2
k+ 12
= f (xi , tk+ 1 ). Define Ωhτ = Ωh × Ωτ , where Ωh = {xi |0 ≤ i ≤ M }, Ωτ = {tk | − n ≤ k ≤ N }, N = 2
[ τ ], Uik = u(xi , tk ), 0 ≤ i ≤ M , −n ≤ k ≤ N. Suppose W = {vik |0 ≤ i ≤ M , −n ≤ k ≤ N } is the grid function space defined on Ωhτ . We define 1 k 1 k+1 1 k+ 21 k+ 12 k+ 21 k+ 21 k+ 12 k+1 k δ t vi δx vi vi+1 − vi vi = (vi + vi ), = (vi − vi ), = , 2 τ h 1 k 1 δx vik+ 1 − δx vik− 1 , Avik = (vi+1 + 10vik + vik−1 ). δx2 vik = h
2
12
2
Lemma 2.1 ([4]). Suppose g (x) ∈ C 6 [xi−1 , xi+1 ]. Then 1 12
[g ′′ (xi−1 ) + 10g ′′ (xi ) + g ′′ (xi+1 )] −
1 h2
[g (xi−1 ) − 2g (xi ) + g (xi+1 )] =
h4 240
g (6) (ωi ),
where ω ∈ (xi−1 , xi+1 ). Considering (1.1) at the point (xi , tk+ 1 ), we have 2
2 ∂ 2u ∂ u ∂ k+ 1 + − s + fi 2 . u xi , tk+ 1 + u xi , tk+ 1 − s = α x , t x , t 1 1 i k+ i k+ 2 2 2 2 2 2 ∂t ∂x ∂x
(2.1)
Taylor expansion yields
∂u τ 2 ∂ 3u k+ 1 (xi , ξik ), xi , tk+ 1 = δt Ui 2 − 2 ∂t 24 ∂ t 3 1 ∂ 2u ∂ 2u τ 2 ∂ 4u ∂ 2u = x , t ( x , t ) + ( x , t ) − (xi , ηik ), 1 i k+ i k i k+1 2 2 2 2 ∂x 2 ∂x ∂x 8 ∂ x2 ∂ t 2 ∂u τ 2 ∂ 3u k+ 1 −n k (xi , ξ i ), xi , tk+ 1 − s = δt Ui 2 − 2 ∂t 24 ∂ t 3 1 ∂ 2u ∂ 2u ∂ 2u τ 2 ∂ 4u = ( x , t − s ) + ( x , t − s ) − (xi , ηki ), x , t 1 − s i k i k + 1 i k + 2 ∂ x2 2 ∂ x2 ∂ x2 8 ∂ x2 ∂ t 2
(2.2) (2.3)
(2.4) (2.5)
k
where ξik , ηik ∈ (tk , tk+1 ), ξ i , ηki ∈ (tk − s, tk+1 − s). Substituting (2.2)–(2.5) into (2.1) and acting with A on both sides of (2.1), we get k+ 12
Aδt Ui
k+ 12 −n
+ Aδt Ui
∂ 2u ∂ 2u = A 2 (xi , tk ) + A 2 (xi , tk+1 ) 2 ∂x ∂x 2 a ∂ u ∂ 2u k+ 1 + A 2 (xi , tk − s) + A 2 (xi , tk+1 − s) + Afi 2 + τ 2 Arik 2 ∂x ∂x a
(2.6)
where rik =
1 ∂ 3u
(xi , ξik ) − 3
α ∂ 4u 1 ∂ 3u α ∂ 4u k k ( x , η ) + ( x , ξ ) − (xi , ηki ). i i i i 8 ∂ x2 ∂ t 2 24 ∂ t 3 8 ∂ x2 ∂ t 2
24 ∂ t According to Lemma 2.1 and Taylor expansion, we have
A
h4 ∂ 6 u k ∂ 2u 2 k ( x , t ) = δ U + (θi , tk ), i k x i ∂ x2 240 ∂ x6
θik ∈ (xi−1 , xi+1 ).
(2.7)
(2.8)
Substituting (2.8) into (2.6), we can get
Aδt
k+ 12
Ui
k+ 12 −n
+ Ui
k+ 1 k + 1 −n k+ 1 = α δx2 Ui 2 + δx2 Ui 2 + Afi 2 + Rki ,
(2.9)
308
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
where Rki = τ 2 Arik +
∂ 6u ∂ 6u k+ 21 k+ 12 −n , t , t θ + θ . 1 1 i k+ 2 k+ 2 −n 480 ∂ x6 ∂ x6 i h4
(2.10)
According to (2.10), there exists a constant c such that
|Rki | ≤ c (τ 2 + h4 ), We omit the small term
Aδt
k+ 21
ui
1 ≤ i ≤ M − 1, 0 ≤ k ≤ N − 1.
Rki ,
k+ 12 −n
replace
+ ui
with
uki
(2.11)
in Eq. (2.9), and construct the compact difference scheme as follows:
k+ 1 k+ 1 −n k+ 1 = α δx2 ui 2 + δx2 ui 2 + Afi 2 ,
uk0 = ua (tk ), ukM = ub (tk ), uki = g (xi , tk ),
Uik
(2.12)
1 ≤ k ≤ N,
0 ≤ i ≤ M , − n ≤ k ≤ 0.
For nonlinear f (u, x, t ), we use the following linearized difference scheme (cf. [4]). Suppose that f (u, x, t ) has a first-order continuous derivative with respect to the first component. Taylor expansion at point (xi , tk ) yields
u xi , tk+ 1
2
τ 2 1 ∂ 2 u τ ∂u 3 1 3 − Uik + Uik−1 = u(xi , tk ) + (xi , tk ) + (xi , tk ) + o(τ 3 ) − Uik 2 2 2 ∂t 2 2! ∂ t 2 2 2 2 1 τ ∂ u ∂ u + Uik − τ (xi , tk ) + (xi , tk ) + o(τ 3 ) 2 ∂t 2 ∂t2 3
=
8
τ2
∂ 2u (xi , tk ) + o(τ 3 ). ∂t2
(2.13)
Similarly, Taylor expansion at point ( 23 Uik − 12 Uik−1 , xi , tk ) yields f
u xi , tk+ 1
2
, xi , tk+ 1
=f
2
+
3 2
Uik −
1 2
Uik−1 , xi , tk+ 1
2
u xi , tk+ 1
2
3
1
2
2
− Uik + Uik−1
∂f k ξi , xi , tk+ 1 , 2 ∂u
(2.14)
where ξik is between 23 Uik − 12 Uik−1 and u(xi , tk+ 1 ). 2
Substituting (2.13) into (2.14), we have f
u xi , tk+ 1 2
f
3 2
Uik
−
1 2
, xi , tk+ 1
2
Uik−1
, xi , tk+ 1
+
3τ 2 ∂ 2 u
2
8 ∂t
(xi , tk ) 2
(2.15)
∂f k ξi , xi , tk+ 1 + o(τ 3 ). 2 ∂u
The local truncation error of scheme (2.12), at present, is also O(τ 2 + h4 ). 3. The solvability and stability Let V = {v|v = (v0 , v1 , . . . , vM ), v0 = vM = 0} be the grid function space defined on Ωh . Define ∥v∥ = v ∈ V.
M −1 h
i =1
(vi )2 ,
Lemma 3.1 ([11]). Let Q be a nonnegative definite matrix. For an arbitrary constant λ > 0, we have
∥(λI − Q )(λI + Q )−1 ∥ = ∥(λI + Q )−1 (λI − Q )∥ ≤ 1. We rewrite the compact scheme (2.12) in vector form: 1
(A + rB)(uk+1 + uk+1−n ) = (A − rB)(uk + uk−n ) + AFk+ 2 + dk , ατ 1 5 where r = 2h 2 > 0, A = tridiag( 12 , 6 , k+ 21
ukN −1 )T , F
k+ 12
= (f1
k+ 12
, f2
k+ 21
1 12
(3.1)
) and B = tridiag(−1, 2, −1) are tridiagonal matrices. uk = (uk1 , uk2 , . . . ,
, . . . , fN −1 )T and dk = (dk1 , dk2 , . . . , dkN −1 )T are column vectors. The components of the vector
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
309
dk are, respectively, dk1 = −
1
(u0k+1 − uk0+1−n + uk0+1−n − uk0−n ) + r (uk0+1 + uk0 + u0k+1−n + uk0−n ) +
12 dki = 0, i = 2, . . . , N − 2, dkN −1 = −
1 12
τ 12
(ukN+1 − ukN+1−n + ukN+1−n − ukN−n ) + r (ukN+1 + ukN + uNk+1−n + uNk−n ) +
k+ 12
f0
τ 12
,
k+ 21
fN
.
Define error vectors ek = uk − vk , −n ≤ k ≤ N, where uk is the exact solution of difference scheme (2.12) and vk is the numerical approximation of uk . Suppose there is no error before time level k = 0, i.e. ek = 0,
−n ≤ k ≤ −1.
(3.2)
Then the error equation corresponding to Eq. (2.12) can be written as
(A + rB)(ek+1 + ek+1−n ) = (A − rB)(ek + ek−n ).
(3.3)
It easily see that the coefficient matrix A + rB of system (3.1) is a positive definite matrix; thus, we have the following results. Theorem 3.2. The compact scheme (2.12) has a unique solution. Proof. Mathematical induction will be used to prove this theorem. When −n ≤ k ≤ 0, uk is determined according to the initial condition of (2.12). Suppose that ul has been determined. Let k = l in (3.1). Then we have get the linear algebraic system with respect to ul+1 . Since the coefficient matrix A + rB of this system is a positive definite matrix, there exists a unique solution ul+1 . By the inductive principle, this completes the proof. Theorem 3.3. The compact scheme (2.12) is unconditionally stable, i.e., there exists a constant C independent of the spatial step h and the temporal step τ such that ∥ek ∥ ≤ C ∥e0 ∥, k = 1, 2, . . . , N. Proof. The error equation (3.3) can be rewritten as
(ek+1 + ek+1−n ) = (A + rB)−1 (A − rB)(ek + ek−n ).
(3.4)
Furthermore,
(A + rB)−1 (A − rB) =
1 r
−1
I + A −1 B
1 r
I − A−1 B .
(3.5)
It is easy to see that A−1 B is a positive semidefinite matrix and 1r > 0. According to Lemma 3.1, we can get ∥(A + rB)−1 (A − rB)∥ ≤ 1. Taking the norm on both sides of (3.5), we have the following recurrence relation:
∥ek+1 + ek+1−n ∥ ≤ ∥ek + ek−n ∥.
(3.6)
For any positive number k, there exist nonnegative integers m and p such that k = mL + p, where 0 ≤ p < L. In the following, we will use mathematical induction to prove that ∥emL+p ∥ ≤ C ∥e0 ∥, where 0 ≤ m ≤ [T /s], 0 ≤ p < L. According to (3.6), we have
∥ek + ek−L ∥ ≤ ∥emL+p + e(m−1)L+p ∥ ≤ · · · ≤ ∥emL + e(m−1)L ∥.
(3.7)
Thus
∥emL+p ∥ = ∥emL+p + e(m−1)L+p − e(m−1)L+p ∥ ≤ ∥emL+p + e(m−1)L+p ∥ + ∥e(m−1)L+p ∥ ≤ ∥emL ∥ + ∥e(m−1)L ∥ + ∥e(m−1)L+p ∥.
(3.8)
When m = 0, for any p = 1, 2, . . . , L, combining (3.2) and (3.8), we have
∥ep ∥ ≤ ∥e0 ∥ + ∥e−L ∥ + ∥e−L+p ∥ ≤ 3∥e0 ∥ ≡ γ0 .
(3.9)
Suppose there exists some constant number γm−1 satisfying the following inequality:
∥e(m−1)L+p ∥ ≤ γm−1 ,
(3.10)
where p = 0, 1, 2, . . . , L. Because of (3.9) and (3.10), for any p = 1, 2, . . . , L, we have
∥emL+p ∥ ≤ ∥emL ∥ + ∥e(m−1)L ∥ + ∥e(m−1)L+p ∥ ≤ 3γm−1 ≡ γm .
(3.11)
310
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
Fig. 4.1. Error surface maps of difference scheme (2.12); Example 1 (left), Example 2 (right).
When p = 0, according to (3.10), we have
∥emL+p ∥ = ∥emL ∥ = ∥e(m−1)L+L ∥ ≤ γm−1 ≤ γm .
(3.12)
Thus ∥emL+p ∥ ≤ γm , where p = 0, 1, 2, . . . , L,
γm = 3γm−1 = 3m γ0 = 3m+1 ∥e0 ∥.
(3.13)
By the inductive principle, for any positive integer m, we have
∥emL+p ∥ ≤ 3m+1 ∥e0 ∥,
0 ≤ p < L.
(3.14) [T /s]+1
As a result, when n = mL + p, 0 ≤ m ≤ [T /s], 0 ≤ p < L, according to (3.14), we have ∥e ∥ ≤ 3 C = 3[T /s]+1 , this completes the proof. n
0
∥e ∥. Taking
4. A numerical experiment We define the discrete L2 -norm and L∞ -norm of the compact scheme (2.12) and Richardson extrapolation algorithms, respectively, by
N M E2 (h, τ ) = |eki |2 , k=0 i=0
N M ET 2 (h, τ ) = |ˆeki |2 ,
E∞ (h, τ ) =
max
0≤i≤M ,0≤k≤N
ET ∞ (h, τ ) =
k=0 i=0
max
|eki |,
0≤i≤M ,0≤k≤N
|ˆeki |,
4 2k u (h, τ2 ) − 31 uki (h, τ ) for linear problems and eˆ ki = u(xi , tk ) − 3 i 1 k u4k (h, τ4 ) − 12 u2k (h, τ2 ) + 21 ui (h, τ ) for nonlinear problems, see e.g. [12–14]. We expect compact scheme (2.12) to 21 i 21 i
where eki = u(xi , tk ) − uki (h, τ ), eˆ ki = u(xi , tk ) −
32
be of second-order accuracy and the extrapolation algorithm to be of fourth-order accuracy.
Example 1. Consider the linear NDPDE with a = 0, b = 1, T = 2, s = 0.1, α = 1/2. The initial condition, boundary condition and f (x, t ) are determined by the theoretical solution t 2 cos π x. Example 2. A nonlinear NDPDE is tested with a = 0, b = 1, T = 1, s = 0.1, α = 2/3, whose theoretical solution of the π system is u(x, t ) = (x2 + x + 1) cos t + . 1 f (u, x, t ) =
π (x2 + x + 1) π (x2 + x + 1) π π sin + sin 2 2 (t + 1) t +1 (t − s + 1)2 t −s+1 4 π π 3 π − cos + cos − (x2 + x + 1)2 cos2 . 3 t +1 t −s+1 2 t +1 3
u2 (x, t ) +
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
311
Fig. 4.2. Error curves of difference scheme (2.12) with different sizes of the time step at the final time step; Example 1 (left), Example 2 (right).
Fig. 4.3. The convergence rates of compact scheme (2.12) and the Richardson extrapolation algorithm with h = τ ; Example 1 (left), Example 2 (right).
As the nonlinear term f (u, x, t ) has been linearized using formula (2.15), the corresponding numerical scheme is linear, we expect our theoretical results to also apply to the nonlinear problem. The error surfaces of numerical solutions are shown in Fig. 4.1 for Example 1 (left) and Example 2 (right), where the step 1 1 , τ = 1600 . Fig. 4.2 denotes error curves of the difference scheme (2.12) with different sizes of the time sizes are h = 40 step at the final time step for Example 1 (left) and Example 2 (right) respectively. The convergence rates of compact scheme (2.12) and the extrapolation algorithm with h = τ are shown in Fig. 4.3. The results suggest that the compact scheme is of second-order accuracy and the extrapolation algorithm is of fourth-order accuracy in both time and space with h = τ . It can be concluded that the numerical results are consistent with the theoretical results. References [1] J.H. Wu, Theory and Applications of Partial Functional Differential Equations, Springer-Verlag, New York, 1996. [2] D.F. Li, C.J. Zhang, Nonlinear stability of discontinuous Galerkin methods for delay differential equations, Appl. Math. Lett. 23 (2010) 457–461. [3] L.Z. Qian, H.B. Gu, High order compact scheme combined with extrapolation technique for solving convection–diffusion equations, J. Shandong Univ. Nat. Sci. 46 (12) (2011) 39–43. [4] Z.Z. Sun, A linearized compact difference scheme for a class of nonlinear delay partial differential equations, Appl. Math. Model. (2012) http://dx.doi.org/10.1016/j.bbr.2011.03.031. [5] Y.R. Chi, A high-precision compact difference scheme for a class of nonlinear delay parabolic partial differential equations, J. Yanbian Univ. Nat. Sci. 36 (4) (2010) 287–290. [6] M.K. Grammatikopoulos, Oscillations of first-order neutral delay differential equations, J. Math. Anal. Appl. 120 (2) (1986) 510–520. [7] X.L. Fu, X.Z. Liu, S. Sivaloganathan, Oscillation criteria for impulsive parabolic differential equations with delay, J. Math. Anal. Appl. 268 (2002) 647–664. [8] E. Hernández, H.R. Henríquez, Existence results for partial neutral functional differential equations with unbounded delay, J. Math. Anal. Appl. 221 (2003) 452–475. [9] M. Adimy, K. Ezzinbi, A class of linear partial neutral functional differential equations with nondense domain, J. Differential Equations 147 (1998) 285–332.
312
Q. Zhang, C. Zhang / Applied Mathematics Letters 26 (2013) 306–312
[10] C.R. Jin, Z.H. Yu, R.N. Qu, An implicit difference scheme for solving the neutral delay parabolic differential equation, J. Shandong Univ. Nat. Sci. 46 (8) (2011) 13–16. [11] G.C. Feng, G.P. Yu, J.F. Gui, Introduction to Matrix Iterative Analysis, Jilin University Press, Jilin, 1991 (in Chinese). [12] J. Zhang, X.Y. Geng, R.X. Dai, Analysis on two approaches for high order accuracy finite difference computation, Appl. Math. Lett. 25 (2012) 2081–2085. [13] H.L. Liao, Z.Z. Sun, Maximum norm error estimates of efficient difference schemes for second-order wave equations, J. Comput. Appl. Math. 235 (2011) 2217–2233. [14] C. Zhang, G. Sun, Nonlinear stability of Runge–Kutta methods applied to infinite-delay-differential equations, Math. Comput. Model. 39 (2004) 495–503.