%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                           %
%  History:                                                 %
%  Ne 10.7.2005 na chalupe dokazana domnenka elegantne      %
%  Ut 12.7.2005 14 kB  8 stran                              %
%  Ct 14.7.2005 16 kB 10 stran                              %
%  Pa 15.7.2005 20 kB 12 stran                              %
%  So 16.7.2005 24 kB 15 stran                              %
%     18.-22.7.2005 Bobok                                   %
%     25.-29.7.2005 Dynamics Days                           %
%      1.-8.8.2005 BCH                                      %
%  So  6.8.2005 24 kB 15 stran opravy                       %
%  Po  8.8.2005 32 kB 19 stran                              %
%  Ut  9.8.2005 42 kB 21 stran                              %
%  St 10.8.2005 44 kB 22 stran                              %
%  Ct 11.8.2005 48 kB 25 stran                              %
%  Pa 12.8.2005 51 kB 26 stran                              %
%  Ne 14.8.2005 51 kB 26 stran Filippov                     %
%  Po 15.8.2005 53 kB pridan bch12f.m SparseArray           %
%  Ut 16.8.2005 55 kB 27 stran Djokovic                     %
%  Ct 18.8.2005 63 kB 31 stran Reinsch, jadro graficky      %
%  Pa 19.8.2005 64 kB 32 stran CD-R                         %
%  Ut 23.8.2005 sparsest ?                                  %
%  St 24.8.2005 67 kB 34 stran                              %
%  Ne 28.8.2005 opraveno #c = 5 -> 4 pro r=6                %
%  Po 29.8.2005 68 kB 76 stran pridany programy             %
%  Ct  1.9.2005 potvrzeno min Ncom = 18 pro r=7             %
%  Ne  4.9.2005 potvrzeno min Ncom = 13 pro r=8             %
%      5.-7.9.2005 Herbertov                                %
%  Ct  8.9.2005 opravy                                      %
%  Pa  9.9.2005 76 kB 38 stran sparsest                     %
%  Ut 13.9.2005 77 kB 39 stran Kolsrud                      %
%  Pa 16.9.2005 holky do Italie                             %
%  So 17.9.2005 79 kB Lojzovy opravy                        %
%  Ne 18.9.2005 85 kB 42 stran supremum?                    %
%  Po 19.9.2005 zacina semestr                              %
%  Ct 22.9.2005 86 kB 43 stran                              %
%  So 24.9.2005 90 kB 45 stran normy ala Robert Israel      %
%                              chem.eng.                    %
%  St 28.9.2005 105 kB 54 stran zig+zag+zug                 %
%                               zig+zag-zig-zag             %
%                               acknowledgement             %
%                               own publications, CV        %
%  Pa 30.9.2005 111 kB 58 stran                             %
%  So 1.10.2005 120 kB 64 stran                             %
%  Ne 2.10.2005 122 kB 65 stran                             %
%  Ne 16.10.2005 114 kB 61 stran odstranena konvergence     %
%  Ut 18.10.2005 117 kB 62 stran odstranena poznamka o toku %
%  Pa 21.10.2005 125 kB 64 stran Nelson Kowalski Gunesch    %
%  Ne 23.10.2005 126 kB 65 stran doplnen Nelson             %
%  Po 24.10.2005 acer 1 GB RAM                              %
%  Po  7.11.2005 126 kB 66 stran Lojzuv komentar            %
%  St  9.11.2005 126 kB 66 stran R^n, ven odkaz na Mma code %
%                                Lojza rekl: Dej to svazat  %
%  Pa 11.11.2005 xw9300 8 GB RAM                            %
%  So 26.11.2005 drobne opravy                              %
%  Po 28.11.2005 po castech: 38 pro r= 9 za     9s on acer  %
%  St 30.11.2005             31 pro r=10 za 13070s on acer  %
%  prosinec, leden: thick diag      r=10 za   649s          %
%  Ut 17. 1.2006 do vazby verze z 30.11.2005 8x109,-Kc      %
%  Ut 24. 1.2006 svazano                                    %
%  St 25. 1.2006 odevzdano na dekanat                       %
%  Ct  2. 2.2006 odevzdano oponentum (Bubenik, Jankovsky)   %
%  Po 20. 2.2006 posudek Jankovskeho                        %
%  Ut  7. 3.2006 posudek Bubenika                           %
%  Pa 17. 3.2006 citace v UK Nonlinear News, Issue 42       %
%  Ut 11. 4.2006 Turzik objevil 4 preklepy                  %
%  St 12. 4.2006 posudek Turzika, skolitele, ustavu         %
%  St  3. 5.2006 dlouhy telefon s I.Markem                  %
%  St 10. 5.2006 obhajeno 6:0                               %
%                Kubicek  -predseda                         %
%                Klic     -skolitel                         %
%                Turzik   -openent                          %
%                Bubenik  -openent                          %
%                Jankovsky-openent-zemrel                   %
%                Janovsky                                   %
%                Vonka                                      %
%                Porubsky                                   %
%                Ivo Marek-zdrzel se                        %
%                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[12pt]{article} 
\usepackage{color}
\usepackage{curves}
\usepackage{times}
\usepackage{graphicx}
\usepackage{verbatim}
\newlength{\figsize}
\setlength{\figsize}{0.7\textwidth}
\newcommand{\U}{\hat U}
\newcommand{\V}{\hat V}
\newcommand{\hatI}{\hat{\cal I}}
\newcommand{\ad}{\mbox{Ad }}
\begin{document} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{ 
\LARGE
 {\bf Institute of Chemical Technology, Prague} \\[0.6em]
 Faculty of Chemical Engineering                \\[0.6em]
 Department of Mathematics                      \\[4em]
 { \bf  Zig+Zag Dynamical Systems 
   \\[1.2em]
  (Dynamick\'{e} syst\'{e}my generovan\'{e}\\
  dv\v{e}ma a v\'{\i}ce vektorov\'{y}mi poli)
}             \\[1em]
  PhD Thesis                                    \\[1.4em]
RNDr. Pavel Pokorn\'y
}
\end{center}

\vspace{\stretch{1}}
{\large
\begin{tabular}{ll}
 Supervisor:   & Prof. RNDr. Alois Kl\'{\i}\v{c}, CSc.\\
 Study program:& Applied Mathematics \\
 Study field:  & Applied Mathematics
\end{tabular}
}
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This thesis was written at the Department of mathematics, 
Institute of Chemical Technology, Prague, in 2005,
summarizing my results of the past ten years.
It was defended on May 10, 2006.

\vspace{\stretch{1}}
For the preparation of this thesis I used exclusively those
sources cited in the list of references.
\\[1em]
Prague,
January 25, 2006,\\
corrected May 11, 2006.
\hfill Pavel Pokorn\'{y}
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Acknowledgment:\\
I would like to thank to my ancestors for giving education to their children,
to my wife Jana for love, patience and understanding, 
to my supervisor prof.\,Kl\'{\i}\v{c} for friendly support and scientific 
guidance, to my colleagues for critical comments at the department 
seminar ``Dynamical Systems", to prof.\,Kowalski for stimulating 
discussion and to posters to {\tt sci.math.research} and
{\tt sci.nonlinear} newsgroups for inspiring comments.

\vspace{\stretch{1}}
Typeset by \LaTeX {} under Linux
(the author has used no Microsoft software
for the preparation of this PhD thesis).
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\tableofcontents


\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{table}[hb]
%\begin{center}
\begin{tabular}{|l|l|}
 \hline   
\multicolumn{2}{|c|}{List of symbols}\\
\hline   
% $N$   & the set of all positive integers \\
  $R$   & the set of all real numbers  \\
  $R^n$ & the set of all $n$-tuples of real numbers \\
  $t $  & real time \\
  $x,x_0,x_1,x_2$ & points in $R^n$ \\
  $\psi   (x,t)$ & switching function \\
  $p$ & half of the switching period $T=2p$\\
  $E$ & identity matrix \\
  $G$ & involution matrix (satisfying $G\cdot G=E$) \\
  $u,v,w,U,V,W$ & vector fields (maps from $R^n$ to $R^n$) \\
  $u'$ & matrix of partial derivatives 
      $u'_{j k}=\frac{\partial u_j}{\partial x_k}$\\
  $I$ & identity map (i.e. $I(x) = x$ for all $x \in R^n$)\\
  $\phi$ & smooth map from $R^n$ to $R^n$ \\
  $\Phi$ & the set of all smooth maps from $R^n$ to $R^n$ \\
  $\hat u,\hat v,\hat w,\U,\V,\hat W$ & operators (maps from $\Phi$ to $\Phi$)\\
  $\hatI$&identity operator (i.e. $\hatI(\phi)=\phi$ for all $\phi\in\Phi$)\\
  $\exp$ & exponential $\exp\hat u=\sum_{k=0}^\infty\frac{\hat u^k}{k!}$\\
  $\log$&logarithm $\log\hat u=\sum_{k=1}^\infty(-1)^{k-1}(\hat u-\hatI)^k/k$\\
  $\U\vdash\!\!\!\dashv\V$ & BCH series $= \log(\exp\U\exp\V) $\\
  ${\cal O}(p^3)$& terms of order 3 and higher in $p$ \\
  $[\U,\V]$ & commutator of two operators $=\U\V-\V\U$ \\
  $[\hat A_1,\hat A_2,\dots,\hat A_r]$ & commutator of $r$ operators 
               $=[\hat A_1,[\hat A_2,[\dots,\hat A_r]]]$\\
  $M_r$ & real $2^r\times 2^r$ matrix \\
  $(\hat w_r)_j$&$j$-th word       of order $r$, e.g. $(\hat w_3)_{6}=\V\U\V$\\
  $(     a_r)_j$& real coefficient at $(\hat w_r)_j$\\
  $(\hat c_r)_j$&$j$-th commutator of order $r$, e.g. $(\hat c_3)_{6}=[\V,\U,\V]$\\
  $(     b_r)_j$& real coefficient at $(\hat c_r)_j$\\
  $N_{tot}$&$=2^r$, number of all possible $r$-th order words \\
  $N_w$&number of $r$-th order words in BCH series\\
  $N_{com}$& number of $r$-th order commutators in BCH series\\
  $\varphi(t,x)$& flow (solution to ODE with initial condition $x$)\\
 \hline
\end{tabular}
%\end{center}
% \caption{ List of symbols.  }
% \label{tabSymbols}
%\end{table}

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Summary}

Keywords:\\
dynamical systems, \\
differential operator, \\
Baker-Campbell-Hausdorff formula,\\
sparsest solution.

The original inspiration to this research comes from chemical engineering.
Experimental studies have shown that some catalytic
processes, which were usually operated in stationary regimes,
can be significantly improved (in terms of selectivity, conversion etc.)
if they are forced to operate in periodic regimes.
A special class of these processes includes reactors with periodic
flow reversal.
The detailed mathematical description of such a reactor
would require a system of partial differential equations
with a large number of parameters. 

To understand the basic types of dynamical behavior of such systems,
we started with models in the form of ordinary differential equations
where the flow reversal is modeled by two different vector 
fields 
$u,v: R^n\rightarrow R^n$
acting one after the other
$$
 \frac{d x}{d t}=f(x,t)
$$
where $f(x,t)$ is periodic in $t$ with period $T=2p$ 
satisfying
$$
 f(x,t)
 =\left\{  % } aby byly v paru
       \begin{array}{r c l}
        u(x)  & \mbox{if} \quad  & 0 \le t < \ p \\
        v(x)  & \mbox{if} \quad  & p \le t < 2p.
       \end{array}
       \right.
$$
We call this type of dynamical systems ``zig+zag dynamical systems".


The output of our research activity in this field consists of three main
parts (besides posters and lectures at 
international conferences):
\begin{itemize}
\item
 \mbox{ [KP] }
 Kl\'{\i}\v{c} A., Pokorn\'y P.:
On Dynamical Systems Generated 
by Two Alternating Vector Fields.
Int.\,J.\,Bif.\,Chaos 
{\bf 6}, 2015 (1996).

\item
 \mbox{ [KPR] }
Kl\'{\i}\v{c} A., Pokorn\'y P., \v{R}eh\'{a}\v{c}ek J.:
Zig-Zag Dynamical Systems and the 
Baker-Campbell-Hausdorff Formula.
Math. Slovaca {\bf 52}, 79-97 (2002).

\item
This PhD thesis which we intend to publish in a scientific journal
after minor modifications.

\end{itemize}

In [KP] Dynamical systems generated by two vector fields are studied
both analytically and numerically.
A special case is considered, namely when the two vector fields are
related by an involution. A map $G$ is called involution if its
second iterate is equal to identity. A linear involution is 
represented by a matrix $G$ satisfying 
$$
G\cdot G = E
$$
where $E$ is the identity matrix.
Two vector fields $u,v$ are $G$-related if 
$$
 v(G\cdot x) = G \cdot u(x).
$$
This is a good starting point to investigate systems
modeling chemical reactors with periodic flow reversal,
when the reactant gases flow for a certain time interval
from one side and then for another time interval they flow
from the opposite side.

Four examples of this type of models were chosen in [KP]
for numerical study: 
\begin{itemize}
\item blinking nodes
\item blinking cycles
\item blinking Lorenz and
\item blinking vortices.
\end{itemize}

Detailed numerical investigation suggests that for small switching
time interval~$p$ the resulting system can be approximated
by the averaged system
$$
 \frac{d x}{d t}= w(x)
$$
where
$$
 w(x) = \frac{u(x)+v(x)}{2}.
$$

In [KPR] this observation was refined in the following way.
Taking the average of the two
vector fields is just the first order approximation of an infinite series
approximation.
The system of two blinking cycles was chosen to test our hypothesis.
By blinking cycles we mean the following:
a 2-dim system given in polar coordinates $(r,\phi)$
$$
 \frac{d r}{d t}= r(a^2-r^2)
$$
$$
 \frac{d \phi}{d t}= 1
$$
has a stable limit cycle with radius $a$ and with the center in the origin.
Writing the system in Cartesian coordinates and shifting the center
one unit to the right we get the vector field $u$ while
shifting the center one unit to the left we get the vector field $v$.

Our conjecture in [KPR] was that the vector field
describing the behavior of the resulting system is given by the 
Baker-Campbell-Hausdorff (BCH) formula.
Taking the averaged system and also taking only the first two or three
terms in the BCH series resulted in a dynamical system
with a qualitatively different phase portrait.
It was our great joy to see that taking the first four terms 
gives a system with the phase portrait indistinguishable from the 
original one.

It took us three years to find a satisfactory proof that 
the BCH series is indeed the right form for the approximation
of the zig+zag dynamical systems. 
This PhD thesis brings the desired proof along with further 
results.

The main results presented in this PhD thesis are as follows.
In section \ref{sec.opera} for a given smooth vector field
$u:R^n\rightarrow R^n$ we introduce a linear first order 
differential operator $\hat u$ acting on smooth maps 
$\phi: R^n \rightarrow R^n$
$$
 \hat u ( \phi) = \phi'\cdot u
$$
where both $\phi$ and $u$ are functions of $x\in R^n$.
Thus in full form it reads
$$
 \hat u ( \phi)(x) = \phi'(x)\cdot u(x).
$$
We show that the solution of ODE
$$
 \dot x = u(x)
$$
with the initial condition
$$
 x(0) = x_0
$$
can be formally written (assuming the series converges) as
$$
 x(t) = 
 \sum_{k=0}^\infty\frac{t^k}{k!}x^{(k)}(0)=
 \sum_{k=0}^\infty\frac{t^k}{k!}\hat u^k I x_0=
 \exp (t\hat u) I x_0.
$$

In section 
\ref{sec.twovf}
we apply this formalism to a dynamical system generated by two
alternating vector fields $u$ and $v$ acting one after the other.
Our aim is to find a shortcut -- a third vector field $w$ that 
moves the point in the state space from the same initial
condition to the same final point in the same time. 
In section 
\ref{sec.shortTaylor}
we use the Taylor expansion to find 
the first few terms of the expansion of the shortcut~$w$.
In section 
 \ref{sec.shortBCH}
we find the full form of the expansion of the shortcut $w$.
To this purpose we first have to replace composition of flows
by composition of operators
$$ 
 \exp(\V)I \circ \exp(\U)I = \exp(\U)\exp(\V) I.
$$
Then we can use the good old Baker--Campbell--Hausdorff formula
that gives 
(for non-commuting $x$ and $y$)
the solution to
$$
 \exp z = \exp x \exp y
$$
for $z$ in the form of an infinite series of commutators of $x$ and $y$.
For commuting $x$ and $y$ (which is the case for real or complex numbers)
we have just
$$
  z = x + y
$$
while for non-commuting $x$ and $y$ 
(which is the case for square matrices or operators)
we have
$$
 z = x + y + \frac{1}{2}[x,y] + \frac{1}{12}([x,x,y]-[y,x,y])
 - \frac{1}{24}[x,y,x,y] +
$$
$$
 + \frac{1}{6!}
(-[x,x,x,x,y]-6[x,x,y,x,y]-2[x,y,y,x,y]+
$$
$$
+2 [y,x,x,x,y]+6[y,x,y,x,y]+[y,y,y,x,y])
 + 
$$
$$
+ \frac{1}{2 \cdot 6!}
([x,y,y,y,x,y]+ [y,x,x,x,x,y]+
$$
$$
+6[y,x,x,y,x,y] -2[y,y,x,x,x,y]) + \dots
$$
Here square brackets denote nested commutator
defined by (\ref{commutator}) and (\ref{ncommutator}).

This form is not unique and in section
  \ref{sec.BCH}
we derive detailed relation between the form of the BCH series
in words and in commutators.
This relation can be expressed conveniently by an infinite series
of square matrices $M_r$, where each $M_r$ is a $r \times r$ matrix
of integer numbers.
Then we apply our results to derive identities between 
commutators. The simplest of them being
$$
 [x,y] + [y,x] = 0.
$$
Each such identity corresponds to one vector of the kernel of the 
matrix $M_r$.
These identities can be used to reduce the number of terms 
of a given order in the BCH series.
We formulate the problem to find the sparsest form of the BCH series
(i.e. the form that has the least number of terms of a given order).
This problem can be solved by finding the sparsest solution
to a given linear under-determined system of algebraic equations.
We give an explicit algorithm to find the sparsest solution
and we present numerical examples for illustration.

In section 
 \ref{sec.zzz}
we generalize our previous results to more than two vector fields.

In section 
  \ref{sec.zzzzs}
we apply our results to the case of a dynamical system 
generated by four vector fields $u,v,-u,-v$
and give a useful interpretation of the commutator $[u,v]$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Introduction}
% \chapter { Introduction }
 \label{sec.intro}

Considerable attention has been devoted to dynamical systems
in the form of a set of~$n$ first order ordinary 
differential equations (ODE)
\begin{equation}
 \frac{d x}{d t}=f(x)
\end{equation}
where $f: R^n \rightarrow R^n$ is a continuous vector field.

The present work deals with a generalization 
of the above dynamical system, namely with dynamical 
systems with discontinuous right hand side
\begin{equation}
 \frac{d x}{d t}=f(x,t)
\end{equation}
where
\begin{equation}
f(x,t)=\left\{  % } aby byly v paru
       \begin{array}{ll}
        u(x) & \mbox{if} \quad  \psi   (x,t) > 0\\
        v(x) & \mbox{if} \quad  \psi   (x,t) < 0
       \end{array}
       \right.
\end{equation}
where $\psi   (x,t)$ (called the switching function)
is a continuous function of both the state
$x$ and of the time~$t$.

We do not deal with the question how to define 
$f(x,t)$ for 
$\psi   (x,t) = 0$ 
here.
This question is important when 
$\psi   (x,t)$ depends on~$x$, see \cite{Filippov}.
Here we can take either one of the single side limits.

We consider a special case when the switching function 
$\psi   (x,t)$
does not depend on~$x$ 
and is periodic in~$t$ with period $T=2p$
and has the special form
\begin{equation}
\psi   (x,t)
\left\{  % } aby byly v paru
       \begin{array}{r c c}
        >0 & \mbox{if} \quad  & 0 < t < \ p \\
        <0 & \mbox{if} \quad  & p < t < 2p.
       \end{array}
       \right.
\end{equation}
So we arrive at the {\bf zig+zag dynamical system}
\begin{equation}
 \frac{d x}{d t}=f(x,t)
