SET ARITHMETIC AND THE ENCLOSING PROBLEM IN DYNAMICS

MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

Abstract. We study the enclosing problem for discrete and continuous dynamical systems in the context of computer assisted proofs. We review and compare the existing methods and emphasise the importance of developing a suitable set arithmetic for efficient algorithms solving the enclosing problem.

1. Introduction.

In the paper we discuss the following enclosing problem, which in general may be formulated as follows. Given a set S ⊂ ℝ^{n}, a natural number d and a function f : ℝ^{n} → ℝ^{n} provide an algorithm constructing a set S̃ such that

1. it may be proved that f^{d}(S) ⊂ S̃

2. the difference S̃ \ f^{d}(S) is possibly small

3. the algorithm is possibly fast. To ensure rigor of numerical computations needed to fulfill postulate 1 interval arithmetic (see Sec. 5) is used in computations. Note that in general it is difficult to fulfill postulates 2 and 3 simultaneously. Hence, depending on applications, the emphasis is either on postulate 2 or 3. Also note that the natural strategy in solving the problem for d > 1 is to iterate a well chosen algorithm for d = 1. Unfortunately, this combined with interval arithmetic leads to the so called wrapping effect (see Sec. 8), which causes that the algorithms are neither fast nor effective.

The enclosing problem is of interest in itself and several classical approaches to it may be found in literature (see [3],[6],[13]). In these

Date: February 14, 2000.

1991 Mathematics Subject Classification. 65G40,65L70.

Key words and phrases. dynamical systems, enclosing problem, interval arithmetic, rigorous numerical analysis.

Research supported by KBN, Grant 2 P03A 029 12.

Research supported by KBN, Grant 2 P03A 029 12 and NSF–NATO grant DGE– 98–04459.

2 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

papers authors propose to fight the wrapping effect by replacing Cartesian products of intervals used to approximate sets by some finer sets. Lohner [6] uses parallelepipeds and Neumaier [13] uses ellipsoids. The right choice of a countable class of sets used to approximate arbitrary sets seems to be the key strategy in developing efficient solutions to the enclosing problem. In other words, a suitable set arithmetic on top of the elementary interval arithmetic must be developed.

The classical approaches considered the case when S was a very small ball or even a point and the emphasis was on minimizing S̃ and not on the speed of the algorithm. Recently the enclosing problem became of vital importance in computer assisted proofs in dynamics based on topological techniques (see [4], [7], [8], [9], [11], [12], [16], [17]). In such applications the map f is usually a generator of a discrete dynamical system, a Poincaré map of a flow or just a t-translation along the trajectories of an ordinary differential equation. In order to verify the assumptions of certain topological criteria for particular dynamical behaviour like the existence of periodic solutions or chaotic invariant sets it is necessary to verify some inclusions of the form f(S) ⊂ T, where S and T are given subsets of ℝ^{n}. However, unlike the classical case the set S is usually large and often irregular in shape. Thus, in order to verify the inclusion one covers the set S with a large number of small, regular in shape sets S_{i} and tries to verify that f(S_{i}) ⊂ T for every i. This reduces the problem to the classical case with the difference that since the set T is given and large when compared to the sets S_{i}, the emphasis is more on the speed of the algorithm than on minimizing the sets S̃_{i} enclosing the sets f(S_{i}).

The aim of this paper is to review and compare various techniques used to solve the enclosing problem in the context of their applicability to computer assisted proofs in dynamics and provide a general framework for future investigations in this area. In particular we want to emphasise the importance of developing a suitable set arithmetic for efficient algorithms solving the enclosing problem.

Special emphasis is given to the method of Lohner [6]. We reformulate his ideas in a language which is more universal and hopefully simpler to understand for people who do not work in numerical analysis. As a simple by-product we obtain a generalization of the Lohner method to discrete dynamical systems.

2. Preliminaries.

Troughout the paper we make a general assumption, that f may be approximated by a rational map, i.e. there exists a rational map g : ℝ^{n} → ℝ^{n} and a small continuous map w : ℝ^{n} → ℝ^{n} of known bound such that in a certain domain D ⊂ ℝ^{n} containing S we have the following decomposition

f(x) = g(x) + w(x) for x ∈ D.  (1)

In some applications f may be just rational and that case is also of our interest but in most applications f is not rational.

We want to emphasize that the case when f comes from a differential equation (t-translation, Poincaré map) may be treated just as a special case of our problem. Let

x′ = V (x), x ∈ ℝ^{n}  (2)

be an ordinary differential equation in ℝ^{n} generating a flow φ : ℝ^{n} × ℝ → ℝ^{n} and let

φ_{t} : ℝ^{n} ∋ x → φ(x, t) ∈ ℝ^{n}

denote the corresponding t-translation operator for t > 0. Taking a small h > 0, a natural number n satisfying t/n < h and a numerical method Φ : ℝ^{n} × ℝ → ℝ^{n} such that

φ(x, s) = Φ(x, s) + ρ(x, s), s ∈ [0, h]

for a small remainder function ρ we reduce the problem of enclosing φ_{t}(S) to the problem of enclosing f^{n}(S), where f := φ_{t/n} = Φ_{t/n} + ρ_{t/n} fulfills the requirement (1). The case of a Poincaré map is similar in spirit though slightly more complicated, because an extra care is needed to estimate the intersection with the Poincaré section.

If X, Y, Z are sets then the function ω : X × Y → Z is called an operation. If (x, y) ∈ dom ω, then we say that xωy is defined and we write xωy := ω(x, y). The three basic arithmetic operations +, −, · are examples of operations in ℝ.

Throughout the paper ℝ and ℤ denote respectively the set of real numbers and the set of integers. We also need the extended set of real numbers

ℝ̄ := ℝ ∪ { ∞, −∞, NAN },

which consists of all real numbers, plus and minus infinity and a special symbol NAN. Here NAN stands for not-a-number. It is assumed that the result of any division by zero, any undefined operation involving infinity as well as any arithmetic operation involving NAN is NAN. This convention enables us to treat all four elementary arithmetic operations +, −, ×, / as defined on ℝ̄ × ℝ̄.

We treat ℝ̄ as a partially ordered set with ℝ ∪ { ∞, −∞ } ordered in the standard way and NAN non-comparable with any other element.

4 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

3. Arithmetic operations on sets.

In the sequel we will frequently deal with sets enclosing numbers, vectors or matrices. We will use bold characters to denote such sets. We will extensively use the standard notation

f(a) := {f(a) | a ∈ a},

where f : X → Y is a map and a ⊂ dom f ⊂ X. In particular we will use the notation in case of operations. Let ⋄ : X × Y → Z. For a × b ⊂ dom ⋄ ⊂ X × Y we write

a ⋄ b := {a ⋄ b | a ∈ a, b ∈ b}

and call this operation an extension of ⋄ to sets.

The notation is convenient in formal calculations. In particular it is convenient to remember the following straightforward propositions.

