On the computation of attractors for delay differential equations

In this work we present a novel framework for the computation of finite dimensional invariant sets of infinite dimensional dynamical systems. It extends a classical subdivision technique [Dellnitz/Hohmann 1997] for the computation of such objects of finite dimensional systems to the infinite dimensional case by utilizing results on embedding techniques for infinite dimensional systems. We show how to implement this approach for the analysis of delay differential equations and illustrate the feasibility of our implementation by computing invariant sets for three different delay differential equations.


Introduction
Over the last two decades so-called set oriented numerical methods have been developed in the context of the numerical treatment of dynamical systems (e. g. [5,6,11,13]). The basic idea is to cover the objects of interest -for instance invariant sets or invariant measures -by outer approximations which are created via multilevel subdivision techniques. These techniques have been used in several different application areas such as molecular dynamics ( [23]), astrodynamics ( [7]) or ocean dynamics ( [12]).
So far the applicability of the subdivision scheme is restricted to finite dimensional dynamical systems, i. e. ordinary differential equations or finite dimensional discrete dynamical systems. In this paper we extend this technique to the infinite dimensional context. More precisely, we develop a set oriented numerical technique which allows us to compute low-dimensional invariant sets for infinite dimensional dynamical systems. Rather than using a straightforward approach based on an appropriate combination of Galerkin expansions and subdivision steps we follow a completely novel path and utilize embedding results in our numerical treatment.
The first result on embeddings in the dynamical systems context is the celebrated Takens Embedding Theorem [25]. Takens has shown that an invariant set -in his context this has to be a compact manifold of dimension d -can generically be reconstructed using the so-called observation map which consists of observations of the dynamical behavior at an appropriate number (at least 2d + 1) of consecutive snapshots in time. This result has been extended by Sauer et al. in [22] to the context of compact invariant sets of box counting dimension d. There it has been shown that the same observation map can be used for the reconstruction of the invariant set as long as more than 2d consecutive snapshots in time are used. Moreover in this work the notion of "genericity" has been replaced by the more intuitive notion of "prevalence". Finally, in [20] Robinson extended the results obtained in [22] to dynamical systems on infinite dimensional Banach spaces. It turns out that also here the same observation map can be used to reconstruct invariant sets of finite dimensional box counting dimension d. However, in addition to the dimension of the set another quantity comes into play, namely the thickness exponent σ. Roughly speaking this exponent measures how well the invariant set can be approximated by finite dimensional subspaces of the underlying Banach space. The lower bound 2d on the number of snapshots has to be replaced by 2(1 + σ)d accordingly.
We remark that there are several further extensions of Takens' theorem. For instance in [24] forced systems are considered, and in [19] one can find a stochastic version of this result.
Our results in this paper are based on Robinson's embedding theorem. In fact, we will combine the reconstruction based on the observation map with the classical subdivision techniques developed in [5]. Assuming that a bound on the box-counting dimension and the thickness exponent of the invariant set are known we use the observation map and its inverse to define a dynamical system ϕ in the embedding space of dimension k > 2(1+σ)d. Then the subdivision scheme is applied to compute the reconstructed invariant set for ϕ. Observe that this way we can always perform the numerical approximation within a finite dimensional space of fixed dimension k, and this is in contrast to Galerkin based approaches where one would have to extend the expansions in order to improve the quality of the approximation.
The numerical approach we are proposing is in principle applicable to arbitrary infinite dimensional dynamical systems. However, here we will restrict our attention to delay differential equations with constant delay in the numerical realization. The applications to e. g. partial differential equations will be done in future work.
A detailed outline of the paper is as follows. In Section 2 we briefly summarize the infinite dimensional embedding theory introduced in [20]. In Section 3 we employ the main result of [20] for the construction of a numerical approach for the computation of compact finite dimensional attractors of infinite dimensional dynamical systems. First we construct a continuous dynamical system ϕ on the embedding space using a generalization of the well-known Tietze extension theorem [9, I.5.3] which is due to Dugundji [8,Theorem 4.1]. Then we extend the results from [5] to the situation where the underlying dynamical system ϕ is just continuous (and not homoemorphic). A numerical realization for the computation of attractors of delay differential equations is given in Section 4. Finally, in Section 5, we illustrate the efficiency of our novel approach for three different delay differential equations.

