BACKWARD EULER TYPE METHODS FOR PARABOLIC

INTEGRO-DIFFERENTIAL EQUATIONS IN BANACH SPACE

N. Yu. Bakaev, S. Larsson, and V. Thomée

Abstract. Time discretization by backward Euler type methods for a parabolic equation with memory is studied. Stability and error estimates are proved under conditions that permit quadrature rules for approximation of the memory term that have reduced storage requirements. The analysis takes place in a Banach space framework, and the results are used to derive error estimates in the L_{2} and maximum norms for piecewise linear finite element discretization in two space dimensions. Résumé. On étudie la discrétisation en temps d’une équation parabolique avec mémoire par des méthodes de type Euler rétrograde. On montre la stabilité et on donne des estimations d’erreur sous des hypothèses qui permettent d’utiliser des formules de quadrature peu exigeantes en stockage pour l’approximation du terme de mémoire. L’analyse est effectuée dans le cadre des espaces de Banach. Appliqués en dimension deux, ces résultats permettent d’obtenir des estimations d’erreur L_{2} et uniforme pour une discrétisation utilisant des éléments finis linéaires par morceaux.

1. Introduction. We consider the initial value problem

(1.1)  u_{t} + Au = ∫_{0}^{t} B(t, s)u(s) ds + f(t), for t ∈ [0, T], with u(0) = v,

in a Banach space X, where A is a closed linear operator with dense domain D(A), and B(t, s) is a smooth linear operator with D(B(t, s)) ⊃ D(A) and such that Q(t, s) = A^{−1}B(t, s) and Q_{t}(t, s) are uniformly bounded for 0 ≤ s ≤ t ≤ T.

We assume that −A generates a bounded analytic semigroup E(t) = e^{−tA}, so that

(1.2)  ∥E(t)∥ + t∥AE(t)∥ ≤ M, for t > 0.

It then follows, by Gronwall’s lemma (see Theorem 2.1 below), that for the solution of (1.1),

(1.3)  ∥u(t)∥ ≤ e^{C(T)M}M(∥v∥ + ∫_{0}^{t} ∥f∥ ds), for t ≤ T.

1991 Mathematics Subject Classification. 65J10. Key words and phrases. Integro-differential equation, parabolic, backward Euler, sparse quad-

rature, finite element method, Banach space, maximum norm.

Published in RAIRO Modél. Math. Anal. Numér. 32 (1998), 85–99. This work was supported by the Swedish Research Council for Engineering Sciences (TFR).

Typeset by AMS-TEX

We shall consider the time discretization of (1.1). Let k be a time step, set t_{n} = nk, n = 0, 1, 2, …, and denote by U^{n} the approximation of u(t_{n}). We shall replace the time derivative by the backward difference quotient ¯∂U^{n} = (U^{n} − U^{n−1})/k and approximate the memory term by a quadrature formula

(1.4)  σ^{n}(φ) = ∑_{j=0}^{n−1} ω_{nj}φ^{j} ≈ ∫_{0}^{t_{n}} φ(t) dt, where φ^{j} = φ(t_{j}).

The backward Euler discretization of (1.1) is then

(1.5)  \begin{aligned}¯∂U^{n} + AU^{n} & = ∑_{j=0}^{n−1} ω_{nj}B(t_{n}, t_{j})U^{j} + f(t_{n}) \\ & = σ^{n}(B_{n}U) + f^{n}, for n ≥ 1, with U^{0} = v. \end{aligned}

Our aim is to extend stability properties such as (1.3) to this discrete problem, and to use these to derive error estimates. In doing so we need to make the following more specific assumptions about the choice of the quadrature formula in (1.4), namely, for some positive integer p,

(1.6) |ϵ^{n}(φ)| ≤ Ck ∫_{0}^{t_{n}} ∑_{l=1}^{p} |φ^{(l)}(t)| dt, where ϵ^{n}(φ) = σ^{n}(φ) − ∫_{0}^{t_{n}} φ(t) dt,

and, for some positive number q,

(1.7)  ∑_{s=0}^{j−1} |ω_{js} − ω_{ns}| ≤ Ck^{q}, for 0 < t_{j} ≤ t_{n} ≤ T.

Under these hypotheses we shall show (Theorem 2.3), that, for small k,

(1.8)  ∥U^{n}∥ ≤ Ce^{C(T)M}M(∥v∥ + k ∑_{j=1}^{n} ∥f^{j}∥), for t_{n} ≤ T.

One example of a quadrature formula satisfying (1.6) (with p = 1) and (1.7) (trivially) is the left-side rectangle rule, corresponding to ω_{nj} = k for j < n. Since this rule requires the storage of all previous U^{j}, sparse quadrature rules have been proposed, e.g., in [4], [6]. A short discussion of such rules is given in Section 5 below, where it is shown that our present assumptions on σ^{n} are different from those made in the earlier work and do not require so called “dominated weights”.

Assume that we want to apply the above result to the case when X = C_{0}(Ω̄) equipped with the maximum norm, where Ω is a smooth domain in R^{2}, and where A is a discrete analogue of the Laplacian −Δ based on piecewise linear finite element spaces S_{h} defined by a family of quasi-uniform triangulations and where B(t, s) is appropriate. In this case it is known, cf. [3], that (1.2) holds with M = C log(1/h). This means that (1.3) and (1.8) contain the stability factor e^{C log(1/h)} = h^{−C}, which is unbounded as h tends to 0, and these stability estimates are therefore of little value.

In order to find a remedy for this we shall assume that, in addition to (1.2), we also have

(1.9)  ∥E(t)∥ + t∥AE(t)∥ ≤ M_{δ}t^{−δ}, for t > 0, for any δ ∈ (0, 1).

This is the case in the above finite element application with M_{δ} independent of h, see Lemma 4.1.

Under these assumptions we shall show that, for the solution of (1.1) (Theorem 2.2),

(1.10)  ∥u(t)∥ ≤ C(T, M_{δ}, δ)M(∥v∥ + ∫_{0}^{t} ∥f∥ ds), for t ≤ T,

and, for the solution of (1.5) (Theorem 2.4), if δ and k are sufficiently small,

(1.11)  ∥U^{n}∥ ≤ C(T, M_{δ}, δ)M(∥v∥ + k ∑_{j=1}^{n} ∥f^{j}∥), for t_{n} ≤ T.

