Calculating the Lyapunov exponents of a piecewise smooth soft impacting system with a time-delayed feedback controller

Lyapunov exponent is a widely used tool for studying dynamical systems. When calculating Lyapunov exponents for piecewise smooth systems with time delayed arguments one faces two difficulties: a high dimension of the discretized state space and a lack of continuity of the variational problem. This paper shows how to build a variational equation for the efficient construction of Jacobians along trajectories of the delayed nonsmooth system. Trajectories of a piecewise smooth system may encounter the so-called grazing events, where the trajectory approaches discontinuity surfaces in the state space in a non-transversal manner. For these events we develop a grazing point estimation algorithm to ensure the accuracy of trajectories for the nonlinear and the variational equations. We show that the eigenvalues of the Jacobian matrices computed by the algorithm converge with an order consistent with the order of the numerical integration method. Finally, we demonstrate the proposed method for a periodically forced impacting oscillator under a time-delayed feedback control, which exhibits grazing and crossing of the impact surface.


Introduction
Analysing grazing bifurcation for nonsmooth systems is a challenging task [1]. In particular, for vibro-impact systems, such as ship mooring interactions [2], bearing looseness [3] and the multi-degreeof-freedom impact oscillators [4], during the process of modelling, their motion could be totally different based on the different types of the grazing impacts [5]. Generally, system may have an abundant of coexisting attractors when grazing occurs, and any tiny differences in modelling will lead to different motion of the system [6] . For example, the motion of an impact oscillator experiences significantly change due to a slight variation of its parameter when grazing bifurcation is encountered [7]. In [8], Nordmark studied the characteristic scaling behaviour in grazing bifurcation, and used the self-similarity under scaling to derive a renomalised mapping. Nordmark [9] presented the grazing bifurcation of a hard impact oscillator which had singular Jacobian of Poincaré map by using the first-order Taylor expansion. It has shown that the more accurate grazing bifurcation to be depicted, the more precise the stability of the oscillator to be studied. This paper will study a new method to improve the accuracy of calculation of grazing bifurcation by estimating the impacting moment. Based on this accurate grazing trajectory, stability analysis of the system can be carried out.
In many applications [10][11][12][13][14], delay is always considered in their differential equations in which the derivative of the unknown functions at a certain time depends on the value of the function at previous time, the so-called delay differential equations (DDEs). For example, Zhang et al. [12] studied a delayed pest control model which was a high-dimensional differential equation with impulsive effects at different fixed impulsive moments. In [13], Carvalho and Pinto used a mathematical model with delay to describe the dynamics of AIDS-related cancers with the treatment of HIV and chemotherapy. In [14], Yan et al. studied the basin of a time-delayed system modelling cutting processes to determine the unsafe cutting zone. The above studies concern about smooth DDEs. The study of nonsmooth DDEs is more challenging due to the lack of an accurate algorithm for the computation of grazing events. Until now, there are very few systematic studies regarding to the nonsmooth DDEs, which is the focus of this paper. This work will study a new algorithm for determining grazing events (improving the computational accuracy) and a new method for calculating Lyapunov exponents (LEs) along trajectories of a piecewise smooth DDE.
LEs are an important tool for studying the stability of dynamical system. If the largest LEs is greater than zero, any small perturbation on the initial condition will result in an exponential divergence of the resulting perturbed trajectory until the distance between perturbed and unperturbed trajectory is no longer small. This senstivity with respect to initial conditions is one of the defining features of chaos. Therefore, the development of an efficient method for calculating the LEs of dynamical system is an active area of research, see e.g. [15][16][17][18][19]. On the one hand, Parker and Chua [15] introduced the systematic methods for estimating LEs for the smooth dynamical systems given in the form of ordinary differential equations or maps. Wolf et al. [16] developed a method for extracting the largest LEs from an experimental time series. For piecewise smooth systems, Jin et al. [19] studied the Poincaré map method for calculating the LEs of a hard impact oscillator, which was constructed using means of local maps to avoid calculating the Jacobian matrices at the points of discontinuity. Lamba and Budd [20] has shown that the largest LE has a discontinuous jump at grazing bifurcation in Filippov systems and scales like 1/| ln ǫ|, where ǫ is the bifurcation parameter. Since DDE is an infinite dimensional system, calculation of LEs for nonsmooth DDEs becomes more complicated than the smooth one. In princple, DDEs could be approximated by ordinary differential equations (ODEs) which can be linearised by numerical integration [21,22], and then the Poincaré map can be constructed for calculating LEs. Studies by Repin [23] and Györi and Turi [24] have shown that DDEs can be analysed using approximating high-dimensional ODEs. However, if delay time is large, in addition, compared with smooth DDEs [25][26][27], calculating the LEs of nonsmooth DDEs needs to store excessive history data points during delay period, such as the datas at the discontinuous moments, and then the global convergence of the system can not be guaranteed. Therefore, this may cause inaccuracy in calculating the eigenvalues of the Jacobian matrix which are used for estimating the LEs of nonsmooth DDEs.
The contribution of this work is the development of a novel method for calculating the LEs of piecewise smooth differential equations with a delay argument, which can provide an improved accuracy in stability analysis of periodic oirbits. If the algorithms do not estimate the points of discontinuity along the trajectory with an accuracy of the same order as the integration method, especially in the grazing event, the expected discontinuous of the coefficients of the variational problem have unexpectedly low accuracy leading to an accumulation of errors. The novelty of our method is that it estimates the points of discontinuity locally along trajectories of the piecewise smooth DDEs, improving the accuracy of computations of system trajectory and LEs. Our proposed method can also be extended to other nonsmooth dynamical systems, such as the hard impact oscillator system with a delay control and stickslip vibrations with a delay term. To demonstrate the reliability of the LEs, we will carry out an error analysis for the nonzero eigenvalues of the Jacobian by comparing it to spectral approximation methods introduced by Chatelin [28] and Breda et al. [25,27] in the special cases where the latter are applicable. Our study indicates that the proposed method can reduce error for the nonzero eigenvalues of the Jacobian by increasing the approximated dimensions of the approximated system (ODEs) slightly, which is generated from linearising the DDEs by the numerical integration.
The rest of this paper is organised as follows. Section 2 introduces a practically relevant example of the type of DDEs the algorithm could be applied to. The example is a mathematical model of a periodically forced mechanical oscillator with a soft impact nonlinearity, which corresponds to a one-sided linear spring. This is followed by some basic relevant definitions and preparations. Section 3 presents the method for constructing the Jacobian of a Poincaré map of a piecewise smooth DDEs. Section 4 describes the estimation method, the so-called grazing estimation algorithm, for determining points of discontinuity accurately. Two cases of grazing events are considered depending on the geometry of the trajectory. Section 5 uses linear operator theory to carry out an error analysis for the eigenvalues of the Jacobian, which can validate the reliability of our proposed method. Section 6 presents the details for computing LEs. Examples and several control scenarios of the system introduced in Section 2 are presented to demonstrate the accuracy of the method in Section 7. Finally, conclusions are drawn in Section 8.

