Computation of attractor dimension and maximal sums of Lyapunov exponents using polynomial optimization (original) (raw)

(1Division of Mathematics, University of Dundee, Dundee, DD1 4HN, United Kingdom
2Department of Mathematics and Statistics, University of Victoria, Victoria, BC, V8P 5C2, Canada)

Abstract

Two approaches are presented for computing upper bounds on Lyapunov exponents and their sums, and on Lyapunov dimension, among all trajectories of a dynamical system governed by ordinary differential equations. The first approach expresses a sum of Lyapunov exponents as a time average in an augmented dynamical system and then applies methods for bounding time averages. This generalizes the work of Oeri & Goluskin (Nonlinearity 36:5378–5400, 2023), who bounded the single leading Lyapunov exponent. The second approach considers a different augmented dynamical system, where bounds on sums of Lyapunov exponents are implied by stability of certain sets, and such stability is verified using Lyapunov function methods. Both of our approaches also can be adapted to directly compute bounds on Lyapunov dimension, which in turn implies bounds on the fractal dimension of a global attractor. For systems of ordinary differential equations with polynomial right-hand sides, all of our bounding formulations lead to polynomial optimization problems with sum-of-squares constraints. These sum-of-squares problems can be solved computationally for any particular system to yield numerical bounds, provided the number of variables and degree of polynomials is not prohibitive. Most of our upper bounds are proven to be sharp under relatively weak assumptions. In the case of the polynomial optimization problems, sharpness means that upper bounds converge to the exact values as polynomial degrees are raised. Computational examples demonstrate upper bounds that are sharp to several digits, including for a six-dimensional dynamical system where sums of Lyapunov exponents are maximized on periodic orbits.

1 Introduction

The detection and characterization of chaos is a major computational challenge for many dynamical systems. Chaotic properties are often quantified via the spectrum of Lyapunov exponents (LEs), which can be computed over individual trajectories to obtain local information about that trajectory and, perhaps, global implications such as existence of hyperchaos or estimates of attractor dimension. In this paper we consider the problem of finding upper bounds on sums of LEs among all possible trajectories, without knowledge of any particular trajectories. Upper bounds on sums of LEs can be used to construct upper bounds on the fractal dimension of the global attractor. Our upper bounds are complementary to typical trajectory-based computations, which give corresponding lower bounds on sums of LEs and on Lyapunov dimension of global attractors.

The LE spectrum quantifies the asymptotic rates at which linearized perturbations in different directions grow or decay along trajectories of a dynamical system. For dynamics in nn dimensions there are nn exponents, defined in decreasing order. Partial sums of the first kk exponents describe the exponential rates of growth or contraction of kk-dimensional volumes in state space, so that the first exponent governs the stretching of lines, the sum of the first two governs the growth of areas, and so on. In dissipative systems, these sums become negative beyond a certain value of kk, reflecting the contraction of higher-dimensional volumes. One way of defining a value at which this contraction begins leads to the Lyapunov dimension, which is non-integer in general, and which provides a natural connection between dynamical stability and the fractal geometry of invariant sets such as strange attractors.

We consider autonomous, continuous-time dynamics governed by an ordinary differential equation (ODE) system,

with F∈C1​(Ω,ℝn)F\in C^{1}(\Omega,\mathbb{R}^{n}), where Ck​(A,B)C^{k}(A,B) denotes the space of functions mapping AA to BB that are kk times continuously differentiable. The state space Ω⊂ℝn\Omega\subset\mathbb{R}^{n} may be all of ℝn\mathbb{R}^{n}, a subdomain of ℝn\mathbb{R}^{n} or a lower-dimensional manifold embedded in ℝn\mathbb{R}^{n}.

In the present work, we apply existing methods for proving global properties of ODE dynamical systems to the problem of bounding LEs and their sums. These existing methods rely on finding auxiliary functions V∈C1​(Ω,ℝ)V\in C^{1}(\Omega,\mathbb{R}) subject to certain pointwise inequalities that, in turn, imply global statements about the dynamics. The best known examples of auxiliary functions are Lyapunov functions, whose inequality conditions imply statements about stability. (Our use of Lyapunov functions to study Lyapunov exponents is not typical, despite both objects bearing the same person’s name.) Among the lesser known types of auxiliary functions are those whose inequality conditions imply bounds on infinite-time averages. For each quantity studied in the present work, we formulate two different upper bounds: a “shifted spectrum” approach based on Lyapunov function conditions for stability, and a “sphere projection” approach based on auxiliary function conditions for bounding time averages.

The shifted spectrum approach that we present in section 2.3 uses conditions on VV that are similar to standard Lyapunov function conditions. In stability theory, the pointwise inequality conditions required of a Lyapunov function may be, for instance,

where D​VDV denotes the gradient of VV with respect to the state space variable XX. Here F​(X)⋅D​V​(X)F(X)\cdot DV(X) is the Lie derivative of VV with respect to the flow since the usual chain rule gives dd​t​V​(X​(t))=V​(X​(t))⋅D​V​(X​(t))\frac{\mathrm{d}}{\mathrm{d}t}V(X(t))=V(X(t))\cdot DV(X(t)). The existence of VV satisfying 2 would imply that the set of X∈ΩX\in\Omega where V​(X)=0V(X)=0 is Lyapunov stable, meaning that trajectories remain arbitrarily close to this set if they start sufficiently close [17]. Until about 25 years ago, there was no broadly applicable method for finding VV subject to pointwise inequalities like those in 2. Now, however, it is often possible to construct such VV computationally using optimization methods based on sum-of-squares (SOS) polynomials [39, 36]. In particular, if the function FF defining the dynamics is polynomial in the components of XX, and if VV is sought from a finite-dimensional space of polynomials, then F⋅D​VF\cdot DV is also a polynomial. Deciding whether the polynomials VV and −F⋅D​V-F\cdot DV are nonnegative on Ω\Omega is prohibitively hard in general [32], but this nonnegativity can be ensured by stronger and more tractable SOS conditions. For instance, simply requiring VV to be SOS, meaning that it can be represented as a sum of squares of other polynomials, would imply V≥0V\geq 0 on all of ℝn\mathbb{R}^{n}. Other SOS conditions can imply nonnegativity on Ω\Omega but not on all of ℝn\mathbb{R}^{n}, as explained in section 2.4.

The sphere projection approach that we present in section 2.2 uses conditions on an auxiliary function VV for bounding time averages. For any Φ∈C0​(Ω,ℝ)\Phi\in C^{0}(\Omega,\mathbb{R}) evaluated along a trajectory X​(t)X(t) with initial condition X0X_{0}, define the infinite-time average Φ¯\overline{\Phi} as

On each X​(t)X(t) that remains bounded forward in time, boundedness of V​(X​(t))V(X(t)) implies dd​t​V¯=0\overline{\frac{\mathrm{d}}{\mathrm{d}t}V}=0, and thus F⋅D​V¯=0\overline{F\cdot DV}=0. Thus, on bounded trajectories,

The right-hand upper bound is useful because it does not involve individual trajectories. Since 4 holds for all bounded trajectories and all V∈C1V\in C^{1}, one can maximize the left-hand side over X0X_{0} and minimize the right-hand side over VV to find

Furthermore, the inequality in 5 is an equality when Ω\Omega is a compact set in which trajectories remain [46, 3], and this equality underlies sharpness of our sphere projection approach. When Φ\Phi, VV and FF are polynomial, so is the expression Φ+F⋅D​V\Phi+F\cdot DV on the right-hand side of 5. In this case the min–max problem can be relaxed into a minimization over polynomial VV subject to an SOS constraint, as explained in section 2.4. This constitutes an SOS optimization problem, where the optimization objective is to minimize the resulting upper bound on Φ¯\overline{\Phi}. These ideas were applied by Oeri and Goluskin [34] to bound the single leading LE. Here we extend this approach to sums of LEs and to Lyapunov dimension.

Our two approaches are the first broadly applicable computational methods that can give convergent upper bounds on global maxima of LE sums or Lyapunov dimension. Unlike methods that rely on pointwise optimization over state space [7], our methods remain sharp when the maximizing orbits are periodic, rather than stationary. Furthermore, our methods give a systematic way to search for the maximizing trajectories themselves.

The rest of the paper is organized as follows. Section 2 defines sums of LEs and Lyapunov dimension and then describes our two general approaches to bounding these quantities. After arriving at more general optimization problems, we strengthen their constraints to obtain SOS optimization problems, then we describe how to exploit symmetries. Section 3 presents computational results from applications of our methods to two examples: the chaotic Duffing oscillator, which can be formulated as an autonomous system in 4 variables, and a 3-degree-of-freedom Hamiltonian system adapted from [35], which is an autonomous system in 6 variables. For all LE sums in both examples, we confirm that our computed bounds are sharp to several digits by finding periodic orbits on which the bounds are approximately saturated. Concluding remarks are given in section 4. For completeness, appendix A summarizes necessary material about exterior products, and appendix B proves a non-standard variant of an existence theorem for Lyapunov functions.

2 Optimization problems and SOS relaxations

This section defines the objects we seek to bound and formulates our methods for doing so. section 2.1 give expressions for LE sums and Lyapunov dimension that are useful for what follows. In section 2.2 we describe approaches for bounding these quantities based on projecting tangent space vectors onto the unit sphere. In section 2.3 we describe approaches that instead shift the spectrum of the Jacobian so that tangent vectors remain bounded. Section 2.4 describes how each formulation is relaxed to into SOS optimization problems that can be implemented numerically. Finally, section 2.5 describes how symmetries may be exploited in the optimization problems.

2.1 Tangent vectors, Lyapunov exponents and Lyapunov dimension

In this section we consider autonomous ODE systems on a domain ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n}, just as in 1 but denoted here using lowercase xx and ff,

with f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}). Uppercase XX and FF and domain Ω\Omega will be reserved for associated ODEs on augmented state spaces that combine the dynamics of xx with dynamics on its tangent space. Trajectories are assumed to remain in ℬ\mathcal{B} for all t≥0t\geq 0. In cases where ℬ\mathcal{B} is unbounded, we further assume that each trajectory remains bounded for t≥0t\geq 0, meaning there is no forward-time blowup. The differentiability of ff and the assumption of no blowup guarantee well-posedness for all positive times, meaning there exists a differentiable flow map ϕ:ℝ+×ℬ→ℬ\phi:\mathbb{R}_{+}\times\mathcal{B}\to\mathcal{B} satisfying x​(t)=ϕt​(x0)x(t)=\phi^{t}(x_{0}) for every initial condition x0∈ℬx_{0}\in\mathcal{B} and all t≥0t\geq 0.

The LEs associated with a trajectory x​(t)x(t) quantify the convergence or divergence of infinitesimally nearby trajectories. Linearization of 6 around each point x​(t)x(t) defines the dynamics of the tangent vector y​(t)y(t), which evolves in the tangent space ℝn\mathbb{R}^{n} according to

where D​f​(x)Df(x) denotes the n×nn\times n Jacobian matrix of ff at xx, and y0y_{0} is the initial tangent vector. Equivalently, y​(t)y(t) is given by the linearized finite-time flow map as

We restrict y0y_{0} to the unit sphere,

| 𝕊n−1={x∈ℝn:|x|=1},\mathbb{S}^{n-1}=\{x\in\mathbb{R}^{n}\,:\,|x|=1\}, | (9) | | ------------------------------------------------------------------------------- | --- |

because the tangent vector dynamics are linear in yy and thus are unaffected by normalization of the initial condition. Considering the separation between trajectories starting at x0x_{0} and at x0+ε​y0x_{0}+\varepsilon y_{0}, Taylor expanding the flow map and using 8 gives

| |ϕt​(x0)−ϕt​(x0+ε​y0)|=ε​|y​(t)|+O​(ε2​|y0|2).\displaystyle|\phi^{t}(x_{0})-\phi^{t}(x_{0}+\varepsilon y_{0})|=\varepsilon|y(t)|+O(\varepsilon^{2}|y_{0}|^{2}). | (10) | | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

The tangent vector typically shrinks or grows exponentially like |y​(t)|≈eμ​t|y(t)|\approx e^{\mu t} for some average rate μ​(x0,y0)\mu(x_{0},y_{0}), in which case 10 suggests

| |ϕt​(x0)−ϕt​(x0+ε​y0)|≈ε​|y0|​eμ​(x0,y0)​t.|\phi^{t}(x_{0})-\phi^{t}(x_{0}+\varepsilon y_{0})|\approx\varepsilon|y_{0}|e^{\mu(x_{0},y_{0})t}. | (11) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

If μ<0\mu<0 this approximation can be valid for all time. If μ>0\mu>0 this approximation can be valid only until the time when separation between trajectories grows too large for 10 to apply, but this time approaches infinity as ε→0\varepsilon\to 0. Motivated by the expectation that |y​(t)|≈eμ​t|y(t)|\approx e^{\mu t}, the corresponding LE can be defined precisely as

| μ​(x0,y0)=lim supt→∞1t​log⁡|y​(t)|.\mu(x_{0},y_{0})=\limsup_{t\to\infty}\frac{1}{t}\log|y(t)|. | (12) | | --------------------------------------------------------------------------------------------------------- | ---- |

For each trajectory x​(t)x(t), which is determined by its initial condition x0x_{0}, the LE μ\mu defined by 12 generally depends on the initial tangent vector direction y0y_{0}. The Oseledets ergodic theorem guarantees that different y0y_{0} give at most nn different LEs [40], the largest of which is called the leading LE. We define this by 2.1 and then give equivalent expressions in sections 2.2 and 2.3 that underlie the methods we introduce there. The literature on LEs uses various definitions that are equivalent in many settings [40].

Definition 2.1.

For x​(t)x(t) evolving by the ODE 6, the leading Lyapunov exponent is

| μ1​(x0)=supy0∈𝕊n−1lim supt→∞1t​log⁡|y​(t)|,\mu_{1}(x_{0})=\sup_{y_{0}\in\mathbb{S}^{n-1}}\limsup_{t\to\infty}\frac{1}{t}\log{|y(t)|}, | (13) | | ------------------------------------------------------------------------------------------------------------------------------------------------------ | ---- |

where the tangent vector y​(t)y(t) evolves by 7.

The leading LE is the average growth rate of lengths in the tangent space, maximized over the initial direction. A natural generalization is to consider growth rates of kk-dimensional volumes in the tangent space for any k≤nk\leq n. From this one can define leading sums of LEs and, in turn, the full spectrum of LEs [45]. It is also possible to define the spectrum of LEs directly using singular values of D​ϕtD\phi^{t} [40], but their sums are easier to define, and the sums are what we study in the present work. Let |y1∧⋯∧yk||y_{1}\wedge\cdots\wedge y_{k}| denote the volume of the parallelepiped whose edges coincide with vectors y1,…,yk∈ℝny_{1},\ldots,y_{k}\in\mathbb{R}^{n}. If each yiy_{i} is a tangent vector that evolves in time according to 7, the volume of the corresponding parallelepiped generally grows or shrinks exponentially in time, so we can define an average exponential rate analogously to 12. Maximizing this rate over all initial directions of the tangent vectors gives [45]

| Mk​(x0)=supy1​(0),…,yk​(0)∈ℝn|y1​(0)∧⋯∧yk​(0)|=1lim supt→∞1t​log⁡|y1​(t)∧⋯∧yk​(t)|,M_{k}(x_{0})=\sup_{\begin{subarray}{c}y_{1}(0),\dots,y_{k}(0)\in\mathbb{R}^{n}\\ | y_{1}(0)\wedge\dots\wedge y_{k}(0)|=1\end{subarray}}\limsup_{t\to\infty}\frac{1}{t}\log{|y_{1}(t)\wedge\dots\wedge y_{k}(t)|}, | (14) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

where each yi​(t)y_{i}(t) evolves as in 7. The value Mk​(x0)M_{k}(x_{0}) is the maximum time-averaged growth rate of kk-dimensional volumes in the tangent space for the trajectory x​(t)x(t) starting at x0x_{0}. This maximum growth rate will be attained by almost every choice of tangent vectors, under the assumptions of certain multiplicative ergodic theorems [40].

The set of all y1∧⋯∧yky_{1}\wedge\cdots\wedge y_{k} form a vector space called the exterior product space of kk-vectors on ℝn\mathbb{R}^{n}. The space of kk-vectors on ℝn\mathbb{R}^{n} has dimension nk=n!k!​(n−k)!n_{k}=\tfrac{n!}{k!(n-k)!}. We can therefore identify each y1∧⋯∧yky_{1}\wedge\cdots\wedge y_{k} with a vector denoted by yk∈ℝnky^{k}\in\mathbb{R}^{n_{k}}, whose Euclidean norm |yk||y^{k}| gives the volume of the parallelepiped with edges y1,…,yky_{1},\ldots,y_{k}. For time-dependent kk-vectors where each yi​(t)y_{i}(t) is a tangent vector satisfying the ODE 7, the kk-vectors satisfy a corresponding ODE (see A.4),

