diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex index 01ecc36..f11d5aa 100644 --- a/docs/nstrophy_doc.tex +++ b/docs/nstrophy_doc.tex @@ -382,62 +382,81 @@ The enstrophy is defined as 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)}) + \hat u(t_{n+1})=\mathfrak F_{t_n}(\hat u(t_n)) . \end{equation} -We then compute +Let $D\mathfrak F_{t_{n}}$ be the Jacobian of this map, in which we split the real and imaginary parts: if \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) + \hat u_k(t_n)=:\rho_k+i\iota_k + ,\quad + \mathfrak F_{t_n}(\hat u(t_n))_k=:\phi_k+i\psi_k \end{equation} -and, by\-~(\ref{T}), +then \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} + (D\mathfrak F_{t_n})_{k,p}:=\left(\begin{array}{cc} + \partial_{\rho_p}\phi_k&\partial_{\iota_p}\phi_k\\ + \partial_{\rho_p}\psi_k&\partial_{\iota_p}\psi_k + \end{array}\right) +\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} - \frac1n\log\mathrm{spec}(\pi_i) + \frac1{t_n}\log\mathrm{spec}(\pi_{t_n}) ,\quad - \pi_i:=\prod_{i=1}^nD\hat u^{(i)} + \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\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. +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} - =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)} - . + =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)}) diff --git a/src/lyapunov.c b/src/lyapunov.c index 8cf2e56..a4782cd 100644 --- a/src/lyapunov.c +++ b/src/lyapunov.c @@ -21,7 +21,8 @@ limitations under the License. #include #include -#define MATSIZE (K1*(2*K2+1)+K2) +#define MATSIZE (2*(K1*(2*K2+1)+K2)) +#define USIZE (K1*(2*K2+1)+K2) // compute Lyapunov exponents int lyapunov( @@ -50,8 +51,8 @@ int lyapunov( ){ double* lyapunov; double* lyapunov_avg; - _Complex double* Du_prod; - _Complex double* Du; + double* Du_prod; + double* Du; _Complex double* u; _Complex double* prevu; _Complex double* tmp1; @@ -63,7 +64,7 @@ int lyapunov( _Complex double* tmp7; _Complex double* tmp8; _Complex double* tmp9; - _Complex double* tmp10; + double* tmp10; double time; fft_vect fft1; fft_vect fft2; @@ -82,9 +83,6 @@ int lyapunov( double step=delta; double next_step=step; - const _Complex double one=1; - const _Complex double zero=0; - int i; // init average for (i=0; i u_{n+1} int ns_D_step( - _Complex double* out, + double* out, _Complex double* un, _Complex double* un_next, int K1, @@ -218,7 +216,7 @@ int ns_D_step( for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){ // derivative in the real direction // finite difference vector - for (i=0;i