Reduced-order Description of Transient Instabilities and Computation of Finite-Time Lyapunov Exponents

TL;DR

This paper introduces a reduced-order method using OTD modes to efficiently detect and analyze finite-time instabilities in high-dimensional chaotic systems, significantly lowering computational costs.

Contribution

It develops a new reduced-order approach leveraging OTD modes for fast convergence to instability directions and efficient FTLE computation.

Findings

OTD modes converge exponentially to eigendirections of the Cauchy-Green tensor.

The reduced-order method drastically reduces computational costs.

Validated on two numerical examples.

Abstract

High-dimensional chaotic dynamical systems can exhibit strongly transient features. These are often associated with instabilities that have finite-time duration. Because of the finite-time character of these transient events, their detection through infinite-time methods, e.g. long term averages, Lyapunov exponents or information about the statistical steady-state, is not possible. Here we utilize a recently developed framework, the Optimally Time-Dependent (OTD) modes, to extract a time-dependent subspace that spans the modes associated with transient features associated with finite-time instabilities. As the main result, we prove that the OTD modes, under appropriate conditions, converge exponentially fast to the eigendirections of the Cauchy--Green tensor associated with the most intense finite-time instabilities. Based on this observation, we develop a reduced-order method for the computation of finite-time Lyapunov exponents (FTLE) and vectors. In high-dimensional systems, the computational cost of the reduced-order method is orders of magnitude lower than the full FTLE computation. We demonstrate the validity of the theoretical findings on two numerical examples.

Equations

$$\dot{z}=f(z,t),z\in {\mathbb{R}}^n,t\in I=[t_0,t_0+T],$$

$$F_{t_0}^t:$$

$$z_0\mapsto z(t;t_0,z_0),$$

$$\dot{v}=L(z(t),t)v,v(t)\in {\mathbb{R}}^n,t\in I=[t_0,t_0+T],$$

$$v(t)=\nabla F_{t_0}^t(z_0)v_0,\text{mbox}witht\in I=[t_0,t_0+T].$$

$$C_{t_0}^t=(\nabla F_{t_0}^t{)}^T\nabla F_{t_0}^t,$$

$$B_{t_0}^t=\nabla F_{t_0}^t(\nabla F_{t_0}^t{)}^T,$$

$$C_{t_0}^t(z_0){\xi }_i(t,t_0,z_0)={\lambda }_i(t,t_0,z_0){\xi }_i(t,t_0,z_0),i=1,\cdots ,n,$$

$$B_{t_0}^t(z_0){\eta }_i(t,t_0,z_0)={\mu }_i(t,t_0,z_0){\eta }_i(t,t_0,z_0),i=1,\cdots ,n.$$

$$\big{<}\xi_{i},\xi_{j}\big{>}=\big{<}\eta_{i},\eta_{j}\big{>}=\delta_{ij},\quad i,j=1,\dots,n,$$

$$B_{t_0}^t\left(\nabla F_{t_0}^t{\xi }_i\right)={\lambda }_i\left(\nabla F_{t_0}^t{\xi }_i\right).$$

$$\nabla F_{t_0}^t{\xi }_i=\sqrt{{\lambda }_i}{\eta }_i,i=1,\dots ,n,$$

$${\lambda }_1>{\lambda }_2>\cdots >{\lambda }_n>0.$$

$${\Lambda }_i(t_0+T,t_0,z_0)=\frac{1}{T}\log \sqrt{{\lambda }_i(t_0+T,t_0,z_0)},i=1,2,\cdots ,n.$$

$$\displaystyle\mathcal{F}=\frac{1}{\left(\delta t\right)^{2}}\sum_{i=1}^{r}\big{\|}u_{i}(t+\delta t)-\nabla F_{t}^{t+\delta t}u_{i}(t)\big{\|}^{2},\quad\quad\delta t\rightarrow 0,$$

$$⟨u_i\left(t\right),u_j\left(t\right)⟩={\delta }_{ij},i,j=1,...,r.$$

$$u_i(t+\delta t)$$

$$\nabla F_t^{t+\delta t}$$

$$\mathcal{F}({\dot{u}}_1,{\dot{u}}_2,\dots ,{\dot{u}}_r)=\sum _{i=1}^r{\Vert \frac{\partial u_i\left(t\right)}{\partial t}-L(z(t),t)u_i(t)\Vert }^2.$$

$$\frac{\partial U}{\partial t}=LU-U{U}^{T}LU,$$

$${\gamma }_{U,V}=\frac{\Vert U^{T}V\Vert }{r^{1/2}},$$

$$1=\Vert v_i{\Vert }^2=\Vert h_i{\Vert }^2+\Vert b_i{\Vert }^2,i=1,2,\cdots ,r.$$

$$\Vert U^{T}V{\Vert }^2=\Vert H{\Vert }^2=\sum _{i=1}^r\Vert h_i{\Vert }^2=r-\sum _{i=1}^r\Vert b_i{\Vert }^2=r-\Vert B{\Vert }^2.$$

$$\dot{T}=L_r(t)T,T(t_0)=T_0,$$

$$\underset{t\to \infty }{\lim }\frac{{\lambda }_i(t,t_0,z_0)}{e^{a_{i}t}}=K_i.$$

$$v_{0_{j}}=\sum_{i=1}^{n}\big{<}v_{0_{j}},\xi_{i}(t,t_{0},z_{0})\big{>}\xi_{i}(t,t_{0},z_{0}),\quad\ j=1,...,r.$$

$$v_{j}(t)=\sum_{i=1}^{n}\sqrt{\lambda_{i}(t,t_{0},z_{0})}\big{<}v_{0_{j}},\xi_{i}(t,t_{0},z_{0})\big{>}\eta_{i}(t,t_{0},z_{0}),\quad\ j=1,...,r.$$

