CC BY
3DIFFERENTIAL EQUATIONS AND
CONTROL PROCESSES N. 1, 2024 Electronic Journal, reg. N $C77-39410 at 15.04.2010 ISSN 1817-2172
http://diffjournal.spbu.ru/ e-mail: jodiff@mail.ru
Delay differential equations
On the recent progress in effective dimension estimates for delay equations
Mikhail Anikushin* Andrey Romanov**
Department of Applied Cybernetics, Faculty of Mathematics and Mechanics,
St Petersburg University
Email: *demolishka@yandex.ru **romanov.andrey.twai@yandex.ru
Abstract. We discuss the recent progress on studying volume contraction for linear cocycles generated by delay equations obtained by the authors. On the basis, we use adapted metrics constructed explicitly or via the Frequency Theorem. In contrast to many existing results, this approach allows to provide effective estimates in terms of the system parameters (including delays). We illustrate the exposed general results by means of the nonautonomous and classical Nicholson blowflies models, where effective dimension estimates for global attractors and robust conditions for the global stability are obtained.
Keywords: Delay equations, Volume contraction, Dimension estimates, Global stability, Nicholson blowflies
1 Introduction
Estimating dimensions of global attractors is one of basic problems in dissipa-tive dynamics. Although the initial interest in the problem was concerned with finite-dimensional reduction (see S. Zelik [45]), further developments revealed a natural generalization of the Bendixson criterion and its applications to the problem of global stability (see R. Smith [41]; M.Y. Li and J.S. Muldowney [34, 35]; G.A. Leonov and V.A. Boichenko [33]). From this perspective, the problem of effective computation or estimation of dimensions is of high interest. In this direction, the main approach is concerned with the study of volume contraction and the related dimensional-like characteristic called the Lyapunov dimension. In [1], it is shown that the Lyapunov dimension theoretically admits exact computation via adapted metrics on exterior products. Although many authors use standard metrics to provide effective estimates (mainly in the case of ODEs and parabolic equations; see R. Temam [44]), sharper estimates or even exact computations require considerations of adapted metrics (see N.V. Kuznetsov and V. Reitmann [25]). This is also the case for some problems in which standard metrics are not appropriate.
This work is a survey devoted to the problem of providing effective dimension estimates for delay equations. On the geometric level, the problem is concerned with studying volume contraction for linear cocycles generated by linearized delay equations. It was clearly posed in the work of the first author [6] based on the analysis of preceding results.
Among them, it should be highlighted the pioneering work of J. Mallet-Paret [40] which utilizes compactness of the cocycle mappings to establish finiteness of dimensions. This method does not provide any way to obtain effective estimates. Starting from [40], most of the works interested in delay equations and dimension estimates follow similar approaches therefore making only qualitative conclusions on the finiteness of dimensions. This is reflected in the classical monographs (for example, J.K. Hale [20]) as well as in relatively recent ones (see A.N. Car-valho, J.A. Langa and J.C. Robinson [13]; I.D. Chueshov [17]). Moreover, even recent works, for example of W. Hu and T. Caraballo [22], proceed to develop this approach. In [22], the resulting estimates are called explicit, but they are not effective because of appealing to objects associated with exponential dichotomies of non-self-adjoint linear problems which cannot be computed or analyzed effectively. In addition, the utilization of exponential dichotomies for the problem of dimension estimates is artificial and therefore it results in estimates with wrong
asymptotics.
It is well-known that one of general ways to provide effective dimension estimates is concerned with the Liouville trace formula and resulting from it upper estimates such as the Douady-Oesterle or Constantin-Foias-Temam estimate (see N.V. Kuznetsov V. Reitmann [25]; R. Temam [44]; S. Zelik [45]). To the best of our knowledge, the first attempt to apply this approach for delay equations was done by J.W-H. So and J. Wu in [43]. However, the authors did not succeed to derive from it effective estimates. In particular, at the end of the introduction they wrote: "The verification of the hypotheses in Theorem 3.10 for reaction-diffusion systems with delays is a non-trivial task and will be reported in a future paper." Since then, none of such declared results appeared.
A reason for such a failure is explained in [6] by means of a simpler example. Namely, to derive from the Liouville formula some uniform upper estimates, one has to compute or estimate the so-called trace numbers. For the computation, one applies the symmetrization procedure which sometimes allows to describe the trace numbers via eigenvalues of self-adjoint operators obtained via the additive symmetrization* (see [1]). In [6], for particular delay operators, it is shown that the trace numbers computed in the standard metric are irrelevant (in particular, they do not depend on the delay). Moreover, for general delay operators, the sym-metrization procedure itself requires additional justification and the false proof (and the false general statement) of Theorem 7 in [6] only confirms this.
In [1], the general symmetrization problem is stated and its investigation is illustrated by means of a fairly general class of delay operators. In particular, it is shown that the resolution of the problem significantly depends on the choice of a metric from the considered class. Moreover, the study allowed to obtain for the first time effective dimension estimates for global attractors in the Mackey-Glass equations and the periodically forced Suarez-Schopf delayed oscillator [5] and the estimates are numerically justified to be asymptotically sharp as the delay tends to infinity. In Section 3, we will discuss some of these results in more details.
As to [43], the authors exclude discrete delays from consideration due to their inability to construct semigroups in Hilbert spaces (this problem was resolved in [6]) and consider only distributed delays. Note that such delays may cause the presence of a continuous spectrum in the additive symmetrization and how to
*For parabolic equations, the symmetrization procedure is trivial and it is usually used implicitly. In the case of reaction-diffusion systems, upper estimates for the trace numbers in standard metrics lead to relevant asymptotics of dimensions [44]. To provide sharper estimates or to study other parabolic equations such as the 2D Navier-Stokes equations, one has to develop additional techniques, for example, the theory of Lieb-Thirring inequalities (see V.V. Chepyzhov and A.A. Ilyin [14]).
deal with this is still an open problem (see [1]).
Another approach to provide effective dimension estimates for delay equations is developed by the first author in [2]. It is based on applications of the recently obtained version of the Frequency Theorem [4] (see also [8]) to study m-fold compound cocycles acting on exterior powers of the Hilbert space. For this, it is required to study their infinitesimal generators along with structural and regularity properties of the associated linear inhomogeneous problems and most of the study [2] is devoted to this. As a result, we obtain frequency conditions for the exponential stability of m-fold compound cocycles. Geometrically, such conditions imply the existence of a constant adapted metric (on the exterior power) which allows to establish the exponential stability.
In [3], the authors developed an approximation scheme to verify some of the arising (in the above study) frequency inequalities for the case of scalar equations*. For m = 2, such conditions (called the Frequency Criterion) are expected to imply global stability in the system due to the generalized Bendixson criterion of M.Y. Li and J.S. Muldowney [35] and robustness of the Frequency Criterion. We will discuss this in Section 4 and illustrate the criterion by means of the classical Nicholson blowflies model in Section 6.
Among preceding works it should be also mentioned the study of J. Mallet-Paret and R.D. Nussbaum [38] concerned with cocycles on injective tensor products generated by a class of scalar delay equations arising after linearization of scalar nonlinear equations with monotone feedback. Such nonlinear equations are known to satisfy the Poincare-Bendixson trichotomy (see J. Mallet-Paret and G.R. Sell [39]). In [38], it is established a comparison principle which allows to compare the Floquet multipliers over periodic orbits and stationary points (considered as periodic with the same period). In [6], this principle is combined with the Ergodic Variational Principle (see [1]) to derive effective dimension estimates*. However, systems of equations or scalar equations without monotone feedback go beyond this context.
Moreover, effective dimension estimates for delay equations can be obtained in the case of small delays by studying vector fields in Rn. This is possible due to the existence of inertial manifolds for such equations (see [4]). In [15, 16], C. Chicone obtained a power series expansion of the vector field on the inertial manifold in terms of the small delay. This series can be used to symmetrize the vector field and provide effective dimension estimates.
^The general scheme works for systems of equations.
*In [5] it is shown that the Frequency Criterion may improve these estimates.
There are results on fractal dimensions of forced almost periodic oscillations obtained by the first author [10] and the authors joint with V. Reitmann [11] for almost periodic ODEs which are interesting to generalize for delay equations. Here a problem arises due to the fact that quadratic functionals constructed by the Frequency Theorem [4] may be not coercive.
This paper is organized as follows. In Section 2, we discuss well-posedness and linearization for a class of nonautonomous nonlinear delay equations in rn posed in a proper Hilbert space h. In Section 3, we consider effective dimension estimates via additive symmetrization of delay operators. In Section 4, we state a frequency criterion guaranteeing volume contraction for m-fold compound cocycles and expose a numerical scheme for its verification in the case m = 2. In Section 5, we obtain effective dimension estimates for global attractors in the nonautonomous Nicholson blowflies model. In Section 6, we apply the Frequency Criterion to provide robust conditions for the global stability in the classical Nicholson blowflies model and compare the result with other methods.
2 Well-posedness and linearization of delay equations in Hilbert spaces
In this section, we are going to discuss the well-posedness and linearization of the following nonautonomous nonlinear delay equations over a semiflow (Q,$) on a complete metric space Q. Over a given q e Q, the system is described by
x(t) = Ax* + BF (^ (q),Cxt) + W (^(q)), (2.1)
where t > 0 is a fixed value (delay), x(-): [—t,T] ^ rn for some T > 0 and xt(0) := x(t + 0) for 0 e [—t, 0] denotes the T-history segment of x(-) at t e [0, T]; A: C([-t, 0]; rn) ^ rn and C: C([-t, 0]; rn) ^ rr2 are bounded linear operators; B is an n x ri-matrix; W: Q ^ rn is a bounded continuous function and F: Q x rr2 ^ rri is a C 1-differentiable in the second argument continuous mapping such that for some A > 0 we have
|F(q,yi) - F(q,y2)|u < A|yi - ^2|m for any q e Q and yi,y2 e m, (2.2)
where u = rri and m = rr2 are endowed with some (not necessarily standard Euclidean) inner products (this will be used in Section 4).
We need to consider (2.1) as an evolutionary equation in the Hilbert space h = rn x L2(—t, 0; rn). In this space, there is an unbounded operator A associated
with A from (2.1) as
(x,0) A (a0,^0J, (2.3)
where (x, 0) belongs to the domain D(A) of A given by
D(A) := {(x,0) e h | 0(0 e W^2(—t,0;rn) and 0(0) = x}, (2.4)
and denotes the derivative in the Sobolev space W 1'2( —t, 0; rn).
We embed the space e = C([—t, 0]; rn) into h via the mapping 0 a (0(0), 0) for 0 e e. It will be convenient to identify the elements of e and their images in h under the embedding and keep the same notations for the induced operators. For example, we put C(x, 0) := C0 for (x, 0) e e for the operator induced by C from (2.1).
With B from (2.1) we associate a bounded linear operator B: u a h given by Bn := (Bn, 0) for n e u. Moreover, let W: Q a h be associated with W from (2.1) as W(q) := (W(q),0) for q e Q. From this we can treat (2.1) as an abstract nonautonomous evolution equation in h over the semiflow (Q,$) which is described by§
v(t) = Av(t) + BF(^(q), Cv(t)) + W(^(q)). (2.5)
Note that A is the generator of a Co-semigroup G in h. By Theorem 1 in [6], for any v0 e h and T > 0 there exists a unique generalized solution v(-) = v(-; q,v0) = (x(-),0(-)) which is a continuous h-valued function on [0,T] satisfying v(0) = v0, 0(t) = xt in L2(—t, 0; rn) and
v(t) = G(t)vo + i G(t — s)[BF(^s(q),Cxs) + W(^s(q))] ds (2.6) 0
for any t e [0,T]. For understanding (2.6), it is important to note that we may interpret the function [0,T] 9 s a Cxs e m for x(-) e L2(—t, T; rn) as a well-defined element of L2(0,T; m) (see [2, 4, 6]). From the variation of constants formula (2.6), one can show that the mappings
0t(q, v0) := v(t; q, v0) for t > 0, q e Q and v0 e h (2.7)
determine a nonlinear cocycle 0 in h over the semiflow (Q,$) (see [1, 2, 3] for precise definitions).
§From what has been said, we have that (2.5) is an abstract form of the transfer equation Jt^>(¿,0) = da0), where t > 0 and 0 e (—t, 0), and (2.1) describes a nonlinear nonlocal boundary condition of Neumann type at 0 = 0.
One can relate the generalized solutions to the classical solutions in e = C([—t, 0]; rn) as follows. Remind that e is identified with its image in h under the embedding 0 ^ (0(0), 0) for 0 e e. By Theorem 1 in [6], the restriction of ^ to e is the cocycle in e generated by classical solutions (see J.K. Hale [20]). More precisely, we have ^(q,v0) = (x(0),xt) for t > 0, where x(-): [—t, ^ rn is continuous and C^differentiate on [0, to) function satisfying (2.1) for t > 0 and x(6) = 0o(6) for 0 e [—t, 0], where v0 = (0(0), 0o). Moreover, (q, h) c e and for any T > t there exists a constant LT > 0 such that the smoothing estimate
||^(q,vi) — #(q,V2)||E < Lt|vi — V2|h for t e [t,T] (2.8)
is valid for any vi,v2 e h and q e Q.
With such we can associate a skew-product semiflow n in Q x h as
nt(q, v) = (^(q), ^(q, v)) for t > 0, q e Q and v e h. (2.9)
From the smoothing property we have that any strictly invariant subset P c Q x h of n, i.e. such that nt(P) = P for any t > 0, must satisfy P c Q x e.
We may linearize (2.5) in the fiber q e Q along the trajectory of v0 e e (for simplicity), as
¿(t) = A£ (t) + BF /(^t(q ),C^(q,v>))C£(t), (2.10)
where F'(q,y) is the derivative of F(q,y) w.r.t. y.
Similarly to the above considerations, we have that for any p = (q, v0) e Qxe and ¿0 e h there exists a unique generalized solution ¿(t) = £(t; q, v0, ¿0) of (2.10) with £(0) = ¿0 and the mappings
^(p,£0) := £(t; q, V0, ¿0) for t > 0, p = (q,V0) e Q x e and ¿0 e h (2.11)
determine a uniformly continuous linear cocycle (in the terminology of [1, 2]) S in h over the skew-product semiflow (Q x e,n). Analogously to (2.8), the cocycle is smoothing in the sense that ST(p, h) c e and for any T > t there is a constant LT > 0 such that
|St(p,£)||e < Lt|£|h for any t e [t,T], p e Q x e and £ e h. (2.12)
Analogously to Theorem 2 in [6], one can deduce from (2.6) that for any q e Q, v0 e e, T > 0 and any bounded subset B of h we have
lim |^t(q,v0 + h£) — #(q,v0) — hst(p,£)|h = 0 (2 13)
h-m0 h , '
where p = (q,v0) and the limit is uniform in t e [0,T] and £ e B. In other words, E^p, •) is the Frechet differential of ^¿(q, •): h —> h at v0 e e, which continuously depends on p = (q, v0) e Q x e.
Let A1(S) > A2(E) > ... be the uniform Lyapunov exponents of E (see [1]). Suppose there exists m > 0 such that ^m=1 Ak(E) > 0 (for m = 0 the sum is zero) and ^m^1 Ak(E) < 0. Then the Lyapunov dimension' of E is defined by the Kaplan-Yorke formula
a (E)
dimL E := m + . (2.14)
1 Am+1(E)|
If there is no such m, we put dimL E := to. Moreover, for Am+1(E) = —to we have dimL E := m.
Let P be a strictly invariant w.r.t. n compact subset, i.e. nt(P) = P for any t > 0. From (2.13) one may relate the Lyapunov dimension of E and the Hausdorff or fractal dimensions of the fibers Pq := PR ({q} x e) (see [14, 25, 44, 45]). Namely, an appropriate extension of the main result from [14] (which is shown for the autonomous case, i.e. for Q being a one-point set) would give the estimate
dimF Pq < dimL E for any q e Q (2.15)
for the fractal dimension dimF Pq of Pq. Note that due to the smoothing estimate (2.8), the Hausdorff and fractal dimensions of Pq in the metrics of h and e coincide.
3 Effective dimension estimates via additive symmetriza-tion of delay operators
Let us rewrite (2.10) as ¿(t) = A(nt(p))^(t), where A(p) = A + BF'(q,Cv0)C for p = (q,v0) e Q x e. It is clear that A(p) has the same form as A and, in particular, it is naturally defined on the same domain D(A) being the generator of a C0-semigroup.
In [1], it is shown that the operator (A*(p)+A(p))/2, where the adjoint A*(p) of A(p) is computed in appropriate metrics in h, may admit a densely-defined self-adjoint extension (called an additive symmetrization of A(p)) being a bounded". operator in h. Under additional conditions (see below), it is also established that
^More precisely, the global or uniform Lyapunov dimension. 'Thus an additive symmetrization of A(p) is unique in this case.
the additive symmetrization is proper in terms of [1]. This allows to characterize trace numbers of A(p) via eigenvalues of the additive symmetrization. From this, the Liouville trace formula can be applied to derive effective dimension estimates. We are going to describe this result.
As a model for A(p), let us consider an operator L: C([—t, 0]; rn) — rn and the corresponding to it operator L as in (2.3), i.e.
(x,0) —— (l^,^), (3.1)
where (x, 0) e D(L) = {(x, 0) e h | 0(0 e W 1'2(-t, 0; rn) and 0(0) = x}. We suppose that L is given by
j
L0 = L00(0) + l—t 0(—t ) + ^ L—Tj 0(—Tj) (3.2)
j=1
for any 0 e C([—t, 0]; rn). Here L0, L—T and L—Tj are n x n-matrices with Tj e (0,t) for j e {1,..., J} being distinct values.
It is convenient to put t0 := 0 and tj+1 := t. Let p(-) be a positive function on [—t, 0] such that** p(-) e C 1([—Tj+1, —t,]; r) for any j e {0,..., J} and infp(0) > 0. With such p, by fixing an inner product (-, •)Rn in rn, we endow h with the inner product
((x,0), (y,0))p := (x,y)rn + ^ p(0)(0(0),0(0))rn(3.3)
where (x, 0), (y, 0) e h.
It is important that p(^) may have discontinuities at —Tj with j e {1,..., J}. Let Aj(p) := p(—Tj + 0) — p(—Tj — 0) be the jump of p at —Tj. Then, by Theorem 5.2 in [1], A(p) admits an additive symmetrization if and only if Aj(p) = 0 for any j e {1,..., J}. Moreover, to admit a proper additive symmetrization we must satisfy Aj(p) > 0 for any j e {1,..., J}.
A useful example of p satisfying the above properties is
p(0) := eKj0, if 0 e (—Tj+1, —Tj] and j e {0,... , J} (3.4)
with k0 < k1 < ... < kj.
Now consider a closed positively invariant w.r.t. n (from (2.9)) subset P C Q x e. Let V: P — r be a bounded scalar function which satisfies
"This means that p(-) g C1 ((—Tj+i, —Tj); r) and it can be naturally extended to [—Tj+i, —Tj] so it becomes C^differentiate on the closed interval.
(V1) For any T > 0, the mapping [0,T] 9 t ^ V(#*(q)) e r is absolutely continuous;
(V2) For each p e P, there exists the right derivative of V along the trajectory of p, i.e.
V(p) := lim V(n'<p» — V<p>. (3.5)
t
With the aid of p as in (3.4) and V as above the following theorem is established in [1] (see Theorem 5.4 therein).
Theorem 1. In the above context, suppose that for any p e P the operators A(p) have the form as L from (3.1) and (3.2) with the corresponding (n x n)-matrices L0 = L0(p), L_T = L-T (p) and L—Tj = L—Tj (p) for j e {1,..., J}. Let k0 < k1 < ... < kj be given and let Ai(p) > ... > An(p) be the eigenvalues (counting multiplicities) of the symmetric matrix (here T denotes the transpose of a matrix^ and In is the (n x n)-identity matrix)
Lc(p) + LT(p) + eKJTL—t(p)L—T(p) + £ ^k?. + /»■ (3.6)
j=i
Let 0 < K(p) < n be the largest number k such that Ak(p) > —k0. For m > 1 consider the value
/ 1 min{m,K (p)} k \
a+(m) := sup I !/(p) + - Y^ Ak(p) - k • max{0, m - K(p)} I . (3.7)
peP y k=i y
Then the cocycle S given by (2.11) over the semiflow (P, n) satisfies
Ai(S) + ... + Am(S) < a+(m). (3.8)
Here the presence of a Lyapunov-like function V(•) in (3.7) is a development of the Leonov method (see N.V. Kuznetsov [27]; N.V. Kuznetsov et al. [26]; G.A. Leonov and V.A. Boichenko [33]). On the geometric level, it is concerned with variations of a constant metric in h in its conformal class via the Lyapunov function. This method led to exact formulas of the Lyapunov dimension for the Henon and Lorenz-like systems (see G.A. Leonov [29, 32]; G.A. Leonov et al. [28]; G.A. Leonov and T. Mokaev [30]; G.A. Leonov, T.A. Alexeeva and N.V. Kuznetsov [31]). In [37], M. Louzeiro et al. developed an algorithm for
^^More generally, one may take the adjoints w.r.t. any fixed inner product in Rn.
optimization in such a conformal class in the case h = rn and V having a particular structure.
From (3.6) with V = 0 it is clear that if we take k0 > 0, then a+(m) is negative for all sufficiently large m provided that the eigenvalues Ak(p) are bounded from above on P (for example, when P is compact).
In [1], Theorem 1 is applied to provide effective dimension estimates for global attractors in the Mackey-Glass equations and periodically forced Suarez-Schopf delayed oscillator (see [5]). Both models are known to be chaotic for certain parameters and this is the first time when effective dimension estimates for chaotic attractors arising in delay equations were obtained. Moreover, the resulting estimates are asymptotically sharp as the delay tends to infinity that is justified by numerical computations of eigenvalues at equilibria. In Section 5, we will apply Theorem 1 to the nonautonomous Nicholson blowflies model.
It is worth mentioning that the estimate (3.8) depends on the spatio-temporal change of variables x(t) — x(k£) and 0(t, 0) — 0(k£, k0) for some k > 0. Although the uniform Lyapunov exponents also depend on the change (they scale by k), the Lyapunov dimension does not that is clear from (2.14). In applications, this means that we may consider a series of estimates for the Lyapunov dimension derived from Theorem 1 which depend on k and then minimize it over k > 0 to get a better result (see Corollary 3 and below for an example).
4 Frequency criterion for the global stability
In this section, we will describe another approach to obtain dimension estimates for invariant sets arising in delay equations. It is developed by the first author in [2] and concerned with applications of the recent version of the Frequency Theorem [4] to study m-fold extensions Sm (called compound cocycles) of the linearization cocycle S (given by (2.11)) to the exterior power ham of the Hilbert space h. In particular, these results allow to obtain frequency conditions for the exponential stability of compound cocycles. As it is shown in [3] by means of the Suarez-Schopf delayed oscillator and Mackey-Glass equations, at least in the case m = 2 (related to the generalized Bendixson criterion and global stability; see below), this approach improves results which can be obtained by Theorem 1 as well as with the aid of some other existing methods in the field. In Section 6, we will compare our approach with a series of works on the global stability of the classical Nicholson blowflies model.
On the infinitesimal level, the m-fold compound cocycle is given by a perturbation of the m-fold additive compound A[Am] of A and one of possible frequency conditions can be described in terms of the transfer operator W(p) := -cm(A[Am] - pI)-iBm as» (here V0 > 0 and A as in (2.2))
suP IIW(-v0 + < A-i. (4.1)
wGm
Here Cm: em ^ mm and Bm: u^, ^ ham are some operators associated with the measurement operator C and the control operator B from (2.10) respectively; mm and um are Hilbert spaces associated with the spaces m and u from (2.2) respectively; and em is a Banach space such that we have the continuous embed-dings
D(A[Am]) c em C ham (4.2)
We refer to [2] or [3] for precise definitions.
In the case (4.1) is satisfied for m = 2 and v0 > 0 such that —v0 > s(A[A2]), where s (AtA2]) is the spectral bound of A[A2], we obtain
Ai(S) + A2(S) < —v0 and, consequently, dimL S < 2. (4.3)
For autonomous systems (2.1) (i.e. when Q is a one-point set) admitting a global attractor, this condition allows to apply the generalized Bendixson criterion of M.Y. Li and J.S. Muldowney [35] to exclude the existence of closed invariant contours on the global attractor. Moreover, since (4.3) is robust in the sense that it is preserved under small perturbations of the cocycle (see [1]), it is expected that the system must be globally stable, i.e. any trajectory must tend to the set of equilibria. In finite dimensions, one utilizes the Ci-Closing Lemma of C.C. Pugh (see [34, 41]) to obtain this. To the best of our knowledge, there are still no variants of the lemma available in infinite dimensions and, in particular, for delay equations. As an intermediate solution, for some delay equations one may construct inertial manifolds (see [4, 6, 7, 9]) and apply the finite-dimensional version of the Closing Lemma. However, we hope that our investigations will stimulate developments of the problem in infinite dimensions.
Note that it is hard to verify (4.1) purely analytically since it requires to solve a first-order PDE on the m-cube (—t, 0)m with boundary conditions on some k-faces (with k < m) adjacent to {0}m involving diagonal derivatives and delays. Moreover, solutions to such PDEs are not usual smooth functions and have only diagonal derivatives (see [2]).
"For convenience, here we omitted mentioning complexifications of operators.
In [3], we developed an approximation scheme to verify (4.1) in the case of scalar equations and avoiding direct examinations of the arising PDEs. For simplicity, we will describe it for m = 2, the system (2.10) with n = r1 = r2 = 1 and the operator C such that for some T0 e [0,t] we have C0 = 0(—T0) for 0 e C([—t, 0]; r). Thus we consider the scalar equation
x(t) = Ax + BFV(p))x(t — T0), (4.4)
where B and F'(p) can be identified with real numbers. As in (2.2), we suppose that there exists A > 0 such that |F'(p)| < A for all p e P.
Let A be the operator in h = rxL2(—t, 0; r) associated with A from (4.4) via (2.3). Let A1(A) > A2 > ... be its eigenvalues arranged by nonincreasing of real parts and according to their multiplicity. For example, if A0 = a0(0) + b0(—t) for 0 e C([—t, 0]; r), then the eigenvalues are given by the solutions p e c to
a + be—Tp — p = 0. (4.5)
By Theorem 3.3 from [2] (see also Proposition 2.1 in [3]), the spectral bound s(A[A2]) of A[A2] is given by Re A1(A) + Re A2(A) (or —to if there are no 2 eigenvalues^).
Let us describe an approximation scheme for verification of the frequency inequality (4.1) in the case (4.4) and m = 2 (for simplicity). It can be easily adapted for the general case r2 = dim m and C: C([—t, 0]; r) — m, although the restrictions n = r1 = 1 are more significant and require many additional constructions. In what follows, we use the trigonometric basis in L2(—t, 0; c) given by 0k(0) = t—1/2eiT-l2nk0 for 0 e [—t, 0] and k e z.
(AS.1) Choose N > 0 and reals T > 0, Q > 0 and v0 > 0 such that —v0 > s(A[A2]) (see below (4.5));
(AS.2) For the linear equation x(t) = Axt compute the fundamental solution x^) on [—t,T] with xTO(0) = 0 for 0 e [—t,0) and xTO(0) = —^2B, and also for each k e {—N,..., N} compute the classical solutions xk(•) on [—t, T] with xk(0) = 0k(0) for 0 e [—t, 0];
(AS.3) For all p = —v0 + iw with w e [—Q, Q] compute the following:
(AS.3.1) For each k e {—N,..., N} and 0 e [—t, 0] compute
fT 1
Mk(0) := / e—pt^ xk(t — T0) • xTO(t + 0) — xk(t + 0) • xTO(t — T0)
dt; (4.6)
§In the case of (4.5), this means b = 0.
(AS.3.2) For all k, l e {—N,..., N} compute
r 0
ck :=V2J M^—(0)d0; (4.7)
(AS.3.3) Compute the largest singular value aT,N(p) of the matrix WT,N(p) with the entries c[ over k, l e {—N,..., N};
(AS.4) For all p = —v0 + iw with w e [—Q, Q] verify the inequality aTN(p) < A-i.
Let us discuss the choice of the parameters T, N and Q and their influence on dynamics of the approximation scheme. Firstly, for a given Q > 0, the approximations aT,N(—v0 + iw) convergence uniformly in w e [—Q, Q] as T, N ^ to to the norm of the transfer operator from (4.1) and they are globally Lipschitz in w e r with a uniform Lipschitz constant. Here, the parameter T can be taken relatively small, since the corresponding to it errors decay exponentially (in an appropriate norm) uniformly in w e r and N as T tends to infinity. For example, in the experiments conducted in [3], the graphs of aT,N(—v0 + iw) for T = 15 and T = 25 are indistinguishable. So, the convergence of aT,N(—v0 + iw) is mostly affected by the choice of N. In the experiments from [3], it is sufficient to take N =10 to obtain the convergence in w e [—Q, Q] for Q = 30 and the graph become indistinguishable with the corresponding graphs for N = 20 and N = 30.
Finally, the choice of Q is related to the conjecture from [3] concerned with that the norm of the operator from (4.1) behave as an aymptotically almost periodic (in the sense of Bohr) function as |w| ^ to. For practical purposes, this means that the inequality in (4.1) can be verified in a finite interval of w. This conjecture is justified by the conducted experiments.
We refer to Section 6 for applications of the Frequency Criterion to the classical Nicholson blowflies model and its comparison with other methods.
5 Nonautonomous Nicholson blowflies model: attractors and dimension estimates
Let (Q,$) be a semiflow on a compact metric space Q. We consider the nonautonomous Nicholson blowflies model which is described over q e Q as
x(t) = -J(^(q))x(t) + P(^(q))x(t - t)e-a(^i(q))x(t-T), (5.1)
where 5(•), P(•) and a(-) are continuous functions on Q satisfying
0 < 5— < 5(q), 0 < P(q) < P + and 0 < a— < a(q) for all q e Q. (5.2)
Let C+ be the cone of nonnegative functions in C([—t, 0]; r) and let BR denote the closed ball of radius R in the supremum norm.
Lemma 5.1. For any R > R0, where
P+
R0 = "T-, (5.3)
solutions x(t) = x(t; 00) to (5.1) with 00 e C+ nBR exist for all t > 0 and satisfy xt e C+ n BR for all t > 0.
Proof. Take x(t; 00) as in the statement for R > R0 and suppose that it exists on some interval [—t, t0) with t0 > 0. From (5.1) and (5.2), it is clearly seen that we must have xt e C+ for any t e [—t, t0).
o 0 _
Let Br be the interior of Br. Suppose 00 e Br. Then in the
case xt leaves
0
BR, we may suppose that t0 is such that x(t0) = R and 0 < x(t) < R for any t e [—t, t0). Note that the maximum of ye—a y over y > 0 is (ea—)—1. Then from (5.1) and (5.2) we have
P+
x(t0) < —5—R +-< 0 (5.4)
ea—
that leads to a contradiction. Since x(t; 00) depend continuously on 00, from this the conclusion of the lemma follows. The proof is finished. □
Similarly to Corollary 5.1 in [1], we obtain the following.
Corollary 1. There is a dissipative skew-product semiflow n in the space Q x C+ given by
nt(q, 00) := (^(q), xt) for t > 0, q e Q and 00 e C+, (5.5)
where x(t) = x(t; 00) is a solution of (5.1) with 00 e C+. Its global attractor A lies in Q x (C+ n BRo), where R0 is given by (5.3).
Now let y0(t) = x(t; 00) be a solution of (5.1) with 00 e C+. Then the linearized over the corresponding trajectory equation reads as
x(t) = —5(^(q))x(t) + P(q)) [1 — a(^t(q))y0(t — t)] e—a(^i(q))yo(t—T)x(t — t).
(5.6)
Thus, in terms of (3.6) for p = (q, 00) e Q x C+ we have
L>(p) = -J(q) and L-:(p) = P(q)[1 - a(q)^(-T)]e-a(q)^0(-T). (5.7)
Consequently, for a fixed k = k0 = kj, the eigenvalue Ai(p) of the 1 x 1-matrix from (3.6) is given by
Ai(p) = —2J(q) + eKT|P(q)|2(1 - a(q)^(-T^e-2^^) + 1. (5.8)
Note that (1 — ay)2e-2ay < 1 for any y > 0 and a > 0. So,
Ai(p) < —2J- + eKT|P+|2 + 1. (5.9)
Let S be the linear cocycle generated by (5.6) in h = r x L2(—t, 0; r) over the semiflow (Q x C+,n).
Theorem 2. Suppose^ that 1 - 2J- + |P + |2 > 0. Then
dimL S < |P +|2ep*+i • t + 1, (5.10)
where p* is the unique root p > —1 of
1 — 9A-
= . (5.11)
Proof. From (5.9) in terms of (3.7), we have
a+(m) < 1 [-2J- + eKT|P+|2 + 1 - k(m - 1)] =: a(m). (5.12)
2
Note that a(m) is a concave function of real m > 1. Then finding d* such that a(d*) = 0 would guarantee that dimL£ < d* (see Remark 3.2 in [1]). Clearly, such d* is given by
d* = ekT|P +|2 + 1 - - +1. (5.13)
*
k
Minimizing it over k > 0, we obtain
1 2J—
ekT(kt - = . (5.14)
where making the change k = t-i(p + 1) with p > —1 leads to (5.11). Note that p = —1 (or k = 0) corresponds to 1 — 2J- + |P+|2 = 0 and d* is well-defined in the limit k ^ 0+. The proof is finished. □
^Otherwise from (5.9) with k = 0 we have Ai(p) = 1 — 2^ + |P +|2 < 0 and, as a consequence, dimL S = 0.
Note that the upper estimate in (5.10) depends on the change*** x(t) ^ x(Kt) (with some k > 0) in (5.1), although the Lyapunov dimension does not. Let us derive from it the following estimate which is invariant w.r.t. such time scaling.
Corollary 2. Suppose that*** P + > 5—. Then
dimL S < +1 • P+t + 1, (5.15)
where p* is the unique root p > — 1 of
—
pep+1 = 1 — 2 _. (5.16)
Proof. Applying the change x(t) ^ x(1) in (5.10) gives the transformation in the coefficients as 5(-) ^ and P(•) ^ ^p^- and t ^ tP+. Then Lemma 2 applied after this change gives the desired. The proof is finished. □
However, (5.15) is not the maximum that can be achieved from Theorem 2 and in general from the change x(t) ^ x(Kt) we have the following.
Corollary 3. Suppose that P+ > 5—. Then
dimL S < inf (Kep*(K)+1 • |P+|2t + 1) , (5.17)
k>0 \ /
where p*(k) is the unique root p > —1 of
Pe"+1 = !—P+—• (5.18)
Let us consider the autonomous version of (5.1) given by
x(t) = —5x(t) + Px(t — t )e—ax(t—T), (5.19)
with some parameters 5 > 0, P > 0 and a > 0. This equation was suggested by W.S.C. Gurney, S.P. Blythe and R.M. Nisbet [19] as a model for laboratory insect populations studied by A.J. Nicholson and it showed a good quantitative agreement with his results.
Although for (5.19) one is usually interested in stable periodic oscillations, the model can demonstrate chaotic behavior. This possibility is not clearly indicated in the literature, although in the survey of Hastings et al. [21] it is noted that
***In terms of (2.10) considered as a transfer equation with the boundary condition (5.1), this is not only the change of time. More precisely, beside x(t) ^ x(«t) we also have ^(t, 0) ^ ^(«t, k0). So, the change is spatio-temporal.
tttFor P + < S- we have dimL S = 0 by the same reasons as above.
the original Nicholson investigations contain chaotic experimental data. On the other hand, in [19], nonperiodic oscillations in the region of linear instability were classified as "formally aperiodic" with a "single clearly marked dominant period".
For (5.19), chaotic behavior can be observed, for example, for a = t = 1, J = 9.5 and P = 250 (such parameters are included into the study from [19]). For such parameters, using the JiTCDDE package for Python to compute Lyapunov exponents (see G. Ansmann [12] and the repository for details), we observed the following values
Af = 0.2176 ± 0.0030, Af = 0.0577 ± 0.0029,
AL = -0.0006 ± 0.0037, Af = -0.1554 ± 0.0053 (5.20)
Af = -0.4462 ± 0.0059.
Thus, there are two positive Lyapunov exponents.
Moreover, for the above parameters, (5.15) gives the estimate dimL S < 883.821 and (5.17) (the minimum is achieved at k « 0.0015) gives dimL S < 671.004. From numerical experiments we have that the zero stationary state
has a 81-dimensional unstable manifold and the Lyapunov dimension at satisfies dimL 00 e (215,216). This contrasts to the Lyapunov dimension at the positive equilibrium 0+ which satisfies dimL 0+ e (16,17). Moreover, for P e [250,10000], the value dimL 00 shows an approximately linear growth as 0.865 • P and the dimension of the unstable manifold of 00 grows approximately as 0.318 • P. Thus (5.15) as well as the resulting estimates for the Hausdorff or fractal dimensions of A seem to be asymptotically sharp as P ^ to since ep*+i < e2. From this, we may conjecture the following.
Conjecture 1. For (5.19), the Eden conjecture at 00 is valid, i.e.
dimL S = dimL 00. (5.21)
6 Classical Nicholson blowflies model: comparison of results on the global stability
In this section, we will discuss results on the global stability of (5.19) in the positive cone C+ and compare them with applications of the Frequency Criterion discussed in Section 4.
It can be shown (see® [42]) that for J > P any solution starting in the cone
iiipor 6 > P, this can be also deduced from the fact that dimL S = 0 (see Corollary 2). Since dimH A <
C+ tends to the zero stationary state = 0. For 5 < P, loses its stability and the positive equilibrium 0+ = a—1 ln P appears in C+.
In [24], M.R.S. Kulenovic, G.F. Ladas and Y.G. Sficas proved that 0+ is globally attracting in the cone C+ provided that
(eST — 1)(P — 0 < 1. (6.1)
Later, their result was sharpened by J.W. So and J.S. Yu to
(eST — 1)ln^p^j < 1 (6.2) and L. Jingwen [23] complemented it with the two new regions given by
(eST — 1) ln fP>) < 1 + 4+, if «0+ > ^^, or \5 / a0+ 2
, st ip1 r VT+iVs — 1
(eST — p\ < i + —, ifa0+ >-2-
(6.3)
All the above methods are based on direct a priori estimates for solutions. In [18], I. Gyori and S. Trofimchuk involved the theory of one-dimensional maps with negative Schwarzian derivative into the problem and complemented the previous results by the condition
(1 — p) ln ^p^j < 1 — p + 2 (1 + \/1 + 4P(1 — P)) , where P = e—ST. (6.4)
It is however does not entirely cover (6.2) and (6.3).
Finally, E. Liz, V. Tkachenko and S. Trofimchuk further developed the method and obtained in [36] the region of global stability whose boundary is very close to the region of linear instability of 0+. Their result significantly improves the previous investigations and cannot be sharpened without taking into account that the delay in (5.19) is independent of time§§§.
It is however interesting to compare the Frequency Criterion (see below (4.3)) with the above mentioned results excluding the result of [36] which is unreachable for rough general methods.
Firstly, note that by the change x(t) ^ ax(Tt) we may assume in (5.19) that a = t = 1 (with the new 5 and P being 5t and Pt respectively in the old terms)
dimL S = 0, the global attractor cannot contain nonstationary points. There is only one stationary point in C+ and so we must have A = j^0}.
§§§It is known that their criterion is sharp in the class of time-dependent delays [36].
and the above inequalities are independent on such a change. Now let us consider the linearized equation along a solution y0 : [—1, ^ r of (5.19) as
±(t) = + p [1 - y0(t - 1)] e-yo(t-1)x(t - 1). (6.5)
Note that for any y G [0, 4] we have
,+ , ^ y , -e-2, if P > 2eJ, , N
1 =: d+ > (1 - y)e-y > d- := I P P> ' (6.6)
1(1 — p)e P, otherwise.
We put k := P(d- + d+)/2. Then (6.5) can be written as
±(t) = -Jx(t) + kx(t - 1) + F'(y0(t - 1))x(t - 1), (6.7)
where |F'(y)| < A := P — k for any y e [0, ef6].
Let A be the global attractor of n in C+ generated by (5.19) (see Corollary 1). Note that A is located in the ball of radius ^6.
Now let us apply the approximation scheme (AS.1)-(AS.4) to (6.7) considered in the context of (4.4) with P := A; A0 := —J0(0) + k0(—1) and C0 := 0(—1) for 0 e C([-1,0];r); B := 1 and the above given F' and A.
We use a numerical realization of the approximation scheme on Python (see the repository for details). Parameters of the scheme are taken as Q = 10, N = 10, T = 15 and v0 = 0.01. With such parameters, we verified the corresponding frequency conditions for (6.7) with J e (0.1, 5) and P e (1.5, 5). Results are collected in Fig. 1. We note that, although the Frequency Criterion seem to completely cover only the result of [24], it also covers a region which is not covered by the other criteria considered^". Moreover, the criterion fails for large P/J since the bound for the global attractor and, consequently, the Lipschitz constant A become rough.
Data availability
The data that support the findings of this study can be generated using the scripts in the repository:
https://gitlab.com/romanov.andrey/nicholson-blowflies-experiments
^^Note that the result of [36] (excluded in Fig. 1) covers all the considered space of parameters.
Figure 1: Regions in the (P, 5)-space with 5 E (0.1, 5) and P E (1.5,5) covered by the Frequency Criterion (gray) verified via the approximation scheme (AS.1)-(AS.4); by (6.1) (blue); by (6.2) (orange); by (6.3) (green); by (6.4) (red). For each criterion, except the Frequency Criterion, we do not color the region covered by the preceding ones.
Funding
The reported study was funded by the Russian Science Foundation (Project 2211-00172).
References
[1] Anikushin M.M. Variational description of uniform Lyapunov exponents via adapted metrics on exterior products. arXiv preprint, arXiv:2304.05713 (2024)
[2] Anikushin M.M. Spectral comparison of compound cocycles generated by delay equations in Hilbert spaces. arXiv preprint, arXiv:2302.02537 (2024)
[3] Anikushin M.M., Romanov A.O. Frequency conditions for the global stability of nonlinear delay equations with several equilibria. arXiv preprint, arXiv:2306.04716 (2024)
[4] Anikushin M.M. Frequency theorem and inertial manifolds for neutral delay equations. J. Evol. Equ, 23, 66 (2023)
[5] Anikushin M.M., Romanov A.O. Hidden and unstable periodic orbits as a result of homoclinic bifurcations in the Suarez-Schopf delayed oscillator and the irregularity of ENSO. Phys. D: Nonlinear Phenom, 445, 133653 (2023)
[6] Anikushin M.M. Nonlinear semigroups for delay equations in Hilbert spaces, inertial manifolds and dimension estimates. Differ. Uravn. Protsessy Upravl., 4, (2022)
[7] Anikushin M.M. Inertial manifolds and foliations for asymptotically compact cocycles in Banach spaces. arXiv preprint, arXiv:2012.03821 (2022)
[8] Anikushin M.M. Frequency theorem for parabolic equations and its relation to inertial manifolds theory. J. Math. Anal. Appl., 505(1), 125454 (2021)
[9] Anikushin M.M. Almost automorphic dynamics in almost periodic cocycles with one-dimensional inertial manifold. Differ. Uravn. Protsessy Upravl., 2, (2021), in Russian
[10] Anikushin M.M. On the Liouville phenomenon in estimates of fractal dimensions of forced quasi-periodic oscillations. Vestn. St. Petersbg. Univ., Math., 52(3), 234-243 (2019)
[11] Anikushin M.M., Reitmann V., Romanov A.O. Analytical and numerical estimates of the fractal dimension of forced quasiperiodic oscillations in control systems. Differ. Uravn. Protsessy Upravl., 2 (2019), in Russian
[12] Ansmann G. Efficiently and easily integrating differential equations with JiT-CODE, JiTCDDE, and JiTCSDE. Chaos, 28(4), 043116 (2018)
[13] Carvalho A.N., Langa J.A., Robinson J.C. Attractors for Infinite-Dimensional Non-Autonomous Dynamical Systems. Springer Science & Business Media (2012)
[14] Chepyzhov V.V., Ilyin A.A. On the fractal dimension of invariant sets; applications to Navier-Stokes equations. Discrete Contin. Dyn. Syst., 10(1&2) 117-136 (2004)
[15] Chicone C. Inertial flows, slow flows, and combinatorial identities for delay equations. J. Dyn. Diff. Equat, 16(3), 805-831 (2004)
[16] Chicone C. Inertial and slow manifolds for delay equations with small delays. J. Differ. Equations, 190(2), 364-406 (2003)
[17] Chueshov I.D. Dynamics of Quasi-stable Dissipative Systems. Berlin: Springer (2015)
[18] Gyori I., Trofimchuk S. Global attractivity in dx/dt = —Jx + pf(x(t — t)). Dynam. Syst. Appl., 8, 197-210 (1999)
[19] Gurney W.S.C., Blythe S.P., Nisbet R.M. Nicholson's blowflies revisited. Nature, 287(5777), 17-21 (1980)
[20] Hale J.K. Theory of Functional Differential Equations. Springer-Verlag, New York (1977)
[21] Hastings A., Hom C.L., Ellner S., Turchin P., Godfray H.C.J. Chaos in ecology: is mother nature a strange attractor? Annu. Rev. Ecol. Evol. Syst., 24(1), 1-33 (1993)
[22] Hu W., Caraballo T. Hausdorff and fractal dimensions of attractors for functional differential equations in Banach spaces. J. Differ. Equ., 385, 395-423 (2024)
[23] Jingwen L. Global attractivity in Nicholson's blowflies. Appl. Math., 11(4), 425-434 (1996)
[24] Kulenovic M.R.S., Ladas G.F., Sficas Y.G. Global attractivity in Nicholson's blowflies. Appl. Anal, 43(1-2), 109-124 (1992)
[25] Kuznetsov N.V., Reitmann V. Attractor Dimension Estimates for Dynamical Systems: Theory and Computation. Switzerland: Springer International Publishing AG (2020)
[26] Kuznetsov N.V., Mokaev T.N., Kuznetsova O.A., Kudryashova E.V. The Lorenz system: hidden boundary of practical stability and the Lyapunov dimension, Nonlinear Dyn., 102, 713-732 (2020)
[27] Kuznetsov N.V. The Lyapunov dimension and its estimation via the Leonov method. Phys. Lett. A, 380(25-26) 2142-2149 (2016)
[28] Leonov G.A., Kuznetsov N.V., Korzhemanova N.A., Kusakin D.V. Lyapunov dimension formula for the global attractor of the Lorenz system. Commun. Nonlinear Sci. Numer. Simul., 41, 84-103 (2016)
[29] Leonov G.A. Lyapunov dimension formulas for Lorenz-like systems. Doklady Mathematics, 93(3) (2016)
[30] Leonov G., Mokaev T. Lyapunov dimension formula for the attractor of the Glukhovsky-Dolzhansky system. Doklady Mathematics, 93(1) (2016)
[31] Leonov G.A., Alexeeva T.A., Kuznetsov N.V. Analytic exact upper bound for the Lyapunov dimension of the Shimizu-Morioka system. Entropy, 17(7), 5101 (2015)
[32] Leonov G.A. Lyapunov dimensions formulas for Henon and Lorenz attrac-tors. Alg. Anal., 13, 155-170 (2001) (Russian); English transl. St. Petersburg Math. J., 13(3), 453-464 (2002)
[33] Leonov G.A., Boichenko V.A. Lyapunov's direct method in the estimation of the Hausdorff dimension of attractors. Acta Appl. Math., 26(1), 1-60 (1992)
[34] Li M.Y., Muldowney J.S. A geometric approach to global-stability problems. SIAM J. Math. Anal., 27(4), 1070-1083 (1996)
[35] Li M.Y., Muldowney J.S. Lower bounds for the Hausdorff dimension of attractors. J. Dyn. Differ. Equ., 7(3), 457-469 (1995)
[36] Liz E., Tkachenko V., Trofimchuk S. A global stability criterion for scalar functional differential equations. SIAM J. Math. Anal., 35(3), 596-622 (2003)
[37] Louzeiro M., Kawan C., Hafstein S., Giesl P., Yuan J. A Projected Subgradient Method for the Computation of Adapted Metrics for Dynamical Systems. SIAM J. Appl. Dyn. Syst, 21(4), 2610-2641 (2022)
[38] Mallet-Paret J., Nussbaum R.D. Tensor products, positive linear operators, and delay-differential equations. J. Dyn. Diff. Equat., 25, 843-905 (2013)
[39] Mallet-Paret J., Sell G.R. The Poincare-Bendixson theorem for monotone cyclic feedback systems with delay. J. Differ. Equations, 125, 441-489 (1996).
[40] Mallet-Paret J. Negatively invariant sets of compact maps and an extension of a theorem of Cartwright. J. Differ. Equations, 22(2), 331-348 (1976)
[41] Smith R.A. Some applications of Hausdorff dimension inequalities for ordinary differential equations. P. Roy. Soc. Edinb. A, 104(3-4), 235-259 (1986)
[42] So J.W., Yu J.S. Global attractivity and uniform persistence in Nicholson's blowflies. Differ. Equ. Dyn. Syst, 2(1), 11-18 (1994)
[43] So J. W-H., Wu J. Topological dimensions of global attractors for semilinear PDE's with delays. Bull. Aust. Math. Soc., 43(3), 407-422 (1991)
[44] Temam R. Infinite-Dimensional Dynamical Systems in Mechanics and Physics. Springer (1997)
[45] Zelik S. Attractors. Then and now. Uspekhi Mat. Nauk, 78(4), 53-198 (2023)
