Effective models and numerical homogenization for wave propagation in heterogeneous media on arbitrary timescales

TL;DR

This paper introduces a family of effective equations for wave propagation in periodic media valid over arbitrary timescales, ensuring well-posedness without regularization and providing a computationally efficient method for calculating effective tensors.

Contribution

The paper develops a new family of effective wave equations for arbitrary timescales, proving their well-posedness and introducing a novel, cost-effective numerical method for computing effective tensors of any order.

Findings

Effective equations are well-posed without regularization.

Effective solutions are close to original wave solutions in a specific norm.

A new relation reduces computational cost for effective tensor calculation.

Abstract

A family of effective equations for wave propagation in periodic media for arbitrary timescales $\mathcal{O}(\varepsilon^{-\alpha})$, where $\varepsilon\ll1$ is the period of the tensor describing the medium, is proposed. The well-posedness of the effective equations of the family is ensured without requiring a regularization process as in previous models [A. Benoit and A. Gloria, arXiv:1701.08600, 2017], [G. Allaire, A. Lamacz, and J. Rauch, arXiv:1803.09455, 2018]. The effective solutions in the family are proved to be $\varepsilon$ close to the original wave in a norm equivalent to the $L^\infty(0,\varepsilon^{-\alpha}T;L^2(\Omega))$ norm. In addition, a numerical procedure for the computation of the effective tensors of arbitrary order is provided. In particular, we present a new relation between the correctors of arbitrary order, which allows to substantially reduce the computational cost of the effective tensors of arbitrary order. This relation is not limited to the effective equations presented in this paper and can be used to compute the effective tensors of alternative effective models.

Equations

$${\partial^{2}_{t}}u^{\varepsilon}(t,x)-\nabla_{\makebox[4.00069pt]{$\scriptstyle x$}}\kern-1.5pt\cdot\big{(}a^{\varepsilon}(x)\nabla_{\makebox[4.00069pt]{$\scriptstyle x$}}\kern-1.5ptu^{\varepsilon}(t,x)\big{)}=f(t,x)\quad\text{in }(0,T]\times\Omega,$$

$${\partial^{2}_{t}}u^{0}(t,x)-\nabla_{\makebox[4.00069pt]{$\scriptstyle x$}}\kern-1.5pt\cdot\big{(}a^{0}(x)\nabla_{\makebox[4.00069pt]{$\scriptstyle x$}}\kern-1.5ptu^{0}(t,x)\big{)}=f(t,x)\quad\text{in }(0,T]\times\Omega,$$

$$a^{\varepsilon}(x)=a\big{(}\tfrac{x}{\varepsilon}\big{)},\qquad\text{where }a(y)\text{ is $Y$-periodic},$$

$${\partial^{2}_{t}}\tilde{u}-a^{0}\nabla^{2}_{x}\tilde{u}-\displaystyle\sum_{r=1}^{\lfloor\alpha/2\rfloor}(-1)^{r}\varepsilon^{2r}\big{(}a^{2r}\nabla^{2r+2}_{x}\tilde{u}-b^{2r}\nabla^{2r}_{x}{\partial^{2}_{t}}\tilde{u}\big{)}=f+\displaystyle\sum_{r=1}^{\lfloor\alpha/2\rfloor}(-1)^{r}\varepsilon^{2r}b^{2r}\nabla^{2r}_{x}f\quad\text{in }(0,\varepsilon^{-\alpha}T]\times\Omega,$$

$$\|w\|_{W}=\inf_{\begin{subarray}{c}w=w_{1}+w_{2}\\ w_{1},w_{2}\in{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)\end{subarray}}\Big{\{}\|w_{1}\|_{{\mathrm{L}^{2}}(\Omega)}+\|\nabla w_{2}\|_{{\mathrm{L}^{2}}(\Omega)}\Big{\}}\quad\forall w\in{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega).$$

$$\big{(}S^{n}(q)\big{)}_{i_{1}\cdots i_{n}}=\frac{1}{n!}\sum_{\sigma\in\mathbb{S}_{n}}q_{i_{\sigma(1)}\cdots i_{\sigma(n)}}.$$

$$p{=}_{S}q\iff S^n(p)=S^n(q).$$

$$q_{i_1\cdots i_n{i}_{n+1}\cdots i_{2n}}=q_{i_{n+1}\cdots i_{2n}{i}_1\cdots i_n}1\le i_1\cdots i_{2n}\le d.$$

$$q\xi :\xi =q_{i_1\cdots i_{2n}}{\xi }_{i_1\cdots i_n}{\xi }_{i_{n+1}\cdots i_{2n}}\ge 0\forall \xi \in {\mathrm{Sym}}^n({\mathbb{R}}^d),$$

$${\otimes }^{s}q=\underset{s\text{ times}}{\underbrace{q\otimes \cdots \otimes q}}.$$

$$q{\nabla }_x^n:=q_{i_1\cdots i_n}{\partial }_{i_1\cdots i_n}^n.$$

$$\lambda |\xi {|}^2\le a(y)\xi \cdot \xi \le \Lambda |\xi {|}^2\forall \xi \in {\mathbb{R}}^d\text{ for a.e. }y\in Y.$$

$$n_i=\frac{{\omega }_i^r-{\omega }_i^l}{{\ell }_i\epsilon }\in {\mathbb{N}}_{>0}\forall i=1,\dots ,d.$$

$$\begin{array}[]{ll}{\partial^{2}_{t}}u^{\varepsilon}(t,x)-\nabla_{\makebox[2.58612pt]{$\scriptstyle x$}}\kern-1.5pt\cdot\big{(}a\big{(}\tfrac{x}{\varepsilon}\big{)}\nabla_{\makebox[2.58612pt]{$\scriptstyle x$}}\kern-1.5ptu^{\varepsilon}(t,x)\big{)}=f(t,x)&\text{in }(0,\varepsilon^{-\alpha}T]\times\Omega,\\[1.0pt] x\mapsto u^{\varepsilon}(t,x)~{}~{}\Omega\text{-periodic}&\text{in }[0,\varepsilon^{-\alpha}T],\\[1.0pt] u^{\varepsilon}(0,x)=u_{0}(x),\quad\partial_{t}u^{\varepsilon}(0,x)=u_{1}(x)&\text{in }\Omega,\end{array}$$

$${\partial^{2}_{t}}\tilde{u}-a^{0}\nabla^{2}_{x}\tilde{u}-\displaystyle\sum_{r=1}^{\lfloor\alpha/2\rfloor}(-1)^{r}\varepsilon^{2r}\Big{(}a^{2r}\nabla^{2r+2}_{x}\tilde{u}-b^{2r}\nabla^{2r}_{x}{\partial^{2}_{t}}\tilde{u}\Big{)}=Qf\quad\text{in }(0,\varepsilon^{-\alpha}T]\times\Omega,$$

$$\displaystyle\big{(}v,w\big{)}_{\mathcal{H}}$$

$$\displaystyle A\big{(}v,w\big{)}$$

$$\mathcal{H}=\big{\{}v\in{{\mathrm{L}^{2}_{0}}(\Omega)}\,:\,(v,v)_{\mathcal{H}}<\infty\big{\}},\qquad\mathcal{V}=\big{\{}v\in{{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)}\,:\,(v,v)_{\mathcal{V}}<\infty\big{\}}.$$

$$\displaystyle\scalebox{1.2}[1.6]{$\int$}\raisebox{-6.0pt}{\makebox[3.5pt]{$\scriptstyle 0$}}\raisebox{10.0pt}{\makebox[6.92812pt]{$\scriptstyle{T^{\varepsilon}}$}}\big{(}\tilde{u}(t),\partial_{t}^{2}v(t)\big{)}_{\mathcal{H}}+A\big{(}\tilde{u}(t),v(t)\big{)}\,\mathrm{d}t=$$

$$u_0\in {\mathrm{W}}_{\mathrm{per}}(\Omega )\cap {\mathrm{H}}^{\lfloor \alpha /2\rfloor +1}(\Omega ),u_1\in {\mathrm{L}}_0^2(\Omega )\cap {\mathrm{H}}^{\lfloor \alpha /2\rfloor }(\Omega ),Qf\in {\mathrm{L}}^2(0,T^{\epsilon };{\mathrm{L}}_0^2(\Omega )),$$

$$\mathcal{B}^{\varepsilon}\tilde{u}(t,x)=\tilde{u}(t,x)+\sum_{k=1}^{\alpha+2}\varepsilon^{k}u^{k}\big{(}t,x,\tfrac{x}{\varepsilon}\big{)},$$

$$\Vert u^{\epsilon }-\tilde{u}{\Vert }_{{\mathrm{L}}^{\infty }(0,{\epsilon }^{-\alpha }T;W)}\le \Vert u^{\epsilon }-{\mathcal{B}}^{\epsilon }\tilde{u}{\Vert }_{{\mathrm{L}}^{\infty }(0,{\epsilon }^{-\alpha }T;W)}+\Vert {\mathcal{B}}^{\epsilon }\tilde{u}-\tilde{u}{\Vert }_{{\mathrm{L}}^{\infty }(0,{\epsilon }^{-\alpha }T;W)},$$

$$r^{\varepsilon}=({\partial^{2}_{t}}+\mathcal{A}^{\varepsilon})(\mathcal{B}^{\varepsilon}\tilde{u}-u^{\varepsilon}),\qquad\mathcal{A}^{\varepsilon}v=-\nabla_{\makebox[4.00069pt]{$\scriptstyle x$}}\kern-1.5pt\cdot\big{(}a\big{(}\tfrac{x}{\varepsilon}\big{)}\nabla_{\makebox[4.00069pt]{$\scriptstyle x$}}\kern-1.5ptv\big{)},$$

$$r^{\epsilon }=$$

$$\displaystyle{+}\,\varepsilon^{0\phantom{-}}\big{(}{{\partial^{2}_{t}}\tilde{u}}{\,+\,}{\mathcal{A}_{yy}u^{2}}{\,+\,}{\mathcal{A}_{xy}u^{1}}{\,+\,}{\mathcal{A}_{xx}\tilde{u}}{\,-\,f\,}\big{)}$$

$$\displaystyle{+}\sum_{k=1}^{\alpha}\varepsilon^{k}\big{(}{\partial^{2}_{t}}u^{k}+\mathcal{A}_{yy}u^{k+2}+\mathcal{A}_{xy}u^{k+1}+\mathcal{A}_{xx}u^{k}\,\big{)}~{}~{}+~{}~{}\mathcal{R}^{\varepsilon}_{\mathrm{ini}},$$

$$\displaystyle\mathcal{A}_{yy}v=-\nabla_{\makebox[3.6831pt]{$\scriptstyle y$}}\kern-1.5pt\cdot\big{(}a(y)\nabla_{\makebox[3.6831pt]{$\scriptstyle y$}}\kern-1.5ptv\big{)},\quad\mathcal{A}_{xy}v=-\nabla_{\makebox[3.6831pt]{$\scriptstyle y$}}\kern-1.5pt\cdot\big{(}a(y)\nabla_{x}v\big{)}-\nabla_{x}\cdot\big{(}a(y)\nabla_{\makebox[3.6831pt]{$\scriptstyle y$}}\kern-1.5ptv\big{)},\quad\mathcal{A}_{xx}v=-\nabla_{x}\cdot\big{(}a(y)\nabla_{x}v\big{)},$$

$$\mathcal{R}^{\varepsilon}_{\mathrm{ini}}=\,\varepsilon^{\alpha+1}\big{(}{\partial^{2}_{t}}u^{\alpha+1}+\mathcal{A}_{xy}u^{\alpha+2}+\mathcal{A}_{xx}u^{\alpha+1}\big{)}\,{+}\,\,\varepsilon^{\alpha+2}\big{(}{\partial^{2}_{t}}u^{\alpha+2}+\mathcal{A}_{xx}u^{\alpha+2}\big{)}.$$

$$u^k(t,x,y)={\chi }^k(y){\nabla }_x^k\tilde{u}(t,x)={\chi }_{i_1\cdot \cdot i_k}^k(y){\partial }_{i_1\cdot \cdot i_k}^k\tilde{u}(t,x),$$

$${\partial^{2}_{t}}\tilde{u}=a^{0}\nabla^{2}_{x}\tilde{u}-\varepsilon^{2}\Big{(}a^{2}\nabla^{4}_{x}\tilde{u}-b^{2}\nabla^{2}_{x}{\partial^{2}_{t}}\tilde{u}\Big{)}.$$

$$\displaystyle{\partial^{2}_{t}}\tilde{u}=a^{0}\nabla^{2}_{x}\tilde{u}-\varepsilon^{2}c^{1}\nabla^{4}_{x}\tilde{u}+\mathcal{R}_{2}^{\varepsilon}\tilde{u},\qquad c^{1}=a^{2}-a^{0}\otimes b^{2},\quad\mathcal{R}_{2}^{\varepsilon}\tilde{u}=\varepsilon^{4}\big{(}a^{2}\otimes b^{2}\nabla^{6}_{x}\tilde{u}-b^{2}\otimes b^{2}\nabla^{4}_{x}{\partial^{2}_{t}}\tilde{u}\big{)}.$$

$${\partial }_t^2\tilde{u}+\epsilon {\partial }_t^2{u}^1+{\epsilon }^2{\partial }_t^2{u}^2$$

Full text

