## Abstract

In this paper, we investigate bifurcation phenomena, such as those of the periodic solutions, for the ‘unperturbed’ nonlinear system , with and *α*>1, *β*>0, when we add the two competing terms , with *f*(*t*) a time-periodic analytic ‘forcing’ function and *γ*>0 the dissipative parameter. The resulting differential equation describes approximately an electronic system known as the saturating inductor circuit. For any periodic orbit of the unperturbed system, we provide conditions which give rise to the appearance of subharmonic solutions. Furthermore, we show that other bifurcation phenomena arise as there is a periodic solution with the same period as the forcing function *f*(*t*) which branches off the origin when the perturbation is switched on. We also show that such a solution, which encircles the origin, attracts the entire phase space when the dissipative parameter becomes large enough. We then compute numerically the basins of attraction of the attractive periodic solutions by choosing specific examples of the forcing function *f*(*t*), which are dictated by experience. We provide evidence showing that all the dynamics of the saturating inductor circuit is organized by the persistent subharmonic solutions and by the periodic solution around the origin.

## 1. Introduction

One of the most important second-order linear systems in practical applications is the non-autonomous simple harmonic oscillator, which is described by(1.1)where *L*, *R* and *C* are constants. This is the standard linear, second-order, constant coefficient ODE with non-zero right-hand side. In an electronics context, the well-known inductor–resistor–capacitor circuit is described exactly by this differential equation, provided that the inductor (*L*), resistor (*R*) and capacitor (*C*) are all independent of *x*, which in this case would be the charge on the capacitor. In many cases, this is an excellent approximation. However, a nonlinear capacitor can easily be implemented in practice, for instance, by using a semiconductor diode, and the resulting nonlinear differential equation has been studied in detail, in the case where *x* does not change sign, in Bartuccelli *et al*. (2006, 2007). Another practical but rather less studied nonlinear version of equation (1.1) is one in which *R* and *C* are constant but *L* is a function of . Model functions are of course based on practical experience. Practical inductors consist of a length of wire wound around a core of magnetizable material, real materials being capable of saturation (so that *L* is a function of , as will become clear) and possibly also displaying hysteresis (so that *L* is a function of and ). Hence, the behaviour of real inductors can depart from the ideal, constant *L*. A simple and not very realistic model of saturation was studied in Chua *et al*. (1982), where an even, piecewise constant if , *L*=*L*_{1} otherwise, with *L*_{0}>*L*_{1}, was used. In a bid for a more faithful model of real inductors, a comparison between experiment and simulation was reported in Deane (1994), where both saturation and hysteresis effects were considered. The resulting system was described by a third-order non-autonomous nonlinear ODE, which was solved numerically and compared with experiment.

In this paper, we study a model which lies between the two extremes described earlier: an inductor whose core saturates, this time smoothly, but which does not display hysteresis.

The rest of the paper is organized as follows. In §2, under appropriate assumptions and approximations, we derive and rescale the equation that we study. Such an equation appears as a perturbation of an integrable equation, whose main properties will be reviewed in §3. We then study the periodic solutions of the full unperturbed equation. First, in §4, we investigate the persistence of periodic solutions of the unperturbed system which are ‘in resonance’ with the forcing term (the so-called subharmonic solutions). Then in §5, we show that the system is characterized by the presence of another periodic solution, which has the same period as the forcing and encircles the origin. More precisely, it branches off the origin in the sense that it reduces to the origin when the forcing is removed. Therefore, such a solution is different in a profound way from the subharmonic solutions, because it arises by bifurcation from the stable equilibrium point and not from the unperturbed periodic solutions. The periodic solution encircling the origin has a non-empty basin of attraction which becomes the entire phase space when the dissipative parameter is large enough; this will be proved in §6. In principle, there could be other solutions relevant for the dynamics. However, there is strong numerical evidence that the full dynamics is organized by the periodic solutions that we have studied analytically. In fact, the numerical simulations performed in §7 suggest that the union of the basins of attraction of those periodic solutions fills the entire phase space; cf. Bartuccelli *et al*. (2001, in press *a*) for similar situations. Finally in §8, we draw some conclusions and discuss some open problems and possible lines of future research.

