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

TL;DR

This paper extends a posteriori error estimation methods for elliptic PDEs with discontinuous coefficients, allowing for more general coefficient approximations in complex geometries, improving error control in numerical simulations.

Contribution

It generalizes existing error estimates to handle coefficients with complex discontinuities by using $L^q$-norm bounds, broadening applicability.

Findings

Extended error estimates to $L^q$-norm bounds for coefficients.

Applicable to problems with curved, complicated interfaces.

Improved error control for numerical solutions of elliptic problems.

Abstract

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.

Equations

$$\Vert \nabla u_{\epsilon }{\Vert }_{p,\Omega }\le C_{\mathrm{reg},A_{\epsilon }}\Vert F{\Vert }_{-1,p,\Omega }$$

$$\Vert \cdot {\Vert }_{p,\Omega }:=\Vert \Vert \cdot {\Vert }_{{\ell }^p}{\Vert }_{p,\Omega }.$$

$$m:=\underset{\zeta \in {\mathbb{R}}^d\backslash \left\{0\right\}}{\sup }\frac{{\Vert M\left(\cdot \right)\zeta \Vert }_{{\ell }^{p^{\prime }}}}{{\Vert \zeta \Vert }_{{\ell }^p}}\text{and}{\left|\left|\left|M\right|\right|\right|}_{p^{\prime \prime },\Omega }={\Vert m\Vert }_{p^{\prime \prime },\Omega }.$$

$${\left|\left|\left|M\right|\right|\right|}_{\infty ,\Omega }=\underset{x\in \Omega }{\mathrm{esssup}}\left(\underset{\zeta \in {\mathbb{R}}^d\backslash \left\{0\right\}}{\sup }\frac{\Vert M\left(x\right)\zeta \Vert }{\Vert \zeta \Vert }\right).$$

$$H(\Omega ,\mathrm{div}):=\{y\in L^2(\Omega ,{\mathbb{R}}^d)\mid \mathrm{div}y\in L^2(\Omega )\},$$

$$(y,z{)}_{\mathrm{div}}:=\left(y,z\right)+\left(\mathrm{div}y,\mathrm{div}z\right)$$

$$\Vert y{\Vert }_{A_0}^2:=\left(A_{0}y,y\right)=\underset{\Omega }{\int }⟨A_{0}y,y⟩\text{and}\Vert y{\Vert }_{A_0^{-1}}^2:=\left(A_0^{-1}y,y\right),$$

$$0<{\alpha }_0:={\left|\left|\left|A_0^{-1}\right|\right|\right|}_{\infty ,\Omega }^{-1}\le {\left|\left|\left|A_0\right|\right|\right|}_{\infty ,\Omega }=:{\beta }_0<\infty .$$

$${\int }_{\Omega }⟨A_0\nabla u,\nabla v⟩={\int }_{\Omega }fv=:F(v)\forall v\in H_0^1\left(\Omega \right).$$

$$0<\underline{\alpha }\le {\alpha }_{\epsilon }:={\left|\left|\left|A_{\epsilon }^{-1}\right|\right|\right|}_{\infty ,\Omega }^{-1}\le {\left|\left|\left|A_{\epsilon }\right|\right|\right|}_{\infty ,\Omega }=:{\beta }_{\epsilon }\le \overline{\beta }<\infty$$

$$\Vert A_0-A_{\epsilon }{\Vert }_{q,\Omega }\le \epsilon .$$

$${\int }_{\Omega }⟨A_{\epsilon }\nabla u_{\epsilon },\nabla v⟩=F(v)\forall v\in H_0^1\left(\Omega \right).$$

$$E_{\mathrm{mod}}^{\epsilon }:=\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}.$$

$$E_{\mathrm{disc}}^{\epsilon ,h}:=\Vert \nabla (u_{\epsilon }-u_{\epsilon ,h}){\Vert }_{A_0}.$$

$$\Vert \nabla (u_{\epsilon }-u_{\epsilon ,h}){\Vert }_{A_{\epsilon }}^2$$

$$:=(1+\gamma )\Vert A_{\epsilon }\nabla u_{\epsilon ,h}-y{\Vert }_{A_{\epsilon }^{-1}}^2+\left(1+\frac{1}{\gamma }\right)C_{\Omega }^2\Vert \mathrm{div}y+f{\Vert }_{2,\Omega }^2.$$

$$C_{F\Omega }:=\underset{w\in H_0^1(\Omega )\backslash \{0\}}{\sup }\frac{\Vert w{\Vert }_{2,\Omega }}{\Vert \nabla w{\Vert }_{2,\Omega }}\le \frac{\mathrm{diam}\Omega }{\sqrt{2}\pi }.$$

