An efficient algorithm for solving piecewise-smooth dynamical systems

This article considers the numerical treatment of piecewise-smooth dynamical systems. Classical solutions as well as sliding modes up to codimension-2 are treated. An algorithm is presented that, in the case of non-uniqueness, selects a solution that is the formal limit solution of a regularized problem. The numerical solution of a regularized differential equation, which creates stiffness and often also high oscillations, is avoided.


Introduction
Piecewise-smooth dynamical systems arise in many applications and they are an active field of recent research. Historically, one of the first examples is Coulomb friction in mechanical systems, where the force of friction is proportional to the sign of velocity (see [5]). Many interesting applications can be found in the monograph [6]: relay control systems, where the control variable admits jump discontinuities; converter circuits, where switching devices lead to a non-smooth dynamics; models in the social and financial sciences, where continuous change can trigger discrete actions. Discontinuity points are also created by the activation/deactivation of inequality constraints in mixed constrained optimization problems. See [24] for a particular application arising in the modelling of atmospheric particles.
For a mathematical formulation of the problem we consider discontinuity hypersurfaces Σ j = {y ∈ R n | α j (y) = 0}, j = 1, . . . , d, (1.1) where α : R n → R d (with d < n) is assumed to be sufficiently differentiable and such that these hyper-surfaces intersect transversally. We denote the discontinuity set by Σ = d j =1 Σ j . The hyper-surfaces Σ j divide the phase space R d \ Σ into 2 d open regions R k = y ∈ R n k j α j (y) > 0 for j = 1, . . . , d , (1.2) where k = (k 1 , . . . , k d ) is a multi-index with k j ∈ {−1, 1}. The discontinuous dynamical system is then given bẏ We assume that the functions f k (y) are defined in a neighbourhood of the closure of R k and that they are sufficiently differentiable. In the discontinuity set Σ the right-hand side of (1.3) is considered to be multi-valued with values from the neighbouring domains. We are thus concerned with a differential inclusion and we adopt a restriction of the approach by Filippov [13,14] for the concept of solutions. Besides classical solutions, which cross the discontinuity surfaces, there are also sliding modes evolving in the discontinuity set Σ. Closely connected to a discontinuous dynamical system is a regularization, where the jump discontinuities are replaced in an ε-neighbourhood by a continuous transition. In this way the differential inclusion is transferred to an ordinary differential equation. It is natural to consider regularizations because, as mentioned in [6, p. 1], ". . . there is strictly speaking no such thing as a piecewise-smooth dynamical system and that in reality all physical systems are smooth". This is precisely what happens in the analysis of gene regulatory networks [12,26], where steep sigmoid-type nonlinearities are approximated by step functions.
Among numerically sound approaches for approximating the solution of (1.3) let us mention the following two: -Algorithm based on event detection. One locates accurately the time instants when the solution enters a new discontinuity surface (or satisfies a criterion for exiting a surface), one stops the integration and investigates the possible solutions leaving the actual point, and then one continues the integration with a new vector field. The disadvantage of this approach is that at the actual point the discontinuous problem can have more than one solution (sometimes even infinitely many), and it may be laborious to follow all of them. -Regularization. One solves numerically the regularized ordinary differential equation, which provides a unique approximation. Here, the difficulty is the choice of the regularization parameter ε > 0. To obtain a good approximation of the solution of (1.3) a very small ε is required. This implies that the regularized differential equation is stiff and sometimes highly oscillatory, so that the numerical integration may become expensive.
Early work on solving piecewise-smooth dynamical systems that is based on detecting, locating, and passing the discontinuity is published in [4,15,27]. For a survey we refer to [10]. There are some recent publications (including Matlab codes for solving piecewise-smooth dynamical systems), like those of [29] and [3], that are reliable and carefully compute the switching points between classical and Filippov solutions. All these publications are restricted to classical solutions and to sliding modes in codimension 1. Our main interest is the situation, where codimension-2 sliding modes can occur.
In the present work we propose an algorithm that combines the advantages of both approaches. With event detection we solve the discontinuous problem (without any ε) but, instead of following all solutions in the case of non-uniqueness, we propose to select the solution which can formally be interpreted as the limit solution (for ε → 0) of a regularized differential equation. This selection is partly done on the basis of the classification in [17].
In Section 2 we recall concepts needed for the understanding of the present article (relation between sliding modes and differential-algebraic equations of index 2, regularization, hidden dynamics, and scaling invariance). The structure of the algorithm for solving the discontinuous system (1.3) is given in Section 3. The main part (Section 4) presents in an algorithmic way the switching between different kinds of solutions at the discontinuity hyper-surfaces. This part is independent of the regularization, in contrast to Section 5, where a justification of the algorithm (based on the hidden dynamics) is given. The article finishes with some comments on the implementation (Section 6) and a conclusion (Section 7).

