Learning high-order spatial discretisations of PDEs with symmetry-preserving iterative algorithms

J. E. Bunder and A. J. Roberts
School of Mathematical Sciences,
University of Adelaide, South Australia \urlmailto:judith.bunder@adelaide.edu.au, orcid:0000-0001-5355-2288\urlmailto:profajroberts@protonmail.com, orcid:0000-0001-8930-1552

Abstract

Common techniques for the spatial discretisation of pdes on a macroscale grid include finite difference, finite elements and finite volume methods. Such methods typically impose assumed microscale structures on the subgrid fields, so without further tailored analysis are not suitable for systems with subgrid-scale heterogeneity or nonlinearities. We provide a new algebraic route to systematically approximate, in principle exactly, the macroscale closure of the spatially-discrete dynamics of a general class of heterogeneous non-autonomous reaction-advection-diffusion pdes. This holistic discretisation approach, developed through rigorous theory and verified with computer algebra, systematically constructs discrete macroscale models through physics informed by the pde out-of-equilibrium dynamics, thus relaxing many assumptions regarding the subgrid structure. The construction is analogous to recent gray-box machine learning techniques in that predictions are directed by iterative layers (as in neural networks), but informed by the subgrid physics (or ‘data’) as expressed in the pdes. A major development of the holistic methodology, presented herein, is novel inter-element coupling between subgrid fields which preserve self-adjointness of the pde after macroscale discretisation, thereby maintaining the spectral structure of the original system. This holistic methodology also encompasses homogenisation of microscale heterogeneous systems, as shown here with the canonical examples of heterogeneous 1D waves and diffusion.

1 Introduction

The observable macroscale dynamics of a multiscale system are typically driven by coherent structures governed by microscopic agents, such as electrons, molecules, or individuals in a population. While such systems are often well-understood at the microscale, large-scale simulations of the full microscale model are usually infeasible, even on the largest high-performance computers, and in particular for nonlinear system with fluctuations that persist and interact across multiple scales [Nasa2018]. Multiscale modelling aims to overcome this issue by efficiently approximating over the multiple scales to provide an accurate numerical model of the coherent dynamics of the system at the macroscale of interest. Numerous computational schemes have been developed for a wide variety of multiscale systems, from composite materials [Dutra2020] to plasmas [Shohet2020]. However, many of these computational schemes are developed for specific systems and are not easily adaptable to other cases. Furthermore, supporting theory is often only rigorously applicable in the limiting case of infinite scale separation (typically characterised as $\epsilon\rightarrow 0$ ). Here we further develop the holistic discretisation methodology for multiscale nonlinear systems—a general purpose approach—supported by analytical theory and straightforwardly adaptable to a wide range of spatio-temporal systems at finite scale separation.

To determine a macroscale description of a multiscale system with microscale heterogeneity we seek an homogenisation, or ‘average’, over the microscale which accounts for the heterogeneity without retaining unnecessary fine-scale details [Geers2107]. For example, for composite materials with periodic microstructures, which are common in nature and are increasingly synthesised for novel industrial applications [[, e.g.,]]NematNasser2011,Bargmann2018, asymptotic homogenisation constructs a power series in scale separation $\epsilon$ with coefficients dependent on a periodic ‘representative volume element’ (e.g., a unit cell) which describes the microscale heterogeneity [Dutra2020, Feppon2020]. However, homogenisation often requires substantial user input in the form of an algebraic analysis of the given microscale system prior to numerical implementation, and theoretical support is often only guaranteed in unphysical circumstances; for example, typically one needs to explicitly identify ‘fast’ and ‘slow’ variables, and also assume an unphysical infinite scale separation between the two [[, e.g.,]]Engquist08. In contrast, the holistic homogenisation methodology, detailed in section 2.1, requires little or no such user input and generally permits an immediate application of the numerical scheme. It is thus more adaptable and accessible than many other common homogenisation methods. Theoretical support for this methodology is provided by centre manifold theory, as discussed in section 2.2, and by the consistency analysis of section 3, and is broadly applicable to many physical scenarios.

The holistic homogenisation methodology was originally developed by [Roberts98a, e.g., ], and has been applied to several different systems [[, e.g.,]]Roberts00a, MacKenzie09b, Roberts2011a, Jarrad2016a. The fundamental idea is to partition the domain of the original system into $N$ disjoint elements, with these $N$ elements defining a macroscale grid, and then to discover appropriate microscale fields that are slaved to a useful set of macroscale variables. In this article, the important new development is the preservation of self-adjointness by appropriate design of the coupling between the $N$ elements, as defined in section 2.1, which ensures that fundamental dynamical features of the multiscale system (specified by the eigenstructure) are maintained—especially for homogenisation as proved in section 2.4. To illustrate the adaptability of this holistic methodology, section 4 considers a pde with periodic heterogeneous diffusion. The problem is firstly embedded in a family of general heterogeneous diffusion pdes, which rephrases the heterogeneity as solely dependent on a ‘phase’ coordinate, where the periodic phase domain is reminiscent of the representative volume element commonly utilised in the asymptotic homogenisation of periodic microstructures. This embedding empowers a straightforward application of the numerical holistic methodology.

In this article, to establish the new homogenisation methodology we mostly consider the spatio-temporal dynamics of a field $u(t,x)$ satisfying non-autonomous reaction-advection-diffusion pdes in the general form

u_{t}=-f(x,u,u_{x})_{x}+\alpha g(t,x,u,u_{x})\,,

(1)

where subscripts $x$ and $t$ denote spatial and temporal derivatives, respectively, and functions $f$ and $g$ are characterised by section 1. The homogenisation derives accurate macroscale spatial discretisations of these nonlinear pdes on a specified macroscale grid. The homogenisation accounts for subgrid (i.e., microscale) dynamics and symmetries, and hence automatically leads to stable discretisations, assuming stability of the original system. We also establish that the same approach generates accurate and stable spatial discretisations of 1D first-order wave pdes (section 3.1) and 1D second-order nonlinear wave pdes obtained by replacing $u_{t}$ in eq. 1 by $u_{tt}$ (section 2.3).

With some adaptations, this methodology is also suitable to multi-dimensional space pdes, $u_{t}=-\text{\boldmath$\nabla$}\cdot{\text{\boldmath$f$}}({\text{\boldmath$x$}},u,\text{\boldmath$\nabla$}u)+\alpha g(t,{\text{\boldmath$x$}},u,\text{\boldmath$\nabla$}u)$ , and is related to the multidimensional homogenisation (albeit non-self-adjoint) of [Roberts2011a], but we leave multi-dimensional space for future research and here we focus the case of 1D space.

Assumption 1.