## 2. The saturating inductor circuit ODE

The circuit under consideration is shown in figure 1 and we now derive the ODE that describes it, briefly explaining the assumptions that we make. Details of the electromagnetism and circuit theory required can be found in Kraus (1991). Kirchhoff's voltage law states that the sum of the voltages across the components must equal the applied voltage at all times, and so *g*(*ωt*)=*v*_{L}+*v*_{R}+*v*_{C}. The definition of capacitance states that *v*_{C}=*q*(*t*)/*C*, where *q*(*t*) is the capacitor charge; also by definition, the current ; and Ohm's law allows us to write *v*_{R}=*i*(*t*)*R*. Faraday's law of induction states that *v*_{L}=d*ϕ*/d*t* where *ϕ* is the magnetic flux, to be defined. To proceed further, we need to make assumptions about the geometry of the inductor core. In order that it will display saturation at relatively low currents, we assume that it is toroidal in shape, with cross-sectional area *A* and mean diameter *D*, and that the inductor itself consists of *n* turns of wire around this core. The magnetic flux density is defined by *B*=*μ*_{0}(*H*+*M*), where *μ*_{0}=4*π*×10^{−7} is the permeability of free space, *H* is the applied magnetic field and *M* is the magnetization. In fact, the quantities *B*, *H* and *M* are all vectors, but under the assumption that they are all parallel to each other and normal to the cross-section of the core, we shall see that we need to consider only their magnitudes. The magnetic flux is then defined byNow, Ampère's circuital law relates *H* to *i* by *H*=*ni*/π*D*, and so(2.1)We now need a model for *M*(*H*). In the absence of hysteresis, and with weak coupling between the magnetic domains, *M*(*H*) is usually described by the Langevin functionand where *M*_{s}>0 and *a*>0 are constants which depend on the core material. Now,whereUsing this in equation (2.1) and Kirchhoff's voltage law allows us to write the differential equation that describes the circuit as(2.2)and all the coefficients are well defined even when . We now rescale time by replacing *ωt* with *t*, and *q* by defining , with *λ* a constant to be defined, giving(2.3)where(2.4)

Finally, we can usefully approximate the function *U* with an expression that is qualitatively similar but easier to handle analytically. From the foregoing we have, for all *λ,* that and . The simplest rational approximation to 1+*U*(*λy*) which has the same limits (independent of *λ*) is *L*(*y*)=(*α*+*y*^{2})/(1+*y*^{2}). Since we are free to choose *λ*, we find the value *λ*=*λ*_{0} that minimizesNumerical computation gives the value *λ*_{0}≈1.99932.

Hence, we set *λ*=*λ*_{0} and using the resulting approximation for the function *U*, the differential equation in its final form is(2.5)

In the rest of the paper, we shall be interested in investigating how the dynamics of the ‘unperturbed’ system , will be influenced by the addition of the two competing terms *f*(*t*)−*γy* under the assumption that they represent a ‘small’ perturbation. We shall also be able to investigate cases far from the perturbation regime, at least when the dissipative parameter is large enough with respect to the forcing term.

## 3. The unperturbed system

Consider the ordinary differential equation (we set *γ*=*ϵC* for notational convenience)(3.1)where *g*(*y*)=*β*(1+*y*^{2})/(*α*+*y*^{2}), with *α*>1, *f*(*t*)=*f*(*t*+2*π*) is an analytic periodic function of time, and *ϵ* (small) and *C* are two real parameters. Without loss we can redefine our *ϵ* so as to absorb the positive constant *β*. For *ϵ*=0, our system is integrable by quadratures; in fact it has an integral of motion, which we identify with its ‘energy’, given by(3.2)

Note that the above system is a particular case of the following class of systems possessing an integral of motion: , , the associated integral being given by

The following result is proved in appendix A.

*Consider the system* (3.1) *with ϵ*=0. *Its phase space is filled with periodic orbits around the origin; the period T*_{0} *of these orbits is strictly decreasing as a function of the energy* (3.2), *and satisfies the inequalities* .