$${\rm Y}(u_{\epsilon ,h},f,\theta (p,t)):={\left(C_{\mathrm{reg},A_{\epsilon }}\Vert f{\Vert }_{t,\Omega }+\Vert \nabla u_{\epsilon ,h}{\Vert }_{t,\Omega }\right)}^{1-\theta (p,t)}{\alpha }_{\epsilon }^{-\frac{1}{2}\theta (p,t)}$$

$$\Theta (u_{\epsilon ,h}):={\left({\mathcal{M}}_{\Omega }(u_{\epsilon ,h},y,\gamma )\right)}^{\theta (p,t)}{\rm Y}(u_{\epsilon ,h},f,\theta (p,t))+\Vert \nabla u_{\epsilon ,h}{\Vert }_{p,\Omega }.$$

$$\Vert \nabla u_{\epsilon }{\Vert }_{t,\Omega }\le C_{\mathrm{reg},A_{\epsilon }}\Vert F{\Vert }_{-1,t,\Omega }.$$

$$\Vert \nabla (u-u_{\epsilon ,h}){\Vert }_{A_0}\le E_{\mathrm{disc}}^{\epsilon ,h}+E_{\mathrm{mod}}^{\epsilon },$$

$$E_{\mathrm{disc}}^{\epsilon ,h}\le {\left|\left|\left|D_{\epsilon }\right|\right|\right|}_{\infty ,\Omega }^{1/2}{\mathcal{M}}_{\Omega }(u_{\epsilon ,h},y,\gamma )$$

$$E_{\mathrm{mod}}^{\epsilon }\le {\left|\left|\left|B_{\epsilon }\right|\right|\right|}_{p^{\prime \prime },\Omega }^{1/2}\Theta (u_{\epsilon ,h}),$$

$$\Vert \nabla (u-u_{\epsilon ,h}){\Vert }_{A_0}\le \Vert \nabla (u_{\epsilon }-u_{\epsilon ,h}){\Vert }_{A_0}+\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}=E_{\mathrm{disc}}^{\epsilon ,h}+E_{\mathrm{mod}}^{\epsilon }.$$

$${\left(E_{\mathrm{disc}}^{\epsilon ,h}\right)}^2$$

$$={\int }_{\Omega }⟨D_{\epsilon }{A}_{\epsilon }^{1/2}\nabla (u_{\epsilon }-u_{\epsilon ,h}),A_{\epsilon }^{1/2}\nabla (u_{\epsilon }-u_{\epsilon ,h})⟩$$

$$\le {\left|\left|\left|D_{\epsilon }\right|\right|\right|}_{\infty ,\Omega }{\int }_{\Omega }⟨A_{\epsilon }\nabla (u_{\epsilon }-u_{\epsilon ,h}),\nabla (u_{\epsilon }-u_{\epsilon ,h})⟩$$

$$={\left|\left|\left|D_{\epsilon }\right|\right|\right|}_{\infty ,\Omega }\Vert \nabla (u_{\epsilon }-u_{\epsilon ,h}){\Vert }_{A_{\epsilon }}^2.$$

$$0=\underset{\Omega }{\int }⟨A_0\nabla (u-u_{\epsilon }),\nabla v⟩+\underset{\Omega }{\int }⟨(A_0-A_{\epsilon })\nabla u_{\epsilon },\nabla v⟩\forall v\in H_0^1(\Omega ).$$

$${\left(E_{\mathrm{mod}}^{\epsilon }\right)}^2=\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}^2$$

$$\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}^2\le {\left({\int }_{\Omega }⟨B_{\epsilon }\nabla u_{\epsilon },\nabla u_{\epsilon }⟩\right)}^{1/2}\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}.$$

$$\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}^2\le {\int }_{\Omega }⟨B_{\epsilon }\nabla u_{\epsilon },\nabla u_{\epsilon }⟩.$$

$$\Vert \nabla (u-u_{\epsilon }){\Vert }_{A_0}$$

Full text

1 Introduction

