A Posteriori Modelling-Discretization Error Estimate for Elliptic Problems with L ^\infty -Coefficients

We consider elliptic problems with complicated, discontinuous diffusion tensor $A_{\scriptscriptstyle 0} $. One of the standard approaches to numerically treat such problems is to simplify the coefficient by some approximation, say $A_{\varepsilon}$, and to use standard finite elements. In \cite{Repin2012} a combined modelling-discretization strategy has been proposed which estimates the discretization and modelling errors by a posteriori estimates of functional type. This strategy allows to balance these two errors in a problem adapted way. However, the estimate of the modelling error is derived under the assumption that the difference $A_{\scriptscriptstyle 0} -A_{\varepsilon}$ is bounded in the $L^{\infty}$-norm, which requires that the approximation of the coefficient matches the discontinuities of the original coefficient. Therefore this theory is not appropriate for applications with discontinuous coefficients along \textit{complicated, curved} interfaces. Based on bounds for $A_{\scriptscriptstyle 0} -A_{\varepsilon}$ in an $L^{q}$-norm with $q<\infty$ we generalize the combined modelling-discretization strategy to a larger class of coefficients.


Introduction
We consider elliptic boundary value problems with complicated, discontinuous diffusion tensor.As a model problem we choose the diffusion equation − div(A 0 ∇u) = f in a two-or three-dimensional bounded domain Ω with homogeneous Dirichlet boundary conditions and A 0 ∈ L ∞ (Ω, R d×d sym ) is symmetric and positive definite.Our emphasis is on diffusion matrices A 0 containing a large number of different scales which we allow to be highly non-uniformly distributed over the domain.
It is well-known that for such problems standard single scale numerical methods such as standard finite element methods are not efficient, since one needs to solve the problem on a sufficiently fine mesh which resolves all the fine-scale behavior of the coefficient.This is usually too costly, especially for three-dimensional problems.Essentially there are two different approaches to overcome this difficulty: One is to design (non-polynomial) generalized finite element methods where the characteristic behavior of the solution is reflected by the shape of the basis functions.This approach has been investigated by many researchers (see, e.g., [3], [4], [2], [15], [21][22][23]).In this paper we follow the second approach which tries to simplify the diffusion coefficient by some approximation and then employs standard finite elements.Standard methods for simplifying the coefficients are based, e.g., on homogenization methods for periodic structures (see.e.g., [13], [10], [8]), or on different upscaling techniques e.g.[11], [17], [19].
In many applications a numerical solution with only moderate guaranteed accuracy is required.For such problems, the combined modelling-discretization strategy has been proposed in [19].This approach consists of two basic steps.In a first step the diffusion coefficient A 0 is replaced by a simpler coefficient A ε and the simplified model is discretized and solved on a rather coarse mesh.In the second step the discretization and modelling errors are controlled using some a posteriori estimates.The total error is bounded by the sum of the discretization and modelling errors which are both explicit and computable.If the total error is larger than a given tolerance, then either the mesh should be refined (if the discretization error dominates) or the coefficient has to be modelled more accurately (if the modelling error dominates).Thus the discretization and modelling errors can be balanced in a problem-adapted way.
The error estimates in [19] are derived by purely functional methods without requiring specific information on the approximating subspace and the numerical method used.Consequently the estimates contain no mesh dependent constants and are valid for any conforming approximation from the respective energy space.
However, the modelling errors arising due to the simplification of the coefficients contain the term |||A 0 − A ε ||| ∞,Ω (cf.[19]) and one has to assume that the approximation A ε matches the discontinuities of A 0 in order to ensure that the term |||A 0 − A ε ||| ∞,Ω becomes small.Since in many applications the discontinuities jump on curved or cracked interfaces, they cannot be captured exactly by the finite element mesh.Hence, this smallness assumption is not suitable for the analysis of numerical methods for problems with discontinuous coefficients along complicated, curved interfaces.
This problem has been addressed in [9] in the context of adaptive finite element methods.Based on a perturbation theory the term |||A 0 − A ε ||| ∞,Ω is replaced by |||A 0 − A ε ||| q,Ω with q := 2p/(p − 2) for some p ≥ 2. The advantage of this approach is that A ε does not have to match the discontinuities of A 0 exactly.However, one needs more regularity on the solution, namely ∇u ∈ L p (Ω) for some p > 2. This also requires additional assumptions on the righthand side f .These requirements are quite mild and are satisfied in many applications.
The goal of our paper is to generalize the modelling-discretization strategy developed in [19] to a larger class of coefficients.Based on the theory presented in [9] which is based on results by [16] we bound the modelling error by a term depending on |||A 0 − A ε ||| q,Ω for some q > 2. Consequently the assumption on A ε can be weakened and the strategy can also be applied to problems where A 0 has discontinuities which are unknown or lie along curves and surfaces.Moreover we have guaranteed, computable upper bounds of the error.
Due to our new approach the bound of the modelling error also depends on a regularity constant C reg,Aε with respect to some Sobolev space W k,p (Ω) as it appears in the inequality for some p > 2, where u ε is the exact solution of the simplified problem and F is a linear functional generated by the right-hand side of the equation.The constant C reg,Aε depends only on p, A ε and Ω.In [9,16], it is shown (by perturbation arguments) that C reg,Aε can be expressed in terms of the constant C reg,I which corresponds to the Laplace operator.Therefore, we need to derive computable upper bounds for C p := C reg,I .
The paper is structured as follows.In Section 2, we first formulate the model problem and the conditions on the coefficient.Then, in Section 3 we present new error estimates for the modelling-discretization strategy introduced in [19].These estimates are based on the theory developed in [9,16] and require that ∇u ε ∈ L p (Ω) for some p > 2. Section 4 is devoted to present some L p -bounds for the gradient of the solution of diffusion problems with L ∞coefficient.These bounds only depend on the size of the jumps in the coefficient.Finally in Section 5 we present an explicit computable estimate of C p for the full space problem.This bound depends on p and the dimension d.

