1 Introduction
Many engineering and scientific structures have important microscale variations, such as lattice materials [Reference Barchiesi, Eugster, Dell’Isola and Hild6, Reference Somnic and Jo62, Reference Somnic and Jo63], the windings in electrical machinery [Reference Romanazzi, Bruna and Howey58], electromagnetism in micro-structured material [Reference Craster, Mityushev and Ruzhansky20, Reference Niyonzima44], food drying [Reference Welsh, Simpson, Khan and Karim70] and slow Stokes flow through porous media [Reference Brown, Popov and Efendiev11]. The engineering challenge is to predict and understand the dynamics on a system scale significantly larger than the microscale details. Homogenization, via the recognition of multiple physical length scales, is the common approach [Reference Abdulle and Pouchon1, Reference Cooper17, Reference Engquist and Souganidis26, Reference Gustafsson and Mossino30]. This article further develops an alternative approach to homogenized modelling of the large scale dynamics via a newly established general rigorous theory [Reference Roberts52, Section 2.5]. This new complementary approach provides a route to systematic refinements of the homogenization for a finite separation of scales and to novel quantification of the remainder error. The significance of this work is not that it effectively reproduces some classic homogenization results, but that it invokes a powerful dynamical systems framework [Reference Bunder and Roberts12, Reference Chen, Roberts and Bunder16, Reference Cox and Roberts19, Reference Hochs and Roberts32, Reference Muncaster42, Reference Watt and Roberts69] which provides complementary new insights and powerful practical extensions of previously known aspects of homogenization. Consequently, a background in modern dynamical systems [Reference Strogatz64, Reference Teschl65] is useful to appreciate the direction developed here for homogenization.
Consider materials with complicated microstructure: we want to model their large-scale dynamics by equations with effective ‘average’ coefficients. Heterogeneous diffusion in one-dimensional (1D) space is the simplest such class of systems: suppose the material, in spatial domain
$[0,L]$
, has structure so that ‘material’ or ‘heat’
$u(t,x)$
diffuses according to
where the heterogeneous diffusion coefficient
$\kappa (x)$
has some complicated multiscale structure. For definiteness, let us initially suppose the macroscale boundary conditions are simply that u is L-periodic in space x (the end of Section 2 discusses other boundary conditions). Our challenge is to derive, support and characterize the macroscale model PDE
for some effective mean field
$U(x,t)$
, some effective macroscale diffusivity K, with some higher order corrections, and possibly an error estimate—both the latter represented by the ellipsis “
$\ldots $
” in (1.2). The problem we address is not to find approximate solutions to PDE (1.1), but instead to provide new justification, construction and interpretation of the simpler PDE (1.2), and cognate higher-order PDEs, as purely macroscale models of the multiscale PDE (1.1). I also comment throughout on the closely connected problem [Reference Abdulle and Pouchon1, Reference Allaire, Lamacz-Keymling and Rauch3, Reference Cornaggia and Guzina18, Reference Duerinckx, Gloria and Shirley25, Reference Shahraki and Guzina61] of wave propagation through heterogeneous material,
and this corresponding effective homogenized PDE model that predicts its macroscale dispersive waves. The homogenization approach for both classes are much the same, so the results here apply to elastic materials and other waves, with some caveats. However, the treatment and theory is clearer and more rigorous for the class of diffusive systems (1.1).
Here we focus on cases where the heterogeneity
$\kappa (x)$
is quasi-periodic in space. For an example application, Sections 4 and 5 focus on homogenization for the family of heterogeneous diffusivities
$\kappa =1/[1+a_1\cos (2\pi x/\ell _1)+a_2\cos (2\pi x/\ell _2)]$
for arbitrary microscale periods
$\ell _1,\ell _2$
with
$\ell _1>\ell _2$
. Sections 2 and 3 develop theory and supporting characteristics of homogenization in the general case where the heterogeneity
$\kappa $
is expressible as a function of two variables, say
$\kappa (x_1,x_2)$
, in which
$\kappa $
is
$\ell _i$
-periodic in variable
$x_i$
. Then the quasi-periodic heterogeneity in (1.1) is
$\kappa (x,x)$
. I emphasize that this work does not address the PDE
$u_t=\partial _x\{\kappa (x/\epsilon ,x/\epsilon )u_x\}$
: there is no
$\epsilon $
in the analysis here; and consequently, and in contrast to extant well established theory of quasi-periodic homogenization (for example, see [Reference Jikov, Kozlov, Oleinik and Yosifian33], [Reference Allaire2, Section 7.4]), the framework here does not invoke the infinite length scale separation limit “
$\epsilon \to 0$
”. Instead, the results herein apply for finite scale separation (quantified by Section 5), and form arbitrarily high-order homogenizations, not just leading order [Reference Jikov, Kozlov, Oleinik and Yosifian33] nor just second-order [Reference Rodrigues Lopes and Andrade Pires57].
Clearly, the approach developed here straightforwardly extends to an arbitrary number of incommensurable periods in the structure of the diffusivity
$\kappa $
, but for simplicity of introducing the key new ideas, techniques and results, we address here the case of just two periods in
$\kappa $
.
Most extant conventional approaches to practically computing a homogenization depend upon the identification of a representative volume element (RVE) [Reference Somnic and Jo62], or a representative cell. Even in the simple case of periodic media, an uncertain issue is an appropriate size of an RVE [Reference Klarmann, Gruttmann and Klinkel34, Reference Matous, Geers, Kouznetsova and Gillman38, Reference Nguyen, Stroeven and Sluys43, Reference Saeb, Steinmann and Javili59, Reference Schneider60, Reference Somnic and Jo63, p. 2, Section 5.1, Section 2.3, Section 1.2.2, p. 2, p. 13, respectively]. A further issue with such RVE-based arguments is that for quasi-periodic heterogeneities, there is no one RVE [Reference Anthoine5, Section 2.1]. When the ratio between microscale lengths
$\ell _1,\ell _2$
is irrational, then all volume elements, no matter what size is chosen, are different: no volume element is representative. Relatedly, when the microscale lengths
$\ell _1,\ell _2$
have a large least common multiple, then any true RVE must be correspondingly large—resulting in onerous analysis or computation. For example, Blanc and Le Bris [Reference Blanc and Le Bris9] comment that for quasi-periodic homogenization, “cell problems must be in general set on the whole space”: herein, we show that this “must” is incorrect for our approach. In contrast to “the whole space”, here, Sections 2 and 3 instead develop a novel rigorous systematic approach that uses an ‘RVE’ cell which is the size of the largest microscale length,
$\max \{\ell _1,\ell _2\}$
, albeit an RVE of more ‘spatial’ dimensions.
I expect straightforward extensions of the approach developed herein to empower cognate analysis and modelling of the following:
-
• cases where the heterogeneity is quasi-periodic in any number of microscale lengths
$\ell _1>\ell _2>\ell _3>\cdots $
; -
• PDEs (1.1) where the field u is itself in a vector space such as the cross-sectional structure of either shear dispersion [Reference Mercer and Roberts40, Reference Mercer and Roberts41], or of elastic beams [Reference Roberts49];
-
• multiple macroscale space dimensions as in the modelling of plates and shells [Reference Roberts and Bunder55];
-
• PDEs like (1.1) but with macroscale variations in the material parameters as well as microscale variations [Reference Anthoine5]—sometimes called functionally graded materials;
-
• PDEs like (1.1) but with nonlinear effects [Reference Matous, Geers, Kouznetsova and Gillman38]—in contrast to the experience of Blanc and Le Bris [Reference Blanc and Le Bris9], in the framework developed here, nonlinear problems do not involve any difficulty other than extra algebraic details.
Further, Section 6 shows that the approach, analysis and results developed herein also include three-scale multiscale systems where notionally
$\ell _2\ll \ell _1\ll L$
, that is, where there is a microscale, a mesoscale and a macroscale [Reference Cruz-González, Ramirez-Torres, Rodriguez-Ramos, Penta and Lebon21, Reference Ramirez-Torres, Penta, Rodriguez-Ramos, Merodio, Sabina, Bravo-Castillero, Guinovart-Diaz, Preziosi and Grillo45].
Other analytic techniques such as mathematical homogenization or the method of multiple scales [Reference Cornaggia and Guzina18, Reference Shahraki and Guzina61] provide a mechanism for deriving effective PDE models such as (1.2). However, such previous underlying theory requires the limit of infinite scale separation between the microscale “cell”-sizes
$\ell _i$
and the macroscale domain size L by invoking
$\epsilon :=\ell _i/L\to 0$
(for example, [Reference Allaire2, Reference Barchiesi, Eugster, Dell’Isola and Hild6, Reference Welsh, Simpson, Khan and Karim70, p. 5, p. 6 and Section 2.4, respectively] and [Reference Altmann, Henning and Peterseim4]). In contrast, and arising from developments in dynamical systems, Sections 2 and 3 respectively use phase-shift embedding and slowly varying theory [Reference Roberts52, Reference Roberts and Bunder55] to provide rigorous support to the modelling at finite scale separation: that is, when the macroscale length scale L of interest is larger than the microscales
$\ell _1$
and
$\ell _2$
, but not infinitely so. In the approach here, not only is there no limit, there is no need for such a defined “
$\epsilon $
”. Such finite scale separation empowers this approach to capture the dispersion in wave homogenization which Abdulle and Pouchon [Reference Abdulle and Pouchon1, p. 3] call for in their discussion stating that “new effective models are required”. The phase-shift embedding developed here further illustrates the general principle that macroscale models of microscale structures are often best phrased as ensemble averages [Reference Birch and Young8, Reference Drugan and Willis23, Reference van Kampen68].
It is well recognized that homogenization is physically effective for reasonably well separated length scales: that the resolved macroscale structures, with length
$\ell _M$
say, satisfy
$\ell _M\gg \ell _1$
. But what does this mean? Much extant homogenization theory requires the mathematical limit
$\ell _1/\ell _M\to 0$
. However, physically, the ratio
$\ell _1/\ell _M$
is always finite. By setting homogenization in a new rigorous framework, Section 5 analyses homogenization for finite
$\ell _1/\ell _M$
to high-order and presents evidence that, in a class of problems, a quantitative limit to the spatial resolution of homogenization is that
$\ell _1/\ell _M\lesssim 1/2$
, equivalently,
$\ell _M\gtrsim 2\ell _1$
. Mathematically, a wider scale separation is not necessary for homogenization.
Of course, practically, we usually prefer a wider scale separation because of many complicating factors in construction and use of practical homogenization models. However, in the sense discussed by Section 5, this quantitative, near-sharp, lower bound on the required scale separation indicates how little scale separation is necessary. Indeed, because it applies at finite scale separation, the analysis here encompasses high-contrast media and so also detects features occurring in this particularly interesting class of problems [Reference Cooper17, p. 2].
The new approach to homogenization that we develop here creates a powerful dynamical systems framework for addressing other modelling issues in the future. For a time-dependent simulation of homogenized dynamics, one must provide some initial conditions for the homogenized PDE (1.2). These initial conditions are surprisingly nontrivial at finite scale separation, but the new framework developed here comes with well-established methods to derive the correct initial conditions for long-term forecasts [Reference Cox and Roberts19, Reference Roberts47, Reference Roberts49, Reference Roberts51]. The projection that provides the correct initial conditions, also provides the correct projection of any applied “forcing” [Reference Roberts47, Section 7], including the cases of control and system uncertainty. System uncertainty is a challenge for other techniques [Reference Altmann, Henning and Peterseim4, Reference Matous, Geers, Kouznetsova and Gillman38, p. 3 and Section 8, respectively]. Further, the homogenized PDE (1.2) needs macroscale boundary conditions—a “highly challenging problem” according to Shahraki and Guzina [Reference Shahraki and Guzina61], and which Cornaggia and Guzina [Reference Cornaggia and Guzina18] called a “complex problem that is still an active research topic”. Such boundary conditions are nontrivial and it is well recognized that consideration of boundary layers must be accounted for. However, again the new framework developed here comes with a rigorous approach to derive correct boundary conditions for (1.2) to empower accurate predictions at finite scale separation [Reference Chen, Roberts, Bunder, Nelson, Hamilton, Jennings and Bunder14–Reference Chen, Roberts and Bunder16, Reference Roberts48, Reference Roberts49]. For an example in a different field but cognate to the macroscale modelling developed here, Mercer and Roberts [Reference Mercer and Roberts41] generalized the advection-diffusion ‘homogenization’ model of shear dispersion in a narrow pipe to arbitrary order and to a quantitative bound on the spatial resolution (its Section 2), with guaranteed accurate initial conditions (its Section 3), and correct macroscale boundary conditions at the pipe’s entry and exit (its Section 4). That is, this dynamical systems approach to homogenization potentially provides a rational complete modelling of not just the PDE at finite scale separation, but also material variations, uncertainty, control, forcing, and initial and boundary conditions. In their review of nonlinear homogenization, Matous et al. [Reference Matous, Geers, Kouznetsova and Gillman38, Section 8] commented that “verification ensures both order of accuracy and consistency … and even mathematical foundations for verification are lacking”—the approach developed here provides the required mathematical foundation.
The new approach developed herein has no need of a variational formulation of the governing equations [Reference Altmann, Henning and Peterseim4, Reference Barchiesi, Eugster, Dell’Isola and Hild6, Section 3 and Section 2.4, respectively], and so applies to a much wider class of problems than many extant homogenization methods.
Throughout this article, in the details of the construction of homogenizations, one may observe that the variables
$x_0,x_1,x_2$
appear in much the same sort of algebraic expressions as traditionally obtained by the method of multiple scales. What this parallelism implies is that in suitable circumstances, the algebraic machinations of the method of multiple scales are indeed sound [Reference Roberts52, Corollary 12]. What is new here is that we place such algebraic details in a new, more rigorous and more extensive framework. The new framework developed here improves the clarity, precision, power and physical interpretability of homogenization.
2 Phase-shift embedding
Let us embed the specific given physical quasi-periodic PDE (1.1) into a family of PDE problems formed by all phase-shifts of the quasi-periodic microscale. This embedding is a novel and rigorous twist to the concept of an RVE.
As indicated by Figure 1, we create the desired embedding by considering a field
${\mathfrak {u}}(t,x_0,x_1,x_2)$
satisfying the PDE
in the domain
${\mathcal D}:=\{(x_0,x_1,x_2) \mid {x_i}\in [0,\ell _i]\},$
for the macroscale
$\ell _0:=L$
, the specified microscale
$\ell _1,\ell _2$
and with
$\ell _i$
-periodic boundary conditions in
${x_i}$
. We assume
$\kappa (x_1,x_2) \geq \kappa _{\min }>0$
. The domain
${\mathcal D}$
in
$x_0x_1x_2$
-space (Figure 1) is often termed cylindrical as it is long and thin. I emphasize that this domain has finite aspect ratio: we do not take any limit involving its aspect ratio tending to zero nor to infinity.
Figure 1 indicates that we consider
$x_0=x$
. The distinction between
$x_0$
and x is that partial derivatives in
$x_0$
are done keeping the other
$x_i$
constant, whereas partial derivatives in x are done keeping phases
$\phi _i$
constant.
Lemma 2.1. For every solution
${\mathfrak {u}}(t,x_0,x_1,x_2)$
of the embedding PDE (2.1) and for every vector of phases 
, the field 
(that is, the field
${\mathfrak {u}}$
evaluated on the solid-blue lines in Figure 1) satisfies the heterogeneous diffusion PDE
Hence, the field
$u(t,x):={\mathfrak {u}}(t,x,x,x)$
satisfies the given heterogeneous PDE (1.1).
Figure 1 Schematic cylindrical domain
${\mathcal D}$
of the multiscale embedding PDE (2.1) for a field
${\mathfrak {u}}(t,x_0,x_1,x_2)$
(here
$L=\ell _0=6$
,
$\ell _1=1.62$
,
$\ell _2=0.72$
). We obtain solutions of the heterogeneous diffusion PDE (1.1), or (2.2), on such blue lines as
$u(t,x)={\mathfrak {u}}(t,x,x+\phi _1,x+\phi _2)$
for every pair of constant phases
$\phi _1$
and
$\phi _2$
(here
$\phi _1=0.82$
and
$\phi _2=0.32$
), and where the last two arguments of
${\mathfrak {u}}$
are modulo
$\ell _1$
and
$\ell _2$
, respectively.
Recall that the most common cell boundary conditions assumed for RVEs are periodic, although in the usual homogenization arguments, other boundary conditions are equally as valid in practice albeit giving slightly different results [Reference Mercer, Mandadapu and Papadopoulos39]. In contrast, here the
$\ell _1,\ell _2$
-periodic cell boundary conditions are not assumed, but arise naturally due to the ensemble of phase-shifts. That is, what previously had to be assumed, here arises naturally.
Proof. Start by considering the left-hand side time derivative in PDE (2.2):
namely the right-hand side of (2.2). Hence, provided PDE (2.1) has boundary conditions of
$\ell _i$
-periodicity in
$x_i$
, every solution of the embedding PDE (2.1) gives a solution of the original PDE (1.1) for every multi-dimensional phase-shift 
of the heterogeneity.
In particular, the field 
satisfies the given heterogeneous PDE (1.1).
Lemma 2.2 (Converse).
Suppose we have a set of solutions 
of the phase-shifted PDE (2.2)—a set parametrized by the phase vector 
—and the set depends smoothly upon 
. Then the field
${\mathfrak {u}}(t,x,x_1,x_2):=u_{(x_1-x,x_2-x)}(t,x)$
satisfies the embedding PDE (2.1).
Proof. First, from the PDE (2.1), consider
Second, since
$\phi _i=x_i-x$
, that is,
$x_i=x+\phi _i$
, then
Third, hence the right-hand-side of PDE (2.1) becomes
the right-hand side of PDE (2.2). Lastly, since 
, it follows that
${\mathfrak {u}}:=u_{(x_1-x,x_2-x)}(t,x)$
satisfies the embedding PDE (2.1).
Corresponding results to Lemmas 2.1 and 2.2 hold for the heterogeneous wave PDE
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
, with 
replaced by 
.
Consequently, PDEs (2.1) and (2.2) are equivalent, and they may provide us with a set of solutions for an ensemble of materials all with the same microscale heterogeneity structure, but with the structural phase of the material shifted through all possibilities. The key difference between PDEs (2.1) and (2.2) is that although PDE (2.2) is heterogeneous in space x (as is (1.1)), the embedding PDE (2.1) is homogeneous in
$x_0$
. Because of this homogeneity, Section 3 is empowered to apply an existing rigorous theory for slow variations in space that leads to the desired homogenized model PDE (1.2).
Discussion on translational symmetry. One should look at a macroscale homogenized model such as PDE (1.2) and wonder from where the continuous x-translation symmetry comes. The symmetry is not in the underlying microscale heterogeneous PDE (1.1). In the case where heterogeneity
$\kappa (x)$
is periodic, the PDE (1.1) has discrete x-translational symmetry, which via the infinite scale separation limit as
$\ell _1\to 0$
, might be viewed as approximately the continuous symmetry in the homogenization. However, in the case of quasi-periodic heterogeneity
$\kappa (x)$
, for incommensurate
$\ell _1$
and
$\ell _2$
, there is no discrete translational symmetry. Instead of such a problematic limit, we here instead create the continuous translational symmetry by considering all together the family of problems formed by all phase shifts of the heterogeneity—the PDEs (2.1) and (2.2). The analysis of the whole family then naturally provides the continuous x-translational symmetry then inherited by the homogenization (1.2).
How does one get predictions correct for a specific heterogeneity? That is, for a given specific phase? The answer is that it is the correct boundary conditions for the homogenization PDE (1.2) that both reflect the specific phase of the heterogeneity at and near a specific boundary, and that also account for boundary layers. Such boundary conditions lie just beyond the scope of this article, but should be derivable from those for the microscale PDE (1.1) by further developing established dynamical systems analysis [Reference Chen, Roberts, Bunder, Nelson, Hamilton, Jennings and Bunder14–Reference Chen, Roberts and Bunder16, Reference Mercer and Roberts41, Reference Roberts48, Reference Roberts49].
3 The slow manifold of any-order homogenization
Let us analyse the embedding PDE (2.1) for a useful slow manifold. This slow manifold expresses and supports the emergence of a precise homogenization of the original heterogeneous PDE (1.1). Since the PDEs herein are linear, the slow manifold is more specifically a slow subspace, but I use the term manifold as the same framework and theory immediately also applies to cognate nonlinear systems.
Rigorous theory [Reference Roberts52], inspired by earlier more formal arguments [Reference Roberts50, Reference Roberts46], establishes how to create a PDE model for the macroscale spatial structure of PDE solutions in cylindrical domains
${\mathcal D}$
like Figure 1. The technique is to base analysis on the case where variations in
$x_0$
are approximately negligible, and then treat slow, macroscale, variations in
$x_0$
as a regular perturbation [Reference Roberts53, Part III]Footnote
1
.
Being “regular”, the perturbation analysis proceeds to any chosen order N in the ‘small’ derivatives in
$x_0$
. For two examples, the usual leading order homogenizations are the case
$N=2$
, and the so-called second-order homogenizations [Reference Anthoine5, Reference Cornaggia and Guzina18, Reference Hii and El Said31] correspond to the higher-order
$N=4$
. Here we consider arbitrary order N. To ensure the required derivatives and operators exist and have requisite properties, we analyse the systems in Sobolov spaces
$\mathbb {H}_D:=W^{N+1,2}(D)$
for various spatial domains D.
To establish the basis of a slow manifold, consider PDE (2.1), in the Sobolov space with
$\mathbb {H}_{\mathcal D}$
, with 
neglected:
called the cell problem, to be solved with
$\ell _i$
-periodic boundary conditions in
${x_i}$
. The basis applies at each and every
$x_0$
. Hence, at every
$x_0$
, we consider PDE (3.1) on the Sobolov space
$\mathbb {H}_{\mathcal C}$
, where
${\mathcal C}$
is the
$x_1x_2$
-cross-section of
${\mathcal D}$
.
Assumption 3.1. Consider the cell eigen-problem for
${\mathfrak {u}}\in \mathbb {H}_{\mathcal {C}}$
,
with
$\ell _i$
-periodic boundary conditions in
${x_i}$
. We assume that
$\kappa (x_1,x_2) \geq \kappa _{\min }>0$
is regular enough that the set of eigenvectors are complete in
$\mathbb {H}_{\mathcal C}$
.
Equilibria. A family of equilibria of (3.1) is
${\mathfrak {u}}(t,x_1,x_2)=U$
, constant in
$x_1,x_2$
, for every U. This family, in
$\mathbb {H}_{\mathcal C}$
for every
$x_0$
, forms a subspace of equilibria
$\mathbb {E}$
for PDE (3.1).
3.1 Spectrum at each equilibrium
We explore the spectrum about every equilibria in
$\mathbb {E}$
as the spectrum is crucial in identifying the existence and properties of invariant manifolds/subspaces. Since the cell PDE (3.1) is linear, the perturbation problem is identical at every U, namely (3.1) itself.
Example 3.2. Consider the case where the heterogeneity
$\kappa $
is constant. In the vector space
$\mathbb {H}_{\mathcal C}$
, and upon defining wavenumbers
$k_i =2\pi /\ell _i$
, a complete set of linearly independent eigenfunctions of the right-hand side of (3.1) are
$\mathrm{e}^{\mathrm{i}(m k_1 x_1+n k_2 x_2)}$
for every integer
$m,n$
. The spectrum of eigenvalues is then
$\lambda _{mn}=-\kappa (mk_1+nk_2)^2$
. This spectrum has one zero-eigenvalue and the rest are negative:
$\lambda _{mn}\leq -\kappa k_1^2 <0$
(recall
$\ell _1>\ell _2$
so
$k_1<k_2$
).
For general heterogeneity
$\kappa $
satisfying Assumption 3.1, we here establish that the cell eigenvalue problem associated with (3.1) is self-adjoint, and hence has only real eigenvalues. Denote the right-hand side operator as
${{\mathcal L}_0{\kern-1pt}:={\kern-1pt}(\partial _{x_1}{\kern-1pt}+{\kern-1pt}\partial _{x_2})\{\kappa (x_1,x_2) (\partial _{x_1}{\kern-1pt}+{\kern-1pt}\partial _{x_2})\}}$
to be solved with
$\ell _i$
-periodic boundary conditions in
${x_i}$
. Also, we use the inner product
$\langle {\mathfrak {u}},{\mathfrak {v}}\rangle :=\iint _{{\mathcal C}} {\mathfrak {u}}\bar {\mathfrak {v}} \,d{x_1}\,d{x_2}$
over the rectangular cross-section
$\mathcal {C}:=[0,\ell _1]\times [0,\ell _2]$
, where an overbar denotes complex conjugation. Then, for every pair of functions
${\mathfrak {u}},{\mathfrak {v}}\in \mathbb {H}_{\mathcal C}$
, consider
$$ \begin{align*} \langle{\mathcal L}_0{\mathfrak{u}},{\mathfrak{v}}\rangle& =\iint_{{\mathcal C}} \bar{\mathfrak{v}}(\partial_{x_1}+\partial_{x_2})\{\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})\} \,d{x_1}\,d{x_2} \\[5pt]& =\iint_{{\mathcal C}} \bar{\mathfrak{v}}\partial_{x_1}\{\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})\} + \bar{\mathfrak{v}}\partial_{x_2}\{\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})\} \,d{x_1}\,d{x_2} \\[5pt]& =\int_0^{\ell_2} \bigg\{ [\bar{\mathfrak{v}}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})]_0^{\ell_1} -\int_0^{\ell_1} \bar{\mathfrak{v}}_{x_1}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2}) \,d{x_1} \bigg\}\,d{x_2} \\[5pt]&\quad{} + \int_0^{\ell_1} \bigg\{ [{\mathfrak{v}}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})]_0^{\ell_2} -\int_0^{\ell_2} \bar{\mathfrak{v}}_{x_2}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2}) \,d{x_2} \bigg\}\,d{x_1} \\[5pt]& =-\iint_{{\mathcal C}} (\bar{\mathfrak{v}}_{x_1}+\bar{\mathfrak{v}}_{x_2})\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2}) \,d{x_1}\,d{x_2}, \end{align*} $$
as the boundary contributions vanish due to the
$\ell _i$
-periodicity. Reverse the above steps with the roles of
${\mathfrak {u}}$
and
${\mathfrak {v}}$
swapped, to then deduce
$\langle {\mathcal L}_0{\mathfrak {u}},{\mathfrak {v}}\rangle =\langle {\mathfrak {u}},{\mathcal L}_0{\mathfrak {v}}\rangle $
. Thus, the operator
${\mathcal L}_0$
is self-adjoint.
Consequently, every eigenvalue is real and the eigenvectors are pairwise orthogonal (or may be chosen to be so).
The next step is to establish that all the nonzero eigenvalues of
${\mathcal L}_0$
are negative and bounded away from zero. Consider the eigen-problem
${\mathcal L}_0{\mathfrak {v}}=\lambda {\mathfrak {v}}$
for cross-sectional structures
${\mathfrak {v}}\in \mathbb {H}_{\mathcal C}$
. Apply
$\langle \cdot ,{\mathfrak {v}}\rangle $
to the eigen-problem to see that
$$ \begin{align*} \lambda\langle{\mathfrak{v}},{\mathfrak{v}}\rangle &=\iint_{\mathcal C} \bar{\mathfrak{v}}(\partial_{x_1}+\partial_{x_2})\{\kappa ({\mathfrak{v}}_{x_1}+{\mathfrak{v}}_{x_2})\}\,d{x_1}\,d{x_2} \\& =\cdots=-\iint_{\mathcal C} \kappa |{\mathfrak{v}}_{x_1}+{\mathfrak{v}}_{x_2}|^2\,d{x_1}\,d{x_2} \\& \leq-\kappa_{\min}\iint_{\mathcal C} |{\mathfrak{v}}_{x_1}+{\mathfrak{v}}_{x_2}|^2\,d{x_1}\,d{x_2} \,, \end{align*} $$
via integration by parts and the
$\ell _i$
-periodic conditions in
$x_i$
. To deduce an upper bound on the nonzero eigenvalues
$\lambda $
, let us minimize
$\iint _{\mathcal C} |{\mathfrak {v}}_{x_1}+{\mathfrak {v}}_{x_2}|^2\,d{x_1}\,d{x_2}$
in
$\mathbb {H}_{\mathcal C}$
subject to
$\langle {\mathfrak {v}},{\mathfrak {v}}\rangle =\iint _{\mathcal C}|{\mathfrak {v}}|^2\,d{x_1}\,d{x_2}=1$
and that
${\mathfrak {v}}({x_1},x_2)$
has zero mean (being orthogonal to the constant eigenvector corresponding to
$\lambda =0$
). Introduce a Lagrange multiplier
$\mu $
, and use calculus of variations to minimize
$\iint _{{\mathcal C}} |{\mathfrak {v}}_{x_1}+{\mathfrak {v}}_{x_2}|^2 \,d{x_1}\,d{x_2} +\mu ( \iint _{{\mathcal C}}|{\mathfrak {v}}|^2 \,d{x_1}\,d{x_2}-1 )$
. Via integration by parts, it follows that a minimizer must satisfy
$(\partial _{x_1}+\partial _{x_2})^2{\mathfrak {v}}=-\mu {\mathfrak {v}}$
over
${\mathfrak {v}}\in \mathbb {H}_{{\mathcal C}}$
. This is the problem Example 3.2 solves, giving potential minimizers as
$\mathrm{e}^{\mathrm{i}(m k_1 x_1+n k_2 x_2)}$
for multiplier
$\mu =(mk_1+nk_2)^2$
. The zero mean condition requires at least one of
$m,n$
nonzero, and since
$\ell _1>\ell _2$
, so minima are thus linear combinations of
$\mathrm{e}^{\mathrm{i} m k_1 x_1}$
for
$m=\pm 1$
. These give
$\iint _{{\mathcal C}} |{\mathfrak {v}}_{x_1}+{\mathfrak {v}}_{x_2}|^2\,d{x_1}\,d{x_2}/\langle {\mathfrak {v}},{\mathfrak {v}}\rangle =k_1^2=4\pi ^2/\ell _1^2$
. Hence, we obtain the upper bound on the nonzero eigenvalues of
$\lambda \leq -\beta _1$
for
$\beta _1=4\pi ^2\kappa _{\min }/\ell _1^2$
.
In the case of the heterogeneous wave propagation problem,
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
, this subsection’s eigenvalue analysis establishes that for the cross-sectional dynamics, there exists a zero eigenvalue separated from pure imaginary eigenvalues of magnitude at least
$\sqrt {\beta _1}=2\pi \sqrt {\kappa _{\min }}/\ell _1$
.
3.2 Results of slowly varying theory
Let us comment on the preconditions for the next Proposition 3.3 as listed in [Reference Roberts52, Assumption 2]. First, Section 3.1 establishes that the spectrum of
${\mathcal L}_0$
is as required for the Hilbert space
$\mathbb {H}_{\mathcal C}$
. Hence, we decompose the space into two closed
${\mathcal L}_0$
-invariant subspaces, the (slow) centre
$\mathbb {E}_c^0$
and the stable
$\mathbb {E}_s^0$
. Since the operator
${\mathcal L}_0$
is self-adjoint, it generates the required strongly continuous semigroups in the Sobolev space
$\mathbb {H}_{\mathcal C}$
[Reference Carr13, Ch. 6]. Hence, the following result holds for homogenization.
Proposition 3.3 (Homogenization exists and emerges).
By [Reference Roberts52, Proposition 6], for every chosen truncation order N, in a domain
$\mathbb {X}$
of ‘slowly varying solutions’ and after the decay of transients 
for
$\beta '\approx \beta _1$
, the mean field
$U(t,x):=(\ell _1\ell _2)^{-1} \iint _{\mathcal C} {\mathfrak {u}}(t,x,x_1,x_2) \,d{x_1}\,d{x_2}$
satisfies a generalized homogeneous PDE
in terms of some deterministic constants
$K_n$
, to a quantifiable error [Reference Roberts52, (9)–(10) and (23)].
The above mentioned error [Reference Roberts52, (23)] involves too many new parameters to meaningfully detail here—but typically, the most important factor at each and every location
$x=x_0$
is the
$(N+1)$
th spatial derivative of the field
${\mathfrak {u}}$
in a spatial neighbourhood of
$x_0$
. The open subset
$\mathbb {X}$
of the domain
$[0,L]$
is then that part of the domain where the error in (3.2) [Reference Roberts52, (23)], dominantly
$\partial _x^{N+1}{\mathfrak {u}}$
, is small enough for the modelling purposes at hand.
Also, using the specific mean field defined in Proposition 3.3—the usual “cell average”—is optional: one may instead parametrize the macroscale slow manifold in terms of any reasonable alternative amplitude [Reference Roberts53, Section 5.3.3].
A generating function argument [Reference Roberts52, Corollaries 12 and 13] establishes that classic, formal, ‘slowly varying’ analyses of the embedding PDE (2.1) are valid. One just needs to formally treat derivatives in the ‘large’ direction, 
, as ‘small’ in a consistent asymptotic sense [Reference Roberts50, Reference Roberts46].
3.3 Wave propagation in heterogeneous media
In the case of the heterogeneous wave propagation problem,
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
, I conjecture that a backwards approach [Reference Hochs and Roberts32, Reference Roberts54] to the slowly varying theory here would support a cognate version of Proposition 3.3 about the generalized homogeneous wave PDE
$U_{tt}=\sum _{n=0}^N K_n\partial _x^nU$
describing its dispersive macroscale waves.
However, in this case of heterogeneous wave propagation in heterogeneous media, the slow manifold is no longer attractive/emergent (unless there is some significant perturbing dissipative mechanism not included in the ideal wave PDE). Instead, we expect that the slow manifold acts as a ‘centre’ for any fast oscillations arising from the initial conditions. Indeed, conjectured backwards theory [Reference Roberts54, Section 2.5] would more precisely assert that there is a system ‘close’ to the given physical PDE
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
(close to some specified order in some sense) that precisely possesses the constructed homogeneous slow manifold for the macroscale, and that precisely acts as a centre for all nearby dynamics. In this way, the backwards theory would straightforwardly establish accurate modelling over “time scales 
” [Reference Abdulle and Pouchon1, Reference Duerinckx and Gloria24], typically for exponent
$\alpha =N$
, the order of truncation.
Although this article addresses homogenization of linear systems, a comment about homogenizing nonlinear wave systems is appropriate here. As geophysical fluid dynamists have understood for decades, for any given wave system, it is common that a nonlinear slow manifold, such as homogenization seeks, cannot actually exist [Reference Lorenz and Krishnamurthy36]. Nonetheless, backwards theory [Reference Hochs and Roberts32, Reference Roberts54] potentially applies to establish that there exists a wave system ‘close’ to the specified nonlinear wave system, where the ‘close’ systems precisely possess the constructed nonlinear slow manifold and the constructed evolution thereon. Thus, the dynamical systems approach developed here provides a methodology to prevail over subtleties in homogenizing nonlinear wave problems.
4 Example construction of any-order homogenized PDEs
Corollary 13 of [Reference Roberts52] established that a procedure [Reference Roberts46], previously viewed as formal, is indeed a rigorous method to construct the slow manifold modelling PDEs (3.2). Here we use the computer algebra system (cas) version of the procedure [Reference Roberts50], as detailed in generality and examples of the book [Reference Roberts53, Part III].
Appendix A lists and documents the computer algebra code. The code constructs the homogenization for the case of heterogeneity,
$$ \begin{align} \kappa(x_1,x_2):=1/[1+a_1\cos k_1x_1+a_2\cos k_2x_2], \end{align} $$
in terms of microscale wavenumbers
$k_i=2\pi /\ell _i$
for the given microscale periodicities
$\ell _1,\ell _2$
. We write approximations to the slow manifold model of the embedding PDE (2.1) in terms of a ‘mean’ field
$U(t,x)$
that evolves according to an homogenized PDE of the form (3.2).
Quickly verify an approximation. First, Appendix A.1 verifies that an approximate field is
such that the mean field U evolves according to the homogenized PDE
(the coefficient
$1$
of the mean diffusivity is the classic harmonic mean of the specific heterogeneous diffusivity (4.1)). To verify, substitute (4.2) into the governing PDE (2.1) and find the PDE’s residual contains only terms in
$U_{xxx}$
and higher derivatives. We denote this by saying that the residual is zero to an error 
. Since the PDE residual is 
, then theory [Reference Roberts52, Proposition 6] assures us that the derived slow manifold and its evolution (4.2) is correct to error 
.
4.1 Iteration systematically constructs approximations
To construct approximations systematically, we repeatedly compute the residual, which then drives corrections, until the residual is zero to any specified order of error. Theory then assures us that the slow manifold is approximated to the same order of error [Reference Roberts52]. For example, the cas iteratively constructs the improved homogenization (the so-called ‘fourth-order’ case
$N=6$
) that
The code of Appendix A.2 constructs this ‘fourth-order’ homogenization (4.3), and executes in approximately 26 seconds. In the case of the heterogeneous wave propagation problem,
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
, we expect the macroscale homogenized dispersive wave PDE to have
$U_{tt}$
equal to the same right-hand side as the above. Most other rigorous homogenizations invoke a limit ‘
$\epsilon \to 0$
’ for some defined scale separation parameter such as
$\epsilon =\ell _1/L$
. Here such a limit corresponds to the scale separation limit
$k_1, k_2\to \infty $
, in which case, the homogenization (4.3) reduces to the classic, simple, diffusion PDE
$U_t=U_{xx}$
. However, I emphasize that our homogenized PDEs like (4.3) are established for the finite scale separation of real physical materials: finite
$\epsilon ,k_1,k_2$
.
In use, any specific truncation of higher-order homogenization, such as (4.3), may need some asymptotically consistent regularization [Reference Benjamin, Bona and Mahony7, Section 2]. For an example of regularization, let us suppose we decide on the
$N=4$
truncation
${U_t=U_{xx}+\alpha ^2U_{xxxx}}$
with derived parameter 
. The fourth derivative term,
$+\alpha ^2U_{xxxx}$
is undesirably destabilizing. However, to the same order of error, namely 
, the chosen truncation is asymptotically equivalent to the regularized equation
${U_t= (1-\alpha ^2\partial _x^2)^{-1} U_{xx}}$
. This regularization may be alternatively written as
$$ \begin{align*} (1-\alpha^2\partial_x^2)U_t=U_{xx}\,, \quad\text{or}\quad U_t =\frac1{2\alpha} \mathrm{e}^{-|x|/\alpha}\star U_{xx}\,. \end{align*} $$
The spatially nonlocal convolution
$\mathrm{e}^{-|x/\alpha |}\star $
, whether explicit or implicit, stabilizes the regularized chosen truncation while preserving asymptotic consistency [Reference Mercer and Roberts41, cf. Section 2.2]. Such nonlocal PDEs have previously been found desirable in various homogenizations of heterogeneous systems [Reference Cooper17, pp. 2–3].
Optional high-order construction. We may set
$a_2=0$
to get simpler results for a single periodic component in the heterogeneity, namely
${\kappa (x_1)=1/[1+a_1\cos k_1x_1]}$
. In this case, it becomes feasible to construct the homogenization to high-order,
$N=34$
, in spatial derivatives. To construct expressions that are exact in heterogeneity amplitude
$a_1$
, we also need to construct to the same order in
$a_1$
because the evidence is that here, the expansions in heterogeneity amplitude
$a_1$
truncate at just less than the same order as the order of derivatives. For this option (Appendix A.2), the algebra of order
$N=34$
takes roughly one minute of compute time.
In general. The code of Appendix A.2 assumes the heterogeneity
$\kappa \approx 1$
, specifically,
$\kappa =1/[1+\kappa ']$
for some
$\kappa '\propto a$
. Then multiplication by
$\kappa $
is realized as multiplication by
$\sum _{n=0}^{N-1}(-\kappa ')^n$
.
The method [Reference Roberts50] is to start the iteration with the leading ‘trivial’ approximation for the field and its evolution: that is,
$u\approx u^{(1)}: = U$
such that 
. Then the nth iteration is given approximations
$u^{(n)}$
and
$g^{(n)}$
, and starts by evaluating, via the flux
$f^{(n)}: =-\kappa (u^{(n)}_x+u^{(n)}_{x_1}+u^{(n)}_{x_2})$
, and then the residual of the embedding PDE (2.1):
$$ \begin{align*} {\operatorname{Res}}^{(n)}:=u^{(n)}_t+f^{(n)}_x+f^{(n)}_{x_1}+f^{(n)}_{x_2}. \end{align*} $$
First update the evolution via the solvability integral that here is the mean over
$x_1,x_2$
:
$g^{(n+1)}:=g^{(n)}+g'$
, where
$g' :=\overline {{\operatorname {Res}}^{(n)}}$
. Second, update the field by a correction
$u'$
through using two steps to solve
$(\partial _{x_1} +\partial _{x_2}) [ \kappa (u^{\prime }_{x_1} +u^{\prime }_{x_2}) ] ={\operatorname {Res}}^{(n)} +g'$
: step 1 integrates to solve
$v_{x_1}+v_{x_2}={\operatorname {Res}}^{(n)} +g'$
with zero mean to avoid unbounded updates; and step 2 integrates again to solve
$u^{\prime }_{x_1}+u^{\prime }_{x_2}=v/\kappa $
with zero mean. This second “zero mean” ensures preservation of our free choice that the macroscale field parameter U is to be the local mean of the microscale field u. Then update the field approximation with
$u^{(n+1)}: =u^{(n)}+u'$
.
The iterative loop terminates when the PDE residual
${\operatorname {Res}}^{(n)}$
is zero to the specified orders of error. Then the constructed homogenization, such as (4.2) and (4.3), is correct to the same order [Reference Roberts52, Proposition 6].
Analogy with machine learning. The recursive iteration used in the above construction was originally developed nearly 30 years ago [Reference Roberts50]. Let us draw an analogy between this iteration and machine learning algorithms where an ai learns the generic form of the macroscale evolution from many thousands of simulations [Reference Frank, Drikakis and Charissis28, Reference González-Garcia, Rico-Martinez and Kevrekidis29, Reference Linot and Graham35]. In the algorithm here, each recursive iteration is analogous to a layer in a deep neural network, evaluating residuals is analogous to a nonlinear neuronal function, and the linear update corrections are analogous to using weighted linear combination of outputs of one layer as the inputs of the next layer. However here, the recursion is a ‘smart neural network’ in that both the neurons and the ‘linear weights’ are crafted to the problem at hand using the physics encoded in the PDE operator—they are ‘physics-informed’. Further, being analytic, a single analysis encompasses all ‘data points’ in the state space’s domain (here all ‘slowly varying’ functions in
$\mathbb {H}_{{\mathcal D}}$
), not just a finite sample as in machine learning (possibly a biased sample). Consequently, via such algorithms, I contend that mathematicians have for decades been using such recursive iteration for implementing smart analogues of machine learning. Such analytic learning empowers the physical interpretation, validation and verification required by modern science [Reference Brenner and Koumoutsakos10].
4.2 Numerics verify these homogenizations
For a computational verification, let us use the specific instance of the heterogeneity (4.1) with
$a_1=a_2=0.4$
, and wavenumbers
$k_1=21$
and
$k_2=34$
, that is, microscale lengths are
$\ell _1=2\pi /k_1\approx 0.30$
and
$\ell _2=2\pi /k_2\approx 0.18$
(all nondimensional). Consequently, let the macroscale then refer to lengths
${}\gtrsim 2\ell _1\approx 0.6$
(Section 5), that is, wavenumbers
$|k|\lesssim 10$
. We compare the macroscale characteristics of the right-hand side of the homogenized (4.3) with the right-hand side of the original heterogeneous PDE (1.1).
Here the microscale heterogeneity (4.1) (see Figure 2), is near-quasi-periodic with
$k_2/k_1=1.6190\approx (\sqrt (5)+1)/2$
, although this heterogeneity does have a large macroscale period
$2\pi $
(as plotted). Consequently, numerically, we solve heterogeneous PDE (1.1) on the domain of length
$L=2\pi $
with a microgrid of
$714$
points, 
. The “numeric” column of Table 1 gives the resultant eigenvalues of the full heterogeneous PDE (1.1) for macroscale wavenumbers k,
$|k|<10$
. These eigenvalues reflect either the decay rate of the diffusion problem (1.1) or the square of frequencies for the corresponding wave problem (1.3). The homogenization (4.3) of this quasi-periodic heterogeneity may be truncated at any chosen order in
$\partial _x$
, so we compare its predictions for three truncations: the usual leading order homogenization with truncation to errors 
; the ‘second-order’ homogenization with errors 
; and the ‘fourth-order’ homogenization with errors 
. We simply substitute the mode
$U=\mathrm{e}^{\lambda t+\mathrm{i}kx}$
into the truncated right-hand side of (4.3) to obtain rates to compare with that of the full problem. Table 1 shows the two expected results: first, for each wavenumber, the accuracy increases with order of truncation; and second, for each order, as the macroscale wavenumber becomes larger (the macroscale length of the mode decreases), the error increases. These expected trends in the errors are shown clearly in Figure 3 that plots the relative errors on log–log axesFootnote
2
. Table 1 and Figure 3 verify our approach to modelling the homogenized dynamics of quasi-periodic heterogeneous problems.
Figure 2 Example (near) quasi-periodic microscale heterogeneous coefficient (4.1) for the case of
${a_1=a_2=0.4}$
and wavenumbers
$k_1=21$
and
$k_2=34$
.
Figure 3 indicates that the homogenization errors become large for wavenumbers
$|k|\gtrsim 10\approx k_1/2$
. This apparent limit of the resolvable wavenumbers matches the estimate obtained in Section 5, the lower bound (5.1), that the macroscale lengths resolved by the homogenization are longer than roughly twice
$\ell _1$
.
Table 1 Macroscale eigenvalues as a function of macroscale wavenumber k: the “numeric” column lists the computed eigenvalues of (1.1) for spatial discretization 
; the last three columns are the rates predicted for
$U=\mathrm{e}^{\mathrm{i} kx}$
by the homogenization (4.3) truncated to errors 
for orders
$n=2,4,6$
, respectively.
5 Spatial resolution of high-order homogenization
An homogenized PDE, such as (3.2) and (4.3), accurately predicts the evolution of the mean field
$U(t,x)$
provided that, on the scale of the heterogeneity, the spatial gradients are sufficiently small, that is, provided the solutions are sufficiently slowly varying in space. Matous et al. [Reference Matous, Geers, Kouznetsova and Gillman38, Section 3.2] phrased the proviso as “the scale of the microstructural fluctuations, [
$\ell _1$
], must be … much smaller than the macroscopic field fluctuations,
$\ell _M$
” (expressed similarly by Nguyen et al. [Reference Nguyen, Stroeven and Sluys43, Section 2.2.2]). The following question arises: what does “sufficiently small”, “sufficiently slowly varying” or “much smaller” actually mean physically? This section provides an indicative quantitative answer.
Here we show that in a prototypical family of homogenization problems, we only need to require that the
$$ \begin{align} \text{macroscale lengths } \ell_M\gtrsim2\ell_1. \end{align} $$
This lower bound on the macroscale lengths is approximate and indicative because the actual performance of an homogenization depends upon the details of the heterogeneity
$\kappa $
, a chosen truncation N of the homogenized PDE, any chosen specific regularization of the PDE, the chosen desired error of predictions and so on. Recall that this lower bound (5.1) agrees with the growth of errors seen in the example of Figure 3. However here, the numerical coefficient
$2$
in (5.1) arises from a well-established quantitative procedure.
Such a result contrasts with the prevailing expectation provided by other homogenization frameworks: for example, Allaire et al. [Reference Allaire, Lamacz-Keymling and Rauch3, p. 2] comment “no convergence is expected”.
The procedure we invoke is that of high-order asymptotic construction followed by Domb–Sykes plots [Reference Domb and Sykes22, Reference van Dyke67] or Mercer–Roberts plots [Reference Mercer and Roberts40, Appendix] that estimate the distance to the nearest, convergence limiting, singularity in the series approximations. We analyse an homogenized PDE (3.2) and (4.3) by considering its corresponding Fourier transform (
$\partial _x\leftrightarrow \mathrm{i} k$
): 
in terms of the Fourier transform
$\tilde U(t,k)$
of the mean field
$U(t,x)$
. Then the coefficient
${\mathcal K}(k):=\sum _{n=0}^\infty (\mathrm{i}^nK_n)k^n$
is a series in wavenumber k. If we can show convergence for
$|k|<k_*$
for some wavenumber bound
$k_*$
, then we deduce that the homogenized PDE is valid for spatial structures of length scale
$\ell _M>2\pi /k_*$
. We know that the radius of convergence
$k_*$
is the distance to the singularity of
${\mathcal K}(k)$
nearest the origin (
$k=0$
). However, we only know N terms in the series for
${\mathcal K}(k)$
, so Domb–Sykes and Mercer–Roberts plots estimate the distance
$k_*$
by extrapolating the ‘ratio test’ from finite n to infinity, in the cases of a single singularity or a complex conjugate pair of singularities, respectively. For example, this technique has been used to quantitatively predict the spatial limit on PDEs modelling shear dispersion in channels and pipes [Reference Mercer and Roberts40, Reference Mercer and Roberts41]. Figure 4 is an example Mercer–Roberts plot for one homogenization case developed here. The extrapolations to
$1/n=0$
(
$n\to \infty $
) in the figure indicate that the convergence is limited by a complex conjugate pair of singularities at a distance in wavenumber space of
$k_*\approx 1/\sqrt {3.6} \approx 0.52$
.
Figure 4 Macroscale wavenumbers k of convergence are estimated from Mercer–Roberts plots [Reference Mercer and Roberts40, Appendix] for the series
${\mathcal K}(k)$
in the case of heterogeneity (5.2), here with
$a=0.975$
. We use the first 17 terms of the series in
$k^2$
. Extrapolation to
$1/n=0$
estimates the location of a pair of complex conjugate singularities that limit the radius of convergence.
To obtain such plots and estimates, we need to compute the homogenization series (3.2) and (4.3) to high-order. For the case of quasi-periodic heterogeneity, the algebraic complexity explodes combinatorially with increasing order. Consequently, as a representative example, we explore the accessible case of homogenization with a single microscale period, specifically the heterogeneity
$$ \begin{align} \kappa(x):=1/[1+a\cos x]. \end{align} $$
The computer algebra code overviewed in Section 4 constructs the homogenization in this case via the option of setting
$a_2=0$
,
$a_1=a$
and
$k_1=1$
. For algebraic simplicity, in this case, the analysis is nondimensionalized on the microscale length so that nondimensionally, the microscale period
$\ell _1=2\pi $
. Executing the code with this option gives that the homogenized PDE (3.2) begins with
As a first step in investigating such high-order homogenization, observe that the
$a^2$
terms are evidently all of the form
$$ \begin{align*}\tfrac12a^2 \partial_x^4 (-4\partial_x^2)^{n/2-2}\end{align*} $$
(this form holds up to at least order
$n=34$
). Hence, in Fourier space (
$\partial _x\leftrightarrow {\mathrm{i}} k$
), the PDE (5.3) becomes the series forFootnote
3
The divisor
$(1-4k^2)$
in this expression indicates that for small heterogeneity amplitude a, the wavenumbers resolvable by such an homogenization series are limited by the pole singularities at
$k=\pm 1/2$
, that is, a bound on the convergence of the series is
$k_*=1/2$
. Hence, for small a, the homogenized PDE (5.3) potentially resolves all wavelengths longer than
$\ell _*:=2\pi /k_*=4\pi =2\ell _1$
Footnote
4
.
I conjecture that the singularities at wavenumber
$k=\pm 1/2$
are connected to, in wave propagation through heterogeneous media, the well-known phenomenon of spectral gaps, namely frequencies at which no wave can propagate through the underlying medium [Reference Cooper17, p. 2].
Figure 5, for small a, confirms this bound on the wavenumbers and wavelengths. The estimates plotted in Figure 5 are obtained from Section 4 constructing the series (5.3) to errors of order
$35$
in
$\partial _x$
, or equivalently to errors 
. This cas construction provides us with the first 17 nonzero coefficients in the series in
$k^2$
for
${\mathcal K}(k)$
.
-
• For heterogeneity amplitude
$a\lesssim 0.4$
, these coefficients have the same sign. Consequently, Domb–Sykes plots [Reference Domb and Sykes22, Reference van Dyke67], akin to the top panel of Figure 4, predict convergence limiting singularities at real
$k=\pm k_*$
as shown in the left-part of Figure 5. -
• For heterogeneity amplitude
$a\gtrsim 0.4$
, the coefficients have a more complicated sign pattern, indicating the convergence limiting singularity is now a complex conjugate pair in the complex
$k^2$
-plane. Mercer–Roberts plots [Reference Mercer and Roberts40, Appendix], such as Figure 4, estimate the radius and angle of these convergence limiting singularities as shown in the right-part of Figure 5.
Figure 5 In complex wavenumber space, plot estimates of the radius
$k_*$
(circles) and angle
$\theta _*$
(crosses) of the convergence limiting singularity of the homogenization (5.3) as a function of the heterogeneity amplitude a in (5.2). The estimates around
$a\approx 0.4$
are unreliable (see text).
The two distinct behaviours summarized in Figure 5 are likely due to singularities moving as a function of heterogeneity parameter a. We conjecture that for
$a\lesssim 0.3$
, there is a second real singularity for wavenumbers larger than
$k_*$
. As a increases, both these singularities move and appear to collide for
$a\approx 0.35$
. After collision, they separate and move out into the complex k plane (
$a\gtrsim 0.4$
). The Domb–Sykes and Mercer–Roberts estimates are both unreliable near such collision, so the estimates for
$a\approx 0.3$
to
$\approx 0.4$
are unreliable. With this caveat, the evidence of Figure 5 indicates that across all
$0<a<1$
, our homogenization can resolve all wavenumbers
$|k|<1/2$
. That is, high-order homogenization can resolve wavelengths longer than just twice the wavelength of the underlying periodicity.
In the case of the heterogeneous wave propagation problem,
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
, we expect the same restrictions on the spatial resolution because we expect the right-hand side of PDE (4.3) to remain the same for the homogenized waves.
The case of quasi-periodic heterogeneity. In the case of the quasi-periodic heterogeneity (4.1), observe in the homogenization PDE (4.3) that the coefficients in
$a_i$
are divided by wavenumbers
$k_i$
, and with increasing powers of
$k_i$
as the order increases. This suggests that higher-order terms in the series are dominated by the smallest wavenumber
$k_i$
, namely
$k_1$
. That is, the higher-order terms should be dominated by the longest microscale heterogeneity, here that of
$a_1\cos (k_1x)$
. Hence, we expect the limiting wavenumber and wavelengths to be the same as those for the single periodic heterogeneity (5.2). Consequently, expect the bound (5.1) to be valid for this quasi-periodic heterogeneity.
Temporal resolution. The above exploration finds that the spatial resolution of the homogenization appears to vary little with heterogeneity strength a. The temporal resolution is quite a different issue, and it may depend significantly upon a. For example, Figure 4 comes from the case
$a=0.975$
, where the microscale heterogeneity varies by a factor of
$80$
across a microscale period. In similar scenarios, the time scale of decay to the homogenized model, roughly
$1/\beta _1=\ell _1^2/(4\pi ^2\kappa _{\min })$
(end of Section 3.1), may be several orders of magnitude longer than for small heterogeneity. For valid predictions by the homogenization, the macroscale time scale of interest needs to be longer than that of this decay.
Of course, all equilibrium problems satisfy this temporal bound.
6 Three-scale homogenization
There is interest in physical systems possessing a microscale, a mesoscale and a macroscale [Reference Cruz-González, Ramirez-Torres, Rodriguez-Ramos, Penta and Lebon21, Reference Ramirez-Torres, Penta, Rodriguez-Ramos, Merodio, Sabina, Bravo-Castillero, Guinovart-Diaz, Preziosi and Grillo45]. Such three-scale systems are encompassed in the framework of Section 2 when
$\ell _2\ll \ell _1\ll L$
. That is, the cylindrical domain of Figure 1 is shaped like a physical ruler in that the domain is significantly thinner vertically than its width, and its width is in turn significantly thinner than its length. Although Section 5 finds that these length ratios need to be at least two, they are not necessarily much bigger.
In this scenario, the arguments of Section 2 still apply: we solve the phase-shift embedding PDE (2.1) to obtain solutions of the heterogeneous PDEs (1.1) and (2.2). It is the slow manifold support and analysis that differs in detail.
Because this rigorous homogenization could, in principle, be exact at finite scale separation, we find here the expected result that a slow manifold homogenization of a ‘meso-slow’ manifold homogenization of the microscale system is the same as the slow manifold homogenization of the microscale (Sections 3 and 4). That is, in principle, this approach to homogenization is transitive.
The same transitivity of homogenization would also hold in the case of the heterogeneous wave propagation problem,
$u_{tt}=\partial _x\{\kappa (x)u_x\}$
.
6.1 Meso-slow homogenization
To analyse the heterogeneous PDEs (1.1) and (2.2) as a three-scale system, we first establish and create the ‘meso-slow’ manifold homogenization. That is, we use that the domain
${\mathcal D}$
, Figure 1, is thin in
$x_2$
, that is,
$\ell _2\ll \ell _1,L$
. Consequently, we consider the embedding PDE (2.1) by assuming the spatial variations in both
$x_0$
and
$x_1$
are on a relatively large length scale. That is, the derivatives of
${\mathfrak {u}}$
and
$\kappa $
in both
$x_0$
and
$x_1$
are treated as small.
This homogenization is an example of that of a functionally graded material [Reference Anthoine5], because the mesoscale variations in heterogeneity in
$x_1$
are on a length scale
$\ell _1\gg \ell _2$
. Our slow manifold approach systematically incorporates all such variations within the analysis: for example, giving the fourth-order homogenization (6.3)–(6.4). At orders higher than the leading order, we systematically find higher-order derivatives of the larger scale material structure: for example, the fourth-order coefficient in the upcoming PDE (6.4) involves the third derivative
${\mathfrak {K}}"'(x_1)$
of the mesoscale variations. This approach automatically incorporates effects obtained by the “second-order homogenization” of Anthoine [Reference Anthoine5].
Akin to Section 3, we here establish the basis of a slow manifold. Consider PDE (2.1), with 
and 
neglected, and in the Sobolov space
$\mathbb {H}_{[0,\ell _2]}$
with
$\ell _2$
-periodic boundary conditions in
${x_2}$
:
The basis established here applies at each and every
$x_0,x_1$
.
Equilibria. For the cell PDE (6.1), and for every
$x_0,x_1$
, in
$\mathbb {H}_{[0,\ell _2]}$
, clearly
${\mathfrak {u}}(t,x_2)={\mathfrak {U}}$
, constant in
$x_2$
, forms a subspace of equilibria
$\mathbb {E}$
.
Spectrum identifies invariant manifolds. Since the cell PDE (6.1) is linear, the perturbation problem is identical at every equilibria in
$\mathbb {E}$
, namely (6.1) itself. We need to characterize the right-hand side operator in (6.1), 
with
$\ell _2$
-periodic boundary conditions, in
$\mathbb {H}_{[0,\ell _2]}$
.
For example, in the case where the heterogeneity is constant with respect to
$x_2$
, a complete set of linearly independent eigenfunctions of
${\mathfrak {L}}_0$
are
$\mathrm{e}^{\mathrm{i} n k_2 x_2}$
for every integer n. The spectrum of eigenvalues at
$(x_0,x_1)$
is then
$\lambda _{n}(x_1):=-\kappa (x_1)k_2^2n^2$
. This spectrum has one zero-eigenvalue and the rest are negative:
$\lambda _{n}(x_1)\leq -\kappa (x_1) k_2^2 <0$
. Hence, for every
$x_0,x_1$
, the nonzero eigenvalues satisfy
$\lambda _n(x_1)\leq -\kappa _{\min }k_2^2<0$
.
For general heterogeneity, similar deductions to those of Section 3.1 establish that the operator
${\mathfrak {L}}_0$
is self-adjoint. Hence, the eigenvalues are real and eigenfunctions orthogonal. By a similar argument to that of Section 3.1, it follows that the nonzero eigenvalues of the operator are bounded above by
$\lambda \leq -\beta _2$
for
${\beta _2:=\kappa _{\min }k_2^2=4\pi ^2\kappa _{\min }/\ell _2^2}$
.
Slowly varying theory. We proceed analogously to Section 3.2, but now with two space dimensions (
$x_0$
and
$x_1$
) in which solutions are slowly varying. Roberts and Bunder [Reference Roberts and Bunder55] developed rigorous slow manifold theory for spatio-temporal systems that have slow variations in multiple space dimensions. The operator
${\mathfrak {L}}_0$
satisfies the preconditions in [Reference Roberts and Bunder55, Assumption 3] with one zero eigenvalue (
$m=1,\alpha =0$
) and the rest negative. By [Reference Roberts and Bunder55, Proposition 1], and analogous to Proposition 3.3, we are assured that in a regime of slowly varying solutions, the mesoscale field
${\mathfrak {U}}(t,x_0,x_1)$
satisfies a PDE akin to (3.2), but in 2D space, to a quantified error [Reference Roberts and Bunder55, (52)], and upon neglecting exponentially decaying transients (decaying faster than
$\mathrm{e}^{-\beta 't}$
for
$\beta '\approx \beta _2$
).
A generating function argument [Reference Roberts and Bunder55, Sections 3.2–3.3] establishes the validity of formal ‘slowly varying’ analysis of the embedding PDE (2.1).
An example mesoscale homogenization. Appendix B lists and documents cas code to construct the mesoscale homogenization for the example case of heterogeneous diffusivity
$$ \begin{align} \kappa(x_1,x_2):=1/[1/{\mathfrak{K}}(x_1)+a_2\cos k_2x_2] \quad\text{for }k_2:=2\pi/\ell_2\,, \end{align} $$
and some given mesoscale heterogeneity function
${\mathfrak {K}}(x_1)$
. The code derives that in terms of the mesoscale mean
${\mathfrak {U}}(t,x_0,x_1)$
, and in terms of
${\mathfrak {D}}_x:=\partial _{x_0}+\partial _{x_1}$
, the microscale field (cf. (4.2a)) is
Simultaneously, the code derives that the mesoscale mean field evolves according to the quasi-diffusion PDE
As to be expected, the leading-order effective diffusivity on the mesoscale,
${\mathfrak {K}}(x_1)$
, is the harmonic mean over
$x_2$
of the microscale diffusivity (6.2). The theory of this subsection assures us that a truncation of the PDE (6.4) models the dynamics of the embedded PDE (2.1) whenever and wherever the spatial gradients in
$x_0$
and
$x_1$
are ‘small’ enough on the microscale
$\ell _2$
. That is, PDE (6.4) is a mesoscale model.
6.2 Homogenization of the mesoscale model
The mesoscale PDE (6.4) is heterogeneous on the spatial length scale
$\ell _1$
of variations in the coefficient
${\mathfrak {K}}(x_1)$
. For the case
$\ell _1\ll L$
, we here derive the macroscale homogenization of this mesoscale model. We have to choose some truncation of PDE (6.4) to analyse further, so let us choose the leading order truncation
${\mathfrak {U}}_t={\mathfrak {D}}_x\{ {\mathfrak {K}}(x_1) {\mathfrak {D}}_x{\mathfrak {U}} \}$
as the mesoscale model.
To compare with the one-step homogenization, let us consider PDE (6.4) in the case where the mesoscale diffusivity is
$$ \begin{align*} {\mathfrak{K}}(x_1):=1/[ 1+a_1\cos k_1x_1] \quad\text{for }k_1:=2\pi/\ell_1\,. \end{align*} $$
The theory and construction of the homogenization follow the same procedure as detailed before, so let us just summarize here. First, we treat derivatives 
as ‘small’. Then the cross-sectional cell PDE is
to be solved with
$\ell _1$
-periodic boundary conditions in
$x_1$
, and in space
${\mathcal H}_{[0,\ell _1]}$
.
Second, a subspace of equilibria of (6.5) is
${\mathfrak {U}}(t,x_1)=U$
, constant in
$x_1$
. The spectrum of the linear operator 
can be straightforwardly shown to be self-adjoint with one zero eigenvalue, and the rest negative and bounded away from zero,
$\lambda \leq -\beta _1$
, for
$\beta _1:=k_1^2{{\mathfrak {K}}}_{\min }=4\pi ^2{{\mathfrak {K}}}_{\min }/\ell _1^2$
.
Third, slowly varying theory for cylindrical domains, akin to Section 3.2, applies to assure us there exists a slow manifold homogenization of the mesoscale model in terms of a macroscale mean field
$U(t,x_0)$
. Further, the slow manifold emerges exponentially quickly on a time scale of
$1/\beta _1$
, and we may approximate the slow manifold via formal ‘slowly varying’ analysis of the mesoscale PDE (6.4).
Fourth, the construction from the mesoscale PDE
${{\mathfrak {U}}}_t={\mathfrak {D}}_x\{ {\mathfrak {K}}(x_1) {{\mathfrak {D}}}_x{\mathfrak {U}} \}$
to the macroscale homogenization is the specific case of Section 4 with
$a_2=0$
and no variable
$x_2$
. Consequently, corresponding results follow immediately, such as the slow manifold is (4.2a) with
$a_2=0$
, and the macroscale homogenization is
$U_t=U_{xx}$
to errors 
.
Here there is no point proceeding to higher order error in
$\partial _x$
, because we chose the leading order mesoscale PDE at the start of this subsection. To obtain correct higher-order macroscale homogenizations, we would chose a correspondingly higher-order mesoscale homogenization.
7 Conclusion
This article further develops a novel theory and practice for mathematical/asymptotic homogenization to illuminate and resolve many issues in homogenizing heterogeneous PDEs such as (1.1). Section 2 introduces the technique of analysing an ensemble of all phase shifts, and extends it to multi-periodic cases, including the case of quasi-periodic. Future research may be able to extend the approach to interesting classes of random heterogeneity. The key to this approach is to be able to recast the ensemble as a system which is mathematically homogeneous in the macroscale.
Given the mathematically homogeneous macroscale ensemble, Section 3 applies recent developments of dynamical systems theory to newly establish the homogenization as a rigorous slow manifold of the given PDE system. No “
$\epsilon $
” is required so the results hold for finite scale separation. No variational principles are required, so the approach applies to a wide range of physical problems, both dissipative and wave-like. The dynamical systems approach provides a framework that future research and applications may exploit to rigorously derive accurate initial conditions, boundary conditions, forcing projections, uncertainty evaluation and error bounds.
This dynamical systems framework provides a very practical way to construct homogenized models, as in the example of Section 4, and as in Appendix A. It illustrates the ease in using computer algebra to perform and check the homogenization (4.2). The method systematically and straightforwardly extends to higher order to derive so-called “second-order homogenization”: indeed, Appendix A.2 simply proceeds further to explicitly construct the “fourth-order homogenization” (4.3).
One advantage in the computer algebra is that, for some classes of heterogeneities, we can proceed to very high order in the analysis. Section 5 computed to 34th-order and uses the resulting generalized homogenization to quantify how ‘much smaller’ the microscale really has to be compared with the macroscale. Remarkably, the evidence is that, ideally, macroscale lengths can be resolved down to just twice the microscale length! Of course, in practice, there are many confounding aspects, but this factor of two is the mathematically ideal lower bound.
The spatio-temporal resolution of homogenizations could be improved, in future research, by multi-modal/multi-zonal homogenization of the heterogeneous microscale system by adapting approaches developed for shear dispersion [Reference Roberts and Strunin56, Reference Watt and Roberts69]. Such multi-modal/multi-zonal homogenization would be akin to Timoshenko–Ehrenfest beam theory [Reference Timoshenko and Woinowsky-Krieger66, 71].
Lastly, Section 6 explores how the approach developed here includes the case when there are three separated length scales in the physical problem: a microscale, mesoscale and macroscale. The exploration shows how one can analyse such a system in one step from the combined micro and meso scales to the macroscale, or one can analyse the system in two steps from micro to meso, and then from meso to macro. Because sound modelling is transitive, the resulting homogenizations are the same. Further, if one’s ultimate aim is a macroscale spatial discretization to compute predictions, that is, the aim is not really the homogenized PDE, then instead of the two-step process of first homogenization and second spatial discretization, future research could do it all in one step from microscale heterogeneous PDE to macroscale discretization using a process like that introduced in a shear dispersion example [Reference MacKenzie, Roberts, Burrage and Sidje37].
A. Appendix: Computer algebra to construct a homogenized PDE
The following code constructs the homogenization (3.2), discussed in Section 4, of the PDEs (1.1) and (2.1) for the family of problems with diffusivity (4.1) for microscale periodicities
$\ell _1,\ell _2$
. Let us use the computer algebra system Reduce (http://www.reduce-algebra.com) as it is free, flexible and fast [Reference Fateman27]. First improve the printed output:
17 on div; off allfac; on revpri;
18 factor d,df;
Let the arguments of trig function heterogeneity be denoted by
$q_i$
and define wavenumbers
$k_i:=2\pi /\ell _i$
, so the diffusivity
$\kappa (x_1,x_2)$
is the following.
24 depend q1,x1,z; depend q2,x2,z;
25 let { df(q1,x1)⇒k1, df(q2,x2)⇒k2 };
26 kappa:=1/(1+a1*cos(q1)+a2*cos(q2));
Write approximations to the slow manifold model of the embedding PDE (2.1) in terms of a “mean” field
$U(t,x)$
, denoted by uu, that evolves according to 
for whatever dudt happens to be.
36 depend uu,x,t;
37 let df(uu,t)⇒dudt;
This code uses x for
$x_0$
as the two are synonymous in practice.
A.1. Quickly verify a leading approximation
First, we guess and then verify the approximate field (4.2a). This is coded, with ordering parameter d that counts the number of x-derivatives in U, as
50 u:=uu+d*(a1/k1*sin(q1)+a2/k2*sin(q2))*df(uu,x)
51 +d
${\verb|^|}$
2*(a1/k1
$\verb|^|$
2*cos(q1)+a2/k2
${\verb|^|}$
2*cos(q2))*df(uu,x,2);
such that U evolves according to the homogenized PDE 
(the coefficient
$1$
of the mean diffusivity is the classic harmonic mean of the specific heterogeneous diffusivity (4.1)):
58 dudt:=d
${\verb|^|}$
2*df(uu,x,2);
To verify, substitute into the governing PDE (2.1) and find the PDE’s residual is zero to an error 
—counting the derivatives with the order parameter d.
65 let d
${\verb|^|}3\Rightarrow 0$
;
66 flux:=-kappa*(d*df(u,x)+df(u,x1)+df(u,x2))$
67 pde:=df(u,t)+d*df(flux,x)+df(flux,x1)+df(flux,x2);
Since the PDE residual is 
, then the slow manifold (4.2a) is correct to error 
[Reference Roberts51].
The next subsection repeatedly computes the residual, to drive corrections, until the residual is zero to a specified order of error, and hence gives a slow manifold to the same order of error.
A.2. Iteration systematically constructs
Second, we iteratively construct the improved homogenization (4.3). This code executes in approximately four seconds.
86 write "
87 Second, Iteratively Construct
88 -----------------------------";
The terms 
appear, at first, impossible to find with exact algebra, so instead we construct in a power series in the amplitude of the heterogeneity. Let a count the number of amplitude factors arising from
$\kappa $
. Choose the following order of errors in both derivatives and
$a:=\sqrt {a_1^2+a_2^2}$
. It eventuates that for this particular heterogeneity
$\kappa $
, the results to error 
appear exact in heterogeneity amplitude a.
104 ord:=7;
Here assume heterogeneity
$\kappa $
is a reciprocal, but could be more general.
109 rkappa:=sub(a1=a*a1,a2=a*a2,1/kappa);
Optional. We may set
$a_2=0$
to get simpler results for a single periodic component, and change truncation to the same order, now high-order. Order 35 takes roughly one minute of compute time. The high-order coefficients of the homogenized PDE resulting from this option are exact, as far as they go, because the evidence is that the expansions in homogeneity amplitude a extend to at most the same order as the order of derivatives.
121 if 0 then begin a2:=0; a1:=k1:=1; ord:=35; end;
In general. Construct heterogeneous diffusivity
$\kappa $
as the reciprocal of rkappa, presuming
$\kappa \approx 1$
:
128 kappa:=1+for n:=1:ord-1 sum (1-rkappa)
${\verb|^|}$
n$
129 res:=(kappa*rkappa-1 where a
${\verb|^|}$
${\verb|~|} \mathrm{p}\ {=}{>}0$
when p>=ord) ;
130 if res neq 0 then rederr("kappa reciprocal error") ;
Start the iteration from the trivial approximation for the field and its evolution.
137 u:=uu;
138 dudt:=0;
Seek solution to the specified orders of errors.
142 for it:=1:999 do begin write "
143 **** ITERATION ", it;
Progressively truncate the order of the order parameter so that we control the residuals better: the bound in this if-statement is the aimed for ultimate order of error.
150 if it<ord then let {d
${\verb|^|}$
(it+1)=>0, a
${\verb|^|}$
(it+1)=>0};
Compute the PDE residual via the flux, trace printing the length of the residual expression:
155 flux:=-kappa*(d*df(u,x)+df(u,x1)+df(u,x2));
156 pde:=df(u,t)+d*df(flux,x)+df(flux,x1)+df(flux,x2);
157 pde:=trigsimp(pde,combine);
158 write lengthpde:=length(pde);
Update the evolution via the solvability condition from the mean over
$x_1,x_2$
:
163 gd:=-(pde where {sin(
${\verb|~|} $
a)=>0,cos(
${\verb|~|} $
a)=>0});
164 if it<6 then write gd:=gd else write lengthgd:=length(gd);
165 dudt:=dudt+gd;
166 rhs:=pde+gd;
Attempt to solve the update in two steps, each step with two integrals. Use the following operator, given we already made qi depend upon dummy z.
173 if it=1 then begin
174 operator intx; linear intx;
175 let { intx(cos(
${\verb|~|} $
q),z) => sin(q)/(df(q,x1)+df(q,x2))
176 , intx(sin(
${\verb|~|} $
q),z) =>-cos(q)/(df(q,x1)+df(q,x2))
177 };
178 end;
-
• First integrate
$v_{x_1}+v_{x_2}=\texttt {rhs}$
with zero mean. If any intx(1,z) appear, then mistake.184 v:=trigsimp( intx(rhs,z));
185 if not freeof(v,intx(1,z)) then rederr("ABORT");
-
• Second, solve
$\kappa u^{\prime }_z =\kappa (u^{\prime }_{x_1} +u^{\prime }_{x_2}) =v$
, using divide by
$\kappa $
and integrate again. Assume any presence of intx(1,z) is an error that eventually sorts itself out in the iteration, so obliterate here! The update then has mean zero.193 ud:=intx(trigsimp(v*rkappa,combine),z);
194 ud:=(ud where intx(1,z)=>0); % obliterate
Add in the update. Perhaps trigsimp is useful here.
199 u:=u+trigsimp(ud);
Exit iterative for-loop when the PDE residual is zero
204 if pde=0 then write "Success: ", it:= it+1 0000;
205 end;
206 showtime;
207 if pde neq 0 then rederr("Iteration failure");
Write out homogenized PDE, optionally its coefficients, and end the script.
213 dudt:=dudt;
214 if ord>7 then begin cs:=coeff(dudt,d);
215 cs:=for n:=3 step 2 until length(cs)
216 collect part(cs,n)/df(uu,x,n-1);
217 on rounded; off nat; write cs:=cs;
218 on nat; off rounded; end;
219 end;
B. Appendix: Meso-scale homogenization code
Construct the asymptotic expansion of the meso-slow homogenization of three-scale diffusion supported by the theory of Section 6.1.
11 on div; off allfac; on revpri;
12 factor d,b,a2,k2;
The diffusivity
$\kappa (x_1,x_2)$
is the following, where b(n) denotes the nth derivative of
$b(x_1)$
. At the end, we recast in
$\texttt {c}={\mathfrak {K}}(x_1):=1/b(x_1)$
.
18 depend kappa,x1,x2;
19 operator b; depend b,x1;
20 kappa:=1/(b(0)+a2*cos(k2*x2));
21 let { df(b(
${\verb|~|} $
n),x1) => b(n+1) } ;
Write approximations to the slow manifold model of the embedding PDE (2.1) in terms of a “mean” field
$U(t,x)$
that evolves according to 
.
29 depend uu,x,x1,t;
30 let df(uu,t)=>dudt;
This code uses x for
$x_0$
as the two are synonymous in practice. Start from the leading approximation for the field and its evolution.
36 u:=uu$
37 dudt:=0$
Iteratively construct to this order of error in slow derivatives.
43 ordd:=5;
44 for it:=1:999 do begin write "
45 **** Iteration ",it;
Progressively truncate the order of the order parameter so that we better control the residuals.
50 if it<ordd then let d
${\verb|^|}$
(it+1)=>0;
Compute the PDE residual via the flux, trace printing the length of the residual expression:
55 flux:=-kappa*(d*df(u,x)+d*df(u,x1)+df(u,x2));
56 pde:=df(u,t)+d*df(flux,x)+d*df(flux,x1)+df(flux,x2);
57 pde:=trigsimp(pde);
58 write lengthpde:=length(pde);
Solvability updates the evolution:
63 dudt:=dudt+(gd:=-int(pde,x2,0,2*pi/k2)/(2*pi/k2));
Update the field via some integrals:
67 ud:=trigsimp(int(int(pde+gd,x2)/kappa,x2));
68 udx:=sub(x2=2*pi/k2,ud)-sub(x2=0,ud);
69 u:=u+ud-int(ud,x2,0,2*pi/k2)/(2*pi/k2)
70 -udx*(x2-pi/k2+a2/k2/b(0)*sin(k2*x2))/(2*pi/k2);
Exit iterative for-loop when the residual is zero.
75 showtime;
76 if pde=0 then write "Success: ", it:=it+1 0000;
77 end;
78 if pde neq 0 then rederr("Iteration failure");
Write out homogenized PDE.
83 dudt:=dudt;
Recast in terms of
$\texttt {c}={\mathfrak {K}}(x_1):=1/b(x_1)$
.
88 operator c; depend c,x1; factor c;
89 let { df(c(
${\verb|~|} $
n),x1)=>c(n+1), b(
${\verb|~|} $
n)=>df(1/c(0),x1,n) };
90 dudt:=dudt;
Check compact form of the evolution PDE and finish the script.
94 procedure dx(a); d*(df(a,x)+df(a,x1))$
95 err:=-(dudt where d
${\verb|^|}$
5=>0) +dx(c(0)*dx(uu))
96 +a2
${\verb|^|}$
2/k2
${\verb|^|}$
2/2*dx(c(0)
${\verb|^|}$
2*dx(dx(c(0)*dx(uu))));
97 if err neq 0 then rederr("simplification failure");
98 end;
Acknowledgements
I thank Arthur Norman and colleagues who maintain the computer algebra software Reduce. The Australian Research Council Discovery Project grants DP200103097 and DP220103156 helped support this research.
which lead to its computations naturally having errors of approximately 
, is an accurate and stable homogenization for all (non-dimensional) wavelengths longer than 
that vanishes as 
Open access
