% Copyright 2017-2025 Ian Jauslin
% 
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
% 
%     http://www.apache.org/licenses/LICENSE-2.0
% 
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

\documentclass{ian}

\usepackage{largearray}
\usepackage{dsfont}

\begin{document}

\hbox{}
\hfil{\bf\LARGE
{\tt nstrophy}
}
\vfill

\tableofcontents

\vfill
\eject

\setcounter{page}1
\pagestyle{plain}

\section{Description of the computation}
\subsection{Irreversible equation}
\indent Consider the incompressible Navier-Stokes equation in 2 dimensions
\begin{equation}
  \partial_tU=\nu\Delta U+G-(U\cdot\nabla)U,\quad
  \nabla\cdot U=0
  \label{ins}
\end{equation}
in which $G$ is the forcing term.
We take periodic boundary conditions, so, at every given time, $U(t,\cdot)$ is a function on the torus $\mathbb T^2:=\mathbb R^2/(L\mathbb Z)^2$. We represent $U(t,\cdot)$ using its Fourier series
\begin{equation}
  \hat U_k(t):=\frac1{L^2}\int_{\mathbb T^2}dx\ e^{i\frac{2\pi}L kx}U(t,x)
\end{equation}
for $k\in\mathbb Z^2$, and rewrite~\-(\ref{ins}) as
\begin{equation}
  \partial_t\hat U_k=
  -\frac{4\pi^2}{L^2}\nu k^2\hat U_k+\hat G_k
  -i\frac{2\pi}L\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
  (q\cdot\hat U_p)\hat U_q
  ,\quad
  k\cdot\hat U_k=0
  \label{ins_k}
\end{equation}
We then reduce the equation to a scalar one, by writing
\begin{equation}
  \hat U_k=\frac{i2\pi k^\perp}{L|k|}\hat u_k\equiv\frac{i2\pi}{L|k|}(-k_y\hat u_k,k_x\hat u_k)
  \label{udef}
\end{equation}
in terms of which, multiplying both sides of the equation by $\frac L{i2\pi}\frac{k^\perp}{|k|}$,
\begin{equation}
  \partial_t\hat u_k=
  -\frac{4\pi^2}{L^2}\nu k^2\hat u_k
  +\hat g_k
  +\frac{4\pi^2}{L^2|k|}\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
  \frac{(q\cdot p^\perp)(k^\perp\cdot q^\perp)}{|q||p|}\hat u_p\hat u_q
  \label{ins_k}
\end{equation}
with
\begin{equation}
  \hat g_k:=\frac{Lk^\perp}{2i\pi|k|}\cdot\hat G_k
  .
  \label{gdef}
\end{equation}
Furthermore
\begin{equation}
  (q\cdot p^\perp)(k^\perp\cdot q^\perp)
  =
  (q\cdot p^\perp)(q^2+p\cdot q)
\end{equation}
and $q\cdot p^\perp$ is antisymmetric under exchange of $q$ and $p$. Therefore,
\begin{equation}
  \partial_t\hat u_k=
  -\frac{4\pi^2}{L^2}\nu k^2\hat u_k+\hat g_k
  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
  =:\mathfrak F_k(\hat u)
  \label{ins_k}
\end{equation}
with
\begin{equation}
  T(\hat u,k):=
  \sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
  \frac{(q\cdot p^\perp)|q|}{|p|}\hat u_p\hat u_q
  .
  \label{T}
\end{equation}
We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k_2|>K_2$. Let
\begin{equation}
  \mathcal K:=\{(k_1,k_2),\ |k_1|\leqslant K_1,\ |k_2|\leqslant K_2\}
  .
\end{equation}
\bigskip

{\bf Remark}:
Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so
\begin{equation}
  \hat u_{-k}=\hat u_k^*
  .
  \label{realu}
\end{equation}
Similarly,
\begin{equation}
  \hat g_{-k}=\hat g_k^*
  .
  \label{realg}
\end{equation}
Thus,
\begin{equation}
  T(\hat u,-k)
  =
  T(\hat u,k)^*
  .
  \label{realT}
\end{equation}
In order to keep the computation as quick as possible, we only compute and store the values for $k_1\geqslant 0$.
(In fact, if we do not enforce the reality conditions, the computation has been found to be unstable.)