Notation and Problem Statement
Throughout the paper it is assumed that Ω is a bounded domain in R d (d = 2, 3) with C 1 boundary.By •, • we denote the usual Euclidean scalar product on R d .For 1 ≤ p ≤ ∞, • ℓ p denotes the discrete ℓ p -norm in R d .If p = 2, then we use • instead of • ℓ 2 .By (•, •) we denote the L 2 -scalar product.The Sobolev space of real-valued functions in L 2 (Ω) with gradients in L 2 (Ω) is denoted by H 1 (Ω) (and • 1,2,Ω is the respective norm).A subspace of H 1 (Ω) containing the functions vanishing on the boundary is denoted by H 1 0 (Ω).For any p ∈ [1, ∞], the adjoint number p ′ is defined by the relation Throughout the paper • p,Ω denotes the norm of L p (Ω).We use the usual notation W 1,p (Ω) for the Sobolev spaces of functions, which generalized derivatives belong to the space L p (Ω).This space is supplied with the standard norm • 1,p,Ω .The space of functions denoted by W 1,p 0 (Ω) is the closure of C ∞ 0 (Ω) with respect to the norm • 1,p,Ω .We also use the space W −1,p (Ω) := (W 1,p ′ 0 (Ω)) ′ endowed with the standard dual norm • −1,p,Ω .For vector and matrix valued functions, we use the same notation for the Lebesgue and Sobolev spaces as well as for the corresponding norms.To explicitly indicate the dimension we write L p Ω, R d and L p Ω, R d×d .For functions in L 2 Ω, R d we set Notice that for p = 2 we have p ′ = 2 and p ′′ = ∞ so that Also we use the space which is a Hilbert space endowed with the scalar product (y, z) div := (y, z) + (div y, div z) and the norm y div := (y, y) 1/2 div .For the functions in L 2 (Ω, R d ), we also introduce two equivalent norms (associated with the energy and complementary energy) where the matrix A 0 ∈ L ∞ Ω, R d×d sym is assumed to be uniformly positive definite, i.e.
Let f be a given function in L 2 (Ω).Consider the following boundary value problem: In view of (2.2), existence and uniqueness of the solution u follows from the Lax-Milgram lemma.
We consider problems with complicated matrix A 0 (which coefficients are complicated functions of x).In this case, direct approximation of u based upon standard numerical approaches may lead to high computational costs.One way to obtain a reasonable approximation of u with minimal expenditures is to consider a simplified problem with a simpler matrix A ε .If the difference between u and the respective solution u ε is explicitly estimated and it is smaller than the desired tolerance, then we can use the simplified problem instead of the problem (2.3).This idea leads to a set of simplified problems generated by uniformly positive definite matrices A ε ∈ L ∞ (Ω, R d×d sym ), where ε is a sequence of positive decreasing numbers (which is either finite, or infinite tending to zero).Henceforth, we assume that there exist positive constants α and β such that for any ε and that the collection of simplified matrices A ε satisfies the condition

