A FINITE ELEMENT MODEL FOR THE

TIME-DEPENDENT JOULE HEATING PROBLEM^{∗}

Charles M. Elliott and Stig Larsson

Abstract. We study a spatially semidiscrete and a completely discrete finite element model for a nonlinear system consisting of an elliptic and a parabolic partial differential equation describing the electric heating of a conducting body. We prove error bounds of optimal order under minimal regularity assumptions when the number of spatial variables d ≤ 3. We establish the existence of solutions with the required regularity over arbitrarily long intervals of time when d ≤ 2.

1. Introduction

In this note we consider the numerical approximation by the finite element method of the following nonlinear elliptic-parabolic system

(1.1)  \begin{matrix}u_{t} − Δu = σ(u)|∇ϕ|^{2}, \\ − ∇ · (σ(u)∇ϕ) = 0, \end{matrix}  x ∈ Ω, t ∈ [0, T],

where u = u(x, t), ϕ = ϕ(x, t), u_{t} = ∂u/∂t, ∇ denotes the gradient with respect to the x-variables and Δ = ∇ · ∇ is the Laplacian. These differential equations are studied for t in a finite interval [0, T] and for x in a bounded convex polygonal domain Ω in R^{d}, d = 1, 2 or 3, together with initial and boundary conditions

(1.2)  \begin{matrix}u(x, t) = 0, ϕ(x, t) = g(x, t), & x ∈ ∂Ω, t ∈ [0, T], \\ u(x, 0) = u_{0}(x), & x ∈ Ω. \end{matrix}

We make the assumption that the function σ ∈ C^{2}(R) and that, for some κ, K > 0 and all s ∈ R,

(1.3)  0 < κ ≤ σ(s) ≤ K,  |σ′(s)| + |σ′′(s)| ≤ K.

This system models the electric heating of a conducting body [5] with u being the temperature, ϕ the electric potential, and σ the temperature-dependent electric conductivity.

Let (·, ·) and ∥ · ∥ denote the inner product and norm in L_{2} = L_{2}(Ω), and H^{1} = H^{1}(Ω) = {u ∈ L_{2} : |∇u| ∈ L_{2}}, H_{0}^{1} = {u ∈ H^{1} : u|_{∂Ω} = 0} be the standard Sobolev spaces. The finite element method is based on the weak formulation of the

Typeset by AMS-TEX

1991 Mathematics Subject Classification. 65N30, 65N15, 35K60.

Key words and phrases. Joule heating, nonlinear, elliptic, parabolic, finite element, backward Euler, existence, regularity.

This research was partially supported by the SERC grant number GR/F85659.

∗Published in Math. Comp. 64 (1995), 1433–1453.

above initial boundary value problem, where we seek u(t) ∈ H_{0}^{1}, ϕ(t) ∈ H^{1} with ϕ(t) − g(t) ∈ H_{0}^{1} such that

(1.4)  \begin{matrix}(u_{t}, χ) + (∇u, ∇χ) = (σ(u)|∇ϕ|^{2}, χ),  ∀χ ∈ H_{0}^{1}, t ∈ [0, T], \\ u(0) = u_{0}, \end{matrix}

and

(1.5)  (σ(u)∇ϕ, ∇χ) = 0,  ∀χ ∈ H_{0}^{1}, t ∈ [0, T].

Let {S_{h}}_{h>0} be a family of approximating subspaces of H^{1}, where each space S_{h} consists of continuous piecewise linear polynomials with respect to a triangulation of Ω with maximum mesh-width h. With each S_{h} we associate the subspace S̊_{h} = {u_{h} ∈ S_{h} : u_{h}|_{∂Ω} = 0}. We assume that the family of triangulations is such that the standard interpolation error estimates [4, Theorem 3.2.1] and inverse estimates [4, Theorem 3.2.6] hold.

We first consider a semidiscrete approximation: find u_{h}(t) ∈ S̊_{h}, ϕ_{h}(t) ∈ S_{h} with ϕ_{h}(t) − π_{h}g(t) ∈ S̊_{h} such that

(1.6)  \begin{matrix}(u_{h,t}, χ) + (∇u_{h}, ∇χ) = (σ(u_{h})|∇ϕ_{h}|^{2}, χ),  ∀χ ∈ S̊_{h}, t ∈ [0, T], \\ u_{h}(0) = u_{h0}, \end{matrix}

and

(1.7)  (σ(u_{h})∇ϕ_{h}, ∇χ) = 0,  ∀χ ∈ S̊_{h}, t ∈ [0, T],

where π_{h} : C(Ω̄) → S_{h} denotes the standard Lagrangian interpolation operator and u_{h0} ∈ S̊_{h} is an appropriate approximation of u_{0}. For this method we prove an error estimate of the form

∥u_{h}(t) − u(t)∥ + ∥ϕ_{h}(t) − ϕ(t)∥ ≤ C(u, ϕ, T)h^{2}, t ∈ [0, T],

(see Theorem 3.1 below) under a certain assumption about the regularity of the exact solutions u and ϕ. This assumption is essentially the same as in the standard error analysis for the corresponding linear elliptic and parabolic problems. The main difficulty here concerns the treatment of the gradient dependent nonlinearity: one has to deal with the expression

σ(u_{h})|∇ϕ_{h}|^{2} − σ(u)|∇ϕ|^{2} = σ(u_{h})∇(ϕ_{h} + ϕ) · ∇(ϕ_{h} − ϕ) + (σ(u_{h}) − σ(u))|∇ϕ|^{2},

where ∇(ϕ_{h} − ϕ) is formally only O(h), and where ∇ϕ and ∇ϕ_{h} enter in a nonlinear way. These difficulties are handled by means of a duality argument and by taking advantage of parabolic smoothing. In particular, we avoid using a maximum norm bound for ∇ϕ_{h}, which would be difficult to obtain.

We also consider a completely discrete scheme based on the backward Euler method with semi-implicit linearization: find U_{n} ∈ S̊_{h}, Φ_{n} ∈ S_{h} with Φ_{n} − π_{h}g(t_{n}) ∈ S̊_{h} such that

(1.8) \begin{matrix}(∂_{n}U_{n}, χ) + (∇U_{n}, ∇χ) = (σ(U_{n−1})|∇Φ_{n−1}|^{2}, χ), ∀χ ∈ S̊_{h}, t_{n} ∈ (0, T], \\ U_{0} = u_{h0}, \end{matrix}

and

(1.9)  (σ(U_{n})∇Φ_{n}, ∇χ) = 0,  ∀χ ∈ S̊_{h}, t_{n} ∈ [0, T).

Here ∂_{n}U_{n} = (U_{n} − U_{n−1})/k, t_{n} = nk, n = 0, 1, 2, …, and k is the timestep. For this scheme we show in Theorem 3.3 that

∥U_{n} − u(t_{n})∥ + ∥Φ_{n} − ϕ(t_{n})∥ ≤ C(u, ϕ, T)(h^{2} + k), t_{n} ∈ [0, T],

again under the same regularity requirement as for linear problems.

We begin the error analysis in Section 2 by recalling some results about linear elliptic and parabolic finite element problems. The nonlinear error analysis is carried out in Section 3 assuming that the number of spatial variables d ≤ 3, and that the exact solutions have minimal regularity. Finally, in Section 4 we prove the global existence of solutions with the required regularity when d ≤ 2. Our argument here builds upon the techniques of Cimatti [5], who showed the existence of weak solutions. We are not aware of any existence and regularity result in the three-dimensional case.

There is a vast literature on finite element methods for nonlinear elliptic and parabolic problems. For example, we mention the work [6], [7] on the porous media equations, which are similar to the Joule heating problem. Roughly speaking, the porous media equations are (1.1) with the term σ(u)|∇ϕ|^{2} replaced by ∇ϕ · ∇u, where u is a concentration, ϕ is the pressure, and ∇ϕ is the velocity. In [6], [7] the equation for ϕ is solved by a mixed method where both ϕ and ∇ϕ are approximated to order O(h^{2}), so that some difficulties that we address here are partly avoided there.

After the present work was finished we became aware of the paper [18], which addresses the same problem as we do, but obtains non-optimal results.

Throughout this work we use the notation ∥u∥_{m,p} = (∑_{|α|≤m} ∥D^{α}u∥_{L_{p}}^{p})^{1/p} for the norm in the standard Sobolev space W_{p}^{m} = W_{p}^{m}(Ω) with the usual modification for p = ∞, and with the exception that ∥ · ∥ and (·, ·) denote the norm and inner product in L_{2}. We also write H^{m} = W_{2}^{m} when p = 2.

2. Linear error analysis

In this section we collect some facts about linear elliptic and parabolic finite element problems that we will need in the sequel. Since ∂Ω is a convex polygon it is wellknown [9] that the Laplacian Δ is an isomorphism from H^{2} ∩ H_{0}^{1} onto L_{2} and we let Δ^{−1} denote its inverse. Let R_{h} : H_{0}^{1} → S̊_{h} be defined by the equation

(2.1)  (∇R_{h}u, ∇χ) = (∇u, ∇χ),  ∀u ∈ H_{0}^{1}, χ ∈ S̊_{h}.

From the standard error analysis [4, Theorems 3.2.2, 3.2.5] for linear elliptic finite element problems we quote the error estimates

(2.2)  ∥(R_{h} − I)u∥ + h∥(R_{h} − I)u∥_{1,2} ≤ Ch^{2}∥u∥_{2,2},  ∀u ∈ H^{2} ∩ H_{0}^{1}.

