NUMERICAL ANALYSIS OF SEMILINEAR PARABOLIC PROBLEMS

STIG LARSSON

Abstract. In these lectures I discuss error analysis techniques for finite element methods for systems of reaction-diffusion equations with applications in dynamical systems theory. The emphasis is on pedagogical aspects and analysis techniques rather than on results. The list of techniques discussed include: analytic semigroup, parabolic smoothing, non-smooth data error estimate, a priori error estimate, a posteriori error estimate, exponential dichotomy, shadowing.

Contents

1. The continuous problem 1 2. Local a priori error estimates 7 2.1. The spatially semidiscrete problem 7 2.2. A completely discrete scheme 10 3. Shadowing—first approach 11 3.1. Linearization 11 3.2. Exponential dichotomies 13 3.3. Shadowing 16 4. A posteriori error estimates 20 4.1. The error equation 20 4.2. Local estimates of the residual 23 4.3. A global error estimate 25 5. Shadowing—second approach 26 References 29

1. The continuous problem

We consider the following initial-boundary value problem for a system of reactiondiffusion equations,

\begin{aligned}u_{t} − ϵΔu & = f̃(u), & x ∈ Ω, t > 0, \\ u & = 0, & x ∈ ∂Ω, t > 0, \\ u(·, 0) & = u_{0}, & x ∈ Ω, \end{aligned} (1.1)

Lectures delivered at the EPSRC VIII-th Summer School in Numerical Analysis, University of Leicester, July 12–17, 1998. Published in “The Graduate Student’s Guide to Numerical Analysis ’98”, editors M. Ainsworth, J. Levesley, and M. Marletta, SSCM volume 26, Springer-Verlag, 1999, pp. 83–117.

where Ω is a bounded domain in R^{d}, d = 1, 2, 3, u = u(x, t) ∈ R^{p}, u_{t} = ∂u/∂t, Δu = ∑_{i=1}^{d} ∂^{2}u/∂x_{i}^{2}, ϵ = diag(ϵ_{1}, …, ϵ_{p}) is a diagonal matrix of constant coefficients ϵ_{j} > 0, and f̃: R^{p} → R^{p} is twice continuously differentiable. If d = 2, 3 we assume, in addition, that the derivatives of f̃ satisfy the growth condition

|f̃^{(l)}(ξ)| ≤ C(1 + |ξ|^{δ+1−l}), ξ ∈ R^{p}, l = 1, 2, (1.2)

where | · | denotes the Euclidean norm on R^{p} and the induced operator norms, and where δ = 2 if d = 3, δ ∈ [1, ∞) if d = 2.

We assume that Ω is a convex polygon, so that we have access to the elliptic regularity theory and so that finite element meshes can be fitted exactly to the domain.

In the sequel we use the Hilbert space H = (L_{2}(Ω))^{p}, with its standard norm and inner product

∥v∥ = ( ∫_{Ω} |v|^{2} dx)^{1/2}, (v, w) = ∫_{Ω} v · w dx. (1.3)

The norms in the Sobolev spaces (H^{m}(Ω))^{p}, m ≥ 0, are denoted by

∥v∥_{m} = ( ∑_{|α|≤m} ∥D^{α}v∥^{2})^{1/2}. (1.4)

The space V = (H_{0}^{1}(Ω))^{p}, with norm ∥·∥_{1}, consists of the functions in (H^{1}(Ω))^{p} that vanish on ∂Ω. V^{∗} = (H^{−1}(Ω))^{p} is the dual space of V with norm

∥v∥_{−1} = sup_{χ∈V} \frac{|(v, χ)|}{∥χ∥_{1}}. (1.5)

If X, Y are Banach spaces, then L(X, Y) denotes the space of bounded linear operators from X into Y, L(X) = L(X, X), and B_{X}(x, R) denotes the closed ball in X with center x and radius R. In particular, we let B_{R} = B_{V} (0, R) denote the the closed ball of radius R in V:

B_{R} = {v ∈ V : ∥v∥_{1} ≤ R}.

We also use the notation

∥v∥_{L_{∞}([0,T],X)} = sup_{t∈[0,T]} ∥v(t)∥_{X}.

We set the problem up in the framework of [6]. We define the unbounded operator A = −ϵΔ on H with domain of definition D(A) = (H^{2}(Ω) ∩ H_{0}^{1}(Ω))^{p}. Then A is a closed, densely defined, and self-adjoint positive definite operator in H with compact inverse. Moreover, our assumption (1.2) guarantees that the mapping f̃ induces a nonlinear operator f : V → H through f(v)(x) = f̃(v(x)), see Lemma 1.1 below. The initial-boundary value problem (1.1) may then be formulated as an initial value problem in V: find u(t) ∈ V such that

u′ + Au = f(u), t > 0; u(0) = u_{0}. (1.6)

The operator −A is the infinitesimal generator of the analytic semigroup E(t) = exp(−tA) defined by

E(t)v = ∑_{j=1}^{∞} e^{−tλ_{j}}(v, φ_{j})φ_{j}, v ∈ H, (1.7)

where λ_{j} and φ_{j} denote the eigenvalues and an orthonormal basis of eigenvectors of A. The semigroup E(t) is the solution operator of the initial value problem for the homogeneous equation,

u′ + Au = 0, t > 0; u(0) = u_{0}.

Its solution is thus given by u(t) = E(t)u_{0}. By Duhamel’s principle it follows that solutions of (1.6) satisfy

u(t) = E(t)u_{0} + ∫_{0}^{t} E(t − s)f(u(s)) ds, t ≥ 0. (1.8)

Conversely, appropriately defined solutions of the nonlinear integral equation (1.8) are solutions of the differential equation (1.6); see Theorem 1.3 below. We shall mainly work with (1.8) and discretized variants of it.

An important ingredient in our framework is that f can be controlled by fractional powers of A. We define the fractional powers of A by means of the spectral theorem and we have

\begin{aligned}∥A^{α}v∥ & = ( ∑_{j=1}^{∞} (λ_{j}^{α}(v, φ_{j}))^{2})^{1/2}, α ∈ R, \\ D(A^{α}) & = \{v : ∥A^{α}v∥ < ∞\}, α ∈ R. \end{aligned} (1.9)

Using also the elliptic regularity estimate,

∥v∥_{2} ≤ C∥Av∥, v ∈ (H^{2}(Ω) ∩ H_{0}^{1}(Ω))^{p}, (1.10)

and the trace inequality,

∥v∥_{L_{2}(∂Ω)} ≤ C∥v∥_{1}, v ∈ (H^{1}(Ω))^{p},

we obtain D(A^{l/2}) = (H^{l}(Ω) ∩ H_{0}^{1}(Ω))^{p}, l = 1, 2, with the equivalence of norms

c∥v∥_{l} ≤ ∥A^{l/2}v∥ ≤ C∥v∥_{l}, v ∈ D(A^{l/2}), l = 1, 2; (1.11)

cf. [12]. A simple exercise using the special case l = 1 of (1.11) shows

c∥v∥_{−1} ≤ ∥A^{−1/2}v∥ ≤ C∥v∥_{−1}; (1.12)

cf. (2.7) below.

The analyticity of the semigroup is reflected in the inequalities (D_{t} = ∂/∂t)

∥D_{t}^{l}E(t)v∥ = ∥A^{l}E(t)v∥ ≤ C_{l}t^{−l}∥v∥, t > 0, v ∈ H, l ≥ 0. (1.13)

These follow easily from (1.7) and Parseval’s identity. Using also the norm equivalences (1.11), (1.12) we obtain the smoothing property

\begin{aligned}∥D_{t}^{l}E(t)v∥_{β} ≤ C_{l}t^{−l−(β−α)/2}∥v∥_{α}, & t > 0, v ∈ D(A^{α/2}), \\ & − 1 ≤ α ≤ β ≤ 2, l = 0, 1. \end{aligned} (1.14)

We need the fact that the operator f : V → H satisfies a local Lipschitz condition. This is contained in our first lemma. We also prove that the Fréchet derivative f′ : V → L(V, H) satisfies a local Lipschitz condition. The proof is based on the assumption (1.2) and the Sobolev inequality (where p = 6 if d ≤ 3, p < ∞ if d = 2, and p = ∞ if d = 1)

∥v∥_{L_{p}} ≤ C∥v∥_{1}, (1.15) and on Hölder’s inequality

∥v^{δ}w∥_{L_{r}} ≤ ∥v∥_{L_{q}}^{δ}∥w∥_{L_{p}}, \frac{δ}{q} + \frac{1}{p} = \frac{1}{r}, δ > 0. (1.16)

Lemma 1.1. For each nonnegative number R there is a constant C(R) such that, for all u, v, w ∈ B_{R}, l = 0, 1,

(1.17)  ∥f′(u)∥_{L(V,H)} ≤ C(R),
(1.18)  ∥f′(u)∥_{L(H,V^{∗})} ≤ C(R),
(1.19)  ∥f(u) − f(v)∥ ≤ C(R)∥u − v∥_{1},
(1.20)  ∥f(u) − f(v)∥_{−1} ≤ C(R)∥u − v∥,
(1.21)  ∥f′(u) − f′(v)∥_{L(V,H)} ≤ C(R)∥u − v∥_{1},
(1.22)  ∥f′(u) − f′(v)∥_{L(H,V^{∗})} ≤ C(R)∥u − v∥_{1},
(1.23) ∥f(u) − f(v) − f′(w)(u − v)∥_{l−1} ≤ C(R)(∥u − w∥_{1} + ∥v − w∥_{1})∥u − v∥_{l}.

If, in addition, u ∈ (H^{2}(Ω))^{p} and z ∈ (H^{2}(Ω))^{p} ∩ V, then

∥A^{1/2}(f′(u)z)∥ ≤ C(R)(∥z∥_{2} + ∥u∥_{2}∥z∥_{1}). (1.24)

Proof. We have in view of (1.2) and the Hölder and Sobolev inequalities, for z ∈ V,

∥f′(u)z∥_{L_{2}} ≤ C(1 + ∥u∥_{L_{q}}^{δ})∥z∥_{L_{p}}
≤ C(1 + ∥u∥_{1}^{δ})∥z∥_{1},

where \frac{1}{p} + \frac{δ}{q} = \frac{1}{2} with p = q = 6 if d = 3, and with arbitrary p ∈ (1, ∞) if d ≤ 2. This proves (1.17) and (1.19) follows. Similarly, for any z ∈ V,

∥(f′(u) − f′(v))z∥_{L_{2}} ≤ C(1 + ∥u∥_{L_{q}}^{δ−1} + ∥v∥_{L_{q}}^{δ−1})∥u − v∥_{L_{p}}∥z∥_{L_{q}}
≤ C(1 + ∥u∥_{1}^{δ−1} + ∥v∥_{1}^{δ−1})∥u − v∥_{1}∥z∥_{1},

with the same p and q as before. This proves (1.21).

Moreover, for any z, χ ∈ V,

(f′(u)z, χ) ≤ C(1 + ∥u∥_{L_{q}}^{δ})∥z∥_{L_{2}}∥χ∥_{L_{p}}
≤ C(1 + ∥u∥_{1}^{δ})∥z∥ ∥χ∥_{1},

and

((f′(u) − f′(v))z, χ) ≤ C(1 + ∥u∥_{L_{q}}^{δ−1} + ∥v∥_{L_{q}}^{δ−1})∥u − v∥_{L_{2}}∥z∥_{L_{q}}∥χ∥_{L_{p}}
≤ C(1 + ∥u∥_{1}^{δ−1} + ∥v∥_{1}^{δ−1})∥u − v∥∥z∥_{1}∥χ∥_{1},

where \frac{δ}{q} + \frac{1}{2} + \frac{1}{p} = 1, i.e., with the same p and q as before. This proves (1.18) and (1.22); (1.20) follows from (1.18).

Next, (1.23) is obtained by applying (1.21) and (1.22) to the identity

f(u) − f(v) − f′(w)(u − v) = ∫_{0}^{1} (f′(su + (1 − s)v) − f′(w)) ds (u − v).

Finally, (1.24) is proved in a similar way, by using the equivalence of norms (1.11) and computing the first order partial derivatives of f′(u)z. Note that this uses the fact that f′(u)z ∈ V, i.e., f′(u)z = 0 on ∂Ω.

We may now prove local existence of solutions of (1.8) and hence of (1.6).