\end{equation}
where $f(x,t)$ is $2p$-periodic in~$t$ satisfying
\begin{equation}
 f(x,t)
 =\left\{  % } aby byly v paru
       \begin{array}{r c l}
        u(x)  & \mbox{if} \quad  & 0 \le t < \ p \\
        v(x)  & \mbox{if} \quad  & p \le t < 2p.
       \end{array}
       \right.
\end{equation}

A note on the term ``zig+zag":
an alternative name is ``zig-zag" where the hyphen (``-")
means simple joining of two words.
We use ``zig+zag" intentionally 
to distinguish between plus (``+") and minus (``-") where minus 
stands for the opposite of plus. 
This convention is useful
when studying systems ``zig+zag-zig-zag", i.e.
systems generated by four vector fields $u,v,-u,-v$
where the third vector field is equal to the first one when 
multiplied by minus one and the fourth vector field
is equal to the second one when multiplied by minus one,
see chapter~\ref{sec.zzzzs}.
Still another name could be ``dynamical systems generated
by two or more alternating vector fields".

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Motivation}

The motivation to investigate zig+zag dynamical systems
comes from chemical engineering, 
see \cite{RussoC}, \cite{RussoI}, \cite{Rehacek1}, \cite{Rehacek2}.
There has been an increasing interest in periodically forced processes
in literature recently. Application of forced unsteady
operating conditions seems to be beneficial in several
heterogeneous catalyzed processes 
(\cite{Bunim}, \cite{Snyder} \cite{Xiao}, \cite{Khinast}).
Experimental studies 
(\cite{Zufle}, \cite{Nieken}) 
have shown that some catalytic
processes, which were usually operated in stationary regimes,
can be significantly improved (in terms of selectivity, conversion etc.)
if they are forced to operate in periodic regimes.
Theoretical and experimental studies have shown that several
processes can be operated auto-thermally avoiding the use of heat
exchanger \cite{Matros85}, \cite{Matros89}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Previous results}

In \cite{KP96} we studied zig+zag dynamical systems,
where the two vector fields are related by a linear involution~$G$ 
(a map is called involution if its second iteration is identity)
i.e. we studied systems of the form
$$
 v(G\cdot x) = G \cdot u(x)
$$
where the matrix~$G$ satisfies 
$$
 G \cdot G = E
$$
here~$E$ is the identity matrix.

%Four examples 
%with different types of dynamics
%of the original system
%$$
% \frac{d x}{d t} = u(x)
%$$
%were presented, namely:
%\begin{itemize}
%\item
% {\bf blinking nodes} 
% when the vector field~$u$ has a single stable stationary point,
%\item
% {\bf blinking cycles}
% when the original system
% has a stable limit cycle,
%\item
% {\bf blinking Lorenz}
% when~$u$ is the vector field 
% from the Lorenz model \cite{Lorenz},
%\item
% {\bf blinking vortices}
% when the original system is a vortex \cite{Ottino89}.
%\end{itemize}

For small switching period~$p$
the method of averaging was used 
to approximate the non-autonomous
zig+zag system by an autonomous system
$$
 \frac{d x}{d t} = \frac{u(x)+v(x)}{2}.
$$

In \cite{Maslo2002} we generalized the method
of averaging by including higher order terms leading to the 
Baker - Campbell - Hausdorff (BCH) 
series for linear systems
and we used the BCH series heuristically for nonlinear system.

The purpose of the present work is to 
give rigorous reasoning for the usage 
of the BCH series for nonlinear systems.
First, we develop the operator formalism to deal with solution
to ODE.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Operator formalism}
 \label{sec.opera}

\subsection {Introduction}

Although we will work only with vector fields on $R^n$ 
that can be given by ordered $n$-tuples of smooth functions
of $n$ variables, it turns out to be useful and necessary 
to consider vector fields as differential operators,
in agreement with \cite{Olver}.
We use hat to stress the operator nature of a vector field, i.e. 
$\hat u$ will denote the vector field on $R^n$ viewed as 
a differential operator and $u$ will be 
its coordinate form in local coordinates.
As we work on $R^n$ only, the coordinate map on $R^n$ 
is the identity map ($Ix=x\ \forall x \in R^n$) and 
$\hat u I = u$ is the coordinate form of the vector field $\hat u$
in the coordinate system on $R^n$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Vector field as an operator}

% For a fixed positive integer~$n$ 
Consider ODE
\begin{equation}
 \label{ODE}
 \frac{d x}{d t} = u(x)
\end{equation}
with a smooth vector field $u:R^n \rightarrow R^n$
(where $n$ is a fixed positive integer)
and with the initial condition 
$$
 x(0) = x_0.
$$
By {\bf smooth} we mean in $C^\infty(R^n)$
i.e. having continuous derivatives of any order.
We shall write~$\dot{x}$ instead of $\frac{d x}{d t}$ where 
possible.
Assuming the solution $x(t) = \varphi(t,x_0)$ can be expressed 
as a Taylor series
\begin{equation}
 \label{Tay}
 x(t) = x(0) + t \dot{x}(0) + \frac{t^2}{2} \ddot{x}(0) + \dots
\end{equation}
for~$t$ in some neighborhood of zero,
we need the time derivatives of~$x$

\begin{eqnarray}
\label{xdots}
 \nonumber
 \dot{x}(t) &=& u(x(t))
 \\
 \nonumber
 \ddot{x}(t)&=& u'(x(t))\cdot u(x(t))
 \\
 \nonumber
 & \vdots & 
\end{eqnarray}
%in short
%\begin{eqnarray}
%\label{xdotsshort}
% \nonumber
% \dot{x}&=& u
% \\
% \nonumber
% \ddot{x}&=& u'\cdot u.
% \\
% \nonumber
% & \vdots & 
%\end{eqnarray}
Here~$u'$ is the matrix of partial derivatives
$$
 u'_{j k}=\frac{\partial u_j}{\partial x_k}
