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

3 MB·

by  Jan Sieber

Convergence of equation-free methods in the case of finite time scale separation
with application to stochastic systems
Jan Sieber∗, Christian Marschler†, and Jens Starke‡
Abstract. A common approach to studying high-dimensional systems with emergent low-dimensional behav-
ior is based on lift-evolve-restrict maps (called equation-free methods): first, a user-defined lifting
operatormapsasetoflow-dimensionalcoordinatesintothehigh-dimensionalphasespace,thenthe
high-dimensional(microscopic)evolutionisappliedforsometime,andfinallyauser-definedrestric-
tionoperatormapsdownintoalow-dimensionalspaceagain. Weproveconvergenceofequation-free
methods for finite time-scale separation with respect to a method parameter, the so-called healing
time.
7
Moreprecisely,ifthehigh-dimensionalsystemhasanattractinginvariantmanifoldwithsmaller
1
0 expansion and attraction rates in the tangential direction than in the transversal direction (normal
2 hyperbolicity), and restriction and lifting satisfy some generic transversality conditions, then the
implicit formulation of the lift-evolve-restrict procedure generates an approximate map that con-
n
verges to the flow on the invariant manifold for healing time going to infinity. In contrast to all
a
J previousresults, ourresultdoesnotrequirethetimescaleseparationtobelarge. Ademonstration
with Michaelis-Menten kinetics shows that the error estimates of our theorem are sharp.
1
Theabilitytoachieveconvergenceevenforfinitetimescaleseparationisespeciallyimportantfor
3
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.
S
Weinvestigatealow-dimensionalstochasticdifferentialequationwheretheratiobetweenthedecay
D
rates of fast and slow variables is 2.
.
h
Key words. implicit equation-free methods, slow-fast systems, stochastic differential equations, Michaelis-
t
a Menten kinetics, dimension reduction
m
[ AMS subject classifications. 65Pxx, 37Mxx, 34E13

1
v 1. Introduction. Dynamical systems with time scale separation have been studied for a
9
long time with different methods. These methods include homogenization methods [27, 31],
9
9 averaging and mean-field methods [35], the slaving principle or adiabatic elimination [15, 16]
8 andtechniquesfromslow-fastdynamicalsystems[11]. Theaimofallmethodsistoreducethe
0
complexity of a high-dimensional (also called microscopic) system to a relatively simple low-
.
1
dimensional(alsocalledmacroscopic)system. Thejustificationforthisreductionissimplestif
0
7 the underlying dynamical system possesses a low-dimensional attracting slow manifold. After
1 reduction, the long-term dynamics of the system can be analyzed on the slow manifold.
:
v Equation-freeframework. Thereisarangeofwell-establishednumericalmethodsthatavoid
i
X the explicit derivation of a macroscopic system by obtaining the required information from
r simulations. The unknown macroscopic dynamics is evaluated by using a wrapper around
a
∗College of Engineering, Mathematics and Physical Sciences, University of Exeter, North Park Road, Exeter
(Devon) EX4 4QF ([e-mail]).
†DepartmentofAppliedMathematicsandComputerScience,TechnicalUniversityofDenmark,Matematiktorvet
303B, DK-2800 Kgs. Lyngby, Denmark.
‡Institute of Mathematics, University of Rostock, Ulmenstraße 69, 18057 Rostock, Germany (jens.starke@
uni-rostock.de).
1
existingmicroscopicsimulatorstoachieveaclosureondemand. Amongthosemethodsarethe
recursive projection method [32], the heterogeneous multiscale method [10] and the equation-
free framework [18, 12, 20]. Here, we focus on the equation-free framework for microscopic
systems with time scale separation, as proposed originally by IG Kevrekidis. The assumption
behind equation-free computations is the existence of a slow low-dimensional description (in
Rd) for some macroscopic quantities of the high-dimensional microscopic system (which is
defined in RD). The framework also relies on the availability of a microscopic time stepper
(a map M(δ;·) : RD (cid:55)→ RD) and two user-defined operators, the lifting L : Rd (cid:55)→ RD and
the restriction R : RD (cid:55)→ Rd, which are maps between the original high-dimensional (RD)
microscopic level and the low-dimensional (Rd) macroscopic level.
The goal is to compose a macroscopic time stepper Φ (δ;·) : Rd (cid:55)→ Rd, which is then
∗
amenable to higher-level tasks. First use cases were macroscopic bifurcation analysis for
microscopic simulations in chemical engineering (see [19] for a review). Recently similar
analysis was performed on stochastic network models of neurons [3] or disease spread [14], or
on agent-based models in ecology [36] and social sciences (for example, for consumer lock-in
[1], for pedestrian flow [25, 24], or for trading [33]). Another example for a high-level task
accessible via equation-free methods is control design [34, 33].
The central building block of equation-free methodology is the “lift-evolve-restrict” map
R◦M(δ;·)◦L: for a given value x ∈ Rd of macroscopic quantities, one first applies the lifting
L to x getting a microscopic state u, then one runs the simulation for time δ starting from u
(applying the microscopic evolution M(δ;u)), and finally one applies the restriction R to the
result M(δ;u).
Figure 1.1. Geometry of lift-evolve-restrict map near slow manifold: macroscopic value x gets lifted to
L(x), then evolved to M(δ;L(x)), then restricted back to R(M(δ;L(x)) ∈ Rd. The aim is to approximate
the slow flow x (cid:55)→ (g ◦L)−1 ◦M(δ;g(L(x))) using this map R◦M(δ;·)◦L, and assuming invertibility of
g◦L:Rd (cid:55)→C.
Recent modifications and improvements to equation-free methods in multi-particle or
agent-based simulations are variance reduction [1, 8], restriction of computations to patches
in space [29, 30, 22] (for which a-priori error estimates can be proven [29, 30]), and automatic,
data-driven selection of the slow variables using diffusion maps [7, 25].
2
Current state of analysis. Analysis of the equation-free framework (based on lift-evolve-
restrict) is still incomplete. 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 C, which is rarely encountered in the
practical applications listed above (see [5] for an exception). Even for this idealized case one
faces the geometric difficulty, illustrated in Fig. 1.1, that in general the lift-evolve-restrict
map will not be compatible with the stable fibers of the slow manifold C. More precisely, after
lifting x ∈ Rd to L(x) ∈ RD, the slow flow is effectively applied to a different point, g◦L(x),
which is the projection of L(x) onto the slow manifold C along the stable fibers illustrated in
Fig. 1.1 (sometimes called isochrones; see Equation (2.3) in §2 for the precise definition of g).
Thus, in the limit of infinite time-scale separation (that is, the derivative of M with respect
to time, ∂ M(t;,x), goes to 0 for x ∈ C) the dynamics of the lift-evolve-restrict map
1
P(t;·) = R◦M(δ;·)◦L
isasmallperturbationofthemapR◦g◦L,whichisindependentoftheslowflow(forexample,
in the geometry shown in Fig. 1.1 with d = 1 and D = 2 the map R◦g◦L has a stable fixed
point). Hence, P(t;·) cannot be a good approximation of the slow flow along the manifold C,
unless R and L are chosen such that R◦g◦L is the identity. Using the coordinate x in the
domain of the lifting L and the map g ◦L : domL (cid:55)→ C, onto the manifold C the slow flow
Φ (δ;·) has the form
∗
y = Φ (δ;x) = (g◦L)−1◦M(δ;·)◦g◦L(x), or, implicitly defining y = Φ (δ;x),
∗ ∗ ∗ ∗
(1.1)
R◦g◦L(y ) = R◦M(δ;·)◦g◦L(x)
∗
(usingthenotation(·)−1 fortheinversemap). Thisdefinitionisimpracticalsincethenonlinear
projection g and the slow manifold C are both unknown in general. There are two ways to
overcome this problem.
Constrained runs. One approach to ensure that R ◦ g ◦ L is close to the identity is to
make sure that the lifting L maps onto the manifold C with sufficient accuracy for all x in its
domain. Usually,thisrequiresanadditionalschemeinvolvingtheiterativeapplicationofL;see
[12, 38, 39]. The a-priori error estimates prove that the lift-evolve-restrict scheme with these
additional iterations has an error of order (d /d )m, where d is the attraction/repulsion
tan tr tan
time scale tangential to the slow manifold C and d is the transversal attraction rate. The
tr
ratio d /d measures the time scale separation. It is assumed to be small when applying
tan tr
contrained runs (and called ε), and convergence is proven in [12, 38, 39] in the limit ε → 0
(which 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 , exploiting that M attracts along the fibers [20, 4]. Marschler et al [23]
skip
show that the healing time t can be justified by introducing an additional shift M(t ;·)
skip skip
and its inverse into (1.1) (note that M(t ;·) is invertible on the slow manifold C):
skip
y = Φ (δ;x) = (g◦L)−1◦M(t ;·)−1◦M(δ+t ;·)◦g◦L(x). (1.2)
∗ ∗ skip skip
Removing the inverses in (1.2) leads to an implicit equation for y = Φ (δ;x) with the healing
∗ ∗
time t as an additional parameter (which has no effect):
skip
R◦M(t ;·)◦g◦L(y ) = R◦M(δ+t ;·)◦g◦L(x) (1.3)
skip ∗ skip
3
However, M(t ;·)◦g−M(t ;·) is small (of order exp(−d t )) such that we may replace
skip skip tr skip
M(t ;·) ◦ g by M(t ;·) in (1.3). This results in a computable approximation y =
skip skip tskip
Φ (δ;x) of y , given implicitly by the equation
t ∗
skip
R(M(t ;L(y ))) = R(M(δ+t ;L(x))). (1.4)
skip tskip skip
Thisapproachwasanalyzedandillustratedinatrafficmodelin[23]andwillalsobestudiedin
this paper. Vandekerckhove et al [37] investigated a similar approach, but applied the healing
timebackwardintimebyfixingtheimageoftherestriction: theysolvex = R(M(t ;L(x )))
skip b
forx firstandthensety = R(M(δ+t ;L(x ))). Thisgivesan(approximate)representation
b skip b
ΦR of the slow flow in the coordinates on the image of the restriction R: ΦR(δ;x) = R ◦
∗ ∗
M(δ;·)◦[R| ]−1(x). Marschler et al [23] proved that the approximation y is exponentially
C tskip
accurate if d /d → 0: (cid:107)y −y (cid:107) ∼ exp(−Kd /d ) (for some constant K depending on
tan tr t ∗ tr tan
skip
t ). Theerrorestimatesin[23]requirethatt d /d and(t +δ)d /d staybounded
skip skip tan tr skip tan tr
from above such that the convergence result is about the limit of infinite time scale separation
d /d → 0, similar to the results for constrained runs schemes [12, 38, 39]. The analysis left
tan tr
open if the error goes to zero for t → ∞ but finite timescale separation: d /d ∈ (0,1).
skip tan tr
Section 2 will introduce the general a-priori error estimate that (cid:107)y −y (cid:107) ∼ exp((d −
t ∗ tan
skip
d )t ) for t → ∞ and fixed d < d under some genericity conditions on R and L. It
tr skip skip tan tr
will also give a convergence result for the derivatives of y with respect to its argument x:
t
skip
(cid:107)∂jy −∂jy (cid:107) ∼ exp(((2j +1)d −d )t ) if (2j +1)d < d .
tskip ∗ tan tr skip tan 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
assumedtohaveafastsubsystem,wherethefasttimescaleconvergesonlyinastatisticalsense
to a stationary measure conditioned on the slow variables. In these cases the microscopic time
stepper M(δ;·) operates on measures (or densities). It may be approximated by Monte Carlo
simulations on ensembles of initial conditions. Barkley et al [4] investigated the behaviour
of the lift-evolve-restrict map P(δ;·) = R◦M(δ;·)◦L where the slow variables were leading
moments (thus, P was called moment map in [4]) on prototype examples from the class of
stochastic problems. The simplest example from [4] is a scalar stochastic differential equation
(SDE), for which the evolution of the probability distribution is governed by the (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 [4] found that the dynamics of
the map P was qualitatively different from the dynamics of the underlying FPE. For example,
P was nonlinear and metastable states were stabilized for certain choices of δ.
Section 4 will demonstrate for two different lifting operators L that the approximation
y , defined by (1.4), behaves exactly as predicted by the convergence theorem presented in
t
skip
Section 2. In particular, it preserves the metastability features and the linearity of the flow
generated by the FPE.
Outlineofresults. Section2statesthepreciseassumptions(timescaleseparationfordecay
rates tangential and transversal to the invariant manifold C (d < d ) and transversality of
tan tr
R and L) for exponential convergence:
∂jy −∂jy ∼ exp(((2j +1)d −d )t ) for t → ∞
tskip ∗ tan tr skip skip
4
(using the convention that ∂0y = y and assuming that the derivatives up to order j+1 exist).
Section 3 demonstrates that the convergence rate is indeed different for the derivatives for a
singularly perturbed ODE modelling the Michaelis-Menten kinetics (which was also used by
[12, 38, 39] for illustration). Section 4 demonstrates convergence of the first two moments
of y to the first two moments of y for the distributions governed by the Fokker-Planck
t ∗
skip
equation of the scalar SDE with a double-well potential drift term studied by Barkley et al
[4]. We demonstrate global convergence for a linear lifting L (a linear combination of initial
lin
Gaussiandistributions). WealsodemonstratelocalconvergenceforthenonlinearliftingL
gauss
used in [4] (the prefactor in the convergence rate depends on the initial point x).
Section 5 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, for example, increasing the number of agents or particles
does not increase the spectral gap (and, thus, the time scale separation). Thus, the results
from Section 2 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 [12, 38, 39] and implicit lifting [25].
2. Convergence in the case of finite time-scale separation. We consider a smooth dy-
namical system
u˙ = f(u), u ∈ RD, (2.1)
where D is large. We assume that the flow M generated by (2.1),
M : R×RD → RD, (t;u) (cid:55)→ M(t;u)
has a d-dimensional compact relatively invariant manifold C (possibly with boundary). That
is, trajectories M(t;u) starting in u ∈ C either stay in C for all times t ∈ R, or they stay in C
until they cross the boundary ∂C of C. We assume that C is at least k times differentiable.
max
For a point u ∈ C, let us denote by N(u) the d-dimensional tangent space to C. The following
assumption states that attraction transversal to the manifold C is faster than attraction or
expansion tangential to C.
Assumption 1 (Hyperbolicity — Separation of time scales and transversal stability).Thereex-
ists an open neighborhood U of the manifold C, a (possibly) nonlinear projection g : U (cid:55)→ C
(the so-called stable fiber projection), a pair of constants (decay rates) 0 < d < d , and a
tan tr
bound C such that the following two conditions hold.
1. (tangentialexpansion/attractionrate) For all points u ∈ C on the manifold, all tangent
vectors v ,...,v ∈ N(u) and all t ∈ R with M(t;u) ∈ C
1 kmax
(cid:107)∂jM(t;u)[v ,...,v ](cid:107) ≤ Cexp(d |t|)(cid:107)v (cid:107)·...·(cid:107)v (cid:107) for all j ∈ {1,...,k }.
2 1 kj tan 1 j max
(2.2)
5
2. (Stability along transversal fiber projections) For all u ∈ U and all t > 0 with
M(t;g(u)) ∈ C
(cid:107)∂jM(t;u)−∂jM(t;g(u))(cid:107) ≤ Cexp(−td )(cid:107)u−g(u)(cid:107) for all j ∈ {0,...,k }.
2 2 tr max
(2.3)
In (2.3) we use the convention that ∂jM is the jth-order partial derivative of M with respect
k
to its kth argument, and that ∂0M (j = 0) is the flow M itself. The norm on the left side
2
of (2.3) is the usual operator norm for the multi-linear operators ∂jM(t,·). The constants
2
C, d and d are assumed to be independent of the point u and the time t. Assumption
tr tan
(2.2) is also made for negative times t such that it is also an assumption about the inverse
of M: M(−t,·) = M−1(t,·). The constant d is the decay rate toward the manifold C, the
tr
constant d is the rate of attraction and expansion along the flow restricted to C. The main
tan
requirement of Assumption 1 is that d > d .
tr tan
Transversality of restriction and lifting. Second, we assume basic compatibility between
R : U ⊂ RD → Rd the restriction operator,
(2.4)
L : domL ⊂ Rd → RD the lifting operator,
and the invariant manifold C: the lifting L should map into the neighborhood U of C in which
the projection g is defined, and the restriction R should be defined on the projection of the
image of L along the stable fibers:
L(domL) ⊂ U, g(L(domL)) ⊂ domR∩C,
In addition to these compatibility conditions, we impose the following two tansversality con-
ditions on lifting L and restriction R.
Assumption 2 (Transversality of R and L).
1. the projection g is a diffeomorphism between rgL = L(domL) and C. In particular,
for all x ∈ domL ⊂ Rd
∂
rank [g(L(x))] = rank[∂g(L(x))◦∂L(x)] = d.
∂x
2. The map R, restricted to C, is a diffeomorphism between C and Rd. In particular, for
all u ∈ C (N(u) is the tangent space to C in u)
dim∂R(u)N(u) = d.
Coordinates on the slow manifold C. The maps R and L create two natural ways to define
local coordinate representations on the invariant manifold C, one by a map from domL to
C, one by a map from C to rgR. Note that domL is not necessarily equal to rgR. For our
presentation we choose the representation in domL coordinates:
g◦L : domL ⊂ Rd (cid:55)→ C ⊂ RD, x (cid:55)→ g(L(x)).
6
The inverse of g ◦L is defined implicitly. Assume that u = g(L(x )) for some x ∈ domL.
0 0 0
Then for u ∈ C near u the pre-image x = (g◦L)−1(u) is found by solving R(u) = R(g(L(x)))
0
for x ≈ x , which has a unique solution by Assumption 2.
0
We can represent the flow M on C as a flow in domL, denoting it by Φ :
∗
Φ : R×domL (cid:55)→ domL, Φ (δ;x) = [(g◦L)−1◦M(δ;·)◦g◦L](x) := y (2.5)
∗ ∗
(for δ ∈ R and x ∈ domL), where y is given implicitly as solution of a d-dimensional system
of nonlinear equations
R(g(L(y))) = R(M(δ;g(L(x)))). (2.6)
Assumption 2 on transversality implies that Φ is well defined for small δ (since y = x is a
∗
regular solution of (2.6) at δ = 0). For larger δ, one can break down the solution into smaller
steps by increasing δ gradually from 0 and tracking the curve y(δ) of solutions of (2.6), which
is well parametrized by δ in every point by Assumption 2. If domL is simply connected then
this continuation approach makes the implicit solution y used in the definition of Φ unique.
∗
Let us define the map
P : R×domL (cid:55)→ Rd, P (t;x) = R(M(t;g(L(x)))). (2.7)
∗ ∗
ThismapP iswelldefinedandinvertibleforallt ∈ Randx ∈ domLsatisfyingM(t;g(L(x))) ∈
∗
C. The implicit definition (2.5) of the flow Φ on C has the following form when expressed
∗
with the help of this map P on domL:
∗
y = Φ (δ;x) if P (0;y) = P (δ;x). (2.8)
∗ ∗ ∗
Since the flow M(δ;·) is a diffeomorphism on C, we can replace the times 0 and δ in the above
implicit definition with t and t +δ for an arbitrary so-called healing time t ∈ R (as
skip skip skip
long as M(t ;g(L(x))) ∈ C). So, equivalent to (2.8):
skip
y = Φ (δ;x) if P (t ;y) = P (t +δ;x). (2.9)
∗ ∗ skip ∗ skip
ConvergenceTheoremforimplicitequation-freecomputationswithfinitetime-scaleseparation.
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 (2.9):
∗
P : [0,∞)×domL (cid:55)→ rgR P(t,x) = R(M(t;L(x))) = [R◦M(t;·)◦L](x) (2.10)
such that we may define correspondingly the approximate flow map
Φ : R×domL (cid:55)→ domL where y = Φ (δ;x) if P(t ;y) = P(t +δ;x). (2.11)
tskip tskip skip skip
Our general convergence theorem, Theorem 2.1, states that Φ is well defined for large t
tskip skip
(that is, the implicit equation defining it has a locally unique solution), and that ∂jΦ is
t
skip
an approximation of ∂jΦ of order exp(((2j+1)d −d )t ) (including j = 0 for the map
∗ tan tr skip
Φ ).
t
skip
7
Theorem 2.1 (Convergence of approximate flow map at finite time-scale separation).
Let us assume that the microscopic flow M satisfies Assumption 1 on time-scale separation,
and that the maps R and L satisfy Assumption 2 on transversality.
Let δ > 0 and x ∈ domL be arbitrary. Let us also assume that x ∈ domL maps into
max
a point under g ◦ L that keeps a positive distance from the boundary ∂C of C for all times
t ≥ −δ under M. That is.
max
dist(M(t;g(L(x))),∂C) ≥ c for all t ≥ −δ and some given c > 0. (2.12)
∂ max ∂
Then there exists a t ≥ δ such that y = Φ (δ;x) is well defined by (2.11) for all
0 max t
skip
δ ∈ [−δ ,δ ] and t > t . The estimate
max max skip 0
(cid:107)∂jΦ (δ;x)−∂jΦ (δ;x)(cid:107) ≤ Cexp(((2j +1)d −d )t ) (2.13)
2 tskip 2 ∗ tan tr skip
holds for all orders j ∈ {0,...,k −1} satisfying (2j+1)d < d . The constant C depends
max tan tr
on δ and x, but not on j or t .
max skip
Assumption (2.12) in Theorem 2.1 is made to permit arbitrarily large t while still
skip
having Assumption 1 and Assumption 2 uniformly satisfied. If one considers x ∈ domL for
which the trajectory t (cid:55)→ M(t;g(L(x))) leaves C (by crossing the boundary ∂C) then one
has to put restrictions on δ and t to avoid crossing ∂C. The theorem permits negative
skip
integration times δ shorter than t ≤ t and positive integration times larger than t as
0 skip skip
long as the factor exp(d |δ|) is of order 1. The statement of Theorem 2.1 does not require
tan
that the time-scale separation d /d goes to zero for convergence of the approximate map.
tan tr
It only requires that d < d , where d is the attraction rate along fibers (see (2.3)) and
tan tr tr
d is the attraction and expansion rate tangential to the invariant manifold C.
tan
Outline of proof of Theorem 2.1. (Existence and error of Φ ) The proof of Theorem 2.1
t
skip
uses the following implicit fixed-point problem for y (denoting the true solution as y =
∗
Φ (δ;x)) :
∗
P (t ;y) = P (t ;y )+[P (t ;y)−P(t ;y)]+[P(t +δ;x)−P (t +δ;x)].
∗ skip ∗ skip ∗ ∗ skip skip skip ∗ skip
(2.14)
Thenormsofbothtermsinbrackets,P (t ;y)−P(t ;y)andP(t +δ;x)−P (t +δ;x),
∗ skip skip skip ∗ skip
are of order exp(−d t ) by Assumption 1, equation (2.3) (transversal stability with rate
tr skip
d of C). For the same reason, the Lipschitz constant of P (t ;y)−P(t ;y) is of order
tr ∗ skip skip
exp(−d t ),too. ByAssumption1,equation(2.2)(tangentialdecayrateinsidethemanifold
tr skip
is less than d ), and Assumption 2 on transversality of L and R the inverse of P (t ;·) has
tan ∗ skip
a local Lipschitz constant of order exp(d t ) near y . These two facts enable us to apply
tan skip ∗
the Banach Contraction Mapping Principle to (2.14) to obtain a unique solution y ≈ y for
∗
large t . More precisly, y−y is of order exp((d −d )t ).
skip ∗ tan tr skip
(Inductive proof of error estimate for derivatives) We differentiate (2.14) with respect to
x in its fixed point y(t ;x) up to j times and then re-arrange the resulting equation for the
skip
jth-order derivatives of y and y into the form
∗
∂P (y)(cid:2)∂jy−∂jy (cid:3) = [∂P (y)−∂P (y )]∂jy +r. (2.15)
∗ ∗ ∗ ∗ ∗ ∗
8
In (2.15) we abbreviated ∂ P (t ;·) = ∂P (·), ∂jy = ∂j(t ;x) and dropped the argument
2 ∗ skip ∗ 2 skip
x from y and the arguments t and x from y. The remainder r is less than Cexp(((2j −
∗ skip
1)d −d )t )forsomeconstantC byinductiveassumption. Theimplicitexpression(2.15)
tan tr skip
for ∂jy−∂jy shows why errors in derivatives of the solution can grow for increasing t and
∗ skip
insufficient time scale separation: the norms of ∂P (y)−∂P (y ) and of [∂ P (t ;y)]−1 are
∗ ∗ ∗ 2 ∗ skip
of order exp(d t ) due to (2.2). The details of the proof are given in Appendix A.
tan skip
3. Example: Michaelis-Menten kinetics. Toillustratetheconsequencesoferrorestimate
(2.13), we look at a model for Michaelis-Menten kinetics with explicit time scale separation
as studied in [26, 12, 38, 39]. The system is given in RD with D = 2 as
x˙ = ε[−x+(x+κ−λ)y],
(3.1)
y˙ = x−(x+κ)y,
where x ∈ R is the slow variable, y ∈ R is the fast variable, and ε measures the time scale
separation. Theparametersκ = 1,λ = 0.5andε = 0.01arekeptfixedthroughoutthesection.
Inthesingularcaseε = 0,system(3.1)hasacriticalmanifoldC . Near-byisatransversally
0
stable invariant manifold, which can be represented as a graph y = h (x) such that d = 1.
ε
The graph h can be expanded in ε for small ε > 0:
ε
x κλx κλx(2κλ−3λx−κx−κ2)
y = + ε+ ε2+O(ε3). (3.2)
x+κ (x+κ)4 (x+κ)7
The dynamics of (3.1) in phase space for different initial conditions is shown in Fig. 3.1(a),
where the fast dynamics in y is observed. After a short transient the trajectories approach
the slow manifold given approximately by (3.2) (shown as the black line in Fig. 3.1(a)).
Wespecifytherestrictionandliftingoperators, RandL, fortheMichaelis-Mentensystem
as
(cid:18) (cid:19) (cid:18) (cid:19)
x x
R : R2 (cid:55)→ R, R = x, and L : R (cid:55)→ R2, L(x) = . (3.3)
y 0.5
The approximate time-δ map Φ (δ;·) on the slow manifold is determined by the root
t
skip
z∗ of
F(z) = R(M(t ;L(z)))−R(M(t +δ;L(x))) (3.4)
skip skip
and setting Φ (δ;x) := z∗ (cf. (2.11)). Note, that F depends on t and δ, which are not
tskip skip
included in the list of arguments to simplify notation. We solve (3.4) using a Newton method
zn+1 = zn−J−1F(zn), n = 0,1,2,... (3.5)
where J is the Jacobian of F with J = ∂F /∂z | . For this example, Jacobians are approx-
ij i j zn
imated numerically using a central finite difference scheme
(cid:12)
∂F (cid:12) F(zn+∆z)−F(zn−∆z)
i(cid:12) ≈ , (3.6)
∂z (cid:12) 2∆z
j(cid:12)
zn
9
0.5 10-2
0.4 10-4
10-6
0.3
E
y0.2 rorre 10-8
10-10 )
d)=dx
0.1
d2)=dx2
10-12 exp(!1:3tskip)
0 0 10 20 30
0 0.05 0.1 0.15 0.2 0.25
t
x skip
(a) (b)
0.5 100
10-2
10-4
E
w 0 ro
rre 10-6
)
10-8 d)=dx
d2)=dx2
−0.5 10-10 exp(!tskip)
0 10 20 30
−2 −1 0 1 t
v skip
(c) (d)
Figure 3.1. Michaelis-Menten dynamics (3.1) for λ=1,κ=0.5,ε=0.01. (a) Phase space representation
ofthedynamicsforvariousinitialconditions(crosses). Thedynamicsisquicklyapproachingtheslowmanifold
(black line, see (3.2)) and evolves on a much slower time scale. Each simulation is performed for t = 100
time steps with a fixed time difference of 0.1. The final states are denoted by a red point. (b) Error (3.9) as a
function of t . (c,d) Same analysis for the rotated system (3.10). The initial large error in the flow of size
skip
1 can be significantly reduced to an order of 10−8 for a larger healing time.
with ∆z = 10−4. The Newton iteration is considered to be converged if both (cid:107)F(zn)(cid:107) and
(cid:107)zn+1−zn(cid:107) are smaller than a given tolerance tol = 10−10.
Higher-order derivatives of the flow are approximated by finite differences
Φ (δ;x+∆x)−Φ (δ;x−∆x)
∂1Φ (δ;x) ≈ tskip tskip , (3.7)
2 tskip 2∆x
Φ (δ;x+2∆x)−2Φ (δ;x)+Φ (δ;x−2∆x)
∂2Φ (δ;x) ≈ tskip tskip tskip , (3.8)
2 tskip 4∆x2
with ∆x = 10−4 .
Sincethe stablefiberprojectionsg containedin theexactflow Φ (see(2.6)) areingeneral
∗
unknown, the error (2.13) cannot be measured directly. Nevertheless, we might analyze the
systematic change in Φ (δ;x) with t . The error is then measured as the distance to the
tskip skip
10