Theorem 1.2. For any R_{0} > 0 there is τ = τ(R_{0}) such that (1.8) has a unique solution u ∈ C([0, τ], V) for any initial value u_{0} ∈ V with ∥u_{0}∥_{1} ≤ R_{0}. Moreover, there is c such that ∥u∥_{L_{∞}([0,τ],V)} ≤ cR_{0}.

Proof. Let u_{0} ∈ B_{R_{0}}, define

S(u)(t) = E(t)u_{0} + ∫_{0}^{t} E(t − s)f(u(s)) ds,

and note that (1.8) is a fixed point equation u = S(u). We shall choose τ and R such that we can apply Banach’s fixed point theorem (the contraction mapping theorem) in the closed ball

B = {u ∈ C([0, τ], V) : ∥u∥_{L_{∞}([0,τ],V)} ≤ R}

in the Banach space C([0, τ], V).

We must show (i) that S maps B into itself, (ii) that S is a contraction on B. In order to prove (i) we take u ∈ B and first note that the Lipschitz condition (1.19) implies that

\begin{aligned}∥f(u(t))∥ & ≤ ∥f(0)∥ + ∥f(u(t)) − f(0)∥ \\ & ≤ ∥f(0)∥ + C(R)∥u(t)∥_{1} \\ & ≤ ∥f(0)∥ + C(R)R, 0 ≤ t ≤ τ. \end{aligned} (1.25)

Hence, using also (1.14), we get

∥S(u)(t)∥_{1} ≤ ∥E(t)u_{0}∥_{1} + ∫_{0}^{t} ∥E(t − s)f(u(s))∥_{1} ds
≤ c_{0}∥u_{0}∥_{1} + c_{1} ∫_{0}^{t}(t − s)^{−1/2}∥f(u(s))∥ ds
≤ c_{0}R_{0} + 2c_{1}τ^{1/2}(∥f(0)∥ + C(R)R), 0 ≤ t ≤ τ.

This implies

∥S(u)∥_{L_{∞}([0,τ],V)} ≤ c_{0}R_{0} + 2c_{1}τ^{1/2}(∥f(0)∥ + C(R)R).

Choose R = 2c_{0}R_{0} and τ = τ(R_{0}) so small that

2c_{1}τ^{1/2}(∥f(0)∥ + C(R)R) ≤ \frac{1}{2}R. (1.26)

Then ∥S(u)∥_{L_{∞}([0,τ],V)} ≤ R and we conclude that S maps B into itself. To show (ii) we take u, v ∈ B and note that

∥f(u(t)) − f(v(t))∥ ≤ C(R)∥u − v∥_{L_{∞}([0,τ],V)}, 0 ≤ t ≤ τ.

Hence

∥S(u)(t) − S(v(t))∥_{1} ≤ ∫_{0}^{t} ∥E(t − s)(f(u(s)) − f(v(s)))∥_{1} ds
≤ c_{1} ∫_{0}^{t}(t − s)^{−1/2}∥f(u(s)) − f(v(s))∥ ds
≤ 2c_{1}τ^{1/2}C(R)∥u − v∥_{L_{∞}([0,τ],V)}, 0 ≤ t ≤ τ,

so that

∥S(u) − S(v)∥_{L_{∞}([0,τ],V)} ≤ 2c_{1}τ^{1/2}C(R)∥u − v∥_{L_{∞}([0,τ],V)}.

It follows from (1.26) that 2c_{1}τ^{1/2}C(R) ≤ \frac{1}{2} and we conclude that S is a contraction on B. Hence S has a unique fixed point u ∈ B.

The initial value problem (1.8) (or (1.6)) thus has a unique local solution for any initial datum u_{0} ∈ V. We denote by S(t, ·) the corresponding (local) solution operator, so that u(t) = S(t, u_{0}) is the solution of (1.8).

The following theorem provides regularity estimates for solutions u of (1.8). These will be used in our error analysis, but the theorem also shows that u′(t) ∈ H and Au(t) ∈ H for t > 0, so that any solution of the integral equation (1.8) is also a solution of the differential equation (1.6).

Theorem 1.3. Let R ≥ 0 and τ > 0 be given and let u ∈ C([0, τ], V) be a solution of (1.8). If ∥u(t)∥_{1} ≤ R for t ∈ [0, τ], then

(1.27)  ∥u(t)∥_{2} ≤ C(R, τ)t^{−1/2},  t ∈ (0, τ],
(1.28)  ∥u_{t}(t)∥_{s} ≤ C(R, τ)t^{−1−(s−1)/2},  t ∈ (0, τ], s = 0, 1, 2.

Proof. In view of (1.19) we have ∥f(u(t))∥ ≤ C(R) for t ∈ [0, τ], cf. (1.25). The proof of (1.27), (1.28) can now obtained by tracing the constants in [6, Theorem 3.5.2]. The argument is based on a generalization of Gronwall’s lemma, Lemma 1.4 below. The case s = 2 of (1.28) requires that we differentiate (1.6) with respect to t and use (1.24).

Lemma 1.4. (Generalized Gronwall lemma.) Let the function φ(t, τ) ≥ 0 be continuous for 0 ≤ τ < t ≤ T. If

φ(t, τ) ≤ A (t − τ)^{−1+α} + B ∫_{τ}^{t}(t − s)^{−1+β}φ(s, τ) ds, 0 ≤ τ < t ≤ T,

for some constants A, B ≥ 0, α, β > 0, then there is a constant C = C(B, T, α, β) such that

φ(t, τ) ≤ CA (t − τ)^{−1+α}, 0 ≤ τ < t ≤ T.

Proof. Iterating the given inequality N − 1 times, using the identity

∫_{τ}^{t}(t − s)^{−1+α}(s − τ)^{−1+β} ds = C(α, β) (t − τ)^{−1+α+β}, α, β > 0, (1.29)

(Abel’s integral) and estimating (t − τ)^{β} by T^{β}, we obtain

φ(t, τ) ≤ C_{1}A (t − τ)^{−1+α} + C_{2} ∫_{τ}^{t}(t − s)^{−1+Nβ}φ(s, τ) ds, 0 ≤ τ < t ≤ T,

where C_{1} = C_{1}(B, T, α, β, N), C_{2} = C_{2}(B, β, N). We now choose the smallest N such that −1 + Nβ ≥ 0, and estimate (t − s)^{−1+Nβ} by T^{−1+Nβ}. If −1 + α ≥ 0 we obtain the desired conclusion by the standard version of Gronwall’s lemma. Otherwise we set ψ(t, τ) = (t − τ)^{1−α}φ(t, τ) to obtain

ψ(t, τ) ≤ C_{1}A + C_{3} ∫_{τ}^{t}(s − τ)^{−1+α}ψ(s, τ) ds, 0 ≤ τ < t ≤ T,

and the standard Gronwall lemma yields ψ(t, τ) ≤ CA for 0 ≤ τ < t ≤ T, which is the desired result.

Note that the constant in Gronwall’s lemma grows exponentially with the length T of the time interval. Hence, results derived by means of this lemma are often useful only for short time intervals. There is also a discrete version of Lemma 1.4; see [3, Lemma 7.1].

We conclude this section with a remark about the abstract framework. It is tempting to start by saying “Let X be a Banach space and A the generator of an analytic semigroup in X.” This approach goes a rather long way with the continuous problem (see [6]) and with pure time discretization. But when it comes to spatial discretization we soon need to know that, for example, X = L_{2}(Ω) and A is an elliptic operator of second order. “PDE theory is not a branch of functional analysis” as L. C. Evans puts it in the preface of his textbook [5].

2. Local a priori error estimates

In this section we first discretize (1.1) with respect to the spatial variables by means of a standard piecewise linear finite element method. We then briefly discuss completely discrete approximation by means of the backward Euler time-stepping. This material is taken from [7]. A general reference is [12].

2.1. The spatially semidiscrete problem. The weak formulation of (1.1) is: find u(t) ∈ V such that

\begin{aligned}(u′, v) + a(u, v) & = (f(u), v), ∀v ∈ V, t > 0, \\ u(0) & = u_{0}, \end{aligned} (2.1)

where a(u, v) = (ϵ∇u, ∇v) = (−ϵΔu, v) = (Au, v) is the bilinear form associated with A.

Let {V_{h}}_{0<h<1} be a family of finite dimensional subspaces of V, where each V_{h} consists of continuous piecewise polynomials of degree ≤ 1 with respect to a triangulation T of Ω with maximal mesh size h, see [2]. The approximate solution u_{h}(t) ∈ V_{h} of (1.1) is defined by

\begin{aligned}(u′_{h}, χ) + a(u_{h}, χ) & = (f(u_{h}), χ), ∀χ ∈ V_{h}, t > 0, \\ u_{h}(0) & = u_{h,0}, \end{aligned} (2.2)

where u_{h,0} ∈ V_{h} is an approximation of u_{0}.

Introducing the linear operator A_{h} : V_{h} → V_{h} and the orthogonal projection P_{h} : H → V_{h}, defined by

(A_{h}ψ, χ) = a(ψ, χ), (P_{h}g, χ) = (g, χ) ∀ψ, χ ∈ V_{h}, g ∈ H, (2.3)

we may write (2.2) as

u′_{h} + A_{h}u_{h} = P_{h}f(u_{h}), t > 0; u_{h}(0) = u_{h,0}. (2.4)

The operator A_{h} is self-adjoint positive definite (uniformly in h); the corresponding semigroup E_{h}(t) = exp(−tA_{h}) : V_{h} → V_{h} therefore satisfies inequalities analogous to (1.13) (uniformly in h):

∥D_{t}^{l}E_{h}(t)v∥ = ∥A_{h}^{l}E_{h}(t)v∥ ≤ C_{l}t^{−l}∥v∥, t > 0, v ∈ V_{h}, l ≥ 0.

Moreover, for the operator A_{h} we have the equivalence of norms (cf. (1.11))

c∥v∥_{1} ≤ ∥A_{h}^{1/2}v∥ = \sqrt{a(v, v)} = ∥A^{1/2}v∥ ≤ C∥v∥_{1}, v ∈ V_{h}. (2.5)

We also have

∥P_{h}f∥ ≤ ∥f∥, f ∈ H,

and (cf. (1.12))

∥A_{h}^{−1/2}P_{h}f∥ ≤ C∥f∥_{−1}, f ∈ H, (2.6) which is obtained by using (2.5) in the calculation

∥A_{h}^{−1/2}P_{h}f∥ = sup_{v_{h}∈V_{h}} \frac{|(A_{h}^{−1/2}P_{h}f, v_{h})|}{∥v_{h}∥} = sup_{v_{h}∈V_{h}} \frac{|(f, A_{h}^{−1/2}v_{h})|}{∥v_{h}∥}

= sup_{w_{h}∈V_{h}} \frac{|(f, w_{h})|}{∥A_{h}^{1/2}w_{h}∥} ≤ C sup_{w_{h}∈V_{h}} \frac{|(f, w_{h})|}{∥w_{h}∥_{1}}

≤ C sup_{w∈V} \frac{|(f, w)|}{∥w∥_{1}} = C∥f∥_{−1}.

(2.7)

Using these inequalities we easily prove the smoothing property of E_{h}(t)P_{h}:

\begin{aligned}∥D_{t}^{l}E_{h}(t)P_{h}f∥_{β} ≤ Ct^{−l−(β−α)/2}∥f∥_{α}, & t > 0, f ∈ D(A^{α/2}), \\ & − 1 ≤ α ≤ β ≤ 1, l = 0, 1. \end{aligned} (2.8)

Note that the upper limit to β is 1, while it is 2 in the continuous case (1.14). The initial-value problem (2.4) is equivalent to the integral equation

u_{h}(t) = E_{h}(t)u_{h,0} + ∫_{0}^{t} E_{h}(t − s)P_{h}f(u_{h}(s)) ds, t ≥ 0. (2.9)

Because of (2.8) the proof of Theorem 1.2 carries over verbatim to the semidiscrete case. We thus have:

Theorem 2.1. For any R_{0} > 0 there is τ = τ(R_{0}) such that (2.9) has a unique solution u_{h} ∈ C([0, τ], V) for any initial value u_{h,0} ∈ V_{h} with ∥u_{h,0}∥_{1} ≤ R_{0}. Moreover, there is c such that ∥u_{h}∥_{L_{∞}([0,τ],V)} ≤ cR_{0}.

