Optimal Dichotomy of Temporal Scales and Boundedness and Stability of Time-Varying Multidimensional Nonlinear Systems

This paper develops a new approach to the estimation of the degree of boundedness or stability of multidimensional nonlinear systems with time-dependent nonperiodic coefficients-an essential task in various engineering and natural science applications. Known approaches to assessing the stability of such systems rest on the utility of Lyapunov functions and Lyapunov first approximation methodologies, typically providing conservative and computationally elaborate criteria for multidimensional systems of this category. Adequate criteria of boundedness of solutions to nonhomogeneous systems of this kind are rare in the contemporary literature. Lately, we develop a new approach to these problems which rests on bounding the evolution of the norms of solutions to initial systems by matching solutions of a scalar auxiliary equation we introduced in [1], [2] and [3]. Still, the technique advanced in [3] rests on the assumption that the average of the linear components of the underlying system is defined by a stable matrix of general position. The current paper substantially amplifies the application domain of this approach. It is merely assumed that the time-dependent linear block of the underlying system can be split into slow and fast varying components by application of any smoothing technique. This dichotomy of temporal scales is determined by the optimal criterion reducing the conservatism of our estimates. In turn, we transform the linear subsystem with slow-varying matrix in a diagonally dominant form by successive applications of the Lyapunov transforms. This prompts the development of novel scalar auxiliary equations embracing the estimation of the norms of solutions to our initial systems. Next, we formulate boundedness or stability criteria and estimate the relevant regions of the underlying systems using analytical and abridged numerical reasoning.


