Stability of a nonideally excited Duffing oscillator

This paper investigates the dynamics of a Duffing oscillator excited by an unbalanced motor. The interaction between motor and vibrating system is considered as nonideal, which means that the excitation provided by the motor can be influenced by the vibrating response, as is the case in general for real systems. This constitutes an important difference with respect to the classical (ideally excited) Duffing oscillator, where the amplitude and frequency of the external forcing are assumed to be known a priori. Starting from pre-resonant initial conditions, we investigate the phenomena of passage through resonance (the system evolves towards a post-resonant state after some transient near-resonant oscillations) and resonant capture (the system gets locked into a near-resonant stationary oscillation). The stability of stationary solutions is analytically studied in detail through averaging procedures, and the obtained results are confirmed by numerical simulations.


Introduction
The analysis of forced vibrations in mechanical systems and structures is generally based on the hypothesis of ideal excitation, which means that the forcing amplitude and frequency are assumed to be known a priori. However, the ideality assumption fails to explain and predict certain kinds of behaviour found in real systems [25,43]. In these situations, the excitation needs to be considered as nonideal, meaning that there is a reciprocal interaction between excitation and vibrating response.
The study of nonideally excited oscillations started with the pioneering work of Sommerfeld in 1904 [43]. The experimental setup consisted in an unbalanced electric motor mounted on an elastically supported table. Sommerfeld measured the rotor speed and the oscillation amplitude, while gradually increasing the input power in order to make the rotor speed pass through the resonance frequency of the vibrating system [30]. As represented in Fig. 1, it was discovered that, with a steady rate of input power increment, the increase in rotor speed was notably slowed down when the system reached the vicinity of the resonance region, while the oscillation amplitude grew considerably. This smooth evolution of the system state was found to be interrupted by an abrupt jump of the rotor speed towards a post-resonant value, together with a sudden decrease in the amplitude of the oscillations. A less pronounced jumped phenomenon occurred when passing through resonance with decreasing input power (see Fig. 1). The Sommerfeld effect is closely related to the more general notion of resonant capture [32,35,38]. This phenomenon occurs when a nonlinear system becomes locked in the vicinity of a particular surface of the phase space, known as a resonance manifold, which usually corresponds to a critical value of one of the frequencies present in the system. The alternative, known as passage through resonance, occurs when the system undergoes a transient motion near the resonance manifold and, after some time, it sets out towards the nonresonant region of the phase space.

Actual Evolution Expected Evolution Natural Frequency of the System
In 1992, Rand et al. [22,35] reported an interesting real case of nonideal excitation. They investigated a specific type of failure found in dual-spin spacecrafts when first placed in orbit, known as precession phase lock. Dual-spin spacecrafts consist of two bodies (platform and rotor) connected by a bearing assembly that allows for relative rotation. An internal motor is used to generate a relative spin motion between platform and rotor. When precession phase lock occurs, the energy provided by the motor is found to be channelled into motions other than the intended rotation. Rand et al. concluded that the phenomenon of precession phase lock is mathematically analogous to the resonant capture of an unbalanced motor attached to an elastic support and driven by a constant torque. Quinn et al. investigated approximate analytical methods to distinguish initial conditions leading to capture or passage through resonance [32]. Quinn also extended the analysis to the case of an unbalanced motor supported by orthogonal, linearly elastic supports, constrained to have a planar motion [33].
Fidlin [14] analysed the near-resonant oscillations produced by an unbalanced motor using two different approaches: a standard second-order averaging and a hierarchic averaging procedure, based on a coordinate transformation suggested by Pechenev [31]. Sanders [38] investigated the same system as Fidlin through first-order averaging, looking for conditions leading to resonant capture or passage through resonance.
The literature mentioned up to this point is particularly connected to the present work, due to the similarity in the considered system and assumptions. (The relation between these references and the results of our investigation will be highlighted in Sect. 6.) Some additional works of relevance are briefly described in the next lines. Cveticanin et al. [7,8] used averaging techniques to investigate different configurations of nonideally excited systems, including cases of variable mass. The behaviour of elastically supported unbalanced motors with fractional damping was numerically analysed in [2,47]. Bharti et al. [3] addressed the appearance of the Sommerfeld effect in the torsional vibrations of a double-Cardan joint driveline. Drozdetskaya and Fidlin used a particular type of averaging method, namely the averaging technique for partially strongly damped systems, to study the dynamics of an unbalanced motor with a passive self-balancing system [11].
As an alternative to averaging, many researchers have exploited the method of Direct Separation of Motions [4] for the analysis of nonideally excited oscillators [9,21,24,39,41,42,49,50]. Of course, together with the analytical and numerical research, experimental studies are critical to understand nonideally excited oscillations. In this regard, it is worth mentioning some recent works from Varanis et al. [45,46] where different setups-a cantilever beam, a portal frame and a 3-DOF shear-building structure-are used for experimental validation. Also relevant is the experimental demonstration given by Kossoski et al. [26] of the use of shape memory alloys to mitigate the Sommerfeld effect.
From the above paragraphs, it is clear that nonideally excited oscillations have attracted wide attention from researchers in the last decades. On the other hand, the Duffing oscillator is surely one of the most studied systems within the field of nonlinear dynamics [44]. The fact that such an apparently simple equation-a linear oscillator with an extra cubic term in the restoring force-can give rise to jumps, chaotic behaviour, bursting oscillations, etc., makes the Duffing oscillator a very useful system to understand the nature of these complex dynamic phenomena [27,34]. In addition, the Duffing equation has great practical importance because it is able to capture the nonlinear behaviour of many physical systems: pendulums [37], cables [48], elastic beams undergoing large displacements [19,29,36], nonlinear electrical circuits [27], etc. The present work is intended to analyse how the Duffing oscillator behaves when, instead of being subjected to the usual harmonic forcing with known frequency and amplitude, it is excited by a nonideal energy source.
The organization of this article is as follows. Section 2 defines the problem under consideration and the underlying assumptions. The analytical treatment of the equations through averaging techniques is presented in Sect. 3. Section 4 shows that an extension of the time scale of validity of the averaged system is required to prove stability of near-resonant solutions. Section 5 presents some numerical results to confirm and illustrate the analytical developments. A discussion and interpretation of the obtained analytical and numerical results is offered in Sect. 6, which also highlights the relation between the present paper and some relevant literature. Section 7 summarizes the work and presents its conclusions.