Here D​f[k]Df^{[k]} denotes the kt​hk^{th} additive compound for the matrix D​fDf, and later D​f(k)Df^{(k)} denotes the kt​hk^{th} multiplicative compound. These compounds are linear transformations acting on the space of kk-vectors, induced by the linear transformation D​fDf on vectors, and they can be computed explicitly as nk×nkn_{k}\times n_{k} matrices. Exterior products and compound matrices are reviewed in appendix A. The expression 14 for MkM_{k} can be equivalently stated using yk​(t)y^{k}(t), resulting in the following 2.2 that we choose for MkM_{k}.

Definition 2.2.

For x​(t)x(t) evolving by the ODE 6 from initial condition x0x_{0}, the Lyapunov exponents are defined such that the sum of the kk leading Lyapunov exponents is

| Mk​(x0)=supy0k∈𝕊nk−1lim supt→∞1t​log⁡|yk​(t)|,M_{k}(x_{0})=\sup_{y^{k}_{0}\in\mathbb{S}^{{n_{k}}-1}}\limsup_{t\to\infty}\frac{1}{t}\log{|y^{k}(t)|}, | (16) | | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

where yk​(t)y^{k}(t) evolves by 15 from initial condition y0ky^{k}_{0}.

Having defined MkM_{k}, we can define the LE spectrum μ1≥μ2≥⋯≥μn\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n} such the MkM_{k} are necessarily LE sums. When k=1k=1 we already have μ1=M1\mu_{1}=M_{1} since 2.2 reduces to 2.1 when yk=yy^{k}=y and D​f[k]=D​fDf^{[k]}=Df. For k≥2k\geq 2, the LE μk\mu_{k} can be defined recursively by

so indeed Mk=∑i=1kμiM_{k}=\sum_{i=1}^{k}\mu_{i}. There are equivalent ways to define μk\mu_{k} as the fundamental object and then define the MkM_{k} as their partial sums, but for our purposes it is simpler to define MkM_{k} directly by 16.

For systems with attractors, a primary use of LE sums is to define and study the dimension of an attractor. 2.3 gives the Kaplan–Yorke formula [16] for the global Lyapunov dimension dLd_{L} of an ODE system; see [18] for equivalent definitions. Hausdorff dimension is perhaps a more natural way to define the dimension of a fractal set in general. For the global attractor of a dynamical system, Hausdorff dimension can be very difficult to estimate directly, but it is bounded above by dLd_{L} [8]. Thus, the upper bounds we produce for dLd_{L} are upper bounds on Hausdorff dimension also.

Definition 2.3.

For x​(t)x(t) evolving in ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n} by the ODE 6, let jj be the smallest integer such that supx0∈ℬMj+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{j+1}(x_{0})<0. The global Lyapunov dimension over ℬ\mathcal{B} is

For particular ODEs, the Kaplan–Yorke formula is usually evaluated numerically along a chaotic trajectory, rather than being maximized over all trajectories. The result indicates a fractal dimension of the chaotic attractor, rather than of the global attractor. Global Lyapunov dimension pertains to the latter, and the so-called Eden conjecture [9, 22] predicts that the maximum in 18 will be attained on predicts that the maximum in 18 is attained on an equilibrium or periodic orbit, not a chaotic orbit. 2.4 states the fact that j≤dL<j+1j\leq d_{L}<j+1, which we prove for completeness, along with a restatement of the definition of dLd_{L} that we will use in sections 2.2 and 2.3.

Lemma 2.4.

The Lyapunov dimension dLd_{L} has the following properties.

1. 1. For the integer jj in 2.3, Lyapunov dimension is bounded by j≤dL<j+1j\leq d_{L}<j+1.
2. 1. The expression 18 for dLd_{L} is equivalent to

      dL=infB∈[0,1)(j+B)s.t.B​[Mj​(x0)−Mj+1​(x0)]−Mj​(x0)≥0​∀x0∈ℬ.d_{L}=\inf_{B\in[0,1)}(j+B)\quad\text{s.t.}\quad B\left[M_{j}(x_{0})-M_{j+1}(x_{0})\right]-M_{j}(x_{0})\geq 0~~\forall~x_{0}\in\mathcal{B}.	(19)

Proof.

The definition of jj implies that there exists ε>0\varepsilon>0 such that

and that supx0∈ℬMj​(x0)≥0\sup_{x_{0}\in\mathcal{B}}M_{j}(x_{0})\geq 0. From these facts it follows that the supremum in the definition 18 of dLd_{L} cannot be negative, which gives the lower bound j≤dLj\leq d_{L} in the first part of the lemma. For the upper bound, note that

which gives dL<j+1d_{L}<j+1.

For the second part of the lemma, note that the supremum in the definition 18 of dLd_{L} is equivalent to

Since −μj+1=Mj−Mj+1>0-\mu_{j+1}=M_{j}-M_{j+1}>0, rearranging the constraint in 22 gives the constraint in 19 as an equivalent form. This infimum is in [0,1)[0,1) and so is unchanged by the B∈[0,1)B\in[0,1) constraint in 19. Thus expression 19 for dLd_{L} is equivalent to expression 18 in the definition. ∎

Remark 2.5.

Combining both parts of 2.4 implies that, if we can find B∈ℝB\in\mathbb{R} and k≥jk\geq j where

for all x0∈ℬx_{0}\in\mathcal{B}, then the Lyapunov dimension is bounded above by dL≤k+Bd_{L}\leq k+B. Note that 23 cannot hold for all x0x_{0} unless B≥0B\geq 0.

2.2 Sphere projection approach

Our first approach to bounding the LE sums MkM_{k} of 2.2 uses the formulation 5 for bounding infinite-time averages. This is a generalization of the approach in [34], which was described only for the k=1k=1 case. First note that the time average in the definition of MkM_{k} can be expressed as a time integral along trajectories in the augmented state space for x​(t)x(t) and yk​(t)y^{k}(t):

| lim supt→∞1t​log⁡|yk​(t)|\displaystyle\limsup_{t\to\infty}\frac{1}{t}\log{|y^{k}(t)|} | =lim supT→∞1T​∫0Tdd​t​log⁡|yk​(t)|​d​t\displaystyle=\limsup_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{\mathrm{d}}{\mathrm{d}t}{\log{|y^{k}(t)|}}\,\mathrm{d}t | (24) | | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- | | =lim supT→∞1T​∫0T(yk)𝖳​dd​t​yk|yk|2​dt\displaystyle=\limsup_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{(y^{k})^{\mathsf{T}}\frac{\mathrm{d}}{\mathrm{d}t}y^{k}}{|y^{k}|^{2}}\,\mathrm{d}t | (25) | | | =lim supT→∞1T​∫0T(yk)𝖳​D​f[k]​(x)​yk|yk|2​dt,\displaystyle=\limsup_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{(y^{k})^{\mathsf{T}}Df^{[k]}(x)y^{k}}{|y^{k}|^{2}}\,\mathrm{d}t, | (26) | |

where (yk)𝖳(y^{k})^{\mathsf{T}} is the transpose of the vector yky^{k}. The formulation 5 for bounding time averages can be applied to the quantity in 26 only if each (x,yk)(x,y^{k}) remains bounded forward in time, but yk​(t)y^{k}(t) will be unbounded when MkM_{k} is positive. The same tangent space dynamics can be captured with bounded variables by projecting yky^{k} onto a sphere, by introducing zk=yk/|yk|z^{k}=y^{k}/|y^{k}|. Expressing the right-hand side of 26 in terms of zkz^{k} and then maximizing both sides over the initial tangent space directions z0kz^{k}_{0} gives

where the overline denotes a time average of

along trajectories in the (x,zk)(x,z^{k}) state space. The ODE governing zk​(t)z^{k}(t) follows from the ODE 15 for yk​(t)y^{k}(t). The x​(t)x(t) dynamics together with the projected dynamics of kk-vectors in its tangent space gives the ODE system

on the augmented state space ℬ×𝕊nk−1\mathcal{B}\times\mathbb{S}^{n_{k}-1}, where

Representations of MkM_{k} as time averages over the projected tangent space dynamics 29 have often appeared in the study of stochastic dynamics, where they are called Furstenberg–Khasminskii formulas [2, chapter 6]. (In the stochastic case, Φk\Phi_{k} is defined with additional terms and is averaged with respect to stationary measures.)

Auxiliary function methods for bounding time averages can be applied to MkM_{k} by using its representation 27 as a time average in the augmented state space for (x,zk)(x,z^{k}). In particular, we directly apply the formula 5 with state variable X=(x,zk)X=(x,z^{k}) on domain ℬ×𝕊nk−1\mathcal{B}\times\mathbb{S}^{n_{k}-1} and with F​(X)F(X) being the right-hand side of 29. The resulting inequality is the first part of proposition 2.6 below, for which no further proof is needed. If F∈C1F\in C^{1} and the domain is compact, then 5 holds with equality, as proved in [46]. This yields the second part of proposition 2.6 since assuming f∈C2f\in C^{2} gives that F∈C1F\in C^{1}. The k=1k=1 special case of proposition 2.6 is proposition 1 in [34].

Proposition 2.6.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}) and dynamics forward invariant on ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n}. Let Φk​(x,zk)\Phi_{k}(x,z^{k}) and ℓk​(x,zk)\ell_{k}(x,z^{k}) be defined as in 28 and 30, respectively.

1. 1. Suppose each trajectory x​(t)x(t) in ℬ\mathcal{B} is bounded for t≥0t\geq 0. For each k∈{1,…,n}k\in\{1,\ldots,n\}, the maximum leading LE sum MkM_{k} among trajectories in ℬ\mathcal{B} is bounded above by

      supx0∈ℬMk​(x0)≤infV∈C1​(ℬ×ℝnk)supx∈ℬzk∈𝕊nk−1[Φk+f⋅Dx​V+ℓk⋅Dzk​V],\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0})\leq\inf_{V\in C^{1}\left(\mathcal{B}\times\mathbb{R}^{n_{k}}\right)}\sup_{\begin{subarray}{c}x\in\mathcal{B}\hskip 13.0pt\\ z^{k}\in\mathbb{S}^{{n_{k}}-1}\end{subarray}}\left[\Phi_{k}+f\cdot D_{x}V+\ell_{k}\cdot D_{z^{k}}V\right],	(31)
      where Dx​VD_{x}V and Dzk​VD_{z^{k}}V denote the gradients of V​(x,zk)V(x,z^{k}) with respect to xx and zkz^{k}, respectively.

2. 1. Suppose ℬ\mathcal{B} is compact and f∈C2​(ℬ,ℝn)f\in C^{2}(\mathcal{B},\mathbb{R}^{n}). Then 31 holds with equality.

The optimization problem on the right-hand side of 31 may be intractable, but it can be relaxed into computationally tractable optimization problems that also yield upper bounds on MkM_{k}. Relaxations using SOS polynomials are given in section 2.4 below.

The Lyapunov dimension 18 can be bounded above using similar ideas. One way is to apply upper bounds Mk≤BkM_{k}\leq B_{k} that have already been found using the right-hand side of 31 or its SOS relaxations. If such bounds are computed with Bk≥0B_{k}\geq 0 but Bk+1<0B_{k+1}<0 for some k∈{1,…,n}k\in\{1,\ldots,n\}, then proposition 2.7 below gives an upper bounds on dLd_{L}. To apply this result, one does not need to know whether kk is the minimal integer jj appearing in the definition 18 of dLd_{L}. It is enough to know that Bk+1<0B_{k+1}<0, which implies k≥jk\geq j.

Proposition 2.7.

Suppose that, for some k∈{1,…,n}k\in\{1,\ldots,n\} and values Bk≥0B_{k}\geq 0 and Bk+1<0B_{k+1}<0,

Then, the Lyapunov dimension 18 is bounded above by

Proof.

The meaning of jj in 2.3 and the fact that Bk+1<0B_{k+1}<0 together imply k≥jk\geq j. According to 2.5, it suffices for the inequality 23 to hold for all x0x_{0} with B=BkBk−Bk+1B=\frac{B_{k}}{B_{k}-B_{k+1}}. For any x0∈ℬx_{0}\in\mathcal{B}, the assumptions give Mk​(x0)≤BkM_{k}(x_{0})\leq B_{k} and 0<−Bk+1≤−Mk+1​(x0)0<-B_{k+1}\leq-M_{k+1}(x_{0}), so

and

Dividing by the positive quantity Bk−Bk+1B_{k}-B_{k+1} and rearranging gives

which is the required inequality 23. This concludes the proof. ∎

The upper bound on Lyapunov dimension given by proposition 2.7 cannot be sharp for certain systems, even if BkB_{k} and Bk+1B_{k+1} are each sharp bounds in 32. The reason is that the inequality 34, when maximized over x0x_{0}, is an equality only if Mk​(x0)M_{k}(x_{0}) and Mk+1​(x0)M_{k+1}(x_{0}) attain their maxima on the same orbits. This is not true in general; see section 3.2 for a counterexample. We thus give a second formulation for upper bounds on dLd_{L} that is generally sharper than proposition 2.7. This formulation directly considers the ratio in the Kaplan–Yorke formula 18, rather than separately bounding the MkM_{k} and Mk+1M_{k+1} terms in that ratio. In particular, we consider the x​(t)x(t) dynamics together with the projected dynamics of kk-vectors and (k+1)(k+1)-vectors in its tangent space,

where ℓk\ell_{k} is defined by 30, and likewise for ℓk+1\ell_{k+1}. The augmented state space of this system is (x,zk,zk+1)∈ℬ×Snk−1×𝕊nk+1−1(x,z^{k},z^{k+1})\in\mathcal{B}\times{S}^{n_{k}-1}\times\mathbb{S}^{n_{k+1}-1}. The constraint in 19 from 2.4 can then be expressed in terms of time averages in the ODE system 37. For this we define

and we consider the condition

where the time average denoted by the overbar is along trajectories of 37. Condition 39 is equivalent to the first condition in 19 of 2.4, as explained below in the proof of proposition 2.8. The condition 39 can be enforced using the general formulation 5 for bounding time averages, applied to the augmented system 37. This yields the first part of proposition 2.8, whose sharpness is addressed by its second part. Proposition 2.8 can be applied only when the LE sum Mk+1M_{k+1} for certain kk is known to be negative on every trajectory. Such negativity can be confirmed by using 31, or its SOS relaxation that we introduce in section 2.4, to obtain a negative upper bound on Mk+1M_{k+1}.

Proposition 2.8.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}) and dynamics forward invariant on ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n}. Let Φk​(x,zk)\Phi_{k}(x,z^{k}) and ℓk​(x,zk)\ell_{k}(x,z^{k}) be as defined by 28 and 30, and likewise for Φk+1​(x,zk+1)\Phi_{k+1}(x,z^{k+1}) and ℓk+1​(x,zk+1)\ell_{k+1}(x,z^{k+1}).

1. 1. Suppose each trajectory x​(t)x(t) in ℬ\mathcal{B} is bounded for t≥0t\geq 0. Let the positive integer k≤nk\leq n be such that supx0∈ℬMk+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{k+1}(x_{0})<0. Then the global Lyapunov dimension dLd_{L} over ℬ\mathcal{B} is bounded above by

      dL≤infB∈[0,1]V∈C1(k+B)s.t.ΨB+f⋅Dx​V+ℓk⋅Dzk​V+ℓk+1⋅Dzk+1​V≤0for allx∈ℬ,zk∈𝕊nk−1,zk+1∈𝕊nk+1−1,d_{L}\leq\inf\limits_{\begin{subarray}{c}B\in[0,1]\\ V\in C^{1}~~\end{subarray}}(k+B)\quad\text{s.t.}\quad\begin{array}[t]{l}\Psi^{B}+f\cdot D_{x}V+\ell_{k}\cdot D_{z^{k}}V+\ell_{k+1}\cdot D_{z^{k+1}}V\leq 0\\[4.0pt] \text{for all}\quad x\in\mathcal{B},~z^{k}\in\mathbb{S}^{n_{k}-1},~z^{k+1}\in\mathbb{S}^{n_{k+1}-1},\end{array}

      (40)

      where V:ℬ×𝕊nk−1×𝕊nk+1−1→ℝV:\mathcal{B}\times\mathbb{S}^{{n_{k}}-1}\times\mathbb{S}^{n_{k+1}-1}\to\mathbb{R}, and Dx​VD_{x}V, Dzk​VD_{z^{k}}V and Dzk+1​VD_{z^{k+1}}V denote the gradients of V​(x,zk,zk+1)V(x,z^{k},z^{k+1}) with respect to each argument.

