Computing stable hierarchies of fiber bundles

Stable fiber bundles are important structures for understanding nonautonomous dynamics. These sets have a hierarchical structure ranging from stable to strong stable fibers. First, we compute corresponding structures for linear systems and prove an error estimate. The spectral concept of choice is the Sacker-Sell spectrum that is based on exponential dichotomies. Secondly, we tackle the nonlinear case and propose an algorithm for the numerical approximation of stable hierarchies in nonautonomous difference equations. This method generalizes the contour algorithm for computing stable fibers from [38, 39]. It is based on Hadamard’s graph transform and approximates fibers of the hierarchy by zero-contours of specific operators. We calculate fiber bundles and illustrate errors involved for several examples, including a nonautonomous Lorenz model.


Introduction
Invariant manifolds in dynamical systems partition the state space into different areas of attraction.They are of great dynamical relevance in various applications like mathematical biology or chaos theory and we refer to [53, Section 1.1] for further applications.Existence results for invariant manifolds of hyperbolic fixed points in autonomous dynamical systems are given, e.g., in [28,34,45,51] and corresponding results for nonautonomous fiber bundles of bounded trajectories have been developed in [2,4].
Finding invariant manifolds analytically is hardly possible for many relevant models and numerical approximations are the only feasible alternatives.The literature in this area is quite vast.Proposed techniques are based on numerical continuation, boundary value problems, Taylor expansions, the parameterization method, fixed point iterations and set oriented methods, see [8,10,17,20,22,23,32,42,49,52], where this list is by no means complete.Some of these methods allow nonautonomous generalizations, e.g.[9,11,27,48,49], but the literature in this area is quite sparse.Finally, we mention techniques for rigorously computing normally hyperbolic manifolds and for proving their existence, see [12,13,24].
The problem gets more involved when asking for hierarchies of invariant manifolds, see [51,Chapter 5].The stable hierarchy ranges from the stable to the strong stable manifold and provides information on regions within the stable manifold, having different rates of attraction.
We even go one step further and aim for computing hierarchies in nonautonomous discrete time systems The family of maps (F n ) n∈ may arise via discrete time modeling or it may result from a discretization of a nonautonomous continuous time system.The reference object for which we define stable hierarchies are bounded trajectories ξ = (ξ n ) n∈ of (1).Fixed points that are frequently chosen as reference objects in autonomous systems, typically do not exist in the nonautonomous context.
Fibers from the stable hierarchy can locally be expressed as graphs over the corresponding stable subspace, see [4,Theorem 4.11] and [46,Theorem 4.6.4].Before we propose an approach for their approximation in Section 3, we specify the precise meaning of a hierarchy of stable subspaces.For this purpose, a linear theory for the so called variational equation is introduced in detail in Section 2. If the underlying system is autonomous and ξ is a fixed point, then (2) turns into a linear, autonomous difference equation Assuming that A is hyperbolic, i.e.A has no eigenvalue on the unit circle, it follows that eigenvalues with corresponding generalized eigenvectors provide complete information on the underlying dynamics.The analysis changes dramatically, when considering the nonautonomous equation (2) or, more generally This setup requires adequate spectral concepts, since eigenvalues are in general meaningless as shown by Vinograd, cf.[14,Section 2.6].The Sacker-Sell spectrum, also called dichotomy spectrum provides the desired stability information.
Of particular interest are corresponding spectral bundles that are the nonautonomous analog of eigenspaces.We review in Section 2 an algorithm from [25, Section 3] for computing spectral bundles and recall a precise estimate of the errors involved from [25,Section 3.1].Its proof that was postponed in [25] to a forthcoming publication is presented here.An example illustrates that this estimate is rather sharp.After establishing the linear theory, we tackle nonlinear approximations of stable fibers in Section 3. The main tool is the contour algorithm that was proposed in [38,39].Note that F n is not assumed to be invertible.The algorithm is based on Hadamard's graph transform [28], but avoids an infinite fixed point iteration.In the end, the problem reduces to the computation of the zero-contour of some function f : Ê r → Ê, where r is the dimension of the stable fiber.For this task, efficient software applies, such as the Matlab commands contour and isosurface in cases r = 2 and r = 3, respectively.This ansatz restricts the space dimension d to d ∈ {2, 3} and requires the unstable fibers to be onedimensional.Overcoming several of these restrictions constitutes a novelty of this paper.Furthermore, an extension of the contour algorithm is introduced that enables approximating the whole stable hierarchy.Error estimates justify the numerical procedure.
We apply the resulting algorithm to various examples of dimension three and four.In particular, error estimates are verified numerically for models with explicitly given stable hierarchies.Numerical strategies for the computation of hierarchies are illustrated in case of non-global graph representations.Finally, we determine the stable hierarchy for a nonautonomous variant of the famous Lorenz system [43].

