Localized bases for finite dimensional homogenization approximations with non-separated scales and high-contrast
aa r X i v : . [ m a t h . NA ] A ug Localized bases for finite dimensional homogenizationapproximations with non-separated scales and high-contrast.
Houman Owhadi ∗ , Lei Zhang † June 30, 2018
Abstract
We construct finite-dimensional approximations of solution spaces of divergenceform operators with L ∞ -coefficients. Our method does not rely on concepts of er-godicity or scale-separation, but on the property that the solution space of theseoperators is compactly embedded in H if source terms are in the unit ball of L in-stead of the unit ball of H − . Approximation spaces are generated by solving ellipticPDEs on localized sub-domains with source terms corresponding to approximationbases for H . The H -error estimates show that O ( h − d )-dimensional spaces withbasis elements localized to sub-domains of diameter O ( h α ln h ) (with α ∈ [ , O ( h − α ) accuracy for elliptic, parabolic and hyperbolic problems. Forhigh-contrast media, the accuracy of the method is preserved provided that local-ized sub-domains contain buffer zones of width O ( h α ln h ) where the contrast of themedium remains bounded. The proposed method can naturally be generalized tovectorial equations (such as elasto-dynamics). Consider the partial differential equation ( − div (cid:16) a ( x ) ∇ u ( x ) (cid:17) = g ( x ) x ∈ Ω; g ∈ L (Ω) , a ( x ) = { a ij ∈ L ∞ (Ω) } u = 0 on ∂ Ω , (1.1)where Ω is a bounded subset of R d with a smooth boundary (e.g., C ) and a is symmetricand uniformly elliptic on Ω. It follows that the eigenvalues of a are uniformly boundedfrom below and above by two strictly positive constants, denoted by λ min ( a ) and λ max ( a ).Precisely, for all ξ ∈ R d and x ∈ Ω, λ min ( a ) | ξ | ≤ ξ T a ( x ) ξ ≤ λ max ( a ) | ξ | . (1.2) ∗ Corresponding author. California Institute of Technology, Computing & Mathematical Sciences ,MC 217-50 Pasadena, CA 91125, [email protected] † University of Oxford, Mathematical Institute
1n this paper, we are interested in the homogenization of (1.1) (and its parabolicand hyperbolic analogues in Sections 4 and 5), but not in the classical sense, i.e., thatof asymptotic analysis [9] or that of G or H -convergence ([47], [57, 32]) in which oneconsiders a sequence of operators − div( a ǫ ∇ ) and seeks to characterize limits of solution.We are interested in the homogenization of (1.1) in the sense of “numerical homogeniza-tion,” i.e., that of the approximation of the solution space of (1.1) by a finite-dimensionalspace.This approximation is not based on concepts of scale separation and/or of ergodicitybut on compactness properties, i.e., the fact that the unit ball of the solution spaceis compactly embedded into H (Ω) if source terms ( g ) are integrable enough. Thishigher integrability condition on g is necessary because if g spans H − (Ω), then thesolution space of (1.1) is H (Ω) (and it is not possible to obtain a finite dimensionalapproximation subspace of H (Ω) with arbitrary accuracy in H -norm). However, if g spans the unit ball of L (Ω), then the solution space of (1.1) shrinks to a compactsubset of H (Ω) that can be approximated to an arbitrary accuracy in H -norm byfinite-dimensional spaces [10] (observe that if a = I d , then the solution space is a closedbounded subset of H ∩ H (Ω), which is known to be compactly embedded into H (Ω)).The identification of localized bases spanning accurate approximation spaces relieson a transfer property obtained in [10]. For the sake of completeness, we will give ashort reminder of that property in Section 2. In Section 3, we will construct localizedapproximation bases with rigorous error estimates (under no further assumptions on a than those given above). In Sub-section 3.4, we will also address the high-contrastscenario in which λ max ( a ) is allowed to be large. In Sections 4 and 5, we will show thatthe approximation spaces obtained by solving localized elliptic PDEs remain accurate forparabolic and hyperbolic time-dependent problems. We refer to Section 6 for numericalexperiments. We refer to Section B of the Appendix for further discussion and a proofof the strong compactness of the solution space when the range of g is a closed boundedsubset of H − ν (Ω) with ν < Recall that the key element in G and H convergence is a notion of “compactness bycompensation” combined with convergence of fluxes. Here, the notion of compactness iscombined with a flux-norm introduced in [10]. The flux-norm.
We will now give a short reminder on the flux-norm and its properties.
Definition 2.1.
For k ∈ ( L (Ω)) d , denote by k pot the potential portion of the Weyl-Helmholtz decomposition of k . Recall that k pot is the orthogonal projection of k onto {∇ f : f ∈ H Ω) } in ( L (Ω)) d . 2 efinition 2.2. For ψ ∈ H (Ω), define k ψ k a -flux := k ( a ∇ ψ ) pot k ( L (Ω)) d . (2.1)We call k ψ k a -flux the flux-norm of Ψ.The following proposition shows that the flux-norm is equivalent to the energy normif λ min ( a ) > λ min ( a ) < ∞ . Proposition 2.1. [Proposition 2.1 of [10]] k . k a -flux is a norm on H (Ω) . Furthermore,for all ψ ∈ H (Ω) λ min ( a ) k∇ ψ k ( L (Ω)) d ≤ k ψ k a -flux ≤ λ max ( a ) k∇ ψ k ( L (Ω)) d . (2.2) Motivations behind the flux-norm:
There are three main motivations behind theintroduction of the flux norm. • The flux-norm allows to obtain approximation error estimates independent fromboth the minimum and maximum eigenvalues of a . In fact, the flux-norm of thesolution of (1.1) is independent from a altogether since k u k a -flux = k∇ ∆ − g k ( L (Ω)) d . (2.3) • The ( · ) pot in the a -flux-norm is explained by the fact that in practice, we are inter-ested in fluxes (of heat, stress, oil, pollutant) entering or exiting a given domain.Furthermore, for a vector field ξ , R ∂ Ω ξ · nds = R Ω div( ξ pot ) dx , which means thatthe flux entering or exiting is determined by the potential part of the vector field. • Classical homogenization is associated with two types of convergence: convergenceof energies (Γ-convergence [33, 15]) and convergence of fluxes ( G or H -convergence[47, 32, 58, 57, 46]). Similarly, one can define an energy norm and a flux-norm. The transfer property.
For V , a finite dimensional linear subspace of H (Ω), wedefine (div a ∇ V ) := span { div( a ∇ v ) : v ∈ V } . (2.4)Note that (div a ∇ V ) is a finite dimensional subspace of H − (Ω). Theorem 2.1. (Transfer property of the flux norm) [Theorem 2.1 of [10]]
Let V ′ and V be finite-dimensional subspaces of H (Ω) . For f ∈ L (Ω) , let u be the solutionof (1.1) with conductivity a and u ′ be the solution of (1.1) with conductivity a ′ . If (div a ∇ V ) = (div a ′ ∇ V ′ ) , then inf v ∈ V k u − v k a -flux k g k L (Ω) = inf v ∈ V ′ k u ′ − v k a ′ -flux k g k L (Ω) . (2.5)3he usefulness of (2.5) can be illustrated by considering a ′ = I so that div a ′ ∇ = ∆.Then, u ′ ∈ H and therefore V ′ can be chosen as, e.g., the standard piecewise linearFEM space, on a regular triangulation of Ω of resolution h , with nodal basis { φ i } . Thespace V is then defined by its basis { θ i } determined by ( div( a ∇ θ i ) = ∆ φ i in Ω θ i = 0 on ∂ Ω . (2.6)Equation (2.5) shows that the approximation error estimate associated with the space V and the problem with arbitrarily rough coefficients is (in a -flux norm) equal to theapproximation error estimate associated with piecewise linear elements and the space H (Ω). More precisely, sup g ∈ L (Ω) inf v ∈ V k u − v k a -flux k g k L (Ω) ≤ Ch, (2.7)where C does not depend on a .We refer to [22], [25] and [11] for recent results on finite element methods for highcontrast ( λ max ( a ) /λ min ( a ) >>
1) but non-degenerate ( λ min ( a ) = O (1)) media underspecific assumptions on the morphology of the (high-contrast) inclusions (in [22], themesh has to be adapted to the morphology of the inclusions). Observe that the pro-posed method remains accurate if the medium is both of high contrast and degenerate( λ min ( a ) << a , at the cost of solving PDEs (2.6)over the whole domain Ω. Remark . We refer to [10] for the optimal constant C in (2.7). This question ofoptimal approximation with respect to a linear finite dimensional space is related to theKolmogorov n-width [54, 44], which measures how accurately a given set of functionscan be approximated by linear spaces of dimension n in a given norm. A surprisingresult of the theory of n-widths is the non-uniqueness of the space realizing the optimalapproximation [54]. Observe also that, as another consequence of the transfer property(2.5), a h k +1 rate of convergence can be achieved in (2.7) by replacing φ i with higher-order basis functions in (2.6), and k g k L with k g k H k in (2.7). Similarly an exponentialrate of convergence can be achieved if the source terms g are analytic. This is thereason behind the near exponential rate of convergence observed in [6] for harmonicfunctions (i.e., with zero source terms, and particular “buffer” solutions computed nearthe boundary) and bounded (non high) contrast media. The elliptic PDEs (2.6) have to be solved on the whole domain Ω. Is it possible tolocalize the computation of the basis elements θ i to a neighborhood of the support ofthe elements φ i ? Observe that the support of each φ i is contained in a ball B ( x i , C h ) ofcenter x i (the node of the coarse mesh associated with x i ) of radius C h . Let 0 < α ≤ φ i ) may, a priori,4ncrease the error estimate in the right hand side of (2.5). This increase can, in fact, belinked to the decay of the Green’s function of the operator − div( a ∇ ). The slower thedecay, the larger the degradation of those approximation error estimates. Inspired bythe strategy used in [35] for controlling cell resonance errors in the computation of theeffective conductivity of periodic or stochastic homogenization (see also [36, 53, 63]), wewill replace the operator − div( a ∇ ) by the operator T − div( a ∇ ) in the left hand sideof (2.6) in order to artificially introduce an exponential decay in the Green’s function.A fine tuning of T is required because although a decrease in T improves the decay ofthe Green function, it also deteriorates the accuracy of the transfer property. In orderto limit this deterioration, we will transfer a vector space with a higher approximationorder than the one associated with piecewise linear elements. Let us now give the mainresult. Let h ∈ (0 , X h be an approximation sub-vector space of H (Ω) such that • X h is spanned by basis functions ( ϕ i ) ≤ i ≤ N (with N = O ( | Ω | /h d )) with supports in B ( x i , C h ) where, the x i are the nodes of a regular triangulation of Ω of resolution h . • X h satisfies the following approximation properties: For all f ∈ H (Ω) ∩ H (Ω)inf v ∈ X h k f − v k H (Ω) ≤ C h k f k H (Ω) , (3.1)and for all f ∈ H (Ω) ∩ H (Ω)inf v ∈ X h k f − v k H (Ω) ≤ C h k f k H (Ω) . (3.2) • For all i , Z Ω |∇ ϕ i | ≤ Ch d − . (3.3) • For all coefficients c i , h d X i c i ≤ C k X i c i ∇ ϕ i k L (Ω) . (3.4) Remark . Examples of such spaces can be found in [17] and constructed using piece-wise quadratic polynomials. From the first bullet point it follows that h can be thoughof as the diameter of the support of the elements ϕ i . The largest parameter h d /C satis-fying (3.4) is the minimal eigenvalue of the stiffness matrix ( R Ω ( ∇ ϕ i ) T ∇ ϕ j ) ≤ i,j ≤ N andCondition (3.4) is obtained from the regularity of the tessellation of Ω. In fact, theproof of Proposition 3.2 shows that Condition (3.4) can be relaxed to the assumption ofexistence of a constant d ϕ > h such that for all coefficients c i h d ϕ X i c i ≤ C k X i c i ∇ ϕ i k L (Ω) . (3.5)5hrough this paper, we will write C any constant that does not depend on h (butmay depend on d , Ω, and the essential supremum and infimum of the maximum andminimum eigenvalues of a over Ω). Let α ∈ (0 ,
1) and C >
0. For each basis element ϕ i of X h let ψ i be the solution of ( h − α ψ i − div( a ∇ ψ i ) = ∆ ϕ i in B ( x i , C h α ln h ) ∩ Ω ψ i = 0 on ∂ (cid:0) B ( x i , C h α ln h ) ∩ Ω (cid:1) . (3.6)Let V h := span( ψ i ) (3.7)be the linear space spanned by the elements ψ i . Theorem 3.1.
For g ∈ L (Ω) , let u be the solution of (1.1) in H (Ω) and u h thesolution of (1.1) in V h . There exists C > such that for C ≥ C , we have k u − u h k H (Ω) k g k L (Ω) ≤ ( Ch if α ∈ (0 , ] Ch − α if α ∈ [ , , (3.8) where the constants C and C depend on a , d , Ω but not on h .Remark . Theorem 3.1 shows the convergence rate in approximation error remainsoptimal (i.e., proportional to h ) after localization if 0 < α ≤ / h − α for ≤ α <
1. In particular, choosing localized domains with radii O ( √ h ln h ) issufficient to obtain the optimal convergence rate O ( h ). Observe that the choice of theconstant α in equation (3.6) is arbitrary. Remark . According to Theorem 3.1, the constant C in (3.6) needs to be chosenlarger than C to achieve the convergence rate h + h − α . The constant C depends on α , d , λ min ( a ) and λ max ( a ). The constant C in the right hand side of (3.8) also dependson α , d , λ min ( a ) and λ max ( a ). It is possible to give an explicit value for C and C bytracking constants in the proof (in particular, as stated in Subsection 3.4, the dependenceon λ max ( a ) can be removed if the elements Ψ i are computed on sub-domains with addedbuffer zones around high-conductivity inclusions). Remark . If one uses piecewise linear basis elements instead of the elements ϕ i (i.e.,in the absence of property (3.2)), then the estimate in the right hand side of (3.8)deteriorates to h − α . The proof of this remark is similar to that of Theorem 3.1. Themain modification lies in replacing h /T by h/T in equations (3.10) and (3.16). Remark . One could use piecewise linear basis elements instead of the elements ϕ i ,and also remove the term h − α ψ i from the transfer property (3.6). In this situation,we numerically observe a rate of convergence of h for periodic, stochastic and low-contrast media after localization of (3.6) to balls of radii O ( h ). In these particular situations (characterized by short range correlations in a ), the term h − α ψ i should beavoided to obtain the optimal convergence rate h after localization to sub-domains of size O ( h ). In that sense, the estimate in the right hand side of (3.8) corresponds to a worst ase scenario with respect to the medium a (characterized by long range correlations),requiring the introduction of the term h − ψ i and a localization to sub-domains of size O ( √ h ln h ) for the optimal convergence rate h . Remark . For the elliptic problem, computational gains result from localization (theelements ψ i are computed on sub-domains Ω i of Ω), parallelization (the elements ψ i canbe computed independently from each other), and the fact that the same basis can beused for different right hand sides g in (1.1). Computational gains are even more signif-icant for time-dependent problems because, once an accurate basis has been determinedfor the elliptic problem, the same basis can be used for the associated (parabolic andhyperbolic) time-dependent problems with the same accuracy (we refer to Sections 4and 5). For the wave equation with rough bulk modulus and density coefficients, theproposed method (based on pre-computing basis elements as solutions of localized el-liptic PDEs) remains accurate, provided that high frequencies are not strongly excited( ∂ t g ∈ L ). On Localization.
We refer to [22], [25] and [6] for recent localization results fordivergence-form elliptic PDEs. The strategy of [22] is to construct triangulations andfinite element bases that are adapted to the shape of high conductivity inclusions viacoefficient dependent boundary conditions for the subgrid problems (assuming a to bepiecewise constant and the number of inclusions bounded). The strategy of [25] is tosolve local eigenvalue problems, observing that only a few eigenvectors are sufficient toobtain a good pre-conditioner. Both [22] and [25] require specific assumptions on themorphology and number of inclusions. The idea of the strategy is to observe that if a is piecewise constant and the number of inclusions bounded, then u is locally H awayfrom the interfaces of the inclusions. The inclusions can then be taken care of by adapt-ing the mesh and the boundary values of localized problems or by observing that thoseinclusions will affect only a finite number of eigenvectors.The strategy of [6] is to construct Generalized Finite Elements by partitioning thecomputational domain into to a collection of preselected subsets and compute optimallocal bases (using the concept of n -widths [55]) for the approximation of harmonic func-tions. Local bases are constructed by solving local eigenvalue problems (correspondingto computing eigenvectors of P ∗ P where P is the restriction of a -harmonic functionsfrom ω ∗ onto ω ⊂ ω ∗ , P ∗ is the adjoint of P , and ω is a sub-domain of Ω surrounded bya larger sub-domain ω ∗ ). The method proposed in [6] achieves a near exponential con-vergence rate (in the number of pre-computed bases functions) for harmonic functions.Non-zero right hand sides ( g ) are then taken care of by solving (for each different g ) par-ticular solutions on preselected subsets with a constant Neumann boundary condition(determined according to the consistency condition).As explained in Remark 2.1, the near exponential rate of convergence observed in [6]is explained by the fact that the source space considered in [6] is more regular than L (since [6] requires the computation particular (local) solutions for each right hand sides g and each non-zero boundary conditions, the basis obtained in [6] is in fact adapted to a -harmonic functions away from the boundary). The strategy proposed here can also7e used to achieve exponential convergence for analytic source terms g by employinghigher-order basis functions ϕ i in (3.6). Furthermore, as shown in sections 4, 5 and 3.4the method proposed here allows for the numerical homogenization of time-dependentproblems (because it does not require the computation of particular solutions for differentsource or boundary terms) and can be extended to high-contrast media. We also notethat the basis functions ψ i are simpler and cheaper to compute (equation (3.6)) thanthe eigenvectors of P ∗ P required by [6]. We refer to page 16 of [6] for a discussion onthe cost of this added complexity. By now, the field of numerical homogenization has become large enough that it is notpossible to give an exhaustive review in this short paper. Therefore, we will restrict ourattention to works directly related to our work.- The multi-scale finite element method [40, 62, 41] can be seen as a numerical gener-alization of this idea of oscillating test functions found in H -convergence. A convergenceanalysis for periodic media revealed a resonance error introduced by the microscopicboundary condition [40, 41]. An over-sampling technique was proposed to reduce theresonance error [40].- Harmonic coordinates play an important role in various homogenization approaches,both theoretical and numerical. These coordinates were introduced in [42] in the con-text of random homogenization. Next, harmonic coordinates have been used in one-dimensional and quasi-one-dimensional divergence form elliptic problems [7, 5], allowingfor efficient finite dimensional approximations. The connection of these coordinates withclassical homogenization is made explicit in [2] in the context of multi-scale finite ele-ment methods. The idea of using particular solutions in numerical homogenization toapproximate the solution space of (1.1) appears to have been first proposed in reservoirmodeling in the 1980s [16], [61] (in which a global scale-up method was introduced basedon generic flow solutions i.e., flows calculated from generic boundary conditions). Itsrigorous mathematical analysis was done only recently [49] and is based on the fact thatsolutions are in fact H -regular with respect to harmonic coordinates (recall that theyare H -regular with respect to Euclidean coordinates). The main message here is that ifthe right hand side of (1.1) is in L , then solutions can be approximated at small scales(in H -norm) by linear combinations of d (linearly independent) particular solutions ( d being the dimension of the space). In that sense, harmonic coordinates are only goodcandidates for being d linearly independent particular solutions.The idea of a global change of coordinates analogous to harmonic coordinates hasbeen implemented numerically in order to up-scale porous media flows [27, 26, 16]. Werefer, in particular, to a recent review article [16] for an overview of some main challengesin reservoir modeling and a description of global scale-up strategies based on genericflows.- In [24, 29], the structure of the medium is numerically decomposed into a micro-scale and a macro-scale (meso-scale) and solutions of cell problems are computed on themicro-scale, providing local homogenized matrices that are transferred (up-scaled) to the8acro-scale grid. This procedure allows one to obtain rigorous homogenization resultswith controlled error estimates for non-periodic media of the form a ( x, xǫ ) (where a ( x, y )is assumed to be smooth in x and periodic or ergodic with specific mixing propertiesin y ). Moreover, it is shown that the numerical algorithms associated with HMM andMsFEM can be implemented for a class of coefficients that is much broader than a ( x, xǫ ).We refer to [34] for convergence results on the Heterogeneous Multiscale Method in theframework of G and Γ-convergence.- More recent work includes an adaptive projection based method [48], which isconsistent with homogenization when there is scale separation, leading to adaptive al-gorithms for solving problems with no clear scale separation; fast and sparse chaos ap-proximations of elliptic problems with stochastic coefficients [60, 37, 23]; finite differenceapproximations of fully nonlinear, uniformly elliptic PDEs with Lipschitz continuousviscosity solutions [19] and operator splitting methods [4, 3].- We refer to [13, 12] (and references therein) for most recent results on homoge-nization of scalar divergence-form elliptic operators with stochastic coefficients. Here,the stochastic coefficients a ( x/ε, ω ) are obtained from stochastic deformations (usingrandom diffeomorphisms) of the periodic and stationary ergodic setting. For each basis element ϕ i of X h , let ψ i,T be the solution of ( T ψ i,T − div( a ∇ ψ i,T ) = ∆ ϕ i in Ω ψ i,T = 0 on ∂ Ω . (3.9)The following Proposition will allow us to control the impact of the introduction of theterm T in the transfer property. Observe that the domain of PDE (3.9) is still Ω (ournext step will be to localize it to Ω i ⊂ Ω).
Proposition 3.1.
For g ∈ L (Ω) let u be the solution of (1.1) in H (Ω) . Then, thereexists v ∈ span( ψ i,T ) such that k u − v k H (Ω) k g k L (Ω) ≤ C (cid:0) h + h T (cid:1) . (3.10) Furthermore, writing v := P i c i ψ i,T we have X i c i ≤ Ch − d (1 + T − ) k g k L (Ω) (3.11) Proof.
Let v = P i c i ψ i,T . We have u − vT − div (cid:0) a ∇ ( u − v ) (cid:1) = g + uT − X i c i ∆ ϕ i . (3.12)9efine a [ v ] to be the energy norm a [ v ] := R Ω ( ∇ v ) T a ∇ v . Multiplying (3.12) by u − v andintegrating by parts, we obtain that k u − v k L (Ω) T + a [ u − v ] = Z Ω ( u − v )( g + uT − X i c i ∆ ϕ i ) . (3.13)Write c i = c i, + c i, and let w and w be the solutions of ∆ w = g − P i c i, ∆ ϕ i and∆ w = uT − P i c i, ∆ ϕ i with Dirichlet boundary conditions on ∂ Ω. Then, we obtain byintegration by parts and the Cauchy-Schwartz inequality that k u − v k L (Ω) T + a [ u − v ] ≤ (cid:13)(cid:13) ∇ ( u − v ) (cid:13)(cid:13) ( L (Ω)) d (cid:0) k∇ w k ( L (Ω)) d + k∇ w k ( L (Ω)) d (cid:1) . (3.14)Using (3.1), we can choose ( c i, ) so that k∇ w k ( L (Ω)) d ≤ Ch k g k L (Ω) . (3.15)Using (3.2), we can choose ( c i, ) so that k∇ w k ( L (Ω)) d ≤ C h T k u k H (Ω) , (3.16)we conclude the proof of the approximation (3.10) by observing that k u k H (Ω) ≤ C k g k L (Ω) .Let us now prove Equation (3.11). First, observe that Equation (3.4) and the triangularinequality imply that (cid:0) X i ( c i ) (cid:1) ≤ Ch − d (cid:16) k X i c i, ∇ ϕ i k L (Ω) + k X i c i, ∇ ϕ i k L (Ω) (cid:17) . (3.17)Next, we obtain from (3.15) and Poincar´e inequality and k X i c i, ∇ ϕ i k L (Ω) ≤ C k g k L (Ω) (3.18)and k X i c i, ∇ ϕ i k L (Ω) ≤ C T k g k L (Ω) (3.19)We conclude by combining equations (3.18) and (3.19) with (3.17).We will now control the error induced by the localization of the elliptic problem (3.9).To this end, for each each basis element ϕ i of X h write S i the intersection of the supportof ϕ i with Ω and let Ω i be a subset of Ω containing S i such that dist( S i , Ω / Ω i ) >
0. Letalso ψ i,T, Ω i be the solution of ( T ψ i,T, Ω i − div( a ∇ ψ i,T, Ω i ) = ∆ ϕ i in Ω i ψ i,T, Ω i = 0 on ∂ Ω i . (3.20)For A, B ⊂ Ω, write d ( A, B ) the Euclidean distance between the sets A and B .10 roposition 3.2. Extending ψ i,T, Ω i by on Ω / Ω i we have (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ Ch d − ( T − + 1) (cid:0) dist( S i , Ω / Ω i ) (cid:1) d +1 exp (cid:16) − dist( S i , Ω / Ω i ) C √ T (cid:17) . (3.21)We refer to Section A of the Appendix for the proof of Proposition 3.2.Taking Ω i := B ( x i , C h α ln h ) ∩ Ω (we use the particular notation C because ourproof of accuracy requires that specific constant to be large enough, i.e., larger thana constant depending on the parameter C appearing in the right hand side of (3.21)and the parameter C describing the balls B ( x i , C h ) containing the support of the basisfunctions ( ϕ i ) ≤ i ≤ N introduced in Subsection 3.1) and T = h α in equation (3.21) ofProposition 3.2, we obtain for C large enough (but independent from h ) that (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ Ch d +1+2 α . (3.22)Let u be the solution of (1.1) in H (Ω). Using Proposition 3.1, we obtain that thereexist coefficients c i such that (cid:13)(cid:13) u − X i c i ψ i,T (cid:13)(cid:13) H (Ω) ≤ C (cid:0) h + h − α (cid:1) k g k L (Ω) . (3.23)and X i c i ≤ Ch − d − α k g k L (Ω) (3.24)Using the triangle inequality, it follows that (cid:13)(cid:13) u − X i c i ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ C (cid:0) h + h − α (cid:1) k g k L (Ω) + X i | c i | (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω) , (3.25)whence, from Cauchy-Schwartz inequality, (cid:13)(cid:13) u − X i c i ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ C (cid:0) h + h − α (cid:1) k g k L (Ω) + (cid:0) X i | c i | (cid:1) (cid:16) X i (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω) (cid:17) . (3.26)Combining (3.26) with (3.24), we obtain that (cid:13)(cid:13) u − X i c i ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ C (cid:0) h + h − α (cid:1) k g k L (Ω) + Ch − d − α k g k L (Ω) (cid:16) X i (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω) (cid:17) . (3.27)Using (3.22) in (3.27), we obtain that (cid:13)(cid:13) u − X i c i ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ C (cid:0) h + h − α (cid:1) k g k L (Ω) . (3.28)Observe that it is the exponential decay in (3.21) that allows us to compensate for thelarge term on the right hand side of (3.27) via (3.22). This concludes the proof ofTheorem 3.1. 11igure 1: Illustrations of the buffer distance. The constant C in the approximation error estimate (3.8) depends, a priori, on thecontrast of a . Is it possible to localize the computation of bases for V h when the contrastof a is high? The purpose of this subsection is to show that the answer is yes providedthat there is a buffer zone between the boundaries of localization sub-domains and thesupports of the elements ϕ i where the contrast of a remains bounded. More precisely,assume that Ω is the disjoint union of Ω bounded and Ω high . Assume that (1.2) holds onlyon Ω bounded , and that on Ω high we have λ min ( a ) | ξ | ≤ ξ T a ( x ) ξ ≤ γ | ξ | . (3.29)where γ can be arbitrarily large. Practical examples include media characterized by abounded contrast background with high conductivity inclusions or channels. Let ψ highi be the solution of ( h − α ψ highi − div( a ∇ ψ highi ) = ∆ ϕ i in Ω i ψ i = 0 on ∂ Ω i . (3.30)Let V highh := span( ψ highi ) (3.31)be the linear space spanned by the elements ψ highi . For each i , define b i to be the largestnumber r such that there exists a subset Ω ′ i such that: the closure of Ω ′ i contains thesupport of ϕ i , (Ω ′ i ) r is a subset of Ω i (where A r are the set of points of Ω that are atdistance at most r for A ), and (Ω ′ i ) r / Ω ′ i is a subset of Ω bounded . If no such subset existswe set b i := 0. b i can be interpreted as the non-high-contrast buffer distance betweenthe support of ϕ i and the boundary of Ω i . We refer to Figure 1 for illustrations of thebuffer distance. 12 heorem 3.2. For g ∈ L (Ω) , let u be the solution of (1.1) in H (Ω) and u h thesolution of (1.1) in V highh . There exists C > such that if for all i , b i ≥ C h α ln h then k u − u h k H (Ω) k g k L (Ω) ≤ ( Ch if α ∈ (0 , ] Ch − α if α ∈ [ , , (3.32) where the constants C and C depend on λ min ( a ) , λ max ( a ) (the bounds on a in Ω bounded ), d , Ω but not on h and γ (The upper bound on a on Ω high ).Remark . Recall that the global basis computed in (2.6) remains accurate if themedium is both of high contrast ( λ max ( a ) >>
1) and degenerate ( λ min ( a ) << C in (3.32) depends on λ min ( a ). Remark . Observe that local solves have to resolve the connected components of highcontrast structures. This is the price to pay for localization with high contrast in themost general case. Recall that in classical homogenization with high contrast the limitof the homogenized operator may be a non-local operator (we refer for instance to [21]).A similar phenomenon is observed here (distant points connected by high conductivitychannels are associated with a low resistance metric and a large coupling coefficient inthe numerically homogenized stiffness matrix).The proof of Theorem 3.2 is similar to that of Theorem 3.1, but it requires a precisetracking of the constants involved. Because of the close similarity we will not includethe proof in this paper but only give its main lines. First, the proof of Proposition3.1 remains unchanged as the constants C in (3.10) and (3.11) do not depend on themaximum eigenvalue of the conductivity a . Only the proof of Proposition 3.2 has tobe adapted and the part of the proof below Proposition 3.2 remains unchanged. Thisrequires an application of the elements of lemmas A.2, A.3, A.4 and A.5 to buffer sub-domains (Ω ′ i ) r / Ω ′ i . The main point is to observe that the decay of the Green’s functionin (Ω ′ i ) r / Ω ′ i can be bounded independently from γ (due to the maximum principle).Observe that the choice of the sub-domain Ω i in (3.30) can be chosen to be the sameas in (3.20) if its intersection with high contrast inclusions is void (i.e., if the maximumeigenvalue of a over Ω i remains bounded independently from γ ); otherwise the choiceof Ω i in (3.30) has to be enlarged (when compared to that associated with (3.20)) tocontain the high-contrast inclusion (plus its buffer). The computational gain of the method proposed in this paper is particularly significantfor time-dependent problems. One such problem is the parabolic equation associatedwith the operator − div( a ∇ ). More precisely, consider the time-dependent partial dif-ferential equation 13 ∂ t u ( x, t ) − div (cid:16) a ( x ) ∇ u ( x, t ) (cid:17) = g ( x, t ) ( x, t ) ∈ Ω T ; g ∈ L (Ω T ) ,u = 0 on ∂ Ω T , (4.1)where a and Ω satisfy the same assumptions as those associated with PDE (1.1), Ω T :=Ω × [0 , T ] for some T > ∂ Ω T := ( ∂ Ω × [0 , T ]) ∪ (Ω × { t = 0 } ).Let V h be the finite-dimensional approximation space defined in (3.7). Let u h be thefinite element solution of (4.1), i.e., u h can be decomposed as u h ( x, t ) = X i c i ( t ) ψ i ( x ) , (4.2)and solves for all j ( ψ j , ∂ t u h ) L (Ω) = − a [ ψ j , u h ] + ( ψ j , g ) L (Ω) , (4.3)with a [ v, w ] := R Ω ( ∇ v ) T a ∇ w . Write k v k L (0 ,T,H (Ω)) := Z T Z Ω |∇ v | ( x, t ) dx dt. (4.4) Theorem 4.1.
We have (cid:13)(cid:13) ( u − u h )( ., T ) (cid:13)(cid:13) L (Ω) + k u − u h k L (0 ,T,H (Ω)) ≤ C k g k L (Ω T ) ( h + h − α ) . (4.5) Proof.
The proof is a generalization of the proof found in [50] (in which approximationspaces are constructed via harmonic coordinates). Let A T be the bilinear form on L (0 , T, H (Ω)) defined by A T [ w , w ] := Z T a [ w , w ] dt. (4.6)Observe that for all v ∈ L (0 , T, V h ), (cid:0) v, ∂ t ( u − u h ) (cid:1) L (Ω T ) + A T [ v, u − u h ] = 0 . (4.7)Writing A T [ v ] := A T [ v, v ], we deduce that for v ∈ L (0 , T, V h ),12 (cid:13)(cid:13) ( u − u h )( ., T ) (cid:13)(cid:13) L (Ω) + A T [ u − u h ] = (cid:0) u − v, ∂ t ( u − u h ) (cid:1) L (Ω T ) + A T [ u − v, u − u h ] . (4.8)Using ∂ t u h in (4.3) and integrating, we obtain that k ∂ t u h k L (Ω T ) + 12 a (cid:2) u h ( ., T ) , u h ( ., T ) (cid:3) = (cid:0) ∂ t u h , g (cid:1) L (Ω T ) . (4.9)14sing Minkowski’s inequality, we deduce that k ∂ t u h k L (Ω T ) + a (cid:2) u h ( ., T ) , u h ( ., T ) (cid:3) ≤ C k g k L (Ω T ) . (4.10)Similarly, k ∂ t u k L (Ω T ) + a (cid:2) u ( ., T ) , u ( ., T ) (cid:3) ≤ C k g k L (Ω T ) . (4.11)Using Cauchy-Schwartz and Minkowski inequalities in (4.8), we obtain that (cid:13)(cid:13) ( u − u h )( ., T ) (cid:13)(cid:13) L (Ω) + A T [ u − u h ] ≤ C k u − v k L (Ω T ) k g k L (Ω T ) + C A T [ u − v ] . (4.12)Take v = R h u to be the projection of u onto L (0 , T, V h ) with respect to the bilinearform A T . Observing that − div( a ∇ u ) = g − ∂ t u with ( g − ∂ t u ) ∈ L (Ω T ), we obtainfrom Theorem 3.1 that (cid:0) A T [ u − R h u ] (cid:1) ≤ C k g k L (Ω T ) ( h + h − α ) . (4.13)Let us now show (using a standard duality argument) that k u − R h u k L (Ω T ) ≤ C ( h + h − α ) k g k L (Ω T ) . (4.14)Choose v ∗ to be the solution of the following linear problem: For all w ∈ L (0 , T, H (Ω)) A T [ w, v ∗ ] = ( w, u − R h u ) L (Ω T ) . (4.15)Taking w = u − R h u in (4.15), we obtain that k u − R h u k L (Ω T ) = A T [ u − R h u, v ∗ − R h v ∗ ] . (4.16)Hence by Cauchy Schwartz inequality and (4.13), k u − R h u k L (Ω T ) ≤ C ( h + h − α ) k g k L (Ω T ) (cid:0) A T [ v ∗ − R h v ∗ ] (cid:1) . (4.17)Using Theorem 3.1 again, we obtain that (cid:0) A T [ v ∗ − R h v ∗ ] (cid:1) ≤ C k u − R h u k L (Ω T ) ( h + h − α ) . (4.18)Combining (4.18) with (4.17) leads to (4.14). Combining (4.12) with v = R h u , (4.14)and (4.13) leads to (cid:13)(cid:13) ( u − u h )( ., T ) (cid:13)(cid:13) L (Ω) + A T [ u − u h ] ≤ C ( h + h − α ) k g k L (Ω T ) , (4.19)which concludes the proof of Theorem 4.1.15 iscretization in time. Let ( t n ) be a discretization of [0 , T ] with time-steps | t n +1 − t n | = ∆ t . Write Z hT , the subspace of L (0 , T, V h ), such that Z hT = ( v ∈ L (0 , T, V h ) : v ( x, t ) = X i c i ( t ) ψ i ( x ) , c i ( t ) are constants on ( t n , t n +1 ] ) . (4.20)Write u h, ∆ t , the solution in Z hT of the following system of implicit weak formulation(such that u h, ∆ t ( x, ≡ n and ψ ∈ V h , (cid:0) ψ, u h, ∆ t ( t n +1 ) (cid:1) L (Ω) = (cid:0) ψ, u h, ∆ t ( t n ) (cid:1) L (Ω) − | ∆ t | a (cid:2) ψ, u h, ∆ t ( t n +1 )] + (cid:0) ψ, Z t n +1 t n g ( t ) dt (cid:1) L (Ω) . (4.21)Then, we have the following theorem Theorem 4.2.
We have (cid:13)(cid:13) ( u − u h, ∆ t )( T ) (cid:13)(cid:13) L (Ω) + k u − u h, ∆ t k L (0 ,T,H (Ω)) ≤ C (cid:0) | ∆ t | + h + h − α (cid:1)(cid:16) k ∂ t g k L (0 ,T,H − (Ω)) + (cid:13)(cid:13) g ( ., (cid:13)(cid:13) L (Ω) + k g k L (Ω T ) (cid:17) . (4.22)The proof of Theorem 4.2 is similar to that of Theorem 1.6 of [50] and will not begiven here. Observe that homogenization in space allows for a discretization in timewith time steps O ( h + h − α ) without compromising the accuracy of the method. Consider the hyperbolic partial differential equation ρ ( x ) ∂ t u ( x, t ) − div (cid:16) a ( x ) ∇ u ( x, t ) (cid:17) = g ( x, t ) ( x, t ) ∈ Ω T ; g ∈ L (Ω T ) ,u = 0 on ∂ Ω T ,∂ t u = 0 on Ω × { t = 0 } , (5.1)where a , Ω, Ω T and ∂ Ω T are defined as in Section 4. In particular, a is assumed to beonly uniformly elliptic and bounded ( a i,j ∈ L ∞ (Ω)). We will further assume that ρ isuniformly bounded from below and above ( ρ ∈ L ∞ (Ω) and essinf ρ ( x ) ≥ ρ min > /h remain weakly excited, because the wavesequation preserves energy and homogenization schemes can not recover energies putinto high frequencies, see [51]). For the sake of conciseness, we will give those resultswith zero boundary conditions. PDE (5.1) corresponds to acoustic wave equations in amedium with density ρ and bulk modulus a − .16et V h be the finite-dimensional approximation space defined in (3.7). Let u h be thefinite element solution of (5.1), i.e., u h can be decomposed as u h ( x, t ) = X i c i ( t ) ψ i ( x ) , (5.2)and solves for all j ( ψ j , ∂ t u h ) L ( ρ, Ω) = − a [ ψ j , u h ] + ( ψ j , g ) L (Ω) , (5.3)where ( v, w ) L ( ρ, Ω) := Z Ω v w ρ. (5.4) Theorem 5.1. If ∂ t g ∈ L (Ω T ) and g ( x, ∈ L (Ω) , then (cid:13)(cid:13) ∂ t ( u − u h )( ., T ) (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) u − u h (cid:13)(cid:13) L (0 ,T,H (Ω)) ≤ C (cid:0) k ∂ t g k L (Ω T ) + k g ( x, k L (Ω) (cid:1) ( h + h − α ) . (5.5) Remark . We refer to [59] for an analysis of the sub-optimal rate of convergenceassociated with finite-difference simulation of wave propagation in discontinuous media(see also [18, 56]). We refer to [51] for an alternative upscaling strategy based on har-monic coordinates. If the medium is locally ergodic with long range correlations [8] andalso characterized by scale separation then we refer to HMM based methods [28, 1].Homogenization based methods require that frequencies larger than 1 /h remain weaklyexcited. For high frequencies, and smooth media (or away from local resonances, e.g.local, nearly resonant cavities), we refer to the sweeping pre-conditioner method [30, 31]. Proof.
Let A T be the bilinear form on L (0 , T, H (Ω)) defined in (4.6). Observe thatfor all v ∈ L (0 , T, V h ), (cid:0) v, ∂ t ( u − u h ) (cid:1) L ( ρ, Ω T ) + A T [ v, u − u h ] = 0 . (5.6)Taking ∂ t u − ∂ t u h − ( ∂ t u − ∂ t v ) as a test function in (5.6) and integrating in time, wededuce that for ∂ t v ∈ L (0 , T, V h ),12 (cid:13)(cid:13) ∂ t ( u − u h )( ., T ) (cid:13)(cid:13) L ( ρ, Ω) + 12 a (cid:2) ( u − u h )( ., T ) (cid:3) = (cid:0) ∂ t ( u − v ) , ∂ t ( u − u h ) (cid:1) L ( ρ, Ω T ) + A T [ ∂ t ( u − v ) , u − u h ] , (5.7)where ( v, w ) L ( ρ, Ω T ) := R T R Ω v w ρ dx dt . Taking the derivative of the hyperbolic equa-tion for u in time, we obtain that ∂ t u − div( a ∇ ∂ t u ) = ∂ t g. (5.8)Integrating (5.8) against the test function ∂ t u and observing that ∂ t u ( x,
0) = g ( x, (cid:13)(cid:13) ∂ t u ( ., T ) (cid:13)(cid:13) L ( ρ, Ω) + a (cid:2) ∂ t u ( ., T ) (cid:3) ≤ C (cid:0) k ∂ t g k L (Ω T ) + k g ( x, k L (Ω) (cid:1) . (5.9)17 L H L ∞ α = 1 / (cid:13)(cid:13) ∂ t u h ( ., T ) (cid:13)(cid:13) L ( ρ, Ω) + a (cid:2) ∂ t u h ( ., T ) (cid:3) ≤ C (cid:0) k ∂ t g k L (Ω T ) + k g ( x, k L (Ω) (cid:1) . (5.10)Take ∂ t v = R h ∂ t u to be the projection of ∂ t u onto L (0 , T, V h ) with respect to thebilinear form A T . Observing that − div( a ∇ ∂ t u ) = ∂ t g − ∂ t u with ( g − ∂ t u ) ∈ L (Ω T ),we obtain from (5.9) and Theorem 3.1 that (cid:0) A T [ u − R h u ] (cid:1) ≤ C (cid:0) k ∂ t g k L (Ω T ) + k g ( x, k L (Ω) (cid:1) ( h + h − α ) . (5.11)Furthermore, using the same duality argument as in the parabolic case, we obtain that k u − R h u k L ( ρ, Ω T ) ≤ C ( h + h − α ) (cid:0) k ∂ t g k L (Ω T ) + k g ( x, k L (Ω) (cid:1) . (5.12)Using Cauchy-Schwartz and Minkowski inequalities and the above estimates in (5.7), weobtain that (cid:13)(cid:13) ∂ t ( u − u h )( ., T ) (cid:13)(cid:13) L ( ρ, Ω) + a (cid:2) ( u − u h )( ., T ) (cid:3) ≤ C ( h + h − α ) (cid:0) A T [ u − u h ] + k ∂ t g k L (Ω T ) + k g ( x, k L (Ω) (cid:1) . (5.13)We conclude using Grownwall’s lemma. We compute the solutions of (1.1) up to time 1 on the fine mesh and in the finite-dimensional approximation space V h defined in (3.7). The physical domain is the square[ − , . Global equations are solved on a fine triangulation with 66049 nodes and 131072triangles.The elements ( ϕ i ) of Sub-section 3.1 are weighted extended B-splines (WEB) [38, 39](obtained by tensorizing one-dimensional elements and using weight function (1 − x )(1 − y ) to enforce the Dirichlet boundary condition). The order of accuracy is not affectedby the choice of weight function given that the boundary is piecewise smooth. Our18 .5 1 1.5 2 2.5−10−9−8−7−6−5−4−3 T= ∞ T=1T=0.5T=0.25T=0.125T=0.0625 (a) L error T= ∞ T=1T=0.5T=0.25T=0.125T=0.0625 (b) H error Figure 2: Example 5 of Section 3 of [49] (percolation at criticality). Logarithm (in base2) of the error with respect to log ( h /h ) (for h = 0 . T used in(3.6).motivation for using WEB elements lies in the fact that, with those elements, (Dirichlet)boundary conditions become simple to enforce. This being said, any finite elementssatisfying the properties (3.1), (3.2), (3.3) and (3.4) would be adequate [14].We write h the size of the coarse mesh. Elements ψ i are obtained by solving (3.6)on localized sub-domains of size h . Table 1 shows errors with α = 1 / a given by(6.1) (Example 1 of Section 3 of [49], trigonometric multi-scale, see also [45]), i.e., for a ( x ) := 16 ( 1 . πx/ǫ )1 . πy/ǫ ) + 1 . πy/ǫ )1 . πx/ǫ ) + 1 . πx/ǫ )1 . πy/ǫ ) +1 . πy/ǫ )1 . πx/ǫ ) + 1 . πx/ǫ )1 . πy/ǫ ) + sin(4 x y ) + 1) , (6.1)where ǫ = , ǫ = , ǫ = , ǫ = , ǫ = .Figure 2 shows the logarithm (in base 2) of the error with respect to log ( h /h ) (for h = 0 . T used in (3.6) for a given by Example 5 of Section 3 of[49] (percolation at criticality, the conductivity of each site is equal to γ or 1 /γ withprobability 1 / γ = 4).Figure 3 shows the logarithm (in base 2) of the error with respect to log ( h /h ) (for h = 0 . T used in (3.6) for a given by Example 3 of Section 3 of [49],i.e., a ( x ) = e h ( x ) , with h ( x ) = P | k |≤ R ( a k sin(2 πk · x ) + b k cos(2 πk · x )), where a k and b k are independent uniformly distributed random variables on [ − . , .
3] and R = 6. Remark . Two factors contribute to the error plots shown in figures 2 and 3: alocalization error which becomes dominant when h /h is small (i.e. the fact (2.6) is notsolved over the whole domain Ω) and the distortion of the transfer property resultingfrom the 1 /T term in (3.9). As expected both figures show that when h /h is large,the error due to the distortion of the transfer property is dominant and is minimizedby a large T . However, when h /h is small, the localization error is dominant and19 .5 1 1.5 2 2.5−8−7−6−5−4−3−2−1 T= ∞ T=1T=0.5T=0.25T=0.125T=0.0625 (a) L error T= ∞ T=1T=0.5T=0.25T=0.125T=0.0625 (b) H error Figure 3: Example 3 of Section 3 of [49] (exponential of a sum of trigonometric functionswith strongly overlapping frequencies). Logarithm (in base 2) of the error with respectto log ( h /h ) (for h = 0 . T used in (3.6).Figure 4: High conductivity channel.is minimized by a small T . The fact that in Figure 3 this error is minimized by thesecond smallest T instead of the smallest T is explained by the fact that the localizationerror remains bounded when h /h is of the order of one whereas the error due to thedistortion of the transfer property blows up as T ↓
0. The fact that in both figures,curves associated with different T interest each other, is indicative of the fact thatfor intermediate values of h /h , the error can be minimized via a fine-tuning of T asexplained in section 3. The differences in the locations of these intersections can beexplained by a larger localization error associated with the example of Figure 3 (dueto longer correlation ranges). In particular, the comparison between figures 2 and 3indicates larger errors for Figure 3. 20 h =C h ln (1/h)h =3h, without bufferh =3h, with buffer (a) L error −4 −3.5 −3 −2.5 −2 −1.5 −1−5.5−5−4.5−4−3.5−3−2.5−2−1.5 h =C h ln (1/h)h =3h, without bufferh =3h, with buffer (b) H error Figure 5: High conductivity channel (Figure 4). The x -axis shows log ( h ), the y -axisshows the log of the error in L and H -norm. The three cases for the localization are h = O ( √ h ln h ) with a buffer around the high conductivity channel (see Sub-section3.4) of size O ( √ h ln h ), h = 3 h with no buffer around the high conductivity channeland h = 3 h with a buffer around the high conductivity channel of size 3 h . In this example, a is characterized by a fine and long-ranged high conductivity channel(Figure 4). We choose a ( x ) = 100, if x is in the channel, and a ( x ) is the percolationmedium, if x is not in the channel (the conductivity of each site, not in channel, is equalto γ or 1 /γ with probability 1 / γ = 4). Figure 5 shows the log of the numericalerror (in L and H norm) versus log ( h ). The three cases for the localization are h = O ( √ h ln h ) with a buffer b i around the high conductivity channel (see Sub-section3.4) of size O ( √ h ln h ), h = 3 h with no buffer around the high conductivity channel and h = 3 h with a buffer b i around the high conductivity channel of size 3 h . The first caseshows that the method of Sub-section 3.4 is converging as expected. The second caseshows that, as expected, taking α = 1, does not guarantee convergence. The third caseshows that adding a buffer around the high conductivity channel improves numericalerrors but is not sufficient to guarantee convergence (as expected, we also need α < h (i.e., log ( h ) = − We compute the solutions of (5.1) up to time 1 on the fine mesh and in the finite-dimensional approximation space V h defined in (3.7). The initial condition is u ( x,
0) = 0and u t ( x,
0) = 0. The boundary condition is u ( x, t ) = 0, for x ∈ ∂ Ω. The density isuniformly equal to one and we choose g = sin( πx ) sin( πy ). Figure 6 shows the fine meshsolutions u and u h at time one, for a given by the trigonometric example (6.1), with h = 0 . h = 3 h and T = h . Figure 6 shows the fine mesh solutions u and u h at time21 a) u (b) u h Figure 6: Wave equation. Trigonometric case, fine mesh solution, h = 0 . h = 3 h , T = h . The L , H and L ∞ relative numerical errors are 0 . . . a given by the high conductivity channel example (Figure 4), with h = 0 . h = 3 h and T = h .We refer to [52] for a list of movies on the numerical homogenization of the waveequation with and without high contrast and with and without buffers (extended buffersin the high contrast case). A Proof of Proposition 3.2.
The proof of Proposition 3.2 is a generalization of the proof of the control of the resonanceerror in periodic medium given in [35].First we need the following lemma, which is the cornerstone of Cacciopoli’s inequality.
Lemma A.1.
Let D be a sub-domain of Ω with piecewise Lipschitz boundary, and let v solve ( vT − div (cid:16) a ( x ) ∇ v ( x ) (cid:17) = f ( x ) x ∈ D ; f ∈ H − ( D ) ,v = 0 on ∂D, (A.1) Let ζ : D → R + be a function of class C such that ζ is identically null on an openneighborhood of the support of f . Then, Z D (cid:12)(cid:12) ∇ ( ζv ) (cid:12)(cid:12) ≤ C Z D v |∇ ζ | , (A.2) where C only depends on the essential supremum and infimum of the maximum andminimum eigenvalues of a over D .Proof. Multiplying (A.1) by ζ v and integrating by parts, we obtain that Z D ζ v T + Z D ∇ ( ζ v ) a ∇ v = 0 . (A.3)22 (a) u −1 −0.5 0 0.5 1−1−0.500.51−0.03−0.02−0.0100.010.020.030.040.050.06 −0.02−0.0100.010.020.030.040.05 (b) u h Figure 7: Wave equation. Channel case, coarse mesh solution, h = 0 . h = 3 h , T = h . The L , H and L ∞ relative numerical errors are 0 . . . Z D ζ v T + Z D ∇ ( ζv ) a ∇ ( ζv ) = Z D v ∇ ζa ∇ ζ, (A.4)which concludes the proof. Lemma A.2.
Let D be a sub-domain of Ω with piecewise Lipschitz boundary. Write G T,D the Green’s function of the operator T − div( a ∇ ) with Dirichlet boundary conditionon ∂D . Then, G T,D ( x, y ) ≤ C | x − y | d − exp (cid:0) − | x − y | C √ T (cid:1) , (A.5) where C only depends on d and the essential supremum and infimum of the maximumand minimum eigenvalues of a over D .Proof. Extending a to R d and using the maximum principle, we obtain that G T,D ( x, y ) ≤ G T, R d ( x, y ) , (A.6)we conclude by using the exponential decay of the Green’s function in R d (we refer toLemma 2 of [35]). Lemma A.3.
Let ψ i,T be the solution of (3.9) and ψ i,T, Ω i the solution of (3.20) . Let Ω ′ i be a sub-domain of Ω i such that S i ⊂ Ω ′ i and dist( S i , Ω i / Ω ′ i ) > . We have (cid:13)(cid:13) ψ i,T (cid:13)(cid:13) H (Ω / Ω ′ i ) ≤ Ch d − (cid:0) dist( S i , Ω / Ω ′ i ) (cid:1) d exp (cid:16) − dist( S i , Ω / Ω ′ i ) C √ T (cid:17) , (A.7) and (cid:13)(cid:13) ψ i,T, Ω i (cid:13)(cid:13) H (Ω i / Ω ′ i ) ≤ Ch d − (cid:0) dist( S i , Ω / Ω ′ i ) (cid:1) d exp (cid:16) − dist( S i , Ω / Ω ′ i ) C √ T (cid:17) . (A.8)23 roof. For A ⊂ Ω, write A r the set of points of Ω that are at distance at most r from A . Let us now use Cacciopoli’s inequality to bound R Ω / Ω ′ i |∇ ψ i,T | . Using Lemma A.1with ζ identically equal to one on Ω / Ω ′ i , zero on (Ω / Ω ′ i ) r with r := dist( S i , Ω / Ω ′ i ) / |∇ ζ | ≤ C/r , we obtain that Z Ω / Ω ′ i |∇ ψ i,T | ≤ Cr Z (Ω / Ω ′ i ) r ψ i,T . (A.9)Next, observe that for x ∈ (Ω / Ω ′ i ) r , ψ i,T ( x ) = − Z S i ∇ G T, Ω ( x, y ) ∇ ϕ i ( y ) dy. (A.10)Hence, (cid:12)(cid:12) ψ i,T ( x ) (cid:12)(cid:12) ≤ k∇ ϕ i k ( L ( S i )) d (cid:13)(cid:13) ∇ G T, Ω ( x, . ) (cid:13)(cid:13) ( L ( S i )) d . (A.11)Another use of Cacciopoli’s inequality leads to (cid:13)(cid:13) ∇ G T, Ω ( x, . ) (cid:13)(cid:13) ( L ( S i )) d ≤ Cr (cid:13)(cid:13) G T, Ω ( x, . ) (cid:13)(cid:13) L ( S ri ) . (A.12)Combining (A.9) with (A.11) with (A.12), we obtain that Z Ω / Ω ′ i |∇ ψ i,T | ≤ k∇ ϕ i k L ( S i )) d Cr Z (Ω / Ω ′ i ) r (cid:13)(cid:13) G T, Ω ( x, . ) (cid:13)(cid:13) L ( S ri ) . (A.13)We conclude the proof of (A.7) using Lemma A.2 and (3.3). The proof of (A.8) is similarobserving that dist( S i , Ω / Ω ′ i ) ≤ dist( S i , Ω i / Ω ′ i ) Lemma A.4.
Let D be a sub-domain of Ω with piecewise Lipschitz boundary. Let ψ ∈ H (Ω) , and let v solve ( vT − div (cid:16) a ( x ) ∇ v ( x ) (cid:17) = 0 x ∈ D,v = ψ on ∂D, (A.14) Write S the intersection of the support of ψ with D . Let D be a sub-domain of D suchthat dist( D , S ) > , then Z D |∇ v | ≤ C (cid:0) dist( D , S ) (cid:1) d ( T − + 1) k ψ k H (Ω) exp (cid:16) − dist( D , S ) C √ T (cid:17) , (A.15) where C does not depend on D, D , S .Proof. Write w := v − ψ . Then, ( wT − div (cid:16) a ( x ) ∇ w ( x ) (cid:17) = − ψT + div( a ∇ ψ ) x ∈ D,v = 0 on ∂D, (A.16)24hus, w ( x ) = − Z D (cid:0) ψ ( y ) T G
T,D ( x, y ) + ∇ ψ ( y ) a ( y ) ∇ G T,D ( x, y ) (cid:1) dy. (A.17)Using Cauchy-Schwartz inequality, we obtain that | w ( x ) | ≤ C k ψ k H (Ω) (cid:16) T (cid:13)(cid:13) G T,D ( x, . ) (cid:13)(cid:13) L ( S ) + (cid:13)(cid:13) ∇ G T,D ( x, . ) (cid:13)(cid:13) ( L ( S )) d (cid:17) . (A.18)For A ⊂ D , write A r the set of points of D that are at distance at most r from A .Let us now use Cacciopoli’s inequality to bound R D |∇ w | . Using Lemma A.1 with ζ identically equal to one on D , zero on D/D r and such that |∇ ζ | ≤ C/r we obtainthat Z D |∇ w | ≤ Cr Z D r w , (A.19)provided that dist( D r , S ) >
0. Hence, for r := dist( D , S ) /
3, we obtain (A.19). Taking r := dist( D , S ) / (cid:13)(cid:13) ∇ G T,D ( x, . ) (cid:13)(cid:13) ( L ( S )) d ≤ Cr (cid:13)(cid:13) G T,D ( x, . ) (cid:13)(cid:13) L ( S r ) . (A.20)Combining (A.19) with (A.18) and (A.20) and observing that w = v on D r we obtainthat Z D |∇ v | ≤ Cr r k ψ k H (Ω) ( T − + 1) Z D r (cid:13)(cid:13) G T,D ( x, . ) (cid:13)(cid:13) L ( S r ) . (A.21)Using Lemma A.2, we deduce that Z D |∇ w | ≤ C | Ω | (dist( D , S )) d k ψ k H (Ω) ( T − + 1) exp (cid:0) − dist( D , S ) C √ T (cid:1) . (A.22)This concludes the proof of Lemma A.4. Lemma A.5.
Let ψ i,T be the solution of (3.9) and ψ i,T, Ω i the solution of (3.20) . Let Ω ′ i be a sub-domain of Ω i such that dist(Ω / Ω i , Ω ′ i ) > . We have (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω ′ i ) ≤ C ( T − + 1) h d − (cid:0) dist(Ω / Ω i , Ω ′ i ) (cid:1) d +1 exp (cid:16) − dist(Ω / Ω i , Ω ′ i ) C √ T (cid:17) . (A.23) Proof.
Lemma A.5 is a direct consequence of Lemma A.4. To this end, we choose D := Ω i , v =: ψ i,T − ψ i,T, Ω i and D := Ω ′ i . We also choose ψ := ηψ i,T where η : Ω → [0 , C , equal to one on Ω / Ω i and 0 on (Ω / Ω i ) r with r := dist(Ω / Ω i , Ω ′ i ) / A r being theset of points in Ω at distance at most r from A ) and |∇ η | ≤ C/r . We obtain fromLemma A.4 that (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω ′ i ) ≤ C ( T − + 1) (cid:0) dist(Ω / Ω i , Ω ′ i ) (cid:1) d k ψ k H (Ω) exp (cid:16) − dist(Ω / Ω i , Ω ′ i ) C √ T (cid:17) . (A.24)We conclude using (3.3) and k ψ k H (Ω) ≤ C dist(Ω / Ω i , Ω ′ i ) k∇ ϕ i k ( L (Ω)) d .25bserving that (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω) ≤ (cid:13)(cid:13) ψ i,T − ψ i,T, Ω i (cid:13)(cid:13) H (Ω ′ i ) + (cid:13)(cid:13) ψ i,T (cid:13)(cid:13) H (Ω / Ω ′ i ) + (cid:13)(cid:13) ψ i,T, Ω i (cid:13)(cid:13) H (Ω i / Ω ′ i ) , (A.25)we conclude the proof of Proposition 3.2 by using Lemma A.5 and Lemma A.3 with Ω ′ i := S ri where S ri are the points in Ω i at distance at most r from S i with r := dist( S i , Ω / Ω i ) / B On the compactness of the solution space.
Although the foundations of classical homogenization [9] were laid down based on as-sumptions of periodicity (or ergodicity) and scale separation, numerical homogenization,as described here, is independent from these concepts and solely relies on the strong com-pactness of the solution space (and the fact that a compact set can be covered with afinite number of balls of arbitrary sizes). Observe that an analogous notion of com-pactness supports the foundations of G and H -convergence ([47], [57, 32]). The maindifference is that G and H -convergence rely on pre-compactness and weak convergenceof fluxes and here, we rely on compactness in the (strong) H -norm, i.e. the followingtheorem.Let W be the range of g in (1.1). Write V := { u ∈ H (Ω) : u solves (1.1) for some g ∈ W } . (B.1) Theorem B.1.
Let ν < . If W is a closed bounded subset of H − ν (Ω) then W is acompact subset of H (Ω) (in the strong H -norm).Proof. We have ( a ∇ u ) pot = −∇ ∆ − g . So using the same notation as in (2.4) we get( a ∇ V ) pot = −∇ ∆ − W . Let u n be a sequence in V then there exists a sequence in W such that − div( a ∇ u n ) = g n . Using the fact that −∇ ∆ − W is a compact subsetof ( L (Ω)) d (we refer, for instance, to the Kondrachov embedding theorem) we getthat there exists g ∗ ∈ W such that k∇ ∆ − g n − ∇ ∆ − g ∗ k L →
0. Writing u ∗ thesolution of − div( a ∇ u ∗ ) = g ∗ and using ( a ∇ ( u n − u ∗ )) pot = −∇ ∆ − ( g n − g ∗ ) we get that k ( a ∇ ( u n − u ∗ )) pot k L →
0. Using the equivalence between the flux norm and the H norm we deduce that k u n − u ∗ k H →
0. This finishes the proof.This notion of compactness of the solution space constitutes a simple but fundamentallink between classical homogenization, numerical homogenization and reduced ordermodeling (or reduced basis modeling [20, 43]) (we also refer to the discussion in Section6 of [10]). This notion is also what allows for atomistic to continuum up-scaling [64], thebasic idea is that if source (force) terms are integrable enough (for instance in L insteadof H − ) then the solution space is no longer H but a sub-space V that is compactlyembedded into H and, hence, it can be approximated by a finite-dimensional space (in H -norm). In other words if these systems are “excited” by “regular” forces or sourceterms (think compact, low dimensional) then the solution space can be approximated bya low dimensional space (of the whole space) and the name of the game becomes “how to26pproximate” this solution space (and this can be done by using local time-independentsolutions). Acknowledgements.
We thank L. Berlyand for stimulating discussions. We alsothank Ivo Babuˇska, John Osborn, George Papanicolaou and Bj¨orn Engquist for pointingus in the direction of the localization problem. The work of H. Owhadi is partially sup-ported by the National Science Foundation under Award Number CMMI-092600 and theDepartment of Energy National Nuclear Security Administration under Award NumberDE-FC52-08NA28613. The work of L. Zhang is partially supported by the EPSRC Sci-ence and Innovation award to the Oxford Centre for Nonlinear PDE (EP/E035027/1).We thank Sydney Garstang for proofreading the manuscript.
References [1] A. Abdulle and M.J. Grote. Finite element heterogeneous multiscale method forthe wave equation.
Submitted , 2010.[2] G. Allaire and R. Brizzi. A multiscale finite element method for numerical homog-enization.
Multiscale Model. Simul. , 4(3):790–812 (electronic), 2005.[3] T. Arbogast and K. J. Boyd. Subgrid upscaling and mixed multiscale finite elements.
SIAM J. Numer. Anal. , 44(3):1150–1171 (electronic), 2006.[4] T. Arbogast, C.-S. Huang, and S.-M. Yang. Improved accuracy for alternating-direction methods for parabolic equations based on regular and mixed finite ele-ments.
Math. Models Methods Appl. Sci. , 17(8):1279–1305, 2007.[5] I. Babuˇska, G. Caloz, and J. E. Osborn. Special finite element methods for a classof second order elliptic problems with rough coefficients.
SIAM J. Numer. Anal. ,31(4):945–981, 1994.[6] I. Babuˇska and R. Lipton. Optimal local approximation spaces for generalized finiteelement methods with application to multiscale problems.
Multiscale Model. Simul. ,9:373–406, 2011.[7] I. Babuˇska and J. E. Osborn. Generalized finite element methods: their performanceand their relation to mixed methods.
SIAM J. Numer. Anal. , 20(3):510–536, 1983.[8] G. Bal and W. Jing. Corrector theory for msfem and hmm in random media. arXiv:1011.5194 , 2010.[9] A. Bensoussan, J. L. Lions, and G. Papanicolaou.
Asymptotic analysis for periodicstructure . North Holland, Amsterdam, 1978.[10] L. Berlyand and H. Owhadi. Flux norm approach to finite dimensional homoge-nization approximations with non-separated scales and high contrast.
Archives forRational Mechanics and Analysis , 198(2):677–721, 2010.2711] C. Bernardi and R. Verf¨urth. Adaptive finite element methods for elliptic equationswith non-smooth coefficients.
Numer. Math. , 85(4):579–608, 2000.[12] X. Blanc, C. Le Bris, and P.-L. Lions. Une variante de la th´eorie del’homog´en´eisation stochastique des op´erateurs elliptiques.
C. R. Math. Acad. Sci.Paris , 343(11-12):717–724, 2006.[13] X. Blanc, C. Le Bris, and P.-L. Lions. Stochastic homogenization and randomlattices.
J. Math. Pures Appl. (9) , 88(1):34–63, 2007.[14] Dietrich Braess.
Finite Elements: Theory, Fast Solvers, and Applications in SolidMechanics . Cambridge Universty Press, 2007.[15] A. Braides. Γ -convergence for beginners , volume 22 of
Oxford Lecture Series inMathematics and its Applications . Oxford University Press, Oxford, 2002.[16] L. V. Branets, S. S. Ghai, L. L., and X.-H. Wu. Challenges and technologies inreservoir modeling.
Commun. Comput. Phys. , 6(1):1–23, 2009.[17] S. C. Brenner and L. R. Scott.
The mathematical theory of finite elements methods ,volume 15 of
Texts in Applied Mathematics . Springer, 2002. Second edition.[18] D. L. Brown. A note on the numerical solution of the wave equation with piecewisesmooth coefficients.
Mathematics of Computation , 42(166):369–391, 1984.[19] L. A. Caffarelli and P. E. Souganidis. A rate of convergence for monotone finitedifference approximations to fully nonlinear, uniformly elliptic PDEs.
Comm. PureAppl. Math. , 61(1):1–17, 2008.[20] Eric Canc`es, Claude Le Bris, Yvon Maday, Ngoc Cuong Nguyen, Anthony T. Patera,and George Shu Heng Pau. Feasibility and competitiveness of a reduced basisapproach for rapid electronic structure calculations in quantum chemistry. In
High-dimensional partial differential equations in science and engineering , volume 41 of
CRM Proc. Lecture Notes , pages 15–47. Amer. Math. Soc., Providence, RI, 2007.[21] K. D. Cherednichenko, V. P. Smyshlyaev, and V. V. Zhikov. Non-local homogenizedlimits for composite media with highly anisotropic periodic fibres.
Proc. Roy. Soc.Edinburgh Sect. A , 136(1):87–114, 2006.[22] C.-C. Chu, I. G. Graham, and T. Y. Hou. A new multiscale finite element methodfor high-contrast elliptic interface problems.
Math. Comp. , 79:1915–1955, 2010.[23] A. Doostan and H. Owhadi. A non-adapted sparse approximation of pdes withstochastic inputs. arXiv:1006.2151 , 2010.[24] W. E, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden. Heterogeneous multi-scale methods: a review.
Commun. Comput. Phys. , 2(3):367–450, 2007.2825] Y. Efendiev, J. Galvis, and X. Wu. Multiscale finite element and domain decompo-sition methods for high-contrast problems using local spectral basis functions. 2009.Submitted.[26] Y. Efendiev, V. Ginting, T. Hou, and R. Ewing. Accurate multiscale finite elementmethods for two-phase flow simulations.
J. Comput. Phys. , 220(1):155–174, 2006.[27] Y. Efendiev and T. Hou. Multiscale finite element methods for porous media flowsand their applications.
Appl. Numer. Math. , 57(5-7):577–596, 2007.[28] B. Engquist, H. Holst, and O. Runborg. Multi-scale methods for wave propagationin heterogeneous media.
To appear in Commun. Math. Sci. , 2010.[29] B. Engquist and P. E. Souganidis. Asymptotic and numerical homogenization.
ActaNumerica , 17:147–190, 2008.[30] B. Engquist and L. Ying. Sweeping preconditioner for the helmholtz equation:Hierarchical matrix representation.
Submitted , 2010.[31] B. Engquist and L. Ying. Sweeping preconditioner for the helmholtz equation:Moving perfectly matched layers.
Submitted , 2010.[32] E. De Giorgi. Sulla convergenza di alcune successioni di integrali del tipo dell’aera.
Rendi Conti di Mat. , 8:277–294, 1975.[33] E. De Giorgi. New problems in Γ-convergence and G -convergence. In Free boundaryproblems, Vol. II (Pavia, 1979) , pages 183–194. Ist. Naz. Alta Mat. Francesco Severi,Rome, 1980.[34] A. Gloria. Analytical framework for the numerical homogenization of elliptic mono-tone operators and quasiconvex energies.
SIAM MMS , 5(3):996–1043, 2006.[35] A. Gloria. Reduction of the resonance error. part 1: Approximation of homogenizedcoefficients.
To appear in M3AS , 2010. HAL: inria-00457159, version 1.[36] A. Gloria and F. Otto. An optimal error estimate in stochastic homogenization ofdiscrete elliptic equations.
To appear , 2010. HAL: inria-00457020, version 1.[37] H. Harbrecht, R. Schneider, and C. Schwab. Sparse second moment analysis forelliptic problems in stochastic domains.
Numer. Math. , 109(3):385–414, 2008.[38] Klaus H¨ollig.
Finite element methods with B-splines , volume 26 of
Frontiers inApplied Mathematics . Society for Industrial and Applied Mathematics (SIAM),Philadelphia, PA, 2003.[39] Klaus H¨ollig, Ulrich Reif, and Joachim Wipper. Weighted extended B-spline approx-imation of Dirichlet problems.
SIAM J. Numer. Anal. , 39(2):442–462 (electronic),2001. 2940] T. Y. Hou and X.-H. Wu. A multiscale finite element method for elliptic problemsin composite materials and porous media.
J. Comput. Phys. , 134(1):169–189, 1997.[41] T. Y. Hou, X.-H. Wu, and Z. Cai. Convergence of a multiscale finite elementmethod for elliptic problems with rapidly oscillating coefficients.
Math. Comp. ,68(227):913–943, 1999.[42] S. M. Kozlov. The averaging of random operators.
Mat. Sb. (N.S.) , 109(151)(2):188–202, 327, 1979.[43] Luc Machiels, Yvon Maday, Ivan B. Oliveira, Anthony T. Patera, and Dimitrios V.Rovas. Output bounds for reduced-basis approximations of symmetric positivedefinite eigenvalue problems.
C. R. Acad. Sci. Paris S´er. I Math. , 331(2):153–158,2000.[44] J. M. Melenk. On n -widths for elliptic problems. J. Math. Anal. Appl. , 247(1):272–289, 2000.[45] Pingbing Ming and Xingye Yue. Numerical methods for multiscale elliptic problems.
J. Comput. Phys. , 214(1):421–445, 2006.[46] F. Murat. Compacit´e par compensation.
Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) ,5(3):489–507, 1978.[47] F. Murat and L. Tartar. H-convergence.
S´eminaire d’Analyse Fonctionnelle etNum´erique de l’Universit´e d’Alger , 1978.[48] J. Nolen, G. Papanicolaou, and O. Pironneau. A framework for adaptive multiscalemethods for elliptic problems.
Multiscale Model. Simul. , 7(1):171–196, 2008.[49] H. Owhadi and L. Zhang. Metric-based upscaling.
Comm. Pure Appl. Math. ,60(5):675–723, 2007.[50] H. Owhadi and L. Zhang. Homogenization of parabolic equations with a continuumof space and time scales.
SIAM J. Numer. Anal. , 46(1):1–36, 2007/08.[51] H. Owhadi and L. Zhang. Homogenization of the acoustic wave equation with acontinuum of scales.
Computer Methods in Applied Mechanics and Engineering ,198(2-4):97–406, 2008. Arxiv math.NA/0604380.[52] H. Owhadi and L. Zhang. Numerical homogenization with localized bases.
Youtube ,2010. .[53] G. C. Papanicolaou and S. R. S. Varadhan. Boundary value problems with rapidlyoscillating random coefficients. In
Random fields, Vol. I, II (Esztergom, 1979) ,volume 27 of
Colloq. Math. Soc. J´anos Bolyai , pages 835–873. North-Holland, Am-sterdam, 1981.[54] A. Pinkus. n-Width in Approximation Theory . Springer, New York, 1985.3055] Allan Pinkus. n -widths in approximation theory , volume 7 of Ergebnisse der Math-ematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)] .Springer-Verlag, Berlin, 1985.[56] A. Sei and W. W. Symes. Error analysis of numerical schemes for the wave equationin heterogeneous media.
Applied Numerical Mathematics , 15(4):465–480, 1994.[57] S. Spagnolo. Sulla convergenza di soluzioni di equazioni paraboliche ed ellittiche.
Ann. Scuola Norm. Sup. Pisa (3) 22 (1968), 571-597; errata, ibid. (3) , 22:673,1968.[58] S. Spagnolo. Convergence in energy for elliptic operators. In
Numerical solution ofpartial differential equations, III (Proc. Third Sympos. (SYNSPADE), Univ. Mary-land, College Park, Md., 1975) , pages 469–498. Academic Press, New York, 1976.[59] W. W. Symes and T. Vdovina. Interface error analysis for numerical wave propa-gation.
Computational Geosciences , 13(3):363–371, 2009.[60] R. A. Todor and C. Schwab. Convergence rates for sparse chaos approximations ofelliptic problems with stochastic coefficients.
IMA J. Numer. Anal. , 27(2):232–261,2007.[61] C. D. White and R. N. Horne. Computing absolute transmissibility in the presenceof finescale heterogeneity.
SPE Symposium on Reservoir Simulation , page 16011,1987.[62] X. H. Wu, Y. Efendiev, and T. Y. Hou. Analysis of upscaling absolute permeability.
Discrete Contin. Dyn. Syst. Ser. B , 2(2):185–204, 2002.[63] V. V. Yurinski˘ı. Averaging of symmetric diffusion in a random medium.
Sibirsk.Mat. Zh. , 27(4):167–180, 215, 1986.[64] L. Zhang, L. Berlyand, M. Federov, and H. Owhadi. Global energy matching methodfor atomistic to continuum modeling of self-assembling biopolymer aggregates. sub-mittedsub-mitted