\subsection{Reversible equation}
\indent The reversible equation is similar to\-~(\ref{ins}) but instead of fixing the viscosity, we fix the enstrophy\-~\cite{Ga22}.
It is defined directly in Fourier space:
\begin{equation}
  \partial_t\hat U_k=
  -\frac{4\pi^2}{L^2}\alpha(\hat U) k^2\hat U_k+\hat G_k
  -i\frac{2\pi}L\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
  (q\cdot\hat U_p)\hat U_q
  ,\quad
  k\cdot\hat U_k=0
\end{equation}
where $\alpha$ is chosen such that the enstrophy is constant.
In terms of $\hat u$\-~(\ref{udef}), (\ref{gdef}), (\ref{T}):
\begin{equation}
  \partial_t\hat u_k=
  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\hat u_k
  +\hat g_k
  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
  .
  \label{rns_k}
\end{equation}
To compute $\alpha$, we use the constancy of the enstrophy:
\begin{equation}
  \sum_{k\in\mathbb Z^2}k^2\hat U_k\cdot\partial_t\hat U_k
  =0
\end{equation}
which, in terms of $\hat u$ is
\begin{equation}
  \sum_{k\in\mathbb Z^2}k^2\hat u_k^*\partial_t\hat u_k
  =0
\end{equation}
that is
\begin{equation}
  \frac{4\pi^2}{L^2}\alpha(\hat u)\sum_{k\in\mathbb Z^2}k^4|\hat u_k|^2
  =
  \sum_{k\in\mathbb Z^2}k^2\hat u_k^*\hat g_k
  +\frac{4\pi^2}{L^2}\sum_{k\in\mathbb Z^2}|k|\hat u_k^*T(\hat u,k)
\end{equation}
and so
\begin{equation}
  \alpha(\hat u)
  =\frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat u_k^*\hat g_k+\sum_k|k|\hat u_k^*T(\hat u,k)}{\sum_kk^4|\hat u_k|^2}
  .
  \label{alpha}
\end{equation}
Note that, by\-~(\ref{realu})-(\ref{realT}),
\begin{equation}
  \alpha(\hat u)\in\mathbb R
  .
\end{equation}

\subsection{Runge-Kutta methods}.
To solve these equations numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$.
These algorithms approximate the solution to an equation of the form $\dot u=f(t;u)$ with
\begin{equation}
  \hat u_k^{(n+1)}=\hat u_k^{(n)}
  +\delta_n\sum_{i=1}^sb_ik_i(t_n;\hat u^{(n)})
  ,\quad
  k_i(t_n;\hat u^{(n)}):=f\left(t_n+c_i\delta_n;\ \hat u+\delta_n\sum_{j=1}^{i-1}a_{i,j}k_j(t_n,\hat u^{(n)})\right)
  .
\end{equation}
The $c_i$ and $a_{i,j}$ are chosen in one of various ways, depending on the desired accuracy.
\bigskip

\indent
{\tt nstrophy} supports the 4th order Runge-Kutta ({\tt RK4}) and 2nd order Runge-Kutta ({\tt RK2}) algorithms.
In addition, several variable step methods are implemented:
\begin{itemize}
  \item the Runge-Kutta-Dormand-Prince method ({\tt RKDP54}), which is of 5th order, and adjusts the step by comparing to a 4th order method;
  \item the Runge-Kutta-Fehlberg method ({\tt RKF45}), which is of 4th order, and adjusts the step by comparing to a 5th order method;
  \item the Runge-Kutta-Bogacki-Shampine method ({\tt RKBS32}), which is of 3d order, and adjusts the step by comparing to a 2nd order method.
\end{itemize}
In these adaptive step methods, two steps are computed at different orders: $\hat u_k^{(n)}$ and $\hat U_k^{(n)}$, the step size is adjusted at every step in such a way that the error is small enough:
\begin{equation}
  D(\hat u^{(n)},\hat U^{(n)})
  <\epsilon_{\mathrm{target}}
\end{equation}
for some given {\it cost function} $D$, and $\epsilon_{\mathrm{target}}$, set using the {\tt adaptive\_tolerance} parameter.
The choice of the cost function matters, and will be discussed below.
If the error is larger than the target, then the step size is decreased.
How this is done depends on the order of algorithm.
If the order is $q$ (here we mean the smaller of the two orders, so 4 for {\tt RKDP54} and {\tt RKF45} and 2 for {\tt RKBS32}), then we expect (as long as the cost function is such that $D(\hat u,\hat u+\varphi)\sim\|\varphi\|$ in some norm)
\begin{equation}
  D(\hat u^{(n)},\hat U^{(n)})=\delta_n^qC_n