Infinite Dimensional Takens Embedding
We start with a short review of the contents of [20]. We consider a dynamical system of the form where Φ : Y → Y is Lipschitz continuous on a Banach space Y . Moreover we assume that Φ has an invariant compact set A, that is Later on we will additionally assume that A is a global attractor in the sense that it attracts all bounded sets within Y as t → ∞.
For the statement of the main result of [20] we need three particular notions: prevalence [22], upper box counting dimension and thickness exponent [15].
Definition 2.2 (Upper box counting dimension). Let Y be a Banach space, and let X ⊂ Y be compact. For ε > 0, denote by N Y (X, ε) the minimal number of balls of radius ε (in the norm of Y ) necessary to cover the set X. Then denotes the upper box-counting dimension of X. is called the thickness exponent of X in Y .
From this definition one sees that essentially, σ(X, Y ) captures how well X can be approximated from within finite dimensional subspaces of Y . In [17], as a more practical expression Kukavica and Robinson prove that where ε Y (X, n) is the minimum distance between X and any n-dimensional linear subspace of Y . These notions are essential in answering the question when a delay embedding technique applied to an invariant subset A ⊂ Y will work generically. More precisely, the results are as follows.
Theorem 2.4 ( [15]). Let Y be a Banach space and X ⊂ Y compact, with upper box counting dimension d(X; Y ) =: d and thickness exponent σ(X, Y ) =: σ. Let N > 2d be an integer, and let α ∈ R with Then for a prevalent set of linear maps L : Note that this implies that a prevalent set of linear maps Y → R N will be one-toone on X. Using this theorem, the following result concerning the delay embedding technique can be proven. Theorem 2.5 ([20]). Suppose that the upper box counting dimension of A is d(A) = d, and that A has a thickness exponent σ. Choose an integer k > 2(1 + σ)d and suppose further that the set A p of p-periodic points of Φ satisfies d(A p ) < p/(2 + 2σ) for p = 1, . . . , k. Then for a prevalent set of Lipschitz maps f : is one-to-one on the invariant set A.
Remark 1. Following an observation already made in [22,Remark 2.9], we note that this result may be generalized to the case where several independent observables are evaluated. More precisely, also for a prevalent set of Lipschitz maps f 1 , . . . , f q :

Computation of Embedded Attractors via Subdivision
3.1. Finite-dimensional Embeddings of Attractors. In this section we employ Theorem 2.5 in order to construct a method for the computation of compact, finite dimensional attractors of infinite dimensional dynamical systems on a Banach space where D k is the map defined in Theorem 2.5 and f is chosen such that D k is one-to-one on A.
We now develop a set oriented numerical technique for the approximation of the set A k . First, we define a dynamical system on R k for which A k is an invariant set, on which the dynamics is conjugate to that of Φ on A. For this we define the map ϕ : R k → R k by where E : R k → Y and R : Y → R k are an embedding and a restriction, respectively, that satisfy More concretely, we define the map R = D k by the right-hand side of (3) (see Figure 1).