Solution concept and regularization
The definition of Filippov solutions for a discontinuous dynamical system (1.3) is ambiguous, because in the intersection of discontinuity hyper-surfaces a convex combination of the adjacent vector fields has too many degrees of freedom. We restrict our study to special convex combinations having m parameters in the intersection of m hyper-surfaces Σ j . Such convex combinations (for m = 2) are called "blending" in [1] and "bilinear interpolation" in [7,8], see also [9,25]. For arbitrary m they are called "convex canopy" in [20]. We consider regularizations that are closely connected to such convex combinations, and we call them "multi-linear interpolation".

Solution concept -classical solutions and sliding modes
For a fixed multi-index k = (k 1 , . . . , k d ) with k j ∈ {−1, 1} the equation (1.3) is a regular ordinary differential equation on the open domain R k , and the standard theory on existence, uniqueness, and continuous dependence on parameters and initial values applies. In this case the solution of (1.3) is called classical.
We next extend the concept of solution to the discontinuity set Σ. For an index vector k = (k 1 , . . . , k d ) with k j ∈ {−1, 0, 1} (note that now k j can also be zero) we consider the set and if at least one component where m counts the number of elements k j being equal to zero. For k = (k 1 , . . . , k d ) we define I k = {j | k j = 0}, and we let which collects the index vectors such that R touches R k . With this notation we consider the differential-algebraic equation (DAE) with algebraic variables λ j , j ∈ I k . In the following we denote the right-hand side of the differential equation in (2. 2) by f k (y, λ k ), where λ k is the vector that collects λ j , j ∈ I k . Differentiating the algebraic constraint of (2.2) with respect to time yields which represents m nonlinear equations in m unknowns λ j , j ∈ I k . We assume that the Implicit Function Theorem can be applied to guarantee that locally λ k can be expressed as function of y. This implies that the DAE has index 2. The special case I k = ∅ includes classical solutions of (1.3), because in this case N k = {k} consists of only one element and the empty product in (2.2) is interpreted as 1.
For λ j ∈ [−1, 1] the vector field in (2.2) is a convex combination of the vector fields f (y) (with ∈ N k ) which are defined on the open domains touching R k . The solution of (2.2) is therefore a Filippov solution.

Definition 2.1
Consider an index vector k with I k = ∅ and let m = |I k | be the cardinality of I k . Then, a solution (y, λ k ) of the differential-algebraic equation (2.2) is called a codimension-m sliding mode in the set R k as long as λ j ∈ [−1, 1] for j ∈ I k .
For a consistent initial value of (2.2), i.e. y(0) ∈ R k and λ k (0) given by (2.3), any technique for the numerical solution of DAE's of index 2 can be applied. Such techniques are explained in detail in the monographs [2,19].

Definition 2.2
A piecewise-smooth, continuous function y : [0, T ] → R n is called a solution of the discontinuous dynamical system (1.3), if there exists a finite partition 0 = t 0 < t 1 < t 2 < . . . < t N = T , such that the following is true: for every subinterval [t i , t i+1 ] there exists k i ∈ {−1, 0, 1} d with m i = |I k i | such that the restriction of y(t) to this interval is a codimension-m i sliding mode in the set R k i (a classical solution if I k i = ∅).

