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∈H01(Ω) is called a weak solution of (1.1) if u is such that b(x,u+w)v∈L1(Ω) for any v∈H01(Ω)∩L∞(Ω) and
(2.1)a(u,v)+∫Ωb(x,u+w)vdx=∫Ωlvdx for allv∈H01(Ω)∩L∞(Ω),
where a(u,v)=∫Ωϵ∇u⋅∇vdx and b(x,z):-k2(x)sinh(z).
Define the functional J:H01(Ω)→ℝ∪{+∞} by the relation
(2.2)J(v):-{∫Ω[ϵ(x)2|∇v|2+k2cosh(v+w)-lv]dxifk2cosh(v+w)∈L1(Ω),+∞ifk2cosh(v+w)∉L1(Ω),
and consider the variational problem:
(2.3)findu∈H01(Ω)such that J(u)=minv∈H01(Ω)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 k2cosh(v+w).
Notice that, for d≤2, the function ev∈L2(Ω) for all v∈H01(Ω) (e.g., see [27, 15]) and, therefore, the set dom(J):-{v∈H01(Ω):J(v)<∞} is a linear subspace of H01(Ω).
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 domJ, it follows that J is convex over H01(Ω) (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)=∫Ωk2cosh(w)dx<∞) 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
J is sequentially weakly lower semicontinuous (s.w.l.s.c.), i.e., J(v)≤lim infn→∞J(vn) for any sequence {vn}n=1∞⊂H01(Ω) weakly converging to v in H01(Ω)) (vn⇀v),
J is coercive, i.e., limn→∞J(vn)=+∞ whenever ∥vn∥H1(Ω)→∞.
To prove that (1) is fulfilled, notice that J is the sum of the functionals
∫Ω(ϵ2|∇v|2-lv)dx and ∫Ωk2(x)cosh(v+w)dx.
The first functional is convex and continuous in H01(Ω) 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 ∫Ωk2(x)cosh(v+w)dx 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 H01(Ω) into L2(Ω).
Let {vn}n=1∞⊂H01(Ω) be a sequence weakly converging in H01(Ω) to a v∈H01(Ω), i.e., vn⇀v.
Since the embedding H01(Ω)↪L2(Ω) is compact, it follows that vn→v (strongly) in L2(Ω).
Therefore, we can extract a subsequence vnm(x)→v(x), which converges almost everywhere in the pointwise sense.
Recall that k2(x)cosh(z(x)+w(x))≥0 for all z∈H01(Ω) and that k2(x)cosh(⋅) is a continuous function for a.e. x∈Ω.
Hence, by Fatou’s lemma, we obtain
(2.4)lim infm→∞∫Ωk2(x)cosh(vnm(x)+w(x))dx≥∫Ωlim infm→∞k2(x)cosh(vnm(x)+w(x))dx=∫Ωk2(x)cosh(v(x)+w(x))dx.
Now it is clear that if {vnm}m=1∞ is an arbitrary subsequence of {vn}n=1∞, then there exists a further subsequence {vnms}s=1∞ for which (2.4) is satisfied.
This means that (2.4) is also satisfied for the whole sequence {vn}n=1∞, and hence ∫Ωk2cosh(v+w)dx is s.l.w.s.c.
The coercivity of J follows from the estimate
J(v)≥12a(v,v)-∥l∥L2(Ω)∥v∥L2(Ω)≥ϵmin∥∇v∥L2(Ω)2-∥l∥L2(Ω)∥v∥H1(Ω)≥ϵmin1+CF2∥v∥H1(Ω)2-∥l∥L2(Ω)∥v∥H1(Ω),
where CF is the constant in the Friedrichs inequality ∥v∥L2(Ω)≤CF∥∇v∥L2(Ω) for all v∈H01(Ω).
Thus conditions (1) and (2) are fulfilled, and the existence of a minimizer u∈H01(Ω) is guaranteed.
Moreover, noting that J is strictly convex, we arrive at the following result.
Theorem 2.1.
There exists a unique minimizer u∈H01(Ω) of problem (2.3).
Remark 2.1.
It is worth noting that dom(J):-{v∈H01(Ω):k2(x)cosh(v+w)∈L1(Ω)} is a linear subspace of H01(Ω) for d≤2 and not a linear subspace of H01(Ω) if d≥3.
In dimension d≤2, from [27, 15], we know that ev∈L2(Ω) for any v∈H01(Ω), and thus eλv1+μv2∈L2(Ω) for any λ,μ∈ℝ and any v1,v2∈H01(Ω).
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=ln1|x|∈H01(B(0,1)).
Since ev=1|x|∈L1(Ω2) and eλv=1|x|λ∉L1(Ω2) for any λ≥d, we find on the one hand that
∫Ωk2cosh(v+w)dx=∫Ω2k2(ev+w+e-v-w)2dx≤12kmax2e∥w∥L∞(Ω2)∫Ω2(ev+e-v)dx≤12kmax2e∥w∥L∞(Ω2)(∫Ω2evdx+|Ω2|)<∞.
On the other hand,
∫Ωk2cosh(λv+w)dx≥12∫Ω2k2eλv+wdx≥12kmin2e-∥w∥L∞(Ω2)∫Ω2eλvdx=∞ 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)⊂H01(Ω) is a convex set.
Indeed, let v1,v2∈dom(J), i.e., k2cosh(v1+w),k2cosh(v2+w)∈L1(Ω).
Since k2cosh(⋅) is convex for almost all x∈Ω and any λ∈[0,1], we have
∫Ωk2cosh(λv1+(1-λ)v2+w)dx≤λ∫Ωk2cosh(v1+w)dx+(1-λ)∫Ωk2cosh(v2+w)dx<+∞.
Hence dom(J) is a convex set.
Remark 2.2.
The functional ∫Ωk2cosh(v+w)dx is not Gateaux differentiable at any u∈H01(Ω)∩L∞(Ω) if d=3 (therefore, J is also not Gateaux differentiable).
In fact, ∫Ωk2cosh(v+w)dx is discontinuous at every u∈H01(Ω)∩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∈H01(Ω) such that ∫Ω2eλvdx=+∞ 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∈H01(Ω), but eλv∉L1(Ω2) for any λ>0 since eλv>|x|-3 for small enough |x|.
In this case, for any u∈H01(Ω)∩L∞(Ω) and any λ>0, we have
∫Ωk2cosh(u+λv+w)dx≥12∫Ω2k2eu+λv+wdx≥kmin2e-∥u+w∥L∞(Ω2)2∫Ω2eλvdx=+∞.
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 k2cosh(u+w)∈L1(Ω).
Since J(u+λv)-J(u)≥0 for all v∈H01(Ω)∩L∞(Ω) and any λ≥0, we have
(2.5)12a(u+λv,u+λv)+∫Ωk2cosh(u+λv+w)dx-∫Ωl(u+λv)dx-12a(u,u)-∫Ωk2cosh(u+w)dx+∫Ωludx≥0.
Making equivalent transformations of (2.5) and dividing by λ>0, we obtain
(2.6)a(u,v)+limλ→0+∫Ωk2(cosh(u+λv+w)-cosh(u+w))λdx-∫Ωlvdx≥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):-k2(x)(cosh(u(x)+w(x)+λv(x))-cosh(u(x)+w(x)))λ→λ→0+k2(x)sinh(u(x)+w(x))v(x) for a.e.x∈Ω.
By the mean value theorem, we obtain
fλ(x)=k2(x)sinh(ς(x))v(x),
where ς(x):-u(x)+w(x)+Θ(x)λv(x) and Θ(x)∈(0,1), a.e. x∈Ω.
Then
|fλ(x)|≤k2(x)(eς(x)+e-ς(x)2)|v(x)|,
from which it follows that
(2.8)|fλ|≤∥v∥L∞(Ω)k2eu+w+e-u-w2e|Θ(x)|λ|v(x)|≤∥v∥L∞(Ω)e∥v∥L∞(Ω)k2cosh(u+w)⏟∈L1(Ω) 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 ∫Ωk2sinh(u+w)vdx, and therefore we obtain
(2.9)a(u,v)+∫Ωb(x,u+w)vdx-∫Ωlvdx≥0 for allv∈H01(Ω)∩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)dx≥0 for allv,z∈H01(Ω)∩L∞(Ω).
If u1,u2∈H01(Ω) are two different solutions of (2.1), then
(2.11)a(u1-u2,v)+∫Ω(b(x,u1+w)-b(x,u2+w))vdx=0 for allv∈H01(Ω)∩L∞(Ω).
Note that the difference u1-u2 is not necessarily in H01(Ω)∩L∞(Ω).
To show that we can test with u1-u2 in (2.11), we apply a property of Sobolev spaces proved in [3].
Theorem 2.2.
Let Ω be an open set in Rd, T∈H-1(Ω)∩Lloc1(Ω), and v∈H01(Ω).
If there exists a function f∈L1(Ω) such that T(x)v(x)≥f(x) a.e in Ω, then Tv∈L1(Ω) and the duality product 〈T,v〉 in H-1(Ω)×H01(Ω) coincides with ∫ΩTvdx.
We have the following situation: a locally summable function g∈Lloc1(Ω) defines a bounded linear functional Tg over the dense subspace C0∞(Ω) of H01(Ω) through the integral formula 〈Tg,φ〉=∫Ωgφdx.
It is clear that the functional Tg is uniquely extendable by continuity to a bounded linear functional T¯g over the whole space H01(Ω).
The question is whether this extension is still representable by the same integral formula for any v∈H01(Ω) (if the integral makes sense at all).
If the function v∈H01(Ω) 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〉=∫Ωgvdx.
Now, applying Theorem 2.2 to the functional Tg defined by
〈Tg,v〉:-〈b(x,u1+w)-b(x,u2+w),v〉 for allv∈H01(Ω)∩L∞(Ω)
and the function v=u1-u2∈H01(Ω), using (2.11), we conclude that
a(u1-u2,u1-u2)+∫Ω(b(x,u1+w)-b(x,u2+w))(u1-u2)dx=0.
Using the monotonicity (2.10) of b(x,⋅) and the coercivity of a(⋅,⋅), we obtain u1=u2.
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 s0≤t<∞ and such that
φ(h)≤Cφ(s)β(h-s)α for allh>s>s0,
where C and α are positive constants and β>1.
If j∈R is defined by jα:-Cφ(s0)β-12αββ-1, then φ(s0+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∥L2(Ω), ϵ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=Gs(u):-sgn(u)max{|u|-s,0},
where s≥∥w∥L∞(Ω2) (we notice that similar test functions Gs have been used in [16, Theorem B.2] in the context of linear elliptic problems).
It is easy to see that Gs(0)=0, this function is Lipschitz continuous and, therefore, Gs(u)∈H01(Ω) (e.g., see [16, 12]).
Next the functional Tb defined by
〈Tb,v〉:-∫Ωb(x,u+w)vdx for allv∈H01(Ω)∩L∞(Ω)
is bounded and linear and b(x,u+w)∈Lloc1(Ω).
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 〈Tb,Gs(u)〉=∫Ωb(x,u+w)Gs(u)dx, it suffices to verify the inequality
(2.13)b(x,u+w)Gs(u)≥f a.e. for somef∈L1(Ω).
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)Gs(u)={b(x,u+w)(u-s)≥0foru>s,0foru∈[-s,s],b(x,u+w)(u+s)≥0foru<-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)Gs(u)dx≥0.
Moreover, using the definition of a(⋅,⋅) and the definition (2.12) of Gs(u), we obtain
(2.16)a(u,Gs(u))=∫Ωϵ∇u⋅∇Gs(u)dx=∫Ωϵ∇Gs(u)⋅∇Gs(u)dx≥ϵmin∥∇Gs(u)∥L2(Ω)2≥ϵminCF2∥Gs(u)∥L2(Ω)2,
where CF is the constant in Friedrichs’ inequality ∥v∥L2(Ω)≤CF∥∇v∥L2(Ω) that holds for all v∈H01(Ω).
Finally, using (2.1), (2.15), and (2.16), we get
∥Gs(u)∥L2(Ω)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,Gs(u)) from below and ∫ΩlGs(u)dx from above using the Sobolev embedding H1(Ω)↪Lq(Ω), where q=∞ for d=1, q<∞ for d=2, and q=2dd-2 for d≥3.
Let q* denote the Hölder conjugate to q.
Then q*=1 for d=1, q*=qq-1>1 for d=2, and q*=2dd+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 CE we denote the embedding constant in the inequality ∥v∥L6(Ω)≤CE∥v∥H1(Ω) for all v∈H1(Ω), which depends only on the domain Ω and d.
For a(u,Gs(u)), we have
(2.17)a(u,Gs(u))=∫Ωϵ∇Gs(u)⋅∇Gs(u)dx≥ϵmin1+CF2∥Gs(u)∥H1(Ω)2
and, for ∫ΩlGs(u)dx, we obtain
(2.18)∫ΩlGs(u)dx=∫A(s)lGs(u)dx≤∥l∥Lq*(A(s))∥Gs(u)∥Lq(Ω)≤CE∥l∥Lq*(A(s))∥Gs(u)∥H1(Ω),
where A(s):-{x∈Ω:|u(x)|>s}.
Combining (2.17), (2.18), (2.15), and (2.1), we arrive at the estimate
(2.19)ϵmin1+CF2∥Gs(u)∥H1(Ω)≤CE∥l∥Lq*(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 H1(Ω)↪Lq(Ω) and Hölder’s inequality yields
(2.20)∥Gs(u)∥H1(Ω)≥1CE(∫Ω|Gs(u)|qdx)1q=1CE(∫A(s)||u|-s|qdx)1q≥1CE(∫A(h)(h-s)qdx)1q=1CE(h-s)|A(h)|1q
and
(2.21)∥l∥Lq*(A(s))≤∥l∥L2(Ω)|A(s)|2-q*2q*.
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)|≤(CE2(1+CF2)ϵmin∥l∥L2(Ω))q|A(s)|q-22(h-s)q for allh>s≥∥w∥L∞(Ω2).
Since q-22=2>1, applying Lemma 2.1, we conclude that there is some j>0 such that
0<jq=(CE2(1+CF2)ϵmin∥l∥L2(Ω))q|A(∥w∥L∞(Ω2))|q-422q(q-2)q-4≤(CE2(1+CF2)ϵmin∥l∥L2(Ω))q|Ω|q-422q(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∈H01(Ω)∩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∈H01(Ω) resulting in a standard weak formulation.
If k2 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, k2 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)findu∈Vsuch that J(u)=infv∈VJ(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 0V 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 allv∈V,y*∈Y*.
The functional J*:V*→ℝ¯ defined by the relation
J*(v*):-supv∈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*)findp*∈Y*such that I*(p*)=supy*∈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(0V) 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)=infv∈VJ(v)=infv∈Vsupy*∈Y*L(v,y*)=supy*∈Y*infv∈VL(v,y*)=supy*∈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 allv∈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*)=DG(Λv,y*)+DF(v,-Λ*y*)-:M⊕2(v,y*),
where
DG(Λv,y*):-G(Λv)+G*(y*)-(y*,Λv),
DF(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 DG and DF 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*)=DG(Λu,y*)+DF(u,-Λ*y*),
(3.2b)J(v)-I*(p*)=M⊕2(v,p*)=DG(Λv,p*)+DF(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)DG(Λv,p*)+DF(v,-Λ*p*)+DG(Λu,y*)+DF(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=H01(Ω), V*=H-1(Ω), Y=[L2(Ω)]d=Y*, G(Λv)=G(∇v)=∫Ω12A∇v⋅∇vdx and F(v)=∫Ω(12v2-lv)dx (where A is a symmetric positive definite matrix with bounded entries), then
(3.5)DG(Λv,p*)=12∫ΩA∇(v-u)⋅∇(v-u)dx,DG(Λu,y*)=12∫ΩA-1(y*-p*)⋅(y*-p*)dx,
DF(v,-Λ*p*)=12∥v-u∥L2(Ω)2,DF(u,-Λ*y*)=12∥div(y*-p*)∥L2(Ω)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)dx+∫ΩA-1(y*-p*)⋅(y*-p*)dx+∥v-u∥L2(Ω)2+∥div(y*-p*)∥L2(Ω)2=|||A∇v-y*|||*2+∥v-divy*-l∥L2(Ω)2=2M⊕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:-H01(Ω), Y:-[L2(Ω)]d, where d=2,3, and Λ the gradient operator ∇:H01(Ω)→[L2(Ω)]d.
We further denote g:Ω×ℝ3→ℝ, g(x,ξ):-ϵ(x)2|ξ|2, and B:Ω×ℝ→ℝ, B(x,ξ):-k2(x)cosh(ξ).
With this notation, we have
G(Λv):-∫Ωg(x,∇v(x))dx=∫Ωϵ2|∇v|2dx,
F(v):-∫ΩB(x,v(x)+w(x))dx-∫Ωlvdx=∫Ωk2cosh(v+w)dx-∫Ωlvdx,
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=log1|x|α, α≥d on the unit ball in ℝd).
However, if d≤2, then exp(v)∈L1(Ω) for all v∈H01(Ω) and F:V→ℝ (see [15]).
Also, F(0V) is obviously finite since w∈L∞(Ω2).
We set V*=H-1(Ω) and Y*=Y=[L2(Ω)]d.
In this case, Λ* coincides with -div considered as an operator from [L2(Ω)]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*)=∫Ω12ϵ(x)|y*(x)|2dx.
For y*∈H(div;Ω) and an arbitrary function z:Ω2→ℝ, we introduce the functional
Iy*(z):-∫Ω2[(divy*+l)z-B(x,z+w)]dx.
Recalling that the nonlinearity B is supported on Ω2, we have
(3.7)F*(-Λ*y*)=supz∈H01(Ω)[〈-Λ*y*,z〉-F(z)]=supz∈H01(Ω)[(-y*,Λz)-F(z)]=supz∈H01(Ω)∫Ω[-y*⋅∇z-B(x,z+w)+lz]dx (ify*∈H(div;Ω))=supz∈H01(Ω)∫Ω[divy*z-B(x,z+w)+lz]dx (finite ifdivy*+l=0inΩ1)=supz∈H01(Ω)Iy*(z)≤∫Ω2supξ∈ℝ[(divy*(x)+l(x))ξ-B(x,ξ+w(x))]dx=∫Ω2[(divy*(x)+l(x))ξ0(x)-B(x,ξ0(x)+w(x))]dx=Iy*(ξ0).
Here ξ0:Ω2→ℝ is computed by the condition
(3.8)ddξ[(divy*(x)+l(x))ξ-B(x,ξ+w(x))]=0 for a.e.x∈Ω2,
which is equivalent to
divy*(x)+l(x)-k2(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*)+ρk2(y*)+1)-w(x)=ln(Θ(ρk(y*)))-w(x),
where ρk(y*):-divy*(x)+l(x)k2(x) and Θ(s):-s+s2+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∈L2(Ω), we see that divp*=k2sinh(u+w)+l∈L2(Ω) and thus p*∈H(div;Ω).
In Proposition 3.1, we will later prove that we have not overestimated the supremum over z∈H01(Ω) 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))=x2+1 for allx∈ℝ,
for any y*∈H(div;Ω)⊂[L2(Ω)]d=Y* with divy*+l=0 in Ω1, we obtain an explicit formula for F*(-Λ*y*):
(3.9)F*(-Λ*y*)=∫Ω2[k2ρk(y*)(ln(Θ(ρk(y*)))-w)-k2ρk2(y*)+1]dx=∫Ω2[k2sinh(S)(S-w)-k2cosh(S)]dx.
Remark 3.1.
Since |ln(t+t2+1)|≤|t| for all t∈ℝ, the function ln(Θ(f(x)))-w(x) belongs to L2(Ω2) for any f∈L2(Ω2), and we conclude that ξ0(x)∈L2(Ω2) if y*∈H(div;Ω).
Therefore, the integral in (3.9) is well defined.
Now our goal is to prove that the inequality supz∈H01(Ω)Iy*(z)≤Iy*(ξ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 divy*+l=0 in Ω1, it holds
supz∈H01(Ω)Iy*(z)=Iy*(ξ0)<∞.
Proof.
The idea is to approximate f=divy*+lk2∈L2(Ω2) and w↾Ω2∈L∞(Ω2) by C0∞(Ω2) functions (in the a.e. sense) and use the Lebesgue dominated convergence theorem.
Let fn∈C0∞(Ω2) and wn∈C0∞(Ω2) be such that fn(x)→f(x) a.e. in Ω2, |fn(x)|≤h(x)∈L2(Ω2) (see [2, Theorem 4.9]), wn(x)→w(x) a.e. in Ω2, |wn(x)|≤m+2, where m:-∥w∥L∞(Ω2).
Then
zn(x):-ln(Θ(fn(x)))-wn(x)→ξ0(x) a.e. inΩ2
and zn∈C0∞(Ω2)⊂H01(Ω2)⊂H01(Ω) (by extending the functions by zero in Ω1).
Since B(x,⋅) is continuous, we have the pointwise a.e. in Ω2 convergence
(divy*(x)+l(x))zn(x)-B(x,zn+w(x))→(divy*(x)+l(x))ξ0(x)-B(x,ξ0(x)+w(x))
Now we search for a function in L1(Ω2) that majorates the function
|(divy*(x)+l(x))zn(x)-B(x,zn+w(x))|:
(3.10)|(divy*(x)+l(x))zn(x)-k2(x)cosh(zn(x)+w(x))|≤|divy*(x)+l(x)||zn(x)|+k2(x)e∥w∥L∞(Ω2)e|zn(x)|.
Our next goal is to bound |zn(x)| in (3.10).
For the first summand, we have
|zn(x)|=|ln(Θ(fn(x)))-wn(x)|≤|fn(x)|+m+2≤h(x)+m+2∈L2(Ω2),
where Remark 3.1 has been used.
However, this bound cannot be used in the second term because eh might not belong even to L1(Ω2).
In order to find an L1-majorant for the second summand in (3.10), we distinguish the following two cases:
In the first case, fn(x)>0.
Then |ln(Θ(fn(x)))|≤|ln(Θ(h(x)))|.
In the second case (fn(x)≤0), we have Θ(fn(x))≤1.
Therefore, 0≥fn(x)≥-h(x).
Since Θ(s) is a monotonically increasing function,
Θ(0)=1≥Θ(fn(x))≥Θ(-h(x))>0.
From here, we obtain
ln(1)=0≥ln(Θ(fn(x)))≥ln(Θ(-h(x))),
and using the relation Θ(-h)=1Θ(h), we conclude that
|ln(Θ(fn(x)))|≤|ln(Θ(-h(x)))|=|ln(Θ(h(x)))|.
Finally, for almost all x∈Ω2, we have
|zn(x)|=|ln(Θ(fn(x)))-wn(x)|≤|ln(Θ(h(x)))|+m+2=ln(Θ(h(x)))+m+2 because h(x)≥0 for a.e.x∈Ω2.
Therefore,
|(divy*(x)+l(x))zn(x)-k2(x)cosh(zn(x)+w(x))|≤|divy*(x)+l(x)|(h(x)+∥w∥L∞(Ω2)+2)+k2(x)e2∥w∥L∞(Ω2)+2Θ(h(x)):-H(x)∈L2(Ω2),
where, in the last line, we used the fact that Θ(h(x))∈L2(Ω2).
All the conditions of the Lebesgue dominated convergence theorem are satisfied, and we see that Iy*(zn)→Iy*(ξ0) and, consequently,
supz∈H01(Ω)Iy*(z)=Iy*(ξ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 divy*+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*)=DG(Λv,y*)+DF(v,-Λ*y*)=G(Λv)+G*(y*)-(y*,Λv)+F(v)+F*(-Λ*y*)+〈Λ*y*,v〉=∫Ωη2(x)dx=12|||ϵ∇v-y*|||*2+DF(v,-Λ*y*),
where
(3.12)η2(x)={12ϵ|ϵ∇v-y*|2,x∈Ω112ϵ|ϵ∇v-y*|2+k2cosh(v+w)-lv+k2ρk(y*)(ln(Θ(ρk(y*)))-w)-k2ρk2(y*)+1-divy*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),divy*(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^*)=12∫Ωϵ|∇(v-u)|2dx-:12|||∇(v-u)|||2,
DG(Λu,y*)=12∫Ω1ϵ|y*-p*|2dx-:12|||y*-p*|||*2.
Now we find explicit expressions for the nonlinear measures DF(v,-Λ*p*) and DF(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)22≤A(s,t)≤(sinh(t)-sinh(s))22,
where A(s,t)=cosh(t)-cosh(s)+ssinh(s)-tsinh(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
DF(v,-Λ*p*)=∫Ω2(k2cosh(v+w)-lv+k2sinh(u+w)u-k2cosh(u+w)-divp*v)dx=∫Ω2k2(cosh(v+w)-cosh(u+w)+usinh(u+w)-vsinh(u+w))dx.
Similarly, DF(u,-Λ*y*)=∫Ω2k2(cosh(T)-cosh(S)+Ssinh(S)-Tsinh(S))dx, where T:-arsinh(ρk(p*)).
The nonlinear quantities DF(v,-Λ*p*) and DF(u,-Λ*y*) measure the error in v and in divy*, respectively.
Using inequality (3.14), we can represent these two measures in a form which resembles the corresponding estimates in the case (3.5) of a quadratic functional F, namely,
(3.15)∫Ω2k22(v-u)2dx≤DF(v,-Λ*p*)≤∫Ω2k22(sinh(v+w)-sinh(u+w))2dx,
(3.16)∫Ω2k22(T-S)2dx≤DF(u,-Λ*y*)≤∫Ω212k2(divp*-divy*)2dx.
Note that, for k≥kmin>0 in Ω, the following equivalences hold:
∫Ωk22(v-u)2dx≂∥v-u∥L2(Ω)2 and ∫Ω12k2(divp*-divy*)2dx≂∥divy*-divp*∥L2(Ω)2.
Moreover, replacing the nonlinear term k2sinh(u+w) with u, inequalities (3.15)
and (3.16) reduce to the equalities for DF(v,-Λ*p*) and DF(u,-Λ*y*) in (3.5) because, in this case, the inverse function of f(x)=x is again f(x).
The functions on the left-hand side, in the middle, and on the right-hand side in inequality (3.14) are depicted on Figure 1.
Further, if v is in a δ1-neighborhood of u in L∞(Ω) norm, then we can find a constant C1(δ1,∥u∥L∞(Ω))>1 such that
(3.17)∫Ω2k22(sinh(v+w)-sinh(u+w))2dx≤C1(δ1,∥u∥L∞(Ω))∫Ω2k22(v-u)2dx.
Analogously, if l∈L∞(Ω2) and ∥div(y*-p*)∥L∞(Ω2)≤δ2 (recall that when l∈L∞(Ω2), divp* is in L∞(Ω2)), then we can find a constant C2(δ2,∥divp*∥L∞(Ω2))<1 such that
(3.18)C2(δ2,∥divp*∥L∞(Ω2))∫Ω212k2(divp*-divy*)2dx≤∫Ω2k22(T-S)2dx.
The constants C1 and C2 are just Lipschitz constants for the locally Lipschitz function sinh.
Notice that if k2≥kmin>0 in Ω, then everywhere in (3.15), (3.16), (3.17), and (3.18), the integrals are taken over the entire domain Ω.
Now the abstract error identity (3.4) takes the form
(3.19)12|||∇(u-v)|||2+12|||p*-y*|||*2+∫Ω2k22(v-u)2dx+C2(δ2,∥divp*∥L∞(Ω))∫Ω212k2(divp*-divy*)2dx≤12|||∇(u-v)|||2+12|||p*-y*|||*2+DF(v,-Λ*p*)+DF(u,-Λ*y*)=M⊕2(v,y*)¯≤12|||∇(u-v)|||2+12|||p*-y*|||*2+C1(δ1,∥u∥L∞(Ω))∫Ω2k22(v-u)2dx+∫Ω212k2(divy*-divp*)2dx,
where we have used p*=ϵΛu=ϵ∇u.
Relation (3.19) shows that the computable majorant M⊕2(v,y*) is bounded from below and above by a multiple of one and the same error norm.
Since DF(v,-Λ*p*)≥0 and DF(u,-Λ*y*)≥0, we also obtain a guaranteed bound on the error in the combined energy norm,
(3.20)|||∇(u-v)|||2+|||p*-y*|||*2≤2M⊕2(v,y*).
From the pointwise equality
(3.21)1ϵ|ϵ∇v-y*|2=1ϵ|ϵ∇(v-u)-(y*-p*)|2=ϵ|∇(v-u)|2+1ϵ|y*-p*|2-2(y*-p*)⋅∇(v-u),
after applying Young’s inequality and integrating over Ω, we obtain a lower bound for the error in combined energy norm,
(3.22)12|||ϵ∇v-y*|||*2≤|||∇(v-u)|||2+|||y*-p*|||*2
Remark 3.2.
Integrating (3.21) over Ω, we obtain the algebraic identity
(3.23)|||ϵ∇v-y*|||*2=|||∇(v-u)|||2+|||y*-p*|||*2-2∫Ω(y*-p*)⋅∇(v-u)dx,
from which the Prager–Synge identity is derived.
Comparing the last relation with (3.19), using the fact that M⊕(v,y*)2=12|||ϵ∇v-y*|||*2+DF(v,-Λ*y*), we arrive at the relation
(3.24)DF(v,-Λ*y*)=DF(v,-Λ*p*)+DF(u,-Λ*y*)+∫Ω(y*-p*)⋅∇(v-u)dx.
From here, it is seen that if the integral on the right-hand side is small compared to the other terms, then the error in v and divy* measured with DF(v,-Λ*p*)+DF(u,-Λ*y*) is controlled mainly by the computable term DF(v,-Λ*y*) in the majorant M⊕2(v,y*).
Moreover, (3.23)
enables us to give a practical estimation of the error in combined energy norm, which is very close to the real error in all of the experiments that we have conducted.
We conclude the section by presenting a near best approximation result.
Contrary to the result in [5, Theorem 6.2], we do not make any restrictive assumptions on the meshes to ensure that the finite element approximations uh are uniformly bounded in L∞ norm.
In our analysis, Vh⊂L∞ is a finite-dimensional subspace of H01, and uh is the minimizer of J over Vh, which is the unique solution of the Galerkin problem:
(3.25)finduh∈Vhsuch that a(uh,v)+∫Ωb(x,uh+w)vdx=(l,v) for allv∈Vh.
Then, using (3.2b) and expression (3.13) for DG(Λv,p*), for any v∈Vh, we can write
|||∇(uh-u)|||2+2DF(uh,-Λ*p*)=2(J(uh)-J(u))≤2(J(v)-J(u))=|||∇(v-u)|||2+2DF(v,-Λ*p*).
Since 2DF(uh,-Λ*p*)≥0, using (3.15), we obtain the following generalization of Cea’s lemma to the case of our nonlinear problem.
Proposition 3.3.
Let Vh⊂L∞(Ω) be a closed subspace of H01(Ω) and uh∈Vh the Galerkin approximation of u defined by (3.25).
Then
(3.26)|||∇(uh-u)|||2≤infv∈Vh{|||∇(v-u)|||2+∫Ω2k2(sinh(v+w)-sinh(u+w))2dx}.
Since we use the finite element method with P1 Lagrange elements, let Vh be the corresponding space, where h refers to the maximum element size.
By Ih(φ), we denote the Lagrange finite element interpolant of φ∈C0(Ω).
Using (3.26), we can show unqualified convergence of the finite element approximations uh to u as h→0.
Let ε>0, and u¯∈C0∞(Ω) is such that ∥∇(u¯-u)∥L2(Ω)≤ε and ∥u¯∥L∞(Ω)≤∥u∥L∞(Ω)+2.
Also, let L be the Lipschitz constant in the inequality
|sinh(s)-sinh(t)|≤L|s-t| for alls,t∈[-2∥w∥L∞(Ω2)-j¯-2,2∥w∥L∞(Ω2)+j¯+2],
where j¯ is the constant from Proposition 2.1.
Then, applying the triangle inequality together with Young’s inequality, we obtain
(3.27)|||∇(uh-u)|||2≤2(|||∇(Ih(u¯)-u¯)|||2+|||∇(u¯-u)|||2)+2(∫Ωk2(sinh(Ih(u¯)+w)-sinh(u¯+w))2dx+∫Ωk2(sinh(u¯+w)-sinh(u+w))2dx).
For the first term in (3.27), assuming mesh regularity, we have
|||∇(Ih(u¯)-u¯)|||2+|||∇(u¯-u)|||2≤ϵmax(C|u¯|22h2+ε2),
where |u¯|2 denotes the H2 seminorm of u¯ and C>0 is a constant depending on the mesh regularity.
Using the fact that ∥Ih(u¯)∥L∞(Ω)≤∥u¯∥L∞(Ω)≤∥u∥L∞(Ω)+2, for the second term in (3.27), we obtain the upper bound
2kmax2L2(∥Ih(u¯)-u¯∥L2(Ω)2+∥u¯-u∥L2(Ω)2)≤2kmax2L2CF2(∥∇(Ih(u¯)-u¯)∥L2(Ω)2+∥∇(u¯-u)∥L2(Ω)2)≤2kmax2L2CF2(C|u¯|22h2+ε2).
This inequality shows that the right-hand side of (3.27) can be made as small as desired provided that we choose ε and h small enough, and therefore |||∇(uh-u)|||→0 when h→0.
Moreover, (3.26) can be also used to obtain qualified convergence of uh in the energy norm under additional assumptions on the interface Γ, the meshes, and the regularity of u.
4 Numerical Results
In this section, we present numerical examples illustrating error identity (3.19) and performance of functional a posteriori error estimates.
All numerical experiments are carried out in FreeFem++ developed and maintained by Frederich Hecht [13], and all pictures are generated in VisIt [6].
We solve adaptively the homogeneous nonlinear problem (1.1) with w:-whref=g-zhref, where zhref is a good Galerkin finite element approximation of the solution z of
-∇⋅(ϵ∇z)=-k2sinh(g)+linΩ1∪Ω2,
[z]Γ=0,
[ϵ∂z∂n]Γ=0,
z=0on∂Ω,
for given functions g and l.
We compare the accuracy of the adaptively computed solution uh of (1.1) for w=whref to the reference solution zhref.
The adaptive mesh refinement is based on the error indicator ∥2η∥L2(Oi) on subdomains Oi, where the function η is defined in (3.12) and η2 is the integrand of the majorant M⊕2(v,y*).
The factor 2 accounts for the factor 2 in (3.20).
More precisely, we find approximations uh to the exact solution u∈H01(Ω) of the problem
∫Ωϵ∇u⋅∇vdx+∫Ωb(x,u+whref)vdx=∫Ωlvdx=0 for allv∈H01(Ω).
In all examples, we use piecewise constant parameters ϵ and k, and for y*∈H(div;Ω), we used a patchwise equilibrated reconstruction of the numerical flux ϵ∇uh based on [1].
More precisely, we find y* in the Raviart–Thomas space RT0 over the same mesh such that its divergence is equal to the L2 orthogonal projection of k2sinh(uh+w)+l onto the space of piecewise constants.
Recall that
M⊕2(v,y*)=M⊕2(v,p*)+M⊕2(u,y*),
where M⊕2(v,y*)=12|||ϵ∇v-y*|||*2+DF(v,-Λ*y*) and M⊕2(v,p*)=J(v)-J(u)=12|||∇(v-u)|||2+DF(v,-Λ*p*) is the primal error, whereas M⊕2(u,y*)=I*(p*)-I*(y*)=12|||y*-p*|||*2+DF(u,-Λ*y*) is the dual error.
Further, we use v for the approximate solution uh and u for the reference solution zhref and define the efficiency index of the lower bound for the error in combined energy norm (3.22) by
IEffCEN,Low:-22|||ϵ∇v-y*|||*|||∇(v-u)|||2+|||y*-p*|||*2.
Similarly,
IEffCEN,Up:-2M⊕2(v,y*)|||∇(v-u)|||2+|||y*-p*|||*2
defines the efficiency index of the upper bound (3.20) for the error in combined energy norm.
Finally,
IEffE:-2M⊕2(v,y*)|||∇(v-u)||| and PrelCEN:-|||ϵ∇v-y*|||*|||∇v|||2+|||y*|||*2
define the efficiency index of the upper bound for the error in energy norm and the practical estimate of the relative error in combined energy norm, respectively.
4.1 Example 1 (2D Problem)
In the first example, the domain Ω is a square with a side 20 with Ω1 being a regular 15-sided polygon with a radius of its circumscribed circle equal to 2.
The coefficients ϵ and k are defined by the relations
ϵ(x)={ϵ1=1,x∈Ω1,ϵ2=100,x∈Ω2, k(x)={k1=0.15,x∈Ω1,k2=0.4,x∈Ω2,
and
g=L(exp(-b1((x1-c1)2σ12-1))-exp(-b2((x2-c2)2σ22-1))),
l=0, where b1=2=b2=2, c1=-1, c2=6, σ1=σ2=1.5, L=0.8.
The reference solution zhref is computed on a multiply refined mesh with 50 086 142 triangles.
Note that k2=0.0225 in Ω1 and k2=0.16 in Ω2.
The mesh adaptation is done with the built-in function “adaptmesh” of FreeFem++.
The localized error indicator ∥2η∥L2(Oi), computed on each vertex patch Oi of the mesh, is compared to its average value over all patches, and the local mesh size is divided by two if this average is smaller than the local value.
Table 1 illustrates the main error identity (3.3) and the convergence of its constituent parts.
Further, it is seen that the dual error 2M⊕2(u,y*) dominates the primal error in this example.
This is due to the fact that the term 2DF(u,-Λ*y*), measuring the error in divy* (cf. (3.16) and (3.18)), is much larger than |||∇(v-u)|||2+DF(v,-Λ*p*), where DF(v,-Λ*p*) behaves like ∥v-u∥L2(Ω2)2 (cf. (3.15) and (3.17)).
As we mentioned earlier, for y*, we use a partially equilibrated reconstruction of the numerical flux ϵ∇v, which is the reason why the integral term in (3.23) is negligible compared to the combined energy norm of the error.
This fact is confirmed by the values of the efficiency index of the lower bound (3.22).
Table 1
Constituent parts of main error identity (3.3) for Example 1 (2D).
| Example 1 (2D): k1=0.15, k2=0.4, ϵ1=1, ϵ2=100 |
| # elts | ∥v-u∥0∥u∥0 [%] | |||∇(v-u)||||||∇u||| [%] | |||y*-p*|||*|||p*|||* [%] | 2M⊕2(v,y*) | 2M⊕2(v,p*) | 2M⊕2(u,y*) |
| 196 | 15.0077 | 51.5582 | 86.1021 | 1778.14 | 66.5980 | 1711.54 |
| 347 | 5.69339 | 30.8534 | 41.7241 | 703.594 | 20.7780 | 682.816 |
| 630 | 4.20384 | 21.7715 | 31.4858 | 217.719 | 10.2201 | 207.498 |
| 1315 | 2.39552 | 15.8532 | 23.1244 | 76.8018 | 5.37574 | 71.4261 |
| 2865 | 1.87075 | 11.7353 | 17.1655 | 33.9310 | 2.94414 | 30.9869 |
| 5938 | 0.64611 | 7.93001 | 11.4692 | 16.0812 | 1.33874 | 14.7425 |
| 12006 | 0.36985 | 5.64786 | 8.23544 | 7.75232 | 0.67872 | 7.07360 |
| 24571 | 0.16023 | 3.94241 | 5.76054 | 3.85268 | 0.33039 | 3.52229 |
| 48483 | 0.08909 | 2.80265 | 4.09366 | 1.90043 | 0.16682 | 1.73361 |
| 97423 | 0.03961 | 1.97875 | 2.88455 | 0.96275 | 0.08304 | 0.87970 |
| 192905 | 0.02230 | 1.39832 | 2.03200 | 0.47524 | 0.04136 | 0.43388 |
| 386185 | 0.01015 | 0.99471 | 1.44616 | 0.24134 | 0.02082 | 0.22052 |
Table 2
Constituent parts of error identity (3.6) for example Example 1 (2D).
| Example 1 (2D): k1=0.15, k2=0.4, ϵ1=1, ϵ2=100 |
| # elts | |||∇(v-u)|||2 | |||y*-p*|||*2 | 2DF(v,-Λ*p*) | 2DF(u,-Λ*y*) |
| 196 | 56.5057 | 157.588 | 10.0923 | 1553.95 |
| 347 | 20.2350 | 37.0058 | 0.54296 | 645.811 |
| 630 | 10.0756 | 21.0729 | 0.14450 | 186.425 |
| 1315 | 5.34235 | 11.3668 | 0.03338 | 60.0593 |
| 2865 | 2.92742 | 6.26338 | 0.01671 | 24.7235 |
| 5938 | 1.33673 | 2.79619 | 0.00200 | 11.9462 |
| 12006 | 0.67805 | 1.44169 | 0.00067 | 5.63191 |
| 24571 | 0.33038 | 0.70538 | 0.00001 | 2.81691 |
| 48483 | 0.16696 | 0.35622 | 0.00000 | 1.37739 |
| 97423 | 0.08323 | 0.17687 | 0.00000 | 0.70283 |
| 192905 | 0.04156 | 0.08777 | 0.00000 | 0.34611 |
| 386185 | 0.02103 | 0.04445 | 0.00000 | 0.17606 |
In Table 3, we can see that IEffCEN,Low is approximately equal to 0.7071≈22.
The value of the efficiency index with respect to the combined energy norm and the value of the ratio DF(v,-Λ*y*)/M⊕2(v,y*) are also coupled in the sense that if we have only one of these two quantities, we can estimate the other one by using the main error equality (3.19).
This estimation is accurate because the integral term in (3.24) is very close to zero, and therefore DF(v,-Λ*y*)≈DF(v-Λ*p*)+DF(u-Λ*y*).
One more consequence of using a partially equilibrated flux is that we obtain a very accurate practical estimate of the absolute and relative error in combined energy norm as illustrated in the last two columns of Table 3.
Table 3
Efficiency indices for Example 1 (2D).
| Example 1 (2D): k1=0.15, k2=0.4, ϵ1=1, ϵ2=100 |
| # elts | DF(v,-Λ*y*)M⊕2(v,y*) [%] | IEffCEN,Low | IEffCEN,Up | IEffE,Up | PrelCEN [%] | True rel. error in CEN [%] |
| 196 | 89.0701 | 0.67371 | 2.88191 | 5.60966 | 74.6973 | 70.9641 |
| 347 | 92.4942 | 0.67919 | 3.50597 | 5.89671 | 36.2638 | 36.6935 |
| 630 | 85.9525 | 0.70066 | 2.64380 | 4.64848 | 27.1574 | 27.0680 |
| 1315 | 78.2616 | 0.70681 | 2.14392 | 3.79158 | 19.9383 | 19.8250 |
| 2865 | 72.8992 | 0.70729 | 1.92142 | 3.40452 | 14.7523 | 14.7032 |
| 5938 | 74.3009 | 0.70708 | 1.97256 | 3.46846 | 9.87419 | 9.85973 |
| 12006 | 72.6473 | 0.70722 | 1.91238 | 3.38130 | 7.06762 | 7.06119 |
| 24571 | 73.1176 | 0.70708 | 1.92864 | 3.41485 | 4.93753 | 4.93591 |
| 48483 | 72.4826 | 0.70694 | 1.90588 | 3.37371 | 3.50789 | 3.50805 |
| 97423 | 73.0084 | 0.70678 | 1.92392 | 3.40108 | 2.47256 | 2.47347 |
| 192905 | 72.8486 | 0.70629 | 1.91692 | 3.38145 | 1.74226 | 1.74418 |
| 386185 | 72.9912 | 0.70546 | 1.91972 | 3.38748 | 1.23829 | 1.24114 |
Figure 2 depicts a mesh that is a part of a sequence of meshes obtained by mesh adaptation using the localized functional error indicator ∥2η∥L2(Oi).
Figure 3 depicts a mesh with approximately the same number of elements but obtained by mesh adaptation using the error indicator |||ϵ∇v-y*|||*(Oi).
The mesh in Figure 2 is refined mainly where the error in divy* is the dominant part of the error M⊕2(v,-Λ*p*)+M⊕2(u,-Λ*y*).
On the other hand, the mesh in Figure 3 is refined most around the extrema of the solution.
Figure 5 depicts the minimal set of elements K of a mesh Th that contains at least 30 % of the total indicated error ∑K∈Th|||ϵ∇v-y*|||*(K) (greedy algorithm with a bulk factor of 0.3), where Th is part of the same sequence as the mesh illustrated in Figure 3.
Figure 6 depicts the elements marked by the greedy algorithm using a bulk factor of 0.5 and employing the true error
2M⊕2(v,p*)+2M⊕2(u,y*)
as indicator.
Figure 7 depicts elements which are marked additionally or fail to be marked by the same greedy algorithm when employing the functional error indicator ∥2η∥L2(Oi) for the same bulk factor.
The ratio of the number of these differently marked elements, that is, elements which are marked by one of the two methods but not by the other one, and the total number of elements is 0.022, and the ratio of the number of differently marked elements to the number of marked elements using the true error is 0.048 (Table 4).
Comparing the indicated error and the true error elementwise, one finds that the error indicator generated by the majorant M⊕2(v,y*) reproduces the local distribution of the error with a very high accuracy.
This is also confirmed by Figure 8, where it can be seen that all error measures are almost identical in both cases of adaptive mesh refinement.
Mesh adaptation based on the functional error indicator ∥2η∥L2(Oi) instead of the error indicator |||ϵ∇v-y*|||*(Oi) (Figure 9) yields approximately twice smaller efficiency indices in energy and combined energy norms and approximately twice smaller values for the full error M⊕2(v,p*)+M⊕2(u,y*) on meshes with a comparable number of elements.
The reason for the higher efficiency indices is that no adaptive control is applied on the nonlinear part of the error measure in (3.19), and consequently, the ratio DF(v,-Λ*y*)/M⊕2(v,y*) is increasing, reaching values close to 100 % on fine meshes.
However, the error in |||∇(v-u)||| and |||y*-p*|||* might be a little higher in the case of the functional error indicator ∥2η∥L2(Oi).
For example, on the mesh from Figure 5 with 24 122 elements, M⊕2(v,p*)+M⊕2(u,y*)=3.8314, |||∇(v-u)|||=0.4674, |||y*-p*|||*=0.6540, whereas on a mesh with 24 571 elements from the sequence adapted with the indicator ∥2η∥L2(Oi), we obtained a value of 1.9263 for M⊕2(v,y*), and 0.574791 and 0.8399 for |||∇(v-u)||| and |||y*-p*|||*, respectively.
This shows that, reducing the error in divy*, the functional error indicator ∥2η∥L2(Oi) provides a better approximation for the primal and dual problem together.
Table 4
Marking based on true error and functional error indicator in Example 1 (2D).
| Example 1 (2D): k1=0.15, k2=0.4, ϵ1=1, ϵ2=100 |
| # elts | # marked elts with true error | # differently marked elts | differently marked elts in % of all mesh elts |
| 196 | 62 | 6 | 3.06122 |
| 347 | 150 | 10 | 2.88184 |
| 630 | 288 | 14 | 2.22222 |
| 1315 | 632 | 39 | 2.96578 |
| 2865 | 1439 | 113 | 3.94415 |
| 5938 | 2949 | 216 | 3.63759 |
| 12006 | 5981 | 534 | 4.44778 |
| 24571 | 12099 | 961 | 3.91111 |
| 48483 | 24194 | 2233 | 4.60574 |
| 97423 | 47784 | 4012 | 4.11812 |
Now we want to demonstrate that flux equilibration is indeed an important subtask to make the proposed error bounds reliable and efficient.
For this purpose, we use a simple global gradient averaging procedure, i.e., project the numerical flux ϵ∇v∈L2(Ω) onto the subspace [Vh]2, where Vh is the finite element space of continuous piecewise linear functions.
Then the problem in Example 1 is solved adaptively once applying the functional error indicator ∥2η∥L2(Oi) and next applying the error indicator |||ϵ∇v-y*|||*(Oi).
Figure 10 shows the adapted mesh with 563 965 elements, which is a part of a sequence of meshes obtained applying the functional error indicator with gradient averaging for y*.
Figure 11 shows a mesh with 444 092 elements, which is part of a sequence of meshes adapted using the second indicator with gradient averaging for y*.
Comparing with the results based on flux equilibration for y*, it can be seen that the mesh in Ω2 close to the interface Γ is refined too much for both error indicators.
Apart from that, the meshes on Figures 11 and 3 look quite similar, unlike the meshes on Figures 10 and 2.
For meshes with a comparable number of elements, applying the indicator |||ϵ∇v-y*|||*(Oi) using gradient averaging instead of flux equilibration, we obtained ∼ 30 % larger values for the error |||∇(v-u)||| and 60 % larger values for the error |||y*-p*|||*.
The difference in the errors when applying the functional indicator ∥2η∥L2(Oi) with gradient averaging for y* instead of flux equilibration resulted in an even more drastic increase of the error, namely, between 40 % and 180 % for |||∇(v-u)||| and between 64 % and 66 % for |||y*-p*|||*, where the meshes had between 21 528 and 563 965 elements.
In both cases, we obtained an increasing sequence of efficiency indices with respect to energy and combined energy norms reaching values of 133 and 107 with the functional error indicator on a mesh with 2 089 022 elements, and 570 and 269 with the error indicator |||ϵ∇v-y*|||*(Oi) on a mesh with 2 954 218 elements.
This is due to the fact that the nonlinear term DF(u,-Λ*y*), which measures the error in divy* (see (3.16) and (3.18)), dominates the other terms in the nonlinear measure M⊕2(v,p*)+M⊕2(u,y*) for the error, reaching more than 99.99 % of it in both cases.
In both experiments with gradient averaging for y*, increasing values of DF(u,-Λ*y*) are in correspondence with increasing error ∥divy*-divp*∥L2(Ω) and increasing efficiency indices.
4.2 Example 2 (2D Problem)
Figures 13 and 15 show how meshes depend on the indicator in another example, where ϵ1=1ϵ2=100, k1=0.2, k2=0.3.
The functions
g=exp(-b1(|x-c1|2σ12-1))-exp(-b2(|x-c2|2σ22-1)) and l=exp(-b3(|x|2σ32-1))sin(x1x24),
where b1=2.2, b2=2.5, b3=6, c1=(-1,0), c2=(5,5), σ1=σ2=2, σ3=10.
The indicator |||ϵ∇v-y*|||*(Oi) correctly approximates elementwise errors in the combined energy norm but does not capture the rest of the error, which results from the nonlinearity k2sinh(u+w) and the right-hand side l in (1.1).
On the other hand, the term DF(v,-Λ*y*) controls the error DF(v,-Λ*p*)+DF(u,-Λ*y*), and this is the reason why the mesh on Figure 13 resembles the wavy features of the function f=-k2sinh(u+w)+l.
The isolines of the reference solution and of the function f are depicted on Figures 14 and 12.
4.3 Example 3 (3D Problem)
Here we consider an example close to a real physical problem.
The computational domain Ω is a cube of side length 20 angstrom with a triangulated water molecule Ω1 in it.
The diameter of the water molecule, which is positioned in the center of the cube, is about 2.75 angstrom.
Its shape is not changed during the mesh adaptation process.
The surface mesh of the water molecule is taken from [28].
Figure 16 illustrates the initial tetrahedral mesh, which consists of 60 222 elements.
It is generated using TetGen [25] and adaptively refined with the help of mmg3d [7].
Using the localized error indicator ∥2η∥L2(Oi) computed on each vertex patch Oi of the mesh, a new local mesh size at each vertex is defined by the formula
hinew=hiold(max{min{AM{∥2η∥L2(Oj)}∥2η∥L2(Oi),1},0.35})
and supplied to mmg3d, where AM{∥2η∥L2(Oj)} is the arithmetic mean of {∥2η∥L2(Oj)} over all vertex patches Oj.
The coefficients ϵ and k for this example are typical for electrostatic computations in biophysics using the PBE and are given by
ϵ(x)={ϵ1=2,x∈Ω1,ϵ2=80,x∈Ω2, k(x)={k1=0,x∈Ω1,k2=0.84,x∈Ω2.
We consider the homogeneous problem, i.e., l=0, and
g=exp(-b1(|x-c1|2σ12-1))-exp(-b2(|x-c2|2σ22-1))+exp(-b3(|x-c3|2σ32-1))+exp(-b4(|x-c4|2σ42-1)),
where
b1=b2=b3=b4=2.3,
c1=(1,1,0), c2=(4,4,0), c2=(0,6,0), c2=(-5,0,0),
σ1=σ2=σ3=σ4=2.
The reference solution zhref is computed on a very fine mesh, obtained after a sequence of adaptive mesh refinements, that contains 79 917 007 tetrahedrons.
Since l=0 in Ω1 is a constant function, the patchwise reconstruction from [1] produces a flux y* with zero divergence in Ω1 and, therefore, the reliability of our majorant is guaranteed.
In this example, we achieved a very tight guaranteed bound on the error in combined energy norm, as well as in energy norm.
The efficiency index IEffCEN,Up settles at around 1.05 and the efficiency index IEffE,Up decreases to 1.30 (Table 6).
This is in a good agreement with the fact that, in this example, the ratio DF(v,-Λ*y*)/M⊕2(v,y*) is well controlled and decreases to around 10 % (column 2 in Table 7).
We also note that, in this example, we obtained very similar results with the error indicator |||ϵ∇v-y*|||*(Oi).
For the efficiency index IEffCEN,Low of the lower bound on the combined energy norm of the error, we obtain values converging to approximately 0.7071, which is the approximate value of 22 (column 3 in Table 6).
This means that the combined energy norm of the error
|||∇(v-u)|||2+|||y*-p*|||*2
is practically equal to |||ϵ∇v-y*|||*.
Table 5
Constituent parts of main error identity (3.3) for Example 3 (3D).
| Example 3 (3D): k1=0, k2=0.84, ϵ1=2, ϵ2=80 |
| # elts | ∥v-u∥0∥u∥0 [%] | |||∇(v-u)||||||∇u||| [%] | |||y*-p*|||*|||p*|||* [%] | 2M⊕2(v,y*) | 2M⊕2(v,p*) | 2M⊕2(u,y*) |
| 60222 | 76.8320 | 108.015 | 167.589 | 425569 | 117373 | 308196 |
| 103236 | 11.9257 | 46.3306 | 55.1210 | 47104.5 | 17845.0 | 29259.5 |
| 222118 | 1.09233 | 17.7353 | 14.9578 | 4484.44 | 2224.69 | 2259.75 |
| 552936 | 0.49820 | 8.67222 | 7.09062 | 965.067 | 513.706 | 451.361 |
| 1.17360e+06 | 0.25609 | 6.58075 | 5.33661 | 539.734 | 295.254 | 244.481 |
| 2.05668e+06 | 0.17094 | 5.37625 | 4.18207 | 350.648 | 197.016 | 153.631 |
| 2.97315e+06 | 0.12317 | 4.73466 | 3.53852 | 265.167 | 152.783 | 112.385 |
| 3.90692e+06 | 0.10071 | 4.32886 | 3.12966 | 216.336 | 127.703 | 88.6336 |
Table 6
Constituent parts of error identity (3.6) for Example 3 (3D).
| Example 3 (3D): k1=0, k2=0.84, ϵ1=2, ϵ2=80 |
| # elts | |||∇(v-u)|||2 | |||y*-p*|||*2 | 2DF(v,-Λ*p*) | 2DF(u,-Λ*y*) | REUp [%] |
| 60222 | 79487.0 | 191346 | 37886.0 | 116850 | — |
| 103236 | 14623.9 | 20699.7 | 3221.12 | 8559.78 | 310.049 |
| 222118 | 2142.92 | 1524.28 | 81.7757 | 735.474 | 33.9219 |
| 552936 | 512.376 | 342.528 | 1.32980 | 108.833 | 13.4714 |
| 1.17360e+06 | 295.039 | 194.026 | 0.21458 | 50.455 | 9.75647 |
| 2.05668e+06 | 196.919 | 119.155 | 0.09743 | 34.4762 | 7.72193 |
| 2.97315e+06 | 152.724 | 85.3044 | 0.05857 | 27.0805 | 6.64970 |
| 3.90692e+06 | 127.666 | 66.7303 | 0.03663 | 21.9033 | 5.96873 |
Table 7
Efficiency indices for Example 3 (3D).
| Example 3 (3D): k1=0, k2=0.84, ϵ1=2, ϵ2=80 |
| # elts | DF(v,-Λ*y*)M⊕2(v,y*) [%] | IEffCEN,Low | IEffCEN,Up | IEffE,Up | PrelCEN [%] | True rel. error in CEN [%] |
| 60222 | 40.0541 | 0.68627 | 1.25353 | 2.31386 | 92.8434 | 140.985 |
| 103236 | 20.4500 | 0.72828 | 1.15478 | 1.79473 | 47.6870 | 50.9159 |
| 222118 | 16.1172 | 0.71615 | 1.10583 | 1.44661 | 16.4040 | 16.4054 |
| 552936 | 11.2249 | 0.70786 | 1.06248 | 1.37241 | 7.90966 | 7.92099 |
| 1.17360e+06 | 9.33477 | 0.70731 | 1.05053 | 1.35254 | 5.98505 | 5.99106 |
| 2.05668e+06 | 9.82289 | 0.70725 | 1.05327 | 1.33442 | 4.81343 | 4.81632 |
| 2.97315e+06 | 10.2057 | 0.70722 | 1.05547 | 1.31767 | 4.17784 | 4.17960 |
| 3.90692e+06 | 10.1194 | 0.70719 | 1.05492 | 1.30175 | 3.77592 | 3.77716 |
Another consequence of this fact is the good accuracy of the practical estimation PrelCEN of the relative error in combined energy norm (columns 6 and 7 in Table 6).
The tight bounds on the error also enable us to compute tight and guaranteed upper bounds on the relative error in energy norm:
(4.1)|||∇(v-u)||||||∇u|||≤2M⊕2(v,y*)|||∇v|||-2M⊕2(v,y*)-:REUp,
where (4.1) is valid if |||∇v|||-2M⊕2(v,y*)>0.
As a remark, we note that the efficiency indices with respect to the energy and combined energy norms of the error can be improved if we use a flux reconstruction in a bigger space, say, RT1, which has better approximation properties.
In this way, the error in divy* will decrease and, as a result, the term DF(v,-Λ*y*) and consequently the dual part of the error 2M⊕2(u,y*)=|||y*-p*|||*2+DF(u,-Λ*y*) will constitute a smaller part of the whole majorant and the error, respectively.
Even better, we can minimize the majorant with respect to y* in a subspace of H(div;Ω) like RT0, possibly on another finer mesh.
Note that, in the limit case, we have infy*∈H(div;Ω)M⊕2(v,y*)=M⊕2(v,p*)=12|||∇(v-u)|||2+DF(v,-Λ*p*), and the dual error completely vanishes.
In this case,
IEffCEN,Up=IEffE=2M⊕2(v,p*)|||∇(v-u)|||=|||∇(v-u)|||2+2DF(v,-Λ*p*)|||∇(v-u)|||,
where the last ratio tends to 1 because, by (3.15) and (3.17), the term DF(v,-Λ*p*)∼∥v-u∥L2(Ω2)2 and has a higher order of convergence than |||∇(v-u)|||2.
In practice, we can minimize the majorant with respect to y* only once on a sufficiently big subspace of H(div;Ω) to find some good approximation y¯* of p* and then reuse this y¯* and obtain guaranteed and tight bounds on the error in energy and combined energy norm.
A key factor that determines the efficiency index is the ratio DF(v,-Λ*y*)M⊕2(v,y*).
Assuming that
DF(v,-Λ*y*)≈DF(v,-Λ*p*)+DF(u,-Λ*y*),
which means that the last term in (3.24) is close to zero, we obtain from (3.19) the estimate
IEffCEN,Up≈11-DF(v,-Λ*y*)M⊕2(v,y*).
From what we have observed, the efficiency index IEffE,Up with respect to the energy norm usually is no more than twice bigger than IEffCEN,Up (assuming we have a good approximation y* to p*).
Therefore, if, during the computations, we detect that this ratio is increasing, we can apply the so-called estimation with one step delay, i.e., compute the value of the majorant M⊕2(v,y*) for the current mesh level with the reconstructed y* from the next level.