We denote by Δ_{h} : S̊_{h} → S̊_{h} the discrete Laplacian defined by

(−Δ_{h}χ, η) = (∇χ, ∇η), ∀χ, η ∈ S̊_{h},

and we let E_{h}(t) = exp(tΔ_{h}) be the analytic semigroup generated by Δ_{h}, and P_{h} : L_{2} → S̊_{h} the orthogonal projector. It is well-known that E_{h}(t)P_{h} satisfies the following bounds:

(2.3) ∥E_{h}(t)P_{h}ψ∥ + t^{1/2}∥E_{h}(t)P_{h}ψ∥_{1,2} + t∥Δ_{h}E_{h}(t)P_{h}ψ∥ ≤ C∥ψ∥,  t > 0,

for ψ ∈ L_{2}, where C is independent of h and t, reflecting the uniform analyticity of the evolution operator. In a similar way, for the discrete evolution operator E_{kh}^{n} = (I − kΔ_{h})^{−n} associated with the backward Euler method, we have

(2.4)  ∥E_{kh}^{n}P_{h}ψ∥ + t_{n}^{1/2}∥E_{kh}^{n}P_{h}ψ∥_{1,2} + t_{n}∥Δ_{h}E_{kh}^{n}P_{h}ψ∥ ≤ C∥ψ∥,  t_{n} > 0.

We may now state and prove error bounds for linear parabolic finite element problems. Such results are common in the literature, but a particular feature of the error bounds presented here is that the regularity requirement is optimal and expressed in a form that is suitable for our regularity analysis in Section 4. Similar results are proved in Chapter 2 of [17] for spatially semidiscrete approximations of the linear homogeneous problem, but are not readily available for the nonhomogeneous problem and completely discrete schemes. Moreover, our proof technique is different from that of [17]; being based on (2.3) and (2.4) this technique will also be used in our nonlinear error analysis below.

Theorem 2.1. (a) Suppose that u(t) ∈ H_{0}^{1} is the solution of the linear heat equation

(2.5)  \begin{matrix}(u_{t}, χ) + (∇u, ∇χ) = (f(t), χ),  ∀χ ∈ H_{0}^{1}, t > 0, \\ u(0) = u_{0}, \end{matrix}

and that u_{h}(t) ∈ S̊_{h} satisfies

\begin{matrix}(u_{h,t}, χ) + (∇u_{h}, ∇χ) = (f(t), χ), ∀χ ∈ S̊_{h}, t > 0, \\ u_{h}(0) = u_{h0}. \end{matrix}

Then, for t > 0, we have

(2.6)  ∥u_{h}(t) − u(t)∥ ≤ C∥u_{h0} − u_{0}∥ + Ch^{2} sup_{0<s≤t}(∥u(s)∥_{2,2} + s∥u_{t}(s)∥_{2,2}),

provided that the solution u has the regularity implied by the norms on the right hand side. (b) If U_{n} ∈ S̊_{h} satisfies

(2.7)  \begin{matrix}(∂_{n}U_{n}, χ) + (∇U_{n}, ∇χ) = (f(t_{n}), χ),  ∀χ ∈ S̊_{h}, t_{n} > 0, \\ U_{0} = u_{h0}, \end{matrix}

then, for t_{n} > 0, we have

\begin{aligned}∥U_{n} − u(t_{n})∥ & ≤ C∥u_{h0} − u_{0}∥ \\ & + Ch^{2} sup_{0<s≤t_{n}}(∥u(s)∥_{2,2} + ∥u_{t}(s)∥ + s∥u_{t}(s)∥_{2,2}) \\ & + Ck sup_{0<s≤t_{n}}(∥Δ^{−1}u_{tt}(s)∥ + s∥u_{tt}(s)∥). \end{aligned} (2.8)

Proof. We prove (2.8) only; (2.6) can be proved in a similar way; indeed, it essentially follows from (2.8) as k → 0. See also [17, Lemma 4 of Chapter 2] for a different proof of (2.6). Let for simplicity

F (u) = sup_{0<s≤t_{n}}(∥u(s)∥_{2,2} + ∥u_{t}(s)∥ + s∥u_{t}(s)∥_{2,2}),
G(u) = sup_{0<s≤t_{n}}(∥Δ^{−1}u_{tt}(s)∥ + s∥u_{tt}(s)∥).

A FINITE ELEMENT MODEL FOR THE JOULE HEATING PROBLEM 5

We write

e_{n} ≡ U_{n} − u(t_{n}) = (U_{n} − R_{h}u(t_{n})) + (R_{h}u(t_{n}) − u(t_{n})) ≡ θ_{n} + ρ_{n},

and the required estimate for ρ_{n} follows from (2.2):

(2.9)  ∥ρ_{n}∥ ≤ Ch^{2}∥u(t_{n})∥_{2,2} ≤ Ch^{2}F (u).

Moreover, we have

(2.10) ∥∂_{n}ρ_{n}∥ ≤ Ch^{2}t_{n}^{−1}F (u), t_{n} > 0,

because

∥∂_{n}ρ_{n}∥ ≤ Ch^{2}t_{n−1}^{−1} max_{t_{n−1}≤t≤t_{n}}(t∥u_{t}(t)∥_{2,2}) ≤ Ch^{2}t_{n}^{−1}F (u),

for t_{n} ≥ t_{2}, and

∥∂_{1}ρ_{1}∥ ≤ Ch^{2}k^{−1} max_{0≤t≤k} ∥u(t)∥_{2,2} ≤ Ch^{2}t_{1}^{−1}F (u).

The remaining term θ_{n} belongs to S̊_{h} and using (2.7), (2.5) and (2.1) we find that it satisfies the equation

∂_{n}θ_{n} − Δ_{h}θ_{n} = P_{h}(−∂_{n}ρ_{n} + ω_{n}),

where ω_{n} = u_{t}(t_{n}) − ∂_{n}u(t_{n}). Hence, by Duhamel’s principle,

θ_{n} = E_{kh}^{n}θ_{0} + k ∑_{j=1}^{n} E_{kh}^{n−j+1}P_{h}(−∂_{j}ρ_{j} + ω_{j}).

Let [n/2] be the integer part of n/2. Summation by parts gives

−k ∑_{j=1}^{[n/2]} E_{kh}^{n−j+1}P_{h}∂_{j}ρ_{j} = E_{kh}^{n}P_{h}ρ_{0} − E_{kh}^{n−[n/2]}P_{h}ρ_{[n/2]} + k ∑_{j=1}^{[n/2]}(∂_{j}E_{kh}^{n−j})P_{h}ρ_{j},

where ∂_{j}E_{kh}^{n−j} = −Δ_{h}E_{kh}^{n−j+1}, so that

θ_{n} = E_{kh}^{n}P_{h}e_{0} − E_{kh}^{n−[n/2]}P_{h}ρ_{[n/2]} − k ∑_{j=1}^{[n/2]} Δ_{h}E_{kh}^{n−j+1}P_{h}ρ_{j}

− k ∑_{j=[n/2]+1}^{n} E_{kh}^{n−j+1}P_{h}∂_{j}ρ_{j} + k ∑_{j=1}^{[n/2]} Δ_{h}E_{kh}^{n−j+1}(Δ_{h}^{−1}P_{h} − P_{h}Δ^{−1})ω_{j}

+ k ∑_{j=1}^{[n/2]} Δ_{h}E_{kh}^{n−j+1}P_{h}Δ^{−1}ω_{j} + k ∑_{j=[n/2]+1}^{n} E_{kh}^{n−j+1}P_{h}ω_{j} ≡ ∑_{i=1}^{7} R_{i}.

We proceed to estimate the seven terms on the right hand side. Using the smoothing property (2.4), the error bounds (2.9) and (2.10), we have

∑_{i=1}^{4} ∥R_{i}∥ ≤ C(∥e_{0}∥ + ∥ρ_{[n/2]}∥) + Ck ∑_{j=1}^{[n/2]} t_{n−j+1}^{−1}∥ρ_{j}∥ + Ck ∑_{j=[n/2]+1}^{n} ∥∂_{j}ρ_{j}∥

≤ C∥e_{0}∥ + Ch^{2}F (u)(1 + k ∑_{j=1}^{[n/2]} t_{n−j+1}^{−1} + k ∑_{j=[n/2]+1}^{n} t_{j}^{−1})

≤ C∥e_{0}∥ + Ch^{2}F (u).

For the fifth term we use the fact that Δ_{h}^{−1}P_{h} − Δ^{−1} = (R_{h} − I)Δ^{−1}, so that by (2.2) and elliptic regularity

\begin{aligned}∥(Δ_{h}^{−1}P_{h} − Δ^{−1})ω_{j}∥ & = ∥(R_{h} − I)Δ^{−1}ω_{j}∥ ≤ Ch^{2}∥Δ^{−1}ω_{j}∥_{2,2} ≤ Ch^{2}∥ω_{j}∥ \\ & ≤ Ch^{2} max_{t_{j−1}≤t≤t_{j}} ∥u_{t}(t)∥ ≤ Ch^{2}F (u). \end{aligned}

Hence

∥R_{5}∥ ≤ Ch^{2}F (u)k ∑_{j=1}^{[n/2]} t_{n−j+1}^{−1} ≤ Ch^{2}F (u).

For the sixth term we note that