We denote by S_{h}(t, ·) the corresponding (local) solution operator, so that u_{h}(t) = S_{h}(t, u_{h,0}) is the solution of (2.9).

We may now estimate the difference between the local solutions u(t) = S(t, u_{0}) and u_{h}(t) = S_{h}(t, u_{h,0}) that we have obtained. We refer to the following result as a local a priori error estimate. It is local because the constant C(R, τ) grows exponentially with the length τ of the time interval and also because it grows with the size R of the solutions. It is a priori because the error is evaluated in terms of derivatives of u, which are estimated a priori in Theorem 1.3.

Theorem 2.2. Let R ≥ 0 and τ > 0 be given. Let u(t) and u_{h}(t) be solutions of (2.1) and (2.2) respectively, such that u(t), u_{h}(t) ∈ B_{R} for t ∈ [0, τ]. Then

(2.10) ∥u_{h}(t) − u(t)∥_{1} ≤ C(R, τ)t^{−1/2}(∥u_{h,0} − P_{h}u_{0}∥ + h),  t ∈ (0, τ],
(2.11)  ∥u_{h}(t) − u(t)∥ ≤ C(R, τ)(∥u_{h,0} − P_{h}u_{0}∥ + h^{2}t^{−1/2}),  t ∈ (0, τ].

Proof. We introduce the elliptic projection operator R_{h} : V → V_{h} defined by

a(R_{h}v, χ) = a(v, χ), ∀χ ∈ V_{h}. (2.12)

Under the usual regularity assumptions on the triangulation and the elliptic regularity estimate (1.10), we have the error bounds (see [2])

∥R_{h}v − v∥ + h∥R_{h}v − v∥_{1} ≤ Ch^{s}∥v∥_{s}, v ∈ V ∩ (H^{s}(Ω))^{p}, s = 1, 2, (2.13)

Following a standard practice we divide the error into two parts:

e(t) ≡ u_{h}(t) − u(t) = (u_{h}(t) − R_{h}u(t)) + (R_{h}u(t) − u(t)) ≡ θ(t) + ρ(t).

In view of (2.13) and (1.27), (1.28) we have, for j = 0, 1 and s = 1, 2,

(2.14)  ∥ρ(t)∥_{j} ≤ Ch^{s−j}∥u(t)∥_{s} ≤ C(R, τ)h^{s−j}t^{−(s−1)/2},  t ∈ (0, τ],
(2.15)  ∥ρ_{t}(t)∥ ≤ Ch^{s}∥u_{t}(t)∥_{s} ≤ C(R, τ)ht^{−1−(s−1)/2},  t ∈ (0, τ].

It remains to estimate θ(t), which belongs to V_{h}. In view of (2.2), the identity A_{h}R_{h} = P_{h}A (which follows easily from (2.3) and (2.12)), and (1.6), we find that

θ_{t} + A_{h}θ = P_{h}(f(u_{h}) − f(u) − ρ_{t}). (2.16)

Hence

θ(t) = E_{h}(t)θ(0) + ∫_{0}^{t} E_{h}(t − σ)P_{h}(f(u_{h}(σ)) − f(u(σ)) − D_{σ}ρ(σ)) dσ, (2.17)

where D_{σ} = ∂/∂σ. Integration by parts yields

\begin{aligned}− ∫_{0}^{t/2} E_{h}(t − σ)P_{h}D_{σ}ρ(σ) dσ & = E_{h}(t)P_{h}ρ(0) − E_{h}(t/2)P_{h}ρ(t/2) \\ & + ∫_{0}^{t/2} (D_{σ}E_{h}(t − σ))P_{h}ρ(σ) dσ. \end{aligned}

Hence

\begin{aligned}θ(t) & = E_{h}(t)P_{h}e(0) − E_{h}(t/2)P_{h}ρ(t/2) + ∫_{0}^{t/2} (D_{σ}E_{h}(t − σ))P_{h}ρ(σ) dσ \\ & − ∫_{t/2}^{t} E_{h}(t − σ)P_{h}D_{σ}ρ(σ) dσ + ∫_{0}^{t} E_{h}(t − σ)P_{h}(f(u_{h}(σ)) − f(u(σ))) dσ. \end{aligned}

Using the smoothing property (2.8) of E_{h}(t)P_{h}, the error estimates (2.14), (2.15) with j = 0, s = 1, and the Lipschitz condition (1.19) we obtain

∥θ(t)∥_{1} ≤ Ct^{−1/2}(∥P_{h}e(0)∥ + ∥ρ(t/2)∥) + C ∫_{0}^{t/2}(t − σ)^{−3/2}∥ρ(σ)∥ dσ
+ C ∫_{t/2}^{t}(t − σ)^{−1/2}∥D_{σ}ρ(σ)∥ dσ

+ C ∫_{0}^{t}(t − σ)^{−1/2}∥f(u_{h}(σ)) − f(u(σ))∥ dσ
≤ C(R, τ)t^{−1/2}(∥P_{h}e(0)∥ + h)
+ C(R, τ)h(∫_{0}^{t/2}(t − σ)^{−3/2} dσ + ∫_{t/2}^{t}(t − σ)^{−1/2}σ^{−1} dσ)

+ C(R) ∫_{0}^{t}(t − σ)^{−1/2}∥e(σ)∥_{1} dσ
≤ C(R, τ)t^{−1/2}(∥P_{h}e(0)∥ + h) + C(R) ∫_{0}^{t}(t − σ)^{−1/2}∥e(σ)∥_{1} dσ,

for t ∈ (0, τ]. Together with (2.14), (2.15) and e = θ + ρ this yields

∥e(t)∥_{1} ≤ C(R, τ)t^{−1/2}(∥P_{h}e(0)∥ + h) + C(R) ∫_{0}^{t}(t − σ)^{−1/2}∥e(σ)∥_{1} dσ, t ∈ (0, τ],

and the desired bound follows by the generalized Gronwall lemma, cf. Lemma 1.4. This proves (2.10), because P_{h}e(0) = u_{h,0} − P_{h}u_{0}.

To prove (2.11) we use the Lipschitz condition (1.20) instead of (1.19). We omit the details.

We now formulate a consequence of the previous theorem, which is the form in which we will apply it later on. (It is convenient to change the notation so that the length of the local interval of existence is now 2τ.)

Theorem 2.3. For any R_{0} ≥ 0 there are numbers τ = τ(R_{0}) > 0 and C(R_{0}) such that if v ∈ B_{R_{0}}, v_{h} ∈ V_{h} ∩ B_{R_{0}}, then S(t, v) and S_{h}(t, v_{h}) exist for t ∈ [0, 2τ]. Moreover, for l = 0, 1,

∥S_{h}(t, v_{h}) − S(t, v)∥_{l} ≤ C(R_{0})t^{−1/2}(∥v_{h} − P_{h}v∥ + h^{2−l}), t ∈ (0, 2τ],

and

∥S_{h}(t, v_{h}) − S(t, v)∥_{l} ≤ C(R_{0})(∥v_{h} − P_{h}v∥ + h^{2−l}), t ∈ [τ, 2τ].

Proof. The local existence theorems 1.2 and 2.1 give S(t, v), S_{h}(t, v_{h}) ∈ B_{R}, t ∈ [0, 2τ], for some τ = τ(R_{0}) and with R = cR_{0}. We now apply Theorem 2.2.

2.2. A completely discrete scheme. In this subsection we show that the above program can be carried through also for a completely discrete scheme based on the backward Euler method. We replace the time derivative in (2.2) by a backward difference quotient ∂_{t}U_{j} = (U_{j} − U_{j−1})/k, where k is a time step and U_{j} is the approximation of u_{j} = u(t_{j}) and t_{j} = jk. The discrete solution U_{j} ∈ V_{h} thus satisfies

∂_{t}U_{j} + A_{h}U_{j} = P_{h}f(U_{j}), t_{j} > 0; U_{0} = u_{h,0}. (2.18)

Duhamel’s principle yields

U_{j} = E_{kh}^{j}u_{h,0} + k ∑_{l=1}^{j} E_{kh}^{j−l−1}P_{h}f(U_{l}), t_{j} ≥ 0, (2.19)

where E_{kh} = (I + kA_{h})^{−1}. Since A_{h} is self-adjoint positive definite (uniformly in h), we have the inequality

∥∂_{t}^{l}E_{kh}^{j}v∥ = ∥A_{h}^{l}E_{kh}^{j}v∥ ≤ C_{l}t_{j}^{−l}∥v∥, t_{j} ≥ t_{l}, v ∈ V_{h}, l ≥ 0,

which is the discrete analogue of (1.13). In view of the inequalities (2.5) and (2.6) this leads to a smoothing property analogous to (2.8):

\begin{aligned}∥∂_{t}^{l}E_{kh}^{j}P_{h}f∥_{β} ≤ Ct_{j}^{−l−(β−α)/2}∥f∥_{α}, & t_{j} > 0, f ∈ D(A^{α/2}), \\ & − 1 ≤ α ≤ β ≤ 1, l = 0, 1. \end{aligned} (2.20)

Again the proof of Theorem 1.2 carries over verbatim to the discrete case. We thus have:

Theorem 2.4. For any R_{0} > 0 there is τ = τ(R_{0}) such that (2.19) has a unique solution U_{j}, t_{j} ∈ [0, τ], for any initial value u_{h,0} ∈ V_{h} with ∥u_{h,0}∥_{1} ≤ R_{0}. Moreover, there is c such that max_{t_{j}∈[0,τ]} ∥U_{j}∥_{1} ≤ cR_{0}.

We denote by S_{hk}(t_{j}, ·) the corresponding (local) solution operator, so that U_{j} = S_{hk}(t_{j}, u_{h,0}) is the solution of (2.19).

We may now estimate the difference between the local solutions u(t) = S(t, u_{0}) and U_{j} = S_{hk}(t_{j}, u_{h,0}) that we have obtained. We refer to the following result as a local a priori error estimate. The proof can be found in [7]. It follows the lines of the proof of Theorem 2.2 but is more technical.

Theorem 2.5. Let R ≥ 0 and τ > 0 be given. Let u(t) and U_{j} be solutions of (2.1) and (2.19) respectively, such that u(t), U_{j} ∈ B_{R} for t, t_{j} ∈ [0, τ]. Then, for k ≤ k_{0}(R), we have

∥U_{j} − u(t_{j})∥_{1} ≤ C(R, τ)(∥u_{h,0} − P_{h}u_{0}∥t_{j}^{−1/2} + ht_{j}^{−1/2} + kt_{j}^{−1}),  t_{j} ∈ (0, τ],
∥U_{j} − u(t_{j})∥ ≤ C(R, τ)(∥u_{h,0} − P_{h}u_{0}∥ + h^{2}t_{j}^{−1/2} + kt_{j}^{−1/2}),  t_{j} ∈ (0, τ].

There is also an analogue of Theorem 2.3.

Theorem 2.6. For any R_{0} ≥ 0 there are numbers τ = τ(R_{0}) > 0, k_{0}(R_{0}) > 0 and C(R_{0}) such that if v ∈ B_{R_{0}}, v_{h} ∈ V_{h} ∩ B_{R_{0}}, then S(t, v) and S_{hk}(t_{j}, v_{h}) exist for t, t_{j} ∈ [0, 2τ]. Moreover, for k ≤ k_{0}(R_{0}), we have

∥S_{hk}(t_{j}, v_{h}) − S(t_{j}, v)∥_{l} ≤ C(R_{0})(∥v_{h} − P_{h}v∥ + h^{2−l} + k), t_{j} ∈ [τ, 2τ], l = 0, 1.

3. Shadowing—first approach

In this section we present a numerical shadowing result. It is taken essentially from [10] but we present it in a slightly more general form.

3.1. Linearization. Let ū ∈ C([0, T], V) with ∥ū(t)∥_{1} ≤ R, t ∈ [0, T], for some T and R. In order to linearize around ū we rewrite the differential equation in (1.6),

u′ + Au + B(t)u = F (t, u), (3.1)

where B(t) is the linear part of f and F (t, ·) is the nonlinear remainder:

\begin{aligned}B(t) & = −f′(ū(t)) ∈ L(V, H), \\ F (t, v) & = f(v) − f′(ū(t))v. \end{aligned} (3.2)

Let L(t, s), 0 ≤ s ≤ t ≤ T, denote the solution operator of the linearized homogeneous problem