Computation of spectral bundles for linear systems
Spectral bundles in nonautonomous linear systems (4) are the analog of eigenspaces in autonomous systems.We review an algorithm from [25] for their numerical approximation and analyze occurring errors.Finally, we illustrate that the derived rates of convergence are rather sharp.
We start with the introduction of our basic notion of hyperbolicity -the so called exponential dichotomy.

Exponential dichotomy
The solution operator Φ of the difference equation ( 4) is given by Let P s,u be two families of complementary, invariant projectors, i.e. ) is used for its inverse.

Definition 1
The difference equation (4) has an exponential dichotomy (ED for short) with data (K, α s,u , P s,u ) on , if there exist constants K, α s , α u > 0 and two families of projectors P s,u that satisfy (5), such that (ii) For n, m ∈ , n ≥ m the following estimates hold:

Sacker-Sell spectrum
The Sacker-Sell spectrum, cf.[50], also known as dichotomy spectrum, is based on considering the scaled equation We define the Sacker-Sell spectum Σ ED := {γ > 0 : (6) has no ED on } and the resolvent set R ED := Ê \ Σ ED .Alternatively, this spectral notion may be derived from operator theory, cf.[47].For n ∈ let (T u) n := A n−1 u n−1 and denote by Σ(T ) the spectrum of this linear operator in suitable spaces of bounded sequences.
Assuming that all matrices A n ∈ Ê d,d are uniformly bounded, it follows from the Spectral Theorem in [5] that Σ ED consists of ℓ ≤ d intervals of the form with two possible cases for I ℓ .In case 1, see Figure 1.
Illustration of spectral and resolvent intervals.
The upper case in Figure 1 may occur, if A n is noninvertible for some n or if Finally, we note that the Sacker-Sell spectrum is closely related to the set of Bohl exponents, see [40,Appendix B.2].

Spectral bundles
For i ∈ {1, . . ., ℓ + 1} and γ ∈ R i we denote by P s,u ,i the dichotomy projectors of the scaled equation (6).Note that dichotomy projectors on are uniquely determined.In particular, dichotomy projectors do not depend in the i-th resolvent interval R i on the specific value of γ ∈ R i .The regularity condition (i) of Definition 1 guarantees that for all γ ∈ R i , the solution operator Φ(•, n) is invertible in R(P u n,i ).Further note that R(P s n,1 ) = Ê d and are well defined, since the dimension of R(P s n,i ) does not depend on n.The value of d s i is strictly decreasing with respect to i. Dichotomy projectors have a hierarchical structure for n ∈ and it follows for all n ∈ and i = 1, . . ., ℓ that We introduce the spectral bundles For i = 1, . . ., ℓ, we note that = {x ∈ Ê d : P u n,i+1 P s n,i x = x} = R(P u n,i+1 P s n,i ).

Spectral bundles and Lyapunov transformations
It is well known that the Sacker-Sell spectrum is invariant w.r.t. a Lyapunov transformation of the equation, i.e. the Sacker-Sell spectrum of (4) coincides with the spectrum of where S n ∈ Ê d,d , n ∈ are invertible, and both S n and S −1 n are uniformly bounded.Corresponding dichotomy projectors are given by S −1 n P s,u n,i S n and consequently, spectral bundles Wi n of (13) satisfy Wi n = S −1 n W i n for i = 1, . . ., ℓ + 1 and n ∈ .
Two systems (4) and ( 13) that are conjugate via a Lyapunov transformation are called kinematically similar.

An algorithm for computing spectral bundles
In the following, we denote by J = [n − , n + ] ∩ and J = [n − , n + − 1] ∩ finite discrete intervals.Our approach for the approximation of spectral bundles W i n is based on techniques for the calculation of dichotomy projectors from [35, Section 2] and [36, Section 2] and on the observation that W i n = R(P s n,i P u n,i+1 ) for i = 1, . . ., ℓ, see (12).The algorithm proposed in [25,Section 3] allows the computation of w = P s N,i P u N,i+1 r for a given vector r ∈ Ê d along the following steps: Approximation of x = P u N,i+1 r: Compute the least squares solution of where δ denotes the Kronecker symbol.Define Approximation of w = P s N,i x: Compute the least squares solution of and define w := u N .

