Cr-Lohner algorithm
Daniel Wilczak1, Piotr Zgliczyński2
Jagiellonian University, Institute of Computer Science,
Nawojki 11, 30–072 Kraków, Poland
e-mail: wilczak@ii.uj.edu.pl, umzglicz@cyf-kr.edu.pl
November 3, 2018
Abstract
We present a Lohner type algorithm for the computation of rigorous
bounds for solutions of ordinary differential equations and its derivatives
with respect to initial conditions up to arbitrary order. As an applica-
tion we prove the existence of multiple invariant tori around some elliptic
periodic orbits for the pendulum equation with periodic forcing and for
Michelson system.
1 Introduction
This paper is a sequel to [Z]. We present here a Lohner-type algorithm for
computation of rigorous enclosures of partial derivatives with respect to initial
conditions up to an arbitrary order r of the flow induced by an autonomous
ODE, hence the name Cr-Lohner algorithm. Let r be a positive integer, then
by Cr-algorithm we will mean the routine which gives rigorous estimates for
partial derivatives with respect to initial conditions up to an order r and Cr-
computations we mean an application of an Cr-algorithm.
Our main motivation for the development of Cr-algorithm was a desire to
provide a tool, which will considerably extend the possibilities of computer
assisted proofs in the dynamics of ODEs. Till now most of such proofs have used
topological conditions (see for example [HZHT, MM, GZ, Z1]) and additionally
conditions on the first derivatives with respect to initial conditions (see for
example [RNS, T, Wi1, WZ, KZ]), hence it required C0- and C1-computations,
respectively. The spectrum of problems treated includes the questions of the
existence of periodic orbits and their local uniqueness, the existence of symbolic
dynamics, the existence of hyperbolic invariants sets, the existence of homo-
and heteroclinic orbits. To treat other phenomena, like bifurcations of periodic
orbits, the route to chaos, invariant tori through KAM theory one needs the
knowledge of partial derivatives with respect to initial conditions of higher order.
In principle, one can think that a good rigorous ODE solver should be en-
ough. Namely, to compute the partial derivatives of the flow induced by
x′ = f(x), x ∈ Rn (1)
1 Research supported by an annual national scholarship for young scientists from the
Foundation for Polish Science
2 Research supported in part by Polish State Ministry of Science and Information Techno-
logy grant N201 024 31/2163
http://arxiv.org/abs/0704.0720v1
it is enough to rigorously integrate a system of variational equations obtained by
a formal differentiation of (1) with respect to the initial conditions. For example
for r = 2 we have the following system
x′ = f(x), (2)
Vij(t) =
(x)Vsj(t) (3)
Hijk(t) =
s,r=1
∂xs∂xr
(x)Vrk(t)Vsj(t) +
(x)Hsjk(x), (4)
with the initial conditions
x(0) = x0, V (0) = Id, Hijk(0) = 0, i, j, k = 1, . . . , n.
It is well known that if by ϕ(t, x0) we denote the (local) flow induced by (1),
(t, x0) = Vi,j(t),
∂xj∂xk
(t, x0) = Hijk(t).
Analogous statements are true for higher order partial derivatives with respect
to initial conditions.
It turns out that a straightforward application of a rigorous ODE solver to
the system of variational Equations (2–4) is very inefficient. Namely, it totally
ignores the structure of the system and leads to a very poor performance and
unnecessary long computation times (see Section 4.1).
Our algorithm is a modification of the Lohner algorithm [Lo], which takes
into account the structure of variational Equations (2–4). Basically it consists of
the Taylor method, a heuristic routine for a priori bounds for solution of (2–4)
during a time step and a Lohner-type control of the wrapping effect, which is
done separately for x and partial derivatives with respect initial conditions (the
variables V and H in (3,4)). The Taylor method is realized using the automa-
tic differentiation [Ra] and the algorithms for computation of compositions of
multivariate Taylor series.
The proposed algorithm has been successfully applied in [HNW] to the Mi-
chelson system [Mi], where a computer assisted proof of the existence of a cocoon
bifurcation was presented. Some parts of this proof required C2-computations.
In the present paper in Section 8 we show an application of our algorithm to
pendulum equation with periodic forcing and the Michelson system. We used
it to compute rigorous bounds for the coefficients of some normal forms up to
order five, which enabled us to prove the existence of invariant tori around some
elliptic periodic orbits in these systems using KAM theorem for twist maps on
the plane. These proofs required C3 and C5 computations.
2 Basic definitions
To effectively deal with the formulas involving partial derivatives we will use ex-
tensively a notation of multiindices, multipointers and submultipointers throu-
ghout the paper.
As an motivation let us consider the formula for the partial derivatives of
the composition of maps. Assume g : Rn → Rn and f : Rn → R are of class C3.
We have
∂3(f ◦ g)
∂xi∂xj∂xc
k,r,s=1
∂xk∂xr∂xs
∂xi∂xj∂xc
k,r=1
∂xk∂xr
∂xi∂xc
∂xj∂xc
∂xi∂xj
To the operator ∂
∂xi1∂xi2∂xi3
we can in a unique way assign a multipointer,
which is a nondecreasing sequence of integers (j1, j2, j3), such that {i1, i2, i3} =
{j1, j2, j3}. A submultipointer is a multipointer, which is a part of a longer mul-
tipointer, for example (i, j, c)(1,3) = (i, c). One observes, that submultipointers
appear at several places in the above formula.
A multiindex is an element of α ∈ Nn. It is another way to represent various
partial derivatives. The coefficient αi tells us how many times to differentiate
a function with respect to the i-th variable. Obviously, we have one-to-one
correspondence between multipointers and multiindices.
2.1 Multiindices
By N we will denote the set of nonnegative integers, i.e. N = {0, 1, 2, . . .}.
Definition 1 An element τ ∈ Nn will be called a multiindex.
For a sequence α = (α1, . . . , αn) ∈ Nn and a vector x = (x1, . . . , xn) ∈ Rn we
1. |α| = α1 + · · ·+ αn
2. α! = α1! · α2! · · ·αn!
3. xα = (xα11 , . . . , x
By eni ∈ Nn we will denote
eni = (0, 0, . . . , 0,
1 , 0, . . . , 0, 0).
We will drop the index n (the dimension) in the symbol eni when it is obvious
from the context.
Put Nnp := {a ∈ Nn : |a| = p}.
For δ = (δ1, . . . , δk) ∈ Nn1 × · · · × Nnk we set
1. |δ| =
i=1 |δi|
2. δ! =
i=1 δi!
Let f = (f1, . . . , fm) : R
n → Rm be sufficiently smooth. For α ∈ Nn we set
1. Dαfi =
∂|α|fi
∂xα11 · · · ∂x
2. Dαf = (Dαf1, D
αf2, . . . , D
For a function f : R× Rn → Rn by Dαfi(t, x) we will denote Dαfi(t, ·)(x) and
similarly
Dαf(t, x) = (Dαf1(t, x), . . . , D
αfn(t, x)).
This convention means that Dα always acts on x-variables.
2.2 Multipointers
For a fixed n > 0 and p > 0 we define
Nnp := {(a1, a2, . . . , ap) ∈ Np : 1 ≤ a1 ≤ · · · ≤ ap ≤ n}
N = Nn :=
Definition 2 An element of Nn will be called a multipointer.
Remark 3 A function
Λ : Nnp ∋ (a1, . . . , ap) →
enai ∈ N
p (5)
is a bijection.
Let f = (f1, . . . , fm) : R
n → Rm be a sufficiently smooth. For a ∈ Nnp we set
1. Dafi :=
∂xa1 . . . ∂xap
2. Daf := (Daf1, . . . , Dafm)
For a function f : R × Rn → Rn by Dafi(t, x) we will denote Dafi(t, ·)(x). In
the light of the above notations Dαf = D
Λ(α)f .
For a = (a1, a2, . . . , an) ∈ Nnp and b = (b1, b2, . . . , bn) ∈ Nnq we define
a+ b = (a1 + b1, . . . , an + bn) ∈ Nnp+q.
For α ∈ Nnp and β ∈ Nnq we define
α+ β = Λ−1 (Λ(α) + Λ(β)) ∈ Nnp+q.
By ≤ we will denote a linear order (lexicographical order) in N defined in
the following way. For a ∈ Nnp and b ∈ Nnq
(a ≤ b) ⇐⇒
either ∃i, i ≤ p, i ≤ q, ai < bi and aj = bj for j < i
or p ≤ q and ai = bi for i = 1, . . . , p.
Definition 4 For k ≤ p we set
N p(k) := {(δ1, . . . , δk) ∈ (N p)k : δ1 ≤ · · · ≤ δk, δ1+ · · ·+δk = (1, 2, . . . , p)} (7)
We will use N p(k) extensively in the next section. Its will be used to label
terms in Dαfi(ϕ(t, x)). Observe that for p > 0
N p(1) = {(1, 2, . . . , p)}
N p(p) = {((1), (2), . . . , (p))}
One can construct all elements of N p(k) using the following recursive procedure.
From the definition of N p(k) it follows that if (δ1, . . . , δm−1) ∈ N p−1(m − 1)
then (δ1, . . . , δm−1, (p)) ∈ N p(m) (notice that order is preserved). Similarly, if
(δ1, . . . , δm) ∈ N p−1(m) then
(δ1, . . . , δs−1, δs + (p), δs+1, . . . , δm) ∈ N p(m)
and again order of elements is preserved. Hence, for p > 2 and 1 < k < p we
have N p(k) = A ∪B where
(δ1, . . . , δk−1, (p)) : (δ1, . . . , δk−1) ∈ N p−1(k − 1)
(δ1, . . . , δs−1, δs + (p), δs+1, . . . , δk) : (δ1, . . . , δk) ∈ N p−1(k)
} (8)
and the sets A and B are disjoint.
Another way to generate all elements of N p(k) can be described as follows
• decompose the set {1, 2, . . . , p} into k nonempty and disjoints sets ∆i,
i = 1, . . . , k
• we sort each ∆i and permute ∆i’s to obtain min(∆1) < min(∆2) < · · · <
min(∆k)
• we define δi to be an ordered set consisting of all elements of ∆i for
i = 1, . . . , k
Definition 5 For an arbitrary a ∈ Nnp and δ ∈ N
k such that k ≤ p we define a
submultipointer aδ ∈ Nnk by (aδ)i = aδi for i = 1, . . . , k, which can be expressed
using Λ as follows
aδ := Λ
enaδi
∈ Nnk
3 Equations for variations
Consider an ODE x′ = f(x) where f is CK+1. Let ϕ : R× Rn−→◦ Rn be a local
dynamical system induced by x′ = f(x). It is well known, that ϕ ∈ CK and one
can derive the equations for partial derivatives of ϕ by differentiating equation
(t, x) = f(ϕ(t, x)) with respect to the initial condition x. As a result we
obtain a system of so-called equations for variations, whose size depends on the
order r of partial derivatives we intend to compute. An example of such system
for r = 2 is given by (2–4) with initial conditions given by (5).
The goal of this section is to write the equations for variations in a compact
form using multipointers and multiindices, which allows us to take into account
the symmetries of partial derivatives,
Lemma 6 Assume f ∈ Cr+1 and let ϕ : R × Rn−→◦ Rn be a local dynamical
system induced by x′ = f(x). Then for a ∈ Nnp such that p ≤ r holds
Daϕi =
i1,...,ik=1
Dei1+···+eik fi
(δ1,...,δk)∈Np(k)
Daδjϕij (9)
for i = 1, . . . , n.
Proof: In the proof the functions Dei1+···+eik fi are always evaluated at
ϕ(t, x), and various partial derivatives of ϕ are always evaluated at (t, x), there-
fore the arguments will be always dropped to simplify formulae. We prove the
lemma by induction on p = |a|. If p = 1 then a = (c) for some c ∈ {1, . . . , n}
and (9) becomes
D(c)ϕi =
Desfi ·D(c)ϕs.
Assume (9) holds true for p − 1, p > 1. Let us fix a ∈ Nnp . We have
a = b + (c), where b = (a1, . . . , ap−1) ∈ Nnp−1 and c = ap. Since (9) is satisfied
for p− 1, therefore we have
Daϕi = D(c)
= D(c)
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδjϕij
i1,...,ik+1=1
β:=ei1+···+eik+1
Dβfi ·D(c)ϕik+1
(δ1,...,δk)∈Np−1(k)
Dbδjϕij
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
For k = 1, . . . , p we set
Tk :=
i1,...,ik=1
Dei1+···+eik fi
(δ1,...,δk)∈Np(k)
Daδjϕij (10)
Now our goal is to prove that:
Daϕi =
Tk (11)
Our strategy of proof is as follows. We will define S1, . . . , Sp, such that
Daϕi =
Sk, Si = Ti, i = 1, . . . , p. (12)
We set
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
k=p−1
i1,...,ik+1=1
β:=ei1+···+eik+1
Dβfi ·D(c)ϕik+1
(δ1,...,δk)∈Np−1(k)
Dbδjϕij .
For m = 2, 3, . . . , p− 1 we set
k=m−1
i1,...,ik+1=1
β:=ei1+···+eik+1
Dβfi ·D(c)ϕik+1
(δ1,...,δk)∈Np−1(k)
Dbδjϕij
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
It remains to show that Si = Ti for i = 1, . . . , p. Consider first i = 1. Recall
that N p−1(1) = {(1, 2, . . . , p− 1)}, hence
Desfi ·Db+(c)ϕs =
Desfi ·Daϕs.
Therefore
S1 = T1. (13)
Consider now i = p. For an arbitrary s > 0 N s(s) contains only one element
((1), (2), . . . , (s)). Therefore we obtain
i1,...,ip=1
Dei1+···+eip fi ·D(c)ϕip
(δ1,...,δp−1)∈Np−1(p−1)
Dbδjϕij
i1,...,ip=1
Dei1+···+eip fi ·D(c)ϕip
Dbjϕij .
Since a = b+ (c), where c = (ap), hence
i1,...,ip=1
Dei1+···+eip fi
Dajϕij
i1,...,ip=1
Dei1+···+eip fi
(δ1,...,δp)∈Np(p)
Daδjϕij = Tp
Consider now m = 2, 3, . . . , p− 1. We have
i1,...,im=1
Dei1+···+eim fi ·D(c)ϕim
(δ1,...,δm−1)∈Np−1(m−1)
Dbδjϕij
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm)∈Np−1(m)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
Using decomposition N p(m) = A ∪B as in (8) we obtain
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm−1,δm=(p))∈A
Daδjϕij
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm)∈B
Daδjϕij
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm)∈Np(m)
Daδjϕij = Tm
We have shown that Ti = Si for i = 1, . . . , p. This finishes the proof.
4 Cr-Lohner algorithm
4.1 Why one needs an Cr-algorithm?
There are several effective algorithms for the computation of rigorous bounds
for solutions of ordinary differential equations, including Lohner method [Lo],
Hermite–Obreschkoff algorithm [NJ] or Taylor models [BM]. For Cr-computa-
tions the number of equations to solve is equal to n
hence, even for
r = 1 direct application of such an algorithms to equations for variations (14)
leads to integration in high dimensional space and is usually inefficient. Let us
recall after [Z, Sec. 6] the basic reason for this. In order to have a good control
over the expansion rate of the set of initial conditions during a time step these
algorithms, while being C0, are C1 ’internally’(or higher for Taylor models),
because they solve non-rigorously equations for (∂ϕ
) - the variational matrix of
the flow. This effectively squares the dimension of phase space of the equation
and impacts heavily the computation time. But as it was observed in [Z] the
equations for partial derivatives of the flow can be seen as non-autonomous and
nonhomogenous linear equations, therefore we do not need additional equations
for variations for them. As a result the dimension of the effective phase space
for our Cr-algorithm is given by n
and not a square of this number.
Another important aspect of the proposed algorithm is the fact that the
Lohner-type control of the wrapping effect is done separately for x-variables
and variables Daϕ. This feature is not present in the blind application of C0
algorithm to the system of variational equations and it turns out that this often
practically switches off the control of the wrapping effect on x-variables, as
various choices used in this control become dominated by the Daϕ-variables.
In [Z] a C1-algorithm has been proposed. Here we present an algorithm for
computation of higher order partial derivatives.
4.2 An outline of the algorithm
Let us fix r ≤ K and consider the following system of differential equations
ϕ = f ◦ ϕ
Daϕ =
i1,...,ik=1
Dei1+···+eik f
(δ1,...,δk)∈Nd(k)
Daδjϕij
for all a ∈ Nnd , d = 1, . . . , r.
Our goal is to present an algorithm for computing a rigorous bound for the
solution of (14) with a set of initial conditions
ϕ(0, x0) ∈ [x0] ⊂ Rn
Dϕ(0, x0) = Id
Daϕ(0, x0) = 0, for a ∈ Nn2 ∪ . . . ∪ Nnr .
In the sequel we will use the following notations:
• if a solution of system (14) is defined for t > 0 and some x0 ∈ Rn, then
for a ∈ N by Va(t, x0) we denote Daϕ(t, x0)
• for [x0] ⊂ Rn by [Va(t, [x0])] we will denote a set for which we have
Va(t, [x0]) ⊂ [Va(t, [x0])]. This set is obtained using an rigorous nume-
rical routine described below.
The Cr-Lohner algorithm is a modification of C1-Lohner algorithm [Z]. One
step of Cr-Lohner is a shift along the trajectory of the system (14) with the
following input and output data
Input data:
• tk - a current time,
• hk - a time step,
• [xk] ⊂ Rn, such that ϕ(tk, [x0]) ⊂ [xk],
• [Vk,a] = [Vk,a(tk, [x0])] ⊂ Rn, such that Daϕ(tk, [x0]) ⊂ [Vk,a] for a ∈
Nn1 ∪ . . . ∪ Nnr .
Output data:
• tk+1 = tk + hk - a new current time,
• [xk+1] ⊂ Rn, such that ϕ(tk+1, [x0]) ⊂ [xk+1],
• [Vk+1,a] = [Vk+1,a(tk+1, [x0])] ⊂ Rn, such that Daϕ(tk+1, [x0]) ⊂ [Vk+1,a]
for a ∈ Nn1 ∪ . . . ∪ Nnr .
We will often skip the arguments of Vk,a when they are obvious from the context.
The values of [xk+1] and [Vk+1,a], a ∈ Nn1 are computed using one step
C1-Lohner algorithm. After it is done, we perform the following operations to
compute [Vk+1,a] for a ∈ Nn2 ∪ . . . ∪ Nnr
1. Find a rough enclosure for Daϕ([0, hk], [xk]).
2. Compute [Vk+1,a], this will also involve some rearrangement computations
to reduce the wrapping effect for V [Mo, Lo].
5 Computation of a rough enclosure for Daϕ
For a fixed multipointer a ∈ Nnd Equation (14) can be written as follows
Daϕ(t, x) = Ba(t, x) +A(t, x)Daϕ(t, x) (16)
where
i1,...,ik=1
Dei1+···+eik f
(δ1,...,δk)∈Nd(k)
Daδjϕij
A = Df ◦ ϕ
The procedure for computing the rough enclosure is based on the notion of
a logarithmic norm, which we give below.
Definition 7 [HNW] For a square matrix A the logarithmic norm µ(A) is de-
fined as a limit
µ(A) = lim sup
‖Id +Ah‖ − 1
where ‖ · ‖ is a given matrix norm.
The formulas for the logarithmic norm of a real matrix in the most frequently
used norms are (see [HNW])
1. for ‖x‖1 =
i |xi|, µ(A) = maxj(ajj +
i6=j |aij |)
2. for m ‖x‖2 =
i |xi|2, µ(A) is equal to the largest eigenvalue of (A +
AT )/2
3. for ‖x‖∞ = maxi |xi|, µ(A) = maxi(aii +
j 6=i |aij |)
In order to find bounds for Daϕ we use the following theorem [HNW, Thm.
I.10.6]
Theorem 8 Let x(t) be a solution of a differential equation
x′(t) = f(t, x(t)), x ∈ Rn (18)
Let ν(t) be a piecewise differentiable function with values in Rn. Assume that
(t, η)
≤ l(t) for η ∈ [x(t), ν(t)]
|ν′(t)− f(t, ν(t))| ≤ δ(t),
where by µ(A), we denote a logarithmic norm of a square matrix A ∈ Rn×n.
Then for t ≥ t0 we have
|x(t) − ν(t)| ≤ eL(t)
|x(t0)− ν(t0)|+
e−L(s)δ(s)ds
, (19)
with L(t) =
l(τ)dτ .
We apply the above theorem to Equation (16) to obtain
Lemma 9 Let us fix x ∈ Rn. Assume that |Ba(t, x)| ≤ δ(t) and µ(A(t, x)) ≤
l(t), then for t > t0
|Daϕ(t, x)| ≤ |Daϕ(t0, x)|eL(t) + eL(t)
e−L(τ)δ(τ)dτ (20)
with L(t) =
l(τ)dτ .
Proof: Consider Equation (16) and a homogenous problem for (16)
w = f(t, w) := A(t, x) · w, w ∈ Rn. (21)
Using Theorem 8 we can estimate the difference between any solution of (21),
w, and a solution of (16), denoted by Daϕ.
|Daϕ(t) − w(t)| ≤ |Daϕ(t0)− w(t0)|eL(t) + eL(t)
e−L(τ)δ(τ)dτ. (22)
After a substitution w(t) = 0, which is a solution of the homogenous equation,
we obtain our assertion.
Usually, we do not have any control over the time dependence of δ and l,
hence we will use the following
Lemma 10 Assume that |Ba(t, x)| ≤ δ and µ(A(t, x)) ≤ l for t ∈ [0, h] then
for t ∈ [0, h] we have
|Daϕ(t, x)| ≤ |Daϕ(0, x)|max(1, ehl) + δ
elt − 1
, if l 6= 0, (23)
|Daϕ(t, x)| ≤ |Daϕ(0, x)|+ δt, when l = 0. (24)
5.1 The procedure for the computation of the rough en-
closure for V .
The procedure for the computing of the rough enclosure is iterative, which
means that given a rough enclosure for ϕ([0, hk], [xk]) and rough enclosures
Daϕ([0, hk], [xk]) for all a ∈ Nn1 ∪ . . . ∪ Nnp we are able to compute the rough
enclosure for Daϕ([0, hk], [xk]) for a ∈ Nnp+1.
The procedures for computation of the rough enclosures of ϕ([0, hk], [xk])
and Daϕ([0, hk], [xk]) for a ∈ Nn1 has been given in [Z]. Below we present an
algorithm for computing [Ea] for a ∈ Nn2 ∪ . . . ∪Nnr .
Input parameters:
• hk - a time step,
• [xk] ⊂ Rn - the current value of x = ϕ(tk, [x0]),
• [E0] ⊂ Rn - a compact and convex such that ϕ([0, hk], [xk]) ⊂ [E0]
• [Ea] ⊂ Rn, a ∈ Nn1 ∪ . . . ∪ Nnp such that Daϕ([0, hk], [xk]) ⊂ [Ea] for
a ∈ Nn1 ∪ . . . ∪ Nnp .
Output:
• [Ea] ⊂ Rn, a ∈ Nnp+1 such that
Daϕ([0, hk], [xk]) ⊂ [Ea]
Before we present an algorithm let us observe that for a fixed a ∈ Nnp+1, Ba
defined in (17) could be seen as a multivariate function of t, x and Vb = Dbϕ for
b ∈ Nn1 ∪. . .∪Nnp . More precisely, putmp := ♯
Nn1 ∪ . . . ∪ Nnp
, where ♯ stands
for number of elements of a set. Recall that, we have defined by (6) a linear
order in Nn. Hence, there is a unique sequence of multipointers b1, . . . , bmp ,
such that bi ∈ Nn1 ∪ . . .∪Nnp for i = 1, . . . ,mp, b1 ≤ b2 ≤ · · · ≤ bmp and bi 6= bj
for i 6= j.
Let us define
B̃a : R× (Rn)mp+1 → Rn,
Fa : R× (Rn)mp+1 → Rn
B̃a(t, x, vb1 , . . . , vbmp ) =
i1,...,ik=1
Dei1+···+eik f(ϕ(t, x))
(δ1,...,δk)∈Np+1(k)
Fa(t, x, vb1 , . . . , vbm) = B̃a(t, x, vb1 , . . . , vbm) +Df(ϕ(t, x))Va(t, x) (26)
Algorithm:
To compute [Ea] for a ∈ Nnp+1 we proceed as follows
1. Find l ≥
maxx∈[E0] µ (Df(x))
2. Compute δa ≥ max ‖B̃a‖, i.e.
δa ≥ max
(x,vb1 ,...,vbmp )∈[E0]×[Eb1 ]×···×[Ebmp ]
∥B̃a(0, x, vb1 , . . . , vbmp )
For example, if a = (j, c) ∈ Nn2 , then δa should be such that
δa ≥ max
x∈[E0],v1∈[E(1)],...,vn∈[E(n)]
r,s=1
∂xr∂xs
(x) (vj)s (vc)r
3. Define [Ea]i = [−1, 1]δa e
, for i = 1, . . . , n, where [Ea]i denotes i-th
coordinate of [Ea].
One can refine the obtained enclosure by
[Ea] :=
[0, hk]Fa
0, [E0], [Eb1 ], . . . , [Ebmp ]
∩ [Ea]
Indeed, for i = 1, . . . , n, t ∈ [0, hk] and x0 ∈ [E0] we have
Daϕi(t, x0) = Daϕi(t, x0)−Daϕ(0, x0)
= t (Fa)i (θi, x0, Db1ϕ(θi, x0), . . . , Dbmpϕ(θi, x0))
= t (Fa)i (0, ϕ(θi, x0), Db1ϕ(θi, x0), . . . , Dbmpϕ(θi, x0))
for some θi ∈ [0, t] ⊂ [0, hk]. In the above we have used the fact that
Fa(t, x, v1, . . . , vmp) = Fa(0, ϕ(t, x), v1, . . . , vmp).
Since ϕ(θi, x0) ∈ [E0] and Dbjϕ(θi, x0) ∈ [Ebj ] for j = 1, . . . ,mp we get
Daϕi(t, x0) ∈ [0, hk] (Fa)i
0, [E0], [Eb1 ], . . . , [Ebmp ]
6 Computation of [Vk+1]
6.1 Composition formulas
For any p-times continuously differentiable functions f, g : Rn → Rn and a ∈ Nnp
we have
Da(f ◦ g) =
i1,...,ik=1
Dei1+···+eik fi
(δ1,...,δk)∈Np(k)
Daδj gij (27)
We can apply the above formula to f = ϕ(hk, ·) and g = ϕ(tk, ·) to obtain
Va(tk + hk, x0) =
i1,...,ik=1
VΛ−1(ei1+...+eik )
(hk, xk)
(δ1,...,δk)∈Np(k)
(tk, x0)
for all x0 ∈ [x0]. Using notations [Vk+1,a] := [Va(tk + hk, [x0])] and [Vk,a] =
[Va(tk, [x0])] we can rewrite the above equation as
[Vk+1,a] =
i1,...,ik=1
VΛ−1(ei1+...+eik )
(hk, [xk])
(δ1,...,δk)∈Np(k)
Vk,aδj
where Λ is defined by (5).
6.2 The procedure for computation of [Vk+1]
We introduce new parameters od - the order of the Taylor method used in
computations of Va for a ∈ Nnd . It makes sense to take o1 ≥ o2 ≥ · · · ≥ or.
Input parameters:
• hk - a time step,
• [xk] ⊂ Rn - the current value of x = ϕ(tk, [x0]),
• [Vk,a] ⊂ Rn - a current value of Vk,a(tk, [x0]), for a ∈ Nn1 ∪ . . . ∪Nnr
• [E0] ⊂ Rn compact and convex, such that ϕ([0, hk], [xk]) ⊂ [E0] - a rough
enclosure for [xk],
• [Ea] ⊂ Rn, compact and convex, such that Daϕ([0, hk], [xk]) ⊂ [Ea], for
a ∈ Nn1 ∪ . . . ∪ Nnr .
Output: [Vk+1,a] ⊂ Rn, such that
Va(tk + hk, x0) ∈ [Vk+1,a] (29)
for x0 ∈ [x0] and a ∈ Nn1 ∪ . . . ∪ Nnr .
Algorithm: We compute [Vk+1] as follows
1. Computation of Va(hk, [xk]) using Taylor method for Equation (14), i.e. for
a ∈ Nnp we compute
[Fa] =
dti−1
Fa(0, [xk], Vb1 , . . . , Vbmp−1 ) (30)
hop+1
(op + 1)!
Fa(0, [E0], [Eb1 ], . . . , [Ebmp−1 ]).
where Vbi = 0 for bi ∈ Nn2 ∪ . . . ∪ Nnp−1 and V(j) = enj for j = 1, . . . , n.
Observe that
Va(hk, [xk]) ⊂ [Fa] (31)
Indeed, using Taylor series expansion we obtain that for xk ∈ [xk] and
j = 1, . . . , n holds
(Va)j(hk, xk) =
dti−1
(Fa)j(0, xk, Vb1(0, xk), . . . , Vbmp−1 (0, xk))
hop+1
(op + 1)!
(Fa)j(θi, xk, Vb1(θi, xk), . . . , Vbmp−1 (θi, xk))
for some θi ∈ [0, hk]. Observe, that
(Fa)j(θi, xk, Vb1(θi, xk), . . . , Vbmp−1 (θi, xk))
(Fa)j(0, ϕ(θi, xk), Vb1(θi, xk), . . . , 0, Vbmp−1 (θi, xk))
Using ϕ(θi, xk) ∈ [E0] and Vbs(θi, xk) ∈ [Ebs ] for s = 1, . . . ,mp−1 we
obtain our assertion.
2. The composition. Put
[Jk] := ([F(1)], . . . , [F(n)])
Using (28) for a ∈ Nnp we have
[Vk+1,a] = [αa] + [Jk] · [Vk,a], (32)
where
[αa] =
i1,...,ik=1
[FΛ−1(ei1+...+eik )
(δ1,...,δk)∈Np(k)
Vk,aδj
In our implementation of the algorithm we use the symbolic differentiation to
obtain formulae for Daf . Next, using the automatic differentiation we compute
Fa(t, x, Vb1 (t, x), . . . , Vbmp−1 (t, x))|t=0 which appear in (30).
6.3 Rearrangement for Va - the evaluation of Equation
It is well know that a direct evaluation of Equation (32) leads to wrapping effect
[Mo, Lo]. To avoid it following the work of Lohner [Lo] we will use the same
scheme as it was proposed in [Z].
Namely, observe that Equation (32) has exactly the same structure as the
propagation equations for C1-method (see [Z, Section 3]). Moreover, all vectors
Vk,a, for a ∈ Nn1 ∪ . . .Nnr ’propagate’ by the same [Jk] as did the variational
part in [Z], hence it makes sense the same approach.
To be more precise, each set [Vk,a], for a ∈ Nn1 ∪ . . . ∪ Nnr is represented in
the following form
[Vk,a] = vk,a + [Bk][rk,a] + Ck[qk,a]
where [Bk] is interval matrix, Ck is point matrix, vk,a is a point vector and
rk,a, qk,a are interval vectors. Observe that [Bk] and Ck are independent of a.
In the sequel we will drop index a. Equation (32) leads to
[Vk+1] = [α] + [Jk](vk + [Bk][rk] + Ck[qk]) (34)
Let m([z]) denotes a center of an interval object, i.e. [z] is interval vector or
interval matrix and ∆([z]) = [z]−m([z]).
Let [Q] be an interval matrix which contains an orthogonal matrix. Usually,
[Q] is computed by the orthonormalisation of the columns of m([Jk])[Bk].
[Z] = m([Jk])Ck
Ck+1 = m([Z])
[Bk+1] = [Q]
Then we rearrange formula (34) as follows
[s] = [α] + [Jk]vk +∆([Jk])[Vk]
vk+1 = m([s])
[qk+1] = [qk]
[rk+1] = [Q
T ](∆([s]) + ∆([Z])[qk]) + ([Q
T ]m([Jk])[Bk])[rk]
Summarizing, we can use the following data structure to represent ϕ(tk, [x0])
and Daϕ(tk, [x0]), for a ∈ Nn1 ∪ . . . ∪ Nnr
type CnSet = record
v0, r0, q0: IntervalVector;
C0, B0, C,B : IntervalMatrix;
{va, ra, qa : IntervalVector}a∈Nn1 ∪...∪Nnr
The set ϕ(tk, [x0]) is represented as v0 + B0r0 + C0q0, the partial derivatives
Daϕ(tk, [x0]) are represented as va+Bra+Cqa. The matrices B,C are common
for all partial derivatives.
Notice, that if we start the Cr computation with an initial condition (15) then
there is no Lipschitz part at the beginning for the partial derivatives. Hence,
the initial values for C and B are set to the identity matrix and the initial values
for qa, ra are set to zero.
If the interval vectors ra become ’thick’ (i.e. theirs diameters are larger than
some threshold value) we can set a new Lipschitz part in our representation (it
must be done simultaneously for all Daϕ) and reset ra in the following way
qa = ra + (B
TC)qa, for a ∈ Nn1 ∪ . . . ∪ Nnr
ra = 0, for a ∈ Nn1 ∪ . . . ∪ Nnr
C = B
B = Id
A similar change of the Lipshitz part may be done when vectors ra become thick
in comparison to qa.
7 Derivatives of Poincaré map
Consider a differential equation
x′ = f(x), x ∈ Rn, f ∈ CK+1 (36)
Let ϕ : R × Rn → Rn be a (local) dynamical system induced by (36). Let
α : Rn → R be C1-map. Put Π = {x | α(x) = C}.
Definition 11 We will say that Π is a local section for the vector field f at
y0 ∈ Π if
〈∇α(y0)|f(y0)〉 6= 0. (37)
Assume x0 ∈ Rn and t0 ∈ R are such that Π is a local section at ϕ(t0, x0).
Consider an implicit equation
α(ϕ(tP (x), x)) = C. (38)
It follows easily from (37) and from the implicit function theorem that there
exists a uniquely defined tP : R
n−→◦ R in a neighborhood of x0, such that
tP (x0) = t0. The function tP is as smooth as the flow ϕ. We will refer to
tP as to the Poincare return time to section Π.
We define a Poincaré map P : Rn ⊃ dom (tP ) → Rn by
P (x) = ϕ(tP (x), x). (39)
Usually the Poincaré map is defined as a map P : Π1−→◦ Π2, where Π1,Π2 are
local sections in Rn. The approach taken here, i.e. treating the Poincaré map
as map P : Rn−→◦ Rn allows us to not to worry about the coordinates on local
section.
In this section we are interested in the partial derivatives of P defined by
(39).
From (39) we can compute ∂Pi
and we obtain
(x) = fi(P (x))
(x) +
(tP (x), x). (40)
We need ∂tP
. We differentiate (38) to obtain
(P (x))
fk(P (x))
(x) +
(tP (x), x)
(∇α(P (x)) · f(P (x))) ∂tP
(x) +
(P (x))
(tP (x), x) = 0. (41)
Hence
(x) = − 1〈∇α(P (x))|f(P (x))〉
(P (x))
(tP (x), x). (42)
7.1 Higher order derivatives of the Poincaré map
To make formulas transparent we will drop arguments of functions in this sec-
tion, but reader should be aware that for tP and its partial derivatives the
argument is x, for ϕ and Daϕ the argument is always the pair (tP (x), x).
From (40) we obtain
D(j,c)P =
ϕD(j)tPD(c)tP +
D(c)ϕD(j)tP +
ϕD(j,c)tP
D(j)ϕD(c)tP +D(j,c)ϕ.
It is easy to see that partial derivatives of high order give rise to quite complex
expressions and it is not entirely obvious how to organize it in some coherent
and programmable way. For this purpose we use the following
Lemma 12 For a multipointer a ∈ Nnp we have
DaP = Daϕ+
(δ1,...,δk)∈Np(k)
j=1Daδj tP
(δ1,...,δk)∈Np(k)
∂tk−1
Daδsϕ
j 6=sDaδj tP
Proof: By induction on p. For p = 1 formula (43) is equivalent to (40), because
the two last sums are taken over empty set. Assume (43) holds true for some
p ≥ 1 and fix a ∈ Nnp+1. Our goal is to show that
DaP = R1 +R2 +R3
where
R1 = Daϕ+
ϕDatP
(δ1,...,δk)∈Np+1(k)
Daδj tP
(δ1,...,δk)∈Np+1(k)
∂tk−1
Daδsϕ
j 6=s
Daδj tP
Write a = β + γ, where β ∈ Nnp and γ = (ap+1) ∈ Nn1 . From the induction
assumption we have
DaP = Dγ
ϕDβtP
(δ1,...,δk)∈Np(k)
j=1Dβδj tP
(δ1,...,δk)∈Np(k)
∂tk−1
Dβδsϕ
j 6=sDβδj tP
i=1 Si
where
S1 = Daϕ+
ϕDatP
DβϕDγtP
ϕDβtPDγtP
DγϕDβtP
(δ1,...,δk)∈Np(k)
j=1Dβδj tP
∂tk+1
ϕDγtP
(δ1,...,δk)∈Np(k)
j=1Dβδj tP
(δ1,...,δk)∈Np(k)
s=1Dβδs+γtP
j 6=s
Dβδj tP
(δ1,...,δk)∈Np(k)
∂tk−1
Dβδs+γϕ
j 6=sDβδj tP
(δ1,...,δk)∈Np(k)
DβδsϕDγtP
j 6=sDβδj tP
S10 =
(δ1,...,δk)∈Np(k)
r 6=s
∂tk−1
DβδsϕDβδr+γtP
j 6=s
j 6=r
Dβδj tP
Obviously R1 = S1. We will show that R2 = S3 + S6 + S7 and R3 = S2 + S4 +
S5 + S8 + S9 + S10.
Denote by Ri,k, i = 2, 3 a part of sum Ri with fixed k = 2, . . . , p + 1.
Similarly, let us denote by Si,k a part of sum Si, i = 5, . . . , 10, for k = 2, . . . , p.
Using decomposition of N p+1(2) as in (8) we obtain that R2,2 = S3 + S7,2.
Similarly, using (8) we observe that R2,k = S6,k−1 + S7,k for k = 3, . . . , p.
Finally, since N p+1(p + 1) = {((1), (2), . . . , (p + 1))} and γ = (ap+1) we find
that R2,p+1 = S6,p. This shows that R2 = S3 + S6 + S7.
It remains to show that R3 = S2 +S4 +S5 + S8 + S9 + S10. We will classify
possible terms by the fact, where p+ 1 appears in δi, i = 1, . . . , k and how this
δi enters in R3 as δs or δj . There are four cases
1. δs = (p+ 1)
2. δj = (p+ 1)
3. p+ 1 ∈ δs, |δs| ≥ 2
4. p+ 1 ∈ δj, |δj | ≥ 2
Let us fix k = 2. Let (δ1, δ2) ∈ N p+1(2). The term for case 1 is S4, for case 2
is S2, case 3 is S8,2 and case 4 is S10,2. Hence, R3,2 = S2 + S4 + S8,2 + S10,2.
For k = 3, . . . , p and fixed (δ1, . . . , δk) ∈ N p+1(k) we have: case 1 is given
by S5,k−1, case 2 by S9,k−1, case 3 by S8,k and case 4 by S10,k Hence, for
k = 3, . . . , p we have R3,k = S5,k−1 + S9,k−1 + S8,k + S10,k.
Finally, for k = p+ 1 we observe, that R3,p+1 = S5,p + S9,p. Indeed, in this
case (δ1, . . . , δp+1) = ((1), (2), . . . , (p + 1)). Hence, either for δs = γ we have
term S5,p and δs 6= γ we have S9,p.
We have showed that R3 = S2 + S4 + S5 + S8 + S9 + S10 and the proof is
finished.
Hence, if we know all the partial derivatives of tP up order p we can compute
the partial derivatives of the Poincaré map up the same order. In next subsection
we show how to compute partial derivatives of tP for affine sections.
7.2 Partial derivatives of tP for affine sections
Assume α : Rn → R is an affine map given by
α(x) = α0 +
αixi.
This is a quite restrictive assumption about sections, but it leads to relatively
simple formulas for DatP and it is sufficient for the applications we have in
mind.
Lemma 13 For a multipointer a ∈ Nnp holds
−DatP
∇α| ∂
= 〈∇α|Daϕ〉
∇α| ∂
(δ1,...,δk)∈Np(k)
j=1Daδj tP
(δ1,...,δk)∈Np(k)
∇α| ∂
∂tk−1
Daδsϕ
j 6=sDaδj tP
Proof: The proof is a direct consequence of Lemma 12 and (38). Since α is
affine, by differentiating of α(P (x)) = C we get 〈∇α|DaP 〉 = 0. Using formula
(43) for DaP we obtain our assertion.
Fix [x] ⊂ Rn and assume we have a rigorous bound for tP ([x]) ∈ [t1, t2] (see
[Z, Section 6] for more details on this). Lemmas 13 and 12 show that given rigo-
rous bounds for the partial derivatives Daϕ([t1, t2], [x]) and
Daϕ([t1, t2], [x])
up to some order p we can compute recursively rigorous bounds for the partial
derivatives of tP ([x]) and P ([x]) up to the same order. Notice, that
Daϕ are
given by Taylor coefficients of the solution of (14) with initial conditions P ([x])
for C0 part and Daϕ(tP (x), [x]) for equations for variations. Hence, these coef-
ficients can be easily computed using the automatic differentiation algorithm.
8 Applications.
One of the typical invariant sets in hamiltonian mechanics are invariant tori.
However, the existence of invariant torus in a given system is often difficult to
prove despite the fact that the theory is quite well developed. Probably the
best work in this direction was done by Celletti and Chercia [CC1, CC2], where
the an effective application (computer assisted proof) of KAM theory to the
restricted three body problem modelling system consisting of Sun, Jupiter and
asteroid 12 Victoria was given. Our aim here is more modest as we focus on the
invariant tori emanating from the elliptic fixed point satisfying suitable twist
condition.
In this section we show that the rigorous computations of partial derivatives
of a dynamical system up to order 3 or 5 can be used to prove that in a particular
system an invariant torus exists around some elliptic periodic orbits. In this
section this will be done for the forced pendulum equation and the Michelson
system.
8.1 Area preserving maps on the plane, normal forms and
KAM theorem
Definition 14 Let f : R2 → R2 be a smooth area preserving map, such that
f(p) = p. Let λ and µ be eigenvalues of df(p). Following [SM] we will call the
point p
• hyperbolic if λ, µ ∈ R and λ 6= µ,
• elliptic if λ = µ and λ 6= µ,
• parabolic if λ = µ.
The following KAM theorem will be the main tool to prove the existence of
invariant tori in this paper.
Theorem 15 [SM, §32] Consider an analytic area preserving map f : R2 →
2, f(r, s) = (r1, s1) where
r1 = r cosα− s sinα+O2l+2
s1 = r sinα+ s cosα+O2l+2 (44)
r2 + s2
and O2l+2 denotes convergent power series in r, s with terms of order greater
than 2l+ 1, only.
If at least one of γ1, . . . , γl is not zero then the origin is a stable fixed point
for map f . Moreover, in any neighborhood U of point 0 there exists an invariant
curve for map f around the origin contained in U .
The next theorem and its proof tells how to bring a planar area preserving
map in the neighborhood of an elliptic fixed point into the form (44).
Theorem 16 [SM, §23] Consider an analytic area preserving map f : R2 → R2
such that f(0) = 0. Let λ, λ̄ be complex eigenvalues of Df(0), such that |λ| =
|λ̄| = 1. If λk 6= 1 for k = 1, . . . , 2l+2, then there is an analytic area preserving
substitution such that in the new coordinates mapping f has form (44).
The proof of the above theorem is constructive, i.e. given the power series
for f at an elliptic fixed point one can construct explicitly an area preserving
substitution and compute the coefficients γ0, . . . , γl in (44). An explicit formula
for the coefficient γ1 in the above normal form is given in Appendix A.
8.2 The existence of invariant tori in forced pendulum.
Consider an equation
θ̈ = − sin(θ) + sin(ωt) (45)
Observe that (45) is hamiltonian.
Let us denote by Pω : R
2−→◦ R2 the Poincaré map for Equation (45) with
a parameter ω, i.e. Pω = ϕ(2π/ω, ·), where ϕ : R × R2−→◦ R2 is a local flow
induced by (45). Observe that (45) is nonautonomous, but it is equivalent to
first order system of autonomous ODE given by
= − sin(θ) + sin(ωt) (46)
In the sequel all rigorous computations for (45) will be in fact performed for the
system (46).
Observe that to any invariant closed curve for Pω corresponds and invariant
2-torus for (45).
Consider a set of parameter values
Ω1 = [2, 2.994], Ω2 = [3, 3.997], Ω3 = [4, 8]
Ω = Ω1 ∪ Ω2 ∪ Ω3
The following lemma was proved with computer assistance
Lemma 17 For all parameter values ω ∈ Ω there exists an elliptic fixed point
xω ∈ R2 for Pω. Moreover, there exists an area-preserving substitution such
that in the new coordinates the map fω(x) = Pω(x+xω)−xω has the form (44)
with l = 1 and γ1 6= 0.
Before we give the proof, let us briefly comment about the choice of the parame-
ter set Ω. For parameter values slightly lower than 2 we observe the parabolic
case, i.e. there exists a parameter value ω1 for which eigenvalues of the deriva-
tive of Pω1 are equal to −1. In two gaps in Ω below 3 and 4 we have resonances
of low order. Namely, we have parameter values with an elliptic fixed with ei-
genvalues to e±2π/3 = −1
i and e±iπ/2 = ±i, respectively. Clearly, in a
computer assisted proof we need to exclude a small interval around those para-
meters. For ω > 4 it seems that the interval Ω3 can be extended much further
to the right without any difficulty.
Proof of Lemma 17: A computer assisted proof consists of the following steps.
We cover the set Ω by 9910 nonequal subintervals ωi. Diameters of ωi’s were
relatively large for values far away from the parabolic cases and very small close
to them. For a fixed subinterval ωi we proceed as follows
1. Let ω̄ denote an approximate center of the interval ωi. We find an approxi-
mate fixed point for Pω̄ using the standard nonrigorous Newton method.
Let us denote such a point by xi.
2. We define a box centered at xi, i.e we set vi := xi + [−εi, εi]2, where
εi > 0 depends on subinterval ωi - the values we used are from the interval
[5 · 10−6, 3 · 10−3], depending on whether xi close to parameter values
corresponding to parabolic cases.
3. Using the C1-Lohner algorithm we compute the Interval Newton operator
[Mo, N, A] Ni := N(Pωi − Id, xi, vi) and verify that Ni ⊂ int vi. This
proves that for all ω ∈ ωi there exists a unique fixed point xω ∈ Ni for
4. Using the C3-Lohner algorithm we compute a rigorous bound for Pωi(Ni)
and DαPωi(Ni), α ∈ N21 ∪N22 ∪N23. Hence, we obtain a rigorous bound for
the coefficients in
fω(x) =
|α|=1
DaP (xω)x
5. We show that an arbitrary matrix M ∈ DPωi(Ni) has a pair of complex
eigenvalues λ, λ̄ which satisfy λk 6= 1 for k = 1, . . . , 4. From Theorem 16
it follows there exists an area-preserving substitution such that in the new
coordinates the map fω for ω ∈ ωi has the form (44) with l = 1.
6. We compute a rigorous bound for γ0 and γ1 which appear in the formula
(44) and verify that for ω ∈ ωi holds γ1 6= 0.
The rigorous bounds for the values of γ1 on Ω are
γ1(Ω1) ⊂ [0.29930416771330087, 30.118260918229566]
γ1(Ω2) ⊂ [0.099747909112924596, 0.56550301088840627]
γ1(Ω3) ⊂ [0.18574835001593507, 0.4129279974577012]
A computer assisted proof of the above took approximately 95 minutes on the
Pentium IV 3GHz processor.
As a straightforward consequence of Lemma 17 and Theorem 15 we obtain
Theorem 18 For all parameter values ω ∈ Ω there exists an elliptic fixed point
xω ∈ R2 for Pω. Moreover, any neighborhood of point xω contains an invariant
curve for Pω around xω.
8.3 Higher order normal forms.
In the previous section it was shown that C3 computations are sufficient to
prove that for (45) a family of invariant tori exists. However, it may happen
that the coefficient γ1 in the normal form vanishes. In this situation we may try
to compute higher order normal form. As an example we consider a pendulum
with a different forcing term,
θ̈ = − sin(θ) + sin(ωt) + sin(2ωt). (47)
Theorem 19 Let Pω be the Poincaré map for (47). For all parameter values
ω ∈ Ω∗ = [2.9957694795, 2.9957694796] there exists an elliptic fixed point xω ∈
2 for Pω. Moreover, any neighbourhood of point xω contains an invariant
curve for Pω around xω.
Proof: The main concept of the proof is the same as in Lemma 17. Using the
nonrigorous Newton method we find an approximate fixed point
x = (−7.7491573604896152 · 10−12,−0.54723831527031352).
We set v = x + 3 · 10−5([−1, 1] × [−1, 1]). Using the C1-Lohner algorithm we
compute the Interval Newton Operator of Pω − Id on v and we obtain that for
all ω ∈ Ω∗, N = N(Pω − Id, center(v), v) ⊂ (N1, N2), where
N1 = [−5.1582932672798325, 5.1582631625020222] · 10−10
N2 = [−0.54723831580217108,−0.54723831470891193]
Since N ⊂ v we conclude that for all ω ∈ Ω∗ there exists a unique fixed point
xω ∈ N for the Poincaré map.
Using C5-Lohner algorithm we compute a rigorous bound for PΩ∗(N) and
DαPΩ∗(N), α ∈ N21 ∪ . . . ∪ N25. Hence, we obtain a rigorous bound for the
coefficients in
fω(x) =
|α|=1
DαP (xω)x
α +O6
We show that an arbitrary matrix M ∈ DPΩ∗(N) has a pair of complex ei-
genvalues λ, λ̄ which satisfy λk 6= 1 for k = 1, . . . , 6. From Theorem 16 it follows
there exists an area-preserving substitution such that in the new coordinates the
map fω for ω ∈ Ω∗ has the form (44) with l = 2.
Next, we compute a rigorous bound for γ1 and γ2 which appear in the formula
(44) and we get
γ1(Ω∗) ⊂ [−5.3924276719042241, 5.381714805052106] · 10−6
γ2(Ω∗) ⊂ [199.95180660157078, 199.99104965939162]
Since for ω ∈ Ω∗, γ2(ω) 6= 0 the assertion follows from Theorem 15.
The main observation which makes this example interesting is that there
exists ω∗ ∈ Ω∗ for which γ1(ω∗) = 0 and we cannot conclude the existence of
invariant tori for all ω ∈ Ω∗ from C3 computations. To be more precise, we
computed the coefficient γ1 for the parameter values ω1 = minΩ∗ and ω2 =
maxΩ∗ and we get
γ1(ω1) ∈ [−2.3559594437885885,−1.3593457220363871] · 10−8
γ1(ω2) ∈ [2.9671154858524365 · 10−9, 1.2819312939263052 · 10−8]
Since γ1 exists for all ω ∈ Ω∗ and depends continuously on ω we conclude, that
γ1(ω∗) = 0 for some ω∗ ∈ Ω∗.
8.4 Application to the Michelson system
The existence of an invariant curve for a planar map f : R2 → R2 can be proven
without assumption that f is measure preserving. The key assumption in the
proof given in [SM] is that any curve γ around an elliptic point intersect its
image under f , i.e. f(γ)∩ γ 6= ∅. Such a situation is also observed in reversible
planar map around an symmetric elliptic fixed points.
Definition 20 An invertible transformation M : Ω −→ Ω is called a reversing
symmetry of a local dynamical system φ : T × Ω −→ Ω, T = R or T = Z if the
following conditions are satisfied
1. if (t, x) ∈ dom (φ) then (−t, S(x)) ∈ dom (φ).
2. S(φ(t, x)) = φ(−t, S(x)))
Remark 21 In the discrete time case, the above two conditions are equivalent
to identity
M ◦ f = f−1 ◦M.
where f = φ(1, ·) is a generator of φ.
Definition 22 Let φ : T×Ω → Ω be a local (discrete or continuous) dynamical
system. For x ∈ Ω put
I(x) = {t ∈ T : (t, x) ∈ dom (φ)}
O (x) = {φ(t, x) ∈ Ω : t ∈ I(x)}
The set O (x) will be called a trajectory of a point x.
Definition 23 Assume S is an reversing symmetry for φ : T × Ω → Ω. An
orbit O (x) is called S-symmetric orbit if O (x) = S(x).
Remark 24 [La] In continuous case the orbit O (x) is S-symmetric if it contains
a point from the set Fix(S) = {y : S(y) = y}.
Remark 25 [Wi2, Lem.3.3] It is easy to see that if Θ ⊂ Ω is a Poincaré section
for a R-reversible flow φ : R × Ω → Ω such that Θ = R(Θ) then the Poincaré
map P : Θ → Θ is R|Θ-reversible.
As we observed at the beginning of this section, an R-reversible planar map
may admit an invariant curve around an R-symmetric elliptic fixed point. In
reversible case a planar map admits the same normal form around symmetric, el-
liptic fixed point as in the area-preserving case and the substitution which tends
the map to the normal form is exactly the same as we described in Appendix A
– for details see [Se, BHS].
Consider an ODE 
ẋ = y
ẏ = z
ż = c2 − y − 1
On one hand, the system (48) is an equation for the steady state solution of
one-dimensional Kuramoto-Sivashinsky PDE and it is known in the literature
as the Michelson system[Mi]. On the other hand, this system appears as a part
of the limit family of the unfolding of the nilpotent singularity of codimension
three (see [DIK1]).
The system (48) is reversible with respect to the symmetry
R : (x, y, z, t) → (−x, y,−z,−t) (49)
and since the divergence vanishes it is also volume preserving.
A dynamical system induced by (48) exhibits several types of dynamics for
different values of parameter. For sufficiently large c there is a simple invariant
set consisting of two equilibria (±c
2, 0, 0) and heteroclinic orbit between them
[MC]. Lau [Lau] numerically observed that when the parameter c decreases a
cascade of cocoon bifurcations occurs and at the limit value c ≈ 1.266232337
a periodic orbit is born through a saddle-node bifurcation. This hypothesis
has been proved in [KWZ]. The computer assisted proof of this fact given in
[KWZ] uses the algorithm presented in this paper in order to compute partial
derivatives up to second order for a certain Poincaré map.
For the parameter value equal to one and slightly smaller than one it was
proven in [DIK2, Wi1, Wi2, Wi3] that the system has rich and complicated
dynamics including symbolic dynamics, heteroclinic solutions, Shilnikov homo-
clinic solutions.
However, as the bifurcations diagram presented by Michelson suggests [Mi,
Fig.1] for all parameter values c ∈ (0, 0.3195) there are at least two elliptic
periodic orbits with large invariant islands around them. In this section we
present a proof that such islands exist for some range of parameter values. The
main idea of the proof is almost the same as in the previous section. There are
two main differences. First, the Poincaré map will not be a time shift. Therefore
computations of the partial derivatives of the Poincaré map require Lemma 12
and Lemma 13. Second difference is: we use the shooting method instead of
the interval Newton method for the proof of the existence of symmetric periodic
orbit.
The aim of this section is to prove the following
Theorem 26 For all parameter values from the set
C = C1 ∪ C2 = [0.1, 0.225]∪ [0.226, 0.25]
there exists a symmetric elliptic periodic orbit for the Michelson system (48).
Moreover, each neighbourhood of such an orbit contains a 2D tori invariant
under the flow generated by the Michelson system.
Let us define the Poincaré section Π := {(0, z, y) : z, y ∈ R}. Let Pc =
(P1, P2) : Π−→◦ Π be the Poincaré map for the system with the parameter value
c. Notice, that Pc is in fact a half Poincaré map, which means that the trajectory
of x crosses Π in opposite directions when passing through x and Pc(x), and
therefore periodic orbits for the Michelson system corresponds to periodic points
for P 2.
Since the section Π is invariant under symmetry (x, y, z) → (−x, y,−z),
from Remark 25 the Poincaré map is also reversible with respect to an involu-
tion R(y, z) = (y,−z). We will use the same letter R to denote the reversing
symmetry of the Poincaré map and the Michelson system.
Let us comment about the choice of the set C. In the gap between intervals
C1 and C2 there is a parameter value c∗ for which the eigenvalues of the Poincaré
map P 2c∗ are ±i. Apparently at this parameter value we have a bifurcation and
four periodic islands are born as it is shown in Fig.1 - see also a movie mpp.mov
available at [Wi4] which presents an animation of the phase portrait of Pc for
the parameter values from the range [0.1, 0.25].
Proof of Theorem 26. The main concept of the proof is quite similar to the
one presented in Lemma 17. We divide the set C of parameter values onto 20800
nonequal parts (smaller when close to the bifurcation parameter c∗ and close to
0.1 and 0.25). For a fixed subinterval ci from the grid we proceed as follows
1. Let c̄ denote a center of the interval ci. We find an approximate fixed point
of P 2c̄ using the standard nonrigorous Newton method. Let us denote this
point by (yi, zi).
2. Since the map Pc is reversible one can prove the existence of the fixed
point for P 2c using the shooting method as follows.
Let Fix(R) = {(y, z) ∈ Π : R(y, z) = (y, z)} = {(y, 0) ∈ Π : y ∈ R}.
Since Pc satisfies (Pc ◦ R)2 = Id whenever the left side is defined, one
can see that if x ∈ Fix(R) and Pc(x) ∈ Fix(R) then P 2c (x) = x. Let us
remark, that we always get an approximate fixed points (yi, zi) resulting
from the nonrigorous Newton method very close to Fix(R). We define
two points u1 = (yi − εi, 0), u2 = (yi + εi, 0) ∈ Fix(R), where εi is a small
number depending on ci and we show that πz(Pci(u1)) · πz(Pci(u2)) < 0,
where πz is a projection onto z coordinate. Hence, if the Pci is defined
on the set Ni = (0, [yi − εi, yi + εi], 0) then for all parameter values c ∈ ci
there is a point uc ∈ N which satisfies πz(Pc(uc)) = 0 and therefore
Figure 1: Phase portrait of the Poincaré map Pc (top) before bifurcation for
c = 0.225 and (bottom) after bifurcation for c = 0.226 with four periodic islands.
Between those parameters resonant case occurs with eigenvalues equal to ±i.
See also auxiliary material [Wi4].
Pc(uc) ∈ Fix(R). This shows that for all c ∈ ci there exists a fixed point
for P 2c inside Ni provided Pc is defined on Ni, which will be discussed
below.
3. Using C3-Lohner algorithm we compute rigorous bounds for P 2ci(Ni) and
DαP 2ci(Ni) for α ∈ N
1 ∪ N22 ∪ N23. This implies also that Ni ⊂ domPci .
4. We show that an arbitrary matrix M ∈ DPci(Ni) has a pair of complex
eigenvalues λ, λ̄ which satisfy λk 6= 1 for k = 1, . . . , 4. From Theorem 16
it follows there exists an area-preserving substitution such that in the new
coordinates the map Pc for c ∈ ci has the form (44) with l = 1.
5. We compute a rigorous bound for γ0 and γ1 which appear in the formula
(44) and verify that for c ∈ ci holds γ1 6= 0.
The rigorous bounds for the values of γ1 on C are
γ1(C1) ⊂ [0.014515898754816965, 157.76639522562903]
γ1(C2) ⊂ [1.1002393483255526, 151.35147664498677]
The computer assisted proof of the above took approximately 7 hours and 50
minutes on the Pentium IV 3GHz processor.
9 Implementation notes.
All the algorithms presented in this paper have been implemented in C++ by
authors and are part of the CAPD library [CAPD]. In particular, the package
implements the computation of partial derivatives of a flow with respect to
initial condition, partial derivatives of Poincaré maps for linear sections and
computations of normal forms for planar maps up to order 5.
The implementation combines the automatic and symbolic differentiation in
order to generate a coefficients in Taylor series for the solutions of the system
(14).
Our tests shows that without difficulty we can compute partial derivatives
up to order 3 for an equation in 8-dimensional phase space (which gives 1320
equations to solve) on a computer with 512MB memory. However, our current
implementation is optimized for lower dimensional problems. All the trees which
represent formulas (14) are stored in the memory of a computer. This speeds
up computations because we do not need to recompute all the multiindices,
multipointers and submultipointers in each step of the algorithm. Unfortunately,
such an implementation is memory-consuming. Therefore, higher dimensional
problems require a computer with huge memory even for C3 or C5 computations.
A Explicit formulas for third order normal forms
for a planar map
The goal of this section is to give some details about the proof of Theorem 16.
We want to present some formulas to give the reader the feeling about the
necessary computations.
Throughout this section we assume that the assumptions of Theorem 16 are
satisfied. In the neighbourhood of 0 f is given by a real, convergent power series
f(x, y) = (x1, y1)
ak−l,lx
lyk−l
bk−l,lx
lyk−l
Denote also by f : C2 → C2 a complex extension of f . Let λ, λ̄ ∈ C be
complex eigenvalues of Df(0) and v, v̄ ∈ C corresponding eigenvectors (here
bar denotes the complex conjugation). Then, using a linear substitution of the
form L = [vT , v̄T ], we can change the coordinate system such that in the new
coordinates the mapping f has the form
f(ξ, η) = (λξ + p(ξ, η), λ̄η + q(ξ, η))
p(ξ, η) =
pl,k−lξ
lηk−l
q(ξ, η) =
ql,k−lξ
lηk−l
pi,j = qj,i for i, j ≥ 0.
The last condition is a consequence of the invariance of R2 ⊂ C2 under the
complex map f . We will refer to it as the reality condition. Namely, the set
2 ⊂ C2 in the new coordinates (ξ, η) is given by ξ = η and the condition
f(R2) ⊂ R2 expressed in coordinates (ξ, η) is equivalent to (50).
Assume now, that λk 6= 1 for k = 1, . . . , 4. Then an analytic area-preserving
substitution satisfying reality condition (50)
(Φ(z, v),Ψ(z, v)) = (z1, v1)
z1 = z +
φl,k−lz
lvk−l + · · ·
v1 = v +
ψl,k−lz
lvk−l + · · ·
where
ψ2,0 = φ0,2 = −λ2p0,2(λ3 − 1)−1
ψ1,1 = φ1,1 = −p1,1(λ− 1)−1
ψ0,2 = φ2,0 = p2,0(λ
2 − λ)−1
ψ3,0 = φ0,3 = −λ3 (p0,3 + p1,1φ0,2 + 2q0,2ψ0,2) (λ4 − 1)−1
ψ2,1 = φ1,2 =
λ2 − 1 (p1,2 + 2p2,0φ0,2 + p1,1φ1,1 + p1,1ψ0,2 + 2p0,2ψ1,1)
ψ1,2 = φ2,1 = −φ2,0ψ0,2 + φ0,2ψ2,0
ψ0,3 = φ3,0 = (p3,0 + 2p2,0φ2,0 + p1,1ψ2,0) (λ
3 − λ)−1
brings f = (f1, f2) to the normal form
(z, v) → (z(α0 + α2zv), v(β0 + β2zv)) +O((zv)2)
β0 = α0 = λ
β2 = α2 = q1,2 + 2q2,0φ0,2 + q1,1φ1,1 + q1,1ψ0,2 + 2q0,2ψ1,1
Finally, let γ0 ∈ R be such that λ = α0 = eiγ0 and we compute coefficient γ1 by
From the proof given in [SM] it follows that γ1 ∈ R and the mapping f in
coordinates (z, v) has the form
f(z, v) =
zei(γ0+γ1zv), ve−i(γ0+γ1zv)
where O4 is a convergent power series with the terms of degree at least 4.
Again, the coefficients of f(z, v) satisfy reality condition (50). In order to
express this normal form in terms of real variables we make a linear substitution
z = r + is, v = r − is
and we obtain the normal form for f
f(r, s) = (r1, s1) +O4
r1 = r cos(γ0 + γ1(r
2 + s2))− s sin(γ0 + γ1(r2 + s2))
s1 = r sin(γ0 + γ1(r
2 + s2)) + s cos(γ0 + γ1(r
2 + s2))
which agrees with (44).
The formulas for higher order terms φi,j , ψi,j (and for γ2, which are not given
here) has been computed in Mathematica.
References
[A] G. Alefeld, Inclusion methods for systems of nonlinear equations - the
interval Newton method and modifications, in Topics in Validated Com-
putations, J. Herzberger (Editor), Elsevier Science B.V., 1994, pages
[BHS] H.W. Broer, G.B. Huitema and M.B. Sevryuk, Quasi-periodicity in
families of dynamical systems: order amidst chaos, Lecture Notes in
Mathematics, Vol. 1645, Springer Verlag, (1996).
[BM] M. Berz, K. Makino, New Methods for High-Dimensional Verified Qua-
drature, Reliable Computing, 5, 13-22 (1999)
[CAPD] CAPD – Computer Assisted Proofs in Dynamics group, a C++ pa-
ckage for rigorous numerics, http://capd.wsb-nlu.edu.pl.
[CC1] A. Celletti, L. Cherchia, KAM Stability for three-body problem of the
Solar System, Z. angw. Math. Phys. 57 (2006) 33-41
[CC2] A. Celletti, L. Cherchia, KAM Stability and Celestial Mechanics, Me-
moirs of the AMS, Vol 187, Num 878 (2007)
[DIK1] F. Dumortier, S. Ibáñez, and H. Kokubu, New aspects in the unfolding
of the nilpotent singularity of codimension three, Dynam. Syst. 16
(2001), 63–95.
[DIK2] F. Dumortier, S. Ibáñez, and H. Kokubu, Cocoon bifurcation in three
dimensional reversible vector fields, Nonlinearity 19 (2006), 305–328.
[GZ] Z. Galias, P. Zgliczyński, Computer assisted proof of chaos in the Lo-
renz system, Physica D, 115, 1998,165–188
[HNW] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential
Equations I, Nonstiff Problems, Springer-Verlag, Berlin Heidelberg
1987.
[HZHT] B. Hassard, J. Zhang, S. Hastings, W. Troy, A computer proof that
the Lorenz equations have ”chaotic” solutions, Appl. Math. Letters, 7
(1994), 79–83
[KZ] T. Kapela, P. Zgliczyński, The existence of simple choreographies for
N-body problem - a computer assisted proof, Nonlinearity, 16 (2003),
1899-1918
[KWZ] H. Kokubu, D.Wilczak, P. Zgliczyński, Rigorous verifi-
cation of cocoon bifurcations in the Michelson system,
http://www.ii.uj.edu.pl/~wilczak, submitted
[La] J.S.W. Lamb, Reversing symmetries in dynamical systems, PhD The-
sis, Amsterdam University, 1994
http://capd.wsb-nlu.edu.pl
http://www.ii.uj.edu.pl/~wilczak
[Lau] Y-T Lau, The “cocoon” bifurcations in three-dimensional systems with
two fixed points, Int. Jour. Bif. Chaos, Vol.2, No.3 (1992) 543-558.
[Lo] 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., Claren-
don Press, Oxford, 1992.
[MC] C.K. McCord, Uniqueness of connecting orbits in the equation Y (3) =
Y 2 − 1, J. Math. Anal. Appl. 114, 584-592.
[Mi] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation,
Physica D, 19, (1986) 89-111.
[Mo] R.E. Moore, Interval Analysis. Prentice Hall, Englewood Cliffs, N.J.,
[MM] K. Mischaikow, M. Mrozek, Chaos in the Lorenz equations: A compu-
ter assisted proof. Part II: Details, Mathematics of Computation, 67,
(1998), 1023–1046
[MZ] M. Mrozek, P. Zgliczyński, Set arithmetic and the enclosing problem
in dynamics, Annales Pol. Math., 2000, 237–259
[NJ] N. S. Nedialkov, K. R. Jackson, An Interval Hermite – Obreschkoff
Method for Computing Rigorous Bounds on the Solution of an Initial
Value Problem for an Ordinary Differential Equation, chapter in the
book Developments in Reliable Computing, editor T. Csendes, 289-310,
Kluwer, Dordrecht, Netherlands, 1999.
[N] A. Neumeier, Interval methods for systems of equations, Cambridge
University Press, 1990.
[RNS] T. Rage, A. Neumaier, C. Schlier, Rigorous verification of chaos in a
molecular model, Phys. Rev. E 50 (1994), 2682–2688
[Ra] L.B. Rall, Automatic Differentiation: Techniques and Applications,
volume 120 of Lecture Notes in Computer Science. Springer Verlag,
Berlin, 1981
[Se] Sevryuk, M. B. Reversible systems, Lecture Notes in Mathematics,
1211. Springer-Verlag, Berlin, 1986
[SM] C.L. Siegel, J.K. Moser, Lectures on Celestial Mechanics, Springer-
Verlag Berlin Heidelberg New York, 1995.
[SK] D. Stoffer, U. Kirchgraber, Possible chaotic motion of comets in the
Sun Jupiter system - an efficient computer-assisted approach, Nonli-
nearity, 17 (2004) 281-300.
[T] W. Tucker, A Rigorous ODE solver and Smale’s 14th Problem, Foun-
dations of Computational Mathematics, (2002), Vol. 2, Num. 1, 53-117
[W] W. Walter, Differential and integral inequalities, Springer-Verlag Ber-
lin Heidelberg New York, 1970
[Wi1] D. Wilczak, The existence of Shilnikov homoclinic orbits in the Mi-
chelson system: a computer assisted proof, Found. Comp. Math. Vol.6,
No.4, 495-535, (2006).
[Wi2] D. Wilczak, Chaos in the Kuramoto–Sivashinsky equations – a compu-
ter assisted proof, J. Diff. Eqns., 194, 433-459 (2003).
[Wi3] D. Wilczak, Symmetric heteroclinic connections in the Michelson sys-
tem – a computer assisted proof, SIAM J. App. Dyn. Sys., Vol.4, No.3,
489-514 (2005).
[Wi4] D. Wilczak, http://www.ii.uj.edu.pl/~wilczak, a refference for
auxiliary materials.
[WZ] D. Wilczak, P. Zgliczyński, Heteroclinic Connections between Periodic
Orbits in Planar Restricted Circular Three Body Problem - A Compu-
ter Assisted Proof, Commun. Math. Phys. 234, 37-75 (2003).
[Z1] P. Zgliczyński, Computer assisted proof of chaos in the Hénon map and
in the Rössler equations, Nonlinearity, 1997, Vol. 10, No. 1, 243–252
[Z] P. Zgliczyński, C1-Lohner algorithm, Foundations of Computational
Mathematics, (2002) 2:429–465
[ZPer] P. Zgliczyński, Lohner Algorithm for perturbations of ODEs and diffe-
rential inclusions, http://www.ii.uj.edu.pl/~zgliczyn
http://www.ii.uj.edu.pl/~wilczak
http://www.ii.uj.edu.pl/~zgliczyn
	Introduction
	Basic definitions
	Multiindices
	Multipointers
	Equations for variations
	Cr-Lohner algorithm
	Why one needs an Cr-algorithm?
	An outline of the algorithm
	Computation of a rough enclosure for Da
	The procedure for the computation of the rough enclosure for V.
	Computation of [Vk+1]
	Composition formulas
	The procedure for computation of [Vk+1]
	Rearrangement for Va - the evaluation of Equation (32)
	Derivatives of Poincaré map
	Higher order derivatives of the Poincaré map
	Partial derivatives of tP for affine sections
	Applications.
	Area preserving maps on the plane, normal forms and KAM theorem 
	The existence of invariant tori in forced pendulum.
	Higher order normal forms.
	Application to the Michelson system
	Implementation notes.
	Explicit formulas for third order normal forms for a planar map

  We present a Lohner type algorithm for the computation of rigorous bounds for
