We extend the theory of spectral submanifolds (SSMs) to general non-autonomous dynamical systems that are either weakly forced or slowly varying. Examples of such systems arise in structural dynamics, fluid–structure interactions, and control problems. The time-dependent SSMs we construct under these assumptions are normally hyperbolic and hence will persist for larger forcing and faster time dependence that are beyond the reach of our precise existence theory. For this reason, we also derive formal asymptotic expansions that, under explicitly verifiable nonresonance conditions, approximate SSMs and their aperiodic anchor trajectories accurately for stronger, faster, or even temporally discontinuous forcing. Reducing the dynamical system to these persisting SSMs provides a mathematically justified model- reduction technique for non-autonomous physical systems whose time dependence is moderate either in magnitude or speed. We illustrate the existence, persistence, and computation of temporally aperiodic SSMs in mechanical examples under chaotic forcing.

Reduced models for complex physical systems are of growing interest in various areas of applied science and engineering. Mathematically justifiable reduction approaches, such as spectral submanifold (or SSM) reduction, seek to identify the internal dynamics on lower-dimensional, attracting invariant sets in the phase space of the system. These reduced dynamics then become viable reduced models for general trajectories that approach the invariant set and synchronize with its inner motions. SSM reduction also allows for periodic or quasiperiodic time dependence in the full system but has been inapplicable to systems with more general time dependence, such as impulsive, chaotic, or discontinuous forcing. This has hindered applications of SSM reduction to a number of problems in structural dynamics. Here, we remove this limitation by extending SSM theory of temporally aperiodic dynamical systems. We obtain exact results for cases of smooth small or smooth slow forcing, but find that our formulas for SSM-reduced dynamics extend to larger and faster forcing in physical examples, including even discontinuous chaotic forcing.

In its simplest form, a spectral submanifold (SSM) of an autonomous dynamical system is an invariant manifold W ( E ) that is tangent to a spectral subspace E of the linearized system at a fixed point (Haller and Ponsioen1). Classic examples of spectral submanifolds are the stable, unstable, and center manifolds tangent to spectral subspaces in which the linearized spectrum has eigenvalues with purely negative, positive, and zero real parts, respectively. The stable and unstable manifolds are well known to be unique and as smooth as the dynamical system, while center manifolds are non-unique and not all of them are guaranteed to be as smooth as the dynamical system (Guckenheimer and Holmes2 and Hirsch et al.3). Near fixed points that only have eigenvalues with negative and zero real parts, center manifolds attract all trajectories, and hence the dynamics restricted to it provides a mathematically exact reduced-order model for the full system (Carr4 and Roberts5).

Dissipative physical systems, however, generically have hyperbolic equilibria and hence admit no center manifolds. Instead, near their stable equilibria, such systems tend to have a set of fastest-decaying modes that die out quickly, leaving a set of slower decaying modes to govern the longer-term dynamics. Slow manifolds (i.e., SSMs constructed over such slower decaying modes) then replace center manifolds as targets for mathematically justified model reduction. Such slow SSMs were first targeted via Taylor expansions as nonlinear normal modes (NNMs) by Shaw and Pierre.6 These insightful calculations were then extended by the same authors to periodically and quasi-periodically forced mechanical systems to approximate forced mechanical response in a number of settings (see, e.g., the reviews by Kerschen et al.,7 Mikhlin and Avramov,8 Touzé et al.,9 and Mikhlin and Avramov10).

Later mathematical analysis of general SSMs yielded precise existence, uniqueness, and smoothness results for these manifolds. Specifically, if the spectral subspace E comprises either only decaying modes or only growing modes (i.e., E is a like-mode spectral subspace) with no integer resonance relationships to the modes outside of E, then the slow SSM family W ( E ) has a unique, primary member, W ( E ), that is as smooth as the full dynamical system (Cabré et al.11 and Haller and Ponsioen1). The remaining secondary (or fractional) members of the SSM family have reduced but precisely known order of differentiability that depends on the ratio of linearized decay rates outside E to those inside E (Haller et al.12). If we further disallow any integer resonance in the full linearized spectrum, then primary and fractional SSMs also exist when E is of mixed-mode type, i.e., spanned by a combination of stable and unstable linear modes (Haller et al.12).

All these results also hold for discrete autonomous dynamical systems, and hence SSM results also cover time-periodic continuous dynamical systems when applied to their Poincaré maps. Based on this fact, numerical implementations of time-periodic SSM calculations for periodically forced finite-element structures have appeared in Ponsioen et al.,13–15 Jain and Haller,16 and Vizzaccaro.17 An open source, equation-driven MATLAB toolbox (SSMTool) with a growing collection of worked problems is also available for mechanical systems with general nonlinearities (Jain et al.18). Data-driven construction of time-periodic SSMs have been developed and applied to numerical and experimental data by Cenedese et al.,19,20 Kaszás et al.21 and Axås et al.,22 with open-source MATLAB implementations (SSMLearn and fastSSM) available from Cenedese et al.23 

Outside autonomous and time-periodic non-autonomous systems, the existence of like-mode SSMs has only been treated under small, time-quasiperiodic perturbations of autonomous dynamical systems (Haro and de la Llave,24 Haller and Ponsioen,1 Opreni et al.,25 and Thurner et al.,26). In such systems, the role of hyperbolic fixed points as anchor points for SSMs is taken over by invariant tori whose dimension is equal to the number of rationally independent frequencies present in the forcing. An SSM in such a case perturbs from the direct product of an unperturbed invariant torus with an underlying spectral subspace E. Such SSMs have been shown to exist for like-mode spectral subspaces E under small quasi-periodic perturbations, provided that the real part of the linearized spectrum within E has no integer relationships with the real part of the spectrum outside E (Haro and de la Llave,24 Haro et al.,27 and Haller and Ponsioen1).

Related work by Fontich et al.28 covers the persistence and smoothness of primary SSMs emanating from an arbitrary attracting orbit of an autonomous dynamical system. While all non-autonomous systems become autonomous when their phase space is extended with the time variable, this theory does not apply under such an extension. The reason is that all attracting orbits lose hyperbolicity in the extended phase space due to the presence of the neutrally stable time direction.

In summary, while available SSM results have proven highly effective in equation-driven and data-driven reduced-order modeling of autonomous, time-periodic, and time-quasiperiodic systems, they offer no theoretical basis or computational scheme for physical systems with general aperiodic time-dependence. Yet the aperiodic setting is clearly of importance in a number of problems, including turbulent fluid–structure interactions, civil engineering structures subject to benchmark aperiodic forcing (such as those mimicking earthquakes) and control of robot motion.

In this paper, we extend available SSMs results to systems with general time dependence. While several powerful linearization results imply the existence of invariant manifolds for such non-autonomous dynamical systems, these results guarantee either no smoothness for SSM-type invariant manifolds (see, e.g., Palmer29) or rely on conditions involving the Lyapunov spectra, Sacker–Sell spectra, or dichotomy spectra of some associated non-autonomous linear systems of ODEs (see, e.g., Yomdin,30 Pötzsche,31 and Cuong et al.32). The latter types of conditions are intuitively clear but not readily verifiable, especially not in a data-driven setting. Here, our objective is to conclude the existence of time-aperiodic SSMs under directly computable conditions that also lead to explicitly computable SSM-reduced models in equation-driven and data-driven applications.

To this end, we consider two settings that arise frequently in practice: weak and additive non-autonomous time dependence and slowly varying (or adiabatic) time dependence. The first setting of weak non-autonomous external forcing is common in structural vibrations, wherein a structure’s steady-state response is of interest under various moderate loading scenarios. So far, related studies have been restricted to temporally periodic or quasiperiodic forcing (see, e.g., Ponsioen et al.,13–15 Li et al.,33,34 Jain and Haller,16 Vizzaccaro et al.,35 and Opreni et al.,25), because the existence and exact form of a steady state and its associated SSM have been unknown for more general forcing profiles. The second setting of slow time dependence arises, for instance, in control applications wherein the intended motion of a structure is generally much slower than the characteristic time scales of its internal vibrations. In those applications, the lack of an adiabatic SSM theory has so far confined model-reduction studies to small-amplitude trajectories along which a single, autonomous SSM computed at a nearby fixed point was used for modeling purposes (Alora et al.36,37).

In both of these non-autonomous settings, we use, modify, or extend prior invariant manifold results and techniques to conclude the existence of weakly aperiodic or adiabatic SSMs in the limit of small enough or slow enough time dependence, respectively. We then derive explicit recursive formulas for the arising non-autonomous SSMs and the aperiodic anchor trajectories to which they are attached. These formulas also cover and extend temporally periodic and aperiodic SSM computations to arbitrarily high order of accuracy. Using simple mechanical examples subjected to chaotic excitation, we illustrate that the new asymptotic formulas yield accurate reduced-order models even for larger and faster forcing.

Consider a non-autonomous dynamical system of the form
(1)
for some integer r 0 and with
(2)
on a compact neighborhood U R n. We consider system (1), a perturbation of the autonomous system,
(3)
that has a fixed point at x = 0 by our assumptions on f 0 ( x ). Both A and f 0 ( x ) may additionally depend on parameters, which we will suppress here for notational simplicity but will point out in the statement of our main results. Note that the time-dependent term f 1 ( x , t ) is also allowed to depend on the phase space variable x. Consequently, in addition to describing purely time-dependent external forcing, f 1 ( x , t ) can also capture what is commonly called parametric forcing in the structural vibration literature, i.e., time-dependence in the internal structure of the system.
We assume that the origin is a hyperbolic fixed point of (3), i.e., all eigenvalues in the spectrum of A,
(4)
listed in the order
(5)
satisfy
(6)
For simplicity of exposition, we assume that A is semisimple and hence has n eigenvectors e 1 , , e n C n corresponding to the eigenvalues listed in (4). The k t h eigenspace E k of A is then the linear span of the real and imaginary parts of the eigenvectors corresponding to the eigenvalue λ k, i.e.,
Note that any such eigenspace in an invariant subspace under the dynamics of the linearized ODE,
(7)
A spectral subspace E is the direct sum of a selected group of eigenspaces, i.e.,
(8)
Two important spectral subspaces are the stable subspace E s and the unstable subspace E u, defined as
(9)
At least one of E s and E u is nonempty due to the hyperbolicity assumption (6). As a consequence of this assumption, there exist also constants K , κ > 0 such that
(10)
This property of A is usually referred to as exponential dichotomy. For κ, we can select any positive number satisfying
In contrasts, the choice of K depends on the eigenvector geometry of A. Specifically, if A is a normal operator and hence has an orthogonal eigenbasis, we can select K = 1.

We first state general results on the fate of the x = 0 fixed point under the non-autonomous forcing term f 1 ( x , t ). Specifically, we give conditions under which this fixed point perturbs for moderately large | f 1 ( x , t ) | into a unique nearby hyperbolic trajectory x ( t ) of (1). This distinguished trajectory remains uniformly bounded for all times and has the same stability type as the x = 0 fixed point of system (3).

(Existence of anchor trajectory for non-autonomous SSMs)

Theorem 1
(Existence of anchor trajectory for non-autonomous SSMs)
Assume that in a ball B δ U of radius δ > 0 around x = 0, the functions f 0 and f 1 are of class C 0 and admit Lipschitz constants L 0 ( δ ) and L 1 ( δ ), respectively, in x for all t R. Assume further that the conditions
(11)
are satisfied with constants K , κ > 0 satisfying (10). Then, the following hold

  1. System (1) has a unique, uniformly bounded trajectory x ( t ) that remains is B δ for all t R and has the same stability type as the x = 0 fixed point of system (3).

  2. The trajectory x ( t ) is as smooth in any parameter as system (1).

Proof.

See Appendix A 1.

By statement (i) of Theorem 1, the anchor trajectory x ( t ) takes over the role of the x = 0 equilibrium of the unperturbed system (3) in the forced system (1): they are both unique uniformly bounded solutions in the ball B δ for their respective systems. Note that if f 0 and f 1 are C 1 in x, then the Lipschitz constants in the inequalities (11) can be chosen as
A unique hyperbolic anchor trajectory x ( t ) may well exist even if the strict bounds listed in (11) are not satisfied. For this reason, in the following theorem, we will simply assume that such a trajectory x ( t ) exists as a perturbation of the x = 0 fixed point and provide a recursively implementable, formal approximation for x ( t ) up to any desired order. These formulas will only assume the uniform-in-time boundedness of f 1 and its derivatives with respect to x at x = 0, without assuming the specific bounds on f 1 listed in Theorem 1. In particular, the non-autonomous term f 1 does not even have to be continuous in time for these formulas to be well-defined. This will enable predictions for anchor trajectories their SSMs in physical systems even under temporally discontinuous forcing terms or under specific realizations of bounded random forcing.
To state our recursively computable formulas for x ( t ) in our upcoming Theorem 2, we will use a matrix T R n × n whose columns comprise the real and imaginary parts of the eigenvectors of A. We order these columns in a way that T block-diagonalizes A into a stable block A s and unstable block A u,
(12)
For later purposes, we also define the time-dependent matrix G ( t ) R n × n as
(13)
Notice that if x = 0 is either asymptotically stable ( E u = { 0 } and A s = T 1 A T) or repelling ( E s = { 0 } and A u = T 1 A T), then formula (13) simplifies to
respectively. Using these quantities, we obtain the following approximation result for the anchor trajectory x ( t ).