$$\big{<}v_{j}(t),\eta_{k}(t,t_{0},z_{0})\big{>}=\sqrt{\lambda_{k}(t,t_{0},z_{0})}\big{<}v_{0_{j}},\xi_{k}(t,t_{0},z_{0})\big{>},\quad\ j,k=1,...,r.$$

$$\big{<}v_{j}(t),v_{j}(t)\big{>}=\sum_{i=1}^{n}{\lambda_{i}(t,t_{0},z_{0})}\big{<}v_{0_{j}},\xi_{i}(t,t_{0},z_{0})\big{>}^{2},\quad\ j=1,...,r.$$

$$\cos{\alpha_{v_{j},\eta_{k}}}=\frac{\big{<}v_{j}(t),\eta_{k}\big{>}}{\|v_{j}(t)\|}=\frac{\sqrt{\lambda_{k}}\big{<}v_{0_{j}},\xi_{k}\big{>}}{\sqrt{\sum_{i=1}^{n}\lambda_{i}\big{<}v_{0_{j}},\xi_{i}\big{>}^{2}}}.$$

$$r\gamma_{v,\eta}^{2}=\sum_{j=1}^{r}\sum_{k=1}^{r}\frac{\lambda_{k}\big{<}v_{0_{j}},\xi_{k}\big{>}^{2}}{\sum_{i=1}^{n}\lambda_{i}\big{<}v_{0_{j}},\xi_{i}\big{>}^{2}}=\sum_{j=1}^{r}\frac{\sum_{k=1}^{r}\lambda_{k}\big{<}v_{0_{j}},\xi_{k}\big{>}^{2}}{\sum_{i=1}^{n}\lambda_{i}\big{<}v_{0_{j}},\xi_{i}\big{>}^{2}}=\sum_{j=1}^{r}\frac{A_{j}}{A_{j}+B_{j}}.$$

Full text

Keywords

Finite-time instabilities; Optimally time-dependent modes; Dynamically orthogonal modes; Transient phenomena; Lagrangian coherent structures; finite-time Lyapunov exponents.

1 Introduction

For a plethora of dynamical systems transient phenomena — usually associated with finite-time instabilities — define the most important aspects of the response. Examples include chaotic fluid systems [3, 15], nonlinear waves [7, 8] and networks [35, 6]. For the analysis of such systems and the formulation of prediction and control algorithms it is critical to identify modes associated with these strongly transient features. A first important challenge towards identifying these modes is the high-dimensionality of the response. This high-dimensionality does not necessarily imply that the transient features are also high-dimensional. In fact it is often the case that the modes associated with transient features are very few. However, because of the broad spectrum of the response it is hard to identify them using energy-based criteria. An additional challenge is associated with the intrinsically time-dependent character of these features which is, in general, non-periodic and makes it essential to consider arbitrary time-dependence for these modes as well.

In the context of dynamical systems numerous approaches have been developed to characterize the transient features associated with non-periodic behavior. The notion of Lagrangian Coherent Structures (LCS) emerged from the fluid dynamics problem of mixing as a general method to characterize dynamical systems with arbitrary time-dependence [20, 18, 24]. The method quantifies finite-time instabilities through the explicit computation of the maximum eigenvalue of the Cauchy–Green tensor on every location of the phase space. For low-dimensional systems it provides a complete and unambiguous description of the finite-time instabilities by identifying parts of the phase space associated with such responses. However, given the fast growth of the computational cost with respect to the phase space dimensionality it is not possible to apply it for very high- or infinite-dimensional systems in order to characterize finite-time instabilities.

For such problems, other approaches focusing more on the description of the system attractor rather than the whole phase space have also been developed in the context of uncertainty quantification and model order reduction. Specifically, the notion of dynamical orthogonality [30, 31, 5, 4] allows for the computation of time-dependent, dynamically orthogonal (DO) modes that lead to reduced-order description of chaotic attractors with low-intrinsic dimensionality but with strong time-dependence [28, 29]. The same notion can be combined with stochastic closures [32, 34] in order to quantify statistics in time-dependent subspaces for turbulent dynamics systems with very broad spectra (or high intrinsic dimensionality) and strongly transient features [33, 25]. In all of these studies the effectiveness of the derived DO modes on capturing transient responses was demonstrated numerically and the very few theoretical studies related to this problem were constrained in very specific setups. One such case is the study of linear parabolic (stable) equations in [26] where it was shown that the DO modes tend to capture the most energetic directions. However, there is no general result indicating where the DO modes converge for an arbitrary (nonlinear) dynamical system.

To circumvent this limitation a minimization principle was introduced in [2] that allows for the derivation of equations that evolve a time-dependent set of modes, the Optimally-Time Dependent (OTD) modes, along a given trajectory of the system. These modes are identical with those obtained from the DO equations in the deterministic limit, i.e. for the case where the stochastic energy goes to zero. To this end, the OTD modes can be seen as the deterministic analog of DO modes. For sufficiently long times where the system reaches an equilibrium it was proven in [2] that the OTD modes converge to the most unstable directions of the system in the asymptotic limit. Moreover, the same modes were utilized in [15] to formulate predictors of extreme events for Navier-Stokes equations and nonlinear waves problems.

Despite their success on numerically capturing transient phenomena associated with finite-time instabilities the exact relation between OTD modes (or DO modes) and finite-time instabilities remains an open problem. In this work we provide a rigorous link between the OTD modes and finite-time instabilities, by showing that under general conditions the OTD modes converge exponentially fast to the eigendirections of the Cauchy–Green tensor associated with the largest eigenvalues, i.e. with the largest finite-time Lyapunov exponents. Thus, we prove that the OTD modes are able to extract over a finite-time interval the modes associated with the most intense finite-time instabilities or the most pronounced transient phenomena. Note, that such convergence does not depend on the dimensionality of the system nor the intrinsic dimensionality of the attractor.