Error estimates
By calculating the unique bounded solutions of ( 15) and ( 17) on , we obtain the vector w ∈ W i N .In practice, we are forced to solve ( 15), ( 17) on finite intervals and numerical errors cannot be avoided.We prove that errors decays exponentially fast by increasing the length of the finite interval J. Furthermore, our estimates indicate optimal choices of γ i and γ i+1 inside corresponding resolvent intervals.
Our main interest is to accurately approximate the subspace W i N rather than the particular vector w ∈ W i N .Thus, approximation errors within W i N are irrelevant and it suffices to discuss the relevant error .
Denote by (K, α s,u j , P s,u ,j ) the corresponding dichotomy data of the scaled equation (6).Let J = [n − , n + ] with N ∈ J. Then the relevant error d N,i of the least squares approximation w = u N of (15), (17) satisfies the inequality with some constant C that does not depend on n ± .
Proof: Denote by Φ j the solution operator of ( 6) with γ = γ j , j ∈ {i, i + 1}.Let ū be the unique bounded solution on of (17).Here, x is defined in ( 16), using the least squares solution v J of (15) w.r.t. the finite interval J.It follows that where We derive estimates for the three terms, starting with d 1 N,i .The unique bounded solution of ( 17) on is given by ūN = P s N,i x, see [36,Theorem 2].Thus, we conclude using ( 9), ( 16) that with least squares solution v J from (15).The vector v N −1 can be decomposed as follows where η v and ξ v are uniformly bounded by some constant L > 0 due to [35,Theorem 4.1].Indeed, η v and ξ v decay exponentially fast to 0 as n ± → ±∞.Inserting ( 21) into ( 20), one gets For an estimate of the remaining terms, we observe from [35,Theorem 4.1] the similar decomposition where η u , ξ u ≤ L. We conclude that Finally, one obtains from ( 9), ( 10), (22) and Combining the three estimates, the proof of ( 19) with C = 3KL is complete.
It follows from that the stable dichotomy rate α s i+1 is maximal, if γ i+1 is chosen close to the right boundary of the resolvent interval R i+1 .On the other hand, a choice of γ i close to the left boundary of R i results in a maximal unstable rate α u i .Pushing these values towards the right respectively left boundary of the corresponding resolvent interval, we obtain in the limit the maximal rates see Figure 2.With these choices, Theorem 2 guarantees optimal rates of convergence of the algorithm from Section 2.5.Using (24), our error estimate (19) reads Note that in case i = 1 there is no unstable direction and the second term in (25) vanishes in accordance with the formal result

A linear example with explicitly given spectral bundles
We apply the algorithm from Section 2.5 to a model, for which spectral bundles are explicitly known.This model is a noninvertible version of [25,Section 5.1].An inspection of approximation errors particularly illustrates that our estimate ( 19) is rather sharp.
We define the diagonal matrix B := diag(σ 1 , . . ., σ d ) and arrive at a nonautonomous system (4) via the Lyapunov transformation Here, X is a sequence of random matrices with independent, uniformly distributed entries in [0, 1].For this model, ℓ = d − 1 and the Sacker-Sell spectrum Σ ED = ℓ i=1 {σ i } consists of isolated point only, see Figure 3.
If follows from (14) that spectral bundles at time N = 0 have the explicit representation , where e i denotes the i-th unit vector in Ê d .
We compute fibers W i 0 by applying the algorithm from Section 2.5.First, one needs an approximation of the Sacker-Sell spectrum in order to set up γ i and 0 Figure 3: Sacker-Sell spectrum and resolvent set of (26).
γ i+1 close to the boundaries of corresponding spectral intervals.In this particular example, we use the explicitly given spectrum.Note that various techniques for calculating spectral intervals have been developed.We refer to [18,19] for SVDand QR-methods that also apply to discrete time systems.An alternative ansatz that is based on solving boundary value problems was proposed in [36].
Let ℓ = 6 and choose γ i+1 = σ i − ε and γ i = σ i + ε with ε = 0.01, i = 1, . . ., 6. Next, the dependence of ( 19) on n + is verified numerically.For this task, we fix n − = −400, resulting in exponentially small errors from the left boundary, which can be neglected.Occurring errors d 0,i , i = 1, . . ., 6, cf. ( 18), of spectral bundles are computed for n + = 10, . . ., 100. Figure 4 illustrates the numerical data.From (19), an error of the form d 0,i ≈ Ce −α u i n + is expected.We compare the exponential slope of the data in Figure 4 with an estimate of α u i that is obtained from the length of the corresponding resolvent interval, see (24).Recall that α u 1 = ∞ and consequently, we expect no approximation error when computing d 0,1 .This is confirmed by the numerical experiment in Figure 4.     Finally, we note that in [25], the algorithm from Section 2.5 is applied to models from particle dynamics and to time-dependent fluid flows.These models fit into our linear framework, after applying a finite rank approximation of the Perron-Frobenius opearator that tracks the evolution of densities over short time intervals, see [26] for details.

