Deciding the Nature of the “Coarse Equation” through Microscopic Simulations: the BabyBathwater Scheme
Abstract
Recent developments in multiscale computation allow the solution of “coarse equations” for the expected macroscopic behavior of microscopically/stochastically evolving particle distributions without ever obtaining these coarse equations in closed form. The closure is obtained “on demand” through appropriately initialized bursts of microscopic simulation. The effective coupling of microscopic simulators with macrosocopic behavior embodied in this approach requires certain decisions about the nature of the unavailable “coarse equation”. Such decisions include (a) the determination of the highest spatial derivative active in the equation, (b) whether the coarse equation satisfies certain conservation laws, and (c) whether the coarse dynamics is Hamiltonian. These decisions affect the number and type of boundary conditions as well as the nature of the algorithms employed in good solution practice. In the absence of an explicit formula for the temporal derivative, we propose, implement and validate a simple scheme for deciding these and other similar questions about the coarse equation using only the microscopic simulator. Microscopic simulations under periodic boundary conditions are carried out for appropriately chosen families of random initial conditions; evaluating the sample variance of certain statistics over the simulation ensemble allows us to infer the highest order of spatial derivatives active in the coarse equation. In the same spirit we show how to determine whether a certain coarse conservation law exists or not, and we discuss plausibility tests for the existence of a coarse Hamiltonian or integrability. We argue that such schemes constitute an important part of the equationfree approach to multiscale computation.
I Introduction
It is often the case that a microscopic or fine description of a physical system is available, while we are interested in its macroscopic or coarse behavior. Consider, as an example, a biased random walk model for which the particle density asymptotically evolves according to a macroscopic law such as the Burgers equation. Typically, the study of macroscopic behavior starts with obtaining a closed PDElevel description (a “coarse equation”) for the time evolution of the expected or ensembleaveraged fields of a few, low order moments of the microstate phase space distribution. For our example, this would be the zeroth moment, the density field. Then, an array of mathematical and computational tools (numerical integration, fixedpoint algorithms etc.) can be brought to bear on the coarse equation.
Over the last few years, we have been developing a class of numerical algorithms which attempt to analyze the coarse behavior without ever obtaining the coarse equation in closed form [1]. The common character of these schemes is to use short, appropriately initialized bursts of microscopic simulations to estimate the quantities which, if the coarse equation was available, we would simply evaluate using the equation itself. Such quantities, estimated ondemand, include the time derivative of the evolving coarse fields, to be used in coarse projective integration [3], or the effect of the timeevolution operator for the implicit coarse Jacobian, to be used in NewtonKrylov type contraction mappings like the Recursive Projection Method (RPM) [4, 5] or in eigenvalue/vector computations. These methods are based on matrixfree large scale scientific computing, and we sometimes collectively refer to them as “equationfree methods”. What makes these computations possible is the assumption of a separation of time scales in the dynamics of the evolving microscopic distribution. Typically, one finds that the hierarchy of coupled equations involving higher cumulants of the microscopic distribution constitutes a singularly perturbed problem: higherorder cumulants become, in the course of simulation, quickly slaved to (become deterministic functionals of) the lowerorder cumulants. The consequences of slaving, realized in the computer as a blackbox, embody the closures that allow us to solve for the coarse behavior. Fundamentally it is no different than if the closures are expressed in closed form first and then evaluated later. This way of thinking and newly developed computing technology can in practice exceed the traditional approach in both accuracy and total cost, especially if the constitutive relation is nonlinear and multidimensional. An example of this type is given in [6]. An additional advantage of such methods is their ability to detect parametric regimes where the (number of moments used in the) present closure model is inadequate and hence appropriate refinements (including higher order moments) are necessary.
In the projective integration method [3] one takes advantage of the slow dynamics of the coarse variables to carry out only bursts of microscopic simulations connected via projections (in effect, extrapolations and/or interpolations) over gaps of time. In the same spirit of exploiting regularity, but now in space, we have developed the socalled gaptooth scheme [7, 8] by evolving the full process only in an array of small spatial boxes (the teeth) separated by empty regions (the gaps). Clearly, the two methods are closely related by the physics of the problem. Indeed we can have a combined gaptoothprojective integration scheme; this is the focus of another paper [8]. Here, we simply want to point out the fact that in the gaptooth method, the teeth communicate with each other via appropriate boundary conditions for the microscopic simulations performed inside them. And here lies the raison d’etre of this paper.
It is a wellknown fact that certain features of a given equation affect the nature of the appropriate numerical solver. A Hamiltonian dynamics problem, for example, is best integrated by a symplectic integrator; often, finite difference solvers of partial differential equations (PDE) are built to respect certain properties of the PDE, such as conservation laws. Most importantly, the highest spatial order of an evolution equation critically affects the types of boundary conditions leading to a wellposed problem.
In a completely analogous manner, the way in which the microscopic model is solved seperately in each tooth in the gaptooth scheme, and the boundary conditions applied to the edges of each tooth, must respect the nature of the unavailable equations and their order. Furthermore, gaptooth algorithms compatible with conservation laws (e.g. using fluxes to estimate temporal derivatives, see for example [9]) are predicated upon knowing that the unavailable equation possesses certain conservation laws.
When the closedform equation is available, some of these questions (e.g. the order of the highest spatial derivative in an evolution equation) can be answered by direct inspection. Other issues (such as the existence of conservation laws, or integrability) may, in the case of closed form equations, be relatively obvious, or may require a lot of work.
What we explore in this paper is the development of computerassisted methodologies to answer the above questions when closed form equations are not available. The idea is that we can probe the consequences of these answers on the dynamics of the unavailable coarse equations using microscopic, particle or agentbased, simulators by trying out large classes of appropriately chosen initial conditions. We will illustrate what we call the BabyBathwater algorithm on examples of particle systems realizing the Burgers and Kortewegde Vries (KdV) equations, for the task of inferring the highest order of spatial derivative on the righthandside, and for answering questions concerning coarse conservation laws.
The paper is organized as follows: In Section II we briefly present our illustrative particlebased example. In Section III we discuss the determination of the highest order spatial derivative active in the “unavailable equation”. In Section IV we explore the possible existence of conservation laws. In the concluding Section we discuss the scope and limitations of the procedure, as well as additional questions that may be addressed through this approach. An interesting “twist” about reverse coarse integration arises in discussing the exploration of possible “coarse Hamiltonianity” of the unavailable coarse equation.
Ii Numerical Experiment Setup
Our illustrative examples will be based on simple numerical experiments. We will first, as a sanity check, demonstrate the approach using a traditional numerical simulator of a known evolution equation as a “black box”. We will then substitute the simulator of the known continuum equation with a particlebased simulator, and repeat the procedure. Our first illustration will be the Burgers equation,
(1) 
as well as a particle based simulator constructed so that the evolution of its density resembles (at the appropriate limit, reproduces) the Burgers evolution. Since one of the issues to be explored is the number and type of boundary conditions in evolving the equation, our simulations must be possible without this a priori knowledge. We therefore use periodic boundary conditions (PBC) enforced on in all our exploratory simulations. One of the attractive features of the Burgers is that, for any initial profile , and even with PBC, the ColeHopf transformation [10] provides an analytical solution. The accuracy of the numerics can thus be checked directly.
A biased random walkerbased particle simulator mimicking the Burgers dynamics was also constructed to demonstrate the direct application of our procedure on microscopic, particle based solvers. A detailed study of the features of this particle model is reported elsewhere [11]. As a reference, the diffusion equation,
(2) 
has the wellknown microscopic realizations of Langevin dynamics or unbiased random walkers. It is not too difficult to conjure up a similar realization motivated by the Burgers equation using random walkers: a unit mass in the coarse description corresponds to walkers, where is a large integer constant. In the simulation, random walkers move on at discrete timesteps . At each step, (a) the walkers’ positions are sorted, (b) each walker checks out the position of the walker places ahead, , and places behind, (properly accounting for PBC, of course). The difference is inversely proportional to the local density of walkers, therefore (c) every walker moves by sampled from , a biased Gaussian distribution. The ’s are then wrapped around to , and the process repeats. This achieves a coarsegrained flux of as motivated by the (1) by assigning each walker a drift speed of . Quantifying the approximation of the Burgers evolution is an interesting subject that we take up separately in [11]; this is not, however, an important issue for this paper. It is only for benchmarking purposes that the relation to a known macroscopic equation is brought up. One can start by presenting a microscopic evolution law, without knowing anything about its corresponding coarse equation, and apply our algorithms on it directly.
Relating the fine with the coarse description requires the use of lifting and restriction operators[5]. Lifting constructs particle distributions conditioned on some of their lower moment fields (here the zeroth moment, or coarse density field); it clearly is a onemany operator, and several microscopic “copies” of a given macroscopic initial condition are often required as discussed in detail in [1, 5]. The restriction (here computing moment fields of a given particle distribution) is a form of projection. Clearly the restriction (back to coarse variables) of a lifting of these coarse variables should be the identity (or close to it, due to noise effects). The lifting and restriction operators we constructed for this work, with interpreted as the coarse density on PBC, are given in Appendix A. Fig. 1 shows a result of the “reversibility test”: we randomly generate a coarse density , lift it to a random walker distribution, then restrict back to , and observe the very good agreement between with . Notice that the only point where this agreement may be less satisfactory is close to local maxima or minima, where the derivative changes sign.
While the proximity of our particle scheme to the Burgers evolution is not the issue in this paper, we briefly illustrate the correspondence of the evolution of an initial profile through the two approaches. Fig. 2 shows the analytical solution obtained through the ColeHopf transformation; it also contains the result of a particle simulation, after the configurations have been processed by the operator of Appendix A to extract the coarse density field estimate. A small value of viscosity is picked to accentuate the behavior of steepening wavefront with time. The microscopic simulation clearly captures the important features of the coarse behavior. Ensemble averaging with the same initial condition in coarse field would reduce the error, but as we can see, even a single microscopic simulation using a reasonable number of particles may still perform quite efficiently.
Lastly, we mention the existence of particle methods to solve partial differential equations (PDEs) such as the Korteweg  de Vries (KdV) equation:
(3) 
which can be formulated as conservation laws [12, 13]. In this paper we use our construction above for the Burgers example. We should highlight once more that our ultimate purpose is not to construct particle solvers of given equations, but rather to decide on features of the unavailable coarse equations for given microscopic schemes.
Iii Identifying The Highest Spatial Order of Coarse Variables
As we mentioned in Section I, system identification lies at the heart of the equationfree approach. Here, we suppose the coarse dynamics follows a certain timeevolution equation of the form:
(4) 
that is unavailable to us. We have already identified (the “coarse variable” for which we believe that a coarse deterministic equation exists in closed form), and constructed the lifting and restriction operators that connect macro/micro descriptions. We seek a general approach to decide qualitative questions, such as (a) what is , and (b) whether can be written as , without having in closed form. One important motivation for this lies in that “production run” simulations of the problem via equationfree computation (for example through the gaptooth scheme) do not require knowledge of , but are affected by the knowledge of (through “teeth” boundary conditions). What we do have is a microscopic simulator embodied in a computer code that can be initialized at will; the physical details of the microscopic code are both extremely important (that is where the “underlying physics” lies) and  for our purposes  irrelevant: we will use the microscopic simulation code as an “inputoutput” (I/O) black box. By probing the coarse I/O response of the black box, the question that we would like to address is whether we can decide on (a) and/or (b).
It may appear initially that we are trying to answer a circular question: in order to probe the coarse inputoutput response of a microscopic simulator we need to run it, and in order to run it we need wellposed boundary conditions, which  among other factors  depend on (a) and (b). To cut the knot, we use (for the decision stage exploratory runs) the Bornvon Karman periodic boundary conditions. We are going to assume that the microscopic simulations can be carried out in PBC, which is an option prevalent among microscopic simulators. This enables us to probe the system’s response to only the initial profile input.
The socalled babybathwater identification scheme works as follows:

Take an integer , starting from .

Pick a random point in the spatial periodic box.

Generate random numbers, designated as , , , .., of .

Generate a conditionally random profile compatible with the PBC and consistent with the above , , , .., requirements. This can almost always be accomplished, for example, by summing sine and cosine harmonics of the PBC:
(5) with . Because we have coefficients, even though there are constraints to satisfy, we still have some random degrees of freedom left in (5). In practice, this initialization can be accomplished by applying conjugate gradient minimization of the dimensional residual norm starting from a random vector.

Lift of (5), run it in the microscopic simulator for time , restrict it back to , and estimate:
(6) Note that here instead of is used in the finite difference. This will cancel some internal noise from the lifting and restriction operations.

Repeat step (v) times to obtain an ensemble averaged to reduce the microscopic noise.

Repeat step (iv) times, collect the estimates:
(7) compute the sample variance .

Repeat step (ii) times, compute the averaged sample variance .

Go back to step (i), . is identified when from to , the averaged sample variance decreases drastically to practically .
Fig. 3 shows such families of constructed initial profiles with progressively more controlled initial derivatives. The basic idea is very simple: even though could have complicated functional dependencies on , if they are all fixed, should have no dispersion even as are varied randomly. The “critical integer order” is identified when the variance at controlled derivatives , , , .., jumps to a finite value; we then have already thrown out the “baby” (the highest relevant spatial derivative ) with the “bathwater” (the higher, nonrelevant ones).
It is important to recognize that the time derivative estimation (6) does not occur instantaneously. A short “healing” period should elapse, during which the higher cumulants of the lifted phase space (microstate) distribution become functionals of the lower order, slow governing cumulants. This separation of time scales, which fundamentally underlies the existence of a deterministic coarse equation closing with the lower cumulants, is discussed in more detail in [1].
As a sanity check, this algorithm is also applied to a traditional continuum PDE timestepper “black box” first. Fig. 4(a) and 4(b) show the results of applying our decision scheme to forward Euler finitedifference PDE solvers of the Burgers and KdV equations, respectively. A spatial mesh of is adopted, and we define
(8)  
(9)  
(10)  
(11) 
We use,
(12) 
to integrate the Burgers equation forward, and
(13) 
to integrate the KdV equation forward. is obtained by cubic spline over , and is evaluated by finite differences, same as in (6). As can be seen in Fig. 4(a) and 4(b), is identified to be using the Burgers PDE timestepper and using the KdV PDE timestepper: the variances drop by more than four decades in both cases when going from to controlled derivatives. To see where the remaining “noise” comes from, note that,
(14) 
and clearly has higherthan spatial derivative dependencies, which are, however, scaled by compared to the leading term. Thus the sample variances should drop by for , which explains the observed magnitude of the fourdecade decrease.
We then apply the identification scheme to the microscopic simulator of Section II, with the lifting and restriction operators constructed in Appendix A. The results are shown in Fig. 5. Under favorable conditions such as and , it takes about minutes of computer time on a single 1GHzCPU personal computer to obtain a reasonably good microscopic noise reduction so the variance drops by about decades going from to . Under unfavorable conditions such as and , it can take up to minutes of computer time to obtain the same decades drop. Compared with the deterministic finitedifference PDE timesteppers, identification of a microscopic simulator is undoubtedly much more computationally intensive, even though fundamentally there is no difference between the two “black boxes”. The problem of microscopic noise reduction is a persistent issue among all “equationfree” methods including bifurcation [1], projective / gaptooth integration [7, 8], and identification, and calls for a unified treatment. This “one time” decision, performed at the beginning of studying a problem, will critically affect subsequent production runs of the microscopic simulator.
Here, one must pay special attention to the rank () of the restriction operator (see Appendix A). As can be seen in Figs. 1 and 2, our proposed restriction operator satisfies the constraint of reversibility and also accurately represents the profile’s longtime evolution. However, these merits do not guarantee automatically good shorttime estimates by finitedifference. Special attention must be paid to the restriction operator : for example, if the highest harmonic in (5) for is , then with , we can get good reversibility test of . Unless we use for restricting the Burgers microscopic dynamics, however, we would not get a good estimate of , because the nonlinear interaction in (1) creates higher harmonics in up to . If is still used, it is equivalent to forcing a leastsquare projection of a vector to a subspace, which may work well enough in the long term, but is too inaccurate for shortterm finite difference estimates. Unless this is taken care of, the estimate using our is found to not even be superior to a crude bincount density estimator with binwidth about , as is the shortest wavelength in .
Lastly, we note that (4) represents a wide category of coarse dynamics; those with higher timederivatives and mixed derivatives can be converted to a multivariate version of (4) and the babybathwater identification scheme will still, in principle, work. A notable exception is the incompressible fluid dynamics case, where the soundspeed is infinite and the pressure plays the role of a global Lagrange multiplier. The incompressible fluid model is but a mathematical idealization of a certain physical limit. It is nonetheless useful and important enough, that the fact that it is not directly amenable to the babybathwater identification is worth mentioning. In general, the babybathwater identification presented here will not work for dynamics with instantaneous remoteaction over macroscopic lengthscales, such as,
(15) 
for which it is easy to show that correlates with infinite number of local spatial derivatives .
Iv Identifying Conservation Laws
In section III above we address the concern of how to identify the highest spatial derivative of an unavailable coarse equation of the type (4). It is natural to try to decide other qualitative questions: for example, whether the coarse dynamics conserve a specific quantity,
(16) 
or not. In the simplest case, we ask whether is conserved. We note that is equivalent to asking whether the RHS of (4) can be written as,
(17) 
or not. Alternatively, we ask whether there exists such that,
(18) 
for arbitrary ,. Whereas in section III we try to identify features of through (4), here we can try to identify consequences of and its features through (18). The process of the babybathwater identification can be carried over; the only difference is that it is going to be a boundary scheme. In one dimension, the boundary sheme reduces to a twopoint scheme as follows:

Take an integer , starting from .

Pick two random points and in the PBC.

Generate random numbers, which are to be designated , , , .., and , , , .., , of .

Generate a conditionally random profile compatible with the PBC that is consistent with the above , , , .., and , , , .., requirements. This can always be done by (5) with . As we discussed above, we have coefficients, even though there are constraints to satisfy, we still have some random degrees of freedom left in .

Lift of (5), run it in the microscopic simulator for time , restrict it back to , estimate:
(19) 
Repeat step (v) times to obtain an ensemble averaged to reduce the microscopic noise.

Repeat step (iv) times, collect the estimates:
(20) compute the sample variance .

Repeat step (ii) times, compute the averaged sample variance .

Go back to step (i), . A conservation law is positively identified when going from to , the averaged sample variance decreases drastically to practically .
Fig. 6 plots families of initial profiles thus constructed with progressively more controlled initial derivatives. Fig. 7(a) and 7(b) show the results of applying the identification scheme to the Burgers finitedifference PDE timestepper (12) and the KdV finitedifference PDE timestepper (13), respectively. is identified to be for (12) and for (13).
Two comments are in order: first, we probe the consequences of conservation (i.e. that boundary fluxes are the only cause of change for the conserved quantity in a domain); second, we obtain (as a sideproduct) the highest spatial derivative of the conserved quantity in the constitutive equation for the flux. It is important to note that if the procedure progressively returns negative answers (e.g., if the sample variance is nonzero for a given number of controlled derivatives), this does not imply that a conservation law does not exist. It only implies that a conservation law of the class encompassed by our equation (4) with spatial derivatives up to the tested order does not exist. In that sense, our procedure provides sufficient confirmation, but its success is not necessary for a conservation law in a different class to prevail. Nevertheless, the class we consider is wide enough to encompass many known examples and problems of interest to applications.
V Discussion
In section IV we proposed methods to check whether a coarse quantity (such as the mass corresponding to the coarse density ) is conserved, without knowledge of the coarse evolution equation. The obvious question that arises is how do we know which to check? The path that we suggest here, in the equationfree setting, is to examine the consequences of conservation laws. For example, consider the conservation of (linear) momentum. An equivalent statement, through Noether’s theorem [14] is the existence of translational invariance. If we numerically establish the latter, then we can claim the former. Let us then consider initial conditions to the available integrator which are shifts of an original profile e.g., , , , etc. Then if we time evolve the problem, using our microscopic timestepper, and the equation is translationally invariant, upon reaching the integration reporting horizon, we can backshift the profile (by the original shift amount). If all backshifts provide an identical profile, we can conclude translational invariance and hence linear momentum conservation. An additional note of caution is that the examination of such consequences is relevant when Noether’s theorem applies, hence when there is an underlying Lagrangian/Hamiltonian structure in the problem (we discuss separately the issue of Hamiltonian nature below). Notice, however the modulo the proviso of “Hamiltonianity”, this methodology can be used to establish additional dynamical invariants, e.g. the invariance with respect to phase of the evolution of a field can be related to norm invariance etc.
Establishing an underlying Hamiltonian structure in a sense proceeds in a similar fashion through its correlation with invariance with respect to time reversal. The crudest way to examine this is by simply running the integrator with a negative timestep (if that option is available). A more refined way to check the same symmetry is by examining computations of the spectrum (e.g., eigenvalues) of linearization of the coarse PDE. In particular, a straightforward consequence of the Hamiltonian nature is that all linearization eigenvalues should come in quartets, namely if is an eigenvalue, then so are , where denotes complex conjugation. It is fortunate that timestepper based numerical analysis techniques for the numerical approximation of the leading spectrum of such a linearization are welldeveloped for the case of large scale continuum simulations (see for example [18, 19, 21, 22, 23, 20]). If an eigenvalue of the linearization is identified, matrixfree eigencomputations with shift can be used to explore the existence of the eigenvalue (in general, real eigenvalues will come in pairs and complex conjugate eigenvalues in quartets).
While coarse timereversibility can be answered by exploring the spectrum of the linearization, it raises the interesting question of how to integrate backward in time with the microscopic code. Consider a molecular dynamics configuration with a certain set of velocities at time zero, and the same molecular configuration with “flipped” velocities. A well known (and testable) consequence of microscopic reversibility and the “molecular chaos” ansatz is that, whether we integrate the molecular dynamics equations forward or backward in time starting from a randomly picked phase point, we will get “the same” forward in time evolution of the coarse macrosocpic observables. It is interesting, however, that the coarse projective integration techniques in an equationfree context can be used to attempt integration of the coarse variables backward in time (under appropriate conditions about the spectrum of the unavailable equation) as follows: Consider the lifting of a particular coarse initial condition to consistent molecular realizations; flipping the molecular velocities for these realizations will not affect the coarse procedure. We then evolve microscopically the molecular configurations (whether with the original or with flipped velocities) forward in time long enough for the higher moments to heal, say for a time . We now estimate the time derivative of the healed coarse variables from the restriction of the “tail end” of the molecular trajectories, and then take a large, macroscopic Euler step backward in time for the coarse variables. We lift again, run molecular dynamics forward or backward microscopically, estimate the coarse forward time derivative, and take another coarse backward step. This procedure can of course be done in a much more sophisticated way as far as the coarse backward time step is concerned  algorithms like RungeKutta or Adams method can be combined with the MD computations to integrate density expectations backward in time for the coarse equations on the “slow manifold”. This “seesaw” forwardbackward coarse integration procedure can also be used on stiff systems of ODEs and even dissipative PDEs under the appropriate conditions to evolve trajectories backwards on a slow manifold. The numerical analysis of these algorithms in the continuum case is an interesting subject in itself, and we are currently pursuing it [24]. It is interesting that the technique, in the molecular dynamics case, can be used to coarsely integrate backward in time on a free energy surface, and thus help molecular simulations escape from free energy minima; we have already confirmed this in the case of Alanine dipeptide folding in water at room temperature through molecular dynamics simulations [25].
Finally, a more complicated question than checking the existence of one (a specific, and hence related to a specific invariance, in accordance with the above discussion) integral of the motion is the one of integrability. The latter necessitates infinite integrals of the motion, normally established by means of identifying Lax pairs and using the inverse scattering transformation machinery [15]. However, one can also use in this case consequences of integrability to establish it. For instance, in recent work [16] it was qualitatively argued (and verified through numerical experiments in different settings) that a feature particular to integrable Hamiltonian systems is the presence of double continuous spectrum eigenvalues, when linearizing around a (coarse PDE) solitary wave under periodic boundary conditions. These as well as other criteria (such as the existence of point spectrum eigenvalues in the spectral gap [17]) can also be (conversely) used to potentially rule out the existence of integrable structure. In short, the spectral properties of the coarse PDE linearization can be used to establish or disprove not only the Hamiltonian (see above), but also potentially the integrable nature of the flow. While these are just initial thoughts towards attempting to decide vital questions about the nature of the unavailable closed equation, it is important to note that, what is computationally involved is a timestepper based identification of facts about the spectrum of the linearization of an operator. This “computational technology” is quite mainstream in the case of large scale continuum simulators, and can be straightfowardly adapted to the case of coarse timesteppers in conjunction with the liftingrestriction steps. Variance reduction will clearly be the most significant step in the wide applicability of these and similarspirited approaches.
Acknowledgements The authors would like to thank Frank Alexander for stimulating discussions. This work was partially supported by the National Science Foundation (IGK,PGK) AFOSR (CWG,IGK) and the Clay Institute (PGK).
Appendix A Coarse Density Lifting / Restriction Operators
For a coarse field , the lifting operator generates a microstate : . Similarly, for a microstate , the restriction operator returns a coarse field estimate : . Both and can be onetoone or onetomany operators, but we demand that asymptotically when the wavelength of is large enough compared to the microscopic length scale [5]. For 1D coarse density field under PBC, we use the following , operators for the sake of definiteness in numerical experiments, even though their construction is not unique.
The lifting operator , :