2. 1. Suppose ℬ\mathcal{B} is compact and f∈C2​(ℬ,ℝn)f\in C^{2}(\mathcal{B},\mathbb{R}^{n}). Let jj be the minimum admissible kk, as in 2.3 for dLd_{L}. Then, for all ε>0\varepsilon>0,
      | dL−ε​|supx0∈ℬMj+1​(x0)|−1≤j+Bε≤dL,d_{L}~-~\varepsilon\,\left|\sup_{x_{0}\in\mathcal{B}}M_{j+1}(x_{0})\right|^{-1}\leq j+B_{\varepsilon}\leq d_{L}, | (41) |
      | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |
      where

      Bε=infB∈[0,1]V∈C1Bs.t.ΨB+f⋅Dx​V+ℓj⋅Dzj​V+ℓj+1⋅Dzj+1​V≤εfor allx∈ℬ,zj∈𝕊nj−1,zj+1∈𝕊nj+1−1.B_{\varepsilon}=\inf\limits_{\begin{subarray}{c}B\in[0,1]\\ V\in C^{1}\hskip 7.0pt\end{subarray}}B\quad\text{s.t.}\quad\begin{array}[t]{l}\Psi^{B}+f\cdot D_{x}V+\ell_{j}\cdot D_{z^{j}}V+\ell_{j+1}\cdot D_{z^{j+1}}V\leq\varepsilon\\[4.0pt] \text{for all}\quad x\in\mathcal{B},~z^{j}\in\mathbb{S}^{n_{j}-1},~z^{j+1}\in\mathbb{S}^{n_{j+1}-1}.\end{array}

      (42)

Proof.

For the first part, the assumption that supx0∈ℬMk+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{k+1}(x_{0})<0 means that k≥jk\geq j where jj is as defined in 2.3. Applying 2.5, we therefore have

The constraint in 43 is equivalent to the Ψ¯B≤0\overline{\Psi}^{B}\leq 0 condition 39. To see this, we rearrange the constraint in 43 and use the expression 27 for MkM_{k} to obtain an equivalent inequality,

This inequality holds if and only if it holds without the suprema for all z0kz_{0}^{k} and z0k+1z_{0}^{k+1}, which is exactly the Ψ¯B≤0\overline{\Psi}^{B}\leq 0 condition 39. Writing this latter condition in place of the equivalent constraint in 43 gives

We want to replace the constraint in 45 with one that does not involve time averages. To do so we upper bound Ψ¯B\overline{\Psi}^{B} using the general formulation 5 for time averages by letting X=(x,zk,zk+1)X=(x,z^{k},z^{k+1}) on domain Ω=ℬ×𝕊nk−1×𝕊nk+1−1\Omega=\mathcal{B}\times\mathbb{S}^{n_{k}-1}\times\mathbb{S}^{n_{k+1}-1}, and letting F​(X)F(X) be the right-hand side of 37. This gives

The constraint of 43, which is equivalent to nonpositivity of the left-hand side of 46, can be strengthened by requiring there to exist V∈C1V\in C^{1} for which the inner supremum on the right-hand side is nonpositive. This gives the first part of the proposition.

For the second part, note that 43 is an equality when k=jk=j, according to 2.4. That equality can be expressed as

where equivalence of the constraints in 43 and 47 is as in the first part. To prove the second inequality in the claim 41, let BB be any value for which the constraint in 47 holds. It suffices to prove existence of VV for which constraint in the definition 42 of BεB_{\varepsilon} holds also. This will imply that the infimum over BB in 42 is at least as small as in 47.

Under the assumptions that ℬ\mathcal{B} is compact and f∈C2f\in C^{2}, the XX domain Ω\Omega is compact, and F​(X)F(X) defined as the right-hand side of 37 is C1C^{1} on an open region containing Ω\Omega. For such XX and FF, the general upper bound 5 is an equality [46], so in particular 46 is an equality. Choosing k=jk=j, equality in 47 implies that the constraint in 47 is equivalent to

This inequality implies that, for any ε>0\varepsilon>0, there exists V∈C1V\in C^{1} such that

(Such VV need not exist with ε=0\varepsilon=0 if the infimum over VV in the constraint of 46 is not attained.) Condition 49 is the same as the constraint in the minimization 42 defining BεB_{\varepsilon}, and this implies Bε≤BB_{\varepsilon}\leq B. We have shown that Bε≤BB_{\varepsilon}\leq B for any BB satisfying the constraint in 47, and thus dL−j≤Bεd_{L}-j\leq B_{\varepsilon}. This establishes the second inequality in 41.

To prove the first inequality in the claim 41, suppose that the inequality constraint in the definition 42 of BεB_{\varepsilon} is satisfied for a pair of BB and VV. Averaging this inequality along any trajectory and then maximizing over the initial directions z0kz_{0}^{k} and z0k+1z_{0}^{k+1} gives

which we divide by −μk+1>0-\mu_{k+1}>0 to find

Maximizing the first term over x0x_{0}, and bounding the second term above using μk+1​(x0)≤supx0∈ℬMk+1<0\mu_{k+1}(x_{0})\leq\sup_{x_{0}\in\mathcal{B}}M_{k+1}<0, we find

| (dL−j)−ε|supx0∈ℬMj+1​(x0)|≤B.(d_{L}-j)-\frac{\varepsilon}{\left|\sup_{x_{0}\in\mathcal{B}}M_{j+1}(x_{0})\right|}\leq B. | (52) | | ------------------------------------------------------------------------------------------------------------------------------------- | ---- |

Since 52 holds at all BB where the constraint set in the definition 42 of BεB_{\varepsilon} is not empty, minimizing over such BB gives 52 with BεB_{\varepsilon} on the right-hand side. This establishes the first inequality in 41, which completes the proof. ∎

Remark 2.9.

The second part of proposition 2.8 implies that

| dL≤j+Bε+ε​|supx0∈ℬMj+1​(x0)|−1d_{L}\leq j+B_{\varepsilon}+\varepsilon\,\left|\sup_{x_{0}\in\mathcal{B}}M_{j+1}(x_{0})\right|^{-1} | (53) | | ------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

at each ε>0\varepsilon>0, and that this upper bound converges to dLd_{L} as ε→0\varepsilon\to 0.

2.3 Shifted spectrum approach

In the preceding subsection, the possibly unbounded kk-vectors yky^{k} in the tangent space are captured by zkz^{k}, which is bounded because it is the projection of yky^{k} onto the unit sphere. In the alternative approach of the present subsection, we capture the dynamics of yky^{k} by letting

where the initial condition w0kw^{k}_{0} can be any kk-vector. This wk​(t)w^{k}(t) is bounded forward in time for sufficiently large BB. The ODE for wkw^{k} is the same as the ODE 15 for yky^{k}, except the spectrum of the Jacobian is shifted by BB. Together with the ODE 6 for x​(t)x(t) this gives the system

where

The following lemma asserts that if BB is large enough to prevent unbounded growth of wkw^{k}, then BB is an upper bound on the LE sum MkM_{k}, and the infimum over such BB is equal to MkM_{k}.

Lemma 2.10.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}) and dynamics forward invariant on ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n}. Let the initial condition x0x_{0} give a trajectory that is bounded for t≥0t\geq 0. Let wk​(t)w^{k}(t) satisfy the ODE 55 with initial condition w0kw_{0}^{k}. For each k∈{1,…,n}k\in\{1,\ldots,n\}:

1. 1. The leading LE sum Mk​(x0)M_{k}(x_{0}) is equivalent to

      Mk​(x0)=infB∈ℝBs.t.for each ​w0k∈ℝnk,wk​(t)​ is bounded for all ​t≥0.M_{k}(x_{0})=\inf_{B\in\mathbb{R}}B\quad\text{s.t.}\quad\text{for each }w^{k}_{0}\in\mathbb{R}^{n_{k}},~w^{k}(t)\text{ is bounded for all }t\geq 0.	(57)

2. 1. For any B>Mk​(x0)B>M_{k}(x_{0}) there exist K​(x0,B)>0K(x_{0},B)>0 and λ​(x0,B)>0\lambda(x_{0},B)>0 such that
      | |wk​(t)|≤K​e−λ​t​|w0k|for all t≥0 and ​w0k∈ℝnk.|w^{k}(t)|\leq Ke^{-\lambda t}|w^{k}_{0}|\quad\text{for all tgeq0t\geq 0tgeq0 and }w^{k}_{0}\in\mathbb{R}^{n_{k}}. | (58) |
      | -------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |
3. 1. If ℬ\mathcal{B} is compact and B>supx0∈ℬMk​(x0)B>\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0}), there exist K​(B)K(B) and λ​(B)\lambda(B) independent of x0x_{0} such that 58 holds for all x0∈ℬx_{0}\in\mathcal{B}.

Proof.

If there exists y0ky^{k}_{0} such that |yk​(t)|→∞|y^{k}(t)|\to\infty in finite time then, according to 54, no BB satisfies the condition in 57. Then both sides of 57 are ∞\infty, and the second part holds vacuously. Henceforth we assume that w​(t)w(t) exists for all t≥0t\geq 0.

We will first show that if BB satisfies the condition of 57 then Mk​(x0)≤BM_{k}(x_{0})\leq B. Substituting yk=eB​t​wky^{k}=e^{Bt}w^{k} into the time average on the right-hand side of 16 gives

| lim supt→∞1t​log⁡|yk​(t)|\displaystyle\limsup_{t\to\infty}\frac{1}{t}\log\left|y^{k}(t)\right| | =lim supt→∞1t​log⁡(|eB​t​wk​(t)|/|w0k|)\displaystyle=\limsup_{t\to\infty}\frac{1}{t}\log\left(\left|e^{Bt}w^{k}(t)\right|/|w^{k}_{0}|\right) | (59) | | --------------------------------------------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- | | =B+lim supt→∞1t​log⁡|wk​(t)|\displaystyle=B+\limsup_{t\to\infty}\frac{1}{t}\log\left|w^{k}(t)\right| | | | | ≤B,\displaystyle\leq B, | | |

where the final inequality follows from the boundedness of wk​(t)w^{k}(t). Taking the supremum of the left-hand side over y0ky^{k}_{0} gives Mk​(x0)≤BM_{k}(x_{0})\leq B.

Next we let B>Mk​(x0)B>M_{k}(x_{0}). Showing that 58 holds will prove the second part of the lemma, and it will complete the proof of the first part by ensuring boundedness of wk​(t)w^{k}(t). Without loss of generality, we assume |w0k|=1|w^{k}_{0}|=1. Since

| Mk​(x0)\displaystyle M_{k}(x_{0}) | ≥lim supt→∞1t​log⁡|yk​(t)|\displaystyle\geq\limsup_{t\to\infty}\frac{1}{t}\log{|y^{k}(t)|} | (60) | | ------------------------------------------------------------------------------------------------------ | --------------------------------------------------------------------------------------------------- | ---- | | =B+lim supt→∞1t​log⁡|wk​(t)|,\displaystyle=B+\limsup_{t\to\infty}\frac{1}{t}\log{|w^{k}(t)|}, | | |

we can choose some λ<(B−Mk​(x0))\lambda<(B-M_{k}(x_{0})) so that lim supt→∞1t​log⁡|wk​(t)|<−λ\limsup_{t\to\infty}\frac{1}{t}\log{|w^{k}(t)|}<-\lambda. Therefore by Grönwall’s inequality there exists a minimal T​(w0k)≥0T(w^{k}_{0})\geq 0 such that |wk​(t)|≤e−λ​t|w^{k}(t)|\leq e^{-\lambda t} for all t≥T​(w0k)t\geq T(w^{k}_{0}). By compactness, we can let

so that |wk​(t)|≤e−λ​t|w^{k}(t)|\leq e^{-\lambda t} for all t≥T∗t\geq T^{*} and all w0k∈𝕊nk−1w^{k}_{0}\in\mathbb{S}^{n_{k}-1}. Finally we let

| K=maxt∈[0,T∗]w0k∈𝕊nk−1⁡eλ​t​|wk​(t)|.K=\max_{\begin{subarray}{c}t\in[0,T^{*}]\\ w^{k}_{0}\in\mathbb{S}^{n_{k}-1}\hskip 1.0pt\end{subarray}}e^{\lambda t}|w^{k}(t)|. | (62) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

which completes the proof of the first two parts of the lemma.

For the third part, assume that Mk​(x0)M_{k}(x_{0}) is bounded above (otherwise the result is vacuously true), so that we can choose λ<B−supx0∈ℬMk​(x0)\lambda<B-\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0}). Let

| K=maxx0∈ℬt∈[0,T†]w0k∈𝕊nk−1⁡eλ​t​|wk​(t)|,K=\max_{\begin{subarray}{c}x_{0}\in\mathcal{B}\\ t\in[0,T^{\dagger}]\\ w^{k}_{0}\in\mathbb{S}^{n_{k}-1}\hskip 1.0pt\end{subarray}}e^{\lambda t}|w^{k}(t)|, | (63) | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

where T†=maxx0∈ℬ⁡T∗​(x0)T^{\dagger}=\max_{x_{0}\in\mathcal{B}}T^{*}(x_{0}) for the T∗​(x0)T^{*}(x_{0}) defined in 61. This maximum is valid because (t,x0,w0k)↦eλ​t​|wk​(t)|(t,x_{0},w^{k}_{0})\mapsto e^{\lambda t}|w^{k}(t)| is continuous on the compact set [0,T†]×ℬ×𝕊nk−1[0,T^{\dagger}]\times\mathcal{B}\times\mathbb{S}^{n_{k}-1}. With this value of λ\lambda and KK, we must then have 58 for every x0x_{0}, which completes the proof. ∎

Practical upper bounds on the LE sum MkM_{k} are found by combining 2.10 with the specialised Lyapunov theorem given as B.1 applied to the system 55. This yields the following proposition.

Proposition 2.11.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}) and dynamics forward invariant on ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n}. Let mkB​(x,wk)m^{B}_{k}(x,w^{k}) be defined as in 56.

1. 1. Suppose each trajectory x​(t)x(t) in ℬ\mathcal{B} is bounded for t≥0t\geq 0. For each k∈{1,…,n}k\in\{1,\ldots,n\}, the maximum leading LE sum MkM_{k} among trajectories in ℬ\mathcal{B} is bounded above by
      | supx0∈ℬMk​(x0)≤infV∈C1​(ℬ×ℝnk)Bs.t.V≥|wk|2f⋅Dx​V+mkB⋅Dwk​V≤0,\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0})\leq\inf_{V\in C^{1}(\mathcal{B}\times\mathbb{R}^{n_{k}})}B\quad\text{s.t.}\quad\begin{array}[t]{rl}V\geq&\hskip-7.0pt\left|w^{k}\right|^{2}\\ f\cdot D_{x}V+m^{B}_{k}\cdot D_{w^{k}}V\leq&\hskip-7.0pt0,\end{array} | (64) |
      | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |
      where the right-hand constraints are imposed for all (x,wk)∈ℬ×𝕊nk(x,w^{k})\in\mathcal{B}\times\mathbb{S}^{n_{k}}.
2. 1. Suppose ℬ\mathcal{B} is compact and f∈C2​(ℬ,ℝn)f\in C^{2}(\mathcal{B},\mathbb{R}^{n}). Then 64 holds with equality between the supremum and the infimum, and the Lyapunov functions may be restricted to the form V​(x,wk)=(wk)𝖳​U​(x)​wkV(x,w^{k})=(w^{k})^{\mathsf{T}}U(x)w^{k} with U∈C1​(ℬ,ℝnk×nk)U\in C^{1}(\mathcal{B},\mathbb{R}^{n_{k}\times n_{k}}).

Proof.

For the first part, assume BB and VV satisfy the constraints on the right-hand side of 64. Letting A​(x)=D​f[k]​(x)−B​InkA(x)=Df^{[k]}(x)-BI_{n_{k}}, the first part of B.1 applies to the system 55 and guarantees boundedness of wk​(t)w^{k}(t) for t≥0t\geq 0. Boundedness of wk​(t)w^{k}(t) means that BB satisfies the constraint on the right-hand side of 2.10, so Mk​(x0)≤BM_{k}(x_{0})\leq B for all x0∈ℬx_{0}\in\mathcal{B}. This implies 64, so the first part of the proposition is proven.

