Discretization, inﬂation and perturbation of attractors ?

. The discretization of attractors for autonomous and nonautonomous systems is considered. Unlike the autonomous case, where most basic issues are now well understood, the nonautonomous case still has many open questions, which will be discussed here.


Introduction
The basic issues concerning the effect of discretization or perturbation on autonomous attractors are now quite well understood. For nonautonomous systems matters are, however, considerably more complicated as solutions now depend explicitly on both the initial and the current time, so limiting objects need not exist in current time or be invariant, the semigroup evolutionary property no longer holds, and the concept of an attractor for autonomous systems is generally too restrictive.
Nonautonomous systems are ubiquitous. They are easily obtained by including time variation in the vector field of an autonomous differential equation and also arise naturally without an underlying autonomous model. Moreover, they cannot be entirely avoided when one is interested primarily in a particular autonomous system, since perturbations and noise terms are more realistically time dependent, while numerical schemes with variable step size are essentially nonautonomous difference equations even when the underlying differential equation is autonomous.
This Chapter begins with a brief review of results for the autonomous case and more recent ideas on inflated autonomous attractors. The cocycle formalism for a nonautonomous system and the concepts of pullback convergence and pullback attractors in such systems are then outlined. Results on the existence of pullback attractors and of Lyapunov functions characterizing pullback attractors are presented, the formulation of a numerical scheme with variable time steps as a discrete time cocycle system is discussed and the comparison of numerical and original pullback attractors considered, at least in special cases, along with the inflation of pullback attractors. Finally, some open questions and desirable future developments are mentioned.