$$
and the dot stands for multiplication of a matrix and a vector
$$
 (u'\cdot u)_j = \sum_{k=1}^{n} u'_{j k} u_k.
$$
We will use multiple primes to denote higher order derivatives
e.g. 
$$
 u''_{j k l}=\frac{\partial^2 u_j}{\partial x_k \partial x_l}.
$$
Using the chain rule for the derivative of the composition
of maps\\$s:R^n\rightarrow R^n$ and $x:R\rightarrow R^n$
$$
 \frac{d}{d t}s(x(t)) = s'(x(t)) \cdot \dot{x}(t)
$$
we can find the time derivative of~$x$ of any order:
\begin{eqnarray}
\label{uuu}
\nonumber
 \dot{x} & = & u \\
\nonumber
 \ddot{x} & = &   u'\cdot u \\
\nonumber
 x^{(3)}  & = &  (u'\cdot u)'\cdot u \\
\nonumber
 x^{(4)}  & = & ((u'\cdot u)'\cdot u)' \cdot u \\
 & \vdots & 
\end{eqnarray}
(each left hand side is evaluated in~$t$ and 
each right hand side is evaluated in~$x(t)$).

Could we introduce multiplication of vector fields
$u,v:R^n\rightarrow R^n$ by 
\begin{equation}
 \label{badmul}
 u v = u'\cdot v
\end{equation}
then the Taylor series~(\ref{Tay}) would become an exponential.
Unfortunately, this multiplication 
would not be associative because for
$u,v,w:R^n\rightarrow R^n$ 
$$
(u v)w = (u'\cdot v)'\cdot w
$$
but on the other hand
$$
u(v w) = u'\cdot(v'\cdot w)
$$
which is different. 
We still want to express the Taylor series (\ref{Tay})
as an exponential. 

To this purpose for any given smooth vector field
$u:R^n\rightarrow R^n$ 
we introduce 
a linear first order differential 
{\bf operator}
$$
 \hat u: \Phi \rightarrow \Phi.
$$
Here $\Phi$ is the set of all smooth maps
$\phi : R^n\rightarrow R^n$.
The action of the operator~$\hat u$ on any 
map $\phi \in \Phi$
is given by 
\begin{equation}
 \label{poleop}
 \hat u \phi = \phi' \cdot u
\end{equation}
i.e.  $\hat u \phi$ 
is a map that assigns to each $x \in R^n$ the point 
$$
 x\mapsto
\hat u\phi x= 
(\hat u\phi)(x)= 
      (\phi'\cdot u)(x)=
      \phi'(x) \cdot u(x).
$$
The alternative grouping, namely 
$\hat u(\phi(x))$ does not make sense, since
the operator can be applied to a map, here~$\phi$,
not to a point~$\phi(x)$.

The action of the operator~$\hat u$ on the identity map~$I$,
i.e. a map satisfying
\begin{equation}
 \label{I}
 I(x) = x \qquad \mbox{ for all } x \in R^n
\end{equation}
is
\begin{equation}
 \label{oppole}
 \hat u I = I' \cdot u = E \cdot u = u
\end{equation}
where~$E$ is the $n\times n$ identity matrix.

% tohle bych tam rad nechal
% Thus, (\ref{poleop}) 
% defines the operator~$\hat u$ when the vector field~$u$ is given,
% while (\ref{oppole}) 
% gives the vector field~$u$ when the operator~$\hat u$ is given.

We also define the identity operator~$\hatI$
by
\begin{equation}
 \label{Ihat}
  \hatI \phi = \phi \qquad \mbox{ for all } \phi\in\Phi. 
\end{equation}
The identity operator $\hatI$ is of order zero
(no derivative of $\phi$ is necessary to evaluate $\hatI \phi = \phi$)
and $\hatI$ does not correspond to any vector field,
because each operator corresponding to a vector field is of order one.

Note: we will apply the operators only to the identity map~$I$ 
or to a map that results from applying another operator 
to the identity map~$I$, which is actually applying the 
composition of operators to the identity map.

Powers of~$\hat u$ can be used to express 
time derivatives of~$x(t)$ in (\ref{uuu}).
Namely
$$
 x^{(k)}(t) = \hat u^k  I x(t)
$$
meaning 
\begin{eqnarray*}
%\label{uhatuhatuhat}
 \nonumber 
 \hat u I &=&u \\ \nonumber 
 \hat u^2I&=&\hat u \hat u I = \hat u u = u'\cdot u \\ \nonumber
 \hat u^3I&=&\hat u u'\cdot u = (u'\cdot u)'\cdot u \\ \nonumber
 \hat u^4I&=&\hat u (u'\cdot u)'\cdot u = ((u'\cdot u)'\cdot u)' \cdot u\\
 & \vdots & 
\end{eqnarray*}
And for $t=0$ we get
$$
 x^{(k)}(0) = \hat u^k  I x_0.
$$
We define (in accord with \cite{Howe}, \cite{Olver} etc.)
the {\bf exponential} of an operator~$\hat u$ as
\begin{equation}
 \label{expop}
 \exp \hat u = \sum_{k=0}^\infty \frac{\hat u ^ k}{k!}.
\end{equation}
%(assuming the series converges in accord with \cite{Olver}, \cite{Howe} etc.).
Note that if~$\hat u$ is a first order differential operator 
(corresponding to a vector field $u:R^n\rightarrow R^n$),
then 
$\hat u^2$ is a second order differential operator,
$\hat u^3$ is a third  order differential operator etc.
and 
$\exp \hat u$
 has not a finite order.
% Neither 
% $\hat u^2$, nor
% $\exp \hat u$
% correspond to a vector field.

Now we can express the Taylor series (\ref{Tay}) 
by the exponential
of the operator ~$\hat u$ corresponding to the vector field~$u$
on the right hand side of the ODE (\ref{ODE})
\begin{equation}
 \label{expopsol}
 x(t) = 
 \varphi(t,x_0)=
 \sum_{k=0}^\infty\frac{t^k}{k!}x^{(k)}(0)=
 \sum_{k=0}^\infty\frac{t^k}{k!}\hat u^k I x_0=
 \exp (t\hat u) I x_0.
\end{equation}

We stress that the exponential in 
(\ref{expopsol}) 
is a formal series, which is equivalent to the Taylor
series of the solution of the given ODE.
The exponential
$$
 \exp (t\hat u) 
$$
is not an operator in the sense that 
when applied to a function it gives another function 
defined on a certain 
subset of $R^n$ independent of $t$.
As we will see in the next example,
(\ref{expopsol}) 
gives the solution of the ODE which is defined on a certain interval
of $t$ that may depend on the initial condition,
which is a common case for ODE.

Note:
some authors (see e.g. \cite{Olver}) 
introduce differential operators acting on real scalar--valued functions alone
and only then they let the operator act com\-po\-nent--wise 
on vector--valued functions. We believe that our approach,
namely to introduce differential operators acting 
on vector--valued functions directly is more natural.
Also, by using the hat notation 
(adopted from quantum mechanics)
we stress the difference between a map (acting on points in $R^n$)
and an operator (acting on maps from $R^n$ to $R^n$).
This difference is often neglected.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Example 1: ODE}

To illustrate the operator formalism just introduced 
consider a differential equation 
$$
 \dot x = x^2
$$
with the initial condition $x(0)=x_0 \in R$.

First, let us solve this ODE by separation of variables:
$$
 x(t) = \frac{x_0}{1-t\ x_0}.
$$
For 
$$
 |t\ x_0| < 1
$$
this solution can be expressed as a Taylor series
$$
 x(t) = x_0 \sum_{k=0}^\infty (t\ x_0)^k.
$$

As the second approach we will use the operator formalism 
introduced above.
For the 1-dim vector field $u(x) = x^2$ we have the corresponding 
operator 
$$
 \hat u (\phi) = \phi' \cdot u
$$
i.e.
$$
 \hat u (\phi)(x) = \phi'(x) \cdot x^2.
$$
This operator when applied to the identity map $I(x)=x$ is
$$
 (\hat u I) (x)  = x^2
$$
and its $k$-th power is
$$
 (\hat u^k I) (x)  = k! \ x^{k+1}.
$$
The exponential of the operator $t\hat u$ is
$$
 \exp(t\hat u) = \sum_{k=0}^\infty \frac{(t\hat u)^k}{k!} 
$$
and this exponential when applied to the identity map $I$ gives in $x$
$$
 (\exp(t\hat u)I) (x) = 
 (\sum_{k=0}^\infty \frac{(t\hat u)^k}{k!}I)(x)=
 \sum_{k=0}^\infty t^k x^{k+1} =
 x \sum_{k=0}^\infty (t x)^k.
$$
Then the solution can be expressed as
$$
 x(t) = (\exp(t\hat u)I) (x_0) = x_0 \sum_{k=0}^\infty (t x_0)^k
$$
which is exactly the same result as we obtained by standard methods.

% It is important to note that the operator formalism introduces
% no further convergence conditions.
The result obtained by operators is equivalent to the result
obtained by Taylor series in both its values and its domain
of convergence.
This example also demonstrates that we do not have to suppose
the existence of a global flow, i.e. a solution defined
for all real $t$. In this example 
for positive initial condition the solution exists for $t<1/x_0$ and 
for negative initial condition the solution exists for $t>1/x_0$ and
for zero     initial condition the solution exists for $t\in R$.

% In this example the domain of convergence of the Taylor series is the
% same as the domain of definition of the solution. 
% In general
The domain of convergence of the Taylor series 
is a subinterval of the domain of definition of the solution.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Two vector fields}
 \label{sec.twovf}

Consider a zig+zag dynamical system generated by two vector 
fields~$u,v$
\begin{equation}
 \label{zzode}
 \dot x=f(x,t)
\end{equation}
where $f(x,t)$ is $2p$-periodic in~$t$ satisfying
\begin{equation}
 \label{zzrhs}
 f(x,t)
 =\left\{  % } aby byly v paru
       \begin{array}{r c l}
        u(x)  & \mbox{if} \quad  & 0 \le t < \ p \\
        v(x)  & \mbox{if} \quad  & p \le t < 2p
       \end{array}
       \right.
\end{equation}
with the initial condition $x(0)=x_0$, see Fig.\,\ref{obrzigzag}.

\begin{figure}[htb]
\begin{center}
\setlength{\unitlength}{0.5pt}
\begin{picture}(500,400)(50,230)
 \curve(378, 566,  230, 510,   87, 417)
 \curve( 87, 417,  227, 438,  384, 422)
 \curve(384, 422,  233, 395,  100, 328)
 \curve(100, 328,  227, 362,  363, 360)
 \curve(363, 360,  217, 334,   92, 276)
 \curve( 92, 276,  206, 314,  335, 322)
 \curve(335, 322,  194, 294,   74, 236)
 \curve( 74, 236,  186, 278,  296, 291)
 \put(161,530){\framebox(110,35){$\dot{x} = u(x)$}}
 \put(241,440){\framebox(100,30){$\dot{x} = v(x)$}}
 \put(380,562){$x_0$}
 \put(  5,400){$x_1=\exp(p\hat u)Ix_0$}
 \put(385,417){$x_2=\exp(p\hat v)Ix_1$}
% \put(370,352){$P^2(x_0)$}
% \put(346,315){$P^3(x_0)$}
% \put(310,280){$P^4(x_0)$}
 \put(170,476){\line(1,1){20}}
 \put(170,476){\line(3,1){26}}
 \put(190,436){\line(-5,1){24}}
 \put(190,436){\line(-3,-1){22}}
 \put(189,378){\line(4,3){20}}
 \put(189,378){\line(5,1){26}}
 \put(200,358){\line(-5,1){24}}
 \put(200,358){\line(-2,-1){22}}
 \put(195,327){\line(4,3){20}}
 \put(195,327){\line(5,0){24}}
 \put(203,313){\line(-5,1){22}}
 \put(203,313){\line(-2,-1){20}}
 \put(195,294){\line(4,3){20}}
 \put(195,294){\line(5,0){24}}
 \put(209,283){\line(-5,1){22}}
 \put(209,283){\line(-2,-1){20}}
\end{picture}
\end{center}
 \caption{
  A typical trajectory of a zig+zag dynamical system 
  (\ref{zzode}--\ref{zzrhs})
  with the initial condition $x(0)=x_0$. In time~$p$ the 
  vector field~$u$ moves 
  the state point from~$x_0$ to~$x_1$ and then in another time~$p$ 
  the vector field~$v$ moves  the state point from~$x_1$ to~$x_2$.
  And then these two steps repeat periodically 
  producing a zig+zag shaped trajectory giving the name to this
  class of dynamical systems.
}
 \label{obrzigzag}
\end{figure}

In the first half of the period, i.e. in time~$p$, the 
vector field~$u$ moves the state point from the initial 
condition $x(0)=x_0$ to the point 
$$
 x_1=\exp(p\hat u) I x_0.
$$
And then in the second half of the period, i.e. in another 
time~$p$, the vector field~$v$ moves the state point from
the point~$x_1$ to the point 
\begin{equation}
 \label{x2}
 x_2=\exp(p\hat v) I (x_1)
    =\exp(p\hat v) I (\exp(p\hat u) I (x_0)).
\end{equation}

Our aim is to find a {\bf shortcut}
(a short way from~$x_0$ to~$x_2$) 
-- a single vector field, let us
call it~$w$, that moves the state point from the same 
initial condition~$x_0$ to the same final point~$x_2$
in the same time~$2p$ (not necessarily going via
the point~$x_1$, actually almost surely not). 
Thus replacing the non-autonomous
dynamical system~(\ref{zzode}) by the autonomous system
$$
 \dot x = w(x).
$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Shortcut by Taylor}
 \label{sec.shortTaylor}

Assuming the solution to the ODE
$$
\dot x = u(x)
$$
can be written as 
$$
 x(t)=x(0)+t u(x(0))+\frac{t^2}{2}u'(x(0))\cdot u(x(0))+{\cal O}(t^3)
$$
we have (omitting the argument if equal to~$x(0)$)
$$
 x_1 = x_0+p u + \frac{p^2}{2} u'\cdot u+{\cal O}(p^3)
$$
and similarly for the ODE
$$
 \dot x = v(x)
$$
assuming we can write the solution in the form of a Taylor series
we have
$$
 x_2=x_1+p v(x_1)+\frac{p^2}{2}v'(x_1)\cdot v(x_1)+{\cal O}(p^3).
$$
Now, we want to express $v(x_1)$ and $v'(x_1)$ 
in terms of $v(x_0)$ and $v'(x_0)$. Using
$$
 x_1 - x_0 = p u + {\cal O}(p^2)
$$
we get
$$
 v(x_1) = v(x_0) + v'(x_0) \cdot p u(x_0) + {\cal O}(p^2)
$$
and
$$
 v'(x_1) = v'(x_0) + {\cal O}(p)
$$
thus (omitting the argument if equal to~$x_0$)
$$
 x_2 = x_1 + p (v + v'\cdot p u) + \frac{p^2}{2}v'\cdot v 
 + {\cal O}(p^3)
 =
$$
$$
 = 
 x_0+p u + \frac{p^2}{2} u'\cdot u
  + p (v + v'\cdot p u) + \frac{p^2}{2}v'\cdot v + {\cal O}(p^3)=
$$
\begin{equation}
 \label{x2uv}
 = 
 x_0+p(u+v)+ \frac{p^2}{2}(u'\cdot u
  +2v'\cdot u+v'\cdot v)+ {\cal O}(p^3).
\end{equation}
We want (\ref{x2uv}) to be equal to the solution of the ODE 
$$
 \dot x = w(x) 
$$
for $t=2p$, namely
$$
 x_2 = x_0 + 2p w + \frac{1}{2}(2p)^2w'\cdot w + {\cal O}(p^3).
$$
We can search for~$w$ in the form of a series in~$p$.
Denoting $w_0$ the leading term and denoting $p w_1$ the 
first order term and neglecting higher order terms we can write
$$
 w = w_0 + p w_1 + {\cal O}(p^2).
$$
Then
$$
x_2=
x_0+2p(w_0+p w_1)+\frac{1}{2}(2p)^2(w_0+p w_1)'\cdot(w_0+p w_1)
+{\cal O}(p^3)=
$$
\begin{equation}
 \label{x2w}
=
x_0+p(2w_0)+
p^2(2w_1+2w_0\cdot w_0)+{\cal O}(p^3).
\end{equation}
Comparing the terms 
of $x_2$ from (\ref{x2uv}) and (\ref{x2w})
linear in~$p$ we have
$$
 w_0 = \frac{u+v}{2}.
$$
In the method called averaging $w$ is approximated by $w_0$ only.

Using this~$w_0$ and 
comparing the terms 
of~$x_2$ from (\ref{x2uv}) and (\ref{x2w})
of order~$p^2$ 
we get
$$
 w_1 = \frac{
 v'\cdot u - u'\cdot v
 }{4}
 .
$$
Together we get
\begin{equation}
 \label{w2Taylor}
 w = w_0 + p w_1 
 + {\cal O}(p^2)
 = 
 \frac{u+v}{2}
 +
 p\frac{
 v'\cdot u - u'\cdot v
 }{4}
 + {\cal O}(p^2).
\end{equation}

This approach can be used to any order. 
However, a more elegant way is to use the BCH  series,
which gives more insight.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Shortcut by BCH}
 \label{sec.shortBCH}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Composition of flows}

First, 
we want to express the composition of maps
$$
 \exp(p\hat v)I \circ \exp(p\hat u)I
$$
from (\ref{x2}) by appropriate 
composition of operators (applied to 
the identity map).

For this purpose consider a map $s:R^n\rightarrow R^n$ 
applied to the solution of the ODE 
$$
 \dot x = u(x)
$$
i.e. applied to the map
$$
 x(t) = \exp(t\hat u)Ix_0.
$$
Let us call the composition~$S$
$$
 S = s \circ x
$$
i.e.
$$
 S(t) = s(x(t)).
$$
Using the chain rule we can find the derivatives of~$S$:
\begin{equation}
 \label{Sdots}
 \begin{array}{c c c c c c}
&  S(t) &=& s(x(t)) &=& \hat s I x(t) \\ 
&  \dot S(t) &=& s'(x(t))\cdot u(x(t)) &=& \hat u \hat s I x(t) \\
\mbox{in short }&  \dot S &=& s'\cdot u&=& \hat u \hat s I x(t) \\
                & \ddot S &=& (s'\cdot u)'\cdot u&=& \hat u^2 \hat s I x(t) \\
& S^{(3)}&=& ((s'\cdot u)'\cdot u)' \cdot u&=& \hat u^3 \hat s I x(t) 
 \end{array}
\end{equation}
in general
$$
 S^{(k)}=\hat u^k \hat s I x(t).
$$
Assuming we can express~$S(t)$ in the form of a Taylor series we have
$$
 s(x(t)) = 
 S(t)=
\sum_{k=0}^\infty \frac{t^k}{k!}S^{(k)}(0)=
\sum_{k=0}^\infty \frac{t^k}{k!} \hat u^k \hat s I x_0=
\exp(t\hat u)\hat s I x_0.
$$
Then consider 
\begin{equation}
 \label{sxt}
 s(x(t)) = \exp(t\hat u) \hat s I x_0
\end{equation}
for a special operator~$\hat s$
\begin{equation}
 \label{sop}
 \hat s = \exp(p\hat v) 
\end{equation}
and its corresponding map (the flow of the ODE $\dot x =v(x)$)
\begin{equation}
 \label{smap}
 s = \hat s I = \exp(p\hat v) I.
\end{equation}
Then after putting (\ref{sop}) and (\ref{smap}) into (\ref{sxt}) we have
$$
 \exp(p\hat v) I (x(t)) = \exp(t \hat u) \exp(p\hat v)Ix_0
$$
and
\begin{equation}
 \label{conjectureg}
 \exp(p\hat v)I \circ \exp(t\hat u)I = \exp(t\hat u) \exp(p\hat v) I.
\end{equation}
This is a key result that allows to treat zig+zag dynamical
systems by operators because we can study composition of flows
by investigating ``multiplication" of exponentials of operators.

For $t=p$ we have
\begin{equation}
 \label{conjecture}
 \exp(p\hat v)I \circ \exp(p\hat u)I = \exp(p\hat u) \exp(p\hat v) I.
\end{equation}
For three vector fields 
$u,v,w:R^n\rightarrow R^n$ 
and for $p>0$ 
we can introduce three new vector fields
$U,V,W:R^n\rightarrow R^n$ by
\begin{eqnarray}
 \label{uU}
 U &=& p u \\  
 \label{vV}
 V &=& p v \\  
 \label{wW}
 W &=&2p w.
\end{eqnarray}
Then the corresponding operators satisfy
\begin{eqnarray}
 \label{uUop}
 \U &=& p \hat u \\  
 \label{vVop}
 \V &=& p \hat v \\  
 \label{wWop}
 \hat W &=&2p \hat w.
\end{eqnarray}
Then (\ref{conjecture}) can be written in a simpler form
\begin{equation}
 \label{sconjecture}
\fbox{$ \exp(\V)I \circ \exp(\U)I = \exp(\U)\exp(\V) I.$}
\end{equation}
%Personal note: It took me more than 6 years to prove this 
%relation in a clear way.
To avoid confusion, if~$\U$ and~$\V$ are operators, so are their
exponentials $\exp(\U)$ and $\exp(\V)$. Then $\exp(\U)I$ and 
$\exp(\V)I$ are these two operators applied to the identity map~$I$
from~$R^n$ to~$R^n$ 
resulting in another two maps
from~$R^n$ to~$R^n$.
On the left we have the composition of these two maps.
On the right the composition of the two 
exponentials 
$\exp(\U)$ and $\exp(\V)$
of operators being another operator 
is applied to the identity map~$I$
from~$R^n$ to~$R^n$
resulting in another map
from~$R^n$ to~$R^n$.

Note the {\bf reversed order} of~$\U$ and~$\V$ on the left and
on the right hand side.
Also note that (\ref{sconjecture}) is far from obvious. E.g.
it is not true when the identity map~$I$ is replaced
by another map~$\phi$, because then even the zero order terms
do not agree: on the left we have $\phi \circ \phi$ while
on the right we have just $\phi$ which is not equal unless
$\phi$ is a projection.

Due to the rescaling (\ref{uUop}-\ref{wWop}),
$\U$ and~$\V$ are proportional to~$p$
(for fixed~$\hat u$ and~$\hat v$).
When~$p$ is small, so are~$\U$ and~$\V$ and
we can work with series expansions.

As an example we give the second order expansion of both 
sides of (\ref{sconjecture}). On the left we have
$$
 \exp(\V)I \circ \exp(\U)I = 
 (\hatI + \V + \frac{1}{2}\V\V)I \circ
 (\hatI + \U + \frac{1}{2}\U\U)I 
 + {\cal O}(p^3)
 =
$$
$$
 = (I + V + \frac{1}{2}V'\cdot V) \circ
 (I + U + \frac{1}{2}U'\cdot U)
 + {\cal O}(p^3)
=
$$
$$
=
 I + U + V 
+ \frac{1}{2}U'\cdot U
+ \frac{1}{2}V'\cdot V
+ V'\cdot U
 + {\cal O}(p^3)
$$
and on the right we have
$$
 \exp(\U)\exp(\V) I
 =
 (\hatI + \U + \frac{1}{2}\U\U)
 (\hatI + \V + \frac{1}{2}\V\V)I
 + {\cal O}(p^3)=
$$
$$
 =
 (\hatI + \U + \V + \frac{1}{2}\U\U + \frac{1}{2}\V\V + \U\V)I
 + {\cal O}(p^3)=
$$
$$
=
 I + U + V 
+ \frac{1}{2}U'\cdot U
+ \frac{1}{2}V'\cdot V
+ V'\cdot U
 + {\cal O}(p^3).
$$

The presence of the term $V'\cdot U$ in the second order expansion
shows that it is not possible to exchange $\U$ and $\V$.
%agrees with the remarkable fact that the order of~$\U$ and~$\V$
%is reversed on the left vs. the right side.

Let us return to our original goal.
The first vector field~$u$ moves the point in the state space
from the initial condition~$x_0$ to the point
$$
 x_1 = \exp(p \hat u)I x_0
$$
in time~$p$
and then the second vector field~$v$ moves the point from~$x_1$ to
$$
 x_2 = \exp(p \hat v)I x_1 = 
 \exp(p \hat v)I  \circ
 \exp(p \hat u)I x_0
$$
in another time~$p$.
We want to find the shortcut -- 
the vector field~$w$ that brings the point in the state space 
from the initial condition~$x_0$ to the final point~$x_2$
in the same time i.e.~$2p$
$$
 x_2 = \exp(2p\hat w)Ix_0.
$$
So we have the equation for~$\hat w$
$$
 \exp(p \hat v)I  \circ
 \exp(p \hat u)I = 
 \exp(2p\hat w)I.
$$
Due to the rescaling (\ref{uUop}-\ref{wWop}) we can write
$$
 \exp(\V)I  \circ
 \exp(\U)I = 
 \exp(\hat W)I.
$$
Using (\ref{sconjecture}) 
to the left hand side
we get 
$$
 \exp \U \exp \V I = \exp \hat W I.
$$
Assuming the operators are applied to the identity map~$I$
only we want to solve
\begin{equation}
 \label{bch}
 \exp \U \exp \V = \exp \hat W
\end{equation}
for~$\hat W$. 

% This solution is not unique.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Non-identity map}

Note however, that 
$$
 \hat a I = \hat b I
$$
does not imply 
$$
 \hat a \phi = \hat b \phi \qquad \mbox{ for all } \phi \in \Phi
$$
or equivalently that
$$
 \hat a I = 0
$$
does not imply 
$$
 \hat a \phi = 0 \qquad \mbox{ for all } \phi \in \Phi.
$$
To see this,
consider the second order operator 
$$
 \hat a = \hat u \hat k 
$$
where~$\hat k$ is a first order operator corresponding 
to a constant vector field~$k$
$$
 k(x) = \mbox{const.}
$$
and~$\hat u$ is another first order operator corresponding to 
a vector field~$u$.
Then 
$$
 \hat a \phi = \hat u \hat k \phi = 
 (\phi'\cdot k)'\cdot u = \phi'' \cdot k \cdot u.
$$
Here $\phi'' \cdot k \cdot u$ is a map which in the point
$x\in R^n$ has the value 
$$
 \hat a \phi (x)  =
 \phi''(x)  \cdot k(x) \cdot u(x)
=\phi''(x)  \cdot \mbox{ const. } \cdot u(x).
$$
Then $\hat a\phi$
is nonzero for quadratic~$\phi$ while zero for linear~$\phi$
and thus for identity.
So we have
$$
 \hat a I = 0
$$
while nonzero result for nonlinear~$\phi$
$$ 
 \hat a \phi = \phi'' \cdot k \cdot u.
$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {BCH series}
  \label{sec.BCH}

Using the series (\ref{expop}) 
for the exponential of an operator 
and using the series 
\begin{equation}
 \label{logop}
 \log \hat u = \sum_{k=1}^\infty (-1)^{k-1} (\hat u - \hatI)^k / k
\end{equation}
for the logarithm
(assuming the series converges in accord with \cite{Olver}, \cite{Howe} etc.)
we get the solution 
of (\ref{bch}) for~$W$
in the form of the Hausdorff series \cite{Bourbaki}
\begin{equation}
 \label{hs}
 \hat W=\log(\exp\U\exp\V).
\end{equation}
% A new symbol $\U\vdash\!\!\!\dashv\V$ 
% (resembling the letter H for Hausdorff)
% is used to denote this series
% $$
%  \hat W=\U\vdash\!\!\!\dashv\V=\log(\exp\U\exp\V).
% $$
In the next section we present a detailed derivation
of a few low order terms of this series.

Reinsch in \cite{Reinsch} uses an interesting alternative approach
to get the terms in the Baker--Campbell--Hausdorff series,
see Appendix \ref{ApReinsch}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Second order expansion}

The second order expansion in~$p$ gives
$$
 \exp \U = \hatI + \U + \frac{1}{2} \U^2 
+ {\cal O}(p^3)
$$
$$
 \exp \V = \hatI + \V + \frac{1}{2} \V^2 + {\cal O}(p^3).
$$
Denoting 
$$
\hat Z =  \exp\U\exp\V-\hatI=
$$
$$
 =
 (\hatI+\U+\frac{1}{2}\U^2)
 (\hatI+\V+\frac{1}{2}\V^2)-\hatI 
 + {\cal O}(p^3)
 =
$$
$$
=
 \U+\V+\U\V+\frac{1}{2}\U^2+\frac{1}{2}\V^2
 + {\cal O}(p^3)
$$
we have
$$
 \hat Z^2=
(
 \U+\V+\U\V+\frac{1}{2}\U^2+\frac{1}{2}\V^2
)^2
 + {\cal O}(p^3)
=
$$
$$
 = 
 \U^2+\U\V+\V\U+\V^2
 + {\cal O}(p^3)
$$
and 
$$
\hat W=
%\U\vdash\!\!\!\dashv\V=
\log(\exp\U\exp\V)=
\hat Z - \frac{1}{2} \hat Z^2 
 + {\cal O}(p^3)
= 
$$
$$
=
 \U+\V+\U\V+\frac{1}{2}\U^2+\frac{1}{2}\V^2
 - \frac{1}{2}
(
 \U^2+\U\V+\V\U+\V^2
)
 + {\cal O}(p^3)
=
$$
$$
=
\U + \V + \frac{1}{2}(\U \V - \V\U)  
 + {\cal O}(p^3).
$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Commutator}

Denoting the {\bf commutator} of two operators
\begin{equation}
 \label{commutator}
 [\U,\V]=\U \V - \V\U
\end{equation}
we get
\begin{equation}
 \label{W2}
 \hat W = \U + \V + \frac{1}{2}[\U,\V]
 + {\cal O}(p^3).
\end{equation}
We say that~$\U$ and~$\V$ commute if 
$[\U,\V]=0$.

Note that 
\begin{equation}
\label{uvplusvu}
 [\U,\V] +  [\V,\U] = 0 
\end{equation}
for any two operators.
This property is called {\bf skew symmetry}.

Going back from 
$\U$, $\V$ and $\hat W$ 
to 
$\hat u$, $\hat v$ and $\hat w$ 
using (\ref{uU}-\ref{wW}) we get
$$
 \hat w = \frac{1}{2p}\hat W = 
\frac{1}{2}(\hat u +\hat v+\frac{p}{2}[\hat u,\hat v] 
+ {\cal O}(p^2)).
$$
Our ultimate goal is not the 
operator~$\hat w$ but rather the
vector field
$$
 w = \hat w I,
$$
i.e. the shortcut from the initial point~$x_0$ 
to the final point~$x_2$
of the zig+zag system.
To this purpose we must apply 
the commutator $[\hat u,\hat v]$
to the identity map~$I$
$$
 [\hat u,\hat v] I = 
 \hat u \hat v I - \hat v \hat u I =
  v'\cdot u-u'\cdot v.
$$
Thus a commutator of two vector fields is defined as
$$
 [u,v] = v'\cdot u-u'\cdot v
$$
and in coordinates
$$
 ([u,v](x))_j =  (([\hat u,\hat v]I)(x))_j = 
 \sum_{k=1}^n 
\frac{\partial v_j(x)}{\partial x_k} u_k(x) -
\frac{\partial u_j(x)}{\partial x_k} v_k(x).
$$
To avoid confusion, if~$\hat u$ and~$\hat v$ are operators,
so is their commutator $[\hat u, \hat v]$. Then 
$ [\hat u,\hat v]I$
is this operator applied to the identity map
from~$R^n$ to~$R^n$,
the result being
another map
from~$R^n$ to~$R^n$.
Then 
$([\hat u,\hat v]I)(x)$
is this map evaluated in the point~$x$, the result being
a point in~$R^n$.
And finally 
$ (([\hat u,\hat v]I)(x))_j  $
is the $j$-th coordinate of this point.

Then
\begin{equation}
 \label{w2BCH}
 w=\hat w I=\frac{u+v}{2}+p\frac{v'\cdot u-u'\cdot v}{4}
 + {\cal O}(p^2).
\end{equation}
This agrees with (\ref{w2Taylor}) completely.
It is always a pleasure to get the same result by 
two different ways.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Commuting vector fields}

We say that two vector fields~$u$, $v$ commute if
$$
 [u,v] = 0.
$$
This occurs if and only if 
their corresponding operators~$\hat u$, $\hat v$ commute.
In one dimension it is easy to find all the vector fields
that commute with a given vector field $u: R \rightarrow R$.
The condition 
$$
 v'(x) u(x) - u'(x) v(x) = 0
$$
is a linear differential equation for the function~$v(x)$
that can be solved by separation of variables to give
$$
 v(x) = k\  u(x) \qquad k \in R.
$$
Thus a one dimensional vector field 
commutes with its multiples only.
Alternatively, this can be shown by considering, for nonzero $f$
$$
 (\frac{g}{f})' = \frac{g' f-f' g}{f^2}
$$
and thus
$$
 [f,g] = f^2 (\frac{g}{f})' .
$$
Then $[f,g]=0$ if 
$$
 (\frac{g}{f})' = 0
$$
and this is for
$$
 \frac{g}{f} = k.
$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Higher order commutators}

We can generalize the notation used for a commutator
of two operators~$\hat A_1$, $\hat A_2$
$$
 [\hat A_1,\hat A_2]=\hat A_1\hat A_2-\hat A_2\hat A_1
$$
to commutators of three or more operators. 
For three operators we have
$$
 [\hat A_1,\hat A_2,\hat A_3]\stackrel{def}{=}
 [\hat A_1, [ \hat A_2,\hat A_3]]=
 \hat A_1 [ \hat A_2,\hat A_3] - 
  [ \hat A_2,\hat A_3] \hat A_1 =
$$
$$ 
 =  \hat A_1 (\hat A_2\hat A_3-\hat A_3\hat A_2) - 
 (\hat A_2\hat A_3-\hat A_3\hat A_2) \hat A_1  =
$$
$$
= 
 \hat A_1 \hat A_2\hat A_3-\hat A_1\hat A_3\hat A_2 - 
 \hat A_2\hat A_3\hat A_1+\hat A_3\hat A_2 \hat A_1
$$
Note that 
\begin{equation}
\label{jacobiI}
 [\hat A_1,\hat A_2, \hat A_3] + 
 [\hat A_2,\hat A_3, \hat A_1] + 
 [\hat A_3,\hat A_1, \hat A_2] = 0
\end{equation}
for any three operators.
This property is called {\bf Jacobi identity}.

Similarly,
we can define 
a commutator of order~$r$
(i.e.  
a commutator of~$r$ operators)
for any positive integer $r$ by
\begin{equation}
 \label{ncommutator}
 [\hat A_1, \hat A_2, \dots , \hat A_r] = 
 [\hat A_1,[\hat A_2,[\dots , \hat A_r]]]
\end{equation}
when we add for the sake of completeness
$$
 [\hat A] = \hat A.
$$
Note that a commutator of order~$r$ is linear in
each one of its arguments.

When each~$\hat A_i$ is either~$\U$ or~$\V$
(the two letters in a two-letter alphabet)
then we can express each commutator 
$ [\hat A_1, \hat A_2, \dots , \hat A_r]$
of order~$r$
as a certain linear combination of $2^r$ 
words 
(ordered $n$-tuples of elements 
that are either the symbol~$\U$ or the symbol~$\V$)
consisting of~$r$ letters.
These linear combinations 
can be conveniently express using
a square $2^r \times 2^r$ matrix~$M_r$.
We will need these matrices later, so we give 
examples of~$M_1$, $M_2$ and $M_3$ and then
we derive a recurrent relation for them.

For $r=1$ we have
\begin{eqnarray}
 \nonumber
 [\U] &=& \U \\
 \nonumber
 [\V] &=& \V 
\end{eqnarray}
and in the matrix form
$$
 \left(
   \begin{array}{c}{}
 {[\U]}
   \\
 {}[\V] 
   \end{array}
   \right)
  =
 \left(
  \begin{array}{r r}
    1&0 \\ 0&1    
  \end{array}
 \right)
 \cdot
 \left(
   \begin{array}{c}
 \U 
   \\
 \V 
   \end{array}
   \right)
  =
  M_1
 \cdot
 \left(
   \begin{array}{c}
 \U 
   \\
 \V 
   \end{array}
   \right).
$$

For $r=2$ we have
\begin{eqnarray}
 {}[\U,\U] &=& 0 \nonumber \\ 
 {}[\U,\V] &=& \U\V - \V\U \nonumber \\ 
 {}[\V,\U] &=& \V\U - \U \V \nonumber \\ 
 {}[\V,\V] &=& 0 \nonumber 
\end{eqnarray}
and in the matrix form
$$
 \left(
   \begin{array}{c}{}
 {}[\U,\U] \\ 
 {}[\U,\V] \\ 
 {}[\V,\U] \\ 
 {}[\V,\V] 
   \end{array}
   \right)
  =
 \left(
  \begin{array}{r r r r}
    0 & 0 & 0 & 0 \\
    0 & 1 &-1 & 0 \\
    0 &-1 & 1 & 0 \\
    0 & 0 & 0 & 0 
  \end{array}
 \right)
 \cdot
 \left(
   \begin{array}{c}{}
 \U\U \\ 
 \U\V \\ 
 \V\U \\ 
 \V\V 
   \end{array}
   \right)
  =
  M_2
 \cdot
 \left(
   \begin{array}{c}{}
 \U\U \\ 
 \U\V \\ 
 \V\U \\ 
 \V\V 
   \end{array}
   \right).
$$

For $r=3$ we have
\begin{eqnarray*}
 {}[\U,\U,\U] &=& 0 \nonumber \\ 
 {}[\U,\U,\V] &=& \U\U\V-2  \U\V\U+\V\U\U \\
 {}[\U,\V,\U] &=& - \U\U\V+2  \U\V\U-\V\U\U \\
 {}[\U,\V,\V] &=& 0 \\
 {}[\V,\U,\U] &=& 0 \\
 {}[\V,\U,\V] &=& - \U\V\V+2 \V\U\V-\V\V\U \\
 {}[\V,\V,\U] &=& \U\V\V-2 \V\U\V+\V\V\U \\
 {}[\V,\V,\V] &=& 0
\end{eqnarray*}
%
% here is how the above code was generated
%
%m[0,x_]=0;
%m[x_,0]=0;
%Clear[c];
%c1 = Tuples [ c[u,v], {3}];
%c[x_,y_,z_]:=c[x,c[y,z]];
%c[x_,y_]:=m[x,y]-m[y,x];
%m[x_,m[y_,z_]]:=m[x,y,z];
%m[m[x_,y_],z_]:=m[x,y,z];
%m[x__,y_+z_]:=m[x,y]+m[x,z];
%m[x__,y_-z_]:=m[x,y]-m[x,z];
%m[x_+y_,z__]:=m[x,z]+m[y,z];
%m[x_-y_,z__]:=m[x,z]-m[y,z];
%s = ToString[TeXForm[ColumnForm[c1]]];
%s2= StringReplace [ s, 
% {"m"->"","("->"",")"->"",","->"","u"->" \\U","v"->"\\V"}];
%Print[s2];
%
and in the matrix form
$$
 \left(
   \begin{array}{c}{}
 {}[\U,\U,\U] \\ 
 {}[\U,\U,\V] \\ 
 {}[\U,\V,\U] \\ 
 {}[\U,\V,\V] \\
 {}[\V,\U,\U] \\ 
 {}[\V,\U,\V] \\ 
 {}[\V,\V,\U] \\ 
 {}[\V,\V,\V] \\
   \end{array}
   \right)
  =
 \left(
 \begin{array}{r r r r r r r r}
  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  0 & 1 & -2 & 0 & 1 & 0 & 0 & 0 \\
  0 & -1 & 2 & 0 & -1 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  0 & 0 & 0 & -1 & 0 & 2 & -1 & 0 \\
  0 & 0 & 0 & 1 & 0 & -2 & 1 & 0 \\
  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
 \end{array}
 \right)
 \cdot
 \left(
   \begin{array}{c}{}
 {}\U\U\U \\ 
 {}\U\U\V \\ 
 {}\U\V\U \\ 
 {}\U\V\V \\
 {}\V\U\U \\ 
 {}\V\U\V \\ 
 {}\V\V\U \\ 
 {}\V\V\V \\
   \end{array}
   \right)
  =
$$
$$
=
  M_3
 \cdot
 \left(
   \begin{array}{c}{}
 {}\U\U\U \\ 
 {}\U\U\V \\ 
 {}\U\V\U \\ 
 {}\U\V\V \\
 {}\V\U\U \\ 
 {}\V\U\V \\ 
 {}\V\V\U \\ 
 {}\V\V\V \\
   \end{array}
   \right).
$$
In general
\begin{equation}
 \label{cMw}
  (\hat c_r)_k = \sum_{j=1}^{2^r} (M_r)_{k,j}\ (\hat w_r)_j
\end{equation}
where $(\hat c_r)_k$
is the $k$-th commutator of order~$r$,
and $(\hat w_r)_j$ is the $j$-th word consisting of $r$ letters.

To work with the matrices $M_r$ it is convenient to derive 
recurrent relations for them. The rule
$$
 [\hat A, \hat X] = \hat A \hat X - \hat X \hat A
$$
where $\hat A$ is either $\U$ or $\V$ and
$\hat X$ is a commutator corresponding to a row of the matrix $M_r$,
implies that each row in the upper half of the matrix $M_{r+1}$
(i.e. a row corresponding to a commutator beginning with~$\U$)
is a difference of two rows:
\begin{itemize}
\item
 a row of matrix~$M_r$ corresponding to the operator~$\hat X$
 padded by zeroes on the right and
 \item
 this row with zeroes inserted to even positions
\end{itemize}
and that each row in the bottom half of the matrix~$M_{r+1}$ 
(i.e. a row corresponding to a commutator beginning with~$\V$)
is a difference of two rows:
\begin{itemize}
\item
 a row of matrix~$M_r$ corresponding to the operator~$\hat X$
 padded by zeroes on the left and
 \item
 this row with zeroes inserted to odd positions
\end{itemize}
namely, 
\begin{equation}
 (M_{r+1})_{i,j}=
 \left\{  % } aby byly v paru
  \begin{array}{r c l l}
(M_r)_{i,j}&-&(M_r)_{i,\frac{j+1}{2}}      &\mbox{ for $i\le 2^r,\ j\le 2^r$,   $j$ odd }\\
(M_r)_{i,j}& &                             &\mbox{ for $i\le 2^r,\ j\le 2^r$,   $j$ even}\\
           &-&(M_r)_{i,\frac{j+1}{2}}      &\mbox{ for $i\le 2^r,\ j >  2^r$,   $j$ odd }\\
           &0&                             &\mbox{ for $i\le 2^r,\ j >  2^r$,   $j$ even}\\
           &0&                             &\mbox{ for $i >  2^r,\ j\le 2^r$,   $j$ odd }\\
           &-&(M_r)_{i-2^r,\frac{j-2^r}{2}}&\mbox{ for $i >  2^r,\ j\le 2^r$,   $j$ even}\\
(M_r)_{i-2^r,j-2^r}& &                             &\mbox{ for $i>2^r,\ j>2^r$, $j$ odd }\\
(M_r)_{i-2^r,j-2^r}&-&(M_r)_{i-2^r,\frac{j-2^r}{2}}&\mbox{ for $i>2^r,\ j>2^r$, $j$ even}\\
       \end{array}
       \right.
\end{equation}
% Readers interested in details can check the algorithm written
% in {\it Mathematica} in the Appendix.

This can be conveniently computed by the 
following {\it Mathematica} code:

{
% \footnotesize
\scriptsize
\begin{verbatim}
m={{1,0},{0,1}};
For[r=2,r<=3,r++,
 n=2^(r-1);
 m= 
 Table[PadRight[m[[i]],2n]-Insert[m[[i]],0,Table[{j+1},{j,n}]],{i,n}]
 ~Join~
 Table[PadLeft[m[[i]],2n]-Insert[m[[i]],0,Table[{j},{j,n}]],{i,n}];
 Print["r=",r," m=",MatrixForm[m]];
];
\end{verbatim}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Higher order expansion in words}

The number $N_w$ of words
in
$$ 
 \hat W=
% \U\vdash\!\!\!\dashv\V=
 \log(\exp\U\exp\V)
$$
of order~$r$ 
is too large to be listed here
even for moderately large~$r$,
see Tab.~\ref{tabNcomm}.
To give a picture of the series we present the expansion
up to order 6 only
$$
 \hat W=
 \hat W_{1}+
 \hat W_{2}+
 \hat W_{3}+
 \hat W_{4}+
 \hat W_{5}+
 \hat W_{6}+
 {\cal O}(p^7)
$$
where
\begin{eqnarray*}
 \hat W_{1}&=&\U + \V
 \\ 
\hat W_{2}&=&\frac{1}{2}(\U\V-\V\U)
 \\ 
 \hat W_{3}&=&\frac{1}{12} ({ \U\U\V}-2 { \U\V\U}+{
   \U\V\V}+{\V\U\U}-2 {\V\U\V}+{\hat
   V\V\U})
 \\ 
 \hat W_{4}&=&\frac{1}{24} ({ \U\U\V\V}-2 { \U\V \hat
   U\V}+2 {\V\U\V\U}-{\V\V\U\U})
 \\ 
 \hat W_{5}&=&\frac{1}{6! } (-
  { \U\U\U\U\V}+ 4 { \U\U\U\V\U}+ 4 { \U\U\U\V\V}- 6 { \U\U\V\U\U}-
 \\ & & 
-6 { \U\U\V\U\V}- 6 { \U\U\V\V\U}+ 4 { \U\U\V\V\V}+4 { \U\V\U\U\U}-
 \\ & & 
-6 { \U\V\U\U\V}+ 24 { \U\V\U\V\U}- 6 { \U\V\U\V\V}-6 { \U\V\V\U\U}-
 \\ & & 
-6 { \U\V\V\U\V}+4 { \U\V\V\V\U}- { \U\V\V\V\V}-{\V\U\U\U\U}+
 \\ & & 
+4 {\V\U\U\U\V}-6 {\V\U\U\V\U}- 6 {\V\U\U\V\V}-6 {\V\U\V\U\U}+
 \\ & & 
+24 {\V\U\V\U\V}-6 {\V\U\V\V\U}+ 4 {\V\U\V\V\V}+4 {\V\V\U\U\U}-
 \\ & & 
-6 {\V\V\U\U\V}-6 {\V\V\U\V\U}- 6 {\V\V\U\V\V}+4 {\V\V\V\U\U}+
 \\ & & 
+4 {\V\V\V\U\V}-{\V\V\V\V\U})
 \\ 
 \hat W_{6}&=&
 \frac{1}{2 \cdot 6!} (
-{ \U\U\U\U\V\V}+4 { \U\U\U\V\U\V}+4 { \U\U\U\V\V\V}-6 { \U\U\V\U\U\V}-
 \\ & & 
-6 { \U\U\V\U\V\V}-6 { \U\U\V\V\U\V}-{ \U\U\V\V\V\V}+4 { \U\V\U\U\U\V}-
 \\ & & 
-6 { \U\V\U\U\V\V}+24 { \U\V \U\V\U\V}+4 { \U\V\U\V\V\V}-6 { \U\V\V\U\U\V}-
 \\ & & 
-6 { \U\V\V\U\V\V}+4 { \U\V\V\V\U\V}-4 {\V\U\U \U\V\U}+6 {\V\U\U\V\U\U}+
 \\ & & 
+6 {\V\U\U\V\V\U}-4 {\V\U\V\U\U\U}-24 {\V\U\V\U\V\U}+6 {\V \U\V\V\U\U}-
 \\ & & 
-4 {\V\U\V\V\V \U}+{\V\V\U\U\U\U}+6 {\V\V\U\U\V\U}+6 {\V\V\U\V\U\U}+
 \\ & & 
+6 {\V\V \U\V\V\U}-4 {\V\V\V\U\U\U}-4 {\V\V\V\U\V\U}+{\V\V\V\V \U\U})
\end{eqnarray*}

% here is how the above LaTeX code was generated by Mathematica
%(*************************************)
%job[r_]:= ( 
%  n = r+1;
%(*************************************)
%m[a_,b_]:=Reap[
% For[p=0,p<n,p++,
%  For[q=0,q<2^p,q++,
%   w = 0;
%   For[i=0,i<=p,i++, 
%    dw = a[[IntegerPart[q/2^(p-i)]  + 2^i]] * 
%         b[[BitAnd[q,2^(p-i)-1] + 2^(p-i)]];
%    w += dw;
%   ];
%   Sow[w];
%  ];
% ];
%]//Rest//First//First;
%(*************************************)
%log[x_]:=(
% xx  = ReplacePart[x,x[[1]]-1,1];
% rr = xx;
% dx = xx;
% For[j=2,j<n,j++,
%  dx /= -j/(j-1);
%  dx = m[dx,xx];
%  rr= rr+ dx;  
% ];
%rr);
%(*************************************)
% ex=Table[{1/i!}~Join~Table[0,{2^i-1}],{i,0,r}]//Flatten;
% ey=Table[Table[0,{2^i-1}]~Join~{1/i!},{i,0,r}]//Flatten;
% exey = m[ex,ey]; 
% z = log[exey];
%(*************************************)
%t[1]={{1,0},{0,1}};
%For[s=2,s<=r,s++,
%    coefw[s] = Take[z,{2^s,2^(s+1)-1}];
%    tu=Tuples[sj["u","v"],s]/.sj->StringJoin//ToExpression;
%    mm[s]=coefw[s].tu//Simplify//TeXForm//ToString;
%    vy=StringReplace[mm[s],{"\\text"->"","u"->" \\U","v"->"\\V"}];
%    vy="$$\n \\hat W_{"<>ToString[s]<>"}="<>vy<>"\n$$";
%    Print[vy];
%];
%);
%(*************************************)
%jobb[r_]:=(
% cpu = First[Timing[job[r]]]/Second//Chop;
%);
%(*************************************)
%jobb[6]; 
%(*************************************)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% order words rank commutators mem real denominator denom/r!
%  2      2     1      1        5                 2   1
%  3      6     2      2        5                12   2
%  4      4     3      1        5                24   1
%  5     30     6      6        5               720   6
%  6     28     9      5        5              1440   2
%  7    126    18     18        5   1   s     30240   6
%  8    124    30     17        5   2.5 s    120960   3
%  9    390    56     55        6   5.6 s   3628800  10
% 10    388    99     54       15  15   s   7257600   2
% 11   2046   186    181       60  88   s 239500800   6
% 12   2044   335    180      180 677   s 958003200   2



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                      bch12f.m            14.8.2005
%
% r=12  mem  = 132 MB
%
% pekrin: real= 3:43 cpu=3:27   r=13 killed after 37 minutes
% acer:   real=10:50 cpu=6:16
% 
%!!a
%(*************************************)
%job[r_]:= ( 
%  n = r+1;
%(*************************************)
%m[a_,b_]:=Reap[
% For[p=0,p<n,p++,
%  For[q=0,q<2^p,q++,
%   w = 0;
%   For[i=0,i<=p,i++, 
%    dw = a[[IntegerPart[q/2^(p-i)]  + 2^i]] * 
%         b[[BitAnd[q,2^(p-i)-1] + 2^(p-i)]];
%    w += dw;
%   ];
%   Sow[w];
%  ];
% ];
%]//Rest//First//First;
%(*************************************)
%log[x_]:=(
% xx  = ReplacePart[x,x[[1]]-1,1];
% rr = xx;
% dx = xx;
% For[j=2,j<n,j++,
%  dx /= -j/(j-1);
%  dx = m[dx,xx];
%  rr= rr+ dx;  
% ];
%rr);
%(*************************************)
% ex=Table[{1/i!}~Join~Table[0,{2^i-1}],{i,0,r}]//Flatten;
% ey=Table[Table[0,{2^i-1}]~Join~{1/i!},{i,0,r}]//Flatten;
% exey = m[ex,ey]; 
% z = log[exey];
%(*************************************)
%t      = SparseArray[{{0,1,-1,0}}];
%zeroes = SparseArray[{0}];
%For[s=3,s<=r,s++,
%    coefw[s] = Take[z,{2^s,2^(s+1)-1}];
%    If[s==r,Clear[z,ex,ey,exey,xx,dx,rr]];
%    coefwnz=Select[coefw[s],#!=0&];
%    llnz = Length[coefwnz];
%    nn=2^(s-1);
%    z1  = Table [{i  },{i,nn}];
%    z2  = Table [{i+1},{i,nn}];
%    t=(PadRight[t,{nn/4,2nn}] - Transpose[Insert[Transpose[t],zeroes,z2]])
%         ~Join~ (
%       PadLeft [t,{nn/4,2nn}] - Transpose[Insert[Transpose[t],zeroes,z1]]);
%    rank = MatrixRank[t];
%    t=Transpose[t];
%    sm = LinearSolve[t,coefw[s]];
%    t=Transpose[t];
%    cc = Tuples[c[u,v],s];
%    cc = Take[cc,{2,-1,4}];
%    bb[s] = sm.cc;
%    hh = Head[bb[s]];
%    lc = -1;
%    If[hh==Plus , lc = Length[bb[s]]];
%    If[hh==Times, lc = 1];
%    bb[s]=Simplify[bb[s]];
%    mem = MaxMemoryUsed[]/1024.0^2;
%Print["bb=",bb[s]];
%Print["s=",s," llnz=",llnz," rank=",rank," lc=",lc," mem=",mem];
%    zeroes = Join[zeroes,zeroes];
%];
%);
%(*************************************)
%jobb[r_]:=(
% cpu = First[Timing[job[r]]]/Second//Chop;
% Print["r=",r," cpu=",cpu];
%);
%(*************************************)
%jobb[12]; 
%(*************************************)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table}[htb]
\begin{center}
\begin{tabular}{|r|r r r r r|}
  \hline   
$r$ &$N_{tot}$&$N_w$ &rank($M_r$) &$N_{com}$ & $d/r!$ \\
  \hline    
  1 &       2 &    2 &   2 &    2 &  1 \\
  2 &       4 &    2 &   1 &    1 &  1 \\
  3 &       8 &    6 &   2 &    2 &  2 \\
  4 &      16 &    4 &   3 &    1 &  1 \\
  5 &      32 &   30 &   6 &    6 &  6 \\  
  6 &      64 &   28 &   9 &    4 &  2 \\
  7 &     128 &  126 &  18 &   18 &  6 \\
  8 &     256 &  124 &  30 &   13 &  3 \\
  9 &     512 &  390 &  56 &        38 & 10 \\
 10 &    1024 &  388 &  99 &        31 &  2 \\
 11 &    2048 & 2046 & 186 & $\le$ 181 &  6 \\
 12 &    4096 & 2044 & 335 & $\le$ 180 &  2 \\
 \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
  \hline
\end{tabular}
\end{center}
 \caption{
 This table shows for each order~$r$ between 1 and 12 
 the total number $N_{tot}=2^r$ of all 
 possible words consisting of $r$ letters ($\U$ or $\V$),
 the number $N_w$ of words appearing 
 in the $r$-th order terms of BCH series, 
 the rank of the matrix~$M_r$ relating the commutators 
 and the words consisting of $r$ letters,
 the least possible number $N_{com}$ of commutators 
 in the $r$-th order terms of BCH,
 and the common denominator $d$ 
 (divided by factorial of $r$) 
 of the terms of order $r$ of the BCH series.
% These results were computed by {\it Mathematica} program
% given in Appendix~\ref{AppM} and by C programs given 
% in Appendices \ref{AppCbch} and \ref{AppCsparsest}.
}
 \label{tabNcomm}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Higher order expansion in commutators}

It is obvious that the second order term can be written using 
a commutator
$$
 \hat W_2 = \frac{1}{2}[\U,\V].
$$
It is still quite easy to find that the third order and the fourth order 
terms can be written using commutators as well
$$
\hat W_{3}=
\frac{1}{12} ([\U,\U,\V]-[\V,\U,\V])
$$
$$
\hat W_{4}=
-\frac{1}{24} [\U,\V,\U,\V].
$$
It is a nontrivial fact that all the terms of 
$ \hat W =  \log (\exp\U \exp\V)$
% \U\vdash\!\!\!\dashv\V$
% of order 2 or more 
can be expressed as linear combinations of commutators only.
% semm
Seven different proofs of this fact are given in
\cite{Baker},
\cite{Cartier},
\cite{Djokovic},
\cite{Dynkin},
\cite{Hochschild},
\cite{Pejas},
\cite{Varadarajan}.
The main idea of the proof by Djokovic 
\cite{Djokovic}
is given in Appendix.
The consequence of this is that 
$ \log(\exp\U\exp\V) = \U + \V$
% $ \U\vdash\!\!\!\dashv\V = \U + \V$
if and only if $\U$ and $\V$ commute.
For a short introduction to the Lie theory see \cite{Howe}.

The formula giving 
$$
 \hat W = \log(\exp\U\exp\V)
$$
in commutators of $\U$ and $\V$ is called the 
Baker -- Campbell -- Hausdorff formula.
It is named after the British mathematician 
Henry Frederick Baker (1866-1956),
after the Irish mathematician John Edward Campbell (1862-1924)
and after the German mathematician Felix Hausdorff (1868-1942).

\begin{table}[htb]
\begin{center}
\begin{tabular}{|l|l|}
  \hline   
\multicolumn{1}{|c|}{words} & \multicolumn{1}{c|}{commutators} 
\\
  \hline   
$(\hat w_1)_1 = \U    $ & $(\hat c_1)_1 = [\U]      $\\
$(\hat w_1)_2 = \V    $ & $(\hat c_1)_2 = [\V]      $\\
$(\hat w_2)_1 = \U\U  $ & $(\hat c_2)_1 = [\U,\U]   $\\
$(\hat w_2)_2 = \U\V  $ & $(\hat c_2)_2 = [\U,\V]   $\\
$(\hat w_2)_3 = \V\U  $ & $(\hat c_2)_3 = [\V,\U]   $\\
$(\hat w_2)_4 = \V\V  $ & $(\hat c_2)_4 = [\V,\V]   $\\
$(\hat w_3)_1 = \U\U\U$ & $(\hat c_3)_1 = [\U,\U,\U]$\\
$(\hat w_3)_2 = \U\U\V$ & $(\hat c_3)_2 = [\U,\U,\V]$\\
$(\hat w_3)_3 = \U\V\U$ & $(\hat c_3)_3 = [\U,\V,\U]$\\
$(\hat w_3)_4 = \U\V\V$ & $(\hat c_3)_4 = [\U,\V,\V]$\\
$(\hat w_3)_5 = \V\U\U$ & $(\hat c_3)_5 = [\V,\U,\U]$\\
$(\hat w_3)_6 = \V\U\V$ & $(\hat c_3)_6 = [\V,\U,\V]$\\
$(\hat w_3)_7 = \V\V\U$ & $(\hat c_3)_7 = [\V,\V,\U]$\\
$(\hat w_3)_8 = \V\V\V$ & $(\hat c_3)_8 = [\V,\V,\V]$\\
\multicolumn{1}{|c|}{\vdots} & \multicolumn{1}{c|}{\vdots}\\ 
  \hline
\end{tabular}
\end{center}
 \caption{ 
    This table explains the symbols for words and commutators 
    used in the text for order $r$ between 1 and 3.
    To avoid confusion, note that $\hat w$ is a word, 
    while $\hat W$ is the BCH series.
     }
 \label{wandc}
\end{table}

To find the explicit form 
of the $r$-th order term in the BCH series
as a linear combination of commutators $(\hat c_r)_k$
with coefficients $(b_r)_k$ (to be found),
 namely
\begin{equation} 
 \label{Wbc}
 \hat W_r = \sum_{k=1}^{2^r} (b_r)_k\ (\hat c_r)_k
\end{equation} 
supposing we know the coefficients~$(a_r)_j$ in
the linear combinations 
of the words 
$(\hat w_r)_j$
giving~$\hat W_r$, namely
\begin{equation} 
 \label{Waw}
 \hat W_r = \sum_{j=1}^{2^r} (a_r)_j\ (\hat w_r)_j
\end{equation} 
all we need to do is to put 
(\ref{cMw}) into (\ref{Wbc}) and to compare it with (\ref{Waw}).
By comparing coefficients in front of individual words
$(\hat w_r)_j$ on the left and on the right we get
$$
 \sum_{k=1}^{2^r} (b_r)_k\ (M_r)_{k,j} = (a_r)_j
$$
which in matrix form is
\begin{equation} 
 \label{Mba}
   M_r^T \cdot \vec b_r = \vec a_r
\end{equation} 
(here $T$ means transpose).
This is the linear equation to give the coefficients 
$(b_r)_k$ in front of commutators in the BCH series.
With the exception of order $r=1$
this linear equation has not a unique solution.
The matrix~$M_r$ is singular, its rank is less than its size~$2^r$,
see Tab.~\ref{tabNcomm}.
But the solution exists.
As an illustration we present the results of order up to eight here
\begin{eqnarray*}
\hat W_{5}&=&
\frac{1}{6! }
(-[\U,\U,\U,\U,\V]-6[\U,\U,\V,\U,\V]-2[\U,\V,\V,\U,\V]+
\\ & & 
+2 [\V,\U,\U,\U,\V]+6[\V,\U,\V,\U,\V]+[\V,\V,\V,\U,\V])
\\
 \hat W_{6}&=&\frac{1}{2 \cdot 6!}
([\U,\V,\V,\V,\U,\V]+ [\V,\U,\U,\U,\U,\V]+\\ & &
+6[\V,\U,\U,\V,\U,\V] -2[\V,\V,\U,\U,\U,\V])
\\
\hat W_{7}&=&
\frac{1}{6\cdot 7!}
([\U,\U,\U,\U,\U,\U,\V]-3[\U,\U,\U,\U,\V,\U,\V]-
\\ & & 
-4[\U,\U,\U,\V,\V,\U,\V]+12[\U,\U,\V,\U,\U,\U,\V]+
\\ & & 
+9[\U,\U,\V,\U,\V,\U,\V]-9[\U,\U,\V,\V,\V,\U,\V]+
\\ & & 
+18[\U,\V,\U,\U,\V,\U,\V]+24 [\U,\V,\U,\V,\V,\U,\V]+
\\ & & 
+3[\U,\V,\V,\V,\V,\U,\V]-3 [\V,\U,\U,\U,\U,\U,\V]-
\\ & & 
-9[\V,\U,\U,\U,\V,\U,\V]+4 [\V,\U,\U,\V,\V,\U,\V]-
\\ & & 
-6 [\V,\U,\V,\U,\U,\U,\V]-30 [\V,\U,\V,\U,\V,\U,\V]-
\\ & & 
-3 [\V,\U,\V,\V,\V,\U,\V]+3 [\V,\V,\U,\U,\V,\U,\V]-
\\ & & 
-6 [\V,\V,\U,\V,\V,\U,\V]-[\V,\V,\V,\V,\V,\U,\V])
\\
\hat W_{8}&=&
\frac{1}{3\cdot 8!}(
6[\U,\U,\V,\V,\V,\V,\U,\V] -24[\U,\V,\V,\U,\V,\V,\U,\V]+
\\ & &
+6[\U,\V,\V,\V,\V,\U,\U,\V] -2[\U,\V,\V,\V,\V,\V,\U,\V]-
\\ & &
-2[\V,\U,\U,\U,\U,\U,\U,\V] -12[\V,\U,\U,\V,\U,\U,\U,\V]+
\\ & &
+6[\V,\U,\U,\V,\U,\V,\U,\V] -6[\V,\U,\V,\U,\U,\U,\U,\V]-
\\ & &
-60[\V,\U,\V,\U,\U,\V,\U,\V]+ 8[\V,\U,\V,\V,\U,\U,\U,\V]+
\\ & &
+6[\V,\V,\U,\U,\U,\U,\U,\V]+ 48[\V,\V,\U,\V,\U,\U,\U,\V]-
\\ & &
-25[\V,\V,\V,\U,\U,\U,\U,\V]
).
\end{eqnarray*}

% the above code was generated by this program
%(*************************************)
%job[r_]:= ( 
%  n = r+1;
%(*************************************)
%m[a_,b_]:=Reap[
% For[p=0,p<n,p++,
%  For[q=0,q<2^p,q++,
%   w = 0;
%   For[i=0,i<=p,i++, 
%    dw = a[[IntegerPart[q/2^(p-i)]  + 2^i]] * 
%         b[[BitAnd[q,2^(p-i)-1] + 2^(p-i)]];
%    w += dw;
%   ];
%   Sow[w];
%  ];
% ];
%]//Rest//First//First;
%(*************************************)
%log[x_]:=(
% xx  = ReplacePart[x,x[[1]]-1,1];
% rr = xx;
% dx = xx;
% For[j=2,j<n,j++,
%  dx /= -j/(j-1);
%  dx = m[dx,xx];
%  rr= rr+ dx;  
% ];
%rr);
%(*************************************)
% ex=Table[{1/i!}~Join~Table[0,{2^i-1}],{i,0,r}]//Flatten;
% ey=Table[Table[0,{2^i-1}]~Join~{1/i!},{i,0,r}]//Flatten;
% exey = m[ex,ey]; 
% z = log[exey];
%(*************************************)
%t[1]={{1,0},{0,1}};
%For[s=2,s<=r,s++,
%    coefw[s] = Take[z,{2^s,2^(s+1)-1}];
%    coefwnz=Select[coefw[s],#!=0&];
%    llnz = Length[coefwnz];
%    nn=2^(s-1);
%    z1 = Table [{i  },{i,nn}];
%    z2 = Table [{i+1},{i,nn}];
%    ax[x_]:= PadRight[x,2nn];
%    xa[x_]:= Insert[x,0,z2];
%    bx[x_]:= PadLeft [x,2nn];
%    xb[x_]:= Insert[x,0,z1];
%    ut[x_]:= Table[ ax[x[[i]]]-xa[x[[i]]],{i,nn}]~Join~
%             Table[ bx[x[[i]]]-xb[x[[i]]],{i,nn}];
%    t[s] = ut[t[s-1]];
%    M[s] = Take[t[s],{2,-1,4}];
%    rank = MatrixRank[M[s]];
%    If[s==r, Clear[t]];
%    sm[s]= LinearSolve[Transpose[M[s]],coefw[s]];
%    Clear[c];
%    cc = Tuples[c[u,v],s];
%    cc = Take[cc,{2,-1,4}];
%    c[x__]:=ff[f1,x,f2];
%    bb[s] = sm[s].cc;
%    bb[s]=Simplify[bb[s]];
%    qq = TeXForm[bb[s]];
%    qq = ToString[qq];
%    qq=StringReplace[qq,{"\\text"->"","f1"->"[", "f2"->"]",
%       "{ff}"->"","u"->"\\U","v"->"\\V","[,"->"[", ",]"->"]"}];
%    qq=StringReplace[qq,{"{[}"->"[", "{]}"-> "]"}];
%    qq=StringReplace[qq,{"[,"->"[", ",]"-> "]"}];
%    
%    q1 = "\\hat W_{"<>ToString[s]<>"}=";
%Print["$$"];
%Print[q1];
%Print[qq];
%Print["$$"];
%Print[];
%];
%);
%(*************************************)
%jobb[r_]:=(
% cpu = First[Timing[job[r]]]/Second//Chop;
%);
%(*************************************)
%jobb[8]; 
%(*************************************)
%

%\begin{equation}
% \label{W4}
% \hat W = 
% \U + \V 
% + \frac{1}{2}[\U,\V]
% + \frac{1}{12} [\U - \V, [\U,\V]]
% - \frac{1}{24}[\U,[\V,[\U,\V]]]
% + {\cal O}(p^5).
%\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Non-uniqueness of results in commutators}

Unlike the results 
for $\hat W_r$ 
in words, the results
for $\hat W_r$ 
in commutators are not unique with the only exception of $r=1$.
This is the consequence of the fact that the rank of 
the matrix $M_r$ that appears in 
(\ref{Mba}) 
is less 
than its size $2^r$ for $r>1$.

The general solution to (\ref{Mba}) can be written as a sum
of a particular solution plus the general solution of the 
corresponding homogeneous equation. In other words we can
add any vector from the kernel of the matrix $M_r^T$ 
(i.e. a vector that is mapped to zero)
and we get another solution of the linear system.
But this another solution corresponds to the same term
in the BCH series.

To illustrate this, let us consider $r=2$.
One possible basis of the kernel of $M_r^T$ is 
\begin{equation}
 \label{threebv}
 (1,0,0,0)^T, \quad 
 (0,1,1,0)^T, \quad 
 (0,0,0,1)^T.
\end{equation}
Then the general solution to (\ref{Mba}) is
$$
 \vec b_2 = 
\left(
  \begin{array}{c}
   0 \\ \frac{1}{2} \\ 0 \\ 0
  \end{array}
 \right)
+k_1\left(\begin{array}{c}1\\0\\0\\0\end{array}\right)
+k_2\left(\begin{array}{c}0\\1\\1\\0\end{array}\right)
+k_3\left(\begin{array}{c}0\\0\\0\\1\end{array}\right),
$$
with $k_1, k_2, k_3 \in R$.
The corresponding second order term in the BCH series is
\begin{eqnarray*}
 \hat W_2 & = & 
   \frac{1}{2}(\hat c_2)_2+k_1(\hat c_2)_1+
   k_2((\hat c_2)_2+(\hat c_2)_3)+k_3(\hat c_2)_4 =
   \\
  &=&  \frac{1}{2}[\U,\V]+k_1[\U,\U]+k_2([\U,\V]+[\V,\U])+k_3[\V,\V].
\end{eqnarray*}
In this case (for $r=2$) 
it is trivial to see that the new terms add nothing new
but they may change the form of the result. 
As an example consider 
\begin{eqnarray*}
  k_1 &=& 0,\\
  k_2 &=& -\frac{1}{2},\\
  k_3 &=& 0.
\end{eqnarray*}
Then we get 
$$
 \hat W_2 = -\frac{1}{2} [\V,\U]
$$
which is equal to the result given above, just in a different form.
The form of the result may be changed in two ways:
\begin{itemize}
\item
 one term may be replaced by another, e.g. 
$[\U,\V]$ may be replaced by $-[\V,\U]$, or
\item
 the number of terms may be changed.
\end{itemize}
We want to find the {\bf sparsest} solution to (\ref{Mba}), 
i.e. the solution that has the least possible number
of nonzero terms.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Identities}
 \label{sec.iden}