We consider elliptic boundary value problems with complicated, discontinuous diffusion tensor. As a model problem we choose the diffusion equation $-\operatorname{div}(A_{\scriptscriptstyle 0}\nabla u)=f$ in a two- or three-dimensional bounded domain $\Omega$ with homogeneous Dirichlet boundary conditions and $A_{\scriptscriptstyle 0}\in L^{\infty}(\Omega,\mathbb{R}_{\operatorname*{sym}}^{d\times d})$ is symmetric and positive definite. Our emphasis is on diffusion matrices $A_{\scriptscriptstyle 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_{\scriptscriptstyle 0}$ is replaced by a simpler coefficient $A_{\varepsilon}$ 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 $\left|\left|\left|A_{\scriptscriptstyle 0}-A_{\varepsilon}\right|\right|\right|_{\infty,\Omega}$ (cf. [19]) and one has to assume that the approximation $A_{\varepsilon}$ matches the discontinuities of $A_{\scriptscriptstyle 0}$ in order to ensure that the term $\left|\left|\left|A_{\scriptscriptstyle 0}-A_{\varepsilon}\right|\right|\right|_{\infty,\Omega}$ 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 $\left|\left|\left|A_{\scriptscriptstyle 0}-A_{\varepsilon}\right|\right|\right|_{\infty,\Omega}$ is replaced by $\left|\left|\left|A_{\scriptscriptstyle 0}-A_{\varepsilon}\right|\right|\right|_{q,\Omega}$ with $q:=2p/(p-2)$ for some $p\geq 2$. The advantage of this approach is that $A_{\varepsilon}$ does not have to match the discontinuities of $A_{\scriptscriptstyle 0}$ exactly. However, one needs more regularity on the solution, namely $\nabla u\in L^{p}(\Omega)$ for some $p>2$. This also requires additional assumptions on the right-hand 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 $\left|\left|\left|A_{\scriptscriptstyle 0}-A_{\varepsilon}\right|\right|\right|_{q,\Omega}$ for some $q>2$. Consequently the assumption on $A_{\varepsilon}$ can be weakened and the strategy can also be applied to problems where $A_{\scriptscriptstyle 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_{\operatorname*{reg},A_{\varepsilon}}$ with respect to some Sobolev space $W^{k,p}\left(\Omega\right)$ as it appears in the inequality

[TABLE]

for some $p>2$, where $u_{\varepsilon}$ 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_{\operatorname*{reg},A_{\varepsilon}}$ depends only on $p$, $A_{\varepsilon}$ and $\Omega$. In [16, 9], it is shown (by perturbation arguments) that $C_{\operatorname*{reg},A_{\varepsilon}}$ can be expressed in terms of the constant $C_{\operatorname*{reg},I}$ which corresponds to the Laplace operator. Therefore, we need to derive computable upper bounds for $C_{p}:=C_{\operatorname*{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 [16, 9] and require that $\nabla u_{\varepsilon}\in L^{p}(\Omega)$ 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^{\infty}$-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$.

2 Notation and Problem Statement

Throughout the paper it is assumed that $\Omega$ is a bounded domain in $\mathbb{R}^{d}$ ($d=2,3$) with $C^{1}$ boundary. By $\langle\cdot,\cdot\rangle$ we denote the usual Euclidean scalar product on $\mathbb{R}^{d}$. For $1\leq p\leq\infty$, $\|\cdot\|_{\ell^{p}}$ denotes the discrete $\ell^{p}$-norm in $\mathbb{R}^{d}$. If $p=2$, then we use $\left\|\cdot\right\|$ instead of $\left\|\cdot\right\|_{\ell^{2}}$. By $\left(\cdot,\cdot\right)$ we denote the $L^{2}$-scalar product. The Sobolev space of real-valued functions in $L^{2}\left(\Omega\right)$ with gradients in $L^{2}(\Omega)$ is denoted by $H^{1}(\Omega)$ (and $\|\cdot\|_{1,2,\Omega}$ is the respective norm). A subspace of $H^{1}(\Omega)$ containing the functions vanishing on the boundary is denoted by $H_{0}^{1}(\Omega)$. For any $p\in[1,\infty]$, the adjoint number $p^{\prime}$ is defined by the relation $\frac{1}{p}+\frac{1}{p^{\prime}}=1$. Analogously, for $p\in\left[2,\infty\right]$, the number $p^{\prime\prime}\in\left[1,\infty\right]$ satisfies the relation $\frac{2}{p}+\frac{1}{p^{\prime\prime}}=1$.

Throughout the paper $\|\cdot\|_{p,\Omega}$ denotes the norm of $L^{p}(\Omega)$. We use the usual notation $W^{1,p}(\Omega)$ for the Sobolev spaces of functions, which generalized derivatives belong to the space $L^{p}(\Omega)$. This space is supplied with the standard norm $\|\cdot\|_{1,p,\Omega}$. The space of functions denoted by $W_{0}^{1,p}(\Omega)$ is the closure of $C_{0}^{\infty}(\Omega)$ with respect to the norm $\|\cdot\|_{1,p,\Omega}$. We also use the space $W^{-1,p}\left(\Omega\right):=(W_{0}^{1,p^{\prime}}\left(\Omega\right))^{\prime}$ endowed with the standard dual norm $\|\cdot\|_{-1,p,\Omega}$. 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}\left(\Omega,\mathbb{R}^{d}\right)$ and $L^{p}\left(\Omega,\mathbb{R}^{d\times d}\right)$. For functions in $L^{2}\left(\Omega,\mathbb{R}^{d}\right)$ we set

[TABLE]

For $M\in L^{\infty}\left(\Omega,\mathbb{R}_{\operatorname{sym}}^{d\times d}\right)$ and $p\geq 2$, we introduce the function $m\in L^{\infty}\left(\Omega\right)$ by

[TABLE]

Notice that for $p=2$ we have $p^{\prime}=2$ and $p^{\prime\prime}=\infty$ so that

[TABLE]

Also we use the space

[TABLE]

which is a Hilbert space endowed with the scalar product

[TABLE]

and the norm $\|y\|_{\operatorname{div}}:=(y,y)_{\operatorname{div}}^{1/2}$. For the functions in $L^{2}(\Omega,\mathbb{R}^{d})$, we also introduce two equivalent norms (associated with the energy and complementary energy)

[TABLE]

where the matrix $A_{\scriptscriptstyle 0}\in L^{\infty}\left(\Omega,\mathbb{R}_{\operatorname*{sym}}^{d\times d}\right)$ is assumed to be uniformly positive definite, i.e.

[TABLE]

Let $f$ be a given function in $L^{2}\left(\Omega\right)$. Consider the following boundary value problem: Find $u\in H_{0}^{1}\left(\Omega\right)$ such that

[TABLE]

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_{\scriptscriptstyle 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_{\varepsilon}$. If the difference between $u$ and the respective solution $u_{\varepsilon}$ 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_{\varepsilon}\in L^{\infty}(\Omega,\mathbb{R}_{\operatorname*{sym}}^{d\times d})$, where $\varepsilon$ is a sequence of positive decreasing numbers (which is either finite, or infinite tending to zero). Henceforth, we assume that there exist positive constants $\underline{\alpha}$ and $\overline{\beta}$ such that for any $\varepsilon$

[TABLE]

and that the collection of simplified matrices $A_{\varepsilon}$ satisfies the condition

[TABLE]

3 Discretization and Combined Error Majorant

The function $u_{\varepsilon}\in H_{0}^{1}\left(\Omega\right)$ (generalized solution of the simplified problem) is defined by the integral identity

[TABLE]

The difference between $u$ and $u_{\varepsilon}$ is the modelling error

[TABLE]

In general $u_{\varepsilon}$ is unknown and instead we use a conforming approximation $u_{\varepsilon,h}\in H_{0}^{1}\left(\Omega\right)$ computed by some numerical method. In view of this, we must also consider the discretization error

[TABLE]

At this point, we do not specify the method by which $u_{\varepsilon,h}$ is found. In the framework of our approach it is not important because the difference between $u_{\varepsilon}$ and $u_{\varepsilon,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

[TABLE]

The majorant $\mathcal{M}_{\Omega}^{2}(u_{\varepsilon,h},y,\gamma)$ contains a vector-valued function $y\in H(\Omega,\operatorname{div})$ and an arbitrary positive parameter $\gamma$. The constant depends on the geometry of $\Omega$, namely, $C_{\Omega}:=C_{F\Omega}/\sqrt{\underline{\alpha}}$, where

[TABLE]

In order to formulate the main result, we introduce the quantities

[TABLE]

and

[TABLE]

Here $\theta(r,t):=\frac{2(t-r)}{r(t-2)}$ (for $2<r<t<\infty$) and $C_{\operatorname*{reg},A_{\varepsilon}}$ is a constant in the inequality

[TABLE]

In Section 4, we show that $C_{\operatorname*{reg},A_{\varepsilon}}$ depends only on the constant $C_{P}:=C_{\operatorname*{reg},I}$ (associated with the Laplace operator) and on the amplitude of the jumps in the coefficient $A_{\varepsilon}$ (cf. Theorem 4.4).

Theorem 3.1**.**

Let $A_{\scriptscriptstyle 0}\in L^{\infty}(\Omega,\mathbb{R}_{\operatorname*{sym}}^{d\times d})$ satisfy (2.2), $f\in L^{P}\left(\Omega\right)$ for some $P\in(2,+\infty)$, and $p\in(2,p^{\ast})$, where the function $p^{\ast}=p^{\ast}\left(\frac{\underline{\alpha}}{\overline{\beta}},P\right)$ is defined by (4.4). Then

[TABLE]

where

[TABLE]

$D_{\varepsilon}:=A_{\varepsilon}^{-1/2}A_{\scriptscriptstyle 0}A_{\varepsilon}^{-1/2}$, $B_{\varepsilon}:=(A_{\varepsilon}-A_{\scriptscriptstyle 0})A_{\scriptscriptstyle 0}^{-1}(A_{\varepsilon}-A_{\scriptscriptstyle 0})$, and $p^{\prime\prime}$ is a number defined in Sect. 2.

Proof.

By the triangle inequality

[TABLE]

First we estimate the discretization error. It holds

[TABLE]

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

[TABLE]

We choose $v=u-u_{\varepsilon}$ and obtain

[TABLE]

Applying the Cauchy-Schwarz inequality yields

[TABLE]

Hence,

[TABLE]

We set $p^{\prime\prime}=p/(p-2)$, apply the triangle and generalized Hölder inequalities, and obtain

[TABLE]

Now, we choose $t\in(p,p^{\ast}(\frac{\underline{\alpha}}{\overline{\beta}},P))$. By Lemma A.1 and (3), we find that

[TABLE]

Now, we apply Theorem 4.4 (where the regularity constant $C_{\operatorname*{reg},A_{\varepsilon}}$ is defined) and arrive at the estimate

[TABLE]

Notice that for $f\in L^{t}(\Omega)$ we have $\|F\|_{-1,t,\Omega}\leq\|f\|_{t,\Omega}$. The combination of (3) and (3.6) yields the desired estimate. ∎

Theorem 3.1 requires several comments. First, from Theorem 3.1 it follows that

[TABLE]

Here the quantity $\Upsilon(u_{\varepsilon,h},f,\theta(p,t))$ is fully computable provided that $C_{\operatorname*{reg},A_{\varepsilon}}$ or a certain upper bound of it is known and norms of $B_{\varepsilon}$ and $D_{\varepsilon}$ 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_{\scriptscriptstyle 0}=A_{\varepsilon}$, then $D_{\varepsilon}=I$ and $B_{\varepsilon}=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_{\varepsilon,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) and Theorem 4.4 we have

[TABLE]

Hence, we obtain another estimate

[TABLE]

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 $\varepsilon$. 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 $\delta$. 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 $\frac{1}{2}\delta$, then we can concentrate on approximations of $u_{\varepsilon}$ instead of $u$. It is natural to await that in many cases $u_{\varepsilon}$ will have better regularity properties than $u$. Then, $u_{\varepsilon,h}$ will converge to $u_{\varepsilon}$ (as $h\rightarrow 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 $\mathcal{M}_{\Omega}$.

Above observations motivate the idea to replace the $L^{\infty}$-norm of $B_{\varepsilon}$ by some $L^{p}$-norm. We are interested to make the norm $\left|\left|\left|B_{\varepsilon}\right|\right|\right|_{p^{\prime\prime},\Omega}$ small and proportional to a certain positive power of $\varepsilon$. For example, let $A_{\scriptscriptstyle 0}=\kappa_{\scriptscriptstyle 0}I$ and $A_{\varepsilon}=\kappa_{\varepsilon}I$, where $\kappa_{\scriptscriptstyle 0}$ is a jump function and $\kappa_{\varepsilon}$ is a continuous piecewise affine function approximating this jump in the $\varepsilon$–strip. Then elementary computations show that $\left|\left|\left|B_{\varepsilon}\right|\right|\right|_{p^{\prime\prime},\Omega}\sim\varepsilon^{1/p^{\prime\prime}}$. Similar a priori analysis is of course possible for more complicated matrices $A_{\scriptscriptstyle 0}$ and $A_{\varepsilon}$.

4 $\boldsymbol{W^{1,p}}$-Regularity Results for Second Order Elliptic

Problems with Rough Coefficients

Now our goal is to derive $L^{p}\left(\Omega\right)$-regularity estimates for the gradient of $\nabla u$ for some $p>2$. We start from the Poisson problem (in this case, $A_{\scriptscriptstyle 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_{\scriptscriptstyle 0}\in L^{\infty}(\Omega,\mathbb{R}_{\operatorname*{sym}}^{d\times d})$. It is important, that our estimates depend only on the amplitude of jumps in the coefficients of $A_{\scriptscriptstyle 0}$.

Consider the following problem: For a given $F\in W^{-1,p}\left(\Omega\right)$ find $\psi\in W_{0}^{1,p}\left(\Omega\right)$ such that

[TABLE]

Theorem 4.1** ([20]).**

Let $1<p<\infty$. Then, for every $F\in W^{-1,p}\left(\Omega\right)$, the problem (4.1) has a unique solution $\psi\in W_{0}^{1,p}\left(\Omega\right)$ which meets the estimate

[TABLE]

with the Laplace $W^{1,p}$-regularity constant $C_{p}$ and

[TABLE]

Remark 4.2**.**

The constant $C_{p}$ is independent of $F$ (and $\psi$) but depends on $\Omega$, $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:=\left\{\left(p,P\right):2<P<\infty,\ 2\leq p\leq P\right\}$ and introduce the function

[TABLE]

as well as the function $p^{\ast}:\left[0,1\right]\times\left(2,\infty\right)\rightarrow\mathbb{R}$

[TABLE]

Remark 4.3**.**

For the function $p^{\ast}\left(\cdot,P\right)$ it holds

[TABLE]

(cf. Figure 1). The function $p^{\ast}$ has the explicit representation

[TABLE]

Writing $A_{\scriptscriptstyle 0}$ as a perturbation of the identity and deducing the $L^{p}$-bound of $\nabla u$ from the $L^{p}$-bound for the solution of the Poisson problem (4.1) we obtain the following result for a uniformly elliptic $L^{\infty}$-coefficient $A_{\scriptscriptstyle 0}$.

Theorem 4.4**.**

Let $A_{\scriptscriptstyle 0}\in L^{\infty}(\Omega,\mathbb{R}_{\operatorname*{sym}}^{d\times d})$ satisfy (2.2), $F\in W^{-1,P}\left(\Omega\right)$ for $2<P<\infty$ and $2\leq p<p^{\ast}\left(\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}},P\right)$. Then the solution of (2.3) exists in $W_{0}^{1,p}\left(\Omega\right)$ and meets the estimate

[TABLE]

with

[TABLE]

For a proof we refer to [16, 9, 24].

Remark 4.5**.**

For the Poisson problem, we have $\alpha_{\scriptscriptstyle 0}=\beta_{\scriptscriptstyle 0}=1$ so that $p^{\ast}\left(\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}},P\right)=P$. The constant in (4.5) simplifies for $p=P$ to $C_{\operatorname*{reg},A_{\scriptscriptstyle 0}}=C_{P}$ and (4.2) is reproduced. 2. 2.

Let $P\geq 2$. For a bounded domain $\Omega$ with $\partial\Omega\in C^{1}$ there exists some $1\leq C_{\operatorname*{L}}=\mathcal{O}\left(1\right)$ such that

[TABLE]

(cf. Remark 5.6).

Theorem 4.6**.**

Let $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$, with $\partial\Omega\in C^{1}$. Let $\alpha_{\scriptscriptstyle 0},\beta_{\scriptscriptstyle 0}$ be as in (2.2) and assume that (4.6) holds.

(small perturbations) For any $P\in\left]2,\infty\right[$, consider $A_{\scriptscriptstyle 0}\in L^{\infty}\left(\Omega,\mathbb{R}_{\operatorname{sym}}^{d\times d}\right)$ in (2.3) with spectral bounds $\alpha_{\scriptscriptstyle 0},\beta_{\scriptscriptstyle 0}$ such that $\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}\geq 1-\frac{1}{2C_{\operatorname*{L}}P^{d+1}}$, and let $F\in W^{-1,P}\left(\Omega\right)$. Then, for any $2\leq p<P$ it holds

[TABLE] 2. 2.

(large perturbations) For any $P\in\left]2,\infty\right[$, consider $A_{\scriptscriptstyle 0}\in L^{\infty}\left(\Omega,\mathbb{R}_{\operatorname{sym}}^{d\times d}\right)$ in (2.3) with spectral bounds $\alpha_{\scriptscriptstyle 0},\beta_{\scriptscriptstyle 0}$ such that $\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}\leq 1-C_{P}^{-1}$. For $c\in\left(0,1\right)$ define111By a Taylor argument it follows that the right-hand side in (4.7) is bounded from below by the expression

$p_{\max}\geq 2+4c\frac{\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}}{\log C_{P}}\left(\frac{1}{2}-\frac{1}{P}\right)$

which shows the qualitative dependence on $c$, $\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}$, $C_{P}$ and $P$ better.

