Equation-free modeling explained

Equation-free modeling is a method for multiscale computation and computer-aided analysis. It is designed for a class of complicated systems in which one observes evolution at a macroscopic, coarse scale of interest, while accurate models are only given at a finely detailed, microscopic, level of description. The framework empowers one to perform macroscopic computational tasks (over large space-time scales) using only appropriately initialized microscopic simulation on short time and small length scales. The methodology eliminates the derivation of explicit macroscopic evolution equations when these equations conceptually exist but are not available in closed form; hence the term equation-free.

Introduction

In a wide range of chemical, physical and biological systems, coherent macroscopic behavior emerges from interactions between microscopic entities themselves (molecules, cells, grains, animals in a population, agents) and with their environment. Sometimes, remarkably, a coarse-scale differential equation model (such as the Navier-Stokes equations for fluid flow, or a reaction–diffusion system) can accurately describe macroscopic behavior. Such macroscale modeling makes use of general principles of conservation (atoms, particles, mass, momentum, energy), and closed into a well-posed system through phenomenological constitutive equations or equations of state. However, one increasingly encounters complex systems that only have known microscopic, fine scale, models. In such cases, although we observe the emergence of coarse-scale, macroscopic behavior, modeling it through explicit closure relations may be impossible or impractical. Non-Newtonian fluid flow, chemotaxis, porous media transport, epidemiology, brain modeling and neuronal systems are some typical examples. Equation-free modeling aims to use such microscale models to predict coarse macroscale emergent phenomena.

Performing coarse-scale computational tasks directly with fine-scale models is often infeasible: direct simulation over the full space-time domain of interest is often computationally prohibitive. Moreover, modeling tasks, such as numerical bifurcation analysis, are often impossible to perform on the fine-scale model directly: a coarse-scale steady state may not imply a steady state for the fine-scale system, since individual molecules or particles do not stop moving when the gas density or pressure become stationary. Equation-free modeling circumvents such problems by using short bursts of appropriately initialized fine-scale simulation, and, in spatial problems, on small well-separated patches of space. A free Matlab/Octave toolbox empowers people to use these equation-free methods.

The coarse time-stepper

Dynamic problems invoke the coarse time-stepper. In essence, short bursts of computational experiments with the fine-scale simulator estimate local time derivatives. Given an initial condition for the coarse variables

U(t _{k)} at time t _{k} , the coarse time-stepper involves four steps:

• Lifting, creates microscale initial conditions

u(t _{k)} , consistent with the macrostate U(t _{k)} ;

• Simulation, uses the microscale simulator to compute the microscale state

u(t) at over a short interval t _{k\leq} t \leq t _{k+\delta} t ;

• Restriction, obtains the macrostate

U(t _{k+\delta} t) from the fine-scale state u(t) ;

• the Time-step, extrapolation of macrostate

U from t _{k} to t _{k+1} = t _{k+\Delta} t predicts the state a macrotime in the future.

Multiple time steps simulates the system into the macro-future.If the microscale model is stochastic, then an ensemble of microscale simulations may be needed to obtain sufficiently good extrapolation in the time step. Such a coarse time-stepper may be used in many algorithms of traditional continuum numerical analysis, such as numerical bifurcation analysis, optimization, control, and even accelerated coarse-scale simulation.For deterministic systems, the Matlab/Octave toolbox provides a user with higher-order accurate time-steppers: a second- and fourth-order Runge--Kutta scheme, and a general interface scheme.

Traditionally, algebraic formulae determine time derivatives of the coarse model. In this approach, the macroscale derivative is estimated by the inner microscale simulator, in effect performing a closure on demand. A reason for the name equation-free is by analogy with matrix-free numerical linear algebra;^[1] the name emphasizes that macro-level equations are never constructed explicitly in closed form.

Restriction

The restriction operator often follows directly from the specific choice of the macroscale variables. For example, when the microscale model evolves an ensemble of many particles, the restriction typically computes the first few moments of the particle distribution (the density, momentum, and energy).

Lifting

The lifting operator is usually much more involved. For example, consider a particle model: we need to define a mapping from a few low order moments of the particle distribution to initial conditions for each particle. The assumption that a relation exists that closes in these low order, coarse, moments, implies that the detailed microscale configurations are functionals of the moments (sometimes referred to as slaving ^[2]). We assume this relationship is established/emerges on time scales that are fast compared to the overall system evolution (see slow manifold theory and applications ^[3]). Unfortunately, the closure (slaving relations) are algebraically unknown (as otherwise the coarse evolution law would be known).

