Abstract
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.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 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 time-scale 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. 2013, 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, 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.
2 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
where \(u(\xi ,t)\) is the scalar-valued displacement of the string, \(t\in \mathbb {R}\) is time and \(\xi \in [0,1]\) is the spatial coordinate along the string; \(\lambda \in [-\,1,1]\) is the force applied at the contact point \(\xi ^{\star }\), prime denotes differentiation with respect to \(\xi \) and dot with respect to t; \(u'(\xi ^{\star }-,t)\), \(u'(\xi ^{\star }+,t)\) denote the left and the right derivative at \(\xi ^{\star }\), respectively. The value of parameter \(\lambda \) is given by
where \(v_{0}\) is the speed of the bow and \(\lambda \) 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 \(\lambda \) is not defined for \(h=0\) by Eq. (2.2). We assume that all values of \(\lambda \in [-\,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 \(\lambda \) using the condition that the functions \(u(\xi ^{\star },\cdot )\) and \(\dot{u}(\xi ^{\star },\cdot )\), that is, the solution of (2.1) and (2.2) evaluated at \(\xi =\xi ^{\star }\) must be continuous in time.
We now consider the case when \(\lambda \) is constant. For constant \(\lambda \) Eq. (2.1) has an equilibrium. The equilibrium shape of the string for \(\lambda =1\) is given by
where H is the Heaviside function. Due to linearity, for a fixed \(\lambda \) the equilibrium is then \(\lambda u_{0}(\xi )\). It is known that free vibrations of a string, that is, the solutions of (2.1) with constant \(\lambda \) can be written as
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 \(\mathcal {M}\). For this set of solutions, we denote the displacement of the string at \(\xi ^{\star }\) by \(y(t)=\lambda y^{\star }+\left( a_{1}\sin \pi t+b_{1}\cos \pi t\right) \sin \pi \xi ^{\star }\), where \(y^{\star }\) is yet to be defined. The value y(t) can be used to recover the displacement of the whole string as
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^{\star }\) gives the same manifold, \(y^{\star }\) only influences the parametrization of the manifold. To fix \(y^{\star }\), we require that the manifold does not move in the tangential direction of \(\mathcal {M}\) when \(\lambda \) changes. This means that the two partial derivatives \(\frac{\partial }{\partial \lambda }u_{\mathcal {M}}(y,\lambda ;\xi )\) and \(\frac{\partial }{\partial y}u_{\mathcal {M}}(y,\lambda ;\xi )\) must be perpendicular, that is,
Solving Eq. (2.6) for \(y^{\star }\), we get
Substituting the immersion (2.5) of manifold \(\mathcal {M}\) into (2.1) and (2.2) while assuming that \(\lambda \) is constant, we get
where \(h=v_{0}-\dot{y}\). 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 \(\lambda \) 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\), \(\lambda \) 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 \(\dot{y}=v_{0}\) on this interval, and consequently, \(\ddot{y}=0\). Substituting the conclusion of this argument into the first line of (2.7) yields
We only allow \(\lambda \in [-\,1,1]\), and therefore if the result of (2.8) is outside of the interval \([-\,1,1]\), \(\lambda \) 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 \(\lambda \). The solid green line on the horizontal red-shaded plane is the stick solution, where the friction force \(\lambda \) grows with a constant rate with respect to t and y until it reaches the limit \(\lambda =\pm \,1\) and slip ensues.
By assuming constant \(\lambda \), we made an error when the relative velocity between the bow and the string becomes zero, i.e., \(\dot{y}=v_{0}\), because contrary to the assumption \(\lambda \) is not constant and can even jump between \(-\,1\) and 1 or from \(\pm \,1\) to the value given 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
where \(w(\xi ,t)\) is a correction term. Assuming that y(t) satisfies (2.7) and substituting (2.9) into (2.1) without assuming constant \(\lambda \), we get the governing equation of the correction term
which describes how the dynamics depart from the invariant manifold when \(\lambda \) varies. Starting from the invariant manifold \(\mathcal {M}\), we have initial conditions
The choice of \(y^{\star }\) dictated by Eq. (2.6) guarantees that the correction term as described by (2.10) does not include vibrations with the first natural frequency of the string (In abstract terms, w is restricted to the invariant normal bundle of \(\mathcal {M}\).) Consequently, the skeleton model (2.7) does not require any correction even when \(\lambda \) varies. Instead, the switching function (2.2) needs to be revised by substituting the corrected solution (2.9). If we use (2.9) in Eq. (2.2), we can write the switching function as
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 \(\lambda \) leads to a jump in the velocity \(\dot{w}(\xi ^{\star },t)\). For this, we assume that \(\lambda (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 \(\dot{w}(\xi ^{\star },t)\) at \(t=0\) has a discontinuity, whose gap is proportional to the jump in \(\lambda \). 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 at \(t=0\). At \(t=0\), we have
which we call the normal discontinuity gap. Due to the linearity of Eq. (2.10), any jump in \(\lambda \) is translated into a discontinuity of the velocity \(\dot{w}(\xi ^{\star },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\le t<\min \left( 2\xi ^{\star },2-2\xi ^{\star }\right) \), before any discontinuity is reflected back to \(\xi =\xi ^{\star }\), we find that the switching function (2.12) in Eq. (2.7) becomes
where \(\kappa \) 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\le t<\min \left( 2\xi ^{\star },2-2\xi ^{\star }\right) \). 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 \(\lambda \), because of the presence of the normal discontinuity gap (2.13). Therefore, when \(h=0\), \(\lambda \) can be solved for, without any closure rule, such as Filippov’s or Utkin’s (see Sect. 3). In our case, solving the equation \(h=0\) for \(\lambda \) yields
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 get
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\) into \(\dot{h}=0\), and therefore (2.18) is valid for initial conditions that satisfy \(h=0\) with \(\lambda \in [-\,1,1]\).
We can now put together the whole model with the correction into a single system
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 \(\mathcal {M}\) and for the time interval \(0\le t<\min \left( 2\xi ^{\star },2-2\xi ^{\star }\right) \). Appendix A shows that the valid time interval can be extended to any length by including delayed values of \(\lambda \) 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 \(\lambda =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 \(\lambda \) 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 \(\lambda \). Due to the higher dimensional dynamics that arise from the inclusion of \(\kappa \) as a dynamic variable and delayed values of \(\lambda \), 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.
3 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_{k}^{j}\) denotes the jth derivative with respect to the kth argument.
Let us consider the vector field
where either
G is a compact and connected subset of \(\mathbb {R}^{n}\) and \(n,p\in \mathbb {N}^{+}\). The function \(h\in C^{p}(G,\mathbb {R})\) is called the switching function, and its zero set defines the switching manifold
A solution \(\varvec{y}:I\rightarrow \mathbb {R}^{n}\) of Eqs. (3.1) and (3.2) is defined on a closed interval of nonzero length \(I\subset \mathbb {R}\). There is no information in Eqs. (3.1) and (3.2) that helps to deduce a value for \(\lambda \) on \(\varSigma \). To explore all possibilities (3.1) can be turned into a differential inclusion
where \(\overline{\mathrm {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 \(\lambda \) on \(\varSigma \).
We have assumed two possibilities, (3.3) or (3.4), for the domain of definition of \(\varvec{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 \(\varSigma \) 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 \(\varvec{y}(0)\in G\) if \(h(\varvec{y}(0))\ne 0\), because \(\varvec{f}\) is a smooth vector field (Coddington and Levinson 1955). However, for \(h(\varvec{y}(0))=0\) the vector field is not defined and one needs to reason how trajectories continue once they reach \(\varSigma \).
The simplest case of a trajectory interacting with \(\varSigma \) is when the trajectory approaches \(\varSigma \) 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 \(\varSigma \) is denoted by \(\varvec{y}^{\star }\), the speeds at which h increases are \({Dh}(\varvec{y}^{\star })\varvec{f}(\varvec{y}^{\star },\pm \,1)\). h changes monotonically through \(h=0\) with a nonzero speed if and only if
In case of (3.7), there is no need to define the dynamics on \(\varSigma \) because the value of \(\lambda \) simply switches between \(\pm \,1\). In our argument, we have used \(\varvec{f}\) for \(\lambda =\pm \,1\) only, therefore, to resolve crossing, it is sufficient to assume (3.3). The subset of \(\varSigma \), where (3.7) holds, is called the crossing region and denoted by \(\varSigma _{cr}\).
Now we discuss the case when
When (3.8) holds, \(\varSigma \) is an attractor in either forward or backward time. This means that trajectories cannot immediately escape \(\varSigma \) once they are on \(\varSigma \). The subset of \(\varSigma \), where (3.8) holds and attracts solutions in forward time is called the sliding region and denoted by \(\varSigma _{sl}\). The repelling subset of \(\varSigma \), where (3.8) holds, is called the escaping region and denoted by \(\varSigma _{\mathrm{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 \(\varSigma \). 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).
3.1 Filippov’s Closure
Filippov’s closure defines a vector field on \(\varSigma \), when condition (3.8) holds, by interpolating between the vector fields \(\varvec{f}(\varvec{y},\pm \,1)\), such that \(\varSigma \) becomes an invariant manifold of the new vector field. The interpolation is carried out as follows. We define a new vector field
where
For \(\varvec{y}\notin \varSigma \) Eq. (3.2) still defines \(\lambda \) and Eq. (3.9) is identical to (3.1) for \(\lambda =\pm \,1\). On \(\varSigma \) and when (3.8) holds we calculate \(\lambda \) from
which stipulates that the vector field (3.9) is tangential to \(\varSigma \). The solution of (3.10) is
Definition 3.2
Assume that (3.8) holds. We call the vector field (3.9), where \(\lambda \) is given by (3.11), Fillipov’s closure.
It can be shown that \(\lambda \in (-\,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 \(\varvec{f}(\varvec{y},\pm \,1)\) so that \(\varvec{r}(\varvec{y})+\varvec{b}(\varvec{y})\lambda \) is tangential to \(\varSigma \). A trajectory at its first point of contact with \(\varSigma \) is continuous, but not continuously differentiable, because \(\lambda \) becomes discontinuous due to (3.11).
When neither (3.7) nor (3.8) holds for \(\varvec{y}^{\star }\in \varSigma \), we have
Equation (3.12) means that one or both of the vector fields \(\varvec{f}(\varvec{y},\pm \,1)\) is tangential to \(\varSigma \), 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 \(\varSigma \).
Trajectories may not have unique continuation when they are tangent to \(\varSigma \). When both vector fields \(\varvec{f}(\varvec{y},\pm \,1)\) are tangential to \(\varSigma \) at a point, \(\lambda \) 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.
3.2 Utkin’s Closure
In this section, we assume that the domain of definition of \(\varvec{f}\) is given by (3.4). In the case of (3.8), we similarly construct the vector field on \(\varSigma \), such that \(\varSigma \) becomes invariant under the vector field. The invariance of \(\varSigma \) is expressed as
Equation (3.13) has at least one solution for some \(\lambda \in [-\,1,1]\), because (3.4) implies that \({Dh}(\varvec{y})\cdot \varvec{f}(\varvec{y},\pm \,1)\) has different signs and due to Bolzano’s theorem there must be a root.
Definition 3.3
Assume that (3.8) holds. We call the vector field (3.9), where \(\lambda \) is given by the solutions of Eq. (3.13), Utkin’s closure.
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 \(\varvec{f}(\varvec{y},1)\) to \(\varvec{f}(\varvec{y},-\,1)\) represents the possible values of the vector field on \(\varSigma \). There is one intersection of this family of vectors with the tangent plane of \(\varSigma \), represented by the thick red arrow, which satisfies Eq. (3.13). Figure 4b shows that there can be multiple intersections of \(\varvec{f}(\varvec{y},[-\,1,1])\) with the tangent plane of \(\varSigma \), 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 \(\lambda \in (-\,1,1)\).
4 Model Reduction
We start with a general continuum model, in the form of
where \(\varvec{x}\) is a function of time \(t\in [s,\infty )\), that is \(\varvec{x}{:}\,[s,\infty )\rightarrow \varvec{X}\) with an initial condition \(\varvec{x}(s)=\varvec{x}_{0}\) and \(\varvec{X}\) is an appropriately defined Banach space. The domain of definition of \(\varvec{F}(\cdot ,\lambda )\), for a fixed \(\lambda \), is denoted by \(\varvec{\mathcal {D}}_{\lambda }(\varvec{F})\subset \varvec{X}\) so that the full domain of definition is \(\varvec{\mathcal {D}}(\varvec{F})=\{(\varvec{\mathcal {D}}_{\lambda }(\varvec{F}),\lambda ):\lambda \in [-\,1,1]\}\) and \(\varvec{F}:\varvec{\mathcal {D}}(\varvec{F})\rightarrow \varvec{X}\). The switching function h is defined on \(\overline{\cup _{\lambda }\varvec{\mathcal {D}}_{\lambda }(\varvec{F})}\) and has values in \(\mathbb {R}\). When \(h=0\), the most general definition of the dynamics is \(\dot{\varvec{x}}\in \overline{\mathrm {co}}\varvec{F}(\varvec{x},[-\,1,1])\). We also require that trajectories are continuous, even when \(h=0\). The smoothness of \(\varvec{F}\) and h is not assumed globally, instead we assume the smoothness of an invariant manifold of \(\varvec{F}\) and related quantities in the next section.
Remark 4.1
The notation of Eq. (4.1) facilitates that \(\lambda \) is an unknown, which needs to be found when \(h(\varvec{x})=0\). Therefore, \(\lambda \) may not be a function of \(\varvec{x}\), but it may become part of the phase space. The solution for \(\lambda \), when \(h(\varvec{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 \(\lambda \).
4.1 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 \(\lambda \) 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 \(\varvec{W}\in C^{p}(\mathbb {R}^{n}\times [-\,1,1],\varvec{X})\), \(p\ge 2\) and a vector field \(\varvec{f}\in C^{p}(G\times [-\,1,1],\mathbb {R}^{n})\), which satisfies the invariance condition
$$\begin{aligned} \varvec{F}(\varvec{W}(\varvec{y},\lambda ),\lambda )=D_{1}\varvec{W}(\varvec{y},\lambda )\varvec{f}(\varvec{y},\lambda ), \end{aligned}$$(4.3)where G is a compact and connected subset of \(\mathbb {R}^{n}\). The invariant manifold is given by
$$\begin{aligned} \mathcal {M}_{\lambda }=\left\{ \varvec{W}(\varvec{y},\lambda ):\varvec{y}\in G,\lambda \in [-\,1,1]\right\} \end{aligned}$$(4.4)and the dynamics of (4.2) on \(\mathcal {M}_{\lambda }\) is described by
$$\begin{aligned} \dot{\varvec{y}}=\varvec{f}(\varvec{y},\lambda ). \end{aligned}$$(4.5)\(\varvec{W}\) is called the immersion of \(\mathcal {M}_{\lambda }\).
- (A2):
-
We assume that for every \(\lambda \in [-\,1,1]\), \(\varvec{F}(\cdot ,\lambda )\) is Frechet differentiable on \(\mathcal {M}_{\lambda }\). This derivative is denoted by
$$\begin{aligned} \varvec{A}_{1}(\varvec{y},\lambda )=D_{1}\varvec{F}(\varvec{W}(\varvec{y},\lambda ),\lambda ). \end{aligned}$$We also assume that the domain of definition of \(\varvec{A}_{1}\), i.e., \(\varvec{\mathcal {D}}(\varvec{A}_{1}(\varvec{y},\lambda ))=\left\{ \varvec{x}\in \varvec{X}:\varvec{A}_{1}(\varvec{y},\lambda )\varvec{x}\in \varvec{X}\right\} \), is independent of \(\varvec{y}\) and \(\lambda \), and we define \(\varvec{Z}=\overline{\varvec{\mathcal {D}}(\varvec{A}_{1}(\varvec{y},\lambda ))}\) (in general, \(\varvec{\mathcal {D}}(\varvec{A}_{1}(\varvec{y},\lambda ))\ne \varvec{\mathcal {D}}_{\lambda }(\varvec{F})\)).
- (A3):
-
Unique continuous solutions. We assume that the abstract Cauchy problem
$$\begin{aligned} \left. \begin{array}{rl} \dot{\varvec{y}} &{} =\varvec{f}(\varvec{y},\lambda )\\ \dot{\varvec{z}} &{} =\varvec{A}_{1}(\varvec{y},\lambda )\varvec{z}-D_{2}\varvec{W}(\varvec{y},\lambda )\dot{\lambda } \end{array}\right\} \end{aligned}$$(4.6)with initial conditions \(\varvec{y}(s)\in G\), \(\varvec{z}(s)\in \varvec{Z}\), \(s\in \mathbb {R}\) and with \(\lambda \in C^{1}([s,\infty ),\mathbb {R})\) has a unique solution \(\left( \varvec{y},\varvec{z}\right) \in C([s,\infty ),G\times \varvec{Z})\), even though we only have \(D_{2}\varvec{W}(\varvec{y}(t),\lambda (t))\in \varvec{X}\). We also assume that the \(\varvec{Z}\) component of the solution can be written as
$$\begin{aligned} \varvec{z}(t)=\varvec{U}(t,s)\varvec{z}(s)-\int _{s}^{t}\varvec{K}(t,\tau )\dot{\lambda }(\tau )\mathrm {d}\tau , \end{aligned}$$(4.7)where \(\varvec{K}\) is bounded for \(\tau \ge s\) and continuous in both variables for \(\tau >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):
-
\(\mathcal {M}_{\lambda }\) is attracting and normally hyperbolic. We assume that there exist two families of projections \(\varPi ^{c}(\varvec{y},\lambda )\) and \(\varPi ^{s}(\varvec{y},\lambda )\), strongly continuous in \(\varvec{y}\) and \(\lambda \) such that
$$\begin{aligned} \varPi ^{c}(\varvec{y},\lambda )+\varPi ^{s}(\varvec{y},\lambda )&=\varvec{I},\nonumber \\ \varPi ^{c}(\varvec{y},\lambda )D_{1}\varvec{W}(\varvec{y},\lambda )&=D_{1}\varvec{W}(\varvec{y},\lambda ),\end{aligned}$$(4.8)$$\begin{aligned} \varvec{U}(t,s)\varPi ^{s}(\varvec{y}(s),\lambda (s))\varvec{z}&=\varPi ^{s}(\varvec{y}(t),\lambda (t))\varvec{U}(t,s)\varvec{z},\quad \forall \varvec{z}\in \varvec{Z},\quad t\ge s. \end{aligned}$$(4.9)Consider the non-autonomous ordinary differential equation \(\dot{\varvec{\eta }}=D_{1}\varvec{f}(\varvec{y},\lambda )\varvec{\eta }\), whose solutions with initial condition \(\varvec{\eta }_{0}\) at \(t=s\) are denoted by \(\varvec{\eta }(t,s,\varvec{\eta }_{0})\). We assume that there exist real numbers \(\sigma _c >0, \sigma _s < - \sigma _c\) ,\(M_{c}>0\) and \(M_{s}>0\) such that
$$\begin{aligned}&\forall (t-s)\in \mathbb {R},\,\varvec{\eta }_{0}\in \mathbb {R}^{n} :\left\| \varvec{\eta }(t,s,\varvec{\eta }_{0})\right\| \le M_{c}\left\| \varvec{\eta }_{0}\right\| \mathrm {e}^{\sigma _{c}\left| t-s\right| },\\&\forall s\le t,\,\varPi ^{s}(\varvec{y}(s),\lambda (s))\varvec{z}=\varvec{z} :\left\| \varvec{U}(t,s)\varvec{z}\right\| \le M_{s}\left\| \varvec{z}\right\| \mathrm {e}^{\sigma _{s}(t-s)}. \end{aligned}$$ - (A5):
-
We assume that for \(t\ge s\) there exists \(0<M<\infty \) and \(\sigma <0\) such that
$$\begin{aligned} \left\| \varvec{K}(t,s)\right\| \le M\mathrm {e}^{\sigma (t-s)}. \end{aligned}$$(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): \(\varvec{z} \left( \varvec{y},\varvec{z}\right) \in C([s,\infty ),G\times \varvec{Z})\) is an F-solution of (4.6) if there exists a sequence \(\varvec{z}_{k}\in C^{1}([s,\infty ),\varvec{Z})\cap C^{1}([s,\infty ),\varvec{Z})\) such that
where \(\left\| \varvec{z}\right\| _{\infty }=\sup _{t\in [s,\infty )}\left\| \varvec{z}(t)\right\| \) and \(\varvec{y}\) satisfies \(\dot{\varvec{y}}=\varvec{f}(\varvec{y},\lambda )\).
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 (1990, 2008), which implies that the convolution integral is
because \(\left( \mu -\varvec{A}_{1}(\varvec{y}(\tau ),\lambda (\tau ))\right) ^{-1}:\varvec{X}\rightarrow \varvec{\mathcal {D}}\). da Prato and Sinestrari (1992) uses a similar technique to approximate the unique solution. The kernel \(\varvec{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 \(\varvec{K}\) for \(t>s\), which allows for a discontinuity at \(t=s\) due to \(D_{2}\varvec{W}(\varvec{y}(\tau ),\lambda (\tau ))\notin \varvec{Z}\).
Remark 4.5
The uniqueness or persistence of \(\mathcal {M}_{\lambda }\) is not addressed by the assumptions. For persistence of \(\mathcal {M}_{\lambda }\) 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 \(\dot{\lambda }\) is bounded. This will be useful later when the reduced-order model is constructed.
Figure 5 shows the invariant manifold \(\mathcal {M}_{\lambda }\) and its intersection with the switching manifold \(\varSigma \). The Banach space \(\varvec{X}\) is represented by two coordinates \(\varvec{x}_{1}\) and \(\varvec{x}_{2}\). The two parts of the invariant manifold (\(\mathcal {M}_{-1}\) and \(\mathcal {M}_{1}\)) do not join up in Fig. 5a. If trajectories cross \(\varSigma \) instantaneously, they are discontinuous. When discontinuity is not allowed, the crossing cannot be instantaneous. In certain cases, however, the disconnected nature of \(\mathcal {M}_{\lambda }\) may be overlooked. For example, when the switching function solely depends on the parameter \(\varvec{y}\) of the immersion \(\varvec{W}\), i.e., \(\varvec{y}=\varvec{x}_{1}\) in Fig. 5. The case when the dynamics is restricted to \(\mathcal {M}_{\lambda }\) is discussed in Sect. 4.2.
Figure 5b shows the extended phase space and how solutions of (4.1) behave about \(\mathcal {M}_{\lambda }\), when instantaneous crossing is not allowed. In Fig. 5b, \(\mathcal {M}_{\lambda }\) is a connected manifold. When a solution of (4.1) arrives at \(\varSigma \), the value of \(\lambda \) must change, so that a trajectory can enter \(\varSigma \). \(\mathcal {M}_{\lambda }\) is only invariant for constant \(\lambda \), and therefore a trajectory (denoted by dotted lines) will not continue on \(\mathcal {M}_{\lambda }\) while also in \(\varSigma \). Once a trajectory has left \(\varSigma \), it will be attracted to \(\mathcal {M}_{\lambda }\) as per Assumption (A4).
In the following sections, we discuss how the departure of a trajectory from \(\mathcal {M}_{\lambda }\) 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 \(\mathcal {M}_{\lambda }\) makes a difference and trajectories can no longer cross \(\varSigma \) instantaneously.
4.2 The Skeleton Model
Having assumed the existence of an invariant manifold \(\mathcal {M}_{\lambda }\), it is natural to consider the dynamics on \(\mathcal {M}_{\lambda }\) in the presence of switching. This can be done by substituting the immersion \(\varvec{W}\) into the full problem (4.1) and disregarding that \(\lambda \) may not be constant on \(\varSigma \). We start with the switching function
In contrast to Sect. 3, the switching function (4.12) depends on \(\lambda \), and therefore the closures described in Sect. 3 may not apply. Using the vector field (4.5) on \(\mathcal {M}_{\lambda }\) and (4.12), we obtain
where \(\varvec{y}(t)\in G\) for all \(t\in [s,s+\varDelta )\), \(\varDelta >0\).
Definition 4.7
Equation (4.13) is called the skeleton model of (4.1) on the invariant manifold \(\mathcal {M}_{\lambda }\).
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 \(\lambda \) 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 \(\lambda \). 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
In the extended state space \(\left( \varvec{y},\lambda \right) \in G\times [-\,1,1]\), \(\varSigma _{0}^{\pm }\) is the boundary of the n-dimensional manifold
\(\varSigma _{0}^{\pm }\) cannot intersect each other in the extended state space.
When \(h_{0}(\varvec{y},\lambda )\ne 0\), the trajectories are described by the vector field
Otherwise, we have an index-1 differential algebraic equation
A unique solution to (4.17) is guaranteed by the implicit function theorem if \(D_{2}h_{0}(\varvec{y},\lambda )\ne 0\), so that there is a unique \(C^{p}\) smooth function \(\lambda (\varvec{y})\) satisfying \(h_{0}(\varvec{y},\lambda (\varvec{y}))=0\). This also implies that \(\lambda (t)=\lambda (\varvec{y}(t))\) is continuous, and hence there is no discontinuity of \(\lambda \) when a trajectory reaches \(\varSigma _{0}^{\pm }\) transversely. Trajectories must spend nonzero time on \(\varSigma _{0}\) in order to keep \(\lambda \) continuous. This short argument highlights a major difference between PWS models described in Sect. 3, where we have \(D_{2}h(\varvec{y},\lambda )=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}(\varvec{y},\lambda )=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 \(\dot{\lambda },\) which yields the differential form of (4.17) for \(\left( \varvec{y},\lambda \right) \in \varSigma _{0}\), that is,
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 \(\varSigma _{0}^{\pm }\) 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 \(\varSigma _{0}^{\pm }\).
Theorem 4.8
Consider a point \(\left( \varvec{y}^{\star },\lambda ^{\star }\right) \in \varSigma _{0}^{\pm }\) and assume that
Further, assume a solution \((\varvec{y}(t),\lambda (t))\), for \(t\in I=(-\,\delta ,0]\) (or \(t\in I=[0,\delta )\)) of either Eq. (4.16) or Eq. (4.20) that reaches \(\left( \varvec{y}^{\star },\lambda ^{\star }\right) \) at \(t=0\). The corresponding trajectory is defined as \(\mathcal {T}=\left\{ (\varvec{y}(t),\lambda (t)):t\in I\right\} \). Trajectory \(\mathcal {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:
-
1.
\(\mathcal {T}\) is not tangent to \(\varSigma _{0}^{\pm }\), i.e., \(D_{1}h_{0}\left( \varvec{y}^{\star },\lambda ^{\star }\right) \varvec{f}\left( \varvec{y}^{\star },\lambda ^{\star }\right) \ne 0\)
-
2.
\(\mathcal {T}\) is tangent to \(\varSigma _{0}^{\pm }\) and the order of the tangency is less than the smoothness order (\(C^{p}\)) of \(h_{0}\). In other words, there exists \(0<\ell \le p\) such that
$$\begin{aligned} \frac{\mathrm {d}^{\ell }}{\mathrm {d}t^{\ell }}h_{0}(\varvec{y}(t),\lambda ^{\star })\vert _{t=0}\ne 0. \end{aligned}$$(4.22)
Proof
The proof can be found in Appendix B. \(\square \)
Remark 4.9
Theorem 4.8 excludes the case \(D_{2}h_{0}(\varvec{y},\lambda )>0\). For \(D_{2}h_{0}(\varvec{y},\lambda )>0\), transverse trajectories (case 1 of Theorem 4.8) cannot cross \(\varSigma _{0}^{\pm }\). Tangential trajectories with even \(\ell \) may have multiple continuation, which is the case of the Teixeira singularity (Colombo and Jeffrey 2011). Tangential trajectories with odd \(\ell \) cannot cross \(\varSigma _{0}^{\pm }\), similar to transverse trajectories. To investigate the case of \(D_{2}h_{0}(\varvec{y},\lambda )>0\) in detail, a definition of how trajectories move along \(\varSigma _{0}^{\pm }\) (with \(\lambda =\pm \,1\)) is also required, which falls outside of the scope of this paper.
4.3 Dynamics About Manifold \(\mathcal {M}_{\lambda }\) Due to Switching
This section describes a correction to the skeleton model (4.13) that resolves the dynamics in the neighborhood of \(\mathcal {M}_{\lambda }\) up to linear order. The correction is necessary, because the \(\dot{\lambda }=0\) assumption does not hold: Eq. (4.20) states that \(\lambda \) varies on \(\varSigma _{0}\). The correction that is introduced here captures trajectories that depart from \(\mathcal {M}_{\lambda }\) when \(h=0\) (see dashed line in Fig. 5b).
Let us suppose that
where \(\varvec{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 \(\varvec{z}\), we substitute (4.23) into (4.1), while taking into account that \(\lambda \) is a function of time. This substitution yields
We assume that \(\varvec{z}\) is a small deviation from \(\mathcal {M}_{\lambda }\) and Taylor expand \(\varvec{F}(\varvec{W}(\varvec{y},\lambda )+\varvec{z},\lambda )\) in \(\varvec{z}\) about \(\varvec{z}=\varvec{0}\), that is,
The expansion (4.25) when substituted into (4.24) yields
We now use the invariance Eq. (4.3) and the dynamics on \(\mathcal {M}_{\lambda }\) as given by (4.5) and notice that two terms cancel in (4.26), so that we get
Combining the skeleton model (4.13) with (4.27) yields the corrected model
where \(\varvec{A}_{1}(\varvec{y},\lambda )=D_{1}\varvec{F}(\varvec{W}(\varvec{y},\lambda ),\lambda )\) is defined in Assumption (A2). A unique solution of (4.28) is assumed in (A3) with a continuously differentiable \(\lambda \). 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
and
When a trajectory is restricted to \(\varSigma \), the solution must satisfy
Similar to the skeleton model, we evaluate how h changes in time and we restrict this change to zero on \(\varSigma \) to find an equation for \(\lambda \) [cf. Eq. (4.18)]. To evaluate Eq. (4.29), we use the following lemma.
Lemma 4.10
Assume (A3) and that \(\lambda \) is continuously differentiable and \(\varvec{y}\), \(\varvec{z}\) satisfy the differential equations
on the interval \(t\in [s,s+\epsilon )\), \(\epsilon >0\) with an initial condition \(\varvec{y}(s)\in G\), \(\varvec{z}(s)\in \varvec{\mathcal {D}}\). Then the right-side derivative of h as a function of time is calculated as
where
Proof
The proof can be found in Appendix C. \(\square \)
Remark 4.11
The quantity \(d^{\pm }(\varvec{y},\varvec{z},\lambda )\) in (4.31) measures the discontinuity of the convolution kernel \(\varvec{K}\) at \(t=s\). A discontinuous \(\varvec{K}\) is possible, because \(D_{2}\varvec{W}(\varvec{y},\lambda )\in \varvec{X}\backslash \varvec{Z}\), and the continuity Assumption (A3) does not apply at \(t=s\). Such a discontinuity allowed us to find a differential equation for \(\lambda \) in Sect. 2.
Definition 4.12
We call the quantity \(d^{\pm }(\varvec{y},\varvec{z},\lambda )\) in Eq. (4.31) the normal discontinuity gap.
We also define two other quantities that will be useful later. These are
and therefore we have the identity \(d^{\pm }(\varvec{y},\varvec{z},\lambda )=d^{+}(\varvec{y},\varvec{z},\lambda )-d^{-}(\varvec{y},\varvec{z},\lambda )\).
We now find the governing equation of the dynamics on \(\varSigma \). We solve equation \(\frac{\mathrm {d}}{\mathrm {d}t^{+}}h=0\), where \(\frac{\mathrm {d}}{\mathrm {d}t^{+}}h\) is given by (4.30) for \(\dot{\lambda }\), which yields
The trajectories of Eq. (4.34) are concatenated with trajectories of
along the boundaries \(\varSigma ^{\pm }\) and form the trajectories of the corrected model (4.28). The following theorem provides a sufficient condition for a unique continuation of trajectories through \(\varSigma ^{\pm }\).
Theorem 4.13
Assume (A1)–(A5). A trajectory \(\mathcal {T}\) of either (4.34) or (4.35) with an end point \(\left( \varvec{y},\varvec{z},\lambda \right) \in \varSigma ^{\pm }\) 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:
-
1.
$$\begin{aligned} d^{\pm }(\varvec{y},\varvec{z},\lambda )>0, \end{aligned}$$(4.36)
-
2.
\({Dh}(\varvec{W}(\varvec{y},\lambda )+\varvec{z})\cdot \varvec{U}(t,s)\varvec{z}\) is continuously differentiable with respect to t for \(t\ge s\) and
-
3.
one of the vector fields, (4.34) or (4.35) is not tangent to \(\varSigma ^{\pm }\), that is,
$$\begin{aligned} {Dh}(\varvec{W}(\varvec{y},\lambda )+\varvec{z})\cdot \left( D_{1}\varvec{W}(\varvec{y},\lambda )\varvec{f}(\varvec{y},\lambda )+\varvec{A}_{1}(\varvec{y},\lambda )\varvec{z}\right) \ne 0. \end{aligned}$$(4.37)
Proof
The proof of Theorem 4.13 can be found in Appendix C. \(\square \)
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 \(\left\| \varvec{z}\right\| \ll 1\), the linear correction is accurate. Because on \({{\mathcal {M}}}_{\lambda }\), we have \(\varvec{z}=\varvec{0}\), when a trajectory enters \(\varSigma \), the rate of change of \(\varvec{z}\) is determined by \(\dot{\lambda }\). The magnitude of \(\dot{\lambda }\) depends on the \(\varvec{f}\) and \(d^{\pm }\). Smaller \(d^{\pm }\) makes \(\lambda \) faster. The value of \(d^{\pm }\) is not necessarily a small parameter, and therefore the deviation from \({{\mathcal {M}}}_{\lambda }\) can stay small. For the linear string \(d^{\pm }=\frac{1}{2}\). In the literature of regularized PWS systems (Jeffrey 2015; Kristiansen and Hogan 2015), to stay close to the skeleton model, fast \(\lambda \) is assumed.
Remark 4.15
If \(d^{\pm }(\varvec{y},\varvec{z},\lambda )=0\), the dynamics about \(\mathcal {M}_{\lambda }\) as captured by variable \(\varvec{z}\) can only have a second-order effect on h due to the nonlinearity of h. Therefore, (4.30) is independent of \(\dot{\lambda }\) and \(\frac{\mathrm {d}}{\mathrm {d}t^{+}}h=0\) cannot be solved for \(\dot{\lambda }\). When \(d^{\pm }(\varvec{y},\varvec{z},\lambda )=0\), the corrected model (4.28) needs a closure, such as Filippov’s or Utkin’s. \(d^{\pm }(\varvec{y},\varvec{z},\lambda )=0\) occurs when \(\varvec{U}\) is strongly continuous on the whole of \(\varvec{X}\), i.e., \(\varvec{Z}=\varvec{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 \(\varvec{N}\). The equation of the exact correction can be written as
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
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 \(\varvec{U}(t,s)\) is only defined on \(\varvec{Z}\), but
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 \(\varvec{N}(\varvec{y},\lambda ;\cdot ):\varvec{\mathcal {D}}\rightarrow \varvec{X}\) be bounded.
4.4 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\le \varepsilon \le 1\) and denote the dependence on \(\varepsilon \) by a subscript, that is \(\varvec{F}_{\varepsilon }\). Here, the \(\varepsilon =0\) limit is represented by the skeleton model (4.13) and \(\varepsilon =1\) refers to the corrected model (4.28). Naturally, the immersion \(\varvec{W}_{\varepsilon }(\varvec{y},\lambda )\) of the invariant manifold also depends on \(\varepsilon \), which implicitly assumes that \({{\mathcal {M}}}_{\lambda }\) persists for \(0\le \varepsilon \le 1\). Whenever we write \(\varvec{F}_{0}\) or \(\varvec{W}_{0}\), we mean the \(\varepsilon =0\) limit.
Let us define the scaled Frechet derivative as
With this notation the corrected model (4.28) becomes
Changing the time scales by introducing \(t=\varepsilon \theta \) we get
where \(\mathring{\;}\) stand for \(\nicefrac {\mathrm {d}}{\mathrm {d}\theta }\). When setting \(\varepsilon =0\) we arrive at the layer system
which stipulates that variable \(\varvec{y}\) is constant along trajectories. We assume the following:
- (\(\overline{\mathbf{A3 }}\)):
-
Assumptions (A3) holds when (4.6) is replaced by (4.42) for all \(\varepsilon \in [0,1]\). The unique solution of (4.42) can be written as
$$\begin{aligned} \varvec{z}(t)=\varvec{U}_{\varepsilon }(t,s)\varvec{z}(s)-\int _{s}^{t}\varvec{K}_{\varepsilon }(t,\tau )\dot{\lambda }(\tau )\mathrm {d}\tau . \end{aligned}$$ - (\(\overline{\mathbf{A4 }}\)):
-
Assumption (A4) holds when (4.6) is replaced by (4.42) and \(\sigma _{s}<-\varepsilon \sigma _{c}\).
- (\(\overline{\mathbf{A5 }}\)):
-
The perturbation \(D_{2}\varvec{W}_{0}(\varvec{y},\lambda )\) acts in the invariant normal bundle of \(\mathcal {M}_{\lambda }\), that is,
$$\begin{aligned} \varPi ^{c}(\varvec{y},\lambda )D_{2}\varvec{W}_{0}(\varvec{y},\lambda )=\varvec{0}. \end{aligned}$$(4.44)
Remark 4.18
As a consequence of (\(\overline{\hbox {A3}}\)) and (\(\overline{\hbox {A4}}\)), \(\varvec{A}_{0}(\varvec{y},\lambda )\) has an n-dimensional kernel spanned by \(D_{1}\varvec{W}_{0}(\varvec{y},\lambda )\), and \(\varPi ^{c}(\varvec{y},\lambda )\varvec{A}_{0}(\varvec{y},\lambda )=\varvec{0}\). Because of (\(\overline{\hbox {A5}}\)) and for \(t\ge s\), we also have
We investigate the non-smooth dynamics for \(\varepsilon =0\). The case of constant \(\lambda \) is trivial, because we have assumed that \({{\mathcal {M}}}_{\lambda }\) is attracting for \(0\le \varepsilon \le 1\). Next we consider the dynamics in \(\varSigma \), which is described by
Any point in \({{\mathcal {M}}}_{\lambda }\cap \varSigma \), i.e., \(\varvec{y}\in G\), \(\varvec{z}=\varvec{0}\), \(\lambda \in [-\,1,1]\) is an equilibrium of (4.45); therefore, \({{\mathcal {M}}}_{\lambda }\) is invariant under all the dynamics for \(\varepsilon =0\). It is however not obvious whether \({{\mathcal {M}}}_{\lambda }\cap \varSigma \) is attracting for \(\varepsilon =0\), which is addressed by the next theorem.
Theorem 4.19
Assume (A1),(A2),(\(\overline{\hbox {A3}}\)),(\(\overline{\hbox {A4}}\)) and \(d^{-}(\varvec{y},\varvec{0},\lambda )\ne 0\). Let \({{\mathcal {M}}}_{\lambda }^{\mathrm{crit}}\) be a compact set, such that
where
Then, \({{\mathcal {M}}}_{\lambda }^{\mathrm{crit}}\) is a normally hyperbolic and attracting critical manifold of Eq. (4.41) for \(\varepsilon =0\).
Proof
In order to calculate whether the critical manifold is attracting, we linearize Eq. (4.45) by using \(\lambda =\lambda _{0}+\alpha \) as a perturbation
The initial conditions \(\alpha (0)\) and \(\varvec{z}(0)\) are linked through Eq. (4.50), such that
It is sufficient to show that \(\alpha \) decays, because by assumptions (\(\overline{\hbox {A3}}\)) (\(\overline{\hbox {A4}}\)) and without forcing, the \(\varvec{z}\) component decays to a constant; if the initial condition satisfies \(\varPi ^{s}(\varvec{y},\lambda )\varvec{z}(0)=\varvec{0}\), \(\varvec{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
The asymptotic properties of \(\alpha (t)\) are determined by the poles of (4.52). The poles of the numerator are already given by the spectrum of \(\varvec{A}_{0}(\varvec{y}_{0},\lambda _{0})\), which is assumed to be in the left half of the complex plane because \({{\mathcal {M}}}_{\lambda }\) is attracting for constant \(\lambda \). Therefore, only the roots of the denominator can cause instability, and hence the condition that \(\Delta (s)\) has roots in the left half of the complex plane is sufficient. \(\square \)
Remark 4.20
The proof can be extended to calculate the initial and final values of \(\alpha \). According to the Laplace final value theorem, we have \(\lim _{t\rightarrow \infty }\alpha (t)=\lim _{s\rightarrow 0}s\alpha (s)\). We observe that
by assumption (\(\overline{\hbox {A5}}\)) and Remark 4.18. Therefore, we have
and if \(\varPi ^{c}(\varvec{y},\lambda )\varvec{z}(0)=\varvec{0}\) we also have \(\lim _{t\rightarrow \infty }\alpha (t)=0\). Applying the Laplace initial value theorem to Eq. (4.52) yields
where \(d^{+}(\varvec{y},\varvec{z},\lambda )\) is given by (4.33). We also have \({Dh}\left( \varvec{W}_{0}(\varvec{y}_{0},\lambda _{0})\right) \cdot \varvec{z}(0)=-d^{-}(\varvec{y}_{0},\varvec{0},\lambda _{0})\alpha (0)\) according to (4.50), and therefore \(\lim _{t\downarrow 0}\alpha (t)=\alpha (0)\), which makes \(\alpha \) continuous at \(t=s\).
Remark 4.21
Similar to Remark 4.5, normal hyperbolicity does not imply the persistence of \(\mathcal {M}_{\lambda }^{\mathrm{crit}}\) under variations in \(\varepsilon \). The theorem of Bates, Lu and Zeng Bates et al. (1998) suggests that the evolution operator \(\varvec{U}\) needs to be differentiable (among other conditions) for \(\mathcal {M}_{\lambda }^{\mathrm{crit}}\) to persist for small \(\varepsilon >0\). Note that the nonlinear string example in Sect. 5 generates such a differentiable \(\varvec{U}\) on \(\varvec{Z}\).
Remark 4.22
When both regions of \(\mathcal {M}_{\lambda }\), that is \(\mathcal {M}_{\lambda }\cap \varSigma \) and \(\mathcal {M}_{\lambda }\backslash \left( \mathcal {M}_{\lambda }\cap \varSigma \right) \), persist for \(\varepsilon >0\), they most likely become discontinuous at the boundaries \(\varSigma ^{\pm }\), and hence as a whole, \(\mathcal {M}_{\lambda }\) does not persist. This is because the vector fields are discontinuous. Therefore, for \(\varepsilon >0\), trajectories that followed one part of \(\mathcal {M}_{\lambda }\) must jump to the other part of \(\mathcal {M}_{\lambda }\), which induces fast transients that we are unable to characterize under general settings.
4.5 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^{\pm }\). To rectify the omission of \(d^{\pm }\), the skeleton model is extended by a scalar variable, which represents the dynamics of the convolution kernel \(\varvec{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
- (A6):
-
\(h(\varvec{x})\) is linear, therefore \(h(\varvec{x})=h(\varvec{0})+{Dh}\cdot \varvec{x}\), where Dh is a constant linear functional.
Assumption (A6) allows us to derive a scalar representation of \(\varvec{z}(t)\) without worrying about a varying \({Dh}(\varvec{x})\). The switching between parts of the state space depends on
In what follows, we approximate the scalar-valued \({Dh}\cdot \varvec{z}\) in (4.53) by a convolution integral. Combining Eqs. (4.7), (4.11) and \(\varvec{z}(0)=\varvec{0}\) yields
In order to proceed, either (A5) or time-scale separation with (\(\overline{\hbox {A5}}\)) can be assumed. When (A5) is assumed, we are restricted to use \(\varepsilon =1\) and if (\(\overline{\hbox {A5}}\)) is assumed we set \(\sigma =\sigma _{s}\). Now we can approximate that
which neglects trajectories in the normal bundle of \(\mathcal {M}_{\lambda }\) that are decaying with exponents smaller than \(\varepsilon ^{-1}\sigma \). Note that
By simply defining \(d^{+}(\varvec{y},\lambda )=d^{+}(\varvec{y},\varvec{0},\lambda )\) we get
After defining \(\kappa ={Dh}\cdot \varvec{z}(t)\), we find that the approximation (4.54) satisfies the differential equation
with initial condition \(\kappa (0)=0\). The switching function (4.53) using the new variable \(\kappa \) becomes
We can also redefine the switching manifolds
With this notation, the skeleton model extended with the approximate normal dynamics becomes
Definition 4.23
We call Eq. (4.57) the reduced-order model of (4.1).
When \(h_{\varepsilon }(\varvec{y},\kappa ,\lambda )\ne 0\), the dynamics of \(\kappa \) is decoupled from the rest of the variables and \(\kappa \) exponentially vanishes, because \(\sigma <0\)—due to Assumption (A5) or (\(\overline{\hbox {A5}}\)). When \(h_{\varepsilon }(\varvec{y},\kappa ,\lambda )=0\), we apply the same technique as in Sect. 4.2 to find a differential equation for \(\lambda \). We express that
in \(\varSigma _{\varepsilon }\), where \(d^{-}(\varvec{y},\lambda )=d^{-}(\varvec{y},\varvec{0},\lambda )\) is defined by (4.32). We drop the \((\varvec{y},\lambda )\) arguments and solve (4.55) and (4.58) for \(\dot{\kappa }\) and \(\dot{\lambda }\) to arrive at
which governs the dynamics on \(\varSigma _{\varepsilon }\).
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 \(\mathcal {T}\) of the reduced-order model (4.57) with an end point at \((\varvec{y}^{\star },\kappa ^{\star },\lambda ^{\star })\in \varSigma _{\varepsilon }^{\pm }\) has a unique continuation through \((\varvec{y}^{\star },\kappa ^{\star },\lambda ^{\star })\) if
-
1.
\(d^{\pm }(\varvec{y}^{\star },0,\lambda ^{\star })>0\) as defined by Eq. (4.31) and
-
2.
when trajectory \(\mathcal {T}\) is not tangent to \(\varSigma ^{\pm }\) or trajectory \(\mathcal {T}\) is tangent to \(\varSigma _{\varepsilon }^{\pm }\) and the of order of the tangency is not greater than the smoothness order (\(C^{p}\)) of \(h_{\varepsilon }\), that is, there exists \(0<\ell \le p\) such that
$$\begin{aligned} \frac{\mathrm {d}^{\ell }}{\mathrm {d}t^{\ell }}h_{\varepsilon }(\varvec{y}(t),\kappa (t),\lambda ^{\star })\vert _{\varvec{y}=\varvec{y}^{\star },\kappa =\kappa ^{\star }}\ne 0. \end{aligned}$$
Proof
The proof is the same as for Theorem 4.8 if we replace \(\left( \varvec{y},\kappa \right) \rightarrow \varvec{y}\) and \(d^{\pm }\rightarrow -D_{2}h_{0}(\varvec{y},\lambda )\). \(\square \)
Proposition 4.25
Let
be a compact set for \(\varepsilon =0\). Then, \({{\mathcal {M}}}_{\lambda }^{\mathrm{crit}}\) is an attracting critical manifold of Eq. (4.59) which persists for a sufficiently small \(\varepsilon >0\). The dynamics on the critical manifold is governed by the skeleton model (4.20).
Proof
First, time is rescaled by \(t=\varepsilon \theta \) in Eq. (4.59) which yields
where \(\mathring{\;}\) stands for \(\nicefrac {\mathrm {d}}{\mathrm {d}\theta }\). Setting \(\varepsilon \rightarrow 0\) yields
Assuming initial conditions \(\kappa (0)=\kappa _{0}\) and \(\lambda (0)=\lambda _{0}\) of (4.61) at \(t=0\), we get \(\lim _{t\rightarrow \infty }\kappa (t)=0\) and \(\lim _{t\rightarrow \infty }\lambda (t)=\lambda _{0}-\kappa _{0}/d^{-}\), if \(\frac{d^{-}\sigma }{d^{\pm }}>0\). This means that the critical manifold is attracting and normally hyperbolic. Therefore, \({{\mathcal {M}}}_{\lambda }^{\mathrm{crit}}\) persists for a sufficiently small \(\varepsilon >0\), according to Fenichel (1972).
While \(\kappa =0\) on the critical manifold, \(\lim _{t\rightarrow \infty }\lim _{\varepsilon \rightarrow 0}\varepsilon ^{-1}\kappa (t)\) may not be zero, that is, the limits \(t\rightarrow \infty \) and \(\varepsilon \rightarrow 0\) do not commute. After introducing \(\varepsilon \gamma =\kappa \), we can write that
Setting \(\varepsilon \rightarrow 0\) and some algebraic manipulation yields
which is the same equation as (4.20) of the skeleton model. \(\square \)
Remark 4.26
We know that \(\sigma <0\), because of Assumption (A5) or (\(\overline{\hbox {A5}}\)). If Proposition 4.24 also holds, \({{\mathcal {M}}}_{\lambda }^{\mathrm{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 \(\varSigma _{\varepsilon }\) the critical manifold is likely to be attracting or repelling under the same conditions. The precise statement is in the following proposition.
Proposition 4.27
Assume that \(d^{\pm }(\varvec{y},\varvec{0},\lambda )>0\) and \({Dh}\cdot \varvec{A}_{0}^{-1}(\varvec{y},\lambda )D_{2}\varvec{W}(\varvec{y},\lambda )>0\) along a smooth curve \(\gamma =\left\{ \left( \varvec{y}(\alpha ),\lambda (\alpha )\right) \in \varSigma _{0}:\alpha \in (-\,\delta ,\delta )\right\} \) with \(\left( \varvec{y},\lambda \right) \in C^{1}\left( (-\,\delta ,\delta ),\varSigma _{0}\right) \) and \(\delta >0\). For \(\varepsilon =0\), the stability of equilibria along \(\gamma \) changes through a zero root (saddle-node bifurcation) at the same value(s) of \(\alpha \in (-\,\delta ,\delta )\) for both systems (4.59) and (4.45).
Proof
Because of the assumption \(d^{\pm }(\varvec{y},\varvec{0},\lambda )>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 \(\Delta (s)\) as given by (4.47), we note that \(\Delta (0)=d^{-}(\varvec{y},\lambda )\). This means that we have a zero root of \(\Delta (s)\) when \(d^{-}=0\). Next, we show that this zero root of \(\Delta (s)\) becomes of the same sign as \(d^{-}\) as \(\varvec{y}(\alpha ),\lambda (\alpha )\) changes along \(\gamma \). Let us now assume that at \(\alpha =0\), we have \(d^{-}(\varvec{y}(0),\lambda (0))=0\) and denote the root of \(\Delta \) that smoothly depends on \(\alpha \) by \(s:(-\,\delta ,\delta )\rightarrow \mathbb {R}\) and for which \(s(0)=0\). We denote the derivative with respect to \(\alpha \) by \(^{\prime }\) and calculate the derivative of s from the definition (4.47), that is,
where we omitted that \(\varvec{y}\), \(\lambda \) are evaluated at \(\alpha =0\). We also calculate the derivative \(d^{-\prime }(\varvec{y},\lambda )=D_{\text {1}}d^{-}(\varvec{y},\lambda )\varvec{y}^{\prime }+D_{\text {2}}d^{-}(\varvec{y},\lambda )\lambda ^{\prime }\) and notice that the derivatives \(s^{\prime }\) and \(d^{-\prime }\) have the same sign when \({Dh}\cdot \varvec{A}_{0}^{-1}(\varvec{y},\lambda )D_{2}\varvec{W}(\varvec{y},\lambda )>0\), which proves the proposition. \(\square \)
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 \(\varepsilon \) 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 \(\varepsilon >0\) the trajectories are still unique, but there is no critical manifold in \(\varSigma \) 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 \(\varepsilon =1\). If this dynamics has a Lyapunov exponent \(\sigma _{0}\) such that \(\sigma _{s}<-\left| \sigma _{0}\right| \) as in Assumption (A4), then this perturbed vector bundle could be attached to \(\mathcal {M}_{\lambda }\), which would become a normally hyperbolic invariant manifold of the corrected model (4.28) in \(\varSigma \).
5 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 \(\mathcal {M}_{\lambda }\) is a linear subspace about an equilibrium that depends nonlinearly on the switching parameter \(\lambda \). The dynamics within the invariant manifold given by \(\varvec{f}(\varvec{y},\lambda )\) and the switching function \(h_{0}(\varvec{y},\lambda )\) 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(\xi ,t)\), where \(\xi \in [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 \(\xi \) and dots to denote differentiation in time. The dimensionless equation under our simplifying assumptions is
where T is the tension within the string and \(\varGamma \) controls the nonlinearity of the string. The boundary conditions are \(u(0,t)=u(1,t)=0\) and
where \(\mu \) is a friction coefficient, \(\xi ^{\star }\) is the position of the contact point with the bow and \(u^{\prime }(\xi ^{\star }-,t)\), \(u^{\prime }(\xi ^{\star }+,t)\) represent the left and right derivative of u with respect to \(\xi \) at \(\xi ^{\star }\), 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 \(\mu \lambda \). The switching parameter \(\lambda \) decides the direction of the friction force and therefore changes sign as the relative velocity \(h=v_{0}-\dot{u}(\xi ^{\star },t)\) between the bow and the string reverses, that is,
To further simplify Eq. (5.1), we use second-order Taylor expansion, that is, \(\sqrt{1+u'^{2}}\approx 1+\frac{1}{2}u'^{2}\), which gives us the equation
with boundary conditions
We require that \(u(\cdot ,t)\in \mathrm {Lip}\left( [0,1],{{\mathbb {R}}}\right) \), i.e., \(u(\cdot ,t)\) is Lipschitz continuous, which allows a finite contact force on the string. We define the operator
The square root of \(\mathfrak {D}^{2}\), can be represented on the series \(u=\sum a_{k}\sin k\pi \xi \) by \(\mathfrak {D}u=\sum k\pi a_{k}\sin k\pi \xi \). Note that \(\mathfrak {D}\) is not producing the first-order derivative. To represent all boundary conditions, we define a restricted \(\mathfrak {D}^{2}\) as
We also introduce damping with a constant damping ratio \(\beta \in [0,1)\) for all vibration modes which transforms Eq. (5.3) into
Let us define \(\varvec{x}_{1}=u(\cdot ),\,\varvec{x}_{2}=\dot{u}(\cdot )\) and \(\varvec{x}=(\varvec{x}_{1},\varvec{x}_{2})^{\mathrm{T}}\), and hence we can write the system (5.3) as the infinite-dimensional dynamical system
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^{\infty }\) 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 \(\lambda \), which is denoted by \(\varvec{x}^{\star }\).
-
2.
Find the smoothest two-dimensional spectral submanifold (Haller and Ponsioen 2016) \(\mathcal {M}_{\lambda }\) about \(\varvec{x}^{\star }\), corresponding to the pair of complex conjugate eigenvalues with the least negative real part. The immersion of the manifold is denoted by \(\varvec{W}:{{\mathbb {R}}}^{2}\times [-\,1,1]\rightarrow \varvec{X}\). Assume a function \(\varvec{y}^{\star }:{{\mathbb {R}}}^{2}\times [-\,1,1]\rightarrow {{\mathbb {R}}}^{2}\), which shifts the parametrization of \(\mathcal {M}_{\lambda }\), such that \(\varvec{W}(\varvec{y},\lambda )=\varvec{W}_{ fix }(\varvec{y}+\varvec{y}^{\star }(\varvec{y},\lambda ),\lambda )\), where \(\varvec{W}_{ fix }\) is just one parametrization of \(\mathcal {M}_{\lambda }\). \(\varvec{y}^{\star }\) is an unknown and will be calculated in step 4.
-
3.
Introduce an artificial parameter \(\varepsilon \) that slows down the dynamics on \(\mathcal {M}_{\lambda }\) to standstill at \(\varepsilon =0\) and has no effect at \(\varepsilon =1\). Then for \(\varepsilon =0\), calculate the invariant normal bundle of \(\mathcal {M}_{\lambda }\), which is formed by the subspace orthogonal to the kernel of the adjoint \(\varvec{A}_{0}^{\star }(\varvec{y},\lambda )\) at each point on \(\mathcal {M}_{\lambda }\).
-
4.
Choose a coordinate shift \(\varvec{y}^{\star }\) so that \(D_{2}\varvec{W}\) falls into the invariant normal bundle of \(\mathcal {M}_{\lambda }\) at \(\varepsilon =0\), that is, \(\varvec{W}\) satisfies assumption (\(\overline{\hbox {A4}}\)). This now fully specifies the immersion \(\varvec{W}\).
-
5.
Obtain the skeleton model by substituting the immersion \(\varvec{W}\) into (5.5).
-
6.
Calculate the normal discontinuity gap from the dynamics in the invariant normal bundle of \(\mathcal {M}_{\lambda }\). Also determine \(\sigma \), the rate of convergence of the trajectory in the normal bundle with initial condition \(D_{2}\varvec{W}\).
Proposition 5.1
Following the six steps above yields the reduced-order model of Eqs. (5.5) and (5.6) in the form of
with switching function
The normal discontinuity gap is
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
where
Proof
These results are proven in Lemmas 5.2–5.4, 5.6 and 5.7. \(\square \)
5.1 The Invariant Manifold and Its Parametrization
We identify the invariant manifold \(\mathcal {M}_{\lambda }\) with the spectral submanifold (Haller and Ponsioen 2016) of the string’s equilibrium corresponding to its first natural frequency. When \(\lambda \) 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.
Lemma 5.2
The immersion of the invariant manifold \(\mathcal {M}_{\lambda }\) about the equilibrium, as specified in steps 1 and 2 of the model reduction process, is
where \(\gamma (\lambda )\) is given by Eq. (5.8). The coordinate shift \(\varvec{y}^{\star }=\left( y_{1}^{\star },y_{2}^{\star }\right) ^{\mathrm{T}}\) is not yet known.
Proof
We choose the representation of the invariant manifold as
where
A substitution of \(\varvec{W}\) into the invariance Eq. (4.3) shows that \(\varvec{W}\) is indeed an immersion of an invariant manifold and corresponds to the first natural frequency. Because \(\varvec{W}\) is linear in \(\varvec{y}\), \(\mathcal {M}_{\lambda }\) is also the smoothest invariant manifold.
The equilibrium \(\varvec{x}^{\star }(\lambda )\) is calculated by setting the time derivative to zero in Eq. (5.3), which yields
In Eq. (5.12), \(\frac{\varGamma }{2}\int _{0}^{1}\varvec{x}_{1}'^{2}\mathrm {d}\xi \) is independent of \(\xi \), and therefore integrating (5.12) twice and applying the boundary conditions, we get
which still needs to be solved for \(\varvec{x}_{1}^{\star }\). We define
which yields
Physically, \(\gamma (\lambda )\xi ^{\star }\left( 1-\xi ^{\star }\right) \) is the displacement of the string at the contact point at the equilibrium. To find the equation for \(\gamma \), we substitute (5.15) into (5.14). We then evaluate the integral in (5.14), that is,
so that Eq. (5.14) becomes
Equation (5.16) can be solved for \(\gamma \) with a single real solution, which is given by Eq. (5.8) that makes the equilibrium fully specified. Figure 7 shows the values of \(\gamma \) 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). \(\square \)
5.2 Linearized Dynamics About the Invariant Manifold
The linearized dynamics about \(\mathcal {M}_{\lambda }\) is characterized by the Frechet derivative \(\varvec{A}_{1}(\varvec{y},\lambda )\) of Eq. (5.5), which is calculated here.
Lemma 5.3
The Frechet derivative of \(\varvec{F}\) evaluated on \(\mathcal {M}_{\lambda }\) is
where
and the instantaneous square of the wave speed on \(\mathcal {M}_{\lambda }\) at the contact point is
and
The domain of definition of \(\varvec{A}_{1}(\varvec{y},\lambda )\) is
and
Proof
The only term in Eq. (5.5) not already linear is
which we now linearize about a general point \((\overline{\varvec{x}}_{1},\overline{\varvec{x}}_{2})\in \mathcal {M}_{\lambda }\). The expression (5.23) is a product, and hence we use the product rule when differentiating it with respect to \(\varvec{x}_{1}\). First, we linearize (5.23) about \(\overline{\varvec{x}}_{1}\) and get
where we have used that \(\varvec{z}_{1}\) must vanish at the boundaries \(\xi =0,1\). Therefore, the first-order Taylor expansion of (5.23) is
The value \(\overline{\varvec{x}}_{1}\) is the first component of the immersion of \(\mathcal {M}_{\lambda }\),
and when applying \(\mathfrak {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 \(\mathcal {M}_{\lambda }\) and at the contact point as
which, after substitution becomes (5.19). Gathering all linear terms, the Frechet derivative on the invariant manifold becomes Eq. (5.17).
The domain of definition of \(\varvec{A}_{1}(\varvec{y},\lambda )\) must now include that \(\varvec{x}_{1}^{\prime }(\xi ^{\star }-)-\varvec{x}_{1}^{\prime }(\xi ^{\star }+)=0\), because the equilibrium is already included in the definition of \(\varvec{W}\). However, \(\varvec{x}_{1}^{\prime }(\xi ^{\star }-)-\varvec{x}_{1}^{\prime }(\xi ^{\star }+)=0\) is already satisfied if \(\varvec{x}_{1}^{\prime \prime }\in L^{\infty }([0,1],\mathbb {R})\), and therefore we arrive at (5.21). To determine the closure \(\overline{\varvec{\mathcal {D}}},\) we first note that \(\varvec{x}_{1}^{\prime }\) must be Lipschitz continuous, and therefore \(\varvec{x}_{1}\) is Lipschitz continuously differentiable. The closure of such functions in the Lipschitz norm is the continuously differentiable functions \(C^{1}([0,1],\mathbb {R})\). Since \(\mathfrak {D}\) is the square root of \(\mathfrak {D}^{2}\), they have the same domain of definition, and therefore \(\varvec{x}_{2}\) is Lipschitz continuously differentiable in \(L^{\infty }\). The closure for this set in the \(L^{\infty }\) norm is the set of continuous functions \(C^{0}([0,1],\mathbb {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 \(\varvec{Z}\ne \varvec{X}\). \(\square \)
5.3 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 \(\mathcal {M}_{\lambda }\). We therefore introduce such a scaling by artificially constructing \(\varvec{A}_{\varepsilon }(\varvec{y},\lambda )\) such that for \(\varepsilon =1\), we recover the original dynamics and for \(\varepsilon =0\) the tangential dynamics becomes \(\dot{\varvec{y}}=\varvec{0}\) when time is rescaled. This allows us to calculate the invariant normal bundle of \(\mathcal {M}_{\lambda }\) at \(\varepsilon =0\) and determine the parametrization \(\mathcal {M}_{\lambda }\) (i.e., the unknown coordinate shift \(\varvec{y}^{\star }(\varvec{y},\lambda )\) in the immersion of \(\mathcal {M}_{\lambda }\)) such that \(D_{2}\varvec{W}(\varvec{y},\lambda )\) is strictly in the invariant normal bundle of \(\mathcal {M}_{\lambda }\).
Lemma 5.4
Applying steps 3 and 4 of the model reduction procedure, we find that the coordinate shift \(\varvec{y}^{\star }=\left( y_{1}^{\star },y_{2}^{\star }\right) ^{\mathrm{T}}\) becomes
and
whose solution is
where
Using the coordinate shift (5.25), the square of the instantaneous wave speed (5.19) becomes (5.7).
Proof
Even though the mode shapes of the nonlinear string are the orthogonal harmonic functions \(\sin k\pi \xi \), the contact force \(\lambda \) at \(\xi ^{\star }\) makes these modes intricately coupled. This coupling is represented by the term \(\gamma (\lambda )\varvec{z}_{1}(\xi ^{\star })\) in Eq. (5.18). Nevertheless, we project \(\varvec{A}_{1}(\varvec{y},\lambda )\) into two subspaces using the projection operators \(\varvec{P}:\varvec{X}\rightarrow T_{\varvec{y}}\mathcal {M}_{\lambda }\),
and \(\varvec{Q}=\varvec{I}-\varvec{P}\). We calculate the projected operators
where we can show that \(\varvec{B}_{21}=\varvec{0}\). Introducing time-scale separation is just a multiplication of matrix \(\varvec{B}_{1}\) by \(\varepsilon \), which represents the rescaled linearized dynamics in the (invariant) tangent bundle of \(\mathcal {M}_{\lambda }\). In the new coordinate system, the scaled linear operator becomes
Given the form of \(\varvec{A}_{\varepsilon }\), the bundle projections assume the form
where \(\varvec{C}\) is an unknown operator. Expanding the bundle invariance Eq. (4.9) yields
which must hold for all \(\varvec{v}\), \(\varvec{P}\varvec{v}=\varvec{0}\). Further expanding (5.30), we get an equation for \(\varvec{C}\) in the form of
The solution is \(\varvec{C}=\varvec{B}_{12}\varvec{B}_{2}^{-1}\) at the critical parameter value \(\varepsilon =0\). We can now introduce another coordinate system in which
We denote this transformation by \(\varvec{x}=\varvec{T}\varvec{a}\), where \(\varvec{a}=\left( \varvec{a}_{1},\varvec{a}_{2}\right) ^{\mathrm{T}}\) with \(\varvec{a}_{\ell }=(a_{\ell 1},a_{\ell 2},\ldots )^{\mathrm{T}}\), \(\ell =1,2\). The projections can be written as
In this new coordinate system, we have the operators
where \(\mathrm {diag}_{k=2}^{\infty }k^{2}\) means an infinite diagonal matrix with elements \(k^{2}\) in the diagonal. The inverse \(\varvec{B}_{2}^{-1}\) is represented by
We can now calculate \(\varvec{C}\) or its representation \(\varvec{C}\varvec{T}\varvec{a}=\varvec{B}_{12}\varvec{B}_{2}^{-1}\varvec{T}\varvec{a}\), which becomes
The derivative of the immersion \(\varvec{W}\), when Fourier expanded is
;hence, the coordinates of \(\varvec{T}^{-1}\varvec{P}D_{2}\varvec{W}(\varvec{y},\lambda )\) and \(\varvec{T}^{-1}\varvec{Q}D_{2}\varvec{W}(\varvec{y},\lambda )\) are
and
, respectively. The constraint (\(\overline{\hbox {A4}}\)), i.e., \(\varPi ^{c}D_{2}\varvec{W}(\varvec{y},\lambda )=\varvec{0}\) gives
which is an equation for \(y_{1}^{\star }(\lambda )\) and \(y_{2}^{\star }(y_{1},\lambda )\). Equation (5.32) is solved for \(D_{2}y_{1}^{\star }\) and \(D_{2}y_{2}^{\star }\), which are then integrated over \(\lambda \). The constant of integration for \(y_{1}^{\star }\) is such that \(y_{1}^{\star }(\varvec{y},0)=0\), which yields (5.25). However, we notice that the constant of integration does not play a role, so we present the simplest formula for \(y_{2}^{\star }(y_{1},\lambda )\), whose derivative is \(D_{2}y_{2}^{\star }(y_{1},\lambda )\) without paying attention to the constant of integration. The result of this integration is (5.27). Evaluating the square of the wave speed with this coordinate shift yields (5.7). \(\square \)
Remark 5.5
For \(\varGamma =0\), we have \(\varvec{C}=\varvec{0}\) and also \(D_{2}y_{2}^{\star }=0\). This implies that for the linear string, the bundle projection is simply \(\varvec{Q}\). The normal bundle is independent of \(\varepsilon \), and there is no need to introduce time-scale separation. Instead of \(\varepsilon \), \(\varGamma \) can be used to track the deformation of the invariant normal bundle, which persists for a sufficiently small \(\varGamma >0\) due to the properties of \(\mathcal {M}_{\lambda }\) (Bates et al. 1998).
5.4 The Vector Field \(\varvec{f}(\varvec{y},\lambda )\) on the Invariant Manifold
Lemma 5.6
The skeleton model of Eq. (5.5) on the invariant manifold specified by (5.10) and with coordinate shifts (5.25) and (5.27) can be written as
After substituting the immersion (5.10), the switching function (5.6) becomes
Proof
The dynamics on the invariant manifold \(\mathcal {M}_{\lambda }\) is given by the invariance condition
This is an equation in the tangent bundle of \(\mathcal {M}_{\lambda }\), and therefore it makes sense to project it using \(\varvec{P}\), as defined by (5.28), to find \(\varvec{f}\). We first calculate that
By inverting the matrix (5.35), we find that the reduced vector field is
Next, we substitute the immersion (5.11) so that the vector field (5.5) on the manifold becomes
Substituting (5.37) into (5.36) yields the reduced vector field (5.33). The switching function, defined by Eq. (5.6) after substituting the immersion, becomes (5.34). \(\square \)
5.5 The Normal Discontinuity Gap \(d^{\pm }\) and Decay Rate \(\sigma \)
The normal discontinuity gap \(d^{\pm }\) measures the discontinuity of the correction about the invariant manifold with initial conditions \(D_{2}\varvec{W}(\varvec{y},\lambda )\) at \(t=0\) and determines the uniqueness of solutions according to Theorem 4.13. We calculate \(d^{\pm }\) for the \(\varepsilon \rightarrow 0\) limit, when the dynamics in the normal bundle of \(\mathcal {M}_{\lambda }\) becomes autonomous. Therefore, it is sufficient to evaluate the limit \(\lim _{t\downarrow 0}{Dh}\cdot \mathrm {e}^{\varvec{A}_{0}(\varvec{y},\lambda )t}D_{2}\varvec{W}(\varvec{y},\lambda )\).
Lemma 5.7
The normal discontinuity gap in the limit \(\varepsilon \rightarrow 0\) is
The rate of decay as defined by (4.10) is
Proof
The calculation is carried out using Fourier series, and hence we write the series expansion
which is calculated from (5.31) by substituting (5.25). Since \(D_{2}\varvec{W}(\varvec{y},\lambda )\) is in the invariant normal bundle of the critical manifold, it is sufficient to restrict the dynamics there. We use the decomposition of \(\mathrm {e}^{\varvec{A}_{0}(\varvec{y},\lambda )t}\) as given by (5.29) to arrive at the expression
where
The relevant component of the solution is
The limit \(d^{+}=-\lim _{t\downarrow 0}\varvec{z}_{2}(t)\vert _{\xi =\xi ^{\star }}\) is calculated as
, 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). \(\square \)
Remark 5.8
The normal discontinuity gap \(d^{\pm }\) is a local property of the string; it depends on material properties and the tension in the string. However, \(d^{\pm }\) is independent of the boundary conditions and the position where the string is forced.
5.6 Spectrum of the Normal Dynamics on \(\varSigma \)
We use Theorem 4.19 to find out whether there exists an attracting critical manifold.
Lemma 5.9
The characteristic function determining the stability of the critical manifold \({{\mathcal {M}}}_{\lambda }^{\mathrm{crit}}\) within \(\varSigma \) is given by
Proof
The characteristic function (4.47), whose roots define stability, is formally written as
It is possible to find a convergent series expansion of \(\Delta (s)\) by using Fourier series. Let us, for the moment, define \(\varvec{x}=\left( s-\varvec{A}_{0}(\varvec{y},\lambda )\right) ^{-1}D_{2}\varvec{W}(\varvec{y},\lambda )\), which is obtained by solving
for \(\varvec{x}\). 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 \pi \xi \) and \(x_{2\ell ,k}\) are the coefficients of \(\sin k\pi \xi \) in the Fourier expansion of \(\varvec{x}\). Now expanding Eq. (5.43) and using the notation (5.44) gives
and
The solution to Eq. (5.46) for the \(k\ge 2\) Fourier coefficients is
and
For the first Fourier coefficients, the solution of Eq. (5.45) is \(sx_{11}=0\) and
The series expansion of the characteristic function \(\Delta (s)\) using notation (5.44) is
and substituting system parameters yields (5.42). \(\square \)
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 \(\Delta (s)\) are perturbations of the eigenvalues of \(\varvec{A}_{0}(\varvec{y},\lambda )\) 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 (\(\varGamma =0\)), the green root near the origin in Fig. 8 remains at the origin, and therefore the normal dynamics is neutrally stable.
Remark 5.10
Continuing from Remark 5.5, we find that for \(\varGamma =0\) and for all \(\varepsilon \in [0,1]\) the characteristic function (5.42) is valid due to \(\varvec{A}_{1}\) being constant. Let us denote the zero root of the characteristic function (5.42) for \(\varGamma =0\) by \(\sigma ^{0}\). The invariant vector bundle corresponding to \(\sigma ^{0}\) is then isomorphic to \(\mathcal {M}_{\lambda }\times \mathbb {R}\). For \(\varGamma >0\), the invariant vector bundle of \(\sigma ^{0}\) is continuously perturbed. The perturbation turns \(\sigma ^{0}\) into a small Lyapunov exponent of the now non-autonomous dynamics within the invariant vector bundle. \(\varGamma >0\) can be chosen such that \(\sigma _{s}<-\varepsilon \left| \sigma ^{0}\right| \), that is, the invariant vector bundle is an attracting normally hyperbolic invariant manifold in \(\left( \mathcal {M}_{\lambda }\times \varvec{X}\right) \cap \varSigma \), that is the phase space of the corrected model in \(\varSigma \). The dynamics in the invariant vector bundle is represented by the reduced-order model on \(\mathcal {M}_{\lambda }\times \mathbb {R}\).
5.7 Equivalent Reduced-Order Model on \(\Sigma _{\varepsilon }\)
We now investigate the dynamics of the string on \(\varSigma _{\varepsilon }\). The dynamics outside \(\varSigma _{\varepsilon }\) is given by \(\dot{\varvec{y}}=\varvec{f}(\varvec{y},\lambda )\), \(\varepsilon \dot{\kappa }=\sigma \kappa \) and \(\dot{\lambda }=0\). The skeleton model on \(\varSigma _{\varepsilon }\) 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}^{\star }(y_{1},\lambda )\) as shown by Eq. (5.27). It is possible to eliminate \(y_{2}^{\star }(y_{1},\lambda )\) using the transformation
The vector field on the invariant manifold now involves \(\dot{\lambda }\) in the form of \(\dot{\overline{\varvec{y}}}=\overline{\varvec{f}}(\overline{\varvec{y}},\lambda ,\dot{\lambda })\), where
The main advantage of this formulation is that the function defining the switching manifold \(\varSigma \) becomes simpler, namely
which is independent of \(\lambda \). The stick dynamics on \(\varSigma \) has the same dependence on \(\dot{\lambda }\) as before, because
Note that \(d^{-}=D_{2}h_{0}(\overline{\varvec{y}},\lambda )=-D_{2}y_{2}^{\star }(\overline{y}_{1},\lambda )\) appears in Eq. (5.50), and remains the coefficient of \(\dot{\lambda }\). The last equation we need is (4.55) that describes \(\kappa \). Note that the transformation (5.47) does not change the values of \(d^{+}\), \(d^{-}\) and \(d^{\pm }\) given by (5.41) and (5.38), because they do not depend on \(y_{2}\). As a result, we have
5.8 Dynamics of the Skeleton Model on \(\Sigma _{0}\)
In this section, we explore the dynamics of the skeleton model, which is the same as the dynamics on the critical manifold, when \(\varepsilon =0\) in Eq. (4.59). The dynamics on the critical manifold can be found by setting \(\overline{y}_{2}=v_{0}/\sin \pi \xi ^{\star }\)and \(\dot{\overline{y}_{2}}=0\) so that \(h=0\) and \(\dot{h}=0\), then solving \(\dot{\overline{\varvec{y}}}=\overline{\varvec{f}}(\overline{\varvec{y}},\lambda ,\dot{\lambda })\) as an algebraic equation with (5.48) on the right side for \(\dot{\lambda }\), that is,
As we noted in Theorem 4.8 in Sect. 4.2, solutions of this model pass through the boundaries \(\varSigma ^{\pm }\) if \(D_{2}y_{2}^{\star }(\overline{y}_{1},\lambda )\sin \pi \xi ^{\star }=-d^{-}>0\). To avoid any problem with having the wrong sign of \(d^{-}\), we rescale time by \(D_{2}y_{2}^{\star }(\overline{y}_{1},\lambda )\) and get
whose forward-time solutions always pass through \(\varSigma ^{\pm }\). This allows for a straightforward numerical solution.
Let us first recall, what Utkin’s closure would produce if we disregarded \(D_{2}y_{2}^{\star }(\overline{y}_{1},\lambda )\). The solution would be given by the equation
which is partly algebraic. Figure 10a shows the phase portrait. The dashed orange lines correspond to \(\lambda \) values jumping between either \(\pm \,1\) or the solution of the algebraic constraint in Eq. (5.54). The continuous green line represents \(\lambda \) values that are admissible by the constraint in Eq. (5.54). Since \(\dot{\overline{y}}_{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}^{\star }(\overline{y}_{1},\lambda )\) 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}^{\star }(\overline{y}_{1},\lambda )<0\). The arrows indicate the correct forward direction of time. Black lines indicate when \(D_{2}y_{2}^{\star }(\overline{y}_{1},\lambda )=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 \(\dot{\overline{y}}_{1}=0\). Another nullcline is shown in green, where \(\dot{\lambda }=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 \(\dot{\lambda }=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 \(\lambda =\pm \,1\) yields sliding solutions. However, this implies that solutions coming from either side of \(\varSigma \) cannot enter \(\varSigma \), while \(\lambda \) stays at \(\pm \,1\). Having a fixed value of \(\lambda \) is not physical in a friction oscillator. The singularities within \(\varSigma \) 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 \(\mathcal {M}_{\lambda }\) as described in Sect. 4.5.
5.9 Dynamics of the Reduced-Order Model on \(\varSigma _{\varepsilon }\)
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 \(\mathcal {M}_{\lambda }\). The dynamics on \(\varSigma _{\varepsilon }\) 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 \(\mathcal {M}_{\lambda }\) for \(\varepsilon =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^{\pm }\) 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 \(\varepsilon =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 \(\kappa \). When trajectories start in the shaded part with \(\lambda =1\), where the critical manifold is repelling, they quickly pass to \(\lambda =-\,1\) without much change in \(y_{1}\), while \(\kappa \) exponentially explodes. Trajectories starting with \(\lambda =-\,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 \(\kappa \). This is illustrated in Fig. 11b. After passing the node, trajectories tend to either \(\lambda =\pm \,1\). It is then likely that trajectories will start a violent oscillation between \(\lambda =\pm \,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 \(\varepsilon \) leads to less violent oscillations between \(\lambda =\pm \,1\), which eventually continues without reaching \(\lambda =\pm \,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 \(\varepsilon =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 \(\varepsilon \) sufficiently large the time scale of the normal dynamics (\(\kappa \) variable) becomes much longer than the dynamics of the rest of the variables, and therefore during a stick phase \(\kappa \) 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 \(\Delta (s)\) is a small perturbation of the zero root, and hence it is easily dominated by other time scales. In fact by removing nonlinearity (\(\varGamma =0\)), this root remains zero, and hence \(\kappa \) simply becomes an integral of other quantities without a dynamics of it own. In our example at \(\varepsilon =1\), \(\kappa \) is almost without its own dynamics. The justification why \(\varepsilon \) can be increased to \(\varepsilon =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^{\pm }\) is positive in all parts of the phase space. For the nonlinear string example, \(d^{\pm }\) is the velocity jump of the contact point due to a unit jump in \(\lambda \), i.e., the contact force. Therefore in light of Newton’s second law, it is understandable why \(d^{\pm }>0\). In case we had found \(d^{\pm }=0\), the reduced-order model (4.57), including an approximation of the normal dynamics about \(\mathcal {M}_{\lambda }\), would not be necessary, the skeleton model would be sufficient.
6 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 reduced-order 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).
References
Bates, P.W., Lu, K., Zeng, C.: Existence and persistence of invariant manifolds for semiflows in Banach space. Memoirs of the American Mathematical Society, vol. 135, No. 645. American Mathematical Society, Providence, Rhode Island (1998)
Benner, P., Cohen, A., Ohlberger, M., Willcox, K.: Model Reduction and Approximation: Theory and Algorithms, Computational Science and Engineering. SIAM, Philadelphia (2017)
Cabré, X., Fontich, E., de la Llave, R.: The parameterization method for invariant manifolds I: manifolds associated to non-resonant subspaces. Indiana Univ. Math. J. 52, 283–328 (2003)
Cardin, P.T., da Silva, P.R., Teixeira, M.A.: On singularly perturbed filippov systems. Eur. J. Appl. Math. 24(6), 835–856 (2013). https://doi.org/10.1017/S0956792513000211
Cardin, P.T., de Moraes, J.R., da Silva, P.R.: Persistence of periodic orbits with sliding or sewing by singular perturbation. J. Math. Anal. Appl. 423(2), 1166–1182 (2015). https://doi.org/10.1016/j.jmaa.2014.10.023
Carr, J.: Applications of Centre Manifold Theory, Applied Mathematical Sciences. Springer, Berlin (1981)
Coddington, E.A., Levinson, N.: Theory of Ordinary Differential Equations. McGraw-Hill, New York (1955)
Colombo, A., Jeffrey, M.: Nondeterministic chaos, and the two-fold singularity in piecewise smooth flows. SIAM J. Appl. Dyn. Syst. 10(2), 423–451 (2011). https://doi.org/10.1137/100801846
Coombes, S., Thul, R., Wedgwood, K.C.A.: Nonsmooth dynamics in spiking neuron models. Phys. D Nonlinear Phenom. 241(22), 2042–2057 (2012). https://doi.org/10.1016/j.physd.2011.05.012
da Prato, G., Sinestrari, E.: Non autonomous evolution operators of hyperbolic type. Semigroup Forum 45(1), 302 (1992). https://doi.org/10.1007/BF03025772
Danca, M.F.: Synchronization of switch dynamical systems. Int. J. Bifurc. Chaos Appl. Sci. Eng. 12(8), 1813–1826 (2002). https://doi.org/10.1142/S0218127402005522
DeLellis, P., di Bernardo, M., Liuzza, D.: Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems. Automatica 56, 1–11 (2015). https://doi.org/10.1016/j.automatica.2015.03.003
di Bernardo, M., Budd, C., Champneys, A.: Grazing, skipping and sliding: analysis of the non-smooth dynamics of the DC/DC buck converter. Nonlinearity 11(4), 859 (1998). https://doi.org/10.1088/0951-7715/11/4/007
di Bernardo, M., Budd, C., Champneys, A.R., Kowalczyk, P.: Piecewise-Smooth Dynamical Systems: Theory and Applications. Springer, Berlin (2008)
Diekmann, O., Walther, H.-O., van Gils, S.A., Lunel, S.M.V.: Delay Equations. Springer, New York (1995)
Fenichel, N.: Persistence and smoothness of invariant manifolds for flows. Indiana Univ. Math. J. 21, 193–226 (1972)
Filippov, A.F.: Differential Equations with Discontinuous Righthand Sides: Control Systems, Mathematics and its Applications. Springer, Berlin (1998)
Fridman, L.M.: Singularly perturbed analysis of chattering in relay control systems. IEEE Trans. Autom. Control 47(12), 2079–2084 (2002). https://doi.org/10.1109/TAC.2002.805672
Glass, L., Kauffman, S.A.: The logical analysis of continuous, non-linear biochemical control networks. J. Theor. Biol. 39(1), 103–129 (1973). https://doi.org/10.1016/0022-5193(73)90208-7
Glendinning, P.A., Jeffrey, M.R.: An Introduction to Piecewise Smooth Dynamics. In Advanced Courses in Mathematics—CRM Barcelona. Birkhäuser, Basel (2017)
Gurtin, M.E., Sternberg, E.: A note on uniqueness in classical elastodynamics. Q. Appl. Math. 19(2), 169–171 (1961)
Haller, G., Ponsioen, S.: Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction. Nonlinear Dyn. 86, 1–42 (2016). https://doi.org/10.1007/s11071-016-2974-z
Haller, G., Ponsioen, S.: Exact model reduction by a slow-fast decomposition of nonlinear mechanical systems. Nonlinear Dyn. 90(1), 617–647 (2017). https://doi.org/10.1007/s11071-017-3685-9
Izhikevich, E.M.: Simple model of spiking neurons. IEEE Trans. Neural Netw. 14(6), 1569–1572 (2003). https://doi.org/10.1109/TNN.2003.820440
Jeffrey, M.R.: Smoothing tautologies, hidden dynamics, and sigmoid asymptotics for piecewise smooth systems. Chaos Interdiscip. J. Nonlinear Sci. 25(10), 103125 (2015). https://doi.org/10.1063/1.4934204
Kanamaru, T.: Van der Pol oscillator 2(1), 2202 (2007). https://doi.org/10.4249/scholarpedia.2202. Revision #138698
Kane, C., Repetto, E.A., Ortiz, M., Marsden, J.E.: Finite element analysis of nonsmooth contact. Comput. Methods Appl. Mech. Eng. 180(1), 1–26 (1999). https://doi.org/10.1016/S0045-7825(99)00034-1
Kausel, E.: Fundamental Solutions in Elastodynamics: A Compendium. Cambridge University Press, Cambridge (2006)
Kristiansen, K.U.: Blowup for flat slow manifolds. Nonlinearity 30(5), 2138 (2017). https://doi.org/10.1088/1361-6544/aa6449
Kristiansen, K.U., Hogan, S.J.: On the use of blowup to study regularizations of singularities of piecewise smooth dynamical systems in \(\mathbb{R}^3\). SIAM J. Appl. Dyn. Syst. 14(1), 382–422 (2015). https://doi.org/10.1137/140980995
Kuehn, C.: Multiple Time Scale Dynamics, Applied Mathematical Sciences. Springer, Berlin (2015)
Kuznetsov, Y.: Elements of Applied Bifurcation Theory, 3rd edn. Springer, Berlin (2004)
Levaggi, L.: Infinite dimensional systems, sliding motions. Eur. J. Control 8(6), 508–516 (2002a). https://doi.org/10.3166/ejc.8.508-516
Levaggi, L.: Sliding modes in Banach spaces. Differ. Integr. Equ. 15(2), 167–189 (2002b)
Londoño, J.M., Serino, G., di Bernardo, M.: Existence and stability of limit cycles in a delayed dry-friction oscillator. Nonlinear Dyn. 67(1), 483–496 (2012). https://doi.org/10.1007/s11071-011-9997-2
Magal, P., Ruan, S.: On semilinear Cauchy problems with non-dense domain. Adv. Differ. Equ. 14(11/12), 1041–1084 (2009)
Martins, J.A.C., Oden, J.T.: Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws. Nonlinear Anal. 11(3), 407–428 (1987). https://doi.org/10.1016/0362-546X(87)90055-1
McIntyre, M.E., Woodhouse, J.: On the fundamentals of bowed-string dynamics. Acustica 43(2), 93–108 (1979)
Mestl, T., Plahte, E., Omholt, S.W.: A mathematical framework for describing and analysing gene regulatory networks. J. Theor. Biol. 176(2), 291–300 (1995). https://doi.org/10.1006/jtbi.1995.0199
Metz, J.A.J., Diekmann, O.: The Dynamics of Physiologically Structured Populations. Lecture Notes in Biomathematics. Springer, Berlin (1986)
Neubrander, F.: Integrated semigroups and their applications to the abstract Cauchy problem. Pac. J. Math. 135(1), 111–155 (1988). https://doi.org/10.2140/pjm.1988.135.111
Oestreich, M., Hinrichs, N., Popp, K.: Dynamics of oscillators with impact and friction. Chaos Solitons Fractals 8(4), 535–558 (1997)
Orlov, Y.: Sliding mode control in Banach space. IFAC Proc. Vol. 28(14), 473–476 (1995). https://doi.org/10.1016/S1474-6670(17)46874-1. 3rd IFAC Symposium on Nonlinear Control Systems Design 1995, Tahoe City, CA, USA, 25–28 June (1995)
Orlov, YuV, Utkin, V.I.: Sliding mode control in indefinite-dimensional systems. Automatica 23(6), 753–757 (1987). https://doi.org/10.1016/0005-1098(87)90032-X
Pazy, A.: Semigroups of Linear Operators and Applications to Partial Differential Equations. Springer, New York (1983)
Rao, S.S.: Mechanical Vibrations. Prentice-Hall, Englewood Cliffs (2016)
Sieber, J.: Dynamics of delayed relay systems. Nonlinearity 19(11), 2489 (2006)
Sieber, J., Kowalczyk, P.: Small-scale instabilities in dynamical systems with sliding. Physica D 239(1–2), 44–57 (2010). https://doi.org/10.1016/j.physd.2009.10.003
Smirnov, G.V.: Introduction to the Theory of Differential Inclusions. Volume. 41 of Graduate Studies in Mathematics. American Mathematical Society, Providence (2002)
Sotomayor, J., Teixeira, M.A.: Regularization of discontinuous vectorfields. In: Equadiff 95, International Conference on Differential Equations, pp. 207–223 (1998)
Szalai, R.: Modelling elastic structures with strong nonlinearities with application to stick–slip friction. Proc. R. Soc. A 470(2161), 1 (2014). https://doi.org/10.1098/rspa.2013.0593
Szalai, R., Jeffrey, M.R.: Nondeterministic dynamics of a mechanical system. Phys. Rev. E 90, 022914 (2014). https://doi.org/10.1103/PhysRevE.90.022914
Szalai, R., Osinga, H.M.: Invariant polygons in systems with grazing-sliding. Chaos Interdiscip. J. Nonlinear Sci. 18(2), 023121 (2008)
Szalai, R., Osinga, H.: Arnol’d tongues arising from a grazing-sliding bifurcation. SIAM J. Appl. Dyn. Syst. 8(4), 1434–1461 (2009). https://doi.org/10.1137/09076235X
Teixeira, M.A.: On topological stability of divergent diagrams of folds. Math. Z. 180(2), 361–371 (1982). https://doi.org/10.1007/BF01214176
Thieme, H.R.: “Integrated semigroups” and integrated solutions to abstract Cauchy problems. J. Math. Anal. Appl. 152(2), 416–447 (1990). https://doi.org/10.1016/0022-247X(90)90074-P
Thieme, H.R.: Differentiability of convolutions, integrated semigroups of bounded semi-variation, and the inhomogeneous Cauchy problem. J. Evol. Equ. 8(2), 283–305 (2008). https://doi.org/10.1007/s00028-007-0355-2
Utkin, V.I.: Sliding Modes in Control and Optimization, Communication and Control Engineering Series. Springer, Berlin (1992)
Wechselberger, M.: Existence and bifurcation of canards in \(\mathbb{R}^3\) in the case of a folded node. SIAM J. Appl. Dyn. Syst. 4(1), 101–139 (2005). https://doi.org/10.1137/030601995
Wechselberger, M.: À propos de canards (Apropos canards). Trans. Am. Math. Soc. 364(6), 3289–3309 (2012)
Weiss, D., Küpper, T., Hosham, H.A.: Invariant manifolds for nonsmooth systems. Phys. D Nonlinear Phenom. 241(22), 1895–1902 (2012). https://doi.org/10.1016/j.physd.2011.07.012
Weiss, D., Küpper, T., Hosham, H.A.: Invariant manifolds for nonsmooth systems with sliding mode. Math. Comput. Simul. 110, 15–32 (2015). https://doi.org/10.1016/j.matcom.2014.02.004
Acknowledgements
The author would like to thank Alan R. Champneys and S. John Hogan for the feedback on the manuscript. He would also like to thank the anonymous reviewers who have helped with the clarity of the text and the accuracy of calculations. The author is especially thankful to Galit Szalai, who has proofread the final draft. Funding was provided by Engineering and Physical Sciences Research Council (GB) (Grant No. EP/K003836/1).
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by George Haller.
Appendices
Solution of the Correction Term in the Introductory Example
This appendix details the solution of equation
with initial condition
and boundary conditions \(w(0,t)=w(1,t)=0\). The solution of (A.1) is then substituted into the switching function
which replaces h in Eq. (2.7) of Sect. 2. We also assume that the solution starts with \(\dot{\lambda }(0)=\ddot{\lambda }(0)=0\), which occurs, for example, when \(h\vert _{t=0}\ne 0\).
The solution of Eqs. (A.1) and (A.2) can be expressed using the variation-of-constants formula. Assume that \(w_{\mathrm{h}}(\xi ,t)\) is the solution of the homogeneous equation \(\ddot{w}_{\mathrm{h}}(\xi ,t)=w_{\mathrm{h}}^{\prime \prime }(\xi ,t)\) with initial and boundary conditions, as in
then the solution of (A.1) for the velocity with zero initial condition (A.2) is
Note that we express the velocity here, because that is what appears in the switching function (A.3). Using integration by parts twice transforms Eq. (A.5) into
We now evaluate Eq. (A.6) at \(\xi =\xi ^{\star }\) and notice that a number of terms vanish. From Eq. (A.1), \(\ddot{\lambda }(0)=0\) and \(w_{\mathrm{h}}(\xi ,0)=0\) it follows that \(\ddot{w}_{\mathrm{h}}(\xi ^{\star },0)=0\), and from Eq. (A.4) we obtain \(\dot{w}_{\mathrm{h}}(\xi ^{\star },0)=y^{\star }-u_{0}(\xi ^{\star })\). Further, using \(\dot{\lambda }(0)=0\) brings (A.6) into
And therefore, the switching function (A.3) becomes
The homogeneous solution \(\dot{w}_{\mathrm{h}}(\xi ,t)\) is found using d’Alembert’s method, which states that there are functions f and g such that \(w_{\mathrm{h}}(\xi ,t)=f(\xi +t)+g(\xi -t)\). Given the initial conditions (A.4), we have \(f(\xi )+g(\xi )=0\) and
It follows that \(g(\xi )=-f(\xi )\), and therefore \(2\dot{f}(\xi )=y^{\star }\frac{\sin \pi \xi }{\sin \pi \xi ^{\star }}-u_{0}(\xi )\). Next, we define \(\varphi (\xi )=\ddot{f}(\xi )\), and we get the solution (for the acceleration) in the form
where
Now we evaluate the boundary conditions, that is \(\ddot{w}_{\mathrm{h}}(0,t)=\ddot{w}_{\mathrm{h}}(1,t)=0\), which gives \(\varphi (\xi )=\varphi (-\,\xi )\) and \(\varphi (1+\xi )=\varphi (1-\xi )\) and by recursion yields
for \(k\in \mathbb {Z}\). Using the initial condition (A.4) and Eq. (A.8), we find that
We can now expand that
Then, the third derivative that appears in the convolution can be written as
where \(\delta \) is the Dirac-delta distribution. Due to the convolution integral (A.5), we are only taking into account past values of \(\lambda \), which yields
where
We can also transform the last remaining integral into a differential equation by introducing
which then gives the initial value problem
Note that the harmonic term in (A.9) is canceled by the homogeneous solution of (A.11), and hence
and the initial condition of (A.11) vanishes
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\le t<\min \left( 2\xi ^{\star },2-2\xi ^{\star }\right) \), it is sufficient to use
which then yields
which is the result we sought.
Proof of Theorem 4.8
The following proof of Theorem 4.8 investigates whether a trajectory approaching \(\varSigma _{0}^{\pm }\) can be continued uniquely after reaching \(\varSigma _{0}^{\pm }\) 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 \(\varSigma _{0}^{\pm }\) 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 \(\varSigma _{0}^{\pm }\). We denote the solution of (4.16) by \((\varvec{\eta }(t),\lambda ^{\star })\), and the solution of (4.20) by \((\varvec{\sigma }(t),\lambda (t))\) either of which can form \(\mathcal {T}\).
First, we prove case 1. The speed of solutions relative to \(\varSigma _{0}^{\pm }\) on the two sides of \(\varSigma _{0}^{\pm }\) is given by \(\frac{\mathrm {d}}{\mathrm {d}t}h_{0}(\varvec{\eta }(t),\lambda ^{\star })\) and \(\dot{\lambda }\), respectively. Trajectories cross \(\varSigma _{0}^{\pm }\) if the signs of these two quantities are the same. We calculate that
and rearrange Eq. (4.18) into
At \(t=0\), the right sides of (B.1) and (B.2) are equal. Using assumption (4.21), we infer that \(\frac{\mathrm {d}}{\mathrm {d}t}h_{0}(\varvec{\eta }(t),\lambda ^{\star })\) and \(\dot{\lambda }\) have the same sign at \(t=0\), and hence \(\mathcal {T}\) has a unique continuation transversely through \(\varSigma _{0}^{\pm }\).
We now show case 2 of the theorem. Assume that \(\mathcal {T}\) is tangent to \(\varSigma _{0}^{\pm }\) to order \(\ell -1\). This means that either
if \(\mathcal {T}\) is a trajectory of (4.20) or
if \(\mathcal {T}\) is a trajectory of (4.16).
We first assume that \(\mathcal {T}\) is a trajectory of (4.20) and (B.3) holds. Let us consider
which is the constraint that keeps the trajectory on \(\varSigma \) [cf. (4.18)]. Any derivative of order j with respect \(\vartheta \) in formula (B.5) includes a \(\frac{\mathrm {d}^{j}}{\mathrm {d}t^{j}}\lambda (t)\) factor. Using (B.3), we can simplify (B.5) to
and
We now show that
The left side of (B.8) is an algebraic expression of \(D_{1}^{j}h_{0}\left( \varvec{y}^{\star },\lambda ^{\star }\right) \) and \(D^{j}\varvec{\sigma }(0)\), \(0<j\le k\). Also, \(D^{j}\varvec{\sigma }(0)\) can be written as
which depends on derivatives of \(\lambda \) up to order \(j-1\) that are all zero according to (B.3). Therefore, none of \(D^{j}\varvec{\sigma }(0)\), \(0<j\le k\) depends on the nonzero \(\frac{\mathrm {d}^{k}}{\mathrm {d}t^{k}}\lambda (0)\). The same holds true for the right side of (B.8) where \(\lambda \) is assumed to be constant, which proves (B.8). Substituting (B.8) into (B.6) implies that (B.4) follows from (B.3). Further, substituting (B.8) into (B.7) we find that
So far, we have shown that if vector field (4.20) is tangent of order \(\ell -1\) to \(\varSigma _{0}^{\pm }\) at \(\left( \varvec{y}^{\star },\lambda ^{\star }\right) \) 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 \(\ell =1\) does not require assumption (B.3). It directly follows from Eq. (B.9) that
Now knowing that \(\frac{\mathrm {d}}{\mathrm {d}t}\lambda (t)=0\) we can apply (B.9) for \(\ell =2\) and conclude that
Repeating this procedure a sufficient number of times shows that (B.4) implies (B.3).
In summary, (B.4) holds if and only if (B.3) holds and we have the equality
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 \(\varSigma _{0}^{\pm }\). If \(\ell \) is odd, trajectory \(\mathcal {T}\) passes through \(\varSigma _{0}^{\pm }\) at the tangency. If \(\ell \) is even, the trajectory continues on the same side of \(\varSigma _{0}^{\pm }\) and there is no joining trajectory from the other side of \(\varSigma _{0}^{\pm }\). Conditions (4.21) and (4.22) also imply that either case 1 or 2 holds for points on \(\varSigma _{0}^{\pm }\) in a sufficiently small open neighborhood of \(\left( \varvec{y}^{\star },\lambda ^{\star }\right) \). This excludes cases where trajectories are forced onto \(\varSigma _{0}^{\pm }\) for a nonzero length of time and implies that there cannot be trajectories joining \(\left( \varvec{y}^{\star },\lambda ^{\star }\right) \) from within \(\varSigma _{0}^{\pm }\). Therefore, there is a unique continuation of \(\mathcal {T}\) for \(t>0\) (or \(t<0\)) sufficiently small. \(\square \)
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 \(\lambda \) from the algebraic constraint
by differentiation, which follows from Lemma 4.10. At this step, we claim continuity of \(\dot{\lambda }\). By investigating the resulting equation, we move to step two and establish that \(\dot{\lambda }\) is indeed continuous in \(\varSigma \). In the final step, we show that trajectories transverse to \(\varSigma ^{\pm }\) cross \(\varSigma ^{\pm }\) when \(d^{\pm }(\varvec{y},\varvec{z},\lambda )>0\), which concludes the proof.
For convenience, we copy here Lemma 4.10.
Lemma C.1
(Lemma 4.10) Assume (A3) and that \(\lambda \) is continuously differentiable and \(\varvec{y}\), \(\varvec{z}\) satisfy the differential equations
on the interval \(t\in [s,s+\epsilon )\), \(\epsilon >0\) with an initial condition \(\varvec{y}(s)\in G\), \(\varvec{z}(s)\in \varvec{\mathcal {D}}\). Then, the right-side derivative of h as a function of time is calculated as
where
Proof
Consider Eq. (C.1) with solution \(\varvec{y}\) and \(\varvec{z}\) on the interval \([s,s+\epsilon )\). We start with the expression
and show that it can be transformed into (C.2). Let us define \(\overline{\varvec{x}}=\varvec{W}(\varvec{y},\lambda )+\varvec{z}\vert _{t=s}\). The only unresolved term
is obtained by taking the derivative of
while assuming constant \(\overline{\varvec{x}}\). The integral in (C.5) is well defined in the Riemann sense, because of Assumption (A5) and because \(\dot{\lambda }\) is continuous. To simplify notation, we define
so that
Assumptions (A5) and (A3) imply that \(\eta (t,s)\) is continuous for \(t>s\), but also allow a discontinuity at \(t=s\), which needs to be taken into account. Differentiating the convolution in (C.7) yields
The integral on the right side of (C.8) is approximated by a Riemann sum
where \(\delta =\left( t-s\right) /m\), \(m>1\) is an integer and \(c_{k}\in (0,\delta )\). We use finite differences to approximate the derivative
The scheme of the finite difference is such that for \(k=m-1\) the second and first argument of \(\eta \) are equal in one of the terms, i.e., \(t-\delta +c_{k}=s+k\delta +c_{k}\), which takes into account the discontinuity of \(\eta \). If this discontinuity is not taken into account, the integral in the limit \(t\downarrow s\) would vanish. In summary, we have the integral
Taking the limit \(t\downarrow s\), is the same as \(\delta \downarrow 0\), and hence we calculate that
which implies that
Using the definition of \(\eta \) and formula (C.8) gives
which in turn is put into (C.4)
This proves Lemma 4.10. \(\square \)
We now prove Theorem 4.13.
Proof
Proof of Theorem 4.13. We want to derive \(\dot{\lambda }\) from the constraint \(h(\overline{\varvec{x}})=0\), where \(\overline{\varvec{x}}=\varvec{W}(\varvec{y},\lambda )+\varvec{z}\). We take the time derivative \(\frac{\mathrm {d}}{\mathrm {d}t^{+}}h(\varvec{W}(\varvec{y},\lambda )+\varvec{z})=0\) and use the expression (C.2) to solve for
This concludes the first step of the proof. We now demonstrate that \(\dot{\lambda }\) is indeed continuous. We only need to recall that the formal solution for any history of \(\lambda \) is
which allows us to write that
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.
Finally, we need to show that if \({Dh}(\overline{\varvec{x}})\cdot \left( D_{1}\varvec{W}(\varvec{y},\lambda )\varvec{f}(\varvec{y},\lambda )+\varvec{A}_{1}(\varvec{y},\lambda )\varvec{z}(t)\right) \ne 0\), trajectories cross \(\varSigma ^{\pm }\). This renders trajectories unique as they pass through \(\varSigma ^{\pm }\). Indeed, if \(\frac{\mathrm {d}}{\mathrm {d}t^{+}}h\) outside of \(\varSigma \) and \(\dot{\lambda }\) inside of \(\varSigma \), but right on the boundary \(\varSigma ^{\pm }\), have the same sign, trajectories cross \(\varSigma ^{\pm }\). We only need to evaluate (4.30) with \(\dot{\lambda }=0\) and compare that to formula (C.9). The two values differ by the factor \(d^{\pm }(\varvec{y},\varvec{z},\lambda )\), and hence if \(d^{\pm }(\varvec{y},\varvec{z},\lambda )>0\) and \({Dh}(\overline{\varvec{x}})\cdot \left( D_{1}\varvec{W}(\varvec{y},\lambda )\varvec{f}(\varvec{y},\lambda )+\varvec{A}_{1}(\varvec{y},\lambda )\varvec{z}(t)\right) \ne 0\), trajectories cross \(\varSigma ^{\pm }\) and solutions are unique. \(\square \)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Szalai, R. Model Reduction of Non-densely Defined Piecewise-Smooth Systems in Banach Spaces. J Nonlinear Sci 29, 897–960 (2019). https://doi.org/10.1007/s00332-018-9508-4
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00332-018-9508-4