Proposition 3.1. If ⋄ : X × X → X is commutative, then its extension to sets is also commutative. If ⋄ : X ×X → X is associative, then its extension to sets is also associative.

Not every property of operations carries over automatically to its extension to sets as the following straightforward proposition shows.

Proposition 3.2. Assume + : X × X → X and · : Y × X → X are operations such that for all y ∈ Y and for all x_{1}, x_{2} ∈ X

y ·(x_{1} + x_{2}) = y · x_{1} + y · x_{2},

whenever both sides make sense. Then for all b ⊂ Y and a_{1}, a_{2} ⊂ X

b ·(a_{1} + a_{2}) ⊂ b · a_{1} + b · a_{2},

whenever both sides make sense. Moreover, the inclusion cannot be replaced by an equality in general.

4. Finite encoding.

We base our considerations concerning computations on the commonly accepted Turing model. There are three important assumptions in this model: computations are performed on a finite sequence of two symbols (0 and 1), there is only a finite number of rules which may be applied to transform the sequence and each rule can read or write only one symbol in the sequence at a time. Obviously, if we want to perform computations on sequences of more than two symbols but the number of symbols is still finite, then a suitable encoding of symbols with sequences of 0 and 1 reduces the computations to the Turing model. A similar relaxation is possible if one wants to let each rule read and write a finite number of symbols at a time.

To perform computations on a finite family of objects which formally are not finite sequences of symbols one has to bijectively encode the objects with such sequences, then perform the computations on the encoding and finally translate the result back to the original objects. For example if one wants to perform computations on the family of subsets of a finite collection of planes in ℝ^{3}, one can first enumerate the planes in some way and then represent a subset with a sequence which has 1 on the i-th place if the i-th plane is present in the subset and 0 otherwise.

A problem arises when we want to perform computations on a set of objects O which is uncountable, because then there is obviously no bijective encoding. In such a situation we usually select a countable subset Ô ⊂ O for which we establish a bijective encoding (we call the elements of Ô representable objects) and we represent arbitrary objects in O via a surjective map ˆ : O → Ô. A drawback of this approach is that we introduce approximations into the computations (an object o ∈ O is approximated by its representation ô ∈ Ô). Since the approximations may take place a large number of times in course of the computation, not only is the final result not exact but usually we do not know how bad the approximation is. Nevertheless we will see in the sequel that sometimes the approach provides exact bounds instead of exact solutions.

As an example consider the common solution used in hardware to perform computations on real numbers. One first chooses a finite subset ℝ̂ ⊂ ℝ̄ such that

0, −∞, +∞, NAN ∈ ℝ̂,  (3)
x ∈ ℝ̂ ⇒ −x ∈ ℝ̂  (4)

and treats it as the set of symbols. Usually ℝ̂ is taken as a set of certain binary fractions to facilitate further encoding into sequences of 0 an 1. Then one takes a map ⟨·⟩ : ℝ̄ → ℝ̂ such that

x ∈ ℝ̂ ⇒ x = ⟨x⟩,  (5)
x ≤ y ⇒ ⟨x⟩ ≤ ⟨y⟩.  (6)

Such a map is obviously a surjection. Any real number x which is not in ℝ̂ is replaced in the course of computations by ⟨x⟩. The process is known as rounding. We call the number ⟨x⟩ the representation of x and the pair (ℝ̂, ⟨·⟩) the representation of ℝ̄. One can easily verify that for any x ∈ ℝ̄ \ {NAN}

↓(x) ≤ ⟨x⟩ ≤ ↑(x),

where

↓(x) := max { y ∈ ℝ̂ | y ≤ x }
↑(x) := min { y ∈ ℝ̂ | y ≥ x }

are the lower and upper representations of x.

As we already mentioned, the price one pays for such an approach is that the computations are approximate and there is no way (other than statistical analysis) to tell how bad the approximations is.

Even if the input of the computation consists of representable numbers, the result of an arithmetic operation on representable numbers may not be a representable number and has to be replaced by its representation. This introduces deviations from the exact results which keep propagating in course of computations. Though usually there is no way to obtain exact results, it is possible to get exact lower and upper bounds. This is achieved by means of the interval arithmetic.

5. Interval arithmetic.

The concept of interval arithmetic goes back to 50’s and has several fathers. Probably the earliest paper on this subject is the work of M. Warmus [14]. Several other early papers are listed in the monograph by R. Moore [10].

The idea of the interval arithmetic consists in replacing every representation ⟨x⟩ on input by the interval [↓(x), ↑(x)] and perform calculations on intervals in such a way that at every stage of calculations the interval contains the exact result.

To be more precise we define the set of intervals over A ⊂ ℝ̄ by

𝕀(A) := {[a^{−}, a^{+}] | a^{−}, a^{+} ∈ A, a^{−} ≤ a^{+}}.

In particular we will write 𝕀, 𝕀̄, 𝕀̂ respectively for 𝕀(ℝ), 𝕀(ℝ̄) and 𝕀(ℝ̂). The elements of 𝕀̂ will be called representable intervals. Note that the only interval in 𝕀̄ and 𝕀̂ containing NAN is [NAN, NAN].

The set of real numbers may be viewed as a subset of 𝕀(ℝ̄) via the embedding

ℝ̄ ∋ x → [x, x] ∈ 𝕀(ℝ̄).

The size of an interval a = [a^{−}, a^{+}] ∈ 𝕀 is defined as the difference a^{+} − a^{−}. We also define the mid point of a by

mid(a) := ⟨\frac{a^{−}+a^{+}}{2} ⟩.

The reader should notice the by the above definition the mid point of an interval is always a representable number.

To perform calculations on intervals we use representable intervals. To encode representable intervals we just store two endpoints of the interval. If a ∈ 𝕀 is an arbitrary interval then we represent a as the smallest representable interval containing a. Such an interval will be called the representation of a and denoted by ↕ (a). One can easily verify that for any a^{−}, a^{+} ∈ ℝ, a^{−} ≤ a^{+}

↕([a^{−}, a^{+}]) = [↓(a^{−}), ↑(a^{+})].

Obviously such an encoding is only surjective.

As in the case of the arithmetic of representable numbers, arithmetic operations on intervals cannot be performed exactly. However in the case of intervals there is an easy way to guarantee that the computed result encloses the exact result. To achieve this, given an arithmetic operation ⋄, one defines the operation ˆ⋄ by

aˆ⋄b := \begin{cases}↕(a ⋄ b) & if NAN ∉ a ⋄ b \\ [NAN, NAN] & otherwise . \end{cases}

for a, b ∈ 𝕀̄. An important remark is that the right hand side of the above formula is useless from the computational point of view, since we cannot compute a ⋄ b (comp. last paragraph of Section 4.). Fortunately, for each concrete arithmetic operation ⋄ ∈ {+, −, ×, /} one can easily find a simple algorithm transforming the encoding (endpoints) of a, b to the encoding of the result, which makes the computation of aˆ⋄b possible.