(Computation of anchor trajectory for non-autonomous SSMs)

Theorem 2
(Computation of anchor trajectory for non-autonomous SSMs)
Assume that a unique, uniformly bounded trajectory x ( t ) U of (1) exists in a compact neighborhood U R n of x = 0, with x ( t ) perturbing from x = 0 under the forcing term f 1 ( x , t ) that satisfies (2). Assume further that f 1 has r 1 continuous derivatives with respect to x at x = 0 that are uniformly bounded in t within U. Then, for any positive integer N r, a formal expansion for x ( t ) exists in the form
(14)
where x ν ( t ) = ( x ν 1 ( t ) , , x ν n ( t ) ) = O ( f 1 U ν ) and
(15)
Here, k j i is the ith component of the integer vector k j N n { 0 } appearing in the index set,
Proof.

See Appendix A 2.

As an example of the approximation provided by Theorem 2 for x ( t ), we evaluate formula (15) up to second order ( N = 2 ) for the case wherein the x = 0 fixed point is attracting ( E u = ) for f 1 ( x , τ ) 0. In that case, we obtain
(16)
where x 2 f 0 ( 0 ) is a three-tensor and refers to the tensor product.

(Applicability to temporally discontinuous forcing)

Remark 1
(Applicability to temporally discontinuous forcing)

The asymptotic approximation (15) only requires f 1 and its x-derivatives at x = 0 to be uniformly bounded in t, as we noted earlier. No derivatives of f 1 ( x , t ) are required to exist with respect to t, and hence temporally discontinuous forcing is also covered by these formal expansions, as we will also see on a specific example with discontinuous chaotic forcing in Sec. IV A 1.

(Simplification for state-independent forcing)

Remark 2
(Simplification for state-independent forcing)

In applications to structural vibrations, the non-autonomous forcing term in system (1) often arises from external forcing that does not depend on x, i.e., f 1 ( x , t ) f 1 ( t ). In that case, the second term in formula (15) is identically zero and the summands in the first term are well defined as long as f 0 ( x ) and its derivatives are bounded at x = 0.

We first recall available SSM results for the autonomous system (3). Let E be a spectral subspace, as defined in Eq. (8). Following the terminology of Haller et al.,12 we call E a like-mode spectral subspaceif the sign of Re λ j is the same for all λ j spect ( A | E ). Otherwise, we call E a mixed-mode spectral subspace.

We call a like-mode spectral subspace E externally nonresonant if
(17)
for any choice of m k N with k = 1 n m k 2 and for any choice of λ j. It turns out that if the condition (17) is satisfied for m k coefficients up to order k = 1 n m k = σ ( E ), where σ ( E ) is the spectral quotient defined by Haller and Ponsioen,1 and then the same condition will also be satisfied for all choices of m k with k = 1 n m k 2 (Cabré et al.11).

If E is an externally nonresonant, like-mode subspace, then E has a unique, smoothest nonlinear continuation in the form of a primary spectral submanifold (or primary SSM), denoted W ( E ) C . This was deduced by Haller and Ponsioen1 in this context from the more abstract, general results of Cabré et al.11 The primary SSM W ( E ) is an invariant manifold of system (3) that is tangent to E at the origin and has the same dimension as E, as shown in Fig. 1 for f 1 ( x , t ) 0 in the extended phase space.

FIG. 1.

Left: the geometry of the autonomous SSM W ( E ) in the extended phase space of the ( x , t ) variables. Shown is the primary (smoothest SSM) tangent to the spectral subspace E at the fixed point at x = 0. Also shown is a fast SSM (blue), defined as the spectral submanifold that is tangent to the direct sum of the eigenspaces outside E. Right: the anchor trajectory x ( t ) and the time-dependent SSM W ( E , t ) in the full, non-autonomous system.

FIG. 1.

Left: the geometry of the autonomous SSM W ( E ) in the extended phase space of the ( x , t ) variables. Shown is the primary (smoothest SSM) tangent to the spectral subspace E at the fixed point at x = 0. Also shown is a fast SSM (blue), defined as the spectral submanifold that is tangent to the direct sum of the eigenspaces outside E. Right: the anchor trajectory x ( t ) and the time-dependent SSM W ( E , t ) in the full, non-autonomous system.

Close modal

In general, an infinite family, W ( E ), of SSMs exists for system (3) with the same properties, but W ( E ) is the unique smoothest among them: the other manifolds in W ( E ) are no more than σ ( E )-times differentiable (see Haller et al.12). A computational algorithm for W ( E ) was developed by Ponsioen et al.,13 then extended to general, finite-element-grade problems in second-order ODE form by Jain and Haller.16 The latest version of the latter algorithm with further extensions is available in the open-source MATLAB live script package SSMTool (see Jain et al.18).

If the spectrum of A is fully nonresonant, i.e.,
(18)
then an invariant manifold family W ( E ) tangent to E at x = 0 exists for any choice of E, whether or not E is of like-mode or mixed-mode type. While W ( E ) may contain just a single manifold [e.g., when E is the stable subspace, E s or the unstable subspace, E u, defined in (9)], there will be infinitely many members in W ( E ) for more general choices of E. In the latter case, generic members of W ( E ) are secondary (or fractional) SSMs, which have a finite order of differentiability consistent with their fractional-powered polynomial representations derived explicitly by Haller et al.12 

Note that the nonresonance conditions (17)–(18) are less restrictive than what is customary in the nonlinear vibration literature. Indeed, conditions (17)–(18) are only violated if both the real and the complex parts of the eigenvalues satisfy simultaneously exactly the same resonance relationship. Also note that 1:1 resonances involving eigenvalues with nonzero real parts do not violate (17)–(18). We will, however, only allow a 1:1 resonance within E to guarantee the normal hyperbolicity of E for the existence of W ( E ). (More specifically, a 1:1 resonance between the spectrum of A within E and outside E would render our upcoming assumption (19) to fail for any ρ > 1.) If a 1:1 resonance arises in a given application, one can simply enlarge E to include all resonant modes. This, in turn, removes any issue with the 1:1 resonance.

We now turn to the existence of non-autonomous SSMs emanating from the uniformly bounded hyperbolic trajectory x ( t ), whose existence and approximation were discussed in Theorems 1 and 2. We will focus on slow SSMs (also called pseudo-unstable manifolds), which are continuations of d-dimensional, ρ-normally attracting (like-mode or mixed-mode) spectral subspaces E of linearized system (7), i.e., can be written as
(19)
for some integer ρ > 1. Therefore, E is spanned by the k modal subspaces carrying the slowest decaying solution families of system (7), including possibly some families that do not even decay but grow. If such unstable modal subspaces are present in the direct sum (19), then E is a mixed-mode spectral subspace which we seek to continue into a mixed-mode non-autonomous SSM in system (1). If, in contrast, only stable modal subspaces are present in the direct sum (19), then E is a like-mode spectral subspace which we seek to continue into a like-mode non-autonomous SSM in system (1).

Note that the third condition in (19) always holds for arbitrary ρ > 1 if Re λ k and Re λ k + 1 have different signs. If Re λ k and Re λ k + 1 have the same sign, then the third condition in (19) holds only if the decay exponents of solutions outside E are at least ρ-times stronger than those of solutions inside E. This ρ will then determine the maximal smoothness that we can a priori guarantee for the SSM emanating from E without further assumptions. Formal expansions for such SSMs will, however, indicate higher degrees of smoothness under appropriate nonresonance conditions.

The following theorem gives our main result on the existence of non-autonomous SSMs associated with a spectral subspace E satisfying (19). We will use the term locally invariant manifold when referring to a manifold carrying trajectories that can only leave the manifold through its boundary (see, e.g., Fenichel38).

(Existence of non-autonomous SSMs)

Theorem 3
(Existence of non-autonomous SSMs)
Assume that conditions (11) of Theorem 1 are satisfied in a ball B δ R n of radius δ > 0 around x = 0. Assume further that f 0 , f 1 C r ( B δ ) for some r > 1 and the uniform bound,
(20)
holds for some finite constant K 3 > 0. Assume finally that E is a ρ-normally hyperbolic spectral subspace with 1 < ρ r that satisfies the nonresonance conditions,
(21)
Then,
  • There exists a non-autonomous spectral submanifold W ( E , t ) B δ of class C ρ that has the same dimension as E, contains x ( t ) and acts as a locally invariant manifold for system (1).

  • The SSM W ( E , t ) is as smooth in any parameter as system (1).

Proof.

See Appendix B 1. We sketch the geometry of the non-autonomous SSM, W ( E , t ), for non-vanishing f 1 ( x , t ) in Fig. 1.

(Uniqueness and smoothness of non-autonomous SSMs)

Remark 3
(Uniqueness and smoothness of non-autonomous SSMs)

The manifolds discussed in Theorem 3 are generally non-unique. In principle, an argument involving the smoothness of the Lyapunov and dichotomy spectra under generic perturbation (see Son39) could be invoked to deduce sharper smoothness results for these manifolds from linearization (see Yomdin30). These results, however, require the identification of the full Lyapunov spectrum of the linear part of the non-autonomous linear differential equation (B1) in our Appendix B 1 in order to exclude resonances. This is generally challenging for equations and unrealistic for datasets. As an alternative, Theorem 4 below will provide unique Taylor expansions for W ( E , t ) under explicitly verifiable nonresonance conditions. These expansions are varied on each persisting non-autonomous SSM up to the degree of smoothness of the SSM. We conjecture that, just as in the time-periodic and time-quasiperiodic case treated by Cabré et al.11 and Haro and de la Llave,24 a unique member of the W ( E , t ) family of manifolds will be smoother than all the others. Our Taylor expansion will then approximate those unique, smoothest manifolds at orders higher than the spectral quotient Σ ( E ) defined in Haller and Ponsioen.1 The remaining manifolds can be constructed via time-dependent versions of the fractional expansions identified in Haller et al.12 

(Non-autonomous SSMs without anchor trajectories)

Remark 4
(Non-autonomous SSMs without anchor trajectories)

For stable hyperbolic fixed points in the f 1 ( x , t ) 0 limit, an alternative to the proof in Appendix B 1 for Theorem 3 can also be given. This alternative uses the theory of non-compact normally hyperbolic invariant manifolds (see Eldering40) coupled with the “wormhole” construct of Eldering et al.41 that enables the handling of inflowing-invariant normally attracting invariant manifolds. Following the steps of our proof in Appendix C 2 for the adiabatic case, this alternative proof yields results similar to those in Theorem 3 but does not rely on a persisting anchor trajectory x ( t ) near the origin. As a consequence, it can also capture non-autonomous SSMs in weakly damped physical systems for higher forcing levels at which x ( t ) is already destroyed. This is because the strength of hyperbolicity of the unforced fixed point at x = 0 (measured by | Re λ 1 | 1) is generally much weaker than the strength of hyperbolicity of the unforced SSM, W ( E ) ( measured by | Re λ k + 1 | | Re λ k | > 1 ) in weakly damped systems. Such persisting SSMs without a stable hyperbolic anchor trajectory have been well-documented in equation- and data-driven studies of time-periodically forced systems. Those SSMs are signaled by overhangs near resonances in the forced response curves at higher forcing levels (see, e.g., Jain and Haller16 and Cenedese et al.19).

As our examples will show, W ( E , t ) will generally persist even under f 1 ( x , t ) perturbations that are significantly larger than those allowed by the rather conservative assumptions of Theorem 1. In addition, W ( E , t ) will also be smoother than C ρ under additional nonresonance conditions. To this end, we will next derive numerically implementable, recursive approximation formulas that are valid for W ( E , t ) as long as it persists.