Remark 2.
Although D k is derived from what would commonly be termed an embedding theorem, we call R a restriction, as it maps from an infinite dimensional domain into a finite-dimensional set, and E an embedding, as it maps (or embeds) a finite-dimensional into an infinite dimensional space.
The map E is obtained in two steps: Firstly, as R is one-to-one on A, the requirement that uniquely defines a mapẼ : A k → A. In a second step, we extend this map to a continuous map E : R k → Y . To do this, we employ a generalization of the well-known Tietze extension theorem [9, I.5.3] found by Dugundji [8, Theorem 4.1], stated here with notation adapted to our needs: Theorem 3.1. Let X be an arbitrary metric space and A ⊂ X be closed, V a locally convex linear space and p : A → V be continuous. Then there is a continuous map P : X → V with P |A = p such that P (X) is contained in the convex hull of p(A).
Using this theorem in the situation introduced above, we obtain the following Proof. By construction, the map R : Y → R k given by (3) is continuous (even Lipschitz) and one-to-one. Thus, restricting R to A we obtain a bijective map R : A → A k . As A is assumed to be compact and A k ⊂ R k is Hausdorff,R is a homeomorphism by a well-known theorem from elementary topology (see e. g. [26,Theorem 17.14]). Thus we obtain a continuous mapẼ : As Y is a normed linear space, it is locally convex. Thus we can apply Dugundji's Theorem with X = R k , A = A k , p =Ẽ and V = Y to see that there is a continuous map E : R k → Y with E |A k =Ẽ. Finally, defining ϕ through (5) we obtain that ϕ is continuous as a composition of continous maps. Now we are in a position to approximate the embedded invariant set A k via the corresponding dynamical system (9) x j+1 = ϕ(x j ) j = 0, 1, . . . .
To this end, we employ a subdivision scheme as defined in [5].
3.2. Subdivision Scheme. We briefly review the classical subdivision procedure. Let Q ⊂ R k be a compact set. We define the global attractor relative to Q by The subdivision procedure allows us to approximate this set. Roughly speaking, the idea of the algorithm is as follows. We start with a finite family of (large) compact subsets of R k which cover the domain in which we want to analyze the dynamical behavior. Then we subdivide each of these sets into smaller ones and throw away subsets which do not contain part of the relative global attractor. Continuing the process with the new collection of (smaller) sets it becomes intuitively clear that this should lead to a successively finer approximation of the relative global attractor.
Let us be more precise. The algorithm generates a sequence B 0 , B 1 , . . . of finite collections of compact subsets of R k such that the diameter converges to zero for → ∞. Given an initial collection B 0 , we inductively obtain B from B −1 for = 1, 2, . . . in two steps.
(1) Subdivision: Construct a new collectionB such that The first step guarantees that the collections B consist of successively finer sets for increasing . In fact, by construction In the second step we remove each subset whose preimage does neither intersect itself nor any other subset inB . As we shall see, this step is responsible for the fact that the unions B∈B B approach the relative global attractor.
Denote by Q the collection of compact subsets obtained after subdivision steps, that is, Moreover let B 0 be a finite collection of closed subsets with Q 0 = B∈B0 B = Q.
Then the main convergence result of [5] states that where h(B, C) is the usual Hausdorff distance between two compact subsets B, C ⊂ R n . However, in that work the authors assume that ϕ is a homeomorphism and not just continuous, as in the situation here. For this reason, in the following we present a proof of convergence for continuous ϕ.

3.3.
Proof of Convergence. Essentially, we will be able to follow the structure of the proof in [5]. However, there are some technical differences, and we will need one additional assumption on A Q . We begin with the following observation: We now can prove our first result.
Proposition 2. Let A Q be the global attractor relative to the compact set Q, and suppose that the embedded attractor A k satisfies A k ⊂ Q. Then Proof. By construction of our dynamical system (see (5)- (7)), we have ϕ(A k ) = A k . Thus, Lemma 3.2 implies that A k ⊂ A Q .
Remark 3. Observe that we can in general not expect that A k = A Q . In fact, by construction A Q may contain several invariant sets and related heteroclinic connections. In this sense A k will be embedded in A Q .
Next observe that the Q 's define a nested sequence of compact sets, that is, Q +1 ⊂ Q . Therefore, for each m, and we may view Q as the limit of the Q 's. Now we will prove the convergence of the subdivision scheme for continuous ϕ. More precisely we will show that Q ∞ = A Q . This will be done in two steps. The first is Proof. We will show that Q ∞ ⊂ ϕ(Q ∞ ). Then the result follows with Lemma 3.2.
Let y ∈ Q ∞ . Then for every ≥ 0 there is a unique B (y) ∈ B with y ∈ B (y). By the selection step of the subdivision scheme (see (13)), there is z ∈ Q with ϕ(z ) ∈ B (y). Choosing a convergent subsequence of (z ), if necessary, we may assume that z = lim →∞ z . By construction, z ∈ Q ∞ , and since lim →∞ diam(B (y)) = 0 we conclude that lim →∞ ϕ(z ) = y. Finally ϕ is continuous, and therefore y = ϕ(z). Hence y ∈ ϕ(Q ∞ ).
For the inverse inclusion, we need to introduce an additional assumption, namely that ϕ −1 (A Q ) ⊂ A Q . This is automatically satisfied in the case where ϕ is a homeomorphism. Moreover if A k is attracting and A k = A Q then A Q is backward invariant. These observations justify this assumption.
This proof is identical to the proof of Lemma 3.2 in [5]. Thus, we will not restate it here.
We summarize the convergence result in the following Proposition 3. Suppose that the relative global attractor

