Numerical homogenization for nonlinear strongly monotone problems

TL;DR

This paper introduces a multiscale finite element method for nonlinear monotone equations that achieves optimal error estimates without requiring structural assumptions, demonstrated through numerical examples including Richards equation.

Contribution

A novel multiscale method for nonlinear monotone problems that constructs problem-adapted spaces via local linear problems, avoiding structural assumptions.

Findings

Method achieves optimal a priori error estimates

Efficient due to linearity and localization of local problems

Validated by numerical examples including Richards equation

Abstract

In this work we introduce and analyze a new multiscale method for strongly nonlinear monotone equations in the spirit of the Localized Orthogonal Decomposition. A problem-adapted multiscale space is constructed by solving linear local fine-scale problems which is then used in a generalized finite element method. The linearity of the fine-scale problems allows their localization and, moreover, makes the method very efficient to use. The new method gives optimal a priori error estimates up to linearization errors. The results neither require structural assumptions on the coefficient such as periodicity or scale separation nor higher regularity of the solution. The effect of different linearization strategies is discussed in theory and practice. Several numerical examples including stationary Richards equation confirm the theory and underline the applicability of the method.

Equations

$$-\nabla\cdot\bigl{(}A(x,\nabla u)\bigr{)}=f,$$

$$\displaystyle-\nabla\cdot\bigl{(}A(x,\nabla u)\bigr{)}$$

$$u$$

$$\mathcal{B}(u;v):=\bigl{(}A(x,\nabla u),\nabla v\bigr{)}_{\Omega}=(f,v)_{\Omega}\qquad\text{for all}\quad v\in H^{1}_{0}(D).$$

$$|\mathcal{B}(v_1;\psi )-\mathcal{B}(v_2;\psi )|\le \Lambda |v_1-v_2{|}_1|\psi {|}_1\text{and}\mathcal{B}(v;v-\psi )-\mathcal{B}(\psi ;v-\psi )\ge \lambda |v-\psi {|}_1^2$$

$$|u{|}_1\le \Lambda /\lambda \Vert f{\Vert }_0,$$

$$\mathcal{B}(u_H;v_H)=(f,v_H{)}_{\Omega }\text{for all}{v}_H\in V_H.$$

$$|u-u_H{|}_1\lesssim \underset{v_H\in V_H}{\inf }|u-v_H{|}_1,$$

$$\Vert u-u_H{\Vert }_1\le CH\Vert u{\Vert }_{H^2(\Omega )}\text{and}\Vert u-u_H{\Vert }_0\le C{H}^2\Vert u{\Vert }_{H^2(\Omega )},$$

$$|I_{H}v{|}_{1,T}$$

$$\Vert v-I_{H}v{\Vert }_{0,T}$$

$$(E_{H}(v))(z)=\sum_{T\in\mathcal{T}_{H},\;z\in T}v|_{T}(z)\bigg{/}\operatorname{card}\{K\in\mathcal{T}_{H}\,,\,z\in K\}.$$

$$\mathcal{B}(u_H^{\mathrm{nl},\mathrm{ms}};v)=(f,v{)}_{\Omega }\text{for all}v\in V_H^{\mathrm{nl},\mathrm{ms}},$$

$$\mathcal{B}(v_H^{\mathrm{nl},\mathrm{ms}};w)=0\text{for all}{v}_H^{\mathrm{nl},\mathrm{ms}}\in V_H^{\mathrm{nl},\mathrm{ms}}\text{and all}w\in W.$$

$$|u-u_H^{\mathrm{nl},\mathrm{ms}}{|}_1\lesssim H\Vert f{\Vert }_0$$

$$\mathcal{B}(v_H-{\mathcal{Q}}^{\mathrm{nl}}{v}_H;w)=0\text{for all}w\in W.$$

$${\lambda }_L|\xi {|}^2\le {\mathbf{A}}_L(x,\nabla u^{\ast })\xi \cdot \xi \le {\Lambda }_L|\xi {|}^2\text{for all}x\in \Omega ,u^{\ast }\in H_0^1(\Omega ),\xi \in {\mathbb{R}}^d.$$

$$\mathcal{A}(v_H-\mathcal{Q}{v}_H,w)=0\text{for all}w\in W,$$

$$\mathcal{A}(v,\psi ):=(\mathfrak{A}(x)\nabla v,\nabla \psi {)}_{\Omega }.$$

$${\int }_{\Omega }{A}_L(x,\nabla u^{\ast },\nabla \mathcal{Q}{v}_H)\cdot \nabla wdx={\int }_{\Omega }{A}_L(x,\nabla u^{\ast },\nabla v_H)\cdot \nabla wdx\text{for all}w\in W.$$

$$\mathcal{B}(u_H^{\mathrm{ms}};v_H^{\mathrm{ms}})=(f,v_H^{\mathrm{ms}}{)}_{\Omega }\text{for all}{v}_H^{\mathrm{ms}}\in V_H^{\mathrm{ms}}.$$

$$\text{N}(T)=\bigcup _{K\in {\mathcal{T}}_H,T\cap K\ne \varnothing }K$$

$$\underset{T\in {\mathcal{T}}_H}{\max }\mathrm{card}\{K\in {\mathcal{T}}_H:K\subset {\text{N}}^m(T)\}\le C_{\mathrm{ol},m}.$$

$${\mathcal{A}}_{{\text{N}}^m(T)}({\mathcal{Q}}_{T,m}{v}_H,w)={\mathcal{A}}_T(v_H,w)\text{for all}w\in W({\text{N}}^m(T)).$$

$${\int }_{{\text{N}}^m(T)}\mathfrak{A}\nabla q_{T,m}^{(j)}\cdot \nabla wdx={\int }_T\mathfrak{A}{e}_j\cdot \nabla wdx\text{for all}w\in W({\text{N}}^m(T)),$$

$$\Bigl{\{}\lambda_{z}-\sum_{T\in\mathcal{T}_{H},\;z\in T}\sum_{j=1}^{d}\Bigl{(}\frac{\partial}{\partial x_{j}}\lambda_{z}|_{T}\Bigr{)}q_{T,m}^{(j)}\Bigr{\}}_{z}.$$

$$\mathcal{B}(u_{H,m};v_{H,m})=(f,v_{H,m}{)}_{\Omega }\text{for all}{v}_{H,m}\in V_{H,m}.$$

$$\mathcal{B}(u_{H,m}^{PG};v_{H,m})=(f,v_{H,m}{)}_{\Omega }\text{for all}{v}_{H,m}\in V_{H,m}.$$

$$|(\mathcal{Q}-{\mathcal{Q}}_m)v_H{|}_1\lesssim C_{\mathrm{ol},m}^{1/2}{\beta }^m|v_H{|}_1.$$

$$|u-u_{H,m}{|}_1\lesssim (H+C_{\mathrm{ol},m}^{1/2}{\beta }^m)\Vert f{\Vert }_0+{\eta }_{\mathrm{lin}}((\mathrm{id}-\mathcal{Q})I_{H}u)$$

$$|u-u_{H,m}{|}_1\lesssim (H+C_{\mathrm{ol},m}^{1/2}{\beta }^m)\Vert f{\Vert }_0+{\eta }_{\mathrm{lin}}(u)$$

Full text

Abstract. In this work we introduce and analyze a new multiscale method for strongly nonlinear monotone equations in the spirit of the Localized Orthogonal Decomposition. A problem-adapted multiscale space is constructed by solving linear local fine-scale problems which is then used in a generalized finite element method. The linearity of the fine-scale problems allows their localization and, moreover, makes the method very efficient to use. The new method gives optimal a priori error estimates up to linearization errors. The results neither require structural assumptions on the coefficient such as periodicity or scale separation nor higher regularity of the solution. The effect of different linearization strategies is discussed in theory and practice. Several numerical examples including stationary Richards equation confirm the theory and underline the applicability of the method.

AMS subject classifications. 65N15, 65N30, 35J60, 74Q15

1 Introduction

