×

Numerical homogenization: survey, new results, and perspectives. (English. French summary) Zbl 1329.65300

Summary: These notes give a state of the art of numerical homogenization methods for linear elliptic equations. The guideline of these notes is analysis. Most of the numerical homogenization methods can be seen as (more or less different) discretizations of the same family of continuous approximate problems, which H-converges to the homogenized problem. Likewise numerical correctors may also be interpreted as approximations of Tartar’s correctors. Hence the convergence analysis of these methods relies on the H-convergence theory. When one is interested in convergence rates, the story is different. In particular one first needs to make additional structure assumptions on the heterogeneities (say periodicity for instance). In that case, a crucial tool is the spectral interpretation of the corrector equation by Papanicolaou and Varadhan. Spectral analysis does not only allow to obtain convergence rates, but also to devise efficient new approximation methods. For both qualitative and quantitative properties, the development and the analysis of numerical homogenization methods rely on seminal concepts of the homogenization theory. These notes contain some new results.

MSC:

65N99 Numerical methods for partial differential equations, boundary value problems
35B27 Homogenization in context of PDEs; PDEs in media with periodic structure
35J25 Boundary value problems for second-order elliptic equations
Full Text: DOI

References:

[1] ij:= \partialj
[2] i- \partiali
[3] j, then the product u{\(\epsilon\)}{\(\cdot\)} v{\(\epsilon\)}converges to u0{\(\cdot\)} v0in the sense of distributions. In view of these results, a natural candidate for the limit of A\ast{\(\epsilon\)},Has {\(\epsilon\)} and H go to zero is A0. In the following section we show how H-convergence can be used to prove the convergence of the self-consistent approach. Throughout the text, we’ll make use of the following notation ●d is the space dimension ; ●D is a bounded open Lipschitz domain of Rd; ●Mddenotes the set of d-dimensional real square matrices ; ESAIM: PROCEEDINGS57 ●M{\(\alpha\)}{\(\beta\)}denotes the set of d-dimensional real square matrices which are {\(\alpha\)}-elliptic and {\(\beta\)}-continuous ;symis the subset of those symmetric matrices of M{\(\alpha\)}{\(\beta\)}; ●M{\(\alpha\)}{\(\beta\)} ●for all {\(\rho\)}, Q{\(\rho\)}= (-{\(\rho\)}/2, {\(\rho\)}/2)d, and we use the short hand notation Q = Q1= (-1/2, 1/2)d; ●for all x \in Rd, Txdenotes the translation by x and for every measurable subset B of Rd, TxB = {x + y : y \in B} ; ●for all x \in Rd, and all {\(\rho\)} > 0, Q{\(\rho\)}(x) := TxQ{\(\rho\)}; ●H1per(Q) denotes the closure of smooth Q-periodic function with zero average in the Hilbert space H1(Q) ; ●for all 1 \leq p \leq \infty, W1,p(D) denotes the Sobolev space of p-integrable functions whose distributional derivatives are p-integrable functions ; ●for all 1 \leq p \leq \infty, W1,p0(D) denotes the subspace of functions W1,p(D) which vanishe on \partialD in the sense of traces ; ●{\(\cdot\)} is the ensemble average, that is the periodic average in the periodic case, and the expectation in the random case ; ●var [{\(\cdot\)}] is the variance associated with the ensemble average ; \bulletandstand for \leq and \geq up to a multiplicative constant which only depends on the dimension d and the coercivity constants (denoted by {\(\alpha\)}, {\(\beta\)} in the text) if not otherwise stated; ●when bothandhold, we simply write \sim; ●we use \gg instead ofwhen the multiplicative constant is (much) larger than 1;d. ●(e1, . . . , ed) denotes the canonical basis of R 2.Analytical framework by H-convergence In this section we present an analytical framework to analyze the convergence of numerical homogenization methods in the case of linear elliptic equations in divergence form. These results are proved using a simplified version of the string of arguments used in [27] to treat the case of general multiple integrals. In addition they cover the case of non symmetric matrices (which was not treated in [27]). 2.1. General framework Let A{\(\epsilon\)}\in M{\(\alpha\)}{\(\beta\)}(D) be a H-convergent sequence whose limit is denoted by Ahom\in M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)}(D). Unlike what we’ve presented in the self-consistent approach, we focus here on a continuous approximation, and shall only later on discretize the equations. We begin with the definition of a local approximation of Ahomon domains of size {\(\rho\)} > 0. Definition 3.For all {\(\rho\)} > 0 and {\(\epsilon\)} > 0, we denote by A{\(\rho\)},{\(\epsilon\)}the element of M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)}(D) defined by: for all i, j \in {1, . . . , d} and for x \in D,
[4] ij:=ej{\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y))dy,(2.1) Q{\(\rho\)}\capT-xD where v{\(\rho\)},{\(\epsilon\)}i(x, {\(\cdot\)}) is the unique weak solution in H10(Q{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y)) = 0in Q{\(\rho\)}\cap T-xD.(2.2) These approximations A{\(\rho\)},{\(\epsilon\)}of Ahomare similar to the coefficients A\ast{\(\epsilon\)},Hof the self-consistent approach. The fact that A{\(\rho\)},{\(\epsilon\)}\in M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)}(D) is proved as follows.{\(\rho\)},{\(\epsilon\)}yields The weak formulation of (2.2) tested with function vi \hat{}\hat{} (ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y))dy =ei{\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y))dy. Since A{\(\epsilon\)}\in M{\(\alpha\)}{\(\beta\)}(D), by Cauchy-Schwarz inequality this turns into {\(\alpha\)}|ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y)|2dy \leq {\(\beta\)}|Q{\(\rho\)}\cap T-xD|1/2Q{\(\rho\)}\capT-xD|ei+ \nablayv{\(\rho\)},{\(\epsilon\)}i(x, y)|2dy,(2.3) \hat{}\hat{}1/2 Q{\(\rho\)}\capT-xD from which the upper bound follows using the defining equation (2.1). We turn to the lower bound. For all {\(\xi\)} \in Rd, we let v{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, {\(\cdot\)}) be the weak solution in H10(Q{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablayv{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)) = 0in Q{\(\rho\)}\cap T-xD. Using the lower bound on A{\(\epsilon\)}and Jensen’s inequality, we have for all {\(\xi\)} \in Rdwith |{\(\xi\)}| = 1 \hat{}{\(\rho\)},{\(\epsilon\)}(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablayv{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y))dy Q{\(\rho\)}\capT-xD{\(\xi\)} ({\(\xi\)} + \nablayv \geq {\(\alpha\)}|{\(\xi\)} + \nablayv{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y))|2dy \hat{} Q{\(\rho\)}\capT-xD \geq {\(\alpha\)}|Q{\(\rho\)}\cap T-xD|,(2.4) which is the desired lower bound since {\(\xi\)} {\(\cdot\)} A{\(\rho\)},{\(\epsilon\)}(x){\(\xi\)} =({\(\xi\)} + \nablayv{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablayv{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y))dy. Q{\(\rho\)}\capT-xD If {A{\(\epsilon\)}} is a family of symmetric matrices, (2.2) is the Euler-Lagrange equation associated with the followingd, equivalent definition of (2.1): for all {\(\xi\)} \in R {\(\xi\)} {\(\cdot\)} A{\(\rho\)},{\(\epsilon\)}(x){\(\xi\)} := inf({\(\xi\)} + \nablav(y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablav(y))dy, v \in H10(Q{\(\rho\)}\cap T-xD) . Q{\(\rho\)}\capT-xD The main result of this section is the following theorem. Theorem 2.Let A{\(\epsilon\)}and A{\(\rho\)},{\(\epsilon\)}be as in Definition 3, then for all {\(\rho\)} > 0 there exists A{\(\rho\)},hom\in M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)}(D) such that for almost every x \in D, lim{\(\epsilon\)}\to0A{\(\rho\)},{\(\epsilon\)}(x)= A{\(\rho\)},hom(x),(2.5) lim{\(\rho\)}\to0A{\(\rho\)},hom(x)= Ahom(x).(2.6) As a direct corollary we have Corollary 1.Let A{\(\epsilon\)}, A{\(\rho\)},{\(\epsilon\)}and A{\(\rho\)},hombe as in Theorem 2, and f \in H-1(D). Then, the weak solution1(D) to u{\(\rho\)},{\(\epsilon\)}\in H0 -\nabla {\(\cdot\)} A{\(\rho\)},{\(\epsilon\)}\nablau{\(\rho\)},{\(\epsilon\)}= f satisfies lim{\(\rho\)}\to0lim{\(\epsilon\)}\to0u{\(\rho\)},{\(\epsilon\)}- uhom H1(D)= 0,(2.7) where uhom\in H10(D) is the weak solution to -\nabla {\(\cdot\)} Ahom\nablauhom= f. ESAIM: PROCEEDINGS59 As a consequence of H-convergence we also have that lim{\(\rho\)}\to0lim{\(\epsilon\)}\to0u{\(\rho\)},{\(\epsilon\)}- u{\(\epsilon\)} L2(D)= 0, Let us point out that without any further assumption on A{\(\epsilon\)}, one cannot get quantitative convergence rates for (2.7). A trivial example is provided by a constant family: A{\(\epsilon\)}:= Ahomfor all {\(\epsilon\)} > 0. In this case, if Ahom: Rd\to Mdis Lipschitz continuous, then the convergence rate in (2.7) is O({\(\rho\)}). Remark 1.Corollary 1 also holds with general Dirichlet, Neumann, mixed Dirichlet-Neumann boundary conditions. In the following subsection we shall complete Corollary 1 with a corrector result in order to approximate correctly \nablau{\(\epsilon\)}in L2(D, Md). The proof of Theorem 2 relies on three ingredients: ●the definition of H-convergence for (2.5), ●the approximate continuity of integrable functions (see (2.9) below), ●the continuous dependence of solutions to linear elliptic problems with respect to the coefficients of the operator stated in the following lemma. Lemma 3.Let A \in M{\(\alpha\)}{\(\beta\)}(D), (A{\(\rho\)}){\(\rho\)}>0\in M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)}(D) and f, (f{\(\rho\)}){\(\rho\)}>0\in H-1(D) be such that A{\(\rho\)}\to A pointwise in D, and f{\(\rho\)}\to f in H-1(D) as {\(\rho\)} goes to zero. Then the unique weak solution u{\(\rho\)}\in H10(D) to -\nabla {\(\cdot\)} A{\(\rho\)}\nablau{\(\rho\)}= f{\(\rho\)} converges in H1(D) to the unique weak solution u in H10(D) to -\nabla {\(\cdot\)} A\nablau = f. We first prove Theorem 2 and Corollary 1, and then turn to the proof of Lemma 3. Proof of Theorem 2. Let x \in D and {\(\rho\)} > 0, and consider problem (2.2). By the locality and definition of H-convergence (see property (2) of Lemma 1 and Definition 1), v{\(\rho\)},{\(\epsilon\)}i(x, {\(\cdot\)}) \rightharpoonup v{\(\rho\)},homi(x, {\(\cdot\)}) in H10(Q{\(\rho\)}\cap T-xD),{\(\rho\)},{\(\epsilon\)}(x, {\(\cdot\)})) \rightharpoonup Ahom(x + {\(\cdot\)})(ei+ \nablayv{\(\rho\)},homi(x, {\(\cdot\)})) in L2(Q{\(\rho\)}\cap T-xD, Rd),(2.8) A{\(\epsilon\)}(x + {\(\cdot\)})(ei+ \nablayvi where v{\(\rho\)},homi(x, {\(\cdot\)}) is the unique solution in H10(Q{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} Ahom(x + y)(ei+ \nablayv{\(\rho\)},homi(x, y)) = 0. Hence, setting
[5] ij:=ej{\(\cdot\)} Ahom(x + y)(ei+ \nablayv{\(\rho\)},homi(x, y))dy, Q{\(\rho\)}\capT-xD (2.8) implies the the claim (2.5). To prove (2.6), we appeal to Lemma 3. To this aim, we note that for all {\(\rho\)} small enough, Q{\(\rho\)}(x) \subset D, so that after a change of variables
[6] ij=ej{\(\cdot\)} Ahom(x + {\(\rho\)}y)(ei+ \nablayv{\(\rho\)},homi(x, {\(\rho\)}y))dy. By the continuity of translations in L1(D) (see [56] or [24] for instance), since Ahom\in L1(D), for all y \in Q = (-1/2, 1/2)dand B \subset\subset D we have \hat{}{\(\rho\)}\to0\to 0. B |Ahom(x + {\(\rho\)}y) - Ahom(x)|dx Integrating over Q and using Fubini’s theorem, one obtains \hat{}\hat{}{\(\rho\)}\to0\to 0.(2.9) BQ |Ahom(x + {\(\rho\)}y) - Ahom(x)|dy dx Consequently, for almost every x \in B, and almost every y \in Q, Ahom(x + {\(\rho\)}y){\(\rho\)}\to0\to Ahom(x).(2.10) Let now x \in B be such a point, and let then w{\(\rho\)}i\in H10(Q) be solutions for i \in {1, . . . , d} to -\nablay{\(\cdot\)} Ahom(x + {\(\rho\)}y)\nablayw{\(\rho\)}i(y) = \nablay{\(\cdot\)} Ahom(x + {\(\rho\)}y)ei. Estimate (2.10) implies that the assumptions of Lemma 3 are satisfied, so that w{\(\rho\)}i\to wiin H1(Q), where wiis10(Q) to the unique weak solution in H -\nablay{\(\cdot\)} Ahom(x)\nablaywi(y) = \nablay{\(\cdot\)} Ahom(x)ei= 0.(2.11) Hence, for all i, j \in {1, . . . , d}
[7] ij=ej{\(\cdot\)} Ahom(x + {\(\rho\)}y)(\nablayw{\(\rho\)}i(y) + ei)dy \hat{} C(0,1) {\(\rho\)}\to0\to\hat{}ej{\(\cdot\)} Ahom(x)(\nablaywi(y) + ei)dy C(0,1) =
[8] ij since wi= 0 is the trivial solution to (2.11). This concludes the proof of the theorem. We now prove Corollary 1. Proof of Corollary 1. Let u{\(\rho\)},hombe the unique weak solution in H10(D) to -\nabla {\(\cdot\)} A{\(\rho\)},hom\nablau{\(\rho\)},hom= f. Due to (2.5), Lemma 3 implies lim{\(\epsilon\)}\to0u{\(\rho\)},hom- u{\(\rho\)},{\(\epsilon\)} H1(D)= 0.(2.12) Similarly, from (2.6) we get lim{\(\rho\)}\to0u{\(\rho\)},hom- uhom H1(D)= 0.(2.13) The claim follows from the combination of (2.12) and (2.13) We conclude this subsection by the proof of Lemma 3. ESAIM: PROCEEDINGS61 Proof of Lemma 3. Let us substract the weak forms of the two equations tested against the admissible testfunction u{\(\rho\)}- u \in H10(D). This yields \hat{} DH-1(D),H10(D), \nabla(u{\(\rho\)}- u) {\(\cdot\)} (A{\(\rho\)}\nablau{\(\rho\)}- A\nablau) = f{\(\rho\)}- f, u{\(\rho\)}- u where {\(\cdot\)}, {\(\cdot\)}H-1(D),H10(D)denotes the duality product between H-1(D) and H10(D), that we rewrite in the form \hat{}\hat{} DDH-1(D),H10(D). \nabla(u{\(\rho\)}- u) {\(\cdot\)} A{\(\rho\)}\nabla(u{\(\rho\)}- u) = -\nabla(u{\(\rho\)}- u) {\(\cdot\)} (A{\(\rho\)}- A)\nablau + f{\(\rho\)}- f, u{\(\rho\)}- u Using the uniform coercivity of A{\(\rho\)}\in M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)}(D) and Cauchy-Schwarz inequality, this turns into {\(\alpha\)} \nabla(u{\(\rho\)}- u)2L2(D)\leq\nablau {\(\cdot\)} (A{\(\rho\)}- A)\nablau\nabla(u{\(\rho\)}- u) {\(\cdot\)} (A{\(\rho\)}- A)\nabla(u{\(\rho\)}- u) \hat{}1/2\hat{}1/2 DD + f{\(\rho\)}- fH-1(D)u{\(\rho\)}- uH10(D). The first factor of the first term of the r. h. s. vanishes as {\(\rho\)} \to 0 by the Lebesgue dominated convergence2/{\(\alpha\)})|\nablau|2. The first factor of the second theorem since A{\(\rho\)}\to A pointwise and 0 \leq \nablau {\(\cdot\)} (A{\(\rho\)}- A)\nablau \leq ({\(\beta\)} + {\(\beta\)} term of the r. h. s. vanishes as {\(\rho\)} \to 0 as well. Since the other terms are bounded using an a priori estimate and Poincar\'{}e’s inequality, the claim follows. 2.2. Numerical corrector Definition 4.Let H > 0, IH\in N, and let {QH,i}i\in[[1,IH]]be a partition of D in disjoint subdomains of diameter2(D) associated with QH,i: for of order H. We define a family (MH) of approximations of the identity on L every w \in L2(D) and H > 0, MH(w) =w 1QH,i. IH i=1QH,i With the notation of Corollary 1, we define the numerical corrector {\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)}associated with u{\(\rho\)},{\(\epsilon\)}on QH,ias the unique weak solution in H10(QH,i) to -\nabla {\(\cdot\)} A{\(\epsilon\)}MH(\nablau{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)}= 0,(2.14) we set \nablauH,i{\(\rho\)},{\(\epsilon\)}:= MH(\nablau{\(\rho\)},{\(\epsilon\)})|QH,i+ \nabla{\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)} for all 1 \leq i \leq IH, and define the corrector as CH{\(\rho\)},{\(\epsilon\)}=\nablauH,i{\(\rho\)},{\(\epsilon\)}1QH,i. IH i=1 We then have the following corrector result: Theorem 3.Under the assumptions of Corollary 1, the corrector of Definition 4 satisfies lim{\(\rho\)},H\to0lim{\(\epsilon\)}\to0\nablau{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}Lp(D)= 0,(2.15) for all exponents p such that ●1 \leq p \leq 2 if A{\(\epsilon\)}is a family of symmetric matrices, ●1 \leq p < 2 if A{\(\epsilon\)}is not a family of symmetric matrices. In addition, if the r. h. s. f \in H-1(D) of equation (1.5) belongs to W-1,q(D) for some q > 2, then one can take p = 2 in (2.15) even if A{\(\epsilon\)}is not symmetric. Remark 2.The order of the limits in H and {\(\rho\)} in (2.15) is not important, and we may take, e.g., H = {\(\rho\)} \to 0. However, we have to first let {\(\epsilon\)} go to zero. The proof of Theorem 3 is rather long and technical. The main idea is to use Tartar’s correctors of on each element QH,iof the partition of D, pass to the limit in {\(\epsilon\)} first, and then in H. Yet this would require us to know \nablauhoma priori – which we don’t. Hence one has to approximate Tartar’s correctors themselves using \nablau{\(\rho\)},{\(\epsilon\)} in place of \nablauhom. As we shall see, the proof relies on two main arguments: ●the div-curl lemma to prove the convergence of Tartar’s correctors (and of their variants), ●the convergence \nablau{\(\rho\)},{\(\epsilon\)}\to \nablauhomin L2(D, Rd) to prove that the approximations of Tartar’s correctors do not spoil the corrector result. Proof. We recall that u{\(\epsilon\)}, uhom, u{\(\rho\)},{\(\epsilon\)}, and u{\(\rho\)},homare solutions in H10(D) to -\nabla {\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}= f,(2.16) -\nabla {\(\cdot\)} Ahom\nablauhom= f,(2.17) -\nabla {\(\cdot\)} A{\(\rho\)},{\(\epsilon\)}\nablau{\(\rho\)},{\(\epsilon\)}= f, -\nabla {\(\cdot\)} A{\(\rho\)},hom\nablau{\(\rho\)},hom= f, that for all 1 \leq i \leq IH, {\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)}is solution in H10(QH,i) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(MH(\nablau{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)}) = 0,(2.18) and that \nablauH,i{\(\rho\)},{\(\epsilon\)}= MH(\nablau{\(\rho\)},{\(\epsilon\)})|QH,i+ \nabla{\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)}. Likewise, for all 1 \leq i \leq IHwe introduce {\(\gamma\)}H,ihom, {\(\gamma\)}H,i{\(\epsilon\)}, and {\(\gamma\)}H,i{\(\rho\)},hom solutions in H10(QH,i) to -\nabla {\(\cdot\)} Ahom(MH(\nablauhom) + \nabla{\(\gamma\)}H,ihom)= 0,(2.19) -\nabla {\(\cdot\)} A{\(\epsilon\)}(MH(\nablau{\(\epsilon\)}) + \nabla{\(\gamma\)}H,i{\(\epsilon\)})= 0,(2.20) -\nabla {\(\cdot\)} Ahom(MH(\nablau{\(\rho\)},hom) + \nabla{\(\gamma\)}H,i{\(\rho\)},hom)= 0;(2.21) and we set \nablauH,ihom=MH(\nablauhom)|QH,i+ \nabla{\(\gamma\)}H,ihom,(2.22) \nablauH,i{\(\epsilon\)}=MH(\nablau{\(\epsilon\)})|QH,i+ \nabla{\(\gamma\)}H,i{\(\epsilon\)},(2.23) \nablauH,i{\(\rho\)},hom=MH(\nablau{\(\rho\)},hom)|QH,i+ \nabla{\(\gamma\)}H,i{\(\rho\)},hom.(2.24) ESAIM: PROCEEDINGS63 We finally define variants of Tartar’s corrector: CH{\(\rho\)},{\(\epsilon\)}=\nablauH,i{\(\rho\)},{\(\epsilon\)}1QH,i, IH i=1 CHhom=\nablauH,ihom1QH,i, IH i=1 CH{\(\epsilon\)}=\nablauH,i{\(\epsilon\)}1QH,i, IH i=1 CH{\(\rho\)},hom=\nablauH,i{\(\rho\)},hom1QH,i. IH i=1 We have by the triangle inequality \hat{}\hat{}\hat{} DDD |\nablau{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|p|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|p+|CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|p, and we shall show that limH\to0lim sup{\(\epsilon\)}\to0D|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|p=0,(2.25) \hat{} lim{\(\rho\)}\to0lim sup{\(\epsilon\)}\to0D|CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|p=0,uniformly in H,(2.26) \hat{} for all 1 \leq p < 2. Step 1. Proof of (2.25). Using the uniform ellipticity of A{\(\epsilon\)}we write {\(\alpha\)}|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|p\leq(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})p/2.(2.27) \hat{}\hat{} DD By H-convergence we know that A{\(\epsilon\)}(\nablau{\(\epsilon\)}-CH{\(\epsilon\)}) and (\nablau{\(\epsilon\)}-CH{\(\epsilon\)}) converge weakly in L2(D, Rd) to Ahom(\nablauhom-H) and (\nablauhom- CHhom), respectively. This does not imply that the limit of the product converges to the Chom product of the limits, and we need to appeal to compensated compactness. We’d like to pass to the limit {\(\epsilon\)} \to 0 in this estimate for p = 2. Unfortunately, in general, the integrand only converges in the sense of distributions, not pointwise – so that one cannot take the characteristic function of D as a test function. Yet the result will hold true for any 1 \leq p < 2 – the proof of which is slightly technical. We let \phi \in C\infty0(D, [0, 1]) be a (non-negative) function such that \phi \in C\infty0(QH,i) for all 1 \leq i \leq IH, and set D\phi:= {x \in D : \phi(x) = 1}. We have by definition of \phi and by H\"{}older’s inequality with exponents (2/p, 2/(2 - p)) for p < 2 {\(\alpha\)}|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|p \hat{} D \leq(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})\phip/2+\hat{}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})p/2 \hat{} DD\phi \leq(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})\phi|D|(2-p)/2 \hat{}p/2 D +(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})|D\phi|(2-p)/2.(2.28) \hat{}p/2 D\phi We begin with the second term of the r. h. s. which we control by \hat{}p/2p+ CH{\(\epsilon\)}pL2(D))|D\phi|(2-p)/2. D\phi{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})|D\phi|(2-p)/2( \nablau{\(\epsilon\)}L2(D) (\nablau{\(\epsilon\)}- CH An a priori estimate combined with Poincar\'{}e’s inequality on H10(D) yields \nablau{\(\epsilon\)} L2(D)fH-1(D). Likewise, for all 1 \leq i \leq IH, by (2.20) & (2.23), \nablauH,i{\(\epsilon\)}2L2(QH,i)QH,i\nablau{\(\epsilon\)}|QH,i| \leqQH,i|\nablau{\(\epsilon\)}|2, 2\hat{} so that CH{\(\epsilon\)}2L2(D)=\nablauH,i{\(\epsilon\)}1QH,i2L2(D)QH,i|\nablau{\(\epsilon\)}|2= \nablau{\(\epsilon\)}2L2(D)\leq f2H-1(D).(2.29) IHIH\hat{} i=1i=1 Hence, \hat{}p/2p|D\phi|(2-p)/2.(2.30) D\phi{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})|D\phi|(2-p)/2fH-1(D) (\nablau{\(\epsilon\)}- CH We now turn to the first term. We apply the div-curl lemma on each QH,i. On the one hand, by Hconvergence, \nablau{\(\epsilon\)}- \nablauH,i{\(\epsilon\)}is curl free and converges weakly in L2(QH,i, Rd) to \nablauhom- \nablauH,ihom. On the otherH,i), and hand, by H-convergence, A{\(\epsilon\)}(\nablau{\(\epsilon\)}- \nablauH,i{\(\epsilon\)}) converges weakly in L2(QH,i, Rd) to Ahom(\nablauhom- \nablauhom-1(QH,i). Hence, the product (\nablau{\(\epsilon\)}- \nablauH,i{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- \nablauH,i{\(\epsilon\)}) its divergence is bounded by 2 fH-1(D)in H converges to (\nablauhom- \nablauH,ihom) {\(\cdot\)} Ahom(\nablauhom- \nablauH,ihom) in the sense of distributions on QH,i. This implies in particular by definition of \phi that \hat{}{\(\epsilon\)}\to0\to\hat{}(\nablauhom- \nablauH,ihom) {\(\cdot\)} Ahom(\nablauhom- \nablauH,ihom)\phi. QH,i{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- \nablauH,i{\(\epsilon\)})\phiQH,i (\nablau{\(\epsilon\)}- \nablauH,i Since the integrand (\nablauhom- \nablauH,ihom) {\(\cdot\)} Ahom(\nablauhom- \nablauH,ihom) is non-negative and \phi \leq 1, \hat{}\hat{} (\nablauhom- \nablauH,i ESAIM: PROCEEDINGS65 so that by (2.28), (2.30), and the definition of CHhom, lim sup|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|pfpH-1(D)|D\phi|(2-p)/2+D(\nablauhom- CHhom) {\(\cdot\)} Ahom(\nablauhom- CHhom). {\(\epsilon\)}\to0D \hat{}\hat{}p/2 Since \phi is arbitrary, |D\phi| can be chosen as small as desired, and the above estimate turns into lim sup|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|pD(\nablauhom- CHhom) {\(\cdot\)} Ahom(\nablauhom- CHhom).(2.31) {\(\epsilon\)}\to0D \hat{}\hat{}p/2 It remains to prove that the r. h. s. goes to zero as H vanishes. To this aim we use the approximation MHof identity. In particular, by the triangle inequality \hat{} D (\nablauhom- CHhom) {\(\cdot\)} Ahom(\nablauhom- CHhom) =(\nablauhom- MH(\nablauhom) + MH(\nablauhom) - CHhom) {\(\cdot\)} Ahom(\nablauhom- MH(\nablauhom) + MH(\nablauhom) - CHhom) \hat{} D \leq\nablauhom- MH(\nablauhom) {\(\cdot\)} Ahom\nablauhom- MH(\nablauhom) \hat{} D +MH(\nablauhom) - CHhom{\(\cdot\)} AhomMH(\nablauhom) - CHhom(2.32) \hat{} D +\nablauhom- MH(\nablauhom) {\(\cdot\)} AhomMH(\nablauhom) - CHhom \hat{} D +MH(\nablauhom) - CHhom{\(\cdot\)} Ahom\nablauhom- MH(\nablauhom) . \hat{} D On the one hand, by the triangle inequality, MHis a contraction on L2(D), so that MH(\nablauhom)L2(D)\leq \nablauhom L2(D)fH-1(D), and on the other hand it follows from (2.19) & (2.22) (the proof is similar to (2.29)) that CHhom L2(D)fH-1(D).(2.33) Hence, the first, third, and last terms of (2.32) vanish as H \to 0. We therefore focus on the second term: \hat{} D MH(\nablauhom) - CHhom{\(\cdot\)} AhomMH(\nablauhom) - CHhom =(MH(\nablauhom) - \nablauH,ihom) {\(\cdot\)} Ahom(MH(\nablauhom) - \nablauH,ihom). IH\hat{} i=1QH,i By (2.19) & (2.22), \hat{} (MH(\nablauhom) - \nablauH,i for all 1 \leq i \leq IHso that \hat{}H,i) {\(\cdot\)} Ahom(MH(\nablauhom) - \nablauH,ihom) QH,ihom (MH(\nablauhom) - \nablau =(MH(\nablauhom) - \nablauH,ihom) {\(\cdot\)} AhomMH(\nablauhom) \hat{} QH,i =MH(\nablauhom) {\(\cdot\)} AhomMH(\nablauhom) -\nablauH,ihom{\(\cdot\)}QH,iAhomMH(\nablauhom) \hat{}\hat{} QH,iQH,i +\nablauH,ihom{\(\cdot\)}QH,iAhom- AhomMH(\nablauhom) \hat{} QH,i =\nablauH,ihom{\(\cdot\)}QH,iAhom- AhomMH(\nablauhom), \hat{} QH,i using thatfflQH,i\nablauH,ihom= MH(\nablauhom)|QH,i. Hence we need to prove that limH\to0DCHhom{\(\cdot\)} MH(Ahom) - AhomMH(\nablauhom) = 0.(2.34) \hat{} Since MHconverges to Id in L2(D), the second factor of the integrand converges to zero in L2(D). The result essentially follows from the Lebesgue dominated convergence theorem, although one needs to take care of the first factor of the integrand which depends on H as well. We conclude as follows. Let {uhom,{\(\lambda\)}}{\(\lambda\)}be a1(D) as {\(\lambda\)} goes to infinity. By definition of MH, sequence of {\(\lambda\)}-Lipschitz functions which converges to uhomin H MH(\nablauhom,{\(\lambda\)}) is essentially bounded by {\(\lambda\)} for all H > 0. We then write \hat{} D CHhom{\(\cdot\)} MH(Ahom) - AhomMH(\nablauhom) =CHhom{\(\cdot\)} MH(Ahom) - AhomMH(\nablauhom,{\(\lambda\)}) \hat{} D +CHhom{\(\cdot\)} MH(Ahom) - AhomMH(\nablauhom- \nablauhom,{\(\lambda\)}). \hat{} D The first term of the r. h. s.converges to zero as H vanishes by the Cauchy-Schwarz inequality since |MH(\nablauhom,{\(\lambda\)})| \leq {\(\lambda\)}. The second term of the r. h. s. is bounded by a constant times fH-1(D)\nablauhom- \nablauhom,{\(\lambda\)} L2(D)using (2.33) and that MHis a contraction on L2(D). Hence it converges to zero uniformly in H as {\(\lambda\)} \to \infty. We have thus proved (2.34), and therefore using (2.32) that limH\to0D(\nablauhom- CHhom) {\(\cdot\)} Ahom(\nablauhom- CHhom) = 0, \hat{} and finally using (2.31) that limH\to0lim sup{\(\epsilon\)}\to0D|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|p= 0, \hat{} as desired. Step 2. Proof of (2.26). By the uniform ellipticity of A{\(\epsilon\)}, {\(\alpha\)}|CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|p\leq(CH{\(\rho\)},{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(CH{\(\rho\)},{\(\epsilon\)}- CH{\(\epsilon\)})p/2.(2.35) \hat{}\hat{} ESAIM: PROCEEDINGS67 The same string of arguments leading to (2.31) in Step 1 (using compensated compactness) allows to pass to the limsup in {\(\epsilon\)} in (2.35), and yields lim sup|CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|pD(CH{\(\rho\)},hom- CHhom) {\(\cdot\)} Ahom(CH{\(\rho\)},hom- CHhom).(2.36) {\(\epsilon\)}\to0D \hat{}\hat{}p/2 Using then equations (2.19) & (2.22) and (2.21) & (2.24) on each QH,iwhich yield \hat{} QH,i (CH{\(\rho\)},hom- CHhom) {\(\cdot\)} Ahom(CH{\(\rho\)},hom- CHhom) =(\nablauH,i{\(\rho\)},hom- \nablauH,ihom) {\(\cdot\)} Ahom(\nablauH,i{\(\rho\)},hom- \nablauH,ihom) \hat{} QH,i =(MH(\nablauhom) - MH(\nablau{\(\rho\)},hom)) {\(\cdot\)} Ahom(\nablauH,i{\(\rho\)},hom- \nablauH,ihom) \hat{} QH,i =(MH(\nablauhom) - MH(\nablau{\(\rho\)},hom)) {\(\cdot\)} Ahom(CH{\(\rho\)},hom- CHhom), \hat{} QH,i and using the a priori estimates CH{\(\rho\)}L2(D), CH{\(\rho\)},hom L2(D)fH-1(D) (whose proofs are similar to (2.29)), (2.36) turns into lim sup|CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|pD(MH(\nablauhom) - MH(\nablau{\(\rho\)},hom)) {\(\cdot\)} Ahom(CHhom- CH{\(\rho\)},hom) {\(\epsilon\)}\to0D \hat{}\hat{}p/2 MH(\nablauhom- \nablau{\(\rho\)},hom)p/2L2(D)( CH{\(\rho\)}p/2L2(D)+ CH{\(\rho\)},homp/2L2(D)) \nablauhom- \nablau{\(\rho\)},homp/2L2(D)fp/2H-1(D), which converges to zero as {\(\rho\)} \to 0 by Corollary 1. Step 3. Extensions. The starting point is as in Step 1: \hat{}\hat{} D{\(\rho\)},{\(\epsilon\)}|2D(\nablau{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}) |\nablau{\(\epsilon\)}- CH \leq 2(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) + 2(CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}). \hat{}\hat{} DD By symmetry we then have \hat{}\hat{}IH\hat{}IH\hat{} DDi=1Di=1QH,i (\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) =\nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}+\nablauH,i{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablauH,i{\(\epsilon\)}- 2\nablauH,i{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}. (2.37) Extending uH,i{\(\epsilon\)}by zero on D QH,i, we obtain a function in H10(D), and the last term of this identity turns into IH\hat{}IH \nablauH,i{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}=f, uH,i{\(\epsilon\)} using the weak form of the defining equation (2.16). Since uH,i{\(\epsilon\)}converges weakly to uH,ihomin H10(QH,i) (and10(D)), we obtain therefore in H lim{\(\epsilon\)}\to0i=1QH,i\nablauH,i{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}=i=1f, uH,ihomH-1(D),H10(D), IH\hat{}IH and therefore lim{\(\epsilon\)}\to0i=1QH,i\nablauH,i{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}=i=1QH,i\nablauH,ihom{\(\cdot\)} Ahom\nablauhom, IH\hat{}IH\hat{} using (2.17). On the other hand, H-convergence implies the convergence of the energy, so that lim{\(\epsilon\)}\to0D\nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}=D\nablauhom{\(\cdot\)} Ahom\nablauhom, \hat{}\hat{} and for all 1 \leq i \leq IH, lim{\(\epsilon\)}\to0QH,i\nablauH,i{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablauH,i{\(\epsilon\)}=QH,i\nablauH,ihom{\(\cdot\)} Ahom\nablauH,ihom \hat{}\hat{} (the proof of which is as above, starting from (2.16) & (2.17) and (2.20) & (2.19)). Putting things back together yields the desired identity lim sup|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|2lim sup{\(\epsilon\)}\to0D(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\epsilon\)}\to0D \hat{}\hat{} =(\nablauhom- CHhom) {\(\cdot\)} Ahom(\nablauhom- CHhom). \hat{} D Likewise, lim sup|CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}|2lim sup{\(\epsilon\)}\to0D(CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(CH{\(\epsilon\)}- CH{\(\rho\)},{\(\epsilon\)}) {\(\epsilon\)}\to0QH,i \hat{}\hat{} =(CHhom- CH{\(\rho\)},hom) {\(\cdot\)} Ahom(CHhom- CH{\(\rho\)},hom). \hat{} D We then finish the proof as in Step 2, which yields the result for p = 2. When the matrix A{\(\epsilon\)}is not symmetric and the r. h. s. f is in W-1,q(D) for some q > 2, we appeal to Meyers’ estimates [49], which yield the higher integrability result for all {\(\rho\)}, {\(\epsilon\)} > 0: \nablauhom Lq(D), \nablau{\(\epsilon\)} Lq(D), \nablau{\(\rho\)},{\(\epsilon\)} Lq(D), \nablau{\(\rho\)},hom Lq(D)fW-1,q(D), provided q - 2 is small enough (this exponent only depends on the ellipticity constants {\(\alpha\)} and {\(\beta\)}). This allows us to take p = 2 in (2.27) and replace (2.28) by {\(\alpha\)}|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|2 \hat{} D \leq(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})\phi +(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) \hat{}\hat{} DD\phi \leq(\nablau{\(\epsilon\)}- CH{\(\epsilon\)}) {\(\cdot\)} A{\(\epsilon\)}(\nablau{\(\epsilon\)}- CH{\(\epsilon\)})\phi +D|\nablau{\(\epsilon\)}- CH{\(\epsilon\)}|q|D\phi|(q-2)/q, \hat{}\hat{}2/q ESAIM: PROCEEDINGS69 the last term of which can be controlled using Meyers’ estimate on the corrector as well. We then conclude the proof as in Step 2, which yields the desired result for p = 2 in the non-symmetric case. Remark 3.There is some flexibility for the choice of the boundary conditions in the definition of the correctors. In particular, the conclusions of Theorem 3 still hold if the homogeneous Dirichlet boundary conditions in (2.14) are replaced either by homogeneous Neumann boundary conditions, or even a zero average condition (that is ffl\nabla{\(\gamma\)}H,i{\(\rho\)},{\(\epsilon\)}= 0). QH,i The proof for nonsymmetric A{\(\epsilon\)}remains essentially unchanged. Only the proof for symmetric A{\(\epsilon\)}has to be slightly adapted to treat the other boundary conditions. In particular, the right way to write the double product in (2.37) is, for Neumann boundary conditions, \hat{}\hat{}\hat{} QH,i{\(\epsilon\)}=QH,i\nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}MH(\nablau{\(\epsilon\)}) =QH,i\nablau{\(\epsilon\)}{\(\cdot\)}QH,iA{\(\epsilon\)}\nablau{\(\epsilon\)} \nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablauH,i {\(\epsilon\)}\to0-\to\nablauhom{\(\cdot\)}\hat{}Ahom\nablauhom=\hat{}\nablauhom{\(\cdot\)} AhomMH(\nablauhom) QH,iQH,iQH,i =\nablauhom{\(\cdot\)} Ahom\nablauH,ihom, \hat{} QH,i and, for the zero average condition, \hat{}\hat{}\hat{} QH,i{\(\epsilon\)}=QH,iQH,i\nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablauH,i{\(\epsilon\)}+QH,i\nablau{\(\epsilon\)}-QH,i\nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablauH,i{\(\epsilon\)} \nablau{\(\epsilon\)}{\(\cdot\)} A{\(\epsilon\)}\nablauH,i =\nablau{\(\epsilon\)}{\(\cdot\)}A{\(\epsilon\)}\nablauH,i{\(\epsilon\)} \hat{} QH,iQH,i {\(\epsilon\)}\to0-\to\nablauhom{\(\cdot\)}\hat{}Ahom\nablauhom=\hat{}\nablauhom{\(\cdot\)} AhomMH(\nablauhom) QH,iQH,iQH,i =\nablauhom{\(\cdot\)} Ahom\nablauH,ihom. \hat{} QH,i 2.3. Direct approach We are now in position to introduce rigorously the direct approach, which consists in approximating u{\(\rho\)},{\(\epsilon\)} for some {\(\rho\)} > {\(\epsilon\)} on the one hand, and then construct a numerical corrector on the other hand. We present the method for a Galerkin approach. Let {VH} be a suitable sequence of finite-dimensional subspaces of H10(D).H{\(\rho\)},{\(\epsilon\)}\in VHthe unique weak solution in VHto We then denote by u -\nabla {\(\cdot\)} A{\(\rho\)},{\(\epsilon\)}\nablauH{\(\rho\)},{\(\epsilon\)}= f, where A{\(\rho\)},{\(\epsilon\)}is defined by (2.1). From Theorem 2, Corollary 1 and standard approximation arguments, we deduce lim suplim{\(\epsilon\)}\to0uH{\(\rho\)},{\(\epsilon\)}- uhom H1(D)= 0,(2.38) H,{\(\rho\)}\to0 where the limits in H and {\(\rho\)} commute. In practice, the matrix A{\(\rho\)},{\(\epsilon\)}is itself approximated. In particular, for all x \in D and h > 0, denoting by Vh(x) a finite dimensional subspace of H10(Q{\(\rho\)}\cap T-xD), we define an approximation Ah{\(\rho\)},{\(\epsilon\)}of A{\(\rho\)},{\(\epsilon\)}at point x by: for all i, j \in {1, . . . , d}
[9] ij:=ej{\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\rho\)},{\(\epsilon\)},hi(x, y))dy, where v{\(\rho\)},{\(\epsilon\)},hi(x, {\(\cdot\)}) is the unique weak solution in Vh(x) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\rho\)},{\(\epsilon\)},hi(x, y)) = 0in Q{\(\rho\)}\cap T-xD. We then define for all {\(\epsilon\)} > 0, {\(\rho\)} > {\(\epsilon\)}, H > 0 and h < {\(\rho\)} the weak solution uH,h{\(\rho\)},{\(\epsilon\)}in VHto -\nabla {\(\cdot\)} Ah{\(\rho\)},{\(\epsilon\)}\nablauH,h{\(\rho\)},{\(\epsilon\)}= f. A further approximation argument then yields limH,{\(\rho\)}\to0lim{\(\epsilon\)}\to0limh\to0uH,h{\(\rho\)},{\(\epsilon\)}- uhom H1(D)= 0.(2.39) Note that the practical implementation of the method makes use of a quadrature rule on D so that Ah{\(\rho\)},{\(\epsilon\)}only has to be calculated at the quadrature points of D. Let us give the argument for (2.39). The first two limits yield by convergence of the Galerkin method and H-convergence lim{\(\epsilon\)}\to0limh\to0uH,h{\(\rho\)},{\(\epsilon\)}- uhom H1(D)= uH{\(\rho\)},hom- uhom H1(D). By the triangle inequality uH{\(\rho\)},hom- uhom H1(D)\leq uH{\(\rho\)},hom- u{\(\rho\)},hom H1(D)+ u{\(\rho\)},hom- uhom H1(D). The first term of the r. h. s. goes to zero as H \to 0 by convergence of the Galerkin approximation. We need to understand how the convergence depends on {\(\rho\)}. From C\'{}ea’s lemma and Poincar\'{}e’s inequality, we have uH{\(\rho\)},hom- u{\(\rho\)},hom H1(D)infvH\inVHvH- u{\(\rho\)},hom H1(D), which, using the triangle inequality, turns into uH{\(\rho\)},hom- u{\(\rho\)},hom H1(D)u{\(\rho\)},hom- uhom H1(D)+ infvH\inVHvH- uhom H1(D). We thus have uH{\(\rho\)},hom- uhom H1(D)u{\(\rho\)},hom- uhom H1(D)+ infvH\inVHvH- uhom H1(D), which vanishes as {\(\rho\)} and H go to zero (independently, as desired). This shows (2.39). We may then turn to the numerical corrector result. As in Definition 4 we let IH\in N, and {QH,i}i\in[[1,IH]]be a partition of D in disjoint subdomains of diameter of order H. For all h > 0 and i \in [[1, IH]] we let VH,i,hbe a1(QH,i). We define the numerical correctors {\(\gamma\)}H,h,i{\(\rho\)},{\(\epsilon\)}associated with uH,h{\(\rho\)},{\(\epsilon\)}on QH,ias the Galerkin subspace of H0 unique weak solution in VH,i,hto -\nabla {\(\cdot\)} A{\(\epsilon\)}MH(\nablauH,h{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}H,h,i{\(\rho\)},{\(\epsilon\)}= 0, we set \nablavH,h,i{\(\rho\)},{\(\epsilon\)}:= MH(\nablauH,h{\(\rho\)},{\(\epsilon\)})|QH,i+ \nabla{\(\gamma\)}H,h,i{\(\rho\)},{\(\epsilon\)} for all 1 \leq i \leq IH, and finally define CH,h{\(\rho\)},{\(\epsilon\)}=\nablavH,h,i{\(\rho\)},{\(\epsilon\)}1QH,i. IH ESAIM: PROCEEDINGS71 We then have the following numerical corrector result: lim{\(\rho\)},H\to0lim sup{\(\epsilon\)}\to0limh\to0\nablau{\(\epsilon\)}- CH,h{\(\rho\)},{\(\epsilon\)}Lp(D)= 0,(2.40) for all exponents p such that ●1 \leq p \leq 2 if A{\(\epsilon\)}is a family of symmetric matrices, ●1 \leq p < 2 if A{\(\epsilon\)}is not a family of symmetric matrices. In addition, if the r. h. s. f \in H-1(D) of equation (1.5) belongs to W-1,q(D) for some q > 2, then one can take p = 2 in (2.40) even if A{\(\epsilon\)}is not symmetric. This result is essentially Theorem 3, although there is an additional approximation argument needed since MH(\nablauH{\(\rho\)},{\(\epsilon\)}) = MH(\nablau{\(\rho\)},{\(\epsilon\)}). It is enough to note that limH,{\(\rho\)}\to0lim sup{\(\epsilon\)}\to0MH(\nablauH{\(\rho\)},{\(\epsilon\)}) - MH(\nablau{\(\rho\)},{\(\epsilon\)})L2(D)\leq lim supH,{\(\rho\)}\to0lim{\(\epsilon\)}\to0\nablauH{\(\rho\)},{\(\epsilon\)}- \nablau{\(\rho\)},{\(\epsilon\)} L2(D)= 0 by (2.38) & (2.39) to conclude. In this subsection we have proved the convergence of the direct approach to numerical homogenization in the framework of H-convergence. This provides a convergence analysis for the so-called Heterogeneous Multiscale Method (HMM) applied to homogenization problems, as introduced by E et. al. in [17, 18]. It also makes rigorous the numerical corrector approach by Arbogast [6]. 2.4. Dual approach As we have already seen, the dual approach consists in approximating u{\(\epsilon\)}in some adapted Galerkin subspace of H10(D) rather than approximating the H-limit of A{\(\epsilon\)}first. The Multiscale Finite Element (MsFEM) basis is constructed as follows. For all H > 0, let THbe a regular mesh of D by tetrahedra of diameter of order H,10(D). We denote by IHand JHthe number of and let VHbe the associated P 1-finite element subspace of H tetrehedra in THand the dimension of VHrespectively, and we let {{\(\psi\)}H,i}1\leqi\leqJHbe the associated hat functions generating VH. For all {\(\epsilon\)} > 0 and all 0 < h < {\(\epsilon\)} we define multiscale hat functions {{\(\psi\)}H,{\(\epsilon\)},h,i}1\leqi\leqJHby their restrictions on the tetrahedra THof TH. In particular, for every tetrahedron THof TH, we let Vh(TH) be a Galerkin subspace of H10(TH) and let {\(\gamma\)}H,{\(\epsilon\)},h,i|THbe the unique weak solution in Vh(TH) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(\nabla{\(\psi\)}H,i|TH+ \nabla{\(\gamma\)}H,{\(\epsilon\)},h,i) = 0in TH. and set {\(\psi\)}H,{\(\epsilon\)},h,i|TH:= ({\(\psi\)}H,i+ {\(\gamma\)}H,{\(\epsilon\)},h,i)|TH. So defined, the multiscale hat functions {{\(\psi\)}H,{\(\epsilon\)},h,i}1\leqi\leqJHbelong to1(D) and for all i \in {1, . . . , JH}, {\(\psi\)}H,{\(\epsilon\)},h,ihas the same support and the same nodal values as {\(\psi\)}H,i. Hence, the H0 multiscale finite element space VH,{\(\epsilon\)},hspanned by the multiscale hat functions {{\(\psi\)}H,{\(\epsilon\)},h,i}1\leqi\leqJHis a subspace of10(D) of dimension JH. H The approximation uH,{\(\epsilon\)},hof u{\(\epsilon\)}is then defined as the unique weak solution in VH,{\(\epsilon\)},hto -\nabla {\(\cdot\)} A{\(\epsilon\)}\nablauH,{\(\epsilon\)},h= f.(2.41) We then have the following convergence result: limH\to0lim{\(\epsilon\)}\to0limh\to0u{\(\epsilon\)}- uH,{\(\epsilon\)},h W1,p(D)= 0(2.42) for all exponents p such that ●1 \leq p \leq 2 if A{\(\epsilon\)}is a family of symmetric matrices, ●1 \leq p < 2 if A{\(\epsilon\)}is not a family of symmetric matrices. In addition, if the r. h. s. f \in H-1(D) of equation (1.5) belongs to W-1,q(D) for some q > 2, then one can take p = 2 in (2.42) even if A{\(\epsilon\)}is not symmetric. The proof of (2.42) consists in two steps. Let us recall there is a one-to-one mapping MH,{\(\epsilon\)},hMsFEMfrom VHtoJH{\(\nu\)}H,i{\(\psi\)}H,i\in VHwe associate the multiscale finite element function VH,{\(\epsilon\)},h. In particular, with every vH=i=1H,{\(\epsilon\)},h(vH) :=JHi=1{\(\nu\)}H,{\(\epsilon\)},h,i{\(\psi\)}H,{\(\epsilon\)},h,i\in VH,{\(\epsilon\)},h(and vice-versa). We may characterize this mappingj|TH vH,{\(\epsilon\)},h= MMsFEM using corrector fields. In particular, for every tetrahedron THof THand every j \in {1, . . . , d} we let {\(\phi\)}H,{\(\epsilon\)},h be the unique weak solution in Vh(TH) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(ej+ \nabla{\(\phi\)}jH,{\(\epsilon\)},h) = 0in TH, and we set {\(\Phi\)}H,{\(\epsilon\)},h:= ({\(\phi\)}1H,{\(\epsilon\)},h, . . . , {\(\phi\)}dH,{\(\epsilon\)},h). By definition, {\(\Phi\)}H,{\(\epsilon\)},h\in H10(D), and we have for all vH\in VH MH,{\(\epsilon\)},hMsFEM(vH) = vH+ \nablavH{\(\cdot\)} {\(\Phi\)}H,{\(\epsilon\)},h. We denote by uH,{\(\epsilon\)},hthe function of VHassociated with the weak solution uH,{\(\epsilon\)},h\in VH,{\(\epsilon\)},hof (2.41) through the one-to-one mapping uH,{\(\epsilon\)},h= MH,{\(\epsilon\)},hMsFEM-1(uH,{\(\epsilon\)},h). We shall first prove that limH\to0lim{\(\epsilon\)}\to0limh\to0uH,{\(\epsilon\)},h- uhom H1(D)= 0.(2.43) To this aim, we write the weak formulation of (2.41) as follows: for all vH,{\(\epsilon\)},h\in VH,{\(\epsilon\)},h, \hat{} DH,{\(\epsilon\)},h H-1(D),H10(D).(2.44) \nablavH,{\(\epsilon\)},h{\(\cdot\)} A{\(\epsilon\)}\nablauH,{\(\epsilon\)},h= f, v Let us focus on the l. h. s. of (2.44), use the characterization of the mapping MH,{\(\epsilon\)},hMsFEMfrom VHto VH,{\(\epsilon\)},handifor all 1 \leq i \leq IH: that functions of VHare locally affine on TH, and that {\(\Phi\)}H,{\(\epsilon\)},hvanishes on \partialTH \hat{} D \nablavH,{\(\epsilon\)},h{\(\cdot\)} A{\(\epsilon\)}\nablauH,{\(\epsilon\)},h =(\nablavH+ \nablavH{\(\cdot\)} \nabla{\(\Phi\)}H,{\(\epsilon\)},h) {\(\cdot\)} A{\(\epsilon\)}(\nablauH,{\(\epsilon\)},h+ \nablauH,{\(\epsilon\)},h{\(\cdot\)} \nabla{\(\Phi\)}H,{\(\epsilon\)},h) \hat{} D =\nablavH{\(\cdot\)} (Id + \nabla{\(\Phi\)}H,{\(\epsilon\)},h)A{\(\epsilon\)}(Id + \nabla{\(\Phi\)}H,{\(\epsilon\)},h) \nablauH,{\(\epsilon\)},h \hat{} D =|TiH|(\nablavH)|TiHTiH(Id + \nabla{\(\Phi\)}H,{\(\epsilon\)},h)A{\(\epsilon\)}(Id + \nabla{\(\Phi\)}H,{\(\epsilon\)},h) (\nablauH,{\(\epsilon\)},h)|TiH IH i=1 =\nablavH{\(\cdot\)} AH,{\(\epsilon\)},h\nablauH,{\(\epsilon\)},h, \hat{} D where AH,{\(\epsilon\)},his the piecewise constant matrix defined by AH,{\(\epsilon\)},h=(Id + \nabla{\(\Phi\)}H,{\(\epsilon\)},h) {\(\cdot\)} A{\(\epsilon\)}(Id + \nabla{\(\Phi\)}H,{\(\epsilon\)},h)1TiH. IH i=1TiH We then focus on the r. h. s. of (2.44), and assume without loss of generality that f \in L\infty(D) (the general case can be dealt with by approximation), so that f, vH,{\(\epsilon\)},h H-1(D),H10(D)=Df vH,{\(\epsilon\)},h=Df vH+D\nablavH{\(\cdot\)} (f{\(\Phi\)}H,{\(\epsilon\)},h). \hat{}\hat{}\hat{} ESAIM: PROCEEDINGS73 We then define fH,{\(\epsilon\)},has fH,{\(\epsilon\)},h:= f - \nabla {\(\cdot\)} (f{\(\Phi\)}H,{\(\epsilon\)},h), and note that by assumption on f we have fH,{\(\epsilon\)},h\in H-1(D), so that the equation for uH,{\(\epsilon\)},h\in VH,{\(\epsilon\)},hturns into an equation for uH,{\(\epsilon\)},h\in VH: for all vH\in VH, \hat{} DH H-1(D),H10(D). \nablavH{\(\cdot\)} AH,{\(\epsilon\)},h\nablauH,{\(\epsilon\)},h= fH,{\(\epsilon\)},h, v By H-convergence, for all H > 0, the sequence {\(\Phi\)}H,{\(\epsilon\)}:= limh\to0{\(\Phi\)}H,{\(\epsilon\)},hconverges weakly to 0 in H1(D) as {\(\epsilon\)} vanishes, so that for all H > 0 lim{\(\epsilon\)}\to0limh\to0{\(\Phi\)}H,{\(\epsilon\)},h L2(D)= 0.(2.45) We are in position to prove that limH\to0lim{\(\epsilon\)}\to0limh\to0uH,{\(\epsilon\)},h- uhom H1(D)= 0. Let denote by uH,homthe weak solution in VHto -\nabla {\(\cdot\)} Ahom\nablauH,hom= f, and by uH,{\(\epsilon\)},hhomthe weak solution in VHto -\nabla {\(\cdot\)} Ahom\nablauH,{\(\epsilon\)},hhom= fH,{\(\epsilon\)},h, Then, by the triangle inequality \nablauH,{\(\epsilon\)},h- \nablauhom L2(D)\leq \nablauH,{\(\epsilon\)},h- \nablauH,{\(\epsilon\)},hhomL2(D) + \nablauH,hom- \nablauhom L2(D)+ \nablauH,{\(\epsilon\)},hhom- \nablauH,hom L2(D).(2.46) We treat the three terms of the r. h. s. separately and start with the first one. The function uH,{\(\epsilon\)},h- uH,{\(\epsilon\)},hhomisthe weak solution in VHto -\nabla {\(\cdot\)} AH,{\(\epsilon\)},h\nabla(uH,{\(\epsilon\)},h- uH,hom) = -\nabla {\(\cdot\)} (Ahom- AH,{\(\epsilon\)},h)\nablauH,{\(\epsilon\)},hhom, so that \nablauH,{\(\epsilon\)},h- \nablauH,{\(\epsilon\)},hhomL2(D)(Ahom- AH,{\(\epsilon\)},h)\nablauH,{\(\epsilon\)},hhomL2(D) \leq (Ahom- AH,{\(\epsilon\)},h)\nablauhom L2(D)+ 2{\(\beta\)} \nablauhom- \nablauH,hom L2(D)+ 2{\(\beta\)} \nablauH,hom- \nablauH,{\(\epsilon\)},hhomL2(D). The first term of the r. h. s. converges to zero as {\(\epsilon\)} and H go to zero by the dominated convergence theorem and using the fact that the following convergence holds pointwise, as in Subsection 2.1, limH\to0lim{\(\epsilon\)}\to0limh\to0AH,{\(\epsilon\)},h= Ahom. The second term coincides with the second term of the r. h. s. of (2.46), and vanishes as H \to 0 by convergence of the Galerkin method for the homogenized equation. We now treat the last and third term, which coincides with the third term of the r. h. s. of (2.46). We recall that uH,hom- uH,{\(\epsilon\)},hhomis the weak solution in VHto -\nabla {\(\cdot\)} Ahom\nabla(uH,hom- uH,{\(\epsilon\)},hhom) = f - fH,{\(\epsilon\)},h= -\nabla {\(\cdot\)} (f{\(\Phi\)}H,{\(\epsilon\)},h), so that \nablauH,hom- \nablauH,{\(\epsilon\)},hhomL2(D)f {\(\Phi\)}H,{\(\epsilon\)},h L2(D) and therefore using (2.45) and the assumption that f \in L\infty(D), for all H > 0, lim{\(\epsilon\)}\to0limh\to0\nablauH,hom- \nablauH,{\(\epsilon\)},hhomL2(D)= 0. We have thus proved (2.43). It remains to note that \nablauH,{\(\epsilon\)},his the corrector associated with \nablauH,{\(\epsilon\)},hand with the partition THof D. Hence, from (2.43) and the same string of arguments as for the direct approach, we deduce that limH\to0lim{\(\epsilon\)}\to0limh\to0\nablau{\(\epsilon\)}- \nablauH,{\(\epsilon\)},h Lp(D)= 0, for all exponents p such that ●1 \leq p \leq 2 if A{\(\epsilon\)}is a family of symmetric matrices, ●1 \leq p < 2 if A{\(\epsilon\)}is not a family of symmetric matrices. This implies the desired convergence result (2.42) by Poincar\'{}e’s inequality for u{\(\epsilon\)}- uH,{\(\epsilon\)},h\in W1,p0(D): Instead of a Galerkin approximation uH,{\(\epsilon\)},hof u{\(\epsilon\)}, we could have considered a Petrov-Galerkin approximation of u{\(\epsilon\)}(in which case the test-functions are in VH, not in VH,{\(\epsilon\)},h). The convergence proof is indeed simpler (one does not need to introduce fH). This variant will be used in the next section. 3.Resonance, windowing, and oversampling In the previous section we have introduced an analytical framework and proved the convergence of some numerical homogenization methods within the framework of H-convergence. Quantitative convergence rates further depend on the class of heterogeneities considered. In this section, we provide convergence rates for the simplest heterogeneities possible, that is we assume the coefficients A{\(\epsilon\)}to be {\(\epsilon\)}-periodic. This allows us to give a complete numerical analysis of the methods, and identify the limiting term in the error. This term is the so-called resonance error. It is related to the boundary conditions used for the corrector. We shall then recall a standard way to reduce the resonance error (windowing and oversampling), check it does indeed reduce the error in the case of periodic structures, and then adapt the analytical framework of Section 2 to include windowing and oversampling. 3.1. Numerical analysis of the periodic case and the resonance error In this subsection we assume that A{\(\epsilon\)}= A({\(\cdot\)}/{\(\epsilon\)}) where A is a symmetric Q = (-1/2, 1/2)d-periodic matrix. In this case, the homogenized matrix Ahomis symmetric, does not depend on the macroscopic space variable,d, and is charaterized by: for all {\(\xi\)} \in R {\(\xi\)} {\(\cdot\)} Ahom{\(\xi\)} =({\(\xi\)} + \nabla{\(\phi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}), \hat{} Q where {\(\phi\)} \in H1#(Q) is the unique Q-periodic weak solution to the corrector equation -\nabla {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}) = 0. Furthermore, we let f \in L\infty(D) and D be smooth enough so that by elliptic regularity, the solution uhom\in H10(D) to the homogenized problem is indeed of class H2(D). ESAIM: PROCEEDINGS75 Then, it is proved in [1, 17, 18] that for the direct approach we have for {\(\rho\)} \sim H, P 1-finite elements for both the macroscopic and the microscopic variables, uH,hH,{\(\epsilon\)}- uhom H1(D)H +{\(\epsilon\)}H+h{\(\epsilon\)}2,(3.1) \nablau{\(\epsilon\)}-\nablavH,h,iH,{\(\epsilon\)}1QH,iL2(D)H +{\(\epsilon\)}H+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)}.(3.2) IH i=1 For the dual approach, it is shown in [42, 43] that uH,{\(\epsilon\)},h- u{\(\epsilon\)} H1(D)H +{\(\epsilon\)}H+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)}.(3.3) We do not display the proofs of these estimates, and refer the reader to the original papers. Note however that: ●the term O(H) comes from the discretization of D, as it is standard for P 1-finite elements, ●the term O(h/{\(\epsilon\)}) in (3.2) and (3.3) comes from the discretizations of the mesh elements THand QH with a meshsize h such that h \ll {\(\epsilon\)},\sqrt{}{\(\epsilon\)} is a theoretical limit due to boundary layers in the neigborhood of \partialD in periodic homog●the term enization, ●the third term in (3.1) may be surprising at a first glance. This part of the error is indeed driven by the finite element error in the computation of the corrector, which is squared when computing the approximation of the homogenized matrix – due to symmetry. The limiting term is however the term involving{\(\epsilon\)}H, which is the inverse of a measure of the number of {\(\epsilon\)}-rescaledd). The largerH{\(\epsilon\)}, the more expensive the method,so that there is a trade-off cost/accuracy between H and {\(\epsilon\)}. This error is even more important at the level of periodic cells contained in QHor TH(which is of order (H/{\(\epsilon\)}) the numerical corrector since its square-root appears in (3.2) and (3.3). There are two sources of this error ●The homogeneous Dirichlet boundary conditions used in (2.2) are not consistent with the periodic boundary conditions of the corrector {\(\phi\)}, ●The average (2.1) defining A{\(\rho\)},{\(\epsilon\)}is not consistent with the average defining Ahomsince Q{\(\rho\)}\cap T-xD is not a multiple of periodic cells in general. A first idea to reduce the error due to boundary conditions on (2.2) consists in imposing the boundary conditions far from the domain of interest, hoping the error is localized on a neighborhood of the boundary. This is indeed the case, as shown on a half plane by Bensoussan, Lions, and Papanicolaou [7]. In particular, this is efficient at the level of the corrector, but not at the level of the homogenized coefficients. To illustrate this, we display the results of two academic series of tests. In the first series of tests, we compare the homogenized coefficients Ahomto two different approximations: AR, which is defined as {\(\xi\)} {\(\cdot\)} AR{\(\xi\)} :=({\(\xi\)} + \nabla{\(\phi\)}R) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}R) QR where QR= (-R/2, R/2)dfor R \in N, and {\(\phi\)}Ris the unique solution in H10(QR) to -\nabla {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}R) = 0, and \tilde{}ARdefined as {\(\xi\)} {\(\cdot\)} \tilde{}AR{\(\xi\)} :=({\(\xi\)} + \nabla{\(\phi\)}R) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}R). Table 1.Error on the approximated homogenized coefficients (performed with [26, FreeFEM ], smooth periodic coefficients, P 2-finite elements, and 100 elements per periodic cell). Number R of periodic cells|Ahom- AR||Ahom- \tilde{}AR| per dimension ErrorRate ofPrefactorErrorRate ofPrefactor convergence(rate=1)convergence(rate=1) 10.157-0.1570.157-0.157 20.08450.8950.1690.02102.900.0420 40.04330.9630.1730.01180.8350.0471 80.02190.9830.1750.005970.9790.0478 120.01461.010.1750.003971.000.0476 160.01100.9650.1760.002990.9850.0478 200.008761.030.1750.002391.000.0478 Table 2.L2-norm of the error on the corrector (performed with [26, FreeFEM ], smooth periodic coefficients, P 2-finite elements, and 100 elements per periodic cell). Number R of periodic cellsE1(R)E2(R) per dimension ErrorRate ofPrefactorErrorRate ofPrefactor convergence(rate=0.5)convergence(rate=1) 10.210-0.2100.210-0.210 20.1560.4250.2210.01160.8930.0232 40.1130.4680.2260.003611.6840.0144 80.08080.4840.2290.001810.9880.0145 120.06620.4910.2290.001211.000.0145 160.05740.4960.2300.0009100.9920.0146 200.05150.4920.2300.0007261.020.0145 In particular, \tilde{}ARis the best approximation of Ahomone can devise using {\(\phi\)}Rsince the center of the computation domain is the place where the effect of the Dirichlet boundary conditions is expected to be the smallest. As can be seen on Table 1, as expected |Ahom- \tilde{}AR| \leq |Ahom- AR|. Yet the rate of convergence in R is -1 in both cases (which corresponds to the term {\(\epsilon\)}/H in (3.1). In Table 2, we do not compare homogenized coefficients, but rather the correctors themselves, and define two errors: E1(R) :=|\nabla{\(\phi\)}R- \nabla{\(\phi\)}|2, 1/2 QR E2(R) :=|\nabla{\(\phi\)}R- \nabla{\(\phi\)}|2. 1/2 Q This time, not only the prefactor of the error changes, but also the convergence rate which passes from 1/2 to 1, which will improve both (3.2) and (3.3), as we shall see below. In order to reduce the second source of error (the fact that the average (2.1) may be calculated on a domain which is not a multiple of the periodic cell), one may use different averaging functions than simply the indicator function which are such that they approximate the mean of a periodic function at a higher order – we call such functions masks, and this approach filtering, see Definition 7 in Section 4. To fix the vocabulary, ESAIM: PROCEEDINGS77 ●windowing amounts to approximating the corrector on a larger domain than needed to reduce the effect of spurious boundary conditions, ●filtering amounts to approximating averages by integrals with a weighted measure (using the so-called masks). 3.2. Windowing and filtering in the direct approach For the direct approach, the use of windowing is straightforward to describe. Let {\(\delta\)} > 1 be fixed, and {\(\eta\)} : Q \to R+be an integrable function of mass one (that is, an averaging mask on Q). For all r > 0, we define+by {\(\eta\)}r(y) = r-d{\(\eta\)}(y/r), which is an averaging mask on Qr. Then, for all x \in D and h > 0,10(Q{\(\delta\)}{\(\rho\)}\cap T-xD), we define an approximation A{\(\delta\)},h{\(\rho\)},{\(\epsilon\)}of A{\(\rho\)},{\(\epsilon\)} {\(\eta\)}r: Qr\to R denoting by Vh(x) a finite dimensional subspace of H at point x by: for all i, j \in {1, . . . , d}
[10] ij:=Q{\(\rho\)}\capT-xD(ej+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)},hj(x, y))) {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)},hi(x, y)){\(\eta\)}{\(\rho\)}(y)dy, \hat{} where for all k \in {1, . . . , d}, v{\(\delta\)}{\(\rho\)},{\(\epsilon\)},hk(x, {\(\cdot\)}) is the unique weak solution in Vh(x) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(x + y)(ek+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)},hk(x, y)) = 0in Q{\(\delta\)}{\(\rho\)}\cap T-xD. We define for all {\(\epsilon\)} > 0, {\(\rho\)} > {\(\epsilon\)}, H > 0 and h < {\(\rho\)} the weak solution u{\(\delta\)},H,h{\(\rho\)},{\(\epsilon\)}in VHto -\nabla {\(\cdot\)} A{\(\delta\)},h{\(\rho\)},{\(\epsilon\)}\nablau{\(\delta\)},H,h{\(\rho\)},{\(\epsilon\)}= f. We may then turn to the numerical corrector. As in Definition 4 we let IH\in N, and {QH,i}i\in[[1,IH]]be a partition of D into disjoint subdomains of diameter of order H. We further set Q{\(\delta\)}H,i:= {x \in D | d(x, QH,i) < {\(\delta\)}H}, which is an enlarged version of QH,i. For all h > 0 and i \in [[1, IH]] we let V{\(\delta\)}H,i,hbe a Galerkin subspace of H10(Q{\(\delta\)}H,i). We define the numerical correctors {\(\gamma\)}{\(\delta\)},H,h,i{\(\rho\)},{\(\epsilon\)}associated with u{\(\delta\)},H,h{\(\rho\)},{\(\epsilon\)}as the unique weak solution in V{\(\delta\)}H,i,hto -\nabla {\(\cdot\)} A{\(\epsilon\)}MH(\nablau{\(\delta\)},H,h{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}{\(\delta\)},H,h,i{\(\rho\)},{\(\epsilon\)}= 0, we set \nablav{\(\delta\)},H,h,i{\(\rho\)},{\(\epsilon\)}:= MH(\nablau{\(\delta\)},H,h{\(\rho\)},{\(\epsilon\)})|QH,i+ (\nabla{\(\gamma\)}{\(\delta\)},H,h,i{\(\rho\)},{\(\epsilon\)})|QH,i for all 1 \leq i \leq IH, and we define C{\(\delta\)},H{\(\rho\)},{\(\epsilon\)}=\nablav{\(\delta\)},H,h,i{\(\rho\)},{\(\epsilon\)}1QH,i. IH i=1 In the case of periodic coefficients, E, Ming and Zhang [18] have essentially proved that (3.1) and (3.2) are replaced by u{\(\delta\)},H,hH,{\(\epsilon\)}- uhom H1(D)H +{\(\epsilon\)}H+h{\(\epsilon\)}2,(3.4) \nablau{\(\epsilon\)}- C{\(\delta\)},H,hH,{\(\epsilon\)}L2(D)H +{\(\epsilon\)}H+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)},(3.5) where the multiplicative constant depends on {\(\delta\)}. This is an improvement for the reconstruction of the fine scale features since{\(\epsilon\)}Hin (3.2) is replaced by{\(\epsilon\)}Hin (3.5).Yet, as already mentioned and further studied in [59], the estimate (3.4) is not better than (3.1) in terms of convergence rates. In addition, there are examples in [29, 59] (using reasonable averaging masks) for which the prefactor in (3.4) is larger than in (3.1). It is therefore not clear whether the use of a mask yields better results than simply taking the average in general. A notable exception is provided in [10], where the mask is not used only as a post-processing to compute the average, but already in the definition of the bilinear form associated with the weak form of the corrector equation. A formal two-scale expansion shows that using this modified corrector equation allows to replace the error corresponding to the second term in (3.1) and (3.4) by {\(\epsilon\)}2– which is definitely better. We do not give more details on this method since we shall present an even Hmore efficient approach in Section 4. 3.3. Oversampling in the dual approach For the dual approach, windowing is usually called oversampling. We construct a new Multiscale Finite Element (MsFEM) basis as follows. Let {\(\delta\)} > 1 be fixed. For all H > 0, let THbe a regular mesh of D by tetrahedra of diameter of order H, and let VHbe the associated P 1-finite element subspace of H10(D). We denote by IHand JHthe number of tetrehedra in THand the dimension of VHrespectively, and we let {{\(\psi\)}H,i}1\leqi\leqJHbe the associated hat functions generating VH. For all {\(\epsilon\)} > 0 and all 0 < h < {\(\epsilon\)} we define multiscale{\(\delta\)}H,{\(\epsilon\)},h,i1\leqi\leqJHby their restrictions on the tetrahedra (Tk)1\leqk\leqIHof TH. In particular, for every hat functions {{\(\psi\)} tetrahedron Tkof TH, we define T{\(\delta\)}k:= {x \in D | d(x, Tk) < ({\(\delta\)} - 1)H}, denote by Vh(T{\(\delta\)}k) a Galerkin subspace of{\(\delta\)},kbe the unique weak solution in Vh(T{\(\delta\)}k) to H10(T{\(\delta\)}k) and let {\(\gamma\)}H,{\(\epsilon\)},h,i -\nabla {\(\cdot\)} A{\(\epsilon\)}(\nabla{\(\psi\)}H,i|Tk+ \nabla{\(\gamma\)}{\(\delta\)},kH,{\(\epsilon\)},h,i) = 0in T{\(\delta\)}k, and set {\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i|Tk:= ({\(\psi\)}H,i+ {\(\gamma\)}{\(\delta\)},kH,{\(\epsilon\)},h,i)|Tk. So defined, the multiscale hat functions {{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i}1\leqi\leqJHdo not belong to H10(D), and we use the notation \nablaH{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,ifor the square integrable function of L2(D) defined on each element Tkby \nablaH{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i|Tk:= (\nabla{\(\psi\)}H,i+ \nabla{\(\gamma\)}{\(\delta\)},kH,{\(\epsilon\)},h,i)|Tk. Note that \nablaH{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,iis not the distributional derivative of {\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i(it is the absolutely continuous part of this derivative with respect to the Lebesgue measure). Hence, the multiscale finite element space V{\(\delta\)}H,{\(\epsilon\)},hspanned by the multiscale hat functions {{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i}1\leqi\leqJHis not a subspace of H10(D) (note that it has dimension JH). The approximation u{\(\delta\)}H,{\(\epsilon\)},h\in V{\(\delta\)}H,{\(\epsilon\)},hof u{\(\epsilon\)}can be defined using either a (discontinuous) Galerkin, or using a Petrov-(discontinuous) Galerkin method (for which the test-functions are in VH, not in V{\(\delta\)}H,{\(\epsilon\)},h). We focus on the Petrov-(discontinuous) Galerkin method, and define u{\(\delta\)}H,{\(\epsilon\)},h\in V{\(\delta\)}H,{\(\epsilon\)},has the unique solution (this statement has to be proved) to: For all vH\in VH, \hat{} DH H-1(D),H10(D).(3.6) \nablavH{\(\cdot\)} A{\(\epsilon\)}\nablaHu{\(\delta\)}H,{\(\epsilon\)},h= f, v In the case of periodic coefficients, both for the (discontinuous) Galerkin and Petrov-Galerkin method (see Efendiev and Hou [22] and Hou, Wu, and Zhang [44]), (3.3) is replaced by the improved estimate u{\(\delta\)}H,{\(\epsilon\)},h- u{\(\epsilon\)} H1H(D)H +{\(\epsilon\)}H+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)},(3.7) where the{\(\cdot\)}H1H(D)is a notation for the broken norm v2H1H(D)= v2L2(D)+ \nablaHv2L2(D,Rd). In [22], two contributions to the resonance error are identified and are proved to be of the same order. In [44], it is showed that one of these two contributions disappears by using the Petrov-Galerkin formulation. Hence, ESAIM: PROCEEDINGS79 the multiplicative constant in (3.7) is expected to be smaller for the Petrov-Galerkin formulation than for the Galerkin formulation – and this is confirmed by numerical tests (the resonance error may even completely disappear in some specific situations, see [44]). So far we’ve seen that windowing and oversampling are efficient in the periodic case to reduce the resonance error for the corrector. In the following paragraph, we show that all the convergence results proved for general H-convergent coefficients hold as well with windowing and oversampling. 3.4. Analytical framework In this paragraph, we shall generalize Theorems 2 and 3 to the case of windowing. The proofs we present here are variants of the ones of [28] (which treat the case of nonlinear operators). We begin with the definition of a local approximation of Ahomon domains of size {\(\rho\)} > 0. Definition 5.Let {\(\delta\)} > 1, and let {\(\eta\)} : Q{\(\delta\)}\to [0, {\(\delta\)}] be a measurable function of mass one (called a mask), such-d{\(\delta\)}] be the rescaled version {\(\eta\)}{\(\rho\)}: y \to {\(\rho\)}-d{\(\eta\)}(y/{\(\rho\)}) of that infQ{\(\eta\)} \geq 1/{\(\delta\)}, and for all {\(\rho\)} > 0, let {\(\eta\)}{\(\rho\)}: C(x, {\(\delta\)}{\(\rho\)}) \to [0, {\(\rho\)} {\(\eta\)}. For all {\(\rho\)} > 0 and {\(\epsilon\)} > 0, we denote by A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}: D \to Mdthe function defined by: for all i, j \in {1, . . . , d} and for x \in D,}
[11] ij:={\(\eta\)}{\(\rho\)}-1\hat{}(ej+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}j(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}i(x, y)){\(\eta\)}{\(\rho\)}(y)dy, (3.8) \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xDQ{\(\delta\)}{\(\rho\)}\capT-xD where for all k \in {1, . . . , k}, v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}k(x, {\(\cdot\)}) is the unique weak solution in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} A{\(\epsilon\)}(x + y)(ek+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}k(x, y)) = 0in Q{\(\delta\)}{\(\rho\)}\cap T-xD.(3.9) Unless {\(\eta\)} is a constant function, it is not clear a priori whether A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}is a coercive matrix. This is indeed the case under mild conditions, which are stated in the main result of this section: Theorem 4.Let D be smooth. Let A{\(\epsilon\)}be a H-convergent sequence, and A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}be as in Definition 5. Then there exist {\(\delta\)} > 1 small enough and {\(\beta\)}’\geq {\(\alpha\)}’> 0 such that for all {\(\rho\)} > 0 and {\(\epsilon\)} > 0, A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}\in M{\(\alpha\)}’{\(\beta\)}’(D). In addition, for all {\(\rho\)} > 0, there exists A{\(\delta\)}{\(\rho\)},hom\in M{\(\alpha\)}’{\(\beta\)}’(D) such that for almost every x \in D, lim{\(\epsilon\)}\to0A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}(x)= A{\(\delta\)}{\(\rho\)},hom(x),(3.10) lim{\(\rho\)}\to0A{\(\delta\)}{\(\rho\)},hom(x)= Ahom(x).(3.11) As a direct corollary we have Corollary 2.Let A{\(\epsilon\)}, A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}and A{\(\delta\)}{\(\rho\)},hombe as in Theorem 4, and f \in H-1(D). Then, the weak solution u{\(\delta\)}{\(\rho\)},{\(\epsilon\)}\in H10(D) to -\nabla {\(\cdot\)} A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}\nablau{\(\delta\)}{\(\rho\)},{\(\epsilon\)}= f satisfies lim{\(\rho\)}\to0lim{\(\epsilon\)}\to0u{\(\delta\)}{\(\rho\)},{\(\epsilon\)}- uhom H1(D)= 0,(3.12) where uhom\in H10(D) is the weak solution to -\nabla {\(\cdot\)} Ahom\nablauhom= f. As a consequence of H-convergence we also have that lim{\(\rho\)}\to0lim{\(\epsilon\)}\to0u{\(\delta\)}{\(\rho\)},{\(\epsilon\)}- u{\(\epsilon\)} L2(D)= 0, Before we proceed with the proofs, let us give the associated corrector result. Definition 6.Let H > 0, IH\in N, and let {QH,i}i\in[[1,IH]]be a partition of D in disjoint subdomains of diameter2(D) associated with QH,i: for of order H. We define a family (MH) of approximations of the identity on L every w \in L2(D) and H > 0, MH(w) =w 1QH,i. IH i=1QH,i Let {\(\delta\)} > 1 be as in Theorem 4, and for all i \in [[1, IH]], set Q{\(\delta\)}H,i:= {x \in D | d(x, QH,i) < ({\(\delta\)} - 1)H}. With the notation of Corollary 2, we define the numerical correctors {\(\gamma\)}{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)}associated with u{\(\delta\)}{\(\rho\)},{\(\epsilon\)}as the unique weak solution in H10(Q{\(\delta\)}H,i) to -\nabla {\(\cdot\)} A{\(\epsilon\)}MH(\nablau{\(\delta\)}{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)}= 0,(3.13) we set \nablau{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)}:= MH(\nablau{\(\delta\)}{\(\rho\)},{\(\epsilon\)})|QH,i+ (\nabla{\(\gamma\)}{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)})|QH,i for all 1 \leq i \leq IH, and we define C{\(\delta\)},H{\(\rho\)},{\(\epsilon\)}:=\nablau{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)}1QH,i. IH i=1 We then have the following corrector result: Theorem 5.Under the assumptions of Corollary 2, the corrector of Definition 6 satisfies lim{\(\rho\)},H\to0lim{\(\epsilon\)}\to0\nablau{\(\epsilon\)}- C{\(\delta\)},H{\(\rho\)},{\(\epsilon\)}Lp(D)= 0,(3.14) for all exponents p such that ●1 \leq p \leq 2 if A{\(\epsilon\)}is a family of symmetric matrices, ●1 \leq p < 2 if A{\(\epsilon\)}is not a family of symmetric matrices. In addition, if the r. h. s. f \in H-1(D) of equation (1.5) belongs to W-1,q(D) for some q > 2, then one can take p = 2 in (3.14) even if A{\(\epsilon\)}is not symmetric. We begin with the proof of Theorem 4, which has the same structure as the proof of Theorem 2. Proof of Theorem 4. We first prove that A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}\in M{\(\alpha\)}’{\(\beta\)}’(D) for some {\(\delta\)} > 1 small enough. By regularity of the domain D, for all 0 < {\(\gamma\)} < 1, there exists {\(\delta\)} > 1 such that for all x \in D, |Q{\(\rho\)}\cap T-xD|\geq {\(\gamma\)}.(3.15) |Q{\(\delta\)}{\(\rho\)}\cap T-xD| By definition of {\(\eta\)}, non-negativity of the integrand, and using that {\(\eta\)}{\(\rho\)}\leq {\(\rho\)}-d{\(\delta\)} and infQ{\(\eta\)} \geq {\(\rho\)}-d/{\(\delta\)} by assumption, one has for all {\(\xi\)} \in Rdsuch that |{\(\xi\)}| = 1 {\(\xi\)} {\(\cdot\)} A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}={\(\eta\)}{\(\rho\)}-1\hat{}({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xDQ{\(\delta\)}{\(\rho\)}\capT-xD \geq1{\(\delta\)}2|Q{\(\delta\)}{\(\rho\)}\cap T-xD|\hat{}Q{\(\rho\)}\capT-xD({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y))dy, ESAIM: PROCEEDINGS81 where v{\(\delta\)}{\(\rho\)},hom{\(\xi\)}(x, {\(\cdot\)}) is the unique solution in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} Ahom(x + y)({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},hom{\(\xi\)}(x, y)) = 0. By Meyers’ estimate and the regularity of D, there exist q > 2 and {\(\sigma\)} > 0 such that for all x \in D and all{\(\delta\)}{\(\rho\)},hom(x, {\(\cdot\)}) \in W1,q0(Q{\(\delta\)}{\(\rho\)}\cap T-xD) and {\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},hom{\(\xi\)}(x, {\(\cdot\)})pLp(Q{\(\delta\)}{\(\rho\)}\capT-xD)\leq {\(\sigma\)}|Q{\(\delta\)}{\(\rho\)}\cap T-xD|. {\(\xi\)} \in Rdwith |{\(\xi\)}| = 1, v{\(\xi\)} Hence, by H\"{}older’s inequality, the fact that |{\(\xi\)}| = 1, the argument leading to (2.4) (coercivity of A{\(\epsilon\)}and Jensen’s inequality), and (3.15), the estimate above turns into {\(\xi\)} {\(\cdot\)} A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}\geq1{\(\delta\)}2|Q{\(\rho\)}\cap T-xD|\hat{}Q{\(\delta\)}{\(\rho\)}\capT-xD({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y))dy -1{\(\delta\)}2|Q{\(\delta\)}{\(\rho\)}\cap T-xD|\hat{}(Q{\(\delta\)}{\(\rho\)}\Q{\(\rho\)})\capT-xD({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablav{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}(x, y))dy \geq{\(\alpha\)}{\(\delta\)}2-{\(\beta\)}{\(\delta\)}2|Q{\(\delta\)}{\(\rho\)}\cap T-xD||(Q{\(\delta\)}{\(\rho\)}Q{\(\rho\)}) \cap T-xD|(q-2)/q({\(\sigma\)}|Q{\(\delta\)}{\(\rho\)}\cap T-xD|)2/q =1{\(\delta\)}2({\(\alpha\)} - {\(\beta\)}(1 - {\(\gamma\)})(q-2)/q{\(\sigma\)}2/q). Since q and {\(\sigma\)} only depend on {\(\alpha\)}, {\(\beta\)} and D, there exists 0 < {\(\gamma\)} < 1 and therefore some {\(\delta\)} > 1 such that {\(\xi\)} {\(\cdot\)} A{\(\delta\)}{\(\rho\)},{\(\epsilon\)}{\(\xi\)}\geq{\(\alpha\)}2, as desired. The upper bound is proved as in (2.3). Let x \in D and {\(\rho\)} > 0, and consider problem (2.2). By the locality and definition of H-convergence (see property (2) of Lemma 1 and Definition 1), for all k \in {1, . . . , d}, v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}k(x, {\(\cdot\)}) \rightharpoonup v{\(\delta\)}{\(\rho\)},homk(x, {\(\cdot\)}) in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD),{\(\delta\)}{\(\rho\)},{\(\epsilon\)}(x, {\(\cdot\)})) \rightharpoonup Ahom(x + {\(\cdot\)})(ek+ \nablayv{\(\delta\)}{\(\rho\)},homk(x, {\(\cdot\)})) in L2(Q{\(\delta\)}{\(\rho\)}\cap T-xD, Rd),(3.16) A{\(\epsilon\)}(x + {\(\cdot\)})(ek+ \nablayvk where v{\(\delta\)}{\(\rho\)},homk(x, {\(\cdot\)}) is the unique solution in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} Ahom(x + y)(ek+ \nablayv{\(\delta\)}{\(\rho\)},homk(x, y)) = 0. Hence, setting
[12] ij:=Q{\(\delta\)}{\(\rho\)}\capT-xD{\(\eta\)}{\(\rho\)}-1\hat{}Q{\(\delta\)}{\(\rho\)}\capT-xD(ej+\nablayv{\(\delta\)}{\(\rho\)},homj(x, y)){\(\cdot\)}Ahom(x+y)(ei+\nablayv{\(\delta\)}{\(\rho\)},homi(x, y)){\(\eta\)}{\(\rho\)}(y)dy, \hat{} (3.10) follows from (3.16) by the div-curl lemma and Meyers’ estimates. Indeed, ej+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}j(x, y) is a gradient, and therefore curl-free, whereas Ahom(x + y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}i(x, y)) is divergence free by definition, so that their product converges in the sense of distributions to the product of the limits as {\(\epsilon\)} \to 0. Yet, by Meyers’ estimates,{\(\delta\)}{\(\rho\)},{\(\epsilon\)}(x, y) is in Lq(Q{\(\delta\)}{\(\rho\)}\cap T-xD) so that y \to (ej+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}j(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}i(x, y)){\(\eta\)}{\(\rho\)}(y) is y \to \nablayvk bounded in Lq/2(Q{\(\delta\)}{\(\rho\)}\cap T-xD) with q/2 > 1, which upgrades the convergence in the sense of distributions to a weak convergence in Lq/2(Q{\(\delta\)}{\(\rho\)}\cap T-xD) (by the Banach-Alaoglu theorem), and yields the desired result (3.10). To prove (3.11), we appeal to Lemma 3. To this aim, we note that for all {\(\rho\)} small enough, Q{\(\delta\)}{\(\rho\)}\subset T-xD, so that\hat{}{\(\delta\)}{\(\rho\)},hom(x, {\(\rho\)}y)) {\(\cdot\)} Ahom(x + {\(\rho\)}y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},homi(x, {\(\rho\)}y)){\(\eta\)}(y)dy.
[13] ij=Q{\(\delta\)}(ej+ \nablayvj By the continuity of translations in L1(D), since Ahom\in L1(D), for all y \in Q{\(\delta\)}and B \subset\subset D we have \hat{}{\(\rho\)}\to0\to 0. B |Ahom(x + {\(\rho\)}y) - Ahom(x)|dx Integrating over Q{\(\delta\)}and using Fubini’s theorem, one obtains \hat{}\hat{}{\(\rho\)}\to0\to 0.(3.17) BQ{\(\delta\)} |Ahom(x + {\(\rho\)}y) - Ahom(x)|dy dx Consequently, for almost every x \in B, and almost every y \in Q{\(\delta\)}, Ahom(x + {\(\rho\)}y){\(\rho\)}\to0\to Ahom(x).(3.18) Let now x \in B be such a point, and let then w{\(\rho\)}k\in H10(Q{\(\delta\)}) be solutions for k \in {1, . . . , d} to -\nablay{\(\cdot\)} Ahom(x + {\(\rho\)}y)\nablayw{\(\rho\)}k(y) = \nablay{\(\cdot\)} Ahom(x + {\(\rho\)}y)ek. Estimate (3.18) implies that the assumptions of Lemma 3 are satisfied, so that w{\(\rho\)}k\to wkin H1(Q{\(\delta\)}), where wk10(Q{\(\delta\)}) to is the unique weak solution in H -\nablay{\(\cdot\)} Ahom(x)\nablaywk(y) = \nablay{\(\cdot\)} Ahom(x)ek= 0.(3.19) Hence, for all i, j \in {1, . . . , d}
[14] ij=Q{\(\delta\)}(ej+ \nablayw{\(\rho\)}j(y)) {\(\cdot\)} Ahom(x + {\(\rho\)}y)(ei+ \nablayw{\(\rho\)}i(y)){\(\eta\)}(y)dy \hat{} {\(\rho\)}\to0\to\hat{}(ej+ \nablaywj(y)) {\(\cdot\)} Ahom(x)(ei+ \nablaywi(y)){\(\eta\)}(y)dy Q{\(\delta\)} =
[15] ij since wk\equiv 0 is the trivial solution to (3.19). This concludes the proof of the theorem. The adaptation of the proof of Theorem 3 to cover Theorem 5 is straightforward, and we leave the details to the reader. To deduce the convergence of the direct approach with windowing and filtering from Theorems 4 and 5, one proceeds as in Subsection 2.3. The proof of convergence of the Petrov-Galerkin version of the dual approach with oversampling is slightly different than in Subsection 2.4. The key observation is that this ”multiscale” Petrov-discontinuous Galerkin method can be interpreted as a simple Galerkin method for an approximate homogenized problem, which makes the analysis much more elementary than expected. Let us recall there is a one-to-one mapping M{\(\delta\)},H,{\(\epsilon\)},hMsFEMfrom VHto V{\(\delta\)}H,{\(\epsilon\)},h. The Petrov-discontinuous Galerlin method is as follows: Find u{\(\delta\)}H,{\(\epsilon\)},h\in V{\(\delta\)}H,{\(\epsilon\)},hsuch that for all vH\in VH, \hat{} DH,{\(\epsilon\)},h= f, vH H-1(D),H10(D). \nablavH{\(\cdot\)} A{\(\epsilon\)}\nablaHu{\(\delta\)} Recall that V{\(\delta\)}H,{\(\epsilon\)},h/\in H1(D), and that for all v{\(\delta\)}H,{\(\epsilon\)},h=JHi=1{\(\nu\)}H,i{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i, \nablaHu{\(\delta\)}H,{\(\epsilon\)},hdenotes the broken gradient JH{\(\nu\)}H,i\nabla{\(\psi\)}{\(\delta\)}H,{\(\epsilon\)},h,i. We now use the one-to-one mapping M{\(\delta\)},H,{\(\epsilon\)},hMsFEMfrom VHto V{\(\delta\)}H,{\(\epsilon\)},h, so that one may write i=1 \nablau{\(\delta\)}H,{\(\epsilon\)},h= M{\(\delta\)},H,{\(\epsilon\)},hMsFEM(u{\(\delta\)}H,{\(\epsilon\)},h) := \nablau{\(\delta\)}H,{\(\epsilon\)},h+ \nablau{\(\delta\)}H,{\(\epsilon\)},h{\(\cdot\)} \nabla{\(\Phi\)}{\(\delta\)}H,{\(\epsilon\)},h. ESAIM: PROCEEDINGS83 for some u{\(\delta\)}H,{\(\epsilon\)},h\in VH. This turns the equation for u{\(\delta\)}H,{\(\epsilon\)},hinto an equation for u{\(\delta\)}H,{\(\epsilon\)},h: Find u{\(\delta\)}H,{\(\epsilon\)},h\in VHsuch that for all vH\in VH,\hat{} DH H-1(D),H10(D), \nablavH{\(\cdot\)} A{\(\delta\)}H,{\(\epsilon\)},h\nablau{\(\delta\)}H,{\(\epsilon\)},h= f, v where A{\(\delta\)}H,{\(\epsilon\)},his defined as A{\(\delta\)}H,{\(\epsilon\)},h:=TiHA{\(\epsilon\)}(Id + \nabla{\(\Phi\)}{\(\delta\)}H,{\(\epsilon\)},h)1TiH. IH i=1 For this equation to be well-posed we need A{\(\delta\)}H,{\(\epsilon\)},hto be coercive. Indeed, there exists some {\(\delta\)} > 1 small enough such that for all H, {\(\epsilon\)}, h > 0, A{\(\delta\)}H,{\(\epsilon\)},his uniformly coercive on D. The proof of this property is similar to the proof of the corresponding result in Theorem 4, and relies on Meyers’ estimates. We leave the details to the reader. We shall now show that limH\to0lim{\(\epsilon\)}\to0limh\to0u{\(\delta\)}H,{\(\epsilon\)},h- uhom H1(D)= 0.(3.20) The limit h \to 0 is standard, and we let A{\(\delta\)}H,{\(\epsilon\)}and u{\(\delta\)}H,{\(\epsilon\)}denote the limits of A{\(\delta\)}H,{\(\epsilon\)},hand u{\(\delta\)}H,{\(\epsilon\)},h. Furthermore we denote by uH,homthe weak solution in VHto -\nabla {\(\cdot\)} Ahom\nablauH,hom= f. Then, by the triangle inequality \nablau{\(\delta\)}H,{\(\epsilon\)}- \nablauhom L2(D)\leq \nablau{\(\delta\)}H,{\(\epsilon\)}- \nablauH,hom L2(D)+ \nablauH,hom- \nablauhom L2(D). We then treat the two terms of the r. h. s. separately and start with the first one. The function u{\(\delta\)}H,{\(\epsilon\)}- uH,hom is the weak solution in VHto -\nabla {\(\cdot\)} A{\(\delta\)}H,{\(\epsilon\)}\nabla(u{\(\delta\)}H,{\(\epsilon\)}- uH,hom) = -\nabla {\(\cdot\)} (Ahom- A{\(\delta\)}H,{\(\epsilon\)})\nablauH,hom, so that \nablau{\(\delta\)}H,{\(\epsilon\)}- \nablauH,hom L2(D)(Ahom- A{\(\delta\)}H,{\(\epsilon\)})\nablauH,hom L2(D) \leq (Ahom- A{\(\delta\)}H,{\(\epsilon\)})\nablauhom L2(D)+ 2{\(\beta\)} \nablauH,hom- \nablauhom L2(D). The first term of the r. h. s. converges to zero as {\(\epsilon\)} and H go to zero by the dominated convergence theorem and using the fact that the following convergence holds pointwise, as above, limH\to0lim{\(\epsilon\)}\to0A{\(\delta\)}H,{\(\epsilon\)}= Ahom. The second term is standard and vanishes due to the convergence of the Galerkin method: limH\to0\nablauH,hom- \nablauhom L2(D)= 0. We have thus proved (3.20). The rest of the proof is as in Subsection 2.4. The case of the (discontinuous) Galerkin version of the dual approach can be dealt with the same way, although care has to be taken for the right hand side (the same way as in Subsection 2.4 where some fHis introduced) Comments are in order. In this subsection we have shown the convergence of the direct and dual approaches of numerical homogenization when combined with windowing, filtering, and oversampling. A crucial question in the analysis (and in practice) is under which conditions the approximate homogenized matrices A{\(\delta\)}H,{\(\epsilon\)},hare coercive. We’ve shown this is the case for {\(\delta\)} close enough to 1. Other conditions may ensure this as well (such as H small and Ahomsmooth). Yet this can be an issue for the numerical practice. It is also worth mentioning that the analysis of convergence of the Petrov-discontinuous Galerkin method does not require any stabilization. This is rather unusual since discontinuous Galerkin methods normally require stabilization to converge to the right solution. The reason for this is that although the method is written in terms of a Petrov-discontinuous Galerkin formulation, it is equivalent (using the one-to-one mapping M{\(\delta\)},H,{\(\epsilon\)},hMsFEM) to a simple conforming Galerkin method on a different equation. The convergence analysis for the latter is then standard, and implies the convergence of the former – this structure is unusual and peculiar to the homogenization problem under consideration. Here we have not considered other boundary conditions for the numerical correctors than homogeneous Dirichlet boundary conditions. This owes to the fact that windowing aims at reducing the effect of boundary conditions, so that the precise choice of the boundary conditions shouldn’t affect the convergence result (indeed, the results hold as well with the boundary conditions of Remark 3). A last and important remark is the following. In (3.8), the quantity which is averaged is the energy density (ej+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}j(x, y)) {\(\cdot\)} A{\(\epsilon\)}(y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}i(x, y)). Yet other choices than this one are possible, and ej{\(\cdot\)} A{\(\epsilon\)}(y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}i(x, y)) would do the job as well. As opposed to the approach without windowing, these two choices give rise to two different approximations, although both converge within the analytical framework (the proof adapt in a straightforward way). Of course, for symmetric diffusion coefficients, it makes ”more sense” to use a symmetric formula. For non-symmetric diffusion coefficients however, there is a more subtle way to generalize this formula. Indeed a property which should be preserved at the level of the approximation is the fact that homogenizing the adjoint problem is equivalent to taking the adjoint of the homogenized problem (in other words, we always have (Ahom)T= (AT)hom). This property is only satisfied at the level of the approximation formula provided the quantity to be averaged is (ej+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}j(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)(ei+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}i(x, y)), where v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}j is the unique weak solution in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD) to -\nabla {\(\cdot\)} AT{\(\epsilon\)}(x + y)(ek+ \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}k(x, y)) = 0in Q{\(\delta\)}{\(\rho\)}\cap T-xD. This property will be important in the next section. 4.Reduction of the resonance error by zero-order regularization We propose now a refinement of the method of windowing to reduce more efficiently the resonance error. The method is based on the introduction of a zero-order term in the corrector equation, and the use of a suitable averaging mask. In this section we describe the method on a prototypical problem, and shall combine it with numerical homogenization methods in the following section. We assume that the homogenized matrix Ahomis given by the asymptotic formula: for all i, j \in {1, . . . , d}, ej{\(\cdot\)} Ahomei= limR\to\inftyQR(ej+ \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}i),(4.1) where for all k \in {1, . . . , d}, {\(\phi\)}kis the weak solution to the corrector equation on Rd -\nabla {\(\cdot\)} A(ek+ \nabla{\(\phi\)}k) = 0.(4.2) If the matrix A is for instance periodic, quasi-periodic, or more generally stationary ergodic, this corrector equation is well-posed (uniqueness follows from the condition of sublinearity at infinity) and the definition of homogenized coefficients makes sense. ESAIM: PROCEEDINGS85 In this form, the numerical homogenization methods presented so far essentially consist in replacing {\(\phi\)}kby {\(\phi\)}kR, unique weak solution in H10(QR) to (4.2) on QRfor some large R > 0, and in approximating Ahomby ej{\(\cdot\)} AR,Lei=(ej+ \nabla{\(\phi\)}jR) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iR){\(\eta\)}L(4.3) \hat{} QR for some averaging mask of support QLof size R/2 \leq L \leq R. As already mentioned, the error in the periodic case is: |Ahom- AR,L|1/R.(4.4) The objective of this section is to obtain a better convergence rate. 4.1. Description of the method, and analysis in the periodic case The starting point of the method is the following observation: solving (4.2) on a bounded domain QR requires to introduce artificial boundary conditions (namely here homogeneous Dirichlet boundary conditions). The associated error then propagates from the boundary \partialQRto the inside of QRdue to the poor decay of the Green’s function associated with the operator -\nabla {\(\cdot\)} A\nabla (which is algebraic), so that windowing ”alone” is not that efficient. We need to find a way to localize the error to a boundary layer in the neighborhood of \partialQR. This can be achieved by adding an absorbing term in the equation, namely a zero-order term. For all T > 0 we consider the ”regularized” corrector equation on Rd: for all k \in {1, . . . , d}, T-1{\(\phi\)}kT- \nabla {\(\cdot\)} A(ek+ \nabla{\(\phi\)}kT) = 0.(4.5) The (unique) weak solution to this equation {\(\phi\)}kTis much easier to approximate than {\(\phi\)}kbecause the Green’s-1- \nabla {\(\cdot\)} A\nabla decays exponentially fast in terms of distance measured in\sqrt{}T . In particular, the difference on QLbetween {\(\phi\)}kTand a solution computed on a bounded domainR-L\sqrt{}T(for R \geq L). In order to deal with non-symmetric matricesefficiently as well, we introduce the adjoint problem and consider the regularized adjoint corrector equation on function associated with the operator T units of QRis essentially of infinite order in terms of Rd: T-1{\(\phi\)}kT- \nabla {\(\cdot\)} AT(ek+ \nabla{\(\phi\)}kT) = 0.(4.6) The use of the adjoint problem indeed follows Tartar’s seminal ideas, whereas the introduction of the zero order term is somehow inspired by the analysis of the corrector equation in the stochastic case. The approximation of Ahomwe shall consider is then given for all i, j \in {1, . . . , d} by ej{\(\cdot\)} AT,R,Lei:=(ej+ \nabla{\(\phi\)}jT,R) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT,R){\(\eta\)}L,(4.7) \hat{} QR where {\(\eta\)}Lis a suitable averaging mask, and {\(\phi\)}kT,R, {\(\phi\)}kT,R\in H10(QR) are the weak solutions in QRto (4.5) and (4.6), respectively. Before we turn to the analysis of the periodic case, let us make the notion of averaging mask more precise. Definition 7.A function {\(\eta\)} : [-1/2, 1/2] \to R+is said to be a filter of order p \geq 0 if (i) {\(\eta\)} \in Cp([-1/2, 1/2]) \cap Wp+1,\infty((-1/2, 1/2)), (ii)\'{}1/2-1/2{\(\eta\)}(x)dx = 1, (iii) {\(\eta\)}(k)(-1/2) = {\(\eta\)}(k)(1/2) = 0 for all k \in {0, . . . , p - 1}. The associated mask {\(\eta\)}L: [-L/2, L/2]d\to R+in dimension d \geq 1 is then defined for all L > 0 by {\(\eta\)}L(x) := L-d{\(\eta\)}(L-1xi), d i=1 where x = (x1, . . . , xd) \in Rd. In the periodic case, we then have the following quantitative convergence result: Theorem 6.Let d \geq 2, A \in M{\(\alpha\)}{\(\beta\)}be Q-periodic, {\(\eta\)} be a filter of order p \geq 0, and Ahomand AT,R,Lbe the2TR, R \geq L \sim R \sim R - L. Then, homogenized matrix and its approximation (4.7) respectively, where R there exists c > 0 depending only on {\(\alpha\)}, {\(\beta\)} and d such that we have |AT,R,L- Ahom|L-(p+1)+ T-2+ T1/4exp -cR - L\sqrt{}T.(4.8) Remark 4.This result is an extension of [29, Theorem 3.1] which allows to cover the case of non-symmetric diffusion coefficients as well. If in (4.7) we had used the (primal) corrector instead of the corrector associated with the adjoint problem, the term T-2would be replaced by T-1for non-symmetric matrices. The starting point of the proof is the decomposition of the error into three contributions: |AT,R,L- Ahom| \leq |AT,R,L- AT,L| + |AT,L- AT| + |AT- Ahom|, where\hat{}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT){\(\eta\)}L, ej{\(\cdot\)} AT,Lei:=(ej+ \nabla{\(\phi\)}T Rd and ej{\(\cdot\)} ATei:=limR\to\inftyQR(ej+ \nabla{\(\phi\)}jT) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT)(4.9) =(ej+ \nabla{\(\phi\)}jT) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT), Q by periodicity of A, {\(\phi\)}iTand {\(\phi\)}jT.k The first contribution measures the error due to the use of boundary conditions to approximate {\(\phi\)}kTand {\(\phi\)}T on a bounded domain. As already mentioned, this term is exponentially small due to the decay of the Green’s function associated with the Helmholtz operators T-1- \nabla {\(\cdot\)} A\nabla and T-1- \nabla {\(\cdot\)} AT\nabla on the whole space (and the maximum principle). The second contribution measures the error made in approximating the average of the periodic function (ej+ \nabla{\(\phi\)}jT) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT) (recall that {\(\phi\)}kT, {\(\phi\)}kT, {\(\phi\)}k, and {\(\phi\)}kare periodic) by the average on QRwith the mask {\(\eta\)}L. By Fourier analysis, one can show that this approximation error decays as in L-(p+1), where p is the order of the filter, so that this contribution is not the limiting one for p sufficiently large. Both errors are analyzed in detail in [29]. The last contribution is the so-called systematic error, which is a consequence of the fact that we have modified the corrector equation by a zero order term, and therefore introduced a bias which is controlled by the parameter T . This is the only place which changes in the proof with respect to the symmetric case treated in [29]. We shall make use of the weak forms of the adjoint corrector and corrector equations: for all v \in H1#(Q) \hat{}j) = 0,(4.10) Q \nablav {\(\cdot\)} AT(ej+ \nabla{\(\phi\)} \hat{} \nablav {\(\cdot\)} A(ei+ \nabla{\(\phi\)}i) = 0.(4.11) ESAIM: PROCEEDINGS87 By definition of ATand Ahomand using (4.11) for v = {\(\phi\)}jand then v = {\(\phi\)}j, we obtain ej{\(\cdot\)} (AT- Ahom)ei=(ej+ \nabla{\(\phi\)}jT) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT) -\hat{}Q(ej+ \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}i) \hat{} Q =(ej+ \nabla{\(\phi\)}jT) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT) -\hat{}Q(ej+ \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}i) \hat{} Q =(\nabla{\(\phi\)}jT- \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT) +\hat{}Q(ej+ \nabla{\(\phi\)}j) {\(\cdot\)} A(\nabla{\(\phi\)}iT- \nabla{\(\phi\)}i) \hat{} Q =(\nabla{\(\phi\)}jT- \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT) +\hat{}Q(\nabla{\(\phi\)}iT- \nabla{\(\phi\)}i) {\(\cdot\)} AT(ej+ \nabla{\(\phi\)}j). \hat{} Q Using then (4.10) with v = {\(\phi\)}iT- {\(\phi\)}i, and (4.11) with v = {\(\phi\)}jT- {\(\phi\)}j, this turns into ej{\(\cdot\)} (AT- Ahom)ei=(\nabla{\(\phi\)}jT- \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iT) -\hat{}(\nabla{\(\phi\)}jT- \nabla{\(\phi\)}j) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}i) \hat{} QQ =(\nabla{\(\phi\)}jT- \nabla{\(\phi\)}j) {\(\cdot\)} A(\nabla{\(\phi\)}iT- \nabla{\(\phi\)}i). \hat{} Q Hence, since A \in M{\(\alpha\)}{\(\beta\)}, this turns into |AT- Ahom|\nabla{\(\phi\)}jT- \nabla{\(\phi\)}jL2(D)\nabla{\(\phi\)}iT- \nabla{\(\phi\)}iL2(D), and we need to estimate the convergence of \nabla{\(\phi\)}iTto \nabla{\(\phi\)}i(the second term is similar). Since {\(\phi\)}iT- {\(\phi\)}iis the1#(Q) to unique weak solution in H T-1({\(\phi\)}iT- {\(\phi\)}i) - \nabla {\(\cdot\)} A\nabla({\(\phi\)}iT- {\(\phi\)}i) = -T-1{\(\phi\)}i, an a priori estimate combined with Poincar\'{}e’s inequality in H1#(Q) allows to conclude that \nabla({\(\phi\)}iT- {\(\phi\)}i)L2(Q)T-1. This finaly implies the desired estimate |AT- Ahom|T-2.(4.12) Let us give a simple application of this estimate: For p \geq 3, the rate in (4.8) is controlled by the last two\sqrt{}T . A possible choice is then given by terms. In particular the last term requires T to be such that L \gg ●T = L2(ln L)-4, ●R = 3L/2, for which (4.8) reads: |AT,R,L- Ahom|R-4ln8R.(4.13) Whereas the convergence rate in estimate (4.4) is of order 1, the regularization approach yields a convergence rate in (4.13) up to order 4-. 4.2. Spectral analysis for symmetric coefficients and consistency in the stationary ergodic case The aim of this subsection is to show that the regularization method introduced in the previous subsection is consistent at the level of the homogenized coefficients for a large class of heterogeneities, and not only for the periodic case. Unlike the consistency results we have proved so far, one cannot consider here any H-convergent sequence. We shall indeed prove consistency for the restricted class of stationary ergodic coefficients (which is already quite large, and includes periodic and quasiperiodic coefficients, as well as standard examples of random inclusions etc.). This is where the second fundamental tool of the homogenization theory pops up: spectral theory. Because of the nature of the tools we shall use (essentially the spectral theorem), the argument we present here only covers the case of symmetric diffusion coefficients (for which we then have {\(\phi\)}T\equiv {\(\phi\)}T), as opposed to Theorem 6. In order for this survey to be as self-contained as possible, we quickly recall standard qualitative results in stochastic homogenization of linear elliptic equations. We refer the reader to the original papers [55] by Papanicolaou and Varadhan, and [47] by Kozlov for details (see also the monography [45]). Let (&#8486;, F, P) be a probability space, and denote by {\(\cdot\)} the associated expectation. We shall say that the family of mappings ({\(\theta\)}z)z\inRdfrom &#8486; to &#8486; is a strongly continuous measure-preserving ergodic translation group if: ●({\(\theta\)}z)z\inRdhas the group property: {\(\theta\)}0= Id (the identity mapping), and for all x, y \in Rd, {\(\theta\)}x+y= {\(\theta\)}x&#9702; {\(\theta\)}y;d, and every measurable set F \in F, {\(\theta\)}xF is measurable ●({\(\theta\)}z)z\inRdpreserves the measure: for all x \in Rand P({\(\theta\)}xF ) = P(F ); ●({\(\theta\)}z)z\inRdis strongly continuous: for any measurable function f on &#8486;, the function ({\(\omega\)}, x) \to f({\(\theta\)}x{\(\omega\)})defined on &#8486; {\(\times\)} Rdis measurable (with the Lebesgue measure on Rd); ●({\(\theta\)}z)z\inRdis ergodic: for all F \in F, if for all x \in Rd, {\(\theta\)}xF \subset F , then P(F ) \in {0, 1}. Let 0 < {\(\alpha\)} \leq {\(\beta\)} < \infty, and let A \in L2(&#8486;, M{\(\alpha\)}{\(\beta\)}). We define a stationary extension of A (still denoted by A) on Rd{\(\times\)} &#8486; as follows: A : (x, {\(\omega\)}) \to A(x, {\(\omega\)}) = A({\(\theta\)}x{\(\omega\)}). Homogenization theory ensures that the solution operator associated with -\nabla {\(\cdot\)} A(x/{\(\epsilon\)}, {\(\omega\)})\nabla converges as {\(\epsilon\)} > 0 vanishes to the solution operator of -\nabla {\(\cdot\)} Ahom\nabla for P-almost every {\(\omega\)}, where Ahomis a deterministic elliptic matrix characterized as follows. For all {\(\xi\)}, {\(\zeta\)} \in Rd, and P-almost every {\(\omega\)}, {\(\xi\)} {\(\cdot\)} Ahom{\(\zeta\)}=limR\to\inftyQR({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}(x, {\(\omega\)})) {\(\cdot\)} A(x, {\(\omega\)})({\(\zeta\)} + \nabla{\(\phi\)}{\(\zeta\)}(x, {\(\omega\)}))dx =({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A({\(\zeta\)} + \nabla{\(\phi\)}{\(\zeta\)}) , where {\(\phi\)}{\(\xi\)}: Rd{\(\times\)} &#8486; \to R is Borel measurable, is such that {\(\phi\)}{\(\xi\)}(0, {\(\cdot\)}) \equiv 0, \nabla{\(\phi\)}{\(\xi\)}is stationary, \nabla{\(\phi\)}{\(\xi\)}= 0, and {\(\phi\)}{\(\xi\)}({\(\cdot\)}, {\(\omega\)}) \in H1loc(Rd) is almost surely a distributional solution to the corrector equation -\nabla {\(\cdot\)} A(x, {\(\omega\)})({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}(x, {\(\omega\)})) = 0 in Rd,(4.14) and likewise for {\(\phi\)}{\(\zeta\)}. The proof of existence and uniqueness of these correctors is obtained by regularization, and we consider for all T > 0 the stationary solution {\(\phi\)}{\(\xi\)}Twith zero expectation to the equation T-1{\(\phi\)}{\(\xi\)}T(x, {\(\omega\)}) - \nabla {\(\cdot\)} A(x, {\(\omega\)})({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T(x, {\(\omega\)})) = 0 in Rd. This equation has an equivalent form in the probability space, to which we can apply the Lax-Milgram theorem. This formulation requires a bit of functional analysis: the stochastic counterpart of \nablai(for i \in {1, . . . , d}) is denoted by Diand defined by Dif ({\(\omega\)}) = limh\to0f ({\(\theta\)}hei{\(\omega\)}) - f({\(\omega\)})h. These are the infinitesimal generators of the d one-parameter strongly continuous unitary groups on L2(&#8486;) defined by the translations in each of the d directions. These operators commute and are closed and densely ESAIM: PROCEEDINGS89 defined on L2(&#8486;). We denote by H(&#8486;) the domain of D = (D1, . . . , Dd). This subset of L2(&#8486;) is a Hilbert space for the norm f2H= |Df|2+ f2. Since the groups are unitary, the operators are skew-adjoint so that we have the ”integration by parts” formula: for all f, g \in H(&#8486;) f Dig = - gDif . The equivalent form of the regularized corrector equation is as follows: T-1{\(\phi\)}{\(\xi\)}T- D {\(\cdot\)} A({\(\xi\)} + D{\(\phi\)}{\(\xi\)}T) = 0,(4.15) which admits a unique weak solution in {\(\phi\)}{\(\xi\)}T\in H(&#8486;), that is such that for all {\(\psi\)} \in H(&#8486;), T-1{\(\phi\)}{\(\xi\)}T{\(\psi\)} + D{\(\psi\)} {\(\cdot\)} A({\(\xi\)} + D{\(\phi\)}{\(\xi\)}T)= 0.(4.16) One may prove using the integration by parts formula that D{\(\phi\)}{\(\xi\)}Tis bounded in L2(&#8486;) and converges weakly (up2(&#8486;) to some solution {\(\Phi\)}{\(\xi\)}, which is a gradient. Using then the spectral representation of the to extraction) in L translation group we may prove uniqueness of the corrector {\(\phi\)}{\(\xi\)}(which is such that \nabla{\(\phi\)}{\(\xi\)}= {\(\Phi\)}{\(\xi\)}), see [55]. Up to here, we have not required A to be symmetric. Let Msym{\(\alpha\)}{\(\beta\)}denote the subset of symmetric matrices of M{\(\alpha\)}{\(\beta\)}, and set Asym{\(\alpha\)}{\(\beta\)}= L\infty(Rd, Msym{\(\alpha\)}{\(\beta\)}).(4.17) In the rest of this section, we shall consider that A \in L2(&#8486;, Msym{\(\alpha\)}{\(\beta\)}) so that one can appeal to spectral theory.sym), and that Ahomis also symmetric. Note that the stationary extension of A belongs to L2(&#8486;, A{\(\alpha\)}{\(\beta\)} Let L = -D {\(\cdot\)} AD be the operator defined on H(&#8486;) as a quadratic form. We still denote by L its Friedrichs extension in L2(&#8486;). This operator is a nonnegative selfadjoint operator. By the spectral theorem it admits a spectral resolution L ={\(\lambda\)}G(d{\(\lambda\)}). \hat{} R+ Recall that this allows one to define a spectral calculus, namely for all suitable function g : R+\to R, one may define the operator g(L) by g(L) :=g({\(\lambda\)})G(d{\(\lambda\)}), \hat{} R+ as one would do for symmetric matrices. We further denote by d := -D {\(\cdot\)} A{\(\xi\)} the local drift in direction {\(\xi\)}, and-1dand shall consider the projection ed= dGd of the spectral measure G onto the local drift. Since {\(\phi\)} = L \nabla{\(\psi\)} {\(\cdot\)} A\nabla{\(\chi\)} = {\(\psi\)}L{\(\chi\)} by integration by parts for all {\(\psi\)}, {\(\chi\)} \in H(&#8486;), one has formally \nabla{\(\phi\)} {\(\cdot\)} A\nabla{\(\phi\)} = (L-1d)L(L-1d) . This identity can be proved using the regularized corrector {\(\phi\)}Tand passing to the limit as T \to \infty. We then apply spectral calculus to L and the function g : {\(\lambda\)} \to1{\(\lambda\)}, which yields the following spectral identity: {\(\xi\)} {\(\cdot\)} Ahom{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}) ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} + \nabla{\(\phi\)} {\(\cdot\)} A\nabla{\(\phi\)} - 2 \nabla{\(\phi\)} {\(\cdot\)} A{\(\xi\)} ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} - \nabla{\(\phi\)} {\(\cdot\)} A\nabla{\(\phi\)} ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} - dg(L)d ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} -1{\(\lambda\)}ded({\(\lambda\)}),(4.18) \hat{} R+ where we have used the weak form of the corrector equation tested with {\(\phi\)} itself \nabla{\(\phi\)} {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}) = 0, which is a consequence of the analysis in [55]. This formula will be crucial for the analysis of the regularization method: With the help of this framework, one may prove the consistency of the regularization approach. Theorem 7.Let d \geq 2, A \in Asym{\(\alpha\)}{\(\beta\)}be a stationary random field, {\(\eta\)} be some mask, and Ahomand AT,R,Lbe the homogenized matrix and its approximation (4.7), respectively. Then, almost surely, limT \to\inftylimR\geqL\to\infty|AT,R,L- Ahom| = 0. Let {\(\xi\)} \in Rdwith |{\(\xi\)}| = 1 be fixed. By the triangle inequality we have: |{\(\xi\)} {\(\cdot\)} AT,R,L{\(\xi\)} - {\(\xi\)} {\(\cdot\)} Ahom{\(\xi\)}| \leq |{\(\xi\)} {\(\cdot\)} AT,R,L{\(\xi\)} - {\(\xi\)} {\(\cdot\)} AT,L{\(\xi\)}| + |{\(\xi\)} {\(\cdot\)} AT,L{\(\xi\)} - {\(\xi\)} {\(\cdot\)} AT{\(\xi\)}| + |{\(\xi\)} {\(\cdot\)} AT{\(\xi\)} - {\(\xi\)} {\(\cdot\)} Ahom{\(\xi\)}|,(4.19) where {\(\xi\)} {\(\cdot\)} AT,L{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T){\(\eta\)}L,(4.20) \hat{} Rd {\(\xi\)} {\(\cdot\)} AT{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) .(4.21) The first term |{\(\xi\)} {\(\cdot\)}AT,R,L{\(\xi\)} -{\(\xi\)} {\(\cdot\)}AT,L{\(\xi\)}| of the r. h. s. of (4.19) is a measure of the error due to boundary conditions and scales likeR - L\sqrt{}T,(4.22)\sqrt{}T \to \infty. |AT,R,L- AT,L|T3/4exp - c see [30, Proposition 2.8 and Remark 2.9]. In particular, it vanishes in the limit (R - L)/ The second term |{\(\xi\)} {\(\cdot\)} AT,L{\(\xi\)} - {\(\xi\)} {\(\cdot\)} AT{\(\xi\)}| is a measure of the error between an average in physical space and the expectation. By the ergodic theorem, this vanishes almost surely as L \to \infty. The last term |{\(\xi\)} {\(\cdot\)} AT{\(\xi\)} - {\(\xi\)} {\(\cdot\)} Ahom{\(\xi\)}| is the systematic error. We shall use spectral theory to show that it converges to zero as T \to \infty. As for (4.18), using the regularized corrector equation (4.16), the spectral theorem indeed yields {\(\xi\)} {\(\cdot\)} AT{\(\xi\)}={\(\xi\)} {\(\cdot\)} A{\(\xi\)} + \nabla{\(\phi\)}T{\(\cdot\)} A\nabla{\(\phi\)}T+ 2 \nabla{\(\phi\)}T{\(\cdot\)} A{\(\xi\)} ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} + \nabla{\(\phi\)}T{\(\cdot\)} A\nabla{\(\phi\)}T- 2 \nabla{\(\phi\)}T{\(\cdot\)} A\nabla{\(\phi\)}T- 2T-1{\(\phi\)}2T ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} -{\(\lambda\)}(T-1+ {\(\lambda\)})2ded({\(\lambda\)}) - 2T-1\hat{}R+1(T-1+ {\(\lambda\)})2ded({\(\lambda\)}) \hat{} R+ ={\(\xi\)} {\(\cdot\)} A{\(\xi\)} -2T-1+ {\(\lambda\)}(T-1+ {\(\lambda\)})2ded({\(\lambda\)}), \hat{} ESAIM: PROCEEDINGS91 so that1ded({\(\lambda\)}),(4.23) {\(\xi\)} {\(\cdot\)} (AT- Ahom){\(\xi\)} = T-2{\(\lambda\)}(T-1+ {\(\lambda\)})2 \hat{} R+ which vanishes as T \to \infty by the dominated convergence theorem (since {\(\lambda\)}-1is integrable for ded): limT \to\infty|AT- Ahom| = 0.(4.24) This concludes the proof of Theorem 7. As a by-product of this analysis we also have consistency of the regularized corrector almost surely: limT \to\inftylimR\geqL\to\inftyQL|\nabla{\(\phi\)}{\(\xi\)}T,R- \nabla{\(\phi\)}{\(\xi\)}|2= 0.(4.25) On the one hand, by ergodicity, limR\to\inftyQR|\nabla{\(\phi\)}T,R- \nabla{\(\phi\)}|2= |\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}|2(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) . On the other hand, using the symmetry of A (in the third identity) and the weak form of the corrector equation for {\(\phi\)}{\(\xi\)}(in the fourth identity), that is for all D{\(\psi\)} \in L2(&#8486;), D{\(\psi\)} {\(\cdot\)} A({\(\xi\)} + D{\(\phi\)}{\(\xi\)}) = 0, we obtain {\(\xi\)} {\(\cdot\)} AT{\(\xi\)} - {\(\xi\)} {\(\cdot\)} Ahom{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) - ({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) =(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) + ({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) =(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}T) + (\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) =({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) - (\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\xi\)}) =(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) {\(\cdot\)} A(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) .(4.26) The combination of this inequality and this identity with (4.24) yields (4.25). Remark 5.With some additional work we can even consider some ”diagonal” extraction in a weaker norm. In particular, in the regime R2TR, R \geq L \sim R \sim R - L, we have limT,R,L\to\infty|AT,R,L- Ahom| = 0. 4.3. Convergence rates in the stochastic case with finite correlation-length In the case of stationary random symmetric coefficients with finite correlation-length cl, that is coefficients A such that for all x, y \in Rd, A(x) and A(y) are independent random fields if |x - y| \geq cl, we may even give quantitative results for the regularization method. This covers for instance the case of a random checkerboard (say when the colors are independent and identically distributed random variables) or the case of random inclusions whose centers are distributed according to a Poisson point process (or hard-core Poisson point process to avoid overlapping inclusions). In this case, Otto and the author proved the following convergence rates for all d = 2, 3, 4 in [37]: |AT,R,L- Ahom|2d = 2 :L-2lnqT,2 < d \leq 4 : L-d,(4.27) for R = 2L, T = L ln2L, and some q > 0 depending only on {\(\alpha\)}, {\(\beta\)}. In particular, these estimates are expected to be optimal (they have the central limit theorem scaling) up to the logarithmic correction for d = 2. They improve the estimates by Bourgeat and Piatnitski in [12] and by E, Ming, and Zhang in [18], which essentially show (for more general statistics on the coefficients however) that there exists some non explicit exponent {\(\gamma\)} > 0 such that |AT,R,L- Ahom|2L-{\(\gamma\)}(4.28) provided d > 2. The core ingredient in their proof is the very insightful contribution by Yurinski\breve{}ı [60]. However, as explained in the introduction of [38], even if all the steps of Yurinski\breve{}ı’s work yielded the expected optimal exponents, the method to pass from the results of [60] to (4.28) does not allow to obtain the optimal convergence rate, so that necessarily {\(\gamma\)} < d. The proof of (4.27) is rather long, and go beyond the scope of this survey. Let us also mention that in the case of discrete elliptic equations on Zdwith random conductivities, the picture is even more complete, and we refer the reader to the series of papers [30, 35, 36, 38, 39]. It is also worth noting that the scaling of (4.27) is better than the scaling of (4.4). 4.4. Improving the convergence rate by Richardson extrapolation In (4.27), the result is only stated up to d = 4 because it does not hold for d > 4. Let us discard the error due to boundary conditions, which is controlled by (4.22) and can be made decay at any order in L (provided\sqrt{}T ). We focus instead on the error term |AT,L- Ahom|2, which splits into two contributions: the R - L \gg random error (AT,Lfluctuates around its expectation AT,L, which by stationarity is nothing but AT) and a systematic error (the difference between the expectation of AT,L= ATand Ahom). As proved in [37], the random error scales as the central limit theorem in any dimension |AT,L- AT|2d = 2: L-2lnqL,d > 2: L-d, for some q > 0 depending only on {\(\alpha\)}, {\(\beta\)}, provided L \leq T (which is compatible with the requirement R-L \gg\sqrt{}T ). The scaling of the systematic error is however different: |AT- Ahom|d = 3 : T-3/2,d = 4 : T-2ln T,(4.29) d = 2 : T-1lnqT d > 4 : T-2, so that the choice T = L ln2L (which is convenient to control the error due to boundary conditions via the exponential decay of the Green’s function of the Helmholtz equation) only yields the central limit theorem scaling in (4.27) up to d = 4. Recall that for stationary random fields A, we have the spectral formula (4.23) for the systematic error: {\(\xi\)} {\(\cdot\)} (AT- Ahom){\(\xi\)} = T-21{\(\lambda\)}(T-1+ {\(\lambda\)})2ded({\(\lambda\)}). \hat{} R+ Hence, the best one can hope for the convergence of ATto Ahomis indeed T-2, which holds provided {\(\lambda\)} \to {\(\lambda\)}-3 is dedintegrable. What (4.29) implies is that {\(\lambda\)} \to {\(\lambda\)}-3is dedintegrable for d > 4 in the case of stationary ESAIM: PROCEEDINGS93 coefficients with finite correlation length. In the periodic case, since there is a spectral gap, there exists {\(\mu\)} > 0 such that ed((0, {\(\mu\)})) = 0, and {\(\lambda\)} \to {\(\lambda\)}-3is dedintegrable in any dimension, so that (4.23) implies (4.12) (provided the matrix is symmetric). In particular, using the approximation ATof Ahom, one cannot benefit from the fact that {\(\lambda\)} \to {\(\lambda\)}-kmay be dedintegrable for some k > 3. In order to benefit from this, one needs to introduce other approximation formulas of Ahomthan AT. Ideally we are looking for approximations AT,kof Ahomwhich are such that |AT,k- Ahom|T-2k1{\(\lambda\)}(T-1+ {\(\lambda\)})2kded({\(\lambda\)}).(4.30) \hat{} R+ In [34], Mourrat and the author defined such a family of approximations of Ahomby induction in terms of their spectral representations. Yet, this approach is only useful if these approximations admit an explicit counterpart in physical space. The formulas introduced in [34] do indeed have the following representation: {\(\xi\)} {\(\cdot\)} AT,k{\(\xi\)} = ({\(\xi\)} + \nabla{\(\phi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}T) + T-1{\(\gamma\)}ij{\(\phi\)}2-iT{\(\phi\)}2-jT, k-1 i,j=0 where the coefficients {\(\gamma\)}ijare defined by induction. To be concrete, the first three approximations of Ahomare given by: {\(\xi\)} {\(\cdot\)} AT,1{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}T) , {\(\xi\)} {\(\cdot\)} AT,2{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}T) - 3T-1{\(\phi\)}2T- 2T-1{\(\phi\)}2T /2+ 5T-1{\(\phi\)}T{\(\phi\)}T /2, {\(\xi\)} {\(\cdot\)} AT,3{\(\xi\)}=({\(\xi\)} + \nabla{\(\phi\)}T) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}T) -559T-1{\(\phi\)}2T- 8T-1{\(\phi\)}2T /2-49T-1{\(\phi\)}2T /4 +413T-1{\(\phi\)}T{\(\phi\)}T /2-229T-1{\(\phi\)}T{\(\phi\)}T /4+103T-1{\(\phi\)}T /2{\(\phi\)}T /4. For general stationary ergodic coefficients, the approximation formulas AT,kare consistent with Ahomin the sense that limT \to\inftyAT,k= Ahom, which can be proved using the Lebesgue dominated convergence theorem as for the convergence of ATin Theorem 7. Note that in order to approximate AT,kin practice one needs to solve the regularized corrector equation for k different zero order terms of magnitude 2k-iT-1for i = 1, . . . k - 1. The associated computable approximations AT,k,R,Lof Ahomare given by {\(\xi\)} {\(\cdot\)} AT,k,R,L{\(\xi\)} :=({\(\xi\)} + \nabla{\(\phi\)}T,R) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}T,R){\(\eta\)}L+ T-1{\(\gamma\)}ij{\(\phi\)}2-iT,R{\(\phi\)}2-jT,R{\(\eta\)}L, \hat{}k-1\hat{} QRi,j=0QR where {\(\eta\)}Lis a suitable averaging mask, and {\(\phi\)}T-i,R\in H10(QR) are the weak solutions in QRto (4.5) with the corresponding magnitude of the zero order term. In the case of (symmetric) periodic coefficients, one may conclude using the spectral gap and (4.30) that for all k \geq 1, |AT,k- Ahom|T-2k. As proved in [34], (4.8) is then replaced for all k \geq 1 by |AT,k,R,L- Ahom|L-(p+1)+ T-2k+ T1/4exp -c2kR - L\sqrt{}T,(4.31) provided {\(\eta\)} is a mask of order p. This allows to drastically reduce the resonance error at the level of the homogenized coefficients. In the case of a discrete elliptic equation on Zdwith i. i. d. bond conductivities, Neukamm, Otto and the author proved in [36] (see also [34]) that for all d \geq 2, there exists kd\geq 1 such that for all k \geq kdwe have |AT,k,R,L- Ahom|2L-d(4.32) provided R = 2L and T = L ln2L – which the central limit theorem scaling in any dimension. It is to be expected that the result will hold for continuous equations as well, which would extend the results (4.27) of [37] to dimensions larger than 4. In this subsection, we have introduced a family of approximation formulas using spectral calculus, which is proved to drastically reduce the resonance error at the level of the homogenized coefficients. Yet, the use of spectral calculus requires the matrix A to be symmetric and it is not clear how to generalize the approximation formulas AT,kto non-symmetric coefficients. To this aim, we give now another interpretation of these approximation formulas, which will allow us to introduce a second class of approximation formulas which readily generalizes to the non-symmetric case. In terms of their spectral representation, AT,kcan be seen as extrapolations of Ahom. In particular, it may make sense to try to extrapolate directly in physical space, and we consider Richardson extrapolations defined by the following induction rule: \tilde{}AT,1= AT, and for all k \geq 1, \tilde{}AT,k+1:=12k- 1(2k\tilde{}A2T,k- \tilde{}AT,k). As shown in [35], in terms of spectral representation, these extrapolation formulas satisfy | \tilde{}AT,k- Ahom|T-(k+1)1{\(\lambda\)}(T-1+ {\(\lambda\)})k+1ded({\(\lambda\)}). \hat{} R+ so that in the case of symmetric periodic matrices, we have using the spectral gap of the elliptic operator: | \tilde{}AT,k- Ahom|T-(k+1).(4.33) The extrapolation formulas obtained by Richardson extrapolation are particularly interesting because it also makes sense to use them in the nonsymmetric case (where ATis defined by the asymptotic formula (4.9) using the corrector of the adjoint problem). As shown in [32], one can bypass the use of spectral theory in the periodic case, and directly prove by PDE arguments that (4.33) holds as well for nonsymmetric matrices. In the stochastic case however, it is not clear how and even whether these convergence rates can be understood without spectral theory. This idea of using Richardson extrapolation is very fruitful and will be applied to the corrector itself in Subsection 5.4. 4.5. Numerical tests In this subsection, we quickly illustrate the efficiency of the regularization method and the sharpness of the analysis. We start with some periodic examples to check the sharpness of the analysis and then turn to more challenging cases such as some quasi-periodic and stochastic cases. 4.5.1. Discrete periodic example To check the validity of Theorem 6 (and of its generalization (4.31)) in the regime of R large, we have considered the example of a discrete elliptic equation: -\nabla\ast{\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}) = 0in Z2,(4.34) ESAIM: PROCEEDINGS95 Figure 1.Periodic cell in the discrete case where for all u : Z2\to R, \nablau(x) :=u(x + e1) - u(x)u(x + e2) - u(x), \nabla\astu(x) :=u(x) - u(x - e1)u(x) - u(x - e2), and A(x) := diag [a(x, x + e1), a(x, x + e2)] . The matrix A is [0, 4)2-periodic, and sketched on a periodic cell on Figure 1. In the example considered, a(x, x + e1) and a(x, x + e2) represent the conductivities 1 or 100 of the horizontal edge [x, x + e1] and the vertical edge [x, x + e2] respectively, according to the colors on Figure 1. The homogenization theory for such discrete elliptic operators is similar to the continuous case (see for instance [58] in two dimensions, and [2] in the general case). Since the periodic pattern is invariant by the rotation R of angle {\(\pi\)}/2, the homogenized matrix satisfies RAhom= Ahom. In dimension d = 2, this implies that Ahomis a multiple of the identity. It can be evaluated numerically (note that we do not make any other error than the machine precision). Its numerical value is Ahom= 26.240099009901 . . .. We have considered the first two approximations formulas AT,1,R,Land AT,2,R,Lof Ahom. In all the cases treated, we’ve taken L = R/3. For the approximation AT,1,R,L, we have tested the following parameters: ●Four values for the zero-order term: T = \infty (no zero-order term), T \sim R, T \sim R3/2, and T \sim R7/4; ●Two different filters: orders p = 0 (no filter) and p = \infty. For the approximation AT,2,R,L, we have tested the following parameters: ●One value of the zero-order term: T \sim R3/2; ●Filter of infinite order p = \infty. The theoretical predictions in terms of convergence rate of AT,k,R,Lto Ahomin function of R are gathered and compared to the results of numerical tests in Table 3. More details are also given on Figures 2–5, where the overall error Error(k, T, R) := |Ahom- AT,k,R,L| is plotted in log scale in function of R. Let us quickly comment on the values of T in Figures 2–5. For the four dependences of T upon R, we have chosen the prefactors so that their values roughly coincide for R = 25 (that Table 3.Order of convergence: predictions and numerical results. T= \inftyT\sim RT\sim R3/2T\sim R7/4 k=1pred.testpred.testpred.testpred.test p = 011111111 p= \infty112233.13.53.4 k=2pred.test p= \infty65.2 is for 25 periodic cells per dimension): T=R/25, T=(4R)3/2/1000, T=(4R)7/4/5000, T=(4R)2/(25 ln4(4R)). The numerical result confirm the analysis, and perfectly illustrate the specific influences of the three parameters k, p and T . 4.5.2. Continuous periodic example We consider the following matrix A : R2\to Msym{\(\alpha\)}{\(\beta\)}, A(x) =2 + 1.8 sin(2{\(\pi\)}x1)2 + 1.8 cos(2{\(\pi\)}x2)+2 + sin(2{\(\pi\)}x2)2 + 1.8 cos(2{\(\pi\)}x1)Id,(4.35) used as benchmark tests in [44]. In this case, {\(\alpha\)} \simeq 0.35, {\(\beta\)} \simeq 20.5, and Ahom\simeq 2.75 Id. We take L = R/3, T = R/10 and a filter of order 2. The global error |AT,1,R,L- Ahom| and the error without zero order term and without filtering are plotted on Figures 6 & 7. Without zero-order term, the convergence rate is R-1as expected, and the use of a filtering method reduces the prefactor but does not change the rate. With the zeroorder term and the filtering method, the apparent convergence rate is R-3(note that the asymptotic theoretical rate R-2is not attained yet), which coincides with the convergence rate associated with filters of order 2. This is in agreement with the tests in the discrete case, and confirms the analysis. 4.5.3. Continuous quasiperiodic example We consider the following coefficients: A(x) =4 + cos(2{\(\pi\)}(x1+ x2)) + cos(2{\(\pi\)}\sqrt{}2(x1+ x2))006 + sin2(2{\(\pi\)}x1) + sin2(2{\(\pi\)}\sqrt{}2x1).(4.36) In this case, the homogenized coefficients are not easy to compute. They can only be extrapolated. We have taken for the approximation of the homogenized coefficients (that we call coefficient of reference) the output of the computation with k = 1, T = R/100 and R = 52. Although this may introduce a bias in favour of the proposed strategy, it can be checked a posteriori: the method without zero-order term and without filtering is expected to converge at a rate R-1. This is effectively what we observe on Figure 8 using this coefficient of reference. Instead, if we use as a reference the output of the computation for R = 52 without zero-order term nor filtering, then we observe a super-linear convergence which is artificial (see Figure 8). With the proposed method, as can be seen on Figure 9, the rate of convergence seems to be much better (the slope of the straight line is -5). The reason for this fact is that there should be a spectral gap in this case as well. ESAIM: PROCEEDINGS97 0.511.522.533.5 -4 -3.5 -3 -2.5 -2 -1.5logError(1,\infty -1 -0.5 0,R 0.5) 1 log 4R Figure 2.Absolute error in log scale without zero order term, no filter (slope -1), infinite order filter (slope -1, better prefactor). 0.511.522.533.5-8 -7 -6 -5 -4 -3gError(1,T -2 -1,R 0 1 log 4R lo ) Figure 3.Absolute error in log scale for T = R/25, no filter (slope -1), infinite order filter (slope -2). 0.511.522.533.5 -10 -8 -6 -4logError(1,T -2 0) 2 log 4R ,R Figure 4.Absolute error in log scale for T = (4R)7/4/5000, no filter (slope -1), infinite order filter (slope -3.4). 0.811.21.41.61.822.22.42.6-12 -10logError(k,T -8 -6), -4,2 -2 0 log 4R ,R k= 1 Figure 5.Absolute error in log scale for T = (4R)3/2/1000, AT,0,R,L(slope -3.1) and AT,1,R,L (slope -5.2), filter of infinite order. ESAIM: PROCEEDINGS99 0.40.60.811.21.41.61.82 -4.5 -4 -3.5 -3 -2.5Error(1,\infty -2 -1.5) -1 log R log ,R Figure 6.Error in log scale for (4.35) in function of the number of cells per dimension R \in
[16] without zero-order term, with and without filtering: Slope -1 in both cases. 0.40.60.811.21.41.61.82 -5.5 -5 -4.5 -4 -3.5 -3logError(1,T -2.5 -2 -1.5,R -1 -0.5 log R ) Figure 7.Error in log scale for (4.35) in function of the number of cells per dimension R \in
[17] with a zero-order term T = R/10, with and without filtering: Slopes -1 and -3. 0.40.50.60.70.80.911.11.21.31.41.51.61.7 -4 -3.8 -3.6gError(1,\infty -3.4 -3.2 -3) -2.8 -2.6 -2.4 -2.2 -2 log R lo ,R Figure 8.Error in log scale for (4.36) in function of the number of cells per dimension R \in [3, 42] without zero-order term and without filtering, for the two different coefficients of reference: Slope -1 and artificial super-linear convergence. 0.40.50.60.70.80.911.11.21.31.41.51.61.7 -7 -6 -5Error(1,T -4 -3) -2 -1 0 1 log R log ,R Figure 9.Error in log scale for (4.36) in function of the number of cells per dimension R \in ESAIM: PROCEEDINGS101 -LL Figure 10.Filter for the discrete stochastic example. 4.5.4. Discrete stochastic example To check the validity of (4.27), we consider the discrete case, in the form of equation (4.34). This time we assume the entries a of the symmetric matrix A to be random, independent from one another, and identically distributed (i. i. d.). More precisely, we consider a Bernoulli law (which amounts to topping a coin) with values {\(\alpha\)} > 0 et {\(\beta\)} \geq {\(\alpha\)} with probability 1/2. In this case the Dykhne formula holds for d = 2 so that the homogenized\sqrt{}{\(\alpha\)}{\(\beta\)}Id (see [30, Appendix A]). For the numerical tests, we take {\(\alpha\)} = 1, {\(\beta\)} = 9 (so matrix is explicit: Ahom= that Ahom= 3Id), and &#63729;T = L + 3 , &#63730;2(L)\sqrt{}L(L + 1) . &#63731;R =1 + 0.1ln The mask is the following: {\(\mu\)}L(x) = \tilde{}{\(\mu\)}L(x1)\tilde{}{\(\mu\)}L(x2), where x = (x1, x2) \in Z2, and \tilde{}{\(\mu\)}L: Z \to R+is given by \tilde{}{\(\mu\)}L(t) = {\(\gamma\)}LL/6 \leq |t| \leq L/2 : 3/2(1 - 2|t|/L), &#63729;L/2 \leq |t|:0, &#63730; &#63731;|t| \leq L/6 :1, and {\(\gamma\)}Lis such that\'{}Z\tilde{}{\(\mu\)}L(t)dt = 1, see Figure 10. For a uniform sampling of log L, L \in [20, 4000], weapproximate the expectation \hat{}2 ZdT,R(x)) {\(\cdot\)} A(x)(e1+ \nabla{\(\phi\)}1T,R(x)){\(\mu\)}L(x)dx - 3 (e1+ \nabla{\(\phi\)}1 by an empirical average over r(L) realizations, and define the error by Error(L) := 1r(L)\hat{}(e1+ \nabla{\(\phi\)}1,jT,R(x)) {\(\cdot\)} Aj(x)(e1+ \nabla{\(\phi\)}1,jT,R(x)){\(\mu\)}L(x)dx - 32, r(L)j=1Zd (4.37) where {{\(\phi\)}1,jT,R} are solutions of the regularized corrector equation on QRfor r(L) different realizations Ajof the coefficients A. In order to be sure that the error computed using the emprirical average ANT,R,Lof AT,R,Lover N independent realization is close to the true error, we have taken r(L) large enough so that the empirical variance of Ar(L)T,R,Lis much smaller than Error(L). 2.533.544.555.566.577.5-3.5 -3 -2.5 -2logError(L) -1.5 -1 -0.5 log(L2) Figure 11.Error (4.37) in log scale in function of L2(slope -1/2). The error is plotted in function of L2in logarithmic scale on Figure 11. The dots (which indicate calculations) are in very good agreement with the straight line of slope -1/2 corresponding to the decay provided by (4.32) (the exponent -1/2 is the central limit theorem scaling): |AT,R,L- Ahom|2 1/2(L2)-1/2. 4.6. Comments on the periodization method A very popular method in numerical homogenization is the so-called periodization method, where homogeneous Dirichlet boundary conditions are replaced by periodic boundary conditions. The numerical tests by Yue and E in [59] tend to show that periodic boundary conditions perform usually better than Dirichlet boundary conditions. The method described there reads: for all L > 0, an approximation of Ahomis given by, for all ij, \in {1, . . . , d}, ej{\(\cdot\)} AL,#ei:=(ej+ \nabla{\(\phi\)}jL,#) {\(\cdot\)} A(ei+ \nabla{\(\phi\)}iL,#), QL where for all k \in {1, . . . , d}, {\(\phi\)}kL,#is the unique weak periodic solution (with zero mean) to -\nabla {\(\cdot\)} A(ek+ \nabla{\(\phi\)}kL,#) = 0in QL. The aim of this subsection is to make two remarks on this numerical observation: ●in the case of random coefficients, the periodization method yields optimal results (without the zero order term regularization) provided the periodization is made at the level of the distribution of the random coefficients and not at the level of a realization, ●in the case of periodic (or quasiperiodic) coefficents or random coefficients with correlations, the periodization method can only be performed at the level of the realizations, and does not always reach optimal convergence rate (unlike the zero order regularization method). ESAIM: PROCEEDINGS103 Let us make these comments more precise by treating a couple of examples. In the case of a discrete elliptic equation with i. i. d. coefficients, it is shown in [35] that |AL,#- Ahom|2L-d, which is the optimal (central limit theorem) scaling. The analysis of [35] indeed yields more details: the error splits into a random error and a systematic error, |AL,#- Ahom|2= var [AL,#] + | AL,#- Ahom|2, whose scalings are for d \geq 2: var [AL,#]L-d, | AL,#- Ahom|2L-2dlndL. In particular, the crucial point in the proof of the estimate of the systematic error is the compatibility of the periodicity and the randomness – which can be easily coupled due to the product structure of the probability space in the i. i. d. case, see also [35] for more general statistics. It is not clear that in the presence of correlations, the right ”periodic approximation” can be constructed, and it is to be expected (although it is not proved) that if this is not done properly the systematic error may scale as | AL,#- Ahom| \sim L-1in any dimension (which is the scaling of a boundary effect, as in (4.4)). The analysis also shows that the scaling of the overall error is the same if the Lddegrees of freedom are distributed over several independent computations. For all N \in N let us define ANL,#as the empirical average of N independent realizations of AL,#, then we have |ANL,#- Ahom|2N-1L-d+ L-2dlndL, so that for the same convergence rate, it may be very advantageous in terms of computational cost to run several realizations on small domains rather than one realization on a large domain (this can indeed be made quantitative using this estimate). This result thus sets on mathematical grounds a method advocated by the mechanical engineering community, see [46]. Let us illustrate how to treat correlations correctly in the continuous case on the example of a homogeneous material perturbed by unit spherical inclusions whose centers are distributed according to a Poisson point process in Rd. The naive ”periodic” approach would be to consider a domain QL, generate a realization of a Poisson point process in QL, construct the associated diffusion matrix by adding spherical inclusions in QL centered at the Poisson points, and solve the corrector equation on QLwith periodic boundary conditions. An intuitive way to see that this periodization is not compatible with the statistics of the random diffusion matrix on Rdis that it creates new inclusion shapes in the picture (inclusions which intersect the boundary \partialQLare simply cut), so that stationarity is broken (the boundary of QLis clearly identified and there is no translation invariance any longer). The right way to proceed is to periodize the underlying point process, see QLas the torus TL, simulate a realization of a Poisson point process on TL(which is actually the same as on the cube), and then construct a diffusion matrix not on the cube but on the torus by adding spherical inclusions centered at the Poisson points. This way, spherical inclusions are not cut, and the statistics of the periodized random diffusion matrix and of the original random diffusion matrix are compatible in the sense that they are both translation invariant. The outputs of these two procedures are sketched on Figure 12. In other terms, there are two operators to be considered: the periodization operator {\(\pi\)}L(seen as acting both on point sets of Rdand on L\infty(Rd, A{\(\alpha\)}{\(\beta\)})) and the operator I which acts on point sets and takes values in L\infty(Rd, A{\(\alpha\)}{\(\beta\)}) (it adds inclusions centered at the points of the point set). What Figure 12 illustrate is that given a Poisson process P, one has {\(\pi\)}L&#9702; I(P) = I &#9702; {\(\pi\)}L(P) – the latter being preferable to obtain optimal error estimates. 000000111111000000111111 Figure 12.Naive periodization (left) versus compatible periodization (right) The worst case for the naive periodization method is in the periodic case itself. Let A be a Q-periodic matrix, and for all L \in R+let AL,#be defined for all {\(\xi\)} \in Rdby {\(\xi\)} {\(\cdot\)} AL,#{\(\xi\)} = inf{({\(\xi\)} + \nabla{\(\phi\)}) {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}), {\(\phi\)} \in H1#(QL)}. QL Then, if L \in N, AL,#= Ahom. Yet we have for all k \in N sup|AL,#- Ahom| \sim 1/k, k\leqL<k+1 which shows that generically the convergence rate of |AL,#- Ahom| to zero is not better than in the case when Dirichlet boundary conditions are used on \partialQL(yet, numerical tests show that the prefactor is better with periodic boundary conditons). In the case of a discrete elliptic equation on Zdwith random conductivities with finite correlation-length, numerical examples in [23] show that the naive periodization indeed yields a systematic error which scales as a surface effect. 5.Numerical homogenization with zero-order regularization In this section, we show how to combine the regularization method with the direct and dual approaches to numerical homogenization. As in the previous sections we begin with the analytical framework. There is a gap in terms of generality between the results we shall prove here and the results of Sections 2 and 3: our analysis does not apply to general H-converging sequences, and we shall limit ourselves to the case of symmetric diffusion matrices which are locally stationary and ergodic. ESAIM: PROCEEDINGS105 This section contains the main new results of this survey, and a more thorough analysis will be given in [32]. Our main achievement here is the introduction of a technique which is consistent in general, and allows to rid numerical homogenization methods of the resonance error completely in the periodic case. 5.1. Analytical framework Let A{\(\epsilon\)}: D {\(\times\)} &#8486; \to Msym{\(\alpha\)}{\(\beta\)}be a family of random symmetric diffusion coefficients parametrized by {\(\epsilon\)} > 0. We make the assumption that A{\(\epsilon\)}is locally stationary and ergodic (and assume some ”cross-regularity”, see below). Hypothesis 1.There exists a random Carath\'{}eodory function (that is continuous in the first variable and measurable in the second variable) \tilde{}A : D {\(\times\)} Rd{\(\times\)} &#8486; \to Msym{\(\alpha\)}{\(\beta\)}, and a constant {\(\kappa\)} > 0 such that ●For all x \in D, the random field \tilde{}A(x, {\(\cdot\)}, {\(\cdot\)}) is stationary on Rd{\(\times\)} &#8486;, and ergodic; ●(cross-regularity) \tilde{}A is {\(\kappa\)}-Lipschitz in the first variable: for all x, y \in D, for almost every z \in Rd, and for almost every realization {\(\omega\)} \in &#8486;, | \tilde{}A(x, z, {\(\omega\)}) - \tilde{}A(y, z, {\(\omega\)})| \leq {\(\kappa\)}|x - y|; ●For all x \in D, for and all {\(\epsilon\)} > 0, and for almost every realization {\(\omega\)} \in &#8486;, A{\(\epsilon\)}(x, {\(\omega\)}) := \tilde{}A(x,x{\(\epsilon\)}, {\(\omega\)}) almost surely. It is not difficult to prove that A{\(\epsilon\)}H-converges on D to some deterministic diffusion function Ahom: D \tosymalmost surely, which is {\(\kappa\)}-Lipschitz and characterized for all x \in D and {\(\xi\)} \in Rdby M{\(\alpha\)},{\(\beta\)}2/{\(\alpha\)} {\(\xi\)} {\(\cdot\)} Ahom(x){\(\xi\)} =({\(\xi\)} + \nabla{\(\phi\)}(x, 0, {\(\cdot\)})) {\(\cdot\)} \tilde{}A(x, 0, {\(\cdot\)})({\(\xi\)} + \nabla{\(\phi\)}(x, 0, {\(\cdot\)})) , where {\(\phi\)}(x, {\(\cdot\)}, {\(\cdot\)}) is the corrector in direction {\(\xi\)} associated with the stationary coefficients \tilde{}A(x, {\(\cdot\)}, {\(\cdot\)}) (x is treated as a parameter). The proof is standard: by H-compacteness, for almost every realization, A{\(\epsilon\)}H-converges up to extraction to some limit A\ast. Using the locality of H-convergence and the fact that \tilde{}A is uniformly Lipschitz in the first variable, one concludes that A\ast(x) = Ahom(x) almost surely. Note that if we weaken the cross-regularity assumption from Lipschitz continuity on D to continuity on D, the result holds true as well – although the proof of the homogenization result is much more subtle, see [11, Theorem 4.1.1]. The combination of the regularization approach with the analytical framework of Section 2 yields the following local approximation of Ahomin the symmetric case. Definition 8.Let {\(\delta\)} > 1, and let {\(\eta\)} : Q{\(\delta\)}\to [0, {\(\delta\)}] be a measurable mask of mass one, such that infQ{\(\eta\)} \geq 1/{\(\delta\)}, and for all {\(\rho\)} > 0, let {\(\eta\)}{\(\rho\)}: Q{\(\delta\)}{\(\rho\)}\to R+be the rescaled version {\(\eta\)}{\(\rho\)}: y \to {\(\rho\)}-d{\(\eta\)}(y/{\(\rho\)}) of {\(\eta\)}. For all {\(\rho\)} > 0, T > 0 and {\(\epsilon\)} > 0, we denote by A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}: D \to Mdthe function defined by: for all i, j \in {1, . . . , d} and for x \in D, {\(\xi\)} {\(\cdot\)} A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x){\(\xi\)} :=Q{\(\delta\)}{\(\rho\)}\capT-xD({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy,(5.1) \hat{} where v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, {\(\cdot\)}) is the unique weak solution in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD) to (T {\(\epsilon\)}2)-1v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y) - \nabla {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) = 0in Q{\(\delta\)}{\(\rho\)}\cap T-xD.(5.2) For an adaptation of this definition to the non-symmetric case (for which the convergence analysis is not yet clear beyond the periodic case) one may use the solution to the adjoint corrector equation, as in (4.7). The scaling of the zero-order term in (5.2) is natural, as can be seen on the periodic case: if A{\(\epsilon\)}(y) = A(y/{\(\epsilon\)}) with A periodic, the change of variable y = {\(\epsilon\)}z turns (5.2) into the equation T-1w(z) - \nabla {\(\cdot\)} A(x + z)({\(\xi\)} + \nablaw(z)) = 0in Q{\(\epsilon\)}-1{\(\delta\)}{\(\rho\)}\cap (T-xD)/{\(\epsilon\)}, which is of the same form as (4.5). Under Hypothesis 1, we have the following convergence result: Theorem 8.Let D be smooth. For all {\(\epsilon\)} > 0, {\(\rho\)} > 0, {\(\delta\)} > 1, and T > 0, let A{\(\epsilon\)}satisfy Hypothesis 1, and A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}’\geq {\(\alpha\)}’> 0 such that for all {\(\epsilon\)} > 0, {\(\rho\)} > 0 and{\(\delta\)}\in Msym{\(\alpha\)}’{\(\beta\)}’(D). In addition, for all x \in D, be as in Definition 8. Then there exist {\(\delta\)} > 1 small enough and {\(\beta\)} T > 0, AT,{\(\rho\)},{\(\epsilon\)} limT \to\infty,{\(\rho\)}\to0lim{\(\epsilon\)}\to0A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x)= Ahom(x),(5.3) almost surely. The limit in T in (5.3) is uniform in {\(\rho\)}. From this theorem, one may directly deduce the convergence of the regularizing method: Corollary 3.Under the assumption of Theorem 8, for all f \in H-1(D), the weak solution u{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}\in H10(D) to -\nabla {\(\cdot\)} A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}\nablau{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}= f satisfies limT \to\infty,{\(\rho\)}\to0lim{\(\epsilon\)}\to0u{\(\delta\)}{\(\rho\)},{\(\epsilon\)}- uhom H1(D)= 0(5.4) almost surely, where uhom\in H10(D) is the weak solution to -\nabla {\(\cdot\)} Ahom\nablauhom= f. As a consequence of H-convergence we also have that limT \to\infty,{\(\rho\)}\to0lim{\(\epsilon\)}\to0u{\(\delta\)}{\(\rho\)},{\(\epsilon\)}- u{\(\epsilon\)} L2(D)= 0 almost surely, where u{\(\epsilon\)}is the unique weak solution in H10(D) to -\nabla {\(\cdot\)} A{\(\epsilon\)}\nablau{\(\epsilon\)}= f Before we proceed with the proofs, let us give the associated corrector result. Definition 9.Let H > 0, IH\in N, and let {QH,i}i\in[[1,IH]]be a partition of D in disjoint subdomains of diameter2(D) associated with QH,i: for2(D) and H > 0, of order H. We define a family (MH) of approximations of the identity on L every w \in L MH(w) =w 1QH,i. IH i=1QH,i Let {\(\delta\)} > 1 be as in Theorem 8, and for all i \in [[1, IH]], set Q{\(\delta\)}H,i:= {x \in D | d(x, QH,i) < ({\(\delta\)} - 1)H}. ESAIM: PROCEEDINGS107 With the notation of Corollary 3, we define the numerical correctors {\(\gamma\)}{\(\delta\)},H,iT,{\(\rho\)},{\(\epsilon\)}associated with u{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}as the unique weak solution in H10(Q{\(\delta\)}H,i) to (T {\(\epsilon\)}2)-1{\(\gamma\)}{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)}- \nabla {\(\cdot\)} A{\(\epsilon\)}MH(\nablau{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)}= 0,(5.5) we set \nablau{\(\delta\)},H,iT,{\(\rho\)},{\(\epsilon\)}:= MH(\nablau{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)})|QH,i+ (\nabla{\(\gamma\)}{\(\delta\)},H,i{\(\rho\)},{\(\epsilon\)})|QH,i\in L2(QH,i) for all 1 \leq i \leq IH, and we define C{\(\delta\)},HT,{\(\rho\)},{\(\epsilon\)}=\nablau{\(\delta\)},H,iT,{\(\rho\)},{\(\epsilon\)}1QH,i. IH i=1 We then have the following corrector result: Theorem 9.Under the assumptions of Corollary 3, the corrector of Definition 9 satisfies limT \to\infty,{\(\rho\)},H\to0lim{\(\epsilon\)}\to0\nablau{\(\epsilon\)}- C{\(\delta\)},HT,{\(\rho\)},{\(\epsilon\)}L2(D)= 0(5.6) almost surely. We only prove Theorem 8, and quickly show how to adapt the proof of Theorem 3 to deal with the regularization. Proof of Theorem 8. The proof that A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}\in M{\(\alpha\)}’{\(\beta\)}’(D) for some {\(\delta\)} > 1 small enough is the same as for Theorem 2 since the argument only relies on Meyers’ estimates, which hold uniformly with respect to the zero-order term. We split the rest of the proof into three steps. Step 1. From locally ergodic to ergodic. In this step we introduce a proxy \tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) for A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) which is uniformly close to A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) in {\(\rho\)}, and for which we can apply our analysis of Section 4. For all x \in D, we define \tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) by: For all {\(\xi\)} \in Rd, {\(\xi\)} {\(\cdot\)} \tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x){\(\xi\)} :=({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}dy, \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD where \tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, {\(\cdot\)}) is the unique weak solution in H10(Q{\(\delta\)}{\(\rho\)}\cap T-xD) to (T {\(\epsilon\)}2)-1\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y) - \nabla {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) = 0in Q{\(\delta\)}{\(\rho\)}\cap T-xD. We shall prove that |A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) - \tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x)|{\(\rho\)}(5.7) uniformly in x, T , {\(\epsilon\)}, and the randomness. We first write the difference as: for all {\(\xi\)} \in Rdwith |{\(\xi\)}| = 1, {\(\xi\)} {\(\cdot\)} (A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) - \tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x)){\(\xi\)} =({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + y)({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD -({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD =({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD -({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD +({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} (A{\(\epsilon\)}(x + y) - \tilde{}A(x,x + y{\(\epsilon\)}))({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy. \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD Using thatx + y)| = | \tilde{}A(x + y,x + y{\(\epsilon\)}) - \tilde{}A(x,x + y{\(\epsilon\)})| \leq {\(\kappa\)}|y|(5.8) |A{\(\epsilon\)}(x + y) - \tilde{}A(x,{\(\epsilon\)} the third term of the r. h. s. is easily estimated: \hat{}x + y))({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy Q{\(\delta\)}{\(\rho\)}\capT-xD{\(\epsilon\)} ({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} (A{\(\epsilon\)}(x + y) - \tilde{}A(x, \leq {\(\kappa\)}{\(\rho\)}|{\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)|2{\(\eta\)}{\(\rho\)}(y)dy{\(\rho\)}(5.9) \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD by definition of {\(\eta\)}{\(\rho\)}and an elementary a priori estimate on v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(using that |{\(\xi\)}| = 1). We now turn to the first two terms, which we write in the form \hat{}x + y)({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy Q{\(\delta\)}{\(\rho\)}\capT-xDT(x, y)) {\(\cdot\)} \tilde{}A(x,{\(\epsilon\)} ({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)} -({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD =(\nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y) - \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})({\(\xi\)} + \nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD +({\(\xi\)} + \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) {\(\cdot\)} \tilde{}A(x,x + y{\(\epsilon\)})(\nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y) - \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)){\(\eta\)}{\(\rho\)}(y)dy \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xD so that using the same a priori estimate as above, it enough to prove that \hat{} Q{\(\delta\)}{\(\rho\)}\capT-xDT(x, y) - \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)|2{\(\eta\)}{\(\rho\)}(y)dy{\(\rho\)}2(5.10) |\nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)} in order to deduce (5.7) from (5.9). This estimate directly follows from the Lipschitz bound (5.8) and an a priori estimate on the equation satisfied by v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, {\(\cdot\)}) - \tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, {\(\cdot\)}): (T {\(\epsilon\)}2)-1(v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y) - \tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) - \nablay\tilde{}A(x,x + y{\(\epsilon\)})\nablay(v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y) - \nablay\tilde{}v{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y)) = \nablay{\(\cdot\)} (A{\(\epsilon\)}(x + y) - \tilde{}A(x,x + y{\(\epsilon\)}))\nablayv{\(\delta\)}{\(\rho\)},{\(\epsilon\)}T(x, y). ESAIM: PROCEEDINGS109 Step 2. Limit as {\(\epsilon\)} \to 0. The change of variables z = y/{\(\epsilon\)} allows us to interprete \tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) as the approximation ”AT,R,L(x)” of Ahom(x) defined in (4.7) (with \tilde{}A(x, {\(\cdot\)}) in place of A({\(\cdot\)})), with the parameters R = {\(\delta\)}{\(\rho\)}{\(\epsilon\)},L ={\(\rho\)}{\(\epsilon\)}. Following Subsection 4.2, we define AT,L(x) and AT(x) by (4.20) and (4.21) (with \tilde{}A(x, {\(\cdot\)}) in place of A({\(\cdot\)})). By the triangle inequality, |AT,R,L(x) - AT(x)| \leq |AT,R,L(x) - AT,L(x)| + |AT,L(x) - AT(x)|. By ergodicity, limL\to\infty|AT,L(x) - AT(x)| = 0 almost surely, so that the almost-sure convergence lim{\(\epsilon\)}\to0\tilde{}A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) =limL,R-L\to\inftyAT,R,L(x) = AT(x)(5.11) follows from (4.22) for all x \in D. Note that the limit does not depend on {\(\rho\)}, as expected by stationarity of\tilde{}A(x, y) in y. Step 3. Limit in T and conclusion. We finally use spectral analysis in the form of (4.24), which yields for all x \in D, limT \to\inftyAT(x) = Ahom(x). We are in position to conclude: the combination of (5.7), (5.11), and (4.24) yields for all x \in D lim suplim sup|A{\(\delta\)}T,{\(\rho\)},{\(\epsilon\)}(x) - Ahom(x)| = 0, T \to\infty,{\(\rho\)}\to0{\(\epsilon\)}\to0 as desired. The proof of Theorem 9 closely follows the proof of Theorem 3. There is yet one slight difference since the numerical corrector of Theorem 9 is constructed using the zero-order regularization of the corrector equation. This additionnal difficulty can be taken care of by using (5.10) and (4.25). We leave the details to the reader. Remark 6.In this section we have only considered the approximation AT,kof Ahomfor k = 1. We could of course use approximations of higher order k > 1. All the convergence results proved above hold true for higher order approximations as well. 5.2. Direct and dual approaches The direct approach associated with the analytical framework is a straightforward modification of the direct approach with oversampling, where the corrector equation is regularized by the zero-order term as in (5.2). Likewise, the Petrov-discontinuous Galerkin version of the dual approach of Subsection 2 can be adapted in a straightforward way by adding a zero-order term, as in (5.2). The extensions of Theorems 8 and 9 to cover the direct and dual approaches can be done as in Section 2, and we leave the details to reader as well. The more interesting question regards the practical use of the method, and therefore the choice of the strength of the zero-order term. The coefficient in front of the zero-order term in (5.2) is of the form (T {\(\epsilon\)}2)-1, and one needs some ”automatic” rule to choose T and {\(\epsilon\)}. The quantity {\(\epsilon\)} represents the local correlation length of A{\(\epsilon\)}(it is the period of A{\(\epsilon\)}in the periodic case, it is the correlation length in the stochastic case with finite correlation length, and in more general cases when there is a scale separation in A{\(\epsilon\)}this length can typically be identified using a wavelet transform). The choice of T is easier since it only depends on the ellipticity ratio {\(\beta\)}/{\(\alpha\)} (recallsym) which drives the value of the coefficient c in the argument of the exponential decay in (4.22). that A{\(\epsilon\)}\in M{\(\alpha\)}{\(\beta\)}sym. This choice of {\(\epsilon\)} and T Hence the same value of T can be used for all the coefficients A{\(\epsilon\)}of some class M{\(\alpha\)}{\(\beta\)} in actual computations may be a subtle issue. To complete the convergence analysis, we provide in the following subsection the numerical analysis of the locally periodic case. 5.3. Numerical analysis of the locally periodic case In this subsection we consider the case of a diffusion function \tilde{}A : D {\(\times\)} Rd\to Msym{\(\alpha\)}{\(\beta\)}such that for all x \in D, \tilde{}A(x, {\(\cdot\)}) is measurable and Q-periodic, and such that \tilde{}A is uniformly Lipschitz in the first variable. Defined for all {\(\epsilon\)} > 0 and x \in D by A{\(\epsilon\)}(x) := \tilde{}A(x, x/{\(\epsilon\)}), A{\(\epsilon\)}satisfies Hypothesis 1, so that the regularized versions of the direct and dual approaches converge. This qualitative convergence result can be turned quantitative, and the combination of the analysis leading to Theorem 6 with [1, 17] and [44] yields for the direct approach (with obvious notation): u{\(\delta\)},H,hT,H,{\(\epsilon\)}- uhom H1(D)H +{\(\epsilon\)}H4-+h{\(\epsilon\)}2,(5.12) \nablau{\(\epsilon\)}- C{\(\delta\)},H,hT,H,{\(\epsilon\)} L2(D)H +{\(\epsilon\)}H2-+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)},(5.13) and for the dual approach (with obvious notation) u{\(\delta\)}T,H,{\(\epsilon\)},h- u{\(\epsilon\)} H1(D)H +{\(\epsilon\)}H2-+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)},(5.14) where 4-(resp. 2-) stands for any exponent strictly less than 4 (resp. 2). These scalings are obtained by taking T :=H{\(\epsilon\)}2ln-2H{\(\epsilon\)} in the estimates. These rates are to be compared to (3.4), (3.5), and (3.7). In the three cases, the regularization method yields better convergence rates. The effect of the regularization on the corrector estimates (5.13) and (5.14) are much less impressive than for (5.12). One can do better. 5.4. Richardson extrapolation for the numerical corrector In this subsection we show how to take advantage of the Richardson extrapolation technique to reduce the resonance error on the corrector itself. Let us go back to the spectral representation (4.18) of the homogenized coefficients in the stationary ergodic case. We have seen in (4.26) that {\(\xi\)} {\(\cdot\)} (AT- Ahom){\(\xi\)} =(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}T) {\(\cdot\)} A(\nabla{\(\phi\)}{\(\xi\)}T- \nabla{\(\phi\)}{\(\xi\)}) , where {\(\phi\)}{\(\xi\)}and {\(\phi\)}{\(\xi\)}Tare the corrector and regularized corrector in direction {\(\xi\)}. In Subsection 4.4, we have also seenthat the approximation ATof Ahomcould be made more accurate by using extrapolation. The same strategy holds for the corrector as well. ESAIM: PROCEEDINGS111 Let A be a symmetric stationary random field, and let {\(\xi\)} \in Rdwith |{\(\xi\)}| = 1 be fixed. We present the first step of the extrapolation, and we define {\(\phi\)}2,Tas: {\(\phi\)}2,T:= 2{\(\phi\)}2T- {\(\phi\)}T, where {\(\phi\)}2Tand {\(\phi\)}Tare the unique stationary solutions with vanishing expectation to {\(\tau\)}-1{\(\phi\)}{\(\tau\)}- \nabla {\(\cdot\)} A({\(\xi\)} + \nabla{\(\phi\)}{\(\tau\)}) = 0, for {\(\tau\)} = 2T and {\(\tau\)} = T , respectively. Recall that edis the projection of the spectral measure of L = -D {\(\cdot\)} AD on the local drift d = -D {\(\cdot\)} A{\(\xi\)}. The spectral representation of the operators (T-1+ L)-1, ((2T )-1+ L)-1, L-1 and L being {\(\lambda\)} \to (T-1+ {\(\lambda\)})-1, {\(\lambda\)} \to ((2T )-1+ {\(\lambda\)})-1, {\(\lambda\)} \to {\(\lambda\)}-1, and {\(\lambda\)} \to {\(\lambda\)}, we have \nabla({\(\phi\)} - {\(\phi\)}2,T) {\(\cdot\)} A\nabla({\(\phi\)} - {\(\phi\)}2,T)={\(\lambda\)}1{\(\lambda\)}- 21T-1/2 + {\(\lambda\)}+1T-1+ {\(\lambda\)}2ded({\(\lambda\)}) \hat{}\infty 0 =T-44\hat{}\infty01{\(\lambda\)}(T-1+ {\(\lambda\)})2(T-1/2 + {\(\lambda\)})2ded({\(\lambda\)}). This implies the consistency of the method for general stationary ergodic random fields A. Indeed, letting T \to \infty and using the Lebesgue dominated convergence theorem, as for (4.24), we have limT \to\infty|\nabla({\(\phi\)} - {\(\phi\)}2,T)|2= 0. In addition, provided {\(\lambda\)} \to {\(\lambda\)}-5is dedintegrable, \nabla({\(\phi\)} - {\(\phi\)}2,T)L2(Q)T-2. This is in particular the case when A is periodic, since then the operator L has a spectral gap (there is a Poincar\'{}e inequality on H1#(Q)), so that the spectral integral is bounded independently of T (the integration interval is-1is ded({\(\lambda\)}) on R+). We have therefore gained one power of T-1with respect indeed isolated from zero, and {\(\lambda\)} to our previous estimate \nabla({\(\phi\)} - {\(\phi\)}T)L2(Q)T-1. We can iterate this procedure and obtain any convergence rate in T for the symmetric periodic case. This approach also generalizes to the nonsymmetric periodic case using only PDE arguments (and not spectral theory), and we refer the reader to [32] for details. Let use this new approximation of the corrector in the direct approach. As in Definition 4 we let IH\in N,{\(\delta\)}:= and {QH,i}i\in[[1,IH]]be a partition of D in disjoint subdomains of diameter of order H. We further set QH,i {x \in D | d(x, QH,i) < ({\(\delta\)} - 1)H}, which is an enlarged version of QH,i. For all h > 0 and i \in [[1, IH]] we let V{\(\delta\)}H,i,h be a Galerkin subspace of H10(Q{\(\delta\)}H,i). For all {\(\tau\)} > 0, we define the numerical correctors {\(\gamma\)}{\(\delta\)},H,h,i{\(\tau\)},{\(\rho\)},{\(\epsilon\)}associated with u{\(\delta\)},H,hT,{\(\rho\)},{\(\epsilon\)}as the unique weak solution in V{\(\delta\)}H,i,hto ({\(\tau\)} {\(\epsilon\)}2)-1{\(\gamma\)}{\(\delta\)},H,h,i{\(\tau\)},{\(\rho\)},{\(\epsilon\)}- \nabla {\(\cdot\)} A{\(\epsilon\)}MH(\nablau{\(\delta\)},H,hT,{\(\rho\)},{\(\epsilon\)}) + \nabla{\(\gamma\)}{\(\delta\)},H,h,i{\(\tau\)},{\(\rho\)},{\(\epsilon\)}= 0, define {\(\gamma\)}{\(\delta\)},H,h,i2,T,{\(\rho\)},{\(\epsilon\)}:= 2{\(\gamma\)}{\(\delta\)},H,h,i2T,{\(\rho\)},{\(\epsilon\)}- {\(\gamma\)}{\(\delta\)},H,h,iT,{\(\rho\)},{\(\epsilon\)}, set \nablav{\(\delta\)},H,h,i2,T,{\(\rho\)},{\(\epsilon\)}:= MH(\nablau{\(\delta\)},H,hT,{\(\rho\)},{\(\epsilon\)})|QH,i+ (\nabla{\(\gamma\)}{\(\delta\)},H,h,i2,T,{\(\rho\)},{\(\epsilon\)})|QH,i for all 1 \leq i \leq IH, and finally consider C{\(\delta\)},H,h2,T,H,{\(\epsilon\)}:=\nablav{\(\delta\)},H,h,i2,T,H,{\(\epsilon\)}1QH,i. IH In the periodic case, we then have: u{\(\delta\)},H,hT,H,{\(\epsilon\)}- uhom H1(D)H +{\(\epsilon\)}H4-+h{\(\epsilon\)}2, \nablau{\(\epsilon\)}- C{\(\delta\)},H,h2,T,H,{\(\epsilon\)} L2(D)H +{\(\epsilon\)}H4-+h{\(\epsilon\)}+\sqrt{}{\(\epsilon\)}, provided T :=H{\(\epsilon\)}2ln-2H{\(\epsilon\)}. Doing so, we’ve upgraded the convergence rate for the corrector. These theoretical results are quite spectacular, and one readily convinces oneself that any finite order of convergence can be reached by iterating the Richardson procedure: the resonance error has disappeared. 6.Other approaches and perspectives 6.1. Other approaches In this survey we have essentially focused on the so-called HMM (direct approach) and so-called MsFEM (dual approach), which are two popular numerical homogenization methods. The convergence analysis of these methods can be made using the same analytical framework – the difference between the two approaches arising during the discretization. Another more global approach has been introduced by Schwab and Hoang in [41], inspired by ●two-scale convergence [4] and periodic unfolding [14], ●sparse tensor-products approximation [40]. The consistency of this approach can be proved using our analytical framework as well. With the notation of Section 2.1, their starting point is indeed the fact that the equation for u{\(\rho\)},{\(\epsilon\)}takes the equivalent form (provided A{\(\epsilon\)}is extended on a neigborhood of D): Find u{\(\rho\)},{\(\epsilon\)}\in H10(D) and U{\(\rho\)},{\(\epsilon\)}\in L2(D, H10(Q)) such that for all v \in H10(D)2(D, H10(Q)), and all V \in L \hat{}\hat{}\hat{} DQD (\nablaxv(x) + \nablayV (x, y)) {\(\cdot\)} A{\(\epsilon\)}(x + {\(\rho\)}y)(\nablaxu{\(\rho\)},{\(\epsilon\)}(x) + \nablayU{\(\rho\)},{\(\epsilon\)}(x, y))dxdy =f (x)v(x)dx. This point of view allows them to look for a efficient approximation of the space L2(D, H10(Q)), and implement sparse-tensor product strategies. Another approach advocated by Owhadi and Zhang in [53], whose origin is not the homogenization theory properly speaking, is based on the notion of harmonic coordinates. The starting point of their work is that although solutions u \in H10(D) to the equation -\nabla {\(\cdot\)} A\nablau = f are usually not more regular than H1in the Euclidian coordinates, they are H2in the so-called harmonic coordinates provided f \in L2(D). These harmonic coordinates are defined as x \to x + {\(\Gamma\)}(x) {\(\cdot\)} x, where {\(\Gamma\)} =10(D) to ({\(\gamma\)}1, . . . , {\(\gamma\)}d) and {\(\gamma\)}kis the unique weak solution in H -\nabla {\(\cdot\)} A(ek+ \nabla{\(\gamma\)}k) = 0 for all k \in {1, . . . , k}. The issue is then to compute (approximations of) these harmonic coordinates. An interesting link between harmonic coordinates and the MsFEM has been made by Allaire and Brizzi in [5], who consider as multiscale finite element space the functions x \to vH(x + {\(\Gamma\)}H(x)), where {\(\Gamma\)}His a local approximation ESAIM: PROCEEDINGS113 of the harmonic coordinates {\(\Gamma\)} on each simplex of the macroscopic tesselation of D. The combination of harmonic coordinates with the regularization method was recently studied by Owhadi and Zhang in [54]. 6.2. Beyond the linear case The reader may wonder what remains of all this when we replace the linear equation by a nonlinear equation, say elliptic monotone. Then part of the results translates quite naturally to the nonlinear case, namely the analytical framework (see [20, 21] for the description of the MsFEM method for nonlinear problems and error estimates in the periodic case, and [27,28] for the extension of the analytical framework of Section 2 to monotone elliptic operators and nonconvex integral functionals). Yet most of the quantitative results break down. It is also not clear how the regularization method behaves in that case. In addition, from a practical point of view, things get much worse. In the linear case, one could have the excuse that the multiscale finite element basis or the local homogenized coefficients and local correctors only had to be computed once, and could be used to solve the equation for a wide variety of boundary conditions and right hand sides. This is definitely an advantage when dealing with optimal control for instance. In the nonlinear case, this does not hold any longer, and the local basis or local values of the operator have to be computed on the fly. This is not doable in practice. For nonlinear problems, even periodic homogenization is not that easy since the object to reconstruct is a homogenized nonlinear map (and not a single matrix). The same fact holds for stochastic homogenization (with the additional difficulty that the corrector equation is posed on the whole space). A possible strategy is to compute approximations of the nonlinear map at some sampling points, and try to reconstruct an analytical approximation of this map by solving an inverse problem. This analytical approximation should at least satisfy similar properties as the homogenized nonlinear map (say convexity, isotropy, etc.). This strategy has been used in [16] to approximate the homogenized energy density associated with a discrete model for rubber introduced in [33] and whose homogenization limit has been etablished in [3]. There, the homogenized energy density is known to be quasiconvex, isotropic, and minimal at identity. An important issue is then to identify a suitable (explicit) subset of such functions in order to approximate the inverse problem. The constraint being nonlinear and the error landscape being quite rough, deterministic optimization methods usually fail to converge. In [16], the use of an evolutionary algorithm has been quite efficient. An intermediate case between nonlinear problems and the linear equations considered throughout this survey is the case of parametrized linear equations. Such an example naturally appears in [31], where a coupled system of elliptic/parabolic equations is studied (for application to nuclear waste storage). Indeed, although the two equations are linear and the coefficients periodic, the nonlinear coupling condition makes the homogenized coefficient of the parabolic equation depend locally on the gradient of the solution of the homogenized elliptic equation. We are thus lead to solve a family of corrector equations of the form: Find {\(\phi\)}{\(\zeta\)}\in H1#(Q) solution to -\nabla {\(\cdot\)} A({\(\zeta\)})({\(\xi\)} + \nabla{\(\phi\)}{\(\zeta\)}) = 0in Q, parametrized by {\(\zeta\)} \in Rd. This is an ideal framework for the reduced basis method [48], whose paradigm is that{\(\zeta\)}, {\(\zeta\)} \in Rd may be a thin subset of H1#(Q), so that it could be well approximated by a suitable low1(Q) (see [9, 15] for the case when A({\(\zeta\)}) = A0+ {\(\zeta\)}A1, with {\(\zeta\)} in some compact set of the family {{\(\phi\)} dimensional subspaces of H# R). This method has already been used in the context of homogenization by Boyaval in [13]. As usual with the reduced basis method, most of the work has to be done on the choice of the estimator, and on the fast-assembly method. In [31] we have used an estimator which is specific to homogenization problems, and appealed to the fast Fourier transform to assemble rigidity matrices efficiently. 6.3. What next ? In the nonlinear case, much remains to be done. The approach which consists in constructing analytical proxies for the homogenized operators raise quite either unaddressed or unsolved challenging questions in approximation theory. Even in the linear case, things are far from complete yet. In particular, as emphasized by Nolen, Papanicolaou, and Pironneau in [52], the numerical homogenization Grail is an adaptive method which would take the best of all worlds: use efficient domain decomposition techniques when needed, and exploit homogenization when possible. References}
[18] A. Abdulle. On a priori error analysis of fully discrete heterogeneous multiscale FEM. Multiscale Model. Simul., 4:447–459, 2005. · Zbl 1092.65093
[19] R. Alicandro and M. Cicalese. A general integral representation result for the continuum limits of discrete energies with superlinear growth. SIAM J. Math. Anal., 36(1):1–37, 2004. · Zbl 1070.49009
[20] R. Alicandro, M. Cicalese, and A. Gloria. Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity. Arch. Ration. Mech. Anal., 200(3):881–943, 2011. · Zbl 1294.74056
[21] G. Allaire. Homogenization and two-scale convergence. SIAM J. Math. Anal., 23:1482–1518, 1992. · Zbl 0770.35005
[22] G. Allaire and R. Brizzi. A multiscale finite element method for numerical homogenization. Multiscale Model. Simul., 4:790–812, 2005. · Zbl 1093.35007
[23] T. Arbogast. Numerical subgrid upscaling of two-phase flow in porous media. In Numerical treatment of multiphase flows in porous media (Beijing, 1999), volume 552 of Lecture Notes in Phys., pages 35–49. Springer, Berlin, 2000. · Zbl 1072.76560
[24] A. Bensoussan, J.-L. Lions, and G. Papanicolaou. Boundary layer analysis in homogenization of diffusion equations with Dirichlet conditions in the half space. Wiley, 1976. · Zbl 0411.60078
[25] A. Bensoussan, J.-L. Lions, and G. Papanicolaou. Asymptotic analysis for periodic structures, volume 5 of Studies in Mathematics and its Applications. North-Holland Publishing Co., Amsterdam, 1978. · Zbl 0404.35001
[26] P. Binev, A. Cohen, W. Dahmen, R. DeVore, G. Petrova, and P. Wojtaszczyk. Convergence rates for greedy algorithms in reduced basis methods. SIAM J. Math. Anal., 43(3):1457–1472, 2011. · Zbl 1229.65193
[27] X. Blanc and C. Le Bris. Improving on computation of homogenized coefficients in the periodic and quasi-periodic settings. Netw. Heterog. Media, 5(1):1–29, 2010. · Zbl 1262.65115
[28] A. Bourgeat, A. Mikeli\'{}c, and S. Wright. Stochastic two-scale convergence in the mean and applications. J. Reine Angew. Math., 456:19–51, 1994. · Zbl 0808.60056
[29] A. Bourgeat and A. Piatnitski. Approximations of effective coefficients in stochastic homogenization. Ann. I. H. Poincar\'{}e Probab. Stat., 40(2):153–165, 2004. · Zbl 1058.35023
[30] S. Boyaval. Reduced-basis approach for homogenization beyond the periodic setting. Multiscale Model. Simul., 7(1):466–494, 2008. · Zbl 1156.74358
[31] D. Cioranescu, A. Damlamian, and G. Griso. Periodic unfolding and homogenization. C. R. Math. Acad. Sci. Paris, 335(1):99– 104, 2002. · Zbl 1001.49016
[32] A. Cohen, R. Devore, and C. Schwab. Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE’s. Anal. Appl. (Singap.), 9(1):11–47, 2011. · Zbl 1219.35379
[33] M. De Buhan, A. Gloria, P. Le Tallec, and M. Vidrascu. Reconstruction of a constitutive law for rubber from in silico experiments. In preparation.
[34] W. E, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden. Heterogeneous multiscale methods: A review. Commun. Comput. Phys., 2:367–450, 2007. · Zbl 1164.65496
[35] W. E, P.B. Ming, and P.W. Zhang. Analysis of the heterogeneous multiscale method for elliptic homogenization problems. J. Amer. Math. Soc., 18:121–156, 2005. · Zbl 1060.65118
[36] Weinan E. Principles of multiscale modeling. Cambridge University Press, Cambridge, 2011. · Zbl 1238.00010
[37] Y. Efendiev and T. Y. Hou. Multiscale finite element methods, volume 4 of Surveys and Tutorials in the Applied Mathematical Sciences. Springer, New York, 2009. Theory and applications. · Zbl 1163.65080
[38] Y.R. Efendiev, T.Y. Hou, and V. Ginting. Multiscale finite element methods for nonlinear problems and their applications. Com. Math. Sc., 2(4):553–589, 2004. · Zbl 1083.65105
[39] Y.R. Efendiev, T.Y. Hou, and X.H. Wu. Convergence of a nonconforming multiscale finite element method. SIAM J. Num. Anal., 37:888–910, 2000. · Zbl 0951.65105
[40] A.-C. Egloffe, A. Gloria, J.-C. Mourrat, and T. N. Nguyen. Theoretical and numerical investigation of convergence rates in stochastic homogenization: corrector equation and random walk in random environment. In prepration. · Zbl 1315.60116
[41] L. C. Evans and R. F. Gariepy. Measure theory and fine properties of functions. Studies in Advanced Mathematics. CRC ESAIM: PROCEEDINGS115 · Zbl 0804.28001
[42] F. F\'{}eyel. Multiscale F E2elastoviscoplastic analysis of composite structures. Comp. Mat. Sci., 16:344–354, 1999.
[43] FreeFEM. http://www.freefem.org/.
[44] A. Gloria. An analytical framework for the numerical homogenization of monotone elliptic operators and quasiconvex energies. Multiscale Model. Simul., 5(3):996–1043, 2006. · Zbl 1119.74038
[45] A. Gloria. An analytical framework for numerical homogenization - Part II: windowing and oversampling. Multiscale Model. Simul., 7(1):275–293, 2008. · Zbl 1156.74362
[46] A. Gloria. Reduction of the resonance error - Part 1: Approximation of homogenized coefficients. Math. Models Methods Appl. Sci., 21(8):1601–1630, 2011. · Zbl 1233.35016
[47] A. Gloria. Numerical approximation of effective coefficients in stochastic homogenization of discrete elliptic equations. M2AN Math. Model. Numer. Anal., 46(1):1–38, 2012. · Zbl 1282.35038
[48] A. Gloria, T. Goudon, and S. Krell. Numerical homogenization of a nonlinearly coupled elliptic-parabolic system, reduced basis method, and application to nuclear waste storage. 2012. Preprint, http://hal.archives-ouvertes.fr/hal-00674519. · Zbl 1330.35027
[49] A. Gloria and Z. Habibi. Reduction of the resonance error - Part 2: Approximation of correctors and spectral theory. In preparation. · Zbl 1335.65086
[50] A. Gloria, P. Le Tallec, and M. Vidrascu. Comparison of network-based models for rubber. 2012. Preprint, http://hal.archivesouvertes.fr/hal-00673406.
[51] A. Gloria and J.-C. Mourrat. Spectral measure and approximation of homogenized coefficients. Probab. Theory. Relat. Fields. DOI 10.1007/s00440-011-0370-7. · Zbl 1264.35021
[52] A. Gloria, S. Neukamm, and F. Otto. Approximation of effective coefficients by periodization in stochastic homogenization. In preparation. · Zbl 1307.35029
[53] A. Gloria, S. Neukamm, and F. Otto. Quantification of ergodicity in stochastic homogenization: optimal bounds via spectral gap on Glauber dynamics. In preparation. · Zbl 1314.39020
[54] A. Gloria and F. Otto. Optimal quantitative estimates in stochastic homogenization of linear elliptic equations. In preparation. · Zbl 1387.35031
[55] A. Gloria and F. Otto. An optimal variance estimate in stochastic homogenization of discrete elliptic equations. Ann. Probab., 39(3):779–856, 2011. · Zbl 1215.35025
[56] A. Gloria and F. Otto. An optimal error estimate in stochastic homogenization of discrete elliptic equations. Ann. Appl. Probab., 22(1):1–28, 2012. · Zbl 1387.35031
[57] M. Griebel and S. Knapek. Optimized tensor-product approximation spaces. Constr. Approx., 16(4):525–540, 2000. · Zbl 0969.65107
[58] V.H. Hoang and Ch. Schwab. High-dimensional finite elements for elliptic problems with multiple scales. Multiscale Model. Simul., 3(1):168–194, 2005. · Zbl 1074.65135
[59] T.Y. Hou and X.H Wu. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys., 134:169–189, 1997. · Zbl 0880.73065
[60] T.Y. Hou, X.H. Wu, and Z.Q. Cai. Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients. Math. Comput., 68:913–943, 1999. · Zbl 0922.65071
[61] T.Y. Hou, X.H. Wu, and Y. Zhang. Removing the cell resonance error in the multiscale finite element method via a PetrovGalerkin formulation. Comm. in Math. Sci., 2(2):185–205, 2004. · Zbl 1085.65109
[62] V.V. Jikov, S.M. Kozlov, and O.A. Oleinik. Homogenization of Differential Operators and Integral Functionals. SpringerVerlag, Berlin, 1994.
[63] T. Kanit, S. Forest, I. Galliet, V. Mounoury, and D. Jeulin. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int. J. Sol. Struct., 40:3647–3679, 2003. · Zbl 1038.74605
[64] S.M. Kozlov. The averaging of random operators. Mat. Sb. (N.S.), 109(151)(2):188–202, 327, 1979.
[65] Y. Maday, A. T. Patera, and G. Turinici. Global a priori convergence theory for reduced-basis approximations of singleparameter symmetric coercive elliptic partial differential equations. C. R. Math. Acad. Sci. Paris, 335(3):289–294, 2002. · Zbl 1009.65066
[66] N. Meyers. An Lp-estimate for the gradient of solutions of second order elliptic divergence equations. Ann. Scuola Norm. Sup. Pisa Cl. Sci. (3), 17(3):189–206, 1963. · Zbl 0127.31904
[67] F. Murat. H-convergence. S\'{}eminaire d’Analyse fonctionnelle et num\'{}erique, Univ. Alger, multigraphi\'{}e, 1978.
[68] F. Murat and L. Tartar. H-convergence. In A.V. Cherkaev and R.V. Kohn, editors, Topics in the Mathematical Modelling of Composites Materials, volume 31 of Progress in nonlinear differential equations and their applications, pages 21–44. Birkh\"{}auser, 1997. · Zbl 0920.35019
[69] J. Nolen, G. Papanicolaou, and O. Pironneau. A framework for adaptive multiscale methods for elliptic problems. Multiscale Model. Simul., 7(1):171–196, 2008. · Zbl 1160.65342
[70] H. Owhadi and L. Zhang. Metric-based upscaling. Comm. Pure Appl. Math., 60(5):675–723, 2007. · Zbl 1190.35070
[71] H. Owhadi and L. Zhang. Localized bases for finite dimensional homogenization approximations with non-separated scales and high-contrast. Multiscale Model. Simul., 9(4):1373–1398, 2011. · Zbl 1244.65140
[72] G.C. Papanicolaou and S.R.S. Varadhan. Boundary value problems with rapidly oscillating random coefficients. In Random fields, Vol. I, II (Esztergom, 1979), volume 27 of Colloq. Math. Soc. J\'{}anos Bolyai, pages 835–873. North-Holland, Amsterdam, 1981.
[73] E. M. Stein. Singular integrals and differentiability properties of functions, volume 30 of Princeton Mathematical Series. Princeton University Press, Princeton, NJ, 1970. · Zbl 0207.13501
[74] L. Tartar. Cours Peccot du Coll‘ege de France. 1977.
[75] M. Vogelius. A homogenization result for planar, polygonal networks. RAIRO Mod\'{}el. Math. Anal. Num\'{}er., 25(4):483–514, 1991. · Zbl 0737.35126
[76] X. Yue and W. E. The local microscale problem in the multiscale modeling of strongly heterogeneous media: effects of boundary conditions and cell size. J. Comput. Phys., 222(2):556–572, 2007. · Zbl 1158.74541
[77] V. V. Yurinski\breve{}ı. Averaging of symmetric diffusion in random medium. Sibirskii Matematicheskii Zhurnal, 27(4):167–180, 1986.
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. In some cases that data have been complemented/enhanced by data from zbMATH Open. This attempts to reflect the references listed in the original paper as accurately as possible without claiming completeness or a perfect matching.