Apart from the fundamental implications of this result on the formulation of reduced-order algorithms for the prediction and control of transient phenomena there is a direct application on the computation of the maximum finite-time Lyapunov exponents for high- or infinite-dimensional systems. Specifically, using the OTD modes we formulate a reduced-order framework that allows, under mild conditions, the computation of finite-time Lyapunov Exponents for general dynamical systems without having to compute the full Cauchy–Green tensor. We demonstrate the derived algorithm in two systems, the Arnold-Beltrami-Childress (ABC) flow and the Six-dimensional Charney-DeVore model with regime transitions. In both cases we examine the effectiveness of the reduced-order method and study limitations and applicability.

2 Preliminaries and Notation

We consider the system of differential equations

[TABLE]

where $f:U\times I\rightarrow\mathbb{R}^{n}$ is a sufficiently smooth vector field. Let $F_{t_{0}}^{t}$ with $t\in I$ denote the associated flow map,

[TABLE]

where $z(t;t_{0},z_{0})$ is a trajectory of system (1) with the initial condition $z_{0}$. The linearized system around the trajectory $z(t)\equiv z(t;t_{0},z_{0})$ satisfies the equation of variations,

[TABLE]

where $L(z,t):=\nabla_{z}f(z,t)$. The deformation gradient $\nabla F^{t}_{t_{0}}$ is the fundamental solution matrix for the linearized dynamics such that the solutions $v(t)\equiv v(t;t_{0},v_{0})$ of the equation of variations satisfy

[TABLE]

2.1 Finite-time Lyapunov exponents

To measure the growth of infinitesimal perturbations in the phase space we use the right Cauchy–Green strain tensor

[TABLE]

and the left Cauchy–Green strain tensor

[TABLE]

where $T$ denote matrix transposition. Let $\xi(t,t_{0},z_{0})\in\mathbb{R}^{n}$ denote the eigenvectors of the right Cauchy–Green strain tensor, so that

[TABLE]

where $\lambda_{i}(t,t_{0},z_{0})$ are the corresponding eigenvalues. Similarly, denote the eigenvectors of the left Cauchy–Green strain tensor by $\eta_{i}(t,t_{0},z_{0})\in\mathbb{R}^{n}$ and their corresponding eigenvalues by $\mu_{i}(t,t_{0},z_{0})$, so that

[TABLE]

When no confusion may arise, we omit the dependence of the eigenvalues and eigenvectors on $(t,t_{0},z_{0})$ for notational simplicity. Since the Cauchy–Green strain tensors are symmetric and positive definite, their eigenvalues are real and positive, $\lambda_{i},\mu_{i}\in\mathbb{R}^{+}$. Furthermore, the eigenvectors are orthogonal, i.e.,

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product.

It is straightforward to show that the right and left tensors have the same eigenvalues, $\lambda_{i}=\mu_{i}$. Specifically, considering the definition of the right Cauchy–Green tensor eigenvalues, $(\nabla F_{t_{0}}^{t})^{T}\nabla F_{t_{0}}^{t}\xi_{i}=\lambda_{i}\xi_{i}$, we multiply the equation with $\nabla F_{t_{0}}^{t}$ from the left to obtain

[TABLE]

Thus, the left and right Cauchy–Green strain tensors have the same set of eigenvalues. Finally, using the singular value decomposition of the deformation gradient, one can show the well-known relation that

[TABLE]

(see, e.g., [21]). In the following, we order the Cauchy–Green eigenvalues in a descending order,

[TABLE]

The finite-time Lyapunov exponents (FTLEs) of system (1) corresponding to the trajectory $z(t;t_{0},z_{0})$ and the time interval $[t_{0},t_{0}+T]$ are defined as

[TABLE]

Note that FTLEs are well-defined for any finite integration time $T$. Under certain assumptions, the limit $T\to\infty$ also exists [27].

2.2 Optimally Time-Dependent (OTD) modes

We give an overview of the OTD modes [2] for general dynamical systems. The OTD modes represent a reduced-order set of $r$ time-dependent orthonormal modes, $U(t)=[u_{1}\left(t\right),u_{2}\left(t\right),...,u_{r}\left(t\right)]$, that minimize the difference between the action of the infinitesimal propagator $\nabla F^{t+\delta t}_{t}$ on $u_{i}(t)$ and its value at time $t+\delta t$. More specifically, we seek to minimize the functional

[TABLE]

with the constraint that the time-dependent basis satisfies the orthonormality conditions

[TABLE]

Using Taylor series expansions, we have

[TABLE]

where $I$ is the identity matrix. Replacing the above equations into the functional (10) results in

[TABLE]

We emphasize that the minimization of the function (14) is considered only with respect to the time-derivative (rate of change) of the basis, $\dot{U}(t),$ instead of the basis $U(t)$ itself. This is because we do not want to optimize the subspace that the operator is acting on, but rather find an optimal set of vectors, $\dot{U}(t),$ that best approximates the linearized dynamics in the subspace $U$. We then solve the resulted equations and compute $U(t)$. We will refer to these modes as the optimally time-dependent (OTD) modes, and the space that these modes span as the OTD subspace. By utilizing the minimization principle and taking into account the orthonormality constraint we obtain the following theorem (proved in [2]).

Theorem 2.1**.**

A one parameter family of vectors $u_{i}(t)\in\mathbb{R}^{n}$, $i=1,2,\cdots,r$, minimizes the functional (14) and satisfies the orthonormality condition (11) if and only if

[TABLE]

where $U(t)=[u_{1}\left(t\right),u_{2}\left(t\right),...,u_{r}\left(t\right)]\in\mathbb{R}^{n\times r}$.

Equation (15) is the evolution equation for the OTD modes. It was proven in [2] that for the case of a time independent operator, $L$, the subspace spanned by the columns of $U(t)$ converges asymptotically to the modes associated with the most unstable directions of the operator $L$.