Problem statement and assumptions
We investigate the dynamics of an unbalanced electric motor attached to a fixed frame through a nonlinear spring-whose force has linear and cubic components-and a linear damper, as represented in Fig. 2. Variable x stands for the linear motion, φ is the angle of the rotor, m 1 is the unbalanced mass with eccentricity r , m 0 is the rest of the vibrating mass, I 0 is the rotor inertia (without including the unbalance), b is the viscous damping coefficient, and k and λ are, respectively, the linear and cubic coefficients of the spring.
The equations of motion for this 2DOF system are [13] mẍ + bẋ + kx + λx 3 = m 1 r (φ 2 cos φ +φ sin φ) where m = m 0 + m 1 , I = I 0 + m 1 r 2 and a dot represents differentiation with respect to time t. Gravity can Fig. 2 Representation of the system under analysis: an unbalanced motor with known torque-speed curve connected to a fixed frame through a Duffing-like spring and a linear damper be shown to have no relevance for the purposes of this work [10], since the corresponding term in the equations becomes zero after averaging. For this reason, no gravity term has been included in Eq. (1). Function L m (φ) represents the motor characteristic. We assume this driving torque to be linearly related to the rotor speed: with C > 0, D < 0 and ω n = √ k/m. It should be pointed out that writing the motor torque as a function of the rotor speed, and not of some internal electrical variables (voltages, currents or magnetic fluxes), corresponds to a quasi-steady state hypothesis. This means that the dynamics of the electrical variables are assumed to be much faster than the dynamics of the mechanical variables. For a thorough discussion on the validity and limitations of this approach, see [6].
It is useful to define the following dimensionless parameters Let us also define dimensionless versions of displacement x and time t: Now system (1) can be rewritten in dimensionless form: where a dot now represents differentiation with respect to dimensionless time τ .
In order to investigate system (5) using perturbation techniques, some assumptions on the order of magnitude of the parameters are needed. We assume the damping, the unbalance, the spring nonlinearity and the driving torque at resonance to be small. This is expressed by making the corresponding parameters proportional to a sufficiently small, positive, dimensionless parameter ε: where parameters with subscript 0 are ε-independent.
It can be shown that the system dynamics and the required perturbation approach are very different depending on the order of magnitude of the motor characteristic [16,17]. We distinguish between the cases of large slope (d = O # (1)) and small slope (d = O # (ε))-symbol O # represents a sharp estimate, see Definition 1.4.4 in [38] or Definition 1.1.4 in [12]. While the case of large slope was investigated by the authors in [16][17][18], the case of small slope is the focus of the present work: A global comparison between the cases of large and small slope can be found in Chapter 7.2 of [15].

Perturbation approach
By introducing Eqs. (6) and (7) into Eq. (5), the system can be written as where subscript 0 has been dropped for convenience. It is useful to transform system (8) according to change of variables Let us also define a new variable for the rotor speed: ≡φ.
It can be shown that system (8), written in terms of the new variables, takes the form (see [17] for details) where A direct inspection of system (11) Then, by expanding the products of sines and cosines in (11), the system can be written aṡ System (14) is in the appropriate form for averaging over the 3 fast angles {ϕ 1 , ϕ 2 , ϕ 3 }. However, this averaging procedure fails in the vicinity of = −1 and = 1-note that ϕ 1 or ϕ 3 is no longer a fast variable when = 1 or = −1, respectively. Since we assume the rotor speed to be always positive, the only relevant singularity occurs at = 1. This can also be expressed by stating that the system has a resonance manifold at = 1. Thus, two different regions of the phase space need to be distinguished: far away from the resonance manifold (outer region), system (14) can be averaged over the 3 fast angles, while in the vicinity of the resonance manifold (inner region) a different averaging approach is required.

Outer region
If the rotor speed is away from 1, we can average system (14) over the three fast angles {ϕ 1 , ϕ 2 , ϕ 3 }. Note that system (14) is T . Its corresponding first-order averaged system can be obtained asẏ This general procedure is described in detail in Section 7.7 of [38]. The resulting averaged system iṡ where System (16) approximates system (14) with O(ε)precision on a time scale 1/ε, as is known from averaging theory [38]. Note that, in the outer region of the phase space, the vibration amplitude and the rotor speed evolve independently of each other. A direct analysis of system (16) yields the conclusion that it has one only fixed point, given by which is globally asymptotically stable as long as d < 0. Note that, according to the assumption that c > 0 and d < 0, the equilibrium point given in Eq. (18) corresponds to a post-resonant regime ( > 1), as represented in Fig. 3.
Two different scenarios can be considered: • If (0) > 1, system (16) is exponentially attracted towards equilibrium (18) without approaching the resonance manifold. It can also be proved that, in this case, the outer averaged system is valid for all τ ∈ [0, ∞) (see Sect. 5.5 in [38]). • If (0) < 1, system (16) is also exponentially attracted towards fixed point (18). However, in its way towards the equilibrium, the considered trajectory will necessarily reach the neighbourhood of the resonance manifold, making system (16) no longer valid. There are, in principle, two options: -The system remains close to the resonance manifold for all subsequent time (resonant capture or locking into resonance). -The system stays near the resonance manifold for some finite time, after which it continues its evolution towards fixed point (18) (passage through resonance).

Inner region
In order to study the system behaviour in the neighbourhood of the resonance manifold, the rotor speed is expanded as Introducing Eq. (19) into system (11) yieldṡ with System (20) contains three slow variables {a, β, σ } and a fast angle φ. It is useful in this case to apply a secondorder averaging procedure, which provides more accurate results than the usual first-order averaging. Observe that system (20) is of the formẋ = √ Its corresponding second-order averaged system can be obtained asẏ = √ εG 1 ( y)+εG where H is a term that depends on G 1 and G 1 . The details of this general methodology can be found in Section 7.9 of [38]. The obtained approximate system iṡ Note that an overbar has been used to emphasize the difference between the original variables {a, β, σ, φ} (solutions of system (20)) and the averaged variables {ā,β,σ ,φ} (solutions of system (22)), whose relation is given by with L an ε-independent constant. The general procedure that allows obtaining the error estimates (23) is described in Section 7.9 of [38]. Despite the fact that the evolution of {ā,β,σ } is independent ofφ-which is precisely the purpose of averaging-system (22) includesφ as a state variable. The reason is that variableφ(τ ) is required to construct the error estimates in (23). However, in order to investigate the dynamics of the averaged system, it is convenient to rewrite it without the fast angle:ȧ A direct analysis of system (24) reveals that, if c > α/2, there are no fixed points. On the other hand, if c < α/2, the system exhibits two fixed points, given by with where R 0 ≡ 1 − a 2 0 and z = ±1. The two possible values of z correspond to the two different equilibrium points.
It is now illustrative to represent the obtained fixed points on a torque-speed plot. This will provide a clear comparison between the cases of large and small slope of the motor characteristic. To this end, consider the conditionσ = 0 applied to Eq. (24), which yields This can clearly be interpreted as a torque balance condition, since the terms on the r.h.s. of the third equation in (24) represent the different torques acting on the rotor. On the other hand, conditionȧ = 0 yields which allows writing Eq. (27) as If we now define functions T m (dimensionless motor torque) and T v (dimensionless vibration torque) as the fixed points are defined by Torque-speed curves in the resonance region for ρ > 0 (for ρ < 0 the curve is inclined to the left). The fixed points of system (24) correspond to the intersections between T m and T v In addition, conditionβ = 0 yields The torque-speed curves can be represented as follows.
First, graph T m versus σ 1 (in this case, since T m ≡ c, a constant function is obtained). Then, plot the parametric curve given by The fact that a is strictly positive comes from its definition as the radius of a polar coordinate transformation-see Eq. (9)-while condition (28) forbids a 0 from being greater than 1. This procedure gives rise to a torque-speed representation such as the one shown in Fig. 4. The reader interested in the relation between the cases of large and small slope of the motor characteristic may find illustrative to compare Fig. 4 with Fig. 4 in [17]. Figure 4 shows the existence of two fixed points, as long as c < α/2. The stability of these equilibria is now investigated. To this end, the Jacobian matrix of system (24) needs to be obtained and evaluated at the equilibrium point of interest: with Now we compute the eigenvalues of matrix J eq . After some algebra, the following result is obtained for the fixed point at the right branch of the torque curve (z = 1): which clearly corresponds to an unstable equilibrium of the saddle type. For the left branch (z = −1) we obtain where i represents the imaginary unit. According to Eq. (37), the fixed point on the left branch is stable as long as d < 0, since all eigenvalues have negative real parts. Note a significant feature of Eqs. (36) and (37): the stability of both fixed points is independent of the cubic nonlinearity. If condition d < 0 is met, then the left equilibrium is stable and the right one is unstable, regardless of the value of parameter ρ. This is an important difference with respect to the classical Duffing oscillator, for which it is well known that the stability of near-resonant solutions is highly dependent on the cubic nonlinearity [44], as represented in Fig. 5.

Analysis of the time scales
This section focuses on a very relevant detail about the stability results obtained in Sect. 3. Note that not only give the eigenvalues information about stability, but also about the time needed by trajectories to be attracted to or repelled by a fixed point.
For z = 1 (right branch of the vibration torque curve), we have Re( where λ 2 is the only eigenvalue with a positive real part-see Eq. (36). Thus, repulsion of trajectories from the fixed point occurs on the time scale 1/ √ ε. Note that this is precisely the time scale on which the inner approximation is valid, as specified in Eq. (23). Consider now a trajectory of the averaged system (24) starting sufficiently close to the considered fixed point. If this trajectory is tracked for the time length on which the averaged system is valid (O(1/ √ ε)), it is found to be repelled out of the equilibrium. Therefore, it can be stated that this fixed point, which is unstable in the averaged system (24), is also unstable in the exact system (20).
, with all three eigenvalues having negative real parts as long as d < 0. Hence, the attraction of trajectories towards the fixed point takes place on the time scale 1/ε. On the other hand, from Eq. (23), the averaged system is known to be valid on the time scale 1/ √ ε. Since 0 < ε 1, we can write With a similar reasoning to the one used for z = 1, we can consider an arbitrary trajectory of the averaged system (24) starting in the neighbourhood of the fixed point. If we follow the trajectory for the time length on which the averaged system is valid (O(1/ √ ε)), it is found to be neither attracted nor repelled by the equilibrium, since the time scale is not long enough. On the other hand, if a longer time scale is considered, solutions of the averaged system may not be good approximations to those of the exact system anymore.
From the above considerations, it is clear that a question remains open at this point of the analysis. The fixed point on the left branch of the vibration torque curve has been shown to be asymptotically stable in the averaged system (24) if d < 0. However, is it also stable in the exact system (20)? We answer this question based upon the following result. Result 4.1 Assuming d < 0, and for initial conditions sufficiently close to the fixed point corresponding to z = −1, solutions of the averaged system (24) and solutions of the exact system (20) satisfy Note by comparing Eq. (39) to Eq. (23) that the time of validity of the averaged approximation has been extended (in the vicinity of the stable fixed point), paying the price of losing an order of magnitude in √ ε in the accuracy of the estimate.
The proof of this result is based on the attraction properties of system (24) in the vicinity of an asymptotically stable fixed point. We note that the use of attraction properties to extend the time of validity of asymptotic approximations is standard in averaging theory (see Chapter 5 in [38]). However, the problem under study presents a particular difficulty in this regard, due to the weakness of the attraction-as has been stressed before in this section, the real part of the eigenvalues in Eq. (37) . This makes necessary a modified version of the general theorem presented in Section 5.5 of [38]. In order to avoid making the main body of the paper too long and cumbersome, the details of the proof can be found in Appendix I. In summary, Result 4.1 allows stating that, if the fixed point corresponding to z = −1 is asymptotically stable in the averaged system (24), the stability is preserved in the exact system (20).

Numerical simulations
This section presents some numerical results intended to confirm and illustrate the analytical developments of previous sections.
The simulations are carried out as follows. After assigning specific values to the system parameters, a set of initial conditions for the original system (8) is chosen with (0) < 1, i.e. in the pre-resonant region of the phase space. Then, differential equations (8) are numerically solved for a time interval [0, τ f ] that is long enough to ascertain whether the system is captured or passes through resonance.
Suppose the system passes through resonance. Looking at the obtained numerical solution, two par-ticular instants τ 1 and τ 2 are defined, at which the system enters and leaves the resonance region, respectively. Although the choice of these two values is somewhat arbitrary, they give an approximation to the limits between the inner and outer solutions. Then, the outer averaged system (16) is solved for τ ∈ [0, τ 1 ] and τ ∈ [τ 2 , τ f ], with the initial conditions obtained as the solution of the original system evaluated at and τ = 0 and τ = τ 2 , respectively. The inner averaged system (22) is solved for τ ∈ [τ 1 , τ 2 ], with the initial conditions corresponding to the solution of the original system particularized at τ = τ 1 . Finally, the solutions of the original, outer and inner systems are represented together, in order to confirm the accuracy of the averaged solutions.
The selected set of dimensionless parameters is with ε = 0.001. This set of parameters gives rise to the torque-speed curves in the resonance region depicted in Fig. 6. Note that condition c < α/2 is fulfilled, which implies that there exist two fixed points in the inner region of the phase space (marked as A and B in Fig. 6).
As discussed in Sects. 3 and 4, the equilibrium on left branch of the vibration torque curve (point A) is  Fig. 6 Torque-speed curves near resonance for parameter values (40) stable-since we have d < 0-while the one on the right branch (point B) is unstable. These fixed points can be readily computed by introducing parameter values (40) into expressions (26):

Fixed points in the inner region
The third equilibrium, in the outer region of the phase space, can be obtained through Eq. (18):

Fixed point in the outer region
This scenario is particularly useful to highlight the difference between the classical (ideally excited) Duffing oscillator and its nonideal counterpart. In order to see this in a clear manner, it is convenient to resort to the amplitude-frequency plot, as usually done when analysing the primary resonance of the classical Duffing oscillator [44]. This can be easily done by representing the last of relations (26) for the parameter values given in Eq. (40), as shown in Fig. 7. Note that the two fixed points A and B are also represented on this curve. There is a straightforward relation between the torquespeed plot in Fig. 6 and the amplitude-frequency (or amplitude-speed) plot of Fig. 7, since the vibration  Fig. 7 Amplitude-frequency relation for the parameters in Eq. (40). This curve is given by the last of relations (26) torque is proportional to the squared amplitude-see Eq. (30). For a classical (ideally excited) Duffing oscillator with an amplitude-frequency curve such as represented in Fig. 7, it is known that point A is unstable and point B is stable [44] (see Fig. 5). However, according to the analytical results of Sects. 3 and 4, the situation is exactly reversed for the nonideally excited Duffing oscillator: point A is stable and point B is unstable, as specified in Eq. (42). In order to confirm this remarkable difference, some simulation results are presented below.
For the first simulation, consider the following set of initial conditions: The obtained numerical results are represented in Fig. 8, with remarkably good accordance between solutions of the original and averaged systems. In this particular case, the outer averaged solution has been depicted for the whole time range, in order to clearly see that this approximation loses all accuracy once the resonance region is reached. Although Fig. 8 certainly displays a scenario of resonant capture, it is not clear from these graphs whether the system is approaching point A or point B. In order to evaluate this, let us first calculate the value of the rotor speed at these two points. By combining rela-    (19), (25), (26) and (42), we can write where A , B , σ 1A and σ 1B represent the values of and σ 1 at fixed points A and B, respectively. Fig. 9 shows a close-up view of the rotor speed evolution plotted in Fig. 8 on which two horizontal lines corresponding to A and B have been added. It is apparent from this figure that variable evolves towards A , which is in complete agreement with the analytical results of Sects. 3 and 4: point A, which is known to be unstable for the classical Duffing oscillator, is found to be stable for the nonideally excited Duffing oscillator.  Consider now a different set of initial conditions: which yields the numerical results displayed in Fig. 10. Unlike in the previous case, the system is now found to pass through resonance, showing again very good agreement between solutions of the original and averaged equations. Once the system has overcome the resonance region, it evolves towards the outer stable equilibrium given by Eq. (43).
It is also illustrative to investigate a scenario with no fixed points in the inner region of the phase space. With this purpose, consider the set of parameters which is exactly the same as (40) except for a larger driving torque at resonance. The corresponding torquespeed curves near resonance are represented in Fig. 11, exhibiting no intersections between the curves.
Clearly, resonant capture cannot occur in this situation, unless an attractor other than a fixed point existed in the inner region. After conducting numerous simulations with different initial conditions, no numerical evi-  Fig. 11 Torque-speed curves near resonance for parameter values (47) dence of such an attractor has been found. This means that the system appears to pass through resonance for any pre-resonant initial condition, leading towards the outer stable equilibrium given by Note that, for the initial conditions given in Eq. (49), the evolution of the rotor speed is nearly unaffected by resonance, while there is a significant effect on the  vibration amplitude. It is interesting that exactly the opposite case is encountered for the initial conditions given in Eq. (50): whereas the displacement is almost unaltered by resonance, the rotor speed undergoes significant oscillations when the system passes through the resonance manifold. Thus, it is clear that the influence of resonance on the system behaviour strongly depends on the initial conditions. In general, we can state that some transient resonant effects can be expected even when the system has no attractors in the resonance region.

Discussion
From the analytical developments of Sects. 3 and 4, confirmed by the numerical results of Sect. 5, we can draw a general picture on the behaviour of the nonideally excited Duffing oscillator.
The system exhibits two equilibrium points in the resonance region as long as condition c < α/2 is met. Fig. 4 represents both equilibria on a torque-speed plot, where the fixed point on the right branch is unstable and the one on the left branch is stable as long as condition d < 0 holds. The existence of a stable fixed point in the inner region justifies the possibility of resonant capture.
The system reaches the resonance manifold whenever the rotor speed is initially below resonance, as  is readily deduced from the second of Eqs. (16). For some sets of initial conditions, the system trajectory in the phase space will enter the basin of attraction of the near-resonant stable fixed point and, therefore, it will remain close to resonance for all subsequent time (resonant capture). Clearly, there may also be sets of initial conditions whose corresponding trajectories do not reach the basin of attraction of the near-resonant stable fixed point. In these cases, the system may leave the resonance region and evolve towards the fixed point given at Eq. (18), in the outer region of the phase space (passage through resonance).
In principle, it would also be possible for trajectories to be attracted by a different object in the inner region, such as a stable limit cycle or a chaotic attractor. This would represent another kind of resonant capture, not due to the presence of a stable fixed point. However, the numerical simulations carried out have not revealed the existence in the inner region of any attractor other than the investigated fixed point. It should be noted that the results presented in Sect. 5 correspond only to a small portion of the total amount of simulations carried out during this research.
It is worth stressing a particularity of the obtained stability results. The cubic term in the classical (ideally excited) Duffing oscillator is known to have a strong influence in the stability of near-resonant motions [44] (see Fig. 5). However, for the nonideally excited Duffing oscillator, we have shown that, although the cubic term modifies the fixed points near resonance (note the inclination of the vibration torque curve in Fig. 4), it does not affect their stability (observe that parameter ρ does not appear in Eqs. (36) and (37)).
It is also interesting to highlight the relation between the presented investigation and some relevant works in the literature. Sanders et al. considered a system analogous to (8) as an illustrative example in Chapters 7 and 8 of their book [38]. Regarding the outer region of the phase space, they conducted a similar analysis to the one presented in Sect. 3.1, averaging over the three fast angles in Eq. (14) and obtaining the equilibrium given in Eq. (18). However, in the inner region they carried out a first-order averaging, in contrast to the second-order averaging addressed in Sect. 3.2 of this paper. This procedure did not allow them to analyse the stability of the fixed points near resonance, since a firstorder averaging is not accurate enough for this purpose. Moreover, although Sanders et al. included a nonlinear term in the restoring force (see term f (x) in Section 7.5.3 of reference [38]), they could not study the effect of this nonlinearity on the near-resonant oscillations of the system, since this also requires a second-order averaging procedure.
Fidlin devoted Chapter 5 of his book [14] to the study of nonideally excited oscillations, taking system (8) (with ρ = 0) as a relevant example. He focused on the resonance region, arriving at a system analogous to (24) after a second-order averaging. One of the main contributions of the present research to Fidlin's work is the analysis of the effect that a cubic nonlinearity in the restoring force has on the near-resonant behaviour of the oscillator. Moreover, by using a torque-speed plot (see Fig. 4), a clear graphical interpretation has been given to the fixed points of the averaged system near resonance, which in turn provides a useful visual comparison between the cases of large and small slope of the motor characteristic. Finally, it has been proved that a stable fixed point near resonance in the averaged system guarantees stability in the original system as well, as discussed in Sect. 4.
We should also mention the thorough investigations by Rand, Quinn and collaborators on the phenomenon of resonant capture [22,23,32,33,35]. These works aimed at unveiling the dynamics of a linear, undamped oscillator driven by an unbalanced motor with a constant torque, which can be written as system (1-2) with b = λ = D = 0. Interestingly, the purpose was to understand and control precession phase lock, an undesired phenomenon in the motion of spacecrafts, which is governed by the same equations as the unbalancedriven oscillator. It can be shown that the analytical developments of Sect. 3.2 are valid under the assumptions of Rand, Quinn, et al. However, recall that the results of Sect. 3.2 guarantee the stability of one of the near-resonant fixed points as long as the slope of the motor characteristic is strictly negative. Under a constant torque, the procedure presented in this paper cannot determine whether there are any stable nearresonant fixed points in the averaged system, as can easily be observed by putting d = 0 in Eq. (37). Hence, different approaches are required to study resonant capture in the case of constant torque, such as the energy method or the invariant manifold approach described in [32].
Regarding future work on the problem addressed in this article, it would be interesting to analyse the robustness of the presented results with respect to parametric errors (i.e. uncertainties on the values of the system parameters) [1]. In addition, some effort should be devoted to determining or estimating the basins of attraction of the different attractors present in the system [20,40]. This would help to quantify the probability of resonant capture and passage through resonance.

Summary and conclusions
In this article we investigated the dynamics of a 1-DOF oscillator, with Duffing-type restoring force and linear damping, excited by an unbalanced motor. The interaction between motor and oscillator was considered as nonideal, meaning that the rotational motion of the rotor is not assumed to be known a priori, but depends on the oscillator response. The behaviour of such a system is known to depend strongly on the order of magnitude of the slope of the motor torque-speed curve [17]. In this regard, two different scenarios are considered: what we have called the case of large slope and the case of small slope. While the former was addressed by the authors in [16,17], the latter is investigated in this paper.
The system behaviour was analytically studied, using averaging techniques, both within the resonance region and far away from resonance, with particular interest in the phenomena of resonant capture and passage through resonance. Under appropriate conditions, two possible stationary motions were found near resonance. A clear graphical interpretation for these solutions was provided in the form of torque-speed curves, which in turn yields a good visual comparison between the cases of large and small slope of the motor curve. A stability analysis of the averaged system revealed that one of these solutions is always unstable, while the other is stable as long as the slope of the motor curve is negative.
It was also analytically proved that the stability properties of the original system coincide with those of the averaged system. Although this step is straightforward for most systems, it required a detailed analysis in the case under study, due to the very weak attraction of the stable solution in the averaged system. It is worth pointing out that justifying the existence of a stable stationary solution near resonance in the original system provides a solid explanation for the possibility of resonant capture.
Regarding the cubic term in the restoring force, it was found that, although this nonlinearity produces an inclination of the peak in the vibration torque curve, it does not affect the stability of near-resonant solutions. This constitutes a very significant difference with respect to the classical (ideally excited) Duffing oscillator. The obtained analytical results were confirmed by numerical simulations.
Author contributions All authors contributed to the study conception and design, material preparation, analytical development and numerical simulations. All authors discussed the results and contributed to the final manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was supported by Grant FPU12/00537 of the Spanish Ministry of Education, Culture and Sport.

Data availability
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or nonfinancial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix I: proof of result 4.1
This appendix shows how the time of validity of the averaged system (24) can be enlarged for trajectories starting close enough to the fixed point corresponding to z = −1, as long as condition d < 0 holds. It is known that attraction properties can be used to expand the time of validity of asymptotic approximations (see Chapter 5 in [38]). The problem in the case of system (24) is that the attraction is very weak, which makes necessary a modified version of the standard result given in Section 5.5 of [38]. The analysis developed in this appendix is based on an example given in [38]-see page 109, The Case n=2, in the mentioned reference.
Among the possible choices for vector and matrix norms, we will use throughout all three appendices the Euclidean norm for vectors, x = x 2 1 + ... + x 2 n , and its associated operator norm for matrices, A = sup{ Ax : x = 1}.
System (20) is of the form For the sake of generality, instead of presenting the proof using the exact system (20) and the averaged system (24), we will use the exact system (51) and its corresponding averaged system. Regarding system (51), let us assume that all functions of x on the r.h.s. are of class C 1 on an open set U ⊂ R n . We assume further that x remains within a connected, bounded, open set D ⊂ U ⊂ R n , while the angle φ is defined on the circle S 1 . The initial conditions can be denoted as x(0) = x 0 , φ(0) = φ 0 . Note that a superscript in square brackets in Eq. (51) denotes a remainder of an expansion in powers of ε (X [3] depends on ε, while X 1 and X 2 do not). Note also the change in notation between Eqs. (20) and (51): √ ε is now written as ε for simplicity. This can be understood as a simple correspondence between two small parameters ( √ ε 1 = ε 2 ). One last remark on notation: dimensional time and dimensionless time have been denoted as t and τ , respectively, throughout the article. In this appendix the distinction is not relevant and a generic time t will be used.
A standard second-order averaging procedure applied to system (51)-see Section 7.9 in [38]-yields the averaged system with initial conditions y(0) = x 0 − εu 1 (x 0 , φ 0 ), ψ(0) = φ 0 . Function u 1 is defined as with constant ψ 0 chosen in such a way that Functions X 1 and X 2 * are defined as where symbol D represents differentiation with respect to y. The error estimate, i.e. the relation between solutions of the averaged system (52) and those of the exact system (51), is given by with L an ε-independent constant. Relation (57) constitutes a standard result for second-order averaging-see Section 7.9 in [38]. We now turn to consider the fixed points of the averaged system. Assume that the averaged system (52), written without the fast angle aṡ has a fixed point y = y eq . Assume further that this fixed point is asymptotically stable, with all eigenvalues of the Jacobian being distinct and where λ i is the eigenvalue of the Jacobian with the greatest real part. This implies that attraction of trajectories towards the equilibrium takes place on the time scale 1/ε 2 . On the other hand, as specified in Eq. (57), the solutions of the averaged system (52) only approximate those of the exact system (51) on the time scale 1/ε. Thus, we can write It is clear from Eq. (60) that, in order to prove that the fixed point under consideration is stable in the exact system, the time of validity needs to be extended somehow. Note that, up to this point in the appendix, we have simply presented a generalized version of the situation found in Sects. 3.2 and 4. Now we turn to the proof itself of Result 4.1.
In order to extend the time of validity of the averaged system by using attraction arguments, the attraction properties of system (58) are investigated. Consider any pair of solutions y 1 (t) and y 2 (t) of system (58) starting in the Poincaré-Lyapunov domain of the fixed point. According to the Poincaré-Lyapunov theorem [38], the difference between these two solutions decreases exponentially: with μ > 0, C ≥ 1. Without loss of generality, we now assume that Eq. (61) holds with C = 1 (see Appendix II). On the other hand, it is known from Eq. (59) that the contraction coefficient μ can be written as with μ 0 an ε-independent constant. Introducing Eq. (62) and assumption C = 1 into Eq. (61) yields In order to apply a contraction argument, we will need to guarantee that, after a time interval O(1/ε), the norm of the difference between y 1 and y 2 has decreased. By introducing t = L/ε into Eq. (63), we have Let x(t), φ(t) denote the solution of system (51) for initial conditions x(0) = x 0 , φ(0) = φ 0 . Likewise, let y(t) denote the solution of system (58) for initial conditions y(0) = x 0 − εu 1 (x 0 , φ 0 ), with y(0) being within the Poincaré-Lyapunov domain of the fixed point y = y eq . Now, consider the following partition of time: On each segment I m = mL ε , (m+1)L ε we define y (m) (t) as the solution of system (58) with initial condition By definingx(t) = x(t) − εu 1 (x(t), φ(t)), we can rewrite Eq. (67) as For all finite m we have from Eq. (57) with c an ε-independent constant and where y (m) (t) has been continued on I m−1 (the existence properties of the solutions permit this). On the other hand, we can use Eq. (64) to write with 0 < k < 1 (observe the notation · I m ≡ sup t∈I m · ). Combining Eqs. (69) and (70) and using the triangle inequality yields We can now combine Eqs. (71) and (69) and apply again the triangle inequality: Using relation (72) recursively, we obtain Taking the limit for m → ∞ yields the following result for t → ∞: Let us now introduce Eq. (65) into Eq. (74): which can be rewritten as with c 1 an ε-independent constant such that c 1 > 2c μ 0 L . Note by comparing Eqs. (76) and (74) that we have lost an order of magnitude in ε in the estimate due to the fact that constant k is close to 1, which is in turn associated to the very weak attraction of the fixed point of interest. Nevertheless, this attraction proves to be sufficient to guarantee the closeness betweenx(t) and y(t) for all t.
Finally, by recalling definitionx(t) = x(t) − εu 1 (x(t), φ(t)) and taking into account that function u 1 is bounded for any x ∈ D and φ ∈ S 1 , we can write which completes the proof. Note that Eq. (39) is the particularization of Eq. (77) to the problem of the nonideally excited Duffing oscillator.

Appendix II: Transformation to real modal form
This appendix is intended to justify the assumption C = 1 in Eq. (61). Let us first summarize the situation of interest: we investigate a system of the forṁ which is assumed to have an asymptotically stable fixed point y = y eq with all eigenvalues of the Jacobian matrix, For any pair of solutions y 1 (t) and y 2 (t) starting within the Poincaré-Lyapunov domain of the fixed point, the Poincaré-Lyapunov theorem states [38] with μ > 0, C ≥ 1. There is, however, one specific scenario where relation (79) holds with C = 1, as proved in Appendix III: the case in which J eq is in real modal form [5]. A real matrix A with distinct eigenvalues is said to be in real modal form if its nonzero elements coincide with its eigenvalues in the following manner. The real eigenvalues are located on the diagonal of A, while the complex conjugate eigenvalues appear as 2 x 2 blocks, with the real part on the diagonal terms and the imaginary part, with alternate signs, on the off-diagonal terms. The rest of the elements of A are zero. For example, a 3 x 3 matrix A in real modal form with one real eigenvalue λ 1 and two complex conjugate eigenvalues λ 2,3 = α ± iβ would be A = Since the proof presented in Appendix I relies upon the assumption that relation (79) holds with C = 1, it is clear that the obtained result-Eq. (77)-is valid if J eq is in real modal form. This appendix will show that the conclusions of Appendix I are also valid in the general case.
Consider again the situation portrayed in the first paragraph of this appendix, with the Jacobian matrix J eq assumed to not be in real modal form. Let us perform a set of coordinate transformations. First, we define a new vector of variables z as In the vicinity of the fixed point of interest, system (78) can be written aṡ We now define T as a matrix whose columns constitute a real eigenbasis [28] of J eq . This means that T is constructed using the eigenvectors of J eq as follows.
The real eigenvectors are directly stored as columns of T , while, for the complex conjugate eigenvectors, the real and imaginary parts are stored in different columns of T . For example, if J eq had one real eigenvector v 1 and a couple of complex conjugate eigenvectors v 2,3 = u ± iw, we would have T = [v 1 |u|w]. Now we define a new vector of variables q as Clearly, system (81), written in terms of q, takes the forṁ with B eq = T −1 J eq T . Due to the way in which matrix B eq is defined, it has some relevant properties, which are well known from the theory of linear dynamical systems [5]: • B eq has the same eigenvalues as J eq . • B eq is in real modal form.
Then, according to the result presented in Appendix III, we can write with μ > 0 and q 1 (t), q 2 (t) being any pair of solutions of system (83) starting within the Poincaré-Lyapunov domain of the fixed point q = 0. The conclusion is that, using variable q instead of y, the Poincaré-Lyapunov theorem can be written with C = 1.
It has been shown that system (78), subjected to the coordinate transformation y → q, can be written as system (83). Analogously, system (51) can also be transformed according to change of variables {x, φ} → { p, φ}, with Then, using variables p and q instead of x and y, we can legitimately follow the steps from Eq. (66) to Eq. (77), which clearly yields with c 3 an ε-independent constant. Now the original variables can be recovered: x(t) − y(t) = T (p(t) − q(t)) ≤ T p(t) − q(t) .
Hence, we have that with c 4 = T c 3 . It has been shown that the conclusion of Appendix I is valid whether or not matrix J eq is in real modal form. The only difference between both cases is the need for an appropriate coordinate transformation. In Appendix I, the Jacobian matrix was assumed to be in real modal form in order to avoid making the proof unnecessarily cumbersome.

Appendix III: The Poincaré-Lyapunov Theorem with the Jacobian matrix being in real modal form
The Poincaré-Lyapunov Theorem is a classic result on the stability of dynamical systems, whose general formulation and proof can be found in any textbook on the subject, such as [38]. In this appendix, we are concerned with the specific form of the theorem when the Jacobian matrix of the system, at the fixed point of interest, is in real modal form: where x, a ∈ R n ; A is a constant n x n matrix in real modal form, with all eigenvalues having strictly negative real part. The vector field g(x) is continuously differentiable with respect to x in a neighbourhood of x = 0, with Then, there exist constants δ, μ > 0 such that if a <δ x(t) ≤ a e −μt for all t ∈ [0, ∞).
The domain a < δ where the attraction is exponential is called the Poincaré-Lyapunov domain of the system at x = 0.
Proof Consider first the linear system associated with system (89): The solution of Eq. (92) can be written as where Φ(t) is the fundamental matrix, i.e. the solution oḟ Φ = AΦ, where I is the n x n-identity matrix. Recall that matrix A is assumed to be in real modal form, which, by definition, means that A is of the form where λ i (i = 1, 2, . . . , m) are the real eigenvalues of A and λ j, j+1 = α j ± iβ j ( j = 1, 2, . . . , r ) are the pairs of complex conjugate eigenvalues of A.
The solution of Eq. (94) is which, for a matrix A such as shown in Eq. (95), takes the form where T j (t) = e α j t cos(β j t) −e α j t sin(β j t) e α j t sin(β j t) e α j t cos(β j t) , j = 1, 2, . . . , r.
On the other hand, a useful lower bound for Φ(t) 2 can be found by choosing an appropriate vector χ as follows: if H = λ i (i ≤ m), then χ = e i ; if H = α j , then χ = e m+2 j (vectors e k represent the standard basis, with the k-th component being 1 and the rest being zero). With such a choice of vector χ , we have that Since it has been shown that Φ(t) 2 ≤ e 2Ht and Φ(t) 2 ≥ e 2Ht , the plain conclusion is Recall that all eigenvalues are assumed to have negative real parts, which means that H < 0. Hence, we can write It should be stressed that, for a system of the form of Eq. (94) with A being a generic real matrix with all eigenvalues having negative real parts, the standard estimate is y(t) ≤ C a e −μt for some constants μ > 0, C ≥ 1 (see Section 5.2 of [38]). What we have justified is that, when A is in real modal form, this inequality holds with C = 1. The next step would be to show that relation (106), which applies to the linear system (94) for any initial condition a, is also valid for the nonlinear system (89) as long as a < δ, for some δ > 0. Since this second part of the proof is standard, meaning that it does not depend on whether or not matrix A is in real modal form, it will not be displayed here. The interested reader can consult it, for example, in Section 5.2 of [38].