Mathematical model and related preparations
The impact oscillator shown in Fig 2.1 represents a mechanical system encountering intermittent socalled soft impacts, which will be studied in the present work. The nondimensional equations of motion of the impact oscillator can be written in a compact form as follows [7,29]: where H(·) stands for the Heaviside step function and x ′ , v ′ denote differentiation with respect to the nondimensional time τ . The discontinuity boundary is fixed at x = e, with e > 0 being the nondimensional gap to the rest point of linear spring. Model (2.1) is derived from the representation in Fig. 2.1 by nondimensionalizing the variables and parameters of the system are nondimensionalized according to: where y 0 > 0 is an arbitrary reference distance, ω n is the natural angular frequency, ω is the frequency ratio, β is the stiffness ratio, ζ is the damping ratio, a is the nondimensionalized forcing amplitude, and e is the nondimensionalized gap between the mass and the secondary spring. In the present work, we will consider a control signal u(τ ), τ ≥ 0, which will be superimposed on the system's external excitation as follows where defines the propotional feedback controller that feedbacks the difference between current measurement of v and a measurement of v from some time τ d ago [30]. In the expression above, k ≥ 0 represents the feedback gain of the controller and τ d > 0 stands for a predefined time delay.
Example (2.3) can be rewritten in the form of a general piecewise continuous DDE with a periodic external excitation as H(y(t), e) = 0, is sufficiently smooth function and p : R + → R d is smooth and periodic with the period T > 0 and τ d > 0. In the present work, we only consider one single delay in the system for simplicity, and assume that for any y,ȳ, y d ,ȳ d ∈ R d , it satisfies where y d = y(t − τ d ) and l 1 , l 2 ≥ 0, and the for any y,ȳ ∈ R d , it has |H(y, e) − H(ȳ, e)| ≤ l 3 |y −ȳ|, where l 3 ≥ 0 and | · | is a given norm on R d . We assume that the initial condition is from a suitable initial function on [t 0 , t 0 − τ d ]. The general form (2.5) belongs to the class of hybrid dynamical systems [1], which consist of a flow (in our case only forward in time), combined with discrete events. Take N ∈ Z + sufficiently large. Define the grid points τ for all t ≥ 0, i = 0, . . . , N , the system (2.5) can be approximated by the d(N + 1) dimensional piecewise smooth abstract Cauchy problem studied in [23]. The details of this approximation approach will be presented in Section 3. This approximation method has aslo been studied by Krasovskii [31], and particularly, the solution of the approximated system uniformly converged to its original DDEs when N → ∞. By using the same approach, Györi and Turi [24] and Banks [32] carried out the convergence analyses for two differential DDEs. Breda et al. [26] studied the characteristic roots of DDEs, and use the Runge-Kutta method to construct a high-dimensional approximating system [25]. The nonzero eigenvalues of evolution operators were computed through a pseudospectral collection, which were used to analyse the asymptotic stability of DDEs. Since system (2.5) is not a smooth DDE, such that its trajectories will encounter discontinuities repeatedly, the methods used for the smooth DDEs are not available for the system (2.5). Our plan is to derive a Poincaré map for discretized (motivated by the periodic forcing of system 2.5) system and study its stability by considering the Jacobian matrix for the map. After reduction to a Poincaré map, we may define LEs for this time-discrete map.
For the piecewise DDE (2.5), we consider a constant phase surface as the Poincaré section defined by P T s := {(y, t) ∈ R d × R + | t = t 0 + kT, k ∈ Z + }, , and the relevant Poincaré map is The LEs of the Poincaré map can be defined as follows.
Definition 2.1. [15] For any initial condition x 0 , let {x m } ∞ m=0 be the corresponding orbit of a ndimensional, discrete-time system P , and let λ whenever the limit exists.

