THE DISCONTINUOUS GALERKIN METHOD

FOR SEMILINEAR PARABOLIC PROBLEMS

Donald Estep and Stig Larsson

Abstract. We prove a priori error estimates for a space-time finite element method for semilinear parabolic problems. The finite element method has basis functions that are continuous in space and discontinuous in time, and variable spatial meshes and time steps are allowed. The effect of numerical quadrature is emphasized. Résumé. Nous montrons des estimations d’erreur a priori pour une méthode des éléments finis en espace et en temps pour des problèmes paraboliques semi-linéaires. La méthode des éléments finis considerée a des fonctions de base continues en espace et discontinues en temps, et admet des maillages spatiales et des pas de temps variables. L’effet de quadrature numérique est accentué.

1. Introduction. In this paper we consider the numerical solution of the semilinear parabolic equation

(1.1)

u_{t} − Δu = f(x, t, u), in Ω × (0, t^{∗}),
u = 0,  on ∂Ω × (0, t^{∗}),
u(·, 0) = u_{0},  in Ω,

by using a finite element method with basis functions that are continuous in space and discontinuous in time, which we refer to as the discontinuous Galerkin method. Here Ω is a bounded convex polygonal domain in R^{2}, t^{∗} > 0, u_{t} = ∂u/∂t, Δu = ∂^{2}u/∂x_{1}^{2} + ∂^{2}u/∂x_{2}^{2} and f is a smooth function satisfying

(1.2)  |f(x, t, s)| + |∂f(x, t, s)/∂s| ≤ M ∀x ∈ Ω, t > 0, s ∈ R.

The discontinuous Galerkin method for linear parabolic equations, i.e., the case of (1.1) when f = f(x, t), has been analyzed by Eriksson, Johnson and Thomee [3] and Eriksson and Johnson [2]. There are two parts to the analysis in [2]: an a posteriori error analysis, which is used to devise a global error control for an adaptive finite element method, and an a priori error analysis, which guarantees convergence. The a priori error bound is of a specific form related to the a posteriori bound and is proved under conditions which allow variable spatial meshes and time steps.

The purpose of the present work is to discuss some ways of handling the nonlinear term in this context, and to prove a priori error bounds in the style of [2]. In particular, we want to allow the possibility that the spatial mesh will change on each time level.

To avoid some technical difficulties we assume that such refinement (un-refinement) is performed by the addition (removal) of nodes to the current mesh. In fact, it is

Typeset by AMS-TEX

1991 Mathematics Subject Classification. 65M15, 65M60, 65M50.

Key words and phrases. Semilinear parabolic problem, discontinuous Galerkin method, finite element method, adaptive mesh, numerical quadrature, error estimate.

Supported by DARPA and the Swedish National Board for Technical Development (STUF).

Published in RAIRO Modél. Math. Anal. Numér. 26 (1992), 35–54.

known that meshes that are changed in an uncontrolled fashion may yield completely false results, see Dupont [1].

The discontinuous Galerkin method involves integrals of the right-hand side of (1.1) with respect to x and t, which we evaluate by numerical quadrature. Our main result concerns the effect of such numerical integration on the a priori error analysis. Eriksson and Johnson do not consider this in [2], since for linear problems this is less important.

We describe the pertinent results of [2] in Section 2. In Section 3 we formulate our numerical method and state our main result, which is proved in Section 4. We conclude in Section 5 by providing an example showing the implementation of one of our algorithms.

2. The linear case. In this section we briefly review the notation and results of Eriksson and Johnson [2], which we will use in our analysis. For the discretization of (1.1) with respect to x = (x_{1}, x_{2}) we let Σ be the class of all finite element discretizations S = (h, T, S) satisfying the following conditions:

(1) h ∈ C^{1}(Ω) is a positive function satisfying

|∇h(x)| ≤ μ ∀x ∈ Ω;

(2) T = {K} is a partition of Ω into triangles K of diameter h_{K} such that

c_{1}h_{K}^{2} ≤ area(K) ∀K ∈ T,

and associated with the function h through the inequalities

c_{2}h_{K} ≤ h(x) ≤ h_{K} ∀x ∈ Ω, K ∈ T;

(3) S is the set of all functions which are continuous on \overline{Ω}, linear on each K ∈ T

and vanish on ∂Ω. For a given domain Ω the positive constants c_{1}, c_{2} and μ characterize the family Σ completely.

For the discretization with respect to t we introduce a partition 0 = t_{0} < t_{1} < … < t_{n−1} < t_{n} < … of R^{+} into subintervals I_{n} = (t_{n−1}, t_{n}) of lengths k_{n} = t_{n} − t_{n−1}, and we associate with each I_{n} a finite element discretization S_{n} = (h_{n}, T_{n}, S_{n}) ∈ Σ. For q = 0 and 1 we define

V_{n}^{q} = {V : V (x, t) = ∑_{j=0}^{q} t^{j}φ_{j}(x), φ_{j} ∈ S_{n}}.

The discontinuous Galerkin method for (1.1) with f = f(x, t) consists in computing a function U such that U|_{Ω×I_{n}} ∈ V_{n}^{q} and (2.1)

∫_{I_{n}}((U_{t}, X) + (∇U, ∇X)) dt + ([U]_{n−1}, X_{n−1}^{+}) = ∫_{I_{n}}(f(·, t), X(t)) dt ∀X ∈ V_{n}^{q},

for n = 1, 2, …, and U_{0}^{−} = u_{0}, where [V]_{n} = V_{n}^{+} − V_{n}^{−}, V_{n}^{±} = lim_{s→0^{±}} V (t_{n} + s) (V_{n}^{−}

should be considered the “nodal value” of V ∈ V_{n}^{q}). Here and below (·, ·) denotes the usual inner product in L_{2} = L_{2}(Ω) and ∥ · ∥ is the corresponding norm. Hence, we are computing a finite element approximation to u on space-time “slabs” Ω × I_{n}. By virtue of the discontinuity of U in time, we can alter the spatial mesh from one time interval to the next.

It turns out that for q = 0 the scheme (2.1) reduces to the following modification of the backward Euler method:

\frac{1}{k_{n}} (U_{n}^{−} − U_{n−1}^{−}, χ) + (∇U_{n}^{−}, ∇χ) = \frac{1}{k_{n}} ∫_{I_{n}}(f(·, t), χ) dt ∀χ ∈ S_{n}.

For q = 1 in the linear homogeneous case (f ≡ 0) the approximation agrees at the nodal points t_{n} with the subdiagonal third order accurate Padé difference approximation. The approximation is second order accurate in the interiors of the intervals I_{n} (see [3]).

We quote the following a priori error bound:

Theorem 2.1. (Eriksson and Johnson [2]) Let u be the solution of (1.1) with f = f(x, t) and U that of (2.1). Suppose that μ is sufficiently small and assume that for each n one of the following assumptions hold:

S_{n} ⊂ S_{n−1} or (max_{x∈Ω} h_{n}(x))^{2} ≤ γk_{n},

with γ sufficiently small and that k_{n} ≤ ck_{n+1} for all n. Then there exists a constant C depending only on c_{1} and c_{2} such that for q = 0, 1 and n = 1, 2, …,

∥U − u∥_{L_{∞}(I_{n},L_{2})} ≤ CL_{n} max_{1≤j≤n} E_{j}^{(q)}(u),
∥U_{n}^{−} − u(t_{n})∥ ≤ CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u),

with L_{n} = \frac{1}{4} \sqrt{1 + log \frac{t_{n}}{k_{n}}} and

E_{j}^{(p)}(u) = ∥h_{j}^{2}D^{2}u∥_{L_{∞}(I_{j},L_{2})} + min_{s≤p+1} k_{j}^{s}∥u_{t}^{(s)}∥_{L_{∞}(I_{j},L_{2})},

where u_{t}^{(1)} = u_{t}, u_{t}^{(2)} = u_{tt}, u_{t}^{(3)} = Δu_{tt} and D^{m}u = (∑_{|α|=m} |D_{x}^{α}u|^{2})^{1/2} (with the usual multi-index notation).

These error bounds imply that the scheme is of second order in x, of order q + 1 uniformly in t, and of order 2q + 1 at the nodes t_{n}.

3. The semilinear case. In this section we discuss some ways of implementing the scheme (2.1) in the semilinear case, i.e., when f = f(x, t, u). (In the sequel we sometimes write f(u) instead of f(x, t, u(x, t)) for compactness of notation. We also use the notation u_{l} = u(·, t_{l}), f(u)_{l} = f(·, t_{l}, u_{l}) and f(V)_{l} = f(·, t_{l}, V_{l}^{−}) for V ∈ V_{l}^{q}.)

In this case the right-hand side of (2.1) is an integral of f(U)X over Ω × I_{n}, which must be evaluated by numerical quadrature if the algorithm is to be employed in a general way. This is the question that we address here. We desire that the completely discrete schemes should retain the order that the discontinuous Galerkin schemes have on linear problems. In particular, for q = 1, in order to retain the nominal third order accuracy of the discretization of the left-hand side of (2.1), we employ a multistep formula based on interpolation of f(U) at the nodes t_{n}, where we expect third order accuracy for U.

For the integral with respect to x we choose the following quadrature rule. Let S = (h, T, S) ∈ Σ be a finite element discretization and define

Q_{K}(f) = \frac{1}{3}area(K) ∑_{j=1}^{3} f(P_{K,j}) ≈ ∫_{K} f(x) dx ∀K ∈ T,

where P_{K,j} are the vertices of the triangle K. We may then define approximations of the inner product and norm of L_{2} by

(χ, ψ)_{S} = ∑_{K∈T} Q_{K}(χψ), ∥χ∥_{S} = (χ, χ)_{S}^{1/2}.

For the approximation of the integral with respect to t we replace f(U) by an interpolant with respect to t and integrate the resulting polynomials exactly. Two possibilities suggest themselves: one is to use extrapolation of an interpolant computed over previous intervals I_{l}, l < n, and the other is to use an interpolant computed over the current and previous intervals I_{l}, l ≤ n. The former process will yield a set of linear equations for U (semi-implicit method), whereas the latter will produce a nonlinear system (nonlinear implicit method).

We use an interpolant of order p = 0 when q = 0 and of order p = 2 when q = 1. Of course, when p = 2 we cannot use this kind of interpolant on the first one or two intervals, necessitating the construction of special interpolants there.

Thus our scheme is of the form: find a function U such that U|_{Ω×I_{n}} ∈ V_{n}^{q}, U_{0}^{−} = u_{0} and

(3.1) ∫_{I_{n}}((U_{t}, X) + (∇U, ∇X)) dt + ([U]_{n−1}, X_{n−1}^{+})

= ∫_{I_{n}}⟨Πf(U)(t), X(t)⟩ dt ∀X ∈ V_{n}^{q}, 0 < t_{n} ≤ t^{∗},

where (except possibly for n = 1 or 2) the integrand of the right-hand side is given by (3.2)

⟨Πf(U)(t), X(t)⟩ = ⟨Π_{n−i}^{p}f(U)(t), X(t)⟩ = ∑_{l=n−i−p}^{n−i} ϕ_{l}(t)(f(·, t_{l}, U_{l}^{−}), X(t))_{S_{l}∪S_{n}}.

Here Π_{n−i}^{p} is the interpolation operator for polynomials of degree p = 2q, computed with respect to the mesh points {t_{l}}_{l=n−i−p}^{n−i} and ϕ_{l} are the corresponding Lagrange basis functions. Note that i = 0 for the nonlinear implicit scheme and i = 1 for the semi-implicit scheme.

In (3.2) we handle the possibility of variable meshes in the following way: we assume that all spatial meshes that occur are refinements of one common coarse mesh, and that each triangulation T_{n} is obtained from its precursor T_{n−1} by adding some nodes and by removing some other nodes. In this way, the union of the mesh points of two triangulations T_{l} and T_{n} form a triangulation whose finite element space is equal to S_{l} ∪ S_{n}. In this situation we define

h_{l,n}(x) = min{h_{l}(x), h_{n}(x)}.

The discrete L_{2} inner product in (3.2) is thus computed over the totality of all mesh points of T_{l} and T_{n}.

