Computational Coarse Graining of a
Randomly Forced 1D Burgers Equation
Abstract
We explore a computational approach to coarse graining the evolution of the largescale features of a randomly forced Burgers equation in one spatial dimension. The long term evolution of the solution energy spectrum appears selfsimilar in time. We demonstrate coarse projective integration and coarse dynamic renormalization as tools that accelerate the extraction of macroscopic information (integration in time, selfsimilar shapes, and nontrivial dynamic exponents) from short bursts of appropriately initialized direct simulation. These procedures solve numerically an effective evolution equation for the energy spectrum without ever deriving this equation in closed form.
pacs:
47.11.St, 47.11.j, 47.27.efI Introduction
The behavior of physical systems is frequently observed, and modeled, at different levels of complexity. Laminar Newtonian fluid flow is a case in point: it is modeled at the atomistic level through molecular dynamics simulators; at a “mesoscopic” level through Lattice Boltzmann models; and (in physical and engineering practice) through continuum macroscopic equations for density, momentum and energy: the Navier–Stokes. What makes the latter, coarse grained description possible is an appropriate closure, in this case Newton’s law of viscosity, which accurately models the effect of higher order, unmodeled quantities (such as the stresses) on the variables (“observables”) in terms of which the model is written (the velocity gradients). Such closures are often known from physical and engineering observation and practice long before their mathematical justification becomes available (in this case, through Chapman–Enskog expansions and kinetic theory). It is, of course, tempting to consider a situation in which the fine scale description is the Navier–Stokes equations themselves, and the coarse grained description is an evolution equation for some observables of turbulent velocity fields, or possibly of ensembles of such fields. Discovering the appropriate observables and deriving such evolution equations has been a longstanding ambition in both physics and mathematics. Our goal here is much more modest: we present a very simple illustration of a problem for which a fine scale description is available, and for which (based on extended observations of direct simulation) we have reason to believe that a coarse grained evolution equation exists – yet it is not available in closed form. We illustrate how (with a good guess of the appropriate coarse grained observables) we circumvent the derivation of a closed form coarse grained model, but still perform scientific computing tasks at the coarse grained level. The numerical tasks we demonstrate are coarse projective integration, which accelerate the (coarse grained) computations of the system evolution, and coarse dynamic renormalization, which (when the coarse grained evolution is selfsimilar, as the case appears to be here) targets the computation of the selfsimilar solution shape and the corresponding exponents. This is an illustration of the socalled EquationFree framework for coarse grained scientific computation in the absence of explicit quantitative closures and the resulting coarse grained evolution equations manifesto (); pnas (); manifesto_short ().
We consider the one dimensional in space, randomly forced Burgers equation,
(1) 
subject to periodic boundary conditions over the domain . We start with random very low energy initial conditions. We are interested in the relatively long term and large scale properties of this system in the inviscid limit, when subject to a random forcing at the small scales. In order to achieve this inviscid limit numerically, we consider a general dissipation term of the form following chekhlov95 (). The coefficients and are chosen so that the dissipation acts only at the highest wavenumbers; here, their values are and . This choice of dissipation is based on the assumption that the universal ranges of the system are not sensitive to a particular choice of the parameters and chekhlov95 (). The white–in–time random force is defined in Fourier space by
(2) 
where and are the spatial and temporal frequencies. The forcing term has a peak at wavenumber and dies off at large wavenumbers; the dissipation term is essentially zero up to around wavenumber and only becomes important at much larger wavenumbers; see Fig. 1.
Before we start our computations, we briefly recall certain known features of the equations and their dynamics. It is quite remarkable that equations (1,2) have nontrivial asymptotic behavior in both the ultraviolet (UV, ) and the infrared (IR, ) limits. Attempts to use (1,2) as a simplified 1–D model for the investigation of small scale features of 3–D turbulence, aimed at recovering in the small scale limit, were based on deep similarities between the two: the total energy in this system is and, as in 3–D turbulence, the inertial range where is the dissipation wave number, is characterized by a constant energy flux , which, in the interval , is reflected in the analytic behavior of the thirdorder structure function
(3) 
Even in decaying high Re–turbulence this equation is correct: is very small. At first glance, due to the shock formation leading to the energy spectrum (rather than with the exponent ), these attempts failed. However, the model revealed a nontrivial biscaling behavior of the structure functions for and for frisch (). The model also generated asymmetric probability densities with algebraically decaying tails for large amplitude fluctuations chekhlov95 (); vyPRL96 (). The properties of this tail attracted the attention of various groups and are reasonably well understood in terms of the dissipation anomaly introduced by polyakov95 (); they have also been analytically obtained in e97 (); e99 (); bec00 (). The Burgers equation stirred by a force with the spectrum does lead to the logarithmically corrected Kolmogorov energy spectrum discovered in chekhlov95b () and does not have multifractality mitraArXiv ().
The IR properties of (1,2) in the range , on which we focus in this paper, are also quite interesting and nontrivial. It was determined, on the basis of the oneloop calculation (also numerically verified in vyPRL88 ()) that in the limit , (1,2) is equivalent to the Burgers equation driven by a random force with –spectral density, which is called the KPZ equation, governing various interfacial phenomena kpz (). The renormalization group methodology applied to the KPZ equation forster77 () led to the equilibrium energy spectrum (here ) and to a nonGaussian dynamic fixed point with nontrivial dynamic scaling exponents,
(4) 
indicating strongly nondiffusive behavior: . These predictions, including dynamic scaling, have been numerically verified in vyPRL88 (). Numerical simulations in the IR range become progressively slower; our computerassisted approach aims at accelerating such simulations (see also kesslerPRE ()).
Our numerical simulator integrates equation (1) in time using a pseudospectral scheme with a Fourier mode spatial discretization. Thus, the smallest scale resolved by our simulator is . The nonlinear term is computed using the fast Fourier transform. The temporal discretization treats the linear term exactly and uses a third order Runge–Kutta–Wray scheme for the nonlinear term. The random forcing term in Fourier space is , where is a random number generated from a Gaussian distribution with zero mean and unit standard deviation, and is the time step of integration. The various parameters chosen are , , , , , and . An ensemble of 16 distinct realizations was computed, and the results averaged to achieve good statistics.
We begin by observing direct numerical simulations of the evolution of equation (1) over time. Fig. 2 shows such a transient, starting from the zero initial condition, , for all . The energy spectra are plotted after , , , and time steps; the insets show the corresponding 1D velocity fields. Visually, we can discern no clear trend in the evolution of the velocity fields themselves, which appear random; yet, we easily see a definite smooth structure in the evolution of the energy spectrum.
For , the spectrum seems to quickly achieve a steady state profile within less than 1000 time steps: the energy provided by the noise is balanced by the strong dissipation at high wavenumbers. It is much more interesting to observe for ; here the spectrum quickly acquires a “cornerlike” shape which appears to evolve smoothly in time, and moves towards lower wavenumbers while maintaining this shape. This “traveling frontlike” evolution progressively slows down as the “corner” evolves towards lower wavenumbers.
Theoretical considerations forster77 (); vyPRL88 () indicate that in the range , the evolving energy spectrum consists of three regimes. For , the spectrum is not universal and depends on the properties of the forcing function. For , (here ); and in the timedependent largescale limit , the energy spectrum is given by a universal relation , which we observe in our simulations. Superficially, the temporal evolution we observe resembles the timedependence of a two dimensional flow driven by a small () random force, eventually leading to the energy spectrum piling up on the smallest (integral) wave number and leading to the IR asymptotics . We stress that the situation we are interested in may be quite different: in 2D, due to enstrophy conservation in the infrared range (), there exists a constant, nonzero, energy flux toward large scales, leading to the so called “inverse cascade”. The thirdorder structure function in a 2D flow is . In the one dimensional flow we are interested in here, the large scale thermodynamic equilibrium is fluxfree and must be equal to zero, provided the system is large enough, that is, .
We do not work at this very long time regime; we rather observe the evolution at times intermediate between this regime and the (relatively) short time it takes for the spectrum to stabilize for . This is the regime in which smooth and, as we argue below, apparently selfsimilar evolution of the energy spectrum prevails.
Based on the study of many transients like those in Fig. 2, we decided to observe the system dynamics in terms of the energy spectrum only, which is defined as
(5) 
where is the spatial Fourier transform of the velocity field , and denotes an ensemble average (here over the 16 realizations). A formal diagrammatic expansion for the evolution of the energy spectrum corresponding to (1,2) can be developed wyld61 (). In the limit , taking into account that all oddorder correlation functions rapidly tend to zero, the solution can be accurately represented in terms of an infinite series involving complicated integral expressions with kernels which are convolutions of the energy spectrum only. Based on this, we assume that an evolution equation of the general form
(6) 
exists and closes for the (ensemble averaged) . We do not know the nature of the operator – it may not be local (in –space), and the equation may not be a partial differential equation – yet the implication is that knowing is enough to predict the (realization ensemble averaged) for later times. If an evolution equation closes with respect to the expected spectrum only, then the remaining degrees of freedom (the phases of the velocity field Fourier coefficients) do not explicitly appear in this equation. This implies that either these variables do not couple to the spectrum evolution, or that their expected dynamics become quickly slaved to the spectrum evolution (with a relaxation time that is small compared to the spectrum evolution time scales). We test this latter ansatz by observing the solution ensemble not only through the spectrum but also through the evolving third order structure function
(7) 
For completeness, the “extremely long term” steady state of and , obtained after integrating for time steps is included in Fig. 3. Fig. 3(a) suggests that reaches a constant value of in the range . The structure function , plotted in Fig. 3(b), is practically zero in the range . The steady state value of is nonzero for , due to the force acting primarily at small scales (corresponding to ) and strong dissipation at the smallest scales.
The remainder of the paper is organized as follows: In Section II we discuss simulations suggesting an apparent dynamic selfsimilarity in the intermediatetime evolution of the spectrum. Section III briefly outlines equation free computational techniques, and demonstrates coarse projective integration of the unavailable equation for the expected spectrum evolution. Section IV then shows how to implement a dynamic renormalization fixed point algorithm that converges on the selfsimilar shape and pinpoints the temporal similarity exponent of its evolution. We conclude with a brief discussion.
Ii Apparent Dynamic Selfsimilarity
Here we argue that the transient evolution of exhibits dynamic selfsimilarity. Fig. 4 clearly shows that the qualitative shape of the spectrum in the region remains effectively unchanged, while gradually “travelling” to the left, slowly filling up the lower wavenumbers. We make this observation more precise as follows: We define a “fulcrum point” as the leftmost end of the apparently steady region of ; here . As shown in Fig. 4(a), we draw straight line segments of various slopes (logarithmic in as well as in ) starting at the fulcrum and intersecting successive computed spectra. Let be the distance of the point of intersection of each such line with the spectrum at time . Our ansatz is that the distance becomes uniformly stretched in (logarithmic) time. This is tested as follows: let the stretching factor along a line of slope be , where is a reference time. Fig. 4(b) is a plot of versus at various times in the order listed in the figure caption; it strongly suggests that the spectrum gets stretched by the same amount along its entire profile (notice that no intersections arise with ). Fig. 4(c) shows the time evolution of the stretching factor (averaged over the lines of different slopes), clearly showing linearity in logarithmic time. Thus, the transient evolution of appears dynamically selfsimilar: the spectrum in the wavenumber region becomes uniformly stretched in logarithmic time, while maintaining its shape. We take advantage of this information in speeding up computations in the next section.
Iii Coarse Projective Integration
Based on our simulation observations above, we now hypothesize that an evolution equation exists and closes at the level of the energy spectrum averaged over our ensemble of realizations. If this equation were explicitly available, a typical numerical integration scheme would start by discretizing it in (for example through finite differences, FD, or finite elements, FE) and then proceed to the derivation of ordinary differential equations (ODEs) in time for the vector of values at the FD discretization points, or of the FE coefficients. These ODEs would be of the general form ; the temporal evolution would be obtained using established numerical integration techniques, of which the simplest is arguably the explicit forward Euler method
(8) 
Clearly, implementing this, as well as other initial value solvers, requires an evaluation of the time derivative (the right hand side of the ODE set) at prescribed system states. In projective integration gk1 (); gk2 () this time derivative is not obtained through such a function evaluation of , since is not available in closed form; instead, it is obtained from processing the results of short bursts of simulation: starting at time , we integrate the original equation for a short time with a time step to obtain a sequence , where . We can then use, for example, the last two points of this sequence to estimate the time derivative, and predict as
(9) 
better estimators (least squares, maximum likelihood) of the time deerivative can clearly be used as well. The procedure provides an approximation of the “projected in time” state , which can then be used to restart the computation, and the process repeats.
When the only available tool is a “black box” simulator of the system evolution equation, projective integration may significantly accelerate computations for systems with gaps in their eigenvalue spectra (separation of time scales). However, in our case, the equation we want to projectively integrate is an equation for energy spectrum evolution, while the available simulator is a simulator for velocity fields; we want to projectively integrate the coarse grained behavior (spectra) through observations of the fine scale behavior (velocities). This now becomes coarse projective integration. From our simulation we have available the full state (the velocity field ensembles) at time ; full direct simulation will then provide full states at each of the intermediate time steps. We restrict the full state to the coarse observables (the discretization of the spectrum); we use these observations to estimate the time derivative of these coarse observables; and we finally pass the time derivative estimates to our integrator of choice (here a simple forward Euler) in order to produce our approximation of the coarse observables at the next (projected) point in time .
Here comes a vital step in our equationfree coarse graining approach: our “black box” simulator requires a fine scale initialization (detailed velocities); yet we only have coarse grained initial conditions (ensembleaveraged spectrum discretizations). For the procedure to continue, we need to somehow construct full finescale states consistent with the desired coarse initial conditions. This construction is the socalled lifting step pnas (); gkt02 (); manifesto (); manifesto_short (). The missing degrees of freedom in our construction are the phases of the velocity fields. Our assumption that an equation exists and closes in terms of only ensemble averaged spectra implies that these remaining degrees of freedom either do not couple with the evolution of the coarse grained variables we retain, or, if they do (as we expect), they quickly get slaved to the coarse variables. “Quickly” here implies a comparison with the time scales of the evolution of the retained coarse variables themselves. Exploiting this “quick slaving” (which is implicit in the assumption that an equation for the retained coarse variables exists and closes) is the main point of equationfree computation: since the fine scale code is available, running it for a short time, even with wrongly initialized “additional variables”, should provide the required slaving on demand. This can be considered a subcase of optimal prediction chorin (). Given an energy spectrum initial condition, we therefore initialize the velocity fields as follows: the amplitudes of the velocity Fourier modes are chosen consistently with the prescribed initial energy spectrum, while the phases are chosen randomly. We then let our direct simulator evolve for some initial “healing” interval (described below) before we start observing the solution energy spectrum, in order to estimate its time derivative; this initial period is allowed so that our “wrong” phase initialization becomes “healed”, that is, so that our phases or their statistics become slaved to the evolution of the spectra. We monitor this slaving through observations of the structure function defined earlier.
Here is a concise description of our coarse (spectrum level) projective (forward Euler) integration. The procedure starts after we have already simulated (starting with zero initial conditions) for a period of approximately 3000 time steps; during this period the “high–” end of the spectrum approaches stationarity.