Constructing the Jacobian Matrix of the Poincaré map
For the nonsmooth system with a delay τ d smaller than the forcing period T (0 < τ d < T ) , the period T can be written as T = nτ d + ∆t, for n ∈ Z + and ∆t ∈ [0, τ d ). For any time interval [t m , t m + τ d ], where t m = t 1 + (m − 1)T , t 1 = t 0 and m ∈ Z + , the solution of system (2.5) can be approximated by N samples at the step h = τ d N by using numerical integration. For example, by using the Modified Euler integration [33], it gives After N + 1 iterations, a local map for the delayed time interval [t m , t m + τ d ] can be defined as Repeating n times of this map and combining them together, we can obtain the following.
Once an arbitrary perturbation δU is applied, the variational equation can be written as In fact, this Eq. (3.6) can be obtained from discretising the continuous variational equation of the system (2.5), and its form can be obtained as with a suitable initial function φ δ satisfying that there exists a sufficiently small ǫ, such that φ δ (t 1 ) = Linearising the equation (3.7) in the interval [t m , t m + nτ d ] by using the modified Euler integration, it gives . . .
By using the map(3.2), the matrix form of the variational equation (3.9) can be rewritten as Since we have n maps, combining all the maps for the interval [t m , t m + T ] gives In addition, the map P △t for the interval [t m + nτ d , t m + T ] can be written as Finally, the variational equation can be obtained as Similarly, for the common system with a large delay time as τ d ≥ T , the solution of system (2.5) can be approximated by N samples at the step h = τ d N by using numerical integration, which can be consider as a special case of the nonsmooth system with a small delay time (0 < τ d < T ) when n = 0. Let N T = T h be the sample number for one period T , construct the map P d , and combine all the linearised maps at the interval [t m , t m + T ]. Finally, we can obtain the same variational equation as Eq. (3.11) and the Jaocbian matrix of the Poincaré map P .