Linear constitutive laws like Hooke’s law in mechanics, Ohm’s law in electromagnetics, or Darcy’s law in fluid flow are very popular, but they are often not accurate enough in practical applications, for instance for high intensities. Instead, nonlinear effects in the constitutive laws have to be taken into account which are often experimentally found and determined, see [41] for a general overview. In this article, we consider as model problem the following nonlinear monotone elliptic equation

[TABLE]

where the exact assumptions as well as boundary conditions are specified later. It is a representative model problem for quasilinear partial differential equations (PDEs) as they occur in mean curvature flow or for non-Newtonian fluids. The transition from linear to nonlinear problems comes with huge additional challenges for the numerical treatment and analysis. As an illustrating example we mention optimal order $L^{2}$-estimates for the finite element method: The classical Aubin-Nitsche trick for linear problems is not applicable, so that, for a long time, only optimal order estimates in the energy norm [10] were known, see [5] and the discussion therein. A similar observation applies to the effect of numerical integration, see [16].

With the view on practical applications such as fluid flow or elasticity, we do not only have to consider nonlinear constitutive laws as discussed above, but also have to consider (spatial) multiscale features in the material coefficients (here, in $A$). For instance, a fluid such as groundwater flows over large distances, while the properties of the soil changes over small distances, see, e.g., [38]. Hence, for applications such as the (quasilinear) porous medium equation, $A$ is subject to rapid variations and/or discontinuities on fine spatial scales or even a cascade of (non-separable) scales. This coincidence of multiscale features and nonlinear material laws makes the problem intractable for standard methods. For example, the finite element method [2, 10, 16] will only give optimal convergence in the asymptotic regime, i.e., if the mesh resolves all features and scales present, which is prohibitively expensive even with today’s computational resources.

In the case of spatially periodic $A$ (with period $\varepsilon\ll 1$), homogenization results using two-scale convergence [7, 30] prove that the solutions of the above model problem converge to the solution of an again monotone elliptic (homogenized) problem for $\varepsilon\to 0$ . The nonlinear effective diffusion tensor can be computed by solving nonlinear so-called cell problems. The (finite element) heterogeneous multiscale method is inspired by this analytical process and it is studied successfully for nonlinear problems in a series of papers [1, 3, 4, 6, 20, 24]. In most cases, the macroscopic nonlinear form involves nonlinear reconstruction operators which require the solution of nonlinear cell problems at each macroscopic quadrature point to incorporate fine-scale information. For parabolic equations, [4] linearizes the macroscopic and cell computations using information from the previous time step. The sparse multiscale FEM [26] tries to reduce the complexity of solving cell problems and a homogenized equation by the introduction of sparse approximations. Another idea is to modify or enrich the standard finite element basis by problem-adapted functions. This is used in the (generalized) multiscale finite element method, for which nonlinear problems are discussed in [9, 11, 12]. Again nonlinear problems have to be solved locally to construct the problem-adapted functions.

The main contribution of this article is the introduction of a new multiscale method for nonlinear strongly monotone problems and its numerical analysis. The idea is to construct a multiscale space by solving local fine-scale problems in the spirit of the Localized Orthogonal Decomposition (LOD) [32, 36]. In contrast to the above discussed methods, the basis construction only requires the solution of linear problems and hence is embarrassingly easy. Moreover, this linearization idea drastically reduces the computational effort for generating a problem-dependent basis and thereby provides a conceptually new view on the treatment of nonlinear multiscale problems. We derive optimal convergence rates (with respect to the mesh size $H$) up to linearization errors without any assumption on the regularity of the exact solution or special properties such as periodicity or scale separation for the coefficient. The occurring linearization errors and resulting possible choices of the linearization are discussed and compared. Extensive numerical experiments show the good performance of the method in agreement with the theoretical estimates. We study periodic as well as completely random multiscale coefficients and also include a model for stationary Richards equation with a high contrast channel. Besides several linear problem classes, the LOD has already been studied for semilinear equations [22] and a nonlinear eigenvalue problem related to the Gross-Pitaevskii equation [23]. These problems, however, are only semilinear and can therefore be handled easier. Yet, we emphasize that these previous works can be re-interpreted in the current framework. We mention the close connections of the LOD to (analytical) homogenization [17, 37], domain decomposition iterative solvers [28, 29, 37], and so-called gamblets [34, 35]. Hence, the current approach can give interesting and useful insights in these areas for nonlinear problems in the future as well.

The article is organized as follows: Section 2 introduces the setting and the standard finite element discretization. We introduce the multiscale method including linearization and localization in Section 3. The arising errors are analyzed in Section 4. Finally, we present several numerical experiments confirming our theory and showing possible applications in Section 5.

2 Problem formulation and discretization

In this section we formulate the considered model problem and introduce necessary finite element prerequisites. We use standard notation on Sobolev spaces. Throughout the whole article, let $\Omega\subset\mathbb{R}^{d}$ be a bounded Lipschitz domain. For a subdomain $D\subset\Omega$, let $\|\cdot\|_{0,D}$, $\|\cdot\|_{1,D}$, and $|\cdot|_{1,D}$ denote the standard $L^{2}(D)$-norm, $H^{1}(D)$-norm, and $H^{1}(D)$-semi norm, respectively. Furthermore, $(\cdot,\cdot)_{D}$ denotes the standard $L^{2}$ scalar product on $D$. We will omit the subscript $D$ if it equals the full domain $\Omega$.

2.1 Model problem

We consider the following nonlinear elliptic problem: Find $u:\Omega\to\mathbb{R}$ such that

[TABLE]

with a right-hand side $f\in L^{2}(\Omega)$. The corresponding weak formulation, with which we will work in the following, reads: Find $u\in H^{1}_{0}(\Omega)$ such that

[TABLE]

For simplicity, we restrict ourselves to homogeneous Dirichlet boundary conditions, but non-homogeneous and Neumann boundary conditions could be treated as well, see [21]. Moreover, we focus on nonlinearities in the highest derivative only, additional (nonlinear) low-order terms can easily be handled as well, cf. [22]. We now specify our assumptions on $A$.

Assumption 2.1**.**

The nonlinearity $A:\Omega\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ satisfies

$A(\cdot,\xi)\in L^{\infty}(\Omega;\mathbb{R}^{d})$* for all $\xi\in\mathbb{R}^{d}$ and $A(x,\cdot)\in C^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$ for almost every $x\in\Omega$;* 2. 2.

there is $\Lambda>0$ such that $|A(x,\xi_{1})-A(x,\xi_{2})|\leq\Lambda|\xi_{1}-\xi_{2}|$ for almost every $x\in\Omega$ and all $\xi_{1},\xi_{2}\in\mathbb{R}^{d}$; 3. 3.

there is $\lambda>0$ such that \bigl{(}A(x,\xi_{1})-A(x,\xi_{2})\bigr{)}\cdot(\xi_{1}-\xi_{2})\geq\lambda|\xi_{1}-\xi_{2}|^{2} for almost every $x\in\Omega$ and all $\xi_{1},\xi_{2}\in\mathbb{R}^{d}$; 4. 4.

$|A(x,0)|\leq C_{0}$* for almost every $x\in\Omega$.*

Assumption 2.1 implies that

[TABLE]

for all $v,v_{1},v_{2},\psi\in H^{1}_{0}(\Omega)$. Therefore, the model problem (2.2) has a unique solution $u\in H^{1}_{0}(\Omega)$, which satisfies

[TABLE]

see [10, Chapter 5]. As discussed in the introduction, we implicitly assume that $A$ is subject to rapid oscillations or discontinuities on a rather fine scale with respect to the spatial variable $x$.

We write $a\lesssim b$ in short for $a\leq Cb$ with a constant C independent of the mesh size $H$ and the oversampling parameter $m$ introduced later. However, $C$ may depend on the monotonicity and Lipschitz constants $\lambda,\Lambda$ of $A$ (cf. Assumption 2.1).

2.2 Finite element discretizations