Firstly, functions $f$ and $g$ on the right-hand side of pde eq. 1 are to be smooth functions of their arguments, and the function $g$ varies relatively slowly in time $t$ . Secondly, in the cases of second-order pdes, the function $f$ , called the flux, is strictly monotonic decreasing in $u_{x}$ : that is, for some positive $\nu$ , $\mathchoice{\frac{\partial f}{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}\leq-\nu<0$ with $f(x,u,0)=0$ . Lastly, for strict support by existing theory, we assume $f$ and $g$ are such that the pde eq. 1, with ‘edge conditions’ eq. 6, satisfies the requirements of the invariant manifold theory by [Hochs2019].

We define smooth (section 1) to mean either infinitely differentiable, $C^{\infty}$ , or differentiable to an order $p$ sufficient for the purposes at hand, $C^{p}$ . This smoothness requirement still permits microscale heterogeneity where the functions $f$ and $g$ vary significantly on a spatial scale of the order of a microscale length $d$ , such as for heterogeneous diffusion (section 4).

We denote the 1D spatial domain of the pde eq. 1 by $\mathbb{X}$ , which in examples is often over the interval $(0,\mathfrak{L})\subset\mathbb{R}$ , and consider spatial dependence of field $u$ in the Hilbert space $\mathbb{H}$ of square integrable, twice differentiable, functions on $\mathbb{X}$ , and often in the subspace $\mathbb{L}$ of all $\mathfrak{L}$ -periodic functions in $\mathbb{H}$ (i.e., we often invoke periodic boundary conditions). We partition the domain $\mathbb{X}$ into $N$ disjoint open interval elements $\mathbb{X}_{j}$ for $j\in\mathbb{J}:=\{1,\ldots,N\}$ and define $\tilde{}\mathbb{X}:=\cup_{j\in\mathbb{J}}\mathbb{X}_{j}$ . Let $u_{j}(t,x)$ denote subgrid fields which satisfy the pde eq. 1 in $\mathbb{X}_{j}$ : that is, $u_{jt}=-f(x,u_{j},u_{jx})_{x}+\alpha g(t,x,u_{j},u_{jx})$ , for $x\in\mathbb{X}_{j}$ . Define macroscale, coarse, ‘grid’ values

U_{j}(t):=\mathcal{U}\big{[}u_{j}(t,x)\big{]}\quad\text{for every }j\in\mathbb{J},

(2)

for some chosen projection $\mathcal{U}:C(\mathbb{X}_{j})\to\mathbb{R}$ , a projection that is some measure of the amplitude of the field $u_{j}$ in the $j$ th element.¹¹1The precise form of $\mathcal{U}$ is largely immaterial because the evolution of states in the slow subspace is independent of how we choose to parametrise the subspace [Roberts2014a, Lemma 5.1]. Hence the macroscale projection $\mathcal{U}$ may be subjectively chosen according to any reasonable measure of the subgrid field $u_{j}$ within each element $\mathbb{X}_{j}$ . Consequently, a spatial discretisation of the pde eq. 1 is a closed set of odes for the vector of macroscale grid variables ${\text{\boldmath$U$}}:=(U_{1},\ldots,U_{N})$ in the form

\mathchoice{\frac{d{\text{\boldmath$U$}}}{dt}}{{d{\text{\boldmath$U$}}}/{dt}}{{d{\text{\boldmath$U$}}}/{dt}}{{d{\text{\boldmath$U$}}}/{dt}}={\text{\boldmath$G$}}(t,{\text{\boldmath$U$}}).

(3)

Four examples discussed herein are the ode systems eqs. 5, 9c, 21b and 36b modelling, respectively, nonlinear wave modulation, spatial diffusion, progressive waves, and homogenised heterogeneous diffusion. section 2.2 uses centre manifold theory (cmt) [[, e.g.,]]Carr81, Haragus2011, Hochs2019 to do three things in support of such a macroscale evolution eq. 3. Firstly, we develop further an approach to proving that in principle there exists an exact macroscale closure eq. 3 to the dynamics of the pde eq. 1. Secondly, cmt establishes that such a closure is emergent from general initial conditions. Thirdly, cmt justifies constructing new systematic approximations to the in-principle closure by learning subgrid field corrections from the pde eq. 1.

Traditional spatial discretisations of pdes eq. 1, whether finite difference, finite element, or finite volume, impose a set of assumed subgrid fields within each element and then derive approximate rules for the macroscale evolution of the parameters of the imposed fit [[, e.g.,]]Efendiev2004, Geers2010. The accuracy of this modelling by discretisation, particularly for nonlinear systems, is strongly dependent on the class of imposed subgrid microstructure and its ability to capture the microscale dynamics which emerge and persist across the macroscale [[, e.g.,]]Aarnes2007, Zhang2020. Our dynamical systems (holistic) approach is to learn the algebraic subgrid fields from the pde eq. 1. Using more general subgrid structures is cognate to the aims of the so-called Generalised/Extended Finite Element Methods [[, e.g.,]]Turner2011b, but crucially we let the pde determine the subgrid structures from the powerful framework of centre manifolds in dynamical systems. The so-called stabilized scheme [[, e.g.,]]Hughes95 appears analogous to the first step of our construction in section 2.2. We learn improved subgrid fields via computer algebra using the governing pde eq. 1, with second and further iterations providing higher-order terms. The resultant holistic discretisation is both accurate and flexible—it naturally captures fine-scale structures which persist across multiple scales, and it does so to a chosen order of accuracy that is routinely achievable via iteration.

A previous alternative holistic approach learns subgrid fields by systematically refining a piecewise constant initial approximation [[, e.g.,]]Roberts98a, Roberts00a, Roberts2011a—an approach that adapts to the multi-scale gap-tooth scheme [[, e.g.,]]Roberts06d, Kevrekidis09a. Our new approach extends and complements previous schemes by systematically preserving self-adjoint symmetries in the learnt macroscale closure.

Example 1 (nls).

As an introductory nontrivial example, consider the nonlinear 1D Schrödinger equation (nls) governing a complex valued field $u(t,x)$ :

\operatorname{i}\partial_{t}u=-\tfrac{1}{2}\partial_{x}^{2}u+\alpha|u|^{2}u\,,

(4)

with periodic boundary conditions on the spatial domain $[-\pi,\pi]$ . The real and imaginary parts of the field $u$ typically evolve to be out-of-phase and each drives persistent oscillations in the other. Common applications include the modulation of light propagating in nonlinear optical fibres and planar waveguides [Kibler2012, e.g.]. The nonlinear wave nature of this example tests our discretisation scheme.

We discretise space into $N$ elements $\mathbb{X}_{j}$ of uniform width $H$ centred about grid points $X_{j}$ . Then, defining macroscale complex variables $U_{j}(t):=u_{j}(t,X_{j})$ (i.e., here the projection $\mathcal{U}$ of eq. 2 is the subgrid field $u_{j}$ at the midpoint $X_{j}$ of the $j$ th element), our holistic homogenisation scheme iteratively derives a macroscale discretisation of eq. 4. The details of the algebraic machinations of the iteration are not important here—computer algebra like that of appendices A, B and C handles the details. The result here is that the macroscale discretisation is, with asterisk superscript denoting complex conjugation,

$\displaystyle\operatorname{i}\partial_{t}U_{j}$	$\displaystyle=\alpha\|U_{j}\|^{2}U_{j}-\gamma^{2}\tfrac{1}{8H^{2}}(U_{j+2}+U_{j-2}-2U_{j})$
	$\displaystyle\quad{}+\alpha\gamma^{2}\tfrac{1}{96}\big{[}2\|U_{j}\|^{2}(U_{j+2}+U_{j-2}-4U_{j})$
	$\displaystyle\qquad{}-4\|U_{j+2}\|^{2}(U_{j+2}-2U_{j})-4\|U_{j-2}\|^{2}(U_{j-2}-2U_{j})$
	$\displaystyle\qquad{}-4U_{j}^{}(U_{j+1}^{2}+U_{j-1}^{2}-U_{j+1}U_{j-1})-4U_{j}(U_{j+1}^{}U_{j-1}+U_{j-1}^{*}U_{j+1})$
	$\displaystyle\qquad{}+U_{j+2}^{}(U_{j}^{2}-2U_{j+1}^{2})+U_{j-2}^{}(U_{j}^{2}-2U_{j-1}^{2})\big{]}+\mathcal{O}\mathchoice{\big{(}\alpha^{2},\gamma^{3}\big{)}}{\big{(}\alpha^{2},\gamma^{3}\big{)}}{(\alpha^{2},\gamma^{3})}{(\alpha^{2},\gamma^{3})},$	(5)

where $\gamma$ is a coupling order parameter labelling each inter-element interaction (section 2.1). Set $\gamma=1$ to obtain the full inter-element coupling needed to model the pde, and then the first line of the macroscale evolution eq. 5 is a standard discretisation of the nls eq. 4 with the spatial derivatives over a step of $2H$ . The subsequent terms in eq. 5 are the leading terms due to subgrid physical interactions, terms learnt by the algebraic methodology developed in this article.

Refer to caption — Figure 1: Simulations of the nls pde compare predictions of field magnitudes: (left) $|u(t,x)|$ computed using a fine-scale discretisation of eq. 4; and (right) $|U(t,x)|$ computed using the coarse scale holistic discretisation eq. 5. The surfaces are coloured according to the arguments of the complex fields, as indicated by the side-bar.

For an example simulation we choose the initial condition $u(0,x)=1-\operatorname{i}\operatorname{sech}x$ which produces Kuznetsov–Ma breather solutions of nonlinear localised oscillating peaks on a non-zero background [Zhao2017]. Breathers and solitons arise in nonlinear pde when there is a balance between the dispersion and the nonlinearity. fig. 1 illustrates the success of the holistic homogenisation by comparing a microscale simulation of the nls eq. 4 (left plot) with the evolution of the macroscale discretisation eq. 5 truncated to errors $\mathcal{O}\mathchoice{\big{(}\alpha^{2},\gamma^{4}\big{)}}{\big{(}\alpha^{2},\gamma^{4}\big{)}}{(\alpha^{2},\gamma^{4})}{(\alpha^{2},\gamma^{4})}$ (right). The microscale simulation of the nls employs a spatial grid of $3300$ points with spacing $h=0.0019$ , whereas the holistic scheme discretises the domain with $N=151$ elements of width $H=0.042$ , which is over twenty times larger and correspondingly less stiff, and so much quicker computationally. At all times we observe fine-scale simulation errors (via Matlab’s ode solver ode113) which do not appear in the simulation with the holistic discretisation eq. 5 (via a fourth order Runge–Kutta).²²2Matlab’s ode113 was the only Matlab ode solver which could simulate the fine-scale model (4) and provide reasonable accuracy (time taken: 8 564 s). In contrast, Matlab ode solvers had no trouble accurately simulating the holistic discretisation (5) (e.g., ode113 time was 246 s), but fourth order Runge–Kutta was substantially faster (time 102 s). The fourth order Runge–Kutta could not accurately simulate the fine-scale model.

Remark 1.

Inter-element coupling is crucial. In order to construct an accurate macroscale discretisation the nature of the coupling between adjacent elements $\mathbb{X}_{j}$ is crucial. For example, in the context of material homogenisation, [Abdulle2020b] in their abstract comment that

A naive treatment of these artificial boundary conditions leads to a first order error … . This error dominates all other errors originating from the discretization of the macro and the micro problems, and its reduction is a main issue in today’s engineering multiscale computations.

The “naive treatment” in this comment corresponds to a low-order approximation in the element coupling order parameter $\gamma$ , but with our iterative algorithm higher order approximations are straightforward. The problematic “boundary conditions” in the comment we term the edge conditions (or coupling conditions)—it is these edge conditions that define the coupling of elements $\mathbb{X}_{j}$ to form the spatially complete system. section 2.1 develops inter-element edge conditions that in a variety of scenarios (e.g., nls of example 1) provably preserve self-adjointness (section 2.4), and have provable accuracy (sections 3.1 and 3.3). Thus this article develops a valuable approach to addressing a key “issue in today’s engineering multiscale computations”.

Remark 2.

Initial conditions. For any macroscale discrete closure such as eq. 3, initial conditions are nontrivial. The paradox is that, despite the definition eq. 2 that the macroscale grid variable $U_{j}(t):=\mathcal{U}[u_{j}(t,x)]$ , in order to obtain accurate forecasts over long times, generally the initial grid value $U_{j}(0)\neq\mathcal{U}[u_{j}(0,x)]$ [Roberts89b, Roberts01a]. In physics, [vanKampen85] termed this phenomena the “initial slip”. There is, in effect, a ‘boundary layer’ in time that accounts for non-trivial rapid transients. Using the geometry of invariant manifolds, [Roberts89b] developed an efficient general method to algebraically learn the correct initial $U_{j}(0)$ for accurate long-time forecasts by reduced dimensional models such as a discretisation of a pde. The method has been applied in various scenarios [[, e.g.,]]Roberts2014a. Further research is needed to form initial conditions for the discretisation framework established herein.

Remark 3.

Analogy with machine learning. The iterative algorithm used to construct the holistic homogenisation of example 1, and the other examples herein, was originally developed by [Roberts96a]. At each iteration the algorithm evaluates residuals of the nonlinear pde within the elements, and improves the resolution of the subgrid fields with a linear correction based upon this subgrid information. We draw an analogy between this iteration and a machine learning algorithm where an AI learns the generic form of the macroscale evolution from many thousands of simulations [[, e.g.,]]GonzalezGarcia98, Frank2020, Linot2020, but here the AI is a ‘smart’ grey-box which is directed by additional algebraic knowledge/data. For example, [BarSinai2018] imposed a conservative form on their learning of a coarse-scale closure to Burgers’ pde, and [Mercer94a] constructed stable schemes to advection-diffusion in pipes that matched the long term evolution to high order. In constructing the holistic discretisation, each iteration is analogous to one layer in a deep neural network: evaluating residuals is analogous to a nonlinear neurone function; and the linear corrections are analogous to using weighted linear combinations of outputs (i.e., activation functions) of one layer as the inputs of the next layer. The analogy requires that, at each execution, the smart neural network is directed by the holistic algorithm to specifically craft the neurones and linear weights to the problem at hand. Being algebraic, this data which directs the algorithm encompasses all points in the state space’s domain, not just sample data as in machine learning. Consequently, we contend that mathematicians have for many decades been doing smart analogues of machine learning. Such algebraic learning empowers the physical interpretation, validation, and verification required by modern science [[, e.g.]]Brenner2021.

2 Self-adjoint preserving coupling

Let the 1D spatial discretisation of the general pde (1) be parametrised by ${\text{\boldmath$U$}}=(U_{1},\ldots,U_{N})$ as defined by some projection (2) and satisfy an evolution ode of the form eq. 3. The discretisation is derived from the subgrid fields $u_{j}(t,x)$ of pde (1) over the $N$ disjoint elements $\mathbb{X}_{j}$ . To construct the spatial discretisation and ode eq. 3 we must specify $u_{j}$ at the two edges of element $\mathbb{X}_{j}$ . Specifying the fields $u_{j}$ at the edges of element $\mathbb{X}_{j}$ defines an inter-element coupling. Such coupling conveys information across space and hence engenders the macroscale dynamics. While there are many possible coupling conditions, here we establish new coupling conditions which both preserve self-adjoint symmetry in the original microscale pde, and also ensure the macroscale dynamics are correctly consistent to an arbitrarily high order of accuracy. As these new coupling conditions constrain fields on the edges of the elements $\mathbb{X}_{j}$ , we refer to them as edge conditions.

We identify the right and left edge-points of the elements by writing the $j$ th element as the open interval $\mathbb{X}_{j}=(L_{j},R_{j})$ . Because the elements abut, the right-edge of the $j$ th element is the left-edge of the $(j+1)$ th element, $R_{j}=L_{j+1}$ , and thus, to recover the original pde (1) over the entire domain, we aim to set $u_{j}(t,R_{j})=u_{j+1}(t,L_{j+1})$ .³³3Strictly, such edge values are the limit as $x$ approaches the edge of the (open) elements. However, to understand and theoretically support the discretisation on the elements, we parameterise a range of inter-element coupling to best control the information flow between elements. We introduce an real-valued inter-element coupling parameter $\gamma$ : when $\gamma=0$ the elements are uncoupled and form a base for some theory; when $\gamma=1$ the elements are fully coupled to recover the original pde problem eq. 1 for field $u(t,x)$ . A solution field $u(t,x)$ then arises as the collection of subgrid fields $u_{j}(t,x)$ . Further, we define the order parameter $\gamma$ so as to ‘label’ each inter-element communication so that a term in $\gamma^{p}$ expresses a composition of interactions between each element and $p$ pairs of its nearest neighbour elements. Consequently, truncating asymptotic expansions to ‘errors’ $o\big{(}\gamma^{p}\big{)}$ naturally and automatically creates local discretisations of stencil width $(2p+1)$ in space.

2.1 Couple the elements

For conciseness, denote the flux in the $j$ th element as $f_{j}(t,x):=f(x,u_{j},u_{jx})$ . Further, use superscripts $R,L$ to denote evaluation at the corresponding right/left edges $x=R_{j},L_{j}$ : for example, $u_{j}^{R}(t):=u_{j}(t,R_{j})$ and $f_{j}^{L}(t):=f_{j}(t,L_{j})$ . In addition to the inter-element coupling parameter $\gamma$ , we introduce a second real-valued parameter $\theta$ that flexibly ‘tunes’ the details of inter-element communication; typically $|\theta|\leq 1$ . For example, for advective processes, $\theta$ controls the upwind character of the discrete model. At every time and every element $j$ , we couple elements with the two edge conditions


$\displaystyle(1-\tfrac{1}{2}\gamma)(u_{j}^{R}-u_{j}^{L})$	$\displaystyle=\tfrac{1}{2}\gamma(u_{j+1}^{L}-u_{j-1}^{R})+\tfrac{1}{2}\gamma\theta(u_{j}^{L}+u_{j}^{R})-\tfrac{1}{2}\gamma\theta(u_{j+1}^{L}+u_{j-1}^{R})\,,$	(6a)
$\displaystyle(1-\tfrac{1}{2}\gamma)(f_{j}^{R}-f_{j}^{L})$	$\displaystyle=\tfrac{1}{2}\gamma(f_{j+1}^{L}-f_{j-1}^{R})-\tfrac{1}{2}\gamma\theta(f_{j}^{L}+f_{j}^{R})+\tfrac{1}{2}\gamma\theta(f_{j+1}^{L}+f_{j-1}^{R}).$	(6b)

When $u(t,x)$ is to be $\mathfrak{L}$ -periodic in space we require the periodicity $u^{e}_{j+N}=u^{e}_{j}$ and $f^{e}_{j+N}=f^{e}_{j}$ for $e\in\{L,R\}$ . section 2.4 proves that these inter-element edge conditions preserve self-adjointness in the pde eq. 1.

•

When parameter $\gamma=0$ , the edge conditions eq. 6 reduce to

$\displaystyle u_{j}^{R}-u_{j}^{L}=0\,,\quad f_{j}^{R}-f_{j}^{L}=0\,.$ (7)

That is, each element is uncoupled from the others when $\gamma=0$ . Consequently we find that $\gamma=0$ forms a useful base for theory. Denoting the element lengths by $H_{j}:=R_{j}-L_{j}$ , the two conditions (7) specify that the subgrid field $u_{j}$ is $H_{j}$ -periodic. For the forthcoming example 2 on the diffusion pde, $u_{t}=u_{xx}$ , such uncoupled elements evolve in time so that $u_{j}\to{}$ constant on the cross-element diffusion time scale of $1/H_{j}^{2}$ . Centre manifold theory then supports the existence and emergence of an exact closure to the spatial discretisation, namely (3), in some finite range of $\gamma$ about $\gamma=0$ (section 2.2).

•

When parameter $\gamma=1$ the edge conditions eq. 6 reduce to

	$\displaystyle\tfrac{1}{2}(1-\theta)(u_{j}^{R}-u_{j+1}^{L})=\tfrac{1}{2}(1+\theta)(u_{j}^{L}-u_{j-1}^{R}),$
	$\displaystyle\tfrac{1}{2}(1+\theta)(f_{j}^{R}-f_{j+1}^{L})=\tfrac{1}{2}(1-\theta)(f_{j}^{L}-f_{j-1}^{R}).$

In matrix notation the $\gamma=1$ edge condition on field $u$ is equivalent to $V\tilde{}{\text{\boldmath$u$}}=\text{\boldmath$0$}$ for vector $\tilde{}{\text{\boldmath$u$}}:=(u_{1}^{R}-u_{2}^{L},\ldots,u_{N}^{R}-u_{1}^{L})$ , and for $N\times N$ circulant matrix

V:=\begin{bmatrix}\tfrac{1}{2}(1-\theta)&0&\cdots&0&\tfrac{1}{2}(1+\theta)\\ \tfrac{1}{2}(1+\theta)&\ddots&&\ddots&0\\ 0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\ 0&\cdots&0&\tfrac{1}{2}(1+\theta)&\tfrac{1}{2}(1-\theta)\end{bmatrix}.

The eigenvalues of this circulant matrix are $\lambda_{k}=\tfrac{1}{2}(1-\theta)+\tfrac{1}{2}(1+\theta)e^{-2\pi\operatorname{i}k(N-1)/N}$ for $k=0,\ldots,N-1$ [Gray2006, Sec. 3.1, Eq. (3.7)], and thus only in the case $k=N/2$ and $\theta=0$ does $V$ have a zero eigenvalue (is singular). Therefore, the only solution of $V\tilde{}{\text{\boldmath$u$}}=\text{\boldmath$0$}$ is $\tilde{}{\text{\boldmath$u$}}=\text{\boldmath$0$}$ (except for the case $\theta=0$ and $N$ even, when there are also nontrivial zigzag solutions in the nullspace). A similar argument holds for the flux condition which is equivalent to $V^{T}\tilde{}{\text{\boldmath$f$}}=\text{\boldmath$0$}$ for vector $\tilde{}{\text{\boldmath$f$}}:=(f_{1}^{L}-f_{N}^{R},f_{2}^{L}-f_{1}^{R},\ldots,f_{N}^{L}-f_{N-1}^{R})$ ; that is, the only solution is $\tilde{}{\text{\boldmath$f$}}=\text{\boldmath$0$}$ (except the case $\theta=0$ and $N$ is even). Therefore, for $\gamma=1$ , edge conditions eq. 6 are equivalent to

u_{j}^{R}-u_{j+1}^{L}=0\,,\quad f_{j}^{R}-f_{j+1}^{L}=0

(8)

(except perhaps for the isolated case $\theta=0$ and $N$ even). Thus $\gamma=1$ restores full inter-element coupling by requiring continuity of the field $u$ and the flux $f$ between elements.

Remark 4.

For the special case of $\theta=0$ and even $N$ , both the continuous edge conditions (8) and a discontinuous zigzag mode, $u_{j}^{R}-u_{j+1}^{L}\propto(-1)^{j}$ and $f_{j}^{R}-f_{j+1}^{L}\propto(-1)^{j}$ , satisfy conditions (6) at $\gamma=1$ . Such a macroscale zigzag mode typically occurs with periodic boundary conditions and an even number of grid points. The zigzag mode is physically deficient because it interpolates from the macro-grid as $\cos[\pi(x-X_{j})/H]+a\sin[\pi(x-X_{j})/H]$ for any arbitrary coefficient $a$ . All other modes have a unique interpolation in resolved wavenumbers $|k|\leq\pi/H$ . Because of this physical deficiency, when necessary we restrict attention to the co-dimension two space without the zigzag modes. Without the zigzag modes, the edge conditions eq. 6 are equivalent to (8).

For a physically reasonable description of the macroscale dynamics, we need to define some measure of the overall size of the fields $u_{j}$ in each element $\mathbb{X}_{j}$ , $j=1,\ldots,N$ , in order to give physical meaning to the macroscale parameters ${\text{\boldmath$U$}}=(U_{1},\ldots,U_{N})$ , that is, we need to define the projection $\mathcal{U}$ invoked by eq. 2. Among many, two possibilities are the mid-element value $\mathcal{U}[u_{j}]:=u_{j}(t,(L_{j}+R_{j})/2)$ (used in example 1), or the element average $\mathcal{U}[u_{j}]:=\frac{1}{H_{j}}\int_{L_{j}}^{R_{j}}u_{j}(t,x)\,dx$ . Any reasonable choice suffices [Roberts2014a, §5.3.3]. The next example implements the element average.

Example 2 (high-order consistency in discretising diffusion).

Consider the simplest case of pde eq. 1, the case of homogeneous diffusion $u_{t}=u_{xx}$ , which arises when $g=0$ and the flux $f=-u_{x}$ . appendix A describes computer algebra code that discretises this diffusion pde via elements, all of size $H$ , with inter-element coupling controlled by edge conditions eq. 6. The computer algebra analyses the pde to learn the emergent slow manifold model as a power series in the coupling parameter $\gamma$ (remark 3). Because of the simplicity of the homogeneous diffusion operator, the algebraically learnt slow manifold is here local polynomials. In terms of subgrid space variable $\xi:=(x-L_{j})/H\in[0,1]$ , tuning parameter $\theta$ , and centred mean $\mu_{j}$ and difference $\delta_{j}$ operators on the macroscale grid variables $U_{j}$ (table 1), the learnt subgrid structures are, when retaining terms up to quadratic order in $\gamma$ ,


	$\displaystyle u_{j}$	$\displaystyle=U_{j}+\gamma(\xi-\tfrac{1}{2})\mu_{j}\delta_{j}U_{j}+\tfrac{1}{2}\gamma^{2}(\xi^{2}-\xi+\tfrac{1}{6})\delta_{j}^{2}U_{j}$
		$\displaystyle\quad{}-\tfrac{1}{4}\gamma^{2}(1+\theta^{2})(\xi-\tfrac{1}{2})\mu_{j}\delta_{j}^{3}U_{j}+\tfrac{1}{8}\gamma^{2}(1-\theta^{2})(\xi^{2}-\xi+\tfrac{1}{6})\delta_{j}^{4}U_{j}$
		$\displaystyle\quad{}+\theta(\xi-\tfrac{1}{2})\big{[}-\tfrac{1}{2}(\gamma-\gamma^{2})\delta_{j}^{2}U_{j}+\tfrac{1}{4}\gamma^{2}\delta_{j}^{4}U_{j}\big{]}+\mathcal{O}\mathchoice{\big{(}\gamma^{3}\big{)}}{\big{(}\gamma^{3}\big{)}}{(\gamma^{3})}{(\gamma^{3})}.$	(9a)
The corresponding learnt lowest-order evolution equation, the closure eq. 3, is

	$\displaystyle\partial_{t}U_{j}$	$\displaystyle=\frac{\gamma^{2}}{H^{2}}\left[(1-\theta^{2})\mu_{j}^{2}+\theta^{2}\right]\delta_{j}^{2}U_{j}+\mathcal{O}\mathchoice{\big{(}\gamma^{3}\big{)}}{\big{(}\gamma^{3}\big{)}}{(\gamma^{3})}{(\gamma^{3})}.$	(9b)
This is a correct leading discretisation of the original system for every $\theta$ —correct as it is consistent with the homogeneous diffusion pde to errors $\mathcal{O}\mathchoice{\big{(}H^{2}\big{)}}{\big{(}H^{2}\big{)}}{(H^{2})}{(H^{2})}$ as $H\to 0$ . But (9b) has two unusual features for holistic discretisation: it only appears at $\mathcal{O}\mathchoice{\big{(}\gamma^{2}\big{)}}{\big{(}\gamma^{2}\big{)}}{(\gamma^{2})}{(\gamma^{2})}$ ; and for tuning $\theta\neq\pm 1$ it invokes the next-nearest neighbours, $U_{j\pm 2}$ (through $\mu_{j}^{2}\delta_{j}^{2}$ ), as well as the two nearest neighbours, $U_{j\pm 1}$ .

Table 1: Useful operator identities [npl61, p.65, e.g.] where

E^{-1}

and

E

represent a shift to the left and right, respectively, with subscript

x

indicating a shift of distance

H

, subscript

j

indicating a shift of one in this element index, and subscript

\xi

indicating a shift of one in this dimensionless subgrid variable. For example,

\delta_{j}^{2}U_{j}=(E_{j}-2+E_{j}^{-1})U_{j}=U_{j+1}-2U_{j}+U_{j-1}

and

\mu_{j}\delta_{j}U_{j}=\frac{1}{2}(E_{j}-E_{j}^{-1})U_{j}=\frac{1}{2}[U_{j+1}-U_{j-1}]

\begin{array}[]{rlrl}\hline\cr&\delta=E^{1/2}-E^{-1/2^{\vphantom{|}}}&&\mu=\tfrac{1}{2}(E^{1/2}+E^{-1/2})\\ &\delta_{x}=2\sinh(H\partial_{x}/2)&&\mu_{x}=\cosh(H\partial_{x}/2)\\ &\delta_{j}=2\sinh(\partial_{j}/2)&&\mu_{j}=\cosh(\partial_{j}/2)\\ &E^{\pm 1}=1\pm\mu\delta+\tfrac{1}{2}\delta^{2}&&\mu^{2}=1+\tfrac{1}{4}\delta^{2}\\[2.15277pt] \hline\cr\end{array}

section 2.2 details how centre manifold theory supports such truncations as an approximation to an in-principle exact slow manifold that exists, and is exponentially quickly attractive for a wide range of initial conditions. This theoretical support also applies to nonlinear modifications of such linear diffusion.

The algebra of appendix A computes to arbitrarily high-order in coupling parameter $\gamma$ , learning more and more subgrid structures of the diffusion pde $u_{t}=u_{xx}$ . Computed to element interaction errors $\mathcal{O}\mathchoice{\big{(}\gamma^{6}\big{)}}{\big{(}\gamma^{6}\big{)}}{(\gamma^{6})}{(\gamma^{6})}$ we find the macroscale closure eq. 3 is

$\displaystyle H^{2}\partial_{t}U_{j}$	$\displaystyle=\theta^{2}\big{[}\gamma^{2}\delta_{j}^{2}-\tfrac{1}{2}\gamma^{3}\delta_{j}^{4}+\tfrac{1}{4}\gamma^{4}(\tfrac{5}{3}\delta_{j}^{4}+\delta_{j}^{6})-\tfrac{1}{4}\gamma^{5}(\tfrac{5}{3}\delta_{j}^{6}+\tfrac{1}{2}\delta_{j}^{8})\big{]}U_{j}$
	$\displaystyle\quad{}+\theta^{2}(1-\theta^{2})\big{[}\tfrac{1}{16}\gamma^{4}(\tfrac{2}{3}+\tfrac{1}{3}\mu_{j}^{2})\delta_{j}^{6}-\tfrac{1}{16}\gamma^{5}(\tfrac{2}{3}+\tfrac{1}{3}\mu_{j}^{2})\delta_{j}^{8}\big{]}U_{j}$
	$\displaystyle\quad{}+(1-\theta^{2})\big{[}\gamma^{2}\mu_{j}^{2}\delta_{j}^{2}-\tfrac{1}{2}\gamma^{3}\mu_{j}^{2}\delta_{j}^{4}+\tfrac{1}{6}\gamma^{4}(1+\tfrac{11}{8}\delta_{j}^{2})\mu_{j}^{2}\delta_{j}^{4}$
	$\displaystyle\qquad{}-\tfrac{1}{6}\gamma^{5}(1+\tfrac{5}{8}\delta_{j}^{2})\mu_{j}^{2}\delta_{j}^{6}\big{]}U_{j}+\mathcal{O}\mathchoice{\big{(}\gamma^{6}\big{)}}{\big{(}\gamma^{6}\big{)}}{(\gamma^{6})}{(\gamma^{6})}.$	(9c)

We test how well such discretisations predict the diffusion dynamics by comparing the pde $u_{t}=u_{xx}$ with the equivalent pde of discretisation eq. 9c. To determine the equivalent pde, the macroscale grid fields $U_{j}$ over all discrete $j$ are interpolated to provide a description of $U_{j}$ over the continuous variable $x$ , so that one step in index $j$ of field $U_{j}$ is the same as one step in $x$ of size $H$ : $\delta_{j}^{2}U_{j}=\delta_{x}^{2}U_{j}$ and $\mu_{j}\delta_{j}U_{j}=\mu_{x}\delta_{x}U_{j}$ , and similarly for higher orders.⁴⁴4Mean operators $\mu_{j}$ and $\mu_{x}$ and difference operators $\delta_{j}$ and $\delta_{x}$ are not equivalent; but they give the same results when operating on macroscale fields $U_{j}$ because the $x$ dependence of the macroscale field is defined by an interpolation of the discrete $j$ field values (remark 5). In contrast, as we show in sections 3.2 and 3.3 (in particular sections 3.2 and 3.3), operators $\delta_{x}$ and $\delta_{j}$ (and similarly $\mu_{x}$ and $\mu_{j}$ ) are distinctly different when operating on microscale fields $u_{j}(t,x)$ . Upon computing to higher-order errors of $\mathcal{O}\mathchoice{\big{(}\gamma^{9}\big{)}}{\big{(}\gamma^{9}\big{)}}{(\gamma^{9})}{(\gamma^{9})}$ , the equivalent pde to the discrete model is, by expanding the operator identity that $\delta_{x}=2\sinh(H\partial_{x}/2)$ (table 1),

$\displaystyle\partial_{t}U$	$\displaystyle=\gamma^{2}\partial_{x}^{2}U+(1-\theta^{2})\left[(\tfrac{1}{3}\gamma^{2}-\tfrac{1}{2}\gamma^{3}+\tfrac{1}{6}\gamma^{4})H^{2}\partial_{x}^{4}U\right.$
	$\displaystyle\quad\left.{}+(\tfrac{2}{45}\gamma^{2}-\tfrac{5}{24}\gamma^{3}+\tfrac{43}{144}\gamma^{4}-\tfrac{1}{6}\gamma^{5}+\tfrac{23}{720}\gamma^{6})H^{4}\partial_{x}^{6}U\right.$
	$\displaystyle\quad\left.{}+(\tfrac{1}{315}-\tfrac{3}{80}\gamma^{3}+\tfrac{61}{480}\gamma^{4}-\tfrac{3}{16}\gamma^{5}+\tfrac{49}{360}\gamma^{6}-\tfrac{23}{480}\gamma^{7}+\tfrac{11}{1680})H^{6}\partial_{x}^{8}U\right]$
	$\displaystyle\quad{}+\theta^{2}\left[(\tfrac{1}{12}\gamma^{2}-\tfrac{1}{2}\gamma^{3}+\tfrac{5}{12}\gamma^{4})H^{2}\partial_{x}^{4}U\right.$
	$\displaystyle\quad\left.{}+(\tfrac{1}{360}\gamma^{2}-\tfrac{1}{12}\gamma^{3}+\tfrac{23}{72}\gamma^{4}-\tfrac{5}{12}\gamma^{5}+\tfrac{8}{45}\gamma^{6})H^{4}\partial_{x}^{6}U\right.$
	$\displaystyle\quad\left.{}+(\tfrac{1}{20160}\gamma^{2}-\tfrac{1}{160}\gamma^{3}+\tfrac{13}{192}\gamma^{4}-\tfrac{11}{48}\gamma^{5}+\tfrac{257}{720}\gamma^{6}-\tfrac{4}{15}\gamma^{7}+\tfrac{13}{168}\gamma^{8})H^{6}\partial_{x}^{8}U\right]$
	$\displaystyle\quad{}+\theta^{2}(1-\theta^{2})\left[(\tfrac{1}{16}\gamma^{4}-\tfrac{1}{16}\gamma^{6})H^{4}\partial_{x}^{6}U\right.$
	$\displaystyle\quad\left.{}+(\tfrac{1}{48}\gamma^{4}-\tfrac{1}{16}\gamma^{5}-\tfrac{1}{192}\gamma^{6}+\tfrac{3}{32}\gamma^{7}-\tfrac{3}{64}\gamma^{8})H^{6}\partial_{x}^{8}U\right]$
	$\displaystyle\quad{}+\theta^{4}(1-\theta^{2})\left[(\tfrac{1}{64}\gamma^{6}-\tfrac{1}{64}\gamma^{8})H^{6}\partial_{x}^{8}U\right]+\mathcal{O}\mathchoice{\big{(}H^{8}\partial_{x}^{9}U\big{)}}{\big{(}H^{8}\partial_{x}^{9}U\big{)}}{(H^{8}\partial_{x}^{9}U)}{(H^{8}\partial_{x}^{9}U)}.$	(9d)

For full coupling $\gamma=1$ and for every tuning $\theta$ , all these coefficients evaluate to zero, except for the leading diffusion coefficient which evaluates to one. Thus the discretisation (9c) is consistent with $u_{t}=u_{xx}$ to error $\mathcal{O}\mathchoice{\big{(}H^{8}\partial_{x}^{9}U\big{)}}{\big{(}H^{8}\partial_{x}^{9}U\big{)}}{(H^{8}\partial_{x}^{9}U)}{(H^{8}\partial_{x}^{9}U)}$ .

Inspection of eq. 9d suggests, and section 3.3 proves, that a construction of the slow manifold to errors $o\big{(}\gamma^{p}\big{)}$ , for every even $p$ , learns a discrete model consistent with the diffusion pde to errors $\mathcal{O}\mathchoice{\big{(}H^{p}\partial_{x}^{p+2}U\big{)}}{\big{(}H^{p}\partial_{x}^{p+2}U\big{)}}{(H^{p}\partial_{x}^{p+2}U)}{(H^{p}\partial_{x}^{p+2}U)}$ . That is, the learnt holistic discretisation is consistent with the diffusion pde to arbitrarily high-order.

In application, this systematic consistency means that one controls and estimates errors in the holistic discretisation by varying the order $p$ of the inter-element coupling.

2.2 Support spatial discretisation with centre manifold theory

We use centre manifold theory [[, e.g.,]]Carr81, Haragus2011, Hochs2019 to establish (theorem 5) the existence and emergence of an exact closure (3) to the spatial discretisation of pdes in the class eq. 1. Centre manifold theory is based upon either an equilibrium or subspace/manifold of equilibria, and primarily follows from the local persistence of a spectral gap in the linearised dynamics [[, e.g.,]Part V]Roberts2014a.

To find useful equilibria we embed the pde eq. 1 into a wider class of problems characterised by parameter $\gamma$ . For definiteness of theory, set the macroscale boundary conditions on $\mathbb{X}$ to be that the field $u(t,x)$ is $\mathfrak{L}$ -periodic in space $x$ . Recall that $u_{j}(t,x)$ denotes solutions of the pde eq. 1 on the elements $\mathbb{X}_{j}$ that partition the domain, so we reserve $u(t,x)$ , over $\tilde{}\mathbb{X}$ , to denote the union over all elements of such solutions. Then the edge-element conditions eq. 6 couple each element together. The information flow between elements is moderated by the coupling parameter $\gamma$ . The parameter $\gamma$ connects the original pde over the whole domain, $\gamma=1$ , to a useful base problem, $\gamma=0$ .

Lemma 2 (equilibria).

The pde eq. 1 with section 1 on elements $\{\mathbb{X}_{j}\}$ with edge conditions eq. 6 possesses an $N$ -dimensional subspace $\mathbb{E}$ of equilibria parametrised by ${\text{\boldmath$U$}}=(U_{1},U_{2},\ldots,U_{N})$ with parameter values $\alpha=\gamma=0$ . Each equilibrium is of a piecewise constant field $u^{*}(x)$ such that for every $j\in\mathbb{J}$ , $u^{*}(x)=U_{j}$ for $x\in\mathbb{X}_{j}$ .

There usually are other equilibria of eq. 1+eq. 6: this lemma only addresses the subspace $\mathbb{E}$ .

Proof.

With advection-reaction coefficient $\alpha=0$ the pde eq. 1 takes the form $u_{t}=-f(x,u,u_{x})_{x}$ . For every piecewise constant field $u^{*}(x)=U_{j}$ in $\mathbb{X}_{j}$ the gradient $u^{*}_{x}=0$ , and by section 1 the flux $f(x,u^{*},u_{x}^{*})=f(x,U_{j},0)=0$ . Thus ${\text{\boldmath$U$}}=(U_{1},U_{2},\ldots,U_{N})$ with $\alpha=0$ are equilibria of the pde on $\tilde{}\mathbb{X}$ . The equilibria also require coupling $\gamma=0$ , for which the right-hand sides of the conditions eq. 6 are zero. The left-hand sides are also zero for fields constant in each element, and so the conditions eq. 6 are satisfied. ∎

We seek solutions $u=u^{*}(x)+\hat{u}(t,x)$ of the general pde eq. 1 where $\hat{u}(t,x)$ denotes a small linearised perturbation of the equilubria in $\mathbb{E}$ . Use $\hat{u}_{j}(t,x)$ as a synonym for $\hat{u}(t,x)$ on the $j$ th element $\mathbb{X}_{j}$ . Then, since parameter $\alpha=0$ and for small $\hat{u}$ flux $f(x,u,u_{x})\approx f(x,u^{*},\hat{u}_{x})$ , the pde eq. 1 linearises to

	$\displaystyle\hat{u}_{t}=-\partial_{x}\left[\mathchoice{\frac{\partial f}{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}(x,u^{*},0)\hat{u}_{x}\right]\quad\text{on }\mathbb{X};$
	$\displaystyle\text{that is,}\quad\hat{u}_{jt}=\partial_{x}\big{[}\kappa_{j}(x)\hat{u}_{jx}\big{]}\quad\text{for }x\in\mathbb{X}_{j}\,,$		(10)

where we define the local diffusivity $\kappa_{j}(x):=-\mathchoice{\frac{\partial f(x,U_{j},0)}{\partial u_{x}}}{{\partial f(x,U_{j},0)}/{\partial u_{x}}}{{\partial f(x,U_{j},0)}/{\partial u_{x}}}{{\partial f(x,U_{j},0)}/{\partial u_{x}}}$ for $x\in\mathbb{X}_{j}$ . From section 1, $\kappa_{j}(x)\geq\nu>0$ . The edge conditions eq. 6 are linear, so they apply here with flux $f_{j}=-\kappa_{j}\hat{u}_{jx}$ , and $u_{j}$ replaced by $\hat{u}_{j}$ —for symbolic simplicity let’s omit the hat hereafter.

A spectrum of eigenvalues arises from the general linearised dynamics (10). We turn to determining this spectrum when the coupling parameter $\gamma=0$ , namely the diffusion eq. 10 with isolating edge conditions eq. 7. In doing this, the following self-adjointness is crucial.

Lemma 3 (self-adjoint).

The diffusion operator

\mathcal{L}:=[\kappa_{j}(x)\partial_{x}\cdot]_{x}

(11)

subject to the $\gamma=0$ edge conditions eq. 7, is self-adjoint in the Hilbert space $\mathbb{L}$ with the usual inner product $\left<v,u\right>:=\sum_{j}\int_{\mathbb{X}_{j}}v_{j}u_{j}\,dx$ .

Proof.

Proceed straightforwardly via integration by parts, where for every $u\in\mathbb{L}$ the flux $f_{j}:=-\kappa_{j}(x)u_{jx}$ , and likewise for $\tilde{u}$ and $\tilde{f}$ :

	$\displaystyle\left<\tilde{u},\mathcal{L}u\right>-\left<\mathcal{L}\tilde{u},u\right>$	$\displaystyle=\sum_{j}\int_{\mathbb{X}_{j}}\tilde{u}_{j}[\kappa_{j}(x)u_{jx}]_{x}-u_{j}[\kappa_{j}(x)\tilde{u}_{jx}]_{x}\,dx$
		$\displaystyle=\sum_{j}\big{[}\tilde{u}_{j}\kappa_{j}(x)u_{jx}-u_{j}\kappa_{j}(x)\tilde{u}_{jx}\big{]}_{L_{j}}^{R_{j}}$
		$\displaystyle=\sum_{j}\big{(}-\tilde{u}_{j}^{R}f_{j}^{R}+u_{j}^{R}\tilde{f}_{j}^{R}+\tilde{u}_{j}^{L}f_{j}^{L}-u_{j}^{L}\tilde{f}_{j}^{L}\big{)}$
		$\displaystyle=\sum_{j}\big{(}-\tilde{u}_{j}^{L}f_{j}^{L}+u_{j}^{L}\tilde{f}_{j}^{L}+\tilde{u}_{j}^{L}f_{j}^{L}-u_{j}^{L}\tilde{f}_{j}^{L}\big{)}\quad(\text{by }\lx@cref{creftype~refnum}{eqsiecc0})$
		$\displaystyle=0\,.$

∎

section 2.4 extends section 2.2 by establishing that the linear operator $\mathcal{L}$ , with inter-element coupling controlled by edge conditions eq. 6 is self-adjoint for every $\gamma$ .

By the self-adjointness of $\mathcal{L}$ (11), we need only seek real eigenvalues in the spectrum of (10). Let’s first discuss the zero eigenvalues, and second, the non-zero eigenvalues. For the diffusion operator eq. 11 any eigenfunction corresponding to an eigenvalue of zero must be linear on each $\mathbb{X}_{j}$ . The element-periodicity conditions eq. 7 for $\gamma=0$ then require that the eigenfunctions must be constant within every element. That is, the zero-eigenspace of $\mathcal{L}$ when $\gamma=0$ is the $N$ -D subspace of equilibria $\mathbb{E}$ (defined in section 2.2). By self-adjointness, there are no generalised eigenfunctions. Hence the operator eq. 11+eq. 7 has $N$ eigenvalues of zero, and $\mathbb{E}$ is the corresponding $N$ -D slow subspace.

Lemma 4 (exponential dichotomy).

Under section 1, the operator eq. 11, $\mathcal{L}=[\kappa_{j}(x)\partial_{x}\cdot]_{x}$ subject to eq. 7, in $\mathbb{L}$ has $N$ zero eigenvalues and all other eigenvalues $\lambda$ are negative and bounded away from zero by $\lambda\leq-\beta$ where a bound $\beta:=\min_{j\in\mathbb{J},\,x\in\mathbb{X}_{j}}\big{[}4\pi^{2}\kappa_{j}(x)/H_{j}^{2}\big{]}>0$ .

Proof.

The precisely $N$ zero eigenvalues are established in the paragraphs preceding section 2.2. Also, section 2.2 establishes all eigenvalues of $\mathcal{L}$ are real. Let $\lambda$ be a non-zero eigenvalue and $v(x)$ be a corresponding eigenfunction. Consequently, and using the usual inner product $\left<v,u\right>:=\sum_{j}\int_{\mathbb{X}_{j}}v_{j}u_{j}\,dx$ , we proceed to derive the inequality

$\displaystyle(-\lambda)\\|v\\|^{2}$	$\displaystyle=-\lambda\left<v,v\right>=\left<v,-\lambda v\right>=\left<v,-\mathcal{L}v\right>=\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}-v_{j}[\kappa_{j}(x)v_{jx}]_{x}\,dx$
	$\displaystyle=\sum_{j\in\mathbb{J}}\big{[}-v_{j}\kappa_{j}v_{jx}\big{]}_{L_{j}}^{R_{j}}+\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}v_{jx}\kappa_{j}v_{jx}\,dx\quad(\text{integrating by parts})$
	$\displaystyle=\sum_{j\in\mathbb{J}}0+\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}\kappa_{j}v_{jx}^{2}\,dx\quad(\text{using \lx@cref{creftype~refnum}{eqsiecc0}})$
	$\displaystyle\geq\frac{\beta}{4\pi^{2}}\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}H_{j}^{2}v_{jx}^{2}\,dx\,,$	(12)