Stable hierarchies of fiber bundles for nonlinear systems and their approximation
We now turn our attention to nonautonomous nonlinear systems.
Recall that the linear system (4) has unique dichotomy projectors, when considered on .As a consequence, spectral bundles W i n -as defined in Section 2.3 -are uniquely determined.
For the system (1), uniqueness of a nonlinear version of W i n cannot be expected.However, strong stable fibers are uniquely determined in nonlinear systems.Assuming 1 ∈ R ED and choosing k such that 1 ∈ Rk, strong stable fibers may be written for a linear system as W(j) := ℓ+1 i=j W i n , j = k, . . ., ℓ + 1.In a two-dimensional system with two disjoint spectral intervals in (0, 1), for example, the case j = 1 yields the stable fiber W(1) = Ê 2 , while W(2) describes the strong stable fiber.
Uniqueness of the strong stable fiber (in black) and of the stable fiber (gray plane) is illustrated in Figure 6.Non-unique objects, like weak stable fibers are shown in green (solid and dashed curves).
The nonlinear generalization of the spaces W(j), j = 1, . . ., ℓ + 1 leads to the so called hierarchy of invariant fiber bundles.Due to [4,Theorem 4.11] and [46,Theorem 4.6.4],these bundles exist and possess a graph representation (under reasonable assumptions), describing fibers as graphs over ranges of dichotomy projectors.For the latter objects, approximation results have been derived in Section 2. Based on these results, we tackle the calculation of hierarchies of fiber bundles in nonlinear systems by using a modified version of the contour algorithm from [38,39].
In particular, we demonstrate that the contour algorithm applies to systems of dimension higher than three, allowing in this way the approximation of twoand three-dimensional projections of high-dimensional fiber bundles.

Nonlinear setup
We consider the nonlinear system (1) and assume that the corresponding solution operator.
First we define a reference solution w.r.t. which we compute a hierarchy of fiber bundles.Since the system (1) typically does not have an n-independent fixed point, the only meaningful nonautonomous analog turns out to be a bounded trajectory ξ of (1): Using Taylor's expansion, we rewrite and suppose that: (A3) For each ε > 0 there exists a δ > 0 such that DQ n (x) ≤ ε for all x − ξ n ≤ δ and all n ∈ .
We require a set-valued Lipschitz condition, see [1, Definition 1.4.5]: (A4) On each compact set K ⊂ Ê d there exists a constant L F,K such that x − y B for all x, y ∈ K and all n ∈ .Here, B denotes the unit ball in Ê d .
We transform the bounded trajectory ξ to zero and rectify subspaces via a kinematic transformation.Let where B(W i n ) denotes a basis of W i n .The transformed system is given as If F satisfies (A1), (A3)-(A5) w.r.t. the bounded solution ξ , then G satisfies these assumptions w.r.t. the zero trajectory x = 0 .Note that (A4) holds for G with a slightly modified Lipschitz constant L G,K .By construction, G n satisfies for n ∈ Note that the upper index of the matrices A i indicates their stability properties from strongly stable (i = ℓ + 1) to strongly unstable (i = 1).The rectified system (27) is introduced for simplicity of presentation only.Our algorithm applies to the original system (1); the kinematic transformation is never computed numerically.Denote by d s,u i , i = 1, . . ., ℓ + 1 the dimension of the stable respectively unstable dichotomy projector in the i-th resolvent interval R i , see (7).We get for k = 1, . . ., ℓ + 1 and all n ∈ , We conclude that g ℓ+1,k n ,k and Φ k−1,1 be the solution operators of and recall that the Sacker-Sell spectrum of the corresponding variational equation is invariant w.r.t. the kinematic transformation (27).The construction of W i n yields that A i n is invertible for all i = 1, . . ., ℓ, but A ℓ+1 n might be singular.We deduce for all n ≥ m that with dichotomy constants −β s k < β u k .The exponential rates β s,u k belong to the unscaled variational equation ( 2) and they may be negative.By considering the scaled equation ( 6) with γ ∈ R k , we obtain an exponential dichotomy with rates α s,u k that have, due to (23), the form