Keywords. homogenization, effective equations, wave equation, heterogeneous media, long-time behavior, dispersive waves, a priori error analysis, multiscale method

AMS subject classifications. 35B27, 74Q10, 74Q15, 35L05, 65M60, 65N30

1 Introduction

The wave equation in heterogeneous media is widely used in many applications such as seismic inversion, medical imaging or the manufacture of composite materials. We consider the following model problem: let $\Omega\subset\mathbb{R}^{d}$ be a hypercube and let $u^{\varepsilon}:[0,T]\times\Omega\to\mathbb{R}$ be the solution of

[TABLE]

where we require $x\mapsto u^{\varepsilon}(t,x)$ to be $\Omega$-periodic and the initial conditions $u^{\varepsilon}(0,x)$ and $\partial_{t}u^{\varepsilon}(0,x)$ are given. As we allow the domain $\Omega$ to be arbitrarily large, (1.1) can be used to model wave propagation in infinite media. We assume here that the tensor $a^{\varepsilon}$ varies at the scale $\varepsilon\ll 1$ while the initial conditions and the source $f$ have wavelength of order $\mathcal{O}(1)$. In such multiscale situations, standard numerical methods such as the finite element (FE) method or the finite difference (FD) method are accurate only if the size of the grid resolves the microscopic scale $\mathcal{O}(\varepsilon)$. Hence, as $\varepsilon\to 0$ or as the domain $\Omega$ grows the computational cost of the method becomes prohibitive and multiscale numerical methods are needed.

Several multiscale methods for the approximation of (1.1) are available in the literature. They can be divided into two groups (see [3] for a review). First, the methods suited when the medium does not have scale separation: [33], [27, 26], [34], and [4]. These methods rely on sophisticated finite element spaces relying on the solutions of localized problems at the fine scale. Second, the methods suited when the medium has scale separation (i.e., a special structure of the medium is required). These methods are built in the framework of the heterogeneous multiscale method (HMM): the FD-HMM [23, 11] and the FE-HMM [1]. In both methods, the effective behavior of the wave is approximated by solving micro problems in small sampling domains.

The FD-HMM and the FE-HMM rely on homogenization theory [14, 36, 12, 28, 17, 32]: they are built to approximate the homogenized equation and thus provide approximations of $u^{\varepsilon}$ in an ${\mathrm{L}^{\infty}}(0,T;{{\mathrm{L}^{2}}(\Omega)})$ sense. The homogenization of the wave equation (1.1) is provided in [15]. For a given sequence of tensors $\{a^{\varepsilon}\}_{\varepsilon>0}$, we have the existence of a subsequence of $\{u^{\varepsilon}\}_{\varepsilon>0}$ that converges weakly∗ in ${\mathrm{L}^{\infty}}(0,T;{{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)})$ to $u^{0}$ as $\varepsilon\to 0$ (definitions of the functional spaces are provided below). The limit $u^{0}$, called the homogenized solution, solves the homogenized equation

[TABLE]

with the same initial conditions as for $u^{\varepsilon}$. The homogenized tensor $a^{0}$ in (1.2) is obtained as the $G$-limit of a subsequence of $\{a^{\varepsilon}\}_{\varepsilon>0}$ (see [38, 19]). In general, $a^{0}$ depends on the choice of the subsequence and thus no formula is available for its computation. In this paper, we consider periodic media, i.e., we assume that the medium is described by

[TABLE]

where $Y$ is a reference cell (typically $Y=(0,1)^{d}$). Under assumption (1.3), $a^{0}$ is proved to be constant and an explicit formula is obtained (see e.g. [14, 12, 28, 17]): it can be computed by means of the first order correctors, which are defined as the solutions of cell problems (i.e., elliptic equations in $Y$ based on $a(y)$ with periodic boundary conditions). Therefore, in the periodic case the homogenized solution $u^{0}$ can be accurately approximated independently of $\varepsilon$.

However, for wave propagation on large timescales, $u^{\varepsilon}$ develops dispersive effects at the macroscopic scale that are not captured by $u^{0}$. Furthermore, if the initial conditions or the source have high spatial frequencies (in between $\mathcal{O}(1)$ and $\mathcal{O}(\varepsilon)$), the dispersion appears at shorter times. Hence, to develop numerical homogenization methods for long-time propagation, or in high frequency regimes, new effective models are required.

The study of this dispersion phenomenon has recently been the subject of considerable interest. Analyses for periodic media and timescales $\mathcal{O}(\varepsilon^{-2})$ are provided in [37, 29, 20, 21, 6, 8, 5] and numerical approaches are studied in [2, 10]. A result for locally periodic media for timescales $\mathcal{O}(\varepsilon^{-2})$ was also obtained in [7]. For arbitrary timescales $\mathcal{O}(\varepsilon^{-\alpha})$, $\alpha\in\mathbb{N}$, effective equations were proposed in [9] and [13]. The well-posedness of these equations is obtained using regularization techniques: in [13], the regularization relies on the tuning of an unknown parameter, which poses problems in practice; in [9] a filtering process is introduced (yet not tested in practice).

In this paper, we present two main results first reported in [35, Chap. 5]. The first main result is the definition of a family of effective equations that approximate $u^{\varepsilon}$ for arbitrary timescales $\mathcal{O}(\varepsilon^{-\alpha})$. The effective equations, derived by generalizing the technique introduced for timescales $\mathcal{O}(\varepsilon^{-2})$ in [5], have the form 111In the whole paper we use the shorthand $q\nabla^{n}_{x}v$ to denote the operator $q_{i_{1}\cdot\cdot i_{n}}\partial^{n}_{i_{1}\cdot\cdot i_{n}}v$, see (1.10).

[TABLE]

where $a^{0}$ is the homogenized tensor and $a^{2r},b^{2r}$ are pairs of nonnegative, symmetric tensors of order $2r+2$ and $2r$, respectively, which satisfy constraints based on high order correctors, solutions of cell problems. Note that the correction of the right-hand side generalizes the one introduced in the case $\alpha=2$ in [8] and discussed in [5]. For all effective solutions $\tilde{u}$ in the family, we prove an error estimate that ensures $\tilde{u}$ to be $\varepsilon$ close to $u^{\varepsilon}$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;W)$ norm (see (1.5)). In contrast to the effective equations proposed in [9] and [13], the well-posedness of (1.4) does not rely on regularization but is naturally ensured by the non-negativity of the tensors. The unregularized versions of the effective equations from [9] and [13] do not belong to the family (1.4) but are closely related: they have the form (1.4) with $a^{2r}=g^{2r}$, $b^{2r}=0$ and these pairs of tensors satisfy the correct constraints. The issue is however that the sign of $g^{2r}$ happens to be negative for some $r$ (this is proved for $g^{2}$ [18]) which results in the ill-posedness of the unregularized version of these effective equations.

The second main result of the paper is an explicit procedure for the computation of the high order effective tensors $\{a^{2r},b^{2r}\}$ in (1.4), for which we provide a new relation between the high order correctors. In particular, while the natural formula to compute $a^{2r},b^{2r}$ requires to solve the cell problems of order $1$ to $2r+1$, this relation ensures that only the cell problems of order $1$ to $r+1$ are in fact necessary. The consequence is a significant reduction of the computational cost needed to compute the effective tensors of arbitrary order. We emphasize that this result can also directly be used to reduce the computational cost for the tensors of the effective equations from [9] and [13].

The paper is organized as follows. In Section 2, we present our first main result: we derive the family of effective equations and state the error estimate. We then compare the obtained effective equations with the ones from [9] and [13]. In Section 3, we construct a numerical procedure to compute the tensors of effective equations. We then present our second main result: a relation between the correctors which allows to reduce the computational cost of the effective tensors. In Section 4, we illustrate our theoretical findings in various numerical experiments. Finally, in Section 5 we provide the proofs of the main results.

Definitions and notation

Let us start by introducing some definitions and notations used in the paper. Let ${\mathrm{H}^{1}_{\mathrm{per}}}(\Omega)$ be the closure of the space $\mathcal{C}^{\infty}_{\mathrm{per}}(\Omega)$ for the ${\mathrm{H}^{1}}$ norm. We denote the quotient spaces ${\mathcal{L}^{2}}(\Omega)={\mathrm{L}^{2}}(\Omega)/\mathbb{R}$ and ${\mathcal{W}_{{\mathrm{per}}}}(\Omega)={\mathrm{H}^{1}_{\mathrm{per}}}(\Omega)/\mathbb{R}$. The space ${\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)$ (resp. ${{\mathrm{L}^{2}_{0}}(\Omega)}$) is composed of the zero mean representatives of the equivalence classes in ${\mathcal{W}_{{\mathrm{per}}}}(\Omega)$ (resp. ${\mathcal{L}^{2}}(\Omega)$). The dual space of ${\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)$ (resp. ${\mathcal{W}_{{\mathrm{per}}}}(\Omega)$) is denoted ${\mathrm{W}_{\!{\mathrm{per}}}^{*}}(\Omega)$ (resp. ${\mathcal{W}_{{\mathrm{per}}}^{*}}(\Omega)$). The integral mean of $v\in{\mathrm{L}^{1}}(\Omega)$ is denoted $\langle{v}\rangle_{\Omega}=\frac{1}{|\Omega|}\int_{\Omega}v$ and $(\cdot,\cdot)_{\Omega}$ denotes the standard inner product in ${\mathrm{L}^{2}}(\Omega)$. We define the following norm on ${\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)$

[TABLE]

Using the Poincaré–Wirtinger inequality, we verify that $\|\cdot\|_{W}$ is equivalent to the ${\mathrm{L}^{2}}$ norm: $\|w\|_{W}\leq\|w\|_{{\mathrm{L}^{2}}(\Omega)}\leq\max\{1,C_{\Omega}\}\|w\|_{W}$ where $C_{\Omega}$ is the Poincaré constant.

We denote ${\mathrm{Ten}}^{n}(\mathbb{R}^{d})$ the vector space of tensors of order $n$. In the whole text, we drop the notation of the sum symbol for the dot product between two tensors and use the convention that repeated indices are summed. The subspace of ${\mathrm{Ten}}^{n}(\mathbb{R}^{d})$ of symmetric tensors is denoted ${\mathrm{Sym}}^{n}(\mathbb{R}^{d})$, i.e., $q\in{\mathrm{Sym}}^{n}(\mathbb{R}^{d})$ iff $q_{i_{1}\cdots i_{n}}=q_{i_{\sigma(1)}\cdots i_{\sigma(n)}}$ for any permutation of order $n$ $\sigma\in\mathbb{S}_{n}$. We define the symmetrization operator $S^{n}:{\mathrm{Ten}}^{n}(\mathbb{R}^{d})\to{\mathrm{Sym}}^{n}(\mathbb{R}^{d})$ as

[TABLE]

The coordinate \big{(}S^{n}(q)\big{)}_{i_{1}\cdots i_{n}} is denoted $S^{n}_{i_{1}\cdots i_{n}}\{q_{i_{1}\cdots i_{n}}\}$. We denote $=_{S}$ an equality holding up to symmetries, i.e., for $p,q\in{\mathrm{Ten}}^{n}(\mathbb{R}^{d})$ we have

[TABLE]

A colon is used to denote the inner product of two tensors in ${\mathrm{Ten}}^{n}(\mathbb{R}^{d})$, $p:q=p_{i_{1}\cdots i_{n}}q_{i_{1}\cdots i_{n}}$. We say that a tensor $q\in{\mathrm{Ten}}^{2n}(\mathbb{R}^{d})$ is major symmetric if it satisfies

[TABLE]

We say that a tensor $q\in{\mathrm{Ten}}^{2n}(\mathbb{R}^{d})$ is positive semidefinite if

[TABLE]

and it is positive definite if the equality in (1.9) holds only for $\xi=0$. The tensor product of $p\in{\mathrm{Ten}}^{m}(\mathbb{R}^{d})$ and $q\in{\mathrm{Ten}}^{n}(\mathbb{R}^{d})$ is the tensor of ${\mathrm{Ten}}^{m+n}(\mathbb{R}^{d})$ defined as $(p\otimes q)_{i_{1}\cdots i_{m+n}}=p_{i_{1}\cdots i_{m}}q_{i_{m+1}\cdots i_{m+n}}$. Note that up to symmetries the tensor product is commutative, i.e., $p\otimes q=_{S}q\otimes p$. We use the shorthand notation

[TABLE]

The derivative with respect to the $i$-th space variable $x_{i}$ is denoted $\partial_{i}$ and the derivation with respect to any other variable is specified. For $q\in{\mathrm{Ten}}^{n}(\mathbb{R}^{d})$, we denote the differential operator

[TABLE]

Settings of the problem

Recall assumption (1.3): a^{\varepsilon}(x)=a\big{(}\tfrac{x}{\varepsilon}\big{)}, where $y\mapsto a(y)$ is a $d\times d$ symmetric, $Y$-periodic tensor. In addition, we assume that $a(y)$ is uniformly elliptic and bounded, i.e., there exists $\lambda,\Lambda>0$ such that

[TABLE]

Without loss of generality, let the reference cell be $Y=(0,\ell_{1})\times\cdots\times(0,\ell_{d})$. We assume that the hypercube $\Omega=(\omega^{l}_{1},\omega^{r}_{1})\times\cdots\times(\omega^{l}_{d},\omega^{r}_{d})$ satisfies

[TABLE]