for the bound $\beta$ defined in terms of the lowest diffusivity of operator $\mathcal{L}$ : $\kappa_{j}(x)\geq\beta H_{j}^{2}/(4\pi^{2})$ for all $x\in\mathbb{X}_{j}$ , for every $j\in\mathbb{J}$ . The bound $\beta>0$ follows from the monotonicity of $f$ , section 1, since $\kappa_{j}(x)=-\mathchoice{\frac{\partial f}{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}{{\partial f}/{\partial u_{x}}}(x,U_{j},0)\geq\nu>0$ for all $j\in\mathbb{J}$ , so $\beta/(4\pi^{2})=\min_{j\in\mathbb{J},\,x\in\mathbb{X}_{j}}\big{[}\kappa_{j}(x)/H_{j}^{2}\big{]}\geq\nu/{\max_{j}H_{j}^{2}}>0$ . A first consequence of inequality (12) is that there are no positive eigenvalues $\lambda$ .

Secondly, we establish that the non-zero eigenvalues are bounded away from zero by $-\beta$ . This bound is established by relating inequality (12) for a general eigenvalue $\lambda$ and associated eigenvector $v$ of some heterogeneous problem to the known lowest magnitude eigenvalue $\lambda_{1}$ of a spatially homogeneous problem. Let $\mathcal{L}_{*}:=H_{j}^{2}\mathchoice{\frac{\partial^{2}}{\partial x^{2}}}{{\partial^{2}}/{\partial x^{2}}}{{\partial^{2}}/{\partial x^{2}}}{{\partial^{2}}/{\partial x^{2}}}$ for $x\in\mathbb{X}_{j}$ with edge conditions eq. 7 (that is, the linear operator (11) for the special homogeneous case of $\kappa_{j}(x)=H_{j}^{2}$ ). Then by the reverse argument to that of the previous paragraph, $\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}H_{j}^{2}v_{x}^{2}\,dx=\cdots=\left<v,-\mathcal{L}_{*}v\right>$ . By the Rayleigh–Ritz theorem, the smallest magnitude, non-zero, eigenvalue $\lambda_{1}$ of $\mathcal{L}_{*}$ satisfies $-\lambda_{1}=\min_{w\perp\mathbb{E}}\left<w,-\mathcal{L}_{*}w\right>/\|w\|^{2}$ , and so $\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}H_{j}^{2}v_{x}^{2}\,dx=\left<v,-\mathcal{L}_{*}v\right>\geq-\lambda_{1}\|v\|^{2}$ . But for $\mathcal{L}_{*}$ , subject to the condition of $H_{j}$ -periodicity eq. 7, it is well-known that the smallest magnitude, non-zero eigenvalue is $\lambda_{1}=-4\pi^{2}$ (corresponding to eigenmodes $\operatorname{e}^{\pm\operatorname{i}2\pi x/H_{j}}$ ). Hence, inequality (12) becomes $(-\lambda)\|v\|^{2}\geq\frac{\beta}{4\pi^{2}}\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}H_{j}^{2}v_{x}^{2}\,dx\geq-\frac{\beta}{4\pi^{2}}\lambda_{1}\|v\|^{2}=\beta\|v\|^{2}$ . That is, $\lambda\leq-\beta$ for every non-zero eigenvalue $\lambda$ . ∎

With section 2.2 proving a spectral gap in the linearised dynamics of eq. 10+eq. 7 with $\alpha=\gamma=0$ , we are now in a position to discuss the corresponding slow manifold of eq. 1+eq. 6. Inspired by the backwards theory of numerical analysis [[, e.g.,]]Grcar2011, we apply backwards theory to establish the existence and emergence of a slow manifold [Hochs2019, Roberts2018a].

For rigorous theory we notionally adjoin the two trivial dynamical equations $\alpha_{t}=\gamma_{t}=0$ to the linearised system eq. 10+eq. 6. Then the equilibria in $\alpha\gamma u$ -space are $(0,0,u^{*}(x))$ since these equilibria are at $\alpha=\gamma=0$ . Hence, strictly, there are two extra zero eigenvalues associated with the trivial $\alpha_{t}=\gamma_{t}=0$ , and the corresponding slow subspace of each equilibria is $(N+2)$ -D. But, except for issues associated with the domain of validity for non-zero $\alpha$ and $\gamma$ , for simplicity, in the following we do not explicitly discuss these two trivial dynamical equations nor their eigenvalues, but consider them implicit.

Theorem 5 (slow manifold).

Consider the generic reaction-advection-diffusion pde eq. 1 with section 1 on domain $\tilde{}\mathbb{X}$ with inter-element edge conditions eq. 6. For every order $p\geq 2$ , there exists a system for $u(t,x)$ , namely the following polynomial map $\upsilon_{p}$ and pde with polynomial nonlinearity $\mathfrak{f}_{p}$

u(t,x)=u^{*}+\upsilon_{p}(t,x,v(t,x))\quad\text{where }v_{t}=\mathcal{L}v+\mathfrak{f}_{p}(t,x,v),

(13)

such that the map+pde eq. 13 is $\mathcal{O}\mathchoice{\big{(}\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p}\big{)}}{\big{(}\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p}\big{)}}{(\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p})}{(\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p})}$ close to eq. 1+eq. 6, and that the map+pde eq. 13 has the following properties.

5(a)

There exists an open domain $\mathcal{E}$ , in $\alpha\gamma u$ -space, containing the subspace $\mathbb{E}$ , such that in $\mathcal{E}$ there exists a constructible slow manifold of eq. 13, polynomial in $U$ , $\alpha$ and $\gamma$ of order $p-1$ :

u=u(t,x,{\text{\boldmath$U$}},\gamma,\alpha)\quad\text{such that}\quad\frac{d{\text{\boldmath$U$}}}{dt}={\text{\boldmath$G$}}(t,{\text{\boldmath$U$}},\gamma,\alpha).

(14)

5(b)

This slow manifold is emergent in the sense that for all solutions $u(t,x)$ of eq. 13 that stay in $\mathcal{E}$ , there exists a solution ${\text{\boldmath$U$}}(t)$ of eq. 14 such that $u(t,x)=u(t,x,{\text{\boldmath$U$}}(t),\gamma,\alpha)+\mathcal{O}\mathchoice{\big{(}e^{-\beta^{\prime}t}\big{)}}{\big{(}e^{-\beta^{\prime}t}\big{)}}{(e^{-\beta^{\prime}t})}{(e^{-\beta^{\prime}t})}$ for some decay rate $0<\beta^{\prime}<\beta$ ( $\beta$ defined by section 2.2).

In this theorem, and subsequently, asserting a quantity $\epsilon=\mathcal{O}\mathchoice{\big{(}v^{p}+\alpha^{q}+\gamma^{r}\big{)}}{\big{(}v^{p}+\alpha^{q}+\gamma^{r}\big{)}}{(v^{p}+\alpha^{q}+\gamma^{r})}{(v^{p}+\alpha^{q}+\gamma^{r})}$ is to mean the asymptotic property that $\epsilon/(v^{p}+\alpha^{q}+\gamma^{r})$ is bounded as $(v,\alpha,\gamma)\to\text{\boldmath$0$}$ [[, e.g.,]]Roberts2014a. Also, the norm $\|u-\mathbb{E}\|$ is a measure of the distance of $u$ from the subspace $\mathbb{E}$ of equilibria. Strictly, this measure of distance is a norm in appropriate graded Frechet spaces that are obtained from intersections of compactly nested sequences of Banach spaces [Hochs2019, Defs. 2.2–2.4] (hereafter denoted HR).

Proof.

For every equilibria $u^{*}\in\mathbb{E}$ , we invoke the backwards Theorems 2.18 and 2.22 of HR. From section 2.2, the eigenfunctions of the operator $\mathcal{L}=[\kappa_{j}(x)\partial_{x}\cdot]_{x}$ subject to eq. 7 form an orthonormal basis of the requisite Hilbert space on $\tilde{}\mathbb{X}:=\cup_{j}\mathbb{X}_{j}$ (Hilbert–Schmidt theorem). Since $f$ and $g$ are smooth (infinitely differentiable) they can be expanded polynomially to any specified order $p$ as required. Hence, with section 1, the conditions of Theorems 2.18 and 2.22 by HR are satisfied. Thus Corollary 2.23 by HR applies.

For every order $p$ , and for $u(t,x)$ defined by the map+pde (13), this corollary asserts that $u(t,x)$ satisfies eq. 1+eq. 6 to a residual $R_{p}=\mathcal{O}\mathchoice{\big{(}\|(v,\alpha,\gamma)\|^{p}\big{)}}{\big{(}\|(v,\alpha,\gamma)\|^{p}\big{)}}{(\|(v,\alpha,\gamma)\|^{p})}{(\|(v,\alpha,\gamma)\|^{p})}=\mathcal{O}\mathchoice{\big{(}\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p}\big{)}}{\big{(}\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p}\big{)}}{(\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p})}{(\|u-\mathbb{E}\|^{p}+\alpha^{p}+\gamma^{p})}$ . In this sense of asymptotic agreement, every such system eq. 13 is close to the system eq. 1+eq. 6.

Also, every system eq. 13 is constructed to have clear invariant subspaces in $v$ , and so clear invariant manifolds in $u$ (Definition 2.11 of HR). Further, from the exponential dichotomy of section 2.2 there exists a domain $D_{\mu}$ , in $\alpha\gamma u$ -space, in which the system eq. 13 has a slow manifold. Since $0\leq\mu<\beta\leq-\lambda$ (Section 2.3.1 of HR) and by Proposition 2.10 of HR, the slow manifold solutions are exponentially quickly attractive, at rate of at least $\beta^{\prime}=\beta-\mu>0$ , to all solutions in the domain $D_{\mu}$ .

These properties hold when based upon each and every equilibrium $u^{*}$ in the slow subspace $\mathbb{E}$ . By smoothness of the system eq. 1+eq. 6, we can ensure we construct the close systems eq. 13 to be smooth as a function of $u^{*}$ . Hence the properties hold smoothly in the union of all the domains $D_{\mu}$ , to form a global domain, denoted $\mathcal{E}$ , in which the properties hold. This establishes theorem 5. ∎