Modifying the algorithm on the discontinuous condition
In this section, we will discuss a special phenomenon of the impact oscillator, the so-called grazing event. Since the system has rich complex dynamics when it experiences grazing [5,34], a careful consideration in calculating this discontinuous moment is required. Another reason is that the global error of our proposed algorithm will depend on the accuracy of the switching function, as the error in the switching boundary may be accumulated leading to undesired large global error. Therefore, during the grazing event, we need to modify our proposed algorithm in the last section by considering the following two grazing cases as illustrated in Fig. 4
. It is worth noting that δu i (t * + δt) can be approximated through linear interpolation, based on the historical datas of the delay time interval, which includes the grazing datas.
Similarly, for the time interval [t * + δt, t * + h], we can obtain

Case 2
If δt is the crossing moment for Case 2, which can be calculated based on Eq. (4.1), δt * is defined as the moment satisfying H max := H(u(t * + δt + δt * )) = max t∈[t * ,t * +h] H(u(t), e), and δt is the moment which satisfies H cr,2 := H(u(t * + δt g )) = 0, where δt g := δt + δt * + δt. The estimation of δt is similar with equation (4.1). From a computational point of view, Case 2 can be triggered either by (1) . (4.6) We also have (4.7) Therefore, for the time interval [t * , t * + δt], its variational equaiton can be written as, For the period [t * + δt, t * + δt graz ], . . . From the discussion above, we can obtain the Jacobian matrix for the Poincaré map (2.6) of the piecewise DDEs (2.5). In the next section, we will discuss the convergence of the eigenvalues of the Jacobian matrix when a perturbation is applied in order to ensure the accuracy of our proposed method.

Properties of the evaluation operator
According to [28,35], the eigenvalues of the Jacobian matrix for the Poincaré map can be considered as the spectral of operator. So, we will study the Poincaré map of the variational Eq. (3.7) and its relevant operator next. For the space C d , assume P := [t 1 , t 1 + ∆T ], which is an bounded interval of R and ∆T < +∞, and 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 can rewrite Eq. (3.7) as where φ δ is defined in Eq. (3.7). Here, we assume δu d (t) = δu N (t), and F can be written as where F j,1 (t) := ∂fj (t,u0(t),u d (t)) ∂u0 and F j,2 (t) := ∂fj (t,u0(t),u d (t)) ∂u d . According to [27], nonautonomous delayed dynamical system can be represented as an evolution operator. So, for any t 1 ∈ P and sufficiently small h > 0, we have where δu 0 (t 1 + h) is the solution of Eq. (5.1) at t 1 + h. For any time t = t 1 + N t h, ∀N t ∈ Z + , δu 0 (t) can be written as Next, we will construct the approximation operator with finite dimensions for the evolution operator U (t 1 + h, t 1 ). In order to simplify our discussion, we define the following operators and | · |, and the space with the map L : P × P + → P * satisfying According to [27], the map L can be divided into two operators L 1 : P → P * and L 2 : P + → P * with where (φ δ , ω) ∈ P × P + , L 1 φ δ = L(φ δ , 0) and L 2 ω = L(0, ω).
In addition, the linear operator Θ : P * → P + satisfies where v, v d ∈ P * and t ∈ [t 1 , t 1 + h], and it has a fixed point ω * ∈ P + which is the solution of the following equation ω * = ΘL(φ δ , ω * ). where I P + is the identity operator for the space P + . Therefore, we can derive the following properties for the operators ΘL 1 and ΘL 2 .

Approximation of the evaluation operator
Since the system (5.1) can be approximated by large finite ODE systems, the approximated operators are constructed thourgh discretization by introducing the relevant discrete space of P and P + along with the following operators. As large finite ODE systems can be obtained from the modified Euler integration, we can adopt linear interpolation to discretize the space P and P + . First of all, based on the time step h, consider the mesh Λ N +1 : In addition, there exists a prolongation operator on the mesh Λ N +1 as, for any is a polynomial with a degree less than or equal to 2.
Let K = 1 ( i.e. h s = h ) and for any given N , the relevant approximated operator U N +1,1 (t 1 + h, t 1 ) : where t ∈ [t 1 , t 1 + h], Φ ∈ P N +1 and Ψ * ∈ P + K+1 , which is the solution of the following equation It is also worth noting that the operatorR hs at the time interval [t 1 , t 1 + h] can be more accurate if the time step h is reduced.