We cover $\Omega$ with a regular mesh $\mathcal{T}_{H}$ consisting of simplices; however, a mesh with quadrilaterals would equally be possible. The mesh is assumed to be shape regular in the sense that the aspect ratio of the elements of $\mathcal{T}_{H}$ is bounded uniformly from below. We introduce the mesh size $H=\max_{T\in\mathcal{T}_{H}}\operatorname{diam}T$ and assume that this is rather coarse, in particular, $\mathcal{T}_{H}$ does not resolve the possible heterogeneities in $A$. We discretize the space $H^{1}_{0}(\Omega)$ with the lowest order Lagrange elements over $\mathcal{T}_{H}$, and denote this space by $V_{H}$. This means that $V_{H}=H^{1}_{0}(\Omega)\cap\mathcal{S}^{1}(\mathcal{T}_{H})$, where $\mathcal{S}^{1}(\mathcal{T}_{H})$ denotes the space of element-wise polynomials of total degree $\leq 1$.

The standard finite element method now seeks a (discrete) solution $u_{H}\in V_{H}$ such that

[TABLE]

This results in a nonlinear system which can be (approximatively) solved via an iteration such as Newton’s method. It is well-known that the properties of $A$ and Galerkin orthogonality imply

[TABLE]

see [10, Chapter 5]. This quasi-optimality by the way holds for any conforming subset $\tilde{V}_{H}\subset H^{1}_{0}(\Omega)$. For the standard finite element method (with linear elements) it is furthermore well-known to have the following error estimates

[TABLE]

see [2]. The higher regularity ($u\in H^{2}(\Omega)$) of the exact solution required in those estimates may not be attainable for nonlinearities $A$ with spatial discontinuities. Even if $u\in H^{2}(\Omega)$, the corresponding norm depends on spatial derivatives of $A$ which behave like $\varepsilon^{-1}$ for coefficients varying on a scale $\varepsilon$. In practice this implies that $H$ needs to be at least $\varepsilon$ in order to observe the linear convergence in the $H^{1}(\Omega)$-norm. In other words, for small $\varepsilon$, there is a large pre-asymptotic region where the error stagnates (at a high level) in practice.

The goal of the multiscale method presented in Section 3 is to circumvent both issues (higher regularity of the solution and dependence on the variations of $A$). At the heart of the method is the choice of a suitable interpolation operator and we now introduce the required properties as well as an appropriate example. Let $I_{H}:H^{1}_{0}(\Omega)\to V_{H}$ denote a bounded local linear projection operator, i.e., $I_{H}\circ I_{H}=I_{H}$, with the following stability and approximation properties for all $v\in H^{1}_{0}(\Omega)$

[TABLE]

where the constants are independent of $H$ and $\textup{N}(T):=\{K\in\mathcal{T}_{H}:K\cap T\neq\emptyset\}$ denotes the neighborhood of an element. A possible choice (which we use in our implementation of the method) is to define $I_{H}:=E_{H}\circ\Pi_{H}$. $\Pi_{H}$ is the $L^{2}$-projection onto the elementwise affine functions $\mathcal{S}^{1}(\mathcal{T}_{H})$, and $E_{H}$ is the averaging operator that maps discontinuous functions in $\mathcal{S}^{1}(\mathcal{T}_{H})$ to $V_{H}$ by assigning to each free vertex the arithmetic mean of the corresponding function values of the neighboring cells, that is, for any $v\in\mathcal{S}^{1}(\mathcal{T}_{H})$ and any vertex $z$ of $\mathcal{T}_{H}$,

[TABLE]

For further details on suitable interpolation operators we refer to [15].

3 Computational multiscale method

In the following, we assume that an interpolation operator $I_{H}:H^{1}_{0}(\Omega)\to V_{H}$ satisfying the projection property as well as (2.5) and (2.6) is at hand. Abbreviating $W:=\ker I_{H}$, we have the splitting $H^{1}_{0}(\Omega)=V_{H}\oplus W$. The main idea of the Localized Orthogonal Decomposition [32, 36] is to make this splitting problem-dependent. In the linear elliptic case the splitting is orthogonalized with respect to the energy scalar product. Below, we discuss how this idea can be transferred to the nonlinear case. We introduce a linearization procedure in the next subsection which makes the computation of a multiscale space in the spirit of the LOD possible. Afterwards, we present the localized computation of the new multiscale basis functions.

3.1 An ideal method and its linearization

Motivated by linear elliptic equations, one could (naively) try to introduce a Galerkin method over a subset $V_{H}^{\mathrm{nl},\mathrm{ms}}\subset H^{1}_{0}(\Omega)$, i.e., we seek $u_{H}^{\mathrm{nl},\mathrm{ms}}$ such that

[TABLE]

where the set $V_{H}^{\mathrm{nl},\mathrm{ms}}$ is defined via

[TABLE]

This is the orthogonalization idea behind the original method, see [32, 36]. Due to the quasi-optimality (2.4) and the properties (2.5) and (2.6) of $I_{H}$, one obtains the a priori error estimate

[TABLE]

with optimal rate in the mesh size, independent of the regularity of the continuous solution $u$. This estimate is derived similar to the linear case [32, 36]. Because of the nonlinearity of $\mathcal{B}$ in its first argument, however, $V_{H}^{\mathrm{nl},\mathrm{ms}}$ is no longer a linear subspace. To be more precise, it holds $V_{H}^{\mathrm{nl},\mathrm{ms}}=(\operatorname{id}-\mathcal{Q}^{\mathrm{nl}})V_{H}$, where $\mathcal{Q}^{\mathrm{nl}}:V_{H}\to W$ solves

[TABLE]

Here, we clearly see that $\mathcal{Q}^{\mathrm{nl}}$ is a nonlinear operator. Therefore, it is by no means clear whether the proposed multiscale method is at all well defined. Even if this is the case, the method is very complicated as it involves two coupled nonlinear problems, where (3.2) is additionally posed on the fine scale.

Here, we propose the following simple yet effective linearization approach. We approximate the nonlinearity $A(x,\nabla v)$ by a function $A_{L}(x,\nabla u^{*},\nabla v)$. Here, $A_{L}:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ is affine in its last argument and we call $u^{*}\in H^{1}_{0}(\Omega)$ the linearization point. We make the following assumption on $A_{L}$.

Assumption 3.1**.**

Write $A_{L}(x,\nabla u^{*},\nabla v)=\mathbf{A}_{L}(x,\nabla u^{*})\nabla v+b_{L}(x,\nabla u^{*})$ with $\mathbf{A}_{L}(x,\nabla u^{*})\in\mathbb{R}^{d\times d}$ and $b_{L}(x,\nabla u^{*})\in\mathbb{R}^{d}$. We assume that

• •

$\mathbf{A}_{L}(\cdot,\nabla u^{*})\in L^{\infty}(\Omega;\mathbb{R}^{d\times d})$* and $b_{L}(\cdot,\nabla u^{*})\in L^{2}(\Omega;\mathbb{R}^{d})$ for all $u^{*}\in H^{1}_{0}(\Omega)$;*

• •

$\mathbf{A}_{L}(x,\nabla u^{*})$* is symmetric for all $u^{*}\in H^{1}_{0}(\Omega)$;*

• •

there exists $0<\lambda_{L}\leq\Lambda_{L}$ such that

[TABLE]

The assumption of symmetry is only made for convenience and to avoid cluttering of notation in the following. Although this linearization model may seem rather abstract, it includes Newton-type as well as Kačanov-type linearizations as illustrated in the example, cf. [13].

Example 3.2**.**

Newton-type linearizations are based on a Taylor expansion up to the first order of the nonlinearity around the linearization point. In particular, we approximate $A(x,\nabla v)\approx A(x,\nabla u^{*})+D_{\xi}A(x,\nabla u^{*})\nabla(v-u^{*})$, where $D_{\xi}A$ denotes the Jacobian of $A$ with respect to the second argument. In the notation of Assumption 3.1, this means that $b_{L}(x,\nabla u^{*})=A(x,\nabla u^{*})-D_{\xi}A(x,\nabla u^{*})\nabla u^{*}$ and $\boldsymbol{A}_{L}(x,\nabla u^{*})=D_{\xi}A(x,\nabla u^{*})$. The assumptions of strict monotonicity and Lipschitz continuity on $A$ (cf. Assumption 2.1) imply that $\boldsymbol{A}_{L}$ indeed satisfies Assumption 3.1, see [27, Lemma 6.5.2].