theorem 5 holds for the pde system eq. 1 on domain $\tilde{}\mathbb{X}$ with edge conditions eq. 6 between elements. The system eq. 1+eq. 6 reverts to the original pde eq. 1 on the domain $\mathbb{X}$ when we evaluate eq. 1+eq. 6 on $\tilde{}\mathbb{X}$ at full coupling $\gamma=1$ . Thus the evolution equation eq. 14 evaluated at full coupling, namely $d{\text{\boldmath$U$}}/dt={\text{\boldmath$G$}}(t,{\text{\boldmath$U$}},1,\alpha)$ , is controllably close to a discretisation of the original pde. Further, this evolution is that on a slow manifold which is emergent in the domain $\mathcal{E}_{1}:=\mathcal{E}|_{\gamma=1}$ (which is not empty as $\mathcal{E}_{1}$ contains at least the equilibria $u={}$ constant and $\alpha=0$ ).

Similar slow manifold properties to those of theorem 5 could be established via the forward Theorems 2.9 and 3.22 of [Haragus2011], together with Proposition 3.6 of [Potzsche2006]. The argument for Theorem 6 by [Jarrad2016a] is an example. However, such a forward approach attempts to prove the existence of an invariant manifold for the nonlinear system eq. 1+eq. 6, and in doing so imposes significant constraints on the functions $f$ and $g$ , constraints that in applications either are often hard to establish, or are not satisfied. The backwards theory of [Hochs2019], which shows the existence of a system with an invariant manifold arbitrarily close to eq. 1+eq. 6, is less restrictive on the functions $f$ and $g$ and so has a wider range of applications.

2.3 The slow manifold of wave-like PDEs

Although this article’s principal scope is the spatial discretisation, or dimensional reduction, of reaction-advection-diffusion pdes eq. 1, much of the approach and theory usefully applies to the spatial discretisation of non-autonomous wave-like pdes in the form

u_{tt}=-f(x,u,u_{x})_{x}+\alpha g(t,x,u,u_{x})

(15)

on a domain $\mathbb{X}$ , with section 1. The only difference between eqs. 15 and 1 is the number of time derivatives on the left-hand side. This subsection comments on the similarities and differences of the theoretical support for such wave systems.

Partition space into $N$ elements as described at the beginning of section 2 with the $j$ th element over the interval $\mathbb{X}_{j}=(L_{j},R_{j})$ , and apply the edge conditions eq. 6. Then, for $\alpha=\gamma=0$ , the subspace $\mathbb{E}$ of piecewise linear equilibria of section 2.2 still exists. Upon linearisation of eq. 15 about each of these equilibria, the spatial differential operator $\mathcal{L}=[\kappa_{j}(x)\partial_{x}\cdot]_{x}$ on $\tilde{}\mathbb{X}=\cup_{j}\mathbb{X}_{j}$ remains self-adjoint. The exponential dichotomy of the operator $\mathcal{L}$ (section 2.2) still applies, namely that there are $N$ eigenvalues of zero, and the others are ${}\leq-\beta$ . So far, the considerations are the same for both diffusion systems and wave systems.

Differences between wave-like pdes eq. 15 and reaction-advection-diffusion pdes eq. 1 start with the eigenvalues of the right-hand side operator $\mathcal{L}$ with its eigenvalues denoted by $\lambda$ . Seeking wave solutions of (15) with frequency $\omega$ , then $\omega=\pm\sqrt{-\lambda}$ and all frequencies $\omega$ are real as every $\lambda\leq 0$ . Here the slow manifold dichotomy is now between slow waves with near zero frequency, separated from fast waves with frequencies ${}\geq\sqrt{\beta}$ (section 2.2). Such subcentre slow manifolds (a class of nonlinear normal modes [Shaw94a, e.g.]) are ubiquitous in both geophysical and elastic engineering applications. However, much less is known rigorously about subcentre slow manifolds: even their existence is problematic [Lorenz87]. Nonetheless, inspired by the backwards theory of [Roberts2018a, Hochs2019] we make the following backwards conjecture for the wave pde eq. 15 that is analogous to theorem 5 for dissipative systems.

Conjecture 6.

For every order $p$ , there exists a (polynomial) coordinate transformation and a (polynomial) pde system in the new variables $({\text{\boldmath$U$}},V(x))$ of the form


	$\displaystyle u=u(t,x,{\text{\boldmath$U$}},V,\gamma,\alpha),\quad$		(16a)
	$\displaystyle\mathchoice{\frac{d^{2}{\text{\boldmath$U$}}}{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}={\text{\boldmath$G$}}(t,{\text{\boldmath$U$}},V,\gamma,\alpha),\quad\mathchoice{\frac{\partial^{2}V}{\partial t^{2}}}{{\partial^{2}V}/{\partial t^{2}}}{{\partial^{2}V}/{\partial t^{2}}}{{\partial^{2}V}/{\partial t^{2}}}=H(t,x,{\text{\boldmath$U$}},V,\gamma,\alpha)V,$		(16b)

such that in the $\alpha\gamma u$ -space the system (16) is the same as the pde eq. 15 to a residual $\mathcal{O}\mathchoice{\big{(}\|u-\mathbb{E}\|^{p}+\gamma^{p}+\alpha^{p}\big{)}}{\big{(}\|u-\mathbb{E}\|^{p}+\gamma^{p}+\alpha^{p}\big{)}}{(\|u-\mathbb{E}\|^{p}+\gamma^{p}+\alpha^{p})}{(\|u-\mathbb{E}\|^{p}+\gamma^{p}+\alpha^{p})}$ , and such that the subspace $\mathbb{E}$ is the tangent space of the manifold $V=0$ at $\gamma=\alpha=0$ .

2.3(a)

Let $\mathcal{E}$ denote a $\alpha\gamma u$ -domain in which the coordinate transform eq. 16a is a diffeomorphism containing $\mathbb{E}$ , then $V=0$ is an exact slow manifold of the dynamics of eq. 16: that is, a slow manifold is

u=u(t,x,{\text{\boldmath$U$}},0,\gamma,\alpha)\quad\text{such that}\quad\mathchoice{\frac{d^{2}{\text{\boldmath$U$}}}{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}={\text{\boldmath$G$}}(t,{\text{\boldmath$U$}},0,\gamma,\alpha).

(17)

2.3(b)

Solutions near, but off the slow manifold with $V\neq 0$ , generally evolve differently (due to wave-wave forcing of mean flow):

\mathchoice{\frac{d^{2}{\text{\boldmath$U$}}}{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}{{d^{2}{\text{\boldmath$U$}}}/{dt^{2}}}={\text{\boldmath$G$}}(t,{\text{\boldmath$U$}},0,\gamma,\alpha)+\mathcal{O}\mathchoice{\big{(}\|V\|^{2}\big{)}}{\big{(}\|V\|^{2}\big{)}}{(\|V\|^{2})}{(\|V\|^{2})}.

Consequently, we contend that the methodology developed here for learning spatially discrete, finite dimensional, accurate, models of dissipative pdes may be also usefully applied to many wave-like pdes eq. 15, provided we accept the potential for microscale wave-wave interactions be unresolved in the model. In order to resolve such wave-wave forcing of the mean flow, one would have to include the amplitude modulation of the fast waves into the list of macroscale variables of interest.

2.4 Preserving self-adjointness

We now return to reaction-advection-diffusion pdes eq. 1 embedded into spatial elements by the coupling conditions eq. 6 (for general $\gamma$ ). theorem 5 proves the existence of systems arbitrarily close to eq. 1+eq. 6—systems that have a constructible slow manifold that forms an emergent discretisation of eq. 1. However, theorem 5 does not address how well the constructed slow manifold discretisation models the original pde eq. 1. This article proves two crucial properties of the slow manifold modelling. Firstly, this section establishes that the particular edge conditions eq. 6 preserve self-adjointness and thus preserves this fundamental symmetry in the discretisation. Secondly, section 3 proves that the slow manifold of eq. 1+eq. 6, the discrete model, is consistent to the pde eq. 1.

Self-adjointness is a property of linear operators, but the pde eq. 1 is nonlinear. Consequently, we establish self-adjointness for linear perturbations about a set of base states of eq. 1, with coefficient $\alpha=0$ to eliminate nonlinearity $g$ from consideration. That is, as before, we linearise pde eq. 1 about the equilibria $u_{j}(t,x)=U_{j}$ constant in each element $\mathbb{X}_{j}$ (section 2.2), to obtain the linearised pde (10). section 2.2 established self-adjointness of the diffusion operator $\mathcal{L}$ of (10) subject to the $\gamma=0$ edge conditions eq. 7, and the following section 2.4 generalises it to edge conditions eq. 6 with $\gamma\neq 0$ .

Lemma 7 (self-adjoint).

The differential operator of the system eqs. 10 and 6, namely $\mathcal{L}=[\kappa_{j}(x)\partial_{x}\cdot]_{x}$ on $\mathbb{L}$ and subject to eq. 6, and in the usual inner product $\left<\tilde{u},u\right>:=\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}\tilde{u}_{j}u_{j}\,dx$ , is self-adjoint for almost every real $\gamma$ and $\theta$ .

The exception in “almost every” is the isolated case when the number of elements is even, $\gamma=1$ , and $\theta=0$ that is associated with the unphysical nature of the zigzag mode as discussed by remark 4.

Proof.

Let’s consider the edge conditions of the $N$ elements as $N\times N$ matrix equations with coupling matrices $C_{\pm}$ . For each of the patch right and left edges, $e\in\{R,L\}$ respectively, define the field vector ${\text{\boldmath$u$}}^{e}:=(u_{1}^{e},\ldots,u_{N}^{e})$ and the flux vector ${\text{\boldmath$f$}}^{e}:=(f_{1}^{e},\ldots,f_{N}^{e})$ . Then, using superscript $T$ to denote matrix transpose, the edge conditions (6) are

	$\displaystyle C_{+}{\text{\boldmath$u$}}^{R}=C_{-}^{T}{\text{\boldmath$u$}}^{L},\quad C_{-}{\text{\boldmath$f$}}^{R}=C_{+}^{T}{\text{\boldmath$f$}}^{L},\quad\text{for circulant matrices}$		(18)
	$\displaystyle C_{\pm}:=\begin{bmatrix}1-\tfrac{\gamma}{2}(1\pm\theta)&0&\cdots&0&\tfrac{\gamma}{2}(1\pm\theta)\\ \tfrac{\gamma}{2}(1\pm\theta)&\ddots&\ddots&&0\\ 0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\ 0&\cdots&0&\tfrac{\gamma}{2}(1\pm\theta)&1-\tfrac{\gamma}{2}(1\pm\theta)\end{bmatrix}.$

We firstly prove that these edge conditions preserve self-adjointness when either $C_{+}$ or $C_{-}$ is invertible, and then secondly when neither is invertible. This proof requires the commutativity

C_{+}^{T}C_{-}=C_{-}C_{+}^{T}.

(19)

This commutativity follows since $C_{\pm}$ , and their transposes, are circulant [Gray2006, Thm. 3.1(1)], and holds for any two circulant matrices.

First consider the case when $C_{+}$ is invertible. Then $C_{+}^{T}$ is also invertible and circulant, so by commutativity (19), $C_{-}C_{+}^{-1,T}=C_{+}^{-1,T}C_{-}$ . Now, for every $u,\tilde{u}\in\mathbb{L}$ ,

	$\displaystyle\left<\tilde{u},\mathcal{L}u\right>-\left<\mathcal{L}\tilde{u},u\right>$
	$\displaystyle=\sum_{j\in\mathbb{J}}\int_{\mathbb{X}_{j}}\tilde{u}_{j}[\kappa_{j}(x)u_{jx}]_{x}-u_{j}[\kappa_{j}(x)\tilde{u}_{jx}]_{x}\,dx$
	$\displaystyle\qquad(\text{then integrate by parts})$
	$\displaystyle=\sum_{j\in\mathbb{J}}\big{[}\tilde{u}_{j}\kappa_{j}(x)u_{jx}-u_{j}\kappa_{j}(x)\tilde{u}_{jx}\big{]}_{L_{j}}^{R_{j}}$
	$\displaystyle\qquad(\text{then as flux }f_{j}:=-\kappa_{j}(x)u_{jx}\text{ and }\tilde{f}_{j}:=-\kappa_{j}(x)\tilde{u}_{jx})$
	$\displaystyle=\sum_{j\in\mathbb{J}}\big{[}-\tilde{u}_{j}^{R}f_{j}^{R}+u_{j}^{R}\tilde{f}_{j}^{R}+\tilde{u}_{j}^{L}f_{j}^{L}-u_{j}^{L}\tilde{f}_{j}^{L}\big{]}$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}{\text{\boldmath$f$}}^{R}+{\text{\boldmath$u$}}^{R,T}\tilde{{\text{\boldmath$f$}}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle\qquad(\text{then, since $C_{+}$ is invertible})$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}C_{+}^{T}C_{+}^{-1,T}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{R,T}C_{+}^{T}C_{+}^{-1,T}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle=-[C_{+}\tilde{{\text{\boldmath$u$}}}^{R}]^{T}C_{+}^{-1,T}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}+[C_{+}{\text{\boldmath$u$}}^{R}]^{T}C_{+}^{-1,T}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle\qquad(\text{apply edge conditions~{}\eqref{eq:LRcc}})$
	$\displaystyle=-[C_{-}^{T}\tilde{{\text{\boldmath$u$}}}^{L}]^{T}C_{+}^{-1,T}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}+[C_{-}^{T}{\text{\boldmath$u$}}^{L}]^{T}C_{+}^{-1,T}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{L,T}C_{-}C_{+}^{-1,T}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{L,T}C_{-}C_{+}^{-1,T}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle\qquad(\text{apply commutativity~{}\eqref{eq:id}})$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{L,T}C_{+}^{-1,T}C_{-}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{L,T}C_{+}^{-1,T}C_{-}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle\qquad(\text{apply edge conditions~{}\eqref{eq:LRcc}})$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{L,T}C_{+}^{-1,T}C_{+}^{T}{\text{\boldmath$f$}}^{L}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{L,T}C_{+}^{-1,T}C_{+}^{T}\tilde{{\text{\boldmath$f$}}}^{L}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle=0\,.$

Thus $\mathcal{L}$ is self-adjoint when $C_{+}$ is invertible.

Second, consider the case when $C_{-}$ is invertible. The detailed argument directly corresponds to the previous case, so we present a summary:

	$\displaystyle\langle\tilde{u},\mathcal{L}u\rangle-\langle\mathcal{L}\tilde{u},u\rangle$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}{\text{\boldmath$f$}}^{R}+{\text{\boldmath$u$}}^{R,T}\tilde{{\text{\boldmath$f$}}}^{R}+\tilde{{\text{\boldmath$u$}}}^{L,T}{\text{\boldmath$f$}}^{L}-{\text{\boldmath$u$}}^{L,T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}{\text{\boldmath$f$}}^{R}+[C_{-}^{T}\tilde{{\text{\boldmath$u$}}}^{L}]^{T}C_{-}^{-1}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{R,T}\tilde{{\text{\boldmath$f$}}}^{R}-[C_{-}^{T}{\text{\boldmath$u$}}^{L,T}]C_{-}^{-1}\tilde{{\text{\boldmath$f$}}}^{L}\$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}{\text{\boldmath$f$}}^{R}+[C_{+}\tilde{{\text{\boldmath$u$}}}^{R}]^{T}C_{-}^{-1}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{R,T}\tilde{{\text{\boldmath$f$}}}^{R}-[C_{+}{\text{\boldmath$u$}}^{R}]^{T}C_{-}^{-1}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{R,T}C_{-}^{-1}C_{+}^{T}{\text{\boldmath$f$}}^{L}+{\text{\boldmath$u$}}^{R,T}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{R,T}C_{-}^{-1}C_{+}^{T}\tilde{{\text{\boldmath$f$}}}^{L}$
	$\displaystyle=-\tilde{{\text{\boldmath$u$}}}^{R,T}{\text{\boldmath$f$}}^{R}+\tilde{{\text{\boldmath$u$}}}^{R,T}C_{-}^{-1}C_{-}{\text{\boldmath$f$}}^{R}+{\text{\boldmath$u$}}^{R,T}\tilde{{\text{\boldmath$f$}}}^{R}-{\text{\boldmath$u$}}^{R,T}C_{-}^{-1}C_{-}\tilde{{\text{\boldmath$f$}}}^{R}$
	$\displaystyle=0\,.$

Thus $\mathcal{L}$ is self-adjoint when $C_{-}$ is invertible.

Now consider the case that neither $C_{+}$ nor $C_{-}$ is invertible, so that their determinants are zero. By a first row expansion, $\det C_{\pm}=\left[1-\tfrac{\gamma}{2}(1\pm\theta)\right]^{N}-(-1)^{N}\left[\tfrac{\gamma}{2}(1\pm\theta)\right]^{N}$ . The determinant is zero only when $\left[1-\tfrac{\gamma}{2}(1\pm\theta)\right]^{N}=\left[-\tfrac{\gamma}{2}(1\pm\theta)\right]^{N}$ . This cannot occur in the case $\gamma(1\pm\theta)=0$ as then the $\textsc{rhs}=0$ and the $\textsc{lhs}=1$ : consequently we can divide by the rhs and seek when $\big{[}1-\frac{2}{\gamma(1\pm\theta)}\big{]}^{N}=1$ . For real $\gamma,\theta$ , the zero determinant thus can only occurs when $N$ is even and $1-\frac{2}{\gamma(1\pm\theta)}=-1$ . That is, when $N$ is even and $\gamma(1\pm\theta)=1$ . Hence the previous arguments which depend on either $C_{+}$ or $C_{-}$ being invertible only fail when $N$ is even and both $\gamma(1\pm\theta)=1$ . This last pair of equations is only satisfied when $\gamma=1$ and $\theta=0$ . Therefore, the operator $\mathcal{L}$ on $\mathbb{L}$ is self-adjoint for every real $\gamma$ and every real $\theta$ , except perhaps for the isolated specific case of even $N$ , $\gamma=1$ and $\theta=0$ . ∎

The proofs of this section hold for periodic boundary conditions for $u(t,x)$ on $\mathbb{X}$ . We expect further research to establish self-adjointness for other common boundary conditions for $u(t,x)$ on $\mathbb{X}$ , such as Dirichlet and Neumann.

3 Prove high-order consistency in 1D space

section 2.2 establishes the in-principle existence and emergence of an exact discrete closure to the dynamics of pdes in the class eq. 1, a closure that also preserves self-adjointness (section 2.4). An issue is whether the systematic approximations, established by theorem 5, are accurate. As evidence for their accuracy, this section proves that the constructed approximate spatially discrete models are consistent to the original pde eq. 1 to any specified order of error.

As a ‘stepping stone’ to the second-order diffusion pde, and because of its interest in its own right, section 3.2 proves consistency of cognate discrete models of a first-order wave pde. section 3.3 then invokes similar arguments to prove consistency of discrete models of diffusion-like pdes. In this section we restrict attention to homogeneous and autonomous pdes, that is, $f$ and $g$ are here assumed independent of $t,x$ .

section 4 uses an embedding to extend the consistency proof to the homogenisation of systems with microscale heterogeneity.

3.1 Example: consistent discretion of the first-order wave PDE

As a first step, corresponding to example 2 for diffusion, here we analogously construct discrete models of the first-order wave pde

u_{t}+cu_{x}=0\,.

(20)

The spatially discrete models are obtained using the elements and coupling introduced in section 2, and are constructed by computer algebra code (appendix B) that applies deep recursive refinement to algebraically learn the spatial discretisations as a slow manifold in a power series in inter-element coupling $\gamma$ . The code also verifies the high-order consistency between the discrete model and the wave pde (20). In addition, this example explores the upwind/downwind character induced by the parameter $\theta$ in the inter-element coupling edge conditions eq. 6.

A major difference between the diffusion operator $\mathcal{L}=\partial_{x}[\kappa_{j}(x)\partial_{x}]$ in eq. 11 and the advection operator $-c\partial_{x}$ in eq. 20 is that (albeit depending upon boundary conditions) the diffusion operator is self-adjoint and the advection operator is not (it is anti-symmetric). Thus the issue of self-adjointness is not discussed here.

For the wave pde (20) with inter-element edge conditions eq. 6a (flux coupling (6b) does not apply here), the subspace $\mathbb{E}$ of piecewise constant equilibria still exists. A proof is a simpler version of that given for section 2.2: for every piecewise constant field $u^{*}(x)=U_{j}$ in $\mathbb{X}_{j}$ the wave pde (20) is simply $u_{t}=0$ and thus ${\text{\boldmath$U$}}=(U_{1},U_{2},\ldots,U_{N})$ are equilibria of the pde on $\tilde{}\mathbb{X}$ with conditions eq. 6a satisfied for $\gamma=0$ .

For reasons corresponding to those discussed by section 2.3, the spatial discretisation of the wave pde (20) arises as a slow subcentre manifold on coupled elements. Since $u_{j}(t,x)$ denotes the field $u(t,x)$ on the $j$ th element, $u_{j}$ is to satisfy the wave pde (20), namely $u_{jt}=-cu_{jx}$ on $\mathbb{X}_{j}$ . Let the elements be coupled with edge conditions eq. 6a with arbitrary coupling parameter $\gamma$ . Recall that the measure of $u_{j}(t,x)$ in each element is almost immaterial [Roberts2014a, Lemma 5.1], so here we choose each macroscale parameter $U_{j}(t)$ to be the element average.