In particular, (1.12) ensures that for a $Y$-periodic function $\gamma$, the map x\mapsto\gamma\big{(}\tfrac{x}{\varepsilon}\big{)} is $\Omega$-periodic ($\gamma$ is extended to $\mathbb{R}^{d}$ by periodicity). Note that the integers $n_{i}$ in (1.12) can be arbitrarily large. In particular $n_{i}$ can be of order $\mathcal{O}(\varepsilon^{-\alpha})$.

Given an integer $\alpha\geq 0$, we consider the wave equation: $u^{\varepsilon}:[0,\varepsilon^{-\alpha}T]\times\Omega\rightarrow\mathbb{R}$ such that

[TABLE]

where $u_{0},u_{1}$ are given initial conditions and $f$ is a source. The well-posedness of (1.13) is proved in [31, 24]: if $u_{0}\in{{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)}$, $u_{1}\in{{\mathrm{L}^{2}_{0}}(\Omega)}$, and $f\in{\mathrm{L}^{2}}(0,{T^{\varepsilon}};{{\mathrm{L}^{2}_{0}}(\Omega)})$, then there exists a unique weak solution $u^{\varepsilon}\in{\mathrm{L}^{\infty}}(0,{T^{\varepsilon}};{{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)})$ with $\partial_{t}u^{\varepsilon}\in{\mathrm{L}^{\infty}}(0,{T^{\varepsilon}};{{\mathrm{L}^{2}_{0}}(\Omega)})$ and ${\partial^{2}_{t}}u^{\varepsilon}\in{\mathrm{L}^{2}}(0,{T^{\varepsilon}};{{\mathrm{W}_{\!{\mathrm{per}}}^{*}}(\Omega)})$.

2 First main result: family of effective equations and a priori error estimate

In this section, we present the family of effective equations and provide the corresponding a priori error estimate. In 2.1, we derive the family in three steps: (i) we discuss the ansatz on the form of the effective equations; (ii) using asymptotic expansion we derive the high order cell problems; (iii) we obtain the constraints on the effective tensors by investigating the well-posedness of the cell problems. In 2.2, we define rigorously the family of effective equations and state the a priori error estimate. Finally, in 2.3 we compare the obtained equations with the other effective equations available in the literature. For the sake of readability we postpone the technical proofs to Section 5.

2.1 Derivation of the family of effective equations

In the whole derivation, we assume that the data are as regular as necessary. The specific requirements are stated in Theorem 2.8. Note that we consider here timescales $\varepsilon^{-\alpha}T$ with $\alpha\geq 2$. For timescales $\varepsilon^{-\alpha}T$ with $\alpha<2$ it can be shown following similar techniques that the standard homogenized equation is a valid effective model (see [35, section 5.1.1]).

Ansatz on the form of the effective equations

We first discuss the ansatz on the form of the effective equations, which has a major importance in the derivation. We assume that the effective equations have the form

[TABLE]

where $a^{0}$ is the homogenized tensor (2.13), $a^{2r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$, $b^{2r}\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$ are tensors to be defined and $Q$ is a differential operator to be defined (the construction of $Qf$ is discussed in Remark 5.2). As discussed in [5], if the set of considered equations is too small, we end up with ill-posed equations. In particular, without the operators $-b^{2r}\nabla^{2r}_{x}{\partial^{2}_{t}}$ in (2.1), our derivation would lead to the same ill-posed effective equations obtained in [9] and [13] (the unregularized versions).

Following the classical Faedo–Galerkin method (see [24, Chap. 7]), we prove the following well-posedness result for (2.1). We define the bilinear forms

[TABLE]

and the associated Banach spaces

[TABLE]

We call a function $\tilde{u}\in{\mathrm{L}^{\infty}}(0,{T^{\varepsilon}};\mathcal{V})$, with $\partial_{t}\tilde{u}\in{\mathrm{L}^{\infty}}(0,{T^{\varepsilon}};\mathcal{H})$, a weak solution of (2.1) if for all test functions $v\in\mathcal{C}^{2}([0,{T^{\varepsilon}}];\mathcal{V})$, with $v({T^{\varepsilon}})=\partial_{t}v({T^{\varepsilon}})=0$, $\tilde{u}$ satisfies

[TABLE]

Theorem 2.1**.**

Assume that the tensors $a^{2r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$, $b^{2r}\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$ are positive semidefinite (1.9) and satisfy the major symmetries (1.8). Furthermore, assume that the data satisfy the regularity

[TABLE]

Then there exists a unique weak solution of (2.1).

Let us provide a short sketch of the proof. We look for successive approximations of a weak solution in the form $u^{m}(t)=\sum_{k=0}^{m}u^{m}_{k}(t)\varphi_{k}$, where $\{\varphi_{k}\}_{k\in\mathbb{N}}$ is a smooth basis of ${{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)}$. For each $m$, $u^{m}(t)$ is obtained as the solution of a well-posed ordinary differential equation. We then prove that the sequence $\{u^{m}\}_{m\geq 0}$ is bounded in ${\mathrm{L}^{\infty}}(0,{T^{\varepsilon}};\mathcal{V})$. In particular, note that the sign assumptions on the tensors ensure that $(v,v)_{\mathcal{H}}\geq\|v\|_{\mathrm{L}^{2}}^{2}$ for any $v\in\mathcal{H}$ and and $A(v,v)\geq\lambda\|\nabla v\|_{\mathrm{L}^{2}}^{2}$ for all $v\in\mathcal{V}$. We thus obtain the existence of a subsequence that weakly∗ converges in ${\mathrm{L}^{\infty}}(0,{T^{\varepsilon}};\mathcal{V})$. We can then prove that the weak∗ limit is the unique weak solution.

Asymptotic expansion, inductive Boussinesq tricks

We make the ansatz that $u^{\varepsilon}$ can be approximated by an adaptation of $\tilde{u}$ of the form

[TABLE]

where $u^{k}$ are to be defined and the map $y\mapsto u^{k}(t,x,y)$ is $Y$-periodic. We split the error as

[TABLE]

and follow the argument presented in [5] based on the error estimate of Lemma 5.8 (see also [35, Section 4.2.2]): for $\|u^{\varepsilon}-\mathcal{B}^{\varepsilon}\tilde{u}\|_{{\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;W)}$ to be of order $\mathcal{O}(\varepsilon)$, we need the terms involving $\tilde{u}$ in the remainder

[TABLE]

to be of order $\mathcal{O}(\varepsilon^{\alpha+1})$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathrm{W}_{\!{\mathrm{per}}}^{*}}(\Omega)})$ norm (it is sufficient that this holds in the stronger ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathrm{L}^{2}}(\Omega)})$ norm). We now expand $r^{\varepsilon}$: using the equation for $u^{\varepsilon}$ (1.13) and the form of the adaptation (2.3), we obtain

[TABLE]

where the operators $\mathcal{A}_{yy},\mathcal{A}_{xy},\mathcal{A}_{xx}$ are defined as

[TABLE]

and the remainder is

[TABLE]

Classical two scale asymptotic expansion [14] advises to look for $u^{k}$ of the form

[TABLE]

where the components of the tensor $\chi^{k}$ are $Y$-periodic functions to be defined.

The next step is the main difficulty of the derivation: we must use the effective equation (2.1) to substitute all the time derivatives in the terms of order $\mathcal{O}(\varepsilon^{0})$ to $\mathcal{O}(\varepsilon^{\alpha})$ in (2.5). Various versions of this process have been used in related works [5, 6, 8, 21, 25, 30] to obtain well-posed effective equations from ill-posed ones. To refer to this kind of manipulation the term Boussinesq trick was coined in [8] (see [16] where such tricks are used for several Boussinesq’s type of equations).

As these inductive Boussinesq tricks represent a technical challenge, let us explain here the case $\alpha=2$ and $f=0$ and postpone the general case to Section 5.1.1 (Theorem 5.1 and Lemma 5.7). For $\alpha=2$ and $f=0$, the effective equation (2.1) can be written as

[TABLE]

We now use the equation to substitute ${\partial^{2}_{t}}\tilde{u}$ in the last term and obtain

[TABLE]

Using (2.7) and the two last equalities, we find that the time derivatives in (2.5) can be written as

[TABLE]

where the remainder $\mathcal{R}^{\varepsilon}\tilde{u}$ has order $\mathcal{O}(\varepsilon^{3})$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathrm{L}^{2}}(\Omega)})$ norm (for $\tilde{u}$ sufficiently regular). Using this expression in (2.5), we obtain the desired development in the case $\alpha=2$ and $f=0$.

This process is generalized in Theorem 5.1 and Lemma 5.7. With these results, we are able to rewrite $r^{\varepsilon}$ in (2.5) without time derivative in the terms of order $\mathcal{O}(\varepsilon^{0})$ to $\mathcal{O}(\varepsilon^{\alpha})$: defining the tensors $c^{r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$ inductively and $p^{k}\in{\mathrm{Ten}}^{k+2}(\mathbb{R}^{d})$ as

[TABLE]

we obtain

[TABLE]

Provided sufficient regularity of $\tilde{u}$, the remainder $\mathcal{R}^{\varepsilon}\tilde{u}$ has order $\mathcal{O}(\varepsilon^{\alpha+1})$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathrm{L}^{2}}(\Omega)})$ norm. Furthermore, provided sufficient regularity of $f$, the remainder $\mathcal{S}^{\varepsilon}f$ has order $\mathcal{O}(\varepsilon)$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathrm{L}^{2}}(\Omega)})$ norm (this is sufficient for $\tilde{u}$ and $u^{\varepsilon}$ to be $\varepsilon$ close, see Theorem 2.8). Using the definition of $u^{k}$ in (2.7) and of $\mathcal{A}_{yy},\mathcal{A}_{xy},\mathcal{A}_{xx}$, we verify that (for $k=0$ use $u^{0}=\tilde{u}$ and $\chi^{0}=1$)

[TABLE]

where $e_{i}$ denotes the $i$-th vector of the canonical basis of $\mathbb{R}^{d}$. Hence, canceling successively the terms of order $\mathcal{O}(\varepsilon^{-1})$ to $\mathcal{O}(\varepsilon^{\alpha})$ in (2.9), we obtain the cell problems: the correctors $\{\chi^{k}_{i_{1}\cdot\cdot i_{k}}\}_{k=1}^{\alpha+2}$ are the functions in ${{\mathrm{W}_{\!{\mathrm{per}}}}(Y)}$ such that

[TABLE]

where the tensors $p^{k}$ are defined in (2.8) and the symmetrization operator is defined in (1.6). Note that the cell problems for $\chi^{1}_{i}$ and $\chi^{2}_{ij}$ are the same as the cell problems obtained with two-scale asymptotic expansion [14].

Remark 2.2*.*

In (2.10), we have chosen symmetric right-hand sides. This is possible thanks to the symmetry of $\nabla^{n}_{x}\tilde{u}$ in (2.9), as a term of the form $q\nabla^{n}_{x}\tilde{u}$ can be rewritten as $S^{n}(q)\nabla^{n}_{x}\tilde{u}$. This choice ensures that the correctors are symmetric tensors functions. In particular, $\chi^{k}$ has only $\binom{k+d-1}{k}$ distinct entries instead of $d^{k}$ if it was not symmetric. When it comes to numerical approximation, each distinct entry of the corrector corresponds to a cell problem to solve and this symmetrization saves computational time. In addition, this choice ensures that the odd order cell problems are well-posed unconditionally (see below).

Constraints on the high order effective tensors

The last step in the derivation of the family of effective equations is to obtain the constraints on the effective tensors by imposing the well-posedness of the cell problems. In order to investigate the solvability of (2.10), let us state the following classical result (obtained with the Fredholm alternative or the Lax–Milgram theorem combined with the characterization of ${{\mathrm{W}_{\!{\mathrm{per}}}^{*}}(Y)}$).

Lemma 2.3**.**

For an elliptic and bounded tensor $a(y)$, consider the following variational problem: find $v\in{{\mathrm{W}_{\!{\mathrm{per}}}}(Y)}$ such that

[TABLE]

where $f^{0},f^{1}_{1},\ldots,f^{1}_{d}$ are given functions. Then (2.11) has a unique solution $v\in{{\mathrm{W}_{\!{\mathrm{per}}}}(Y)}$ if and only if

[TABLE]

We now proceed to the two following tasks: first, we verify that the odd order cell problems satisfy unconditionally the solvability condition (2.12); second, we impose the solvability condition (2.12) on the right-hand sides of the even order cell problems to obtain the constraints on the effective tensors.

First note that (2.10a) is well-posed as its right-hand side unconditionally satisfies (2.12). Next, we verify that the cell problem for $\chi^{2}$ is well-posed as the homogenized tensor satisfies

[TABLE]

which ensures the solvability condition (2.12) to hold. We then continue this process to derive the constraints on the higher order effective tensors imposed by the well-posedness of the higher order cell problems.

We assume that the cell problems are well-posed up to order $2r$, $r\geq 1$. Consider the cell problem for $\chi^{2r+1}$. Recalling that $p^{2r-1}=0$ and as the correctors $\chi^{1},\ldots,\chi^{2r}$ have zero mean, we verify that the solvability condition (2.12) is equivalent to the following relation between the correctors.

Lemma 2.4**.**

For any $1\leq r\leq\lfloor\alpha/2\rfloor$, the correctors $\chi^{2r}$ and $\chi^{2r-1}$ satisfy the equality