To state these approximations, we first introduce a small perturbation parameter ϵ 0 with which we rescale the non-autonomous term in (1) as
(22)
in order to focus on moderate values of f 1 ( x , t ). As a consequence of this scaling, the expansion (14) for the anchor trajectory is also rescaled to
(23)
given that the form of the coefficients x ν ( t ) in the formulas (15) yields
(24)
Second, we let P = [ e 1 , , e n ] C n contain the complex unit eigenvectors corresponding to the ordered eigenvalues (5) of A. Then, under a coordinate change x ( u , v ) C d × C n d defined as
(25)
we obtain a complex system of ODEs,
(26)
where
(27)
and
(28)
In system (26), the fixed point ( u , v ) = ( 0 , 0 ) corresponds to the anchor trajectory x ϵ ( t ) and the u coordinate space is aligned with the spectral subspace E.
Third, using the integer multi-index ( k , p ) N d × N with | ( k , p ) | 0, we define a (now ϵ-dependent) order- | ( k , p ) | approximation to x ϵ ( t ) and the evaluation of f ^ ( u , v , ϵ ; t ) on this approximation as
(29)
respectively, where
(30)
In short, x ϵ ( t ; k , p ) is an approximation of x ( t ) in the scaled variable ϵ x up to order | k | + p 1. Accordingly, f ^ ( u , v , ϵ ; t ; k , p ) is an approximation of f ^ ( u , v , ϵ ; t ) that uses x ϵ ( t ; k , p ) instead of x ( t ).
Finally, for any multi-index k N d, we define the time-dependent diagonal matrix G k ( t ) C ( n d ) × ( n d ) as
(31)
for = d + 1 , , n. Using these quantities, we can state the following results.

(Computation of non-autonomous SSMs and their reduced dynamics)

Theorem 4
(Computation of non-autonomous SSMs and their reduced dynamics)
Assume that after the rescaling (22), an ϵ-dependent SSM, W ϵ ( E , t ), of the type described in Theorem 3 exists in the coordinates ( u , v ) for system (26) for all ϵ [ 0 , ϵ ]. Assume further that W ϵ ( E , t ) is N-times continuously differentiable for any fixed t and the nonresonance conditions,
(32)
hold for all j = k + 1 , , n and for all m k N with
(33)
Then, for all ϵ [ 0 , ϵ ],:
  1. The SSM W ϵ ( E , t ) admits a formal asymptotic expansion,
    (34)
    The uniformly bounded h k p ( t ) coefficients in this expansion can be computed recursively from their initial conditions,
    (35)
    via the formula,
    (36)
    where the functions M k p are defined as
    (37)
  2. The reduced dynamics on W ϵ ( E , t ) is obtained by restricting the u-component of system (26) to W ϵ ( E , t ), which yields
    (38)
    Here, Q u C d × n is a matrix whose jth row is e ^ j / ( e ^ j e j ) for j = 1 , d, where e ^ j is the jth unit left eigenvector of P corresponding to its unit right eigenvector e j. Equivalently, Q u is composed of the first d rows of P 1.
  3. Let ( ξ , η ) T = P 1 x denote coordinates aligned with the subspace E and the direct sum of the remaining eigenspaces, respectively, emanating from the original x = 0 fixed point of system (1) for ϵ = 0. In these coordinates, the reduced dynamics (38) becomes
    (39)
Proof.

See Appendix B 2.

(Related results)

Remark 5
(Related results)

Pötzsche and Rassmussen42 derive Taylor approximations for a broad class of invariant manifolds emanating from fixed points of non-autonomous ODE. Due to the generality of that setting, the resulting formulas are less explicit than ours, assume global boundedness on the nonlinear terms, and also assume that the origin remains a fixed point for all times (and hence disallow additive external forcing). Nevertheless, under further assumptions and appropriate reformulations, those formulas should ultimately yield results equivalent to ours when applied to SSMs. Although it does not directly overlap with our results here, we also mention related work by Pötzsche and Rassmussen43 on the computation of classic and strong stable/unstable manifolds, as well as center-stable/unstable manifolds, for non-autonomous difference equations.

(Extension to discontinuous forcing)

Remark 6
(Extension to discontinuous forcing)

The derivatives in the definition of (37) can be evaluated using the same higher-order chain-rule formula of Constantine and Savits44 that we applied in the proof of Theorem 2. Using these formulas reveal that h k p ( t ) can be formally computed as long as f 1 and its x-derivatives at x = 0 remain uniformly bounded in t. Therefore, the formal expansion we have derived for the SSM W ϵ ( E , t ) is valid for small enough ϵ even if f 1 is discontinuous in the time t. We illustrate the continued accuracy of these expansions for discontinuous forcing on an example in Sec. IV A 1.

(Relation to the case of periodic or quasi-periodic forcing)

Remark 7
(Relation to the case of periodic or quasi-periodic forcing)

The recursively defined coefficient vectors M k p ( s , h j m ( s ) ) defined in Eq. (37) are explicitly time dependent but otherwise satisfy the same formulas as the terms on the right-hand side of the invariance equations arising from autonomous SSM calculations. Consequently, the recursive formulas originally implemented in SSMTool by Jain et al.18 for autonomous SSM calculations apply here without modification. The difference is that those autonomous M k coefficients are used in solving linear algebraic systems of equations for h k, as opposed to the non-autonomous M k p coefficients are used in linear ODEs for h k p ( t ).

(Accuracy of asymptotic SSM formulas)

Remark 8
(Accuracy of asymptotic SSM formulas)

The accuracy of the reduced-order dynamics (38) increases with increasing | ( k , p ) |. There is no general convergence result under N in the reduced dynamics (38), but such convergence is known if f ( x , t ) is analytic and has periodic or quasi-periodic time dependence (Cabré et al.11 and Haro and de la Llave24). For those types of time dependencies, the O ( ϵ ) term in the reduced dynamics (38) is the same as that used for time-periodic and time-quasiperiodic SSMs in earlier publications (see, e.g., Breunung and Haller45 and Cenedese et al.20).

By formula (39), in the ξ coordinate aligned with E and emanating from x = 0, the reduced dynamics up to first order in ϵ is of the form
(40)
where we used formulas (16) and (24) in evaluating x ~ 1 ( t ), the first coefficient in the expansion of x ϵ ( t ) in Eq. (23).
In prior treatments of time-periodically and time-quasi-periodically forced SSMs, the second line on the right-hand side of Eq. (40) has been justifiably dropped in leading-order approximations, because the factor D f 0 ( P ( ξ h 0 ( ξ ) ) ) is small close to the origin (see, e.g., Breunung and Haller,45 Ponsioen et al.,15 Jain and Haller,16 and Cenedese et al.,19). In that case, the leading-order correction to the unforced reduced dynamics is simply the third line in (40), which is just the projection of the restriction of the forcing term to the autonomous SSM onto the spectral subspace E. If, in addition, the forcing term only depends on time (i.e., x f ~ 1 0), as is often the case in structural vibrations, then a leading-order approximation for small-amplitude motions on the SSM W ϵ ( E , t ) becomes
(41)
This level of approximation of the reduced dynamics was shown to be accurate for small ϵ and | ξ | for periodic and quasi-periodic forcing by the references cited above. Equation (41) now establishes the same approximation formula for general, x-independent forcing f 1 ( t ) = ϵ f ~ 1 ( t ) for trajectories that stay close to the origin. Under general forcing, however, trajectories may well stray away from the origin, given that W ϵ ( E , t ) itself will generally move substantially away from the origin. In that case, the next level of approximation obtained from Eq. (40) is
(42)
Higher-order approximation can be systematically developed by Taylor—expanding the reduced dynamics (39) in ϵ and utilizing the terms from the expansions for x ϵ ( t ) and W ϵ ( E , t ) using Theorems 2 and 4.

Alternative parametrizations of W ϵ ( E , t ) beyond the graph-style parametrization worked out here in detail are also possible (see Haller and Ponsioen,1 Haro et al.,27 and Vizzaccaro et al.35). One option is the normal-form style parameterization which simultaneously brings the reduced dynamics on W ϵ ( E , t ) to normal form. Normal form transformations can also be applied separately once W ϵ ( E , t ) is located (see, e.g., Cenedese et al.19).

We now consider slowly varying non-autonomous dynamical systems of the form
(43)
for some integer r 1. We obtain such a system, for instance, when we replace the forcing function f 1 ( x , t ) in system (1) with its slowly varying counterpart, f 1 ( x , ϵ t ), in which case we have f ( x , ϵ t ) = A x + f 0 ( x ) + f 1 ( x , ϵ t ) . The form in (43) is, however, more general than this specific example.
Introducing the slow variable α = ϵ t, we rewrite system (43) as
(44)
For ϵ = 0, we assume that the x-component of system (44) has an α-dependent, uniformly bounded, and uniformly hyperbolic fixed point x 0 ( α ), i.e., with the notation,
(45)
we have
(46)
for some constant K > 0. We further assume that x 0 ( α ) is a class C r function of α, in which case, by the first equation in (46), we have
(47)
By the hyperbolicity assumption (46), we can order the spectrum of A ( α ) as
so that it satisfies
(48)
for some j 1. Conditions (46) imply that for ϵ = 0, the slow-fast system (44) has a one-dimensional, normally attracting, non-compact (but uniformly bounded) invariant manifold of the form
(49)
The manifold L 0 is composed of fixed points of system (44) for ϵ = 0 and hence is called a critical manifold in the language of geometric singular perturbation theory (Fenichel46).

As in the non-autonomous case treated in Sec. II, we first discuss the continued existence of an anchor trajectory that acts as a continuation of the critical manifold L 0.

(Existence and computation of anchor trajectory for adiabatic SSMs)

Theorem 5
(Existence and computation of anchor trajectory for adiabatic SSMs)
For ϵ > 0 small enough, there exists a unique, attracting, one-dimensional slow manifold L ϵ, composed of a slow trajectory x ϵ ( α ) that is uniformly bounded in α R. The trajectory x ϵ ( α ) is O ( ϵ ) C 1-close and C r diffeomorphic to the critical manifold L 0 defined in (49). Finally, for any non-negative integer N r, x ϵ ( α ) can be approximated as
(50)
with the recursively defined coefficients,
(51)
where the index set p s ( ν , γ ) is defined as
Proof.

See Appendix C 1.

Specifically, for N = 2, formula (51) gives an approximation for the slow anchor trajectory x ϵ ( α ) in the form
(52)
We seek to construct an adiabatic SSM anchored along the slow invariant manifold L ϵ that is spanned by x ϵ ( α ) for small ϵ > 0. To this end, we assume that the first k eigenvalues λ 1 ( α ) , , λ k ( α ) in the list (48) have negative real parts and there is a nonzero spectral gap between Re λ k ( α ) and Re λ k + 1 ( α ),
(53)
We also assume that the real spectral subspace E ( α ) formed by the corresponding first k eigenvalues varies smoothly in α and hence has a constant dimension,
We then define the integer part ρ of the spectral gap associated with E ( α ) as
(54)
Note that ρ 1 holds by assumption (53), with ρ measuring the minimum of how many times the attraction rates toward E ( α ) overpower the contraction rates within E ( α ).
For ϵ = 0 and for each α R, system (44) has an SSM W 0 ( E ( α ) ) tangent to E ( α ) at x 0 ( α ) (see Haller and Ponsioen1). As a consequence,
is a ( d + 1 )-dimensional, ρ-normally attracting invariant manifold for system (44) for ϵ = 0 in the sense of Fenichel38 and Eldering et al.41 More specifically, M 0 is a ρ-normally attracting, non-compact, inflowing-invariant manifold with a boundary, as shown in Fig. 2 for ϵ = 0 .
FIG. 2.

Left: the geometry of the critical manifold L 0 (spanned by the anchor trajectory x 0 ( α )), the attracting invariant manifold M 0, and the fast stable manifold N 0 (blue) of L 0. Two trajectories off these manifolds are shown in red. Right: the persistent slow manifold L ϵ (spanned by x ϵ ( α )) and and the adiabatic SSM M ϵ for small enough ϵ > 0.

FIG. 2.

Left: the geometry of the critical manifold L 0 (spanned by the anchor trajectory x 0 ( α )), the attracting invariant manifold M 0, and the fast stable manifold N 0 (blue) of L 0. Two trajectories off these manifolds are shown in red. Right: the persistent slow manifold L ϵ (spanned by x ϵ ( α )) and and the adiabatic SSM M ϵ for small enough ϵ > 0.

Close modal

We have the following main results on the persistence of the manifold M 0 in the form of an adiabatic SSM for small enough ϵ.

(Existence of adiabatic SSMs)

Theorem 6
(Existence of adiabatic SSMs)

