COMPUTATIONAL TECHNIQUES TO LOCATE CROSSING/SLIDING REGIONS AND THEIR SETS OF ATTRACTION IN NON-SMOOTH DYNAMICAL SYSTEMS

. Several models in the applied sciences are characterized by instan-taneous changes in the solutions or discontinuities in the vector ﬁeld. Knowl- edge of the geometry of interaction of the ﬂow with the discontinuities can give new insights on the behaviour of such systems. Here, we focus on the class of the piecewise smooth systems of Filippov type. We describe some numerical techniques to locate crossing and sliding regions on the discontinuity surface, to compute the sets of attraction of these regions together with the mathematical form of the separatrices of such sets. Some numerical tests will illustrate our approach.


1.
Introduction. Many nonsmooth dynamical systems can be described by piecewise smooth (PWS) systems of ordinary differential equations (ODEs), whose theory has been extensively described in [28]. Examples include mechanical, robotic, electrical, biological, climate, and even financial models (see for instance [12,18,23,38,39]). Recently, standard techniques to study discontinuous systems of ODEs have been employed in the simulation of reaction diffusion systems with discontinuous forcing terms (see for instance [4,2,3]). The solutions of such systems can be simulated using suitable numerical methods, see for instance the book [1] or the papers [10,11,37].
While simulations can provide valuable insight into a system's dynamics, they have limited value when parameters or initial conditions are highly uncertain. In these cases, more robust information is gathered by looking at the structure of the system's flow as a whole (see [14,15,24,20,29]). By identifying specific features that steer whole families of trajectories towards one or another behaviour, one can obtain information that is robust with respect to small perturbation of the initial conditions as well as small modelling errors.
Such an analysis could be particularly valuable in a discontinuous system. Discontinuities in dynamical systems often represent physically relevant transitions, for example, switches of a mechanical or electronic relay, or crossings of a triggering threshold in chemical or biological systems. Thus, by identifying groups of trajectories that encounter different sequences of discontinuities, we can classify regions of state space that exhibit measurably different behaviour, often with important implications for the system's function.
As discussed in this paper, trajectories can be grouped by their switching sequences by looking at the set of tangencies of the flow with the discontinuity surfaces, and the family of orbits crossing these sets. In the following, after laying out some mildly restrictive hypotheses, we discuss some numerical techniques which can be used to compute such sets, and to represent them with implicit functions so as to efficiently classify initial conditions in terms of the future behaviour of their trajectory. In particular, in Section 2, crossing and attractive/repulsive sliding regions on the discontinuity surface are defined together with the notion of singular points and singular sets of a PWS system; Section 3 briefly outlines how one can compute the singular sets via standard computational continuation techniques. In Section 4, definition of separatrices and sets of attraction of a sliding region is provided; while Section 5 presents some numerical tools for the computation of separatrices and their representation in implicit function form and exemplifies their performance through several numerical simulations.
The main limitation of the techniques we propose is that, being based on the computation of separatrices, they are only practical in scenarios where such sets are at most two-dimensional, and therefore in systems of up to three dimensions.
2. Crossing/sliding regions and singular sets. For the sake of simplicity, we assume that the state space is split in only two subregions, R 1 and R 2 , separated by a switching surface (or discontinuity surface) Σ, implicitly defined by a smooth scalar function h : R n → R, so that and We assume that h ∈ C k , k ≥ 1, and that ∇h(x) = 0 for all x ∈ Σ where ∇h(x) is the gradient vector of h with respect to the x variables. Under these assumptions, we consider an n-dimensional systems of the following form:ẋ where f 1 and f 2 are smooth functions in their domain of definition. The Lie derivative of a smooth function h : R n → R along a smooth vector field f : R n → R n is defined as w(x) := ∇ h(x)f (x) and represents the rate of change of h(x) along the flow of f (x) or, equivalently, the projection of f (x) onto the normal direction to Σ at x. Similarly, the second Lie derivative w 2 (x) := ∇ w(x)f (x) represents the rate of change of w(x) along the flow of f (x). The dynamics of (3) inside the region R 1 (respectively R 2 ) is governed by the smooth vector field f 1 (respectively f 2 ). To continue solutions across the switching surface, we define solutions near Σ as any solution of the differential inclusion [28] x ∈ F (x), and The subset of points on Σ where a trajectory, coming from the subregion R 1 (respectively R 2 ), crosses Σ and enters R 2 (respectively R 1 ) is called crossing region; the subset of points on Σ where the vector fields f 1 and f 2 are both pointing towards Σ is called attractive sliding region; while the subset of points on Σ where the vector fields f 1 and f 2 are both pointing away from Σ is called repulsive sliding region.
According to the definitions above, a solution reaching Σ in the crossing region crosses it moving on the other side of Σ, while one reaching Σ in the attractive sliding region is bound to slide on Σ. In both cases, the forward-time solution is uniquely defined, which is not the case for solutions on the repulsive sliding region. The above considerations can be precisely formulated in terms of the Lie derivatives w i (x) = ∇ h(x)f i (x), with i = 1, 2, as follows.
Definition 2.1. The negative and the positive crossing regions R C 1 and R C 2 , where trajectories cross from R 1 to R 2 or from R 2 to R 1 are respectively: The crossing region, Definition 2.2. The attractive and the repulsive sliding regions on the discontinuity surface, are respectively: The sliding regionof the PWS system, Clearly, R C and R S are open disjoint (and possibly empty) sets. In case of attractive sliding motion, the forward-time solution is the integral of the unique element of f F (x) that is tangent to Σ (that is, ∇ h(x)f F (x) = 0), which yields where Note that equation (8) is not defined when w 1 (x) and w 2 (x) are simultaneously zero, that is, at points where both flows f 1 and f 2 are tangent to Σ. Dynamics at these points is known to be extremely complex [16,17]. We now turn our attention to points on Σ where the Lie derivatives vanish: Sets E 1 and E 2 are singular sets while points belonging to E i aresingular points [7]). Singular sets separate crossing region from sliding one.
During an attractive sliding motion, when a solution x(t) reaches a pointx = x(t) on E 1 \ E I or E 2 \ E I , it leaves the sliding region if w 2 1 (x) = 0 or w 2 2 (x) = 0, respectively. Whereas w 2 1 (x) = 0 or w 2 2 (x) = 0, the solution leaves the sliding region when the first non-zero Lie derivative is of even order; it remains in the sliding region otherwise. The points at which the solution leaves the sliding region are said exit points. The subset of singular points of E 1 and E 2 is called the singular set. When a solution, during an attractive sliding motion, reaches a point in E I , its behaviour and uniqueness depend on the higher order Lie derivatives of the fields f 1 and f 2 .
Finally, we define the singular entry set as the set of singular points for the system integrated backward in time. This may be considered a symmetric set with respect to the singular set.
2.1. Example. Let us consider the following two-dimensional PWS system: , and the discontinuity surface Σ defined by the zero set of h(x) = x 2 − 0.2. The gradient of Σ is given by the constant vector ∇h = [0 1] , thus on Σ we have ∇h f 1 (x) = −x 1 +1, and ∇h f 2 (x) = −x 1 −1 , which means that the singular sets are There is an attractive sliding segment, on the line x 2 = 0.2, which is given by R S A = {x ∈ R 2 : x 1 ∈ (−1, 1), x 2 = 0.2}, while the Filippov sliding vector field is given by 2}, are two crossing regions. The set E I = E 1 ∩ E 2 is empty. Remark 1. An example of taxonomy of the singular points for planar discontinuity surfaces has been proposed in [5,6]: for three dimensional discontinuous systems of Filippov type several different kinds of singular points can be defined on the boundary between the sliding and crossing region and classified using an integrationfree method based on the evaluation of the vector fields.
3. Continuation of boundaries of crossing/sliding regions. From now on, we focus our analysis on systems in 2 and 3 dimensions. It should be pointed out that numerical methods described in the following sections can be extended -in principle-to higher dimensions, but at considerable cost in terms of numerical and algorithmic complexity.
Singular sets E 1 and E 2 constitute the boundaries of the crossing/sliding regions on Σ and are described by means of nonlinear equations. When we deal with PWS systems in two and three dimensions, these sets reduce to isolated points and curves on Σ, respectively. Hence, standard root-finding algorithms can be used to compute singular sets (points) in the two-dimensional case, while standard continuation algorithms can be adopted to numerically reconstruct such sets in three-dimensional systems. In this latter case, it can be useful to write E 1 and E 2 as the zero-sets of the following system: where G i : R n → R 2 is for any x ∈ R n . It should be observed that the Jacobian matrix of G i has full rank almost everywhere on E I . Proposition 1. (See [31]). Let G : R n → R n−k (n > k ≥ 1) be a smooth function and suppose that there exists x ∈ R n where G(x) = 0 and its Jacobian matrix J G has full rank n − k. If p ∈ R n is in the kernel of J G (x), then there exists a curve x = c(t) in R n such that c(0) = x,ċ(0) = p, which lies locally in the zero set of G, that is {x ∈ R n : G(x) = 0}.
Based on the above result, available software can be used to compute the boundaries of crossing and sliding regions in systems of up to 3 dimensions. In our examples we used CL MATCONT 1 ( [21,22]) to continue the singular sets starting from points on the switching manifold by means of equations (12). Additional MATLAB code has been developed to graphically represent the different regions on the switching manifold once the singular sets are obtained. This code provides a graphical representations of the singular manifolds in the appropriate parameter space and is based on the function ezimplot3.m 2 .
3.1. Example. (See [5]). As a first example, let us consider the following PWS systemẋ where where the discontinuity surface is the zero set of h(x) = x 3 , and λ is a real parameter. Figure  , of the two singular sets of (13) for a negative value of the parameter, that is λ = −10. In this case, an attractive sliding region (which is a strip) appears on the discontinuity surface. All points on the discontinuity surface located outside the boundaries of this sliding region (depicted as red and blue lines) belong to the positive or negative crossing region. The width of the sliding strip increases when the absolute value of the parameter λ increases. Moreover, a change in the sign of λ affects the attractiveness of the sliding region, particularly it is repulsive for negative values of λ and attractive otherwise.
The top plot in Figure 2 illustrates the effect of continuous variation of the parameter λ on the singular sets; the graph shows the shape of the singular sets in the (x 1 , x 2 , λ) space. The red and blue singular manifolds correspond to the red and blue lines illustrated in Figure 1; the singular sets for λ = −5 are also reported. The bottom plot in Figure 2, instead, depicts the projection of the singular sets in the λx 1 -plane; this plot provides information about the qualitative behavior of the system for different values of the parameter λ.  [37]). As a second example, consider the undamped dryfriction oscillator with one degree of freedom