[TABLE]

Then, for any $2\leq p<p_{\max}$ it holds

[TABLE]

Proof.

@1 The condition $\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}\geq 1-\frac{1}{2C_{\operatorname*{L}}P^{d+1}}$ implies $\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}\geq 1-\frac{1}{2C_{P}}>1-\frac{1}{C_{P}}$ so that (cf. (4.3))

[TABLE]

Hence, we may choose any $2\leq p<P$ for the following. It is easy to see that then $\eta\left(p,P\right)\leq 1$ and (4.5) gives us

[TABLE]

This leads to the first assertion.

@2 The condition $0\leq\frac{\alpha_{\scriptscriptstyle 0}}{\beta_{\scriptscriptstyle 0}}\leq 1-C_{P}^{-1}$ implies that for any $c\in(0,1)$ the number $p_{\max}$ satisfies

[TABLE]

and Theorem 4.4 implies $\nabla u\in L^{p}\left(\Omega\right)$ for $2\leq p<p_{\max}$. Note that then $\eta\left(p,P\right)\leq\eta\left(p_{\max},P\right)=c\frac{\log\left(\frac{1}{1-\frac{\alpha_{0}}{\beta_{\scriptscriptstyle 0}}}\right)}{\log C_{P}}$ so that

[TABLE]