For the case of equi-sized elements, $H_{j}=H$ , the code of appendix B constructs the slow manifold as a power series in inter-element coupling parameter $\gamma$ . The tuning parameter $\theta$ provides a range of alternative discrete models for the macroscale dynamics. The code finds the wave pde (20) with inter-element coupling eq. 6a has a slow manifold with evolution governed by the odes


	$\displaystyle\dot{U}_{j}$	$\displaystyle=-\frac{\gamma c}{2H}(U_{j+1}-U_{j-1})+\theta\frac{\gamma c}{2H}(U_{j+1}-2U_{j}+U_{j-1})+\mathcal{O}\mathchoice{\big{(}\gamma^{2}\big{)}}{\big{(}\gamma^{2}\big{)}}{(\gamma^{2})}{(\gamma^{2})},$	(21a)
in terms of the element average values $U_{j}$ . At full coupling $\gamma=1$ , the parameter $\theta$ includes the following alternatives: $\theta=0$ gives the anti-symmetric, centred difference, form $\dot{U}_{j}\approx-c(U_{j+1}-U_{j-1})/(2H)$ ; $\theta=1$ gives the upwind (when $c>0$ ) difference $\dot{U}_{j}\approx-c(U_{j}-U_{j-1})/H$ ; whereas $\theta=-1$ gives the downwind (when $c>0$ ) difference $\dot{U}_{j}\approx-c(U_{j+1}-U_{j})/H$ (and vice versa when $c<0$ ).

Such discrete models of the macroscale evolution are learnt in conjunction with the subgrid field structure of the slow manifold. The learnt subgrid field corresponding to the low-order (21a) is the linear $u_{j}=U_{j}+\gamma(\xi-\tfrac{1}{2})(\mu_{j}\delta_{j}-\tfrac{1}{2}\theta_{j}\delta_{j}^{2})U_{j}+\mathcal{O}\mathchoice{\big{(}\gamma^{2}\big{)}}{\big{(}\gamma^{2}\big{)}}{(\gamma^{2})}{(\gamma^{2})}$ , written in terms of the local subgrid variable $\xi:=(x-L_{j})/H$ , and centred mean $\mu_{j}$ and difference $\delta_{j}$ operators (table 1). Higher-order models have more detailed subgrid fields representing more effects of both inter-element communication and subgrid scale physics.

The computer algebra of appendix B easily computes to higher-order in $\gamma$ . Deep iterative refinement learns that the macroscale evolution (21a) on the slow manifold is refined to the following, where $U_{j}$ is omitted for brevity:

$\displaystyle\tfrac{H}{c}\partial_{t}$	$\displaystyle=-\gamma\mu_{j}\delta_{j}+\tfrac{1}{2}\theta\gamma(1-\gamma)\delta_{j}^{2}+\tfrac{1}{4}\gamma^{2}\big{[}1-\tfrac{1}{3}\gamma+\theta^{2}(1-\gamma)\big{]}\mu_{j}\delta_{j}^{3}$
	$\displaystyle\quad{}-\tfrac{1}{8}\theta\gamma^{2}(1-\gamma)(2-\gamma-\theta^{2}\gamma)\delta_{j}^{4}-\tfrac{1}{16}\gamma^{3}\big{[}\tfrac{4}{3}-\gamma+\theta^{2}(4-6\gamma)-\theta^{4}\gamma\big{]}\mu_{j}\delta_{j}^{5}$
	$\displaystyle\quad{}+\mathcal{O}\mathchoice{\big{(}\gamma^{5},\delta_{j}^{6}\big{)}}{\big{(}\gamma^{5},\delta_{j}^{6}\big{)}}{(\gamma^{5},\delta_{j}^{6})}{(\gamma^{5},\delta_{j}^{6})}.$	(21b)

The $\gamma$ -dependence has the appealing feature that when evaluated at full coupling a derived discrete model such as (21b) is independent of the tuning $\theta$ up to some order in $\delta$ : here, for example, at $\gamma=1$ the evolution (21b) becomes simply $\tfrac{H}{c}\partial_{t}=-\mu_{j}\delta_{j}+\tfrac{1}{6}\mu_{j}\delta_{j}^{3}+\mathcal{O}\mathchoice{\big{(}\delta_{j}^{5}\big{)}}{\big{(}\delta_{j}^{5}\big{)}}{(\delta_{j}^{5})}{(\delta_{j}^{5})}$ . That such a discrete model is consistent with the wave pde (20) is seen by the equivalent pde to the discrete model. For example, constructing (21b) to next order, errors $\mathcal{O}\mathchoice{\big{(}\gamma^{6}\big{)}}{\big{(}\gamma^{6}\big{)}}{(\gamma^{6})}{(\gamma^{6})}$ , and substituting expansions for $\delta_{j}=\delta_{x}=2\sinh(H\partial_{x}/2)$ and $\mu_{j}=\mu_{x}=\cosh(H\partial_{x}/2)$ (table 1; when operating on macroscale grid field $U_{j}$ ), the equivalent pde of the discrete model (21b) is

$\displaystyle\frac{1}{c}\partial_{t}U$	$\displaystyle=-\gamma\partial_{x}U+\tfrac{1}{2}\gamma(1-\gamma)\theta H\partial_{x}^{2}U-\tfrac{1}{12}\gamma(1-\gamma)(2-\gamma-3\theta^{2}\gamma)H^{2}\partial_{x}^{3}U$
	$\displaystyle\quad{}+\tfrac{1}{24}\gamma(1-\gamma)\theta\big{[}1-6\gamma+3(1+\theta^{2})\gamma^{2}\big{]}H^{3}\partial_{x}^{4}U$
	$\displaystyle\quad{}-\tfrac{1}{240}\gamma(1-\gamma)\big{[}2-(13+15\theta^{2})\gamma+12(1+5\theta^{2})\gamma^{2}$
	$\displaystyle\qquad{}-3(1+10\theta^{2}+15\theta^{4})\gamma^{3}\big{]}H^{4}\partial_{x}^{5}U\quad{}+\mathcal{O}\mathchoice{\big{(}H^{5}\partial_{x}^{6}U\big{)}}{\big{(}H^{5}\partial_{x}^{6}U\big{)}}{(H^{5}\partial_{x}^{6}U)}{(H^{5}\partial_{x}^{6}U)}.$	(21c)

At full coupling $\gamma=1$ , due to the factors $(1-\gamma)$ , this equivalent pde reduces to the required wave pde (20) for every value of $\theta$ . The order of error in this consistency appears proportional to the chosen order of coupling in $\gamma$ .

Perturbations

One can perturb the wave pde (20) and the computer algebra still constructs a discrete model that is consistent to high-order to the perturbed pde. The code of appendix B caters for perturbations to the first-order wave pde (20) such as

\mathchoice{\frac{\partial u}{\partial t}}{{\partial u}/{\partial t}}{{\partial u}/{\partial t}}{{\partial u}/{\partial t}}+c\mathchoice{\frac{\partial u}{\partial x}}{{\partial u}/{\partial x}}{{\partial u}/{\partial x}}{{\partial u}/{\partial x}}=\alpha\left[c_{0}u+c_{2}\mathchoice{\frac{\partial^{2}u}{\partial x^{2}}}{{\partial^{2}u}/{\partial x^{2}}}{{\partial^{2}u}/{\partial x^{2}}}{{\partial^{2}u}/{\partial x^{2}}}+c_{3}\mathchoice{\frac{\partial^{3}u}{\partial x^{3}}}{{\partial^{3}u}/{\partial x^{3}}}{{\partial^{3}u}/{\partial x^{3}}}{{\partial^{3}u}/{\partial x^{3}}}+c_{4}\mathchoice{\frac{\partial^{4}u}{\partial x^{4}}}{{\partial^{4}u}/{\partial x^{4}}}{{\partial^{4}u}/{\partial x^{4}}}{{\partial^{4}u}/{\partial x^{4}}}\right].

appendix B also codes provision for a microscale lattice version of the advection, namely $u_{t}=-c\big{[}u(t,x+d)-u(t,x-d)\big{]}/(2dH)$ . The computer analysis in both cases learns discrete macroscale models with high-order consistency to the corresponding given perturbed wave system.

Remark 5.

The consistency results of sections 3.2 and 3.3 use deductions expressed in the operator algebra of spatial operators $\partial_{x}$ , $\delta_{x}$ and $\mu_{x}$ . The application of these operators evaluates fields at specific spatial points, but exactly how depends on the scale at which they are operating. At the microscale $u_{j}$ is governed by a given pde (such as (1) in the general case) and this pde justifies the evaluation of such operators. For example, in the context of the wave pde, $\partial_{t}u_{j}=-c\partial_{x}u_{j}$ , effectively $\partial_{x}=-\tfrac{1}{c}\partial_{t}$ and $\delta_{x}=2\sinh[H\partial_{x}/2]=-2\sinh[H\partial_{t}/(2c)]$ (table 1). Thus such spatial operators can be interpreted as transformed to time derivatives via the pde, evaluated at a point, and then transformed back via the pde to spatial operators. Often, when operating at the microscale, we convert to the subgrid variable $\xi=(x-L_{i})/H$ with $\partial_{\xi}=\partial_{x}/H$ , $\delta_{\xi}=\delta_{x}$ and $\mu_{\xi}=\mu_{x}$ . In contrast, the macroscale grid fields $U_{j}$ are discrete so there is no formal spatial derivative, but the main aim of this article is to construct an approximate macroscale pde from a discretised evolution equation of macroscale grid fields on elements $j$ , and to this end we define the spatial derivative at the macroscale in terms of finite-difference equations of element indices $j$ . For example, to obtain (21c) from (21b) we substitute $\delta_{j}=\delta_{x}=2\sinh(H\partial_{x}/2)$ .

3.2 Prove consistency to the first-order wave PDE

As indicated by its equivalent pde (21c), we here prove that the holistic discretisation eq. 21b is consistent to the first-order wave pde eq. 20, for every order of analysis. sections 3.2 and 3.2 establish two key identities for the subsequent consistency theorems 11 and 12, but first we prove coupling identities for derivatives of the subgrid fields which apply in subsequent lemmas.

Lemma 8.

Let the subgrid fields $u_{j}$ be smooth, satisfy the wave pde eq. 20, and be coupled by eq. 6a. Then the edge conditions eq. 6a hold for all spatial derivatives of $u_{j}$ .

Proof.

Denote the $n$ th derivative $u_{j}^{(n)}:=\partial_{x}^{n}u_{j}$ . Proceed via induction. For every $n\geq 1$ consider, recalling superscripts $R,L$ denote evaluation at the right/left edge of an element,

	$\displaystyle c(1-\tfrac{1}{2}\gamma)(u_{j}^{(n)R}-u_{j}^{(n)L})-\tfrac{c}{2}\gamma(u_{j+1}^{(n)L}-u_{j-1}^{(n)R})$
	$\displaystyle{}-\tfrac{c}{2}\gamma\theta(u_{j}^{(n)L}+u_{j}^{(n)R})+\tfrac{c}{2}\gamma\theta(u_{j+1}^{(n)L}+u_{j-1}^{(n)R})$
	$\displaystyle\quad(\text{due to the identity }u_{j}^{(n)}=\partial_{x}u_{j}^{(n-1)})$
	$\displaystyle=c(1-\tfrac{1}{2}\gamma)(\partial_{x}u_{j}^{(n-1)R}-\partial_{x}u_{j}^{(n-1)L})-\tfrac{c}{2}\gamma(\partial_{x}u_{j+1}^{(n-1)L}-\partial_{x}u_{j-1}^{(n-1)R})$
	$\displaystyle{}-\tfrac{c}{2}\gamma\theta(\partial_{x}u_{j}^{(n-1)L}+\partial_{x}^{2}u_{j}^{(n-1)R})+\tfrac{c}{2}\gamma\theta(\partial_{x}u_{j+1}^{(n-1)L}+\partial_{x}u_{j-1}^{(n-1)R})$
	$\displaystyle\quad(\text{then using the {pde}\ as in \lx@cref{creftype~refnum}{remcom}})$
	$\displaystyle=-(1-\tfrac{1}{2}\gamma)(\partial_{t}u_{j}^{(n-1)R}-\partial_{t}u_{j}^{(n-1)L})+\tfrac{1}{2}\gamma(\partial_{t}u_{j+1}^{(n-1)L}-\partial_{t}u_{j-1}^{(n-1)R})$
	$\displaystyle{}+\tfrac{1}{2}\gamma\theta(\partial_{t}u_{j}^{(n-1)L}+\partial_{t}u_{j}^{(n-1)R})-\tfrac{1}{2}\gamma\theta(\partial_{t}u_{j+1}^{(n-1)L}+\partial_{t}u_{j-1}^{(n-1)R})$
	$\displaystyle=-\partial_{t}\left\{(1-\tfrac{1}{2}\gamma)(u_{j}^{(n-1)R}-u_{j}^{(n-1)L})-\tfrac{1}{2}\gamma(u_{j+1}^{(n-1)L}-u_{j-1}^{(n-1)R})\right.$
	$\displaystyle{}\left.-\tfrac{1}{2}\gamma\theta(u_{j}^{(n-1)L}+u_{j}^{(n-1)R})+\tfrac{1}{2}\gamma\theta(u_{j+1}^{(n-1)L}+u_{j-1}^{(n-1)R})\right\}$
	$\displaystyle=\partial_{t}\left\{0\right\}=0\,,$

where the penultimate equality above holds provided $u^{(n-1)}_{j}$ satisfies eq. 6a. By induction, since $u^{(0)}$ satisfies eq. 6a then so does $u^{(n)}$ for every $n\geq 0$ . ∎

Now we recast the coupling condition so the spatial difference is on the left-hand side and all coupling terms are on the right-hand side.

Lemma 9.

Let all elements be of length $H_{j}=H$ . Define the mid-element locations $X_{j}:=(R_{j}+L_{j})/2$ . Then the inter-element edge condition eq. 6a is equivalent to the coupling condition

\delta_{x}u_{j}=\frac{\gamma(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\gamma\theta\mu_{j}\delta_{j}+\tfrac{1}{4}\gamma^{2}(\theta^{2}-1)\delta_{j}^{2}}}u_{j}\quad\text{at }x=X_{j}\,.

(22)

Further, this identity also holds for every spatial derivative $\partial_{x}^{n}u_{j}$ evaluated at the mid-element point $x=X_{j}$ .

Proof.

We first rewrite the edge conditions eq. 6a in terms of mean and difference operators $\mu$ and $\delta$ . Identities herein involving $u_{j}$ are evaluated “at $x=X_{j}$ ”—equivalently “at $\xi=1/2$ ” in terms of the subgrid space variable $\xi:=(x-L_{j})/H$ (as introduced in the subgrid fields (9)) because in this proof we regard $u_{j}$ as a function of $\xi$ , not directly $x$ . Using the spatial shift operators $E_{j}$ and $E_{\xi}$ (table 1), $u_{j}^{R}=E_{\xi}^{1/2}u_{j}$ , $u_{j}^{L}=E_{\xi}^{-1/2}u_{j}$ , $u_{j-1}^{R}=E_{j}^{-1}E_{\xi}^{1/2}u_{j}$ , $u_{j+1}^{L}=E_{j}E_{\xi}^{-1/2}u_{j}$ , which we substitute into (6a), and then replace the shift operators with mean and difference operators to obtain after some rearrangement that (6a) is the same as

(1+\tfrac{1}{4}\gamma\delta_{j}^{2}-\tfrac{1}{2}\gamma\theta\mu_{j}\delta_{j})\delta_{\xi}u_{j}=\gamma(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}\mu_{\xi}u_{j}\,.

(23)

On squaring the operators⁵⁵5 Equation eq. 23 is an operation on $u_{j}$ evaluated at $\xi=\frac{1}{2}$ to define the coupling conditions at element edges $\xi=0,1$ . When squaring the operators we must smoothly extend $u_{j}$ over $\xi\in[-\frac{1}{2},0),(1,\frac{3}{2}]$ (as defined by the pde) with the new ‘squared operator’ equation defining new coupling conditions at $\xi=-\frac{1}{2},\frac{3}{2}$ (or equivalently, at $x=X_{j-1},X_{j+1}$ but still in the smooth extension of element $j$ ). For example, $\delta^{2}_{\xi}u_{j}(\xi=\frac{1}{2})=u_{j}(\frac{3}{2})-u_{j}(-\frac{1}{2})$ and we remain in a smooth extension of element $j$ .

\displaystyle(1+\tfrac{1}{4}\gamma\delta_{j}^{2}-\tfrac{1}{2}\gamma\theta\mu_{j}\delta_{j})^{2}\delta_{\xi}^{2}u_{j}=\gamma^{2}(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})^{2}\delta_{j}^{2}\mu_{\xi}^{2}u_{j}\,.

Expand both sides, and on the left substitute $\mu_{j}^{2}=1+\frac{1}{4}\delta_{j}^{2}$ while on the right substitute $\mu_{\xi}^{2}=1+\frac{1}{4}\delta_{\xi}^{2}$ (table 1), to obtain

\big{[}1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\gamma\theta\mu_{j}\delta_{j}+\tfrac{1}{4}\gamma^{2}(\theta^{2}-1)\delta_{j}^{2}\big{]}\delta_{\xi}^{2}u_{j}={\gamma^{2}(\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2})^{2}\delta_{j}^{2}}u_{j}\,.

Since $\delta_{\xi}=\delta_{x}$ on the microscale, dividing the above by the $[\cdot]$ factor, taking the square root, and recalling that this operator applies to $u_{j}$ at $\xi=\tfrac{1}{2}$ (that is, $x=X_{j}$ ), we obtain (22). ∎

Lemma 10.

The first-order wave pde eq. 20 on elements with coupling eq. 22 has macroscale dynamics such that

\partial_{t}U_{j}=-\frac{2c}{H}\operatorname{asinh}\left[\frac{\tfrac{1}{2}\gamma(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\gamma\theta\mu_{j}\delta_{j}+\tfrac{1}{4}\gamma^{2}(\theta^{2}-1)\delta_{j}^{2}}}\right]U_{j}\,,

(24)

where macroscale grid field $U_{j}$ is a local average of $u_{j}(x)$ about the centre of the $j$ th element; that is, $U_{j}:=\tfrac{1}{2\ell}\int_{X_{j}-\ell}^{X_{j}+\ell}w(x-X_{j})u_{j}(x)\,dx$ for some $\ell\leq H/2$ and some weight function $w$ (the macroscale field reduces to the mid-element value $U_{j}:=u_{j}(X_{j})$ as $\ell\rightarrow 0$ and for weight $w(0)=1$ ).

Knowing how spatial differences of pde solutions transform now translates into the corresponding evolution of the macroscale variables.

Proof.

Consider the wave pde eq. 20 in the $j$ th element and higher order spatial derivatives (and using table 1): for every $p=0,1,2,\ldots$ , $\partial_{t}\partial_{x}^{p}u_{j}=-c\partial_{x}\partial_{x}^{p}u_{j}=-\tfrac{c}{H}2\operatorname{asinh}(\delta_{x}/2)\partial_{x}^{p}u_{j}$ . Inverting the operator $\operatorname{asinh}(\delta_{x}/2)$ , this identity is equivalent to $2\sinh(-H\partial_{t}/2c)\partial_{x}^{p}u_{j}=\delta_{x}\partial_{x}^{p}u_{j}$ . Evaluating this at the element mid-point and by the coupling eq. 22,

	$\displaystyle 2\sinh(-H\partial_{t}/2c)\partial_{x}^{p}u_{j}\big{\|}_{X_{j}}$	$\displaystyle=\delta_{x}\partial_{x}^{p}u_{j}\big{\|}_{X_{j}}$
		$\displaystyle=\frac{\gamma(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\gamma\theta\mu_{j}\delta_{j}+\tfrac{1}{4}\gamma^{2}(\theta^{2}-1)\delta_{j}^{2}}}\partial_{x}^{p}u_{j}\big{\|}_{X_{j}}\,.$

Premultiplying by $(x-X_{j})^{p}/p!$ provides the $p$ th term in a Taylor expansion of $u_{j}(x)$ about $x=X_{j}$ (the term $(x-X_{j})^{p}=H^{p}(\xi-\frac{1}{2})^{p}$ commutes with $\partial_{t}$ , $\delta_{j}$ and $\mu_{j}$ since it is independent of $t$ and $j$ , but it does not commute with $\delta_{x}=\delta_{\xi}$ when $p>0$ ) and on summing over all $p=0,1,2,\ldots$ we obtain

2\sinh(-H\partial_{t}/2c)u_{j}(x)=\frac{\gamma(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\gamma\theta\mu_{j}\delta_{j}+\tfrac{1}{4}\gamma^{2}(\theta^{2}-1)\delta_{j}^{2}}}u_{j}(x)

(25)

in some open interval containing $X_{j}$ . Then multiply by the weight function $w(x-X_{j})$ and integrate over the centre of the $j$ th element from $x=X_{j}-\ell$ to $X_{j}+\ell$ and finally revert the $\sinh(-H\partial_{t}/2c)$ operator to deduce eq. 24; or, for the averaging width $\ell\rightarrow 0$ , set $x=X_{j}$ and revert the $\sinh$ operator to deduce eq. 24. ∎