v′ + Av + B(t)v = 0, t > s; v(s) = ϕ. (3.3)

Thus, v(t) = L(t, s)ϕ is the solution of (3.3). The operator w ↦ −B(t)w = f′(ū(t))w satisfies a local Lipschitz condition from V into H (uniformly in t); cf. Lemma 1.1. We therefore obtain local existence in C([0, τ], V) by the same argument as for (1.6). The following lemma provides a global a priori bound, which implies that any local solution can be extended to a global solution on the whole interval [0, T]. The lemma will also be useful in our analysis below. The bound grows exponentially with the length T of the time interval; this is realistic when the linearized problem is unstable, otherwise it is an overestimate. Note, by the way, that the linearized operator A + B(t) is not self-adjoint in general even though A is.

Lemma 3.1. Assume that ū ∈ C([0, T], V) with ∥ū(t)∥_{1} ≤ R, t ∈ [0, T], for some T and R. Then any solution v(t) = L(t, s)ϕ of (3.3) satisfies

∥L(t, s)ϕ∥_{m} ≤ C(R, T)(t − s)^{−(m−l)/2}∥ϕ∥_{l}, 0 ≤ s < t ≤ T, 0 ≤ l ≤ m ≤ 1.

Proof. By Duhamel’s principle we have

v(t) = E(t − s)ϕ + ∫_{s}^{t} E(t − σ)f′(ū(σ))v(σ) dσ,

so that, in view of (1.14), (1.17) and (1.18),

\begin{aligned}∥v(t)∥_{m} & ≤ C(t − s)^{−(m−l)/2}∥ϕ∥_{l} + C ∫_{s}^{t}(t − σ)^{−1/2}∥f′(ū(σ))v(σ)∥_{m−1} dσ \\ & ≤ C(t − s)^{−(m−l)/2}∥ϕ∥_{l} + C(R) ∫_{s}^{t}(t − σ)^{−1/2}∥v(σ)∥_{m} dσ. \end{aligned}

The desired bound now follows by Gronwall’s lemma 1.4.

In order to conveniently take advantage of the smoothing properties of our problems, as reflected in the error bounds of Theorems 2.3 and 2.6, we replace the continuous dynamical system u(t) = S(t, u_{0}), t ≥ 0, by a discrete dynamical system u_{n} = u(nτ) = S(τ, u_{n−1}), n = 1, 2, …, with an appropriately chosen step τ.

Let ∥ū(t)∥_{1} ≤ R, t ∈ [0, T], and let τ = τ(R+1) be half the time of local existence as in Theorems 2.3 and 2.6, adjusted so that Nτ = T for some positive integer N. This implies that the operator S = S(τ, ·) is defined on the ball D = B_{R+1}. We may thus think of S as a nonlinear operator S : D ⊂ X → X with both X = V and X = H; its domain D = B_{R+1} is a ball in V in both cases.

Let T_{n} = nτ and L_{n} = L(T_{n+1}, T_{n}). Duhamel’s principle applied to (3.1) gives

u(t) = S(t − T_{n}, v) = L(t, T_{n})v + ∫_{T_{n}}^{t} L(t, s)F (s, u(s)) ds. (3.4)

With t = T_{n+1} we get a decomposition of S = S(τ, ·),

S = L_{n} + N_{n},

where

N_{n}(v) = ∫_{T_{n}}^{T_{n+1}} L(t, s)F (s, S(s − T_{n}, v)) ds.

We think of L_{n} as the linear part of the nonlinear time-τ evolution S at time level T_{n}. In the following lemma we estimate the nonlinear remainder N_{n}.

Lemma 3.2. Let ρ ≤ 1. (a) If

∥S(t − T_{n}, v_{i}) − ū(t)∥_{1} ≤ ρ, t ∈ [T_{n}, T_{n+1}], i = 1, 2, (3.5)

then

∥N_{n}(v_{1}) − N_{n}(v_{2})∥_{1} ≤ C(R, τ)ρ∥v_{1} − v_{2}∥_{l}, l = 0, 1. (3.6)

(b) If ū is a solution of (1.6) and ∥v_{i} − ū(T_{n})∥_{1} ≤ ρ, i = 1, 2, then (3.6) holds.

Proof. Let u_{i}(t) = S(t − T_{n}, v_{i}) and w_{i}(t) = ∫_{T_{n}}^{t} L(t, s)F (s, u_{i}(s)) ds for i = 1, 2. Then N_{n}(v_{i}) = w_{i}(T_{n+1}). (Note that u_{i}(t) exists for t ∈ [T_{n}, T_{n+1}] because v_{i} ∈ D = B_{R+1}.) Applying Lemma 3.1 over the (short!) time interval [T_{n}, T_{n+1}] we get, with e = w_{1} − w_{2},

\begin{aligned}∥e(t)∥_{1} & ≤ ∫_{T_{n}}^{t} ∥L(t, s)(F (s, u_{1}(s)) − F (s, u_{2}(s)))∥_{1} ds \\ & ≤ C(R, τ) ∫_{T_{n}}^{t}(t − s)^{−1/2}∥F (s, u_{1}(s)) − F (s, u_{2}(s))∥ ds. \end{aligned}

In view of (1.23) and (3.5) we have here (note that u_{i}(t) ∈ B_{R+1})

∥F (s, u_{1}) − F (s, u_{2})∥ ≤ C(R) ∑_{i=1}^{2} ∥u_{i} − ū∥_{1}∥u_{1} − u_{2}∥_{1} ≤ C(R)ρ∥u_{1} − u_{2}∥_{1},

so that

∥e(t)∥_{1} ≤ C(R, τ)ρ ∫_{T_{n}}^{t}(t − s)^{−1/2}∥u_{1}(s) − u_{2}(s)∥_{1} ds.

Here, by (3.4) and Lemma 3.1,

\begin{aligned}∥u_{1}(s) − u_{2}(s)∥_{1} & = ∥L(s, T_{n})(v_{1} − v_{2}) + w_{1}(s) − w_{2}(s)∥_{1} \\ & ≤ C(R, τ)(s − T_{n})^{−(1−l)/2}∥v_{1} − v_{2}∥_{l} + ∥e(s)∥_{1}. \end{aligned}

Hence, since

∫_{T_{n}}^{t}(t − s)^{−1/2}(s − T_{n})^{−(1−l)/2} ds ≤ Cτ^{l/2},

see (1.29), we get

∥e(t)∥_{1} ≤ C(R, τ)ρ∥v_{1} − v_{2}∥_{l} + C(R, τ)ρ ∫_{T_{n}}^{t}(t − s)^{−1/2}∥e(s)∥_{1} ds,

and by Gronwall’s lemma 1.4

∥e(T_{n+1})∥_{1} ≤ C(R, τ)ρ∥v_{1} − v_{2}∥_{l},

which is (3.6) and part (a) is proved. If ū is a solution of (1.6), then a similar Gronwall argument gives

∥u_{i}(t) − ū(t)∥_{1} ≤ C(R, τ)∥v_{i} − ū(T_{n})∥_{1}, t ∈ [T_{n}, T_{n+1}].

Combining this with the result of part (a) we obtain the second statement of the lemma.

We have now decomposed the operator S as S = L_{n} + N_{n}, where, according to Lemma 3.1, the L_{n} are bounded linear operators on both V and H,

∥L_{n}∥_{L(X)} ≤ C(R, τ), n = 0, …, N − 1, X = V, H,

and where the nonlinear remainder N_{n} : D ⊂ X → V has a Lipschitz constant that can be rendered arbitrarily small (uniformly in n) by restricting to a small neighborhood of ū; see Lemma 3.2.

Note that L_{n} is the linearization of S at ū(T_{n}). In fact, the mapping S : D ⊂ V → V is Fréchet differentiable with L_{n} = S′(ū(T_{n})) ∈ L(V). However, the mapping S : D ⊂ H → H is not differentiable, because D_{n} = B_{V} (ū(T_{n}), ρ) is not a neighborhood of ū(T_{n}) with respect to the topology of H.

3.2. Exponential dichotomies. We need a condition under which we can prove existence of global solutions of (1.6). Let X denote either V or H. We assume that the linear evolution operator L(t, s) has an exponential dichotomy in X on the interval J = [0, T]; see [6]. This means that there are projections P (t) ∈ L(X), t ∈ J, and constants M ≥ 1, β > 0 such that the following conditions are satisfied for s, t ∈ J, t ≥ s.

1. L(t, s)P (s) = P (t)L(t, s).

2. The restriction L(t, s)|_{R(I−P (s))} : R(I − P (s)) → R(I − P (t)) is an isomorphism. We define L(s, t) : R(I − P (t)) → R(I − P (s)) to be its inverse.

3. ∥L(t, s)P (s)∥_{L(X)} ≤ Me^{−β(t−s)}.

4. ∥L(s, t)(I − P (t))∥_{L(X)} ≤ Me^{−β(t−s)}.

The ranges R(P (t)) of the projections P (t) are called stable subspaces; the complementary spaces R(I − P (t)) = N (P (t)) (null spaces of P (t)) are called unstable subspaces. Note that, due to parabolic smoothing, L(t, s) is not invertible. Condition (2) means that L(t, s) can be inverted on the unstable subspaces, so that we can evolve the solution backwards there. Condition (1) implies that the stable/unstable subspaces are invariant. Conditions (3) and (4) mean that solutions decay exponentially: forwards in time on the stable subspaces, backwards in time on the unstable subspaces. The projections are not orthogonal in general.

Example 3.1. (Hyperbolic fixed point.) Let ū ∈ V be a stationary solution of (1.6); Aū = f(ū). Assume, for simplicity, that B = −f′(ū) is self-adjoint so that the linearized operator A + B is self-adjoint, too, and has real spectrum. Let μ_{j}, ϕ_{j} be the eigenvalues and an orthonormal basis of eigenvectors of A + B. The evolution operator L(t, s) = exp(−(t − s)(A + B)) is then given by

e^{−t(A+B)}v = ∑_{j=1}^{∞} e^{−tμ_{j}}(v, ϕ_{j})ϕ_{j}.

The fixed point ū is said to be hyperbolic if the spectrum does not intersect the imaginary axis, i.e., in this case, if all eigenvalues μ_{j} are non-zero. We then have an unbounded sequence

μ_{1} ≤ ⋯ ≤ μ_{q} < 0 < μ_{q+1} ≤ ⋯ < ∞.

It is now easy to show that we have a constant exponential dichotomy with constant stable and unstable projections

P v = ∑_{j=q+1}^{∞}(v, ϕ_{j})ϕ_{j}, I − P = ∑_{j=1}^{q}(v, ϕ_{j})ϕ_{j}.

When we linearize around a time-dependent solution, ū(t), then we cannot use spectral theory in this way. This is why we need the more general concept of exponential dichotomy.

If L(t, s) has an exponential dichotomy on the interval J, then the inhomogeneous boundary value problem, with u_{0} ∈ R(P_{0}), u_{T} ∈ R(I − P_{N}),

\begin{matrix}u & ′ + Au + B(t)u = g(t), t ∈ J, \\ P & (0)u(0) = u_{0}, (I − P (T))u(T) = u_{T} , \end{matrix} (3.7)

is well-posed. We leave this as an exercise; see Exercise 3.1 below. Instead we switch the discussion to the discrete case, where we shall prove an analogous result.

Recall the sequence of operators L_{n} = L(T_{n+1}, T_{n}), n = 0, …, N − 1. Define

L_{mn} = L_{m−1} … L_{n}, m > n; L_{nn} = I

(note L_{n} = L_{n+1,n}). If L(t, s) has an exponential dichotomy on J, then {L_{n}}_{n=0}^{N−1} has a discrete exponential dichotomy; see [6]. This means that there are projections P_{n}, n = 0, …, N, and constants M ≥ 1 and θ ∈ (0, 1) such that the following conditions are satisfied for n = 0, …, N − 1.

1. L_{n}P_{n} = P_{n+1}L_{n}.

2. The restriction L_{n}|_{R(I−P_{n})} : R(I − P_{n}) → R(I − P_{n+1}) is an isomorphism. We define L_{n,n+1} : R(I − P_{n+1}) → R(I − P_{n}) to be its inverse.

3. ∥L_{mn}P_{n}∥_{L(X)} ≤ Mθ^{m−n}, m ≥ n.