The OTD modes were also used to capture the transient non-normal growth of instabilities. In figure 1, we recall a system from [2] exhibiting a non-normal growth (cf. [2] for the details). A typical trajectory of the system, shown in the left panel, has an almost periodic behavior where each cycle exhibits three distinct regimes: (A) an exponential growth in the $z_{3}$ direction, (B) a non-normal growth in the $z_{1}-z_{2}$ plane with simultaneous decay in the $z_{3}$ direction and, (C) exponential decay in the $z_{1}-z_{2}$ plane to the origin. This configuration allows for the repeated occurrence of exponential (regime A) and non-normal (regime B) instabilities.

Due to the exponential instability close to the origin the system undergoes chaotic transitions between positive and negative values of $z_{3}$. We compute the instantaneous growth rate corresponding to the direction of a single OTD mode. The OTD mode initially captures the severe exponential growth and subsequently captures the non-normal growth. On the other hand the real part of the eigenvalues of the full linearized operator can only capture the exponential growth, even in regimes where it is not relevant, while they completely miss the non-normal growth. The maximum largest singular value $\sigma$ of the matrix $L$ exhibits a similar behavior.

The OTD modes share some fundamental characteristics of other stability measures. In the next section, we show the relation between the OTD modes and the finite-time Lyapunov vectors. In Appendix A, we also discuss the connection between them and the Dynamically Orthogonal (DO) modes introduced in [30]. Their connection with covariant Lyapunov vectors [16] remains to be explored.

3 Alignment of the OTD modes with the dominant Cauchy–Green strain eigenvectors

Here we prove that under proper assumptions, the OTD modes not only capture the unstable directions of steady linear operators, but also the transient instabilities associated with time-dependent operators. More specifically, we prove that the OTD subspace converges exponentially fast to the eigenspace associated with the dominant eigenvalues of the Cauchy–Green strain tensor. The latter describes the finite-time growth of infinitesimal perturbations in the phase space. For what follows we will need to measure the distance between two subspaces of the same dimensionality. To this end we need the following definitions.

Definition 3.1**.**

Two $r$-dimensional linear subspaces of $\mathbb{R}^{n}$ spanned by the columns of $U\in\mathbb{R}^{n\times r}$ and $V\in\mathbb{R}^{n\times r}$ are equivalent if there is an invertible matrix $R\in\mathbb{R}^{r\times r}$ such that $U=VR$.

The next definition provides a quantity for measuring the ‘angle’ between two subspaces.

Definition 3.2**.**

For two $r$-dimensional linear subspaces of $\mathbb{R}^{n}$ spanned by the columns of $U=[u_{1},u_{2},...,u_{r}]\in\mathbb{R}^{n\times r}$ and $V=[v_{1},v_{2},...,v_{r}]\in\mathbb{R}^{n\times r}$, we define the distance function

[TABLE]

where $\|\cdot\|$ denotes the Frobenius norm.

Note that the entries of $U^{T}V$ are the inner products between $u_{i}$’s and $v_{j}$’s, i.e., (U^{T}V)_{ij}=\big{<}u_{i},v_{j}\big{>}. The following lemma shows that the equivalence of two subspaces can be deduced from the distance $\gamma_{U,V}$.

Lemma 3.1**.**

Given a subspace defined by the orthonormal columns of $U\in\mathbb{R}^{n\times r}$, and a subspace defined by the unit-length (but not necessarily orthonormal) columns of $V\in\mathbb{R}^{n\times r}$, then we always have $\gamma_{U,V}\leq 1$. Moreover, the two subspaces spanned by the columns of $U$ and $V$ are equivalent if and only if $\gamma_{U,V}=1$.

Proof.

By the orthogonal projection theorem, there are $h_{i}\in\mathbb{R}^{r}$ and $b_{i}\in\mathbb{R}^{n}$ such that $v_{i}=Uh_{i}+b_{i}$ and $U^{T}b_{i}=0$ for $i=1,2,\cdots,r$. Note that since $v_{i}$ are unit length, we have

[TABLE]

Defining $H=[h_{1}|\cdots|h_{r}]$ and $B=[b_{1}|\cdots|b_{r}]$, we have $V=UH+B$. This implies $U^{T}V=H$ and hence

[TABLE]

Therefore, $\gamma_{U,V}^{2}=1-\|B\|^{2}/r$. Note that the subspaces $\mathrm{col}(U)$ and $\mathrm{col}(V)$ are equal if and only if $\|B\|=0$. Therefore, $\mathrm{col}(U)=\mathrm{col}(V)$ if and only if $\gamma_{U,V}=1$. We also note that if the two subspaces are not equivalent $\gamma_{U,V}<1.$ ∎

To prove the main theorem it is more convenient to work with the subspace spanned by the linearized flow (3), instead of the OTD subspace. The following result states that the two subspaces are equivalent.

Theorem 3.1**.**

Let $V(t)\in\mathbb{R}^{n\times r}$ solve the equation of variations (3) and $U(t)\in\mathbb{R}^{n\times r}$ with $U^{T}U=I$ solve the OTD equation (15). Assume that the two subspaces spanned by the columns of $U$ and $V$ are initially equivalent, i.e. there is an invertible matrix $T_{0}\in\mathbb{R}^{r\times r}$ such that $V_{0}=U_{0}T_{0}$. Then there exists a one-parameter family of linear transformations $T(t)$ such that $V(t)=U(t)T(t)$ for all $t$ and $T(t)$ satisfies the differential equation

[TABLE]

where $L_{r}=U^{T}LU$ is the orthogonal projection of $L$ to the linear subspace spanned by the columns of $U$.

Proof.

See Theorem 2.4 in [2]. ∎

Based on this equivalence, it is sufficient to study the behavior of the subspace evolved under the time-dependent linearized dynamics (3) in order to understand the properties of the OTD subspace. Specifically, we use this lemma to prove the main result of this paper, namely, that under a spectral gap condition, the OTD subspace will converge to the dominant directions of the Cauchy–Green strain tensor.

