Exponential decay of the resonance error in numerical homogenization via parabolic and elliptic cell problems

TL;DR

This paper introduces two innovative methods using parabolic and elliptic cell problems to compute homogenized coefficients in multiscale elliptic PDEs, achieving exponential decay of resonance error and improving accuracy over traditional techniques.

Contribution

The paper proposes novel approaches based on parabolic and elliptic cell problems that significantly reduce resonance error in homogenization computations.

Findings

Resonance error decays exponentially with the new methods.

New approaches outperform standard techniques in accuracy.

Methods are applicable to multiscale elliptic PDEs.

Abstract

This paper presents two new approaches for finding the homogenized coefficients of multiscale elliptic PDEs. Standard approaches for computing the homogenized coefficients suffer from the so-called resonance error, originating from a mismatch between the true and the computational boundary conditions. Our new methods, based on solutions of parabolic and elliptic cell-problems, result in an exponential decay of the resonance error.

Equations

$$\{\begin{array}{rlrl}-\nabla \cdot \left(a^{\epsilon }(x)\nabla u^{\epsilon }\right)& =f& & \text{in }\Omega \subset {\mathbb{R}}^d\\ u^{\epsilon }& =0& & \text{on }\partial \Omega ,\end{array}$$

$$a_{ij}^0=\frac{1}{\left|K\right|}{\int }_K\left({\mathbf{e}}_i+\nabla {\chi }^i(y)\right)\cdot a(y)\left({\mathbf{e}}_j+\nabla {\chi }^j(y)\right)dy,$$

$$-\nabla\cdot\Big{(}a(y)\left(\nabla\chi^{i}+\mathbf{e}_{i}\right)\Big{)}=0\quad\text{in }K.$$

$$a^{0,R}_{ij}=\frac{1}{\left\lvert K_{R}\right\rvert}\int_{K_{R}}\big{(}\mathbf{e}_{i}+\nabla\psi^{i}_{R}(y)\big{)}\cdot a(y)\big{(}\mathbf{e}_{j}+\nabla\psi^{j}_{R}(y)\big{)}\,dy,$$

$$\left\{\begin{aligned} -\nabla\cdot\Big{(}a(y)\left(\nabla\psi^{i}_{R}+\mathbf{e}_{i}\right)\Big{)}&=0&\quad&\text{in }K_{R}\\ \psi^{i}_{R}&=0&\quad&\text{on }\partial K_{R}.\end{aligned}\right.$$

$$\{\begin{array}{rlrl}& \frac{\partial u_R^i}{\partial t}-\nabla \cdot (a(y)\nabla u_R^i)=0& & \text{in }{K}_R\times (0,T]\\ & u_R^i=0& & \text{on }\partial K_R\times (0,T]\\ & u_R^i(y,0)=\nabla \cdot (a(y){\mathbf{e}}_i)& & \text{in }{K}_R.\end{array}$$

$$a_{ij}^{0,R,L,T}:={\int }_{K_L}{\mathbf{e}}_i\cdot a(y){\mathbf{e}}_j{\mu }_L(y)dy-2{\int }_0^T{\int }_{K_L}{u}_R^i(y,t)u_R^j(y,t){\mu }_L(y)dydt.$$

$$\left\{\begin{aligned} -\nabla\cdot\Big{(}a(y)\left(\nabla\chi^{i}_{R,T,N}+\mathbf{e}_{i}\right)\Big{)}+[e^{-A_{N}T}g^{i}](y)&=0&\quad&\text{in }K_{R}\\ \chi^{i}_{R,T,N}&=0&\quad&\text{on }\partial K_{R},\end{aligned}\right.$$

$$[e^{-A_{N}T}{g}^i](y):=\sum _{k=1}^N{e}^{-{\lambda }_{k}T}{g}_k^j{\phi }_k(y),$$

$$b_{ij}^{0,R,L,T,N}={\int }_{K_L}\left(a_{ij}(y)+\sum _{k=1}^d{a}_{ik}(y){\partial }_k{\chi }_{R,T,N}^j(y)\right){\mu }_L(y)dy.$$

$$X_0(R_{+},\Omega ):=\left\{v\in L^2\left(R_{+};H_0^1(\Omega )\right),{\partial }_{t}v\in L^2\left(R_{+};H^{-1}(\Omega )\right)\right\}.$$

$${\psi }_R^k(y)={\int }_0^{+\infty }{u}_R^k(y,t)dt,$$

$$\frac{1}{2}{\int }_{K_R}\nabla {\psi }_R^k(y)\cdot a(y)\nabla {\psi }_R^j(y)dy={\int }_0^{+\infty }{\int }_{K_R}{u}_R^k(y,t)u_R^j(y,t)dydt.$$

$$\{\begin{array}{rl}& \frac{d{u}_R^k}{dt}+A{u}_R^k=0\\ & u_R^k(0)=g^k,g^k(y)=\nabla \cdot \left(a(y){\mathbf{e}}_k\right)\text{ in }{L}^2(K_R).\end{array}$$

$$S_{\omega }=\left\{z\in \mathbb{C}:\left|\text{arg }z\right|<\omega ,0<\omega <\frac{\pi }{2}\right\}.$$

$$e^{-tA}=\frac{1}{2\pi i}{\int }_{\Gamma }{e}^{-tz}(zI-A{)}^{-1}dz$$

$${\int }_0^{+\infty }{u}_R^k(t)dt$$

$$=\frac{1}{2\pi i}{\int }_{\Gamma }{\int }_0^{+\infty }{e}^{-tz}dt{\left(zI-A\right)}^{-1}{g}^{k}dz$$

$$=\frac{1}{2\pi i}{\int }_{\Gamma }\frac{1}{z}{\left(zI-A\right)}^{-1}{g}^{k}dz=A^{-1}{g}^k.$$

$${\int }_{K_R}\nabla {\psi }_R^j\cdot a(y)\nabla {\psi }_R^{k}dy$$

$$={\int }_0^{+\infty }{\left(\nabla \cdot \left(a{\mathbf{e}}_k\right),u_R^j\right)}_{L^2(K_R)}dt.$$

$${\int }_{K_R}\nabla {\psi }_R^k(y)\cdot a(y)\nabla {\psi }_R^j(y)dy={\int }_0^{+\infty }{\left(u_R^k(\cdot ,t/2),u_R^j(\cdot ,t/2)\right)}_{L^2(K_R)}dt,$$

$$a_{ij}^{0,R}=\frac{1}{\left|K_R\right|}{\int }_{K_R}{\mathbf{e}}_i\cdot a(y){\mathbf{e}}_{j}dy-\frac{1}{\left|K_R\right|}{\int }_{K_R}\nabla {\psi }_R^i(y)\cdot a(y)\nabla {\psi }_R^j(y)dy.$$

$$a^{0,R,R,+\infty }=a^{0,R}.$$

$$\Vert a^{0,R,R,+\infty }-a^0{\Vert }_F\le \frac{C}{R}.$$

$$\Vert a^{0,R,L,T}-a^0{\Vert }_F\le C\left(L^{-(q+1)}+e^{-\alpha {\pi }^{2}T}+{\left(\frac{R}{\sqrt{T}}+1\right)}^{d-1}\frac{e^{-c\frac{{|R-L|}^2}{T}}}{\left|R-L\right|}+\frac{T^{d+1}}{{\left|R-L\right|}^{2d}}{e}^{-2c\frac{{|R-L|}^2}{T}}\right),$$

$$\Vert a^{0,R,L,T}-a^0{\Vert }_F\le C\left[R^{-(q+1)}+\gamma (R)e^{-\zeta R}\right],$$

$$\zeta =\frac{\pi k_o}{2\sqrt{\beta /\alpha }}\text{ and }\gamma (R)=1+\frac{{\left(\sqrt{R}+1\right)}^{d-1}}{R}+\frac{e^{-\frac{\pi k_o}{2\sqrt{\beta /\alpha }}R}}{R^{d-1}}.$$

$$\Vert b^{0,R,L,T,N}-a^0{\Vert }_F\le C\left(\sqrt{T}{L}^{-(q+1)}+e^{-\frac{\alpha {\pi }^2}{2}T}+\frac{R^{d-1}{T}^{\frac{5-d}{2}}}{{\left|R-L\right|}^3}{e}^{-c\frac{{|R-L|}^2}{T}}+\frac{R^2}{L^{d/2}}{e}^{-\frac{C_d{N}^{2/d}T}{R^2}}\right),$$

$$\Vert b^{0,R,L,T,N}-a^0{\Vert }_F\le C\left[R^{-q-\frac{1}{2}}+\gamma (R)e^{-\zeta R}\right],$$

$$\zeta =\frac{\pi k_o}{\sqrt{8\beta /\alpha }}\text{ and }\gamma (R)=R^{2-d/2}+R^{\frac{d-3}{2}}+1.$$

$$a(y)=\left(\begin{array}{cc}{\left(3+\frac{2\sqrt{17}}{8\sin (2\pi y_1)+9}\right)}^{-1}& 0\\ 0& {\left(\frac{1}{20}+\frac{2\sqrt{17}}{8\cos (2\pi y_2)+9}\right)}^{-1}\end{array}\right).$$

Full text

Résumé Décroissance exponentielle de l’erreur de résonance en homogénéisation numerique via des problèmes de cellules paraboliques et elliptiques. Cette note présent deux nouvelles approches pour trouver les coefficients homogénéisés des EDP elliptiques multi-échelles. Les approches standard pour calculer les coefficients homogénéisés souffrent de ce que l’on appelle l’erreur de résonance, qui découle d’une inadéquation entre les vraies conditions aux limites et celles computationelles. Nos nouvelles méthodes, basées sur des solutions des problèmes de cellules paraboliques et elliptiques, entraînent une décroissance exponentielle de l’erreur de résonance.

1 Introduction

We consider the numerical homogenization of multiscale elliptic partial differential equations (PDEs) of the form

[TABLE]

where $a^{\varepsilon}\in\left[L^{\infty}(\Omega)\right]^{d\times d}$ is symmetric, uniformly elliptic and bounded, and ${\varepsilon\ll|\Omega|^{1/d}=O(1)}$ is the wave-length of the small scale variations in the medium. A direct numerical approximation of $u^{\varepsilon}$ by standard finite element/difference methods is prohibitively expensive as the $\varepsilon$-scale variations need to be resolved on the whole computational domain $\Omega$. Homogenization theory aims at finding an effective coefficient $a^{0}$ (or solution $u^{0}$) such that $-\nabla\cdot\left(a^{0}(x)\nabla u^{0}\right)=f$ describes the coarse-scale behaviour of (1). The coefficient $a^{0}$ (and hence the solution) is no more oscillatory, and a standard solver may be directly applied to the homogenized system once $a^{0}$ is determined. Explicit representations for $a^{0}$ are available only in a few cases, such as periodic microstructures or stationary ergodic random materials. For example, when the medium is such that $a^{\varepsilon}(x)=a(x/\varepsilon)$, and $a$ is a ${K:=[-1/2,1/2]^{d}}$-periodic function, then $a^{0}$ is given by

[TABLE]

where $\chi^{i}$ is a $K$-periodic solution of the so-called corrector problem, see [6, 9, 13]:

[TABLE]

In several situations of interest, e.g. when the period of $a$ is not known or when $a$ is quasi-periodic or stochastic, equations (2) and (3) have to be posed over the whole $\mathbb{R}^{d}$. In this case, we write ${a^{0}=\lim_{R\to+\infty}a^{0,R}}$, where

[TABLE]

the domain $K_{R}:=\left[-\frac{R}{2},\frac{R}{2}\right]^{d}$, and $\psi^{i}_{R}$ solves the Dirichlet problem

[TABLE]

In practice, the value of $a^{0,R}$ can be computed only for finite values of $R$, and an error occurs due to the mismatch on the boundary $\partial K_{R}$ between the values of $\psi^{i}_{R}$ and $\chi^{i}$. This error will then propagate into the domain $K_{R}$ and deteriorate the accuracy of the approximation $a^{0,R}$. It is well-known that if $a$ is $K$-periodic and $R$ is not integer, then ${\left\lVert a^{0,R}-a^{0}\right\rVert_{F}\leq CR^{-1}}$, see [4, 10]. A similar result exists also for stationary ergodic random coefficients, both in continuous [8] and discrete [12] settings. This first order resonance error dominates all other discretization errors in modern multiscale methods and, therefore, better approximation techniques with reduced resonance errors are needed.

In order to reduce the resonance error, previous approaches improved the prefactor (but not the convergence rate) [15], or gave second order rates in $1/R$ [7], or fourth order in the asymptotic limit for large values of $R$ [11]. Another strategy results in arbitrary orders in $1/R$, but at the cost of solving a computationally expensive wave equation [5].

This paper, inspired by [14], presents two strategies based on parabolic and elliptic corrector problems which have exponentially decaying boundary errors at a cost comparable to the one of solving the classical elliptic model.

2 New algorithms for computing the homogenized tensor

It can be seen that the computation of homogenized coefficients is linked to the average of oscillatory functions, as formula (4) shows. A naive averaging of a $K$-periodic function $f$ over $K_{R}$ converges to the mean value ${\frac{1}{|K|}\int_{K}f(y)\;dy}$ with a first order accuracy in $1/R$. To improve this accuracy to arbitrarily high rates, a set of smooth averaging filters can be used, see [11].

Definition 1**.**

We say that a function $\mu:[-1/2,1/2]\mapsto\operatorname{\mathbb{R}}_{+}$ belongs to the space $\mathbb{F}_{q}$, $q\geq 1$, if:

1. i)

$\displaystyle\mu\in C^{q}([-1/2,1/2])\cap W^{q+1,\infty}((-1/2,1/2))$; 2. ii)

$\displaystyle\mu^{(k)}(-1/2)=\mu^{(k)}(1/2)=0,\,\forall k\in\left\{0,\dots,q-1\right\}$; 3. iii)

$\displaystyle\int_{-\frac{1}{2}}^{\frac{1}{2}}\mu\left(y\right)dy=1$.

For $q=0$, we define $\mu\in\mathbb{F}_{0}$ as ${\mu(y)=\mathds{1}_{[-1/2,1/2]}}$, where $\mathds{1}_{I}$ is the characteristic function on the interval $I$.

We say that a function ${\mu_{L}:K_{L}:=[-L/2,L/2]^{d}\subset\mathbb{R}^{d}\to\mathbb{R}_{+}}$, with $L>0$, belongs to the space $\mathbb{F}_{q}(K_{L})$ if $\mu_{L}(y)=\frac{1}{L^{d}}\prod_{i=1}^{d}\mu\left(\frac{y_{i}}{L}\right)$, where $\mu\in\mathbb{F}_{q}$ and $y_{i}$ is the coordinate along the $i$-th direction.

2.1 A parabolic approach

In this section we introduce a numerical homogenization scheme based on the solution of parabolic differential equations, as proposed in the discrete setting [14], since the parabolic Green’s function decays exponentially in space. This yields to a reduced influence of mismatching boundary values on the corrector functions. The new cell problems are defined by

[TABLE]

Then, the homogenized coefficient is approximated by

[TABLE]

As it will be shown in Section 3, the $T$ parameter is crucial to obtain an exponential convergence of the resonance error.

2.2 A modified elliptic approach

The second approach that we propose can be viewed as adding a correction term to the elliptic cell problem (5) to reduce the boundary effect in the interior region $K_{L}$. The new elliptic cell problems are given by

[TABLE]

where

[TABLE]

and $\{\lambda_{k},\varphi_{k}\}_{k=1}^{N}$ are the first $N$ dominant eigenvalues and eigenfunctions of the operator $A:=-\nabla\cdot\left(a(\cdot)\nabla\right)$, equipped with Dirichlet boundary conditions. Moreover, $g^{i}(y):=\nabla\cdot a(y)\mathbf{e}_{i}$, and $g^{i}_{k}:=\langle g^{i},\varphi_{k}\rangle_{L^{2}(K_{R})}$. The homogenized coefficient can then be approximated by

[TABLE]

It is worth mentioning that the correction term $[e^{-A_{N}T}g^{i}](y)$ is an approximation to $[e^{-AT}g^{i}](y)$, which corresponds to the solution of the parabolic PDE (6) at time $T$. However, due to the exponential decay of the semigroup $e^{-AT}$ with respect to the eigenvalues of the operator $A$, one can approximate this correction term with exponential accuracy by computing a few dominant eigenmodes of $A$, instead of solving the full parabolic PDE (6).

3 Main results

3.1 Equivalence between the standard elliptic and the parabolic formulations

In this section, we give a proof of the equivalence between elliptic and parabolic equations, thus legitimating the use of (7) and (9), in place of (4), as upscaling model. In the statement of the results, we will refer to the space $\mathcal{M}(\alpha,\beta,\Omega)$, which consists of symmetric matrices $a\in[L^{\infty}(\Omega)]^{d\times d}$ such that ${\alpha|\zeta|^{2}\leq\zeta\cdot a(y)\zeta\leq\beta|\zeta|^{2}}$, ${\forall\zeta\in\mathbb{R}^{d}}$, ${\text{ a.e. }y\in\Omega\subset\mathbb{R}^{d}}$. We will also use the notation

[TABLE]

Theorem 1**.**

Let $a\in\mathcal{M}(\alpha,\beta,K_{R})$ and let $\nabla\cdot\left(a\mathbf{e}_{k}\right)\in L^{2}(K_{R})$, for $k=1,\dots,d$. Let $u_{R}^{k}\in X_{0}(\operatorname{\mathbb{R}}_{+},K_{R})$ be the unique weak solution of (6) and $\psi^{k}_{R}\in H^{1}_{0}(K_{R})$ be the unique weak solution of (5). Then, for $1\leq j,k\leq d$ the following identities hold

[TABLE]

[TABLE]

Proof.

We reformulate problem (6) as the abstract Cauchy problem in $L^{2}(K_{R})$

[TABLE]

Here, the operator $A:H_{0}^{1}(K_{R})\to H^{-1}(K_{R})$ is defined as $Au:=-\nabla\cdot\left(a\nabla u\right)$. Then, $u^{k}_{R}(t)=e^{-tA}g^{k}$. We know that $\sigma(A)$, the spectrum of $A$, is contained in an open sectorial domain $\alpha+S_{\omega}$, where $\alpha\in\operatorname{\mathbb{R}}$, $\alpha>0$ and

[TABLE]

Then, the Dunford integral representation

[TABLE]

holds, where $\Gamma$ is an infinite curve lying in $\rho(A):=\mathbb{C}\setminus\sigma(A)$ and surrounding $\sigma(A)$ counterclockwise. Then, integrating in time we obtain

[TABLE]

The first equality is given by the Dunford integral formula. The second equality is obtained by Fubini’s theorem. The third equality is true because the double integral is bounded, $\lim_{t\rightarrow+\infty}e^{-tz}=0$ since $Re(z)>0$ on $\Gamma$. The last equality follows from the fact that the function $f(z)=1/z$ is holomorphic in the interior of $\alpha+S_{\omega}$. Since $A$ is an isomorphism and ${\psi^{k}_{R}}$ is the weak solution of ${A\psi^{k}_{R}=g^{k}}$, we have that $A^{-1}g^{k}=\psi^{k}_{R}$ and (10) is proved.

To prove (11), we write the weak formulation of (5) and choose ${\psi^{j}_{R}=\int_{0}^{+\infty}u^{j}_{R}\;dt}$ as test function:

[TABLE]

Using the semigroup property of $e^{-tA}$ and the self-adjointness of $A$ we obtain

[TABLE]

and conclude the proof by the change of variable $t/2\mapsto t$. ∎

By the symmetry of $a$ and the weak form of (5), we can rewrite (4) as:

[TABLE]

Theorem 1 provides an equivalent expression, based on the solutions $u^{i}_{R}$ of the parabolic cell problems (6) over infinite time domain, for the second integral in the above expression. This result is summarized as a corollary below.

Corollary 1**.**

Let $a$ satisfy the assumptions of Theorem 1, $a^{0,R}$ be defined by (4) and $a^{0,R,R,+\infty}$ be defined by (7) with $\mu_{R}\in\mathbb{F}_{0}(K_{R})$ (note that $L=R$). Then

[TABLE]

Hence, using the classical result stated in Section 1, there exist a constant $C>0$ independent of $R$ such that

[TABLE]

From this analysis, we can immediately see that, when ${T=+\infty}$, the parabolic approach does not result in any gain in comparison to the standard cell-problem (5), as the two strategies are equivalent and have first order convergence rates in $1/R$.

3.2 Exponential convergence of the parabolic approach (7)

The following Theorem 2 shows that an exponential convergence rate for the boundary error can be attained when the parameter $T$ is sufficiently small and an appropriate filter is used, as fully proved in [2]. If the coefficients $a_{ij}(y)$ are $K$-periodic and (6) is solved with periodic boundary conditions and integer $R$, then Theorem 1 still holds true by substituing functions $\psi^{i}_{R}$ in (10) and (11) with $\chi^{j}$ defined in (3). Thus, it is possible to find an equivalent formula for the exact homogenized coefficients, which is based on a parabolic model with periodic boundary conditions. The proof of Theorem 2 is based on such an equivalence result and on a decomposition of the resonance error in several terms, respectively accounting for the averaging error of periodic functions, the boundary mismatch between the two problems and the truncation in time.

Theorem 2**.**

Let $a\in\mathcal{M}(\alpha,\beta,K_{R})$ be $K$-periodic, $\nabla\cdot\left(a\mathbf{e}_{i}\right)\in L^{2}(K_{R})$ for any $i=1,\dots,d$, $\mu_{L}\in\mathbb{F}_{q}(K_{L})$, for $0<L<\lfloor R\rfloor$, $R>1$ and $T>0$. Then

[TABLE]

where $C>0$ is a constant independent of $R,L,T$ and $c=1/4\beta$. Moreover, the choices $L=(1-k_{o})R$ and $T=k_{T}R$, with $0<k_{o}<1$ and $k_{T}=\frac{k_{o}}{\pi\sqrt{4\beta\alpha}}$ result in the following convergence rate in terms of $R$:

[TABLE]

with

[TABLE]

The term $L^{-(q+1)}$ is the averaging error induced by using a filter function $\mu_{L}\in\mathbb{F}_{q}(K_{L})$, and it can be made arbitrarily small by taking higher values for $q$. The term $e^{-\alpha\pi^{2}T}$ originates from using a finite $T$ for the parabolic cell problem (6). The remaining terms are the errors due to the boundary conditions, which decay exponentially provided $T<|R-L|^{2}$. Moreover, the quasi-optimal scaling of $L$ and $T$ in terms of $R$ are found by equating the exponents of the truncation and boundary errors. Note that bound (12) is similar to the one obtained in [11], except for the term $T^{-2}$ that accounts for the effect of using a biased model equation.

3.3 Exponential convergence of the modified elliptic approach (9)

The boundary error associated with formula (9) relies on how the parameter $T$ is tuned, similarly to the parabolic case. If $T=+\infty$, then the term $e^{-A_{N}T}g^{j}$ vanishes, and the standard Dirichlet cell problem (5) is recovered. Hence, no improvement over the first order convergence rate will be observed. In the following theorem, we specify the precise exponential upper bound for the modified elliptic approach. The proof is based on the equality ${\lim\limits_{N\to\infty}\chi^{i}_{R,T,N}=\int_{0}^{T}u^{i}_{R}\,dt}$ and it is developed in [3].

Theorem 3**.**

Let $a\in\mathcal{M}(\alpha,\beta,K_{R})$ be $K$-periodic, $\nabla\cdot\left(a\mathbf{e}_{i}\right)\in L^{2}(K_{R})$ for any $i=1,\dots,d$, $\mu_{L}\in\mathbb{F}_{q}(K_{L})$, for $0<L<\lfloor R\rfloor$, $R>1$ and $T>0$. Then

[TABLE]

where $C>0$ is a constant independent of $R,L$ and $T$ but may depend on $a$ and $\mu$, $c=1/4\beta$ and $c_{d}>0$ is a constant (independent of $R,T,L$) that may depend on the dimension $d$, $\alpha$ and $\beta$. Moreover, the choices $L=(1-k_{o})R$, $T=k_{T}R$ and $N=\left\lfloor{\left(\frac{\alpha\pi^{2}}{2c_{d}}\right)^{d/2}R^{d}}\right\rfloor$, with $0<k_{o}<1$ and $k_{T}=\frac{k_{o}}{\pi\sqrt{2\beta\alpha}}$ result in the following convergence rate in terms of $R$:

[TABLE]

with

[TABLE]

The upper bounds in Theorem 3 have a similar character to those in Theorem 2, except the error $e^{-\frac{c_{d}N^{2/d}T}{R^{2}}}$ which comes from the spectral truncation. In particular, Theorem 3 shows that an exponential convergence (for the spectral error) will be achieved if the number of modes scales as ${N=O(R^{d})}$ and ${T=O(R)}$. In practice, small values of $R$, e.g. $R=10$, are preferred for simulations. This makes the elliptic approach very favourable both from a computation and accuracy point of view.

4 Numerical validation

Here, we show the results of numerical tests performed using a two-dimensional periodic tensor for validating the convergence rates of Theorems 2 and 3 (note that both theorems assume periodicity of the coefficient). In particular, we consider the following $2\times 2$ tensor

[TABLE]

We compute a numerical approximation of $a^{0,R,L,T}$ and $b^{0,R,L,T,N}$ for many values of $R$ and with the optimal values for $L$ and $T$ (as expressed in Theorems 2 and 3). The error between the numerical approximations and the exact value $a^{0}$ is then plotted against $R$, see Figure 1. The reference value for $a^{0}$ is computed by solving the standard elliptic corrector problems (3) with $R=1$, periodic boundary conditions and by using formula (2). Numerical approximations for the parabolic formulation (6) are computed by a P1-Finite Elements discretization with meshsize $h=1/100$ in space, while a Rosenbrock formula of order 2 with tolerance $tol=10^{-5}$ and adaptive stepping scheme is used in time. For the modified elliptic approach we used a meshsize $h=1/160$ and $N=60$ eigenmodes for approximating the right-hand side $e^{-A_{N}T}g^{j}$, in disregard of the size $R$.

The decay of the overall upscaling error is pictured in Figure 1. In particular, for relatively low to moderate values of $R$, e.g., $1<R\leq 10$, the exponentially decaying boundary error is negligible in comparison to the averaging error when $q=1$ and $q=3$.

Besides the improved convergence rate, it is desirable that the computational cost of the proposed methods is comparable to the one of the classical model, that equals the cost of solving $d$ linear systems. Those can be solved by different numerical schemes, from LU decomposition, more suitable for smaller and two dimensional problems, to iterative methods like GMRES or CG, that are more indicated for large systems coming from three dimensional models. The first of the two proposed approaches, the parabolic model, can be efficiently solved by stabilized explicit ODE solver (such as RKC2 [16] or ROCK4 [1]), whose algorithms perform cheap matrix-vector multiplications iteratively. The number of iterations depends on the number of time steps and stages, but not on the dimension of the system. Lastly, the modified elliptic method has the same cost as the classical model, with the additional expense of accurately reconstructing the modified right-hand side by eigenfunction decomposition of the operator $A$. The eigenmodes computation can be done, for instance, by Krylov-Schur decomposition. A full analysis of the computational cost for the two methods will be addressed in future studies.