Approximation of Attracting Sets.
Observe that by construction of ϕ, the set A Q defined in Section 3.2 contains the one-to-one image A k of the invariant set A of Φ. We now show that by using sufficiently high powers of Φ we can actually approximate a one-to-one image of A if A is attracting. Therefore, we now assume that A is an attracting set, that is, A attracts all bounded sets within a neighborhood U of A. Moreover we assume that the set Q is chosen in such a way that (17) A k ⊂ Q and E(Q) ⊂ U.
Hence, for every x ∈ Q, Φ j (E(x)) will eventually approach the attracting set A for j → ∞. However, this alone does not guarantee that A k is also an attracting set for the dynamical system ϕ. For instance, it may be the case that for a certainx ∈ Q one has a "spurious fixed point" in the sense that although Φ(E(x)) = E(x) may be closer to A than E(x).
In order to overcome this problem we now define for m ≥ 1 the continuous maps and denote the corresponding relative global attractors by A m Q . Remark 4. Observe that A is an invariant set for Φ m for every m and therefore we can still use R as the restriction in our construction of the dynamical system ϕ m . (17)

), and Lemma 3.2 implies that
Set V = E(Q) ⊂ U by assumption (see (17)). Since A is attracting and V is a compact set within U we can find an m ≥ 1 such that where h is the Hausdorff distance. By our choice of δ it follows that Thus, yielding a contradiction.

Remark 5.
Roughly speaking Proposition 4 states that it will be possible to approximate an attracting set for Φ if we perform the computations with appropriately high iterates of Φ.

Numerical Realization for Delay Differential Equations
As one important setting where finite dimensional dynamical phenomena occur in infinite dimensional Banach spaces, we consider delay differential equations with constant time delay τ > 0. More precisely we consider equations of the form where y(t) ∈ R n and g : R n × R n → R n is a smooth map. Following [14], we denote by C = C([−τ, 0], R n ) the (infinite dimensional) state space of the dynamical system (19). Equipped with the maximum norm, C is a Banach space.
Let y u (t) be the trajectory generated by (19) with the initial condition u ∈ C. Then the flow Φ s : C → C of (19) is given by Next we fix ω > 0 and consider the corresponding time-ω-map Φ ω : C → C as our dynamical system. That is, we set (20) Φ = Φ ω and Y = C in our abstract dynamical system (1). In order to numerically realize the construction of the map ϕ = R•Φ•E described in Section 3, we have to work on three tasks: the implementation of E, of R, and of Φ ω respectively. For the latter we will rely on standard methods for forward time integration of DDEs [2]. The map R will be realized on the basis of Theorem 2.5 and Remark 1 by an appropriate choice of observables. For the numerical construction of the embedding E we will employ a bootstrapping method that re-uses results of previous computations. This way we will in particular guarantee that the identities in (6) are at least approximately satisfied.
From now on we assume that upper bounds for both the box counting dimension d and the thickness exponent σ are available. This allows us to fix k > 2(1 + σ)d according to Theorem 2.5.

Numerical Realization of R.
For the definition of R we have to specify the time span ω and appropriate corresponding observables. In the case of a scalar equation (n = 1) we choose the observable f to be Thus, in this case the restriction R is simply given by The time span ω (see (20)) is defined to be a natural fraction of τ , that is Remark 6.
(a) Observe that a natural choice for K in (21) would be K = k−1 for k > 1. That is, for each evaluation of R the observable would be applied to a function u : [−τ, 0] → R at k equally distributed time steps within the interval [−τ, 0]. (b) As described in Section 3.4 (see Remark 4) we will frequently replace Φ by Φ m (m > 1) in order to speed up the convergence towards the invariant sets A resp. A k . For an illustration of this procedure see Figure 2. Figure 2. Numerical realization of the restriction R for n = 1, K = 2 and m = 6K.
For the numerical analysis of systems of delay differential equations (n > 1) we make use of Remark 1 as follows: For each component u j of u we define a separate observable f j (j = 1, . . . , n) by (22) f j (u) = u j (ν j ) for a ν j ∈ [−τ, 0], and choose different time spans accordingly. Finally, we note that also more general constructions for the restriction R : C → R k can be employed. In fact, by virtue of Theorem 2.4, for any k that is sufficiently large for the delay embedding construction, an arbitrary linear map C → R k will generically be one-to-one on A. Therefore almost every linear combination of trajectory points computed during forward integration can be used for the construction of the map R.