Regularization
We are interested in solutions of (1.3) in the sense of Definition 2.2 that can be considered as the formal limit of a regularized differential equation, where jump discontinuities in the vector field are smoothed out. For this we consider a transition function π(u), which is assumed to be continuous, piecewise-smooth, and satisfies π(u) = −1 for u ≤ 1 and π(u) = 1 for u ≥ 1. We also assume that π (u) > 0 for u ∈ (−1, 1), and that π(u) is centrally symmetric. A typical example is π(u) = u for |u| ≤ 1 (see Fig. 1).
For a discontinuous dynamical system (1.3) we consider the regularizatioṅ where u j = α j (y)/ε. We denote the right-hand side of this regularized differential equation by f y, π(u 1 ), . . . , π(u d ) . The complete phase space (including the discontinuity set Σ) is the union of 3 d sets where k = (k 1 , . . . , k d ) with k j ∈ {−1, 0, 1}. For the case that all k j = 0, we have that R k ε ⊂ R k , and = k is the only vector for which the product in (2.4) is non-zero. Therefore, on the set R k ε the regularization coincides with the differential equatioṅ y = f k (y) of the un-regularized problem.
For k with I k = ∅ the set R k ε approximates R k . On the set R k ε only the vectors ∈ N k give rise to a non-vanishing product in (2.4). Since j π(u j ) = k j π(u j ) = 1 for ∈ N k and j ∈ I k , the regularized differential (2.4) becomeṡ which is in complete analogy to (2.2). If m denotes the cardinality of I k , then for m = 1 the sum in (2.6) consists of two terms (linear interpolation), for m = 2 it consists of four terms (bilinear interpolation), and in general it consists of 2 m terms.

Hidden dynamics
A justification of our algorithm is based on the study of the solution of the regularized differential equation, when it is close to an intersection of discontinuity surfaces. In the region R k ε it follows from (2.6) that u i = α i (y)/ε satisfies which is a singularly perturbed differential equation. Close to a point y * ∈ R k of the discontinuity manifold it can be studied by separating a transient part from the smooth solution. For this we introduce the fast time τ = t/ε, we denote the derivative with respect to τ by a prime, and we substitute the constant vector y * for y. This yields which is a regular dynamical system for u i , i ∈ I k . It is called hidden dynamics (a term coined in [21]). We expect that this system credibly describes the transient behaviour of the solution of the regularized differential equation. 1

Special case of two intersecting surfaces
We assume that only two components of k are zero, say, k 1 = k 2 = 0. We then have I k = {1, 2} and N k consists of four elements. The differential (2.8) of the hidden dynamics is then given by, for i = 1, 2, where in the notation 1]. We also denote de right-hand side of (2.9) by g i π(u 1 ), π(u 2 ) . This 2-dimensional system has been discussed in detail in [17,Section 5]: how initial values are determined by the incoming solution, how the behaviour of the solution for τ → ∞ determines which kind of solution (classical or sliding) will be followed by the regularized equation, how a geometric study of the flow is possible, etc. We note that the right-hand side of (2.9) is a quadratic polynomial in π(u 1 ), π(u 2 ), which vanishes on a hyperbola in the π(u 1 ), π(u 2 ) -space. Throughout the present work we consider the transition function of Fig. 1.
The study of the hidden dynamics is an essential tool for designing the algorithm proposed in the present paper. A whole monograph [22] is devoted to this topic. Let us also mention the work [23], which concentrates on the 2-dimensional system (2.9). On the basis of singular perturbation theory it discusses stability of sliding, and it shows that there exists at most one stable sliding vector field.

Scaling invariance
A substitution α j (y) → κ j α j (y) with κ j ≥ 1 neither changes the discontinuous hyper-surfaces and the open regions R k nor the solution of the discontinuous dynamical system (1.3). However, it changes the regularization (2.4) (u j = α j (y)/ε is replaced by κ j u j ) and therefore also the solution of the regularized differential equation. Consequently, also in the hidden dynamics the expression π(u j ) is replaced by π(κ j u j ).
One of our aims is to design an algorithm for the numerical solution of (1.3) that is invariant with respect to such a scaling.

