A toolbox of Equation-Free functions in Matlab\Octave for efficient system level simulation
NNoname manuscript No. (will be inserted by the editor)
A toolbox of Equation-Free functions inMatlab/Octave for efficient system level simulation
John Maclean · J. E. Bunder · A. J. Roberts
April 8, 2020
Abstract
The ‘equation-free toolbox’ empowers the computer-assisted analysisof complex, multiscale systems. Its aim is to enable you to immediately usemicroscopic simulators to perform macro-scale system level tasks and analysis,because micro-scale simulations are often the best available description of asystem. The methodology bypasses the derivation of macroscopic evolutionequations by computing the micro-scale simulator only over short bursts in timeon small patches in space, with bursts and patches well-separated in time and
John MacleanSchool of Mathematical Sciences, University of Adelaide, South Australia.
J. E. BunderSchool of Mathematical Sciences, University of Adelaide, South Australia. mailto:[email protected] , http://orcid.org/0000-0001-5355-2288 · A. J. RobertsSchool of Mathematical Sciences, University of Adelaide, South Australia. , http://orcid.org/0000-0001-8930-1552 a r X i v : . [ c s . M S ] A p r John Maclean et al. space respectively. We introduce the suite of coded equation-free functions inan accessible way, link to more detailed descriptions, discuss their mathematicalsupport, and introduce a novel and efficient algorithm for Projective Integration.Some facets of toolbox development of equation-free functions are then detailed.Download the toolbox functions and use to empower efficient and accuratesimulation in a wide range of your science and engineering problems. Keywords
Multiscale methods · code toolbox · numerical algorithms Contents
Suppose that you have a detailed and trustworthy computational simulationof some problem of interest. When the detailed computation is too expensiveto simulate all the times of interest over all the space of interest, then the‘Equation-Free Methodologies’ aim to accurately empower long-time simulationand system level analysis [e.g., Kevrekidis and Samaey, 2009, Kevrekidis et al., https://github.com/uoa1184615/EquationFreeGit toolbox of Equation-Free functions 3 Figure 1
Snapshots in time of the patch scheme simulation of the nonlinear diffusivemicro-scale system (1).(a) t = 0 (b) t = 0 . − − . space x space y u ( x , y ) − − . space x space y u ( x , y ) (c) t = 1 (d) t = 3 − − . space x space y u ( x , y ) − − . space x space y u ( x , y ) Matlab /Octave functions.1.1 Simulation on only small patches of spaceFigure 1 illustrates an ‘equation-free’ computation on only small well-separatedpatches of the spatial domain. The micro-scale simulations within each patch,here the nonlinear diffusive system (1), are craftily coupled to neighbouringpatches and thus interact to provide accurate macro-scale predictions over the
April 8, 2020
John Maclean et al. whole spatial domain [e.g., Roberts et al., 2014]. We have proved that thepatches may be tiny, and still the scheme makes accurate macro-scale predic-tions [Roberts and Kevrekidis, 2007]. Thus the computational savings may beenormous, especially when combined with projective integration (Section 1.2).The example system illustrated in Figure 1 is a nonlinear discrete diffusionsystem inspired by the lubrication flow of a thin layer of fluid, namely ∂u∂t = ∇ · (3 u ∇ u ) , equivalently ∂u∂t = ∇ ( u ) , which on a 2D micro-scale lattice x i,j with tiny micro-scale spacing d is herediscretised simply to du i,j dt = u i +1 ,j + u i − ,j + u i,j +1 + u i,j − − u i,j d . (1)We want to predict the dynamics of this spatial micro-scale lattice systemon the macro-scale spatial domain [ − , × [ − , , but suppose full directcomputation is too expensive. Instead, the micro-scale simulation illustratedby Figure 1 was performed only on about 20% of the domain (it could bemuch less)—the small patches of space in Figure 1. The key to an accuratemacro-scale prediction is that each patch is coupled to nearby patches, at everycomputed time, by appropriate macro-scale interpolation that gives the edgevalues for every patch [e.g., Roberts et al., 2014].The patch scheme is most useful in applications where there is no knownmacro-scale closure. Then the patch scheme automatically achieves a compu-tational macro-scale closure, without the need for any analytic constructionoften invoked in numerical/computational homogenization [e.g., Saeb et al., toolbox of Equation-Free functions 5 x i are coordinatesof a micro-scale lattice, for potentially exhaustingly many lattice points in-dexed by i ; for example, a full atmospheric simulation. And suppose yourdetailed and trustworthy simulation is coded in terms of micro-field variablevalues u i ( t ) ∈ R p at lattice point x i at time t . When a detailed computa-tional simulation is prohibitively expensive over all the desired spatial domain, x ∈ X ⊂ R d , our toolbox provides functions that empower you to use yourmicro-scale code as a ‘black-box’ inside only small, well-separated, patchesof space by appropriately coupling across un-simulated space between thepatches (Section 2.2). The toolbox functions have many options including bothnewly developed spectral coupling, and new symmetry preserving coupling.Section 3.2 gives an introductory tutorial.1.2 Projective Integration skips through timeSimulation over time is a complementary dynamic problem. The ‘equation-free’approach is to simulate for only short bursts of time, and then to extrapolateover un-simulated time into the future, or into the past, or perform systemlevel analysis [e.g., Gear and Kevrekidis, 2003b, Rico-Martinez et al., 2004,Erban et al., 2006, Givon et al., 2006]. Figure 2 plots one example where thegaps in time show the un-computed times between bursts of computation. April 8, 2020
John Maclean et al.
Figure 2
Projective Integration by the new function
PIG of the example multiscale system (2).The macro-scale solution U ( t ) is represented by the blue circles ( ◦ ). The black dotted line,underneath the PI solution, shows an accurate micro-scale simulation over the whole timedomain. . . . . . time t PI: U ( t ) = u ( t ) PI: u ( t ) micro burstsPI: u ( t ) micro burstsSimulation with ode45() In general such a simulation is coded in terms of detailed (micro-scale)variable values u ( t ) , in R p for some p (typically large in applications), andevolving in time t . The details u could represent particles, agents, or states ofa system. Both forward and backward in time computations may be performedby Projective Integration (PI) with provable accuracy [Gear and Kevrekidis,2003b, Givon et al., 2006, Maclean and Gottwald, 2015]. For efficient simulationon long times, Section 2.1 describes how to provide your micro-scale detailed Matlab /Octave code as a ‘black box’ to the novel Projective Integrationfunctions in the toolbox. Section 3.1 gives a user-friendly introductory tutorial. toolbox of Equation-Free functions 7
The example simulation of Figure 2 is that of a toy system that nonethelesshas challenging qualities of the multiscale phenomena that Projective Integra-tion resolves. Here the micro-scale simulation is a pair of coupled slow-fast ode s for u ( t ) = ( u ( t ) , u ( t )) : du dt = cos( u ) sin( u ) cos( t ) , (2a) du dt =10 (cid:2) cos( u ) − u (cid:3) . (2b)Using simple integration schemes, numerical solutions can be rapidly computedon micro-times of O (cid:0) − (cid:1) , but solutions over O (cid:0) (cid:1) times are computationallyprohibitive—except by stiff integrators. Section 3.1.4 discusses scenarios wherestiff integrators cannot be used or are relatively expensive, but where projectiveintegration is effective.The system (2) represents the realities of, for example, molecular dynamicssimulations with rapid modes represented by u ( t ) and slow macro-scale statevariables (like temperature) represented by u ( t ) . In applications the dynam-ics of the slow modes are usually not known, instead they emerge over themicro-scale simulation bursts [e.g., Cisternas et al., 2004, Setayeshgar et al.,2005, Erban et al., 2006]. Consequently, although in this toy system the slowvariable u and fast variable u are obvious, here we compute only with thefull system (2)—it is a ‘black-box’ for which we do not necessarily know what‘variables’ are fast or slow. An alternative is to invoke algebraic analysis toconstruct the slow manifold of slow-fast systems like (2) [e.g., Roberts, 2015,Ch. 4–5]: however, in many scenarios such analysis is not feasible. Since oneoften only measures macro-scale state variables, here we suppose the micro- April 8, 2020
John Maclean et al. simulator only outputs U ( t ) = u ( t ) , called a restriction. Projective Integration(PI) uses only short bursts of the ‘black-box’ simulation, and then invokes anappropriate extrapolation to accurately predict large ‘projective’ steps. After aprojective step the PI algorithm appropriately re-initialises u ( t ) , called lifting,for another burst of the micro-simulator.Figure 2 shows the PI computation and the associated micro-scale bursts,compared to expensive simulation with ode45() . The micro-scale simulationswith u ( t ) show the characteristic fast transients. The absolute error betweenthe Projective Integration simulation and the trusted ode45() simulationis · − .The PI simulation of Figure 2 was performed by a novel function, called PIG , introduced in Section 2.1 and included in the toolbox. The advantage ofthe Projective Integration simulation over
Matlab ’s ode45() is that PI usesonly . of the number of evaluations of the ode s (2). Section 3.1 discusseshow this new PI scheme and function may be constructed through recursiveuse of ode45() which inherits all of Matlab ’s adaptive error control—an errorcontrol active on both the micro-scale and the macro-scale.
This section outlines the key Projective Integration and Patch Dynamicsalgorithms implemented in the toolbox [Roberts et al., 2020]. Pseudo-codehighlights the essential features of each algorithm and the accompanying toolbox of Equation-Free functions 9 discussion demonstrates the extended capabilities. Theoretical support for thealgorithms is also discussed.2.1 Algorithms for projective integration in timeOur efficient simulation of the stiff dissipative system (2) is due to the classof algorithms called Projective Integrators. Small time steps with the micro-simulator are alternated with long time steps comprised of ‘projective’ extrap-olations. Sometimes PI is done with the macrostep being a Forward Eulermethod [Siettos et al., 2003, Chuang et al., 2015, e.g.], and so incurs globalerror proportional to the size of the projective time steps. Such low accuracymotivated the development of PI algorithms of provably higher accuracy [Rico-Martinez et al., 2004, Lee and Gear, 2007]. Here we additionally extend themethodology to a projective integration that invokes adaptive integrators.The toolbox provides
PIG() —denoting a Projective Integration with Gen-eral macro-integrator and micro-simulator. The user specified macro-integratoris any integrator suitable for time stepping on the macro-scale, and includesadaptive integrators such as
Matlab ’s ode45 as used for Figure 2. Then PIG() provides to that integrator accurate time derivatives at the requiredtimes—time derivatives estimated from appropriate relatively short bursts ofthe user-defined micro-simulator microSim() . The effect is to do PI with anyintegrator, explicit or implicit, taking the projective steps. Error analyses sug-gested that schemes of this sort usually incur significant errors proportional tothe duration of the micro-simulator burst [E, 2003, Maclean and Gottwald, 2015,
April 8, 2020
Maclean, 2015]. However, such errors are avoided by
PIG() through a novelimplementation of a ‘constraint-defined manifold computing’ scheme based ona methodology originally proposed by Gear et al. [2005], Gear and Kevrekidis[2005]. Algorithm 1 provides pseudo-code for
PIG() that details these steps.This package is the first time such functionality for Projective Integration hasbeen developed into a general function, tested, and made available.Algorithm 1 outlines the essence of
PIG() , but additional features may beinvoked by a user. In particular, sometimes the projective steps are performedon a few macro-scale variables only. That is, the PI is done in a space ofreduced dimension to that of the micro-simulator [Frederix et al., 2007, Boldet al., 2012, Sieber et al., 2018, e.g.]. In such cases the user provides ‘lifting’and ‘restriction’ operators to convert between the micro- and macro-simulationspaces [Roose et al., 2009].In addition to
PIG() , the toolbox provides efficient integrators
PIRK2() and
PIRK4() , which are Projective Integrators similar to
PIG() but with the user-defined macro-integrator replaced with second and fourth order, respectively,Runge–Kutta macro-integrators that take user specified macro-scale time-steps.2.2 Algorithms to simulate on patches of spaceThe spatial multiscale ode s (1) are simulated only on a fraction of space byemploying a patch scheme [also known as the gap-tooth scheme, Samaey et al.,2010]. This scheme applies to dynamic systems evolving in time with some‘spatial’ structure (‘space’ could be some other type of domain), and is useful toolbox of Equation-Free functions 11
Algorithm 1
The
PIG () function uses short bursts of a user provided micro-scale process to simulate from time T to T , with initial condition u , and viaa user-specified macro-scale integrator macroInt (). PIG tells the specified macro-scale integrator to use time derivatives estimated by a constraint-defined manifold function.1: function
PIG ( macroInt () , microSim () , constrDeriv () , { T , T } , u )2: return { T n , u n } Nn =0 := macroInt ( constrDeriv () , { T , T } , u ); end function Here, the toolbox’s constraint-defined manifold function uses two applications of microSim and a backwards projective step in order to provide the time derivative at the time specifiedby macroInt ().1: function constrDeriv ( t , u , δ )2: { t m , u m } Mm =0 := microSim ( t , u , δ ); (burst of simulation)3: (cid:16) ∆ u ∆t (cid:17) M := u M − u M − t M − t M − ; (derivative estimate at t + δ ) v := u m − δ (cid:16) ∆ u ∆t (cid:17) M ; (PI backwards to t − δ ) { s m , v m } Mm =0 := microSim ( t − δ , v , δ ); (burst of simulation)6: return (cid:16) ∆ v ∆s (cid:17) M := v M − v M − s M − s M − ; (derivative estimate at t ) end function A user provides a micro-scale simulator with the following input and output.1: function microSim ( t , u , δ )2: compute a burst with end time t M := t + δ ;3: return { t m , u m } Mm =0 end function April 8, 2020
Algorithm 2
Patch scheme when a user has a function fun ( t, u ) that computesa time-step/time-derivative on a micro-scale lattice in space. Suppose thatdirect simulation is prohibitively expensive with fun by some time integrator,say macroInt ( fun , { T , T } , u ) . A user then invokes the patch scheme withthe following two steps.
1: Invoke configPatches1 ( fun , { a , b } , . . . ) to configure the patch scheme on the spatialdomain [ a , b ] of interest (Algorithm 3).2: Then executing macroInt ( patchSmooth1 , { T , T } , u ) computes the solution on patchesin space (Algorithm 4). when direct simulation over the spatial domain of interest is infeasible due tothe overwhelming computational cost. The toolbox assumes that space on themicro-scale is represented as a rectangular lattice network with grid points x i —one example being the spatial discretization (1) of a pde . The toolbox functionrequires the user to provide a function that computes a micro-scale time-stepor time derivative for variables u i at each of the lattice nodes.The patch scheme computes the detailed micro-scale at only a (small) subsetof lattice nodes in space—the patches as illustrated by Figure 1. The schemecraftily couples the patches together to ensure macro-scale predictions areprovably accurate [Roberts and Kevrekidis, 2007, Roberts et al., 2014, Bunderet al., 2017, Cao and Roberts, 2016]. For simplicity, this article discusses specif-ically the case of patches in 1D space, but the scheme and the toolbox extendto higher dimensional space [Roberts et al., 2014], as in the 2D simulationsof Figure 1. The essence of the toolbox patch algorithms are outlined, for aspecific case in 1D space, by Algorithms 2 to 4. toolbox of Equation-Free functions 13 Algorithm 3 configPatches1 initialises the basic patch structures for Algo-rithm 2. It creates nPatch equi-spaced patches inside the 1D interval [ a , b ] . Thepatches are of size proportional to ratio with each having nSubP micro-scalelattice points, and coupled by interpolation of order ordCC . All defined variablesare saved as global for subsequent use in simulation. function configPatches1 ( fun , {a,b}, nPatch , ratio , ordCC , nSubP )2: H := ( b − a ) /N ; (distance between macro-scale mid-patches)3: X j := a + ( j − ) H for j := 1 : nPatch ; (macro-scale mid-patches)4: i0 := ( nSubP + 1) / ; (index for the centre point of a patch)5: d := ( ratio × H ) / ( i0 − ; (micro-scale lattice spacing)6: x i,j := X j + d ( i − i0 ) for j := 1 : nPatch and i := 1 : nSubP ; (micro-scale latticepoints in the j th patch)7: Calculate coefficients of interpolation to determine the edge-values of a patchfrom mid-patch values, for interpolation from ordCC / patches to either side ofa patch.8: end function Algorithm 2 introduces the necessary code for using the toolbox to imple-ment the patch scheme on a problem of interest. In the specific case of a macro-scale 1D domain [ a, b ] , the toolbox function configPatches1 (Algorithm 3)constructs nPatch equi-spaced patches, each of width proportional to ratio and each containing nSubP micro-scale lattice points. The integer ordCC deter-mines the form and order of the inter-patch coupling that ensures macro-scaleaccuracy: – even ordCC , the usual, and as described in Algorithm 3, invokes clas-sic Lagrange interpolation, of order ordCC , from mid-patch values to the April 8, 2020
Algorithm 4
Evolve in time on patches via function patchSmooth1 whichinterfaces between a user’s macroInt and fun . Function patchSmooth1 calcu-lates the time-derivative/time-step in the interior of all patches via function patchEdgeInt1 which couples the patches with a macro-scale interpolation. function patchSmooth1 ( t , u )2: u := patchEdgeInt1 ( u ) ;3: dudt interior := fun ( t , u , x ) ; (evaluate time-step inside every patch)4: dudt edges ,j := 0 for j := 1 : nPatch ;5: dudt := reshape ( dudt , nSubP × nPatch , ; (make column vector)6: return dudt end function function patchEdgeInt1 ( u )2: u := reshape ( u , nSubP , nPatch ) ; (reshape so each column is a patch)3: u edges ,j := Lagrange interpolation of mid-patch values u i0 ,j from ordCC / patches to either side of each patch.4: return u end function patch edge-values to generally give a macro-scale that is consistent toerrors O (cid:0) H ordCC (cid:1) [Roberts and Kevrekidis, 2007]; – the special case of ordCC = 0 , as used in 2D for the example of Figure 1,invokes a new spectral interpolation that recent numerical experimentsindicate has consistency errors exponentially small in H ; – odd ordCC creates a scheme with a staggered grid of patches suitable formany wave systems [Cao and Roberts, 2013], staggered in the sense thatmid-patch values for odd/even patches use order ordCC interpolation todetermine edge-patch values of even/odd patches, respectively, and that toolbox of Equation-Free functions 15 typically has macro-scales consistent to errors O (cid:0) H ordCC +1 (cid:1) [Cao andRoberts, 2016]; – and the special case of ordCC = − invokes a new staggered spectralinterpolation that recent numerical experiments indicate has consistencyerrors exponentially small in H .After constructing the patches, a user-specified integrator macroInt (step 2 ofAlgorithm 2) simulates the user-defined fun over the time interval [ T , T ] withinitial condition u . Algorithm 4 overviews the toolbox function patchSmooth1 which interfaces between the patch coupling and the user’s function fun thatcomputes the micro-scale time-derivative/time-step within all the patches.The toolbox currently implements corresponding functionality for prob-lems in 2D space, as Figure 1 shows, via toolbox functions configPatches2 , patchSmooth2 , and patchEdgeInt2 .This patch scheme is an example of so-called computational homogenization [e.g., Geers et al., 2010, Saeb et al., 2016, Geers et al., 2017] and is relatedto numerical homogenization [e.g., Craster, 2015, Owhadi, 2015, Peterseim,2019, Maier and Peterseim, 2019]. The three main distinguishing features of thepatch scheme are: that a user need not perform any analysis of the micro-scalestructures; that computations are done only on a (small) subset of the spatialdomain; and for a wide class of systems the scheme is proved to be accurate to auser specified order [Roberts and Kevrekidis, 2007, Roberts et al., 2014, Cao andRoberts, 2016]. However, in application to micro-scale heterogeneous media—the main interest of computational/numerical homogenization—more research April 8, 2020 needs to be done. Bunder et al. [2017] started exploring the patch scheme inheterogeneous media and established it generally has small macro-scale errorsfor diffusion in random media. Further, they found, analogous to that foundin some numerical homogenization [Peterseim, 2019], that when the patchhalf-width is an integral number of periods of the micro-scale heterogeneity,then the macro-scale predictions are accurate to errors O (cid:0) H ordCC (cid:1) , as before.A recent innovation in the toolbox (by setting parameter patches.EdgyInt=1 )is the capability to couple patches by interpolating next-to-edge values to theedge-patch values, but to the opposite edge. This coupling is the subject of anarticle currently in preparation which will discuss how the coupling usefullypreserves symmetry in many applications of interest, and is of controllablemacro-scale accuracy for micro-scale heterogeneous media. Users need to download the toolbox via GitHub . Place the folder of thistoolbox in a path searched by your Matlab /Octave. The toolbox providesboth a user’s and a developer’s manual: start by looking at the User’s Manual.Many of the main toolbox functions, when invoked without any arguments,will simulate an example. Executing the command
PIG() reproduces thePro-jective Integration example presented in Section 1.2. Similarly, the nonlineardiffusion/lubrication-like example of Section 1.1 is reproduced by executing https://github.com/uoa1184615/EquationFreeGit toolbox of Equation-Free functions 17 configPatches2() . The following Sections 3.1 and 3.2 explain the code forboth of these introductory examples as templates to adapt for other problems.3.1 Invoking Projective Integration in GeneralThis subsection discusses some key factors when constructing a ProjectiveIntegration (PI) simulation. Suppose the slow dynamics of your system occurs at rate/frequency of magni-tude about α ; and the rate of decay of your fast modes are higher than thelower bound β (e.g., if three fast modes decay roughly like e − t , e − t , e − t ,then β ≈ ). The PI must be able to stably project the (damped) fast modes,and so the duration δ of the micro-scale burst must be sufficiently long. Thefast modes decay like e − βδ over the micro-burst, and then grow like β∆ on aprojective step of length ∆ . For stability, the product of these effects must beless than one; that is, β∆e − βδ (cid:46) . Rearranging requires the burst length δ (cid:38) β log | β∆ | (3)Figure 3 plots this lower bound as a function of β∆ by writing δ/∆ (cid:38) β∆ log | β∆ | .Once stability of the fast modes is assured by satisfying (3), the accuracy ofthe macro-scale time-step must be considered. The second and fourth order PIfunctions PIRK2() and
PIRK4() have global errors proportional to ∆ and ∆ April 8, 2020
Figure 3
Minimum value (3) for δ/∆ , the ratio of burst length for the micro-simulator tothe macro-time-step, in order to achieve stability in Projective Integration. The horizontalaxis β∆ is the product of the minimum decay rate for the fast modes and the macro-time-step. − − β ∆ ( δ / ∆ ) m i n respectively, for a projective step of length ∆ . More specifically, for a desiredaccuracy of ε and recalling that the slow dynamics occurs at rate/frequency ofmagnitude about α : with PIRK2 we suggest α∆ (cid:46) (6 ε ) / ; whereas with PIRK4 try α∆ (cid:46) ε / . The accuracy of PIG() is controlled by adjusting whatever theaccuracy parameters are provided to the user-specified macro-integrator.
PIG tutorial
We now discuss details of the simulation of the multiscale, slow-fast, ode s (2)shown in Figure 2. First we code the right-hand side function of the system(often people would phrase it in terms of the small parameter (cid:15) = 1 /β ), toolbox of Equation-Free functions 19 beta = 1e5;dxdt=@(t,x) [ cos(x(1))*sin(x(2))*cos(t)beta*( cos(x(1))-x(2) ) ]; Second, we code micro-scale bursts, here using the standard ode45() . Wechoose a burst length (2 /β ) log β as the fast rate of decay is approximately β .Because we do not know the macro-scale time-step invoked by the adaptivefunction to be specified for macroInt (), so we blithely assume ∆ (cid:46) and thendouble the formula (3) for safety. bT = 2/beta*log(beta); Then define the micro-scale burst from state xb0 at time tb0 to be an integrationby the adaptive ode45 of the coded ode s (2), over a burst-time bT . microBurst = @(tb0, xb0) feval(’ode45’,dxdt,[tb0 tb0+bT],xb0); Third, code functions to convert between macro-scale state u and micro-scalestate ( u , u ) . restrict = @(x) x(1);lift = @(X,xApprox) [X; xApprox(2)]; The (optional) restrict() and lift() functions are detailed in the followingSection 3.1.3. Fourth, invoke
PIG to use
Matlab ’s ode45 on the macro-scaleslow evolution. Integrate the micro-bursts over ≤ t ≤ from the initialcondition x (0) = (1 , . Tspan = [0 6];x0 = [1;0];
April 8, 2020 macroInt=’ode45’;[Ts,Xs,tms,xms] = PIG(macroInt,microBurst,Tspan,x0,restrict,lift);
We pause this example to discuss the available outputs from
PIG() . Betweenzero and five outputs may be requested from
PIG() . Most often you wouldstore the first two output results, via say [Ts,Xs] = PIG(...) . – T , an L -vector of times at which macroInt () produced results. – X , an L × N array of the computed solution: the i th row of X , X(i,:) , isto be the macro-state vector X ( t i ) at time t i = T(i) .However, micro-scale details of the underlying Projective Integration com-putations may be helpful, and so
PIG() provides some optional outputs of themicro-scale bursts, via [Ts,Xs,tms,xms] = PIG(...) – tms , optional, is an (cid:96) -dimensional column vector containing micro-scaletimes with bursts, each burst separated by NaN ; – xms , optional, is an (cid:96) × n array of the corresponding micro-scale states.In some contexts it may be helpful to see directly how Projective Integrationapproximates a reduced slow vector field, via [T,X,tms,xms,svf] = PIG(...) in which – svf is a struct containing the Projective Integration estimates of the slowvector field. – svf.T is a ˆ L -dimensional column vector containing all times at whichthe micro-scale simulation data is extrapolated to form an estimate of d x /dt in macroInt (). toolbox of Equation-Free functions 21 – svf.dX is a ˆ L × N array containing the estimated slow vector field.If macroInt () is, for example, the forward Euler method (or the Runge–Kuttamethod), then ˆ L = L (or ˆ L = 4 L ).Returning to the example one remarkable feature of PIG() is revealed: thisPI simulation, which as mentioned in Section 1.2 uses only . as manyevaluations of (2) as direct simulation with ode45() , is accomplished with arecursive call to ode45() . The standard Matlab integrator is used to bothcompute the micro-scale bursts, and also compute the projective steps. All theusual machinery of adaptive time stepping and error control is used to regulatethe micro-scale simulation, and also to regulate the macro-scale projectivetime-steps.
As described in the penultimate paragraph of Section 2.1, the user may requirean intricate or application-specific process to convert between micro-scalevariables u and macro-scale variables U . Provide these bespoke functions tothe toolbox PI functions with the optional inputs restrict() and lift() .The user-provided function restrict(u) should map a micro-scale state u , at some fixed time, to a lower-dimensional macro-scale state U at the sametime. The reverse effect is accomplished by the function lift(U,uApprox) ,which takes two inputs—as lifting is often a non-unique process. The first input U is the macro-scale state at the (present) time at which a micro-scale state April 8, 2020 is desired, and the second input uApprox is the last micro-scale state outputfrom the micro-simulator. This uApprox is typically a micro-scale state froman earlier time, but nonetheless should be useful in initialising a consistentmicro-scale state.A guiding principle in the restriction and lifting functions is that we cannotanticipate a user’s every need; therefore, the functions are straightforward toedit. In particular, their inputs may readily be expanded as needed.
The 2D example of Section 1.2 may be efficiently and simply simulated bystandard stiff integrators, e.g. ode15s() in Matlab . We here demonstratethe advantage of PI over such integrators as the model dimension increases.Consider linear systems of the form d u dt = A u + b , where (10 + N ) × (10 + N ) matrix A is randomly generated so that it has ten eigenvalues with real partwithin [ − . , . , corresponding to ten slow variables, and N eigenvalues withreal part within [ −
20 000 , −
10 000] , corresponding to N fast variables. Thevector b ∈ R N is randomly generated with variance one.We generate such a system at values of N between and , and ateach one simulate to final time from randomly chosen initial conditionswith each of PIRK4() , PIG() and ode15s() . After each simulation we recordthe time elapsed over each simulation as well as the relative error betweenthe simulation estimates and the exact solution. This procedure is repeateda further eleven times (the toolbox script pirk4mance.m details code for this toolbox of Equation-Free functions 23
Figure 4
Performance of Projective Integration compared to a standard stiff integrator,while scaling the dimension of the numerical model described in Section 3.1.4. Solid linesshow the median, and dotted lines show the 25th and 75th percentiles. The stiff integrator isfast and reasonably accurate at low system dimension, but performs poorly at dimension 60and above. − − − System dimension T i m e t a k e n t o s i m u l a t e ode15sPIRK4PIG − − − System dimension E rr o r a t fin a l t i m e ode15sPIRK4PIG experiment): Figure 4 displays the recorded statistics. It appears that forsystem dimension larger than about , the cost by ode15s in setting up andmanaging the Jacobian is too high. For these larger systems, the ProjectiveIntegration functions appear the better choice.The relative performance of ode15s may be improved by providing it withmore information. By providing both the Jacobian A explicitly to the integrator, April 8, 2020 and specifying the initial conditions on the attracting slow manifold, ode15s becomes competitive with the Projective Integration functions (up to systemdimension ). But knowing and coding both the system Jacobian and slowmanifold is usually too hard in practice.3.2 Patch dynamics tutorialThe code here shows one way to get started with the patch scheme. A user’sscript may have the following three steps (arrows indicate function recursion).1. configPatches22. ode integrator ↔ patchSmooth2 ↔ user’s micro-scale code3. process resultsWe reproduce the simulation (Figure 1) of the lattice spatial discretization (1)of a nonlinear diffusion pde . As a micro-scale discretization of the pde weuse the following function to compute the time derivatives: an array variable u(i,j,:,:) refers to the ( i, j ) th point in each and every patch as the thirdand forth indices index the 2D array of patches. function ut = nonDiffPDE(t,u,x,y)dx = diff(x(1:2)); dy = diff(y(1:2)); % microgrid spacingi = 2:size(u,1)-1; j = 2:size(u,2)-1; % interior patch pointsut = nan(size(u)); % preallocate storageut(i,j,:,:) = diff(u(:,j,:,:).^3,2,1)/dx^2 ...+diff(u(i,:,:,:).^3,2,2)/dy^2;end toolbox of Equation-Free functions 25 To use this micro-scale code on patches we need to establish global data,in struct patches , to characterise the patches. Here we aim to simulate thenonlinear diffusive (1) on a macro-scale × -periodic domain with a 2D array of × patches. Choose spectral interpolation ( ordCC = 0 ) to couple the patches,and set each patch of half-size ratio . (relatively large for visualization), andwith × micro-scale lattice points within each patch. Roberts et al. [2014]established that such a patch scheme is consistent with the discretisation (1),as the patch spacing H decreases, and hence is consistent to the original pde . global patchesnSubP = 5;configPatches2(@nonDiffPDE,[-3 3 -2 2],nan,[9 7],0,0.25,nSubP); The third argument to configPatches2() is intended for boundary conditionson the macro-scale simulation, not yet implemented in the toolbox. The inputsare otherwise 2D analogues of the inputs to configPatches1 in Algorithm 2.For an initial condition of the simulation, here set a perturbed-Gaussian.The reshape functions give the x and y lattice coordinates as two 4D arrays ofsize × × × and × × × respectively. Then auto-replication fills u0 with values from the Gaussian u = e − x − y . x = reshape(patches.x,nSubP,1,[],1);y = reshape(patches.y,1,nSubP,1,[]);u0 = exp(-x.^2-y.^2);u0 = u0.*(0.9+0.1*rand(size(u0))); April 8, 2020
Initiate a plot of the simulation using only the micro-scale values interior tothe patches: set x and y -edges to nan to leave, in plots, gaps between patches. figure(1), clfx = patches.x; y = patches.y;x([1 end],:) = nan; y([1 end],:) = nan; Start by showing the initial conditions—the top-left panel of Figure 1. u = reshape(permute(u0,[1 3 2 4]), [numel(x) numel(y)]);hsurf = surf(x(:),y(:),u’); drawnow
Integrate in time using a standard function, or alternatively via projectiveintegration. [ts,us] = ode15s( @patchSmooth2, [0 4], u0(:));
Animate the computed simulation to finish with the bottom-right picture inFigure 1. Use patchEdgeInt2 to interpolate patch-edge values (even if notdrawn) as it reconstitutes the row orientated output from ode15s into arraysthat reflect the 2D patch structure. for i = 1:length(ts)u = patchEdgeInt2(us(i,:));u = reshape(permute(u,[1 3 2 4]), [numel(x) numel(y)]);set(hsurf,’ZData’, u’);legend([’time = ’ num2str(ts(i),2)])pause(0.1)end toolbox of Equation-Free functions 27
This project is collectively developing a
Matlab /Octave toolbox of equation-free algorithms [Roberts et al., 2020]. The algorithms currently implement auseful functionality, but much more is desirable so the plan is to subsequentlydevelop more capability.
Matlab /Octave appears a good choice for a firstversion since it is widespread, efficient, supports various parallel modes, anddevelopment costs are reasonably low. Further, it is built on blas and lapack so the cache and superscalar cpu are potentially well utilised.Projective Integration and Patch Scheme functions are designed to generallyenable users to invoke ‘equation-free’ algorithms in a wide variety of applications.In addition to simulations analogous to those in previous sections, the toolboxmay empower users in other scenarios such as the following list. In the UserManual [Roberts et al., 2020] each toolbox function and application is presentedin its own section, so here we refer to sections of the User Manual by using thename of the appropriate function. – To projectively integrate in time a multiscale, slow-fast, system of ode swith a simple
PI macro-integrator, you could use
PIRK2() , or
PIRK4() for higher-order accuracy. Perhaps adapt the Michaelis–Menten examplepresented in the User Manual at the beginning of
PIRK2.m . – One may use short forward bursts of micro-scale simulation in order tostably predict the slow dynamics backward in time using PI [Gear andKevrekidis, 2003a], as in egPIMM.m . April 8, 2020 – The lifting and restriction functions may be utilised in more complicatedapplications like those others have previously addressed [e.g., Frederix et al.,2007, Roose et al., 2009, Bold et al., 2012, Sieber et al., 2018]. The liftingfunction has two inputs: the macro-scale state (after a macro-scale timestep); and the micro-scale state at the end of the last micro-simulation. Fulldetails for constructing these functions are in
PIG.m . – For the patch scheme in 1D adapt the code at the beginning of configPatches1.m for Burgers’ pde , or the staggered patches of 1D waterwave equations in waterWaveExample.m . – For the patch scheme in 2D adapt the code at the beginning of configPatches2.m for nonlinear diffusion, or the regular patches of the 2Dwave equation of wave2D.m . – The previous two examples are for systems that have smooth spatial struc-tures on the micro-scale: when the micro-scale is ‘rough’ with a knownheterogeneous period (so far only in 1D), then adapt the example of
HomogenisationExample.m . – Employ an ensemble of patch dynamics simulations, averaging appropriately,by adapting the example of ensembleAverageExample.m . – Combine the projective integration and patch functions to simulate a systemwith both time and spatial scale separation. In this case use a PI functionas the integrator for a patch scheme, as done in the second portion of ensembleAverageExample.m . toolbox of Equation-Free functions 29 We encourage adaptation and further development of the toolbox algorithms,and are keen to include new collaborators in future versions of the toolbox.In particular, as well as developing multi-D functions further, corresponding1D coded functionality, we need to code and prove capability for generalmacro-scale boundaries. We also plan to develop projective integration foroscillatory/stochastic systems, perhaps via Dynamic Mode Decomposition [e.g.,Kutz et al., 2016, 2018].
Acknowledgements
This project is supported by the Australian Research Coun-cil via grants DP150102385 and DP180100050.
References
Katherine A. Bold, Karthikeyan Rajendran, Balazs Rath, and Ioannis G.Kevrekidis. An equation-free approach to coarse-graining the dynamics ofnetworks. Technical report, [ http://arxiv.org/abs/1202.5618v1 ], 2012.J. E. Bunder, A. J. Roberts, and I. G. Kevrekidis. Good coupling for themultiscale patch scheme on systems with microscale heterogeneity.
J. Com-putational Physics , 337:154–174, may 2017. doi:10.1016/j.jcp.2017.02.004.Meng Cao and A. J. Roberts. Multiscale modelling couples patches of wave-like simulations. In Scott McCue, Tim Moroney, Dann Mallet, and JudithBunder, editors,
Proceedings of the 16th Biennial Computational Techniquesand Applications Conference, CTAC-2012 , volume 54 of
ANZIAM J. , pagesC153–C170, May 2013. doi:10.21914/anziamj.v54i0.6137.
April 8, 2020
Meng Cao and A. J. Roberts. Multiscale modelling couples patches of non-linear wave-like simulations.
IMA J. Applied Maths. , 81(2):228–254, 2016.doi:10.1093/imamat/hxv034.Claire Y. Chuang, Sang M. Han, Luis A. Zepeda-Ruiz, and Talid Sinno. Oncoarse projective integration for atomic deposition in amorphous systems.
The Journal of Chemical Physics , 143(13):134703, 2015. ISSN 0021-9606.doi:10.1063/1.4931991. URL https://aip.scitation.org/doi/full/10.1063/1.4931991 .J. Cisternas, C. W. Gear, S. Levin, and I. G. Kevrekidis. Equation-free mod-eling of evolving diseases: Coarse-grained computations with individual-based models.
Proc. R. Soc. Lond. A , 460:2761–2779, October 2004.doi:10.1098/rspa.2004.1300 1471-2946.Richard V. Craster.
Dynamic Homogenization , volume 116 of
SpringerProceedings in Mathematics and Statistics , pages 41–50. Springer, 2015.doi:10.1007/978-3-319-12148-2_3.Weinan E. Analysis of the heterogeneous multiscale method for ordinarydifferential equations.
Communications in Mathematical Sciences , 1(3):423–436, 2003. doi:10.4310/CMS.2003.v1.n3.a3.Radek Erban, Ioannis G. Kevrekidis, and Hans G. Othmer. Anequation-free computational approach for extracting population-level behavior from individual-based models of biological dis-persal.
Physica D: Nonlinear Phenomena , 215(1):1–24, 2006.ISSN 0167-2789. doi:10.1016/j.physd.2006.01.008. URL http: toolbox of Equation-Free functions 31 .Yves Frederix, Giovanni Samaey, Christophe Vandekerckhove, and DirkRoose. Equation-free methods for molecular dynamics: a liftingprocedure.
Proc. Appl. Meth. Mech. , 7:20100003–20100004, 2007.doi:10.1002/pamm.200700025.C. W. Gear and I. G. Kevrekidis. Computing in the past with forward integra-tion.
Phys. Lett. A , 321:335–343, 2003a. doi:10.1016/j.physleta.2003.12.041.C. W. Gear and Ioannis G. Kevrekidis. Projective methods for stiffdifferential equations: Problems with gaps in their eigenvalue spec-trum.
SIAM Journal on Scientific Computing , 24(4):1091–1106, 2003b.doi:10.1137/S1064827501388157. URL http://link.aip.org/link/?SCE/24/1091/1 .C. W. Gear, T. J. Kaper, I. G. Kevrekidis, and A. Zagaris. Projecting to a slowmanifold: singularly perturbed systems and legacy codes.
SIAM J. AppliedDynamical Systems , 4(3):711–732, 2005. doi:10.1137/040608295.C William Gear and Ioannis G Kevrekidis. Constraint-defined manifolds: alegacy code approach to low-dimensional computation.
Journal of ScientificComputing , 25(1):17–28, 2005. doi:10.1007/s10915-004-4630-x.M. G. D. Geers, V. G. Kouznetsova, and W. A. M. Brekelmans. Multi-scale computational homogenization: Trends and challenges.
Journal ofComputational and Applied Mathematics , 234(7):2175–2182, August 2010.ISSN 0377-0427. doi:10.1016/j.cam.2009.08.077.
April 8, 2020
Marc G. D. Geers, Varvara G. Kouznetsova, Karel Matouš, and JulienYvonnet. Homogenization methods and multiscale modeling: Non-linear problems. In
Encyclopedia of Computational Mechanics, Sec-ond Edition , pages 1–34. Wiley, 2017. ISBN 978-1-119-17681-7.doi:10.1002/9781119176817.ecm2107. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119176817.ecm2107 .Dror Givon, Ioannis G. Kevrekidis, and Raz Kupferman. Strong con-vergence of projective integration schemes for singularly perturbedstochastic differential systems.
Comm. Math. Sci. , 4(4):707–729, 2006.doi:10.4310/CMS.2006.v4.n4.a2.I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg, andK. Theodoropoulos. Equation-free, coarse-grained multiscale computation:enabling microscopic simulators to perform system level tasks.
Comm. Math.Sciences , 1:715–762, 2003. doi:10.4310/CMS.2003.v1.n4.a5.I. G. Kevrekidis, C. W. Gear, and G. Hummer. Equation-free: the computer-assisted analysis of complex, multiscale systems.
A. I. Ch. E. Journal , 50:1346–1354, 2004. doi:10.1002/aic.10106.Ioannis G. Kevrekidis and Giovanni Samaey. Equation-free multiscale compu-tation: Algorithms and applications.
Annu. Rev. Phys. Chem. , 60:321—44,2009. doi:10.1146/annurev.physchem.59.032607.093610.J. Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua L. Proctor.
Dynamic Mode Decomposition: Data-driven Modeling of Complex Systems .Number 149 in Other titles in applied mathematics. SIAM, Philadelphia, toolbox of Equation-Free functions 33
Complexity , (6010634):1–16, 2018. doi:10.1155/2018/6010634.Steven L Lee and C William Gear. Second-order accurate projective integratorsfor multiscale problems.
Journal of Computational and Applied Mathematics ,201(1):258–274, 2007.John Maclean. A note on implementations of the boosting algorithm andheterogeneous multiscale methods.
SIAM Journal on Numerical Analysis ,53(5):2472–2487, 2015. doi:10.1137/140982374.John Maclean and Georg A. Gottwald. On convergence of higher order schemesfor the projective integration method for stiff ordinary differential equa-tions.
Journal of Computational and Applied Mathematics , 288:44–69, 2015.doi:http://dx.doi.org/10.1016/j.cam.2015.04.004.Roland Maier and Daniel Peterseim. Explicit computational wave propagationin micro-heterogeneous media.
BIT Numerical Mathematics , 59(2):443–462,June 2019. ISSN 1572-9125. doi:10.1007/s10543-018-0735-8.Houman Owhadi. Bayesian numerical homogenization.
Multiscale Modeling& Simulation , 13(3):812–828, 2015. doi:10.1137/140974596. URL http://dx.doi.org/10.1137/140974596 .Daniel Peterseim. Numerical homogenization beyond scale separation and peri-odicity. Technical report, AMSI Winter School on Computational Modelingof Heterogeneous Media, https://ws.amsi.org.au/wp-content/uploads/
April 8, 2020 sites/70/2019/06/numhomamsi2019.pdf , 6 2019.R. Rico-Martinez, C. W. Gear, and Ioannis G. Kevrekidis. Coarseprojective kmc integration: forward/reverse initial and bound-ary value problems.
Journal of Computational Physics , 196(2):474–489, 2004. doi:10.1016/j.jcp.2003.11.005. URL .A. J. Roberts.
Model emergent dynamics in complex systems . SIAM, Philadel-phia, jan 2015. ISBN 9781611973556. URL http://bookstore.siam.org/mm20/ .A. J. Roberts and I. G. Kevrekidis. General tooth boundary conditions forequation free modelling.
SIAM J. Scientific Computing , 29(4):1495–1510,2007. doi:10.1137/060654554.A. J. Roberts, Tony MacKenzie, and Judith E. Bunder. A dynamical systemsapproach to simulating macroscale spatial dynamics in multiple dimensions.
J. Engineering Mathematics , 86(1):175–207, 2014. doi:10.1007/s10665-013-9653-6.A. J. Roberts, John Maclean, and J. E. Bunder. Equation-free function toolboxfor matlab/octave. Technical report, [ https://github.com/uoa1184615/EquationFreeGit ], 2020.Dirk Roose, Erik Nies, Ting Li, Christophe Vandekerckhove, Giovanni Samaey,and Yves Frederix. Lifting in equation-free methods for molecular dynamicssimulations of dense fluids.
Discrete and Continuous Dynamical Systems— toolbox of Equation-Free functions 35
Series B , 11(4):855–874, April 2009. doi:10.3934/dcdsb.2009.11.855.Saba Saeb, Paul Steinmann, and Ali Javili. Aspects of Computational Homog-enization at Finite Deformations: A Unifying Review From Reuss’ to Voigt’sBound.
Applied Mechanics Reviews , 68(5), September 2016. ISSN 0003-6900.doi:10.1115/1.4034024.G. Samaey, A. J. Roberts, and I. G. Kevrekidis. Equation-free computation:an overview of patch dynamics. In Jacob Fish, editor,
Multiscale methods:bridging the scales in science and engineering , chapter 8, pages 216–246.Oxford University Press, 2010. ISBN 978-0-19-923385-4.S. Setayeshgar, C. W. Gear, H. G. Othmer, and I. G. Kevrekidis. Application ofcoarse integration to bacterial chemotaxis.
SIAM J. Mathematical Modelingand Simulation , 4:307–327, 2005. http://epubs.siam.org/sam-bin/dbq/article/60087 .J. Sieber, C. Marschler, and J. Starke. Convergence of Equation-Free Methodsin the Case of Finite Time Scale Separation with Application to Deterministicand Stochastic Systems.
SIAM Journal on Applied Dynamical Systems , 17(4):2574–2614, January 2018. doi:10.1137/17M1126084.C. I. Siettos, M. D. Graham, and I. G. Kevrekidis. Coarse brownian dynamicsfor nematic liquid crystals: Bifurcation, projective integration, and controlvia stochastic simulation.
J. Chemical Physics , 118(22):10149–10156, 2003.doi:10.1063/1.1572456., 118(22):10149–10156, 2003.doi:10.1063/1.1572456.