## Abstract

This study is concerned with a new generalized mathematical model for single degree-of-freedom bistable oscillators with harmonic excitation of low-frequency, linear viscous damping and a restoring force that contains a negative linear term and a positive nonlinear term which is a power-form function of the generalized coordinate. Comprehensive numerical mapping of the range of bifurcatory behaviour shows that such non-autonomous systems can experience mixed-mode oscillations, including bursting oscillations (fast flow oscillations around the outer curves of a slow flow), and relaxation oscillations like a classical (autonomous) van der Pol oscillator. After studying the global system dynamics the focus of the investigations is on cubic oscillators of this type. Approximate techniques are presented to quantify their response, i.e. to determine approximations for both the slow and fast flows. In addition, a clear analogy between the behaviour of two archetypical oscillators—the non-autonomous bistable oscillator operating at low frequency and the strongly damped autonomous van der Pol oscillator—is established for the first time.

## 1. Introduction

This study is concerned with the dynamics of bistable single degree-of-freedom oscillators with low-frequency harmonic excitation, linear viscous damping and a restoring force that contains a linear term and a nonlinear term expressed as a general power-form function, whose mathematical model is given by
*a*and
*b*where *x* is a non-dimensional variable that is henceforth referred to as the displacement; overdots denote differentiation with respect to non-dimensional time *t*; *δ* is a non-dimensional damping coefficient; *α* is the power of the nonlinear part of the restoring force that can take any real values higher than unity (to ensure that the restoring force is odd for all values of *α*, including even numbers or non-integer values, noting that the sign and absolute value functions are used in equation (1.1*a*)); *β* is the coefficient of the nonlinear term and is a positive real number; *f*_{0} is the magnitude of the applied harmonic excitation and is also assumed to be positive; *ω* is the angular frequency of the harmonic excitation, and this is considered to be low, i.e. it is of order *ε*.

The dynamics of cubic oscillators modelled by equations (1.1*a*,*b*) with *α*=3 and with *ω* not restricted to low values has been widely examined (e.g. [1] and references therein). However, the general model given by equations (1.1*a*,*b*) has not attracted the same level of attention, despite the fact that it has been associated with, or inspired by, systems appearing in different fields, e.g. dipteran flight [2,3], neurodynamics [4] and mechanical energy harvesting [4]. In order to describe the kinematics of the wing beat in dipteran flight, Thomson & Thompson [2] modelled this as a bistable mechanism which ‘clicks’. They used the fold catastrophe with two control parameters to explain behaviour of the dipteran click mechanism. Brennan *et al*. [3] conducted a dynamic analysis of a simplified model of such a mechanism, with an externally excited bistable Duffing oscillator as the corresponding mathematical model. Making the comparison between this oscillator and a linear resonant one with no click mechanism, by considering the peak kinetic energy of the wing root normalized by power input to the mechanism, they demonstrated a significant advantage in having a click mechanism in the flight motor in the case when the system is highly damped and driven harmonically well below its resonant frequency. Han & Bi [4] presented an analysis of bursting oscillations in damped bistable Duffing oscillators with low-frequency excitation, considering the magnitude of the excitation as a control parameter and determining its threshold for the occurrence of bursting. Cohen *et al.* [5] have recently investigated a bistable vibration-based energy harvesting device by employing numerical slow–fast signal decomposition, showing the merits of this for the purpose of energy harvesting through the conversion of slow into fast oscillations.

The aim of this study is to provide deeper insights into the global behaviour of a wide class of oscillators modelled by equations (1.1*a*,*b*), and to answer the question as to when and how they exhibit mixed-mode oscillations: bursting and relaxation. Section 2 contains basic dynamical characteristics of the system under consideration in terms of fixed points and fold bifurcation points. Section 3 provides a general overview of the types of the system response, both via parametric planes as well as by time-history diagrams and phase portraits. Section 4 is concerned with quantitative techniques for obtaining the solution for bursting oscillations, i.e. the slow and fast flow for cubic oscillators, which have not been obtained so far. Section 5 is concerned with relaxation oscillations and the design of non-autonomous bistable oscillators operating at low frequency with the same jump points and period as autonomous van der Pol oscillators. Section 6 contains conclusion and proposes some directions for further potential use of the results presented here.