For the second part, let B>supx0Mk​(x0)B>\sup_{x_{0}}M_{k}(x_{0}). By the third part of 2.10, there exist K>0K>0 and λ>0\lambda>0such that 135 holds for all x0∈ℬx_{0}\in\mathcal{B}. Observe also that A∈C1​(ℬ,ℝnk×nk)A\in C^{1}(\mathcal{B},\mathbb{R}^{n_{k}\times n_{k}}) since f∈C2f\in C^{2}, and that A​(x)A(x) is bounded uniformly in xx since ℬ\mathcal{B} is compact. Thus the second part of B.1 applies, so there exists V​(x,wk)=(wk)𝖳​U​(x)​wkV(x,w^{k})=(w^{k})^{\mathsf{T}}U(x)w^{k} satisfying the constraint in 64. This is true for any B>supx0Mk​(x0)B>\sup_{x_{0}}M_{k}(x_{0}), so 64 is an equality. ∎

This method gives bounds for the different sums of LEs, which can then be used with proposition 2.7 to bound the Lyapunov dimension. Analogously to proposition 2.8, a direct bound for the dimension is also possible with a modified shifted spectrum approach, using the augmented system

where we have defined

| SkB​(x,yk+1,w)=D​f[k]​(x)​w+B1−B​(yk+1)𝖳​D​f[k+1]​(x)​yk+1|yk+1|2​w.S_{k}^{B}(x,y^{k+1},w)=Df^{[k]}(x)w+\frac{B}{1-B}\frac{(y^{k+1})^{\mathsf{T}}Df^{[k+1]}(x)y^{k+1}}{|y^{k+1}|^{2}}w. | (66) | | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---- |

No shifting is needed in the yk+1y^{k+1} dynamics because this quantity decays exponentially on average when Mk+1<0M_{k+1}<0. The ww dynamics 66 shift the D​f[k]Df^{[k]} spectrum not by a negative constant, as in 56, but by a factor proportional to the instantaneous growth rate of yk+1y^{k+1}, which is negative on average.

Proposition 2.12.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}) and dynamics forward invariant on ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n} with each trajectory x​(t)x(t) bounded for t≥0t\geq 0. Let the positive integer k≤nk\leq n be such that supx0∈ℬMk+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{k+1}(x_{0})<0. Then the global Lyapunov dimension dLd_{L} over ℬ\mathcal{B} is bounded above by

| dL≤infV∈C1(k+B)s.t.V​(x,yk+1,w)≥|w|2f⋅Dx​V+D​f[k+1]​yk+1⋅Dyk​V+SkB⋅Dw​V≤0,d_{L}\leq\inf_{V\in C^{1}}(k+B)\quad\text{s.t.}\quad\begin{array}[t]{rl}V(x,y^{k+1},w)\geq&\hskip-7.0pt\left|w\right|^{2}\\ f\cdot D_{x}V+Df^{[k+1]}y^{k+1}\cdot D_{y^{k}}V+S_{k}^{B}\cdot D_{w}V\leq&\hskip-7.0pt0,\end{array} | (67) | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

where the right-hand constraints are imposed for all (x,yk+1,w)∈ℬ×ℝnk+1×ℝnk(x,y^{k+1},w)\in\mathcal{B}\times\mathbb{R}^{n_{k+1}}\times\mathbb{R}^{n_{k}}, and Dx​VD_{x}V, Dzk​VD_{z^{k}}V and Dzk+1​VD_{z^{k+1}}V denote the gradients of V​(x,yk+1,w)V(x,y^{k+1},w) with respect to each argument.

Proof.

We will prove that, for any BB and VV satisfying the constraints in 67, dL≤k+Bd_{L}\leq k+B. We will use 2.5 to do this. First note that supx0∈ℬMk+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{k+1}(x_{0})<0 implies that k≥jk\geq j where jj is the minimal integer of 2.3. Observe that the dynamics of ww in the system 65 are of the form d​wd​t=A​(x,yk+1)​w\frac{\mathrm{d}w}{\mathrm{d}t}=A(x,y^{k+1})w for some matrix AA, so we may apply B.1 to this system. This tells us that given any BB and VV satisfying the constraints, w​(t)w(t) is bounded for any initial conditions. We now suppose that we have such a VV and BB and show that 23 holds for all x0x_{0}.

Consider any initial conditions for 65. Let yk=|yk+1|−B1−B​wy^{k}=|y^{k+1}|^{\frac{-B}{1-B}}w, which combined with 66 implies that yky^{k} satisfies the usual tangent dynamics 15. We then have

| w𝖳​D​f[k]​w|w|2=(yk)𝖳​D​f[k]​yk|yk|2.\frac{w^{\mathsf{T}}Df^{[k]}w}{|w|^{2}}=\frac{(y^{k})^{\mathsf{T}}Df^{[k]}y^{k}}{|y^{k}|^{2}}. | (68) | | --------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

Since |w​(t)||w(t)| is bounded,

| 0\displaystyle 0 | ≥lim supT→∞1T​log⁡|w​(T)||w​(0)|\displaystyle\geq\limsup_{T\to\infty}\frac{1}{T}\log\frac{|w(T)|}{|w(0)|} | (69) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------- | ---- | | =limT→∞1T​∫0Tdd​t​log⁡|w​(t)||w​(0)|​d​t\displaystyle=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{\mathrm{d}}{\mathrm{d}t}\log\frac{|w(t)|}{|w(0)|}\mathrm{d}t | (70) | | | =lim supT→∞1T​∫0T12​(d​|w|2d​t|w|2)​dt\displaystyle=\limsup_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{1}{2}\left(\frac{\frac{\mathrm{d}|w|^{2}}{\mathrm{d}t}}{|w|^{2}}\right)\mathrm{d}t | (71) | | | =lim supT→∞1T​∫0T(w𝖳​(D​f)[k]​w|w|2+B1−B​(yk+1)𝖳​(D​f)[k+1]​yk+1|yk+1|2)​dt.\displaystyle=\limsup_{T\to\infty}\frac{1}{T}\int_{0}^{T}\left(\frac{w^{\mathsf{T}}(Df)^{[k]}w}{|w|^{2}}+\frac{B}{1-B}\frac{(y^{k+1})^{\mathsf{T}}(Df)^{[k+1]}y^{k+1}}{|y^{k+1}|^{2}}\right)\mathrm{d}t. | (72) | |

Using 68 and 26, this gives

| 0≤lim supt→∞1t​log⁡|yk​(t)|+B1−B​lim supt→∞1t​log⁡|yk+1​(t)|.0\leq\limsup_{t\to\infty}\frac{1}{t}\log{|y^{k}(t)|}+\frac{B}{1-B}\limsup_{t\to\infty}\frac{1}{t}\log{|y^{k+1}(t)|}. | (73) | | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---- |

Since this holds for all yk+1y^{k+1} and yky^{k}, by 27 we must have

which rearranges to

2.5 finishes the proof. ∎

2.4 Sum-of-squares relaxations

In sections 2.2 and 2.3 we have formulated upper bounds on LE sums and on Lyapunov dimension that involve pointwise inequality constraints on functions over sets or, equivalently, pointwise suprema of functions over sets. All of these bounds can be expressed as solutions to optimization problems with pointwise nonnegativity constraints. When all expressions are polynomial, nonnegativity can be strengthened into polynomial SOS conditions. This approach is also the basis of the previous paper [34], and we briefly recall the main concepts here.

In order for a function S:ℝn→ℝS:\mathbb{R}^{n}\to\mathbb{R} to satisfy S​(x)≥0S(x)\geq 0 for all X∈ℝnX\in\mathbb{R}^{n}, it suffices that S∈ΣS\in\Sigma, where Σ\Sigma denotes the set of polynomials that are representable as a sum of squares of other polynomials. The set of nonnegative polynomials strictly contains Σ\Sigma, and the relationship between these sets has been extensively studied [43]. Often one wants to enforce S​(x)≥0S(x)\geq 0 on a subset of ℝn\mathbb{R}^{n} without requiring nonnegativity on all of ℝn\mathbb{R}^{n}. For instance, one may be interested in a set Ω\Omega with a ‘semiagebraic’ definition in terms of polynomial equalities and inequalities,

where the gi∈ℝ​[X]g_{i}\in\mathbb{R}[X] and hi∈ℝ​[X]h_{i}\in\mathbb{R}[X] are given. Here ℝ​[X]\mathbb{R}[X] denotes the ring of polynomials in XX with real coefficients. There are various ways to define a set ΣΩ\Sigma^{\Omega} of polynomials whose nonnegativity on Ω\Omega is implied by SOS conditions. Here we choose a set of polynomials that is called the quadratic module generated by the polynomials that define Ω\Omega,

where all ℝ​[X]\mathbb{R}[X] denotes the set of all polynomial functions of XX with real coefficients. Any polynomial in ΣΩ\Sigma^{\Omega} must be nonnegative on Ω\Omega because σ0​(X)≥0\sigma_{0}(X)\geq 0 at every X∈ℝnX\in\mathbb{R}^{n}, while each σi​(X)​gi​(X)≥0\sigma_{i}(X)g_{i}(X)\geq 0 and each ρj​(X)​hj​(X)=0\rho_{j}(X)h_{j}(X)=0 at every X∈ΩX\in\Omega. When Ω=ℝn\Omega=\mathbb{R}^{n}, there are no gig_{i} or ρj\rho_{j}, so ΣΩ=Σ\Sigma^{\Omega}=\Sigma. In the preceding paper [34], the condition S∈ΣΩS\in\Sigma^{\Omega} was expressed in a different but equivalent way: by saying there exist σi\sigma_{i} and ρj\rho_{j} such that S−∑σi​gi−∑ρj​hjS-\sum\sigma_{i}g_{i}-\sum\rho_{j}h_{j} is SOS.

Computational formulations enforce membership in ΣΩ\Sigma^{\Omega} via a stronger constraint of membership in a finite-dimensional subspace. The simplest choice is to enforce membership in ΣνΩ\Sigma^{\Omega}_{\nu}, which the truncation of the quadratic module 77 where each σi\sigma_{i} and ρj\rho_{j} has polynomial degree no larger than ν\nu. This gives a hierarchy of computational formulations, where raising ν\nu cannot worsen results–and typically improves them–but makes computations more expensive. In practice, it is often useful to choose the spaces for σi\sigma_{i} and ρi\rho_{i} more carefully. In our computations of section 3, for instance, we impose symmetries and allow different maximum degrees in different variables. For simplicity in the present subsection, however, we write formulations in terms of the simple truncation ΣνΩ\Sigma^{\Omega}_{\nu}.

The upper bounds on LE sums and Lyapunov dimension formulated in sections 2.2 and 2.3 are stated as, or can be reformulated as, minimizations subject to pointwise inequalities on VV. When the ODE right-hand side f​(x)f(x) is polynomial, and when we restrict to auxiliary functions VV in some finite-dimensional space of polynomials, then the pointwise inequalities can be enforced using SOS conditions. This gives SOS optimization problems that are computationally tractable when they are not too large, and whose minima give upper bounds on the desired quantities. For the shifted spectrum approach of section 2.3, the SOS relaxations are immediate after choosing a polynomial space for VV and enforcing nonnegativity on some set Ω\Omega via membership in ΣΩ\Sigma^{\Omega}. For the sphere projection approach of section 2.2, the relevant minmax problems can be expressed as constrained minimizations, after which the SOS relaxations are simple.

2.4.1 Sphere projection

The sphere projection approach of section 2.2 is based on the general upper bound 5 on time averages, which can be expressed as

An SOS relaxation of the right-hand side gives an upper bound, stated in the first part of proposition 2.13 below. Furthermore, this upper bound is an equality if Ω\Omega is forward invariant and its semialgebraic specification 76 has the Archimedean property, as stated in the second part of proposition 2.13 below. The Archimedean property implies compactness of Ω\Omega and is only slightly stronger; any specification of compact Ω\Omega that is not Archimedean can be given this property by adding the polynomial gi​(X)=R2−|X|2g_{i}(X)=R^{2}-|X|^{2} with a constant RR that is large enough to not change Ω\Omega itself. A proof of the second part of proposition 2.13 can be found in [19, Theorem 1]. Briefly, since Ω\Omega is compact and forward invariant, equality in 78 is given by the main theorem of [46], whose applicability on compact manifolds is addressed after (8) in [34]. Then, the Archimedean property guarantees that ΣΩ\Sigma^{\Omega} contains B−Φ−F⋅D​VB-\Phi-F\cdot DV for each suboptimal BB, by Putinar’s Positivstellensatz [41, Lemma 4.1].

Proposition 2.13.

Let dd​t​X​(t)=F​(X​(t))\frac{\mathrm{d}}{\mathrm{d}t}X(t)=F(X(t)) with F∈ℝn​[X]F\in\mathbb{R}^{n}[X], and let Φ∈ℝ​[X]\Phi\in\mathbb{R}[X]. Let the dynamics be forward invariant on a semialgebraic set Ω⊂ℝn\Omega\subset\mathbb{R}^{n}, whose specification defines the quadratic module ΣΩ\Sigma^{\Omega}.

1. 1. If each trajectory X​(t)X(t) in Ω\Omega is bounded for t≥0t\geq 0, then

      maxX0∈Ω⁡Φ¯≤infV∈ℝ​[X]Bs.t.B−Φ−F⋅D​V∈ΣΩ.\max_{X_{0}\in\Omega}\overline{\Phi}\leq\inf_{V\in\mathbb{R}[X]}B\quad\text{s.t.}\quad B-\Phi-F\cdot DV\in\Sigma^{\Omega}.	(79)

2. 1. If Ω\Omega is compact and its semialgebraic specification is Archimedean, then,

      maxX0∈Ω⁡Φ¯=infV∈ℝ​[X]Bs.t.B−Φ−F⋅D​V∈ΣΩ.\max_{X_{0}\in\Omega}\overline{\Phi}=\inf_{V\in\mathbb{R}[X]}B\quad\text{s.t.}\quad B-\Phi-F\cdot DV\in\Sigma^{\Omega}.	(80)

All of the upper bounds from section 2.2 can be relaxed into SOS optimization problems by the same reasoning that leads from 5 to 79. Proposition 2.15 below is an SOS relaxation of proposition 2.6 for bounding LE sums. Computations giving such bounds also imply bounds on Lyapunov dimension via proposition 2.7, but potentially sharper bounds on dimension can be computed directly using proposition 2.16, which is an SOS relxation of proposition 2.8. In the proposition statements, ℝn​[x]\mathbb{R}^{n}[x] denotes the space of polynomials in xx taking vector values in ℝn\mathbb{R}^{n}, and the x​(t)x(t) domain is a semialgebraic set,

This is understood to include the ℬ=ℝn\mathcal{B}=\mathbb{R}^{n} case when there are no gig_{i} or hjh_{j}. The k=1k=1 special case of proposition 2.15 is essentially Proposition 2 of [34], although the SOS conditions are more explicit in [34], and the restriction to finite-dimensional polynomial spaces is more explicit here.

Remark 2.14.

In the following propositions 2.15, 2.16, 2.17 and 2.18, for fixed k∈{1,…,n}k\in\{1,\ldots,n\} and polynomial degrees ν\nu and ν′\nu^{\prime}, the computation of the resulting upper bound Ci​(k,ν,ν′)C_{i}(k,\nu,\nu^{\prime}) is an SOS program that can be solved using standard software (see section 3). For certain combinations (k,ν,ν′)(k,\nu,\nu^{\prime}), there will not exist any VV in the chosen space of polynomials that satisfy the SOS constraints, in which case we understand the infimum over bounds to be +∞+\infty, meaning that we obtain no upper bound. When bounding Lyapunov dimension, choosing k>dL−1k>d_{L}-1 is necessary for the constraint set to not be empty. In all cases, raising the polynomial degrees ν\nu and ν′\nu^{\prime} cannot shrink the constraint set and will often enlarge it.

Proposition 2.15.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈ℝn​[x]f\in\mathbb{R}^{n}[x], for which a semialgebraic set ℬ\mathcal{B} is forward invariant. Let ΣΩ\Sigma^{\Omega} be the quadratic module 77 for the set Ω=ℬ×𝕊nk\Omega=\mathcal{B}\times\mathbb{S}^{n_{k}} whose semialgebraic specification includes the polynomial constraints on xx in 81 and the equality |zk|−1=0\big|z^{k}\big|-1=0. Let Φk​(x,zk)\Phi_{k}(x,z^{k}) and ℓk​(x,zk)\ell_{k}(x,z^{k}) be defined as in 28 and 30, respectively.

1. 1. Suppose that each trajectory x​(t)x(t) in ℬ\mathcal{B} is bounded for t≥0t\geq 0. For each k∈{1,…,n}k\in\{1,\ldots,n\}, and each pair of polynomial degrees (ν,ν′)(\nu,\nu^{\prime}), the maximum leading LE sum MkM_{k} among trajectories in ℬ\mathcal{B} is bounded above by

      supx0∈ℬMk​(x0)≤C1​(k,ν,ν′),\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0})\leq C_{1}(k,\nu,\nu^{\prime}),	(82)
      where
      C1​(k,ν,ν′)=infV∈ℝ​[x,zk]νBs.t.B−Φk−f⋅Dx​V−ℓk⋅Dzk​V∈Σν′Ω.C_{1}(k,\nu,\nu^{\prime})=\inf_{V\in\mathbb{R}[x,z^{k}]_{\nu}}B\quad\text{s.t.}\quad B-\Phi_{k}-f\cdot D_{x}V-\ell_{k}\cdot D_{z^{k}}V\in\Sigma^{\Omega}_{\nu^{\prime}}.	(83)
      -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------	----