Example (See
where ω is the forcing frequency and F is the size of the Coloumb friction force. Equation (14) is equivalent to the three dimensional PWS system (3) with vector fields and discontinuity surface being the zero set of h(x) = x 2 .   points inside the region bounded by these two sinusoids constitute the sliding region on the discontinuity surface. The sliding region is attractive for positive values of the parameter F , while it is repulsive when F is negative. Points on the discontinuity surface located outside the sliding region belong to the positive or negative crossing regions. Figure 4 portrays the two singular sets as a function of parameter F , in the space (x 1 , x 3 , F ). The red and blue singular sets correspond to the red and blue curves  in Figure 3. Figure 4 highlights the configuration the system would assume in the state space when the parameter value is F = 0 (up) and F = 1 (down). It should be observed that the width of the sliding region increases as the absolute value of the parameter F increases. To answer this question, definitions of subsets containing all initial conditions which give rise to solutions which encounter the discontinuity surface in a given crossing or sliding region are provided. A further classification of such sets will help to determine -for any arbitrary initial condition-its possibly infinite sequence of sliding and crossing behaviour.  Note that, by the above definition, a separatrix is a family of orbit segments entirely contained in the closure of one of the regions R i , for i = 1, 2.

5.
Computation and representation of a separatrix. Interpolation or continuation methods can be useful to numerically compute separatrices of sets of attraction for the crossing/sliding regions. In this section, we present in some detail the techniques adopted to construct separatrix curves and surfaces for a given PWS system in 3D, and we discuss their advantages and drawbacks.

Computation of a separatrix by backward integration
is the exact solution of the PWS system in R 1 , at t > 0, starting from x 0 . An orbit, in a separatrix of R 1 , can be computed by integrating the exact flow backward in time starting from an initial condition . This property ensures that the backward integration of an orbit from a point x 0 ∈ E 1 results in an orbit that -integrated forward in the same flow-reaches the same point x 0 .
In general, numerical integration methods are not symmetric; there are however examples of symmetric implicit second order schemes, such as the implicit midpoint rule and trapezoidal rule ( [32]). Hence, to compute an orbit on a separatrix of R 1 starting from any given arbitrary point in E 1 , the symmetric implicit midpoint rule with step size −τ or the symmetric trapezoidal rule can be adopted. During backward integration, the time at which a trajectory reaches the discontinuity surface Σ (as shown in Example (6.1)) needs to be exactly located since, as stated in Definition (4.1), we are only interested in the first intersection of the trajectory with a crossing or sliding region. Hence, schemes (16) or (17) must be implemented as event driven methods (see [1,37]). In particular, if h(x k )h(x k+1 ) < 0, the event pointx (such that h(x) = 0) can be found using a zero finding routine or by means of a direct event location procedure which, after a reparametrization of the time variable, computesx in a fixed number of time steps ( [25,34]). Finally, to compute a set of orbits on a separatrix of R 1 , a symmetric scheme can be applied to a sequence of initial conditions in E 2 . Similarly, sets of orbits on a separatrices on R 2 can be approximated using symmetric backward integration starting from a sequence of initial conditions in E 1 .

Remark 2.
It should be observed that the approach we have previously described requires to integrate a given trajectory backwards for a given T . In general, a separatrix needs not be defined over the full (backward) time interval [0, −T ], since trajectories may intersect further the discontinuity surface before time −T . Given our definition of set of attraction and separatrix, it is sufficient to take the intersection of the computed set with the region R i , for i = 1, 2 to obtain the separatrix. The complete integration backward in time of the trajectories composing a separatrix, taking into consideration the intersection with further points of discontinuity is, in general, a complex numerical and theoretical problem only partially addressed in [15].

5.2.
Computation of a separatrix through 2D-continuation. The full 2dimensional separatrix can be computed by mesh reconstruction. Using the notation indroduced in Section 5.1, let ϕ t (x 0 ) be the solution of the PWS system at time t having initial condition x 0 for t = 0. Without loss of generality, we restrict ourselves to a portion of the separatrix lying in R 1 . The main idea is to use a continuation technique to numerically compute the surface implicitly defined as the zero set of the map where τ = τ (x) > 0 is the time at which the solution ϕ τ (x) first hits the discontinuity surface Σ. Recall that w 2 has been defined in Section 2. We point out that F 1 is well defined at points x ∈ R 1 for which such a value of τ (x) exists (i.e. initial conditions for which the solution does eventually reach Σ), and that it is smooth as long as f 1 (ϕ τ (x)) is transverse to Σ. Therefore the object of interest, hereafter denoted by S 1 = {x ∈ R 1 : F 1 (x) = 0}, is a two dimensional manifold (a surface) when the PWS system is three dimensional. The continuation can be performed via a predictor-corrector algorithm such as the one presented in [30], where it was designed and employed in the context of rigorous computation of solutions' manifolds for infinite dimensional dynamical systems. Here we only sketch the main features of the algorithm, and refer to [30] for more details.
The algorithm is essentially an advancing front method. Assuming that an initial point x (0) on the separatrix S 1 is given, it proceeds to explore the surface from there. First, six predictions x (1) pred , . . . , x (6) pred are constructed on the tangent space T x (0) to S 1 at x (0) , and then projected down on the surface via a stationary Gauss-Newton method: where J F1 (x (0) ) is the Jacobian of F 1 at x (0) and J F1 (x (0) ) + is its Moore-Penrose pseudoinverse. The six predictions lie at the vertices of a regular hexagon, whose size is chosen by the user. See Figure 5.
The actual implementation of each projection step follows closely the one adopted in [30], but is simplified by the fact that F 1 is real-valued, so that its Jacobian matrix is actually a gradient vector. Let be a QR factorization of J F1 (x (0) ) T , with q 1 ∈ R 3×1 and the columns of Q 2 ∈ R 3×2 being mutually orthonormal, and r being the Euclidean norm of J F1 (x (0) ). Then projection of the six predictions proceeds as follows: for k=1, 2, . . . , 6 do Note that the columns of Q 2 form an orthonormal basis for ker(J F1 (x (0) )), which spans T x0 . Therefore, predictions are formed through taking linear combinations of the columns of Q 2 and the iterates above proceed in the direction perpendicular to T x (0) . At this stage, six new points x (1) , . . . , x (6) on S 1 and six new simplices (triangles) have been computed, which form what is called a patch centered ad x (0) . The algorithm then proceeds by completing the patches centered ad x (1) , x (2) , . . . in the same predictor-corrector fashion. More precisely, given a point x (k) on the front of the mesh, it identifies the external "gap angle" that needs to be filled to complete the patch, builds new predictions x (k+1) pred , x (k+2) pred , . . . on the tangent plane T x (k) towards the unexplored portion of the manifold, and projects them down on S 1 similarly to how it is described above. See Figure 6. The number of predictions is chosen aiming at dividing the gap into angles which are as close as possible to π/3, while the length of the predictors is chosen adaptively in order to keep a balance between quick convergence of the projection steps, quality of the numerical representation, and overall cost of the algorithm.
The algorithm stops if an exception occurs (e.g. the stepsize reaches a user defined minimum threshold, projection on the manifold fails) or the portion of manifold sought by the user has been successfully approximated. Upon completion, the algorithm produces a database of nodes (points) and simplices (triangles) that provides a faithful numerical representation of the separatrix S 1 .
Remark 3. We point out that evaluation of F 1 relies on integration of an ODE, and is therefore performed numerically. The Jacobian for F 1 (required by both prediction and correction steps) is approximated via centered finite differences, for gap angle   Remark 4. Herebelow we summarize the main advantages of the 2D continuation method: 1.
Stepsize is chosen adaptively based on the corrector's length, so that the selection of points of S 1 is influenced by the local curvature of the manifold rather than by the dynamics of the PWS system; 2. Continuation of each portion of the separatrix starts off at a single point x (0) ; no knowledge of the boundary between crossing and sliding region on the discontinuity surface Σ is required; 3. With this method, the presence on the separatrix of invariant sets (e.g. equilibrium points) which attract orbits in backward time does not pose any concern since integration is only performed in forward time (towards Σ) and no two nodes on the separatrix need to be connected by an orbit.
Remark 5. Notice that for a 3D system, as an alternative to the method in this section, the separatrix can be computed using standard continuation software (e.g. Matconto or AUTO), as the continuation of a backward-time trajectory rooted on the 1-dimensional singular set, with the final condition as continuation parameter. A similar approach was illustrated in [14].

5.3.
Representation of the separatrix as an implicit function. Once a collection of points on a separatrix is provided by any of the means outlined before, the surface which implicitly describes the separatrix can be defined in implicit function form using suitable interpolation techniques. Deriving a mathematical expression of this surface is useful to determine the relative position of points with respect to the separatrices and to establish the behaviour (reaching the crossing or sliding region) of trajectories starting from a selected point. Similar techniques have been used in for approximating the basins of attraction for equilibrium points of dynamical systems (see for example [13,26]).
In this section, we illustrate the main phases of a technique based on radial basis function interpolation which can be used to reconstruct the separatrix surface in a given domain of the phase space.
Given a set of sampled data points, X K = {x i ∈ R n : i = 1, . . . , K} on or near an unknown surface M of R n , the goal is to determine a surface M * which is a reasonable approximation to M and is the zero iso-surface of a smooth function m : R n → R (i.e., for all points x ∈ R n , it satisfies the implicit equation m(x) = 0). Non-standard interpolation techniques based on radial basis functions (RBF), ψ : [0, ∞] → R, centered at x i , can approximate M * , so that being · the Euclidean norm and c i real coefficients. In n = 2 or n = 3 dimensions, the RBFs functions most frequently used are ψ(r) = r 2 log r and ψ(r) = r 3 , r ∈ R ( [9,40,36]), while compactly supported RBFs (such as φ(r) = (1 − r) + ) or shape parameter depending RBFs (as ψ(r) = exp[(− r) 2 ], being a given parameter) have been recently proposed (see [42,41]). Since, imposing only the interpolation conditions m(x i ) = 0, for i = 1, . . . , K, leads to the trivial solution (i.e., the identically zero function), an auxiliary set of offsurface points needed to be defined to add supplementary significant interpolation conditions [41,27]. Assume that a set of oriented normals to M at each x i is also provided, namely, n i ∈ R n , i = 1, . . . , K. These normals allow to construct off and in surface points for each x i ∈ X K where δ > 0 is a small step size whose specific magnitude should be adequately chosen for a good surface fitting ( [40,27]). Note that if normals are not explicitly given at each point, they can be approximated as described in Section 5.3.1.
The union of X +δ = {x K+1 , . . . , x 2K }, X −δ = {x 2K+1 , . . . , x 3K }, and X K constitutes the overall set of points on which the interpolation conditions are assigned. The interpolant is computed by determining a function P(x) whose zero contour interpolates the points x i and whose in and out offset contours interpolate the offsurface points in X + δ and X − δ . Particularly, the complete interpolation conditions are where the values of +1 and −1 are arbitrary. Conditions (22) lead to solve the 3K × 3K linear system where P = [p ki ], with p ki = ψ( x k − x i ), for k, i = 1, . . . , 3K, c = [c 1 , . . . , c 3K ] is the unknown vector and y is the vector whose entries are 0, 1, −1. The matrix P is symmetric and full rank if RBFs with infinity support are considered, while it is a sparse matrix for RBFs with compact support.
Remark 6. System (23) has large dimension 3K and can suffer from ill-conditioning of the interpolation matrices P which is directly connected to the order of smoothness of the basis function and to the node distribution. Hence its numerical solution has to be carefully computed using for instance partition of unit interpolation method proposed in [35,13].
5.3.1. Normals estimation. When normals n i at x i are not known, they can be estimated adopting the algorithm proposed in [33]. This algorithm firstly approximates a tangent plane at each sample point x i and then computes the normal to each obtained plane using Principal Component Analysis (PCA). In particular, the tangent plane is obtained by applying the least square method to thek-neighborhood of each x i . Thek nearest neighbors of each x i are thek closest points (in the Euclidean distance) surrounding x i (wherek is a user-specified parameter). Tangent planes are represented by their center point (centroids) and by their normal vector. The estimated normal to the tangent plane can be obtained performing an analysis of the eigenvectors and eigenvalues (or Principal Component Analysis) of a covariance matrix created from the nearest neighbors of the query point. More specifically, for each x i , the covariance matrix C is computed: wherek is the number of point neighbors considered and x is the centroid of the nearest neighbors. Being λ j the j-th eigenvalue of the covariance matrix, ordered in a decreasing way, and v j the j-th eigenvector (such that Cv j = λ j v j ), then the normal n i to x i can be approximated by thek normalized eigenvector vk, corresponding to the smallest eigenvalue λk.
In general, the normals computed via PCA could be not consistently oriented over the entire set of points. In fact, if the tangent planes of two sufficiently closed points x i and x j are consistently oriented, then n i n j = +1, otherwise it is n i n j = −1 and either n i or n j should be flipped. A consistent global orientation is difficult to obtain because of the difficulty of specifying "sufficiently closeness" condition. However, the problem of selecting a consistent orientations for the tangent planes can be modeled as a graph optimization task. Particularly, a propagation order of normal signs can be achieved by traversing the minimal spanning tree of the weighted unoriented graph whose nodes are the tangent planes of each point x i , edge (i, j) between nodes i and j exists if the tangent plane centers are sufficiently closed and the cost of an edge (i, j) is computed as 1 − |n i n j | (more details can be found in [33]).

Remark 7.
The normal approximation procedure has been implemented in MAT-LAB2014 R routines. Particularly, given a collection of points on a separatrix surface, for each point thek = 10 nearest neighbors have been obtained using knnsearch routing with Euclidean distance and kdtree research method. Then, the minimal spanning tree of the weighted unoriented graph has been obtained by the Kruskal algorithm. As initial orientation for the normals, we consider the unit normal of the tree root and then we consistently propagate it. The breath-first search algorithm has been used to traverse the tree, assigning each normal an orientation that is consistent with that of its parent. Particularly, if during traversal, the current node has been assigned a normal n i and j is the next node to be visited, then n j is replaced with −n j , if the orientation is not consistent, that is n i n j < 0.

Remark 8.
Here, we summarize the main features of reconstruction method based on RBF approximation.
1. The reconstruction of the unknown manifold is influenced by its local curvature. Particularly, the algorithm can fail to correctly reconstruct the surface when it presents sharp edges. In this case, in fact, the method to propagate normal orientation can generate ambiguous situations sharp edges. 2. The neighborhood of each data point x i depends on itsk nearest neighbors, wherek has been currently assumed to be a fixed input parameter. It should be observed that for general curve reconstruction from unorganized points the choice of this parameter could influence the stability of the algorithm. In fact, when data noise tends to dominate, the eigenvectors obtained using PCA do not reveal the correct surface tangent plane for small values ofk. However, in the separatrix reconstruction case, data points are obtained using backward integration so they contain little or no noise and small values ofk have been empirically observed to produce good reconstruction results. 3. The reconstruction of the unknown manifold is based on the availability of a set of sampled data points near or on the manifold and of the distribution of these points. These points can be obtained using symmetric backward integrators to a sequence of initial conditions in E 1 or E 1 and then randomly sampling the obtained numerical orbits. However, the presence on the separatrix of invariant sets (e.g. equilibrium points) attracting orbits in backward time can generate points with a peculiar distribution in the variable space restricting the reconstruction of the manifold only to this domain.
6. Numerical illustrations of separatrix reconstructions. The following two examples illustrate the reconstruction of the separatrix curve/surface obtained using the approach described in Section 5.3.
Different behaviours can be obtained when starting from a point in R 11 ∪R 12 . All trajectories (forward in time) will reach the segment [B, C] on Σ: since w 1 (x) > 0, a trajectory starting from R 11 will meet the sliding segment; while a trajectory starting from R 12 , firstly crosses W 1 and then reaches Σ. Since the system has no limit cycle in R 1 and the equilibrium point of f 1 , i.e. x * = (1/1.2, 0), is a repulsive point, then any trajectory starting in a neighborhood of x * approaches Σ. In Figure  7 a trajectory starting in a neighborhood of x * , crossing W 1 and then reaching the sliding attractive segment is also reported. In the region R 2 , we have w 2 (x) = −x 1 − 1 0.8+x2 ; thus W 2 = {(x 1 , x 2 ) ∈ R 2 : Starting with e 1 = (1, 0.2) (the singular set is E 1 = {e 1 }), the numerical solution of x = f 2 (x), in reverse time, will meet the curve W 2 at a finite time T 2 ∼ = 1.33 (i.e. w 2 (φ −T2 (e 1 )) = 0) and Σ at T * 2 ∼ = 1.75 (i.e. h(φ −T * 2 (e 1 )) = 0). The trajectory Γ 2 = {φ −t (e 1 ) : t ∈ [0, T * 2 ]} ⊂ R 2 can be used to devide R 2 in different regions. In particular, the trajectory (in forward time) starting from a point in R 2 but outside the region R 21 ∪ R 22 could reach the sliding segment only after crossing Σ at a point on the segment [C, D]. A trajectory starts from a point inside the region R 21 ∪ R 22 will reach the segment [A, B] of Σ. Moreover, no limit cycle is present in R 2 and the virtual equilibrium point of f 2 , x * = (−1/0.8, 0) (which is outside R 2 ) cannot be reached starting from any point in R 2 .
6.2. 3D example. Consider the three dimensional PWS system (15), with F = 1 and ω = 1, whose singular curves are sinusoids on the discontinuity plane x 2 = 0, as depicted in Figure 3. Randomly sampling points on the singular curves (blue and red dots in Figure 8) and solving backward in time the PWS system, we obtained a collection of points on the separatrix surfaces which were used to reconstruct the separatrices via RBF interpolation method. The separatrix surfaces reconstructed in the domain [−2, 2] × [−1, 1] × [−5, 25] are depicted in Figure 8 and can be useful to discuss the behaviour of trajectories starting in the chosen domain.
Particularly, the separatrix surface m 2 in R 2 (the blue surface in Figure 8) identifies two subregions R 21 = {x ∈ R 2 : m 2 (x) < 0} and R 22 = {x ∈ R 2 : m 2 (x) > 0} showing a different behaviour. Any trajectories starting from initial conditions in R 21 will move forward the crossing region on the discontinuity surface x 2 = 0, Figure 8. Reconstruction of the separatrices surfaces from some collection of points (blue and red dots on the surfaces). The black line illustrates a piece of a trajectory starting from the initial point [1,1,5] . The chosen initial condition belongs to the subregion of points in R 2 starting from which any trajectory (forward in time) reaches the crossing region on the discontinuity surface x 2 = 0.
while any trajectories starting from initial conditions in R 22 will reach the attractive sliding region. Analogously, the separatrix surface m 1 in R 1 (the red surface in Figure 8) identifies two subregions R 11 = {x ∈ R 1 : m 1 (x) < 0} and R 12 = {x ∈ R 1 : m 1 (x) > 0}. Any trajectory starting from point in R 11 will reach the attractive sliding region on the discontinuity surface, while points in R 12 give rise to trajectories which will reach the crossing region of the discontinuity surface.
The black line also depicted in Figure 8 corresponds to a part of a sample trajectory obtained integrating forward in time by an event location method (see [8,19]) from the initial point x 0 = [1, 1, 5] ∈ R 2 . As it can be observed, this initial point generates trajectories reaching the crossing region on x 2 = 0. Details on the behaviour of this sample trajectory are reported also in Figures 9 where the projections of on the x 1 x 3 , x 1 x 2 and x 3 x 2 planes, respectively; the singular curves on the plane x 2 = 0 are also depicted as red and blue sinusoid. Table 1 reports the values w.r.t. the interpolated surface obtained for some points randomly chosen into the interpolation domain. This piece of quantitative information allows to a priori know the behaviour of the trajectories starting from these points. Moreover, Figure 10 illustrates the separatrix surfaces m 1 in R 1 and m 2 in R 2 reconstructed (as described in Section 5.3) from some collections of points (red and blue dots on the surfaces) randomly chosen on trajectories obtained integrating the PWS system backward in time for a fixed T f in or until the discontinuity surface is reached for the first time. Figure 11 reports, instead, only the separatrix m 2 in R 2 together with the details of the two exit curves on Σ. These latter pictures can  give an idea of the qualitative behaviour of the system inside the selected domain: particularly, each trajectory starting from initial point external to the colored surfaces and pointing toward the discontinuity surface will hit Σ on a crossing point. On the other hand, initial condition located into the "inside" of the surface w.r.t. Σ will produce trajectories sliding on it. Figure 10. Reconstruction of some portions of the separatrix surfaces m 1 in R 1 (red surface) and m 2 in R 2 (blue surface) obtained interpolating a collection of points (red and blue dots on the surface, respectively) randomly chosen from trajectories obtained integrating -backward in time-the PWS system Figure 11. Separatrix surface m 2 in R 2 and exit curves on the discontinuity surface Σ.
Hereunder, we show how the algorithm based on the 2D continuation discussed in Section 5.2 performs in computing portions of the separatrix for the same problem.
First we have computed two portions of the separatrix on both sides of Σ. We have started the computation from two nodes obtained via successive bisections along two lines transversal to the separatrix, and have obtained a triangulation made of 1899 nodes and 3455 triangles. The computation has been halted once all the triangles inside the bounding box B = [−2, 0] × [−1, 0] × [− π 2 , 9 2 π] have been successfully computed. The triangulation is shown in Figure 12. Then a second portion of the separatrix inside the region R 1 (x 2 < 0) has been computed relaxing the constraints on x 1 , x 2 and strengthening the constraint on x 3 . In this case, a triangulation made of 8220 nodes and 16083 triangles was obtained. The triangulation is shown in Figure 13. It should be observed that, allowing for larger values of x 2 , the separatrix curves back towards Σ. Conclusions. In this paper we proposed some computational approaches to approximate the sets of attraction for the crossing/sliding regions of a nonsmooth dynamical system of Filippov type. While the singular sets (that is the boundaries of crossing/sliding regions) on the discontinuity surface Σ can be computed by Figure 13. 3D example via 2D continuation, second portion standard numerical continuation, the separatrices of the sets of attraction of these crossing/sliding regions required to carefully combina different numerical tools. For such sets we have proposed two different approaches: interpolation of points on the separatrix computed by backward integration starting from initial points on the singular sets, or genuine 2D-continuation of the separatrix starting from one given point on it.
The set of tools presented in this paper is limited in practice to the analysis of systems in 2 and 3 dimensions. However, when applicable, it can provide extremely valuable information. In the future, we hope to develop an effective user oriented Matlab tool for 2D/3D problems based on the described computational approaches.