\documentclass{TEMA}

\usepackage[latin1]{inputenc} 
\usepackage{amsmath,amssymb,amsthm,mathrsfs,amsfonts}
\usepackage[latin1]{inputenc}
%\usepackage{dsfont}
\usepackage{stmaryrd}
\usepackage{bm}
\usepackage{graphicx}
\usepackage{float}
\usepackage[english]{babel}
\usepackage{epsfig}
%\usepackage{refcheck}
\begin{document}

%********************************************************
\title
    {Homogenization of a continuously microperiodic multidimensional medium\thanks{Financial support by project CAPES No. 88881.030424/2013-01 is gratefully acknowledged; Paper presented at the National Congress of Applied and Computational Mathematics 2016.}}

\author
   {M. P. Lima%
     \thanks{marcos.lima@ufrgs.br; MSc in Mathematical Modeling}\,,
     Institute of Mathematics and Statistics,
     UFRGS - Universidade Federal do Rio Grande do Sul,  91509-900 Porto Alegre, RS, Brazil
     \\ \\
     L. D. Pérez-Fernández%
     \thanks{leslie.fernandez@ufpel.edu.br; PhD in Mathematics}\,,
     Department of Mathematics and Statistics, UFPel - Universidade Federal de Pelotas,
96160-000 Capão do Leão, RS, Brazil.
	\\	\\
J. Bravo-Castillero%
\thanks{jbravo@matcom.uh.cu; PhD in Mathematics}\,,
     Facultad de Matemática y Computación, UH - Universidad de la Habana, CP 10400 Ciudad de la Habana, Cuba.	
	}




\criartitulo
%\runningheads {Lima, Pérez-Fernández and Bravo-Castillero}{Homogenization of a Microperiodic Medium}

\begin{abstract}
{\bf Resumo}. The asymptotic homogenization method is applied to obtain formal asymptotic solution and the homogenized solution of a Dirichlet boundary-value problem for an elliptic equation with rapidly oscillating coefficients. The proximity of the formal asymptotic solution and the homogenized solution to the exact solution is proved, which provides the mathematical justification of the homogenization process. Preservation of the symmetry and positive-definiteness of the effective coefficient in the homogenized problem is also proved. An example is presented in order to illustrate the theoretical results.

{\bf Palavras-chave}. Microperiodic medium, effective behavior, asymptotic homogenization method.
\end{abstract}



 
\newsec{Introduction}

The heterogeneous medium is formed by distributions of occupied domains: (a) by different homogeneous materials called phase, thus constituting a composite; or (b) of the same material in different states, such as a polycrystal \cite{Torquato2002}, or a functionally graded material \cite{Sadd2005}. Heterogeneous medium abound in nature and in manufactured products. For example: among the natural medium that have heterogeneity are the bone, atmosphere, soil, sandstone, wood, lungs, vegetable and animal tissues, cell aggregates, tumors; (solid, granular or particulate, fibrous, and combinations thereof), and cell solids, gels, foams, metal alloys, microemulsions, ceramic and block copolymers \cite{Torquato2002}. The theoretical prediction of mechanical, electromagnetic, and transport properties of heterogeneous materials has a long and venerable history, attracting the attention of science icons, including Maxwell \cite{Maxwell1873}, Rayleigh \cite{Rayleigh1892} and Einstein \cite{Einstein1906}. Generally, the physical phenomena of interest associated with such properties occur in the "microscale", called generically so, because it can be in the order of tenths of nanometers (gels) up to the order of meters (geological processes). Such medium exhibit separation of structural scales, which is characterized by the small parameter $\varepsilon = l/L$, $0<\varepsilon\ll1$, where $l$ and $L$ are, respectively, the characteristic lengths of the micro- and macroscale. Also, these medium are assumed to satisfy the continuum hypothesis, that is, $l$ is much greater than the characteristic length of the molecular scale. With such assumptions, micro-heterogeneous medium are regarded as a continuum at the microscale, so they can be characterized via effective properties \cite{Torquato2002}. More precisely, the hypothesis of equivalent homogeneity states that, at the macroscale, an heterogeneous medium is physically equivalent to an ideally homogeneous medium in such a way that the effective properties of the former are the properties of the latter \cite{Bakhvalov1989}. The process of obtaining the effective behavior of micro-heterogeneous medium is called homogenization. The periodicity is characterized by periodic replication of a recurrent element $Y$ called basic cell. 

Phenomena occurring within a micro-heterogeneous medium are described by initial/boundary-value problems whose differential equations have rapidly oscillating coefficients, whereas the equations of the problems for the equivalent homogeneous medium have constant coefficients \cite{Bakhvalov1989}. Thus, the hypothesis of equivalent homogeneity will be valid if the solution $u^\varepsilon$ of a problem for the heterogeneous medium is $\varepsilon$-near of the solution $u_0$ of the problem for the equivalent homogeneous medium with respect to some norm, that is, $\left\|u^\varepsilon - u_0\right\| = \mathcal{O}(\varepsilon)$.

We can find various applications of homogenization theory in literature, for example: obtaining properties coupled nonexistent in the constituents (magnetoelectric
\cite{Benveniste1995}, pyroelectric and 
pyromagnetic \cite{Bravo2015}); topology optimization \cite{Bendsoe2003}; optimal design of heterogeneous materials \cite{Torquato2010}; biomechanics of bone \cite{Parnell2009}, prediction of structural failures \cite{Perez2014}; propagation of seismic waves \cite{Capdeville2010}; physics of nuclear reactors \cite{Allaire1999}; transport of a chemical species \cite{NG2006}.

