Correctly deal with complex values in lyapunov exponents documentation

docs/nstrophy_doc.tex

@@ -103,9 +103,97 @@ We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k
   \mathcal K:=\{(k_1,k_2),\ |k_1|\leqslant K_1,\ |k_2|\leqslant K_2\}
   .
 \end{equation}
+\bigskip
 
-\subsubsection{Runge-Kutta methods}.
+{\bf Remark}:
-To solve the equation numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$.
+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}
@@ -190,30 +278,8 @@ It can be made by specifying the parameter {\tt adaptive\_cost}.
   These cost functions are selected by choosing {\tt adaptive\_cost=k3} and {\tt adaptive\_cost=k32} respectively.
 \end{itemize}
 
-\subsubsection{Reality}.
-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$.
 
-\subsubsection{FFT}. We compute T using a fast Fourier transform, defined as
+\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}
@@ -235,7 +301,7 @@ in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase
   \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_q\hat u_p
+  (q\cdot p^\perp)\frac{|q|}{|p|}\hat u_p\hat u_q
 \end{equation}
 provided
 \begin{equation}
@@ -257,6 +323,11 @@ Therefore,
   \right)(k)
 \end{equation}
 
+\subsection{Observables}
+\indent
+We define the following observables.
+\bigskip
+
 \subsubsection{Energy}.
 We define the energy as
 \begin{equation}
@@ -367,49 +438,124 @@ The enstrophy is defined as
   .
 \end{equation}
 
-\subsubsection{Lyapunov exponents}
+\subsection{Lyapunov exponents}
 \indent
-To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$.
+The Lyapunov are defined from the {\it tangent flow} of the dynamics.
-This map is always of Runge-Kutta type, that is,
+Consider an equation of the form
 \begin{equation}
-  \hat u(t_{n+1})=\mathfrak F_{t_n}(\hat u(t_n))
+  \dot u=f(t;u)
   .
 \end{equation}
-Let $D\mathfrak F_{t_{n}}$ be the Jacobian of this map, in which we split the real and imaginary parts: if
+Now, the flow may not be complex-differentiable, so the tangent flow should be computed on the real and imaginary parts.
+Let
 \begin{equation}
-  \hat u_k(t_n)=:\rho_k+i\iota_k
+  u=\zeta+i\xi
   ,\quad
-  \mathfrak F_{t_n}(\hat u(t_n))_k=:\phi_k+i\psi_k
+  f(t;u)=\theta(t;\zeta,\xi)+i\psi(t;\zeta,\xi)
-\end{equation}
+  .
-then
+\end{equation}
-\begin{equation}
+The tangent flow is given by
-  (D\mathfrak F_{t_n})_{k,p}:=\left(\begin{array}{cc}
+\begin{equation}
-    \partial_{\rho_p}\phi_k&\partial_{\iota_p}\phi_k\\
+  \dot\delta=\left(\begin{array}{cc}
-    \partial_{\rho_p}\psi_k&\partial_{\iota_p}\psi_k
+    D_\zeta\theta&D_\xi\theta\\
-  \end{array}\right)
+    D_\zeta\psi&D_\xi\psi
-\end{equation}
+  \end{array}\right)\delta
-We compute this Jacobian numerically using a finite difference, by computing
+\end{equation}
-\begin{equation}
+where $D_\zeta\theta$ is the Jacobian of $\theta$ with respect to $\zeta$ and so forth...
-  (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
+The flow of this equation is denoted by
-  \left(\begin{array}{cc}
+\begin{equation}
-    \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
+  \varphi_{t_0,t_1}(\delta_0)
-    \psi_k(\hat u+\epsilon\delta_p)-\psi_k(\hat u)&\psi_k(\hat u+i\epsilon\delta_p)-\psi_k(\hat u)
+\end{equation}
-  \end{array}\right)
+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.
+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}
-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)
@@ -428,88 +574,97 @@ The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
 %  \right)\hat u_{k-\ell}
 %  .
 %\end{equation}
-\bigskip
+For the reversible equation,
-
-\indent
-The Lyapunov exponents are then defined as
 \begin{equation}
-  \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
+  f(\hat u)=
-  ,\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}
-
-\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)^*$.
-After imposing this condition, the algorithm no longer blows up, but it is still unstable (for instance, increasing $K_1$ or $K_2$ leads to very different results).
-
-\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:
+so
 \begin{equation}
-  \sum_{k\in\mathbb Z^2}k^2\hat U_k\cdot\partial_t\hat U_k
+  ((D f(\hat u))\delta)_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}\alpha(\hat u) k^2\delta_k
-  +\frac{4\pi^2}{L^2}\sum_{k\in\mathbb Z^2}|k|\hat u_k^*T(\hat u,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}
-and so
+where
 \begin{equation}
-  \alpha(\hat u)
+  D\alpha(\hat u)\delta
-  =\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}
+  =
-  .
+  \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}
-  \label{alpha}
+  -\alpha(\hat u)\frac{2\sum_kk^4\delta_k^*\hat u_k}{\sum_kk^4|\hat u_k|^2}
-\end{equation}
-Note that, by\-~(\ref{realu})-(\ref{realT}),
-\begin{equation}
-  \alpha(\hat u)\in\mathbb R
   .
 \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