For any order $r$ 
each vector in the kernel of the matrix $M_r^T$
corresponds to a certain identity.
Three vectors in the kernel of $M_2^T$ given in (\ref{threebv})
correspond 
to three identities
$$
  (\hat c_2)_1 = [\U,\U] = 0
$$
$$
 (\hat c_2)_2 +  (\hat c_2)_3 = [\U,\V] + [\V,\U]  = 0
$$
and 
$$
 (\hat c_2)_4 =  [\V,\V] = 0.
$$
This is a direct consequence of the skew symmetry of the commutator.

Oteo in \cite{Oteo} finds some of these identities by 
comparison of the BCH series computed by various methods.
Our approach is more systematic.
Based on the kernel of the matrices $M_r$ 
we are able to give both the total number 
and the form of these identities.

Each such identity says that 
a certain linear combination of a certain number
(we call this number the length $L$ of the identity) of commutators 
is equal to zero.
Table~\ref{tabLR} shows the number of independent identities
and the maximal length of an non-decomposable 
identity for order $r$ between $1$ and $12$.
We say that an identity is non-decomposable if none of its proper parts
is an identity itself.
E.g. the identity
$$
 [\U,\V] + [\V,\U]  = 0
$$
is non-decomposable, while the identity
$$
  [\U,\U] + [\V,\V] = 0