∥Δ^{−1}ω_{j}∥ = ∥k^{−1} ∫_{t_{j−1}}^{t_{j}}(t − t_{j−1})Δ^{−1}u_{tt}(t) dt∥ ≤ k max_{t_{j−1}≤t≤t_{j}} ∥Δ^{−1}u_{tt}(t)∥ ≤ kG(u),

so that

∥R_{6}∥ ≤ CkG(u)k ∑_{j=1}^{[n/2]} t_{n−j+1}^{−1} ≤ CkG(u).

Finally, we have

∥ω_{j}∥ ≤ Ckt_{j}^{−1}G(u), t_{j} > 0,

because

∥ω_{j}∥ = ∥k^{−1} ∫_{t_{j−1}}^{t_{j}}(t − t_{j−1})u_{tt}(t) dt∥ ≤ kt_{j−1}^{−1} max_{t_{j−1}≤t≤t_{j}}(t∥u_{tt}(t)∥) ≤ Ckt_{j}^{−1}G(u),

for t_{j} ≥ t_{2}, and

∥ω_{1}∥ = ∥k^{−1} ∫_{0}^{k} tu_{tt}(t) dt∥ ≤ max_{0<t≤k}(t∥u_{tt}(t)∥) ≤ G(u) = kt_{1}^{−1}G(u).

Hence

∥R_{7}∥ ≤ Ck ∑_{j=[n/2]+1}^{n} ∥ω_{j}∥ ≤ CkG(u)k ∑_{j=[n/2]+1}^{n} t_{j}^{−1} ≤ CkG(u).

Taken together these estimates prove (2.8). □

3. Nonlinear error analysis

3.1 The semidiscrete case. Let u_{h}, ϕ_{h} be the semidiscrete finite element approximations of the solutions u, ϕ of the nonlinear problem (1.1). In this section we estimate the errors u(t) − u_{h}(t) and ϕ(t) − ϕ_{h}(t) uniformly over a finite time interval 0 ≤ t ≤ T under minimal assumptions about the regularity of u and ϕ. The error analysis is carried out assuming that the number of spatial variables d ≤ 3; the regularity assumptions, however, have only been verified for d ≤ 2, see Section 4 below. The result is presented in the following theorem.

Theorem 3.1. Let u, ϕ and u_{h}, ϕ_{h} be solutions of (1.4)–(1.5) and (1.6)–(1.7), respectively, with u_{h0} chosen so that

(3.1)  ∥u_{0} − u_{h0}∥ ≤ M_{1}h^{2}.

Assume further that d ≤ 3 and that

sup_{0<t≤T}(∥u(t)∥_{2,2} + t∥u_{t}(t)∥_{2,2}) ≤ M_{2}, (3.2)

sup_{0<t≤T}(∥g(t)∥_{H^{2}(∂Ω)} + ∥ϕ(t)∥_{2,2} + ∥ϕ(t)∥_{1,∞}) ≤ M_{3}, (3.3)

for some positive numbers T and M_{i}, i = 1, …, 3. Then there is a constant C = C(κ, K, M_{1}, M_{2}, M_{3}, T) such that

(3.4)  ∥u(t) − u_{h}(t)∥ + ∥ϕ(t) − ϕ_{h}(t)∥ ≤ Ch^{2},  t ∈ [0, T].

Here ∥ · ∥_{H^{2}(∂Ω)} is defined by summation over the flat parts of the polygon ∂Ω. In the remainder of this section we let C denote various quantities that may depend on the data of our problem as in the statement of Theorem 3.1. All estimates that are derived hold uniformly with respect to t ∈ [0, T]. We prepare for the proof of Theorem 3.1 by proving some preliminary bounds for ϕ(t) − ϕ_{h}(t).

Lemma 3.2. Under the assumptions of Theorem 3.1 we have

∥∇(ϕ(t) − ϕ_{h}(t))∥ ≤ C(h + ∥u(t) − u_{h}(t)∥), (3.5)

∥ϕ(t) − ϕ_{h}(t)∥ ≤ C(h^{2} + ∥u(t) − u_{h}(t)∥ + h^{−d/6}∥u(t) − u_{h}(t)∥^{2}). (3.6)

Proof. Since all results below are uniform in t we do not make the t-dependence explicit. Let e_{ϕ} = ϕ − ϕ_{h}. It follows from (1.5) and (1.7) that

(3.7)  (σ(u_{h})∇e_{ϕ}, ∇χ) = ([σ(u_{h}) − σ(u)]∇ϕ, ∇χ),  ∀χ ∈ S̊_{h},

and

(3.8)  (σ(u)∇e_{ϕ}, ∇χ) = ([σ(u_{h}) − σ(u)]∇ϕ_{h}, ∇χ),  ∀χ ∈ S̊_{h}.

Since π_{h}ϕ − ϕ_{h} ∈ S̊_{h}, we have from (3.7) that

\begin{aligned}(σ(u_{h})∇e_{ϕ}, ∇e_{ϕ}) & = (σ(u_{h})∇e_{ϕ}, ∇(ϕ − π_{h}ϕ)) + (σ(u_{h})∇e_{ϕ}, ∇(π_{h}ϕ − ϕ_{h})) \\ & = (σ(u_{h})∇e_{ϕ}, ∇(ϕ − π_{h}ϕ)) + ([σ(u_{h}) − σ(u)]∇ϕ, ∇(π_{h}ϕ − ϕ_{h})). \end{aligned}

It now follows from (1.3) and (3.3) that

∥∇e_{ϕ}∥^{2} ≤ C(∥∇e_{ϕ}∥ ∥∇(ϕ − π_{h}ϕ)∥ + ∥u_{h} − u∥ ∥ϕ∥_{1,∞}∥∇(ϕ_{h} − π_{h}ϕ)∥)
≤ C(∥∇e_{ϕ}∥ ∥∇(ϕ − π_{h}ϕ)∥ + ∥u_{h} − u∥(∥∇(ϕ − π_{h}ϕ)∥ + ∥∇e_{ϕ}∥))
≤ C(∥∇(ϕ − π_{h}ϕ)∥^{2} + ∥u_{h} − u∥^{2}) + \frac{1}{2}∥∇e_{ϕ}∥^{2}.

Hence

∥∇e_{ϕ}∥ ≤ C(∥∇(ϕ − π_{h}ϕ)∥ + ∥u_{h} − u∥),

which immediately yields (3.5) in view of an interpolation error estimate and (3.3).

The L_{2} estimate of e_{ϕ} is obtained by the standard duality argument. Let ψ be the unique solution of the Dirichlet problem

(3.9)  −∇ · (σ(u)∇ψ) = e_{ϕ} in Ω;  ψ = 0 on ∂Ω.

Then the standard regularity estimate yields (recall that Ω is a convex polygon)

∥ψ∥_{2,2} ≤ C∥Δψ∥ = C∥\frac{σ′(u)∇u · ∇ψ + e_{ϕ}}{σ(u)} ∥.

Using (1.3), Hölder’s inequality, an interpolation inequality (see (4.2) below), the trivial estimate ∥ψ∥_{1,2} ≤ C∥e_{ϕ}∥, and ∥u∥_{2,2} ≤ M_{2} from (3.2) we obtain

∥ψ∥_{2,2} ≤ C(∥∇u∥_{0,6} ∥∇ψ∥_{0,3} + ∥e_{ϕ}∥)
≤ C(∥u∥_{2,2} ∥ψ∥_{1,2}^{1−d/6} ∥ψ∥_{2,2}^{d/6} + ∥e_{ϕ}∥)
≤ C(∥u∥_{2,2}^{1/(1−d/6)}∥ψ∥_{1,2} + ∥e_{ϕ}∥) + \frac{1}{2}∥ψ∥_{2,2}
≤ C∥e_{ϕ}∥ + \frac{1}{2}∥ψ∥_{2,2}.

Hence

(3.10)  ∥ψ∥_{2,2} ≤ C∥e_{ϕ}∥.

From (3.9) and Green’s formula it follows that

\begin{aligned}∥e_{ϕ}∥^{2} & = −(∇ · (σ(u)∇ψ), e_{ϕ}) \\ & = (σ(u)∇ψ, ∇e_{ϕ}) − ⟨σ(u)∇ψ · ν, e_{ϕ}⟩ \\ & = (σ(u)∇(ψ − ψ_{h}), ∇e_{ϕ}) + (σ(u)∇ψ_{h}, ∇e_{ϕ}) \\ & − ⟨σ(u)∇ψ · ν, g − π_{h}g⟩ \\ & ≡ T_{1} + T_{2} + T_{3}, \end{aligned} (3.11)

where ψ_{h} = π_{h}ψ ∈ S̊_{h}, and where we have used ⟨·, ·⟩ to denote the inner product in L_{2}(∂Ω), and ν is the unit outward normal vector on ∂Ω.

The first term on the right of (3.11) is easily estimated by means of an interpolation error bound, (3.5) and (3.10):

\begin{aligned}|T_{1}| & ≤ C∥∇(ψ − ψ_{h})∥ ∥∇e_{ϕ}∥ ≤ Ch∥ψ∥_{2,2} C(h + ∥u − u_{h}∥) \\ & ≤ C(h^{2} + ∥u − u_{h}∥)∥e_{ϕ}∥. \end{aligned}

Using (3.8) we have