Numerical Realization of E.
In the application of the subdivision scheme for the computation of the relative global attractor A k described in Section 3.2 one has to perform the selection step (13)). Numerically this is realized as follows: At first ϕ is evaluated for a large number of test points z k j ∈ B k for each box B k ∈B . Then a box B m is kept in the collection B if it is hit by (at least) one of the images ϕ(zk j ).

Remark 7.
In practice the test points z k j ∈ B k can be chosen according to several different strategies: In low dimensional problems one can choose them from a regular grid within each box B k . Alternatively one can select the test points from the boundaries of the boxes. In our computations we have sampled a fixed number of test points from each box at random with respect to a uniform distribution.
For the evaluation of ϕ = R • Φ • E at a test point z we need to define the image E(z), that is, we need to generate adequate initial conditions for the forward integration of the DDE (19). In the first step of the subdivision procedure, when no information on A is available, we proceed as follows. In the case of a scalar delay equation, that is n = 1, we construct a piecewise linear function u = E(z), where (24) u(t i ) = z i , Observe that by this choice of E and R the condition R • E(z) = z is satisfied for each test point z (see (6) and Remark 6 (a)).
For n > 1 we proceed analogously and distribute the components of z ∈ R k to the components u j of u = E(z) ∈ R n according to (22) and (23). Also in this case the condition R • E(z) = z still holds.
In the following steps of the subdivision procedure we proceed as follows: Note that if B ∈ B , then, by the selection step, there must have been aB ∈ B −1 such that R(Φ ω (E(ẑ))) ∈ B for at least one test pointẑ ∈B. Therefore, we can use the information from the computation of Φ ω (E(ẑ)) to construct an appropriate E(z) for each test point z ∈ B.
More concretely, in every step of the subdivision procedure, for every set B ∈ B we keep additional information about the points Φ ω (E(ẑ)) that were mapped into B by R in the previous step. In the simplest case, we store k i ≥ 1 additional equally distributed function values for each interval (−τ + (i − 1)ω, −τ + iω) for i = 1, . . . , k − 1. When ϕ(B) is to be computed using test points from B, we first use the points in B for which additional information is available and generate the corresponding initial value functions via spline interpolation. Note that the more information we store, the smaller the error Φ ω (E(ẑ)) − E(z) becomes for z = R(Φ ω (E(ẑ))). That is, we enforce an approximation of the identity E•R(u) = u for all u ∈ A (see (6)).
If the additional information is available only for a few points in B, we generate new test points in B at random and construct corresponding trajectories by piecewise linear interpolation.

Numerical Results
In this section we present results of computations carried out for three different delay differential equations. In each case, u(t) is scalar, although for the DDE considered in Section 5.2 the problem is recast into a three-dimensional form in order to obtain a first-order equation.

The Modified Wright Equation.
As the first example, we consider a modification of the Wright equation, In [14] it has been shown that the stationary solution u 0 (t) ≡ 0 undergoes a supercritical Hopf bifurcation at α = π/2. Thus, (25) possesses a stable periodic solution for α > π/2 -at least locally. In our computations we set α = 2, choose the embedding dimension k = 5, and approximate the relative global attractor A Q ⊂ R 5 for Q = [−2, 2] 5 , see Figure 3. Here the set A Q consists of a reconstruction of the two-dimensional unstable manifold of u 0 ≡ 0 which accumulates on a stable periodic orbit at its boundary.