Theorem 3.2**.**

Let $\xi_{i}(t,t_{0},z_{0})$ and $\eta_{i}(t,t_{0},z_{0})$ $(i=1,2,\dots,n)$ denote the eigenvectors of the right and left Cauchy–Green strain tensors, respectively, with corresponding eigenvalues $\lambda_{1}(t,t_{0},z_{0})\geq\ \lambda_{2}(t,t_{0},z_{0})\geq\dots\geq\lambda_{n}(t,t_{0},z_{0})$. Moreover, let $U(t)\in\mathbb{R}^{n\times r}$ solve the OTD equation (15) along the trajectory $z(t;t_{0},z_{0})$. Assume that

i) The subspace $U(t)$ is initiated so that each basis element $u_{i}(t_{0})$ has a non-zero projection on at least one of the eigenvectors $\xi_{i}(t,t_{0},z_{0})$, $i=1,2,\cdots,r$ for all times $t$.

ii) The trajectory $z(t;t_{0},z_{0})$ is hyperbolic in the sense that there exist constants $a_{1}\geq a_{2}\geq\cdots\geq a_{r}>a_{r+1}\geq\cdots\geq a_{n}$ and $K_{i}>0$ $(i=1,2,\cdots,n)$ such that

[TABLE]

Then the subspace $U(t)$ aligns with the $r$ most dominant left Cauchy–Green strain eigenvectors, $\eta_{i}(t,t_{0},z_{0})$, $i=1,...,r$, exponentially fast as $t\to\infty$.

Proof.

Based on Lemma 3.1 the subspace spanned by the OTD modes is equivalent to the subspace that we obtain if we evolve $V(t)=[v_{1}(t),...,v_{r}(t)]$ using the linearized dynamics (3). We represent the initial condition for the equation of variations (3) as $V(t_{0})=[v_{0_{1}},v_{0_{2}},...,v_{0_{r}}]$ and express it in terms of the orthonormal strain eigenbasis $\{\xi_{i}(t,t_{0},z_{0})\}$ as

[TABLE]

Then, equations (4) and (7) applied to equation (19) yield

[TABLE]

We then have

[TABLE]

Moreover,

[TABLE]

Based on this last equation, the cosine of the angle $\alpha_{v_{j},\eta_{k}}(t,t_{0},z_{0})$ between the vector $v_{j}(t)$ and the eigenvector $\eta_{k}(t,t_{0},z_{0})$ of $B_{t_{0}}^{t}(z_{0})$ obeys the relation

[TABLE]

Identity (20) and definition 3.2 yield

[TABLE]

where

[TABLE]

Note that assumption (ii) implies $\lim_{t\to\infty}\lambda_{r+1}(t,t_{0},z_{0})/\lambda_{r}(t,t_{0},z_{0})=0$ since

[TABLE]

As a consequence, we obtain from the last inequality the estimate

[TABLE]

Note that \sum_{k=1}^{r}\big{<}v_{0_{j}},\xi_{k}\big{>}^{2}>0 from the first assumption. Therefore, for the asymptotic limit $t\to\infty$, we have

[TABLE]

Note that the way we have defined the two subspaces satisfy the assumptions of Lemma 3.1 and therefore we always have $\gamma_{v,\eta}^{2}\leq 1$. Thus, $\gamma_{v,\eta}^{2}=1$ which implies that the two subspaces are equivalent. This completes the proof. ∎

We point out that the hyperbolicity condition (ii) is not a necessary condition. In fact, as long as the ratio $\lambda_{r+1}(t)/\lambda_{r}(t)$ tends to zero asymptotically, the alignment takes place (see equation (22)).

4 Eigenvalue crossing and OTD modes

Theorem 3.2 shows that under appropriate assumptions the OTD modes converge exponentially fast to the eigenspace of the left Cauchy–Green strain tensor associated with the largest Lyapunov exponents. This convergence, however, may be interrupted when the smallest eigenvalue spanned by the OTD modes crosses with the one that is not being spanned. In figure 2, we illustrate this situation. The green curve denotes the trajectory of the system, while the ellipses indicate the principal axes of the left Cauchy–Green strain tensor. Assuming that we are resolving only one OTD mode, after sufficient time this must have converged to the correct principal direction (red arrow). As we go through the critical time instant where the ellipsoid becomes a circle due to eigenvalue crossing, we observe that the principal directions are instantaneously ill-defined. Immediately after the eigenvalue crossing, the dominant principal direction undergoes a $90$ degrees internal rotation compared to immediately before the crossing. The OTD modes, on the other hand, evolve smoothly and are unable to capture such an instantaneously unbounded transition.

In the context of shearless Lagrangian coherent structures, these eigenvalue crossings are referred to as the Cauchy–Green singularities and play an important role in the computation of jet cores [11]. In a more general application, the significance of these singularities is pointed out by Lancaster [22] who derived the rate of change of the eigenvectors of an arbitrary symmetric matrix, proving the following theorem.

Theorem 4.1**.**

The rate of change of the eigenvectors $R(t)\in\mathbb{R}^{n\times n}$ of a symmetric matrix $G(t)\in\mathbb{R}^{n\times n}$ is described by the equation

[TABLE]

where $K(t)$ is a skew-symmetric matrix given by

[TABLE]

with $\widetilde{G}(t)=R^{T}(t)\dot{G}(t)R(t)$ and $\lambda_{i}(t)$ being the eigenvalues of $G(t)$.

Proof.

See Appendix B or [22]. ∎

As $\lambda_{i}$ approaches $\lambda_{j}$, $\dot{R}$ becomes unbounded, consistent with the description given previously. This is not a problem if both of the eigen-directions corresponding to $\lambda_{i}$ and $\lambda_{j}$ are being spanned/resolved by the OTD subspace. The corresponding OTD subspace will just evolve smoothly in this case. However, if the eigenvalue crossing occur with an eigendirection that is not being spanned by the OTD subspace, the latter will have to instantaneously ‘steer’ towards a direction that is orthogonal to it.