## 2. Basic dynamical characteristics

Let us first analyse free ( *f*_{0}=0) undamped (*δ*=0) oscillators (1.1*a*,*b*). The restoring force and the potential energy are, respectively, given by
*a*and
*b*The potential energy is a double-well potential with two local minima (centres C_{−} and C_{+}) at *x*_{C-,+}=±(1/*β*)^{1/α−1} and one local maximum (saddle S) at *x*_{S}=0. Thus, the oscillators are bistable. In the case where there is small damping in the system the saddle S remains a saddle, but the centres C_{−} and C_{+} either become stable foci or stable nodes.

We now focus on the force resulting from the restoring force and the excitation force:
*f* treated as a constant. A geometric presentation of equation (2.2) in the *f*–*x* plane is shown in figure 1 and represents the so-called S-curve (or nullcline). Two vertical tangents are defined by d*F*_{R}/d*x*=0 and exist at points *x*_{c1,2}, corresponding to a fold bifurcation. Together with *f*_{c1,2}=*F*_{R}(*x*_{c1,2}), they give
*a*

and
*b*Two additional points *x*_{c3,4}, which correspond to the solution of *F*_{R}(*f*_{c1,2}),

The values of these characteristic points for cubic oscillators are of importance for the subsequent investigation. They are as follows:
*a*
*b*
*c*It can be seen that the coordinate of the centre of free undamped oscillations C_{+} is between the points *x*_{c2} and *x*_{c3}. These identified characteristic points are of interest for the subsequent numerical and analytical analyses. Firstly, the aim is to determine numerically how the system with low-frequency excitation behaves when the magnitude of the force takes values close to *f*_{c1} and, especially, how the corresponding time-history diagrams look when they are superimposed on the nullcline shown in figure 1.

## 3. Insight into mixed-mode behaviour: numerical investigations

To gain an insight into the regular global behaviour of the oscillators (1.1*a*,*b*), the equations of motion are solved numerically. It is initially assumed that the parameters are *α*=3 and *β*=1. Consequently, equation (2.3*b*) gives the critical value *f*_{0} used in the numerical calculations is chosen to be close to this critical value but slightly higher in order to encompass the influence of damping, i.e. this parameter is fixed at *f*_{0}=0.4. The parameters *ω* and *δ* are varied, and qualitatively different resulting behaviours are captured. Depending on the combination of *ω* and *δ*, three characteristic motions are detected and depicted by dots of different colour (red, blue and green) in the parametric plane *ω*–*δ* (figure 2), dividing this plane into three regions. Note that each of these three regions is characterized by a single attractor of regular behaviour, so there is no need to show basins of attraction in terms of initial conditions.

To illustrate the responses which they represent, three points (Points 1, 2 and 3) are labelled in figure 2 (note that they correspond to the same value of *ω*=0.015). Their responses, in reverse order, are presented in figures 3–5. Figure 3 shows the behaviour typical of Point 3, which dominates in figure 2. In figure 3*a*, the numerical solution of the equations of motion (1.1*a*,*b*) is superimposed onto the nullcline. As can be seen, the motion is associated with only one branch of the nullcline and there is no jump from one branch to another. This motion is periodic about a non-zero value and the corresponding time-history diagram in plotted in figure 3*b*, while the phase trajectory is given in figure 3*c*.

The behaviour of Point 2 which is characteristic for the blue region in figure 2 is presented in figure 4*a*–*c*. Figure 4*a* again shows the numerical solution of the equations of motion superimposed onto the nullcline: the motion consists of stretching along the two branches of the nullcline and jumps at the points *x*_{c1} and *x*_{c2}, approximately. The arrows plotted indicate the system behaviour. Note that when capturing this behaviour numerically, a condition of no evidence of at least one full oscillation around the branches is set. The corresponding time-history diagram in plotted in figure 4*b* and has the form of relaxation oscillations, consisting of slow motions along the outer curves and fast changes in the amplitude (jumps). These relaxation oscillations and the phase trajectory shown in figure 4*c* indicate that this behaviour resembles that of the classical van der Pol oscillator [6–8]. However, the system (1.1*a*,*b*) studied here is non-autonomous, while the classical van der Pol oscillator is autonomous. The links between these two oscillators are established in §5.

The behaviour of Point 1 that is characteristic for the green region in figure 3, is shown in figure 5*a*–*c*. Figure 5*a* contains the numerical solution of the equations of motion superimposed on the nullcline. The motion is such that it stretches along the lower branch of the nullcline then jumps up from the point close to *x*_{c1} to the upper branch, oscillating about it as a damped fast flow, stretching along it, jumping down from the point close to *x*_{c2} to the lower branch, oscillating around it and repeating this behaviour periodically (the arrows plotted give a clear indication of the system behaviour). The corresponding time-history diagram in figure 5*b* shows fast oscillations around the outer curves, slow motion along it, and then sudden changes of the amplitude. The associated phase trajectory is plotted in figure 5*c*.

The parametric chart *ω*−*δ* with the three distinctive regions presented in figure 2 is plotted for one specific value of *f*_{0}. To determine how these regions change if the value of *f*_{0} decreases or increases, two more charts have been produced and are given in figure 6. In the first case, the regions of relaxation oscillations and bursting oscillations shrink, whereas in the second case the opposite is true. It is also notable that relaxation oscillations appear for higher values of *δ* as *ω* decreases. To provide more detailed insights into the change of these regions for other values of *f*_{0}, an animation has been produced (see the electronic supplementary material, animation S1). When *f*_{0} is small (smaller than the critical value), the regions with oscillations around a non-zero value dominate. As *f*_{0} increases the regions of relaxation oscillations and bursting oscillations widen. To illustrate how this parametric plane appears for different values of the power of the nonlinearity and for changeable *f*_{0}, two more animations were created: electronic supplementary material, animation S2 represents the oscillator with a fractional-order nonlinearity defined by *α*=3/2 and the electronic supplementary material, animation S3 is for the quintic oscillator where *α*=5. The general conclusions reached for the system behaviour illustrated in the electronic supplementary material, animation S1 hold for these two cases too, but the sizes and shapes of the regions are different.

## 4. Approximate analytical solution for bursting oscillations

In has been demonstrated in the previous section that the system under consideration can exhibit mixed-mode oscillations, one type of which shown in figure 5*b* includes fast flow oscillations around periodical slow flows. The aim of the following analysis is to provide approximate analytical modelling of these slow and fast flows, which are subsequently labelled *x*_{s} and *X*, respectively. Thus, after representing the overall solution as *x*=*x*_{s}+*X*, substituting it into equations (1.1*a*,*b*) with *α*=3, one can separate this equation into the equation of motion for the slow flow
*τ*=*ωt*, and also the equation of motion for the fast flow

### (a) Approximate analytical solution for the slow flow

To find the approximate solution of equation (4.1) with *α*=3, one can assume a series expansion in the small parameter *ω* for the solution
*ω*:
*x*_{s0} has been obtained from equation (4.4), as the components *x*_{s1} and *x*_{s2} depend on this as well as on its first and second derivative.

Equation (4.4) can be considered as a cubic equation and solved by using the trigonometric method (the corresponding details are given in appendix A). Its solution can be represented in the form:
*p*=−1/*β*, *n*=0.) Several remarks need to be made about the solution given in equation (4.7). First of all, this solution has different forms depending whether it models the outer curve corresponding to *x*_{s}>0 (this will be named the ‘Upper outer curve’ and labelled by ‘U’) or *x*_{s}<0 (this will be named the ‘Lower outer curve’ and labelled by ‘L’). Based on equations (2.5*a*–*c*), this solution for the Upper outer curve can be expressed as
*f*_{0}/*f*_{c1}>1, the argument of the arccos function can be higher than unity for some values of *τ* and, consequently, this function then leads to a complex solution. For example, considering the expression for *x*_{s0U} then from *τ*_{0U}=0 until *Mathematica*, recognizes this issue as it uses the cosine function defined as an entire function). Furthermore, from *τ*_{1U} until *τ*_{1U}+*π* the arccos function gives a solution with a real part and a non-vanishing imaginary part, so that *x*_{s0U} is not purely real anymore. Thus, the solution for *x*_{s0U} given by equations (4.7) and (4.8) is valid until *τ*_{2U}.

This zeroth-order solution is plotted in figure 7 with the time scale *τ* redefined back to the original quantity *t* (it is depicted by the black dotted line). This solution is presented together with the steady-state numerical solution of equations (1.1*a*,*b*), depicted by the solid line, corresponding to the 60th period (note that this moment of time is taken as the origin for the time scale, and that the system parameters are the same as those used in figure 5). It is seen that even this zeroth-order solution agrees well with the Upper outer curve obtained numerically.

Based on equation (4.7), one can express the zeroth-order solution corresponding to the Lower outer curve as follows:
*τ*_{0L}=*π* until *τ*_{2L} the solution for *x*_{s0L} yields a non-vanishing imaginary part.

Now, by differentiating equations (4.7) and (4.9), and then substituting them into equations (4.5) and (4.6), the solutions *x*_{s1}and *x*_{s2} can be obtained. The overall solution (4.3) is not shown here due to its complexity, but is also plotted in figure 7 as magenta dots. Note that the analysis of the expressions for *x*_{s1} and *x*_{s2} in equations (4.5) and (4.6) implies that they become unbounded for *τ*_{2U} and *τ*_{2L} discussed above. As can be seen the difference between the overall and the zeroth-order solution is very small. It should be noted that figure 7 illustrates the accuracy of the approximations for the slow flow time only for the parameters used in figure 5. A few more examples are provided in §4c.

### (b) Approximate analytical solution for the fast flow

In order now to solve the equation of motion for the fast flow (4.2), we treat *x*_{s} as a constant and make it equal to *x*_{c3} for the fast flow around the Upper outer curve and to *x*_{c4} for the fast flow around the Lower outer curve (figures 1 and 5*a*). Hence, equation (4.2) becomes
*a*
*b*
*c*Equation (4.11) contains quadratic and cubic nonlinearity and belongs to the class of the so-called Helmholtz–Duffing equations [1]. Note that the quadratic term is positive for fast flow around the Upper outer curve and negative for fast flow around the Lower outer curve. In addition, the coefficients of the quadratic and cubic terms are not small so this equation is strongly nonlinear. To find the corresponding approximate solutions by using a perturbation approach this equation needs to be transformed into the form conventionally used to represent what one would consider to be a weakly nonlinear equation, and we return to this assumption later in the conclusion. Thus, *X* is rescaled in accordance with *ε*≪1. In addition, to capture appropriately the influence of damping, the damping coefficient is also rescaled by using *δ* =*ε*^{2} *δ*_{1}. Introducing these two substitutions into equation (4.11) and omitting the tildes, one derives the following:
*a*and
*b*where *T*_{N}=*ε*^{N} *t*, *N*=0, 1, 2,… and

Substituting equations (4.14) and (4.15*a*,*b*) into equation (4.13) and collecting terms of the same perturbation order, leads to
*X*_{0} is assumed as
*a* and phase *α* are unknown functions of *T*_{1} and *T*_{2}.

The right-hand side (RHS) of equation (4.17) now becomes
*a*/∂*T*_{1}=0 and ∂*α*/∂*T*_{1}=0, which both imply that the amplitude *a* and phase *α* do not depend on *T*_{1}, but only on *T*_{2}, i.e.
*a*and
*b*The remaining terms in equation (4.20) are considered in order to obtain the following particular integral of equation (4.17):
*a*,*b*) and (4.22), it can be seen that the RHS of equation (4.18) becomes
*a*and
*b*where *a*_{0} and *α*_{0} are constants defined by initial conditions.

