$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_{*}$ has the form (44) with $l=2$.

Next, we compute a rigorous bound for $\gamma_{1}$ and $\gamma_{2}$ which appear in the formula (44) and we get

[TABLE]

Since for $\omega\in\Omega_{*}$, $\gamma_{2}(\omega)\neq 0$ the assertion follows from Theorem 15.

The main observation which makes this example interesting is that there exists $\omega_{*}\in\Omega_{*}$ for which $\gamma_{1}(\omega_{*})=0$ and we cannot conclude the existence of invariant tori for all $\omega\in\Omega_{*}$ from $\mathcal{C}^{3}$ computations. To be more precise, we computed the coefficient $\gamma_{1}$ for the parameter values $\omega_{1}=\min\Omega_{*}$ and $\omega_{2}=\max\Omega_{*}$ and we get

[TABLE]

Since $\gamma_{1}$ exists for all $\omega\in\Omega_{*}$ and depends continuously on $\omega$ we conclude, that $\gamma_{1}(\omega_{*})=0$ for some $\omega_{*}\in\Omega_{*}$.

8.4 Application to the Michelson system

The existence of an invariant curve for a planar map $f:\mathbb{R}^{2}\to\mathbb{R}^{2}$ can be proven without assumption that $f$ is measure preserving. The key assumption in the proof given in [SM] is that any curve $\gamma$ around an elliptic point intersect its image under $f$, i.e. $f(\gamma)\cap\gamma\neq\emptyset$. Such a situation is also observed in reversible planar map around an symmetric elliptic fixed points.

Definition 20

An invertible transformation $M:\Omega\longrightarrow\Omega$ is called a reversing symmetry of a local dynamical system $\phi:\mathbb{T}\times\Omega\longrightarrow\Omega$, $\mathbb{T}=\mathbb{R}$ or $\mathbb{T}=\mathbb{Z}$ if the following conditions are satisfied

if $(t,x)\in\mathrm{dom}\,(\phi)$ then $(-t,S(x))\in\mathrm{dom}\,(\phi)$. 2. 2.

$S(\phi(t,x))=\phi(-t,S(x)))$**

Remark 21

In the discrete time case, the above two conditions are equivalent to identity

[TABLE]

where $f=\phi(1,\cdot)$ is a generator of $\phi$.

Definition 22

Let $\phi:\mathbb{T}\times\Omega\to\Omega$ be a local (discrete or continuous) dynamical system. For $x\in\Omega$ put

[TABLE]

The set $\mathcal{O}\left(x\right)$ will be called a trajectory of a point $x$.

Definition 23

Assume $S$ is an reversing symmetry for $\phi:\mathbb{T}\times\Omega\to\Omega$. An orbit $\mathcal{O}\left(x\right)$ is called $S$-symmetric orbit if $\mathcal{O}\left(x\right)=S(x)$.

Remark 24

[La]** In continuous case the orbit $\mathcal{O}\left(x\right)$ is $S$-symmetric if it contains a point from the set $\mathrm{Fix}(S)=\{y:S(y)=y\}$.

Remark 25

[Wi2, Lem.3.3]** It is easy to see that if $\Theta\subset\Omega$ is a Poincaré section for a $R$-reversible flow $\phi:\mathbb{R}\times\Omega\to\Omega$ such that $\Theta=R(\Theta)$ then the Poincaré map $P:\Theta\to\Theta$ is $R|_{\Theta}$-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, elliptic 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

[TABLE]

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