Introduction
Assessment of the boundedness and stability of solutions to nonlinear systems with variable and nonperiodic coefficients is a long-standing problem that finds its origin in diverse arrays of vital applied fields that, for instance, are concerned with the design of observers and controllers. Sufficient conditions of asymptotic stability for such systems were firstly endorsed by Lyapunov [1] via the utility of the Lyapunov exponents which were widely used subsequently. In fact, it was shown in [1] that under some natural conditions on the underlying homogeneous system, its trivial solution is asymptotically stable if the linearized at zero systems are regular and its maximal Lyapunov exponent is negative, see the contemporary review and historical perspectives on this subject, for instance, in [2] and [3]. It turns out that the validation of the regularity condition is a challenging problem and the computation of Lyapunov exponents requires hefty numerical simulations for multidimensional systems. Subsequently, Perron [4] confirmed that the regularity condition is essential for the asymptotic stability of such nonlinear systems and observed that the Lyapunov exponents of a linearized system can be sensitive to small perturbations.
A more robust but conservative methodology for stability analysis of nonautonomous nonlinear systems was introduced in [5], where the concept of generalized exponents was developed. It turns out that the upper generalized exponent is equal to or larger than the maximal Lyapunov exponent, which extends the conservatism of this approach. Yet, the estimation of generalized exponents for nonautonomous and nonlinear systems presents a challenging problem as well. More relaxed stability conditions resting on the application of this methodology were given in [6].
The concept of the exponential dichotomy of solutions to time-varying systems was used in qualitative analysis and characterization of the stability of some of these systems [5]. However, practical applications of such criteria typically require comprehension of the fundamental set of solutions to the relevant linearized system.
It appears that there are no necessary and sufficient criteria of local stability of the trivial solution to a sufficiently broad class of time-varying nonlinear systems, see, e.g., [2,3,[5][6][7][8][9].
Several numerical techniques for the estimation of the stability regions of autonomous nonlinear systems were developed in the last few decades [21][22][23], but the extension of these techniques to time-varying systems presently remains elusive.
To our knowledge, the problem of estimating the trapping regions for a wide class of nonhomogeneous nonlinear systems with variable coefficients has not been ascribed in the current literature yet.
An alternative approach to the analysis of the stability of a nonlinear system was developed in [24]. This approach rests on the analysis of convergence of nearby trajectories of the initial system. Yet, to our knowledge, the application of this technique to complex nonlinear systems is still limited.
In [25], we developed a novel approach for assessing the boundedness/stability of time-varying nonlinear systems and applied it to the estimation of the trapping/stability regions of these systems. This approach rests on the development of a scalar auxiliary equation with solutions bounding from the above the norm of matching solutions to the original system. The obtained estimates appear to be sufficiently accurate if the solutions emanated from the central part of the trapping/stability regions but turned out to be more conservative near the boundaries of such regions. Subsequently, [26] combines this technique with the relevant successive approximations granting estimation of the boundaries of the trapping/stability regions and the norms of adjacent solutions with increasing precision.
Yet, the auxiliary equation developed in [25] contains the condition number of the transition matrix of a linear counterpart of the initial system. Such functions can approach infinity on large time intervals, which would make our inferences overconservative for a number of practically important systems. Pinsky [27] escapes this limitation through developing a technique casting the auxiliary equation in a modified form under the assumption that the average of the time-dependent matrix of a linear block of the given system is a stable matrix of general position.
The current paper develops a novel methodology that essentially extends the application domain of the approach developed in [27]. It merely assumes that the linear block of the given system can be split into slow-and fast-varying components by application of, e.g., the moving averaging technique. Such temporal dichotomy is defined by the optimal criterion reducing the conservatism of our subsequent estimates. In turn, we apply recursively the Lyapunov transform to map the slow linear subsystem to its diagonally dominant form which ultimately aids the development of novel scalar auxiliary equations. Solutions to these scalar equations are bound from the above time histories of the norms of solutions for a broad class of nonautonomous nonlinear systems and unfold a more inclusive assessment of their boundedness/stability properties. Lastly, we authenticate our inferences in representative simulations.

Definitions and preliminaries
This paper studies the behavior of solutions to a time-varying nonlinear system with marked linearization which can be written in the following form, where functions,F * : [t 0 , ∞) → R n , and matrix B : 1, F 0 ≥ 0, · stands for induced 2-norm of a matrix or 2-norm of a vector, and x(t, t 0 , x 0 ) : [t 0 , ∞) × T × H → R n is a solution to (1). We will frequently write below that x(t, t 0 , x 0 ) ≡ x(t, x 0 ) to simplify notation and assume that (1) possesses a unique solution ∀x 0 ∈ H and ∀t ≥ t 0 .
The last assumption together with the continuity of the right side of (1) implies that x(t, x 0 ) is a continuous and continuously differential function that is bounded ∀t ∈ [t 0 , t * ], t * < ∞ which let us focus this paper on the behavior of x(t, x 0 ) for t → ∞.
To simplify referencing of the formulas, we also write a homogeneous counterpart to (1) as,ẋ We will study the stability of the trivial solution to this equation. The next useful equation, which simplifies our further referencing, is a linear counterpart to (2) which takes the following form,ẋ At this point, we assume that Eqs. (2) and (3) possess unique solutions for ∀x 0 ∈ H and ∀t ≥ t 0 .
Clearly, the solution to (3) takes the following form, and w(t) are the transition and fundamental matrices for (3), respectively, which, due to continuity of B(t), are continuous and continuously differentiable in t.
Next, we mention that the Lyapunov exponents of solutions to (2), which are widely used to assess the exponential growth/decay of solutions to differential equations, can be defined as follows [2,5], After that we acknowledge the definition of the Lyapunov transform, v : y which can possess some desirable properties.
In order, let us bring the definition of the comparison principle [8] which is repeatedly used below. Consider a scalar differential equation, where the function g(t, u 1 ) is continuous in t and locally Lipschitz in u 1 for ∀u 1 ∈ ℘ ⊂ R. Suppose that the solution to this equation, u 1 t, u 1,0 ∈ ℘, ∀t ≥ t 0 . Next, consider a differential inequality, where D + u 2 denotes the upper right-hand derivative in t of u 2 t, u 2,0 and u 2 t, u 2,0 ∈ ℘, ∀t ≥ t 0 . Then,u 1 t, u 1,0 ≥ u 2 t, u 2,0 , ∀t ≥ t 0 .
Our subsequent inferences rely on the application of the techniques that were developed in our former papers [25,27], which are abbreviated in this section to make this paper more inclusive.
In [25] under the normalizing condition w(t 0 ) 1 we developed a scalar equation,Ẋ X : [t 0 , ∞) → R ≥0 (R ≥0 is a set of nonnegative real numbers) and show that the norm of the solution to (1) is bounded from above by the solution of (4) with matching initial conditions, i.e., x(t, x 0 ) ≤X (t, X 0 ), X 0 w −1 (t 0 )x 0 , ∀t ≥ t 0 . Note that p(t) is a continuous function and c(t) is a continuously differentiable function in t due to continuity of B(t).
To represent (4) in a tractable form, we also developed in [25] a nonlinear extension of Lipschitz continuity condition which was written as follows, where L : [t 0 , ∞)×R ≥0 → R ≥0 is continuous function in t and x , L(t, 0) 0, and x is a bounded neighborhood of x ≡ 0. Note that L(t, x ) can be readily defined in a closed-form if f * is either a piece-wise polynomial in x or can be approximated by such function with bounded in x error term [25]. In the former case, x ≡ R n and in the latter case, x ≡ R n if the error term is also globally bounded for ∀x ∈ R n . Thus, the condition, x ≡ R n is readily met in numerous nonlinear systems that are appeared in applications.
Let us show how to define L(t, x ) for a simple system which can be readily generalized to a piecewise polynomial vector-field. Assume that x [x 1 x 2 ] T , f : where we used that x n i ≤ x n ,i 1, 2, n ∈ N, N is a set of positive integers. Note that additional and more complex examples of this kind are provided in Sect. 7 as well as in [25][26][27]. With (7) at hand, (4) can be written as follows, Let us assume that (8) possesses a unique solution for ∀t ≥ t 0 and ∀z 0 ≥ 0 which implies that z(t, z 0 ) is a continuous and continuously differentiable function ∀t ≥ t 0 due to the last assumption and continuity of the right side of (8).
Clearly, due to the comparison principle, (8). Hence, the last inequality prompts the boundedness/stability criteria for solutions to the multidimensional Eqs. (1) and (2) through the relevant properties of a scalar Eq. (8).
Yet, these estimates become over conservative for relatively large values of t if lim t→∞ c(t) ∞ which frequently takes place even if A const is a stable matrix.
Next, we write (1) as follows, where G * (t) B(t) − A is a zero-mean matrix and, subsequently, rewrite the last equation in the eigenbasis of the matrix A as, where Afterward, we rearranged the last equation as follows, and use in (5) and (6) that w(t) exp(ηI + iβ(t))(t − t 0 ), which returns that p(w(t)) η and c(w(t)) 1. Then, the auxiliary Eq. (8) for (11) can be written as, where y(t, y 0 ) ≤ z(t, z 0 ), ∀t ≥ t 0 ,y ∈ y ⊂ C n , y is a neighborhood of y ≡ 0 and function L(t, z) : [t 0 , ∞) × R ≥0 → R ≥0 is developed through application of the analog of inequality (7) to function, f (t, y) v −1 f * (t, vy), which returns that.
see an illustrative example above in this section as well as more inclusive examples in Sect. 7. Note that for piecewise polynomial vector fields (13) holds for ∀y ⊂ C n . In the sequel, η was selected through the application of the condition, min η (η + α − ηI ) yielding that η α n , (η + α − ηI ) |η α n α 1 , which brings (12) into the following form, Lastly, we present a less conservative counterpart to the prior equation as follows, Afterward, we proved in [27] that where x(t, t 0 ) is a solution to (1) or (2) and z(t, z 0 ) is the matching solution to the scalar Eqs. (14) or (15).
Clearly, this approach fails if A(t 0 ) 0 and becomes over conservative if A(t 0 ) is an unstable matrix. Moreover, the conservatism of the above technique increases under the assumption that α 1 < 0, but |α 1 | is sufficiently small and sup ∀t≥t 0 Ĝ (t) is relatively large.
Next, we present for convenience some conventional definitions of trapping/stability regions, see also [27], where such definitions were also adopted.

Definition 1 A connected and compact set of all initial vectors
Clearly, this definition acknowledges that 1 is the invariant set of (1).

Definition 2
A connected and open set of all initial vectors, 2 (t 0 ), that includes zero-vector is called a region of stability of the trivial solution to (2) if the condition

Definition 3
A connected and open set of all initial vectors, x 0 ∈ 3 (t 0 ), that includes zero-vector is called a region of asymptotic stability of the trivial solution to (2) if the condition,x 0 ∈ 3 (t 0 ) implies that lim t→∞ x(t, x 0 ) 0.

Generalized auxiliary equation
This section develops a novel framework deriving a sequence of scalar auxiliary equations for our initial Eqs. (1) or (2). Our approach splits matrix B(t) into slow-and fast-varying components via the application of moving averages to the entries of this matrix. Next, we apply successively the Lyapunov transforms to (1) to map a subsystem defined by the slow-varying matrix of this equation to its simplified form. This aids the development of more general auxiliary equations escaping the limitations of our prior technique [27].
Our approach is sequenced below in a few consecutive steps where our assumptions are outlined in an ad hoc form. Then, we present a statement encapsulating our hypotheses and outcomes.
3.1. Let us split matrix B(t) into slow-and fast-varying components, B(t) a 0 (t, δ) + g 0 (t, δ), where a slow component, a 0 (t, δ) can be seized, for instance, by application to every entry of B(t) of some moving averages with fixed windows, δ { δ i j }, δ ∈ R n×n ≥0 . The components of the matrix δ are defined later in this section via the application of an optimal criterion. Clearly, the fast matrix g 0 (t, δ) B(t) − a 0 (t, δ). Note that the applications of different smoothing techniques aiding the selection of slow-varying components, like, e.g., wavelets, etc., do not alter the underlying structure of our methodology.
Next, we assume that a slow matrix a 0 (t) is a matrix of general position [28] which is continuously differentiable for an appropriate number of times that is detailed subsequently. Such matrices possess simple eigenvalues and are diagonalizable for all but possibly some isolated values of t. Such isolated values of t, for which a 0 (t) is not diagonalizable, we call irregular in this paper. Complementary values of t correspond to simple eigenvalues of a 0 (t) and are called regular. Thus, the eigenvalues and eigenvectors of a 0 (t) are continuously differentiable for a relevant number of times at regular values of t, see, e.g., [29,30].
Revision of our technique for irregular values of t is undertaken in Sect. 6, whereas other sections of this paper assume that all values of t ≥ t 0 are regular.
Under the assumption that the matrix a 0 (t) is continuously differentiable, the Lyapunov transform x v 1 (t)y 1 (t) maps (1) into the following equation, is a continuous matrix due to our assumption on a 0 (t). Hence, 1 (t), v 1 (t) andv 1 (t) can be considered as slow matrices, which implies that a 1 (t) is a slow and diagonally dominant matrix since sup t 0 ≤t v 1 (t) is a relatively small value.
To develop a consecutive approximation, we assume that matrix a 0 (t) is twice continuously differentiable, a matrix a 1 (t) possesses simple eigenvalues for ∀t ≥ t 0 and, hence, is diagonalizable by application of Lyapunov transform This conveys the following equation, In the sequel, we assume that matrix a 0 (t) is k-times continuously differentiable, matrix a k−1 (t) possesses simple eigenvalues for ∀t ≥ t 0 , and v k (t) is eigenvector matrix of a k−1 (t). Then, the application of Lyapunov transform y k−1 (t) v k (t)y k (t), k ≥ 2 brings the following equation, is a continuous matrix. Lastly, we rewrite the former equation as follows,

At this point, let us develop a scalar auxiliary equation for Eq. (17) by extending a technique that was abbreviated in the previous section
Firstly, we write (17) alike (11) as follows, , and a continuous function η k−1 : [t 0 , ∞) → R that will be defined subsequently. Next, we apply to the last equation a technique that derives the auxiliary equation for (11). For this sake, we set that a diagonal matrix (5) and (6) Then, as prior, we select η k−1 (t) through the application of the following condition Lastly, this comprises the following scalar auxiliary equation, where functions L k : [t 0 , ∞) × R ≥0 → R ≥0 are continuous in both variables. Note that L k is defined through the application of inequality (13) to a continuous function f k (t, y k ) which yields that where as prior, L k (t, y k ) can be defined in the closed-form ∀y k ∈ C n if, e.g., f k (t, y k ) are piecewise polynomials in y k , see illustrative example in Sect. 2 and a more inclusive example in Sect. 7 as well as [25,27].
Next, as previously, we write a more efficient counterpart to Eq. (15) as follows, To simplify further referencing, we present a homogeneous counterpart to (19) as follows,ż 3.3. Subsequently, we recall that all entries of the right sides of Eqs. (18)-(20) depend upon entries of matrix δ controlling separation of temporal scales of the matrix B(t). While the above equations are defined for ∀δ ∈ R n×n ≥0 , the efficacy of the estimates they provide can be enhanced by appropriate choice of δ. To this end, to abate the influence of the nonlinear component in the right sides of our auxiliary equations in a sufficiently small neighborhood of z k ≡ 0, we assume that and scalars ς k,2 > 1. Next, we define δ to maximize the degrees of stability of either equation For this sake, we define Lyapunov exponents of solutions to these equations as follows, and select δ to conform one of the following relations inf 0≤δ≤ φ K ,1 (δ) or inf 0≤δ≤ φ K ,2 (δ), where entries of matrix ⊂ R n×n ≥0 assign some practically feasible intervals of variation of components of matrix δ. The efficacy of this approach is assessed in the simulations presented in Sect. 7.
Then, the right sides of Eqs. (18)(19)(20) are continuous functions in the relevant variables and the following inequality holds, where x(t, x 0 ) is a solution to either (1) or (2) and z k t, z k,0 is a solution to one of the corresponding scalar Eqs. (18), (19), or (20), respectively.
Proof The proof of this statement directly follows from our previous derivation. In fact, the implied conditions assure that the matrices a k (t), k ≥ 1 are diagonalizable and their eigenvalues and eigenvectors are continuously differentiable [29,30]. This implies that α k (t) and, hence, α k,max (t), and other terms in the right side of (18)(19)(20) are at least continuous functions.
Next, diagonalization of matrices a k (t) by the appropriate Lyapunov transforms leads to (17) and consequently forges the design of auxiliary Eqs. (18) or (19) ∀y k ∈ C n , which prompts inequality (21) .
Apparently, there is no warranty that inequality (21) becomes sharper for every successive value of k which reminisced alike behavior of asymptotic series. Still, to prove the boundedness/stability of solutions to our initial multidimensional equations, one should attest a required property for at least one scalar auxiliary equation with k K .
Our approach is essentially simplified if B(t) is already a slow matrix. In this case, it let to conform and extend some criteria of stability of the trivial solution to Eq. (2) which were developed via the utility of nonlinear versions of so known frozen coefficients approach, see, e.g., [19,20]. Moreover, our methodology also embraces the boundedness criteria for nonhomogeneous systems of this kind which practically were overlooked in the current literature.
Let us note as well that Eqs. (18) and (19) are more general than our prior auxiliary Eqs. (14) and (15) and approach their former counterparts as δ i j → ∞.
Local analysis of the boundedness and stability of solutions to Eqs. (19) and (20) is abridged and can be carried out in nearly closed-form under the following conditions, where functions l k ; [t 0 , ∞) × R ≥0 → R ≥0 are continuous in t and bounded inẑ k . Clearly, (22) can be interpreted, for instance, as an application of Lipschitz continuity condition to a scalar nonnegative function L k (t, z k ). More efficient definition of functions l k t,ẑ k can be readily given if L k (t, z k ) are piecewise polynomial in z k , see, e.g., [27] for further details. The next two sections develop novel boundedness/stability criteria and estimate the relevant regions for Eqs. (1) and (2), respectively, under the presumption that all values of t ≥ t 0 are regular. Our expositions in these two sections are similar to ones made in [27]. To underscore these similarities, we employ alike notations. Yet, the current derivations employ more general auxiliary Eqs. (19) and, thus, deliver essentially more general results.

Local stability and boundedness criteria
This section utilizes (22) to derive simplified stability and boundedness criteria and abridge the estimations of boundedness/stability regions.
Application of (22) to (19) yields a sequence of scalar linear equations which can be written as follows, where functions μ k t,ẑ k α k,max (t) + G k (t) + l k t,ẑ k are continuous in t and bounded inẑ k , and functions F k (t) are continuous. The solutions to (23) take the form, Note that below we frequently adopt a reduced notation, i.e., Z k,h t, t 0 ,ẑ k Z k,h t,ẑ k and Z k,F t, t 0 ,ẑ Z k,F t,ẑ . Clearly, Z k,h t,ẑ k > 0, ∀t ≥ t 0 , ∀ẑ k > 0 and θ k t, τ ,ẑ > 0, ∀t > τ, ∀t ≥ t 0 , ∀ẑ k > 0, and the last inequality yields that Z k,F t,ẑ > 0, ∀t ≥ t 0 , ∀ẑ k > 0.

Stability criteria for time-varying homogeneous nonlinear systems
Let us write a homogeneous counterpart to (23) with the following solutions, Z k t, z k,0 ,ẑ k z k.0 Z k,h t,ẑ k for further referencing. As is known, the stability/asymptotic stability of Eq. (25) can be described through conditions set on the fundamental solution matrices for these equations, see, e.g., [31] and additional references therein. Application of these conditions yields that Eq. (25) with k K ≥ 1 is stable if and only if, and is asymptotically stable if and only if whereẑ K ,b andẑ K ,B are the maximal values ofẑ K for which either (26) or (27) hold.
Next lemma provides sufficient conditions prompting the linearization of (20) by (22).

Lemma Assume that ∃k
K ≥ 1 such that a function μ K t,ẑ K is continuous in t and bounded inẑ k ∀t ≥ t 0 , ∀ẑ k > 0 and (25) is stable ∀ẑ K ∈ 0,ẑ K ,max , whereẑ K ,max is the maximal value ofẑ K for which (25) is stable that can be infinity. Also assume that z K ,0 <ẑ (25) is stable, which yields that sup t≥t 0 Z K t, t 0 , z K ,0 ,ẑ K Z K ,h,s t 0 ,ẑ K z K ,0 < ∞ and, consequently, (28) holds .
Note that condition z K ,0 <ẑ K /Z K ,h,s t 0 ,ẑ K is automatically assured if z K ,0 is sufficiently small,ẑ K > 0 and (25) is stable.
Theorem 2 Assume that ∃k K ≥ 1 such that assumptions of Theorem 1 and Lemma are met, and Eq. (25) is stable/ asymptotically stable ∀ẑ K ∈ 0,ẑ K ,max . Then, the trivial solution to (2) is stable/ asymptotically stable, respectively, and in both cases inequality (21) takes the form, where Z K t, z K ,0 ,ẑ K , z K t, z K ,0 and x(t, x 0 ) are solutions to Eqs. (25), (20) and (2), respectively. Furthermore, the region of stability or asymptotic stability of the trivial solution to (2) includes the following set of initial vectors, Clearly, (30) defines the interior of the relevant ellipsoid.
Proof Firstly, let us show that (25) is stable. Pretend in contrary that t * > t 0 is the smallest value of t such that z K t * , z K ,0 ẑ K under prior conditions. Then, (22) and, thus, Eq. (25) holds for ∀t ∈ [t 0 , t * ] which, due to the comparison principle and (28), implies that z K t * , z K ,0 ≤ Z t * , z K ,0 ,ẑ K <ẑ K . Thus,z K t, z K ,0 <ẑ K , ∀t ≥ t 0 , ∀z K ,0 <ẑ K /Z K ,s which enables the application of the linear Eq. (25). Then, the application of comparison principal prompts that z K t, z K ,0 ≤ Z t, z K ,0 ,ẑ K , ∀t ≥ t 0 , z K ,0 <ẑ K /Z K ,h,s ,ẑ K ∈ 0,ẑ K ,max which implies that the trivial solution to (20) is stable.
Due to the last inequality, (21) can be written as (29) which implies stability of the trivial solution to (2) under the conditions of this theorem.
Let us subsequently estimate the regions of stability of the trivial solutions to Eqs. (20) and, in turn, (2). Clearly, a solution to (20) is stable if z K ,0 < supˆz K ∈(0,ẑ K ,max) ẑ K /Z K ,s t 0 ,ẑ K which permits linearization of (20). This implies that the region of stability of the trivial solution to (2) is defined by (30) In the sequel, we assume that (25) is asymptotically stable and other conditions of this statement hold. Then, under this more conservative condition, our prior inferences and (29) hold, which implies that lim t→∞ Z K t, z K ,0 ,ẑ K 0 and, subsequently,lim t→∞ x(t, x 0 ) 0 . The above statement prompts some apparent sufficient stability conditions of the trivial solution to (2) which can be readily checked.
Proof Really, under the above conditions, (25) is either stable or asymptotically stable, which, due to Theorem 2, assures this statement .
Proof In fact, the condition of this statement assumes that the Lyapunov exponent of solutions to Eq. (25) with k K is negative which yields that the trivial solution to this equation is asymptotically stable [2] and, due to Theorem 2, (29) and (30) hold .

Boundedness criteria for solutions of nonhomogeneous time-varying nonlinear systems
Now we are going to develop some boundedness criteria of solutions to (1) and estimate the relevant regions of the initial data. These conditions frequently assure either the stability or asymptotic stability of the forced solutions to (1) and, thus, are important in applications. Our inferences have some correspondence with the concept of inputto-state stability developed by Sontag [32,33] under more restrictive conditions. A solution to (23) Z K t, t 0 , z K ,0 ,ẑ K < ∞, ∀t ≥ t 0 ,ẑ K ∈ 0,ẑ K ,1 if the trivial solution to (25) is either stable or asymptotically stable and.
whereẑ K ,F is the maximal value ofẑ K for which (32) holds, andẑ K ,1 min ẑ K ,max ,ẑ K ,F . This leads us to.

Theorem 3 Assume that ∃k
K ≥ 1 such that assumptions of Theorem 1 hold, Eq. (25) is stable ∀ẑ K ∈ 0,ẑ K ,max , function μ K t,ẑ K is continuous in t and bounded inẑ K , and (32) and inequality hold. Then, (I) x(t, x 0 ) ≤ ∞, ∀t ≥ t 0 , where x(t, x 0 ) are solutions to (1) stemming from the interior of the ellipsoid in R n which is defined as follows, (II) inequality (21) takes the form, where Z t, z K ,0 ,ẑ K and x(t, x 0 ) are solutions to (23) and (1), respectively, and (III) if ∃k K such that (25) is asymptotically stable, other conditions of this statement hold, and lim t→∞ sup V K (t) V K < ∞; then, Proof Firstly, we show that under the condition of this theorem, sup t≥t 0 As is mentioned prior, Z K ,h,s t 0 ,ẑ K < ∞ since we assume that for k K (25) is stable. Then, (25) is stable and z K ,0 < z K ,F . In fact, pretend, in the contrary, that t * > t 0 is the smallest value of t such that z K t * , z K ,0 ẑ K . Then, (22) and, thus, Eq. (23) holds for ∀t ∈ [t 0 , t * ]. Hence, the comparison principle and (37) yield that z K t * , z K ,0 ≤ Z t * , z K ,0 ,ẑ K <ẑ K . This contradiction authenticates that z K t, z k,0 <ẑ K , ∀t ≥ t 0 , ∀z K ,0 < z K ,F ,ẑ K ∈ 0,ẑ K ,1 which enables the application of Eq. (23) in our inferences.
Consequently, application of the comparison principal returns that Clearly, the last inequality let us rewrite (21) as (35), which shows that the norms of solutions to (1) are bounded if they are outgoing from the region specified by (34).
Next, two corollaries present simplified sufficient conditions under which the presumptions of Theorem 4 are met.

Corollary 3 Assume that ∃k
K such that conditions of Theorem 1are met, the inequalities (31), (32), and (33) hold, and function μ K t,ẑ K is continuous in t and bounded inẑ K . Then, solutions to (1) are bounded in the region that is defined by (34) and inequalities (35) and (36) hold ifẑ K ,1 min ẑ K ,v ,ẑ K ,F . Proof In fact, (31) implies asymptotic stability of (25) for ∀ẑ K ∈ 0,ẑ K ,v which, together with other assumptions, assures that the conditions of theorem 3 hold .
In the sequel, let us assume that the Lyapunov exponent of solutions to (25), i.e., φ K ẑ K < 0,ẑ K ∈ 0,ẑ K ,φ , see Corollary 2 for more details. This This leads to.

Proof In fact, in this case
Note that the above statements grant boundedness/stability conditions of solutions to Eqs. (1) or (2) in abridged form.

Nonlocal boundedness/stability criteria
Under the assumption that that all values of t ≥ t 0 are regular, this section devises the boundedness or stability criteria of the initial Eqs. (1) and (2) 19) and (14), foster a formal similarity between our current and former criteria that were developed in [27]. Yet, the current criteria are substantially broader since they rely on the applications of the more general auxiliary equations.
The next theorem accentuates that under some natural conditions, the solutions to (19) monotonically increase in the initial value for ∀t ≥ t 0 which abridges simulations of this equation.