For the regularity constant we obtain

[TABLE]

and the assertion follows.∎

5 Analysis of the Laplace $\boldsymbol{W^{1,p}}$-Regularity Constant

for the Full Space Problem

Let $f$ be a measurable function on a domain $\Omega$ (bounded or unbounded) in $\mathbb{R}^{d}$. The distribution function $\mu_{f}$ is defined by 222For a measurable subset $M\subset\mathbb{R}^{d}$ we set $|M|:=\int_{M}1$.

[TABLE]

for $t>0$ and measures the relative size of $f$. A basic property of $\mu$ is given by the following lemma which is proved in [12, Lemma 9.7].

Lemma 5.1**.**

For any $p>0$ and $|f|^{p}\in L^{1}(\Omega)$, we have

[TABLE]

For $f\in L^{p}(\Omega)$, $1<p<\infty$, the Newtonian potential of $f$ is defined as

[TABLE]

where $G$ is the fundamental solution of Laplace’s equation which is given by

[TABLE]

where $\Gamma\left(\cdot\right)$ denotes the Gamma function. For fixed $i,j$ we define the linear operator $T:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ by

[TABLE]

Further for $d\geq 2$ and $1<p<2$ we define the constants

[TABLE]

and

[TABLE]

Notation 5.2**.**

For $r>0$ and $x\in\mathbb{R}^{d}$, we denote the ball with radius $r$ around $x$ by $B_{r}\left(x\right):=\left\{y\in\mathbb{R}^{d}:\left\|y-x\right\|<r\right\}$. The dimension $d$ is clear from the argument $x$. We write short $B_{r}$ for $B_{r}\left(0\right)$.