The above proof also shows that eq. 22 holds for $u_{j}$ evaluated at all $x\in(X_{j}-\ell,X_{j}+\ell)$ , not just $x=X_{j}$ . To prove this, in eq. 25 substitute $2\sinh(-H\partial_{t}/2c)=2\sinh(H\partial_{x}/2)=\delta_{x}$ .

Computer algebra readily learns the dynamics (24) in the form of a Taylor series of the coupling parameter $\gamma$ . That is, the computer algebra recovers precisely the series eq. 21b to a chosen order of $\gamma$ . This recovery of the series verifies sections 3.2, 3.2 and 3.2.

Theorem 11.

At full coupling $\gamma=1$ , and for every $\theta$ , the holistic discretisation (24), equivalently (21b), is consistent with the first-order wave pde (20).

Proof.

Evaluating eq. 24 at full-coupling $\gamma=1$ gives (table 1)

	$\displaystyle\partial_{t}U_{j}$	$\displaystyle=-\frac{2c}{H}\operatorname{asinh}\left[\frac{\tfrac{1}{2}(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{1+\tfrac{1}{2}\delta_{j}^{2}-\theta\mu_{j}\delta_{j}+\tfrac{1}{4}(\theta^{2}-1)\delta_{j}^{2}}}\right]U_{j}$
		$\displaystyle=-\frac{2c}{H}\operatorname{asinh}\left[\frac{\tfrac{1}{2}(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{\mu_{j}^{2}-\theta\mu_{j}\delta_{j}+\tfrac{1}{4}\theta^{2}\delta_{j}^{2}}}\right]U_{j}$
		$\displaystyle=-\frac{2c}{H}\operatorname{asinh}\left[\frac{\tfrac{1}{2}(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}}{\sqrt{(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})^{2}}}\right]U_{j}$
		$\displaystyle=-\frac{2c}{H}\operatorname{asinh}\left[\delta_{j}/2\right]U_{j}=-\frac{2c}{H}\operatorname{asinh}\left[\delta_{x}/2\right]U_{j}\quad\text{(\lx@cref{creftype~refnum}{remcom})}$
		$\displaystyle=-\frac{c}{H}\partial_{j}U_{j}=-c\partial_{x}U_{j}\,.$

That $\partial_{t}U_{j}=-c\partial_{x}U_{j}$ confirms that computing on elements with full coupling recovers macroscale simulations and predictions consistent with the macroscale dynamics of the first-order wave pde eq. 20. ∎

Substituting $\gamma=1$ reduces the series eq. 21b, obtained from computer algebra, to the first-order wave pde eq. 20 and verifies theorem 11. However, a practical derivation of the series eq. 21b necessarily truncates to some finite order in the coupling, say to errors $o\big{(}\gamma^{p}\big{)}$ . Two important questions concerning this error in $\gamma$ and the practical application of this holistic discretisation are: how is this error reinterpreted in terms of a stencil of nearest-element coupling? and how does it affect the consistency of the series eq. 21b with the wave pde eq. 20?

Theorem 12.

Consider constructing, to errors $o\big{(}\gamma^{p}\big{)}$ , the slow manifold of the first-order wave pde eq. 20 on elements coupled by eq. 6a. The resulting holistic discretisation (that is, eq. 21b truncated to errors $o\big{(}\gamma^{p}\big{)}$ ) is

•

a scheme of stencil width $(2p-1)$ on the macroscale grid, and
•

consistent with the wave pde to errors $o\big{(}\partial_{x}^{p}\big{)}$ .

Proof.

Consider macroscale solutions for which the differences $\delta_{j}$ are ‘small’ (i.e., $U_{j}$ varies slowly with changes in $j$ ). For such solutions, the Taylor series in $\delta_{j}$ of the holistic discretisation eq. 24 is appropriate. But we only need the pattern of powers of $\gamma$ , $\mu_{j}$ and $\delta_{j}$ , so in the derivation we omit almost all details of the coefficients, and use $\sim$ instead of $=$ to denote such omission.

The $\operatorname{asinh}$ term in eq. 24 is

\operatorname{asinh}\left[\frac{\tfrac{1}{2}(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})(\gamma\delta_{j})}{\sqrt{1+\tfrac{1}{2}\delta_{j}(\gamma\delta_{j})-\theta\mu_{j}(\gamma\delta_{j})+\tfrac{1}{4}(\theta^{2}-1)(\gamma\delta_{j})^{2}}}\right],

(26)

where we pair each instance of $\gamma$ with an instance of $\delta_{j}$ . As all $\gamma$ are paired with a $\delta_{j}$ , but not all $\delta_{j}$ (or $\mu_{j}$ ) are paired with a $\gamma$ , we conclude that in a Taylor expansion in small $\delta_{j}$ , up to a given power $\delta_{j}^{p}$ , every $\gamma$ exponent is $\leq p$ (as in eq. 21b). Furthermore, truncating construction to errors $o\big{(}\gamma^{p}\big{)}$ ensures that the terms are all complete and correct for every power of $\delta_{j}$ less than $p$ ; that is, the error is $o\big{(}\delta_{j}^{p}\big{)}$ . Since $\partial_{x}\sim\delta_{j}$ when operating on macroscale grid fields, in such a construction the consistency error is $o\big{(}\partial_{x}^{p}\big{)}$ as asserted in the theorem.

To determine the width of the stencil it is sufficient to expand

	$\displaystyle\operatorname{asinh}[\cdot]$	$\displaystyle\sim\operatorname{asinh}\left[\frac{(\mu_{j}+\delta_{j})(\gamma\delta_{j})}{\sqrt{1+(\mu_{j}+\delta_{j})(\gamma\delta_{j})+(\gamma\delta_{j})^{2}}}\right]$
		$\displaystyle\sim\sum_{m=0}^{\infty}(\mu_{j}+\delta_{j})^{2m+1}(\gamma\delta_{j})^{2m+1}[1+(\mu_{j}+\delta_{j})(\gamma\delta_{j})+(\gamma\delta_{j})^{2}]^{-m-1/2}$
		$\displaystyle\sim\sum_{m=0}^{\infty}(\mu_{j}+\delta_{j})^{2m+1}(\gamma\delta_{j})^{2m+1}\sum_{n=0}^{\infty}(\mu_{j}+\delta_{j})^{0:n}(\gamma\delta_{j})^{n:2n}$
		$\displaystyle\sim\sum_{m,n=0}^{\infty}(\mu_{j}+\delta_{j})^{2m+1:2m+n+1}(\gamma\delta_{j})^{2m+n+1:2m+2n+1}\,.$

On expanding to errors $o\big{(}\gamma^{p}\big{)}$ we retain up to $\gamma^{p-1}$ . Then, the lowest order $\gamma$ which obtains the highest order $E_{j}^{\pm 1/2}\sim\mu_{j},\delta_{j}$ is when $p-1=2m+n+1$ , so that in the above expansion we have up to $E_{j}^{[(2m+n+1)+(2m+n+1)]/2}\sim E_{j}^{\pm(p-1)}$ . So the stencil is of width $(2p-1)$ . ∎

This completes the proof of the consistency of these fixed-bandwidth holistic discretisations, elements coupled by (6a), to the first-order wave pde eq. 20.

3.3 Holistic discretisation of generalised diffusion is consistent

Recall that example 2 shows that in the case of the diffusion pde, $u_{t}=u_{xx}$ , the holistic discretisation eq. 9c is consistent with the pde to all computed orders of accuracy. The computer algebra (appendix A) deriving the results eq. 9 is easily modified to analyse more general pdes (such as advection-diffusion, $u_{t}=-cu_{x}+u_{xx}$ , or lattice diffusion, $u_{t}=u|_{x+d}-2u+u|_{x-d}$ ) and in all cases investigated we found corresponding high-order consistency. Consequently, here we prove consistency for a generalised class of pdes. The proofs in this section are analogous to those in section 3.2 for wave pdes.

The proof of consistency is cognate to an earlier proof [Roberts2011a, §6.1] of high-order consistency emerging from a distinctly different inter-element coupling. However, the new inter-element coupling (6), implemented on element edges in order to preserve important symmetries, requires new proof.

The analysis here explores, for some Hilbert space $\mathbb{U}$ , the evolution in time $t$ of a $\mathbb{U}$ -valued field $u(t,x)$ over space $x\in\mathbb{X}$ . We restrict attention to smooth $\mathbb{U}$ -valued fields $u$ in the vector space $\mathbb{C}:=C^{\omega}(\mathbb{X})$ , and correspondingly $\mathbb{C}_{j}:=C^{\omega}(\mathbb{X}_{j})$ on the spatial elements. The field $u$ is to satisfy the isotropic, homogenous, generalised diffusion-like pde, the evolution equation,

u_{t}=\mathcal{K}_{0}u+\mathcal{K}(\partial_{x}^{2})u\,,

(27)

for some linear operators $\mathcal{K}_{k}:\mathbb{U}\to\mathbb{U}$ and isotropic invertible $\mathcal{K}:\mathbb{C}\to\mathbb{C}$ , with inverse $\mathcal{K}^{-1}$ , where $\mathcal{K}$ is formally expandable as $\mathcal{K}=\sum_{k=1}^{\infty}\mathcal{K}_{k}\partial_{x}^{2k}$ with $\mathcal{K}_{1}$ invertible.

The most basic example of (27) is simple diffusion, $u_{t}=u_{xx}$ , obtained with $\mathbb{U}=\mathbb{R}$ , zero $\mathcal{K}_{0}$ , and identity $\mathcal{K}$ . For a second example, with microscale space step $d$ the difference equation $u_{t}=\delta^{2}u:=u|_{x+d}-2u+u|_{x-d}$ corresponds to (table 1) setting $\mathcal{K}_{0}:=0$ and $\mathcal{K}(\partial_{x}^{2}):=4\sinh^{2}(d\partial_{x}/2)=d^{2}\partial_{x}^{2}+\tfrac{1}{12}d^{4}\partial_{x}^{4}+\tfrac{1}{360}d^{6}\partial_{x}^{6}+\cdots$ and so $\mathcal{K}_{1}=d^{2}$ . Further in this second example, the invertibility gives that since $\delta^{2}=\mathcal{K}(\partial_{x}^{2})$ then $\partial_{x}^{2}=\mathcal{K}^{-1}(\delta^{2})=(4/d^{2})\operatorname{asinh}^{2}(\delta/2)=\tfrac{1}{d^{2}}\big{[}\delta^{2}-\tfrac{1}{12}\delta^{4}+\tfrac{1}{90}\delta^{6}-\cdots\big{]}$ .

sections 3.3, 3.3, 3.3 and 3.3 establish key identities for the subsequent consistency theorems 17 and 18.

Lemma 13.

Define the transformed time-derivative operator $\mathcal{T}:=\mathcal{K}^{-1}(\partial_{t}-\mathcal{K}_{0})$ . Then the evolution equation eq. 27 is equivalent to

\mathcal{T}u=u_{xx}\,.

(28)

Proof.

As in remark 5, and since $\mathcal{K}$ is invertible, we rearrange the operator equation $\partial_{t}=\mathcal{K}_{0}+\mathcal{K}(\partial_{x}^{2})$ to $\mathcal{K}^{-1}(\partial_{t}-\mathcal{K}_{0})=\partial_{x}^{2}$ . Hence, (28) is equivalent to the evolution equation (27). ∎

Lemma 14.

Let the subgrid fields $u_{j}\in\mathbb{C}_{j}$ satisfy the evolution equation eq. 27, equivalently (28), and be coupled by edge conditions eq. 6 with flux $f_{j}:=-u_{jx}$ . Then the edge condition eq. 6a holds for all even spatial derivatives of $u_{j}$ , and the edge condition eq. 6b holds for all odd spatial derivatives of $u_{j}$ .

Proof.

Denote the $n$ th derivative $u_{j}^{(n)}:=\partial_{x}^{n}u_{j}$ . Proceed via induction. From the definition of the flux, multiplying the flux edge condition eq. 6b by $-1$ is the edge condition for the first derivative of $u_{j}$ , and this coupling is equivalent to eq. 6a with $\theta\mapsto-\theta$ . Then for every $n\geq 2$ consider

	$\displaystyle(1-\tfrac{1}{2}\gamma)(u_{j}^{(n)R}-u_{j}^{(n)L})-\tfrac{1}{2}\gamma(u_{j+1}^{(n)L}-u_{j-1}^{(n)R})$
	$\displaystyle\quad{}-\tfrac{1}{2}\gamma(-1)^{n}\theta(u_{j}^{(n)L}+u_{j}^{(n)R})+\tfrac{1}{2}\gamma(-1)^{n}\theta(u_{j+1}^{(n)L}+u_{j-1}^{(n)R})$
	$\displaystyle\quad(\text{due to the identity }u_{j}^{(n)}=\partial_{x}^{2}u_{j}^{(n-2)})$
	$\displaystyle=(1-\tfrac{1}{2}\gamma)(\partial_{x}^{2}u_{j}^{(n-2)R}-\partial_{x}^{2}u_{j}^{(n-2)L})-\tfrac{1}{2}\gamma(\partial_{x}^{2}u_{j+1}^{(n-2)L}-\partial_{x}^{2}u_{j-1}^{(n-2)R})$
	$\displaystyle\quad{}-\tfrac{1}{2}\gamma(-1)^{n}\theta(\partial_{x}^{2}u_{j}^{(n-2)L}+\partial_{x}^{2}u_{j}^{(n-2)R})+\tfrac{1}{2}\gamma(-1)^{n}\theta(\partial_{x}^{2}u_{j+1}^{(n-2)L}+\partial_{x}^{2}u_{j-1}^{(n-2)R})$
	$\displaystyle\quad(\text{then using \eqref{Xeqrpde}})$
	$\displaystyle=(1-\tfrac{1}{2}\gamma)(\mathcal{T}u_{j}^{(n-2)R}-\mathcal{T}u_{j}^{(n-2)L})-\tfrac{1}{2}\gamma(\mathcal{T}u_{j+1}^{(n-2)L}-\mathcal{T}u_{j-1}^{(n-2)R})$
	$\displaystyle\quad{}-\tfrac{1}{2}\gamma(-1)^{n}\theta(\mathcal{T}u_{j}^{(n-2)L}+\mathcal{T}u_{j}^{(n-2)R})+\tfrac{1}{2}\gamma(-1)^{n}\theta(\mathcal{T}u_{j+1}^{(n-2)L}+\mathcal{T}u_{j-1}^{(n-2)R})$
	$\displaystyle=\mathcal{T}\left\{(1-\tfrac{1}{2}\gamma)(u_{j}^{(n-2)R}-u_{j}^{(n-2)L})-\tfrac{1}{2}\gamma(u_{j+1}^{(n-2)L}-u_{j-1}^{(n-2)R})\right.$
	$\displaystyle\left.\quad{}-\tfrac{1}{2}\gamma(-1)^{n}\theta(u_{j}^{(n-2)L}+u_{j}^{(n-2)R})+\tfrac{1}{2}\gamma(-1)^{n}\theta(u_{j+1}^{(n-2)L}+u_{j-1}^{(n-2)R})\right\}$
	$\displaystyle=\mathcal{T}\left\{0\right\}=0\,,$

where the penultimate equality above holds provided $u^{(n-2)}_{j}$ satisfies eq. 6a when $n$ is even and eq. 6b when $n$ is odd. By induction, since $u^{(0)}$ satisfies eq. 6a then so does $u^{(n)}$ for every even $n\geq 0$ , and since $u^{(1)}$ satisfies eq. 6b then so does $u^{(n)}$ for every odd $n\geq 0$ . ∎

In the case of simple diffusion, the code of appendix A verifies using computer algebra that for the constructed slow manifold discretisation (section 2.2), the field $u_{j}(t,x)$ satisfies the following identity eq. 29. We now prove this identity more generally. The proof is similar to that of section 3.2, but a little more complex as we now have edge conditions for both field and flux.

Lemma 15.

Recall the mid-element locations $X_{j}:=(R_{j}+L_{j})/2$ . For $u_{j}(t,x)$ as specified in section 3.3, and $u_{j}$ smoothly extended to being over $[X_{j-1},X_{j+1}]$ ,

\delta_{x}^{2}u_{j}=\frac{\gamma^{2}(\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2})\delta_{j}^{2}}{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\tfrac{1}{4}\gamma^{2}(\theta^{2}+1)\delta_{j}^{2}}u_{j}\quad\text{at }x=X_{j}\,.

(29)

Further, this identity also holds for every spatial derivative $\partial_{x}^{n}u_{j}$ .

Proof.

The field and flux edge conditions eq. 6, and their derivatives (section 3.3), imply that, in terms of mean and difference operators, the following holds at $x=X_{j}$ , for every integer $k=0,1,\ldots$ ,

	$\displaystyle[\delta_{x}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{x}-\mu_{j}\mu_{x})\delta_{j}]\partial_{x}^{2k}u_{j}\quad$	$\displaystyle=\phantom{+}\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{x}-\delta_{j}\mu_{x})\delta_{j}\partial_{x}^{2k}u_{j}\,,$
	$\displaystyle[\delta_{x}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{x}-\mu_{j}\mu_{x})\delta_{j}]\partial_{x}^{2k+1}u_{j}$	$\displaystyle=-\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{x}-\delta_{j}\mu_{x})\delta_{j}\partial_{x}^{2k+1}u_{j}\,.$		(30)

Upon multiplying by appropriate factors, and changing to subgrid variable $\xi=(x-L_{j})/H$ , the above set of edge conditions is equivalent to the following set: at $\xi=\tfrac{1}{2}$ , for every $k=0,1,\ldots$


$\displaystyle\frac{[\delta_{\xi}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{\xi}-\mu_{j}\mu_{\xi})\delta_{j}]\ \partial_{\xi}^{2k}u_{j}}{2^{2k}(2k)!}\quad$	$\displaystyle=\phantom{+}\frac{\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{\xi}-\delta_{j}\mu_{\xi})\delta_{j}\ \partial_{\xi}^{2k}u_{j}}{2^{2k}(2k)!}\,,$	(31a)
$\displaystyle\frac{[\delta_{\xi}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{\xi}-\mu_{j}\mu_{\xi})\delta_{j}]\ \partial_{\xi}^{2k+1}u_{j}}{2^{2k+1}(2k+1)!}$	$\displaystyle=-\frac{\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{\xi}-\delta_{j}\mu_{\xi})\delta_{j}\ \partial_{\xi}^{2k+1}u_{j}}{2^{2k+1}(2k+1)!}\,.$	(31b)

Summing each of (31) over $k=0,1,2,\ldots$ , and for simplicity here omitting “ $u_{j}$ at $\xi=\tfrac{1}{2}$ ”, gives

	$\displaystyle[\delta_{\xi}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{\xi}-\mu_{j}\mu_{\xi})\delta_{j}]\cosh(\tfrac{1}{2}\partial_{\xi})$	$\displaystyle=\phantom{+}\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{\xi}-\delta_{j}\mu_{\xi})\delta_{j}\cosh(\tfrac{1}{2}\partial_{\xi}),$
	$\displaystyle[\delta_{\xi}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{\xi}-\mu_{j}\mu_{\xi})\delta_{j}]\sinh(\tfrac{1}{2}\partial_{\xi})$	$\displaystyle=-\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{\xi}-\delta_{j}\mu_{\xi})\delta_{j}\sinh(\tfrac{1}{2}\partial_{\xi}).$

Recall that $\cosh(\tfrac{1}{2}\partial_{\xi})=\mu_{\xi}$ and $\sinh(\tfrac{1}{2}\partial_{\xi})=\tfrac{1}{2}\delta_{\xi}$ (table 1). Consequently the above pair of operator equations are

	$\displaystyle[\delta_{\xi}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{\xi}-\mu_{j}\mu_{\xi})\delta_{j}]\mu_{\xi}$	$\displaystyle=\phantom{+}\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{\xi}-\delta_{j}\mu_{\xi})\delta_{j}\mu_{\xi}\,,$
	$\displaystyle[\delta_{\xi}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{\xi}-\mu_{j}\mu_{\xi})\delta_{j}]\delta_{\xi}$	$\displaystyle=-\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{\xi}-\delta_{j}\mu_{\xi})\delta_{j}\delta_{\xi}\,.$

Each of these equations rearrange, respectively, to

	$\displaystyle(1+\tfrac{1}{4}\gamma\delta_{j}^{2}-\tfrac{1}{2}\gamma\theta\mu_{j}\delta_{j})\mu_{\xi}\delta_{\xi}$	$\displaystyle=\gamma(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}\mu_{\xi}^{2}\,,$
	$\displaystyle(1+\tfrac{1}{4}\gamma\delta_{j}^{2}+\tfrac{1}{2}\gamma\theta\mu_{j}\delta_{j})\delta_{\xi}^{2}\quad$	$\displaystyle=\gamma(\mu_{j}+\tfrac{1}{2}\theta\delta_{j})\delta_{j}\mu_{\xi}\delta_{\xi}\,.$

Multiply both sides of the first equation by $\gamma(\mu_{j}+\tfrac{1}{2}\theta\delta_{j})\delta_{j}$ , then substitute the second to obtain⁶⁶6Here we smoothly extend the $j$ th element ( $\xi\in[0,1]$ ) to also include $\xi\in[-\frac{1}{2},0),(1,\frac{3}{2}]$ and the following equation is an additional coupling condition for element $j$ at $\xi=-\frac{1}{2},\frac{3}{2}$ , or equivalently $x=X_{j-1},X_{j+1}$ (as in the proof of section 3.2).

	$\displaystyle(1+\tfrac{1}{4}\gamma\delta_{j}^{2}-\tfrac{1}{2}\gamma\theta\mu_{j}\delta_{j})(1+\tfrac{1}{4}\gamma\delta_{j}^{2}+\tfrac{1}{2}\gamma\theta\mu_{j}\delta_{j})\delta_{\xi}^{2}$
	$\displaystyle=\gamma^{2}(\mu_{j}+\tfrac{1}{2}\theta\delta_{j})(\mu_{j}-\tfrac{1}{2}\theta\delta_{j})\delta_{j}^{2}\,.$

