Choose interelement coupling to preserve self-adjoint dynamics in multiscale modelling and computation
aa r X i v : . [ n li n . C G ] N ov Choose interelement coupling to preserveself-adjoint dynamics in multiscale modellingand computation
A. J. Roberts ∗ November 13, 2018
Abstract
Consider the macroscale modelling of microscale spatiotemporaldynamics. Here we develop a new approach to ensure coarse scalediscrete models preserve important self-adjoint properties of the finescale dynamics. The first part explores the discretisation of microscalecontinuum dynamics. The second addresses how dynamics on a finelattice are mapped to lattice a factor of two coarser (as in multigrids).Such mapping of discrete lattice dynamics may be iterated to empowerus in future research to explore scale dependent emergent phenomena.The support of dynamical systems, centre manifold, theory ensuresthat the coarse scale modelling applies with a finite spectral gap, in afinite domain, and for all time. The accuracy of the models is limitedby the asymptotic resolution of subgrid coarse scale processes, and iscontrolled by the level of truncation. As given examples demonstrate,the novel feature of the approach developed here is that it ensures thepreservation of important conservation properties of the microscaledynamics.
Contents ∗ School of Mathematical Sciences, University of Adelaide, South Australia 5005. mailto:[email protected] Self-adjoint preserving coupling conditions 53 Centre manifold theory supports discretisation 8
Dynamical systems theory gives new assurances about the quality of finitedifference and finite element models of nonlinear spatiotemporal systems. Ef-forts to construct approximations to the long-term, low-dimensional dynam-ics of dissipative partial differential equations ( pde s), on its inertial mani-fold [43, e.g.], have largely been based upon the global nonlinear Galerkinmethod [29, 14, 25, e.g.], and its variants [19, 15, e.g.]. In contrast, a ‘holisticdiscretisation’ [30] developed further here is based purely upon the local dy-namics on finite elements while maintaining, as do inertial manifolds, fidelitywith the solutions of the original pde .2o generate a macroscale model of a dissipative pde , we divide spaceinto finite elements, then specially crafted coupling conditions empower usto support a macroscale model with centre manifold theory [8, 22, e.g.]. Suchmacroscale models are used for computations and for analytic understandingof the dynamics. Crucially, the theory supports the existence, exponentiallyquick attractiveness, and approximate construction of a slow manifold of thesystem dynamics, both deterministic [30, e.g.] and stochastic [36]. Earlierresearch developed coupling conditions that both had centre manifold sup-port on finite elements and ensured consistency as the element size becamesmall [33, 24, e.g.]. However, albeit effective for many systems, these couplingconditions fail to preserve at each level of approximation any self-adjoint sym-metry in the underlying spatial dynamics [31]. The innovations introduced inSection 2 are interelement coupling conditions that not only engender centremanifold support, Section 3, and assure consistency for vanishing elementsize, but also preserve self-adjoint symmetry in each approximation . For justone example, Section 4.2, as deduced in the model (21), argues that a partic-ular one-dimensional, nonlinear, continuum diffusion is soundly mapped asfollows to nonlinear dynamics on a macroscale grid (spacing h ): ∂u∂t = ∂∂x (cid:20) u ∂u∂x (cid:21) → dU j dt ≈ (cid:0) U + − + U + (cid:1) , (1)where the grid values U j = u ( X j , t ) . Why this discretisation instead ofothers? Because the new coupling conditions as well as having dynamicalsystems support and classic consistency, also preserve in the model the self-adjoint symmetry of material conservation that is present in the original pde .In many science and engineering applications it is essential to preserve theconservative form of the governing equations. I argue that interelement cou-pling rules akin to those of Section 2 will automatically preserve conservativeforms.Most methods for modelling dynamics posit just two time scales: a fastand a slow scale [12, 13, 28, e.g.]. Indeed, Sections 2–4 implicitly separatethe dynamics of dissipative pde s into the ‘uninteresting’ fast subgrid dynam-ics, and the relevant slow dynamics of macroscale evolution resolved by thediscretisation. But many applications possess a wide variety of interestingspace-time scales [5, 12, e.g.]. Recent research developed a methodology withrigorous support for changing the resolved spatial grid scale by just a factorof two [35]. Homogenisation, in Section 7.2, is one example: the derivation of3quation (50) recommends that the evolution of discrete diffusion on a grid,with spatially varying diffusivity, is mapped to a coarser grid as du i dt = κ i − u i − − ( κ i − + κ i + ) u i + κ i + u i + → dU j dt ≈ (cid:8) K j − U j − − ( K j − + K j + ) U j + K j + U j + (cid:9) , (2)where the coarse grid index j = and the coarser scale diffusivity K j ≈ ( κ − + κ − + κ + + κ + ) is an average over the microgrid diffusivi-ties. The mapping of dynamics from a finer grid to a coarser grid, via finiteelements formed from a small number of fine grid points, may then be iter-ated to generate a hierarchy of models across a wide range of spatial scales,with the theory of centre manifolds to support across the whole hierarchy.This approach promises to empower us with great flexibility in modellingcomplex dynamics over multiple scales. Sections 5–7 further develop thismodelling transformation of discrete dynamics on grids by exploring cou-pling conditions which result in the coarse grid dynamics also preserving theself-adjoint symmetries of the fine grid dynamics.Most two scale modelling methods can also be applied over many scales.However, most established methods require each application to be based upona large ‘spectral gap’: a parameter such as ǫ measures the scale separation,and invoking “as ǫ → ” provides the extreme scale separation. In contrast,multigrid iteration for solving linear equations transforms between lengthscales that are different by (usually) a factor of two [7, 32, 6, e.g.]. Anal-ogously, Sections 5–7 explore modelling dynamics on a hierarchy of lengthscales that differ by a factor of two and hence the ‘spectral gap’ is finite, notinfinite as required by popular extant, non-multigrid, methods for modellingdynamics. Section 6 describes how to divide fine grid lattice dynamics intosmall finite elements, develops interelement coupling rules that preserve self-adjoint symmetries, and then establishes the centre manifold support for theresulting coarse grid models. Section 7 outlines three applications includingthe homogenisation (2).The methodology proposed here uses dynamical systems theory to sup-port and construct accurate coarse grid models of fine scale dynamics, bothcontinuum and discrete. I expect that the basis developed here for deter-ministic dynamics in one spatial dimension can be extended to both higherdimensions and stochastic dynamics. For non-dissipative dynamics, althoughthe formal construction of coarse models follows analogously, current theory4 xu ( x, t ) X j − X j − X j X j + X j + ✲ t t t t t U j − U j − U j U j + U j + ❞ ❞ ❞❞ ❞ ❞❞ ❞ ❞ u j − ( x, t ) u j ( x, t ) u j + ( x, t ) Figure 1: to discretise continuum dynamics, bottom, to the grid, top, rewritethe continuum dynamics of u ( x, t ) as the dynamics of u j ( x, t ) on overlappingelements. This figure plots three consecutive elements (blue, black and ma-genta): the j th element stretches from X j − to X j + . Then analysis supportsand generates rules for the evolution of grid values U j ( t ) = u j ( X j , t ) .gives little support for the relevance of the resulting sub-centre slow manifoldmodels [41, 2, e.g.]. This article is confined to one dimensional dissipativedynamics.I also conjecture an implication for the equation-free modelling methodol-ogy of Kevrekidis et al. [20, 23, 17, 40, e.g.]. When coupling sparse patches ofmicrosimulators, previous research established that an analogous dynamicalsystems approach provided coupling of patches with high order accuracy onthe macroscale [37, 38]. For dynamical systems where we want to automat-ically preserve in the macroscale any microscale symmetries, the couplingcondition (5), shown to be required here, suggests that each patch shouldhave a point source at mid-patch in proportion to fluxes extrapolated fromneighbouring patches. Further research will tell. The next three sections explore the dynamics of a continuum field u ( x, t ) in one spatial dimension. Figure 1 shows a part of an L -periodic domain, u ( x + L, t ) = u ( x, t ) , which is divided into m overlapping elements E j = { x | X j − ≤ x ≤ X j + } for m grid points X j . These grid points need not be equallyspaced. Typically u j ( x ) and v j ( x ) denote fields in the j th element E j . Define5he inner product h u, v i = m X j = Z E j uv dx = m X j = (cid:14)Z X − j X j − uv dx + Z X j + X + j uv dx (cid:15) . (3)These integrals in the inner product are split over the two halves of eachelement to emphasise that perhaps unexpected contributions come from thecentral grid point X j . Theorem 1 (self-adjoint coupling)
The operator L in the linear systemof coupled pde s ∂u j ∂t = L u j = ∂∂x (cid:20) f j ( x ) ∂u j ∂x (cid:21) + αg j ( x ) u j , x ∈ E j , (4) is self-adjoint when coupled with the interelement coupling conditions f j ( X + j ) u jx ( X + j ) − f j ( X − j ) u jx ( X − j )+ γ (cid:2) f j + ( X j ) u j + ( X j ) − f j − ( X j ) u j − ( X j ) (cid:3) − γ ′ (cid:2) f j ( X j + ) u jx ( X j + ) − f j ( X j − ) u jx ( X j − ) (cid:3) = (5) u j ( X j ± ) = γ ′ u j ( X j ) + γu j ± ( X j ± ) and u j ( X − j ) = u j ( X + j ) , (6) in which subscript x denotes spatial differentiation. Usually I link the coupling parameters by γ + γ ′ = ; however, thisconstraint is not necessary in this theorem.In many applications the coefficient functions f j and g j do not dependupon the element j . However, in nonlinear systems we may need the reas-surance of the theorem when applied with nonlinear f ( u, u x ) whence f j = f ( u j , u jx ) will be element dependent. Any dependence upon time and otherparameters are suppressed for clarity. Proof:
Undergraduate algebra proves the theorem. Integration by partsgives hL u, v i = X j (cid:14) (cid:2) f j u jx v j − f j u j v jx (cid:3) X − j X j − + (cid:2) f j u jx v j − f j u j v jx (cid:3) X j + X + j + Z E j u j L v j dx (cid:15) = h u, L v i + X j (cid:8) f j ( X − j ) u jx ( X − j ) v j ( X j ) − f j ( X − j ) u j ( X j ) v jx ( X − j )− f j ( X j − ) u jx ( X j − ) (cid:2) γ ′ v j ( X j ) + γv j − ( X j − ) (cid:3) + f j ( X j − ) (cid:2) γ ′ u j ( X j ) + γu j − ( X j − ) (cid:3) v jx ( X j − )+ f j ( X j + ) u jx ( X j + ) (cid:2) γ ′ v j ( X j ) + γv j + ( X j + ) (cid:3) − f j ( X j + ) (cid:2) γ ′ u j ( X j ) + γu j + ( X j + ) (cid:3) v jx ( X j + )− f j ( X + j ) u jx ( X + j ) v j ( X j ) + f j ( X + j ) u j ( X j ) v jx ( X + j ) (cid:9) (upon renumbering u j ± ( X j ± ) by j ′ = j ± becomes) = h u, L v i + X j (cid:8) u j ( X j ) (cid:2) − f j ( X − j ) v jx ( X − j ) + γ ′ f j ( X j − ) v jx ( X j − )+ γf j + ( X j ) v j + ( X j ) − γ ′ f j ( X j + ) v jx ( X j + )− γf j − ( X j ) v j − ( X j ) + f j ( X + j ) v jx ( X + j ) (cid:3) + v j ( X j ) (cid:2) f j ( X − j ) u jx ( X − j ) − γ ′ f j ( X j − ) u jx ( X j − )− γf j + ( X j ) u j + ( X j ) + γ ′ f j ( X j + ) u jx ( X j + )+ γf j − ( X j ) u j − ( X j ) − f j ( X + j ) u j,x ( X + j ) (cid:3) (cid:9) = h u, L v i + by (5).Hence L with coupling conditions (5)–(6) is self-adjoint. This proof also ap-plies, with only minor modifications, to multi-component systems where thefield u j in each element is a vector and f j and g j are symmetric matrices. ♠ Remark
The virtue of such coupling conditions for overlapping domainsis that when finding the adjoint, the algebra involves quantities that arealready involved, namely the grid values. If one tries to find the adjointfor disjoint elements, then the field values at the element boundaries alsoare involved, thus leading to too many free variables after the integrationby parts. The overlap of finite elements used herein are increasingly beingused in multiscale modelling: examples include the ‘border regions’ of theheterogeneous multiscale method [13, e.g.], the ‘buffers’ of the gap-toothscheme [40, e.g.], and the overlapping domain decomposition that improvesconvergence in waveform relaxation of parabolic pde s [16, e.g.].7
Centre manifold theory supports discreti-sation
Based upon the equilibria and spectra about equilibria, centre manifold the-ory rigorously supports the existence, relevance and construction of low di-mensional models of dynamical systems [8, 22, e.g.]. This section exploresthe theoretical support for forming discretisations of the class of self-adjoint pde s in the form ∂u∂t = ∂∂x (cid:20) f ( x, u, u x ) ∂u∂x (cid:21) + αg ( x, u, u x ) , (7)using the self-adjoint preserving coupling conditions (5)–(6) with γ ′ = − γ .Embed the L -periodic dynamics of the pde (7) into the overlapping ele-ments of Figure 1 by introducing the subgrid field u j ( x, t ) in each element E j .Define these subgrid fields to satisfy the pde (7) in each element, namely ∂u j ∂t = ∂∂x (cid:20) f ( x, u j , u jx ) ∂u j ∂x (cid:21) + αg ( x, u j , u jx ) , j =
1, . . . , m , (8)and coupled by the conditions (5)–(6).In the extended state space E = ( u , . . . , u m , γ, α ) there exist a sub-space E of equilibria with parameters γ = α = ( γ ′ = ) and subgridfields constant in each element, u j ( x, t ) = U j ; that is, fields which are piece-wise constant over the domain are equilibria when γ = α = . For eachelement define f j ( x ) = f ( x, U j , 0 ) then the differential equation for pertur-bations u ′ j ( x, t ) , linearised about the piecewise constant equilibria, is simplythe diffusion equation ∂u ′ j ∂t = L u ′ j = ∂∂x (cid:20) f j ( x ) ∂u ′ j ∂x (cid:21) , j =
1, . . . , m , (9)with, as γ = ( γ ′ = ), the ‘insulating’ version of the coupling condi-tions (5)–(6), namely for j =
1, . . . , mf j ( X + j ) u ′ jx ( X + j ) − f j ( X − j ) u ′ jx ( X − j )− f j ( X j + ) u ′ jx ( X j + ) + f j ( X j − ) u ′ jx ( X j − ) = (10) u ′ j ( X j ± ) = u ′ j ( X − j ) = u ′ j ( X + j ) . (11)8he centre manifold support for discrete models of the nonlinear sys-tem (7) with coupled elements rests upon the eigenstructure of the linearisedsystem (9)–(11). As for all self-adjoint systems, following the proof of orthog-onality of eigenmodes, Section 3.2, I proceed to prove that the eigenvalues arereal and they are all non-positive. The proofs are elementary undergraduatealgebra. Then Section 3.3 uses these to prove the existence and relevance ofa slow manifold, ‘holistic’ discretisation for the nonlinear system (7). As one important example, this subsection considers the special case wherethe grid is equispaced and the ‘diffusivity’ f only depends upon space throughits dependence upon the field u , that is, f = f ( u, u x ) . Then the lineariseddynamics about the piecewise constant solutions are described by the pde (9)with coefficients f j which are constant on each element. Assume all f j arebounded above zero. This subsection then finds the spectrum and eigenmodesof the linearised pde (9).Seek eigenvalues λ in the j th element such that f j u ′ xx = λu ′ (droppingsubscript j for simplicity). As the eigenvalues must be real (see Theorem 3),set λ = − f j k for some k ( ≥ ) to be determined. For an equi-spaced grid, ∆X j = h , solutions must be of the form u ′ = A cos k ( x − X j ) + B sin k ( x − X j ) + C sin k | x − X j | , (12)upon using the continuity at x = X j . Then the insulating boundary condi-tions (11) of u ′ ( X j ± ) = u ′ ( X j ) require A ( cos kh − ) + ( C ± B ) sin kh = .The difference of these two conditions is sin kh = and hence either B = or kh = nπ . Explore both in turn. kh = nπ Then cos kh = (− ) n and the two conditions reduce to simply A (cid:2) (− ) n − (cid:3) = for which two cases arise depending upon even orodd values for n . even n Then A , B and C are unrestrained except possibly by thederivative condition in (10). However, all modes independentlysatisfy the derivative condition (10). Hence the three orthogonalmodes u ′ = cos k ( x − X j ) , u ′ = sin k ( x − X j ) and u ′ = sin k | x − X j | form a basis for the three dimensional eigenspace correspondingto eigenvalue λ = − f j k . 9 dd n Then A = . Furthermore, sin k ( x − X j ) satisfies the derivativecondition in (10), but sin k | x − X j | does not, C = . Hence u ′ = sin k ( x − X j ) is the only eigenmode.One may like to view the modes sin k | x − X j | for even n withwavenumber k = nπ/h as serving in place of the usual Fouriermodes cos k ( x − X j ) for odd n . kh = nπ In this case, necessarily B = ; then the condition A ( cos kh − ) + C sin kh = and the derivative condition (10) become (cid:20) cos kh − sin kh − sin kh cos kh − (cid:21) (cid:20) AC (cid:21) = . The determinant ( cos kh − ) + sin kh = ( − cos kh ) is non-zeroexcept for the outlawed cases kh = . Hence, A = C = also, andso there are no further eigenmodes.The spectrum of the linearised dynamics (9)–(11) (and multiplicity) onthe j th element is thus { − f j π /h , − j π /h (triple) , − j π /h , − j π /h (triple) , . . . } . (13)This set, with one zero eigenvalue and the rest negative, provided ‘diffusivity’ f j > 0 , are as required for the support of centre manifold theory. This subsection proves, using elementary undergraduate algebra, that crucialproperties of the spectrum (13) hold for the more general system (9)–(11)for possibly non-uniform elements. Alternatively, one may view the firsttwo of these properties as a consequence of the Self-adjoint Theorem 1, andthe third as a natural consequence of the dissipation of diffusion. A readercomfortable with such a view may proceed direct to the next Section 3.3.Throughout this subsection we only address a generic j th element in iso-lation, the interelement coupling parameter γ = . This isolation is toexplore desirable properties of the dynamics in each isolated element. Thenext Section 3.3 then uses centre manifold theory, based upon these results,to support the modelling of fully coupled dynamics.10 heorem 2 (orthogonal eigenmodes) Consider the operator L on theright-hand side of the pde (9) with boundary conditions (10) – (11) ; eigen-modes corresponding to distinct eigenvalues are orthogonal.Proof: Adapt the classical proof in many undergraduate texts. Let u and v be eigenmodes of L in the j th element corresponding to eigenvalues λ and µ , respectively. For slight simplicity use f to denote f j . Since L u = λu and L v = µv , = h v, L u − λu i − h u, L v − µv i (integrating by parts) = (cid:2) fvu x − fuv x (cid:3) X j + X + j + (cid:2) fvu x − fuv x (cid:3) X − j X j − + ( µ − λ ) h u, v i = f ( X j + ) v ( X j + ) u x ( X j + ) − f ( X j + ) u ( X j + ) v x ( X j + )− f ( X + j ) v ( X + j ) u x ( X + j ) + f ( X + j ) u ( X + j ) v x ( X + j )+ f ( X − j ) v ( X − j ) u x ( X − j ) − f ( X − j ) u ( X − j ) v x ( X − j )− f ( X j − ) v ( X j − ) u x ( X j − ) + f ( X j − ) u ( X j − ) v x ( X j − )+ ( µ − λ ) h u, v i (by (11), u ( X ± j ) is continuous and u ( X j ± ) = u ( X j ) ,and similarly for v ) = v ( X j ) (cid:2) f ( X j + ) u x ( X j + ) − f ( X + j ) u x ( X + j )+ f ( X − j ) u x ( X − j ) − f ( X j − ) u x ( X j − ) (cid:3) + u ( X j ) (cid:2) − f ( X j + ) v x ( X j + ) + f ( X + j ) v x ( X + j )− f ( X − j ) v x ( X − j ) + f ( X j − ) v x ( X j − ) (cid:3) + ( µ − λ ) h u, v i (by the boundary condition (10) on gradients) = ( µ − λ ) h u, v i . Hence, for distinct eigenvalues, λ = µ , the corresponding eigenmodes mustbe orthogonal, h u, v i = . ♠ Theorem 3 (reality)
Consider the operator L in pde (9) with boundaryconditions (10) – (11) ; its eigenvalues are all real.Proof: Again adapt the classic proof. Let u denote any eigenmodecorresponding to any eigenvalue λ on the j th element; they are potentially11omplex. Let v and µ denote the complex conjugate of u and λ , respectively.Since the pde (9) and boundary conditions (10)–(11) have all real coefficients,these complex conjugates must also be an eigenmode/eigenvalue pair. Thederivation within the proof of Theorem 2 then establishes that ( µ − λ ) h u, v i = Here h u, v i = R E j | u | dx which is necessarily non-zero and hence µ = λ . But µ and λ are complex conjugates, so any eigenvalue λ must be real. ♠ Theorem 4 (spectral gap)
Consider the operator L in pde (9) with bound-ary conditions (10) – (11) : when the coefficient function f j ( x ) is bounded abovezero, f j, min = min x ∈ E j f j ( x ) > 0 , then there is a zero eigenvalue and all othereigenvalues are negative and bounded away from zero by an amount propor-tional to f j, min / ( X j + − X j − ) .Proof: Let u denote any eigenmode corresponding to any eigenvalue λ on the j th element. Thus λu = [ f j u x ] x . Multiply this ode by u and integrate: λ Z E j u dx = Z E j u [ f j u x ] x dx = (cid:2) f j uu x (cid:3) X j + X + j + (cid:2) f j uu x (cid:3) X − j X j − − Z E j f j u dx = f j ( X j + ) u ( X j + ) u x ( X j + ) − f j ( X + j ) u ( X + j ) u x ( X + j )+ f j ( X − j ) u ( X − j ) u x ( X − j ) − f j ( X j − ) u ( X j − ) u x ( X j − )− Z E j f j u dx (by (11), u ( X ± j ) is continuous and u ( X j ± ) = u ( X j ) ) = u ( X j ) (cid:2) f j ( X j + ) u x ( X j + ) − f j ( X + j ) u x ( X + j )+ f j ( X − j ) u x ( X − j ) − f j ( X j − ) u x ( X j − ) (cid:3) − Z E j f j u dx (by the boundary condition (10) on gradients) = − Z E j f j u dx . f j bounded above zero, the right-hand sideis non-positive and hence so must all the eigenvalues λ . The eigenvalue zerocorresponds only to solutions that are constant on the j th element, u x = .Now prove that all other eigenvalues are bounded away from zero byusing a bound from the constant coefficient case of Section 3.1. Constrainthe magnitude of the eigenmodes by R E j u dx = . Then use the aboveidentity to bound the eigenvalue λ = − Z E j f j u dx ≤ − f j, min Z E j u dx . (14)Thus minimise R E j u dx subject to R E j u dx = and R E j u dx = (as othereigenmodes are necessarily orthogonal to the constant eigenmode) and theboundary conditions (10)–(11). Using Lagrange multipliers µ and ν , stan-dard Calculus of Variations asserts this minimum occurs for functions u ( x ) satisfying the Euler–Lagrange equation ν + − xx = . For simplicity,let ξ = x − X j and X j ± − X j = ± ¯ h + ˜ h . Three cases arise: µ = + k > 0 ⇒ u = ν2µ + A cosh kξ + B sinh kξ + C sinh k | ξ | ; µ = ⇒ u = νx + A + Bξ + C | ξ | ; µ = − k < 0 ⇒ u = ν2µ + A cos kξ + B sin kξ + C sin k | ξ | . In the first case (hyperbolic), substituting u into the boundary conditions (10)–(11) and the orthogonality condition gives four linear equations for A , B , C and ν which have nontrivial values only when the determinant − ¯ hk sinh k ¯ h (cid:2) cosh k ¯ h − cosh k ˜ h (cid:3) = This has no real solutions for k . In the second case, substituting leads to thedeterminant ¯ h ( ¯ h + ˜ h )( ¯ h − ˜ h ) , which also cannot be zero for non-degenerategrids. In the third case (trigonometric), A , B , C and ν have nontrivial valuesonly when the determinant ¯ hk sin k ¯ h (cid:2) cos k ¯ h − cos k ˜ h (cid:3) = This only has solutions for finite k —for non-degenerate grids the smallest k = π/ ¯ h —hence R E j u dx is bounded away from zero by an amount propor-tional to ¯ h − ∝ ( X j + − X j − ) − , and thus so are the negative eigenvalues. ♠ j th element is qualitatively like that of the constant coefficient case (13). Inparticular, there always exists a spectral gap between the zero and the othereigenvalues. Recall we aim to rigorously support discretisation of the nonlinear dynam-ics of the self-adjoint nonlinear general reaction diffusion equation (7). Fordefiniteness we seek spatially periodic solutions, u ( x + L, t ) = u ( x, t ) forsome period L , and divide the domain into m overlapping elements as shownschematically in Figure 1 (with X m − X = L ). The subgrid fields u j ( x, t ) ineach element satisfy the pde (8) and are coupled by the conditions (5)–(6).Then Theorem 4 proves that linearised about any of the piecewise constantsolutions of the subspace E there is a spectral gap in the dynamics of thereaction diffusion pde (7). Consequently centre manifold theory [8, 22, e.g.]asserts the following. Corollary 5 (slow manifold)
For sufficiently smooth reaction g and diffu-sivity f in some neighbourhood of the subset of E for which f j, min are boundedabove zero:1. there exists a ( m + ) dimensional slow manifold M of the subgrid pde (8) coupled by (5) – (6) , one dimension for each element, and onedimension each for parameters γ and α ;2. the slow manifold M may be parametrised by any reasonable mea-sure U j of the field in each element, that is, the slow manifold and theevolution thereon may be written, for some u j and g j , j =
1, . . . , m , as u j = u j ( U , x, α, γ ) such that ˙ U j = dU j dt = g j ( U , α, γ ) ; (15)
3. the dynamics on M is ‘asymptotically complete’ [39] in that for allsolutions of the subgrid pde (8) coupled by (5) – (6) from initial con-ditions in some neighbourhood of M , there exists a solution of theslow manifold model (15) that is approached exponentially quickly intime—roughly at a rate min j { f j, min / ( X j + − X j − ) } ; . the order of error of an approximation to the slow manifold M and itsevolution, (15) , is the same as the order of the residuals of the governing pde (8) and coupling (5) – (6) when evaluated at the approximation. The evolution (15) on the slow manifold M forms the discrete modelof the dynamics of the pde (8) coupled by (5)–(6). In principle, such amodel is an exact closure for the discretisation in that the model tracks theevolution from general initial conditions [Corollary 5.3]. However, we hardlyever can construct an exact slow manifold. Nonetheless, computer algebrareadily constructs the slow manifold and its evolution to a controllable orderof accuracy [Corollary 5.4]. Then evaluating the model for full coupling, γ = , generates a model for the dynamics of the physical reaction diffusion pde (7). Section 4.3 provides evidence that γ = is within the finite domainof validity of the slow manifold M . This section briefly describes three interesting applications of the precedingCorollary 5. Section 4.1 investigate the simplest case of linear diffusion inorder to show how the slow manifold varies with coupling parameter γ . Theinterest is to see how the approach generates classic interpolation for thesubgrid fields and classic finite difference rules for the evolution. Section 4.2explores nonlinear diffusion to demonstrate that this approach supports spe-cific discretisations for nonlinear problems and that, through using the cou-pling conditions (5)–(6), the discretisation preserves the self-adjointness ofthe original physical system. Lastly, Section 4.3 illustrates one way to ac-count for domain boundary conditions other than periodic, and also verifiesconvergence in the coupling γ in one particular example. The simplest application of the Slow Manifold Corollary 5 is to formingdiscrete models of linear homogeneous diffusion as governed by the pde ∂u∂t = ∂ u∂x . (16)As described in Section 3, embed the L -periodic spatial domain into dynamicson m overlapping elements. Corollary 5 guarantees there exists a relevant15iscrete, slow manifold, model of the diffusion (16).Computer algebra [34, §
2] readily constructs approximations to the slowmanifold; one may check the approximation by confirming (17)–(18) satisfiesthe governing equations, (16) and (5)–(6), to the specified order of error. Letthe grid be uniform, ∆X j = h , and define the grid values U j ( t ) = u j ( X j , t ) .Computer algebra finds the subgrid, intraelement field is u j = (cid:2) + γ ( ξµδ + | ξ | δ ) + γ ( ξ δ − | ξ | δ ) (cid:3) U j + O (cid:0) γ (cid:1) . (17)in terms of the subgrid variable ξ = ( x − X j ) /h , and the centred mean µ and difference δ operators, µU j = ( U j + + U j − ) /2 and δU j = U j + − U j − . Reassuringly, when evaluated at full coupling γ = the terms linearand quadratic in the coupling parameter γ form classic linear and quadraticinterpolation, respectively, between the grid values U j . The correspondingevolution on the slow manifold is the discretisation˙ U j = γ δ U j − − γ δ U j + − + γ δ U j − − + − γ δ U j + O (cid:0) γ , δ (cid:1) . (18)Evaluated at full physical coupling γ = , this model recovers the classiccentred finite difference formula for the discretisation. Computing to higherorders in coupling γ , gives more and more terms in the classic formula.This approach recovers classic formula in such simple linear dynamics.However, the theoretical support is different: centre manifold theory appliesat finite element size, to guarantee a relevant model for all initial conditionsin some finite domain, and, after initial transients decay, for all times. Theonly approximation is the error incurred by the truncation of the descriptionof the slow manifold in coupling parameter γ . Nonlinear diffusion has many applications and is of continuing interest [21,45, e.g.] As a specific example application, consider nonlinear diffusion gov-erned by the pde ∂u∂t = ∂∂x (cid:18) u ∂u∂x (cid:19) . (19)16s described in Section 3, embed the L -periodic domain into dynamics on m overlapping elements, with the coupling conditions on the gradient im-plemented with nonlinear diffusivity f ( x, u, u x ) = u . Corollary 5 guaranteesthere exists a relevant discrete, slow manifold, model of the nonlinear diffu-sion (19).The rate of attraction to the slow manifold is no longer a constant: insteadit is now proportional to some measure of the smallest amplitude of theinitial field. Consequently, high order asymptotic approximations to the slowmanifold are replete with divisions by grid values U j . Such divisors reflectthe nonlinear diffusion and suggest that the domain of attraction of the slowmanifold model reduces for fields u of small magnitude.Computer algebra [34, §
3] readily constructs the subgrid field of the slowmanifold: u j = (cid:2) + γξµδ + γ ( − γ ) | ξ | δ ) + γ ξ δ (cid:3) U j + γ (cid:2) ξ ( − | ξ | )( µδU j + j µδU ) + ( | ξ | − ξ )( δ U j − j δ U ) (cid:3) + O (cid:0) γ (cid:1) . (20)The first line is identical to that for linear diffusion, (17); thus the secondline is due to the nonlinearity in the diffusion. Unlike usual finite differencesor finite elements methodology which imposes an interpolation between gridvalues, part of the value of this ‘holistic discretisation’ [33, e.g.] is thatthe subgrid field is constructed to satisfy the pde (19) and thus generatesaccurate closures for the discrete model. The evolution on the slow manifoldgives the discrete model˙ U j = γ δ U − − γ δ U + γ δ U + O (cid:0) γ (cid:1) . (21)Truncated to errors O (cid:0) γ (cid:1) then evaluated at full coupling, γ = , this isthe introductory model (1). To errors of order O (cid:0) γ (cid:1) , the discretisation (21)of the nonlinear diffusion is the same as the discretisation (18) of lineardiffusion, but applied to U instead of to U j . This is reasonable since thiscontinuum nonlinear diffusion operator ( uu x ) x = ( u ) xx . This nontrivialcorrespondence, not imposed at all but instead a natural closure from theintricate subgrid scale interactions, confirms this approach to discretisationis sound.However, such a very close correspondence breaks down at O (cid:0) γ (cid:1) whendivisions by U j ± start invading the slow evolution (21) of the nonlinear17iffusion. Subtleties in such higher order subgrid scale interactions result insuch more complicated discretisations.Nonetheless, the equivalent differential equation of (21) is ∂U∂t = γ ∂∂x (cid:18) U ∂U∂x (cid:19) + h γ ( − γ ) − (cid:18) U∂x + U ∂ U∂x (cid:19) + O (cid:0) h , γ (cid:1) , (22)which is in conservative form no matter what order we truncate the anal-ysis in coupling parameter γ . The coupling conditions (5)–(6) do preserveconservation. In every other section we explore dynamics far away from physical boundariesby assuming spatial periodicity. Conversely, this section explores the extremecase of precisely one element between physical boundaries forming that oneelement, say the element is non-dimenisonalised to − . This ex-treme case empowers us to compute to high order and show convergence inthe ‘coupling’ parameter γ , as well as illustrating the ease of incorporatingphysical boundary conditions on the global domain rather than assumingperiodic conditions for the m elements as used elsewhere.For definiteness, suppose the physical boundary conditions for the exam-ple pde of nonlinear diffusion (19) are the Dirichlet conditions that u = at x = ± . Implement these Dirichlet conditions on the one element [−
1, 1 ] via the adaptation of the coupling conditions (5)–(6) to u ( + , t ) u x ( + , t ) − u ( − , t ) u x ( − , t )= γ ′ (cid:2) u (
1, t ) u x (
1, t ) − u (−
1, t ) u x (−
1, t ) (cid:3) , (23) u ( ±
1, t ) = γ ′ u (
0, t ) and u ( − , t ) = u ( + , t ) . (24)When coupling parameter γ = ( γ ′ = ) the element is isolated fromthe physical boundaries and the spectrum of the dynamics on the one el-ement are as for the previous Section 4.2. Thus there exists a slow manifoldparametrised by γ and the mid-element value U ( t ) = u (
0, t ) : the slow man-ifold may be described by u ( x, t ) = u ( x, U , γ ) such that ˙ U = g ( U , γ ) . Itis exponentially quickly attractive, roughly at a rate ∝ U .18traightforward adaptions of the computer algebra [34] construct the slowmanifold and the evolution thereon to high order in the coupling parameter γ .The slow manifold is u = U (cid:2) − γ | x | + γ ( | x | − x ) + γ ( | x | + − | x | ) + O (cid:0) γ (cid:1)(cid:3) . (25)Observe the γ and γ subgrid structures give classic linear and quadraticinterpolation when evaluated at the physical γ = . It is the O (cid:0) γ (cid:1) termsthat begins to account for the nonlinear nature of the diffusion, hence pro-vide correct subgrid structures, and consequently provide a sound closure forthe macroscale, one dimensional model on the element. The correspondingevolution on the slow manifold is˙ U = − U (cid:2) γ + γ − γ − γ − γ + γ + γ + O (cid:0) γ (cid:1)(cid:3) . (26)This predicts that from all nearby initial conditions, and apart from expo-nentially quick transients, solutions decay algebraically in proportion to as ‘material’ u diffuses through the physical boundaries at x = ± where thediffusivity is zero but the flux u ∂u∂x is not.With just one ‘amplitude’ U computer algebra generates approxima-tions to O (cid:0) γ (cid:1) within a minute. Then generalised Domb–Sykes plots [26,Appendix] estimate the location of the convergence limiting singularity. Fig-ure 2 strongly suggests convergence in coupling parameter γ with radius ofconvergence due to a complex conjugate pair of (logarithmic) singulari-ties in complex γ at an angle of ◦ to the real γ axis. These plots providegood evidence that evaluation at γ = of the power series’ is convergent forat least this simple case. Research into spatio-temporal dynamics commonly invokes a discrete lat-tice [9, 10, 27, 42, 18, e.g.]. Section 6 introduces rigorous support for atransformation from fine scale lattice dynamics to a coarser scale lattice dy-namics. Similar to the continuum case of Sections 2–4, the transformationrequires embedding the dynamics in a system of twice the dimensionality,before reducing the dimension by a factor of four, to achieve an overall halv-ing of the dimensionality. The resultant model then resolves dynamics on alattice with twice the grid spacing. This section argues that it is difficult,19 / r (a) c o s a n g l e (b) Figure 2: power series to errors O (cid:0) γ (cid:1) generate these generalised Domb–Sykes plots [26, Appendix] that, when extrapolated to = , strongly sug-gest the power series in coupling parameter γ has radius of convergence .20 t t ❞❞ t t ❞ u u u u u u v v v v v v v v Figure 3: rewrite the dynamics of u as the dynamics of v on two elements byrenaming variables: the solid discs correspond to differential equations, andthe open circles correspond to algebraic coupling equations.if not impossible, to achieve such net halving of the dimensionality withoutthe initial doubling via the embedding.As the most basic, but key, example of self-adjoint dynamics on a lattice,consider the discrete diffusion equation nondimensionalised to˙ u i = u i − − i + u i + , (27)on equi-spaced grid points x i = ih . This section seeks to construct a soundcoarse scale model for these dynamics, but fails because of interesting reasonsthat inspire the next sections.For almost extreme simplicity, suppose the discrete diffusion (27) applieson a small finite domain at just four lattice points, i =
1, 2, 3, 4 , with insu-lating Neumann-like boundary conditions provided by u − u = u − u = (28)We seek a model in just two dynamical variables, for this system with fourdynamical variables. Because the dynamics are so low dimensional, andlinear, we reasonably explore all options.Divide the domain into two elements, the first containing u and u andthe second containing u and u . As shown in Figure 3, to connect thetwo elements introduce new names v = v = u and v = v = u forthe middle two dynamical variables. For convenience, rename the othersvariables v = u , v = u , v = u and v = u . Then write the discretediffusion (27) and insulating boundary conditions (28), together with theinterelement coupling identities, as the differential-algebraic system D ˙ v = L v , (29)21here v = ( v , . . . , v ) , D = diag (
0, 1, 1, 0, 0, 1, 1, 0 ) and the matrix L = − − − − − − − − , (30)where I explain the import of the italic entries shortly. The spectrum of thelinear discrete diffusion (29), from det ( L − λD ) , is { − + √ − − − √ } .We construct a long term model from the two slow modes corresponding tothe two eigenvalues nearest zero. Thus we seek a model with just two dy-namical variables that systematically track the amplitude of the two slowestmodes.Centre manifold theory provides rigorous support for such low dimen-sional modelling [8, 22, e.g.]. But, analogously to the continuum analysis ofSections 2–4, we need to artificially modify the linear operator L to have twoeigenvalues of zero, then implement a homotopy (a smooth path) to recover L and the original dynamics [33, 35, e.g.]. I argue this is impossible, and hencewe need the more complicated embedding of the next section.To see what freedom we have available, first identify the aspects of L thatcannot be changed. We posit that the evolution equations, corresponding tothe second, third, sixth and seventh lines of L cannot be changed as they areto encode the microscale dynamics (27): if we need to modify the microscaledynamics, then any embeddings we find are almost certainly problem specificand thus not of general power. To maintain the self-adjoint symmetry of L ,this then fixes L , L , L and L to be one, and many other entries tobe zero. The entries to vary are those in italics in (30). Second, of these, set L = L = to avoid excessively nonlocal equations. Third, preservingthe zero eigenvalue of u i = constant, the conservation mode, we need eachrow sum to remain zero. Lastly, as well as self-adjoint symmetry, we requireisotropy, left-right symmetry. These four requirements result in three degrees22f freedom spanned by the three matrices Y = − − ,Y = − − −
10 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 − ,Y = − − −
10 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 1 − . Thus the (nearly) most general self-adjoint, isotropic, conservative, matrix ofthe discrete diffusion dynamics is L + y Y + y Y + y Y for some constants y , y and y . Its spectrum is given by the zeros of the characteristic polynomialdet ( L + y Y + y Y + y Y − λD ) = c λ + c λ + c λ + c λ , (31)for some coefficients c k depending upon y , y and y : for example, c = −( − )( + + y − + y ) . Recall that we seek to create a homotopy from dynamics with a double zeroeigenvalue on two decoupled elements, y = and y = (the ‘green line’ in23 −1.0−0.6−0.20.20.61.0 y −1.0−0.6−0.20.2 0.6 1.0y1−1.0−0.6−0.20.20.61.0 y2 Figure 4: red: iso-surface of ‘infinite’ eigenvalues from the coefficient c = of the characteristic polynomial (31). Blue blob: the origin indicatingthe original diffusion dynamics (30). Green line: two decoupled elements.There appears no route from green to blue without encountering ‘infinite’eigenvalues, red. 24igure 4), to the original dynamics with y = y = y = (the ‘blue blob’in Figure 4). Figure 4 indicates that it is impossible to create a homotopyfrom the two decoupled elements to the original dynamics without crossingthe red surfaces. The red surfaces in Figure 4 are the surfaces of c = in the charateristic polynomial (31), and hence represent neighbourhoodswhere eigenvalues become infinitely large, both positive and negative. Wecannot create a smooth homotopy from two decoupled elements to the cou-pled original dynamics through such neighbourhoods. This failure with justtwo elements suggests that, within the class of self-adjoint, isotropic, diffusiondynamics, we cannot artificially divide a domain into disjoint elements.Thus the next section proceeds to explore embedding lattice dynamicsonto overlapping elements analogous to the overlapping elements used for thecontinuum dynamics of Sections 2–4 and for other multiscale approaches [13,40, 16, e.g.]. This section explores how to transform the discrete dynamics of variables u i ( t ) on a fine grid of spacing h , into discrete dynamics of variables U j ( t ) on acoarser grid of spacing H = .Figure 5 schematically shows that the j th element stretches from a neigh-bourhood of X j − to a neighbourhood of X j + . Figure 5 also shows eachelement is divided into two halves, and the variables duplicated so that, no-tionally, u = v j − = v j,2 = v j,4 = v j + and u + = v j − = v j,3 = v j,5 = v j + . Embed the fine grid dynamics in these overlapping elements in a spaceof double the dimensionality by treating v j − and v j,2 , and v j,5 and v j + asindependently evolving variables.In the embedding space we use the inner product h v , w i = P j,i v j,i w j,i . Analogous to the continuum dynamics of Sections 2–4, we need to find rulesto couple neighbouring overlapping elements that are sufficiently generic thatthey are useful for modelling a wide variety of self-adjoint lattice dynamics.Here we focus on the basic case of the coarse grid modelling of the discrete25 t t t t t t t t t t u − u − u − u − u − u + u + u + u + u + ✲ t t t t t U j − U j − U j U j + U j + ❞ t t ❞❞ t t ❞❞ t t ❞❞ t t ❞❞ t t ❞❞ t t ❞ j th element (cid:12) ( j + ) th element (cid:12) (cid:13) ( j − ) th element v j,0 v j,1 v j,2 v j,3 v j,4 v j,5 v j,6 v j,7 v j − v j − v j − v j − v j − v j − v j − v j − v j + v j + v j + v j + v j + v j + v j + v j + Figure 5: to transform dynamics from the fine grid, bottom, to the coarsegrid, top, rewrite the dynamics of u i as the dynamics of v j,i on overlappingelements, here see three consecutive elements (blue, black and magenta), byduplicating and renaming variables: the solid discs correspond to differentialequations showing the number of dynamic variables is doubled, and the circlescorrespond to algebraic coupling equations. Each element has two halves,also shown separated for clarity.diffusion dynamics (27): equation (27) applies at the internal dynamic vari-ables v j,1 , v j,2 , v j,5 and v j,6 of each element, the discs in Figure 5. Now takeup the challenge of using the variables v j,0 , v j,4 , v j,5 and v j,7 , the open circlesin Figure 5, to couple together not only the elements but also the two halvesof each element.Again analogous to the continuum dynamics of Sections 2–4, where rig-orous theorems are invoked we require that the lattice dynamics are periodicin the fine grid with period in i , and thus periodic on the coarse gridwith period m in j .Let us identify the possible domain of evolution and coupling rules. Fol-lowing Section 5, one evolution rule within each element is D ˙ v j = L v j (32)where v j = ( v j,0 , . . . , v j,7 ) , D = diag (
0, 1, 1, 0, 0, 1, 1, 0 ) and equation (30)gives the matrix L . As in Section 5, to preserve the self-adjoint isotropicconservation of the dynamics the upright entries in the matrix L are fixed; wecan only modify the italic entries. As in Section 5, one set of possible changes26o L is spanned by Y , Y and Y . This set changes only internal interactions.There are also three basic matrices that couple elements together preservingself-adjoint isotropic conservation, with the additional constraint that thecoupling has to be ‘local’ on the grid . Let ¯ ε ± denote the shift operators fromone element to its neighbours: define ¯ ε ± v j,i = v j ± . Then write three basismatrices for nearest neighbour coupling changes to L as Y = − ¯ ε − − ¯ ε + ,Y = − ¯ ε − ¯ ε − ¯ ε + ¯ ε + − ,Y = − ¯ ε − − ¯ ε + − ¯ ε − − ¯ ε + . For example, the non-zero elements of Y v j are v j,0 − v j − and v j,7 − v j + which connects end values of neighbouring overlapping elements. Thus weexplore the self-adjoint isotropic conservative dynamics of D ˙ v j = L + X l = y l Y l ! v j , (33)27or some coefficients y l that we are free to choose.When fully coupled, the element dynamics (33) must reduce to that ofthe discrete diffusion (27). These are most compactly expressed in termsof the fine grid shift operator ε ± defined as ε ± v j,i = v j,i ± . (When usedappropriately, the coarse grid shift ¯ ε ± = ε ± .) Then the fine grid diffusion,expressed as ˙ u i = ( ε + − + ε − ) u i , has ‘eigenvalue’ ( ε + − + ε − ) correspond-ing to ‘eigenvector’ u = ( . . . , ε − , ε − , 1, ε + , ε + , . . . ) where all are interpretedin an operator sense. For the element dynamics (33) to reproduce these dy-namics exactly, the dynamics must have the same ‘eigenvalue’ ( ε + − + ε − ) but corresponding to the ‘eigenvector’ v j = ¯ ε j + (
1, ε + , ε + , ε + , ε + , ε + , ε + , ε + ) .Substituting these into (33), equating coefficients of ε ± in all componentsgives a set of equations that uniquely determine y = y = y = and y = y = y = . Thus the unique operator that gives the correct evo-lution of the diffusion (27) with the correct subgrid microstructure for thediffusion is L = − − ¯ ε − ¯ ε − − − − ¯ ε + − ¯ ε − ¯ ε + − − ¯ ε − − − ¯ ε + − ¯ ε + − . (34)In a homotopy from a useful base for centre manifold theory, we must end thehomotopy at this operator for the fully coupled dynamics on the elements. Now seek a base of decoupled elements from which a slow manifold may beconstructed to model the coarse grid dynamics. Perhaps the closest operatorto the fully coupled L of (34) is simply to change the inter-element coupling28hift operators ¯ ε ± to be ones: that is, define L = L + Y + Y = − − − − − − − −
10 0 0 0 1 − − − − . (35)Then straightforward algebra derives that the spectrum of D ˙ v j = L v j is { − − − } for each of the m decoupled elements. The zero eigenvaluewith all the rest negative implies, Section 6.3, that there exists a relevantslow manifold model [8, 22, e.g.] which we can construct with one dimensionfor each element—the dynamics on the m -dimensional slow manifold formsthe coarse grid model.But is L a good choice? and, is it the only choice? The spectrum of thehomotopy answers. Create a general homotopy from decoupled elements tofully coupled by defining the convex combination L γ = ( − γ ) [ L + y Y + y Y + y Y ] + γL , where parameter γ morphs the operator from the decoupled case, γ = , tothe fully coupled case, γ = . The characteristic polynomial, det ( L γ − λD ) ,is too hideous to record here, but computer algebra derives it easily. Thencomputer algebra iteration finds asymptotic approximations to the eigenvaluenear zero: λ = γ ( ¯ ε − − + ¯ ε + ) (cid:2) + ( y − − ) (cid:3) + O (cid:0) γ + | y | (cid:1) . The factor ( ¯ ε − − + ¯ ε + ) , although unexpectedly involving double shiftsover the coarse grid, is exactly what we want as to leading order it is thecentred second difference δ of the fine grid diffusion dynamics (27). However,the term linear in y l ruins this identity and so we choose y = + . Similar computer algebra to higher order in | y | indicates we also need y = − , and analysis to higher order in coupling parameter γ indicatesthat then we need y = . To best match the dynamics of the discretediffusion (27) we choose y = y = y = . Consequently, accuracy of29he coarse grid model for just simple diffusion requires us to use the convexcombination operator L γ = γ ′ L + γL , where normally γ ′ = − γ , (36)as the homotopy to smoothly connect the decoupled matrix (35) with thefully coupled operator (34). Inhomogeneity requires generalisation
Linearisation in Section 6.3 leadsus to consider the class of inhomogeneous, discrete, reaction-diffusion equa-tions on the fine grid of˙ u i = δ ( f i δu i ) + αg i u i = f i − u i − − ( f i − + f i + ) u i + f i + u i + − αg i u i , (37)for some spatially varying ‘reactions’ g i and ‘diffusivities’ f i ± governing thedynamic exchange between u i and u i ± . Equation (37) describes relativelygeneral self-adjoint dynamics on the fine grid.Embed the fine grid dynamics (37) into the dynamics on the finite ele-ments shown in Figure 5 and generalising the self-adjoint, consistent opera-tor (36) by considering the dynamics D ˙ v j = F j v j + α g j , (38)where the linear ‘reaction’ g j = (cid:0)
0, g − v j,1 , g − v j,2 , 0, 0, g + v j,5 , g + v j,6 , 0 (cid:1) , and where the diffusivity matrix F j = (cid:20) F F F F (cid:21) for the four sub-blocks (39) F = − f − + f − −( γ ′ + γ ¯ ε − ) f + f − − f − − f − + f − + f − − f − − f + f − f ( γ ′ + γ ¯ ε + ) + f − f ,F = − f + f − f ( γ ′ + γ ¯ ε − )+ f − f − f + + f + + f + − f + − f + + f + −( γ ′ + γ ¯ ε + ) f + f + − f + , = diag [( γ ′ + γ ¯ ε − ) f , 0, 0, f ( γ ′ + γ ¯ ε − )] ,F = diag [ f ( γ ′ + γ ¯ ε + ) , 0, 0, ( γ ′ + γ ¯ ε + ) f ] , recalling that ¯ ε ± f v j,i = f ± v j ± . The operator on the right-hand sideof (38) is self-adjoint: the only subtlety arises via the operators couplingelements. For example, suppose the operator ¯ ε ± f occurs at element ( i, k ) in F j : then the inner product h w , F v i has the components X j w j,i ¯ ε ± ( f v j,k ) = X j w j,i f ± v j ± = X j w j ∓ f v j,k = X j f ( ¯ ε ∓ w j,i ) v j,k which suitably corresponds to the ( k, i ) elements of F j in the inner prod-uct h F w , v i . Thus the dynamics of the embedded, coupled, finite element,system (38) preserves the self-adjointness in the fine grid inhomogeneousdynamics (37) in a manner consistent with the required homotopy (36) ofhomogeneous dynamics. We have arrived at a separation of the self-adjoint, isotropic, diffusion dy-namics (27) into elements, drawn schematically in Figure 5, with couplingbetween the elements that ranges from decoupled to fully coupled, and pre-serves the self-adjoint nature of the linear dynamics. We now invoke centremanifold theory [8, 22, e.g.] to rigorously support coarse grid modelling of nonlinear lattice dynamics. The nonlinear dynamics are introduced, embed-ded into elements, and then linearised to connect to the preceding results.In analogy with the continuum dynamics of (7), here we consider the classof fine grid dynamics expressible by the general nonlinear, local interaction,lattice rule ˙ u i = δ (cid:2) f ( i, µu i , δu i ) δu i (cid:3) + αg ( i, u i , µδu i ) , (40)for sufficiently smooth ‘diffusivities’ f and ‘reactions’ g . For definite theoret-ical statements, suppose the fine grid lattice and the dynamics (40) on it are -periodic in i . 31mbed the fine grid, -periodic, nonlinear reaction-diffusion dynam-ics (40) into the higher dimensional dynamics of v j within the corresponding m elements of Figure 5 by a nonlinear version of (38). To help write the mapfrom the fine grid to the overlapping elements let i ′ = ( i − ) − sign ( i − ) ,then the nonlinear fine grid dynamics (40) in each element are˙ v j,i = δ (cid:2) f + i ′ δv j,i (cid:3) + αg + i ′ , i =
1, 2, 5, 6, (41)where f + i ′ = f ( + i ′ , µv j,i , δv j,i ) and g + i ′ = g ( + i ′ , v j,i , µδv j,i ) . Controlthe coupling of these to neighbouring elements by the conditions f ± δv j,i + ( γ ′ + γ ¯ ε ± ) (cid:2) f δv j, (cid:3) = = (cid:10) , (42) f ( γ ′ + γ ¯ ε + ) v j,0 + f δv j,i − f ( γ ′ + γ ¯ ε − ) v j,7 = = , . (43)The embedded element dynamics (41)–(43) then has a useful m dimensionalsubspace E of equilibria: α = γ = and v j,i = U j constant in each of the m elements.Linearise the dynamics about each equilibria in E to obtain (38) withoutreaction, α = , and with diffusivites f + i = f ( + i, U j , 0 ) . As the couplingparameter γ = , each element is isolated. When the diffusivities f i areconstant, each element has spectrum proportional to { − − − } —thezero eigenvalue corresponds to v j,i being constant in each element. Elemen-tary algebra shows that provided for all jf j > 0 and ( f − + f + ) − f − f + > 0 , (44)then the spectrum remains as one zero eigenvalue with the other three beingnegative. Consequently, centre manifold theory assures us of the followingcorollary [8, 22, e.g.]. Corollary 6 (slow manifold)
Provided (44) , in some finite neighbourhoodof the subspace E :1. there exists a ( m + ) dimensional slow manifold M of the nonlin-ear, fine grid, element dynamics (41) – (43) , one dimension for eachelement, one for the inter-element coupling parameter γ , and one forthe amplitude parameter α of the reaction; . the slow manifold M may be parametrised by any reasonable mea-sure U j of v j,i in each element, that is, the slow manifold and the evo-lution thereon may be written, for some v j and g j , j =
1, . . . , m , as v j = v j ( U , γ, α ) such that ˙ U j = dU j dt = g j ( U , γ, α ) ; (45)
3. the dynamics on M is ‘asymptotically complete’ [39] in that from allinitial conditions in some neighbourhood of M , there exists a solu-tion of (45) approached exponentially quickly in time by the solutionof (41) – (43) ;4. the order of error of an approximation to the slow manifold M and its evolution, (45) , is the same as the order of the residuals of the gov-erning dynamics (41) – (43) , in the coupling parameter γ and reactionparameter α , when evaluated at the approximation. Basic example of linear diffusion
The simplest example of the classof lattice dynamics to which Corollary 6 applies is the linear discrete dif-fusion (27). Computer algebra [34, §
4] readily iterates to asymptoticallyapproximate the slow manifold model. Here we choose to parametrise theslow manifold of coarse scale dynamics in terms of the coarse variables U j = ( v j,2 + v j,3 + v j,4 + v j,5 ) , (46)which Figure 5 shows to be estimates of the mid-element values of the finegrid variables. Executing the computer algebra deduces that the fine grid,intraelement structure is v j = h (
1, 1, 1, 1, 1, 1, 1, 1 )+ γ4 (− − −
1, 1, −
1, 1, 3, 5 ) ¯ µ ¯ δ + γ (
3, 1, 0, 0, 0, 0, 1, 3 )( ¯ δ + ¯ δ )+ γ (
13, 7, 1, −
5, 5, − − − ) ¯ µ ¯ δ i U j + O (cid:0) γ (cid:1) , in terms of the centred mean ¯ µ and difference ¯ δ operators on the coarse grid.The first line gives the piecewise constant basis for the slow subspace E
33f equilibria. The second line gives a linear variation between neighbouringelements, the third line quadratic, and so on. Corollary 6 assures us that nomatter what the initial conditions for the fine grid diffusion (27), exponen-tially quickly the system will settle onto the slow manifold with the abovestructure on the fine grid.The corresponding evolution of the coarse grid variables U j is˙ U j = h γ (cid:0) ¯ δ + ¯ δ (cid:1) + γ (cid:0) − ¯ δ − ¯ δ (cid:1) + γ (cid:0) ¯ δ + ¯ δ + ¯ δ (cid:1) i U j + O (cid:0) γ (cid:1) . (47)Obtain the lowest accuracy model from the first line, by neglecting terms O (cid:0) γ (cid:1) , and then evaluating at the physical fully coupled case γ = : themodel is ˙ U j = ( U j − − j + U j + ) which is appropriate although sur-prisingly only involves every second coarse grid value. This ‘surprise’ wasforecast in Section 6.2 by the leading operator eigenvalue of the operator L γ .Higher orders in coupling parameter γ give coarse grid models with a widerstencil, and of more accuracy. For example, we see the accuracy through the equivalent fine grid expression of the coarse grid model (47) obtained by theoperator identity that ¯ δ = + δ : in terms of fine grid differences δ , thecoarse model (47) is˙ U j = (cid:2) γ δ + γ ( − γ )( − ) δ (cid:3) U j + O (cid:0) δ , γ (cid:1) . This equivalent fine grid expression demonstrates that when evaluated forfull coupling, γ = , the fourth order differences vanish to leave the correctequivalent model ˙ U j ≈ δ U j . Computer algebra [34, §
4] to high order incoupling γ confirms that higher order differences in the equivalent fine gridexpression similarly vanish. Thus the coarse grid model (47) is an accurateclosure for the coarse scales of the fine grid dynamics.As introduced in earlier work with non-self-adjoint coupling [35], sucha mapping of lattice dynamics from fine grid scale to the coarser scale canbe iterated and renormalised to cover step-by-step the wide range of lengthand time scales on a multigrid [7, 44, e.g.]. Such step-by-step dynamicaltransformations could empower us, in future research, to carefully explorethe development of emergent phenomena in nonlinear and stochastic systemsover multigrids. 34 Three further applications indicates range
The previous section concluded by modelling the dynamics of discrete dif-fusion. By itself, discrete diffusion is well understood. The value of thepreceding section is that it empowers us to model more complicated dynam-ics from the same base with the same powerful centre manifold support. Thissection introduces three interesting applications.
Centre manifold theory, as recorded in Corollary 6, applies to nonlinear dy-namics such as the coarse grid modelling of the discrete reaction-diffusionequation [46, e.g.] ˙ u i = δ u i + α ( u i − u ) , (48)where α ( u i − u ) is some example nonlinear reaction.Simple modifications of earlier code [34, §
5] gives computer algebra thatconstructs the slow manifold model of the coarse grid evolution. Executingthe code derives, for example, the coarse grid model˙ U j = γ ( U j + − j + U j − ) + α ( U j − U )+ αγ (cid:0) − j U + − + − j U j + + j + U j + − + − j U j − + j − U j + − − − j U j − + j − U j − − − (cid:1) + O (cid:0) γ , α (cid:1) . (49)The first line, at full coupling γ = , is a classic model of the reaction-diffusion (48). The second and subsequent lines account for sub-elementreaction and diffusion, interacting together and with neighbouring elements.Such terms are required for an accurate closure of the nonlinear fine griddynamics on the coarse grid. A critical issue in material science is the effective large scale properties of acomposite material with significant microscopic structure. A canonical prob-lem is the effective large scale diffusion through a domain with microscopicvariations in diffusion coefficient [4, 40, 1, e.g.]. Here we transform diffusion35n a fine grid, with fine grid variations in coefficient, into diffusion onto acoarser grid. Here we provide a new and powerful view of the classic resultthat, to leading order, the coarse grid diffusion is a local average of the finegrid diffusion.Consider the fine grid dynamics (40), without reaction α = , where thediffusion coefficient governing flux between fine grid points x i − and x i + is f i = + ǫκ i . The parameter ǫ controls the overall size of the variations inthe diffusivity. For simplicity in construction and interpretation I treat thevariations in diffusivity as small; that is, we construct the slow manifold as apower series in ǫ . Computer algebra then finds the slow manifold as a powerseries in the strength ǫ . Straightforward modifications to earlier computeralgebra [34, §
6] finds, for the example (2) mentioned in the Introduction,˙ U j = γ (cid:8) K j − U j − − ( K j − + K j + ) U j + K j + U j + (cid:9) + O (cid:0) γ (cid:1) , (50)where the coarse grid effective diffusivity coefficients K j = + ǫ (cid:2) κ − + κ − + κ + + κ + (cid:3) + ǫ (cid:2) − ( κ − + κ − − κ + − κ + ) − ( κ − − κ − ) − ( κ + − κ + ) (cid:3) + O (cid:0) ǫ (cid:1) . Because the element coupling preserves self-adjoint symmetry we find thatthe coarse grid, slow manifold model (50) is indeed self-adjoint with theseparticular effective diffusivities.Even more beautiful is that the effective diffusivity on the coarse gridis local : the diffusivity K j ± governing the flux between coarse grid points X j and X j ± depends only upon the fine scale diffusivities between X j and X j ± ,namely between κ and κ ± . This beautiful feature is not built into theapproach, but appears as a natural consequence of this scheme to preserveself-adjoint properties. Zigzag microstructure
The specific zigzag microstructure κ i = (− ) i is straightforward to analyse to higher order [34, § U j = D U j , the coarse grid diffusion operator D = γ ( − ǫ )( + ¯ δ ) ¯ δ (cid:10) − γ ( + ǫ ) ¯ δ γ (cid:2) + + + ( + + ) ¯ δ (cid:3) ¯ δ (cid:11) + O (cid:0) γ , ǫ (cid:1) . (51)As it should, the − ǫ factor shows that the coarse scale diffusion is de-pressed by the microstructure, and indeed vanishes for ǫ = ± reflecting thevanishing of the microscale diffusivity between every second fine scale gridpoint that occurs at ǫ = ± . An outstanding issue in modelling is the direct discrete, macroscale modellingof pattern evolution [11, 3, e.g.]. Currently, in fields such as fluid convectionand in some reaction-diffusion equations, we first model the evolution of rolls,spirals and spots by variants of the Ginzburg–Landau pde derived assumingamplitudes and phases of the structures vary slowly in space. Second, wethen discretise the Ginzburg–Landau pde for numerical simulations of themodulation of the pattern over large scales. The challenge is to analyse theoriginal system dynamics and directly generate such a macroscale discreti-sation in one step . My first attempt to do this showed potential [31], butthe coupling conditions used therein did not preserve the self-adjoint natureof the dynamics and so the discretisation unsatisfactorily did not preserverequired symmetries. Now with self-adjoint coupling conditions we return tomodelling pattern evolution.
Linear dynamics
On a lattice, the simplest microscale pattern is perhapsthe long lasting, zigzag mode of the discrete equation˙ u i = − u i = − u i − − i − u i + . (52)Now adapt the self-adjoint coupling conditions and analysis of Section 6 toderive a model of these zigzag mode dynamics on the coarser grid. Embedthe dynamics on the elements shown in Figure 5 with the evolution rule D ˙ v j = Z γ v j , (53)37here the ‘zigzag’ operator Z γ = − − ′ + γ ¯ ε − γ ′ + γ ¯ ε − − − − − − − ′ + γ ¯ ε + − − ′ + γ ¯ ε − γ ′ + γ ¯ ε + − − ′ + γ ¯ ε − − − − − − −
10 0 0 γ ′ + γ ¯ ε + γ ′ + γ ¯ ε + − − . Then based about the decoupled case of γ = , centre manifold theory [8, 22,e.g.] similarly guarantees the existence and relevance of a slow manifoldparametrised by the amplitude of the zigzag mode local to each element.Computer algebra [34, § (− ) i to account for the zigzag neutral mode, constructs the slowmanifold model. The corresponding evolution of the coarse grid order pa-rameters, the amplitudes U j , is exactly (47) discussed before. The differenceis that here the equation governs the local amplitude of the zigzag mode,rather than the local mean field. Nonlinear amplitude modulation
Even more interesting is nonlinearlattice dynamics which centre manifold theory also supports. Modify thediscrete equation (52) by a local quadratic nonlinearity to˙ u i = − u i + au = − u i − − i − u i + + au . (54)What is the corresponding coarser grid equation governing the local am-plitude of the zigzag mode? In particular, does the quadratic nonlinearitystabilise or destabilise the origin?Simple modifications to the computer algebra then derives the coarserscale, modulation equation˙ U j = γ ( ¯ δ + ¯ δ ) U j − aγ ( ¯ δ ¯ µU j )( ¯ δ U j ) + a U + O (cid:0) γ + a (cid:1) . (55)Evaluated at γ = this coarse grid model predicts that the quadratic nonlin-earity in (54) is destabilising as it generates the cubic growth term a U . Inother problems, the coarse grid model will similarly form a discrete version ofthe Ginzburg–Landau equation. One might derive such a Ginzburg–Landau38quation by traditional modulation theory, but here analysis of the fine scalegrid dynamics generates the novel term aγ ( ¯ δ ¯ µU j )( ¯ δ U j ) which accountsfor interactions over the scale of a few grid points and thus forms a morecomprehensive closure.This analysis empowers us to derive discrete models of nonlinear latticepattern dynamics without having to assume the infinite scale separation re-quired by traditional modulation theory. This article uses centre manifold theory to further develop a novel approachto high quality coarse scale discrete models of nonlinear spatiotemporal sys-tems. The method is to divide space into overlapping finite elements withspecially crafted coupling conditions. The innovative coupling conditions ofSection 2 not only engender centre manifold support, Section 3, and assureconsistency for vanishing element size, but also preserve the self-adjoint sym-metry which is often so important in applications.The first half of this article extracts accurate finite scale discrete latticemodels from analysis of the infinitesimal scale modelling of dissipative pde s.A companion problem is to extract an accurate coarse scale lattice modelfrom a fine scale lattice model. The second half of this article does this forthe possibly extreme case when the coarse scale lattice is just a factor oftwo coarser than the fine scale lattice, and thus connects this approach tomultigrid methods [7]. The approach is again based upon dividing the finelattice into overlapping elements with specially crafted coupling conditions,Section 6. To supplement an earlier approach [35], here the coupling condi-tions preserve self-adjointness in the dynamics as seen in the three exampleapplications of Section 7.Importantly, centre manifold theory [8, 22, e.g.] supports modelling forfinite spectral gap, not just the infinite spectral gap required by other meth-ods. Furthermore, the support applies, in principle, to a finite domain andfor all time; the only approximation is in the approximation of the slow man-ifold. The approach developed here for deterministic dynamics in one spatialdimension should be extendable to both higher dimensions and stochasticdynamics. 39 cknowledgement
I thank the Australian Research Council for supportvia grants DP0774311 and DP0988738.
References [1] Todd Arbogast and Kirsten J. Boyd. Subgrid upscaling and mixed mul-tiscale finite elements.
SIAM J. Numer. Anal. , 44:1150–1171, 2006.doi:10.1137/050631811.[2] O. Bokhove and T. G. Shepherd. On the Hamiltonian balanced dynamicsand the slowest invariant manifold.
J. Atmos Sci , 53(2):276–297, 1996.doi:10.1175/1520-0469(1996)053¡0276:OHBDAT¿2.0.CO;2.[3] H. R. Brand. Phase dynamics—a review and a perspective.
PropagationIn Systems Far From Equilibrium , pages 206–224, 1988.[4] Achi Brandt. General highly accurate algebraic coars-ening.
Elect. Trans. Num. Anal. , 10:1–20, 2000. .[5] Achi Brandt. Multiscale scientific computation: review 2001. InT. F. Chan T. J. Barth and R. Haimes, editors,
Multiscale and Multires-olution Methods: Theory and Applications , pages 1–96. Springer–Verlag,Heidelberg, 2001.[6] Achi Brandt. Methods of systematic upscaling. Technical re-port, Department of Computer Science and Applied Math-ematics, The Weizmann Institute of Science, March 2006. .[7] William L. Briggs, Van Emden Henson, and Steve F. McCormick.
Amultigrid tutorial, second edition . SIAM, 2nd edition, 2001.[8] J. Carr.
Applications of centre manifold theory , volume 35 of
AppliedMath. Sci.
Springer–Verlag, 1981.[9] Shiyi Chen and G. D. Doolen. Lattice Boltzmann methodfor fluid flows.
Annu. Rev. Fluid Mechanics , 30:329–364, 1998.doi:10.1146/annurev.fluid.30.1.329.4010] J. Cisternas, C. W. Gear, S. Levin, and I. G. Kevrekidis.Equation-free modeling of evolving diseases: Coarse-grainedcomputations with individual-based models. Technical report,[ http://arXiv.org/abs/nlin.AO/0310011 ], 2003.[11] M. C. Cross and P. C. Hohenberg. Pattern formation out-side of equilibrium.
Rev. Mod. Phys. , 65(3):851–1112, 1993.doi:10.1103/RevModPhys.65.851.[12] J. Dolbow, M. A. Khaleel, and J. Mitchell. Multiscale mathematicsinitiative: A roadmap. Report from the 3rd DoE workshop on mul-tiscale mathematics. Technical report, Department of Energy, USA, , December 2004.[13] Weinan E, Bjorn Engquist, Xiantao Li, Weiqing Ren, and Eric Vanden-Eijnden. The heterogeneous multiscale method: A review. Technicalreport, ,2004.[14] C. Foias, M. S. Jolly, I. G. Kevrekidis, G. R. Sell, and E. S. Titi. Onthe computation of inertial manifolds.
Phys. Lett. A , 131:433–436, 1988.doi:10.1016/0375-9601(88)90295-2.[15] C. Foias and E. S. Titi. Determining nodes, finite differenceschemes and inertial manifolds.
Nonlinearity , 4:135–153, 1991.doi:10.1088/0951-7715/4/1/009.[16] Martin J. Gander and Andrew M. Stuart. Space-time continuous anal-ysis of waveform relaxation for the heat equation.
SIAM Journal onScientific Computing , 19(6):2014–2031, 1998.[17] C. W. Gear, Ju Li, and I. G. Kevrekidis. The gap-toothmethod in particle simulations.
Phys. Lett. A , 316:190–195, 2003.doi:10.1016/j.physleta.2003.07.004.[18] Johannes Giannoulis, Michael Herrmann, and Alexander Mielke. La-grangian and hamiltonian two-scale reduction. Technical report, http://arxiv.org/abs/0802.2820v1 , 2008.4119] M. S. Jolly, I. G. Kevrekidis, and E. S. Titi. Approximate inertial man-ifolds for the Kuramoto–Sivashinsky equation: analysis and computa-tions.
Physica D , 44:38–60, 1990. doi:10.1016/0167-2789(90)90046-R.[20] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Run-borg, and K. Theodoropoulos. Equation-free, coarse-grained multiscalecomputation: enabling microscopic simulators to perform system leveltasks.
Comm. Math. Sciences , 1:715–762, 2003.[21] Alexander Kurganov and Eitan Tadmor. New high-resolutioncentral schemes for nonlinear conservation laws and convection-diffusion equations.
J. Computational Physics , 160:241–282, 2000.doi:10.1006/jcph.2000.6459.[22] Y. A. Kuznetsov.
Elements of applied bifurcation theory , volume 112 of
Applied Mathematical Sciences . Springer–Verlag, 1995.[23] J. Li, C. W. Gear P. G. Kevrekidis, and I. G. Kevrekidis.Deciding the nature of the coarse equation through mi-croscopic simulation: an augmented Lagrangian approach.
SIAM Multiscale Modeling and Simulation , 1:391–407, 2003. http://epubs.siam.org/sam-bin/dbq/article/41916 .[24] T. MacKenzie and A. J. Roberts. Accurately model the Kuramoto–Sivashinsky dynamics with holistic discretisation.
SIAM J. Ap-plied Dynamical Systems , 5(3):365–402, 2006. doi:10.1137/050627733 http://epubs.siam.org/SIADS/volume-05/art_62773.html .[25] M. Marion and R. Temam. Nonlinear Galerkin meth-ods.
SIAM J. Numer. Anal. , 26(5):1139–1157, 1989. http://locus.siam.org/SINUM/volume-26/art_0726063.html .[26] G. N. Mercer and A. J. Roberts. A centre manifold de-scription of contaminant dispersion in channels with varyingflow properties.
SIAM J. Appl. Math. , 50:1547–1565, 1990. http://link.aip.org/link/?SMM/50/1547/1 .[27] Mauro Mobilia, Ivan T. Georgiev, and Uwe C. Tauber.Spatial stochastic predator-prey models. Technical report, http://arxiv.org/abs/q-bio.PE/0609039 , 2006.4228] G. A. Pavliotis and A. M. Stuart.
An introduction to multiscalemethods . University of Warwick and Imperial College London, 2006. .[29] A. J. Roberts. The utility of an invariant manifold description of theevolution of a dynamical system.
SIAM J. Math. Anal. , 20:1447–1458,1989. http://locus.siam.org/SIMA/volume-20/art_0520094.html .[30] A. J. Roberts. Holistic discretisation ensures fidelity to Burg-ers’ equation.
Applied Numerical Modelling , 37:371–396, 2001.doi:10.1016/S0168-9274(00)00053-2.[31] A. J. Roberts. Holistically discretise the Swift-Hohenberg equa-tion on a scale larger than its spatial pattern. Technical report, http://arXiv.org/abs/math.NA/0110153 , 2001.[32] A. J. Roberts. Simple and fast multigrid solution of Poisson’s equationusing diagonally oriented grids.
ANZIAM J. , 43(E):E1–E36, July 2001. http://anziamj.austms.org.au/V43/E025 .[33] A. J. Roberts. A holistic finite difference approach models linear dy-namics consistently.
Mathematics of Computation , 72:247–262, 2002. .[34] A. J. Roberts. Computer algebra derives discretisationsvia self-adjoint multiscale modelling. Technical report, http://eprints.usq.edu.au/4275/ , 2008.[35] A. J. Roberts. Model dynamics on a multigrid across multiple lengthand time scales. Technical report, http://arxiv.org/abs/0802.1098 ,2008.[36] A. J. Roberts. Resolve subgrid microscale interactions to discretisestochastic partial differential equations. preprint, January 2008.[37] A. J. Roberts and I. G. Kevrekidis. Higher order accuracy inthe gap-tooth scheme for large-scale dynamics using microscopicsimulators. In Rob May and A. J. Roberts, editors,
Proc. of12th Computational Techniques and Applications Conference CTAC-2004 , volume 46 of
ANZIAM J. , pages C637–C657, July 2005. http://anziamj.austms.org.au/V46/CTAC2004/Robe [July 20, 2005].4338] A. J. Roberts and I. G. Kevrekidis. General tooth boundary conditionsfor equation free modelling.
SIAM J. Scientific Computing , 29(4):1495–1510, 2007. doi:10.1137/060654554.[39] J. C. Robinson. The asymptotic completeness of in-ertial manifolds.
Nonlinearity , 9:1325–1340, 1996. .[40] G. Samaey, I. G. Kevrekidis, and D. Roose. The gap-tooth scheme forhomogenization problems.
SIAM Multiscale Modeling and Simulation ,4:278–306, 2005. doi:10.1137/030602046.[41] J. Sijbrand. Properties of centre manifolds.
Trans. Amer. Math. Soc. ,289:431–469, 1985.[42] Matthew J. Simpson, Alistair Merrifield, Kerry A. Landman, andBarry D. Hughes. Simulating invasion with cellular automata: Con-necting cell-scale and population-scale properties.
Physical Review E ,76:021918, 2007. doi:10.1103/PhysRevE.76.021918.[43] R. Temam. Inertial manifolds.
Mathematical Intelligencer , 12:68–74,1990.[44] James L. Thomas, Boris Diskin, and Achi Brandt. Textbook multigridefficiency for fluid simulations.
Annu. Rev. Fluid Mechanics , 35:317–340,2003. doi:10.1146/annurev.fluid.35.101101.161209.[45] T. P. Witelski and A. J. Bernoff. Self-similarasymptotics for linear and nonlinear diffusion equa-tions.
Studies in Applied Maths. , 100:153–193, 1998. .[46] D. J. Wollkind, V. S. Manoranjan, and L. Zhang. Weakly nonlinearanalyses of prototype reaction-diffusion model equations.