T_{2} = (σ(u)∇e_{ϕ}, ∇ψ_{h}) = ([σ(u_{h}) − σ(u)]∇ϕ_{h}, ∇ψ_{h})
= ([σ(u_{h}) − σ(u)]∇ϕ, ∇ψ_{h}) + ([σ(u_{h}) − σ(u)]∇(ϕ_{h} − ϕ), ∇ψ)
+ ([σ(u_{h}) − σ(u)]∇(ϕ_{h} − ϕ), ∇(ψ_{h} − ψ)).

Hence

|T_{2}| ≤ C∥u − u_{h}∥(∥ϕ∥_{1,∞}∥ψ_{h}∥_{1,2} + ∥∇e_{ϕ}∥_{0,3}∥ψ∥_{1,6}) + C∥∇e_{ϕ}∥ ∥∇(ψ_{h} − ψ)∥.

Using the facts that ∥ϕ∥_{1,∞} ≤ M_{3}, ∥ψ_{h}∥_{1,2}+∥ψ∥_{1,6} ≤ C∥ψ∥_{2,2} ≤ C∥e_{ϕ}∥ and ∥∇(ψ_{h}− ψ)∥ ≤ Ch∥ψ∥_{2,2} ≤ Ch∥e_{ϕ}∥ we obtain

|T_{2}| ≤ C(∥u − u_{h}∥(1 + ∥∇e_{ϕ}∥_{0,3})+h∥∇e_{ϕ}∥)∥e_{ϕ}∥.

By interpolation error estimates, an inverse estimate, (3.3) and (3.5) we have here

∥∇e_{ϕ}∥_{0,3} ≤ ∥∇(ϕ − π_{h}ϕ)∥_{0,3} + ∥∇(π_{h}ϕ − ϕ_{h})∥_{0,3}
≤ Ch^{1−d/6}∥ϕ∥_{2,2} + Ch^{−d/6}∥∇(π_{h}ϕ − ϕ_{h})∥
≤ Ch^{1−d/6}∥ϕ∥_{2,2} + Ch^{−d/6}(Ch∥ϕ∥_{2,2} + ∥∇e_{ϕ}∥)
≤ C(1 + h^{−d/6}∥u − u_{h}∥),

so that

|T_{2}| ≤ C(h^{2} + ∥u − u_{h}∥ + h^{−d/6}∥u − u_{h}∥^{2})∥e_{ϕ}∥.

Finally we have

|T_{3}| ≤ C∥g − π_{h}g∥_{L_{2}(∂Ω)}∥ψ∥_{2,2} ≤ Ch^{2}∥g∥_{H^{2}(∂Ω)}∥e_{ϕ}∥ ≤ Ch^{2}∥e_{ϕ}∥,

where a trace inequality, an interpolation error estimate and (3.10) have been used. Together the above estimates prove (3.6). □

Proof of Theorem 3.1. It is convenient to split the error into two parts: u_{h} − u = (u_{h} − ũ_{h}) + (ũ_{h} − u), where ũ_{h} : [0, T] → S̊_{h} is uniquely defined by

(3.12)  \begin{matrix}(ũ_{h,t}, χ) + (∇ũ_{h}, ∇χ) = (F (u, ϕ), χ),  ∀χ ∈ S̊_{h}, t ∈ [0, T], \\ ũ_{h}(0) = u_{h0}, \end{matrix}

with F (u, ϕ) = σ(u)|∇ϕ|^{2}. Applying the known error analysis for linear parabolic equations we obtain

(3.13)  ∥ũ_{h}(t) − u(t)∥ ≤ Ch^{2},

where C depends on M_{1} and M_{2}, see Theorem 2.1.

Forming the difference between (1.6) and (3.12) we have for ζ = u_{h} − ũ_{h} that

ζ_{t} − Δ_{h}ζ = P_{h}(F (u_{h}, ϕ_{h}) − F (u, ϕ)), t ∈ [0, T]; ζ(0) = 0.

Hence the variation of constants formula implies that

(3.14)  ∥ζ(t)∥ ≤ ∫_{0}^{t}∥E_{h}(t − s)P_{h}(F (u_{h}(s), ϕ_{h}(s)) − F (u(s), ϕ(s)))∥ ds.

We proceed to estimate ∥ζ(t)∥ by bounding the right side in various ways. In doing so we shall need several bounds for the operator E_{h}(t)P_{h}. In addition to (2.3) we quote from [11, Lemma 5.2] a bound of the norm of E_{h}(t)P_{h} considered as an operator from L_{2} into L_{∞}, namely, for any ϵ > 0 there is C_{ϵ} > 0 such that

(3.15)  ∥E_{h}(t)P_{h}ψ∥_{0,∞} ≤ C_{ϵ}t^{−d/4−ϵ}∥ψ∥,  t > 0.

By duality we also have the same bound for the norm of E_{h}(t)P_{h} : L_{1} → L_{2}, that is,

(3.16)  ∥E_{h}(t)P_{h}ψ∥ ≤ C_{ϵ}t^{−d/4−ϵ}∥ψ∥_{0,1},  t > 0, ϵ > 0.

In fact,

(3.17)  ∥E_{h}(t)P_{h}ψ∥ = sup_{χ∈L_{2}} \frac{|(E_{h}(t)P_{h}ψ, χ)|}{∥χ∥} = sup_{χ∈L_{2}} \frac{|(ψ, E_{h}(t)P_{h}χ)|}{∥χ∥},

since E_{h}(t)P_{h} is selfadjoint, so that (3.16) follows from (3.15).

We begin by deriving a preliminary low order estimate of ∥ζ(t)∥. We have

(3.18) \begin{matrix}∥F (u_{h}, ϕ_{h}) − F (u, ϕ)∥_{0,1} \\ ≤ ∥σ(u_{h})∇(ϕ_{h} + ϕ) · ∇(ϕ_{h} − ϕ)∥_{0,1} + ∥(σ(u_{h}) − σ(u))|∇ϕ|^{2}∥_{0,1} \\ ≤ C(∥∇ϕ_{h}∥ + ∥∇ϕ∥)∥∇(ϕ_{h} − ϕ)∥ + C∥u_{h} − u∥ ∥ϕ∥_{1,∞}^{2} \\ ≤ C(h + ∥u_{h} − u∥), \end{matrix}

where we have used the easily proved fact that ∥∇ϕ_{h}(t)∥+∥∇ϕ(t)∥ ≤ C, the assumption ∥ϕ(t)∥_{1,∞} ≤ M_{3}, and the error bound (3.5). Hence by (3.14), (3.16) and (3.13) we have

∥u_{h}(t) − u(t)∥ ≤ ∥ũ_{h}(t) − u(t)∥ + ∥ζ(t)∥
≤ Ch^{2} + C ∫_{0}^{t}(t − s)^{−α}∥F (u_{h}(s), ϕ_{h}(s)) − F (u(s), ϕ(s))∥_{0,1} ds
≤ Ch^{2} + C ∫_{0}^{t}(t − s)^{−α}(h + ∥u_{h}(s) − u(s)∥) ds
≤ Ch + C ∫_{0}^{t}(t − s)^{−α}∥u_{h}(s) − u(s)∥ ds,

where it is possible to choose α ∈ (3/4, 1), since d ≤ 3. Hence a variant of Gronwall’s lemma (see, e.g., [14, Lemma 5.6.7], [10, Lemma 7.1.1] or Lemma 3.4 below) yields the preliminary bound

∥u_{h}(t) − u(t)∥ ≤ Ch.

Inserted into (3.5) and (3.6) this gives

(3.19)  \begin{matrix}∥ϕ_{h}(t) − ϕ(t)∥_{1,2} ≤ Ch, \\ ∥ϕ_{h}(t) − ϕ(t)∥ ≤ C(h^{2} + ∥u_{h}(t) − u(t)∥). \end{matrix}

The reason for the suboptimality of the preliminary bound is that we estimated F (u_{h}, ϕ_{h}) − F (u, ϕ) in terms of ∇(ϕ_{h} − ϕ), which is only O(h). In order to obtain an estimate of second order we shall use a duality argument to remove the gradient from the latter term. This argument requires a more accurate expansion of F (u_{h}, ϕ_{h}) − F (u, ϕ), namely,

\begin{aligned}F (u_{h}, ϕ_{h}) − F (u, ϕ) & = [σ(u_{h}) − σ(u)]|∇ϕ|^{2} \\ & + 2σ(u)∇ϕ · ∇(ϕ_{h} − ϕ) \\ & + 2[σ(u_{h}) − σ(u)]∇ϕ · ∇(ϕ_{h} − ϕ) \\ & + σ(u_{h})|∇(ϕ_{h} − ϕ)|^{2} \\ & ≡ R_{1} + R_{2} + R_{3} + R_{4}. \end{aligned} (3.20)

Using (2.3) and (3.16) we shall estimate each of the terms ∥E_{h}(t − s)P_{h}R_{i}(s)∥ and substitute the result into the right hand side of (3.14).

Omitting the dependence on t − s and s in intermediate steps we obtain for the first term

\begin{aligned}∥E_{h}(t − s)P_{h}R_{1}(s)∥ & = ∥E_{h}P_{h}([σ(u_{h}) − σ(u)]|∇ϕ|^{2})∥ ≤ C∥u_{h} − u∥ ∥ϕ∥_{1,∞}^{2} \\ & ≤ C∥u_{h}(s) − u(s)∥. \end{aligned}

For the second term we use a duality argument (cf. (3.17)): for χ ∈ L_{2} we have