Starting with a “current” velocity field initial condition, we evolve the randomly forced Burgers simulator for a short interval , and restrict the velocity fields at successive instances within this interval to ensembleaveraged energy spectra. Fig. 5(a) shows the restricted spectra at times 3500, 4500, 5500, and 7000.

Using the procedure described in section II, we compute the stretching factors along various straight lines passing through the fulcrum and intersecting these spectra. Our discretization of the spectrum does not come in the form for a number of mesh points (or ), but rather in the form for a number of “discretization slopes” ; the ensembleaveraged spectrum evolution appears smooth enough, so that approximately 20 such discretization angles and simple interpolation is sufficient. We now compute the rate of change of the stretching factors in logarithmic time, and use the same to project the energy spectrum to a future time . Here, the stretching factors are computed along straight lines of slopes , and their rate of change in logarithmic time is computed using a least squares fit as
(10) where and are the coefficients of the least squares fit. We then use equation (10) to project the stretching factor to a time , and compute the projected spectrum using
(11) An estimate of the spectrum at using those at times 3500, 4500, 5500, and 7000 is shown in Fig. 5(a).

We now lift the projected spectrum to an ensemble of consistent velocity fields in order to repeat the process. This is done by randomizing the phases for all the wavenumbers, while retaining the amplitudes corresponding to the projected spectrum. More precisely, we generate 16 different initial velocity fields as , where is a random number between and chosen from a uniform distribution.