[TABLE]

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 $(\pm c\sqrt{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\approx 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 Poincaré 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 homoclinic solutions.

However, as the bifurcations diagram presented by Michelson suggests [Mi, Fig.1] for all parameter values $c\in(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 Poincaré map will not be a time shift. Therefore computations of the partial derivatives of the Poincaré 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

[TABLE]

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 Poincaré section $\Pi:=\{(0,z,y):z,y\in\mathbb{R}\}$. Let $P_{c}=(P_{1},P_{2}):\Pi{\mathrel{\longrightarrow\mkern-25.0mu\circ}\;\;}\Pi$ be the Poincaré map for the system with the parameter value $c$. Notice, that $P_{c}$ is in fact a half Poincaré map, which means that the trajectory of $x$ crosses $\Pi$ in opposite directions when passing through $x$ and $P_{c}(x)$, and therefore periodic orbits for the Michelson system corresponds to periodic points for $P^{2}$.

Since the section $\Pi$ is invariant under symmetry $(x,y,z)\to(-x,y,-z)$, from Remark 25 the Poincaré map is also reversible with respect to an involution $R(y,z)=(y,-z)$. We will use the same letter $R$ to denote the reversing symmetry of the Poincaré map and the Michelson system.

Let us comment about the choice of the set $C$. In the gap between intervals $C_{1}$ and $C_{2}$ there is a parameter value $c_{*}$ for which the eigenvalues of the Poincaré map $P^{2}_{c_{*}}$ are $\pm 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 $P_{c}$ 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 $c_{i}$ from the grid we proceed as follows

Let $\bar{c}$ denote a center of the interval $c_{i}$. We find an approximate fixed point of $P^{2}_{\bar{c}}$ using the standard nonrigorous Newton method. Let us denote this point by $(y_{i},z_{i})$. 2. 2.

Since the map $P_{c}$ is reversible one can prove the existence of the fixed point for $P_{c}^{2}$ using the shooting method as follows.

Let $\mathrm{Fix}(R)=\{(y,z)\in\Pi:R(y,z)=(y,z)\}=\{(y,0)\in\Pi:y\in\mathbb{R}\}$. Since $P_{c}$ satisfies $(P_{c}\circ R)^{2}=\mathrm{Id}$ whenever the left side is defined, one can see that if $x\in\mathrm{Fix}(R)$ and $P_{c}(x)\in\mathrm{Fix}(R)$ then $P^{2}_{c}(x)=x$. Let us remark, that we always get an approximate fixed points $(y_{i},z_{i})$ resulting from the nonrigorous Newton method very close to $\mathrm{Fix}(R)$. We define two points $u_{1}=(y_{i}-\varepsilon_{i},0),u_{2}=(y_{i}+\varepsilon_{i},0)\in\mathrm{Fix}(R)$, where $\varepsilon_{i}$ is a small number depending on $c_{i}$ and we show that $\pi_{z}(P_{c_{i}}(u_{1}))\cdot\pi_{z}(P_{c_{i}}(u_{2}))<0$, where $\pi_{z}$ is a projection onto $z$ coordinate. Hence, if the $P_{c_{i}}$ is defined on the set $N_{i}=(0,[y_{i}-\varepsilon_{i},y_{i}+\varepsilon_{i}],0)$ then for all parameter values $c\in c_{i}$ there is a point $u_{c}\in N$ which satisfies $\pi_{z}(P_{c}(u_{c}))=0$ and therefore $P_{c}(u_{c})\in\mathrm{Fix}(R)$. This shows that for all $c\in c_{i}$ there exists a fixed point for $P^{2}_{c}$ inside $N_{i}$ provided $P_{c}$ is defined on $N_{i}$, which will be discussed below. 3. 3.

Using $\mathcal{C}^{3}$-Lohner algorithm we compute rigorous bounds for $P^{2}_{c_{i}}(N_{i})$ and $D^{\alpha}P^{2}_{c_{i}}(N_{i})$ for $\alpha\in\mathbb{N}^{2}_{1}\cup\mathbb{N}^{2}_{2}\cup\mathbb{N}^{2}_{3}$. This implies also that $N_{i}\subset\mathrm{dom}\,P_{c_{i}}$. 4. 4.

We show that an arbitrary matrix $M\in DP_{c_{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 $P_{c}$ for $c\in c_{i}$ has the form (44) with $l=1$. 5. 5.

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

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

[TABLE]

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 Poincaré 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 $\mathcal{C}^{3}$ or $\mathcal{C}^{5}$ computations.

Appendix 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 [math] $f$ is given by a real, convergent power series

[TABLE]

Denote also by $f:\mathbb{C}^{2}\to\mathbb{C}^{2}$ a complex extension of $f$. Let $\lambda,\bar{\lambda}\in\mathbb{C}$ be complex eigenvalues of $Df(0)$ and $v,\bar{v}\in\mathbb{C}$ corresponding eigenvectors (here bar denotes the complex conjugation). Then, using a linear substitution of the form $L=[v^{T},\bar{v}^{T}]$, we can change the coordinate system such that in the new coordinates the mapping $f$ has the form

[TABLE]

The last condition is a consequence of the invariance of $\mathbb{R}^{2}\subset\mathbb{C}^{2}$ under the complex map $f$. We will refer to it as the reality condition. Namely, the set $\mathbb{R}^{2}\subset\mathbb{C}^{2}$ in the new coordinates $(\xi,\eta)$ is given by $\xi=\overline{\eta}$ and the condition $f(\mathbb{R}^{2})\subset\mathbb{R}^{2}$ expressed in coordinates $(\xi,\eta)$ is equivalent to (A).

Assume now, that $\lambda^{k}\neq 1$ for $k=1,\ldots,4$. Then an analytic area-preserving substitution satisfying reality condition (A)

[TABLE]

where

[TABLE]

brings $f=(f_{1},f_{2})$ to the normal form

[TABLE]

with

[TABLE]

Finally, let $\gamma_{0}\in\mathbb{R}$ be such that $\lambda=\alpha_{0}=e^{i\gamma_{0}}$ and we compute coefficient $\gamma_{1}$ by

[TABLE]

From the proof given in [SM] it follows that $\gamma_{1}\in\mathbb{R}$ and the mapping $f$ in coordinates $(z,v)$ has the form

[TABLE]

where $O_{4}$ is a convergent power series with the terms of degree at least $4$.

Again, the coefficients of $f(z,v)$ satisfy reality condition (A). In order to express this normal form in terms of real variables we make a linear substitution

[TABLE]

and we obtain the normal form for $f$

[TABLE]

which agrees with (44).

The formulas for higher order terms $\phi_{i,j},\psi_{i,j}$ (and for $\gamma_{2}$, which are not given here) has been computed in Mathematica.