Initializing the unknown microscale modes randomly introduces a lifting error: we rely on the separation of macro and micro time scales to ensure a quick relaxation to functionals of the coarse macrostates (healing). A preparatory step may be required, possibly involving microscale simulations constrained to keep the macrostates fixed.^[4] When the system has a unique fixed point for the unknown microscale details conditioned upon the coarse macrostates, a constrained runs algorithm may perform this preparatory step using only the microscale time-stepper.^[5]

An illustrative example

A toy problem illustrates the basic concepts. For example, consider the differential equation system for two variables

(U(t),u(t)) :(dU)/(dt) =- U - u +2 , (du)/(dt) =100 (U ^{3-u).} Capital U denotes the presumed macroscale variable, and lowercase u the microscale variable. This classification means that we assume a coarse model of the form dU/dt = G(U) exists, although we do not necessarily know what it is. Arbitrarily define the lifting from any given macrostate U as (U,u) = (U, 0 . 5 ) . A simulation using this lifting and the coarse time-stepper is shown in the figure.u ≈ U ^{3} for any initial data. The coarse time-stepper solution would agree better with the full solution when the 100 factor is increased. The graph shows the lifted solution (blue solid line) (U,u) . At times t = n\delta t , the solution is restricted and then lifted again, which here is simply setting u(n\delta t) =0 . 5 . The slow manifold is shown as a red line. The right plot shows the time derivative of the restricted solution as a function of time (blue curve), as well as the time derivative dU/dt (the coarse time derivative), as observed from a full simulation (red curve).

On application to concrete multiscale problems

The equation-free approach has been applied to many examples. The examples illustrate the various ways to construct and assemble the algorithmic building blocks. Numerical analysis establishes the accuracy and efficiency of this approach. Additional numerical analysis on other methods of this type has also been done.^[6]

Applying the equation-free paradigm to a real problem requires considerable care, especially defining the lifting and restriction operators, and the appropriate outer solver.

• The first challenge is to identify the macroscale observables. They need to be complete enough so that the unknown microscale variables can be reliably reconstructed (lifted). Physical arguments often identify the macroscale observables. Almost always one invokes densities, but there are some surprisingly simple examples where correlation functions are essential macroscale variables.^[7] If not resorting to physical arguments, then modern data-mining or manifold learning techniques, such as Isomap or diffusion maps, may obtain the macroscale variables from microscale simulation.^[8]
• There should be a clear separation between the time-scales of the macroscale observables and the time-scales of the remaining microscale modes quasi-equilibrate given any macrostate.
• Knowing the macroscale observables might not be sufficient. One strategy to obtain such information is the baby-bathwater scheme which uses only appropriately initialized simulations.^[9]

Coarse bifurcation analysis

The recursive projection method^[10] enables the computation of bifurcation diagrams using legacy simulation code. It also empowers the coarse time-stepper to perform equation-free bifurcation computations. Consider the coarse time stepper in its effective form

\vec U ^{n+1} = \vec S(\vec U ^{n,λ;} \delta t), which includes explicit dependence upon one or more parameters λ . Bifurcation analysis computes equilibria or periodic orbits, their stability and dependence upon parameter λ .

Compute a coarse equilibrium as a fixed point of the coarse time stepper

\vec U - \vec S(\vec U, λ ; \delta t) = \vec 0 . In the equation-free context, the recursive projection method is the outer solver of this equation, and the coarse time-stepper enables this method to be performed using fine scale dynamics.

Additionally, for problems where the macroscale has continuous symmetries, one can use a template based approach ^[11] to compute coarse self-similar or travelling wave solutions as fixed points of a coarse time-stepper that also encodes appropriate rescaling and/or shifting of space-time and/or solution.For example, self-similar diffusion solutions may be found as the probability density function of detailed molecular dynamics.^[12]

An alternative to the recursive projection method is to use Newton—Krylov methods.^[13]

Coarse projective integration

The coarse time stepper accelerates simulation over large macroscale times. In the scheme described above, let the large macro-time-step