4. ∥L_{nm}(I − P_{m})∥_{L(X)} ≤ Mθ^{m−n}, m ≥ n. Note that (1) implies L_{mn}P_{n} = P_{m}L_{mn}, m ≥ n.

Let X_{0}^{S} = R(P_{0}) and X_{N}^{U} = R(I − P_{N}). Consider the following linear boundary value problem with data {y_{n}}_{n=0}^{N+1}, where y_{0} ∈ X_{0}^{S}, y_{N+1} ∈ X_{N}^{U}, and y_{n} ∈ X for n = 1, …, N.

\begin{matrix}x & _{n+1} = L_{n}x_{n} + y_{n+1}, n = 0, …, N − 1, \\ P & _{0}x_{0} = y_{0}, (I − P_{N})x_{N} = y_{N+1}. \end{matrix} (3.8)

A simple calculation gives (Duhamel’s principle)

x_{n} = L_{n0}x_{0} + ∑_{j=1}^{n} L_{nj}y_{j}.

Hence, by condition (1),

P_{n}x_{n} = L_{n0}P_{0}x_{0} + ∑_{j=1}^{n} L_{nj}P_{j}y_{j}.

Similarly, evolving from T_{n} to T_{N},

(I − P_{N})x_{N} = L_{Nn}(I − P_{n})x_{n} + ∑_{j=n+1}^{N} L_{Nj}(I − P_{j})y_{j},

where, by condition (2), L_{Nn}|_{R(I−P_{n})} is invertible, so that

(I − P_{n})x_{n} = L_{nN}(I − P_{N})x_{N} − ∑_{j=n+1}^{N} L_{nj}(I − P_{j})y_{j}.

Hence, using also the boundary conditions, we get

\begin{aligned}x_{n} & = P_{n}x_{n} + (I − P_{n})x_{n} \\ & = L_{n0}P_{0}y_{0} + ∑_{j=1}^{n} L_{nj}P_{j}y_{j} \\ & + L_{nN}(I − P_{N})y_{N+1} − ∑_{j=n+1}^{N} L_{nj}(I − P_{j})y_{j}. \end{aligned} (3.9)

Using (3) and (4) we get

∥x_{n}∥_{X} ≤ M(θ^{n} + ∑_{j=1}^{n} θ^{n−j} + θ^{N−n} + ∑_{j=n+1}^{N} θ^{j−n}) max_{j} ∥y_{j}∥_{X}

≤ M(1 + ∑_{j=n+1}^{N} θ^{j−n} + θ^{n} + ∑_{j=1}^{n} θ^{n−j}) max_{j} ∥y_{j}∥_{X}

≤ M( ∑_{j=n}^{∞} θ^{j−n} + ∑_{j=−∞}^{n} θ^{n−j}) max_{j} ∥y_{j}∥_{X} = \frac{2M}{1 − θ} max_{j} ∥y_{j}∥_{X}.

It follows that (3.8) is uniquely solvable and the solution satisfies

max_{0≤n≤N} ∥x_{n}∥_{X} ≤ \frac{2M}{1 − θ} max_{0≤n≤N+1} ∥y_{n}∥_{X}. (3.10)

We formulate this as a lemma.

Lemma 3.3. Assume that {L_{n}}_{n=0}^{N−1} has a discrete exponential dichotomy in the Banach space X with projections P_{n} and constants M ≥ 1 and θ ∈ (0, 1). Set X = X^{N+1}, Y = X_{0}^{S} × X^{N} × X_{N}^{U}, and define a linear operator A : X → Y by A : (x_{0}, …, x_{N}) ↦ (y_{0}, …, y_{N+1}), where

\begin{matrix}y & _{n+1} = x_{n+1} − L_{n}x_{n}, n = 0, …, N − 1, \\ y & _{0} = P_{0}x_{0}, y_{N+1} = (I − P_{N})x_{N}. \end{matrix}

Then A is invertible and

∥A^{−1}∥_{L(Y,X)} ≤ \frac{2M}{1 − θ} .

For future reference we note that the solution formula (3.9) may be written

x_{n} = G_{n0}y_{0} − G_{nN}y_{N+1} + ∑_{j=1}^{N} G_{nj}y_{j}, (3.11)

by means of the Green’s operator

G_{nj} = \begin{cases}L_{nj}P_{j}, & 0 ≤ j ≤ n, \\ −L_{nj}(I − P_{j}), & n < j ≤ N. \end{cases} (3.12)

Exercise 3.1. Show that the solution of (3.7) may be written

u(t) = G(t, 0)u_{0} − G(t, T)u_{T} + ∫_{0}^{T} G(t, s)g(s) ds, (3.13)

where

G(t, s) = \begin{cases}L(t, s)P (s), & 0 ≤ s ≤ t, \\ −L(t, s)(I − P (s)), & t < s ≤ T. \end{cases} (3.14)

Use this to derive an estimate of ∥u(t)∥_{X} similar to (3.10).

3.3. Shadowing. Let ū ∈ C([0, T], V) be a solution of (1.6) with ∥ū(t)∥_{1} ≤ R, t ∈ [0, T], for some T and R. Assume that we have a numerical solution u_{h}(t) = S_{h}(t, u_{h,0}), t ∈ [0, T], which belongs to a neighborhood of ū in the sense that, for some ρ_{0} ≤ 1,

∥u_{h}(t) − ū(t)∥_{1} ≤ ρ_{0}, t ∈ [0, T]. (3.15)

With T_{n} = nτ chosen as in §3.1 we form the sequence u_{h,n} = S_{h}(nτ, u_{h,0}). This means that u_{h,n} is an orbit of the nonlinear mapping S_{h} = S_{h}(τ, ·), i.e., u_{h,n+1} = S_{h}(u_{h,n}). Comparing with the mapping S = S(τ, ·) we get, in view of Theorem 2.3,

\begin{aligned}∥u_{h,n+1} − S(u_{h,n})∥_{l} & = ∥S_{h}(τ, u_{h,n}) − S(τ, u_{h,n})∥_{l} \\ & ≤ C(R + 1)h^{2−l} =: δ, n = 0, …, N − 1, l = 0, 1. \end{aligned} (3.16)

(A similar result holds for a completely discrete sequence via Theorem 2.6.) With the terminology used in shadowing theory this means that the sequence {u_{h,n}}_{n=0}^{N} is a δ-pseudo-orbit of S with δ = C(R+1)h^{2−l}. Our aim is to find an orbit {u_{n}}_{n=0}^{N} of S, i.e., u_{n+1} = S(u_{n}), such that, for some σ,

max_{0≤n≤N} ∥u_{n} − u_{h,n}∥_{l} ≤ σδ =: ϵ, l = 0, 1.

This means that {u_{h,n}}_{n=0}^{N} is ϵ-shadowed by the true orbit {u_{n}}_{n=0}^{N}.

In order to achieve this we shall prove the following theorem which is set in a more general context.

Theorem 3.4. (Shadowing theorem.) Consider a nonlinear operator S : D ⊂ X → X, where X is a Banach space and D a nonempty subset of X. Assume that S can be decomposed in the form

S = L_{n} + N_{n} (3.17)

where {L_{n}}_{n=0}^{N−1} has a discrete exponential dichotomy on X with projections P_{n} and constants M ≥ 1, θ ∈ (0, 1), and where the operators N_{n} : D_{n} ⊂ X → X are Lipschitz continuous with Lipschitz constants that satisfy

Lip(N_{n}) ≤ \frac{1 − θ}{4M} , n = 0, …, N − 1. (3.18)

Set

σ = \frac{4M}{1 − θ} . (3.19)

(a) Let {x̃_{n}}_{n=0}^{N} be a sequence with x̃_{n} ∈ D_{n}, n = 0, …, N −1. If {x_{n}}_{n=0}^{N} satisfies x_{n} ∈ D_{n}, n = 0, …, N − 1, and

\begin{matrix}x & _{n+1} = S(x_{n}), n = 0, …, N − 1, \\ P & _{0}x_{0} = P_{0}x̃_{0}, (I − P_{N})x_{N} = (I − P_{N})x̃_{N}, \end{matrix} (3.20)

then

max_{0≤n≤N} ∥x_{n} − x̃_{n}∥_{X} ≤ σ max_{0≤n≤N−1} ∥x̃_{n+1} − S(x̃_{n})∥_{X}. (3.21)

(b) Assume that D_{n} contains a closed ball B_{X}(z_{n}, ρ) for n = 0, …, N − 1, where {z_{n}}_{n=0}^{N} is a sequence that satisfies

∥z_{n+1} − S(z_{n})∥_{X} ≤ ρ/σ, n = 0, …, N − 1. (3.22)

Set ρ_{0} = ρ/(σM). Then, for any data x̃_{0} ∈ B_{X}(z_{0}, ρ_{0}), x̃_{N} ∈ B_{X}(z_{N}, ρ_{0}), the boundary value problem (3.20) has unique solution {x_{n}}_{n=0}^{N} with x_{n} ∈ B_{X}(z_{n}, ρ).

The proof is based on the following simple lemmas.

Lemma 3.5. Let X and Y be Banach spaces and let A : X → Y be a linear bijection with bounded inverse A^{−1}. Assume that the mapping N : D ⊂ X → Y, defined in a nonempty subset D of X, is Lipschitz continuous with

α := ∥A^{−1}∥_{L(Y,X)} Lip(N) < 1.

Define F = A + N and set σ := ∥A^{−1}∥_{L(Y,X)}/(1 − α). Then

(3.23) ∥x_{1} − x_{2}∥_{X} ≤ σ∥F(x_{1}) − F(x_{2})∥_{Y}, for all x_{1}, x_{2} ∈ D.

Proof. The bound (3.23) follows readily from the identity

x_{1} − x_{2} = A^{−1}(F(x_{1}) − F(x_{2})) − A^{−1}(N(x_{1}) − N(x_{2})).

Lemma 3.6. In addition to the hypotheses of Lemma 3.5, assume that the domain D of N contains a closed ball B_{X}(z, ρ). Then, for each y ∈ Y with

∥y − F(z)∥_{Y} ≤ ρ/σ, (3.24)

the equation F(x) = y has a unique solution x ∈ B_{X}(z, ρ).

Proof. This is a simple consequence of the contraction mapping theorem. In fact, if we define T(x) = A^{−1}(y − N(x)), then Lip(T) = α < 1. Moreover, the identity

T(x) − z = A^{−1}(y − F(z)) − A^{−1}(N(x) − N(z))

and (3.24) imply that T maps B_{X}(z, ρ) into itself.

Proof of Theorem 3.4. Set X_{0}^{S} = R(P_{0}), X_{N}^{U} = R(I − P_{N}), X = X^{N+1}, Y = X_{0}^{S} × X^{N} × X_{N}^{U}, D = D_{0} × ⋯ × D_{N−1} × X, and define y ∈ Y, A : X → Y, N : D ⊂ X → Y by

\begin{aligned}y & = [\begin{matrix}y_{0} \\ y_{n+1} \\ y_{N+1} \end{matrix}] := [\begin{matrix}P_{0}x̃_{0} \\ 0 \\ (I − P_{N})x̃_{N} \end{matrix} ] , \\ Ax & := [\begin{matrix}P_{0}x_{0} \\ x_{n+1} − L_{n}x_{n} \\ (I − P_{N})x_{N} \end{matrix}] , N(x) := [\begin{matrix}0 \\ −N_{n}(x_{n}) \\ 0 \end{matrix}] , \end{aligned}

where n = 0, …, N − 1, and set F = A + N. The boundary value problem (3.20) then becomes

y = F(x).

The assumption (3.18) implies that Lip(N) ≤ (1 − θ)/(4M). We may therefore apply Lemma 3.5 with α = 1/2; this yields a value of σ = ∥A^{−1}∥_{L(Y,X)}/(1 − α) = 4M/(1 − θ) that coincides with (3.19). The error bound (3.21) follows from (3.23) since

F(x̃) − F(x) = [\begin{matrix}0 \\ x̃_{n+1} − S(x̃_{n}) \\ 0 \end{matrix}] .

This proves part (a). To prove part (b) we apply Lemma 3.6. Set z = (z_{0}, …, z_{N}). Our assumption implies that B_{X}(z, ρ) ⊂ D. Moreover, with y as above,

y − F(z) = [\begin{matrix}P_{0}(x̃_{0} − z_{0}) \\ −(z_{n+1} − S(z_{n})) \\ (I − P_{N})(x̃_{N} − z_{N}) \end{matrix}] .