Stable hierarchy of invariant fiber bundles
We define a hierarchy of locally invariant fibers under the assumptions (A1)-(A5).For simplicity, we start with the transformed system (27).By [4,Theorem 4.11] and [46,Theorem 4.6.4]there exists a hierarchy of invariant fiber bundles that are given for k = 1, . . ., ℓ + 1 and n ∈ as These fiber bundles are locally invariant in the following sense.Let y = x h k n (x) Furthermore, L k n is a nested sequence of locally invariant fiber bundles Assuming (A5), there exists an index 1 ≤ k ≤ ℓ + 1 such that 1 ∈ Rk.Then L k n denotes the stable n-fiber of ( 27) and defines the hierarchy of local stable fibers.
Returning to the original system (1), the hierarchy of locally invariant fibers reads for k = 1, . . ., ℓ + 1 and n ∈ where U denotes a sufficiently small neighborhood of 0. Let U n = U + ξ n for n ∈ .We introduce a global continuation of the local stable fibers (29) for k = k, . . ., ℓ + 1 ) Finally, the global stable hierarchy is given by the nested sequence In autonomous systems, this hierarchy ranges from the strong stable to the stable manifolds.We refer to [51, Chapter 5, Appendix III] in case the reference object is a fixed point of a diffeomorphism.For autonomous models, lower indices in (31) are skipped, since these sets are time independent.

Approximation of bounded trajectories
Before we start with the computation of fiber bundles, an accurate approximation on a interval [n − , n + ] of the reference object, i.e. the bounded trajectory ξ from (A2) is needed.For this task, we choose a buffer interval [m − , m + ], where m − < n − < n + < m + .Then, we solve the periodic boundary value problem using Newton's method.Note that due to its sparse structure, the linear system that occurs in each Newton step is efficiently solvable.It follows from (A5) that the variational equation ( 2) has an ED on .Exploiting this hyperbolicity condition yields that any fixed tolerance ∆ is reached on [n − , n + ], if the buffer interval [m − , m + ] is sufficiently large, see [37,Theorem 8].We neglect the influence of this small error on the output of the contour algorithm in the following sections.

Approximation results for stable hierarchies
In practical applications, the graph representation of L k n is typically not explicitly given.We approximate stable fiber via the following sets: where U is a sufficiently small neighborhood of ξ .These sets are nested for all n ∈ and all m, p ∈ AE. (34) Unless the graph representation L k n is global for all n, it is essential for the proof of Theorem 3 and for defining (30) that starting points from H k n converge forward in time towards the bounded trajectory ξ .This convergence holds true if and only if k ≥ k and consequently, the approximation result cannot be extended to unstable fibers.
The parameter m in (32) controls the number of steps inside the neighborhood U .Using the estimate (34), we conclude that all elements of S k n (m, p) lie arbitrarily close to H k n , if p is fixed and m is chosen to be large.The parameter p in (32) controls the number of steps outside the neighborhood U and thus determines the length of the computed fiber.Outside the local neighborhood U , convergence cannot be concluded and errors may grow at most at an exponential rate, resulting in the term e ωp in (34).However, an increase of m compensates for this divergent behavior.

Core of the algorithm
Our ansatz for the approximation of invariant fiber bundles is based on computing the set S k n (m, p) numerically.The core of this approach is the calculation of the set Z k n (m, p), i.e. of the kernel of P u n+m+p,k • Ψ(n + m + p, n).Thus, one first has to determine dichotomy projectors, for which good algorithms are available as described in Section 2.5.
The second task is the detection of the zero-contour.It turns out that standard tools like Newton's method and pseudo arclength continuation for finding and continuing zeros of P u n+m+p,k • Ψ(n + m + p, n) are hardly applicable.This is due to the fact that a large number (m+p) of iteration steps, results in extremely steep gradients.Further difficulties occur, if fiber bundles have the form of several components that in addition may intersect.Various noninvertible examples show these characteristics.Finding seed points in each component -as starting points for a numerical continuation -is a nontrivial task.
We avoid these problems by using contour algorithms in Matlab that are based on the mean-value theorem.They apply to functions f : Ê r → Ê, where r ∈ {2, 3}.Their application requires the evaluation of f (x) for x ∈ Γ, where Γ ⊂ Ê r is a grid in some compact area of interest.The resulting discrete data are used for detecting zero-contours.More precisely, the Matlab command contour approximates zero-curves in case r = 2 and for r = 3, isosurface computes zero-surfaces.Note that S k n (m, p) ⊂ Z k n (m, p) and thus, one additionally has to verify for x ∈ Z k n (m, p), whether Ψ(n + i, n)x ∈ U n+i for i ≥ p. Fortunately, this extra test can be neglected in several applications, see the discussion in [39, Section 3.1].
In the next sections, we overcome the restriction that the proposed method applies to two-and three-dimensional models, only.Furthermore, we introduce an extended version, allowing the calculation of hierarchies of stable fiber bundles.