In the case that $A$ takes the form $A(x,\nabla u)=\alpha(x,|\nabla u|^{2})\nabla u$, Kačanov-type linearizations are very popular, which “freeze the nonlinearity”. In the language of Assumption 3.1, one sets $A_{L}(x,\nabla u^{*},\nabla v):=\alpha(x,|\nabla u^{*}|^{2})\nabla v$, and, hence, $\boldsymbol{A}_{L}(x,\nabla u^{*})=\alpha(x,|\nabla u^{*}|^{2})$ and $b_{L}(x,\nabla u^{*})=0$.

Let us come back to the linearization of (3.2). We pick and fix a function $u^{*}\in H^{1}_{0}(\Omega)$ and set $\mathfrak{A}:=\mathbf{A}_{L}(\cdot,\nabla u^{*})\in L^{\infty}(\Omega;\mathbb{R}^{d\times d})$. We define the linear correction operator $\mathcal{Q}:V_{H}\to W$ via

[TABLE]

where the bilinear form $\mathcal{A}$ is defined as

[TABLE]

Due to Assumption 3.1, the linear corrector problem (3.3) has a unique solution. Further, note that (3.3) can equivalently be written as

[TABLE]

Having linearized $\mathcal{Q}$, we define the linear multiscale space $V_{H}^{\mathrm{ms}}:=(\operatorname{id}-\mathcal{Q})V_{H}$. This problem-adapted space is used in a Galerkin method to seek $u_{H}^{\mathrm{ms}}\in V_{H}^{\mathrm{ms}}$ as solution of the nonlinear problem

[TABLE]

Let $\{\lambda_{z}\}_{z}$ be the nodal basis of $V_{H}$ (i.e., the standard hat functions). Then $\{\lambda_{z}-\mathcal{Q}\lambda_{z}\}_{z}$ forms a basis of $V_{H}^{\operatorname{ms}}$. Note that this requires only to solve linear problems. After this basis has been (pre-)computed, the nonlinear problem (3.4) can be solved with any suitable iteration method such as Newton’s method, which is rather cheap since the dimension of the multiscale space is small (note $\operatorname{dim}(V_{H}^{\mathrm{ms}})=\operatorname{dim}(V_{H})$). We emphasize that although we linked nonlinear iterative methods and linearization strategies in Example 3.2, the iterative method chosen to solve (3.4) does not need to correspond to the linearization strategy chosen for (3.3).

3.2 Localization of the basis generation

As in the linear case, the corrector problems (3.3) are global fine-scale problems, which are as expensive to solve as the solution of a (linear) multiscale model problem on a fine-scale mesh. However since (3.3) is a standard elliptic problem (cf. Assumption 3.1), we can localize these corrector problems in the well-known way for the linear case. To this end, we define the neighborhood

[TABLE]

associated with an element $T\in\mathcal{T}_{H}$. Thereby, for any $m\in\mathbb{N}_{0}$, the $m$-layer patches are defined inductively via $\textup{N}^{m+1}(T)=\textup{N}(\textup{N}^{m}(T))$ with $\textup{N}^{0}(T):=T$. The shape regularity implies that there is a bound $C_{\operatorname{ol},m}$ (depending only on $m$) of the number of the elements in the $m$-layer patch, i.e.,

[TABLE]

Throughout this article, we assume that $\mathcal{T}_{H}$ is quasi-uniform, which implies that $C_{\operatorname{ol},m}$ grows at most polynomially with $m$.

We then define the truncated correction operator $\mathcal{Q}_{m}:V_{H}\to W$ as $\mathcal{Q}_{m}=\sum_{T\in\mathcal{T}_{H}}\mathcal{Q}_{T,m}$, where for any $v_{H}\in V_{H}$ the truncated element corrector $\mathcal{Q}_{T,m}v_{H}\in W(\textup{N}^{m}(T)):=\{w\in W:w=0\text{ in }\Omega\setminus\textup{N}^{m}(T)\}$ solves

[TABLE]

Here, $\mathcal{A}_{D}$ denotes the restriction of the bilinear form $\mathcal{A}$ to the subdomain $D\subset\Omega$. We then set up the multiscale space $V_{H,m}:=(\operatorname{id}-\mathcal{Q}_{m})V_{H}$. For each element $T\in\mathcal{T}_{H}$, we only have to solve $d$ problems of type (3.6) with $v_{H}|_{T}=x_{j}$, $j=1,\ldots,d$, or precisely, the following cell problems: Find $q_{T,m}^{(j)}\in W(\textup{N}^{m}(T))$, $j=1,\ldots,d$, such that

[TABLE]

where $e_{j}$ denotes the $j$th canonical unit vector. Denoting by $\{\lambda_{z}\}_{z}$ the standard hat functions, a basis of $V_{H,m}$ is hence given by

[TABLE]

The localized multiscale method consists of replacing $V_{H}^{\mathrm{ms}}$ by $V_{H,m}$ in (3.4). More precisely, we seek (in a Galerkin method) $u_{H,m}\in V_{H,m}$ such that

[TABLE]

Again, the nonlinear problem (3.7) is solved with an iterative method, where the multiscale basis can be pre-computed. This requires the storage of all correctors, which can be very memory consuming since $q_{T,m}^{(j)}$ includes fine-scale features. Instead the correctors could also be computed on the fly inside each Newton iteration. Note that $\mathcal{Q}_{m}$ and therefore also the solution $u_{H,m}$ depends on $\mathfrak{A}$ and thereby, they implicitly depend on (i) the chosen linearization model and (ii) the chosen linearization point $u^{*}$ in Section 3.1. The choice of $u^{*}$ and its consequences will be discussed in Section 4.2 below.

Remark 3.3**.**

To avoid communication between the correctors, one can also consider the Petrov-Galerkin method to seek $u_{H,m}^{PG}\in V_{H}$ such that

[TABLE]

In the Petrov-Galerkin method, $q_{T,j}$ and $q_{T^{\prime},j}$ for $T,T^{\prime}\in\mathcal{T}_{H}$ with $T\neq T^{\prime}$ are never needed at the same time. Hence, these correctors can immediately be discarded once the contributions of element $T$ to the linear system (in each Newton iteration) are assembled. Hence, when memory becomes the limiting factor, the Petrov-Galerkin variant is preferable, see [14, 15]. Note, however, that $u_{H,m}^{PG}$ only contains information on the coarse scale $H$. The following error analysis will be restricted to the Galerkin case for simplicity.

Remark 3.4**.**

The present method is still semi-discrete since the corrector problems (3.6) are infinite-dimensional. The discretization procedure for them is equivalent to the case of linear elliptic equations: We introduce a second (fine) simplicial mesh $\mathcal{T}_{h}$ of $\Omega$ which resolves all features of $A$. Denoting by $V_{h}=H^{1}_{0}(\Omega)\cap\mathcal{S}^{1}(\mathcal{T}_{h})$ the corresponding lowest order Lagrange finite element space, we set $W_{h}(\textup{N}^{m}(T)):=\{w_{h}\in V_{h}:w_{h}=0\text{ in }\Omega\setminus\textup{N}^{m}(T),I_{H}w_{h}=0\}$ and discretize (3.6) by solving over the space $W_{h}(\textup{N}^{m}(T))$ instead of $W(\textup{N}^{m}(T))$. In the following, we work with the semi-discrete version and emphasize that similar error estimates (with respect to a reference solution $u_{h}\in V_{h}$) can be shown in the fully discrete variant, as illustrated for the linear case, see [21, 25, 32].

4 Error analysis

In this section, let the linear model $\mathfrak{A}$ be fixed, i.e., the linearization model $A_{L}$ and the linearization point $u^{*}\in H^{1}_{0}(\Omega)$ are fixed. Occasionally, we will also use $\mathfrak{A}^{u}:=\mathbf{A}_{L}(\cdot,\nabla u)$ with the above fixed linearization model $A_{L}$ and the exact solution $u$ to (2.2). $\mathcal{Q}_{m}^{u}$ then denotes the corrector associated with $\mathfrak{A}^{u}$, i.e., $\mathcal{Q}_{m}^{u}$ is defined via (3.6) with $\mathfrak{A}$ replaced by $\mathfrak{A}^{u}$. We have the following result on the error $\mathcal{Q}-\mathcal{Q}_{m}$.