For *ϵ*=0, the system (3.1), in the limits of small and large energy, becomes ‘essentially’ the harmonic oscillator with frequencies and , respectively. We observe that for any fixed value of the energy, we have formally an infinite number of unperturbed periodic solutions parametrized by a phase shift *t*_{0}. This phase shift is irrelevant in the autonomous case (*ϵ*=0), but becomes very important when we add the time-dependent perturbation (*ϵ*≠0), because in the extended phase space (*x*, *y*, *t*), each ‘initial’ phase will give rise to a different orbit. Therefore, the unperturbed solution to be continued for *ϵ*≠0 can be written in terms of the initial phase, in the form *x*_{0}(*t*)=*X*_{0}(*t*−*t*_{0}), where *X*_{0}(*t*) is an unperturbed solution with the phase appropriately fixed (for instance with ) in order to impose some constraints on the perturbed system (see below). For practical purposes, it will be more convenient to change the origin of time, by writing the unperturbed solution as *x*_{0}(*t*)=*X*_{0}(*t*) and the periodic forcing function in the form *f*(*t*+*t*_{0}).

*Consider the system* (3.1) *with ϵ*=0. *The energy function* (3.2) *satisfies the identities E*(±*x*, ±*y*)=*E*(*x*, *y*). *The solutions x*_{0}(*t*) *can be chosen to be even in t, and satisfy x*_{0}(*t*)=−*x*_{0}(*t*−*T*_{0}/2), *if T*_{0} *is the corresponding period. Furthermore, they possess only odd components in their Fourier series expansion*.

The equation (3.2) shows that the energy function is even in both variables (*x*, *y*), namely *E*(±*x*, ±*y*)=*E*(*x*, *y*). Therefore, the solutions have reflectional symmetry with respect to both axes (*x*, *y*). Hence, giving initial conditions on the *x*-axis, one obtains solutions *x*_{0}(*t*) that are even in time (similarly if one chooses initial conditions on the *y*-axis then the solutions are odd in time). Using this property, it follows that any solution of period *T*_{0} can be chosen to be even in time, and moreover it must have the property *x*_{0}(*t*)=−*x*_{0}(*t*−*T*_{0}/2).

To show the validity of the part concerning the Fourier spectrum, consider the Fourier series of *x*_{0}(*t*), namely(3.3)where *ω*=2*π*/*T*_{0}. The Fourier coefficients can be written asTherefore, if *ν* is an even number (including zero) , while if *ν* is an odd number . ▪

We now investigate the dynamics of the system when *ϵ*≠0 (and small). In particular, we are interested in studying which periodic orbits of the unperturbed system persist under the combined effect of the periodic perturbation *f*(*t*+*t*_{0})=*f*(*t*+*t*_{0}+2*π*) and the dissipative term ; cf. Melnikov (1963), Hale & Táboas (1978), Chow & Hale (1982), Guckenheimer & Holmes (1990), Gentile *et al*. (2007), Bartuccelli *et al*. (in press *a*) for related results. Note that in order to have a periodic solution with period commensurate with 2*π*, we must impose a so-called resonance condition given by *ω*=2*π*/*T*_{0}=*p*/*q*. If *ω* satisfies the above condition, then the perturbed periodic solution will have the period *T*=*pT*_{0}, during which the periodic forcing function *f*(*t*+*t*_{0}) has ‘rotated’ *q* times. Note that the phase . It is more convenient to consider the system(3.4)defined in the extended phase space (*x*, *y*, *t*). Then for *ϵ*=0, the motion of the variables (*x*, *y*, *t*) is in general quasi-periodic and reduces to a periodic motion whenever *T*_{0} becomes commensurate with 2*π*. Thus, in the extended phase space (*x*, *y*, *t*), the solution runs on an invariant torus, which is uniquely determined by the energy *E*; if *ω*=*p*/*q*, we say that the torus is *resonant*. In general, the non-resonant tori disappear under the effect of the perturbation, unless the system belongs to the class that can be analysed by using the KAM theory. In addition, the resonant tori disappear, but some remnants are left: indeed usually a finite number of periodic orbits, called *subharmonic solutions*, lying on the unperturbed torus, can survive under the effect of the perturbation. If *ω*=*p*/*q*, we shall call *q*/*p* the *order* of the corresponding subharmonic solutions. The subharmonic solutions will be studied in §4.