[TABLE]

Remark 2.5*.*

A proof of Lemma 2.4 can be found in [35, Lemma 5.2.5]. A similar result is also known in the context of Bloch wave theory (see e.g. [20], [13] and the references therein). As discussed, Lemma 2.4 guarantees that the odd order cell problems are well-posed unconditionally. Pursuing the reasoning, this result ensures that no operator of odd order is needed in the effective equations. It is also the reason why no additional correction is required in the effective equation when increasing the timescale from an even integer $\alpha$ to $\alpha+1$.

Finally, consider the cell problem for $\chi^{2r+2}$. Imposing the solvability condition (2.12) on the right-hand side, we obtain the following constraint on the tensor $p^{2r}$:

[TABLE]

This constraint can be rewritten in terms of the effective tensors $a^{2r},b^{2r}$ using the definition of $c^{r}$ and $p^{2r}$ in (2.8):

[TABLE]

The constraints (2.15) characterize the family of effective equations. Indeed, if the effective tensors $a^{2r},b^{2r}$ satisfy (2.15), the cell problems for $\chi^{1}$ to $\chi^{\alpha+2}$ are well-posed because their right-hand sides satisfy the solvability condition (2.12). Therefore, the adaptation $\mathcal{B}^{\varepsilon}\tilde{u}$ given by (2.3) and (2.7) is well defined and can be used to prove an error estimate between $u^{\varepsilon}-\tilde{u}$. This result is presented in the next section and rigorously proved in Section 5.2.

2.2 A priori error estimate for the family of effective equations

In the previous section, we derived the constraints on the effective tensors for the effective equations to approximate $u^{\varepsilon}$. We present here an error estimate that ensures that the solutions of the derived effective equations are $\varepsilon$ close to $u^{\varepsilon}$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;W)$ norm.

Let us first rigorously define the family of effective equations derived in Section 2.1.

Definition 2.6**.**

The family $\mathcal{E}$ of effective equations is the set of equations (2.1), where $Qf$ is defined as

[TABLE]

and for $1\leq r\leq\lfloor\alpha/2\rfloor$ the tensors $a^{2r},b^{2r}$ satisfy the following requirements:

1. (i)

$a^{2r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$, $b^{2r}\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$; 2. (ii)

$a^{2r}$ and $b^{2r}$ are positive semidefinite (see (1.9)); 3. (iii)

$a^{2r}$ and $b^{2r}$ satisfy the major symmetries (1.8); 4. (iv)

$a^{2r}$ and $b^{2r}$ satisfy the constraints (2.15).

Remark 2.7*.*

In the case $\alpha=2$, the correction $Qf$ of the right-hand side in the effective equations was introduced in [8] and discussed in [5]. We verify that the definition of $Qf$ in (2.16) leads to a lower constant multiplying $\|f\|_{{\mathrm{L}^{1}}(0,\varepsilon^{-\alpha}T;{\mathrm{H}^{r(\alpha)}(\Omega)})}$ in estimate (2.17). More details on the origin of (2.16) are given in Remarks 5.2 and 5.10.

For the effective solutions in family $\mathcal{E}$, we prove in Section 5.2 the following error estimate.

Theorem 2.8**.**

Let $d\leq 3$, $\alpha\geq 2$ and let $\tilde{u}$ belong to the family of effective equations $\mathcal{E}$ (Definition 2.6). Furthermore, assume that $a(y)\in{\mathrm{L}^{\infty}}(Y)$ and that the data and $\tilde{u}$ satisfy the following regularity

[TABLE]

where $r(\alpha)=\alpha+2\lfloor\alpha/2\rfloor+2$. Then the following error estimate holds

[TABLE]

where $C$ depends only on $T$, $\lambda$, $\Lambda$, $\{|a^{2r}|_{\infty},|b^{2r}|_{\infty}\}_{r=1}^{\lfloor\alpha/2\rfloor}$, and $Y$, and the norm $\|{\cdot}\|_{W}$ is defined in (1.5).

Remark 2.9*.*

For $d\geq 4$, the result of Theorem 2.8 holds provided a higher regularity on the effective solution and $f$ are assumed. Specifically, assuming that $m$ is a sufficiently large integer for the embedding $\mathrm{H}^{m}_{\mathrm{per}}(\Omega)\hookrightarrow\mathcal{C}^{0}(\Omega)$ to hold, the statement of Theorem 2.8 is true for $r(\alpha)=\alpha+2\lfloor\alpha/2\rfloor+m$.

Remark 2.10*.*

From the derivation in Section 2.1, we may hope that $\mathcal{B}^{\varepsilon}\tilde{u}$ is a better approximation of $u^{\varepsilon}$ than $\tilde{u}$. For example, in the elliptic case the error between $u^{\varepsilon}$ and an adaptation can be estimated in the energy norm (see e.g. [28], note that under (1.12) there is no boundary layer). However, in the case of the wave equation, the answer is not as simple. In fact, an error estimate in the energy norm can be proved only in some specific settings: for example if $f=0$ and for well prepared initial position $u_{0}$ of the form $\mathcal{B}^{\varepsilon}\bar{u}_{0}$ for some $\bar{u}_{0}$. Note that this issue is related to the lack of convergence of the energy of the fine scale wave toward the energy of the homogenized wave (see [15]).

2.3 Comparison with other effective equations in the literature

In [13], effective equations for arbitrary timescales are derived. The settings are slightly different as the wave equation (1.13) is considered in the whole space $\mathbb{R}^{d}$ and with a tensor that can be of a more general nature: periodic, almost periodic, quasiperiodic and random (we refer to [13] for the specific definitions). In these circumstances, [13] proposes effective equations that have the form

[TABLE]

The effective tensors in this equation are indeed verified to match the tensors $g^{2r}$ defined in (2.14) (see [35, section 5.2.6]). Under some weak regularity assumption on the initial condition, an error estimate between $u^{\varepsilon}$ and $\bar{u}$ is proved in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{\mathrm{L}^{2}}(\mathbb{R}^{d}))$ norm. This estimate is a strong theoretical result as it holds in the norm we expect in the context of homogenization of the wave equation. However, the use of the effective equation (2.18) in practice is problematic as no procedure for the computation of the regularization parameter $\gamma$ is available. In fact, numerical tests indicate that the range of acceptable values for $\gamma$ is narrow: if $\gamma$ is too small, the equation is ill-posed and if $\gamma$ is too large, the solution $\bar{u}$ of (2.18) does not describe $u^{\varepsilon}$ accurately (see [35, section 5.4.3]).

In [9], another effective equation for arbitrary timescales is proposed. The settings of this result are the following: the wave equation (1.13) is considered in the whole space $\mathbb{R}^{d}$ with a periodic tensor, it includes an oscillating density and admits nonzero source and trivial initial conditions. In the particular case of a constant density (as considered in the present paper), the effective equation proposed in [9] reads

[TABLE]

where $S^{\varepsilon}_{1}$, $S^{\varepsilon}_{2}$ are filtering differential operators that ensure the equation to be well-posed. The main result is an error estimate in energy norm between $u^{\varepsilon}$ and an adaptation of $\hat{u}$. In particular, the results enable the adaptation to approximate $u^{\varepsilon}$ as accurately as one wants by increasing $\alpha$ accordingly. Nevertheless, the filtering process used in (2.19) to obtain a well-posed equation has yet to be tested in practice.

Despite the inherent difference between (2.18) and (2.19) and effective equations (2.1) in family $\mathcal{E}$, their respective effective tensors are tied through the relation (see (2.15)):

[TABLE]

Hence, if we let $b^{2r}=0$ for all $r$ in the effective equations in $\mathcal{E}$, relation (2.20) reads $g^{2r}=(-1)^{r}a^{2r}$ and we end up with the unregularized versions of (2.18) and (2.19) (i.e., $\gamma=0$ and $S_{1}^{\varepsilon}=S_{2}^{\varepsilon}=\mathrm{Id}$). Similarly, equations (2.18) and (2.19) without their regularization devices satisfy requirement (iv) in Definition 2.6. However, while the well-posednesses of (2.18) and (2.19) rely on their respective regularization process, the well-posedness of the effective equations in $\mathcal{E}$ is guaranteed by requirements (ii) and (iii). To fulfill these requirements, we have an explicit and constructive algorithm described in Section 3.1.

Finally, before approximating any of the effective equations (2.1), (2.18) or (2.19), the effective tensors $g^{2r}$ must be computed. For this calculation, a substantial gain of computational time is obtained by using the formula provided by our second main result in Theorem 3.5, Section 3.2.

3 Second main result: computation of effective tensors and reduction of the computational cost

In this section, we provide a numerical procedure for the computation of the tensors of some effective equations in the family $\mathcal{E}$. In particular, in 3.2 we present a new relation between the correctors that allows to reduce significantly the computational cost for the effective tensors. We emphasize that this result is not limited to the family $\mathcal{E}$ and can be used to compute the effective tensors of the alternative effective models from [9] and [13] discussed in 2.3. The final algorithm is provided in 3.4.

3.1 Construction of high order effective tensors

Recall that the tensors $\{a^{2r},b^{2r}\}_{r=1}^{\lfloor\alpha/2\rfloor}$ of an effective equation in family $\mathcal{E}$ are characterized by the following properties (Definition 2.6): for $1\leq r\leq\lfloor\alpha/2\rfloor$

1. (i)