We now repeat the process: we evolve the fine scale simulator (the randomly forced Burgers) with the initial conditions thus obtained; but before we start using the resulting spectra to estimate a new coarse time derivative and make a new forward Euler projection, we test the slaving of the injected degrees of freedom (the random phases) to our coarse observables (the spectrum). This slaving is monitored through the evolution of the newly initialized structure function . Fig. 5(b) shows the structure function from a simulation that has been running for some time; if we restrict to energy spectra, and then randomize the phases (that is, lift back to velocity fields) this phase information is lost. When we restart the simulation from the lifted velocity fields, however, we clearly observe that the phases “heal”, that is, acquires its original shape quickly compared to the evolution time scales of the spectrum (in less than 500 time steps, which is much smaller than the projection horizon of time steps). We continue now the “short burst” of fine scale simulation to collect data (restrictions to spectra) for the next projection step; the procedure (steps 1–3) then repeats.
Fig. 5 shows an instantiation of such a coarse projective step. The restriction of an initial simulation used to estimate the coarse time derivative are obtained at 3500, 4500, 5500, and 7000; the projection step (which here is performed in logarithmic time) brings us to ; the projected spectrum (black dashed line) is compared to the result of continuing the full simulation until time (blue). More importantly, the result of lifting from the projected spectrum and continuing the evolution till (black dashed line) is also seen to coincide with the result of a full simulation till (blue). Also, the computational savings in the first step of coarse projective integration can now be quantified. Using coarse projective integration, the randomly forced Burgers equation is simulated for 7000 time steps and then a projection step is used to obtain the spectrum at , thus resulting in a saving of 23000 time steps, or a speedup factor of 4.3. The factor is actually a little lower because of the computational effort required for the projection step. Since the dynamics evolve in logarithmic time, this factor will increase for the later projective integration steps.
The choice of performing projective integration (in effect, Taylor series expansion of the coarse solution in the independent variable) in logarithmic rather than linear time is motivated by the observation that our transient is close to a dynamically selfsimilar solution. A detailed discussion of the usefulness of of (coarse) projective integration in problems with continuous symmetries, such as scale invariance, having (possible approximately) selfsimilar solutions can be found in mihalis ().
Iv Coarse Dynamic Renormalization
Direct simulations starting at showed that the spectrum quickly evolves to a “steady state” for ; after this initial period, its low wavenumber portion acquires a “constant shape” which is then observed, over time, to stretch gradually and rather smoothly towards . This stretching becomes increasingly slower in linear time. In what follows we do not concern ourselves with the very fast initial transient (the establishment of a steady state for ); nor will we study the very long term dynamics that occur when the “corner” of the spectrum finally approaches . We focus on the “intermediate time” regime during which the spectrum acquires its “constant shape” which then appears to travel (in logarithmic and logarithmic time) towards . This latter observation suggests an even simpler caricature of the selfsimilar evolution; see Fig. 6(a): the spectrum to the left of acquires quickly the appropriate “corner like” shape, which we approximate by two straight lines meeting at a point with wavenumber . The energy of the wavenumbers to the right of remains more or less unchanged, anchored at the fulcrum point. The portion of the spectrum to the left of subsequently simply moves towards lower with constant shape. A constant shift of this part of the spectrum in corresponds to a stretching by a constant factor in , and clearly suggests that for our apparently selfsimilar solution
(12) 
and an appropriate choice of constants and . For the data shown in Fig. 4, initialized with a zero spectrum at and taken at times , the values of the exponent and the constant are computationally found to be and . These values can be extracted from Fig. 6(b), which shows plots of the wavenumbers versus and versus on a logarithmic scale.
This appears consistent with the numerical simulations of Yakhot and She vyPRL88 () as well as with dynamic renormalization group results suggesting a nontrivial dynamic scaling of the temporal frequency (or ) as opposed to or . The same exponent is obtained using different spectrum “caricatures” based on the same lowwavenumber straightline approximations (not shown here). We would like to stress that computation of dynamic scaling exponents is an extremely difficult and computationally expensive task. The fact that the nontrivial exponent (4) emerged from our relatively fast calculation seems remarkable. Moreover, it serves as an independent test of our coarse graining procedure.
Given these observations, it is clear that the intermediate time behavior can be recovered in an even more economical fashion than coarse projective integration: it is enough to find the “right shape” (the slope of the leftmost part of the spectrum) as well as the correct stretching rate (encoded in the exponent); this information suffices to approximate the spectra at all intermediate times. This is analogous to computing the shape and the (constant) speed of a traveling wave in a problem with translational invariance. Knowing the shape and speed we can construct the solution at any future time; scale invariance here takes the place of translational invariance for traveling, and the stretching rate (quantified by the similarity exponents) is analogous to the traveling speed. The problem therefore reduces to finding the “right” spectrum shape (for ): a spectrum which, when evolved for some time to and rescaled in –space, remains unchanged. This gives rise to the fixed point problem:
(13) 
Solving such coarse fixed point problems using coarse time steppers has been discussed and illustrated in several contexts chenJNNFM (); zouPRE (); chenPRE (); kesslerPRE (). The simplest fixed point algorithm is successive substitution: we start with an initial guess of the selfsimilar spectrum shape which we evolve for some time using the direct randomly forced Burgers simulator, including ensemble averaging. We then rescale the resulting spectrum by a constant factor, and repeat the process; for a discussion of a socalled template based approach to choosing the appropriate scaling factor at each iteration, see aronson01 (); rowley03 (), and our choice for this problem will be discussed below. If the selfsimilar spectrum is stable (as is the case here), this iteration converges to its scaleinvariant shape. We remark, before proceeding, that other fixed point algorithms (in particular, matrixfree Newton–Krylov–GMRES, kelley04 ()) have been used for solving this type of problem; such Newton–based algorithms can converge to selfsimilar solutions even if they are not stable. In what follows, we demonstrate that successive substitution using a dynamically rescaled coarse time stepper does indeed converge to the right (invariant) spectrum shape.
Long direct simulation will also eventually converge to this invariant shape; yet, while approaching the right shape, the transient will also move towards lower wavenumbers (will “stretch”) and this makes the evolution increasingly slower. We therefore have the potential to realize a computational advantage by rescaling the spectrum after partial evolution to a constant reference scale, a scale in whose neighborhood the dynamics are as fast as possible in physical time; this allows us to observe the approach to the stable invariant shape without this shape stretching.
We now consider another caricature of the (ensemble averaged) spectra during our intermediate–time simulation, approximated by two straight lines in the region as shown in Fig. 7(a). These two straight lines intersect at a “corner”; we call this corner wavenumber . As we discussed earlier in this section, we observe that the slope of the straight line approximating the region of the spectrum to the left of the corner remains (approximately) constant () during selfsimilar evolution; finding the right shape becomes then equivalent to finding this slope. In order to obtain the selfsimilar shape, we initialize our spectrum as follows: The portion of the spectrum to the right of the fulcrum is initialized at its (already) steady shape. For , a straight line approximation is used, the slope of which is computed by a local fitting to the spectrum obtained from the initial direct simulation. Finally, for the region , we initialize the spectrum with a straight line of a “guess” slope . We now pin our reference scale by selecting a particular, fixed ; this means that our procedure will converge on the representative of the oneparameter family of selfsimilar shapes that has its corner fixed there.
Our coarse renormalized time stepper involves lifting from this guess spectrum to full Burgers velocity fields, (ensemble) evolution of these fields by the direct randomly forced Burgers simulator for a time , restriction back to spectra, and ensemble averaging, after this time. A further restriction step is now available by our choice of caricature, in the form of approximating the spectrum by two piecewise linear segments. This gives us a new corner and a new slope . We now rescale the resulting spectrum by shifting back to the old but keeping the new ; this allows us to persistently observe the system at a convenient scale. Iterations of this “liftevolverestrictrescale” procedure is shown in Fig. 7(a); Fig. 7(b) shows a sequence of several successive renormalization substitutions converging to the invariant shape; here it takes approximately six iterations to convergence if the originally guessed slope is less steep than the true one; steeper initial slopes practically converge within one iteration.
From the last iteration (upon convergence) we obtain the shape of the leftmost part of the spectrum (so ; see forster77 ()). We can also extract the scaling exponent by doing a similar least squares fit to a plot of versus . The constants in the expression (12) are and (here we started iterating at of the original simulation). The larger we can make the corner (the closer to we can practically choose it), the fastest the approach to the selfsimilar shape (the dynamics slow down significantly at progressively lower wavenumbers).
V Summary and Discussion
Using a particular type of randomly forced Burgers equation in one spatial dimension as our illustrative example, and based on direct simulations, we postulated that the system dynamics could, under certain conditions, be coarse grained to an effective evolution equation for the system energy spectrum. The direct simulations also strongly suggest an effectively selfsimilar evolution of this energy spectrum; it is important to reiterate that this ansatz is only made for a particular (intermediate) wavenumber regime, and only for initial conditions that contain no appreciable energy in the large scales. Based on these observations we demonstrated the implementation and use of the “equationfree” computational methodology: the design of short bursts of appropriately initialized computational experiments with the direct simulator that accelerate the study of the system behavior – as compared to full direct simulation only. In particular, we demonstrated coarse projective integration, as well as coarse dynamic renormalization.
Our computations are predicated upon advance knowledge of the right set of coarse observables; these are the quantities in terms of which the (unavailable) evolution equation can in principle be written. Our coarse observables here consisted of (a discretization of) the (ensemble averaged) energy spectrum; the phase variables were assumed (and observed, through the use of the thirdorder structure function ) to become relatively quickly slaved to the energy spectrum dynamics. Algorithms that attempt to detect such observables from datamining simulation results, such as the diffusion map approach coifman (); nadler (), may prove valuable in extending the applicability of the computational tools we presented here. The computational construction of synthetic turbulence fields consistent with such observables (see, for example, the recent paper meneveau06 ()) can be considered as a type of “lifting” scheme in our framework. Beyond the design of (numerical) experiments to accelerate system modeling, it is also possible to use such experiments to answer questions about the nature of the unavailable equation (for example, whether it is a conservation law). Our approach made no assumption about the local or nonlocal nature of the (explicitly unavailable) evolution equation: it may be a partial differential equation or a nonlocal, integrodifferential one. It would be particularly interesting to devise computational experiments with the direct simulation code to explore whether the righthandside of our unavailable evolution equation appears to be local or not in –wavenumber space.
The method proposed here may prove useful for the efficient extraction of dynamic scaling exponents in models of various physical phenomena for which we do not know the governing equations in closed form.
Acknowledgements. This work was partially supported by DARPA and the US DOE(CMPD); the use of computer resources at the Boston University Scientific Computing and Visualization Group (SCV) is gratefully acknowledged.
References
 (1) D. G. Aronson, S. I. Betelu, and I. G. Kevrekidis, Going with the flow: a Lagrangian approach to selfsimilar dynamics and its consequences, http://arxiv.org/abs/nlin.AO/0111055, 2001.
 (2) J. Bec and U. Frisch, Probability distribution functions of derivatives and increments for decaying Burgers turbulence, Phys. Rev. E, 61(2), 13951402, 2000.
 (3) A. Chekhlov and V. Yakhot, Kolmogorov turbulence in a randomforcedriven Burgers equation, Phys. Rev. E, 51(4), 273942, 1995.
 (4) A. Chekhlov and V. Yakhot, Kolmogorov turbulence in a randomforcedriven Burgers equation: Anomalous scaling and probability density functions, Phys. Rev. E, 52(5), 56814, 1995.
 (5) L. Chen, P. G. Debenedetti, C. W. Gear, and I. G. Kevrekidis, From molecular dynamics to coarse selfsimilar solutions: a simple example using equationfree computation, J. Non Newt. Fluid Mech., 120, 215223, 2004.
 (6) L. Chen, I. G. Kevrekidis, and P. G. Kevrekidis, Equationfree dynamic renormalization in a glassy compaction model, Phys. Rev. E, 74, 016702, 2006; condmat/0412773 at arXiv.org.
 (7) A. J. Chorin, O. H. Hald, and R. Kupferman, Optimal prediction and the MoriZwanzig representation of irreversible processes, Proc. Nat. Acad. Sci., 97(7), 296873, 2000.
 (8) R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner and S. Zucker, Geometric diffusion as a tool for harmonic analysis and structure definition of data: diffusion maps, Proc. Nat. Acad. Sci., 102(21), 742631, 2005; Geometric diffusion as a tool for harmonic analysis and structure definition of data: multiscale methods, Proc. Nat. Acad. Sci., 102(21), 742631, 2005.
 (9) W. E, K. Khanin, A. Mazel, and Y. Sinai, Probability distribution functions for the random forced Burgers equation, Phys. Rev. Lett., 78(10), 19047, 1997.
 (10) W. E and E. Vanden Eijnden, Asymptotic theory for the probability density functions in Burgers turbulence, Phys. Rev. Lett., 83(13), 25725, 1999.
 (11) D. Forster, D. Nelson, and M. Stephen, Largedistance and longtime properties of a randomly stirred fluid, Phys. Rev. A, 16(2), 732749, 1977.
 (12) U. Frisch, Turbulence, Cambridge University Press, January 1996.
 (13) C. W. Gear and I. G. Kevrekidis, Projective methods for stiff differential equations: problems with gaps in their eigenvalue spectrum, SIAM J. Sci. Comp., 24(4), 10911106, 2003.
 (14) C. W. Gear and I. G. Kevrekidis, Telescopic projective integrators for stiff differential equations, J. Comp. Phys., 187(1), 95109, 2003.
 (15) C. W. Gear and I. G. Kevrekidis, Constraintdefined manifolds: a legacy code approach to lowdimensional computation, J. Sci. Comp., 25(1), 1728, 2005.
 (16) C. W. Gear, I. G. Kevrekidis, and C. Theodoropoulos, “Coarse” integration/bifurcation analysis via microscopic simulators: microGalerkin methods, Comp. Chem. Eng., 26, 941963, 2002.
 (17) M. Kardar, G. Parisi, and Y.C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett., 58(9), 882892, 1986.
 (18) M. E. Kavousanakis, R. Erban, A. G. Boudouvis, C. W. Gear, and I. G. Kevrekidis, Projective and coarse projective integration for problems with continuous symmetries, in press, J. Comp. Phy., 2006; also available at http://arxiv.org/math.DS/0608122.
 (19) C.T. Kelley, I.G. Kevrekidis, and L. Qiao, NewtonKrylov solvers for timesteppers, available as math.DS/0404374 at arxiv.org, 2004.
 (20) D. A. Kessler, I. G. Kevrekidis, and L. Chen, Equationfree dynamic renormalization of a KPZtype Equation, Phys. Rev. E., 73, 036703(6pp), 2006; available as condmat/0505758 at arxiv.org.
 (21) I. G. Kevrekidis, C. W. Gear, and G. Hummer, Equationfree: the computerassisted analysis of complex, multiscale systems, A.I.Ch.E Journal, 50(7), 13461354, 2004.
 (22) I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg and K. Theodoropoulos, Equationfree coarsegrained multiscale computation: enabling microscopic simulators to perform systemlevel tasks, Comm. Math. Sciences, 1(4), 715762, 2003; also at http://arXiv.org/physics/0209043.
 (23) D. Mitra, J. Bec, R. Pandit, and U. Frisch, Is multiscaling an artifact in the stochastically forced Burgers equation?, http://arxiv.org/abs/nlin.CD/0406049, 2005.
 (24) B. Nadler, S. Lafon, R. C. Coifman and I. G. Kevrekidis, Diffusion maps, spectral clustering and the reaction coordinates of dynamical systems, Appl. Comp. Harm. Anal., 21, 113127, 2006; also at http://arxiv.org/math/0503445.
 (25) A. Polyakov, Turbulence without pressure, Phys. Rev. E, 52(6), 618388, 1995.
 (26) C. Rosales and C. Meneveau, A minimal multiscale Lagrangian map approach to synthesize nonGaussian turbulent vector fields, Physics of Fluids, 18, 075104, 2006.
 (27) C. W. Rowley, I. G. Kevrekidis, J. E. Marsden, and K. Lust, Reduction and reconstruction for selfsimilar dynamical systems, Nonlinearity, 16, 12571275, 2003.
 (28) L. Smith and V. Yakhot, Bose condensation and smallscale structure generation in a random force driven 2D turbulence, Phy. Rev. Letters, 71, 352355, 1993.
 (29) C. Theodoropoulos, YH. Qian, and I. G. Kevrekidis, “Coarse” stability and bifurcation analysis using timesteppers: A reactiondiffusion example, Proc. Nat. Acad. Sci., 97, 98409843, 2000.
 (30) H. W. Wyld, Formulation of the theory of turbulence in an incompressible fluid, Ann. Phys., 14, 143165, 1961.
 (31) V. Yakhot, Meanfield approximation and a small parameter in turbulence theory, Phys. Rev. E, 63(2), 026307, 2001.
 (32) V. Yakhot and A. Chekhlov, Algebraic tails of probability density functions in the randomforcedriven Burgers turbulence, Phys. Rev. Lett., 77(15), 3118  3121, 1996.
 (33) V. Yakhot and Z.S. She, Longtime, largescale properties of the randomforcedriven Burgers equation, Phys. Rev. Lett., 60(18), 184043, 1988.
 (34) Y. Zou, I. G. Kevrekidis, and R. Ghanem, Equationfree dynamic renormalization: selfsimilarity in multidimensional particle system dynamics, Phys. Rev. E, 72, 046702 (7pp.), 2005; also at http://arxiv.org/math.DS/0505358.