The Arneodo
System with Delay. The second example is a modification of the Arneodo system [1] where a delay is introduced in the first order derivative of u, This equation has been introduced and analyzed in [21]. In our computations we use the equivalent reformulation as a first-order systeṁ  The undelayed system (i. e. (26) with τ = 0) has been studied extensively. It possesses the equilibria O 1 = (0, 0, 0) and O 2 = (α, 0, 0), the latter is asymptotically stable for α < 2. At α = 2 the equilibrium O 2 undergoes a supercritical Hopf bifurcation (cf. [16]). For values of α which are slightly larger than two, points on the two-dimensional unstable manifold of O 2 converge to the corresponding limit cycle on the branch of periodic solutions. That is, topologically speaking, we have the same situation as in Figure 3.
For the delayed (i. e. τ > 0) equation, the Hopf bifurcation occurs at decreasing values of α for increasing values of τ . For fixed α = 2.5, the amplitude of the limit cycle grows with increasing values of τ and loses its stability in a perioddoubling bifurcation at τ ≈ 0.11 [21]. Our purpose is to investigate the structure of the relative global attractor right after the occurrence of the period-doubling bifurcation. Concretely we set α = 2.5, τ = 0.13, choose the embedding dimension k = 5, and approximate the relative global attractor 4,4]. This way we compute a reconstruction of the twodimensional unstable manifold of the origin which accumulates on a period-doubled limit cycle.
In this example we have made use of Remark 1 in our numerical realization. Concretely we have chosen ω = τ /2 and the following three observables (see (22) and (23)) Thus, the restriction R can be written as Observe that R : C → R 5 is linear and therefore also Theorem 2.4 could be used in order to justify this construction. The corresponding reconstructions of the relative global attractor are shown in Figure 5.
Observe that after the period doubling bifurcation the relative global attractor contains a Moebius strip with the period-doubled periodic solution at its boundary. Thus, in the course of the period doubling bifurcation there has to occur a significant change of the geometry of the unstable manifold at its boundary so that it can accommodate the period-doubled solution. In fact, the corresponding mechanism has been analyzed analytically already in 1984 by Crawford and Omohundro [3]. It turns out that at the period doubling the unstable manifold starts to wrap itself "infinitesimally" around the unstable periodic solution. In a corresponding Poincaré section this becomes a spiraling behavior with very sharp curvature, and we analyze this behavior at the reconstruction in Figure 6 (see also Figure 16 in [3]). However, we expect that one would have to choose a much higher resolution (i. e. higher number of subdivisions) in order to reveal this dynamical behavior more clearly.

The Mackey-Glass Equation.
Our final example is the well-known delay differential equation introduced by Mackey and Glass in 1977 [18], namely where we choose β = 2, γ = 1, η = 9.65, and τ = 2. This equation is a model of blood production, where u(t) represents the concentration of blood at time t,u(t) represents production at time t and u(t − τ ) is the concentration at an earlier time, when the request for more blood is made. Direct numerical simulations indicate that the dimension of the corresponding attracting set is approximately d = 2. Thus, we choose the embedding dimension k = 7, and approximate the relative global attractor A Q for Q = [0, 1.5] 7 ⊂ R 7 . In Figure 7 we show projections of the coverings obtained after = 28, 42 and 63 subdivision steps as well as a direct simulation.
Finally we conclude this section with an outlook and show a corresponding invariant measure for the reconstructed Mackey-Glass attractor in Figure 8. This measure has been computed with the software package GAIO [4] which is based on the techniques developed in [6]. However, a detailed investigation on how to  approximate invariant measures in infinite dimensional problems efficiently using embedding theory will be done in future work.

Conclusion
In contrast to the situation for finite dimensional dynamical systems, for which there exists a wide range of advanced numerical tools, there are currently only few options besides direct simulation for the numerical computation of attractors of infinite dimensional dynamical systems generated e. g. by DDEs. In this paper we develop a general methodology for the computation of finite dimensional compact attractors of infinite dimensional dynamical systems, and illustrate its application to several DDEs.
Combining the delay embedding technique with a set-oriented method for the computation of attractors, we obtain a flexible method for the analysis of infinite dimensional dynamical systems. More concretely, we use standard techniques for short-time simulation of the system in question to approximate the infinite dimensional dynamics, and employ the delay embedding technique on simulation results in order to obtain a representation of the dynamics by a continuous map on a (c) Zoom into the upper left corner of (b) (d) Zoom into the lower right corner of (b) Figure 6. Detailed illustration of a Poincaré section through the relative global attractor for (26).
moderately sized space. This map, in turn, can be analyzed using the subdivision scheme in order to compute a covering of an attractor. We show that in this way one obtains sets that are in one-to-one correspondence with the infinite dimensional attractor. The method proposed shares vital characteristics with its counterparts for finite dimensional systems. Most importantly, the numerical effort essentially depends on the dimension of the object to be computed, and not on the dimension of ambient space used for computations, that is in the case of this paper, the dimension of the space used for the delay embedding. In three examples, we have illustrated the suitability for the computation of two-dimensional attractors. (In the case of the Mackey-Glass equation, the attractor probably even has a box counting dimension larger than d = 2.6, see e.g. [10].) Furthermore, due to the set oriented nature of the underlying subdivision algorithm, the method does not depend on any special geometric properties (besides compactness) of the attractor.