Lemma 5.3**.**

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain, $f\in L^{2}(\Omega)$ and $t>0$. Then we have the estimate

[TABLE]

with $C(d)$ as in (5.1).

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 $\Omega$ and fix a cube $K_{0}\supset\Omega$ so that for fixed $t>0$ we have

[TABLE]

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

[TABLE]

are similarly subdivided and the process is repeated indefinitely. In this way we obtain a sequence of parallel subcubes $(K_{\ell})_{\ell=1}^{\infty}$ such that

[TABLE]

and

[TABLE]

The function $f$ is now split into a “good part” $g$ defined by

[TABLE]

and a “bad part” $b=f-g$. Clearly,

[TABLE]

Since $T$ is linear we have $Tf=Tg+Tb$ and thus

[TABLE]

From [12, p. 232] we know that

[TABLE]

In a next step we want to estimate $\mu_{Tb}(t/2)$. Writing

[TABLE]

we have

[TABLE]

Let us now fix some $\ell$ and a sequence $b_{\ell m}\subset C_{0}^{\infty}(K_{\ell})$ converging to $b_{\ell}$ in $L^{2}(\Omega)$ and satisfying

[TABLE]

Then for $x\notin K_{\ell}$ we have the relation

[TABLE]

where $\bar{y}_{\ell}$ denotes the center of $K_{\ell}$. Further for $x\notin K_{\ell}$ and $y\in K_{\ell}$ it holds