2. 1. Suppose further that the invariant set ℬ\mathcal{B} is compact, and its semialgebraic specification is Archimedean. Then, for each kk,

      supx0∈ℬMk​(x0)=supν,ν′∈ℤ≥C1​(k,ν,ν′).\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0})=\sup_{\nu,\nu^{\prime}\in\mathbb{Z}_{\geq}}C_{1}(k,\nu,\nu^{\prime}).	(84)

Proof.

The first part is an SOS relaxation of the upper bound 31 on MkM_{k} for each kk. Equivalently, it is a particular case of 79 with XX, FF and Φ\Phi as described above in proposition 2.6, where VV is further constrained by having degree no larger than ν\nu, and where σi\sigma_{i} and ρj\rho_{j} in the quadratic module representation 77 have degrees no larger than ν′\nu^{\prime}. The second part is the analogous special case of 80. This equality is applicable because the ℬ\mathcal{B} is specified as an Archimedean subset of ℝn\mathbb{R}^{n}, and so adding |zk|−1=0|z^{k}|-1=0 gives an Archimedean specification of Ω\Omega in ℝn+nk\mathbb{R}^{n+n_{k}}. All polynomial degrees are admitted in 80, which is equivalent to the supremum over ν\nu and ν′\nu^{\prime} in 84. ∎

Proposition 2.16.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈ℝn​[x]f\in\mathbb{R}^{n}[x], for which a semialgebraic set ℬ\mathcal{B} is forward invariant, and let the integer k≤nk\leq n be such that supx0∈ℬMk+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{k+1}(x_{0})<0. Let ΣΩ\Sigma^{\Omega} be the quadratic module 77 for the set Ω=ℬ×𝕊nk−1×𝕊nk+1−1\Omega=\mathcal{B}\times\mathbb{S}^{n_{k}-1}\times\mathbb{S}^{n_{k+1}-1} whose semialgebraic specification includes the polynomial constraints on xx in 81, and the equalities |zk|−1=0\big|z^{k}\big|-1=0 and |zk+1|−1=0\big|z^{k+1}\big|-1=0. Let Φk​(x,zk)\Phi_{k}(x,z^{k}) and ℓk​(x,zk)\ell_{k}(x,z^{k}) be defined as in 28 and 30, respectively, and likewise for Φk+1​(x,zk+1)\Phi_{k+1}(x,z^{k+1}) and ℓk+1​(x,zk+1)\ell_{k+1}(x,z^{k+1}).

1. 1. Suppose that each trajectory x​(t)x(t) in ℬ\mathcal{B} is bounded for t≥0t\geq 0. For each pair of polynomial degrees (ν,ν′)(\nu,\nu^{\prime}), the the maximum Lyapunov dimension dLd_{L} among trajectories in ℬ\mathcal{B} is bounded above by

      dL≤k+C2​(k,ν,ν′,0)d_{L}\leq k+C_{2}(k,\nu,\nu^{\prime},0)	(85)
      where
      C2​(k,ν,ν′,ε)=infB∈[0,1]V∈ℝ​[x,zk,zk+1]νBs.t.ε−[ΨB+f⋅Dx​V+ℓk⋅Dzk​V+ℓk+1⋅Dzk+1​V]∈Σν′Ω.C_{2}(k,\nu,\nu^{\prime},\varepsilon)=\hskip-6.0pt\inf\limits_{\begin{subarray}{c}B\in[0,1]\\ V\in\mathbb{R}[x,z^{k},z^{k+1}]_{\nu}\end{subarray}}\hskip-10.0ptB\quad\text{s.t.}\quad\varepsilon-[\Psi^{B}+f\cdot D_{x}V+\ell_{k}\cdot D_{z^{k}}V+\ell_{k+1}\cdot D_{z^{k+1}}V]\in\Sigma^{\Omega}_{\nu^{\prime}}.	(86)
      ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------	----

2. 1. Suppose further that the invariant set ℬ\mathcal{B} is compact, and its semialgebraic specification is Archimedean. Let jj be the minimum admissible kk, as in 2.3 for dLd_{L}. Then,

      dL=j+infν,ν′∈ℤ≥ε>0C2​(j,ν,ν′,ε).d_{L}=j~+\inf_{\begin{subarray}{c}\nu,\nu^{\prime}\in\mathbb{Z}{{}_{\geq}}\\ \varepsilon>0\end{subarray}}C_{2}(j,\nu,\nu^{\prime},\varepsilon).	(87)

Proof.

The first part follows from the upper bound 40 on dLd_{L} because the SOS constraint in 86 with ε=0\varepsilon=0 implies the nonpositivity constraint on the right-hand side of 40. For the second part, we first establish upper and lower bounds on C2C_{2} in terms of dLd_{L}. For the lower bound, note that the constraint in 86 implies the constraint in the definition 42 of BεB_{\varepsilon}. This implies Bε≤C2​(j,ν,ν′,ε)B_{\varepsilon}\leq C_{2}(j,\nu,\nu^{\prime},\varepsilon), with which the first inequality in 41 gives

| dL−ε​|supx0∈ℬMj+1​(x0)|−1≤j+C2​(j,ν,ν′,ε)d_{L}~-~\varepsilon\,\left|\sup_{x_{0}\in\mathcal{B}}M_{j+1}(x_{0})\right|^{-1}\leq j+C_{2}(j,\nu,\nu^{\prime},\varepsilon) | (88) | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

for all ε>0\varepsilon>0 and ν,ν′∈ℤ≥\nu,\nu^{\prime}\in\mathbb{Z}_{\geq}. Minimizing over these parameters gives

The infimum will be approached or attained as ν,ν′→∞\nu,\nu^{\prime}\to\infty and ε→0+\varepsilon\to 0^{+}.

For the upper bound on C2C_{2}, note that since ℬ\mathcal{B} has an Archimedean specification, so does Ω\Omega. Thus we can apply the general SOS upper bounds with equality 80, in particular with F​(X)F(X) being the right-hand side of the ODE system 37, with Φ​(X)\Phi(X) being ΨB\Psi^{B}, and with k=jk=j. This gives

Note that the variable called BB in 80 is ε\varepsilon in 90, where BB has a different meaning. Let BB be such that the constraint holds in expression 47 for dLd_{L}. This constraint is equivalent to both sides of 90 being nonpositive, and nonpositivity of the right-hand infimum in 90 implies that there exists polynomial VV satisfying the SOS constraint for any ε>0\varepsilon>0. Since such VV exists for any B>dL−jB>d_{L}-j ,

for all ε>0\varepsilon>0. Restricting the right-hand minimization to VV of degree ν\nu and weights σi\sigma_{i} and ρj\rho_{j} of degree ν′\nu^{\prime} gives the definition of C2​(j,ν,ν′,ε)C_{2}(j,\nu,\nu^{\prime},\varepsilon). Thus, for all ε>0\varepsilon>0,

(This inequality need not hold with finite ν\nu and ν′\nu^{\prime}.) This upper bound of dLd_{L}, combined with the lower bound of dLd_{L} 89, established the claim 87 of the second part and completes the proof. ∎

The minimization problems defining C1C_{1} and C2C_{2} for each (k,ν,ν′)(k,\nu,\nu^{\prime}) in 82 of proposition 2.15 are SOS programs. That is, polynomials of finite degree are constrained to be SOS, and these polynomials as well as the optimization objective are affine in the optimization parameters. Here the optimization parameters include BB as well as the coefficients of VV in whatever basis it is represented. Such SOS programs are computationally tractable so long as the dimension of XX and the degree of polynomials in the SOS constraints are not too large. Increasing ν\nu and ν′\nu^{\prime} gives a hierarchy of SOS programs with increasing computational cost whose minima cannot increase. Typically these minima decrease as ν\nu and ν′\nu^{\prime} are raised, and the second part of each theorem guarantees convergence to the sharp bound under mild conditions. More generally, a hierarchy of SOS programs can be defined by choosing different sequences of finite-dimensional polynomial spaces, rather than simply truncating by maximum degrees of ν\nu and ν′\nu^{\prime}. This often includes imposing symmetries on each polynomial subspace due to symmetries of f​(x)f(x), as explained in section 2.5 below.

2.4.2 Shifted spectrum

Next we give the SOS formulations based on the shifted spectrum approach. Proposition 2.17 below is an SOS relaxation of proposition 2.11 for bounding LE sums. Computations giving such bounds also imply bounds on Lyapunov dimension via proposition 2.7, but potentially sharper bounds on dimension can be computed directly using proposition 2.18, which is an SOS relaxation of proposition 2.12. Unlike our other methods, we do not prove that proposition 2.18 gives sharp results in general.

Proposition 2.17.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈ℝn​[x]f\in\mathbb{R}^{n}[x], for which a semialgebraic set ℬ\mathcal{B} is forward invariant. Let mkB​(x,wk)m^{B}_{k}(x,w^{k}) be defined as in 56. Let ΣΩ\Sigma^{\Omega} be the quadratic module 77 for the set Ω=ℬ×ℝnk\Omega=\mathcal{B}\times\mathbb{R}^{n_{k}} whose semialgebraic specification includes the same constraints 81 defining ℬ\mathcal{B}.

1. 1. For each k∈{1,…,n}k\in\{1,\ldots,n\}, and each pair of polynomial degrees (ν,ν′)(\nu,\nu^{\prime}), the maximum leading LE sum MkM_{k} among trajectories in ℬ\mathcal{B} is bounded above by

      supx0∈ℬMk​(x0)≤C3​(k,ν,ν′),\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0})\leq C_{3}(k,\nu,\nu^{\prime}),	(93)
      where
      C3​(k,ν,ν′)=infV∈ℝ​[x,wk]νBs.t.V−|wk	2∈Σν′Ω−[f⋅Dx​V+mkB⋅Dwk​V]∈Σν′Ω.C_{3}(k,\nu,\nu^{\prime})=\inf_{V\in\mathbb{R}[x,w^{k}]_{\nu}}B\quad\text{s.t.}\quad\begin{array}[t]{r}V-
      -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------	----

2. 1. Suppose that the invariant set ℬ\mathcal{B} is compact, and that its semialgebraic specification is Archimedean. Then, for each kk,

      supx0∈ℬMk​(x0)=infν,ν′∈ℤ≥C3​(k,ν,ν′),\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0})=\inf_{\nu,\nu^{\prime}\in\mathbb{Z}_{\geq}}C_{3}(k,\nu,\nu^{\prime}),	(95)
      where the infimum in 94 is over V​(x,wk)=(wk)𝖳​U​(x)​wkV(x,w^{k})=(w^{k})^{\mathsf{T}}U(x)w^{k} with U∈ℝnk×nk​[x]U\in\mathbb{R}^{n_{k}\times n_{k}}[x] a symmetric polynomial matrix, and where σi​(x,wk)\sigma_{i}(x,w^{k}) and ρj​(x,wk)\rho_{j}(x,w^{k}) in the SOS representation 77 are likewise quadratic in wkw^{k}.

Proof.

The first part follows from the upper bound 64 on MkM_{k} because the SOS constraints in 94 imply the nonnegativity constraints in 64 from proposition 2.11. For the second part, let B1B_{1} and B2B_{2} be arbitrary with B1>B2>supx0∈ℬMk​(x0)B_{1}>B_{2}>\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0}). It suffices to show that there exists a polynomial matrix U​(x)U(x) such that V​(x,wk)=(wk)𝖳​U​(x)​wkV(x,w^{k})=(w^{k})^{\mathsf{T}}U(x)w^{k} satisfies the SOS constraints of 94 for B=B1B=B_{1} with no restrictions on polynomial degree. Since B2>supx0∈ℬMk​(x0)B_{2}>\sup_{x_{0}\in\mathcal{B}}M_{k}(x_{0}), the second part of proposition 2.11 guarantees the existence of a symmetric U∈C1​(ℝn,ℝnk×nk)U\in C^{1}(\mathbb{R}^{n},\mathbb{R}^{n_{k}\times n_{k}}) for which VV satisfies the pointwise inequalities in 64 with B=B2B=B_{2}. We will show that this implies existence of the desired polynomial UU under the present assumptions.

The pointwise inequalities in 64, with B=B2B=B_{2} in the V=(wk)𝖳​U​wkV=(w^{k})^{\mathsf{T}}Uw^{k} case, can be written in the form (wk)𝖳​Z^i​(x)​wk≥0(w^{k})^{\mathsf{T}}\hat{Z}_{i}(x)w^{k}\geq 0, where

with InkI_{n_{k}} denoting the nk×nkn_{k}\times n_{k} identity matrix. The pointwise inequalities on ℬ×ℝnk\mathcal{B}\times\mathbb{R}^{n_{k}} are equivalent to the condition that both Z^i​(x)⪰0\hat{Z}_{i}(x)\succeq 0 for all x∈ℬx\in\mathcal{B}. In fact, by scaling UU we can satisfy both of these equations with the additional restriction that Z^1≻0\hat{Z}_{1}\succ 0 strictly on ℬ\mathcal{B}. Let

which are the same definitions with B2B_{2} replaced by B1B_{1}. Since B1>B2B_{1}>B_{2}, we now have both Zi​(x)≻0Z_{i}(x)\succ 0 strictly for all x∈ℬx\in\mathcal{B}.

Because there exists U∈C1U\in C^{1} such that both Zi​(x)≻0Z_{i}(x)\succ 0 on ℬ\mathcal{B}, and ℬ\mathcal{B} is compact, there also exists a polynomial matrix U∈ℝnk×nk​[x]U\in\mathbb{R}^{n_{k}\times n_{k}}[x] for which both Zi​(x)≻0Z_{i}(x)\succ 0 on ℬ\mathcal{B}. This follows from a strengthened version of the Stone-Weierstrass approximation theorem [25]. This polynomial UU satisfies the SOS constraints of 94 due to ℬ\mathcal{B} being Archimedean and ZiZ_{i} being strictly positive definite on ℬ\mathcal{B}. Recalling that ℬ\mathcal{B} is the set where all gi​(x)≥0g_{i}(x)\geq 0 and hj​(x)=0h_{j}(x)=0, applying a matrix Positivstellensatz [44, theorem 2] guarantees existence of polynomial matrices Si​(x)S_{i}(x) and Ti​(x)T_{i}(x), with the Si​(x)S_{i}(x) being squares of other polynomial matrices, such that

and likewise for Z2Z_{2}. The existence of such representations implies that both (wk)𝖳​Zi​wk(w^{k})^{\mathsf{T}}Z_{i}w^{k} belong to the quadratic module ΣΩ\Sigma^{\Omega}. Therefore V=(wk)𝖳​U​wkV=(w^{k})^{\mathsf{T}}Uw^{k} satisfies the constraints of 94 with B=B1B=B_{1}, and the second part of the proposition is proven. ∎

Proposition 2.18.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f∈ℝn​[x]f\in\mathbb{R}^{n}[x], for which a semialgebraic set ℬ\mathcal{B} is forward invariant on ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n} with each trajectory x​(t)x(t) bounded for t≥0t\geq 0. Let the positive integer k≤nk\leq n be such that supx0∈ℬMk+1​(x0)<0\sup_{x_{0}\in\mathcal{B}}M_{k+1}(x_{0})<0. Let ΣΩ\Sigma^{\Omega} be the quadratic module 77 for the set Ω=ℬ×ℝnk+1×ℝnk\Omega=\mathcal{B}\times\mathbb{R}^{n_{k+1}}\times\mathbb{R}^{n_{k}} whose semialgebraic specification includes the polynomial constraints on xx in 81. Let each trajectory x​(t)x(t) in ℬ\mathcal{B} be bounded for t≥0t\geq 0. For each k∈{1,…,n}k\in\{1,\ldots,n\}, and each pair of polynomial degrees (ν,ν′)(\nu,\nu^{\prime}), the the maximum Lyapunov dimension dLd_{L} among trajectories in ℬ\mathcal{B} is bounded above by

where