Here, by conditions (3) and (4) in the definition of discrete exponential dichotomy, we have ∥P_{0}∥_{L(X)} ≤ M, ∥I − P_{N}∥_{L(X)} ≤ M, so that (3.24) follows by our assumptions.

We now return to the problem of shadowing a numerical solution u_{h}(t) in a neighborhood of ū.

Theorem 3.7. (Numerical shadowing theorem 1.) Let ū ∈ C([0, T], V) be a solution of (1.6) with

max_{t∈[0,T]} ∥ū(t)∥_{1} ≤ R, (3.25)

for some T and R. Assume that the solution operator L(t, s) of the linearized problem (3.3) has an exponential dichotomy in V and H on the interval [0, T].

Then there are numbers ρ_{0} and C such that, for each solution u_{h}(t) = S_{h}(t, u_{h,0}), t ∈ [0, T], of (2.2) with

max_{t∈[0,T]} ∥u_{h}(t) − ū(t)∥_{1} ≤ ρ_{0}, (3.26)

there is a solution u of (1.6) such that

∥u_{h}(t) − u(t)∥_{l} ≤ C(1 + t^{−1/2})h^{2−l}, t ∈ (0, T], l = 0, 1. (3.27)

Proof. Let u_{h,n} = u_{h}(T_{n}), T_{n} = nτ, and S = L_{n} + N_{n} be as above. We apply Theorem 3.4 with x̃_{n} = u_{h,n}. Set D = B_{R+1} = B_{V} (0, R + 1), z_{n} = ū(T_{n}), and

D_{n} = {v ∈ V : ∥v − ū(T_{n})∥_{1} ≤ ρ} = B_{V} (z_{n}, ρ).

Lemma 3.2 (b) shows that N_{n} : D_{n} → X has Lipschitz constant

Lip(N_{n}) ≤ C(R)ρ

for both X = V and X = H. We fix ρ so that (3.18) holds and set ρ_{0} = ρ/(σM).

In order to prove existence we use X = V. Then D_{n} is a closed ball in X. Moreover, z_{n+1} − S(z_{n}) = 0 so that (3.22) holds, and (3.26) implies that u_{h,n} = x̃_{n} ∈ B_{V} (z_{n}, ρ_{0}) ⊂ D_{n}, n = 0, …, N. Theorem 3.4 (b) now gives the existence of an orbit u_{n}, n = 0, …, N, of S, satisfying (3.20), and hence, by part (a), (3.21) and (3.16) give the special case l = 1 of the inequality

max_{0≤n≤N} ∥u_{h,n} − u_{n}∥_{l} ≤ σ max_{0≤n≤N−1} ∥u_{h,n+1} − S(u_{h,n})∥_{l} ≤ Ch^{2−l}, l = 0, 1. (3.28)

Another application of Theorem 3.4 (a), now with X = H, proves the case l = 0 of (3.28). Note that D_{n} does not contain any ball defined using the norm of X = H so part (b) does not apply with X = H.

Using the sequence u_{n} we define u(t) = S(t − T_{n}, u_{n}) for t ∈ [T_{n}, T_{n+1}]. By existence and uniqueness of solutions (Theorem 2.3) this is a solution of (1.6). Error bounds at intermediate times are obtained by combining (3.28) with Theorem 2.3 as follows. For t ∈ [0, T_{1}] we have

∥u_{h}(t) − u(t)∥_{1} ≤ C(R)t^{−1/2}(∥u_{h,0} − P_{h}u_{0}∥_{1} + h) ≤ Ct^{−1/2}h.

For t ∈ [T_{n+1}, T_{n+2}], n ≥ 0, we have

∥u_{h}(t) − u(t)∥_{1} ≤ C(R)(∥u_{h,n} − u_{n}∥_{1} + h) ≤ Ch.

This proves the special case l = 1 of (3.27). The case l = 0 is obtained similarly.

Exercise 3.2. State and prove a similar theorem without the assumption that ū is a solution. (Hint: (3.22) means that z_{n} needs only be a pseudo-orbit of S.)

Exercise 3.3. State and prove a similar theorem for the backward Euler discretization in §2.2.

Note that the shadowing property is a property of the evolution operator S; properties of S_{h} (or S_{hk}) enter only through (3.16).

Shadowing in a neighborhood of a hyperbolic fixed point (Example 3.1) was proved in [9] and [10] by the above technique. There we also proved the converse: every exact solution near the fixed point is shadowed by a numerical solution. The proof involves showing that the exponential dichotomy of the evolution operator L carries over to its numerical counterpart and then letting S and S_{h} change roles.

This work was inspired by the paper [1], where similar results were proved in the context of ordinary differential equations.

In [8] we obtained shadowing in a neighborhood of the whole attractor of S. Here we assume that the attractor has the Morse-Smale property, which means, among other things, that S has a finite number of fixed points, all of them are hyperbolic, and the attractor consists of the union of their unstable manifolds. After reduction to finite dimension by means of an inertial manifold, we apply a more complicated, global shadowing theorem of Pilyugin. Again properties of S_{h} (or S_{hk}) enter mainly through (3.16).

4. A posteriori error estimates

In this section we express the error in terms of an appropriately defined residual of the computed solution and the solution of a linearized adjoint problem. The residual is estimated in terms of a posteriori computable residuals. The presentation here is inspired by [4].

4.1. The error equation. We will derive an equation for the error by means of a duality argument. We therefore begin by recalling the weak form of (1.1): find u(t) ∈ V such that

\begin{matrix}( & u′, v) + a(u, v) − (f(u), v) = 0, ∀v ∈ V, t > 0, \\ u(0) = u_{0}, \end{matrix} (4.1)

where a(u, v) = (ϵ∇u, ∇v) = (−ϵΔu, v) = (Au, v) is the bilinear form associated with A. To allow for time-discretization by means of discontinuous piecewise polynomials in t, we introduce a mesh on the time interval [0, T],

0 = t_{0} < ⋯ < t_{j−1} < t_{j} < ⋯ < t_{N} = T,

with time intervals I_{j} = (t_{j−1}, t_{j}) of lengths k_{j} = t_{j} − t_{j−1}, and we define function spaces

\begin{aligned}V & = C([0, T], V), \\ Ṽ & = \{v : v ∈ C(Ī_{j}, V) ∩ C^{1}(Ī_{j}, V^{∗}), j = 1, …, N\}. \end{aligned} (4.2)

Functions in Ṽ are piecewise smooth with respect to the mesh, but may be discontinuous at the nodes t_{j}.

Let U ∈ Ṽ be arbitrary. (In the end U will be a numerical approximation of u, but this is not essential right now.) We define the local residual of U,

\begin{aligned}r_{t_{j−1},t}(U, v) & = ∫_{t_{j−1}}^{t} ((U′, v) + a(U, v) − (f(U), v)) dσ \\ & + ([U]_{t_{j−1}}, v(t_{j−1}^{+})), v ∈ Ṽ, t ∈ (t_{j−1}, t_{j}], \end{aligned} (4.3)

where [v]_{t_{j}} = v(t_{j}^{+}) − v(t_{j}^{−}) denotes the jump at t_{j}.

Let s, t ∈ [0, T] with s < t. We define a punctured interval and two sets of nodes: if s ∈ [t_{m}, t_{m+1}), t ∈ (t_{n−1}, t_{n}], then

\begin{aligned}J_{s,t} & = (s, t_{m+1}) ∪ I_{m+2} ∪ … I_{n−1} ∪ (t_{n−1}, t), \\ J_{s,t}^{−} & = \{t_{j} : s ≤ t_{j} < t\}, J_{s,t}^{+} = \{t_{j} : s < t_{j} ≤ t\}. \end{aligned}

We can now define a global residual

\begin{aligned}R_{s,t}(U, v) & = ∫_{J_{s,t}} ((U′, v) + a(U, v) − (f(U), v)) dσ \\ & + ∑_{J_{s,t}^{−}}([U]_{t_{j}}, v(t_{j}^{+})), v ∈ Ṽ. \end{aligned} (4.4)

Any solution u of (4.1) satisfies

R_{s,t}(u, v) = 0, ∀v ∈ Ṽ. (4.5)

Hence

R_{s,t}(U, v) = R_{s,t}(U, v) − R_{s,t}(u, v)

= ∫_{J_{s,t}} ((U′ − u′, v) + a(U − u, v) − (f(U) − f(u), v)) dσ

+ ∑_{J_{s,t}^{−}}([U − u]_{t_{j}}, v(t_{j}^{+})), ∀v ∈ Ṽ.

We linearize around U. Thus, let

\begin{aligned}e & := U − u, \\ B(t) & := −f′(U(t)) ∈ L(V, H), \\ η = η(U, u) & := f(U) − f(u) − f′(U)(U − u); \end{aligned} (4.6)

cf. (3.2) with ū replaced by U. Then

\begin{aligned}R_{s,t}(U, v) & = ∫_{J_{s,t}} ((e′, v) + a(e, v) + (Be, v) − (η, v)) dσ \\ & + ∑_{J_{s,t}^{−}}([e]_{t_{j}}, v(t_{j}^{+})), ∀v ∈ Ṽ. \end{aligned} (4.7)

We define the bilinear form

\begin{aligned}B_{s,t}(w, v) & = ∫_{J_{s,t}} ((w′, v) + a(w, v) + (Bw, v)) dσ \\ & + ∑_{J_{s,t}^{−}}([w]_{t_{j}}, v(t_{j}^{+})) + (w(s^{−}), v(s^{+})) \\ & = ∫_{J_{s,t}} ((w, −v′) + a(w, v) + (w, B^{∗}v)) dσ \\ & + ∑_{J_{s,t}^{+}}(w(t_{j}^{−}), −[v]_{t_{j}}) + (w(t^{−}), v(t^{+})), w, v ∈ Ṽ, \end{aligned} (4.8)

where the second variant is obtained by integration by parts. When restricted to smooth test functions v ∈ C_{0}^{∞} the relevant terms in (4.8) coincide with the definition of the derivative w′ in the sense of distributions. The point is that we extend this definition to piecewise smooth test functions v ∈ Ṽ.

Now (4.7) becomes

e ∈ Ṽ; B_{s,t}(e, v) = (e(s^{−}), v(s^{+})) + R_{s,t}(U, v) + ∫_{s}^{t}(η, v) dσ, ∀v ∈ Ṽ. (4.9)

We argue by duality. We therefore introduce the linearized adjoint problem with data ψ ∈ H:

−z′(s) + Az(s) + B^{∗}(s)z(s) = 0, s < t; z(t) = ψ. (4.10)

This is the adjoint of the linearized problem (3.3) (but linearized around U(t) instead of ū(t)). Using the second form of (4.8) we write (4.10) as:

z ∈ V; B_{s,t}(w, z) = (w(t^{−}), ψ), ∀w ∈ Ṽ. (4.11)

Let K(s, t) be the solution operator of (4.10), so that its solution is z(s) = K(s, t)ψ. Global existence on [0, T] is obtained in the same way as for the evolution operator L(t, s) in §3.1. In fact, the solution u(t) = L(t, s)ϕ of (3.3) with data ϕ ∈ H satisfies

u ∈ V; B_{s,t}(u, v) = (ϕ, v(s^{+})), ∀v ∈ Ṽ. (4.12)

With v = z in (4.12) and w = u in (4.11) we get

(L(t, s)ϕ, ψ) = B(u, z) = (ϕ, K(s, t)ψ).

We conclude that K(s, t) = L(t, s)^{∗}.

We now return to the equation (4.9) for e. With v = z in (4.9) and w = e in (4.11) we get

(e(t^{−}), ψ) = (e(s^{−}), z(s^{+})) + R_{s,t}(U, z) + ∫_{s}^{t}(η, z) dσ, ∀ψ ∈ H, (4.13)

or, with z(s) = K(s, t)ψ,

\begin{aligned}(e(t^{−}), ψ) & = (e(s^{−}), K(s, t)ψ) \\ & + ∫_{J_{s,t}} ((U′ − f(U), K(σ, t)ψ) + a(U, K(σ, t)ψ)) dσ \\ & + ∑_{J_{s,t}^{−}}([U]_{t_{j}}, K(t_{j}, t)ψ) + ∫_{s}^{t}(η, K(σ, t)ψ) dσ, ∀ψ ∈ H. \end{aligned} (4.14)