Sometimes it is convenient to decompose an interval as a sum of a representable number and a representable interval containing zero. For this end we introduce the notation

𝕀̂_{0} := { a ∈ 𝕀̂ | 0 ∈ a }.

and define the decomposition of an interval a ∈ 𝕀 as the pair

(mid(a), aˆ-mid(a)) ∈ ℝ̂ × 𝕀̂_{0}.

Note the following straightforward proposition.

Proposition 5.1. If (c, c_{0}) ∈ ℝ̂ × 𝕀̂_{0} is the decomposition of c then c ⊂ c + c_{0}.

The elements of 𝕀^{n} will be called n-dimensional interval vectors and elements of 𝕀^{n×n} - n-dimensional interval matrices.

The encoding of interval vectors and matrices is componentwise. We just represent interval vectors and interval matrices as sequences of representations of their entries.

8 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

We define the size of an interval vector or an interval matrix as the maximum of sizes of the components, the mid point of an interval vector or an interval matrix as the interval vector or matrix consisting of the mid points of the components. Similarly we define the decomposition of an interval vector or matrix.

Apart from interval vectors and matrices we also define interval sets. An interval set is a Cartesian product of intervals. An interval representable set is a Cartesian product of representable intervals. There is a one-to-one correspondence between interval vectors and interval sets given by

(a_{1}, a_{2}, …, a_{n}) ↔ a_{1} × a_{2} × … × a_{n}

for a_{1}, a_{2}, … a_{n} ∈ 𝕀. Since we can treat matrices as n × m-vectors, a similar correspondence applies to matrices. In the sequel we will freely pass from one interpretation to the other.

If s ⊂ ℝ^{n} is a set, then [s] will stand for the minimal representable interval set which contains a.

We define the sum of interval vectors or interval matrices, the product of an interval and an interval vector or an interval matrix, the product of interval matrices, the scalar product of interval vectors and the product of an interval matrix by an interval vector exactly as in the case of real vectors and matrices.

6. Rational expressions and rational functions.

Practically speaking the numeric evaluation of a function at a given argument is somehow reduced to calculations with the four elementary arithmetic operations +, −, ·, /. Therefore what is actually computed is a rational function. Since we will evaluate rational functions with interval arguments where the arithmetic operations are neither associative nor distributive, we need a slightly modified definition of a rational function.

We first recall the formalism of arithmetic expressions. For this purpose we fix a countably infinite set of symbols V = { X_{1}, X_{2}, … }, called variables, and define a word to be an arbitrary finite sequence of symbols in V ∪ ℝ̂ ∪ { ( , ) , +, −, ·, / }. If ω is a word then |ω| will stand for the length of ω, i.e. the number of elements in ω. The set of (arithmetic) expressions W is defined as the smallest set of words satisfying the following two properties

x ∈ V ∪ ℝ̂ ⇒ x ∈ W

ω_{1}, ω_{2} ∈ W ⇒ (ω_{1} + ω_{2}), (ω_{1} − ω_{2}), (ω_{1} · ω_{2}), (ω_{1}/ω_{2}) ∈ W The following proposition is an easy exercise.

Proposition 6.1. If ω ∈ W then exactly one of the following possibilities holds

(i) ω = X_{i} for some X_{i} ∈ V (ii) ω = s for some s ∈ ℝ̂ (iii) there exist ω_{1}, ω_{2} ∈ W such that ω = (ω_{1} + ω_{2}) or ω =

(ω_{1} − ω_{2}) or ω = (ω_{1} · ω_{2}) or ω = (ω_{1}/ω_{2}).

For instance ((X_{1} +X_{2})+X_{3}) is an arithmetic expression but +++, +(X_{1}X_{2}), X_{1} + X_{2} + X_{3}, X_{1} + (X_{2} + X_{3}) are not.

We define recursively Var(ω), the set of variables in ω by

Var(ω) := \begin{cases}∅ & if ω = s, where s ∈ ℝ̂ \\ \{X_{i}\} & if ω = X_{i}, where X_{i} ∈ V \\ Var(ω_{1}) ∪ Var(ω_{2}) & if ω = (ω_{1} ⋄ ω_{2}) \end{cases} (7)

where ⋄ ∈ { +, −, ·, / }.

Let ω be an expression. Let n ∈ ℕ be large enough so that Var(ω) ⊂ {X_{1}, … X_{n}}. Let x_{1}, …, x_{n} be real numbers. We define recursively ω(x_{1}, … x_{n}), the evaluation of ω at (x_{1}, …, x_{n}) ∈ ℝ^{n} by

ω(x_{1}, … x_{n}) :=

\begin{cases}s & if ω = s, where s ∈ ℝ̂ \\ x_{i} & if ω = X_{i}, where X_{i} ∈ V \\ ω_{1}(x_{1}, … x_{n}) ⋄ ω_{2}(x_{1}, … x_{n}) & if ω = (ω_{1} ⋄ ω_{2}) \end{cases} (8)

where ⋄ ∈ { +, −, ·, / }.

Similarly we define recursively ⟨ω⟩(x_{1}, …, x_{n}), the interval evaluation of ω at (x_{1}, …, x_{n}) ∈ 𝕀̂^{n} by

⟨ω⟩(x_{1}, … x_{n}) :=

\begin{cases}s & if ω = s, where s ∈ ℝ̂ \\ x_{i} & if ω = X_{i}, where X_{i} ∈ V \\ ⟨ω_{1}⟩(x_{1}, … x_{n})⟨⋄⟩⟨ω_{2}⟩(x_{1}, … x_{n}) & if ω = (ω_{1} ⋄ ω_{2}). \end{cases} (9)

Let ω_{1} := ((X_{1} +X_{2})+X_{3}) and ω_{2} := (X_{1} +(X_{2} +X_{3})) be arithmetic expressions. Obviously ω_{1} ≠ ω_{2} (the sequences of symbols are different) but for any (x_{1}, x_{2}, x_{3}) ∈ ℝ^{3}

ω_{1}(x_{1}, x_{2}, x_{3}) = ω_{2}(x_{1}, x_{2}, x_{3}),

because of the associativity of the addition operation in ℝ. However, the interval evaluations ⟨ω_{1}⟩(x_{1}, x_{2}, x_{3}) and ⟨ω_{2}⟩(x_{1}, x_{2}, x_{3}) need not be the same, because ⟨+⟩ is not associative in general. This explains why we set up the definition of the arithmetic expression in such a way that there is no ambiguity in the order of evaluation.

10 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

Assume ω = (ω_{i})_{i=1,... ,m} is a vector of arithmetic expressions such that Var(ω_{i}) ⊂ {X_{1}, … X_{n}}. There is then an associated function f^{ω} = (f_{1}^{ω}, … f_{m}^{ω}) : ℝ^{n} → ℝ^{m} given by