Theorem 4 Assume that Eqs. (19)/ (20) admit unique solutions ∀z
Also, assume that z k t, z k0 and z k t, z k0 are solutions to (19)/ (20) Proof In fact, this theorem directly follows from the uniqueness property of solutions to the scalar Eqs. (19)/ (20) .
The prior statement leads to the boundedness/stability criteria for our initial equations which can be tested via reduced numerical simulations of (19)/ (20). To formulate these criteria, let us first define a set of centered at zero concentric ellipsoids, E(z) ⊂ R n as follows, Also, we assume that ∂ E(z) ⊂ R n−1 defines the boundaries of these ellipsoids and . This prompts the following.

Theorem 5 1. Assume that ∃k
K ≥ 1 such that the trivial solution to the relevant Eq. (20) is either stable or asymptotically stable. Then the trivial solution to (2) is stable/asymptotically stable as well. Furthermore, let the interval [0, z K ) defines the region of stability/asymptotic stability of the trivial solution to (20). Then, the set E − (z K ) is enclosed in the stability region of (1).
Proof Really, the proof of both statements follows from inequality (21), where it is presumed that x(t, x 0 ) and z(t, z K 0 ) are solutions to either Eqs. (2) and (20) or (1) and (19), respectively .
Thus, to estimate the trapping/stability regions of multidimensional Eqs. (1) or (2) about x ≡ 0 we, firstly, shall distinguish in simulations the threshold value z K which defines the relevant intervals of the initial values for Eqs. (19) or (20). This task is markedly streamlined since z K (t, z K 0 ) is monotonically increases in z K 0 , ∀t ≥ t 0 due to theorem 4. Then the ellipsoids estimating the trapping/stability regions for Eqs. (1) or (2) are defined by the following formula, z K V −1 (t 0 )x 0 . To get further insight into the behavior of solutions to our nonlinear auxiliary equations, we can bound from the above and below the right sides of (19) or (20) by their time-invariant counterparts which introduce two scalars, autonomous and, thus, integrable equations. Solutions to these equations define bilateral bounds on the relevant solutions of (19) or (20) which evoke an analytical assessment of the behavior of the solutions to our initial equations. This approach was realized in [27] using more restrained auxiliary equations but can be immediately extended to a wider class of systems corresponding to the auxiliary equations developed in this paper.