The remaining term in equation (4.23) yields the following particular integral for equation (4.18):
*a*,*b*):
*a*and
*b*Owing to the transformation previously introduced the constants *a*_{0} and *α*_{0} should be found from the initial conditions. These are

Numerical solutions for the rescaled equation of the fast flow (4.13) are presented in figure 8 for the Upper curve (*μ*≡*μ*_{U}) and in figure 9 for the Lower curve (*μ*≡*μ*_{L}), together with the approximations defined by equations (4.26) and (4.27*a*,*b*) for the parameters from figure 5 and for *ε*=0.1.

It is seen that the approximate analytical solutions successfully capture the character of the response. Their accuracy is reasonably good given the fact that this is initially large amplitude motion.

### (c) Fast flow superimposed on slow flow

Now, the solutions for the slow and fast flow can be summed but it is important to note that the solutions from figures 8 and 9 need to be rescaled as they actually represent *a*,*b*) and presented in figure 10 for the parameters from figure 5, while the influence of a different value of the damping parameter is shown in figure 11 (*δ*=0.1). These two figures confirm that the approximate analytical solutions capture the response very well, both qualitatively and quantitatively. Note that in both figures the slow flow is approximated by the zeroth-order solution.

## 5. Relaxation oscillations in non-autonomous bistable oscillators and autonomous van der Pol oscillators