Besides the subharmonic solutions, there are other periodic solutions which are relevant for the dynamics: indeed there is a periodic solution branching off the origin. This is a feature characteristic of forced systems in the presence of damping (Gentile *et al*. 2005, 2006; Bartuccelli *et al*. 2007, in press *b*). Existence and properties of such a solution will be studied in §§5 and 6. In particular, we shall show that this periodic solution becomes a global attractor if the dissipative parameter is large enough.

## 4. Subharmonic solutions

We consider the system (3.1) with *ϵ*≠0 and small, and we look for subharmonic solutions which are analytic in *ϵ*. First, we formally define power series in *ϵ* of the form(4.1)where, for all *k*∈, the functions *x*^{(k)}(*t*) and *y*^{(k)}(*t*) are periodic with period *T*=*pT*_{0}=2*πq*, with . We shall see that the functions *x*^{(k)}(*t*) and *y*^{(k)}(*t*) can be determined provided the parameters *C*_{k} are chosen to be the appropriate functions of the initial phase *t*_{0}.

If we introduce the decompositions (4.1) into (3.1) and denote with *W*(*t*), the Wronskian matrix for the unperturbed linearized system, we obtain(4.2)where are corrections to the initial conditions, and(4.3)where . Note that by construction *F*^{(k)}(*t*) depends only on the coefficients *y*^{(k′)}(*t*) and with .

We shall denote by (*x*_{0}, *y*_{0}) the solution of the unperturbed system , , with *x*_{0} even (hence *y*_{0} odd): this is possible by lemma 3.1. The Wronskian matrix appearing in (4.2) is a solution of the unperturbed linear system(4.4)so that the matrix *W*(*t*) can be written as(4.5)where *c*_{1} and *c*_{2} are suitable constants such that *W*(0)= (one can easily check that *c*_{1}*c*_{2}=−*α*).

By making explicit the dependence of the unperturbed periodic solutions on *E*, we can write(4.6)so that we can writewhereNote in particular that *h*(*t*) and *b*(*t*) are, respectively, odd and even in *t*, so that is even in time. Hence we havewhere *μ*≠0 because the system is anysochronous and the derivative of the period d*T*_{0}(*E*)/d*E*<0 as shown in lemma 3.1; also the functions are periodic with period *T*_{0}. Note that the lower two components of the Wronskian matrix are the time derivatives of the upper two. Hence we can express the Wronskian matrix as follows(4.7)By introducing (4.7) in (4.2) and using thatwe obtain for all *k*≥1where ; a similar expression can be written for the component . After some manipulation, we obtainFor notational convenience, let us put , and . Then we obtain a periodic solution of period *T*=*pT*_{0}=2*πq* if, to any order *k*∈, one has(4.8)where, given any *T*-periodic function *H* we denote its mean by 〈*H*〉. In addition, recall that and .

The parameters are left undetermined, and we can fix them arbitrarily, as the initial phase *t*_{0} is still a free parameter. For instance, we can set for all *k*∈ or else we can define for *k*∈, with the values to be fixed in the way which turns out to be most convenient (Gentile *et al*. 2007).

Thus, under the restrictions (4.8), we have that for all *k*≥1(4.9)is a periodic function of time. Hence, we must prove that the first condition in (4.8) can be fulfilled at all orders *k*≥1; the second condition just gives a prescription how to fix the parameters .

