Convergence of equation-free methods in the case of finite time scale separation with application to deterministic and stochastic systems

TL;DR

This paper proves the convergence of equation-free methods for systems with finite time-scale separation, including stochastic systems, justifying their use in high-dimensional bifurcation analysis without requiring large scale separation.

Contribution

It provides the first convergence proof for equation-free methods at finite time-scale separation, applicable to both deterministic and stochastic systems, with sharp error estimates.

Findings

Convergence of equation-free methods is proven for finite healing time.

The results apply to systems with normal hyperbolicity without large time-scale separation.

Demonstration with Michaelis-Menten kinetics confirms sharpness of error estimates.

Abstract

A common approach to studying high-dimensional systems with emergent low-dimensional behavior is based on lift-evolve-restrict maps (called equation-free methods): first, a user-defined lifting operator maps a set of low-dimensional coordinates into the high-dimensional phase space, then the high-dimensional (microscopic) evolution is applied for some time, and finally a user-defined restriction operator maps down into a low-dimensional space again. We prove convergence of equation-free methods for finite time-scale separation with respect to a method parameter, the so-called healing time. Our convergence result justifies equation-free methods as a tool for performing high-level tasks such as bifurcation analysis on high-dimensional systems. More precisely, if the high-dimensional system has an attracting invariant manifold with smaller expansion and attraction rates in the tangential direction than in the transversal direction (normal hyperbolicity), and restriction and lifting satisfy some generic transversality conditions, then an implicit formulation of the lift-evolve-restrict procedure generates an approximate map that converges to the flow on the invariant manifold for healing time going to infinity. In contrast to all previous results, our result does not require the time scale separation to be large. A demonstration with Michaelis-Menten kinetics shows that the error estimates of our theorem are sharp. The ability to achieve convergence even for finite time scale separation is especially important for applications involving stochastic systems, where the evolution occurs at the level of distributions, governed by the Fokker-Planck equation. In these applications the spectral gap is typically finite. We investigate a low-dimensional stochastic differential equation where the ratio between the decay rates of fast and slow variables is 2.

Equations

$$P:[0,\infty )\times \mathrm{dom}L\ni (t,x_L)\mapsto R(M(t;L(x_L)))\in \mathrm{rg}R\text{mbox}$$

$$P(t_{\mathrm{skip}};y)=P(t_{\mathrm{skip}}+\delta ;x)$$

$${\Phi }_{\ast }(\delta ;\cdot )$$

$${\Phi }_{\ast }:$$

$$R(g(L(y_{\ast })))$$

$$y_{\ast }={\Phi }_{\ast }(\delta ;x)=(g\circ L{)}^{-1}\circ M(t_{\mathrm{skip}};\cdot {)}^{-1}\circ M(\delta +t_{\mathrm{skip}};\cdot )\circ g\circ L(x)\text{mbox}.$$

$$R\circ M(t_{\mathrm{skip}};\cdot )\circ g\circ L(y_{\ast })=R\circ M(\delta +t_{\mathrm{skip}};\cdot )\circ g\circ L(x)$$

$$R(M(t_{\mathrm{skip}};L(y_{t_{\mathrm{skip}}})))=R(M(\delta +t_{\mathrm{skip}};L(x)))\text{mbox}.$$

$$\partial^{j}y_{t_{\mathrm{skip}}}-\partial^{j}y_{*}\sim\exp(((2j+1)d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})\quad\mbox{for ${t_{\mathrm{skip}}}\to\infty$}$$

$$\dot{u}=f(u),u\in {\mathbb{R}}^D\text{mbox},$$

$$M:\mathbb{R}\times {\mathbb{R}}^D\to {\mathbb{R}}^D\text{mbox},(t;u)\mapsto M(t;u)$$

$$\Vert {\partial }_2^{j}M(t;u)[v_1,\dots ,v_{k_j}]\Vert$$

$$\Vert {\partial }_2^{j}M(t;u)-{\partial }_2^{j}M(t;g(u))\Vert$$

$$R$$

$$L$$

$$L(\mathrm{dom}L)$$

$$\mathrm{rank}\frac{\partial }{\partial x}\left[g(L(x))\right]=\mathrm{rank}\left[\partial g(L(x))\circ \partial L(x)\right]=d\text{mbox}.$$

$$\dim \partial R(u)\mathcal{N}(u)=d\text{mbox}.$$

$$g\circ L$$

$${\Phi }_{\ast }$$

$$R(g(L(y)))$$

$$P_{\ast }$$

$$y={\Phi }_{\ast }(\delta ;x)\text{mbox}if{P}_{\ast }(0;y)=P_{\ast }(\delta ;x)\text{mbox}.$$

$$y={\Phi }_{\ast }(\delta ;x)\text{mbox}if{P}_{\ast }(t_{\mathrm{skip}};y)=P_{\ast }(t_{\mathrm{skip}}+\delta ;x)\text{mbox}.$$

$$P$$

$$\displaystyle\Phi_{t_{\mathrm{skip}}}:\mathbb{R}\times\operatorname{dom}\operatorname{\mathcal{L}}\ni(\delta,x)\mapsto y\in\operatorname{dom}\operatorname{\mathcal{L}}\mbox{,\ where $y$ solves\ }P({t_{\mathrm{skip}}};y)=P({t_{\mathrm{skip}}}+\delta;x)\mbox{}$$

$$\operatorname{dist}(M(t;g(\operatorname{\mathcal{L}}(x))),\partial{\cal C})\geq c_{\partial}\mbox{\quad for all $t\geq-\delta_{\max}$ and some given $c_{\partial}>0$.}$$

$$\Vert {\partial }_2^j{\Phi }_{t_{\mathrm{skip}}}(\delta ;x)-{\partial }_2^j{\Phi }_{\ast }(\delta ;x)\Vert \le C\exp (((2j+1)d_{\tan }-d_{\mathrm{tr}})t_{\mathrm{skip}})$$

$$R(M(t_{\mathrm{skip}};L(y)))$$

$$R(M(t_{\mathrm{skip}};g(L(y_{\ast }))))$$

$$P_{\ast }(t_{\mathrm{skip}};y)=P_{\ast }(t_{\mathrm{skip}};y_{\ast })+[P_{\ast }(t_{\mathrm{skip}};y)-P(t_{\mathrm{skip}};y)]+[P(t_{\mathrm{skip}}+\delta ;x)-P_{\ast }(t_{\mathrm{skip}}+\delta ;x)]\text{mbox}.$$

$$\partial P_{\ast }(y)\left[{\partial }^{j}y-{\partial }^j{y}_{\ast }\right]=\left[\partial P_{\ast }(y)-\partial P_{\ast }(y_{\ast })\right]{\partial }^j{y}_{\ast }+r\text{mbox}.$$

Full text

More precisely, if the high-dimensional system has an attracting invariant manifold with smaller expansion and attraction rates in the tangential direction than in the transversal direction (normal hyperbolicity), and restriction and lifting satisfy some generic transversality conditions, then an implicit formulation of the lift-evolve-restrict procedure generates an approximate map that converges to the flow on the invariant manifold for healing time going to infinity. In contrast to all previous results, our result does not require the time scale separation to be large. A demonstration with Michaelis-Menten kinetics shows that the error estimates of our theorem are sharp.

The ability to achieve convergence even for finite time scale separation is especially important for applications involving stochastic systems, where the evolution occurs at the level of distributions, governed by the Fokker-Planck equation. In these applications the spectral gap is typically finite. We investigate a low-dimensional stochastic differential equation where the ratio between the decay rates of fast and slow variables is $2$.

{AMS}

65Pxx, 37Mxx, 34E13

1 Introduction

High-dimensional dynamical systems with time scale separation have under certain assumptions the potential to be studied and understood through a reduction to low-dimensional systems. In most cases these reduction methods are applied directly to the high-dimensional systems of equations [39]. The most common approaches are referred to as averaging and mean-field approximation [44], the slaving principle or adiabatic elimination [19] in the physics literature. The aim of these methods is to reduce the complexity of a high-dimensional (here also called microscopic) system to a relatively simple low-dimensional (here also called macroscopic) system. After reduction, the long-term dynamics of the system can be analyzed by studying the low-dimensional macroscopic system, using techniques that may only be available for low-dimensional deterministic systems (e.g., detailed bifurcation analysis). The underlying assumption is that a trajectory of the microscopic system will rapidly relax onto a low-dimensional manifold, which it will then track on a longer time scale, following the slower macroscopic equations. Thus, one speaks of slow variables, which are the coordinates on the slow low-dimensional manifold, and fast variables transversal to the slow manifold. The notion that the fast variables are “slaved” by the slow variables describes that over long time the microscopic trajectories track the slow manifold.

The justification for this reduction is simplest and strongest if the underlying microscopic dynamical system possesses a low-dimensional attracting invariant manifold. In these cases mathematical theorems on persistence of invariant manifolds can be applied. Proofs were given by Fenichel [15] and Hirsch et al. [20] for finite-dimensional smooth dynamical systems such as ordinary differential equation (ODEs) and maps and by Bates et al. [6, 7] for general semiflows (covering certain classes of partial differential equations). Certain cases of averaging (such as periodic and quasi-periodic forcing) may also be reduced to invariant manifold persistence.

The case for model reduction is more subtle if the microscopic system is stochastic (or more generally, ergodic), for example, if the model is given by a multi-particle or agent-based simulation. The time-scale separation for these systems occurs if, for example, the number of particles is large. A model case for stochastic systems is the reduction of a high-dimensional system of stochastic differential equations (SDE) to a low-dimensional SDE acting on a slower time scale (smaller drift terms and smaller noise amplitude than the microscopic system). In this case, the arguments for model reduction look formally similar to the case of attracting manifolds in deterministic systems [36, 17]. However, the underlying mathematical convergence results are not as strong. Two aspects in which the deterministic results are stronger than the stochastic results will have implications on convergence results for computational methods:

Validity for finite time-scale separation: For invariant manifolds in a deterministic ODE the time scale separation (let us call it $\varepsilon$) is measured as the ratio between the rate of attraction along directions tangential to the manifold ($d_{\mathrm{tan}}$) and transversal to the manifold ($d_{\mathrm{tr}}$, so $\varepsilon=d_{\mathrm{tan}}/d_{\mathrm{tr}}$). As long as this ratio $\varepsilon$ is less than unity, the manifold persists. Let us call this persistent low-dimensional manifold ${\cal C}$. Persistence implies that, even for a finite $\varepsilon$, a reduced model on this manifold ${\cal C}$ exists, describing some trajectories of the microscopic system with perfect accuracy (those that lie on ${\cal C}$). In practice the dynamics on the slow manifold ${\cal C}$ is often approximated by an expansion in $\varepsilon$.

Shadowing: Even more, every point $u$ from an open neighborhood of ${\cal C}$ has a shadowing point $g(u)\in{\cal C}$. The difference between the trajectories starting from $u$ and $g(u)$ goes to zero in time with a rate close to $d_{\mathrm{tr}}$. This means that the reduced model describes all nearby trajectories even for positive $\varepsilon$ with perfect accuracy except for rapidly decaying terms. The nonlinear projection $u\mapsto g(u)\in{\cal C}$ is called the stable fiber projection.

Compared to the above, the precise mathematical convergence statements in [36, 17] for stochastic systems with time scale separation are weaker. They are concerned with the limit $\varepsilon\to 0$ and prove that moments of the slow coordinates of the microscopic trajectory and of the trajectories of the slow model, derived by a formal expansion in $\varepsilon$, converge to each other for $\varepsilon\to 0$ [36].

On-demand computation of slow flow — Equation-free

framework

The above mathematical theorems underpin the derivation of approximate low-dimensional models for large-dimensional systems. However, they also provide guidance for the convergence analysis of computational methods that avoid the explicit derivation of a low-dimensional model, but merely assume its existence. A general framework for analysing slow-time scale behaviour of systems with time scale separation was proposed by Kevrekidis et al. under the name “equation-free computations” [23, 16, 25]. The assumption behind equation-free computations is the existence of a slow low-dimensional description (in $\mathbb{R}^{d}$) for some macroscopic quantities of the high-dimensional microscopic system (which is defined in $\mathbb{R}^{D}$) that contains a $d$-dimensional invariant manifold ${\cal C}$, which we will call the slow manifold. We do not append a subscript $\varepsilon$ to ${\cal C}$, since our main result will only assume existence and smoothness of the invariant manifold ${\cal C}$ and its stable fiber projection, but not consider the limit of time scale separation $\varepsilon\to 0$. The framework, illustrated in Fig. 1, only relies on the availability of a microscopic time stepper (a map $M(t;\cdot):\mathbb{R}^{D}\mapsto\mathbb{R}^{D}$ for $t\geq 0$) that can be called at selected microscopic initial values $u\in\mathbb{R}^{D}$. The goal is to compose a macroscopic time stepper $\Phi_{*}(\delta;\cdot):\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ for $\delta\in\mathbb{R}$ (possibly including $\delta<0$) in some coordinates for the slow manifold ${\cal C}$, which is then amenable to higher-level tasks such as bifurcation analysis.

For equation-free computations the user also has to choose two operators, the lifting $\operatorname{\mathcal{L}}:\mathbb{R}^{d}\mapsto\mathbb{R}^{D}$ and the restriction $\operatorname{\mathcal{R}}:\mathbb{R}^{D}\mapsto\mathbb{R}^{d}$, which are maps between the original high-dimensional ($\mathbb{R}^{D}$) microscopic level and the low-dimensional ($\mathbb{R}^{d}$) macroscopic level. The user-defined lifting ${\cal L}$ and restriction ${\cal R}$, together with the time stepper $M(t;\cdot)$ define the central building block of equation-free methodology, the “lift-evolve-restrict” map,

[TABLE]

(see Fig. 1). For a given value $x_{L}\in\mathbb{R}^{d}$ of macroscopic quantities, one first applies the lifting $\operatorname{\mathcal{L}}$ to $x_{L}$ getting a microscopic state $u$, then one runs the microscopic simulation for time $t$ starting from $u$ (applying the microscopic evolution $M(t;u)$), and finally one applies the restriction $\operatorname{\mathcal{R}}$ to the result $M(t;u)$.

The use of the lift-evolve-restrict map $P$ assumes that the trajectory $t\mapsto M(t;u)$ of the time stepper will be close to the slow manifold ${\cal C}$ most of the time. Assuming this, equation-free methods aim to extract information about the slow flow along ${\cal C}$ by calling the lift-evolve-restrict map $P$ judiciously. The simplest approach would be to use $P(t;\cdot)$ as an approximation for $M(t;\cdot)$ restricted to ${\cal C}$ (called explicit equation-free computation in [33]).

When using equation-free methods one faces several challenges, both analytical and in terms of implementation. First, as the slow manifold ${\cal C}$ cannot be assumed to be known to the user, the method cannot assume that the user provides a lifting operator $\operatorname{\mathcal{L}}$ that maps onto ${\cal C}$. This leads to initial fast transients in the trajectory that will also change the supposedly slow variables, unless the stable fiber projection $g:\mathbb{R}^{D}\to{\cal C}$ keeps the restriction constant (the criterion would be $\operatorname{\mathcal{R}}\circ g\circ\operatorname{\mathcal{L}}\approx I$). Since the projection $g$ cannot be assumed to be known either, this implies that an unknown nonlinear transformation is applied to the variables in $\operatorname{dom}\operatorname{\mathcal{L}}$ before the slow dynamics start. A detailed illustration of this problem is given in Fig. 2 and its description in Section 2.

Second, the justification for equation-free methods relies on the stronger results for classical attracting invariant manifolds of deterministic systems (including persistence of the slow manifold for finite time-scale separation and its shadowing properties via stable fiber projection). However, the methods are commonly applied to stochastic or deterministic chaotic systems with time scale separation, for which convergence results are weaker. In the stochastic case the microscopic time stepper $M(t;\cdot)$ applies to densities not single trajectories. Finally, for applications with stochastic microscopic systems the additional difficulty of low computational accuracy in the evaluation of $M(t;\cdot)$ and possibly $\operatorname{\mathcal{L}}$ and $\operatorname{\mathcal{R}}$ may impose practical limitations.

This paper addresses the first challenge, the unknown slow manifold and fiber projection. It proves convergence of the implicit approximation $y=\Phi_{t_{\mathrm{skip}}}(\delta;x)$ for the slow flow $\Phi_{*}(\delta;x)$, given by the solution $y$ of the $d$-dimensional nonlinear system

[TABLE]

for sufficiently large healing time ${t_{\mathrm{skip}}}$ and a fixed finite time scale separation for the scenario of an attracting $d$-dimensional invariant manifold in $\mathbb{R}^{D}$ (strong reduction results are available in this scenario). We also give a demonstration how the implicit equation-free formulation behaves when it is applied to moments of distributions in a stochastic system. Starting from this demonstration, we outline in our subsequent discussion how convergence statements for stochastic systems may have to be formulated.

Applications and recent practical improvements