An interesting question now concerns the possibility of designing the non-autonomous system (1.1*a*,*b*) to perform relaxation oscillations (figure 4*b*) with the same characteristics (four jumping points and the period) as the autonomous van der Pol oscillator. To that end, one can use the results from §2 and relate them to the characteristics of the van der Pol oscillators, both the classical form [6–8] and generalized forms [9–11].

### (a) Links with the classical van der Pol oscillator

The classical van der Pol oscillator [6–8] is represented by
*ε* in equation (5.1) is large (figure 12). It is known that its jump-up *x*_{ju} and jump-down points *x*_{jd} (figure 12) occur at [7,10]

and we note that the starting points of the outer curves *x*_{u} and *x*_{d} are at
*a*), one obtains *β*=1/3. For this value of *β*, *x*_{c3,4} from equation (2.5*c*) coincides with *x*_{u} and *x*_{d} from equation (5.3). The corresponding value for *f*_{0} is then defined by the value for |*f*_{c1,2}|, equation (2.5*b*), i.e. *f*_{0}=2/3 (as in the previous sections, due to the influence of damping, a slightly higher value is chosen at *f*_{0}=0.72). Finally, this equivalence requires two more parameters to be related to each other, namely the coefficient *ε* in (5.1) and the angular frequency *ω* in equation (1.1*a*,*b*). Firstly, it is known that relaxation oscillations in the autonomous van der Pol oscillator have the period [7] given by
*a*,*b*) has a period defined by
*ω*=0.01, one can then calculate that *ε*=389.293. By using these values both equations (1.1*a*,*b*) and (5.1) can be solved numerically and the results obtained are shown in figure 13 (the dashed lines define the solution for the van der Pol oscillator (5.1) and the solid lines represent the oscillator (1.1*a*,*b*) with *δ*=2.5). Three periods before the 60th period are shown. Note that the time axis is shifted by −260 after solving (5.1) in order to overlap the responses and make the visual comparison easier.

