Computing and Visualization in Science manuscript No. (will be inserted by the editor)

An adaptive subdivision technique for the approximation of attractors and invariant measures

Michael Dellnitz^{⋆}, Oliver Junge

Mathematisches Institut, Universität Bayreuth, D-95440 Bayreuth, Germany

(http://www.uni-bayreuth.de/departments/math/~mdellnitz; http://www.uni-bayreuth.de/departments/math/~ojunge)

Received: 13 January 1997 / Accepted: 29 September 1997

Communicated by: P. Deuflhard

Abstract. Recently subdivision techniques have been introduced in the numerical investigation of the temporal behavior of dynamical systems. In this article we intertwine the subdivision process with the computation of invariant measures and propose an adaptive scheme for the box refinement which is based on the combination of these methods. Using this new algorithm the numerical effort for the computation of box coverings is in general significantly reduced, and we illustrate this fact by several numerical examples.

1 Introduction

Recently subdivision methods have been successfully applied to the numerical analysis of complex dynamical behavior (e.g. [7,3,2,4,5]). These methods can be used for two essentially different purposes: the first is to understand the geometric structure of an underlying attractor. Secondly the goal may be to approximate the observable dynamical behavior of the underlying system in a specific region of state space by the computation of invariant measures. This paper concerns the second possibility, and we propose an adaptive scheme incorporated into the subdivision technique which allows to reduce the numerical effort significantly.

Obviously the dynamical behavior just needs to be approximated on the support of a certain invariant measure. Indeed, the idea for the adaptive principle stated here is to intertwine the subdivision techniques with the computation of a natural invariant measure, an SBRmeasure, say. Roughly speaking, the size of the covering boxes is reduced in those parts of state space where the natural invariant measure μ is concentrated, and, on the other hand, boxes are not subdivided in areas which have μ-measure zero.

The main goal of this article is to illustrate the efficiency of the new method by numerical examples. For that purpose we consider several dynamical systems for which the SBR-measure is known analytically since this allows us to compare the numerical results obtained by the adaptive subdivision algorithm to those obtained by the standard subdivision procedure. The adaptive algorithm is essentially based on the combination of two existing methods for which convergence results are known. However, this fact does not immediately imply convergence of the adaptive method as well. Rather this theoretical but relevant problem is currently under investigation, and the results will be published elsewhere.

An outline of the paper is as follows: in Sect. 2 we recall the standard subdivision technique from [3]. The numerical method for the approximation of SBR-measures is described in Sect. 3. Then, in Sect. 4, we present our adaptive subdivision technique, and the efficiency of this method is illustrated by several examples in Sect. 5.

2 The standard subdivision algorithm

The purpose is to approximate invariant sets of discrete dynamical systems of the form

x_{j+1} = f(x_{j}), j = 0, 1, …,

where f is a continuous mapping on ℝ^{n}. The central object which is approximated by the subdivision algorithm developed in [3] is the so-called relative global attractor,

A_{Q} = ⋂_{j≥0} f^{j}(Q), (1)

where Q ⊂ ℝ^{n} is a compact subset. Roughly speaking, the set A_{Q} should be viewed as the union of unstable manifolds of invariant objects inside Q. In particular, A_{Q} may contain subsets of Q which cannot be approximated by direct simulation.

The subdivision algorithm for the approximation of A_{Q} generates a sequence B_{0}, B_{1}, B_{2}, … of finite collections of boxes with the property that for all integers k the set Q_{k} = ⋃_{B∈B_{k}} B is a covering of the relative global attractor under consideration. Moreover the sequence of coverings is constructed in such a way that the diameter of the boxes,

⋆ Research of the authors is partly supported by the Deutsche Forschungsgemeinschaft under Grant De 448/5-2 and by the Konrad-Zuse-Zentrum für Informationstechnik Berlin.

diam(B_{k}) = max_{B∈B_{k}} diam(B)

converges to zero for k → ∞.

Given an initial collection B_{0}, one inductively obtains B_{k} from B_{k−1} for k = 1, 2, … in two steps.

(i) Subdivision: Construct a new collection B̂_{k} such that

⋃_{B∈B̂_{k}} B = ⋃_{B∈B_{k−1}} B

and diam(B̂_{k}) ≤ θ diam(B_{k−1})

for some 0 < θ < 1. (ii) Selection: Define the new collection B_{k} by

B_{k} = {B ∈ B̂_{k}:f^{−1}(B)∩B̂ ≠∅ for some B̂ ∈B̂_{k}} .

The following proposition establishes a general convergence property of this algorithm.

Proposition 2.1 ([3]). Let A_{Q} be the global attractor relative to the compact set Q, and let B_{0} be a finite collection of closed subsets with Q_{0} = Q. Then

lim_{k→∞} h (A_{Q}, Q_{k}) = 0,

where we denote by h(B, C) the usual Hausdorff distance between two compact subsets B, C ⊂ ℝ^{n}.

3 Approximation of SBR-measures

Recently it has been shown in [4] how to compute numerically approximations of an SBR-measure supported on a hyperbolic invariant set. Since we want to use this method in our adaptive scheme we now sketch the main ingredients of this algorithm. To make the ideas more transparent we simplify the description drastically by avoiding all technical details concerning the underlying mathematical foundation in Ergodic Theory.

The crucial observation is that the calculation of invariant measures can be viewed as a fixed point problem. Let M be the set of probability measures on ℝ^{n}. Then μ ∈ M is invariant if and only if it is a fixed point of the Frobenius-Perron operator P : M → M,

(P μ)(B) = μ(f^{−1}(B)) for all measurable B ⊂ ℝ^{n}. (1)

For a discretization of the operator P : M → M we replace M by a finite dimensional set M_{k}: let B_{i} ∈ B_{k}, i = 1, …, N, denote the boxes in the covering obtained after k steps in the subdivision algorithm and set as before Q_{k} = ⋃_{B∈B_{k}} B. We choose M_{k} to be the set of “discrete probability measures” on B_{k}, that is,

M_{k} = {u : B_{k} → [0, 1] | ∑_{i=1}^{N} u(B_{i}) = 1} .

Assuming that f(Q_{k}) ⊂ Q_{k} the discretized Frobenius-Perron operator P_{k} : M_{k} → M_{k} is given by

v = P_{k}u, v(B_{i}) = ∑_{j=1}^{N} \frac{m(f^{−1}(B_{i}) ∩ B_{j})}{m(B_{j})}u(B_{j}),  (2)
i = 1, …, N,

where m denotes Lebesgue measure. Now a fixed point u = P_{k}u of P_{k} provides an approximation to an invariant measure of f. Remark 3.1. For the mathematically precise statement on the convergence of this method one would have to introduce the concept of small random perturbations. The reason is that this allows one to use a result of Yu. Kifer on the convergence of invariant measures in the perturbed systems to the SBR-measure ([9]). However, the purpose of this article is to develop and to test an adaptive scheme for the box refinement in the subdivision algorithm rather than to explain the theoretical background concerning the computation of SBR-measures. Therefore the reader is referred to [4] for the rigorous mathematical treatment.

4 The adaptive subdivision algorithm

As mentioned above, the standard subdivision algorithm may approximate a part of the global attractor which is dynamically irrelevant in the sense that no invariant measure has support on this subset. The reason is that each box is subdivided in a step of the subdivision algorithm regardless of any information on the dynamical behavior. In particular, also those subsets of the relative global attractor corresponding to unstable or transient dynamical behavior are approximated by the standard procedure.

On the other hand, if one is mainly interested in the approximation of the support of the (natural) invariant measure rather than in the precise geometric structure of the global attractor then this strategy may lead to unnecessary high storage and computation requirements. In the following we present a modified subdivision strategy which avoids this drawback: roughly speaking,

– in the subdivision step we use the information on the actual approximation of the invariant measure to decide whether or not a box should be subdivided;

– in the selection step we keep only those boxes which have a nonempty intersection with the support of the invariant measure. To be more precise, let {δ_{k}} be a sequence of positive

real numbers such that δ_{k} → 0 for k → ∞. The algorithm generates a sequence of pairs (B_{0}, u_{0}), (B_{1}, u_{1}), (B_{2}, u_{2}), … where the B_{k}’s are finite collections of compact subsets of ℝ^{n} and the discrete measures u_{k} : B_{k} → [0, 1] can be interpreted as approximations to the SBR-measure μ_{SBR}: u_{k}(B) ≈ μ_{SBR}(B) for all B ∈ B_{k}.

Given an initial pair (B_{0}, u_{0}), one inductively obtains (B_{k}, u_{k}) from (B_{k−1}, u_{k−1}) for k = 1, 2, … in three steps:

(i) Subdivision: Define

B_{k−1}^{−} = {B ∈ B_{k−1} : u_{k−1}(B) < δ_{k−1}}

and

B_{k−1}^{+} = B_{k−1}\B_{k−1}^{−}.

Construct a new (sub-)collection B̂_{k}^{+} such that

⋃_{B∈B̂_{k}^{+}} B = ⋃_{B∈B_{k−1}^{+}} B

where diam(B̂_{k}^{+}) ≤ θ diam(B_{k−1}^{+}) for some 0 < θ < 1.

(ii) Calculation of the invariant measure: Set

B̂_{k} = B_{k−1}^{−} ∪ B̂_{k}^{+}.

For the collection B̂_{k} calculate the approximating invariant measure as the fixed point û_{k} of the discretized Frobenius-Perron operator defined by (2).

(iii) Selection: Set

B_{k} = {B ∈ B̂_{k} : û_{k}(B) > 0} and u_{k} = û_{k}|_{B_{k}}.

Remark 4.1. (a) In the realization of the algorithm we typically subdivide the boxes in the collection B_{k}^{+} by bisection. This guarantees that the number of boxes is not growing too fast. For the details concerning the implementation the reader is again referred to [3,4].

(b) In principle there is some freedom in choosing the sequence {δ_{k}} of positive numbers used in the subdivision step. Note however that this sequence determines the number of boxes which will be subdivided and hence it has a significant influence on the storage requirement. In the computations it turned out to be quite efficient to choose the average

δ_{k} = \frac{1}{N_{k}} ∑_{B∈B_{k}} u_{k}(B) = \frac{1}{N_{k}} ,

where N_{k} is the number of boxes in B_{k}.

(c) In the numerical realization of the selection step (iii) we check whether û_{k}(B) > ϵ where ϵ > 0 is chosen sufficiently small with respect to the machine precision.

5 Numerical examples

In this section we illustrate the efficiency of the adaptive scheme by several numerical examples. First we consider three one-dimensional mappings for which the SBRmeasures are known analytically. For these cases we will see that, as expected, the new technique is particularly useful if the underlying invariant density has singularities. Additionally we consider the Hénon map as a twodimensional example and show the box refinement produced by the adaptive subdivision algorithm at a certain step.

Before proceeding let us indicate some details concerning the implementation of the adaptive subdivision algorithm:

Table 1. Comparison between the standard and the adaptive subdivision algorithm for the Logistic Map. The minimal box volume in each row is 2^{−ℓ}

ℓ number of boxes L^{1}-error
standard adaptive standard adaptive
6 64  17  0.1390  0.1431
8 256  34  0.0670  0.0765
12 4096  151  0.0210  0.0258
16 65536  679  0.0064  0.0073
20 2^{20}  2831  -  0.0021

(a) The subdivision is always done by bisection and the boxes are stored in a binary tree. This way we keep the storage requirement at a low level.

(b) For the computation of the transition probabilities m(f^{−1}(B_{i}) ∩ B_{j}) in (2) we use an exhaustion technique as described in [8]. This method is particularly useful when – as in our examples – local Lipschitz constants are available for the underlying dynamical system.

(c) The computation of the discrete measures is done by an inverse power method. In the solution of the corresponding linear systems the fact is taken into account that the discretized Frobenius-Perron operator is extremely sparse.

The adaptive algorithm is integrated into the C++ code GAIO (Global Analysis of Invariant Objects). A link to a detailed description of GAIO can be found on the homepages of the authors.

Three one-dimensional examples

Motivated by the numerical investigations in [6] we apply the adaptive subdivision algorithm to three different one-dimensional dynamical systems on the interval [0, 1]. In each case we have chosen the initial collection B_{0} = {[0, 1]}.

1. As a first example we consider the Logistic Map f_{1} : [0, 1] → [0, 1], f_{1}(x) = λx(1 − x) for λ = 4. The unique absolutely continuous invariant measure μ of f_{1} has the density

h_{1}(x) = \frac{1}{π\sqrt{x(1 − x)}}

(see e.g. [10]). Using the standard and the adaptive algorithm we have approximated this density on several levels and the results are shown in Table 1. We remark that even the computation for ℓ = 20 only takes about 50 sec on an MIPS R4400 cpu. In Fig. 5 we illustrate the fact that the size of the boxes is much smaller for those which are close to zero or one. Indeed, this is what we expect since the density has singularities in these points.

2. We consider the map f_{2} : [0, 1] → [0, 1],

f_{2}(x) = \{ \begin{matrix}\frac{2x}{1 − x^{2}} & for 0 ≤ x < \sqrt{2} − 1, \\ \frac{1 − x^{2}}{2x} & for \sqrt{2} − 1 ≤ x ≤ 1. \end{matrix}

Fig. 1. Illustration of the relation between the density and the actual box refinement produced by the adaptive subdivision algorithm for ℓ = 10: a the density h_{1}; b the radii versus the midpoints of boxes

Table 2. Comparison of the numerical results for f_{2} (ℓ as in Table 1)

ℓ  number of boxes  L^{1}-error
standard adaptive standard adaptive
6 64  45  0.0027  0.0053
8 256  187  6.5 · 10^{−4} 0.0013
10 1024  759  1.7 · 10^{−4} 3.2 · 10^{−4}
12 4096  3047  4.6 · 10^{−5} 7.8 · 10^{−5}

Its invariant density is

h_{2}(x) = \frac{4}{π(1 + x^{2})} ,

see Fig. 5. In Table 2 we present the numerical results for this case. As expected the application of the adaptive subdivision algorithm is not more efficient than the standard one since the invariant measure is quite close to Lebesgue measure.

3. Finally we consider the map f_{3} : [0, 1] → [0, 1],

f_{3}(x) = (\frac{1}{8} − 2 |x − \frac{1}{2}|^{3})^{\frac{1}{3}} + \frac{1}{2},

with the invariant density

h_{3}(x) = 12 (x − \frac{1}{2})^{2} .

The graphs of f_{3} and h_{3} are shown in Fig. 5. We now discuss the numerical results presented in Table 3. Note that the derivative of f_{3} has singularities at two points inside [0, 1]. This is the reason why several boxes get lost in the realization of the selection step in the standard subdivision algorithm. Consequently the computation of the invariant measure does not lead to satisfying results. In contrast to this no boxes are lost in the application of the adaptive subdivision technique, and accordingly the L^{1}-error is decreasing with an increasing number of subdivision steps.

Table 3. Comparison of the numerical results for f_{3} (ℓ as in Table 1)

ℓ  number of boxes  L^{1}-error
standard adaptive standard adaptive
6 63  11  0.0260  0.2931
8 249  30  0.0105  0.2583
10 993  189  0.0110  0.0435
12 3967  816  0.0133  0.0065

The Hénon map

We apply the adaptive subdivision algorithm to a twodimensional example, namely a scaled version of the well known Hénon map

f : ℝ^{2} → ℝ^{2}, f(x, y) = (1 − ax^{2} + y/5, 5bx).

In the computations we have fixed the parameters by a = 1.2, b = 0.2, and we have chosen B_{0} = {[−2, 2]^{2}}.

In Fig. 4 we present a tiling of the square [−2, 2]^{2} obtained by the adaptive subdivision algorithm after several subdivision steps. The resulting box-collection B consists of the grey boxes shown in part (a) of this figure. We expect that due to the numerical approximation some boxes have positive discrete measure although they do not intersect the support of the real natural invariant measure. Having this in mind we neglect those boxes with very small discrete measure and show in Fig. 4b a subcollection B̃ ⊂ B with the property that

∑_{B∈B̃} u(B) ≈ 0.99 (3)

(see also Remark 4.1(c)).

Remark 5.1. For our choice of the parameter values we cannot explicitly write down a natural invariant measure. Hence it is impossible to compare our numerical results using analytical ones. Moreover, it is not even known for an arbitrary choice of the parameter values whether or not the Hénon map possesses an SBR-measure. However, recently it was proved by M. Benedicks and L.-S. Young that the Hénon map indeed has an SBRmeasure for a “large” set of parameter values, see [1].

Finally we apply the numerical techniques described in [4] to determine the essential dynamical behavior of the Hénon map for this choice of parameter values. An

Fig. 2. a The map f_{2}; and b its invariant density h_{2}

Fig. 3. a The map f_{3}; and b its invariant density h_{3}

Fig. 4. a A tiling of the square [−2, 2]^{2} obtained by the adaptive subdivision algorithm; and b the subcollection B̃ of boxes with discrete density bigger than 0.35 (see also (3))

approximation of the (natural) invariant measure obtained by the adaptive subdivision algorithm is shown in Fig. 5. This computation is based on the total number of 1514 boxes inside the square [−2, 2]^{2}, whereas the support of the invariant measure is covered by 1442 boxes. The results indicate that the Hénon map exhibits complicated dynamical behavior.

Moreover it can be shown that the areas where the density is colored red resp. blue are permuted cyclically by the mapping. Hence altogether we may conclude that for these parameter values the Hénon map exhibits a two-cycle (the “macro-dynamics”) in addition to unpredictable (chaotic) behavior. This fact is also demonstrated by a Java-animation for which a link can be found on the homepages of the authors.

References

1. M. Benedicks, L.-S. Young. Sinai-Bowen-Ruelle measures for certain Henon maps. Invent. math., 112:541–576, 1993

2. M. Dellnitz, A. Hohmann. The computation of unstable manifolds using subdivision and continuation. In H.W. Broer, S.A. van Gils, I. Hoveijn, F. Takens (eds), Nonlinear Dynamical Systems and Chaos, pages 449–459. Birkhäuser, PNLDE 19, 1996

3. M. Dellnitz, A. Hohmann. A subdivision algorithm for the computation of unstable manifolds and global attractors. Numerische Mathematik, 75:293–317, 1997

4. M. Dellnitz, O. Junge. On the approximation of complicated dynamical behavior. To appear in SIAM J. on Numerical Analysis, 1998

5. M. Dellnitz, O. Junge. Almost invariant sets in Chua’s circuit. To appear in Int. J. Bif. and Chaos, 1997 Fig. 5. Illustration of the (natural) invariant measure for the Hénon map. The picture shows the density of the discrete measure on B̃, see (3)

6. J. Ding, Q. Du,, T. Y. Li. High order approximation of the Frobenius-Perron operator. Appl. Math. Comp., 53:151–171, 1993

7. M. Eidenschink. Exploring Global Dynamics: A Numerical Algorithm Based on the Conley Index Theory. PhD Thesis, Georgia Institute of Technology, 1995 8. R. Guder, M. Dellnitz,, E. Kreuzer. An adaptive method

for the approximation of the generalized cell mapping. Chaos, Solitons and Fractals, 8(4):525–534, 1997

9. Yu. Kifer. General random perturbations of hyperbolic and expanding transformations. J. Analyse Math., 47:111–150, 1986

10. A. Lasota, M.C. Mackey. Chaos, Fractals and Noise. Springer, 1994