For *k*=1, the condition on the zero average, in terms of the Melnikov integral (note that gives *M*(*t*_{0}, *C*_{0})=0, and we can choose *C*_{0}=*C*_{0}(*t*_{0}) so that this holds becausefor all *t*_{0}∈[0, 2*πq*/*p*).

For *k*=2, the condition reads (recall that )where the primes denote differentiation with respect to *y*. For higher values of *k*, in order to satisfy the first condition in (4.8), we must choose(4.10)where *Γ*^{(k)} can be proved, by induction, to be a well-defined periodic function depending on the coefficients . Therefore, we conclude that if we set *C*_{0}=*C*_{0}(*t*_{0}) and, for all *k*≥1, we choose , according to the second condition in (4.8) and *C*_{k}=*C*_{k}(*t*_{0}) according to (4.10), we obtain that in the series expansions (4.1) the coefficients *x*^{(k)}(*t*) and *y*^{(k)}(*t*) are well-defined periodic functions of period *T*. All of this of course makes sense if the series expansions (4.1) converge; this can be seen by showing that the coefficients admit bounds of the form |*x*^{(k)}(*t*)|<*δ*^{k} for some positive constant *δ*; hence by taking *ϵ* small enough and *C*=*C*(*ϵ*, *t*_{0}) as a function of *t*_{0} according to (4.10), it follows that the series converges to a periodic function with period *T*=*T*_{0}*p*, analytic in *t* and *ϵ*. The strategy for showing the convergence of the series is essentially the same as that explained in Gallavotti & Gentile (2002) and Gentile *et al*. (2007), and all the details can be found in these references. Hence, the formal power series converges for |*ϵ*|<*ϵ*_{0}=*δ*^{−1}. For fixed *ϵ*∈(−*ϵ*_{0}, *ϵ*_{0}), we shall find the range allowed for *C* by computing the supremum and the infimum, for *t*_{0}∈[0, 2*π*), of the function *t*_{0}→*C*(*ϵ*, *t*_{0}). The bifurcation curves will be defined in terms of the function *C*(*ϵ*, *t*_{0}) as(4.11)These bifurcation curves divide the parameter plane (*ϵ*, *γ*) into two disjoint sets such that only inside the region delimited by the bifurcation curves are there subharmonic solutions; the region of existence contain the real *ϵ*-axis. For more details on this point see Gentile *et al*. (2007).

We are now in a position to summarize the above discussion by stating theorem 4.1 which is an adaptation of theorem 5 of Gentile *et al*. (2007).

*Consider the system* (3.1). *For all p, q relatively prime positive integers such that* , *there exist ϵ*_{0}>0 *and two continuous functions γ*_{1}(*ϵ*) and *γ*_{2}(*ϵ*), *with γ*_{2}(0)=*γ*_{2}(0), *γ*_{1}(*ϵ*)≥0≥*γ*_{2}(*ϵ*) *for ϵ*>0 *and γ*_{1}(*ϵ*)≤0≤*γ*_{2}(*ϵ*) *for ϵ*<0, *such that* (3.1) *has at least one subharmonic solution of order q*/*p for γ*_{2}(*ϵ*)≤*ϵC*≤*γ*_{1}(*ϵ*) *when ϵ*∈(0, *ϵ*_{0}) *and for γ*_{1}(*ϵ*)≤*ϵC*≤*γ*_{2}(*ϵ*) *when ϵ*∈(−*ϵ*_{0}, 0).

The above theorem applies for any analytic 2*π*-periodic forcing function of time; so, provided the integers *p*, *q* are relatively prime, it follows that generically any resonant torus of the system (3.1) possesses subharmonic solutions of order *q*/*p*. Naturally if we choose a particular forcing function with special properties, then in turn this will influence the set of possible subharmonic solutions. For example, if we take the function *f*(*t*)=sin *t*, then for dissipation large enough—that is for *γ*=*O*(*ϵ*)—the only existing subharmonic solutions are those with frequency *ω*=1/(2*k*+1) with *k* any positive integer. This can be seen as follows. The solutions of the unperturbed system have odd components only in their Fourier representation; moreover we know that to first order the Melnikov condition *M*(*t*_{0}, *C*_{0})=0 must be satisfied. Hence we can choose(4.12)which is well defined because the denominator is non-zero; hence by computing the integral in the numerator of (4.12) by using the Fourier representations of and *f*(*t*)=sin *t* and by taking the zero Fourier component of their product, one can see that the only possible frequencies are of the type *ω*=1/(2*k*+1); see Bartuccelli *et al*. (in press) for details. We are implicitly assuming that the integral is non-zero; this must be checked (either analytically or numerically), because if it happens to be zero then we have to go to higher orders until we find a *k*>1 such that(4.13)Other subharmonic solutions appear only for much smaller values of the dissipative parameter, that is for *γ*=*O*(*ϵ*^{p}), with *p*>1. Again we refer to Gentile *et al*. (2007) and Bartuccelli *et al*. (in press *a*) for further details.

The subharmonic solutions turn out to be attractive. The numerical simulations performed in §7 give evidence that, at least for small values of the perturbation parameter *ϵ*, all trajectories in phase space are eventually attracted either to some subharmonic solution or to the periodic solution, which will be investigated in §5. Note that this situation is different from the kind of problem studied in Stoker (1950): first of all our unperturbed system is strongly nonlinear; moreover, the limit cycles described by the attractive periodic solutions are generated by the perturbation itself, so that the problem does not consist in studying the persistence of the self-sustained limit cycle.

## 5. Periodic solution branching off the origin

In this section, we wish to investigate the periodic solution which originates from the zero solution in the presence of the periodic forcing function *f*(*t*)=*f*(*t*+2*π*). For this it is more convenient to write our system as follows:(5.1)where , and again the constant *β* has been absorbed into the constants *ϵ* and *γ*.

The following result will be proved.

*Consider the system* (5.1), *and assume that* . *There exists ϵ*_{0}∝max{*γ*, *Q*_{0}} *such that for* |*ϵ*|<*ϵ*_{0}, *the system* (5.1) *admits a periodic solution x*_{0}(*t*) *with the same period as the forcing, which describes a closed curve* *of diameter O*(*ϵ*/*ϵ*_{0}) *around the origin. Moreover, for any value of ϵ there exists γ*_{0}>0 *such that for γ*>*γ*_{0}, *the periodic solution x*_{0}(*t*) *describes a global attractor*.

The idea is to expand the periodic solution around the origin and then to use a similar strategy to the one of §4. Thus we have(5.2)note that the dissipative parameter *γ* is now included in the ‘unperturbed’ system and therefore it is (in general) not necessarily small. We now introduce the decompositions (5.2) into (5.1) and we perform a similar analysis to that of §5. To first order in *ϵ*, we obtain (remember that the zero-order solution is the origin)(5.3)where *g*(0)=*β*/*α*. This equation is linear non-homogeneous and its time-asymptotic solution (obtained by setting to zero the transient non-periodic part which decays exponentially in time) has Fourier coefficients(5.4)In particular it has the same period as *f*. Note that the first-order solution *x*^{(1)}(*t*) has the same average as the forcing function *f*(*t*) (as it should). It is now instructive to compute the second order in *ϵ* and then the *k*th order, because at second order one can already see the general properties of the structure of the solution at any order. By equating powers of *ϵ*^{2}, we obtain(5.5)which can be easily solved in Fourier space.

By inspecting the structure of equation (5.1), one can see that at any order *k* the linearized part is the same—compare for instance (5.3) with (5.5)—and the perturbation is formed of terms that are periodic functions having the same period as the forcing term *f*(*t*); hence the structure of the equation at order *k* has the form(5.6)where the function *P*^{(k)} is a periodic function having the same period as the forcing term *f*(*t*); hence, by writing(5.7)we find that the Fourier coefficients are given by(5.8)according to (5.6). Note that depends on the coefficients , with *k*′<*k* and , so that it is a well-defined periodic function (as can be checked by induction). The explicit expression for can be found in appendix B.

The convergence of the perturbation series (5.2) can be discussed as in the references Gallavotti & Gentile (2002) and Gentile *et al*. (2007) quoted in §4. The presence of the factor *γ* appearing in the denominator of (5.8) implies that the corresponding radius of convergence in *ϵ* grows linearly in *γ*. More precisely, by iterating (5.8) for every *k* down to *k*=1, one finds that, uniformly in *t*, , where *A* and *B* are constants which depend upon the functions *g*(*y*) and *f*(*t*), but not on *γ*, and *Q*(*γ*)=max{*g*(0)*γ*, *Q*_{0}}, with ; cf. appendix B. Therefore by looking at the power series expansion (5.2), one sees that by choosing *ϵ* judiciously the series converges. The radius of convergence becomes smaller with decreasing *γ* up to a threshold value proportional to , below which it remains constant. For *γ*=0, we can choose *Q*_{0}/*B* as the radius of convergence of the power series, provided *n*^{2}−*β*/*α*≠0 (absence of resonances). Moreover, under this hypothesis, for any fixed *ϵ*—not necessarily small—we have that for *γ* large enough, say for *γ*>*γ*_{0}, there is a periodic solution *x*_{0}(*t*) encircling the origin with the same period as the forcing. We shall see in §6 that there exists *γ*_{1}≥*γ*_{0} such that for *γ*>*γ*_{1} the solution *x*_{0}(*t*) describes a global attractor.

## 6. Global attraction to the periodic solution encircling the origin

In this section, we wish to prove that if the dissipative parameter *γ*>0 is large enough, then every initial condition will evolve so as to tend asymptotically to the periodic solution encircling the origin, whose existence was proved in §5. We start with the perturbed equation in the form(6.1)where , and the parameter *ϵ* (arbitrary, cf. the comments at the end of §5) has been absorbed into the forcing function *f*(*t*). The strategy we follow is to split the solution *x*(*t*)=*x*_{0}(*t*)+*ξ*(*t*), where *x*_{0}(*t*) satisfies (6.1). Thus we must prove that *ξ*(*t*)→0 as time goes to infinity. By introducing *x*(*t*)=*x*_{0}(*t*)+*ξ*(*t*) into (6.1) and simplifying the terms satisfied by the periodic solution, we obtain the equation(6.2)We write it as(6.3)where as usual . Hence, we have to prove the following result.

*Consider the equation* (6.3) *with x*_{0}(*t*), *the periodic solution encircling the origin. Then if the dissipative parameter γ*>0 *is large enough, all trajectories in phase space converge towards the origin asymptotically in time*.

By using the Liouville-like transformation (Bartuccelli *et al*. 2004)we first remove the time dependence from the coefficient of *ξ*. We obtain and , where the prime denotes the derivative with respect to the variable *τ*. Substituting and simplifying we get(6.4)By defining the Hamiltonian for the harmonic oscillator by *H*=(*ξ*^{2}+*y*^{2})/2, one can see that *H*′=−*Γy*^{2}−*Fy*, where and . Now one can see that we can estimate with *L* a positive constant of order 1/*γ*. The last property follows, when explicitly estimating *L*, from the fact that are of order 1/*γ* (see (5.4)). Then one finds *H*′≤−*y*^{2}(*Γ*−*L*/*Γ*)≤0 for every (*ξ*,*y*). Moreover, *H*′=0 only on the *y*-axis provided that *γ* is large enough. Hence, we are in the hypotheses of the theorem of Barbashin–Krasovsky (Krasovsky 1963) and therefore we can conclude that if *γ* is large enough so as to have and *Γ*−*L*/*Γ*>0, it follows that the origin is asymptotically stable and attracts all the solutions for every initial condition in phase space. This completes the proof of the lemma. ▪

## 7. Numerical simulations

We now report on some numerical computations which illustrate theorems 4.1 and 5.1 by considering a special case of equation (2.5) with *f*(*t*)=*A* sin *t*, and choosing *A*=1.2, *α*=40, *β*=1 and *γ*=0.003. Note that (2.5) corresponds to (3.1) with *ϵ*=1. Technically theorems 4.1 and 5.1, for fixed *γ*, have been proved for values of *ϵ* relatively small so as to be sure of the convergence of the power series that arise. However, it is very time consuming to obtain basin of attraction pictures for small values of the dissipative parameter *γ*=*ϵC*. Hence, we have compromised by verifying, within the limitations of finite precision and time numerical computations, that only attracting periodic solutions having periods 2*πn* for *n*=1, 3 and 5 exist, for (3.1) with *ϵ*=1, *ϵ*=0.1 and *ϵ*=0.01. Only in the case *ϵ*=1, we have computed the actual basins of attraction of each of these solutions, figure 2. The figure shows solutions with periods 2*πn* for figure 2*a*, *n*=1; figure 2*b*, *n*=3 and figure 2*c*, *n*=5. Below there are the corresponding basins of attraction. These were computed in the obvious way—by numerically integrating the differential equation starting from a set of initial conditions on a 200×200 grid, and in each case, after allowing the transient to decay sufficiently, deciding which solution has been reached. Of the 40 000 pixels in the figure, 17 335 correspond to the *n*=1 solution, 18 718 to *n*=3 and 3947 to *n*=5. Hence, all 40 000 initial conditions considered lead to one of only these three solutions.

This is in accordance with the analytical results of §§4 to 6. Lemma 3.1 implies that subharmonic solutions are possible corresponding only to orders *q*/*p* with . Furthermore, in the light of the discussion after theorem 4.1, one has *p*=1, and *q* is an odd integer owing to lemma 3.2, so for these parameters, *q* can only be 3 or 5. We also need to check that *γ* is below the maximum of *C*_{0}(*t*_{0}) as defined in equation (4.12). For *f*(*t*)=*A* sin *t*, one finds *C*_{0}(*t*_{0})=*R*(*q*)cos *t*_{0}, with *R*(*q*) a constant depending on *q*. Carrying out the implied integrations numerically, we find that *R*(3)≈–0.02113 and *R*(5)≈0.05174, while (for *ϵ*=1) one has *C*=*γ*=0.003. A similar scenario has been checked to arise for *γ*=0.015 (and the same values of the other parameters).