### (b) Links with generalized van der Pol oscillators

This approach can be extended to design generalized bistable oscillators (1*a*–*c*) to have the same period and jump points as generalized van der Pol oscillators [10]
*α*_{1}>1 and *β*_{1}=*α*_{1}−1 (other possible cases existing in [10,11] will be considered later in a separate publication).

To fulfil the requirements, the results from [10] are used: the jump-up *x*_{ju} and jump-down points *x*_{jd} of the oscillator (5.7) are defined by equation (5.2); and the coordinate *x*_{u} is defined implicitly by the following expression:
*a*) with equation (5.2), it follows that *α*=1/*β*. If, for example, the bistable oscillators is quintic, one then has *α*=5 and *β*=1/5. Equation (2.3*b*) gives *f*_{c1}=4/5 (consequently, *f*_{0} is chosen to be slightly higher: *f*_{0}=0.9). Furthermore, equation (2.4) now yields the value *x*_{c3}=1.650629. The next requirement for the equivalence of the oscillators is *x*_{c3}=*x*_{u}, which due to equation (5.8) gives *β*_{1}=4. Given the condition below equation (5.7), one finds that *α*_{1}=5. The period given by equation (5.9) emerges as *T*=0.569668*ε*. By equating this with equation (5.5), it can be shown that *ε*≈11.0295/*ω*. Taking *ω*=0.01 again it can be calculated that *ε*=1102.95. Figure 14 contains the steady-state responses of these two oscillators as in figure 13 (the time shift for the van der Pol oscillator is shifted by −270 and the damping coefficient for equation (1.1*a*,*b*) is *δ*=3.5). The time responses shown confirm that the requirements for equivalent characteristics for the two oscillators have thus been achieved.

## 6. Conclusions

This study has investigated a generalized phenomenological model for a system that exhibits mixed-mode oscillations in the forms of bursting and relaxation motions. The system considered offers the distinct advantages of mathematical simplicity and elegance and exhibits mixed-mode oscillations that are driven by low-frequency external forcing.

The basic dynamical characteristics have been analysed in terms of characteristic values of the magnitude of the force and the displacements. Comprehensive numerical mapping of the range of bifurcatory behaviour in terms of the magnitude of the force, the damping coefficient, the forcing angular frequency and the power of the nonlinearity has been provided. It has been shown that the system under consideration can exhibit oscillations around a non-trivial point, and also bursting and relaxation oscillations. Illustrations of these responses with respect to the nullcline are also provided.

Then, the focus of the investigation was concentrated on cubic oscillators of this type. A quantitative treatment of bursting oscillations in such systems has been carried out, separating these into a periodic upper and lower slow flow and fast flow oscillations around the upper and lower curves. The expressions for the upper and lower slow flows at the zeroth order and second order have been derived and compared with numerical solutions, for which very good agreement has been found. The fast flow oscillations were modelled by a nonlinear damped Helmholtz–Duffing equation whose second-order solutions were derived and compared also with appropriate numerical solutions. It should be noted that the fast flow dynamics typically have the characteristics of an underdamped transient and therefore the normal assumptions needed to employ the method of multiple scales are reasonably adhered to for this part of the solution. The solutions for the slow and fast flows were then summed and it was verified numerically that the resulting overall approximation captures the response very well, in both qualitative and quantitative aspects.

Finally, as the non-autonomous bistable system under consideration can readily exhibit relaxation oscillations, the question was then posed as to how one could design it to have the same response characteristics as a strongly damped autonomous van der Pol oscillator in terms of jump points and the period. As far as the authors are aware, these are the first links to have been established between these two archetypical oscillators to date. In addition because mechanical manifestations of the van der Pol oscillator are rarely found then the results presented in this paper potentially open up new possibilities for the design of analogous mechanical systems, further investigations into such systems and the determination of beneficial applications for the resulting dynamic behaviour.

## Author' contributions

I.K. designed the study, coordinated it, made the concept for obtaining approximate solutions, established the links between two types of oscillators and drafted the manuscript; M.C. participated in the design of the study, supervised all their steps, carried out an approximate technique for the fast flow analyses, and critically revised the draft; M.Z. carried out numerical simulations for behavioural mapping, producing all the animations and the corresponding figures. All authors gave final approval for publication.

## Competing interests

The authors have no competing interests.

## Funding

The involvement of I.K. and M.C. was financially supported by the UK’s Royal Society under the International Exchanges Scheme, grant number IE130406 entitled: ‘Slow-fast nonlinear dynamic phenomena in neural systems’, while M.Z. received financial support from the Ministry of Education, Science and Technological Development of the Republic of Serbia (III 41007).

## Appendix A. Solution of the cubic equation by a trigonometric method

A general cubic equation of the form
*u*^{3} leads to

This leads to *x* in terms of *p* and *q* emerges:

- Received September 8, 2015.
- Accepted November 18, 2015.

- © 2015 The Author(s)