Model Reduction of Non-densely Defined Piecewise-Smooth Systems in Banach Spaces

In this paper, a model reduction technique is introduced for piecewise-smooth (PWS) vector fields, whose trajectories fall into a Banach space, but the domain of definition of the vector fields is a non-dense subset of the Banach space. The vector fields depend on a parameter that can assume different discrete values in two parts of the phase space and a continuous family of values on the boundary that separates the two parts of the phase space. In essence, the parameter parametrizes the possible vector fields on the boundary. The problem is to find one or more values of the parameter so that the solution of the PWS system on the boundary satisfies certain requirements. In this paper, we require continuous solutions. Motivated by the properties of applications, we assume that when the parameter is forced to switch between the two discrete values, trajectories become discontinuous. Discontinuous trajectories exist in systems whose domain of definition is non-dense. It is shown that under our assumptions the trajectories of such PWS systems have unique forward-time continuation when the parameter of the system switches. A finite-dimensional reduced-order model is constructed, which accounts for the discontinuous trajectories. It is shown that this model retains uniqueness of solutions and other properties of the original PWS system. The model reduction technique is illustrated on a nonlinear bowed string model.


Introduction
The purpose of model reduction is to extract the essence of a complex model, disregarding details that are irrelevant to a specific application. Depending on the question asked from the model, different kinds of model reduction are required. In many cases, only qualitative predictions are needed, where low-order analytically solvable models, such as normal forms used in bifurcation theory (Kuznetsov 2004), are useful. In other cases, the reduced-order model has to be solved numerically with a specified accuracy using constrained computational resources (Benner et al. 2017). Similar to model reduction, any numerical scheme that solves a continuum problem, such as finite elements, spectral collocation or finite differences, turns an infinite-dimensional continuous-time problem into a finite-dimensional problem. A numerical scheme, however, tends to emphasize quantitative accuracy, which might miss some qualitative features, such as differentiability of solutions. In this paper, we focus on the qualitative properties of solutions of piecewise-smooth (PWS) systems, with applications to numerical schemes and reduced-order models in mind.
For smooth systems, there are rigorous ways to obtain reduced-order models. Center manifold reduction (Carr 1981) about an invariant set, such as an equilibrium or periodic orbit, captures the slowest dynamics and can be used to study bifurcations, regardless of the dimensionality of the system (Kuznetsov 2004). In multiple timescale systems (Kuehn 2015) attracting slow manifolds that contain dynamics much slower than the rest of the system can be used to obtain reduced-order models. This paper discusses model reduction for infinite-dimensional systems that are piecewise smooth. The theory of PWS systems is summarized in Filippov (1998), which contains the basic definitions and results on existence of solutions in finite dimensions. There are numerous applications of PWS systems, where discontinuities are essential to the model or where rapid variations of the vector field over small regions of the phase space naturally lead to discontinuous approximations. Some applications in finite dimensions include neuron models with resetting (Coombes et al. 2012;Izhikevich 2003), DC-DC converters (di Bernardo et al. 1998), network dynamics (Danca 2002;DeLellis et al. 2015), friction oscillators (Oestreich et al. 1997;Szalai and Osinga 2008), gene regulatory networks (Glass and Kauffman 1973;Mestl et al. 1995) and so on. We consider the special case of differential equations that are discontinuous along a codimension-one hypersurface of their phase space, called the switching manifold. We assume that the phase space is a Banach space and that the domain of definition of the differential equation is not dense.
In contrast to smooth systems, center manifolds or slow manifolds that continue through switching manifolds do not exist for PWS systems. In general, the dynamics of singularly perturbed PWS systems cannot be reduced to an invariant manifold, because small-scale instabilities persist as the perturbation vanishes (Sieber and Kowalczyk 2010). For a special class of PWS systems, slow manifolds with similar properties to smooth systems exist (Fridman 2002;Cardin et al. 2013Cardin et al. , 2015. It is also possible to find equivalents of invariant manifolds which allow model reduction by considering the dynamics on the invariant manifold. Invariant cones can be found in systems with equilibria on the switching manifold (Weiss et al. 2012(Weiss et al. , 2015. Invariant polygons may also appear when an unstable focus-type periodic orbit interacts with discontinuities of the vector field (Szalai and Osinga 2008), which leads to periodic or chaotic dynamics (Szalai and Osinga 2009).
In infinite dimensions, the theory of PWS systems is focused on sliding mode control (Orlov and Utkin 1987) and PWS delay equations (Sieber 2006;Londoño et al. 2012). Sliding mode control applies a discontinuous control signal to a plant, in order to restrict the system onto an engineered hypersurface with a prescribed dynamics. The main objective of sliding mode control is to establish conditions that guarantee the prescribed dynamics. The results in this area concern systems that are densely defined on reflexive Banach spaces (Levaggi 2002a, b), which suggests that these systems are similar to finite-dimensional PWS systems.
In this paper, we relax the assumption of a dense domain of definition and not surprisingly we find different dynamics to what has been studied before. For this class of systems, we are able to prove uniqueness of solutions and also construct a reduced-order model. One consequence of the non-dense domain is the existence of discontinuous solutions, which is just the inverse of the Hille-Yosida theorem (Pazy 1983): trajectories of a linear autonomous system (as described by a semigroup) are strongly (or weakly) continuous if and only if the infinitesimal generator is closed and densely defined. The relevant mathematics describing our class of systems is the non-autonomous generalization of integrated semigroups (Neubrander 1988;da Prato and Sinestrari 1992). To illustrate that our class of systems are necessary to describe physical phenomena, we refer to McIntyre and Woodhouse (1979). In McIntyre and Woodhouse (1979), the authors have noticed that the measured impulse response function of a string has a discontinuity in the velocity component, which is manifest of the non-dense domain and that the initial condition is outside of the closure of the domain. This is shown later in the paper for the relevant mathematical model. Crucially, accounting for the discontinuity of the impulse response explains the observed asymmetric hysteresis of the stick-slip motion that causes 'flattening' of notes when the string is bowed in a certain way. The discontinuity of the impulse response is exactly the property that allows us to carry out model reduction.
The outline of the paper is as follows. We first carry out model reduction on a simple linear example to illustrate each step of the process, but without a rigorous justification of the steps. In Sect. 3, we review basic classes of PWS systems and highlight some cases where solutions may be non-unique. Section 4 describes the model reduction process in a general setting and shows that uniqueness of solutions and some other properties carry over to the reduced-order model. Section 5 describes a nonlinear example, the classical example of the bowed string, which highlights the significance that nonlinearity plays in the reduction process and uncovers some possibly surprising results that were not known about friction oscillators.

The Reduction Procedure Through an Example
To provide a straightforward template for the model reduction procedure, we take an idealized linear bowed string model with a single contact point and systematically apply our abstract procedure without rigorous justification. The list of steps is found at the start of Sect. 5. Let us consider the equation of motion of a linear bowed string u(ξ, t) = u (ξ, t), u(0, t) = u(1, t) = 0, u (ξ −, t) − u (ξ +, t) = λ, (2.1) where u(ξ, t) is the scalar-valued displacement of the string, t ∈ R is time and ξ ∈ [0, 1] is the spatial coordinate along the string; λ ∈ [− 1, 1] is the force applied at the contact point ξ , prime denotes differentiation with respect to ξ and dot with respect to t; u (ξ −, t), u (ξ +, t) denote the left and the right derivative at ξ , respectively. The value of parameter λ is given by where v 0 is the speed of the bow and λ represents the friction force between the bow and the string. In general, h is a smooth scalar-valued function of the state variables and is called the switching function. The equation h = 0 implicitly defines a surface in the phase space of (2.1) which is called the switching manifold (di Bernardo et al. 2008). Note that λ is not defined for h = 0 by Eq. (2.2). We assume that all values of λ ∈ [− 1, 1] are possible when h = 0, and therefore the model is a differential inclusion (Smirnov 2002). Later on, we will find a unique value for λ using the condition that the functions u(ξ , ·) andu(ξ , ·), that is, the solution of (2.1) and (2.2) evaluated at ξ = ξ must be continuous in time.
We now consider the case when λ is constant. For constant λ Eq. (2.1) has an equilibrium. The equilibrium shape of the string for λ = 1 is given by where H is the Heaviside function. Due to linearity, for a fixed λ the equilibrium is then λu 0 (ξ ). It is known that free vibrations of a string, that is, the solutions of (2.1) with constant λ can be written as u(ξ, t) = λu 0 (ξ ) + ∞ k=1 sin kπξ (a k sin kπ t + b k cos kπ t) , (2.4) where a k and b k are determined from initial conditions (Rao 2016, Sect. 8.2). We now consider solutions for which a k = b k = 0 for k > 1 in (2.4). The remaining two parameters a 1 , b 1 describe a set of solutions that are restricted to a two-dimensional invariant manifold, which we denote by M. For this set of solutions, we denote the displacement of the string at ξ by y(t) = λy + (a 1 sin π t + b 1 cos π t) sin πξ , where y is yet to be defined. The value y(t) can be used to recover the displacement of the whole string as u M (y(t), λ; ξ) = λu 0 (ξ ) + y(t) − λy sin πξ sin πξ .
(2.5) Expression (2.5) is the immersion of the invariant manifold into the configuration space, but not into the full phase space, which owing to the second-order time derivative in (2.1) should also contain velocities. Note that any value of y gives the same manifold, y only influences the parametrization of the manifold. To fix y , we require that the manifold does not move in the tangential direction of M when λ changes. This means that the two partial derivatives ∂ ∂λ u M (y, λ; ξ) and ∂ ∂ y u M (y, λ; ξ) must be perpendicular, that is, 1 0 ∂ ∂λ u M (y, λ; ξ) ∂ ∂ y u M (y, λ; ξ)dξ = 0.
Substituting the immersion (2.5) of manifold M into (2.1) and (2.2) while assuming that λ is constant, we getÿ where h = v 0 −ẏ. Equation (2.7) has the form of a common PWS system, which is widely used as a reduced-order model of (2.1). However, the assumption that λ is constant does not hold when h = 0, and therefore we consider (2.7) a skeleton of a more accurate description and call (2.7) the skeleton model. The skeleton model (2.7) is a typical friction oscillator and therefore can be solved using techniques known from mechanics. At stick, when h = 0, λ is not explicitly defined by the skeleton model (2.7), instead we need to argue the following. If h = 0 on some interval of time, thenẏ = v 0 on this interval, and consequently,ÿ = 0. Substituting the conclusion of this argument into the first line of (2.7) yields λ = y y . (2.8) We only allow λ ∈ [− 1, 1], and therefore if the result of (2.8) is outside of the interval [− 1, 1], λ simply swaps from 1 to − 1 or vice versa. The argument made to find (2.8) is equivalent to Filippov's closure, which is summarized in Sect. 3. The phase portrait of the skeleton model (2.7) can be seen in Fig. 1. We focus on the dynamics at stick, which occurs on the switching manifold, highlighted by the horizontal red-shaded plane. The dashed red lines correspond to discontinuities in λ. The solid green line on the horizontal red-shaded plane is the stick solution, where the friction force λ grows with a constant rate with respect to t and y until it reaches the limit λ = ± 1 and slip ensues.
By assuming constant λ, we made an error when the relative velocity between the bow and the string becomes zero, i.e.,ẏ = v 0 , because contrary to the assumption λ is not constant and can even jump between − 1 and 1 or from ± 1 to the value given . Solid lines represent solutions where λ is continuous, and dashed lines represent solutions where λ switches between two values. The horizontal plane withẏ = v 0 contains sticking solutions. There is only one sticking solution represented by the solid line, which is also marked with an arrowhead denoting the direction of time by (2.8). The desire to correct for this error is the subject of the paper, because this is the source of the qualitative discrepancy between solutions of Eqs. (2.1), (2.2) and the skeleton model (2.7). To account for the error made, the solution of Eqs. (2.1) and (2.2) is now written as (2.9) where w(ξ, t) is a correction term. Assuming that y(t) satisfies (2.7) and substituting (2.9) into (2.1) without assuming constant λ, we get the governing equation of the correction termẅ (2.12) The difficulty of evaluating Eq. (2.12) lies with solving (2.10), which we carry out in Appendix A in detail. We now illustrate how a discontinuity of λ leads to a jump in the velocityẇ(ξ , t). For this, we assume that λ(t) = H (t) in Eq. (2.10). The solution of Eq. (2.10) can be seen in Fig. 2, with initial conditions (2.11). By comparing Fig. 2b, d, it can be seen that the velocityẇ(ξ , t) at t = 0 has a discontinuity, whose gap is proportional to the jump in λ. The time history of this velocity in Fig. 2j has further discontinuities. Discontinuities for t > 0 are due to reflections at the boundaries, and they are specific to this example that lacks damping. Typically, wave dispersion or damping that is present in other mechanical systems would destroy discontinuities for t > 0, but not which we call the normal discontinuity gap. Due to the linearity of Eq. (2.10), any jump in λ is translated into a discontinuity of the velocityẇ(ξ , t). This velocity jump also appears in the switching function (2.12), which makes a qualitative difference between the dynamics of the original model and the skeleton model (2.7) on the switching manifold as we show below.
After solving Eq. (2.10) on the interval 0 ≤ t < min (2ξ , 2 − 2ξ ), before any discontinuity is reflected back to ξ = ξ , we find that the switching function (2.12) in Eq. (2.7) becomes where κ is a variable satisfying the initial value problem Note that by using h as defined by (2.14) in Eq. (2.7), we get an exact representation of the dynamics of the original problem (2.1) and (2.2) on the time interval 0 ≤ t < min (2ξ , 2 − 2ξ ). The valid time interval can be extended to any length by considering the full solution of (2.10) derived in Appendix A. The switching function (2.14) depends on λ, because of the presence of the normal discontinuity gap (2.13). Therefore, when h = 0, λ can be solved for, without any closure rule, such as Filippov's or Utkin's (see Sect. 3). In our case, solving the equation A nonzero normal discontinuity gap, as calculated in Eq. (2.13), turns the skeleton model (2.7) at stick into an index-1 differential algebraic equation. When gathering all dynamic equations at stick, we geẗ (2.17) By definition, an index-1 differential algebraic equation can be turned into an ordinary differential equation by differentiating the algebraic Eq. (2.16) once, which for Eq.
(2.17) of the bowed string at stick givesÿ Note that the differentiation of (2.16) also transforms the stick constraint h = 0 intȯ h = 0, and therefore (2.18) is valid for initial conditions that satisfy h = 0 with λ ∈ [− 1, 1]. We can now put together the whole model with the correction into a single system Phase portrait of the full model (2.1) and (2.2). The simulation was carried out using an exact representation of the switching function h, a further extension of (2.14) to t ≥ 0, which is derived in Appendix A We call Eq. (2.19) the reduced-order model of the initial problem (2.1) and (2.2), because it exactly reproduces the dynamics for initial conditions in M and for the time interval 0 ≤ t < min (2ξ , 2 − 2ξ ). Appendix A shows that the valid time interval can be extended to any length by including delayed values of λ in the switching function. It is noteworthy that instead of the two-phase space regions defined by (2.2), the reduced-order model (2.19) has three regions, where the dynamics is defined. The additional phase space region corresponds to the stick phase of motion, which has its own regular dynamics. This dynamics follows from the assumption that the velocity y 2 is continuous in time and there is a normal discontinuity gap, i.e., the correction term w is discontinuous at t = 0, when λ = H (t).
The phase portrait of the reduced-order model (2.19) can be seen in Fig. 3. In the simulation, we have used the result of Appendix A to extend the valid time interval to an appropriate length. In comparison with the skeleton model (2.7) shown in Fig. 1 the dynamics at stick becomes more complicated. The dynamics in slip is the same, because λ is constant and decoupled from the rest of the variables due to the choice of immersion (2.6). The stick dynamics is now described by the differential Eq. (2.18), and therefore there is no discontinuity of λ. Due to the higher dimensional dynamics that arise from the inclusion of κ as a dynamic variable and delayed values of λ, the dynamics depicted in Fig. 3 is only a projection. Regardless of the differences, the phase portrait in Fig. 3 appears as a smoothed version of the same dynamics in Fig. 1, even though no smoothing or regularization was applied. Furthermore, to solve the reduced-order model (2.19) we did not need an arbitrary closure, such as Filippov's to define the dynamics at stick, instead the solution followed straight from the initial problem (2.1) and (2.2).
In Sect. 4, we explore a generalization of Eqs. (2.1) and (2.2). We consider models whose solutions may be norm-discontinuous as illustrated by the linear string model. Before we embark on the general theory, we recall basic definitions and properties of PWS models.

Finite-Dimensional PWS Models
In this section, we summarize two commonly used closures of PWS systems. As it turns out, these common PWS systems are special forms of the skeleton model to be defined in Sect. 4.2. An introduction to the state of the art can be found in Glendinning and Jeffrey (2017); however, the book of Filippov (1998) contains the most general definitions of PWS systems. Below, we review the cases defined in (Filippov 1998, Chapter 2, §4), which are used most commonly in applications. We avoid cases where the vector field is a set-valued function (Filippov 1998, Chapter 2, §5, §6). We also limit the description to the bi-modal case, where the discontinuity occurs along a single implicitly defined manifold in the phase space.
Note 3.1 In addition to various notation for derivatives, in what follows D is also used to denote the Frechet derivative of a function; a subscript of D, such as D k denotes the partial derivative with respect to the kth argument of a function and a superscript such as D j k denotes the jth derivative with respect to the kth argument. Let us consider the vector fielḋ G is a compact and connected subset of R n and n, p ∈ N + . The function h ∈ C p (G, R) is called the switching function, and its zero set defines the switching manifold A solution y : I → R n of Eqs. (3.1) and (3.2) is defined on a closed interval of nonzero length I ⊂ R. There is no information in Eqs. (3.1) and (3.2) that helps to deduce a value for λ on Σ. To explore all possibilities (3.1) can be turned into a differential inclusionẏ ∈ co f ( y, where co denotes the closure of the convex hull of a set. The existence of solutions of (3.6) is investigated in Filippov (1998, Chapter 2, §7). In this section, we review different definitions of λ on Σ.
We have assumed two possibilities, (3.3) or (3.4), for the domain of definition of f . The case of (3.3) is the minimum necessary to make Eqs. (3.1) and (3.2) consistent. In many applications, such as mechanics, the larger domain of definition (3.4) is naturally given, which is useful to define the solutions of (3.1) and (3.2) on Σ as we show later in this section.
The system (3.1) and (3.2) has a unique solution on an interval of nonzero length for initial condition y(0) ∈ G if h( y(0)) = 0, because f is a smooth vector field (Coddington and Levinson 1955). However, for h( y(0)) = 0 the vector field is not defined and one needs to reason how trajectories continue once they reach Σ.
The simplest case of a trajectory interacting with Σ is when the trajectory approaches Σ transversely from one side and continues on the other side; this is called crossing. The exact condition for crossing is that the value of h must change monotonically through h = 0 with a nonzero speed along the trajectory at h = 0. If the common point of the trajectory with Σ is denoted by y , the speeds at which h increases are Dh( y ) f ( y , ± 1). h changes monotonically through h = 0 with a nonzero speed if and only if (3.7) In case of (3.7), there is no need to define the dynamics on Σ because the value of λ simply switches between ± 1. In our argument, we have used f for λ = ± 1 only, therefore, to resolve crossing, it is sufficient to assume (3.3). The subset of Σ, where (3.7) holds, is called the crossing region and denoted by Σ cr . Now we discuss the case when When (3.8) holds, Σ is an attractor in either forward or backward time. This means that trajectories cannot immediately escape Σ once they are on Σ. The subset of Σ, where (3.8) holds and attracts solutions in forward time is called the sliding region and denoted by Σ sl . The repelling subset of Σ, where (3.8) holds, is called the escaping region and denoted by Σ esc . Equations. (3.8), (3.1) and (3.2) are not sufficient to define a solution, and an assumption is required that specifies how a trajectory continues on Σ. In this paper, we call such an assumption a closure, because it completes (3.1) and (3.2). In the following, we discuss two commonly used closures. The first closure is attributed to Filippov (1998, Chapter 2, §4, 2.a), the second closure is due to Utkin (1992), and also explored in Filippov's book (Filippov 1998, Chapter 2, §4, 2.b). We note that there is no common terminology in the literature for closures, various closures have their own name. For example, Filippov's closure is commonly called Filippov's method and Utkin's closure is called the equivalent control method, due to its origin in control theory. There are many possibilities to define a closure, for example, Filippov explores systems where the closure is not explicitly defined, but only certain constraints are placed on it Filippov (1998, Chapter 4
It can be shown that λ ∈ (− 1, 1) when condition (3.8) holds (Filippov 1998). Fillipov's closure is illustrated in Fig. 4c, which shows that the vector field given by (3.9) and (3.11) is chosen from all convex combinations of f ( y, ± 1) so that r( y) + b( y)λ is tangential to Σ. A trajectory at its first point of contact with Σ is continuous, but not continuously differentiable, because λ becomes discontinuous due to (3.11). When neither (3.7) nor (3.8) holds for y ∈ Σ, we have Equation (3.12) means that one or both of the vector fields f ( y, ± 1) is tangential to Σ, which we call a tangency. The boundaries of crossing, sliding and escaping regions are formed by tangencies, which generally occur as codimension-one surfaces of Σ. Trajectories may not have unique continuation when they are tangent to Σ. When both vector fields f ( y, ± 1) are tangential to Σ at a point, λ is not uniquely defined by (3.11) as both the numerator and denominator of (3.11) vanish. Consequently, the forward-time solution of (3.1), (3.2) and (3.11) is not unique. This case is illustrated in Fig. 4d, which shows a set of possible directions that a solution can follow. A particular case of this double tangency is the Teixeira singularity, where an open set of initial conditions generate trajectories that go through the double tangency. The Teixeira singularity (Teixeira 1982) was studied extensively (Colombo and Jeffrey 2011;Filippov 1998;Kristiansen and Hogan 2015;Szalai and Jeffrey 2014) in various contexts.

Utkin's Closure
In this section, we assume that the domain of definition of f is given by (3.4). In the case of (3.8), we similarly construct the vector field on Σ, such that Σ becomes invariant under the vector field. The invariance of Σ is expressed as (3.13) Equation (3.13) has at least one solution for some λ ∈ [− 1, 1], because (3.4) implies that Dh( y) · f ( y, ± 1) has different signs and due to Bolzano's theorem there must be a root.
The root of (3.13) may not be unique, which renders the solution of (3.1) and (3.2) non-unique. We also note that (3.13) can have a solution even when (3.7) holds in the crossing region. A simple case of Utkin's closure is illustrated in Fig. 4a. The green curve connecting f ( y, 1) to f ( y, − 1) represents the possible values of the vector field on Σ. There is one intersection of this family of vectors with the tangent plane of Σ, represented by the thick red arrow, which satisfies Eq. (3.13). Figure 4b shows that there can be multiple intersections of f ( y, [− 1, 1]) with the tangent plane of Σ, that then yields multiple solutions. Note that in the case of Fig. 4b, Filippov's closure yields a unique solution. The contrary, when Utkin's closure predicts a unique solution and Filippov's closure predicts a family of solutions, is also possible. For example, when the convex hull represented by the green dashed line in Fig. 4d is deformed slightly, the possible number of solutions can be reduced to three. Out of these three solutions, there is only one with λ ∈ (− 1, 1).

Model Reduction
We start with a general continuum model, in the form oḟ The switching function h is defined on ∪ λ D λ (F) and has values in R. When h = 0, the most general definition of the dynamics isẋ ∈ coF(x, [− 1, 1]). We also require that trajectories are continuous, even when h = 0. The smoothness of F and h is not assumed globally, instead we assume the smoothness of an invariant manifold of F and related quantities in the next section.
Remark 4.1 The notation of Eq. (4.1) facilitates that λ is an unknown, which needs to be found when h(x) = 0. Therefore, λ may not be a function of x, but it may become part of the phase space. The solution for λ, when h(x) = 0 is defined in Sect. 4.3. This is a similar setting to Sect. 3, except that the phase space is now infinite dimensional, and therefore a different kind of solution is required for λ.

The Invariant Manifold
For PWS systems, such as Eq. (4.1), differentiable invariant manifolds that extend over switching boundaries do not exist. This fact makes model reduction more complicated than for smooth models. It is however possible to find invariant manifolds for constant λ of the PWS system (4.1). Our approach is therefore to first consider the smooth systemẋ We make the following initial assumptions and definitions: (A1) Existence of an invariant manifold. We assume that there exists a function which satisfies the invariance condition where G is a compact and connected subset of R n . The invariant manifold is given by and the dynamics of (4.2) on M λ is described bẏ This derivative is denoted by We also assume that the domain of definition of A 1 , i.e., D( A 1 ( y, λ)) = {x ∈ X : A 1 ( y, λ)x ∈ X}, is independent of y and λ, and we define Z = D( A 1 ( y, λ)) (in general, D( A 1 ( y, λ)) = D λ (F)). (A3) Unique continuous solutions. We assume that the abstract Cauchy probleṁ We also assume that the Z component of the solution can be written as where K is bounded for τ ≥ s and continuous in both variables for τ > s. The underlying conditions of existence of unique solutions can be found in da Prato and Sinestrari (1992). For discussion, see Remarks 4.3 and 4.4.
(A4) M λ is attracting and normally hyperbolic. We assume that there exist two families of projections Π c ( y, λ) and Π s ( y, λ), strongly continuous in y and λ such that Consider the non-autonomous ordinary differential equationη = D 1 f ( y, λ)η, whose solutions with initial condition η 0 at t = s are denoted by η(t, s, η 0 ). We assume that there exist real numbers σ c > 0, σ s < −σ c ,M c > 0 and M s > 0 such that (4.10)

Remark 4.2
For systems with an equilibrium it is natural to consider spectral submanifolds (Haller and Ponsioen 2016), that are the smoothest invariant manifolds tangent to an invariant linear subspace of the variational problem about the equilibrium. The uniqueness and existence of such manifolds are established in Cabré et al. (2003). In order to be meaningful, these manifolds need to contain the slowest dynamics within the system to capture long-time behavior. This requirement is outlined in points R1 and R2 of Haller and Ponsioen (2017).

Remark 4.3
We do not fully specify the definition of a solution of (4.1) and (4.6) apart from the solution being continuous. The results of this paper only depend on the form of the solution as given by (4.7) and not how it is obtained. However, it might be helpful to think of F-solutions of (4.6) as defined by da Prato and Sinestrari (1992):

Remark 4.4
The existence of unique F-solutions of Eq. (4.6) is established in da Prato and Sinestrari (1992) in Theorem 5.1. However, for many examples, e.g., elastodynamics (Gurtin and Sternberg 1961;Martins and Oden 1987) and delay equations (Diekmann et al. 1995), existence and uniqueness results are already known and it is not necessary to check the conditions listed in da Prato and Sinestrari (1992). The existence and regularity of a convolution kernel for non-autonomous problems are not discussed in the literature. However, the autonomous problem is discussed in Thieme (1990Thieme ( , 2008, which implies that the convolution integral is because (μ − A 1 ( y(τ ), λ(τ ))) −1 : X → D. da Prato and Sinestrari (1992) uses a similar technique to approximate the unique solution. The kernel K , however, has two parameters, and therefore its smoothness properties are not trivial even if we know that the convolution (4.11) is continuous in t. We have therefore assumed the continuity of K for t > s, which allows for a discontinuity at t = s due to D 2 W ( y(τ ), λ(τ )) / ∈ Z.

Remark 4.5
The uniqueness or persistence of M λ is not addressed by the assumptions. For persistence of M λ under a perturbation, additional smoothness conditions on the solutions of (4.2) have to hold, which can be found in Bates et al. (1998).

Remark 4.6
The condition (4.10) implies that the convolution in (4.7) remains bounded whenλ is bounded. This will be useful later when the reduced-order model is constructed. Figure 5 shows the invariant manifold M λ and its intersection with the switching manifold Σ. The Banach space X is represented by two coordinates x 1 and x 2 . The two parts of the invariant manifold (M −1 and M 1 ) do not join up in Fig. 5a. If trajectories cross Σ instantaneously, they are discontinuous. When discontinuity is not allowed, the crossing cannot be instantaneous. In certain cases, however, the disconnected nature of M λ may be overlooked. For example, when the switching function solely depends on the parameter y of the immersion W , i.e., y = x 1 in Fig. 5. The case when the dynamics is restricted to M λ is discussed in Sect. 4.2. Figure 5b shows the extended phase space and how solutions of (4.1) behave about M λ , when instantaneous crossing is not allowed. In Fig. 5b, M λ is a connected manifold. When a solution of (4.1) arrives at Σ, the value of λ must change, so that a trajectory can enter Σ. M λ is only invariant for constant λ, and therefore a trajectory (denoted by dotted lines) will not continue on M λ while also in Σ. Once a trajectory has left Σ, it will be attracted to M λ as per Assumption (A4).
In the following sections, we discuss how the departure of a trajectory from M λ can be captured and whether or not capturing this dynamics makes a qualitative difference in the predictions of the model. In Sect. 2, we have already seen that including a correction that captures the departure from M λ makes a difference and trajectories can no longer cross Σ instantaneously.

The Skeleton Model
Having assumed the existence of an invariant manifold M λ , it is natural to consider the dynamics on M λ in the presence of switching. This can be done by substituting the immersion W into the full problem (4.1) and disregarding that λ may not be constant on Σ. We start with the switching function (4.12) In contrast to Sect. 3, the switching function (4.12) depends on λ, and therefore the closures described in Sect. 3 may not apply. Using the vector field (4.5) on M λ and (4.12), we obtainẏ Definition 4.7 Equation (4.13) is called the skeleton model of (4.1) on the invariant manifold M λ .
Definition 4.7 alludes to what follows next. We will use Eq. (4.13) to build upon and not consider it as an end result. Equation (4.13) is inaccurate when λ varies and that causes solutions to become non-unique, even if they were unique in the full problem (4.1). Nevertheless, we highlight some properties of the skeleton model that carry over to the reduced-order model. We note that already in the skeleton model the switching function h 0 can become dependent on λ. This means that the dynamics when h 0 = 0 may be defined as an index-1 differential algebraic equation. To describe such dynamics, in the introductory example in Eq. (2.19) we needed to separate the switching manifold into two components. Here, we formalize this splitting and define two new switching manifolds (4.14) In the extended state space ( y, λ) ∈ G × [− 1, 1], Σ ± 0 is the boundary of the ndimensional manifold (4.15) Σ ± 0 cannot intersect each other in the extended state space.
When h 0 ( y, λ) = 0, the trajectories are described by the vector fielḋ Otherwise, we have an index-1 differential algebraic equatioṅ ( 4.17) A unique solution to (4.17) is guaranteed by the implicit function theorem if D 2 h 0 ( y, λ) = 0, so that there is a unique C p smooth function λ( y) satisfying h 0 ( y, λ( y)) = 0. This also implies that λ(t) = λ( y(t)) is continuous, and hence there is no discontinuity of λ when a trajectory reaches Σ ± 0 transversely. Trajectories must spend nonzero time on Σ 0 in order to keep λ continuous. This short argument highlights a major difference between PWS models described in Sect. 3, where we have D 2 h( y, λ) = 0 and where the Implicit Function Theorem does not apply. Models in Sect. 3 are special cases of the skeleton model.
We can also write the index-1 differential algebraic Eq. (4.17) in a differential form by differentiating the constraint h 0 ( y, λ) = 0, that is, As discussed, whether solutions are well defined, depends on the term If (4.19) is nonzero, Eq. (4.18) can be solved forλ, which yields the differential form of (4.17) for ( y, λ) ∈ Σ 0 , that is, (4.20) The continuous concatenation of solutions of Eqs. (4.16) and (4.20) gives the full solution of Eq. (4.13). This concatenation is a PWS problem, where Σ ± 0 are now separating the phase space into three regions. The following theorem looks at the case when there is no need to define the dynamics on Σ ± 0 .
Theorem 4.8 Consider a point ( y , λ ) ∈ Σ ± 0 and assume that Trajectory T has a unique continuation for t > 0 (or t < 0) sufficiently small as a solution of the skeleton model (4.13) if one of the following conditions holds: 0 and the order of the tangency is less than the smoothness order (C p ) of h 0 . In other words, there exists 0 < ≤ p such that Proof The proof can be found in Appendix B.
Remark 4.9 Theorem 4.8 excludes the case D 2 h 0 ( y, λ) > 0. For D 2 h 0 ( y, λ) > 0, transverse trajectories (case 1 of Theorem 4.8) cannot cross Σ ± 0 . Tangential trajectories with even may have multiple continuation, which is the case of the Teixeira singularity (Colombo and Jeffrey 2011). Tangential trajectories with odd cannot cross Σ ± 0 , similar to transverse trajectories. To investigate the case of D 2 h 0 ( y, λ) > 0 in detail, a definition of how trajectories move along Σ ± 0 (with λ = ± 1) is also required, which falls outside of the scope of this paper.

Dynamics About Manifold M Due to Switching
This section describes a correction to the skeleton model (4.13) that resolves the dynamics in the neighborhood of M λ up to linear order. The correction is necessary, because theλ = 0 assumption does not hold: Eq. (4.20) states that λ varies on Σ 0 . The correction that is introduced here captures trajectories that depart from M λ when h = 0 (see dashed line in Fig. 5b).
Let us suppose that where z represents the difference between the trajectories of the full model (4.1) and the skeleton model (4.13). This setup is illustrated in Fig. 5a. To derive an equation for z, we substitute (4.23) into (4.1), while taking into account that λ is a function of time. This substitution yieldṡ We assume that z is a small deviation from M λ and Taylor expand F(W ( y, λ) + z, λ) in z about z = 0, that is, The expansion (4.25) when substituted into (4.24) yields (4.26) We now use the invariance Eq. (4.3) and the dynamics on M λ as given by (4.5) and notice that two terms cancel in (4.26), so that we get ( 4.27) Combining the skeleton model (4.13) with (4.27) yields the corrected model . A unique solution of (4.28) is assumed in (A3) with a continuously differentiable λ. In this paper, we do not investigate whether the corrected model (4.28) is a faithful representation of the fully nonlinear system (4.1); for some discussion, see Remark 4.14.
We define the switching manifolds as When a trajectory is restricted to Σ, the solution must satisfy Similar to the skeleton model, we evaluate how h changes in time and we restrict this change to zero on Σ to find an equation for λ [cf. Eq. (4.18)]. To evaluate Eq. (4.29), we use the following lemma.
Lemma 4.10 Assume (A3) and that λ is continuously differentiable and y, z satisfy the differential equations Then the right-side derivative of h as a function of time is calculated as The proof can be found in Appendix C.

Remark 4.11
The quantity d ± ( y, z, λ) in (4.31) measures the discontinuity of the convolution kernel K at t = s. A discontinuous K is possible, because D 2 W ( y, λ) ∈ X\Z, and the continuity Assumption (A3) does not apply at t = s. Such a discontinuity allowed us to find a differential equation for λ in Sect. 2.
We also define two other quantities that will be useful later. These are and therefore we have the identity d ± ( y, z, We now find the governing equation of the dynamics on Σ. We solve equation (4.34) The trajectories of Eq. (4.34) are concatenated with trajectories oḟ along the boundaries Σ ± and form the trajectories of the corrected model (4.28).
The following theorem provides a sufficient condition for a unique continuation of trajectories through Σ ± .
Theorem 4.13 Assume (A1)-(A5). A trajectory T of either (4.34) or (4.35) with an end point ( y, z, λ) ∈ Σ ± at t = s has a unique continuation for t > s with t −s sufficiently small, as a solution of the corrected model (4.28), if the following conditions hold: s)z is continuously differentiable with respect to t for t ≥ s and 3. one of the vector fields, (4.34) or (4.35) is not tangent to Σ ± , that is, (4.37) Proof The proof of Theorem 4.13 can be found in Appendix C.
Remark 4.14 The linear correction about the invariant manifold is carried out here without an assessment whether trajectories of the corrected model (4.28) and the full model (4.1) are qualitatively the same. If z 1, the linear correction is accurate. Because on M λ , we have z = 0, when a trajectory enters Σ, the rate of change of z is determined byλ. The magnitude ofλ depends on the f and d ± . Smaller d ± makes λ faster. The value of d ± is not necessarily a small parameter, and therefore the deviation from M λ can stay small. For the linear string d ± = 1 2 . In the literature of regularized PWS systems (Jeffrey 2015;Kristiansen and Hogan 2015), to stay close to the skeleton model, fast λ is assumed.

Remark 4.15
If d ± ( y, z, λ) = 0, the dynamics about M λ as captured by variable z can only have a second-order effect on h due to the nonlinearity of h. Therefore, (4.30) is independent ofλ and d dt + h = 0 cannot be solved forλ. When d ± ( y, z, λ) = 0, the corrected model (4.28) needs a closure, such as Filippov's or Utkin's. d ± ( y, z, λ) = 0 occurs when U is strongly continuous on the whole of X, i.e., Z = X. This case for linear systems is explored in Orlov (1995) and Levaggi (2002a, b).

Remark 4.16
The transversality condition (4.37) is the equivalent of case 1 of Theorem 4.8. The equivalent of case 2 of Theorem 4.8 is not proven here, but a similar argument can be made while carefully accounting for the infinite-dimensional nature of the problem.

Remark 4.17
It is possible to consider a nonlinear correction, so that (4.23) becomes exact. Let us define the nonlinear term without discussing the constraints on N. The equation of the exact correction can be written asẏ (4.38) Equation (4.38) is a semi-linear abstract Cauchy problem, which is frequently analyzed in the mathematical literature. The solutions of (4.38) are formally obtained from the integral equation (4.39) In general, existence and uniqueness of solutions of (4.39) are established using a contraction mapping argument. However, under our assumptions the convolution is not justified because U(t, s) is only defined on Z, but (4.40) Regardless of (4.40), the autonomous case (Thieme 2008;Magal and Ruan 2009) has unique solutions under appropriate conditions. The author is confident that a similar argument can be made to establish unique solutions (4.39) although that might require that the nonlinearity N( y, λ; ·) : D → X be bounded.

Time-Scale Separation
We already have some indication that switching has a great influence on the normal dynamics. For example, ignoring the normal dynamics as in the skeleton model (4.13) leads to a different uniqueness condition than for the corrected model (4.28). In this section, we restrict the analysis to the simplest case where there is a separation of time scales. We assume a parameter 0 ≤ ε ≤ 1 and denote the dependence on ε by a subscript, that is F ε . Here, the ε = 0 limit is represented by the skeleton model (4.13) and ε = 1 refers to the corrected model (4.28). Naturally, the immersion W ε ( y, λ) of the invariant manifold also depends on ε, which implicitly assumes that M λ persists for 0 ≤ ε ≤ 1. Whenever we write F 0 or W 0 , we mean the ε = 0 limit. Let us define the scaled Frechet derivative as A ε ( y, λ) = ε D 1 F ε (W ε ( y, λ), λ).
With this notation the corrected model (4.28) becomeṡ Changing the time scales by introducing t = εθ we get where˚stand for d /dθ. When setting ε = 0 we arrive at the layer system which stipulates that variable y is constant along trajectories. We assume the following: (A3) Assumptions (A3) holds when (4.6) is replaced by (4.42) for all ε ∈ [0, 1]. The unique solution of (4.42) can be written as (A4) Assumption (A4) holds when (4.6) is replaced by (4.42) and (4.44)

Remark 4.18
As a consequence of (A3) and (A4), A 0 ( y, λ) has an n-dimensional kernel spanned by D 1 W 0 ( y, λ), and Π c ( y, λ)A 0 ( y, λ) = 0. Because of (A5) and for t ≥ s, we also have We investigate the non-smooth dynamics for ε = 0. The case of constant λ is trivial, because we have assumed that M λ is attracting for 0 ≤ ε ≤ 1. Next we consider the dynamics in Σ, which is described bẙ (4.45) Any point in M λ ∩ Σ, i.e., y ∈ G, z = 0, λ ∈ [− 1, 1] is an equilibrium of (4.45); therefore, M λ is invariant under all the dynamics for ε = 0. It is however not obvious whether M λ ∩ Σ is attracting for ε = 0, which is addressed by the next theorem. Proof In order to calculate whether the critical manifold is attracting, we linearize Eq. (4.45) by using λ = λ 0 + α as a perturbation The initial conditions α(0) and z(0) are linked through Eq. (4.50), such that It is sufficient to show that α decays, because by assumptions (A3) (A4) and without forcing, the z component decays to a constant; if the initial condition satisfies Π s ( y, λ)z(0) = 0, z decays to zero. Applying the Laplace transform to (4.49), we find that where s is the Laplace parameter. By substituting (4.51) into (4.50), we find which can be rearranged into (4.52) The asymptotic properties of α(t) are determined by the poles of (4.52). The poles of the numerator are already given by the spectrum of A 0 ( y 0 , λ 0 ), which is assumed to be in the left half of the complex plane because M λ is attracting for constant λ. Therefore, only the roots of the denominator can cause instability, and hence the condition that (s) has roots in the left half of the complex plane is sufficient.

Remark 4.20
The proof can be extended to calculate the initial and final values of α.

Remark 4.21
Similar to Remark 4.5, normal hyperbolicity does not imply the persistence of M crit λ under variations in ε. The theorem of Bates, Lu and Zeng Bates et al. (1998) suggests that the evolution operator U needs to be differentiable (among other conditions) for M crit λ to persist for small ε > 0. Note that the nonlinear string example in Sect. 5 generates such a differentiable U on Z.

Remark 4.22
When both regions of M λ , that is M λ ∩ Σ and M λ \ (M λ ∩ Σ), persist for ε > 0, they most likely become discontinuous at the boundaries Σ ± , and hence as a whole, M λ does not persist. This is because the vector fields are discontinuous. Therefore, for ε > 0, trajectories that followed one part of M λ must jump to the other part of M λ , which induces fast transients that we are unable to characterize under general settings.

Qualitative Approximation of Normal Dynamics and the Reduced-Order Model
A key difference between the skeleton model (4.13) and the corrected model (4.28) is that they have unique solutions under different conditions. This difference is caused by the fact that the skeleton model does not take into account the normal discontinuity gap d ± . To rectify the omission of d ± , the skeleton model is extended by a scalar variable, which represents the dynamics of the convolution kernel K in Eq. (4.7). We call this extension the reduced-order model. It is then shown that the reduced-order model reproduces uniqueness of solutions and the existence of a critical manifold under equivalent conditions to those of Theorems 4.13 and 4.19.
To simplify the ensuing analysis, we assume that where Dh is a constant linear functional.
Assumption (A6) allows us to derive a scalar representation of z(t) without worrying about a varying Dh(x). The switching between parts of the state space depends on In what follows, we approximate the scalar-valued Dh · z in (4.53) by a convolution integral. Combining Eqs. (4.7), (4.11) and z(0) = 0 yields In order to proceed, either (A5) or time-scale separation with (A5) can be assumed. When (A5) is assumed, we are restricted to use ε = 1 and if (A5) is assumed we set σ = σ s . Now we can approximate that which neglects trajectories in the normal bundle of M λ that are decaying with exponents smaller than ε −1 σ . Note that By simply defining d + ( y, λ) = d + ( y, 0, λ) we get (4.54) After defining κ = Dh · z(t), we find that the approximation (4.54) satisfies the differential equationκ with initial condition κ(0) = 0. The switching function (4.53) using the new variable κ becomes We can also redefine the switching manifolds With this notation, the skeleton model extended with the approximate normal dynamics becomesẏ (4.57) Definition 4.23 We call Eq. (4.57) the reduced-order model of (4.1).
When h ε ( y, κ, λ) = 0, the dynamics of κ is decoupled from the rest of the variables and κ exponentially vanishes, because σ < 0-due to Assumption (A5) or (A5). When h ε ( y, κ, λ) = 0, we apply the same technique as in Sect. 4.2 to find a differential equation for λ. We express that is defined by (4.32). We drop the ( y, λ) arguments and solve (4.55) and (4.58) forκ andλ to arrive aṫ which governs the dynamics on Σ ε . We can now check that the reduced-order model (4.57) has the same key properties as the corrected model (4.28). In what follows, we outline the equivalents of Theorems 4.13 and 4.19 for the reduced-order model (4.57).

Proposition 4.24 A trajectory T of the reduced-order model (4.57) with an end point at ( y , κ , λ ) ∈ Σ ±
ε has a unique continuation through ( y , κ , λ ) if 1. d ± ( y , 0, λ ) > 0 as defined by Eq. (4.31) and 2. when trajectory T is not tangent to Σ ± or trajectory T is tangent to Σ ± ε and the of order of the tangency is not greater than the smoothness order (C p ) of h ε , that is, there exists 0 < ≤ p such that d dt h ε ( y(t), κ(t), λ )| y= y ,κ=κ = 0.
Proof The proof is the same as for Theorem 4.8 if we replace ( y, κ) → y and d ± → −D 2 h 0 ( y, λ).

Proposition 4.25 Let
be a compact set for ε = 0. Then, M crit λ is an attracting critical manifold of Eq. (4.59) which persists for a sufficiently small ε > 0. The dynamics on the critical manifold is governed by the skeleton model (4.20).
While κ = 0 on the critical manifold, lim t→∞ lim ε→0 ε −1 κ(t) may not be zero, that is, the limits t → ∞ and ε → 0 do not commute. After introducing εγ = κ, we can write thatẏ Setting ε → 0 and some algebraic manipulation yieldṡ which is the same equation as (4.20) of the skeleton model.

Remark 4.26
We know that σ < 0, because of Assumption (A5) or (A5). If Proposition 4.24 also holds, M crit λ is attracting when d − < 0. This is the same condition under which solutions of the skeleton model (4.13) are unique due to Theorem 4.8.
Next, we investigate in what sense the reduced-order model (4.57) is similar to the corrected model (4.41) with time-scale separation. It turns out that on Σ ε the critical manifold is likely to be attracting or repelling under the same conditions. The precise statement is in the following proposition.
Proof Because of the assumption d ± ( y, 0, λ) > 0, the stability of an equilibrium of (4.61) purely depends on d − , i.e., the equilibrium is attracting when d − < 0. On the other hand, substituting s = 0 into (s) as given by (4.47), we note that (0) = d − ( y, λ). This means that we have a zero root of (s) when d − = 0. Next, we show that this zero root of (s) becomes of the same sign as d − as y(α), λ(α) changes along γ . Let us now assume that at α = 0, we have d − ( y(0), λ(0)) = 0 and denote the root of that smoothly depends on α by s : (− δ, δ) → R and for which s(0) = 0. We denote the derivative with respect to α by and calculate the derivative of s from the definition (4.47), that is, where we omitted that y, λ are evaluated at α = 0. We also calculate the derivative d − ( y, λ) = D 1 d − ( y, λ) y + D 2 d − ( y, λ)λ and notice that the derivatives s and d − have the same sign when Dh · A −1 0 ( y, λ)D 2 W ( y, λ) > 0, which proves the proposition.

Remark 4.28
In the next section, for the example of the nonlinear string, d − is a small parameter, which measures how well the equilibrium shape of the string is approximated by a truncated Fourier series. The error gets smaller with increasing number of terms in the truncated series, and therefore d − also gets smaller. Without damping or nonlinearity, d − entirely vanishes, as was the case in Szalai (2014). When both d − and ε vanish, we arrive at a system that is subject to Utkin's closure in Sect. 3.2. If d − vanishes, but we have d + > 0 then for ε > 0 the trajectories are still unique, but there is no critical manifold in Σ that is being perturbed.

Remark 4.29
A more rigorous analysis would inspect the dynamics in the perturbed vector bundle corresponding to the near zero root of (4.47) for ε = 1. If this dynamics has a Lyapunov exponent σ 0 such that σ s < − |σ 0 | as in Assumption (A4), then this perturbed vector bundle could be attached to M λ , which would become a normally hyperbolic invariant manifold of the corrected model (4.28) in Σ. Fig. 6 Schematic of the nonlinear bowed string model. The continuous line represents the deformation of the string under vibration, and the dashed line represents the equilibrium shape of the string. μλ represents the friction force between the bow and the string, which acts at the contact point ξ = ξ

A Bowed Nonlinear String Model Reduced to Single Degree of Freedom
In this section, we illustrate the theory through a non-trivial example. In this example, the invariant manifold M λ is a linear subspace about an equilibrium that depends nonlinearly on the switching parameter λ. The dynamics within the invariant manifold given by f ( y, λ) and the switching function h 0 ( y, λ) are also nonlinear, which yields neither a Filippov-nor an Utkin-type model, but the skeleton model described in Sect. 4.2. In addition to the nonlinearity, we also include damping to make the invariant manifold attracting. We consider a nonlinear string with both ends rigidly held as illustrated in Fig. 6. The string has no resistance to bending, and any motion that occurs is due to the tension within the string. Whenever lateral deformation occurs, the string becomes stretched, which in turn causes an increase in tension and makes the model nonlinear. The tension is uniform along the length of the string. We denote the lateral deformation of the string by u(ξ, t), where ξ ∈ [0, 1] represents the distance along the string and t represents time. Moreover, we assume that this deformation occurs within a fixed plane so that u is a scalar-valued function. We also ignore any gravitational effect. We use primes to denote differentiation with respect to ξ and dots to denote differentiation in time. The dimensionless equation under our simplifying assumptions is where T is the tension within the string and Γ controls the nonlinearity of the string. The boundary conditions are u(0, t) = u(1, t) = 0 and where μ is a friction coefficient, ξ is the position of the contact point with the bow and u (ξ −, t), u (ξ +, t) represent the left and right derivative of u with respect to ξ at ξ , respectively. The boundary condition (5.2) reflects the equilibrium of forces at the contact point. The slope of the string together with the tension forms a force vector on both sides of the contact point. Since the string at the contact point is not smooth, the two force vectors do not cancel, and therefore to reach equilibrium an external force is necessary, supplied by the friction force μλ. The switching parameter λ decides the direction of the friction force and therefore changes sign as the relative velocity h = v 0 −u(ξ , t) between the bow and the string reverses, that is, To further simplify Eq. (5.1), we use second-order Taylor expansion, that is, √ 1 + u 2 ≈ 1 + 1 2 u 2 , which gives us the equation with boundary conditions We require that u(·, t) ∈ Lip ([0, 1], R), i.e., u(·, t) is Lipschitz continuous, which allows a finite contact force on the string. We define the operator The square root of D 2 , can be represented on the series u = a k sin kπξ by Du = kπa k sin kπξ . Note that D is not producing the first-order derivative. To represent all boundary conditions, we define a restricted D 2 as We also introduce damping with a constant damping ratio β ∈ [0, 1) for all vibration modes which transforms Eq. (5.3) intö Let us define x 1 = u(·), x 2 =u(·) and x = (x 1 , x 2 ) T , and hence we can write the system (5.3) as the infinite-dimensional dynamical systeṁ and the switching function is In order to represent solutions that were encountered in Sect. 2, we chose for the phase space of (5.5), where L ∞ stands for the space of bounded functions. The domain of definition is In what follows, we carry out a number of steps to arrive at the reduced-order model. These steps are applicable to systems where the invariant manifold is a spectral submanifold of an equilibrium. The steps are 1. Calculate the equilibrium of (5.5) as a function of λ, which is denoted by x . 2. Find the smoothest two-dimensional spectral submanifold (Haller and Ponsioen 2016) M λ about x , corresponding to the pair of complex conjugate eigenvalues with the least negative real part. The immersion of the manifold is denoted by where W fix is just one parametrization of M λ . y is an unknown and will be calculated in step 4. 3. Introduce an artificial parameter ε that slows down the dynamics on M λ to standstill at ε = 0 and has no effect at ε = 1. Then for ε = 0, calculate the invariant normal bundle of M λ , which is formed by the subspace orthogonal to the kernel of the adjoint A 0 ( y, λ) at each point on M λ . 4. Choose a coordinate shift y so that D 2 W falls into the invariant normal bundle of M λ at ε = 0, that is, W satisfies assumption (A4). This now fully specifies the immersion W . 5. Obtain the skeleton model by substituting the immersion W into (5.5). 6. Calculate the normal discontinuity gap from the dynamics in the invariant normal bundle of M λ . Also determine σ , the rate of convergence of the trajectory in the normal bundle with initial condition D 2 W . y 2 + y 2 (y 1 , λ) − c 2 (y 1 , λ) π 2 y 1 − 2γ (λ) sin πξ +2βπ y 2 + y 2 (y 1 , λ) − D 1 y 2 (y 1 , λ) y 2 + y 2 (y 1 , λ) with switching function h ε (y 1 , y 2 , κ, λ) = v 0 − y 2 + y 2 (y 1 , λ) sin πξ − εκ.
The normal discontinuity gap is λ) π c 2 (y 1 , λ) − β 2 and d + (y 1 , λ) The coordinate shift on the manifold in the velocity coordinate is any function that satisfies the differential equation The instantaneous square of the wave speed at the contact point is Proof These results are proven in Lemmas 5.2-5.4, 5.6 and 5.7.

The Invariant Manifold and Its Parametrization
We identify the invariant manifold M λ with the spectral submanifold (Haller and Ponsioen 2016) of the string's equilibrium corresponding to its first natural frequency. When λ is constant, the string has an equilibrium. We choose the smoothest invariant manifold about the equilibrium corresponding to the first natural frequency of the string, which is a unique two-dimensional linear subspace. We note that the theory of Cabré et al. (2003) does not apply, because damping makes backward-time solutions non-unique.  y 1 ( y, λ)) sin πξ (y 2 + y 2 ( y, λ)) sin πξ , (5.10) where γ (λ) is given by Eq. (5.8). The coordinate shift y = y 1 , y 2 T is not yet known.

Lemma 5.2 The immersion of the invariant manifold M λ about the equilibrium, as specified in steps 1 and 2 of the model reduction process, is
Proof We choose the representation of the invariant manifold as A substitution of W into the invariance Eq. (4.3) shows that W is indeed an immersion of an invariant manifold and corresponds to the first natural frequency. Because W is linear in y, M λ is also the smoothest invariant manifold.
The equilibrium x (λ) is calculated by setting the time derivative to zero in Eq. (5.3), which yields (5.12) In Eq. (5.12), Γ 2 1 0 x 2 1 dξ is independent of ξ , and therefore integrating (5.12) twice and applying the boundary conditions, we get (5.13) which still needs to be solved for x 1 . We define (5.14) which yields Physically, γ (λ)ξ (1 − ξ ) is the displacement of the string at the contact point at the equilibrium. To find the equation for γ , we substitute (5.15) into (5.14). We then evaluate the integral in (5.14), that is, so that Eq. (5.14) becomes (5.16) Equation (5.16) can be solved for γ with a single real solution, which is given by Eq. (5.8) that makes the equilibrium fully specified. Figure 7 shows the values of γ for various levels of nonlinearity. The nonlinearity is hardening, because the string deforms less than it would under the same force with the linear model. Substituting the equilibrium into (5.11) yields Eq. (5.10).

Linearized Dynamics About the Invariant Manifold
The linearized dynamics about M λ is characterized by the Frechet derivative A 1 ( y, λ) of Eq. (5.5), which is calculated here. (5.18) and the instantaneous square of the wave speed on M λ at the contact point is

Lemma 5.3 The Frechet derivative of F evaluated on M λ is
and D 2 x 1 = π 2 y 1 + y 1 ( y, λ) sin πξ. (5.20) The domain of definition of A 1 ( y, λ) is The only term in Eq. (5.5) not already linear is which we now linearize about a general point (x 1 , x 2 ) ∈ M λ . The expression (5.23) is a product, and hence we use the product rule when differentiating it with respect to x 1 . First, we linearize (5.23) about x 1 and get where we have used that z 1 must vanish at the boundaries ξ = 0, 1. Therefore, the first-order Taylor expansion of (5.23) is The value x 1 is the first component of the immersion of M λ , and when applying D 2 , we get (5.20). The remaining term in the Taylor approximation (5.24) is We also define the square of the instantaneous wave speed on M λ and at the contact point as The domain of definition of A 1 ( y, λ) must now include that x 1 (ξ −) − x 1 (ξ +) = 0, because the equilibrium is already included in the definition of W . However, x 1 (ξ −) − x 1 (ξ +) = 0 is already satisfied if x 1 ∈ L ∞ ([0, 1], R), and therefore we arrive at (5.21). To determine the closure D, we first note that x 1 must be Lipschitz continuous, and therefore x 1 is Lipschitz continuously differentiable. The closure of such functions in the Lipschitz norm is the continuously differentiable functions C 1 ([0, 1], R). Since D is the square root of D 2 , they have the same domain of definition, and therefore x 2 is Lipschitz continuously differentiable in L ∞ . The closure for this set in the L ∞ norm is the set of continuous functions C 0 ([0, 1], R). Summing up this argument, we have found (5.22). Lipschitz functions are not dense in C 1 , and continuous functions are not dense in the set of bounded functions either, and therefore Z = X.

Invariant Normal Bundle and Time-Scale Separation
There is no small parameter in Eq. (5.5) that controls the spectral gap between the tangential and normal dynamics about M λ . We therefore introduce such a scaling by artificially constructing A ε ( y, λ) such that for ε = 1, we recover the original dynamics and for ε = 0 the tangential dynamics becomesẏ = 0 when time is rescaled. This allows us to calculate the invariant normal bundle of M λ at ε = 0 and determine the parametrization M λ (i.e., the unknown coordinate shift y ( y, λ) in the immersion of M λ ) such that D 2 W ( y, λ) is strictly in the invariant normal bundle of M λ .
Proof Even though the mode shapes of the nonlinear string are the orthogonal harmonic functions sin kπξ , the contact force λ at ξ makes these modes intricately coupled. This coupling is represented by the term γ (λ)z 1 (ξ ) in Eq. (5.18). Nevertheless, we project A 1 ( y, λ) into two subspaces using the projection operators P : X → T y M λ , P x = 2 sin πξ 1 0 x 1 (Ξ ) sin πΞ dΞ 2 sin πξ 1 0 x 2 (Ξ ) sin πΞ dΞ (5.28) and Q = I − P. We calculate the projected operators where we can show that B 21 = 0. Introducing time-scale separation is just a multiplication of matrix B 1 by ε, which represents the rescaled linearized dynamics in the (invariant) tangent bundle of M λ . In the new coordinate system, the scaled linear operator becomes Given the form of A ε , the bundle projections assume the form where C is an unknown operator. Expanding the bundle invariance Eq. (4.9) yields which must hold for all v, Pv = 0. Further expanding (5.30), we get an equation for C in the form of The solution is C = B 12 B −1 2 at the critical parameter value ε = 0. We can now introduce another coordinate system in which x = ∞ k=1 a k sin kπξ, = 1, 2.

Remark 5.5
For Γ = 0, we have C = 0 and also D 2 y 2 = 0. This implies that for the linear string, the bundle projection is simply Q. The normal bundle is independent of ε, and there is no need to introduce time-scale separation. Instead of ε, Γ can be used to track the deformation of the invariant normal bundle, which persists for a sufficiently small Γ > 0 due to the properties of M λ (Bates et al. 1998).

Proof
The dynamics on the invariant manifold M λ is given by the invariance condition This is an equation in the tangent bundle of M λ , and therefore it makes sense to project it using P, as defined by (5.28), to find f . We first calculate that By inverting the matrix (5.35), we find that the reduced vector field is f ( y, λ) = 1 0 −D 1 y 2 (y 1 , λ) 1 P F (W ( y, λ), λ). (5.36) Next, we substitute the immersion (5.11) so that the vector field (5.5) on the manifold becomes

The Normal Discontinuity Gap d ± and Decay Rate
The normal discontinuity gap d ± measures the discontinuity of the correction about the invariant manifold with initial conditions D 2 W ( y, λ) at t = 0 and determines the uniqueness of solutions according to Theorem 4.13. We calculate d ± for the ε → 0 limit, when the dynamics in the normal bundle of M λ becomes autonomous. Therefore, it is sufficient to evaluate the limit lim t↓0 Dh · e A 0 ( y,λ)t D 2 W ( y, λ).

Lemma 5.7
The normal discontinuity gap in the limit ε → 0 is The rate of decay as defined by (4.10) is The calculation is carried out using Fourier series, and hence we write the series expansion sin kπξ k 2 π 2 sin kπξ D 2 y 2 (y 1 , λ) sin πξ , which is calculated from (5.31) by substituting (5.25). Since D 2 W ( y, λ) is in the invariant normal bundle of the critical manifold, it is sufficient to restrict the dynamics there. We use the decomposition of e A 0 ( y,λ)t as given by (5.29) to arrive at the expression sin kπξ k 2 π 2 sin kπξ 0 .

The relevant component of the solution is
sin kπξ k 2 π 2 sin kπξ. (5.40) The limit d + = − lim t↓0 z 2 (t)| ξ =ξ is calculated as Dh · z(t) = γ (λ) c 2 cos −1 β /c π c 2 − β 2 − D 2 y 2 (y 1 , λ) sin πξ (5.41) , and therefore we have shown (5.38). The calculation of (5.41) involves lengthy algebraic manipulations, converting the product of exponentials and trigonometric functions in (5.40) into sums of pure exponential expressions, which yields a sum of series with exponential terms. Each of the sub-series converge to logarithms of exponential polynomials. It then turns out that the result has discontinuities due to branch cuts of the complex logarithm and the limit at the branch cut brings the result. The detailed calculation (with slightly different notations) can be found in Sect. 2 of the electronic supplementary material of Szalai (2014). The decay rate (5.39) is found by reading off the smallest exponent from formula (5.40).

Remark 5.8
The normal discontinuity gap d ± is a local property of the string; it depends on material properties and the tension in the string. However, d ± is independent of the boundary conditions and the position where the string is forced.

Spectrum of the Normal Dynamics onẆ
e use Theorem 4.19 to find out whether there exists an attracting critical manifold.

Proof
The characteristic function (4.47), whose roots define stability, is formally written as It is possible to find a convergent series expansion of (s) by using Fourier series. Let us, for the moment, define We separate the solution into the first Fourier coefficient and the rest, such that where x 11 and x 12 are the coefficients of sin πξ and x 2 ,k are the coefficients of sin kπξ in the Fourier expansion of x. Now expanding Eq. (5.43) and using the notation (5.44) gives The solution to Eq. (5.46) for the k ≥ 2 Fourier coefficients is x 22,k = −2γ (λ) c 2 (y 1 , λ) sin kπξ s 2 + 2sβπk + c 2 (y 1 , λ)π 2 k 2 and For the first Fourier coefficients, the solution of Eq. (5.45) is sx 11 = 0 and 1 k 2 π 2 − c 2 (y 1 , λ) s 2 + 2sβπk + c 2 (y 1 , λ)π 2 k 2 sin 2 kπξ .
The series expansion of the characteristic function (s) using notation (5.44) is x 22,k sin kπξ + D 2 y 2 (y 1 , λ) sin πξ and substituting system parameters yields (5.42). Figure 8 shows the roots of (5.42) at an attracting point of the critical manifold. It can be seen that there is a real root near zero, while other roots are well inside the left complex half space. It seems that roots of (s) are perturbations of the eigenvalues of A 0 ( y, λ) apart from the rightmost root, that appears due to switching.
The plot of this rightmost root is shown in Fig. 9 in orange (solid lines), which indicates that the critical manifold is partly attracting (negative root), partly repelling (positive root). This might be surprising because the system dissipates energy as a whole. However, the constraint h = 0 and nonlinearity couple the dynamics in the tangent and normal bundles and energy is exchanged between them causing instability. In contrast, there is no such coupling in the linear string (Γ = 0), the green root near the origin in Fig. 8 remains at the origin, and therefore the normal dynamics is neutrally stable. Fig. 9 The rightmost root of the characteristic function (5.42) (orange, solid lines) and 3 × (− σ d − /d ± ) in (4.60) (green, dashed lines) determine whether the critical manifold is attracting at any given value of y and λ. Parameters are Γ = 10, μ = 1, β = 0.1 and ξ = √ 2/2 (Color figure online) Remark 5.10 Continuing from Remark 5.5, we find that for Γ = 0 and for all ε ∈ [0, 1] the characteristic function (5.42) is valid due to A 1 being constant. Let us denote the zero root of the characteristic function (5.42) for Γ = 0 by σ 0 . The invariant vector bundle corresponding to σ 0 is then isomorphic to M λ × R. For Γ > 0, the invariant vector bundle of σ 0 is continuously perturbed. The perturbation turns σ 0 into a small Lyapunov exponent of the now non-autonomous dynamics within the invariant vector bundle. Γ > 0 can be chosen such that σ s < −ε σ 0 , that is, the invariant vector bundle is an attracting normally hyperbolic invariant manifold in (M λ × X) ∩ Σ, that is the phase space of the corrected model in Σ. The dynamics in the invariant vector bundle is represented by the reduced-order model on M λ × R.

Equivalent Reduced-Order Model on 6 "
We now investigate the dynamics of the string on Σ ε . The dynamics outside Σ ε is given byẏ = f ( y, λ), εκ = σ κ andλ = 0. The skeleton model on Σ ε is formally given by Eq. (4.20), while the reduced-order model including a qualitative approximation of the normal dynamics is given by (4.59). The complication with Eqs. (4.20) and (4.59) is that they involve the lengthy term y 2 (y 1 , λ) as shown by Eq. (5.27). It is possible to eliminate y 2 (y 1 , λ) using the transformation y = y 1 , y 2 + y 2 (y 1 , λ) T .

Dynamics of the Skeleton Model on 6 0
In this section, we explore the dynamics of the skeleton model, which is the same as the dynamics on the critical manifold, when ε = 0 in Eq. (4.59). The dynamics on the critical manifold can be found by setting y 2 = v 0 / sin πξ andẏ 2 = 0 so that h = 0 andḣ = 0, then solvingẏ = f ( y, λ,λ) as an algebraic equation with (5.48) on the right side forλ, that is, As we noted in Theorem 4.8 in Sect. 4.2, solutions of this model pass through the boundaries Σ ± if D 2 y 2 (y 1 , λ) sin πξ = −d − > 0. To avoid any problem with having the wrong sign of d − , we rescale time by D 2 y 2 (y 1 , λ) and geṫ y 1 = D 2 y 2 (y 1 , λ)v 0 / sin πξ λ = c 2 (y 1 , λ) π 2 y 1 − 2γ (λ) sin πξ + 2βπv 0 / sin πξ , (5.53) whose forward-time solutions always pass through Σ ± . This allows for a straightforward numerical solution. Let us first recall, what Utkin's closure would produce if we disregarded D 2 y 2 (y 1 , λ). The solution would be given by the equation which is partly algebraic. Figure 10a shows the phase portrait. The dashed orange lines correspond to λ values jumping between either ± 1 or the solution of the algebraic constraint in Eq. (5.54). The continuous green line represents λ values that are admissible by the constraint in Eq. (5.54). Sinceẏ 1 is a positive constant, solutions can only move in one direction along the green line. This is typical of friction oscillators, and it is consistent with rigid body dynamics, where forces are allowed to be discontinuous. A different picture emerges when D 2 y 2 (y 1 , λ) is considered. Figure 10b shows a typical phase portrait of (5.52). Parts of trajectories are denoted by dashed lines when D 2 y 2 (y 1 , λ) < 0. The arrows indicate the correct forward direction of time. Black lines indicate when D 2 y 2 (y 1 , λ) = 0, and therefore Eq. (5.52) is singular and the direction of time changes in Eq. (5.53). The black lines also form a set of nullclines of Eq. (5.53), because at these pointsẏ 1 = 0. Another nullcline is shown in green, wherė λ = 0. At the intersection of the green and black lines, Eq. (5.53) has an equilibrium, which is a node. The weak stable manifold of this equilibrium is close to the green nullcline ofλ = 0.
The phase portrait of Fig. 10b is not typical for a friction oscillator. Yet, the skeleton model is obtained through a careful reduction to an invariant manifold, where we made sure any perturbation due to the discontinuities would only affect the invariant normal bundle. Applying Filippov's closure at the boundaries λ = ± 1 yields sliding solutions. However, this implies that solutions coming from either side of Σ cannot enter Σ, while λ stays at ± 1. Having a fixed value of λ is not physical in a friction oscillator. The singularities within Σ are reached infinitely fast. This resembles the dynamics of the van der Pol oscillator at the fold point of its critical manifold (Kanamaru 2017) or in general the dynamics of singularly perturbed systems (Kuehn 2015). For example, the equilibrium of Eq. (5.53), formed by the intersection of two nullclines, resembles folded-node singularities (Wechselberger 2012;Kristiansen 2017). Therefore, our best chance to gain more insight is to consider the reduced-order model (5.52) that includes a representation of the dynamics in the invariant normal bundle of M λ as described in Sect. 4.5.

Dynamics of the Reduced-Order Model on "
In this section, we investigate the reduced-order model (4.57), which is the extension of the skeleton model by a single variable representing the dynamics in the normal bundle of M λ . The dynamics on Σ ε is given by Eq. (5.51) with parameters derived in Sect. 5.5. Proposition 4.27 shows that the reduced-order model captures the stability of M λ for ε = 0 well. Figure 9 confirms this: In the illustrated part of the phase space, the stability of the critical manifold of the corrected model and the reduced-order model is the same. The critical manifold is repelling where d − > 0 as per Proposition 4.25. The skeleton model does not capture the repelled trajectories and also displays singular dynamics for d − > 0 as shown in Fig. 10b. Here, we illustrate that the positive value of d ± for the reduced-order model resolves the singularities that occur in the skeleton model according to Proposition 4.24.
We first choose a small parameter value ε = 10 −7 to show the qualitative differences between Eqs. (5.52) and (5.51). Figure 11a shows the two-dimensional projection of the phase portrait ignoring variable κ. When trajectories start in the shaded part with λ = 1, where the critical manifold is repelling, they quickly pass to λ = − 1 without much change in y 1 , while κ exponentially explodes. Trajectories starting with λ = − 1, in the region where the critical manifold is attracting, follow the manifold while being attracted to the stable node of (5.52) at the intersection of the green and black lines. At the node, the stability of the critical manifold changes and trajectories are again repelled with growing magnitude of κ. This is illustrated in Fig. 11b. After passing the node, trajectories tend to either λ = ± 1. It is then likely that trajectories will start a violent oscillation between λ = ± 1, because they interact with the two repelling parts of the critical manifold. This dynamics has some resemblance to Fig. 10b except that there is no need to rescale time, since there is no division by d − .
Increasing ε leads to less violent oscillations between λ = ± 1, which eventually continues without reaching λ = ± 1. Such a case is shown in Fig. 11d, where the oscillation is reduced to a single loop about the line where the critical manifold becomes repelling. For ε = 1, the dynamics becomes relatively slow for all variables and resembles that of typical friction oscillators with well-defined stick and slip phases. This phase portrait is shown in Fig. 11c. For ε sufficiently large the time scale of the normal dynamics (κ variable) becomes much longer than the dynamics of the rest of the variables, and therefore during a stick phase κ does not change much, which also means that the instability of the critical manifold loses its influence on the dynamics. Indeed, the leading characteristic root of (s) is a small perturbation of the zero root, and hence it is easily dominated by other time scales. In fact by removing nonlinearity (Γ = 0), this root remains zero, and hence κ simply becomes an integral of other quantities without a dynamics of it own. In our example at ε = 1, κ is almost without its own dynamics. The justification why ε can be increased to ε = 1 can be found in Remark 5.10.
The conclusion from the analysis is that simply applying reduction to an invariant manifold is not sufficient, one needs to take into account at least a qualitative approximation of the normal dynamics. This is because the skeleton model (4.13) over-emphasizes instabilities and turns them into singularities. The main component that makes the reduced-order model (4.57) well behaved is that d ± is positive in all parts of the phase space. For the nonlinear string example, d ± is the velocity jump of the contact point due to a unit jump in λ, i.e., the contact force. Therefore in light of Newton's second law, it is understandable why d ± > 0. In case we had found d ± = 0, the reduced-order model (4.57), including an approximation of the normal dynamics about M λ , would not be necessary, the skeleton model would be sufficient.

Conclusion
In this paper, we have investigated PWS systems on Banach spaces with non-dense domain of definition. Specific application areas that satisfy this assumption are the elastodynamics equations (Kausel 2006), delay equations (Diekmann et al. 1995) or age-dependent population dynamics (Metz and Diekmann 1986). Such systems are different from other classes of PWS systems, because they can have unique solutions under general conditions. We were also able to construct a finite-dimensional reducedorder model that inherits key properties of an infinite-dimensional model. Non-dense domain of definition can arise if the phase space is non-reflexive, for example, when the phase space consists of continuous, bounded or Lipschitz continuous functions. In some cases, boundary conditions can make the domain non-dense (Neubrander 1988).
The key quantity that decides uniqueness of solutions is the normal discontinuity gap, which is due to discontinuous trajectories that systems with non-dense domains have. For the linear and nonlinear string, the normal discontinuity gap represents the velocity jump of a contact point in response to a unit jump in force. The presence of the normal discontinuity gap allows the dynamics inside the switching manifold to become smooth. As a result, two new discontinuity boundaries arise, where trajectories can enter or leave the switching manifold. If the normal discontinuity gap is positive, trajectories cross the new discontinuity boundaries under general conditions.
Despite uniqueness of solutions, invariant manifolds that extend over the switching manifold do not exist. We have assumed the existence of an invariant manifold when the switching parameter of the vector field is constant. This invariant manifold does not persist when the switching parameter varies, but we have found that pieces of this manifold do persist, while discontinuities between the persisting pieces develop along the two new discontinuity boundaries. We have also shown that switching can make the invariant manifold repelling. However, in the example of the nonlinear string the invariant manifold is repelling only in a single direction, which can be captured by a scalar variable. We have constructed a reduced-order model that captures this instability. It remains to be shown under what conditions there is a spectral gap between the reduced model and the rest of the dynamics, so that the reduced-order model captures all the essential dynamics. We have only shown that the invariant manifold becomes repelling within the reduced-order model and within the infinite-dimensional system under the same conditions through a real root (see Proposition 4.27).
While the theory presented is incomplete, we hope that the results in this paper will find applications in simulating PWS continuum systems. Using the reduced-order model instead of the skeleton model eliminates singularities and allows for a unique solution. This allows well-conditioned numerical schemes that lead to robust solutions unlike what is currently possible (Kane et al. 1999). While it is not proven that the reduced-order model fully captures all dynamics, we expect that this will be shown in the future either in general or under further conditions.
We have demonstrated the model reduction procedure on a bowed nonlinear string model. In this example, we have found that the skeleton model has non-physical singularities, where the friction force between the bow and the string remains at its maximal limit. The skeleton model also has a singularity that resembles a folded node of singularly perturbed systems (Wechselberger 2005;Kristiansen 2017). After correcting the skeleton model with the dynamics that arises in the eliminated parts of the system due to switching, the pictures becomes clearer. It turns out that the correction is a largely decaying motion with the possibility of an instability along a one-dimensional subspace. When this possible instability is taken into account, the model becomes free of singularities and the dynamics resembles what a friction oscillator would exhibit when the friction force is regularized (Sotomayor and Teixeira 1998).
Note that the harmonic term in (A.9) is canceled by the homogeneous solution of (A.11), and hence .12) and the initial condition of (A.11) vanishes κ = π 2 y λ − κ , κ(0) = 0,κ(0) = 0. (A.13) The switching function (A.10) with (A.12) and (A.13) takes into account the full perturbation (A.1) exactly. If we are seeking to solve for a finite time interval, infinitely long delays in (A.12) can be neglected. In case of very short simulation on the interval 0 ≤ t < min (2ξ , 2 − 2ξ ), it is sufficient to use which then yields which is the result we sought.

B Proof of Theorem 4.8
The following proof of Theorem 4.8 investigates whether a trajectory approaching Σ ± 0 can be continued uniquely after reaching Σ ± 0 in the two cases set out by the theorem. Proof Both of the Eqs. (4.16) and (4.20) that govern the dynamics on the two sides of Σ ± 0 already have unique solutions. We need to exclude the possibility that a trajectory can be continued by both Eqs. (4.16) and (4.20) simultaneously and also exclude the existence of a sliding trajectory on Σ ± 0 . We denote the solution of (4.16) by (η(t), λ ), and the solution of (4.20) by (σ (t), λ(t)) either of which can form T .
(B.9) So far, we have shown that if vector field (4.20) is tangent of order − 1 to Σ ± 0 at ( y , λ ) then so is (4.16) and the orientation of the tangencies are the same. We now show that this sufficient condition is also necessary.
Using Eq. (B.9) for = 1 does not require assumption (B.3). It directly follows from Eq. (B.9) that Now knowing that d dt λ(t) = 0 we can apply (B.9) for = 2 and conclude that d 2 dt 2 h 0 (η(t), λ )| t=0 = 0 ⇒ Using assumption (4.21) and Eq. (B.10), we conclude that the order and orientation of the tangency are the same on both sides of Σ ± 0 . If is odd, trajectory T passes through Σ ± 0 at the tangency. If is even, the trajectory continues on the same side of Σ ± 0 and there is no joining trajectory from the other side of Σ ± 0 . Conditions (4.21) and (4.22) also imply that either case 1 or 2 holds for points on Σ ± 0 in a sufficiently small open neighborhood of ( y , λ ). This excludes cases where trajectories are forced onto Σ ± 0 for a nonzero length of time and implies that there cannot be trajectories joining ( y , λ ) from within Σ ± 0 . Therefore, there is a unique continuation of T for t > 0 (or t < 0) sufficiently small.

C Proofs of Lemma 4.10 and Theorem 4.13
There are three steps to the proof of Theorem 4.13. First, we derive a differential equation for λ from the algebraic constraint h(W ( y, λ) + z) = 0 by differentiation, which follows from Lemma 4.10. At this step, we claim continuity ofλ. By investigating the resulting equation, we move to step two and establish thaṫ λ is indeed continuous in Σ. In the final step, we show that trajectories transverse to Σ ± cross Σ ± when d ± ( y, z, λ) > 0, which concludes the proof.
For convenience, we copy here Lemma 4.10.
The scheme of the finite difference is such that for k = m − 1 the second and first argument of η are equal in one of the terms, i.e., t − δ + c k = s + kδ + c k , which takes into account the discontinuity of η. If this discontinuity is not taken into account, the integral in the limit t ↓ s would vanish. In summary, we have the integral (η(t + c k , s + kδ + c k ) − η(t − δ + c k , s + kδ + c k ))λ(s + kδ + c k ).
We now prove Theorem 4.13.

(C.9)
This concludes the first step of the proof. We now demonstrate thatλ is indeed continuous. We only need to recall that the formal solution for any history of λ is which allows us to write that Dh(x) · A 1 ( y, λ)z(t) = d dt Dh(x) · U(t, s)z(s) − t s d dt Dh(x) · K (t, τ )λ(τ )dτ.
(C.10) Under the assumptions of Theorem 4.13, the expression (C.10) is continuous in t and so is (C.9), which concludes the second part of the proof.