Proposition 4.1**.**

Let Assumptions 2.1 and 3.1 be fulfilled. Let $\mathcal{Q}$ be defined in (3.3) and $\mathcal{Q}_{m}$ defined via (3.6). There exists $0<\beta<1$ such that for any $v_{H}\in V_{H}$

[TABLE]

Proposition 4.1 follows from the linear elliptic case in [21, 32, 36]. The main idea is that the ideal element corrector $\mathcal{Q}_{T}$, which is defined via (3.6) with $\textup{N}^{m}(T)=\Omega$, decays exponentially fast (measured in $m$) away from $T$. With a slightly different localization strategy, the procedure can also be interpreted in the spirit of an iterative domain decomposition solver, see [28, 29]. Since Assumption 3.1 holds for all $u^{*}\in H^{1}_{0}(\Omega)$, Proposition 4.1 is still valid if we replace $\mathcal{Q}$ and $\mathcal{Q}_{m}$ by $\mathcal{Q}^{u}$ and $\mathcal{Q}_{m}^{u}$, respectively.

We first discuss estimates for the Galerkin method (3.7) in Section 4.1. (Additional) error terms arise from the linearization, which we discuss separately in Section 4.2 together with the choice of $u^{*}$.

4.1 A priori error estimates

Since $V_{H,m}$ is a linear subspace of $H^{1}_{0}(\Omega)$, the Galerkin method (3.4) is automatically well defined, i.e., there exists a unique solution $u_{H,m}$ and its satisfies the following error estimate.

Theorem 4.2**.**

Let Assumptions 2.1 and 3.1 be fulfilled. Let $u$ be the solution to (2.2) and $u_{H,m}$ the solution to (3.7). Then it holds that

[TABLE]

as well as

[TABLE]

with, for any $v\in H^{1}_{0}(\Omega)$, the linearization error

[TABLE]

Further, if $A_{L}(x,\nabla u,\nabla u)=A(x,\nabla u)$, it also holds that

[TABLE]

Note that the assumption $A_{L}(x,\nabla u,\nabla u)=A(x,\nabla u)$ is satisfied for the Newton-type and the Kačanov-type linearization introduced in Example 3.2.

Proof of Theorem 4.2.

Since (3.7) defines a Galerkin method, the quasi-optimality (2.4) leads to

[TABLE]

Proof of (4.1): We choose $v_{H,m}=(\operatorname{id}-\mathcal{Q}_{m})I_{H}u=(\operatorname{id}-\mathcal{Q})I_{H}u+(\mathcal{Q}-\mathcal{Q}_{m})I_{H}u$ and observe that the second term can directly be estimated using Proposition 4.1, the stability of $I_{H}$ and (2.3). Note that by definition $u-(\operatorname{id}-\mathcal{Q})I_{H}u\in W$. Hence we obtain with the strong monotonicity of $A$ (cf. Assumption 2.1), the approximation property (2.6), and the definition of $\mathcal{Q}$ in (3.3) that

[TABLE]

Combination of this estimate with Proposition 4.1 leads to (4.1).

Proof of (4.2): We again choose $v_{H,m}=(\operatorname{id}-\mathcal{Q}_{m})I_{H}u=(\operatorname{id}-\mathcal{Q})I_{H}u+(\mathcal{Q}-\mathcal{Q}_{m})I_{H}u$, but we treat the first term in a slightly different manner. Namely, employing Assumption 3.1 and the definition of $\mathcal{Q}$ in (3.3), we deduce

[TABLE]

Combination of this estimate with Proposition 4.1 leads to (4.2).

Proof of (4.3): Again, we choose $v_{H,m}:=(\operatorname{id}-\mathcal{Q}_{m})I_{H}u$, but we split it in a different way this time. We write

[TABLE]

where we recall that $\mathcal{Q}^{u}$ and $\mathcal{Q}_{m}^{u}$ are the solutions to (3.3) and (3.6), respectively, with coefficient $\mathfrak{A}^{u}:=\mathbf{A}_{L}(x,\nabla u)$. The last term $(\mathcal{Q}_{m}^{u}-\mathcal{Q}_{m})I_{H}u$ is directly included in (4.3), while the second term $(\mathcal{Q}^{u}-\mathcal{Q}_{m}^{u})I_{H}u$ is again estimated with Proposition 4.1. For the first term $(\operatorname{id}-\mathcal{Q}^{u})I_{H}u$, we observe due to Assumption 3.1 that

[TABLE]

Since $u-(\operatorname{id}-\mathcal{Q}^{u})I_{H}u\in W$, the definition of $\mathcal{Q}^{u}$ in (3.3) implies that

[TABLE]

Hence, employing the assumption $A_{L}(x,\nabla u,\nabla u)=A(x,\nabla u)$ and (2.6), we deduce

[TABLE]

Up to the linearization errors, all variants of the previous theorem are identical to the linear elliptic case [32]. In particular, if we choose $m\approx|\log(H)|$, we have linear convergence without any assumptions on the regularity of $u$ or the variations of $A$. By Friedrich’s inequality, the same estimate also holds for the $L^{2}$-norm. In contrast to the linear case, the Aubin-Nitsche trick cannot be applied so that higher order convergence for nonlinear problems is rather difficult to achieve, see the discussion in [2]. Using the idea of the elliptic projection in [2], we obtain an $L^{2}$-estimate in Theorem A.1 in the Appendix. Roughly speaking, it yields quadratic convergence (up to (new) linearization errors) for the choice $m\approx|\log(H)|$.

By the stability of $I_{H}$ we deduce an estimate for the error to $I_{H}u_{H,m}$, which describes the finite element part of the Galerkin solution.

Corollary 4.3**.**

Let Assumptions 2.1 and 3.1 be fulfilled. Let $u$ be the solution to (2.2) and $u_{H,m}$ the solution to (3.7). Then it holds that

[TABLE]

Proof.

With the triangle inequality we split

[TABLE]

which finishes the proof together with the properties (stability, approximation, and projection) of $I_{H}$. ∎

Note that the two last terms in (4.4) can be estimated via Theorem 4.2. For the $L^{2}$-norm we also have the estimates from Theorem A.1 in the appendix. For $m\approx|\log(H)|$, the error $\|u-I_{H}u_{H,m}\|_{0}$ is at least of order $H$ and might be up to order $H^{2}$ if $u$ is sufficiently regular. If we assume $L^{2}$-stability of $I_{H}$, i.e., $\|I_{H}v\|_{0}\lesssim\|v\|_{0}$ for all $v\in H^{1}_{0}(\Omega)$, the estimate in Corollary 4.3 simplifies to

[TABLE]

Together with Theorem A.1 this implies that the error in the finite element part of $u_{H,m}$ is dominated by the $L^{2}$-best-approximation error in the finite element space.

4.2 Linearization errors and choice of linearization points

In this section, we discuss estimates for the linearization errors of Theorem 4.2, namely $\eta_{\mathrm{lin}}(v)$ and $|(\mathcal{Q}_{m}-\mathcal{Q}_{m}^{u})I_{H}u|_{1}$, and their implication on the choice of the linearization point $u^{*}$.

Linearization error $\eta_{\mathrm{lin}}(v)$.

First, we note that due to Assumptions 2.1 and 3.1, $\eta_{\mathrm{lin}}(v)$ can always be bounded as follows

[TABLE]

Hence, the stabilities of $\mathcal{Q}$ and $I_{H}$ as well as (2.3) directly imply that

[TABLE]

with $C(u^{*})=C_{0}+\|b_{L}(x,\nabla u^{*})\|_{0}$. Considering Example 3.2, we have $C(u^{*})=C_{0}$ for Kačanov-type linearizations and $C(u^{*})\leq 2C_{0}+(\Lambda+\Lambda_{L})|u^{*}|_{1}$ for Newton-type lineariaztions so that in both cases the linearization error is bounded by the data $f$ and $u^{*}$. As a consequence, if the nonlinearity and $\|f\|_{0}$ are small, the linearization error is negligible. This of course is a rather restrictive assumption since it basically means that we are still in the almost linear case with only a small nonlinear perturbation. Note that the above bound on $\eta_{\mathrm{lin}}(v)$ can be simplified in the case $v=(\operatorname{id}-\mathcal{Q})I_{H}u$ because due to (3.3)