Expand both sides, and on the left substitute $\mu_{j}^{2}=1+\frac{1}{4}\delta_{j}^{2}$ whereas on the right substitute $\mu_{\xi}^{2}=1+\frac{1}{4}\delta_{\xi}^{2}$ (table 1), to obtain

\big{[}1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\tfrac{1}{4}\gamma^{2}(\theta^{2}+1)\delta_{j}^{2}\big{]}\delta_{\xi}^{2}={\gamma^{2}(\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2})\delta_{j}^{2}}\,.

Since $\delta_{\xi}=\delta_{x}$ , dividing the above by the $[\cdot]$ factor, and now explicitly applying to $u_{j}$ at $\xi=\tfrac{1}{2}$ (that is, $x=X_{j}$ ) gives (29).

Furthermore, to prove that (29) holds for all spatial derivatives $\partial_{x}^{n}u_{j}$ , adapt (30) for $n$ even, so that $n+2k$ is even:

	$\displaystyle[\delta_{x}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{x}-\mu_{j}\mu_{x})\delta_{j}]\partial_{x}^{2k}\partial_{x}^{n}u_{j}\quad$	$\displaystyle=\phantom{+}\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{x}-\delta_{j}\mu_{x})\delta_{j}\partial_{x}^{2k}\partial_{x}^{n}u_{j}\,,$
	$\displaystyle[\delta_{x}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{x}-\mu_{j}\mu_{x})\delta_{j}]\partial_{x}^{2k+1}\partial_{x}^{n}u_{j}$	$\displaystyle=-\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{x}-\delta_{j}\mu_{x})\delta_{j}\partial_{x}^{2k+1}\partial_{x}^{n}u_{j}\,,$

and for $n$ odd, so that $n+2k$ is odd:

	$\displaystyle[\delta_{x}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{x}-\mu_{j}\mu_{x})\delta_{j}]\partial_{x}^{2k}\partial_{x}^{n}u_{j}\quad$	$\displaystyle=-\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{x}-\delta_{j}\mu_{x})\delta_{j}\partial_{x}^{2k}\partial_{x}^{n}\,,$
	$\displaystyle[\delta_{x}+\gamma(\tfrac{1}{4}\delta_{j}\delta_{x}-\mu_{j}\mu_{x})\delta_{j}]\partial_{x}^{2k+1}\partial_{x}^{n}u_{j}$	$\displaystyle=\phantom{+}\tfrac{1}{2}\gamma\theta(\mu_{j}\delta_{x}-\delta_{j}\mu_{x})\delta_{j}\partial_{x}^{2k+1}\partial_{x}^{n}u_{j}\,.$

Then, in the proof of (29) for $u_{j}$ we replace $u_{j}$ with $\partial_{x}^{n}u_{j}$ and $\theta$ with $(-1)^{n}\theta$ . Since only $\theta^{2}$ terms arise in (29), it thus also holds for $\partial_{x}^{n}u_{j}$ . ∎

Lemma 16.

The evolution equation (27), equivalently (28), on equi-sized elements for $u_{j}\in\mathbb{C}_{j}$ coupled by eq. 6, has dynamics satisfying both the following two spatially discrete equations:


	$\displaystyle\mathcal{T}U_{j}=\frac{4}{H^{2}}\operatorname{asinh}^{2}\left[\sqrt{\frac{\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2}}{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\tfrac{1}{4}\gamma^{2}(\theta^{2}+1)\delta_{j}^{2}}}\frac{\gamma\delta_{j}}{2}\right]U_{j}\,,$		(32a)
	$\displaystyle\partial_{t}U_{j}=\mathcal{K}_{0}U_{j}+\mathcal{K}\left\{\frac{4}{H^{2}}\operatorname{asinh}^{2}\left[\sqrt{\frac{\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2}}{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\tfrac{1}{4}\gamma^{2}(\theta^{2}+1)\delta_{j}^{2}}}\frac{\gamma\delta_{j}}{2}\right]\right\}U_{j}\,,$		(32b)

where macroscale grid field $U_{j}(t)\in\mathbb{U}$ is a local average of $u_{j}(t,x)$ about the centre of the $j$ th element; that is, $U_{j}(t):=\tfrac{1}{2\ell}\int_{X_{j}-\ell}^{X_{j}+\ell}w(x-X_{j})u_{j}(t,x)\,dx$ for some $\ell\leq H/2$ and some weight function $w$ (the macroscale field reduces to the mid-element value $U_{j}(t):=u_{j}(t,X_{j})$ as $\ell\rightarrow 0$ and for weight $w(0)=1$ ).

Proof.

We adapt the proof of section 3.2. Given (29) of section 3.3, first consider the evolution equation (28) in the elements (and smoothly extended from the elements): for every $p=0,1,2,\ldots$ , $\mathcal{T}\partial_{x}^{p}u_{j}=\partial_{x}^{2}\partial_{x}^{p}u_{j}=\tfrac{4}{H^{2}}\operatorname{asinh}^{2}(\delta_{x}/2)\partial_{x}^{p}u_{j}=\tfrac{4}{H^{2}}\operatorname{asinh}^{2}(\sqrt{\delta_{x}^{2}/4})u_{j}$ . Upon inverting the operator $\operatorname{asinh}^{2}(\sqrt{\cdot/4})$ , this evolution equation is equivalent to $4\sinh^{2}(H\sqrt{\mathcal{T}}/2)\partial_{x}^{p}u_{j}=\delta_{x}^{2}\partial_{x}^{p}u_{j}$ . Evaluating this form of the evolution equation at the element mid-point, and using eq. 29, gives

4\sinh^{2}(H\sqrt{\mathcal{T}}/2)\partial_{x}^{p}u_{j}\big{|}_{X_{j}}=\delta_{x}^{2}\partial_{x}^{p}u_{j}\big{|}_{X_{j}}=\frac{\gamma^{2}(\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2})\delta_{j}^{2}}{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\tfrac{1}{4}\gamma^{2}(\theta^{2}+1)\delta_{j}^{2}}\partial_{x}^{p}u_{j}\big{|}_{X_{j}}\,.

Premultiplying by $(x-X_{j})^{p}/p!$ provides the $p$ th term in a Taylor expansion of the smooth $u_{j}(t,x)$ about $x=X_{j}$ (the term $(x-X_{j})^{p}$ commutes with $\delta_{t}$ , $\delta_{j}$ and $\mu_{j}$ , but not with $\delta_{x}$ when $p>0$ ). Given $u_{j}\in\mathbb{C}_{j}$ , the Taylor series converges in some open interval containing $X_{j}$ , and so summing over all $p$ we obtain

4\sinh^{2}(H\sqrt{\mathcal{T}}/2)u_{j}(t,x)=\frac{\gamma^{2}(\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2})\delta_{j}^{2}}{1+\tfrac{1}{2}\gamma\delta_{j}^{2}-\tfrac{1}{4}\gamma^{2}(\theta^{2}+1)\delta_{j}^{2}}u_{j}(t,x)

(33)

in that interval. Then multiply by the weight function $w(x-X_{j})$ and integrate over the middle of the $j$ th element from $x=X_{j}-\ell$ to $X_{j}+\ell$ , and revert the $\sinh(-H\mathcal{T}/2c)$ operator, to deduce eq. 32a; or for the averaging width $\ell\rightarrow 0$ , set $x=X_{j}$ and revert the $\sinh$ operator to also deduce eq. 32a.

Finally, from eq. 32a revert the operator $\mathcal{T}=\mathcal{K}^{-1}(\partial_{t}-\mathcal{K}_{0})$ to deduce (32b). ∎

The above proof also shows that identity eq. 29 holds for $u_{j}$ evaluated at all $x\in(X_{j}-\ell,X_{j}+\ell)$ , not just $x=X_{j}$ . To prove this, in eq. 33 substitute $4\sinh^{2}(H\sqrt{\mathcal{T}}/2)=\delta_{x}^{2}$ . Similarly, the code of appendix A shows that the identity holds for all $\xi$ , not just $\xi=1/2$ .

With these lemmas we now prove the consistency of the slow manifold discretisation.

Theorem 17.

At full coupling $\gamma=1$ , and for every $\theta$ , the holistic discretisation (32) is consistent with the evolution equation (27).

For example, the holistic discretisation eq. 9c is consistent with the simple diffusion pde $u_{t}=u_{xx}$ .

Proof.

Evaluating (32a) at full coupling $\gamma=1$ gives (table 1)

	$\displaystyle\mathcal{T}U_{j}$	$\displaystyle=\frac{4}{H^{2}}\operatorname{asinh}^{2}\left[\sqrt{\frac{\mu_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2}}{1+\tfrac{1}{4}\delta_{j}^{2}-\tfrac{1}{4}\theta^{2}\delta_{j}^{2}}}\frac{\delta_{j}}{2}\right]U_{j}$
		$\displaystyle=\frac{4}{H^{2}}\operatorname{asinh}^{2}\left[{\tfrac{1}{2}\delta_{j}}\right]U_{j}=\frac{1}{H^{2}}\partial_{j}^{2}U_{j}\,.$

Let $U(t,x)$ be a smooth interpolation of $U_{j}(t)$ , then $\frac{1}{H}\partial_{j}U_{j}=\partial_{x}U$ , and so the above reduces to $\mathcal{T}U=\partial_{x}^{2}U$ , that is, $\mathcal{K}^{-1}(\partial_{t}-\mathcal{K}_{0})U=\partial_{x}^{2}U$ . By reverting the operator $\mathcal{K}^{-1}$ , it follows that $\partial_{t}U=\mathcal{K}_{0}U+\mathcal{K}(\partial_{x}^{2})U$ . Thus at $\gamma=1$ , both equations in (32) are consistent to the evolution equation (27) to all orders. ∎

Of more practical interest is what happens when the inter-element coupling is analysed to some finite order in $\gamma$ in order to construct an holistic discretisation of some finite stencil width in space.

Theorem 18.

Consider constructing, to errors $o\big{(}\gamma^{p}\big{)}$ , the spatial discretisation of the evolution equation (27) on elements coupled by eq. 6. The resulting holistic discretisation (e.g., eq. 9c truncated to errors $o\big{(}\gamma^{p}\big{)}$ ) is

•

a scheme of stencil width $(2p-1)$ on the macroscale grid, and
•

consistent with the general pde (27) to errors $o\big{(}\partial_{x}^{p}\big{)}$ .

Proof.

By section 3.3 the spatial discretisation of the evolution (27) satisfies eq. 32. Consider macroscale solutions of (32) for which the differences $\delta_{j}$ are ‘small’: for such solutions, the Taylor series in $\delta_{j}$ is appropriate. We expand $\operatorname{asinh}$ in eq. 32 but, similarly to the proof of theorem 12, only retain the powers of $\gamma$ , $\mu_{j}$ and $\delta_{j}$ , omitting all other details:

	$\displaystyle\operatorname{asinh}[\cdot]$	$\displaystyle\sim\operatorname{asinh}\left[\frac{\sqrt{(\mu_{j}^{2}-\delta_{j}^{2})}(\gamma\delta_{j})}{\sqrt{1+\delta_{j}(\gamma\delta_{j})+(\gamma\delta_{j})^{2}}}\right]$
		$\displaystyle\sim\sum_{m=0}^{\infty}(\mu_{j}^{2}-\delta_{j}^{2})^{m+1/2}(\gamma\delta_{j})^{2m+1}[1+\delta_{j}(\gamma\delta_{j})+(\gamma\delta_{j})^{2}]^{-m-1/2}$
		$\displaystyle\sim\sum_{m=0}^{\infty}(\mu_{j}^{2}-\delta_{j}^{2})^{m+1/2}(\gamma\delta_{j})^{2m+1}\sum_{n=0}^{\infty}\delta_{j}^{0:n}(\gamma\delta_{j})^{n:2n}$
		$\displaystyle\sim\sum_{m,n=0}^{\infty}(\mu_{j}^{2}-\delta_{j}^{2})^{m+1/2}\delta_{j}^{0:n}(\gamma\delta_{j})^{2m+n+1:2m+2n+1}\,.$

Thus, since $\mu_{j}\sim\delta_{j}$ the $\operatorname{asinh}^{2}$ in eq. 32 expands as

\operatorname{asinh}^{2}[\cdot]\sim\delta_{j}^{2}(\gamma\delta_{j})^{2}\left[\sum_{m,n=0}^{\infty}\delta_{j}^{2m:2m+n}(\gamma\delta_{j})^{2m+n:2m+2n}\right]^{2}.

First establish the theorem’s properties for (32a). For a given power $\delta_{j}^{p}$ , every $\gamma$ exponent is $\leq p$ . Hence, truncating the analysis to errors $o\big{(}\gamma^{p}\big{)}$ ensures that the terms are all complete and correct for every power of $\delta_{j}$ less than $p$ ; that is, the error is $o\big{(}\delta_{j}^{p}\big{)}$ . Since $\partial_{x}\sim\delta_{j}$ for smooth fields, the consistency error in such analysis is $o\big{(}\partial_{x}^{p}\big{)}$ .

The highest power of $\gamma$ retained is $p-1$ which could appear in the above sum as $\gamma^{2+2(2m+n)}$ . For such a power of $\gamma$ , the highest power of the difference is $\delta_{j}^{2+2+2(2m+n)+2(2m+n)}=\delta_{j}^{2(p-1)}$ , which involves shifts of $(E_{j}^{\pm 1/2})^{2(p-1)}=E_{j}^{\pm(p-1)}$ . That is, a construction to errors $o\big{(}\gamma^{p}\big{)}$ results in a macroscale discretisation with stencil width $(2p-1)$ .

Second, establish the theorem’s properties for (32b). Denote the the right-hand side operator of (32a) as $\mathcal{R}=\frac{4}{H^{2}}\operatorname{asinh}^{2}(\cdots)$ so that (32b) is $\partial_{t}U_{j}=\mathcal{K}_{0}U_{j}+\mathcal{K}\{\mathcal{R}\}U_{j}$ , and recall that $\mathcal{K}(\mathcal{R})=\sum_{k=1}^{\infty}\mathcal{K}_{k}\mathcal{R}^{k}$ . From the above analysis we have

\mathcal{R}^{k}\sim\operatorname{asinh}^{2k}[\cdot]\sim\delta_{j}^{2k}(\gamma\delta_{j})^{2k}\left[\sum_{m,n=0}^{\infty}\delta_{j}^{2m:2m+n}(\gamma\delta_{j})^{2m+n:2m+2n}\right]^{2k}.

As before, for a given power $\delta_{j}^{p}$ , every $\gamma$ exponent is $\leq p$ so truncating the to errors $o\big{(}\gamma^{p}\big{)}$ ensures the error is $o\big{(}\delta_{j}^{p}\big{)}$ and the consistency error is $o\big{(}\partial_{x}^{p}\big{)}$ .

When the highest power of $\gamma$ retained is $p-1$ , then this power could appear in $\mathcal{K}$ as $\gamma^{2k+2k(2m+n)}$ , corresponding to the highest power of the difference operator $\delta_{j}^{2k+2k+2k(2m+n)+2k(2m+n)}=\delta_{j}^{2(p-1)}$ . As before, this involves shifts of $(E_{j}^{\pm 1/2})^{2(p-1)}=E_{j}^{\pm(p-1)}$ , and so we conclude that a construction to errors $o\big{(}\gamma^{p}\big{)}$ results in a macroscale discretisation with stencil width $(2p-1)$ . ∎

4 Application to accurate numerical homogenisation

Consider a general heterogeneous diffusion pde, for field $u(t,x)$ in some spatial domain $\mathbb{X}$ ,

\mathchoice{\frac{\partial u}{\partial t}}{{\partial u}/{\partial t}}{{\partial u}/{\partial t}}{{\partial u}/{\partial t}}=\mathchoice{\frac{\partial}{\partial x}}{{\partial}/{\partial x}}{{\partial}/{\partial x}}{{\partial}/{\partial x}}\left[\kappa(x)\mathchoice{\frac{\partial u}{\partial x}}{{\partial u}/{\partial x}}{{\partial u}/{\partial x}}{{\partial u}/{\partial x}}\right],

(34)

where the diffusivity $\kappa(x)$ is $d$ -periodic in $x$ . The microscale length $d$ is fixed. We do not invoke the limit $d\to 0$ : $d$ is some constant that happens to be relatively small compared to the size $\mathfrak{L}$ of the the domain of interest. Our approach and results here apply to the physically relevant case of a finite scale separation ratio $d/\mathfrak{L}$ . The results are not restricted to the mathematical limit $d/\mathfrak{L}\to 0$ .

After some general considerations, we focus on the specific example case when the diffusivity is $\kappa=1/(1+a\cos kx)$ for some amplitude $a$ of the heterogeneity. For $d$ -periodic heterogeneity, here the wavenumber of the microscale is $k:=2\pi/d$ .

Figure 2: ‘cylindrical’ domain of the embedding pde eq. 35 for field

{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}(t,{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}},z)

. Obtain solutions of the heterogeneous diffusion pde eq. 34 on the blue line as

u_{\phi}(t,x)={\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}(t,x,x+\phi)

for every constant phase

\phi

To make rigorous progress we embed the general heterogeneous diffusion pde eq. 34 in the family of all phase shifts of the $d$ -periodic diffusion [[, e.g.,]]Roberts2013a, Roberts2016a. Consider the following pde for a field ${\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}(t,{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}},z)$ in the ‘cylindrical’ domain $\mathbb{X}\times[0,d)$ (fig. 2):

	$\displaystyle\partial_{t}{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}$	$\displaystyle=(\partial_{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}+\partial_{z})\big{[}\kappa(z)({\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}+{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{z})\big{]}$
		$\displaystyle=\partial_{z}[\kappa\partial_{z}]{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}+[2\kappa\partial_{z}+\kappa^{\prime}]{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}+\kappa{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}}\,,$		(35)

with a boundary condition of $d$ -periodicity in $z$ . For every field ${\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}(t,{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}},z)$ that satisfies pde (35), define the field $u_{\phi}(t,x):={\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}(t,x,x+\phi)$ where the $z$ -argument, $x+\phi$ , is taken modulo $d$ as indicated in fig. 2. Then from the pde eq. 35, the field $u_{\phi}$ satisfies the heterogeneous diffusion pde $\partial_{t}u_{\phi}=\partial_{x}[\kappa(x+\phi)\partial_{x}u_{\phi}]$ . Hence solutions ${\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}(t,{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}},z)$ of the embedding pde eq. 35 provide solutions of the original heterogeneous pde eq. 34 for every phase shift $\phi$ of the heterogeneity. For the specific phase shift $\phi=0$ , the field $u_{0}(t,x)$ solves the specific original pde eq. 34.

The crucial property of the embedding pde eq. 35 is that the pde is homogeneous in the longitudinal coordinate ${\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}$ . The heterogeneity in eq. 35 appears only in the cross-section coordinate $z$ . Consequently, theory and techniques developed for the case of systems homogeneous in $x$ apply to the ${\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}$ -dependence of the embedding pde eq. 35, and thence to the heterogeneous diffusion eq. 34. This phase-shift embedding also extends to nonlinear systems [[, e.g., §3.3]]Roberts2013a, but here we confine attention to linear homogenisation.

Specifically, theorem 18 indicates that a holistic discretisation constructed via inter-element coupling controlled by edge conditions (6), and constructed to errors $o\big{(}\gamma^{p}\big{)}$ , is a spatially discrete scheme with a stencil width of $(2p-1)$ that is consistent to the macroscale dynamics of the heterogeneous diffusion eq. 34 to errors $o\big{(}\partial_{x}^{p}\big{)}$ . The theorem only “indicates” because in order to apply, the theorem needs to address fields $u(t,x)$ in some Hilbert space of the $z$ -dependence.

Now consider the specific case of homogenising the specific heterogeneous diffusion $\kappa(x):=1/(1+a\cos kx)$ . With inter-element edge conditions (6) and tuning $\theta=0$ , appendix C lists code that constructs the holistic discretisation, the numerical homogenisation, of the heterogeneous diffusion eq. 34 in this case. The code implicitly assumes $d$ divides evenly into the element length $H$ as otherwise various coded means are not correct. For simplicity, we seek the results as a power series in both the coupling $\gamma$ and the amplitude $a$ of the heterogeneity.