| C4​(k,ν,ν′)=infB∈[0,1]V∈ℝ​[x,yk+1,w]νBs.t.V−|w|2∈Σν′Ω(1−B)​|yk+1|2​SkB∈Σν′Ω.C_{4}(k,\nu,\nu^{\prime})=\inf_{\begin{subarray}{c}B\in[0,1]\\ V\in\mathbb{R}[x,y^{k+1},w]_{\nu}\end{subarray}}B\quad\text{s.t.}\quad\begin{array}[t]{r}V-\left|w\right|^{2}\in\Sigma^{\Omega}_{\nu^{\prime}}\phantom{,}\\ (1-B)\left|y^{k+1}\right|^{2}S_{k}^{B}\in\Sigma^{\Omega}_{\nu^{\prime}}.\end{array} | (102) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ----- |

2.5 Symmetry exploitation

If an ODE being studied using auxiliary function methods has symmetry, then in many cases the same symmetry can be imposed on the auxiliary functions, leading to symmetries in SOS programs that can be exploited computationally by block diagonalizing the corresponding semidefinite programs (SDPs) [11]. More precisely, an ODE system 1 with right-hand side F​(X)F(X) is invariant under a map Λ:ℝn→ℝn\Lambda:\mathbb{R}^{n}\to\mathbb{R}^{n} if and only if FF is Λ\Lambda-equivariant, meaning that F​(Λ​X)=Λ​F​(X)F(\Lambda X)=\Lambda F(X) for all XX. In such cases, X​(t)X(t) solves the ODE if and only if Λ​X​(t)\Lambda X(t) does. The trajectory itself may be Λ\Lambda-invariant, meaning that X​(t)=Λ​X​(t)X(t)=\Lambda X(t), or it may not be. For any ODE, the set of such transformations forms a group. We focus on the common situation where each Λ\Lambda is an orthogonal linear transformation on ℝn\mathbb{R}^{n}, meaning the group of such Λ\Lambda is a subgroup of the orthogonal group O​(n)O(n). (In some contexts the set of Λ\Lambda would be called the linear representation of a group, rather than a group itself, but we do not draw this distinction.) In many formulations of auxiliary function methods for ODEs, if the ODE is invariant under Λ\Lambda then the results are provably unchanged if the optimization is over the restricted set of VV which are Λ\Lambda-invariant, meaning V​(X)=V​(Λ​X)V(X)=V(\Lambda X) for all XX. In the context of bounding time averages, see the appendix of [34] for a proof that V​(X)V(X) can be taken as Λ\Lambda-invariant under certain conditions. For our shifted spectrum methods, the Lyapunov function in B.1 can be taken as Λ\Lambda-invariant also. We omit the proof, which closely follows the argument in [34].

In order to apply such results in the context of bounding LE sums and Lyapunov dimension, we must determine the symmetries of ODEs on the augmented state spaces for (x,zk)(x,z^{k}) and (x,wk)(x,w^{k}), governed by 29 and 55. First, regardless of whether the xx ODE has any symmetry, the augmented ODEs have a sign symmetry due to linearity of the yky^{k} dynamics. These sign symmetries are summarized by proposition 2.19. Second, if the right-hand side f​(x)f(x) of the xx ODE is equivariant under an orthogonal transformation, this induces equivariance on the right-hand sides of the augmented ODEs. This is summarized by proposition 2.20, which invokes the kt​hk^{th} multiplicative compound Λ(k)\Lambda^{(k)} of Λ\Lambda (cf. appendix A). Note that InI_{n} denotes the n×nn\times n identity matrix, and diag(A,…,B)\operatorname*{diag}(A,\dots,B) is a block diagonal matrix with square blocks A,…,BA,\dots,B.

Proposition 2.19.

Let dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f:ℬ→ℝnf:\mathcal{B}\to\mathbb{R}^{n}. The augmented ODEs are invariant under the following orthogonal transformations.

1. 1. The ODE system 29 is invariant under diag(In,−Ink)\operatorname*{diag}(I_{n},-I_{n_{k}}).
2. 1. The ODE system 37 is invariant under diag(In,±Ink,∓Ink+1)\operatorname*{diag}(I_{n},\pm I_{n_{k}},\mp I_{n_{k+1}}).
3. 1. The ODE system 55 is invariant under diag(In,−Ink)\operatorname*{diag}(I_{n},-I_{n_{k}}).
4. 1. The ODE system 65 is invariant under diag(In,±Ink+1,∓Ink)\operatorname*{diag}(I_{n},\pm I_{n_{k+1}},\mp I_{n_{k}}).

Proof.

All four parts have very similar proofs, so we prove only the last. The claim amounts to the right-hand side of 65 being equivariant under negation of yk+1y^{k+1}. For f​(x)f(x) this statement is trivial. The requirement that D​f[k+1]​(x)​(−yk+1)=−D​f[k+1]​(x)​yk+1Df^{[k+1]}(x)(-y^{k+1})=-Df^{[k+1]}(x)y^{k+1} is also immediate. The final requirement that SkB​(x,−yk+1,w)=SkB​(x,yk+1,w)S_{k}^{B}(x,-y^{k+1},w)=S_{k}^{B}(x,y^{k+1},w) holds because yk+1y^{k+1} appears quadratically in both the numerator and denominator in 66. ∎

Proposition 2.20.

Let ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n} be invariant under Λ∈O​(n)\Lambda\in O(n). If dd​t​x​(t)=f​(x​(t))\frac{\mathrm{d}}{\mathrm{d}t}x(t)=f(x(t)) with f:ℬ→ℝnf:\mathcal{B}\to\mathbb{R}^{n} is invariant under Λ\Lambda, then the augmented ODEs are invariant under the following orthogonal transformations.

1. 1. The ODE system 29 is invariant under diag(Λ,Λ(k))\operatorname*{diag}(\Lambda,\Lambda^{(k)}).
2. 1. The ODE system 37 is invariant under diag(Λ,Λ(k),Λ(k+1))\operatorname*{diag}(\Lambda,\Lambda^{(k)},\Lambda^{(k+1)}).
3. 1. The ODE system 55 is invariant under diag(Λ,Λ(k))\operatorname*{diag}(\Lambda,\Lambda^{(k)}).
4. 1. The ODE system 65 is invariant under diag(Λ,Λ(k+1),Λ(k))\operatorname*{diag}(\Lambda,\Lambda^{(k+1)},\Lambda^{(k)}).

Proof.

By assumption f​(x)f(x) is equivariant under the orthogonal transformation Λ\Lambda, meaning f​(Λ​x)=Λ​f​(x)f(\Lambda x)=\Lambda f(x). The multiplicative compound Λ(k)\Lambda^{(k)} is orthogonal by proposition A.2, thus all transformations in the proposition are orthogonal. The claim is proved by showing that the right-hand sides of the augmented ODEs are equivariant under the claimed transformations. All four parts are very similar, so we prove only the first. For this we must show that ℓk​(Λ​x,Λ(k)​zk)=Λ(k)​ℓk​(x,zk)\ell_{k}(\Lambda x,\Lambda^{(k)}z^{k})=\Lambda^{(k)}\ell_{k}(x,z^{k}). It is sufficient to show that

so that

by orthogonality of Λ(k)\Lambda^{(k)}. Oeri and Goluskin [34] showed thatD​f​(Λ​x)=Λ​D​f​(x)​Λ𝖳Df(\Lambda x)=\Lambda Df(x)\Lambda^{\mathsf{T}} . Taking the kkth additive compound of this equation gives, via proposition A.3,

Orthogonality then gives 103. The proofs of the other parts follow similarly, noting that 103 remains valid with yky^{k} or wkw^{k} in place of zkz^{k}. ∎

Combining both propositions gives any finite symmetry group for the augmented systems. For example, if ℬ\mathcal{B} is invariant and ff is equivariant under a group generated by {Λ1,…,Λj}\{\Lambda_{1},\dots,\Lambda_{j}\} then the right-hand sides of the augmented systems 29 and 55 are equivariant under the group generated by

and we can restrict to VV that are invariant under the same transformations. This can be enforced by representing VV in a symmetry-invariant basis. For sign symmetries this can be done with monomial bases, but general rotations require more sophisticated choices of symmetry-invariant polynomial basis. Note that for the shifted spectrum method proposition 2.17 it is sufficient to consider VV that are is quadratic in wkw^{k}, which automatically gives invariance of VV with respect to the symmetry of proposition 2.19.

3 Computational examples

We now present computational examples in which our methods based on SOS optimization are applied to particular ODEs. It has been conjectured that, for all systems with equilibria embedded in chaotic attractors, the supremum in the formula for the Lyapunov dimension 18 is attained on an equilibrium [9, 22]. This indeed is the case for many simple systems, including Lorenz’s 1963 model [23]. In order to demonstrate the power of our methods, we focus on two example systems for which we know there are no equilibria.

For all of our examples, the code to replicate the results, implemented in Julia, is available at https://github.com/jeremypparker/Parker_Goluskin_LEs. We used SumOfSquares.jl [20] to reformulate SOS optimization problems as SDPs, along with SymbolicWedderburn.jl [15] for symmetry exploitation. To solve the resulting SDPs, we primarily used Mosek [1]. Certain SDPs were also solved using Clarabel.jl [14] and Hypatia.jl [5], which gave very similar results with longer runtimes. In applications of the sphere projection approach, the bounds BB can be minimized as the objective of the SDPs, and the bounds we report are the numerical optima returned by Mosek. In applications of the shifted spectrum approach, the value of BB must be fixed during each SDP solution because it multiplies coefficients of VV in the SOS constraints, and the coefficients to be optimized over can appear only linearly. In such computations, if there exists a polynomial VV satisfying the SOS constraints, then the SDP will have a nonempty constraint set, and the SDP solver should report feasibility. We thus perform a bisection search in BB to find the boundary between large BB values where the SDPs are feasible and small BB values where they are infeasible. Close to this boundary the numerical conditioning of the SDPs becomes increasingly poor, making it difficult to precisely identify the infimum over admissible BB. Confirming that near-optimal bounds are valid calls for rigorous and/or multiple-precision computations [12, 37], but for simplicity we have used only double precision arithmetic in our examples. We call a bound ‘feasible’ if Mosek reports a dual objective value of less than 10−610^{-6}, ‘infeasible’ if the final dual objective value is greater than 11 and ‘indeterminate’ otherwise, but we stress that these computations are not rigorously validated.

To confirm that accuracy of our computed bounds, we computed various period orbits and the LE sums along each of these orbits. We used the Tsit5 solver from the OrdinaryDiffEq.jl package [42] to solve integrate the ODE systems, find periodic orbits and compute LEs. For a given periodic orbit of period TT, the LEs μk\mu_{k} were found using the eigenvalues of the Jacobian D​ϕTD\phi^{T}. To compute the sums MkM_{k}, we simply sum the μk\mu_{k} in descending order.

3.1 The Duffing Oscillator

The Duffing oscillator is a second-order nonautonomous ODE governing amplitude of oscillations that are damped and sinusoidally forced. Different versions exist, but we choose

with parameters (α,β,γ,δ,ω)=(1,−1,0.3,0.2,1)(\alpha,\beta,\gamma,\delta,\omega)=(1,-1,0.3,0.2,1). We numerically estimate the Lyapunov exponents on the chaotic attractor to be approximately 0.1580.158 and −0.358-0.358, which gives a local Lyapunov dimension of dL≈1.441d_{L}\approx 1.441.Kuznetsov and Reitmann [18] derive a rigorous bound on the Lyapunov dimension of the global attractor in the Duffing oscillator which, for these parameter values, gives a bound dL≤1.9910d_{L}\leq 1.9910. Note that we have not defined LEs or Lyapunov dimension for nonautonomous systems. This is straightforward to do, but instead we convert 107 to an autonomous system.

The ODE 107 is not a polynomial system, as required for the formulations presented in section 2, but it can be transformed into a polynomial system by a standard trick (cf. Parker et al. [38]). We let x1=Ax_{1}=A and x2=A˙x_{2}=\dot{A}, and we introduce x3=sin⁡ω​tx_{3}=\sin{\omega t} and x4=cos⁡ω​tx_{4}=\cos{\omega t}. Then the x​(t)x(t) vector evolves according to

Recalling that x3x_{3} and x4x_{4} are sine and cosine of the same argument, we are interested only in dynamics on the manifold

| ℬ={x∈ℝ4|x32+x42=1}.\mathcal{B}=\{x\in\mathbb{R}^{4}\,|\,x_{3}^{2}+x_{4}^{2}=1\}. | (109) | | -------------------------------------------------------------------------------------------- | ----- |

Since only odd powers of xx appear in 108, the right-hand side is equivariant under the simple sign symmetry Λ=−I\Lambda=-I, and ℬ\mathcal{B} is invariant, so we can also impose that our polynomials VV are invariant under the same symmetry. This is exploited in our code using the methods described in section 2.5.

For the system 108, the Jacobian is

From this, we compute the additive compound matrices

In this case, the fourth additive compound matrix—the trace of the Jacobian—is independent of xx, and therefore either of our methods gives the sharp result that M4=−δM_{4}=-\delta on any orbit. This implies that the dimension dLd_{L} is optimized on the same orbit as the optimizer for M3M_{3} and, given a sharp bound for M3M_{3}, proposition 2.7 will give a sharp bound for dLd_{L}. Furthermore, since there are no equilibria, every trajectory must have a zero LE associated with the flow, and an additional vanishing LE corresponding to amplitude change in the x3x_{3} and x4x_{4} variables. Therefore M1=M2=M3M_{1}=M_{2}=M_{3} for all orbits of interest, and it suffices to compute bounds on M1M_{1} only in order to deduce bounds for M2M_{2}, M3M_{3} and dLd_{L}. Nevertheless, to test our methods we computed bounds using our formulations for M1M_{1}, M2M_{2}, M3M_{3} and dLd_{L}.

To formulate SOS optimization problems, we must choose finite dimensional spaces for the polynomials to be optimized over. In all cases we choose a maximum degree of VV and then construct a basis from symmetry-invariant monomials. For each SOS constraint, the restriction that hj​(x)=0h_{j}(x)=0, where hj=x32+x42−1h_{j}=x_{3}^{2}+x_{4}^{2}-1, is enforced via additional unknown polynomials ρj\rho_{j}, whose maximum degrees are chosen so that hj​ρjh_{j}\rho_{j} matches the maximum degree of other terms in the expression. See the accompanying code for full details.

For the sphere projection methods proposition 2.15 and proposition 2.16, computational results are reported in tables 1 and 2 respectively. In both cases, clearly converged bounds are achieved with VV of degree 4 in xx and 2 in zkz^{k}. For the shifted spectrum methods proposition 2.17 and proposition 2.18, numerical errors are large when the degree of VV in xx is 4 or smaller, so we let VV have degree 66 in xx. In this case, we find feasible bounds M1≤0.887898M_{1}\leq 0.887898 and dL≤3.8162d_{L}\leq 3.8162, and we find infeasible bounds when each final digit is decreased by 1. These bounds perfectly match those found via the sphere projection methods.

Table 1: Upper bounds on the sum MkM_{k} of the kk leading LEs, computed by solving the SOS optimization problem in proposition 2.15 for the Duffing oscillator in the variables of 108. Polynomial degrees are at most dxd_{x} in xx and dzd_{z} in zkz^{k}. Tabulated values are rounded up to the precision shown, and the underlines indicate the digits that agree (after rounding) with the values computed by numerical integration for the periodic orbit of figure 1. On this orbit M4=0.2M_{4}=0.2 and M1=M2=M3≈0.88785M_{1}=M_{2}=M_{3}\approx 0.88785.

Table 2: Upper bounds on the Lyapunov dimension dLd_{L}, computed by solving the SOS optimization problem in proposition 2.16 with k=3k=3 for the Duffing oscillator in the variables of 108. Polynomial degrees are at most dxd_{x} in xx and dzd_{z} in zkz^{k} and zk+1z^{k+1}. Tabulated values are rounded up to the precision shown, and underlines indicate digits that agree (after rounding) with the value dL≈3.81615d_{L}\approx 3.81615 computed by numerical integration for the periodic orbit of figure 1.

Figure 1: A chaotic trajectory (grey) in the system 108, and the unstable periodic orbit of period 2​π2\pi discussed in the text (blue).

To find an orbit that maximises the leading LE M1M_{1}, and consequently dLd_{L} also, we follow a strategy similar to that of Tobasco et al. [46]. We restrict to a Poincaré section of ℬ\mathcal{B} with x3=0x_{3}=0 and x4=1x_{4}=1, then we attempt to find the global minimum of the expression

in the constraint of 83, where BB and VV are the approximately optimal values returned by the SDP solver for the sphere projection method. This expression must vanish on average over the optimal orbit, and it is everywhere non-negative, so close to the orbit it must be very small in value. The minimization strategy is to plot, at each x1x_{1} and x2x_{2} over the approximate range of the chaotic attractor, the minimum of this expression over all zk∈𝕊nk−1z^{k}\in\mathbb{S}^{n_{k}-1}, which is approximated by random sampling on the unit sphere. This minimum value is then passed to a Newton-shooting algorithm to converge the true periodic orbit, which approximately intersects the point