This is not possible since it can be easily seen that for any bounded operator $L,$ the rate of change of the corresponding OTD modes (eq. (15)) is bounded:

[TABLE]

where $Q$ is a projection operator (therefore bounded) and $U$ is also bounded. If the eigenvalue crossing involves modes already contained in $U$ then there is no need for infinite rate of change of the modes, because existing modes just exchange roles (or indices). However, if there is eigenvalue crossing between a mode that is contained in $U$ and another one that is not contained then it is essential to have an evolution of the direction in $U$ to an orthogonal direction (not contained in $U$) and this cannot happen instantaneously.

Therefore, we see that the convergence of the OTD modes can be compromised if there are eigenvalue crossings involving directions not spanned by the OTD subspace. Note however that such eigenvalue crossing will only influence the convergence of the OTD mode associated with the smallest finite-time Lyapunov exponent captured within the OTD subspace. Assuming that there are no subsequent eigenvalue crossings with larger eigenvalues, one expects that the OTD subspace would capture the directions associated with the largest finite-time Lyapunov exponents.

Thus, for systems where eigenvalue crossing occurs, a larger OTD subspace guarantees that the principal axes associated with the largest Lyapunov exponents will be captured by the OTD modes. Note that such convergence does not depend on the physical dimension of the system or the intrinsic dimension of the underlying attractor. This is particularly important for systems where instabilities of isolated modes emerge out of high-dimensional stochastic backgrounds, such as intermittent dissipation bursts in turbulence [15] or focusing instabilities occurring in nonlinear stochastic waves [8, 15].

As we demonstrate in the next section, for many practical applications it is sufficient to compute just the maximum finite-time Lyapunov exponent so what is important is to have the most dominant direction of the Cauchy–Green tensor captured by the OTD subspace. The benefits of such approach are, of course, more important as the dimension $n$ of the underlying system increases.

5 Reduced-order computation of finite-time Lyapunov exponents

As shown above, the OTD modes are able to extract the transiently most unstable directions independently of the attractor dimensionality. Consequently, the OTD basis can be used for the reduced-order approximation of finite-time Lyapunov exponents (FLTE) of high-dimensional systems. Specifically, we apply the following approximation scheme for the computation of a FTLE corresponding to a time interval of length $T$.

Advect the trajectory $z(t;t_{0},z_{0})$ for an interval $t\in[t_{0},t_{0}+T]$ where $z_{0}$ is the initial point. 2. 2.

Compute the $r$-dimensional OTD subspace $U(t)$ corresponding to this trajectory. 3. 3.

Compute the low-dimensional fundamental solution matrix $\Phi_{t_{0}}^{t}\in\mathbb{R}^{r\times r}$ using the reduced linear dynamical system

[TABLE]

where $L_{r}(z,t)=U(t)^{T}L(z,t)U(t)$ is the projection of the full linearized operator $L$ onto the OTD subspace, which contains the most unstable directions. 4. 4.

Compute the reduced-order finite-time Lyapunov exponents,

[TABLE]

where $\gamma_{i}$ denotes the eigenvalues of the reduced-order right Cauchy–Green strain tensor $\left(\Phi_{t_{0}}^{t_{0}+T}\right)^{\top}\Phi_{t_{0}}^{t_{0}+T}$, ordered such that $\gamma_{1}\geq\gamma_{2}\geq\cdots\geq\gamma_{n}$.

The described algorithm allows for the computation of FTLE by solving for each initial point, $z_{0}$, one $n$-dimensional equation, as well as additional $r$ equations which are $n$-dimensional for the OTD subspace. This cost has to be summed together with the solution of reduced-order linear system. To this end we have the total cost for the reduced-order and full algorithms:

• •

Full-order calculation cost: $(n+1)\times n$ consisting of $n+1$ equations of dimension $n$ as follows:

[TABLE]

• •

Reduced-order calculation cost: $n\times(r+1)+r^{2}$ consisting of $r+1$ equations of dimension $n$ and $r$ equations of dimension $r$ as follows:

[TABLE]

We observe that the cost of computing FTLE using the full system is $\mathcal{O}(n^{2})$. The computational cost of the reduced-order FTLE, however, is $\mathcal{O}(n)$ for fixed $r$. In principle, the dimension $r$ of the reduced system can be as low as $1$. But, as we will see in what follows, a one-dimensional OTD subspace may lead to false troughs in the FTLE field. A higher dimensional OTD subspace returns more accurate estimation of the FTLE field and often circumvents the false troughs.

Figure 3 shows the computational cost of the full FTLE calculation versus the reduced-order FTLE calculation. Note that for low-dimensional systems (small $n$) the reduction in the computational cost is not significant. However, as the system dimension increases ($n\gg r$), the computational cost of the reduced-order FTLE can be orders of magnitude lower than the full FTLE computations.

5.1 The ABC Flow

As the first example we consider the ABC (Arnold–Beltrami–Childress) flow with three dimensional velocity field

[TABLE]

which is an exact solution of the Euler equation for the ideal incompressible fluids [1]. The ABC flow has served as a prototype example for testing numerical methods for computing Lagrangian coherent structures [17, 13, 14]. Here, we set $A=\sqrt{3}$, $B=\sqrt{2}$ and $C=1$ which allows for the existence of chaotic Lagrangian trajectories [17].