Application to higher-dimensional problems
For higher-dimensional mappings f : Ê r → Ê, zero-hyper-surfaces are difficult to illustrate and their approximation turns out to be rather expensive due to the curse of dimension.As we will see, the intersection of this hyper-surface with an s-dimensional subspace can be visualized easily for s = 2, r ≥ 2 and for s = 3, r ≥ 3. Setting s = r results in cases that have been discussed in Section 3.5.

The Lorenz manifold in a two-dimensional subspace
We illustrate this idea for s = {2, 3}, r = 3 by computing the Lorenz manifold, i.e. the stable manifold of the fixed point 0 of the Lorenz system, cf.[43 (35) Let Ψ be the solution operator of a discretization of ( 35) with the classical Runge-Kutta scheme of order 4 and step size h = 0.02.The fixed point 0 has two stable and one unstable eigenvalue.Stable dichotomy projectors in the resolvent interval R 2 possess the stable respectively unstable subspace of DG(0) as range and kernel.For s = 3 we define f (x) = (Ψ(m, 0)x) |R(P u m,2 ) .The left diagram in Figure 7 shows the zero-surface for m = 100 that we compute, using the Matlab command isosurface.The right diagram shows an approximation of the intersection of the Lorenz manifold with the (x 2 , x 3 )-plane in high resolution.
This task corresponds to the case s = 2 and we define the function f : .
The zero-contour in Figure 7 (right) is calculated with m = 200 iteration steps of each grid point and by applying the Matlab command contour.For generating Figure 7, we use a 400 × 400 × 400 respectively a 2000 × 2000 grid.

A four-dimensional map
Our second example is a four-dimensional autonomous model with an explicitly given three-dimensional stable manifold: Here, {λ 1 , λ 2 } defines the Sacker-Sell spectrum w.r.t. the fixed point ξ = 0.In this example, the stable subspace is not aligned to the coordinate axes and we obtain the graph representation of the stable manifold, see (29), as , where Fix λ 1 = 1 2 , λ 2 = 2.We consider the case s = 3, r = 4 and calculate in Figure 8 intersections of The case s = 2, r = 4 is shown in Figure 9 by intersecting H 2 with selected coordinate planes.All diagrams are computed with the algorithm

Computing hierarchies of stable fibers
In the previous sections, the calculation of zero-contours of mappings f : Ê r → Ê has been considered.Standard tools like the Matlab commands contour and isosurface solve this problem in case r = 2 and r = 3, respectively.For higherdimensional systems, we discussed in Section 3.6 the computation of two-and three-dimensional projections.Approximation techniques for hierarchies of stable fibers require the calculation of zero-contours of mappings f : Ê r → Ê j for j = 1, . . . ,r−1.These contours are nested and the case j = 1 is efficiently solvable.Calculating these contours one by one results in an algorithm for computing the whole stable hierarchy.Note that efficient software is currently not available for the direct computation of contours and isosurfaces of f : Ê m+j → Ê m if m ≥ 2 and j ∈ {1, 2}.
The ranks of two consecutive projectors from (8) may vary by more than one.In this case, we artificially insert further projectors into the sequence (8).This results in a new nested sequence of projectors for n ∈ In case d = ℓ we obtain P u n,j = P u n,j and the rank of these projectors equals j − 1 as required.
We define for j = 1, . . ., d + 1 the family of functions as well as their zero-contours that are nested, cf.(33).
Note that N 1 = Ê d and N 2 is the zero-contour of f 2 : Ê d → Ê. Inductively, we proceed as follows.Assume that N j has already been computed.This set In our applications, we find the parametrization g j via numerical interpolation.Using R( P u n,j+1 ) ∩ R(I − P u n,j ) = R( P u n,j+1 − P u n,j ), we conclude that where ) .This set is numerically accessible, since f j,j+1 •g j : Ê d−j+1 → Ê.In higher-dimensional models, we apply the approach from Section 3.6 throughout these computations.In a nutshell, the nested sequence and Z k n (m, p) respectively S k n (m, p) (see the discussion at the end of Section 3.5) approximate the stable fiber H k n for k ≥ k thanks to Theorem 3. Thus, the sets N ζ(k) , k = k, . . ., ℓ + 1 establish the desired approximation of the stable hierarchy.