Assume that x 0 ( α ) and f ( x , α ) have r continuous derivatives that are uniformly bounded in α R in a small, closed neighborhood of M 0. Let m = min ( r , ρ ). Then, for ϵ > 0 small enough, there exists a persistent invariant manifold (adiabatic SSM anchored along L ϵ), denoted M ϵ, that is C m diffeomorphic to M 0, has m uniformly bounded derivatives and is O ( ε ) C 1 -close to M 0. Furthermore, L ϵ M ϵ holds.

Proof.

See Appendix B 2.

We show the geometry of the adiabatic SSM, M ϵ, in Fig. 2 for ϵ > 0.

(Applicability to mixed-mode adiabatic SSMs)

Remark 9
(Applicability to mixed-mode adiabatic SSMs)

For the purposes of constructing M ϵ, we had to assume in Eq. (53) that the spectrum of A ( α ) lies in the negative complex half plane, i.e., E ( α ) contains only stable directions. In the terminology of Haller et al.,12 this means that W ( E ( α ) ) has to be a like-mode SSM for all values of α for Theorem 6 to apply. This enables us to select an inflowing-invariant, normally attracting manifold M ~ 0 in our proof to which the wormhole construct in Proposition B1 of Eldering et al.41 is applicable. There is every indication that M ϵ also exists under the weaker assumption (46), which only requires A ( α ) to have no eigenvalues on the imaginary axes, and hence allows W ( E ( α ) ) to be a mixed-mode SSM that contains both stable and unstable directions. In this case, however, an M ~ 0 with an inflowing or overflowing boundary cannot be chosen and hence Proposition B1 of Eldering et al.41 would need a technical extension that we will not pursue here, although it appears doable. We simply note that the asymptotic formulas we will derive for M ϵ are formally valid under the weaker assumption (46) as well [i.e., for mixed-mode W ( E ( α ) )], as long as E ( α ) is normally hyperbolic, i.e., attracts solutions at a rate that is uniformly stronger than the contraction rates inside E ( α ). We will evaluate and confirm these formulas in an example with a mixed-mode adiabatic SSMs in Sec. IV B 2.

As in the general non-autonomous case treated in Sec. II, the adiabatic M ϵ will generally persist for larger values of ϵ and will also be smoother than C m under additional nonresonance conditions. In the following, we will provide a numerically implementable recursive scheme for computing M ϵ, assuming that it exists and is smooth enough.

To state these recursive approximation results, we first introduce the new coordinates,
(55)
where P ( α ) C n diagonalizes the matrix A ( α ). The coordinates ( u , v ) C d × C n d align with the d-dimensional stable spectral subspace family E ( α ) and with the direct sum of the remaining eigenspaces of A ( α ), respectively. In these new coordinates, system (44) becomes
(56)
where
(57)
and
(58)
By existing SSM theory, any α = c o n s t . slice of M 0 can be written in the form of a Taylor expansion in u with coefficients depending smoothly on α. Therefore, the whole of M 0 can be written as
(59)
We can now state the following results on the computation of the adiabatic SSM, M ϵ, for ϵ 0.

(Computation of adiabatic SSMs)

Theorem 7
(Computation of adiabatic SSMs)
Assume that an ϵ-dependent adiabatic SSM, M ϵ, of the type described in Theorem 6 exists in the transformed system (56) for all ϵ [ 0 , ϵ ] . Assume further that M ϵ is N-times continuously differentiable and the nonresonance conditions,
(60)
hold for all j = k + 1 , , n and for all m k N with
(61)

Then, for all ϵ [ 0 , ϵ ].

  1. The adiabatic SSM, M ϵ, admits a formal asymptotic expansion
    (62)
    The functions h k p ( α ) are uniformly bounded in α R for all ( k , p ) and can be computed recursively from their initial conditions,
    (63)
    via the formula
    (64)
    where
    (65)
    Specifically, for ϵ = 0, the coefficients in expansion (59) for M 0 can be computed as
    (66)

  1. The reduced dynamics on M ϵ is given by
    (67)
    Here, Q u ( α ) C d × n is a matrix whose jth row is e ^ j ( α ) / ( e ^ j ( α ) e j ( α ) ) for j = 1 , d, where e ^ j ( α ) is the jth unit left eigenvector of P ( α ) corresponding to its unit right eigenvector e j ( α ). Equivalently, Q u ( α ) is composed of the first d rows of P 1 ( α ).

Proof.

See Appendix C 3.

A leading-order approximation for the reduced dynamics on the adiabatic SSM, M ϵ, can be obtained from formula (67) as
(68)
where we have used the expression for x 1 ( α ) from expansion (52).

We now assume that the slow time dependence in system (43) arises from slow (but generally not small) additive forcing. This is a reasonable assumption, for instance, in the setting of external control forces acting on a highly damped structure, whose unforced decay rate to an equilibrium is much faster than the rate at which the forcing changes in time. Another relevant setting is very slow (quasi-static) forcing in structural dynamics.

In such cases, system (43) can be rewritten as
(69)
or, equivalently,
(70)
For ϵ = 0, the invariant manifold L 0 of fixed points satisfies
(71)
and we have
(72)
Note that we have changed our notation slightly for the matrix A to point out that it depends solely on x 0 ( α ).

These simplifications imply that the matrix A ( x 0 ( α ) ), the column matrix P ( x 0 ( α ) ) of the right eigenvectors of A ( x 0 ( α ) ) and the row matrix Q u ( x 0 ( α ) ) of the first d, appropriately normalized left eigenvectors of A ( x 0 ( α ) ) can now be determined solely from the unforced part f 0 ( x ) of the right-hand side of system (69) along any parametrized path x 0 ( α ) defined implicitly by Eq. (71). Similarly, an inspection of the formulas (66) defining the slow SSM, M 0, reveals that the graph h 0 ( u , x 0 ( α ) ) of M 0 can be computed along the path x 0 ( α ) solely based on the knowledge of f 0 ( x ); only the path itself depends on the specific forcing term f 1 ( α ).

From Eq. (67), we obtain that the leading-order reduced dynamics on M ϵ are now
(73)
Along any envisioned path x 0 ( α ), one can a priori compute the parametric family of frozen-time reduced models,
(74)
at all points p within an open set U in the phase space that is expected to contain x 0 ( α ) for possible forcing functions f 1 ( α ) of interest. Then, for any specific forcing f 1 ( α ), we can determine the path x 0 ( α ) from Eq. (71) and use the appropriate elements of the model family (74) with p = x 0 ( α ) in Eq. (73) to obtain the non-autonomous leading-order model,
(75)
In practice, one would only precompute (74) at a discrete set of points { p k } k = 1 K U to obtain a finite model family from (74) and then interpolate from these points to approximate the full leading-order model (73) along x 0 ( α ).
A specific feature arising in applications to forced mechanical systems is that the phase space variable x = ( q , v ) is composed of positions q R n / 2 and their corresponding velocities v = q ˙ R n / 2, where n is an even number denoting twice the number of degree of freedom of the mechanical system. In that case, forcing only appears in the v component of the ODE (69), implying
As a consequence, Eq. (71) takes the more specific form
(76)
leading to the single equation,
for the path x 0 ( α ) = ( q ( α ) , 0 ) . This means that the reduced-order model family (74) can be constructed a priori along different spatial locations q of the configuration space under the application of static forces f 1 v ( α ) that keep the mechanical system at equilibrium ( v = 0) at those q ( α ) locations. This is a significant simplification over the general setting of Eq. (74), because it only requires the construction of the model family (74) over a codimension- n subset U q = { ( q , v ) R n : v = 0 } of the full phase space R n of system (69).

This simplified reduced model construction is expected to provide further improvement over current autonomous-SSM-based control strategies that use only a single member of the model family (74), computed at fixed point x 0 of the unforced system (see Alora et al.36,37). This single model is then used along the full path x 0 ( α ) to approximate the leading-order model (73). Clearly, this approximation will generally only be valid when the instantaneous equilibrium path x 0 ( α ) remains close to the unperturbed fixed point x 0 . We will explore the application of the ideas described in this section to robot control in other upcoming publications.

Here, we consider two main mechanical examples to illustrate the construction of anchor trajectories and the corresponding non-autonomous SSMs under weak or adiabatic forcing. The governing equations of these mechanical systems are of the general form
(77)
where x is the vector of generalized coordinates; M = M T is the positive definite mass matrix; and C = C T is the positive definite damping matrix. The stiffness matrix K = K T is positive definite in our first example and indefinite in our second example; f ( x ) is the vector of geometric nonlinearities; and F ( t ) is the vector of external forcing. We can rewrite the second-order system (77) in the first-order form we have used in this paper by letting
With this notation, system (77) will take the form of the first-order system (1) or that of (69), depending on whether F ( t ) can be considered small or slow.
Both of our examples treated below admit a hyperbolic fixed point at x = 0 in the absence of forcing. To add general aperiodic forcing to these systems, we will construct the vector F ( t ) from a trajectory on the chaotic attractor of the classic Lorenz system,
(78)
with parameter values ρ = 28 , σ = 10 , and β = 8 3 . We solve this system over the time interval [ 0 , 500 ], starting from the initial condition ( x 0 , y 0 , z 0 ) = ( 0 , 0.3 , 0.5 ) to generate weak forcing and over the time interval [ 15 , 20 ] from ( x 0 , y 0 , z 0 ) = ( 0.8 , 0.3 , 0.2 ) to generate slow forcing. Specifically, we use the x ( t ) component of this solution as forcing profile after appropriate scaling in both cases. Outside the time interval [ 0 , 500 ], we select the forcing to be identically zero for the weak forcing case. To emulate slow forcing, we use a slowed-down version of the chaotic signal by rescaling time as t α = ϵ t. These choices of the forcing ensure its uniform boundedness, which was one of our fundamental assumptions in deriving our results for weak forcing. In the slow forcing case, the signal is smooth enough to ensure accurate numerical differentiability in the α-interval [ 0 , 6 ] we consider.

We formulated our results in Secs. II and III for fully non-dimensionalized systems for which smallness or slowness can simply be imposed by selecting a single perturbation parameter ϵ small enough. In specific physical examples, such a non-dimensionalization can be done in multiple ways and is often cumbersome to carry out for large systems. One nevertheless needs to assess whether the external forcing is de facto smaller or slower than the magnitude or the speed of variation of internal forces along trajectories. While such a consideration has been largely ignored in applications of perturbation results to physical problems, here we will propose and apply heuristic physical measures that help assessing the magnitude and the slowness of the perturbation.

To assess smallness or slowness in a systematic and non-dimensional fashion without non-dimensionalization of the full mechanical system, we first express the mechanical system in the form
(79)
with the subscripts referring to the autonomous internal forces and non-autonomous external forces, respectively. Using these forces, we define the non-dimensional forcing weakness r w and forcing speed r s as
(80)
where the over-bar represents averaging with respect to the initial conditions ( x 0 , x ˙ 0 ) of the unperturbed (i.e., ϵ = 0 ) trajectories [ x ( t ) , x ˙ ( t ) ] of the system over an open domain of interest in the phase space.

In a practical setting, we deem a specific external forcing weak if r w 1 or slow if r s 1. This is the range in which the perturbation-based invariant manifold techniques used in proving our relevant theorems can reasonably expected to apply. We have noted, however, that the SSMs obtained from these techniques are normally hyperbolic and hence persist under further increases in r w and r s. Indeed, we will show that our formal asymptotic SSM expansions tend to remain valid for r w > 1 and and r s > 1. In some cases, we also find these formulas to have predictive power for r w , r s 1.

All MATLAB scripts used in producing the results for the examples below can be publicly accessed at Kaundinya and Haller.47 

1. Weak forcing

We consider the two-degree-of-freedom (2 DOF) mechanical system from Haller and Ponsioen,1 placed now on a moving cart subject to external chaotic shaking, as shown in Fig. 3(a). In this example, the fixed point of the unforced system is asymptotically stable. As a result, under external forcing, the SSM attached to a unique nearby anchor trajectory will be of like-mode type.

FIG. 3.

(a) The physical setup for the forced cart problem. (b) The chaotic forcing signal g ( t ), generated as the x ( t ) component of the numerically solved Lorenz system (78) with initial condition ( 0 , 0.3 , 0.5 ), computed over the time interval [ 0 , 500 ] measured in seconds. For times below t < 0, we set g ( t ) 0.

FIG. 3.

(a) The physical setup for the forced cart problem. (b) The chaotic forcing signal g ( t ), generated as the x ( t ) component of the numerically solved Lorenz system (78) with initial condition ( 0 , 0.3 , 0.5 ), computed over the time interval [ 0 , 500 ] measured in seconds. For times below t < 0, we set g ( t ) 0.

