STABILITY OF IMPLICIT–EXPLICIT BACKWARD DIFFERENCE FORMULAS FOR NONLINEAR PARABOLIC EQUATIONS GEORGIOS AKRIVIS Abstract. We analyze stability properties of BDF methods of order up to 5 for linear parabolic equations as well as of implicit–explicit BDF methods for nonlinear parabolic equations by energy techniques; time dependent norms play also a key role in the analysis.
1. Introduction Let T > 0, u0 ∈ H, and consider two abstract initial value problems, one for a linear parabolic equation, { ′ u (t) + A(t)u(t) = 0, 0 < t < T, (1.1) u(0) = u0 , and one for a possibly nonlinear parabolic equation, { ′ u (t) + A(t)u(t) = B(t, u(t)), 0 < t < T, (1.2) u(0) = u0 , in a usual triple of separable complex Hilbert spaces V ⊂ H = H ′ ⊂ V ′ , with V densely and continuously embedded in H. Here A(t) : V → V ′ are linear operators, while the operators B(t, ·) : V → V ′ may be nonlinear. We denote by (·, ·) both the inner product in H and the antiduality pairing between V ′ and V, and by | · | and ∥ · ∥ the norms in H and V, respectively. The space V ′ may be considered the completion of H with respect to the dual norm ∥ · ∥⋆ , ∀v ∈ V ′
∥v∥⋆ :=
|(v, w)| = sup |(v, w)|. w∈V w∈V \{0} ∥w∥
sup
∥w∥=1
Date: January 29, 2015. 2010 Mathematics Subject Classification. Primary 65M12, 65M60; Secondary 65L06. Key words and phrases. Parabolic equations, BDF methods, stability, energy technique, time dependent norms. Work partially supported by GSRT-ESET “Excellence” grant 1456. 1
2
GEORGIOS AKRIVIS
For q = 1, . . . , 6, consider the implicit q−step BDF method (α, β) and the explicit q−step method (α, γ) described by the polynomials α, β and γ, q q ∑ ∑ 1 q−j j ζ (ζ − 1) = αi ζ i , β(ζ) = ζ q , α(ζ) = j i=0 j=1 (1.3) q−1 ∑ q q q γi ζ i . γ(ζ) = β(ζ) − βq (ζ − 1) = ζ − (ζ − 1) = i=0
The BDF methods are A−stable for q = 1 and q = 2, i.e., A(ϑq )−stable with ϑ1 = ϑ2 = 90◦ , and A(ϑq )−stable for q = 3, . . . , 6 with ϑ3 = 86.03◦ , ϑ4 = 73.35◦ , ϑ5 = 51.84◦ and ϑ6 = 17.84◦ ; see [10, Section V.2]. Their order is q. For a given α, the scheme (α, γ) is the unique explicit q−step scheme of order q; the order of all other explicit q−step schemes (α, γ˜ ) is at most q − 1. Let N ∈ N, N ≥ q, and consider a uniform partition tn := nk, n = 0, . . . , N, of the interval [0, T ], with time step k := T /N. Assuming we are given starting approximations U 0 , . . . , U q−1 ∈ V, we discretize (1.1) in time by the q−step BDF method, i.e., we define approximations U m ∈ V to the nodal values um := u(tm ) of the exact solution as follows q ∑ (1.4) αi U n+i + kA(tn+q )U n+q = 0, i=0
n = 0, . . . , N − q. With the same notation, we discretize (1.2) in time by the implicit– explicit q−step (α, β, γ)−scheme,
(1.5)
q ∑
αi U n+i + kA(tn+q )U n+q = k
i=0
q−1 ∑
γi B(tn+i , U n+i ),
i=0
n = 0, . . . , N − q. The scheme (1.5) is referred to as the implicit–explicit q−step BDF method. The unknown U n+q appears only on the left-hand side of (1.5); therefore, to advance in time, we only need to solve one linear equation, which reduces to a linear system if we discretize also in space, at each time level. The stability results for the schemes (1.4) and (1.5), respectively, combined with the consistency of the methods for the underlying equations, lead to optimal order a priori error estimates for the initial value problems (1.1) for the (inhomogeneous) linear equation and for (1.2), respectively.
1.1. Abstract setting. Natural conditions for the parabolicity of the abstract equation in (1.1) are coercivity and boundedness of the operators A(t) : V → V ′ , i.e., (1.6)
Re(A(t)v, v) ≥ κ(t)∥v∥2
∀v ∈ V,
and (1.7)
∥A(t)v∥⋆ ≤ ν(t)∥v∥
∀v ∈ V,
respectively, with two smooth positive functions κ, ν : [0, T ] → R.
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
3
In the stability analysis of the implicit–explicit scheme (1.5) we assume, in addition, that B(t, ·) satisfies the following local Lipschitz condition in a ball Bu(t) := {v ∈ V : ∥v−u(t)∥ ≤ 1}, centered at the value u(t) of the solution u at time t, and, for simplicity, defined here in terms of the norm of V, ˜ ∥B(t, v) − B(t, v˜)∥⋆ ≤ λ(t)∥v − v˜∥ + µ ˜|v − v˜| ∀v, v˜ ∈ Bu(t) , (1.8) ˜ : [0, T ] → R and an arbitrary for all t ∈ [0, T ], with a smooth nonnegative function λ constant µ˜. Using (1.6) and (1.7), existence and uniqueness of the approximations U q , . . . , U N can be easily established by the Lax–Milgram lemma.
1.2. An auxiliary result by Nevanlinna & Odeh. Based on Dahlquist’s G−stability theory, Nevanlinna and Odeh [15] proved the following result for BDF methods of order up to five; this result allows us to establish stability by the energy method. Lemma 1.1. ([15]) Let α ∈ Pq , q ≤ 5, be the generating polynomial of the q−step BDF method; see (1.3). Let (·, ·) be an inner product with associated norm | · |. Then, there exist 0 ≤ ηq < 1, a positive definite symmetric matrix G = (gij ) ∈ Rq,q and reals δ0 , . . . , δq such that for v 0 , . . . , v q in the inner product space, q q q q ∑ 2 (∑ ) ∑ ∑ Re αi v i , v q − ηq v q−1 = gij (v i , v j ) − gij (v i−1 , v j−1 ) + δi v i . i=0
i,j=1
i,j=1
i=0
The smallest possible values of ηq are η1 = η2 = 0, η3 = 0.0836, η4 = 0.2878, η5 = 0.8160.
□
1.3. Stability results. We establish stability of the BDF scheme (1.4) and local stability of the implicit–explicit BDF scheme (1.5), for q = 1, . . . , 5, under the sufficient stability conditions (1.9)
∀t ∈ [0, T ]
ηq ν(t) < κ(t)
and (1.10)
∀t ∈ [0, T ]
˜ < κ(t), ηq ν(t) + (2q − 1)(1 + ηq )λ(t)
respectively. Using time-dependent norms, we relax these stability conditions to (2.29) and (3.30), respectively. For early, influential work on implicit multistep methods for nonlinear stiff differential equations and for linear parabolic equations, respectively, we refer to [13] and [16]; in [16] a rather complete analysis of strongly A(ϑ)−stable multistep schemes is presented. Implicit–explicit multistep schemes, for linear parabolic equations, were introduced and analyzed in [8]; such methods for nonlinear equations are studied in [3] and [2]. The stability analysis in [3, 2] is based on spectral and Fourier techniques, and led to sharp stability conditions; the analysis was motivated by similar techniques used in [13] and [16].
4
GEORGIOS AKRIVIS
Our analysis here is based on energy techniques and is motivated by the auxiliary Lemma 1.1; this Lemma was recently used first in [14] for the analysis of implicit BDF methods for a class of linear parabolic equations and subsequently in [4] both for BDF methods and some computationally less expensive variants for quasi-linear parabolic equations. In contrast to [4], here we restrict our attention to the case of linear operators A(t), i.e., operators independent of the solution u, and establish stability of the methods under less stringent assumptions. The advantage of the energy approach is that it is elementary and the proofs are quite short; the drawback is that it does not always lead to sharp results and does not apply to all strongly A(ϑ)−stable multistep methods; having said this, let us emphasize that this approach allows us to substantially improve some results of [2], as we will see. For the analysis of BDF methods in the case of linear parabolic equations with time independent, positive definite self-adjoint operator A we refer to [17, Chapter 10]. Energy techniques are employed to A−stable multistep schemes in [18] and to the three-step BDF method for the Navier–Stokes equations in [5]. Parabolic equations arise, for instance, in time dependent diffusion problems, such as the transient flow of heat according to Fourier’s law of heat conduction. There are implicit Runge–Kutta methods, and the closely related Galerkin time stepping schemes, that combine excellent stability properties for parabolic equations, such as A−stability and B−stability, with high order of accuracy; however, these methods are computationally time consuming per time step and, when applied to parabolic equations, suffer, in general, from the so-called order reduction phenomenon; cf. [17, Chapters 8 and 12] and the references therein. High order linear multistep methods, on the other hand, do not have so good stability properties as they can not be A−stable, for instance, according to the famous second Dahlquist barrier, and their stability properties are sensitive to variations of the time step; these methods do not suffer from order reduction when applied to parabolic equations, require at every time level only the solution of an equation (system) of the form of the implicit Euler method and are, thus, computationally efficient, provided, of course, that they are stable for the underlying equation. To make implicit schemes of either class implementable for nonlinear differential equations, we need to linearize at some stage of the discretization process. Implicit–explicit multistep methods are designed to discretize the linear part of the equations implicitly and the nonlinear part explicitly and can thus have good stability properties for some classes of equations; they are very efficient computationally, since they only require solving one linear equation (system) of the form of the implicit Euler method for the corresponding linear equation at every time level. Exponential integrators with Krylov subspace implementations, cf. [11, 12], as well as Chebyshev Runge–Kutta methods, cf. [1] and the references therein, are competing methods for parabolic equations. The relative merits and disadvantages of each of the methods depends on the particular equation and the available fast numerical linear algebra, such as multigrid preconditioners for linear systems.
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
5
The time stepping schemes (1.4) and (1.5) can not be implemented in this form, since they require solving a linear elliptic equation at every time level. If we want to really compute approximations, we need to combine the time stepping schemes with discretization in space, for instance, by the finite element method; then, the corresponding elliptic equations reduce to linear algebraic systems of equations and can be solved. Our approach extends easily to the fully discrete case; cf., e.g., [3]. For the mathematical theory of finite element methods we refer the reader to [7]; in particular, the conditioning of finite element equations is discussed in [7, §9.6–9.8]. Numerical methods for linear algebraic systems are analyzed in [9, 6]. The local Lipschitz condition (1.8) is usually satisfied in the applications in balls Bu(t) , centered at the value u(t) of the solution u at time t, defined in terms of L∞ −based Sobolev norms, often different for each argument, rather than in terms of L2 −based Sobolev norms. In such cases, our analysis does not directly apply if we only consider the discretization in time, since, in contrast to the present framework (see Theorem 3.2 in the sequel), it can not ensure that the approximations are sufficiently close to the exact solution in the required norm; it does, however, apply, usually under mild mesh-conditions, in the fully discrete case, i.e., if we combine the time stepping schemes with discretization in space; cf., e.g., [3]. An outline of the paper is as follows: Section 2 is of somewhat preparatory nature but the results may be of independent interest: we study stability properties of the implicit BDF methods for the linear equation (1.1). In Section 3 we present our main results: we establish stability properties of the implicit–explicit BDF schemes for the nonlinear equation (1.2) and derive optimal order error estimates. 2. Implicit BDF methods for linear equations This section is devoted to the analysis of stability properties of the implicit BDF methods (1.4) for the linear equation (1.1). It is well known that without additional assumptions on the functions κ and ν one can show stability only for A−stable methods. The three-, four- and five-step BDF methods do not have this property; the point then is which conditions on the functions κ and ν, or more precisely, on their ratio λ(t) := ν(t)/κ(t), t ∈ [0, T ], are suitable to allow us to establish stability of the schemes (1.4); the stability condition will depend on the specific method, i.e., on q ∈ {3, 4, 5}. We mention in passing that the ratio λ(t) does depend on the norm ∥ · ∥; if this norm is replaced by an equivalent norm, the ratio is multiplied by a constant factor; this factor is 1 only in the case the two norms are multiples of each other. This is the main reason why we will also rely on time dependent norms to get by by less stringent stability conditions. 2.1. Necessary stability condition. Here we will use the von Neumann criterion to show that the condition ν(t) 1 (2.1) max λ(t) = max ≤ 0≤t≤T 0≤t≤T κ(t) cos ϑ
6
GEORGIOS AKRIVIS
on the ratio λ(t) is necessary for the stability of an A(ϑ)−stable method, with ϑ < π/2, when applied to (1.1). We let φ ∈ (− π2 , π2 ) and consider the initial value problem for the parabolic equation (2.2)
˜ = − cos φAu ˜ − i sin φAu, ˜ ut = −Au = −eiφ Au
t ∈ (0, T ],
with A˜ : V → V ′ a positive definite self-adjoint bounded operator. The eigenvalues of eiφ A˜ are of the form reiφ , with a positive number r, i.e., they lie on the half-line in the complex plane starting at the origin and forming angle φ with the positive real half-axis. For |φ| > ϑ, this half-line is not contained in the stability sector Sϑ := {z ∈ C : z = reiφ , r ≥ 0, |φ| ≤ ϑ} of A(ϑ)−stable methods; thus, according to the von Neumann criterion, such a method can not be stable for this equation. Let us now ˜ v)1/2 . Then, the determine λ(t). The most suitable norm in V is ∥v∥ := |A˜1/2 v| = (Av, dual norm ∥ · ∥⋆ in V ′ is ∥v∥⋆ = |A˜−1/2 v| = (v, A˜−1 v)1/2 . Now, for v ∈ V, we have ˜ v) = cos φ(Av, ˜ v) + i sin φ(Av, ˜ v), (eiφ Av, ˜ v) = (cos φ)∥v∥2 ; we infer that κ(t) = cos φ. Furthermore, obviwhence Re(eiφ Av, ously, ˜ L(V,V ′ ) = |eiφ | ∥A∥ ˜ L(V,V ′ ) = ∥A∥ ˜ L(V,V ′ ) = 1, ∥eiφ A∥
whence ν(t) = 1. Therefore, λ(t) = 1/ cos φ. Since a necessary stability condition is |φ| ≤ ϑ, or, equivalently, cos φ ≥ cos ϑ, we infer that a necessary condition for the stability of A(ϑ)−stable methods for all equations satisfying the assumptions (1.6) and (1.7), expressed in terms of ϑ and the ratio λ(t) = ν(t)/κ(t), is (2.1). Let us close this subsection by mentioning that Savaré [16] imposed the condition that all complex numbers (A(t)v, v), for all v ∈ V and for all t ∈ [0, T ], are contained in a sector Sϑ˜, for some ϑ˜ < ϑ, and established stability of all strongly A(ϑ)−stable multistep schemes, avoiding any explicit condition on the ratio ν(t)/κ(t). 2.2. First sufficient stability condition. Here we will derive a sufficient stability condition, expressed in terms of the ratio λ(t) = ν(t)/κ(t), for the scheme (1.4). Motivated by the approach in [15], [14] and [4], we will use an energy method with suitable weight to establish stability of (1.4). Proposition 2.1 (Stability of the implicit BDF scheme (1.4)). Assume (1.6) and (1.7). Then, for q ∈ {1, . . . , 5}, under the stability condition (2.3)
κ(t) − ηq ν(t) ≥ ρ > 0,
the BDF method (1.4) is stable in the sense that, for k sufficiently small, ∑ 1 ∑ ℓ 2 cq |U | + ρk ∥U ∥ ≤ Cq |U j |2 + ηq ck∥U q−1 ∥2 , 2 ℓ=q j=0 n
(2.4)
q−1
n 2
for n = q, . . . , N, with cq , Cq positive constants depending only on q, and c a constant depending only on the maximum of ν.
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
7
Proof. We take in (1.4) the inner product with U n+q − ηq U n+q−1 , and then real parts to obtain q (∑ ) (2.5) Re αi U n+i , U n+q − ηq U n+q−1 + kIn+q = 0 i=0
with
( ) In+q := Re A(tn+q )U n+q , U n+q − ηq U n+q−1 .
(2.6)
The first term on the left-hand side of (2.5) can be taken care of exactly as in [15], [14] and [4]: With the notation U n := (U n−q+1 , . . . , U n )T and the norm |U n |G given by q ∑ n 2 gij (U n−q+i , U n−q+j ), |U |G = i,j=1
from Lemma 1.1 we have q (∑ ) Re αi U n+i , U n+q − ηq U n+q−1 ≥ |U n+q |2G − |U n+q−1 |2G . i=0
Thus, (2.5) yields (2.7)
|U n+q |2G − |U n+q−1 |2G + kIn+q ≤ 0.
It now remains to estimate In+q from below in a suitable way. First, we have ( ) ( ) In+q = Re A(tn+q )U n+q , U n+q − ηq Re A(tn+q )U n+q , U n+q−1 , and thus, in view of the coercivity condition (1.6), ( ) (2.8) In+q ≥ κ(tn+q )∥U n+q ∥2 − ηq Re A(tn+q )U n+q , U n+q−1 . To estimate the second term on the right-hand side of (2.8), we notice that ( ) Re A(tn+q )U n+q , U n+q−1 ≤ ∥A(tn+q )U n+q ∥⋆ ∥U n+q−1 ∥, whence, in view of the boundedness condition (1.7), ( ) Re A(tn+q )U n+q , U n+q−1 ≤ ν(tn+q )∥U n+q ∥ ∥U n+q−1 ∥, and hence
] ( ) ν(tn+q ) [ n+q 2 ∥U ∥ + ∥U n+q−1 ∥2 . Re A(tn+q )U n+q , U n+q−1 ≤ 2 In view of (2.9), estimate (2.8) leads to ] [ 1 1 n+q n+q (2.10) In+q ≥ κ(t ) − ηq ν(t ) ∥U n+q ∥2 − ηq ν(tn+q )∥U n+q−1 ∥2 . 2 2 From (2.7) and (2.10), we obtain
(2.9)
|U n+q |2G − |U n+q−1 |2G [ ] 1 1 + k κ(tn+q )− ηq ν(tn+q ) ∥U n+q ∥2 − k ηq ν(tn+q )∥U n+q−1 ∥2 ≤ 0; 2 2
8
GEORGIOS AKRIVIS
therefore, in view also of the stability condition (2.3), |U n+q |2G − |U n+q−1 |2G + ρk∥U n+q ∥2 (2.11) [ ] 1 + ηq kν(tn+q ) ∥U n+q ∥2 − ∥U n+q−1 ∥2 ≤ 0. 2 m+1 m e with L e the Lipschitz constant of ν, whence Now, |ν(t ) − ν(t )| ≤ Lk, e ν(tn+q ) ≤ ν(tn+q−1 ) + Lk;
(2.12) thus, estimate (2.11) yields
|U n+q |2G − |U n+q−1 |2G + ρk∥U n+q ∥2 [ ] 1 (2.13) 1 e 2 ∥U n+q−1 ∥2 . + ηq k ν(tn+q )∥U n+q ∥2 − ν(tn+q−1 )∥U n+q−1 ∥2 ≤ ηq Lk 2 2 Summing here from n = 0 to n = m − q, we obtain m ∑ 1 |U m |2G + ρk ∥U n ∥2 + ηq kν(tm )∥U m ∥2 ≤ |U q−1 |2G 2 n=q ∑ [ ] 1 e ∥U q−1 ∥2 + 1 ηq Lk e 2 + ηq k ν(tq−1 ) + Lk ∥U n ∥2 . 2 2 n=q m−1
e ≤ ρ, the last term on the right-hand side can For k sufficiently small such that ηq Lk be absorbed into the second term on the left-hand side, and we get m [ ] ρ ∑ n 2 1 m 2 |U |G + k ∥U ∥ ≤ |U q−1 |2G + ηq k ν(tq−1 ) + 1 ∥U q−1 ∥2 . 2 n=q 2
Using now the lower bound |U m |2G ≥ cq |U m |2 with cq the( smallest eigenvalue of ) the q−1 2 0 2 q−1 2 matrix G as well as the obvious estimate |U |G ≤ Cq |U | + · · · + |U | , we □ obtain the desired stability estimate (2.4). Notice that the sufficient stability condition (2.3) is void for q = 1, 2, and takes the form 1 (2.14) max λ(t) < 0≤t≤T ηq for q = 3, 4, 5. Remark 2.1 (Discrepancy between the sufficient and the necessary stability conditions). In the case of the q−step BDF methods, the values of the denominators η˜q := cos ϑq on the right-hand side of the necessary stability condition (2.1) are (2.15)
η˜3 = 0.0692,
η˜4 = 0.2865,
η˜5 = 0.6139,
η˜6 = 0.9524.
The ratios rq between the bounds of the sufficient and necessary stability conditions (2.14) and (2.1) are as follows: η˜3 0.0692 0.2865 0.6139 (2.16) r3 = = = 0.8277, r4 = = 0.9955, r5 = = 0.7523. η3 0.0836 0.2878 0.8160
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
9
In other words, the bounds of our sufficient stability conditions can be improved at most by 17.23%, 0.45% and 24.77%, for the three-, four- and five-step BDF methods, respectively. The result is particularly satisfactory for the four-step BDF scheme. □ Remark 2.2 (The case of a Gårding inequality). A straightforward modification of our stability proof yields stability under the sufficient stability condition (2.14) also in the case the coercivity condition (1.6) is replaced by a Gårding inequality, i.e., a positive multiple of |v|2 is subtracted from the right-hand side in (1.6), and/or a constant multiple of |v| is added to the right-hand side of (1.7). □ 2.3. Second sufficient stability condition: by means of time-dependent norms. In Remark 2.1, we discussed how much the constant ηq could be possibly decreased in the stability condition (2.14). Now, we focus on the function λ(t) in (2.14): since it depends on the choice of the norm of V, our effort here is to specify a suitable norm, such that λ(t) takes on smaller values. Notice that, for q = 3, 4, 5, the sufficient stability condition (2.14) may not be satisfied even in the case of positive definite selfadjoint operators A(t); the sufficient stability condition we will derive in this section is always satisfied for self-adjoint operators and is usually more favourable than (2.14) in the general case. Motivated by the approach in [14] and [4], where time-dependent norms were used in the case of self-adjoint operators, to get by by less stringent stability conditions we rely here on time dependent norms: We decompose the operators A(t) in their self-adjoint and anti-self-adjoint parts As (t) and Aa (t), respectively, ] ] 1[ 1[ As (t) := A(t) + A(t)⋆ , Aa (t) := A(t) − A(t)⋆ , 2 2 and introduce in V the time-dependent norm ∥ · ∥t by ∥v∥t := (As (t)v, v)1/2
∀v ∈ V.
We denote by ∥ · ∥⋆,t the corresponding dual norm on V ′ , ∀v ∈ V ′
∥v∥⋆,t :=
|(v, w)| = sup |(v, w)|. w∈V w∈V \{0} ∥w∥t
sup
∥w∥t =1
It follows easily from (1.6) and (1.7) that the norms ∥ · ∥t and ∥ · ∥ are equivalent, √ √ (2.17) κ(t) ∥v∥ ≤ ∥v∥t ≤ ν(t) ∥v∥ ∀v ∈ V. We denote by λa (t) : [0, T ] → [1, ∞) a smooth function such that (2.18)
∥A(t)v∥⋆,t ≤ λa (t)∥v∥t
∀v ∈ V.
An obvious consequence of (1.6) and (1.7) is that (2.18) is valid with λa (t) = λ(t) = ν(t)/κ(t). In general, however, (2.18) may be satisfied with λa (t) much smaller than λ(t); see Example 2.1 in the sequel. In the case of positive definite self-adjoint operators A(t), the estimate (2.18) holds as an equality with λa (t) = 1. The difference λa (t) − 1 may be viewed as a measure of the anti-self-adjoint part Aa (t) of A(t), or, in
10
GEORGIOS AKRIVIS
other words, as a measure of the deviation of A(t) from a positive definite self-adjoint operator. We will also use a mild Lipschitz condition on As (t), with respect to t, namely ( ) (2.19) ∥ As (t) − As (t˜) v∥⋆ ≤ L|t − t˜| ∥v∥ ∀t, t˜ ∈ [0, T ] ∀v ∈ V, with a Lipschitz constant L, and our goal is to prove stability under a sufficient stability condition expressed in terms of λa (t). Proposition 2.2 (Stability of the implicit BDF scheme (1.4)). Assume (1.6), (2.18) and (2.19). Then, for q ∈ {3, 4, 5}, under the stability condition (2.20)
1 − ηq λa (t) ≥ ρ > 0,
the BDF method (1.4) is stable in the sense that, for k sufficiently small, ∑ ∑ 1 cq |U | + ρκ⋆ k ∥U ℓ ∥2 ≤ Cq |U j |2 + ηq ck∥U q−1 ∥2 , 2 j=0 ℓ=q n
(2.21)
q−1
n 2
for n = q, . . . , N, with κ⋆ := min0≤t≤T κ(t), cq , Cq positive constants depending only on q, and c a constant depending only on the maximum of λa and ν. Proof. Our starting point is again (2.7). The analogue of (2.8) is in this case ( ) (2.22) In+q = ∥U n+q ∥2tn+q − ηq Re A(tn+q )U n+q , U n+q−1 . To estimate the second term on the right-hand side of (2.22), we notice that ( ) Re A(tn+q )U n+q , U n+q−1 ≤ ∥A(tn+q )U n+q ∥⋆,tn+q ∥U n+q−1 ∥tn+q ; therefore, in view of the boundedness condition (2.18), ( ) Re A(tn+q )U n+q , U n+q−1 ≤ λa (tn+q )∥U n+q ∥tn+q ∥U n+q−1 ∥tn+q , and hence
( ) λa (tn+q ) [ n+q 2 ] Re A(tn+q )U n+q , U n+q−1 ≤ ∥U ∥tn+q + ∥U n+q−1 ∥2tn+q . 2 In view of (2.23), relation (2.22) leads to [ ] 1 1 (2.24) In+q ≥ 1 − ηq λa (tn+q ) ∥U n+q ∥2tn+q − ηq λa (tn+q )∥U n+q−1 ∥2tn+q . 2 2 n+q−1 2 n+q−1 2 Next, we need to relate ∥U ∥tn+q back to ∥U ∥tn+q−1 . Since ∥v∥2t = ∥v∥2t˜ + ((As (t) − As (t˜))v, v), using the Lipschitz condition (2.19), we have
(2.23)
∥v∥2tn+q ≤ ∥v∥2tn+q−1 + Lk∥v∥2 ,
whence, in view of the first inequality in the equivalence of norms (2.17), ( ) L ∥v∥2tn+q ≤ 1 + k ∥v∥2tn+q−1 , n+q−1 κ(t ) and we obtain the desired estimate (2.25)
∥U n+q−1 ∥2tn+q ≤ (1 + ck)∥U n+q−1 ∥2tn+q−1
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
11
with c := L/ min0≤t≤T κ(t). Therefore, (2.24) yields [ ] 1 1 (2.26) In+q ≥ 1 − ηq λa (tn+q ) ∥U n+q ∥2tn+q − ηq λa (tn+q )(1 + ck)∥U n+q−1 ∥2tn+q−1 . 2 2 From (2.7) and (2.26), in view of the stability condition (2.20), we obtain |U n+q |2G − |U n+q−1 |2G + ρk∥U n+q ∥2tn+q (2.27) ] [ 1 + ηq kλa (tn+q ) ∥U n+q ∥2tn+q − (1 + ck)∥U n+q−1 ∥2tn+q−1 ≤ 0. 2 n+q Since λa (t ) ≤ λa (tn+q−1 )(1 + cˆk), with a suitable constant cˆ, the estimate (2.27) yields |U n+q |2G − |U n+q−1 |2G + ρk∥U n+q ∥2tn+q
] [ 1 n+q−1 n+q−1 2 n+q n+q 2 + ηq k λa (t )∥U ∥tn+q − λa (t )(1 + c˜k)∥U ∥tn+q−1 ≤ 0. 2 Summing here from n = 0 to n = m − q, we obtain
(2.28)
|U m |2G + ρk
m ∑
1 ∥U n ∥2tn + ηq kλa (tm )∥U m ∥2tm ≤ |U q−1 |2G 2 n=q
∑ 1 1 + ηq k(1 + c˜k)λa (tq−1 )∥U q−1 ∥2tq−1 + ηq c˜k 2 λa (tn )∥U n ∥2tn . 2 2 n=q m−1
Now, for k sufficiently small such that ηq c˜k max0≤t≤T λa (t) ≤ ρ, the last term on the right-hand side can be absorbed in the second term on the left-hand side; using then the equivalence of norms (2.17), we obtain m ∑ 1 1 e a (tq−1 )ν(tq−1 )∥U q−1 ∥2 , |U m |2G + ρκ⋆ k ∥U n ∥2 ≤ |U q−1 |2G + ηq k Cλ 2 2 n=q which is the desired stability estimate.
□
The sufficient stability condition (2.20) takes the form 1 (2.29) max λa (t) < 0≤t≤T ηq for q = 3, 4, 5; cf. (2.14). Remark 2.3 (Necessary stability conditions). It immediately follows from the example used in subsection 2.1 that the following two conditions, the first on λa (t) and the second on the ratio of the norms of Aa (t) and As (t), 1 , (2.30) max λa (t) ≤ 0≤t≤T cos ϑ and ∥Aa (t)∥L(V,V ′ ) (2.31) max ≤ tan ϑ, 0≤t≤T ∥As (t)∥L(V,V ′ )
12
GEORGIOS AKRIVIS
are necessary for the stability of an A(ϑ)−stable method, with ϑ < π/2, when applied to (1.1). □ Example 2.1 (A time-dependent second order elliptic operator). Here we demonstrate the advantages of the use of time-dependent norms with a simple example; we will see that the ratio λ(t) may take on much larger values than λa (t). We let ¯ × [0, T ] → R be Ω ⊂ Rd be a bounded domain with smooth boundary ∂Ω and a, b : Ω smooth functions, with a positive, and consider )a family of second order elliptic op( erators A(t), A(t)v := −∇ [a(x, t) + ib(x, t)]∇v , t ∈ [0, T ], subject to homogeneous Dirichlet boundary conditions. Obviously, A(t) is an operator from V := H01 (Ω) to V ′ = H −1 (Ω). A first choice of time-independent norm ∥ · ∥ in V is the standard H01 (Ω)−seminorm, ∥v∥ := |∇v|, with | · | the L2 (Ω)−norm. In this case the best possible choices of functions κ and ν satisfying (1.6) and (1.7) are κ(t) = min a(x, t), ¯ x∈Ω
ν(t) = max |a(x, t) + ib(x, t)|. ¯ x∈Ω
A second, more natural and, in general, better choice of a time-independent norm ∥ · ∥ in V is in this case (∫ )1/2 ∥v∥ := a(x, t⋆ )|∇v(x)|2 dx , Ω
for some fixed t ∈ [0, T ]. It is then easily seen that the best possible choices of functions κ and ν satisfying (1.6) and (1.7) are ⋆
κ(t) = min ¯ x∈Ω
a(x, t) , a(x, t⋆ )
ν(t) = max ¯ x∈Ω
|a(x, t) + ib(x, t)| . a(x, t⋆ )
Therefore, for instance, for a(x, t) = ef (t)g(x) and b = 0, the best possible choice of κ and ν leads to the ratio λ(t) = ν(t)/κ(t), [ ] |f (t) − f (t⋆ )| max g(x) − min g(x) ¯ ¯ x∈Ω x∈Ω λ(t) = e , which may take on large values. The use of time-dependent norms, on the other hand, leads to |a(x, t) + ib(x, t)| λa (t) = max ; ¯ a(x, t) x∈Ω notice that λa (t) is at most equal to and, in general, much smaller than λ(t) = ν(t)/κ(t). In particular, in the case of self-adjoint A(t), i.e., b = 0, if we use the time-dependent norm, then we have λa (t) = 1, t ∈ [0, T ]. □ 2.4. Comparison with results from [2]. Let us first recall a stability result from [2]: We decompose the operators A(t) in the form A(t) = A1 + A2 (t) with A1 a positive definite self-adjoint bounded operator from V to V ′ . Then, an A(ϑ)−stable multistep method is stable for (1.1), provided (2.32)
max ∥A2 (t)∥L(V,V ′ ) < sin ϑ ∥A1 ∥L(V,V ′ ) .
0≤t≤T
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
13
(Actually, this result was formulated in [2] in the case V is endowed with the norm ∥ · 1/2 ∥, ∥v∥ := |A1 v|; then, ∥A1 ∥L(V,V ′ ) = 1.) The ratio max0≤t≤T ∥A2 (t)∥L(V,V ′ ) /∥A1 ∥L(V,V ′ ) in (2.32) can be viewed as a measure of the deviation of the family of operators A(t), t ∈ [0, T ], from a time-independent, positive definite self-adjoint operator A1 . We shall compare the sufficient stability conditions (2.14) and (2.32) for the BDF methods in the case [ ] (2.33) A(t) = 1 + z(t) A1 = A1 + A2 (t) with A1 as above and z : [0, T ] → C a complex-valued function with real-part larger than −1, Re z(t) > −1, for all t ∈ [0, T ], such that the corresponding evolution equation be parabolic. We restrict our attention to the three-, four- and five-step BDF methods for two reasons: first, the sufficient stability condition (2.3) is void (and, in particular, optimal) for q = 1 and q = 2, and, second, for the one- and two-step BDF schemes stability results by the energy technique were established also in [2]. Before we proceed, let us note that the sufficient stability conditions (2.14) and (2.29) coincide for operators of the form (2.33). For the operators (2.33) condition (2.32) reads max |z(t)| < sin ϑ,
0≤t≤T
which in the case of the q−step BDF methods means (2.34)
max |z(t)| < λq ,
0≤t≤T
with
λq := sin ϑq .
Thus z(t) belongs to the open disc of radius sin ϑq in the complex plane centered at the origin. In other words, 1 + z(t) belongs to the open disc of radius sin ϑq centered at 1; see Figure 2.1, middle. Since 1 + z(t) must belong to the sector Sϑq , if we want the method to be stable, the radius sin ϑq on the right-hand side of (2.34) is optimal, it can not be increased. The values of the constants λq are (2.35)
λ3 = 0.9976, λ4 = 0.9581, λ5 = 0.7863, λ6 = 0.3063.
Let us now see what the new sufficient stability condition (2.14) means for the operators (2.33). For convenience, we use the norm ∥ · ∥, ∥v∥ := |A1/2 1 v|, v ∈ V. For v ∈ V, we have [ ] [ ] Re(A(t)v, v) = Re 1 + z(t) (A1 v, v) = Re 1 + z(t) ∥v∥2 and
] −1/2 [ 1/2 ∥A(t)v∥⋆ = |A1 1 + z(t) A1 v| = |1 + z(t)| |A1 v| = |1 + z(t)| ∥v∥, [ ] whence κ(t) = Re 1 + z(t) and ν(t) = |1 + z(t)|. Therefore, condition (2.14) reads
(2.36)
max
0≤t≤T
|1 + z(t)| 1 [ ]< . ηq Re 1 + z(t)
This means that 1 + z(t) belongs to the interior of a sector Sϑˆq , with ϑˆq < ϑq such that cos ϑˆq = ηq ; see Figure 2.1, right. Notice that the disc of radius sin ϑq centered at 1,
14
GEORGIOS AKRIVIS
see (2.34), is not entirely contained in the sector Sϑˆq , since ϑˆq < ϑq , i.e., the stability condition (2.36) is not always more favourable than (2.34). Notice also that, in contrast to (2.32), the stability condition (2.14) leads again to (2.36), if the operators A(t) in (2.33) are multiplied by a positive function σ(t). y C
y C
Sϑ
ϑ ϑ
x
y
r = sin ϑq r
Sϑ q
1
C
x
Sϑˆq
x
Figure 2.1. The stability sector Sϑ of A(ϑ)−stable methods in the complex plane, left; illustration in light blue of the stability conditions (2.34), middle, and (2.36), right, for the values of 1 + z(t).
3. Implicit–explicit BDF methods for nonlinear equations In this section we present the main results of the paper; we prove stability of the implicit–explicit BDF methods (1.5), of order up to five, for the nonlinear equation (1.2) and establish optimal order a priori error estimates. 3.1. First sufficient stability condition. Since the differential equation is in general nonlinear in this case, besides the approximations U n ∈ Bu(tn ) satisfying (1.5), we consider implicit–explicit BDF approximations V n ∈ Bu(tn ) such that (3.1)
q ∑
αi V n+i + kA(tn+q )V n+q = k
i=0
q−1 ∑
γi B(tn+i , V n+i ),
i=0
n = 0, . . . , N − q.
Theorem 3.1 (Stability of the implicit–explicit BDF scheme (1.5)). Assume (1.6), (1.7) and (1.8). Then, for q ∈ {1, . . . , 5}, under the stability condition (3.2)
∀t ∈ [0, T ]
˜ ≥ ρ > 0, κ(t) − ηq ν(t) − (2q − 1)(1 + ηq )λ(t)
the implicit–explicit BDF method (1.5) is locally stable in the sense that, with ϑm := U m − V m , for k sufficiently small, ∑( ) 1 ∑ ℓ 2 cq |ϑ | + ρk ∥ϑ ∥ ≤ C |ϑj |2 + k∥ϑj ∥2 , 2 ℓ=q j=0 n
(3.3)
q−1
n 2
for n = q, . . . , N, with cq a positive constant depending only on q, and C a constant independent of ρ, k, n and the approximations.
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
15
Proof. Letting bm := B(tm , U m ) − B(tm , V m ) and subtracting (3.1) from (1.5), we obtain q ∑
(3.4)
αi ϑn+i + kA(tn+q )ϑn+q = k
i=0
q−1 ∑
γi bn+i ,
i=0
n = 0, . . . , N − q. As in Section 2, we take in (3.4) the inner product with ϑn+q − ηq ϑn+q−1 , and take real parts to obtain
(3.5)
Re
q (∑
) αi ϑn+i , ϑn+q − ηq ϑn+q−1 + kIn+q = kJn+q
i=0
with (3.6)
( ) In+q := Re A(tn+q )ϑn+q , ϑn+q − ηq ϑn+q−1
and (3.7)
Jn+q := Re
q−1 (∑
γi b
n+i
,ϑ
n+q
− ηq ϑ
n+q−1
) .
i=0
With the notation Θn := (ϑn−q+1 , . . . , ϑn )T and the norm |Θn |G given by |Θn |2G
=
q ∑
gij (ϑn−q+i , ϑn−q+j ),
i,j=1
in view of Lemma 1.1, relation (3.5) yields the estimate (3.8)
|Θn+q |2G − |Θn+q−1 |2G + kIn+q ≤ kJn+q .
Furthermore, In+q can be estimated from below exactly as in the case of the implicit BDF scheme, (3.9)
[ ] 1 1 In+q ≥ κ(tn+q ) − ηq ν(tn+q ) ∥ϑn+q ∥2 − ηq ν(tn+q )∥ϑn+q−1 ∥2 ; 2 2
see (2.10). Therefore, all that remains to be done, is to estimate Jn+q from above in a suitable way. For simplicity of presentation, we assume µ˜ = 0 in the following; the general case can be treated similarly via a straightforward use of the discrete Gronwall inequality at the end of the proof. First, we have Jn+q ≤
q−1 ∑ i=0
( ) |γi | ∥bn+i ∥⋆ ∥ϑn+q ∥ + ηq ∥ϑn+q−1 ∥ ,
16
GEORGIOS AKRIVIS
whence, in view of the local Lipschitz condition (1.8), Jn+q ≤
q−1 ∑
( ) ˜ n+i )∥ϑn+i ∥ ∥ϑn+q ∥ + ηq ∥ϑn+q−1 ∥ |γi |λ(t
i=0
( ) 1∑ ˜ n+i ) ∥ϑn+i ∥2 + ∥ϑn+q ∥2 + ηq ∥ϑn+i ∥2 + ηq ∥ϑn+q−1 ∥2 |γi |λ(t 2 i=0 q−1
≤
q−1 q−1 )( (∑ ∑ ) 1 1 n+i n+i 2 ˜ ˜ n+i ) ∥ϑn+q ∥2 + ηq ∥ϑn+q−1 ∥2 . = (1 + ηq ) |γi |λ(t )∥ϑ ∥ + |γi |λ(t 2 2 i=0 i=0
˜ n+i ) ≤ λ(t ˜ n+q−j )+ Lk, b i = 0, . . . , q −1, j = 0, 1, and |γ0 |+· · ·+|γq−1 | = Now, since λ(t |γ(−1)| = 2q − 1, we easily see that q−1 ∑
˜ n+i ) ≤ (2q − 1)λ(t ˜ n+q−j ) + Ck, b |γi |λ(t
j = 0, 1.
i=0
Therefore, the above estimate for Jn+q yields q−1 ∑ [ ] 1 ˜ n+i )∥ϑn+i ∥2 + Ck b ∥ϑn+q ∥2 + ηq ∥ϑn+q−1 ∥2 Jn+q ≤ (1 + ηq ) |γi |λ(t 2 i=0 (3.10) [ n+q ] 1 q ˜ ˜ n+q−1 )∥ϑn+q−1 ∥2 . + (2 − 1) λ(t )∥ϑn+q ∥2 + ηq λ(t 2 In view of the stability assumption (3.2), from (3.8), (3.9) and (3.10) we infer that ( ) 1 |Θn+q |2G − |Θn+q−1 |2G + ρk∥ϑn+q ∥2 + ηq kν(tn+q ) ∥ϑn+q ∥2 − ∥ϑn+q−1 ∥2 2 q−1 ∑ 1 q n+q n+q 2 ˜ ˜ n+i )∥ϑn+i ∥2 + (2 − 1)(1 + ηq )k λ(t )∥ϑ ∥ ≤ (1 + ηq )k |γi |λ(t 2 (3.11) i=0 [ n+q ] 1 q ˜ ˜ n+q−1 )∥ϑn+q−1 ∥2 + (2 − 1)k λ(t )∥ϑn+q ∥2 + ηq λ(t 2 [ ] b 2 ∥ϑn+q ∥2 + ηq ∥ϑn+q−1 ∥2 . + Ck
Estimating the coefficient ν(tn+q ) of ∥ϑn+q−1 ∥2 on the left-hand side of (3.11) as in (2.12), proceeding as in the proof of Proposition 2.1 and using the fact that |γ0 | + · · · + |γq−1 | = 2q − 1, we easily arrive at the desired stability estimate (3.3), provided k is sufficiently small. □ The sufficient stability condition (3.2) can also be written in the form (3.12)
∀t ∈ [0, T ]
˜ < κ(t) ηq ν(t) + (2q − 1)(1 + ηq )λ(t)
˜ vanishes. How much the coefficient ηq of ν(t) in and reduces to (2.14) in case λ (3.12) could be possibly decreased, can be seen from the discussion in Remark 2.1. ˜ Concerning the coefficient of λ(t), it follows from the analysis in [3], where the case of time-independent, positive definite self-adjoint operator A(t) is considered, that the
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
17
best one can hope for is to get rid of ηq , i.e., (2q − 1)(1 + ηq ) can be at most decreased to 2q − 1. A priori error estimates are usually established by combining stability and consistency of the numerical method. As is typical for multistep methods, it is very easy to prove consistency of the scheme (1.5). All (local) stability results we present in this paper can be used to derive optimal order error estimates. As an example, we will next derive error estimates using the local stability result of Theorem 3.1; see, also, e.g., [3, 2, 4]. 3.1.1. Consistency. The order of the q−step methods (α, β) and (α, γ) is q, i.e., q ∑
(3.13)
ℓ
i αi = ℓq
ℓ−1
=ℓ
q−1 ∑
i=0
iℓ−1 γi ,
ℓ = 0, 1, . . . , q.
i=0
The consistency error E n of the scheme (1.5) for the solution u of (1.2), i.e., the amount by which the exact solution misses satisfying (1.5), is given by (3.14)
n
kE =
q ∑
αi u
n+i
+ kA(t
n+q
n+q
)u
−k
i=0
q−1 ∑
γi B(tn+i , un+i ),
i=0
n = 0, . . . , N − q. Here, um := u(tm ) denote the nodal values of the exact solution u(t). Letting
(3.15)
E1n
:=
q ∑
αi u
n+i
′
−ku (t
n+q
),
E2n
:= kB(t
n+q
,u
i=0
n+q
)−k
q−1 ∑
γi B(tn+i , un+i ),
i=0
and using the differential equation in (1.2), we infer that (3.16)
kE n = E1n + E2n .
Now, by Taylor expanding about tn and using the order conditions of the implicit (α, β)−scheme, i.e., the first equality in (3.13), and the second equality in (3.13), respectively, we obtain ] [ q ∫ tn+q ∑ ∫ tn+i 1 (tn+i − s)q u(q+1) (s) ds − qk (tn+q − s)q−1 u(q+1) (s) ds , αi En = 1 q! i=0 tn tn [ ∫ n+q ] ∫ tn+i q t ∑ k n n+q e (q) (s) ds − e (q) (s) ds , − s)q−1 B γi (tn+i − s)q−1 B E2 = (q − 1)! tn (t tn i=0 e with B(t) := B(t, u(t)), t ∈ [0, T ]. Thus, under obvious regularity requirements, we obtain the desired optimal order consistency estimate
(3.17)
max ∥E n ∥⋆ ≤ Ck q .
0≤n≤N −q
18
GEORGIOS AKRIVIS
3.1.2. Error estimates. Combining local stability and consistency we derive optimal order error estimates: Theorem 3.2 (Error estimate). Assume that the stability condition (3.2) is satisfied, that the solution u of (1.2) is sufficiently smooth such that the consistency estimate (3.17) be valid and that we are given starting approximations U 0 , U 1 , . . . , U q−1 ∈ V to u0 , . . . , uq−1 such that ( ) (3.18) max |uj − U j | + k 1/2 ∥uj − U j ∥ ≤ Ck q . 0≤j≤q−1
Let U q , . . . , U N ∈ V be recursively defined by (1.5), and en := un − U n , n = 0, . . . , N. Then, there exists a constant C, independent of k and m, such that, for k sufficiently small, q−1 m−q m {∑ } ∑ ∑ ( j2 ) m 2 ℓ 2 j 2 (3.19) |e | + k |e | + k∥e ∥ + k ∥E ℓ ∥2⋆ , ∥e ∥ ≤ C j=0
ℓ=0
ℓ=0
m = q − 1, . . . , N, and
max |u(tn ) − U n | ≤ Ck q .
(3.20)
0≤n≤N
Proof. According to (3.17) and (3.18), there exists a constant C⋆ such that the righthand side of (3.19) can be estimated by C⋆2 k 2q , q−1 N −q } {∑ ∑ j 2 j 2 ∥E ℓ ∥2⋆ ≤ C⋆2 k 2q . (3.21) C (|e | + k∥e ∥ ) + k j=0
ℓ=0
Now, obviously, (3.20) is a consequence of (3.19) and (3.21). Thus, it remains to prove the stability estimate (3.19). Subtracting (1.5) from (3.14), we have q ∑
αi e
n+i
+ kA(t
n+q
n+q
)e
q−1 ∑ [ ] =k γi B(tn+i , un+i ) − B(tn+i , U n+i ) + kE n .
i=0
i=0
If we take here the inner product with en+q − ηq en+q−1 , proceed exactly as in the proof of Theorem 3.1, and assume for the time being that U j ∈ Bu(tj ) , j = 0, . . . , n + q − 1, we easily arrive at ∑ ∑( ) 1 ∑ ℓ 2 cq |en+q |2 + ρk ∥e ∥ ≤ C Re(E ℓ , eℓ+q − ηq eℓ+q−1 ); |ej |2 + k∥ej ∥2 + k 2 ℓ=q j=0 ℓ=0 n+q
q−1
n
cf. (3.3). Now, bounding Re(E ℓ , eℓ+q − ηq eℓ+q−1 ) ≤
ρ ρηq 1 + ηq ℓ 2 ∥eℓ+q ∥2 + ∥eℓ+q−1 ∥2 + 2 ∥E ∥⋆ 4(1 + ηq ) 4(1 + ηq ) ρ
and summing up, we obtain ∑( ) 1 ∑ ℓ 2 1 + ηq ∑ ℓ 2 | + ρk k ∥E ∥⋆ , ∥e ∥ ≤ C |ej |2 +k∥ej ∥2 +kηq c∥eq−1 ∥2 +2 4 ℓ=q ρ j=0 ℓ=0 n+q
cq |e
n+q 2
q−1
n
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
19
and infer that (3.19) holds true for m = n + q. Now, the estimate (3.19) is obviously valid for m = q − 1. Assume inductively that it holds for m = q − 1, . . . , n + q − 1, 0 ≤ n ≤ N − q. Then, according to (3.21) and the induction hypothesis, we have, for k small enough, max
(3.22)
0≤j≤n+q−1
∥ej ∥ ≤ C⋆ k q−1/2 ≤ 1,
and thus U j ∈ Bu(tj ) , j = 0, . . . , n + q − 1. Therefore, as we proved above, (3.19) holds indeed for m = n + q as well, and the proof is complete. □ 3.2. Second sufficient stability condition: by means of time-dependent norms. As in subsection 2.3, we use here the time-dependent norms ∥ · ∥t and ∥ · ∥⋆,t to study stability properties of the implicit–explicit BDF methods (1.5) of order up to 5. In analogy to (1.8), we assume here that the operators B satisfy the local Lipschitz condition (3.23)
˜ b (t)∥v − v˜∥t + µ ∥B(t, v) − B(t, v˜)∥⋆,t ≤ λ ˜b |v − v˜|
∀v, v˜ ∈ Bu(t) ,
˜ b : [0, T ] → R and an arbitrary for all t ∈ [0, T ], with a smooth nonnegative function λ ˜ constant µ˜b . It follows easily from √ (1.8) and (1.6) that (3.23) is valid with λb (t) = ˜ λ(t)/κ(t) and µ˜b = µ˜/ min0≤t≤T κ(t). In general, however, (3.23) may be satisfied ˜ ˜ with λb (t) much smaller than λ(t)/κ(t); see Example 3.1 in the sequel.
Theorem 3.3 (Stability of the implicit–explicit BDF scheme (1.5)). Assume (1.6), (2.18), (2.19) and (3.23). Then, for q ∈ {1, . . . , 5}, under the stability condition (3.24)
∀t ∈ [0, T ]
˜ b (t) ≥ ρ > 0, 1 − ηq λa (t) − (2q − 1)(1 + ηq )λ
the implicit–explicit BDF method (1.5) is locally stable in the sense that, with ϑm := U m − V m , for k sufficiently small, ∑( ∑ ) 1 cq |ϑ | + ρκ⋆ k ∥ϑℓ ∥2 ≤ C |ϑj |2 + k∥ϑj ∥2 , 2 j=0 ℓ=q n
(3.25)
q−1
n 2
for n = q, . . . , N, with κ⋆ := min0≤t≤T κ(t), cq a positive constant depending only on q, and C a constant independent of ρ, k, n and the approximations. Proof. Our starting point is the estimate (3.26)
|Θn+q |2G − |Θn+q−1 |2G + kIn+q ≤ kJn+q
with In+q and Jn+q given in (3.6) and (3.7), respectively; see (3.8). First, In+q can be estimated from below as in the case of the implicit BDF method for the linear equation in the form [ ] 1 1 (3.27) In+q ≥ 1 − ηq λa (tn+q ) ∥ϑn+q ∥2tn+q − ηq λa (tn+q )(1 + ck)∥ϑn+q−1 ∥2tn+q−1 ; 2 2
20
GEORGIOS AKRIVIS
cf. (2.26). Therefore, it only remains to estimate Jn+q from above in a suitable way. As in the proof of Theorem 3.1, for simplicity of presentation, we again assume µ˜b = 0. Now, we have Jn+q ≤
q−1 ∑
( ) |γi | ∥bn+i ∥⋆,tn+i ∥ϑn+q ∥tn+i + ηq ∥ϑn+q−1 ∥tn+i ,
i=0
whence, in view of the local Lipschitz condition (3.23), proceeding as in the derivation of the corresponding estimation in the proof of Theorem 3.1, we obtain ∑ 1 ˜ b (tn+i )∥ϑn+i ∥2n+i |γi |λ ≤ (1 + ηq ) t 2 i=0 q−1
Jn+q
+
q−1 ( ) 1∑ ˜ b (tn+i ) ∥ϑn+q ∥2n+i + ηq ∥ϑn+q−1 ∥2n+i . |γi |λ t t 2 i=0
˜ b (tn+i ) ≤ λ ˜ b (tn+q−j ) + Lk, b i = 0, . . . , q − 1, j = 0, 1, Therefore, using the estimates λ q the fact that |γ0 | + · · · + |γq−1 | = 2 − 1, as well as estimates analogous to (2.25), we end up with an estimate of Jn+q of the form ∑ 1 ˜ b (tn+i )∥ϑn+i ∥2n+i Jn+q ≤ (1 + ηq ) |γi |λ t 2 i=0 [ ] 1 ˜ b (tn+q )∥ϑn+q ∥2n+q + ηq λ ˜ b (tn+q−1 )∥ϑn+q−1 ∥2n+q−1 + (2q − 1) λ t t 2 [ n+q 2 ] b ∥ϑ ∥ n+q + ηq ∥ϑn+q−1 ∥2n+q−1 , + Ck t t q−1
(3.28)
b with an appropriate constant C. Combining (3.26), (3.27) and (3.28), and proceeding as in the proof of Theorem 3.1, we easily infer that ∑( ) 1 ∑ n 2 e cq |ϑ | + ρk ∥ϑ ∥tn ≤ C |ϑj |2 + k∥ϑj ∥2tj , 2 j=q j=0 n
(3.29)
q−1
n 2
for n = q, . . . , N. The desired stability estimate (3.25) follows then immediately from (3.29) in view of the equivalence of the norms ∥ · ∥t and ∥ · ∥; see (2.17). □ The stability result of Theorem 3.3, combined with the consistency result (3.17), along the lines of Theorem 3.2, yields easily optimal order a priori error estimates, under the stability condition (3.24). The sufficient stability condition (3.24) can also be written in the form (3.30)
∀t ∈ [0, T ]
˜ b (t) < 1. ηq λa (t) + (2q − 1)(1 + ηq )λ
Remark 3.1 (Sufficient stability condition in the case of self-adjoint A(t)). In the case of positive definite self-adjoint operators A(t), we have λa (t) = 1 and the sufficient
STABILITY PROPERTIES OF IMPLICIT–EXPLICIT BDF METHODS
21
stability condition (3.30) takes the form (3.31)
∀t ∈ [0, T ]
˜ b (t) < 1 − ηq =: r˜q . (2q − 1)λ 1 + ηq
The values of r˜q are (3.32)
r˜1 = r˜2 = 1, r˜3 = 0.8457, r˜4 = 0.5530, r˜5 = 0.1013.
It is known that even in the case of time-independent, positive definite self-adjoint operator A, the constant on the right-hand side of (3.31) can not be replaced by a constant larger than 1, if we want the method to be stable under our conditions; see [3]. □ Remark 3.2 (The implicit–explicit one- and two-step BDF methods). In the case of the implicit–explicit one- and two-step BDF methods, the sufficient stability conditions (3.12) and (3.30) reduce to (3.33)
∀t ∈ [0, T ]
˜ < κ(t) (2q − 1)λ(t)
∀t ∈ [0, T ]
˜ b (t) < 1, (2q − 1)λ
and (3.34)
respectively, since η1 = η2 = 0, with q = 1, 2. Even in the case of time-independent, positive definite self-adjoint operator A, the coefficient 2q − 1 on the left-hand sides of (3.33) and (3.34) can not be replaced by a smaller one, if we want the method to be stable under our conditions; see [3]. In other words, the sufficient stability conditions (3.12) and (3.30) are sharp for the implicit–explicit one- and two-step BDF methods. □ ˜ b (t) and λ(t)/κ(t) ˜ Example 3.1 (Comparison between λ ). Here we demonstrate the advantages of the use of time-dependent norms with a simple example; we will see that ˜ ˜ b (t). With the notation of Exthe ratio λ(t)/κ(t) may take on much larger values than λ ample 2.1, we(slightly modify the definition of the operator A(t) by letting b = 0, i.e., ) A(t)v := −∇ a(x, t)∇v , t ∈ [0, T ], assuming again homogeneous Dirichlet boundary conditions. We have chosen A(t) self-adjoint here, since any anti-self-adjoint part ˜ ˜ b . Furtheradded to A(t) is irrelevant for κ(t), λ(t), the time-dependent(norms and for)λ more, we let B(t, ·) : V → V ′ be given by B(t, v) := −∇ c(x, t)d(v)∇v , t ∈ [0, T ], with c : Ω¯ × [0, T ] → R a smooth function and d : R → R√a smooth, bounded function. If we choose the time-independent norm ∥v∥ := | a(·, t⋆ ) ∇v|, with a fixed t⋆ ∈ [0, T ], in the real space V as in Example 2.1, then the best possible choices of ˜ satisfying (1.6) and (1.8) are functions κ and λ κ(t) = min ¯ x∈Ω
a(x, t) , a(x, t⋆ )
˜ = max |c(x, t)| sup |d(y)|. λ(t) ¯ a(x, t⋆ ) y∈R x∈Ω
22
GEORGIOS AKRIVIS
˜ As demonstrated in Example 2.1, the ratio λ(t)/κ(t) may take on large values. The use of time-dependent norms, on the other hand, leads to ˜ b (t) = max |c(x, t)| sup |d(y)|, λ ¯ a(x, t) y∈R x∈Ω ˜ ˜ b (t) is much smaller than which is at most equal to the ratio λ(t)/κ(t); in general, λ this ratio. □
3.3. Comparison with results from [2]. Let us first recall a stability result from [2]: As in §2.4, we decompose the linear operators A(t) in the form A(t) = A1 + A2 (t) with A1 a bounded, positive definite self-adjoint linear operator from V to V ′ . For convenience, we use the norm ∥ · ∥, ∥v∥ := |A1/2 1 v|, v ∈ V. Then, the implicit– explicit multistep scheme (1.5) is locally stable for (1.2), with A(t) as above and B(t, ·) satisfying the local Lipschitz condition (1.8), provided (3.35)
1 ˜ < 1. max ∥A2 (t)∥L(V,V ′ ) + (2q − 1) max λ(t) 0≤t≤T sin ϑq 0≤t≤T
Also, if we restrict our attention to linear constraints in the maxima of ∥A2 (t)∥ and ˜ λ(t), as in (3.35), then this condition is sharp, in the sense that none of the coefficients 1/ sin ϑq and 2q − 1 can be replaced by smaller constants; see [2]. We refer to [2] also for a necessary (nonlinear) stability condition on these quantities and its discrepancy to (3.35). Now, as in §2.4, we shall compare the sufficient stability conditions (3.12) and (3.35) in the case [ ] (3.36) A(t) = 1 + z(t) A1 = A1 + A2 (t) with A1 as above and z : [0, T ] → C a complex-valued function with real-part larger ′ than −1, Re z(t) > [ −1, for]all t ∈ [0, T ]; cf. (2.33). In this case we have ∥A2 (t)∥L(V,V ) = |z(t)|, κ(t) = Re 1 + z(t) and ν(t) = |1 + z(t)| (see §2.4). Therefore, the sufficient stability conditions (3.35) and (3.12) take the form (3.37)
1 ˜