In this contribution, we study the following elliptic problem: Find $u^\varepsilon\in C^2(\overline{\Omega})$, $\overline{\Omega} = [0,1]^d$, such that
\begin{equation}\label{probmult}
%\left\{
%\begin{array}{l}
%\displaystyle
\mathcal{L}^\varepsilon u^\varepsilon \equiv \frac{\partial}{\partial x_j}\left(a^\varepsilon_{jl}(x)\frac{\partial u^\varepsilon}{\partial x_l}\right) = f(x), \mbox{ } x\in\Omega, \quad 
u^\varepsilon = 0, \mbox{ } x\in\partial\Omega,
%\end{array}
%\right.
\end{equation}
where, $j,l=1,2,...,d$, $a^\varepsilon_{jl}\in C^1(\overline{\Omega})$ are $\varepsilon Y$-periodic, $Y = [0,1]^d$, symmetric and positive-definite, that is, $a^\varepsilon_{jl}(x) = a^\varepsilon_{lj}(x)$ for all $x\in\Omega$ (symmetry) and $\forall \eta\in\mathbb{R}^d, \exists c>0 : \forall x\in\mathbb{R}^d, a^\varepsilon_{jl}(x)\eta_j\eta_l\geq c\eta_l\eta_l$ (positive definiteness), respectively. In (\ref{probmult}), and throughout the work, Einstein's notation for the sum over repeated indexes is adopted. 

Problem (\ref{probmult}) models, for instance, the distribution of a stationary temperature field $u^\varepsilon$ over a periodic micro-heterogeneous conductive medium with thermal conductivity $a^\varepsilon_{jl}(x) = a_{jl}(x/\varepsilon)$ in the presence of a distributed heat source $f$. Notice that this problem is, in general, hard of to solve analytically, while a direct numerical approach requires a very fine discretization of the domain in order to capture the rapidly oscillating behavior of the coefficients $a^\varepsilon_{jl}(x)$, which increases considerably the computational cost and compromises the convergence of the adopted numeric method. An effective and mathematically rigorous alternative is the Asymptotic Homogenization Method (AHM) \cite{Bakhvalov1989}. The AHM proposes a formal asymptotic solution (f.a.s.) of the original problem (\ref{probmult}), that is, a two-scale asymptotic series of powers of $\varepsilon$,
\begin{equation}\label{multsaf0}
u^{(\infty)}(x,\varepsilon)=\sum_{i=0}^{\infty}\varepsilon^i u_i(x,y), \quad y=\frac{x}{\varepsilon}.
\end{equation}
Note that, as $\mathcal{L}^\varepsilon$ is a linear operator, the asymptotic series (\ref{multsaf0}) is also an asymptotic expansion of the exact solution of the problem: $u^\varepsilon \sim u^{(\infty)}$. The unknown coefficients $u_i\in C^2(\overline{\Omega}\times\mathbb{R}^d)$ of the powers of $\varepsilon$, which are $Y$-periodic functions with respect to the so-called local variable $y$, are obtained as the solutions of a recurrent sequence of problems which results from subsituting (\ref{multsaf0}) into the original problem (\ref{probmult}).

The work is structured as follows: in section 2, we develop the AHM in detail to build the f.a.s.; in section 3, we mathematically justify the AHM, i.e., we prove that $\|u^\varepsilon-u_0\|_{H^1_0(\Omega)}=\mathcal{O}(\sqrt{\varepsilon)}$, where $H^1_0(\cdot)$ is the space of null-trace square-integrable functions whose generalized first-order derivatives are square-integrable; in sections 4 and 5, we show some properties preserved on apply of the AHM; an example is solved in section 6, showing the development  analytical (separation of variables) and numerical (finite difference method) applied; finally, in section 7 we describes the conclusions about development realized in the previous sections. 

\newsec{Aplication of the AHM}

Following \cite{Bakhvalov1989}, consider the f.a.s 
\begin{equation}\label{multsaf}
u^{(2)}(x,\varepsilon)=u_0(x,y)+\varepsilon u_1(x,y)+\varepsilon^2u_2(x,y),\;\;y=\frac{x}{\varepsilon},
\end{equation}
for problem (\ref{probmult}). By substituting (\ref{multsaf}) into the equation of (\ref{probmult}), applying the chain rule 
\begin{equation*}\left.\frac{\partial \Psi^\varepsilon}{\partial x_j}=\frac{\partial \Psi}{\partial x_j}+\frac{1}{\varepsilon}\frac{\partial \Psi}{\partial y_j}\right\rvert_{y=x/\varepsilon},
\end{equation*}
where $\displaystyle \Psi^\varepsilon = \Psi\left(x,x\varepsilon\right)$, and grouping by power $\varepsilon$, we have
\begin{multline}\label{expr1}
\mathcal{L}^\varepsilon
u^{(2)}-f(x) = \varepsilon^{-2}\mathcal{L}_{yy}u_0+\varepsilon^{-1}\left(\mathcal{L}_{xy}u_0+\mathcal{L}_{yx}u_0+\mathcal{L}_{yy}u_1\right)\\+\varepsilon^0\left(\mathcal{L}_{xx}u_0
+\mathcal{L}_{xy}u_1+\mathcal{L}_{yx}u_1+\mathcal{L}_{yy}u_2-f(x)\right)+\mathcal{O}(\varepsilon),
\end{multline}
where the differential operators $\mathcal{L}_{\alpha\beta}$, $\alpha,\beta\in\{x,y\}$, are defined as
\begin{equation*}
\mathcal{L}_{\alpha\beta} \equiv \frac{\partial}{\partial \alpha_j}\left(a_{jl}(y)\frac{\partial }{\partial \beta_l}\right).
\end{equation*}


Thus, in order to $u^{(2)}$ given by (\ref{multsaf}) be a f.a.s of problem (\ref{probmult}), that is, for the right-hand side of (\ref{expr1}) be assimptotically null, its unknown coefficients $u_i$ must satisfy the following recurrent equations:
\begin{eqnarray}
\varepsilon^{-2}&:&\mathcal{L}_{yy}u_0=0 \label{u0mult}\\
\varepsilon^{-1}&:&\mathcal{L}_{yy}u_1=-\mathcal{L}_{xy}u_0-\mathcal{L}_{yx}u_0\label{u1mult}\\
\varepsilon^0&:&\mathcal{L}_{yy}u_2=-\mathcal{L}_{xx}u_0
-\mathcal{L}_{xy}u_1-\mathcal{L}_{yx}u_1+f(x).\label{u2mult}
\end{eqnarray} 

The following Lemma \cite{Bakhvalov1989} guarantees the existence of $Y$-periodic solutions $u_i$ of (\ref{u0mult})-(\ref{u2mult}): 
\begin{lemmaTEMA}\label{lemageral}
Let $a_{jl}(y)$ and $F(y)$ be $Y$-periodic differentiable functions, such that $a_{jl}(y)$ is symmetric and positive definite. Then, a necessary and sufficient condition for a $Y$-periodic solution of the equation
\begin{equation}\label{eqLema}
\mathcal{L}_{yy} N \equiv\frac{\partial }{\partial y_j}\left(a_{jl}(y)\frac{\partial N}{\partial y_l}\right)=F(y),
\end{equation}
to exist is that $F(y)$ has null mean value, that is,
\begin{equation*}
\langle F(y)\rangle\equiv\int_Y F(y)dy=0,
\end{equation*}
where $\langle \cdot \rangle$ is the local averaging operator over the periodic cell.
\end{lemmaTEMA}

\begin{remTEMAi}\label{obs}
Note that the $Y$-periodic solution $N$ of (\ref{eqLema}) is unique up to an additive constant, that is, $N(y)=\widetilde{N}(y)+C$, where $\widetilde{N}(y)$ is the null-average $Y$-periodic solution of (\ref{eqLema}), $\langle \widetilde{N}\rangle=0$, and $C$ is an arbitrary constant.
\end{remTEMAi}

Remark \ref{obs} of Lemma \ref{lemageral} applied to (\ref{u0mult}) implies that $u_0$ does not depend on $y$:
\begin{equation}\label{u0r}
u_0(x,y)=u_0(x).
\end{equation}

By substituting (\ref{u0r}) into (\ref{u1mult}), we have
\begin{equation}\label{u1r}
\mathcal{L}_{yy}u_1=-\mathcal{L}_{yx}u_0,
\end{equation}
which, by applying Lemma \ref{lemageral}, implies that there exists $u_1$, $Y$-periodic solution of (\ref{u1r}), as $a_{jl}(y)$ is $Y$-periodic, so $\langle\partial a_{jl}/\partial y_j\rangle = 0$. In order to solve (\ref{u1r}), we seek its solution $u_1$ by separation of variables as 
\begin{equation}\label{sep}
u_1(x,y)=N_p(y)\frac{\partial u_0}{\partial x_p}.
\end{equation}
Substitution of (\ref{sep}) into (\ref{u1r}) yields
\begin{equation}\label{u1r0}
\frac{\partial}{\partial y_j}\left(a_{jl}(y)\frac{\partial N_p}{\partial y_l}\right)\frac{\partial u_0}{\partial x_p}=-\frac{\partial a_{jl}}{\partial y_j}\frac{\partial u_0}{\partial x_l},
\end{equation}
which, by replacing the index $l$ by $p$ in the right-hand side, leads to
\begin{equation}\label{u1r1}
\frac{\partial}{\partial y_j}\left(a_{jp}(y)+a_{jl}(y)\frac{\partial N_p}{\partial y_l}\right)\frac{\partial u_0}{\partial x_p}=0.
\end{equation}
Assume that $\partial u_0/\partial x_p\neq0$. Then, (\ref{u1r1}) is satisfied by $N_p(y)$, $p=1,\ldots,d$, the $Y$-periodic solutions of the so-called local problems:
\begin{equation}\label{problocalmult}
\frac{\partial}{\partial y_j}\left(a_{jp}(y)+a_{jl}(y)\frac{\partial N_p}{\partial y_l}\right)=0, \mbox{ } y\in Y, \quad \langle N_p(y)\rangle = 0.
\end{equation}
Lemma \ref{lemageral} guarantees the existence of $N_p(y)$, $Y$-periodic solutions of $\mathcal{L}_{yy}N_p=F$ with $F=-\partial a_{jp}/\partial y_j$, and their uniqueness is ensured by condition $\langle N_p(y)\rangle=0$.

Application of Lemma \ref{lemageral} provides a necessary and sufficient condition for the existence of $u_2$, $Y$-periodic solution of (\ref{u2mult}):
\begin{equation}\label{aplicLema}
\langle\mathcal{L}_{xx}u_0 + \mathcal{L}_{xy}u_1 + \mathcal{L}_{yx}u_1\rangle - f(x)=0.
\end{equation}
Substitution of (\ref{sep}) into (\ref{aplicLema}) followed by the appropriate index replacement yields
\begin{equation}\label{aplicLema1}
\left\langle\frac{\partial}{\partial y_l}\left(a_{lj}(y)N_p(y)\right)+a_{jl}(y)\frac{\partial N_p}{\partial y_l}+a_{jp}(y)\right\rangle\frac{\partial^2 u_0}{\partial x_j \partial x_p}=f(x).
\end{equation}
Notice that, $\left\langle\partial\left(a_{lj}(y)N_p(y)\right)/\partial y_l\right\rangle=0$, as $a_{jp}(y)$ and $N_p(y)$ are $Y$-periodic. Thus, (\ref{aplicLema1}) can be rewritten as
\begin{equation}\label{eqhomomult}
\widehat{a}_{jp}\frac{\partial^2 u_0}{\partial x_j \partial x_p}=f(x),
\end{equation}
which is the equation of the so-called homogenized problem, and
\begin{equation}\label{coefmult}
\widehat{a}_{jp}=\left\langle a_{jl}(y)\frac{\partial N_p}{\partial y_l}+a_{jp}(y)\right\rangle
\end{equation}
are the so-called effective coefficients. The homogenized problem is defined as: find $u_0\in C^2(\overline{\Omega})$ such that
\begin{equation}\label{probhomomult}
\widehat{a}_{jp}\frac{\partial^2 u_0}{\partial x_j \partial x_p} = f(x), \mbox{ } x\in\Omega, \quad u_0(x) = 0, \mbox{ } x\in\partial\Omega
\end{equation}
where $\widehat{a}_{jp}$ are the effective coefficients given by (\ref{coefmult}). Also, notice that the existence of $u_0(x)$ solution of (\ref{probhomomult}) guarantees the existence of $u_2(x,y)$ $Y$-periodic solution of (\ref{u2mult}). By substituting (\ref{sep}) and (\ref{eqhomomult}) into (\ref{u2mult}), we have 
\begin{equation*}\label{conprox}
\mathcal{L}_{yy}u_2=-\left(\frac{\partial}{\partial y_l}\left(a_{lj}(y)N_p(y)\right)+a_{jl}(y)\frac{\partial N_p}{\partial y_l}+a_{jp}(y)-\widehat{a}_{jp}\right)\frac{\partial^2 u_0}{\partial x_j \partial x_p}.
\end{equation*}
This show that the solution $u_2$ can be sought as
\begin{equation}\label{u2sep}
u_2(x,y)=N_{pq}(y)\frac{\partial^2 u_0}{\partial x_p\partial x_q},
\end{equation}
where $N_{pq}(y)$ is the $Y$-periodic solution of the so-called second local problem, which is obtained by substituting (\ref{eqhomomult}) and (\ref{u2sep}) into (\ref{u2mult}) and following the same steps as for obtaining $u_1$. The second local problem is
\begin{equation}\label{problocalmult2}
\mathcal{L}_{yy}N_{pq}=\widehat{a}_{pq}-T_{pq}, \mbox{ } y\in Y, \quad \langle N_{pq}\rangle=0,
\end{equation}
where
\begin{equation*}
T_{pq}(y)=a_{pq}(y)+a_{ql}(y)\frac{\partial N_p}{\partial y_l}+\frac{\partial}{\partial y_l}\left(a_{lq}(y)N_p(y)\right).
\end{equation*}
The existence of the $Y$-periodic solution $N_{pq}$ of problem (\ref{problocalmult2}) is guaranteed by Lemma \ref{lemageral}.

Finally, the second-order f.a.s. $u^{(2)}(x,\varepsilon)$ of the original problem (\ref{probmult}) is
\begin{equation}\label{saf2}
u^{(2)}(x,\varepsilon)=u_0(x)+\varepsilon N_p\left(\frac{x}{\varepsilon}\right)\frac{\partial u_0}{\partial x_p}+\varepsilon^2 N_{pq}\left(\frac{x}{\varepsilon}\right)\frac{\partial^2 u_0}{\partial x_p
\partial x_q}.
\end{equation}
On the other hand, we show in the next section that the first-order f.a.s $u^{(1)}(x,\varepsilon)$ consisting of the two first terms of $u^{(2)}(x,\varepsilon)$ in (\ref{saf2}) yields a good approximation of the exact solution $u^\varepsilon(x)$ of the original problem (\ref{probmult}). 


\newsec{Proximity}

In general, the f.a.s $u^{(2)}$, constructed to approximate the exact solution $u^\varepsilon$ of the original problem (\ref{probmult}), does not satisfy the boundary condition $u^\varepsilon|_{\partial \Omega}=0$. As the solution $u_0$ of the homogenized problem (\ref{probhomomult}) satisfies the boundary condition $u_0|_{\partial \Omega}=0$, it follows that the error of $u^{(2)}|_{\partial \Omega}$ is nonzero of order $\varepsilon$. In order to overcome such a situation, we multiply $\varepsilon u_1+\varepsilon^2 u_2$ by a function $\chi(x)$ to force the f.a.s. to satisfy exactly the boundary condition, but this adds some new terms to the error. However, these terms are evaluated as being of order $\sqrt{\varepsilon}$ in the norm of $H^{-1}(\Omega)$, that is, the error in form $f_0+\partial f_i/\partial x_i$, with $\|f_0\|_{L^2(\Omega)}$ and $\|f_i\|_{L^2(\Omega)}$ of order $\sqrt{\varepsilon}$, where $L^2(\cdot)$ is space of square-integrable functions, so we can apply an estimate for the solution in $H^1_0(\Omega)$.

Specifically, if $N_p$ are solutions of the local problems (\ref{problocalmult}), $N_{pq}$ are solutions of the second local problems (\ref{problocalmult2}), $u_0$ is the solution of the homogenized problem (\ref{probhomomult}), and $u_0(x)\in C^2 (\overline{\Omega})$, then the f.a.s (\ref{saf2}) is solution of the equation
\begin{equation*}
\mathcal{L}^\varepsilon u^{(2)}-f(x)=\varepsilon r^{(2)}(x,\varepsilon)
\end{equation*}
where $\varepsilon r^{(2)}(x,\varepsilon)$ is the error of taking $u^{(2)}$ as the solution of the original problem:
\begin{equation*}
\varepsilon r^{(2)}(x,\varepsilon)=\varepsilon\left(\mathcal{L}_{xy}u_2+\mathcal{L}_{yx}u_2+\mathcal{L}_{xx}u_1\right)+\varepsilon^2
\mathcal{L}_{xx}u_2,\;y=\frac{x}{\varepsilon}.
\end{equation*}
Notice that, by assuming that $u_0\in\mathcal{C}^4(\overline{\Omega})$, the application of Weierstrass theorem \cite{Kudriavtsev1983} ensures that $\|r^{(2)}(x,\varepsilon)\|_{L^2(\Omega)}=\mathcal{O}(1)$. However, the error of $u^{(2)}$ in the boundary $\partial \Omega$ is of order $\varepsilon$:
\begin{equation*}
u^{(2)}|_{\partial \Omega}=\varepsilon u_1|_{\partial \Omega}+\varepsilon^2 u_2|_{\partial \Omega}.
\end{equation*}

Let $\chi(x)$ be an infinitely differentiable function with support contained inside the $\varepsilon$-neighbourhood of $\partial\Omega$ and such that $\chi|_{\partial \Omega}=1$, $\left|\chi\right|\leq 1$, and $\| \varepsilon \partial \chi/\partial x\|_{C(\overline{\Omega})}\leq c_1$, where $c_1$ is a constant independent of $\varepsilon$.

Consider the modified second-order f.a.s
\begin{align}
\widetilde{u}^{(2)}(x,\varepsilon) = &u^{(2)}(x,\varepsilon)-\chi(x)\left(\varepsilon u_1(x,y)+\varepsilon^2 u_2(x,y)\right)\nonumber\\
=&u_0(x)+\varepsilon\left(1-\chi(x)\right)u_1(x,y)+\varepsilon^2\left(1-\chi(x)\right)u_2(x,y), \mbox{ } y=\frac{x}{\varepsilon}\label{aproxcc}
\end{align}
that satisfies exactly the boundary conditions of the original problem (\ref{probmult}). In order to assess how good $\widetilde{u}^{(2)}$ is as an approximation of $u^\varepsilon$, we have to estimate the error given by the right-hand side of
\begin{equation}\label{probexatmult}
\mathcal{L}^\varepsilon \widetilde{u}^{(2)}-f(x)=\varepsilon r^{(2)}(x,\varepsilon)-\frac{\partial}{\partial x_j}\Phi_j(x,\varepsilon),
\end{equation}
where
\begin{equation}\label{novoerro}
\Phi_j(x,\varepsilon)=a_{jl}\left(\frac{x}{\varepsilon}\right)\left[\varepsilon\frac{\partial (\chi u_1)}{\partial x_l}+\varepsilon^2 \frac{\partial (\chi u_2)}{\partial x_l}\right].
\end{equation}

First, by the hypotheses on $\chi$, we have that $\Phi_j(x, \varepsilon)$ is such that $\left| \Phi_j(x,\varepsilon) \right|\leq c_2$ for every $x\in \overline{\Omega}$, where $c_2$ is a constant independent of $\varepsilon$, and its support is contained inside the $\varepsilon$-neighbourhood of $\partial\Omega$.
For instance, for $d=3$, we have that
\begin{equation}\label{defchi}
\Phi_j(x,\varepsilon)=0,\; x\in [\varepsilon,1-\varepsilon]^3,
\end{equation}
and $\|\Phi_j(x,\varepsilon)\|_{L^2(\Omega)}=\mathcal{O}(\sqrt{\varepsilon})$. Indeed, by (\ref{defchi}) we have that
\begin{equation}\label{novest}
\|\Phi_j(x,\varepsilon)\|^2_{L^2(\Omega)}=\int_{\Omega} \Phi^2_j(x,\varepsilon)dx= \int_{\Omega\setminus[\varepsilon,1-\varepsilon]^3} \Phi^2_j(x,\varepsilon)dx.
\end{equation}
Consider the domains $\Omega^*=\Omega\setminus[\varepsilon,1-\varepsilon]^3\ni x$ and $\Omega^*_\varepsilon\ni x/\varepsilon$ such that $\Omega^*=\varepsilon\Omega^*_\varepsilon$. By substituting (\ref{sep}) and (\ref{u2sep}) into (\ref{novest}) with (\ref{novoerro}), and differentiating, we obtain
\begin{multline*}
\int_{\Omega^*}{\left[a_{jl}\left(\frac{x}{\varepsilon}\right)\left(\varepsilon\frac{\partial (\chi u_1)}{\partial x_l}+\varepsilon^2 \frac{\partial (\chi u_2) }{\partial x_l}\right)\right]^2}dx\\
=\int_{\Omega^*}\left[a_{jl}\left(\frac{x}{\varepsilon}\right)\left(\varepsilon N_p\left(\frac{x}{\varepsilon}\right)\left(\frac{\partial \chi}{\partial x_l}\frac{\partial u_0}{\partial x_p} + \chi(x)\frac{\partial^2 u_0 }{\partial x_p \partial x_l} \right) \right.\right.\\
\left.\left.+\varepsilon^2 N_{pq}\left(\frac{x}{\varepsilon}\right) \left(\frac{\partial \chi}{\partial x_l}\frac{\partial^2 u_0}{\partial x_p \partial x_q}+\chi(x) \frac{\partial^3 u_0}{\partial x_p \partial x_q \partial x_l}\right)\right)\right]^2 dx.
\end{multline*}
By Weierstrass theorem, there are constants $B_0$, $B_1$ e $B_2$ such that $|a_{jl}|\leq B_0$, $|N_p|\leq B_1$ and $|N_{pq}|\leq B_2$ for every $x/\varepsilon\in \Omega^*_\varepsilon$. By taking $B=\max\{B_0,B_1,B_2\}$, we have that
\begin{multline*}
\int_{\Omega^*}\left[a_{jl}\left(\frac{x}{\varepsilon}\right)\left(\varepsilon\frac{\partial (\chi u_1)}{\partial x_l}+\varepsilon^2 \frac{\partial (\chi u_2) }{\partial x_l}\right)\right]^2dx\\
\leq B^4 \int_{\Omega^*} \left(\left|\varepsilon \frac{\partial \chi}{\partial x_l}\right|\left|\frac{\partial u_0}{\partial x_p}\right|+\varepsilon \left|\chi\right|\left|\frac{\partial^2 u_0}{\partial x_p\partial x_l}\right|\right.\\
\left.+\varepsilon\left|\varepsilon \frac{\partial \chi}{\partial x_l}\right|\left|\frac{\partial^2 u_0}{\partial x_p \partial x_q}\right|+\varepsilon^2 |\chi|\left|\frac{\partial^3 u_0}{\partial x_p \partial x_q \partial x_l}\right|\right)^2dx.
\end{multline*}
By Weierstrass theorem, there are constants $C_0$, $C_1$, $C_2$ and $C_3$ such that 
\begin{equation*}
\left|\varepsilon \frac{\partial \chi}{\partial x_l}\right|\leq C_0,\;\;\left|\frac{\partial u_0}{\partial x_p}\right|\leq C_1,\;\;\left|\frac{\partial^2 u_0}{\partial x_p \partial x_l}\right|\leq
C_2,\;\;\left|\frac{\partial^3 u_0}{\partial x_p \partial x_q \partial x_l}\right|\leq C_3,\;\;\forall x\in \Omega^*.
\end{equation*}
By recalling that $|\chi|\leq 1$ and taking $C=\max\{1,C_0,C_1,C_2,C_3\}$, we have that
\begin{equation*}
\int_{\Omega^*}\left[a_{jl}\left(\frac{x}{\varepsilon}\right)\left(\varepsilon\frac{\partial (\chi u_1)}{\partial x_l}+\varepsilon^2 \frac{\partial (\chi u_2) }{\partial x_l}\right)\right]^2dx\leq B^4 C^4 (1+\varepsilon)^4\int_{\Omega^*}dx.
\end{equation*}
Then, from (\ref{novest}) it follows that
\begin{equation*}
\|\Phi_j(x,\varepsilon)\|^2_{L^2(\Omega)}\leq 2B^4 C^4 (1+\varepsilon)^4 |\Omega^*|=2B^4 C^4 (1+\varepsilon)^4 D(\varepsilon),
\end{equation*}
where $D(\varepsilon)=3\varepsilon-6\varepsilon^2+4\varepsilon^3$ is of order $\mathcal{O}(\varepsilon)$. Therefore,
\begin{equation*}
\|\Phi_j(x,\varepsilon)\|_{L^2(\Omega)}\leq \sqrt{2}B^2 C^2 (1+\varepsilon)^2 \sqrt{D(\varepsilon)},
\end{equation*}
which proves that $ \|\Phi_j(x,\varepsilon)\|_{L^2(\Omega)}=\mathcal{O}(\sqrt{\varepsilon})$.

Now, subtracting of the equation of problem (\ref{probmult}) from (\ref{probexatmult}) yields
\begin{equation*}
\frac{\partial}{\partial x_j}\left(a^\varepsilon_{jl}(x)\frac{\partial}{\partial x_l}(\widetilde{u}^{(2)}-u^\varepsilon)\right)=\varepsilon r^{(2)}(x,\varepsilon)-\frac{\partial}{\partial x_j}\Phi_j(x,\varepsilon),
\end{equation*}
from which, recalling that $(\widetilde{u}^{(2)}-u^\varepsilon)|_{\partial \Omega}=0$, the maximum principle for elliptic equations \cite{Bakhvalov1989} provides the estimate
\begin{equation}\label{modsafdif}
\|\widetilde{u}^{(2)}-u^\varepsilon\|_{H^1_0(\Omega)}\leq c_3\left(\varepsilon\|r^{(2)}(x,\varepsilon)\|_{L^2(\Omega)}+\sum_{j=1}^3\|\Phi_j(x,\varepsilon)\|_{L^2(\Omega)}\right),
\end{equation}
where constant $c_3$ is independent of $\varepsilon$. Recall that, for $u_0\in C^4(\overline{\Omega})$, we have $\|r^{(2)}(x,\varepsilon)\|_{L^2(\Omega)}=
\mathcal{O}(1)$. Then, (\ref{modsafdif}) implies that the modified second-order f.a.s (\ref{aproxcc}) satisfies the relation
\begin{equation*}
\|\widetilde{u}^{(2)}-u^\varepsilon\|_{H^1_0(\Omega)}=\mathcal{O}(\sqrt{\varepsilon})
\end{equation*}
and, in particular, 
\begin{equation*}
\|\widetilde{u}^{(1)}-u^\varepsilon\|_{H^1_0(\Omega)}=\mathcal{O}(\sqrt{\varepsilon}),
\end{equation*}
where $\widetilde{u}^{(1)}(x,\varepsilon)=u^{(1)}-\varepsilon\chi u_1$ is the related first-order modified f.a.s

Finally, in order to show the proximity between $u_0$ and $u^\varepsilon$, it suffices to prove, following the procedure described above, that
\begin{equation*}
\|\widetilde{u}^{(1)}-u_0\|_{H^1_0(\Omega)}=\mathcal{O}(\sqrt{\varepsilon}).
\end{equation*}
Therefore, from the Minkowski inequality's \cite{EVANS}
\begin{equation*}
\|u^\varepsilon-u_0\|_{H^1_0(\Omega)}\leq\|u^\varepsilon-\widetilde{u}^{(1)}\|_{H^1_0(\Omega)}+\|\widetilde{u}^{(1)}-u_0\|_{H^1_0(\Omega)},
\end{equation*}
it follows that
\begin{equation*}
\|u^\varepsilon-u_0\|_{H^1_0(\Omega)}=\mathcal{O}(\sqrt{\varepsilon}).
\end{equation*}

An alternative approach to estimate $u^\varepsilon -\tilde{u}^{(1)}$ and  $\tilde{u}^{(1)}-u_0$ is through an estimative for the solution of non homogeneous Dirichlet problem provided in \cite{Cioranescu}. In that book, such an approach was applied and similar results were obtained, that is, the $\mathcal{O}(\sqrt{\varepsilon})$ proximity for homogeneous boundary condition. Other approaches can be applied to know to proximity between the $u^\varepsilon$ and $u_0$ solutions, for instance, proposes asymptotic expansions of which its coefficients are convolution forms of Green function and the source term of the elliptic equation, so from the representation of the elliptic solution by the Green function, we estimate the proximity of the expansion truncated and the exact solution \cite{Choe2003}. 

Other two important results of the AHM in multidimensional problems is that the effective coefficients preserve the symmetry and the positive-definiteness of the original coefficients, which will be proved next in detail.

\newsec{Preservation of Symmetry}

First, observe that the effective coefficients (\ref{coefmult}) can be rewritten as
\begin{equation}\label{simet0}
\widehat{a}_{jp}=\left\langle a_{kl}(y)\frac{\partial M_j}{\partial y_k}\frac{\partial M_p}{\partial y_l}\right\rangle-\left\langle a_{kl}(y)\frac{\partial N_j}{\partial y_k}\frac{\partial M_p}{\partial y_l}\right\rangle,
\end{equation}
where $M_p=N_p+y_p$.
We will show that the first term on the right-hand side of (\ref{simet0}) is symmetric with respect to indexes $j$ and $p$, and that the second term is null. 

First, we show that the second term of (\ref{simet0}) is null. Consider the equation of the local problem (\ref{problocalmult}) rewritten as
\begin{equation*}
\frac{\partial}{\partial y_j}\left(a_{jl}(y)\frac{\partial M_p}{\partial y_l}\right)=0,\;y\in Y.
\end{equation*}
Let $\phi=\phi(y)$ be $Y$-periodic and continuously differentiable. Then,
\begin{equation}\label{simet1}
\frac{\partial}{\partial y_j}\left(a_{jl}(y)\frac{\partial M_p}{\partial y_l}\right)\phi(y)=0,\;y\in Y,
\end{equation}
so
\begin{equation*}
\frac{\partial}{\partial y_j}\left(a_{jl}(y)\frac{\partial M_p}{\partial y_l}\phi(y)\right) = a_{jl}(y)\frac{\partial M_p}{\partial y_l}\frac{\partial \phi}{\partial y_j}.
\end{equation*}
By applying the average operator, we get
\begin{equation}\label{simet3}
\left\langle\frac{\partial}{\partial y_j}\left(a_{jl}(y)\frac{\partial M_p}{\partial y_l}\phi(y)\right)\right\rangle=\left\langle a_{jl}(y)\frac{\partial M_p}{\partial y_l}\frac{\partial \phi}{\partial y_j}\right\rangle.
\end{equation}
Note that $M_p$ is not $Y$-periodic. However, $\partial M_p/\partial y_l=\partial N_p/\partial y_l+\delta_{pl}$ is $Y$-periodic. So, $\phi(y)a_{jl}(y)\partial M_p/\partial y_l$ is $Y$-periodic, which implies that the term on the left-hand side of (\ref{simet3}) is null. Therefore, we have
\begin{equation*}
\left\langle a_{jl}(y)\frac{\partial M_p}{\partial y_l}\frac{\partial
\phi}{\partial y_j}\right\rangle=0,
\end{equation*}
for every $Y$-periodic $\phi\in C^1(Y)$. In particular, by taking $\phi=N_k$, it follows that
\begin{equation}\label{simet2}
\left\langle a_{jl}(y)\frac{\partial M_p}{\partial y_l}\frac{\partial
N_k}{\partial y_j}\right\rangle=0.
\end{equation}
Thus, by putting (\ref{simet2}) into (\ref{simet0}), we have
\begin{equation}\label{formredcoef}
\widehat{a}_{jp}=\left\langle a_{kl}(y)\frac{\partial M_j}{\partial y_k}\frac{\partial M_p}{\partial y_l}\right\rangle.
\end{equation}
Finally, the symmetry of $\widehat{a}_{jp}$ follows from interchanging the indexes $k$ and $l$, and considering the symmetry of $a_{lk}$:
\begin{equation*}
\widehat{a}_{jp}=\left\langle a_{lk}(y)\frac{\partial M_j}{\partial y_l}\frac{\partial M_p}{\partial y_k}\right\rangle =\left\langle a_{lk}(y)\frac{\partial M_p}{\partial y_k}\frac{\partial M_j}{\partial y_l}\right\rangle=\left\langle a_{kl}(y)\frac{\partial M_p}{\partial y_k}\frac{\partial M_j}{\partial y_l}\right\rangle=\widehat{a}_{pj}.
\end{equation*}
Therefore, $\widehat{a}_{jp}=\widehat{a}_{pj}$, which provides a way to control calculations of the effective coefficients performed via other analytical or numerical techniques.

\newsec{Preservation of Positive-Definiteness}

Preservation of positive-definiteness of the effective coefficients implies that the equation of the homogenized problem is elliptic. Recall from the formulation of the original problem (\ref{probmult}) that the coefficient $a_{jl}$ is positive definite, that is, $\exists c>0,\;\forall\eta\in\mathbb{R}^d\;|\;a_{jl}(y)\eta_j\eta_l\geq c\eta_l\eta_l,\;\forall y\in Y$. In order to show the positive-definiteness of $\widehat{a}_{jp}$, a similar relation must be proved. Indeed, it follows from (\ref{formredcoef}) that
\begin{equation}\label{posit}
\widehat{a}_{jp}\eta_j\eta_p=\left\langle a_{kl}(y)\frac{\partial M_p}{\partial y_k}\frac{\partial M_j}{\partial y_l}\right\rangle\eta_j\eta_p =\left\langle a_{kl}\eta_p\frac{\partial M_p}{\partial y_k}\eta_j\frac{\partial M_j}{\partial y_l}\right\rangle.
\end{equation}
As $a_{kl}$ is positive-definite, we have
\begin{equation}\label{posit3}
a_{kl}(y)\eta_p\frac{\partial M_p}{\partial y_k}\eta_j\frac{\partial M_j}{\partial y_l}\geq c \sum_{l=1}^d \left(\eta_j\frac{\partial M_j}{\partial y_l}\right)^2\mbox{, }\forall y\in Y.
\end{equation}
Then, by substitution (\ref{posit3}) into (\ref{posit}), we have
\begin{equation*}\label{posit1}
\widehat{a}_{jp}\eta_j\eta_p\geq c \sum_{l=1}^d\left\langle\left(\eta_j\frac{\partial M_j}{\partial y_l}\right)^2\right\rangle.
\end{equation*}
From the Cauchy-Buniakovski inequality \cite{Kolmogorov1975}, we have %Cauchy-Schwarz na versão em inglês
\begin{equation*}
\widehat{a}_{jp}\eta_j\eta_p\geq c\sum_{l=1}^d\left\langle\eta_j\frac{\partial M_j}{\partial y_l}\right\rangle^2,
\end{equation*}
which can be rewritten as
\begin{equation}\label{posit2}
\widehat{a}_{jp}\eta_j\eta_p\geq c\sum_{l=1}^d\left[\left\langle\eta_j\frac{\partial N_j}{\partial y_l}\right\rangle+\left\langle\eta_j\frac{\partial y_j}{\partial y_l}\right\rangle\right]^2,
\end{equation}
where
\begin{eqnarray*}
\left\langle\eta_j\frac{\partial N_j}{\partial y_l}\right\rangle=0&\mbox{ and }&\left\langle\eta_j\frac{\partial y_j}{\partial y_l}\right\rangle=\eta_j\langle\delta_{jl}\rangle=\eta_l.
\end{eqnarray*}
So, it follows from (\ref{posit2}) that
\begin{equation*}
\widehat{a}_{jp}\eta_j\eta_p\geq c\sum_{l=1}^d \eta_{l}^2=c\eta_l\eta_l,
\end{equation*}
that is, the effective coefficient $\widehat{a}_{jp}$ is positive-definite, which provides another way to control calculations of $\widehat{a}_{jp}$ performed via other analytical or numerical techniques.

\section{Example}
%Let $d=3$. Consider a thermally-isotropic micro-heterogeneous and $\varepsilon Y$-periodic solid, which occupies the region $\overline{\Omega}=\Omega\cup\partial\Omega=[0,1]^3$. 
Consider problem (\ref{probmult}) with the source term $f(x)=-1$ and isotropic coefficients $a_{jl}^\varepsilon(x)=a^\varepsilon(x)\delta_{jl}$, where $a^\varepsilon(x)=1+0.25\sin(2\pi x_1/\varepsilon)\sin(2\pi x_2/\varepsilon)$. Figure \ref{coefrapmult} shows the behavior of function $a^\varepsilon(x)$, for different values of $\varepsilon$:
\begin{figure}[H]
\centering
%\epsfig{file=coefPB,width=\textwidth}
\includegraphics[width=\textwidth]{coefPB}
\caption{Continuously differentiable, positive, bounded and rapidly oscillating coefficients for the three-dimensional elliptic linear problem.}\label{coefrapmult}
\end{figure}

For this problem, we solve the local problem (\ref{problocalmult}), obtain the effective coefficients (\ref{coefmult}) and solve the homogenized problem (\ref{probhomomult}). Notice that the
first of these problems is two-dimensional from the point of view of the homogenization, as $a_{jl}(y)$ depend only on $y_1$ and $y_2$ because $a^\varepsilon(x)$ depend only on $x_1$ and $x_2$. Therefore, from (\ref{problocalmult}), we solve the following local problem:
\begin{equation}\label{problocalmultexemplo}
 \frac{\partial}{\partial y_1}\left(a(y)\frac{\partial N_p}{\partial y_1}\right)+\frac{\partial}{\partial y_2}\left(a(y)\frac{\partial N_p}{\partial y_2}\right)=-\frac{\partial a}{\partial y_p}, \mbox{  }y\in Y, \quad \langle N_p(y)\rangle=0,
\end{equation}
for each $p=1,\;2,\;3$. The equation in (\ref{problocalmultexemplo}) is a Poisson equation's for $p=1,\;2$ and a Laplace equation's for $p=3$. To solve it, the finite difference method was used  \cite{Samarskii2007,Burden2010}, which comprises the following steps:
\begin{itemize}
\item A mesh of $(N+2)^2$ nodes for the periodic cell $Y$ in the $y_1y_2$-plane is defined. In particular, the local functions $N_p$ are taken to be null on the nodes of $\partial Y$.
\item The discretization of the equation in (\ref{problocalmultexemplo}) on an arbitrary point $(y_{1_i},y_{2_j})$ is
\begin{equation*}
c_{i,j}N^p_{i+1,j}+d_{i,j}N^p_{i-1,j}+e_{i,j}N^p_{i,j+1}+f_{i,j}N^p_{i,j-1}+g_{i,j}N^p_{i,j}=F^p_{i,j},
\end{equation*}
where
\begin{equation*}
c_{i,j}=\frac{a^{(1)}_{i,j}}{2h}+\frac{a_{i,j}}{h^2},\; d_{i,j}=\frac{a_{i,j}}{h^2}-\frac{a^{(1)}_{i,j}}{2h},\; e_{i,j}=\frac{a^{(2)}_{i,j}}{2h}+\frac{a_{i,j}}{h^2},\;  f_{i,j}=\frac{a_{i,j}}{h^2}-\frac{a^{(2)}_{i,j}}{2h},
\end{equation*}
\begin{equation*}
g_{i,j}=-\frac{a_{i,j}}{h^2},\; F^1_{i,j}=-a^{(1)}_{i,j},\; F^2_{i,j}=-a^{(2)}_{i,j},\; F^3_{i,j}=0,
\end{equation*}
\begin{equation*}
a^{(1)}_{i,j}=\frac{\partial a}{\partial y_1}(y_{1_i},y_{2_j}),\; \displaystyle a^{(2)}_{i,j}=\frac{\partial a}{\partial y_2}(y_{1_i},y_{2_j}),
\end{equation*}
and $h=1/(N+1)$ is the spacing between the mesh points in $y_1$ and $y_2$, being $N$ an even number that define the number of points inside of the region $Y$ with $y_3$ fixed.

\item A system of linear equations $\mathbf{A}\mathbf{b}=\mathbf{F}^p$ is assembled as follows. The $N^2\times N^2$ matrix $\mathbf{A}$ is structured in blocks as:
\begin{equation}\label{matrizcoefDF}
\mathbf{A}=\left[\begin{array}{llllll}
S_1&R_1&0&\cdots&\cdots&0\\
T_2&S_2&R_2&\ddots&\ddots&\vdots\\
0&T_3&S_3&R_3&\ddots&\vdots\\
\vdots&\ddots&\ddots&\ddots&\ddots&0\\
\vdots&\ddots&\ddots&T_{N-1}&S_{N-1}&R_{N-1}\\
0&\cdots&\cdots&0&T_N&S_N
\end{array}\right]_{N^2\times N^2},
\end{equation}
where $\{S_j\}_{1\leq j\leq N}$, $\{R_j\}_{1\leq j\leq N-1}$ and $\{T_j\}_{2\leq j \leq N}$, are $N\times N$ blocks given by:
\begin{equation*}
S_{j}=\left[\begin{array}{llllll}
c_{1,j}&d_{1,j}&0&\cdots&\cdots&0\\
e_{2,j}&c_{2,j}&d_{2,j}&\ddots&\ddots&\vdots\\
0&e_{3,j}&c_{3,j}&d_{3,j}&\ddots&\vdots\\
\vdots&\ddots&\ddots&\ddots&\ddots&0\\
\vdots&\ddots&\ddots&e_{N-1,j}&c_{N-1,j}&d_{N-1,j}\\
0&\cdots&\cdots&0&e_{N,j}&c_{N,j}
\end{array}\right]_{N\times N},
\end{equation*}
\begin{equation*}
R_{j}=\left[\begin{array}{llllll}
f_{1,j}&0&0&\cdots&\cdots&0\\
0&f_{2,j}&0&\ddots&\ddots&\vdots\\
0&0&f_{3,j}&0&\ddots&\vdots\\
\vdots&\ddots&\ddots&\ddots&\ddots&0\\
\vdots&\ddots&\ddots&0&f_{N-1,j}&0\\
0&\cdots&\cdots&0&0&f_{N,j}
\end{array}\right]_{N\times N},
\end{equation*}
and
\begin{equation*}
T_{j}=\left[\begin{array}{llllll}
g_{1,j}&0&0&\cdots&\cdots&0\\
0&g_{2,j}&0&\ddots&\ddots&\vdots\\
0&0&g_{3,j}&0&\ddots&\vdots\\
\vdots&\ddots&\ddots&\ddots&\ddots&0\\
\vdots&\ddots&\ddots&0&g_{N-1,j}&0\\
0&\cdots&\cdots&0&0&g_{N,j}
\end{array}\right]_{N\times N}.
\end{equation*}
The vector $\mathbf{b}$ of unknowns in the linear system is defined as:
\begin{equation}\label{vetorsolDF}
\mathbf{b}=\left[w_1\;\;w_2\;...\;w_{N^2}\right]_{N^2\times 1}^T,
\end{equation}
where $w_k$, $k = 1,\ldots,N^2$ are defined as $N^p_{i,j}$ for $i,j = 1,\ldots,N$. Finally, vector $\mathbf{F}^p$ is defined as:
\begin{equation}\label{vetorfonteDF}
\mathbf{F}^p=\left[F^p_{1,1}\;F^p_{1,2}\;...\;F^p_{1,N}\;F^p_{2,1}\;F^p_{2,2}\;...\;F^p_{2,N}\;...\;F^p_{N,1}\;F^p_{N,2}\;...\;F^p_{N,N}\right]_{N^2\times 1}^T.
\end{equation}
%
%Therefore, we get a linear equation system of the form
%\begin{equation*}
%\mathbf{A}\mathbf{b}=\mathbf{F}^p,
%\end{equation*}
%which determine values of the vector (\ref{vetorsolDF}) from coefficient matrix  (\ref{matrizcoefDF}) and of the vector (\ref{vetorfonteDF}).
\end{itemize}

The following results were obtained for a mesh of 10000 nodes, for each $p$. The full time of this simulation was approximately 100 seconds. Figure \ref{local2D} shows the solutions $N_p$ of the local problems, for $p=1,2$, as $N_3(y)=0$.
\begin{figure}[H]
\centering
%\epsfig{file=LocaisPB,width=\textwidth}
\includegraphics[width=\textwidth]{LocaisPB}
\caption{Solutions of the local problems $N_1$ and $N_2$, being
$N_3=0$.}\label{local2D}
\end{figure}

Notice that the local solutions $N_p$ were calculated assuming homogeneous Dirichlet conditions on the boundary of the periodic cell, that is, the condition $\langle N_p \rangle=0$ in local problem (\ref{problocalmult}) was replaced by $N_p|_{\partial Y}=0$, in order to ensure that the first-order f.a.s $u^{(1)}$ satisfy exactly the boundary conditions of the original problem. On the other hand, notice that Lemma \ref{lemageral} uses the null-average condition for uniqueness. In general, different choices of uniqueness conditions on the local solutions will affect the behavior of $u^{(1)}$ only in the neighborhood of the boundary. However, the difference between the local solutions obtained from these different conditions is an additive constant and, as the effective coefficients
depend only on the derivatives of the local solutions, such a difference does not affect its values. So, in order to assess the consequence of using such an alternative uniqueness condition, we calculated the averages of the local solutions obtained with homogeneous conditions, which were found to be very close to zero. Specifically, $\langle N_1\rangle=2.365\times10^{-5}$,
$\langle N_2\rangle=-2.309\times10^{-5}$ e $\langle N_3\rangle=N_3=0$. Furthermore, these values tend to zero as the step of the discretization tends to zero, $h\rightarrow0$, which means that the exact local solutions of this example satisfy both uniqueness conditions.

Now, in order to calculate the effective coefficients $\widehat{a}_{jl}$, we start by simplifying (\ref{coefmult}), as it requires to obtain the derivative of the solutions $N_p$, which would require extra steps in programming. In addition, the numerical calculation of the derivatives would increase the simulation time with the consequent increase of computational cost and numerical error. So, in order to simplify (\ref{coefmult}), we rewrite the term that contains the derivative of $N_p$ as
\begin{equation}\label{relGreen}
a_{jl}(y)\frac{\partial N_p}{\partial y_l}=\frac{\partial}{\partial y_l}\left(N_p(y) a_{jl}(y)\right)-N_p(y)\frac{\partial a_{jl}}{\partial y_l}.
\end{equation}
By substituting (\ref{relGreen}) into (\ref{coefmult}), we get
\begin{equation*}
\widehat{a}_{jp}=\left\langle a_{jp}(y)\right\rangle-\left\langle N_p(y)\frac{\partial a_{jl}}{\partial y_l}\right\rangle,
\end{equation*}
as $N_p(y)$ and $a_{jl}(y)$ are $Y$-periodic. As the coefficient in this example is isotropic, we have
\begin{equation*}
\widehat{a}_{jp}=\left\langle a(y)\right\rangle\delta_{jp}-\left\langle N_p(y)\frac{\partial a}{\partial y_l}\right\rangle\delta_{jl},
\end{equation*}
that is,
\begin{equation*}
(\widehat{a}_{jp})_{1\leq j,p\leq3}=\left[
\begin{array}{ccc}
\displaystyle\langle a(y)\rangle-\left\langle N_1(y)\frac{\partial a}{\partial y_1}\right\rangle&\displaystyle-\left\langle N_2(y)\frac{\partial a}{\partial y_1}\right\rangle&0\\
\displaystyle-\left\langle N_1(y)\frac{\partial a}{\partial y_2}\right\rangle&\displaystyle\langle a(y)\rangle-\left\langle N_2(y)\frac{\partial a}{\partial y_2}\right\rangle&0\\
0&0&\langle a(y)\rangle
\end{array}
\right].
\end{equation*}
We calculated each element of the effective coefficient matrix by the Simpson's Method \cite{Burden2010}, which yields
\begin{equation*}
(\widehat{a}_{jp})_{1\leq j,p \leq3}=
\left[\begin{array}{rrr} 0.994587&0.000103&0\\
 -0.000084&0.994607&0\\
0&0&1
\end{array}\right].
\end{equation*}

Notice that this matrix is positive-definite, and that the elements off the main diagonal are very close to zero, that is, we can to assume that $\widehat{a}_{jp}$ satisfy the symmetry condition. Indeed, finer meshes would produce values closer to $1$ in the main diagonal and to $0$ outside it, respectively, but the computational cost would be also much higher, which might compromise convergence of the numerical method. Thus, when the step of the discretization tends to zero, $h\rightarrow0$, we have that $\widehat{a}_{jp}=\delta_{jp}$, which is, obviously, symmetric and positive definite, as expected. Also, this result can be controlled via variational bounds for the effective coefficients \cite{Bakhvalov1989}
\begin{equation*}\label{cotascoefsmult}
\langle a^{-1}(y)\rangle^{-1}\delta_{jp}\leq \widehat{a}_{jp} \leq \langle
a(y)\rangle\delta_{jp},
\end{equation*}
where $\langle a(y)\rangle=1$ and $\langle a^{-1}(y)\rangle^{-1}=0.984055$.

Using the fact that $\widehat{a}_{jp}=\delta_{jp}$ as $h\rightarrow0$, the homogenized problem (\ref{probhomomult}) is given by the following Poisson equation with homogeneous Dirichlet conditions:
\begin{equation}\label{probhomomult2}
\frac{\partial^2 u_0}{\partial x_1^2}+\frac{\partial^2 u_0}{\partial x_2^2}+\frac{\partial^2 u_0}{\partial x_3^2} = -1, \mbox{ } x\in\Omega \quad u_0 = 0, \mbox{ } x\in\partial\Omega.
\end{equation}
We solve problem (\ref{probhomomult2}) via separation of variables as in \cite{Weinberger1995}, which leads to a Fourier series with orthogonal basis $\{\sin(n\pi x_i)\}_{(n,i)\in\mathbb{N}\times\{1,2,3\}}$. Then, the solution $u_0$ is
\begin{equation*}
u_0(x)=\frac{64}{\pi^5}\sum_{i,j,k=1}^\infty \frac{\sin(i\pi x_1)\sin(j\pi x_2)\sin(k\pi x_3)}{\left((2i-1)^2+(2j-1)^2+(2k-1)^2\right)(2i-1)(2j-1)(2k-1)}.
\end{equation*}

Figures \ref{surfsol} and \ref{erros} show the first-order f.a.s $u^{(1)}$ and the homogenized solution $u_0$, and the absolute differences $|u^{(1)}-u_0|$ for
$\varepsilon=0.5$, $0.25$ and $0.125$, respectively:
\begin{figure}[H]
\centering
%\epsfig{file=solPB,width=\textwidth}
\includegraphics[width=\textwidth]{solPB}
\caption{First-order f.a.s $u^{(1)}$ (left) and homogenized solution $u_0$ (right).}
\label{surfsol}
\end{figure}

\begin{figure}[H]
\centering
%\epsfig{file=erroPB,width=\textwidth}
\includegraphics[width=\textwidth]{erroPB}
\caption{Absolute differences between the first-order f.a.s $u^{(1)}$ and the homogenized solution $u_0$ for $\varepsilon=1/2,1/4,1/8$.}
\label{erros}
\end{figure}

In the three-dimensional visualization corresponding to Figure \ref{surfsol}, it is not possible to distinguish the slight differences between the first-order f.a.s $u^{(1)}$ and the homogenized solution $u_0$. However, such differences are observable in the two-dimensional representation. Specifically, the left plot in Figure \ref{surfsol} shows a slight asymmetry, whereas the right plot is perfectly symmetrical. However, both solutions are very close, as shown in Figure \ref{erros}. 

\newsec{Conclusions}

The example presented in this work illustrates the fact that the AHM is an efficient alternative approach to the direct resolution of problems with rapidly oscillating coefficients. Indeed, a direct approach via a traditional method, such as the finite difference method used here, would require an extremely fine three-dimensional discretization of the domain in order to capture the rapid variation of the coefficients. This implies a high computational cost and compromises the convergence of the numerical method, which in fact happened when we attempted to solve the original problem of this example by the finite difference method. Thus, AHM provided a useful alternative to obtain an approximate solution of this problem in the form of a formal asymptotic solution, which was proved to be very close to the exact solution. Moreover, from a practical point of view, the problems generated from the application of AHM, that is, the homogenized and local problems, are much easier to  solve. In this case, local problems are two-dimensional and their coefficients, although non-constant, they do not oscillate rapidly, which allowed the finite difference method to be successfully applied. Also, as the homogenized problem has constant coefficients, it was solved analytically via the Fourier method's. On the other hand, the search of the exact solution of the original problem would require non-traditional approaches, such as multiscale numerical methods.The example is important for those who will go study the AHM, so know the practical process after applying the method, such: sequence of solve the problems decoupled in the process of homogenization, possible analytical and numerical methods to be used, how to control numerical results for the effective property. This example is a particular case (isotropic property, constant source, homogeneous Dirichlet boundary conditions and continuous coefficients) of much more general and complex problems, but it illustrates essentially the sequence of resolutions to obtain an exact solution approximation, indeed, an approximation of $\mathcal{O}(\sqrt{\varepsilon})$.

%\begin{thebibliography}{8}

%\bibitem{Bakhvalov1989}
%N.S. Bakhvalov, G.P. Panasenko, ``Homogenisation: Averaging Processes in Periodic Media'', Kluwer Academic Publishers, Dordrecht, 1989.
%
%\bibitem{Burden2008}
%R.L. Burden, J.D. Faires, ``Análise Numérica'', Cengage Learning, São Paulo, 2008.
%
%\bibitem{Cioranescu}
%D. Cioranescu, P. Donato, ``An Introduction to Homogenization'', Oxford University Press, New York, 1999.
%
%\bibitem{Cunha2000}
%C. Cunha, ``Métodos Numéricos'', Editora da Unicamp, São Paulo, 2000.
%
%\bibitem{Kolmogorov1975}
%A.N. Kolmogorov, S.V. Fomin, ``Elementos de la Teoria de Funciones y del Análisis Funcional'', Mir, Moscou, 1975.
%
%\bibitem{Torquato2002}
%S. Torquato, ``Random Heterogeneous Materials: Microstructure and Macroscopic Properties'', Springer-Verlag, New York, 2002.
%
%\bibitem{Weinberger1995}
%H.F. Weinberger, ``A First Course in Partial Differential Equations with Complex Variables and Transform Methods'', Dover, New York, 1995.
%
%\bibitem{Choe2003}
%H.J. Choe, K.-B. Kong, C.-O. Lee, ``Convergence in $L^p$ space for the homogenization problems of elliptic and parabolic equations in the plane'', J. Math. Anal. Appl., 287, 321-336, 2003.

%\end{thebibliography}

%\appendix
%
%\section{Voigt-Reuss's bounds}
%
%In the variational setting, the effective behavior of medium micro-heterogeneous is characterized by effective energy density $\widehat{\epsilon}(\overline{\sigma})$ given by the minimum energy principle
%\begin{equation}\label{PME1}
%\widehat{W}(\overline{\epsilon})=\inf_{\epsilon\in S_W(\overline{\epsilon})}\langle W(\epsilon)\rangle,
%\end{equation}
%where $S_W(\overline{\epsilon})$ is the admissibility set in $Y$-periodic compatible gradiente-type fields $\epsilon$ with components in $L^2(Y)$ and mean value $\langle \epsilon \rangle= \overline{\epsilon}$, and $W(\epsilon)$ is the convex constitutive energy. The effective complementary energy density $\widehat{U}(\overline{\sigma})=\widehat{W}^*(\overline{\sigma})$ is given by the minimum complementary energy principle
%\begin{equation}\label{PME2}
%\widehat{U}(\overline{\sigma})=\inf_{\sigma\in S_U(\overline{\sigma})}\langle U(\sigma)\rangle,
%\end{equation}
%where $S_U(\overline{\sigma})$ is the admissibility set of 
%$Y$-periodic divergence-free flux-type fields $\sigma$ with components in $L^2(Y)$ and mean value $\langle \sigma \rangle= \overline{\sigma}$, and $U(\sigma)=W^*(\sigma)$ with the asterisk denoting the Legendre transform, that is,
%\begin{displaymath}
%W^*(\sigma)=\sup_{\epsilon}\{\sigma\cdot\epsilon -W(\epsilon)\},
%\end{displaymath}
%where the dot stands for a suitable inner product. In this context, the constitutive relation and the effective law are given by
%\begin{displaymath}
%\sigma=\frac{\partial W}{\partial \epsilon}(\epsilon)
%\end{displaymath}
%and
%\begin{displaymath}
%\overline{\sigma}=\frac{\partial \widehat{W}}{\partial \overline{\epsilon}}(\overline{\epsilon}),
%\end{displaymath}
%respectively. Convexity of $W(\epsilon)$ guarantees that $\widehat{U}^*(\overline{\sigma})=\widehat{W}(\overline{\epsilon})$, so it suffices to take constant fields $\sigma=\overline{\sigma}$ and $\epsilon=\overline{\epsilon}$ in both minimum principles (\ref{PME1}) and (\ref{PME2}) to obtain the elementary bounds
%\begin{displaymath}
%\langle U \rangle^*(\overline{\epsilon})\leq \widehat{W}(\epsilon)\leq \langle W \rangle (\overline{\epsilon}).
%\end{displaymath}
%By linearity of behavior, the lower bound $\langle U \rangle^*(\overline{\epsilon})$ reduces to bound $\langle W^{-1}\rangle^{-1}(\overline{\epsilon})$. Thus, in this case, for a linear behavior, the constitutive and effective energy functions $W(x,\epsilon)=(1/2)a_{ij}(x)\epsilon_i\epsilon_j$ and $\widehat{W}(\overline{\epsilon})=(1/2)\widehat{a}_{ij}\overline{\epsilon}_i\overline{\epsilon}_j$ result
%\begin{displaymath}
%\frac{1}{2}\langle a^{-1}_{ij}(x) \rangle^{-1}\overline{\epsilon}_i\overline{\epsilon}_j\leq \frac{1}{2}\widehat{a}_{ij}\overline{\epsilon}_i\overline{\epsilon}_j\leq \frac{1}{2}\langle a_{ij}(x)\rangle \overline{\epsilon}_i\overline{\epsilon}_j,
%\end{displaymath}
%from where the classical limits upper Voigt \cite{Voigt1889} and lower Reuss \cite{Reuss1929} are obtained for the effective property.
%
%\section{Generalized Maximum Principle for Elliptic Equations}
%
%Let $u\in H_0^1(\Omega)$, $\Omega\subset\mathbb{R}^d$, generalized solution of the problem
%\begin{equation*}
%\left\{\begin{array}{l}
%\displaystyle \frac{\partial}{\partial x_i}\left(a_{ij}(x)\frac{\partial u}{\partial x_j}\right)+\frac{\partial}{\partial x_i}\left(a_{i}(x)u\right)+b_{i}(x)\frac{\partial u}{\partial
%x_i}+a(x)u=f(x)+\frac{\partial f_i(x)}{\partial
%x_i},\;\;x\in\Omega \setminus \Gamma,\\%\bigcup_{i=1}^N \partial\Omega_i
%u|_{\partial\Omega}=0,
%\end{array}\right.
%\end{equation*}
%where $\Gamma=\cup_{i=1}^N \partial \Omega_i$, and $\partial \Omega_i$ are functions that split $\Omega$ into a number finite of disjoint domains $\widetilde{\Omega}_q$, and the coefficients $a_{ij}$, $a_i$, $b_i$, $a$, as well as $f$ and $f_i$ are infinitely differentiable in each of the $\overline{\widetilde{\Omega}}_q$.
%The coefficients $a_{ij}$ satisfy symmetry condition, $a_{ij}(x)=a_{ji}(x)$, and positivity $a_{ij}(x)\eta_i\eta_j\geq\kappa \eta_i\eta_i,\;\;\forall\eta\in\mathbb{R}^d$,
%where $\kappa$ is a positive constant.
%
%For this solution the following estimate is valid:
%\begin{equation}
%\|u\|_{H_0^1(\Omega)}\leq c\left(\|f\|_{L^2(\Omega)}+\sum_{i=1}^s
%\|f_i\|_{L^2(\Omega)}\right),
%\end{equation}
%where $c>0$ is a constant \cite{Bakhvalov1989}.

\bibliographystyle{ieeetr}
\bibliography{TEMA}

\end{document}
\newpage
$ \  \  $  \thispagestyle{myheadings}  \markboth{      }{   }