IMPLICIT-EXPLICIT MULTISTEP FINITE ELEMENT METHODS FOR ...
IMPLICIT-EXPLICIT MULTISTEP FINITE ELEMENT METHODS FOR NONLINEAR PARABOLIC PROBLEMS GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS Abstract. We approximate the solution of initial boundary value problems for nonlinear parabolic equations. In space we discretize by finite element methods. The discretization in time is based on linear multistep schemes. One part of the equation is discretized implicitly and the other explicitly. The resulting schemes are stable, consistent and very efficient, since their implementation requires at each time step the solution of a linear system with the same matrix for all time levels. We derive optimal order error estimates. The abstract results are applied to the Kuramoto-Sivashinsky and the Cahn-Hilliard equations in one dimension, as well as to a class of reaction diffusion equations in Rν , ν = 2, 3.
1. Introduction In this paper we construct and analyze implicit-explicit multistep schemes for the time discretization of a class of nonlinear parabolic problems of the form: Given T > 0 and u0 ∈ H, find u : [0, T ] → D(A) such that (1.1)
u′ (t) + Au(t) = B(u(t)),
0 < t < T,
u(0) = u0 ,
where A is a linear, selfadjoint, positive definite operator on a Hilbert space (H, (·, ·)) with domain D(A) dense in H, and B : D(A) → H is a (possibly nonlinear) differentiable operator. Considering a finite dimensional subspace Vh of V, V := D(A1/2 ), we are led to a semidiscrete problem approximating (1.1): We seek a function uh , uh (t) ∈ Vh , defined by (1.2)
u′h (t) + Ah uh (t) = Bh (uh (t)),
0 < t < T,
uh (0) = u0h ;
here u0h ∈ Vh is a given approximation to u0 , and Ah , Bh are appropriate operators on Vh with Ah linear, selfadjoint and positive definite. Following an idea of Crouzeix, [3], for the time discretization of parabolic equations with time dependent coefficients, we combine implicit and explicit multistep schemes to discretize (1.2) in time: Implicit schemes are used for discretizing the left-hand side of the o.d.e. in (1.2), and explicit schemes for the nonlinear right-hand side. Thus, letting 1991 Mathematics Subject Classification. Primary 65M60, 65M12; Secondary 65L06. The work of the first and third authors was supported in part by a research grant from the University of Crete. 1
2
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
k be a (constant) time step, tn = nk, n = 0, . . . , N, T = Nk, we define a sequence of approximations U n , U n ∈ Vh , to un = u(tn ), by (1.3)
q X
αi U
n+i
+k
q X
βi Ah U
n+i
=k
γi Bh (U n+i ).
i=0
i=0
i=0
q−1 X
A scheme of this form is characterized by three polynomials α, β and γ, α(ζ) =
q X i=0
i
αi ζ ,
β(ζ) =
q X
i
βi ζ ,
i=0
γ(ζ) =
q−1 X
γi ζ i .
i=0
We note that, when Bh vanishes, the scheme (1.3) reduces to the implicit linear multistep method (α, β) (or (ρ, σ) in the notation of Dahlquist) for the equation u′h (t)+Ah uh (t) = 0; similarly if Ah vanishes, the scheme (1.3) reduces to the linear multistep method (α, γ), which is explicit since γq = 0, for the equation u′h (t) = Bh (uh (t)). The scheme (1.3) is a combination of the methods (α, β) and (α, γ); it is linearly implicit and nonlinearly explicit. We shall refer to it as the (α, β, γ) scheme. We shall assume in the sequel that the method (α, β) is strongly A(0)−stable; this implies, in particular, that αq βq is positive, which in turn ensures invertibility of the operator αq I + kβq Ah . Thus, given U 0 , . . . , U q−1 in Vh , U q , . . . , U N are well defined by the (α, β, γ) scheme. We will assume in the sequel, without loss of generality, that both αq and βq are positive. For the analysis of the scheme (1.3) we will need some additional assumptions for the operators A and B as well as for the finite dimensional spaces Vh ; the operators Ah and Bh will then be appropriately defined. Let, thus, | · | denote the norm of H, and introduce in V the norm k · k by kvk = (A1/2 v, A1/2 v)1/2 . We identify H with its dual, and denote by V ′ the dual of V , again by (·, ·) the duality pairing on V ′ and V, and by k · k⋆ the dual norm on V ′ . We assume that B can be extended to an operator from V into V ′ , which is differentiable, and an estimate of the form (1.4)
|(B ′ (v)w, ω)| ≤ λkwkkωk + µ(v)|w| |ω|
∀v, w, ω ∈ V
holds with a constant λ < 1 and a functional µ(v) which is bounded for v bounded in V. Indeed, depending on the particular (α, β, γ) scheme, we shall need to assume that the constant λ is appropriately small. Further, the assumption that µ(v) is bounded for v bounded in V will suffice to derive our results under some mild meshconditions; these conditions can be weakened if µ(v) is bounded for v bounded in some appropriate weaker norms, and even avoided if µ(v) is bounded for v bounded in H. We will assume in the sequel that (1.1) possesses a solution which is sufficiently regular for our results to hold. Uniqueness of smooth solutions follows easily in view of (1.4). For the space discretization we use a family Vh , 0 < h < 1, of finite dimensional subspaces of V. In the sequel the following discrete operators will play an essential role: Define Po : V ′ → Vh , Rh : V → Vh , Ah : V → Vh , and the nonlinear operator
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
3
Bh : V → Vh by (Po v, χ) = (v, χ) ∀χ ∈ Vh ,
(ARh v, χ) = (Av, χ) ∀χ ∈ Vh ,
(Ah ϕ, χ) = (Aϕ, χ) ∀χ ∈ Vh ,
(Bh (ϕ), χ) = (B(ϕ), χ) ∀χ ∈ Vh . Then, obviously, Ah Rh = Po A and Bh = Po B . In the error analysis we shall use the approximation properties of the “elliptic projection” Rh . We assume that Rh v is an approximation to v of order r, provided that v is sufficiently regular, (1.5)
|v − Rh v| + hd/2 kv − Rh vk ≤ M(v)hr ,
where r and d are two integers, 2 ≤ d ≤ r, and M(v) is bounded if an appropriate norm of v is bounded. We further assume that (1.6)
kB(v) − B(Rh v)k⋆ ≤ M(v)hr .
We emphasize that the condition (1.6) serves consistency purposes only. It is needed to prove consistency of the (α, β, γ) scheme—and, in fact, already of the semidiscrete problem (1.2)—for Rh u of optimal order with respect to h. It is somehow restrictive though, since it essentially means that, if A is a differential operator of order d, then B may contain derivatives of up to order d/2 only. For some concrete differential equations, however, one can get by with a less stringent condition by taking into account in the definition of Rh the terms of B of order higher than d/2. The scheme (1.3) is very efficient since its implementation requires at every time step solving a linear system with the same matrix for all time levels. If the order of both the implicit and the explicit scheme is p, then under some mild meshconditions and for appropriately small λ, and appropriate starting values U 0 , . . . , U q−1 , we derive the optimal order error estimate max |u(tn ) − U n | ≤ c(k p + hr ).
0≤n≤N
An outline of the paper is as follows: Section 2 is devoted to the analysis of a simple one step semiexplicit scheme of first order accuracy; its purpose is to motivate the analysis, in section 3, of more general multistep schemes of higher accuracy. Finally, in the last section, we apply our abstract results to three examples, namely the KuramotoSivashinsky equation and the Cahn-Hilliard equation in one space dimension, and to a class of reaction diffusion equations in Rν , ν = 2, 3. 2. An implicit-explicit one step scheme As a motivation for the analysis of multistep schemes, we study in this section the simplest implicit-explicit scheme which is a combination of the backward Euler and the Euler method.
4
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
Let U 0 ∈ Vh be given. We define fully discrete approximations U n ∈ Vh to u(tn ), recursively by U n+1 − U n + Ah U n+1 = Bh (U n ), n = 0, . . . , N − 1. k Our main concern in this section is to analyze the approximation properties of the sequence {U n }. As an intermediate step, we shall show consistency of the scheme (2.1) for the elliptic projection of the solution u of (1.1). Let (2.1)
W (t) = Rh u(t),
W (t) ∈ Vh , 0 ≤ t ≤ T .
We note for later use that, in view of the definition of Rh , kW (t)k is obviously bounded by ku(t)k, and thus bounded uniformly in h and t, (2.2)
sup sup kW (t)k ≤ C . h
t
Consistency. The consistency error E n of the scheme (2.1) for W is given by (2.3)
kE n = (W n+1 − W n ) + kAh W n+1 − kBh (W n ),
n = 0, 1, . . . , N − 1.
Using the relation Ah Rh = Po A and the definition of Bh , we rewrite E n as kE n = Rh (un+1 − un ) + kPo Aun+1 − kPo B(Rh un ),
and thus, in view of (1.1), E n = E1n + E2n , where (2.4)
kE1n =(Rh − Po )(un+1 − un ) + Po (un+1 − un − ku′ (tn+1 ) )
and
+ kPo (B(un+1) − B(un ) ),
E2n = Po (B(un ) − B(Rh un ) ).
(2.5) Obviously
max |E1n | ≤ C(k + hr ).
(2.6)
0≤n≤N −1
Further, in view of (1.6), (2.7)
max kE2n k⋆ ≤ Chr .
0≤n≤N −1
Next, we show optimal order rate of convergence of our approximations to the (sufficiently regular) solution of (1.1), provided that the initial approximation U 0 ∈ Vh satisfies (2.8)
|u0 − U 0 | + hd/2 ku0 − U 0 k ≤ c(u)hr .
We first introduce some more notation: For v ∈ V and b > 0, let 1/2 |||v|||b := |v|2 + bkkvk2 , and set |||v||| := |||v|||1 . Further, let
(2.9)
m := sup{µ(v) : kvk ≤ sup ku(t)k + 1} t
with µ(v) as in (1.4).
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
5
The main result in this section is given in the following theorem: Theorem 2.1. Assume that k and h2r k −1 are sufficiently small. Then, we have the local stability estimate
(2.10)
n−1 n X |E1j | |||W n − U n |||λ ≤ emtn |||W 0 − U 0 |||λ + k j=0
1
+p k 2(1 − λ)
and the error estimate
n−1 X j=0
kE2j k2⋆
1/2 o
,
max |u(tn ) − U n | ≤ C(k + hr ).
(2.11)
0≤n≤N
Proof. Let ρn = un − W n and ϑn = W n − U n , n = 0, . . . , N. Then, according to (1.5), max |ρn | ≤ Chr .
(2.12)
0≤n≤N
Further, in view of (1.5) and (2.8), |||ϑ0 |||λ ≤ |ϑ0 | + k 1/2 kϑ0 k
≤ Chr + Ck 1/2 hr−d/2 ,
i.e., since d ≤ r, |||ϑ0 |||λ ≤ C(k + hr ) .
(2.13)
Now, if we assume that (2.10) holds, using (2.13), (2.6) and (2.7), we obtain max |||ϑn |||λ ≤ C(k + hr ) ,
(2.14)
0≤n≤N
and then (2.11) follows in view of (2.12). Thus, it remains to prove (2.10). To this end, we proceed as follows: From (2.1) and (2.3) we obtain an error equation for ϑn , ϑn+1 + kAh ϑn+1 =ϑn + k(Bh (W n ) − Bh (U n )) + kE n Z 1 n =ϑ + k Bh′ (W n − sϑn ) ϑn ds + kE n , 0
n = 0, . . . , N − 1.
According to (2.13), (2.6) and (2.7), there exists a constant C⋆ such that mT
(2.15)
e
n
0
0
|||W − U |||λ + k
N −1 X j=0
|E1j |
X j 1/2 o 1 k kE2 k2⋆ +p 2(1 − λ) j=0 N −1
≤ C⋆ (k + hr ) .
Next, we split the error ϑn as ϑn = ϑn1 +ϑn2 . Here ϑ01 = ϑ0 , ϑ02 = 0, and ϑni , n = 1, . . . , N, are recursively defined by Z 1 n+1 n+1 n (2.16) ϑi + kAh ϑi = ϑi + k Bh′ (W n − sϑn ) ϑni ds + kEin , i = 1, 2. 0
6
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
We shall show inductively that, for n = 0, . . . , N, % n−1 o n X 0 n mtn |||ϑ1 |||λ + k |E1j | (2.17) |||ϑ1 |||λ ≤ e j=0
and |||ϑn2 |||λ
(2.18)
mtn
≤e
1
p k 2(1 − λ)
n−1 X j=0
kE2j k2⋆
1/2
.
Obviously, (2.10) follows from (2.17) and (2.18). The estimates (2.17) and (2.18) are of course valid for n = 0. Assume that they hold for some n, 0 ≤ n ≤ N − 1. Then, according to (2.15) and to (2.10), which also holds for n, we have, for k and k −1 h2r small enough, 1 kϑn k ≤ √ C⋆ (k 1/2 + k −1/2 hr ) ≤ 1, λ i.e., sup µ (W n − sϑn ) ≤ m.
(2.19)
0<s 0, j = 1, . . . , s .
We also consider an explicit multistep scheme (α, γ) and we assume that both methods (α, β) and (α, γ) are of order p, i.e., (with the convention 00 = 1) q X
(iii)
ℓ
i αi = ℓ
q X
ℓ−1
i
i=0
i=0
βi = ℓ
q−1 X
iℓ−1 γi ,
ℓ = 0, 1, . . . , p .
i=0
An example of a class of (α, β, γ) methods satisfying the above assumptions with the order p = q is given by the polynomials q X 1 q−j α(ζ) = ζ (ζ − 1)j , j j=1
β(ζ) = ζ q ,
and γ(ζ) = ζ q − (ζ − 1)q .
The corresponding implicit (α, β) schemes are the well-known B.D.F. methods which are strongly A(0)−stable for 1 ≤ q ≤ 6. Remark 3.1. The hypothesis (iii) may be written in the equivalent form α(ex ) = xβ(ex ) + O(xp+1 ) = xγ(ex ) + O(xp+1 ),
as x → 0,
which implies β(y) − γ(y) = O((y − 1)p ),
(3.1)
as y → 1.
Since β −γ is a polynomial of degree q, we necessarily have p ≤ q ; recall that, cf. Cryer [5], the strong A(0)−stability of the (α, β) scheme implies also that p ≤ q . In the same paper Cryer (see also Grigorieff and Schroll [8]) shows that for any q there exists an (α, β) q-step A(0)−stable method of order q. Following the proof given in Hairer and Wanner, [9, Thm. 2.2, p. 270], it can be seen that these methods can be chosen to be strongly A(0)−stable. On the other hand, given an (α, β) method of order p = q, and since the degree of γ is ≤ q − 1, we deduce from (3.1) that the (α, γ) scheme will be of order q if and only if γ(ζ) = β(ζ) − βq (ζ − 1)q . In the sequel assume that we are given approximations U 0 , U 1 , . . . , U q−1 ∈ Vh to u0 , . . . , uq−1 such that (3.2)
q−1 X j=0
|W j − U j | + k 1/2 kW j − U j k ≤ c(k p + hr ).
We define U n ∈ Vh , n = q, . . . , N, recursively by the (α, β, γ) scheme (1.3). We shall prove in this section that the method (1.3) is stable, and we shall show convergence of
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
9
the approximations U n to u(tn ), as h and k tend to zero. In particular, under a mild meshcondition and for λ sufficiently small, we derive the optimal order error estimate max |u(tn ) − U n | ≤ c(k p + hr ).
0≤n≤N
As in the previous section, we shall first show consistency of the scheme (1.3) for the projection W, W (t) = Rh u(t), 0 ≤ t ≤ T . Consistency. The consistency error E n of the scheme ((1.3) for W is given by n
(3.3)
kE =
q X
αi W
n+i
+k
q X
βi Ah W
n+i
i=0
i=0
−k
q−1 X
γi Bh (W n+i ),
i=0
n = 0, . . . , N − q. Using the relations Ah Rh = Po A and Bh = Po B, we can rewrite (3.3) as n
kE = Rh
q X
αi u
n+i
+ kPo
q X
βi Au
n+i
i=0
i=0
− kPo
q−1 X
γi B(Rh un+i ),
i=0
and thus, in view of (1.1), we split E n as E n = E1n + E2n , where kE1n
=(Rh − Po )
(3.4i)
q X
αi u
n+i
+ Po
i=0
i=0
X
q X
q
+ kPo
βi B(u
n+i
i=0
)−
q−1 X
αi un+i − kβi u′ (tn+i ) γi B(un+i )
i=0
and (3.4ii)
E2n
= Po
q−1 X i=0
γi B(un+i ) − B(Rh un+i ) .
First, we will estimate E1n . Using (1.5) and the fact that α1 + · · · + αq = 0, we easily see that (3.5i)
q X αi un+i ≤ Ckhr . (Rh − Po ) i=0
Further, in view of the consistency properties of (α, β), (3.5ii)
q X n+i ′ α u − kβ u (t ) ≤ Ck p+1. i i n+i i=0
10
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
Moreover, q X
βi B(u
n+i
i=0
)−
q−1 X
γi B(un+i )
i=0
p−1 q−1 X X X (ik)ℓ ∂ ℓ (ik)ℓ ∂ ℓ γ B(u)(t ) − B(u)(tn ) + ϕn i n ℓ ℓ ℓ! ∂t ℓ! ∂t i=0 ℓ=0 ℓ=0 p−1
q
=
X
=
q X kℓ X
βi
i=0
p−1
ℓ=0 n
ℓ!
i=0
ℓ
βi i −
q−1 X
γi iℓ
i=0
∂ℓ
∂tℓ
B(u)(tn ) + ϕn
=ϕ , where the last equality holds in view of the consistency properties of (α, β) and (α, γ), and, obviously, |ϕn | ≤ ck p . This relation and (3.5i,ii) yield max |E1n | ≤ C(k p + hr ).
(3.6)
0≤n≤N −q
Finally, using (1.6), we obtain max kE2n k⋆ ≤ Chr ,
(3.7)
0≤n≤N −q
which completes the estimate of E n . Convergence of the multistep scheme. Let ϑn = W n − U n , n = 0, . . . , N. Then (3.3) and (1.3) yield the error equation for ϑn (3.8)
q X
n+i
αi ϑ
+k
q X
n+i
βi Ah ϑ
=k
i=0
i=0
i=0
q−1 X
γi {Bh (W n+i ) − Bh (U n+i )} + kE n .
According to the splitting of E n , we introduce ϑj1 and ϑj2 , cf. section 2, by
(3.9)
q X i=0
αi ϑjn+i
+k
q X
βi Ah ϑjn+i
=k
q−1 X i=0
i=0
γi
Z
1 0
Bh′ (W n+i − sϑn+i ) ds ϑjn+i
+ kEjn , j = 1, 2, n = 0, . . . , N − q, with initial values ϑj1 = ϑj and ϑj2 = 0 for j = 0, . . . , q−1. Then, obviously, ϑn = ϑn1 +ϑn2 . In the sequel we shall use the notation n n+q−1 Ej ϑj 0 αi + βi x , Θjn = ... , Ejn = .. , δi (x) = − αq + βq x . n ϑj 0 ∆i = δi (kAh ),
Γni
= γi
Z
0
1
Bh′ (W n+i − sϑn+i ) ds,
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
∆q−1 ∆q−2 . . . ∆0 I 0 0 Λ = Λ(kAh ) = , . . . . . . 0 I 0
and
11
q−1 Γn . . . Γn0 0 ... 0 Γn = . .. , . . . 0 ... 0
(αq + kβq Ah )ϑn+q−1 j .. (αq + kβq Ah ) Θjn = . . n (αq + kβq Ah )ϑj Equation (3.9) can then be written in the form
(3.10)
(αq + kβq Ah )Θjn+1 = (αq + kβq Ah ) ΛΘjn + k Γn Θjn + k Ejn ,
j = 1, 2.
We quote the following result from Crouzeix, [3]: Lemma 3.1. There exist a constant η, with 0 ≤ η < 1, and a continuous map H : ¯ + → Cq×q such that for all x ≥ 0 the matrix H(x) is invertible and the Euclidean R norm k · k2 of the matrix L(x), L(x) =
αq + βq x H(x)−1 Λ(x)H(x), αq + ηβq x
is less or equal to one for all x ≥ 0. Let H = H(kAh ),
and
L = L(kAh ),
Yjn = H−1 Θjn , Γen = H−1 Γn , Eejn = H−1 Ejn ;
then, we can rewrite (3.10) as (3.11)
(αq + kβq Ah )Yjn+1 = (αq + kηβq Ah ) LYjn + k Γen HYjn + k Eejn ,
j = 1, 2.
In view of the fact that kH(x)k2 , kH(x)−1 k2 are uniformly bounded, see relations (3.27) and (3.28) in [3], it suffices to estimate Y n . We adjust in this section the definition of ||| · ||| to the scheme under consideration by setting |||v||| := (αq |v|2 + βq kkvk2 )1/2 ,
v ∈ V.
Further, for V = (v1 , . . . , vq )T and W = (w1 , . . . , wq )T in H q or in V q we shall use the notation q q X 1/2 X (vi , wi ), |V | := |vi |2 (V, W ) := , i=1
X q
kV k :=
i=1
kvi k
2
1/2
i=1
X q
,
|||V ||| :=
i=1
|||vi |||
2
1/2
,
kV k⋆ :=
q X i=1
and, for a linear operator M : H q → H q , we set |M| := supV ∈H q ,V 6=0 The main result in this paper is given in the following theorem:
kvi k2⋆
|M V | . |V |
1/2
,
12
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
Theorem 3.1. Assume that the constant λ in (1.4) is appropriately small (depending on the particular scheme) and that k and h2r k −1 are sufficiently small. Then, we have the local stability estimate |ϑn | + k 1/2 kϑn k ≤ CecmT (3.12)
q−1 nX j=0
n−q X |E1j | |ϑj | + k 1/2 kϑj k + k j=0
n−q
+ k
X
kE2j k2⋆
j=0
1/2 o
,
n = q − 1, . . . , N,
and the error estimate max |u(tn ) − U n | ≤ C(k p + hr ).
(3.13)
0≤n≤N
Proof. Let ρn = un − W n , n = 0, . . . , N. Then, according to ((1.5), max |ρn | ≤ Chr .
(3.14)
0≤n≤N
Now, if we assume that (3.12) holds, using (3.2), (3.6) and (3.7), we obtain max |ϑn | ≤ C(k p + hr ) ,
(3.15)
0≤n≤N
and (3.13) follows immediately from (3.14) and (3.15). Thus, it remains to prove (3.12). According to (3.2), (3.6) and (3.7)), there exists a constant C⋆ such that the right-hand side of (3.12) can be estimated by C⋆ (k p + hr ), cmT
(3.16)
Ce
q−1 nX j=0
j
|ϑ | + k
1/2
j
kϑ k + k
n−q X j=0
|E1j |
+ k
N −q
X j=0
kE2j k2⋆
1/2 o
≤ C⋆ (k p + hr ) . We will estimate ϑn by estimating Yjn . In fact, we shall show that for some ε = ε(λ, (α, β, γ) ), 0 < ε < (1 − η 2 )βq with η as in Lemma 3.1, (3.17)
n−1 n k X ej o |||Y1n ||| ≤ ecmtn |||Y10 ||| + √ |E | αq j=0 1
and n−1
(3.18)
|||Y2n |||
cmtn
≤e
X j 1/2 1 √ k kEe2 k2⋆ . ε j=0
Then, (3.12) follows, and the proof will be complete. We shall use induction: The estimates (3.17) and (3.18) are valid for n = 0. Assume that they hold for 0, . . . , n, 0 ≤ n ≤ N −q. Then, according to (3.16) and (3.12), which is then valid for 0, . . . , n+ q −1, we have, for k and k −1 h2r small enough, max
0≤j≤n+q−1
kϑj k ≤ C⋆ (k p−1/2 + k −1/2 hr ) ≤ 1,
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
13
i.e., (3.19)
sup
max
0<s 0 we seek a realvalued function u defined on R × [0, T ], 1-periodic in the space variable and satisfying (4.1)
ut + uux + uxx + νuxxxx = 0 in R × [0, T ]
and (4.2)
u(·, 0) = u0
in R,
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
15
where u0 is a given, smooth 1-periodic function. The periodic initial value problem (4.1)–(4.2) is well-posed, cf. [13], [17]. For numerical methods to this problem see, e.g., [1] and the references therein. The KS equation is related to turbulence phenomena in chemistry and combustion, and arises also in plasma physics and in two-phase flows in cylindrical geometries, see [18] and [14]. s For s ∈ N0 , let Hper denote the periodic Sobolev space of order s, consisting of the s . The 1-periodic elements of H sloc (R), and let k · kH s be the norm over a period in Hper 2 0 inner product in H := Lper = Hper is denoted by (·, ·), and the induced norm by | · |. 2 Let A : H 4per → H be defined by Av = ν(vxxxx + v). Then V := D(A1/2 ) = Hper , and 1/2 2 2 1/2 the norm in V is given by kvk = ν (|vxx | + |v| ) . Let B : V → H be given by B(v) = −vxx − vvx + νv. Then, B ′ (v)w = −wxx − (vw)x + νw , and thus by periodicity (B ′ (v)w, ω) = (wx , ωx ) + (v, wωx ) + ν(w, ω) . 1 , Therefore, in view of the inequality kwkL∞ ≤ |w| + |w ′|, w ∈ Hper
|(B ′ (v)w, ω)| ≤|wx | |ωx| + |v|kwkL∞ |ωx | + ν|w| |ω|
≤(1 + |v|)|wx | |ωx| + ν|w| |ω| + |v| |w| |ωx|,
and thus, using the inequality |ux |2 ≤ |u| |uxx|, u ∈ V, we easily see that the condition (1.4) is satisfied for any λ > 0, (4.3)
|(B ′ (v)w, ω)| ≤ λkwkkωk + µ(v)|w| |ω|,
∀v, w, ω ∈ V,
1 [ 1 + 2λν 2 + 2|v|(1 + |v|) ] . We note that, since λ can be taken arbiwith µ(v) := 2λν trarily small, our results hold for this problem for all implicit-explicit schemes satisfying our stability and consistency assumptions; further, since µ is bounded for v bounded in H, the meshcondition of Theorems 2.1 and 3.1 is not needed here, cf. Remark 2.1. For the space discretization, we let 0 = x0 < x1 < · · · < xJ = 1 be a partition of [0, 1], and h := maxj (xj+1 − xj ). Setting xjJ+s := j + xs , j ∈ Z, s = 0, . . . , J − 1, this partition is periodically extended to a partition of R. For integer r ≥ 4, let Vh denote a space of at least once continuously differentiable, 1-periodic splines of degree r − 1, in which approximations to the solution u(·, t) of (4.1)–(4.2) will be sought for 0 ≤ t ≤ T . The following approximation property of the family {Vh }0 0 we seek a real-valued function u, defined on ¯ × [0, T ], satisfying Ω ut − ∆u = f (u) in Ω × [0, T ], u=0
(4.18)
on ∂Ω × [0, T ],
u(·, 0) = u0
in Ω,
where u0 is a given smooth function and f : R → R is a smooth function. We are interested in approximating smooth solutions of this problem, and assume therefore that the data are smooth and compatible such that (4.18) gives rise to sufficiently regular solutions. We shall distinguish two cases. Assuming that f satisfies a polynomial growth condition of order at most ρ, see (4.19) below, and provided that ρ ≤ 4 for ν = 3, we show that the abstract theory of sections 2 and 3 is directly applicable and yields, without any additional assumptions on the discretization spaces, optimal order error estimates for all (α, β, γ) schemes considered in this paper. For general f our results apply as well, provided that meshconditions stronger than those of Theorem 3.1 are fulfilled. Reaction diffusion equations model various physical phenomena related, for instance, to phase transitions, chemical reactions, pattern formation, cf., e.g., [2], [7], [10] and their references. The Allen-Cahn model (ρ = 3), which in the limit describes evolution by mean curvature, [2], [7], and generalized Ginzburg-Landau equations, [10], are reaction diffusion equations of particular interest with polynomial nonlinearity. Let H s = H s (Ω) be the usual Sobolev space of order s, and k · kH s be the norm of s H . The inner product in H := L2 (Ω) is denoted by (·, ·), and the induced norm by | · |; the norm of Ls (Q), 1 ≤ s ≤ ∞, is denoted by k · kLs (Q) and simply by k · kLs for Q = Ω. Obviously, V = H01 (Ω) = H01 and the norm in V is equivalent to the H 1 norm. We now consider the case that f satisfies the following growth condition: there exists L ∈ R such that |f ′ (ξ)| ≤ L(1 + |ξ|ρ−1) ∀ξ ∈ R .
(4.19)
In the sequel, we shall use the Sobolev inequality (4.20)
kvkLs ≤ CkvkH 1 ,
1 ≤ s < ∞ for ν = 2,
and 1 ≤ s ≤ 6 for ν = 3 ,
as well as its consequence (4.21)
1−ν/4
kvkL4 ≤ CkvkL2
ν/4
kvkH 1 ,
ν = 2, 3.
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
19
Let B : V → V ′ , B(v) = f (v). First, we note that B is well defined. Indeed, for v, w ∈ V, Z 1 (f (v), w) = f ′ (sv)v ds, w + (f (0), w), 0
and, therefore, by (4.19) and H¨older’s inequality, i hZ 5/6 ρ6/5 (4.22) |(f (v), w)| ≤ C |v| dx kwkL6 + (1 + kvkL2 ) kwkL2 , Ω
and thus, in view of (4.20), we see that f (v) ∈ V ′ , provided ρ ≤ 5 for ν = 3. Further, by (4.19), for v, w, ω ∈ V, Z ′ |(f (v)w, ω)| ≤ C |v|ρ−1|w||ω|dx + CkwkL2 kωkL2 Ω Z 1/2 2(ρ−1) ≤C |v| dx kwkL4 kωkL4 + CkwkL2 kωkL2 , Ω
i.e., in view of (4.21), (4.23)
|w| |(f ′ (v)w, ω)| ≤ Ckvkρ−1 L2(ρ−1)
4−ν 4
|ω|
4−ν 4
ν
ν
kwk 4 kωk 4 + C|w| |ω| .
Thus, B is differentiable. Furthermore, (4.23) and Young’s inequality (ab ≤ qε aq + ′
ε−q /q q ′ 1 b , q q′
(4.24)
+
1 q′
= 1 ) yield
|(B ′ (v)w, ω)| ≤ C ν4 εkwk kωk 4(ρ−1) ν 4−ν − 4−ν 4−ν +C 1+ 4 ε |w| |ω| . kvkL2(ρ−1)
Therefore, by Sobolev’s inequality, we see that, for ρ ≤ 4 when ν = 3 (and for any ρ when ν = 2), the assumption (1.4) is satisfied with λ = C ν4 ε, µ(v) = C + C
4−ν 4
ν
ε− 4−ν kvkL2(ρ−1)
4(ρ−1) 4−ν
,
and µ(v) is bounded for v bounded in V. Again, λ can be taken arbitrarily small and thus our theory applies to all (α, β, γ) schemes satisfying our stability and consistency assumptions. We further analyze the case that f satisfies the growth condition (4.19), with ρ ≤ 4 for ν = 3, by introducing the discretization spaces Vh ; the case of general f shall be discussed at the end of this section. For simplicity, let Vh be the subspace of V defined on a finite element partition Th of Ω, and consisting of piecewise polynomial functions of degree at most r − 1, r ≥ 2. Let hK denote the diameter of an element K ∈ Th , and h = maxK∈Th hK . We define the elliptic projection operator Rh , Rh : V → Vh , by (∇Rh v, ∇χ) = (∇v, ∇χ) ∀χ ∈ Vh . We assume that (we do not attempt here to deal with problems that may arise in the case of a curved boundary ∂Ω concerning the requirement Vh ⊂ V, cf. Remark 2.2) (4.25)
|v − Rh v| + hkv − Rh vk ≤ Chr kvkH r ,
v ∈ H r ∩ H01 ;
20
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
then, in particular, the estimate (1.5)) will hold in this case with d = 2. Next, we will verify (1.6). We shall further assume that (4.26)
sup{kRh vkL∞ : 0 < h < 1} ≤ CkvkH r ,
For v ∈ H r ∩ H01 , we have
B(v) − B(Rh v) =
Z
v ∈ H r ∩ H01 .
1
f ′ (v − τ (v − Rh v))dτ (v − Rh v);
0
using here the Sobolev inequality kvkL∞ ≤ CkvkH r (r ≥ 2) and (4.26), we conclude, in view of (4.25), ˜ |B(v) − B(Rh v)| ≤ Chr M(v)kvk Hr ,
(4.27)
˜ with M(v) bounded for v ∈ H r ∩ H01 bounded in H r . Thus, (1.6) is satisfied. Now, let W (t) := Rh u(t), and assume that we are given approximations U 0 , . . . , U q−1 ∈ Vh to u0 , . . . , uq−1 such that q−1 X
(4.28)
j=0
|W j − U j | + k 1/2 kW j − U j k ≤ c(k p + hr ).
Then, we define U n ∈ Vh , n = q, . . . , N, recursively by the (α, β, γ) scheme
(4.29)
q X
αi (U n+i ,χ) + k
q X i=0
i=0
βi (∇U n+i , ∇χ)
q−1
=k
X i=0
γi (f (U n+i ), χ) ∀χ ∈ Vh ,
n = 0, . . . , N − q,
where (α, β) and (α, γ) are multistep schemes of order p, and (α, β) is strongly A(0)−stable. Then, Theorem 3.1 yields, for sufficiently small k and k −1 h2r , the error estimate (4.30)
max |un − U n | ≤ c(k p + hr ). n
Remark 4.2. By (4.24), we see that (1.4) holds for any λ > 0 with µ(v) bounded for v bounded in L2(ρ−1) . Using this fact, Sobolev’s inequality kvkLs ≤ C|v|1−akvka ,
a = ν s−2 , ν = 2, 3, (with s ≤ 6, of course, for ν = 3) and Remark 2.1, it is easily seen 2s that the meshcondition k −1 h2r ≤ c, under which (4.30) holds, can be weakened. We shall not dwell on this. We close this subsection by briefly considering the case of a general smooth function f, as well as the case that f satisfies (4.19) for ν = 3 but with ρ > 4. First, we note that in our analysis it suffices to assume that B is well defined and differentiable on a subspace Ve of V ∩ L∞ containing Vh , for all h. By tracing back the proof of (4.24) we see that in this case µ(v) is bounded, provided that Z |f ′ (v(x))|2 dx Ω
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
21
is bounded. Note that the assumption k −1 h2r ≤ c, for appropriate c, of Theorem 3.1 is only used to show that kϑn k ≤ 1, which implies (3.19), i.e., that µ(W j − sϑj ), s ∈ (0, 1), is bounded by a constant. In the case under investigation, by using appropriate inverse inequalities, we show that if stronger meshconditions are R ′ satisfied, then sup0<s 4, and ν = 3 and general f. i. ν = 2 and general f. First, we note that ∀χ ∈ Vh , kχkL∞ ≤ C| log(h)|1/2 kχkH 1 R with h = minK∈Th hK , cf. [19, p. 67]. Obviously, Ω |f ′ (χ(x))|2 dx is bounded if kχkL∞ is bounded. Now, max
0≤j≤n+q−1
kϑj kL∞ ≤ C| log(h)|1/2
max
0≤j≤n+q−1
kϑj k,
and thus, according to (3.16), max
0≤j≤n+q−1
kϑj kL∞ ≤ C⋆ C| log(h)|1/2 (k p−1/2 + k −1/2 hr ).
Therefore, if k and h are chosen such that | log(h)|k 2p−1 and | log(h)|k −1 h2r are sufficiently small, then µ(W j − sϑj ) will be bounded, and the convergence results hold. ii. ν = 3 and f satisfies (4.19) with ρ > 4. If s ≥ 6, we have kχkLs (Ω) ≤ Ch−
(4.31)
s−6 2s kχkH 1 (Ω)
∀χ ∈ Vh .
Indeed, employing standard homogeneity arguments, one can show that for an arbitrary element K ∈ Th , s−6 2s
−
kχkLs (K) ≤ ChK
k∇χkL2 (K) + kχkL6 (K) ,
and (4.31) follows in view of (4.20). Hence, max
0≤j≤n+q−1
ρ−4 − 2(ρ−1)
kϑj kL2(ρ−1) ≤ Ch −
Therefore, if k and h are such that h then, as before, max
0≤j≤n+q−1
ρ−4 ρ−1 k 2p−1
max
0≤j≤n+q−1 −
and k −1 h
kϑj k.
ρ−4 ρ−1 h2r
are sufficiently small,
ρ−4
kϑj kL2(ρ−1) ≤ C⋆ h− 2(ρ−1) (k p−1/2 + k −1/2 hr ) ≤ 1,
µ(W j − sϑj ) will be bounded in view of (4.24) and (4.26), and our convergence results hold. iii. ν = 3 and general f. In this case, kχkL∞ ≤ Ch−1/2 kχkH 1
∀χ ∈ Vh ,
22
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
as one can see by modifying the proof of the two dimensional case. By the same arguments as for ν = 2, the convergence results of this paper hold, provided that h−1 k 2p−1 and k −1 h−1 h2r are sufficiently small.
References 1. G. Akrivis, High-order finite element methods for the Kuramoto-Sivashinsky equation, RAIRO Mod´el. Math. Anal. Num´er. 30 (1996), 157–183. MR 97e:65095 2. S.M. Allen and J.W. Cahn, A macroscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall. 27 (1979), 1085-1095. 3. M. Crouzeix, Une m´ethode multipas implicite-explicite pour l’approximation des ´equations d’´evolution paraboliques, Numer. Math. 35 (1980), 257–276. MR 82b:65084 4. M. Crouzeix and P.-A. Raviart, Approximation des ´equations d’´evolution lin´eares par des m´ethodes a pas multiples, C.R. Acad. Sc. Paris, S´erie A 283 (1976), 367–370. MR 54:14377 ` 5. C. W. Cryer, A new class of highly stable methods, BIT 13 (1973), 153–159. MR 48:1469 6. C.M. Elliott and S.-M. Zheng, On the Cahn-Hilliard equation, Arch. Rational Mech. Anal. 96 (1986), 339–357. MR 87k:80007 7. L.C. Evans, H.M. Soner and P.E. Souganidis, Phase transitions and generalized motion by mean curvature, Comm. Pure Appl. Math. 45 (1992), 1097–1123. MR 93g:35064 ¨ 8. R.D. Grigorieff and J. Schroll, Uber A(α)-stabile Verfahren hoher Konsistenzordnung, Computing 20 (1978), 343–350. MR 83b:65086 9. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, Springer-Verlag, Berlin Heidelberg, Springer Series in Computational Mathematics v. 14, 1991. MR 92a:65016 10. Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer-Verlag, New York, Series in Synthetics, v. 19, 1984. MR 87e:92054 11. M.N. Le Roux, Semi-discr´etisation en temps pour les ´equations d’´evolution paraboliques lorsque l’op´erateur d´epend du temps, RAIRO 13 (1979), 119–137. MR 80e:65056 12. W.R. McKinney, Optimal error estimates for high order Runge-Kutta methods applied to evolutionary equations, Ph.D. thesis, University of Tennessee, Knoxville, 1989. 13. B. Nicolaenko and B. Scheurer, Remarks on the Kuramoto-Sivashinsky equation, Physica 12 D (1984), 391–395. MR 86d:80007 14. D.T. Papageorgiou, C. Maldarelli and D.S. Rumschitzki, Nonlinear interfacial stability of coneannular film flow, Phys. Fluids A2 (1990), 340–352. MR 91b:76067 15. G. Savar´e, A(Θ)-stable approximations of abstract Cauchy problems, Numer. Math. 65 (1993), 319–336. MR 94h:65062 16. L.L. Schumaker, Spline Functions: Basic Theory, Wiley, New York, 1981. MR 82j:41001 17. E. Tadmor, The well-posedness of the Kuramoto-Sivashinsky equation, SIAM J. Math. Anal. 17 (1986), 884–893. MR 87g:35117 18. R. Temam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics, Springer-Verlag, New York, 1988. MR 89m:58056 19. V. Thom´ee, Galerkin finite element methods for parabolic problems, Springer-Verlag, Lecture Notes in Mathematics v. 1054, 1984. MR 86k:65006 20. M. Zl´amal, Finite element multistep discretizations of parabolic boundary value problems, Math. Comp. 29 (1975), 350–359. MR 51:7326 21. M. Zl´amal, Finite element methods for nonlinear parabolic equations, RAIRO 11 (1977), 93–107. MR 58:19239
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
23
Department of Computer Science, University of Ioannina, 451 10 Ioannina, Greece E-mail address: [email protected] IRMAR, Universit´ e de Rennes I, Campus de Beaulieu, F-35042 Rennes, France E-mail address: [email protected] Department of Mathematics, University of Crete, 714 09 Heraklion, Crete, Greece, and IACM, Foundation for Research and Technology - Hellas, 711 10 Heraklion, Crete, Greece E-mail address: [email protected]
1. Introduction In this paper we construct and analyze implicit-explicit multistep schemes for the time discretization of a class of nonlinear parabolic problems of the form: Given T > 0 and u0 ∈ H, find u : [0, T ] → D(A) such that (1.1)
u′ (t) + Au(t) = B(u(t)),
0 < t < T,
u(0) = u0 ,
where A is a linear, selfadjoint, positive definite operator on a Hilbert space (H, (·, ·)) with domain D(A) dense in H, and B : D(A) → H is a (possibly nonlinear) differentiable operator. Considering a finite dimensional subspace Vh of V, V := D(A1/2 ), we are led to a semidiscrete problem approximating (1.1): We seek a function uh , uh (t) ∈ Vh , defined by (1.2)
u′h (t) + Ah uh (t) = Bh (uh (t)),
0 < t < T,
uh (0) = u0h ;
here u0h ∈ Vh is a given approximation to u0 , and Ah , Bh are appropriate operators on Vh with Ah linear, selfadjoint and positive definite. Following an idea of Crouzeix, [3], for the time discretization of parabolic equations with time dependent coefficients, we combine implicit and explicit multistep schemes to discretize (1.2) in time: Implicit schemes are used for discretizing the left-hand side of the o.d.e. in (1.2), and explicit schemes for the nonlinear right-hand side. Thus, letting 1991 Mathematics Subject Classification. Primary 65M60, 65M12; Secondary 65L06. The work of the first and third authors was supported in part by a research grant from the University of Crete. 1
2
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
k be a (constant) time step, tn = nk, n = 0, . . . , N, T = Nk, we define a sequence of approximations U n , U n ∈ Vh , to un = u(tn ), by (1.3)
q X
αi U
n+i
+k
q X
βi Ah U
n+i
=k
γi Bh (U n+i ).
i=0
i=0
i=0
q−1 X
A scheme of this form is characterized by three polynomials α, β and γ, α(ζ) =
q X i=0
i
αi ζ ,
β(ζ) =
q X
i
βi ζ ,
i=0
γ(ζ) =
q−1 X
γi ζ i .
i=0
We note that, when Bh vanishes, the scheme (1.3) reduces to the implicit linear multistep method (α, β) (or (ρ, σ) in the notation of Dahlquist) for the equation u′h (t)+Ah uh (t) = 0; similarly if Ah vanishes, the scheme (1.3) reduces to the linear multistep method (α, γ), which is explicit since γq = 0, for the equation u′h (t) = Bh (uh (t)). The scheme (1.3) is a combination of the methods (α, β) and (α, γ); it is linearly implicit and nonlinearly explicit. We shall refer to it as the (α, β, γ) scheme. We shall assume in the sequel that the method (α, β) is strongly A(0)−stable; this implies, in particular, that αq βq is positive, which in turn ensures invertibility of the operator αq I + kβq Ah . Thus, given U 0 , . . . , U q−1 in Vh , U q , . . . , U N are well defined by the (α, β, γ) scheme. We will assume in the sequel, without loss of generality, that both αq and βq are positive. For the analysis of the scheme (1.3) we will need some additional assumptions for the operators A and B as well as for the finite dimensional spaces Vh ; the operators Ah and Bh will then be appropriately defined. Let, thus, | · | denote the norm of H, and introduce in V the norm k · k by kvk = (A1/2 v, A1/2 v)1/2 . We identify H with its dual, and denote by V ′ the dual of V , again by (·, ·) the duality pairing on V ′ and V, and by k · k⋆ the dual norm on V ′ . We assume that B can be extended to an operator from V into V ′ , which is differentiable, and an estimate of the form (1.4)
|(B ′ (v)w, ω)| ≤ λkwkkωk + µ(v)|w| |ω|
∀v, w, ω ∈ V
holds with a constant λ < 1 and a functional µ(v) which is bounded for v bounded in V. Indeed, depending on the particular (α, β, γ) scheme, we shall need to assume that the constant λ is appropriately small. Further, the assumption that µ(v) is bounded for v bounded in V will suffice to derive our results under some mild meshconditions; these conditions can be weakened if µ(v) is bounded for v bounded in some appropriate weaker norms, and even avoided if µ(v) is bounded for v bounded in H. We will assume in the sequel that (1.1) possesses a solution which is sufficiently regular for our results to hold. Uniqueness of smooth solutions follows easily in view of (1.4). For the space discretization we use a family Vh , 0 < h < 1, of finite dimensional subspaces of V. In the sequel the following discrete operators will play an essential role: Define Po : V ′ → Vh , Rh : V → Vh , Ah : V → Vh , and the nonlinear operator
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
3
Bh : V → Vh by (Po v, χ) = (v, χ) ∀χ ∈ Vh ,
(ARh v, χ) = (Av, χ) ∀χ ∈ Vh ,
(Ah ϕ, χ) = (Aϕ, χ) ∀χ ∈ Vh ,
(Bh (ϕ), χ) = (B(ϕ), χ) ∀χ ∈ Vh . Then, obviously, Ah Rh = Po A and Bh = Po B . In the error analysis we shall use the approximation properties of the “elliptic projection” Rh . We assume that Rh v is an approximation to v of order r, provided that v is sufficiently regular, (1.5)
|v − Rh v| + hd/2 kv − Rh vk ≤ M(v)hr ,
where r and d are two integers, 2 ≤ d ≤ r, and M(v) is bounded if an appropriate norm of v is bounded. We further assume that (1.6)
kB(v) − B(Rh v)k⋆ ≤ M(v)hr .
We emphasize that the condition (1.6) serves consistency purposes only. It is needed to prove consistency of the (α, β, γ) scheme—and, in fact, already of the semidiscrete problem (1.2)—for Rh u of optimal order with respect to h. It is somehow restrictive though, since it essentially means that, if A is a differential operator of order d, then B may contain derivatives of up to order d/2 only. For some concrete differential equations, however, one can get by with a less stringent condition by taking into account in the definition of Rh the terms of B of order higher than d/2. The scheme (1.3) is very efficient since its implementation requires at every time step solving a linear system with the same matrix for all time levels. If the order of both the implicit and the explicit scheme is p, then under some mild meshconditions and for appropriately small λ, and appropriate starting values U 0 , . . . , U q−1 , we derive the optimal order error estimate max |u(tn ) − U n | ≤ c(k p + hr ).
0≤n≤N
An outline of the paper is as follows: Section 2 is devoted to the analysis of a simple one step semiexplicit scheme of first order accuracy; its purpose is to motivate the analysis, in section 3, of more general multistep schemes of higher accuracy. Finally, in the last section, we apply our abstract results to three examples, namely the KuramotoSivashinsky equation and the Cahn-Hilliard equation in one space dimension, and to a class of reaction diffusion equations in Rν , ν = 2, 3. 2. An implicit-explicit one step scheme As a motivation for the analysis of multistep schemes, we study in this section the simplest implicit-explicit scheme which is a combination of the backward Euler and the Euler method.
4
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
Let U 0 ∈ Vh be given. We define fully discrete approximations U n ∈ Vh to u(tn ), recursively by U n+1 − U n + Ah U n+1 = Bh (U n ), n = 0, . . . , N − 1. k Our main concern in this section is to analyze the approximation properties of the sequence {U n }. As an intermediate step, we shall show consistency of the scheme (2.1) for the elliptic projection of the solution u of (1.1). Let (2.1)
W (t) = Rh u(t),
W (t) ∈ Vh , 0 ≤ t ≤ T .
We note for later use that, in view of the definition of Rh , kW (t)k is obviously bounded by ku(t)k, and thus bounded uniformly in h and t, (2.2)
sup sup kW (t)k ≤ C . h
t
Consistency. The consistency error E n of the scheme (2.1) for W is given by (2.3)
kE n = (W n+1 − W n ) + kAh W n+1 − kBh (W n ),
n = 0, 1, . . . , N − 1.
Using the relation Ah Rh = Po A and the definition of Bh , we rewrite E n as kE n = Rh (un+1 − un ) + kPo Aun+1 − kPo B(Rh un ),
and thus, in view of (1.1), E n = E1n + E2n , where (2.4)
kE1n =(Rh − Po )(un+1 − un ) + Po (un+1 − un − ku′ (tn+1 ) )
and
+ kPo (B(un+1) − B(un ) ),
E2n = Po (B(un ) − B(Rh un ) ).
(2.5) Obviously
max |E1n | ≤ C(k + hr ).
(2.6)
0≤n≤N −1
Further, in view of (1.6), (2.7)
max kE2n k⋆ ≤ Chr .
0≤n≤N −1
Next, we show optimal order rate of convergence of our approximations to the (sufficiently regular) solution of (1.1), provided that the initial approximation U 0 ∈ Vh satisfies (2.8)
|u0 − U 0 | + hd/2 ku0 − U 0 k ≤ c(u)hr .
We first introduce some more notation: For v ∈ V and b > 0, let 1/2 |||v|||b := |v|2 + bkkvk2 , and set |||v||| := |||v|||1 . Further, let
(2.9)
m := sup{µ(v) : kvk ≤ sup ku(t)k + 1} t
with µ(v) as in (1.4).
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
5
The main result in this section is given in the following theorem: Theorem 2.1. Assume that k and h2r k −1 are sufficiently small. Then, we have the local stability estimate
(2.10)
n−1 n X |E1j | |||W n − U n |||λ ≤ emtn |||W 0 − U 0 |||λ + k j=0
1
+p k 2(1 − λ)
and the error estimate
n−1 X j=0
kE2j k2⋆
1/2 o
,
max |u(tn ) − U n | ≤ C(k + hr ).
(2.11)
0≤n≤N
Proof. Let ρn = un − W n and ϑn = W n − U n , n = 0, . . . , N. Then, according to (1.5), max |ρn | ≤ Chr .
(2.12)
0≤n≤N
Further, in view of (1.5) and (2.8), |||ϑ0 |||λ ≤ |ϑ0 | + k 1/2 kϑ0 k
≤ Chr + Ck 1/2 hr−d/2 ,
i.e., since d ≤ r, |||ϑ0 |||λ ≤ C(k + hr ) .
(2.13)
Now, if we assume that (2.10) holds, using (2.13), (2.6) and (2.7), we obtain max |||ϑn |||λ ≤ C(k + hr ) ,
(2.14)
0≤n≤N
and then (2.11) follows in view of (2.12). Thus, it remains to prove (2.10). To this end, we proceed as follows: From (2.1) and (2.3) we obtain an error equation for ϑn , ϑn+1 + kAh ϑn+1 =ϑn + k(Bh (W n ) − Bh (U n )) + kE n Z 1 n =ϑ + k Bh′ (W n − sϑn ) ϑn ds + kE n , 0
n = 0, . . . , N − 1.
According to (2.13), (2.6) and (2.7), there exists a constant C⋆ such that mT
(2.15)
e
n
0
0
|||W − U |||λ + k
N −1 X j=0
|E1j |
X j 1/2 o 1 k kE2 k2⋆ +p 2(1 − λ) j=0 N −1
≤ C⋆ (k + hr ) .
Next, we split the error ϑn as ϑn = ϑn1 +ϑn2 . Here ϑ01 = ϑ0 , ϑ02 = 0, and ϑni , n = 1, . . . , N, are recursively defined by Z 1 n+1 n+1 n (2.16) ϑi + kAh ϑi = ϑi + k Bh′ (W n − sϑn ) ϑni ds + kEin , i = 1, 2. 0
6
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
We shall show inductively that, for n = 0, . . . , N, % n−1 o n X 0 n mtn |||ϑ1 |||λ + k |E1j | (2.17) |||ϑ1 |||λ ≤ e j=0
and |||ϑn2 |||λ
(2.18)
mtn
≤e
1
p k 2(1 − λ)
n−1 X j=0
kE2j k2⋆
1/2
.
Obviously, (2.10) follows from (2.17) and (2.18). The estimates (2.17) and (2.18) are of course valid for n = 0. Assume that they hold for some n, 0 ≤ n ≤ N − 1. Then, according to (2.15) and to (2.10), which also holds for n, we have, for k and k −1 h2r small enough, 1 kϑn k ≤ √ C⋆ (k 1/2 + k −1/2 hr ) ≤ 1, λ i.e., sup µ (W n − sϑn ) ≤ m.
(2.19)
0<s 0, j = 1, . . . , s .
We also consider an explicit multistep scheme (α, γ) and we assume that both methods (α, β) and (α, γ) are of order p, i.e., (with the convention 00 = 1) q X
(iii)
ℓ
i αi = ℓ
q X
ℓ−1
i
i=0
i=0
βi = ℓ
q−1 X
iℓ−1 γi ,
ℓ = 0, 1, . . . , p .
i=0
An example of a class of (α, β, γ) methods satisfying the above assumptions with the order p = q is given by the polynomials q X 1 q−j α(ζ) = ζ (ζ − 1)j , j j=1
β(ζ) = ζ q ,
and γ(ζ) = ζ q − (ζ − 1)q .
The corresponding implicit (α, β) schemes are the well-known B.D.F. methods which are strongly A(0)−stable for 1 ≤ q ≤ 6. Remark 3.1. The hypothesis (iii) may be written in the equivalent form α(ex ) = xβ(ex ) + O(xp+1 ) = xγ(ex ) + O(xp+1 ),
as x → 0,
which implies β(y) − γ(y) = O((y − 1)p ),
(3.1)
as y → 1.
Since β −γ is a polynomial of degree q, we necessarily have p ≤ q ; recall that, cf. Cryer [5], the strong A(0)−stability of the (α, β) scheme implies also that p ≤ q . In the same paper Cryer (see also Grigorieff and Schroll [8]) shows that for any q there exists an (α, β) q-step A(0)−stable method of order q. Following the proof given in Hairer and Wanner, [9, Thm. 2.2, p. 270], it can be seen that these methods can be chosen to be strongly A(0)−stable. On the other hand, given an (α, β) method of order p = q, and since the degree of γ is ≤ q − 1, we deduce from (3.1) that the (α, γ) scheme will be of order q if and only if γ(ζ) = β(ζ) − βq (ζ − 1)q . In the sequel assume that we are given approximations U 0 , U 1 , . . . , U q−1 ∈ Vh to u0 , . . . , uq−1 such that (3.2)
q−1 X j=0
|W j − U j | + k 1/2 kW j − U j k ≤ c(k p + hr ).
We define U n ∈ Vh , n = q, . . . , N, recursively by the (α, β, γ) scheme (1.3). We shall prove in this section that the method (1.3) is stable, and we shall show convergence of
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
9
the approximations U n to u(tn ), as h and k tend to zero. In particular, under a mild meshcondition and for λ sufficiently small, we derive the optimal order error estimate max |u(tn ) − U n | ≤ c(k p + hr ).
0≤n≤N
As in the previous section, we shall first show consistency of the scheme (1.3) for the projection W, W (t) = Rh u(t), 0 ≤ t ≤ T . Consistency. The consistency error E n of the scheme ((1.3) for W is given by n
(3.3)
kE =
q X
αi W
n+i
+k
q X
βi Ah W
n+i
i=0
i=0
−k
q−1 X
γi Bh (W n+i ),
i=0
n = 0, . . . , N − q. Using the relations Ah Rh = Po A and Bh = Po B, we can rewrite (3.3) as n
kE = Rh
q X
αi u
n+i
+ kPo
q X
βi Au
n+i
i=0
i=0
− kPo
q−1 X
γi B(Rh un+i ),
i=0
and thus, in view of (1.1), we split E n as E n = E1n + E2n , where kE1n
=(Rh − Po )
(3.4i)
q X
αi u
n+i
+ Po
i=0
i=0
X
q X
q
+ kPo
βi B(u
n+i
i=0
)−
q−1 X
αi un+i − kβi u′ (tn+i ) γi B(un+i )
i=0
and (3.4ii)
E2n
= Po
q−1 X i=0
γi B(un+i ) − B(Rh un+i ) .
First, we will estimate E1n . Using (1.5) and the fact that α1 + · · · + αq = 0, we easily see that (3.5i)
q X αi un+i ≤ Ckhr . (Rh − Po ) i=0
Further, in view of the consistency properties of (α, β), (3.5ii)
q X n+i ′ α u − kβ u (t ) ≤ Ck p+1. i i n+i i=0
10
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
Moreover, q X
βi B(u
n+i
i=0
)−
q−1 X
γi B(un+i )
i=0
p−1 q−1 X X X (ik)ℓ ∂ ℓ (ik)ℓ ∂ ℓ γ B(u)(t ) − B(u)(tn ) + ϕn i n ℓ ℓ ℓ! ∂t ℓ! ∂t i=0 ℓ=0 ℓ=0 p−1
q
=
X
=
q X kℓ X
βi
i=0
p−1
ℓ=0 n
ℓ!
i=0
ℓ
βi i −
q−1 X
γi iℓ
i=0
∂ℓ
∂tℓ
B(u)(tn ) + ϕn
=ϕ , where the last equality holds in view of the consistency properties of (α, β) and (α, γ), and, obviously, |ϕn | ≤ ck p . This relation and (3.5i,ii) yield max |E1n | ≤ C(k p + hr ).
(3.6)
0≤n≤N −q
Finally, using (1.6), we obtain max kE2n k⋆ ≤ Chr ,
(3.7)
0≤n≤N −q
which completes the estimate of E n . Convergence of the multistep scheme. Let ϑn = W n − U n , n = 0, . . . , N. Then (3.3) and (1.3) yield the error equation for ϑn (3.8)
q X
n+i
αi ϑ
+k
q X
n+i
βi Ah ϑ
=k
i=0
i=0
i=0
q−1 X
γi {Bh (W n+i ) − Bh (U n+i )} + kE n .
According to the splitting of E n , we introduce ϑj1 and ϑj2 , cf. section 2, by
(3.9)
q X i=0
αi ϑjn+i
+k
q X
βi Ah ϑjn+i
=k
q−1 X i=0
i=0
γi
Z
1 0
Bh′ (W n+i − sϑn+i ) ds ϑjn+i
+ kEjn , j = 1, 2, n = 0, . . . , N − q, with initial values ϑj1 = ϑj and ϑj2 = 0 for j = 0, . . . , q−1. Then, obviously, ϑn = ϑn1 +ϑn2 . In the sequel we shall use the notation n n+q−1 Ej ϑj 0 αi + βi x , Θjn = ... , Ejn = .. , δi (x) = − αq + βq x . n ϑj 0 ∆i = δi (kAh ),
Γni
= γi
Z
0
1
Bh′ (W n+i − sϑn+i ) ds,
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
∆q−1 ∆q−2 . . . ∆0 I 0 0 Λ = Λ(kAh ) = , . . . . . . 0 I 0
and
11
q−1 Γn . . . Γn0 0 ... 0 Γn = . .. , . . . 0 ... 0
(αq + kβq Ah )ϑn+q−1 j .. (αq + kβq Ah ) Θjn = . . n (αq + kβq Ah )ϑj Equation (3.9) can then be written in the form
(3.10)
(αq + kβq Ah )Θjn+1 = (αq + kβq Ah ) ΛΘjn + k Γn Θjn + k Ejn ,
j = 1, 2.
We quote the following result from Crouzeix, [3]: Lemma 3.1. There exist a constant η, with 0 ≤ η < 1, and a continuous map H : ¯ + → Cq×q such that for all x ≥ 0 the matrix H(x) is invertible and the Euclidean R norm k · k2 of the matrix L(x), L(x) =
αq + βq x H(x)−1 Λ(x)H(x), αq + ηβq x
is less or equal to one for all x ≥ 0. Let H = H(kAh ),
and
L = L(kAh ),
Yjn = H−1 Θjn , Γen = H−1 Γn , Eejn = H−1 Ejn ;
then, we can rewrite (3.10) as (3.11)
(αq + kβq Ah )Yjn+1 = (αq + kηβq Ah ) LYjn + k Γen HYjn + k Eejn ,
j = 1, 2.
In view of the fact that kH(x)k2 , kH(x)−1 k2 are uniformly bounded, see relations (3.27) and (3.28) in [3], it suffices to estimate Y n . We adjust in this section the definition of ||| · ||| to the scheme under consideration by setting |||v||| := (αq |v|2 + βq kkvk2 )1/2 ,
v ∈ V.
Further, for V = (v1 , . . . , vq )T and W = (w1 , . . . , wq )T in H q or in V q we shall use the notation q q X 1/2 X (vi , wi ), |V | := |vi |2 (V, W ) := , i=1
X q
kV k :=
i=1
kvi k
2
1/2
i=1
X q
,
|||V ||| :=
i=1
|||vi |||
2
1/2
,
kV k⋆ :=
q X i=1
and, for a linear operator M : H q → H q , we set |M| := supV ∈H q ,V 6=0 The main result in this paper is given in the following theorem:
kvi k2⋆
|M V | . |V |
1/2
,
12
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
Theorem 3.1. Assume that the constant λ in (1.4) is appropriately small (depending on the particular scheme) and that k and h2r k −1 are sufficiently small. Then, we have the local stability estimate |ϑn | + k 1/2 kϑn k ≤ CecmT (3.12)
q−1 nX j=0
n−q X |E1j | |ϑj | + k 1/2 kϑj k + k j=0
n−q
+ k
X
kE2j k2⋆
j=0
1/2 o
,
n = q − 1, . . . , N,
and the error estimate max |u(tn ) − U n | ≤ C(k p + hr ).
(3.13)
0≤n≤N
Proof. Let ρn = un − W n , n = 0, . . . , N. Then, according to ((1.5), max |ρn | ≤ Chr .
(3.14)
0≤n≤N
Now, if we assume that (3.12) holds, using (3.2), (3.6) and (3.7), we obtain max |ϑn | ≤ C(k p + hr ) ,
(3.15)
0≤n≤N
and (3.13) follows immediately from (3.14) and (3.15). Thus, it remains to prove (3.12). According to (3.2), (3.6) and (3.7)), there exists a constant C⋆ such that the right-hand side of (3.12) can be estimated by C⋆ (k p + hr ), cmT
(3.16)
Ce
q−1 nX j=0
j
|ϑ | + k
1/2
j
kϑ k + k
n−q X j=0
|E1j |
+ k
N −q
X j=0
kE2j k2⋆
1/2 o
≤ C⋆ (k p + hr ) . We will estimate ϑn by estimating Yjn . In fact, we shall show that for some ε = ε(λ, (α, β, γ) ), 0 < ε < (1 − η 2 )βq with η as in Lemma 3.1, (3.17)
n−1 n k X ej o |||Y1n ||| ≤ ecmtn |||Y10 ||| + √ |E | αq j=0 1
and n−1
(3.18)
|||Y2n |||
cmtn
≤e
X j 1/2 1 √ k kEe2 k2⋆ . ε j=0
Then, (3.12) follows, and the proof will be complete. We shall use induction: The estimates (3.17) and (3.18) are valid for n = 0. Assume that they hold for 0, . . . , n, 0 ≤ n ≤ N −q. Then, according to (3.16) and (3.12), which is then valid for 0, . . . , n+ q −1, we have, for k and k −1 h2r small enough, max
0≤j≤n+q−1
kϑj k ≤ C⋆ (k p−1/2 + k −1/2 hr ) ≤ 1,
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
13
i.e., (3.19)
sup
max
0<s 0 we seek a realvalued function u defined on R × [0, T ], 1-periodic in the space variable and satisfying (4.1)
ut + uux + uxx + νuxxxx = 0 in R × [0, T ]
and (4.2)
u(·, 0) = u0
in R,
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
15
where u0 is a given, smooth 1-periodic function. The periodic initial value problem (4.1)–(4.2) is well-posed, cf. [13], [17]. For numerical methods to this problem see, e.g., [1] and the references therein. The KS equation is related to turbulence phenomena in chemistry and combustion, and arises also in plasma physics and in two-phase flows in cylindrical geometries, see [18] and [14]. s For s ∈ N0 , let Hper denote the periodic Sobolev space of order s, consisting of the s . The 1-periodic elements of H sloc (R), and let k · kH s be the norm over a period in Hper 2 0 inner product in H := Lper = Hper is denoted by (·, ·), and the induced norm by | · |. 2 Let A : H 4per → H be defined by Av = ν(vxxxx + v). Then V := D(A1/2 ) = Hper , and 1/2 2 2 1/2 the norm in V is given by kvk = ν (|vxx | + |v| ) . Let B : V → H be given by B(v) = −vxx − vvx + νv. Then, B ′ (v)w = −wxx − (vw)x + νw , and thus by periodicity (B ′ (v)w, ω) = (wx , ωx ) + (v, wωx ) + ν(w, ω) . 1 , Therefore, in view of the inequality kwkL∞ ≤ |w| + |w ′|, w ∈ Hper
|(B ′ (v)w, ω)| ≤|wx | |ωx| + |v|kwkL∞ |ωx | + ν|w| |ω|
≤(1 + |v|)|wx | |ωx| + ν|w| |ω| + |v| |w| |ωx|,
and thus, using the inequality |ux |2 ≤ |u| |uxx|, u ∈ V, we easily see that the condition (1.4) is satisfied for any λ > 0, (4.3)
|(B ′ (v)w, ω)| ≤ λkwkkωk + µ(v)|w| |ω|,
∀v, w, ω ∈ V,
1 [ 1 + 2λν 2 + 2|v|(1 + |v|) ] . We note that, since λ can be taken arbiwith µ(v) := 2λν trarily small, our results hold for this problem for all implicit-explicit schemes satisfying our stability and consistency assumptions; further, since µ is bounded for v bounded in H, the meshcondition of Theorems 2.1 and 3.1 is not needed here, cf. Remark 2.1. For the space discretization, we let 0 = x0 < x1 < · · · < xJ = 1 be a partition of [0, 1], and h := maxj (xj+1 − xj ). Setting xjJ+s := j + xs , j ∈ Z, s = 0, . . . , J − 1, this partition is periodically extended to a partition of R. For integer r ≥ 4, let Vh denote a space of at least once continuously differentiable, 1-periodic splines of degree r − 1, in which approximations to the solution u(·, t) of (4.1)–(4.2) will be sought for 0 ≤ t ≤ T . The following approximation property of the family {Vh }0 0 we seek a real-valued function u, defined on ¯ × [0, T ], satisfying Ω ut − ∆u = f (u) in Ω × [0, T ], u=0
(4.18)
on ∂Ω × [0, T ],
u(·, 0) = u0
in Ω,
where u0 is a given smooth function and f : R → R is a smooth function. We are interested in approximating smooth solutions of this problem, and assume therefore that the data are smooth and compatible such that (4.18) gives rise to sufficiently regular solutions. We shall distinguish two cases. Assuming that f satisfies a polynomial growth condition of order at most ρ, see (4.19) below, and provided that ρ ≤ 4 for ν = 3, we show that the abstract theory of sections 2 and 3 is directly applicable and yields, without any additional assumptions on the discretization spaces, optimal order error estimates for all (α, β, γ) schemes considered in this paper. For general f our results apply as well, provided that meshconditions stronger than those of Theorem 3.1 are fulfilled. Reaction diffusion equations model various physical phenomena related, for instance, to phase transitions, chemical reactions, pattern formation, cf., e.g., [2], [7], [10] and their references. The Allen-Cahn model (ρ = 3), which in the limit describes evolution by mean curvature, [2], [7], and generalized Ginzburg-Landau equations, [10], are reaction diffusion equations of particular interest with polynomial nonlinearity. Let H s = H s (Ω) be the usual Sobolev space of order s, and k · kH s be the norm of s H . The inner product in H := L2 (Ω) is denoted by (·, ·), and the induced norm by | · |; the norm of Ls (Q), 1 ≤ s ≤ ∞, is denoted by k · kLs (Q) and simply by k · kLs for Q = Ω. Obviously, V = H01 (Ω) = H01 and the norm in V is equivalent to the H 1 norm. We now consider the case that f satisfies the following growth condition: there exists L ∈ R such that |f ′ (ξ)| ≤ L(1 + |ξ|ρ−1) ∀ξ ∈ R .
(4.19)
In the sequel, we shall use the Sobolev inequality (4.20)
kvkLs ≤ CkvkH 1 ,
1 ≤ s < ∞ for ν = 2,
and 1 ≤ s ≤ 6 for ν = 3 ,
as well as its consequence (4.21)
1−ν/4
kvkL4 ≤ CkvkL2
ν/4
kvkH 1 ,
ν = 2, 3.
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
19
Let B : V → V ′ , B(v) = f (v). First, we note that B is well defined. Indeed, for v, w ∈ V, Z 1 (f (v), w) = f ′ (sv)v ds, w + (f (0), w), 0
and, therefore, by (4.19) and H¨older’s inequality, i hZ 5/6 ρ6/5 (4.22) |(f (v), w)| ≤ C |v| dx kwkL6 + (1 + kvkL2 ) kwkL2 , Ω
and thus, in view of (4.20), we see that f (v) ∈ V ′ , provided ρ ≤ 5 for ν = 3. Further, by (4.19), for v, w, ω ∈ V, Z ′ |(f (v)w, ω)| ≤ C |v|ρ−1|w||ω|dx + CkwkL2 kωkL2 Ω Z 1/2 2(ρ−1) ≤C |v| dx kwkL4 kωkL4 + CkwkL2 kωkL2 , Ω
i.e., in view of (4.21), (4.23)
|w| |(f ′ (v)w, ω)| ≤ Ckvkρ−1 L2(ρ−1)
4−ν 4
|ω|
4−ν 4
ν
ν
kwk 4 kωk 4 + C|w| |ω| .
Thus, B is differentiable. Furthermore, (4.23) and Young’s inequality (ab ≤ qε aq + ′
ε−q /q q ′ 1 b , q q′
(4.24)
+
1 q′
= 1 ) yield
|(B ′ (v)w, ω)| ≤ C ν4 εkwk kωk 4(ρ−1) ν 4−ν − 4−ν 4−ν +C 1+ 4 ε |w| |ω| . kvkL2(ρ−1)
Therefore, by Sobolev’s inequality, we see that, for ρ ≤ 4 when ν = 3 (and for any ρ when ν = 2), the assumption (1.4) is satisfied with λ = C ν4 ε, µ(v) = C + C
4−ν 4
ν
ε− 4−ν kvkL2(ρ−1)
4(ρ−1) 4−ν
,
and µ(v) is bounded for v bounded in V. Again, λ can be taken arbitrarily small and thus our theory applies to all (α, β, γ) schemes satisfying our stability and consistency assumptions. We further analyze the case that f satisfies the growth condition (4.19), with ρ ≤ 4 for ν = 3, by introducing the discretization spaces Vh ; the case of general f shall be discussed at the end of this section. For simplicity, let Vh be the subspace of V defined on a finite element partition Th of Ω, and consisting of piecewise polynomial functions of degree at most r − 1, r ≥ 2. Let hK denote the diameter of an element K ∈ Th , and h = maxK∈Th hK . We define the elliptic projection operator Rh , Rh : V → Vh , by (∇Rh v, ∇χ) = (∇v, ∇χ) ∀χ ∈ Vh . We assume that (we do not attempt here to deal with problems that may arise in the case of a curved boundary ∂Ω concerning the requirement Vh ⊂ V, cf. Remark 2.2) (4.25)
|v − Rh v| + hkv − Rh vk ≤ Chr kvkH r ,
v ∈ H r ∩ H01 ;
20
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
then, in particular, the estimate (1.5)) will hold in this case with d = 2. Next, we will verify (1.6). We shall further assume that (4.26)
sup{kRh vkL∞ : 0 < h < 1} ≤ CkvkH r ,
For v ∈ H r ∩ H01 , we have
B(v) − B(Rh v) =
Z
v ∈ H r ∩ H01 .
1
f ′ (v − τ (v − Rh v))dτ (v − Rh v);
0
using here the Sobolev inequality kvkL∞ ≤ CkvkH r (r ≥ 2) and (4.26), we conclude, in view of (4.25), ˜ |B(v) − B(Rh v)| ≤ Chr M(v)kvk Hr ,
(4.27)
˜ with M(v) bounded for v ∈ H r ∩ H01 bounded in H r . Thus, (1.6) is satisfied. Now, let W (t) := Rh u(t), and assume that we are given approximations U 0 , . . . , U q−1 ∈ Vh to u0 , . . . , uq−1 such that q−1 X
(4.28)
j=0
|W j − U j | + k 1/2 kW j − U j k ≤ c(k p + hr ).
Then, we define U n ∈ Vh , n = q, . . . , N, recursively by the (α, β, γ) scheme
(4.29)
q X
αi (U n+i ,χ) + k
q X i=0
i=0
βi (∇U n+i , ∇χ)
q−1
=k
X i=0
γi (f (U n+i ), χ) ∀χ ∈ Vh ,
n = 0, . . . , N − q,
where (α, β) and (α, γ) are multistep schemes of order p, and (α, β) is strongly A(0)−stable. Then, Theorem 3.1 yields, for sufficiently small k and k −1 h2r , the error estimate (4.30)
max |un − U n | ≤ c(k p + hr ). n
Remark 4.2. By (4.24), we see that (1.4) holds for any λ > 0 with µ(v) bounded for v bounded in L2(ρ−1) . Using this fact, Sobolev’s inequality kvkLs ≤ C|v|1−akvka ,
a = ν s−2 , ν = 2, 3, (with s ≤ 6, of course, for ν = 3) and Remark 2.1, it is easily seen 2s that the meshcondition k −1 h2r ≤ c, under which (4.30) holds, can be weakened. We shall not dwell on this. We close this subsection by briefly considering the case of a general smooth function f, as well as the case that f satisfies (4.19) for ν = 3 but with ρ > 4. First, we note that in our analysis it suffices to assume that B is well defined and differentiable on a subspace Ve of V ∩ L∞ containing Vh , for all h. By tracing back the proof of (4.24) we see that in this case µ(v) is bounded, provided that Z |f ′ (v(x))|2 dx Ω
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
21
is bounded. Note that the assumption k −1 h2r ≤ c, for appropriate c, of Theorem 3.1 is only used to show that kϑn k ≤ 1, which implies (3.19), i.e., that µ(W j − sϑj ), s ∈ (0, 1), is bounded by a constant. In the case under investigation, by using appropriate inverse inequalities, we show that if stronger meshconditions are R ′ satisfied, then sup0<s 4, and ν = 3 and general f. i. ν = 2 and general f. First, we note that ∀χ ∈ Vh , kχkL∞ ≤ C| log(h)|1/2 kχkH 1 R with h = minK∈Th hK , cf. [19, p. 67]. Obviously, Ω |f ′ (χ(x))|2 dx is bounded if kχkL∞ is bounded. Now, max
0≤j≤n+q−1
kϑj kL∞ ≤ C| log(h)|1/2
max
0≤j≤n+q−1
kϑj k,
and thus, according to (3.16), max
0≤j≤n+q−1
kϑj kL∞ ≤ C⋆ C| log(h)|1/2 (k p−1/2 + k −1/2 hr ).
Therefore, if k and h are chosen such that | log(h)|k 2p−1 and | log(h)|k −1 h2r are sufficiently small, then µ(W j − sϑj ) will be bounded, and the convergence results hold. ii. ν = 3 and f satisfies (4.19) with ρ > 4. If s ≥ 6, we have kχkLs (Ω) ≤ Ch−
(4.31)
s−6 2s kχkH 1 (Ω)
∀χ ∈ Vh .
Indeed, employing standard homogeneity arguments, one can show that for an arbitrary element K ∈ Th , s−6 2s
−
kχkLs (K) ≤ ChK
k∇χkL2 (K) + kχkL6 (K) ,
and (4.31) follows in view of (4.20). Hence, max
0≤j≤n+q−1
ρ−4 − 2(ρ−1)
kϑj kL2(ρ−1) ≤ Ch −
Therefore, if k and h are such that h then, as before, max
0≤j≤n+q−1
ρ−4 ρ−1 k 2p−1
max
0≤j≤n+q−1 −
and k −1 h
kϑj k.
ρ−4 ρ−1 h2r
are sufficiently small,
ρ−4
kϑj kL2(ρ−1) ≤ C⋆ h− 2(ρ−1) (k p−1/2 + k −1/2 hr ) ≤ 1,
µ(W j − sϑj ) will be bounded in view of (4.24) and (4.26), and our convergence results hold. iii. ν = 3 and general f. In this case, kχkL∞ ≤ Ch−1/2 kχkH 1
∀χ ∈ Vh ,
22
GEORGIOS AKRIVIS, MICHEL CROUZEIX, AND CHARALAMBOS MAKRIDAKIS
as one can see by modifying the proof of the two dimensional case. By the same arguments as for ν = 2, the convergence results of this paper hold, provided that h−1 k 2p−1 and k −1 h−1 h2r are sufficiently small.
References 1. G. Akrivis, High-order finite element methods for the Kuramoto-Sivashinsky equation, RAIRO Mod´el. Math. Anal. Num´er. 30 (1996), 157–183. MR 97e:65095 2. S.M. Allen and J.W. Cahn, A macroscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall. 27 (1979), 1085-1095. 3. M. Crouzeix, Une m´ethode multipas implicite-explicite pour l’approximation des ´equations d’´evolution paraboliques, Numer. Math. 35 (1980), 257–276. MR 82b:65084 4. M. Crouzeix and P.-A. Raviart, Approximation des ´equations d’´evolution lin´eares par des m´ethodes a pas multiples, C.R. Acad. Sc. Paris, S´erie A 283 (1976), 367–370. MR 54:14377 ` 5. C. W. Cryer, A new class of highly stable methods, BIT 13 (1973), 153–159. MR 48:1469 6. C.M. Elliott and S.-M. Zheng, On the Cahn-Hilliard equation, Arch. Rational Mech. Anal. 96 (1986), 339–357. MR 87k:80007 7. L.C. Evans, H.M. Soner and P.E. Souganidis, Phase transitions and generalized motion by mean curvature, Comm. Pure Appl. Math. 45 (1992), 1097–1123. MR 93g:35064 ¨ 8. R.D. Grigorieff and J. Schroll, Uber A(α)-stabile Verfahren hoher Konsistenzordnung, Computing 20 (1978), 343–350. MR 83b:65086 9. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, Springer-Verlag, Berlin Heidelberg, Springer Series in Computational Mathematics v. 14, 1991. MR 92a:65016 10. Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer-Verlag, New York, Series in Synthetics, v. 19, 1984. MR 87e:92054 11. M.N. Le Roux, Semi-discr´etisation en temps pour les ´equations d’´evolution paraboliques lorsque l’op´erateur d´epend du temps, RAIRO 13 (1979), 119–137. MR 80e:65056 12. W.R. McKinney, Optimal error estimates for high order Runge-Kutta methods applied to evolutionary equations, Ph.D. thesis, University of Tennessee, Knoxville, 1989. 13. B. Nicolaenko and B. Scheurer, Remarks on the Kuramoto-Sivashinsky equation, Physica 12 D (1984), 391–395. MR 86d:80007 14. D.T. Papageorgiou, C. Maldarelli and D.S. Rumschitzki, Nonlinear interfacial stability of coneannular film flow, Phys. Fluids A2 (1990), 340–352. MR 91b:76067 15. G. Savar´e, A(Θ)-stable approximations of abstract Cauchy problems, Numer. Math. 65 (1993), 319–336. MR 94h:65062 16. L.L. Schumaker, Spline Functions: Basic Theory, Wiley, New York, 1981. MR 82j:41001 17. E. Tadmor, The well-posedness of the Kuramoto-Sivashinsky equation, SIAM J. Math. Anal. 17 (1986), 884–893. MR 87g:35117 18. R. Temam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics, Springer-Verlag, New York, 1988. MR 89m:58056 19. V. Thom´ee, Galerkin finite element methods for parabolic problems, Springer-Verlag, Lecture Notes in Mathematics v. 1054, 1984. MR 86k:65006 20. M. Zl´amal, Finite element multistep discretizations of parabolic boundary value problems, Math. Comp. 29 (1975), 350–359. MR 51:7326 21. M. Zl´amal, Finite element methods for nonlinear parabolic equations, RAIRO 11 (1977), 93–107. MR 58:19239
MULTISTEP METHODS FOR NONLINEAR PARABOLIC EQUATIONS
23
Department of Computer Science, University of Ioannina, 451 10 Ioannina, Greece E-mail address: [email protected] IRMAR, Universit´ e de Rennes I, Campus de Beaulieu, F-35042 Rennes, France E-mail address: [email protected] Department of Mathematics, University of Crete, 714 09 Heraklion, Crete, Greece, and IACM, Foundation for Research and Technology - Hellas, 711 10 Heraklion, Crete, Greece E-mail address: [email protected]