$a^{2r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$, $b^{2r}\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$; 2. (ii)

$a^{2r}$ and $b^{2r}$ are positive semidefinite, i.e., (see (1.9))

[TABLE] 3. (iii)

$a^{2r}$ and $b^{2r}$ satisfy the major symmetries, i.e.,

[TABLE] 4. (iv)

$a^{2r}$ and $b^{2r}$ satisfy the constraints

[TABLE]

with the tensors $\check{q}^{r}\in{\mathrm{Sym}}^{2r+2}(\mathbb{R}^{d})$ defined as

[TABLE]

where the tensor $c^{1}\in{\mathrm{Ten}}^{4}(\mathbb{R}^{d}),\ldots,c^{r-1}\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$, defined in (2.8), are computed with $a^{0}$, $\{a^{2s},b^{2s}\}_{s=1}^{r-1}$ and the tensor $g^{2r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$, defined in (2.14), can be computed from the correctors $\chi^{2r}$ and $\chi^{2r+1}$ (a much cheaper alternative is to compute only $S^{2r+2}(g^{2r})$ with formula (3.8), see Section 3.2).

Our goal is now to construct tensors $\{a^{2r},b^{2r}\}_{r=1}^{\lfloor\alpha/2\rfloor}$ satisfying the requirements (i) to (iv). Note that if $\check{q}^{r}$ was positive semidefinite, the pair $a^{2r}=\check{q}^{r}$ and $b^{2r}=0$ would trivially satisfy (i) to (iv). However, although the sign of $\check{q}^{r}$ is unknown for $r\geq 2$, it is known that $\check{q}^{1}$ is negative definite (see [18, 8, 5]). Hence, the main challenge of the construction lies in the sign of the tensors and to build valid tensors, we need two basics from the tensor world. First, we use the following result, proved in Section 5.3.

Lemma 3.1**.**

For any $n\geq 1$, the tensor $S^{2n}(\otimes^{n}a^{0})$ is positive definite.

Second, we use a “matricization” operator, which linearly maps a symmetric tensor $q\in{\mathrm{Sym}}^{n}(\mathbb{R}^{d})$ to a symmetric matrix $M(q)$ whose sign is the same as $q$, i.e.,

[TABLE]

where $\nu$ is a simple tensor-to-vector transformation. One construction for such an operator is provided in Section 3.3.

With Lemma 3.1 and the operator $M$, we are able to build tensors that satisfy (i) to (iv).

Lemma 3.2**.**

For $r\geq 1$, assume that $a^{0}$ and $\{a^{2s},b^{2s}\}_{s=1}^{r-1}$ have already been computed and let $\check{q}^{r}$ be the tensor defined in (3.2). Define the tensor $R\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$ as

[TABLE]

where $\lambda_{\min}(\cdot)$ denotes the minimal eigenvalue of a symmetric matrix and $\{\cdot\}_{+}=\max\{0,\cdot\}$. Then the tensors

[TABLE]

satisfy the requirements (i) to (iv) of Definition 2.6.

Proof.

First, note that the orders of the tensors in (i) are correct. Second, as $a^{2r},b^{2r}$ are fully symmetric tensors, they trivially satisfy the major symmetries and (iii) is verified. Next, we verify (iv), i.e., $a^{2r},b^{2r}$ satisfy (3.1) (recall the meaning of $=_{S}$ in (1.7)):

[TABLE]

We are left with (ii): verifying that $a^{2r},b^{2r}$ are positive semidefinite. The positive semidefiniteness of $b^{2r}=R$ follows directly from the non-negativity of $\delta^{*}$ and Lemma 3.1. Finally, let us verify that $a^{2r}$ is positive semidefinite. Using (3.3) and the definitions in (3.4), we have

[TABLE]

where we denoted $v_{\xi}=\nu(\xi)$ and $|\cdot|$ the Euclidian norm. The definition of $\delta^{*}$ in (3.4) ensures that the right-hand side is nonnegative, proving that $a^{2r}$ is positive semidefinite. That concludes the verification of (i) to (iv) and the proof of the lemma is complete. ∎

Lemma 3.2 allows to compute the tensors of one effective equations in the family $\mathcal{E}$. We emphasize that this is one possible construction among others. To obtain other effective equations, we have many options. For example, note that replacing $\delta^{*}$ in (3.4) with any $\delta\geq\delta^{*}$ provides other valid pairs of tensors. Alternatively, we can replace $S^{2r}(\otimes^{r}a^{0})$ in the definition of $R$ with another positive definite tensor. Finally, let us note the following alternative (used in [5] in the case $\alpha=2$).

Remark 3.3*.*

Assume that the tensor $\check{q}^{r}$ in (3.2) can be decomposed as $\check{q}^{r}=_{S}\check{q}^{r,1}-\check{q}^{r,2}\!\otimes\!a^{0}$, where $\check{q}^{r,2}\in{\mathrm{Ten}}^{2r}(\mathbb{R}^{d})$ is positive semidefinite. Then we verify that the pairs (3.5), where $R$ is defined as

[TABLE]

define effective equations in the family.

3.2 A new remarkable relation between the correctors to reduce the cost of computation of the effective tensors

Let us discuss the cost of the procedure described in the previous section. For an integer $\alpha$, assume that we want to construct the tensors of an effective equation for a timescale $\mathcal{O}(\varepsilon^{-\alpha})$: $a^{0},\{a^{2r},b^{2r}\}_{r=1}^{s}$, where $s=\lfloor\alpha/2\rfloor$. The main computational cost of this construction lies in the calculation of the tensors $\{\check{q}^{r}\}_{r=1}^{s}$ in (3.2), which involves the tensors $\{S^{2r+2}(g^{2r})\}_{r=1}^{s}$, where

[TABLE]

Following this natural—but naive—formula, we thus need approximations of the correctors $\chi^{1}$ to $\chi^{2s+1}$. In this section, we present a new relation between the correctors ensuring that in fact the correctors $\chi^{1}$ to $\chi^{s+1}$ are sufficient (Theorem 3.5).

Let us quantify the computational gain achieved thanks to this result. As each $\chi^{r}$ is a symmetric tensor valued function, it has $\binom{r+d-1}{r}$ distinct components (see Remark 2.2). The number of cell problems to solve to have all the distinct components of $\chi^{1}$ to $\chi^{k}$ is thus

[TABLE]

Hence, to compute $\{\check{q}^{r}\}_{r=1}^{s}$ Theorem 3.5 allows to avoid the approximation of $\chi^{s+2},\ldots,\chi^{2s+1}$, i.e., we spare the approximation of

[TABLE]

cell problems. To fully appreciate this gain, assume that we want to compute the effective tensors of an effective equation for a timescale $\mathcal{O}(\varepsilon^{-6})$, i.e. $s=\lfloor\alpha/2\rfloor=3$: for $d=2$, we spare $21$ cell problems ($14$ cell problems to solve instead of $35$) and if $d=3$, we spare $85$ cell problems ($34$ cell problems to solve instead of $119$).

Remark 3.4*.*

It is important to note that this gain is not specific to the effective equations presented here. Indeed, formula (3.8) can directly be used to compute the (sufficient) symmetric part of the high order tensors in the effective models from [9] and [13] discussed in Section 2.3.

The new relation between the correctors is presented in the following result, proved in Section 5.4. Note that this result was obtained independently in [22, Remark 2.5] in the context of stochastic homogenization.

Theorem 3.5**.**

Let $\{\chi^{k}_{i_{1}\cdot\cdot i_{k}}\}_{k=1}^{2r+1}$ be the zero mean correctors defined in (2.10). Then, for $1\leq r\leq\lfloor\alpha/2\rfloor$, the tensor $g^{2r}$ defined by (3.7) satisfies the decomposition

[TABLE]

where the tensors $k^{r},h^{r}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$ are defined as

[TABLE]

where the second double sum in the definition of $h^{r}$ vanishes for $r=1$.

Observe that the tensor $h^{r}$ only depends on the correctors $\chi^{1}$ to $\chi^{r}$ so that decomposition (3.8) indeed guarantees that $S^{2r+2}(g^{2r})$ can be computed from $\chi^{1}$ to $\chi^{r+1}$.

Remark 3.6*.*

It can be verified that the homogenized tensor (2.13) satisfies $a^{0}=g^{0}=k^{0}$. Hence, decomposition (3.8) also holds in the case $r=0$ with $h^{0}=0$ (recall that $\chi^{0}=1$).

3.3 Matrix associated to a symmetric tensor of even order

In this section, we construct the “matricization” operator used in the construction of Section 3.1. This operator maps a given symmetric tensor of even order $q$ to a matrix whose sign is the same as $q$.

We consider the bilinear map

[TABLE]

Denote $I(d,n)$ the set of multiindices of the distinct entries of a tensor in ${\mathrm{Sym}}^{n}(\mathbb{R}^{d})$, i.e.,

[TABLE]

We verify that the cardinality of $I(d,n)$ is $N(d,n)=|I(d,n)|=\binom{d+n-1}{n}.$ We denote $J(d,n)=\{1,\ldots,N(d,n)\}$ and let $\ell:J(d,n)\to I(d,n)$ be a bijection. We define then the bijective mapping

[TABLE]

For $i\in I(d,n)$, let $z(i)$ be the number of multiindices in $\{1,\ldots,d\}^{n}$ that are equal to $i$ up to symmetries, i.e.,

[TABLE]

With these notations, we rewrite the map defined in (3.10) as

[TABLE]

Finally, we define the matrix associated to a tensor as

[TABLE]

We verify that for any $\xi,\eta\in{\mathrm{Sym}}^{n}(\mathbb{R}^{d})$, $q\xi:\eta=M(q)\nu(\xi)\cdot\nu(\eta)$. Hence, $q$ is positive definite (resp. semidefinite) if and only if $M(q)$ is positive definite (resp. semidefinite). In particular, property (3.3) holds.

3.4 Algorithm for the computation of the high order effective tensors

We present here the algorithm for the computation of the effective tensors of one equation in the family $\mathcal{E}$. The procedure relies on the construction explained in Lemma 3.2, the formula provided by Theorem 3.5, and the “matricization” operator $M$ defined in (3.12). The set of index $I(d,n)$ is defined in (3.11).

Algorithm 7 Compute the tensors of an effective equation (2.1) in the family $\mathcal{E}$.

1:$Y$-periodic, $d\times d$ symmetric tensor $a(y)$ satisfying (1.11); timescale $\alpha$.

2:effective tensors $a^{0},\{a^{2r},b^{2r}\}_{r=1}^{\lfloor\alpha/2\rfloor}$.

3:for all $i\in I(d,1)$ solve the cell problem for $\chi^{1}_{i}$ in (2.10a) with \big{\langle}{\chi^{1}_{i}}\big{\rangle}_{Y}=0

4:for all $i\in I(d,2)$ c^{0}_{i_{1}i_{2}}=a^{0}_{i_{1}i_{2}}=-\big{\langle}{a\nabla_{\makebox[3.6831pt]{\scriptstyle y}}\kern-1.5pt\chi_{i_{2}}^{1}\cdot\nabla_{\makebox[3.6831pt]{\scriptstyle y}}\kern-1.5pt\chi_{i_{1}}^{1}}\big{\rangle}_{Y}+\big{\langle}{ae_{i_{2}}\cdot e_{i_{1}}}\big{\rangle}_{Y}

5:for $r=1,\ldots,\lfloor\alpha/2\rfloor$ do

6:      for all $i\in I(d,r+1)$ solve the cell problems for $\chi^{r+1}_{i_{1}\cdot\cdot i_{r+1}}$ in (2.10b) with \big{\langle}{\chi^{r+1}_{i_{1}\cdot\cdot i_{r+1}}}\big{\rangle}_{Y}=0

7:      compute the tensor S^{2r+2}\big{(}g^{2r}\big{)} with formula (3.8)

8:      compute the tensor $\check{q}^{r}$ with formula (3.2)

9:      build the matrices $M(\check{q}^{r})$, A^{r}=M\big{(}S^{2r+2}(\otimes^{r+1}a^{0})\big{)} and compute \delta^{*}=\bigg{\{}-\displaystyle\frac{\lambda_{\min}(M(\check{q}^{r}))}{\lambda_{\min}(A^{r})}\bigg{\}}_{+}

10:      compute a^{2r}=\check{q}^{2r}+\delta^{*}S^{2r+2}\big{(}\otimes^{r+1}a^{0}\big{)},  b^{2r}=\delta^{*}S^{2r}\big{(}\otimes^{r}a^{0}\big{)}, $c^{r}=a^{2r}-\displaystyle\sum_{\ell=0}^{r-1}c^{\ell}\!\otimes\!b^{2(r-\ell)}$

11:end for

4 Numerical experiments

In this section, we present numerical experiments to illustrate the result of Theorem 2.8: the effective equations in the family capture the long-time behavior of $u^{\varepsilon}$. Note that reporting numerical error in the approximation of $u^{\varepsilon}$ in a pseudo-infinite medium for very large timescales is not conceivable as computing a reference solution is out of reach even for one-dimensional problems. We can however consider a small periodic domain as this setting is covered by our theory. In addition, we will illustrate that high order effective models are also useful when we deal with high frequency regimes.

We consider the one-dimensional model problem (1.13) given by the data $u_{0}(x)=e^{-4{x^{2}}}$, $u_{1}=f=0$, $a(y)=\sqrt{2}-\cos(2\pi y)$ (we verify that $a^{0}=1$) with $\varepsilon=1/10$, and the periodic domain $\Omega=(-L,L)$, $L=84$. Let us give some insight on the macroscopic evolution of $u^{\varepsilon}$: the central pulse $u_{0}$ separates into left- and right-going waves packets with speed $a^{0}=\pm 1$; these packets meet at $x=L$ (equivalently $x=-L$) for $t=L+2kL$, $k\in\mathbb{N}$, and at $x=0$ for $t=2kL$, $k\in\mathbb{N}$. As time increases, dispersion appears in each packet. As the domain has a periodic boundary, at some point in time the packets start to superpose. In our settings, $t=\varepsilon^{-4}$ is the larger timescale until the bulk of the dispersion spans approximately half of the domain $\Omega$. For $\alpha=0,2,4$, we compare $u^{\varepsilon}$ with the corresponding effective solutions of order $s=\lfloor{\alpha/2}\rfloor$, denoted $\tilde{u}^{\lfloor{\alpha/2}\rfloor}$ ($\tilde{u}^{0}$ is the homogenized solution). We approximate $u^{\varepsilon}$ using a spectral method grid of size $h=\varepsilon/16$ and a leap frog scheme for the time integration with $\Delta t=h/16$. The effective solutions $\{\tilde{u}^{s}\}_{s=0}^{2}$ are approximated with a Fourier method (the effective coefficients are computed using Algorithm 3.4) on a grid of size $h=\varepsilon/16$ (no time integration). Details on the numerical methods can be found in [35, section 5.3].

The results are displayed in Figure 1. In the top-left plot, we compare $u^{\varepsilon}$, $\tilde{u}^{0}$, $\tilde{u}^{1}$ at $t=\varepsilon^{-2}=10^{2}$. We observe that $\tilde{u}^{0}$ does not capture the macroscopic dispersion developed by $u^{\varepsilon}$, while $\tilde{u}^{1}$ accurately describes it. In the top-right plot, we compare $u^{\varepsilon}$, $\{\tilde{u}^{s}\}_{s=0}^{2}$, at $t=\varepsilon^{-4}=10^{4}$ (observe that right- and left- going wave packets are in fact superposed). At this timescale, the first order effective solution $\tilde{u}^{1}$ does not capture accurately the dispersion anymore, as expected, but the second order $\tilde{u}^{2}$ does. In the bottom plot, we compare the normalized errors $\mathrm{err}(\tilde{u}^{s})(t)=\|(u^{\varepsilon}-\tilde{u}^{s})(t)\|_{{\mathrm{L}^{2}}(\Omega)}/\|u^{\varepsilon}(t)\|_{{\mathrm{L}^{2}}(\Omega)}$ for the effective solutions $\{\tilde{u}^{s}\}_{s=0}^{2}$ on the time interval $[0,\varepsilon^{-4}]$ (the $x$-axis is in log scale). We observe that $\tilde{u}^{0}$ is accurate up to $\varepsilon^{-1}$ and then deteriorates. The first order $\tilde{u}^{1}$ is accurate up to $\varepsilon^{-3}$, then deteriorates. Finally, the second order $\tilde{u}^{2}$ has a satisfying accuracy over the whole time interval $(0,\varepsilon^{-4})$. The observations of this experiment corroborate the result of Theorem 2.8: the effective solution $\tilde{u}^{\lfloor{\alpha/2}\rfloor}$ accurately describe $u^{\varepsilon}$ up to $\mathcal{O}(\varepsilon^{-\alpha})$ timescales.

In two dimensions, computing a reference solution for the previous experiment has a huge computational cost. However, we can illustrate an interesting fact: high order effective models are useful in high frequency regimes. Let us first give some insights on this fact. We consider the following specific settings: let a\big{(}\tfrac{x}{\varepsilon}\big{)} be a given tensor with $\varepsilon>0$ fixed, $g(x)=e^{\beta^{2}|x|^{2}}$ be a Gaussian with $\beta=\mathcal{O}(1)$ and $\nu>0$ a scaling parameter. In (1.13), we let $f=0$ and the initial conditions $u_{0}(x)=g(\nu x)$ and $u_{1}=0$ (we assume that $\Omega$ is arbitrarily large). To specify the dependence of $u^{\varepsilon}$ on the parameter $\nu$, we denote it as $u^{\varepsilon}_{\nu}$. In the equation for $u^{\varepsilon}_{\nu}$, making the changes of variables $\hat{x}=\nu x$ and $\hat{t}=\nu t$ and introducing $\hat{u}^{\varepsilon}_{\nu}(\hat{t},\hat{x})=u^{\varepsilon}_{\nu}(\hat{t}/\nu,\hat{x}/\nu)$, we find

[TABLE]

and we conclude that $u^{\varepsilon}_{\nu}(\hat{t}/\nu,\hat{x}/\nu)=\hat{u}^{\varepsilon}_{\nu}(\hat{t},\hat{x})=u^{\nu\varepsilon}_{1}(\hat{t},\hat{x})$. In other words, for $\nu>1$ (i.e., an increase in the frequencies of the initial wave) the long-time effects of $u^{\nu\varepsilon}_{1}$ can be observed at a shorter time in $u^{\varepsilon}_{\nu}$ (modulo a contraction of space). However, for high values of $\nu$ we meet situations where Theorem 2.8 does not provide a satisfactory error estimate: on the one hand, in the estimate for $u^{\varepsilon}_{\nu}$, the increase of $\nu$ deteriorates the error constant in Theorem 2.8; one the other hand, in the estimate for $u^{\nu\varepsilon}_{1}$, the constant is good, but the period of the tensor $\varepsilon^{\prime}=\nu\varepsilon$ is to far from the asymptotic regime for homogenization to be meaningful. In practice, we observe that the increase of the frequencies of the initial wave (increase of $\nu$) leads to additional dispersive effects in $u^{\varepsilon}_{\nu}$ (or $u^{\nu\varepsilon}_{1}$). Furthermore, the use of higher order effective models allows to capture these additional effects.

To illustrate this, we consider the two dimensional model problem given by the data $f=u_{1}=0$ and

[TABLE]

We compute the effective tensors using Algorithm 3.4 and approximate $u^{\varepsilon}_{\nu}$ and $\{\tilde{u}^{s}_{\nu}\}_{s=1}^{3}$ at time $t=20$. As above, $u^{\varepsilon}$ is approximated with a spectral method on a grid of size $h=\varepsilon/16$ and leap frog scheme with $\Delta t=h/100$ and $\{\tilde{u}^{s}\}_{s=0}^{2}$ are computed with a Fourier method on a grid of size $h=\varepsilon/16$ (see [35, section 5.3]). In Figure 2, we compare the solution in the periodic medium $u^{\varepsilon}_{\nu}$ (top-left) with the effective solutions of order $1$, $\tilde{u}^{1}_{\nu}$ (top-right), and order $2$, $\tilde{u}^{2}_{\nu}$ (bottom-left). We observe that $\tilde{u}^{2}_{\nu}$ captures more accurately the dispersion developed by $u^{\varepsilon}_{\nu}$ than $\tilde{u}^{1}_{\nu}$. This is even better seen in the 1d cut at $\{x_{1}=0\}$ in the bottom-right plot of Figure 2. Furthermore, even though distinguishing the higher order dispersion from the $\varepsilon$-scale oscillation is not easy in this regime, in the zoom we can guess that the model of order $3$ is better than the order $2$.

5 Proofs of the main results

In this section, we provide the proofs of the main results of the paper. In 5.1, we prove the inductive Boussinesq tricks: Theorem 5.1 and Lemma 5.7. In 5.2, the a priori error estimate for the family of effective equations of Theorem 2.8 is proved. In 5.3, we prove the result on positive definite tensors from Lemma 3.1 and in 5.4, we prove the new relation between the correctors from Theorem 3.5.

5.1 Inductive Boussinesq tricks for the derivation of the family of effective equations

In this section, we present the technical task that was postponed in the derivation of the family of effective equations in Section 2.1. Specifically, we provide the result allowing to substitute the time derivatives in the terms of order $\mathcal{O}(\varepsilon^{0})$ to $\mathcal{O}(\varepsilon^{\alpha})$ in $r^{\varepsilon}$ (2.5). To that end, the main challenge is to proceed to inductive Boussinesq tricks (Theorem 5.1). This result is then used in Lemma 5.7, which provides the specific relation used in Section 2.1.

5.1.1 Inductive Boussinesq tricks

The following theorem is the key result of this section.

Theorem 5.1**.**

Let $\tilde{u}$ be the solution of (2.1) and let $N$ be an even integer such that $0\leq N\leq 2\lfloor\alpha/2\rfloor$. If the right-hand side of (2.1) is defined as

[TABLE]

then $\tilde{u}$ satisfies

[TABLE]

where $c^{r}$ is the tensor defined inductively as

[TABLE]

and the remainders $\mathcal{R}^{\varepsilon}_{N}\tilde{u}$ and $\mathcal{S}^{\varepsilon}_{N}f$ are defined in (5.22) and satisfy the following estimates: for any integer $n$ such that $N+n\leq\alpha$

[TABLE]

where $r(\alpha,N,n)=2\lfloor\alpha/2\rfloor+N+n$ and the constants depend only on the tensors $a^{0},\{a^{2r},b^{2r}\}_{r=1}^{\lfloor\alpha/2\rfloor}$.

Remark 5.2*.*

We verify that the definition (5.1) of $Qf$ ensures the remainder $\mathcal{S}^{\varepsilon}_{N}f$ to have maximal order in terms of $\varepsilon$ in the second estimate in (5.4). More specifically, among the right-hand side of the form $Qf=\sum_{s=0}^{\lfloor\alpha/2\rfloor}(-1)^{s}\varepsilon^{2s}q^{s}\nabla^{2s}_{x}f$, $q^{s}\in{\mathrm{Ten}}^{2s}(\mathbb{R}^{d})$, (5.1) is the only one that ensures $\|\mathcal{S}^{\varepsilon}_{N}f\|=\mathcal{O}(\varepsilon^{N+2})$ for all even $N$ such that $0\leq N\leq 2\lfloor\alpha/2\rfloor$. This affirmation is proved in Remark 5.6.

The result of Theorem 5.1 relies on inductive Boussinesq tricks. Let us summarize this process. We start from the expression of ${\partial^{2}_{t}}\tilde{u}$ given by the effective equation (2.1), which we rewrite as

[TABLE]

where $Qf$ is defined in (5.1) (see Remark 5.2). The result for $N=0$ trivially follows (5.5) with $\mathcal{S}^{\varepsilon}_{0}f=G^{0}(0)$ and $\mathcal{R}^{\varepsilon}_{0}\tilde{u}=T^{0}(0)$ (defined in (5.10) below). Let us then assume that $N$ is an even number such that $2\leq N\leq\alpha$. In the right-hand side of (5.5), we inductively substitute ${\partial^{2}_{t}}\tilde{u}$ (using (5.5) itself) in the terms of order $\mathcal{O}(\varepsilon^{2})$ to $\mathcal{O}(\varepsilon^{N})$. To understand this technical process at best, let us define inductively some quantities and functions. Note that these definitions arise in the result of Lemma 5.3, which exhibits the decomposition process of one Boussinesq trick.

We define the tensors

[TABLE]

and

[TABLE]

Associated with the tensors $A^{r}(j)$ and $B^{r}(j)$, we define for the given $N$ the functions

[TABLE]

With these notations, we verify that (5.5) reads (recall the definition of $Qf$ in (5.1))

[TABLE]

where the remainders are $T^{N}(0)=G^{N}(0)=0$ if $N=\alpha$ and

[TABLE]

otherwise.

The proof is divided into three steps. In the first step (Lemma 5.3), we apply one Boussinesq trick: in $S^{N}(j)$, we use (5.5) to substitute ${\partial^{2}_{t}}\tilde{u}$ and decompose the result into terms of interest plus remainders. In the second step (Lemma 5.4), the decomposition is used inductively to obtain a first version of Theorem 5.1, which involves the tensor $\tilde{c}^{r}$ instead of $c^{r}$. The final step is to prove that $\tilde{c}^{r}$ and $c^{r}$ are equal (Lemma 5.5).

Lemma 5.3**.**

The functions defined in (5.8) satisfy the relations ((5.11a) and (5.11b) hold only if $N\geq 4$)

[TABLE]

where the remainders $\{T^{N}(j)\}_{j=1}^{N/2}$ and $\{G^{N}(j)\}_{j=1}^{N/2}$ are defined in (5.16), (5.17), (5.19), and (5.20).

Lemma 5.4**.**

The effective solution of (2.1) $\tilde{u}$ satisfies the equality

[TABLE]

where $\tilde{c}^{r}$ is the tensor defined as

[TABLE]

and the remainders $\mathcal{R}^{\varepsilon}_{N}\tilde{u}$ and $\mathcal{S}^{\varepsilon}_{N}f$ are defined in (5.22) and satisfy the estimates (5.4).

Lemma 5.5**.**

For $1\leq r\leq\lfloor\alpha/2\rfloor$, the tensor $\tilde{c}^{r}$, defined in (5.13), equals the tensor $c^{r}$, defined in (5.3).

Remark 5.6*.*

Let us explain how the definition of $Qf$ is obtained. Assuming that $Qf$ is unknown, (5.9) reads $\partial_{t}\tilde{u}=Qf+R^{N}(0)+S^{N}(0)+T^{N}(0)$. Using inductively (5.11a) and (5.11c), which hold independently of $Qf$, we obtain (5.12) with $\mathcal{S}^{\varepsilon}_{N}f$ defined as

[TABLE]

We want to build $Qf$ so that $\mathcal{S}^{\varepsilon}_{N}f$ has order $\mathcal{O}(\varepsilon^{N+2})$ for all even $0\leq N\leq 2\lfloor\alpha/2\rfloor$. Inserting the ansatz $Qf=\sum_{s=0}^{\lfloor\alpha/2\rfloor}(-1)^{s}\varepsilon^{2s}q^{s}\nabla^{2s}_{x}f$, where $q^{s}\in{\mathrm{Ten}}^{2s}(\mathbb{R}^{d})$ are unknown, we obtain (after some work)

[TABLE]

Canceling the terms of order $\mathcal{O}(1)$ to $\mathcal{O}(\varepsilon^{N})$, we obtain $q^{0}=1$ and an inductive definition for $q^{m}$, $m\geq 1$. We then verify by induction that $q^{m}=b^{2m}$ for $1\leq m\leq\lfloor\alpha/2\rfloor$, i.e., $Qf$ is defined as (5.1).

*Proof of Lemma 5.3. * To slightly simplify the notation, Let us denote $\bar{\alpha}=2\lfloor\alpha/2\rfloor$. We first prove (5.11a). Using (5.5), we substitute ${\partial^{2}_{t}}\tilde{u}$ in $S^{N}(j)$ and obtain

[TABLE]

where the remainder is $T^{N}_{1}(j+1)=0$ if $N=\bar{\alpha}$ and

[TABLE]

otherwise. In order to decompose the terms in (5.14), we study the two index sets (denoting $\ell=N/2$)

[TABLE]

where we recall that $0\leq j\leq\ell-2$. We verify that these sets can be written as

[TABLE]

Hence, we rewrite (5.14) as

[TABLE]

where the remainder is

[TABLE]

Using the definitions of $A^{r}(j)$ in (5.7) and $B^{r}(j)$ in (5.6), we verify that the two sums in (5.15) match $R^{N}(j+1)$ and $S^{N}(j+1)$, respectively. That proves (5.11a).

Next, let us prove (5.11b). Using the definition of $Qf$ in (5.1) and recalling (5.8c) and (5.8d), we have

[TABLE]

As done above, we decompose the sum and use the expression of $J^{\ell}_{m}(j)$ to obtain

[TABLE]

where the remainder is

[TABLE]

Using the definitions of $B^{m}(j+1)$ in (5.6) and $F^{N}(j+1)$ in (5.8c), we obtain (5.11b).

Next, we prove (5.11c). From the definition of $S^{N}(j)$ in (5.8b), we use (5.5) to get

[TABLE]

where the remainder is

[TABLE]

The first term of the right-hand side of (5.18) matches the definition of $R^{N}(\bar{\alpha}/2)$ in (5.8a). That proves (5.11c).

Finally, let us prove (5.11d). From the definition of $\tilde{F}^{N}(j)$, we verify that (5.11d) holds with $G^{N}(N/2)$ defined as

[TABLE]

That concludes the proof of Lemma 5.3. $\square$

*Proof of Lemma 5.4. * Starting from (5.9), we use inductively the decompositions in (5.11) and obtain

[TABLE]

We define the remainders as

[TABLE]

where $\{T^{N}(j)\}_{j=0}^{N/2}$ and $\{G^{N}(j)\}_{j=0}^{N/2}$ are defined in (5.10), (5.16), (5.17), (5.19), and (5.20). From the definitions of $T^{N}(j)$ and $G^{N}(j)$, we verify that $\mathcal{R}^{\varepsilon}_{N}\tilde{u}$ and $\mathcal{S}^{\varepsilon}_{N}f$ satisfy the estimates in (5.4). Next, we have to develop $\sum_{j=0}^{N/2}R^{N}(j)$ in (5.21). Using the definition of $R^{N}(j)$ in (5.8a) and exchanging the sums, we find that

[TABLE]

Using the definition of $A^{r}(j)$ in (5.7), we verify that for $r=0$, $\sum_{j=0}^{r}A^{r}(j)=a^{0}=\tilde{c}^{0}$, and for $1\leq r\leq\alpha/2$,

[TABLE]

Combining (5.21), (5.22), (5.23), and (5.24), we obtain (5.12) and the proof of Lemma 5.4 is complete. $\square$

*Proof of Lemma 5.5. * We prove by induction on $r$ that the tensor $\tilde{c}^{r}$, defined in (5.13), equals the tensor $c^{r}$, defined in (5.3). The base case is trivially verified as $\tilde{c}^{0}=a^{0}=c^{0}$. Then, for $r\geq 2$, we assume that $\tilde{c}^{s}=c^{s}$ for $1\leq s\leq r-1$ and we have to verify that the tensor

[TABLE]

equals $\tilde{c}^{r}$, defined in (5.13). Using the induction assumption and (5.13), we have

[TABLE]

Let us denote the triple sum $T$ and its summand $x^{r}_{\ell,j,s}=a^{2s}\!\otimes\!B^{\ell-s}(j-1)\!\otimes\!b^{2(r-\ell)}$. Applying the change of indices $m=r-\ell$ and exchange the sums twice to get

[TABLE]

Using the definition of $B^{r}(j)$ in (5.6), we then have

[TABLE]

We change the index $k=j+1$ in this equality to obtain from (5.25) (recall that $B^{r}(0)=-b^{2r}$)

[TABLE]

This expression matches the definition of $\tilde{c}^{r}$ in (5.13). Hence, we have proved that $\tilde{c}^{r}=c^{r}$ for all $0\leq r\leq\lfloor\alpha/2\rfloor$ and the proof of Lemma 5.5 is complete. $\square$

5.1.2 Use of the inductive Boussinesq tricks

In the asymptotic expansion in Section 2.1, we need to substitute the terms involving ${\partial^{2}_{t}}\tilde{u}$ in the development of $r^{\varepsilon}$ in (2.5). This task is performed in the following result, obtained thanks to Theorem 5.1.

Lemma 5.7**.**

Define tensors $p^{k}\in{\mathrm{Ten}}^{k}(\mathbb{R}^{d})$ as

[TABLE]

where $c^{r}$ is defined in (5.3) and let $\chi^{0}=1$ and $\{\chi^{k}\}_{k=1}^{\alpha}$ be the tensor functions in (2.7). Then the following equality holds

[TABLE]

where the remainders $\mathcal{S}^{\varepsilon}f$ and $\mathcal{R}^{\varepsilon}\tilde{u}$ are defined in (5.29), (5.30). In particular, provided sufficient regularity of $\tilde{u}$, $\mathcal{R}^{\varepsilon}\tilde{u}$ has order $\mathcal{O}(\varepsilon^{\alpha+1})$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathrm{L}^{2}}(\Omega)})$ norm.