Close modal
Using the coordinate vector x = ( q 1 , q 2 , x c ), we obtain the equations of motion in the form (77) with
(81)
where M T = M f + m 1 + m 2. We further set m 1 = m 2 = 1 ( kg ), M f = 4 ( kg ), k = k f = 1 ( N / m ), c = c f = 0.3 ( Nm / s ), and γ = 0.5 ( N / m 3 ).

We will consider two different cases of forcing by scaling the magnitude of the chaotic signal g ( t ) in two different ways. In the first case, we will have max | F ( t ) | = 0.06 ( N ), which gives the non-dimensional forcing weakness r w = 0.16 from the first formula in Eq. (80). This forcing scheme can, thus, be considered weak albeit not very small. Our subsequent choice of max | F ( t ) | = 3 ( N ) gives r w = 7.8, which is definitely outside the range of small forcing.

As the unforced x = 0 fixed point is asymptotically stable and the forcing in uniformly bounded, this fixed point perturbs for small enough | M T g ( t ) | into a nearby attracting anchor trajectory x ( t ) by Theorem (1). From the asymptotic formula (14) listed in the theorem, a 5th-order approximation for x ( t ) has the terms,
(82)
with the external forcing f 1 ( t ) = ( 0 , 0 , 0 , 0 , 0 , g ( t ) ) T appearing in the first-order version of the equation motion.

Utilizing that the forcing signal g ( t ) vanishes for t < 0, we evaluate the integrals in (82) via the trapezoidal integration scheme of MATLAB. For sufficient accuracy, we use forcing values evaluated at 10 4 equally spaced times to ensure accurate numerical integration. The terms x ν ( t ) in the expansion for x ( t ) are then spline-interpolated in time to obtain an expression for x ( t ) for arbitrary t [ 0 , 500 ].

To compute the non-autonomous SSM coefficients, we change coordinates such that x ( t ) serves as the origin of the transformed system (26). As a result, the matrices defined in Eq. (27) are of the form
with λ 1 , 2 = 0.0227 ± 0.3956 i, λ 3 , 4 = 0.1234 ± 1.2700 i, and λ 5 , 6 = 0.3788 ± 1.6707 i. These eigenvalues satisfy the non-resonance condition (32), and hence by Theorem 4, there exists a 2D autonomous SSM, W ϵ ( E , t ), that admits a truncated Taylor expansion,
(83)
for any N 2 in system (26). Here, we choose the dimensional book-keeping parameter,
in computing expansion (83). As already noted, the actual ratio between external forces to internal forces along trajectories is better represented by the dimensionless forcing weakness parameter r w.
Up to the order of truncation in Eq. (83), the two-dimensional SSM-reduced order model (38) can be written in the complex coordinate u 1 as
(84)
After computing the trajectory from formula (82) up to fifth order, we use the recursive formulas in statement (i) of Theorem 4 to compute the coefficients h j m ( t ) in Eq. (84) up to the same order recursively. Details for the coefficients in the SSM-reduced model (84) can be found in the code available from Kaundinya and Haller.47 

To assess the performance of the truncated SSM-reduced model (84) for N = 5, we first plot the normalized mean trajectory error computed from ten different initial conditions for max | | F | | = 0.06 ( N ) in Fig. 4(a). We compute the dependence of this error on max | | F | | in Fig. 4(b). As expected, the errors grow with max | | F | | but remains an order of magnitude less than max | | F | |.

FIG. 4.

Error distribution of the SSM-reduced modeling error for the chaotically shaken cart example. (a) Ensemble average (black) of the mean trajectory errors (grey, normalized by the maximum of the maximal trajectory amplitude) computed from ten initial conditions under forcing with max | | F | | = 0.06 ( N ). (b) The mean trajectory error under varying max | | F | | on trajectories starting from the initial condition u 1 = u 2 = 1.2 in all cases.

FIG. 4.

Error distribution of the SSM-reduced modeling error for the chaotically shaken cart example. (a) Ensemble average (black) of the mean trajectory errors (grey, normalized by the maximum of the maximal trajectory amplitude) computed from ten initial conditions under forcing with max | | F | | = 0.06 ( N ). (b) The mean trajectory error under varying max | | F | | on trajectories starting from the initial condition u 1 = u 2 = 1.2 in all cases.

Close modal

Probing larger forcing amplitudes, we still obtain accurate reduced-order models up to max | | F | | = 3 ( N ). As an illustration, Fig. 5 uses the center-of-mass coordinate x c to show our 5th order asymptotic approximation using Theorem 1 for the anchor trajectory of the chaotic SSM (blue), a simulation from the full model for a trajectory in the 2D SSM converging to the anchor point (black), and a prediction from the 2D SSM-reduced model obtained from Theorem 4 for the same trajectory (red).

FIG. 5.

SSM-reduced model predictions and their verification in the chaotically forced cart example (see the text for details). The non-dimensional forcing weakness measure gives r w = 7.8, which is clearly outside the small forcing regime. We also zoom in to show the initial transient phase to highlight how the reduced order model captures the dynamics of the full-order model.

FIG. 5.

SSM-reduced model predictions and their verification in the chaotically forced cart example (see the text for details). The non-dimensional forcing weakness measure gives r w = 7.8, which is clearly outside the small forcing regime. We also zoom in to show the initial transient phase to highlight how the reduced order model captures the dynamics of the full-order model.

Close modal

In Fig. 6 (Multimedia available online), we also plot snapshots of the chaotic 2D SSM, its anchor trajectory, a trajectory from the full-order model, and the prediction for this trajectory from the SSM-reduced order model at three different times. Note how both the full trajectory and the predicted model trajectory synchronize with the anchor trajectory of the SSM (see Fig. 6). Finally, we would like to illustrate that the improper integrals in our asymptotic formulas for anchor trajectories and SSMs indeed converge and give correct results even for discontinuous forcing, as noted in Remark 6. To this end, we make the forcing discontinuous at ten random points in time, as shown in Fig. 7(b). The same formulas for the anchor trajectory, its attached SSM, and the SSM-reduced order model remain formally applicable and continue to give accurate approximations even for the relatively large forcing amplitude of max | | F | | = 3 ( N ), as seen in Fig. 7.

FIG. 6.

(a)–(c) Snapshots of the evolution of the 2D SSM, its anchor trajectory, and a trajectory of the SSM-reduced model shown in the coordinates used in Theorem 4. The non-dimensional forcing weakness measure is again r w = 7.8, which is no longer small. See Movie 1 in the supplementary material for the full evolution. Multimedia available online.

FIG. 6.

(a)–(c) Snapshots of the evolution of the 2D SSM, its anchor trajectory, and a trajectory of the SSM-reduced model shown in the coordinates used in Theorem 4. The non-dimensional forcing weakness measure is again r w = 7.8, which is no longer small. See Movie 1 in the supplementary material for the full evolution. Multimedia available online.

Close modal
FIG. 7.

(a) Discontinuous chaotic forcing profile for the cart example with non-dimensional forcing weakness r w = 7.8. We set g ( t ) 0 for t < 0. (b) Same as Fig. 5 but for the discontinuous chaotic forcing profile shown in the subplot (a).

FIG. 7.

(a) Discontinuous chaotic forcing profile for the cart example with non-dimensional forcing weakness r w = 7.8. We set g ( t ) 0 for t < 0. (b) Same as Fig. 5 but for the discontinuous chaotic forcing profile shown in the subplot (a).

Close modal
We close this subsection with a slight modification of our example that adds more nonlinearity to the problem and underlines the utility of the higher-order approximations we have developed for the anchor trajectory (or generalized steady state) and the SSM attached to it. For brevity, we will only illustrate the need for higher-order approximations to the anchor trajectory by replacing the localized nonlinearity in Eq. (81) with the more general nonlinearity,
(85)
This makes the originally linear right-most spring between the cart and the wall in Fig. 3(a) nonlinear with cubic nonlinear stiffness coefficient γ f.

First, Fig. 8(a) shows for γ f = 0 ( N / m 3 ) (i.e., for the linear limit of the right-most spring), the O ( 1 ) and O ( 11 ) asymptotic approximations to the anchor trajectory of the chaotic SSM using the results of Theorem 1. Already in this case, the O ( 1 ) approximation is noticeably improved by the O ( 11 ) approximation, but the O ( 1 ) approximation is also close to the actual generalized steady state.

FIG. 8.

Leading-order and higher-order approximation of the generalized steady state of the cart system under weak chaotic forcing. (a) For max | | F | | = 3 ( N ) with non-dimensional forcing weakness r w = 15.82 and a linear right-most spring [ γ f = 0 ( N / m 3 )]. (b) For one-fifth of the forcing level in the case (a) [ max | | F | | = 0.6 ( N ), r w = 6] but for a nonlinear right-most spring [ γ f = 0.5 ( N / m 3 )].

FIG. 8.

Leading-order and higher-order approximation of the generalized steady state of the cart system under weak chaotic forcing. (a) For max | | F | | = 3 ( N ) with non-dimensional forcing weakness r w = 15.82 and a linear right-most spring [ γ f = 0 ( N / m 3 )]. (b) For one-fifth of the forcing level in the case (a) [ max | | F | | = 0.6 ( N ), r w = 6] but for a nonlinear right-most spring [ γ f = 0.5 ( N / m 3 )].

Close modal

In contrast, Fig. 8(b) shows a case in which much smaller forcing is applied but the right-most spring is nonlinear with γ f = 0.5 ( N / m 3 ). In this case, more significant error develops between the actual anchor trajectory and its O ( 1 ), while the O ( 11 ) approximation remains indistinguishably close to the actual anchor trajectory. This example, therefore, underlines the need for higher approximations to the anchor point in the presence of stronger nonlinearities. As the SSM is attached to the anchor trajectory, the accuracy of an SSM-reduced model also depends critically on this refined approximation.

2. Slow forcing

For the same cart system shown in Fig. 3(a), we select the physical parameters m 1 = m 2 = 1 ( kg ), M f = 2 ( kg ), k = k f = 1 ( N / m ), and c = c f = 0.3 ( Ns / m ). We now slow down the previously applied chaotic forcing by replacing the external forcing vector in Eq. (81) with F ( α ) = ( 0 , 0 , N s g ( α ) ) T. We select this forcing to be of large amplitude by letting N s = 10 and consider the phase range α [ 0 , 6 ]. In our analysis, we will start with the minimal forcing speed ϵ = α ˙ = 0.001, for which obtain the non-dimensional forcing speed r s = 0.72 from the second formula in (80). This qualifies as moderately slow external forcing relative to the speed of variation of internal forces along unperturbed trajectories. In contrast, for the maximal forcing speed ϵ = α ˙ = 0.008, we will consider, we obtain r s = 7.26, which can no longer be considered slow. Again, of particular interest will be how our asymptotic formulas for adiabatic SSMs and their anchor points perform in the latter case, which is clearly faster than what is normally considered slowly varying in perturbation theory.

We start by solving for the fixed point x 0 ( α ) of the system under static forcing, as defined by the algebraic equation in Eq. (46). In the present example, this amounts to solving the equation,
(86)
which is obtained from the general forced equation of motion (79). We solve this nonlinear algebraic equation for the selected range of α values by Newton iteration.

Once x 0 ( α ) is available numerically, we evaluate the recursive formula (51) up to order N = 3 to obtain

Having computed this approximation for the anchor trajectory, we change to the coordinates defined in Eq. (55). In those coordinates, we compute the α-dependent matrices defined in Eq. (57) numerically. Their general form in the present examples is
To perturb from this frozen-time limit, we set α = ϵ Ω t and Ω = 1 Hz to make the small parameter ϵ non-dimensional.

We set N = 3 and solve for all the α dependent coefficients in the expansion of the 2D slow chaotic SSM by solving a recursive linear algebraic equations (64) order by order. We find that all 28 coefficients that describe the SSM up to cubic order are nonzero. An important numerical step used in this task is to perform differentiation with respect to α. For this purpose, we first need to re-orient the unit eigenvectors originally returned by MATLAB for specific values of α to make these eigenvectors smooth functions of α. We then perform a central finite differencing which includes four adjacent points in the α direction. We finally employ the Savitzky–Golay filtering function of MATLAB to smoothen the results obtained from finite differencing an already finite-differenced signal.

To illustrate the ultimate accuracy of these numerical procedures, Table I shows the normalized mean trajectory error for predictions from two different SSM-reduced models. We see that under increasing ϵ (i.e., faster forcing), the errors still remain an order-of-magnitude smaller than ϵ.

TABLE I.

Normalized mean trajectory error for predictions from two different SSM-reduced models under different ϵ values. The corresponding non-dimensional forcing speed rs is also shown for each listed value of ϵ.