f_{i}^{ω}(x_{1}, … x_{n}) = ω_{i}(x_{1}, … x_{n}) for i = 1, 2, … m

and an associated interval function ⟨f⟩^{ω} = (⟨f⟩_{1}^{ω}, … ⟨f⟩_{m}^{ω}) : 𝕀̂^{n} → 𝕀̂^{m} given by

⟨f⟩_{i}^{ω}(x_{1}, … x_{n}) = ⟨ω_{i}⟩(x_{1}, … x_{n}) for i = 1, 2, … m.

We say that the function g : ℝ^{n} → ℝ^{m} is a rational function if g = f^{ω} for a vector ω of arithmetic expressions. The function ⟨f⟩^{ω} is then called an interval extension of g. Note that for reasons explained earlier there may be two different arithmetic expressions ω_{1} and ω_{2} such that f^{ω_{1}} = f^{ω_{2}}. Consequently the interval extension of g need not be unique. Nevertheless we will write ⟨g⟩ meaning any possible interval extension of g.

The following theorem, which may be proved easily by induction, is crucial for interval evaluations of rational functions.

Theorem 6.2. Assume f : ℝ^{n} → ℝ^{m} is a rational function. Then for any x_{1}, … x_{n} ∈ 𝕀̂^{n} we have

f(x_{1}, … x_{n}) ⊂ ⟨f⟩(x_{1}, … x_{n}).  (10)

Here there is one important comment in order. The left-hand side of (10) more formally should be written as f(x_{1}×⋯×x_{n}) and denotes the image of the set x_{1} × ⋯ × x_{n} by f. Therefore it is a pure notation and does not contribute towards a method of its computation. The righthand side may be computed easily from its definition thus providing a method of enclosing the left-hand side.

We will say that an operation ⋄ : ℝ^{m} × ℝ^{n} → ℝ^{p} is rational if it is a rational function. Note that matrix product provides a simple example of a rational operation.

We can rewrite the last theorem in terms of an operation as in the following corollary.

Corollary 6.3. If ⋄ : ℝ^{m} × ℝ^{n} → ℝ^{p} is a rational operation and v ∈ ℝ^{m}, w ∈ ℝ^{n} then

v ⋄ w ⊂ v⟨⋄⟩w.

7. Enclosures via simple representations of sets.

The first question we have to ask when discussing the problem of estimation of f^{d}(S) is how we should represent sets. Actually the important question is what the family of representable sets should be, because once the family Ŝ of representable sets is fixed, it is very natural to define Ŝ, the representation of S as the smallest representable set containing S. Since then f^{d}(S) ⊂ f^{d}(Ŝ), the next question we have to address is how to construct an estimate of f^{d}(S) when S is a representable set. Obviously the two questions are strongly interconnected.

The choice of the class of representable sets depends on several factors. The first factor is the dimension of the space. In dimension one intervals serve well the purpose of representing compact, connected sets. However in higher dimensions we should also consider what our sets look like. If they are as simple as cuboids, balls or parallelepipeds the strategy will be simpler then in the case when they are relatively big and irregular in shape.

There are a few classes of simple sets which can be easily parametrized by a few numbers or vectors, which enables easy encoding of such sets, when the parameters are representable.

7.1. Bounds via balls. Assume a norm ||.|| in ℝ^{n} is given. Recall that the ball centered at x ∈ ℝ^{n} of radius r > 0 is the set B(x, r) := {y ∈ ℝ^{n} | ||x − y|| ≤ r}. We say that the ball B(x, r) is representable if both x and r are representable. We encode such a ball as the pair (x, r) ∈ ℝ̂^{2}.

The simplest strategy for finding an enclosure for f(S) requires that f is Lipschitz with Lipschitz constant L ∈ ℝ̂, S = B(x, r) for some representable vector x ∈ ℝ^{n} and some representable number r ∈ ℝ^{+}. It also requires that an algorithm is known which for a given representable ϵ > 0 constructs a vector y ∈ B(f(x), ϵ). One then easily verifies the following proposition. Proposition 7.1.

f(B(x, r)) ⊂ B(y, Lˆ·r ˆ+ ϵ).

7.2. Logarithmic norms. The applicability of the last proposition depends crucially on the knowledge of a possibly small Lipschitz constant. In case the map f = φ_{t} is the t-translation of the flow generated by the diffrential equation (2), a Lipschitz constant for f may be obtained from the classical differential inequalities theory as e^{Kt}, where

K := sup { ||d_{x}V || | x ∈ φ(S, [0, t]) }

12 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

and ||.|| denotes a norm in ℝ^{n}. Unfortunately this is usually a very bad Lipschitz constant. In particular it is always greater than one, even if the flow is contracting.

Much better Lipschitz constants for f may be obtained by means of the so called logarithmic norm. If Q is a square matrix, then its logarithmic norm μ(Q) is defined by

μ(Q) := lim_{h→0,h>0} \frac{||I + hQ|| − 1}{h}.

It turns out that the differential inequalities theory has a counterpart in which the norm is replaced by the function μ (for a very detailed case concerning non-autonomous differential equations see Theorem 10.6 in [5]). For our purposes a very special case of that theorem may be formulated as follows.

Theorem 7.2. If

K_{μ} := sup { μ(d_{x}V) | x ∈ φ(S, [0, t]) }

then e^{K_{μ}t} is a Lipschitz constant for f = φ_{t}.

7.3. Bounds via interval sets. If f is a rational function and S is a representable interval set, one can use the interval evaluation of f as an enclosure of f(S):

Proposition 7.3. Assume f is a rational function and S = x ∈ 𝕀̂^{n} is an interval set. Then f(x) ⊂ ⟨f⟩(x).

If f is not rational, we still may be able to approximate f by a rational function g and a small function w from the decomposition (1). We assume that we know an interval set w ∈ 𝕀̂^{n} such that w(x) ∈ w for x ∈ D. We do not assume how big the bound w is, but the bound for f(S) we are looking for will crucially depend on the size of w.

Proposition 7.4.

f(x) ⊂ ⟨g⟩(x) ˆ+ w.

7.4. Set valued mean value theorem. If f is C^{1}, then we define the Jacobian matrix of f at x by

Jxf := ([\frac{∂f_{i}}{∂x_{j}} (x)])_{i,j=1,n}

We will first prove the set valued mean value theorem.

Theorem 7.5. For any h ∈ 𝕀_{0}^{n} we have

f(x + h) ⊂ f(x) + J_{x+h}f · h.

Proof: Take h ∈ h. Using the standard mean value theorem for the components f_{i} of f we get for some θ_{i} ∈ [0, 1] and i = 1, 2, … n

f_{i}(x + h) = f_{i}(x) + ∑_{j=1}^{n} \frac{∂f_{i}}{∂x_{j}} (x + θ_{i}h)h_{j}.

Thus it follows that

