Slowly varying, macroscale models emerge from microscale dynamics over multiscale domains
SSlowly varying, macroscale models emerge frommicroscale dynamics over multiscale domains
A. J. Roberts ∗ J. E. Bunder † December 15, 2016
Abstract
Many physical systems are well described on domains which are relativelylarge in some directions but relatively thin in other directions. In this sce-nario we typically expect the system to have emergent structures that varyslowly over the large dimensions. For practical mathematical modelling ofsuch systems we require efficient and accurate methodologies for reducingthe dimension of the original system and extracting the emergent dynamics.Common mathematical approximations for determining the emergent dy-namics often rely on self-consistency arguments or limits as the aspect ratioof the ‘large’ and ‘thin’ dimensions becomes unphysically infinite. Here webuild on a new approach, previously establish for systems which are largein only one dimension, which analyses the dynamics at each cross-section ofthe domain with a rigorous multivariate Taylor series. Then centre mani-fold theory supports the local modelling of the system’s emergent dynamicswith coupling to neighbouring cross-sections treated as a non-autonomousforcing. The union over all cross-sections then provides powerful support forthe existence and emergence of a centre manifold model global in the largefinite domain. Quantitative error estimates are determined from the inter-actions between the cross-section coupling and both fast and slow dynamics.Two examples provide practical details of our methodology. The approachdeveloped here may be used to quantify the accuracy of known approxima-tions, to extend such approximations to mixed order modelling, and to openpreviously intractable modelling issues to new tools and insights.
Contents ∗ School of Mathematical Sciences, University of Adelaide, South Australia. mailto:[email protected] orcid: † School of Mathematical Sciences, University of Adelaide, South Australia. mailto:[email protected] orcid: a r X i v : . [ m a t h . D S ] D ec A PDE models interior plate dynamics 11
B.1 Transform PDEs to ODEs for Taylor coefficients . . . . . . . . . . . 35B.2 Initialise the construction of a transform . . . . . . . . . . . . . . . 36B.3 Iterate to separate slow-fast coordinates . . . . . . . . . . . . . . . 37B.4 Write out the transform . . . . . . . . . . . . . . . . . . . . . . . . 37B.5 Check the generating multinomial form . . . . . . . . . . . . . . . . 37B.6 Write transform in LaTeX . . . . . . . . . . . . . . . . . . . . . . . 38B.7 Transform back to slowly varying . . . . . . . . . . . . . . . . . . . 39
System of large spatial extent in some directions and relatively thin extent in otherdimensions are important in engineering and physics. Examples include thin fluidfilms, flood and tsunami modelling (Noakes, King, and Riley 2006; Bedient andHuber 1988; LeVeque, George, and Berger 2011, e.g.), pattern formation in sys-tems near onset (Newell and Whitehead 1969; Cross and Hohenberg 1993; Westra,Binks, and Water 2003, e.g.), wave interactions (Nayfeh and Hassan 1971; Griffiths,Grimshaw, and Khusnutdinova 2005, e.g.), elastic shells (Naghdi 1972; Mielke 1988;Lall, Krysl, and Marsden 2003, e.g.), and microstructured materials (Romanazzi,Bruna, and Howey 2016, e.g.). There are many formal approaches to mathemati-cally describe, by means of modulation or amplitude equations, the relatively longtime and space evolution of these systems (Dyke 1987, e.g.). This article devel-ops a general approach to illuminate and enhance such practical approximations.Roberts (2015a) originally developed this approach for systems which have onlyone large dimension, and any number of significantly smaller dimensions. Here weconsider the general case where the system has (finite) number of large dimensionsand any number of thin dimensions.The approach is to examine the dynamics in the locale around any cross-section.We find that a truncated Taylor series—a Taylor multinomial—for local spatialstructures is only coupled to neighbouring locales via the highest order resolvedderivative. This coupling as is treated as an ‘uncertain forcing’ of the local dy-namics, and with Assumption 3, we apply non-autonomous centre manifold theory(Potzsche and Rasmussen 2006; Haragus and Iooss 2011, e.g.) to prove the exis-tence and emergence of a slowly varying local model. The theoretical supportprovided by centre manifold theory applies for all cross-sections and so establishesexistence and emergence of a centre manifold model globally over the spatial do-main (Proposition 1) to form an ‘infinite’ dimensional centre manifold (Gallay1993; Aulbach and Wanner 1996; Aulbach and Wanner 2000, e.g.). Section 3 de-velops the methodology for a general linear pde system defined on a general do-2ain consisting of both ‘thin’ and ‘large’ dimensions. In addition to rigorousproofs, Section 3 also establishes a practical construction procedure based upon amultinomial generating function.The new approach derives a novel quantitative estimate of the leading error,equation (52), obtained from the final term of the exact Taylor multinomial (18c).Thus our approach not only provides new theoretical support for established meth-ods such as the method of multiple scales, it extends the methods and providesnew error estimates of the slowly varying model. Interestingly, the theory is stillvalid in boundary layers and shocks, it is just that then the error terms are solarge that the analysis in terms of centre-stable space variations is inappropriate.However, we restrict the current analysis to linear pde systems and will considernonlinear systems in future work.Two examples illustrate the general method and theory of Section 3. First,Section 2 introduces the new approach with the example of a random walker whowalks over a large plane but randomly changes direction among three directions.We construct the emergent mean dynamics of the random walker across the plane,as a Fokker–Planck pde , with a quantifiable error which agrees with the generalform of the error (52). The computer algebra code of Appendix B implementsthe practical construction algorithm for this example and confirms that the mod-elling extends to arbitrary order to rigorously derive a generalised Kramers–Moyalexpansion for the random walker (Pawula 1967, e.g.).Second, Section 4 discusses the more complex example of a two dimensionalheterogeneous diffusion problem—one with a spatially varying diffusivity in a cellu-lar pattern. Via an ensemble of cellular phase-shifted problems, the heterogeneousdiffusion problem is embedded in a family of problems which are homogeneous inthe large dimensions and heterogeneous only in the thin dimensions of the cellularensemble. Then the local dynamics within a cellular cross-section leads to theensemble averaged homogenisation of the original pde . This approach should un-derlie future rigorous modelling of pattern formation problems in multiple spacedimensions.The new methodology developed herein is (cid:15) -free. Although the analysis isbased upon a fixed reference equilibrium, crucially the subspace and centre mani-fold theorems guarantee the existence and emergence of models in a finite domainabout this reference equilibrium. Sometimes such a finite domain of applicabil-ity is large. The only epsilons in this article appear in comparisons with othermethodologies.
This section introduces the novel approach in perhaps the simplest example sys-tem of the effective drift of a vacillating random walker. Section 3 develops theapproach for general linear pde s on general large but thin domains.Consider the example of walker located somewhere in large domain in the xy -plane. At any time the walker steps in one of three directions, as illustratedby different colours in Figure 1. But the walker randomly changes directions,as also illustrated in Figure 1. We model the probability that the walker is atposition ( x , y ) and walking in direction j at time t .Let’s explore the evolution of the probability density function ( pdf ) for therandom walker. Let the pdf be denoted by p j ( x , y , t ) for direction states j =
1, 2, 3 . To avoid some symmetry in the system, the walker cannot move directlybetween the first and third states, but must go indirectly via the second state.3 ( x , y , t ) p ( x , y , t ) p ( x , y , t ) Figure 1: Depending upon state of mind, a walker steps in one of three directions,indicated by different colours and indices j =
1, 2, 3 , and directions indicated bythe intra-plane arrows. However, the walker randomly in time changes state ofmind, indicated by the inter-plane double-headed arrows.Choose constant non-dimensional walking velocities in the three states of v =(
1, 1 ) , v = (−
1, 0 ) and v = ( − ) . The principle of conservation of probabilitygives that the three governing non-dimensional Fokker–Planck pde s are ∂p ∂t = − ∂p ∂x − ∂p ∂y + ( p − p ) , (1a) ∂p ∂t = + ∂p ∂x + ( p − p + p ) , (1b) ∂p ∂t = − ∂p ∂x + ∂p ∂y + ( p − p ) , (1c)for ( x , y ) in some large spatial domain X . For simplicity, in this section we assumethe domain X is convex—a restriction removed in the general analysis of Section 3.An equivalent interpretation of the system of pde s (1) is that for a sort ofheat exchanger when we view p j ( x , y , t ) as the temperature field in three adjacentplates, which is carried by some fluid flow at velocity v j in plate j and diffusesbetween plates j , with only plate j = j =
1, 3 . Ourchallenge then is to find a description of the large time, emergent, heat distribution.The equivalent challenge for the random walker is to describe his/her emergentprobability distribution.We focus upon the emergent solution of the pde s (1) in the interior of thedomain X . This section ultimately finds that the mean pdf , u ( x , y , t ) = ( p + p + p ) , satisfies the anisotropic Fokker–Planck/advection-diffusion pde ∂u ∂t ≈ − ∂u ∂x + ∂ u ∂x + ∂ u ∂y for ( x , y ) ∈ X . (2)Many extant mathematical methods, such as homogenisation and multiple scales(Engquist and Souganidis 2008; Pavliotis and Stuart 2008, e.g.), would derive suchan pde . The main results of this section are to rigorously derive this pde froma ‘local’ analysis, complete with a novel quantitative error formula, and with aninnovative proof that the pde arises as a naturally emergent model from a widevariety of initial and boundary conditions.The pde s (1) would have some boundary conditions specified on the bound-ary ∂ X . A future challenge is to determine the corresponding boundary conditionson ∂ X to be used with the model pde (2) in order to correctly predict the meanfield u (Roberts 1992; Chen, Roberts, and Bunder 2016, e.g.).The analysis is clearer in cross-state spectral modes. Thus transform to fields u ( x , y , t ) := ( p + p + p ) ( the mean field ) ,4 x , y ) = ( X , Y ) p p p Figure 2: schematic diagram of the random walker indicating that we focus onmodelling the pdf dynamics in the locale of a fixed station ( x , y ) = ( X , Y ) . Thevertical cylinder illustrates an example of the locale of the station ( x , y ) = ( X , Y ) . u ( x , y , t ) := ( p − p ) , u ( x , y , t ) := ( p − p + p ) . (3)Equivalently, p = u + u + u , p = u − u and p = u − u + u . Then thegoverning pde s (1) become the separated ‘slow-fast’ system ∂u ∂t = − ∂u ∂x − ∂u ∂y − ∂u ∂x , (4a) ∂u ∂t = − u − ∂u ∂y − ∂u ∂x − ∂u ∂y , (4b) ∂u ∂t = − u − ∂u ∂x − ∂u ∂y + ∂u ∂x . (4c)This form highlights that the difference fields u and u tend to decay exponentiallyquickly, but that interaction between gradients of the mean and difference fieldsgenerates other effects, effects that are crucial in deriving the macroscale Fokker–Planck/advection-diffusion pde (2). Our approach expands the fields in their local spatial structure based around anystation ( x , y ) = ( X , Y ) , and then the results apply to all stations. As commentedearlier, this approach is (cid:15) -free.Fix upon a station in X at ( x , y ) = ( X , Y ) as shown in Figure 2. Consider the pdf fields in the vicinity of ( x , y ) = ( X , Y ) , and denote the collective as vector u ( x , y , t ) := ( u , u , u ) (an unsubscripted u denotes such a vector in R ). For all ( x , y ) ∈ X (convex) invoke a multivariate Taylor’s Remainder Theorem to expressthe fields exactly: u = u ( X , Y , t ) + u ( X , Y , t )( x − X ) + u ( X , Y , t )( y − Y )+ u ( X , Y , x , y , t ) ( x − X ) + u ( X , Y , x , y , t )( x − X )( y − Y )+ u ( X , Y , x , y , t ) ( y − Y )
2! , (5a)where, assuming u is twice differentiable in X , the coefficient ‘derivatives’ are • when m + n < u mn ( X , Y , t ) := ∂ m + n u∂x m ∂y n (cid:12)(cid:12)(cid:12)(cid:12) ( x , y )=( X , Y ) , (5b)5 whereas when m + n = u mn ( X , Y , x , y , t ) := (cid:90) ( − s ) ∂ m + n u∂x m ∂y n (cid:12)(cid:12)(cid:12)(cid:12) ( X , Y )+ s ( x − X , y − Y ) ds . (5c)For definiteness and to avoid a combinatorial explosion of equations, this sec-tion truncates the Taylor expansion of u to the lowest order of interest, namely thequadratic approximation N =
2; Appendix B lists computer algebra code that notonly derives the results summarised here, but also derives corresponding resultsfor arbitrarily specified truncation order N . Local ODEs
Substituting the Taylor expansion (5a) into the governing pde s (4)leads to a set of equations which are exact everywhere. But in some places (namelynear the station ( X , Y ) ) the equations are useful in that remainder terms are neg-ligibly small. We derive a set of linearly independent equations for the coefficientfunctions u mnj in (5) simply by differentiation and evaluation at ( x , y ) = ( X , Y ) (Appendix B): this process is almost the same as equating coefficients of termsin ( x − X ) m ( y − Y ) n , but with care to maintain exactness one finds exact remainderterms. To express the remainder terms, and because of the subsequent evaluationat the station ( x , y ) = ( X , Y ) , the symbols u mnj such that m + n = u mnj ( X , Y , X , Y , t ) ; further, the x and y subscripts in the symbols u mnjx and u mnjy for m + n = u mnjx := ∂u mnj ∂x (cid:12)(cid:12)(cid:12)(cid:12) ( x , y )=( X , Y ) and u mnjy := ∂u mnj ∂y (cid:12)(cid:12)(cid:12)(cid:12) ( x , y )=( X , Y ) . (6)Now, after substituting (5a), the various derivatives of (4a) evaluated at ( x , y ) =( X , Y ) lead to six ode s for the six u coefficients:˙ u = − u − u − u , (7a)˙ u = − u − u − u , (7b)˙ u = − u − u − u , (7c)˙ u = − u x − u x − u x − u y , (7d)˙ u = − u y − u y − u y − u x − u x − u x , (7e)˙ u = − u y − u y − u x − u y − u x . (7f)Similarly, after substituting (5a), the various derivatives of (4b) evaluated at ( x , y ) = ( X , Y ) lead to six ode s for the six u coefficients,˙ u = − u − u − u − u , (7g)˙ u = − u − u − u − u , (7h)˙ u = − u − u − u − u , (7i)˙ u = − u − u y − u y − u x − u x − u x ; (7j)˙ u = − u − u y − u y − u y − u x − u x − u x , (7k)˙ u = − u − u y − u y − u y − u x , (7l)and the various derivatives of (4c) evaluated at ( x , y ) = ( X , Y ) lead to six ode sfor the six u coefficients,˙ u = − u − u − u + u , (7m) The ‘uncertain’ derivatives u mnjx and u mnjy might appear to be simple third-order derivatives,but they are more subtle because of the integral (5c). u = − u − u − u + u , (7n)˙ u = − u − u − u + u , (7o)˙ u = − u − u y − u x − u x + u x . (7p)˙ u = − u − u y − u y + u y − u x − u x + u x , (7q)˙ u = − u − u y − u y + u y − u x + u x , (7r)The functions u mnjx and u mnjy for m + n = u mnjx and u mnjy couplethe dynamics at a station ( X , Y ) with the dynamics at neighbouring stations. Itis by treating terms in u mnjx and u mnjy as ‘uncertain’, time-dependent, inputs intothe local dynamics that we notionally make the vast simplification in apparently reducing the problem from one of an infinite dimensional dynamical system to atractable eighteen dimensional system. For a dynamical system approach to modelling the local dynamics, define the statevector u ( X , Y , t ) = ( u , u , u , u , u , u ) and group the eighteen ode s (7)into the matrix-vector system, of the form d u /dt = L u + r ( t ) , d u dt = L L L L L L L L L L L L (cid:124) (cid:123)(cid:122) (cid:125) L u + r r r (cid:124) (cid:123)(cid:122) (cid:125) r ( t ) (8a)where blanks are zero-blocks, and where L := − − , L := − − − − + , L := − − − − ,(8b)and where r m , n ( t ) in (8a) contain the appropriate terms from (7) involving thedefinite but ‘uncertain’ u mnjx and u mnjy .The hierarchical structure of the matrix in (8) directly reflects the approach:the blocks L m , n directly encode the various derivatives appearing in the originalgoverning pde s (4); whereas the structure of these blocks in L are a consequenceof the multivariate Taylor expansion (5a) and its derivatives. Local slow subspace
The system (8) appears in the form of a ‘forced’ linearsystem, so our first task is to understand the corresponding linear homogeneoussystem obtained by omitting the ‘forcing’ (although here the the ‘forcing’ is theuncertain coupling with neighbouring locales). The corresponding homogeneoussystem is upper triangular, so its eigenvalues are those of L , namely 0, − − O (cid:0) e − t (cid:1) , the system evolves on the 6D slowsubspace of the six eigenvalues 0.Let’s construct this 6D slow subspace. Two eigenvectors corresponding to thezero eigenvalue are found immediately, namely v := (
1, 0, . . . , 0 ) , (9a)7 = ( −
1, 0, 0, 0, 0, 1, 0 . . . , 0 ) . (9b)A further four linearly independent vectors are generalised eigenvectors: v := (
0, 0, − , 1, 0, . . . , 0 ) , (9c) v := (
0, 0, − , 0, 0, − , 0, 0, 0, 1, 0, . . . , 0 ) , (9d) v = ( , 0, 0, −
1, 0, 0, 0, − , 0, 0, 0, 1, 0, . . . , 0 ) , (9e) v := (
0, 0, , 0, 0, 0, 0, −
1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 ) . (9f)Setting the matrix V := (cid:2) v · · · v (cid:3) ∈ R × , the slow subspace is then u = V u where we conveniently choose to use u := ( u , . . . , u ) ∈ R to directlyparametrise the slow subspace because of the form chosen for the eigenvectors v mn .On this slow subspace, from the eigenvectors via u = V u , the difference variables u := (cid:0) − u + u , − u , − u , 0, 0, 0 (cid:1) , u := (cid:0) u − u − u , − u , − u , 0, 0, 0 (cid:1) .Further, on this slow subspace the evolution is determined by the upper triangularmatrix in d u dt = A u = − − −
00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 u (10)The large number of zeros in the ode system (10) reflects the zero eigenvalues, thezero-mean of the y velocities for the three layers, and a useful upper triangularnature. The non-zero elements in the first row give ˙ u = − u + u + u which directly corresponds to the macroscale model pde (2). The novelty of ourapproach is that we now go beyond such a basic model to quantitatively determineremainders—these remainders are quantifiable errors in the model. Near identity coordinate transforms underpin modelling dynamics. In particular,time-dependent coordinate transforms empower understanding of the modelling ofnon-autonomous, and stochastic, dynamical systems (Aulbach and Wanner 1999;Arnold and Imkeller 1998; Roberts 2008, e.g.). This section analogously uses atime-dependent coordinate transformation to separate exactly the slow and fastmodes of the system (8) in the presence of the ‘uncertain forcing’ that couples thedynamics to neighbouring stations. This is the first time the effects of such couplinghave been quantified in multiscale problems with multiple large dimensions.The coordinate transform introduces new dependent variables U ( t ) . In somesense, the new variables u ≈ U so the coordinate transform is ‘near identity’. Let’schoose to parametrise precisely the slow subspace of the system (8) by the vari-ables U : that is, on the subspace where the new stable variables U = U = , weinsist on the exact identity u = U . This choice simplifies subsequent constructionof slowly varying models such as (2).In the coordinate transform, the effects of the uncertain forcings appear asintegrals over their past history. For any µ > e − µt (cid:63) w ( t ) := (cid:90) t e µ ( τ − t ) w ( τ ) dτ . (11)8onsequently, the time derivative d ( e − µt (cid:63) w ) /dt = − µe − µt (cid:63) w + w which is akey property in the upcoming analysis.Well established iteration described elsewhere (Roberts 2008, e.g.) constructsthe coordinate transformation. The details are not significant here, all we needare the results. The computer algebra code of Appendix B, for the case N = U (cid:55)→ u : wherethe ellipses represent many terms we choose not to present so that you can moreeasily appreciate the overall structure, u = U + U − U − U + U + U , (12a) u = U + U + U , (12b) u = U + U + U , (12c) u = U , (12d) u = U , (12e) u = U , (12f) u = U + U − U + e − t (cid:63) u y + · · · − e − t (cid:63) e − t (cid:63) e − t (cid:63) u x , (12g) u = U − U − e − t (cid:63) u y + · · · + e − t (cid:63) e − t (cid:63) u x , (12h) u = U − U − e − t (cid:63) u y + · · · + e − t (cid:63) e − t (cid:63) u x , (12i) u = U − e − t (cid:63) u y − e − t (cid:63) u x − e − t (cid:63) u x − e − t (cid:63) u y − e − t (cid:63) u x , (12j) u = U − e − t (cid:63) u y − · · · − e − t (cid:63) u x , (12k) u = U − e − t (cid:63) u y − e − t (cid:63) u y − e − t (cid:63) u x − e − t (cid:63) u y , (12l) u = U − U − U + U − e − t (cid:63) u x + · · · − e − t (cid:63) e − t (cid:63) e − t (cid:63) u x , (12m) u = U − U − e − t (cid:63) u x · · · − e − t (cid:63) e − t (cid:63) u x , (12n) u = U − U − e − t (cid:63) u y · · · − e − t (cid:63) e − t (cid:63) u x , (12o) u = U + e − t (cid:63) u x − e − t (cid:63) u y − e − t (cid:63) u x − e − t (cid:63) u x , (12p) u = U + e − t (cid:63) u y + · · · − e − t (cid:63) u x , (12q) u = U + e − t (cid:63) u y + · · · − e − t (cid:63) u x . (12r)In these new variables U the original system (8) is exactly the separated system˙ U = U − U + U + e − t (cid:63) u x · · · − e − t (cid:63) e − t (cid:63) u x , (13a)˙ U = − U − e − t (cid:63) u x · · · + e − t (cid:63) u x , (13b)˙ U = − U − e − t (cid:63) u y − · · · + e − t (cid:63) u y , (13c)˙ U = − u x − u y − u x − u x , (13d)˙ U = − u y − u x − u y − u x − u y − u x , (13e)˙ U = − u y − u x − u y − u y − u x , (13f)˙ U = − U − U − U − U − U , (13g)˙ U = − U − U − U , (13h)˙ U = − U − U − U , (13i)˙ U = − U , (13j)˙ U = − U , (13k)9 U = − U , (13l)˙ U = − U + U − U − U − U , (13m)˙ U = U − U − U , (13n)˙ U = U − U − U , (13o)˙ U = − U , (13p)˙ U = − U , (13q)˙ U = − U . (13r)In these new variables, the ode s (13g)–(13r) in this separated system immediatelyshow that all the new stable variables U mn , U mn → t → ∞ . Moreover,they decay exponentially quickly, U mn , U mn = O (cid:0) e − γt (cid:1) for any chosen rate 0 <γ <
1. That is, U = U = is the exact slow subspace for the ‘forced’ system (8).Further, and usefully, by the construction of (12a)–(12f), on this slow subspace u mn = U mn , exactly. Recall the exact Taylor multinomial (5a). Given the exact coordinate trans-form (12), and that U mn , U mn = O (cid:0) e − γt (cid:1) , the Taylor multinomial (5a) establishesthat, based upon the station ( X , Y ) , the mean field u ( x , y , t ) = u ( X , Y , t ) + u ( X , Y , t )( x − X ) + u ( X , Y , t )( y − Y )+ u ( X , Y , x , y , t ) ( x − X ) + u ( X , Y , x , y , t )( x − X )( y − Y )+ u ( X , Y , x , y , t ) ( y − Y ) + O (cid:0) e − γt (cid:1) . (14)Crucially, the left-hand side is independent of the station ( X , Y ) . If the right-hand side was just a local approximation, then the field it generates would dependupon the station ( X , Y ) . But the right-hand side is exact (with its unknown butexponentially quickly decaying transients). This exactness is maintained becausewe keep the remainder terms in the analysis. Consequently, the mean field givenby expression (14) is independent of the station ( X , Y ) despite ( X , Y ) appearing onthe right-hand side.To obtain an exact pde of the slow variations in the mean field u , take the timederivative of (14) and evaluate at ( x , y ) = ( X , Y ) . Remembering the derivative ofthe history convolution, that d/dt ( e − µt (cid:63) w ) = − µe − t (cid:63) w + w , the ode (13a)implies ∂u ∂t = ∂U ∂t + O (cid:0) e − γt (cid:1) = − U + U + U + e − t (cid:63) u x + · · · − e − t (cid:63) e − t (cid:63) u x + O (cid:0) e − γt (cid:1) = − u + u + u + e − t (cid:63) (cid:0) u y + u x − u y − u x + u y − u y − u x (cid:1) + e − t (cid:63) (cid:0) u x + u y + u x + u y + u x + u y + u x + u y + u x (cid:1) + e − t (cid:63) e − t (cid:63) (cid:0) − u y − u x − u y − u x − u y − u x (cid:1) + e − t (cid:63) e − t (cid:63) (cid:0) − u x + u y + u x + u x (cid:1) + O (cid:0) e − γt (cid:1) = − ∂u ∂x + ∂ u ∂x + ∂ u ∂y + r + O (cid:0) e − γt (cid:1) ,10or remainder r := e − t (cid:63) (cid:0) u y + u x − u y − u x + u y − u y − u x (cid:1) + e − t (cid:63) (cid:0) u x + u y + u x + u y + u x + u y + u x + u y + u x (cid:1) + e − t (cid:63) e − t (cid:63) (cid:0) − u y − u x − u y − u x − u y − u x (cid:1) + e − t (cid:63) e − t (cid:63) (cid:0) − u x + u y + u x + u x (cid:1) ,that couples the local dynamics to its neighbouring locales. Consequently, an exactstatement of the mean field u is thus ∂u ∂t = − ∂u ∂x + ∂ u ∂x + ∂ u ∂y + r + O (cid:0) e − γt (cid:1) . (15)In principle, equation (15) is an exact integro-differential equation for the system:the integral part coming from the history convolutions hidden within the coupling r with other locales. The macroscale approximation is to neglect both the transientsand the integral coupling. The rigorous, macroscale, slowly varying, model is thenthe Fokker–Planck/advection-diffusion pde (2) obtained from (15) with O (cid:0) e − γt (cid:1) neglected as a quickly decaying transient, and the uncertain coupling r neglectedas its error.Since the remainder coupling r is a linear combination of the subtle spatialderivatives (6) of second derivatives, here the magnitude of the neglected couplingis that of third order spatial derivatives.The analysis detailed in this section is for the chosen truncation N = N = ∂u ∂t = − ∂u ∂x + ∂ u ∂x + ∂ u ∂y + ∂ u ∂x − ∂ u ∂x∂y + r + O (cid:0) e − γt (cid:1) for some even more voluminous coupling remainder r . Both such pde s are exam-ples of so-called mixed order models, as are the other pde s derived with largertruncation N . The extant mathematical methodologies of homogenisation andmultiple scales promote an aversion to such mixed order models. But our ap-proach rigorously justifies such pde s as mathematical models with the derivedremainder r being the quantifiable error. Inspired by the successful exact modelling of the random walker (heat exchanger)in Section 2, this section establishes analogous exact modelling in the very wideclass (16) of linear systems of pde s in multiple space dimensions. Figure 3 schemat-ically overviews this section.This section develops a rigorous approach to the modelling of fields in multiscaledomains X × Y . We suppose that X is an open domain in R M of large ‘macroscale’extent, called the plate, whereas Y is a relatively small ‘microscale’ cross-sectionaldomain in some Hilbert space. We consider the dynamics of some field u ( x , y , t ) in In some disciplines this pde would be viewed as providing the next term in a Kramers–Moyalexpansion of the mean random walk pdf (Pawula 1967, e.g.). For example, in the case of the modelling of an elastic plate, Y represents the thickness of theplate, and X represents the relatively large width of the plate. Alternatively, if we were analysingthe probability density function of some stochastic system, then space Y could be unboundedso long as there exists a quasi-stationary distribution in Y with operator L having a suitablespectral gap (Roberts 2015b, § pde (16) for u ( x , y , t ) , u ∈ U Lagranges Remainder Theorem (18) m N de s (23) for u ( n ) ( X , y , t ) , u ( n ) ∈ U Generating function G (31)local pde (25) for ˜ u ( ξ , X , y , t ) , ˜ u ∈ U N Linear algebra on U N , §§ m pde s (17) for U ( x , t ) with closure error (52)Figure 3: Flow chart describing the modelling scheme, from the original pde for u ( x , y , t ) to the macroscale model pde for the emergent ‘mean field’ vari-ables U ( x , t ) . Blue rectangles describe different stages of the modelling and reddiamonds describe the method for obtaining one model from another.12 given Hilbert space U (finite or infinite dimensional), where u : X × Y × R → U isa function of M -dimensional position x ∈ X ⊂ R M , cross-sectional position y ∈ Y ,and time t ∈ R (in § y denoted a ‘large’ dimension variable which in thisgeneral theory are all gathered into x , whereas hereafter y denotes position in the‘thin’ dimensions). We suppose the field u ( x , y , t ) satisfies a specified pde in thelinear class ∂u∂t = L u + L ( ) ∂u∂x + · · · + L ( ) ∂u∂x M + · · · = ∞ (cid:88) | k | = L k ∂ kx u , (16)where the M -dimensional (mixed) derivative of order | k | = k + k + · · · + k M is ∂ kx = ∂ | k | ∂x k ∂x k · · · ∂x k M M ,and where the infinite sum is notionally written as being over all possible multi-indices k ∈ N M (as usual, the set of natural numbers N := {
0, 1, 2, . . . } ). However,in application the sum of terms in the pde (16) typically truncate at the first orsecond order derivatives (as in the pde (1)). Nonetheless, our analysis caters forarbitrarily high order pde s. Such sums truncate as we do assume that only a finitenumber of operators L k are non-zero. Thus, such ‘infinite’ sums truncate at somedefinite but unspecified order of | k | .The plate domain X (open) may be finite, or infinite, or multi-periodic. Theplate domain X generally excludes boundary layers and any internal ‘shocks’ or‘cracks’—it need not even be connected. In application, the domain X will bechosen a finite distance from boundary layers or shocks so that boundary layerstructures have decayed enough for the remainder coupling effects (e.g., (15)) tobe below some chosen threshold error on ∂ X . The ‘microscale’ cross-section Y maybe as simple as the index set {
1, 2, 3 } as for the random walker/heat exchanger (2),or the whole of R Y . The operators L k are assumed autonomous and independentof in-plate position x ; they only operate in the ‘microscale’ cross-section Y . Manyproblems which at first appear inhomogeneous in x may be embedded into anensemble of phase-shifted problems such that the resultant ensemble has opera-tors L k independent of position x —see the example of Section 4.For pde s in the general form (16), assume the field u is smooth enough to havecontinuous 2 N derivatives in x , u ∈ C N ( X × Y × R , U ) , for some pre-specifiedTaylor series truncation N . The choice of truncation N is only constrained by thedifferentiability class of field u , and not by any truncation of the infinite k sumin (16).This section, through to the very end of section 3.4, establishes the followingproposition. Proposition 1 (slowly varying pde ) . Let u ( x , y , t ) be governed by a linear pde of the form (16) satisfying Assumption 3. Define the ‘macroscale field’ U ( x , t ) := (cid:104) Z ( y ) , u ( x , y , t ) (cid:105) , U : X × R → R m , for both Z ( y ) and inner product of Defi-nition 4. Then, in the regime of ‘slowly varying solutions’ the macroscale field U satisfies the pde ∂U∂t = N (cid:88) | n | = A n ∂ nx U , x ∈ X , (17) in terms of m × m matrices A n given by (37a) , to a closure error quantified by (52) ,and upon neglecting transients decaying exponentially quickly in time. To cater for the intricacies of this problem we use a lot of symbols. Appendix Asummarises the notation for convenient reference.13 .1 Rewrite the local field
Choose a cross-section at plate station X ∈ X , as shown schematically in Figure 2.Then, generalising (5), invoke the multivariate Lagrange’s Remainder Theorem towrite the field u in terms of an N th degree local Taylor multinomial about thecross-section x = X : u ( x , y , t ) = N − (cid:88) | n | = u ( n ) ( X , y , t ) ( x − X ) n n ! + (cid:88) | n | = N u ( n ) ( X , x , y , t ) ( x − X ) n n ! , (18a)with the multi-index factorial n ! := n ! n ! · · · n M ! , the multi-index power x n := x n x n · · · x n M M , and where • in the first sum, for | n | < N , u ( n ) ( X , y , t ) := ∂ nx u (cid:12)(cid:12) x = X , (18b)and u ( n ) : X × Y × R → U . • whereas for the cases in the second sum of | n | = N , by Lagrange’s RemainderTheorem for multivariate Taylor series (e.g., Taylor 2011, Eq. (1.27)), u ( n ) ( X , x , y , t ) := N (cid:90) ( − s ) N − ∂ nx u (cid:12)(cid:12) X + s ( x − X ) ds , (18c)and u ( n ) : X × X × Y × R → U .As the domain X is open, Lagrange’s Remainder Theorem (18) holds for all x ∈ χ ( X ) ⊆ X where χ ( X ) is any open subset of X such that for all points x ∈ χ ( X ) the convex combination X + s ( x − X ) ∈ χ ( X ) for every 0 (cid:54) s (cid:54) u ( n ) reflects that the functions (18b)–(18c) are fundamentally derivativesof the field u , even though they appear as independent variables in the analysisof this section (and as first posited by Roberts (1988)). Although at the highestorder | n | = N the coefficient functions u ( n ) ( X , x , y , t ) in principle are well defined,in our macroscale modelling the spatial derivatives of these u ( n ) appear as anuncertain part of the modelling closure. The reason for the uncertainty is that,at the highest order | n | = N , the derivatives of these u ( n ) depend upon thefield u ( x , y , t ) between the station X and a point of interest x . Derive exact local DEs
Generalising (7), let’s derive some exact de s (either ode s or pde s depending upon Y ) for the the evolution of the coefficients u ( n ) ( X , y , t ) and u ( n ) ( X , x , y , t ) . The specified pde (16) invokes various spatial derivatives ofthe field u . Since the domain X is open and the multivariate Taylor multino-mial (18a) applies in the open neighbourhood χ ( X ) of each station X , and u ∈ C N ,the Taylor multinomial can be differentiated (cid:96) times, | (cid:96) | (cid:54) N . The multivariateTaylor multinomial (18a) gives, after some rearrangement ∂ (cid:96)x u = N − | (cid:96) | − (cid:88) | n | = u ( n + (cid:96) ) ( x − X ) n n ! + (cid:88) | n | = N n (cid:88) m =( n − (cid:96) ) ⊕ (cid:18) (cid:96)n − m (cid:19) ∂ m + (cid:96) − nx u ( n ) ( x − X ) m m ! , (19)14here in the range of some sums, the notation ( k ) ⊕ means the index vector whose i th component is max ( k i , 0 ) . In deriving expansion (19), the partial derivativeswith respect to x maintain constant X , y and t . Similarly, when taking partialderivatives of any of X , x , y or t , the other variables in this group remain constant.Take the (cid:96) th spatial derivative of the specified pde (16), ∂ (cid:96)x (cid:18) ∂u∂t (cid:19) = ∂ (cid:96)x ∞ (cid:88) | k | = L k ∂ kx u = ⇒ ∂ (cid:0) ∂ (cid:96)x u (cid:1) ∂t = ∞ (cid:88) | k | = L k ∂ (cid:96) + kx u . (20)Then substitute (19) for the spatial derivatives of field u (replacing (cid:96) with (cid:96) + k where appropriate), N − | (cid:96) | − (cid:88) | n | = ∂u ( n + (cid:96) ) ∂t ( x − X ) n n ! + (cid:88) | n | = N n (cid:88) m =( n − (cid:96) ) ⊕ (cid:18) (cid:96)n − m (cid:19) ∂ m + (cid:96) − nx ∂u ( n ) ∂t ( x − X ) m m ! = ∞ (cid:88) | k | = L k N − | (cid:96) + k | − (cid:88) | n | = u ( n + (cid:96) + k ) ( x − X ) n n ! + ∞ (cid:88) | k | = L k (cid:88) | n | = N n (cid:88) m =( n − (cid:96) − k ) ⊕ (cid:18) (cid:96) + kn − m (cid:19) ∂ m + (cid:96) + k − nx u ( n ) ( x − X ) m m ! , (21)This equation is exact for every x ∈ χ ( X ) , for all stations X ∈ X , as the multivariateTaylor multinomial (18a) is exact. However, from the last line of (21), regions ofrapid spatial variation will have large ‘uncertain’ terms involving ∂ kx u ( n ) for | k | > | n | = N .Now set x = X in equation (21) so that all terms containing factors of ( x − X ) vanish. Unless otherwise specified, hereafter u ( n ) denotes u ( n ) ( X , y , t ) when | n | 1) determines theremainder (22b). For all indices | n | = 0, . . . , N , the u ( n ) terms in equation (22a)are evaluated at station X , but the remainder term r n implicitly contains effectsdue to variations along the line joining fixed station X to variable position x viathe integral (18c).Now rewrite the ode s (22a) of the local field derivatives u ( n ) in a form corre-sponding to (8a): ∂u ( n ) ∂t = (cid:88) | k | (cid:54) N , k (cid:62) n L k − n u ( k ) + r n , for | n | (cid:54) N , (23)and with the remainder r n playing the role of time dependent forcing. Equa-tion (23) forms a large system of N := (cid:18) N + MM (cid:19) (24)15 e s in U since there are N possible values for M dimensional n constrained by | n | (cid:54) N ( N = R for the N = ode s).The system of N de s (23) is an exact statement of the dynamics in the lo-cale χ ( X ) of every station X . The system (23) might appear closed, but it iscoupled to the dynamics of neighbouring stations by the derivatives within theremainder (22b): ∂ kx u ( (cid:96) ) for | k | (cid:62) | (cid:96) | = N . Thus system (23) is two faced:when viewed globally as the union over all stations X ∈ X it is a deterministicautonomous system; but when viewed locally at any one station X ∈ X , the inter-station coupling implicit in the remainder r n appears as time dependent ‘forcing’.Our plan is to treat the remainders as ‘uncertainties’ and derive models wherethe effects of the uncertain remainders can be bounded into the precise errorstatement (52) for the models. Roughly, since the remainder is linear in ∂ kx u ( (cid:96) ) ∝ ∂ k + (cid:96)x u when | (cid:96) | = N for slowly varying fields u these high derivatives are smalland so the errors due to the uncertain remainder will be small. If the field u hasany localised internal or boundary layers, then in these locales the errors due tothe uncertain remainder will be correspondingly large and so such locales shouldbe excluded from X . From the specified pde (16) the previous sections dissected out a possibly com-binatorially large system of N de s (23) describing localised dynamics near anystation. This section uses a multinomial generating function to pack all these de sback together into one equation. The theoretical utility is that we then compactlyhandle the many de s all together. The practical utility is that the symbologyconnects with and validates extant, widely used, heuristic methodologies.This section establishes the following proposition. Proposition 2 (linear equivalence) . Let u ( x , y , t ) be governed by the specifiedlinear pde (16) . Then the dynamics at all locales X ∈ X are equivalently governedby the pde ∂ ˜ u∂t = N (cid:88) | k | = L k ∂ kξ ˜ u + ˜ r [ u ] , (25) for the generating function multinomial ˜ u ( X , t ) defined in (26) , and for the ‘un-certain’ coupling term ˜ r [ u ] given by (30) . For every station X ∈ X and time t consider the field u in terms of a local Taylormultinomial (18a) about the cross-section x = X . In terms of the indeterminate ξ ∈ R M , define the generating multinomial˜ u ( X , t ) := N − (cid:88) | n | = ξ n n ! u ( n ) ( X , y , t ) + (cid:88) | n | = N ξ n n ! u ( n ) ( X , X , y , t ) , (26)where this generating multinomial ˜ u , through its range denoted by U N , is implic-itly a function of the indeterminate ξ , and the cross-sectional variable y . Thisgenerating multinomial ˜ u : X × R → U N for the vector space U N := U ⊗ G N where G N := { multinomials in ξ of degree (cid:54) N } , (27) The k th derivative of (18c) with respect to x = X + ξ evaluated at | ξ | = ∂ k + (cid:96)x u ( X , y , t ) = (cid:0) N + | k | N (cid:1) ∂ kx u ( X , X , y , t ) . ⊗ represents the vector space tensor product. In the heat exchangerexample of section 2.2, ( u , u , . . . , u ) ∈ R ⊗ R = R whereas the equivalentgenerating multinomial is ˜ u ∈ U where U = R ⊗ G with ξ = ( ξ , ξ ) . Impor-tantly, the pde (25) for ˜ u is symbolically the same as the specified pde (16) with u ↔ ˜ u , x ↔ ξ , and the addition of a remainder term ˜ r [ u ] . But the derivatives ∂ ξ in the pde (25) are considerably simpler than the derivatives ∂ x in the pde (16)because in pde (25) the derivatives act only on G N , the space of multinomials in ξ of degree at most N . Simpler because, although derivatives are often confound-ingly unbounded in mathematical analysis, here the derivatives ∂ ξ are boundedin G N .Our first task is to show that the generating multinomial (26) satisfies the pde (25) and to derive an expression for the coupling term ˜ r [ u ] . Our second taskis to prove that systematic modelling of the specified pde (16), via pde (25), isequivalent to well-known heuristic procedures expressed in terms of the generatingmultinomial, to some error which is now determined.To find the remainder ˜ r [ u ] , first observe that k derivatives with respect to theindeterminate ξ of the generating multinomial (26) lead to the identity ∂ kξ ˜ u = N − | k | (cid:88) | n | = ξ n n ! u ( n + k ) . (28)Second, take the time derivative of (26) and replace ∂u ( n ) /∂t using (22a): ∂ ˜ u∂t = N (cid:88) | n | = ξ n n ! N − | n | (cid:88) | k | = L k u ( n + k ) + N (cid:88) | n | = ξ n n ! r n = N (cid:88) | k | = L k N − | k | (cid:88) | n | = ξ n n ! u ( n + k ) + N (cid:88) | n | = ξ n n ! r n = N (cid:88) | k | = L k ∂ kξ ˜ u + N (cid:88) | n | = ξ n n ! r n . (29)This is precisely pde (25) of Proposition 2 with coupling remainder˜ r [ u ] = N (cid:88) | n | = ξ n n ! r n = (cid:88) | k | (cid:62) L k N (cid:88) | n | = ξ n n ! (cid:88) | (cid:96) | = N (cid:96) < n + k (cid:18) k + n(cid:96) (cid:19) (cid:2) ∂ k + n − (cid:96)x u ( (cid:96) ) ( X , x , y , t ) (cid:3) x = X (30)upon using expression (22b) for the remainders r n . Expression (30) gives theremainder term appearing in Proposition 2.We now turn to the second key task which is to relate fields in physical spacewith their corresponding field in the generating multinomial space U N . Define theoperator G := N (cid:88) | n | = ξ n n ! ∂ nx x = X , (31)where these subscripted brackets denote evaluation. This operator is denoted by G to signify it determines the generating multinomial corresponding to a given field.17or example, it is straightforward to deduce G [( x − X ) n ] = ξ n when | n | (cid:54) N , and G (cid:2) ( x − X ) n u ( n ) ( X , x , y , t ) (cid:3) = ξ n u ( n ) ( X , X , y , t ) when | n | = N . (32)Thus G operating on the Taylor multinomial expansion (18a) gives G u ( x , y , t ) = ˜ u ( X , t ) ∈ U N = U ⊗ G N . (33)But to use operator G on some general function f ( x ) ∈ C N + , observe that theTaylor expansion of f ( X + ξ ) about x = X gives (cid:2) f ( x ) (cid:3) x = X + ξ = f ( X + ξ ) = N (cid:88) | n | = ξ n n ! ∂ nx f ( X ) + R N ( f ) = G f ( x ) + O (cid:0) | ξ | N + (cid:1) , (34)where R N ( f ) is the appropriate Lagrange remainder term. Thus, G f ( x ) evalu-ates f ( X + ξ ) to a difference O (cid:0) | ξ | N + (cid:1) .Terms of the form “ O (cid:0) | ξ | N + (cid:1) ” are not to be viewed as errors, instead theyrepresent differences. Such differences arise through terms that are of no interestor relevance in this context. The reason is that we are only interested in opera-tions and identities in the multinomial space G N of degree N multinomials. Forexpressions such as f ( X + ξ ) that typically are off G N , it is only the projectiononto G N that is relevant as we only address relations of the components in G N . Aterm “ O (cid:0) | ξ | N + (cid:1) ” reflects such a projection. Establish Proposition 2 The equivalence of indeterminate ξ and space x , towithin the quantified difference, is the key to the equivalence between our rigorousapproach to modelling and the well established heuristics of slow scaling of thespace variables. From its definition (31), the operator G commutes with ∂/∂t andthus, from (33), G ∂u ( x , y , t ) ∂t = ∂ ˜ u ( X , t ) ∂t . (35)However, from (34) we also know that G ∂u ( x , y , t ) /∂t is the Taylor expansionof ∂u ( x , y , t ) /∂t at x = X + ξ with a difference O (cid:0) | ξ | N + (cid:1) . Thus, to a differ-ence O (cid:0) | ξ | N + (cid:1) , the pde (25) of the generating multinomial ˜ u ( X , t ) is equivalentto the pde (16) of original field u ( x , y , t ) when evaluated at x = X + ξ . At x = X (i.e., ξ = ) the pde (25) of ˜ u ( X , t ) and the pde (16) of u ( x , y , t ) are identical.This completes the proof of Proposition 2. Having transformed the physical pde (16) to the equivalent generating functionform (25) we now analyse this later form in order to establish the Slowly Varying pde Proposition 1. The generating function form (25) is symbolically the sameas the physical pde (16) but has two crucial differences: firstly, in G N the deriva-tives ∂ nξ are bounded and finite-D, whereas in X the derivatives ∂ nx are unboundedand ‘infinite’-D; and secondly, it is the presence of the ‘uncertain forcing’ term ˜ r [ u ] that couples the local, approximately finite-D, dynamics to the ‘infinite’-D dynam-ics over the whole domain X .To analyse the ‘uncertainly forced’ system (25) we must first understand theautonomous local system ∂ ˜ u∂t = ˜ L ˜ u with ˜ L := N (cid:88) | k | = L k ∂ kξ ( ˜ L : U N → U N ) . (36)18he invariant subspaces of ˜ L in U N are a key part of our understanding of theautonomous system (36). This subsection establishes Proposition 7 that solutionson the centre subspace of ˜ L satisfy the local, linear, N th order, pde (40).Operator ˜ L is effectively block upper-triangular. The upper triangular nature isdue to the derivatives in the definition (36) of ˜ L . Decompose operator ˜ L = L + ˜ K ,the sum of its ‘diagonal’ part L and its ‘off-diagonal’ part ˜ K := (cid:80) N | k | = L k ∂ kξ .Terms of ˜ K ˜ u in ξ n depend only upon terms of ˜ u in ξ k for | k | > | n | ( | k | (cid:54) N ) , andso ˜ K is zero for all | k | (cid:54) | n | . The diagonal block L of the ‘block upper-triangular’operator ˜ L ensures that the spectrum of operator ˜ L is that of L , but repeated N times (once for each power of ξ in G N ), for the combinatorially large N definedby (24). Assumption 3. We assume the following for the primary case of purely centre-stable dynamics.1. The Hilbert space U is the direct sum of two closed L -invariant subspaces, E c and E s , and the corresponding restrictions of L generate strongly continuoussemigroups (Gallay 1993; Aulbach and Wanner 1996, e.g.).2. The operator L has a discrete spectrum of eigenvalues λ , λ , . . . (repeatedaccording to multiplicity) with corresponding linearly independent (gener-alised) eigenvectors v , v , . . . that are complete ( U = span { v , v , . . . } ).3. The first m eigenvalues λ , . . . , λ m of L all have real part satisfying |(cid:60) λ j | (cid:54) α and hence the m -dimensional centre subspace E c = span { v , . . . , v m } (Chicone2006, Chap. 4, e.g.).4. All other eigenvalues λ m + , λ m + , . . . have real part negative and well sep-arated from the centre eigenvalues, namely (cid:60) λ j (cid:54) − β < − Nα for j = m + m + 2, . . . , and so the stable subspace E s = span { v m + , v m + , . . . } .Although we almost entirely address the case when the Hilbert space U decom-poses into only a centre and a stable subspace, much of the following derivationand discussion applies to other cases that may be of interest in other circum-stances. One may be interested in a centre subspace among both stable andunstable modes, or in a slow subspace corresponding to pure zero eigenvalues (asin the random walker/heat exchanger example of Section 2), or in some other ‘nor-mal mode’ subcentre subspace (Lamarque, Touz´e, and Thomas 2012, e.g.), or inthe centre-unstable subspace , and so on. We primarily focus on the centre subspaceamong otherwise decaying modes as then the centre subspace contains the longterm emergent dynamics from quite general initial conditions (Robinson (1996)called it asymptotically complete). Definition 4. Recall Assumption 3 identifies a subset of m eigenvectors of L which span the centre subspace E c ( ⊂ U ) . • With these eigenvectors define V := (cid:2) v v · · · v m (cid:3) ∈ U × m . • Since the centre subspace is an invariant space of L , define A ∈ R m × m tobe such that L V = V A (often A will be in Jordan form, but it is notnecessarily so). 19 Use (cid:104)· , ·(cid:105) to also denote the inner product on the Hilbert space U , (cid:104)· , ·(cid:105) : U × U → R .Interpret this inner product when acting on two matrices/vectors with el-ements in U as the matrix/vector of the corresponding elementwise innerproducts. For example, for Z , V ∈ U × m , (cid:104) Z , V (cid:105) ∈ R m × m . • Also define Z to have m columns, linearly independent, which both spanthe corresponding centre subspace of the adjoint L † (in the chosen innerproduct), and also are normalised such that (cid:104) Z , V (cid:105) = I m .Developing from this eigen-decomposition of L , we require m N basis vectorsof ˜ L to span the centre subspace of ˜ L , denoted E Nc , in the generating multinomialspace U N . We now establish the existence of suitable basis vectors of ˜ L from thoseof L by considering suitable multinomials in ξ . These basis vectors are typicallygeneralised eigenvectors of ˜ L in the multinomial space U N —they are derived from,but are very different to, the eigenvectors of L Recursively define generalised eigenvectors After solving the basic eigen-problem for A , V and Z , Definition 4, now recursively solve the following se-quence of problems for A n ∈ R m × m and V n ∈ U × m , 0 < | n | (cid:54) N , A n := (cid:88) < | k | , k (cid:54) n (cid:10) Z , L k V n − k (cid:11) , (37a) L V n − V n A = − (cid:88) < | k | , k (cid:54) n L k V n − k + (cid:88) < | k | , k (cid:54) n V n − k A k , (37b) (cid:10) Z , V n (cid:11) = m . (37c)In applications, the m columns of each of these V n contain information about theinteractions between plate-wise gradients of the field u , as felt through the mech-anisms encoded in the L k for | k | > L . Lemma 5. The recursive equations (37) are solvable for A n and V n for all < | n | (cid:54) N .Proof. Consider (cid:104) Z , (37b) (cid:105) . The left-hand side, using the choice (37a) and theorthogonality (37c), becomes (cid:10) Z , L V n (cid:11) − (cid:10) Z , V n A (cid:11) = (cid:10) L † Z , V n (cid:11) − (cid:10) Z , V n (cid:11) A = (cid:10) Z A † , V n (cid:11) − m A = A (cid:10) Z , V n (cid:11) = A m = m .Whereas the right-hand side of (cid:104) Z , (37b) (cid:105) , also using the orthogonality (37c),becomes − (cid:88) < | k | , k (cid:54) n (cid:10) Z , L k V n − k (cid:11) + (cid:88) < | k | < | n | , k (cid:54) n (cid:10) Z , V n − k A k (cid:11) + (cid:10) Z , V A n (cid:11) = − (cid:88) < | k | , k (cid:54) n (cid:10) Z , L k V n − k (cid:11) + (cid:88) < | k | < | n | , k (cid:54) n (cid:10) Z , V n − k (cid:11) A k + (cid:10) Z , V (cid:11) A n = − (cid:88) < | k | , k (cid:54) n (cid:10) Z , L k V n − k (cid:11) + (cid:88) < | k | < | n | , k (cid:54) n m A k + I m A n = − (cid:88) < | k | , k (cid:54) n (cid:10) Z , L k V n − k (cid:11) + A n = m 20y the choice (37a). Consequently the right-hand side of (37b) is in the rangeof the left-hand side. Also, the left-hand side operator has a null space spannedby the m columns of V and so there is enough freedom to impose the normalitycondition (37c) at every step. Lemma 6. For the homogeneous system (36) in the multinomial space U N , a basisfor the centre subspace of ˜ L is the collective columns of ˜ V n := n (cid:88) k = V n − k ξ k k ! , 0 (cid:54) | n | (cid:54) N . (38) Let ˜ V denote the collection of columns of ˜ V n over all n , (cid:54) | n | (cid:54) N , orderedwithin the partial ordering of | n | . Denote the centre subspace of ˜ L , spanned bycolumns of ˜ V , by E Nc .Proof. First prove the space spanned by (38) is invariant: for all n , 0 (cid:54) | n | (cid:54) N ,consider ˜ L ˜ V n = N (cid:88) | (cid:96) | = L (cid:96) ∂ (cid:96)ξ ˜ V n = N (cid:88) | (cid:96) | = L (cid:96) ∂ (cid:96)ξ n (cid:88) k = V n − k ξ k k ! = n (cid:88) k = k (cid:88) (cid:96) = L (cid:96) V n − k ξ ( k − (cid:96) ) ( k − (cid:96) ) ! = n (cid:88) k = k (cid:88) (cid:96) = L (cid:96) V n − k ξ ( k − (cid:96) ) ( k − (cid:96) ) ! = n (cid:88) k = k (cid:88) (cid:96) = L k − (cid:96) V n − k ξ (cid:96) (cid:96) ! = n (cid:88) (cid:96) = ξ (cid:96) (cid:96) ! n (cid:88) k = (cid:96) L k − (cid:96) V n − k = n (cid:88) (cid:96) = ξ (cid:96) (cid:96) ! n − (cid:96) (cid:88) k = L k V n − (cid:96) − k = n (cid:88) (cid:96) = ξ (cid:96) (cid:96) ! n − (cid:96) (cid:88) k = V n − (cid:96) − k A k = n (cid:88) k = (cid:34) n − (cid:96) (cid:88) (cid:96) = ξ (cid:96) (cid:96) ! V n − k − (cid:96) (cid:35) A k = n (cid:88) k = ˜ V n − k A k , (39)which is a linear combination of { ˜ V n } , the columns of ˜ V . Since this identity holdsfor all n , there exists an m N × m N matrix A , composed of N × N blocks of ( m × m ) A n and 0 m , such that ˜ L ˜ V = ˜ V A . (The random walker/heat exchanger matrix A in (10) furnishes an example of such a block upper-triangular matrix A —in a casewith m = L ˜ V n = ˜ V n A + ( lower order ˜ V (cid:96) ) , where ‘lowerorder’ means | (cid:96) | < | n | . Consequently, matrix A is upper-triangular. Further,the N diagonal blocks of A are A . Thus the eigenvalues corresponding to theeigenspace spanned by ˜ V are the centre eigenvalues in A repeated N times. Thesefully account for the centre eigenvalues of ˜ L (counted according to multiplicity).Third, the linear independence of both { ξ n n ! } in G N and the columns of V in U imply, via definition (38), that the columns of ˜ V are linearly independent in U N = U ⊗ G N .Our aim is to show the evolution on the centre subspace has a physicallyappealing, compact, representation directly corresponding to the physical space pde (17). Further, the representation sets up connections to other establishedmethodologies. The next proposition is this desired result.21 roposition 7. Let { U n ∈ R m : n ∈ N M , 0 (cid:54) | n | (cid:54) N } parametrise the centresubspace E Nc via ˜ u = (cid:80) N | n | = ˜ V n U n . For these parameters, define the generatingmultinomial ˜ U ( ξ , t ) := (cid:80) N | n | = ξ n n ! U n ( t ) . Then on the centre subspace E Nc of ˜ L theevolution of the autonomous (36) is governed by the pde ∂ ˜ U∂t = N (cid:88) | n | = A n ∂ nξ ˜ U , (40) for matrices A n constructed by (37a) .Proof. Consider the autonomous (36), ∂ ˜ u/∂t = ˜ L ˜ u : by the parametrisation itsleft-hand side ∂ ˜ u/∂t = (cid:80) N | n | = ˜ V n ∂U n /∂t ; whereas the right-hand side˜ L ˜ u = ˜ L N (cid:88) | n | = ˜ V n U n = N (cid:88) | k | = ˜ L ˜ V k U k = N (cid:88) | k | = k (cid:88) n = ˜ V k − n A n U k ( by (39) )= N (cid:88) | k | = k (cid:88) n = ˜ V n A k − n U k = N (cid:88) | n | = ˜ V n (cid:88) k (cid:62) n , | k | (cid:54) N A k − n U k = N (cid:88) | n | = ˜ V n N − | n | (cid:88) | (cid:96) | = A (cid:96) U (cid:96) + n .By the linear independence of the columns of ˜ V , these two sides are equal iff ∂U n ∂t = N − | n | (cid:88) | (cid:96) | = A (cid:96) U (cid:96) + n for all | n | (cid:54) N . (41)Second, consider the time derivative of the generating multinomial ˜ U : ∂ ˜ U∂t = ∂∂t N (cid:88) | n | = ξ n n ! U n = N (cid:88) | n | = ξ n n ! ∂U n ∂t = N (cid:88) | n | = ξ n n ! N − | n | (cid:88) | (cid:96) | = A (cid:96) U (cid:96) + n = N (cid:88) | n | = N − | n | (cid:88) | (cid:96) | = A (cid:96) ξ n n ! U (cid:96) + n = N (cid:88) | n | = N − | n | (cid:88) | (cid:96) | = A (cid:96) ∂ (cid:96)ξ ξ ( n + (cid:96) ) ( n + (cid:96) ) ! U (cid:96) + n = N (cid:88) | (cid:96) | = A (cid:96) ∂ (cid:96)ξ N − | (cid:96) | (cid:88) | n | = ξ ( n + (cid:96) ) ( n + (cid:96) ) ! U (cid:96) + n = N (cid:88) | (cid:96) | = A (cid:96) ∂ (cid:96)ξ N (cid:88) | k | = ξ k k ! U k = N (cid:88) | (cid:96) | = A (cid:96) ∂ (cid:96)ξ ˜ U .Thus the pde (40) holds on the autonomous centre subspace of ˜ L in U N .This result establishes that locally, and in the absence of coupling with nearbylocales, a macroscale field together with its slowly-varying derivatives exist thatevolve consistently with the expected macroscale pde .As an example we return to the random walker of Section 2 and show how toconstruct all ˜ V n for | n | (cid:54) N = 6. We use initial basis vectors V = v = ( 1, 0, 0 ) (the eigenvector of L with zero eigenvalue) and Z = ( 1, 0, 0 ) (the eigenvectorof adjoint L † with zero eigenvalue). Either, we first recursively calculate all A n V n from (37) and use definition (38) to construct all ˜ V n ; or, wecalculate A n and ˜ V n directly with A n := (cid:88) < | k | , k (cid:54) n (cid:10) Z , L k ˜ V n − k (cid:11) ξ = , (42a) L ˜ V n − ˜ V n A = − (cid:88) < | k | , k (cid:54) n L k ˜ V n − k + (cid:88) < | k | , k (cid:54) n ˜ V n − k A k , (42b) (cid:10) Z , ˜ V n (cid:11) = ξ n n ! , (42c)obtained from combining (37) and (38). In either case, A , A , A = A =− , A = and A = (in agreement with matrix A in (10)) and˜ V = V = ( 1, 0, 0 ) ,˜ V = V + V ξ = ( ξ , 0, − ) ,˜ V = V + V ξ = ( ξ , − 1, 0 ) ,˜ V = V + V ξ + V 00 12 ξ = ( ξ , 0, − − ξ ) ,˜ V = V + V ξ + V ξ + V ξ ξ = ( ξ ξ , − ξ , − ξ ) ,˜ V = V + V ξ + V 00 12 ξ = ( ξ , − ξ , ) .These three dimensional generating basis vectors of order N = N = 18 dimensional eigenvectors v mn in equation (9).Straightforward substitution confirms that A n and basis vectors ˜ V n satisfy iden-tity (39). Our aim is not to model the autonomous (36), but the exact system (25) withits uncertain coupling ˜ r [ u ] to neighbouring locales. Let’s proceed to project theuncertain coupling by treating it as an arbitrary , time dependent, forcing. Theresult of this subsection completes the proof of Proposition 1. Change basis to centre and stable variables Write ˜ u = ˜ VU + ˜ WS wherethe centre variables U = (cid:2) U n (cid:3) parametrise the centre subspace, and the sta-ble variables S parametrise the stable subspace. Here, ˜ W is a differential ma-trix operator containing eigenvectors that form a basis for the multinomial stablesubspace E Ns and is analogous to ˜ V , the differential operator containing eigen-vectors forming a basis fo the multinomial centre subspace E Nc (Lemma 6). Asdetailed in Assumption 3, the eigenvectors of the stable subspace have eigenval-ues λ j where (cid:60) λ j (cid:54) − β < − Nα , indicating rapid exponential decay of thesemodes and the emergence of the centre subspace over long time scales (with itseigenvalues |(cid:60) λ j | (cid:54) α ). Despite the rapid decay of the stable modes, when forcedby coupling with neighbouring locales, their influence on the dynamics of the sys-tem is not negligible in general. Here we derive a slow macroscale model whichquantifies the effects of the coupling.Analogous to ˜ V , ˜ W is associated with the following properties: • ˜ W forms a basis for the multinomial stable subspace E Ns of ˜ L ; and • there exists an operator B : E Ns → E Ns such that ˜ L ˜ W = ˜ WB and all eigen-values of B have real part (cid:54) − β (the eigenvalues of B are λ m + , λ m + , . . .).23e separate the ‘forcing’ in system (25) into components in E Nc and E Ns : ˜ r = ˜ V ˜ r c + ˜ W ˜ r s . Then by the linear independence of the complete basis { ˜ V , ˜ W } , the‘forced’ system (25) separates to ∂ U ∂t = A U + ˜ r c , (43a) ∂ S ∂t = BS + ˜ r s . (43b)Now consider the stable variables S . Since ˜ L generates a continuous semigroup(Assumption 3), so does its restriction B , and so we rewrite (43b) in the integralequation form S ( t ) = e B t S ( ) + (cid:90) t e B ( t − τ ) ˜ r s ( τ ) dτ = e B t S ( ) + e B t (cid:63) ˜ r s , (44)as convolutions f ( t ) (cid:63) g ( t ) = (cid:82) t f ( t − τ ) g ( τ ) dτ . Since all eigenvalues of B havereal part (cid:54) − β , then for some decay rate γ ∈ ( α , β ) S ( t ) = e B t (cid:63) ˜ r s + O (cid:0) e − γt (cid:1) , written S ( t ) (cid:39) e B t (cid:63) ˜ r s (45)upon invoking the following definition. Definition 8. Define f ( t ) (cid:39) g ( t ) to mean f − g = O (cid:0) e − γt (cid:1) as t → ∞ for someexponential decay rate γ such that α < γ < β .Consequently, equation (45) quantifies how the local stable variables S areforced by coupling with neighbouring locales via the remainder effects in ˜ r s . The centre subspace dynamics with remainder As invoked in Proposi-tion 1, define the macroscale amplitude field of slowly varying solutions by theprojection U ( x , t ) := (cid:104) Z , u ( x , y , t ) (cid:105) , (46)which as yet is distinct from the multinomial local centre variables U . In orderto discover how the amplitude field U ( x , t ) evolves, our task is to now relate thefield U ( x , t ) to the local centre subspace variables U .As a preliminary step, and since Z is independent of ξ , for any index (cid:96) simplify (cid:10) Z , ∂ (cid:96)ξ ˜ V (cid:11) ξ = = (cid:68) Z , (cid:2) ∂ (cid:96)ξ ˜ V n (cid:3) ξ = (cid:69) ( by Lemma 6, with matrix index n )= (cid:10) Z , (cid:2) V n − (cid:96) for n (cid:62) (cid:96) , else 0 (cid:3)(cid:11) ( by (38) )= (cid:2) (cid:104) Z , V n − (cid:96) (cid:105) for n (cid:62) (cid:96) , else 0 m (cid:3) ( by Defn 4 )= (cid:2) I m for n = (cid:96) , else 0 m (cid:3) ( by Defn 4 and (37c) ) . (47)Consequently, with (cid:96) = , (cid:104) Z , ˜ V (cid:105) ξ = A is the first row of blocks of A , namely the m × m blocks A n in appropriate order. Recalling that U = (cid:2) U n (cid:3) , an identity tobe used shortly is then that (cid:104) Z , ˜ V (cid:105) ξ = A U = N (cid:88) | n | = A n U n . (48)Recall that the multinomial ˜ u ( X , t ) is the Taylor expansion in ξ of u ( X + ξ , y , t ) to differences O (cid:0) | ξ | N + (cid:1) . Hence the field u ( X , y , t ) = ˜ u ( X , t ) (cid:12)(cid:12) ξ = and so24acroscale amplitude U ( X , t ) = (cid:104) Z , ˜ u ( X , t ) (cid:12)(cid:12) ξ = (cid:105) . Take the temporal derivativeof U ( x , t ) at cross-section x = X , ∂U ( X , t ) ∂t = (cid:28) Z , ∂ ˜ u ( X , t ) ∂t (cid:29) ξ = = (cid:28) Z , ˜ V ∂ U ∂t (cid:29) ξ = + (cid:28) Z , ˜ W ∂ S ∂t (cid:29) ξ = ( using separated system (43) )= (cid:104) Z , ˜ V A U (cid:105) ξ = + (cid:104) Z , ˜ V ˜ r c (cid:105) ξ = + (cid:104) Z , ˜ WBS (cid:105) ξ = + (cid:104) Z , ˜ W ˜ r s (cid:105) ξ = ( using convolution solution (45) )= (cid:104) Z , ˜ V (cid:105) ξ = A U + (cid:104) Z , ˜ V ˜ r c (cid:105) ξ = + (cid:104) Z , ˜ WB e B t (cid:63) ˜ r s (cid:105) ξ = + (cid:104) Z , ˜ W ˜ r s (cid:105) ξ = + O (cid:0) e − γt (cid:1) ( using identity (48) and Defn 8 ) (cid:39) N (cid:88) | n | = A n U n + (cid:104) Z , ˜ V ˜ r c (cid:105) ξ = + (cid:104) Z , ˜ WB e B t (cid:63) ˜ r s (cid:105) ξ = + (cid:104) Z , ˜ W ˜ r s (cid:105) ξ = . (49)For this result to form a pde for the macroscale field U we need to write thecentre subspace parameters U n in terms of spatial derivatives of U . Identity (28)ensures ∂ kξ ˜ u (cid:12)(cid:12) ξ = = u ( k ) = ∂ kx u (cid:12)(cid:12) x = X for | k | (cid:54) N . Then the (cid:96) th spatial derivativeof U ( x , t ) at x = X is ∂ (cid:96)x U (cid:12)(cid:12) x = X = ∂ (cid:96)x (cid:104) Z , u (cid:105) (cid:12)(cid:12) x = X = (cid:104) Z , ∂ (cid:96)x u | x = X (cid:105) = (cid:104) Z , ∂ (cid:96)ξ ˜ u (cid:105) ξ = ( using ˜ u = ˜ VU + ˜ WS )= (cid:104) Z , ∂ (cid:96)ξ ˜ V (cid:105) ξ = U + (cid:104) Z , ∂ (cid:96)ξ ˜ WS (cid:105) ξ = ( using (47) and the solution (45) )= (cid:2) I m for n = (cid:96) , else 0 m (cid:3) U + (cid:104) Z , ∂ (cid:96)ξ ˜ W e B t (cid:63) ˜ r s (cid:105) ξ = + O (cid:0) e − γt (cid:1) (cid:39) U (cid:96) + (cid:104) Z , ∂ (cid:96)ξ ˜ W e B t (cid:63) ˜ r s (cid:105) ξ = . (50)The above shows that, discounting exponential transients, U (cid:96) is the (cid:96) th spatialderivative of the amplitude field U ( x , t ) evaluated at x = X , with a remainderterm determined from the forcing coupling.Substituting equation (50) into (49) gives ∂U∂t (cid:12)(cid:12)(cid:12)(cid:12) x = X (cid:39) N (cid:88) | n | = A n ∂ nx U | x = X + ρ , (51)where the remainder ρ = (cid:104) Z , ˜ V ˜ r c (cid:105) ξ = + (cid:104) Z , ˜ WB e B t (cid:63) ˜ r s (cid:105) ξ = + (cid:104) Z , ˜ W ˜ r s (cid:105) ξ = − N (cid:88) | n | = A n (cid:104) Z , ∂ nξ ˜ W e B t (cid:63) ˜ r s (cid:105) ξ = . (52)The pde (51) applies for every station X in the domain X . Strictly, the‘ pde ’ (51) is a coupled differential-integral equation: the dynamics at each sta-tion X being coupled by the gradients and their history convolution integrals oc-curring within the remainder (52). But when the uncertain remainder term is neg-ligible, as in slowly varying regimes where the remainder ρ is O (cid:0) (cid:80) | n | = N + | ∂ nx u | (cid:1) ,25 x (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109) (cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:109)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102)(cid:102) Figure 4: microscale periodic cells of size h × h in 2D are represented by a spatiallyvarying, doubly h -periodic, diffusion coefficient K ( x ) in the pde (53) over somelarge domain.then equation (51) reduces to the plate pde closure (17). This completes theargument that establishes Proposition 1. Many engineering structures have microscale structure, such as the windings inelectrical machinery (Romanazzi, Bruna, and Howey 2016, e.g.), electromagnetismin micro-structured material (Craster 2015; Niyonzima 2014, e.g.), and slow Stokesflow through porous media (Brown, Popov, and Efendiev 2011, e.g.). The engi-neering challenge is to understand the dynamics on a scale significantly largerthan the micro-structure. Homogenization, via multiple length scales, is the com-mon approach (Gustafsson and Mossino 2003; Engquist and Souganidis 2008, e.g.).Building on the 1D case (Roberts 2015a, § x = ( x , x ) be the spatial coordinates, then we seek to model on the macroscale—that is, across many cells—the diffusion in time t of the temperature field u ( x , t ) satisfying the diffusion pde ∂ u ∂t = ∂∂ x (cid:20) K ( x ) ∂ u ∂ x (cid:21) + ∂∂ x (cid:20) K ( x ) ∂ u ∂ x (cid:21) . (53)The spatially varying diffusion coefficient K ( x ) is here assumed to be doubly pe-riodic over a length h as illustrated schematically by Figure 4; that is, K ( x + ph , x + qh ) = K ( x , x ) for integer p , q . Here we use upright roman letters forfield u and space x to distinguish these direct physical quantities from those of themathematical analysis which use the maths font u and x for closely related butdifferent quantities.Ensemble averages provides our route to rigorous modelling. Let’s embed thespecific problem in the family of problems of all phase shifted versions of the26 − 20 20 h x x ~ y h h ~ x y y Figure 5: schematic diagrams both illustrating that at every macroscale posi-tion x ∈ X , the pde (55) has a small square cross-section Y = [ h ] of onecell. Solutions of pde (54) are the values of field u on any 2D slice through thesepictures of y = ( x + φ ) mod h .material. That is, for all microscale phase shifts φ , 0 (cid:54) φ , φ < h , seek the fieldu ( x , φ , t ) that satisfies the pde ∂ u ∂t = ∂∂ x (cid:20) K ( x + φ ) ∂ u ∂ x (cid:21) + ∂∂ x (cid:20) K ( x + φ ) ∂ u ∂ x (cid:21) . (54)The original pde (53) and its solution is included in this family as the phase shift φ = version of pde (54) and its solution. The second step in the embeddingis to write the solution field u in terms of a new field u ( x , y , t ) that is a functionof macroscale coordinates x ∈ X , microscale ‘cell’ coordinates y ∈ Y = [ h ] ,and time t . Hereafter let’s consider the diffusivity K to be only a function of themicroscale cell coordinate y , that is, K ( y ) . Then consider solutions u ( x , y , t ) tothe pde ∂ t u = ∂ y (cid:2) K ( y ) ∂ y (cid:3) u + ∂ y (cid:2) K ( y ) ∂ y (cid:3) u + (cid:2) K y + K∂ y (cid:3) u x + (cid:2) K y + K∂ y (cid:3) u x + K ( y ) u x x + K ( y ) u x x . (55)Elementary algebra shows that solutions of this pde (55), via u ( x, φ , t ) = u ( x, x + φ , t ) , give solutions to the pde (54)—and hence (53)—and vice-versa.We need to consider boundary conditions for both domains X and Y . First, theslowly varying theory of Section 3 applies usefully in a macroscale domain X fromwhich boundary layers and shocks have been excised, and addresses the generalsolutions in such a domain X . Thus macroscale boundary conditions on ∂ X are notneeded to apply the theory. Second, on the microscale Y , the diffusion coefficient is h -periodic in y , so we correspondingly require the field u ( x , y , t ) to be h -periodicin y . Then identities such as y = x + φ are implicitly to be interpreted modulo h in both components. This periodicity in the cell structure is sufficient to form awell-posed problem in the cross-sectional, cell, microstructure.The embedding pde (55) is in the linear class (16) of Section 3 with M = L = ∂ y (cid:2) K ( y ) ∂ y (cid:3) + ∂ y (cid:2) K ( y ) ∂ y (cid:3) , L ( ) = (cid:2) K y ( y ) + K ( y ) ∂ y (cid:3) L ( ) = (cid:2) K y ( y ) + K ( y ) ∂ y (cid:3) L ( ) = L ( ) = K ( y ) ,and all other L k zero. Other crucial requisites are those of Assumption 3 upon theproperties of L with boundary conditions of double h -periodicity in y ∈ Y ; that27s, upon the basic properties of the fundamental cell problem. The self-adjoint celleigen-problem L v = ∂ y (cid:2) K ( y ) v y (cid:3) + ∂ y (cid:2) K ( y ) v y (cid:3) = λv is well known and for diffusion K ( y ) (cid:62) K min > m = v ( y ) = constant over the cell Y , and the othereigenvalues are real and negative, λ j (cid:54) − β for β = π K min /h . The correspondingset of eigenfunctions are orthogonal and complete in the Hilbert space of smooth L functions on Y . Consequently the requisite semigroups exist (Carr 1981, Ch. 6,e.g.). Hence Proposition 1 applies. Set Z ( y ) = /h so that the macroscalefield U ( x , t ) := h (cid:82)(cid:82) Y u ( x , y , t ) d y is a cell-mean. Then, for any chosen order oftruncation N , the field U satisfies the pde ∂U∂t = N (cid:88) | n | = A n ∂ n U∂ x n , x ∈ X , (56)for some 2 × A n depending upon K ( y ) , to a closure remainder errorquantified by (52), and upon neglecting transients decaying roughly like O (cid:0) e − βt (cid:1) as t → ∞ .The pde (56) is the homogenisation of the original pde (53), generalised to anyorder of truncation. Further, the novel expression (52) would be used to quantifythe remainder error in any large scale modelling. Vibrations of an inhomogeneous plate If, instead of the inhomogeneousdiffusion (53), suppose we wanted to model the vibrations of an inhomogeneousplate satisfying the corresponding, second order in time, pde ∂ u ∂t = ∂∂ x (cid:20) K ( x ) ∂ u ∂ x (cid:21) + ∂∂ x (cid:20) K ( x ) ∂ u ∂ x (cid:21) . (57)Then all of the above algebra would be effectively the same except for two as-pects. First, the exponential emergence of the model from all initial conditions,expressed in terms “ O (cid:0) e − γt (cid:1) ”, would be replaced by long-lived oscillations e iω j t of relatively high frequency. Second, also the exponentially decaying convolutionsin the remainder term r would turn into tricky convolutions with oscillating fac-tors e iω j t , potentially causing an algebraic growth in the error of the model. Bothof these mechanisms would be ameliorated by any small viscous damping or smallradiative damping that is physically present but omitted from the mathematical pde (57). Consequently, in practice and with care the approach here should helpmodel the vibration of inhomogeneous plates. This article develops a general theoretical approach to supporting the much in-voked practical approximation of slow variations in space. The key idea, suggestedby Roberts (2015a), is to examine the dynamics in the locale around any cross-section. We find that a Taylor series approximation to the dynamics is only cou-pled to neighbouring locales via the highest order resolved derivative. Treating One subtlety in the result is that in any given physical realisation, the ensemble of initialconditions has to be chosen so that the remainder error becomes small. If an ensemble of initialconditions were chosen poorly, then the remainder error would stay large over space-time. X ⊂ R .This approach opens much for future research. It may be able to illuminatethe thorny issue of providing boundary conditions to slowly varying models ofproblems such as shells, plates and Turing patterns (Segel 1969; Roberts 1992;Mielke 1992, e.g). Acknowledgement The Australian Research Council Discovery Project grantDP120104260 and DP150102385 helped support this research. We thank ArthurNorman and colleagues who maintain the Reduce software. References Arnold, Ludwig and Peter Imkeller (1998). “Normal forms for stochasticdifferential equations”. In: Probab. Theory Relat. Fields doi : (cit. on pp. 8, 29).Aulbach, Bernd and Thomas Wanner (1996). “Integral manifolds forCaratheodory type differential equations in Banach spaces”. In: Six Lectureson Dynamical Systems . Ed. by B. Aulbach and F. Colonius. World Scientific,Singapore, pp. 45–119 (cit. on pp. 2, 19).— (1999). “Invariant foliations for Caratheodory type differential equations inBanach spaces”. In: Advances of Stability Theory at the End of XX Century .Ed. by V. Lakshmikantham and A. A. Martynyuk. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.45.5229&rep=rep1&type=pdf .Gordon & Breach Publishers (cit. on p. 8).— (2000). “The Hartman–Grobman theorem for Caratheodory-type differentialequations in Banach spaces”. In: Nonlinear Analysis doi : (cit. on p. 2).Bedient, P. B. and W. C. Huber (1988). Hydrology and floodplain analysis .Addison–Wesley (cit. on p. 2). 29ridges, Thomas J. and Daniel J. Ratliff (June 2016). “Double criticality and thetwo-way Boussinesq equation in stratified shallow water hydrodynamics”. In: Physics of Fluids doi : .Bridges, Thomas, Jonathan Pennant, and Sergey Zelik (July 2014). “DegenerateHyperbolic Conservation Laws with Dissipation: Reduction to and Validity ofa Class of Burgers-Type Equations”. In: Archive for Rational Mechanics andAnalysis doi : .Brown, Donald, Peter Popov, and Yalchin Efendiev (2011). “On homogenizationof stokes flow in slowly varying media with applications to fluid–structureinteraction”. In: GEM - International Journal on Geomathematics , pp. 1–25. issn : 1869-2672. doi : (cit. on p. 26).Carr, J. (1981). Applications of centre manifold theory . Vol. 35. Applied Math.Sci. Springer–Verlag. url : http://books.google.com.au/books?id=93BdN7btysoC (cit. on p. 28).Chen, Chen, A. J. Roberts, and J. E. Bunder (2016). Boundary conditions formacroscale waves in an elastic system with microscale heterogeneity .Tech. rep. [ http://arxiv.org/abs/1603.06686 ] (cit. on p. 4).Chicone, Carmen (2006). Ordinary Differential Equations with Applications .Ed. by J. E. Marsden, L. Sirovich, and S. S. Antman. Vol. 34. Texts inApplied Mathematics. Springer (cit. on p. 19).Craster, Richard V. (2015). “Dynamic Homogenization”. In: ed. byV. V. Mityushev and M. Ruzhansky. Vol. 116. Springer Proceedings inMathematics and Statistics. Springer, pp. 41–50. doi : (cit. on p. 26).Cross, M. C. and P. C. Hohenberg (1993). “Pattern formation outside ofequilibrium”. In: Rev. Mod. Phys. doi : (cit. on p. 2).Dyke, M. van (1987). “Slow Variations In Continuum Mechanics”. In: AdvApplied Mech 25, pp. 1–45 (cit. on p. 2).Engquist, B. and P. E. Souganidis (2008). “Asymptotic and numericalhomogenization”. In: Acta Numerica 17. doi:10.1017/S0962492906360011,pp. 147–190. doi : (cit. on pp. 4, 26).Gallay, Th. (1993). “A center-stable manifold theorem for differential equationsin Banach spaces”. In: Commun. Math. Phys Modulational instability of two pairs of counter-propagating waves and energyexchange in two-component media . Tech. rep.[ http://arXiv.org/abs/nlin.PS/0503047 ] (cit. on p. 2).Gustafsson, Bj¨orn and Jacqueline Mossino (2003). “Non-periodic explicithomogenization and reduction of dimension: the linear case”. In: IMAJournal of Applied Mathematics 68, pp. 269–298. doi : (cit. on p. 26).Haragus, Mariana and Gerard Iooss (2011). Local Bifurcations, CenterManifolds, and Normal Forms in Infinite-Dimensional Dynamical Systems .Springer. doi : (cit. on p. 2).Lall, Sanjay, Petr Krysl, and Jerrold E. Marsden (2003). “Structure-preservingmodel reduction for mechanical systems”. In: Physica D: NonlinearPhenomena doi : (cit. on p. 2). 30amarque, Claude-Henri, Cyril Touz´e, and Olivier Thomas (2012). “An upperbound for validity limits of asymptotic analytical approaches based on normalform theory”. In: Nonlinear Dynamics , pp. 1–19. doi : (cit. on p. 19).LeVeque, Randall J., David L. George, and Marsha J. Berger (2011). “Tsunamimodelling with adaptively refined finite volume methods”. In: Acta Numerica 20. doi:10.1017/S0962492911000043, pp. 211–289. doi : (cit. on p. 2).MacCullum, M. and F. Wright (1991). Algebraic computing with REDUCE .Oxford Science Pub (cit. on p. 35).Mercer, G. N. and A. J. Roberts (1990). “A centre manifold description ofcontaminant dispersion in channels with varying flow properties”. In: SIAMJ. Appl. Math. 50, pp. 1547–1565. doi : .Mielke, A. (1986). “A Reduction Principle For Non-autonomous Systems InInfinite Dimensional Spaces”. In: J. Diff Equat 65, pp. 68–88.— (1988). “On Saint-venant’s Problem For An Elastic Strip”. In: Proc Roy SocEdin A A. angew MathPhys doi : (cit. on p. 29).Naghdi, P. M. (1972). The Theory Of Plates And Shells . Vol. 6. Handbuch DerPhysik, pp. 425–640 (cit. on p. 2).Nayfeh, A. H. and S. D. Hassan (1971). “The Method Of Multiple Scales AndNonlinear Dispersive Wave”. In: J. Fluid Mech. 48, p. 463 (cit. on p. 2).Newell, A. C. and J. A. Whitehead (1969). “Finite amplitude, finite bandwidthconvection”. In: J. Fluid Mech. 38, pp. 279–303 (cit. on p. 2).Niyonzima, Innocent (2014). “Multiscale Finite Element Modeling of NonlinearQuasistatic Electromagnetic Problems”. PhD thesis. Universit´e de Li`ege,Facult´e des Sciences Appliqu´ees (cit. on p. 26).Noakes, C. J., J. R. King, and D. S. Riley (2006). “On the development ofrational approximations incorporating inertial effects in coating and rimmingflows: a multiple-scales approach”. In: Q. J. Mechanics Appl Math doi : (cit. on p. 2).Pavliotis, G. A. and A. M. Stuart (2008). Multiscale methods: averaging andhomogenization . Vol. 53. Texts in Applied Mathematics. Springer (cit. onp. 4).Pawula, R. F. (1967). “Approximation of the Linear Boltzmann Equation by theFokker-Planck Equation”. In: Phys. Rev. doi : (cit. on pp. 3, 11).Potzsche, Christian and Martin Rasmussen (2006). “Taylor Approximation ofIntegral Manifolds”. In: Journal of Dynamics and Differential Equations doi : (cit. on p. 2).Roberts, A. J. (1988). “The application of centre manifold theory to theevolution of systems which vary slowly in space”. In: J. Austral. Math. Soc. B 29, pp. 480–500. doi : (cit. on p. 14).— (1992). “Boundary conditions for approximate differential equations”. In: J. Austral. Math. Soc. B 34. doi:10.1017/S0334270000007384, pp. 54–80. doi : (cit. on pp. 4, 29).31oberts, A. J. (2007). Computer algebra derives normal forms of stochasticdifferential equations . Tech. rep. http://eprints.usq.edu.au/archive/00001873 (cit. on p. 36).— (2008). “Normal form transforms separate slow and fast modes in stochasticdynamical systems”. In: Physica A doi : (cit. on pp. 8, 9, 29, 36).— (2015a). “Macroscale, slowly varying, models emerge from the microscaledynamics in long thin domains”. In: IMA Journal of Applied Mathematics ,pp. 1–27. doi : (cit. on pp. 2, 26, 28, 29).— (2015b). Model emergent dynamics in complex systems . SIAM, Philadelphia. isbn : 9781611973556. url : http://bookstore.siam.org/mm20/ (cit. onp. 11).Robinson, J. C. (1996). “The asymptotic completeness of inertial manifolds”. In: Nonlinearity ,pp. 1325–1340. doi : (cit. on p. 19).Romanazzi, Pietro, Maria Bruna, and David Howey (Aug. 2016). “Thermalhomogenisation of electrical machine windings applying the multiple-scalesmethod”. In: Journal of Heat Transfer . doi : (cit. onpp. 2, 26).Segel, L. A. (1969). “Distant Side Walls Cause Slow Amplitude Modulation OfCellular Convection”. In: J. Fluid Mech 38, pp. 203–224. doi : (cit. on p. 29).Taylor, G. I. (1953). “Dispersion of soluble matter in solvent flowing slowlythrough a tube”. In: Proc. Roy. Soc. Lond. A Partial Differential Equations I . Applied MathematicalSciences. doi:10.1007/978-1-4419-7055-8. Springer (cit. on p. 14).Westra, Mark-Tiele, Doug J. Binks, and Willem van de Water (2003). “Patternsof Faraday waves”. In: J. Fluid Mech. doi : (cit. on p. 2). A Notation This appendix summarises a lot of the symbols used in the general theory. Thethird column in the table of notation gives the case corresponding to the randomwalker/heat exchanger of Section 2.Symbol Meaning Random walker x ∈ X ‘large’ open spatial domain in R M of the plate ( x , y ) ∈ R y ∈ Y space of the ‘cross-section’ of the plate Y = { 0, 1, 2 } U Hilbert space of field values R u ( x , y , t ) function of field values, in C N ( u , u , u ) u ( n ) ( X , y , t ) for | n | < N , the n th derivative in x evaluatedat x = X u ( n ) ( X , x , y , t ) for | n | = N , weighted integral, between X and x , of the n th derivative in x N some chosen order of truncation of the multi-variable Taylor series N = L k maps U → U , operator coefficients of the k th derivative in x in the original pde × N { 0, 1, 2, . . . } , the non-negative integers32ymbol Meaning Random walker k , n , . . . M -dimensional integer multi-indices in N M , | k | (cid:54) N k = ( k , k ) , . . . U ( x , t ) macroscale emergent variables, in R m ; in gen-eral U := (cid:104) Z , u ( x , y , t ) (cid:105) U ( x , y , t ) A n m × m matrix coefficient of the n th spatialderivative of U in the macroscale pde scalars 0, − , , (cid:80) b | k | = a denotes a sum over all indices k ∈ N M satis-fying a (cid:54) | k | (cid:54) b . (cid:80) ∞ | k | = a denotes a sum over all indices k ∈ N M satis-fying | k | (cid:62) a , but the sum truncates as thereare a finite number of non-zero L k . (cid:80) condition denotes a sum over all variable indices thatsatisfy the specified condition. (cid:80) bk = a denotes a sum over indices k such that a (cid:54) k (cid:54) b , that is, a (cid:96) (cid:54) k (cid:96) (cid:54) b (cid:96) for all (cid:96) = 1, 2, . . . , M . r , r n complicated remainder terms for various ex-pressions, functionals of the field u N = (cid:0) N + MM (cid:1) is the number of de s, each in U ,of the dynamics for any locale in X N = ξ ∈ R M generating multinomial variable—is effec-tively a local space variable, x = X + ξ G N space of multinomials in ξ of degree (cid:54) N U N = U ⊗ G N the space of multinomial fields˜ u ( X ,t) the generating multinomial with coeffi-cients u ( n ) , depends implicitly on ξ and y ˜ r generating multinomial of remainder terms, acomplicated functional of field u G operator to give the generating multinomial,at X , of any given field˜ L (cid:80) N | k | = L k ∂ kξ differential operator on multino-mials, corresponds to given pde E c , E s subspaces of U , invariant under L , centreand stable respectvely λ j , v j complete eigenvalues and (generalised) eigen-vectors (in U ) of L eigenvalues0, − − V ∈ U × m = (cid:2) v v · · · v m (cid:3) , m columns are a basisfor centre subspace E c ( 1, 0, 0 ) Z ∈ U × m m columns are a basis for the centre subspaceof the adjoint L † , also (cid:104) Z , V (cid:105) = I m ( 1, 0, 0 ) E Nc ⊂ U N multinomial centre subspace of ˜ L V n ∈ U × m each are components of centre eigenvectorsof ˜ L , derived recursively˜ V n ∈ U × mN m columns are linearly independent centreeigenvectors of ˜ L ,˜ V ∈ U × m N N = [ ˜ V n ] , its m N columns are a basis for E Nc A matrix of N × N blocks of A n and 0 m , asappropriate. the 6 × U n ∈ R m each are m parameters of the multinomialcentre subspace E Nc , m N in total U ( m , n ) = U mn U ∈ G mN corresponding generating multinomial of thecentre subspace parameters U ∈ R m N = (cid:2) U n (cid:3) , all parameters of the multinomialcentre subspace of ˜ L U = (cid:2) U mn (cid:3) E Ns ⊂ U N multinomial stable subspace of ˜ L ˜ W a basis for the multinomial stable sub-space E Ns of ˜ L (formal) S parameters of the multinomial stable sub-space of ˜ L S = (cid:2) U mnj (cid:3) , j = 1, 2 B restriction of multinomial operator ˜ L to E Ns ˜ r c , ˜ r s components of the remainder term ˜ r in thesubspaces E Nc , E Ns respectively34 Computer algebra models the random walker This section lists and describes computer algebra code to analyse the Taylor seriesapproach to the slowly varying modelling of the heat exchanger (1) of Figure 1.The code is written in the free computer algebra package Reduce (MacCullum andWright 1991, e.g.). Analogous code will work for other computer algebra packages.Make the printing appears nicer. on div; on revpri; off allfac; B.1 Transform PDEs to ODEs for Taylor coefficients Define the x and y components of the velocities in each plane. The matrices vx and vy are diagonal with v j = ( vx jj , vy jj ) . Also define the mixing operator lop which describes the rate of the walker changing among the three directions. nn:=2; vx:=mat((1,0,0),(0,-1,0),(0,0,+1)); vy:=mat((1,0,0),(0, 0,0),(0,0,-1)); lop:=mat((-1,1,0),(1,-2,1),(0,1,-1))$ Define the model in terms of the heat flows c j ( x , y , t ) . array resc(3); operator ct; depend ct,t,x,y; for j:=1:3 do resc(j):=-df(ct(j),t)+(for k:=1:3 sum ( lop(j,k)*ct(k)-vx(j,k)*df(ct(k),x)-vy(j,k)*df(ct(k),y) ))$ Map from the c fields to the slow and fast u heat fields u defined by (3) The u ( x , y , t ) field is slow with eigenvalue λ = u ( x , y , t ) and u ( x , y , t ) are fast with eigenvalues λ = − λ = − array resu(2); operator ut; depend ut,t,x,y; ct(1):=(ut(0)+ut(1)+ut(2))$ ct(2):=(ut(0)-2*ut(2))$ ct(3):=(ut(0)-ut(1)+ut(2))$ write resu(0):=(resc(1)+resc(2)+resc(3))/3; write resu(1):=(resc(1)-resc(3))/2; write resu(2):=(resc(1)-2*resc(2)+resc(3))/6; array evl(2); evl(0):=0$ evl(1):=-1$ evl(2):=-3$ Extract coefficient matrices of these modal ode s for later use. matrix ll00(3,3),ll10(3,3),ll01(3,3); for i:=0:2 do ll00(i+1,i+1):=evl(i); for i:=0:2 do for j:=0:2 do ll10(i+1,j+1):=df(resu(i),df(ut(j),x)); for i:=0:2 do for j:=0:2 do ll01(i+1,j+1):=df(resu(i),df(ut(j),y)); Construct the Taylor expansion of the u fields to order N using Lagrange’sremainder theorem with Taylor coefficients u jmn where j = 0, 1, 2, and m , n = 0, 1, . . . , N with | ( m , n ) | (cid:54) N . Each coefficient u jmn is a function of X , Y and t ,but those with | ( m , n ) | = N are also functions of x and y . operator u; depend u,xx,yy,t; for j:=0:2 do for k:=0:nn do depend u(j,k,nn-k), x,y$ for j:=0:2 do ut(j):=(for m:=0:nn sum for n:=0:nn-m sum gives full information about Reduce. u(j,m,n)*(x-xx)^m/factorial(m)*(y-yy)^n/factorial(n))$ Obtain the ode for each Taylor coefficient. array odeu(2,nn,nn); for j:=0:2 do for m:=0:nn do for n:=0:nn-m do write odeu(j,m,n):=sub(y=yy,x=xx,df(df(resu(j),x,m),y,n))$ The uncertain ‘forcing’ terms are ∂ x u ( m , n ) j and ∂ y u ( m , n ) j with order m + n = N .Rename these w ( j , m , n , z ) where the subscript refers to the derivative ∂ z with z = x , y . (Because sub() is still active, we have to replace xx by x and yy by y.) inclw:=0; % =1 to include, =0 to exclude operator w; depend w,tt; for j:=0:2 do for m:=0:nn do for n:=0:nn-m do write odeu(j,m,n):=((odeu(j,m,n) where df(u(~k,~l,~p),~z)=>w(k,l,p,z)*inclw when z neq t) where {xx=>x, yy=>y}); depend tt,t; B.2 Initialise the construction of a transform We want to determine a new set of fields U jmn for which the evolution ˙ U jmn sepa-rates fast and slow fields. Store the transform for u ( m , n ) j in array ux(j,m,n) , andthe right-hand side of the corresponding ode for U ( m , n ) j in array duudt(j,m,n) operator uu; depend uu,xx,yy,t; let { df(uu(~p,~q,~r),t)=>duudt(p,q,r) , u(~p,~q,~r)=>ux(p,q,r) }; array ux(2,nn,nn), duudt(2,nn,nn); As a first approximation the transform u (cid:55)→ U is the identity and ˙ U jmn = λ j U jmn . for m:=0:nn do for n:=0:nn-m do for j:=0:2 do begin ux(j,m,n):=uu(j,m,n); duudt(j,m,n):=evl(j)*uu(j,m,n); end; Need to express the uncertain remainders as history integrals so use well es-tablished operators (Roberts 2008; Roberts 2007, e.g.). operator z; linear z; let { df(z(~f,tt,~mu),t)=>-sign(mu)*f+mu*z(f,tt,mu) , z(1,tt,~mu)=>1/abs(mu) , z(z(~r,tt,~nu),tt,~mu) => (z(r,tt,mu)+z(r,tt,nu))/abs(mu-nu) when (mu*nu<0) , z(z(~r,tt,~nu),tt,~mu) => -sign(mu)*(z(r,tt,mu)-z(r,tt,nu))/(mu-nu) when (mu*nu>0)and(mu neq nu) }; Define an operator to separate out terms in fast stable variables. operator fast; linear fast; let { fast(uu(0,~m,~n),t)=>0 , fast(uu(~j,~m,~n),t)=>uu(j,m,n) when j>0 , fast(w(~a,~b,~c,~d),t)=>0 , fast(z(~a,tt,~b),t)=>0 }; The above properties are critical : they must be correct for the results to becorrect.Determine the effect on the slow manifold of the fast forcing.36 operator slow; linear slow; let { slow(uu(~j,~m,~n),uu)=>uu(j,m,n)/evl(j) }; B.3 Iterate to separate slow-fast coordinates Iterate to obtain the transform u (cid:55)→ U and the evolution ˙ U jmn . The evolutionof the slow fields U mn should only contain slow fields, whereas the evolution ofthe fast fields, U mn and U mn , should only contain fast fields. For j = 1, 2 , allfast fields in the residue of the ode of fast field u jmn are placed in the evolution˙ U jmn and all remaining terms are placed in the transform of u jmn . All fast fieldsin the residue of the ode of slow field u mn are placed in the transform of u mn and all remaining terms are placed in the evolution ˙ U mn . The coupling terms w may have both slow and fast components but are designated slow. for iter:=1:9 do begin ok:=1; lengthRess:={}; for m:=0:nn do for n:=0:nn-m do begin res:=odeu(0,m,n); % slow modes lengthRess:=length(res).lengthRess; ux(0,m,n):=ux(0,m,n)+slow(gd:=fast(res,t),uu); duudt(0,m,n):=duudt(0,m,n)+(res-gd); if res neq 0 then ok:=0; for j:=1:2 do begin% fast modes res:=odeu(j,m,n); lengthRess:=length(res).lengthRess; duudt(j,m,n):=duudt(j,m,n)+(gd:=fast(res,t)); ux(j,m,n):=ux(j,m,n)+z(res-gd,tt,evl(j)); if res neq 0 then ok:=0; end; end; write lengthRess:=lengthRess; showtime; if ok then write iter:=iter+10000; end; Check the iteration converged to the specified order. if not ok then rederr("The iteration failed to converge"); B.4 Write out the transform On completing the iteration, write all transforms and evolutions. for j:=0:2 do for n:=0:nn do for m:=0:n do write ux(j,n-m,m):=ux(j,n-m,m); for j:=0:2 do for n:=0:nn do for m:=0:n do write duudt(j,n-m,m):=duudt(j,n-m,m); B.5 Check the generating multinomial form First find the slow subspace eigenvectors in the multinomial form: they come fromthe various coefficients of each of the uu(0,m,n) amplitudes. array vt(nn,nn); for n:=0:nn do for m:=0:n do begin vt(n-m,m):=tp mat((0,0,0)); for p:=0:nn do for q:=0:nn-p do vt(n-m,m):=vt(n-m,m)+df( tp mat((ux(0,p,q),ux(1,p,q),ux(2,p,q))) ,uu(0,n-m,m))*xi^p*yi^q/factorial(p)/factorial(q); end; Find that ∂ v ( m , n ) /∂ξ = v ( m − n ) and ∂ v ( m , n ) /∂ξ = v ( m , n − ) . This patternmust be useful.Then put these eigen-multinomials into one whole (works for any basis vectorsat all). The second version below uses dx and dy to notionally symbolise differentialoperators in some manner yet to be decided. p:=-1$ factor zz; vtt:=for n:=0:nn sum for m:=0:n sum vt(n-m,m)*zz^(p:=p+1); vtt:=for n:=0:nn sum for m:=0:n sum vt(n-m,m)*dx^(n-m)*dy^m; Now pre-multiply by ˜ L obtained from the basic modal ode s: lltvtt:=ll00*vtt+ll10*df(vtt,xi)+ll01*df(vtt,yi); Is this the same as the following? matrix avtt(3,1); for i:=1:3 do avtt(i,1):= (duudt(0,0,0) where uu(0,~m,~n)=>df(vtt(i,1),xi,m,yi,n)); errorInOps:=avtt-lltvtt; Yes it is! So, post-multiplying by the eigen-matrix is equivalent to premultiplyingby the slow evolution operator. And this happens automatically.Evaluating avtt or lttvtt at ξ = ζ = B.6 Write transform in LaTeX Write out in pretty LaTeX. The definition of the L A TEX command is a bit dodgyas convolutions of convolutions are not printed in the correct order; however,convolutions commute so it does not matter. load_package rlfi; mathstyle math; defindex ux(down,up,up); defindex duudt(down,up,up); defindex uu(down,up,up); defindex w(down,up,up,down); defid uu,name="U"; defid ux,name="u"; defid duudt,name="\dot U"; defid tt,name="}{"; defid w,name="u"; defid z,name="\z"; Change name to get braces, not left-right parentheses. deflist(’((!( !{) (!) !}) ),’name)$ Force all fractions (coded in Reduce as quotient ) to use \frac command so wecan change how it appears. put(’quotient,’laprifn,’prinfrac); glmsmvs.red for later reading. Prepend the ex-pressions with an instruction to write a heading, and surround the heading withanti-math mode to cancel the math environment that rlfi puts in. out "glmsmvs.red"$ for j:=0:2 do for m:=0:nn do for n:=0:m do write "ux(",j,",",m-n,",",n,"):=ux(",j,",",m-n,",",n,");"; for j:=0:2 do for m:=0:nn do for n:=0:m do write "duudt(",j,",",m-n,",",n,"):=duudt(",j,",",m-n,",",n,");"; write "end;"; shut "glmsmvs.red"; Now write the LaTeX: out "glmsmvs.ltx"$ on latex; write "% need to delete commands by hand: /.*;.*// % and delete commas: /,// \newcommand{\z}[2]{e^{\ifnum "; in "glmsmvs.red"$ off latex; shut "glmsmvs.ltx"$ end;end;end;%%%%%%% overkill reduce B.7 Transform back to slowly varying We now treat the slow fields U mn as the Taylor coefficients of some slow fieldfunction U ( x , y , t ) and construct a slow field pde . operator uufun; depend uufun,x,y,t; All field coefficients U mn are functions of X , Y and t , but those field coefficientswith | ( m , n ) | = N are also functions of x and y . for j:=0:2 do for k:=0:nn do depend uu(j,k,nn-k),x,y$ Define the Taylor series expansion of U ( x , y , t ) with coefficients U c mn = U mn . operator uut; depend uut,x,y,t; operator uuc; depend uuc,t; for j:=0:2 do uut(j):=(for m:=0:nn sum for n:=0:(nn-m) sum uuc(j,m,n)*(x-xx)^m/factorial(m)*(y-yy)^n/factorial(n)); Construct a pde for U ( x , y , t ) which, like the original u ( x , y , t ) pde , is firstorder in time. For clarity, we place w x ( u jmn ) and w y ( u jmn ) for j = 1, 2 and allconvolutions in the function ‘force’. resuu:=(-df(uufun(0),t)+(df(uut(0),t)) where {uuc(~k,~m,~n)=>uu(k,m,n)})$ forced:={z0(~f,t)=>0, z1(~f,t)=>0, z2(~f,t)=>0, wx(1,~k,nn-k)=>0, wy(1,~k,nn-k)=>0, wx(2,~k,nn-k)=>0, wy(2,~k,nn-k)=>0}$ force:=resuu-(resuu where forced)$ resuu:=resuu-force$ for m:=0:nn do resuu:=sub( {wy(0,m,nn-m)=df(ux(0,m,nn-m),y) ,wx(0,m,nn-m)=df(ux(0,m,nn-m),x)} ,resuu); U mn with derivatives of U ( x , y , t ) and its Taylorseries: U mn = ∂ mx ∂ ny U ( x , y , t ) − [ ∂ mx ∂ ny U Taylor ( x , y , t ) − U c mn ] . (58) for l:=0:nn do for k:=0:l do begin resuu:=(resuu where {uu(0,l-k,k)=> df(uufun(0),x,l-k,y,k)-(df(uut(0),x,l-k,y,k)-uuc(0,l-k,k))})$ resuu:=(resuu where {uuc(0,~m,~n)=>uu(0,m,n)}); end; In the pde of U ( x , y , t ) the forcing terms are defined by the function ‘force’. factor df; resuu:=resuu+forcing; End the if-statement that chooses whether to execute the code of this appendix. end;end;