Autonomous dynamical systems
The solution x(t) = x(t; x 0 ) with initial value x(0; x 0 ) = x 0 of an autonomous differential equationẋ generates a continuous time semigroup φ = {φ t } t∈I R + on IR d defined by φ t (x 0 ) := x(t; x 0 ) for each t ≥ 0 and x 0 ∈ IR d under assumptions on the vector field F that ensure the existence, uniqueness and global extendability of all such solutions. In particular, to simplify the exposition, it will be assumed here that F satisfies a uniform Lipschitz condition on IR d with Lipschitz constant K.
Recall that the Hausdorff separation of nonempty compact subsets A and B of IR d is defined by The long term dynamical behaviour of a semidynamical system φ often occurs in or near its maximal attractor, that is, a nonempty compact subset A 0 of IR d which is φ-invariant, i.e. with φ t (A 0 ) = A 0 for all t ≥ 0, and attracting, i.e. with lim t→∞ H * (φ t (D), A 0 ) = 0 for any bounded subset D ⊂ R d . ( The existence of a maximal attractor follows from that of geometrically simpler and more easily found absorbing sets. A positively invariant compact subset B of R d , i.e. with φ t (B) ⊆ B for all t ≥ 0, is called an absorbing set for the semidynamical system φ on IR d if for every bounded subset D of IR d there exists a t D ∈ IR + such that φ t (D) ⊂ B for all t ≥ t D . The maximal attractor is then given uniquely by A maximal attractor is uniformly asymptotically stable [29] and, as shown by Yoshizawa [31], there then exists a Lyapunov function V : 2. there exist continuous strictly increasing functions α, β : with α(0) = β(0) = 0 and 0 < α(r) < β(r) for all r > 0 such that 3. there exists a constant c > 0 such that Such Lyapunov functions are a very convenient tool for discretization and perturbation investigations as they do not require explicit knowledge of the solutions of the differential equation. For example, the inequality is satisfied [19] by a pth-order one-step numerical scheme (possibly implicit) with constant step size h > 0 applied to the differential equation (1), where is the local discretization error with constant C p . Similarly, the inequality is satisfied by a solution y(t; y 0 ) of the perturbed differential equation with uniformly bounded continuously differentiable perturbations g satisfying g(y) ≤ 1 on IR d . These Lyapunov inequalities can be used to show the existence of absorbing sets for the the discrete time semidynamical system generated by the numerical scheme (6) and for the continuous time semidynamical system generated by the perturbed differential equation (7). From this follows the existence of a maximal numerical attractor A h num and maximal perturbed attractor A h pert , which converge upper semicontinuously to and similarly for the perturbed attractor A h pert [12,19,29,24].
In general, the Hausdorff separation H * above cannot be replaced by Hausdorff metric H, so the numerical attractor A h num or the perturbed attractor A h pert may approximate in the near limit only a proper subset of the original attractor A 0 , which represents a collapse of the original attractor A 0 under discretization or perturbation. For example, the closed unit disc A 0 = {(x, y) ∈ IR 2 : x 2 + y 2 ≤ 1} is the maximal attractor of the two-dimensional system and a disc of radius slightly larger than 1 is the maximal attractor of the explicit Euler scheme applied to (9), while the singleton set A h num = {(0, 0)} is the maximal attractor of the corresponding implicit Euler scheme when the step size h is sufficiently small.

Inflation of the attractor
The totality of possible elements of such discretized or perturbed attractors can be determined [17] by inflating the vector field of the differential equation (1) to form a differential inclusion or setvalued differential equation The set F (x) here is nonempty, compact and convex, and depends continuously on , while the mapping x → F (x) satisfies a uniform Lipschitz condition on IR d with the same Lipschitz constant K as the function F . These properties ensure the existence [2] of an absolutely continuous solution with initial value x(0) = x 0 satisfying Moreoever, the setvalued mapping (t, is the attainability set formed by all such solutions, is continuous with respect to the Hausdorff metric, while Φ t (x 0 ) is a nonempty compact connected subset of where φ t (x 0 ) is the solution of the singlevalued differential equation (1). The solutions can also be shown to satisfy a Lyapunov inequality like (7) with the parameter h replaced by , which can then be used to construct an absorbing set and hence to establish the existence of a maximal attractor A infl for the setvalued semidynamical system [30] generated by the Φ t on IR d . The attractor A infl was called the -inflated attractor [17] of the original singlevalued semidynamical system φ. By the construction, A infl contains A 0 and converges continuously rather than just upper semicontinuously to A 0 , i.e.
If the step size or perturbation parameter h in the numerical scheme (6) or perturbed differential equation (8) is chosen small enough compared with , then the numerical and the perturbed dynamics will be contained within and carried along by the inflated setvalued dynamics Φ t (x 0 ) and hence the numerical attractor A h num and the perturbed attractor A h pert will be contained in the -inflated attractor A infl . The effects of roundoff error, which usually vary from step to step so the actual numerical dynamical system generated within the computer will be nonautonomous even if a constant time step h is used, will similarly be contained in the inflated attractor A infl provided is larger than the machine precision. The inflated attractor A infl is thus the smallest set containing all possible limiting behaviour or approximate autonomous attractors or nonautonomous attractor components (to be defined later) resulting from all possible perturbations and approximations of appropriate magnitude of the original semidynamical system φ. In particular, there is no loss of information in the inflated attractor about the original asymptotic dynamics as may occur with certain approximate systems for which the approximate attractors converge only upper semicontinuously to the original maximal attractor A 0 .

Convergence rates
The theorems used above giving the upper semi continuous convergence of the numerical and perturbed attractors to the original one have a rate of the form α −1 (h p ), where α the strictly increasing function that bounds the Lyapunov function V from below and is usually not known explicitly in practice. To be able to say something more specific about the convergence rate, one needs to know or assume something more about the attractor A 0 .
For example, the -inflated attractor A infl = [− 1/ρ , 1/ρ ] of the scalar differential equationẋ = −x|x| ρ−1 , where ρ ≥ 1, converges to the maximal attractor A 0 = {0} with order 1/ρ. Essentially, the rate of convergence here depends on how fast the unperturbed attractor attracts its neighbourhoods. This is, in fact, typical of the general situation, as was shown in [9] (see also [10]) using a different kind of perturbation that is, however, equivalent to the inflated dynamics of [17].
Let A 0 be the maximal attractor of the semidynamical system generated by (1). A family of forward invariant compact sets {B µ , µ ≥ 0} that depend continuously on µ with respect to the Hausdorff metric H and satisfy A 0 ⊂ int B 0 is called a contracting family of neighbourhoods if there exist a T > 0 with The existence of a contracting family of neighbourhoods with rate of contraction γ is both necessary and sufficient for the rate of convergence γ of the inflated attractor [9].

Theorem 1.
Let B be an absorbing set for a maximal attractor A 0 for which A 0 ⊂ int B. Then A 0 admits a contracting family of neighbourhoods B µ with B 0 = B and contraction rate γ if and only if there is an * > 0 such that there exists an inflated attractor for some constant K.

Nonautonomous dynamical systems
Suppose that a unique solution x(t) = x(t; t 0 , x 0 ) of a nonautonomous differential equationẋ with initial value x(t 0 ; t 0 , x 0 ) = x 0 at time t 0 exists for all x 0 ∈ IR d and t ≥ t 0 ∈ IR. The semigroup property of solutions of an autonomous differential equation now becomes for all x 0 ∈ IR d and all t 0 ≤ t 1 ≤ t 2 in IR, which is called a cocycle property. An abstract nonautonomous dynamical system that is sometimes called a process [12] can be defined in terms of the solution mapping (t, t 0 , x 0 ) → x(t; t 0 , x 0 ) with this generalized semigroup property together with initial condition and continuity properties. An alternative formulation due to Sell [28] that retains the semigroup representation is somewhat more abstract but includes more information about how they evolve in time. It is based on the fact that whenever x(t) is a solution of the differential equation (11), then the x τ (t) := x(τ + t) with fixed τ satisfies the nonautonomous differential equation Denote by F a set of functions F : IR × IR d → IR d such that F τ (·, ·) := F (τ + ·, ·) ∈ F for all τ ∈ IR; for example, F is a compact metric space for almost periodic differential equations. Then introduce a group of shift operators θ τ : F → F by θ τ F := F τ for each τ ∈ IR, define X = IR d × F and write x(t; x 0 , F ) for the solution of (11) with initial value x 0 at initial time Then the family of mappings {Ψ t , t ∈ IR} is a continuous-time semigroup on the state space X and with an appropriate topology on F so that (t, which is also a cocycle property.

Cocycle formalism
The shift operators in the cocycle property (13) can be considered as a driving mechanism that indicates how the dynamics of the nonautonomous system changes with time. This motivates the following definition of an abstract nonautonomous dynamical system. Let θ = {θ t , t ∈ IR} be a group of mappings on a nonempty parameter set P , that is, θ t : P → P with θ 0 = id. and θ t • θ s = θ t+s for all t, s ∈ IR, and write θ t p for θ t (p).

Definition 2.
A family of mappings φ (t,p) : IR d → IR d for t ∈ IR + and p ∈ P is called a cocycle on IR d with respect to a group θ of mappings on P if for all t, s ∈ IR + and p ∈ P .
The use of a general parameter set P here may seem an unnecessary abstraction, but in fact allows for broader applicability and richer dynamical behaviour, particularly when P is a compact metric space. In the skew-product formalism above P is the function space F, while for a periodic differential equation (i.e. with F (t + T, x) = F (t, x) in (11)) a circle S 1 ∼ = IR (mod T ) representing the fundamental periodic interval can be used as P . A control systemẋ(t) = f(x(t), u(t)) can be formulated as a cocyle with a compact metric space of all measurable control functions taking values in a given compact and convex set as the parameter set P . A general nonautonomous differential equation (11) can also included in this new formalism with P being the set IR of initial times and the shift operators θ t by θ t t 0 := t 0 + t, but the parameter space P is now no longer compact. The parameter space P may not even be a topological space, as happens with random dynamical systems for which a canonical probabilistic sample space is used as the parameter space [1,16].

Nonautonomous attraction and attractors
A nonautonomous differential equation (11) can sometimes have an attractor as defined for autonomous systems. For example,φ(t) ≡ 0 is a solution of (11) if the vector field satisfies F (t, 0) = 0 for all t ∈ IR and could be asymptotically stable in the sense of Lyapunov [31]. However, even in this simple case, the rate of attraction and absorbing sets need not be uniform in time, as can be seen from the exampleẋ = −2tx, which has the asymptotically stable solution u(t) ≡ 0 and general solutions x(t; t 0 , x 0 ) = x 0 e −t 2 +t 2 0 . The situation is more complicated for a nonzero asymptotically stable solutionφ. Of course, the time varying change of coordinates z(t) = x(t) −φ(t) will convert this to the preceding situation providedφ is known explicitly. If not, how can a specific pointφ(t) ∈ {φ(s), s ∈ IR} be determined analytically or numerically for a given finite t ∈ IR? The exampleẋ = −x + g(t) for a continuous function g : IR → IR, which has the general solution gives some insight here. Holding t fixed and letting the initial time t 0 → −∞ gives the limit provided that the improper integrals here exist and are finite for each t ∈ IR; see Figure 1. Note thatφ is a solution of the differential equation here that exists for all t ∈ IR and although the geometric trajectory set {φ(s), s ∈ IR} is not invariant under the nonautonomous dynamics, the solution satisfies a dynamical invariance of the form The convergence in (14) can be rewritten as and is called pullback convergence to distinguish it from the usual forwards convergence given by The idea of pullback convergence has been used in other contexts for many years, for example by Mark Krasnosel'skii [25] in the 1960s to establish the existence of solutions of (11) that remain bounded for all t ∈ IR. To help understand what it means, recall that in an autonomous system convergence with time t → ∞ gives the same result as convergence with the elapsed time t − t 0 → ∞ with t fixed and t 0 → −∞, since autonomous dynamics depend only on the elapsed time and the attractor or limit set exists for all time and is invariant. In a nonautonomous system pullback convergence involves essentially t − t 0 → ∞ with t fixed and t 0 → −∞, and thus differs from the usual forward convergence with t − t 0 → ∞ for fixed t 0 . In general, pullback convergence and forwards convergence are independent concepts in nonautonomous systems, as the examplesẋ = −2tx andẋ = 2tx show, since, as the Figures 2 and 3 indicate, the first is forwards but not pullback convergent, whereas the latter is pullback but not forwards convergent.

Pullback attractors
The above observations suggest that a nonautonomous attractor could be defined in terms of pullback convergence, with such a pullback attractor consisting of a family of compact subsets that are mapped into each other under the forward action of the cocycle mappings. A family A = {A p , p ∈ P } of compact subsets of IR d is called a pullback attractor of a cocycle {φ (t,p) , t ∈ IR + , p ∈ P } on IR d if it is invariant in the sense that and pullback attracting in the sense that For example, in (14) above with P = IR, p = t 0 and θ t t 0 = t 0 + t, the component sets A t0 = {φ(t 0 )} for each t 0 ∈ IR form a pullback attractor. The definition also includes the usual autonomous attractor by representing the autonomous semigroup as a cocycle with respect to a singleton parameter set P = {p}. The existence of a pullback attractor also follows from that of more easily found absorbing sets, but these are now defined in terms of the pullback action and families of parametrized sets are used to allow for nonuniformities that are ubiquitous in nonautonomous systems. A family B = {B p , p ∈ P } of compact subsets of IR d is called a pullback absorbing family for a cocycle {φ (t,p) , t ∈ IR + , p ∈ P } on IR d if for each p ∈ P and every bounded subset D of IR d there exists a t D (p) ∈ IR + such that B is said to be uniformly absorbing if the t D (p) here do not depend on p.

Theorem 4.
Let {φ (t,p) , t ∈ IR + , p ∈ P } be a cocycle of continuous mappings on IR d with a pullback absorbing family B = {B(p), p ∈ P }. Then there exists a pullback attractor A = {A p , p ∈ P } with components uniquely determined by Proofs of various versions of this theorem in a number of different contexts can be found in [6,8,20,21,24,27]. Although pullback convergence does not in general imply forwards convergence, additonal continuity and compactness assumptions allow one to conclude that all forward limiting behaviour is contained in the union of all of the pullback attractor component sets.

Corollary 5.
Suppose in addition to the assumptions of Theorem 4 that the mappings φ (t,·) (·) : P × IR d → IR d are continuous, P is a compact metric space, the θ t are continuous and B is uniformly absorbing. Then Furthermore [3], if B consists of just a single absorbing set for all p ∈ P , then the autonomous skew-product flow Ψ t (x 0 , p) := φ (t,p) (x 0 ), θ t p) on the extended state space IR d × P has a maximal autonomous attractor A in IR d × P with the sectional structure A = p∈P A p × {p}.

Lyapunov functions for pullback attractors
A pullback attractor can also be characterized by a Lyapunov function [14]. Suppose that the cocycle dynamical system (φ, θ) is generated by a nonau-tonomous differential equatioṅ i.e. with d dt where for simplicity (p, x) → f(p, x) is assumed to be continuous in (p, x) ∈ P × IR d , x → f(p, x) to be uniformly Lipschitz continuous on IR d with Lipschitz constant L(p) for each p ∈ P , and (t, p) → θ t p to be continuous. Let A be a pullback attractor for (φ, θ). Then there exists a pullback neighbourhood system B with A p ⊂ intB p for each p ∈ P such that the function V : where T p,t := t+ t 0 L(θ −s p) ds with T p,0 := 0, satifies the following properties: (1) For each p ∈ P there exists a function a(p, ·) : IR + → IR + with a(p, 0) = 0 and a(p, r) > 0 for all r > 0 which is monotonic increasing in r such that for all x 0 ∈ IR d ; (2) V is uniformly Lipschitz on IR d with Lipschitz constant 1 for all p ∈ P ; (3) For all p ∈ P and any bounded set D in IR d In addition, it can be shown that there exists a family N = {N p , p ∈ P } of nonempty compact sets of IR d which are positively invariant w.r.t. φ in the sense that φ t,p (N p ) ⊆ N θtp for all t ≥ 0, p ∈ P , and satisfying A p ⊂ intN p for each p ∈ P such that for all x 0 ∈ N p and t ≥ 0, which in turn implies that However, this does not imply Lyapunov stability or asymptotic stability, since there is no guarantee (without additional assumptions) that inf j≥0 a(θ t p, r) > 0 for r > 0, so dist(φ t,p (x 0 ), A θtp ) need not become small as t → ∞.
This is in fact what happens with the differential equationẋ = 2tx with solutions x(t; t 0 , x 0 ) = x 0 e t 2 −t 2 0 .
The pullback attractor here has components A t0 = {0} for each t 0 ∈ P = IR and a Lyapunov function meeting the above requirements is given by and Property (2) are immediate, while Property (3) follows from In addition, V satisfies inequality (23), since although from Figure 3 the zero solution is clearly not Lyapunov stable .

Approximation of pullback attractors
A nonautonomous dynamical system arises if variable stepsizes h n are used in the numerical scheme (6) or a timedependent perturbation g(t, y) in the perturbed differential equation (8), even though the original dynamical system generated by the differential equation (1) is autonomous. However, in both cases the Lyapunov inequalites (5) and (7) remain valid and can be used to construct a single uniform absorbing set about the autonomous maximal attractor A 0 for each of the resulting nonautonomous dynamical systems. The nonautonomously perturbed ordinary differential equation obviously generates a cocycle with respect to the parameter set P = IR of initial times t 0 , for which there thus exists a pullback attractor A h pert = {A h pert,t0 , t 0 ∈ IR} and the individual component sets converge upper semi continuously to A 0 [24], i.e.
An analogous result holds for the numerical scheme with variable stepsizes [21,23], but the formulation of such a numerical scheme as a discrete time cocycle is not as obvious.

Numerical schemes as discrete time cocycles
Consider an explicit one-step numerical scheme (6) with variable time-steps, applied to the autonomous differential equation (1). Define H δ to be the set of all two sided sequences h = {h n } n∈I Z satisfying for δ > 0 (the particular factor 1/2 here is chosen just for convenience) and define the shift operatorθ : H δ → H δ byθh =θ{h n } n∈I Z := {h n+1 } n∈I Z . The set H δ is a compact metric space with the metric n and the shift operatorθ is a homeomorphism on this metric space, so its iterations form a discrete time group. It then follows that the numerical scheme (24) with variable time steps generates a discrete time cocycle ψ on IR d with the parameter space H δ and shift operator group defined by for any n ∈ IN , x 0 ∈ IR d and h = {h n } n∈I Z ∈ H δ . As mentioned above, it then has a pullback attractor A h num = {A δ num,h , h ∈ H δ } for which the components converge upper semicontinuously to the autonomous maximal attractor A 0 uniformly in the sense that The situation is somewhat more complicated for the discretization of a nonautonomous differential equation of the form (21) that generates a cocycle on IR d with respect to the given parameter space P and group θ. An explicit one-step numerical scheme with variable step size applied to (20) now takes the form where the times t n are related to a sequence of stepsizes h ∈ H δ by t 0 = 0 and define t n = t n (h) := n−1 j=0 h j and t −n = t −n (h) := − n j=1 h −j for n ≥ 1. Define a mapping ψ : where Q δ := H δ × P and x n is the nth iterate of the numerical scheme (26) with initial value x 0 ∈ IR d , initial parameter p ∈ P and stepsize sequence h ∈ H δ . Then ψ is a dsicrete time cocycle on IR d with the extended parameter space Q δ and the group Θ = {Θ n } n∈I Z on H δ × P with Θ n :Q δ → Q δ for n ∈ I Z defined by iteration of the component shift operators, A numerical pullback attractor now has the form if it exists. The existence of both continuous time and discrete time numerical pullback attractors were established in [3] under very strongly uniform structural assumptions on the vector field of the nonautonomous differential equation (21). Here the upper semicontinuous convergence reads A practical complication here is that a numerical scheme (26) applied to a differential equation of the form (21) may have have a lower order than the scheme on which it is based (e.g. a Runge-Kutta scheme) since the mapping t → F (θ t p, x) may not be sufficiently smooth to justify the usual error estimations. The original higher order may still be retained if one first averages the vector field over each discretization subinterval with an appropriately chosen sampling step [11].

Inflated pullback attractors
Inflating the vector field of the differential equation (21) leads to a nonautonomous differential inclusion or setvalued differential equation of the form where the use of a family := { p , p ∈ P } of inflation parameters is to handle nonuniformities in the nonautonomous vector field. Solutions of this equation are interpreted as absolutely continuous functions x(t) satisfying which requires that the mappings t → θtp must satisfy some kind of continuity property to ensure that the resulting attainability sets Φ t,p (x 0 ) generate a setvalued cocycle mapping [23].
Consider the uniform inflation of the differential equationẋ = 2tx, that is, the differential inclusion which generates the setvalued cocycle t0 e −s 2 ds over the parameter set P = IR with shift θ t t 0 = t 0 + t. The setvalued orinflated pullback attractor A infl = A infl,t0 , t 0 ∈ IR here has components for t 0 ∈ IR, but requires a restriction on the regions of pullback attraction to subsets D t0 = {x 0 ∈ IR 1 : |x 0 | ≤ e t 2 0 √ π} for each t 0 ∈ IR, see Figure  4, where the positive part of the pullback attractor is given by the shaded region and the upper curve indicates the upper bound on D t0 . Such a restriction on the regions of pullback attraction is typical in many examples and the theory of pullback attractors has been extended to handle it [16]. The component sets of any perturbed or numerical pullback attractor for sufficiently close perturbation or numerical approximations will lie within the corresponding component of the inflated pullback attractor. Their regions of pullback attraction will, in general, also need to be parameter dependent or the magnitude of the error need to be made increasingly smaller with increasing t 0 .
6 Scope of the project and future work This article surveys the progress that has been made as well as sketching the background and basic issues involved with a project that has been funded as part of the DANSE research program during the final two years of its six year existence. The original broad aim of the project was to investigate the effects of discretization and perturbation on attractors of nonautonomous dynamical systems and, more specifically, to generalize the 1986 result of Kloeden and Lorenz [19] on the discretization of autonomous attractors, for which the Lyapunov inequality (5) was used to construct an absorbing set for the discretized system. The characterization of pullback attractors by a Lyapunov function was thus seen as the crucial initial step in the project. The longer term motivation for the project was to understand the effects of discretization and perturbation on random dynamical systems, which are intrinsically nonautonomous and, moreover, highly nonuniform.
Results on the discretization of attractors are essentially perturbation results, if for rather atypical types of perturbations, for which the methods that are traditionally applied usually require for some kind of uniformity in the assumed behaviour under consideration. For cocycle systems with compact parameter sets and other nice topological and structural properties, as is assumed for the skew-product flow formalism of nonautonomous deterministic differential eqautions, reasonable progress can be expected and has been made [3]. The example of an inflated pullback attractor in the previous section is only uniform over the past up to any present time rather than for all times so just how much of uniformity with respect to the parameter space is required to give the sought result remains to be seen. Convergence rates for approximations of pullback and inflated pullback attractors also still need to be carefully investigated as does the apparently strong connection between the pullback attractors and controllability properties of the control systems that are representable as cocycle systems.
Measurability rather than continuity is the dominate characteristic of random dynamical systems, at least with respect to the parameter in the cocycle mapping, so some very deep and challenging theoretical analysis seems to be required. At present only simple special cases have been investigated, e.g. [16]. Results of numerical simulations of random dynamical systems, again only for special cases, using the subdivision algorithm of Dellnitz and Hohmann [7] reveal very interesting dynamical behaviour and suggest that such an investigation will be well worth the effort [26].