In the above finite element application the bound now contains a single factor log(1/h).

We note that, in the case of one space dimension, (1.2) holds in maximum norm for finite elements of any order with a constant independent of h, cf. [1], so that (1.8) shows a uniform bound for T bounded.

The above stability estimates are proved in Section 2 below, and in Section 3 we give corresponding error estimates. In Section 4 we discuss the application to the finite element case in more detail, and Section 5 is concerned with sparse quadrature rules.

2. Stability estimates. We begin with the basic stability result (1.3) for the continuous equation.

Theorem 2.1. Assume that (1.2) holds. Then, for the solution of (1.1),

∥u(t)∥ ≤ e^{C(T)M}M(∥v∥ + ∫_{0}^{t} ∥f∥ ds), for t ≤ T.

Proof. Using Duhamel’s principle, we have

(2.1) \begin{aligned}u(t) & = E(t)v + ∫_{0}^{t} E(t − s)f(s) ds + ∫_{0}^{t} E(t − y) ∫_{0}^{y} B(y, s)u(s) ds dy \\ & = F (t) + ∫_{0}^{t} G(t, s)u(s) ds, \end{aligned}

where F (t) = E(t)v + ∫_{0}^{t} E(t − s)f(s) ds, and, since AE(t) = −E′(t),

\begin{aligned}G(t, s) & = ∫_{s}^{t} AE(t − y)Q(y, s) dy \\ & = (I − E(t − s))Q(t, s) + ∫_{s}^{t} AE(t − y)(Q(y, s) − Q(t, s)) dy, \end{aligned}

where Q(t, s) = A^{−1}B(t, s). Since by our assumption

(2.2)  ∥Q(t, s)∥ + ∥Q_{t}(t, s)∥ ≤ C, 0 ≤ s ≤ t ≤ T,

we may use (1.2) to conclude that

(2.3) ∥G(t, s)∥ ≤ C(1 + ∥E(t − s)∥) + C ∫_{s}^{t}(t − y)∥AE(t − y)∥ dy ≤ C(T)M,

and we hence obtain from (2.1)

∥u(t)∥ ≤ M(∥v∥ + ∫_{0}^{t} ∥f∥ ds) + C(T)M ∫_{0}^{t} ∥u(s)∥ ds, for t ≤ T,

from which the desired result follows by Gronwall’s lemma. □

We now turn to the modified stability estimate (1.10) for the continuous equation.

Theorem 2.2. Assume that (1.2) and (1.9) hold. Then for the solution of (1.1) we have, for any δ ∈ (0, 1),

∥u(t)∥ ≤ C(T, M_{δ}, δ)M(∥v∥ + ∫_{0}^{t} ∥f∥ ds), for t ≤ T.

Proof. This time we use (1.9) in the first inequality of (2.3) to obtain

∥G(t, s)∥ ≤ C + CM_{δ}(t − s)^{−δ} + CM_{δ} ∫_{s}^{t}(t − y)^{−δ} dy ≤ C(T, M_{δ})(t − s)^{−δ}.

Hence, by (2.1),

(2.4)  ∥u(t)∥ ≤ M(∥v∥ + ∫_{0}^{t} ∥f∥ ds) + C(T, M_{δ}) ∫_{0}^{t}(t − s)^{−δ}∥u(s)∥ ds,

and our result follows by a variant of Gronwall’s lemma (cf. [2, Lemma 5.6.7]; it also follows easily from the time-continuous version of Lemma 2.2 below). □

We note that under assumption (1.9) one may also show for the solution of (1.1)

(2.5)  ∥u(t)∥ ≤ C(T, M_{δ}, δ)(t^{−δ}∥v∥ + ∫_{0}^{t}(t − s)^{−δ}∥f(s)∥ ds), for t ≤ T.

In fact, instead of (2.4) one has, with C = C(M_{δ}),

∥u(t)∥ ≤ C(t^{−δ}∥v∥ + ∫_{0}^{t}(t − s)^{−δ}∥f∥ ds) + C ∫_{0}^{t}(t − s)^{−δ}∥u(s)∥ ds,

from which (2.5) follows by a variant of Gronwall’s lemma.

In order to prove the discrete analogues (1.8) and (1.11) we introduce the backward Euler one step evolution operator E_{k} = (I+kA)^{−1} and first show the following discrete analogues of (1.2) and (1.9) (with the former contained as a special case with M_{0} = M).

Lemma 2.1. Assume that, for some δ ∈ [0, 1),

∥E(t)∥ + t∥AE(t)∥ ≤ M_{δ}t^{−δ}, for t > 0.

Then

∥E_{k}^{n}∥ + t_{n}∥AE_{k}^{n}∥ ≤ CM_{δ}t_{n}^{−δ}, for t_{n} > 0.

Proof. We have, cf., e.g., [2, p. 21],

E_{k}^{n} = (I + kA)^{−n} = \frac{1}{(n − 1)!} ∫_{0}^{∞} t^{n−1}e^{−t}E(kt) dt, for n ≥ 1.

Hence, by our assumption on E(t), and since Γ(n − δ)/Γ(n) ≤ Cn^{−δ},

∥E_{k}^{n}∥ ≤ M_{δ} \frac{k^{−δ}}{(n − 1)!} ∫_{0}^{∞} t^{n−1−δ}e^{−t} dt = M_{δ}k^{−δ} \frac{Γ(n − δ)}{Γ(n)} ≤ CM_{δ}t_{n}^{−δ}.

The estimation of AE_{k}^{n} is similar for n ≥ 2; for n = 1 we use kAE_{k} = I − E_{k}. □

Our proofs in the discrete case will depend on the following version of Gronwall’s lemma. We remark that, in the earlier work [4], [6], instead of (1.7) the analysis was based on the assumption that the weights ω_{ns} are “dominated” in the sense that there are weights ω_{s}, independent of n and such that ω_{ns} ≤ ω_{s} for 0 ≤ t_{s} < t_{n} ≤ T with ∑_{s=0}^{n−1} ω_{s} ≤ C, which makes it possible to apply a more standard Gronwall lemma.

Lemma 2.2. Let u_{n}, n = 0, 1, …, be nonnegative numbers such that