Solving piecewise-smooth dynamical systems
Typically, a numerical algorithm for solving piecewise-smooth dynamical systems (1.3) is composed of three parts: -Computation. Use any code that permits to solve the index-2 differentialalgebraic equation (2.2) starting at consistent initial values. Techniques and codes are well documented in text books like [19] and [2]. In the beginning one is usually concerned with a classical solution, for which I k = ∅, so that all y ∈ R k are consistent. At a transition point t i the initial value is determined by continuity. -Event location. The code has to be equipped with an event location algorithm that stops the integration either (i) when the solution enters a new discontinuity surface or (ii) when one of the Lagrange multipliers λ j leaves the interval [−1, 1] or (iii) when the solution λ k of the algebraic system (2.3) ceases to exist in the unit cube or becomes unstable. Since all λ j are functions of y, each of the conditions gives raise to an algebraic relation g y(t) = 0. Event detection strategies are made for finding such points. This defines a new grid point t i . Algorithms for event location are discussed in [28] (based on the BDF code DASSL) and in [16,24] (based on the implicit Runge-Kutta code RADAU5). -Switching. As soon as a new transition point t i is detected, one can check all possible multi-indices k which, for the present solution value, give raise to a meaningful solution in the sense of Definition 2.1. In the case of non-uniqueness one can follow all possible solutions (which may be laborious and inefficient) or one can select one of them -but which one? The present article is devoted to a theoretically founded switching criterion based partially on the classification of [17]. It provides a solution that can be considered as the formal limit of a regularized problem.
The present work focuses on the switching algorithm. It is neither our intention to give details on the numerical computation of differential-algebraic equations of index 2 nor to discuss techniques for event location.

Switching algorithm
The idea is to select a solution that can be considered as the limit solution of a regularized differential equation. For the case that a solution enters a codimension-2 discontinuity the algorithm is based on the classification of [17]. The treatment of the more challenging situation of exiting a codimension-2 discontinuity is new. Depending on whether, on the interval [t i−1 , t i ], the solution is a classical solution or a sliding mode, the switching algorithm at t i is discussed in the following subsections: for a classical solution in Section 4.1; for a codimension-1 sliding mode in Section 4.2; for a codimension-2 sliding mode in Section 4.3.
for an accumulation of grid points and spiraling solutions in Section 4.4.

Classical solution
We consider a classical solution of the differential equation (1.3) for t ≥ t i−1 until it enters a discontinuity surface at time t i . Without loss of generality we assume that the discontinuity surface is Σ 1 = {y | α 1 (y) = 0}. Removing irrelevant indices from the vector k (for notational convenience), we assume the classical solution to be in R −1 = {y | α 1 (y) < 0} with vector field f −1 (y). On the opposite side of Σ 1 the vector field is f 1 (y). We assume that the solution enters transversally the surface Σ 1 , so that f −1 1 := α 1 (y)f −1 (y) > 0 at the entry point. We then distinguish the cases, where f 1 1 := α 1 (y)f 1 (y) is positive or negative (see Fig. 2). We do not consider the non generic situation, where this expression vanishes.
If f 1 1 > 0, the only possible solution is classical in the region R 1 . If f 1 1 < 0, there is no classical solution leaving the solution point in Σ 1 . The solution in R 0 = Σ 1 is defined by the DAE (2.2), i.e.
Differentiating the constraint with respect to time yields which determines λ as function of y, namely, The initial value for (4.1) is defined by continuity for y, and for λ it satisfies −1 < λ < 1, because f −1 1 > 0 and f 1 1 < 0.

Codimension-1 sliding mode
Suppose that we are concerned with a codimension-1 sliding mode along Σ 1 for t ≥ t i−1 . The type of solution can change at some t i > t i−1 , if either it exits the discontinuity surface Σ 1 or it enters an additional discontinuity surface, say Σ 2 . Both situations are discussed in the next two subsections.