Convergence analysis for the nonzero eigenvalues of the Jacobian Matrix
In this section, we will present the convergence analysis for 0 < τ d < T only, and the proof for τ d ≥ is similar, so will be neglected. In order to ensure the only solution for initial problem (5.1), we introduce the subspace P + Lip of P + with the norm 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 the convergence analysis for the eigenvalues of Jacobian matrix of the Poincaré map (2.6), the following lemmas are given.
Lemma 5.2. If the operatorŪ 2 (t 1 + h, t 1 ) is defined as Eq. (5.18), we have where c 3 is a positive constant.
Proof. According the result of Lemma 5.2, we assume that there are two positive constants M 1 and M 2 such that and ||N −1 Therefore, Combined the inequality (5.20), with results in [28,35] and Theorem 4.6 and 4.7 of [27], the following results can be obtained. It should be noted thatŪ N +1,2 andŪ 2 have the same nonzero eigenvalues, geometric and partial multiplicities and eigenvectors.So, the following theorem can be deduced. For any interval [t m , t m + T ], the result (5.24) always exists. From the above study, we can ensure that our proposed calculation method has a good convergence rate on the nonzero characteristic multipliers of the system (5.1). So, our approximation for the Jacobian matrix of the Poincaré map (2.6) is reliable. It is also worth noting that by adopting a high-order integration method (eg. the Runge-kutta method) with a sufficiently small time step h, the approximated operator could be more accurate O(h 4 ). If the local estimation in Section 4 is not used, the convergence of the approximated operator can not be guaranteed as the same with the order of the numerical integration, and if the system witness sufficiently large number of times of grazing, the convergence rate also can be lower than O(h 2 ), which is resulted from the grazing situation. This is one of the major contributions of this work.

Calculation of the Lyapunov Exponents
Since the entire motion of the system (2.5) can be represented by the Poincaré map (2.6) as Y m+1,0 = P m (Y 1,0 ) = P • · · · P • P (Y 1,0 ), (6.1) where the Jaobian matrix of P m is m i=1 M i . According to Definition 2.1, LEs can be calculated as is the ith eigenvalues of m i=1 M i . However, calculating LEs by using Eq. (6.2) may introduce an overflow problem. In detail, some elements of the Jacobian matrix could be very large for chaotic attractors, while, some of them could be very small for periodic attractors, which causes unreliable calculation. On the other hand, calculating LEs from the Jacobian matrix directly is time-consuming as time-delayed dynamical system is highdimensional. To overcome tthese issues , LEs can be computed according to the average exponent divergence rate between the basis orbit started from Y 1 (0) and its neighborhood orbit along the direction where ||δY m,0 || is the norm of δY m,0 and m ∈ Z + . Next, choose Y 1,0 ∈ R d(N +1) , and its related linearly independent initial perturbed vector (δY 1 1,0 , δY 2 1,0 , · · · , δY ), and the Gram-Schmidt ortho-normalization [33] is applied to normalise the second vectors, which gives the new vector (δv 1 2,0 , δv 2 2,0 , · · · , δv ). At the next iteration, the second vectors will be treated as the initial vectors to be substituted into Eq. (6.1). Similarly, repeating m times for this process gives the mth vector (δY 1 m,0 , δY 2 m,0 , · · · , δY d(N +1) m,0 ). The steps of the Gram-Schmidt ortho-normalization are given as following where ||V i m,0 || is the norm of V i m,0 , δY i m,0 , δvī m,0 (i,ī = 1, 2, · · · , d(N + 1)) is a standard scalar product. Finally, LEs can be calculated by using ln ||V i ̺ (0)||. (6.5)