(2.6)  u_{n} ≤ K + ∑_{s=0}^{n−1} μ_{ns}u_{s}, for t_{n} = nk ≥ 0.

Assume that the coefficients μ_{ns} = μ_{ns}(k) are nonnegative and that there exist positive numbers γ, τ and k_{0}, independent of k, such that

(2.7)  ∑_{s=m}^{l−1} μ_{ns} ≤ γ < 1, for t_{l−1} − t_{m} ≤ τ, 0 ≤ t_{m} < t_{l} ≤ t_{n}, k ≤ k_{0}.

Then, with C = log ((1 − γ)^{−1}),

(2.8)  u_{n} ≤ Ke^{C(t_{n}/τ+1)}, for t_{n} ≥ 0, k ≤ k_{0}.

Proof. Let I_{j} = [(j − 1)τ, jτ]. Then, by (2.6) and (2.7),

max_{t_{n}∈I_{j}} u_{n} ≤ K + γ ∑_{l=1}^{j} max_{t_{n}∈I_{l}} u_{n},

from which we easily conclude that max_{t_{n}∈I_{j}} u_{n} ≤ K(1 − γ)^{−j}, which implies (2.8) since j ≤ t_{n}/τ + 1 for t_{n} ∈ I_{j}. □

In the following we let T be fixed (but arbitrary) and we let C denote various constants that may depend on T.

We now show the following stability result for (1.5). For the purpose of our error estimates below it is phrased in a more general way than (1.8), which latter is contained for g^{j} = 0.

Theorem 2.3. Assume (1.2) and let {g^{j}} ⊂ X be arbitrary. Then we have, for the solution of (1.5), for t_{n} ≤ T, k ≤ k_{0}(M),

∥U^{n}∥ ≤ Ce^{CM}M(∥v∥ + k ∑_{j=1}^{n} ∥f^{j} − g^{j}∥ + ∥A^{−1}g^{1}∥ + k ∑_{j=2}^{n} ∥A^{−1} ¯∂g^{j}∥).

Proof. The proof is modeled on the proof of Theorem 2.1. We have by Duhamel’s principle and a change of the order of summation in the double sum, cf. (2.1),

(2.9) U^{n} = E_{k}^{n}v+k ∑_{j=1}^{n} E_{k}^{n−j+1}f^{j}+k ∑_{j=1}^{n} E_{k}^{n−j+1} ∑_{s=0}^{j−1} ω_{js}B_{js}U^{s} = F^{n}+∑_{s=0}^{n−1} G_{ns}U^{s},

where, since kAE_{k}^{n} = E_{k}^{n−1} − E_{k}^{n}, we may write

F^{n} = E_{k}^{n}v + k ∑_{j=1}^{n} E_{k}^{n−j+1}(f^{j} − g^{j}) + A^{−1}g^{n} − E_{k}^{n}A^{−1}g^{1} − k ∑_{j=2}^{n} E_{k}^{n−j+1}A^{−1} ¯∂g^{j},

and, with Q_{js} = A^{−1}B_{js}, we have G_{ns} = k ∑_{j=s+1}^{n} AE_{k}^{n−j+1}ω_{js}Q_{js}. Here, by (1.2) and Lemma 2.1, for small k and t_{n} ≤ T,