O ( ϵ , j = 1 , | m | = 3 ) O ( ϵ 3 , j + | m | = 3 )
ϵ = 0.001 (rs = 0.72)  9.8 × 10−5  7.2 × 10−6 
ϵ = 0.008(rs = 7.26)  6.4 × 10−3  7.8 × 10−4 
ϵ = 0.010 (rs = 8.66)  1.1 × 10−2  2.9 × 10−3 
O ( ϵ , j = 1 , | m | = 3 ) O ( ϵ 3 , j + | m | = 3 )
ϵ = 0.001 (rs = 0.72)  9.8 × 10−5  7.2 × 10−6 
ϵ = 0.008(rs = 7.26)  6.4 × 10−3  7.8 × 10−4 
ϵ = 0.010 (rs = 8.66)  1.1 × 10−2  2.9 × 10−3 

As we did for weak forcing, we track the center of mass coordinate x c to test the accuracy of our asymptotic formulas up to order N = 3 for ϵ = 0.008 in Fig. 9. Note how this cubic-order approximation of the SSM-reduced model already tracks the full solution closely. Finally, as in the weakly forced case, we also plot in Fig. 10 snapshots of the evolution of the slow anchor trajectory, the attached chaotic SSM, a simulated full-order trajectory, and its prediction from the same initial condition based on the SSM-reduced dynamics (Multimedia available online). These snapshots further confirm the accuracy of our analytic approximation formulas, now in several coordinate directions.

FIG. 9.

(a) Slow chaotic forcing signal g ( ϵ t ) for the cart example for ϵ = 0.008, which amounts to r s = 7.26. (b) Assessment of the accuracy of the anchor trajectory and the reduced dynamics of the slow chaotic SSM in the cart example, based on the history of the x c ( t ) coordinate of the center of mass.

FIG. 9.

(a) Slow chaotic forcing signal g ( ϵ t ) for the cart example for ϵ = 0.008, which amounts to r s = 7.26. (b) Assessment of the accuracy of the anchor trajectory and the reduced dynamics of the slow chaotic SSM in the cart example, based on the history of the x c ( t ) coordinate of the center of mass.

Close modal
FIG. 10.

(a)–(c) Snapshots of the evolution of the slow chaotic SSM, its anchor trajectory, the true solution of the full system, and the SSM-reduced model prediction for the true solution of the chaotically shaken cart, all represented in the physical coordinates [ q 2 , x c , x ˙ c ]. The non-dimensional forcing speed is r s = 7.26. See Movie 2 in the supplementary material for the full evolution. Multimedia available online.

FIG. 10.

(a)–(c) Snapshots of the evolution of the slow chaotic SSM, its anchor trajectory, the true solution of the full system, and the SSM-reduced model prediction for the true solution of the chaotically shaken cart, all represented in the physical coordinates [ q 2 , x c , x ˙ c ]. The non-dimensional forcing speed is r s = 7.26. See Movie 2 in the supplementary material for the full evolution. Multimedia available online.

Close modal

1. Weak forcing

The physical setup of our second example, a particle moving on a shaken rail with a bump in the middle, is shown in Fig. 11(a). A major difference from our previous mechanical example is that the origin of the unforced system is now an unstable, saddle-focus-type fixed point. Therefore, under weak chaotic forcing, a nearby chaotic SSM of mixed-mode type is expected to arise attached to a saddle-focus-type chaotic anchor trajectory.

FIG. 11.

(a) Setup for the chaotically forced bumpy real example. (b) Plots of 16 initial conditions on local mixed-mode non-autonomous 2D SSM, released at four different times. The left inset shows how the transient behavior of the full system matches those of the anchor trajectories calculated from our theory. The right inset shows fast convergence of the full solution to the anchor trajectories together with local predictions of the reduced order model. The non-dimensional forcing weakness measure in this example is r w = 44.6.

FIG. 11.

(a) Setup for the chaotically forced bumpy real example. (b) Plots of 16 initial conditions on local mixed-mode non-autonomous 2D SSM, released at four different times. The left inset shows how the transient behavior of the full system matches those of the anchor trajectories calculated from our theory. The right inset shows fast convergence of the full solution to the anchor trajectories together with local predictions of the reduced order model. The non-dimensional forcing weakness measure in this example is r w = 44.6.

Close modal
We use the center of mass position x c and the relative position x of the smaller mass m as generalized coordinates. In terms of these coordinates, the quantities featured in the general equation of motion (77) take the form

with for the coordinate vector x = ( x c , x ) T. The parameter β controls the height of the bumpy rail and a the distance to the wells. The parameter c models linear damping of the mass m due to the rail. The velocity-dependent nonlinearity f ( x , x ˙ ) is the result of expanding the exact equations of motion with respect to the height of the rail and keeping only the leading-order nonlinearities.

We set m = 1 ( kg ), M f = 4 ( kg ), k f = 1 ( N / m ), c f = c = 0.3 ( Ns / m ), a = 0.3 ( m ), β = 1 5 a 3, and g = 9.8 ( m / s 2 ). The forcing signal g ( t ) will be the same as in the cart example. We scale the chaotic forcing signal in a way that the non-dimensional forcing weakness measure defined in formula (80) gives r w = 44.6 when evaluated on the phase space region containing the trajectories used in our analysis. Therefore, the external force here is large relative to the internal forces of the system.

This forcing magnitude is outside the small range in which Theorem 4 strictly guarantees the existence of a saddle-type anchor trajectory for a mixed-mode chaotic SSM in this problem. We nevertheless evaluate our asymptotic expansions for the anchor trajectory, as those expansions hold as long as the SSM smoothly persist under increasing forcing. Specifically, a cubic-order approximation x ( t ) ν = 1 3 x ν ( t ) for the anchor trajectory x ( t ) from formula (14) requires the functions,
with Green’s function G ( t ) defined in formula (13) for saddle-type anchor trajectories. The external forcing appears now as f 1 ( t ) = ( 0 , g ( t ) , 0 , 0 ) T in the first-order formulation of the equation of motions. We calculate these integrals using the same numerical methodology as in our first cart example. In this example, the calculation of G ( t τ ) from (13) involves the matrices
where the eigenvalues of the unforced saddle-focus-type fixed point at the origin are λ 1 = 5.5196, λ 2 = 5.9094, and λ 3 , 4 = 0.0301 ± 0.4465 i. The 2D non-autonomous SSM can be constructed as a graph over the mixed-mode tangent space corresponding to the real eigenvalues λ 1 and λ 2.

Without listing the details here, we also perform a similar calculation for the two attracting anchor trajectories created by the chaotic forcing from the two asymptotically stable fixed points of the unperturbed system that lie at the two bottom points of the bumpy rail. These two stable anchor trajectories will represent two chaotic attractors for the forced system. Trajectories initialized near the unstable anchor trajectory on its attached 2D mixed-mode SSM will converge to one of these chaotic attractors. This is indeed observable in Fig. 11(b), in which we launch several trajectories near the unstable anchor trajectory (red) within the 2D mixed-mode SSM approximated up to cubic order. These trajectories then converge to one of the two predicted attracting chaotic anchor trajectories (blue) which we have computed up to cubic order of accuracy from our formulas.

The instability of the anchor trajectory perturbing from the origin makes it challenging to verify the accuracy of our local expansion (64) for the 2D mixed-mode SSM attached to this trajectory. We can nevertheless verify the local accuracy for these formulas by tracking nearby trajectories launched from the predicted 2D SSM and confirm in Fig. 12 (Multimedia available online) that those trajectories evolve together with the predicted SSM until ejected from a neighborhood of the unstable anchor trajectory.

FIG. 12.

(a) and (b) SSM dynamics in the (not so) weakly forced bumpy rail example. Shown are snapshots of the predicted 2D mixed-mode SSM (grey), its anchor trajectory (red), trajectories predicted by the SSM-reduced model (green), and full system simulations of trajectories launched from the same initial conditions (black). The non-dimensional forcing weakness measure is again r w = 44.6. See Movie 3 in the supplementary material for the full evolution. Multimedia available online.

FIG. 12.

(a) and (b) SSM dynamics in the (not so) weakly forced bumpy rail example. Shown are snapshots of the predicted 2D mixed-mode SSM (grey), its anchor trajectory (red), trajectories predicted by the SSM-reduced model (green), and full system simulations of trajectories launched from the same initial conditions (black). The non-dimensional forcing weakness measure is again r w = 44.6. See Movie 3 in the supplementary material for the full evolution. Multimedia available online.

Close modal

2. Slow forcing

In the same bumpy rail problem but now with height β = 1 10 3, we rescale the forcing in time to make it slowly varying. The slow forcing vector is given by F ( α ) = ( N s g ( α ) 0 ) with N s = 3 and the range of interest is α [ 0 , 6 ] . We will perform simulations for the forcing speed ϵ = α ˙ = 0.01, which yields the non-dimensional forcing speed measure r s = 1.22, signaling moderately fast forcing.

As in our first example, we compute the zeroth order term x 0 ( α ) in the expansion of the slow anchor trajectory for the SSM by solving an algebraic equation of the form (86) by Newton iteration,
Up to cubic order, expansion (51) for the slow, unstable anchor trajectory perturbing from the origin under the slow forcing can be written as
With this approximation for the anchor trajectory at hand, we can change to the coordinates defined in Eq. (55). In those coordinates, we compute the α-dependent matrices defined in Eq. (57) now arise as
where we again set α = ϵ Ω t with Ω = 1 Hz, as in our first example.

As in the weakly forced version of this example, we plot the predictions from the 2D SSM-reduced order model and compare it to simulations of the full system in Fig. 13(a) (Multimedia available online). We again see close agreement of the reduced dynamics with the true solution in the vicinity of the slow, unstable anchor trajectory, prior to the convergence of the trajectories to the two stable, chaotic anchor trajectories arising from the two stable unforced equilibria along the rail. Finally, as in our first example, we show snapshots of the projected 2D slow, mixed-mode SSM along with some trajectories of its restricted dynamics. These again stay close locally to trajectories of the full system that are launched from the same initial conditions, thereby confirming the expectations put forward in Remark 9.

FIG. 13.

(a)–(c) Same as Fig. 12, but now for moderately slow chaotic forcing with non-dimensional forcing speed measure r s = 1.22. See Movie 4 in the supplementary material for the full evolution. Multimedia available online.

FIG. 13.

(a)–(c) Same as Fig. 12, but now for moderately slow chaotic forcing with non-dimensional forcing speed measure r s = 1.22. See Movie 4 in the supplementary material for the full evolution. Multimedia available online.

Close modal

In this paper, we have derived existence results for spectral submanifolds (SSMs) in non-autonomous dynamical systems that are either weakly non-autonomous or slowly varying (adiabatic). While our exact proofs for these SSMs results only cover weak enough or slow enough time dependence, the invariant manifolds we obtain are structurally stable and hence are persistent away from these limits. Under further nonresonance conditions, the manifolds and their reduced dynamics also admit asymptotic expansions up to any finite order for which we present recursive formulas. These formulas suggest that one of the non-autonomous SSMs is as smooth as the original dynamical system, just as primary SSMs are known to be in the autonomous case reviewed in Sec. I. MATLAB implementations of our recursive formulas for the examples treated here are available from Kaundinya and Haller.47 

In our results, weakly forced SSMs are guaranteed to exist under uniformly bounded forcing. This restriction is expected because unbounded perturbations to an autonomous limit of a dynamical system will generally wipe out any structure identified in that limit. For adiabatic SSMs, the uniformity of the forcing magnitude is replaced by the uniformity of the strength of hyperbolicity along the anchor trajectory of the SSM in the limit of frozen time.

In the weakly forced case, the SSMs emanate from unique, uniformly bounded hyperbolic anchor trajectories. Our expansions for these anchor trajectories are of independent interest as they provide formal approximations of the generalized steady state response of any nonlinear system subject to moderate but otherwise arbitrary time-dependent forcing. In finite-element simulations under aperiodic forcing, such steady states have been assumed to exist but their computation tends to involve lengthy simulations of randomly chosen initial conditions until they have reached a (somewhat vaguely defined) statistical steady state.

Even with today’s computational power, such simulations are still prohibitively long due to the small damping of typical structural materials, the high degrees of freedom of the models used, and the costly evaluations of nonlinearities arising from complex geometries and multi-physics. With the explicit recursive formulas we obtained here, one can now directly approximate these generalized steady states without lengthy simulations.