Estimate . Pick that is “safely” greater than , for example .

Define . Create particles with each independently drawn from uniform distribution on .

Go to each particle , randomly decimate it with probability . Count the total number of surviving particles .

Compute quadrature,
(21) randomly round to or such that . Randomly pick particles out of the survivors and decimate them. We now have a set of particles , totally numbered either or .
The restriction operator , :

Define microscopic density function,
(22) and corresponding cumulant function,
(23) Clearly, at the the first, second, third particle positions , , , , , , etc. And we have , .

Define a residual function ,
(24) which is the difference between and the cumulant of a homogenized particle gas background. The idea is that , so it is a periodic function and can be approximated by,
(25) In fact, a sound strategy is to leastsquare fit (its , coefficients) to at ’s, where is the sorted list of . can be easily evaluated at ’s, noting the last sentence of step 1.

The coarse density estimate can be obtained by taking the derivative of ,
(26)
It is worth noting that although the constructed depends on , it satisfies the particle number conservation exactly because the finite harmonics all integrate to zero and only the background contribution remains. In fact, also satisfies exact particle number conservation to the original under probabilistic average. Further, one can show exactly for in the first harmonics subspace.
References
 [1] I.G. Kevrekidis, C.W. Gear, J.M. Hyman, P.G. Kevrekidis, O. Runborg, C. Theodoropoulos, “EquationFree Multiscale Computation: enabling microscopic simulators to perform systemlevel tasks,” submitted to Communications in the Mathematical Sciences;
 [2] A.G. Makeev, D. Maroudas, I.G. Kevrekidis, “Coarse stability and bifurcation analysis using stochastic simulators: Kinetic Monte Carlo examples,” J. Chem. Phys. 116 (2002) 1008391; A.G. Makeev, D. Maroudas, A.Z. Panagiotopoulos, I.G. Kevrekidis, “Coarse bifurcation analysis of kinetic Monte Carlo simulations: A latticegas model with lateral interactions,” J. Chem. Phys. 117 (2002) 822940.
 [3] C.W. Gear, I.G. Kevrekidis, “Projective Methods for Stiff Differential Equations: problems with gaps in their eigenvalue spectrum,” SIAM Journal on Scientific Computing, in press; can also be obtained as NECITR 2001–29 at http://www.neci.nj.nec.com/homepages/cwg
 [4] G.M. Shroff, H.B. Keller, “Stabilization of Unstable Procedures  the Recursive Projection Method,” SIAM Journal on Numerical Analysis 30 (1993) 10991120.
 [5] C.W. Gear, I.G. Kevrekidis, C. Theodoropoulos, “Coarse Integration/Bifurcation Analysis via Microscopic Simulators: microGalerkin methods,” Computers and Chemical Engineering 26 (2002) 941963.
 [6] J. Li, K.J. Van Vliet, T. Zhu, S. Yip, S. Suresh, “Atomistic mechanisms governing elastic limit and incipient plasticity in crystals,” Nature 418 (2002) 307310. The quasicontinuum method is a superset of IPFEM. See, for instance, R. Phillips, D. Rodney, V. Shenoy, E. Tadmor and M. Ortiz, Model. Simul. Mater. Sci. Eng. 7, 769 (1999).
 [7] I.G. Kevrekidis, “Coarse Bifurcation Studies of Alternative Microscopic/Hybrid Simulators,” Plenary Lecture, CAST Division, AIChE Annual Meeting, Los Angeles, 2000. Slides can be obtained at http://arnold.princeton.edu/yannis/
 [8] C.W. Gear, J. Li, I.G. Kevrekidis, “Patch Dynamics on the Burgers Equation,” to be submitted to Physical Review Letters.
 [9] W. E, B. Engquist, “The Heterogeneous MultiScale Methods,” Communications in Mathematical Sciences, in press (2003).
 [10] J.D. Cole, “On a quasi linear parabolic equation occurring in aerodynamics”, Quart. Appl. Math. 9 (1951) 225236; E. Hopf, “The partial differential equation ”, Comm. Pure Appl. Math. 3 (1950) 201230.
 [11] J. Li, C.W. Gear, I.G. Kevrekidis, “Particle solution of the Burgers equation and microscopic correlations,” to be submitted (2003).
 [12] A. Chertock, D. Levy, “Particle methods for dispersive equations,” Journal of Computational Physics 171 (2001) 70830.
 [13] A. Chertock, D. Levy, “A particle method for the KdV equation,” Journal of Scientific Computing 17 (2002) 491499.
 [14] V.I. Arnold, Mathematical Methods of Classical Mechanics, SpringerVerlag (New York, 1989).
 [15] M.J. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, SIAM (Philadelphia, 1981).
 [16] P.G. Kevrekidis, N.R. Quintero, “Using the Continuous Spectrum to Feel Integrability, the Effect of Boundary Conditions,” (preprint).
 [17] P.G. Kevrekidis, “Integrability Revisited: A Necessary Condition”, Phys. Lett. A 285, 383 (2001).
 [18] K.N. Christodoulou, L.E. Scriven, “Finding leading modes of a viscous free surface flow: an asymmetric generalized eigenproblem,” J. Sci. Comput. 3 (1988) 355406.
 [19] I. Goldhirsch, S.A. Orszag, B.K. Maulik, “An efficient method for computing leading eigenvalues and eigenvectors of large asymmetric matrices,” J. Sci. Comput. 2 (1987) 3358.
 [20] Z. Bai, J.W. Demmel, J. Dongarra, A. Ruhe, H. van der Vorst eds. Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide (SIAM publications, Philadelphia, 2000).
 [21] W.E. Arnoldi, “The principle of minimized iterations in the solution of the matrix eigenvalue problem,” Q. Appl. Math. 9 (1951) 1729.
 [22] C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators,” J. Res. Natl. Bur. Stand. 45 (1950) 255282.
 [23] L.S. Tuckerman, D. Barkley, “Bifurcation analysis for timesteppers,” in Numerical Methods for Bifurcation Problems and LargeScale Dynamical Systems, IMA Volumes in Mathematics and its Applications 119, p.453466, E. Doedel and L.S. Tuckerman eds., (Springer, New York, 2000).
 [24] C.W. Gear, I.G. Kevrekidis “Reverse Projective Integration” in preparation.
 [25] G. Hummer, I. G. Kevrekidis, “Coarse Molecular Dynamics of a Peptide Fragment: Free Energy, Kinetics and LongTime Dynamics Computations,” submitted to J.Chem.Phys (2002).