(2.10) \begin{aligned}∥F^{n}∥ & ≤ CM(∥v∥ + k ∑_{j=1}^{n} ∥f^{j} − g^{j}∥ + max_{j=1,n} ∥A^{−1}g^{j}∥ + k ∑_{j=2}^{n} ∥A^{−1} ¯∂g^{j}∥ \\ & ≤ CM(∥v∥ + k ∑_{j=1}^{n} ∥f^{j} − g^{j}∥ + ∥A^{−1}g^{1}∥ + k ∑_{j=2}^{n} ∥A^{−1} ¯∂g^{j}∥). \end{aligned}

In view of Lemma 2.2, the result therefore follows once we have shown that, for small k,

(2.11)  ∑_{s=m}^{l−1} ∥G_{ns}∥ ≤ CM(t_{l−m−1} + k^{q} log \frac{1}{k} ), for 0 ≤ t_{m} < t_{l} ≤ t_{n} ≤ T,

so that (2.7) holds for μ_{ns} = ∥G_{ns}∥, with γ = \frac{1}{2}, τ = (4CM)^{−1} and k_{0} = k_{0}(M). In order to prove (2.11) we write

\begin{aligned}G_{ns} & = k ∑_{j=s+1}^{n} AE_{k}^{n−j+1}ω_{ns}Q_{ns} + k ∑_{j=s+1}^{n} AE_{k}^{n−j+1}ω_{ns}(Q_{js} − Q_{ns}) \\ & + k ∑_{j=s+1}^{n} AE_{k}^{n−j+1}(ω_{js} − ω_{ns})Q_{js} ≡ G_{ns}^{1} + G_{ns}^{2} + G_{ns}^{3}. \end{aligned}

Here G_{ns}^{1} = (I − E_{k}^{n−s})ω_{ns}Q_{ns}, so that by (1.2) and Lemma 2.1 (cf. (2.3))

∥G_{ns}^{1} + G_{ns}^{2}∥ ≤ Cω_{ns}(1 + k ∑_{j=s+1}^{n} t_{n−j}∥AE_{k}^{n−j+1}∥) ≤ CMω_{ns}.

Also

∥G_{ns}^{3}∥ ≤ CMk ∑_{j=s+1}^{n} t_{n−j+1}^{−1}|ω_{js} − ω_{ns}|.

We shall show

(2.12)  ∑_{s=m}^{l−1} ω_{ns} ≤ t_{l−m} + Ck^{q},

and, for small k and t_{n} ≤ T,

(2.13)  k ∑_{s=m}^{l−1} ∑_{j=s+1}^{n} t_{n−j+1}^{−1}|ω_{js} − ω_{ns}| ≤ Ck^{q} log \frac{1}{k} ,

which together show (2.11). In order to prove (2.12) we note that the quadrature formula is exact for constants (cf. (1.6)) and use (1.7) to get

∑_{s=m}^{l−1} ω_{ns} = ∑_{s=0}^{l−1} ω_{ls} − ∑_{s=0}^{m−1} ω_{ms} − ∑_{s=0}^{l−1}(ω_{ls} − ω_{ns}) + ∑_{s=0}^{m−1}(ω_{ms} − ω_{ns}) ≤ t_{l−m} + Ck^{q}.

For the proof of (2.13) we use (1.7) as follows:

k ∑_{s=m}^{l−1} ∑_{j=s+1}^{n} t_{n−j+1}^{−1}|ω_{js} − ω_{ns}| = k ∑_{j=m+1}^{n} t_{n−j+1}^{−1} ∑_{s=m}^{min (l−1,j−1)} |ω_{js} − ω_{ns}|

≤ Ck^{q+1} ∑_{j=m+1}^{n} t_{n−j+1}^{−1} ≤ Ck^{q}(1 + log \frac{t_{n} − t_{m}}{k} ) ≤ Ck^{q} log \frac{t_{n}}{k} . □

The following is a modified stability result for the time stepping method.

Theorem 2.4. Assume (1.2) and (1.9), let {g^{j}} ⊂ X and r ∈ (1, ∞] be arbitrary. Then we have, for the solution of (1.5), for some δ ∈ (0, 1), for t_{n} ≤ T, k ≤ k_{0}(M_{δ}),

∥U^{n}∥ ≤ C(T, M_{δ}, δ){M(∥v∥ + k ∑_{j=1}^{n} ∥f^{j} − g^{j}∥) + M_{δ}(k ∑_{j=1}^{n} ∥g^{j}∥^{r})^{1/r}}.

Proof. Again we have (2.9) and instead of (2.10) we use (1.9) and Lemma 2.1 to get, for δ sufficiently small,

∥k ∑_{j=1}^{n} E_{k}^{n−j+1}g^{j}∥ ≤ CM_{δ}k ∑_{j=1}^{n} t_{n−j+1}^{−δ}∥g^{j}∥ ≤ C(T, δ)M_{δ}(k ∑_{j=1}^{n} ∥g^{j}∥^{r})^{1/r}.

It remains to bound ∑_{s=m}^{l−1} ∥G_{ns}∥, so that Lemma 2.2 may be applied with μ_{ns} = ∥G_{ns}∥ and with τ depending on M_{δ} instead of M. For this purpose we use (1.9) and Lemma 2.1 to get, similarly to the proof of Theorem 2.3, (this is a correction of an error in the version the paper printed in RAIRO M^{2}AN)

∥G_{ns}^{1} + G_{ns}^{2}∥ ≤ Cω_{ns}(1 + ∥E_{k}^{n−s}∥ + k ∑_{j=s+1}^{n} t_{n−j+1}∥AE_{k}^{n−j+1}∥) ≤ CM_{δ}t_{n−s}^{−δ}ω_{ns},

and

∥G_{ns}^{3}∥ ≤ CM_{δ}k ∑_{j=s+1}^{n} t_{n−j+1}^{−1−δ}|ω_{js} − ω_{ns}|.

Estimating ∑_{s=m}^{l−1} ∥G_{ns}^{3}∥ we have, uniformly in δ, (cf. the proof of (2.13))

\begin{aligned}k ∑_{s=m}^{l−1} ∑_{j=s+1}^{n} t_{n−j+1}^{−1−δ}|ω_{js} − ω_{ns}| & ≤ k ∑_{j=m+1}^{n} t_{n−j+1}^{−1−δ} ∑_{s=0}^{j−1} |ω_{js} − ω_{ns}| \\ & ≤ Ck^{q−δ+1} ∑_{j=m+1}^{n} t_{n−j+1}^{−1} ≤ Ck^{q−δ} log \frac{1}{k} . \end{aligned}

In order to bound ∑_{s=m}^{l−1} ∥G_{ns}^{1} + G_{ns}^{2}∥ in terms of t_{l−m} and k, we argue as in the proof of (2.12). We have

∑_{s=m}^{l−1} t_{n−s}^{−δ}ω_{ns} = ∑_{s=0}^{l−1} t_{n−s}^{−δ}ω_{ls} − ∑_{s=0}^{m−1} t_{n−s}^{−δ}ω_{ms}

− ∑_{s=0}^{l−1} t_{n−s}^{−δ}(ω_{ls} − ω_{ns}) + ∑_{s=0}^{m−1} t_{n−s}^{−δ}(ω_{ms} − ω_{ns})

= ∫_{t_{m}}^{t_{l}}(t_{n} − s)^{−δ} ds + ϵ^{l} − ϵ^{m} − η^{l} + η^{m},

where, uniformly for small δ, and using the elementary inequality (x+y)^{γ} ≤ x^{γ} +y^{γ} for x, y ≥ 0, γ ∈ (0, 1),

∫_{t_{m}}^{t_{l}}(t_{n} − s)^{−δ} ds = \frac{(t_{n} − t_{m})^{1−δ} − (t_{n} − t_{l})^{1−δ}}{1 − δ} ≤ \frac{t_{l−m}^{1−δ}}{1 − δ} ≤ Ct_{l−m}^{1−δ}.

Moreover, according to (1.6) we have

|ϵ^{i}| = |ϵ^{i}((t_{n} − ·)^{−δ})| ≤ Ck ∫_{0}^{t_{i}} ∑_{l=0}^{p} |D_{t}^{l}(t_{n} − t)^{−δ}| dt ≤ Ckt_{n−i}^{−p},

again uniformly for small δ, where for simplicity we have replaced δ by 1 in the final step. Also, according to (1.7),

|η^{i}| = | ∑_{s=0}^{i−1} t_{n−s}^{−δ}(ω_{is} − ω_{ns})| ≤ k^{−δ} ∑_{s=0}^{i−1} |ω_{is} − ω_{ns}| ≤ Ck^{q−δ}.

Together these estimates show

(2.14)  ∑_{s=m}^{l−1} t_{n−s}^{−δ}ω_{ns} ≤ Ct_{l−m}^{1−δ} + Ckt_{n−l}^{−p} + Ck^{q−δ}.

Since our estimate of ϵ^{l} may be as large as Ck^{1−p} (when t_{l} = t_{n−1}), we have to make a refined estimation. Let ϵ > 0 be arbitrary. The contribution of the terms in ∑_{s=m}^{l−1} t_{n−s}^{−δ}ω_{ns} with t_{n−s} ≥ (k/ϵ)^{1/p} is then bounded by Ct_{l−m}^{1−δ} + Cϵ + Ck^{q−δ} according to (2.14). On the other hand, the contribution of the terms with t_{n−s} < (k/ϵ)^{1/p} can be bounded, using (2.12), by

k^{−δ} ∑_{t_{n−s}<(k/ϵ)^{1/p}} ω_{ns} ≤ k^{−δ}((k/ϵ)^{1/p} + Ck^{q}).

Thus, putting these bounds together we have (cf. (2.11))

∑_{s=m}^{l−1} ∥G_{ns}∥ ≤ CM_{δ}(t_{l−m}^{1−δ} + ϵ + k^{−δ+1/p}ϵ^{−1/p} + k^{q−δ} log \frac{1}{k} ),

so that (2.7) holds for μ_{ns} = ∥G_{ns}∥, if 0 < δ < min(q, 1/p), by choosing ϵ = k^{(1−pδ)/(1+p)}, with k ≤ k_{0}(M_{δ}) sufficiently small, and τ = (4CM_{δ})^{−1/(1−δ)}. □

In analogy with (2.5) one may also show that under assumption (1.9) we have for the solution of (1.5)

∥U^{n}∥ ≤ C(T, M_{δ}, δ)(t_{n}^{−δ}∥v∥ + k ∑_{j=1}^{n} t_{n−j+1}^{−δ}∥f^{j}∥), for t_{n} ≤ T.

Since this will not be used below, we refrain from the details.

3. Error estimates. We show the following error estimate in our abstract framework.

Theorem 3.1. Assume that (1.2) holds. Then we have, for the solutions of (1.5) and (1.1), with C = C(T), for t_{n} ≤ T

∥U^{n} − u(t_{n})∥ ≤ Ce^{CM}Mk ∫_{0}^{t_{n}} (∥u_{tt}∥ + ∑_{l=1}^{p} ∥Au^{(l)}∥) dt.

If also (1.9) holds, then the result holds with the stability factor C(T, M_{δ})M, for some (sufficiently small) δ > 0.

Proof. Set e^{n} = U^{n} − u^{n}. Then

¯∂e^{n} + Ae^{n} = σ(B_{n}e) + τ_{1}^{n} + τ_{2}^{n},

where τ_{1}^{n} = −(¯∂u^{n} − u_{t}(t_{n})) and τ_{2}^{n} = ϵ^{n}(B_{n}u). Here, ∥τ_{1}^{n}∥ ≤ ∫_{t_{n−1}}^{t_{n}} ∥u_{tt}∥ dt, and using (1.6),

∥τ_{2}^{n}∥ ≤ Ck ∫_{0}^{t_{n}} ∑_{l=1}^{p} ∥Au^{(l)}∥ dt.

(Note that, by duality, the assumption (1.6) for scalar functions implies the corresponding statement for vector valued functions.) By the stability result of Theorem 2.3, we now have

∥e^{n}∥ ≤ Ce^{CM}Mk ∑_{j=1}^{n} (∥τ_{1}^{j}∥ + ∥τ_{2}^{j}∥),

from which the first result follows. The proof of the second is analogous, using Theorem 2.4 instead of Theorem 2.3. □

4. Application to piecewise linear finite elements. In this section we consider the case of the initial value problem (1.1) when A is a self-adjoint positive definite elliptic operator, and B(t, s) is a second order partial differential operator with smooth coefficients, in a plane convex domain Ω and with Dirichlet boundary conditions. Together with this problem we shall consider its spatial discretization in piecewise linear finite element spaces S_{h}. With (·, ·) the standard L_{2}-inner product on Ω, the spatially discrete analogue of (1.1) is to find u_{h}(t) ∈ S_{h} for t ≥ 0 such that

(4.1) \begin{aligned}(u_{h,t}, χ) + A(u_{h}, χ) & = ∫_{0}^{t} B(t, s; u_{h}(s), χ) ds + (f, χ), ∀χ ∈ S_{h}, t ∈ [0, T], \\ u_{h}(0) & = v_{h}, \end{aligned}

where A(·, ·) and B(t, s; ·, ·) are the standard bilinear forms associated with A and B(t, s). Introducing the discrete operators A_{h} and B_{h}(t, s) : S_{h} → S_{h} by

(A_{h}ψ, χ) = A(ψ, χ), (B_{h}(t, s)ψ, χ) = B(t, s; ψ, χ), ∀ψ, χ ∈ S_{h},

the problem (4.1) may be expressed as (1.1), with A and B(t, s) replaced by their discrete analogues A_{h} and B_{h}(t, s) in S_{h}. It is to this spatially discrete problem that we now apply the backward Euler discretization (1.5), which yields the completely discrete problem to find U^{n} ∈ S_{h} such that, with B_{n}(s; ·, ·) = B(t_{n}, s; ·, ·),

(4.2)  \begin{aligned}(¯∂U^{n}, χ) + A(U^{n}, χ) & = σ^{n}(B_{n}(U, χ)) + (f^{n}, χ), ∀χ ∈ S_{h}, n ≥ 1, \\ U^{0} & = v_{h}. \end{aligned}

We shall begin by considering this problem in the Hilbert space L_{2}(Ω). In order to apply the above theory to this problem, we recall the well known fact that for E_{h}(t) = e^{−A_{h}t}, we have, with respect to the L_{2}-norm,

∥E_{h}(t)∥ + t∥A_{h}E_{h}(t)∥ ≤ C, for t > 0,

so that (1.2) is valid. We also need to assume that, uniformly in h,

(4.3)  ∥A_{h}^{−1}B_{h}(t, s)∥ + ∥A_{h}^{−1}B_{h,t}(t, s)∥ + ∥B_{h}(t, s)A_{h}^{−1}∥ ≤ C

(cf. (2.2)). This is the case, e.g., if the triangulation underlying the definition of S_{h} is quasi-uniform, or if the principal part of B(t, s) equals a scalar function b(t, s) times the principal part of A, see [4], [5], [6].

Under these assumptions we have the following. We assume for simplicity that the discrete initial value is R_{h}v, where R_{h} : H_{0}^{1}(Ω) → S_{h} denotes the Ritz projection, i.e., the orthogonal projection with respect to the inner product A(·, ·).

BACKWARD EULER TYPE METHODS 11

Theorem 4.1. We have, for the solutions of (1.1) and (4.2) with v_{h} = R_{h}v,

∥U^{n} − u(t_{n})∥ ≤ C(T, u)(h^{2} + k), for t_{n} ≤ T.

Proof. We write

e^{n} = U^{n} − u^{n} = (U^{n} − R_{h}u^{n}) + (R_{h}u^{n} − u^{n}) = θ^{n} + ρ^{n}.

It is well known that

(4.4)  ∥ρ^{n}∥ ≤ Ch^{2}∥u^{n}∥_{H^{2}}.

For θ^{n} ∈ S_{h} we note that, with B_{h,n}(s) = B_{h}(t_{n}, s),

(4.5) ¯∂θ^{n} + A_{h}θ^{n} = σ^{n}(B_{h,n}θ) + τ^{n},

where, with ϵ^{n} defined in (1.6),

\begin{aligned}τ^{n} & = −¯∂(R_{h}u^{n} − u^{n}) + (¯∂u^{n} − u_{t}^{n}) + ϵ^{n}(B_{h,n}R_{h}u) + ∫_{0}^{t_{n}}(B_{h,n}R_{h} − P_{h}B_{n})u dt \\ & = τ_{1}^{n} + τ_{2}^{n} + τ_{3}^{n} + τ_{4}^{n}. \end{aligned}

Here, by standard estimates,

∥τ_{1}^{n}∥ + ∥τ_{2}^{n}∥ ≤ Ch^{2} \frac{1}{k} ∫_{t_{n−1}}^{t_{n}} ∥u_{t}^{n}∥_{H^{2}} dt + ∫_{t_{n−1}}^{t_{n}} ∥u_{tt}∥ dt ≤ C(u)(h^{2} + k).

Further, since ∥B_{h,n}A_{h}^{−1}∥ ≤ C and R_{h} = A_{h}^{−1}P_{h}A,

∥τ_{3}^{n}∥ = ∥ϵ^{n}(B_{h,n}R_{h}u)∥ ≤ Ck ∫_{0}^{t_{n}} ∑_{l=1}^{p} ∥Au^{(l)}∥ dt ≤ C(u)k.

For τ_{4}^{n} we have, with Z_{h,n} = B_{h,n}R_{h} − P_{h}B_{n}, k ¯∂τ_{4}^{n} = ∫_{0}^{t_{n}} Z_{h,n}u dt − ∫_{0}^{t_{n−1}} Z_{h,n−1}u dt = ∫_{t_{n−1}}^{t_{n}} Z_{h,n}u dt + k ∫_{0}^{t_{n−1}} ¯∂Z_{h,n}u dt. We shall show presently that

(4.6)  ∥A_{h}^{−1}Z_{h,n}u∥ + ∥A_{h}^{−1} ¯∂Z_{h,n}u∥ ≤ C(u)h^{2}.

We now apply Theorem 2.3 to (4.5), with f^{n} = ∑_{i=1}^{4} τ_{i}^{n}, g^{n} = τ_{4}^{n}, to obtain, since θ^{0} = 0,

\begin{aligned}∥θ^{n}∥ & ≤ C(T)(k ∑_{j=1}^{n}(∥τ_{1}^{j}∥ + ∥τ_{2}^{j}∥ + ∥τ_{3}^{j}∥) + ∥A_{h}^{−1}τ_{4}^{1}∥ + k ∑_{j=2}^{n} ∥A_{h}^{−1} ¯∂τ_{4}^{j}∥) \\ & ≤ C(T, u)(h^{2} + k). \end{aligned}

Together with (4.4) this completes the proof.

It remains thus to show (4.6). We have, taking the supremum over all χ ∈ S_{h} with ∥χ∥ = 1,

∥A_{h}^{−1}Z_{h,n}u∥ = sup_{χ}(A_{h}^{−1}Z_{h,n}u, χ) = sup_{χ}(A_{h}^{−1}(B_{h,n}R_{h} − P_{h}B_{n})u, χ)

= sup_{χ} B_{n}(ρ, A_{h}^{−1}χ) ≤ sup_{χ} B_{n}(ρ, (A_{h}^{−1} − A^{−1})χ) + sup_{χ} B_{n}(ρ, A^{−1}χ)

≤ C∥ρ∥_{H^{1}}h + C∥ρ∥ sup_{χ∈S_{h}} ∥B_{n}^{∗}A^{−1}χ∥ ≤ C(u)h^{2},

which completes the proof of the first part. The second part follows similarly from

∥A_{h}^{−1} ¯∂Z_{h,n}u∥ ≤ Ck^{−1} sup_{χ} (B_{n}(ρ, A_{h}^{−1}χ) − B_{n−1}(ρ, A_{h}^{−1}χ)) ≤ C(u)h^{2}. □

We now turn to a discussion of the above problem in the Banach space C_{0}(Ω̄) and throughout the rest of this section we now use the maximum norm ∥v∥ = sup_{x∈Ω} |v(x)|. We note that in the general case A^{−1}B(t, s) is then not a bounded operator and (4.3) cannot be expected to hold. We therefore now restrict the considerations to the case that A = −Δ, B = −b(t, s)Δ, where Δ is the Laplacian and b(t, s) is a smooth scalar function. In this case A(·, ·) = (∇·, ∇·) and B(t, s; ·, ·) = b(t, s)(∇·, ∇·), and the discrete analogues of A and B(t, s) are defined by A_{h} = −Δ_{h}, B_{h}(t, s) = −b(t, s)Δ_{h}, where

−(Δ_{h}ψ, χ) = (∇ψ, ∇χ), ∀ψ, χ ∈ S_{h}.

In order to apply our abstract theory in this case, we note that A_{h}^{−1}B_{h}(t, s) = b(t, s)I is bounded together with its derivatives. We also need to know to what extent the assumptions (1.2) and (1.9) are satisfied for A_{h} = −Δ_{h}. We first recall from [3] that, if the family of triangulations underlying the definition of S_{h} is quasiuniform, then for E_{h}(t) = e^{−A_{h}t} we have, with respect to the maximum norm,

(4.7)  ∥E_{h}(t)∥ + t∥A_{h}E_{h}(t)∥ ≤ C log \frac{1}{h}, for t > 0,

so that (1.2) is satisfied with M = C log(1/h). We also want to show that (1.9) is satisfied with M_{δ} independent of h.

Lemma 4.1. Under the present assumptions we have, for any δ ∈ (0, 1) and for h ≤ h_{δ},

∥E_{h}(t)∥ + t∥A_{h}E_{h}(t)∥ ≤ C_{δ}t^{−δ}, for t > 0.

Proof. We use techniques from [3, Theorem 3.3]. It is easy to show that ∥E_{h}(t)∥ and ∥A_{h}E_{h}(t)∥ are bounded for t ≥ t_{0} > 0 and decay exponentially (uniformly in h) as t → ∞, so it suffices to consider 0 < t ≤ 1, say. By the maximum principle, ∥E(t)∥ ≤ 1. We shall show that, with P_{h} the L_{2}-projection onto S_{h},

(4.8)  ∥E_{h}(t)P_{h}v − E(t)v∥ ≤ C_{δ}t^{−δ}∥v∥, for 0 < t ≤ 1, h ≤ h_{δ},

which, applied with v ∈ S_{h} shows the desired estimate for E_{h}(t).

From [3] we quote that, for any ϵ > 0,

(4.9)  ∥E_{h}(t)P_{h} − E(t)∥ ≤ C_{ϵ}h^{2−3ϵ}t^{−1+ϵ},

which implies (4.8) for any δ ≥ \frac{1}{3}. To consider smaller δ, we use (4.9) with ϵ = \frac{1}{2}, say, together with (4.7), to obtain, for 0 < δ ≤ \frac{1}{2} and h ≤ h_{δ},

∥E_{h}(t)P_{h} − E(t)∥ ≤ (C log \frac{1}{h})^{1−2δ}(Ch^{1/2}t^{−1/2})^{2δ} ≤ C_{δ}t^{−δ}.

We now turn to the estimate for A_{h}E_{h}(t) = −E′_{h}(t). As before the statement is valid for the continuous analogue E′(t), so it suffices to show

t∥E′_{h}(t)P_{h}v − E′(t)v∥ ≤ C_{δ}t^{−δ}∥v∥.

With u_{h}(t) = E_{h}(t)P_{h}v and u(t) = E(t)v we write

u_{h}(t) − u(t) = (u_{h}(t) − P_{h}u(t)) + (P_{h}u(t) − u(t)) = η + ζ.

Here, since P_{h} is bounded in maximum norm,

t∥ζ_{t}(t)∥ = t∥P_{h}u_{t}(t) − u_{t}(t)∥ ≤ Ct∥u_{t}(t)∥ ≤ C∥v∥ ≤ C_{δ}t^{−δ}∥v∥.

For η we note, with ρ = R_{h}u − u,

η_{t} + A_{h}η = A_{h}(R_{h} − P_{h})u = A_{h}P_{h}ρ.

Differentiating and setting ω = tη_{t} we have

ω_{t} + A_{h}ω = A_{h}P_{h}(tρ_{t}) + η_{t},

so that, since ω(0) = 0 and η(0) = 0,

ω(t) = ∫_{0}^{t} E_{h}(t − s)(A_{h}P_{h}(sρ_{t}(s)) + η_{t}(s)) ds
= η(t) + ∫_{0}^{t} A_{h}E_{h}(t − s)(P_{h}(sρ_{t}(s)) − η(s)) ds.

We recall from [3] that, since (α + β)^{−1} ≤ α^{−1+γ}β^{−γ}, for α, β > 0, 0 < γ < 1, we have

∥A_{h}E_{h}(t)∥ ≤ C log \frac{1}{h}(t + h^{2})^{−1} ≤ Ch^{−2(1−γ)} log \frac{1}{h}t^{−γ}.

Further, using the stability of P_{h}, the logarithmic stability of R_{h}, and the analyticity of E(t) in L_{p}(Ω), we have for a suitable ν < ∞,

s∥P_{h}ρ_{t}(s)∥ ≤ Cs∥ρ_{t}(s)∥ ≤ C log \frac{1}{h}s inf_{χ∈S_{h}} ∥u_{t}(s) − χ∥

≤ C log \frac{1}{h}h^{2−3ϵ}s∥u_{t}(s)∥_{W_{∞}^{2−3ϵ}} ≤ C log \frac{1}{h}h^{2−3ϵ}s∥u_{t}(s)∥_{W_{ν}^{2−2ϵ}}

≤ C log \frac{1}{h}h^{2−3ϵ}s^{−1+ϵ}∥v∥_{L_{ν}} ≤ C log \frac{1}{h}h^{2−3ϵ}s^{−1+ϵ}∥v∥.

It is easy to show η(t) ≤ C_{δ}t^{−δ} and, moreover, from [3, (3.17)] we have

∥η(t)∥ ≤ Ch^{2−3ϵ}s^{−1+ϵ}∥v∥,

and we conclude

\begin{aligned}∥ω(t)∥ & ≤ C_{δ}t^{−δ}∥v∥ + Ch^{2γ−3ϵ}( log \frac{1}{h})^{2} ∫_{0}^{t}(t − s)^{−γ}s^{−1+ϵ} ds∥v∥ \\ & = C_{δ}t^{−δ}∥v∥ + Ch^{2γ−4ϵ}t^{ϵ−γ}∥v∥ ≤ C_{δ}t^{−δ}∥v∥, \end{aligned}

by the choice ϵ = δ, γ = 2δ. This completes the proof. □

We now show the analogue of Theorem 4.1 in the case of the maximum norm.

Theorem 4.2. With respect to the maximum norm we have under the above assumptions, for the solutions of (1.1) and (4.2) with v_{h} = R_{h}v, for h ≤ h_{0}, k ≤ k_{0},

∥U^{n} − u(t_{n})∥ ≤ C(T, u) log \frac{1}{h}(h^{2} + k), for t_{n} ≤ T.

Proof. We follow the lines of the proof of Theorem 4.1. In maximum norm we have

∥ρ^{n}∥ ≤ Ch^{2} log \frac{1}{h}∥u^{n}∥_{W_{∞}^{2}}.

Similarly to above we now get

∥τ_{1}^{n}∥ ≤ C(u)h^{2} log \frac{1}{h}, ∥τ_{2}^{n}∥ ≤ C(u)k.

Further, since Δ_{h}R_{h} = P_{h}Δ,

∥τ_{3}^{n}∥ = ∥ϵ^{n}(b_{n}Δ_{h}R_{h}u)∥ ≤ Ck ∫_{0}^{t_{n}} ∑_{l=1}^{p} ∥Δu^{(l)}∥ dt ≤ C(u)k.

In this case Z_{h,n} = −b(t_{n}, ·)(Δ_{h}R_{h} − P_{h}Δ) = 0 and hence τ_{4}^{n} = 0. Theorem 2.4 with f^{n} = ∑_{i=1}^{4} τ_{i}^{n}, g^{n} = τ_{1}^{n}, therefore shows (recall that now M = C log 1/h)

∥θ^{n}∥ ≤ C(T, u) log \frac{1}{h}(h^{2} + k). □

5. Quadrature rules. We now give some examples of quadrature rules satisfying our assumptions. As mentioned in the introduction, the most obvious choice is the rectangle rule, which corresponds to taking all ω_{ns} = k, for s < n. Clearly then (1.6) holds with p = 1 and the sum in (1.7) vanishes.

A drawback of this method is that all the previously computed values of the solution enter into the equation (1.5), so that all of these have to be stored for future use. Following the philosophy of [4] and [6] we shall now turn to some sparse rules, that reduce the storage requirement.

We begin with a quadrature rule based on the trapezoidal rule on intervals of lengths O(k^{1/2}), with a slight modification near t_{n}. Let m = [k^{−1/2}], set k_{1} = mk and t̄_{j} = jk_{1}, and let j_{n} be the largest integer with t̄_{j_{n}} < t_{n}. For the interval (0, t_{n}) we then apply the composite trapezoidal rule with stepsize k_{1} on (0, t̄_{j_{n}}), then the one-interval trapezoidal rule on (t̄_{j_{n}}, t_{n−1}), and finally the left side rectangle rule on (t_{n−1}, t_{n}). Thus

σ^{n}(φ) = \frac{k_{1}}{2} ∑_{j=1}^{j_{n}}(φ(t̄_{j}) + φ(t̄_{j−1})) + \frac{1}{2}(t_{n−1} − t̄_{j_{n}})(φ(t_{n−1}) + φ(t̄_{j_{n}})) + kφ(t_{n−1}).

Since the rule is second order in k_{1} over (0, t̄_{j_{n}}) and (t̄_{j_{n}}, t_{n−1}), and first order on (t_{n−1}, t_{n}), (1.6) holds with p = 2. Here ω_{ns} ≤ k^{1/2} for s < n, and it is easy to see that (1.7) holds with q = 1/2. The number of time levels that enter the computation is of order O(k^{−1/2}) for this rule, as compared with O(k^{−1}) for the rectangle rule.

In [4] a similar quadrature rule was used, with the difference that the left side composite rectangle rule with stepsize k was used on (t̄_{j_{n}}, t_{n}). Again (1.6) and (1.7) hold with p = 2, q = 1/2. The reason for using the rectangle rule on (t̄_{j_{n}}, t_{n}) in this way was that here the ω_{ns} are “dominated weights” in the sense described before Lemma 2.2. In this example, we may take ω_{s} = O(k_{1}) = O(k^{1/2}) for s divisible by m and ω_{s} = O(k) for all other s. Even though the storage requirement is increased compared to the above method, it is still of the same order O(k^{−1/2}). We remark that the first sparse rule described above does not have dominated quadrature weights, since ω_{s,s−1} ≥ ck^{1/2} with c > 0 for all t_{s} in the right hand halves (say) of the intervals (t̄_{j−1}, t̄_{j}), so that ω_{s} ≥ ω_{ns} implies ∑_{s=0}^{n−1} ω_{s} ≥ \frac{1}{2}ck^{−1/2}.

Going one step further with the idea of reducing the storage requirement, we may set m = [k^{−1/4}] and k_{2} = m^{3}k = O(k^{1/4}), and do the following. We first use Simpson’s rule on as many intervals of lengths 2k_{2} that can be fitted into [0, t_{n−1}), and then, on the remaining interval, which is of length at most O(k^{1/2}), the composite trapezoidal rule on as many intervals of lengths k_{1} = m^{2}k = O(k^{1/2}) as fit in, thus reaching t̄_{j_{n}}, then the one-interval trapezoidal rule on the interval (t̄_{j_{n}}, t_{n−1}), and finally the left rectangle rule on (t_{n−1}, t_{n}). Similarly to above, (1.6) and (1.7) hold with p = 4, q = 1/4, and the number of time-levels that need to be stored per unit time is now O(k_{2}^{−1}) + O(k_{2}k_{1}^{−1}) + 1 = O(k^{−1/4}). This rule does not have dominated weights. Thus, our present assumptions allow some advantageous rules that were not covered in [4] or [6].

References

1. M. Crouzeix, S. Larsson, and V. Thomée, Resolvent estimates for elliptic finite element operators in one dimension, Math. Comp. 63 (1994), 121–140.

2. A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Springer-Verlag, Berlin and New York, 1983.

3. A. H. Schatz, V. Thomée and L. B. Wahlbin, Maximum norm stability and error estimates in parabolic finite element equations, Comm. Pure Appl. Math. 33 (1980), 265–304.

4. I. H. Sloan and V. Thomée, Time discretization of an integro-differential equation of parabolic type, SIAM J. Numer. Anal. 23 (1986), 1052–1061.

5. V. Thomée and N.-Y. Zhang, Error estimates for semidiscrete finite element methods for parabolic integro-differential equations, Math. Comp. 53 (1989), 121–139.

6. N.-Y. Zhang, On fully discrete Galerkin approximations for partial integro-differential equations of parabolic type, Math. Comp. 60 (1993), 133–166.

N. Yu. Bakaev, Department of Mathematics, Air Force Engineering Academy, Leningradskii Prospekt 40, Moscow 125190, Russia

   E-mail address: bakaev@postman.ru

S. Larsson and V. Thomée, Department of Mathematics, Chalmers University of Technology and Göteborg University, SE–412 96 Göteborg, Sweden

   E-mail address: stig@math.chalmers.se, thomee@math.chalmers.se