Formally, since K(s, t) = L(t, s)^{∗}, we may write

\begin{aligned}e(t^{−}) & = L(t, s)e(s^{−}) + ∫_{s}^{t} L(t, σ)(U′ + AU − f(U)) dσ \\ & + ∫_{s}^{t} L(t, σ)η dσ. \end{aligned} (4.15)

This is the variation of constants formula (Duhamel’s principle) for the equation

e′ + Ae + B(t)e = U′ + AU − f(U) + η.

The first term in the right side of (4.15) describes the evolution of the initial error e(s^{−}); the second term is the contribution to e(t^{−}) from the residual U′+AU −f(U). The third term is the remainder from linearization. However, the derivatives U′(σ) and AU(σ) = −ϵΔU(σ), where U(σ) is only in V, have to be understood in the sense of distributions. This is why we argue by duality and use the weak formulation (4.13).

We finish this subsection by stating an estimate of the remainder η defined in (4.6). It is an immediate consequence of (1.23).

Lemma 4.1. Assume that ∥U(t)∥_{1} ≤ R, ∥u(t)∥_{1} ≤ R, for t ∈ [0, T]. Then, for s, t ∈ [0, T], l = 0, 1,

| ∫_{s}^{t}(η, v) dσ| ≤ C(R)∥e∥_{L_{∞}([s,t],H^{1})}∥e∥_{L_{∞}([s,t],H^{l})}∥v∥_{L_{1}([s,t],H^{1−l})}.

4.2. Local estimates of the residual. We now assume that U is defined by the discontinuous Galerkin finite element method; see [12]. We define a space of functions that are discontinuous, piecewise constant in time, and piecewise linear in space:

Ṽ_{h} = {v ∈ Ṽ : v|_{I_{j}}(t) = v(t_{j}^{−}) ∈ V_{h}, j = 1, …, N}.

Then U ∈ Ṽ_{h} is defined by

∫_{I_{j}} ((U′, v) + a(U, v) − (f(U), v)) dσ + ([U]_{t_{j−1}}, v(t_{j−1}^{+})) = 0, ∀v ∈ Ṽ_{h}, (4.16)

for j = 1, …, N. An initial value U(t_{0}^{−}) = u_{h,0} must be provided. Here, of course, U′ = 0, since U is constant on I_{j}. Setting U_{j} = U(t_{j}^{−}), v(t) = χ ∈ V_{h} on I_{j}, we see that (4.16) becomes

k_{j}a(U_{j}, χ) − k_{j}(f(U_{j}), χ) + (U_{j} − U_{j−1}, χ) = 0, ∀χ ∈ V_{h}, (4.17)

or

A_{h}U_{j} − P_{h}f(U_{j}) + k_{j}^{−1}(U_{j} − U_{j−1}) = 0,

which is the backward Euler method (2.18).

The equation (4.16) means that

r_{t_{j−1},t_{j}}(U, v) = 0, ∀v ∈ Ṽ_{h}. (4.18)

As a consequence of this “orthogonality” we have the following a posteriori estimate of the local residual.

Lemma 4.2. Let U ∈ Ṽ_{h} be a solution of (4.16) and assume that v ∈ V is sufficiently smooth. Then, for l = 0, 1, m = 1, 2, Ĩ_{j} = (t_{j−1}, t) ⊂ I_{j},

\begin{aligned}|r_{t_{j−1},t}(U, v)| & ≤ C∥kR_{t}∥_{L_{∞}(Ĩ_{j},H^{l})} min (∥v′∥_{L_{1}(Ĩ_{j},H^{−l})}, ∥v∥_{L_{∞}(Ĩ_{j},H^{−l})}) \\ & + C∥h^{m}R_{x}∥_{L_{∞}(Ĩ_{j},H)}∥∫Ĩ_{j} v dσ∥_{m}, \end{aligned}

where the functions h = h(x), k = k(t), R_{t} = R_{t}(x, t), R_{x} = R_{x}(x, t), are defined piecewise by

\begin{aligned}h|_{K} & = h_{K} = diam(K), k|_{I_{j}} = k_{j}, \\ R_{t}|_{K×I_{j}} & = k_{j}^{−1}[U]_{t_{j−1}} = k_{j}^{−1}(U_{j} − U_{j−1}), \\ R_{x}|_{K×I_{j}} & = |k_{j}^{−1}[U]_{t_{j−1}} − f(U_{j})| + |K|^{−1/2}h_{K}^{−1/2}∥ϵ[∂U_{j}/∂n]∥_{L_{2}(∂K)}. \end{aligned} (4.19)

Here [∂U_{j}/∂n] is the jump across ∂K in the exterior normal derivative ∂U_{j}/∂n.

Proof. Since U is constant on I_{j} we have by (4.17) and (4.3), for any χ ∈ V_{h},

r_{t_{j−1},t}(U, v) = a(U, V (t)) + (U′ − f(U), V (t)) + ([U]_{t_{j−1}}, v(t_{j−1}^{+}))
= a(U_{j}, V (t)) + (k_{j}^{−1}[U]_{t_{j−1}} − f(U_{j}), V (t))
+ ([U]_{t_{j−1}}, v(t_{j−1}^{+}) − k_{j}^{−1}V (t))
= (a(U_{j}, V (t) − χ) + (k_{j}^{−1}[U]_{t_{j−1}} − f(U_{j}), V (t) − χ))
+ ([U]_{t_{j−1}}, v(t_{j−1}^{+}) − k_{j}^{−1}V (t)) ≡ I + II,

where V (t) = ∫_{t_{j−1}}^{t} v dσ = ∫Ĩ_{j} v dσ. Here, by the definition of R_{t} and a simple estimation,

\begin{aligned}|II| & ≤ ∥k_{j}R_{t}∥_{l}∥v(t_{j−1}^{+}) − k_{j}^{−1}V (t)∥_{−l} \\ & ≤ ∥k_{j}R_{t}∥_{l} min (∥v′∥_{L_{1}(Ĩ_{j},H^{−l})}, 2∥v∥_{L_{∞}(Ĩ_{j},H^{−l})}). \end{aligned}

We now estimate I. We write g = −f(U_{j}) + k_{j}^{−1}[U]_{t_{j−1}} for short. Elementwise integration by parts gives, since ΔU_{j} = 0 on K,

a(U_{j}, w) + (g, w) = ∑_{K} ((ϵ∇U_{j}, ∇w)_{K} + (g, w)_{K})

= ∑_{K} ((−ϵΔU_{j} + g, w)_{K} + (ϵ∂U_{j}/∂n, w)_{∂K})

= ∑_{K} ((g, w)_{K} − \frac{1}{2}(ϵ[∂U_{j}/∂n], w)_{∂K\∂Ω}), ∀w ∈ V.

With w = V (t) − χ we get

\begin{aligned}|I| & ≤ ∑_{K} (∥g∥_{L_{2}(K)}∥V (t) − χ∥_{L_{2}(K)} \\ & + \frac{1}{2}∥ϵ[∂U_{j}/∂n]∥_{L_{2}(∂K\∂Ω)}∥V (t) − χ∥_{L_{2}(∂K)}). \end{aligned}

We choose χ = ΠV (t), where Π : V → V_{h} is an interpolation operator such that

∥D^{l}(w − Πw)∥_{L_{2}(K)} ≤ Ch_{K}^{m−l}∥D^{m}w∥_{L_{2}(S_{K})}, l = 0, 1, m = 1, 2, (4.20)

where S_{K} is the union of all simplices adjacent to K; see [11]. Using also the (scaled) trace inequality

∥w∥_{L_{2}(∂K)} ≤ C(h_{K}^{−1/2}∥w∥_{L_{2}(K)} + h_{K}^{1/2}∥Dw∥_{L_{2}(K)}), (4.21)

we obtain

∥V (t) − ΠV (t)∥_{L_{2}(∂K)} ≤ Ch_{K}^{m−1/2}∥D^{m}V (t)∥_{L_{2}(S_{K})}, m = 1, 2.

Hence

|I| ≤ C ∑_{K} h_{K}^{m}(∥g∥_{L_{2}(K)}
+ h_{K}^{−1/2}∥ϵ[∂U_{j}/∂n]∥_{L_{2}(∂K\∂Ω)})∥D^{m}V (t)∥_{L_{2}(S_{K})}
≤ C∥h^{m}R_{x}∥_{L_{2}} ∥D^{m}V (t)∥_{L_{2}} ≤ C∥h^{m}R_{x}∥_{L_{2}} ∥∫Ĩ_{j} v dσ∥_{m},

where the constant depends only on the constants in (4.20) and (4.21).

Note that the quantities R_{t} and R_{x}, in the previous theorem are computable a posteriori from the numerical solution U.

We now show that the local residuals tend to zero with h and k by using local a priori estimates from §2. Essentially, the idea is to compare U with an auxiliary local solution ũ(t) = S(t − T_{n}, U(T_{n})) (we use the notation T_{n} = nτ from §3). In order to present the idea with a minimum of technical details we make a very strong assumption about the mesh. For a more general result we refer to [4].

Theorem 4.3. Assume that the family of spatial meshes is quasi-uniform and that the the time step is constant with h^{2}k^{−1} ≤ C, where now h = max_{K} h_{K}, k = k_{j}. Let U ∈ Ṽ_{h} be a solution of (4.16) with ∥U(t)∥_{1} ≤ R for t ∈ [0, T]. Let τ be as in Theorems 2.3 and 2.6 and assume that t_{j−1} ≥ τ. Then, for l = 0, 1, we have

∥kR_{t}∥_{L_{∞}(I_{j},H^{l})} + ∥h^{2−l}R_{x}∥_{L_{∞}(I_{j},H)} ≤ C(R)(h^{2−l} + k).

Proof. Since t_{j−1} ≥ τ we have I_{j} ⊂ [T_{n+1}, T_{n+2}] for some n ≥ 0 (we use the notation T_{n} = nτ from §3). Let ũ(t) = S(t − T_{n}, U(T_{n})) be the local exact solution starting with ũ(T_{n}) = U(T_{n}). Theorem 2.6 gives for the error ẽ = U − ũ the bound

∥ẽ∥_{L_{∞}((T_{n+1},T_{n+2}),H^{l})} ≤ C(R)(h^{2−l} + k).

Hence

\begin{aligned}∥kR_{t}∥_{L_{∞}(I_{j},H^{l})} & = ∥[U]_{t_{j−1}}∥_{l} = ∥[ẽ]_{t_{j−1}}∥_{l} ≤ 2∥ẽ∥_{L_{∞}((T_{n+1},T_{n+2}),H^{l})} \\ & ≤ C(R)(h^{2−l} + k). \end{aligned}

Similarly, for the first term in R_{x},

∥h^{2−l}k_{j}^{−1}[U]_{t_{j−1}}∥_{L_{∞}(I_{j},H)} ≤ C(R)h^{2−l}k^{−1}(h^{2−l} + k) ≤ C(R)(h^{2−l} + k).

Also, by Lemma 1.1,

∥h^{2−l}f(U)∥_{L_{∞}(I_{j},H)} ≤ C(R)h^{2−l}.

Finally, for the last contribution to R_{x}, we use (4.21) and Theorem 1.3,

h^{2−l}( ∑_{K} h^{−1}∥ϵ[∂U_{j}/∂n]∥_{L_{2}(∂K)}^{2})^{1/2} = h^{2−l}( ∑_{K} h^{−1}∥ϵ[∂ẽ/∂n]∥_{L_{2}(∂K)}^{2})^{1/2}

≤ Ch^{2−l}( ∑_{K} h^{−1}∥Dẽ∥_{L_{2}(∂K)}^{2})^{1/2}

≤ Ch^{2−l}( ∑_{K} h^{−1}(h^{−1/2}∥Dẽ∥_{L_{2}(K)} + h^{1/2}∥D^{2}ẽ∥_{L_{2}(K)})^{2})^{1/2}

≤ Ch^{2−l}(h^{−1}∥Dẽ∥_{L_{∞}((T_{n+1},T_{n+2}),H)} + ∥D^{2}ẽ∥_{L_{∞}((T_{n+1},T_{n+2}),H)})

≤ C(h^{1−l}∥ẽ∥_{L_{∞}((T_{n+1},T_{n+2}),H^{1})} + h^{2−l}∥D^{2}ũ∥_{L_{∞}((T_{n+1},T_{n+2}),H)})