and has a period of 2​π2\pi. This is the small, approximately circular periodic orbit depicted in figure 1. We computed LEs on this orbit to be approximately

giving

The Kaplan-Yorke formula 18 then gives dL=3.81615d_{L}=3.81615, given that the supremum must be attained on this orbit. These values perfectly match the bounds given by C1C_{1}, C2C_{2}, C3C_{3} and C4C_{4}. In the limit γ→0\gamma\to 0 of 107, the orbit becomes an equilibrium point at the origin. This matches previous observations that global Lyapunov dimension is attained on equilibria in systems that have such points. Finally, we note that our bounds on dLd_{L} given by C2C_{2} and C4C_{4} of 108 are equivalent to dL≤1.8162d_{L}\leq 1.8162 for 107, a significant improvement over the bound of Kuznetsov and Reitmann [18], and that this result must be sharp as an equivalent periodic orbit exists within the nonautonomous system 107.

3.2 A three-degree-of-freedom Hamiltonian system

For an example where LE sums for different kk are maximized on different orbits, we consider a higher-dimensional system whose large number of symmetries keeps the computational cost tractable. In particular, we consider a Hamiltonian system with Hamiltonian function

| H​(x)=12​|x|2+x12​x22+x22​x32+x32​x12,H(x)=\frac{1}{2}\left|x\right|^{2}+x_{1}^{2}x_{2}^{2}+x_{2}^{2}x_{3}^{2}+x_{3}^{2}x_{1}^{2}, | (118) | | --------------------------------------------------------------------------------------------------------------------------------------------- | ----- |

the first three components of xx being the position variables, and the final three the corresponding momenta. This is a specific form of a triaxially symmetric potential studied by Palacián et al. [35], which has applications in galactic dynamics [6, 4]. The Hamiltonian 118 gives rise to a system of ODEs in 6 variables,

The right-hand side of 119 is equivariant under x↦Λ​xx\mapsto\Lambda x where Λ\Lambda is any matrix in the subgroup of O​(6)O(6) generated by

We consider the fixed energy level H=1H=1. At this value of HH, there are no equilibria, so μ3=μ4=0\mu_{3}=\mu_{4}=0 on all orbits. We consider only orbits whose position variables remain in the region where 12​(x12+x22+x32)+(x12​x22+x22​x32+x32​x12)≤1\frac{1}{2}\left(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}\right)+\left(x_{1}^{2}x_{2}^{2}+x_{2}^{2}x_{3}^{2}+x_{3}^{2}x_{1}^{2}\right)\leq 1, depicted in figure 2, and whose full state vector xx remains in the compact sphere 12​|x|2≤1\frac{1}{2}|x|^{2}\leq 1. Two families of particularly short periodic orbits, which are also depicted in figure 2, were found by a brute-force search over random initial conditions. The first family, which we call S1S_{1}, has a period of approximately 5.278 and passes near

or near points that map to this xx under a symmetry. The S1S_{1} orbits have LEs of approximately

corresponding to LE sums of

The second family of orbits, which we call S2S_{2}, has a period of approximately 4.857 and passes near

or near a symmetry-related point. These S2S_{2} orbits have LEs of approximately

corresponding to sums of

Figure 2: Plots of the position variables (x1,x2,x3)(x_{1},x_{2},x_{3}) for the Hamiltonian system 118 at the energy level H=1H=1. Left: The family of periodic orbits S1S_{1}. Right: The family of periodic orbits S2S_{2}. Trajectories are confined to stay within the region 12​(x12+x22+x32)+(x12​x22+x22​x32+x32​x12)≤1\frac{1}{2}\left(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}\right)+\left(x_{1}^{2}x_{2}^{2}+x_{2}^{2}x_{3}^{2}+x_{3}^{2}x_{1}^{2}\right)\leq 1, which is bounded by the shaded surface. On S1S_{1}, the LE sums are M1=M2=M3=M4=M5=0.42913M_{1}=M_{2}=M_{3}=M_{4}=M_{5}=0.42913 and on S2S_{2}, M1=M5=0.25919M_{1}=M_{5}=0.25919 and M2=M3=M4=0.51837M_{2}=M_{3}=M_{4}=0.51837.

Table 3: Upper bounds on the sum MkM_{k} of the kk leading LEs, computed by solving the SOS optimization problem in proposition 2.15 for the Hamiltonian system 119. Polynomial degrees are at most dxd_{x} in xx and dzd_{z} in zkz^{k}. Tabulated values are rounded up to the precision shown, and the underlines indicate the digits that agree (after rounding) with the values computed by numerical integration for the periodic orbits S1S_{1} (k=1,5,6k=1,5,6) or S2S_{2} (k=2,3,4k=2,3,4) of figure 2. Values are missing in cases where computations were prohibitively memory-intensive.

Table 4: Upper bounds on the sum MkM_{k} of the kk leading LEs, computed by solving the SOS optimization problem in proposition 2.17 for the Hamiltonian system 119. Polynomial degrees are at most dxd_{x} in xx and 2 in wkw^{k}. The values in this table are not rounded. The underlines indicate the digits that agree (after rounding) with the values computed by numerical integration for the periodic orbits S1S_{1} (k=1,5,6k=1,5,6) or S2S_{2} (k=2,3,4k=2,3,4) of figure 2.

Because this system is Hamiltonian, volumes in state space are conserved, and our methods for bounding Lyapunov dimension are not relevant. We carried out SOS computations to bound LE sums using the sphere projection and shifted spectrum approaches, and tables 3 and 4 give the results. In both cases we add the explicit constraint |x|2≤2|x|^{2}\leq 2 along with H​(x)=1H(x)=1 in the definition of the set ℬ\mathcal{B} because this improves the numerical results. Using the shifted spectrum approach, we find sharp bounds for M1M_{1}, M2M_{2}, M4M_{4}, M5M_{5} and M6M_{6}. Using the sphere projection approach, we found sharp bounds on M1M_{1}, M5M_{5} and M6M_{6}. Sharp bounds for the other kk values seemed to require larger polynomial degrees that were beyond our computational resources.

Given the Hamiltonian structure of the system, bounding M5M_{5} is directly equivalent to bounding M1M_{1}, bounding M4M_{4} is equivalent to bounding M2M_{2}, and bounding M6M_{6} is trivial, since this is always zero. These facts are borne out by the results. Bounding the sum M3M_{3} is not a priori identical to any other case, but the lack of equilibria means that we expect the true value to be equal to M2M_{2} and M4M_{4}. Furthermore, bounds for M3M_{3} require a particularly high-dimensional augmented system of dimension n+nk=26n+n_{k}=26, which is computationally challenging at even modest polynomial degree despite exploiting all symmetries of the system. This is especially true for the sphere projection method, which has higher polynomial degree in the tangent space dynamics. In table 3, the missing values were not possible to compute when solving SDPs with Mosek due to insufficient memory, despite having 3TB of memory. This could be surmounted by using a first-order SDP solver at the expense of very long runtimes and/or lower precision. In table 4, the M3M_{3} result is not sharp at the maximum degree dx=6d_{x}=6 that was computationally feasible.

4 Conclusion

We have described two distinct approaches to bounding sums of Lyapunov exponents as well as Lyapunov dimension in ODE systems. The general versions of our formulations require optimizing over auxiliary functions in the class C1C^{1}. We prove that most of these formulations give sharp bounds under weak assumptions, but the optimization problems are not tractable in general. In the case of ODEs with polynomial right-hand sides, our bounding formulations can be relaxed using SOS constraints, which we again prove are convergent, and which are often computationally tractable using SOS optimization. Carrying out such SOS computations for particular ODEs, we obtained nearly perfect bounds in challenging examples. For one of the examples, the sphere projection method gave much better conditioned SDPs in practice. For the other example, the shifted spectrum method was more appropriate due to onerous memory requirements of the sphere projection method. In section 3.1, we demonstrated how the results of our bounding methods could be used to find the optimal, potentially previously undetected, periodic orbits which maximize the bounds.

All of our methods could be used to construct fully rigorous bounds via a computer assisted proof, making use of interval arithmetic or rational projection [12, 37]. In both of our examples, the results of the SOS programs allowed us to find periodic orbits that attained our bounds to within several digits, thus confirming their sharpness. In principle, computer-assisted proofs based on our framework could also give rigorous upper bounds on LEs in partial differential equations, by combining our methods with some in Goluskin and Fantuzzi [13], but such computations would be challenging.

In the thesis of Oeri [33], an analogous approach to the sphere projection method was formulated for bounding the leading LE of discrete-time dynamical systems, with an SOS implementation for polynomial maps. This can be directly adapted to sums of LEs by using additive (rather than multiplicative) compound matrices. An approach analogous to our shifted spectrum method is also possible for bounding sums of LEs for discrete-time systems, and methods can additionally be derived for the Lyapunov dimension. In general, these are computationally difficult because the equivalent of the Lie derivative is a composition of the dynamics with an auxiliary function, leading to very high polynomial degrees. This may be the focus of a future publication.

Data availability statement

Acknowledgements

The authors thank Samuel Punshon-Smith and Anthony Quas for useful discussions. During this work DG was supported by the NSERC Discovery Grants Program (awards RGPIN-2018-04263 and RGPIN-2025-06823). Computational resources were provided by the Digital Research Alliance of Canada.

Appendix A Exterior products and compound matrices

The exterior (or wedge or Grassmann) product is a multilinear map of vectors, defined by the property that x2∧x1=−x1∧x2x_{2}\wedge x_{1}=-x_{1}\wedge x_{2}, and x1∧x2=0x_{1}\wedge x_{2}=0 only if x1=α​x2x_{1}=\alpha x_{2} for some α\alpha, or x2=0x_{2}=0. An exterior product of kk vectors in a vector space of dimension n≥kn\geq k is called a kk-vector. The properties of the exterior product imply that permuting the xix_{i} in x1∧⋯∧xkx_{1}\wedge\dots\wedge x_{k} either leaves the kk-vector unchanged or negates it, and that the kk-vector is nonzero if and only if the xix_{i} are linearly independent. The set of kk-vectors in ℝn\mathbb{R}^{n} is therefore a linear space of dimension nk=n!(n−k)!n_{k}=\tfrac{n!}{(n-k)!}. We identify this space with ℝnk\mathbb{R}^{n_{k}} and choose as a basis the set of kk-vectors of the form

where i1<…<iki_{1}<\ldots<i_{k} and eie_{i} is the it​hi^{th} standard basis vector in ℝnk\mathbb{R}^{n_{k}}. The standard Euclidean norm on this space satisfies the property that

| |x1∧⋯∧xk||x_{1}\wedge\dots\wedge x_{k}| | (121) | | --------------------------------------------- | ----- |

is the volume of the parallepiped spanned by x1,…,xk∈ℝnx_{1},\dots,x_{k}\in\mathbb{R}^{n} [47].

Linear maps on vectors in ℝn\mathbb{R}^{n} induce linear maps on corresponding kk-vectors in ℝnk\mathbb{R}^{n_{k}}. Likewise, linear ODEs on ℝn\mathbb{R}^{n} induce linear ODEs on kk-vectors. As defined below, the map on kk-vectors induced by a map AA is captured by the multiplicative compound A(k)A^{(k)}, and the ODE for kk-vectors induced by an ODE with right-hand side AA is captured by the additive compound A[k]A^{[k]}. The multiplicative compound is common in many different areas of mathematics, and is known variously as the compound representation, exterior power representation, wedge power representation or simply the induced action of a linear map. Compound matrices have been applied previously in the calculation and bounding of LEs [26, 29, 30, 18] and also stability of dynamical systems [31]. Here we state key definitions and theorems; a more detailed introduction can be found in Marshall et al. [28, Chapter 19F], in more generality for complex non-square matrices.

Definition A.1.

For any real matrix A∈ℝn×nA\in\mathbb{R}^{n\times n} and integer k∈{1,…,n}k\in\{1,\ldots,n\}:

1. 1. The kt​hk^{th} multiplicative compound of AA is the real matrix A(k)∈ℝnk×nkA^{(k)}\in\mathbb{R}^{{n_{k}}\times n_{k}} whose entries are the k×kk\times k subdeterminants of AA, taken in lexicographic order.
2. 1. The kt​hk^{th} additive compound of AA is the real matrix A[k]∈ℝnk×nkA^{[k]}\in\mathbb{R}^{{n_{k}}\times n_{k}} such that
      | A[k]=dd​δ​(I+δ​A)(k)|δ=0.A^{[k]}=\frac{\mathrm{d}}{\mathrm{d}\delta}\left.\left(I+\delta A\right)^{(k)}\right|_{\delta=0}. | (122) |
      | -------------------------------------------------------------------------------------------------------------------------------------------- | ----- |

The lexicographic order in A.1 refers to the ordering of all nkn_{k} permutations of kk increasing integers chosen from {1,…,n}\{1,\ldots,n\}. As an example, the lexicographic order of k=3k=3 integers in the case n=4n=4 is

Therefore, the (2,3)(2,3) entry of A(3)A^{(3)} is

where ap,qa_{p,q} denotes the (p,q)(p,q) element of AA. Note that the row and column indices come from the second and third permutations, respectively. It follows from the definitions that A(1)=A[1]=AA^{(1)}=A^{[1]}=A, and A(n)=detAA^{(n)}=\det{A} whereas A[n]=tr​AA^{[n]}=\mathrm{tr}A. In fact, the eigenvalues of multiplicative and additive compounds are respectively the products and sums of the eigenvalues of the original matrix, which gives rise to their names. An explicit example of additive compound matrices is given in 111.

The next two propositions list additional properties that will be useful to us for multiplicative and additive compounds, respectively. In proposition A.2, part 1 is a special case of the Cauchy–Binet theorem. Parts 2–4 follow quickly from the definition, and the last part is theorem 6.2.10 of Fiedler [10].

Proposition A.2 (Properties of multiplicative compound matrices).

For any real square matrices AA and BB:

1. 1. (A​B)(k)=A(k)​B(k)(AB)^{(k)}=A^{(k)}B^{(k)}.
2. 1. (A(k))𝖳=(A𝖳)(k)\left(A^{(k)}\right)^{\mathsf{T}}=\left(A^{\mathsf{T}}\right)^{(k)}.
3. 1. If AA is invertible, then (A−1)(k)=(A(k))−1\left(A^{-1}\right)^{(k)}=\left(A^{(k)}\right)^{-1}.
4. 1. AA is orthogonal if and only if A(k)A^{(k)} is orthogonal.
5. 1. For kk-vectors identified with ℝnk\mathbb{R}^{n_{k}} and for A∈ℝn×nA\in\mathbb{R}^{n\times n},

      A(k)​(x1∧⋯∧xk)=(A​x1)∧⋯∧(A​xk).A^{(k)}\left(x_{1}\wedge\dots\wedge x_{k}\right)=(Ax_{1})\wedge\dots\wedge(Ax_{k}).	(125)

Proposition A.3 (Properties of additive compound matrices).

For any real square matrices:

1. 1. (A+B)[k]=A[k]+B[k]\left(A+B\right)^{[k]}=A^{[k]}+B^{[k]}.
2. 1. A[k]​(x1∧⋯∧xk)=(A​x1∧x2∧⋯∧xk)+(x1∧A​x2∧⋯∧xk)+⋯+(x1∧x2∧⋯∧A​xk).A^{[k]}\left(x_{1}\wedge\dots\wedge x_{k}\right)=(Ax_{1}\wedge x_{2}\wedge\dots\wedge x_{k})+(x_{1}\wedge Ax_{2}\wedge\dots\wedge x_{k})+\dots+(x_{1}\wedge x_{2}\wedge\dots\wedge Ax_{k}).
3. 1. If QQ is invertible, (Q​A​Q−1)[k]=Q(k)​A[k]​(Q(k))−1\left(QAQ^{-1}\right)^{[k]}=Q^{(k)}A^{[k]}\left(Q^{(k)}\right)^{-1}.

Proof.

Part 1 is immediate from the definition. For part 2, see theorem 2f of [24]. Part 3 follows by calculation,

| (Q​A​Q−1)[k]\displaystyle\left(QAQ^{-1}\right)^{[k]} | =dd​δ​(I+δ​Q​A​Q−1)(k)|δ=0\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\delta}\left.\left(I+\delta QAQ^{-1}\right)^{(k)}\right|_{\delta=0} | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------ | | =dd​δ​(Q​(I+δ​A)​Q−1)(k)|δ=0\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\delta}\left.\left(Q\left(I+\delta A\right)Q^{-1}\right)^{(k)}\right|_{\delta=0} | | | =dd​δ​Q(k)​(I+δ​A)(k)​(Q(k))−1|δ=0\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\delta}\left.Q^{(k)}\left(I+\delta A\right)^{(k)}\left(Q^{(k)}\right)^{-1}\right|_{\delta=0} | |