$$
is decomposable, because 
it can be decomposed into two identities, namely
$$
  [\U,\U] = 0
$$
and
$$
  [\V,\V] = 0.
$$

\begin{table}[htb]
\begin{center}
\begin{tabular}{|r|r r |}
  \hline   
$r$ & $N_{rules}$ & $L_{max}$ \\
  \hline    
  1 &       0 &       \\
  2 &       3 &    2  \\
  3 &       6 &    2  \\
  4 &      13 &    2  \\
  5 &      26 &    2  \\  
  6 &      55 &    4  \\
  7 &     110 &    4  \\
  8 &     226 &    8  \\
  9 &     456 &   15  \\
 10 &     925 &   26  \\
 11 &    1862 &   43  \\
 12 &    3761 &   76  \\
 \vdots & \vdots & \vdots \\
  \hline
\end{tabular}
\end{center}
 \caption{
 This table shows for each order~$r$ between 1 and 12 
 the number $N_{rules}$ of independent rules (identities) of order $r$ 
 (which is the dimension of the kernel of the matrix $M_r^T$,
  cf. Table \ref{tabNcomm})
 and the maximal length $L_{max}$ of a non-decomposable rule of order $r$.
 The field of $L_{max}$ for $r=1$ is blank because there are no rules
for $r=1$, the kernel of $M_1$ is just a zero vector. 
}
 \label{tabLR}
\end{table}

For higher order~$r$ it may be less obvious whether two forms
of the result are equal or not.

It was a surprise for me that the kernel of the matrix $M_4^T$
contains the vector
$$
 (0,0,0,0,0,1,0,0,0,-1,0,0,0,0,0,0)
$$
which corresponds to the identity
$$
 (\hat c_4)_6  - (\hat c_4)_{10}=0
$$
i.e.
$$
 [\U,\V,\U,\V] - [\V,\U,\U,\V] = 0.
$$
I expected the plus sign instead of the minus sign,
in analogy with 
$$
 [\U,\V] + [\V,\U] = 0.
$$
But the above identity is indeed correct, as can be shown
by direct evaluation of
$$
 [\U,\V,\U,\V] = [\U,[\V,[\U,\V]]] = 
 [\U,[\V,\U\V-\V\U]] = 
$$
$$
 =[\U,\V\U\V-\V\V\U-\U\V\V+\V\U\V]=
$$
$$
=\U\V\U\V-\U\V\V\U-\U\U\V\V+\U\V\U\V-\V\U\V\U+\V\V\U\U+\U\V\V\U-\V\U\V\U = 
$$
$$
=-\U\U\V\V+2\U\V\U\V-2\V\U\V\U+\V\V\U\U
$$
and similarly for 
$$
 [\V,\U,\U,\V] = \dots = -\U\U\V\V+2\U\V\U\V-2\V\U\V\U+\V\V\U\U.
$$
As a consequence both
$$
\hat W_{4}=
-\frac{1}{24} ([\U,\V,\U,\V])
$$
and
$$
\hat W_{4}=
-\frac{1}{24} ([\V,\U,\U,\V])
$$
are correct.
This changes the form of the result but not its sparsity. 
Each form has exactly one term.

As an illustration in 
Fig.\,\ref{obrBJ} 
we present a basis of the kernel of $M_6^T$ shown
as a $55 \times 64 $ matrix. 
Here positive elements are shown as red squares, negative 
elements are shown as blue squares (we adopt this convention
from electronics) and zero elements as blank squares.
The higher the absolute value of an element, the higher the 
saturation of the color.
As an example, the $13$-th vector in this basis is
$$
(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -3, 
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, $$ $$ 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 
 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)^T
$$
corresponding to the identity
$$
   [\U, \U, \V, \V, \U, \V] - 
 3 [\U, \V, \U, \V, \U, \V] + 
$$
$$
 +3 [\V, \U, \U, \V, \U, \V] + 
   [\V, \V, \U, \U, \V, \U] = 0 .
