Equation free modeling
Equation-free modeling is a collective term referring to a paradigm for multiscale computation and computer-aided analysis. It is designed for a class of complex/multiscale problems in which one observes evolution at a macroscopic, coarse scale of interest, while accurate models are only given at a more detailed (fine-scale, microscopic, possibly atomistic) level of description. The framework enables one to perform macroscopic computational tasks (over extended spatiotemporal scales) using only appropriately initialized microscopic simulation on short time and length scales. The methodology bypasses the derivation of explicit macroscopic evolution equations when these equations conceptually exist but are not available in closed form – hence the term equation-free.
In a wide range of chemical, physical and biological systems, macroscopic, coherent behavior emerges from interactions between microscopic entities (molecules, cells, individuals in a population) among themselves and with their environment. Sometimes, remarkably, a coarse-scale model (such as the Navier-Stokes equations for fluid flow, or a reaction-diffusion system) can accurately describe behavior at this level, making use of general principles of conservation (species, mass, momentum, energy) closed through phenomenological constitutive equations (reaction rates as functions of concentrations, viscous stresses as functionals of velocity gradients). However, one increasingly encounters complex systems that can only be modeled with sufficient accuracy at a microscopic, fine scale. In such cases, although one observes the emergence of coarse-scale, macroscopic behavior in practice, modeling it through explicit closure relations may be impossible or impractical without simplifying assumptions that are hard to justify. Non-Newtonian fluid flow, chemotaxis, porous media transport, and neuronal systems are but a few typical examples in which there is a trend of moving away from empirical closures towards more fine-scale models.
Performing coarse-scale computational tasks with fine-scale models is often infeasible: direct simulation over the full spatiotemporal domain of interest can be computationally prohibitive. Moreover, additional modeling tasks, such as numerical bifurcation analysis, are often impossible to perform on the fine-scale model directly: a coarse steady state may not imply a steady state for the fine-scale system, since (for instance) individual molecules do not stop moving when the gas density or pressure become stationary.
The goal of equation-free modeling is to circumvent the explicit derivation of coarse equations by using short bursts of appropriately initialized fine-scale simulation. Thus, a (deterministic or stochastic) fine-scale model <math>f</math> and an associated time-stepper <math>s</math> with time step <math>\textrm{d} t\ ,</math>
- <math>\label{fine}
are given in terms of fine-scale variables <math>u(t)\ ,</math> while a coarse model <math>F</math> and its associated time-stepper <math>S</math> with time-step <math>\delta t\ ,</math>
- <math>\label{coarse}
in terms of some (known) coarse variables are assumed to exist, but are unavailable in closed form. The article mainly treats the case when this coarse-scale model is deterministic. Some comments on stochastic models at the coarse scale are given towards the end.
The idea was first introduced in (Theodoropoulos et al, 2000) and (Gear et al., 2002); see (Kevrekidis et al., 2003) for an early review. A more recent review, on which this entry is based, is (Kevrekidis and Samaey, 2009). The reviews also contain ample references on how the equation-free framework relates to, and borrows ideas from, a long history of (computational) multiscale approaches. (See also "Further reading" for a limited list of references).
A key tool is the coarse time-stepper: a short computational experiment with the fine-scale simulator. Given an initial condition for the coarse variables <math>U(t^*)</math> at <math>t^*\ ,</math> the coarse time-stepper involves:
- Lifting. Create fine-scale initial conditions <math>u(t^*)\ ,</math> consistent with <math>U(t^*)\ ;</math>
- Simulation. Use the fine-scale simulator \eqref{fine} to compute the fine-scale state <math>u(t)</math> at <math>t \in [t^*,t^*+\delta]\ ;</math>
- Restriction. Obtain the coarse state <math>U(t^*+\delta t)</math> from the fine-scale state <math>u(t)\ .</math>
- <math>\label{coarse-ts}
where the superscript on <math>s</math> denotes the <math>k</math> fine-scale time steps <math>\textrm{d} t\ ,</math> such that <math>\delta t = k \textrm{d} t\ ,</math> <math>\mu</math> is the lifting operator, <math>\mathcal{M}</math> is the restriction operator, and the overbars are introduced to emphasize that these are (time-discrete) approximations. If the fine-scale model is stochastic, multiple replica simulations, using an ensemble of fine-scale initial conditions, may be needed to obtain sufficiently low-variance results.
Once a coarse time-stepper is available, a direct bridge can be built between fine-scale simulation and algorithms of traditional continuum numerical analysis, such as numerical bifurcation analysis, optimization, control, and even accelerated coarse-scale simulation. Traditionally, a coarse-level solver (here called outer solver) obtains the numerical quantities it requires for computation (time derivatives, the action of Jacobians) using explicit formulae from the coarse model. In our approach, the outer solver designs, executes, and processes the results of computational experiments with an inner fine-scale simulator, in effect performing a closure on demand. The analogy with matrix-free numerical linear algebra (Kelley, 1995) provides another reason for the name equation-free; it emphasizes that the coarse-level equations are never constructed explicitly in closed form.
The restriction operator can often be determined as soon as the coarse variables are chosen. For instance, when the fine-scale model evolves an ensemble of many particles, the restriction typically computes the first few moments of the particle distribution (density, momentum, energy), which can be represented by evaluation on a coarse grid, in a finite element or spectral basis, or even using empirical basis functions (EOFs, KL or POD modes).
The construction of the lifting operator is usually much more involved. Again considering a particle model example, 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 an equation exists that closes at the level of these low order moments, implies that higher order moments become functionals of (slaved to) the low order ones on time scales that are fast compared to the overall system evolution. (If this were not true, a description in terms of the chosen low order moments would contain memory terms, as can be seen from the Mori-Zwanzig formalism (Mori, 1965) and (Zwanzig, 1973).) The same principle also underpins quasi-steady state approximations. Unfortunately, the slaving relations (closures) are unknown (otherwise the coarse evolution law would be explicitly known).
Initializing the unknown degrees of freedom randomly introduces a lifting error, and one then relies on separation of time scales to ensure their quick relaxation to functionals of the coarse observables (healing). A preparatory step, possibly involving fine-scale simulations constrained to keep the coarse observables fixed, may be required, e.g. (Ryckaert et al., 1977). When the system has a unique fixed point for the unknown degrees of freedom conditioned upon the coarse observables, it is possible to construct a constrained runs algorithm to perform this preparatory step using only the fine-scale time-stepper, see e.g. (Gear et al.,2005).
The basic concepts are most easily illustrated on a toy problem: a singularly perturbed system of ordinary differential equations. Consider for example
- <math>\label{eq-example}
as the fine-scale equation. As the coarse variable, we consider <math>U=\mathcal{M}(u_1,u_2)=u_1\ ,</math> implying that we assume a coarse model of the form <math>\dot{U}=G(U)</math> to exist. We define the lifting as <math>u=\mu(U) = (U,1)\ .</math> A simulation using the coarse time-stepper is shown in <figref>Equation-free_modeling_example.png</figref>.
The solution of equation \eqref{eq-example} rapidly moves to the slow manifold <math>u_2=u_1^3</math> for any initial data when <math>\varepsilon\ll 1\ .</math> The coarse time-stepper solution agrees better with the full solution (actually, its restriction) when <math>\varepsilon</math> is decreased (for fixed <math>\delta t</math>). The left figure shows the lifted solution (blue solid line) <math>u(n\delta t +s)\ ,</math> <math>s\in [0,\delta]\ ,</math> plotted in the <math>(u_1,u_2)</math> phase space. At times <math>t=n{\delta t}\ ,</math> the solution is restricted and then lifted again, which here amounts to setting <math>u_2(n\delta t)=1\ .</math> The slow manifold is shown as a dashed red line. The right figure shows the time derivative of the restricted solution <math>\mathcal{M}(u(n\delta t+s))</math> as a function of time <math>t=n\delta t+s</math> (blue solid line), as well as the time derivative of <math>u_1</math> (the coarse time derivative), as observed from a full simulation of \eqref{eq-example}. (red dashed line).
The approach has been applied to many concrete examples that illustrate how the algorithmic building blocks can be constructed and assembled. Moreover, some numerical analysis has been performed that discusses accuracy and efficiency of these methods. Extensive additional numerical analysis results on methods of this type can be found in the literature on heterogeneous multiscale methods (E and Engquist, 2003).
However, applying the equation-free paradigm to a concrete problem requires considerable care in the selection of all building blocks. One should judiciously define lifting and restriction operators, as well as the appropriate outer solver. Several remarks are in order:
- One has to know which coarse observables appear in the unavailable coarse equation, or, equivalently, from which coarse observables the unknown fine-scale degrees of freedom can be reliably reconstructed (lifted). If physical arguments are not sufficient to obtain these coarse observables, one might use modern data-mining/manifold learning techniques, such as Isomap or diffusion maps (Coifman et al., 2005) to obtain the coarse variables themselves from fine-scale simulation.
- For the methodology to be effective, there should be a clear separation between the time-scales at which the coarse observables evolve and the time-scales at which the remaining degrees of freedom equilibrate with respect to these observables.
- If one has knowledge on the dependence of the unknown fine-scale degrees of freedom on the coarse observables, this can and should be included when designing a lifting procedure.
- To decide on an appropriate coarse-scale solver (see later), knowing the coarse observables might not be sufficient, and one might need more information on the nature of the coarse model, such as the order or character (parabolic, hyperbolic) of a postulated unavailable partial differential equation. A strategy to obtain such information, also using only appropriately initialized simulations with the fine-scale model, is the baby-bathwater scheme (Li et al., 2003).
Inspired by the recursive projection method (RPM) (Schroff and Keller, 1993), a computational superstructure (wrapper) that enables the computation of bifurcation diagrams using a legacy simulation code, the coarse time-stepper was first used to perform equation-free bifurcation computations. Consider the coarse equation, and its coarse time-stepper
- <math>\label{eq-coarse-bif}
which now include explicit dependence on one or more parameters <math>\lambda\ ,</math> and suppose one wants to compute asymptotic solutions (steady states, periodic orbits), their stability and dependence on <math>\lambda\ .</math>
A coarse steady state can be computed as a fixed point of the coarse time-stepper,
- <math>\label{eq-fixed}
RPM combines Newton-Raphson iterations in a low-dimensional slow subspace with Picard iterations in its orthogonal complement. In the equation-free context, RPM is the outer solver; the coarse time-stepper enables the RPM computation to be performed using a fine-scale inner solver. The coarse time-stepper immediately enables continuation/bifurcation computations at the coarse level, even when parameters of the fine-scale model are varied, whose influence on the coarse level is difficult to assess directly. The procedure is illustrated in <figref>Equation-free_modeling_equation-free-main.png</figref>.
Additionally, for problems with continuous symmetries, one can use a template-based approach (Rowley and Marsden, 2000) to compute coarse self-similar solutions, such as traveling waves/exploding solutions, as fixed points of a coarse time-stepper that incorporates an appropriate shifting/rescaling of space, time and/or the solution. Coarse traveling speeds/similarity exponents are a byproduct of the procedure upon convergence. An illustrative example of outer coarse self-similar diffusion dynamics based on inner molecular dynamics (MD) (Chen et al., 2004), is depicted in <figref>Equation-free_modeling_self-similar.png</figref>.
As alternative to RPM, one can also use a Newton-Krylov method. One can estimate matrix-vector products with the Jacobian as
- <math>\left(I-\bar{S}_U(\bar{U},\Delta t)\right)v\approx v -
The coarse time-stepper can also be used to accelerate the simulation of equation \eqref{coarse} over large (macroscopic) time intervals. Let <math>\Delta t\gg \delta t</math> be a large time step (commensurate with the slow coarse dynamics). Denote by <math>\bar{U}^n\approx \bar{U}(n\Delta t)</math> the numerical approximation of the coarse solution <math>U(t)</math> and by <math>\bar{U}^{k,n}</math> the iterates of the coarse time-stepper at <math>t_{k,n}=n\Delta t+k\delta t</math> (<math>\bar{U}^{0,n}\equiv\bar{U}^n</math>). (A slight abuse of notation: it is only necessary to lift the coarse solution for <math>k=0\ ,</math> since the corresponding fine-scale state for <math>k>0</math> is seamlessly available inside the coarse time-stepper.) One can then approximate <math>\bar{U}^{n+1}</math> via extrapolation, for instance,
- <math>\label{pfe}
Here <math>\bar{F}</math> approximates <math>F\ ,</math> e.g. via finite differencing
- <math>\label{pfe2}
The equation-free algorithm in \eqref{pfe} and \eqref{pfe2} is called coarse projective forward Euler, and it is the simplest instantiation of this class of coarse integration methods. The procedure is illustrated in<figref>Equation-free_modeling_projective.png</figref> for the example \eqref{eq-example}.
On the right the results are shown for coarse projective
integration with <math>\Delta t = 4 \delta t\ .</math> Intervals of full
simulation are in blue. Magenta dotted lines represent projections
in time, and the red dashed line is the full fine-scale simulation.
A phase space view is seen on the left. Note that, while <math>u_1</math> is
projected in time corresponding to the estimated time derivative,
<math>u_2</math> is reset during the lifting, so that the new initial
condition is off the slow manifold.
The first <math>k</math> steps with the coarse time-stepper serve two purposes. First, they reflect that, at the fine scale, we must allow for healing of the lifting errors (that perturb the fine-scale higher order moments) before we can start estimating coarse time derivatives. Second, in \eqref{pfe}, extrapolation is performed starting from <math>\bar{U}^{k,n}\ .</math> This is particularly useful if the coarse equation contains both fast and slow modes; in that case, if the number <math>k</math> of inner steps is chosen large enough, the size of the projective step is practically limited by stability restrictions on the slow modes only (Gear and Kevrekidis, 2003).
Higher order versions of \eqref{pfe} can be constructed in several ways. A straightforward idea is to use polynomial extrapolation (Gear et al., 2002), resembling a Taylor method. Implementations of coarse projective integration using Adams-Bashforth or Runge-Kutta as the outer solver are also possible, as are implicit versions, partially discussed in (Gear et al., 2002).
Projective integration is especially suited for problems with a large gap in their eigenvalue spectrum. However, extensions to deal with parabolic problems, have been proposed with properties similar to Runge-Kutta-Chebychev methods. Telescopic projective integration (Gear and Kevrekidis, 2003) is a recursive version of the method: the coarse projective integrator \eqref{pfe} itself is now considered as the inner integrator, around which a new projective integrator is wrapped. Another idea is to design a multistep state extrapolation method using the last points of a number of sequences of inner steps (Vandekerckhove et al., 2007).
In coarse projective integration, short time simulations by the fine-scale inner solver are used to explore long-time intervals at the coarse level. By analogy, spatially localized simulations performed in a number of small domains (teeth) separated by gaps, can be appropriately linked through spatial interpolation to enable the exploration of spatially extended systems at the coarse level.
A discussion of the resulting gap-tooth scheme can be found in (Kevrekidis et al., 2003), (Samaey et al., 2005),(Roberts and Kevrekidis, 2004); the combination of gap-tooth with coarse integration is called patch dynamics.
For simplicity, consider the unavailable coarse model \eqref{coarse} to be a partial differential equation (PDE) in one space dimension
- <math>\label{eq-coarse-pde}
where <math>\partial_x^k</math> denotes the k-th spatial derivative and the order <math>d</math> of \eqref{eq-coarse-pde} is known.
Given equation \eqref{eq-coarse-pde}, its finite difference method-of-lines discretization on a regular mesh can be written as
- <math>\label{eq-comparison}
where <math>U_i(t)\approx U(x_i,t)\ ,</math> <math></math>x_i\in\Pi(\Delta x):=\{0=x_0<x_1=x_0+\delta <math>\bar{U}^n="(\bar&#123;U&#125;_0(t_n),\ldots," t)="\frac&#123;\bar&#123;U&#125;_i^&#123;n+\delta&#125;&#45;\bar&#123;U&#125;^n&#125;&#123;\delta" t} </math>
=" Initial" <math>i="0,\ldots,N\" ,</math>
- <math>\label{eq-local-taylor}
- <math>\frac{1}{h}\int_{x_i-h/2}^{x_i+h/2}
="Boundary" conditions=" An" <math>H="O(\delta" center manifold theory.
="Open" directions=" *"
=" *" 3(3):1924.
="" reading=""
="See" links=" <a href="http://www.princeton.edu/che/people/faculty/kevrekidis/" target="_blank">website of Yannis Kevrekidis</a> <a href="http://www.cs.kuleuven.be/~giovanni" target="_blank">website of Giovanni Samaey</a>" also=" Category:Dynamical_Systems"></x_1=x_0+\delta> __AUTOLINKER{burst,bursting,memory} Category:Multiple_Curators