For these asymptotic formulas to converge, the forcing does not actually have to be uniformly bounded: it may grow temporally unbounded in a compact neighborhood of the origin as long as the improper integrals in the statements of Theorems 2 and 4 converge. To the best of our knowledge, such explicit, readily computable asymptotic expansions for time-dependent steady states have been unavailable in the literature even for time-periodic or time-quasiperiodic forcing, let alone time-aperiodic forcing. Weak or slow time-periodic and time-quasi-periodic forcing is also covered by these formulas as special cases.

Our existence results for SSMs assume temporal smoothness for the dynamical system and hence do not strictly cover the case of random forcing with continuous paths. Our approximation formulas, however, are formally applicable (i.e., the improper integrals in them converge) for much broader forcing types, including mildly exponentially growing or temporally discontinuous forcing. Therefore, the results in this paper provide formal nonlinear model-reduction formulas that can be implemented as approximations for all forcing types arising in practice.

The simple examples considered here already illustrate the ability of our asymptotic formulas to work under various forcing types, including chaotic and discontinuous forcing. Our weakly non-autonomous SSM theory is expected to be helpful in model reduction for structural vibration problems in which so far time-periodic and time-quasiperiodic SSMs have been constructed. An application of our adiabatic SSM theory is already under way in the model-predictive control of soft robots, wherein the target trajectories of the robot generally have a much slower time scale than the highly damped robot itself.

Going from the simple, illustrative examples treated here to finite-element problems with temporally aperiodic forcing will require no further theoretical development. Indeed, the existence and smoothness results derived in this paper are valid in arbitrarily high (finite) dimensions. However, a computational reformulation will be required to make these results directly applicable to second-order mechanical systems without the need to convert them to first-order systems with diagonalized linear parts. The development of such a reformulation using the approach of Jain and Haller16 is currently under way.

We are grateful to Christian Pötzsche for very helpful and detailed explanations and references on the continuity of the dichotomy spectrum for non-autonomous differential equations. We are also grateful to Martin Rasmussen and Peter Kloeden for helpful references on non-autonomous invariant manifolds.

The authors have no conflicts to disclose.

George Haller: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Methodology (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). Roshan S. Kaundinya: Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (supporting); Writing – review & editing (equal).

The codes required to generate these numerical results are openly available at https://github.com/haller-group/Aperiodic-SSMs, Ref. 47.

1. Proof of Theorem 1: Existence and regularity of a hyperbolic anchor trajectory

We first recall a result of Palmer,29 reformulated from Hartman,48 for a general non-autonomous ODE of the form
(A1)
where A ( t ) and f ( x , t , p ) are C 0 functions of their arguments; | f ( x , t , p ) | μ is assumed to hold for a constant μ > 0 and for all x R n and t R; f ( , t , p ) admits a global Lipschitz constant L > 0 for all times t; and p R m is a parameter vector.
Assume that the linearization
(A2)
of system (A1) has an exponential dichotomy, i.e., there exists a fundamental matrix solution Φ ( t ) of (A2), constants K , κ > 0 and a projection P (with P 2 = P ) such that
(A3)
This condition guarantees that the x = 0 solution of the linearized system (A2) is hyperbolic, i.e., has well-defined stable and/or unstable subbundles (no neutrally stable directions). Finally, assume that
(A4)
Under these assumptions, there exists a unique, globally bounded solution x ( t ; p ) of the nonlinear system (A1) that satisfies
(A5)
as shown in Lemma 1 of Palmer.29 Furthermore, x ( t ; p ) is continuous in the parameter p. [If f is smooth in its arguments, then the continuity of x ( t ; p ) in p can be strengthened to smooth dependence on p using more general results for the persistence of non-compact, normally hyperbolic invariant manifolds; see Eldering.40 Such a strengthened version, however, would not be specific enough about the norm of | x ( t ; p ) | to the extent given in (A5).] Note that x ( t ; p ) takes over the role of x = 0 as a unique, uniformly bounded solution under nonzero f ( x , t , p ). The result of Palmer29 quoted above makes no assumption on f ( x , t , p ) containing only nonlinear terms, and hence f ( x , t , p ) is allowed to be an arbitrary perturbation to the linear ODE (A2), as long as it is uniformly bounded and uniformly Lipschitz.

The uniform boundedness conditions (A4)(A5) will not be satisfied in realistic applications. Nevertheless, one can still use the above results in such applications using smooth cutoff functions (see, e.g., Fenichel38) that make f ( x , t , p ) vanish for all t B δ outside a small ball B δ R n. The latter classic tool from invariant manifold theory is, however, only applicable here if the predicted hyperbolic trajectory x ( t ; p ) lies inside B δ, where the cutoff right-hand side and the original right-hand side of (A1) still coincide.

To ensure this, we consider a δ-ball B δ R n around x and define the uniform bound μ ( δ ) and uniform Lipschitz constant L ( δ ) as
(A6)
Assume further that these bounds satisfy
(A7)
which assures that assumption (A4) holds and also that the predicted x ( t ; p ) falls in B δ by the estimate (A5). Note that the inequalities in (A7) will always hold for δ > 0 small enough if f ( x , t , p ) is only composed of terms that are nonlinear in x. Indeed, in that case, we can select μ ( δ ) = O ( δ 2 ) and L ( δ ) = O ( δ ) for all x B δ, and hence the inequalities in (A7) are satisfied for small enough δ, given that K , κ > 0 do not depend on δ.
We now apply the above refined persistence result in the setting of the perturbed nonlinear system (1) by letting
(A8)
and assume a separate uniform Lipschitz constant L 1 ( δ ) for f 1 ( x , t ) satisfying
By our hyperbolicity assumption on the x = 0 fixed point, we can select the constant κ > 0 so that
(A9)
We can also bound within B δ the nonlinear term f 0 and its spatial derivative as
(A10)
where K 1 , K 2 > 0 are appropriate constants. Additionally, we can specifically choose K 2 = max x B δ | x f 0 ( x ) | by the mean value theorem. Given these bounds, the constants in formulas (A6) can be estimated as
(A11)
Based on the inequalities (A11), we see that assumptions (A7) are satisfied if
or, equivalently, if
(A12)
We stress again that these conditions do not have to hold globally for all x, only for x B δ, but they must hold uniformly within B δ for all times. This follows because we can apply the C cutoff procedure mentioned above to the whole right-hand side of (A1), including f 1 ( x , t ). Based on these considerations and using the inequalities in Eq. (A10), we obtain statement (i) of Theorem 1.

Statement (ii) of Theorem 1 follows from the smooth dependence of x ( t ; p ) on the parameter vector p, as shown by Palmer.29 

2. Proof of Theorem 2: Approximation of the anchor trajectory x*(t)

First, we introduce a perturbation parameter ϵ 0 and rewrite the full non-autonomous system (1) as
(A13)
By the smooth dependence of the uniformly bounded hyperbolic solution x ( t ) on parameters, we can seek this solution in the form of a Taylor expansion
(A14)
with Taylor coefficients ξ ν ( t ) that are uniformly bounded in time. By statement (ii) of Theorem (1), such a formal Taylor expansion in ϵ is justified up to any finite order, although may not necessarily converge.
Substitution of the expansion (A14) into Eq. (A13) gives
(A15)
Equating equal powers of ϵ on both sides yields the system of differential equations
(A16)
Note that the term formally containing ξ ν ( t ) in D ϵ ν f ( x ϵ ( t ) , t ; ϵ ) | ϵ = 0 is
because f 0 ( x ) = O ( | x | 2 ) and x 0 ( t ) 0. Therefore, the right-hand side of Eq. (A16) only depends on ξ 1 ( t ) , , ξ ν 1 ( t ), which makes the whole system of equations a recursively solvable sequence of inhomogeneous, linear, constant-coefficient system of ODEs. The recursive solutions are of the form
(A17)
Assume first, for simplicity, that the x = 0 fixed point is asymptotically stable for ϵ = 0, and hence
(A18)
In that case, by the uniform boundedness of ξ 1 ( t 0 ) for all t 0 R, we can take the t 0 limit in (A17) to obtain
(A19)
To obtain a more explicit recursive formula for ξ ν ( t ) for general ν that is suitable for direct numerical implementation, we will use the multi-variate Faá di Bruno formula of Constantine and Savits44 for higher-order derivatives of composite functions. To recall the general form of this formula, we first consider a general composite function H : R p R, defined as
(A20)
and introduce the non-negative multi-index ν and related notation as
We also introduce an ordering relation on N p for arbitrary μ , ν N p such that ν μ holds provided one of the following is satisfied:

  • | μ | < | ν |,

  • | μ | = | ν | and μ 1 < ν 1, or

  • | μ | = | ν | and μ 1 = ν 1 , μ k = ν k , μ k + 1 < ν k + 1 for some 1 k < p.

Finally, for any μ , γ N p and s N +, we define the index set
With this notation, Constantine and Savits44 prove the following multi-variate version of the Faá di Bruno formula:
(A21)
where k j i denotes the ith element of the multi-index k j N m { 0 }.
Of relevance to us is the case p = 1, wherein we have H = f q, m = n , g i = j 1 ϵ j ξ j i ( τ ), i = 1 , , n; x = ϵ R , x 0 = 0 R, ν = ν N, i = i N, and y = x R n. In that case, we can write
(A22)
Note, however, that H ( ϵ ) also has explicit dependence on ϵ and hence is not exactly of the form (A20). To address this issue, we also observe that
and hence H 0 ( ϵ ) and H 1 ( ϵ ) are individually of the form (A20). Applied to these two functions, the formula (A21) simplifies to
Substitution of these expressions into formula (A19) then gives the final recursive formulas
(A23)
where k j i is the ith component of the integer vector k j N n { 0 } appearing in the index set
We have also used the notational convection | γ | x 1 γ 1 x n γ n = I for γ = 0.
Substitution of (A23) into the expansion (A14) gives
(A24)
These results cover anchor trajectories of like-mode SSMs of stable hyperbolic fixed points but do not cover anchor trajectories for mixed-mode SSMs of such fixed points or anchor trajectories of like-mode SSMs of unstable hyperbolic fixed points.
We now extend formula (A23) to mixed-mode SSMs by weakening the stability assumption (A18) back to our general hyperbolicity assumption (6). In this case, the matrix A has an exponential dichotomy, as described by the inequalities (10). The dichotomy exponent κ > 0 can be selected as in (A9). We use the matrix T defined in (12) in the coordinate transformation
which brings system system (1) to the form
with ( f s ( y , t ) , f u ( y , t ) ) T = T 1 f ( T y , t ) and
(A25)
In these new coordinates, with the notation
the system of ODEs (A16) becomes
(A26)
whose solutions can be written as
(A27)
Based on the spectral properties of A s and A u listed in (A25) and the uniform boundedness of ξ ^ ν s , u ( t 0 s , u ), we can take the limits t 0 s and t 0 u + to obtain
(A28)
Therefore, introducing Green’s function (13), we can write the solution ξ ν ( t ) in the y coordinates as
This then leads to the final formula (15) for a general hyperbolic trajectory x ( t ), in analogy with formula (A24) for an asymptotically stable hyperbolic trajectory.

1. Proof of Theorem 3: Existence of non-autonomous spectral submanifolds

Under the conditions of Theorem 1, we can shift coordinates by letting
which transforms system (1) to the form
(B1)
where
(B2)
Next, we introduce a perturbation parameter 0 ϵ < δ and the scalings
(B3)
These scalings imply
With these expressions, we can rewrite Eq. (B1) as
We can further rewrite these equations as
(B4)
(B5)
While linearization results do not apply to Eq. (B4) due to the non-hyperbolicity of the origin, the invariant manifold results of Kloeden and Rassmussen49 do apply. Their main condition stated in our context is
(B6)
of which the first one is already satisfied and the second one requires the local uniform boundedness of the first derivatives of G ( Y , t ) in a neighborhood of Y = 0. By inspection of formula (B5), we see that the second condition in (B6) is satisfied if the following four conditions are fulfilled:

  • | x f 0 ( x ( t ) ) | is uniformly bounded,

  • | x f 1 ( x , t ) | and | x 2 f 1 ( x , t ) | are uniformly bounded in the B δ neighborhood of x = 0,

  • | g ( y , t ) | is uniformly bounded in a neighborhood of y = 0, and

  • y | g ( y , t ) | is uniformly bounded in a neighborhood of y = 0.

Given the definition of g ( y , t ) in (B2), conditions (c)–(d) are satisfied if | x f ( x ( t ) , t ) | and
(B7)
are uniformly bounded in B δ.