$$
Note that there is just one minus sign and three plus signs
in this identity, as seen also in
Fig. \ref{obrBJ} where the $13$-th row contains one blue square
and 3 red squares.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=\figsize]{bj.ps}
\end{center}
 \caption{
  One possible basis of the kernel of the matrix $M_6^T$.
  The matrix $M_6^T$ is a $64\times 64$ matrix. Its rank is 9,
  so the dimension of the kernel is $55$. 
  Each vector in the kernel has 64 components.
  Thus there are 55 rows 
  each consisting of 64 squares in this picture.
  Each row depicts one vector in the basis of the kernel.
  Red squares represent positive components while blue squares
  represent negative ones. The saturation of the color shows
  the absolute value of the corresponding element.
  Zero components are shown as blank squares (most of the components
  are zero).
 }\label{obrBJ}
\end{figure}
% the figure was made by this Mathematica code:
%ap[x_]:=(
% ArrayPlot[x,
%  ColorFunctionScaling -> True,
%  ColorFunction->(If[#==1/2,Hue[1,0,1],
%         If[#>1/2, Hue[1,2Abs[#-1/2],1],  Hue[0.6,1,1]]]&)
% ]);
%r = 6;
%n = r+1;
%t = {{1,0},{0,1}};
%For[s=2,s<=r,s++,
%  nn=2^(s-1);
%  z1 = Table [{i  },{i,nn}];
%  z2 = Table [{i+1},{i,nn}];
%  ax[x_]:= PadRight[x,2nn];
%  xa[x_]:= Insert[x,0,z2];
%  bx[x_]:= PadLeft [x,2nn];
%  xb[x_]:= Insert[x,0,z1];
%  ut[x_]:= Table[ ax[x[[i]]]-xa[x[[i]]],{i,nn}]~Join~
%           Table[ bx[x[[i]]]-xb[x[[i]]],{i,nn}];
%  t    = ut[t];
%  If[s==r,
%    cc  = Tuples[c[u,v],s];
%    tt   = Transpose[t];
%    nstt = NullSpace [tt];
%    o = cc . nstt[[13]];
%    Print["o=",o];
%    rank = MatrixRank[M[s]];
%    Print["s=",s," max t=",Max[Flatten[t]]];
%    ap[nstt]; 
% ];
%];
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Our observations and conjectures about BCH series}

Based on results given in Tab.~\ref{tabNcomm}
we formulate the following observations and conjectures:
\begin{itemize}

\item
Observation 1:\\
If the order~$r$ is odd and between 1 and 11
then the number~$N_w(r)$ of words of order~$r$
in the BCH series as a function of~$r$
forms a strictly
increasing sequence 2, 6, 30, 126, 390, 2046.

\item
Conjecture 1:\\
If the order~$r$ is odd 
then the number~$N_w(r)$ of words of order~$r$
in the BCH series as a function of~$r$
forms a strictly
increasing sequence.

\item
Observation 2:\\
If the order~$r$ is even and between 4 and 12
then the number~$N_w(r)$ of words of order~$r$
in the BCH series is 
$$
 N_w(r)= N_w(r-1).
$$

\item
Conjecture 2:\\
If the order~$r$ is even and more than 2
then the number~$N_w(r)$ of words of order~$r$
in the BCH series is 
$$
 N_w(r)= N_w(r-1).
$$

\item
Observation 3:\\
If the order~$r$ is between 2 and 12
then the minimal number~$N_{com}(r)$ of commutators of order~$r$
is strictly less than
then the number~$N_w(r)$ of words of order~$r$
in the BCH series.

Thus it is more convenient to write the BCH series 
in commutators than in words, because then the results are shorter.
There is one even more important reason to write the terms
using commutators: a commutator of operators corresponding
to vector fields corresponds to a vector field as well.

\item
Conjecture 3:\\
If the order~$r$ is more than 1
then the minimal number~$N_{com}(r)$ of commutators of order~$r$
is strictly less than
then the number~$N_w(r)$ of words of order~$r$
in the BCH series.

\item
Observation 4:\\
For certain~$r$ (e.g. 2, 3, 5, 7, 11) 
the number $N_w$ of words reaches
a remarkably high value. Considering that two words, 
namely 
the word consisting of $\U$'s only and 
the word consisting of $\V$'s only
cannot be generated as a linear
combination of commutators, 
the value $N_w=N_{tot}-2$ is the highest possible result 
consistent with the theory.
%
%\item
%The results given in 
%Tab.~\ref{tabNcomm}
%allow to formulate a conjecture:
%For even~$r$ 
%(let us write it as $r=2k$) 
%the number $(N_{com})_{2k}$ of commutators of order~$r$
%is less by one than the number of commutators of order~$r-1$.
%$$
% (N_{com})_{2k} = (N_{com})_{2k-1} - 1.
%$$
%
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Sparsest form of BCH series}

Using (\ref{hs}) with exp given by (\ref{expop}) 
and log given by (\ref{logop}) we get the series to solve (\ref{bch})
in words. Then we can solve (\ref{Mba}) to get the series in
commutators.
As mentioned earlier, the system of linear algebraic equations (\ref{Mba})
is under-determined, it has more unknowns (the coefficients $(b_r)_k$)
than independent equations. It has infinitely many solutions.
We can search for the sparsest solution, i.e. the solution
that has the least possible number of nonzero elements.

The sparsest solution of an under-determined system
is not only elegant, it is also practical, because 
it contains the information in the most compressed way.
This is desirable when we want to store or to transmit the information.
This problem has been discussed intensively recently, see \cite{Donoho}
and the references there.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Finite forms of BCH series}

In this section we present two special cases when the BCH series
has only a finite number of terms.

If $\U$ and $\V$ commute, i.e. if 
$$
 [\U,\V] = 0
$$
then 
$$
 \hat W = \log (\exp \U \exp\V)
$$
has finitely many terms, namely
$$
 \hat W = \U + \V.
$$
This is also the case when $\U$ and $\V$ are replaced by (complex) numbers.

Another special case of finite form of the BCH series is this one:
If 
$$
 [\U,\V] \neq 0
$$
$$
 [\U,\U,\V] = 0
$$
$$
 [\V,\U,\V] = 0
$$
then 
$$
 \hat W = \U + \V + \frac{1}{2} [\U,\V].
$$
As a consequence we can write
$$
 \exp\U\exp\V=\exp(\U+\V)\exp(\frac{1}{2}[\U,\V])=\exp\V\exp\U\exp[\U,\V].
$$
Proof:
since $[\U,\V]$ commutes with both $\U$ and $\V$ 
and thus also with $(\U+\V)$
(even though $\U$ and $\V$ do not commute)
we can write
$$
 \exp\U\exp\V=\exp(\U+\V+\frac{1}{2}[\U,\V])=
 \exp(\U+\V)\exp(\frac{1}{2}[\U,\V]).
$$
This proves the first statement.
And we can also write
$$
 \exp\V\exp\U\exp[\U,\V]=
 \exp(\U+\V)\exp(-\frac{1}{2}[\U,\V])\exp[\U,\V]=
$$
$$
  = \exp(\U+\V)\exp(\frac{1}{2}[\U,\V]).
$$
This proves the second statement.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Disprove of Kolsrud conjecture}

Kolsrud in \cite{Kolsrud} computes the sparsest form of the BCH
series up to the order 9 by intuitive applications
of simplification rules.
He observes that for the order $r$ between 1 and 9 
if the order $r$ is a prime number then 
the minimal number of commutators to express the terms of the BCH series
of order $r$ is equal to the rank of the matrix $M_r$ (in our
symbols).
Unfortunately, we have found that for $r=11$ which is prime,
there is a solution having 181 terms (see Table \ref{tabNcomm}).
This number may not be the least possible value, but it is 
definitely less than the rank of the matrix $M_r$ which is
186 for $r=11$. This result disproves the Kolsrud conjecture.

In the next subsection we describe our systematic algorithm
to find the sparsest form of the BCH series.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Our algorithm to find the sparsest solution}

We use a systematic approach to find the sparsest solution
to (\ref{Mba}), when the right hand side $\vec a_r$ is already known.
First, we reduce the size of the problem.
For the order $r$ there are originally
$2^r$ equations for $2^r$ unknowns. 
Since 
$$
 [\U,\U] =0 
$$
any commutator ending in two $\U$'s is equal to zero
(and the same is true for $\V$) and we can put to zero each 
coefficient $(b_r)_k$ in front of a commutator ending in two 
$\U$'s or in two $\V$'s. 
In this way we reduce the number of unknowns to one half 
to get $2^r/2$ unknowns.

Further, since
$$
 [\U,\V] = - [\V,\U]
$$
any commutator ending in $\dots,\V,\U$ is equal to minus one times
the same commutator with the last two symbols exchanged.
Thus we can put to zero each coefficient $(b_r)_k$ in front of a commutator
ending in $\dots,\V,\U$. 
In this way we reduce the number of unknowns to one half once more
to get $2^r/4$ unknowns.

We can look at the system of linear algebraic equations
as of finding the sparsest linear combination of columns of the matrix
to get the right hand side vector.
Thus we can remove any column that is a multiple of another column
and we put the unknown corresponding to the removed column to zero.
Unfortunately we cannot remove a column that is a linear combination
of other columns because then the solution to the new system
might have worse sparsity (more nonzero components).
But we can remove any row (any equation) that is a linear
combination of other rows.

Performing these reductions along with the Gauss elimination
the system can be reduced considerably ending with $m$ rows (equations)
and $n$ unknowns.
After this reduction we try to find the sparsest solution.

First, we test whether there is a solution containing only
one nonzero component. To do this we remove all but one columns
of the matrix $M_r^T$. This can be done in $n$ ways.
We try to solve this system in all these $n$ cases.
If the rank of the matrix of coefficients is equal to the rank
of the augmented matrix 
(the matrix with the right hand side column added) 
then the system has solution 
(according to the Frobenius theorem), 
if the ranks do not agree,
the system has no solution.
If we find a solution we are done.

If not, we test whether there is a solution containing two
nonzero components. To do this we remove all but two columns
of the matrix $M_r^T$. This can be done in
$$
 \left(\begin{array}{c} n \\ 2 \end{array} \right)
$$
ways.
In all these cases we again test whether there is a solution 
using the Frobenius theorem. If so, we are done.

If not, we test whether there is a solution containing three
nonzero components. To do this we remove all but three columns
of the matrix $M_r^T$. This can be done in
$$
 \left(\begin{array}{c} n \\ 3 \end{array} \right)
$$
ways.
In all these cases we again test whether there is a solution 
using the Frobenius theorem. If so, we are done.
If not we test whether there is a solution containing 4, then 5,
then 6 and so on solutions. This algorithm stops after a finite 
(though possibly large)
number of steps.

As shown in Table \ref{tabNcomm} we have found the sparsest 
solution using the algorithm described above for order up to 10
and the minimal number $N_{com}$ of the commutators 
in the terms of order $r$ in the BCH series is given 
in the table.
For order 11 and 12 we have found a solution 
using the Gauss elimination 
without removing columns
but this solution
is not guaranteed to be the sparsest one 
because the algorithm described above takes too long for large~$r$.
So the numbers in Table \ref{tabNcomm}
give only an upper estimate for $N_{com}$ for $r$ 11 and 12.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Example 2: Zig+Zag system}

To illustrate the above results consider the following 1-dim example
of a zig+zag system with 
$$
 u(x) = x^2
$$
acting for time $p$
and
$$
 v(x) = x^3
$$
acting for another time $p$
with the initial condition 
$$
 x(0) = x_0 = 1.
$$
We first find $x_1$ and $x_2$ by separation of variables
and we find the Taylor series of $x_2$ as a function of the
parameter $p$. Then we compare this result with that 
obtained by the BCH series.

The ODE 
$$
 \dot x = x^2
$$
has the solution 
$$
 x(t) = \frac{x(0)}{1-t x(0)}
$$
thus
$$
 x_1 = \frac{1}{1-p}.
$$

The ODE 
$$
 \dot x = x^3
$$
has the solution 
$$
 x(t) = \frac{x(0)}{\sqrt{1-2t x(0)^2}}
$$
thus
$$
 x_2 = \frac{x_1}{\sqrt{1-2 p x_1^2}}=
 \frac{  \frac{1}{1-p}}   {\sqrt{1-\frac{2p}{(1-p)^2}}}=
 \frac{1}{\sqrt{1-4p+p^2}}
$$
for 
$$
 0 \le p < 2-\sqrt{3}.
$$
The Taylor series of $x_2$  is
$$
 x_2 =
 1+2 p + \frac{11 p^2}{2} + 17 p^3 + \frac{443 p^4}{8} +
 {\cal O}\left(p^{5}\right).
$$

Let us now use the BCH series to find $x_2$ and to compare it with 
the above result.
We first prepare the commutator 
$$
 [u,v] = v'\cdot u - u'\cdot v
$$
of two power functions $x^m$ and $x^k$
$$
 [x^m,x^k] = k x^{k-1} x^m - m x^{m-1} x^k = (k-m) x^{m+k-1}.
$$
Thus
$$
 [u,v] = [x^2,x^3] = x^4,
$$ 
$$
 [u,u,v] = [x^2,x^4] = 2 x^5,
$$
$$
 [v,u,v] = [x^3,x^4] = x^6, 
$$
$$
 [u,v,u,v] = [x^2,x^6] = 4 x^7
$$
etc.

The shortcut vector field $w$ 
that moves the same initial condition $x_0$ to the same 
final point $x_2$ 
(not necessarily, actually almost never, via the point $x_1$)
in the same time $t=2p$ 
can be found by the BCH formula, recalling (\ref{uU})-(\ref{wW})
$$
 w =  \frac{1}{2p} W =
$$
$$
 = 
 \frac{1}{2p}(U+V+\frac{1}{2}[U,V]+\frac{1}{12}([U,U,V]-[V,U,V])-
 \frac{1}{24}[U,V,U,V]+{\cal O}(p^5) ) =
$$
$$
 = 
 \frac{u+v}{2} +\frac{p}{4}[u,v]+\frac{p^2}{24}([u,u,v]-[v,u,v])-
 \frac{p^3}{48}[u,v,u,v]
+{\cal O}(p^4)=
$$
$$
=
\frac{x^2+x^3}{2}+\frac{p x^4}{4}+\frac{p^2}{24}(2x^5-x^6)-\frac{p^3x^7}{12}
+{\cal O}(p^4).
$$
Then the solution of the ODE
$$
 \dot x = w(x)
$$
with the initial condition 
$$
 x(0) = 1
$$
in time $t=2p$ can be found by Taylor expansion
$$
 x(t)=x(0)+t\dot{x}(0) + 
 \frac{t^2}{2} \ddot{x}(0) + 
 \frac{t^3}{3!} x^{(3)}(0) + 
 \frac{t^4}{4!} x^{(4)}(0) 
 +{\cal O}(t^5)
$$
$$
 x_2 = x(2p)=
 1+2 p + \frac{11 p^2}{2} + 17 p^3 + \frac{443 p^4}{8} +
 {\cal O}\left(p^{5}\right)
$$
which agrees completely with the result obtained without BCH.

% Note, however, that the solution 
% \begin{equation}
%  \label{x2ex}
%  x_2 = \frac{1}{\sqrt{ (\frac{1}{x_0}-p)^2 - 2p }}
% \end{equation}
% is not of the form
% $$
%  x_2 = \varphi(2p,x_0)
% $$
% with $\varphi$ being a flow.
% By flow $\varphi(t,x_0)$ we mean a diffeomorphism
%  $\varphi: R \times R^n \rightarrow R^n$
% % solution to ODE with the initial condition $x_0$ 
% satisfying
% \begin{equation}
%  \label{compflow}
%  \varphi(s, \varphi(t,x_0)) = \varphi(t+s,x_0).
% \end{equation}
% Assuming (\ref{x2ex}) is a flow satisfying (\ref{compflow}) we get
% $$
%  \frac{1}{ \sqrt{ ( \sqrt{(\frac{1}{x_0}-t)^2-2t}-s)^2-2s } }=
%  \frac{1}{ \sqrt{ ( \frac{1}{x_0}-(t+s))^2-2(t+s)}}
% $$
%  which can be reduced to
% $$
%  s \sqrt{(\frac{1}{x_0}-t)^2-2t} = s(\frac{1}{x_0} -t).
% $$
% This is satisfied only for $t=0$ or $s=0$.
% The consequence is that the shortcut~$w$ that we have found
% is only an approximation useful for small $p$, 
% not the exact replacement of the original dynamics.
% And this is true in a general case too.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Application to linear ODE}

Suppose the two vector fields $U$ and $V$ are linear, i.e.
there are two $n\times n$ matrices $\cal U$ and $\cal V$ such that
$$
 U(x) = {\cal U} \cdot x 
$$
and
$$
 V(x) = {\cal V} \cdot x .
$$

Note to symbols: names of matrices are usually ordinary capital 
letter. We want to distinguish between a vector field and
a matrix and still use related symbols, so we choose
calligraphic capital letters for matrices.

Then the first order term in the expansion of the operator~$\hat W$
$$
 \hat W_1 = \U + \V
$$
gives the first order term in the expansion of the vector field~$W$
$$
 W_1 = \hat W_1 I = (\U + \V) I = U + V.
$$
The second order term in the expansion of the operator
$$
 \hat W_2 = \frac{1}{2}[\U,\V]
$$
gives
$$
 W_2 = \hat W_2 I  = 
 \frac{1}{2}[\U,\V]I = 
 \frac{1}{2}(\U\V-\V\U)I = 
 \frac{1}{2}(V'\cdot U-U'\cdot V).
$$
Since 
$$
 U'(x) = \cal U 
$$
and
$$
 V'(x) = \cal V
$$
we have
$$
 W_2(x) = 
 \frac{1}{2}(
  {\cal V} \cdot {\cal U} \cdot x-
  {\cal U} \cdot {\cal V} \cdot x)
$$
so we can write
$$
 W_2 = 
 \frac{1}{2}(
  {\cal V} \cdot {\cal U} -
  {\cal U} \cdot {\cal V} )
 =
 \frac{1}{2}[{\cal V},{\cal U}] =
 -\frac{1}{2}[{\cal U},{\cal V}].
$$
As a result, for linear systems
we can use the BCH  series directly for matrices,
but then we have to multiply by $-1$ each even order term.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Ambiguity of matrix commutator}

This uncomfortable result leads various authors to define 
the commutator of two matrices in different ways.
While Rossmann in \cite{Rossmann} p.14 defines the commutator
of two matrices $\cal U$ and $\cal V$
$$
 [{\cal U},{\cal V}]_{Rossmann}={\cal U}\cdot{\cal V}-{\cal V}\cdot{\cal U},
$$
Olver in \cite{Olver} p.44 defines the same commutator as
$$
 [{\cal U},{\cal V}]_{Olver}={\cal V}\cdot{\cal U}-{\cal U}\cdot{\cal V}.
$$
%This is a dangerous situation.
We follow Rossmann, 
because this notation is more common in literature.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Zig+Zag+Zug Systems}
 \label{sec.zzz}

In the preceding chapters we studied dynamical systems
generated by two vector fields $u,v:R^n\rightarrow R^n$
acting one after the other in time.
In this chapter we investigate a generalization to the case
when there are more than two vector fields. As 
the first step in this generalization is the case with three
vector fields, we choose to call them ``zig+zag+zug" systems.
To be more precise we study dynamical systems of the form

\begin{equation}
 \frac{d x}{d t}=f(x,t)
\end{equation}
where $f(x,t)$ is $3p$-periodic in~$t$ satisfying
\begin{equation}
 f(x,t)
 =\left\{  % } aby byly v paru
       \begin{array}{r c l}
        u_1(x)  & \mbox{if} \quad  & 0 \le t < \ p \\
        u_2(x)  & \mbox{if} \quad  & p \le t < \ 2p \\
        u_3(x)  & \mbox{if} \quad  & 2p \le t < 3p.
       \end{array}
       \right.
\end{equation}
In analogy with (\ref{uU}) we introduce rescaled vector fields
$$
 U_k = p\ u_k \qquad \mbox{ for } k=1,\dots,3
$$
and the corresponding operators $\hat U_k$.

We again want to find a shortcut -- a single vector field $W$
that moves the point in the state space from the same initial
condition $x_0=x(0)$ to the same final point $x_3=x(3p)$.

To use our previous results it is convenient to introduce 
a new symbol $\vdash\!\!\!\dashv$ resembling the letter H for
Hausdorff (see \cite{Bourbaki}) by 
$$
 \U\vdash\!\!\!\dashv\V = \log(\exp\U\exp\V)
$$
or equivalently
\begin{equation}
 \label{eH}
 \exp (\U\vdash\!\!\!\dashv\V) = \exp \U \exp \V.