When q = 0 the equation (3.1) becomes

\frac{1}{k_{n}} (U_{n}^{−}−U_{n−1}^{−}, χ)+(∇U_{n}^{−}, ∇χ) = (f(·, t_{n−i}, U_{n−i}^{−}), χ)_{S_{n−i}∪S_{n}} ∀χ ∈ S_{n}, 0 < t_{n} ≤ t^{∗},

which is a variant of the backward Euler method. In Section 5 we give details of the implementation of a scheme with q = 1.

It remains to devise a starting procedure for the case q = 1, p = 2, i.e., to define the right-hand side ∫_{I_{n}}⟨Πf(U), X⟩ dt for 1 ≤ n ≤ 1 + i. For the nonlinear implicit scheme (i = 0) we simply take

⟨Πf(U), X⟩ = ⟨Π_{1}^{1}f(U), X⟩ on I_{1},

i.e., we use a linear interpolant in (3.2). For the semi-implicit method (i = 1) we similarly define

⟨Πf(U), X⟩ = ⟨Π_{1}^{1}f(U), X⟩ on I_{2}.

For n = 1 we use a prediction-correction procedure in order to obtain the correct accuracy. For the predicted value of U_{1}^{−} we define Ũ|_{Ω×I_{1}} ∈ V_{1}^{q} such that Ũ_{0}^{−} = u_{0} and

∫_{I_{1}}((Ũ_{t}, X) + (∇Ũ, ∇X)) dt + ([Ũ]_{0}, X_{0}^{+}) = ∫_{I_{1}}⟨Π_{0}^{0}f(Ũ), X⟩ dt ∀X ∈ V_{1}^{q},

that is,

∫_{I_{1}}((Ũ_{t}, X) + (∇Ũ, ∇X)) dt + ([Ũ]_{0}, X_{0}^{+}) = ∫_{I_{1}}(f(·, 0, u_{0}), X(t))_{S_{0}∪S_{1}} dt ∀X ∈ V_{1}^{q}.

Then we compute U|_{Ω×I_{1}} ∈ V_{1}^{q} such that U_{0}^{−} = u_{0} and

∫_{I_{1}}((U_{t}, X) + (∇U, ∇X)) dt + ([U]_{0}, X_{0}^{+}) = ∫_{I_{1}}⟨Π_{1}^{1}f(Ũ), X⟩ dt ∀X ∈ V_{1}^{q}.

For the scheme described above we have the following a priori error bound. We need all the assumptions of this and the previous sections and, in addition, we assume that, for some c_{3}, c_{4} > 0,

(3.3)  c_{3} ≤ \frac{k_{n}}{k_{n−1}} ≤ c_{4}, 0 < t_{n} ≤ t^{∗}.

It is also convenient to define

(3.4)  m = max\{i + p − 1, 0\} = \{ \begin{matrix}0, & if q = 0, \\ 1 + i, & if q = 1, \end{matrix}

to distinguish between the general case of (3.2) (t_{n} ≥ t_{m+1}) and the special initialization procedure (0 < t_{n} ≤ t_{m}) needed when q = 1.

Theorem 3.1. Let u be the solution of (1.1) and U be that of (3.1). Then there are constants k and C depending on M, t^{∗} and c_{l}, l = 1, …, 4, such that, if k_{n} ≤ k for 0 < t_{n} ≤ t^{∗}, then

∥U − u∥_{L_{∞}(I_{n},L_{2})} ≤ CL_{n} max_{1≤j≤n} E_{j}^{(q)}(u) + C max_{1≤j≤n} F_{j}^{(2q)}(u), 0 < t_{n} ≤ t^{∗},
∥U_{n}^{−} − u(t_{n})∥ ≤ CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u) + C max_{1≤j≤n} F_{j}^{(2q)}(u), 0 < t_{n} ≤ t^{∗}.

For the semi-implicit schemes (i = 1) the condition k_{n} ≤ k is not needed. Here L_{n} and E_{j}^{(p)}(u) are as in Theorem 2.1 and

F_{j}^{(p)}(u) = ∑_{l=j−i−p}^{j−i}(∥h_{l,j}^{2}D_{x}^{2}u_{l}∥ + ∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥)

+ min_{s≤p+1}(k_{j}^{s}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{j},L_{2})}), t_{m+1} ≤ t_{j} ≤ t^{∗},

and for q = 1 we have, in addition, if i = 0,

F_{1}^{(2)}(u) = ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}u_{l}∥ + ∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥)

+ min_{s≤2}{k_{1}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})}},

and, if i = 1,

F_{1}^{(2)}(u) = ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}u_{l}∥ + ∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥)

+ min_{s≤2}{k_{1}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})}} + min_{s≤1}{k_{1}^{s+2}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})}},

and

F_{2}^{(2)}(u) = ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}u_{l}∥ + ∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥)

+ max_{1≤l≤2} min_{s≤2}{k_{l}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{l},L_{2})}}.

Remark. If one would settle for second order convergence (when q = 1), then Πf(U) in (3.2) could be taken to be a linear interpolant, which could even be computed at interior points of the interval I_{n−1} or I_{n}, with corresponding simplifications of the analysis and implementation.

4. Proof of Theorem 3.1. For the proof of Theorem 3.1 we shall compare the solution U of (3.1) with the solution V of the discrete linear problem (2.1) with f(x, t) replaced by f(x, t, u(x, t)). That is, V |_{Ω×I_{n}} ∈ V_{n}^{q}, 0 < t_{n} ≤ t^{∗} is defined by V_{0}^{−} = u_{0} and

(4.1) ∫_{I_{n}}((V_{t}, X) + (∇V, ∇X)) dt + ([V]_{n−1}, X_{n−1}^{+}) = ∫_{I_{n}}(f(u), X) dt ∀X ∈ V_{n}^{q}.

Theorem 2.1 then shows that

(4.2)  ∥V_{n}^{−} − u_{n}∥ ≤ CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u), 0 < t_{n} ≤ t^{∗},

and the proof of Theorem 3.1 will be accomplished once we have shown that for θ = U − V we have

(4.3)  ∥θ∥_{L_{∞}(I_{n},L_{2})} ≤ CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u) + C max_{1≤j≤n} F_{j}^{(2q)}(u), 0 < t_{n} ≤ t^{∗}.