[TABLE]

Next, we will show that $\eta_{\mathrm{lin}}(v)$ is small for the linearizations of Example 3.2 if $u^{*}$ is close to the exact solution $u$.

Lemma 4.4** (Linearization error for Kačanov-type linearization).**

Let $A(x,\nabla v)=\alpha(x,|\nabla v|^{2})\nabla v$ and set $\mathbf{A}_{L}(x,\nabla u^{*})=\alpha(x,|\nabla u^{*}|^{2})$ and $b_{L}(x,\nabla u^{*})=0$. Then,

[TABLE]

Proof.

Assumptions 2.1 and 3.1 directly yield

[TABLE]

Lemma 4.5** (Linearization error for Newton-type linearization).**

Set $\mathbf{A}_{L}(x,\nabla u^{*})=D_{\xi}A(x,\nabla u^{*})$ and $b_{L}(x,\nabla u^{*})=A(x,\nabla u^{*})-D_{\xi}A(x,\nabla u^{*})\nabla u^{*}$. Assume that $D_{\xi}A$ is Lipschitz continuous in its last argument, i.e., there is $L_{A}>0$ such that

[TABLE]

Then,

[TABLE]

Proof.

We perform a Taylor expansion of $A$ around $u^{*}$ and obtain that for any $w\in H^{1}_{0}(\Omega)$ it holds that

[TABLE]

The assumed Lipschitz continuity of $D_{\xi}A$ (4.5) finishes the proof. ∎

Note that we did not exploit $w\in W$ in the estimates for $\eta_{\mathrm{lin}}(v)$. Lemmas 4.4 and 4.5 underline that $u^{*}$ should be close to $u$, which can be seen as a local convergence result and motivates the discussion of different linearization points further below in this section.

Linearization error $|(\mathcal{Q}_{m}-\mathcal{Q}_{m}^{u})I_{H}u|_{1}$.

In case of a Kačanoc-type linearization from Example 3.2, we observe that $\eta_{\mathrm{lin}}(u)$ in (4.2) can be bounded as $\eta_{\mathrm{lin}}(u)\leq\|\mathfrak{A}^{u}-\mathfrak{A}\|_{L^{\infty}(\Omega)}|u|_{1}$. Hence, this linearization error can be expressed as the $L^{\infty}(\Omega)$-error between the “true” coefficient $\mathfrak{A}^{u}$ and the “reference” coefficient $\mathfrak{A}$ used in practice for computing the corrector $\mathcal{Q}_{m}$. The error $|(\mathcal{Q}_{m}-\mathcal{Q}_{m}^{u})I_{H}u|_{1}$ in (4.3) is obviously also closely related to the error between the coefficients $\mathfrak{A}^{u}$ and $\mathfrak{A}$. However, estimating directly the error between the correction operators allows for a refined estimate. Before we present this in detail, let us first note that, obviously,

[TABLE]

so that we will estimate $(\mathcal{Q}_{m}-\mathcal{Q}_{m}^{u})v_{H}$ for an arbitrary $v_{H}\in V_{H}$ in the following. Further, by the definition of $\mathcal{Q}_{m}$ and $\mathcal{Q}_{m}^{u}$ as sums of element correctors with support only in $\textup{N}^{m}(T)$ we have

[TABLE]

We now have the following estimate for the error between the correction operators, the proof is postponed to the Appendix B.

Proposition 4.6**.**

Fix $T\in\mathcal{T}_{H}$ and $v_{H}\in V_{H}$ and recall the notation $\textup{N}^{m}(T)$ from Section 3.2. For any $a\in L^{\infty}(\textup{N}^{m}(T),\mathbb{R}^{d\times d})$ denote by $\overline{a}$ the average of the trace, i.e., $\overline{a}:=\fint_{\textup{N}^{m}(T)}\mathrm{tr}(a)\,dx$, and denote by $\widehat{a}\in L^{\infty}(\textup{N}^{m}(T),\mathbb{R}^{d\times d})$ the scaled coefficient, i.e., \widehat{a}:=a\big{/}\overline{a}. Recall $\mathfrak{A}=\boldsymbol{A}_{L}(\cdot,\nabla u^{*})$ with associated corrector $\mathcal{Q}_{T,m}$ and $\mathfrak{A}^{u}=\boldsymbol{A}_{L}(\cdot,\nabla u)$ with associated corrector $\mathcal{Q}_{T,m}^{u}$, cf. Assumption 3.1. Define

[TABLE]

Then,

[TABLE]

We emphasize that Proposition 4.6 relates the error between the element correctors to the error between the coefficients $\mathfrak{A}$ and $\mathfrak{A}^{u}$, but (i) locally on each element or element patch, respectively, and (ii) only the $L^{\infty}$-error between the scaled coefficients is relevant. The latter point is possible because (3.6) can be multiplied by the scalar-valued constant $\overline{\mathfrak{A}}$ without changing the element corrector. Hence, if a coefficient is (locally) multiplied by a constant, the error between the associated correction operators is zero while the error between the (unscaled) coefficients themselves can be very large.

Error indicators similar to $E_{\mathcal{Q},T}$ in Proposition 4.6 were already presented in [19, 18] in the context of time-changing or perturbed diffusion coefficients with the following differences. First, the idea of studying the scaled coefficients is new in Proposition 4.6. Second, since $\mathfrak{A}^{u}$ is not available in practice in our case, Proposition 4.6 is an a priori result and aims at linking the error between the correction operators to the error between the coefficients. This is in sharp contrast to the previous works, which use $E_{\mathcal{Q},T}$ as a practical indicator to steer re-computation of element correctors. For this reason, some terms in [18, 19] are (slightly) different to obtain contrast-independent bounds. The result of Proposition 4.6 can probably be refined in this spirit as well, but we omit this for simplicity.

Choice of linearization points $u^{*}$.

The previous discussion has shown that the choice of $u^{*}$ is crucial for the performance of the method and that, ideally, $u^{*}$ should be chosen close the exact solution $u$. With this in mind, let us compare some possible choices of $u^{*}$ also used in the numerical experiments.

We can select for $u^{*}$ the initial value of the nonlinear iteration, for instance $u^{*}=0$. This results in a multiscale space $V_{H,m}$ that can be computed a priori and yields a very cheap method. The overall error of course will only be dominated by the discretization error of order $H+\beta^{m}$ if $u^{*}$ is already close to the exact solution. For instance, a similar assumption is also needed to guarantee convergence of Newton’s iteration method for the nonlinear problem. In particular the choice $u^{*}=0$ should yield satisfactory results for problems without steep gradients. 2. 2.

Another option is to compute the finite element solution $u_{H}$ using $V_{H}$ and to use it as $u^{*}$. Since $V_{H}$ is associated with a coarse mesh, the computation of $u_{H}$ is rather cheap. Due to the finescale features of $A$, however, $u_{H}$ alone will not be a satisfactory approximation of $u$ in general, not even in the $L^{2}(\Omega)$-norm. The finescale features of $A$ are then taken into account in the ensuing LOD solution. To estimate the linearization error in this case, we can combine Lemmas 4.4 or 4.5, respectively, with standard a priori estimates for the finite element solution, see Section 2.2 for estimates in the $L^{2}(\Omega)$-norm and $H^{1}(\Omega)$-semi norm. In case of the Newton-type linearization, Lemma 4.5 also requires a $W^{1,\infty}(\Omega)$-estimate, where we obtain from [8, Chapter 8]

[TABLE]

Note that for all finite element error estimates, higher regularity of $u$ is required and that the constants will depend on the finescale parameter $\varepsilon$. Computing the LOD solution $u_{H,m}$ with $u^{*}=u_{H}$ is notably different from the two-grid approach [40] where the FE solution $u_{H}$ is “post-processed” by another (linear) FE solve on a finer mesh. 3. 3.

