Johannes Kraus, Svetoslav Nakov and Sergey I. Repin

Reliable Numerical Solution of a Class of Nonlinear Elliptic Problems Generated by the Poisson–Boltzmann Equation

Accessible
De Gruyter | Published online: June 13, 2019

Abstract

We consider a class of nonlinear elliptic problems associated with models in biophysics, which are described by the Poisson–Boltzmann equation (PBE). We prove mathematical correctness of the problem, study a suitable class of approximations, and deduce guaranteed and fully computable bounds of approximation errors. The latter goal is achieved by means of the approach suggested in [19] for convex variational problems. Moreover, we establish the error identity, which defines the error measure natural for the considered class of problems and show that it yields computable majorants and minorants of the global error as well as indicators of local errors that provide efficient adaptation of meshes. Theoretical results are confirmed by a collection of numerical tests that includes problems on 2D and 3D Lipschitz domains.

1 Introduction

1.1 Classical Statement of the Problem

Let Ω d , d = 2 , 3 be a bounded domain with Lipschitz boundary Ω . We assume that Ω contains an interior subdomain Ω 1 with Lipschitz boundary Γ. In general, Ω 1 may consist of several disconnected parts (in this case, all of them are assumed to have Lipschitz continuous boundaries). We consider a class of nonlinear elliptic equations motivated by the Poisson–Boltzmann equation (PBE), which is widely used for computation of electrostatic interactions in a system of biomolecules in ionic solution [23, 10, 11],

(1.1)

(1.1a) - ( ϵ u ) + k 2 sinh ( u + w ) = l in Ω 1 Ω 2 ,
(1.1b) [ u ] Γ = 0 ,
(1.1c) [ ϵ u n ] Γ = 0 ,
(1.1d) u = 0 on Ω ,

where Ω 2 :- Ω ( Ω 1 Γ ) , the coefficients ϵ , k L ( Ω ) , ϵ max ϵ ϵ min > 0 , w is measurable, and l L 2 ( Ω ) . Typically, in biophysical applications, Ω 1 is occupied by one or more macromolecules, and Ω 2 is occupied by a solution of water and moving ions. The coefficients ϵ and k represent the dielectric constant and the modified Debye–Hückel parameter, and u is the dimensionless electrostatic potential. Concerning the given functions k and w, we can identify three main cases:

  1. (a)

    k max k ( x ) k min > 0 in Ω and w L ( Ω ) ,

  2. (b)

    k ( x ) 0 in Ω 1 , k max k ( x ) k min > 0 in Ω 2 and w L ( Ω 2 ) ,

  3. (c)

    k ( x ) 0 in Ω 2 , k max k ( x ) k min > 0 in Ω 1 and w L ( Ω 1 ) .

Throughout the paper, the major attention is paid to case b, which arises when solving the PBE and which is the most interesting from the practical point of view. Cases a and c can be studied analogously (with some rather obvious modifications). The case with nonhomogeneous Dirichlet boundary condition u = g on Ω can also be treated in this framework provided that the boundary condition is defined as the trace of a function g such that g H 1 ( Ω ) L ( Ω ) and g L s ( Ω ) with s > max { 2 , d } .

The ability to find reliable and efficient solutions of the nonlinear Poisson–Boltzmann equation (PBE) for complex geometries of the interior domain Ω 1 (with Lipschitz boundary) and piecewise constant dielectrics is important for applications in biophysics and biochemistry, e.g., in modeling the effects of water and ion screening on the potentials in and around soluble proteins, nucleic acids, membranes, and polar molecules and ions; see [23] and the references therein. Although the solution of the linearized PBE (as in the linear Debye–Hückel theory) often yields accurate approximations [22], certain mathematical models are valid only if they are based on the nonlinear PBE.

