- Email: [email protected]
Convergence and asymptotic stability of Galerkin methods for linear parabolic equations with delays
Applied Mathematics and Computation 264 (2015) 160–178
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Convergence and asymptotic stability of Galerkin methods for linear parabolic equations with delays ✩ Hui Liang∗ School of Mathematical Sciences, Heilongjiang University, Harbin, Heilongjiang 150080, China
a r t i c l e
i n f o
Keywords: Linear parabolic equations Delay Convergence Asymptotic stability Galerkin methods
a b s t r a c t This paper is concerned with the convergence and asymptotic stability of semidiscrete and full discrete schemes for linear parabolic equations with delay. These full discrete numerical processes include forward Euler, backward Euler and Crank–Nicolson schemes. The optimal convergence orders are consistent with those of the original parabolic equation. It is proved that the semidiscrete scheme, backward Euler and Crank–Nicolson full discrete schemes can unconditionally preserve the delay-independent asymptotic stability, but some additional restrictions on time and spatial stepsizes of the forward Euler full discrete scheme is needed to preserve the delay-independent asymptotic stability. Numerical experiments illustrate the theoretical results. © 2015 Elsevier Inc. All rights reserved.
1. Introduction Since there is sometimes a memory or delayed effect, partial functional differential equations are used frequently for describing and modeling mathematically many real life phenomena as demonstrated [18]. The following initial-boundary value problem for linear parabolic differential equations with a delay term:
⎧ ⎨ut (x, t) = auxx (x, t) + buxx (x, t − τ ), u(0, t) = u(π , t) = 0, ⎩ u(x, t) = φ(x, t),
t > 0, x ∈ := [0, π ], t ≥ −τ , −τ ≤ t ≤ 0, x ∈ ,
(1.1)
with a > 0, b ࣙ 0 and τ > 0 has been studied by Houwen et al. [8]. Uniqueness, local existence, and global continuation are straightforward by Green and Stech [6]. In particular, Tian [16] has analyzed the asymptotic stability of finite difference methods applied to (1.1). Applying the process of semidiscretization with respect to the spatial variable x to partial functional differential equations (1.1) results in a system of delay differential equations (DDEs). There are lots of references discussed the numerical methods and numerical stability for DDEs as demonstrated [2,3,9,13], among which we refer the reader to the monograph by Bellen and Zennaro [2]. And there are many references studied the convergence of difference methods for partial functional differential equations as demonstrated [1,7,10]. However, up to now, there are few papers published with respect to the numerical stability
✩ This work is supported by the National Nature Science Foundation of China (no. 11101130), the Research Fund of the Heilongjiang Provincial Education Department for the academic backbone of the excellent young people (no. 1254G044), Science and Technology Innovation Team in Higher Education Institutions of Heilongjiang Province (no. 2014TD005) and the Heilongjiang University Science Funds for Distinguished Young Scholars (no. JCL201303). ∗ Tel.: +86 13936648074. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.amc.2015.04.104 0096-3003/© 2015 Elsevier Inc. All rights reserved.
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
161
analysis for partial functional differential equations. The paper by Tian [16] gives a complete analysis of the numerical stability of finite difference methods, which is easy to carry out, but has more limitations comparing with Galerkin methods; the paper by Liang, Shi and Lv [12] presents a complete analysis of the numerical stability of Galerkin methods for a partial differential equation with piecewise constant argument, which is a partial differential equation with a special delay term, and Liang, Liu and Lv [11] also investigate the numerical stability of θ -methods for this equation. The method of lines and its stability for another type of linear parabolic differential equations with a delay term has been discussed by Zubik-Kowal [19,20]. The delay-dependent stability for a two-dimensional parabolic functional-differential equation with both fixed and distributed delays has been analyzed by Huang and Vandewalle [9]. The Galerkin method for elliptic, parabolic, and hyperbolic problems has been studied in [4,14,17], etc. In this paper, the optimal convergence order and the asymptotic stability of the Galerkin method applied to the process of semidiscretization with respect to the spatial variable x to (1.1)are investigated. Then it is discrete with respect to the time variable t by using forward Euler, backward Euler and Crank–Nicolson schemes, and the optimal convergence orders and asymptotic stability of these three full discretization schemes are studied respectively. Some notations: Denote · and · r as the norms in L2 = L2 () and the Sobolev space Hr = Hr () = W2r (), respectively, so that for real-valued function v,
v = vL2 :=
|v(x)|2 dx
1/2
,
and for a positive integer r ,
⎞1/2
di v(x) := ⎝ 2 ⎠ . dxi ⎛
vr = vHr
i≤r
Let v(x), w(x) (x ࢠ ) be real-valued functions. Throughout this paper,
(v, w) :=
v(x)w(x)dx,
∇ v, ∇ w :=
∂ v(x) ∂ w(x) dx, ∂x ∂x
and C will denote a positive constant not necessarily the same at different occurrences, which may depends on a, b and t of (1.1), but is independent of h and t (the stepsize in t-direction). Assumptions: In this paper, assume that u(t) := u(·, t), ut (t) := ut (·, t), utt (t) := utt (·, t), uttt (t) := uttt (·, t) ∈ H01 (), and at t = nτ (n = −1, 0, 1, . . . ), the derivative of t is denoted as the left derivative. Definition 1.1. The zero solution of (1.1) is called asymptotically stable if the solution u(x, t) of (1.1) corresponding to any sufficiently differential function φ (x, t) with φ (0, t) = φ (π , t) = 0 satisfies
lim u(x, t) = 0,
t→∞
x ∈ .
Theorem 1.2. [16] The zero solution of (1.1) is delay-independently asymptotically stable for a > b. 2. Semidiscretization of Galerkin methods For the preliminary step of defining a spatially semidiscrete approximation solution to (1.1), we thus write the problem in weak form: multiply the equation by a smooth function v(x) which vanishes on ࢚ (or v ∈ H01 ()), then integrate over , and apply Green’s formula to the second and third terms for the first equation of (1.1), yields
(ut (·, t), v) + a(∇ u(·, t), ∇ v) + b(∇ u(·, t − τ ), ∇ v) = 0, ∀v ∈ H01 (), t > 0.
(2.1)
Let N > 0 be a given positive integer number, h := πN be a stepsize in x-direction and
Ih := {xn := nh, n = 0, 1, . . . , N} be a given mesh on [0, π ]. Set In ࣎ [xn − 1 , xn ] (1 ࣘ n ࣘ N). Define
Shr := {χ ∈ C () : χ |In ∈ P (r−1)(In ), 1 ≤ n ≤ N,
χ (0) = χ (π ) = 0},
where P(r − 1) denotes the space of all (real) polynomials of degree not exceeding r − 1. We may then pose the approximate problem to find uh (t) ࣎ uh ( ·, t), belonging to Shr for each t, such that
(uh,t (t), χ) + a(∇ uh (t), ∇χ )) + b(∇ uh (t − τ ), ∇χ )) = 0, ∀χ ∈ Shr , t > 0, uh (x, t) = φh (x, t), −τ ≤ t ≤ 0,
(2.2)
where φ h ( ·, t) is some an approximation of φ ( ·, t) in Shr . N
In terms of the basis { j }1 h for Shr , our semidiscrete problem may be stated: find the coefficients α j (t) in uh (x, t) := such that Nh
j=1
αj (t)( j , k )) + a
Nh
j=1
αj (t)(∇ j , ∇ k )) + b
Nh
j=1
αj (t − τ )(∇ j , ∇ k )) = 0,
Nh j=1
αj (t) j (x)
162
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
k = 1, . . . , Nh , t > 0, and α j (t) = γ j (t) for j = 1, . . . , Nh and − τ ࣘ t ࣘ 0, with γ j (t) being the components of the given initial approximation φ h (x, t). In matrix notation this may be expressed as
Aα (t) + aBα(t) + bBα(t − τ ) = 0, for t > 0, α(t) = γ (t), for − τ ≤ t ≤ 0,
where A = (ajk ) is the mass matrix with elements ajk = ( j , k ), B = (bjk ) is the stiffness matrix with bjk = (࢟ j , ࢟ k ), α (t) and γ (t) are the vectors of α j and γ j respectively. The dimension of all these items equals Nh , the dimension of Shr . Since like the stiffness matrix B the mass matrix A is a Gram matrix, and thus in particular positive definite and invertible (see [15]), the above system of ordinary differential equations may be written as
α (t) + aA−1 Bα(t) + bA−1 Bα(t − τ ) = 0, for t > 0, α(t) = γ (t), for − τ ≤ t ≤ 0,
(2.3)
and hence obviously has a unique solution for positive t. 2.1. The convergence analysis Introduce the elliptic or Ritz projection Rh onto Shr as the orthogonal projection with respect to the inner product (࢟v, ࢟w), so that
(∇ Rh v, ∇χ) = (∇ v, ∇χ ), ∀χ ∈ Shr , for v ∈ H01 ().
(2.4)
Lemma 2.1. [15] Assume that for any v ∈ Hs () ∩ H01 (),
inf {v − χ + h∇(v − χ )} ≤ Chs vs ,
χ ∈Shr
for 1 ≤ s ≤ r
holds. Then, with Rh defined by (2.4) we have
Rh v − v + h∇(Rh v − v) ≤ Chs vs , for any v ∈ Hs () ∩ H01 (), 1 ≤ s ≤ r. Let u(t) ࣎ u( ·, t) and u : [0, +∞) → H01 (). Define Dh : H01 () → Shr by
a(∇ Dh u(t) − ∇ u(t), ∇χ ) + b(∇ Dh u(t − τ ) − ∇ u(t − τ ), ∇χ ) = 0, Dh u(t) = Rh u(t) = Rh φ(t),
for t > 0, χ ∈ Shr ,
for − τ ≤ t ≤ 0.
Theorem 2.2. Let u be the solution of (2.1). Assume that u(t) − Rh u(t)−l ࣘ that Dh u(t) ∈ Shr is the above projection of u(t), we have
(2.6) Chr + l u(t)r ,
− 1 ࣘ l ࣘ r − 2, −τ ࣘ t ࣘ 0 holds and
Dh u(t) − u(t)−l ≤ Chr+l , −1 ≤ l ≤ r − 2, t > 0. Proof. Denote
ρ u − Dh u = (u − Rh u) + (Rh u − Dh u) ξ + η, so for any χ ∈ Shr ,
a(∇η(t), ∇χ) = a(∇ Rh u(t) − ∇ Dh u(t), ∇χ ) = a(∇ u(t) − ∇ Dh u(t), ∇χ ) = b(∇ Dh u(t − τ ) − ∇ u(t − τ ), ∇χ ) = b(∇ Dh u(t − τ ) − ∇ Rh u(t − τ ) + ∇ Rh u(t − τ ) − ∇ u(t − τ ), ∇χ ) = −b(∇η(t − τ ) + ∇ξ (t − τ ), ∇χ ). Without loss of generality, let t ∈ ((k − 1)τ , kτ ], k ∈ N. Taking χ = η(t), then
η(t)1 ≤ C (η(t − τ )1 + ξ (t − τ )1 ) ≤ C (η(t − 2τ )1 + ξ (t − τ )1 + ξ (t − 2τ )1 ) ≤ · · · ≤ C (η(t − kτ )1 +
k
ξ (t − iτ )1 ),
i=1
by (2.6),
η(t − kτ ) = Rh u(t − kτ ) − Dh u(t − kτ ) = 0, therefore by using Lemma 2.1 with s = r,
η(t)1 ≤ C
k
i=1
ξ (t − iτ )1 = C
k
i=1
u(t − iτ ) − Rh u(t − iτ )1
(2.5)
(2.7)
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
≤ Chr−1
k
163
u(t − iτ )r .
i=1
By the definition of ρ and Lemma 2.1,
ρ(t)1 ≤ ξ (t)1 + η(t)1 ≤ Chr−1
k
u(t − iτ )r .
(2.8)
i=0
Thus the estimate in H1 () follows. To get other estimates, construct an auxiliary problem −w = ϕ in , with ϕ ∈ H01 (). Define wh = Rh w. By (2.5) and Lemma 2.1, there is
(ρ(t), ϕ) = (ρ(t), −w) = (∇ρ(t), ∇ w) = (∇ρ(t), ∇ w(t) − ∇ wh ) + (∇ρ(t), ∇ wh ) =
b a
b a
(∇ρ(t), ∇ w − ∇ wh ) − (∇ρ(t − τ ), ∇ wh − ∇ w) − (∇ρ(t − τ ), ∇ w)
≤ C ((ρ(t)1 + ρ(t − τ )1 )wh − w1 + ρ(t − τ )−l wl+2 ) ≤ C (hl+1 (ρ(t)1 + ρ(t − τ )1 ) + ρ(t − τ )−l )ϕl By using the definition of the norm ρ−l ,
ρ(t)−l ≤ C (hl+1 (ρ(t)1 + ρ(t − τ )1 ) + ρ(t − τ )−l ) ≤ C (hl+1 (ρ(t)1 + ρ(t − τ )1 + ρ(t − 2τ )1 ) + ρ(t − 2τ )−l ) ⎛
≤ · · · ≤ C ⎝h
l+1
k
⎞
ρ(t − iτ )1 + ρ(t − kτ )−l ⎠ ,
i=1
moreover, by the assumption for − τ ࣘ t − kτ ࣘ 0,
ρ(t − kτ )−l = u(t − kτ ) − Dh u(t − kτ )−l = u(t − kτ ) − Rh u(t − kτ )−l ≤ Chr+l u(t − kτ )r , which together with (2.8) complete the proof. The following theorem can be proved by a similar way to Theorem 2.2. Theorem 2.3. Let u be the solution of (2.1). Assume that ut (t) − Rh ut (t)i = ut (t) − (Rh u)t (t)i ࣘ Chr + i ut (t)r , i = 1, 2, −τ ࣘ t ࣘ 0 holds and ρ (t) ࣎ u(t) − Dh u(t). Then
ρt (t) + hρt (t)1 ≤ Chr , for t > 0.
(2.9)
Theorem 2.4. Let u and uh be the solutions of (2.1) and (2.2) respectively. Assume that u(t) − Rh u(t)i ࣘ Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Chr + i u(t)r ,
ut (t) −
uh (t) − u(t) ≤ Chr , for t > 0. Proof. Denote
uh (t) − u(t) = θ (t) + ρ(t), where θ (t) = uh (t) − Dh u(t),
ρ(t) = Dh u(t) − u(t).
(2.10)
The second term can be easily bounded by Theorem 2.2 and obvious there exists the following estimate:
ρ(t) ≤ Chr .
(2.11)
In order to estimate θ (t), note that by (2.5) and our definitions,
(θt (t), χ) + a(∇θ (t), ∇χ ) + b(∇θ (t − τ ), ∇χ ) = (uh,t , χ) − (Dh ut (t), χ ) + a(∇ uh (t), ∇χ ) − a(∇ Dh u(t), ∇χ ) + b(∇ uh (t − τ ), ∇χ ) − b(∇ Dh u(t − τ ), ∇χ ) = (uh,t , χ) − (Dh ut (t), χ ) + a(∇ uh (t), ∇χ ) − a(∇ u(t), ∇χ ) + b(∇ uh (t − τ ), ∇χ ) − b(∇ u(t − τ ), ∇χ ) = (ut (t) − Dh ut (t), χ ) = −(ρt (t), χ ), for ∀χ ∈ Shr . (2.12) Taking χ = θ (t),
1 d θ (t)2 + aθ (t)21 ≤ C (θ (t − τ )1 θ (t)1 + ρt (t)θ (t)). 2 dt So
a 1 d θ (t)2 + θ (t)21 ≤ C (θ (t − τ )21 + ρt (t)θ (t)). 2 dt 2
(2.13)
164
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
Without loss of generality, assume that t ∈ ((k − 1)τ , kτ ], k ∈ N. Integrating both sides of the above inequality from (k − 1)τ to t yields
θ (t)2 +
t
(k−1)τ
θ (s)21 ds ≤ C (θ ((k − 1)τ )2 +
Note that
t
(k−1)τ
θ (s)21 ds ≤ C (θ ((k − 1)τ )2 + = C (θ ((k − 1)τ )2 +
t
(k−1)τ t−τ
(k−2)τ
t
(k−1)τ
θ (s − τ )21 ds +
θ (s)21 ds +
t
ρt (s)θ (s)ds +
(k−1)τ
⎛ k−1
≤ C⎝ θ (iτ )2 + i=0
therefore
θ (t)
⎛ 2
≤ C⎝
k−1
θ (iτ ) +
i=0
So
⎛
θ (kτ ) ≤ C ⎝ 2
2
k−1
t−kτ
−τ
θ (iτ ) + 2
i=0
t−kτ −τ
t−τ
(k−2)τ
−τ
θ (s)
θ (s)
t
(k−3)τ
k−1
i=0
2 1 ds
+
k−1
2 1 ds
+
t−(k−i−1)τ iτ
(k−1)τ
ρt (s)θ (s)ds).
ρt (s)θ (s)ds)
ρt (s)θ (s)ds)
θ (s)21 ds
t−(k−i−1)τ
iτ
⎞
ρt (s)θ (s)ds⎠ , ⎞
ρt (s)θ (s)ds⎠ . ⎞
k−1 (i+1)τ
i=0
(k−1)τ
t
ρt (s)θ (s)ds) ≤ · · ·
θ (s)21 ds +
i=0
0
t
(k−1)τ t−2τ
≤ C (θ ((k − 1)τ )2 + θ ((k − 2)τ )2 + +
θ (s − τ )21 ds +
ρt (s)θ (s)ds⎠ ,
iτ
by the discrete Gronwall inequality (see [3,5]),
θ (kτ )2 ≤ C (θ (0)2 + therefore,
θ (t)2 ≤ C (θ (0)2 +
0 −τ
0
−τ
θ (s)21 ds +
θ (s)21 ds +
k−1 (i+1)τ
i=0
t 0
iτ
ρt (s)θ (s)ds),
ρt (s)2 ds +
0
t
θ (s)2 ds),
then θ (t) ࣘ Chr follows by Gronwall inequality (see [3]) and the assumption of the theorem, which together with (2.11) complete the proof. 2.2. The stability analysis Definition 2.5. If the solution uh (x, t) of (2.2) corresponding to any sufficiently differentiable function φ h (x, t) with φ h (0, t) = φ h (π , t) = 0 satisfies
lim uh (x, t) = 0,
t→∞
x ∈ [0, π ],
then the zero solution of (2.2) is called asymptotically stable. Theorem 2.6. The zero solution of (2.2) is asymptotically stable for a > b. Proof. Let α (t) = eλt C1 , where C1 is a constant vector. The characteristic equation of (2.3) is:
|λI + aA−1 B + bA−1 Be−λτ | = 0,
(2.14)
which is equivalent to
λ = −(a + be−λτ )λA−1 B , where λA−1 B denotes the corresponding eigenvalue of A−1 B. Moreover, because A and B are positive definite, so λA−1 B > 0. In fact, since B is a positive definite matrix, so there exists an invertible matrix C, such that B = CT C, and
CA−1 BC −1 = CA−1 C T CC −1 = CA−1 C T , which is positive definite and similar to A−1 B.
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
165
Let λ = λ1 + iλ2 , then
λ1 = −(a + be−λ1 τ cos λ2 τ )λA−1 B , λ2 = be−λ1 τ sin λ2 τ λA−1 B .
If λ1 ࣙ 0, then a + be−λ1 τ cos λ2 τ ≥ a − b and λ1 = −(a + be−λ1 τ cos λ2 τ )λA−1 B ≤ −(a − b)λA−1 B < 0, which is a contradiction. So if a > b, then λ1 < 0, i.e., the zero solution of (2.2) is asymptotically stable. 3. Fully discretization for Galerkin methods In this section, we shall focus our attention on numerical stability of backward Euler, Crank–Nicolson and forward Euler– Galerkin methods applied to (2.2). 3.1. The backward Euler–Galerkin method Let t = τ /m be a given stepsize with integer m ࣙ 1, the gridpoints tn be defined by tn = nt (n = 0, 1, 2, . . .), and Un be the approximation in Shr of u(t) at t = tn = nt, then the backward Euler–Galerkin method to (2.2) leads to a numerical process of the following type
U n − U n−1 ,χ t
+ a(∇ U n , ∇χ ) + b(∇ U n−m , ∇χ ) = 0,
for n > 0,
χ ∈ Shr ,
(3.1)
where Un ( · ) ࣎ φ h ( ·, tn ) for − m ࣘ n ࣘ 0. Nh Let U n (x) := i (x)αin . Then (3.1) can be reduced to i=1
Nh Nh Nh
1
1 n αi ( i , j ) − αin−1 ( i , j ) = −a αin (∇ i , ∇ j ) − b αin−m (∇ i , ∇ j ). t t Nh
i=1
i=1
i=1
(3.2)
i=1
With the notation as in the semidiscrete situation, this may be written as
(A + atB)α n = Aα n−1 − btBα n−m , for n > 0, αn = γ n, for − m ≤ n ≤ 0,
with γ n being the corresponding vector of components of the given initial approximation φ h (tn ), and α n := (α1 , . . . , αNh )T , and hence obviously has a unique solution for positive n. 3.1.1. The convergence analysis We shall prove the following error estimate: Theorem 3.1. Let u and Un be the solutions of (2.1) and (3.1), respectively. Assume that u(t) − Rh u(t)i ࣘ Chr + i u(t)r , ut (t) − Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Un − u(tn ) ≤ C (hr + t), for n = 1, 2, . . . . Proof. Denote
U n − u(tn ) = (U n − Dh u(tn )) + (Dh u(tn ) − u(tn )) = θ n + ρ n , and ρ n = ρ (tn ) is bounded as claimed in (2.11). This time, a calculation corresponding to (2.12) yields
θ n − θ n−1 t
,χ
+ a(∇θ n , ∇χ ) + b(∇θ n−m , ∇χ ) = −(wn , χ ),
∀χ ∈ Shr ,
where
Dh u(tn ) − Dh u(tn−1 ) − ut (tn ) = Dh − I ∂¯ u(tn ) + (∂¯ u(tn ) − ut (tn )) t u(tn ) − u(tn−1 ) = : wn1 + wn2 , ∂¯ u(tn ) := . t
wn =
Choosing χ = θ n , yields
θ n − θ n−1 t
, θ n + aθ n 21 + b(∇θ n−m , ∇θ n ) = −(wn , θ n ).
By using Schwarz inequality,
θ n − θ n−1 t
, θ n + θ n 21 ≤ C θ n−m 21 + wn θ n .
(3.3)
166
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
So
θ n 2 + tθ n 21 ≤ C (θ n−1 2 + tθ n−m 21 + (t)2 wn 2 ). Without loss of generality, assume that n ∈ ((k − 1)m, km], k ∈ N. Then
tθ n 21 ≤ C (θ n−1 2 + tθ n−m 21 + (t)2 wn 2 ) ≤ C (θ n−1 2 + θ n−m−1 2 + tθ n−2m 21 + (t)2 wn 2 + (t)2 wn−m 2 ) ⎞ ⎛ k−1 k−1
≤ ··· ≤ C⎝ θ n−im−1 2 + tθ n−km 21 + (t)2 wn−im 2 ⎠ . i=0
So
⎛
θ n 2 ≤ C ⎝
k−1
i=0
θ n−im−1 2 + tθ n−km 21 + (t)2
i=0
k−1
⎞
wn−im 2 ⎠ ,
i=0
therefore by using the discrete Gronwall inequality (see [3,5]),
⎛
θ ≤ C ⎝θ + tθ n 2
0 2
+ (t)
n−km 2 1
2
k−1
w
⎞
n−im 2 ⎠
.
(3.4)
i=1
Write
wn1 = (Dh − I)∂¯ u(tn ) = t−1
tn tn−1
d/dt (Dh − I)u(t)dt,
so
(t)2
k−1 k−1
n−im 2 = w1 i=0
i=0
≤
tn−im tn−im−1
k−1
tn−im tn−im−1
i=0
d/dt (Dh − I)u(t)dt Chr ut (t)r dt
2
2
2r
≤ Ch .
(3.5)
Further
twi2 = u(ti ) − u(ti−1 ) − tut (ti ) = −
ti ti−1
(t − ti−1 )utt (t)dt,
so that
(t)2
k−1
wn−im 2 = 2
i=0
k−1
i=0
≤ (t)2
tn−im tn−im−1
2
(t − tn−im−1 )utt (t)dt
k−1 tn−im
i=0
tn−im−1
2 utt (t)dt
≤ C (t) , 2
(3.6)
which together with (3.4), (3.5) and the assumption of the theorem complete the proof. 3.1.2. The stability analysis Definition 3.2. If the solution Un of (3.1) corresponding to any sufficiently differentiable function φ h (x, t) with φ h (0, t) = φ h (π , t) = 0 satisfies
lim U n = 0,
n→∞
x ∈ [0, 1],
then the zero solution of (3.1) is called asymptotically stable. Let K ࣎ [xi , xi + 1 ] be an element of the finite element, and Kˆ := [−1, 1] be the reference element in ξ -plane. Then
K
i j dx =
h 2
Kˆ
ˆ i ˆ j dξ ,
K
∇ i ∇ j dx =
2 h
Kˆ
ˆ i∇ ˆ j dξ . ∇
By (3.2),
αn =
h 2
−1 −1 2bt h h 2at 2at Bˆα n−m Aˆα n−1 − Aˆ + Bˆ Aˆ + Bˆ 2 h h 2 h
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
4at −1 −1 n−1 4bt 4at −1 −1 −1 n−m ˆ Bˆ ˆ Bˆ I+ α − , A A Aˆ Bˆα I + h2 h2 h2
=
where Aˆ and Bˆ is the corresponding mass and stiffness matrices. Let α n = λn C1 . The characteristic equation of (3.7) is:
λm − 1 +
4at λ ˆ−1 ˆ h2 A B
−1
λm−1 +
4bt h2
1+
4at λ ˆ−1 ˆ h2 A B
−1
167
(3.7)
λAˆ−1 Bˆ = 0,
(3.8)
ˆ Therefore the backward Euler scheme is delay-independently where λAˆ−1 Bˆ denotes the corresponding eigenvalue of Aˆ−1 B. asymptotically stable if and only if (3.8) is Schur polynomial for all m ࣙ 1(that is, the modulus of all zeros of the characteristic polynomial is less than 1). In order to prove that the characteristic polynomial is a Schur polynomial, the following lemma is needed. Lemma 3.3. [13] Let γ m (z) = α (z)zm − β (z) be a polynomial, where α (z) and β (z) are polynomial of constant order. Then the polynomial γ m (z) is a Schur polynomial for any m ࣙ 1 if and only if the following conditions hold (a) α (z) = 0⇒|z| < 1, (b) |β(z)| ≤ |α(z)|, ∀z ∈ C, |z| = 1 and (c) γm (z) = 0, ∀z ∈ C, |z| = 1, ∀m ≥ 1. Theorem 3.4. Suppose that 0 ࣘ b < a. Then the zero solution of the backward Euler scheme is delay-independently asymptotically stable. Proof. Denote
4at λ ˆ−1 ˆ h2 A B
α(λ) = λ − 1 + and
β(λ) = −
4bt h2
−1 4at λ λAˆ−1 Bˆ . 1+ ˆ−1 ˆ h2 A B
(a) If α (λ) = 0, then |λ| = (1 + (b) For ∀λ ∈ C, |λ| = 1,
|α(λ)|
−1
4at λAˆ−1 Bˆ )−1 h2
< 1.
−1 4at λAˆ−1 Bˆ 4at h2 = λ − 1 + λ ˆ−1 ˆ = λ − 1 + h2 A B 1 + 4at λAˆ−1 Bˆ h2 ≥
λ
1
4at Aˆ−1 Bˆ h2 + 4at Aˆ−1 Bˆ h2
λ
>
λ
1
4bt Aˆ−1 Bˆ h2 + 4at Aˆ−1 Bˆ h2
λ
= |β(λ)|.
(c) By (b), it is trivial. The proof is completed. 3.2. The Crank––Nicolson–Galerkin method The Crank––Nicolson–Galerkin method to (2.2) corresponding to (3.1) for the backward Euler methods leads to a numerical process of the following type
U n − U n−1 ,χ t
+a
∇ Un + ∇ Un−1 2
, ∇χ
+b
∇ Un−m + ∇ Un−m−1 2
, ∇χ
= 0,
(3.9)
for n > 0, where Un ( · ) = φ h ( ·, tn ) for − m ࣘ n ࣘ 0. Nh Let U n (x) := i (x)αin . Then (3.9) can be reduced to i=1
Nh 1
1 n αi ( i , j ) − αin−1 ( i , j ) t t Nh
i=1
i=1
Nh a n b
=− αi + αin−1 (∇ i , ∇ j ) − αin−m + αin−m−1 (∇ i , ∇ j ) for n > 0. 2 2 Nh
i=1
i=1
With the notation as in the semidiscrete situation, this may be written as
A + 2a tB α n = A − 2a tB α n−1 − t 2b Bα n−m − t 2b Bα n−m−1 , n > 0,
αn = γ n,
and hence obviously has a unique solution for positive n.
(3.10)
168
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
3.2.1. The convergence analysis We shall prove the following error estimate: Theorem 3.5. Let u and Un be the solutions of (2.1) and (3.9), respectively. Assume that u(t) − Rh u(t)i ࣘ Chr + i u(t)r , ut (t) − Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Un − u(tn ) ≤ C (hr + (t)2 ), for n = 1, 2, . . . . Proof. Denote
U n − u(tn ) = (U n − Dh u(tn )) + (Dh u(tn ) − u(tn )) = θ n + ρ n , and ρ n = ρ (tn ) is bounded as claimed in (2.11). This time, a calculation corresponding to (2.12) yields
θ n − θ n−1 t
,χ
+a
∇θ n + ∇θ n−1 2
, ∇χ
+b
∇θ n−m + ∇θ n−m−1 2
, ∇χ
= −(wn , χ), ∀χ ∈ Shr ,
(3.11)
where
Dh u(tn ) − Dh u(tn−1 ) ut (tn ) + ut (tn−1 ) − t 2 u ( t t n ) + ut (tn−1 ) = (Dh − I)∂¯ u(tn ) + ∂¯ u(tn ) − =: wn1 + wn2 . 2
wn =
Choosing χ =
θ n + θ n−1
, yields
2
n n−m θ + θ n−1 2 + ∇θ n−m−1 + b ∇θ , + a t 2 2 2 1 θ n + θ n−1 = − wn , . 2
θ n − θ n−1 θ n + θ n−1
,
By using Schwarz inequality,
θ n − θ n−1 θ n + θ n−1 t
So
,
2
∇θ n + ∇θ n−1
2
n n n−1 n−m + θ n−m−1 2 θ + θ n−1 2 ≤ C θ + w n θ + θ . + 2 2 2 1 1
n n−m 2 2 θ + θ n−1 2 θ + θ n−m−1 n−1 2 n 2 ≤C θ + t + ( t ) w . 2 2
θ n 2 + t
1
1
Without loss of generality, assume that n ∈ ((k − 1)m, km], k ∈ N. Then
2 n n−m θ + θ n−1 2 θ + θ n−m−1 n−1 2 2 n 2 t + t ≤ C θ + (t) w 2 2 1 1 n−2m 2 θ + θ n−2m−1 + (t)2 (wn 2 + wn−m 2 ) ≤ C θ n−1 2 + θ n−m−1 2 + t 2 1 ⎛ ⎞ k−1 k−1 θ n−km + θ n−km−1 2
≤ ··· ≤ C⎝ θ n−im−1 2 + t wn−im 2 ⎠ . + (t)2 2 i=0
Therefore
⎛
θ n 2 ≤ C ⎝
k−1
i=0
1
i=0
⎞ k−1 θ n−km + θ n−km−1 2
θ n−im−1 2 + t wn−im 2 ⎠ . + (t)2 2 i=0
1
By the assumption of the theorem and by using the discrete Gronwall inequality (see [3,5]),
⎛
θ ≤ C ⎝θ n 2
0 2
⎞ k−1 θ n−km + θ n−km−1 2
2 n−im 2 ⎠ + t w . + (t) 2 1
Write
wn1 = (Dh − I)∂¯ u(tn ) = t−1
tn tn−1
(Dh − I)ut (t)dt,
i=0
(3.12)
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
169
so
(t)2
k−1 k−1
n−im 2
w1 ≤ i=1
i=1
tn−im tn−im−1
Chr ut (t)r dt
2 ≤ Ch2r .
(3.13)
Further
ti ut (ti ) + ut (ti−1 ) 2 ( t ) − u ( t ) − t ( t ) u ≤ C uttt (t) dt, twi2 = i−1 i 2 ti−1
so that
(t)2
k−1 k−1
n−im 2 w2 ≤ C (t)4 i=1
tn−im
tn−im−1
i=1
uttt (s) dt
2
≤ C (t)4 ,
(3.14)
which together with (3.12), (3.13) and the assumption of the theorem complete the proof. 3.2.2. The stability analysis By (3.10),
−1 −1 h bt h h at at at Bˆ(α n−m + α n−m−1 ) Aˆ + Bˆ Aˆ − Bˆ α n−1 − Aˆ + Bˆ 2 h 2 h h 2 h 2at −1 −1 2bt 2at −1 2at −1 −1 −1 ˆ Bˆ ˆ Bˆ α n−1 − ˆ Bˆ = I+ A A A Aˆ Bˆ(α n−m + α n−m−1 ). I − I + h2 h2 h2 h2
αn =
(3.15)
Let α n = λn C1 . The characteristic equation of (3.15) is:
λm −
1− 1+
2at h2 2at h2
2bt λAˆ−1 Bˆ m−1 λAˆ−1 Bˆ h2 λ − (λ + 1) = 0. 2at λAˆ−1 Bˆ 1 + h2 λAˆ−1 Bˆ
(3.16)
Theorem 3.6. Suppose that 0 ࣘ b < a. Then the zero solution of the Crank–Nicolson scheme is delay-independently asymptotically stable. Proof. Denote α(λ) = λ −
1− 2at λAˆ−1 Bˆ 2 h 1+ 2at λAˆ−1 Bˆ h2
(a) If α (λ) = 0, then |λ| =
and β(λ) =
1− 2at λAˆ−1 Bˆ 2 h
1+ 2at λAˆ−1 Bˆ 2
2bt λAˆ−1 Bˆ h2 1+ 2at λAˆ−1 Bˆ 2 h
λ+1 .
< 1.
h
(b) For ∀λ ∈ C, |λ| = 1, denote λ = cos θ + isin θ , we have
2i sin θ λ − 1 cos θ − 1 + i sin θ = = . λ + 1 cos θ + 1 + i sin θ 2 + 2 cos θ So we have
2at λAˆ−1 Bˆ (λ − 1 ) −1 B 2 λA ˆ ˆ h = +
2at 2at λ+1 1 + h2 λAˆ−1 Bˆ λ + 1 1 + 2 λAˆ−1 Bˆ h 2at 2bt β(λ) λAˆ−1 Bˆ λAˆ−1 Bˆ h2 h2 . > = ≥ 2at 2at λ + 1 1 + 2 λ ˆ−1 ˆ 1 + 2 λ ˆ−1 ˆ
α(λ) λ − λ + 1 =
1− 2at 2 λAˆ−1 Bˆ h 1+ 2at h2
h
A
B
h
A
B
(c) By (b), it is trivial. The proof is completed. 3.3. The forward Euler–Galerkin method The forward Euler–Galerkin method to (2.2) corresponding to (3.1) for the backward Euler method leads to a numerical process of the following type
U n − U n−1 ,χ t
+ a(∇ U n−1 ,
∇χ ) + b(∇ Un−m−1 , ∇χ ) = 0, for n > 0.
where Un ( · ) = vh ( ·, tn ) for − m ࣘ n ࣘ 0.
(3.17)
170
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π . Fig. 1. The semidiscretization Galerkin method with h = 10
π . Fig. 2. The semidiscretization Galerkin method with h = 20
Let U n (x) :=
Nh i=1
i (x)αin . Then (3.1) can be reduced to
Nh Nh Nh Nh
1
1
αin ( i , j ) − αin−1 ( i , j ) = −a αin−1 (∇ i , ∇ j ) − b αin−m−1 (∇ i , ∇ j ). t t i=1
i=1
i=1
(3.18)
i=1
With the notation as in the semidiscrete situation, this may be written as
n Aα = (A − atB)α n−1 − btBα n−m−1 , for n > 0, αn = γ n, for − m ≤ n ≤ 0, and hence obviously has a unique solution for positive n.
3.3.1. The convergence analysis By using the same method of Section 3.1.1 as the proof of the convergence for the backward Euler–Galerkin method, the following theorem can be easily obtained. Theorem 3.7. Let u and Un be the solutions of (2.1) and (3.1). Assume that u(t) − Rh u(t)i ࣘ Chr + i u(t)r , ut (t) − Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Un − u(tn ) ≤ C (hr + t), for n = 1, 2, · · · .
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π . Fig. 3. The semidiscretization Galerkin method with h = 40
π , t = τ . Fig. 4. The backward Euler–Galerkin method with h = 10 100
π , t = τ . Fig. 5. The backward Euler–Galerkin method with h = 20 400
171
172
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = Fig. 6. The backward Euler–Galerkin method with h = 40
τ
1600
.
π , t = τ . Fig. 7. The Crank–Nicolson–Galerkin method with h = 10 100
3.3.2. The stability analysis By (3.18),
2 −1 h 4bt −1 n−m 2at Aˆ Bˆα Aˆ Aˆ − Bˆ α n−1 − h 2 h h2 4at −1 4bt −1 n−m ˆ Bˆ α n−1 − = I− . A Aˆ Bˆα h2 h2
αn =
(3.19)
Let α n = λn C1 . The characteristic equation of (3.7) is:
λm − 1 −
4at λ ˆ−1 ˆ h2 A B
λm−1 +
4bt λ ˆ−1 ˆ = 0. h2 A B
(3.20)
Theorem 3.8. Suppose that 0 ࣘ b < a. Then the zero solution of the backward Euler scheme is delay-independently asymptotically stable if and only if t ≤ 4aλ1 . h2 Aˆ−1 Bˆ
Proof. Denote α(λ) = λ − (1 −
4at λAˆ−1 Bˆ ) and h2
(a) If α (λ) = 0, then |λ| = 1 −
4at λAˆ−1 Bˆ h2
β(λ) = − 4bt λAˆ−1 Bˆ . h2
< 1 if and only if
t h2
1 . 2aλAˆ−1 Bˆ 1 if t ≤ , 4aλAˆ−1 Bˆ h2
<
(b) For ∀λ ∈ C, |λ| = 1, we assume that λ ࣎ cos θ + isin θ , then
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
173
π , t = τ . Fig. 8. The Crank–Nicolson–Galerkin method with h = 20 200
π , t = τ . Fig. 9. The Crank–Nicolson–Galerkin method with h = 40 400
Table 1 The 2-norm of the errors for the semidiscrezation Galerkin method at t = 1.
N
m
u(tn ) − Un (h = π /N)
u(tn ) − Un (h = π /N) u(tn ) − Un (h = π /2N)
10 20 40
100 400 1600
4.4903e-06 2.8574e-07 1.7939e-08
– 1.5715e+01 1.5928e+01
4at
2
|α(λ)|2 = λ − 1 + 2 λAˆ−1 Bˆ h
2 4at 4at = 2(1 − cos θ ) 1 − λ λ + ˆ−1 ˆ ˆ−1 ˆ h2 A B h2 A B 2 2 4at 4bt ≥ λ ˆ−1 ˆ > λ ˆ−1 ˆ = |β(λ)|2 ; h2 A B h2 A B
Otherwise, if
t h2
>
1
4aλAˆ−1 Bˆ
, then
|α(λ)|2
<
( 4t λ )2 a2 . h2 Aˆ−1 Bˆ
Choose
b 4t λ = |β(λ)|. h2 Aˆ−1 Bˆ Therefore,
t h2
≤
1 4aλAˆ−1 Bˆ
is a sufficient and necessary condition.
b2
=
|α(λ)|2 a2 + 4t ( 2 λ ˆ−1 ˆ )2 h
A
2
B
, so
|α(λ)|
4t λ −1 h2 Aˆ Bˆ
< b < a and |α(λ)| <
174
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = τ . Fig. 10. The forward Euler–Galerkin method with h = 10 100
π , t = τ . Fig. 11. The forward Euler–Galerkin method with h = 20 400
Table 2 The 2-norm of the errors for the backward Euler–Galerkin method at t = 1. N
m
u(tn ) − Un
u(tn ) − Un (h = π /N, t = τ /m) u(tn ) − Un (h = π /2N, t = τ /4m)
10 20 40
100 400 1600
4.1615e–06 2.5857e–07 1.6136e–08
– 1.6094e+01 1.6024e+01
Table 3 The 2-norm of the errors for the Crank–Nicolson–Galerkin method at t = 1. N
m
u(tn ) − Un (h = π /N, t = τ /m)
u(tn ) − Un (h = π /N, t = τ /m) u(tn ) − Un (h = π /2N, t = τ /2m)
10 20 40
100 200 400
4.5255e–06 2.8796e–07 1.8078e–08
– 1.5716e+01 1.5929e+01
(c) By (b), it is trivial. The proof is completed. Remark 3.9. In Theorem 3.8, λAˆ−1 Bˆ is an implicit expression, which will change as r (the degree of the polynomial) and Nh change. We list the following three cases for r = 2: • For Nh = 10, max λAˆ−1 Bˆ ≈ 2.7900e + 00;
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = Fig. 12. The forward Euler–Galerkin method with h = 40
τ
1600
175
.
π . Fig. 13. The semidiscrezation Galerkin method with h = 10
Table 4 The 2-norm of the errors for the forward Euler–Galerkin method at t = 1. N
m
u(tn ) − Un (h = π /N, t = τ /m)
u(tn ) − Un (h = π /N, t = τ /m) u(tn ) − Un (h = π /2N, t = τ /4m)
10 20 40
100 400 1600
3.2208e–05 2.0297e–06 1.2712e–07
– 1.5868e+01 1.5967e+01
• For Nh = 20, max λAˆ−1 Bˆ ≈ 2.9453e + 00; • For Nh = 40, max λAˆ−1 Bˆ ≈ 2.9862e + 00. 4. Numerical experiments In this section, the same parameters as [16] are used: φ (x, t) = sin x, τ = 1, a = 1.5, b = 1. We take r = 2. In Figs. 1–12, the errors of the semidiscretization, backward Euler, Crank–Nicolson and forward Euler Galerkin methods with the stepsize π , h = π and h = π for the semidiscretization Galerkin method, h = π , t = τ ; h = π , t = τ ; h = π , t = τ for h = 10 20 40 10 100 20 400 40 1600 π , t = τ ; h = π , t = τ ; h = π , t = τ for the Crank–– the backward and forward Euler–Galerkin methods and h = 10 100 20 200 40 400 Nicolson–Galerkin method are drawn. From these figures, it can be observed that the convergence orders obey the theoretical analysis. In Tables 1–4, the norms of the errors for the semidiscretization Galerkin and the three full discrete Galerkin methods at t = 1, u(1)−u (1)(h=π /N) )−Un (h=π /N,t=τ /m) and the ratios of u(t)−u h(t)(h=π /2N) for the semidiscretization Galerkin methods, uu(t(tn)−U n (h=π /2N,t=τ /4m) for the backward and h
n
176
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = τ . Fig. 14. The backward Euler–Galerkin method with h = 10 40
π , t = τ . Fig. 15. The Crank–Nicolson–Galerkin method with h = 10 40
)−Un (h=π /N,t=τ /m) forward Euler–Galerkin methods and uu(t(tn)−U n (h=π /2N,t=τ /2m) for the Crank––Nicolson–Galerkin method are given, respectively. n It can be seen that the convergence orders obey our theoretical analysis. π , t = τ , τ respectively for the forward Euler–Galerkin method; the In Figs. 13– 16, the same stepsize as [16] is taken: h = 10 40 50 π for the semidiscrezation Galerkin methods and h = π , t = τ for backward Euler and the Crank––Nicolson–Galerkin h = 10 10 40 methods. From Fig. 13, it can be seen that the semidiscrezation Galerkin method is asymptotically stable. From Figs. 14 and 15, it can 1 are asymptotically stable. be observed that the backward Euler and the Crank–Nicolson Galerkin methods with stepsize t = 40 1 1 From Figs. 16 and 17, it can be seen that for the two stepsizes t = 40 and t = 50 , the forward Euler methods are unstable in [16], the finite difference method with the first stepsize is unstable, but which with the second stepsize is asymptotically stable). 1 for the forward Euler–Galerkin method. It can be observed that the numerical solution with this In Fig. 18, we take t = 500 stepsize is asymptotically stable. Therefore, the mesh ratio for the forward Euler–Galerkin method is stricter than which for the forward Euler finite difference method. Conclusion. In practice, the stability analysis of a numerical method is very important. However, up to now, there are few papers published with respect to the numerical stability analysis for partial functional differential equations, no other than the numerical stability analysis of Galerkin methods for partial functional differential equations. In this paper, a complete analysis of the asymptotic stability of the semidiscrete (Galerkin methods) scheme and the full (forward, backward Euler and Crank– Nicolson methods) schemes for linear parabolic differential equations with delay are given: the semidiscrete, backward Euler and Crank–Nicolson full discrete schemes can unconditionally preserve the delay-independent asymptotic stability, but an additional restriction on time and spatial stepsize of the forward Euler full discrete scheme is needed to preserve the delay-independent asymptotic stability. These results for stability are consistent with the corresponding results of the finite difference methods. Moreover, numerical experiment has shown that in general, the addition restriction for the forward Euler–Galerkin full discrete
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = τ . Fig. 16. The forward Euler–Galerkin method with h = 10 40
π , t = τ . Fig. 17. The forward Euler–Galerkin method with h = 10 50
π , t = τ . Fig. 18. The forward Euler–Galerkin method with h = 10 500
177
178
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
scheme is stricter than which for the forward Euler finite difference method. In addition, it is proved that the convergence order of the semidiscrete scheme and the full discrete schemes are consistent with those of the original parabolic equation. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]
A. Ashyralyev, D. Agirseven, On convergence of difference schemes for delay parabolic equations, Comput. Math. Appl. 66 (2013) 1232–1244. A. Bellen, M. Zennaro, Numerical Methods for Delay Differential Equations, Clarendon Press, Oxford, 2003. H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, Cambridge University Press, Cambridge, 2004. H. Brunner, K. Mustapha, H. Mustapha, D. Schoetzau, An HP-version discontinuous Galerkin method for integro-differential equations of parabolic type, SIAM J. Numer. Anal. 49 (2011) 1369–1396. C.M. Chen, T. Shih, Finite Element Methods for Integrodifferential Equations, World Scientific, Singapore, 1998. D. Green, H.W. Stech, Diffusion and hereditary effects in a class of population models, in: S.N. Busenberg, K.L. Cooke (Eds.), Differential Equations and Applications in Ecology, Epidemics, and Population Problems, Academic Press, New York, 1981. W. Gu, P. Wang, A Crank–Nicolson difference scheme for solving a type of variable coefficient delay partial differential equations, J. Appl. Math. 2014 (2014) 1–6. P.J.V.d. Houwen, B.P. Sommeijer, C.T.H. Baker, On the stability of predictor-corrector methods for parabolic equations with delay, IMA J. Numer. Anal. 6 (1986) 1–23. C.M. Huang, S. Vandewalle, An analysis of delay-dependent stability for ordinary and partial differential equations with fixed and distributed delays, SIAM J. Sci.Comput. 25 (2004) 1608–1632. Y.F. Jin, J.X. Jiang, C.M. Hou, D.H. Guan, New difference scheme for general delay parabolic equations, J. Inform. Comput. Sci. 18 (2012) 5579–5586. H. Liang, M.Z. Liu, W.J. Lv, Stability of θ -schemes in the numerical solution of a partial differential equation with piecewise continuous arguments, Appl. Math. Lett. 23 (2010) 198–206. H. Liang, D.Y. Shi, W.J. Lv, Convergence and asymptotic stability of Galerkin methods for a partial differential equation with piecewise constant argument, Appl. Math. Comput. 217 (2010) 854–860. M.Z. Liu, M.N. Spijker, The stability of the θ -methods in the numerical solution of delay differential equations, IMA J. Numer. Anal. 10 (1990) 31–48. D.Y. Shi, W. Gong, J.C. Ren, A new stable second order nonconforming mixed finite element scheme for the stationary Stokes and Navier–Stokes equations, Math. Comput. Model. 53 (2011) 1956–1969. V. Thomée, Galerkin Finite Element Methods for Parabolic Problem, Springer, New York, 1997. H.J. Tian, Asymptotic stability of numerical methods for linear delay parabolic differential equations, Comput. Math.Appl. 56 (2008) 1758–1765. J.P. Wang, R. Zhang, Maximum principles for p1-conforming finite element approximations of quasi-linear second order elliptic equations, SIAM J. Numer. Anal. 50 (2012) 626–642. J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer-Verlag, New York, 1996. B. Zubik-Kowal, The method of lines for parabolic differential-functional equations, IMA J. Numer. Anal. 17 (1997) 103–123. B. Zubik-Kowal, Stability in the numerical solution of linear parabolic equations with delay term, BIT 41 (2001) 191–206.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Convergence and asymptotic stability of Galerkin methods for linear parabolic equations with delays ✩ Hui Liang∗ School of Mathematical Sciences, Heilongjiang University, Harbin, Heilongjiang 150080, China
a r t i c l e
i n f o
Keywords: Linear parabolic equations Delay Convergence Asymptotic stability Galerkin methods
a b s t r a c t This paper is concerned with the convergence and asymptotic stability of semidiscrete and full discrete schemes for linear parabolic equations with delay. These full discrete numerical processes include forward Euler, backward Euler and Crank–Nicolson schemes. The optimal convergence orders are consistent with those of the original parabolic equation. It is proved that the semidiscrete scheme, backward Euler and Crank–Nicolson full discrete schemes can unconditionally preserve the delay-independent asymptotic stability, but some additional restrictions on time and spatial stepsizes of the forward Euler full discrete scheme is needed to preserve the delay-independent asymptotic stability. Numerical experiments illustrate the theoretical results. © 2015 Elsevier Inc. All rights reserved.
1. Introduction Since there is sometimes a memory or delayed effect, partial functional differential equations are used frequently for describing and modeling mathematically many real life phenomena as demonstrated [18]. The following initial-boundary value problem for linear parabolic differential equations with a delay term:
⎧ ⎨ut (x, t) = auxx (x, t) + buxx (x, t − τ ), u(0, t) = u(π , t) = 0, ⎩ u(x, t) = φ(x, t),
t > 0, x ∈ := [0, π ], t ≥ −τ , −τ ≤ t ≤ 0, x ∈ ,
(1.1)
with a > 0, b ࣙ 0 and τ > 0 has been studied by Houwen et al. [8]. Uniqueness, local existence, and global continuation are straightforward by Green and Stech [6]. In particular, Tian [16] has analyzed the asymptotic stability of finite difference methods applied to (1.1). Applying the process of semidiscretization with respect to the spatial variable x to partial functional differential equations (1.1) results in a system of delay differential equations (DDEs). There are lots of references discussed the numerical methods and numerical stability for DDEs as demonstrated [2,3,9,13], among which we refer the reader to the monograph by Bellen and Zennaro [2]. And there are many references studied the convergence of difference methods for partial functional differential equations as demonstrated [1,7,10]. However, up to now, there are few papers published with respect to the numerical stability
✩ This work is supported by the National Nature Science Foundation of China (no. 11101130), the Research Fund of the Heilongjiang Provincial Education Department for the academic backbone of the excellent young people (no. 1254G044), Science and Technology Innovation Team in Higher Education Institutions of Heilongjiang Province (no. 2014TD005) and the Heilongjiang University Science Funds for Distinguished Young Scholars (no. JCL201303). ∗ Tel.: +86 13936648074. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.amc.2015.04.104 0096-3003/© 2015 Elsevier Inc. All rights reserved.
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
161
analysis for partial functional differential equations. The paper by Tian [16] gives a complete analysis of the numerical stability of finite difference methods, which is easy to carry out, but has more limitations comparing with Galerkin methods; the paper by Liang, Shi and Lv [12] presents a complete analysis of the numerical stability of Galerkin methods for a partial differential equation with piecewise constant argument, which is a partial differential equation with a special delay term, and Liang, Liu and Lv [11] also investigate the numerical stability of θ -methods for this equation. The method of lines and its stability for another type of linear parabolic differential equations with a delay term has been discussed by Zubik-Kowal [19,20]. The delay-dependent stability for a two-dimensional parabolic functional-differential equation with both fixed and distributed delays has been analyzed by Huang and Vandewalle [9]. The Galerkin method for elliptic, parabolic, and hyperbolic problems has been studied in [4,14,17], etc. In this paper, the optimal convergence order and the asymptotic stability of the Galerkin method applied to the process of semidiscretization with respect to the spatial variable x to (1.1)are investigated. Then it is discrete with respect to the time variable t by using forward Euler, backward Euler and Crank–Nicolson schemes, and the optimal convergence orders and asymptotic stability of these three full discretization schemes are studied respectively. Some notations: Denote · and · r as the norms in L2 = L2 () and the Sobolev space Hr = Hr () = W2r (), respectively, so that for real-valued function v,
v = vL2 :=
|v(x)|2 dx
1/2
,
and for a positive integer r ,
⎞1/2
di v(x) := ⎝ 2 ⎠ . dxi ⎛
vr = vHr
i≤r
Let v(x), w(x) (x ࢠ ) be real-valued functions. Throughout this paper,
(v, w) :=
v(x)w(x)dx,
∇ v, ∇ w :=
∂ v(x) ∂ w(x) dx, ∂x ∂x
and C will denote a positive constant not necessarily the same at different occurrences, which may depends on a, b and t of (1.1), but is independent of h and t (the stepsize in t-direction). Assumptions: In this paper, assume that u(t) := u(·, t), ut (t) := ut (·, t), utt (t) := utt (·, t), uttt (t) := uttt (·, t) ∈ H01 (), and at t = nτ (n = −1, 0, 1, . . . ), the derivative of t is denoted as the left derivative. Definition 1.1. The zero solution of (1.1) is called asymptotically stable if the solution u(x, t) of (1.1) corresponding to any sufficiently differential function φ (x, t) with φ (0, t) = φ (π , t) = 0 satisfies
lim u(x, t) = 0,
t→∞
x ∈ .
Theorem 1.2. [16] The zero solution of (1.1) is delay-independently asymptotically stable for a > b. 2. Semidiscretization of Galerkin methods For the preliminary step of defining a spatially semidiscrete approximation solution to (1.1), we thus write the problem in weak form: multiply the equation by a smooth function v(x) which vanishes on ࢚ (or v ∈ H01 ()), then integrate over , and apply Green’s formula to the second and third terms for the first equation of (1.1), yields
(ut (·, t), v) + a(∇ u(·, t), ∇ v) + b(∇ u(·, t − τ ), ∇ v) = 0, ∀v ∈ H01 (), t > 0.
(2.1)
Let N > 0 be a given positive integer number, h := πN be a stepsize in x-direction and
Ih := {xn := nh, n = 0, 1, . . . , N} be a given mesh on [0, π ]. Set In ࣎ [xn − 1 , xn ] (1 ࣘ n ࣘ N). Define
Shr := {χ ∈ C () : χ |In ∈ P (r−1)(In ), 1 ≤ n ≤ N,
χ (0) = χ (π ) = 0},
where P(r − 1) denotes the space of all (real) polynomials of degree not exceeding r − 1. We may then pose the approximate problem to find uh (t) ࣎ uh ( ·, t), belonging to Shr for each t, such that
(uh,t (t), χ) + a(∇ uh (t), ∇χ )) + b(∇ uh (t − τ ), ∇χ )) = 0, ∀χ ∈ Shr , t > 0, uh (x, t) = φh (x, t), −τ ≤ t ≤ 0,
(2.2)
where φ h ( ·, t) is some an approximation of φ ( ·, t) in Shr . N
In terms of the basis { j }1 h for Shr , our semidiscrete problem may be stated: find the coefficients α j (t) in uh (x, t) := such that Nh
j=1
αj (t)( j , k )) + a
Nh
j=1
αj (t)(∇ j , ∇ k )) + b
Nh
j=1
αj (t − τ )(∇ j , ∇ k )) = 0,
Nh j=1
αj (t) j (x)
162
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
k = 1, . . . , Nh , t > 0, and α j (t) = γ j (t) for j = 1, . . . , Nh and − τ ࣘ t ࣘ 0, with γ j (t) being the components of the given initial approximation φ h (x, t). In matrix notation this may be expressed as
Aα (t) + aBα(t) + bBα(t − τ ) = 0, for t > 0, α(t) = γ (t), for − τ ≤ t ≤ 0,
where A = (ajk ) is the mass matrix with elements ajk = ( j , k ), B = (bjk ) is the stiffness matrix with bjk = (࢟ j , ࢟ k ), α (t) and γ (t) are the vectors of α j and γ j respectively. The dimension of all these items equals Nh , the dimension of Shr . Since like the stiffness matrix B the mass matrix A is a Gram matrix, and thus in particular positive definite and invertible (see [15]), the above system of ordinary differential equations may be written as
α (t) + aA−1 Bα(t) + bA−1 Bα(t − τ ) = 0, for t > 0, α(t) = γ (t), for − τ ≤ t ≤ 0,
(2.3)
and hence obviously has a unique solution for positive t. 2.1. The convergence analysis Introduce the elliptic or Ritz projection Rh onto Shr as the orthogonal projection with respect to the inner product (࢟v, ࢟w), so that
(∇ Rh v, ∇χ) = (∇ v, ∇χ ), ∀χ ∈ Shr , for v ∈ H01 ().
(2.4)
Lemma 2.1. [15] Assume that for any v ∈ Hs () ∩ H01 (),
inf {v − χ + h∇(v − χ )} ≤ Chs vs ,
χ ∈Shr
for 1 ≤ s ≤ r
holds. Then, with Rh defined by (2.4) we have
Rh v − v + h∇(Rh v − v) ≤ Chs vs , for any v ∈ Hs () ∩ H01 (), 1 ≤ s ≤ r. Let u(t) ࣎ u( ·, t) and u : [0, +∞) → H01 (). Define Dh : H01 () → Shr by
a(∇ Dh u(t) − ∇ u(t), ∇χ ) + b(∇ Dh u(t − τ ) − ∇ u(t − τ ), ∇χ ) = 0, Dh u(t) = Rh u(t) = Rh φ(t),
for t > 0, χ ∈ Shr ,
for − τ ≤ t ≤ 0.
Theorem 2.2. Let u be the solution of (2.1). Assume that u(t) − Rh u(t)−l ࣘ that Dh u(t) ∈ Shr is the above projection of u(t), we have
(2.6) Chr + l u(t)r ,
− 1 ࣘ l ࣘ r − 2, −τ ࣘ t ࣘ 0 holds and
Dh u(t) − u(t)−l ≤ Chr+l , −1 ≤ l ≤ r − 2, t > 0. Proof. Denote
ρ u − Dh u = (u − Rh u) + (Rh u − Dh u) ξ + η, so for any χ ∈ Shr ,
a(∇η(t), ∇χ) = a(∇ Rh u(t) − ∇ Dh u(t), ∇χ ) = a(∇ u(t) − ∇ Dh u(t), ∇χ ) = b(∇ Dh u(t − τ ) − ∇ u(t − τ ), ∇χ ) = b(∇ Dh u(t − τ ) − ∇ Rh u(t − τ ) + ∇ Rh u(t − τ ) − ∇ u(t − τ ), ∇χ ) = −b(∇η(t − τ ) + ∇ξ (t − τ ), ∇χ ). Without loss of generality, let t ∈ ((k − 1)τ , kτ ], k ∈ N. Taking χ = η(t), then
η(t)1 ≤ C (η(t − τ )1 + ξ (t − τ )1 ) ≤ C (η(t − 2τ )1 + ξ (t − τ )1 + ξ (t − 2τ )1 ) ≤ · · · ≤ C (η(t − kτ )1 +
k
ξ (t − iτ )1 ),
i=1
by (2.6),
η(t − kτ ) = Rh u(t − kτ ) − Dh u(t − kτ ) = 0, therefore by using Lemma 2.1 with s = r,
η(t)1 ≤ C
k
i=1
ξ (t − iτ )1 = C
k
i=1
u(t − iτ ) − Rh u(t − iτ )1
(2.5)
(2.7)
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
≤ Chr−1
k
163
u(t − iτ )r .
i=1
By the definition of ρ and Lemma 2.1,
ρ(t)1 ≤ ξ (t)1 + η(t)1 ≤ Chr−1
k
u(t − iτ )r .
(2.8)
i=0
Thus the estimate in H1 () follows. To get other estimates, construct an auxiliary problem −w = ϕ in , with ϕ ∈ H01 (). Define wh = Rh w. By (2.5) and Lemma 2.1, there is
(ρ(t), ϕ) = (ρ(t), −w) = (∇ρ(t), ∇ w) = (∇ρ(t), ∇ w(t) − ∇ wh ) + (∇ρ(t), ∇ wh ) =
b a
b a
(∇ρ(t), ∇ w − ∇ wh ) − (∇ρ(t − τ ), ∇ wh − ∇ w) − (∇ρ(t − τ ), ∇ w)
≤ C ((ρ(t)1 + ρ(t − τ )1 )wh − w1 + ρ(t − τ )−l wl+2 ) ≤ C (hl+1 (ρ(t)1 + ρ(t − τ )1 ) + ρ(t − τ )−l )ϕl By using the definition of the norm ρ−l ,
ρ(t)−l ≤ C (hl+1 (ρ(t)1 + ρ(t − τ )1 ) + ρ(t − τ )−l ) ≤ C (hl+1 (ρ(t)1 + ρ(t − τ )1 + ρ(t − 2τ )1 ) + ρ(t − 2τ )−l ) ⎛
≤ · · · ≤ C ⎝h
l+1
k
⎞
ρ(t − iτ )1 + ρ(t − kτ )−l ⎠ ,
i=1
moreover, by the assumption for − τ ࣘ t − kτ ࣘ 0,
ρ(t − kτ )−l = u(t − kτ ) − Dh u(t − kτ )−l = u(t − kτ ) − Rh u(t − kτ )−l ≤ Chr+l u(t − kτ )r , which together with (2.8) complete the proof. The following theorem can be proved by a similar way to Theorem 2.2. Theorem 2.3. Let u be the solution of (2.1). Assume that ut (t) − Rh ut (t)i = ut (t) − (Rh u)t (t)i ࣘ Chr + i ut (t)r , i = 1, 2, −τ ࣘ t ࣘ 0 holds and ρ (t) ࣎ u(t) − Dh u(t). Then
ρt (t) + hρt (t)1 ≤ Chr , for t > 0.
(2.9)
Theorem 2.4. Let u and uh be the solutions of (2.1) and (2.2) respectively. Assume that u(t) − Rh u(t)i ࣘ Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Chr + i u(t)r ,
ut (t) −
uh (t) − u(t) ≤ Chr , for t > 0. Proof. Denote
uh (t) − u(t) = θ (t) + ρ(t), where θ (t) = uh (t) − Dh u(t),
ρ(t) = Dh u(t) − u(t).
(2.10)
The second term can be easily bounded by Theorem 2.2 and obvious there exists the following estimate:
ρ(t) ≤ Chr .
(2.11)
In order to estimate θ (t), note that by (2.5) and our definitions,
(θt (t), χ) + a(∇θ (t), ∇χ ) + b(∇θ (t − τ ), ∇χ ) = (uh,t , χ) − (Dh ut (t), χ ) + a(∇ uh (t), ∇χ ) − a(∇ Dh u(t), ∇χ ) + b(∇ uh (t − τ ), ∇χ ) − b(∇ Dh u(t − τ ), ∇χ ) = (uh,t , χ) − (Dh ut (t), χ ) + a(∇ uh (t), ∇χ ) − a(∇ u(t), ∇χ ) + b(∇ uh (t − τ ), ∇χ ) − b(∇ u(t − τ ), ∇χ ) = (ut (t) − Dh ut (t), χ ) = −(ρt (t), χ ), for ∀χ ∈ Shr . (2.12) Taking χ = θ (t),
1 d θ (t)2 + aθ (t)21 ≤ C (θ (t − τ )1 θ (t)1 + ρt (t)θ (t)). 2 dt So
a 1 d θ (t)2 + θ (t)21 ≤ C (θ (t − τ )21 + ρt (t)θ (t)). 2 dt 2
(2.13)
164
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
Without loss of generality, assume that t ∈ ((k − 1)τ , kτ ], k ∈ N. Integrating both sides of the above inequality from (k − 1)τ to t yields
θ (t)2 +
t
(k−1)τ
θ (s)21 ds ≤ C (θ ((k − 1)τ )2 +
Note that
t
(k−1)τ
θ (s)21 ds ≤ C (θ ((k − 1)τ )2 + = C (θ ((k − 1)τ )2 +
t
(k−1)τ t−τ
(k−2)τ
t
(k−1)τ
θ (s − τ )21 ds +
θ (s)21 ds +
t
ρt (s)θ (s)ds +
(k−1)τ
⎛ k−1
≤ C⎝ θ (iτ )2 + i=0
therefore
θ (t)
⎛ 2
≤ C⎝
k−1
θ (iτ ) +
i=0
So
⎛
θ (kτ ) ≤ C ⎝ 2
2
k−1
t−kτ
−τ
θ (iτ ) + 2
i=0
t−kτ −τ
t−τ
(k−2)τ
−τ
θ (s)
θ (s)
t
(k−3)τ
k−1
i=0
2 1 ds
+
k−1
2 1 ds
+
t−(k−i−1)τ iτ
(k−1)τ
ρt (s)θ (s)ds).
ρt (s)θ (s)ds)
ρt (s)θ (s)ds)
θ (s)21 ds
t−(k−i−1)τ
iτ
⎞
ρt (s)θ (s)ds⎠ , ⎞
ρt (s)θ (s)ds⎠ . ⎞
k−1 (i+1)τ
i=0
(k−1)τ
t
ρt (s)θ (s)ds) ≤ · · ·
θ (s)21 ds +
i=0
0
t
(k−1)τ t−2τ
≤ C (θ ((k − 1)τ )2 + θ ((k − 2)τ )2 + +
θ (s − τ )21 ds +
ρt (s)θ (s)ds⎠ ,
iτ
by the discrete Gronwall inequality (see [3,5]),
θ (kτ )2 ≤ C (θ (0)2 + therefore,
θ (t)2 ≤ C (θ (0)2 +
0 −τ
0
−τ
θ (s)21 ds +
θ (s)21 ds +
k−1 (i+1)τ
i=0
t 0
iτ
ρt (s)θ (s)ds),
ρt (s)2 ds +
0
t
θ (s)2 ds),
then θ (t) ࣘ Chr follows by Gronwall inequality (see [3]) and the assumption of the theorem, which together with (2.11) complete the proof. 2.2. The stability analysis Definition 2.5. If the solution uh (x, t) of (2.2) corresponding to any sufficiently differentiable function φ h (x, t) with φ h (0, t) = φ h (π , t) = 0 satisfies
lim uh (x, t) = 0,
t→∞
x ∈ [0, π ],
then the zero solution of (2.2) is called asymptotically stable. Theorem 2.6. The zero solution of (2.2) is asymptotically stable for a > b. Proof. Let α (t) = eλt C1 , where C1 is a constant vector. The characteristic equation of (2.3) is:
|λI + aA−1 B + bA−1 Be−λτ | = 0,
(2.14)
which is equivalent to
λ = −(a + be−λτ )λA−1 B , where λA−1 B denotes the corresponding eigenvalue of A−1 B. Moreover, because A and B are positive definite, so λA−1 B > 0. In fact, since B is a positive definite matrix, so there exists an invertible matrix C, such that B = CT C, and
CA−1 BC −1 = CA−1 C T CC −1 = CA−1 C T , which is positive definite and similar to A−1 B.
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
165
Let λ = λ1 + iλ2 , then
λ1 = −(a + be−λ1 τ cos λ2 τ )λA−1 B , λ2 = be−λ1 τ sin λ2 τ λA−1 B .
If λ1 ࣙ 0, then a + be−λ1 τ cos λ2 τ ≥ a − b and λ1 = −(a + be−λ1 τ cos λ2 τ )λA−1 B ≤ −(a − b)λA−1 B < 0, which is a contradiction. So if a > b, then λ1 < 0, i.e., the zero solution of (2.2) is asymptotically stable. 3. Fully discretization for Galerkin methods In this section, we shall focus our attention on numerical stability of backward Euler, Crank–Nicolson and forward Euler– Galerkin methods applied to (2.2). 3.1. The backward Euler–Galerkin method Let t = τ /m be a given stepsize with integer m ࣙ 1, the gridpoints tn be defined by tn = nt (n = 0, 1, 2, . . .), and Un be the approximation in Shr of u(t) at t = tn = nt, then the backward Euler–Galerkin method to (2.2) leads to a numerical process of the following type
U n − U n−1 ,χ t
+ a(∇ U n , ∇χ ) + b(∇ U n−m , ∇χ ) = 0,
for n > 0,
χ ∈ Shr ,
(3.1)
where Un ( · ) ࣎ φ h ( ·, tn ) for − m ࣘ n ࣘ 0. Nh Let U n (x) := i (x)αin . Then (3.1) can be reduced to i=1
Nh Nh Nh
1
1 n αi ( i , j ) − αin−1 ( i , j ) = −a αin (∇ i , ∇ j ) − b αin−m (∇ i , ∇ j ). t t Nh
i=1
i=1
i=1
(3.2)
i=1
With the notation as in the semidiscrete situation, this may be written as
(A + atB)α n = Aα n−1 − btBα n−m , for n > 0, αn = γ n, for − m ≤ n ≤ 0,
with γ n being the corresponding vector of components of the given initial approximation φ h (tn ), and α n := (α1 , . . . , αNh )T , and hence obviously has a unique solution for positive n. 3.1.1. The convergence analysis We shall prove the following error estimate: Theorem 3.1. Let u and Un be the solutions of (2.1) and (3.1), respectively. Assume that u(t) − Rh u(t)i ࣘ Chr + i u(t)r , ut (t) − Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Un − u(tn ) ≤ C (hr + t), for n = 1, 2, . . . . Proof. Denote
U n − u(tn ) = (U n − Dh u(tn )) + (Dh u(tn ) − u(tn )) = θ n + ρ n , and ρ n = ρ (tn ) is bounded as claimed in (2.11). This time, a calculation corresponding to (2.12) yields
θ n − θ n−1 t
,χ
+ a(∇θ n , ∇χ ) + b(∇θ n−m , ∇χ ) = −(wn , χ ),
∀χ ∈ Shr ,
where
Dh u(tn ) − Dh u(tn−1 ) − ut (tn ) = Dh − I ∂¯ u(tn ) + (∂¯ u(tn ) − ut (tn )) t u(tn ) − u(tn−1 ) = : wn1 + wn2 , ∂¯ u(tn ) := . t
wn =
Choosing χ = θ n , yields
θ n − θ n−1 t
, θ n + aθ n 21 + b(∇θ n−m , ∇θ n ) = −(wn , θ n ).
By using Schwarz inequality,
θ n − θ n−1 t
, θ n + θ n 21 ≤ C θ n−m 21 + wn θ n .
(3.3)
166
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
So
θ n 2 + tθ n 21 ≤ C (θ n−1 2 + tθ n−m 21 + (t)2 wn 2 ). Without loss of generality, assume that n ∈ ((k − 1)m, km], k ∈ N. Then
tθ n 21 ≤ C (θ n−1 2 + tθ n−m 21 + (t)2 wn 2 ) ≤ C (θ n−1 2 + θ n−m−1 2 + tθ n−2m 21 + (t)2 wn 2 + (t)2 wn−m 2 ) ⎞ ⎛ k−1 k−1
≤ ··· ≤ C⎝ θ n−im−1 2 + tθ n−km 21 + (t)2 wn−im 2 ⎠ . i=0
So
⎛
θ n 2 ≤ C ⎝
k−1
i=0
θ n−im−1 2 + tθ n−km 21 + (t)2
i=0
k−1
⎞
wn−im 2 ⎠ ,
i=0
therefore by using the discrete Gronwall inequality (see [3,5]),
⎛
θ ≤ C ⎝θ + tθ n 2
0 2
+ (t)
n−km 2 1
2
k−1
w
⎞
n−im 2 ⎠
.
(3.4)
i=1
Write
wn1 = (Dh − I)∂¯ u(tn ) = t−1
tn tn−1
d/dt (Dh − I)u(t)dt,
so
(t)2
k−1 k−1
n−im 2 = w1 i=0
i=0
≤
tn−im tn−im−1
k−1
tn−im tn−im−1
i=0
d/dt (Dh − I)u(t)dt Chr ut (t)r dt
2
2
2r
≤ Ch .
(3.5)
Further
twi2 = u(ti ) − u(ti−1 ) − tut (ti ) = −
ti ti−1
(t − ti−1 )utt (t)dt,
so that
(t)2
k−1
wn−im 2 = 2
i=0
k−1
i=0
≤ (t)2
tn−im tn−im−1
2
(t − tn−im−1 )utt (t)dt
k−1 tn−im
i=0
tn−im−1
2 utt (t)dt
≤ C (t) , 2
(3.6)
which together with (3.4), (3.5) and the assumption of the theorem complete the proof. 3.1.2. The stability analysis Definition 3.2. If the solution Un of (3.1) corresponding to any sufficiently differentiable function φ h (x, t) with φ h (0, t) = φ h (π , t) = 0 satisfies
lim U n = 0,
n→∞
x ∈ [0, 1],
then the zero solution of (3.1) is called asymptotically stable. Let K ࣎ [xi , xi + 1 ] be an element of the finite element, and Kˆ := [−1, 1] be the reference element in ξ -plane. Then
K
i j dx =
h 2
Kˆ
ˆ i ˆ j dξ ,
K
∇ i ∇ j dx =
2 h
Kˆ
ˆ i∇ ˆ j dξ . ∇
By (3.2),
αn =
h 2
−1 −1 2bt h h 2at 2at Bˆα n−m Aˆα n−1 − Aˆ + Bˆ Aˆ + Bˆ 2 h h 2 h
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
4at −1 −1 n−1 4bt 4at −1 −1 −1 n−m ˆ Bˆ ˆ Bˆ I+ α − , A A Aˆ Bˆα I + h2 h2 h2
=
where Aˆ and Bˆ is the corresponding mass and stiffness matrices. Let α n = λn C1 . The characteristic equation of (3.7) is:
λm − 1 +
4at λ ˆ−1 ˆ h2 A B
−1
λm−1 +
4bt h2
1+
4at λ ˆ−1 ˆ h2 A B
−1
167
(3.7)
λAˆ−1 Bˆ = 0,
(3.8)
ˆ Therefore the backward Euler scheme is delay-independently where λAˆ−1 Bˆ denotes the corresponding eigenvalue of Aˆ−1 B. asymptotically stable if and only if (3.8) is Schur polynomial for all m ࣙ 1(that is, the modulus of all zeros of the characteristic polynomial is less than 1). In order to prove that the characteristic polynomial is a Schur polynomial, the following lemma is needed. Lemma 3.3. [13] Let γ m (z) = α (z)zm − β (z) be a polynomial, where α (z) and β (z) are polynomial of constant order. Then the polynomial γ m (z) is a Schur polynomial for any m ࣙ 1 if and only if the following conditions hold (a) α (z) = 0⇒|z| < 1, (b) |β(z)| ≤ |α(z)|, ∀z ∈ C, |z| = 1 and (c) γm (z) = 0, ∀z ∈ C, |z| = 1, ∀m ≥ 1. Theorem 3.4. Suppose that 0 ࣘ b < a. Then the zero solution of the backward Euler scheme is delay-independently asymptotically stable. Proof. Denote
4at λ ˆ−1 ˆ h2 A B
α(λ) = λ − 1 + and
β(λ) = −
4bt h2
−1 4at λ λAˆ−1 Bˆ . 1+ ˆ−1 ˆ h2 A B
(a) If α (λ) = 0, then |λ| = (1 + (b) For ∀λ ∈ C, |λ| = 1,
|α(λ)|
−1
4at λAˆ−1 Bˆ )−1 h2
< 1.
−1 4at λAˆ−1 Bˆ 4at h2 = λ − 1 + λ ˆ−1 ˆ = λ − 1 + h2 A B 1 + 4at λAˆ−1 Bˆ h2 ≥
λ
1
4at Aˆ−1 Bˆ h2 + 4at Aˆ−1 Bˆ h2
λ
>
λ
1
4bt Aˆ−1 Bˆ h2 + 4at Aˆ−1 Bˆ h2
λ
= |β(λ)|.
(c) By (b), it is trivial. The proof is completed. 3.2. The Crank––Nicolson–Galerkin method The Crank––Nicolson–Galerkin method to (2.2) corresponding to (3.1) for the backward Euler methods leads to a numerical process of the following type
U n − U n−1 ,χ t
+a
∇ Un + ∇ Un−1 2
, ∇χ
+b
∇ Un−m + ∇ Un−m−1 2
, ∇χ
= 0,
(3.9)
for n > 0, where Un ( · ) = φ h ( ·, tn ) for − m ࣘ n ࣘ 0. Nh Let U n (x) := i (x)αin . Then (3.9) can be reduced to i=1
Nh 1
1 n αi ( i , j ) − αin−1 ( i , j ) t t Nh
i=1
i=1
Nh a n b
=− αi + αin−1 (∇ i , ∇ j ) − αin−m + αin−m−1 (∇ i , ∇ j ) for n > 0. 2 2 Nh
i=1
i=1
With the notation as in the semidiscrete situation, this may be written as
A + 2a tB α n = A − 2a tB α n−1 − t 2b Bα n−m − t 2b Bα n−m−1 , n > 0,
αn = γ n,
and hence obviously has a unique solution for positive n.
(3.10)
168
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
3.2.1. The convergence analysis We shall prove the following error estimate: Theorem 3.5. Let u and Un be the solutions of (2.1) and (3.9), respectively. Assume that u(t) − Rh u(t)i ࣘ Chr + i u(t)r , ut (t) − Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Un − u(tn ) ≤ C (hr + (t)2 ), for n = 1, 2, . . . . Proof. Denote
U n − u(tn ) = (U n − Dh u(tn )) + (Dh u(tn ) − u(tn )) = θ n + ρ n , and ρ n = ρ (tn ) is bounded as claimed in (2.11). This time, a calculation corresponding to (2.12) yields
θ n − θ n−1 t
,χ
+a
∇θ n + ∇θ n−1 2
, ∇χ
+b
∇θ n−m + ∇θ n−m−1 2
, ∇χ
= −(wn , χ), ∀χ ∈ Shr ,
(3.11)
where
Dh u(tn ) − Dh u(tn−1 ) ut (tn ) + ut (tn−1 ) − t 2 u ( t t n ) + ut (tn−1 ) = (Dh − I)∂¯ u(tn ) + ∂¯ u(tn ) − =: wn1 + wn2 . 2
wn =
Choosing χ =
θ n + θ n−1
, yields
2
n n−m θ + θ n−1 2 + ∇θ n−m−1 + b ∇θ , + a t 2 2 2 1 θ n + θ n−1 = − wn , . 2
θ n − θ n−1 θ n + θ n−1
,
By using Schwarz inequality,
θ n − θ n−1 θ n + θ n−1 t
So
,
2
∇θ n + ∇θ n−1
2
n n n−1 n−m + θ n−m−1 2 θ + θ n−1 2 ≤ C θ + w n θ + θ . + 2 2 2 1 1
n n−m 2 2 θ + θ n−1 2 θ + θ n−m−1 n−1 2 n 2 ≤C θ + t + ( t ) w . 2 2
θ n 2 + t
1
1
Without loss of generality, assume that n ∈ ((k − 1)m, km], k ∈ N. Then
2 n n−m θ + θ n−1 2 θ + θ n−m−1 n−1 2 2 n 2 t + t ≤ C θ + (t) w 2 2 1 1 n−2m 2 θ + θ n−2m−1 + (t)2 (wn 2 + wn−m 2 ) ≤ C θ n−1 2 + θ n−m−1 2 + t 2 1 ⎛ ⎞ k−1 k−1 θ n−km + θ n−km−1 2
≤ ··· ≤ C⎝ θ n−im−1 2 + t wn−im 2 ⎠ . + (t)2 2 i=0
Therefore
⎛
θ n 2 ≤ C ⎝
k−1
i=0
1
i=0
⎞ k−1 θ n−km + θ n−km−1 2
θ n−im−1 2 + t wn−im 2 ⎠ . + (t)2 2 i=0
1
By the assumption of the theorem and by using the discrete Gronwall inequality (see [3,5]),
⎛
θ ≤ C ⎝θ n 2
0 2
⎞ k−1 θ n−km + θ n−km−1 2
2 n−im 2 ⎠ + t w . + (t) 2 1
Write
wn1 = (Dh − I)∂¯ u(tn ) = t−1
tn tn−1
(Dh − I)ut (t)dt,
i=0
(3.12)
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
169
so
(t)2
k−1 k−1
n−im 2
w1 ≤ i=1
i=1
tn−im tn−im−1
Chr ut (t)r dt
2 ≤ Ch2r .
(3.13)
Further
ti ut (ti ) + ut (ti−1 ) 2 ( t ) − u ( t ) − t ( t ) u ≤ C uttt (t) dt, twi2 = i−1 i 2 ti−1
so that
(t)2
k−1 k−1
n−im 2 w2 ≤ C (t)4 i=1
tn−im
tn−im−1
i=1
uttt (s) dt
2
≤ C (t)4 ,
(3.14)
which together with (3.12), (3.13) and the assumption of the theorem complete the proof. 3.2.2. The stability analysis By (3.10),
−1 −1 h bt h h at at at Bˆ(α n−m + α n−m−1 ) Aˆ + Bˆ Aˆ − Bˆ α n−1 − Aˆ + Bˆ 2 h 2 h h 2 h 2at −1 −1 2bt 2at −1 2at −1 −1 −1 ˆ Bˆ ˆ Bˆ α n−1 − ˆ Bˆ = I+ A A A Aˆ Bˆ(α n−m + α n−m−1 ). I − I + h2 h2 h2 h2
αn =
(3.15)
Let α n = λn C1 . The characteristic equation of (3.15) is:
λm −
1− 1+
2at h2 2at h2
2bt λAˆ−1 Bˆ m−1 λAˆ−1 Bˆ h2 λ − (λ + 1) = 0. 2at λAˆ−1 Bˆ 1 + h2 λAˆ−1 Bˆ
(3.16)
Theorem 3.6. Suppose that 0 ࣘ b < a. Then the zero solution of the Crank–Nicolson scheme is delay-independently asymptotically stable. Proof. Denote α(λ) = λ −
1− 2at λAˆ−1 Bˆ 2 h 1+ 2at λAˆ−1 Bˆ h2
(a) If α (λ) = 0, then |λ| =
and β(λ) =
1− 2at λAˆ−1 Bˆ 2 h
1+ 2at λAˆ−1 Bˆ 2
2bt λAˆ−1 Bˆ h2 1+ 2at λAˆ−1 Bˆ 2 h
λ+1 .
< 1.
h
(b) For ∀λ ∈ C, |λ| = 1, denote λ = cos θ + isin θ , we have
2i sin θ λ − 1 cos θ − 1 + i sin θ = = . λ + 1 cos θ + 1 + i sin θ 2 + 2 cos θ So we have
2at λAˆ−1 Bˆ (λ − 1 ) −1 B 2 λA ˆ ˆ h = +
2at 2at λ+1 1 + h2 λAˆ−1 Bˆ λ + 1 1 + 2 λAˆ−1 Bˆ h 2at 2bt β(λ) λAˆ−1 Bˆ λAˆ−1 Bˆ h2 h2 . > = ≥ 2at 2at λ + 1 1 + 2 λ ˆ−1 ˆ 1 + 2 λ ˆ−1 ˆ
α(λ) λ − λ + 1 =
1− 2at 2 λAˆ−1 Bˆ h 1+ 2at h2
h
A
B
h
A
B
(c) By (b), it is trivial. The proof is completed. 3.3. The forward Euler–Galerkin method The forward Euler–Galerkin method to (2.2) corresponding to (3.1) for the backward Euler method leads to a numerical process of the following type
U n − U n−1 ,χ t
+ a(∇ U n−1 ,
∇χ ) + b(∇ Un−m−1 , ∇χ ) = 0, for n > 0.
where Un ( · ) = vh ( ·, tn ) for − m ࣘ n ࣘ 0.
(3.17)
170
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π . Fig. 1. The semidiscretization Galerkin method with h = 10
π . Fig. 2. The semidiscretization Galerkin method with h = 20
Let U n (x) :=
Nh i=1
i (x)αin . Then (3.1) can be reduced to
Nh Nh Nh Nh
1
1
αin ( i , j ) − αin−1 ( i , j ) = −a αin−1 (∇ i , ∇ j ) − b αin−m−1 (∇ i , ∇ j ). t t i=1
i=1
i=1
(3.18)
i=1
With the notation as in the semidiscrete situation, this may be written as
n Aα = (A − atB)α n−1 − btBα n−m−1 , for n > 0, αn = γ n, for − m ≤ n ≤ 0, and hence obviously has a unique solution for positive n.
3.3.1. The convergence analysis By using the same method of Section 3.1.1 as the proof of the convergence for the backward Euler–Galerkin method, the following theorem can be easily obtained. Theorem 3.7. Let u and Un be the solutions of (2.1) and (3.1). Assume that u(t) − Rh u(t)i ࣘ Chr + i u(t)r , ut (t) − Rh ut (t)i ࣘ Chr + i ut (t)r , − τ ࣘ t ࣘ 0 and φ h (t) − φ (t)i ࣘ Chr , i = 0, 1, then
Un − u(tn ) ≤ C (hr + t), for n = 1, 2, · · · .
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π . Fig. 3. The semidiscretization Galerkin method with h = 40
π , t = τ . Fig. 4. The backward Euler–Galerkin method with h = 10 100
π , t = τ . Fig. 5. The backward Euler–Galerkin method with h = 20 400
171
172
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = Fig. 6. The backward Euler–Galerkin method with h = 40
τ
1600
.
π , t = τ . Fig. 7. The Crank–Nicolson–Galerkin method with h = 10 100
3.3.2. The stability analysis By (3.18),
2 −1 h 4bt −1 n−m 2at Aˆ Bˆα Aˆ Aˆ − Bˆ α n−1 − h 2 h h2 4at −1 4bt −1 n−m ˆ Bˆ α n−1 − = I− . A Aˆ Bˆα h2 h2
αn =
(3.19)
Let α n = λn C1 . The characteristic equation of (3.7) is:
λm − 1 −
4at λ ˆ−1 ˆ h2 A B
λm−1 +
4bt λ ˆ−1 ˆ = 0. h2 A B
(3.20)
Theorem 3.8. Suppose that 0 ࣘ b < a. Then the zero solution of the backward Euler scheme is delay-independently asymptotically stable if and only if t ≤ 4aλ1 . h2 Aˆ−1 Bˆ
Proof. Denote α(λ) = λ − (1 −
4at λAˆ−1 Bˆ ) and h2
(a) If α (λ) = 0, then |λ| = 1 −
4at λAˆ−1 Bˆ h2
β(λ) = − 4bt λAˆ−1 Bˆ . h2
< 1 if and only if
t h2
1 . 2aλAˆ−1 Bˆ 1 if t ≤ , 4aλAˆ−1 Bˆ h2
<
(b) For ∀λ ∈ C, |λ| = 1, we assume that λ ࣎ cos θ + isin θ , then
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
173
π , t = τ . Fig. 8. The Crank–Nicolson–Galerkin method with h = 20 200
π , t = τ . Fig. 9. The Crank–Nicolson–Galerkin method with h = 40 400
Table 1 The 2-norm of the errors for the semidiscrezation Galerkin method at t = 1.
N
m
u(tn ) − Un (h = π /N)
u(tn ) − Un (h = π /N) u(tn ) − Un (h = π /2N)
10 20 40
100 400 1600
4.4903e-06 2.8574e-07 1.7939e-08
– 1.5715e+01 1.5928e+01
4at
2
|α(λ)|2 = λ − 1 + 2 λAˆ−1 Bˆ h
2 4at 4at = 2(1 − cos θ ) 1 − λ λ + ˆ−1 ˆ ˆ−1 ˆ h2 A B h2 A B 2 2 4at 4bt ≥ λ ˆ−1 ˆ > λ ˆ−1 ˆ = |β(λ)|2 ; h2 A B h2 A B
Otherwise, if
t h2
>
1
4aλAˆ−1 Bˆ
, then
|α(λ)|2
<
( 4t λ )2 a2 . h2 Aˆ−1 Bˆ
Choose
b 4t λ = |β(λ)|. h2 Aˆ−1 Bˆ Therefore,
t h2
≤
1 4aλAˆ−1 Bˆ
is a sufficient and necessary condition.
b2
=
|α(λ)|2 a2 + 4t ( 2 λ ˆ−1 ˆ )2 h
A
2
B
, so
|α(λ)|
4t λ −1 h2 Aˆ Bˆ
< b < a and |α(λ)| <
174
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = τ . Fig. 10. The forward Euler–Galerkin method with h = 10 100
π , t = τ . Fig. 11. The forward Euler–Galerkin method with h = 20 400
Table 2 The 2-norm of the errors for the backward Euler–Galerkin method at t = 1. N
m
u(tn ) − Un
u(tn ) − Un (h = π /N, t = τ /m) u(tn ) − Un (h = π /2N, t = τ /4m)
10 20 40
100 400 1600
4.1615e–06 2.5857e–07 1.6136e–08
– 1.6094e+01 1.6024e+01
Table 3 The 2-norm of the errors for the Crank–Nicolson–Galerkin method at t = 1. N
m
u(tn ) − Un (h = π /N, t = τ /m)
u(tn ) − Un (h = π /N, t = τ /m) u(tn ) − Un (h = π /2N, t = τ /2m)
10 20 40
100 200 400
4.5255e–06 2.8796e–07 1.8078e–08
– 1.5716e+01 1.5929e+01
(c) By (b), it is trivial. The proof is completed. Remark 3.9. In Theorem 3.8, λAˆ−1 Bˆ is an implicit expression, which will change as r (the degree of the polynomial) and Nh change. We list the following three cases for r = 2: • For Nh = 10, max λAˆ−1 Bˆ ≈ 2.7900e + 00;
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = Fig. 12. The forward Euler–Galerkin method with h = 40
τ
1600
175
.
π . Fig. 13. The semidiscrezation Galerkin method with h = 10
Table 4 The 2-norm of the errors for the forward Euler–Galerkin method at t = 1. N
m
u(tn ) − Un (h = π /N, t = τ /m)
u(tn ) − Un (h = π /N, t = τ /m) u(tn ) − Un (h = π /2N, t = τ /4m)
10 20 40
100 400 1600
3.2208e–05 2.0297e–06 1.2712e–07
– 1.5868e+01 1.5967e+01
• For Nh = 20, max λAˆ−1 Bˆ ≈ 2.9453e + 00; • For Nh = 40, max λAˆ−1 Bˆ ≈ 2.9862e + 00. 4. Numerical experiments In this section, the same parameters as [16] are used: φ (x, t) = sin x, τ = 1, a = 1.5, b = 1. We take r = 2. In Figs. 1–12, the errors of the semidiscretization, backward Euler, Crank–Nicolson and forward Euler Galerkin methods with the stepsize π , h = π and h = π for the semidiscretization Galerkin method, h = π , t = τ ; h = π , t = τ ; h = π , t = τ for h = 10 20 40 10 100 20 400 40 1600 π , t = τ ; h = π , t = τ ; h = π , t = τ for the Crank–– the backward and forward Euler–Galerkin methods and h = 10 100 20 200 40 400 Nicolson–Galerkin method are drawn. From these figures, it can be observed that the convergence orders obey the theoretical analysis. In Tables 1–4, the norms of the errors for the semidiscretization Galerkin and the three full discrete Galerkin methods at t = 1, u(1)−u (1)(h=π /N) )−Un (h=π /N,t=τ /m) and the ratios of u(t)−u h(t)(h=π /2N) for the semidiscretization Galerkin methods, uu(t(tn)−U n (h=π /2N,t=τ /4m) for the backward and h
n
176
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = τ . Fig. 14. The backward Euler–Galerkin method with h = 10 40
π , t = τ . Fig. 15. The Crank–Nicolson–Galerkin method with h = 10 40
)−Un (h=π /N,t=τ /m) forward Euler–Galerkin methods and uu(t(tn)−U n (h=π /2N,t=τ /2m) for the Crank––Nicolson–Galerkin method are given, respectively. n It can be seen that the convergence orders obey our theoretical analysis. π , t = τ , τ respectively for the forward Euler–Galerkin method; the In Figs. 13– 16, the same stepsize as [16] is taken: h = 10 40 50 π for the semidiscrezation Galerkin methods and h = π , t = τ for backward Euler and the Crank––Nicolson–Galerkin h = 10 10 40 methods. From Fig. 13, it can be seen that the semidiscrezation Galerkin method is asymptotically stable. From Figs. 14 and 15, it can 1 are asymptotically stable. be observed that the backward Euler and the Crank–Nicolson Galerkin methods with stepsize t = 40 1 1 From Figs. 16 and 17, it can be seen that for the two stepsizes t = 40 and t = 50 , the forward Euler methods are unstable in [16], the finite difference method with the first stepsize is unstable, but which with the second stepsize is asymptotically stable). 1 for the forward Euler–Galerkin method. It can be observed that the numerical solution with this In Fig. 18, we take t = 500 stepsize is asymptotically stable. Therefore, the mesh ratio for the forward Euler–Galerkin method is stricter than which for the forward Euler finite difference method. Conclusion. In practice, the stability analysis of a numerical method is very important. However, up to now, there are few papers published with respect to the numerical stability analysis for partial functional differential equations, no other than the numerical stability analysis of Galerkin methods for partial functional differential equations. In this paper, a complete analysis of the asymptotic stability of the semidiscrete (Galerkin methods) scheme and the full (forward, backward Euler and Crank– Nicolson methods) schemes for linear parabolic differential equations with delay are given: the semidiscrete, backward Euler and Crank–Nicolson full discrete schemes can unconditionally preserve the delay-independent asymptotic stability, but an additional restriction on time and spatial stepsize of the forward Euler full discrete scheme is needed to preserve the delay-independent asymptotic stability. These results for stability are consistent with the corresponding results of the finite difference methods. Moreover, numerical experiment has shown that in general, the addition restriction for the forward Euler–Galerkin full discrete
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
π , t = τ . Fig. 16. The forward Euler–Galerkin method with h = 10 40
π , t = τ . Fig. 17. The forward Euler–Galerkin method with h = 10 50
π , t = τ . Fig. 18. The forward Euler–Galerkin method with h = 10 500
177
178
H. Liang / Applied Mathematics and Computation 264 (2015) 160–178
scheme is stricter than which for the forward Euler finite difference method. In addition, it is proved that the convergence order of the semidiscrete scheme and the full discrete schemes are consistent with those of the original parabolic equation. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]
A. Ashyralyev, D. Agirseven, On convergence of difference schemes for delay parabolic equations, Comput. Math. Appl. 66 (2013) 1232–1244. A. Bellen, M. Zennaro, Numerical Methods for Delay Differential Equations, Clarendon Press, Oxford, 2003. H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, Cambridge University Press, Cambridge, 2004. H. Brunner, K. Mustapha, H. Mustapha, D. Schoetzau, An HP-version discontinuous Galerkin method for integro-differential equations of parabolic type, SIAM J. Numer. Anal. 49 (2011) 1369–1396. C.M. Chen, T. Shih, Finite Element Methods for Integrodifferential Equations, World Scientific, Singapore, 1998. D. Green, H.W. Stech, Diffusion and hereditary effects in a class of population models, in: S.N. Busenberg, K.L. Cooke (Eds.), Differential Equations and Applications in Ecology, Epidemics, and Population Problems, Academic Press, New York, 1981. W. Gu, P. Wang, A Crank–Nicolson difference scheme for solving a type of variable coefficient delay partial differential equations, J. Appl. Math. 2014 (2014) 1–6. P.J.V.d. Houwen, B.P. Sommeijer, C.T.H. Baker, On the stability of predictor-corrector methods for parabolic equations with delay, IMA J. Numer. Anal. 6 (1986) 1–23. C.M. Huang, S. Vandewalle, An analysis of delay-dependent stability for ordinary and partial differential equations with fixed and distributed delays, SIAM J. Sci.Comput. 25 (2004) 1608–1632. Y.F. Jin, J.X. Jiang, C.M. Hou, D.H. Guan, New difference scheme for general delay parabolic equations, J. Inform. Comput. Sci. 18 (2012) 5579–5586. H. Liang, M.Z. Liu, W.J. Lv, Stability of θ -schemes in the numerical solution of a partial differential equation with piecewise continuous arguments, Appl. Math. Lett. 23 (2010) 198–206. H. Liang, D.Y. Shi, W.J. Lv, Convergence and asymptotic stability of Galerkin methods for a partial differential equation with piecewise constant argument, Appl. Math. Comput. 217 (2010) 854–860. M.Z. Liu, M.N. Spijker, The stability of the θ -methods in the numerical solution of delay differential equations, IMA J. Numer. Anal. 10 (1990) 31–48. D.Y. Shi, W. Gong, J.C. Ren, A new stable second order nonconforming mixed finite element scheme for the stationary Stokes and Navier–Stokes equations, Math. Comput. Model. 53 (2011) 1956–1969. V. Thomée, Galerkin Finite Element Methods for Parabolic Problem, Springer, New York, 1997. H.J. Tian, Asymptotic stability of numerical methods for linear delay parabolic differential equations, Comput. Math.Appl. 56 (2008) 1758–1765. J.P. Wang, R. Zhang, Maximum principles for p1-conforming finite element approximations of quasi-linear second order elliptic equations, SIAM J. Numer. Anal. 50 (2012) 626–642. J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer-Verlag, New York, 1996. B. Zubik-Kowal, The method of lines for parabolic differential-functional equations, IMA J. Numer. Anal. 17 (1997) 103–123. B. Zubik-Kowal, Stability in the numerical solution of linear parabolic equations with delay term, BIT 41 (2001) 191–206.