A motivation for using the equation-free framework is that it extends methods which are otherwise only applicable to low-dimensional dynamical systems directly to simulations of high-dimensional complex systems. Classical applications of equation-free methods were macroscopic bifurcation analysis for microscopic simulations in chemical engineering (see [24] for a review). Recently similar analysis was performed on stochastic network models of neurons [4, 30] or disease spread [18], or on agent-based models in ecology [45] and social sciences (for example, for consumer lock-in [3], for pedestrian flow [33, 32], or for trading [42]). Another example for a high-level task accessible via equation-free methods is control design [43, 42].

Recent modifications and improvements to equation-free methods in multi-particle or agent-based simulations are variance reduction [35, 3], restriction of computations to patches in space [40, 41, 28] (for which a-priori error estimates can be proven [40, 41]), and data-driven selection of the slow variables using diffusion maps [12, 33]. Debrabant et al. [13] construct an acceleration scheme for Monte-Carlo simulations of high-dimensional SDEs based on moments of densities (the macroscopic variables), and prove its convergence as the number of moments goes to infinity.

2 Current state of analysis

Geometry of the idealized case of an attracting slow manifold

Analysis of the equation-free framework (based on lift-evolve-restrict) is still ongoing. Convergence analysis with general a-priori error estimates has been performed mostly for the idealized case where the $D$-dimensional microscopic problem has a $d$-dimensional attracting invariant slow manifold ${\cal C}$, which is rarely encountered in the practical applications listed above. Exceptions are, for example, a study of bursting neurons [8] and the application of implicit equation-free computations to generalize an algorithm for growing stable manifolds of fixed points of two-dimensional maps a delay-differential equation with an unknown two-dimensional slow manifold [38]. Even for this idealized case one faces the geometric difficulty illustrated in Fig. 2.

The geometry shows an example scenario where the microscopic system is two-dimensional, and the slow manifold ${\cal C}$ is horizontal (and, thus, the slow motion is purely horizontal, drifting to the left). Here we choose a lifting $\operatorname{\mathcal{L}}$ that maps also onto a horizontal line $\operatorname{rg}\operatorname{\mathcal{L}}$. However, $\operatorname{rg}\operatorname{\mathcal{L}}$ is at a distance to ${\cal C}$, because the precise location of ${\cal C}$ is in practice unknown. The restriction $\operatorname{\mathcal{R}}$ is the horizontal component of any point $u\in\mathbb{R}^{2}$. The spaces $\operatorname{dom}\operatorname{\mathcal{L}}$ and $\operatorname{rg}\operatorname{\mathcal{R}}$ (both one-dimensional) are drawn separately for clarity in Fig. 2, but they may be identical in examples. The fast motion of $M(t;\cdot)$ is not perfectly vertical, but has a significant horizontal component. Figure 2 also shows how the map $P$ acts on a typical point $x_{L}$, showing its image $\operatorname{\mathcal{L}}(x_{L})$, the result of the evolution, $M(t;\operatorname{\mathcal{L}}(x_{L}))$, and the result of the restriction $P(t;x_{L})=x_{R}=\operatorname{\mathcal{R}}(M(t;\operatorname{\mathcal{L}}(x_{L})))$.

The point $x_{s}\in{\cal C}$ in Fig. 2 is defined as the unique point $x_{s}$ on ${\cal C}$ such that $M(t;\operatorname{\mathcal{L}}(x_{L}))-M(t;x_{s})$ converges at an exponential rate $d_{\mathrm{tr}}$ that is larger than the maximal rate of contraction $d_{\mathrm{tan}}$ tangential to ${\cal C}$ (which is horizontal). As mentioned in the introduction as shadowing, this mapping is defined for every point $u$ in the neighborhood of ${\cal C}$: for every $u$ near ${\cal C}$ there exists a point $g(u)\in{\cal C}$ such that $M(t;u)-M(t;g(u))\sim\exp(-d_{\mathrm{tr}}t)$ (in the illustration $u=\operatorname{\mathcal{L}}(x_{L})$, $g(u)=x_{s}$). This point $g(u)$ is called the stable fiber projection of $u$. The map $g$ is known to have the same regularity as ${\cal C}$ [15, 20]. For $d_{\mathrm{tan}}\ll d_{\mathrm{tr}}$ the map can be expanded in orders of $\varepsilon=d_{\mathrm{tan}}/d_{\mathrm{tr}}$. The thin grey lines (called stable fibers or isochrones) in Fig. 2 indicate how points in the plane are projected onto ${\cal C}$ under the nonlinear projection $g$ for the illustrative example.

Figure 2 makes clear that the dynamics of the map $P(t;x)$ is qualitatively different from the dynamics of $M(t;\cdot)$ restricted to ${\cal C}$, $M(t;\cdot)|_{\cal C}$. For the particular geometry shown in Fig. 2 $P(t;\cdot)$ has a unique stable fixed point if the horizontal attraction/expansion rate $d_{\mathrm{tan}}$ of $M(t;\cdot)$ on ${\cal C}$ is sufficiently small compared to the attraction rate $d_{\mathrm{tr}}$ transversal to ${\cal C}$. This fixed point is nearly independent of the dynamics of $M(t;\cdot)$ on ${\cal C}$.

More generally, if the lifting operator $\operatorname{\mathcal{L}}$ does not map $x_{L}$ into the low-dimensional slow manifold ${\cal C}$ then the initial part of the trajectory $t\mapsto M(t;\operatorname{\mathcal{L}}(x_{L}))$, which is computed as part of the lift-evolve-restrict map $P$, is a rapidly changing transient toward the slow manifold ${\cal C}$, which will generically also change the resulting $x_{R}$.

In the limit of infinite time-scale separation (that is, the derivative of $M$ with respect to time, $\partial_{1}M(t;,u)$, goes to [math] for $u\in{\cal C}$) the dynamics of the lift-evolve-restrict map $P$ is a small perturbation of the map $\operatorname{\mathcal{R}}\circ g\circ\operatorname{\mathcal{L}}$. Unless this limit map equals the identity, $P(t;\cdot)$ cannot be a good approximation of the slow flow along the manifold ${\cal C}$. Using $x$ in the domain of the lifting $\operatorname{\mathcal{L}}$ and the map $g\circ\operatorname{\mathcal{L}}:\operatorname{dom}\operatorname{\mathcal{L}}\mapsto{\cal C}$ onto the manifold ${\cal C}$ as the coordinate map, the slow flow $\Phi_{*}$ has the form

[TABLE]

(using the notation $(\cdot)^{-1}$ for the inverse map). This definition is not directly computable since the nonlinear projection $g$ is unknown in general.

Feasible approaches to construct an accurate approximation of $M(t;\cdot)$ restricted to ${\cal C}$ are constrained runs, as discussed by Gear, Zagaris et al. [16, 48, 49], or the introduction of a healing time ${t_{\mathrm{skip}}}$. The latter approach is studied in this paper.

Constrained runs

The approach of [16, 48, 49] to ensuring that $\operatorname{\mathcal{R}}\circ g\circ\operatorname{\mathcal{L}}$ is close to the identity is to enforce that the lifting $\operatorname{\mathcal{L}}$ maps onto the manifold ${\cal C}$ with sufficient accuracy for all $x$ in its domain. Usually, this requires an additional scheme involving the iterative application of $\operatorname{\mathcal{L}}$ and $M$; see [16, 48, 49]. The a-priori error estimates prove that the lift-evolve-restrict scheme with these additional iterations has an error of order $(d_{\mathrm{tan}}/d_{\mathrm{tr}})^{m}$ if the constrained runs scheme is of order $m$, where $d_{\mathrm{tan}}$ is the attraction/repulsion time scale tangential to the slow manifold ${\cal C}$ and $d_{\mathrm{tr}}$ is the transversal attraction rate. The ratio $d_{\mathrm{tan}}/d_{\mathrm{tr}}$ measures the time scale separation. It is assumed to be small when applying constrained runs (and called $\varepsilon$), and $O(\varepsilon^{m})$ convergence is proven in [16, 48, 49] in the limit $\varepsilon\to 0$. This limit will not be required in our proof, later on.

Implicit formulation with healing time

A second, alternative, approach is to introduce a healing time ${t_{\mathrm{skip}}}$, exploiting that $M$ attracts along the fibers [25, 5].

Marschler et al. [31] show that the healing time ${t_{\mathrm{skip}}}$ can be motivated by introducing an additional shift $M({t_{\mathrm{skip}}};\cdot)$ and its inverse into Eq. 1 (note that $M({t_{\mathrm{skip}}};\cdot)$ is invertible on the slow manifold ${\cal C})$:

[TABLE]

Removing the inverses in Eq. 2 leads to an implicit equation for $y_{*}=\Phi_{*}(\delta;x)$ with the healing time ${t_{\mathrm{skip}}}$ as an additional parameter:

[TABLE]

In Eq. 3 the parameter ${t_{\mathrm{skip}}}$ has no effect since $M({t_{\mathrm{skip}}};\cdot)$ is invertible on the slow manifold. However, the difference $M({t_{\mathrm{skip}}};\cdot)\circ g-M({t_{\mathrm{skip}}};\cdot)$ decreases with ${t_{\mathrm{skip}}}$ (at rate $\sim\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$). In Fig. 3 the distance between points along the trajectory starting from $\operatorname{\mathcal{L}}(x)$ (in red) and their projections $g\circ\operatorname{\mathcal{L}}(x)$ (white) illustrates this convergence. Thus, we may approximate $M({t_{\mathrm{skip}}};\cdot)\circ g$ by $M({t_{\mathrm{skip}}};\cdot)$ in Eq. 3. This results in a computable approximation $y_{t_{\mathrm{skip}}}=\Phi_{t_{\mathrm{skip}}}(\delta;x)$ of $y_{*}$, given implicitly by the equation

[TABLE]

Figure 3 illustrates the effect of increasing healing time ${t_{\mathrm{skip}}}$ in the scenario introduced in Fig. 2. The points $y_{1}$ and $y_{2}$ are the solutions of Eq. 4 for two different healing times $t_{\mathrm{skip},1}<t_{\mathrm{skip},2}$. Equation (4) means that the points $y_{j}$ are defined as those elements of $\operatorname{dom}\operatorname{\mathcal{L}}$ for which the trajectory starting from $\operatorname{\mathcal{L}}(y_{j})$ has the same horizontal component (restriction $\operatorname{\mathcal{R}}$) as $M(t_{\mathrm{skip},j}+\delta;\operatorname{\mathcal{L}}(x))$ after time $t_{\mathrm{skip},j}$.

The implicit approach was analyzed and illustrated in a traffic model in [31] and will also be studied in this paper. Vandekerckhove et al. [46] proposed and demonstrated a similar approach, but applied the healing time backward in time by fixing the image of the restriction: they solve $x=\operatorname{\mathcal{R}}(M({t_{\mathrm{skip}}};\operatorname{\mathcal{L}}(x_{b})))$ for $x_{b}$ first and then set $y=\operatorname{\mathcal{R}}(M(\delta+{t_{\mathrm{skip}}};\operatorname{\mathcal{L}}(x_{b})))$. This gives an (approximate) representation $\Phi_{*}^{\operatorname{\mathcal{R}}}$ of the slow flow in the coordinates on the image of the restriction $\operatorname{\mathcal{R}}$: $\Phi_{*}^{\operatorname{\mathcal{R}}}(\delta;x)=\operatorname{\mathcal{R}}\circ M(\delta;\cdot)\circ[\operatorname{\mathcal{R}}|_{\cal C}]^{-1}(x)$.

The coordinates for the flow on the slow manifold ${\cal C}$ are somewhat arbitrary as the difference between the expressions used by Vandekerckhove et al. [46] and the implicit expression Eq. 4 for $\Phi_{t_{\mathrm{skip}}}$ shows. For the coordinates in the space $\operatorname{dom}\operatorname{\mathcal{L}}$ the diffeomorphism between $\operatorname{dom}\operatorname{\mathcal{L}}$ and ${\cal C}$ is $g\circ\operatorname{\mathcal{L}}$, where $g$ is the stable fiber projection, as implied by Eq. 2. The diffeomorphism can be approximately computed by solving $\operatorname{\mathcal{R}}(M(2{t_{\mathrm{skip}}};\operatorname{\mathcal{L}}(x_{g})))=\operatorname{\mathcal{R}}(M({t_{\mathrm{skip}}};\operatorname{\mathcal{L}}(x)))$ for $x_{g}$ and then using $M({t_{\mathrm{skip}}};\operatorname{\mathcal{L}}(x_{g}))$ as the approximation for $[g\circ\operatorname{\mathcal{L}}](x)$. The approximate diffeomorphism for the expression of Vandekerckhove et al. [46] is $\left[\operatorname{\mathcal{R}}|_{\cal C}\right]^{-1}:x\mapsto M({t_{\mathrm{skip}}};\operatorname{\mathcal{L}}(x_{b}))$.

Marschler et al. [31] proved that the approximation $y_{t_{\mathrm{skip}}}$ is exponentially accurate if $d_{\mathrm{tan}}/d_{\mathrm{tr}}\to 0$: $\|y_{t_{\mathrm{skip}}}-y_{*}\|\sim\exp(-Kd_{\mathrm{tr}}/d_{\mathrm{tan}})$ (for some constant $K$ depending on ${t_{\mathrm{skip}}}$). The error estimates in [31] require that ${t_{\mathrm{skip}}}d_{\mathrm{tan}}/d_{\mathrm{tr}}$ and $({t_{\mathrm{skip}}}+\delta)d_{\mathrm{tan}}/d_{\mathrm{tr}}$ stay bounded from above such that the convergence result is valid in the limit of infinite time scale separation $d_{\mathrm{tan}}/d_{\mathrm{tr}}\to 0$. This means that the assumptions of [31] are similar to those required by schemes involving constrained runs [16, 48, 49]. The analysis left open if the error goes to zero for ${t_{\mathrm{skip}}}\to\infty$ but the time scale separation stays finite: $d_{\mathrm{tan}}/d_{\mathrm{tr}}\in(0,1)$.