First, note that condition (a) is satisfied because x ( t ) is uniformly bounded. Second, since x ( t ) stays in the B δ ball under the conditions (11), we obtain that | x f 0 ( x ( t ) + y ) x f 0 ( x ( t ) ) | in uniformly bounded in a compact neighborhood of y = 0. Therefore, for condition (B7) to hold [and hence to ensure that (c)–(d) to hold], it remains to require that | x f 1 ( x ( t ) + y ) x f 1 ( x ( t ) ) | remain uniformly bounded in a neighborhood of y = 0. By the mean value theorem, this holds if | x 2 f 1 ( x , t ) | remains uniformly bounded in B δ, which is just condition (b) above. Therefore, in addition to our assumptions in (11), we need to assume condition (20) for (a)–(d) to hold. In summary, under the conditions (11) and (20), the local invariant manifold results of Kloeden and Rassmussen49 are applicable near the Y = 0 fixed point of system (B4)(B5).

Specifically, under conditions (11) and (20), the results in Sec. 6.3 of Kloeden and Rassmussen49 apply to general non-autonomous systems of differential equations of the form (B4). These results in turn build on classic result of Sacker and Sell50 that establish the existence of a dichotomy spectrum Σ for A ( t ) that consists up to n + 1 disjoint closed intervals of the form
(B8)
(By the definition of the dichotomy spectrum, the linear system of ODE Y ˙ = ( A ( t ) λ I ) Y admits no exponential dichotomy for any λ Σ.) Assuming m 2, one can therefore select constants κ j + , κ j R for j = 1 , , m 1 from the gaps among the closed spectral subintervals in (B8) as follows:
such that for appropriate constants K > 0 and projection maps P ± j ( t 0 ) R ( n + 1 ) × ( n + 1 ) with P ± j ( t 0 ) P ± j ( t 0 ) = P ± j ( t 0 ), the normalized fundamental matrix solution Φ ( t , t 0 ) of Y ˙ = A ( t ) Y satisfies
(B9)
for all j = 1 , , m 1.
Then, by Theorem 6.10 and Remark 6.6 of Kloeden and Rassmussen,49 for each j = 1 , , m 1, there exist two non-autonomous invariant manifolds W j ± ( t ) for system (B4) such that
  • W j ± ( t ) contain the Y = 0 fixed point of system (B4).

  • In a neighborhood U R n + 1 , the manifolds W j ± ( t ) can described as graphs with the help of C 0 functions w j ± : U × R R n + 1 such that
    (B10)
  • w j ± ( Z , t ) is uniformly o ( | Z | ), i.e., lim u 0 w j ± ( Z , t ) Z = 0, uniformly in t.

  • If
    (B11)
    holds for a positive integer m j +, then W j + ( t ) is of class C m j +. Similarly, if
    (B12)
    holds for a positive integer m j , then W j ( t ) is of class C m j .
  • If κ j + < 0 then for all γ > κ j +, we have sup t 0 Y ( t , t 0 , Y 0 ) e γ t < for all Y 0 W j + ( t 0 ) U in a small enough neighborhood U of Y = 0. Similarly, if κ j > 0 then for all γ < κ j , we have sup t 0 Y ( t , t 0 , Y 0 ) e γ t < for all Y 0 W j + ( t 0 ) U in a small enough neighborhood U of Y = 0.

Based on these results, further non-autonomous invariant manifolds can be obtained by letting
(B13)
Note that Theorem 6.10 and Remark 6.6 of Kloeden and Rassmussen49 assume no hyperbolicity for the Y = 0 fixed point of system (B4), which makes them applicable to the Y = 0 fixed point of the extended system (B4). Specifically, the dichotomy spectrum Σ of the matrix A ( t ) defined in Eq. (B5) is discrete and given by
As Kloeden and Rassmussen49 point out, the definition (B13) yields a hierarchy of invariant manifolds, which in turn implies the hierarchy of invariant manifolds shown in Table II for the original system (1) arising from this construct for the extended system.
TABLE II.

Hierarchy of invariant manifolds for the extended system (B4).

W 1 + ( t )  ⊂  W 2 + ( t )  ⊂  ⋅⋅ ⋅  ⊂  W n 1 + ( t )  ⊂  R n + 1 
                 
    W 2 , 1 ( t )  ⊂  ⋅⋅ ⋅  ⊂  W n 1 , 1 ( t )  ⊂  W 1 ( t ) 
                 
            ⋮    ⋮ 
            W n 1 , n 2 ( t )  ⊂  W n 2 ( t ) 
                 
                W n 1 ( t ) 
W 1 + ( t )  ⊂  W 2 + ( t )  ⊂  ⋅⋅ ⋅  ⊂  W n 1 + ( t )  ⊂  R n + 1 
                 
    W 2 , 1 ( t )  ⊂  ⋅⋅ ⋅  ⊂  W n 1 , 1 ( t )  ⊂  W 1 ( t ) 
                 
            ⋮    ⋮ 
            W n 1 , n 2 ( t )  ⊂  W n 2 ( t ) 
                 
                W n 1 ( t ) 
For any index j c, all invariant manifolds W i , j + ( t ) defined in formula (B13) are graphs over a spectral subspace that contains the μ c = 0 center direction and hence they are as smooth in ϵ as they are in the other variables along that spectral subspace. This smoothness property is in turn shared by all invariant manifolds of the form

The ϵ =constant slices of W i + ( t ), W j ( t ), and W i , j ( t ) then provide similar invariant manifolds W i + ( t ; ϵ ), W j ( t ; ϵ ), and W i , j ( t ; ϵ ) of one less dimension along the x ( t ; ϵ ) trajectory of the original system (1) for small enough with ϵ > 0 under the assumed scaling (B3). The degree of smoothness of W i + ( t ; ϵ ), W j ( t ; ϵ ), and W i , j ( t ; ϵ ) in x and ϵ coincides with the general degree of smoothness that can be inferred from the above construct for their extended counterparts, W j + ( t ), W j ( t ), and W i , j ( t ). Specifically, the following non-autonomous invariant manifolds can be inferred from Table II for system (1):

  1. W i + ( t ; ϵ ): the time-dependent smooth continuations of the fastest decaying i modes for any 1 i n 1. These manifolds satisfy (B11) for arbitrary large m and hence are as smooth as the original dynamical system. They are generally C 0 in ϵ which can be improved to C r for pseudo-stable manifolds, i.e., for all W i + ( t ; ϵ ) with μ i > 0.

  2. W j ( t ; ϵ ): the time-dependent smooth continuations of the slowest decaying (or even growing) n 1 j modes. If μ j 1 / μ j > m j for some positive integer m j , then these manifolds satisfy (B12) with m j 1 and hence are of smoothness class C m j . For all j with μ j < 0, W j ( t ; ϵ ) is also C m j smooth in ϵ. Any manifold W j ( t ; ϵ ) with μ j > 0 (fast unstable manifolds) can only be concluded to be C 0 in ϵ from this argument.

  3. W i , j ( t ; ϵ ); the continuation of any (like-mode or mixed-mode) spectral subspace, that is, r i , j-normally hyperbolic. In that case, W i , j ( t ; ϵ ) is also of smoothness class C r i , j in ϵ.

2. Proof of Theorem 4: Approximation of the SSM, W (E, t), anchored at x ϵ ( t )

a. Invariance equation
Next, we derive an approximation for the non-autonomous SSMs anchored to the hyperbolic trajectory x ( t ). We perform a linear change of coordinates
where P = [ e 1 , , e n ] C n contains the complex eigenvectors corresponding to the ordered eigenvalues (5) of A, and ( u , v ) C d × C n d. Rewriting the scaled version (A13) of system (1) in these complex coordinates, we obtain
(B14)
where
(B15)
Note that under the scaling (A13), the anchor trajectory x ( t ) becomes ϵ-dependent, which is reflected by our modified notation x ϵ ( t ) for the same trajectory.
By definition, we have
therefore,
(B16)
By the definition of f ^ in (B15) and by formula (B16) with p = 1 , we also have
(B17)
For ϵ = 0, under the nonresonance conditions (17), we have a unique, primary spectral submanifold
(B18)
where h 0 ( u ) = O ( | u | 2 ) C r defines the primary SSM as a smooth graph over E with a quadratic tangency to E at x = 0 (see Haller and Ponsioen1). For ϵ > 0 small, under the conditions of statement (i) of Theorem 3, this primary SSM persists in the form of an (generally non-unique) invariant manifold
(B19)
with coefficients h k p ( t ) that are uniformly bounded in t within a B δ ball around the origin x = 0 for 0 ϵ δ, and with the notation
Note that this expansion is only valid up to the order of smoothness of the SSM. For general members of the W ϵ ( E , t ), this smoothness can only guaranteed to be class C ρ, but exceptional members may admit higher-order expansions, just as primary SSMs do in the case of autonomous, time-periodic, and time-quasiperiodic SSMs (see Haller and Ponsioen1).

Also note that for the smooth persistence of W 0 ( E ) as W ϵ ( E , t ), a possible 1:1 resonance between the an eigenvalue inside E and another one outside E has to be excluded in order to secure the normal hyperbolicity of W 0 ( E ). This is reflected by the lowering of the lower index in the summation in Eq. (33) from 2 to 1 relative to what one requires for the existence of W 0 ( E ).

Differentiating the definition (B19) of the invariant manifold W ϵ ( E , t ) in time and using the system of ODEs (B14), we obtain
(B20)
At the same time, substitution into (B14) gives
(B21)
Comparing (B20) and (B21) give the invariance PDE satisfied by W ϵ ( E , t ),
(B22)
As the origin ( u , v ) = ( 0 , 0 ) is a fixed point of system (B14) for all t and ϵ, we must have h ( 0 , t ; ϵ ) 0 for all ϵ 0 small and all t R. This, in turn, implies,
(B23)
Furthermore, for ϵ = 0, the non-autonomous SSM W ϵ ( E , t ) becomes the autonomous SSM W ( E ), and hence we must have
(B24)
with the time-independent Taylor series coefficients h k in the expansion for W 0 ( E ) in Eq. (B18).
b. Structure of the invariance equation
Substitution of (B19) into the invariance equation (B22) gives
(B25)
We observe that
(B26)
Therefore, the invariance equation (B25) can be rewritten as
(B27)
where
(B28)
Recall from Eq. (B17) that f ^ v and f ^ u vanish for u , v , ϵ = 0 and have no O ( | v | ) terms. Therefore, when | k | , p 0 h k p ( t ) u k ϵ p is substituted into the terms in the Taylor expansion of f ^ v ( u , v , ϵ ; t ), then it is multiplied in each case by at least the first power of u or at least the first power of ϵ . (When | k | , p 0 h k p ( t ) u k ϵ p is substituted into terms of O ( | v | 2 ) then it is multiplied by itself, but h 0 0 ( t ) 0, and hence the lowest order term multiplying | k | , p 0 h k p ( t ) u k ϵ p will be again at least the first power of u or at least the first power of ϵ .) As a result, M k p ( h j m , t ) will only depend on h j m that are lower in order, i.e.,
Equating coefficients of equal powers of u in Eq. (B27), we therefore obtain the recursively solvable linear system of inhomogeneous linear ODEs
(B29)
with M k p ( t , h j m ) defined in formula (37). So far, we know from Eqs. (B23) and (B24) that
(B30)
which is a homogeneous system of ODEs in (B29) for | k | > 0.
c. Solution of the invariance equation
A first observation about solving the family (B29) of linear ODEs is that for arbitrary | k | and p = 0, we obtain from Eqs. (B24) and (B29) the algebraic equations
We know from classic SSM theory that this algebraic system of equations is a recursively solvable linear system of equations, because | j | < | k | holds for all M k 0 ( h j 0 ), and hence
(B31)

Next, we note that the nonresonance assumption (18) implies that the coefficient matrix A k of the homogeneous part of Eq. (B29) is nonsingular. We additionally require now this homogenous part to admit a hyperbolic fixed point, which leads to the stronger nonresonance condition (32) of Theorem 3. This is the same non-resonance condition that arises in the work of Haro and de la Llave24 for quasiperiodic forcing, which is a subset of the general forcing class we are considering here.

Under the nonresonance condition (32), all diagonal elements of A k have nonzero real parts. We can therefore uniquely split A k into the sum of a diagonal matrix A k , which contains the stable eigenvalues of A k as well as zeros, and a diagonal matrix A k +, which contains the unstable eigenvalues of A k as well as zeros,
We define the corresponding splitting for the vectors h k p ( t ) and M k p ( t , h j m ( t ) ) as
We then then split the expression for the solution of system (B29) as
(B32)
If h k p ( t ) is a uniformly bounded solution of system (B29) for all forward and backward times, then using the signs of the nonzero diagonal entries of the matrices A k ±, we obtain from the formulas (B32) that