solutions of ordinary differential equations and its derivatives with respect
to initial conditions up to arbitrary order. As an application we prove the
existence of multiple invariant tori around some elliptic periodic orbits for
the pendulum equation with periodic forcing and for Michelson system.

Introduction
This paper is a sequel to [Z]. We present here a Lohner-type algorithm for
computation of rigorous enclosures of partial derivatives with respect to initial
conditions up to an arbitrary order r of the flow induced by an autonomous
ODE, hence the name Cr-Lohner algorithm. Let r be a positive integer, then
by Cr-algorithm we will mean the routine which gives rigorous estimates for
partial derivatives with respect to initial conditions up to an order r and Cr-
computations we mean an application of an Cr-algorithm.
Our main motivation for the development of Cr-algorithm was a desire to
provide a tool, which will considerably extend the possibilities of computer
assisted proofs in the dynamics of ODEs. Till now most of such proofs have used
topological conditions (see for example [HZHT, MM, GZ, Z1]) and additionally
conditions on the first derivatives with respect to initial conditions (see for
example [RNS, T, Wi1, WZ, KZ]), hence it required C0- and C1-computations,
respectively. The spectrum of problems treated includes the questions of the
existence of periodic orbits and their local uniqueness, the existence of symbolic
dynamics, the existence of hyperbolic invariants sets, the existence of homo-
and heteroclinic orbits. To treat other phenomena, like bifurcations of periodic
orbits, the route to chaos, invariant tori through KAM theory one needs the
knowledge of partial derivatives with respect to initial conditions of higher order.
In principle, one can think that a good rigorous ODE solver should be en-
ough. Namely, to compute the partial derivatives of the flow induced by
x′ = f(x), x ∈ Rn (1)
1 Research supported by an annual national scholarship for young scientists from the
Foundation for Polish Science
2 Research supported in part by Polish State Ministry of Science and Information Techno-
logy grant N201 024 31/2163
http://arxiv.org/abs/0704.0720v1
it is enough to rigorously integrate a system of variational equations obtained by
a formal differentiation of (1) with respect to the initial conditions. For example
for r = 2 we have the following system
x′ = f(x), (2)
Vij(t) =
(x)Vsj(t) (3)
Hijk(t) =
s,r=1
∂xs∂xr
(x)Vrk(t)Vsj(t) +
(x)Hsjk(x), (4)
with the initial conditions
x(0) = x0, V (0) = Id, Hijk(0) = 0, i, j, k = 1, . . . , n.
It is well known that if by ϕ(t, x0) we denote the (local) flow induced by (1),
(t, x0) = Vi,j(t),
∂xj∂xk
(t, x0) = Hijk(t).
Analogous statements are true for higher order partial derivatives with respect
to initial conditions.
It turns out that a straightforward application of a rigorous ODE solver to
the system of variational Equations (2–4) is very inefficient. Namely, it totally
ignores the structure of the system and leads to a very poor performance and
unnecessary long computation times (see Section 4.1).
Our algorithm is a modification of the Lohner algorithm [Lo], which takes
into account the structure of variational Equations (2–4). Basically it consists of
the Taylor method, a heuristic routine for a priori bounds for solution of (2–4)
during a time step and a Lohner-type control of the wrapping effect, which is
done separately for x and partial derivatives with respect initial conditions (the
variables V and H in (3,4)). The Taylor method is realized using the automa-
tic differentiation [Ra] and the algorithms for computation of compositions of
multivariate Taylor series.
The proposed algorithm has been successfully applied in [HNW] to the Mi-
chelson system [Mi], where a computer assisted proof of the existence of a cocoon
bifurcation was presented. Some parts of this proof required C2-computations.
In the present paper in Section 8 we show an application of our algorithm to
pendulum equation with periodic forcing and the Michelson system. We used
it to compute rigorous bounds for the coefficients of some normal forms up to
order five, which enabled us to prove the existence of invariant tori around some
elliptic periodic orbits in these systems using KAM theorem for twist maps on
the plane. These proofs required C3 and C5 computations.
2 Basic definitions
To effectively deal with the formulas involving partial derivatives we will use ex-
tensively a notation of multiindices, multipointers and submultipointers throu-
ghout the paper.
As an motivation let us consider the formula for the partial derivatives of
the composition of maps. Assume g : Rn → Rn and f : Rn → R are of class C3.
We have
∂3(f ◦ g)
∂xi∂xj∂xc
k,r,s=1
∂xk∂xr∂xs
∂xi∂xj∂xc
k,r=1
∂xk∂xr
∂xi∂xc
∂xj∂xc
∂xi∂xj
To the operator ∂
∂xi1∂xi2∂xi3
we can in a unique way assign a multipointer,
which is a nondecreasing sequence of integers (j1, j2, j3), such that {i1, i2, i3} =
{j1, j2, j3}. A submultipointer is a multipointer, which is a part of a longer mul-
tipointer, for example (i, j, c)(1,3) = (i, c). One observes, that submultipointers
appear at several places in the above formula.
A multiindex is an element of α ∈ Nn. It is another way to represent various
partial derivatives. The coefficient αi tells us how many times to differentiate
a function with respect to the i-th variable. Obviously, we have one-to-one
correspondence between multipointers and multiindices.
2.1 Multiindices
By N we will denote the set of nonnegative integers, i.e. N = {0, 1, 2, . . .}.
Definition 1 An element τ ∈ Nn will be called a multiindex.
For a sequence α = (α1, . . . , αn) ∈ Nn and a vector x = (x1, . . . , xn) ∈ Rn we
1. |α| = α1 + · · ·+ αn
2. α! = α1! · α2! · · ·αn!
3. xα = (xα11 , . . . , x
By eni ∈ Nn we will denote
eni = (0, 0, . . . , 0,
1 , 0, . . . , 0, 0).
We will drop the index n (the dimension) in the symbol eni when it is obvious
from the context.
Put Nnp := {a ∈ Nn : |a| = p}.
For δ = (δ1, . . . , δk) ∈ Nn1 × · · · × Nnk we set
1. |δ| =
i=1 |δi|
2. δ! =
i=1 δi!
Let f = (f1, . . . , fm) : R
n → Rm be sufficiently smooth. For α ∈ Nn we set
1. Dαfi =
∂|α|fi
∂xα11 · · · ∂x
2. Dαf = (Dαf1, D
αf2, . . . , D
For a function f : R× Rn → Rn by Dαfi(t, x) we will denote Dαfi(t, ·)(x) and
similarly
Dαf(t, x) = (Dαf1(t, x), . . . , D
αfn(t, x)).
This convention means that Dα always acts on x-variables.
2.2 Multipointers
For a fixed n > 0 and p > 0 we define
Nnp := {(a1, a2, . . . , ap) ∈ Np : 1 ≤ a1 ≤ · · · ≤ ap ≤ n}
N = Nn :=
Definition 2 An element of Nn will be called a multipointer.
Remark 3 A function
Λ : Nnp ∋ (a1, . . . , ap) →
enai ∈ N
p (5)
is a bijection.
Let f = (f1, . . . , fm) : R
n → Rm be a sufficiently smooth. For a ∈ Nnp we set
1. Dafi :=
∂xa1 . . . ∂xap
2. Daf := (Daf1, . . . , Dafm)
For a function f : R × Rn → Rn by Dafi(t, x) we will denote Dafi(t, ·)(x). In
the light of the above notations Dαf = D
Λ(α)f .
For a = (a1, a2, . . . , an) ∈ Nnp and b = (b1, b2, . . . , bn) ∈ Nnq we define
a+ b = (a1 + b1, . . . , an + bn) ∈ Nnp+q.
For α ∈ Nnp and β ∈ Nnq we define
α+ β = Λ−1 (Λ(α) + Λ(β)) ∈ Nnp+q.
By ≤ we will denote a linear order (lexicographical order) in N defined in
the following way. For a ∈ Nnp and b ∈ Nnq
(a ≤ b) ⇐⇒
either ∃i, i ≤ p, i ≤ q, ai < bi and aj = bj for j < i
or p ≤ q and ai = bi for i = 1, . . . , p.
Definition 4 For k ≤ p we set
N p(k) := {(δ1, . . . , δk) ∈ (N p)k : δ1 ≤ · · · ≤ δk, δ1+ · · ·+δk = (1, 2, . . . , p)} (7)
We will use N p(k) extensively in the next section. Its will be used to label
terms in Dαfi(ϕ(t, x)). Observe that for p > 0
N p(1) = {(1, 2, . . . , p)}
N p(p) = {((1), (2), . . . , (p))}
One can construct all elements of N p(k) using the following recursive procedure.
From the definition of N p(k) it follows that if (δ1, . . . , δm−1) ∈ N p−1(m − 1)
then (δ1, . . . , δm−1, (p)) ∈ N p(m) (notice that order is preserved). Similarly, if
(δ1, . . . , δm) ∈ N p−1(m) then
(δ1, . . . , δs−1, δs + (p), δs+1, . . . , δm) ∈ N p(m)
and again order of elements is preserved. Hence, for p > 2 and 1 < k < p we
have N p(k) = A ∪B where
(δ1, . . . , δk−1, (p)) : (δ1, . . . , δk−1) ∈ N p−1(k − 1)
(δ1, . . . , δs−1, δs + (p), δs+1, . . . , δk) : (δ1, . . . , δk) ∈ N p−1(k)
} (8)
and the sets A and B are disjoint.
Another way to generate all elements of N p(k) can be described as follows
• decompose the set {1, 2, . . . , p} into k nonempty and disjoints sets ∆i,
i = 1, . . . , k
• we sort each ∆i and permute ∆i’s to obtain min(∆1) < min(∆2) < · · · <
min(∆k)
• we define δi to be an ordered set consisting of all elements of ∆i for
i = 1, . . . , k
Definition 5 For an arbitrary a ∈ Nnp and δ ∈ N
k such that k ≤ p we define a
submultipointer aδ ∈ Nnk by (aδ)i = aδi for i = 1, . . . , k, which can be expressed
using Λ as follows
aδ := Λ
enaδi
∈ Nnk
3 Equations for variations
Consider an ODE x′ = f(x) where f is CK+1. Let ϕ : R× Rn−→◦ Rn be a local
dynamical system induced by x′ = f(x). It is well known, that ϕ ∈ CK and one
can derive the equations for partial derivatives of ϕ by differentiating equation
(t, x) = f(ϕ(t, x)) with respect to the initial condition x. As a result we
obtain a system of so-called equations for variations, whose size depends on the
order r of partial derivatives we intend to compute. An example of such system
for r = 2 is given by (2–4) with initial conditions given by (5).
The goal of this section is to write the equations for variations in a compact
form using multipointers and multiindices, which allows us to take into account
the symmetries of partial derivatives,
Lemma 6 Assume f ∈ Cr+1 and let ϕ : R × Rn−→◦ Rn be a local dynamical
system induced by x′ = f(x). Then for a ∈ Nnp such that p ≤ r holds
Daϕi =
i1,...,ik=1
Dei1+···+eik fi
(δ1,...,δk)∈Np(k)
Daδjϕij (9)
for i = 1, . . . , n.
Proof: In the proof the functions Dei1+···+eik fi are always evaluated at
ϕ(t, x), and various partial derivatives of ϕ are always evaluated at (t, x), there-
fore the arguments will be always dropped to simplify formulae. We prove the
lemma by induction on p = |a|. If p = 1 then a = (c) for some c ∈ {1, . . . , n}
and (9) becomes
D(c)ϕi =
Desfi ·D(c)ϕs.
Assume (9) holds true for p − 1, p > 1. Let us fix a ∈ Nnp . We have
a = b + (c), where b = (a1, . . . , ap−1) ∈ Nnp−1 and c = ap. Since (9) is satisfied
for p− 1, therefore we have
Daϕi = D(c)
= D(c)
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδjϕij
i1,...,ik+1=1
β:=ei1+···+eik+1
Dβfi ·D(c)ϕik+1
(δ1,...,δk)∈Np−1(k)
Dbδjϕij
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
For k = 1, . . . , p we set
Tk :=
i1,...,ik=1
Dei1+···+eik fi
(δ1,...,δk)∈Np(k)
Daδjϕij (10)
Now our goal is to prove that:
Daϕi =
Tk (11)
Our strategy of proof is as follows. We will define S1, . . . , Sp, such that
Daϕi =
Sk, Si = Ti, i = 1, . . . , p. (12)
We set
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
k=p−1
i1,...,ik+1=1
β:=ei1+···+eik+1
Dβfi ·D(c)ϕik+1
(δ1,...,δk)∈Np−1(k)
Dbδjϕij .
For m = 2, 3, . . . , p− 1 we set
k=m−1
i1,...,ik+1=1
β:=ei1+···+eik+1
Dβfi ·D(c)ϕik+1
(δ1,...,δk)∈Np−1(k)
Dbδjϕij
i1,...,ik=1
β:=ei1+···+eik
(δ1,...,δk)∈Np−1(k)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
It remains to show that Si = Ti for i = 1, . . . , p. Consider first i = 1. Recall
that N p−1(1) = {(1, 2, . . . , p− 1)}, hence
Desfi ·Db+(c)ϕs =
Desfi ·Daϕs.
Therefore
S1 = T1. (13)
Consider now i = p. For an arbitrary s > 0 N s(s) contains only one element
((1), (2), . . . , (s)). Therefore we obtain
i1,...,ip=1
Dei1+···+eip fi ·D(c)ϕip
(δ1,...,δp−1)∈Np−1(p−1)
Dbδjϕij
i1,...,ip=1
Dei1+···+eip fi ·D(c)ϕip
Dbjϕij .
Since a = b+ (c), where c = (ap), hence
i1,...,ip=1
Dei1+···+eip fi
Dajϕij
i1,...,ip=1
Dei1+···+eip fi
(δ1,...,δp)∈Np(p)
Daδjϕij = Tp
Consider now m = 2, 3, . . . , p− 1. We have
i1,...,im=1
Dei1+···+eim fi ·D(c)ϕim
(δ1,...,δm−1)∈Np−1(m−1)
Dbδjϕij
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm)∈Np−1(m)
Dbδs+(c)ϕis
j 6=s
Dbδjϕij
Using decomposition N p(m) = A ∪B as in (8) we obtain
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm−1,δm=(p))∈A
Daδjϕij
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm)∈B
Daδjϕij
i1,...,im=1
Dei1+···+eim fi
(δ1,...,δm)∈Np(m)
Daδjϕij = Tm
We have shown that Ti = Si for i = 1, . . . , p. This finishes the proof.
4 Cr-Lohner algorithm
4.1 Why one needs an Cr-algorithm?
There are several effective algorithms for the computation of rigorous bounds
for solutions of ordinary differential equations, including Lohner method [Lo],
Hermite–Obreschkoff algorithm [NJ] or Taylor models [BM]. For Cr-computa-
tions the number of equations to solve is equal to n
hence, even for
r = 1 direct application of such an algorithms to equations for variations (14)
leads to integration in high dimensional space and is usually inefficient. Let us
recall after [Z, Sec. 6] the basic reason for this. In order to have a good control
over the expansion rate of the set of initial conditions during a time step these
algorithms, while being C0, are C1 ’internally’(or higher for Taylor models),
because they solve non-rigorously equations for (∂ϕ
) - the variational matrix of
the flow. This effectively squares the dimension of phase space of the equation
and impacts heavily the computation time. But as it was observed in [Z] the
equations for partial derivatives of the flow can be seen as non-autonomous and
nonhomogenous linear equations, therefore we do not need additional equations
for variations for them. As a result the dimension of the effective phase space
for our Cr-algorithm is given by n
and not a square of this number.
Another important aspect of the proposed algorithm is the fact that the
Lohner-type control of the wrapping effect is done separately for x-variables
and variables Daϕ. This feature is not present in the blind application of C0
algorithm to the system of variational equations and it turns out that this often
practically switches off the control of the wrapping effect on x-variables, as
various choices used in this control become dominated by the Daϕ-variables.
In [Z] a C1-algorithm has been proposed. Here we present an algorithm for
computation of higher order partial derivatives.
4.2 An outline of the algorithm
Let us fix r ≤ K and consider the following system of differential equations
ϕ = f ◦ ϕ
Daϕ =
i1,...,ik=1
Dei1+···+eik f
(δ1,...,δk)∈Nd(k)
Daδjϕij
for all a ∈ Nnd , d = 1, . . . , r.
Our goal is to present an algorithm for computing a rigorous bound for the
solution of (14) with a set of initial conditions
ϕ(0, x0) ∈ [x0] ⊂ Rn
Dϕ(0, x0) = Id
Daϕ(0, x0) = 0, for a ∈ Nn2 ∪ . . . ∪ Nnr .
In the sequel we will use the following notations:
• if a solution of system (14) is defined for t > 0 and some x0 ∈ Rn, then
for a ∈ N by Va(t, x0) we denote Daϕ(t, x0)
• for [x0] ⊂ Rn by [Va(t, [x0])] we will denote a set for which we have
Va(t, [x0]) ⊂ [Va(t, [x0])]. This set is obtained using an rigorous nume-
rical routine described below.
The Cr-Lohner algorithm is a modification of C1-Lohner algorithm [Z]. One
step of Cr-Lohner is a shift along the trajectory of the system (14) with the
following input and output data
Input data:
• tk - a current time,
• hk - a time step,
• [xk] ⊂ Rn, such that ϕ(tk, [x0]) ⊂ [xk],
• [Vk,a] = [Vk,a(tk, [x0])] ⊂ Rn, such that Daϕ(tk, [x0]) ⊂ [Vk,a] for a ∈
Nn1 ∪ . . . ∪ Nnr .
Output data:
• tk+1 = tk + hk - a new current time,
• [xk+1] ⊂ Rn, such that ϕ(tk+1, [x0]) ⊂ [xk+1],
• [Vk+1,a] = [Vk+1,a(tk+1, [x0])] ⊂ Rn, such that Daϕ(tk+1, [x0]) ⊂ [Vk+1,a]
for a ∈ Nn1 ∪ . . . ∪ Nnr .
We will often skip the arguments of Vk,a when they are obvious from the context.
The values of [xk+1] and [Vk+1,a], a ∈ Nn1 are computed using one step
C1-Lohner algorithm. After it is done, we perform the following operations to
compute [Vk+1,a] for a ∈ Nn2 ∪ . . . ∪ Nnr
1. Find a rough enclosure for Daϕ([0, hk], [xk]).
2. Compute [Vk+1,a], this will also involve some rearrangement computations
to reduce the wrapping effect for V [Mo, Lo].
5 Computation of a rough enclosure for Daϕ
For a fixed multipointer a ∈ Nnd Equation (14) can be written as follows
Daϕ(t, x) = Ba(t, x) +A(t, x)Daϕ(t, x) (16)
where
i1,...,ik=1
Dei1+···+eik f
(δ1,...,δk)∈Nd(k)
Daδjϕij
A = Df ◦ ϕ
The procedure for computing the rough enclosure is based on the notion of
a logarithmic norm, which we give below.
Definition 7 [HNW] For a square matrix A the logarithmic norm µ(A) is de-
fined as a limit
µ(A) = lim sup
‖Id +Ah‖ − 1
where ‖ · ‖ is a given matrix norm.
The formulas for the logarithmic norm of a real matrix in the most frequently
used norms are (see [HNW])
1. for ‖x‖1 =
i |xi|, µ(A) = maxj(ajj +
i6=j |aij |)
2. for m ‖x‖2 =
i |xi|2, µ(A) is equal to the largest eigenvalue of (A +
AT )/2
3. for ‖x‖∞ = maxi |xi|, µ(A) = maxi(aii +
j 6=i |aij |)
In order to find bounds for Daϕ we use the following theorem [HNW, Thm.
I.10.6]
Theorem 8 Let x(t) be a solution of a differential equation
x′(t) = f(t, x(t)), x ∈ Rn (18)
Let ν(t) be a piecewise differentiable function with values in Rn. Assume that
(t, η)
≤ l(t) for η ∈ [x(t), ν(t)]
|ν′(t)− f(t, ν(t))| ≤ δ(t),
where by µ(A), we denote a logarithmic norm of a square matrix A ∈ Rn×n.
Then for t ≥ t0 we have
|x(t) − ν(t)| ≤ eL(t)
|x(t0)− ν(t0)|+
e−L(s)δ(s)ds
, (19)
with L(t) =
l(τ)dτ .
We apply the above theorem to Equation (16) to obtain
Lemma 9 Let us fix x ∈ Rn. Assume that |Ba(t, x)| ≤ δ(t) and µ(A(t, x)) ≤
l(t), then for t > t0
|Daϕ(t, x)| ≤ |Daϕ(t0, x)|eL(t) + eL(t)
e−L(τ)δ(τ)dτ (20)
with L(t) =
l(τ)dτ .
Proof: Consider Equation (16) and a homogenous problem for (16)
w = f(t, w) := A(t, x) · w, w ∈ Rn. (21)
Using Theorem 8 we can estimate the difference between any solution of (21),
w, and a solution of (16), denoted by Daϕ.
|Daϕ(t) − w(t)| ≤ |Daϕ(t0)− w(t0)|eL(t) + eL(t)
e−L(τ)δ(τ)dτ. (22)
After a substitution w(t) = 0, which is a solution of the homogenous equation,
we obtain our assertion.
Usually, we do not have any control over the time dependence of δ and l,
hence we will use the following
Lemma 10 Assume that |Ba(t, x)| ≤ δ and µ(A(t, x)) ≤ l for t ∈ [0, h] then
for t ∈ [0, h] we have
|Daϕ(t, x)| ≤ |Daϕ(0, x)|max(1, ehl) + δ
elt − 1
, if l 6= 0, (23)
|Daϕ(t, x)| ≤ |Daϕ(0, x)|+ δt, when l = 0. (24)
5.1 The procedure for the computation of the rough en-
closure for V .
The procedure for the computing of the rough enclosure is iterative, which
means that given a rough enclosure for ϕ([0, hk], [xk]) and rough enclosures
Daϕ([0, hk], [xk]) for all a ∈ Nn1 ∪ . . . ∪ Nnp we are able to compute the rough
enclosure for Daϕ([0, hk], [xk]) for a ∈ Nnp+1.
The procedures for computation of the rough enclosures of ϕ([0, hk], [xk])
and Daϕ([0, hk], [xk]) for a ∈ Nn1 has been given in [Z]. Below we present an
algorithm for computing [Ea] for a ∈ Nn2 ∪ . . . ∪Nnr .
Input parameters:
• hk - a time step,
• [xk] ⊂ Rn - the current value of x = ϕ(tk, [x0]),
• [E0] ⊂ Rn - a compact and convex such that ϕ([0, hk], [xk]) ⊂ [E0]
• [Ea] ⊂ Rn, a ∈ Nn1 ∪ . . . ∪ Nnp such that Daϕ([0, hk], [xk]) ⊂ [Ea] for
a ∈ Nn1 ∪ . . . ∪ Nnp .
Output:
• [Ea] ⊂ Rn, a ∈ Nnp+1 such that
Daϕ([0, hk], [xk]) ⊂ [Ea]
Before we present an algorithm let us observe that for a fixed a ∈ Nnp+1, Ba
defined in (17) could be seen as a multivariate function of t, x and Vb = Dbϕ for
b ∈ Nn1 ∪. . .∪Nnp . More precisely, putmp := ♯
Nn1 ∪ . . . ∪ Nnp
, where ♯ stands
for number of elements of a set. Recall that, we have defined by (6) a linear
order in Nn. Hence, there is a unique sequence of multipointers b1, . . . , bmp ,
such that bi ∈ Nn1 ∪ . . .∪Nnp for i = 1, . . . ,mp, b1 ≤ b2 ≤ · · · ≤ bmp and bi 6= bj
for i 6= j.
Let us define
B̃a : R× (Rn)mp+1 → Rn,
Fa : R× (Rn)mp+1 → Rn
B̃a(t, x, vb1 , . . . , vbmp ) =
i1,...,ik=1
Dei1+···+eik f(ϕ(t, x))
(δ1,...,δk)∈Np+1(k)
Fa(t, x, vb1 , . . . , vbm) = B̃a(t, x, vb1 , . . . , vbm) +Df(ϕ(t, x))Va(t, x) (26)
Algorithm:
To compute [Ea] for a ∈ Nnp+1 we proceed as follows
1. Find l ≥
maxx∈[E0] µ (Df(x))
2. Compute δa ≥ max ‖B̃a‖, i.e.
δa ≥ max
(x,vb1 ,...,vbmp )∈[E0]×[Eb1 ]×···×[Ebmp ]
∥B̃a(0, x, vb1 , . . . , vbmp )
For example, if a = (j, c) ∈ Nn2 , then δa should be such that
δa ≥ max
x∈[E0],v1∈[E(1)],...,vn∈[E(n)]
r,s=1
∂xr∂xs
(x) (vj)s (vc)r
3. Define [Ea]i = [−1, 1]δa e
, for i = 1, . . . , n, where [Ea]i denotes i-th
coordinate of [Ea].
One can refine the obtained enclosure by
[Ea] :=
[0, hk]Fa
0, [E0], [Eb1 ], . . . , [Ebmp ]
∩ [Ea]
Indeed, for i = 1, . . . , n, t ∈ [0, hk] and x0 ∈ [E0] we have
Daϕi(t, x0) = Daϕi(t, x0)−Daϕ(0, x0)
= t (Fa)i (θi, x0, Db1ϕ(θi, x0), . . . , Dbmpϕ(θi, x0))
= t (Fa)i (0, ϕ(θi, x0), Db1ϕ(θi, x0), . . . , Dbmpϕ(θi, x0))
for some θi ∈ [0, t] ⊂ [0, hk]. In the above we have used the fact that
Fa(t, x, v1, . . . , vmp) = Fa(0, ϕ(t, x), v1, . . . , vmp).
Since ϕ(θi, x0) ∈ [E0] and Dbjϕ(θi, x0) ∈ [Ebj ] for j = 1, . . . ,mp we get
Daϕi(t, x0) ∈ [0, hk] (Fa)i
0, [E0], [Eb1 ], . . . , [Ebmp ]
6 Computation of [Vk+1]
6.1 Composition formulas
For any p-times continuously differentiable functions f, g : Rn → Rn and a ∈ Nnp
we have
Da(f ◦ g) =
i1,...,ik=1
Dei1+···+eik fi
(δ1,...,δk)∈Np(k)
Daδj gij (27)
We can apply the above formula to f = ϕ(hk, ·) and g = ϕ(tk, ·) to obtain
Va(tk + hk, x0) =
i1,...,ik=1
VΛ−1(ei1+...+eik )
(hk, xk)
(δ1,...,δk)∈Np(k)
(tk, x0)
for all x0 ∈ [x0]. Using notations [Vk+1,a] := [Va(tk + hk, [x0])] and [Vk,a] =
[Va(tk, [x0])] we can rewrite the above equation as
[Vk+1,a] =
i1,...,ik=1
VΛ−1(ei1+...+eik )
(hk, [xk])
(δ1,...,δk)∈Np(k)
Vk,aδj
where Λ is defined by (5).
6.2 The procedure for computation of [Vk+1]
We introduce new parameters od - the order of the Taylor method used in
computations of Va for a ∈ Nnd . It makes sense to take o1 ≥ o2 ≥ · · · ≥ or.
Input parameters:
• hk - a time step,
• [xk] ⊂ Rn - the current value of x = ϕ(tk, [x0]),
• [Vk,a] ⊂ Rn - a current value of Vk,a(tk, [x0]), for a ∈ Nn1 ∪ . . . ∪Nnr
• [E0] ⊂ Rn compact and convex, such that ϕ([0, hk], [xk]) ⊂ [E0] - a rough
enclosure for [xk],
• [Ea] ⊂ Rn, compact and convex, such that Daϕ([0, hk], [xk]) ⊂ [Ea], for
a ∈ Nn1 ∪ . . . ∪ Nnr .
Output: [Vk+1,a] ⊂ Rn, such that
Va(tk + hk, x0) ∈ [Vk+1,a] (29)
for x0 ∈ [x0] and a ∈ Nn1 ∪ . . . ∪ Nnr .
Algorithm: We compute [Vk+1] as follows
1. Computation of Va(hk, [xk]) using Taylor method for Equation (14), i.e. for
a ∈ Nnp we compute
[Fa] =
dti−1
Fa(0, [xk], Vb1 , . . . , Vbmp−1 ) (30)
hop+1
(op + 1)!
Fa(0, [E0], [Eb1 ], . . . , [Ebmp−1 ]).
where Vbi = 0 for bi ∈ Nn2 ∪ . . . ∪ Nnp−1 and V(j) = enj for j = 1, . . . , n.
Observe that
Va(hk, [xk]) ⊂ [Fa] (31)
Indeed, using Taylor series expansion we obtain that for xk ∈ [xk] and
j = 1, . . . , n holds
(Va)j(hk, xk) =
dti−1
(Fa)j(0, xk, Vb1(0, xk), . . . , Vbmp−1 (0, xk))
hop+1
(op + 1)!
(Fa)j(θi, xk, Vb1(θi, xk), . . . , Vbmp−1 (θi, xk))
for some θi ∈ [0, hk]. Observe, that
(Fa)j(θi, xk, Vb1(θi, xk), . . . , Vbmp−1 (θi, xk))
(Fa)j(0, ϕ(θi, xk), Vb1(θi, xk), . . . , 0, Vbmp−1 (θi, xk))
Using ϕ(θi, xk) ∈ [E0] and Vbs(θi, xk) ∈ [Ebs ] for s = 1, . . . ,mp−1 we
obtain our assertion.
2. The composition. Put
[Jk] := ([F(1)], . . . , [F(n)])
Using (28) for a ∈ Nnp we have
[Vk+1,a] = [αa] + [Jk] · [Vk,a], (32)
where
[αa] =
i1,...,ik=1
[FΛ−1(ei1+...+eik )
(δ1,...,δk)∈Np(k)
Vk,aδj
In our implementation of the algorithm we use the symbolic differentiation to
obtain formulae for Daf . Next, using the automatic differentiation we compute
Fa(t, x, Vb1 (t, x), . . . , Vbmp−1 (t, x))|t=0 which appear in (30).
6.3 Rearrangement for Va - the evaluation of Equation
It is well know that a direct evaluation of Equation (32) leads to wrapping effect
[Mo, Lo]. To avoid it following the work of Lohner [Lo] we will use the same
scheme as it was proposed in [Z].
Namely, observe that Equation (32) has exactly the same structure as the
propagation equations for C1-method (see [Z, Section 3]). Moreover, all vectors
Vk,a, for a ∈ Nn1 ∪ . . .Nnr ’propagate’ by the same [Jk] as did the variational
part in [Z], hence it makes sense the same approach.
To be more precise, each set [Vk,a], for a ∈ Nn1 ∪ . . . ∪ Nnr is represented in
the following form
[Vk,a] = vk,a + [Bk][rk,a] + Ck[qk,a]
where [Bk] is interval matrix, Ck is point matrix, vk,a is a point vector and
rk,a, qk,a are interval vectors. Observe that [Bk] and Ck are independent of a.
In the sequel we will drop index a. Equation (32) leads to
[Vk+1] = [α] + [Jk](vk + [Bk][rk] + Ck[qk]) (34)
Let m([z]) denotes a center of an interval object, i.e. [z] is interval vector or
interval matrix and ∆([z]) = [z]−m([z]).
Let [Q] be an interval matrix which contains an orthogonal matrix. Usually,
[Q] is computed by the orthonormalisation of the columns of m([Jk])[Bk].
[Z] = m([Jk])Ck
Ck+1 = m([Z])
[Bk+1] = [Q]
Then we rearrange formula (34) as follows
[s] = [α] + [Jk]vk +∆([Jk])[Vk]
vk+1 = m([s])
[qk+1] = [qk]
[rk+1] = [Q
T ](∆([s]) + ∆([Z])[qk]) + ([Q
T ]m([Jk])[Bk])[rk]
Summarizing, we can use the following data structure to represent ϕ(tk, [x0])
and Daϕ(tk, [x0]), for a ∈ Nn1 ∪ . . . ∪ Nnr
type CnSet = record
v0, r0, q0: IntervalVector;
C0, B0, C,B : IntervalMatrix;
{va, ra, qa : IntervalVector}a∈Nn1 ∪...∪Nnr
The set ϕ(tk, [x0]) is represented as v0 + B0r0 + C0q0, the partial derivatives
Daϕ(tk, [x0]) are represented as va+Bra+Cqa. The matrices B,C are common
for all partial derivatives.
Notice, that if we start the Cr computation with an initial condition (15) then
there is no Lipschitz part at the beginning for the partial derivatives. Hence,
the initial values for C and B are set to the identity matrix and the initial values
for qa, ra are set to zero.
If the interval vectors ra become ’thick’ (i.e. theirs diameters are larger than
some threshold value) we can set a new Lipschitz part in our representation (it
must be done simultaneously for all Daϕ) and reset ra in the following way
qa = ra + (B
TC)qa, for a ∈ Nn1 ∪ . . . ∪ Nnr
ra = 0, for a ∈ Nn1 ∪ . . . ∪ Nnr
C = B
B = Id
A similar change of the Lipshitz part may be done when vectors ra become thick
in comparison to qa.
7 Derivatives of Poincaré map
Consider a differential equation
x′ = f(x), x ∈ Rn, f ∈ CK+1 (36)
Let ϕ : R × Rn → Rn be a (local) dynamical system induced by (36). Let
α : Rn → R be C1-map. Put Π = {x | α(x) = C}.
Definition 11 We will say that Π is a local section for the vector field f at
y0 ∈ Π if
〈∇α(y0)|f(y0)〉 6= 0. (37)
Assume x0 ∈ Rn and t0 ∈ R are such that Π is a local section at ϕ(t0, x0).
Consider an implicit equation
α(ϕ(tP (x), x)) = C. (38)
It follows easily from (37) and from the implicit function theorem that there
exists a uniquely defined tP : R
n−→◦ R in a neighborhood of x0, such that
tP (x0) = t0. The function tP is as smooth as the flow ϕ. We will refer to
tP as to the Poincare return time to section Π.
We define a Poincaré map P : Rn ⊃ dom (tP ) → Rn by
P (x) = ϕ(tP (x), x). (39)
Usually the Poincaré map is defined as a map P : Π1−→◦ Π2, where Π1,Π2 are
local sections in Rn. The approach taken here, i.e. treating the Poincaré map
as map P : Rn−→◦ Rn allows us to not to worry about the coordinates on local
section.
In this section we are interested in the partial derivatives of P defined by
(39).
From (39) we can compute ∂Pi
and we obtain
(x) = fi(P (x))
(x) +
(tP (x), x). (40)
We need ∂tP
. We differentiate (38) to obtain
(P (x))
fk(P (x))
(x) +
(tP (x), x)
(∇α(P (x)) · f(P (x))) ∂tP
(x) +
(P (x))
(tP (x), x) = 0. (41)
Hence
(x) = − 1〈∇α(P (x))|f(P (x))〉
(P (x))
(tP (x), x). (42)
7.1 Higher order derivatives of the Poincaré map
To make formulas transparent we will drop arguments of functions in this sec-
tion, but reader should be aware that for tP and its partial derivatives the
argument is x, for ϕ and Daϕ the argument is always the pair (tP (x), x).
From (40) we obtain
D(j,c)P =
ϕD(j)tPD(c)tP +
D(c)ϕD(j)tP +
ϕD(j,c)tP
D(j)ϕD(c)tP +D(j,c)ϕ.
It is easy to see that partial derivatives of high order give rise to quite complex
expressions and it is not entirely obvious how to organize it in some coherent
and programmable way. For this purpose we use the following
Lemma 12 For a multipointer a ∈ Nnp we have
DaP = Daϕ+
(δ1,...,δk)∈Np(k)
j=1Daδj tP
(δ1,...,δk)∈Np(k)
∂tk−1
Daδsϕ
j 6=sDaδj tP
Proof: By induction on p. For p = 1 formula (43) is equivalent to (40), because
the two last sums are taken over empty set. Assume (43) holds true for some
p ≥ 1 and fix a ∈ Nnp+1. Our goal is to show that
DaP = R1 +R2 +R3
where
R1 = Daϕ+
ϕDatP
(δ1,...,δk)∈Np+1(k)
Daδj tP
(δ1,...,δk)∈Np+1(k)
∂tk−1
Daδsϕ
j 6=s
Daδj tP
Write a = β + γ, where β ∈ Nnp and γ = (ap+1) ∈ Nn1 . From the induction
assumption we have
DaP = Dγ
ϕDβtP
(δ1,...,δk)∈Np(k)
j=1Dβδj tP
(δ1,...,δk)∈Np(k)
∂tk−1
Dβδsϕ
j 6=sDβδj tP
i=1 Si
where
S1 = Daϕ+
ϕDatP
DβϕDγtP
ϕDβtPDγtP
DγϕDβtP
(δ1,...,δk)∈Np(k)
j=1Dβδj tP
∂tk+1
ϕDγtP
(δ1,...,δk)∈Np(k)
j=1Dβδj tP
(δ1,...,δk)∈Np(k)
s=1Dβδs+γtP
j 6=s
Dβδj tP
(δ1,...,δk)∈Np(k)
∂tk−1
Dβδs+γϕ
j 6=sDβδj tP
(δ1,...,δk)∈Np(k)
DβδsϕDγtP
j 6=sDβδj tP
S10 =
(δ1,...,δk)∈Np(k)
r 6=s
∂tk−1
DβδsϕDβδr+γtP
j 6=s
j 6=r
Dβδj tP
Obviously R1 = S1. We will show that R2 = S3 + S6 + S7 and R3 = S2 + S4 +
S5 + S8 + S9 + S10.
Denote by Ri,k, i = 2, 3 a part of sum Ri with fixed k = 2, . . . , p + 1.
Similarly, let us denote by Si,k a part of sum Si, i = 5, . . . , 10, for k = 2, . . . , p.
Using decomposition of N p+1(2) as in (8) we obtain that R2,2 = S3 + S7,2.
Similarly, using (8) we observe that R2,k = S6,k−1 + S7,k for k = 3, . . . , p.
Finally, since N p+1(p + 1) = {((1), (2), . . . , (p + 1))} and γ = (ap+1) we find
that R2,p+1 = S6,p. This shows that R2 = S3 + S6 + S7.
It remains to show that R3 = S2 +S4 +S5 + S8 + S9 + S10. We will classify
possible terms by the fact, where p+ 1 appears in δi, i = 1, . . . , k and how this
δi enters in R3 as δs or δj . There are four cases
1. δs = (p+ 1)
2. δj = (p+ 1)
3. p+ 1 ∈ δs, |δs| ≥ 2
4. p+ 1 ∈ δj, |δj | ≥ 2
Let us fix k = 2. Let (δ1, δ2) ∈ N p+1(2). The term for case 1 is S4, for case 2
is S2, case 3 is S8,2 and case 4 is S10,2. Hence, R3,2 = S2 + S4 + S8,2 + S10,2.
For k = 3, . . . , p and fixed (δ1, . . . , δk) ∈ N p+1(k) we have: case 1 is given
by S5,k−1, case 2 by S9,k−1, case 3 by S8,k and case 4 by S10,k Hence, for
k = 3, . . . , p we have R3,k = S5,k−1 + S9,k−1 + S8,k + S10,k.
Finally, for k = p+ 1 we observe, that R3,p+1 = S5,p + S9,p. Indeed, in this
case (δ1, . . . , δp+1) = ((1), (2), . . . , (p + 1)). Hence, either for δs = γ we have
term S5,p and δs 6= γ we have S9,p.
We have showed that R3 = S2 + S4 + S5 + S8 + S9 + S10 and the proof is
finished.
Hence, if we know all the partial derivatives of tP up order p we can compute
the partial derivatives of the Poincaré map up the same order. In next subsection
we show how to compute partial derivatives of tP for affine sections.
7.2 Partial derivatives of tP for affine sections
Assume α : Rn → R is an affine map given by
α(x) = α0 +
αixi.
This is a quite restrictive assumption about sections, but it leads to relatively
simple formulas for DatP and it is sufficient for the applications we have in
mind.
Lemma 13 For a multipointer a ∈ Nnp holds
−DatP
∇α| ∂
= 〈∇α|Daϕ〉
∇α| ∂
(δ1,...,δk)∈Np(k)
j=1Daδj tP
(δ1,...,δk)∈Np(k)
∇α| ∂
∂tk−1
Daδsϕ
j 6=sDaδj tP
Proof: The proof is a direct consequence of Lemma 12 and (38). Since α is
affine, by differentiating of α(P (x)) = C we get 〈∇α|DaP 〉 = 0. Using formula
(43) for DaP we obtain our assertion.
Fix [x] ⊂ Rn and assume we have a rigorous bound for tP ([x]) ∈ [t1, t2] (see
[Z, Section 6] for more details on this). Lemmas 13 and 12 show that given rigo-
rous bounds for the partial derivatives Daϕ([t1, t2], [x]) and
Daϕ([t1, t2], [x])
up to some order p we can compute recursively rigorous bounds for the partial
derivatives of tP ([x]) and P ([x]) up to the same order. Notice, that
Daϕ are
given by Taylor coefficients of the solution of (14) with initial conditions P ([x])
for C0 part and Daϕ(tP (x), [x]) for equations for variations. Hence, these coef-
ficients can be easily computed using the automatic differentiation algorithm.
8 Applications.
One of the typical invariant sets in hamiltonian mechanics are invariant tori.
However, the existence of invariant torus in a given system is often difficult to
prove despite the fact that the theory is quite well developed. Probably the
best work in this direction was done by Celletti and Chercia [CC1, CC2], where
the an effective application (computer assisted proof) of KAM theory to the
restricted three body problem modelling system consisting of Sun, Jupiter and
asteroid 12 Victoria was given. Our aim here is more modest as we focus on the
invariant tori emanating from the elliptic fixed point satisfying suitable twist
condition.
In this section we show that the rigorous computations of partial derivatives
of a dynamical system up to order 3 or 5 can be used to prove that in a particular
system an invariant torus exists around some elliptic periodic orbits. In this
section this will be done for the forced pendulum equation and the Michelson
system.
8.1 Area preserving maps on the plane, normal forms and
KAM theorem
Definition 14 Let f : R2 → R2 be a smooth area preserving map, such that
f(p) = p. Let λ and µ be eigenvalues of df(p). Following [SM] we will call the
point p
• hyperbolic if λ, µ ∈ R and λ 6= µ,
• elliptic if λ = µ and λ 6= µ,
• parabolic if λ = µ.
The following KAM theorem will be the main tool to prove the existence of
invariant tori in this paper.
Theorem 15 [SM, §32] Consider an analytic area preserving map f : R2 →
2, f(r, s) = (r1, s1) where
r1 = r cosα− s sinα+O2l+2
s1 = r sinα+ s cosα+O2l+2 (44)
r2 + s2
and O2l+2 denotes convergent power series in r, s with terms of order greater
than 2l+ 1, only.
If at least one of γ1, . . . , γl is not zero then the origin is a stable fixed point
for map f . Moreover, in any neighborhood U of point 0 there exists an invariant
curve for map f around the origin contained in U .
The next theorem and its proof tells how to bring a planar area preserving
map in the neighborhood of an elliptic fixed point into the form (44).
Theorem 16 [SM, §23] Consider an analytic area preserving map f : R2 → R2
such that f(0) = 0. Let λ, λ̄ be complex eigenvalues of Df(0), such that |λ| =
|λ̄| = 1. If λk 6= 1 for k = 1, . . . , 2l+2, then there is an analytic area preserving
substitution such that in the new coordinates mapping f has form (44).
The proof of the above theorem is constructive, i.e. given the power series
for f at an elliptic fixed point one can construct explicitly an area preserving
substitution and compute the coefficients γ0, . . . , γl in (44). An explicit formula
for the coefficient γ1 in the above normal form is given in Appendix A.
8.2 The existence of invariant tori in forced pendulum.
Consider an equation
θ̈ = − sin(θ) + sin(ωt) (45)
Observe that (45) is hamiltonian.
Let us denote by Pω : R
2−→◦ R2 the Poincaré map for Equation (45) with
a parameter ω, i.e. Pω = ϕ(2π/ω, ·), where ϕ : R × R2−→◦ R2 is a local flow
induced by (45). Observe that (45) is nonautonomous, but it is equivalent to
first order system of autonomous ODE given by
= − sin(θ) + sin(ωt) (46)
In the sequel all rigorous computations for (45) will be in fact performed for the
system (46).
Observe that to any invariant closed curve for Pω corresponds and invariant
2-torus for (45).
Consider a set of parameter values
Ω1 = [2, 2.994], Ω2 = [3, 3.997], Ω3 = [4, 8]
Ω = Ω1 ∪ Ω2 ∪ Ω3
The following lemma was proved with computer assistance
Lemma 17 For all parameter values ω ∈ Ω there exists an elliptic fixed point
xω ∈ R2 for Pω. Moreover, there exists an area-preserving substitution such
that in the new coordinates the map fω(x) = Pω(x+xω)−xω has the form (44)
with l = 1 and γ1 6= 0.
Before we give the proof, let us briefly comment about the choice of the parame-
ter set Ω. For parameter values slightly lower than 2 we observe the parabolic
case, i.e. there exists a parameter value ω1 for which eigenvalues of the deriva-
tive of Pω1 are equal to −1. In two gaps in Ω below 3 and 4 we have resonances
of low order. Namely, we have parameter values with an elliptic fixed with ei-
genvalues to e±2π/3 = −1
i and e±iπ/2 = ±i, respectively. Clearly, in a
computer assisted proof we need to exclude a small interval around those para-
meters. For ω > 4 it seems that the interval Ω3 can be extended much further
to the right without any difficulty.
Proof of Lemma 17: A computer assisted proof consists of the following steps.
We cover the set Ω by 9910 nonequal subintervals ωi. Diameters of ωi’s were
relatively large for values far away from the parabolic cases and very small close
to them. For a fixed subinterval ωi we proceed as follows
1. Let ω̄ denote an approximate center of the interval ωi. We find an approxi-
mate fixed point for Pω̄ using the standard nonrigorous Newton method.
Let us denote such a point by xi.
2. We define a box centered at xi, i.e we set vi := xi + [−εi, εi]2, where
εi > 0 depends on subinterval ωi - the values we used are from the interval
[5 · 10−6, 3 · 10−3], depending on whether xi close to parameter values
corresponding to parabolic cases.
3. Using the C1-Lohner algorithm we compute the Interval Newton operator
[Mo, N, A] Ni := N(Pωi − Id, xi, vi) and verify that Ni ⊂ int vi. This
proves that for all ω ∈ ωi there exists a unique fixed point xω ∈ Ni for
4. Using the C3-Lohner algorithm we compute a rigorous bound for Pωi(Ni)
and DαPωi(Ni), α ∈ N21 ∪N22 ∪N23. Hence, we obtain a rigorous bound for
the coefficients in
fω(x) =
|α|=1
DaP (xω)x
5. We show that an arbitrary matrix M ∈ DPωi(Ni) has a pair of complex
eigenvalues λ, λ̄ which satisfy λk 6= 1 for k = 1, . . . , 4. From Theorem 16
it follows there exists an area-preserving substitution such that in the new
coordinates the map fω for ω ∈ ωi has the form (44) with l = 1.
6. We compute a rigorous bound for γ0 and γ1 which appear in the formula
(44) and verify that for ω ∈ ωi holds γ1 6= 0.
The rigorous bounds for the values of γ1 on Ω are
γ1(Ω1) ⊂ [0.29930416771330087, 30.118260918229566]
γ1(Ω2) ⊂ [0.099747909112924596, 0.56550301088840627]
γ1(Ω3) ⊂ [0.18574835001593507, 0.4129279974577012]
A computer assisted proof of the above took approximately 95 minutes on the
Pentium IV 3GHz processor.
As a straightforward consequence of Lemma 17 and Theorem 15 we obtain
Theorem 18 For all parameter values ω ∈ Ω there exists an elliptic fixed point
xω ∈ R2 for Pω. Moreover, any neighborhood of point xω contains an invariant
curve for Pω around xω.
8.3 Higher order normal forms.
In the previous section it was shown that C3 computations are sufficient to
prove that for (45) a family of invariant tori exists. However, it may happen
that the coefficient γ1 in the normal form vanishes. In this situation we may try
to compute higher order normal form. As an example we consider a pendulum
with a different forcing term,
θ̈ = − sin(θ) + sin(ωt) + sin(2ωt). (47)
Theorem 19 Let Pω be the Poincaré map for (47). For all parameter values
ω ∈ Ω∗ = [2.9957694795, 2.9957694796] there exists an elliptic fixed point xω ∈
2 for Pω. Moreover, any neighbourhood of point xω contains an invariant
curve for Pω around xω.
Proof: The main concept of the proof is the same as in Lemma 17. Using the
nonrigorous Newton method we find an approximate fixed point
x = (−7.7491573604896152 · 10−12,−0.54723831527031352).
We set v = x + 3 · 10−5([−1, 1] × [−1, 1]). Using the C1-Lohner algorithm we
compute the Interval Newton Operator of Pω − Id on v and we obtain that for
all ω ∈ Ω∗, N = N(Pω − Id, center(v), v) ⊂ (N1, N2), where
N1 = [−5.1582932672798325, 5.1582631625020222] · 10−10
N2 = [−0.54723831580217108,−0.54723831470891193]
Since N ⊂ v we conclude that for all ω ∈ Ω∗ there exists a unique fixed point
xω ∈ N for the Poincaré map.
Using C5-Lohner algorithm we compute a rigorous bound for PΩ∗(N) and
DαPΩ∗(N), α ∈ N21 ∪ . . . ∪ N25. Hence, we obtain a rigorous bound for the
coefficients in
fω(x) =
|α|=1
DαP (xω)x
α +O6
We show that an arbitrary matrix M ∈ DPΩ∗(N) has a pair of complex ei-
genvalues λ, λ̄ which satisfy λk 6= 1 for k = 1, . . . , 6. From Theorem 16 it follows
there exists an area-preserving substitution such that in the new coordinates the
map fω for ω ∈ Ω∗ has the form (44) with l = 2.
Next, we compute a rigorous bound for γ1 and γ2 which appear in the formula
(44) and we get
γ1(Ω∗) ⊂ [−5.3924276719042241, 5.381714805052106] · 10−6
γ2(Ω∗) ⊂ [199.95180660157078, 199.99104965939162]
Since for ω ∈ Ω∗, γ2(ω) 6= 0 the assertion follows from Theorem 15.
The main observation which makes this example interesting is that there
exists ω∗ ∈ Ω∗ for which γ1(ω∗) = 0 and we cannot conclude the existence of
invariant tori for all ω ∈ Ω∗ from C3 computations. To be more precise, we
computed the coefficient γ1 for the parameter values ω1 = minΩ∗ and ω2 =
maxΩ∗ and we get
γ1(ω1) ∈ [−2.3559594437885885,−1.3593457220363871] · 10−8
γ1(ω2) ∈ [2.9671154858524365 · 10−9, 1.2819312939263052 · 10−8]
Since γ1 exists for all ω ∈ Ω∗ and depends continuously on ω we conclude, that
γ1(ω∗) = 0 for some ω∗ ∈ Ω∗.
8.4 Application to the Michelson system
The existence of an invariant curve for a planar map f : R2 → R2 can be proven
without assumption that f is measure preserving. The key assumption in the
proof given in [SM] is that any curve γ around an elliptic point intersect its
image under f , i.e. f(γ)∩ γ 6= ∅. Such a situation is also observed in reversible
planar map around an symmetric elliptic fixed points.
Definition 20 An invertible transformation M : Ω −→ Ω is called a reversing
symmetry of a local dynamical system φ : T × Ω −→ Ω, T = R or T = Z if the
following conditions are satisfied
1. if (t, x) ∈ dom (φ) then (−t, S(x)) ∈ dom (φ).
2. S(φ(t, x)) = φ(−t, S(x)))
Remark 21 In the discrete time case, the above two conditions are equivalent
to identity
M ◦ f = f−1 ◦M.
where f = φ(1, ·) is a generator of φ.
Definition 22 Let φ : T×Ω → Ω be a local (discrete or continuous) dynamical
system. For x ∈ Ω put
I(x) = {t ∈ T : (t, x) ∈ dom (φ)}
O (x) = {φ(t, x) ∈ Ω : t ∈ I(x)}
The set O (x) will be called a trajectory of a point x.
Definition 23 Assume S is an reversing symmetry for φ : T × Ω → Ω. An
orbit O (x) is called S-symmetric orbit if O (x) = S(x).
Remark 24 [La] In continuous case the orbit O (x) is S-symmetric if it contains
a point from the set Fix(S) = {y : S(y) = y}.
Remark 25 [Wi2, Lem.3.3] It is easy to see that if Θ ⊂ Ω is a Poincaré section
for a R-reversible flow φ : R × Ω → Ω such that Θ = R(Θ) then the Poincaré
map P : Θ → Θ is R|Θ-reversible.
As we observed at the beginning of this section, an R-reversible planar map
may admit an invariant curve around an R-symmetric elliptic fixed point. In
reversible case a planar map admits the same normal form around symmetric, el-
liptic fixed point as in the area-preserving case and the substitution which tends
the map to the normal form is exactly the same as we described in Appendix A
– for details see [Se, BHS].
Consider an ODE 
ẋ = y
ẏ = z
ż = c2 − y − 1
On one hand, the system (48) is an equation for the steady state solution of
one-dimensional Kuramoto-Sivashinsky PDE and it is known in the literature
as the Michelson system[Mi]. On the other hand, this system appears as a part
of the limit family of the unfolding of the nilpotent singularity of codimension
three (see [DIK1]).
The system (48) is reversible with respect to the symmetry
R : (x, y, z, t) → (−x, y,−z,−t) (49)
and since the divergence vanishes it is also volume preserving.
A dynamical system induced by (48) exhibits several types of dynamics for
different values of parameter. For sufficiently large c there is a simple invariant
set consisting of two equilibria (±c
2, 0, 0) and heteroclinic orbit between them
[MC]. Lau [Lau] numerically observed that when the parameter c decreases a
cascade of cocoon bifurcations occurs and at the limit value c ≈ 1.266232337
a periodic orbit is born through a saddle-node bifurcation. This hypothesis
has been proved in [KWZ]. The computer assisted proof of this fact given in
[KWZ] uses the algorithm presented in this paper in order to compute partial
derivatives up to second order for a certain Poincaré map.
For the parameter value equal to one and slightly smaller than one it was
proven in [DIK2, Wi1, Wi2, Wi3] that the system has rich and complicated
dynamics including symbolic dynamics, heteroclinic solutions, Shilnikov homo-
clinic solutions.
However, as the bifurcations diagram presented by Michelson suggests [Mi,
Fig.1] for all parameter values c ∈ (0, 0.3195) there are at least two elliptic
periodic orbits with large invariant islands around them. In this section we
present a proof that such islands exist for some range of parameter values. The
main idea of the proof is almost the same as in the previous section. There are
two main differences. First, the Poincaré map will not be a time shift. Therefore
computations of the partial derivatives of the Poincaré map require Lemma 12
and Lemma 13. Second difference is: we use the shooting method instead of
the interval Newton method for the proof of the existence of symmetric periodic
orbit.
The aim of this section is to prove the following
Theorem 26 For all parameter values from the set
C = C1 ∪ C2 = [0.1, 0.225]∪ [0.226, 0.25]
there exists a symmetric elliptic periodic orbit for the Michelson system (48).
Moreover, each neighbourhood of such an orbit contains a 2D tori invariant
under the flow generated by the Michelson system.
Let us define the Poincaré section Π := {(0, z, y) : z, y ∈ R}. Let Pc =
(P1, P2) : Π−→◦ Π be the Poincaré map for the system with the parameter value
c. Notice, that Pc is in fact a half Poincaré map, which means that the trajectory
of x crosses Π in opposite directions when passing through x and Pc(x), and
therefore periodic orbits for the Michelson system corresponds to periodic points
for P 2.
Since the section Π is invariant under symmetry (x, y, z) → (−x, y,−z),
from Remark 25 the Poincaré map is also reversible with respect to an involu-
tion R(y, z) = (y,−z). We will use the same letter R to denote the reversing
symmetry of the Poincaré map and the Michelson system.
Let us comment about the choice of the set C. In the gap between intervals
C1 and C2 there is a parameter value c∗ for which the eigenvalues of the Poincaré
map P 2c∗ are ±i. Apparently at this parameter value we have a bifurcation and
four periodic islands are born as it is shown in Fig.1 - see also a movie mpp.mov
available at [Wi4] which presents an animation of the phase portrait of Pc for
the parameter values from the range [0.1, 0.25].
Proof of Theorem 26. The main concept of the proof is quite similar to the
one presented in Lemma 17. We divide the set C of parameter values onto 20800
nonequal parts (smaller when close to the bifurcation parameter c∗ and close to
0.1 and 0.25). For a fixed subinterval ci from the grid we proceed as follows
1. Let c̄ denote a center of the interval ci. We find an approximate fixed point
of P 2c̄ using the standard nonrigorous Newton method. Let us denote this
point by (yi, zi).
2. Since the map Pc is reversible one can prove the existence of the fixed
point for P 2c using the shooting method as follows.
Let Fix(R) = {(y, z) ∈ Π : R(y, z) = (y, z)} = {(y, 0) ∈ Π : y ∈ R}.
Since Pc satisfies (Pc ◦ R)2 = Id whenever the left side is defined, one
can see that if x ∈ Fix(R) and Pc(x) ∈ Fix(R) then P 2c (x) = x. Let us
remark, that we always get an approximate fixed points (yi, zi) resulting
from the nonrigorous Newton method very close to Fix(R). We define
two points u1 = (yi − εi, 0), u2 = (yi + εi, 0) ∈ Fix(R), where εi is a small
number depending on ci and we show that πz(Pci(u1)) · πz(Pci(u2)) < 0,
where πz is a projection onto z coordinate. Hence, if the Pci is defined
on the set Ni = (0, [yi − εi, yi + εi], 0) then for all parameter values c ∈ ci
there is a point uc ∈ N which satisfies πz(Pc(uc)) = 0 and therefore
Figure 1: Phase portrait of the Poincaré map Pc (top) before bifurcation for
c = 0.225 and (bottom) after bifurcation for c = 0.226 with four periodic islands.
Between those parameters resonant case occurs with eigenvalues equal to ±i.
See also auxiliary material [Wi4].
Pc(uc) ∈ Fix(R). This shows that for all c ∈ ci there exists a fixed point
for P 2c inside Ni provided Pc is defined on Ni, which will be discussed
below.
3. Using C3-Lohner algorithm we compute rigorous bounds for P 2ci(Ni) and
DαP 2ci(Ni) for α ∈ N
1 ∪ N22 ∪ N23. This implies also that Ni ⊂ domPci .
4. We show that an arbitrary matrix M ∈ DPci(Ni) has a pair of complex
eigenvalues λ, λ̄ which satisfy λk 6= 1 for k = 1, . . . , 4. From Theorem 16
it follows there exists an area-preserving substitution such that in the new
coordinates the map Pc for c ∈ ci has the form (44) with l = 1.
5. We compute a rigorous bound for γ0 and γ1 which appear in the formula
(44) and verify that for c ∈ ci holds γ1 6= 0.
The rigorous bounds for the values of γ1 on C are
γ1(C1) ⊂ [0.014515898754816965, 157.76639522562903]
γ1(C2) ⊂ [1.1002393483255526, 151.35147664498677]
The computer assisted proof of the above took approximately 7 hours and 50
minutes on the Pentium IV 3GHz processor.
9 Implementation notes.
All the algorithms presented in this paper have been implemented in C++ by
authors and are part of the CAPD library [CAPD]. In particular, the package
implements the computation of partial derivatives of a flow with respect to
initial condition, partial derivatives of Poincaré maps for linear sections and
computations of normal forms for planar maps up to order 5.
The implementation combines the automatic and symbolic differentiation in
order to generate a coefficients in Taylor series for the solutions of the system
(14).
Our tests shows that without difficulty we can compute partial derivatives
up to order 3 for an equation in 8-dimensional phase space (which gives 1320
equations to solve) on a computer with 512MB memory. However, our current
implementation is optimized for lower dimensional problems. All the trees which
represent formulas (14) are stored in the memory of a computer. This speeds
up computations because we do not need to recompute all the multiindices,
multipointers and submultipointers in each step of the algorithm. Unfortunately,
such an implementation is memory-consuming. Therefore, higher dimensional
problems require a computer with huge memory even for C3 or C5 computations.
A Explicit formulas for third order normal forms
for a planar map
The goal of this section is to give some details about the proof of Theorem 16.
We want to present some formulas to give the reader the feeling about the
necessary computations.
Throughout this section we assume that the assumptions of Theorem 16 are
satisfied. In the neighbourhood of 0 f is given by a real, convergent power series
f(x, y) = (x1, y1)
ak−l,lx
lyk−l
bk−l,lx
lyk−l
Denote also by f : C2 → C2 a complex extension of f . Let λ, λ̄ ∈ C be
complex eigenvalues of Df(0) and v, v̄ ∈ C corresponding eigenvectors (here
bar denotes the complex conjugation). Then, using a linear substitution of the
form L = [vT , v̄T ], we can change the coordinate system such that in the new
coordinates the mapping f has the form
f(ξ, η) = (λξ + p(ξ, η), λ̄η + q(ξ, η))
p(ξ, η) =
pl,k−lξ
lηk−l
q(ξ, η) =
ql,k−lξ
lηk−l
pi,j = qj,i for i, j ≥ 0.
The last condition is a consequence of the invariance of R2 ⊂ C2 under the
complex map f . We will refer to it as the reality condition. Namely, the set
2 ⊂ C2 in the new coordinates (ξ, η) is given by ξ = η and the condition
f(R2) ⊂ R2 expressed in coordinates (ξ, η) is equivalent to (50).
Assume now, that λk 6= 1 for k = 1, . . . , 4. Then an analytic area-preserving
substitution satisfying reality condition (50)
(Φ(z, v),Ψ(z, v)) = (z1, v1)
z1 = z +
φl,k−lz
lvk−l + · · ·
v1 = v +
ψl,k−lz
lvk−l + · · ·
where
ψ2,0 = φ0,2 = −λ2p0,2(λ3 − 1)−1
ψ1,1 = φ1,1 = −p1,1(λ− 1)−1
ψ0,2 = φ2,0 = p2,0(λ
2 − λ)−1
ψ3,0 = φ0,3 = −λ3 (p0,3 + p1,1φ0,2 + 2q0,2ψ0,2) (λ4 − 1)−1
ψ2,1 = φ1,2 =
λ2 − 1 (p1,2 + 2p2,0φ0,2 + p1,1φ1,1 + p1,1ψ0,2 + 2p0,2ψ1,1)
ψ1,2 = φ2,1 = −φ2,0ψ0,2 + φ0,2ψ2,0
ψ0,3 = φ3,0 = (p3,0 + 2p2,0φ2,0 + p1,1ψ2,0) (λ
3 − λ)−1
brings f = (f1, f2) to the normal form
(z, v) → (z(α0 + α2zv), v(β0 + β2zv)) +O((zv)2)
β0 = α0 = λ
β2 = α2 = q1,2 + 2q2,0φ0,2 + q1,1φ1,1 + q1,1ψ0,2 + 2q0,2ψ1,1
Finally, let γ0 ∈ R be such that λ = α0 = eiγ0 and we compute coefficient γ1 by
From the proof given in [SM] it follows that γ1 ∈ R and the mapping f in
coordinates (z, v) has the form
f(z, v) =
zei(γ0+γ1zv), ve−i(γ0+γ1zv)
where O4 is a convergent power series with the terms of degree at least 4.
Again, the coefficients of f(z, v) satisfy reality condition (50). In order to
express this normal form in terms of real variables we make a linear substitution
z = r + is, v = r − is
and we obtain the normal form for f
f(r, s) = (r1, s1) +O4
r1 = r cos(γ0 + γ1(r
2 + s2))− s sin(γ0 + γ1(r2 + s2))
s1 = r sin(γ0 + γ1(r
2 + s2)) + s cos(γ0 + γ1(r
2 + s2))
which agrees with (44).
The formulas for higher order terms φi,j , ψi,j (and for γ2, which are not given
here) has been computed in Mathematica.
References
[A] G. Alefeld, Inclusion methods for systems of nonlinear equations - the
interval Newton method and modifications, in Topics in Validated Com-
putations, J. Herzberger (Editor), Elsevier Science B.V., 1994, pages
[BHS] H.W. Broer, G.B. Huitema and M.B. Sevryuk, Quasi-periodicity in
families of dynamical systems: order amidst chaos, Lecture Notes in
Mathematics, Vol. 1645, Springer Verlag, (1996).
[BM] M. Berz, K. Makino, New Methods for High-Dimensional Verified Qua-
drature, Reliable Computing, 5, 13-22 (1999)
[CAPD] CAPD – Computer Assisted Proofs in Dynamics group, a C++ pa-
ckage for rigorous numerics, http://capd.wsb-nlu.edu.pl.
[CC1] A. Celletti, L. Cherchia, KAM Stability for three-body problem of the
Solar System, Z. angw. Math. Phys. 57 (2006) 33-41
[CC2] A. Celletti, L. Cherchia, KAM Stability and Celestial Mechanics, Me-
moirs of the AMS, Vol 187, Num 878 (2007)
[DIK1] F. Dumortier, S. Ibáñez, and H. Kokubu, New aspects in the unfolding
of the nilpotent singularity of codimension three, Dynam. Syst. 16
(2001), 63–95.
[DIK2] F. Dumortier, S. Ibáñez, and H. Kokubu, Cocoon bifurcation in three
dimensional reversible vector fields, Nonlinearity 19 (2006), 305–328.
[GZ] Z. Galias, P. Zgliczyński, Computer assisted proof of chaos in the Lo-
renz system, Physica D, 115, 1998,165–188
[HNW] E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential
Equations I, Nonstiff Problems, Springer-Verlag, Berlin Heidelberg
1987.
[HZHT] B. Hassard, J. Zhang, S. Hastings, W. Troy, A computer proof that
the Lorenz equations have ”chaotic” solutions, Appl. Math. Letters, 7
(1994), 79–83
[KZ] T. Kapela, P. Zgliczyński, The existence of simple choreographies for
N-body problem - a computer assisted proof, Nonlinearity, 16 (2003),
1899-1918
[KWZ] H. Kokubu, D.Wilczak, P. Zgliczyński, Rigorous verifi-
cation of cocoon bifurcations in the Michelson system,
http://www.ii.uj.edu.pl/~wilczak, submitted
[La] J.S.W. Lamb, Reversing symmetries in dynamical systems, PhD The-
sis, Amsterdam University, 1994
http://capd.wsb-nlu.edu.pl
http://www.ii.uj.edu.pl/~wilczak
[Lau] Y-T Lau, The “cocoon” bifurcations in three-dimensional systems with
two fixed points, Int. Jour. Bif. Chaos, Vol.2, No.3 (1992) 543-558.
[Lo] 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., Claren-
don Press, Oxford, 1992.
[MC] C.K. McCord, Uniqueness of connecting orbits in the equation Y (3) =
Y 2 − 1, J. Math. Anal. Appl. 114, 584-592.
[Mi] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation,
Physica D, 19, (1986) 89-111.
[Mo] R.E. Moore, Interval Analysis. Prentice Hall, Englewood Cliffs, N.J.,
[MM] K. Mischaikow, M. Mrozek, Chaos in the Lorenz equations: A compu-
ter assisted proof. Part II: Details, Mathematics of Computation, 67,
(1998), 1023–1046
[MZ] M. Mrozek, P. Zgliczyński, Set arithmetic and the enclosing problem
in dynamics, Annales Pol. Math., 2000, 237–259
[NJ] N. S. Nedialkov, K. R. Jackson, An Interval Hermite – Obreschkoff
Method for Computing Rigorous Bounds on the Solution of an Initial
Value Problem for an Ordinary Differential Equation, chapter in the
book Developments in Reliable Computing, editor T. Csendes, 289-310,
Kluwer, Dordrecht, Netherlands, 1999.
[N] A. Neumeier, Interval methods for systems of equations, Cambridge
University Press, 1990.
[RNS] T. Rage, A. Neumaier, C. Schlier, Rigorous verification of chaos in a
molecular model, Phys. Rev. E 50 (1994), 2682–2688
[Ra] L.B. Rall, Automatic Differentiation: Techniques and Applications,
volume 120 of Lecture Notes in Computer Science. Springer Verlag,
Berlin, 1981
[Se] Sevryuk, M. B. Reversible systems, Lecture Notes in Mathematics,
1211. Springer-Verlag, Berlin, 1986
[SM] C.L. Siegel, J.K. Moser, Lectures on Celestial Mechanics, Springer-
Verlag Berlin Heidelberg New York, 1995.
[SK] D. Stoffer, U. Kirchgraber, Possible chaotic motion of comets in the
Sun Jupiter system - an efficient computer-assisted approach, Nonli-
nearity, 17 (2004) 281-300.
[T] W. Tucker, A Rigorous ODE solver and Smale’s 14th Problem, Foun-
dations of Computational Mathematics, (2002), Vol. 2, Num. 1, 53-117
[W] W. Walter, Differential and integral inequalities, Springer-Verlag Ber-
lin Heidelberg New York, 1970
[Wi1] D. Wilczak, The existence of Shilnikov homoclinic orbits in the Mi-
chelson system: a computer assisted proof, Found. Comp. Math. Vol.6,
No.4, 495-535, (2006).
[Wi2] D. Wilczak, Chaos in the Kuramoto–Sivashinsky equations – a compu-
ter assisted proof, J. Diff. Eqns., 194, 433-459 (2003).
[Wi3] D. Wilczak, Symmetric heteroclinic connections in the Michelson sys-
tem – a computer assisted proof, SIAM J. App. Dyn. Sys., Vol.4, No.3,
489-514 (2005).
[Wi4] D. Wilczak, http://www.ii.uj.edu.pl/~wilczak, a refference for
auxiliary materials.
[WZ] D. Wilczak, P. Zgliczyński, Heteroclinic Connections between Periodic
Orbits in Planar Restricted Circular Three Body Problem - A Compu-
ter Assisted Proof, Commun. Math. Phys. 234, 37-75 (2003).
[Z1] P. Zgliczyński, Computer assisted proof of chaos in the Hénon map and
in the Rössler equations, Nonlinearity, 1997, Vol. 10, No. 1, 243–252
[Z] P. Zgliczyński, C1-Lohner algorithm, Foundations of Computational
Mathematics, (2002) 2:429–465
[ZPer] P. Zgliczyński, Lohner Algorithm for perturbations of ODEs and diffe-
rential inclusions, http://www.ii.uj.edu.pl/~zgliczyn
http://www.ii.uj.edu.pl/~wilczak
http://www.ii.uj.edu.pl/~zgliczyn
	Introduction
	Basic definitions
	Multiindices
	Multipointers
	Equations for variations
	Cr-Lohner algorithm
	Why one needs an Cr-algorithm?
	An outline of the algorithm
	Computation of a rough enclosure for Da
	The procedure for the computation of the rough enclosure for V.
	Computation of [Vk+1]
	Composition formulas
	The procedure for computation of [Vk+1]
	Rearrangement for Va - the evaluation of Equation (32)
	Derivatives of Poincaré map
	Higher order derivatives of the Poincaré map
	Partial derivatives of tP for affine sections
	Applications.
	Area preserving maps on the plane, normal forms and KAM theorem 
	The existence of invariant tori in forced pendulum.
	Higher order normal forms.
	Application to the Michelson system
	Implementation notes.
	Explicit formulas for third order normal forms for a planar map