f(x + h) ∈ f(x) + (∑_{j=1}^{n}[\frac{∂f_{i}}{∂x_{j}} (x + h)]h_{j})_{i=1,n} ⊂

f(x) + J_{x+h} · h ⊂ f(x) + J_{x+h} · h.

7.5. Improved bounds via interval sets. We again assume that S is an interval set, i.e. S = x for a certain interval vector x. We rewrite x as the sum x + r, where (x, r) is the decomposition of x and take the pair (x, r) as the representation of S.

Lemma 7.6. Assume (x_{1}, x_{1}^{0}) is a decomposition of ⟨g⟩(x) and put

r_{1} := Jxgˆ·r ˆ+ w ˆ+ x_{1}^{0},

then f(x + r) ⊂ x_{1} + r_{1}.

Proof: Let u ∈ x + r. Then u = x + r for some r ∈ r. Thus we have

f(u) = g(u) + w(u) ∈ g(x) + Jxg · r + w(u) ⊂
⟨g⟩(x) + Jxg · r + w ⊂ x_{1} + Jxg · r + w + x_{1}^{0} ⊂
x_{1} + Jxgˆ·r ˆ+ w ˆ+ x_{1}^{0} = x_{1} + r_{1}.

The classes of representable sets we considered in this section as well as classes we will consider in the next section are gathered in the following table. The table contains the set, their geometric form and their encoding.

Basic sets

Description Geometric Encoding Comments

Form

a ball B(c, r) (c, r) c ∈ ℝ̂^{n}, r ∈ ℝ̂,

r > 0

the geometric form depends on norm used

an interval c + r (c, r) c ∈ ℝ̂^{n},
set  r ∈ 𝕀̂_{0}^{n}
a paralle- c + B · r (c, B, r) c ∈ ℝ̂^{n}, B ∈ 𝕀̂^{n×n},
lepiped  r ∈ 𝕀̂_{0}^{n}
a cuboid c + Q · r (c, Q, r) c ∈ ℝ̂^{n}, Q ∈ 𝕀̂^{n×n},

r ∈ 𝕀̂_{0}^{n} and Q contains an orthogonal matrix

a doubleton c + C · r_{0} + s (c, C, r_{0})  c ∈ ℝ̂^{n}, C ∈ ℝ̂^{n×n},
followed by the r_{0} ∈ 𝕀̂_{0}^{n}
encoding of s s - another basic set
an elliptic c + L · B(0, r) (c, L, r)  c ∈ ℝ̂^{n}, L ∈ ℝ̂^{n×n}
set  r ∈ ℝ̂

8. Shape sensitive enclosures.

At first it may seem that the enclosures presented in the previous chapter are perfect if our set S is an interval set or a ball. Actually the choices are usually satisfactory if we need to estimate only the first iterate of f. However if we try to iterate this process, we often encounter problems. The reason is that the image of an interval or a ball need not resemble the shape of the argument. Thus enclosing the image in an interval or a ball we may introduce a strong overestimation which will keep propagating under iteration in an exponential way. This process is known as the wrapping effect.

Lohner (see [6]) proposed a solution to this problem in the case f is a translation map along trajectories of a differential equation. We present a version of his method adapted to the more general discrete case in the Lemmas 8.3, 8.4 and 8.5.

A good way to understand the key idea of Lohner approach is to consider first the case when f is a linear map. Even in this simple case the wrapping effect will take place, because linear maps neither transform interval sets to interval sets nor balls to balls. To avoid wrapping in this case we need a class of representable sets which is closed under linear maps. A possible choice are parallelepipeds in ℝ^{n}. A parallelepiped in ℝ^{n} may be written as p := c + B · r, where c ∈ ℝ^{n}, B ∈ ℝ^{n×n} and r ∈ 𝕀^{n}. We will say that p is representable if c, B and r are representable. The encoding of a representable parallelepiped is obviously (c, B, r).

Using parallelepipeds with linear maps we would compute the images exactly if not the rounding errors. But we can expect a very good upper bound.

8.1. Rearrangement. One of the important problems concerning computations with a given class of representable sets is enclosing sets which appear in course of computations by a representable set. We would like to consider here a relatively simple problem of enclosing with a parallelepiped the sum of a parallelepiped and a small interval set. Thus a parallelepiped c + P · r and an interval v are given and we want to find a possibly small parallelepiped c + P_{1} · r_{1} such that

c + P · r + v ⊂ c + P_{1} · r_{1}.

We will consider a slightly more general problem, where P and P_{1} are replaced by interval matrices P and P_{1}. Obviously it is sufficient to find a matrix P_{1} and an interval vector r_{1} such that

P · r + v ⊂ P_{1} · r_{1}.

The simplest choice is to fix arbitrarily P_{1} and search only for r_{1}. To do so we introduce the concept of a weak inverse of an interval matrix. We say that an interval matrix P^{∗} ∈ 𝕀^{n×n} is a weak inverse of an interval matrix P ∈ 𝕀^{n×n} if there exists a non-singular matrix P ∈ P such that P^{−1} ∈ P^{∗}. Obviously the weak inverse is not unique. Note that to find a weak inverse of a given interval matrix P it is enough to find an enclosure of the inverse matrix of a non-singular selector of P. In particular, if P contains an orthogonal matrix then P^{T} is its weak inverse.

For P, P_{1} ∈ 𝕀^{n×n} and r, v ∈ 𝕀^{n} put

Rearr_{1}(P, P_{1}, r, v) := (P_{1}^{∗} ˆ·P)ˆ·r ˆ+ P_{1}^{∗} ˆ·v  (11)

We have the following lemma.

Lemma 8.1. Assume P_{1} is an interval matrix and P_{1}^{∗} is its week inverse. Then for r_{1} := Rearr_{1}(P, P_{1}, r, v) we have the inclusion

P · r + v ⊂ P_{1} · r_{1}.

16 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

Proof: We have

P · r + v ⊂ P_{1} · P_{1}^{∗} ·(P · r + v) ⊂
P_{1} ·(P_{1}^{∗} · P · r + P_{1}^{∗} · v) ⊂ P_{1} ·(P_{1}^{∗} ˆ·Pˆ·r ˆ+ P_{1}^{∗} ˆ·v) ⊂
P_{1} · r_{1}.

Obviously in practice we try to choose P_{1} in such a way that P_{1} · r_{1} is as small as possible. Note that we put explicit parentheses around the matrix P_{1}^{∗} ˆ·P in the formula (11) requiring that the computation of the product precedes the multiplication by r. The reason is as follows. Matrices in the product P_{1}^{∗} ˆ·P in general are of relatively small size compared to the size of r, therefore the potential wrapping effect in this computation is also relatively small. The biggest wrapping will be introduced in the multiplication by r. Hence it is worth to put maximum effort to minimize this wrapping. There is no wrapping in the multiplication by a matrix, when the off diagonal elements of the matrix are non-negative (see Anguelov, [1],[2] ), in particular when the matrix is the identity matrix. Thus one can expect small wrapping when the matrix is close to identity. When P_{1} is a selector of P, the matrix P_{1}^{∗} ˆ·P contains identity. Thus, since the size of this matrix is usually small, it is close to identity and one can expect small wrapping.

