$C^r$-Lohner algorithm

TL;DR

This paper introduces a rigorous computational algorithm for solving ordinary differential equations and their derivatives, enabling the proof of complex dynamical phenomena like invariant tori in specific systems.

Contribution

The paper develops a $C^r$-Lohner algorithm that computes bounds for solutions and derivatives of ODEs, extending previous methods to arbitrary order derivatives.

Findings

Proves existence of multiple invariant tori in pendulum with periodic forcing.

Establishes invariant tori in Michelson system.

Demonstrates the algorithm's effectiveness in complex dynamical systems.

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 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.

Figures

Equations

$$x^{\prime }=f(x),x\in {\mathbb{R}}^n$$

$$x^{\prime }$$

$$\frac{d}{dt}{V}_{ij}(t)$$

$$\frac{d}{dt}{H}_{ijk}(t)$$

$$x(0)=x_0,V(0)=Id,H_{ijk}(0)=0,i,j,k=1,\dots ,n.$$

$$\frac{\partial {\phi }_i}{\partial x_j}(t,x_0)$$

$$\frac{{\partial }^2{\phi }_i}{\partial x_j\partial x_k}(t,x_0)$$

$$\frac{{\partial }^3(f\circ g)}{\partial x_i\partial x_j\partial x_c}=\sum _{k,r,s=1}^n\frac{{\partial }^{3}f}{\partial x_k\partial x_r\partial x_s}\frac{\partial g_k}{\partial x_i}\frac{\partial g_r}{\partial x_j}\frac{\partial g_s}{\partial x_c}+\sum _{k=1}^n\frac{\partial f}{\partial x_k}\frac{{\partial }^3{g}_k}{\partial x_i\partial x_j\partial x_c}$$

$$+\sum _{k,r=1}^n\frac{{\partial }^{2}f}{\partial x_k\partial x_r}\left(\frac{{\partial }^2{g}_k}{\partial x_i\partial x_c}\frac{\partial g_r}{\partial x_j}+\frac{\partial g_k}{\partial x_i}\frac{{\partial }^2{g}_r}{\partial x_j\partial x_c}+\frac{{\partial }^2{g}_k}{\partial x_i\partial x_j}\frac{\partial g_r}{\partial x_c}\right)$$

$$e_i^n=(0,0,\dots ,0,\stackrel{i}{\overbrace{1}},0,\dots ,0,0).$$

$$D^{\alpha }f(t,x)=(D^{\alpha }{f}_1(t,x),\dots ,D^{\alpha }{f}_n(t,x)).$$

$${\mathcal{N}}_p^n$$

$$\mathcal{N}={\mathcal{N}}^n$$

$$\Lambda :{\mathcal{N}}_p^n\ni (a_1,\dots ,a_p)\to \sum _{i=1}^p{e}_{a_i}^n\in {\mathbb{N}}_p^n$$

$$a+b=(a_1+b_1,\dots ,a_n+b_n)\in {\mathbb{N}}_{p+q}^n.$$

$$\alpha +\beta ={\Lambda }^{-1}\left(\Lambda (\alpha )+\Lambda (\beta )\right)\in {\mathcal{N}}_{p+q}^n.$$