[TABLE]

for some point $\zeta_{y}$ between $y$ and $\bar{y}_{\ell}$. Some computations show that

[TABLE]

where $\delta_{ij}$ denotes the Kronecker symbol. This leads to the estimate

[TABLE]

Note that

[TABLE]

so that

[TABLE]

where $\delta_{\ell}:=\operatorname*{diam}K_{\ell}$. The combination of (5.6) and (5.7) yields

[TABLE]

It is easy to see that for $x\notin B_{\delta_{\ell}}\left(\bar{y}_{\ell}\right)$ it holds

[TABLE]

Thus,

[TABLE]

We set $F^{\ast}:=\bigcup_{\ell}B_{\delta_{\ell}}\left(\bar{y}_{\ell}\right)$ and $G^{\ast}:=K_{0}\backslash F^{\ast}$. Letting $m\rightarrow\infty$ and summing over $\ell$ we get

[TABLE]

By Lemma 5.1 and using the last estimate we obtain

[TABLE]

Moreover by [12, p. 234] and (5.3) we have with $F:=\bigcup_{\ell}K_{\ell}$

[TABLE]

Finally we get by (5.4), (5.5), (5.8) and (5.9)

[TABLE]

∎

Theorem 5.4** (Calderon-Zygmund estimate).**

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain and $f\in L^{p}(\Omega)$, $1<p<\infty$. Then it holds