Besides these subharmonic solutions, there is also the periodic solution branching off the origin. Note that in the case investigated, one has , so that we are far from resonances, and theorem 5.1 applies.

The numerical results show that these three periodic solutions attract almost all trajectories (at least for initial data not too far from the origin).

## 8. Discussion and conclusion

In this paper, we have discussed various properties of the solutions of a nonlinear ODE with periodic forcing, in which the nonlinearity consists, in a mechanical analogy, of a velocity-dependent mass, the restoring force and dissipation being linear. The ODE arises as a model of an electronic circuit containing an inductor whose core saturates. The importance of this circuit is that it is very simple, consisting of only three electronic components, and yet for certain values of the parameters can display chaotic behaviour. Previous studies made use of various simplified models of the saturating inductor: for instance, by using a piecewise constant approximation to the function *g*(*y*) and neglecting hysteresis effects (Chua *et al*. 1982); or by using a more realistic model including hysteresis effects (so that *g* depends on *y* and ) (Deane 1994). The latter results in a third-order non-autonomous ODE (Deane 1994). The ODE studied in this paper is in between these two models, as it still neglects hysteresis but uses a smooth—more realistic—function for the magnetization.

In our analysis, we considered first the unperturbed version of the system in which dissipation and forcing were absent, and showed that here, only periodic solutions are displayed and their period is a monotonically decreasing function of the energy. When the perturbation is present, we have proved the existence of subharmonic solutions—that is solutions whose periods are rational multiples of the period of the forcing—and additionally, of a periodic solution with the same period as the forcing, the latter becoming a global attractor in the presence of sufficient dissipation. Numerical computations have been carried out to illustrate this behaviour and show the relevance of these solutions for the dynamics.

The main limitation in our model is the fact that we neglect hysteresis. On the other hand, as mentioned earlier, taking it into account leads to a much more complicated system, and application of the analysis techniques used here, at best, would require non-trivial generalization. It would be interesting to investigate in a real circuit whether hysteresis can be really neglected for certain inductor core materials. In general, the results presented in Deane (1994) show that hysteresis produces a distortion of the periodic solutions studied in this paper. Moreover, for larger values of the perturbation parameter new bifurcation phenomena appear (such as period doublings) and also chaos.

An analytical description of this behaviour would be highly desirable, but would require dealing with a non-perturbative regime: the mathematics becomes much more involved and new ideas are necessary. Note that our system, in the absence of perturbation (*ϵ*=0), does not possess homoclinic or heteroclinic orbits: indeed, all the unperturbed solutions are periodic. This situation is different from that of pendulum-like systems, where Smale horseshoes transient chaos can manifest itself through homoclinic tangles (Holmes & Marsden 1981/82; Greenspan & Holmes 1984; Guckenheimer & Holmes 1990).

## Acknowledgments

It is a pleasure to acknowledge fruitful discussions with Alexei Ivanov.

## Footnotes

- Received July 12, 2006.
- Accepted May 30, 2007.

- © 2007 The Royal Society