Auxiliary equation in a neighborhood of an irregular time moment
Currently, let us presume that a k (t), k ≥ 0 are matrices of general position which are not diagonalizable for some isolated irregular values of t where their normal forms possess two-dimensional Jordan blocks corresponding to repeated eigenvalues, see pp. 219-235 [28].
To simplify further reasoning, let us pretend that matrices a k (t), k k 1 possess simple eigenvalues ∀t ≥ t 0 but a matrix a k 1 (t) is not diagonalizable at t t k 1 and is a small interval about this isolated irregular time moment. Thus, our prior technique fails to provide the auxiliary equation at t t k 1 ,∀k > k 1 . Nonetheless, under some additional conditions, we can define the auxiliary equation in a neighborhood of an isolated irregular time moment using the mean value theorem for matrices which we adopt in the following form where a small real value d 0 and a k 1 (t) is the derivative of the matrix a k 1 (t).
Next, ∀k > k 1 and ∀t ∈ k 1 , we set that V k 1 k (t, y k , d), ∀y k ∈ C n , ∀t ∈ k 1 and also recall that the Lyapunov transforms introduced in Sect. 3, i.e.,y k−1 (t) v k (t)y k (t) remain intact in the current section ∀t / ∈ k 1 and k ≥ 1 which leads to the following, Theorem 6 Suppose that the matrix a k 1 t k 1 , t k 1 > t 0 is not diagonalizable but all other assumptions of Theorem 1 are met and the matrix a 0 (t) is k+1-times continuously differentiable ∀t ∈ k 1 . Also assume that matrices a k 1 (t + d) and a k 1 ,k (t + d) , k > k 1 are matrices of general position which possess simple eigenvalues ∀t ∈ k 1 and some small d 0 and that functions L k 1 ,k (t, z k , d), ∀t ∈ k 1 , z k ≥ 0, k > k 1 are continuous in t, z k and d.
Next, using (39), we write the underlying equation for k k 1 and ∀t ∈ k 1 as follows, Then, subsequent applications to the last equation of the transforms y k (t + d) v k+1 (t + d)y k+1 (t + d),∀t ∈ k 1 , ∀k ≥ k 1 lead to the following equations, In turn, we alter the previous equation as follows Consequently, the application of the technique used in Sect. 3 to the last equation, brings the second equation in (40) which assures the application of inequality (41)∀t ≤ t k 1 + t k 1 .
In turn, the inequality x 0 holds due to the comparison principle, which assures (41)∀t ≥ t 0 .
Note that immediate extension of the prior reasoning assures Theorem 6 under more relax conditions that some matrices a k 1 (t), …,a k J (t),k 1 < ... < k J are not diagonalizable at the same isolated time moment and also if some of matrices a k (t) are not diagonalizable at a finite number of isolated time moments. However, such generalizations of the prior theorem are left beyond the scope of this paper.
Consequently, we extend the application of Theorem 2 to the case concerned in this section. For this sake, we restrain inequality (22) on the pertaining time interval as follows, where a function l ↔k 1 ,k t,ẑ k , d is continuous in t and bounded inẑ k and d. Clearly, we employed in (22) and (42) the same value ofẑ k for simplicity which, in some cases, might extend the conservatism of the follow-up inferences. The applications of (22) and (42) yield a set of three linear equations with matching initial conditionṡ where functions μ ←k , μ →k , F←k and F→k are developed by restraining of their relevant counterparts, that are defined in Sect. 4, on the appropriate time intervals and Z→k(t), ∀t >t k 1 + t k 1 , ∀k > k 1 due to the comparison principle.
Next, let us define piecewise continuous functions Then, alike (24), a continuous general solution to (43) can be written as follows . Assume subsequently that ∃k K ≥ k 1 such that a homogeneous counterpart to (43) is stable ∀ẑ K ∈ 0,ẑ K ,Max . Then, which let us reframe Theorem 2 as follows.
Theorem 7 Assume that the assumptions of Theorem 6 are met and ∃k K > k 1 such that function l K t,ẑ K is continuous in t and bounded and inẑ K , function l ↔k 1 ,k (t,ẑ k , d) is continuous in t and bounded inẑ K and d, and the homogeneous counterpart of (43) is stable or asymptotically stable ∀ẑ K ∈ 0,ẑ K ,Max . Then, the trivial solution to Eq. (2) is stable or asymptotically stable, respectively, and in both cases, inequality (29) takes the form, (II) inequality (21) takes the form, where Z K t, z K ,0 ,ẑ K , d and x(t, x 0 ) are solutions to (43) and (1), respectively, and.
(III) if ∃k K > k 1 such that homogeneous counterpart to (43) is asymptotically stable, other conditions of this statement hold, and lim t→∞ sup V K (t) Proof. The proof of this statement can be attained directly by the utility of (44) and the apparent adjustment of the proof of Theorem 3 .