Three examples with explicitly given hierarchies
For an illustration of the numerical procedure from Section 3.7, we construct three autonomous models in Ê 3 with explicitly given hierarchies w.r.t. the fixed point 0.
The first one is invertible and has a three-dimensional stable manifold, i.e. k = 1, with k from Section 3.2.In the second model, k = 2 and the Jacobian at the fixed point has a zero-eigenvalue.The third model possesses a two-dimensional unstable manifold, i.e. k = 3.

An invertible system with k = 1
We consider the system Since this system is autonomous, fiber bundles are invariant manifolds that do not depend on time.We aim for an approximation of the hierarchy of stable manifolds of the fixed point 0. The variational equation exhibits the Sacker-Sell spectrum Σ ED = { 1 4 , 1 2 , 3 4 } and the resolvent intervals are 2 ) and R 4 = (0, 1  4 ).Here, ℓ = 3 and k = 1.Note that the ranges of dichotomy projectors -which are eigenspaces of DF (0) in this autonomous model -are not aligned to the coordinate axes.Eigenvectors w.r.t. the eigenvalues 1  4 , 1 2 , 3  4 are All manifolds from the stable hierarchy have global graph representations (29) that are given as with These sets are nested, i.e. (31) holds.Furthermore, we find for H 2 and H 3 the alternative representations Computing the zero-contour N 2 of f 2 , defined in (37), we obtain an approximation of H 2 .In Figure 10 (left) we choose a 300 × 300 × 300 grid and m = 10 iteration steps within the contour algorithm.Then, we interpolate the numerical data w.r.t. the tangent space span{v 1 , v 2 } and obtain g2 as an approximation of g 2 .The graph of g2 is shown in the right panel of Figure 10.
For calculating an approximation of H 3 , we use the proposed method from Section 3.7.We choose a 2000 × 2000 grid in [−4, 4] × [−4, 4].Then we compute the set N 3 as defined in (38) (right panel of Figure 10).For this task, the Matlab command contour applies.
Finally, we discuss approximation errors.Starting with H 2 , we plot in the left diagram of Figure 11 pointwise errors in case m = 10.It turns out that an increase of m does not result in smaller errors as suggested by our error estimate (34).Observed errors are caused by linear interpolation that is used for both, computing the zero-contour and for interpolating the data to obtain g2 .As a consequence, errors are of magnitude O(∆ 2 ), where ∆ is the spacing of the grid, used for the contour algorithm.
Interestingly, the computation of H 3 does not show this problem.Interpolation errors turn out to be negligible when calculating the zero-contour N proj 3 := {y ∈ Ê 2 : f 2,3 (g 2 (y)) = 0} ⊂ Ê 2 , cf. (38).Here, it is essential not to map the resulting points onto H 2 by computing N 3 ≈ g2 (N proj 3 ), since this operation will  Dichotomy projectors within the resolvent interval R ℓ+1 = R 3 are nontrivial for the next model.Note that the resolvent interval R 3 corresponds to the eigenvalue 0, see Figure 1 (lower case).We consider The manifolds H 2 and H 3 have the explicit representations We numerically determine the stable manifold H 2 as well as the strong stable manifold H 3 as described in Section 3.8.1 and analyze approximation errors of H 3 in Figure 12.Due to the zero-eigenvalue of DF (0), we conclude that α s ℓ+1 = ∞ as discussed at the end of Section 2.6.As a consequence, we expect from (34) that no (relevant) errors occur, when calculating H 3 numerically.This confirms the error plot in Figure 12.The next model illustrates that one-dimensional stable manifolds in three-dimensional systems are numerically accessible.Let The fiber L 2 and the stable manifold H 3 have the explicit representations We apply the contour algorithm as described in Section 3.8.1 for obtaining approximations of L 2 and H 3 in Figure 13.Here, it is important to note that the graph representation of L 2 contains a weak unstable manifold, which is not unique.By considering this weak unstable manifold as a quasibounded trajectory forward in time, see