Finally, one can compute a whole cascade of LOD solutions by starting with some $u^{*}$, computing the solution $u_{H,m}$ to (3.7), and using this solution as new $u^{*}$. For this cascade, we can inductively combine (4.2) and Lemmas 4.4 or 4.5, respectively, to obtain an a priori error estimate, which consists of the discretization error of order $H+\beta^{m}$ and the linearization error for the initial choice of $u^{*}$. This straightforward error estimate, however, does not seem to be optimal because the initial linearization error will always remain and one does not exploit that the LOD solution should be closer to the exact solution with each step. For this cascade of LOD solutions, an error indicator in the spirit of Proposition 4.6 – with $\mathcal{Q}_{m}^{u}$ replaced by $\mathcal{Q}_{m}^{u_{H,m}}$ with the current LOD solution – can be used to locally determine which element correctors should be recomputed with the new linearization point. Since (3.7) is solved by an iterative method as well, one can also imagine to already update some element correctors during this nonlinear iteration. Such an iterative linearized approach is studied for the nonlinear Helmholtz equation in [31], but the application to quasilinear problems is a future research topic of its own.

5 Numerical experiments

We present the results of several numerical experiments, subject to different multiscale coefficients and nonlinearities. In all cases, the computational domain is $\Omega=[0,1]^{2}$. Since no exact solution is known in any of the examples, we compute a reference solution $u_{h}\in V_{h}$ using a standard finite element method on a fine mesh $\mathcal{T}_{h}$ that resolves all multiscale features. Specifically, we fix $h=2^{-8}$ and vary $H=2^{-2},2^{-3},\ldots 2^{-6}$. We present results for oversampling parameters $m=1,2,3$ and already point out that $m=2,3$ is a sufficient choice in most experiments. To solve the nonlinear problems, we use Newton’s method with tolerance $10^{-11}$ for the residual as stopping criterion. Sections 5.1 and 5.2 consider nonlinearities of the type $A(x,\nabla u)$, where Assumption 2.1 is only satisfied in Section 5.1. Further, we consider a model for the stationary Richards equation with a quasilinear coefficient of the form $a(x,u)\nabla u$ in Section 5.3. Note that the resulting nonlinear form $\mathcal{B}$ is no longer monotone such that the above proof techniques do not directly transfer, see [5, 6]. Nevertheless, we use the presented multiscale method with the obvious modifications for (3.7) and $\mathfrak{A}:=a(x,0)$ in the corrector problems. The code is available at Zenodo with doi 10.5281/zenodo.4311614 and is based on a preliminary implementation of the LOD for linear problems developed at the Chair of Computational Mathematics, University of Augsburg cf. [33].

5.1 Periodic coefficient

We choose a model problem similar to [2] with nonlinearity

[TABLE]

with $\varepsilon=2^{-5}$ and sources $f_{1}(x)=10\exp(-0.1|x-x_{0}|^{2})$ or $f_{2}(x)=100\exp(-0.1|x-x_{0}|^{2})$ with $x_{0}=(0.45,0.5)^{T}$. The coefficient and the reference solutions $u_{h}$ are depicted in Figure 5.1. We will consider the two right-hand sides $f_{1}$ and $f_{2}$ to study the influence of higher values of the solution $u$ and its gradient on the errors. Note that one can hope for higher regularity of the exact solution $u$ in this experiment.

We use the Galerkin method (3.7) and obtain a solution $u_{H,m}\in V_{H,m}$, where the correctors are computed using the Newton-type linearization at $u^{*}$, which is specified below. We focus on two (relative) errors in the following: The so-called (relative) upscaled error

[TABLE]

for which we expect a linear convergence rate (cf. Theorem 4.2); and the so-called (relative) macroscopic error

[TABLE]

for which we expect the same behavior as the $L^{2}$-best approximation in $V_{H}$ (cf. Corollary 4.3).

For the simple choice $u^{*}=0$, these two errors are depicted for the two right-hand sides in Figure 5.2. We note that the (relative) macroscopic errors $e_{H}$ in the left column closely follow the error of the (relative) $L^{2}$-best approximation in the space $V_{H}$ for $m=2,3$ (cf. the discussion after Corollary 4.3). Since $I_{H}u_{H,m}$ lies in the same space, we cannot hope for anything better. In particular, the reduced convergence rate for $H$ between $\sqrt{\varepsilon}$ and $\varepsilon$ is no defect of the method, but intrinsic to the problem, see also the discussion of this so-called resonance effect in [17]. We emphasize that in the pre-asymptotic range $\varepsilon\ll H$, the standard finite element method shows no convergence rates in contrast to the multiscale method. For approximations in the $H^{1}$-semi norm, the (coarse-scale) space $V_{H}$ is no longer sufficient. Therefore, we consider the (relative) upscaled error $e_{\operatorname{LOD}}$ in the right column of Figure 5.2. This error overall converges linearly as expected from Theorem 4.2. All in all, the experiment clearly confirms the predicted convergence rates of Theorem 4.2 and Corollary 4.3.

When comparing the top and bottom row of Figure 5.2, we also observe that the larger gradients of $u_{h}$ caused by $f_{2}$ influence the performance of the method. In particular $e_{\mathrm{LOD}}$ shows some deviation from the optimal linear convergence. Therefore, we compare the behavior of $e_{H}$ and $e_{\mathrm{LOD}}$ for different linearization points with fixed $m=3$ in Figure 5.3. We consider the following choices of $u^{*}$ as discussed in Section 4.2: $u^{*}=0$, $u^{*}=u_{H}$ with the (coarse) FE solution $u_{H}\in V_{H}$ and $u^{*}=u_{H,m}$ as well as $u^{*}=I_{H}u_{H,m}$, where $u_{H,m}$ is the LOD solution with linearization at zero. All these solutions show a similar qualitative behavior, but in the quantitative errors we observe clear differences, especially for $e_{\mathrm{LOD}}$. The choices $u^{*}=u_{H}$ and $u^{*}=I_{H}u_{H,m}$ result in almost identical results because they lie in the same space $V_{H}$. A bit surprisingly, those choices even perform slightly worse than $u^{*}=0$. A possible explanation is that the gradients are relevant for the nonlinearity, where the space $V_{H}$ does not provide sufficient approximations. Both for $e_{H}$ and $e_{\mathrm{LOD}}$, the choice $u^{*}=u_{H,m}$ performs best, which motivates to study iterative LOD approximations as discussed in Section 4.2 in future research.

5.2 Random coefficient

We choose the nonlinear coefficient as

[TABLE]

where $c(x)$ is piece-wise constant on a quadrilateral mesh $\mathcal{T}_{\varepsilon}$ with $\varepsilon=2^{-6}$ and the values are random numbers in $[0.1,1]$. The right-hand side is

[TABLE]

[[, see]]Henn11phdhmm for the nonlinearity and a similar right-hand side. Note that we clearly cannot expect higher regularity for the exact solution because of the spatial discontinuities in $A$. Further, $A$ is only locally Lipschitz constant in its second argument. Thus, Assumption 2.1 is violated and we even have $|A|\to\infty$ for $|\xi|\to\infty$. The coefficient and corresponding reference solution $u_{h}$ are depicted in Figure 5.4.

For this example, we study the Petrov-Galerkin LOD as briefly discussed in Remark 3.3. Precisely, we compute $u_{H,m}^{PG}\in V_{H}$ as the solution of

[TABLE]

where $V_{H,m}$ is defined as in the Galerkin case. Hence, the solution $u_{H,m}^{PG}$ lies in the FE space and we expect the relative macroscopic error $e_{H}:=\frac{\|u_{h}-u_{H,m}\|_{0}}{\|u_{h}\|_{0}}$ to follow the $L^{2}$-best approximation, cf. Corollary 4.3.