where the last equality uses properties 1 and 2 from proposition A.2. ∎

Multiplicative compound matrices are straightforward to compute directly, at least when nn is small. Additive compound matrices can be computed with symbolic manipulation. In our computational examples we compute additive compounds using Julia with the DynamicPolynomials.jl package [21]. Explicit formulas for the elements of additive compounds also exist [27].

We now have the ingredients to prove 15:

Lemma A.4.

If

for i=1,…,ki=1,\dots,k, and we define

then

Proof.

We have

so, using proposition A.2,

Observe that ϕ0​(x0)=x0\phi^{0}(x_{0})=x_{0}, and ϕt+δ​(x0)=ϕδ​(ϕt​(x0))\phi^{t+\delta}(x_{0})=\phi^{\delta}(\phi^{t}(x_{0})), so that

Then by proposition A.2,

so by the usual definition of the derivative,

Furthermore, we have

and so

by A.1. ∎

Appendix B A specialised Lyapunov theorem

Here we present a Lyapunov theorem for the stability of one variable w∈ℝmw\in\mathbb{R}^{m} whose linear dynamics is controlled by a second variable xx. This is applicable to both proposition 2.11 and proposition 2.12. The second part, whose proof follows that of similar results in Khalil [17], is the so-called converse Lyapunov theorem, which states that under stronger conditions, it is always possible to find a Lyapunov function and that this VV is quadratic in ww.

Lemma B.1.

For a domain ℬ⊂ℝn\mathcal{B}\subset\mathbb{R}^{n}, let x​(t)x(t) and w​(t)w(t) solve the ODE system

on ℬ×ℝm\mathcal{B}\times\mathbb{R}^{m}, with f∈C0​(ℬ,ℝn)f\in C^{0}(\mathcal{B},\mathbb{R}^{n}) and A∈C0​(ℬ,ℝm×m)A\in C^{0}(\mathcal{B},\mathbb{R}^{m\times m}).

1. 1. If there exists V∈C1​(ℬ×ℝm)V\in C^{1}\left(\mathcal{B}\times\mathbb{R}^{m}\right) such that
      | V​(x,w)≥|w|2 and f​(x)⋅Dx​V​(x,w)+[A​(x)​w]⋅Dw​V​(x,w)≤0V(x,w)\geq|w|^{2}\quad\text{ and }\quad f(x)\cdot D_{x}V(x,w)+\left[A(x)w\right]\cdot D_{w}V(x,w)\leq 0 | (134) |
      | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ----- |
      for all (x,w)∈ℬ×ℝm(x,w)\in\mathcal{B}\times\mathbb{R}^{m}, then w​(t)w(t) is bounded forward in time for each w​(0)∈ℝmw(0)\in\mathbb{R}^{m}.
2. 1. If A∈C1​(ℬ,ℝm×m)A\in C^{1}(\mathcal{B},\mathbb{R}^{m\times m}) is bounded and f∈C1​(ℬ,ℝn)f\in C^{1}(\mathcal{B},\mathbb{R}^{n}), and if there exist K>0K>0 and λ>0\lambda>0 such that
      | |w​(t)|≤K​e−λ​t​|w​(0)||w(t)|\leq Ke^{-\lambda t}|w(0)| | (135) |
      | ---------------------------------------------------------- | ----- |
      for all solutions to 133, then there exists U∈C1​(ℬ,ℝm×m)U\in C^{1}\left(\mathcal{B},\mathbb{R}^{m\times m}\right) such that 134 is satisfied with V​(x,w)=w𝖳​U​(x)​wV(x,w)=w^{\mathsf{T}}U(x)w.

Proof.

For the first part, let V∈C1V\in C^{1} satisfy the two inequality conditions. For contradiction, suppose that there is a solution with |w​(t)||w(t)| unbounded. The first VV condition then implies that V​(x​(t),w​(t))V(x(t),w(t)) is unbounded forward in time. The second VV condition implies

This expression is continuous in tt, so we can integrate to show that VV remains bounded for all time, which is a contradiction.

For the second part, fix a time T>0T>0, to be chosen later. Define

| V​(x0,w0)=2​L1−e−2​L​T​∫0T|w​(t)|2​𝑑t,V(x_{0},w_{0})=\frac{2L}{1-e^{-2LT}}\int_{0}^{T}|w(t)|^{2}\,dt, | (137) | | -------------------------------------------------------------------------------------------------------------- | ----- |

where (x​(t),w​(t))(x(t),w(t)) evolves under (133) from initial condition (x​(0),w​(0))=(x0,w0)(x(0),w(0))=(x_{0},w_{0}), and LL is an upper bound for the operator norm ‖A​(x)‖\|A(x)\| over x∈ℬx\in\mathcal{B}. Observe that, by linearity of w​(t)w(t) with respect to w0w_{0}, this is quadratic in w0w_{0} and thus can be written

for some symmetric matrix U​(x0)U(x_{0}). Since A∈C1​(ℬ,ℝm×m)A\in C^{1}(\mathcal{B},\mathbb{R}^{m\times m}), it follows that U∈C1​(ℬ,ℝm×m)U\in C^{1}(\mathcal{B},\mathbb{R}^{m\times m}).

Note that

| |dd​t​|w|2|\displaystyle\left|\frac{\mathrm{d}}{\mathrm{d}t}|w|^{2}\right| | ≤2​|w|​|A​(x)​w|\displaystyle\leq 2|w|\left|A(x)w\right| | (139) | | --------------------------------------------------------------------------------- | ------------------------------------------------------------- | ----- | | ≤2​‖A​(x)‖​|w|2\displaystyle\leq 2\|A(x)\||w|^{2} | | | | ≤2​L​|w|2\displaystyle\leq 2L|w|^{2} | | |

i.e. dd​t​|w|2≥−2​L​|w|2\frac{\mathrm{d}}{\mathrm{d}t}|w|^{2}\geq-2L|w|^{2}. Thus Grönwall’s inequality gives us

| |w​(t)|2≥e−2​L​t​|w0|2|w(t)|^{2}\geq e^{-2Lt}|w_{0}|^{2} | (140) | | ----------------------------------------------------------- | ----- |

and so

| V​(x0,w0)≥2​L1−e−2​L​T​∫0Te−2​L​t​|w0|2​𝑑t=|w0|2,V(x_{0},w_{0})\geq\frac{2L}{1-e^{-2LT}}\int_{0}^{T}e^{-2Lt}|w_{0}|^{2}dt=|w_{0}|^{2}, | (141) | | ------------------------------------------------------------------------------------------------------------------------------------------------- | ----- |

satisfying the first condition of 134.

For the second condition of 134, we define

| v​(t)=V​(x​(t),w​(t))=2​L1−e−2​L​T​∫tt+T|w​(s)|2​ds.v(t)=V(x(t),w(t))=\frac{2L}{1-e^{-2LT}}\int_{t}^{t+T}|w(s)|^{2}\,{\rm d}s. | (142) | | -------------------------------------------------------------------------------------------------------------------------------------- | ----- |

Then,

| f​(x​(t))⋅Dx​V​(x​(t),w​(t))+[A​(x​(t))​w​(t)]⋅Dw​V​(x​(t),w​(t))=d​vd​t=2​L1−e−2​L​T​(|w​(t+T)|2−|w​(t)|2)≤2​L1−e−2​L​T​(K2​e−2​λ​T−1)​|w​(t)|2,\begin{split}f(x(t))&\cdot D_{x}V(x(t),w(t))+\left[A(x(t))w(t)\right]\cdot D_{w}V(x(t),w(t))\\ &=\frac{\mathrm{d}v}{\mathrm{d}t}\\ &=\frac{2L}{1-e^{-2LT}}\left(|w(t+T)|^{2}-|w(t)|^{2}\right)\\ &\leq\frac{2L}{1-e^{-2LT}}\left(K^{2}e^{-2\lambda T}-1\right)|w(t)|^{2},\end{split} | (143) | | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ----- |

so if T≥1λ​log⁡KT\geq\frac{1}{\lambda}\log{K} the condition is satisfied. ∎

References

• ApS [2025] M. ApS. MOSEK Optimizer API for Julia, 2025.
• Arnold [1998] L. Arnold. Random dynamical systems. Springer Monographs in Mathematics. Springer-Verlag, 1998.
• Bochi [2018] J. Bochi. Ergodic optimization of Birkhoff averages and Lyapunov exponents. In Proc. Int. Congr. Math., pages 1825–1846, 2018.
• Caranicolas [1994] N. Caranicolas. 1:1:1 resonant periodic orbits in 3-dimensional galactic-type Hamiltonians. Astronomy and Astrophysics, 282:34–36, 1994.
• Coey et al. [2022] C. Coey, L. Kapelevich, and J. P. Vielma. Solving natural conic formulations with Hypatia.jl. INFORMS Journal on Computing, 34(5):2686–2699, 2022.
• de Zeeuw [1985] T. de Zeeuw. Motion in the core of a triaxial potential. Monthly Notices of the Royal Astronomical Society, 215(4):731–760, 1985.
• Doering and Gibbon [1995] C. R. Doering and J. Gibbon. On the shape and dimension of the Lorenz attractor. Dynamics and Stability of Systems, 10(3):255–268, 1995.
• Eden et al. [1991] A. Eden, C. Foias, and R. Temam. Local and global Lyapunov exponents. Journal of Dynamics and Differential Equations, 3:133–177, 1991.
• Eden [1989] A. O. Eden. An abstract theory of L-exponents with applications to dimension analysis. Indiana University, 1989.
• Fiedler [2008] M. Fiedler. Special matrices and their applications in numerical mathematics. Courier Corporation, 2008.
• Gatermann and Parrilo [2004] K. Gatermann and P. A. Parrilo. Symmetry groups, semidefinite programs, and sums of squares. Journal of Pure and Applied Algebra, 192(1-3):95–128, 2004.
• Goluskin [2018] D. Goluskin. Bounding averages rigorously using semidefinite programming: mean moments of the Lorenz system. Journal of Nonlinear Science, 28:621–651, 2018.
• Goluskin and Fantuzzi [2019] D. Goluskin and G. Fantuzzi. Bounds on mean energy in the Kuramoto–Sivashinsky equation computed using semidefinite programming. Nonlinearity, 32(5):1705, 2019.
• Goulart and Chen [2024] P. J. Goulart and Y. Chen. Clarabel: An interior-point solver for conic programs with quadratic objectives, 2024.
• Kaluba et al. [2021] M. Kaluba, D. Kielak, and P. W. Nowak. On property (T)(T) for A​u​t​(Fn)Aut(F_{n}) and S​Ln​(Z)SL_{n}(Z). Annals of Mathematics, 193(2):539–562, 2021.
• Kaplan and Yorke [1979] J. Kaplan and J. Yorke. Chaotic behavior of multidimensional difference equations. Functional Differential Equations and Approximation of Fixed Points, pages 204–227, 1979.
• Khalil [2002] H. Khalil. Nonlinear systems. Prentice Hall, 2002.
• Kuznetsov and Reitmann [2020] N. Kuznetsov and V. Reitmann. Attractor dimension estimates for dynamical systems: theory and computation. Springer, 2020.
• Lakshmi et al. [2020] M. Lakshmi, G. Fantuzzi, J. Fernández-Caballero, Y. Hwang, and S. Chernyshenko. Finding extremal periodic orbits with polynomial optimisation, with application to a nine-mode model of shear flow. SIAM J. Appl. Dyn. Syst., 19:763–787, 2020.
• Legat et al. [2017] B. Legat, C. Coey, R. Deits, J. Huchette, and A. Perry. Sum-of-squares optimization in Julia. In The First Annual JuMP-dev Workshop, 2017.
• Legat et al. [2024] B. Legat, S. Timme, T. Weisser, L. Kapelevich, N. Sajko, Manuel, A. Demin, C. Hansknecht, C. Rackauckas, J. TagBot, M. Forets, and T. Holy. JuliaAlgebra/DynamicPolynomials.jl: v0.5.7, May 2024.
• Leonov and Lyashko [1993] G. A. Leonov and S. Lyashko. Eden’s hypothesis for a Lorenz system. Vestnik St. Petersburg University: Mathematics, 26(3):15–18, 1993.
• Leonov et al. [2016] G. A. Leonov, N. V. Kuznetsov, N. Korzhemanova, and D. Kusakin. Lyapunov dimension formula for the global attractor of the Lorenz system. Communications in Nonlinear Science and Numerical Simulation, 41:84–103, 2016.
• London [1976] D. London. On derivations arising in differential equations. Linear and Multilinear Algebra, 4(3):179–189, 1976.
• Lorentz [1986] G. G. Lorentz. Approximation of Functions. Chelsea Publishing Company, 2nd edition, 1986.
• Manika [2013] D. Manika. Application of the compound matrix theory for the computation of Lyapunov exponents of autonomous Hamiltonian systems. Master’s thesis, Aristotle University of Thessaloniki, 2013.
• Margaliot and Sontag [2019] M. Margaliot and E. D. Sontag. Revisiting totally positive differential systems: A tutorial and new results. Automatica, 101:1–14, 2019.
• Marshall et al. [2010] A. Marshall, I. Olkin, and B. Arnold. Inequalities: Theory of Majorization and Its Applications. Springer Series in Statistics. Springer New York, 2010.
• Martini et al. [2022] D. Martini, D. Angeli, G. Innocenti, and A. Tesi. Ruling out positive Lyapunov exponents by using the Jacobian’s second additive compound matrix. IEEE Control Systems Letters, 6:2924–2928, 2022.
• Martini et al. [2023] D. Martini, D. Angeli, G. Innocenti, and A. Tesi. Bounding Lyapunov exponents through second additive compound matrices: Case studies and application to systems with first integral. International Journal of Bifurcation and Chaos, 33(10):2350114, 2023.
• Muldowney [1990] J. S. Muldowney. Compound matrices and ordinary differential equations. The Rocky Mountain Journal of Mathematics, pages 857–872, 1990.
• Murty and Kabadi [1987] K. G. Murty and S. N. Kabadi. Some NP-complete problems in quadratic and nonlinear programming. Mathematical Programming: Series A and B, 39(2):117–129, 1987.
• Oeri [2023] H. Oeri. Convex optimization methods for bounding Lyapunov exponents. PhD thesis, University of Victoria, 2023.
• Oeri and Goluskin [2023] H. Oeri and D. Goluskin. Convex computation of maximal Lyapunov exponents. Nonlinearity, 36(10):5378, 2023.
• Palacián et al. [2017] J. F. Palacián, C. Vidal, J. Vidarte, and P. Yanguas. Periodic solutions and KAM tori in a triaxial potential. SIAM Journal on Applied Dynamical Systems, 16(1):159–187, 2017.
• Papachristodoulou and Prajna [2002] A. Papachristodoulou and S. Prajna. On the construction of Lyapunov functions using the sum of squares decomposition. In Proceedings of the 41st IEEE Conference on Decision and Control, 2002., volume 3, pages 3482–3487. IEEE, 2002.
• Parker [2024] J. P. Parker. The Lorenz system as a gradient-like system. Nonlinearity, 37(9):095022, 2024.
• Parker et al. [2021] J. P. Parker, D. Goluskin, and G. Vasil. A study of the double pendulum using polynomial optimization. Chaos, 31(10):103102, 2021.
• Parrilo [2000] P. Parrilo. Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization. PhD thesis, California Institute of Technology, 2000.
• Pikovsky and Politi [2016] A. Pikovsky and A. Politi. Lyapunov exponents: a tool to explore complex dynamics. Cambridge University Press, 2016.
• Putinar [1993] M. Putinar. Positive polynomials on semi-algebraic sets. Indiana University Mathematics Journal, 452:969–984, 1993.
• Rackauckas and Nie [2017] C. Rackauckas and Q. Nie. Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia. Journal of Open Research Software, 5(1):15, 2017.
• Reznick [2000] B. Reznick. Some concrete aspects of hilbert’s 17th problem. Contemporary Mathematics, 253(251-272), 2000.
• Scherer and Hol [2006] C. W. Scherer and C. W. J. Hol. Matrix sum of squares relaxations for robust semidefinite programs. Mathematical Programming, 107:189–211, 2006.
• Temam [2012] R. Temam. Infinite-dimensional dynamical systems in mechanics and physics, volume 68. Springer Science & Business Media, 2012.
• Tobasco et al. [2018] I. Tobasco, D. Goluskin, and C. R. Doering. Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems. Physics Letters A, 382(6):382–386, 2018.
• Winitzki [2009] S. Winitzki. Linear algebra via exterior products. Sergei Winitzki, 2009.