[TABLE]

with

[TABLE]

and $C(d,p)$ as in (5.2).

Proof.

For the case $p=2$ we refer to [12, Theorem 9.9].

By Lemma 5.1 and since $\|Tf\|_{2,\Omega}=\|f\|_{2,\Omega}$ we know that

[TABLE]

for all $t>0$ and all $f\in L^{2}(\Omega)$.

Further by Lemma 5.3 we have

[TABLE]

for all $t>0$ and all $f\in L^{2}(\Omega)$. Thus it follows by the Marcinkiewicz interpolation theorem (cf. Theorem A.2 with $q=1$ and $r=2$) that

[TABLE]

with $C(d,p)$ as in (5.2). By a duality argument (cf. [12, Theorem 9.9]) we obtain that

[TABLE]

Next we employ the Riesz-Thorin interpolation theorem to remove the singular behavior of $C\left(d,p\right)$ as $p\rightarrow 2$. Note that $T:L^{p}\left(\Omega\right)\rightarrow L^{p}\left(\Omega\right)$ is continuous for $p\in\left\{p_{0},p_{1}\right\}$ with $p_{0}=3/2$ and $p_{1}=2$ and $T:L^{p_{0}}\left(\Omega\right)+L^{p_{1}}\left(\Omega\right)\rightarrow L^{p_{0}}\left(\Omega\right)+L^{p_{1}}\left(\Omega\right)$ (observe that $L^{p_{0}}\left(\Omega\right)+L^{p_{1}}\left(\Omega\right)=L^{p_{0}}\left(\Omega\right)$ since $\Omega$ is bounded). We know that

[TABLE]

Let

[TABLE]

Then the Riesz-Thorin interpolation theorem implies that $T:L^{p}\left(\Omega\right)\rightarrow L^{p}\left(\Omega\right)$ is bounded and

[TABLE]

for $\frac{3}{2}<p<2$. Note that $p_{0}^{\prime}=3$ and $p_{1}^{\prime}=2$. Let

[TABLE]

Applying again the Riesz-Thorin interpolation theorem yields

[TABLE]

for $2<p<3$. In total we have proved that

[TABLE]

with $C_{1}(d,p)$ as in (5.10). ∎

Remark 5.5**.**

Note that $C_{1}(d,p)=C_{1}(d,p^{\prime})$ for $1<p<\infty$. Further observe that $C_{1}(d,p)\geq 1$. 2. 2.

Since for $p\geq 3$ we have

[TABLE]

the constant $C_{1}(d,p)$ is of order $\mathcal{O}(p)$ as $p\rightarrow\infty$.

Remark 5.6**.**

*Observe that $T$ can be defined as a bounded operator on $L^{p}(\Omega)$ even if $\Omega$ is unbounded. In this case Theorem 5.4 still holds provided that $d\geq 3$ (cf. [12]). The constant $C_{1}(d,p)$ in the Calderon-Zygmund estimate is equal to the Laplace $W^{1,p}$-regularity constant for the full space problem (cf. [20, Chapter II, Lem. 2.1]).

The analysis of the Laplace $W^{1,p}$-regularity constant for bounded domains with $C^{1}$-boundary can be found in [20, Chapter II]. However, the usual localization techniques such as using cutoff functions on overlapping balls does not lead to sharp estimates.

From the analysis in [20] it follows that for a bounded domain $\Omega$ with $\partial\Omega\in C^{1}$ and $2\leq p<\infty$ there exists some $1\leq C_{\operatorname*{L}}=\mathcal{O}\left(1\right)$ such that*

[TABLE]

with $C_{p}$ as in (4.2).

Appendix A Interpolation Estimates

Lemma A.1**.**

Let $2<r<t<\infty$ and $\theta\in(0,1)$ be such that

[TABLE]

Then if $u\in L^{t}(\Omega)$, we have the estimate

[TABLE]

A proof can be found in [1]. Note that inequality (A.1) also holds for vector-valued functions.

Theorem A.2** (Marcinkiewicz interpolation theorem).**

Let $S$ be a linear mapping from $L^{q}(\Omega)\cap L^{r}(\Omega)$ into itself, $1\leq q<r<\infty$ and suppose that there are constants $S_{1}$ and $S_{2}$ such that

[TABLE]

for all $f\in L^{q}(\Omega)\cap L^{r}(\Omega)$ and $t>0$. Then $S$ extends as a bounded linear mapping from $L^{p}(\Omega)$ into itself for any $p$ such that $q<p<r$ and

[TABLE]

for all $f\in L^{q}(\Omega)\cap L^{r}(\Omega)$, where

[TABLE]

For a proof we refer to [12, Theorem 9.8].