Discretization and Combined Error Majorant
The function u ε ∈ H 1 0 (Ω) (generalized solution of the simplified problem) is defined by the integral identity The difference between u and u ε is the modelling error In general u ε is unknown and instead we use a conforming approximation u ε,h ∈ H 1 0 (Ω) computed by some numerical method.In view of this, we must also consider the discretization error At this point, we do not specify the method by which u ε,h is found.In the framework of our approach it is not important because the difference between u ε and u ε,h can be estimated within the framework of a unified method that follows from a posteriori error estimates of the functional type (see e.g.[18] and the references therein).In our case, the respective estimate has the form The majorant M 2 Ω (u ε,h , y, γ) contains a vector-valued function y ∈ H(Ω, div) and an arbitrary positive parameter γ.The constant depends on the geometry of Ω, namely, where In order to formulate the main result, we introduce the quantities Here θ(r, t) := 2(t−r) r(t−2) (for 2 < r < t < ∞) and C reg,Aε is a constant in the inequality In Section 4, we show that C reg,Aε depends only on the constant C P := C reg,I (associated with the Laplace operator) and on the amplitude of the jumps in the coefficient A ε (cf.Theorem 4.4).
2), f ∈ L P (Ω) for some P ∈ (2, +∞), and p ∈ (2, p * ), where the function p * = p * α β , P is defined by (4.4).Then where , and p ′′ is a number defined in Sect. 2. Proof.By the triangle inequality First we estimate the discretization error.It holds The last norm can be estimated by the error majorant and thus we obtain (3.3).Next we will estimate the modelling part of the error.Observe that We choose v = u − u ε and obtain Applying the Cauchy-Schwarz inequality yields Hence, We set p ′′ = p/(p − 2), apply the triangle and generalized Hölder inequalities, and obtain Now, we choose t ∈ (p, p * ( α β , P )).By Lemma A.1 and (3.2), we find that t) .Now, we apply Theorem 4.4 (where the regularity constant C reg,Aε is defined) and arrive at the estimate Notice that for f ∈ L t (Ω) we have F −1,t,Ω ≤ f t,Ω .The combination of (3.5) and (3.6) yields the desired estimate.
Theorem 3.1 requires several comments.First, from Theorem 3.1 it follows that Here the quantity Υ(u ε,h , f, θ(p, t)) is fully computable provided that C reg,Aε or a certain upper bound of it is known and norms of B ε and D ε can be computed a priori.Since these matrices have low order (typically d = 2 or 3), the computation of norms is reduced to well-known algebraic procedures.
It is easy to see that if A 0 = A ε , then D ε = I and B ε = 0.In this case, the second and the third terms in the right-hand side vanish and the total error is completely determined by the discretization error encompassed in the first term.Another limit case arises if u ε,h coincides with the exact solution of the problem (3.1).Then, the first two terms can be made arbitrary small and the overall error is determined by the modelling error encompassed in the last term.In other words, the first two terms can be made (at least theoretically) arbitrary small if h tends to zero.
Next, it is worth noticing that a somewhat different modus operandi leads to a simpler upper bound of the error.Indeed, in view of (3.5) and Theorem 4.4 we have Hence, we obtain another estimate which contains only two terms associated with the discretization and modelling errors, respectively.Here, the second term associated with the modelling error depends only on ε.This majorant may be coarser than (3.7), but it allows us to make an a priori estimation of the modelling error and decide whether or not the problem (3.1) could be used for getting an approximation with the desired accuracy δ.This fact suggests the strategy of solving a simple (smoothened, averaged) problem (3.1) instead of the complicated problem (2.3).For example, if we know that the modelling error is smaller than 1 2 δ, then we can concentrate on approximations of u ε instead of u.It is natural to await that in many cases u ε will have better regularity properties than u.Then, u ε,h will converge to u ε (as h → 0) much faster than analogous approximations u h in (2.3) will converge to u (it is known that this convergence for problems with complicated coefficients may be arbitrary slow, see [5]).Hence, we obtain an efficient method of getting an approximation with the required accuracy.Moreover, in this way proper reconstructions of the vector function y can be performed by various "gradient recovery" or "gradient averaging" methods, which has been investigated by many researchers (see, e.g., [6,7,14,25]).Thus, for small h the approximation errors can be controlled by the majorant M Ω .
Above observations motivate the idea to replace the L ∞ -norm of B ε by some L p -norm.We are interested to make the norm |||B ε ||| p ′′ ,Ω small and proportional to a certain positive power of ε.For example, let A 0 = κ 0 I and A ε = κ ε I, where κ 0 is a jump function and κ ε is a continuous piecewise affine function approximating this jump in the ε-strip.Then elementary computations show that |||B ε ||| p ′′ ,Ω ∼ ε 1/p ′′ .Similar a priori analysis is of course possible for more complicated matrices A 0 and A ε .