Simulations
This section applies the developed methodology to evaluate the degree of stability/ boundedness and estimate the evolutions of the norms of solutions to two nonlinear systems with time-dependent nonperiodic coefficients which are frequently appeared in various applications [34]. These systems comprise two coupled Van der Pol-like or Duffing-like oscillators with variable coefficients. The simulations authenticate our theoretical inferences and contrast the efficacy of the current and prior techniques [27] in their common application domains.
Note that the technique developed in [27] works under the hypothesis that A is a stable and X(t) is zero-mean matrices, whereas our current approach does not rest on these assumptions.
To simplify our task, we adopt in this section the appropriate formulas from Sect. 3 where we set that k 1.
In turn, let us recall that v 1 (t) v 1,i j (t) , i, j 1, ..., 4 and 1 (t) diag λ 1 (t), λ * 1 (t), λ 2 (t), λ * 2 (t) are the eigenvectors and eigenvalues matrices of a 0 (t), λ k (t) and Finally, we select G 1 G 1 − Im(diag(G 1 )) and set that where formulas for e i (t), i 1, 2 were derived in [27] under the assumption that v 1 const; still these derivations remain intact if v 1 v 1 (t). Below we briefly replicate some of these steps for the sake of completeness.
Let  3 which, finally, returns (55). Afterward, let us assess in simulations the efficacy of our optimal criterion which yields separation of temporal scales of matrix B(t) via maximizing the degrees of stability of the linearized auxiliary Eq. (20) with k 1. This also let us to contrast the methodologies developed in [27] and in the current paper since the former -relays on averaging of B(t) on infinite interval t 0 ∞ .
To reduce the scope of simulations, we set that δ i j δ, i, j 1, 2 which embraces a suboptimal separation of temporal scales in B(t). Figure 1a-f examines the dependence of the Lyapunov exponents φ 11 and φ 12 upon δ, see Sect. 3 of this paper. For this sake, we used the following sets of parameters: in Fig. 1a, a 11 a 12  by the values of the core frequencies of the variable coefficients used in simulations. The curves in Fig. 1a and d as well as 1b and 1e are fairly similar. Figure 1a and d displays that the current approach is somewhat superior to the one developed in [27] which averages B(t) over interval [t 0 , ∞). Figure 1c shows that the technique developed in [27], which averages B(t) on infinite interval [t 0 , ∞), is practically optimal for the set of parameters used in these simulations. Yet, Fig. 1b, e, and f displays that the current methodology frequently selects the optimal moving average with a window of finite length, which is determined by the core frequencies of the variable coefficients included in the relevant simulations. Figure 2a-d displays in solid and dashed lines time histories of slow, i.e., α 1,max t, δ , and fast, G 1 t, δ , time-varying components of the linear part of Eq. (20) which are simulated for the optimal value of δ and the sets of parameters used in Fig. 1a-c and f, respectively. Clearly, the proposed criterion provides a quite effective separation of the temporal scales of the linear block of the auxiliary equation despite that we evoke here the most restricted and, hence, suboptimal version of this approach.
Note that all subsequent simulations of the auxiliary equations presented in this section include suboptimal values of δ that were preliminarily defined for the relevant sets of parameters.  20) and (18) with F 1 0, respectively. These simulations assumed the following set of parameters, ε 1 ε 2 0.5, a 11 a 12 a 21 a 22 0.5, ω 11 20, ω 12 1, ω 21 10, ω 22 0.5. As is expected, solutions to Eq. (20) provide a more efficient estimation of the norms of solutions to Eq. (2) than its more conservative counterpart. Figure 3b displays in solid and dashed lines the time histories of the norms of the solutions to Eqs. (2) and (20) for the following set of parameters, ε 1 ε 2 0.5, a 11 a 12 0, a 21 a 22 0.1 ,ω 21 10, ω 22 0.5. Clearly, the reduction in the amplitudes of variable coefficients improves the accuracy of our estimates. Figure 3c and d contrasts in the solid and dashed lines the time histories of the norms of solutions to the nonhomogeneous Eq. (1) and its scalar counterpart (19), where we set that F 2 0 and either F 1 0.025 or F 1 0.05, respectively. Note that the other parameters in these two figures and Fig. 3b are the same. Clearly, the precisions of our estimates are correlated inversely with F .  (19), respectively. In this case, we set that F 1 0.025, F 2 0. As prior, the latter simulations are used to estimate the threshold value of z 1 t, z 1,0 , i.e., z 1 . Clearly, the developed technique provides rather conservatively estimates of the boundary of the actual trapping region which, nonetheless, authenticate our theoretical inferences. Note that additional simulations provide alike estimates for projections on complementary coordinate frames.