(E_{h}(t − s)P_{h}R_{2}(s), χ) = 2(∇(ϕ_{h} − ϕ), σ(u)∇ϕE_{h}P_{h}χ)
= −2(ϕ_{h} − ϕ, ∇ · [σ(u)∇ϕE_{h}P_{h}χ])
= −2(ϕ_{h} − ϕ, σ(u)∇ϕ · ∇[E_{h}P_{h}χ]),

since ∇ · (σ(u)∇ϕ) = 0. Hence, by (2.3),

\begin{aligned}|(E_{h}(t − s)P_{h}R_{2}(s), χ)| & ≤ 2∥ϕ_{h} − ϕ∥ ∥ϕ∥_{1,∞} ∥E_{h}P_{h}χ∥_{1,2} \\ & ≤ C(t − s)^{−1/2}∥ϕ_{h}(s) − ϕ(s)∥ ∥χ∥, \end{aligned}

so that, in view of the second estimate in (3.19),

∥E_{h}(t − s)P_{h}R_{2}(s)∥ ≤ C(t − s)^{−1/2}(h^{2} + ∥u_{h}(s) − u(s)∥).

By (3.16) there is α ∈ (3/4, 1) such that

∥E_{h}(t − s)P_{h}R_{3}(s)∥ ≤ C(t − s)^{−α}∥[σ(u_{h}) − σ(u)]∇ϕ · ∇(ϕ_{h} − ϕ)∥_{0,1}
≤ C(t − s)^{−α}∥u_{h} − u∥ ∥ϕ∥_{1,∞}(∥∇ϕ_{h}∥ + ∥∇ϕ∥)
≤ C(t − s)^{−α}∥u_{h}(s) − u(s)∥.

Similarly, by the first estimate in (3.19),

\begin{aligned}∥E_{h}(t − s)P_{h}R_{4}(s)∥ & ≤ C(t − s)^{−α}∥σ(u_{h})|∇(ϕ_{h} − ϕ)|^{2}∥_{0,1} \\ & ≤ C(t − s)^{−α}∥ϕ_{h}(s) − ϕ(s)∥_{1,2}^{2} ≤ C(t − s)^{−α}h^{2}. \end{aligned}

Combining these bounds with (3.14) and (3.13) we obtain

(3.21)  ∥u_{h}(t) − u(t)∥ ≤ Ch^{2} + C ∫_{0}^{t}(t − s)^{−α}(h^{2} + ∥u_{h}(s) − u(s)∥) ds.

Gronwall’s lemma now yields the desired bound for ∥u_{h}(t)−u(t)∥ in (3.4), and hence, in view of (3.19), also the bound for ∥ϕ_{h}(t) − ϕ(t)∥. □

3.2 The completely discrete case. We now turn to the completely discrete scheme. Our result is the following.

Theorem 3.3. Let u, ϕ and U_{n}, Φ_{n} be solutions of (1.4)–(1.5) and (1.8)–(1.9), respectively, with u_{h0} chosen so that

∥u_{0} − u_{h0}∥ ≤ M_{1}h^{2}.

Assume further that d ≤ 3 and that

sup_{0<t≤T}(∥u(t)∥_{2,2} + ∥u_{t}(t)∥ + ∥Δ^{−1}u_{tt}(t)∥ + t∥u_{t}(t)∥_{2,2} + t∥u_{tt}(t)∥) ≤ M_{2},

sup_{0<t≤T}(∥g(t)∥_{H^{2}(∂Ω)} + ∥ϕ(t)∥_{2,2} + ∥ϕ(t)∥_{1,∞} + ∥ϕ_{t}(t)∥_{1,2}) ≤ M_{3},

and that k ≤ M_{4}h^{d/6} for some positive numbers T and M_{i}, i = 1, …, 4. Then there is C = C(κ, K, M_{1}, M_{2}, M_{3}, M_{4}, T) such that

∥u(t_{n}) − U_{n}∥ + ∥ϕ(t_{n}) − Φ_{n}∥ ≤ C(h^{2} + k), t_{n} ∈ [0, T].

We will need a discrete version of the generalized Gronwall lemma that we referred to in the previous proof. We formulate this in the following lemma, where we use the convention that a sum is considered to be empty if its upper limit is smaller than its lower limit, that is, ∑_{l=n}^{m} a_{l} = 0 if m < n.

Lemma 3.4. Assume that the sequence φ_{n} satisfies

0 ≤ φ_{n} ≤ A + Bk ∑_{l=0}^{n−1} t_{n−l}^{−1+β}φ_{l}, for t_{n} ∈ [0, T],

where A, B, T are positive numbers and 0 < β < 1. Then there is a constant C = C(β, B, T) such that φ_{n} ≤ CA for t_{n} ∈ [0, T].

Proof. Iterating the given inequality once, using the inequalities

k ∑_{l=0}^{n−1} t_{n−l}^{−1+β} ≤ ∫_{0}^{t_{n}}(t_{n} − s)^{−1+β} ds = C_{β}t_{n}^{β},

k ∑_{l=j+1}^{n−1} t_{n−l}^{−1+β}t_{l−j}^{−1+γ} ≤ C_{β} ∫_{t_{j}}^{t_{n}}(t_{n} − s)^{−1+β}(s − t_{j})^{−1+γ} ds = C_{β,γ}t_{n−j}^{−1+β+γ},

valid for 0 < β, γ < 1, we obtain

\begin{aligned}φ_{n} & ≤ A + ABk ∑_{l=0}^{n−1} t_{n−l}^{−1+β} + B^{2}k ∑_{j=0}^{n−2}(k ∑_{l=j+1}^{n−1} t_{n−l}^{−1+β}t_{l−j}^{−1+β})φ_{j} \\ & ≤ C_{1}(β, B, T)A + C_{2}(β, B)k ∑_{l=0}^{n−2} t_{n−l}^{−1+2β}φ_{l}. \end{aligned}

After N − 1 iterations, where N is the smallest integer such that −1 + Nβ ≥ 0, we have

\begin{aligned}φ_{n} & ≤ C_{1}(β, B, T)A + C_{2}(β, B)k ∑_{l=0}^{n−N} t_{n−l}^{−1+Nβ}φ_{l} \\ & ≤ C_{1}(β, B, T)A + C_{2}(β, B, T)k ∑_{l=0}^{n−N} φ_{l}, \end{aligned}

and the desired conclusion follows by the standard Gronwall lemma. □

Proof of Theorem 3.3. The proof is a generalization of the proof Theorem 3.1. We begin by splitting the error into two parts: U_{n} − u_{n} = (U_{n} − Ũ_{n}) + (Ũ_{n} − u_{n}), where Ũ_{n} ∈ S̊_{h} is uniquely defined by

(3.22)  \begin{matrix}(∂_{n}Ũ_{n}, χ) + (∇Ũ_{n}, ∇χ) = (F (u_{n}, ϕ_{n}), χ),  ∀χ ∈ S̊_{h}, t_{n} ∈ [0, T], \\ Ũ_{0} = u_{h0}, \end{matrix}

with u_{n} = u(t_{n}), ϕ_{n} = ϕ(t_{n}) and F (u_{n}, ϕ_{n}) = σ(u_{n})|∇ϕ_{n}|^{2}. Applying the known error analysis for linear parabolic equations we obtain

(3.23)  ∥Ũ_{n} − u_{n}∥ ≤ C(h^{2} + k),

where C depends on M_{1} and M_{2}, see Theorem 2.1.

Forming the difference between (1.8) and (3.22) and applying the variation of constants formula we have for ζ_{n} = U_{n} − Ũ_{n} that