In order to prove this we shall show below that

(4.4)  ∥θ_{n}^{−}∥ ≤ CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u) + C max_{1≤j≤n} F_{j}^{(2q)}(u), 0 < t_{n} ≤ t^{∗}.

It is convenient to begin by demonstrating that (4.3) follows from (4.4). From (3.1) and (4.1) it follows that θ|_{Ω×I_{n}} ∈ V_{n}^{q}, 0 < t_{n} ≤ t^{∗}, satisfies θ_{0}^{−} = 0 and

(4.5) ∫_{I_{n}}((θ_{t}, X) + (∇θ, ∇X)) dt + ([θ]_{n−1}, X_{n−1}^{+})

= ∫_{I_{n}}(⟨Πf(U), X⟩ − (f(u), X)) dt ∀X ∈ V_{n}^{q}.

We split the integrand on the right-hand side into three terms,

\begin{aligned}⟨Πf(U), X⟩ − (f(u), X) & = ⟨Π[f(U) − f(u)], X⟩ \\ & + (⟨Πf(u), X⟩ − (Πf(u), X)) \\ & + ([Π − I]f(u), X). \end{aligned} (4.6)

These are estimated in the following three lemmas. The first term requires a uniform Lipschitz condition on the nonlinearity f, which holds in view of assumption (1.2).

Lemma 4.1. Let U|_{Ω×I_{l}} ∈ V_{l}^{q} for j − i − p ≤ l ≤ j − i and X ∈ V_{j}^{q}. Then

|⟨Π_{j−i}^{p}[f(U)−f(u)](t), X(t)⟩| ≤ CM ∑_{l=j−i−p}^{j−i}(∥U_{l}^{−}−u_{l}∥+∥h_{l,j}^{2}D_{x}^{2}u_{l}∥)∥X(t)∥, t ∈ I_{j}.

Proof. Using the uniform boundedness of the Lagrangian basis functions,

(4.7)  |ϕ_{l}(t)| ≤ C, t ∈ I_{j},

under the assumption (3.3), and employing the Lipschitz condition for f, we have

\begin{aligned}|⟨Π_{j−i}^{p}[f(U) − f(u)](t), X(t)⟩| & = | ∑_{l=j−i−p}^{j−i} ϕ_{l}(t) ([f(U)_{l} − f(u)_{l}], X(t))_{S_{l}∪S_{j}} | \\ & ≤ CM( ∑_{l=j−i−p}^{j−i} ∥U_{l}^{−} − u_{l}∥_{S_{l}∪S_{j}})∥X(t)∥_{S_{l}∪S_{j}}. \end{aligned}

But ∥ · ∥_{S_{l}∪S_{j}} and ∥ · ∥ are uniformly equivalent norms on S_{l} ∪ S_{j}, and X(t) ∈ S_{j} ⊂ S_{l} ∪ S_{j}, so that ∥X(t)∥_{S_{l}∪S_{j}} ≤ C∥X(t)∥. Similarly, with J_{S_{l}∪S_{j}} : C(\overline{Ω}) → S_{l} ∪ S_{j} the Lagrange interpolation operator, we have

\begin{aligned}∥U_{l}^{−} − u_{l}∥_{S_{l}∪S_{j}} & = ∥J_{S_{l}∪S_{j}}(U_{l}^{−} − u_{l})∥_{S_{l}∪S_{j}} ≤ C∥U_{l}^{−} − J_{S_{l}∪S_{j}}u_{l}∥ \\ & ≤ C (∥U_{l}^{−} − u_{l}∥ + ∥u_{l} − J_{S_{l}∪S_{j}}u_{l}∥), \end{aligned}

since U_{l}^{−} ∈ S_{l} ⊂ S_{l} ∪ S_{j}. In view of a well-known error bound for J_{S_{l}∪S_{j}}, this proves the lemma. □

The second term on the right of (4.6) involves the error in spatial quadrature.

Lemma 4.2. Let X ∈ V_{j}^{q}. Then

|⟨[Π_{j−i}^{p}f(u)](t), X(t)⟩ − ([Π_{j−i}^{p}f(u)](t), X(t))|

≤ C ∑_{l=j−i−p}^{j−i}(∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥)∥X(t)∥_{H^{1}}, t ∈ I_{j}.

Proof. Using (4.7) we have

\begin{aligned}|⟨[Π_{j−i}^{p}f(u)](t), X(t)⟩ − ([Π_{j−i}^{p}f(u)](t), X(t))| & = | ∑_{l=j−i−p}^{j−i} ϕ_{l}(t) ϵ_{l,j}(f(u)_{l}, X(t))| \\ & ≤ C ∑_{l=j−i−p}^{j−i} |ϵ_{l,j}(f(u)_{l}, X(t))|, \end{aligned}

where ϵ_{l,j}(·, ·) = (·, ·)_{S_{l}∪S_{j}} − (·, ·) is the quadrature error, for which we have

|ϵ_{l,j}(f, χ)| ≤ C(∥h_{l,j}^{2}D_{x}^{2}f∥ + ∥h_{l,j}^{2}D_{x}f∥)∥χ∥_{H^{1}}, f ∈ H^{2}(Ω), χ ∈ S_{l} ∪ S_{j}.

This follows by a simple modification of the proof of Lemma 2.3 of [4] or Lemma 1 on p. 170 of [5]. Since X(t) ∈ S_{j} ⊂ S_{l} ∪ S_{j} this proves the lemma. □

The third term in (4.6) involves the error in interpolation with respect to t. The proof of the following lemma is well-known and we omit it. The assumption (3.3) is used here again.

Lemma 4.3. Let m = max{i + p − 1, 0}. We have

∥(Π_{j−i}^{p} − I)f(u)∥_{L_{∞}(I_{j},L_{2})} ≤ C max_{j−m≤l≤j} min_{s≤p+1}{k_{l}^{s}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{l},L_{2})}}.

We can now estimate ∥θ∥_{L_{∞}(I_{n},L_{2})} by applying a simple energy argument to (4.5). We formulate this as a lemma.