In the ensemble of fig. 2, the flux in the ${\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}$ -direction (as distinct from the flux in the $x$ -direction) from the ensemble pde eq. 35 is $f=-2\kappa{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{z}-\kappa^{\prime}{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}-\kappa{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}}$ . Since ${\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}$ satisfies the edge condition eq. 6a, then so does its derivative in the transverse direction ${\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{z}$ . Hence, the flux edge condition eq. 6b, upon dividing by $-\kappa(z)$ , is satisfied by the derivative ${\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{\cal\mathchoice{\scriptstyle X}{\scriptstyle X}{\scriptscriptstyle X}{\scriptscriptstyle X}}$ . This is the second of the two edge conditions coded in appendix C.

Executing the code of appendix C finds that, to low-orders, the subgrid field

	$\displaystyle{\cal\mathchoice{\scriptstyle U}{\scriptstyle U}{\scriptscriptstyle U}{\scriptscriptstyle U}}_{j}$	$\displaystyle=U_{j}+\gamma\left\{\xi-\tfrac{1}{2}+\frac{a}{kH}\sin kz\right\}\mu_{j}\delta_{j}U_{j}\quad{}+\gamma^{2}\left\{\vphantom{\frac{1}{2}}(\tfrac{1}{12}-\tfrac{1}{2}\xi+\tfrac{1}{2}\xi^{2})\delta_{j}^{2}\right.$
		$\displaystyle\quad\left.{}+(\tfrac{1}{8}-\tfrac{1}{4}\xi)\mu_{j}\delta_{j}^{3}+(\tfrac{1}{48}-\tfrac{1}{8}\xi+\tfrac{1}{8}\xi^{2})\delta_{j}^{4}+\frac{a}{k^{2}H^{2}}\cos kz\,(\delta_{j}^{2}+\tfrac{1}{4}\delta_{j}^{4})\right.$
		$\displaystyle\quad\left.{}-\frac{a}{kH}\sin kz\,\left[(\tfrac{1}{2}-\xi)\delta_{j}^{2}+\tfrac{1}{4}\mu_{j}\delta_{j}^{3}+(\tfrac{1}{8}-\tfrac{1}{4}\xi)\delta_{j}^{4}\right]\right\}U_{j}+\mathcal{O}\mathchoice{\big{(}\gamma^{3},a^{3}\big{)}}{\big{(}\gamma^{3},a^{3}\big{)}}{(\gamma^{3},a^{3})}{(\gamma^{3},a^{3})}$	(36a)
This shows that the sub-element field has smooth macroscale structures in space through its $\xi$ -dependence, structures that are modified by the heterogeneity via the microscale diffusivity-variations represented by terms in $\sin kz$ and $\cos kz$ . Since their coefficients involve divisions by wavenumber $k$ , which are proportional to multiplication by the small microscale periodicity $d$ , these trigonometric modifications to the field are relatively small. Higher-order terms involve harmonics of these trigonometric functions.

Executing the code of appendix C to higher-order errors we find that the corresponding evolution modifies the homogeneous case eq. 9c to

\displaystyle H^{2}\partial_{t}U_{j}

\displaystyle=\cdots\lx@cref{creftype~refnum}{eqgopdiff}\cdots+\frac{a^{2}}{2k^{2}H^{2}}(\gamma^{4}-\gamma^{5}\delta_{j}^{2})\mu_{j}^{4}\delta_{j}^{4}U_{j}+\mathcal{O}\mathchoice{\big{(}\gamma^{6},a^{7}\big{)}}{\big{(}\gamma^{6},a^{7}\big{)}}{(\gamma^{6},a^{7})}{(\gamma^{6},a^{7})},

(36b)

where the microscale heterogeneity only affects the macroscale evolution through the magnitude $a$ and the wavenumber $k$ of the heterogeneity. Alternatively, executing the code to errors $\mathcal{O}\mathchoice{\big{(}\gamma^{7},a^{3}\big{)}}{\big{(}\gamma^{7},a^{3}\big{)}}{(\gamma^{7},a^{3})}{(\gamma^{7},a^{3})}$ , and evaluated at full coupling $\gamma=1$ , we compute that this evolution has equivalent pde

\displaystyle\mathchoice{\frac{\partial U}{\partial t}}{{\partial U}/{\partial t}}{{\partial U}/{\partial t}}{{\partial U}/{\partial t}}

\displaystyle=\mathchoice{\frac{\partial^{2}U}{\partial x^{2}}}{{\partial^{2}U}/{\partial x^{2}}}{{\partial^{2}U}/{\partial x^{2}}}{{\partial^{2}U}/{\partial x^{2}}}+\frac{a^{2}}{2k^{2}}\mathchoice{\frac{\partial^{4}U}{\partial x^{4}}}{{\partial^{4}U}/{\partial x^{4}}}{{\partial^{4}U}/{\partial x^{4}}}{{\partial^{4}U}/{\partial x^{4}}}-\frac{2a^{2}}{k^{4}}\mathchoice{\frac{\partial^{6}U}{\partial x^{6}}}{{\partial^{6}U}/{\partial x^{6}}}{{\partial^{6}U}/{\partial x^{6}}}{{\partial^{6}U}/{\partial x^{6}}}+\mathcal{O}\mathchoice{\big{(}a^{2}H^{2}d^{4}\partial_{x}^{8}U\big{)}}{\big{(}a^{2}H^{2}d^{4}\partial_{x}^{8}U\big{)}}{(a^{2}H^{2}d^{4}\partial_{x}^{8}U)}{(a^{2}H^{2}d^{4}\partial_{x}^{8}U)}.

(36c)

The leading-order homogenisation, the diffusion term $\mathchoice{\frac{\partial^{2}U}{\partial x^{2}}}{{\partial^{2}U}/{\partial x^{2}}}{{\partial^{2}U}/{\partial x^{2}}}{{\partial^{2}U}/{\partial x^{2}}}$ , is exactly the well-known correct harmonic mean of the diffusivity $\kappa(x)$ . The higher-order terms, $\mathchoice{\frac{\partial^{4}U}{\partial x^{4}}}{{\partial^{4}U}/{\partial x^{4}}}{{\partial^{4}U}/{\partial x^{4}}}{{\partial^{4}U}/{\partial x^{4}}}$ , being divided by powers of $k=2\pi/d$ , vanish in the usual theoretical limit of $d\to 0$ . However, here our analysis is rigorous for finite $d$ and so these fourth and higher order terms quantify effects due to the physical finite scale separation of the macroscale from a finite sized microscale. These higher-order derivatives depend upon $a^{2}$ (recall that $a$ is the magnitude of the heterogeneity), as appropriate by symmetry in $a$ .

Consequently, the holistic discretisation (36b) is an accurate, analytically learnt, numerical homogenisation for the heterogeneous diffusion (34).

Numerical homogenisation for waves

Wave propagation (section 2.3) through heterogeneous media can be analysed almost identically to the heterogeneous diffusion of this section. The resultant numerical homogenisation would be (36b) but with $\partial_{t}^{2}U_{j}$ on the left-hand side. Its equivalent pde would be (36c) but with $\mathchoice{\frac{\partial^{2}U}{\partial t^{2}}}{{\partial^{2}U}/{\partial t^{2}}}{{\partial^{2}U}/{\partial t^{2}}}{{\partial^{2}U}/{\partial t^{2}}}$ on the left. Consequently, the numerical homogenisation would predict, through the fourth- and sixth-order differences/derivatives in (36), a wave speed dependence upon wavelength (wave dispersion) caused by the microscale heterogeneity at finite scale separation.

Spatial boundaries

The homogenisation (36b) is constructed with periodic diffusivity and is independent of the phase $\phi$ of the microscale heterogeneity. Macroscale effects of the phase $\phi$ only arise via correctly determined influences of the boundary conditions on the macroscale domain. The development of homogenisation for general boundary conditions is the subject of ongoing research. For spatial discretisations, such as (36b), one could develop correct discretisations near a boundary by adapting the arguments of [Roberts01b, MacKenzie03].

5 Conclusion

Fine-scale heterogeneity or nonlinearities both complicate the derivation of macroscale discretisation of pdes, with many common methods failing to accurately account for the effects at the macroscale of subgrid physics. Holistic discretisations systematically construct macroscale closures which can accurately account for subgrid structures. In particular, for an accurate homogenisation of a pde, the discretisation must maintain symmetries of pde. Here we developed a new holistic discretisation which is guaranteed to maintain the self-adjointness of a pde through carefully crafted inter-element coupling conditions, thus ensuring that the general spectral structure of the pde is captured by the discretisation. For both generic 1D waves and diffusion, we show that the homogenisation is consistent with the original pde.

In this article we focus on 1D systems, but simulations of 2D systems indicate that analogous discretisations are possible in 2D. Such 2D systems will be explored and theory developed in future research.

Acknowledgement

This research was funded by the Australian Research Council under grants DP150102385 and DP200103097. We thank Peter Hochs for his comments on earlier versions of this article.

Appendix A High-order consistency for holistic discretisation of advection-diffusion

This is script diffAdvecHolistic.tex— Learns the holistic discretisation of advection-diffusion on elements with self-adjoint preserving coupling. Find on slow manifold both $u$ and $dU/dt$ linear in advection $c$ , although the code must truncate in $c$ in order to converge. All code is written in the computer algebra package Reduce.⁷⁷7Reduce [http://reduce-algebra.com/] is a free, fast, general purpose, computer algebra system.

1 
2 on div; off allfac; on revpri;
3 factor gamma,hh,c;

Subgrid structures are functions of sub-element variable $\xi=(x-L_{j})/H$

4 
5 depend xi,x;  let df(xi,x)=>1/hh;

Operator $u(1)=u(0)$ and mean zero.

6 
7 operator linv; linear linv;
8 let { linv(xi^~~p,xi)=>(xi^(p+2)-xi+1/2-1/(p+3))/(p+1)/(p+2)
9     , linv(1,xi)=>(xi^2-xi+1/2-1/3)/2 };

Operator to compute mean over an element

10 
11 operator mean; linear mean;
12 let { mean(xi^~~p,xi)=>1/(p+1)
13     , mean(1,xi)=>1 };

Parametrise discretisation slow manifold with evolving order parameters

14 
15 operator uu;  depend uu,t;
16 let df(uu(~k),t)=>sub(j=k,gj);

Initial approximation in $j$ th element

17 
18 uj:=uu(j);  gj:=0;

Here specify required orders of errors in the result. Intermediate working needs $c$ truncated to some order

19 
20 let { gamma^3=>0, c^3=>0};

Deep iterative refinement learns emergent slow manifold

21 
22 R:= xi=1;  L:= xi=0;
23 for iter:=1:99 do begin

Compute residuals for selected one of possible PDEs, the two coupling conditions, and the definition of the order-parameter macroscale variable.

24 
25    pde:= -df(uj,t) + (if 1
26      then -c*df(uj,x)+df(uj,x,2)
27      else (sub(xi=xi+d,uj)-2*uj+sub(xi=xi-d,uj))/(d*hh)^2 );
28    ucc:= -(1-gamma/2)*(sub(R,uj)-sub(L,uj))
29      +gamma/2*(sub({L,j=j+1},uj)-sub({R,j=j-1},uj))
30      +gamma*theta/2*(sub(R,uj)+sub(L,uj))
31      -gamma*theta/2*(sub({L,j=j+1},uj)+sub({R,j=j-1},uj));
32    ux:=df(uj,x);
33    udc:= -(1-gamma/2)*(sub(R,ux)-sub(L,ux))
34      +gamma/2*(sub({L,j=j+1},ux)-sub({R,j=j-1},ux))
35      -gamma*theta/2*(sub(R,ux)+sub(L,ux))
36      +gamma*theta/2*(sub({L,j=j+1},ux)+sub({R,j=j-1},ux));
37    amp:=mean(uj,xi)-uu(j);

Trace write lengths of residuals.

38 
39    write lengthress:=map(length(~a),{pde,ucc,udc,amp});

Update approximations from the current residuals

40 
41    gj:=gj+(gd:=udc/hh-mean(pde,xi));
42    uj:=uj-hh^2*linv(pde-gd,xi)+(xi-1/2)*ucc;

Exit loop when residuals are zero to specified order of error

43 
44    if {pde,ucc,udc,amp}={0,0,0,0}
45    then write "Success: ",iter:=100000+iter;
46 end;%for
47 if {pde,ucc,udc,amp}neq{0,0,0,0} then rederr("iteration fail");

Compute and report equivalent pde of discrete

48 
49 in_tex "../convert2EquivalentDE.tex"$

Optionally check some identities, and exit script.

50 
51 in "diffTestProof.txt"$
52 end;

A.1 Convert evolution to equivalent PDE

All algorithms invoke this conversion script

53 
54 write "Convert evolution dUjdt=gj to equivalent PDE for
55 U(x,t) via operator form.  AJR, from a long time ago.";

Convert to central difference operator form

56 
57 rules:={ mu^2=>1+delta^2/4, uu(j)=>1
58     , uu(j+~p)=>(1+sign(p)*mu*delta+delta^2/2)^abs(p)}$
59 gop:=(gj where rules);
60 uop:=(uj where rules)$

Convert to equivalent pde using Taylor expansion

61 
62 remfac gamma; factor df;
63 let hh^10=>0;
64 depend uu,x;
65 rules:={uu(j)=>uu, uu(j+~p)=>uu+(for n:=1:10 sum
66                df(uu,x,n)*(hh*p)^n/factorial(n)) }$
67 duujdt:=(gj where rules);
68 duujdt1:=sub(gamma=1,duujdt);

69 
70 end;

Appendix B High-order consistency for holistic discretisation of simple wave PDE

This is script Holistic discretisation of modified uni-directional wave PDE, with speed c, on elements with ‘self-adjoint’ coupling. The sub-element field is independent of speed c. For error $\mathcal{O}\mathchoice{\big{(}\gamma^{p}\big{)}}{\big{(}\gamma^{p}\big{)}}{(\gamma^{p})}{(\gamma^{p})}$ gives consistency $\mathcal{O}\mathchoice{\big{(}H^{p-1}\big{)}}{\big{(}H^{p-1}\big{)}}{(H^{p-1})}{(H^{p-1})}$ . For linear modifications measured by alpha, the slow manifold is linear in alpha, but must truncate in alpha for iteration to terminate!

71 
72 on div; off allfac; on revpri;
73 factor gamma,hh,c,alpha,c0,c2,c3,c4;

Subgrid structures are functions of sub-element variable $\xi=(x-L_{j})/H$

74 
75 depend xi,x;  let df(xi,x)=>1/hh;

Operator

76 
77 operator linv; linear linv;
78 let { linv(xi^~~p,xi)=>(xi^(p+1)-1/(p+2))/(p+1)
79     , linv(1,xi)=>(xi-1/2) };

Operator to compute mean over an element

80 
81 operator mean; linear mean;
82 let { mean(xi^~~p,xi)=>1/(p+1)
83     , mean(1,xi)=>1 };

Parametrise discretisation slow manifold with evolving order parameters

84 
85 operator uu;  depend uu,t;
86 let df(uu(~k),t)=>sub(j=k,gj);

Initial approximation in $j$ th element

87 
88 uj:=uu(j);  gj:=0;

Here specify required orders of errors in the result.

89 
90 let { gamma^6=>0, alpha=>0 };

Deep iterative refinement learns emergent slow manifold

91 
92 R:= xi=1;  L:= xi=0;
93 for iter:=1:99 do begin

Compute residuals for selected one of possible PDEs, the coupling condition, and the definition of the order-parameter macroscale variable.

94 
95    pde:= -df(uj,t) + (if 1
96      then -c*df(uj,x)
97        +alpha*(c0*uj+c2*df(uj,x,2)+c3*df(uj,x,3)+c4*df(uj,x,4))
98      else -c*(sub(xi=xi+d,uj)-sub(xi=xi-d,uj))/(2*d*hh) );
99    fx:=c*uj; the ’flux’
100    ucc:=(1+theta)/2*( -sub(L,fx)
101            +(1-gamma)*sub(R,fx) +gamma*sub({R,j=j-1},fx) )
102        -(1-theta)/2*(-sub(R,fx)
103            +(1-gamma)*sub(L,fx) +gamma*sub({L,j=j+1},fx) );
104    amp:=mean(uj,xi)-uu(j);

Trace write lengths of residuals.

105 
106    write lengthress:=map(length(~a),{pde,ucc,amp});

Update approximations from the current residuals

107 
108    gj:=gj+(gd:=ucc/hh+0*mean(pde,xi));
109    uj:=uj+hh/c*linv(pde-gd,xi);

Exit loop when residuals are zero to specified order of error

110 
111    if {pde,ucc,amp}={0,0,0}
112    then write "Success: ",iter:=100000+iter;
113 end;%for
114 if {pde,ucc,amp}neq{0,0,0} then rederr("iteration fail");

Compute and report equivalent pde, then exit script.

115 
116 in_tex "../convert2EquivalentDE.tex"$
117 end;

Appendix C Holistic discretisation for heterogeneous diffusion and its consistency

This is script Holistic discretisation of heterogeneous diffusion on elements with self-adjoint preserving coupling. Here seek discretisation as series in $a$ , the amplitude of the heterogeneity.

118 
119 on div; off allfac; on revpri;
120 factor gamma,hh,a,k;

Subgrid structures are functions of sub-element variable $\xi=(x-L_{j})/H$

121 
122 depend xi,x;  let df(xi,x)=>1/hh;
123 depend xi,xz; depend z,xz;

Operator $u(1)=u(0)$ , periodic in $z$ , and mean zero.

124 
125 operator linv; linear linv;
126 let { linv(xi^~~p,xz)=>hh^2*(xi^(p+2)-xi+1/2-1/(p+3))/(p+1)/(p+2)
127     , linv(1,xz)=>hh^2*(xi^2-xi+1/2-1/3)/2
128     , linv(sin(~~q*z),xz)=>-1/(q)^2*sin(q*z)
129     , linv(cos(~~q*z),xz)=>-1/(q)^2*cos(q*z)
130     , linv(sin(~~q*z)*xi^~~p,xz)=>-(xi^p-xi)/(q)^2*sin(q*z)
131     , linv(cos(~~q*z)*xi^~~p,xz)=>-(xi^p-xi)/(q)^2*cos(q*z)
132     };

Operator to compute mean over an element, these assume trig functions fit in period

133 
134 operator mean; linear mean;
135 let { mean(1,xz)=>1 , mean(xi^~~p,xz)=>1/(p+1)
136     , mean(cos(~~q*z),xz)=>0 , mean(xi^~~p*cos(~~q*z),xz)=>0
137     , mean(sin(~~q*z),xz)=>0 , mean(xi^~~p*sin(~~q*z),xz)=>0
138     };

Parametrise discretisation slow manifold with evolving order parameters

139 
140 operator uu;  depend uu,t;
141 let df(uu(~k),t)=>sub(j=k,gj);

Initial approximation in jth element

142 
143 uj:=uu(j);  gj:=0;

Iterative refinement to specified order of error, and using Taylor series of $\kappa(x)$ in coefficient $a$

144 
145 let { gamma^7=>0, a^3=>0 };
146 kappa:=for n:=0:deg((1+a)^9,a) sum (-a*cos(k*z))^n;
147 for iter:=1:99 do begin

Compute residuals for the PDE, the two coupling conditions, and the definition of the order-parameter macroscale variable.

148 
149    flux:=kappa*(df(uj,x)+df(uj,z));
150    pde:=trigsimp( -df(uj,t)+df(flux,x)+df(flux,z)
151        ,combine);
152    ucc:=(1-gamma/2)*(-sub(xi=+1,uj)+sub(xi=0,uj))
153    +gamma/2*(sub({xi=0,j=j+1},uj)-sub({xi=1,j=j-1},uj));
154    ux:=df(uj,x);
155    udc:=(1-gamma/2)*(-sub(xi=+1,ux)+sub(xi=0,ux))
156    +gamma/2*(sub({xi=0,j=j+1},ux)-sub({xi=1,j=j-1},ux));
157    amp:=mean(uj,xz)-uu(j);

Trace write lengths of residuals.

158 
159    write lengthress:=map(length(~a),{pde,ucc,udc,amp});

Update approximations from the current residuals

160 
161    gj:=gj+(gd:=mean(udc,xz)/hh+mean(pde,xz));
162    uj:=uj-linv(pde-gd,xz)+(xi-1/2)*ucc;

Exit loop when residuals are zero to specified order of error

163 
164    if {pde,ucc,udc,amp}={0,0,0,0}
165    then write "Success: ",iter:=100000+iter;
166 end;
167 if {pde,ucc,udc,amp}neq{0,0,0,0} then rederr("iteration fail");

Compute and report equivalent pde, then exit script.

168 
169 in_tex "../convert2EquivalentDE.tex"$
170 end;

$\displaystyle\operatorname{i}\partial_{t}U_{j}$	$\displaystyle=\alpha\|U_{j}\|^{2}U_{j}-\gamma^{2}\tfrac{1}{8H^{2}}(U_{j+2}+U_{j-2}-2U_{j})$
	$\displaystyle\quad{}+\alpha\gamma^{2}\tfrac{1}{96}\big{[}2\|U_{j}\|^{2}(U_{j+2}+U_{j-2}-4U_{j})$
	$\displaystyle\qquad{}-4\|U_{j+2}\|^{2}(U_{j+2}-2U_{j})-4\|U_{j-2}\|^{2}(U_{j-2}-2U_{j})$
	$\displaystyle\qquad{}-4U_{j}^{}(U_{j+1}^{2}+U_{j-1}^{2}-U_{j+1}U_{j-1})-4U_{j}(U_{j+1}^{}U_{j-1}+U_{j-1}^{*}U_{j+1})$
	$\displaystyle\qquad{}+U_{j+2}^{}(U_{j}^{2}-2U_{j+1}^{2})+U_{j-2}^{}(U_{j}^{2}-2U_{j-1}^{2})\big{]}+\mathcal{O}\mathchoice{\big{(}\alpha^{2},\gamma^{3}\big{)}}{\big{(}\alpha^{2},\gamma^{3}\big{)}}{(\alpha^{2},\gamma^{3})}{(\alpha^{2},\gamma^{3})},$	(5)