W 1,p -Regularity Results for Second Order Elliptic Problems with Rough Coefficients
Now our goal is to derive L p (Ω)-regularity estimates for the gradient of ∇u for some p > 2.
We start from the Poisson problem (in this case, A 0 is equal to the identity matrix I).Then, we employ perturbation arguments in order to get the desired estimates for a uniformly elliptic matrix A 0 ∈ L ∞ (Ω, R d×d sym ).It is important, that our estimates depend only on the amplitude of jumps in the coefficients of A 0 .
Consider the following problem: For a given 20]).Let 1 < p < ∞.Then, for every F ∈ W −1,p (Ω), the problem (4.1) has a unique solution ψ ∈ W 1,p 0 (Ω) which meets the estimate with the Laplace W 1,p -regularity constant C p and Remark 4.2.The constant C p is independent of F (and ψ) but depends on Ω, d and p.We have C 2 = 1 and, for p > 2, C p is non-decreasing and continuous in p (cf. [16]).
Let T := {(p, P ) : 2 < P < ∞, 2 ≤ p ≤ P } and introduce the function  (cf. Figure 1).The function p * has the explicit representation Writing A 0 as a perturbation of the identity and deducing the L p -bound of ∇u from the L p -bound for the solution of the Poisson problem (4.1) we obtain the following result for a uniformly elliptic 2), F ∈ W −1,P (Ω) for 2 < P < ∞ and 2 ≤ p < p * α 0 β 0 , P .Then the solution of (2.3) exists in W 1,p 0 (Ω) and meets the estimate For a proof we refer to [9,16,24].