$$\left(a\le b\right)⟺\{\begin{array}{l}\text{either }\exists i,i\le p,i\le q,a_i<b_i\text{ and }{a}_j=b_j\text{ for }j<i\\ \text{or }p\le q\text{ and }{a}_i=b_i\text{ for }i=1,\dots ,p.\end{array}$$

$${\mathcal{N}}^p(k):=\{({\delta }_1,\dots ,{\delta }_k)\in ({\mathcal{N}}^p{)}^k:{\delta }_1\le \cdots \le {\delta }_k,{\delta }_1+\cdots +{\delta }_k=(1,2,\dots ,p)\}$$

$${\mathcal{N}}^p(1)=\{(1,2,\dots ,p)\}$$

$${\mathcal{N}}^p(p)=\{((1),(2),\dots ,(p))\}$$

$$({\delta }_1,\dots ,{\delta }_{s-1},{\delta }_s+(p),{\delta }_{s+1},\dots ,{\delta }_m)\in {\mathcal{N}}^p(m)$$

$$\begin{array}{rl}A& =\left\{({\delta }_1,\dots ,{\delta }_{k-1},(p)):({\delta }_1,\dots ,{\delta }_{k-1})\in {\mathcal{N}}^{p-1}(k-1)\right\}\\ B& =\bigcup _{s=1}^k\left\{({\delta }_1,\dots ,{\delta }_{s-1},{\delta }_s+(p),{\delta }_{s+1},\dots ,{\delta }_k):({\delta }_1,\dots ,{\delta }_k)\in {\mathcal{N}}^{p-1}(k)\right\}\end{array}$$

$$a_{\delta }:={\Lambda }^{-1}\left(\sum _{i=1}^k{e}_{a_{{\delta }_i}}^n\right)\in {\mathcal{N}}_k^n$$

$$\frac{d}{dt}{D}_a{\phi }_i=\sum _{k=1}^p\sum _{i_1,\dots ,i_k=1}^n\left(D^{e_{i_1}+\cdots +e_{i_k}}{f}_i\right)\circ \phi \sum _{({\delta }_1,\dots ,{\delta }_k)\in {\mathcal{N}}^p(k)}\prod _{j=1}^k{D}_{a_{{\delta }_j}}{\phi }_{i_j}$$

$$\frac{d}{dt}{D}_{(c)}{\phi }_i=\frac{d}{dt}\frac{\partial {\phi }_i}{\partial x_c}=\sum _{s=1}^n\frac{\partial f_i}{\partial x_s}\frac{\partial {\phi }_s}{\partial x_c}=\sum _{s=1}^n{D}^{e_s}{f}_i\cdot D_{(c)}{\phi }_s.$$

$$\frac{d}{dt}{D}_a{\phi }_i$$

$$=D_{(c)}\left(\sum _{k=1}^{p-1}\sum _{\begin{array}{c}{i}_1,\dots ,i_k=1\\ \beta :=e_{i_1}+\cdots +e_{i_k}\end{array}}^n{D}^{\beta }{f}_i\sum _{({\delta }_1,\dots ,{\delta }_k)\in {\mathcal{N}}^{p-1}(k)}\prod _{j=1}^k{D}_{b_{{\delta }_j}}{\phi }_{i_j}\right)$$

$$=\sum _{k=1}^{p-1}\sum _{\begin{array}{c}{i}_1,\dots ,i_{k+1}=1\\ \beta :=e_{i_1}+\cdots +e_{i_{k+1}}\end{array}}^n{D}^{\beta }{f}_i\cdot D_{(c)}{\phi }_{i_{k+1}}\sum _{({\delta }_1,\dots ,{\delta }_k)\in {\mathcal{N}}^{p-1}(k)}\prod _{j=1}^k{D}_{b_{{\delta }_j}}{\phi }_{i_j}$$

$$+\sum _{k=1}^{p-1}\sum _{\begin{array}{c}{i}_1,\dots ,i_k=1\\ \beta :=e_{i_1}+\cdots +e_{i_k}\end{array}}^n{D}^{\beta }{f}_i\sum _{({\delta }_1,\dots ,{\delta }_k)\in {\mathcal{N}}^{p-1}(k)}\sum _{s=1}^k{D}_{b_{{\delta }_s}+(c)}{\phi }_{i_s}\prod _{\begin{array}{c}j=1,\\ j\ne s\end{array}}^k{D}_{b_{{\delta }_j}}{\phi }_{i_j}$$

$$T_k:=\sum _{i_1,\dots ,i_k=1}^n{D}^{e_{i_1}+\cdots +e_{i_k}}{f}_i\sum _{({\delta }_1,\dots ,{\delta }_k)\in {\mathcal{N}}^p(k)}\prod _{j=1}^k{D}_{a_{{\delta }_j}}{\phi }_{i_j}$$

$$\frac{d}{dt}{D}_a{\phi }_i=\sum _{k=1}^p{T}_k$$

$$\frac{d}{dt}{D}_a{\phi }_i=\sum _{k=1}^p{S}_k,S_i=T_i,i=1,\dots ,p.$$

$$S_1$$

$$S_p$$

$$S_m$$

Full text

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 $\mathcal{C}^{r}$-Lohner algorithm. Let $r$ be a positive integer, then by $\mathcal{C}^{r}$-algorithm we will mean the routine which gives rigorous estimates for partial derivatives with respect to initial conditions up to an order $r$ and $\mathcal{C}^{r}$-computations we mean an application of an $\mathcal{C}^{r}$-algorithm.

Our main motivation for the development of $\mathcal{C}^{r}$-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 $\mathcal{C}^{0}$- and $\mathcal{C}^{1}$-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 enough. Namely, to compute the partial derivatives of the flow induced by

[TABLE]

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

[TABLE]

with the initial conditions

[TABLE]

It is well known that if by $\varphi(t,x_{0})$ we denote the (local) flow induced by (1), then

[TABLE]

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 automatic differentiation [Ra] and the algorithms for computation of compositions of multivariate Taylor series.

The proposed algorithm has been successfully applied in [HNW] to the Michelson system [Mi], where a computer assisted proof of the existence of a cocoon bifurcation was presented. Some parts of this proof required $\mathcal{C}^{2}$-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 $C^{3}$ and $C^{5}$ computations.

2 Basic definitions

To effectively deal with the formulas involving partial derivatives we will use extensively a notation of multiindices, multipointers and submultipointers throughout the paper.

As an motivation let us consider the formula for the partial derivatives of the composition of maps. Assume $g:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $f:\mathbb{R}^{n}\to\mathbb{R}$ are of class $\mathcal{C}^{3}$. We have

[TABLE]

To the operator $\frac{\partial^{3}}{\partial x_{i_{1}}\partial x_{i_{2}}\partial x_{i_{3}}}$ we can in a unique way assign a multipointer, which is a nondecreasing sequence of integers $(j_{1},j_{2},j_{3})$, such that $\{i_{1},i_{2},i_{3}\}=\{j_{1},j_{2},j_{3}\}$. A submultipointer is a multipointer, which is a part of a longer multipointer, 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 $\alpha\in\mathbb{N}^{n}$. It is another way to represent various partial derivatives. The coefficient $\alpha_{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 $\mathbb{N}$ we will denote the set of nonnegative integers, i.e. $\mathbb{N}=\{0,1,2,\ldots\}$.

Definition 1

An element $\tau\in\mathbb{N}^{n}$ will be called a multiindex.

For a sequence $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{N}^{n}$ and a vector $x=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$ we set

$|\alpha|=\alpha_{1}+\cdots+\alpha_{n}$ 2. 2.

$\alpha!=\alpha_{1}!\cdot\alpha_{2}!\cdots\alpha_{n}!$ 3. 3.

$x^{\alpha}=(x_{1}^{\alpha_{1}},\ldots,x_{n}^{\alpha_{n}})$

By $e_{i}^{n}\in\mathbb{N}^{n}$ we will denote

[TABLE]

We will drop the index $n$ (the dimension) in the symbol $e^{n}_{i}$ when it is obvious from the context.

Put $\mathbb{N}^{n}_{p}:=\{a\in\mathbb{N}^{n}:|a|=p\}$.

For $\delta=(\delta_{1},\ldots,\delta_{k})\in\mathbb{N}^{n_{1}}\times\cdots\times\mathbb{N}^{n_{k}}$ we set

$|\delta|=\sum_{i=1}^{k}|\delta_{i}|$ 2. 2.

$\delta!=\prod_{i=1}^{k}\delta_{i}!$

Let $f=(f_{1},\ldots,f_{m}):\mathbb{R}^{n}\to\mathbb{R}^{m}$ be sufficiently smooth. For $\alpha\in\mathbb{N}^{n}$ we set

$\displaystyle D^{\alpha}f_{i}=\frac{\partial^{|\alpha|}f_{i}}{\partial x_{1}^{\alpha_{1}}\cdots\partial x_{n}^{\alpha_{n}}}$ 2. 2.

$D^{\alpha}f=(D^{\alpha}f_{1},D^{\alpha}f_{2},\ldots,D^{\alpha}f_{m})$

For a function $f:\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ by $D^{\alpha}f_{i}(t,x)$ we will denote $D^{\alpha}f_{i}(t,\cdot)(x)$ and similarly

[TABLE]

This convention means that $D^{\alpha}$ always acts on $x$-variables.

2.2 Multipointers

For a fixed $n>0$ and $p>0$ we define

[TABLE]

Definition 2

An element of $\mathcal{N}^{n}$ will be called a multipointer.

Remark 3

A function

[TABLE]

is a bijection.

Let $f=(f_{1},\ldots,f_{m}):\mathbb{R}^{n}\to\mathbb{R}^{m}$ be a sufficiently smooth. For $a\in\mathcal{N}_{p}^{n}$ we set

$\displaystyle D_{a}f_{i}:=\frac{\partial^{p}f_{i}}{\partial x_{a_{1}}\ldots\partial x_{a_{p}}}$ 2. 2.

$\displaystyle D_{a}f:=(D_{a}f_{1},\ldots,D_{a}f_{m})$

For a function $f:\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ by $D_{a}f_{i}(t,x)$ we will denote $D_{a}f_{i}(t,\cdot)(x)$. In the light of the above notations $D_{\alpha}f=D^{\Lambda(\alpha)}f$.

For $a=(a_{1},a_{2},\ldots,a_{n})\in\mathbb{N}^{n}_{p}$ and $b=(b_{1},b_{2},\ldots,b_{n})\in\mathbb{N}^{n}_{q}$ we define

[TABLE]

For $\alpha\in\mathcal{N}^{n}_{p}$ and $\beta\in\mathcal{N}^{n}_{q}$ we define

[TABLE]

By $\leq$ we will denote a linear order (lexicographical order) in $\mathcal{N}$ defined in the following way. For $a\in\mathcal{N}^{n}_{p}$ and $b\in\mathcal{N}^{n}_{q}$

[TABLE]

Definition 4

For $k\leq p$ we set

[TABLE]

We will use $\mathcal{N}^{p}(k)$ extensively in the next section. Its will be used to label terms in $D^{\alpha}f_{i}(\varphi(t,x))$. Observe that for $p>0$

[TABLE]

One can construct all elements of $\mathcal{N}^{p}(k)$ using the following recursive procedure. From the definition of $\mathcal{N}^{p}(k)$ it follows that if $(\delta_{1},\ldots,\delta_{m-1})\in\mathcal{N}^{p-1}(m-1)$ then $(\delta_{1},\ldots,\delta_{m-1},(p))\in\mathcal{N}^{p}(m)$ (notice that order is preserved). Similarly, if $(\delta_{1},\ldots,\delta_{m})\in\mathcal{N}^{p-1}(m)$ then

[TABLE]

and again order of elements is preserved. Hence, for $p>2$ and $1<k<p$ we have $\mathcal{N}^{p}(k)=A\cup B$ where

[TABLE]

and the sets $A$ and $B$ are disjoint.

Another way to generate all elements of $\mathcal{N}^{p}(k)$ can be described as follows

• •

decompose the set $\{1,2,\dots,p\}$ into $k$ nonempty and disjoints sets $\Delta_{i}$, $i=1,\dots,k$

• •

we sort each $\Delta_{i}$ and permute $\Delta_{i}$’s to obtain $\min(\Delta_{1})<\min(\Delta_{2})<\dots<\min(\Delta_{k})$

• •

we define $\delta_{i}$ to be an ordered set consisting of all elements of $\Delta_{i}$ for $i=1,\dots,k$

Definition 5

For an arbitrary $a\in\mathcal{N}_{p}^{n}$ and $\delta\in\mathcal{N}^{p}_{k}$ such that $k\leq p$ we define a submultipointer $a_{\delta}\in\mathcal{N}^{n}_{k}$ by $(a_{\delta})_{i}=a_{\delta_{i}}$ for $i=1,\ldots,k$, which can be expressed using $\Lambda$ as follows

[TABLE]

3 Equations for variations

Consider an ODE $x^{\prime}=f(x)$ where $f$ is $\mathcal{C}^{K+1}$. Let $\varphi:\mathbb{R}\times\mathbb{R}^{n}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\mathbb{R}^{n}$ be a local dynamical system induced by $x^{\prime}=f(x)$. It is well known, that $\varphi\in\mathcal{C}^{K}$ and one can derive the equations for partial derivatives of $\varphi$ by differentiating equation $\frac{\partial\varphi}{\partial t}(t,x)=f(\varphi(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 (1).

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\in\mathcal{C}^{r+1}$ and let $\varphi:\mathbb{R}\times\mathbb{R}^{n}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\mathbb{R}^{n}$ be a local dynamical system induced by $x^{\prime}=f(x)$. Then for $a\in\mathcal{N}^{n}_{p}$ such that $p\leq r$ holds

[TABLE]

for $i=1,\ldots,n$.

Proof: In the proof the functions $D^{e_{i_{1}}+\cdots+e_{i_{k}}}f_{i}$ are always evaluated at $\varphi(t,x)$, and various partial derivatives of $\varphi$ are always evaluated at $(t,x)$, therefore 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\in\{1,\ldots,n\}$ and (9) becomes

[TABLE]

Assume (9) holds true for $p-1$, $p>1$. Let us fix $a\in\mathcal{N}^{n}_{p}$. We have $a=b+(c)$, where $b=(a_{1},\ldots,a_{p-1})\in\mathcal{N}^{n}_{p-1}$ and $c=a_{p}$. Since (9) is satisfied for $p-1$, therefore we have

[TABLE]

For $k=1,\dots,p$ we set

[TABLE]

Now our goal is to prove that:

[TABLE]

Our strategy of proof is as follows. We will define $S_{1},\dots,S_{p}$, such that

[TABLE]

We set

[TABLE]

For $m=2,3,\ldots,p-1$ we set

[TABLE]

It remains to show that $S_{i}=T_{i}$ for $i=1,\dots,p$. Consider first $i=1$. Recall that $\mathcal{N}^{p-1}(1)=\{(1,2,\dots,p-1)\}$, hence

[TABLE]

Therefore

[TABLE]

Consider now $i=p$. For an arbitrary $s>0$ $\mathcal{N}^{s}(s)$ contains only one element $((1),(2),\ldots,(s))$. Therefore we obtain

[TABLE]

Since $a=b+(c)$, where $c=(a_{p})$, hence

[TABLE]

Consider now $m=2,3,\ldots,p-1$. We have

[TABLE]

Using decomposition $\mathcal{N}^{p}(m)=A\cup B$ as in (8) we obtain

[TABLE]

We have shown that $T_{i}=S_{i}$ for $i=1,\dots,p$. This finishes the proof.

4 $\mathcal{C}^{r}$-Lohner algorithm

4.1 Why one needs an $\mathcal{C}^{r}$-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 $\mathcal{C}^{r}$-computations the number of equations to solve is equal to $n\begin{pmatrix}n+r\\ n\end{pmatrix}$ 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 $\mathcal{C}^{0}$, are $\mathcal{C}^{1}$ ’internally’(or higher for Taylor models), because they solve non-rigorously equations for ($\frac{\partial\varphi}{\partial x}$) - 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 $\mathcal{C}^{r}$-algorithm is given by $n\begin{pmatrix}n+r\\ n\end{pmatrix}$ 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 $D_{a}\varphi$. This feature is not present in the blind application of $\mathcal{C}^{0}$ 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 $D_{a}\varphi$-variables.

In [Z] a $\mathcal{C}^{1}$-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\leq K$ and consider the following system of differential equations

[TABLE]

for all $a\in\mathcal{N}^{n}_{d}$, $d=1,\ldots,r$.

Our goal is to present an algorithm for computing a rigorous bound for the solution of (14) with a set of initial conditions

[TABLE]

In the sequel we will use the following notations:

• •

if a solution of system (14) is defined for $t>0$ and some $x_{0}\in\mathbb{R}^{n}$, then for $a\in\mathcal{N}$ by $V_{a}(t,x_{0})$ we denote $D_{a}\varphi(t,x_{0})$

• •

for $[x_{0}]\subset\mathbb{R}^{n}$ by $[V_{a}(t,[x_{0}])]$ we will denote a set for which we have $V_{a}(t,[x_{0}])\subset[V_{a}(t,[x_{0}])]$. This set is obtained using an rigorous numerical routine described below.

The $\mathcal{C}^{r}$-Lohner algorithm is a modification of $\mathcal{C}^{1}$-Lohner algorithm [Z]. One step of $\mathcal{C}^{r}$-Lohner is a shift along the trajectory of the system (14) with the following input and output data

Input data:

• •

$t_{k}$ - a current time,

• •

$h_{k}$ - a time step,

• •

$[x_{k}]\subset{\hbox{{\bbc R}}}^{n}$, such that $\varphi(t_{k},[x_{0}])\subset[x_{k}]$,

• •

$[V_{k,a}]=[V_{k,a}(t_{k},[x_{0}])]\ \subset{\hbox{{\bbc R}}}^{n}$, such that $D_{a}\varphi(t_{k},[x_{0}])\subset[V_{k,a}]$ for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$.

Output data:

• •

$t_{k+1}=t_{k}+h_{k}$ - a new current time,

• •

$[x_{k+1}]\subset{\hbox{{\bbc R}}}^{n}$, such that $\varphi(t_{k+1},[x_{0}])\subset[x_{k+1}]$,

• •

$[V_{k+1,a}]=[V_{k+1,a}(t_{k+1},[x_{0}])]\ \subset{\hbox{{\bbc R}}}^{n}$, such that $D_{a}\varphi(t_{k+1},[x_{0}])\subset[V_{k+1,a}]$ for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$.

We will often skip the arguments of $V_{k,a}$ when they are obvious from the context.

The values of $[x_{k+1}]$ and $[V_{k+1,a}]$, $a\in\mathcal{N}^{n}_{1}$ are computed using one step $\mathcal{C}^{1}$-Lohner algorithm. After it is done, we perform the following operations to compute $[V_{k+1,a}]$ for $a\in\mathcal{N}^{n}_{2}\cup\ldots\cup\mathcal{N}^{n}_{r}$

1.

Find a rough enclosure for $D_{a}\varphi([0,h_{k}],[x_{k}])$.

2.

Compute $[V_{k+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 $D_{a}\varphi$

For a fixed multipointer $a\in\mathcal{N}^{n}_{d}$ Equation (14) can be written as follows

[TABLE]

where

[TABLE]

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 $\mu(A)$ is defined as a limit

[TABLE]

where $\|\cdot\|$ is a given matrix norm.

The formulas for the logarithmic norm of a real matrix in the most frequently used norms are (see [HNW])

for $\|x\|_{1}=\sum_{i}|x_{i}|$, $\mu(A)=\max_{j}(a_{jj}+\sum_{i\neq j}|a_{ij}|)$ 2. 2.

for m $\|x\|_{2}=\sqrt{\sum_{i}|x_{i}|^{2}}$, $\mu(A)$ is equal to the largest eigenvalue of $(A+A^{T})/2$ 3. 3.

for $\|x\|_{\infty}=\max_{i}|x_{i}|$, $\mu(A)=\max_{i}(a_{ii}+\sum_{j\neq i}|a_{ij}|)$

In order to find bounds for $D_{a}\varphi$ we use the following theorem [HNW, Thm. I.10.6]

Theorem 8

Let $x(t)$ be a solution of a differential equation

[TABLE]

Let $\nu(t)$ be a piecewise differentiable function with values in ${\hbox{{\bbc R}}}^{n}$. Assume that

[TABLE]

where by $\mu(A)$, we denote a logarithmic norm of a square matrix $A\in{\hbox{{\bbc R}}}^{n\times n}$.

Then for $t\geq t_{0}$ we have

[TABLE]

with $L(t)=\int_{t_{0}}^{t}l(\tau)d\tau$.

We apply the above theorem to Equation (16) to obtain

Lemma 9

Let us fix $x\in\hbox{{\bbc R}}^{n}$. Assume that $|B_{a}(t,x)|\leq\delta(t)$ and $\mu(A(t,x))\leq l(t)$, then for $t>t_{0}$

[TABLE]

with $L(t)=\int_{t_{0}}^{t}l(\tau)d\tau$.

Proof: Consider Equation (16) and a homogenous problem for (16)

[TABLE]

Using Theorem 8 we can estimate the difference between any solution of (21), $w$, and a solution of (16), denoted by $D_{a}\varphi$.

[TABLE]

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 $\delta$ and $l$, hence we will use the following

Lemma 10

Assume that $|B_{a}(t,x)|\leq\delta$ and $\mu(A(t,x))\leq l$ for $t\in[0,h]$ then for $t\in[0,h]$ we have

[TABLE]

or

[TABLE]

5.1 The procedure for the computation of the rough enclosure for $V$.

The procedure for the computing of the rough enclosure is iterative, which means that given a rough enclosure for $\varphi([0,h_{k}],[x_{k}])$ and rough enclosures $D_{a}\varphi([0,h_{k}],[x_{k}])$ for all $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{p}$ we are able to compute the rough enclosure for $D_{a}\varphi([0,h_{k}],[x_{k}])$ for $a\in\mathcal{N}^{n}_{p+1}$.

The procedures for computation of the rough enclosures of $\varphi([0,h_{k}],[x_{k}])$ and $D_{a}\varphi([0,h_{k}],[x_{k}])$ for $a\in\mathcal{N}^{n}_{1}$ has been given in [Z]. Below we present an algorithm for computing $[E_{a}]$ for $a\in\mathcal{N}^{n}_{2}\cup\ldots\cup\mathcal{N}^{n}_{r}$.

Input parameters:

• •

$h_{k}$ - a time step,

• •

$[x_{k}]\subset{\hbox{{\bbc R}}}^{n}$ - the current value of $x=\varphi(t_{k},[x_{0}])$,

• •

$[E_{0}]\subset{\hbox{{\bbc R}}}^{n}$ - a compact and convex such that $\varphi([0,h_{k}],[x_{k}])\subset[E_{0}]$

• •

$[E_{a}]\subset{\hbox{{\bbc R}}}^{n}$, $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{p}$ such that $D_{a}\varphi([0,h_{k}],[x_{k}])\subset[E_{a}]$ for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{p}$.

Output:

• •

$[E_{a}]\subset{\hbox{{\bbc R}}}^{n}$, $a\in\mathcal{N}^{n}_{p+1}$ such that

[TABLE]

Before we present an algorithm let us observe that for a fixed $a\in\mathcal{N}^{n}_{p+1}$, $B_{a}$ defined in (17) could be seen as a multivariate function of $t$, $x$ and $V_{b}=D_{b}\varphi$ for $b\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{p}$. More precisely, put $m_{p}:=\sharp\left\{\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{p}\right\}$, where $\sharp$ stands for number of elements of a set. Recall that, we have defined by (6) a linear order in $\mathcal{N}^{n}$. Hence, there is a unique sequence of multipointers $b_{1},\ldots,b_{m_{p}}$, such that $b_{i}\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{p}$ for $i=1,\ldots,m_{p}$, $b_{1}\leq b_{2}\leq\cdots\leq b_{m_{p}}$ and $b_{i}\neq b_{j}$ for $i\neq j$.

Let us define

[TABLE]

by

[TABLE]

and

[TABLE]

Algorithm:

To compute $[E_{a}]$ for $a\in\mathcal{N}^{n}_{p+1}$ we proceed as follows

1.

Find $l\geq\left(\max_{x\in[E_{0}]}\mu\left(Df(x)\right)\right)$.

2.

Compute $\delta_{a}\geq\max\|\tilde{B}_{a}\|$, i.e.

[TABLE]

For example, if $a=(j,c)\in\mathcal{N}^{n}_{2}$, then $\delta_{a}$ should be such that

[TABLE]

3.

Define $[E_{a}]_{i}=[-1,1]\delta_{a}\frac{e^{lt}-1}{l}$, for $i=1,\ldots,n$, where $[E_{a}]_{i}$ denotes $i$-th coordinate of $[E_{a}]$.

One can refine the obtained enclosure by

[TABLE]

Indeed, for $i=1,\ldots,n$, $t\in[0,h_{k}]$ and $x_{0}\in[E_{0}]$ we have

[TABLE]

for some $\theta_{i}\in[0,t]\subset[0,h_{k}]$. In the above we have used the fact that

[TABLE]

Since $\varphi(\theta_{i},x_{0})\in[E_{0}]$ and $D_{b_{j}}\varphi(\theta_{i},x_{0})\in[E_{b_{j}}]$ for $j=1,\ldots,m_{p}$ we get

[TABLE]

6 Computation of $[V_{k+1}]$

6.1 Composition formulas

For any $p$-times continuously differentiable functions $f,g:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $a\in\mathcal{N}^{n}_{p}$ we have

[TABLE]

We can apply the above formula to $f=\varphi(h_{k},\cdot)$ and $g=\varphi(t_{k},\cdot)$ to obtain

[TABLE]

for all $x_{0}\in[x_{0}]$. Using notations $[V_{k+1,a}]:=[V_{a}(t_{k}+h_{k},[x_{0}])]$ and $[V_{k,a}]=[V_{a}(t_{k},[x_{0}])]$ we can rewrite the above equation as

[TABLE]

where $\Lambda$ is defined by (5).

6.2 The procedure for computation of $[V_{k+1}]$

We introduce new parameters $o_{d}$ - the order of the Taylor method used in computations of $V_{a}$ for $a\in\mathcal{N}^{n}_{d}$. It makes sense to take $o_{1}\geq o_{2}\geq\cdots\geq o_{r}$.

Input parameters:

• •

$h_{k}$ - a time step,

• •

$[x_{k}]\subset{\hbox{{\bbc R}}}^{n}$ - the current value of $x=\varphi(t_{k},[x_{0}])$,

• •

$[V_{k,a}]\subset{\hbox{{\bbc R}}}^{n}$ - a current value of $V_{k,a}(t_{k},[x_{0}])$, for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$

• •

$[E_{0}]\subset{\hbox{{\bbc R}}}^{n}$ compact and convex, such that $\varphi([0,h_{k}],[x_{k}])\subset[E_{0}]$ - a rough enclosure for $[x_{k}]$,

• •

$[E_{a}]\subset{\hbox{{\bbc R}}}^{n}$, compact and convex, such that $D_{a}\varphi([0,h_{k}],[x_{k}])\subset[E_{a}]$, for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$.

Output: $[V_{k+1,a}]\subset{\hbox{{\bbc R}}}^{n}$, such that

[TABLE]

for $x_{0}\in[x_{0}]$ and $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$.

Algorithm: We compute $[V_{k+1}]$ as follows

1.

Computation of $V_{a}(h_{k},[x_{k}])$ using Taylor method for Equation (14), i.e. for $a\in\mathcal{N}^{n}_{p}$ we compute

[TABLE]

where $V_{b_{i}}=0$ for $b_{i}\in\mathcal{N}^{n}_{2}\cup\ldots\cup\mathcal{N}^{n}_{p-1}$ and $V_{(j)}=e^{n}_{j}$ for $j=1,\ldots,n$. Observe that

[TABLE]

Indeed, using Taylor series expansion we obtain that for $x_{k}\in[x_{k}]$ and $j=1,\ldots,n$ holds

[TABLE]

for some $\theta_{i}\in[0,h_{k}]$. Observe, that

[TABLE]

Using $\varphi(\theta_{i},x_{k})\in[E_{0}]$ and $V_{b_{s}}(\theta_{i},x_{k})\in[E_{b_{s}}]$ for $s=1,\ldots,m_{p-1}$ we obtain our assertion.

2.

The composition. Put

[TABLE]

Using (28) for $a\in\mathcal{N}^{n}_{p}$ we have

[TABLE]

where

[TABLE]

In our implementation of the algorithm we use the symbolic differentiation to obtain formulae for $D_{a}f$. Next, using the automatic differentiation we compute $\frac{d^{i}}{dt^{i}}F_{a}(t,x,V_{b_{1}}(t,x),\ldots,V_{b_{m_{p-1}}}(t,x))|_{t=0}$ which appear in (1.).

6.3 Rearrangement for $V_{a}$ - the evaluation of Equation

(32)

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 $\mathcal{C}^{1}$-method (see [Z, Section 3]). Moreover, all vectors $V_{k,a}$, for $a\in\mathcal{N}^{n}_{1}\cup\ldots\mathcal{N}^{n}_{r}$ ’propagate’ by the same $[J_{k}]$ as did the variational part in [Z], hence it makes sense the same approach.

To be more precise, each set $[V_{k,a}]$, for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$ is represented in the following form

[TABLE]

where $[B_{k}]$ is interval matrix, $C_{k}$ is point matrix, $v_{k,a}$ is a point vector and $r_{k,a},q_{k,a}$ are interval vectors. Observe that $[B_{k}]$ and $C_{k}$ are independent of $a$.

In the sequel we will drop index $a$. Equation (32) leads to

[TABLE]

Let $m([z])$ denotes a center of an interval object, i.e. $[z]$ is interval vector or interval matrix and $\Delta([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([J_{k}])[B_{k}]$.

Let

[TABLE]

Then we rearrange formula (34) as follows

[TABLE]

Summarizing, we can use the following data structure to represent $\varphi(t_{k},[x_{0}])$ and $D_{a}\varphi(t_{k},[x_{0}])$, for $a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}$

type CnSet = record

$v_{0},r_{0},q_{0}$: IntervalVector;

$C_{0},B_{0},C,B$ : IntervalMatrix;

$\{v_{a},r_{a},q_{a}:\text{ IntervalVector}\}_{a\in\mathcal{N}^{n}_{1}\cup\ldots\cup\mathcal{N}^{n}_{r}}$

end;

The set $\varphi(t_{k},[x_{0}])$ is represented as $v_{0}+B_{0}r_{0}+C_{0}q_{0}$, the partial derivatives $D_{a}\varphi(t_{k},[x_{0}])$ are represented as $v_{a}+Br_{a}+Cq_{a}$. The matrices $B,C$ are common for all partial derivatives.

Notice, that if we start the $\mathcal{C}^{r}$ 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 $q_{a},r_{a}$ are set to zero.

If the interval vectors $r_{a}$ 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 $D_{a}\varphi$) and reset $r_{a}$ in the following way

[TABLE]

A similar change of the Lipshitz part may be done when vectors $r_{a}$ become thick in comparison to $q_{a}$.

7 Derivatives of Poincaré map

Consider a differential equation

[TABLE]

Let $\varphi:{\hbox{{\bbc R}}}\times{\hbox{{\bbc R}}}^{n}\to{\hbox{{\bbc R}}}^{n}$ be a (local) dynamical system induced by (36). Let $\alpha:{\hbox{{\bbc R}}}^{n}\to{\hbox{{\bbc R}}}$ be $\mathcal{C}^{1}$-map. Put $\Pi=\{x\ |\ \alpha(x)=C\}$.

Definition 11

We will say that $\Pi$ is a local section for the vector field $f$ at $y_{0}\in\Pi$ if

[TABLE]

Assume $x_{0}\in\mathbb{R}^{n}$ and $t_{0}\in\mathbb{R}$ are such that $\Pi$ is a local section at $\varphi(t_{0},x_{0})$. Consider an implicit equation

[TABLE]

It follows easily from (37) and from the implicit function theorem that there exists a uniquely defined $t_{P}:\mathbb{R}^{n}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\mathbb{R}$ in a neighborhood of $x_{0}$, such that $t_{P}(x_{0})=t_{0}$. The function $t_{P}$ is as smooth as the flow $\varphi$. We will refer to $t_{P}$ as to the Poincare return time to section $\Pi$.

We define a Poincaré map $P:\mathbb{R}^{n}\supset\mathrm{dom}\,(t_{P})\to\mathbb{R}^{n}$ by

[TABLE]

Usually the Poincaré map is defined as a map $P:\Pi_{1}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\Pi_{2}$, where $\Pi_{1},\Pi_{2}$ are local sections in $\mathbb{R}^{n}$. The approach taken here, i.e. treating the Poincaré map as map $P:\mathbb{R}^{n}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\mathbb{R}^{n}$ 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 $\frac{\partial P_{i}}{\partial x_{j}}$ and we obtain

[TABLE]

We need $\frac{\partial t_{P}}{\partial x_{j}}$. We differentiate (38) to obtain

[TABLE]

Hence

[TABLE]

7.1 Higher order derivatives of the Poincaré map

To make formulas transparent we will drop arguments of functions in this section, but reader should be aware that for $t_{P}$ and its partial derivatives the argument is $x$, for $\varphi$ and $D_{a}\varphi$ the argument is always the pair $(t_{P}(x),x)$.

From (40) we obtain

[TABLE]

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\in\mathcal{N}^{n}_{p}$ we have

[TABLE]

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\geq 1$ and fix $a\in\mathcal{N}^{n}_{p+1}$. Our goal is to show that

[TABLE]

where

[TABLE]

Write $a=\beta+\gamma$, where $\beta\in\mathcal{N}^{n}_{p}$ and $\gamma=(a_{p+1})\in\mathcal{N}^{n}_{1}$. From the induction assumption we have

[TABLE]

where

[TABLE]

Obviously $R_{1}=S_{1}$. We will show that $R_{2}=S_{3}+S_{6}+S_{7}$ and $R_{3}=S_{2}+S_{4}+S_{5}+S_{8}+S_{9}+S_{10}$.

Denote by $R_{i,k}$, $i=2,3$ a part of sum $R_{i}$ with fixed $k=2,\ldots,p+1$. Similarly, let us denote by $S_{i,k}$ a part of sum $S_{i}$, $i=5,\ldots,10$, for $k=2,\ldots,p$.

Using decomposition of $\mathcal{N}^{p+1}(2)$ as in (8) we obtain that $R_{2,2}=S_{3}+S_{7,2}$. Similarly, using (8) we observe that $R_{2,k}=S_{6,k-1}+S_{7,k}$ for $k=3,\ldots,p$. Finally, since $\mathcal{N}^{p+1}(p+1)=\{((1),(2),\ldots,(p+1))\}$ and $\gamma=(a_{p+1})$ we find that $R_{2,p+1}=S_{6,p}$. This shows that $R_{2}=S_{3}+S_{6}+S_{7}$.

It remains to show that $R_{3}=S_{2}+S_{4}+S_{5}+S_{8}+S_{9}+S_{10}$. We will classify possible terms by the fact, where $p+1$ appears in $\delta_{i}$, $i=1,\dots,k$ and how this $\delta_{i}$ enters in $R_{3}$ as $\delta_{s}$ or $\delta_{j}$. There are four cases

$\delta_{s}=(p+1)$ 2. 2.

$\delta_{j}=(p+1)$ 3. 3.

$p+1\in\delta_{s}$, $|\delta_{s}|\geq 2$ 4. 4.

$p+1\in\delta_{j}$, $|\delta_{j}|\geq 2$

Let us fix $k=2$. Let $(\delta_{1},\delta_{2})\in\mathcal{N}^{p+1}(2)$. The term for case 1 is $S_{4}$, for case 2 is $S_{2}$, case 3 is $S_{8,2}$ and case 4 is $S_{10,2}$. Hence, $R_{3,2}=S_{2}+S_{4}+S_{8,2}+S_{10,2}$.

For $k=3,\ldots,p$ and fixed $(\delta_{1},\ldots,\delta_{k})\in\mathcal{N}^{p+1}(k)$ we have: case 1 is given by $S_{5,k-1}$, case 2 by $S_{9,k-1}$, case 3 by $S_{8,k}$ and case 4 by $S_{10,k}$ Hence, for $k=3,\ldots,p$ we have $R_{3,k}=S_{5,k-1}+S_{9,k-1}+S_{8,k}+S_{10,k}$.

Finally, for $k=p+1$ we observe, that $R_{3,p+1}=S_{5,p}+S_{9,p}$. Indeed, in this case $(\delta_{1},\ldots,\delta_{p+1})=((1),(2),\ldots,(p+1))$. Hence, either for $\delta_{s}=\gamma$ we have term $S_{5,p}$ and $\delta_{s}\neq\gamma$ we have $S_{9,p}$.

We have showed that $R_{3}=S_{2}+S_{4}+S_{5}+S_{8}+S_{9}+S_{10}$ and the proof is finished.

Hence, if we know all the partial derivatives of $t_{P}$ up order $p$ we can compute the partial derivatives of the Poincaré map up the same order. In next subsection we show how to compute partial derivatives of $t_{P}$ for affine sections.

7.2 Partial derivatives of $t_{P}$ for affine sections

Assume $\alpha:\mathbb{R}^{n}\to\mathbb{R}$ is an affine map given by

[TABLE]

This is a quite restrictive assumption about sections, but it leads to relatively simple formulas for $D_{a}t_{P}$ and it is sufficient for the applications we have in mind.

Lemma 13

For a multipointer $a\in\mathcal{N}^{n}_{p}$ holds

[TABLE]

Proof: The proof is a direct consequence of Lemma 12 and (38). Since $\alpha$ is affine, by differentiating of $\alpha(P(x))=C$ we get $\left\langle\nabla\alpha|D_{a}P\right\rangle=0$. Using formula (43) for $D_{a}P$ we obtain our assertion.

Fix $[x]\subset\mathbb{R}^{n}$ and assume we have a rigorous bound for $t_{P}([x])\in[t_{1},t_{2}]$ (see [Z, Section 6] for more details on this). Lemmas 13 and 12 show that given rigorous bounds for the partial derivatives $D_{a}\varphi([t_{1},t_{2}],[x])$ and $\frac{\partial^{k}}{\partial t^{k}}D_{a}\varphi([t_{1},t_{2}],[x])$ up to some order $p$ we can compute recursively rigorous bounds for the partial derivatives of $t_{P}([x])$ and $P([x])$ up to the same order. Notice, that $\frac{\partial^{k}}{\partial t^{k}}D_{a}\varphi$ are given by Taylor coefficients of the solution of (14) with initial conditions $P([x])$ for $\mathcal{C}^{0}$ part and $D_{a}\varphi(t_{P}(x),[x])$ for equations for variations. Hence, these coefficients 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:\mathbb{R}^{2}\to\mathbb{R}^{2}$ be a smooth area preserving map, such that $f(p)=p$. Let $\lambda$ and $\mu$ be eigenvalues of $df(p)$. Following [SM] we will call the point $p$

• •

hyperbolic* if $\lambda,\mu\in\mathbb{R}$ and $\lambda\neq\mu$,*

• •

elliptic* if $\lambda=\overline{\mu}$ and $\lambda\neq\mu$,*

• •

parabolic* if $\lambda=\mu$.*

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:\mathbb{R}^{2}\to\mathbb{R}^{2}$, $f(r,s)=(r_{1},s_{1})$ where

[TABLE]

and $O_{2l+2}$ denotes convergent power series in $r,s$ with terms of order greater than $2l+1$, only.

If at least one of $\gamma_{1},\ldots,\gamma_{l}$ is not zero then the origin is a stable fixed point for map $f$. Moreover, in any neighborhood $U$ of point [math] 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:\mathbb{R}^{2}\to\mathbb{R}^{2}$ such that $f(0)=0$. Let $\lambda,\bar{\lambda}$ be complex eigenvalues of $Df(0)$, such that $|\lambda|=|\bar{\lambda}|=1$. If $\lambda^{k}\neq 1$ for $k=1,\ldots,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 $\gamma_{0},\ldots,\gamma_{l}$ in (44). An explicit formula for the coefficient $\gamma_{1}$ in the above normal form is given in Appendix A.

8.2 The existence of invariant tori in forced pendulum.

Consider an equation

[TABLE]

Observe that (45) is hamiltonian.

Let us denote by $P_{\omega}:\mathbb{R}^{2}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\mathbb{R}^{2}$ the Poincaré map for Equation (45) with a parameter $\omega$, i.e. $P_{\omega}=\varphi(2\pi/\omega,\cdot)$, where $\varphi:\mathbb{R}\times\mathbb{R}^{2}{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\mathbb{R}^{2}$ 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

[TABLE]

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_{\omega}$ corresponds and invariant 2-torus for (45).

Consider a set of parameter values

[TABLE]

The following lemma was proved with computer assistance

Lemma 17

For all parameter values $\omega\in\Omega$ there exists an elliptic fixed point $x_{\omega}\in\mathbb{R}^{2}$ for $P_{\omega}$. Moreover, there exists an area-preserving substitution such that in the new coordinates the map $f_{\omega}(x)=P_{\omega}(x+x_{\omega})-x_{\omega}$ has the form (44) with $l=1$ and $\gamma_{1}\neq 0$.

Before we give the proof, let us briefly comment about the choice of the parameter set $\Omega$. For parameter values slightly lower than $2$ we observe the parabolic case, i.e. there exists a parameter value $\omega_{1}$ for which eigenvalues of the derivative of $P_{\omega_{1}}$ are equal to $-1$. In two gaps in $\Omega$ below $3$ and $4$ we have resonances of low order. Namely, we have parameter values with an elliptic fixed with eigenvalues to $e^{\pm 2\pi/3}=\frac{-1}{2}\pm\frac{\sqrt{3}}{2}i$ and $e^{\pm i\pi/2}=\pm i$, respectively. Clearly, in a computer assisted proof we need to exclude a small interval around those parameters. For $\omega>4$ it seems that the interval $\Omega_{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 $\Omega$ by $9910$ nonequal subintervals $\omega_{i}$. Diameters of $\omega_{i}$’s were relatively large for values far away from the parabolic cases and very small close to them. For a fixed subinterval $\omega_{i}$ we proceed as follows

Let $\bar{\omega}$ denote an approximate center of the interval $\omega_{i}$. We find an approximate fixed point for $P_{\bar{\omega}}$ using the standard nonrigorous Newton method. Let us denote such a point by $x_{i}$. 2. 2.

We define a box centered at $x_{i}$, i.e we set $v_{i}:=x_{i}+[-\varepsilon_{i},\varepsilon_{i}]^{2}$, where $\varepsilon_{i}>0$ depends on subinterval $\omega_{i}$ - the values we used are from the interval $[5\cdot 10^{-6},3\cdot 10^{-3}]$, depending on whether $x_{i}$ close to parameter values corresponding to parabolic cases. 3. 3.

Using the $\mathcal{C}^{1}$-Lohner algorithm we compute the Interval Newton operator [Mo, N, A] $N_{i}:=N(P_{\omega_{i}}-\mathrm{Id},x_{i},v_{i})$ and verify that $N_{i}\subset{\rm int}\,v_{i}$. This proves that for all $\omega\in\omega_{i}$ there exists a unique fixed point $x_{\omega}\in N_{i}$ for $P_{\omega}$. 4. 4.

Using the $\mathcal{C}^{3}$-Lohner algorithm we compute a rigorous bound for $P_{\omega_{i}}(N_{i})$ and $D^{\alpha}P_{\omega_{i}}(N_{i})$, $\alpha\in\mathbb{N}^{2}_{1}\cup\mathbb{N}^{2}_{2}\cup\mathbb{N}^{2}_{3}$. Hence, we obtain a rigorous bound for the coefficients in

[TABLE] 5. 5.

We show that an arbitrary matrix $M\in DP_{\omega_{i}}(N_{i})$ has a pair of complex eigenvalues $\lambda,\bar{\lambda}$ which satisfy $\lambda^{k}\neq 1$ for $k=1,\ldots,4$. From Theorem 16 it follows there exists an area-preserving substitution such that in the new coordinates the map $f_{\omega}$ for $\omega\in\omega_{i}$ has the form (44) with $l=1$. 6. 6.

We compute a rigorous bound for $\gamma_{0}$ and $\gamma_{1}$ which appear in the formula (44) and verify that for $\omega\in\omega_{i}$ holds $\gamma_{1}\neq 0$.

The rigorous bounds for the values of $\gamma_{1}$ on $\Omega$ are

[TABLE]

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 $\omega\in\Omega$ there exists an elliptic fixed point $x_{\omega}\in\mathbb{R}^{2}$ for $P_{\omega}$. Moreover, any neighborhood of point $x_{\omega}$ contains an invariant curve for $P_{\omega}$ around $x_{\omega}$.

8.3 Higher order normal forms.

In the previous section it was shown that $\mathcal{C}^{3}$ computations are sufficient to prove that for (45) a family of invariant tori exists. However, it may happen that the coefficient $\gamma_{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,

[TABLE]

Theorem 19

Let $P_{\omega}$ be the Poincaré map for (47). For all parameter values $\omega\in\Omega_{*}=[2.9957694795,2.9957694796]$ there exists an elliptic fixed point $x_{\omega}\in\mathbb{R}^{2}$ for $P_{\omega}$. Moreover, any neighbourhood of point $x_{\omega}$ contains an invariant curve for $P_{\omega}$ around $x_{\omega}$.

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

[TABLE]

We set $v=x+3\cdot 10^{-5}([-1,1]\times[-1,1])$. Using the $\mathcal{C}^{1}$-Lohner algorithm we compute the Interval Newton Operator of $P_{\omega}-\mathrm{Id}$ on $v$ and we obtain that for all $\omega\in\Omega_{*}$, $N=N(P_{\omega}-\mathrm{Id},\mathrm{center}(v),v)\subset(N_{1},N_{2})$, where

[TABLE]

Since $N\subset v$ we conclude that for all $\omega\in\Omega_{*}$ there exists a unique fixed point $x_{\omega}\in N$ for the Poincaré map.

Using $\mathcal{C}^{5}$-Lohner algorithm we compute a rigorous bound for $P_{\Omega_{*}}(N)$ and $D^{\alpha}P_{\Omega_{*}}(N)$, $\alpha\in\mathbb{N}^{2}_{1}\cup\ldots\cup\mathbb{N}^{2}_{5}$. Hence, we obtain a rigorous bound for the coefficients in

[TABLE]

We show that an arbitrary matrix $M\in DP_{\Omega_{*}}(N)$ has a pair of complex eigenvalues $\lambda,\bar{\lambda}$ which satisfy $\lambda^{k}\neq 1$ for $k=1,\ldots,6$. From Theorem 16 it follows there exists an area-preserving substitution such that in the new coordinates the map $f_{\omega}$ for $\omega\in\Omega_{*}$