Numerical studies
In this section, we will show the effectiveness of the proposed algorithm by studying a soft impacting system with a delayed feedback controller. Since the system has many coexisting attractors when grazing is encountered [29], our control objective here is to drive the system from its current attractor to a desired one. Calculating the LEs of the system could allow us to monitor the stability of the delayed feedback controller and its effective parametric regime. Now, we consider the following parameters for the impacting system,   . 7.2 (a), the largest LEs are all greater than 0 for k ∈ [0, 0.04] and the system presents a chaotic motion as shown in FIG. 5(b). The phase trajectory of the chaotic motion for k = 0.02 is presented in FIG. 7.2 (c). For k ∈ (0.04, 0.055), the largest LE decreases and suddenly increases to the neighbourhood of zero at k = 0.055 indicating a period doubling of the system. Similarly, at k = 0.065, such a fluctuation is observed again. Thereafter, the largest LE decreases dramatically, and then increases gradually from k = 0.07. For k ∈ [0.07, 1.4], both LEs are below zero, and the system has period-1 response which is demonstrated in FIG. 7.2 (d) and (e).
A critical issue of computation of the nonsmooth dynamical systems is that the accumulated computational error in grazing boundary, could lead to incorrect simulation. FIG.7.3 compares the computations of the impacting system for e = 1.2609 controlled from a chaotic response to a period-1 response by using the delayed feedback control with and without the grazing estimation algorithm. The time histories of times of impacts without (black line) and with (orange line) the grazing estimation algorithm were presented in FIG.7.3 (a) which were counted from t = 9722, and the phase trajectories from chaotic (grey line) to period-1 (red line) response were shown in FIG.7.3 (b). It can be seen from the figure that the accumulated error was built up in the times of impacts, and obvious different can be identified from t = 10411. The cause of the difference can be found from FIG. 7.3 (c) and (d), at where the time histories of displacement of the impacting system were shown. As can been seen from the figures, the system with the grazing estimation algorithm was stabilised quicker than the one without the algorithm. As the modified Euler integration method was used in this work, the difference could be significant if a higher order integration method such as the Rung-kutta integration method is adopted.

Case 0 < τ d < T
For the case of a small time delay, we present the example for τ d = T /2 in FIG. 7.4 under ζ = 0.01, e = 1.26, a = 0.7, β = 28 and ω = 0.802. It can be seen from the figures that the system has chaotic motion for k ∈ [0, 0.007] and its largest LEs are all greater than zero (green line). For k ∈ (0.007, 0.015], the system experiences transient periodic motions, and the relevant largest LEs are smaller than zero which is consistent with the result in FIG. 7.4 (b) indicating the alternation between chaotic and periodic motions. At k = 0.016, the system has a very small chaotic regime and bifurcates into a non-impact period-1 response immediately lasting until k = 0.0425 at where another chaotic regime is encountered. For k ∈ [0.0425, 0.045], the system has chaotic response in most of the regime, but has small window of period-3 response in k ∈ [0.044, 0.04475]. After k = 0.045, the non-impact period-1 response presents again as the control parameter k increases. To compare FIG. 7.4 (a) and (b), the evolution of the calculated LEs is consistent with system's bifurcations, and phase trajectory from each regime has been illustrated in FIG. 7.4 (c-f).

Conclusions
This paper studies a numerical method for calculating the LEs of time-delayed piecewise-smooth systems by using a soft impacting system under the delayed feedback control with a particular focus on its near-grazing dynamics.
The main tasks are to build an effective variational equation and obtain the Jacobian matrix for the delayed impacting system. As the delayed system is infinite dimensional, the delayed impacting system can be approximated by finite dimensional systems, which can be linearised by the modified Euler integration method at each time step. Then the system was discretised by constructing the Poincaré map, and perturbation was introduced to obtain the variational equation. Finally the Jacobian matrix was obtained through combining all the approximated systems linearised from the variational equation at each time step in the one period. In order to increase the convergence rate and improve computational accuracy, a grazing estimation algorithm was studied. In addition, the convergence rate of the eigenvalues of the Jacobian matrix was studied by the spectral theory of the evolutionary operator. In particular, the delayed impacting system was described as the evolutionary operator whose convergence was the same as the relevant nonzero eigenvalues of the Jacobian matrix, therefore guaranteeing the reliability of the proposed numerical method.
Our numerical studies considered two scenarios for the delay time in the impacting system which were larger (τ d ≥ T ) and smaller (0 < τ d < T ) than the period of excitation. Both cases showed that the calculated LEs were consistent with the bifurcations of the system, and the grazing estimation algorithm was more accurate and effective for calculating the nonsmooth dynamical systems.