Lemma 4.4. Let U ∈ V_{l}^{q} for j − i − p ≤ l ≤ j − i and assume that θ ∈ V_{j}^{q} satisfies

(4.8) ∫_{I_{j}}((θ_{t}, X) + (∇θ, ∇X)) dt + ([θ]_{j−1}, X_{j−1}^{+})

= ∫_{I_{j}}(⟨Π_{j−i}^{p}f(U), X⟩ − (f(u), X)) dt ∀X ∈ V_{j}^{q}.

Then with m = max{i + p − 1, 0} we have

∥θ∥_{L_{∞}(I_{j},L_{2})} ≤ C∥θ_{j−1}^{−}∥ + CMk_{j} ∑_{l=j−i−p}^{j−i}(∥U_{l}^{−} − u_{l}∥ + ∥h_{l,j}^{2}D_{x}^{2}u_{l}∥)

+ C\sqrt{k_{j}} ∑_{l=j−i−p}^{j−i}(∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥)

+ Ck_{j} max_{j−m≤l≤j} min_{s≤p+1}(k_{l}^{s}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{l},L_{2})}).

Proof. Since θ(t) is a polynomial in t of degree 0 or 1, we have

∥θ∥_{L_{∞}(I_{j},L_{2})} ≤ max{∥θ_{j}^{−}∥, ∥θ_{j−1}^{+}∥},

so it suffices to estimate the two terms on the right. Taking X = θ in (4.8), we have

∫_{I_{j}}((θ_{t}, θ) + (∇θ, ∇θ)) dt + ([θ]_{j−1}, θ_{j−1}^{+}) = ∫_{I_{j}}(⟨Π_{j−i}^{p}f(U), θ⟩ − (f(u), θ)) dt,

whence,

∥θ_{j}^{−}∥^{2} + ∥θ_{j−1}^{+}∥^{2} + ∥θ∥_{L_{2}(I_{j},H^{1})}^{2} ≤ C∥θ_{j−1}^{−}∥^{2} + C ∫_{I_{j}}(⟨Π_{j−i}^{p}f(U), θ⟩ − (f(u), θ)) dt.

In view of Lemmas 4.1, 4.2 and 4.3, we have here

∫_{I_{j}}(⟨Π_{j−i}^{p}f(U), X⟩ − (f(u), θ)) dt = ∫_{I_{j}}⟨Π_{j−i}^{p}[f(U) − f(u)], θ⟩ dt
+ ∫_{I_{j}}(⟨Π_{j−i}^{p}f(u), θ⟩ − (Π_{j−i}^{p}f(u), θ)) dt + ∫_{I_{j}}([Π_{j−i}^{p} − I]f(u), θ) dt
≤ Ck_{j}R_{1}∥θ∥_{L_{∞}(I_{j},L_{2})} + C\sqrt{k_{j}}R_{2}∥θ∥_{L_{2}(I_{j},H^{1})} + Ck_{j}R_{3}∥θ∥_{L_{∞}(I_{j},L_{2})}
≤ Cϵ^{−1}k_{j}^{2}(R_{1}^{2} + R_{3}^{2}) + Cϵ^{−1}k_{j}R_{2}^{2} + ϵ∥θ∥_{L_{∞}(I_{j},L_{2})}^{2} + ϵ∥θ∥_{L_{2}(I_{j},H^{1})}^{2},

where ϵ > 0 and

R_{1} = M ∑_{l=j−i−p}^{j−i}(∥U_{l}^{−} − u_{l}∥ + ∥h_{l,j}^{2}D_{x}^{2}u_{l}∥),

R_{2} = ∑_{l=j−i−p}^{j−i}(∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥),

R_{3} = max_{j−m≤l≤j} min_{s≤p+1}(k_{l}^{s}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{l},L_{2})}).

Hence, choosing ϵ small enough we obtain the desired result. □

We are now in a position to finish the proof that (4.4) implies (4.3). We have to distinguish between the general case and the startup terms. In the general case, i.e., when n ≥ m + 1, where m is defined in (3.4), we have Π = Π_{n−i}^{2q} in (4.5) and Lemma 4.4, (4.2) and (4.4) show

∥θ∥_{L_{∞}(I_{n},L_{2})} ≤ C max_{1≤j≤n}(∥θ_{j}^{−}∥ + ∥V_{j}^{−} − u_{j}∥ + ∑_{l=j−i−2q}^{j−i}(∥h_{l,j}^{2}D_{x}^{2}u_{l}∥

+ ∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥)

+ min_{s≤2q+1}(k_{j}^{s}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{j},L_{2})}))

≤ CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u) + C max_{1≤j≤n} F_{j}^{(2q)}(u), t_{m+1} ≤ t_{n} ≤ t^{∗}.

This is the desired result in the general case.

If q = 1 we also have to consider the startup terms. If q = 1, i = 0, and n = 1, then we have Π = Π_{1}^{1} in (4.5) and Lemma 4.4 shows

\begin{aligned}∥θ∥_{L_{∞}(I_{1},L_{2})} & ≤ C(∥θ_{1}^{−}∥ + ∥V_{1}^{−} − u_{1}∥ + ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}u_{l}∥ \\ & + ∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥)) \\ & + C min_{s≤2}\{k_{1}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})}\} \\ & ≤ CL_{1}E_{1}^{(2)}(u) + CF_{1}^{(2)}(u). \end{aligned} (4.9)

If q = 1, i = 1, then we first apply Lemma 4.4 with θ replaced by θ̃ = Ũ − V and with Π = Π_{0}^{0}. This gives

\begin{aligned}∥θ̃∥_{L_{∞}(I_{1},L_{2})} & ≤ C(∥h_{0,1}^{2}D_{x}^{2}u_{0}∥ + ∥h_{0,1}^{2}D_{x}^{2}f(u)(t_{0})∥ + ∥h_{0,1}^{2}D_{x}f(u)(t_{0})∥ \\ & + min_{s≤1}\{k_{1}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})}\}). \end{aligned}

Next we apply the lemma with θ = U − V and with Π_{1}^{1}f(Ũ) − f(u) on the right-hand side. This gives

∥θ∥_{L_{∞}(I_{1},L_{2})} ≤ C(k_{1}∥θ̃_{1}^{−}∥ + ∥V_{1}^{−} − u_{1}∥ + ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}u_{l}∥

+ ∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥))

+ C min_{s≤2}{k_{1}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})}}.

Using the above bound for θ̃ = Ũ − V we conclude

∥θ∥_{L_{∞}(I_{1},L_{2})} ≤ CL_{1}E_{1}^{(2)}(u) + CF_{1}^{(2)}(u).

We then continue to the interval I_{2} and use Lemma 4.4 with θ = U − V, Π = Π_{1}^{1}:

∥θ∥_{L_{∞}(I_{2},L_{2})} ≤ C(∥θ_{1}^{−}∥ + ∥V_{1}^{−} − u_{1}∥ + ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}u_{l}∥

+ ∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥))

+ C max_{1≤l≤2} min_{s≤2}{k_{l}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{l},L_{2})}}

≤ CL_{1}E_{1}^{(2)}(u) + CF_{2}^{(2)}(u).

Assumption (3.3) was also employed here. This completes the proof that (4.4) implies (4.3).

We now proceed to prove (4.4). Following [2] we represent ∥θ_{n}^{−}∥ by duality. For this purpose we define V^{q} = {V : V ∈ V_{n}^{q}, 0 < t_{n} ≤ t^{∗}}. We then note that, by summation of equation (4.5), we have

(4.10)  B(θ, X) = ∑_{j=1}^{n} ∫_{I_{j}}⟨Πf(U), X⟩ dt − ∫_{0}^{t_{n}}(f(u), X) dt ∀X ∈ V^{q},

where

\begin{aligned}B(V, W) & = ∑_{j=1}^{n} ∫_{I_{j}}((V_{t}, W) + (∇V, ∇W)) dt + ∑_{j=1}^{n−1}([V]_{j}, W_{j}^{+}) + (V_{0}^{+}, W_{0}^{+}) \\ & = ∑_{j=1}^{n} ∫_{I_{j}}(−(V, W_{t}) + (∇V, ∇W)) dt − ∑_{j=1}^{n−1}(V_{j}^{−}, [W]_{j}) + (V_{n}^{−}, W_{n}^{−}). \end{aligned}

Next we consider the discrete analog of the “backward problem”

−z_{t} − Δz = 0, 0 < t < t; z(t_{n}) = θ_{n}^{−}.

In view of the second form of B(·, ·) above, it is clear that the corresponding discrete problem consists in finding Z ∈ V^{q} such that

(4.11)  B(X, Z) = (X_{n}^{−}, θ_{n}^{−}) ∀X ∈ V^{q}.

The following stability bound is proved in [3]:

(4.12)  ∥Z∥_{L_{∞}([0,t_{n}],L_{2})} + ∥Z∥_{L_{2}([0,t_{n}],H^{1})} ≤ C∥θ_{n}^{−}∥.

Taking X = θ in (4.11) and using (4.10), we obtain

∥θ_{n}^{−}∥^{2} = B(θ, Z) = (∑_{j=1}^{m} + ∑_{j=m+1}^{n}) ∫_{I_{j}}(⟨Πf(U), Z⟩ − (f(u), Z)) dt,

where m is as in (3.4) and the first sum is empty if m = 0. Splitting the integrand on the right as in (4.6) and using Lemmas 4.1, 4.2 and 4.3, we obtain

| ∑_{j=m+1}^{n} ∫_{I_{j}}(⟨Πf(U), Z⟩ − (f(u), Z)) dt|

≤ C ∑_{j=m+1}^{n} k_{j} ∑_{l=j−i−2q}^{j−i}(∥θ_{l}^{−}∥ + ∥V_{l}^{−} − u_{l}∥ + ∥h_{l,j}^{2}D_{x}^{2}u_{l}∥)∥Z∥_{L_{∞}([0,t_{n}],L_{2})}

+ C\sqrt{t_{n}} max_{1≤j≤n} ∑_{l=j−i−2q}^{j−i}(∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥)∥Z∥_{L_{2}([0,t_{n}],H^{1})}

+ Ct_{n} max_{1≤j≤n} min_{s≤2q+1}(k_{j}^{s}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{j},L_{2})})∥Z∥_{L_{∞}([0,t_{n}],L_{2})}

≤ C((1 − i)k_{n}∥θ_{n}^{−}∥ + ∑_{j=1}^{n−1} k_{j}∥θ_{j}^{−}∥ + L_{n} max_{1≤j≤n} E_{j}^{(2q)}(u) + max_{1≤j≤n} F_{j}^{(2q)}(u)))∥θ_{n}^{−}∥.

Here we have also used (4.12), (4.2) and the fact that θ_{0}^{−} = 0. The factor 1 − i in the first term on the right means that this term vanishes for the lagged schemes (i = 1).

If m = 1, i.e., q = 1, i = 0, we have also the term

|∫_{I_{1}}(⟨Π_{1}^{1}f(U), Z⟩ − (f(u), Z)) dt|

≤ Ck_{1}(∥θ_{1}^{−}∥ + ∥V_{1}^{−} − u_{1}∥ + ∑_{l=0}^{1} ∥h_{l,1}^{2}D_{x}^{2}u_{l}∥)∥Z∥_{L_{∞}(I_{1},L_{2})}

+ C\sqrt{k_{1}} ∑_{l=0}^{1}(∥h_{l,1}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,1}^{2}D_{x}f(u)_{l}∥)∥Z∥_{L_{2}(I_{1},H^{1})}

+ C min_{s≤2}(k_{1}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{1},L_{2})})∥Z∥_{L_{∞}(I_{1},L_{2})}

≤ (CL_{1}E_{1}^{(2)}(u) + CF_{1}^{(2)}(u))∥θ_{n}^{−}∥,

where we have used the estimate of ∥θ_{1}^{−}∥ from (4.9).

If m = 2, i.e., q = 1, i = 1, we have the terms

|∫_{I_{1}}(⟨Π_{1}^{1}f(Ũ), Z⟩ − (f(u), Z)) dt| + |∫_{I_{2}}(⟨Π_{1}^{1}f(U), Z⟩ − (f(u), Z)) dt|