*Proof of Lemma 5.7. * Let us define $T\coloneqq\sum_{k=0}^{\alpha}\varepsilon^{k}\chi^{k}\nabla^{k}_{x}{\partial^{2}_{t}}\tilde{u}$. Using Theorem 5.1 and the definition of $p^{k}$, we obtain for any choice of even integers $0\leq N(k)\leq\alpha$:

[TABLE]

where $z^{j}_{k}$ and the remainders are defined as

[TABLE]

with $\mathcal{S}^{\varepsilon}_{N(k)}f$ and $\mathcal{R}^{\varepsilon}_{N(k)}\tilde{u}$ the remainders defined in (5.22). The first estimate in (5.4) implies that

[TABLE]

Hence, as our goal is for this quantity to have order $\mathcal{O}(\varepsilon^{\alpha+1})$ (see the discussion in Section 2.1), we need to set $N(k)$ such that

[TABLE]

Dealing with even and odd indices separately, we verify that the smallest even integer satisfying the above inequality is

[TABLE]

We now rewrite the double sum in the right-hand side of (5.28) separately for even and odd index $k$. For the sum on even index \big{\{}k=2\ell\,:\,0\leq\ell\leq\lfloor\alpha/2\rfloor\big{\}}, we have

[TABLE]

where in the first equality we changed the index $m=r+\ell$ and in the second we exchanged the order of summation and used that $z^{2s+1}_{k}=0$ for any $s,k$. Similarly, for the sum on odd index \big{\{}k=2\ell-1\,:\,1\leq\ell\leq\lceil\alpha/2\rceil\big{\}}, we have