Coupled Duffing-like system
In a coupled Duffing-like system, the linear components are still defined by (53), whereas the definition of nonlinear components should be adjusted as follows, which prompts the relevant correction in (55), i.e., v 2k → v 1k and v 4k → v 3k . respectively. As for the Van der Pol-like system, the solutions to Eq. (20) provide more efficient estimates than the ones defined by (18) with F 1 0. These simulations assumed the following set of parameters a i j 0.5, i, j 1, 2, ω 11 0.94, ω 12 2.12, ω 21 20,ω 22 0.3, ε 1 ε 2 0.5. Figure 4b displays in solid and dashed lines time histories of the norms of the solutions to Eq. (2) and its upper estimates which are developed in simulations of (20) with matched initial values for the following set of parameters, ε 1 ε 2 0.5, a 11 a 12 0, a 21 a 22 0.1 , ω 21 10, ω 22 0.5. As in Fig. 3b, the amplitudes of the variable coefficients are scaled down in this figure. Nonetheless, the accuracy of our estimates improves marginally in this case. Figure 4c contrasts in dashed and solid lines projections on x 1 × x 2 -plane of the boundaries of the stability region of the trivial solution to (2) which were developed in simulations of this equation and its scalar counterpart (20), respectively. Note that in Fig. 4c-f, we set a i j 0.1, i, j 1, 2 whereas other parameters adopted in these simulations are the same that the ones used in Fig. 4a.
Furthermore, Fig. 4d and e contrasts in solid and dashed line time histories of the norms of solutions to (1) and its scalar counterpart (19). In these figures, we set that F 1 0.01, F 2 0 and F 1 0.025, F 2 0, respectively. As is expected, enlargement of F decreases the accuracy of our estimates. Figure 4f contrasts in dashed and solid lines projections on x 1 × x 2 -plane of the boundaries of the trapping region of (1) which are developed in simulations of this equation and its scalar counterpart (19), respectively. In this case F 1 0.01, F 2 0.
Apparently, the developed estimates of the norms of solutions to the coupled Van der Pol-like and Duffing-like systems provide adequate accuracy if the solutions emanated from a central part of the stability/ trapping regions. Yet, the boundaries of these regions are estimated rather conservatively by our current methodology. However, our estimates become more plausible if the original systems include some unknown but bounded in norm components -a standard presumption in control theory. Still, all the above simulations affirmed our theoretical inferences.

