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  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.1 Classical Statement of the Problem
Let be a bounded domain with Lipschitz boundary . We assume that Ω contains an interior subdomain with Lipschitz boundary Γ. In general, 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],
where , the coefficients , , w is measurable, and . Typically, in biophysical applications, is occupied by one or more macromolecules, and 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:
in Ω and ,
in , in and ,
in , in and .
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 on can also be treated in this framework provided that the boundary condition is defined as the trace of a function g such that and with .
The ability to find reliable and efficient solutions of the nonlinear Poisson–Boltzmann equation (PBE) for complex geometries of the interior domain (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  and the references therein. Although the solution of the linearized PBE (as in the linear Debye–Hückel theory) often yields accurate approximations , 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 ). 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 to its dual . For this (more general) weak formulation, we can guarantee existence of a solution and prove its uniqueness using a result of Brezis and Browder . Additionally, in Proposition 2.1, we show that the solution is bounded (here  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 and a nonlinear measure. In the case of a linear elliptic equation of the form , this nonlinear measure reduces to , where v and are approximations to the exact solution u and the exact flux . One advantage of the presented error estimate is that it is valid for any conforming approximations of u and 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 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 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
A function is called a weak solution of (1.1) if u is such that for any and
where and .
Define the functional by the relation
and consider the variational problem:
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 . Notice that, for , the function for all (e.g., see [27, 15]) and, therefore, the set is a linear subspace of . However, if , then is a convex set but not a linear subspace (Remark 2.1). Since is convex and obviously J is convex over , it follows that J is convex over (e.g., see ). Next we note that J is a proper functional, i.e., J is not identically equal to (e.g., ) and does not take the value ( is nonnegative). Therefore, existence of the minimizer u is guaranteed by known theorems of the calculus of variations (e.g., see ) if we will show that
J is sequentially weakly lower semicontinuous (s.w.l.s.c.), i.e., for any sequence weakly converging to v in ) ( ),
J is coercive, i.e., whenever .
To prove that (1) is fulfilled, notice that J is the sum of the functionals
The first functional is convex and continuous in and, therefore, it is s.w.l.s.c. (sequentially weakly lower semicontinuous). The second functional is convex and, for , 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 , the functional 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 into .
Let be a sequence weakly converging in to a , i.e., . Since the embedding is compact, it follows that (strongly) in . Therefore, we can extract a subsequence , which converges almost everywhere in the pointwise sense. Recall that for all and that is a continuous function for a.e. . Hence, by Fatou’s lemma, we obtain
Now it is clear that if is an arbitrary subsequence of , then there exists a further subsequence for which (2.4) is satisfied. This means that (2.4) is also satisfied for the whole sequence , and hence is s.l.w.s.c.
The coercivity of J follows from the estimate
where is the constant in the Friedrichs inequality for all .
Thus conditions (1) and (2) are fulfilled, and the existence of a minimizer is guaranteed. Moreover, noting that J is strictly convex, we arrive at the following result.
There exists a unique minimizer of problem (2.3).
On the other hand, if , let , and let the inner domain be given by for some , where denotes the ball in with radius r centered at zero. Consider the function . Since and for any , we find on the one hand that
On the other hand,
Hence , but for any and, therefore, is not a linear subspace. However, is a convex set. Indeed, let , i.e., . Since is convex for almost all and any , we have
Hence is a convex set.
The functional is not Gateaux differentiable at any if (therefore, J is also not Gateaux differentiable). In fact, is discontinuous at every . This fact is easy to see by the following example. Let be the ball centered at 0 with radius 2. There exists a function such that for any . In particular, we can set , where ϕ is a smooth function equal to 1 in and 0 in . Then , but for any since for small enough . In this case, for any and any , we have
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 . Since for all and any , we have
Making equivalent transformations of (2.5) and dividing by , we obtain
To compute the limit in the second term of (2.6), we will apply the Lebesgue dominated convergence theorem. We have
By the mean value theorem, we obtain
where and , a.e. . Then
from which it follows that
2.2 Uniqueness of the Solution to (2.1)
Uniqueness of the solution of (2.1) follows from the monotonicity of :
If are two different solutions of (2.1), then
Let Ω be an open set in , , and . If there exists a function such that a.e in Ω, then and the duality product in coincides with .
We have the following situation: a locally summable function defines a bounded linear functional over the dense subspace of through the integral formula . It is clear that the functional is uniquely extendable by continuity to a bounded linear functional over the whole space . The question is whether this extension is still representable by the same integral formula for any (if the integral makes sense at all). If the function is fixed, then Theorem 2.2 gives a sufficient condition for gv to be summable and for the extension evaluated at v to be representable with the same integral formula as above, i.e., .
Now, applying Theorem 2.2 to the functional defined by
and the function , using (2.11), we conclude that
Using the monotonicity (2.10) of and the coercivity of , we obtain .
2.3 Boundedness of the Minimizer
Let be a nonnegative function, which is nonincreasing for and such that
where C and α are positive constants and . If is defined by , then .
Now we present a main result of this section.
The unique weak solution u to problem (2.1) belongs to . Moreover, there exists a positive constant , depending only on d, Ω, , , such that . If , then the constant is equal to zero.
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
where (we notice that similar test functions have been used in [16, Theorem B.2] in the context of linear elliptic problems).
Choosing , using the monotonicity of , and the fact that , we obtain
which shows that assumption (2.13) holds for .
Now we are ready to prove that . First we consider the case . From (2.14), it follows that
Moreover, using the definition of and the definition (2.12) of , we obtain
where is the constant in Friedrichs’ inequality that holds for all . Finally, using (2.1), (2.15), and (2.16), we get for all . Consequently, almost everywhere and for all . In the case where l is not identically zero in Ω, we further estimate from below and from above using the Sobolev embedding , where for , for , and for . Let denote the Hölder conjugate to q. Then for , for , and for . In order to treat both cases in which we are interested simultaneously, namely, , we can take and . By we denote the embedding constant in the inequality for all , which depends only on the domain Ω and d. For , we have
and, for , we obtain
The final step before applying Lemma 2.1 is to estimate the left-hand side of (2.19) from below in terms of for and the right-hand side of (2.19) from above in terms of . Again, using the Sobolev embedding and Hölder’s inequality yields
Since , applying Lemma 2.1, we conclude that there is some such that
and . Hence . ∎
The results of this section are summarized in the following theorem.
Since in , and , we conclude that (2.1) holds for all resulting in a standard weak formulation. If is uniformly positive in the whole domain Ω and , then . On the other hand, if in , is uniformly positive in , and , we have .
3 A Posteriori Error Estimates
3.1 Abstract Framework
Here are reflexive Banach spaces with the norms and , respectively, , are convex and proper functionals, and is a bounded linear operator. By 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 and , respectively. They are endowed with the norms and . Henceforth, denotes the duality product of and . Analogously, is the duality product of and , and is the operator adjoint to Λ. It is defined by the relation
The functional defined by the relation
where and are the functionals conjugate to G and F, respectively. Problems (P) and (P*) are generated by the Lagrangian defined by the relation . If we additionally assume that is coercive and that is finite, then it is well known that problems (P) and (P*) have unique solutions and and that strong duality relations hold (see , or the book [8, Proposition 2.3, Remark 2.3, and Proposition 1.2 in Chapter VI]):
Furthermore, the pair is a saddle point of the Lagrangian L, i.e.,
and u and satisfy the relations
are the compound functionals for G and F, respectively . A compound functional is nonnegative by the definition. Equality (3.1) shows that and can vanish simultaneously if and only if and . Moreover, setting and in (3.1), we obtain analogous identities for the primal and dual parts of the error:
Notice that depends on the approximations v and only and, therefore, is fully computable. The right-hand side of (3.3) can be viewed as a certain measure of the distance between and , which vanishes if and only if and . Hence the relation
establishes the equality of the computable term 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 , , , and (where A is a symmetric positive definite matrix with bounded entries), then
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 , , where , and Λ the gradient operator . We further denote , , and , . With this notation, we have
and the functional J, defined by (2.2), can be written in the form . For any the functional is finite, while may take the value for some if (e.g., , on the unit ball in ). However, if , then for all and (see ). Also, is obviously finite since . We set and . In this case, coincides with considered as an operator from to . 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 (or equivalently ) which in some cases may dominate the usual energy norm terms. We start by deriving explicit expressions for , 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 . For and an arbitrary function , we introduce the functional
Recalling that the nonlinearity B is supported on , we have
Here is computed by the condition
which is equivalent to
We notice that (3.8) is a necessary condition for a maximum which is also sufficient since is convex. The solution of the last equation exists, is unique, and is given by
where and for . Note that the exact solution of dual problem (P*) also satisfies the relation because, for any , it holds . Moreover, since , , and , we see that and thus . In Proposition 3.1, we will later prove that we have not overestimated the supremum over in (3.7) and that we actually have equalities everywhere. Denoting and using the expression for and the formula
for any with in , we obtain an explicit formula for :
Since for all , the function belongs to for any , and we conclude that if . Therefore, the integral in (3.9) is well defined.
Now our goal is to prove that the inequality holds as the equality. In other words, we want to prove that the error estimate remains sharp and that the computed majorant will be indeed zero if approximations coincide with the exact solution .
For any with in , it holds
The idea is to approximate and by functions (in the a.e. sense) and use the Lebesgue dominated convergence theorem. Let and be such that a.e. in , (see [2, Theorem 4.9]), a.e. in , , where . Then
and (by extending the functions by zero in ). Since is continuous, we have the pointwise a.e. in convergence
Now we search for a function in that majorates the function :
Our next goal is to bound in (3.10). For the first summand, we have
where Remark 3.1 has been used. However, this bound cannot be used in the second term because might not belong even to . In order to find an -majorant for the second summand in (3.10), we distinguish the following two cases: In the first case, . Then .
In the second case ( ), we have . Therefore, . Since is a monotonically increasing function, . From here, we obtain
and using the relation , we conclude that
Finally, for almost all , we have
where, in the last line, we used the fact that . All the conditions of the Lebesgue dominated convergence theorem are satisfied, and we see that and, consequently,
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 with in , the quantity is fully computable and is given by the relation
and we have used the expressions for and . It is clear that since it is the sum of the compound functionals generated by and evaluated at and , respectively. It therefore qualifies as an error indicator, provided that is chosen appropriately, which we demonstrate with numerical experiments in the next section. Using the expression for , we obtain
Now we find explicit expressions for the nonlinear measures and similar to the ones for the case of quadratic F in (3.5) for the linear elliptic equation . We will need the following assertion, which is easy to prove.
For all , it holds
Since, for the exact solution u, we have and for a.e. , we find that
Similarly, , where