\end{equation}
for some number $C_n$.
We wish to set $\delta_{n+1}$ so that
\begin{equation}
  \delta_{n+1}^qC_n=\epsilon_{\mathrm{target}}
\end{equation}
so
\begin{equation}
  \delta_{n+1}
  =\left(\frac{\epsilon_{\mathrm{target}}}{C_n}\right)^{\frac1q}
  =\delta_n\left(\frac{\epsilon_{\mathrm{target}}}{D(\hat u^{(n)},\hat U^{(n)})}\right)^{\frac1q}
  .
  \label{adaptive_delta}
\end{equation}
(Actually, to be safe and ensure that $\delta$ decreases sufficiently, we multiply this by a safety factor that can be set using the {\tt adaptive\_factor} parameter.)
If the error is smaller than the target, we increase $\delta$ using\-~(\ref{adaptive_delta}) (without the safety factor).
To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} parameter.
\bigskip

\indent
The choice of the cost function $D$ matters.
It can be made by specifying the parameter {\tt adaptive\_cost}.
\begin{itemize}
  \item
  For computations where the main focus is the enstrophy\-~(\ref{enstrophy}), one may want to set the cost function to the relative difference of the enstrophies:
  \begin{equation}
    D(\hat u,\hat U):=\frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)}
    .
  \end{equation}
  This cost function is selected by choosing {\tt adaptive\_cost=enstrophy}.

  \item
  For computations where the main focus is the value of $\alpha$\-~(\ref{alpha}), one may want to set the cost function to the relative difference of $\alpha$:
  \begin{equation}
    D(\hat u,\hat U):=\frac{|\alpha(\hat u)-\alpha(\hat U)|}{|\alpha(\hat u)|}
    .
  \end{equation}
  This cost function is selected by choosing {\tt adaptive\_cost=alpha}.

  \item Alternatively, one my take $D$ to be the normalized $L_1$ norm:
  \begin{equation}
    D(\hat u,\hat U):=
    \frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|
    ,\quad
    \mathcal N:=\sum_k|\hat u_k|
    .
  \end{equation}
  This function is selected by choosing {\tt adaptive\_cost=L1}.

  \item Empirically, we have found that $|\hat u-\hat U|$ behaves like $k^{-3}$ for {\tt RKDP54} and {\tt RKF45}, and like $k^{-\frac32}$ for {\tt RKBS32}, so a cost function of the form
  \begin{equation}
    D(\hat u,\hat U):=\frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|k^{-3}
    ,\quad
    \mathcal N:=\sum_k|\hat u_k|k^{-3}
  \end{equation}
  or
  \begin{equation}
    D(\hat u,\hat U):=\frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|k^{-\frac32}
    ,\quad
    \mathcal N:=\sum_k|\hat u_k|k^{-\frac32}
  \end{equation}
  are sensible choices.
  These cost functions are selected by choosing {\tt adaptive\_cost=k3} and {\tt adaptive\_cost=k32} respectively.
\end{itemize}


\subsection{Computation of $T$: FFT}. We compute T using a fast Fourier transform, defined as
\begin{equation}
  \mathcal F(f)(n):=\sum_{m\in\mathcal N}e^{-\frac{2i\pi}{N_1}m_1n_1-\frac{2i\pi}{N_2}m_2n_2}f(m_1,m_2)
\end{equation}
where
\begin{equation}
  \mathcal N:=\{(n_1,n_2),\ 0\leqslant n_1< N_1,\ 0\leqslant n_2< N_2\}
\end{equation}
for some fixed $N_1,N_2$. The transform is inverted by
\begin{equation}
  \frac1{N_1N_2}\mathcal F^*(\mathcal F(f))(n)=f(n)
\end{equation}
in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase.
\bigskip

\indent The condition $p+q=k$ can be rewritten as
\begin{equation}
  T(\hat u,k)
  =
  \sum_{p,q\in\mathcal K}
  \frac1{N_1N_2}
  \sum_{n\in\mathcal N}e^{-\frac{2i\pi}{N_1}n_1(p_1+q_1-k_1)-\frac{2i\pi}{N_2}n_2(p_2+q_2-k_2)}
  (q\cdot p^\perp)\frac{|q|}{|p|}\hat u_p\hat u_q