Conclusions and future studies
This paper develops a novel methodology for estimation of the degree of boundedness/stability of nonlinear systems with time-dependent nonperiodic coefficients, which builds upon our recent studies in this area [25][26][27]. In [25], we develop a scalar auxiliary equation with solutions bounding from the above the norm of solutions to multidimensional nonlinear equations with time-varying coefficients, which let us to enhance the boundedness/stability criteria and estimate the corresponding regions for the underlying systems. However, some underlying assumptions constrain the broad application of this methodology, which, consequently, was recast in [27] under the condition that the average of time-dependent matrix of the linear block of our system is a stable matrix.
In contrast, the current approach is substantially more flexible. It merely requires splitting the time-dependent matrix of the linear block into slow and fast time-varying components via the application of moving averages and the optimal criterion developed in this paper. Next, the subsystem with a slow matrix is simplified using the relevant Lyapunov transform. Repeated applications of such transformations yield a sequence of abridged systems which leads to the development of their scalar auxiliary counterparts. Ultimately, the windows of moving averages are defined by maximizing the degree of stability of the linearized homogeneous auxiliary equation. This minimizes the conservatism of our estimates and aids the estimation of the behavior of solutions to the original systems through the assessment of the matching solutions to scalar auxiliary equations. Using this reasoning, we developed novel boundedness/stability criteria and estimated the boundedness/stability regions of the original systems.
Abridged estimates of the solutions to scalar nonlinear auxiliary equations were developed using linear upper bounds for these equations. This lets us simplify the assessment of the evolutions of the norms of solutions to the original systems stemming from their trapping/stability regions and estimate the boundaries of these regions. The inferences developed by the application of this technique should be useful in appropriate control tasks.
The current methodology is simplified if the underlying system possesses only slowvarying coefficients. In this case, it can be used to authenticate and extend practically plausible frozen coefficients or gain scheduling techniques that were frequently used in applications without appropriate justification.
Our inferences were examined in inclusive simulations that were partly discussed in this paper. These simulations contrast our current and prior methodologies and show that solutions to our scalar auxiliary equations deliver sound upper bounds of the norms of solutions to the original systems if they are stemmed from the central parts of their trapping/stability regions. However, the efficacy of our estimates is declined if the solutions originated near the boundaries of these regions, which, in turn, are estimated rather conservatively. Nonetheless, all simulations reinforce our theoretical inferences and turn out to be more plausible for systems under uncertainty. Yet, integration of the methodologies developed in this paper and in [2] should substantially enhance the estimation of the boundaries of the trapping/stability regions as well.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.