Calculating the Lyapunov exponents of a piecewise-smooth soft impacting system with a time-delayed feedback controller
aa r X i v : . [ n li n . C D ] J un Calculating the Lyapunov exponents of a piecewise-smooth softimpacting system with a time-delayed feedback controller
Zhi Zhang, Yang Liu ∗ , Jan Sieber a College of Engineering Mathematics and Physical Sciences, University of Exeter, Harrison Building, North Park Road,Exeter EX4 4QF, UK
Abstract
Lyapunov exponents are a widely used tool for studying dynamical systems. When calculating Lyapunovexponents for piecewise-smooth systems with time-delayed arguments one faces a lack of continuity in thevariational problem. This paper studies how to build a variational equation for the efficient constructionof Jacobians along trajectories of the delayed nonsmooth system. Trajectories of the piecewise-smoothsystem may encounter a so-called grazing event where the trajectory approaches a discontinuity surfacein the state space in a non-transversal manner. For this event we develop a grazing point estimationalgorithm to ensure the accuracy of trajectories for the nonlinear and the variational equations. We showthat the eigenvalues of the Jacobian matrix computed by the algorithm converge with an order consistentwith the order of the numerical integration method, therefore guaranteeing the reliability of our proposednumerical method. Finally, the method is demonstrated on a periodically forced impacting oscillatorunder the time-delayed feedback control.
Keywords:
Lyapunov exponents; Piecewise-smooth dynamical system; Delay differential equation;Grazing; Impact oscillator.
1. Introduction
Analysing grazing events for nonsmooth systems is a challenging task [1]. In general, vibro-impactsystems, such as ship mooring interactions [2], bearing looseness [3] and the multi-degree-of-freedomimpact oscillators [4], may have abundant coexisting attractors when grazing occurs. Tiny differences inmodelling will lead to different motion of the system [5–7]. For example, the motion of an impact oscillatorexperiences significant change due to a slight variation on its parameter when a grazing bifurcation isencountered [8]. In [9], Nordmark studied the characteristic scaling behaviour near grazing bifurcations,and used the self-similarity under scaling to derive a renormalised mapping. Nordmark [10] presentedthe grazing bifurcation of a hard impact oscillator, where the Poincar´e map has a singular Jacobian,by using a first-order Taylor expansion. It has shown that the stability of the oscillator can be studiedmore precisely if its grazing events are computed more accurately. This paper will study a new methodto improve the accuracy of calculating the grazing events by estimating the impacting moment of thesystem. Based on this accurate grazing trajectory, stability analysis of the system can be carried out.In many applications [11–17] arise differential equations in which the derivative of the unknownfunctions at a certain time depends on the value of the function at previous time. These are so-called
Email addresses: [email protected] (Zhi Zhang), [email protected] (Yang Liu ∗ ), [email protected] (JanSieber) Preprint submitted for review June 26, 2020 elay differential equations (DDEs). For example, Zhang et al . [13] studied a delayed pest control modelwhich was a high-dimensional differential equation with impulsive effects at different fixed impulsivemoments. In [14], Carvalho and Pinto used a mathematical model with delay to describe the dynamics ofAIDS- related cancers with the treatment of HIV and chemotherapy. In [15], Yan et al . used the basin ofa time-delayed system modelling cutting process to determine the unsafe cutting zone. The above studiesare concerned with smooth DDEs. The analysis of nonsmooth DDEs is more challenging due to the lackof an accurate algorithm for computing the grazing events. Until now, there are very few systematicstudies regarding to nonsmooth DDEs, which is the focus of this paper. The present work will study anew algorithm to determine the occurrence of grazing for improving computational accuracy and a newmethod for calculating Lyapunov exponents (LEs) along the trajectories of a nonsmooth DDE.The LE of a trajectory is a quantity that characterises the rate of separation of infinitesimally near-by trajectories [18]. It determines a notion of sensitivity of this trajectory to perturbations in initialconditions. If the largest LE, which is referred to the maximal LE, is greater than zero, any smallperturbation of the initial condition will result in an exponential divergence of the resulting perturbedtrajectory until the distance between the perturbed and unperturbed trajectories is no longer small.This sensitivity with respect to initial condition is one of the defining features of chaos. If the LEare identical for typical trajectories of an attractor in a dynamical system, one speaks of the LE forthis attractor (or this dynamical system). The LE indicate predictability (or lack of it) for dynamicalsystems, such that they are considered as an important tool for studying the stability of dynamicalsystems. Therefore, the development of an efficient method for calculating the LEs of dynamical systemis an active area of research, see e.g. [19–26]. For finite-dimensional dynamical systems Benettin et al .[20] introduced a systematic method for estimating the LEs of smooth dynamical systems. Wolf et al .[21] developed a method for extracting the largest LE from an experimental time series. For nonsmoothsystems, M¨uller [24] developed a model-based algorithm to calculate the LEs of nonlinear dynamicalsystems with discontinuities. They found that the required linearised equations must be supplementedby certain transition conditions when crossing the discontinuities. In [25], Dellago et al . generalisedBenettin’s classical algorithm and applied it to the case of dynamical systems where smooth streamingwas interrupted by a differentiable map at discrete times. Lamba and Budd [27] have shown that thelargest LE has a discontinuous jump at grazing bifurcations in Filippov systems and scales like 1 / | ln ǫ | ,where ǫ is the bifurcation parameter. In contrast to ordinary differential equatrions (ODEs), DDEsare infinite dimensional systems. such that the computation of LEs for nonsmooth DDEs combinesdifficulties from discontinuities and high dimensionality. In principle, a DDE could be approximated bya high-dimensional ODE, which can be linearised along trajectories obtained by numerical integration[28, 29], such that the LEs can be constructed for the Poincar´e map. Studies by Repin [30] and Gy¨ori andTuri [31] have shown that DDEs can be analysed using approximating high-dimensional ODEs. However,if the delay time is large, calculating the LEs of nonsmooth DDEs needs to store excessive history datapoints during delay period compared to smooth DDEs [32–34], e.g. the data at past encounters of thediscontinuity. In this case the global convergence of the system cannot be guaranteed. Therefore, it maycause inaccuracy in calculating the eigenvalues of Jacobian matrix which is used for estimating the LEsof nonsmooth DDEs.The contribution of the present work is the development of a novel method for precisely calculatingthe LEs of piecewise-smooth differential equations with a delayed argument, which can provide improvedaccuracy for stability analysis of periodic orbits. In detail, if an algorithm cannot estimate the point ofdiscontinuity along trajectory with an accuracy of the same order as its integration method, especially2n the grazing event, the expected discontinuous coefficients of the variational problem will have unex-pectedly low accuracy leading to an accumulation of errors. Similar work was reported by M¨uller [24]who studied a method for constructing the map of systems with discontinuity, and combined it withthe map obtained along the differentiable parts of the trajectories to generate a composition of Jacobianmatrices for calculating LEs. However, M¨uller’s approach is difficult to implement for piecewise-smoothDDE due to its high dimension and complex dynamics, which could cause a high computational costand an accumulation of computational errors at discontinuous moments. We address this issue in thepresent work, demonstrating our approach for the delayed piecewise-smooth oscillator. We construct aPoincar´e map that consists of many local maps for each small time step, which are linearised for the LEcomputation. As the linearised Poincar´e map requires accurate information about the time of crossingor grazing of a discontinuity (when impact occurs), we will introduce a grazing estimation algorithm toobtain an accurate Jacobian matrix for the oscillator. The novelty of our proposed method is that it canestimate the point of discontinuity locally along trajectories of piecewise-smooth DDEs, improving theaccuracy of computations of the system trajectory and of the LEs. The proposed method can also beextended to other nonsmooth dynamical systems, such as the hard impact oscillator with a time-delayedcontroller or stick-slip vibrations with a delay term. To demonstrate the reliability of the method, wewill carry out an error analysis for the nonzero eigenvalues of the Jacobian by adopting the spectralapproximation methods introduced by Chatelin [35] and Breda et al. [32, 34]. Our study indicates thatthe proposed method can reduce the error for the nonzero eigenvalues of the Jacobian by increasing thedimensions of the system of ODEs approximating the DDE slightly, which is generated by linearising theDDEs along trajectories obtained by numerical integration.The rest of this paper is organised as follows. Section 2 introduces the mathematical model of aperiodically forced mechanical oscillator subjected to a one-sided soft impact. This is followed by somebasic relevant definitions and preparations. Section 3 presents the method for constructing the Jacobianof Poincar´e map of piecewise-smooth DDEs. However, such a construction is inaccurate due to thenonsmoothness of the considered system. Thus, Section 4 studies an estimation method for determiningthe points of discontinuity accurately. Here, two cases of grazing events are considered based on thegeometry of the trajectory. Section 5 uses linear operator theory to carry out an error analysis for theeigenvalues of the Jacobian, which can validate the reliability of our proposed method. In Section 6,the steps for computing LEs are detailed. Examples and several control scenarios of the oscillator arepresented in Section 7 to demonstrate the accuracy of the method. Finally, some concluding remarks aredrawn in Section 8.
2. Mathematical model and relevant preparations
The impact oscillator shown in Fig. 2.1 represents a mechanical system encountering intermittent so-called soft impacts, which will be studied in the present work. Soft impacts occur in mechanical systemswhen an object hits an obstacle of negligible mass but non-negligible stiffness. In Fig. 2.1 the object ismodelled by the block of mass m and the obstacle is modelled by the spring with stiffness k (a backlash spring). The collision occurs when the distance g between block and spring reaches 0. Since at impactthe spring is relaxed, the forces in the system depend continuously on g (and, hence, on the position y ofthe block), but the spring constants exerted by the backlash spring are discontinuous: 0 for g > k for g = 0. Systems with soft impacts are common to a broad range of engineering applications, e.g. [36–39],where the repeated collision of mechanical parts is unavoidable [40]. The vibro-impact capsule system[41, 42] is a typical two-degrees-of-freedom dynamical system experiencing soft impacts and nonlinear3riction. Any small variations in friction or system parameters (e.g. the stiffness of the backlash spring)may lead to a qualitative change of the dynamics of the system [43, 44]. Thus, accurate prediction ofits collision is crucial to fully understand the dynamics of the system, in particular, in the presence oftime-delay effects [45].The nondimensional equations of motion of the impact oscillator can be written in a compact formas below [8], x ′ ( τ ) = v ( τ ) ,v ′ ( τ ) = aω sin( ωτ ) − ζv ( τ ) − x ( τ ) − β ( x ( τ ) − e ) H ( x ( τ ) − e ) , (2.1)where H ( · ) stands for the Heaviside step function and x ′ , v ′ denote differentiation with respect to thenondimensional time τ . The discontinuity boundary is fixed at x = e , with e > ω n = r k m , τ = ω n t, ω = Ω ω n , ζ = c mω n ,x = yy , e = gy , a = Ay , β = k k , where y > ω n is the natural angular frequency of the mass-springsystem ( m , k in Fig. 2.1), ω is the ratio between forcing forcing frequency and natural frequency, β isthe stiffness ratio, ζ is the damping ratio, and a is the nondimensionalised forcing amplitude. k c m Asin( t) W g k y Figure 2.1: Physical model of the soft impact oscillator [40].
In the present work, we will consider a control signal u ( τ ) , τ ≥
0, which will be superimposed on thesystem’s external excitation as follows x ′ ( τ ) = v ( τ ) ,v ′ ( τ ) = (cid:0) aω sin( ωτ ) + u ( τ ) (cid:1) − ζv ( τ ) − x ( τ ) − β ( x ( τ ) − e ) H ( x ( τ ) − e ) , (2.2)where u ( τ ) = k (cid:0) v ( τ − τ d ) − v ( τ ) (cid:1) , τ ≥ , (2.3)defines the proportional feedback controller that feedbacks the difference between the current measure-ment of v and a measurement of v from some time τ d ago [46]. In the expression above, k ≥ τ d > τ d = 2 π/ω and if (2.3) successfully stabilises a period-1 motion. This isthe case even if we do not know the precise time profile of this period-1 motion, which is in contrast tostandard linear feedback control u ( τ ) = k ( v ref ( τ ) − v ( τ ). The asymptotically vanishing control signal isattractive in applications where energy consumption is a critical issue, e.g. [47].Eq. (2.2) can be rewritten in the form of a general piecewise continuous DDE with a periodic externalexcitation as ˙ y ( t ) = f ( y ( t ) , y ( t − τ d )) + p ( t ) , for H ( y ( t ) , e ) > , ˙ y ( t ) = f ( y ( t ) , y ( t − τ d )) + p ( t ) , for H ( y ( t ) , e ) < ,y ( t + ) = y ( t − ) , for H ( y ( t ) , e ) = 0 , (2.4)where f , : R d × R d → R d , H : R d → R are sufficiently smooth functions and p : R + → R d is smoothand periodic with the period T >
0. The delay τ d is assumed to be positive but may be different fromthe period in general. In the present work, we only consider one single delay in the system for simplicity,and assume that for any y, ¯ y, y d , ¯ y d ∈ R d , f , f and H satisfy the Lipschitz condition | f , ( y, y d ) − f , (¯ y, ¯ y d ) | ≤ l | y − ¯ y | + l | y d − ¯ y d | , | H ( y, e ) − H (¯ y, e ) | ≤ l | y − ¯ y | ,where l , l , l ≥ | · | is a norm on R d . We assume that the initial condition is a suitable initialfunction on [ t , t − τ d ]. The general form (2.4) belongs to the class of hybrid dynamical systems [1],which consists of a flow (in our case only forward in time), combined with discrete events.Take N ∈ Z + sufficiently large, and define the discretisation grid points τ id := i τ d N , i = 0 , . . . , N ,and u i ( t ) := y ( t − τ id ) for all t ≥ i = 0 , . . . , N . Eq. (2.4) can be approximated by a d ( N + 1)dimensional piecewise-smooth discretised problem studied in [30], which will be presented in Section 3.This approximation method has also been studied by Krasovskii [48], finding that the solution of theapproximating system uniformly converges to the solution of the original DDEs when N → ∞ . Byusing the same approach, Gy¨ori and Turi [31] and Banks [49] carried out convergence analyses for twoDDEs. Breda et al. [33] studied the characteristic roots of linear DDEs, and used a Runge-Kutta methodto construct a high-dimensional approximating system. The nonzero eigenvalues of evolution operatorswere computed through a pseudospectral collection, which was used to analyse the asymptotic stability ofDDEs. Since Eq. (2.4) is a piecewise-smooth DDE whose trajectories can encounter discontinuities, themethods used for smooth DDEs are not suitable, or, at least, converge with lower-than-expected order.Therefore, motivated by the periodic forcing of Eq. (2.4), our plan here is to derive a Poincar´e map (alsocalled stroboscopic map) for discretising the system and study linear stability of its orbits by consideringthe Jacobian matrix of the map in these orbits. After such a reduction to the Poincar´e map, we will beable to define LEs for this time-discrete map.For the piecewise DDE (2.4), we consider a constant phase surface as the Poincar´e section defined by P Ts := { ( y, t ) ∈ C ([ − τ d , , R d ) × R + | t = t + kT, k ∈ Z + } . For the corresponding Poincar´e map P : P Ts → P Ts (2.5)the LEs can be defined as follows. Definition 2.1. [19] For any initial condition x ∈ P Ts , let { x m } ∞ m =0 be the corresponding orbit of5he map P , and let λ m , · · · , λ mn be the n largest in modulus eigenvalues of DP m ( x ), sorted such that | λ m | ≥ . . . ≥ | λ mn | . The Lyapunov exponents of x are ϑ i := lim m →∞ ln | λ mi | m , i = 1 , . . . , n (2.6)whenever the limit exists for x and for all i ≤ n .The above definition is applicable to our map P acting on the infinite dimensional space P Ts , since P is differentiable and its linearisation is bounded and has a spectrum consisting only of a sequence (finiteor infinite) of eigenvalues of finite multiplicity converging to 0 and zero. The expression in the limit (2.6)is not a practical recipe for computation since λ mi may be very large or very small.
3. Constructing the Jacobian matrix of the Poincar´e map
For the nonsmooth system with a delay τ d smaller than its forcing period T , i.e. 0 < τ d < T , theperiod T can be written as T = nτ d + ∆ t , for some n ∈ Z + and ∆ t ∈ [0 , τ d ). For any time interval[ t m , t m + τ d ], where t m = t + ( m − T , t = t and m ∈ Z + , the solution of system (2.4) can beapproximated by N steps of size h = τ d N by using numerical integration. The expression derived in thissection initially ignore grazing of the discontinuity surface { H = 0 } . Section 4 will explain how theexpressions will be modified at the respective events. The modified Euler integration formula [50] givesfor a single step of size h = τ d /Nu ( t m + h ) = u ( t m ) + h (cid:2) f j ( u ( t m ) , u ( t m − hN ))+ f j ( u ( t m + h ) , u ( t m − h ( N − (cid:3) + h (cid:2) p ( t m ) + p ( t m + h ) (cid:3) , (3.1)(here written only for the first step at t m ) where j = 1 , if H ( u ( t m ) , e ) > ,j = 2 , if H ( u ( t m ) , e )) < ,u ( t + m ) = u ( t − m ) , if H ( u ( t m ) , e ) = 0 . Iterating this map N + 1 times gives a discretised map for the delay-time interval [ t m , t m + τ d ], whichwe call P d : R d ( N +1) → R d ( N +1) . It satisfies U m, = P d ( U m, ) , (3.2)where U m, := ( u TN ( t m ) , · · · , u T ( t m ) , u T ( t m )) T ∈ R d ( N +1) and U m, := ( u TN ( t m + τ d ) , · · · , u T ( t m + τ d ) , u T ( t m + τ d )) T ∈ R d ( N +1) , and we use the general convention that u i ( t ) = u ( t − ( i/N ) τ d ) forarbitrary i ∈ { , . . . , N } and t . Iterating the map P d n times, we can obtain a map P nd from U at time t m to U at time t m + nτ d , U m,n = P d ◦ · · · ◦ P d ( U m, ) = P nd ( U m, ) , (3.3)where U m,i := ( u TN ( t m + ihN ) , · · · , u T ( t m + ihN )) T ∈ R d ( N +1) . Finally the discretised map for the time∆ t is defined as P ∆ t : R d ( N +1) → R d ( N +1) , which can be represented as U m,n +∆ N = P ∆ t ( U m,n ) , (3.4)6here U m,n +∆ N := ( u TN ( t m + h ( nN + ∆ N )) T , · · · , u T ( t m + h ( nN + ∆ N )) T ∈ R d ( N +1) and ∆ N := ∆ th .Thus combining Eqs. (3.3) and (3.4) we can construct map P disc as the discretised Poincar´e map P advancing by time T U m,n +∆ N = P disc ( U m +1 , ) = P ∆ t ◦ P nd ( U m, ) , (3.5)which can then iterate further by setting U m +1 , = U m,n +∆ N . For an arbitrary perturbation δU isapplied, the variational equation for P disc can be written as δU m +1 , = N +1 X i =1 ∂P disc ( U m, ) ∂u i − ( t m ) δu i − ( t m ) , (3.6)where δU m, := ( δu TN ( t m ) , · · · , δu T ( t m ) , δu T ( t m )) T ∈ R d ( N +1) , and we use again the convention that δu i ( t ) := δu ( t − τ id ) , i = 0 , · · · , N . In fact, Eq. (3.6) can be obtained from discretising the continuousvariational equation of system (2.4), and its form can be obtained asdd t δu ( t ) = ∂f j ( t, u ( t ) , u N ( t )) ∂u δu ( t ) + ∂f j ( t, u ( t ) , u N ( t )) ∂u N δu N ( t ) , (3.7)where j = 1 , if H ( u ( t ) , e ) > ,j = 2 , if H ( u ( t ) , e ) < ,u ( t + ) = u ( t − ) , if H ( u ( t ) , e ) = 0 . An example initial function φ δ for (3.7) is of the form φ δ ( t ) = ( ǫ, , · · · , T ∈ R d and φ δ ( t ) =(0 , · · · , T ∈ R d for t ∈ [ t − τ d , t ), and sufficiently small ǫ . Discretising Eq. (3.7) in the interval[ t m , t m + nτ d ] by using the modified Euler integration gives δu ( t m + lh ) = δu ( t m + ( l − h )+ h (cid:2) A m,l δu ( t m + ( l − h + B m,l δu ( t m − ( N − l + 1) h )] (3.8)+ h [ A m,l +1 δu ( t m + lh ) + B m,l +1 δu ( t m − ( N − l ) h ) (cid:3) , where l = 1 , · · · , N, · · · , nN + ∆ N , A m,l = ∂f j ( u ( t ) ,u N ( t )) ∂u | t = t m + h ( l − , B m,l = ∂f j ( u ( t ) ,u N ( t )) ∂u N | t = t m + h ( l − and m ∈ Z + . Rewriting Eq. (3.8) in a matrix form gives δu N ( t m + lh )... δu ( t m + lh ) δu ( t m + lh ) = M m,l δu N ( t m + ( l − h )... δu ( t m + ( l − h ) δu ( t m + ( l − h ) , (3.9)where M m,l = ˆ M m,l ˜ M m,l , M m,l = I · · · · · · I − h B m,l +1 · · · I − h A m,l +1 − , and ˜ M m,l = I · · · · · · I h B m,l · · · I + h A m,l . By using the map (3.2), the matrix form of the variational equation (3.9) can be rewritten as δU m,N = M m,N ◦ · · · ◦ M m, δU m, . Since we have n maps, combining all the maps for the interval [ t m , t m + T ] gives δU m,n = M m,nN ◦ · · · ◦ M m, ◦ M m, δU m, . In addition, the map P △ t for the interval [ t m + nτ d , t m + T ] can be written as δU m +1 , = M m,nN +∆ N ◦ · · · ◦ M m,nN δU m,n . (3.10)Finally, the overall variational equation can be obtained as δU m +1 , = M m δU m, , (3.11)where M m = M m,nN +∆ N ◦ · · · ◦ M m,nN ◦ · · · ◦ M m, is the approximation of Jacobian matrix of thePoincar´e map P .Similarly, for the system with a large delay time, e.g. τ d ≥ T , the solution of system (2.4) can beapproximated by N steps f size h = τ d N by using numerical integration, which can be considered as aspecial case of the nonsmooth system with a small delay time (0 < τ d < T ) when n = 0. Let N T = Th bethe sample number for one period T , construct the map P d , and combine all the linearised maps at theinterval [ t m , t m + T ]. Finally, we can obtain the same variational equation as Eq. (3.11) and the Jaocbianmatrix of the Poincar´e map P .
4. Modifying the algorithm at the discontinuity
In this section, we will discuss a special phenomenon of the impact oscillator, the so-called crossingand grazing events. Since the system has rich complex dynamics when it experiences grazing [5, 51], acareful consideration in calculating this discontinuous moment is required. In addition, the global errorof our proposed algorithm will depend on how accurately we capture the effect of switching, as the errormade at the switching boundary could accumulate, leading to unexpected large global error. Therefore,during the grazing event, we need to modify our proposed algorithm from Section 3 by considering thetwo grazing cases illustrated in Fig. 4.1. 8 v x=e (x(t ),v(t )) * * (x(t +h),v(t +h)) * * (a) (b) (x(t +h),v(t +h)) * * (x(t ),v(t )) * * x=e x v Figure 4.1: (a) Case 1: for t = t ∗ >
0, such that H := H ( u ( t ∗ ) , e ) < H := H ( u ( t ∗ + h ) , e ) > H > H < t = t ∗ >
0, and there exists δt ∈ (0 , h ), such that H := H ( u ( t ∗ ) , e ) < H := H ( u ( t ∗ + h )) < H cr, := H ( u ( t ∗ + δt ) , e ) = 0 (or H > H > H cr, = 0). For Case 1, we assume that for time step l ∗ ∈ Z + at time t ∗ := t m + ( l ∗ − h the switching function H changes sign: H := H ( u ( t ∗ ) , e ) < H := H ( u ( t ∗ + h ) , e ) >
0, or H > H <
0. Thus, weexpect that for some time δt ∈ (0 , h ), the switching fucntion is zero: H cr, := H ( u ( t ∗ + δt ) , e ) = 0. Inorder to guarantee the order of convergence of our proposed algorithm to O ( h ), the crossing time δt needs to be estimated first. Since δt < h , the condition H cr, = 0 can be linearised as H cr, ≈ H ( u ( t ∗ ) + ˙ u ( t ∗ ) δt, e ) ≈ H + ddu H [ ˙ u ( t ∗ ) δt ] = 0 , such that δt = − H ddu H [ ˙ u ( t ∗ )] . (4.1)Once δt is calculated, the switching time t ∗ + δt can be obtained, and the variational equation at thestep crossing the switching can be written as δu ( t ∗ + δt ) = δu ( t ∗ ) + δt [ A m,l ∗ δu ( t ∗ ) + B m,l ∗ δu N ( t ∗ )]+ δt [ A δtm,l ∗ δu ( t ∗ + δt ) + B δtm,l ∗ δu N ( t ∗ + δt )] , (4.2)where A δtm,l ∗ = ∂f j ( u ( t ) ,u N ( t )) ∂u | t = t − m + h ( l ∗ − δt , B δtm,l ∗ = ∂f j ( u ( t ) ,u N ( t )) ∂u N | t = t − m + h ( l ∗ − δt and l ∗ = 1 , · · · , N, · · · , nN + ∆ N . Thus the discretied map from t ∗ to t ∗ + δt can be written as δu N ( t ∗ + δt )... δu ( t ∗ + δt ) δu ( t ∗ + δt ) = M δtm,l ∗ δu N ( t ∗ )... δu ( t ∗ ) δu ( t ∗ ) , (4.3)where M δtm,l ∗ = ˆ M δtm,l ∗ ˜ M δtm,l ∗ , M δtm,l ∗ := I · · · · · · I − δt B δtm,l ∗ · · · I − δt A δtm,l ∗ − , ˜ M δtm,l ∗ := I · · · · · · I δt B m,l ∗ · · · I + δt A m,l ∗ ,A m,l ∗ = ∂f j ( u ( t ) ,u N ( t )) ∂u | t = t m + h ( l ∗ − and B m,l ∗ = ∂f j ( u ( t ) ,u N ( t )) ∂u N | t = t m + h ( l ∗ − . It is worth noting that δu i ( t ∗ + δt ) can be approximated through linear interpolation based on the historical data obtained fromthe delayed time interval which also includes the grazing data.Similarly, for the time interval [ t ∗ + δt, t ∗ + h ], we can obtain δu N ( t ∗ + h )... δu ( t ∗ + h ) δu ( t ∗ + h ) = ¯ M hm,l ∗ δu N ( t ∗ + δt )... δu ( t ∗ + δt ) δu ( t ∗ + δt ) , (4.4)where ¯ M hm,l ∗ := ˆ M hm,l ∗ ˜ M hm,l ∗ ,˜ M hm,l ∗ := I · · · · · · I h − δt B δtm,l ∗ · · · I + h − δt A δtm,l ∗ and ˆ M hm,l ∗ := I · · · · · · I − h − δt B m,l ∗ +1 · · · I − h − δt A m,l ∗ +1 − . Finally, we have δu N ( t ∗ + h )... δu ( t ∗ + h ) δu ( t ∗ + h ) = ¯ M hm,l ∗ M δtm,l ∗ δu N ( t ∗ )... δu ( t ∗ ) δu ( t ∗ ) . (4.5)Therefore, when Case 1 occurs, ¯ M hm,l ∗ M δtm,l ∗ should be inserted between M m,l ∗ +1 and M m,l ∗ for the time10nterval [ t ∗ , t ∗ + h ] in Eq. (3.10). The expressions for crossing events from H >
H < f are reversed. Let δt be the first crossing time for Case 2, which can be calculated based on Eq. (4.1). We define δt ∗ as the time where H is maximal, such that H max := H ( u ( t ∗ + δt + δt ∗ )) = max t ∈ [ t ∗ ,t ∗ + h ] H ( u ( t ) , e ), and δ ¯ t as the time where H changes sign back, such that H cr, := H ( u ( t ∗ + δt g )) = 0, where δt g := δt + δt ∗ + δ ¯ t .The estimate of δt follows Eq. (4.1). From a computational point of view, Case 2 can be triggered eitherby (i) H < H > ddt H > ddt H > < δt g < h , or (ii) H > H < ddt H < ddt H < < δt g < h .Since ddt H ( t ∗ + δt + t ) | t = δt ∗ ≈ ddu H cr, [ ˙ u ( t ∗ + δt ) + ¨ u ( t ∗ + δt ) δt ∗ ] + d du H cr, [ ˙ u ( t ∗ + δt ) δt ∗ ] = 0 , we have δt ∗ = − ddu H cr, [ ˙ u ( t ∗ + δt )] ddu H cr, [¨ u ( t ∗ + δt )] + d du H cr, [ ˙ u ( t ∗ + δt )] . (4.6)For δ ¯ t we have H cr, ≈ H max + ddu H max [ ˙ u ( t ∗ + δt + δt ∗ )] δ ¯ t ≈ H cros , + ddu H cros , [ ˙ u ( t ∗ + δt )] δt ∗ + (cid:2) ddu H cros , + d du H cros , [ ˙ u ( t ∗ + δt ) δt ∗ ] (cid:3)(cid:2) ˙ u ( t ∗ + δt ) + ¨ u ( t ∗ + δt ) δt ∗ (cid:3) δ ¯ t = 0 , which gives δ ¯ t = − (cid:2) H cr, + ddu H cr, [ ˙ u ( t ∗ + δt )] δt ∗ (cid:3)(cid:2) d du H cr, + d du H cr, [ ˙ u ( t ∗ + δt ) δt ∗ ] (cid:3)(cid:2) ˙ u ( t ∗ + δt ) + ¨ u ( t ∗ + δt ) δt ∗ (cid:3) − . (4.7)Therefore, for the step from t ∗ to t ∗ + δt , the variational equation can be written as δu N ( t ∗ + δt )... δu ( t ∗ + δt ) δu ( t ∗ + δt ) = M δtm,l ∗ δu N ( t ∗ )... δu ( t ∗ ) δu ( t ∗ ) . (4.8)For the step from t ∗ + δt to t ∗ + δt graz ] we have δu N ( t ∗ + δt graz )... δu ( t ∗ + δt graz ) δu ( t ∗ + δt graz ) = M δt graz m,l ∗ δu N ( t ∗ + δt )... δu ( t ∗ + δt ) δu ( t ∗ + δt ) , (4.9)11here M δt graz m,l ∗ := ˆ M δt graz m,l ∗ ˜ M δt graz m,l ∗ , ˆ M δt graz m,l ∗ := I · · · · · · I − δt ∗ + ¯ δt B δt graz m,l ∗ · · · I − δt ∗ + ¯ δt A δt graz m,l ∗ − , ˜ M δt graz m,l ∗ := I · · · · · · I δt ∗ + ¯ δt B δtm,l ∗ · · · I + δt ∗ + ¯ δt A δtm,l ∗ ,A δt graz m,l ∗ = ∂f j ( u ( t ) ,u N ( t )) ∂u | t = t − m + h ( l ∗ − δt graz and B δt graz m,l ∗ = ∂f j ( u ( t ) ,u N ( t )) ∂u N | t = t − m + h ( l ∗ − δt graz . For the pe-riod [ t ∗ + δt graz , t ∗ + h ], δu N ( t ∗ + h )... δu ( t ∗ + h ) δu ( t ∗ + h ) = ¯ M hm,l ∗ δu N ( t ∗ + δt graz )... δu ( t ∗ + δt graz ) δu ( t ∗ + δt graz ) , (4.10)where ¯ M hm,l ∗ = ˆ M hm,l ∗ ˜ M hm,l ∗ , ˆ M hm,l ∗ := I · · · · · · I − h − δt graz B m,l ∗ +1 · · · I − h − δt graz A m,l ∗ +1 − and ˜ M hm,l ∗ := I · · · · · · I h − δt graz B δt graz m,l ∗ · · · I + h − δt graz A δt graz m,l ∗ . δu N ( t ∗ + h )... δu ( t ∗ + h ) δu ( t ∗ + h ) = ¯ M hm,l ∗ M δt graz m,l ∗ M δtm,l ∗ δu N ( t ∗ )... δu ( t ∗ ) δu ( t ∗ ) , (4.11)Thus, once Case 2 is encountered, ¯ M hm,l ∗ M δt graz m,l ∗ M δtm,l ∗ should be inserted between M m,l ∗ +1 and M m,l ∗ inEq. (3.10) for the step from t ∗ to t ∗ + h .From the discussion above, we can obtain an accurate Jacobian matrix for the Poincar´e map (2.5). Inthe next section, we will discuss the convergence of eigenvalues of the Jacobian matrix when a perturbationis introduced in order to ensure the accuracy of our proposed method.
5. Convergence analysis
According to [35, 52], the spectrum of the Jacobian for the Poincar´e map consists of eigenvalues and0. So we will study the Poincar´e map of Eq. (3.7) and its relevant Jacobian.For the space C d , assume P := [ t , t + ∆ T ], which is an bounded interval of R and ∆ T < + ∞ . C ( P , C d ) denotes the Banach space with all bounded continuous functions from P to C d with the norm || u || C = max t ∈ P | u ( t ) | , where u ∈ C ( P , C d ) and | · | is a given norm on C d .Now, we rewrite Eq. (3.7) as dd t δu ( t ) = F ( t, δu ( t ) , δu N ( t )) , where t ∈ P and F : P × C d × C d → C d ,δu ( t ) = φ δ ( t ) , where t ∈ [ t − τ d , t ] and φ δ ∈ C ([ t − τ d , t ] , C d ) , (5.1)where φ δ is defined in Eq. (3.7). Here, we assume δu d ( t ) = δu N ( t ), and F can be written as F ( t, δu ( t ) , δu d ( t )) = F j, ( t ) δu ( t ) + F j, ( t ) δu d ( t ) , (5.2)where j = 1 , if H ( u ( t ) , e ) > ,j = 2 , if H ( u ( t ) , e ) < ,F ( t − , δu ( t − ) , δu d ( t − )) = F ( t + , δu ( t + ) , δu d ( t + )) , if H ( u ( t ) , e ) = 0 ,F j, ( t ) := ∂f j ( t,u ( t ) ,u d ( t )) ∂u , and F j, ( t ) := ∂f j ( t,u ( t ) ,u d ( t )) ∂u d .According to [34], nonautonomous delayed dynamical system can be represented as an evolutionoperator. So, for any t ∈ P and sufficiently small h >
0, we have U ( t + h, t ) φ δ = δu ( t + h ) , (5.3)where δu ( t + h ) is the solution of Eq. (5.1) at t = t + h . For any time t = t + N t h, ∀ N t ∈ Z + , δu ( t )can be written as δu ( t ) = U ( t + hN t , t + h ( N t − · · · U ( t + 2 h, t + h ) U ( t + h, t ) φ δ . U ( t + h, t ). In order to simplify our discussion, we define the following spaces P := C ([ t − τ d , t ] , C d ) , and P + := C ([ t , t + h ] , C d ) , their relevant norms || · || := max t ∈ [ t − τ d ,t ] | · | , and || · || + := max t ∈ [ t ,t + h ] | · | , and the space P ∗ := C ([ t − τ d , t + h ] , C d ) , with the map L : P × P + → P ∗ satisfying L ( φ δ , z )( η ) = φ δ ( t ) + R ηt z ( θ ) dθ, if η ∈ [ t , t + h ] ,φ δ ( η ) , if η ∈ [ t − τ d , t ] . According to [34], the map L can be divided into two operators L : P → P ∗ and L : P + → P ∗ with L ( φ δ , ω ) = L φ δ + L ω, (5.4)where ( φ δ , ω ) ∈ P × P + , L φ δ = L ( φ δ ,
0) and L ω = L (0 , ω ).In addition, we define the linear operator Θ : P ∗ → P + via[Θ v ]( t ) = F ( t, v ( t ) , v d ( t )) , (5.5)where v ∈ P ∗ , t ∈ [ t , t + h ] and v d ( t ) = v ( t − τ d ). The fixed point problem ω ∗ = Θ L ( φ δ , ω ∗ ) . (5.6)has a fixed point ω ∗ ∈ P + if the original problem (5.1) has a solution in [ t , t + h ]. So ω ∗ satisfies U ( t + h, t ) φ δ = L ( φ δ , ω ∗ ) . (5.7)According to Eq. (5.4), Eq. (5.6) can be rewritten as( I P + − Θ L ) ω ∗ = Θ L φ δ (5.8)where I P + is the identity operator for the space P + . Therefore, we can derive the following propertiesfor the operators Θ L and Θ L . Proposition 1.
If the operator Θ is defined by Eq. (5.5), it is a bounded linear operator with v ∈ P ∗ . Proposition 2. If L and L are defined by Eq. (5.4), then Θ L : P → P + and Θ L : P + → P + are14ounded linear operators with regard to ω ∈ P + . Since system (5.1) can be approximated by large finite ODE systems, the approximated operators areconstructed through discretisation by introducing the relevant discrete space of P and P + along with thefollowing operators. As large finite ODE systems can be obtained from the modified Euler integration,we can adopt linear interpolation to discretise the space P and P + .First of all, based on the time step h , consider the mesh Λ N +1 := ( t − N h, · · · , t − h, t ) in [ t − τ d , t ].We construct a restriction operator r h : P → P N +1 := C d ( N +1) on Λ N +1 , such that r h φ δ ∈ P N +1 ,where [ r h φ δ ] i = φ δ ( t − ( N + 1 − i ) h ) ∈ C d . In addition, there exists a prolongation operator onthe mesh Λ N +1 such that for any ̟ N +1 := ( ̟ T ( t − N h ) , · · · , ̟ T ( t )) T ∈ P N +1 , where ̟ ∈ P ,¯ r h : t ∈ [ t − τ d , t ] → ¯ r h ( t ) ∈ C × d ( N +1) , ¯ r h ( t − ( N +1 − i ) h ) ̟ N +1 = ̟ ( t − ( N +1 − i ) h ) , i ∈ Z [1 , N +1] , and ¯ r h ( t ) ̟ N +1 is a polynomial with a degree less than or equal to 2.Similarly, consider the mesh Λ K +1 := ( t , t + h s , · · · , t + Kh s ) in [ t , t + h ], where 0 < h s < h , K = h/h s , the space P + can be discretised by the restriction operator R h s : P + → P + K +1 := C d ( K +1) on themesh Λ K +1 such that R h s ψ ∈ P + K +1 , where R h s ψ i = ψ ( t + ( i − h s ) ∈ C d . For mesh Λ K +1 we constructa relevant prolongation operator as follows. For any ̟ K +1 := ( ̟ T ( t ) , · · · , ̟ T ( t + Kh s )) ∈ P + K +1 , where ̟ ∈ P + , ¯ R h s : t ∈ [ t , t + h ] → ¯ R h s ( t ) ∈ C d ( K +1) , such that ¯ R h s ( t + ( i − h s ) ̟ K +1 = ̟ ( t + ( i − h s ), i ∈ Z [1 , K + 1], and ¯ R h s ( t ) ̟ K +1 is a polynomial with degree less than or equal to K + 1. Here, theoperator L := ¯ R h s ( t ) R h s is a Lagrange operator [53].Let K = 1 ( i.e. h s = h ) and for any given N , the relevant approximated operator U N +1 , ( t + h, t ) : P N +1 → P N +1 satisfies U N +1 , ( t + h, t )Φ = r h L (¯ r h ( t − τ d )Φ , ¯ R h s ( t )Ψ ∗ ) , (5.9)where t ∈ [ t , t + h ], Φ ∈ P N +1 and Ψ ∗ ∈ P + K +1 , which is the solution of the following equationΨ ∗ = R h s Θ L (¯ r h ( t − τ d )Φ , ¯ R h s ( t )Ψ ∗ ) . (5.10)It is worth noting that the operator ¯ R h s at the time interval [ t , t + h ] can be more accurate if the timestep h is reduced. In this section, we will present the convergence analysis for 0 < τ d < T only. The proof for τ d ≥ T is similar, so will be omitted here. In order to ensure a unique solution for the initial problem (5.1), weintroduce the subspace P + Lip of P + with the norm || ψ || + Lip = l ( ψ ) + || ψ || + , ψ ∈ P + Lip , where l ( ψ ) is the Lipschitz constant of ψ , and the subspace P Lip of P with the norm as || ψ || Lip = l ( ψ ) + || ψ || , ψ ∈ P Lip . To carry out convergence analysis for the eigenvalues of Jacobian of the Poincar´e map (2.5), thefollowing lemmas are given based on [34]. 15 emma 5.1.
For any σ ∗ , σ ∗ ∈ P + , σ ∗ = L Θ L ( φ δ , σ ∗ ) , φ δ ∈ P , (5.11)and σ ∗ = Θ L ( φ δ , σ ∗ ) , φ δ ∈ P , for sufficiently small h , and we have || σ ∗ − σ ∗ || + ≤ c h , (5.12)where c is a positive constant.Based on Eq. (5.9), a new operator in the interval [ t , t + h ] can be introduced as¯ U N +1 , ( t + h, t ) = ¯ r h U N +1 , ( t + h, t ) r h : P → P , (5.13)which has the same geometric and partial multiplicities as the operator U N +1 , ( t + h, t ) in Eq. (5.9).Therefore, there exists a map ¯ U ( t + h, t ) : P → P such that¯ U ( t + h, t ) φ δ = L ( φ δ , σ ∗ ) , φ δ ∈ P , (5.14)where σ ∗ ∈ P + is the solution of Eq. (5.11), and ¯ U N +1 , ( t + h, t ) can be written as¯ U N +1 , ( t + h, t ) = L ¯ U ( t + h, t ) L . Lemma 5.2.
If the operator ¯ U ( t + h, t ) is defined as Eq. (5.14), we have || ¯ U ( t + h, t ) − U ( t + h, t ) || ≤ c h , (5.15)where c is a positive constant.It is worth noting that the evolution operator ¯ U ( t + ih, t + ( i − h ), where i = 1 , · · · , ¯ N and¯ N := N + n + △ N + 1, must have the same properties as the operator U ( t + ih, t + ( i − h ) in theinequality (5.15). Thus, the Poincar´e map can be obtained by combining all the evolution operators U ( t + ih, t + ( i − h ) over the entire time interval [ t , t + T ]. As a result, the convergence problem isequivalent to studying the convergence of the operator Q ¯ Ni =1 ¯ U ( t + ih, t + ( i − h ) to U ( t , t + T ). Lemma 5.3.
For the entire interval [ t , t + T ] and a sufficiently small time step h , we can obtain || U ( t + T, t ) − ¯ N Y i =1 ¯ U ( t + ih, t + ( i − h ) || ≤ c h , (5.16)where i = 1 , , · · · , ¯ N , ¯ N := N + n + △ N + 1, and c is a positive constant.Combining the inequality (5.16) with the results in [35, 52] and Theorem 4.6 and 4.7 in [34], thefollowing lemma can be obtained. Lemma 5.4.
Let λ ∈ C \ { } be an isolated eigenvalue for the operator U ( t + T, t ) with the finitealgebraic multiplicity m a and ascent (length of longest Jordan chain) κ , and Γ be a neighborhood of λ λ of U on the time interval [ t , t + T ]. For a sufficiently small h , ¯ U ( t + T, t ) has m eigenvalues λ ,ι , where ι = 1 , . . . . , m a , and we have max ι =1 ,...,m a | λ − λ ,ι | ≤ c h κ , (5.17)where c is a positive constant.It should be noted that ¯ U N +1 , and ¯ U have the same nonzero eigenvalues, geometric and partialmultiplicities and eigenvectors. This leads to the following theorem. Theorem 5.5.
Let λ ∈ C \ { } be an isolated eigenvalue for the operator U ( t + T, t ) with the finitealgebraic multiplicity m a and the ascent κ , and let Γ be a neighborhood of λ for the time interval[ t , t + T ]. For a sufficiently small h , ¯ U N +1 , ( t + T, t ) has m eigenvalues λ N +1 , ,ι , where ι = 1 , . . . . , m a and we have max ι =1 ,...,m a | λ − λ N +1 , ,ι | ≤ c h κ , (5.18)where c is a positive constant.The inequality (5.18) holds for any interval [ t m , t m + T ]. From the above study, we can ensure thatour proposed approximation method has the expected convergence rate on the nonzero characteristicmultipliers of the system (5.1). So our approximation for the Jacobian of the Poincar´e map (2.5) isreliable. It is also worth noting that by adopting a high-order integration method (e.g. Runge-Kuttamethod) with a sufficiently small time step h , the approximated operator could be more accurate O ( h ).However, this would also require higher-order corrections at the crossing and grazing events for the termsderived in Section 4. Without these corrections the convergence of the approximated operator cannotbe guaranteed as the same with the order of the numerical integration. Furthermore, if the systemencounters sufficiently large number grazing events, the convergence rate will be lower than O ( h ) dueto these grazing events.
6. Calculation of the Lyapunov exponents
The dynamics of system (2.4) can be represented by the Poincar´e map (2.5) as Y m +1 , = P m ( Y , ) = P ◦ · · · P ◦ P ( Y , ) , (6.1)where the Jaobian matrix of P m is Q mi =1 M i . According to Definition 2.1, LEs can be calculated as ϑ i = lim m →∞ m ln | λ mi | , i = 1 , · · · , d ( N + 1) , (6.2)where λ mi is the i th eigenvalues of Q mi =1 M i .However, calculating LEs by using Eq. (6.2) will introduce an overflow problem. Specifically, someelements of the Jacobian matrix will be very large for chaotic attractors, and some of them could be verysmall for periodic attractors, which may cause inaccuracies. On the other hand, calculating LEs from theJacobian matrix directly is time-consuming as the time-delayed dynamical system is high-dimensional.To overcome these issues , LEs can be computed according to the average exponential divergence ratebetween the basis orbit started from Y (0) and its neighborhood orbit along the direction of v , = Y , || Y , || ϑ ( Y , , v , ) = lim m →∞ m ln || δY m, |||| δY , || , (6.3)where || δY m, || is the norm of δY m, and m ∈ Z + .Next, choose Y , ∈ R d ( N +1) , and its related linearly independent initial perturbed vector ( δY , , δY , , · · · , δY d ( N +1)1 , ) can be normalised as( δv , , δv , , · · · , δv d ( N +1)1 , ) = ( δY , || δY , || , δY , || δY , || , · · · , δY d ( N +1)1 , || δY d ( N +1)1 , || ) . (6.4)Substituting the vector (6.4) to Eq. (6.1) obtains the second vector ( δY , , δY , , · · · , δY d ( N +1)2 , ), andGram-Schmidt orthonormalization [50] can be applied to normalise the second vector, which gives a newvector ( δv , , δv , , · · · , δv d ( N +1)2 , ). For the next iteration, the second vector will be used as the initialvector to be substituted into Eq. (6.1). Likewise, repeating m times for this process gives the m th vector( δY m, , δY m, , · · · , δY d ( N +1) m, ). The steps of Gram-Schmidt orthonormalization are given as follows V m, = δY m, ,δv m, = V m, || V m, || ,V m, = δY m, − < δY m, , δv m, > δv m, ,δv m, = V K (0) || V m, || , ... V d ( N +1) m, = δY N +1) m, − < δY N +1) m, , δv m, > δv m, − · · ·− < δY d ( N +1) m, , δv d ( N +1) − m, > δv d ( N +1) − m, ,δv d ( N +1) m, = V d ( N +1) m, || V d ( N +1) m, || , where || V im, || is the norm of V im, , h δY im, , δv ¯ im, i ( i, ¯ i = 1 , , · · · , d ( N + 1)) is a standard scalar product.Finally, LEs can be calculated by using ϑ i ≈ m ln m Y ̺ =1 || V i̺ (0) || = 1 m m X ̺ =1 ln || V i̺ (0) || . (6.5) Remark.
Based on the above analysis, a guideline for the implementation of the algorithm is presentedas follows.
Step 1 : Calculate the Jacobian matrix according to the relevant trajectory at the time step after thesystem is stabilised by the time-delayed feedback controller;
Step 2 : If the trajectory approaches to grazing, calculate its relevant Jacobian using Eq. (4.5) orEq. (4.11), and then insert it to the matrix M m in Eq. (3.11) at the grazing moment; Step 3 : Choose appropriate initial perturbed unit vectors, and calculate the Floquet Multipliers ofeach Poincar´e map using Gram-Schmidt orthonormalization;18 tep 4 : Calculate the LEs using Eq. (6.5) after several evolutions of Poincar´e map.
7. Numerical studies
In this section, we will show the effectiveness of our proposed method by studying the soft impactingsystem with a delayed feedback controller presented in Fig. 2.1. Since the system has many coexistingattractors when grazing is encountered [40], our control objective here is to drive the system from itscurrent attractor to a desired one. Calculating the LEs of the system allows us to monitor the stabilityof the delayed feedback controller and its effective parametric regime.We choose the following parameters for the impacting system, ζ = 0 . , e = 1 . , a = 0 . , β = 28 and ω = 0 . . For these parameters a grazing event is encountered, and a chaotic and a period-5 attractors coexist asshown in Fig. 7.1.
Figure 7.1: Basin of attraction of the impacting system computed for ζ = 0 . , e = 1 . , a = 0 . , β = 28 and ω = 0 . τ d ≥ T Fig. 7.2 presents the first example of using the delayed feedback controller (2.3) for which a largedelayed time (i.e. τ d ≥ T ) was considered, and the control parameter k was varied from 0 to 1 .
4. As canbe seen from Fig. 7.2(a), the largest LEs are all greater than 0 for k ∈ [0 , .
04] and the system presentsa chaotic motion as shown in Fig. 7.2(b). The phase trajectory of the chaotic motion for k = 0 .
02 ispresented in Fig. 7.2(c). For k ∈ (0 . , . k = 0 .
055 indicating a period doubling of the system. Similarly, at k = 0 . k = 0 .
07. For k ∈ [0 . , . igure 7.2: (a) LEs and (b) displacement of the impacting system under the delayed feedback controller as functions of thecontrol parameter k . Black, red and green lines denote the two largest LEs and the zero line, respectively. Additional panelsshow the phase trajectories of the system calculated for (c) k = 0 .
02, (d) k = 0 .
55 and (e) k = 1. Black dots represent thePoincar´e sections, and blue lines represent the impact boundary. A critical issue for computing nonsmooth dynamical systems is that the accumulated computationalerror from the impact boundary due to grazing event could lead to inaccurate simulation. Fig. 7.3compares the computations of the impacting system for e = 1 . t = 9722, and thephase trajectories from chaotic (grey line) to period-1 (red line) response are shown in Fig. 7.3(b). It canbe seen from the figure that the accumulated error was built up in the number of impacts, and a cleardifference can be observed from t = 10411. The cause of such a difference can be found from Figs. 7.3(c)and (d), where the time histories of displacement of the impacting system are shown. As can be seenfrom the figures, the system with the grazing estimation algorithm was stabilised quicker than the onewithout the algorithm. < τ d < T For the case of a small time delay (i.e. 0 < τ d < T ), we present the example for τ d = T / k ∈ [0 , . k ∈ (0 . , . k = 0 . k = 0 . k ∈ [0 . , . k ∈ [0 . , . igure 7.3: (a) Number of impacts as a function of time without (black line) and with (orange line) the grazing estimationalgorithm based on the discontinuous condition calculated for ζ = 0 . e = 1 . a = 0 . β = 28, ω = 0 .
802 and k = 1 . k = 0 . k increases. Tocompare Figs. 7.4(a) and (b), the evolution of the calculated LEs is consistent with system’s bifurcation,which is also demonstrated by the phase trajectories presented in Figs. 7.4(c)-(f).
8. Conclusions
This paper studies a numerical method for calculating the LEs of time-delayed piecewise-smoothsystems by using a soft impacting system under the delayed feedback control with a particular focus onits near-grazing dynamics. The main feature of the proposed method is that it can provide improvedaccuracy for the stability analysis of periodic orbits by estimating the point of discontinuity locally alongtrajectories of piecewise-smooth DDEs with an accuracy of the same order as its integration method.In addition, the method can also be applied to the other nonsmooth dynamical systems with a delayedargument, such that it can be used as a generic computational tool for stability analysis.The main tasks were to build an effective variational equation and obtain the Jacobian for the delayedimpacting system. As the delayed impacting system is infinite dimensional, it was approximated by finitedimensional systems, which were discretised by the modified Euler integration method at each timestep. Then the DDE system converted to a time-discrete map by constructing a Poincar´e map, and itslinearisation was introduced to obtain its variational equation. Then the Jacobian of the map was obtainedby combining all the approximating systems linearised from the variational equation at each time step inone period of external excitation. In order to increase the convergence rate and improve computationalaccuracy, a grazing estimation algorithm was introduced. The convergence rate of eigenvalues of theJacobian matrix was studied by using the spectral theory of the evolutionary operator. In particular, thedelayed impacting system was described as an evolutionary operator with the expected convergence rate21 igure 7.4: (a) LEs and (b) displacement of the impacting system under the delayed feedback controller as functions of thecontrol parameter k . Black, red and green lines denote the two largest LEs and the zero line, respectively. Phase trajectoriesof the system calculated for (c) k = 0 .
01, (d) k = 0 .
03, (e) k = 0 .
043 and (f) k = 0 .
052 are shown. Black dots represent thePoincar´e sections, and blue lines indicate the nonsmooth boundary. for the relevant nonzero eigenvalues of the Jacobian, therefore guaranteeing the reliability of the proposednumerical method.Our numerical studies considered two scenarios of delay time in the system, a larger ( τ d ≥ T ) anda smaller (0 < τ d < T ) delay than the period of excitation. Both cases showed that the calculated LEswere consistent with the bifurcation of the system, and the grazing estimation algorithm had improvedaccuracy for simulating nonsmooth dynamical systems. Acknowledgements
This work has been supported by EPSRC under Grant No. EP/P023983/1. Mr Zhi Zhang would liketo acknowledge the financial support from the University of Exeter for his Exeter International ExcellenceScholarship. JS’ research is supported by EPSRC Fellowship EP/N023544/1 and the European UnionsHorizon 2020 research and innovation programme under grant agreement number 820970, project TiPES.
Appendix.
Proof of Proposition 1 : Let v , v ∈ P ∗ , whereΘ v ( t ) = F ( t, v ( t ) , v ( t − τ d )) , v ( t ) = F ( t, v ( t ) , v ( t − τ d )) . Then we can obtain || Θ( v + v ) || ≤ | F ¯ j, ( t ) v ( t ) + F ¯ j, ( t ) v ( t − τ d ) | + | F j, ( t ) v ( t ) + F j, ( t ) v ( t − τ d ) | = || Θ v || + || Θ v || In addition, according to Eq. (5.2), there must exist a positive constant B Θ satisfying that, for any v ∈ P ∗ , || Θ v || ≤ B Θ || v || . Therefore, the operator Θ is bounded and linear in the space. Proof of Proposition 2 : For ∀ φ δ , there exist ω , ω , ω ∈ P + such that Θ L ( φ δ , ω ) = ω and ω = ω + ω . So we have Θ L ( φ δ , ω ) = Θ[ L ( φ δ , ω ) + L ω ]= Θ L ( φ δ ) + Θ L ω + Θ L ω = Θ L ( φ δ ) + Θ L ( ω + ω ) . According to the Eqs. (5.5)-(5.7), if Θ L (0 , ω ) = ω (where ω ∈ P + ) holds, L (0 , ω ) must be the solution ofthe following system dd t δu ( t ) = F ( t, δu ( t ) , δu d ( t )) ,δu ( s ) = 0 , (8.1)where F ∈ C ( P , R d ) and s ∈ [ t − τ d , t ]. Then for any ω ∈ P + , it gives || Θ L ω || = || ω || , so Θ L isbounded.Let φ δ, , φ δ, ∈ P , and for φ δ, + φ δ, , there exists ω ∈ P + such that Θ L ( φ δ, + φ δ, , ω ) = ω . Also,there exists ω , ω ∈ P + , such that ω = ω + ω . Then we have L ( φ δ, + φ δ, , ω ) = L ( φ δ, , ω ) + L ( φ δ, , ω )= L φ δ, + L ω + L φ δ, + L ω = L ( φ δ, + φ δ, ) + L ω Since Θ[ L φ δ, + L ω + L φ δ, + L ω ] = Θ L φ δ, + Θ L φ δ, + Θ L ω and Θ[ L ( φ δ, + φ δ, ) + L ω ] = Θ L ( φ δ, + φ δ, ) + Θ L ω, Θ L is a bounded linear operator. Proof of Lemma 5.1 : Based on Theorem 3.3 in [34] and let σ ∗ = σ ∗ + ρ ∗ , we have || ρ ∗ || := || σ ∗ − σ ∗ || + = || ( I P + − L Θ L ) − || || ( I P + − L ) || + || σ ∗ || + Lip , (8.2)For sufficiently small h , || ( I P + − L ) || + is the global error from the modified Euler integration, which23atisfies || ( I P + − L ) || + ≤ c h , where c is a positive constant. Since I P + − L Θ L = ( I P + − Θ L ) + ( I P + − L )Θ L , and Θ L is bounded, if h →
0, ( I P + − L Θ L ) − = ( I P + − Θ L ) . In addition, as σ ∗ = ( I P + Lip − Θ L ) − Θ L φ δ , (8.3)and || σ ∗ || + Lip ≤ || ( I P + Lip − Θ L ) − || || Θ L || || φ δ || Lip , (8.4) || σ ∗ || + Lip is bounded. Thus, there must exist a positive constant c for Eq. (8.2) satisfying || ρ ∗ || ≤ c h . Proof of Lemma 5.2 : For ( φ δ , ω ∗ ) , ( φ δ , ω ∗ ) ∈ P + Lip × P + , based on Eq. (5.4), we have || ¯ U ( t + h, t ) − U ( t + h, t ) || = || L ( φ δ , ω ∗ ) − L ( φ δ , ω ∗ ) || = || L ( ω ∗ − ω ∗ ) || , where ω ∗ = ¯ I P + Θ L ( φ δ , ω ∗ )and σ ∗ = Θ L ( φ δ , ω ∗ ) . So || ¯ U ( t + h, t ) − U ( t + h, t ) || = || L ( ω ∗ − ω ∗ ) || = || Z t + ht ( ω ∗ − ω ∗ )( t )d t || = || ( ω ∗ − ω ∗ ) || + h. According to Eqs. (8.2) and (8.3) and the inequality (8.4), we have || ¯ U ( t + h, t ) − U ( t + h, t ) || ≤ c h . Proof of Lemma 5.3 : According to Lemma 5.2, we assume that there are two positive constants M and M such that || ¯ N − Y i =2 U ( t + ih, t + ( i − h ) || ≤ M , (8.5)and || ¯ N − Y j =1 ¯ U ( t + jh, t + ( j − h ) || ≤ M . (8.6)24herefore, || U ( t + T, t ) − ¯ N Y i =1 ¯ U ( t + ih, t + ( i − h ) ||≤ ¯ N ¯ N − Y i =2 U ( t + ih, t + ( i − h ) ¯ N − Y j =1 ¯ U ( t + jh, t + ( j − h ) c h ≤ ¯ N M M c h = Th M M c h = c h . Compliance with ethical standardsConflict of interest
The authors declare that they have no conflict of interest concerning the publication of this manuscript.
Data accessibility
The datasets generated and analysed during the current study are available from the correspondingauthor on reasonable request.