\end{equation}
provided
\begin{equation}
  N_i>3K_i.
\end{equation}
Indeed, $\sum_{n_i=0}^{N_i}e^{-\frac{2i\pi}{N_i}n_im_i}$ vanishes unless $m_i=0\%N_i$ (in which $\%N_i$ means `modulo $N_i$'), and, if $p,q,k\in\mathcal K$, then $|p_i+q_i-k_i|\leqslant3K_i$, so, as long as $N_i>3K_i$, then $(p_i+q_i-k_i)=0\%N_i$ implies $p_i+q_i=k_i$.
Therefore,
\begin{equation}
  T(\hat u,k)
  =
  \textstyle
  \frac1{N_1N_2}
  \mathcal F^*\left(
    \mathcal F\left(\frac{p_x\hat u_p}{|p|}\right)(n)
    \mathcal F\left(q_y|q|\hat u_q\right)(n)
    -
    \mathcal F\left(\frac{p_y\hat u_p}{|p|}\right)(n)
    \mathcal F\left(q_x|q|\hat u_q\right)(n)
  \right)(k)
\end{equation}

\subsection{Observables}
\indent
We define the following observables.
\bigskip

\subsubsection{Energy}.
We define the energy as
\begin{equation}
  E(t)=\frac12\int\frac{dx}{L^2}\ U^2(t,x)=\frac12\sum_{k\in\mathbb Z^2}|\hat U_k|^2
  .
\end{equation}
We have
\begin{equation}
  \partial_t E=\int\frac{dx}{L^2}\ U\partial tU
  =
  \nu\int\frac{dx}{L^2}\ U\Delta U
  +\int\frac{dx}{L^2}\ UG
  -\int\frac{dx}{L^2}\ U(U\cdot\nabla)U
  .
\end{equation}
Since we have periodic boundary conditions,
\begin{equation}
  \int dx\ U\Delta U=-\int dx\ |\nabla U|^2
  .
\end{equation}
Furthermore,
\begin{equation}
  I:=\int dx\ U(U\cdot\nabla)U
  =\sum_{i,j=1,2}\int dx\ U_iU_j\partial_jU_i
  =
  -\sum_{i,j=1,2}\int dx\ (\partial_jU_i)U_jU_i
  -\sum_{i,j=1,2}\int dx\ U_i(\partial_jU_j)U_i
\end{equation}
and since $\nabla\cdot U=0$,
\begin{equation}
  I
  =
  -I
\end{equation}
and so $I=0$.
Thus,
\begin{equation}
  \partial_t E=
  \int\frac{dx}{L^2}\ \left(-\nu|\nabla U|^2+UG\right)
  =
  \sum_{k\in\mathbb Z^2}\left(-\frac{4\pi^2}{L^2}\nu k^2|\hat U_k|^2+\hat U_{-k}\hat G_k\right)
  .
\end{equation}
Furthermore,
\begin{equation}
  \sum_{k\in\mathbb Z^2}k^2|\hat U_k|^2\geqslant
  \sum_{k\in\mathbb Z^2}|\hat U_k|^2-|\hat U_0|^2
  =2E-|\hat U_0|^2
\end{equation}
so
\begin{equation}
  \partial_t E\leqslant -\frac{8\pi^2}{L^2}\nu E+\frac{4\pi^2}{L^2}\nu\hat U_0^2+\sum_{k\in\mathbb Z^2}\hat U_{-k}\hat G_k
  \leqslant
  -\frac{8\pi^2}{L^2}\nu E+\frac{4\pi^2}{L^2}\nu\hat U_0^2+
  \|\hat G\|_2\sqrt{2E}
  .
\end{equation}
In particular, if $\hat U_0=0$ (which corresponds to keeping the center of mass fixed),
\begin{equation}
  \partial_t E\leqslant -\frac{8\pi^2}{L^2}\nu E+\|\hat G\|_2\sqrt{2E}
  .
\end{equation}
Now, if $\frac{8\pi^2}{L^2}\nu\sqrt E<\sqrt2\|\hat G\|_2$, then
\begin{equation}
  \frac{\partial_t E}{-\frac{8\pi^2}{L^2}\nu E+\|\hat G\|_2\sqrt{2E}}\leqslant1
\end{equation}
and so
\begin{equation}
  \frac{\log(1-\frac{8\pi^2\nu}{L^2\sqrt2\|\hat G\|_2}\sqrt{E(t)})}{-\frac{4\pi^2}{L^2}\nu}\leqslant t+
  \frac{\log(1-\frac{8\pi^2\nu}{L^2\sqrt2\|\hat G\|_2}\sqrt{E(0)})}{-\frac{4\pi^2}{L^2}\nu}
\end{equation}
and
\begin{equation}
  E(t)
  \leqslant
  \left(
    \frac{L^2\sqrt2\|\hat G\|_2}{8\pi^2\nu}(1-e^{-\frac{4\pi^2}{L^2}\nu t})
    +e^{-\frac{4\pi^2}{L^2}\nu t}\sqrt{E(0)}
  \right)^2
  .
\end{equation}
If $\frac{8\pi^2}{L^2}\nu\sqrt E>\sqrt2\|\hat G\|_2$,
\begin{equation}
  \frac{\partial_t E}{-\frac{8\pi^2}{L^2}\nu E+\|\hat G\|_2\sqrt{2E}}\geqslant1
\end{equation}
and so
\begin{equation}
  \frac{\log(\frac{8\pi^2\nu}{L^2\sqrt2\|\hat G\|_2}\sqrt{E(t)}-1)}{-\frac{4\pi^2}{L^2}\nu}\geqslant t+
  \frac{\log(\frac{8\pi^2\nu}{L^2\sqrt2\|\hat G\|_2}\sqrt{E(0)})-1}{-\frac{4\pi^2}{L^2}\nu}
\end{equation}
and
\begin{equation}
  E(t)
  \leqslant
  \left(
    \frac{L^2\sqrt2\|\hat G\|_2}{8\pi^2\nu}(1-e^{-\frac{4\pi^2}{L^2}\nu t})
    +e^{-\frac{4\pi^2}{L^2}\nu t}\sqrt{E(0)}
  \right)^2
  .
  \label{enstrophy}
\end{equation}

\subsubsection{Enstrophy}.
The enstrophy is defined as
\begin{equation}
  \mathcal En(t)=\int\frac{dx}{L^2}\ |\nabla U|^2
  =\frac{4\pi^2}{L^2}\sum_{k\in\mathbb Z^2}k^2|\hat U_k|^2
  .
\end{equation}

\subsection{Lyapunov exponents}
\indent
The Lyapunov are defined from the {\it tangent flow} of the dynamics.
Consider an equation of the form
\begin{equation}
  \dot u=f(t;u)
  .
\end{equation}
Now, the flow may not be complex-differentiable, so the tangent flow should be computed on the real and imaginary parts.
Let
\begin{equation}
  u=\zeta+i\xi
  ,\quad
  f(t;u)=\theta(t;\zeta,\xi)+i\psi(t;\zeta,\xi)
  .
\end{equation}
The tangent flow is given by
\begin{equation}
  \dot\delta=\left(\begin{array}{cc}
    D_\zeta\theta&D_\xi\theta\\
    D_\zeta\psi&D_\xi\psi
  \end{array}\right)\delta
\end{equation}
where $D_\zeta\theta$ is the Jacobian of $\theta$ with respect to $\zeta$ and so forth...
The flow of this equation is denoted by
\begin{equation}
  \varphi_{t_0,t_1}(\delta_0)
\end{equation}
and defined by
\begin{equation}
  \frac d{dt}\varphi_{t_0,t}(\delta_0)=\left(\begin{array}{cc}
    D_\zeta\theta(t;\zeta,\xi)&D_\xi\theta(t;\zeta,\xi)\\
    D_\zeta\psi(t;\zeta,\xi)&D_\xi\psi(t;\zeta,\xi)
  \end{array}\right)\varphi_{t_0,t}(\delta_0)
  ,\quad
  \varphi_{t_0,t_0}(\delta_0)=\delta_0
  .
\end{equation}
The flow $\varphi_{t_0,t}(\delta_0)$ is linear in $\delta_0$, and so can be represented as a matrix.
The Lyapnuov exponents are defined from the $QR$ decomposition of $\varphi_{0,t}$: if
\begin{equation}
  \varphi_{0,t}=Q_tR_t
\end{equation}
where $Q_t$ is orthogonal, and $R_t$ is triangular with positive diagonal entries $(r_t^{(j)})_j$, then the Lyapunov exponents are
\begin{equation}
  \lim_{t\to\infty}\frac 1t\log|r_t^{(j)}|
  .
\end{equation}
\bigskip

\indent
One problem to compute the Lyapunov exponents numerically is that the spectrum depends exponentially on time, and so has a tendency to blow up or shrink very rapidly, which leads to large truncation errors.
To avoid these, we compute the tangent flow {\it in parts}: we split time into intervals:
\begin{equation}
  [0,t)=\bigcup_{i=0}^{N-1}[t_i,t_{i+1})
  ,\quad
  \varphi_{0,t}=\prod_{i=N-1}^0\varphi_{t_{i},t_{i+1}}
  .
\end{equation}
At the thresholds between these intervals, we perform a QR decomposition:
\begin{equation}
  Q_0R_0:=\varphi_{0,t_0}
  ,\quad
  Q_iR_i:=\varphi_{t_{i-1},t_i}Q_{i-1}
\end{equation}
where $Q_i$ is orthogonal and $R_i$ is upper triangular with $\geqslant 0$ diagonal entries (in addition, we assume these are not $0$).
Doing so, we find
\begin{equation}
  \varphi_{0,t}=Q_{N-1}\prod_{i=N-1}^0R_i
\end{equation}
and since the product of triangular matrices is triangular and the diagonal elements multiply, we find that the Lyapunov exponents are given by
\begin{equation}
  \lim_{t_{N-1}\to\infty}\frac1{t_{N-1}}\sum_{i=N-1}^0\log|r_i^{(j)}|
  .
\end{equation}
To do the computation numerically, we drop the limit, and compute the logarithm of the product of the diagonal entries of $R_i$.
\bigskip

\indent
In practice, we approximate $\varphi_{t_{i-1},t_i}$ by running a Runge-Kutta algorithm for the tangent flow equation.
Because the tangent flow equation depends on $u(t)$, we must run it at times for which $u$ has been computed, so the Runge-Kutta algorithm for the tangent flow cannot be an adaptive step method, and should be one of {\tt RK4} for the fourth order algorithm or {\tt RK2} for the second order algorithm.
To obtain the full matrix, we consider every element of the canonical basis as an initial condition $\delta_0$.
We then iterate the Runge-Kutta algorithm until the time $t_0$ (chosen in one of two ways, see below), at which point we perform a QR decomposition, save the diagonal entries of $R$, replace the family of initial conditions with the columns of $Q$, and continue the flow from there.
The choice of the times $t_i$ can be done either by fixed-length intervals, specified with the option {\tt lyapunov\_reset}, or the QR decomposition can be triggered whenever $\|\delta\|_1$ exceeds a threshold, specified in {\tt lyapunov\_maxu} (after all, the intervals are used to prevent $\delta$ from becoming too large). 
\bigskip

\indent
To compute the Lyapunov exponents, we thus need the Jacobians of $\theta$ and $\psi$.
Note that, by the linearity of the tangent flow equation,
\begin{equation}
  ((D \theta(\hat u))\delta)_{k}
  =
  \mathcal Re(Df(\hat u)(\delta_{\mathrm r}+i\delta_{\mathrm i}))
  ,\quad
  ((D \psi(\hat u))\delta)_{k}
  =
  \mathcal Im(Df(\hat u)(\delta_{\mathrm r}+i\delta_{\mathrm i}))
  .
\end{equation}
For the irreversible equation,
\begin{equation}
  f(\hat u)=
  -\frac{4\pi^2}{L^2}\nu k^2\hat u_k+\hat g_k
  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
\end{equation}
and
\begin{equation}
  ((D f(\hat u))\delta)_{k}
  =
  -\frac{4\pi^2}{L^2}\nu k^2\delta_{k}
  +\frac{4\pi^2}{L^2|k|}DT(\hat u,k)\delta
\end{equation}
\begin{equation}
  DT(\hat u,k)\delta
  =
  \sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
  \left(\frac{(q\cdot p^\perp)|q|}{|p|}+\frac{(p\cdot q^\perp)|p|}{|q|}\right)\hat u_p(\delta_{q,\mathrm r}+i\delta_{q,\mathrm i})
  .
\end{equation}
%and, by\-~(\ref{T}),
%\begin{equation}
%  \partial_{\hat u_\ell}T(\hat u,k)
%  =
%  \sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}}
%  \left(
%    \frac{(q\cdot \ell^\perp)|q|}{|\ell|}
%    +
%    \frac{(\ell\cdot q^\perp)|\ell|}{|q|}
%  \right)\hat u_q
%  =
%  (k\cdot \ell^\perp)\left(
%    \frac{|k-\ell|}{|\ell|}
%    -
%    \frac{|\ell|}{|k-\ell|}
%  \right)\hat u_{k-\ell}
%  .
%\end{equation}
For the reversible equation,
\begin{equation}
  f(\hat u)=
  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\hat u_k
  +\hat g_k
  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
\end{equation}
so
\begin{equation}
  ((D f(\hat u))\delta)_k
  =
  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\delta_k
  -\frac{4\pi^2}{L^2}k^2\hat u_k D\alpha(\hat u)\delta
  +\frac{4\pi^2}{L^2|k|}DT(\hat u,k)\delta
\end{equation}
where
\begin{equation}
  D\alpha(\hat u)\delta
  =
  \frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat \delta_k^*\hat g_k+\sum_k|k|\hat (\delta_k^*T(\hat u,k)+\hat u_k^*DT(\hat u,k)\delta)}{\sum_kk^4|\hat u_k|^2}
  -\alpha(\hat u)\frac{2\sum_kk^4\delta_k^*\hat u_k}{\sum_kk^4|\hat u_k|^2}
  .
\end{equation}


%We compute this Jacobian numerically using a finite difference, by computing
%\begin{equation}
%  (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
%  \left(\begin{array}{cc}
%    \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
%    \psi_k(\hat u+\epsilon\delta_p)-\psi_k(\hat u)&\psi_k(\hat u+i\epsilon\delta_p)-\psi_k(\hat u)
%  \end{array}\right)
%  .
%\end{equation}
%The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
%%, so
%%\begin{equation}
%%  D\hat u^{(n+1)}=\mathds 1+\delta\sum_{i=1}^q w_iD\mathfrak F(\hat u^{(n)})
%%  .
%%\end{equation}
%%We then compute
%%\begin{equation}
%%  (D\mathfrak F(\hat u))_{k,\ell}
%%  =
%%  -\frac{4\pi^2}{L^2}\nu k^2\delta_{k,\ell}
%%  +\frac{4\pi^2}{L^2|k|}\partial_{\hat u_\ell}T(\hat u,k)
%%\end{equation}
%%and, by\-~(\ref{T}),
%%\begin{equation}
%%  \partial_{\hat u_\ell}T(\hat u,k)
%%  =
%%  \sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}}
%%  \left(
%%    \frac{(q\cdot \ell^\perp)|q|}{|\ell|}
%%    +
%%    \frac{(\ell\cdot q^\perp)|\ell|}{|q|}
%%  \right)\hat u_q
%%  =
%%  (k\cdot \ell^\perp)\left(
%%    \frac{|k-\ell|}{|\ell|}
%%    -
%%    \frac{|\ell|}{|k-\ell|}
%%  \right)\hat u_{k-\ell}
%%  .
%%\end{equation}
%\bigskip
%
%\indent
%The Lyapunov exponents are then defined as
%\begin{equation}
%  \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
%  ,\quad
%  \pi_{t_n}:=\prod_{i=1}^nD\hat u(t_i)
%  .
%\end{equation}
%However, the product of $D\hat u$ may become quite large or quite small (if the exponents are not all 1).
%To avoid this, we periodically rescale the product.
%We set $\mathfrak L_r>0$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $t_n$ crosses a multiple of $\mathfrak L_r$, we rescale the eigenvalues of $\pi_i$ to 1.
%To do so, we perform a $QR$ decomposition:
%\begin{equation}
%  \pi_{\alpha\mathfrak L_r}
%  =R^{(\alpha)}Q^{(\alpha)}
%\end{equation}
%where $Q^{(\alpha)}$ is orthogonal and $R^{(\alpha)}$ is a diagonal matrix, and we divide by $R^{(\alpha)}$ (thus only keeping $Q^{(\alpha)}$).
%The Lyapunov exponents at time $\alpha\mathfrak L_r$ are then
%\begin{equation}
%  \frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})
%  .
%\end{equation}



\vfill
\eject

\begin{thebibliography}{WWW99}
\small
\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}
\end{thebibliography}

\end{document}