[TABLE]

Using the two last equalities in (5.28), we gather the sums to find

[TABLE]

Replacing $z^{k-j}_{j}$ by its definition in (5.29), we obtain (5.27) and that concludes the proof of Lemma 5.7. $\square$

5.2 Proof of the a priori error estimate for the family of effective equations (Theorem 2.8)

In this section, we prove Theorem 2.8. Note that the main ingredient of the proof is an adaptation operator based on the adaptation constructed in Section 2.1.

Let us first introduce some notations. We use a bracket \scalebox{2.1}[1.0]{[}{v}\scalebox{2.1}[1.0]{]} to denote the equivalence class of $v\in{\mathrm{L}^{2}}(\Omega)$ in ${\mathcal{L}^{2}}(\Omega)$, and a bold face letter $\boldsymbol{v}$ to denote elements of ${\mathcal{W}_{{\mathrm{per}}}}(\Omega)$. Note that the Hilbert space ${\mathcal{L}^{2}}(\Omega)$ is equipped with the inner product \big{(}\scalebox{2.1}[1.0]{[}{v}\scalebox{2.1}[1.0]{]},\scalebox{2.1}[1.0]{[}{w}\scalebox{2.1}[1.0]{]}\big{)}_{{\mathcal{L}^{2}}(\Omega)}=\big{(}v-\langle{v}\rangle_{\Omega},w-\langle{w}\rangle_{\Omega}\big{)}_{{\mathrm{L}^{2}}(\Omega)}. We define the following norm on ${\mathcal{W}_{{\mathrm{per}}}}(\Omega)$

[TABLE]

In particular, we verify that a function $w\in{{\mathrm{W}_{\!{\mathrm{per}}}}(\Omega)}$ satisfies \|w\|_{W}=\|\scalebox{2.1}[1.0]{[}{w}\scalebox{2.1}[1.0]{]}\|_{\mathcal{W}}, where $\|{\cdot}\|_{W}$ is defined in (1.5).

The proof of Theorem 2.8 is structured as follows. First, based on the adaptation $\mathcal{B}^{\varepsilon}\tilde{u}$ defined by (2.3) and (2.7), we define the adaptation operator $\boldsymbol{\mathcal{B}}^{\varepsilon}\tilde{u}$. Then, using the triangle inequality, we split the error as

[TABLE]

and estimate the two terms of the right-hand side separately.

Let us first discuss the consequences of the assumptions made in the theorem. The fact that the effective equation belongs to the family $\mathcal{E}$ (Definition 2.6) has two major implications: first, the equation is well-posed; second, the tensors $\{a^{2r},b^{2r}\}_{r=1}^{\lfloor\alpha/2\rfloor}$ satisfy the constraints (2.15) and thus the cell problems (2.10) are well-posed. Hence, we have the existence and uniqueness of the effective solution $\tilde{u}$ and of the correctors $\chi^{1},\ldots,\chi^{\alpha+2}$. Inductively, we can show that the correctors satisfy the following bound

[TABLE]

where the constant $C$ depends only on the ellipticity constant $\lambda$ and the reference cell $Y$. Let us then investigate the regularity of $\tilde{u}$ and $f$. As we assume $d\leq 3$, the embedding $\mathrm{H}^{2}_{\mathrm{per}}(\Omega)\hookrightarrow\mathcal{C}^{0}_{\mathrm{per}}(\bar{\Omega})$ is continuous and the regularity assumption ensures that

[TABLE]

where we recall that $r(\alpha)=\alpha+2\lfloor\alpha/2\rfloor+2$ (as $\alpha\geq 2$ we have $r(\alpha)\geq\alpha+4$).

We are now able to define the adaptation operator. Thanks to the regularity of $f$ and the correctors, $\mathcal{S}^{\varepsilon}f$ defined in (5.29) belongs to ${\mathrm{L}^{2}}(0,\varepsilon^{-\alpha}T;{{\mathcal{L}^{2}}(\Omega)})$. We can thus define $\boldsymbol{\varphi}$ as the unique solution in ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathcal{W}_{{\mathrm{per}}}}(\Omega)})$ of

[TABLE]

where $\mathcal{A}^{\varepsilon}$ is defined in (2.4). The adaptation operator is then defined as

[TABLE]

where $\mathcal{B}^{\varepsilon}\tilde{u}$ is the adaptation constructed in Section 2.1 (see (2.3) and (2.7)):

[TABLE]

Assumption (1.12) ensures that x\mapsto\chi^{k}\big{(}\tfrac{x}{\varepsilon}\big{)}\nabla^{k}_{x}\tilde{u}(t,x) is $\Omega$-periodic. Using the regularity of the correctors and $\tilde{u}$, we verify that the adaptation

[TABLE]

(notice that the regularity of $\boldsymbol{\varphi}$ and its derivatives are weaker than that of $\mathcal{B}^{\varepsilon}\tilde{u}$).

Define then \boldsymbol{\eta}^{\varepsilon}(t)=\boldsymbol{\mathcal{B}}^{\varepsilon}\tilde{u}(t)-\scalebox{2.1}[1.0]{[}{u^{\varepsilon}(t)}\scalebox{2.1}[1.0]{]} and recall that we developed $r^{\varepsilon}=({\partial^{2}_{t}}+\mathcal{A}^{\varepsilon})(\mathcal{B}^{\varepsilon}\tilde{u}-u^{\varepsilon})$ in Section 2.1. The development in (2.9) together with the cell problems (2.10) ensure that

[TABLE]

where $\mathcal{R}_{\mathrm{ini}}^{\varepsilon}\tilde{u}(t)$ is defined in (2.6) and $\mathcal{R}^{\varepsilon}\tilde{u}$ in (5.29). The following error estimate can then be used to quantify $\boldsymbol{\eta}^{\varepsilon}$ (see [5, Corollary 2.2] or [35, Corollary 4.2.2]).

Lemma 5.8**.**

Assume that $\boldsymbol{\eta}\in{\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathcal{W}_{{\mathrm{per}}}}(\Omega)})$, with $\partial_{t}\boldsymbol{\eta}\in{\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathcal{L}^{2}}(\Omega)})$, ${\partial^{2}_{t}}\boldsymbol{\eta}\in{\mathrm{L}^{2}}(0,\varepsilon^{-\alpha}T;{{\mathcal{W}_{{\mathrm{per}}}^{*}}(\Omega)})$ satisfies

[TABLE]

where $\boldsymbol{\eta}^{0}\in{{\mathcal{W}_{{\mathrm{per}}}}(\Omega)}$, $\boldsymbol{\eta}^{1}\in{{\mathcal{L}^{2}}(\Omega)}$, and $\boldsymbol{r}\in{\mathrm{L}^{1}}(0,\varepsilon^{-\alpha}T;{{\mathcal{L}^{2}}(\Omega)})$. Then the following estimate holds

[TABLE]

where $C(\lambda)$ depends only on the ellipticity constant $\lambda$. If in addition, $\boldsymbol{r}\in{\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;{{\mathcal{L}^{2}}(\Omega)})$, then

[TABLE]

In order to estimate the remainder terms in (5.32), we also need the following result.

Lemma 5.9**.**

Let $\gamma\in\mathrm{L}^{2}_{\mathrm{per}}(Y)$ and $v\in\mathrm{H}^{2}_{\mathrm{per}}(\Omega)$. Then the following estimate holds

[TABLE]

where the constant $C$ depends only on $Y$ and $d$.

Proof.

Recall that $Y=(0,\ell_{1})\times\cdots\times(0,\ell_{d})$ and $\Omega=(\omega_{1}^{l},\omega_{1}^{r})\times\cdots\times(\omega^{l}_{d},\omega^{r}_{d})$. As $\Omega$ satisfies (1.12), the numbers $N_{i}=\frac{\omega_{i}^{l}-\omega_{i}^{r}}{\ell_{i}\varepsilon}$ are integers and the cells constituting $\Omega$ belongs to the set $\{\varepsilon(n\cdot\ell+Y):0\leq n_{i}\leq N_{i}-1\}$. Denoting $\Xi=\{\xi=n\cdot\ell:0\leq n_{i}\leq N_{i}-1\}$, the domain $\Omega$ satisfies

[TABLE]

Hence, almost every $x\in\Omega$ can be written as $x=\varepsilon(\xi+y)$ for some $\xi\in\Xi,y\in Y$. For such triplet $(x,\xi,y)$, the $Y$-periodic function $\gamma$ satisfies \gamma\big{(}\tfrac{x}{\varepsilon}\big{)}=\gamma(\xi+y)=\gamma(y). Let $Z\subset\mathbb{R}^{d}$ be an open set with a $\mathcal{C}^{1}$ boundary, that contains $Y$ and is contained in its neighborhood, i.e.,

[TABLE]