We compute the FTLE on a $251\times 251$ grid of initial conditions on each face of the cube $[0,2\pi]^{3}$ with an integration time of length $T=8$. We use finite differences to approximate the deformation gradient $\nabla F_{t_{0}}^{t_{0}+T}$. To increase the finite difference accuracy, we use six auxiliary grids $(z_{0,1}\pm h,z_{0,2}\pm h,z_{0,3}\pm h)$ around each grid point $z_{0}=(z_{0,1},z_{0,2},z_{0,3})$. The parameter $h$ controls the accuracy of the finite differences and is set to $h=10^{-8}$ in the following. We refer to [23, 12] for further details on the approximation of the deformation gradient. Once the deformation gradient is approximated, computing the Cauchy–Green strain tensor and consequently the FTLE is straightforward. Note that the deformation gradient can also be obtained by numerically integrating the equation of variations (3).

To compute the reduced-order FTLE, we numerically integrate system (29c) from the initial grid points $z_{0}$ with an integration time of length $T$. The initial condition $U(0)$ for the OTD equation (29b) is the $r$ most dominant eigenvectors of the symmetric linear operator at $t=0$, i.e. $L^{s}=(L(z_{0},0)+L(z_{0},0)^{T})/2$. This choice of the OTD initial condition is made since these eigenvectors are the instantaneously most unstable directions [19].

Figure 4 shows the reduced-order FTLE fields using OTD reduction of sizes $r=1$ (left column) and $r=2$ (middle column). The true FTLE field (right column) is also shown for comparison. Each row of the figure shows a different cross section of the field. For $r=2$, the reduced-order FTLE is close to the true FTLE field. This is also true qualitatively for $r=1$. However, for $r=1$, the reduced FTLE exhibits additional troughs that do not exist in the true FTLE field.

We demonstrate in figure 5 that these false troughs occur when there is a crossing (or near crossing) between the first and second most dominant eigenvalues of the right Cauchy–Green tensor, i.e., $\lambda_{1}$ and $\lambda_{2}$. Recall that this eigenvalue crossing violates condition (ii) of Theorem 3.2 and hence the discrepancy between the reduced-order FTLE and the true FTLE is expected. The false troughs are eliminated as we increase the dimension of the OTD reduction from $r=1$ to $r=2$. This is due to the fact that there is not eigenvalue crossing between $\lambda_{2}$ and $\lambda_{3}$.

The above observation is demonstrated in figure 5 where the reduced-order FTLEs at two select points are compared against the true FTLE as the integration time $T$ varies. In figure 5(a), the initial point $z_{0}=(2.0,0.6,0.0)$ is chosen to be at one of the false troughs. An eigenvalue crossing occurs at approximately $T=0.2$. For advection times before $T=0.2$, both one- and two-dimensional reductions follow the largest FTLE very closely. However, near the eigenvalue crossing the most dominant eigenvector of the Cauchy–Green strain tensor undergoes rotation with an unbounded rate. As demonstrated in inequality 25, the one-dimensional OTD reduction cannot follow such rapid rotation. This is confirmed in figure 5(a) where it can be seen that the FTLE obtained from the one-dimensional OTD reduction, after $T=0.2$, follows the second most dominant FTLE instead of undergoing a rapid change that is required to follow the most dominant FTLE. Near $T=1$, the FTLE of the one-dimensional reduction departs from non-dominant FTLEs and starts increasing. This behavior is to be expected since in the remaining advection time no eignevalue crossing occurs and as a result the one-dimensional reduction continues to recover to converge to the most dominant eigenvector – however slowly.

The two-dimensional OTD reduction, however, recovers much faster compared to the one-dimensional OTD reduction after $T=0.2$. We observe that the two eigenvalues of OTD ($r=2$) approach each other nearly at $T=0.5$ causing in a relatively rapid repulsion of these two eigen-directions. This repulsion results in a rapid rotation of the dominant OTD vector towards the dominant eigen-direction of the Cauchy–Green tensor. This mechanism of rapid correction of the OTD vectors is entirely absent in the one-dimensional OTD reduction.

Figure 5(b) shows the evolution of the same quantities along a different ABC trajectory. Here, eigenvalue crossings are absent. As a result both one-dimensional and two-dimensional OTD reductions closely follow the most dominant FTLE and in the case of $r=2$ both most and second most dominant FTLEs. This demonstrates that, in the absence of eigenvalue crossing, one-dimensional OTD reduction accurately approximates the most dominant FTLE.

5.2 Charney–DeVore model with regime transitions

In this section, we apply the OTD reduction to the six-dimensional truncation of the equations for barotropic flow in a plane channel with orography that are known as the Charney–DeVore (CDV) model. The truncated system is given by

[TABLE]

The model coefficients are given by

[TABLE]

Following [10], we set $(z_{1}^{*},z_{4}^{*},C,\beta,\gamma,b)=(0.95,-0.76095,0.1,1.25,0.2,0.5)$. For these parameters the CDV model generates regime transitions due to the interaction of barotropic and topographic instabilities. It was shown in [9] that a Proper Orthogonal Decomposition (POD) reduction of this model to three leading POD modes resolves 97% of cumulative variance but cannot capture the chaotic regime transitions present in the six-dimensional model. These highly transient instabilities render this model an appropriate test case for evaluating the performance of the OTD reduction in computing the FTLE. Moreover, using this model we asses the performance of the OTD reduction for a higher dimensional problem.

We compute the true FTLE using the finite difference method analogous to the ABC flow problem. Similarly, the initial condition $U(0)$ for the OTD equation (29b) are the $r$ most dominant eigenvectors of the symmetric linearized operator. The fourth-order Runge-Kutta scheme with time step size $\Delta t=0.4$ (days) is used for the numerical integration of equations (28b) and (29c). Using the smaller time step $\Delta t=0.2$ (days) did not change the numerical results.

Figure 6 shows the reduced FTLE fields obtained from the one-dimensional reduction (left column) and two-dimensional reduction (middle column) as well as the full FTLE (right column). Each row shows a two-dimensional cross section of the phase space. As in the ABC flow, the reduced FTLE fields agree qualitatively with the fulll FTLE field. However, for $r=1$ reduction, some false troughs are observed. Away from the false troughs, the one-dimensional reduction estimates the FTLE accurately. For $r=2$ reduction, there are no false troughs.