\Delta t\gg\delta t , and be on the time scale of the slow coarse dynamics. Let the computed U ^{n ≈} U(n\Delta t) in terms of the coarse variable, and let the microscale simulation compute U ^{n,k} ≈ U(n\Delta t + k\delta t) from a local time simulation with initial condition that the coarse variable U(n\Delta t) = U ^{n=U} ^{n,0} . Then we approximate U((n +1 )\Delta t) via extrapolating over a gap byU ^{n+1} = U ^{n,k} + (\Delta t - k\delta t)F(\vec U ^{n)} where, for example, simple linear extrapolation would be F(\vec U ^{n)=(U} ^{n,k} - U ^{n,k-1} )/\delta t This scheme is the called coarse projective forward Euler, and is the simplest in the class.

The

k steps taken before the extrapolation reflect that we must allow the system to settle onto a quasi-equilibrium (from the microscale point of view), so that we can make a reliable extrapolation of the slow dynamics. Then the size of the projective integration step is limited by stability of the slow modes.^[14]

Higher order versions of coarse projective integration can be formed, analogous to Adams–Bashforth or Runge–Kutta.^[15] Higher order schemes for systems where the microscale noise is still apparent on the macroscale time step are more problematic.^[16]

Patch dynamics

The spatial analogue of projective integration is the gap-tooth scheme.The idea of the gap-tooth scheme is to perform simulations of small patches of space, the teeth, separated by unsimulated space, the gaps.By appropriately coupling the small patches of simulations we create a large scale, coarse level, simulation of the spatially extended system.When the microscale simulator is computationally expensive the gap-tooth scheme empowers efficient large scale prediction.Furthermore, it does this without us ever having to identify an algebraic closure for a large scale model.^{[17] [18] [19]} The Matlab/Octave toolbox provides support to users to implement simulations on an rectangular grid of patches in 1D or 2D space.

The combination of the gap-tooth scheme with coarse projective integration is called patch dynamics.

Coupling boundary conditions

The key to the gap-tooth and patch scheme is the coupling of the small patches across unsimulated space.Surprisingly, the generic answer is to simply use classic Lagrange interpolation, whether in one dimension^[19] or multiple dimensions.^[20] This answer is related to the coupling in holistic discretization and theoretical support provided by the theory of slow manifolds.The interpolation provides value or flux boundary conditions as required by the microscale simulator. High order consistency between the macroscale gap-tooth/patch scheme and the microscale simulation is achieved through high order Lagrange interpolation.

However, commonly the microscale is a noisy particle based or agent-based model.In such cases the relevant macroscale variables are averages such as mass and momentum density. Then one generally has to form averages over a core of each tooth/patch, and apply the coupling condition over a finite action region on the edges of each tooth/patch.The provisional recommendation is to make these regions as big as half the tooth/patch.^[21] That is, for efficiency one makes the microscale tooth/patch as small as possible, but limited by the need to fit in action and core regions big enough to form accurate enough averages.

Lifting

Patch dynamics is the combination of the gap-tooth scheme and coarse projective integration. Just as for normal projective integration, at the start of each burst of microscale simulation, one has to create an initial condition for each patch that is consistent with the local macroscale variables, and the macroscale gradients from neighboring interpolated patches. The same techniques suffice.

Open problems and future directions

Assumptions and choices about the macroscale evolution are crucial in the equation-free scheme. The key assumption is that the variables we choose for the macroscale coupling must effectively close on the chosen macroscale. If the chosen macroscale length is too small then more coarse scale variables may be needed: for example, in fluid dynamics we conventionally close the PDEs for density, momentum and energy; yet in high speed flow especially at lower densities we need to resolve modes of molecular vibration because they have not equilibrated on the time scales of the fluid flow. Qualitatively the same considerations apply to the equation-free approach.

For many systems appropriate coarse variables are more-or-less known by experience. However, in complex situations there is a need to automatically detect the appropriate coarse variables, and then use them in the macroscale evolution. This needs much more research utilizing techniques from data mining and manifold learning. In some problems it could be that as well as densities, the appropriate coarse variables also need to include spatial correlations, as in the so-called Brownian bugs.^[22]

The macroscale may have to be treated as a stochastic system, but then the errors are likely to be much larger and the closures more uncertain.

Notes and References