Over the recent years, adaptive finite element methods have proved to be an adequate technique in the numerical solution of elliptic problems with sharp local features such as point sources, heterogeneous coefficients or nonsmooth boundaries or interfaces (e.g., see [4, 9] and also have been successfully used to solve the nonlinear PBE [5, 14]. Adaptivity heavily relies on reliable and efficient error indicators, which are typically developed in the framework of a posteriori error control methods. While the theory of a posteriori error estimates for linear elliptic partial differential equations is already well established and understood, it is far less developed for nonlinear problems. A posteriori error analysis based on functional estimates has already been successfully applied to variational nonlinear problems including obstacle problems in [20, 21]. The accuracy verification approach taken in this work is also based on arguments that are commonly used in duality theory and convex analysis and can be found, e.g., in [8, 17]. Fast solution methods for systems of nonlinear differential equations is another important issue that affects efficiency of computer simulation methods. Multigrid methods may provide optimal or nearly optimal algorithms (in terms of complexity) to perform this task (e.g., see [18]). However, a systematic discussion of this topic is beyond the framework of this paper.

The main questions studied in the paper are related to the well-posedness of problem (1.1) and a posteriori error estimation of its numerical solution. We use a suitable weak formulation (Definition 2.1), where the nonlinearity does not satisfy any polynomial growth condition, and consequently it does not induce a bounded mapping from H 0 1 ( Ω ) to its dual H - 1 ( Ω ) . For this (more general) weak formulation, we can guarantee existence of a solution and prove its uniqueness using a result of Brezis and Browder [3]. Additionally, in Proposition 2.1, we show that the solution is bounded (here [3] is used, again together with special test functions suggested in [26, 16]). Boundedness of the solution is important and later used in the derivation of functional a posteriori error estimates. Applying the general approach from [19, 17], we derive guaranteed and computable bounds of the difference between the exact solution and any function from the respective energy class in terms of the energy and combined energy norms (equation (3.20)). Moreover, we obtain an error identity (3.19) with respect to a certain measure for the error which is the sum of the usual combined energy norm | | | ( v - u ) | | | 2 + | | | y * - p * | | | * 2 and a nonlinear measure. In the case of a linear elliptic equation of the form - div ( ϵ u ) + u = l , this nonlinear measure reduces to v - u L 2 ( Ω ) 2 + div ( y * - p * ) L 2 ( Ω ) 2 , where v and y * are approximations to the exact solution u and the exact flux p * = ϵ u . One advantage of the presented error estimate is that it is valid for any conforming approximations of u and ϵ u and that it does not rely on Galerkin orthogonality or properties specific to the used numerical method. Another advantage is that only the mathematical structure of the problem is exploited, and therefore no mesh dependent constants are present in the estimate. Majorants of the error not only give guaranteed bounds of global (energy) error norms but also generate efficient error indicators (cf. (1.1a), Figures 12 and 13). Also, we derive a simple, but efficient lower bound for the error in the combined energy norm. Using only the error majorant, we obtain an analog of Cea’s lemma which forms a basis for the a priori convergence analysis of finite element approximations for this class of semilinear problems. Finally, we present three numerical examples that verify the accuracy of error majorants and minorants and confirm efficiency of the error indicator in mesh adaptive procedures.

The outline of the paper is as follows. In Section 2, we discuss correctness of problem (1.1) and prove an a priori L ( Ω ) estimate for the solution u. In Section 3, first we recall some facts from the duality theory and general a posteriori error estimation method for convex variational problems. Then we apply this abstract framework and derive explicit forms of all the respective terms. A special attention is paid to the general error identity that defines a combined error measure natural for the considered class of problems. At the end of Section 3, we prove convergence of the conforming finite element method based on P 1 Lagrange elements. In Section 4, we consider numerical examples associated with 2D and 3D problems and make a systematic comparison of numerical solutions computed by adaptive mesh refinements based on different indicators. The last section includes a summary of the presented results.

2 Variational Form of the Problem

From now on, we assume that the functions k and w fall in case b of Section 1.1.

Definition 2.1.

A function u H 0 1 ( Ω ) is called a weak solution of (1.1) if u is such that b ( x , u + w ) v L 1 ( Ω ) for any v H 0 1 ( Ω ) L ( Ω ) and

(2.1) a ( u , v ) + Ω b ( x , u + w ) v d x = Ω l v d x for all v H 0 1 ( Ω ) L ( Ω ) ,

where a ( u , v ) = Ω ϵ u v d x and b ( x , z ) :- k 2 ( x ) sinh ( z ) .

Define the functional J : H 0 1 ( Ω ) { + } by the relation

(2.2) J ( v ) :- { Ω [ ϵ ( x ) 2 | v | 2 + k 2 cosh ( v + w ) - l v ] d x if k 2 cosh ( v + w ) L 1 ( Ω ) , + if k 2 cosh ( v + w ) L 1 ( Ω ) ,

and consider the variational problem:

(2.3) find u H 0 1 ( Ω ) such that J ( u ) = min v H 0 1 ( Ω ) J ( v ) .

2.1 Existence of a Minimizer

We begin with proving that the variational problem is well-posed. First it is necessary to make some comments on specific features of the above defined variational problem associated with the term k 2 cosh ( v + w ) . Notice that, for d 2 , the function e v L 2 ( Ω ) for all v H 0 1 ( Ω ) (e.g., see [27, 15]) and, therefore, the set dom ( J ) :- { v H 0 1 ( Ω ) : J ( v ) < } is a linear subspace of H 0 1 ( Ω ) . However, if d = 3 , then dom ( J ) is a convex set but not a linear subspace (Remark 2.1). Since dom ( J ) is convex and obviously J is convex over dom J , it follows that J is convex over H 0 1 ( Ω ) (e.g., see [8]). Next we note that J is a proper functional, i.e., J is not identically equal to + (e.g., J ( 0 ) = Ω k 2 cosh ( w ) d x < ) and does not take the value - ( J ( v ) is nonnegative). Therefore, existence of the minimizer u is guaranteed by known theorems of the calculus of variations (e.g., see [8]) if we will show that

  1. (1)

    J is sequentially weakly lower semicontinuous (s.w.l.s.c.), i.e., J ( v ) lim inf n J ( v n ) for any sequence { v n } n = 1 H 0 1 ( Ω ) weakly converging to v in H 0 1 ( Ω ) ) ( v n v ),

  2. (2)

    J is coercive, i.e., lim n J ( v n ) = + whenever v n H 1 ( Ω ) .

To prove that (1) is fulfilled, notice that J is the sum of the functionals

Ω ( ϵ 2 | v | 2 - l v ) d x and Ω k 2 ( x ) cosh ( v + w ) d x .

The first functional is convex and continuous in H 0 1 ( Ω ) and, therefore, it is s.w.l.s.c. (sequentially weakly lower semicontinuous). The second functional is convex and, for d = 2 , it is Gateaux differentiable, which implies that it is also s.w.l.s.c. (the proof of this implication can be found in [24, Corollary 2.4]). However, for d = 3 , the functional Ω k 2 ( x ) cosh ( v + w ) d x is not Gateaux differentiable (Remark 2.2). Nevertheless, one can show that it is also s.w.l.s.c. For this purpose, we use Fatou’s lemma and compact embedding of H 0 1 ( Ω ) into L 2 ( Ω ) .

Let { v n } n = 1 H 0 1 ( Ω ) be a sequence weakly converging in H 0 1 ( Ω ) to a v H 0 1 ( Ω ) , i.e., v n v . Since the embedding H 0 1 ( Ω ) L 2 ( Ω ) is compact, it follows that v n v (strongly) in L 2 ( Ω ) . Therefore, we can extract a subsequence v n m ( x ) v ( x ) , which converges almost everywhere in the pointwise sense. Recall that k 2 ( x ) cosh ( z ( x ) + w ( x ) ) 0 for all z H 0 1 ( Ω ) and that k 2 ( x ) cosh ( ) is a continuous function for a.e. x Ω . Hence, by Fatou’s lemma, we obtain

(2.4) lim inf m Ω k 2 ( x ) cosh ( v n m ( x ) + w ( x ) ) d x Ω lim inf m k 2 ( x ) cosh ( v n m ( x ) + w ( x ) ) d x = Ω k 2 ( x ) cosh ( v ( x ) + w ( x ) ) d x .

Now it is clear that if { v n m } m = 1 is an arbitrary subsequence of { v n } n = 1 , then there exists a further subsequence { v n m s } s = 1 for which (2.4) is satisfied. This means that (2.4) is also satisfied for the whole sequence { v n } n = 1 , and hence Ω k 2 cosh ( v + w ) d x is s.l.w.s.c.

The coercivity of J follows from the estimate

J ( v ) 1 2 a ( v , v ) - l L 2 ( Ω ) v L 2 ( Ω ) ϵ min v L 2 ( Ω ) 2 - l L 2 ( Ω ) v H 1 ( Ω ) ϵ min 1 + C F 2 v H 1 ( Ω ) 2 - l L 2 ( Ω ) v H 1 ( Ω ) ,

where C F is the constant in the Friedrichs inequality v L 2 ( Ω ) C F v L 2 ( Ω ) for all v H 0 1 ( Ω ) .

Thus conditions (1) and (2) are fulfilled, and the existence of a minimizer u H 0 1 ( Ω ) is guaranteed. Moreover, noting that J is strictly convex, we arrive at the following result.

Theorem 2.1.

There exists a unique minimizer u H 0 1 ( Ω ) of problem (2.3).

Remark 2.1.

It is worth noting that dom ( J ) :- { v H 0 1 ( Ω ) : k 2 ( x ) cosh ( v + w ) L 1 ( Ω ) } is a linear subspace of H 0 1 ( Ω ) for d 2 and not a linear subspace of H 0 1 ( Ω ) if d 3 . In dimension d 2 , from [27, 15], we know that e v L 2 ( Ω ) for any v H 0 1 ( Ω ) , and thus e λ v 1 + μ v 2 L 2 ( Ω ) for any λ , μ and any v 1 , v 2 H 0 1 ( Ω ) .

On the other hand, if d 3 , let Ω = B ( 0 , 1 ) , and let the inner domain Ω 2 be given by Ω 2 = B ( 0 , r ¯ ) for some r ¯ < 1 , where B ( 0 , r ) denotes the ball in d with radius r centered at zero. Consider the function v = ln 1 | x | H 0 1 ( B ( 0 , 1 ) ) . Since e v = 1 | x | L 1 ( Ω 2 ) and e λ v = 1 | x | λ L 1 ( Ω 2 ) for any λ d ,[1] we find on the one hand that

Ω k 2 cosh ( v + w ) d x = Ω 2 k 2 ( e v + w + e - v - w ) 2 d x 1 2 k max 2 e w L ( Ω 2 ) Ω 2 ( e v + e - v ) d x 1 2 k max 2 e w L ( Ω 2 ) ( Ω 2 e v d x + | Ω 2 | ) < .

On the other hand,

Ω k 2 cosh ( λ v + w ) d x 1 2 Ω 2 k 2 e λ v + w d x 1 2 k min 2 e - w L ( Ω 2 ) Ω 2 e λ v d x = for any λ > d .

Hence v dom ( J ) , but λ v dom ( J ) for any λ d and, therefore, dom ( J ) is not a linear subspace. However, dom ( J ) H 0 1 ( Ω ) is a convex set. Indeed, let v 1 , v 2 dom ( J ) , i.e., k 2 cosh ( v 1 + w ) , k 2 cosh ( v 2 + w ) L 1 ( Ω ) . Since k 2 cosh ( ) is convex for almost all x Ω and any λ [ 0 , 1 ] , we have

Ω k 2 cosh ( λ v 1 + ( 1 - λ ) v 2 + w ) d x λ Ω k 2 cosh ( v 1 + w ) d x + ( 1 - λ ) Ω k 2 cosh ( v 2 + w ) d x < + .

Hence dom ( J ) is a convex set.

Remark 2.2.

The functional Ω k 2 cosh ( v + w ) d x is not Gateaux differentiable at any u H 0 1 ( Ω ) L ( Ω ) if d = 3 (therefore, J is also not Gateaux differentiable). In fact, Ω k 2 cosh ( v + w ) d x is discontinuous at every u H 0 1 ( Ω ) L ( Ω ) . This fact is easy to see by the following example. Let Ω 2 = B ( 0 , 2 ) Ω be the ball centered at 0 with radius 2. There exists a function v H 0 1 ( Ω ) such that Ω 2 e λ v d x = + for any λ > 0 . In particular, we can set v = ϕ | x | - 1 / 3 , where ϕ is a smooth function equal to 1 in B ( 0 , 1 ) and 0 in 3 B ( 0 , 2 ) . Then v H 0 1 ( Ω ) , but e λ v L 1 ( Ω 2 ) for any λ > 0 since e λ v > | x | - 3 for small enough | x | . In this case, for any u H 0 1 ( Ω ) L ( Ω ) and any λ > 0 , we have

Ω k 2 cosh ( u + λ v + w ) d x 1 2 Ω 2 k 2 e u + λ v + w d x k min 2 e - u + w L ( Ω 2 ) 2 Ω 2 e λ v d x = + .

Now our goal is to show that the minimizer u is a solution of (2.1). To prove this, we use the Lebesgue dominated convergence theorem and the fact that, at the unique minimizer u of J, we have k 2 cosh ( u + w ) L 1 ( Ω ) . Since J ( u + λ v ) - J ( u ) 0 for all v H 0 1 ( Ω ) L ( Ω ) and any λ 0 , we have

(2.5) 1 2 a ( u + λ v , u + λ v ) + Ω k 2 cosh ( u + λ v + w ) d x - Ω l ( u + λ v ) d x - 1 2 a ( u , u ) - Ω k 2 cosh ( u + w ) d x + Ω l u d x 0 .

Making equivalent transformations of (2.5) and dividing by λ > 0 , we obtain

(2.6) a ( u , v ) + lim λ 0 + Ω k 2 ( cosh ( u + λ v + w ) - cosh ( u + w ) ) λ d x - Ω l v d x 0 .

To compute the limit in the second term of (2.6), we will apply the Lebesgue dominated convergence theorem. We have

(2.7) f λ ( x ) :- k 2 ( x ) ( cosh ( u ( x ) + w ( x ) + λ v ( x ) ) - cosh ( u ( x ) + w ( x ) ) ) λ λ 0 + k 2 ( x ) sinh ( u ( x ) + w ( x ) ) v ( x ) for a.e. x Ω .

By the mean value theorem, we obtain

f λ ( x ) = k 2 ( x ) sinh ( ς ( x ) ) v ( x ) ,

where ς ( x ) :- u ( x ) + w ( x ) + Θ ( x ) λ v ( x ) and Θ ( x ) ( 0 , 1 ) , a.e. x Ω . Then

| f λ ( x ) | k 2 ( x ) ( e ς ( x ) + e - ς ( x ) 2 ) | v ( x ) | ,

from which it follows that

(2.8) | f λ | v L ( Ω ) k 2 e u + w + e - u - w 2 e | Θ ( x ) | λ | v ( x ) | v L ( Ω ) e v L ( Ω ) k 2 cosh ( u + w ) L 1 ( Ω ) for all λ 1 .

From the Lebesgue dominated convergence theorem, (2.7) and (2.8), it follows that the limit in (2.6) is equal to Ω k 2 sinh ( u + w ) v d x , and therefore we obtain

(2.9) a ( u , v ) + Ω b ( x , u + w ) v d x - Ω l v d x 0 for all v H 0 1 ( Ω ) L ( Ω ) .

Since the test functions belong to a linear manifold, (2.9) is equivalent to the weak formulation (2.1).

2.2 Uniqueness of the Solution to (2.1)

Uniqueness of the solution of (2.1) follows from the monotonicity of b ( x , ) :

(2.10) Ω ( b ( x , v + w ) - b ( x , z + w ) ) ( v - z ) d x 0 for all v , z H 0 1 ( Ω ) L ( Ω ) .

If u 1 , u 2 H 0 1 ( Ω ) are two different solutions of (2.1), then

(2.11) a ( u 1 - u 2 , v ) + Ω ( b ( x , u 1 + w ) - b ( x , u 2 + w ) ) v d x = 0 for all v H 0 1 ( Ω ) L ( Ω ) .

Note that the difference u 1 - u 2 is not necessarily in H 0 1 ( Ω ) L ( Ω ) . To show that we can test with u 1 - u 2 in (2.11), we apply a property of Sobolev spaces proved in [3].

Theorem 2.2.

Let Ω be an open set in R d , T H - 1 ( Ω ) L loc 1 ( Ω ) , and v H 0 1 ( Ω ) . If there exists a function f L 1 ( Ω ) such that T ( x ) v ( x ) f ( x ) a.e in Ω, then T v L 1 ( Ω ) and the duality product T , v in H - 1 ( Ω ) × H 0 1 ( Ω ) coincides with Ω T v d x .

We have the following situation: a locally summable function g L loc 1 ( Ω ) defines a bounded linear functional T g over the dense subspace C 0 ( Ω ) of H 0 1 ( Ω ) through the integral formula T g , φ = Ω g φ d x . It is clear that the functional T g is uniquely extendable by continuity to a bounded linear functional T ¯ g over the whole space H 0 1 ( Ω ) . The question is whether this extension is still representable by the same integral formula for any v H 0 1 ( Ω ) (if the integral makes sense at all). If the function v H 0 1 ( Ω ) is fixed, then Theorem 2.2 gives a sufficient condition for gv to be summable and for the extension T ¯ g evaluated at v to be representable with the same integral formula as above, i.e., T ¯ g , v = Ω g v d x .

Now, applying Theorem 2.2 to the functional T g defined by

T g , v :- b ( x , u 1 + w ) - b ( x , u 2 + w ) , v for all v H 0 1 ( Ω ) L ( Ω )

and the function v = u 1 - u 2 H 0 1 ( Ω ) , using (2.11), we conclude that

a ( u 1 - u 2 , u 1 - u 2 ) + Ω ( b ( x , u 1 + w ) - b ( x , u 2 + w ) ) ( u 1 - u 2 ) d x = 0 .

Using the monotonicity (2.10) of b ( x , ) and the coercivity of a ( , ) , we obtain u 1 = u 2 .

2.3 Boundedness of the Minimizer

Next we show that the solution to problem (2.1) is essentially bounded. To prove this, we need the following lemma [16].

Lemma 2.1.

Let φ ( t ) be a nonnegative function, which is nonincreasing for s 0 t < and such that

φ ( h ) C φ ( s ) β ( h - s ) α for all h > s > s 0 ,

where C and α are positive constants and β > 1 . If j R is defined by j α :- C φ ( s 0 ) β - 1 2 α β β - 1 , then φ ( s 0 + j ) = 0 .

Now we present a main result of this section.

Proposition 2.1.

The unique weak solution u to problem (2.1) belongs to L ( Ω ) . Moreover, there exists a positive constant j ¯ > 0 , depending only on d, Ω, l L 2 ( Ω ) , ϵ min , such that u L ( Ω ) w L ( Ω 2 ) + j ¯ . If l = 0 , then the constant j ¯ is equal to zero.

Proof.

To prove that u is bounded, we apply Theorem 2.2 once again.

The first step is to show that (2.1) holds for the test function

(2.12) v = G s ( u ) :- sgn ( u ) max { | u | - s , 0 } ,

where s w L ( Ω 2 ) (we notice that similar test functions G s have been used in [16, Theorem B.2] in the context of linear elliptic problems).

It is easy to see that G s ( 0 ) = 0 , this function is Lipschitz continuous and, therefore, G s ( u ) H 0 1 ( Ω ) (e.g., see [16, 12]). Next the functional T b defined by

T b , v :- Ω b ( x , u + w ) v d x for all v H 0 1 ( Ω ) L ( Ω )

is bounded and linear and b ( x , u + w ) L loc 1 ( Ω ) . This follows from (2.1) and from the fact that the functionals a ( u , ) and ( l , ) belong to H - 1 ( Ω ) . In view of Theorem 2.2, to show that T b , G s ( u ) = Ω b ( x , u + w ) G s ( u ) d x , it suffices to verify the inequality

(2.13) b ( x , u + w ) G s ( u ) f a.e. for some f L 1 ( Ω ) .

Choosing s w L ( Ω 2 ) , using the monotonicity of b ( x , ) , and the fact that b ( x , 0 ) = 0 , we obtain

(2.14) b ( x , u + w ) G s ( u ) = { b ( x , u + w ) ( u - s ) 0 for u > s , 0 for u [ - s , s ] , b ( x , u + w ) ( u + s ) 0 for u < - s ,

which shows that assumption (2.13) holds for f = 0 .

Now we are ready to prove that u L ( Ω ) . First we consider the case l = 0 . From (2.14), it follows that

(2.15) Ω b ( x , u + w ) G s ( u ) d x 0 .

Moreover, using the definition of a ( , ) and the definition (2.12) of G s ( u ) , we obtain

(2.16) a ( u , G s ( u ) ) = Ω ϵ u G s ( u ) d x = Ω ϵ G s ( u ) G s ( u ) d x ϵ min G s ( u ) L 2 ( Ω ) 2 ϵ min C F 2 G s ( u ) L 2 ( Ω ) 2 ,

where C F is the constant in Friedrichs’ inequality v L 2 ( Ω ) C F v L 2 ( Ω ) that holds for all v H 0 1 ( Ω ) . Finally, using (2.1), (2.15), and (2.16), we get G s ( u ) L 2 ( Ω ) 2 0 for all s w L ( Ω 2 ) . Consequently, | u | s almost everywhere and for all s w L ( Ω 2 ) . In the case where l is not identically zero in Ω, we further estimate a ( u , G s ( u ) ) from below and Ω l G s ( u ) d x from above using the Sobolev embedding H 1 ( Ω ) L q ( Ω ) , where q = for d = 1 , q < for d = 2 , and q = 2 d d - 2 for d 3 . Let q * denote the Hölder conjugate to q. Then q * = 1 for d = 1 , q * = q q - 1 > 1 for d = 2 , and q * = 2 d d + 2 for d > 2 . In order to treat both cases in which we are interested simultaneously, namely, d = 2 , 3 , we can take q = 6 and q * = 6 / 5 . By C E we denote the embedding constant in the inequality v L 6 ( Ω ) C E v H 1 ( Ω ) for all v H 1 ( Ω ) , which depends only on the domain Ω and d. For a ( u , G s ( u ) ) , we have

(2.17) a ( u , G s ( u ) ) = Ω ϵ G s ( u ) G s ( u ) d x ϵ min 1 + C F 2 G s ( u ) H 1 ( Ω ) 2

and, for Ω l G s ( u ) d x , we obtain

(2.18) Ω l G s ( u ) d x = A ( s ) l G s ( u ) d x l L q * ( A ( s ) ) G s ( u ) L q ( Ω ) C E l L q * ( A ( s ) ) G s ( u ) H 1 ( Ω ) ,

where A ( s ) :- { x Ω : | u ( x ) | > s } . Combining (2.17), (2.18), (2.15), and (2.1), we arrive at the estimate

(2.19) ϵ min 1 + C F 2 G s ( u ) H 1 ( Ω ) C E l L q * ( A ( s ) ) .

The final step before applying Lemma 2.1 is to estimate the left-hand side of (2.19) from below in terms of | A ( h ) | for h > s w L ( Ω 2 ) and the right-hand side of (2.19) from above in terms of | A ( s ) | . Again, using the Sobolev embedding H 1 ( Ω ) L q ( Ω ) and Hölder’s inequality yields

(2.20) G s ( u ) H 1 ( Ω ) 1 C E ( Ω | G s ( u ) | q d x ) 1 q = 1 C E ( A ( s ) | | u | - s | q d x ) 1 q 1 C E ( A ( h ) ( h - s ) q d x ) 1 q = 1 C E ( h - s ) | A ( h ) | 1 q

and

(2.21) l L q * ( A ( s ) ) l L 2 ( Ω ) | A ( s ) | 2 - q * 2 q * .

Combining (2.20), (2.21), and (2.19), we obtain the following inequality for the nonnegative and non-increasing function φ ( t ) :- | A ( t ) | :

| A ( h ) | ( C E 2 ( 1 + C F 2 ) ϵ min l L 2 ( Ω ) ) q | A ( s ) | q - 2 2 ( h - s ) q for all h > s w L ( Ω 2 ) .

Since q - 2 2 = 2 > 1 , applying Lemma 2.1, we conclude that there is some j > 0 such that

0 < j q = ( C E 2 ( 1 + C F 2 ) ϵ min l L 2 ( Ω ) ) q | A ( w L ( Ω 2 ) ) | q - 4 2 2 q ( q - 2 ) q - 4 ( C E 2 ( 1 + C F 2 ) ϵ min l L 2 ( Ω ) ) q | Ω | q - 4 2 2 q ( q - 2 ) q - 4 -: j ¯ q

and | A ( w L ( Ω 2 ) + j ¯ ) | = 0 . Hence u L ( Ω ) w L ( Ω 2 ) + j ¯ . ∎

The results of this section are summarized in the following theorem.

Theorem 2.3.

Problem (2.1) has a unique solution u H 0 1 ( Ω ) L ( Ω ) , which is the unique minimizer of variational problem (2.3).

Remark 2.3.

Since k = 0 in Ω 1 , w L ( Ω 2 ) and u L ( Ω ) , we conclude that (2.1) holds for all v H 0 1 ( Ω ) resulting in a standard weak formulation. If k 2 is uniformly positive in the whole domain Ω and w L ( Ω ) , then u L ( Ω ) w L ( Ω ) + j ¯ . On the other hand, if k = 0 in Ω 2 , k 2 is uniformly positive in Ω 1 , and w L ( Ω 1 ) , we have u L ( Ω ) w L ( Ω 1 ) + j ¯ .

3 A Posteriori Error Estimates

3.1 Abstract Framework

First we briefly recall some results from the duality theory [17, 8]. Consider a class of variational problems having the following common form:

(P) find u V such that J ( u ) = inf v V J ( v ) , where J ( v ) = G ( Λ v ) + F ( v ) .

Here V , Y are reflexive Banach spaces with the norms V and Y , respectively, F : V ¯ , G : Y ¯ are convex and proper functionals, and Λ : V Y is a bounded linear operator. By 0 V we denote the zero element in V. It is assumed that J is coercive and lower semicontinuous. In this case, problem (P) has a solution u, which is unique if J is strictly convex.

The spaces topologically dual to V and Y are denoted by V * and Y * , respectively. They are endowed with the norms V * and Y * . Henceforth, v * , v denotes the duality product of v * V * and v V . Analogously, ( y * , y ) is the duality product of y * Y * and y Y , and Λ * : Y * V * is the operator adjoint to Λ. It is defined by the relation

Λ * y * , v = ( y * , Λ v ) for all v V , y * Y * .

The functional J * : V * ¯ defined by the relation

J * ( v * ) :- sup v V { v * , v - J ( v ) }

is called dual (or Fenchel conjugate) to J (see, e.g., [8]). In accordance with the general duality theory of the calculus of variations, the primal problem (P) has a dual counterpart:

(P*) find p * Y * such that I * ( p * ) = sup y * Y * I * ( y * ) , where I * ( y * ) :- - G * ( y * ) - F * ( - Λ * y * ) ,

where G * and F * are the functionals conjugate to G and F, respectively. Problems (P) and (P*) are generated by the Lagrangian L : V × Y * ¯ defined by the relation L ( v , y * ) = ( y * , Λ v ) - G * ( y * ) + F ( v ) . If we additionally assume that G * is coercive and that F ( 0 V ) is finite, then it is well known that problems (P) and (P*) have unique solutions u V and p * Y * and that strong duality relations hold (see [17], or the book [8, Proposition 2.3, Remark 2.3, and Proposition 1.2 in Chapter VI]):

J ( u ) = inf v V J ( v ) = inf v V sup y * Y * L ( v , y * ) = sup y * Y * inf v V L ( v , y * ) = sup y * Y * I * ( y * ) = I * ( p * ) .

Furthermore, the pair ( u , p * ) is a saddle point of the Lagrangian L, i.e.,

L ( u , y * ) L ( u , p * ) L ( v , p * ) for all v V , y * Y * ,

and u and p * satisfy the relations

Λ u G * ( p * ) , p * G ( Λ u ) .

We have

(3.1) J ( v ) - I * ( y * ) = G ( Λ v ) + F ( v ) + G * ( y * ) + F * ( - Λ * y * ) = D G ( Λ v , y * ) + D F ( v , - Λ * y * ) -: M 2 ( v , y * ) ,

where

D G ( Λ v , y * ) :- G ( Λ v ) + G * ( y * ) - ( y * , Λ v ) ,
D F ( v , - Λ * y * ) :- F ( v ) + F * ( - Λ * y * ) + Λ * y * , v

are the compound functionals for G and F, respectively [17]. A compound functional is nonnegative by the definition. Equality (3.1) shows that D G and D F can vanish simultaneously if and only if v = u and y * = p * . Moreover, setting v :- u and y * :- p * in (3.1), we obtain analogous identities for the primal and dual parts of the error:

(3.2)

(3.2a) J(u)-I^*(y^*) = M 2 ( u , y * ) = D G ( Λ u , y * ) + D F ( u , - Λ * y * ) ,
(3.2b) J ( v ) - I * ( p * ) = M 2 ( v , p * ) = D G ( Λ v , p * ) + D F ( v , - Λ * p * ) .

Using the fact that J ( u ) = I * ( p * ) and that the above equalities (3.2) hold, we obtain another important relation (see [17])

(3.3) M 2 ( v , y * ) = J ( v ) - I * ( y * ) = J ( v ) - I * ( p * ) + J ( u ) - I * ( y * ) = M 2 ( v , p * ) + M 2 ( u , y * ) .

Notice that M 2 ( v , y * ) depends on the approximations v and y * only and, therefore, is fully computable. The right-hand side of (3.3) can be viewed as a certain measure of the distance between ( u , p * ) and ( v , y * ) , which vanishes if and only if v = u and y * = p * . Hence the relation

(3.4) D G ( Λ v , p * ) + D F ( v , - Λ * p * ) + D G ( Λ u , y * ) + D F ( u , - Λ * y * ) = M 2 ( v , y * )

establishes the equality of the computable term M 2 ( v , y * ) and an error measure natural for this class of variational problems.

It is worth noting that identity (3.4) can be represented in terms of norms if G and F are quadratic functionals. For example, if V = H 0 1 ( Ω ) , V * = H - 1 ( Ω ) , Y = [ L 2 ( Ω ) ] d = Y * , G ( Λ v ) = G ( v ) = Ω 1 2 A v v d x and F ( v ) = Ω ( 1 2 v 2 - l v ) d x (where A is a symmetric positive definite matrix with bounded entries), then

(3.5)

D G ( Λ v , p * ) = 1 2 Ω A ( v - u ) ( v - u ) d x , D G ( Λ u , y * ) = 1 2 Ω A - 1 ( y * - p * ) ( y * - p * ) d x ,
D F ( v , - Λ * p * ) = 1 2 v - u L 2 ( Ω ) 2 , D F ( u , - Λ * y * ) = 1 2 div ( y * - p * ) L 2 ( Ω ) 2 .

In this case, the minimizer of (P) solves the linear elliptic problem - div ( A u ) + u = l in Ω, and (3.4) is reduced to the error identity

(3.6) Ω A ( v - u ) ( v - u ) d x + Ω A - 1 ( y * - p * ) ( y * - p * ) d x + v - u L 2 ( Ω ) 2 + div ( y * - p * ) L 2 ( Ω ) 2 = | | | A v - y * | | | * 2 + v - div y * - l L 2 ( Ω ) 2 = 2 M 2 ( v , y * ) .

The sum of the first and the third term in (3.6) represents the primal, the sum of the second and fourth term the dual error.

We set V :- H 0 1 ( Ω ) , Y :- [ L 2 ( Ω ) ] d , where d = 2 , 3 , and Λ the gradient operator : H 0 1 ( Ω ) [ L 2 ( Ω ) ] d . We further denote g : Ω × 3 , g ( x , ξ ) :- ϵ ( x ) 2 | ξ | 2 , and B : Ω × , B ( x , ξ ) :- k 2 ( x ) cosh ( ξ ) . With this notation, we have

G ( Λ v ) :- Ω g ( x , v ( x ) ) d x = Ω ϵ 2 | v | 2 d x ,
F ( v ) :- Ω B ( x , v ( x ) + w ( x ) ) d x - Ω l v d x = Ω k 2 cosh ( v + w ) d x - Ω l v d x ,

and the functional J, defined by (2.2), can be written in the form J ( v ) = G ( Λ v ) + F ( v ) . For any v V the functional G ( Λ v ) is finite, while F : V { + } may take the value + for some v V if d 3 (e.g., v = log 1 | x | α , α d on the unit ball in d ). However, if d 2 , then exp ( v ) L 1 ( Ω ) for all v H 0 1 ( Ω ) and F : V (see [15]). Also, F ( 0 V ) is obviously finite since w L ( Ω 2 ) . We set V * = H - 1 ( Ω ) and Y * = Y = [ L 2 ( Ω ) ] d . In this case, Λ * coincides with - div considered as an operator from [ L 2 ( Ω ) ] d to H - 1 ( Ω ) . We will present the particular form of error equality (3.4) where the error is measured in a special “nonlinear norm”. This measure contains the usual combined energy norm terms, i.e., the sum of the energy norms of the errors for the primal and dual problem, but also two additional nonnegative terms due to the nonlinearity B ( x , ξ ) (or equivalently b ( x , ξ ) ) which in some cases may dominate the usual energy norm terms. We start by deriving explicit expressions for G * , F * , and then we will use these expressions to get an explicit form of the abstract error equality (3.4).

3.2 Fenchel Conjugates of the Functionals G and F

It is easy to find that G * ( y * ) = Ω 1 2 ϵ ( x ) | y * ( x ) | 2 d x . For y * H ( div ; Ω ) and an arbitrary function z : Ω 2 , we introduce the functional

I y * ( z ) :- Ω 2 [ ( div y * + l ) z - B ( x , z + w ) ] d x .

Recalling that the nonlinearity B is supported on Ω 2 , we have

(3.7) F * ( - Λ * y * ) = sup z H 0 1 ( Ω ) [ - Λ * y * , z - F ( z ) ] = sup z H 0 1 ( Ω ) [ ( - y * , Λ z ) - F ( z ) ] = sup z H 0 1 ( Ω ) Ω [ - y * z - B ( x , z + w ) + l z ] d x ( if y * H ( div ; Ω ) ) = sup z H 0 1 ( Ω ) Ω [ div y * z - B ( x , z + w ) + l z ] d x ( finite if div y * + l = 0 in Ω 1 ) = sup z H 0 1 ( Ω ) I y * ( z ) Ω 2 sup ξ [ ( div y * ( x ) + l ( x ) ) ξ - B ( x , ξ + w ( x ) ) ] d x = Ω 2 [ ( div y * ( x ) + l ( x ) ) ξ 0 ( x ) - B ( x , ξ 0 ( x ) + w ( x ) ) ] d x = I y * ( ξ 0 ) .

Here ξ 0 : Ω 2 is computed by the condition

(3.8) d d ξ [ ( div y * ( x ) + l ( x ) ) ξ - B ( x , ξ + w ( x ) ) ] = 0 for a.e. x Ω 2 ,

which is equivalent to

div y * ( x ) + l ( x ) - k 2 ( x ) sinh ( ξ + w ( x ) ) = 0 for a.e. x Ω 2 .

We notice that (3.8) is a necessary condition for a maximum which is also sufficient since B ( x , ) is convex. The solution of the last equation exists, is unique, and is given by

ξ 0 ( x ) = arsinh ( ρ k ( y * ) ) - w ( x ) = ln ( ρ k ( y * ) + ρ k 2 ( y * ) + 1 ) - w ( x ) = ln ( Θ ( ρ k ( y * ) ) ) - w ( x ) ,

where ρ k ( y * ) :- div y * ( x ) + l ( x ) k 2 ( x ) and Θ ( s ) :- s + s 2 + 1 for s . Note that the exact solution p * = ϵ u of dual problem (P*) also satisfies the relation div ( ϵ u ) + l = 0 because, for any x Ω 1 , it holds k ( x ) = 0 . Moreover, since u L ( Ω ) , w L ( Ω 2 ) , and l L 2 ( Ω ) , we see that div p * = k 2 sinh ( u + w ) + l L 2 ( Ω ) and thus p * H ( div ; Ω ) . In Proposition 3.1, we will later prove that we have not overestimated the supremum over z H 0 1 ( Ω ) in (3.7) and that we actually have equalities everywhere. Denoting S :- arsinh ( ρ k ( y * ) ) and using the expression for ξ 0 ( x ) and the formula

cosh ( arsinh ( x ) ) = x 2 + 1 for all x ,

for any y * H ( div ; Ω ) [ L 2 ( Ω ) ] d = Y * with div y * + l = 0 in Ω 1 , we obtain an explicit formula for F * ( - Λ * y * ) :

(3.9) F * ( - Λ * y * ) = Ω 2 [ k 2 ρ k ( y * ) ( ln ( Θ ( ρ k ( y * ) ) ) - w ) - k 2 ρ k 2 ( y * ) + 1 ] d x = Ω 2 [ k 2 sinh ( S ) ( S - w ) - k 2 cosh ( S ) ] d x .

Remark 3.1.

Since | ln ( t + t 2 + 1 ) | | t | for all t , the function ln ( Θ ( f ( x ) ) ) - w ( x ) belongs to L 2 ( Ω 2 ) for any f L 2 ( Ω 2 ) , and we conclude that ξ 0 ( x ) L 2 ( Ω 2 ) if y * H ( div ; Ω ) . Therefore, the integral in (3.9) is well defined.

Now our goal is to prove that the inequality sup z H 0 1 ( Ω ) I y * ( z ) I y * ( ξ 0 ) holds as the equality. In other words, we want to prove that the error estimate remains sharp and that the computed majorant M 2 ( v , y * ) will be indeed zero if approximations ( v , y * ) coincide with the exact solution ( u , p * ) .

Proposition 3.1.

For any y * H ( div ; Ω ) with div y * + l = 0 in Ω 1 , it holds

sup z H 0 1 ( Ω ) I y * ( z ) = I y * ( ξ 0 ) < .

Proof.

The idea is to approximate f = div y * + l k 2 L 2 ( Ω 2 ) and w Ω 2 L ( Ω 2 ) by C 0 ( Ω 2 ) functions (in the a.e. sense) and use the Lebesgue dominated convergence theorem. Let f n C 0 ( Ω 2 ) and w n C 0 ( Ω 2 ) be such that f n ( x ) f ( x ) a.e. in Ω 2 , | f n ( x ) | h ( x ) L 2 ( Ω 2 ) (see [2, Theorem 4.9]), w n ( x ) w ( x ) a.e. in Ω 2 , | w n ( x ) | m + 2 , where m :- w L ( Ω 2 ) . Then

z n ( x ) :- ln ( Θ ( f n ( x ) ) ) - w n ( x ) ξ 0 ( x ) a.e. in Ω 2

and z n C 0 ( Ω 2 ) H 0 1 ( Ω 2 ) H 0 1 ( Ω ) (by extending the functions by zero in Ω 1 ). Since B ( x , ) is continuous, we have the pointwise a.e. in Ω 2 convergence

( div y * ( x ) + l ( x ) ) z n ( x ) - B ( x , z n + w ( x ) ) ( div y * ( x ) + l ( x ) ) ξ 0 ( x ) - B ( x , ξ 0 ( x ) + w ( x ) )

Now we search for a function in L 1 ( Ω 2 ) that majorates the function | ( div y * ( x ) + l ( x ) ) z n ( x ) - B ( x , z n + w ( x ) ) | :

(3.10) | ( div y * ( x ) + l ( x ) ) z n ( x ) - k 2 ( x ) cosh ( z n ( x ) + w ( x ) ) | | div y * ( x ) + l ( x ) | | z n ( x ) | + k 2 ( x ) e w L ( Ω 2 ) e | z n ( x ) | .

Our next goal is to bound | z n ( x ) | in (3.10). For the first summand, we have

| z n ( x ) | = | ln ( Θ ( f n ( x ) ) ) - w n ( x ) | | f n ( x ) | + m + 2 h ( x ) + m + 2 L 2 ( Ω 2 ) ,

where Remark 3.1 has been used. However, this bound cannot be used in the second term because e h might not belong even to L 1 ( Ω 2 ) . In order to find an L 1 -majorant for the second summand in (3.10), we distinguish the following two cases: In the first case, f n ( x ) > 0 . Then | ln ( Θ ( f n ( x ) ) ) | | ln ( Θ ( h ( x ) ) ) | .

In the second case ( f n ( x ) 0 ), we have Θ ( f n ( x ) ) 1 . Therefore, 0 f n ( x ) - h ( x ) . Since Θ ( s ) is a monotonically increasing function, Θ ( 0 ) = 1 Θ ( f n ( x ) ) Θ ( - h ( x ) ) > 0 . From here, we obtain

ln ( 1 ) = 0 ln ( Θ ( f n ( x ) ) ) ln ( Θ ( - h ( x ) ) ) ,

and using the relation Θ ( - h ) = 1 Θ ( h ) , we conclude that

| ln ( Θ ( f n ( x ) ) ) | | ln ( Θ ( - h ( x ) ) ) | = | ln ( Θ ( h ( x ) ) ) | .

Finally, for almost all x Ω 2 , we have

| z n ( x ) | = | ln ( Θ ( f n ( x ) ) ) - w n ( x ) | | ln ( Θ ( h ( x ) ) ) | + m + 2 = ln ( Θ ( h ( x ) ) ) + m + 2 because h ( x ) 0 for a.e. x Ω 2 .

Therefore,

| ( div y * ( x ) + l ( x ) ) z n ( x ) - k 2 ( x ) cosh ( z n ( x ) + w ( x ) ) | | div y * ( x ) + l ( x ) | ( h ( x ) + w L ( Ω 2 ) + 2 ) + k 2 ( x ) e 2 w L ( Ω 2 ) + 2 Θ ( h ( x ) ) :- H ( x ) L 2 ( Ω 2 ) ,

where, in the last line, we used the fact that Θ ( h ( x ) ) L 2 ( Ω 2 ) . All the conditions of the Lebesgue dominated convergence theorem are satisfied, and we see that I y * ( z n ) I y * ( ξ 0 ) and, consequently,

sup z H 0 1 ( Ω ) I y * ( z ) = I y * ( ξ 0 ) .

3.3 Error Measures

In this section, we apply the abstract framework from Section 3.1 and derive an explicit form of relation (3.4) adapted to our problem. For any y * H ( div ; Ω ) with div y * + l = 0 in Ω 1 , the quantity M 2 ( v , y * ) is fully computable and is given by the relation

(3.11) M 2 ( v , y * ) = D G ( Λ v , y * ) + D F ( v , - Λ * y * ) = G ( Λ v ) + G * ( y * ) - ( y * , Λ v ) + F ( v ) + F * ( - Λ * y * ) + Λ * y * , v = Ω η 2 ( x ) d x = 1 2 | | | ϵ v - y * | | | * 2 + D F ( v , - Λ * y * ) ,

where

(3.12) η 2 ( x ) = { 1 2 ϵ | ϵ v - y * | 2 , x Ω 1 1 2 ϵ | ϵ v - y * | 2 + k 2 cosh ( v + w ) - l v + k 2 ρ k ( y * ) ( ln ( Θ ( ρ k ( y * ) ) ) - w ) - k 2 ρ k 2 ( y * ) + 1 - div y * v , x Ω 2

and we have used the expressions for G * and F * . It is clear that η 2 ( x ) 0 since it is the sum of the compound functionals generated by g ~ x ( s ) :- g ( x , s ) and B ~ x ( s ) - l ( x ) s = B ( x , s + w ( x ) ) - l ( x ) s evaluated at ( v ( x ) , y * ( x ) ) and ( v ( x ) , div y * ( x ) ) , respectively. It therefore qualifies as an error indicator, provided that y * is chosen appropriately, which we demonstrate with numerical experiments in the next section. Using the expression for G * , we obtain

(3.13)  D_G( v,p^*) = 1 2 Ω ϵ | ( v - u ) | 2 d x -: 1 2 | | | ( v - u ) | | | 2 ,
D G ( Λ u , y * ) = 1 2 Ω 1 ϵ | y * - p * | 2 d x -: 1 2 | | | y * - p * | | | * 2 .

Now we find explicit expressions for the nonlinear measures D F ( v , - Λ * p * ) and D F ( u , - Λ * y * ) similar to the ones for the case of quadratic F in (3.5) for the linear elliptic equation - div ( A u ) + u = l . We will need the following assertion, which is easy to prove.

Proposition 3.2.

For all s , t R , it holds

(3.14) ( t - s ) 2 2 A ( s , t ) ( sinh ( t ) - sinh ( s ) ) 2 2 ,

where A ( s , t ) = cosh ( t ) - cosh ( s ) + s sinh ( s ) - t sinh ( s ) .

Since, for the exact solution u, we have ρ k ( p * ) = sinh ( u + w ) and u = arsinh ( ρ k ( p * ) ) - w for a.e. x Ω 2 , we find that

D F ( v , - Λ * p * ) = Ω 2 ( k 2 cosh ( v + w ) - l v + k 2 sinh ( u + w ) u - k 2 cosh ( u + w ) - div p * v ) d x = Ω 2 k 2 ( cosh ( v + w ) - cosh ( u + w ) + u sinh ( u + w ) - v sinh ( u + w ) ) d x .

Similarly, D F ( u , - Λ * y * ) = Ω 2 k 2 ( cosh ( T ) - cosh ( S ) + S sinh ( S ) - T sinh ( S ) ) d x , where T :- arsinh ( ρ k