A well known fact in numerical analysis is that inverting a matrix numerically is not a good procedure if the goal is to solve a linear equation. It turns out that finding numerically an inverse of P_{1} may be replaced in the rearrengment problem by solving a suitable interval linear equation. To explain this assume that P ∈ 𝕀^{n×n}, P_{1} ∈ ℝ^{n×n} is an invertible selector of P, R := Pˆ-P_{1} and r, v ∈ 𝕀^{n}. Let t be an interval vector obtained by applying the Gauss elimination to the linear interval equation

P_{1} ˆ·x = Rˆ·r ˆ+ v.

Note that t in particular encloses the set of solutions of all linear equations

P_{1} · x = b, where x ∈ ℝ^{n} and b ∈ Rˆ·r ˆ+ v.

Put

Rearr_{2}(P, P_{1}, r, v) := r ˆ+ t,  (12)

We have the following lemma.

Lemma 8.2. Assume P_{1} is an invertible selector of P. Then for r_{1} := Rearr_{2}(P, P_{1}, r, v) we have the inclusion

P · r + v ⊂ P_{1} · r_{1}.

Proof: Take u ∈ P · r+v. Then u = P · r +v for some P ∈ P, r ∈ r and v ∈ v. Let R := P − P_{1}. We have

u = P_{1} · P_{1}^{−1} ·(P · r + v) = P_{1} · (P_{1}^{−1} ·(P_{1} + R) · r + P_{1}^{−1} · v) =

P_{1} ·(r + P_{1}^{−1} · R · r + P_{1}^{−1} · v) =
P_{1} · (r + P_{1}^{−1} ·(R · r + v)) ∈ P_{1} ·(r + t) ⊂
P_{1} · r_{1}.

8.2. Parallelepipeds. We now turn to the Lohner method of parallelepipeds in enclosing f(S).

We recall that we consider a function f : ℝ^{n} → ℝ^{n} such that for a certain rational map g : ℝ^{n} → ℝ^{n} and a small continuous map w : ℝ^{n} → ℝ^{n} the decomposition (1) holds.

Let x := x + B · r ⊂ D be a fuzzy parallelepiped.

Lemma 8.3. Assume B_{1} is an interval matrix, B_{1}^{∗} is a weak inverse of B_{1}, (x_{1}, x_{1}^{0}) is a decomposition of ⟨g⟩(x) and

r_{1} := Rearr_{1}(Jxgˆ·B, B_{1}, r, w ˆ+ x_{1}^{0})

then f(x + B · r) ⊂ x_{1} + B_{1} · r_{1}.

Proof: Let u ∈ x + B · r. Then u = x + B · r for some B ∈ B and r ∈ r. Thus we have

f(u) = g(u) + w(u) ∈ g(x) + Jxg · B · r + w(u) ⊂
⟨g⟩(x) + Jxg · B · r + w ⊂
x_{1} + Jxg · B · r + w + x_{1}^{0} ⊂
x_{1} + (Jxgˆ·B) · r + (w ˆ+ x_{1}^{0})

and the conclusion follows from Lemma 8.1.

Note that an analogous lemma may be formulated in terms of Rearr_{2}. The above lemma guarantees the correctness of the enclosure for any

choice of B_{1} but obviously we still try to choose B_{1} in such a way that the the resulting enclosure is possibly small. Unless there is a special reason to do otherwise we take B_{1} as a singleton B_{1}. Also, as we explained in the preceding paragraph a good idea is to take B_{1} as a selector of Jxgˆ·B.

However there is a problem with the described method. It requires constructing the weak inverse of B_{1}. This may lead to very large interval matrices if B_{1} is ill-conditioned.

18 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

8.3. Cuboids. To avoid the problem with ill-conditioned matrices we can take cuboids instead of parallepipeds as our representable sets. A cuboid is a special case of a parallelepiped q := c+Q · r, where c ∈ ℝ^{n}, r ∈ 𝕀^{n} and Q is an interval matrix which contains an orthogonal matrix. We have the following lemma.

Let x := x + Q · r ⊂ D be a cuboid.

Lemma 8.4. If Q_{1} be an interval matrix such that Q_{1} contains an orthogonal matrix, (x_{1}, x_{1}^{0}) is a decomposition of ⟨g⟩(x) and

r_{1} := Rearr_{1}(Jxgˆ·Q, Q_{1}, r, w ˆ+ x_{1}^{0}),

then f(x + Q · r) ⊂ x_{1} + Q_{1} · r_{1}.

Proof: Let Q_{1} be an orthogonal selector of Q_{1}. Then id = Q_{1} · Q_{1}^{T} ∈ Q_{1} · Q_{1}^{T}, therefore Q_{1}^{T} is a weak inverse of Q_{1} and the result follows from the preceding lemma.

To keep wrapping small Lohner proposes to take as Q_{1} the orthogonal factor in the orthogonal-triangular decomposition of a selector of Jxgˆ·B.

8.4. Doubletons. The need for inverting matrices in the Lohner method arises from the fact that in course of computation there appear some extra interval vectors (for instance the nonlinear terms sitting in the function w), which have to be incorporated into a parallelepiped enclosure as in Lemma 8.1. However, if these interval vectors are small with respect to the size of the set S, then it is worth to keep the extra terms separately, so that the matrix inversion problem will be solved only with respect to small terms. In other words we use the sets of the form x + C · r_{0} + s, where x is a vector, C is a matrix, and r_{0} is an interval vector and s is a small set from a certain class of representable sets. We call such sets doubletons (Lohner calls the corresponding method an inner enclosure). The encoding of a representable doubleton consists of (x, C, r_{0}) and the encoding of s. Formally a doubleton is the algebraic sum of a parallelepiped and a representable set s and we want to keep them in this separated way. The parallelepiped, x + C · r_{0} is used to store the main part of the image of the set S as if there were no additional terms present in the course of the computation. Here we avoid multiplying inverted matrices by a relatively large interval set r_{0}. The set s collects all the remainders which appear during the computation. In practice we make take s in the form of an interval set, another parallelepiped, a cuboid or any other representable set in a class for which we have an algorithm encl which for a given set s in this class and a given A ∈ 𝕀^{n×n}, v ∈ 𝕀^{n} constructs a set encl(A, v, s) in this class such that

A · s + v ⊂ encl(A, v, s).

Thus let x := x + C · r_{0} + s ⊂ D be a doubleton.

Lemma 8.5. Assume B_{1} is an interval matrix, B_{1}^{∗} is a weak inverse of B_{1}, C_{1}, C_{1}^{0} is a decomposition of Jxgˆ·C, (x_{1}, x_{1}^{0}) is a decomposition of ⟨g⟩(x). Then