Entering the intersection Σ 1 ∩ Σ 2
As before we disregard irrelevant indices from the index vector k, and we keep only those corresponding to Σ 1 and Σ 2 . We consider a codimension-1 sliding along R 0,−1 = {y | α 1 (y) = 0, α 2 (y) < 0}, and we generically assume that where we use the notation f k j = α j (y)f k (y) for j ∈ {1, 2} and k = (k 1 , k 2 ) (all vector fields are evaluated at the entry point). The first two inequalities of (4. The algorithm presented in the two figures needs some more explanations. In addition to classical and codimension-1 solutions we have to consider codimension-2 solutions. They are defined bẏ subject to the algebraic constraints α 1 (y) = 0 and α 2 (y) = 0. The right-hand side of (4.5) is denoted by f 0,0 (y, λ 1 , λ 2 ). Differentiating the algebraic relations with for j ∈ {1, 2}. We also use the notation g j (y, λ 1 , λ 2 ) = 0 for this equation. For the existence of a locally unique solution (λ 1 , λ 2 ) of the system (4.6), we assume that the Implicit Function Theorem can be applied, which means that the 2 × 2 matrix is invertible. For a fixed value of y the equation (4.6) represents a hyperbola with vertical and horizontal asymptotes in the (λ 1 , λ 2 )-space. We are only interested in values (λ 1 , λ 2 ) lying in the square [−1, 1] × [−1, 1] (which we sometimes call unit square). When, in the algorithms of Figs. 5 and 6, we write that a solution (λ 1 , λ 2 ) of (4.6) exists (or not) in the unit square, we mean only solutions on the branch of the first hyperbola (j = 1) that crosses the bottom side of the square. The expression f 1,0 2 (λ 2 ), appearing in the switching algorithm, is defined by f 1,0 2 (λ 2 ) = α 2 (y)f 1,0 (y, λ 2 ), and f 1,0 (y, λ 2 ) is the vector field of (2.2) for k = (1, 0).

Remark 4.1
In the situation, where the algorithm of Figs. 5 and 6 proposes a codimension-2 sliding, the discussion of [17] shows that in the beginning of the sliding the solution (λ 1 , λ 2 ) of (4.6) is such that the determinant of G(y, λ 1 , λ 2 ) is positive and at least one of its diagonal elements is negative. This is important for the strategy in Section 4.3.3.

Codimension-2 sliding mode
Suppose that there is a codimension-2 sliding mode along Σ 1 ∩ Σ 2 for t ≥ t i−1 . It is characterized by the existence of λ 1 (t), λ 2 (t) ∈ (−1, 1) 2 satisfying the polynomial system (4.6). The type of solution can change at some t i > t i−1 , if either (a) the pair λ 1 (t), λ 2 (t) leaves the unit square (−1, 1) 2 , (b) the solution λ 1 (t), λ 2 (t) of (4.6) becomes double and ceases to exist, (c) the matrix (4.7) changes stability (see Remark 4.1), (d) the sliding mode enters an additional discontinuity surface, say Σ 3 . The algorithms for the situations (a), (b), and (c) are presented in the following subsections, their justification is discussed in Section 5. In this paper we do not consider the situation (d), because we are not aware of results on the limit solution of the regularized differential equation close to a codimension-3 discontinuity surface.

Exiting Σ 1 ∩ Σ 2 from a codimension-2 sliding -type (a)
We assume that at t = t i the pair λ 1 (t), λ 2 (t) of the system (4.6) leaves the unit square at one side (generically, we can exclude the corners). Without loss of generality we can assume that hold. The proposed switching, based on [18, Theorems 2 and 3], is shown in Fig. 7. All vector fields are evaluated at the exit point. According to the decision tree the solution of the discontinuous problem either continues, beyond the exit point, as a codimension-1 sliding mode or as a classical solution.
Note that all situations of Fig. 7 admit further solutions (classical or codimension-1). The proposed algorithm chooses the solution that can be realized as the limit of a regularized differential equation.

Exiting Σ 1 ∩ Σ 2 from a codimension-2 sliding -type (b)
We assume that the two (real) solutions of the system (4.6) coalesce in the unit square and disappear at t = t i . This implies that at this time instant the two hyperbolas in the (λ 1 , λ 2 )-space are tangential at λ 1 (t i ), λ 2 (t i ) . Without loss of generality we assume that the hyperbolas have positive slope which, expressed in terms of the vector fields, is (see Lemma 6.3 of [17]) (otherwise we reflect the picture at the vertical axis, i.e. change the sign of λ 1 ). Moreover, we assume that the hyperbola corresponding to α 2 (y) lies above that of α 1 (y) (4.10) (otherwise we exchange α 1 and α 2 ). Denoting the left-hand expression in (4.6) by g j (λ 1 , λ 2 ) and the derivative with respect to λ i by ∂ i , we distinguish cases according to the signs of ∂ i g j . Figure 8 gives a complete characterization of all possible situations (we shall explain later in Section 5.2 that the apparantly missing situation ∂ 2 g 1 > 0, ∂ 1 g 2 < 0 cannot arise at a vanishing stationary point). As in previous  There is only one situation (blue box) with a unique solution. In the case of nonuniqueness our algorithm selects the solution which can be interpreted as the formal limit of the solution of a regularized problem.

Exiting Σ 1 ∩ Σ 2 from a codimension-2 sliding -type (c)
A stationary point of the hidden dynamics (2.9) (corresponding to a solution (λ 1 , λ 2 ) of (4.6)) is asymptotically stable if both eigenvalues of (4.7) have negative real part. This is equivalent to det G(y, λ 1 , λ 2 ) > 0 and trace G(y, λ 1 , λ 2 ) < 0. (4.11) As explained in [17,Section 8] the trace of the matrix G is not scaling invariant. It is shown that, if at least one of the diagonal elements of G is negative, there exists a scaling that makes the stationary point asymptotically stable. The condition for a stationary point to be asymptotically stable after a suitable scaling of the constraints therefore becomes where G i,j stands for the elements of the matrix G of (4.7); see also Remark 4.1. Our (scaling invariant) strategy is to exit a codimension-2 sliding, if one of the two conditions in (4.12) becomes violated.
(1) Assume first that at time t = t i we have det G(y * , λ * 1 , λ * 2 ) = 0 (i.e. det G(y, λ 1 , λ 2 ) changes from positive to negative), while the second condition of (4.12) still holds. Due to the special structure of (4.6) both of its solutions coalesce at t = t i and, while det G(y, λ 1 , λ 2 ) changes from positive to negative for the actual solution, it changes from negative to positive for the other solution. In this situation, we propose to continue with a codimension-2 sliding, and we take for (λ 1 , λ 2 ) the solution of (4.6) for which the determinant of G(y, λ 1 , λ 2 ) is positive.

Accumulation of grid points, entering Σ 1 ∩ Σ 2 through spiraling
We consider the situation, where a solution of (1.3) enters the intersection at y ∈ Σ 1 ∩ Σ 2 by spiraling inwards. This can be clockwise or counterclockwise. Assuming the second, this is the case, if the vector fields, evaluated at y, satisfy and if the contractivity condition holds. Under these two assumptions the solution of the discontinuous system (1.3) converges to Σ 1 ∩Σ 2 in finite time. It spirals around Σ 1 ∩Σ 2 and produces an infinity of grid points that converge geometrically to the entry point. From there on we have a codimension-2 sliding.

Justification of the algorithm
In the situation of The philosophy of the presented algorithm is that in the situation of nonuniqueness we choose a solution that can be interpreted as the limit for ε → 0 of the solution of a regularization. The combined system (2.6)-(2.7) is a singularly perturbed differential equation which is typically studied by asymptotic expansions in powers of ε and by separating slow and fast dynamics. This, however, is not always possible in the present context. The experiment of [18,Section 4.2] even demonstrates the lack of an expansion in integer powers of ε.
On the other hand, close to a value y * in the discontinuity manifold, it is expected that the hidden dynamics (2.8) reproduces well the bahaviour of the solution of the discontinuous equation. For the case that α j (y) is an affine function of y and that the vector fields are constant in a neighbourhood of y * , the solution of (2.8) describes the functions u i = α i (y)/ε without any error.
We thus trust the hidden dynamics and we propose to select the solution after a switching point according to the behaviour the hidden dynamics. This implies that we have a transition to (c.f. [

Justification of the algorithms of Sections 4.2.2 and 4.3.1
The algorithms of Figs. 5 and 6 are just a transcription of Theorems 6.1 and 6.2 of [17], where the conditions are written in terms of the four vector fields rather than in terms of the vector field of the hidden dynamics.
The assumption (4.3) is equivalent to (6.1) of [17]. By Lemma 6.3 of [17] the condition ∂ 2 g 1 (u 1,0 , −1) < 0 (left turning situation) in [17] is equivalent to Fig. 5). Item (a) of Theorem 6.1 in [17] corresponds to the existence of a solution of (4.6) in the unit square. In the notation of the present work the expression g β (1, u * ) of [17, Theorem 6.1] is equal to g 2 (1, λ 2 ) = f 1,0 2 (λ 2 ), which appears in Fig. 5. The situation in [17,Theorem 6.1], where the solution of the hidden dynamics approaches a limit cycle around a stationary point, corresponds to high oscillations of amplitude O(ε) in the solution of the regularized differential equation, and to a codimension-2 sliding in the discontinuous system. Items (b) and (c) of Theorem 6.1 in [17] correspond to the part in Fig. 6, where the system (4.6) does not have a solution in the unit square that lies on the branch of the hyperbola crossing the bottom line of the square.

Justification of the algorithm of Section 4.3.2
The first seven pictures of Fig. 11 show the vector fields corresponding to the seven situations of Fig. 8. The hyperbolas g 1 = 0 and g 2 = 0 are drawn in blue. By (4.9) both hyperbolas have positive slope and by (4.10) the hyperbola for g 2 = 0 lies above that for g 1 = 0. We mark the hyperbolas with an arrow which make them to an oriented path. To the left of the hyperbola g 1 = 0 we have by convention g 1 > 0 so that the vector field points to the right, and on the other side the vector field points to the left. Similarly, to the left of g 2 = 0 the vector field points upwards, and downwards on the other side. We start with the assumption ∂ 2 g 1 > 0, ∂ 1 g 2 > 0. The branches of the hyperbolas passing through the vanishing stationary point are directed upwards. This can be observed in the pictures (a1) and (a2). In the picture (a1), where f −− 1 > 0, the solution starting at the vanishing stationary point leaves the unit square at the bottom side. This gives rise to a codimension-1 sliding in R 0,−1 . In this situation (again picture (a1)) there is also a classical solution leaving the intersection into the region R −1,1 , which however cannot be realized as the limit of a regularization. If f −− 1 < 0 (picture (a2)), the solution leaves the unit square at the lower left corner, which gives rise to a classical solution in R −1,−1 . If the branch of the hyperbola g 2 = 0 would leave the unit square at the right side, there would be another classical solution in R 1,1 , which however is irrelevant.
All other situations, namely (b1), (b2), (b3), (c1), (c2), can be explained similarly by looking at the vector fields in Fig. 11. All of them, with the exception of (b2), admit a second (non relevant) solution by slightly modifying the hyperbolas. For example, in the situation (b1) one can change the hyperbola g 1 = 0 such that it Fig. 11 Vector field of the hidden dynamics in the situation where both stationary points coalesce. The first seven pictures correspond to the colored boxes of Fig. 8 in the same order. Red arrows represent the solution after the switching enters at the left side and such that its second branch surrounds the lower right corner. Consequently, there is also a classical solution in R 1,−1 . Slight modifications of the hyperbolas permit to produce a classical solution in R −1,1 for (b3), a classical solution in R 1,1 for (c1), and a classical solution in R 1,−1 for (c2).
The last two pictures of Fig. 11 treat the situation ∂ 2 g 1 > 0, ∂ 1 g 2 < 0. The picture (d1) gives the impression that more than one solution, starting at the vanishing stationary point, are possible. However, when looking at the situation just before the stationary points vanish (picture (d2)), one sees that both stationary points are unstable and that there is no bounded limit cycle in the unit square. Therefore, this situation cannot occur at the end of a codimension-2 sliding.

Justification of the algorithm of Section 4.3.3
In the situation of the first three pictures of Figs. 9 and 10 there is exactly one solution exiting the codimension-2 hyper-surface. Let us nevertheless briefly discuss the hidden dynamics in these situations. The condition G 1,1 (y * , λ * 1 , λ * 2 ) = 0, which in terms of (2.9) reads ∂ 1 g 1 (u * 1 , u * 2 ) = 0, implies that the hyperbola g 1 (u 1 , u 2 ) = 0 is degenerate. It is the union of the horizontal and vertical asymptote. The assumption (4.13) implies that at the left side of the unit square and in a neighbourhood of the stationary point (u * 1 , u * 2 ) the hyperbola g 1 (u 1 , u 2 ) = 0 is oriented to the right. The positivity of det G and of G 2,2 imply that the hyperbola g 2 (u 1 , u 2 ) = 0 crosses the stationary point (u * 1 , u * 2 ) from bottom left to top right. The additional assumption f 1,1 1 > 0 in Fig. 9 implies that the vertical asymptote of g 1 (u 1 , u 2 ) = 0 is outside the unit square. The upper pictures of Fig. 12 illustrate the four situations of Fig. 9. The vector field of the hidden dynamics is shown on the unit square and on a neighbourhood of it. The oriented hyperbolas are indicated in blue, and 20 solutions with randomly chosen initial values 10 −6 -close to the stationary Fig. 12 Vector field of the hidden dynamics for the situations discussed in Fig. 9 (upper pictures) and in Fig. 10 (lower pictures). Twenty solutions corresponding to initial values that are random perturbations of the equilibrium (u * 1 , u * 2 ) are plotted in red point are plotted in red. An inspection of the vector field shows that the solutions spiral outwards in the first picture. Apparently, there is an infinity of solutions starting at t = t i at the point y * . In the second and third situations the solutions all tend to a classical solution in R 1,1 and R −1,−1 , respectively. The fourth picture indicates the existence of two classical solutions, one in R 1,1 and the other in R −1,−1 . Figure 13 shows the basin of attraction of the two solutions. For any initial value in the grey region the solution converges to the classical solution in R −1,−1 , and for initial values in the white region it converges to R 1,1 . This clearly shows that there is non-uniqueness of the solution of the discontinuous problem.
The assumption f 1,1 1 < 0 in Fig. 10 implies that the vertical asymptote of g 1 (u 1 , u 2 ) = 0 is inside the unit square. The lower pictures of Fig. 12 illustrate the four situations of Fig. 10. The second, third, and fourth pictures are similar as before with the exception that the classical solution in R 1,1 is now a codimension-1 sliding along R 0,1 . This is because the solution of the hidden dynamics cannot cross the vertical asymptote.
An explanation for the first picture of the lower row in Fig. 12 is more tricky. The solution spirals outwards around the stationary point, but remains to the left of the asymptote. On the other hand, the only stable solution leaving the square is classical in the region R 1,−1 , which is to the right of the asymptote. How can we reach this solution? The reason is that G 1,1 (y, λ 1 , λ 2 ) is negative before t = t i , but positive after it. Hence, immediately after t = t i the hyperbola g 1 (u 1 , u 2 ) = 0 is not degenerate and the vertical asymptote is no longer a separation of solutions of the hidden dynamics. After a few spirals around the stationary point, the solution can escape to the right and follow the classical solution in R 1,−1 .

Some details for an implementation
Every code for solving ordinary differential equations having an option for event location is suitable for using the algorithm of the present work. -Exiting the intersection Σ 1 ∩ Σ 2 . In the right pictures of Fig. 12 we also have non-uniqueness. This is independent of the scaling, because in any case the basin of attraction of both solutions are non empty (Fig. 13). -Exiting Σ 1 ∩Σ 2 through spiraling. In the upper left picture of Fig. 12 the solution exits a codimension-2 sliding through spiraling. As a consequence of our scalinginvariant criterion (4.12) we have a one-parameter family of exiting solutions. With the usual criterion (4.11), we would have a two-parameter family of exiting solution, because the exit point is not unique.