≤ C(R)(h^{2−l} + k).

4.3. A global error estimate. We now combine the results of §4.1 and §4.2 to obtain a global a posteriori error estimate. We estimate both the H and V norms of the error (l = 0, 1). To achieve this we use a duality argument:

∥e∥_{l} ≤ C∥A^{l/2}e∥ = C sup_{∥φ∥=1}(A^{l/2}e, φ) = C sup_{∥φ∥=1}(e, A^{l/2}φ).

We thus take ψ = A^{l/2}φ, ∥φ∥ = 1, in (4.13) and (4.10), so that z(s) = K(s, t)A^{l/2}φ. Using also Lemma 4.2 with m = 2 − l and Lemma 4.1 we get the following.

Theorem 4.4. Let U and u be solutions of (4.16) and (4.1), respectively. Assume that ∥U(t)∥_{1} ≤ R, ∥u(t)∥_{1} ≤ R, for t ∈ [0, T]. Then, for t ∈ [0, T], l = 0, 1, with n such that t ∈ [t_{n}, t_{n+1}),

∥e(t^{−})∥_{l} ≤ CS_{l}^{0}(0, t)∥e(0^{−})∥_{l}
+ CS_{l}^{1}(0, t)∥kR_{t}∥_{L_{∞}((0,t),H^{l})}
+ CS_{l}^{2}(0, t)∥h^{2−l}R_{x}∥_{L_{∞}((0,t),H)}
+ C(R)S_{l}^{3}(0, t)∥e∥_{L_{∞}((0,t),H^{1})} ∥e∥_{L_{∞}((0,t),H^{l})},

where the so called stability factors are defined by

S_{l}^{0}(0, t) = sup_{∥φ∥=1} max_{τ∈[0,t]} ∥A^{−l/2}K(0, τ)A^{l/2}φ∥,

S_{l}^{1}(0, t) = sup_{∥φ∥=1} ( ∫_{0}^{t_{n−1}} ∥A^{−l/2}D_{s}K(s, t)A^{l/2}φ∥ ds

+ max_{s∈[t_{n−1},t]} ∥A^{−l/2}K(s, t)A^{l/2}φ∥),

S_{l}^{2}(0, t) = sup_{∥φ∥=1} ( ∫_{0}^{t_{n−1}} ∥A^{(2−l)/2}K(s, t)A^{l/2}φ∥ ds

+ ∥∫_{t_{n−1}}^{t} A^{(2−l)/2}K(s, t)A^{l/2}φ ds∥),

S_{l}^{3}(0, t) = sup_{∥φ∥=1} ∫_{0}^{t} ∥A^{(1−l)/2}K(s, t)A^{l/2}φ∥ ds.

Note in this connection that

∥A^{−l/2}K(s, t)A^{l/2}∥_{L(H)} = ∥A^{l/2}L(t, s)A^{−l/2}∥_{L(H)}.

The special treatment of the last interval (t_{n−1}, t) is necessitated by the expected singular behavior (cf. Lemma 3.1)

∥A^{(2−l)/2}K(s, t)A^{l/2}φ∥ ∼ ∥A^{−l/2}D_{s}K(s, t)A^{l/2}φ∥ ∼ (t − s)^{−1}, as s → t.

The stability factors grow fast with t if the linearized problem is unstable.

Of course, the error estimate of the the previous theorem is not “closed” because the square of the error occurs on the right hand side. However, this square remainder can be removed by means of the contraction mapping principle, if we compute with a sufficiently small tolerance. We demonstrate this in a slightly different context in the next section.

5. Shadowing—second approach

We finish these lectures by describing another approach to shadowing based on a posteriori error analysis. This is work in progress together with M. G. Larson, and we only sketch a simplified version.

The error estimate in Theorem 4.4 is based on the stability of the adjoint (4.10) of the linearized initial value problem (3.3). We now use the stability of the boundary value problem (3.7). Let t ∈ (0, T), ψ ∈ H, and let z be a solution of

−z′ + Az + B^{∗}(s)z = 0, s ∈ (0, T), s ≠ t; [z]_{t} = −ψ. (5.1)

Thus z has a jump at s = t. This equation does not determine z uniquely; we must provide appropriate boundary conditions at s = 0, T. But for the time being we simply assume that we have some solution operator, denoted z(s) = G(t, s)^{∗}ψ, of (5.1).

Multiplying (from the left) by a test function w ∈ Ṽ and integrating by parts as in (4.8), taking into account that z is discontinuous at t, we get

0 = ∫_{0}^{T} (w, −z′ + Az + B^{∗}z) dσ

= ∫_{J_{0,T}} ((w′, z) + a(w, z) + (Bw, z)) dσ

+ ∑_{J_{0,T}^{−}}([w]_{t_{j}}, z(t_{j})) + (w(0^{−}), z(0)) − (w(T^{−}), z(T)) + (w(t^{−}), [z]_{t}).

Thus z satisfies

B̃_{0,T} (w, z) = (w(t^{−}), ψ), ∀w ∈ Ṽ, (5.2)

where now, instead of (4.8),

\begin{aligned}B̃_{0,T} (w, v) & = ∫_{J_{0,T}} ((w′, v) + a(w, v) + (Bw, v)) dσ \\ & + ∑_{J_{0,T}^{−}}([w]_{t_{j}}, v(t_{j}^{+})) + (w(0^{−}), v(0^{+})) − (w(T^{−}), v(T^{+})). \end{aligned} (5.3)

Combining this definition with (4.7) yields, instead of (4.9),

\begin{aligned}e ∈ Ṽ; B̃_{0,T} (e, v) & = (e(0^{−}), v(0^{+})) − (e(T^{−}), v(T^{+})) \\ & + R_{0,T} (U, v) + ∫_{0}^{T} (η, v) dσ, ∀v ∈ Ṽ. \end{aligned} (5.4)

Together with (5.2) this gives, instead of (4.13),

\begin{aligned}(e(t^{−}), ψ) & = (e(0^{−}), z(0)) − (e(T^{−}), z(T)) \\ & + R_{0,T} (U, z) + ∫_{0}^{T} (η, z) dσ, ∀ψ ∈ H. \end{aligned} (5.5)

Note in this connection that z(s) may have an additional jump at s = t in which case z ∉ Ṽ. However, we may still use v = z in (5.4) because R_{0,T} (U, v) does not involve v′; see (4.4).

The error representation (5.5) is similar to (3.13). In fact, if the linearized problem has an exponential dichotomy, then we choose z(s) = G(t, s)^{∗}ψ, where G(t, s) is the Green’s operator in (3.14).

With U given we want to solve equation (5.5) for e. We then obtain a shadow solution u = U −e. In order to achieve this we set the boundary data (e(0^{−}), z(0)) = (e(T^{−}), z(T)) = 0. With U fixed we consider η = η(U, u) = η(U, U −e) as a function of e and define the nonlinear operator

N_{0,T} (e, v) = ∫_{0}^{T} (η(U, U − e), v) dσ.

We then write (5.5) as a fixed point equation

(e(t^{−}), ψ) = R_{0,T} (U, z) + N_{0,T} (e, z), ∀ψ ∈ H. (5.6)

As before we take ψ = A^{l/2}φ, ∥φ∥ = 1, in (5.6), so that z(s) = G(t, s)^{∗}A^{l/2}φ. Define the ball

B_{ρ} = {e ∈ L_{∞}((0, T), V) : ∥e∥_{L_{∞}((0,T),H^{1})} ≤ ρ}.

Similarly to Lemma 4.1 it follows that N_{0,T} has a small Lipschitz constant in B_{ρ}: if ψ = A^{1/2}φ, ∥φ∥ = 1, then

|N_{0,T} (e_{1}, z) − N_{0,T} (e_{2}, z)| ≤ S_{1}^{3}C(R)ρ∥e_{1} − e_{2}∥_{L_{∞}((0,T),H^{1})}, e_{1}, e_{2} ∈ B_{ρ}. (5.7)

By Lemma 4.2 we get an error bound similar to that in Theorem 4.4. It is of the form

|R_{0,T} (U, z)| ≤ CS_{1}^{1}∥kR_{t}∥_{L_{∞}((0,T),H^{1})} + CS_{1}^{2}∥hR_{x}∥_{L_{∞}((0,T),H)} =: R(U). (5.8)

Here the stability factors S_{1}^{i} = S_{1}^{i}(0, T), i = 1, 2, 3, are defined as in Theorem 4.4, but using G(t, s)^{∗} instead of K(s, t) = L(t, s)^{∗} and giving special treatment to an interval near s = t instead of the endpoint s = T.

We can now apply the contraction mapping principle to the equation (5.6) in the ball B_{ρ} and obtain the following.

Theorem 5.1. (Numerical shadowing theorem 2.) Assume that

∥U∥_{L_{∞}((0,T),H^{1})} ≤ R.

Let ρ be so small that

S_{1}^{3}C(R)ρ ≤ \frac{1}{2}. (5.9)

Let h and k be so small that

R(U) ≤ \frac{1}{2}ρ. (5.10)

Then U has a shadow u with

∥U − u∥_{L_{∞}((0,T),H^{1})} ≤ ρ. (5.11)

Proof. The assumption (5.9) together with (5.7) guarantees that the right side of (5.6) defines a contraction on B_{ρ}. The assumption (5.10) together with (5.8) guarantees that B_{ρ} is mapped into itself. Hence there is a unique solution e ∈ B_{ρ} of (5.6). The shadow solution is given by u = U − e.

As in Theorem 4.3 it follows from local a priori estimates that the quantity R(U) tends to zero with h and k, so that (5.10) can be achieved if the stability factors occurring in its definition are not too large. This would be the case if the linearized problem has an exponential dichotomy on the interval (0, T).

Note that, in addition to the existence of u, we also get an error bound (5.11) in the V norm. By an additional application of the contraction property, now with respect to H, we get an error bound also in H, cf. the proof of Theorem 3.7.

The radius ρ plays the role of a computational tolerance; see (5.10) and (5.8). Thus, we have shadowing if we compute with a sufficiently small tolerance.

The fixed point iteration that we employ in the proof is a quasi-Newton iteration for the equation

u′ + Au − f(u) = 0,

namely

e′_{n+1} + Ae_{n+1} − f′(U(t))e_{n+1} = −(u′_{n} + Au_{n} − f(u_{n})), u_{n} = U − e_{n}.

However, the derivatives that occur here must be understood in a weak sense. That is why we argue by duality as in (5.6).

References

1. W.-J. Beyn, On the numerical approximation of phase portraits near stationary points, SIAM J. Numer. Anal. 24 (1987), 1095–1113.

2. P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.

3. C. M. Elliott and S. Larsson, Error estimates with smooth and nonsmooth data for a finite element method for the Cahn-Hilliard equation, Math. Comp. 58 (1992), 603–630, S33–S36.

4. D. Estep, M. G. Larson, and R. D. Williams, Estimating the error of numerical solutions of systems of nonlinear reaction-diffusion equations, Mem. Amer. Math. Soc., to appear.

5. L. C. Evans, Partial Differential Equations, Graduate Studies in Mathematics, vol. 19, American Mathematical Society, Providence, RI, 1998.

6. D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, vol. 840, Springer-Verlag, 1981.

7. S. Larsson, Nonsmooth data error estimates with applications to the study of the long-time behavior of finite element solutions of semilinear parabolic problems, preprint 1992–36, Department of Mathematics, Chalmers University of Technology, 1992.

8. S. Larsson and S. Yu. Pilyugin, Numerical shadowing near the global attractor for a semilinear parabolic equation, preprint 1998–21, Department of Mathematics, Chalmers University of Technology and Göteborg University, 1998.

9. S. Larsson and J.-M. Sanz-Serna, The behavior of finite element solutions of semilinear parabolic problems near stationary points, SIAM J. Numer. Anal. 31 (1994), 1000–1018.

10. , A shadowing result with applications to finite element approximation of reactiondiffusion equations, Math. Comp. 68 (1999), 55–72.

11. L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (1990), 483–493.

12. V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Series in Computational Mathematics, vol. 25, Springer-Verlag, Berlin, 1997.

Department of Mathematics, Chalmers University of Technology and Göteborg University, SE–412 96 Göteborg, Sweden

E-mail address: stig@math.chalmers.se, http://www.math.chalmers.se/~stig