As $Z$ has a $\mathcal{C}^{1}$ boundary and $d\leq 3$, Sobolev embedding theorem ensures that the embedding $\mathrm{H}^{2}(Z)\hookrightarrow\mathcal{C}^{0}(\bar{Z})$ is continuous. Hence, there exists a constant $C_{Y}$, depending only on $Y$, such that

[TABLE]

We now prove (5.33). Using (5.34), we have

[TABLE]

where we made the change of variables $x=\varepsilon(\xi+y)$. As $v\in\mathrm{H}^{2}_{\mathrm{per}}(\Omega)\hookrightarrow\mathcal{C}^{0}_{\mathrm{per}}(\bar{\Omega})$, we have

[TABLE]

where $v_{\xi,\varepsilon}$ is the function of $\mathcal{C}^{0}(\bar{Y})$ defined by v_{\xi,\varepsilon}(y)=v\big{(}\varepsilon(\xi+y)\big{)}. Using (5.35) gives $\|v_{\xi,\varepsilon}\|_{\mathcal{C}^{0}(\bar{Y})}\leq C_{Y}\|v_{\xi,\varepsilon}\|_{\mathrm{H}^{2}(N_{Y})}$. Furthermore, we have

[TABLE]

As $\partial_{y_{i}}v_{\xi,\varepsilon}=\varepsilon\partial_{x_{i}}v$ and $\partial^{2}_{y_{ij}}v_{\xi,\varepsilon}=\varepsilon^{2}\partial^{2}_{x_{ij}}v$, the change of variable $x=\varepsilon(\xi+y)$ leads to

[TABLE]

where we used that every cell $\varepsilon(\xi+Y)$ belongs to the neighborhoods of $(6d-3)$ cells (including itself). This proves (5.33) and the proof of the lemma is complete. ∎

With Lemma 5.9, we can estimate the remainders in (5.32). Using (5.29) and estimate (5.4), we obtain

[TABLE]

From (5.30), we verify that $\alpha-1\leq N(k)+k\leq\alpha$ and thus

[TABLE]

Similarly, using (5.4) to estimate $\mathcal{R}^{\varepsilon}_{\mathrm{ini}}$ (2.6) and $\mathcal{S}^{\varepsilon}f$ (5.29), we verify that

[TABLE]

Remark 5.10*.*

From (5.29), we verify that

[TABLE]

Referring to Remark 5.2, the definition of $Qf$ ensures that the second term has order $\mathcal{O}(\varepsilon^{\alpha+1})$. Hence, $Qf$ ensures a lower constant in the second estimate in (5.38).

With these estimates, we are able to prove Theorem 2.8.

*Proof of Theorem 2.8. * We have to estimate both error terms the right-hand side of (5.31). Using Lemma 5.8 and (5.38), we verify that $\boldsymbol{\varphi}$ satisfies the estimate

[TABLE]

where $r(\alpha)=\alpha+2\lfloor\alpha/2\rfloor+2$. Hence, from the definition of $\boldsymbol{\mathcal{B}}^{\varepsilon}\tilde{u}$, we deduce that

[TABLE]

Next, applying Lemma 5.8 with \boldsymbol{\eta}^{\varepsilon}(t)=\boldsymbol{\mathcal{B}}^{\varepsilon}\tilde{u}(t)-\scalebox{2.1}[1.0]{[}{u^{\varepsilon}(t)}\scalebox{2.1}[1.0]{]} (see (5.38)), we obtain

[TABLE]

where we used that for $\alpha\geq 2$, $r(\alpha)\geq\alpha+4$. Using (5.39) and (5.40) in (5.31) proves estimate (2.17) and the proof of Theorem 2.8 is complete. $\square$

5.3 A symmetrized tensor product of symmetric positive definite matrices is positive definite

(proof of Lemma 3.1)

Lemma 3.1 states that the tensor $S^{2n}(\otimes^{n}a^{0})$ is positive definite, where $a^{0}$ is the homogenized tensor. We prove here that this property is true for any symmetric, positive definite matrix.

A first important result is the following.

Lemma 5.11**.**

Let $R\in{\mathrm{Ten}}^{2n}(\mathbb{R}^{d})$ be a positive definite tensor and let $A\in{\mathrm{Sym}}^{2}(\mathbb{R}^{d})$ be a symmetric, positive definite matrix. Then the tensor of ${\mathrm{Ten}}^{2n+2}(\mathbb{R}^{d})$ defined by $A_{i_{1}i_{2n+2}}R_{i_{2}\cdot\cdot i_{2n+1}}$ is positive definite.

Proof.

As $A$ is symmetric positive definite, the Cholesky factorization provides an invertible matrix $H$ such that $A=H^{T}H$. For $\xi\in{\mathrm{Sym}}^{n+1}(\mathbb{R}^{d})$, we thus have

[TABLE]

As $R$ is positive definite, the equality holds if and only if $H_{rj}\xi_{ji_{2}\cdot\cdot i_{n+1}}=0$ for all $r,i_{2},\ldots,i_{n+1}\in\{1,\ldots,d\}$. Let $i_{2},\cdots,i_{n+1}$ be arbitrarily fixed and denote $v_{j}=\xi_{ji_{2}\cdot\cdot i_{n+1}}$. Hence, we have $H_{rj}v_{j}=0$ for all $r$, which is equivalent to $H^{T}v=0$. As $H^{T}$ is regular, we obtain that $v=0$. We have proved that the equality in (5.41) implies $\xi=0$. Hence the tensor is positive definite and the proof of the lemma is complete. ∎

With Lemma 5.11 at hand, we are able to prove the following result, which implies Lemma 3.1.

Lemma 5.12**.**

If $A\in{\mathrm{Sym}}^{2}(\mathbb{R}^{d})$ is a symmetric, positive definite matrix, then the tensor $S^{2n}(\otimes^{n}A)\in{\mathrm{Sym}}^{2n}(\mathbb{R}^{d})$ is positive definite.

Proof.

We proceed by induction on $n$. The case $n=1$ is ensured by [5, Lemma 4.1] (it can also be deduced from Lemma 5.11). We assume that the result holds for $1,\ldots,n-1$ and prove it for $n$. Let $\xi\in{\mathrm{Sym}}^{n}(\mathbb{R}^{d})\backslash\{0\}$. First, assume that $n$ is odd. Then, the product $S^{2n}(\otimes^{n}A)\xi:\xi$ is composed of terms of the form

[TABLE]

i.e., one of the factor $A_{i_{r}i_{s}}$ share on index with both entities of $\xi$. The induction hypothesis combined with Lemma 5.11 ensure that all these terms are strictly positive and thus $S^{2n}(\otimes^{n}A)$ is positive definite. Second, we assume that $n$ is even. Then, the product $S^{2n}(\otimes^{n}A)\xi:\xi$ is composed of terms of two forms. First, there are terms of the form (5.42). By the same induction argument as before, they are strictly positive. Second, terms of the form

[TABLE]

Altogether, we verify that $S^{2n}(\otimes^{n}A)\xi:\xi>0$ and the proof of the lemma is complete. ∎

5.4 Proof of the new relation between the correctors (Theorem 3.5)

We prove the result for $r\geq 1$. Note that in the case $r=1$, we adopt the convention that empty sums vanish, i.e., $\sum_{k=1}^{0}x_{k}=0$. For the sake of clarity, we assume that $|Y|=1$ so that $(\cdot,1)_{Y}=\langle{\cdot}\rangle_{Y}$.

In Section 2.1, we explained that the cell problems (2.10) are well-posed if and only if $p^{j}=_{S}g^{j}$. Let us then replace $p^{j}$ by $g^{j}$ in the expression of the cell problems (2.10). For $r\geq 1$, let $\chi_{1},\ldots,\chi^{2r+1}$ be the $2r+1$ first zero mean correctors. We define the tensors $A^{k},B^{k},C^{k}\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$ as

[TABLE]

Note that the symmetry of $a$ ensures the following symmetry relations for $A^{k}$ and $C^{k}$:

[TABLE]

Furthermore, using the test function $w=\chi^{2r+2-k}_{i_{k+1}\cdot\cdot i_{2r+2}}$ in the cell problem for $\chi^{k}_{i_{1}\cdot\cdot i_{k}}$ in (2.10) and the symmetry of $a$, we obtain the following relations

[TABLE]

Define then the tensor $T\in{\mathrm{Ten}}^{2r+2}(\mathbb{R}^{d})$ as

[TABLE]

Using the definition of $\sigma^{k}$ and (5.43), we verify that

[TABLE]

where in the second equality we changed the index $k=2r-m$. Using (5.44), we decompose the tensor $T$ as

[TABLE]

where the tensors $U^{1}$, $U^{2}$ and $U^{3}$ are

[TABLE]

The rest of the proof relies on the following result.

Lemma 5.13**.**

The tensors $U^{i}$ defined above satisfy the following relations:

[TABLE]

where $h^{r}$ is the tensor defined in (3.9).

While the proof of (i) and (ii) is direct, the proof of (iii) requires preliminary work, done in the two following lemmas.

Lemma 5.14**.**

The tensors $D^{k}$ satisfy

[TABLE]

where for $r=1$ the sum on the right-hand side vanishes.

Lemma 5.15**.**

The tensor $H$ defined in (5.48) satisfies $H=_{S}h^{r}$, where $h^{r}$ is the tensor defined in (3.9).

*Proof of Lemma 5.14. * Let us denote z^{j,k}=\big{\langle}{g^{k-j}\otimes\chi^{j}\otimes\chi^{2r-k}}\big{\rangle}_{Y}. As $\chi^{0}=1$ and the correctors have zero mean, we verify that $z^{0,k}=z^{j,2r}=0$. Hence, $D^{k}$ can be written as $D^{k}=\sum_{j=1}^{k}z^{j,k}$. Splitting the sum and making the change of index $\ell=2r-j$, we write for $1\leq k\leq r-1$

[TABLE]

As we verify that $z^{2r-\ell,2r-k}=_{S}z^{k,\ell}$, we can write

[TABLE]

As $g^{s}$ is nonzero only for even index (see (5.26)), we have for all $s$, $(-1)^{s}g^{s}=g^{s}$. In particular, for all $s=k-\ell$, we have $(-1)^{k}g^{k-\ell}=(-1)^{\ell}g^{k-\ell}$ which ensures that $(-1)^{k}z^{k,\ell}=(-1)^{\ell}z^{k,\ell}$. Hence, changing the summation order, we find

[TABLE]

Furthermore, we verify that

[TABLE]

Summing the two last equalities we obtain (5.48) and the proof is complete. $\square$

*Proof of Lemma 5.15. * Lemma 2.4 implies that $g^{2s+1}=_{S}0$. Hence, the only nonvanishing terms in the double sum

[TABLE]

are $\{1\leq m,n\leq r:m+n\text{ is even}\}$. These terms are exactly given by

[TABLE]

We verify that the expression of $h^{r}$ in (3.9) corresponds precisely to the sum on these two sets of index, i.e., $H=_{S}h^{r}$. $\square$

We now prove Lemma 5.13.

*Proof of Lemma 5.13. * Let us start by proving (i). Using the definition of $\sigma^{k}$, we have

[TABLE]

where in the second equality we changed the index $k=2r+1-m$. As $(-1)^{2r+1-k}=(-1)^{k+1}$, we obtain equality (i).

Next, we prove (ii). Using again the definition of $\sigma^{k}$, we find

[TABLE]

where we changed the index $k=2r+2-m$. The symmetry (5.43) implies that $C^{2r+1-k}=_{S}C^{k-1}$. Hence, noting that $(-1)^{2r+2-k}=(-1)^{k}$, we obtain (ii).

Finally, we prove (iii). With the definition of $\sigma^{k}$, we find

[TABLE]

where we made the changes of index $k=n-1$ and $k=2r+1-m$, respectively. As $(-1)^{2r+2-k}=(-1)^{k}$, using Lemma 5.14 we obtain

[TABLE]

As $\chi^{0}=1$ and the correctors have zero mean, we verify that $D^{0}=0$. Using Lemma 5.15, we obtain $U^{3}=-h^{r}$ and (iii) is proved. $\square$

With Lemma 5.13 at hand we can prove Theorem 3.5.

*Proof of Theorem 3.5. * Combining (5.46), (5.47) and Lemma 5.13 yields

[TABLE]

We verify that the tensors $g^{2r}$ and $k^{r}$, defined in (3.7) and (3.9), respectively, can be written as

[TABLE]

Hence, we deduce from (5.49) that

[TABLE]

This equality matches the decomposition (3.8) and the proof of Theorem 3.5 is complete. $\square$

6 Conclusion

In this paper, we presented a family of effective equations for wave propagation in periodic media for arbitrary timescales. In particular, for any given $\alpha\geq 0$, our main result (Theorem 2.8) ensures the effective solutions to be close to $u^{\varepsilon}$ in the ${\mathrm{L}^{\infty}}(0,\varepsilon^{-\alpha}T;W)$ norm. As emphasized, the effective equations are well-posed without requiring any regularization process. In addition, we described a numerical procedure to compute the effective tensors of equations in the family. We showed that the computational cost of this procedure can be significantly reduced by using a new relation between the correctors. This relation should also be used to compute the tensors of the alternative effective model available in the literature.

One question that is raised is how to find the best effective equation in the family? More precisely, can we find a criterion to find an optimal equation in the family and can we build such an equation explicitly? These interrogations also concern the effective models from [9] and [13]. Answers to these questions are exciting topics left for future research.