\end{equation}
We want to find $\hat W$ satisfying 
$$
 \exp \hat W = \exp \hat U_1 \exp \hat U_2 \exp \hat U_3.
$$
Using (\ref{eH}) twice we get
$$
 \exp \hat W = 
 \exp (\hat U_1 \vdash\!\!\!\dashv \hat U_2) \exp \hat U_3 = 
 \exp (\hat U_1 \vdash\!\!\!\dashv \hat U_2 \vdash\!\!\!\dashv \hat U_3 )
$$
and thus
$$
 \hat W = \hat U_1 \vdash\!\!\!\dashv \hat U_2 \vdash\!\!\!\dashv \hat U_3.
$$
It is easy to show that our previous result 
$$
 \U\vdash\!\!\!\dashv\V = \U + \V + \frac{1}{2}[\U,\V] 
 + {\cal O}(p^3)
$$
implies
$$
 \hat U_1 \vdash\!\!\!\dashv \hat U_2 \vdash\!\!\!\dashv \hat U_3  = 
 \hat U_1 + \hat U_2 + \hat U_3 
 + \frac{1}{2} (
 [\hat U_1,\hat U_2] +
 [\hat U_1,\hat U_3] +
 [\hat U_2,\hat U_3])
 + {\cal O}(p^3).
$$
It is easy to generalize this result from 3 to $m$ vector fields 
$ \hat U_1,\dots,\hat U_m$ acting one after another. 
By induction we find 
$$
 \hat W = \sum_{k=1}^m \hat U_k  + 
 \frac{1}{2} 
 \sum_{j=1}^{m-1} \sum_{k=j+1}^m 
 [\hat U_j,\hat U_k]
 + {\cal O}(p^3).
$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Zig+Zag-Zig-Zag System}
  \label{sec.zzzzs}

Consider the dynamical system generated by four vector fields
$u,v,-u,-v$ acting one after the other. Using the result of 
the previous chapter for $m=4$ we have
$$
 \hat W = \sum_{k=1}^4 \hat U_k  + 
 \frac{1}{2} 
 \sum_{j=1}^{3} \sum_{k=j+1}^4 
 [\hat U_j,\hat U_k]
 + {\cal O}(p^3) =
$$
$$
 = \hat U_1 + 
   \hat U_2 + 
   \hat U_3 + 
   \hat U_4 + 
  \frac{1}{2} ( 
   [\hat U_1,\hat U_2] +
   [\hat U_1,\hat U_3] +
$$
$$ 
  +[\hat U_1,\hat U_4] +
   [\hat U_2,\hat U_3] +
   [\hat U_2,\hat U_4] +
   [\hat U_3,\hat U_4])
 + {\cal O}(p^3).
$$
As
\begin{eqnarray*}
 \hat U_1 &=&  \U \\
 \hat U_2 &=&  \V \\
 \hat U_3 &=& -\U \\
 \hat U_4 &=& -\V
\end{eqnarray*}
we have
$$
 \hat W = \U + \V - \U - \V +
  [ \U, \V] + 
  [ \U,-\U] + 
  [ \U,-\V] + 
$$
$$ 
 +[ \V,-\U] + 
  [ \V,-\V] + 
  [-\U,-\V] 
 + {\cal O}(p^3)
= [\U,\V]
 + {\cal O}(p^3).
$$
This gives an important interpretation of the commutator of two
vector fields: it is a vector field which is 
tangent to the curve with parametric equation
$$
 x_4 = x_4 (p)
$$ 
of the ``zig+zag-zig-zag" system generated by four vector fields,
where the third vector field is the opposite to the first one and 
the fourth vector field is the opposite to the second one.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\appendix

\section{Appendices}

\subsection{Proof of the BCH formula}
 \label{ApDjokovic}

The main idea of the proof by Djokovic in \cite{Djokovic}
that $\hat W$ can be expressed by commutators of $\U$ and $\V$ only
is as follows.
He denotes 
$$
 h = \log (\exp(t\U) \exp(t\V)) = 
  \sum_{n=1}^\infty 
   h_n(\U,\V)\ t^n
$$
and
$$
 \ad a\ (x) = [a,x] = a x - x a 
$$
and $D$ the formal differentiation with respect to~$t$.
Then he derives step by step 
% by simple manipulation
two different results of 
$(D\exp(h))\exp(-h)$, namely
$$
 (D\exp(h))\exp(-h)=
 \U+\exp(\ad h)(\V)
$$
and
$$
 (D\exp(h))\exp(-h)=
 f(\ad h) (D h)
$$
where 
$$
 f(x) = 
 \frac{\exp(x)-1}{x}=
 \sum_{n=1}^\infty \frac{x^{n-1}}{n!}=
 1+\frac{x}{2}+\frac{x^2}{6}+\frac{x^3}{24}+\frac{x^4}{120}+\dots
$$
Thus 
$$
 f(\ad h) (D h) = \U+\exp(\ad h)(\V).
$$
Introducing 
$$
 g(x) = \frac{1}{f(x)}=
 \frac{x}{\exp(x)-1}=
 1-\frac{x}{2}+\frac{x^2}{12}-\frac{x^4}{6! }+\frac{x^6}{6\cdot 7!}+\dots
$$
$D h$ can be expressed as
$$
 D h = g(\ad h)(\U+\exp(\ad h)(\V)).
$$
This is a differential equation for~$h$ that 
provides a recursion formulas for~$h_n$.
The coefficient of~$t^n$ on the left is 
$$
 (n+1) h_{n+1}
$$
and that on the right is a linear combination of expressions of the form 
$$
 \ad h_{m_1} \dots \ad h_{m_k} (z)
$$
where $z=\U$ or $\V$ and $m_1 + \dots +m_k=n$.
Thus if $h_j$ are Lie elements (i.e. can be expressed by commutators
of $\U$ and $\V$ only) for $j\le n$ so is~$h_{n+1}$.
Since $h_0 = 0$ (and also $h_1 = \U + \V$), $h$ is a Lie element.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Alternative way to find the terms in the BCH series}
 \label{ApReinsch}

Reinsch in \cite{Reinsch} in order to find the $r$-th order term
in the BCH series uses two $(r+1)\times(r+1)$ matrices
$$
 M_{i j} = \delta_{i+1,j}
$$
and 
$$
 N_{i j} = \delta_{i+1,j}\ \sigma_i
$$
i.e.
$$ 
 M=
 \left(
  \begin{array}{c c c c c c c}
    0 & 1 & 0 & \dots & & &  \\ 
      & 0 & 1 &  0 & \dots & &   \\ 
      &   & . &    &       & &   \\ 
      &   &   & .  &       & &   \\ 
      &   &   &    & .     & &   \\ 
      &   &   &    &       & 0 & 1 \\ 
      &   &   &    &       &   & 0 
  \end{array}
 \right)
$$ 
and
$$ 
 N=
 \left(
  \begin{array}{c c c c c c c}
    0 & \sigma_1 & 0 & \dots & & &  \\ 
      & 0 & \sigma_2 &  0 & \dots & &   \\ 
      &   & . &    &       & &   \\ 
      &   &   & .  &       & &   \\ 
      &   &   &    & .     & &   \\ 
      &   &   &    &       & 0 & \sigma_r \\ 
      &   &   &    &       &   & 0 
  \end{array}
 \right).
$$ 
When these matrices multiply a $(r+1)$-dim vector 
$$
 (0,0,\dots,1)^T
$$
from the left, then the only nonzero component of the vector
moves one position up with each multiplication.
Multiplication by the matrix $M$ leaves this nonzero component
unchanged, while multiplication by the matrix $N$ multiplies
this nonzero component by one more $\sigma_i$, thus leaving 
a mark to be used later by the transformation rule $T$.

Then the $r$-th order term of 
$$
 z=\log(\exp(M)\exp(N))
$$ is
$$
 z_r=T(\log(\exp(M)\exp(N)))_{1,r+1}
$$
where $T$ stands for a transformation that puts 
$y$ at $i$-th position when $\sigma_i$ is found
(this is the mark left by multiplication by the matrix $N$)
and it puts $x$ at $i$-th position when $\sigma_i$ is not found.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection {Example 3: Parking problem}
% Zig+Zag-Zig-Zag system 

 \label{Apcar}

To illustrate the Zig+Zag-Zig-Zag system consider 
the following example from \cite{Nelson} and \cite{Rossmann}:
how to park your car into a space 
that is just a little longer than your car?

A car can be driven forward or backward (along an arc 
depending on the position of the driving wheel).
In a situation shown in Fig. \ref{Fparking}
where each rectangle represents a car parked in a street,
the driver may wish to be able to drive his car to the right instead.
Surprising enough, this can be done as a limiting case of a Zig+Zag-Zig-Zag
dynamical system.

\begin{figure}[htb]
\begin{center}
\setlength{\unitlength}{3.0pt}
\begin{picture}(104,44)
 \put(0,0){\line(1,0){104}}
 \put(0,44){\line(1,0){104}}
 \put(10, 2){\framebox(24,10){}}
 {
 \color{red}
 \put(40,12){\framebox(24,10){}}
 }
 \put(70, 2){\framebox(24,10){}}
 \put(10,32){\framebox(24,10){}}
 \put(40,32){\framebox(24,10){}}
 \put(70,32){\framebox(24,10){}}
 \put(52,17){\vector(0,-1){10}}
 \put(55, 5){\large ?}
\end{picture}
\end{center}
 \caption{
  Parking problem: how to park your car into a space 
  that is just a little longer than your car?
  Each black rectangle represents a parked car in a street,
  the red rectangle is your car.
}
 \label{Fparking}
\end{figure}

\begin{figure}[htb]
\begin{center}
\setlength{\unitlength}{0.6pt}
\begin{picture}(300,300)
 {
 \color{red}
 \qbezier(248.589,166.742)(239.229,221.855)(229.87,276.967)
 }
 \qbezier(150,150)(199.294,158.371)(248.589,166.742)
 \qbezier(150,150)(189.935,213.484)(229.87,276.967)
 \qbezier(229.87,276.967)(198.128,296.935)(166.386,316.903)
 \qbezier(229.87,276.967)(225.191,304.524)(220.511,332.08)
 \qbezier(248.589,166.742)(268.589,166.742)(288.589,166.742)
 \qbezier(300.,150.)(300.,198.738)(271.353,238.168)
 \qbezier(271.353,238.168)(242.705,277.598)(196.353,292.658)
 \qbezier(196.353,292.658)(150.,307.719)(103.647,292.658)
 \qbezier(103.647,292.658)(57.2949,277.598)(28.6475,238.168)
 \qbezier(28.6475,238.168)(0,198.738)(0,150.)
 \qbezier(0,150.)(0,101.262)(28.6475,61.8322)
 \qbezier(28.6475,61.8322)(57.2949,22.4024)(103.647,7.34152)
 \qbezier(103.647,7.34152)(150.,-7.71933)(196.353,7.34152)
 \qbezier(196.353,7.34152)(242.705,22.4024)(271.353,61.8322)
 \qbezier(271.353,61.8322)(300.,101.262)(300.,150.)
 \qbezier(250.,150.)(250.,182.492)(230.902,208.779)
 \qbezier(230.902,208.779)(211.803,235.065)(180.902,245.106)
 \qbezier(180.902,245.106)(150.,255.146)(119.098,245.106)
 \qbezier(119.098,245.106)(88.1966,235.065)(69.0983,208.779)
 \qbezier(69.0983,208.779)(50.,182.492)(50.,150.)
 \qbezier(50.,150.)(50.,117.508)(69.0983,91.2215)
 \qbezier(69.0983,91.2215)(88.1966,64.9349)(119.098,54.8943)
 \qbezier(119.098,54.8943)(150.,44.8538)(180.902,54.8943)
 \qbezier(180.902,54.8943)(211.803,64.9349)(230.902,91.2215)
 \qbezier(230.902,91.2215)(250.,117.508)(250.,150.)
 \put(164.576,159.733){$\beta$}
 \put(210.972,294.763){$\beta$}
 \put(231.866,145.758){$A$}
 \put(257.212,175.691){$\alpha$}
 \put(236.26,287.125){$B$}
 \put(247.379,221.651){$L$}
 \put(197.422,169.394){$R_1$}
 \put(159.61,212.87){$R_2$}
\end{picture}
\end{center}
 \caption{
  A driving car with a fixed orientation of the driving wheel.
  $A$ (the center of the rear axle) moves on a circle with radius $R_1$,
  $B$ (the center of the front axle) moves on a circle with radius $R_2$,
  $L$ is the length of the car,
  $\alpha$ is the orientation of the car (relative to a chosen direction),
  $\beta$ is the orientation of the driving wheel of the car 
   (relative to the axis of the car).
}
 \label{Fcar}
\end{figure}

The configuration of a car in a plane can be given by four real values
(see Fig.\,\ref{Fcar}):
\begin{itemize}
\item
 the position of the center $B$ of the front axle in Cartesian
 coordinates $(x,y)$,
\item
 the angle $\alpha$ specifying the orientation 
 of the car relative to a chosen direction
\item
 the angle $\beta$ specifying the orientation 
 of the driving wheel of the car relative to the axis of the car.
\end{itemize}
 
The driver can perform two actions: steering (turning the driving wheel)
and driving. Let us denote the vector field describing steering
by Steer and the vector field describing driving by Drive.
Then 
\begin{equation}
 \label{Steer}
 \mbox{Steer} = (0,0,0,1)^T
\end{equation}
meaning 
\begin{eqnarray}
 \dot x &=& 0 \nonumber \\
 \dot y &=& 0 \nonumber \\
 \dot \alpha &=& 0 \nonumber \\
 \dot \beta  &=& 1.\nonumber
\end{eqnarray}
Using Fig.\,\ref{Fcar} we find that the center $B$ of the front axle
has coordinates
$$ x = R_2 \cos (\alpha+\beta-\frac{\pi}{2}) = R_2 \sin (\alpha+\beta) $$
$$ y = R_2 \sin (\alpha+\beta-\frac{\pi}{2}) =-R_2 \cos (\alpha+\beta) $$
and thus
$$ \dot x = \dot \alpha R_2 \cos (\alpha+\beta) = v \cos (\alpha+\beta) $$
$$ \dot y = \dot \alpha R_2 \sin (\alpha+\beta) = v \sin (\alpha+\beta) $$
where the velocity $v$ is 
$$
 v = \dot \alpha  R_2 .
$$
Since
$$
 \sin \beta = \frac{L}{R_2}
$$
we have
$$
 \dot \alpha = \frac{v}{R_2} = \frac{v}{L}\sin \beta.
$$
Choosing the length scale so that the length of the car is $L=1$
and choosing the time scale so that the velocity is $v=1$ we get
the vector field corresponding to driving
\begin{equation}
 \label{Drive}
 \mbox{Drive} = (
  \cos(\alpha+\beta), 
  \sin(\alpha+\beta),
  \sin\beta ,
  0)^T
\end{equation}
meaning 
\begin{eqnarray}
 \dot x      &=& \cos(\alpha+\beta) \nonumber \\
 \dot y      &=& \sin(\alpha+\beta) \nonumber \\
 \dot \alpha &=& \sin(\beta) \nonumber \\
 \dot \beta  &=& 0\nonumber .
\end{eqnarray}
Let us compute the commutator
$$
 [u,v] = v'\cdot u - u' \cdot v
$$
of the two vector fields Steer and Drive.
As
$$
 \mbox{Steer}' = 
\left( \begin{array}{l l l l}
                  0 & 0 & 0 & 0 \\
                  0 & 0 & 0 & 0 \\
                  0 & 0 & 0 & 0 \\
                  0 & 0 & 0 & 0
                 \end{array}
                 \right)
$$
and 
$$
 \mbox{Drive}' = 
\left( \begin{array}{c c c c}
                  0 & 0 & -\sin (\alpha+\beta) & -\sin (\alpha+\beta) \\
                  0 & 0 & \cos (\alpha+\beta) & \cos (\alpha+\beta) \\
                  0 & 0 & 0 & \cos (\beta) \\
                  0 & 0 & 0 & 0
                 \end{array}
                 \right)
$$
we have
\begin{equation}
 \label{wriggle}
 [\mbox{Steer},\mbox{Drive}] = 
(-\sin(\alpha+\beta), \cos(\alpha+\beta), \cos(\beta),0)^T
\end{equation}
meaning 
\begin{eqnarray}
 \dot x      &=& -\sin (\alpha+\beta) \nonumber \\
 \dot y      &=& \cos (\alpha+\beta)  \nonumber \\
 \dot \alpha &=& \cos (\beta) \nonumber \\
 \dot \beta  &=& 0\nonumber .
\end{eqnarray}
Note that the angle $\alpha$ gives the orientation of the car
(relative to a chosen direction),
the angle $\beta$ gives the orientation of the driving wheel
(relative to the car)
and the angle $\alpha+\beta$ gives the direction at which 
the center of the front axle moves when driving.
Then the vector field (\ref{wriggle}) describes a motion of the car
when the car moves in a direction perpendicular to the direction
it would drive; rotating with a constant angular velocity 
 $\dot \alpha = \cos (\beta)$; with fixed position $\beta$ of the driving 
wheel. This motion can be called Wriggle (\v{c}esky: vrt\v{e}t se) defined by
$$
 \mbox{Wriggle} = [\mbox{Steer},\mbox{Drive}].
$$
It is easy to show that 
$$
 [ \mbox{Wriggle} , \mbox{Steer}] = \mbox{Drive}
$$
thus giving nothing new. But 
$$
 [ \mbox{Wriggle} , \mbox{Drive}] = 
(-\sin (\alpha),\cos (\alpha),0,0)^T.
$$
This motion can be called Slide, thus
$$
 \mbox{Slide} = [ \mbox{Wriggle} , \mbox{Drive}]
$$
because the car moves in a direction perpendicular to its axis
and does not rotate! This is the very motion needed
to park a car to a space that is just a little
longer than the length of the car 
(or to come out of such a difficult position).

To give the complete picture we add that 
$$
 [\mbox{Slide},\mbox{Steer}] = (0,0,0,0)^T,
$$
$$
 [\mbox{Slide},\mbox{Drive}] = \sin\beta\ (\cos\alpha,\sin\alpha,0,0)^T,
$$
and
$$
 [\mbox{Slide},\mbox{Wriggle}] = \cos\beta\ (\cos\alpha,\sin\alpha,0,0)^T.
$$
The last two vector fields are proportional to the vector field
$$
 (\cos\alpha,\sin\alpha,0,0)^T
$$
which, in turn, is equal to Drive for $\beta=0$, thus giving nothing new.

This can be conveniently computed by the following 
{\it Mathematica} program
\begin{verbatim}
Steer = {0,0,0,1};
Drive = {Cos[a+b],Sin[a+b],Sin[b],0};
d[f_]:= Outer[D,f,{x,y,a,b}];
com[f_,g_] := d[g].f - d[f].g;
Wriggle = com[Steer,Drive];
Slide   = com[Wriggle,Drive]//Simplify;
SlideSteer = com[Slide,Steer];
SlideDrive = com[Slide,Drive];
SlideWriggle = com[Slide,Wriggle];

Print["Slide=",Slide];
Print["[Slide,Steer]=",SlideSteer];
Print["[Slide,Drive]=",SlideDrive];
Print["[Slide,Wriggle]=",SlideWriggle];
\end{verbatim}

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Curriculum Vitae}

Born: 26.12.1962 in Pardubice, Czech Republic.\\
Father: Doc. MUDr. \v{C}estm\'{\i}r Pokorn\'y, CSc.\\
Mother: MUDr. Irena Pokorn\'a.  \\
Married, two children.  \\[0.8em]
1986: RNDr. in Physics, 
      Faculty of Mathematics and Physics,
      Charles University Prague. \\
1986 - 1987: Military service. \\
1987 - 1994: Dept.\,of Chemical Engineering, 
Institute of Chemical Technology, Prague.\\
1.7.-30.9.1991 Dept.\,of Physical Chemistry, 
Free University of Brussels, Belgium. \\
1994 - today: Dept.\,of Mathematics, 
Institute of Chemical Technology, Prague.\\[0.8em]
Teaching experience:\\
since 1994: Calculus I and Calculus II in Czech.\\
since 1996: Calculus I and Calculus II in English.\\
since 1997: Fourier Transform in Czech.
\\[0.7em]
Language experience:\\
state language examination in English, German and Russian,
partially French.\\
Member of Mensa.