≤ Ck_{1}(∥θ̃_{1}^{−}∥ + ∥V_{1}^{−} − u_{1}∥ + ∑_{l=0}^{1} ∥h_{l,1}^{2}D_{x}^{2}u_{l}∥)∥Z∥_{L_{∞}(I_{1},L_{2})}

+ Ck_{2}(∥θ_{1}^{−}∥ + ∥V_{1}^{−} − u_{1}∥ + ∑_{l=0}^{1} ∥h_{l,1}^{2}D_{x}^{2}u_{l}∥)∥Z∥_{L_{∞}(I_{2},L_{2})}

+ C\sqrt{t_{2}} max_{1≤j≤2} ∑_{l=0}^{1}(∥h_{l,j}^{2}D_{x}^{2}f(u)_{l}∥ + ∥h_{l,j}^{2}D_{x}f(u)_{l}∥)∥Z∥_{L_{2}([0,t_{2}],H^{1})}

+ C max_{1≤l≤2} min_{s≤2}(k_{l}^{s+1}∥D_{t}^{s}f(u)∥_{L_{∞}(I_{l},L_{2})})∥Z∥_{L_{∞}([0,t_{2}],L_{2})}

≤ (CL_{1}E_{1}^{(2)}(u) + C max_{1≤j≤2} F_{j}^{(2)}(u))∥θ_{n}^{−}∥,

where we have used the already proven bounds for ∥θ̃_{1}^{−}∥ and ∥θ_{1}^{−}∥.

Altogether we now have

∥θ_{n}^{−}∥ ≤ (1 − i)Ck_{n}∥θ_{n}^{−}∥ + C ∑_{j=1}^{n−1} k_{j}∥θ_{j}^{−}∥ + CL_{n} max_{1≤j≤n} E_{j}^{(2q)}(u) + C max_{1≤j≤n} F_{j}^{(2q)}(u),

and (4.4) follows by Gronwall’s Lemma. If i = 0 we first have to eliminate the first term on the right by taking k_{n} small. This completes the proof of Theorem 3.1.

5. Example. We conclude this paper by providing an example showing the implementation of one of our algorithms in the case q = 1. We solve the equation u_{t} − Δu = 10(u − u^{3}) with homogeneous Dirichlet boundary conditions and Ω = (0, 1) × (0, 1). For the spatial discretization we use the piecewise linear finite element method computed on the fixed mesh shown in Figure 1 with mesh spacing h = 1/m. We number the nodes as indicated so that there are m = 1/h nodes on a side and m^{2} nodes in total.

Figure 1. Mesh for the sample computation.

Associated with this mesh are the finite element space S, the mass matrix A, A_{ij} = (φ_{i}, φ_{j}) and stiffness matrix B, B_{ij} = (∇φ_{i}, ∇φ_{j}), where {φ_{i}}_{i=1}^{m^{2}} is the standard computational basis for S. We will also need the lumped mass matrix Ā, Ā_{ij} = (φ_{i}, φ_{j})_{S}, see [4] or [5].

We derive the equations for the semi-implicit q = 1, i = 1 scheme, which in some sense is the most complicated of our algorithms. We also have experience computing with q = 0, i = 0, 1 and q = 1, i = 0, 1 in one and two dimensions using fixed spatial meshes, and with q = 1, i = 0 in the fully adaptive two dimensional code TRANSI.

Throughout the startup procedure and the rest of the steps the same discrete system results from the left-hand side of (3.1). We let U ∈ V_{n}^{q} be given by

U(t) = U_{n−1}^{+} \frac{t − t_{n}}{−k_{n}} + U_{n}^{−} \frac{t − t_{n−1}}{k_{n}} ,

where U_{n−1}^{+} and U_{n}^{−} now denote vectors of values at the spatial nodes, that is, U_{n}^{±} = ((U_{n}^{±})_{1}, (U_{n}^{±})_{2}, …, (U_{n}^{±})_{m^{2}})^{T}, corresponding to the basis of S. (3.1) yields the pair of equations

A\frac{U_{n}^{−} − U_{n−1}^{+}}{k_{n}} ∫_{I_{n}} \frac{t − t_{n}}{−k_{n}} dt + BU_{n−1}^{+} ∫_{I_{n}} \frac{(t − t_{n})^{2}}{k_{n}^{2}} dt

+ BU_{n}^{−} ∫_{I_{n}} \frac{(t − t_{n})(t − t_{n−1})}{−k_{n}^{2}} dt + AU_{n−1}^{+} = AU_{n−1}^{−} + Ā ∫_{I_{n}} Πf(U)\frac{t − t_{n}}{−k_{n}} dt,

and

A\frac{U_{n}^{−} − U_{n−1}^{+}}{k_{n}} ∫_{I_{n}} \frac{t − t_{n−1}}{k_{n}} dt + BU_{n−1}^{+} ∫_{I_{n}} \frac{(t − t_{n})(t − t_{n−1})}{−k_{n}^{2}} dt

+ BU_{n}^{−} ∫_{I_{n}} \frac{(t − t_{n})^{2}}{k_{n}^{2}} dt = Ā ∫_{I_{n}} Πf(U)\frac{t − t_{n−1}}{k_{n}} dt,

where f(U) = (f((U)_{1}), f((U)_{2}), …, f((U)_{m^{2}}))^{T}. This is equivalent to

[ \begin{matrix}\frac{1}{2}A + \frac{1}{3}k_{n}B & \frac{1}{2}A + \frac{1}{6}k_{n}B \\ −\frac{1}{2}A + \frac{1}{6}k_{n}B & \frac{1}{2}A + \frac{1}{3}k_{n}B \end{matrix} ] [ \begin{matrix}U_{n−1}^{+} \\ U_{n}^{−} \end{matrix} ] = [ \begin{matrix}0 & A \\ 0 & 0 \end{matrix} ] [ \begin{matrix}U_{n−2}^{+} \\ U_{n−1}^{−} \end{matrix} ]

+ [ \begin{matrix}Ā ∫_{I_{n}} Πf(U)\frac{t−t_{n}}{−k_{n}} dt \\ Ā ∫_{I_{n}} Πf(U)\frac{t−t_{n−1}}{k_{n}} dt \end{matrix} ] ,

