An ergodic averaging method to differentiate covariant Lyapunov vectors
NNoname manuscript No. (will be inserted by the editor)
An ergodic averaging method to differentiatecovariant Lyapunov vectors
Computing the curvature of one-dimensional unstablemanifolds of strange attractors
Nisha Chandramoorthy · Qiqi Wang
Submitted version
Abstract
Covariant Lyapunov vectors or CLVs span the expanding and contract-ing directions of perturbations along trajectories in a chaotic dynamical system.Due to efficient algorithms to compute them that only utilize trajectory infor-mation, they have been widely applied across scientific disciplines, principallyfor sensitivity analysis and predictions under uncertainty. In this paper, we de-velop a numerical method to compute the directional derivatives of CLVs alongtheir own directions. Similar to the computation of CLVs, the present methodfor their derivatives is iterative and analogously uses the second-order derivativeof the chaotic map along trajectories, in addition to the Jacobian. We validatethe new method on a super-contracting Smale-Williams Solenoid attractor. Wealso demonstrate the algorithm on several other examples including smoothly per-turbed Arnold Cat maps, and the Lorenz’63 attractor, obtaining visualizations ofthe curvature of each attractor. Furthermore, we reveal a fundamental connectionof the CLV self-derivatives with a statistical linear response formula.
Keywords chaotic dynamics · Lyapunov vectors · uniform hyperbolicity Covariant Lyapunov Vectors (CLVs) [17] are specific bases for tangent spaces alonga trajectory, characterized by Lyapunov exponents. Ginelli e t al.’s [11] efficient al-gorithm to compute CLVs has led to several applications of Lyapunov analysis in This work was supported by Air Force Office of Scientific Research Grant No. FA8650-19-C-2207.N. ChandramoorthyMechanical Engineering and Center for Computational Science and Engineering77 Massachusetts Avenue, Cambridge MA E-mail: [email protected]. WangDepartment of Aeronautics and Astronautics and Center for Computational Science and En-gineering77 Massachusetts Avenue, Cambridge MAE-mail: [email protected] a r X i v : . [ n li n . C D ] J u l Nisha Chandramoorthy, Qiqi Wang engineering, in both deterministic and stochastic chaotic systems. These applica-tions include uncertainty quantification, data assimilation and forecasting, acrossa range of disciplines such as numerical weather prediction and aerospace engineer-ing ([5], [21], [23], [15]; see [6] for a survey of applications of Lyapunov analysis). Inthis work, we propose a new method to computationally estimate the directionalderivatives of CLVs along their own respective directions. We shall refer to thesein this paper as
CLV self-derivatives .In the case of a one-dimensional unstable manifold, the CLV corresponding tothe largest Lyapunov exponent is the unit tangent vector field along the unstablemanifold. The norm of this CLV self-derivative is hence also the curvature of theunstable manifold. Thus, a primary contribution of this work is the computation ofthe curvature of one-dimensional unstable manifolds. Furthermore, we draw uponthe connection between statistical linear response, or the change in the invariantprobability distribution on the attractor due to system parameter perturbations,and the CLV self-derivatives, to reveal the latter’s fundamental nature. This con-nection is not limited to one-dimensional unstable manifolds. In particular, asshown in previous work [9], a tractable computation of the derivative of the in-variant probability measure with respect to perturbations along an unstable CLV,involves the computation of the CLV self-derivatives. The CLV self-derivatives area key ingredient that will enable sensitivity computation in chaotic systems, whichis an active area of research (see [27], [1], [20] for some of the recently developedsensitivity computation methods for chaotic systems).The outline of the subsequent sections is as follows. In section 2, we brieflysummarize the theory of covariant Lyapunov vectors and establish the settingwe derive our results in: uniformly hyperbolic chaos. The numerical method tocompute CLV self-derivatives, henceforth known as differential CLV method isderived in section 3; while the main steps are in section 3.3, notational setup andthe intuition for the steps are developed in the prior subsections. We validate themethod using a super-contracting Solenoid map in section 4.1. Further numericalexperiments demonstrating the method on the Lorenz’63 attractor, a volume-preserving perturbed Cat map, a dissipative perturbed Cat map, and the H´enonmap are in sections 4.2, 4.3, 4.4 and 4.5 respectively. The implication of the methodfor the computation of linear response is discussed in section 5. We summarize ourresults and conclude in section 6.
The dynamical system studied in this paper is the iterative application of a smooth( ∈ C k , for k ≥
3) self-map ϕ : M → M of a domain M . We write ϕ n to denotean n -time composition of ϕ . For simplicity, we take throughout, the domain M tobe a compact subset of R d . The iterates under ϕ , or the points along trajectoriesof the dynamical system, are represented using the following subscript notation: if p ∈ M , p n := ϕ n p ; p is simply written as p , which we use to denote an arbitraryphase point. A similar notation is also adopted for scalar or vector-valued functionsor observables. If f is an observable, f n := f ◦ ϕ n . The derivative with respectto the state is denoted as D and the partial derivative operators, with respectto the Euclidean coordinate functions x , x , · · · , x d are written as ∂ , ∂ , · · · , ∂ d , n ergodic averaging method to differentiate covariant Lyapunov vectors 3 respectively. For instance, if f : M → R is a scalar-valued observable, the derivative Df evaluated at p is given by Df ( p ) = [ ∂ f ( p ) , · · · , ∂ d f ( p )] T . Using the notationintroduced, an application of the chain rule would be as follows:( Df n ) T = (( Df ) n ) T Dϕ n . Finally, we assume the existence of an ergodic, physical, invariant measure for ϕ ,known as the SRB measure, and denote it as µ . As a result, ergodic (Birkhoff)averages of observables in L ( µ ) converge to their expectations with respect to µ :lim N →∞ (1 /N ) (cid:80) N − n =0 f n ( p ) = (cid:104) f, µ (cid:105) for Lebesgue-a.e. p. Note that such a measureis guaranteed to exist [28] in the uniformly hyperbolic setting, which we discuss insection 2.2.2.1 Tangent dynamicsIn order to introduce covariant Lyapunov vectors (CLVs), whose derivatives arethe subject of this paper, we briefly discuss the asymptotic behavior of tangentdynamics in chaotic systems. We refer to as tangent dynamics the linear evolutionof perturbations under the Jacobian matrix, Dϕ . Denoting the tangent space at p as T p M , Dϕ n ( p ) is a map from T p M to T p n M . Given a tangent vector v ∈ T p M ,we denote its iterate under the tangent dynamics at time n as v n ∈ T p n M . Thatis, v n = Dϕ n ( p ) v . Intuitively, if a perturbation of norm O ( (cid:15) ) is applied to p along v , up to first order in (cid:15) , the deviation from the original trajectory, after time n ,is along v n . In other words, v n = lim (cid:15) → ϕ n ( p + (cid:15)v ) − p n (cid:15) = Dϕ n ( p ) v . (1)In practice, the above equation for the tangent dynamics is solved iteratively, sinceusing the chain rule, Dϕ n ( p ) = Dϕ ( p n − ) · · · Dϕ ( p ) , and hence v n +1 = Dϕ ( p n ) v n . A classical result in nonlinear dynamics, known as the Oseledets multiplicativeergodic theorem (OMET) [4] deals with the asymptotic behavior of v n as n → ∞ ,in ergodic systems. The OMET implies the following: at µ -a.e. p ∈ M , the tangentspace splits as a direct sum of Dϕ -invariant subspaces, as T p M = ⊕ ai =1 E i ( p ), with a ≤ d . This splitting is based on the asymptotic, exponential growth/decay ratesof tangent dynamics in the subspaces E i ( p ). More precisely, if v i ∈ E i ( p ), its normunder the tangent dynamics grows/decays exponentially at a rate that convergesto a constant. The limits λ i := lim n →∞ log (cid:13)(cid:13)(cid:13)(cid:13) v in v i (cid:13)(cid:13)(cid:13)(cid:13) , (2)1 ≤ i ≤ a , are known as the Lyapunov exponents (LEs). In the case of invert-ible, ergodic, chaotic maps, which is our setting, a) the LEs are almost every-where independent of p , and b) there is at least one positive LE. Arranging the LEs in descending order as λ ≥ λ ≥ · · · ≥ λ a , λ > , in a chaotic sys-tem. Let d u be the number of positive LEs, and d s = a − d u be the numberof negative LEs. Then, E u ( p ) := ⊕ d u i =1 E i ( p ) is called the unstable subspace of T p M . In other words, the unstable subspace E u ( p ) is the set of tangent vectorsthat asymptotically grow exponentially in norm under the tangent dynamics; by Nisha Chandramoorthy, Qiqi Wang definition, the unstable subspaces at points on a chaotic orbit are non-empty.Similarly, the set of tangent vectors that asymptotically decay exponentially innorm under the tangent dynamics, make up the stable subspace, denoted using E s ( p ) := ⊕ ai = d u +1 E i ( p ) = T p M \ E u ( p ). If each E i is one-dimensional and a = d ,the covariant Lyapunov vectors or CLVs, denoted as V i in this paper, are unitvector fields along E i . That is, CLVs satisfy the following properties: – The covariance property: Dϕ ( p ) V i ( p ) ∈ E i ( p ) . (3)Since by definition V i ( p ) is a unit vector, we introduce a scalar function z i : M → R + defined as z i ( p ) = (cid:13)(cid:13) Dϕ ( p ) V i ( p ) (cid:13)(cid:13) , to indicate the local stretching orcontraction factor of the i th CLV. Hence, the covariance property of the i thCLV can be expressed as Dϕ ( p ) V i ( p ) = z i V i ( p ) . (4) – The i th CLV grows/decays asymptotically on an exponential scale, at the rate λ i , and, in addition, is invariant under time-reversal: λ i := lim n →±∞ n log (cid:13) (cid:13)(cid:13) Dϕ n ( p ) V i ( p ) (cid:13) (cid:13)(cid:13) . (5) c > λ ∈ (0 ,
1) such that, at every point p ∈ M , i) every stable tangent vector v ∈ E s ( p ) satisfies: (cid:107) Dϕ n ( p ) v (cid:107) ≤ c λ n (cid:107) v (cid:107) , and ii) every unstable tangent vector v ∈ E u ( p ) satisfies: (cid:107) Dϕ − n ( p ) v (cid:107) ≤ c λ n (cid:107) v (cid:107) , for all n ∈ N . As a result, in thesesystems, there exist upper (lower) bounds that are independent of the base point p , on the slowest stretching (contracting) factors among z i ( p ). In particular, defin-ing C := cµ , we have z i ( p ) ≥ (1 /C ), 1 ≤ i ≤ d u and z i ( p ) ≤ C , d u + 1 ≤ i ≤ d. From the definition of the LEs (Eq. 5), it is also clear that they are the ergodic(Birkhoff) averages of the stretching/contraction factors: (cid:104) z i , µ (cid:105) := lim N →∞ N N − (cid:88) n =0 log z i ( p n ) = λ i , p ∈ M µ − a . e . (6)2.3 Examples A simple example of a uniformly hyperbolic system is Arnold’s Cat map, a smoothself-map of the surface of the torus ( T ≡ R / Z ): ϕ ([ x , x ] T ) = (cid:20) (cid:21) (cid:20) x x (cid:21) mod 1 . (7) n ergodic averaging method to differentiate covariant Lyapunov vectors 5 This is a linear hyperbolic system, i.e., the Jacobian matrix of the map is a constantin phase space and has eigenvalues other than 1. In this simple example, the CLVsand the stretching/contracting factors, are also independent of the phase point.The logarithm of the eigenvalues of the constant Jacobian matrix, are the LEs ofthis map: λ = log | (3 + √ / | and λ = log | (3 − √ / | . It is also clear that E = E u and E = E s are one-dimensional subspaces spanned by V and V ,the eigenvectors of the Jacobian matrix at eigenvalues of e λ and e λ respectively.Moreover, z and z are also constant on R / Z : z = e λ , and z = e λ . Further,the SRB measure for this map is the Lebesgue measure on R / Z . Since the Jacobian matrix is symmetric, the CLVs V and V are everywhereorthogonal to each other, but it is worth noting that this is a special case. Ina generic uniformly hyperbolic system, it is only true that the angle between theCLVs is uniformly bounded away from zero. The perturbed Cat maps treated laterhave additive perturbations to the Cat map above, that are smooth functions onthe torus. Two types of smooth perturbations are considered later, both designedto produce non-uniform behavior of the CLVs. Both perturbed Cat maps are stilluniformly hyperbolic, and differ in whether or not the resulting maps are area-preserving, in order to represent the two distinct cases of conservative (symplectic)and dissipative chaos.2.4 Lack of differentiability of E u and E s On hyperbolic sets, it is known that E u and E s are H¨older-continuous functionsof phase space, in a sense clarified in the Appendix section A. When the H¨olderexponent α , from Appendix section A, equals 1, we have Lipschitz continuity, butthis is indeed rare. Several examples (see [13] and references therein) have beenconstructed in which α is made to be arbitrarily small at almost all phase points,even in C ∞ maps. In rare cases, E u and E s are continuously differentiable whena certain bunching condition ([13], or section 19.1 of [16]) is satisfied by the LEs.Revisiting the examples, the perturbed Cat maps discussed above belong to therare category of maps with continuously differentiable stable/unstable subspaces.In fact, it can be shown that all uniformly hyperbolic maps on compact sets ofdimension 2, belong to this category (see Corollary 19.1.11 of [16]). While it wouldbe typical of a higher-dimensional map, even when uniformly hyperbolic, to shownon-smoothness of the stable and unstable subspaces, we have chosen to work withtwo-dimensional examples in this paper for easy visualization of the subspaces,which are lines in these maps.2.5 Derivatives of CLVs in their own directionsWhile the CLVs may lack differentiability on M , they have directional derivatives,in their own directions. In fact, it can be shown that these directional deriva- tives, which we refer to here as CLV self-derivatives, are themselves H¨older con-tinuous with the same exponent α (see Remark in the proof of Theorem 19.1.6of [16]). To wit, in two-dimensional uniformly hyperbolic systems, examples ofwhich are considered in this paper, both partial derivatives (along coordinate di-rections) of the CLVs exist, and hence the CLVs have directional derivatives in all Nisha Chandramoorthy, Qiqi Wang directions. The purpose of this paper, however, is to numerically compute direc-tional derivatives of CLVs along their respective directions in a general uniformlyhyperbolic system, regardless of their differentiability in phase space. Thus, wecompute the CLV self-derivatives, without using the partial derivatives along co-ordinate directions, which may not exist. The CLV self-derivatives are denoted by W i ( p ) ∈ T p T p M ≡ R d . They are defined using curves C p,i : [ − (cid:15) p , (cid:15) p ] → M withthe properties: i) C p,i (0) = p , ii) ( d C p,i /dt )( t ) = V i ( C p,i ( t )) , ∀ t ∈ [ − (cid:15) p , (cid:15) p ], as W i ( p ) := ∂ V i V i ( p ) := lim t → V i ( C p,i ( t )) − V i ( p ) t . (8)For example, in the case of a 1-dimensional unstable manifold, the curve C p, ,coincides with a local unstable manifold at p. Further discussion on the definitionof W i based on these curves, is postponed until section 3.1. Here we explain theexistence of these curves. The vector fields V i , 1 ≤ i ≤ d u are infinitely smoothon an open set in a local unstable manifold, and likewise, V i , for d u + 1 ≤ i ≤ d are infinite smooth on an open set in a local stable manifold. As a result, due tothe existence and uniqueness theorem, the flow of vector field V i , denoted by thecurve C · ,i exists and is uniquely defined, for some (cid:15) · > , justifying the definitionin Eq. 8. For simplicity, given T p M ≡ R d , we write all the tangent vectors in Euclidean coordinates. The output of the numerical method to be developed, W i , are d -dimensional vector fields consisting of component-wise directional derivatives ( ∂ V i )of V i . W i , being self-derivatives of CLVs, are naturally defined along trajectories, just like the CLVs.Thus, we seek a trajectory-based iterative procedure, in order to compute them. Weassume as input to the method the map, its Jacobian and second-order derivative,all computed along a long, µ -typical trajectory. The CLVs that need to be differ-entiated are also assumed as input, along the trajectory. To compute the CLVs,a standard algorithm such as Ginelli e t al.’s algorithm [11] can be used. This isan iterative procedure involving repeated QR factorizations of nearby subspacesto the one that is spanned by the required CLVs. For Ginelli e t al.’s algorithm,the reader is referred to [11] and [22] for its convergence with respect to trajectorylength, and for other algorithms that involve LU factorizations instead of QR, to[17].Besides using the computed CLVs as input, the differential CLV method wedevelop here for W i does not follow Ginelli e t al.’s or other algorithms for thecomputation of CLVs, primarily because the vector fields W i do not satisfy the covariance property. But the method resembles the latter algorithms in beingiterative and trajectory-based. One advantage of trajectory-based computation isthat we exploit for fast convergence (this aspect again being similar to the CLVcomputation algorithms) the hyperbolic splitting of the tangent space. This willbe clear at the end of the next section in which we give a step-by-step derivation. n ergodic averaging method to differentiate covariant Lyapunov vectors 7 In this section, we derive a numerical method to determine the quantity of interest, W i , which is defined in Eq. 8. In particular, fixing a reference trajectory p, p , · · · , we develop an iterative scheme that converges asymptotically to vectors W in := W i ( p n ), starting from an arbitrary guess for W i := W i ( p ) ∈ R d . The derivationresults in the following iteration, valid for 1 ≤ i ≤ d u , n ∈ Z + : W in +1 = (cid:16) I − V in +1 ( V in +1 ) T (cid:17) D ϕ ( p n ) : V in V in + Dϕ ( p n ) W in ( z in ) . (9)The iteration mainly uses the chain rule and the covariance property of V i , ina convenient set of coordinate systems centered along each µ -typical trajectory.These trajectory-based coordinates help us uncover each term on the right handside of Eq. 9.3.1 Change of coordinates and associated notationFix a µ -typical point p ∈ M and consider again the curves C p,i , 1 ≤ i ≤ d, whichwere introduced to define W i in Eq. 8. To reiterate, the curves C p,i : [ − (cid:15) p , (cid:15) p ] → M are such that i) C p,i (0) = p and ii) ( C p,i ) (cid:48) ( t ) = V i ( C p,i ( t )) , for all t ∈ [ − (cid:15) p , (cid:15) p ] . There exists a measurable function p → (cid:15) p that defines the extent of the curves sothat such a coordinate change, from [ − (cid:15) p , (cid:15) p ] d to a neighborhood of p , exists andis additionally continuous at each p ∈ M . This partly follows from an assertionproved in standard stable-unstable manifold theory: a closed (cid:15) p Euclidean ballaround the origin in R d u ( R d s ) has an embedding into a local unstable (stable)manifold at p . These pointwise coordinate systems are referred to as Lyapunovcharts or adapted coordinates in the theoretical literature ([16] Ch. 6, [18]).In writing Eq. 8, we made a particular choice of adapted coordinates. Wechose adapted coordinates that are adapted specifically to the CLVs, as opposedto any other basis of T p M , in the following sense. At each p , the image of the i thEuclidean basis vector e i under the differential of the coordinate change, is V i .More intuitively, we have chosen adapted coordinates such that the i th Euclideancoordinate, corresponds, under these coordinate changes, to points that are per-turbations along V i , at the phase point corresponding to the origin. Thus, ourquantity of interest, can be written by definition of CLV-adapted coordinates, as W i ( p ) = ( ∂ V i V i )( p ) = ddt ( V i ◦ C p,i )(0) . (10)3.2 The map in adapted coordinates Now we introduce the transformation of the CLV-adapted coordinates on R d , dueto iteration of the map, ϕ : M → M . To do that, we fix an i ≤ d u and focus on therelationship between the curves C p ,i : [ − (cid:15) p , (cid:15) p ] → M and ϕ ◦C p,i : [ − (cid:15) p , (cid:15) p ] → M . Define f p,i := ( C p ,i ) − ◦ ϕ ◦ C p,i , noting that this definition makes sense at apoint t ∈ [ − (cid:15) p , (cid:15) p ] whenever ϕ ( C p,i ( t )) lies in the image of C p ,i . In addition to
Nisha Chandramoorthy, Qiqi Wang continuity of adapted coordinates, the requirement that the definition of f p,i isvalid on [ − (cid:15) p , (cid:15) p ] imposes another constraint on the function p → (cid:15) p . This function, p → (cid:15) p , which determines the length of the curves at each p , must be such thatorbits of f p,in := f p n ,i ◦ · · · ◦ f p ,i ◦ f p,i , n ∈ Z + are possible, for 1 ≤ i ≤ d u . Clearly, 0 is a fixed point of f p,in for all n ∈ Z + , andproduces the orbit p, p , p , · · · , . Intuitively, if orbits of f p,in excluding the fixedpoint exist, say t n := f p,in ( t ), it means that it is possible to control the sizes (cid:15) p n along the orbit of p so that C p n ,i ( t n ) lies in a local unstable manifold of p n , at each n. This constraint can indeed be satisfied based on an assertion in stable-unstablemanifold theory (in particular, see Lemma 2.2.2 of [18]).To summarize, we make a specific choice of p → (cid:15) p such that the curves C p n ,i ateach n lie inside a local unstable manifold at p n , and are tangent to V in := V i ( p n ) . This allows us to obtain expressions for the derivative of a CLV V in +1 with respectto V in , which will in turn enter into the computation of W i . In particular, usingCLV-adapted coordinates, a suitable p → (cid:15) p , as described above, and the definitionof f p,i , ( df p,i /dt )(0) = z i ( p ) . (11)Now we usefully relate the iterates through ϕ of the differential operator on M : ∂ V i , and its analog on R : d/dt , along the trajectory lying in the unstable manifoldof p n . In particular, for the function V i , when combined with Eq. 10, and Eq. 11, d (cid:0) V i ◦ ϕ ◦ C p,i (cid:1) dt (0) = d (cid:0) V i ◦ C p ,i ◦ f p,i (cid:1) dt (0)= z i ( p ) W i ( p ) . (12)3.3 Computation of unstable CLV self-derivativesStarting from Eq. 12, and by definition of CLVs (Eq. 4) W i ( p ) = 1 z i ( p ) ddt (cid:18) Dϕ ◦ C p,i V i ◦ C p,i z i ◦ C p,i (cid:19) (0) (13)= 1( z i ( p )) ddt (cid:16) Dϕ ◦ C p,i (cid:17) (0) V i ( p ) + 1( z i ( p )) Dϕ ( p ) ddt ( V i ◦ C p,i )(0)+ V i ( p ) ddt (cid:18) z i ◦ C p,i (cid:19) (0) (14)By Eq. 10, we can write the second term above as (1 /z i ( p )) Dϕ ( p ) W i ( p ) . The firstterm can be written using the chain rule in terms of the d × d × d second-orderderivative of ϕ , which is denoted as D ϕ. Let the elements of the second-orderderivative of the map be indexed such that D ϕ [ i, j, k ] = ∂ k ∂ j ( ϕ ) i , and let D ϕ : b indicate a d × d matrix resulting from taking the dot product of the last axis of D ϕ and the vector b . Then, Eq. 14 becomes W i ( p ) = 1( z i ( p )) D ϕ ( p ) : V i ( p ) V i ( p ) + 1( z i ( p )) Dϕ ( p ) W i ( p )+ V i ( p ) ddt (cid:18) z i ◦ C p,i (cid:19) (0) (15) n ergodic averaging method to differentiate covariant Lyapunov vectors 9 ddt (cid:18) z i ◦ C p,i (cid:19) (0) = − z i ( p )) ddt (cid:32)(cid:16) ( Dϕ V i ) T Dϕ V i (cid:17) ◦ C p,i (cid:33) (0)= − ( Dϕ V i ) T ( p )( z i ( p )) (cid:16) ( D ϕ : V i )( p ) V i ( p )+ Dϕ ( p ) W i ( p ) (cid:17) = − ( V i ( p )) T ( z i ( p )) (cid:16) (( D ϕ : V i ) V i )( p ) + ( Dϕ W i )( p ) (cid:17) . (16)Substituting Eq. 16 into Eq. 15, we see that Eq. 15 simply projects out the V i ( p )direction. That is, W i ( p ) = (cid:16) I − V i ( p )( V i ( p )) T (cid:17)(cid:32) D ϕ ( p ) : V i ( p ) V i ( p ) + Dϕ ( p ) W i ( p )( z i ( p )) (cid:33) , (17)where I is the d × d Identity matrix. Now, Eq. 17 can be marched forward in timerecursively by replacing W i ( p ) with W i ( p ), and W i ( p ) with W i ( p ). Use thesubscript notation, e.g. W in := W i ( p n ) , and start from a random initial vector ∈ R d as a guess for W i := W i ( p ) . The following iteration is proposed as thedifferential CLV method to obtain W in , n ∈ Z + , ≤ i ≤ d u W in +1 = (cid:16) I − V in +1 ( V in +1 ) T (cid:17)(cid:32) ( D ϕ ) n : V in V in + ( Dϕ ) n W in ( z in ) (cid:33) . (18)In Appendix section D, we show that the above equation converges asympoticallyat an exponential rate, under certain conditions on the LEs. Finally, note thatthe entire procedure above was derived for the unstable CLV self-derivatives. Forthe stable ones, we must apply the same procedure with time reversal since thestable and unstable CLVs are the same, except their roles are exchanged upontime reversal. More precisely, when d u + 1 ≤ i ≤ d , i.e. we must apply the aboveiterative procedure (Eq. 18) by replacing ϕ with the inverse map, ϕ − . In this section, we implement the differential CLV algorithm discussed in the previous section to several examples of low-dimensional chaotic attractors, someof which were introduced in section 2. In every example, the unstable subspace isone-dimensional (a line) and numerical estimates of W are shown. The Pythoncode for the implementation, along with the files needed to generate the plots inthis section, can be found in [7]. Fig. 1
Comparison of the x , x , x components of V computed analytically (shown in or-ange) and numerically (in blue), for the super-contracting Solenoid map ϕ ([ r, t, z ] T ) = s + ( r − s ) /s + (cos t ) / tz/s + (sin t ) / (19)Clearly the parameter s is a contraction factor along the ˆ r and ˆ z directions.In the limit s → ∞ , the attractor of the map, henceforth referred to as thesuper-contracting Solenoid attractor, becomes a space curve. It is described bythe following curve parameterized by the coordinate t , expressed in Cartesiancoordinates: γ ( t ) := x ,n +1 x ,n +1 x ,n +1 = (cid:18) s + cos t (cid:19) cos 2 t (cid:18) s + cos t (cid:19) sin 2 t sin t , (20)where t = arctan( x ,n /x ,n ). As an aside, note that in the ˆ t direction, the map issimply a linear expanding map and hence the ˆ t component of the state vector hasa uniform probability distribution in [0 , π ) . We fix s at 1 throughout. The one-dimensional unstable manifold is given by the curve γ ( t ) defined in Eq. 20. Then, the tangent vector field to the curve, γ (cid:48) ( t ), must be along V ( γ ( t )). This is verifiednumerically in Figure 1, where the numerically computed vector field V agreesclosely with the unit tangent vector field γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) : in each of the subfigures,the components of the two vector fields lie superimposed on each other. Conse-quently, the acceleration along the curve γ ( t ), ∂ γ (cid:48) ( t ) γ (cid:48) ( t ) must be in the direction n ergodic averaging method to differentiate covariant Lyapunov vectors 11 Fig. 2
Comparison of the x , x , x components of W computed analytically (shown inorange) and numerically (in blue), for the super-contracting Solenoid map. Fig. 3
The vector field V is shown for the Solenoid map. The color represents (cid:107) W (cid:107) , whichis the curvature of the attractor. of W ( γ ( t )). In particular, the acceleration in the direction of the unit tangentvector, ∂ γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) ( γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) ), must match W ( γ ( t )). This is also clearly seennumerically. In Figure 2, each component of the two vector fields ∂ γ (cid:48) / (cid:107) γ (cid:48) (cid:107) ( γ (cid:48) / (cid:107) γ (cid:48) (cid:107) ),computed analytically, and W , computed numerically using Eq. 18, are seen to coincide. Thus, the norms of the two vector fields are of course, in close agree-ment as well, as can be seen in Figure 3. Both the analytically computed norm (cid:107) ∂ γ (cid:48) / (cid:107) γ (cid:48) (cid:107) ( γ (cid:48) / (cid:107) γ (cid:48) (cid:107) ) (cid:107) , and the numerically computed (cid:107) W (cid:107) are shown as a colormapon the vector field V = γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) . The plots in Figure 3 are a visualizationof the curvature of the one-dimensional unstable manifold γ ( t ), by definition of Fig. 4
Trajectories of the Lorenz system shown on the x - x plane, after T = 18 (left) and T = 20 (right). The initial conditions were 10001 equi-spaced points on the short line segmentjoining (-0.01,0,1) and (0.01,0,1). curvature. The final results of the analytical curvature calculations are providedin Appendix section B.4.2 Numerical verification of the curvature of the Lorenz attractorNext we consider the well-known Lorenz’63 system, given by the following systemof ODEs: ddt x x x = F ([ x , x , x ] T ) := x − x ) x (28 − x ) − x x x − x . (21)The map ϕ here, is defined to be a time-discretized form of the above system ofODEs. In particular, we use a second-order Runge-Kutta scheme with a timestepof δt = 0 .
01. The map ϕ ( p ) = p , is the time-integrated solution after time δt ,starting from p := [ x , x , x ] T ∈ R . The Lorenz’63 map defined this way has thefollowing Lyapunov exponents: λ ≈ . λ ≈ λ ≈ − . . The unstablemanifold, which corresponds to λ is one-dimensional. There is a one-dimensionalcenter manifold tangent to the right hand side of the ODE, F . This correspondsto λ ≈
0, i.e., since clearly F ( p ) ≈ Dϕ ( p ) F ( p ) , the tangent vector F ( p ) ∈ T p R does not show exponential growth or decay under the tangent dynamics. Thus, thismap is not uniformly hyperbolic as per the description in section 2.2. Rather, it isa partially hyperbolic system –a generalization of a uniformly hyperbolic system that allows a center direction – in which the center-unstable manifold is two-dimensional and tangent to F ⊕ E u . The Lorenz attractor nevertheless mimics thestatistical behavior of a uniformly hyperbolic attractor. For instance, the centrallimit theorem holds for H¨older continuous observables and an SRB-type invariantdistribution exists [3]. n ergodic averaging method to differentiate covariant Lyapunov vectors 13 Fig. 5
Comparison between V from an iteration of the tangent dynamics (shown in orange)and V from finite difference of the primal trajectories (in blue). The first column shows thecomponents of V at time T = 18 and the second column at T = 20 . The first, second andthird rows show the x , x , x components of V respectively. In Figure 4, we numerically calculate the one-dimensional unstable manifoldat p := (0 , ,
1) of the Lorenz attractor. We populate the small line segment ([-0.01,0,1], [0.01,0,1]) with 10001 equi-spaced initial conditions. In Figure 4, thesepoints are shown after time evolution for time T = 18 or n = 1800 steps (on theleft) and T = 20 or n = 2000 steps (on the right). These trajectory points atboth times are colored according to the distance between them and the reference trajectory starting at p , at the respective times. The points that are within adistance of 0.1, are considered iterates (at the indicated times) from within a localunstable manifold of p. These points selected approximately from a local unstable manifold of p arealso used to compute the first CLV V along their trajectories, using the following Fig. 6
Comparison between W from the differential CLV method (shown in orange) and W from finite difference (in blue). The first column shows the components of W at time T = 18 and the second column at T = 20 . The first, second and third rows show the x , x , x components of W respectively. finite difference approximation of V : V ( q n ) ≈ p n − q n (cid:107) p n − q n (cid:107) . (22)The 3 components V x i , i = 1 , , V ( q n ) uses onlythe trajectory q, q , · · · , q n and the tangent dynamics along this trajectory, and works as follows: randomly initialize v ( q ) and propagate the tangent dynamicswith repeated normalization. v ( q n +1 ) = Dϕ ( q n ) v ( q n ) , (23) v ( q n +1 ) ←− v ( q n +1 ) / (cid:107) v ( q n +1 ) (cid:107) . (24) n ergodic averaging method to differentiate covariant Lyapunov vectors 15 Carrying this out for n ∈ Z + , similar to a power iteration method for the com-putation of the dominant eigenvector of a matrix, yields a unit vector v ( q n ) thataligns with V ( q n ) . Clearly, this procedure is equivalent to the above-mentionedfinite difference procedure, as long as q n is in a small neighborhood of p n , for thelength of the trajectory considered.Having visualized V along trajectories, we now compute W using our differ-ential CLV method in section 3. To test its correctness, we also compute W usinga finite difference method as follows. As usual, let the reference trajectory alongwhich we require to compute W be p, p , · · · , p N , and assume that we know theCLVs V ( p ) , V ( p ) , · · · , V ( p N ) . Let q, q , · · · , q N and r, r , · · · , r N be two othertrajectories that are at most a distance of O (1) away from the reference trajectory,at each of the N timesteps. Then, according to our preceding discussion, V ( q n ) ≈ − V ( p n ) ≈ p n − q n (cid:107) p n − q n (cid:107) . (25)At each n , we rescale q n and r n along V ( p n ) to obtain the two points i) ˜ q n = p n + (cid:15) qn V ( q n ), ii) ˜ r n = p n + (cid:15) rn V ( r n ). Then, we can approximately compute W ( p n ) as W ( p n ) ≈ (˜ r n − p n ) /(cid:15) rn − (˜ q n − p n ) /(cid:15) qn (cid:107) ˜ r n − ˜ q n (cid:107) . (26)In Figure 6, we plot the three components of W : W x , W x , W x computed usingthe above procedure in blue and the same quantity computed using the differentialCLV algorithm in section 3 in orange. The closeness of the two results indicatesthe correctness of our algorithm. It is also a numerical verification of the fact that V is differentiable along itself in this system, even though it is only partiallyhyperbolic.4.3 Qualitative verification on a perturbed cat mapWe consider a smoothly perturbed Cat map (PCM) (see section 2.3) due toSlipantschuk et al. [26]. The PCM [26] was designed to be an analytic, area-preserving, uniformly hyperbolic map of the torus, whose spectral properties canbe computed analytically. The PCM is given by ϕ ([ x , x ]) = (cid:20) (cid:21) (cid:20) x x (cid:21) + (cid:20) Ψ s ,s ( x ) Ψ s ,s ( x ) (cid:21) , (27)where Ψ s ,s ( x ) := (1 /π ) arctan (cid:16) s sin(2 πx − s ) / (1 − s cos(2 πx − s )) (cid:17) , is a perturbation whose maximum magnitude is controlled by the parameter s and the location of the maximum, by s . Clearly, the original Cat map is recovered at s = 0 . As in the Cat map, the sum of the LEs is 0 but their values are sensitive tothe parameters, with lesser sensitivity to s when compared to s . Unlike the Catmap, the CLVs are no longer uniform in phase space and are also not orthogonal toeach other. In Figure 7, we show the vector fields V and V computed at s = 0 . Fig. 7
The vector fields V (left) and V (right) are shown for the PCM at s = 0 . , s = 0 . Fig. 8
The vector field V is shown for the PCM at s = 0 . , s = 0 .
2. The color representsthe values of (cid:107) W × V (cid:107) , which equals the norm of (cid:107) W (cid:107) multiplied by a sign representingthe orientation with respect to V . and s = 0 . . Notably, non-zero values of s create a curvature in the CLVs, which is again non-uniform in space. We compute the self-derivative of the unstableCLV using our differential CLV method in section 3. By construction, the methodproduces a vector field W that is orthogonal to V . The norm of the computedvectors, (cid:107) W (cid:107) , is shown signed according to its orientation with respect to V . Inparticular, in Figure 8, we plot (cid:107) W × V (cid:107) as a colormap on the vector field V . n ergodic averaging method to differentiate covariant Lyapunov vectors 17 Figure 8 is a qualitative representation of the fact that (cid:13)(cid:13) W (cid:13)(cid:13) is the curvature ofthe unstable manifold, which is everywhere tangent to the plotted vector field V . The V self-derivative W is the acceleration of a particle moving with the velocityfield V . This intuitive picture is mirrored by Figure 8, in which (cid:13)(cid:13) W (cid:13)(cid:13) is higher inregions of velocity changes than where the velocity appears rather uniform (e.g. ina thin strip around the diagonal of the square). The regions of similar magnitudeof acceleration but of opposite sign, reflect the symmetry in the velocity field V about x = x , and moreover indicate the opposite directions of the turns madein those regions by traveling particles.4.4 Qualitative verification on the volume-decreasing perturbed CatWhile the PCM was an example of a symplectic, uniformly hyperbolic system,now we consider a dissipative, uniformly hyperbolic map. We introduce anotherperturbed Cat map, with smooth nonlinear perturbations that cause the resultingmap to be volume-decreasing. The norm of the perturbations is controlled by a setof four parameters s = [ s , s , s , s ] T and the unperturbed Cat map (the originalAnosov Cat) is recovered at s = [0 , , , ϕ ([ x , x ] T ) = (cid:20) (cid:21) (cid:20) x x (cid:21) + (cid:18) s (cid:20) v v (cid:21) + s (cid:20) v v (cid:21)(cid:19) sin(2 π ˜ V · x ) /c + (cid:18) s (cid:20) v v (cid:21) + s (cid:20) v v (cid:21)(cid:19) sin(2 π ˜ V · x ) /c (28)where ˜ V := [ v , v ] T = [5 , − T ∈ R is a rational approximation of the stableCLV of the unperturbed Cat map. Similarly, ˜ V := [ v , v ] T = [8 , T ∈ R isa rational approximation of the unstable CLV of the unperturbed Cat map. Theconstant c serves to normalize the perturbations and is set to c = 2 π ( v + v ) . Thefour parameters together determine the norm and direction of the perturbation. InFigure 9, V in plotted in each case of turning on just one of the four parameters, inorder to isolate their effects. Each subfigure reflects the effect of a single parameter,on V , in comparison to the unperturbed Cat map (in which V is roughly parallelto the line ˜ V ). For instance, when s = [1 , , , T , a perturbation is applied alongthe direction ˜ V , which is approximately along the stable direction of the DCM.The norm of this perturbation varies sinusoidally with the orientation along theapproximately stable direction, ˜ V . As can be seen in the top-left of Figure 9, theCLV V is rather uniform in its own direction but shows a striated pattern in theperpendicular direction, roughly along ˜ V . As another example, the bottom-leftsubfigure shows V at s = [0 , , , T . From Eq. 28, we know that s being non-zero introduces a perturbation, along ˜ V , whose norm varies in the approximatelyunstable direction, ˜ V . This is portrayed in the figure, wherein V appears aswaves, which are seen traveling approximately along ˜ V but the amplitudes of the waves clearly vary in the perpendicular, approximately unstable direction.With this understanding of the effect of each parameter, we expect that V would show a smaller sensitivity, in its own direction, when the norm of the per-turbation is uniform along ˜ V . This is the case when s , s are set to 0. Thisintuition is confirmed by the numerical results obtained on using the differential Fig. 9
The vector field V is shown for the DCM at different parameter choices. The param-eters not indicated are set to 0 in each case. CLV method. As shown in Figure 10, when either s = 1 or s = 1, and the other3 parameters are set to 0, we see that the numerically computed W has a smallernorm, when compared to the other cases.On the bottom row in Figure 10 are the vector fields W when either s or s are set to 1 and the rest to 0. In these cases, the norm of the perturbation variesalong the approximately unstable direction, and this is clearly reflected in thehigher (when compared to the other two cases) magnitudes of W . In addition, thevariation in W itself, which gives information about the second-order derivativeof V , is also consistent with our expectations. For instance, W shows a markedvariation along V when s = 1 (bottom-left of Figure 10). This can be explained by the applied perturbations being sinusoidal in the direction of ˜ V , giving riseto a harmonic functions for the higher-order derivatives along V as well. Finally,when s = 1, (bottom-right of Figure 10), it is easy to observe that, qualitatively,the density of the lines V is reflected in the magnitudes of W . This is not acoincidence, as we shall see in section 5. There, we describe that W partly givesthe variation in the density of the SRB measure on the unstable manifold, due to n ergodic averaging method to differentiate covariant Lyapunov vectors 19 Fig. 10
The vector field V is shown for the DCM, colored according to (cid:13)(cid:13) W × V (cid:13)(cid:13) . Theparameters not indicated as 1 are set to zero in each case. perturbations along V . Now we can see that especially the s = 1 case provides avisualization consistent with this theoretical insight. Particularly, the pronouncedvariation in the unstable direction (bottom-right, Figure 10), mirrors the changesin probability density on the unstable manifold, which is qualitatively measuredby the closeness of the V lines in Figure 9.4.5 Numerical results on the Henon mapAs our final example, we consider the classical H`enon attractor. The Henon mapis the canonical form for a two-dimensional area-decreasing quadratic map [14]: ϕ ([ x , x ] T ) = (cid:20) x + 1 − ax bx (cid:21) . (29)Taking the parameters a and b at their standard values of a = 1 . b = 0 .
3, weobtain the Henon attractor, on which the CLVs are shown in Figure 11. At theseparameter values, the H´enon attractor is nonhyperbolic due to the presence oftangencies between the stable and unstable manifolds [2]. On this map, we applythe differential CLV method we derived in section 3, and the resulting W is shownin Figure 12. The CLVs may not be differentiable everywhere, as seen by the largemagnitudes of the numerically computed W at the sharp turns in the attractor. In Figure 13, we dissect the derivatives further to investigate the issue of dif-ferentiability numerically. In each subfigure, the vector field V is plotted coloredaccording to (cid:107) W (cid:107) ; the points at which (cid:107) W (cid:107) is not in the range indicated by thecolormap, are excluded. From the top row of Figure 13, it is clear that (cid:107) W (cid:107) < . Fig. 11
The CLV V on the henon attractor. Inset is the CLV field in a neighborhood of thefixed point ≈ (0 . , . . Fig. 12
The vector field V is shown for the Henon map. The color represents the V self-derivative norm, (cid:107) W (cid:107) . Fig. 13
The vector field V is shown for the Henon map. The color represents (cid:107) W (cid:107) , thecurvature of the unstable manifold.n ergodic averaging method to differentiate covariant Lyapunov vectors 21 curved side of the attractor, still have a curvature less than 1. On the bottom row,the more rounded portions of the attractor, as expected, have a higher curvaturewhen compared to the previous cases. On the bottom-right, we see that only thecorners and turns appear to have (cid:107) W (cid:107) higher than 100. Among these points, thevariation in the curvature, (cid:107) W (cid:107) , is over six orders of magnitude, with the sharpcorners, having the highest curvatures. In this case, our numerical method for W acts as an indicator for the lack of differentiability at some points. At least in twodimensions, this also turns out to be a detector for uniform hyperbolicity, basedon our discussion in section 2.4. A landmark result in the theory of uniformly hyperbolic systems due to Ruelle([24][25] ; [12] contains a modern proof of the result) is the smooth response of theirstatistics to perturbations. Here we briefly describe this result, called the linearresponse formula, and draw a connection between the formula and the unstableCLV self-derivatives W i , i ≤ d u . Consider a family of uniformly hyperbolic maps ϕ t ∈ C ( M ), where t is a smallparameter around 0. Let the reference map ϕ be written simply as ϕ , and V bea smooth vector field such that ϕ t = ϕ + tV up to first order in t . Let the SRBmeasure associated to ϕ t be µ t : that is, µ t is a ϕ t -invariant probability distributionon M such that for any scalar observable f ∈ L ( M ), the ergodic average startingfrom a p ∈ M Lebesgue-a.e., lim N →∞ (1 /N ) (cid:80) N − n =0 f ( ϕ t ( p )) = (cid:104) f, µ (cid:105) . Ruelle’s linear response theory [24][25] proves the existence of the statistical re-sponse, (cid:104) f, ∂ t µ t (cid:105) , including expressing this quantity as an exponentially convergingseries. The quantity (cid:104) f, ∂ t | µ t (cid:105) represents the derivative with respect to t of er-godic averages or equivalently ensemble averages of observables with respect to theSRB measure, and is of immense interest in practical applications. The statisticalsensitivity (cid:104) f, ∂ t | µ t (cid:105) is useful for sensitivity analysis, uncertainty quantification,model selection etc, in every scientific discipline from climate studies [23,19] toaerodynamic fluid flows [21,15]. The linear response formula [24,25] is as follows: (cid:104) f, ∂ t (cid:12)(cid:12) µ t (cid:105) = ∞ (cid:88) n =0 (cid:104) Df n · V, µ (cid:105) . (30)Although the above series is exponentially converging, previous works [8][10] sug-gest that it is computationally infeasible to calculate the series in its original formwhen V ∈ E u , especially in high-dimensional practical systems. If each term inthe series in regularized by an integration by parts, the resulting form of the linearresponse formula is more amenable to computation. For a simple illustration, wefix the perturbation field V = V i , ≤ i ≤ d u , along an unstable CLV, instead of ageneral unstable perturbation, which would be a linear combination of the unsta-ble CLVs. Applying integration by parts to Eq. 30, and using the fact that ergodic averages converge to ensemble averages (for a full derivation of the following fromEq. 30, see section 5 of [9]), (cid:104) f, ∂ t | µ t (cid:105) = − ∞ (cid:88) m =0 lim N →∞ N N − (cid:88) n =0 f m (cid:32) W i ( q n ) + ( ∂ V i ρ )( q n ) ρ ( q n ) (cid:33) (31) where ρ is the density of the distribution µ along the unstable manifolds. While W i can be computed using the differential CLV method, since the density ρ is unknown, the above formula cannot still be used to compute linear response.However, a recent reformulation of the formula has been derived in [9]. In thatreformulation, the distribution g := W i + ∂ V i ρ /ρ is arrived at using an iterativetrajectory-based algorithm known as the space-split sensitivity or S3 algorithm (seesection 4 of [9]). For completion, we give the final expression for g in AppendixC. We note that the trajectory-based computation of the expression involves abyproduct of the differentital CLV method: ∂ V i (1 /z i ); this can be computed us-ing Eq. 16. Thus, the differential CLV method, with and without the orthogonalprojection step (Eq. 18) can supply both the terms – ∂ V i (1 /z i ) and W i – that areneeded to compute linear response using Eq. 31. In this work, we have derived a numerical method, called the differential CLVmethod, to compute the derivatives of Covariant Lyapunov Vectors along their owndirections: the CLV self-derivatives. These directional derivatives exist in smoothuniformly hyperbolic systems with compact attractors. We demonstrate the appli-cation of the differential CLV method on a variety of low-dimensional attractorsincluding a quasi-hyperbolic attractor (Lorenz’63) and a non-hyperbolic attractor(H´enon). In the two-dimensional uniformly hyperbolic systems considered, includ-ing perturbations of the Cat map, our method provides rich visualizations of thecurvature of the one-dimensional unstable manifold. The CLV self-derivatives arefundamentally linked to the statistical linear response of a chaotic attractor. Thelink is through their utility to compute the divergence of perturbations on theunstable manifold, with respect to the SRB-type measure. We draw a concreteconnection between the computed derivatives and an efficient method to differen-tiate statistics with respect to system parameters in uniformly hyperbolic systems.We hope that the differential CLV method to obtain the CLV self-derivatives andthe numerical experiments in this work will invite applications to linear response,and beyond, of these fundamental objects.
Acknowledgments:
We offer our sincere thanks to Dr. Jizhou Li for com-ments on this manuscript. n ergodic averaging method to differentiate covariant Lyapunov vectors 23
References
1. Abramov, R.V., Majda, A.J.: New approximations and tests of linear fluctuation-responsefor chaotic nonlinear forced-dissipative dynamical systems. Journal of Nonlinear Science (3), 303–341 (2008)2. Arai, Z.: On hyperbolic plateaus of the hnon map. Experimental Mathematics (2),181–188 (2007). DOI 10.1080/10586458.2007.10128992. URL
3. Arajo, V., Melbourne, I., Varandas, P.: Rapid Mixing for the Lorenz Attractor and Sta-tistical Limit Laws for Their Time-1 Maps. Commun. Math. Phys. (3), 901–938(2015). DOI 10.1007/s00220-015-2471-0. URL http://link.springer.com/10.1007/s00220-015-2471-0
4. Arnold, L.: The Multiplicative Ergodic Theorem on Bundles and Manifolds, pp. 163–199.Springer Berlin Heidelberg, Berlin, Heidelberg (1998). DOI 10.1007/978-3-662-12878-7 4.URL https://doi.org/10.1007/978-3-662-12878-7_4
5. Beeson, R., Sri Namachchivaya, N.: Particle filtering for chaotic dynamical systems usingfuture right-singular vectors. Nonlinear Dyn (2020). DOI 10.1007/s11071-020-05727-y.URL http://link.springer.com/10.1007/s11071-020-05727-y
6. Cencini, M., Ginelli, F.: Lyapunov analysis: from dynamical systems theory to applica-tions. Journal of Physics A: Mathematical and Theoretical (25), 250301 (2013). DOI10.1088/1751-8113/46/25/250301. URL https://doi.org/10.1088%2F1751-8113%2F46%2F25%2F250301
7. Chandramoorthy, N.: nishachandramoorthy/s3: Differential clv method. 1.0.0 (2020). DOI10.5281/zenodo.3941678. URL https://doi.org/10.5281/zenodo.3941678
8. Chandramoorthy, N., Fernandez, P., Talnikar, C., Wang, Q.: Feasibility analysis of ensem-ble sensitivity computation inturbulent flows. AIAA Journal (10), 4514–4526 (2019).DOI 10.2514/1.J058127. URL https://doi.org/10.2514/1.J058127
9. Chandramoorthy, N., Wang, Q.: A computable realization of Ruelle’s formula for linearresponse of statistics in chaotic systems. arXiv e-prints arXiv:2002.04117 (2020)10. Eyink, G., Haine, T., Lea, D.: Ruelle’s linear response formula, ensemble adjoint schemesand l´evy flights. Nonlinearity , 1867 (2004). DOI 10.1088/0951-7715/17/5/01611. Ginelli, F., Chat´e, H., Livi, R., Politi, A.: Covariant lyapunov vectors. Journal of PhysicsA: Mathematical and Theoretical (25), 254005 (2013). DOI 10.1088/1751-8113/46/25/254005. URL https://doi.org/10.1088%2F1751-8113%2F46%2F25%2F254005
12. Gou¨ezel, S., Liverani, C., et al.: Compact locally maximal hyperbolic sets for smoothmaps: fine statistical properties. Journal of Differential Geometry (3), 433–477 (2008).DOI 10.4310/jdg/121379818413. Hasselblatt, B., Wilkinson, A.: Prevalence of non-Lipschitz Anosov foliations. ErgodicTheory and Dynamical Systems (3), 643–656 (1999). DOI 10.1017/S0143385799133868.Publisher: Cambridge University Press14. H´enon, M.: A Two-dimensional Mapping with a Strange Attractor, pp. 94–102. SpringerNew York, New York, NY (2004)15. Huhn, F., Magri, L.: Optimisation of chaotically perturbed acoustic limit cycles. NonlinearDyn (2), 1641–1657 (2020). DOI 10.1007/s11071-020-05582-x. URL http://link.springer.com/10.1007/s11071-020-05582-x
16. Katok, A., Hasselblatt, B.: Introduction to the modern theory of dynamical systems,vol. 54. Cambridge university press (1997). DOI 10.1017/CBO978051180918717. Kuptsov, P.V., Parlitz, U.: Theory and Computation of Covariant Lyapunov Vectors.J Nonlinear Sci (5), 727–762 (2012). DOI 10.1007/s00332-012-9126-5. URL https://doi.org/10.1007/s00332-012-9126-5
18. Ledrappier, F., Young, L.S.: The metric entropy of diffeomorphisms: Part i: Character-ization of measures satisfying pesin’s entropy formula. Annals of Mathematics (3),509–539 (1985). URL
19. Lucarini, V.: Response operators for markov processes in a finite state space: Radius ofconvergence and link to the response theory for axiom a systems. Journal of StatisticalPhysics (2), 312–333 (2016). DOI 10.1007/s10955-015-1409-4. URL https://doi.org/10.1007/s10955-015-1409-4
20. Ni, A.: Sensitivity analysis on chaotic dynamical systems by non-intrusive least squaresadjoint shadowing (nilsas). arXiv preprint arXiv:1801.08674 (2018)21. Ni, A.: Hyperbolicity, shadowing directions and sensitivity analysis of a turbulent three-dimensional flow. Journal of Fluid Mechanics , 644669 (2019). DOI 10.1017/jfm.2018.9864 Nisha Chandramoorthy, Qiqi Wang22. Noethen, F.: A projector-based convergence proof of the ginelli algorithm for covariantlyapunov vectors. Physica D: Nonlinear Phenomena , 18 – 34 (2019). DOI https://doi.org/10.1016/j.physd.2019.02.012. URL
23. Ragone, F., Lucarini, V., Lunkeit, F.: A new framework for climate sensitivity andprediction: a modelling perspective. Climate Dynamics , 1459–1471 (2016). DOI10.1007/s00382-015-2657-324. Ruelle, D.: Differentiation of srb states. Communications in Mathematical Physics ,227–241 (1997). DOI 10.1007/s00220005013425. Ruelle, D.: Differentiation of srb states: correction and complements. Communications inmathematical physics , 185–190 (2003). DOI 10.1007/s00220-002-0779-z26. Slipantschuk, J., Bandtlow, O.F., Just, W.: Complete spectral data for analytic Anosovmaps of the torus. Nonlinearity (7), 2667–2686 (2017). DOI 10.1088/1361-6544/aa700f.URL https://iopscience.iop.org/article/10.1088/1361-6544/aa700f
27. Wang, Q.: Convergence of the least squares shadowing method for computing derivativeof ergodic averages. SIAM Journal on Numerical Analysis , 156–170 (2014). DOI10.1137/13091706528. Young, L.S.: Statistical properties of dynamical systems with some hyperbolicity. Annalsof Mathematics , 585–650 (1998). DOI 10.2307/120960 A The lack of differentiability of CLVs
In general, we say that a subspace E is H¨older continuous on M if there exist constants K, δ > α ∈ (0 ,
1] such that (cid:107) E ( p ) − E ( q ) (cid:107) ∗ ≤ K (cid:107) p − q (cid:107) α , whenever p, q ∈ M are suchthat (cid:107) p − q (cid:107) ≤ δ. As mentioned in section 2.4, the subspaces E u , E s are H¨older continuousspaces with an α that is rarely equal to 1. The reader is referred to classical texts such as [16](Chapter 19) or [13] for a detailed exposition on H¨older structures on hyperbolic sets.There, the norm (cid:107)·(cid:107) ∗ uses the adapted coordinate system introduced in section 3.1. Theset of H¨older continuous functions themselves, is independent of the coordinate system how-ever. The norm (cid:107)·(cid:107) ∗ used in the above references (e.g. in Theorem 19.1.6 of [16]), for ourparticular choice of adapted coordinates introduced in 3.1, results in the following definitions,which are exactly what one might expect. Suppose (cid:107) p − q (cid:107) ≤ δ , and Q ( p ) , Q ( q ) are ma-trix representations of the CLV basis whose i th columns respectively are V i ( p ) , V i ( q ). Then, (cid:107) E u ( p ) − E u ( q ) (cid:107) ∗ := (cid:107) Q ( p )[: , d u ] − Q ( q )[: , d u ] (cid:107) where the norm on the right hand sideis a matrix norm on R d × d u , say the induced 2-norm. Here we have again used programmaticnotation: given a matrix A , A [: , i : j ] refers to the columns of A from i to j , limits included. Sim-ilarly, for E s , (cid:107) E s ( p ) − E s ( q ) (cid:107) ∗ := (cid:107) Q ( p )[: , d u + 1 : d ] − Q ( q )[: , d u + 1 : d ] (cid:107) . Consistent withthese definitions, for a one-dimensional E i , we have (cid:13)(cid:13) E i ( p ) − E i ( q ) (cid:13)(cid:13) := (cid:13)(cid:13) V i ( p ) − V i ( q ) (cid:13)(cid:13) , which is simply the 2-norm on R d . B Computations on the super-contracting Solenoid attractor
The super-contracting Solenoid attractor is the curve γ : [0 , π ] → R (defined in Eq. 20)parameterized by a single parameter t . Since we have a closed form expression for the one-dimensional attractor, we can compute its tangent vector field, as: dγdt = − r ( t ) sin 2 t − (sin t cos 2 t ) / r ( t ) cos 2 t − (sin t sin 2 t ) / t , (32)where r ( t ) = (cid:18) s + cos t (cid:19) . n ergodic averaging method to differentiate covariant Lyapunov vectors 25As explained in section 4.1, V ( t ) = γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) . Further, we analytically calculate that ∂ γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) (cid:0) γ (cid:48) ( t ) / (cid:107) γ (cid:48) ( t ) (cid:107) (cid:1) = 12 − (193 cos t + 392 cos 2 t + 267 cos 3 t + 68 cos 4 t + 6 cos 5 t + 36) /c − (189 sin t + 392 sin 2 t + 267 sin 3 t + 68 sin 4 t + 6 sin 5 t ) /c − (19 sin t + 8 sin t cos t + 2 sin t cos 2 t − t cos t ) / ( c / (cid:2) − sin 2 t/r cos 2 t/r (cid:3) γ (cid:48) ( t ) (cid:107) γ (cid:48) ( t ) (cid:107) (33)where c := 2(16 cos t + 2 cos 2 t + 19) / . In Figures 3 and 2, we observe that the vector field W computed using the differential CLVmethod (Eq. 18), matches almost exactly against the above expression in Eq. 33. C Computation of linear response
For a perturbation vector field V = V i , ≤ i ≤ d u , the linear response of the statistic (cid:104) f, µ t (cid:105) ,as derived in [9] is (cid:104) f, ∂ t | µ t (cid:105) = ∞ (cid:88) n =0 (cid:104) Df n · V i , µ (cid:105) (34)= ∞ (cid:88) n =0 (cid:104) f n g i , µ (cid:105) . (35)The distribution g i is derived, using the space-split sensitivity algorithm in [9], to be: g i := ∞ (cid:88) k =1 z i − k (cid:81) k − j =0 z ij − k ∂ V i − k (1 /z i − k )) − k . (36)As in [9], we could use the following scalar equation to approximate the distribution g i : g i ( p n +1 ) = g i ( p n ) /z i ( p n ) − (cid:0) ∂ V i (1 /z i ) (cid:1) ( p n ) . (37)in which the term ∂ V i (1 /z i ) must be evaluated along a trajectory. Note that through ouralgorithm for W i , via Eq. 18, we have also essentially obtained ∂ V i (1 /z i ) as a byproduct,through Eq. 16. D Convergence of the differential CLV method
In this section, we show that when 1 ≤ i ≤ d u , iteration via Eq. 18 converges to W i ( q n ).Moreover, the asymptotic convergence is exponentially fast under some conditions. Fix a ref-erence trajectory q, q , · · · , , and use the notation f n to denote f ( q n ). Let W i , W i , · · · and˜ W i , ˜ W i , · · · be two sequences of vectors generated by iterating Eq. 18. Then, from Eq. 18, (cid:13)(cid:13)(cid:13) W in − ˜ W in (cid:13)(cid:13)(cid:13) = 1 (cid:81) n − m =0 ( z im ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) n − (cid:89) m =0 (cid:16) ( I − V im +1 ( V im +1 ) T )( Dϕ ) m (cid:17) ( W i − ˜ W i ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) . (38)We can apply Oseledets MET to the cocycle Coc( q m , n ) = (cid:81) n − k =0 ( I − V im + k +1 ( V im + k +1 ) T )( Dϕ ) m + k , and to the Jacobian cocycle to obtain the following asymptotic inequality. In particular, usingthe relationship Eq. 6, we get that for every (cid:15) >
0, there exists an N ∈ N such that for all n ≥ N , (cid:13)(cid:13)(cid:13) W in − ˜ W in (cid:13)(cid:13)(cid:13) = 1 (cid:81) n − m =0 ( z im ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) n − (cid:89) m =0 (cid:16) ( I − V im +1 ( V im +1 ) T )( Dϕ ) m (cid:17) ( W i − ˜ W i ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ e − n ( λ i − (cid:15) ) e n ( ω i + (cid:15) ) (cid:13)(cid:13)(cid:13) W i − ˜ W i (cid:13)(cid:13)(cid:13) . (39)6 Nisha Chandramoorthy, Qiqi WangIn the above inequality 39, ω i := max j (cid:54) = i, ≤ j ≤ d u λ j . Thus, asymptotic exponential convergenceis guaranteed whenever 2 λ i ≥ ω i , which is of course true when i = 1= 1