(3.24)  ∥ζ_{n}∥ ≤ k ∑_{l=0}^{n−1}∥E_{kh}^{n−l}P_{h}((F (U_{l}, Φ_{l}) − F (u_{l+1}, ϕ_{l+1}))∥.

In addition to (2.4) the discrete evolution operator E_{kh}^{n}P_{h} satisfies the bound

(3.25)  ∥E_{kh}^{n}P_{h}ψ∥ ≤ C_{ϵ}t_{n}^{−d/4−ϵ}∥ψ∥_{0,1},  t_{n} > 0, ϵ > 0.

The proof of this is analogous to that of (3.16). Indeed, by inspection of the proof of [11, Lemma 5.2] it is clear that a completely discrete analog of (3.15) holds and hence (3.25) follows by duality.

Lemma 3.2 is directly applicable to the equation for Φ_{n} and gives

(3.26)  \begin{matrix}∥∇(Φ_{n} − ϕ_{n})∥ ≤ C(h + ∥U_{n} − u_{n}∥), \\ ∥Φ_{n} − ϕ_{n}∥ ≤ C(h^{2} + ∥U_{n} − u_{n}∥ + h^{−d/6}∥U_{n} − u_{n}∥^{2}). \end{matrix}

Hence,

(3.27)  ∥F (U_{n−1}, Φ_{n−1}) − F (u_{n−1}, ϕ_{n−1})∥_{0,1} ≤ C(h + ∥U_{n−1} − u_{n−1}∥),

cf. (3.18). In a similar way we have

\begin{matrix}∥F (u_{n−1}, ϕ_{n−1}) − F (u_{n}, ϕ_{n})∥_{0,1} \\ ≤ ∥σ(u_{n−1})∇(ϕ_{n−1} + ϕ_{n}) · ∇(ϕ_{n−1} − ϕ_{n})∥_{0,1} \\ + ∥(σ(u_{n−1}) − σ(u_{n}))|∇ϕ_{n}|^{2}∥_{0,1} \\ ≤ C(∥∇ϕ_{n−1}∥ + ∥∇ϕ_{n}∥)∥ϕ_{n−1} − ϕ_{n}∥_{1,2} \\ + C∥u_{n−1} − u_{n}∥ ∥ϕ_{n}∥_{1,∞}^{2} \\ ≤ C(∥ϕ_{n−1} − ϕ_{n}∥_{1,2} + ∥u_{n−1} − u_{n}∥) \\ ≤ Ck sup_{0≤t≤T}(∥ϕ_{t}(t)∥_{1,2} + ∥u_{t}(t)∥) ≤ Ck. \end{matrix} (3.28)

Therefore

∥F (U_{n−1}, Φ_{n−1}) − F (u_{n}, ϕ_{n})∥_{0,1} ≤ C(h + k + ∥U_{n−1} − u_{n−1}∥),

so that, by (3.24), (3.23) and (3.25),

∥U_{n} − u_{n}∥ ≤ ∥Ũ_{n} − u_{n}∥ + ∥ζ_{n}∥

≤ C(h^{2} + k) + Ck ∑_{l=0}^{n−1} t_{n−l}^{−α}∥F (U_{l}, Φ_{l}) − F (u_{l+1}, ϕ_{l+1})∥_{0,1}

≤ C(h^{2} + k) + Ck ∑_{l=0}^{n−1} t_{n−l}^{−α}(h + k + ∥U_{l} − u_{l}∥)

≤ C(h + k) + Ck ∑_{l=0}^{n−1} t_{n−l}^{−α}∥U_{l} − u_{l}∥,

where we have chosen α ∈ (3/4, 1). Hence Lemma 3.4 yields the preliminary bound

∥U_{n} − u_{n}∥ ≤ C(h + k),

which, in view of (3.26), leads to

(3.29)  \begin{matrix}∥Φ_{n} − ϕ_{n}∥_{1,2} ≤ C(h + k), \\ ∥Φ_{n} − ϕ_{n}∥ ≤ C(h + k + ∥U_{n} − u_{n}∥), \end{matrix}

where we have also used the assumption that kh^{−d/6} ≤ M_{4}. In order to complete the proof we repeat the steps leading to (3.21) using (3.28) and replacing (3.27) by a more accurate expansion as in (3.20). This gives

∥U_{n} − u_{n}∥ ≤ C(h^{2} + k) + Ck ∑_{l=0}^{n−1} t_{n−l}^{−α}∥U_{l} − u_{l}∥,

so that ∥U_{n} − u_{n}∥ ≤ C(h^{2} + k) follows by the Gronwall argument of Lemma 3.4, and hence, in view of (3.26), we also obtain ∥Φ_{n} − ϕ_{n}∥ ≤ C(h^{2} + k). □

4. Existence and regularity

In this section we study the solvability of the system (1.1)–(1.2). The existence of global weak solutions in two space dimensions (Ω ⊂ R^{2}) was shown by Cimatti [5], see also Rodrigues [15] and Allegretto and Xie [2] for existence results for related problems. The regularity of these solutions, however, is insufficient for the purpose of proving error bounds of optimal order, cf. Theorems 3.1 and 3.3. Building upon the techniques of [5] we obtain global strong solutions with the required regularity in one and two space dimensions. The three-dimensional case remains an open problem.

Theorem 4.1. Assume that Ω ⊂ R^{d}, d ≤ 2, is a bounded domain whose boundary is either smooth or a convex polygon. Let T > 0, r > 2 and assume that u_{0} ∈ H^{2} ∩ H_{0}^{1}, g ∈ L_{∞}([0, T], W_{r}^{2}), g_{t} ∈ L_{∞}([0, T], H^{2}) and g_{tt} ∈ L_{∞}([0, T], H^{1}). Then (1.1)–(1.2) has a unique solution u ∈ C^{1}([0, T], L_{2})∩C([0, T], H^{2}), ϕ ∈ L_{∞}([0, T], H^{2}). Moreover, there is a constant C, depending on T, r, u_{0}, g, Ω and on σ through the constants κ, K in (1.3), such that for t ∈ [0, T] we have

∥u(t)∥_{2,2} + ∥u_{t}(t)∥ + ∥Δ^{−1}u_{tt}(t)∥ + t∥u_{t}(t)∥_{2,2} + t∥u_{tt}(t)∥ + ∥ϕ(t)∥_{2,2} + ∥ϕ(t)∥_{1,∞} + ∥ϕ_{t}(t)∥_{1,2} ≤ C.

In order to prepare the way for the proof we recall some facts that we shall need. The assumption about Ω guarantees that for any p ∈ [2, ∞) there is C such that

(4.1)  ∥u∥_{2,p} ≤ C∥Δu∥_{0,p},  ∀u ∈ W_{p}^{2} ∩ H_{0}^{1},

see [9]. Under even weaker assumptions about Ω ⊂ R^{d} we have the following interpolation result [1]: let 1 ≤ p ≤ ∞, m ≥ 1 and assume that u ∈ W_{p}^{m}. Then there is a constant C = C(m, p, q, d, Ω) such that the inequality

(4.2)  ∥u∥_{0,q} ≤ C∥u∥_{0,p}^{1−θ}∥u∥_{m,p}^{θ},  where θ = \frac{d}{m}(\frac{1}{p} − \frac{1}{q}),

holds for q ∈ [p, ∞] if m − d/p > 0, for q ∈ [p, ∞) if m − d/p = 0, and for q ∈ [p, −d/(m − d/p)] if m − d/p < 0. Note that 0 ≤ θ ≤ 1. We shall also use a theorem of Meyers [13], which we quote here in a special case suitable for our purpose. We use the standard notation W̊_{q}^{1} = {u ∈ W_{q}^{1} : u|_{∂Ω} = 0} and W_{q}^{−1} is the dual space of W̊_{q′}^{1}, where q and q′ are conjugate exponents.

Theorem 4.2. (Meyers [13]) Assume that Ω ⊂ R^{d} has the property that, for some q ∈ (2, ∞) and L ≥ 1, the Laplacian Δ is an isomorphism from W̊_{q}^{1} onto W_{q}^{−1} with ∥Δ^{−1}∥ ≤ L, and let the function a satisfy the inequalities 0 < κ ≤ a(x) ≤ K for all x ∈ Ω. Then there are p ∈ (2, q) and C > 0 depending only on q, κ, K, L such that the following holds true: let f ∈ (L_{p})^{d} be a vector field and let u ∈ H_{0}^{1} be the unique solution of

(a∇u, ∇χ) = (f, ∇χ), ∀χ ∈ H_{0}^{1}.

Then u ∈ W̊_{p}^{1} and ∥∇u∥_{0,p} ≤ C∥f∥_{0,p}.

The assumption in Meyers’ theorem is satisfied, for example, if Ω ⊂ R^{2} is bounded and ∂Ω is either smooth [12] or a polygon [8]. See also [3] for a modern presentation of Meyers’ theorem, and [16] for a finite element version.

Proof of Theorem 4.1. Let V_{m} be the eigenspace corresponding to the m smallest eigenvalues of the operator −Δ with domain of definition H^{2} ∩ H_{0}^{1}. We consider the initial value problem

(4.3)  \begin{matrix}U(t) ∈ V_{m}, \\ (U_{t}, χ) + (∇U, ∇χ) = (σ(U)|∇Φ|^{2}, χ), & ∀χ ∈ V_{m}, t > 0, \\ (U(0), χ) = (u_{0}, χ), & ∀χ ∈ V_{m}, \end{matrix}

where Φ(t) is determined by the linear elliptic boundary value problem

(4.4)  \begin{matrix}Φ(t) ∈ H^{1}, Φ(t) − g(t) ∈ H_{0}^{1}, \\ (σ(U)∇Φ, ∇χ) = 0,  ∀χ ∈ H_{0}^{1}, t > 0. \end{matrix}

Clearly, given U(t) and g(t) there is a unique solution Φ(t) of (4.4). It follows that (4.3) is an initial value problem for a finite dimensional system of ordinary differential equations for the Fourier coefficients of U. Hence there exists a unique solution to (4.3) on a time interval [0, T_{m}], where T_{m} depends on u_{0} and m. We proceed to derive a priori bounds, which show that (4.3) has a solution on the prescribed time interval [0, T]. These estimates will also allow the passage to the limit in U and Φ as m → ∞, yielding the existence of a solution u, ϕ to (1.1)–(1.2) with the desired regularity. This passage to the limit is rather standard and we omit the details (cf. [5]).

Throughout this proof we let C denote various quantities that may depend on the data T, r, u_{0}, g, Ω and on σ through the constants κ, K in (1.3), but not on m and t. All estimates that we derive below hold uniformly with respect to t ∈ [0, T].

Step 1. We begin by showing some preliminary estimates of Φ. The starting point is the maximum principle, which yields

(4.5)  ∥Φ(t)∥_{0,∞} ≤ ∥g(t)∥_{0,∞} ≤ C.

Next we apply Meyers’ theorem 4.2 to equation (4.4), which implies that there is p > 2 such that ∥∇Φ∥_{0,p} ≤ C∥∇g∥_{0,p}. The constant C depends only on Ω and σ through the bounds in (1.3). The optimal value of p is unknown; for simplicity we assume that 2 < p ≤ r. Together with (4.5) this shows that

(4.6)  ∥Φ(t)∥_{1,p} ≤ C.

Further estimates of Φ depend on derivatives of σ(U), and we shall take this carefully into account.

First we note that equation (4.4) implies −σ(U)ΔΦ − ∇σ(U) · ∇Φ = 0, so that, by Hölder’s inequality,

∥ΔΦ∥_{0,p} = ∥\frac{σ′(U)}{σ(U)}∇U · ∇Φ∥_{0,p} ≤ C∥∇U∥_{0,q′}∥∇Φ∥_{0,q},

for any q, q′ satisfying 1/q + 1/q′ = 1/p. Applying the regularity estimate (4.1) to the function Φ − g ∈ H_{0}^{1}, we thus obtain

∥Φ∥_{2,p} ≤ C(∥ΔΦ∥_{0,p} + ∥g∥_{2,p}) ≤ C(1 + ∥∇U∥_{0,q′}∥∇Φ∥_{0,q}).

Hence using also the interpolation inequality (4.2) and (4.6) we have

∥Φ∥_{2,p} ≤ C(1 + ∥U∥_{2,2}∥Φ∥_{1,p}^{α}∥Φ∥_{2,p}^{1−α})
≤ C(1 + ∥U∥_{2,2}∥Φ∥_{2,p}^{1−α})
≤ C(1 + ∥U∥_{2,2}^{1/α}) + \frac{1}{2}∥Φ∥_{2,p},

where α = 1 − d/p + d/q = 1 − d/q′. In the last step we also used Young’s inequality

(4.7)  ab ≤ ϵ^{1−1/θ}a^{1/θ} + ϵb^{1/(1−θ)},  ϵ > 0, 0 < θ < 1, a, b ≥ 0.

For the above estimate of ∥∇U∥_{0,q′} to hold it is required that q′ < ∞, which in its turn is equivalent to α < 1. We have thus proved the preliminary estimate

(4.8)  ∥Φ(t)∥_{2,p} ≤ C(1 + ∥U(t)∥_{2,2}^{1/α}),  for α ∈ [1 − d/p, 1),

where C is independent of α. We next proceed to show that there are β < 1 and C > 0 such that

(4.9)  ∥Φ(t)∥_{1,∞} + ∥Φ(t)∥_{1,4}^{2} ≤ C(1 + ∥U(t)∥_{2,2}^{β}).

In fact, arguing as above using (4.6) and (4.8) we have

∥∇Φ∥_{0,∞} ≤ C∥Φ∥_{1,p}^{1−γ}∥Φ∥_{2,p}^{γ} ≤ C(1 + ∥U∥_{2,2}^{γ/α}),

where γ = d/p < 1. By taking α sufficiently near 1, and in view of the maximum norm estimate in (4.5), we obtain the desired sublinear estimate of ∥Φ∥_{1,∞} in (4.9). Similarly,

∥∇Φ∥_{0,4}^{2} ≤ (C∥Φ∥_{1,p}^{1−γ}∥Φ∥_{2,p}^{γ})^{2} ≤ C(1 + ∥U∥_{2,2}^{2γ/α}),

where now γ = d/p − d/4. Since 2γ = (4/p − 1)d/2 < 1, the bound for ∥Φ∥_{1,4}^{2} in (4.9) follows by taking α sufficiently near 1.

Remark. This is where the restriction to two space variables occurs: if d = 3, then we must have q′ ≤ 6 and α ≤ 1/2, so that we can only guarantee that β < 3 in (4.9).

Step 2. We now estimate ∥U∥_{2,2} and ∥U_{t}∥. We begin by noting that it suffices to estimate ∥U_{t}∥. Indeed, equation (4.3) implies that U_{t} − ΔU = P_{m}(σ(U)|∇Φ|^{2}), where P_{m} denotes the orthogonal projection onto V_{m}. Hence using the regularity estimate (4.1), (4.9) and (4.7), we obtain

∥U∥_{2,2} ≤ C∥ΔU∥ ≤ C(∥U_{t}∥ + ∥σ(U)|∇Φ|^{2}∥)
≤ C(∥U_{t}∥ + ∥Φ∥_{1,4}^{2}) ≤ C(1 + ∥U_{t}∥ + ∥U∥_{2,2}^{β})
≤ C(1 + ∥U_{t}∥) + \frac{1}{2}∥U∥_{2,2},

since β < 1. This shows that

(4.10)  ∥U(t)∥_{2,2} ≤ C(1 + ∥U_{t}(t)∥).

Taking χ = U_{t} in (4.3) we obtain

\begin{aligned}∥U_{t}∥^{2} + \frac{1}{2} \frac{d}{dt}∥∇U∥^{2} & = (σ(U)|∇Φ|^{2}, U_{t}) ≤ C∥Φ∥_{1,4}^{2}∥U_{t}∥ \\ & ≤ C(1 + ∥U_{t}∥^{β})∥U_{t}∥ ≤ C + \frac{1}{2}∥U_{t}∥^{2}, \end{aligned}

since β < 1, where we have employed (4.9), (4.10) and (4.7). Integration with respect to t then yields

(4.11)  ∫_{0}^{t} ∥U_{t}∥^{2} ds + ∥∇U(t)∥^{2} ≤ C∥∇(P_{m}u_{0})∥^{2} + Ct ≤ C,

since U(0) = P_{m}u_{0} where P_{m} is bounded independently of m with respect to the norm ∥∇ · ∥.

In order to obtain further estimates of U_{t} we differentiate equations (4.3) and (4.4) with respect to t. Beginning with (4.4) we have

(4.12)  (σ(U)∇Φ_{t}, ∇χ) = −(σ(U)_{t}∇Φ, ∇χ),  ∀χ ∈ H_{0}^{1}.

Straightforward estimates based on taking χ = Φ_{t} − g_{t} give

∥Φ_{t}∥_{1,2} ≤ C(∥g_{t}∥_{1,2} + ∥σ′(U)U_{t}∇Φ∥) ≤ C(1 + ∥U_{t}∥ ∥Φ∥_{1,∞}),

so that in view of (4.9) and (4.10)

(4.13)  ∥Φ_{t}(t)∥_{1,2} ≤ C(1 + ∥U_{t}(t)∥^{2}).

Next we note that the source term in (4.3) may be transformed as follows using Green’s formula and equation (4.4):

(σ(U)|∇Φ|^{2}, χ) = −(σ(U)Φ∇Φ, ∇χ), ∀χ ∈ V_{m}.

Differentiating equation (4.3) with respect to t we thus have

(4.14)  (U_{tt}, χ) + (∇U_{t}, ∇χ) = −((σ(U)Φ∇Φ)_{t}, ∇χ),  ∀χ ∈ V_{m}.

With χ = U_{t} this leads to

(4.15)  \frac{d}{dt}∥U_{t}∥^{2} + ∥∇U_{t}∥^{2} ≤ ∥(σ(U)Φ∇Φ)_{t}∥^{2}.

Here we have

∥(σ(U)Φ∇Φ)_{t}∥ ≤ ∥σ′(U)U_{t}Φ∇Φ∥ + ∥σ(U)Φ_{t}∇Φ∥ + ∥σ(U)Φ∇Φ_{t}∥
≤ C(∥U_{t}∥ ∥Φ∥_{0,∞}∥Φ∥_{1,∞} + ∥Φ_{t}∥_{0,q}∥Φ∥_{1,p} + ∥Φ∥_{0,∞}∥Φ_{t}∥_{1,2}),

where 1/q + 1/p = 1/2. Using Sobolev’s inequality ∥Φ_{t}∥_{0,q} ≤ C∥Φ_{t}∥_{1,2} and known bounds for Φ and Φ_{t} in (4.5), (4.9), (4.13) and (4.10), we arrive at

∥(σ(U)Φ∇Φ)_{t}∥ ≤ C(1 + ∥U_{t}∥^{2}).

Using the fact that U(0) = P_{m}u_{0}, so that ∥U(0)∥_{2,2} ≤ C∥ΔU(0)∥ ≤ C∥Δu_{0}∥ ≤ C, and hence by (4.9),

\begin{aligned}∥U_{t}(0)∥ & ≤ ∥ΔU(0)∥ + ∥P_{m}(σ(U(0))|∇Φ(0)|^{2})∥ \\ & ≤ ∥ΔU(0)∥ + C∥Φ(0)∥_{1,4}^{2} ≤ C(1 + ∥U(0)∥_{2,2}) ≤ C, \end{aligned}

we integrate (4.15) to get

∥U_{t}(t)∥^{2} + ∫_{0}^{t} ∥∇U_{t}∥^{2} ds ≤ C + C ∫_{0}^{t} ∥U_{t}∥^{4} ds.

Applying Gronwall’s lemma together with (4.11) we obtain

(4.16)  ∥U_{t}(t)∥^{2} + ∫_{0}^{t} ∥∇U_{t}∥^{2} ds ≤ C exp (∫_{0}^{t} ∥U_{t}∥^{2} ds) ≤ C.

Substituting this result into (4.10), (4.9), (4.13) and (4.8) we may conclude

(4.17)  ∥U(t)∥_{2,2} + ∥U_{t}(t)∥ + ∥Φ(t)∥_{2,2} + ∥Φ(t)∥_{1,∞} + ∥Φ_{t}(t)∥_{1,2} ≤ C.

Step 3. It remains to bound ∥Δ^{−1}U_{tt}(t)∥, t∥U_{t}(t)∥_{2,2} and t∥U_{tt}(t)∥. We begin by noting that U_{tt} − ΔU_{t} = P_{m}(σ(U)|∇Φ|^{2})_{t}, where, in view of (4.17),

\begin{aligned}∥(σ(U)|∇Φ|^{2})_{t}∥ & ≤ ∥σ′(U)U_{t}|∇Φ|^{2}∥ + 2∥σ(U)∇Φ · ∇Φ_{t}∥ \\ & ≤ C(∥U_{t}∥ ∥Φ∥_{1,∞}^{2} + ∥Φ∥_{1,∞}∥Φ_{t}∥_{1,2}) ≤ C, \end{aligned} (4.18)

so that

(4.19)  ∥U_{t}(t)∥_{2,2} ≤ C(1 + ∥U_{tt}(t)∥),

and also

∥Δ^{−1}U_{tt}(t)∥ ≤ ∥U_{t}∥ + ∥Δ^{−1}(σ(U)|∇Φ|^{2})_{t}∥ ≤ C.

It now only remains to estimate t∥U_{tt}(t)∥. In order to do so we differentiate equation (4.3) with respect to t and substitute χ = U_{tt}, which after some simple manipulations gives

∥U_{tt}∥^{2} + \frac{d}{dt}∥∇U_{t}∥^{2} ≤ ∥(σ(U)|∇Φ|^{2})_{t}∥^{2} ≤ C,

where we employed (4.18) in the last step. Multiplication by t and integration now yields

(4.20)  ∫_{0}^{t} s∥U_{tt}∥^{2} ds + t∥∇U_{t}(t)∥^{2} ≤ C + ∫_{0}^{t} ∥∇U_{t}∥^{2} ds ≤ C,

in view of (4.16). In the next step we differentiate equation (4.14) with respect to t and substitute χ = U_{tt} to obtain

(4.21)  \frac{d}{dt}∥U_{tt}∥^{2} + ∥∇U_{tt}∥^{2} ≤ ∥(σ(U)Φ∇Φ)_{tt}∥^{2}.

Here we have

(σ(U)Φ∇Φ)_{tt} = σ(U)_{tt}Φ∇Φ + σ(U)Φ_{tt}∇Φ
+ σ(U)Φ∇Φ_{tt} + 2σ(U)_{t}Φ_{t}∇Φ
+ 2σ(U)_{t}Φ∇Φ_{t} + 2σ(U)Φ_{t}∇Φ_{t}.

(4.22)

In order to estimate the terms on the right we shall repeatedly use bounds from (4.17) together with

∥U_{t}∥_{0,∞} ≤ C∥U_{t}∥_{2,2} ≤ C(1 + ∥U_{tt}∥),

which follows from Sobolev’s inequality and (4.19). Thus,

∥σ(U)_{tt}Φ∇Φ∥ = ∥(σ′(U)U_{tt} + σ′′(U)U_{t}^{2})Φ∇Φ∥
≤ C(∥U_{tt}∥ + ∥U_{t}∥ ∥U_{t}∥_{0,∞})∥Φ∥_{1,∞}^{2}
≤ C(1 + ∥U_{tt}∥),

and

∥σ(U)_{t}Φ_{t}∇Φ∥ + ∥σ(U)_{t}Φ∇Φ_{t}∥ ≤ C∥U_{t}∥_{0,∞}∥Φ_{t}∥_{1,2}∥Φ∥_{1,∞} ≤ C(1 + ∥U_{tt}∥).

Similarly

∥σ(U)Φ_{t}∇Φ_{t}∥ ≤ C∥Φ_{t}∥_{0,∞}∥Φ_{t}∥_{1,2} ≤ C∥Φ_{t}∥_{0,∞}.

Application of Meyers’ theorem to equation (4.12) shows that

\begin{aligned}∥∇Φ_{t}∥_{0,p} & ≤ C(∥σ′(U)U_{t}∇Φ∥_{0,p} + ∥σ(U)∇g_{t}∥_{0,p}) \\ & ≤ C(1 + ∥U_{t}∥_{0,∞}∥Φ∥_{1,p}) ≤ C(1 + ∥U_{tt}∥), \end{aligned}

where p ∈ (2, r] is the same as before, and since ∥Φ_{t}∥_{0,∞} ≤ C∥Φ_{t}∥_{1,p} we may conclude that

∥σ(U)Φ_{t}∇Φ_{t}∥ ≤ C(1 + ∥U_{tt}∥).

Finally, for the remaining two terms in (4.22) we have

∥σ(U)Φ_{tt}∇Φ∥ + ∥σ(U)Φ∇Φ_{tt}∥ ≤ C∥Φ∥_{1,∞}∥Φ_{tt}∥_{1,2} ≤ C∥Φ_{tt}∥_{1,2}.

In order to bound ∥Φ_{tt}∥_{1,2} we differentiate equation (4.12) with respect to t to get

(σ(U)∇Φ_{tt}, ∇χ) = −(σ(U)_{tt}∇Φ + 2σ(U)_{t}∇Φ_{t}, ∇χ), ∀χ ∈ H_{0}^{1}.

With χ = Φ_{tt} − g_{tt} this gives

\begin{aligned}∥Φ_{tt}∥_{1,2} & ≤ C(∥g_{tt}∥_{1,2} + (∥U_{tt}∥ + ∥U_{t}∥ ∥U_{t}∥_{0,∞})∥Φ∥_{1,∞} + ∥U_{t}∥_{0,∞}∥Φ_{t}∥_{1,2}) \\ & ≤ C(1 + ∥U_{tt}∥). \end{aligned}

Together these estimates show that

∥(σ(U)Φ∇Φ)_{tt}∥ ≤ C(1 + ∥U_{tt}∥).

If we substitute this into (4.21), multiply by t^{2} and integrate, then we get

t^{2}∥U_{tt}(t)∥^{2} + ∫_{0}^{t} s^{2}∥∇U_{tt}∥^{2} ds ≤ C(1 + ∫_{0}^{t} s∥U_{tt}∥^{2} ds) ≤ C,

in view of (4.20), and the proof of the a priori bounds is complete.

Step 4. Finally, in order to prove uniqueness we let u_{1}, ϕ_{1} and u_{2}, ϕ_{2} be two solutions of (1.1)–(1.2). Using the a priori bounds ∥ϕ_{i}∥_{1,∞} ≤ C, i = 1, 2, it is straightforward to show that ∥ϕ_{1} − ϕ_{2}∥_{1,2} ≤ C∥u_{1} − u_{2}∥, and \frac{d}{dt}∥u_{1} − u_{2}∥^{2} ≤ C∥u_{1} − u_{2}∥^{2}, so that uniqueness follows. □

References

1. R. A. Adams and J. Fournier, Cone conditions and properties of Sobolev spaces, J. Math. Anal. Appl. 61 (1977), 713–734.

2. W. Allegretto and H. Xie, Existence of solutions for the time-dependent thermistor equations, IMA J. Appl. Math. 48 (1992), 271–281.

3. A. Bensoussan, J. L. Lions and G. Papanicolaou, Asymptotic Analysis for Periodic Structures, North-Holland, Amsterdam, 1978.

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

5. G. Cimatti, Existence of weak solutions for the nonstationary problem of the Joule heating of a conductor, Ann. Mat. Pura Appl. (4) 162 (1992), 33–42.

6. J. Douglas, Jr., R. E. Ewing and M. F. Wheeler, The approximation of the pressure by a mixed method in the simulation of miscible displacement, RAIRO Anal. Numér. 17 (1983), 17–33.

7. J. Douglas, Jr., R. E. Ewing and M. F. Wheeler, A time-discretization procedure for a mixed finite element approximation of miscible displacement in porous media, RAIRO Anal. Numér. 17 (1983), 249–265.

8. P. Grisvard, Le problème de Dirichlet dans l’espace W_{p}^{1}, Portugal. Math. 43 (1986), 393–398.

9. P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, New York, 1985.

10. D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Math., vol. 840, Springer-Verlag, Berlin, Heidelberg, New York, 1981.

11. C. Johnson, S. Larsson, V. Thomée and L. B. Wahlbin, Error estimates for spatially discrete approximations of semilinear parabolic equations with nonsmooth initial data, Math. Comp. 49 (1987), 331–357.

12. J. L. Lions and E. Magenes, Problèmes aux limites non homogènes III, Ann. Scuola Norm. Sup. Pisa 15 (1961), 39–101.

13. N. G. Meyers, An L^{p} estimate for the gradient of solutions of second order elliptic divergence equations, Ann. Scuola Norm. Sup. Pisa 17 (1963), 189–206.

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

15. J.-F. Rodrigues, A nonlinear parabolic system arising in thermomechanics and in thermomagnetism, Math. Models Methods Appl. Sci. 2 (1992), 271–281.

16. P. Saavedra and L. R. Scott, Variational formulation of a model free-boundary problem, Math. Comp. 57 (1991), 451–475.

17. V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Lecture Notes in Math., vol. 1054, Springer-Verlag, Berlin and New York, 1984.

18. X.-Y. Yue, Numerical analysis of nonstationary thermistor problem, J. Comput. Math. 12 (1994), 213–223.

Centre for Mathematical Analysis and Its Applications, School of Mathematical and Physical Sciences, University of Sussex, Brighton BN1 9QH, England

E-mail address: C.M.Elliott@sussex.ac.uk

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

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