From 97a127a8db787c57c5bae093e0b625ca392019ac Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@rutgers.edu>
Date: Mon, 19 Feb 2024 11:41:35 -0500
Subject: [PATCH] Document Lyapunov exponents

---
 docs/nstrophy_doc.tex | 70 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 70 insertions(+)

diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex
index 9092dea..d7ad7b3 100644
--- a/docs/nstrophy_doc.tex
+++ b/docs/nstrophy_doc.tex
@@ -15,6 +15,7 @@
 \documentclass{ian}
 
 \usepackage{largearray}
+\usepackage{dsfont}
 
 \begin{document}
 
@@ -86,6 +87,7 @@ and $q\cdot p^\perp$ is antisymmetric under exchange of $q$ and $p$. Therefore,
   \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
@@ -375,6 +377,74 @@ The enstrophy is defined as
   .
 \end{equation}
 
+\subsubsection{Lyapunov exponents}
+\indent
+To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$.
+This map is always of Runge-Kutta type, that is,
+\begin{equation}
+  \hat u^{(n+1)}=\hat u^{(n)}+\delta\sum_{i=1}^q w_i\mathfrak F(\hat u^{(n)})
+\end{equation}
+(see\-~(\ref{ins_k})), 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
+  =
+  \sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}}
+  (q\cdot \ell^\perp)\left(
+    \frac{|q|}{|\ell|}
+    -
+    \frac{|\ell|}{|q|}
+  \right)\hat u_q
+  .
+\end{equation}
+\bigskip
+
+\indent
+The Lyapunov exponents are then defined as
+\begin{equation}
+  \frac1n\log\mathrm{spec}(\pi_i)
+  ,\quad
+  \pi_i:=\prod_{i=1}^nD\hat u^{(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\in\mathbb N_*$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $n$ is 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}
+  =Q^{(\alpha)}R^{(\alpha)}
+\end{equation}
+where $Q^{(\alpha)}$ is diagonal and $R^{(\alpha)}$ is an orthogonal matrix, and we divide by $Q^{(\alpha)}$ (thus only keeping $R^{(\alpha)}$).
+Thus, we replace
+\begin{equation}
+  \pi_{\alpha\mathfrak L_r+i}\mapsto R^{(\alpha)}\prod_{j=1}^iD\hat u^{(j)}
+  .
+\end{equation}
+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}
+
 \subsubsection{Numerical instability}.
 In order to prevent the algorithm from blowing up, it is necessary to impose the reality of $u(x)$ by hand, otherwise, truncation errors build up, and lead to divergences.
 It is sufficient to ensure that the convolution term $T(\hat u,k)$ satisfies $T(\hat u,-k)=T(\hat u,k)^*$.