In figure 7, the FTLE values are plotted as a function of the integration time. These correspond to the point $z_{0}=(1.14,0,0,-0.91,0,0)$, marked by a cross symbol in the inset, which lies on a false trough in the $z_{1}-z_{4}$ section. An eigenvalue crossing occurs at approximately $T=9$. Analogues to the ABC flow, we observe that before $T=9$ both one- and two-dimensional OTD reductions capture the most dominant eigenvalues of the Cauchy–Green strain tensor. After $T=9$, the one-dimensional OTD reduction continuously follows the second most dominant eigen-direction of the Cauchy–Green strain tensor. The inability of a single OTD vector to undergo a dramatic rotation to follow the most dominant eigen-direction of the Cauchy–Green tensor leads to the false trough that is observed in Figure 6(a). The two-dimensional OTD subspace, however, converges to the subspace spanned by the two most dominant eigenvectors of the Cauchy–Green strain tensor. This convergence is guaranteed by Theorem 3.2 since the second and the third Cauchy–Green eigenvalues $(\lambda_{2},\lambda_{3})$ do not cross (see figure 7). As a result, the two-dimensional ($r=2$) OTD reduction closely approximates the two larges FTLEs for all times.

The above results demonstrate that, as long as $\lambda_{2}$ and $\lambda_{3}$ do not coincide, a reduced FTLE with two OTD modes reliably approximates the dominant FTLEs of the full system regardless of the system’s dimension $n$. This amounts to a significant reduction in the computational cost of the FTLE evaluations, and the gain in the computational speed is larger as the dimension of the dynamical system increases. Investigating the numerical performance of the OTD reduction for infinite-dimensional systems is the subject of future study.

6 Conclusions

We have examined the properties of the optimally time-dependent (OTD) modes for general time-dependent dynamical systems. Specifically, we have shown that under mild conditions, related to the spectrum of the Cauchy–Green tensor, the OTD modes converge exponentially fast to the eigendirections of the Cauchy–Green tensor associated with the dominant eigenvalues (largest finite-time Lyapunov exponents). Therefore, the OTD modes can be employed to extract time-dependent subspaces that encode information associated with transient dynamics (finite-time instabilities).

We have applied the derived result on the formulation of a reduced order algorithm for the computation of the maximum finite-time Lyapunov exponent, a measure that has been used for the quantification of Lagrangian Coherent Structures. We demonstrated the derived results through two specific examples, the three-dimensional ABC flow and the six-dimensional Charney-DeVore model. In both case we thoroughly analyzed the limitations due to dimensionality reduction in combination with eigenvalue crossing. Apart of its value as a computational method for finite-time Lyapunov exponents the presented result paves the way for the development of efficient control and prediction strategies for chaotic dynamical systems exhibiting transient features, such as extreme events and off-equilibrium dynamics.

Acknowledgments

T.P.S. has been supported through the ARO grant 66710-EG-YIP, the AFOSR grant FA9550-16-1-0231, the ONR grant N00014-15-1-2381 and the DARPA grant no. HR0011-14-1-0060. H.B. and M.F. have been supported through the first, second and fourth grants as postdoctoral associates.

Appendix A: The OTD modes and the dynamically orthogonal modes

The set of OTD equations have the same form as the Dynamically Orthogonal (DO) field equations [28, 30]. In fact, because the minimization of the function $\mathcal{F}$ is performed only over the rate of change of the basis elements, $\dot{u}_{i}(t)$, (and not the basis elements $u_{i}(t)$) one can follow exactly the same steps as in Theorem 2.1 in [2] to show that the evolution equations (15) can be derived for the general case of a nonlinear operator $\mathcal{N}(U)$. More specifically, by minimization of the functional:

[TABLE]

we can obtain the evolution equations.

[TABLE]

For the case of the DO equations the nonlinear operator takes the form

[TABLE]

where we have used the notation in [28], i.e. $\mathcal{L}$ is the right hand side of the stochastic PDE, $\bar{u}$ is the trajectory of the mean, $Y_{i}(t)$ are the stochastic coefficients and $C_{Y_{i}Y_{j}}(t)$ is their covariance matrix. The operation $E^{\omega}$ denotes the average with respect to an appropriate probability measure. Conversely, one can obtain the OTD mode equations from the DO equations simply by considering their deterministic limit, i.e. by taking the limit $C_{Y_{i}Y_{j}}(t)\rightarrow 0$.

Appendix B - Proof of Theorem 4.1

We consider the rate of change of the eigenvectors of a general one-parameter family of symmetric matrices, $G(t)\in\mathbb{R}^{n\times n}$. For such a matrix we have the eigenvalue problem

[TABLE]

where the columns of $R(t)\in\mathbb{R}^{n\times n}$ are the eigenvectors of $G(t)$ and $\Lambda(t)\in\mathbb{R}^{n\times n}$ is the diagonal matrix of corresponding eigenvalues. Note that, since $G$ is symmetric, its eigenvectors are orthogonal and therefore $RR^{T}=I$. Differentiating with respect to time, we obtain

[TABLE]

Multiplying the above equation by $R^{T}$ from the left yields

[TABLE]

Since $R^{T}R=I$, we have:

[TABLE]

Let $K=R^{T}\dot{R}$. From the above relation we observe that $K$ is a skew-symmetric matrix. Denoting $\widetilde{G}=R^{T}\dot{G}R$, we have

[TABLE]

For off-diagonal terms $i\neq j$, we have $\dot{\Lambda}_{ij}=0$. Therefore,

[TABLE]

Transposing the above relation

[TABLE]

Using the skew-symmetric property of $K$, i.e. $K_{ij}=-K_{ji}$, and summing up the above two equations yields

[TABLE]

This results in the following closed-form evolution equations for $\Lambda$ and $R$,

[TABLE]