Our paper will prove the general a-priori error estimate that $\|y_{t_{\mathrm{skip}}}-y_{*}\|\sim\exp((d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ for ${t_{\mathrm{skip}}}\to\infty$ and fixed $d_{\mathrm{tan}}<d_{\mathrm{tr}}$ under some genericity conditions on $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$. It will also give a convergence result for the derivatives of $y_{t_{\mathrm{skip}}}$ with respect to its argument $x$: $\|\partial^{j}y_{t_{\mathrm{skip}}}-\partial^{j}y_{*}\|\sim\exp(((2j+1)d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ if $(2j+1)d_{\mathrm{tan}}<d_{\mathrm{tr}}$.

Analysis beyond attracting manifolds in slow-fast systems

As mentioned above, equation-free analysis based on lift-evolve-restrict maps is more commonly applied to problems that are assumed to have a fast subsystem, where the fast time scale converges only in a statistical sense to a stationary measure conditioned on the slow variables. In these cases the microscopic time stepper $M(\delta;\cdot)$ operates on measures (or densities). It may be approximated by Monte Carlo simulations on ensembles of initial conditions. Barkley et al. [5] investigated the behaviour of the lift-evolve-restrict map $P(\delta;\cdot)=\operatorname{\mathcal{R}}\circ M(\delta;\cdot)\circ\operatorname{\mathcal{L}}$ where the slow variables were leading moments (thus, $P$ was called moment map in [5]) on prototype examples from the class of stochastic problems. The simplest example from [5] is a scalar stochastic differential equation (SDE), for which the evolution of the probability distribution is governed by a (linear) Fokker-Planck equation (FPE). Hence, the measure of time-scale separation is the size of the spectral gap in the right-hand side of the FPE. The analysis in [5] found that the dynamics of the map $P$ was qualitatively different from the dynamics of the underlying linear FPE. For example, $P$ was nonlinear and had several coexisting fixed points for certain choices of time $\delta$.

Our paper will demonstrate for two different lifting operators $\operatorname{\mathcal{L}}$ that the approximation $y_{t_{\mathrm{skip}}}$, defined by Eq. 4, behaves exactly as predicted by our convergence theorem. In particular, it preserves the metastability features and the linearity of the flow generated by the FPE, thus, addressing the problems highlighted in [5].

2.1 Outline of results

Section 3 states the precise assumptions (time scale separation for decay rates tangential and transversal to the invariant manifold ${\cal C}$ ($d_{\mathrm{tan}}<d_{\mathrm{tr}}$) and transversality of $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$) for exponential convergence:

[TABLE]

(using the convention that $\partial^{0}y=y$ and assuming that the derivatives up to order $j+1$ exist). Estimate Eq. 5 predicts that convergence in ${t_{\mathrm{skip}}}$ is slower for derivatives of higher order. Section 4 demonstrates the convergence rates in ${t_{\mathrm{skip}}}$ for $y_{t_{\mathrm{skip}}}$ and its first two derivatives with respect to $x$ for a singularly perturbed ODE modelling the Michaelis-Menten kinetics (which was also used by [16, 48, 49] for illustration). Section 5 studies the evolution of densities under a scalar SDE with a double-well potential drift term also considered by Barkley et al. [5]. We demonstrate global convergence of implicit equation-free methods for a linear lifting $\operatorname{\mathcal{L}}_{\mathrm{lin}}$. We also demonstrate local convergence for the nonlinear lifting $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$ used in [5].

Section 6 discusses differences between observations of the behaviour in the SDE and the predictions from the theoretical result. These are caused by the numerical errors in the evaluations of lifting, evolution and restriction and their growth along trajectories.

We conclude with an outlook on possible consequences of the results on application of equation-free methods to Monte-Carlo simulations of multi-particle or agent-based systems. One important observation is that in some cases increasing the number of agents or particles does not increase the spectral gap (and, thus, the time scale separation). Didactic examples where the finiteness of the spectral gap is apparent are the dynamic networks as considered by Gross and Kevrekidis [18]. The slow system is an ODE derived from the pair-wise interaction approximation, cutting off an infinite series of ODEs of higher-order interaction terms. The spectral gap between pair-wise interaction terms and triplet interaction terms is finite even in the limit of infinitely large networks.

Thus, the results from Section 3 are potentially applicable to equation-free analysis of stochastic multi-particle systems, where distributions of microscopic initializations are studied. This is in contrast to previous convergence results on constrained runs [16, 48, 49] and implicit lifting [31], which only apply in the limit of infinite time scale separation.

3 Convergence in the case of finite time-scale

separation

We consider a smooth dynamical system

[TABLE]

where $D$ is large. We assume that the flow $M$ generated by Eq. 6,

[TABLE]

has a $d$-dimensional compact relatively invariant manifold ${\cal C}$ (possibly with boundary). That is, trajectories $M(t;u)$ starting in $u\in{\cal C}$ either stay in ${\cal C}$ for all times $t\in\mathbb{R}$, or they stay in ${\cal C}$ until they cross the boundary $\partial{\cal C}$ of ${\cal C}$. We assume that ${\cal C}$ is at least $k_{\max}$ times differentiable. For a point $u\in{\cal C}$, let us denote by ${\cal N}(u)$ the $d$-dimensional tangent space to ${\cal C}$. The following assumption states that attraction transversal to the manifold ${\cal C}$ is faster than attraction or expansion tangential to ${\cal C}$. {assumption}[Hyperbolicity — Separation of time scales and transversal stability] There exists an open neighborhood ${\cal U}$ of the manifold ${\cal C}$, a (possibly nonlinear) projection $g:{\cal U}\mapsto{\cal C}$ (the so-called stable fiber projection), a pair of constants (decay rates) $0<d_{\mathrm{tan}}<d_{\mathrm{tr}}$, and a bound $C$ such that the following two conditions hold.

(tangential expansion/attraction rate) For all points $u\in{\cal C}$ on the manifold, all tangent vectors $v_{1},\ldots,v_{k_{\max}}\in{\cal N}(u)$ and all $t\in\mathbb{R}$ with $M([0,t];u)\subset{\cal C}$

[TABLE]

for all $j\in\{1,\ldots,k_{\max}\}$. 2. 2.

(Stability along transversal fiber projections) For all $u\in{\cal U}$ and all $t>0$ with $M([0,t];g(u))\in{\cal C}$

[TABLE]

for all $j\in\{0,\ldots,k_{\max}\}$.

In Eq. 7 and Eq. 8 we use the convention that $\partial^{j}_{k}M$ is the $j$th-order partial derivative of $M$ with respect to its $k$th argument, and that $\partial_{2}^{0}M$ ($j=0$) is the flow $M$ itself. The norm on the left side of Eq. 8 is the usual operator norm for the multi-linear operators $\partial_{2}^{j}M(t,\cdot)$. The constants $C$, $d_{\mathrm{tr}}$ and $d_{\mathrm{tan}}$ are assumed to be independent of the point $u$ and the time $t$. Assumption Eq. 7 is also made for negative times $t$ (using the convention that $M([0,t];u)$ means $M([t,0];u)$ for $t<0$) such that it is also an assumption about the inverse of $M$, when restricted to ${\cal C}$: $M(-t,\cdot)=M^{-1}(t,\cdot)$. The constant $d_{\mathrm{tr}}$ is the decay rate toward the manifold ${\cal C}$, the constant $d_{\mathrm{tan}}$ is the rate of attraction and expansion along the flow restricted to ${\cal C}$. The main requirement of Section 3 is that $d_{\mathrm{tr}}>d_{\mathrm{tan}}$.

Transversality of restriction and lifting

Second, we assume basic compatibility between

[TABLE]

and the invariant manifold ${\cal C}$: the lifting $\operatorname{\mathcal{L}}$ should map into the neighborhood ${\cal U}$ of ${\cal C}$ in which the stable fiber projection $g$ is defined, and the restriction $\operatorname{\mathcal{R}}$ should be defined on the projection $g$ of the image of $\operatorname{\mathcal{L}}$ along the stable fibers:

[TABLE]

In addition to these compatibility conditions, we impose the following two transversality conditions on lifting $\operatorname{\mathcal{L}}$ and restriction $\operatorname{\mathcal{R}}$. {assumption}[Transversality of $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$]

the projection $g$ is a diffeomorphism between $\operatorname{rg}\operatorname{\mathcal{L}}=\operatorname{\mathcal{L}}(\operatorname{dom}\operatorname{\mathcal{L}})$ and ${\cal C}$. In particular, for all $x\in\operatorname{dom}\operatorname{\mathcal{L}}\subset\mathbb{R}^{d}$

[TABLE] 2. 2.

The map $\operatorname{\mathcal{R}}$, restricted to ${\cal C}$, is a diffeomorphism between ${\cal C}$ and $\mathbb{R}^{d}$. In particular, for all $u\in{\cal C}$ (${\cal N}(u)$ is the tangent space to ${\cal C}$ in $u$)

[TABLE]

Coordinates on the slow manifold ${\cal C}$

The maps $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$ create two natural ways to define local coordinate representations on the invariant manifold ${\cal C}$, one by a map from $\operatorname{dom}\operatorname{\mathcal{L}}$ to ${\cal C}$, one by a map from ${\cal C}$ to $\operatorname{rg}\operatorname{\mathcal{R}}$. For our presentation we choose the representation in $\operatorname{dom}\operatorname{\mathcal{L}}$ coordinates:

[TABLE]

The inverse of $g\circ\operatorname{\mathcal{L}}$ is defined implicitly. Assume that $u_{0}=g(\operatorname{\mathcal{L}}(x_{0}))$ for some $x_{0}\in\operatorname{dom}\operatorname{\mathcal{L}}$. Then for $u\in{\cal C}$ near $u_{0}$ the pre-image $x=(g\circ\operatorname{\mathcal{L}})^{-1}(u)$ is found by solving $\operatorname{\mathcal{R}}(u)=\operatorname{\mathcal{R}}(g(\operatorname{\mathcal{L}}(x)))$ for $x\approx x_{0}$, which has a locally unique solution by Section 3.

We can represent the flow $M$ on ${\cal C}$ as a flow in $\operatorname{dom}\operatorname{\mathcal{L}}$, denoting it by $\Phi_{*}$:

[TABLE]

(for $\delta\in\mathbb{R}$ and $x\in\operatorname{dom}\operatorname{\mathcal{L}}$), where $y$ is given implicitly as solution of a $d$-dimensional system of nonlinear equations

[TABLE]

Section 3 on transversality implies that $\Phi_{*}$ is well defined for small $\delta$ (since $y=x$ is a regular solution of Eq. 11 at $\delta=0$). For larger $\delta$, one can break down the solution into smaller steps by increasing $\delta$ gradually from [math] and tracking the curve $y(\delta)$ of solutions of Eq. 11, which is well parametrized by $\delta$ in every point by Section 3. If $\operatorname{dom}\operatorname{\mathcal{L}}$ is simply connected then this continuation approach makes the implicit solution $y$ used in the definition of $\Phi_{*}$ unique. Let us define the map

[TABLE]

This map $P_{*}$ is well defined and invertible for all $t\in\mathbb{R}$ and $x\in\operatorname{dom}\operatorname{\mathcal{L}}$ for which the trajectory $s\mapsto M(s;g(\operatorname{\mathcal{L}}(x)))$ stays in ${\cal C}$ for all $s$ between [math] and $t$. The implicit definition Eq. 10 of the flow $\Phi_{*}$ on ${\cal C}$ has the following form when expressed with the help of this map $P_{*}$ on $\operatorname{dom}\operatorname{\mathcal{L}}$:

[TABLE]

Since the flow $M(\delta;\cdot)$ is a diffeomorphism on ${\cal C}$, we can replace the times [math] and $\delta$ in the above implicit definition with ${t_{\mathrm{skip}}}$ and ${t_{\mathrm{skip}}}+\delta$ for an arbitrary so-called healing time ${t_{\mathrm{skip}}}\in\mathbb{R}$ (as long as $M([0,{t_{\mathrm{skip}}}];g(\operatorname{\mathcal{L}}(x)))\subset{\cal C}$). So, equivalent to Eq. 13, we have for ${t_{\mathrm{skip}}}>0$ with $M([0,{t_{\mathrm{skip}}}];g(\operatorname{\mathcal{L}}(x)))\subset{\cal C}$

[TABLE]

Convergence Theorem for implicit equation-free

computations with finite time-scale separation

The stable fiber projection $g$ (which is part of the definition of $P_{*}$) is not known in most practical applications. Thus, implicit equation-free computations use the explicit macroscopic time-$t$ map $P$ instead of $P_{*}$ in the equation defining $y$ in Eq. 14:

[TABLE]

such that we may define the approximate flow map

[TABLE]

implicitly in a similar way to Eq. 14. Our general convergence theorem, the following Theorem 3.1, states that $\Phi_{t_{\mathrm{skip}}}$ is well defined for large ${t_{\mathrm{skip}}}$ (that is, the equation in Eq. 16, defining $\Phi_{t_{\mathrm{skip}}}$ implicitly, has a locally unique solution), and that $\partial^{j}\Phi_{t_{\mathrm{skip}}}$ is an approximation of $\partial^{j}\Phi_{*}$ of order $\exp(((2j+1)d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ (including $j=0$ for the map $\Phi_{t_{\mathrm{skip}}}$).

Theorem 3.1** (Convergence of approximate flow map at finite

time-scale separation).**

Let us assume that the microscopic flow $M$ satisfies Section 3 on time-scale separation, and that the maps $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$ satisfy Section 3 on transversality.

Let $\delta_{\max}>0$ and $x\in\operatorname{dom}\operatorname{\mathcal{L}}$ be arbitrary. Let us also assume that $x\in\operatorname{dom}\operatorname{\mathcal{L}}$ maps to a point under $g\circ\operatorname{\mathcal{L}}$ that keeps a positive distance from the boundary $\partial{\cal C}$ of ${\cal C}$ for all times $t\geq-\delta_{\max}$ under $M$. That is,

[TABLE]

Then there exists a $t_{0}\geq\delta_{\max}$ such that $y=\Phi_{t_{\mathrm{skip}}}(\delta;x)$ is well defined by Eq. 16 for all $\delta\in[-\delta_{\max},\delta_{\max}]$ and ${t_{\mathrm{skip}}}>t_{0}$. The estimate

[TABLE]

*holds for all orders $j\in\{0,\ldots,k_{\max}-1\}$ satisfying $(2j+1)d_{\mathrm{tan}}<d_{\mathrm{tr}}$. The constant $C$ depends on $\delta_{\max}$ and $x$, but not on ${t_{\mathrm{skip}}}$. *

Assumption Eq. 17 in Theorem 3.1 is made to permit arbitrarily large ${t_{\mathrm{skip}}}$ while still having Section 3 and Section 3 uniformly satisfied. If one considers $x\in\operatorname{dom}\operatorname{\mathcal{L}}$ for which the trajectory $t\mapsto M(t;g(\operatorname{\mathcal{L}}(x)))$ leaves ${\cal C}$ (by crossing the boundary $\partial{\cal C}$) then one has to put restrictions on $\delta$ and ${t_{\mathrm{skip}}}$ to avoid crossing $\partial{\cal C}$. The theorem permits negative integration times $\delta$ shorter than $t_{0}\leq{t_{\mathrm{skip}}}$ and positive integration times larger than ${t_{\mathrm{skip}}}$ as long as the factor $\exp(d_{\mathrm{tan}}|\delta|)$ is of order $1$. Since $1/d_{\mathrm{tan}}$ and $1/d_{\mathrm{tr}}$ are the time scales of the dynamics inside the invariant manifold ${\cal C}$ and transversal to it, the theorem covers time steps of length $\delta$ of order $1$ in the slow ($1/d_{\mathrm{tan}}$) time scale. The statement of Theorem 3.1 does not require that the time-scale separation $d_{\mathrm{tan}}/d_{\mathrm{tr}}$ goes to zero for convergence of the approximate map. It only requires that $d_{\mathrm{tan}}<d_{\mathrm{tr}}$, where $d_{\mathrm{tr}}$ is the attraction rate along fibers (see Eq. 8) and $d_{\mathrm{tan}}$ is the attraction and expansion rate tangential to the invariant manifold ${\cal C}$. Since the constant $C$ in Eq. 18 is independent of ${t_{\mathrm{skip}}}$ it can be chosen uniformly for compact domains $\operatorname{dom}\operatorname{\mathcal{L}}$.

Outline of proof of Theorem 3.1

Existence and error of $\Phi_{t_{\mathrm{skip}}}$: For the proof of Theorem 3.1 we have to analyze the difference between the two defining equations for the approximate solution $y$ and the true solution $y_{*}=\Phi_{*}(\delta;x)$, both depending on $x$ as a parameter:

[TABLE]

Rearranging the difference between Eq. 19 and Eq. 20, we obtain an implicit fixed-point problem for $y$ (recall that $P_{*}(t;\cdot)=\operatorname{\mathcal{R}}\circ M(t;\cdot)\circ g\circ\operatorname{\mathcal{L}}$ and $P(t;\cdot)=\operatorname{\mathcal{R}}\circ M(t;\cdot)\circ\operatorname{\mathcal{L}}$):

[TABLE]

The norms of both terms in square brackets, $P_{*}({t_{\mathrm{skip}}};y)-P({t_{\mathrm{skip}}};y)$ and $P({t_{\mathrm{skip}}}+\delta;x)-P_{*}({t_{\mathrm{skip}}}+\delta;x)$, are of order $\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$ by Section 3, equation Eq. 8 (transversal stability with rate $d_{\mathrm{tr}}$ of ${\cal C}$). For the same reason, the Lipschitz constant of $P_{*}({t_{\mathrm{skip}}};y)-P({t_{\mathrm{skip}}};y)$ is of order $\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$, too. By Section 3, equation Eq. 7 (tangential decay rate inside the manifold is less than $d_{\mathrm{tan}}$), and Section 3 on transversality of $\operatorname{\mathcal{L}}$ and $\operatorname{\mathcal{R}}$, the inverse of $P_{*}({t_{\mathrm{skip}}};\cdot)$ has a local Lipschitz constant of order $\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$ near $y_{*}$. These two facts enable us to apply the Banach Contraction Mapping Principle to Eq. 21 to obtain a unique solution $y\approx y_{*}$ for large ${t_{\mathrm{skip}}}$. More precisely, $y-y_{*}$ is of order $\exp((d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$.

Inductive proof of error estimate for derivatives: We differentiate Eq. 21 with respect to $x$ in its fixed point $y({t_{\mathrm{skip}}};x)$ up to $j$ times and then re-arrange the resulting equation for the $j$th-order derivatives of $y$ and $y_{*}$ into the form

[TABLE]

In Eq. 22 we abbreviated $\partial_{2}P_{*}({t_{\mathrm{skip}}};\cdot)=\partial P_{*}(\cdot)$, $\partial^{j}y=\partial^{j}_{2}y({t_{\mathrm{skip}}};x)$ and dropped the argument $x$ from $y_{*}$ and the arguments ${t_{\mathrm{skip}}}$ and $x$ from $y$. The remainder $r$ is less than $C\exp(((2j-1)d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ for some constant $C$ by induction hypothesis. The implicit expression Eq. 22 for $\partial^{j}y-\partial^{j}y_{*}$ shows why errors in derivatives of the solution can grow for increasing ${t_{\mathrm{skip}}}$ and insufficient time scale separation: the norms of $\partial P_{*}(y)-\partial P_{*}(y_{*})$ and of $\left[\partial P_{*}(y)\right]^{-1}$ are of order $\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$ due to Eq. 7. The details of the proof are given in Appendix A.

4 Example: Michaelis-Menten kinetics

To illustrate the consequences of error estimate Eq. 18, we look at a model for Michaelis-Menten kinetics with explicit time scale separation as studied in [34, 16, 48, 49]. The system is given in $\mathbb{R}^{D}$ with $D=2$ as

[TABLE]

where $x\in\mathbb{R}$ is the slow variable, $y\in\mathbb{R}$ is the fast variable, and $\varepsilon$ measures the time scale separation. The parameters $\kappa=1,\lambda=0.5$ and $\varepsilon=0.01$ are kept fixed throughout this section.

In the singular case $\varepsilon=0$, system Eq. 23 has a critical manifold $\mathcal{C}_{0}$ of equilibria, given by the graph ${\cal C}_{0}=\{(x,y):y=x/(x+\kappa)\}$. For positive $\varepsilon$, the system has a transversally stable invariant manifold, which can be represented as a graph ${\cal C}_{\varepsilon}=\{(x,y):y=h_{\varepsilon}(x)\}$, such that $d=1$. In this section we put a subscript $\varepsilon$ on quantities to indicate their dependence on the parameter $\varepsilon$ (so, writing, for example, ${\cal C}_{\varepsilon}$ instead of ${\cal C}$). The graph $h_{\varepsilon}$ can be expanded in $\varepsilon$ for small $\varepsilon>0$:

[TABLE]

We plan to compare an equation-free approximate flow $\Phi_{t_{\mathrm{skip}}}(\delta;\cdot)$, which is constructed below, to the true flow $\Phi_{*}(\delta;\cdot)$. For this simple example we may approximate the true flow $\Phi_{*}(\delta;\cdot)$ by obtaining an approximation of the stable fiber projection $g_{\varepsilon}$. For $\varepsilon\to 0$ in Eq. 23, $g_{\varepsilon}$ has the limit $g_{\varepsilon}(x,y)\to(x,h_{0}(x))$. Thus, every point in phase space is approximately projected along vertical lines, as shown in Fig. 4(a). For positive $\varepsilon$ this stable fiber projection persists and is perturbed by terms of order $\varepsilon$. A general approximation algorithm for stable fibers in slow-fast systems was provided by Kristiansen et al. [26]. However, we need the stable fibers only to a degree of accuracy that permits us to compare $\Phi_{*}(\delta;\cdot)$ to $\Phi_{t_{\mathrm{skip}}}(\delta;\cdot)$ for demonstration purposes. Thus, we expand $g_{\varepsilon}=(g_{x,\varepsilon}(x,y),g_{y,\varepsilon}(x,y))$ in $\varepsilon$. Since $g_{\varepsilon}$ projects onto the manifold ${\cal C}_{\varepsilon}$, we know that $g_{y,\varepsilon}(x,y)=h_{\varepsilon}(g_{x,\varepsilon}(x,y))$. The first-order expansion of $g_{x,\varepsilon}$ is

[TABLE]

Figure 4 was produced using a third-order expansion of $h_{\varepsilon}$ and $g_{x,\varepsilon}$. The supplementary material provides Matlab code which reproduces the graphs in Fig. 4 and computes the expansion coefficients for $h_{\varepsilon}$ and $g_{x,\varepsilon}$ to third order (see also Appendix B). Figures 4(a) and 4(c) show the phase space geometry. The slow manifold ${\cal C}_{\varepsilon}=\{(x,y):x=h_{\varepsilon}(x)\}$ is shown in green, the stable fibers (pre-images of selected points $(x_{j},h_{\varepsilon}(x_{j}))$ on the slow manifold under the projection $g_{\varepsilon}$) are the almost straight grey lines, and some sample trajectories of Eq. 23 are shown in purple. After a rapid transient all trajectories approach the slow manifold ${\cal C}_{\varepsilon}$, given approximately by Eq. 24. Furthermore, Fig. 4(c) shows in more detail how initial conditions on the same stable fiber, defined as the pre-image of $x_{0}$ under $g_{x,\varepsilon}$, $G_{\mathrm{pre}}=\{(x,y):g_{x,\varepsilon}(x,y)=x_{0}\}$, collapse onto the same slow limiting trajectory (trajectories shown in red in Fig. 4(c)). In contrast, initial conditions with the same $y$-component do so only up to an error of order $\varepsilon$ (trajectories shown in purple in Fig. 4(c)).

To define the approximate flow $\Phi_{t_{\mathrm{skip}}}$, we specify the restriction and lifting operators, $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$, for the Michaelis-Menten system as

[TABLE]

The approximate time-$\delta$ map $\Phi_{t_{\mathrm{skip}}}(\delta;\cdot)$ on the slow manifold is determined by the root $z_{t_{\mathrm{skip}}}$ of

[TABLE]

and setting $\Phi_{t_{\mathrm{skip}}}(\delta;x):=z_{t_{\mathrm{skip}}}$ (cf. Eq. 16). We compare this to the true solution (or, rather, the alternative approximation by expansion) $\Phi_{*}(\delta;x)$, determined by the root $z_{*}$ of

[TABLE]

setting $\Phi_{*}(\delta;x):=z_{*}$. Note, that $F$ depends on ${t_{\mathrm{skip}}}$ and $\delta$, and $F_{*}$ depends on $\delta$, which are not included in the list of arguments to simplify notation. We solve Eq. 26 and Eq. 27 using a Newton iteration with tolerance $10^{-12}$, where we approximate $M$ with the fifth-order component of the DOPRI45 Runge-Kutta scheme with fixed step size $0.1$ for the ODE Eq. 23. We approximate the first two derivatives of $\Phi_{*}$ and $\Phi_{t_{\mathrm{skip}}}$ (and the Jacobians needed inside the Newton iteration) by central finite differences with step size $\Delta_{z}=10^{-4}$. The supplementary material contains didactic implementations of $\Phi_{*}$ and $\Phi_{t_{\mathrm{skip}}}$ for this example in the form of matlab code and its published output. The error

[TABLE]

for $x_{0}=-0.1$ and $\delta=25$ is shown in Fig. 4(e) for a range of healing times ${t_{\mathrm{skip}}}\in[0;30]$. Since $\varepsilon=10^{-2}$ and $d_{\mathrm{tan}}\sim\varepsilon$, the quantity $\exp(d_{\mathrm{tan}}\delta)$ is of order $1$, as required by Theorem 3.1.

The error plot Fig. 4(e) shows that $\Phi_{t_{\mathrm{skip}}}$, $\partial_{2}\Phi_{t_{\mathrm{skip}}}$ and $\partial_{2}^{2}\Phi_{t_{\mathrm{skip}}}$ approach a limit at an exponential rate in ${t_{\mathrm{skip}}}$ up to an accuracy determined by the accuracy of the asymptotic expansion of $\Phi_{*}$ ($\sim\varepsilon^{4}$), round-off errors in the finite difference approximations of the derivatives ($\sim 10^{-4}$ for $\partial_{2}^{2}\Phi_{t_{\mathrm{skip}}}$), and the tolerance of the Newton iteration. Furthermore, the convergence rate is indeed lower for higher orders of the derivative of $\Phi_{t_{\mathrm{skip}}}$ as the estimate Eq. 18 in Theorem 3.1 suggests. The slope of $\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$ is included as a lower bound for the error for comparison.

We also observe that the error of the flow and its derivatives is acceptably small ($\approx 10^{-3}$) even for the minimal value ${t_{\mathrm{skip}}}=0$. Since $\operatorname{\mathcal{R}}\circ\operatorname{\mathcal{L}}$ equals the identity (see Eq. 25), the implicit equation-free method turns into an explicit formulation if ${t_{\mathrm{skip}}}=0$. However, the geometry of system Eq. 23 is not generic. The system Eq. 23 is given in an explicit slow-fast form with one fast and one slow variable. This leads with our choice of lifting and restriction to the degenerate situation that the stable fiber projection $g_{\varepsilon}$ is aligned with lifting and restriction to first order: $\operatorname{\mathcal{R}}\circ g_{0}\circ\operatorname{\mathcal{L}}=I$ such that $\operatorname{\mathcal{R}}\circ g_{\varepsilon}\circ\operatorname{\mathcal{L}}$ is a small (order $\varepsilon$) perturbation of the identity. In this case the explicit equation-free method without healing time ($y=\operatorname{\mathcal{R}}(M(0;\operatorname{\mathcal{L}}(y)))=\operatorname{\mathcal{R}}(M(\delta;\operatorname{\mathcal{L}}(x)))$, such that ${t_{\mathrm{skip}}}=0$) is accurate up to order $\varepsilon$. To create a situation with a generic arrangement of the stable fiber projection $g_{\varepsilon}$, we study a rotated system of the Michaelis-Menten dynamics (which was also used by [16]).

We apply the rotation matrix $R$ to the system in order to obtain the dynamics in the new coordinates $(v,w)^{T}\in\mathbb{R}^{2}$ by

[TABLE]

is the microscopic simulator in the new $(v,w)^{T}$ coordinates. In this rotated system the time scale separation is no longer visible between $v$ and $w$, since the slow and fast variables are mixed. Figure 4(b) shows the phase space geometry with slow manifold (green), stable fibers (grey) and sample trajectories. The initial transients are no longer following a straight line parallel to a coordinate axis such that both, $v$ and $w$, change rapidly during transients. This situation is expected in a generic situation when one applies equation-free methods without precise knowledge about the slow and fast variables. We use the same restriction and lifting operators as defined in Eq. 25 (but in the new coordinates $(v,w)$: $\operatorname{\mathcal{L}}(x)=(x,0.5)^{T}$, $\operatorname{\mathcal{R}}(v,w)=v$). All parameter values are as in the unrotated system Eq. 23, otherwise. The error $E^{j}$, defined in Eq. 28, is now much larger: it is of order $1$ for ${t_{\mathrm{skip}}}=0$; see Fig. 4(f). In the implicit framework the error decreases with increasing healing time ${t_{\mathrm{skip}}}$ down to $10^{-8}$ for ${t_{\mathrm{skip}}}=20$. Note again that the slope of the curves is smaller for higher-order derivatives as predicted by the estimate Eq. 18 in Theorem 3.1. The error for the flow is bounded from below again by the accuracy of the asymptotic expansion for $\Phi_{*}$, the accuracy of the Newton iteration and round-off error caused by the finite difference approximations of $\partial_{2}^{j}\Phi_{*}$ and $\partial_{2}^{j}\Phi_{t_{\mathrm{skip}}}$.

5 Application: stochastic dynamics

A common area where equation-free methods are applied are multi-particle systems where slow dynamics emerges for macroscopic (typically averaged) quantities, e. g. [29, 5]. More precisely, the macroscopic quantities are assumed to satisfy a low-dimensional stochastic differential equation (SDE). For example, the SDE could be assumed to be of the form $\mathop{}\!\mathrm{d}x=f(x)\mathop{}\!\mathrm{d}t+\sigma\mathop{}\!\mathrm{d}W_{t}$, where the noise term $\sigma\mathop{}\!\mathrm{d}W_{t}$ approximates the microscopic fluctuation as white noise and the deterministic part $f(x)$ is the systematic average drift of the macroscopic quantities. Givon et al. [17] review rigorous results concerning dimension reduction of SDEs.

Typically, a stochastic simulation is performed not just once, but for an ensemble of initial conditions and realizations (as part of a Monte Carlo simulation). At the level of an SDE, an ensemble of initial conditions corresponds to (a sampling of) an initial distribution density $\rho(x)$. In this section we restrict ourselves to the study of a scalar SDE of the form

[TABLE]

where $W_{t}$ is a Wiener process, an example for which explicit equation-free methods have been thoroughly analyzed by Barkley et al. [5]. As in [5], we set the noise strength $\sigma$ equal to $1$ in Eq. 30 without loss of generality. The potential

[TABLE]

forms for $\mu>0$ a double well with two local minima $Q_{\pm}$ and a local maximum $Q_{s}$ (see lower panel of Fig. 5 for a graph of $V$). The parameters $\mu$ and $\nu$ determine the depth and the asymmetry of the double-well potential, respectively. We use $\mu=6$, $\nu=0.3$ such that $Q_{-}<Q_{s}<Q_{+}$ and the well around $Q_{-}$ is deeper than the well around $Q_{+}$. The microscopic simulation is a Monte Carlo simulation of Eq. 30 starting from initial (ensemble) density $\rho_{0}(Q)$ of initial conditions. Thus, the phase space is the space of possible initial distributions in $Q$, which has dimension $D$ equal to infinity. Strictly, the infinite-dimensional case is outside of the scope of Theorem 3.1. However, the observations to follow agree with the convergence predicted by the theorem for reasons that will be discussed after defining lifting, evolution and restriction. We will make the connection to multi-particle systems or high-dimensional SDEs in Section 6.

5.1 Lifting, evolution and restriction for distributions

The evolution of the probability density function (pdf) $\rho(Q,t)$ for the realization of Eq. 30 is determined by the Fokker-Planck equation with $\sigma=1$,

[TABLE]

The right-hand side of Eq. 32 is linear, of the form

[TABLE]

where the operator $L:\mathbb{H}^{2}_{1}(\mathbb{R};\mathbb{R})\mapsto\mathbb{L}^{2}_{1}(\mathbb{R};\mathbb{R})$ is self-adjoint with respect to the scalar product

[TABLE]

The space $\mathbb{L}^{2}_{1}(\mathbb{R};\mathbb{R})$ is in our case the space of all measurable functions $u:\mathbb{R}\to\mathbb{R}$ with $\int_{\mathbb{R}}u^{2}(x)/\varphi_{1}(x)\mathop{}\!\mathrm{d}x<\infty$ (a subset of $\mathbb{L}^{2}(\mathbb{R};\mathbb{R})$, which has the scalar product $\langle\rho_{1},\rho_{2}\rangle=\int_{\mathbb{R}}\rho_{1}(Q)\rho_{2}(Q)\mathop{}\!\mathrm{d}Q$). The space $\mathbb{H}^{\ell}_{1}(\mathbb{R};\mathbb{R})$ is the space of all $u\in\mathbb{L}^{2}_{1}(\mathbb{R};\mathbb{R})$ with $u^{(j)}\in\mathbb{L}^{2}_{1}(\mathbb{R};\mathbb{R})$ for all $j\leq\ell$. The spectrum of $L$ is real and consists of point spectrum only. It has the form $0=\lambda_{1}>\lambda_{2}>\ldots$ with eigenvectors $\varphi_{j}(Q)\in\mathbb{H}^{2}_{1}(\mathbb{R};\mathbb{R})$ that can be orthonormalized with respect to $\langle\cdot,\cdot\rangle_{1}$. The function $\varphi_{1}$ is the eigenvector for the trivial eigenvalue $\lambda_{1}=0$ (which is present due to the preservation of total probability $\int_{\mathbb{R}}\rho(Q,t)\mathop{}\!\mathrm{d}Q$ along trajectories). The spectrum and the corresponding eigenfunctions $\varphi_{j}$ are shown in Fig. 6 (left panel), together with the $\mathbb{L}^{2}$-adjoint eigenfunctions $\varphi_{j}/\varphi_{1}$ (right panel). A solution of the Fokker-Planck equation Eq. 32 can be expanded in the eigenfunctions of $L$ with time-dependent coefficients $a_{j}(t)$:

[TABLE]

The coefficients satisfy $\dot{a}_{j}(t)=\lambda_{j}a_{j}(t)$ for all $j$, and the series $\sum_{j=1}^{\infty}a_{j}^{2}$ converges for all $t>0$. The orthonormality of the basis $\{\varphi_{j}:j\geq 1\}$ with respect to $\langle\cdot,\cdot\rangle_{1}$, defined in Eq. 34, implies that

[TABLE]

Since $\lambda_{1}=0$, $a_{1}(t)$ equals $a_{1}(0)$ for all times $t\geq 0$. One usually chooses $a_{1}(0)=1$ such that $\rho(Q,t)$ converges to the stationary density $\varphi_{1}(Q)$ for $t\to\infty$. While Theorem 3.1 was only formulated for flows in $\mathbb{R}^{D}$, the linearity of $L$ implies that statements identical to Theorem 3.1 can be made for the PDE Eq. 32. Instead of Fenichel’s Theorem on invariant manifolds in ODEs [15] (persistence and regularity of invariant manifolds and fiber projections) we rely on the spectral mapping properties for the linear operator $L$. For any chosen dimension $d$ of the slow variables, the slow manifold ${\cal C}$ is the subspace spanned by $\varphi_{1},\ldots,\varphi_{d}$. Instead of the stable fiber projection in the finite-dimensional case, we have the linear spectral (for $L$) projection $g:\mathbb{L}^{2}_{1}(\mathbb{R};\mathbb{R})\mapsto{\cal C}$ ($M$ is also linear, such that we write $M(t)\rho$ and $g\rho$) which is explicitly known in terms of the eigenvectors of $L$:

[TABLE]

With this definition of ${\cal C}$ and $g$ the decay and growth properties of the evolution map $M$ replacing Section 3 are

[TABLE]

and some constant $C$, such that $d_{\mathrm{tan}}=-\lambda_{d}$, $d_{\mathrm{tr}}=-\lambda_{d+1}$. The approximation statement of Theorem 3.1 then follows immediately from error estimates for finite-dimensional matrices and will be derived after the definition of the restriction and lifting operators. The lifting and restriction operators are chosen to map from a macroscopic description of $\rho$, for example, by moments, to the full density $\rho$ and vice versa. In particular, we will investigate the behaviour of implicit equation-free methods for $d=3$ using the following restriction and two different lifting operators:

[TABLE]

Thus, $\operatorname{\mathcal{R}}$ projects a density onto its first $d$ moments (counting from the zeroth moment, which is preserved by $M$ since $\lambda_{1}=0$). In a Monte-Carlo simulation the zeroth moment would correspond to the (possibly scaled) number of realizations. The functions $\rho_{j}$ in the definition Eq. 41 of $\operatorname{\mathcal{L}}_{\mathrm{lin}}$ are arbitrary in $\mathbb{L}^{2}_{1}(\mathbb{R};\mathbb{R})$ with $\int_{\mathbb{R}}\rho_{j}(Q)\mathop{}\!\mathrm{d}Q=1$, which ensures that $\int_{\mathbb{R}}\operatorname{\mathcal{L}}_{\mathrm{lin}}(x)(Q)\mathop{}\!\mathrm{d}Q=\sum_{j=1}^{d}x_{j}$ is conserved under $M(t)$. For $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$, the $x_{1}$ component is preserved under $M(t)$ and becomes the first component of $\operatorname{\mathcal{R}}$ such that always $[\operatorname{\mathcal{R}}M(t)\operatorname{\mathcal{L}}_{\mathrm{Gauss}}x]_{1}=x_{1}$.

For the combination of $\operatorname{\mathcal{L}}_{\mathrm{lin}}$ and $\operatorname{\mathcal{R}}$ all components of the lift-evolve-restrict map $P(t;\cdot)=\operatorname{\mathcal{R}}\circ M(t;\cdot)\circ\operatorname{\mathcal{L}}$ and its exact counterpart $P_{*}(t;\cdot)=\operatorname{\mathcal{R}}\circ M(t;\cdot)\circ g\circ\operatorname{\mathcal{L}}$ from Section 3 are linear such that we can reduce the study of convergence for arbitrary coordinates $x$ to convergence estimates for matrices.

The combination of $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$ and $\operatorname{\mathcal{R}}$ was studied in detail in [5] for explicit equation-free methods, where the authors observed that the nonlinearity of $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$ introduced a nonlinearity in the moment map and that the resulting flow depended qualitatively on the choice of the healing time ${t_{\mathrm{skip}}}$. We will demonstrate that for $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$ the implicitly defined flow $\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}$ converges to a nonlinear transformation of the linear flow $M(t)|_{\cal C}$. Since the $x_{1}$ component does not change under $P(t;\cdot)$ and $P_{*}(t;\cdot)$ for $\operatorname{\mathcal{L}}=\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$, it can be ignored, making the choice of $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$ and $\operatorname{\mathcal{R}}$ identical to the situation studied in [5].

We use the MATLAB [21] package chebfun [10, 14] to numerically compute the spectrum and eigenfunctions of $L$, the flow $M$, the projection $g$, restriction and lifting for the example potential $V$ given in Eq. 31. The package chebfun uses Chebyshev polynomials of adaptive degree to approximate arbitrary functions on finite intervals to optimal precision. For a typical result, as shown in Fig. 6 the degree is larger than $100$ ($394$ for the left panel, $941$ for the right panel). The numerically computed spectrum of $L$ is

[TABLE]

Note that $\lambda_{1}=0$ is the correct value for the first eigenvalue on an infinite domain. In numerical computations we choose a bounded domain $[-10,10]$ with Dirichlet boundary conditions, leading to a small probability of escape from the domain. The spectrum and the corresponding eigenfunctions $\varphi_{j}$ are shown in Fig. 6. The eigenvector $\varphi_{1}$ corresponds to the stationary solution of the Fokker-Planck equation and $\varphi_{2}$ is the mode representing escape from one well to another.

5.2 Convergence for the linear lifting operator $\operatorname{\mathcal{L}}_{\mathrm{lin}}$ with

$d=3$

We express the maps $P_{*}(t;\cdot)$ and $P(t;\cdot)$ in terms of $M$, the eigenvectors $\varphi_{j}$ and the scalar product $\langle\cdot,\cdot\rangle_{1}$, initially for a general dimension $d$. The exact macroscopic flow $\Phi_{*}$ is defined using the map $P_{*}(t;\cdot)$ in Eq. 14, and the approximate macroscopic flow $\Phi_{t_{\mathrm{skip}}}$ is defined using the map $P(t;\cdot)$ in Eq. 16. The definitions Eq. 40 for $\operatorname{\mathcal{R}}$ and Eq. 41 for $\operatorname{\mathcal{L}}_{\mathrm{lin}}$ imply

[TABLE]

where $k=1,\ldots,d$. Using the $d\times d$ matrices ($k,\ell,j=1,\ldots,d$)

[TABLE]

we can express the map $P_{{\mathrm{lin}},*}(t)$ and the exact slow flow $\Phi_{{\mathrm{lin}},*}(\delta)$ in the form

[TABLE]

General Section 3 on transversality of $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}$ for Theorem 3.1, when applied to the SDE Eq. 30 and $\operatorname{\mathcal{R}}$ and $\operatorname{\mathcal{L}}_{\mathrm{lin}}$, demands the regularity of the matrices $R_{d}$ and $T_{\mathrm{lin}}$. If both matrices are indeed regular then $P_{{\mathrm{lin}},*}$ is invertible: $P_{{\mathrm{lin}},*}(t)^{-1}=T_{\mathrm{lin}}^{-1}M_{d}(-t)R_{d}^{-1}$. Thus, the claim of Theorem 3.1 can be simplified to a statement about perturbations of matrices using the quantities

[TABLE]

The estimate for $r(t)$ follows from Eqs. 45 and 46:

[TABLE]

The integrand contains the spectral projection $\rho\mapsto\rho-\sum_{\ell=0}^{d}\langle\varphi_{\ell},\rho_{j}\rangle\varphi_{\ell}$ onto the complement of the space spanned by $\varphi_{1},\ldots,\varphi_{d}$. On the complement of ${\cal L}(\{\varphi_{1},\ldots,\varphi_{d}\})$ the evolution operator $M(t)$ decays exponentially with rate $\lambda_{d+1}$ in time. Together with the boundedness of $\operatorname{\mathcal{R}}$ and the spectral projection, this decay of $M(t)$ implies estimate Eq. 51. Since $P_{\mathrm{lin}}(t;\cdot)$ is linear, the approximate flow map $\Phi_{{\mathrm{lin}},{t_{\mathrm{skip}}}}(\delta)$ is given by

[TABLE]

assuming the inverse of $P_{\mathrm{lin}}({t_{\mathrm{skip}}})$ exists (all involved matrices have dimension $d\times d$). The linear expressions for the exact flow $\Phi_{{\mathrm{lin}},*}$, Eq. 49, and $\Phi_{{\mathrm{lin}},{t_{\mathrm{skip}}}}$, Eq. 52, imply

[TABLE]

The exponential estimates $r(t)\leq C\exp(\lambda_{d+1}t)$ in Eq. 51 and $n(t)\leq\exp(-\lambda_{d}t)$ in Eq. 50 immediately imply the statement of Theorem 3.1 (given that $\Phi_{{\mathrm{lin}},*}$ is globally bounded).

The semilogarithmic plot in Fig. 7(a) shows the difference between $\Phi_{{\mathrm{lin}},{t_{\mathrm{skip}}}}(\delta)$ and $\Phi_{{\mathrm{lin}},*}(\delta)$ (blue line with circles) for $d=3$, for a linear basis of three Gaussians $\rho_{j}$ (with variance $1$ and means $-1.5$, $-0.5$ and $1$), $\delta=0.1$ and the double-well potential well $V$ with parameters $\mu=6$, $\nu=0.3$. The decay rate inside the slow manifold ${\cal C}=\operatorname{span}(\varphi_{1},\varphi_{2},\varphi_{3})$ is $d_{\mathrm{tan}}=-\lambda_{3}\approx 5.71$ and the attraction rate toward ${\cal C}$ is $d_{\mathrm{tr}}=-\lambda_{4}\approx 10.3$. Figure 7(a) also shows the two components of the error $\Phi_{{\mathrm{lin}},{t_{\mathrm{skip}}}}(\delta)-\Phi_{{\mathrm{lin}},*}(\delta)$ and their theoretical estimates:

(In yellow) The difference $r({t_{\mathrm{skip}}})$ between $P_{\mathrm{lin}}({t_{\mathrm{skip}}})=\operatorname{\mathcal{R}}\circ M({t_{\mathrm{skip}}})\circ\operatorname{\mathcal{L}}_{\mathrm{lin}}$ and $P_{{\mathrm{lin}},*}({t_{\mathrm{skip}}})=\operatorname{\mathcal{R}}\circ M({t_{\mathrm{skip}}})\circ g\circ\operatorname{\mathcal{L}}_{\mathrm{lin}}=R_{d}M_{d}({t_{\mathrm{skip}}})T_{\mathrm{lin}}$, which decays according to the attraction toward ${\cal C}$ until it reaches the limits of numerical accuracy of chebfun ($\sim 10^{-8}$): $\|P_{\mathrm{lin}}({t_{\mathrm{skip}}})-P_{{\mathrm{lin}},*}({t_{\mathrm{skip}}})\|\sim\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$. 2. 2.

(In red) The norm $n({t_{\mathrm{skip}}})$ of the inverse of $P_{\mathrm{lin}}({t_{\mathrm{skip}}})$, which grows like $\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$. Figure 7(a) shows the inverse (the minimal singular value).

The overall error Eq. 53 is approximately the product of these two components, which is proportional to $\exp((d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ (shown as a blue dashed line in Fig. 7(a)). In particular, the combination of $\|P_{{\mathrm{lin}},*}({t_{\mathrm{skip}}})^{-1}\|\sim\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$, $\|P_{\mathrm{lin}}({t_{\mathrm{skip}}})-P_{{\mathrm{lin}},*}({t_{\mathrm{skip}}})\|\sim\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$ and $d_{\mathrm{tan}}<d_{\mathrm{tr}}$ implies that $P_{\mathrm{lin}}({t_{\mathrm{skip}}})$ is invertible for sufficiently large ${t_{\mathrm{skip}}}$ and that $\|P_{\mathrm{lin}}({t_{\mathrm{skip}}})^{-1}\|\sim\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$.

Figure 7(b) shows a phase portrait of the exact flow in the coordinates in $\operatorname{dom}\operatorname{\mathcal{L}}_{\mathrm{lin}}$. Since $\Phi_{{\mathrm{lin}},*}$ and $\Phi_{{\mathrm{lin}},{t_{\mathrm{skip}}}}$ both preserve the quantity $\sum_{j=1}^{d}x_{j}$ (which corresponds to $\int_{\mathbb{R}}\operatorname{\mathcal{L}}_{\mathrm{lin}}(x)(Q)\mathop{}\!\mathrm{d}Q$), we set $x_{3}=1-x_{1}-x_{2}$ in the initial values for the sample trajectories, keeping $\sum_{j=1}^{d}x_{j}=1$ along trajectories without loss of generality. This leads to an affine flow in the $(x_{1},x_{2})$-plane with a non-trivial fixed point (shown in black in Fig. 7(b)). The coloring along the sample trajectories illustrates the extreme difference in the time scale along the directions corresponding to $\lambda_{2}$ ($\approx-10^{-7}$; escape between wells), mostly evolving on time scales $\gg 10^{4}$ (dark red in Fig. 7(b)), and $\lambda_{3}$ ($\approx-5.71$; relaxation into the nearest well), mostly decaying on time scale of order $1$ and less (blue and light blue in Fig. 7(b)).

Remark: Densities with sign changes in

Section 5.2

The phase portrait Fig. 7(b) of the exact flow $\Phi_{{\mathrm{lin}},*}$ includes coordinates $x=(x_{1},x_{2},x_{3})$ where the lifted initial density $\operatorname{\mathcal{L}}_{\mathrm{lin}}(x)$ has sign changes. This is not unphysical. If one performs Monte Carlo simulations with ensembles on the example with the lifting operator $\operatorname{\mathcal{L}}_{\mathrm{lin}}(x)$, one would run a Monte Carlo simulation on an ensemble for each of the three initial densities $\rho_{j}$. Then one would sum the densities at the end of the simulation with the weights $x_{j}$ ($j=1,\ldots,3$). These weights can be negative to get a combined density.

5.3 Convergence for the nonlinear lifting operator $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$

The exact and approximate lift-evolve-restrict maps for lifting with a Gaussian distribution of mass $x_{1}$, mean $x_{2}$ and variance $x_{3}$, of the form $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}(x)=q\mapsto x_{1}\exp(-(q-x_{2})^{2}/(2x_{3}))/\sqrt{2\pi x_{3}}$, are given by

[TABLE]

where $k=1,\ldots,d$ ($d=3$). The flow $M(t)$ preserves the integral of the initial distribution such that $P_{\mathrm{Gauss}}(t;x)_{1}=x_{1}$ and $P_{{\mathrm{Gauss}},*}(t;x)_{1}=x_{1}$. Thus, we can fix $x_{1}=1$ without loss of generality and focus on the dynamics in the $(x_{2},x_{3})$-plane in $\operatorname{dom}\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$.

Phase portrait of the exact flow $\Phi_{{\mathrm{Gauss}},*}$

The exact flow $\Phi_{{\mathrm{Gauss}},*}$ on ${\cal C}$ in the coordinates of $\operatorname{dom}\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$ is a nonlinear transformation of the linear map $M_{d}(t)=\operatorname{diag}\left[\exp(\lambda_{\ell}t)_{\ell=1}^{\ell=d}\right]:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$, defined in Eq. 47 (with $d=3$). We call the nonlinear transformation

[TABLE]

In particular $T_{\mathrm{Gauss}}(x)_{1}=x_{1}$ by construction. Using $T_{\mathrm{Gauss}}$, $M_{d}$ and the matrix $R_{d}$ (defined in Eq. 47), the map $P_{{\mathrm{Gauss}},*}(t;x)$, and the exact flow $\Phi_{{\mathrm{Gauss}},*}$ are given by (using the notation $T_{\mathrm{Gauss}}^{-1}$ for the inverse of the nonlinear map $T_{\mathrm{Gauss}}$)

[TABLE]

where all involved quantities are maps from $\mathbb{R}^{3}$ to $\mathbb{R}^{3}$. Since the map $T_{\mathrm{Gauss}}$ is nonlinear, it is not clear if the inverse exists for all $x\in\mathbb{R}^{3}$, or if it is unique where it exists. Figure 8(a) shows the contours of $T_{\mathrm{Gauss}}(x)_{2}$ (in black) and $T_{\mathrm{Gauss}}(x)_{3}$ (in blue; remember that $T_{\mathrm{Gauss}}(x)_{1}=x_{1}$), and the norm of $[\partial T_{\mathrm{Gauss}}(x)]^{-1}$ as color shading (in logarithmic scale). Since the difference in time scale between motion along $\varphi_{2}$ and motion along $\varphi_{3}$ is large ($0>\lambda_{2}\gg\lambda_{3}$), the flow $\Phi_{{\mathrm{Gauss}},*}$ follows the black curves in the direction of the arrow until it reaches the zero-level of $T_{\mathrm{Gauss}}(x)_{3}$ (slightly wider blue curve, only visible close to the bottom of Fig. 8(a)).

Near-singularity of $T_{\mathrm{Gauss}}$

The zero curve $\{x:T_{\mathrm{Gauss}}(x)_{3}=0\}$ in the $(x_{2},x_{3})$ plane (wide blue in Fig. 8(a)) is given by $\int_{\mathbb{R}}\operatorname{\mathcal{L}}_{\mathrm{Gauss}}(x)(Q)\varphi_{3}(Q)/\varphi_{1}(Q)\mathop{}\!\mathrm{d}Q=0$, where $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}(x)$ is a Gaussian of mean $x_{2}$ and variance $x_{3}$ and $\varphi_{3}(Q)/\varphi_{1}(Q)$ is shown in Fig. 8(b,c). From the profile of $\varphi_{3}/\varphi_{1}$ it is clear that the zero-level forms a single curve connecting the two pieces of the wide blue curve $\{x:T_{\mathrm{Gauss}}(x)_{3}=0\}$ visible in Fig. 8(a). However, this curve has a large radius (passing through the region $x_{3}\gg 1$). For example, there exists a Gaussian $u=\operatorname{\mathcal{L}}_{\mathrm{Gauss}}(x)$ with mean $x_{2}=0$ and large variance $x_{3}$ such that $T_{\mathrm{Gauss}}(x)_{3}=0$, because $\varphi_{3}/\varphi_{1}$ is negative everywhere outside its peak, but the negative values have small modulus (note the scaling of the vertical axis in the zoom of $\varphi_{3}/\varphi_{1}$ in Fig. 8(c)). The fixed point of $x\mapsto\Phi_{{\mathrm{Gauss}},*}(\delta;x)$ (assuming $x_{1}=1$) is the intersection of the two zero-level curves (not visible in Fig. 8(a) as it has large $x_{3}$). The color shading in Fig. 8(a) indicates that the nonlinear transformation $T_{\mathrm{Gauss}}$ is nearly singular close to the line $x_{3}=0$, because the $\mathbb{L}^{2}$-adjoint modes $\varphi_{2}/\varphi_{1}$ and $\varphi_{3}/\varphi_{1}$ are both nearly constant away from the region around $Q\in[-2,2]$ (see Fig. 6, right panel) such that, when inverting $T_{\mathrm{Gauss}}$, the mean $x_{2}$ is very sensitive for small changes in the coefficients for the $\mathbb{L}^{2}$-adjoint modes $\varphi_{2}/\varphi_{1}$ and $\varphi_{3}/\varphi_{1}$ .

Components of error $\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}-\Phi_{{\mathrm{Gauss}},*}$

We perform a detailed convergence analysis along the example trajectory of the exact flow $\Phi_{{\mathrm{Gauss}},*}$ shown in Fig. 8(a): $y_{*}=\Phi_{{\mathrm{Gauss}},*}(\delta;x)$, where $x=(1,0.5,2)^{T}$ and $\delta=0.1$ (thus, $y_{*}\approx(1,0.8459,6.4556)^{T}$).

To understand the factors entering the practically achievable lower limit of the error $\|y_{t_{\mathrm{skip}}}-y_{*}\|=\|\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}(\delta;x)-\Phi_{{\mathrm{Gauss}},*}(\delta;x)\|$, we consider again the identity Eq. 21 used in the proof of Theorem 3.1:

[TABLE]

but re-arrange it using the concrete expressions for $P_{{\mathrm{Gauss}},*}$ and $P_{\mathrm{Gauss}}$:

[TABLE]

Since the matrices $R_{d}$ and $M_{d}({t_{\mathrm{skip}}})$ are invertible, we can apply their inverses to both sides in Eq. 59. For a general distribution $\rho$, the composition of $R_{d}^{-1}$ and $\operatorname{\mathcal{R}}$

[TABLE]

is a projection onto the slow manifold ${\cal C}$ in the coordinates $(\varphi_{1},\varphi_{2},\varphi_{3})$. Furthermore, the nonlinear map $T_{\mathrm{Gauss}}$ is locally invertible in $y_{*}$ (and, hence, also in $y_{t_{\mathrm{skip}}}$, if $y_{t_{\mathrm{skip}}}$ is near $y_{*}$). Its Jacobian is invertible in $y_{*}$ with a moderate norm of its inverse $\|[\partial T_{\mathrm{Gauss}}(y_{*})]^{-1}\|\approx 10$ for the chosen $y_{*}$. Hence, the identity Eq. 59 can be written in the form

[TABLE]

The two residual terms on the right-hand side, labelled $\operatorname{res}(y_{t_{\mathrm{skip}}})$ and $\operatorname{res}_{\delta}(x)$, are the two contributions to the error, before it gets amplified by a moderate factor ($\|[\partial T_{\mathrm{Gauss}}(y_{*})]^{-1}\|\approx 10$) when inverting $T_{\mathrm{Gauss}}$. The spectral properties of the flow $M$ ensure that

[TABLE]

where $d_{\mathrm{tr}}=-\lambda_{4}\approx 10.3$. Applying this estimate to $\eta=y_{t_{\mathrm{skip}}}$ and $t={t_{\mathrm{skip}}}$, and to $\eta=x$ and $t={t_{\mathrm{skip}}}+\delta$ gives the asymptotics $\sim\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$ in ${t_{\mathrm{skip}}}$ for $M_{d}({t_{\mathrm{skip}}})\operatorname{res}(y_{t_{\mathrm{skip}}})$ and $M_{d}({t_{\mathrm{skip}}})\operatorname{res}_{\delta}(x)$, shown in Fig. 9 (red and blue curves with circles). The healed residuals $M_{d}({t_{\mathrm{skip}}})\operatorname{res}(y_{t_{\mathrm{skip}}})$ and $M_{d}({t_{\mathrm{skip}}})\operatorname{res}_{\delta}(x)$ indeed decay with rate $d_{\mathrm{tr}}$ until computational errors for computing the distributions dominate (in this case $10^{-8}$). The matrix $M_{d}({t_{\mathrm{skip}}})^{-1}=M_{d}(-{t_{\mathrm{skip}}})$ has norm of order $\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$ (where $d_{\mathrm{tan}}=-\lambda_{3}\approx 5.71$; see grey dashed line sloping upward in Fig. 9) such that the residuals $\operatorname{res}(y_{t_{\mathrm{skip}}})$ and $\operatorname{res}_{\delta}(x)$ are of order $\sim\exp((d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ (blue and red curves with $+$ marks in Fig. 9). The residuals indeed decrease with rate $d_{\mathrm{tr}}-d_{\mathrm{tan}}$ for increasing ${t_{\mathrm{skip}}}$ until the amplification of the computational errors by $\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}})$ starts to dominate (at ${t_{\mathrm{skip}}}\approx 2$). The true error $y_{t_{\mathrm{skip}}}-y_{*}$ (shown in black in Fig. 9) is then amplified approximately by the norm of $\|[\partial T_{\mathrm{Gauss}}(y_{*})]^{-1}\|\approx 10$, because the residuals $\operatorname{res}$ and $\operatorname{res}_{\delta}$ occur on the manifold ${\cal C}$ (in the coordinates $(\varphi_{1},\varphi_{2},\varphi_{3})$), while the error $y_{t_{\mathrm{skip}}}-y_{*}$ is defined in $\operatorname{dom}\operatorname{\mathcal{L}}$. The relation between the error $y_{t_{\mathrm{skip}}}-y_{*}$ and the residual errors is independent of ${t_{\mathrm{skip}}}$. Overall, the error $y_{t_{\mathrm{skip}}}-y_{*}$ decays with rate $d_{\mathrm{tr}}-d_{\mathrm{tan}}$ asymptotically for increasing ${t_{\mathrm{skip}}}$, but the computational error grows with rate $d_{\mathrm{tan}}$. The optimal healing time ${t_{\mathrm{skip}}}$ is when both errors are of the same order of magnitude.

The identity Eq. 60 becomes a nonlinear fixed-point problem after applying $T_{\mathrm{Gauss}}^{-1}$, for which the right-hand side is a contraction for sufficiently large ${t_{\mathrm{skip}}}$ (see the proof of Theorem 3.1). For Fig. 9 we applied this fixed-point iteration. The final fixed-point iteration correction (shown as a yellow curve in Fig. 9 is always smaller than the error $y_{t_{\mathrm{skip}}}-y_{*}$.

5.4 The size of

computational errors in ensemble computations

The results shown in Figs. 7(a) and 9 show the qualitative behaviour of implicit lifting for increasing ${t_{\mathrm{skip}}}$. Two sources contribute to the overall error. One source is the mismatch between the trajectory started from the lifted point and the projected (along the stable fiber) trajectory on the slow manifold. The size of this contribution is estimated in Theorem 3.1 as decaying with rate $d_{\mathrm{tr}}-d_{\mathrm{tan}}$ with increasing ${t_{\mathrm{skip}}}$ (also observed in Figs. 7(a) and 9). The other source is the limited accuracy in the computations of the lifting $\operatorname{\mathcal{L}}$, the microscopic flow $M$ and the restriction $\operatorname{\mathcal{R}}$. Errors introduced from this limited accuracy grow with rate $d_{\mathrm{tan}}$ for increasing ${t_{\mathrm{skip}}}$. The analysis in Figs. 7(a) and 9 illustrates the trade-off between these two sources of error when the computational error is small ($\approx 10^{-8}$, using chebfun [10, 14]).

If the microscopic flow $M$ describes a multi-particle or high-dimensional stochastic system and is estimated using ensembles of realizations then the computational error of the flow estimate (and, possibly, the computation of $\operatorname{\mathcal{L}}$ and $\operatorname{\mathcal{R}}$) is determined by the ensemble size $N$. This error decreases asymptotically like $1/\sqrt{N}$ for increasing $N$, unless one is able to apply variance reduction techniques (see, for example, [3] for a technique to reduce noise in the computations of Jacobians needed to solve nonlinear systems). In this section we demonstrate that the error behavior can be expected to be qualitatively the same as in Figs. 7(a) and 9, but with stricter limitations on ${t_{\mathrm{skip}}}$ due to larger computational errors in $\operatorname{\mathcal{L}}$, $\operatorname{\mathcal{R}}$, and the flow $M$. To keep the computations simple and comparable to the previous subsection, we perform a Monte-Carlo simulation directly for the SDE Eq. 30.

Figure 10 shows the overall behaviour of the error when performing computations based on random ensembles of finite size $N$, using the lifting operator $\operatorname{\mathcal{L}}_{\mathrm{Gauss}}$, based on Gaussians. For an ensemble size $N$, mean $\bar{Q}$ and we create a random set of initial conditions

[TABLE]

where $\eta\thicksim\mathcal{N}(0,1)$ is a random variable drawn from a standard normal distribution for each $n$. An ensemble of $N$ realizations at positions $Q_{n}$ is restricted according to

[TABLE]

Similar, to the definitions Eq. 40 and Eq. 42, the first component of the argument to $\operatorname{\mathcal{L}}$ and of the output of $\operatorname{\mathcal{R}}$ is the number of realizations, which is preserved. In order to solve Eq. 30 numerically, we use the Euler-Maruyama scheme

[TABLE]

where $h=0.01$ is the step size, $f=-\partial_{Q}V$ and $\xi_{n}\thicksim\mathcal{N}(0,1)$ is standard normal random noise that is uncorrelated, i. e.  $\langle\xi_{n}(t)\xi_{n}(t^{\prime})\rangle=\delta(t-t^{\prime})$.

The error for each ${t_{\mathrm{skip}}}$ in Fig. 10 was estimated by comparing the value of $\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}(\delta;x)$ to the value $\Phi_{{\mathrm{Gauss}},t_{\max}}(\delta;x)$ for the largest ${t_{\mathrm{skip}}}$ (called $t_{\max}$, equalling $1$). Thus, the value of ${t_{\mathrm{skip}}}$ at which the error starts to grow and the growth rate may not have been captured accurately. However, we observe an exponential decay with increasing ${t_{\mathrm{skip}}}$ over approximately two orders of magnitude and the more stringent limitation on ${t_{\mathrm{skip}}}$, as the error stops decreasing at ${t_{\mathrm{skip}}}\approx 0.6$.

Two problems limit the computational accuracy of function evaluations.

In Monte Carlo simulations with ensemble size $N$ the evaluation of the macroscopic lift-evolve-restrict map $P(t;\cdot)$ of the dynamics is noisy in $(\bar{Q},\operatorname{var}Q)$. This is due to the inherent noise in Eq. 30 and due to the noise in the lifting procedure Eq. 61. Hence the evaluation of $P$ with the same input parameters might yield different outputs. The result of $P$ is a random variable with an ensemble-dependent distribution (see Fig. 10, where the distribution of the second component (the mean) of $P(1;(N,-0.5,0.2))$ is shown for a range of $N$). The standard deviation of $P$ decreases with the ensemble size like $\thicksim 1/\sqrt{N}$. 2. 2.

Function evaluations for large $\operatorname{var}Q$ become computationally difficult since a large $\operatorname{var}Q$ implies sampling of trajectories far away from the minima of the potential. Since the potential is steep away from the minima, the drift forces $V^{\prime}$ become large, which results in stability problems of the numerical scheme Eq. 63 for a fixed step size $h$.

When solving $P({t_{\mathrm{skip}}};y)=P({t_{\mathrm{skip}}}+\delta;x)$ for $y$ in the analysis in Fig. 10 we use a Newton iteration with damping $\gamma=0.5$ on the macroscopic level with tolerance $\texttt{tol}=5\cdot 10^{-2}$ where Jacobians are computed by a central finite-difference scheme with $\Delta\bar{Q}=\Delta\operatorname{var}Q=5\cdot 10^{-2}$. The ensemble size is $N=10^{7}$. The level of the minimal error is limited by the finite ensemble size $N$ and the accuracy of function evaluations and approximations of the Jacobian in the Newton iterations (see [3] how the accuracy of the Jacobians can be improved).

6 Discussion

6.1 General estimate for the influence of evaluation errors

While the theoretical convergence result in Theorem 3.1 appears to suggest that a larger ${t_{\mathrm{skip}}}$ always leads to a smaller error, the demonstrations for the Michaelis-Menten kinetics model in Section 4 and the SDE in Section 5 illustrate that there is a trade-off and, hence, an optimal value for ${t_{\mathrm{skip}}}$ in practice. One source for the difference between the estimates of Theorem 3.1 and numerical observations are numerical errors in the evaluation of lifting $\operatorname{\mathcal{L}}$, evolution $M(t;\cdot)$ and restriction $\operatorname{\mathcal{R}}$. The effect of these errors grow along trajectories inside the slow manifold ${\cal C}$ if the vector field tangent to ${\cal C}$ has non-zero expansion rates forward or backward in time. This becomes clear when looking at the arguments in the proof of Theorem 3.1. The approximate solution $y_{t_{\mathrm{skip}}}$ is the fixed point of the map (see Equation Eq. 82)

[TABLE]

According to Theorem 3.1, $y_{t_{\mathrm{skip}}}-y_{*}\sim\exp((d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$, where $d_{\mathrm{tan}}$ is defined as $\max\{d_{\mathrm{tan}}^{+},d_{\mathrm{tan}}^{-}\}$, the maximum of the forward ($d_{\mathrm{tan}}^{+}$) and backward ($d_{\mathrm{tan}}^{-}$) expansion rate of the flow $M|_{\cal C}$ tangential to ${\cal C}$, and $d_{\mathrm{tr}}$ is the rate of attraction transversal to ${\cal C}$. However, if we take into account evaluation errors, we have to distinguish between the exact and approximate operators. That is, $P_{*}(t;\cdot)$ equals $\operatorname{\mathcal{R}}\circ M(t;\cdot)|_{\cal C}\circ g\circ\operatorname{\mathcal{L}}$ (recall that $g$ is the stable fiber projection) and $P(t;\cdot)$ equals $\operatorname{\mathcal{R}}_{\Delta}\circ M_{\Delta}(t;\cdot)\circ\operatorname{\mathcal{L}}_{\Delta}$ where we use the subscript $\Delta$ to indicate that the operator is affected by small errors. For $\operatorname{\mathcal{L}}_{\Delta}$ and $\operatorname{\mathcal{R}}_{\Delta}$ this means simply that they are perturbations of $\operatorname{\mathcal{L}}$ and $\operatorname{\mathcal{R}}$ of size $\Delta$. The evaluation error in $M$ along trajectories in ${\cal C}$ causes errors of size

[TABLE]

These errors in $\operatorname{\mathcal{L}}_{\Delta}$, $\operatorname{\mathcal{R}}_{\Delta}$ and $M_{\Delta}(t;\cdot)$ are all part of the term $P({t_{\mathrm{skip}}}+\delta;x)$ in Eq. 64) such that the error grows for increasing ${t_{\mathrm{skip}}}$ at the rate

[TABLE]

which gets then amplified by the expansion rate of $M(-{t_{\mathrm{skip}}};\cdot)|_{\cal C}$ when applying $P_{*}({t_{\mathrm{skip}}};\cdot)^{-1}$. Thus, there will be an error between the exact fixed point of the map Eq. 64 and the fixed point with evaluation errors. This error is of order $\Delta\exp((d_{\mathrm{tan}}^{+}+d_{\mathrm{tan}}^{-}){t_{\mathrm{skip}}})$, which is growing exponentially in ${t_{\mathrm{skip}}}$. This is visible in all computational results:

• •

In the Michaelis-Menten kinetics model in Section 4 the error $\Delta$ is of order $10^{-10}$ and $d_{\mathrm{tan}}$ is of order $\varepsilon$ (which is $10^{-2}$) such that the growth of the error with ${t_{\mathrm{skip}}}$ is not visible in the range of ${t_{\mathrm{skip}}}$ between [math] and $30$ in Figs. 4(e) and 4(f).

• •

For the stochastic differential equation in Section 5, $d_{\mathrm{tan}}^{+}$ is zero and $d_{\mathrm{tan}}^{-}=-\lambda_{3}\approx 5.71$. For Figs. 7(a) and 9 we computed the evolution of densities directly using the Fokker-Plank equation and chebfun such that the evaluation error $\Delta$ is of the order $10^{-8}$ (visible as the lower bound on the residuals $\|P_{\mathrm{lin}}({t_{\mathrm{skip}}})-P_{{\mathrm{lin}},*}({t_{\mathrm{skip}}})\|$ in Fig. 7(a) and in the residuals after healing in Fig. 9). Thus, the overall influence of the evaluation error is of order $\Delta\exp({t_{\mathrm{skip}}}d_{\mathrm{tan}}^{-})$. The amplification factor reaches $\sim 10^{5}$ for ${t_{\mathrm{skip}}}=2$. In Fig. 7(a) evaluation errors dominate only from ${t_{\mathrm{skip}}}\approx 3$, while in Fig. 9 they dominate from ${t_{\mathrm{skip}}}\approx 2.5$.

• •

In Fig. 10 the growth rate of the evaluation error is the same as in Fig. 9, but the basic evaluation error of a single time step of $M_{\Delta}(t,\cdot)$ and the lifting $\operatorname{\mathcal{L}}_{\Delta}$ is larger (as they are generated from ensembles): $\Delta\sim 10^{-3.5}$ for ensemble size $N=10^{7}$. Thus, the effects of evaluation error start to dominate already for ${t_{\mathrm{skip}}}\approx 0.7$. With smaller, more realistic, ensemble sizes the restriction on ${t_{\mathrm{skip}}}$ posed by evaluation errors will be even more severe. Since the necessary length of ${t_{\mathrm{skip}}}$ to reduce the projection error $y_{t_{\mathrm{skip}}}-y_{*}$ (from Theorem 3.1) is dictated by $d_{\mathrm{tan}}^{-}-d_{\mathrm{tr}}$, we have a general approximate optimal healing time for positive evaluation errors $\Delta$ of the order

[TABLE]

resulting in an optimal error of the order

[TABLE]

In the limit of large time scale separation ($d_{\mathrm{tan}}^{\pm}/d_{\mathrm{tr}}\to 0$) the power $p$ of the error reaches $1$ and the optimal ${t_{\mathrm{skip}}}$ is of order $-\log\Delta/d_{\mathrm{tr}}$.

6.2 Consequences for equation-free analysis of stochastic

systems

The lift-evolve-restrict map $P_{\mathrm{Gauss}}(t;\cdot)$ in Section 5 reduced the SDE $\mathop{}\!\mathrm{d}Q=-V^{\prime}(Q)+\sigma\mathop{}\!\mathrm{d}W_{t}$ (or, more precisely, its Fokker-Planck equation) to the slow manifold (a linear subspace) spanned by its first $3$ modes. Barkley et al [5] observed that the map $P_{\mathrm{Gauss}}(t;\cdot)$ (called moment map in [5]) is nonlinear and, hence, suspected that the nonlinearity of $P_{\mathrm{Gauss}}$ may be the object of interest for nonlinear analysis (such as finding multiple equilibria, bifurcations under parameter changes, etc). However, as equation Eq. 57 shows, the exact flow map $\Phi_{{\mathrm{Gauss}},*}(\delta;\cdot)$ of the low-order moments is still a nonlinear transformation (by $T_{\mathrm{Gauss}}$) of a linear map such that there is no nonlinear dynamic behaviour present. More precisely, the phase portrait of the exact flow map $\Phi_{{\mathrm{Gauss}},*}(\delta;\cdot)$ is topologically conjugate to the phase portrait of a linear system. Since the approximate flow $\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}(\delta;\cdot)$, computed with $P_{\mathrm{Gauss}}$, converges to $\Phi_{{\mathrm{Gauss}},*}(\delta;\cdot)$ for ${t_{\mathrm{skip}}}\to\infty$ we do not expect nonlinear behavior for $\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}$ either.

This raises the question what the natural nonlinearity of the underlying system is in the case of equation-free methods applied to stochastic systems.

6.2.1 Artificial nonlinearity

Since the Fokker-Planck equation is linear, the apparent nonlinear dynamics arises only due to artificial projections of nonlinearly transformed phase portraits of the linear Fokker-Planck equation when the healing time ${t_{\mathrm{skip}}}$ is not sufficiently large. For example, let us consider again the SDE with lifting to a Gaussian distribution from Section 5.3. What happens if we choose a moment map for only the zeroth and first moment but an insufficiently large ${t_{\mathrm{skip}}}$ (which would have to be $\sim 1/\lambda_{2}\approx 10^{6}$ to make Theorem 3.1 applicable)? For illustration we choose a lifting to near-delta Gaussian distributions, similar to [5]. In the notation from Section 5.3 this means that we keep $x_{1}$ equal to $1$ (mass), vary $x_{2}$ (mean) between $-3$ and $3$, and keep $x_{3}\ll 1$ (variance) fixed ($x_{3}=0.04$ for the illustration in Fig. 11). The restriction is then the projection on the zeroth and first moment. If the healing time ${t_{\mathrm{skip}}}$ satisfies $1/\lambda_{3}\ll{t_{\mathrm{skip}}}\ll 1/\lambda_{2}$ (instead of ${t_{\mathrm{skip}}}\sim 1/\lambda_{2}$), then we obtain for the approximate flow $\Phi_{{\mathrm{Gauss}},{t_{\mathrm{skip}}}}$ a projection of the phase portrait Fig. 8(a) onto the line with $x_{3}=0.04$. Figure 11 shows this projected phase portrait (arrows on the $x$-axis) and the associated right-hand side (in blue). It resembles a phase portrait of a scalar ODE with two coexisting stable fixed points, separated by an unstable fixed point. Of course, this nonlinearity is created artificially by projecting the accurate nonlinearly transformed two-dimensional phase portrait of a linear system onto an arbitrarily chosen line in $\mathbb{R}^{2}$.

6.2.2 Reduction of high-dimensional SDEs

While in high-dimensional SDEs there is at first sight no obvious nonlinearity present in the evolution of densities (see Fokker-Planck equation Eq. 32), the reduction to low-order moments of a multi-particle system with randomness still gives a valid dimension reduction procedure. We give an informal outline of the argument for a particularly simple case in which dimension reduction is in theory possible according to Givon et al. [17] (see also textbook [36]).

Let us assume that the simulation (say, an agent-based simulation) can be modelled by a high-dimensional SDE (which is the microscopic model)

[TABLE]

where $u\in\mathbb{R}^{n_{u}}$ and (to keep the argument simple) $\sigma_{u}$ is constant and regular, and $W_{u,t}$ are $n_{u}$ independent instances of Brownian motion. Let us also assume that there exist coordinates $(x,y)\in\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{y}}$ ($n_{x}+n_{y}=n_{u}$) for $u$ such that in these coordinates we have a time scale separation:

[TABLE]

and that for each $x$ the random variable $y$ converges to its stationary density with rate of order $1$ (fast). Let $v_{0}(x,y)$ be the nullvector of the Fokker-Planck operator of the fast subsystem of Eq. 66, $p\mapsto L_{0}p=\partial_{y}[\frac{1}{2}\sigma_{y}^{T}\sigma_{y}\partial_{y}p-gp]$, with $\int v_{0}(x,y)\mathop{}\!\mathrm{d}y=1$. Any function of the form $v_{0}(x,y)p_{x}(x)$ is also a nullvector of $L_{0}$. If $(\varepsilon\lambda,p)$ (with $O(\lambda)=1$) is an eigenpair of the Fokker-Planck operator $L_{0}+\varepsilon L_{1}$ with $L_{1}p=\partial_{x}[\frac{1}{2}\sigma_{x}^{T}\sigma_{x}\partial_{x}p-fp]$ for the combined system Eq. 66 in $(x,y)$ coordinates, then $\lambda=\lambda_{0}+O(\varepsilon)$, $p(x,y)=v_{0}(x,y)p_{x}(x)+O(\varepsilon)$, where $(\varepsilon\lambda_{0},p_{x})$ is an eigenpair of the of the right-hand side of the Fokker-Planck equation for the reduced SDE

[TABLE]

In Eq. 67 $\tilde{f}(x)=\int f(x,y)v_{0}(x,y)\mathop{}\!\mathrm{d}y$ is the conditional expectation with respect to $x$ of the drift in $x$ and $\tilde{\sigma}_{x}(x)=\sigma_{x}\left[\int v_{0}(x,y)\mathop{}\!\mathrm{d}y/2\right]^{1/2}$ is the standard deviation of $x$ in the stationary distribution of $y$. Consequently, performing equation-free analysis on the high-dimensional SDE Eq. 66 using a small number $d$ of variables gives the same results as equation-free analysis on the reduced system Eq. 67 (up to order $\varepsilon^{2}$).

Givon et al. [17] discuss dimension reduction more generally (independent of explicit spatial coordinates $x$ and $y$) for Fokker-Planck operators of the form $L_{0}+\varepsilon L_{1}$, assuming that the linear operator $L_{0}$ has a non-trivial kernel (dimension greater than $1$, implying that $\varepsilon$ is a singular perturbation parameter). Hence, equation-free-analysis based on implicit lifting and sufficiently large healing times can be used to perform closure-on-demand, as described in [25], rigorously. Convergence of the approximate system created by lift-evolve-restrict maps to the Fokker-Planck operator of the reduced system Eq. 67 occurs in the sense of classical singular perturbation theory toward an attracting low-dimensional linear invariant subspace of densities in the domain of definition of $L_{0}+\varepsilon L_{1}$, as ensured by Theorem 3.1 for sufficiently large healing times ${t_{\mathrm{skip}}}$.

For the case that the high-dimensional SDE consists of a large number $N$ of random variables (for example, describing agents) our analysis in Section 5 and the above discussion raise an important point. Applying the equation-free procedure to initial densities and the Fokker-Planck operator $L_{0}+\varepsilon L_{1}$ does not reduce the high-dimensional SDE to a low-dimensional SDE, but it reduces the high-dimensional SDE to a low-dimensional linear ODE for the coefficients of the leading modes of the Fokker-Planck equation. Hence, increasing the number of variables $N$ (e.g., agents) does not increase the spectral gap or the time scale separation. This is obvious for the simple SDE example in Section 5: decreasing the noise level will let $\lambda_{2}/\lambda_{3}$ converge to [math] (the time scale for escape from one well to the other), but $\lambda_{3}/\lambda_{4}$ will remain approximately $1/2$. Hence, we need the convergence result for finite time scale separation to prove validity of the model reduction. Results for sufficiently large time scale separation such as those by Zagaris et al. [16, 48, 49] (using, for example, constrained runs) and Marschler et al. [31] are not applicable to equation-free methods operating on Fokker-Planck equations, if the aim is to extract the decay rate or shape of the dominant modes of the Fokker-Planck equation.

In summary, one possible work flow for analysing a high-dimensional SDE with generator splittable as $L_{0}+\varepsilon L_{1}$ with equation-free methods is: (1) use the equation-free moment map to determine properties of the leading $d$ eigenmodes $\varphi_{j}$ and eigenvalues $\lambda_{j}$ of $L_{0}+\varepsilon L_{1}$; (2) if these $\varphi_{j}$ and $\lambda_{j}$ are also the leading eigenmodes and eigenvalues to an operator $L_{1}$ for a Fokker-Planck equation of a low-dimensional SDE, identify the properties of $L_{1}$ from the modes (for example, singular points of the potential).

7 Outlook

The arguments in Section 5, studying the simple scalar SDE $\mathop{}\!\mathrm{d}Q=-V^{\prime}(Q)\mathop{}\!\mathrm{d}t+\mathop{}\!\mathrm{d}W_{t}$, and the discussion in Section 6.2 treat SDEs as linear evolution equations for densities. The sections below outline how one may have to modify the arguments of Theorem 3.1 for other tasks of equation-free analysis, which are beyond the scope of this paper.

7.1 Bifurcation analysis for the drift of the

reduced system

Assume that we have access to a simulator of a system that can be modelled by a high-dimensional SDEs of type Eq. 65,

[TABLE]

with time-scale separation as in Eq. 66. A sensible object for nonlinear equation-free analysis is a bifurcation analysis of the deterministic part $\dot{x}=\tilde{f}(x)$ of the reduced SDE Eq. 67,

[TABLE]

assuming a reduction as discussed in Section 6.2.2 is possible. For example, one may want to determine its phase portraits and their parameter dependence. If one had a direct simulator of the low-dimensional reduced SDE Eq. 69, one could approximate $\tilde{f}$ in any given $x_{0}\in\mathbb{R}^{n_{x}}$ via

[TABLE]

where $X_{\delta}$ (a random variable in $\mathbb{R}^{n_{x}}$) is the solution of the SDE Eq. 69 at time $\delta$ starting from the deterministic $x_{0}$, and $EX_{\delta}\in\mathbb{R}^{n_{x}}$ is its expectation.

Equation-free analysis based on a lift-evolve-restrict map $P$ with healing time provides an approximation for Eq. 70 if only a simulation of the high-dimensional SDE Eq. 68 is available. The healing time permits the fast variable $y$ to settle to its stationary density $v_{0}(x,y)$ before one measures $\tilde{f}$. Since the slow-fast coordinate split of $u$ into $x$ and $y$ is unknown, one has to define a lifting $\operatorname{\mathcal{L}}$ and a restriction $\operatorname{\mathcal{R}}$ between $\mathbb{R}^{n_{x}}$ and the space of random variables $U$ in $\mathbb{R}^{n_{u}}$.

Let us assume that the lift $\operatorname{\mathcal{L}}(x_{L})$ of $x_{L}\in\mathbb{R}^{n_{x}}$ is a random variable $U_{0}$ in $\mathbb{R}^{n_{u}}$ with density $p_{0}$ on $\mathbb{R}^{n_{u}}$. The SDE Eq. 68 creates a Markov process $t\mapsto U_{t}$ for $t\geq 0$. Let us consider a restriction $\operatorname{\mathcal{R}}$ of a random variable $U_{t}$ that is the expectation $ER(U_{t})$ of a map $R:\mathbb{R}^{n_{u}}\mapsto\mathbb{R}^{n_{x}}$. Thus, the lift-evolve-restrict map $P:\mathbb{R}\times\mathbb{R}^{n_{x}}\mapsto\mathbb{R}^{n_{x}}$ is $P(t;x_{L})=E[R(U_{t})|U_{0}=\operatorname{\mathcal{L}}(x_{L})]$. A good approximation of the deterministic part of the slow flow (in $x_{L}$ coordinates) would not be $(y-x_{0})/\delta$ where $y$ is the solution of $P({t_{\mathrm{skip}}}+\delta;x_{0})=P({t_{\mathrm{skip}}};y)$. Rather, a possible construction is to define $x_{R}=P({t_{\mathrm{skip}}};x_{L})$ and then compute

[TABLE]

This means that one first solves the SDE for the healing time ${t_{\mathrm{skip}}}$, then increases time to ${t_{\mathrm{skip}}}+\delta$, and uses the conditional expectation of $R(U_{{t_{\mathrm{skip}}}+\delta})$, with the condition that $R(U_{t_{\mathrm{skip}}})=x_{R}$. This conditional expectation enters the difference quotient for $\tilde{f}_{L}(x_{L})$, which is otherwise similar to Eq. 70. Constructions of the form Eq. 71 do not fit into the framework of Theorem 3.1. Still, we conjecture that the function $\tilde{f}_{L}$ approximates $\tilde{f}$ (up to a coordinate change from $x_{L}$ to $x$) for sufficiently small $\delta$ and large ${t_{\mathrm{skip}}}$. The approximation will become accurate only in the limit of large time scale separation for a set of $n_{x}$ slow variables (in contrast to Theorem 3.1), but we need only genericity conditions on $\operatorname{\mathcal{L}}$ and $\operatorname{\mathcal{R}}$.

7.2 Averaging deterministic high-dimensional systems

There is still another gap to applications for multi-particle systems, which are commonly deterministic at the microscopic level. For example, Barkley et al. [5] used the scalar SDE Eq. 30 as a simple model for a heat bath problem where the position $Q$ of a heavy particle of mass $M$ and generalized coordinates $(Q,P)$ is coupled to a heat bath of $N$ smaller particles of masses $m_{i}$ and generalized coordinates $(q_{i},p_{i})$ for $i=1,\ldots,N$. The full system in [5] was described by the Hamiltonian

[TABLE]

where the number $N$ of particles is large and the masses $m_{i}$ and spring coupling constants $k_{i}$ are small (with particular $N$-dependent distributions, see [5], eq. (2.2)). The necessary assumption to enable treatment of a fast deterministic subsystem as a stochastic system is some form of ergodicity: any distribution of initial conditions of the fast subsystem converges rapidly to a unique stationary distribution (conditioned on the slow variables). This condition is hard to verify (even empirically) for any particular system. In particular, it is not true for Eq. 72 if one treats the coordinates $(Q,P)$ as the slow variables since the small masses are only coupled through the heavy particle. Convergence to an SDE is only guaranteed for the system with Hamiltonian Eq. 72 if the initial conditions for $q_{i}$ and $p_{i}$ are set according to the stationary measure conditioned on $P$ and $Q$ (which was done in [5], see [5, 27, 37] for background results). Hence, the introduction of a healing time ${t_{\mathrm{skip}}}$ will not have an improving effect for equation-free analysis of the heat bath problem Eq. 72.

7.3 Approximation of stochastic slow manifolds

As mentioned already in the introduction, our convergence result for finite time scale separation relies on a result about model reduction that is valid for finite time scale separation, namely the persistence of normally hyperbolic invariant manifolds and their stable fibers. While the model reduction results for stochastic systems in [17, 36] provide only statements for the limit of infinite time scale separation, stronger results are available for stochastic systems, if one is able to fix the noise realization (for example, the Brownian path) [1, 2]. In this case, the microscopic map $M$ has, for the example of an SDE of the form $\mathop{}\!\mathrm{d}u=F(u)\mathop{}\!\mathrm{d}t+\sigma_{u}\mathop{}\!\mathrm{d}W_{u,t}$, the form $M(t;u,\omega)$, where $\omega\in C([0,\infty);\mathbb{R}^{D})$ is a realization of the Wiener process $W_{u,t}$ and $M$ satisfies the invariance relation $M(t+s;x,\omega)=M(t;M(s;x,\omega),\omega(s+\cdot))$.

Invariant stochastic manifolds ${\cal C}$ are then invariant objects depending on the realization (one may write ${\cal C}(\omega)$). Their persistence and attraction properties have been proven for some cases such as finite-dimensional SDEs [9, 47] and SPDEs [11]. For these cases, an implicit equation-free scheme $y=\Phi_{t_{\mathrm{skip}}}(\delta;x,\omega)$ defined implicitly via

[TABLE]

may converge in a similar way as claimed in Theorem 3.1. However, the stochastic invariant manifold results and the implementation of Eq. 73 depend on the ability to use the same realization $\omega$ throughout the computation, as was done in [22] (for example, for different arguments $y$ during a Newton iteration for Eq. 73). While fixing the realization is possible for SDEs, for many of the applications for equation-free analysis [18, 24, 30, 32, 33, 42, 45] it is not clear how to do that.

8 Conclusion

This paper proves convergence of equation-free methods, based on lift-evolve-restrict maps $P(t;\cdot)=\operatorname{\mathcal{R}}\cdot M(t;\cdot)\circ\operatorname{\mathcal{L}}$. Our convergence proof does not assume that the time scale separation becomes large, in contrast to previous results [49, 31]. Rather, convergence is achieved for finite time scale separation, but in the limit of large healing time ${t_{\mathrm{skip}}}$ and an implicit approximation of the slow flow $\Phi_{*}(t;x)$: $P({t_{\mathrm{skip}}};y)=P(t+{t_{\mathrm{skip}}};x)$ defines the approximation $\Phi_{t_{\mathrm{skip}}}(t;x):=y$. The original explicit equation-free framework, as proposed by Kevrekidis et al., corresponds to the case where ${t_{\mathrm{skip}}}=0$ and $\operatorname{\mathcal{R}}\circ\operatorname{\mathcal{L}}=I$. The analysis is performed for attracting slow manifolds in deterministic systems. However, we demonstrate on a simple SDE that our result may also be useful for stochastic systems, where the time scale separation is in the spectrum of the Fokker-Planck equation and is often only of order $1$. In particular, for the prototype example investigated by [5] the implicit flow approximation $\Phi_{t_{\mathrm{skip}}}$ converges to the true solution $\Phi_{*}$ of the linear Fokker-Planck equation for large healing times ${t_{\mathrm{skip}}}$.

Acknowledgements

J. Sieber’s research was supported by funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement number 643073, by the EPSRC Centre for Predictive Modelling in Healthcare (Grant Number EP/N014391/1) and by the EPSRC Fellowship EP/N023544/1.

C. Marschler and J. Starke would like to thank Civilingeniør Frederik Christiansens Almennyttige Fond for financial support. J. Starke would also like to thank the Villum Fonden (VKR-Centre of Excellence Ocean Life), the Technical University of Denmark and Queen Mary University of London for financial support.

Appendix A Proof of Convergence Theorem 3.1

For the proof of Theorem 3.1 we have to analyze the two equations (for $y$ and $y_{*}$ respectively)

[TABLE]

In both equations $x\in\mathbb{R}^{d}$ enters as a parameter. Section 3 ensures that the solution $y_{*}$ of Eq. 75 is unique and independent of ${t_{\mathrm{skip}}}$. For equation Eq. 74 we have to prove the existence of a solution $y$, and prove that it is close to $y_{*}$ for sufficiently large ${t_{\mathrm{skip}}}$. Throughout this appendix we will use the notations

[TABLE]

to describe that $\|u(t)\exp(-\alpha t)\|$ is bounded uniformly for all $t\geq 0$, and that the function $v(t)\exp(-\alpha t)$ tends to zero for $t\to\infty$. For the special case $\alpha=0$ we write $O(1)$ and $o(1)$. If the quantity depends also on other parameters (say, $y\in\operatorname{dom}\operatorname{\mathcal{L}}$) then the expression implies uniformity (for example, for $y$ close to $y_{*}$) unless stated explicitly otherwise.

Using the definitions Eq. 12 of $P_{*}(t;x)=\operatorname{\mathcal{R}}(M(t;g(\operatorname{\mathcal{L}}(x))))$ and Eq. 15 for the map $P(t;x)=\operatorname{\mathcal{R}}(M(t;\operatorname{\mathcal{L}}(x)))$, equation Eq. 74 can be written in the form (using Eq. 75)

[TABLE]

The operator $P_{*}$ and the newly introduced $G$ and $H$ satisfy the following conditions on their derivatives by Section 3, Eq. 7 and Eq. 8 on separation of time scales for the flow $M$:

[TABLE]

for all $j\in\{0,\ldots,k_{\max}\}$ and all $y$ in a neighborhood of $y_{*}$. In the case of $H$ the bound is also uniform for $\delta\in[-\delta_{\max},\delta_{\max}]$. Thus, the parameter $\delta$ has been dropped from the list of arguments in $H$. Combining the separation of time scales in Section 3, Eq. 7, with Section 3 on the uniform invertibility of $\operatorname{\mathcal{R}}|_{\cal C}$ and $g\circ\operatorname{\mathcal{L}}:\operatorname{dom}\operatorname{\mathcal{L}}\mapsto{\cal C}$, we have a Lipschitz constant ($C$ is independent of $y_{1}$, $y_{2}$ and ${t_{\mathrm{skip}}}$)

[TABLE]

when inverting $P_{*}({t_{\mathrm{skip}}};\cdot)$ for all $y_{1},y_{2}$ in a neighborhood of $y_{*}$ and all ${t_{\mathrm{skip}}}\geq 0$. We also note that

[TABLE]

Specifically, these derivatives depend only on $\delta\in[-\delta_{\min},\delta_{\max}]$. Thus, $\partial^{j}y_{*}(x)$ are uniformly bounded due to Eq. 7, and because we required $\exp(d_{\mathrm{tan}}\delta_{\max})=O(1)$.

Abbreviating notation In the following all derivatives of the functions $P_{*}$, $G$ and $H$ are with respect to their second argument ($y$ or $x$). The argument ${t_{\mathrm{skip}}}$ enters the functions $P_{*}$, $G$ and $H$ as a parameter that we will drop in our notation such that we will write, for example, $\partial^{3}P_{*}(y_{*})[\partial y_{*}]^{2}[\partial^{2}y_{*}]$ for $\partial_{2}^{3}P_{*}({t_{\mathrm{skip}}};y_{*})[\partial y_{*}/\partial x]^{2}[\partial^{2}y_{*}/(\partial x)^{2}]$. The parameter ${t_{\mathrm{skip}}}$ enters estimates via the bounds Eq. 77–Eq. 81 for $P_{*}$, $G$ and $H$.

The properties Eq. 77–Eq. 81 make Banach’s contraction mapping principle applicable to Equation Eq. 76 in a sufficiently small neighborhood of $y_{*}$ and for sufficiently large ${t_{\mathrm{skip}}}$ (as shown in the paragraph that follows). We then estimate the error of the derivatives of $y$ with respect to $x$.

Existence of solution $y$ and its error

We apply the Banach contraction mapping principle to the map

[TABLE]

($P_{*}^{-1}(\cdot)$ is the inverse of the diffeomorphism $P_{*}:U(y_{*})\mapsto U(P_{*}(y_{*}))$). Let $B$ be a closed ball around $y_{*}$ of radius $R$ in which all estimates Eq. 77–Eq. 80 on $P_{*}$, $G$ and $H$ hold. Combining the estimate Eq. 80 for the Lipschitz constant of $P_{*}^{-1}$ with $y_{1}=y$ and $y_{2}=y_{*}$, and the bound on the derivatives for $G$ (w.r.t. $y$) gives an estimate for the difference of $N(y)$ from $y_{*}$:

[TABLE]

Thus, choosing ${t_{\mathrm{skip}}}$ sufficiently large, we can ensure that $N$ maps $B$ back into itself (since $d_{\mathrm{tan}}<d_{\mathrm{tr}}$). Similarly, the Lipschitz constant of $N$ in $B$ can be estimated by

[TABLE]

where $C\exp((d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})\max_{B}\|\partial G\|$ is smaller than unity for sufficiently large ${t_{\mathrm{skip}}}$. Consequently, $N$ has a unique fixed point $y$ in $B$, which solves the perturbed problem Eq. 74. Moreover, the difference $y-y_{*}$ satisfies

[TABLE]

Error of derivatives

The smoothness of the coefficients in Eq. 76 ensures that $y$ is also differentiable as a function of $x$ up to order $k_{\max}$. We want to prove that for $\ell$ satisfying $\ell\leq k_{\max}-1$ (where $k_{\max}$ is the order of differentiability of the coefficients in Eq. 76) and $(2\ell+1)d_{\mathrm{tan}}<d_{\mathrm{tr}}$ the bound on the error is

[TABLE]

We prove this by induction starting from $\ell=1$, which we check first using the previous paragraph’s results.

Assume that the bound Eq. 84 holds for all derivatives up to $\ell-1$. This implies, in combination with Eq. 81, that $y$, $\partial y$, …, $\partial^{\ell-1}y$ are bounded uniformly for all ${t_{\mathrm{skip}}}\geq 0$ (just like $\partial^{\ell}y_{*}$ for $\ell=1\ldots k_{\max}$ by Eq. 81). In order to estimate the difference $\partial^{\ell}y-\partial^{\ell}y_{*}$, we return to Eq. 76 and differentiate each of the terms $\ell$ times with respect to $x$ (noting that $y_{*}$ and $y$ are also functions of $x$):

[TABLE]

The term $\partial^{\ell}H(x)$ is $O(1)$ for all ${t_{\mathrm{skip}}}\geq 0$ by Eq. 79. In the term $\partial^{\ell}/(\partial x^{\ell})[G(y)]$ we extract the highest-order derivative of $y$ by writing it in the form

[TABLE]

For Appendix A the boundedness of the $O(1)$ terms follows from the boundedness of all their parts: the derivatives of $G$ are bounded by Eq. 78, $\partial^{\ell}y_{*}$ is bounded by Eq. 81, and $y$, $\partial y$,…, $\partial^{\ell-1}y$ are bounded by induction hypothesis. The pre-factor $\partial G(y)$ of $\partial^{\ell}y-\partial^{\ell}y_{*}$ is also bounded uniformly for all ${t_{\mathrm{skip}}}\geq 0$.

Inserting the right-hand side of Appendix A into the right-hand side of Eq. 85, we obtain

[TABLE]

Expanding the left-hand side of the above equation using the chain rule, we get a sequence of differences with equal powers of derivatives of $P_{*}$, $y$ and $y_{*}$. From this sequence of differences we extract the difference between derivatives involving $\partial^{\ell}y$ and $\partial^{\ell}y_{*}$ and collect all other terms in a remainder $r$ (which is present only for $\ell>1$ and will later turn out to be of order $O(\exp((2\ell d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$):

[TABLE]

From the difference with the highest-order derivatives of $y$ and $y_{*}$ we extract the difference $\partial^{\ell}y-\partial^{\ell}y_{*}$ by adding zeroes. Using the notational convention

[TABLE]

for the mean between two points of a single-argument function $F$ in the following,

[TABLE]

The order $O(\exp((2d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$ of the second term follows from the bounds on $y-y_{*}$ (given in Eq. 83), $\partial^{2}P_{*}$ (given in Eq. 77) and the boundedness of $\partial^{\ell}y_{*}$ (given in Eq. 81). This immediately implies the estimate for the case $\ell=1$: inserting Eq. 90 into Eq. 87, we have for $\ell=1$

[TABLE]

In Eq. 91 we have collected the bounded terms with pre-factors $\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}})$ and $\exp((2d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$ using the larger pre-factor $\exp((2d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}})$. Since (by Eq. 80) the inverse of $\partial P_{*}(y)$ satisfies $\partial P_{*}(y)^{-1}=O(\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}}))$ we can rearrange Eq. 91 to isolate $\partial y-\partial y_{*}$ for large ${t_{\mathrm{skip}}}$, giving the estimate (note that $\partial G(y)=O(1)$)

[TABLE]

which is what we had to prove for $\ell=1$.

Error of higher-order derivatives

Let us assume that the assumptions of the theorem are satisfied for all $j<\ell$ with $\ell\geq 2$. By the conditions of the theorem we assume that $(2\ell+1)d_{\mathrm{tan}}<d_{\mathrm{tr}}$ and the conditions Eq. 77–Eq. 81 are satisfied for $j\leq\ell$ (including existence of the corresponding derivatives).

For $\ell>1$ we have to include the remainder $r$ from Eq. 88 in our consideration. This remainder is a sum of expressions $a_{\nu}$ of the form

[TABLE]

where $2\leq j\leq\ell$, and $\nu$ is a $j$-tuple of integers $\nu_{i}\in\{1,\ldots,\ell-1\}$ with $\sum_{i=1}^{j}\nu_{i}=\ell$. All factors $\partial^{\nu_{i}}y$ and $\partial^{\nu_{i}}y_{*}$ are of order $O(1)$ with respect to ${t_{\mathrm{skip}}}$ according to Eq. 81 and induction hypothesis. The terms $\partial^{j}P_{*}(y)$ and $\partial^{j}P_{*}(y_{*})$ are of order $O(\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}}))$ according to Eq. 77. The difference in Eq. 93 can be expressed as a sum of $j+1$ differences involving $\partial^{i}y-\partial^{i}y_{*}$ for some $i\in\{0\ldots,\ell-1\}$ by adding $j+1$ zeros:

[TABLE]

The right-hand side in Eq. 94 is of order $O(\exp((2d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$. The $i$th term in the sum in Eq. 95 is of order $O(\exp((d_{\mathrm{tan}}(1+(2\nu_{i}+1))-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$. So, since $\nu_{i}\leq\ell-1$ and $\ell>1$, all terms in the sum for $a_{\nu}$ are at most of order $O(\exp((2\ell d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$. Consequently,

[TABLE]

Inserting this estimate in combination with Eq. 88 and Eq. 90 into Eq. 87, we obtain

[TABLE]

In Eq. 97 we have included the smaller error terms $O(\exp((2d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$ and $O(\exp(-d_{\mathrm{tr}}{t_{\mathrm{skip}}}))$ into the (for $\ell>1$) larger $O(\exp((2\ell d_{\mathrm{tan}}-d_{\mathrm{tr}}){t_{\mathrm{skip}}}))$. Since, $d_{\mathrm{tr}}<d_{\mathrm{tan}}$, $\partial G(y)=O(1)$ and $\partial P_{*}(y)=O(\exp(d_{\mathrm{tan}}{t_{\mathrm{skip}}}))$, we can isolate $\partial^{\ell}y-\partial^{\ell}y_{*}$ in Eq. 97. This results in the asymptotic estimate claimed in Theorem 3.1:

[TABLE]

Appendix B Brief description of supplementary material

The supplementary material contains Matlab/octave scripts and functions that reproduce Fig. 4 from Section 4. The provided zip file unpacks into folder \seqsplitdemo_Michaelis_Menten/. The main script is \seqsplitdemo_Michaelis_Menten.m, which will reproduce Fig. 4, showing phase space geometry of the Michaelis-Menten kinetics Eq. 23 with explicit time scale separation as also studied by Gear et al. and others [34, 16, 48, 49].

• •

Folder \seqsplitdemo_Michaelis_Menten/rotated/ contains the published html output from the script for the rotated coordinate system Eq. 29 in file \seqsplitdemo_Michaelis_Menten.html.

• •

Folder demo_Michaelis_Menten/unrotated/ contains the published html output from the script for the coordinate system with explicit time scale separation Eq. 23 in file demo_Michaelis_Menten.html.

• •

Folder tools/ contains some auxiliary functions called in the script (a simple Newton iteration ScSolve.m, an explicit initial-value-problem solver using the Dormand-Prince scheme and fixed step size ScIVP.m, and a function for approximating the Jacobian with finite differences ScJacobian.m.