\begin{center}
 \includegraphics[width=\figsize]{rod.ps}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \section {{\it Mathematica} program}
% \label{AppM}
%
% Here is the {\it Mathematica} code to generate the 
% coefficients in the BCH series in commutators.
%
% \verbatiminput{bch12ff.m}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \section {C program for words}
% \label{AppCbch}
%
% Here is the C code to generate the coefficients in the BCH series in words.
%
% \verbatiminput{bch.c}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \section {C program for the sparsest solution}
% \label{AppCsparsest}
%
% Here is the C code to find the sparsest form of 
% the terms in the BCH series in commutators.
%
% \verbatiminput{sparsest.c}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage

\begin{thebibliography}{99}

\bibitem{Baker}
Baker H.\,F.: Alternants and continuous groups. 
Proc. London Math. Soc., II. Ser. {\bf 3}, 24-27 (1905).

\bibitem{Bourbaki}
Bourbaki N.: {\it Lie Groups and Lie Algebras.} 
Springer (1989).

\bibitem{Bunim}
Bunimovich G.\,A. et.al.:
SO2 oxidation in reverse flow reactor: influence of a vanadium catalytic
dynamic properties. 
Chem. Eng. Sci. {\bf 50}, 565 (1995).

\bibitem{Cartier}
Cartier P.: Demonstration algebrique 
de la formule de Hausdorff.
Bull. Soc. math. France {\bf 84}, 241-249 (1956).

\bibitem{Djokovic}
Djokovic D.\,Z.:
An Elementary Proof of the 
Baker--Campbell--Hausdorff--Dynkin Formula.
Math. Z. {\bf 143}, 209-211 (1975).

\bibitem{Donoho}
Donoho D.\,L.: 
For Most Large Under-determined Systems of Linear Equations 
the Minimal $\ell^1$-norm Solution is also the Sparsest Solution.
Preprint {\verb|http://www-stat.stanford.edu/~donoho|}

\bibitem{Dynkin}
Dynkin E.\,B.:
On the representation of the series log($e^x e^y$) with 
non-commuting $x$ and $y$ by commutators. 
Mat. Sbornik, n. Ser. {\bf 25}, 155-162 (1949).

\bibitem{Filippov}
Filippov A.\,F.: 
{ \it Differential Equations with Discontinuous Righthand Sides.}
Kluwer Academic Publishers Dordrecht (1988).

\bibitem{Hochschild}
Hochschild G.: {\it The structure of Lie groups.}, 
San Francisco--London--Amsterdam: Holden--Day, 1965.

\bibitem{Howe}
Howe R.: Very basic Lie theory. Amer. Math. Monthly {\bf 90} (1983).

\bibitem{Israel}
Israel R.: discussion in the {\tt sci.math.research} newsgroup.
% Department of Mathematics, University of British Columbia, 
% Vancouver, BC, Canada,

\bibitem{Khinast}
Khinast J., Luss D.:
Efficient bifurcation analysis of periodically--forced
distributed parameters systems.
Comp. Chem. Eng. {\bf 24}, 139 (2000).

\bibitem{Kolsrud}
Kolsrud M.: Maximal reductions in the Baker--Hausdorff formula.
J. Math. Phys. {\bf 34}, 270 (1993).

\bibitem{Lorenz}
Lorenz E.\,N.: Deterministic Nonperiodic Flow.
J. Atmos. Sci., Vol. {\bf 20}, (1963).

\bibitem{Matros85}
Matros Y.\,Sh.:
{\it Unsteady processes in catalytic reactor.}
Elsevier, Amsterdam (1985).

\bibitem{Matros89}
Matros Y.\,Sh.:
{\it Catalytic process under unsteady--state conditions.}
Elsevier, Amsterdam (1989).

\bibitem{Nelson}
Nelson E.:
{\it Tensor Analysis.}
Princeton University Press, Princeton (1967).

\bibitem{Nieken}
Nieken U., Kolios G., Eigenberger G.:
Control of ignited steady state in autothermal fixed bed
reactors for catalytic combustion.
Chem. Eng. Sci. {\bf 49}, 5507 (1994).

\bibitem{Olver}
Olver P.\,J.:
{\it Applications of Lie Groups to Differential Equations.}
Springer (1993).

\bibitem{Oteo}
Oteo J.\,A.: 
The Baker Campbell Hausdorff formula and nested commutator identities.
J. Math. Phys. {\bf 32}, 419 (1991).

\bibitem{Ottino89}
Ottino J.\,M.:
 {\it The kinematics of mixing, stretching, 
chaos and transportations}.
  Cambridge University Press (1989).

\bibitem{Pejas}
Pejas W.: 
Ein Beweis der qualitativen Aussage 
der Campbell--Hausdorff Formel
fur analytische Gruppen. Arch. der Math. {\bf 19}, 453-456 (1968).

\bibitem{Rehacek1}
\v{R}eh\'{a}\v{c}ek J., Kub\'{\i}\v{c}ek M., Marek M.:
Modeling of a tubular catalytic reactor with flow reversal.
Chem. Eng. Sci. {\bf 47}, 2897 (1992).

\bibitem{Rehacek2}
\v{R}eh\'{a}\v{c}ek J., Kub\'{\i}\v{c}ek M., Marek M.:
Periodic, quasi--periodic and chaotic spatio--temporal patterns
in a tubular catalytic reactor with periodic flow reversal.
Comp. Chem. Eng. {\bf 22}, 283 (1998).

\bibitem{Reinsch}
Reinsch M.\,W.: A simple expression for the terms in the 
Baker--Campbell--Hausdorff series.
J. Math. Phys. {\bf 41}, 2434 (2000).
%{\tt arXiv:math-ph/9905012 v2 13 Jan 2000}

\bibitem{Rossmann}
Rossmann W.:
{\it Lie Groups An Introduction Through Linear Groups.}
Oxford University Press (2002).

\bibitem{RussoC}
Russo L. et.al.:
Symmetry properties and bifurcational analysis of a class
of periodically forced chemical reactors. 
Chem. Eng. Sci. (2002).

\bibitem{RussoI}
Russo L. et.al.:
Nonlinear Analysis of a network of three continuous stirred tank
reactors with periodic feed switching: symmetry and symmetry breaking.
Int.\,J.\,Bif.\,Chaos {\bf 14}, 1325 (2004).

\bibitem{Snyder}
Snyder J.\,D., Subramanian S.:
A novel reverse flow strategy for ethylbenzene
dehydrogenation in a packed bed reactor.
Chem. Eng. Sci. {\bf 49}, 5585 (1994).

\bibitem{Varadarajan}
Varadarajan V.\,S.: {\it Lie groups, 
Lie algebras, and their representations.}
Englewood Cliffs; Prentice Hall 1974.

\bibitem{Xiao}
Xiao W.\,D., Wang H., Yuan W.\,K.:
Unsteady state converters with inter--stage heat removal.
Chem. Eng. Sci. {\bf 54}, 4629 (1999).

\bibitem{Zufle}
Zufle H., Turek T.:
Catalytic combustion in a reactor with periodic flow reversal.
Chem. Eng. Proc. {\bf 36}, 327 (1997).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{\bf Own publications:}
\end{center}

\bibitem{KP96}
Kl\'{\i}\v{c} A., Pokorn\'y P.:
On Dynamical Systems Generated 
by Two Alternating Vector Fields.
Int.\,J.\,Bif.\,Chaos 
{\bf 6}, 2015 (1996).

\bibitem{Maslo2002}
Kl\'{\i}\v{c} A., Pokorn\'y P., \v{R}eh\'{a}\v{c}ek J.:
Zig-Zag Dynamical Systems and the 
Baker-Campbell-Hausdorff Formula.
Math. Slovaca {\bf 52}, 79-97 (2002).

\bibitem{w1}
 Pokorn\'y P., Marek M.:
     Origin and Dimension of Strange Non-chaotic Attractors.
     Series on Nonlinear Science E, Vol. 4 (Eds. Kapitaniak T. and Brindley J.)
     World Scientific Singapore (1994).

\bibitem{w2}
 Pokorn\'y P., Schreiber I., Marek M.
     On the Route to Strangeness without Chaos in the Quasi-periodically Forced 
     van der Pol Oscillator.
     Chaos, Solitons \& Fractals Vol. 7 No. 3 (1996).

\bibitem{w3}
 Dynamical Systems Generated by Two Alternating Vector Fields 
     and the Baker-Campbell-Hausdorff Formula. In Proceedings
     Int. Symposium of Nonlinear Theory and its Applications,
     Nolta 2000, 17.-21.9.2000 Dresden, Germany. IEICE. Eds. S. Oishi and 
     W. Schwarz. ISBN 3-933592-84-4, pp.827-830.

\bibitem{w4}
 Excitable dynamical systems. In 12th Czech-Slovak-Polish Optical
     Conference on Wave and Quantum Aspects of Contemporary Optics.
     Pe\v{r}ina, Hrabovsk\'{y}, K\v{r}epelka, Editors. 
     Proceedings of SPIE Vol. 4356, 
     4356-60 (2001). ISBN 0-8194-4047-7.

\bibitem{w5}
 Turz\'{\i}k, Kl\'{\i}\v{c}, Dubcov\'a, Pokorn\'y:
     Stability of Spatially Homogeneous Solutions in Lattice Dynamical 
     Systems.  Proceedings of 3rd Scientific Colloquium, Math. Dept. 
     Prague Institute of Chemical Technology, 2001, pp. 38-46.

\bibitem{w6}
 Kl\'{\i}\v{c}, Turz\'{\i}k, Dubcov\'a, Pokorn\'y:
     A Simple Example of Lattice Dynamical System with Spatial Chaos.
     Proceedings of 3rd Scientific Colloquium, Math. Dept. 
     Prague Institute of Chemical Technology, 2001, pp. 47-54.

\bibitem{w7}
 Pokorn\'y, Turz\'{\i}k, Dubcov\'a, Kl\'{\i}\v{c}: 
     On Spectra of Certain Class of Linear Operators.
     Proceedings of 3rd Scientific Colloquium, Math. Dept. 
     Prague Institute of Chemical Technology, 2001, pp. 55-63.

\bibitem{w8}
 Turz\'{\i}k, Dubcov\'a, Kl\'{\i}\v{c}, Pokorn\'y:
     Application of Bernoulli Shift to Investigation of Stability of
     Spatially Periodic Solutions in Lattice Dynamical Systems.
     Dynamical Systems: An Int. Journal Vol.18, No.1 (2003) 23-33. 
     ISSN 1468-9367.

\bibitem{DKPT}
 Dubcov\'a, Kl\'{\i}\v{c}, Pokorn\'y, Turz\'{\i}k:
     Stability of Spatially Periodic Solutions in Coupled Map Lattices.
     Int.\,J.\,Bif.\,Chaos 
     Vol.13 No.2 (2003) 343-356 ISSN 0218-1274.

\bibitem{w10}
 Pro\v{c} vr\v{z}ou dve\v{r}e, 
 aneb Hopfova bifurkace ve va\v{s}\'{\i} lo\v{z}nici. 
 Ve sborn\'{\i}ku Matematika na vysok\'{y}ch \v{s}kol\'{a}ch, 
 Herbertov 1.-3.9.2003 ed. L. Herrmann.
     
\bibitem{w11}
 \'{U}vod do teorie neline\'{a}rn\'{\i}ch dynamick\'{y}ch 
 syst\'{e}m\accent23u{} a deterministick\'{e}ho chaosu. 
 Sborn\'{\i}k 6. ro\v{c}n\'{\i}ku semin\'{a}\v{r}e 
 Matematika na vysok\'{y}ch \v{s}kol\'{a}ch
 (Determinismus a n\'{a}hoda) Herbertov 5.-7.9.2005 ed. L. Herrmann.

\bibitem{w12}
 \'{U}skal\'{\i} \v{c}\'{\i}slicov\'{y}ch po\v{c}\'{\i}ta\v{c}\accent23u{}
  p\v{r}i simulaci neline\'{a}rn\'{\i}ch dynamick\'{y}ch 
     syst\'{e}m\accent23u{}.  
 Sborn\'{\i}k 6. ro\v{c}n\'{\i}ku semin\'{a}\v{r}e 
 Matematika na vysok\'{y}ch \v{s}kol\'{a}ch
 (Determinismus a n\'{a}hoda) Herbertov 5.-7.9.2005 ed. L. Herrmann.

\bibitem{w13}
 Modelov\'{a}n\'{\i} chaotick\'{y}ch syst\'{e}m\accent23u{} 
 kone\v{c}n\'{y}mi automaty.
 Sborn\'{\i}k 6. ro\v{c}n\'{\i}ku semin\'{a}\v{r}e 
 Matematika na vysok\'{y}ch \v{s}kol\'{a}ch
 (Determinismus a n\'{a}hoda) Herbertov 5.-7.9.2005 ed. L. Herrmann.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{\bf Invited lectures:}
\end{center}

\bibitem{i1}
Blinking Dynamical Systems and Campbell-Baker-Hausdorff formula.
     9.Int. Colloquium on Differential Equations, Plovdiv Bulgaria 18-23.8.1998.

\bibitem{i2}
 Nonlinear Dynamical Systems and Prediction.
     Nostradamus Prediction Workshop Zl\'{\i}n, Czech Republic 22-23.9.1998.

\bibitem{i3}
 Nonlinear Dynamics and Prediction.
     Nostradamus Prediction Conference 7-8.10.1999 Zl\'{\i}n, Czech Republic.

\bibitem{i4}
 History and Future of Prediction.
     Nostradamus Int. Conference, 
     25.-26. 9. 2001 UTB Zl\'{\i}n, Czech Republic.

\bibitem{i5}
 Chaos je kr\'{a}sn\'{y}.
     21.5.2002, KMOP MFF UK.

\bibitem{i6}
  Deterministick\'{y} chaos a matematick\'{e}, fyzik\'{a}ln\'{\i} 
  a filozofick\'{e} souvislosti. J\v{C}MF na MFF UK Praha, 19.5.2004.

\bibitem{i7}
 \'{U}vod do teorie neline\'{a}rn\'{\i}ch dynamick\'{y}ch 
 syst\'{e}m\accent23u{} a deterministick\'{e}ho chaosu. 
 6. ro\v{c}n\'{\i}k semin\'{a}\v{r}e 
 Matematika na vysok\'{y}ch \v{s}kol\'{a}ch
 (Determinismus a n\'{a}hoda) Herbertov 5.-7.9.2005.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{\bf Contributed lectures:}
\end{center}

\bibitem{c1}
  Coupled Map Lattices. 
     Summer School Mathematical Modeling in Chemical Engineering.
     5.-9.9.1988 Se\v{c}, Czech Republic.

\bibitem{c2}
  Cellular Automata. 
     Summer School Mathematical Modeling in Chemical Engineering.
     4.-8.9.1989 Ka\v{s}persk\'{e} Hory, Czech Republic.

\bibitem{c3}
  Excitable Dynamical Systems - Experiment and Modeling.
     Summer School Mathematical Modeling in Chemical Engineering.
     10.-14.9.1990 Harrachov, Czech Republic.

\bibitem{c4}
  Excitability in Chemical and Biological Systems.
     Department of Physical Chemistry, Free University of Brussels, Belgium.
     September 1991.

\bibitem{c5}
  Reduction of Dimensionality in Excitable Dynamical Systems.
     Tempus Summer School Mathematical Modeling in Chemical Engineering.
     24.-31.8.1992 Prague, Czech Republic.

\bibitem{c6}
  Origin and Dimension of Strange Non-chaotic Attractors.
     Euromech Colloquium 20.-25.9.1993 Spala, Poland.

\bibitem{c7}
  Non-chaotic Strange Attractors in Quasi-periodically Forced Excitable 
     Systems.
     12.-15.6.1994 Woudschotten, Netherlands.

\bibitem{c8}
 On Dynamical Systems Generated by Two Alternating Vector Fields.
     Second Int. Conference on Difference Equations and Applications,
     7.-11.8.1995, Veszprem, Hungary.

\bibitem{c9}
 Initiation and Control of Chaos by Blinking - Results.
     Int. Conference on Nonlinear Dynamics, Chaotic and Complex Systems,
     7.-12.11.1995, Zakopane, Poland.

\bibitem{c10}
 Evaluation of Time Series.
     Dynamics Days Edinburgh 28.6.-1.7.1998.

\bibitem{c11}
 Dynamical Systems with Blinking.
     Dynamics Days Edinburgh 28. 6.-1. 7. 1998.

\bibitem{c12}
 Zig-Zag Dynamical Systems.
     Dynamics Days, Como, Italy 20-23.6.1999.

\bibitem{c13}
 Mathematical Modeling of Chemical and Biological Systems.
     Bio-ma\-the\-ma\-ti\-cal Conference AV\v{C}R, Srn\'{\i}, 
     Czech Republic 8.-9.4.2000.

\bibitem{c14}
 Dynamical Systems Generated by Two Alternating Vector Fields and the 
     Baker-Campbell-Hausdorff formula.
     Dynamics Days, University of Surrey, England, 25.-29.6.2000.

\bibitem{c15}
 Excitable Dynamical Systems.
     XII Czech-Slovak-Polish Optical Conference Wave and Quantum Aspects of 
     Contemporary Optics, Velke Losiny, Czech Republic, 12.-15.9. 2000.

\bibitem{c16}
 Dynamical Systems and the Baker-Campbell-Hausdorff Formula.
     Conference Nolta 2000, 17.-21.9.2000 Dresden, Germany.

\bibitem{c17}
 Mathematical Modeling of Multi-Patch Predator-Prey Model.
     Bio-ma\-the\-ma\-ti\-cal Seminar 30.3.-1.4.2001 Cikhaj, Czech Republic.

\bibitem{c18}
 On Dynamical Systems Generated by Two Vector Fields.
     3rd Scientific Colloquium, Math. Dept. 
     Prague Institute of Chemical Technology, 26.-28.6.2001.

\bibitem{c19}
 Zig-zag Dynamical Systems and the Baker-Campbell-Hausdorff formula.
     The Czechoslovak Int. Conference on Differential
     Equations and Their Applications EQUADIFF 10, Prague, 28.-31.8.2001.

\bibitem{c20}
 Pru\v{z}n\'{e} kyvadlo a deterministick\'{y} chaos.
 Matematika na vy\-so\-k\'{y}ch \v{s}kol\'{a}ch, Herbertov 1.-3.9.2003.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{\bf My software packages published under the GNU General Public License:}
\end{center}

\bibitem{xpplot} {\tt xpplot} graphical data plotting tool for Unix.
\bibitem{easynum} {\tt  easynum} numerical library in C.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{\bf Publications citing my work:}
\end{center}

\bibitem{C1}
\v{S}t\v{e}p\'anek F., Kub\'{\i}\v{c}ek M., Marek M.:
Modeling and Optimization of Continuous Adsorption in a desiccant Rotor.
Ind. Eng. Chem. Res. {\bf 37}, 1435 (1998) 
cites \cite{KP96}.

\bibitem{C2}
Galan J., Freire E.:
Chaos in a Mean Field Model of Coupled Quantum Wells; 
Bifurcations of Periodic Orbits in a Symmetric Hamiltonian System.
Rep. Math. Phys. {\bf 44}, 87 (1999)
cites \cite{easynum}.

\bibitem{C3}
Russo L., Mancusi E., Maffettone P.L., Crescitelli S.:
Symmetry Properties and Bifurcation Analysis of a Class of 
Periodically Forced Chemical Reactors.
Chem. Eng. Sci. {\bf 57}, 5065 (2002)
cites \cite{KP96}.

\bibitem{C4}
Almeida J., Peralta-Salas D., Romera M.:
Can Two Chaotic Systems Give Rise to Order?
Physica {\bf D 200}, 124 (2005)
cites \cite{KP96}, \cite{Maslo2002}.

\bibitem{C5}
Afraimovich V.:
Some Topological Properties of Lattice Dynamical Systems.
In Chazottes J.\,R., Fernandez B. (Eds.): 
{\it Dynamics of Coupled Map Lattices and of Related Spatially
Extended Systems}. 
Series: Lecture Notes in Physics, {\bf 671}, Springer-Verlag (2005)
cites \cite{DKPT}.

\bibitem{C6}
Greve G. C. F.: Brave GNU World, Issue 5 (1999)
cites \cite{xpplot}.

\begin{center}
{\bf This PhD thesis has been announced in:}
\end{center}

\bibitem{C7}
Brindley J. et.\,al. (Eds.): UK Nonlinear News, Issue 42 (2006).

\end{thebibliography}

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