for n = 1, 2, …, with U_{0}^{+} := 0. We let

U_{n} := [ \begin{matrix}U_{n−1}^{+} \\ U_{n}^{−} \end{matrix} ] , B_{n} := [ \begin{matrix}\frac{1}{2}A + \frac{1}{3}k_{n}B & \frac{1}{2}A + \frac{1}{6}k_{n}B \\ −\frac{1}{2}A + \frac{1}{6}k_{n}B & \frac{1}{2}A + \frac{1}{3}k_{n}B \end{matrix} ] , A := [ \begin{matrix}0 & A \\ 0 & 0 \end{matrix} ] .

Then, in order of application, the equations are

B_{1}Ũ_{1} = AU_{0} + [ \begin{matrix}0 & \frac{1}{2}k_{1}Ā \\ 0 & \frac{1}{2}k_{1}Ā \end{matrix} ] [ \begin{matrix}0 \\ f(U_{0}^{−}) \end{matrix} ] ,

B_{1}U_{1} = AU_{0} + [ \begin{matrix}\frac{1}{3}k_{1}Ā & \frac{1}{6}k_{1}Ā \\ \frac{1}{6}k_{1}Ā & \frac{1}{3}k_{1}Ā \end{matrix} ] [ \begin{matrix}f(U_{0}^{−}) \\ f(Ũ_{1}^{−}) \end{matrix} ] ,

B_{2}U_{2} = AU_{1} + [ \begin{matrix}−\frac{k_{2}^{2}}{6k_{1}} Ā & \frac{k_{2}}{k_{1}} (\frac{k_{2}}{6} + \frac{k_{1}}{2} )Ā \\ −\frac{k_{2}^{2}}{3k_{1}} Ā & \frac{k_{2}}{k_{1}} (\frac{k_{2}}{3} + \frac{k_{1}}{2} )Ā \end{matrix} ] [ \begin{matrix}f(U_{0}^{−}) \\ f(U_{1}^{−}) \end{matrix} ] ,

and finally, in general,

B_{n}U_{n} = AU_{n−1}

+ [ \begin{matrix}\frac{k_{n}^{2}}{k_{n−2}(k_{n−2}+k_{n−1})}(\frac{k_{n}}{12} + \frac{k_{n−1}}{6} )Ā & −\frac{k_{n}^{2}}{k_{n−2}k_{n−1}} (\frac{k_{n}}{12} + \frac{k_{n−1}}{6} + \frac{k_{n−2}}{6} )Ā \\ \frac{k_{n}^{2}}{k_{n−2}(k_{n−2}+k_{n−1})}(\frac{k_{n}}{4} + \frac{k_{n−1}}{3} )Ā & −\frac{k_{n}^{2}}{k_{n−2}k_{n−1}} (\frac{k_{n}}{4} + \frac{k_{n−1}}{3} + \frac{k_{n−2}}{3} )Ā \end{matrix}

\begin{matrix}\frac{k_{n}}{k_{n−1}(k_{n−2}+k_{n−1})}(\frac{k_{n}^{2}}{12} + \frac{k_{n}k_{n−1}}{3} + \frac{k_{n}k_{n−2}}{6} + \frac{k_{n−1}^{2}}{2} + \frac{k_{n−1}k_{n−2}}{2} )Ā \\ \frac{k_{n}}{k_{n−1}(k_{n−2}+k_{n−1})}(\frac{k_{n}^{2}}{4} + \frac{2k_{n}k_{n−1}}{3} + \frac{k_{n}k_{n−2}}{3} + \frac{k_{n−1}^{2}}{2} + \frac{k_{n−1}k_{n−2}}{2} )Ā \end{matrix} ] [ \begin{matrix}f(U_{n−3}^{−}) \\ f(U_{n−2}^{−}) \\ f(U_{n−1}^{−}) \end{matrix} ] .

Implementing these in code is now straightforward. There is one point that should be discussed, namely how to solve the associated systems of equations. In particular, the q = 1 formulas yield nonsymmetric linear systems which can be expensive to solve. (This difficulty is common in systems arising from implicit high-order Runge-Kutta schemes, which are related to these formulas.) Both B and A have a “principal” tridiagonal band as well as two other bands located O(1/h) above and below the diagonal. Thus, for the q = 0 scheme (backward Euler), using just this band structure without further rearrangement yields an operations count of O(1/h^{4}) for solving the system once. In our case B_{n} has a block 2 × 2 band structure. Solving this system directly however leads to an operations count of O(1/h^{6}) (as if the matrix were full) because of fill-in. Rearranging the matrix can alleviate this load. For example, it is straightforward to rearrange the matrix to achieve an operations count of O(1/h^{4}) once again. First permute the rows of B_{n} into the order 1, m^{2} + 1, 2, m^{2} + 2, ..., m^{2}, 2m^{2} and then do the same by columns. The resulting matrix has a band of width 7 down the diagonal and bands of width 5 located O(1/h) above and below the diagonal.

The following computation was made on a 35 × 35 mesh with initial data equal to 1 inside the disc (x − \frac{1}{2})^{2} + (y − \frac{1}{2})^{2} ≤ \frac{1}{16} and equal to 0 outside. The solution eventually becomes all 0 at a fairly slow rate. See Figure 2.

References

1. T. Dupont, Mesh modification for evolution equations, Math. Comp. 39 (1982), 85–107.

2. K. Eriksson and C. Johnson, Adaptive finite element methods for parabolic problems I: a linear model problem, SIAM J. Numer. Anal. 28 (1991), 43–77.

3. K. Eriksson, C. Johnson and V. Thomée, Time discretization of parabolic problems by the discontinuous Galerkin method, M^{2}AN 19 (1985), 611–643.

4. Y.-Y. Nie and V. Thomée, A lumped mass finite-element method with quadrature for a non-linear parabolic problem, IMA J. Numer. Anal. 5 (371–396).

5. V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Lecture Notes in Mathematics, vol. 1054, Springer-Verlag, 1984.

Figure 2. The discontinuous Galerkin q = 1, i = 0 approximation for the equation u_{t} − Δu = 10(u − u^{3}).

School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332, USA E-mail address: estep@math.gatech.edu

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

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