1. C. T. Kelley. Iterative Methods for linear and nonlinear equations, SIAM, Philadelphia, 1995.
2. H. Haken. Slaving principle revisited. Physica D, 97:95–103, 1996.
3. A. J. Roberts. Effectively model dynamics, deterministic and stochastic, across multiple space and time scales. In J. G. Hartnett and P. C. Abbott, editors, Frontiers of Fundamental and Computational Physics: 10th International Symposium, volume 1246, pages 75–87. AIP, 2010.
4. J. P. Ryckaert, G. Ciccotti, and H. Berendsen. Numerical integration of the Cartesian equation of motion of a system with constraints: molecular dynamics of N-alkanes. J. Comput. Phys., 23:237, 1977.
5. C. W. Gear, T. J. Kaper, I. G. Kevrekidis, and A. Zagaris. Projecting to a Slow Manifold: Singularly Perturbed Systems and Legacy Codes. SIAM Journal on Applied Dynamical Systems 4(3):711–732, 2005.
6. W. E and B. Engquist (2003). The heterogeneous multiscale methods Comm. Math. Sciences 1(1):87–132.
7. W. R. Young, A. J. Roberts, and G. Stuhne. Reproductive pair correlations and the clustering of organisms. Nature, 412:328–331, 2001.
8. R. R. Coifman et al. (2005). Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps Proceedings of the National Academy of Sciences 102(21):7426–7431.
9. J. Li, P. G. Kevrekidis, C. W. Gear and I. G. Kevrekidis (2003). Deciding the nature of the coarse equation through microscopic simulations: the baby-bathwater scheme SIAM Multiscale Modeling and Simulation 1(3):391–407.
10. G.M. Schroff and H.B. Keller (1993). Stabilization of unstable procedures: the recursive projection method SIAM Journal on Numerical Analysis 30: 1099–1120.
11. C. Rowley and J. Marsden (2000). Reconstruction equations and the Karhunen–Loève expansion for systems with symmetry Physica D: Nonlinear Phenomena 142: 1–19.
12. L. Chen, P. Debenedetti, C.W. Gear, and I.G. Kevrekidis (2004). From molecular dynamics to coarse self-similar solutions: a simple example using equation-free computation Journal of Non-Newtonian Fluid Mechanics 120: 215–223.
13. C.T. Kelley (1995). Iterative Methods for linear and nonlinear equations SIAM, Philadelphia.
14. C.W. Gear and I.G. Kevrekidis. Projective methods for stiff differential equations: problems with gaps in their eigenvalue spectrum . SIAM Journal on Scientific Computing 24(4):1091–1106, 2003.
15. C.W. Gear; I.G. Kevrekidis and Theodoropoulos. Coarse integration/bifurcation analysis via microscopic simulators: micro-Galerkin methods Computers and Chemical Engineering 26: 941–963, 2002.
16. X. Chen, A. J. Roberts, and I. G. Kevrekidis. Projective integration of expensive multiscale stochastic simulation. In W. McLean and A. J. Roberts, editors, Proceedings of the 15th Biennial Computational Techniques and Applications Conference, CTAC-2010, volume 52 of ANZIAM J., pages C661–C677, Aug. 2011. http://journal.austms.org.au/ojs/ index.php/ANZIAMJ/article/view/3764
17. Kevrekidis, I.G. et al. (2003). Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level tasks Comm. Math. Sciences 1(4): 715–762.
18. Samaey, G.; Roose, D. and Kevrekidis, I.G. (2005). The gap-tooth scheme for homogenization problems SIAM Multiscale Modeling and Simulation 4: 278–306.
19. Roberts, A.J. and Kevrekidis, I.G. (2007). General tooth boundary conditions for equation free modelling SIAM J. Scientific Computing 29(4): 1495–1510.
20. A. J. Roberts, T. MacKenzie, and J. Bunder. A dynamical systems approach to simulating macroscale spatial dynamics in multiple dimensions. J. Engineering Mathematics, 86(1):175–207, 2014.
21. Bunder, J. E., A. J. Roberts, and I. G. Kevrekidis (2017). “Good coupling for the multiscale patch scheme on systems with microscale heterogeneity”. In: J. Computational Physics 337, pp. 154–174. http://dx.doi.org/10.1016/j.jcp.2017.02.004
22. W. R. Young, A. J. Roberts, and G. Stuhne. Reproductive pair correlations and the clustering of organisms. Nature, 412:328–331, 2001.