f(x + C · r_{0} + s) ⊂ x_{1} + C_{1} · r_{0} + encl(Jxg, w ˆ+ x_{1}^{0} ˆ+ C_{1}^{0} ˆ·r_{0}, s).

Proof: Let u ∈ x + C · r_{0} + s. Then u = x + C · r_{0} + s for some r_{0} ∈ r_{0} and s ∈ s. Thus we have

f(u) = g(u) + w(u) ∈
g(x) + Jxg · C · r_{0} + Jxg · s + w(u) ⊂
⟨g⟩(x) + Jxg · C · r_{0} + Jxg · s + w ⊂
x_{1} + C_{1} · r_{0} + Jxg · s + (w + x_{1}^{0} + C_{1}^{0} · r_{0}) ⊂
x_{1} + C_{1} · r_{0} + Jxg · s + (w ˆ+ x_{1}^{0} ˆ+ C_{1}^{0} ˆ·r_{0}) ⊂
x_{1} + C_{1} · r_{0} + encl(Jxg, w ˆ+ x_{1}^{0} ˆ+ C_{1}^{0} ˆ·r_{0}, s).

8.5. Elliptic sets. An elliptic set (an ellipsoid) in ℝ^{n} may be written as e := c + L · b, where c ∈ ℝ^{n}, L ∈ ℝ^{n×n} and b is a Euclidean ball. The elliptic set e is representable if c, L and b are representable.

Elliptic sets as a tool in the enclosing problem were introduced by A. Neumaier [13] and we refer the reader there for a very clear presentation of an enclosing problem algorithm based on elliptic sets.

8.6. Large or irregular sets. When the set S is very large or has an irregular shape, a covering strategy may be used to obtain a good upper bound. We start with a family C of certain representable sets which forms a locally finite minimal covering of ℝ^{n}. We call the sets in C elementary representable sets. Possible choices of elementary representable set are for example: a family of small representable intervals in ℝ^{n} forming a grid in ℝ^{n} or a collection of balls in ℝ^{n} of fixed representable radius and centers at a grid of representable vectors in ℝ^{n}. We then define representable sets as sets which are finite unions of elements in C. The encoding for such sets is just the list in a certain fixed order of the elements of the union. If S is represented by the union of s_{1}, s_{2}, … s_{n} then f(S) ⊂ f(s_{1}) ∪ f(s_{2}) ∪ … ∪ f(s_{n}) and the problem of enclosing f(S) is reduced to enclosing the sets f(s_{i}) for i = 1, 2, …, n.

20 MARIAN MROZEK AND PIOTR ZGLICZYŃSKI

9. Numerical examples.

We tested the methods discussed in this paper on dynamical systems obtained from:

1 the classical Rössler equations

2 the classical Lorenz equations

3 a three dimensional ODE derived form the Kuramoto-Sivashinsky

equations 4 the Hénon map In the first three cases we took as f the h-translation along the

trajectories of the flow and as its rational approximation g the 5-th order Taylor method. In case of the Hénon map we took the Hénon map itself as both f and g.

We performed two kinds of tests: blowUp and Poincare map.

BlowUp test We iterated the process of enclosing a small initial set until the ratio of the size of enclosure to the size of the initial set passed above a prescribed value (called later blowUp limit).

Poincare map test. We compute a Poincare map P for a flow induced by an ODE on some section Θ, which was relevant for the computer assisted proofs [4], [7], [16]. We are interested in two quantities. First one is L, which is given as the ratio of the size of enclosure of P (x) to the size of the initial set x. We will call L a computed Lipschitz constant.

Another quantity is a local computation cost for a given method denoted by C_{l}. We define C_{l} as follws: Let c be a time (cost) of computation of x. We set

C_{l} = Lc  (13)

We justify the above definition as follows. Our goal is to compute P (S) for some set S with a prescribed error ϵ > 0. Let S ⊂ ⋃_{i=1}^{k} x_{i}. C a total cost of computation of P (S) is given by

C = kc

Observe that the size of x_{i} is given by the following error equation

\frac{diam(x_{i})L}{2} < ϵ  (14)

From this we see that we have the following relation

k ≈ \frac{diam(S)^{d}}{diam(x)^{d}} ≥ L^{d} (\frac{diam(S)}{2ϵ})^{d}  (15)

where d is a dimension of S. Hence we see that

C ∼ L^{d}c  (16)

In out test examples we have d = 1.

Observe that in (14) we silently assumed that the only source of error in computation P (S) is related to space discretization. In real computations we have also round–off errors and errors related to time discretization. For discussion of these errors and their dependence upon various numerical methods the reader is referred to [18].

9.1. Rössler equations. Rössler equations are given by

x′ = −y − z
y′ = x + by  (17)
z′ = b + y(x − a)

where a = 5.7 and b = 0.2.

Let P be a Poincaré map for a section Θ = {x = 0, y < 0, −y − z > 0}. This is a map studied in [17]. The observed Poincaré return time belongs to (4.5, 6.5).

Tests are performed with time step h = 0.0025 for interval sets v_{1} = (0, −10.2, 0) + {0} × [−5 · 10^{−6}, 5 · 10^{−6}]^{2}, v_{2} = (0, −3.8, 0) + {0} × [−0.001, 0.001]^{2}. For blowUp test we set blowUpLimit=100 for v_{1} and blowUpLimit=10 for v_{2}.

The results of tests are in contained in tables 1 and 2. Symbol ∞ in Tab. 2 means that due to blow up of the size of the set during computations we were unable to compute Poincaré map.

The results of blowUp tests show that we have significant wrapping effect. No blow up was observed for interval and cuboid doubletons. The failure of performing so well for the parallelepiped doubleton is related to problems with inversion of ill-conditioned matrices.

In Poincaré map test we get a qualitatively different results for v_{1} and v_{2}. ¿From the point of view of rigorous numerics these sets differ significantly. v_{2} represents an example of an easy set, whereas v_{1} is an difficult set. This distinction as justified as follows: the wrapping effect during computation of P (v_{2}) is small, whereas is very big for P (v_{1}). As a result we see that for v_{2} the winner (minimal C_{l}) are euclidean logarithmic balls, but interval and parallelepiped doubletons are performing almost as good.

But for a difficult set v_{1} the doubletons cuboid and interval are nearly 10^{3}-times better then euclidean logarithmic balls.

In [17] Poincare map was computed using euclidean logarithmic balls and the total computation time was around 50 hours. ¿From the table 2 it follows that using doubletons we can redo these computation in just few minutes.

Sets  v_{1}  v_{2}

logarithmic balls:

euclidean  1.67 7.77
max  1.404 1.153
sum  1.449 1.92

balls:

euclidean  1.359 0.47
max  1.313 0.48
sum  1.215 0.39

single sets

interval  1.725 1.98
parallelepiped 4.595 3.73
cuboid  1.842 8.125
ellipsoid  4.36 3.39