We first consider the convergence history of $e_{H}$ for $u^{*}=0$ and different choices of $m$ in Figure 5.5 (left). The multiscale method performs obviously better than the standard FEM. The error of the Petrov-Galerkin LOD is following the $L^{2}$-best approximation as expected up to a saturation or stagnation for the last considered mesh. Since this effect occurs for all $m=1,2,3$ and was also observed for $m=4,5$ (data not shown), the linearization error most probably starts to dominate. As in the previous example, we also consider $e_{H}$ for fixed $m$ and different $u^{*}$ in Figure 5.5 (right). The considered choices are $u^{*}=0$, $u^{*}=u_{H}$ with the FE solution $u_{H}$ on the coarse mesh as before and $u^{*}=u_{H,m}^{PG}$ as well as $u^{*}=(\operatorname{id}-\mathcal{Q}_{m})u_{H,m}^{PG}$, where $u_{H,m}^{PG}$ is the Petrov-Galerkin LOD solution with $u^{*}=0$ and $\mathcal{Q}_{m}$ is the corrector computed with $u^{*}=0$. The latter two linearization points play roles comparable to $I_{H}u_{H,m}$ and $u_{H,m}$ in the Galerkin setting discussed in the previous experiment. Only $u^{*}=0$ shows the saturation at the end, but nevertheless, this choice of linearization performs best. We emphasize that the impractical choice of $u^{*}=u_{h}$ would lead to an $e_{H}$ completely following the $L^{2}$-best approximation error (data not shown). This illustrates the validity of our error estimates also for the Petrov-Galerkin LOD, but also stresses the influence of the chosen linearization. Moreover, this example confirms and underlines that the present multiscale method does not rely on assumptions such as periodicity or scale separation.

5.3 Stationary Richards equation

We now test the applicability of our method to quasilinear non-monotone problems, which, for instance, are frequently encountered in (unsaturated) groundwater flow and can be modeled by the (stationary) Richards equation. Here, we consider a quasilinear coefficient of the form $A(x,u,\nabla u)=c(x)k(u)\nabla u$. For $c$ we choose a spatial multiscale model with a channel as depicted in Figure 5.6 left, which is often present in geophysical applications. Note that this can easily be extended to the case of several channels. For the nonlinearity $k$, we consider the following so-called van Genuchten model [39]

[TABLE]

with $\alpha=0.005$, [[, see also]]EGKL12nonlinflowitsolver. The right-hand side is

[TABLE]

The reference solution for the present setting is depicted in Figure 5.6 right. Note the influence of the channel, which analytically manifests itself in a low regularity of the solution.

As a consequence, we observe exactly the (worst-case) convergence rates as predicted by our theory (for the monotone case), but not more. In particular, $e_{H}$ follows the best-approximation error, which in this case is “only” linear as discussed after Corollary 4.3. Furthermore, we see the expected linear convergence of the upscaled error $e_{\operatorname{LOD}}$ up to a slight saturation for small $H$. This is most probably caused again by a dominating linearization error, which could be cured by a different choice of the linearization. This experiment clearly underlines the applicability of the approach also beyond the strictly monotone case and indicates that similar convergence rates as in Theorem 4.2 and Corollary 4.3 can be expected.

Conclusion

We presented a multiscale method for nonlinear monotone elliptic problems with spatial multiscale features. A problem-adapted multiscale basis is constructed by solving local linear fine-scale problems for each coarse-scale mesh element. Numerical analysis shows optimal error estimates up to linearization errors and we discussed choices of the linearization. Several numerical experiments underline and confirm the applicability of the method as well as the expected convergence rates. We also numerically compared the influence of the chosen linearization on the performance of the method. As mentioned, iterative (adaptive) LOD methods where either a cascade of LOD solutions to different linearization points is computed or the correctors are (partially) updated during the nonlinear iteration are interesting extensions of the presented method and will be subject of future research. Also the numerical analysis for the quasilinear non-monotone problems will be studied in the future.

Acknowledgments

We are grateful to D. Peterseim for fruitful discussion on the subject and for providing a preliminary implementation of the LOD for linear problems. We thank the anonymous reviewers for their valuable remarks.

Appendix A $L^{2}$-error estimate for the Galerkin method

In this appendix, we prove an $L^{2}$-estimate for the Galerkin method (3.4).

Theorem A.1**.**

Let Assumptions 2.1 and 3.1 be fulfilled and suppose that $A(x,\nabla u)=A_{L}(x,\nabla u,\nabla u)$. Let $u$ be the solution to (2.2) and $u_{H,m}$ the solution to (3.7). Then it holds that

[TABLE]

with the linearization error $\eta_{\mathrm{lin},1}$ defined as

[TABLE]

and where $\widetilde{u}_{H,m}$ is defined via (A.1) below.

The first term of the above error estimate presumably is of $O(H^{2})$ for the choice $m\approx|\log H|$, see Section 4.1. The linearization error $\eta_{\mathrm{lin},1}(u,v)$ can be estimated in the same way as in Lemmas 4.4 and 4.5 by replacing $u^{*}$ with $u$. We then need to estimate $|u-\widetilde{u}_{H,m}|_{1}$ and $\|\nabla(u-\widetilde{u}_{H,m})\|_{L^{\infty}(\Omega)}$. As shown below, $|u-\widetilde{u}_{H,m}|_{1}$ is of the same order as $|u-u_{H,m}|_{1}$ discussed in Theorem 4.2. However, it is not clear whether $W^{1,\infty}$-estimates for the multiscale space $V_{H,m}$ can be established in the spirit of [8, Chapter 8] to obtain a quadratic rate for $\eta_{\mathrm{lin},1}$ in case of a Newton-type linearization.

Proof.

Let $\widetilde{u}_{H,m}\in V_{H,m}$ be the unique solution to

[TABLE]

cf. [2]. Note that this can equivalently be written as

[TABLE]

We split the error $u-u_{H,m}$ into $u-\widetilde{u}_{H,m}$ and $\widetilde{u}_{H,m}-u_{H,m}$ and estimate both parts separately.

First step: Estimate of $\widetilde{u}_{H,m}-u_{H,m}$: By the monotonicity of $A$, Galerkin orthogonality, the assumption $A(x,\nabla u)=A_{L}(x,\nabla u,\nabla u)$ and (A.1), we deduce

[TABLE]

The $L^{2}$-estimate is obtained by applying Friedrich’s inequality.

Second step: Estimate of $u-\widetilde{u}_{H,m}$: This is in principle an $L^{2}$-estimate for a Galerkin LOD for an elliptic diffusion problem. The main issue, however, is that the multiscale space $V_{H,m}$ is not built with respect to the diffusion tensor $\mathfrak{A}^{u}=\mathbf{A}_{L}(x,\nabla u)$ but with respect to $\mathfrak{A}=\mathbf{A}_{L}(x,\nabla u^{*})$. We first of all note that due to the projection, we have

[TABLE]

Let now $z\in H^{1}_{0}(\Omega)$ and $z_{H,m}\in V_{H,m}$ be the solutions of the following dual problems

[TABLE]

Assumption 3.1 and Galerkin orthogonality imply

[TABLE]

where the second term is estimated with Proposition 4.1. Employing (2.5) and the a priori (stability) estimate for $z$, the last term yields

[TABLE]

For the first term, we obtain with the ellipticity of $\mathfrak{A}^{u}$ as well as the definitions of $\mathcal{Q}^{u}$ and $z$ that

[TABLE]

Combining the foregoing estimates, we conclude

[TABLE]

Finally, the definition of the dual problems yields

[TABLE]

which in combination with the already derived estimates finishes the proof. ∎

Appendix B Proof of Proposition 4.6

The proof of Proposition 4.6 simply relies on the definition of the element correctors and is similar to the results in [18, 19].

Proof of Proposition 4.6.

Let $T\in\mathcal{T}_{H}$ and $v_{H}\in V_{H}$ be fixed. Abbreviate $w:=(\mathcal{Q}_{T,m}-\mathcal{Q}_{T,m}^{u})v_{H}$ and note that $w\in W(\textup{N}^{m}(T))$. We deduce by the definition of $\mathcal{Q}_{T,m}$ and $\mathcal{Q}_{T,m}^{u}$ that

[TABLE]

We then proceed as follows

[TABLE]

which finishes the proof. Note that only in the very first step we hide the lower spectral bound of $\widehat{\mathfrak{A}^{u}}|_{\textup{N}^{m}(T)}$ in the notation $\lesssim$, all other estimates are constant-free. ∎