An noninvertible example
We consider the noninvertible autonomous system The stable manifold H 2 of the fixed point 0 does not possess a global graph representation, see Figure 16.For its calculation, we apply the contour algorithm on a 400 × 400 × 400 grid in [−4, 4] × [−2, 2] × [−2, 0.5] and iterate each grid point for 6 steps.Next, we aim for an approximation of the strong stable manifold H 3 .The calculation of (38) requires local parametrizations of H 2 .For constructing them, we first compute the primary part of H 2 , using a cutoff version of F The resulting function exhibits almost global parametrizations.We parametrize With respect to these parametrizations, we determine zero-contours as proposed in (38).Carrying out the contour algorithm, each point on a 2000 × 2000 grid in the blue respectively green rectangle from Figure 14 is mapped to H 2 , using the corresponding parametrization.Then these points are iterated for 4 steps.The resulting data are used to compute the zero-contours in Figure 15.
These contours are mapped to H 2 in Figure 16.Exploiting the symmetry property we obtain data for both, positive and negative values of x 2 from the right contour in Figure 15., since all elements from H 3,1 converge to 0 with the exponential rate log 1 4 and F (x) ∈ H 3,1 for all x ∈ H 3,2 , see Figure 17.Note that Figure 16 contains further parts of the strong stable manifold that are mapped to H 3,2 .
x 1 x 1 x 2 x 2 x 3 x 3 To this ODE-model, we apply the classical Runge-Kutta scheme of order 4 with step size 0.02.Approximations of the stable fiber bundle H 2 t of the fixed point 0 are computed on a 300 × 300 × 300 grid in [−30, 30] 3 .The contour algorithm iterates each grid point for 100 Runge-Kutta steps.
For finding an approximation of the strong stable fiber H 3 t , we proceed as in Section 3.9 and calculate the principle part of the stable fiber using cutoff techniques.Then, we evaluate (38), map the resulting contours to the approximation of H 2 t and obtain the hierarchy of fiber bundles, shown in Figures 18, 19.In Figure 20, we recompute the fibers H 2 0 and H 3 0 on a larger domain.

Conclusion
The appropriate spectral concept for nonautonomous linear systems turns out to be the Sacker-Sell spectrum.Dichotomy projectors and thus, spectral bundles are numerically accessible.Corresponding methods and error estimates are given, including their proofs.
Appropriate reference objects for nonlinear systems are bounded trajectories and we aim for an approximation of their stable hierarchies.The proposed variant of the contour algorithm solves this task by computing fibers from the hierarchy one after another.In principle, restrictions on the dimension of the system and on the dimension of the stable fiber are not required.
If the system (1) is invertible, then our approach also applies to the inverted system and provides the stable hierarchy of (46), which is the unstable hierarchy of the original system (1).
Hierarchies of fiber bundles for ODE models can be computed by applying our method to a time h-map that we obtain, for example, using a one-step discretization scheme.In particular, a negative step size results in an approximation of the unstable hierarchy.
It is an interesting problem to consider the difference equation ( 1) on a finite time interval only, resulting in a finite time dynamical system.Corresponding linear and nonlinear results as well as finite time versions of invariant fibers have been developed, see for example [7,6,21,29,30,31].However, fibers in finite time are fat objects.One may interpret approximations that are computed with the contour algorithm as parts of fat finite time fibers.But it is still an open question whether the whole finite time object is accessible in this way.
P u n = I − P s n , P s,u n Φ(n, m) = Φ(n, m)P s,u m for all n, m ∈ , n ≥ m. (5) Denote by R(P ) the range of a projector P .If the operator Φ(n, m) |R(P u m ) : R(P u m ) → R(P u n ), n ≥ m is invertible, then the same symbol Φ(m, n) : R(P u n ) → R(P u m

Figure 6 :
Figure 6: Illustration of strong and weak stable fibers in nonlinear systems.
p) for all m, p ∈ AE and all n ∈ .Along the lines of [38, Theorem 3.3] we conclude upper-semicontinuity w.r.t. the Hausdorff semi-distance dist(A, B) = sup a∈A inf b∈B a − b .Theorem 3 Assume (A1)-(A5) and fix k ≥ k.Denote by β s,u k the dichotomy rates w.r.t. the resolvent interval R k , cf. (28).For any −β s k < − βs k < βu k < β u k there exist a neighborhood U of ξ , a C > 0 and ω ∈ Ê such that dist(S k n (m, p), H k n ) ≤ Ce ωp e −( βs k + βu k )m

Figure 7 :
Figure 7: Approximation of the Lorenz manifold (left) and of its intersection with the (x 2 , x 3 ) plane (right).

Figure 8 :
Figure 8: Intersection of the stable manifold of (36) with three-dimensional cubes.

4 Figure 9 :
Figure 9: Intersection of the stable manifold of (36) with two-dimensional coordinate planes.

3 Figure 16 : 3
Figure 16: Approximation of the stable manifold H 2 of (43) w.r.t. the fixed point 0 (red ball).The red lines are parts of the strong stable manifold H 3 .

Figure 17 : 2 − x 1 ) 3 
Figure 17: Illustration of parts of the strong stable manifold of the fixed point 0.

Figure 18 :Figure 19 : 3 Figure 20 :
Figure 18: Approximation of the stable and strong stable fiber (in red) of the fixed point 0 (red ball) for the nonautonomous Lorenz system (45) at time t = 0.