doubletons

interval  > 10 > 10
parallelepiped 4.595 3.73
cuboid  > 10 > 10

Table 1. BlowUp times for Rössler system.

9.2. Tests for Lorenz equation, Kuramoto-Shivashinski equations and Hénon map. Lorenz equations are given by

x′ = s(y − x)
y′ = (r − z)x − y  (18)
z′ = xy − qz

where s = 10, r = 28 and q = 8/3. We set Θ = {z = r − 1}. We made tests for the set x = (0.137, 0.378, 27) + [−0.0005, 0.0005]^{2} × {0}. The time step h was equal to 0.003. This is an example of difficult point for computations from [4]. For blowUp test we set blowUpLimit equal to 100.

We made test also for the following three dimensional ODE derived from Kuramoto-Sivashinsky equations (see [12])

x′ = y
y′ = z  (19)
z′ = −d^{2}ly − x^{2}/2 + d^{6}

where d = 1 and l = 1.

Sets  L(v_{1})  C_{l}(v_{1})  L(v_{2}) C_{l}(v_{2})

logarithmic balls:

euclidean  1.33 · 10^{4} 1  3.53 1
max  ∞  ∞  ∞  ∞
sum  ∞  ∞  ∞  ∞

balls:

euclidean  ∞  ∞  ∞  ∞
max  ∞  ∞  ∞  ∞
sum  ∞  ∞  ∞  ∞

single sets:

interval  4440  0.71  1530 950
parallelepiped ∞  ∞  ∞  ∞
cuboid  435.7  0.073  3.23 2.04
ellipsoid  ∞  ∞  ∞  ∞

doubletons:

interval  4.99  8.1 · 10^{−4} 2.26 1.39
parallelepiped ∞  ∞  ∞  ∞
cuboid  4.96  8.4 · 10^{−4} 2.21 1.41

Table 2. Results of Poincare map test for Rössler system

Hénon map is given by

H(x, y) = (1 − ax^{2} + y, bx)  (20)

where a = 1.4 and b = 0.3. We do test for x = (1, 1) + [−0.5 · 10^{−5}, 0.5 · 10^{−5}]^{2}. For blowUp test we set blowUpLimit equal to 100.

The table 3 collects the results of blowUp test. ¿From the contents of tables 3 and 4 one can see that the winner for

Lorenz equation is doubleton cuboid , but it is only three times better than euclidean logarithmic balls, which were used in [4].

References

[1] R. Anguelov, S. Markov, Wrapping effect and wrapping function, Reliable Computing 4(1998), 1-20.

[2] R. Anguelov, Wrapping function of the initial value problem for ODE: Applications, Reliable Computing, to appear.

[3] George F. Corliss, Robert Rihm, Validating an a priori enclosure using highorder Taylor series, Scientific computing and validated numerics (Wuppertal, 1995), 228–238, Math. Res., 90, Akademie Verlag, Berlin, 1996.

[4] Z. Ga1ias and P. Zg1iczyński, Computer assisted proof of chaos in the Lorenz system, Physica D, 115, 1998,165-188.

Set  Lorenz K-S Hénon
interval  0.59 2.28  13
parallelepiped  0.42 2.71  7
cuboid  1.14 2.63  13
ellipsoid  0.44 2.84  8

doubletons:

interval  1.23 2.87  13
parallelepiped  0.42 2.71  7
cuboid  2.54 2.86  13

balls:

euclidean  0.28 1.71  8
max  0.23 1.56  5
sum  0.25 1.68  8

logarithmic balls:

euclidean  0.61 2.14
max  0.41 1.56
sum  0.45 1.68

Table 3. BlowUp test results

[5] E. Hairer, S.P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems, Springer-Verlag, Berlin Heidelberg 1987.

[6] R.J. Lohner, Computation of Guaranteed Enclosures for the Solutions of Ordinary Initial and Boundary Value Problems, in: Computational Ordinary Differential Equations, J.R. Cash, I. Gladwell Eds., Clarendon Press, Oxford, 1992.

[7] K. Mischaikow and M. Mrozek, Chaos in Lorenz equations: a computer assisted proof, Bull. Amer. Math. Soc. (N.S.), 32(1995), 66-72.

[8] K. Mischaikow and M. Mrozek, Chaos in the Lorenz Equations: A Computer Assisted Proof. Part II: Details, Mathematics of Computation, 67(1998), 1023-1046.

[9] K. Mischaikow, M. Mrozek and A. Szymczak, Chaos in the Lorenz equations: a computer assisted proof. Part III: the classical parameter values, submitted.

[10] R.E. Moore, Interval analysis, Prentice-Hall, Inc., Englewood Cliffs, N.J., 1966.

[11] M. Mrozek, Topological invariants, multivalued maps and computer assisted proofs, Computers & Mathematics, 32(1996),83-104.

[12] M. Mrozek and M. Żelawski, Heteroclinic connections in the Kuramoto-Sivashinsky equation, Reliable Computing, 3(1997), 277-285.

[13] A. Neumaier, The wrapping effect, ellipsoid arithmetic, stability and confidence regions, Computing Supplementum, 9 (1993), 175-190.

[14] M. Warmus, Calculus of approximations, Bull. de l’Académie Polonaise des Sciences Cl. III, 4(1956), 253-259.

Set  Lorenz K-S

logarithmic balls:

euclidean  1.00 1.00
max  ∞ 2.14
sum  ∞ 3.39

sets:

interval  62.9 5.08
parallelepiped  ∞ 2.86
cuboid  0.59 3.66
ellipsoid  ∞ 2.42

doubletons:

interval  0.39 2.82
parallelepiped  ∞ 2.90
cuboid  0.31 2.93

balls:

euclidean  ∞ 2.46
max  ∞ 2.17
sum  ∞ 3.47

Table 4. Results for Poincaré map tests, local computation cost comparison

[15] M. Warmus, Approximation and Inequalities in the Calculus of Approximations. Classification of Approximate Numbers, Bull. de l’Académie Polonaise des Sciences, Série des sciences math., astr. et phys., 9, (1961), 241-245.

[16] P. Zgliczyński, Rigorous verification of chaos in the Rössler equations, Scientific Computing and Validated Numerics, G. A1efe1d, A. Frommer and B. Lang editors, Akademie Ver1ag, Ber1in, 1996, 287-292

[17] P. Zgliczyński, Computer assisted proof of chaos in the Hénon map and in the Rössler equations, Nonlinearity, 10(1997), 243-252.

[18] P. Zgliczyński, Remarks on computer assisted proof of chaotice behavior in ODE’s, in preparation.

Instytut Informatyki, Uniwersytet Jagielloński, Kraków, Poland

E-mail address: mrozek@ii.uj.edu.pl

School of Mathematics, Georgia Institute of Technology, Atlanta, Georgia, USA.

E-mail address: piotrz@math.gatech.edu