(large perturbations) For any
Then, for any 2 ≤ p < p max it holds Proof.@1 The condition Hence, we may choose any 2 ≤ p < P for the following.It is easy to see that then η (p, P ) ≤ 1 and (4.5) gives us This leads to the first assertion.@2 The condition 0 implies that for any c ∈ (0, 1) the number p max satisfies p max ≤ p * α 0 β 0 , P , and Theorem 4.4 implies ∇u ∈ L p (Ω) for 2 ≤ p < p max .Note that then η (p, P ) ≤ For the regularity constant we obtain and the assertion follows.
1 By a Taylor argument it follows that the right-hand side in (4.7) is bounded from below by the expression which shows the qualitative dependence on c, α0 β0 , C P and P better.
5 Analysis of the Laplace W For f ∈ L p (Ω), 1 < p < ∞, the Newtonian potential of f is defined as where G is the fundamental solution of Laplace's equation which is given by where Γ (•) denotes the Gamma function.For fixed i, j we define the linear operator T : Proof.We follow the proof of [12,Theorem 9.9] and track the dependence of the constants on p and d.We first extend f to vanish outside Ω and fix a cube K 0 ⊃ Ω so that for fixed t > 0 we have By bisection of the edges of K 0 , we subdivide K 0 into 2 d congruent subcubes with disjoint interiors.Those subcubes K which satisfy K f ≤ t|K| are similarly subdivided and the process is repeated indefinitely.In this way we obtain a sequence of parallel subcubes (K ℓ ) ∞ ℓ=1 such that The function f is now split into a "good part" g defined by From [12, p. 232] we know that Let us now fix some ℓ and a sequence b ℓm ⊂ C ∞ 0 (K ℓ ) converging to b ℓ in L 2 (Ω) and satisfying Then for x / ∈ K ℓ we have the relation where ȳℓ denotes the center of K ℓ .Further for x / ∈ K ℓ and y ∈ K ℓ it holds for some point ζ y between y and ȳℓ .Some computations show that where δ ij denotes the Kronecker symbol.This leads to the estimate where δ ℓ := diam K ℓ .The combination of (5.6) and (5.7) yields It is easy to see that for Thus, We set F * := ℓ B δ ℓ (ȳ ℓ ) and G * := K 0 \F * .Letting m → ∞ and summing over ℓ we get By Lemma 5.1 and using the last estimate we obtain  Next we employ the Riesz-Thorin interpolation theorem to remove the singular behavior of C (d, p) as p → 2. Note that T : L p (Ω) → L p (Ω) is continuous for p ∈ {p 0 , p 1 } with p 0 = 3/2 and p 1 = 2 and T : L p 0 (Ω) + L p 1 (Ω) → L p 0 (Ω) + L p 1 (Ω) (observe that L p 0 (Ω) + L p 1 (Ω) = L p 0 (Ω) since Ω is bounded).We know that T f p,Ω ≤ C (d, p) f p,Ω for p ∈ {p 0 , p 1 } .

Further for d ≥ 2 Notation 5 . 2 .
and 1 < p < 2 we define the constants C(d) := 2 d+2 + 2 d+1 d (d + 5) + 2π dFor r > 0 and x ∈ R d , we denote the ball with radius r around x by B r (x) := y ∈ R d : y − x < r .The dimension d is clear from the argument x.We write short B r for B r (0).Lemma 5.3.Let Ω ⊂ R d be a bounded domain, f ∈ L 2 (Ω) and t > 0. Then we have the estimate µ T f (t) ≤ C(d) f 1,Ω t with C(d) as in (5.1).
and a "bad part" b = f − g.Clearly, |g| ≤ 2 d t a.e., b(x) = 0 for x ∈ G and K ℓ b = 0 for ℓ = 1, 2, . . .Since T is linear we have T f = T g + T b and thus

Figure 2 :
Figure 2: The constant C 1 (d, p) as a function of p for d = 2 (left) and for d = 3 (right).