Self-Oscillations in a Generalized Van Der Pol and Rayleigh Oscillator

Borys I. Basok, Volodymyr V. Gotsulenko

Institute of Engineering Thermophysics of NAS of Ukraine, Kyiv, Ukraine.

DOI: 10.4236/jamp.2025.139165                

Abstract

This paper investigates the emergence and transformation of self-oscillatory regimes in a nonlinear dynamic system based on a generalized van der Pol oscillator model. Using asymptotic and bifurcation analysis, conditions for the existence of limit cycles are derived and their stability is examined. Special attention is given to both global and local bifurcations, including the formation of homoclinic loops and Andronov-Hopf bifurcations. The influence of wave resistance and equilibrium position on the qualitative structure of the phase portrait is also analyzed.

Keywords

Self-Oscillations, Limit Cycle, Bifurcation, Wave Resistance

1. Introduction

Self-oscillating systems represent an important class of nonlinear dynamical models widely used in various fields of science and engineering, from radio engineering and mechanics to biophysics and thermophysics. Unlike forced oscillations, self-oscillations arise due to internal mechanisms of self-excitation and are maintained without external periodic forcing. Understanding the nature of such oscillations and their mathematical modeling is of fundamental importance for analyzing stability, synchronization, bifurcations, and transitions to chaos.

Among the most well-known and studied models of self-oscillations are the Van der Pol and Rayleigh equations. The classical Van der Pol equation, proposed in the 1920s to describe electrical oscillations in nonlinear circuits, has the form [1] [2]:

$\ddot{x}-\mu \left(1-x^2\right)\dot{x}+x=0,$

where μ is a parameter indicating the nonlinearity and the strength of the damping. This model exhibits a stable limit cycle resulting from the balance between negative linear and positive nonlinear damping.

The Rayleigh equation, proposed earlier to describe acoustic oscillations in pipes, is written as [3]:

$\ddot{x}+\mu \left(\frac{1}{3}{\dot{x}}^3-\dot{x}\right)+x=0.$

It also describes self-oscillations but with a different type of nonlinear resistance depending only on velocity. Despite differences in physical assumptions, both models share similar qualitative properties and serve as a basis for constructing more complex generalized models.

Both systems can be conveniently represented as first-order systems. For the Van der Pol equation:

$\{\begin{array}{l}\dot{x}=-y,\hfill \\ \dot{y}=\mu \left(1-x^2\right)y+x,\hfill \end{array}$

and for the Rayleigh equation:

$\{\begin{array}{l}\dot{x}=-y,\hfill \\ \dot{y}=-\mu \left(\frac{1}{3}{y}^3-y\right)+x.\hfill \end{array}$

The Rayleigh equation can be transformed into the form of the Van der Pol equation by a change of variables. Let: $x=\dot{z}$ , then $\dot{x}=\ddot{z}$ . Substituting into the Rayleigh equation yields:

$\ddot{z}+\mu \left(\frac{1}{3}{\dot{z}}^3-\dot{z}\right)+z=0,$

which structurally coincides with the Van der Pol equation if we denote $z=x$ , $\dot{z}=y$ . Thus, the Rayleigh equation can be viewed as the Van der Pol equation expressed in variables where $x$ corresponds to velocity rather than position. This observation highlights the structural similarity between the two models and allows the use of common analytical methods.

It is also worth noting that when the nonlinearity parameter $\mu$ is set to zero, both the Van der Pol and Rayleigh equations reduce to the classical harmonic oscillator equation:

$\ddot{x}+x=0,$

which describes simple sinusoidal oscillations with constant amplitude and frequency. This highlights that the nonlinear terms in each model are responsible for the emergence of self-sustained oscillations and complex dynamical behavior.

It is also known that for small positive values of $\mu$ , both the Van der Pol and Rayleigh equations admit orbitally asymptotically stable periodic solutions corresponding to a soft regime of self-excited oscillations. In this regime, the system undergoes a supercritical Andronov-Hopf bifurcation on the phase plane, leading to the birth of a stable limit cycle. This bifurcation scenario is typical for systems with weak nonlinear damping and provides a smooth transition from equilibrium to sustained oscillatory behavior as the nonlinearity parameter increases. Both models describe self-oscillatory processes that arise due to the internal nonlinear structure of the system, without any external excitation. Studying these models allows us to identify the conditions for the emergence of stable limit cycles, bifurcations, and transitions to complex dynamics.

Of particular interest is the regime of hard excitation of self-oscillations, in which a stable equilibrium coexists with a stable limit cycle. This regime is characteristic of a subcritical Andronov-Hopf bifurcation, where the system loses stability not through a smooth emergence of small-amplitude oscillations, but through an abrupt transition to finite-amplitude oscillations. This is accompanied by hysteresis and multistability, which are crucial for applications requiring sharp switching between oscillatory and stationary modes. The relevance of this study stems from the need for more accurate descriptions of complex oscillatory processes occurring in modern engineering and physical systems, where simple models are insufficient. Generalized equations that incorporate multiple types of nonlinearities expand the applicability of classical models and deepen our understanding of self-organization mechanisms in nonlinear media.

2. Formalization of the Dynamic System of the Considered Problem

In this work, we consider a generalized model that combines key features of both the Van der Pol and Rayleigh equations [4] [5]:

$\frac{\text{d}x}{\text{d}t}=\alpha \left(H\left(x\right)-y\right),\frac{\text{d}y}{\text{d}t}=\beta \left(x-\phi \left(y\right)\right),$ (1)

where $H\left(x\right)=\eta -\left[\frac{k_1\epsilon }{{\left(x-x_0\right)}^{\epsilon }}+k_2\epsilon {\left(x-x_0\right)}^{\nu }\right]+k_3\epsilon {\left(x-\xi \right)}^{\mu }$ ; $\eta >0$ ; $\nu >\mu \ge 1$ ; $k_i>0$ $\left(i=\overline{1;3}\right)$ ; $\epsilon \ge 0$ , $\epsilon \to 0$ ; $\alpha >0$ , $\beta >0$ ; $\phi \left(y\right)=\xi {\left(\frac{y-H_0}{\eta -H_0}\right)}^a$ , $0<x_0<\xi$ , $H_0<\eta$ , $a=a_0\epsilon +O\left({\epsilon }^2\right)$ as $\epsilon \to 0$ , $a_0\ge 0$ .

It is convenient to switch to new variables in system (1), assuming

$X=\sqrt{\beta }\left(x-\xi \right),Y=\sqrt{\alpha }\left(y-\eta \right),\tau ={\omega }_{0}t,\text{where}{\omega }_0=\sqrt{\alpha \beta }.$ (2)

In variables (2), the dynamic system (1) takes the following form

$\frac{\text{d}X}{\text{d}\tau }=-Y+\sqrt{\alpha }\left[H\left(\xi +\frac{X}{\sqrt{\beta }}\right)-\eta \right],\frac{\text{d}Y}{\text{d}\tau }=X+\sqrt{\beta }\left[\xi -\phi \left(\eta +\frac{Y}{\sqrt{\alpha }}\right)\right].$ (3)

When $\epsilon =0$ , system (3) is conservative. It is known that certain perturbations can transform a conservative system into a self-oscillating one, with a limit cycle close to one of the closed phase trajectories of the original unperturbed conservative system. In this regard, the following result holds [6].

Proposition 1. As $\epsilon \to 0$ , consider the following dynamic system

$\frac{\text{d}x}{\text{d}t}=-y+\epsilon f_1\left(x,y\right),\frac{\text{d}y}{\text{d}t}=x+\epsilon f_2\left(x,y\right).$

Let ${\ell }_A=\left\{\left(x,y\right):x=A\cos \left(t\right),y=A\sin \left(t\right),0\le t<2\pi \right\}$ , where ${\ell }_A$ is a circle centered at the origin with radius $A$ . Also, let $F\left(A\right)={\displaystyle {\oint }_{{\ell }_A}{f}_1\left(x,y\right)\text{d}y-f_2\left(x,y\right)\text{d}x}$ . Then, if the function $F\left(A\right)$ has a simple positive root $A^{\text{*}}$ , the system has a limit cycle ${\text{Γ}}_{\epsilon }\simeq {\ell }_{A^{\text{*}}}$ for small $\epsilon >0$ , which is stable if ${\frac{\text{d}F}{\text{d}A}|}_{A=A^{\text{*}}}<0$ and unstable if ${\frac{\text{d}F}{\text{d}A}|}_{A=A^{\text{*}}}>0$ .

It is easy to verify that, up to terms of order $O\left({\epsilon }^2\right)$ as $\epsilon \to 0$ , the following equalities hold:

$\begin{array}{l}\sqrt{\alpha }\left[H\left(\xi +\frac{X}{\sqrt{\beta }}\right)-\eta \right]\\ =\epsilon \sqrt{\frac{\alpha }{{\beta }^{\mu }}}\left[k_3{X}^{\mu }-k_2{\left(\xi +\frac{X}{\sqrt{\beta }}-x_0\right)}^{\nu }-k_1\right]+O\left({\epsilon }^2\right),\end{array}$

$\begin{array}{c}\sqrt{\beta }\left[\xi -\phi \left(\eta +\frac{Y}{\sqrt{\alpha }}\right)\right]=\xi \sqrt{\beta }\left[1-{\left(1+\frac{Y}{\sqrt{\alpha }\left(\eta -H_0\right)}\right)}^a\right]\\ =-\epsilon a_0\xi \sqrt{\beta }\ln \left(1+\frac{Y}{\sqrt{\alpha }\left(\eta -H_0\right)}\right)+O\left({\epsilon }^2\right).\end{array}$

In the considered problem:

$f_1\left(X,Y\right)=\sqrt{\frac{\alpha }{{\beta }^{\mu }}}\left[k_3{X}^{\mu }-k_2{\left(\xi +\frac{X}{\sqrt{\beta }}-x_0\right)}^{\nu }-k_1\right],$ (4)

$f_2\left(X,Y\right)=-a_0\xi \sqrt{\beta }\ln \left(1+\frac{Y}{\sqrt{\alpha }\left(\eta -H_0\right)}\right).$

Thus, according to (4), for system (3), we obtain

$\begin{array}{c}F\left(A\right)={\displaystyle {\oint }_{{\ell }_A}{f}_1\left(X,Y\right)\text{d}Y}-f_2\left(X,Y\right)\text{d}X\\ ={\displaystyle {\oint }_{{\ell }_A}\sqrt{\frac{\alpha }{{\beta }^{\mu }}}\left[k_3{X}^{\mu }-k_2{\left(\xi +\frac{X}{\sqrt{\beta }}-x_0\right)}^{\nu }-k_1\right]\text{d}Y}\\ +a_0\xi \sqrt{\beta }\ln \left(1+\frac{Y}{\sqrt{\alpha }\left(\eta -H_0\right)}\right)\text{d}X.\end{array}$ (5)

To compute the line integral (5), we need the following preliminary relation: $\forall p\in \mathbb{N}\cup \left\{0\right\}$

$I_p={\displaystyle {\int }_0^{2\pi }{\cos }^p\left(t\right)\text{d}t}=2\sqrt{\pi }\frac{{\left(-1\right)}^p+1}{p}\frac{\Gamma \left(\frac{p+1}{2}\right)}{\Gamma \left(\frac{p}{2}\right)},$ (6)

where $\text{Γ}\left(x\right)$ is the Euler gamma function. Taking into account relation (6), we obtain

$\begin{array}{l}{\displaystyle {\oint }_{{\ell }_A}{f}_1\left(X,Y\right)\text{d}Y}\\ =\sqrt{\frac{\alpha }{{\beta }^{\mu }}}A\left[k_3{A}^{\mu }{I}_{\mu +1}-k_2{\displaystyle {\int }_0^{2\pi }\cos \left(t\right){\left(\xi -x_0+\frac{A}{\sqrt{\beta }}\cos \left(t\right)\right)}^{\nu }\text{d}t}\right]\\ =k_3{I}_{\mu +1}\sqrt{\frac{\alpha }{{\beta }^{\mu }}}{A}^{\mu +1}-k_2\sqrt{\frac{\alpha }{{\beta }^{\mu }}}A{\displaystyle \sum }_{k=0}^{\nu }{\left(\xi -x_0\right)}^{\nu -k}{I}_{k+1}\frac{A^k}{\sqrt{{\beta }^k}}.\end{array}$ (7)

Further,

${\displaystyle {\oint }_{{\ell }_A}{f}_2\left(X,Y\right)\text{d}X}=-a_0\xi \sqrt{\beta }{\displaystyle {\oint }_{{\ell }_A}\ln \left(1+\frac{Y}{\sqrt{\alpha }\left(\eta -H_0\right)}\right)\text{d}X}.$ (8)

Let us compute the integral depending on the parameter

$I\left(r\right)={\displaystyle {\int }_0^{2\pi }\ln \left(1+r\sin \tau \right)\sin \tau \text{d}\tau },$

using Leibniz’s rule for differentiation under the integral sign. We have

$I\left(0\right)=0,\frac{\text{d}I}{\text{d}r}={\displaystyle {\int }_0^{2\pi }\frac{{\sin }^2\tau }{1+r\sin \tau }\text{d}\tau }=\frac{2\pi }{r^2}\left[\frac{1}{\sqrt{1-r^2}}-1\right]\Rightarrow$

$I\left(r\right)={\displaystyle {\int }_0^r\frac{\text{d}I}{\text{d}r}\text{d}r}=2\pi {\displaystyle {\int }_0^r\frac{1}{r^2}\left[\frac{1}{\sqrt{1-r^2}}-1\right]\text{d}r}=\frac{2\pi }{r}\left[1-\sqrt{1-r^2}\right].$

Thus, the expression for the line integral (8) can be written as

${\displaystyle {\oint }_{{\ell }_A}{f}_2\left(X,Y\right)\text{d}X}=2\pi a_0{\omega }_0\xi \left(\eta -H_0\right)\left(1-\sqrt{1-\frac{A^2}{\alpha {\left(\eta -H_0\right)}^2}}\right).$

Taking into account the obtained relations, the expression for the function $F\left(A\right)$ in the considered problem can finally be written in the form

$F\left(A\right)=q_1{A}^{\mu +1}-q_{2}A{\displaystyle \sum }_{k=0}^{\nu }{q}_{3,k}{A}^k-q_4\left(1-\sqrt{1-q_5{A}^2}\right),$ (9)

where $q_1=k_3{I}_{\mu +1}\sqrt{\frac{\alpha }{{\beta }^{\mu }}}$ , $q_2=k_2\sqrt{\frac{\alpha }{{\beta }^{\mu }}}$ , $q_{3,k}=I_{k+1}\frac{{\left(\xi -x_0\right)}^{\nu -k}}{\sqrt{{\beta }^k}}$ for $k=0,\cdots ,\nu$ , $q_4=2\pi a_0{\omega }_0\xi \left(\eta -H_0\right)$ , $q_5=\frac{1}{\alpha {\left(\eta -H_0\right)}^2}$ .

The proof of the existence of positive roots of function (9) in the general case is a separate and rather complex problem. However, in the case $\mu =1$ , the proof can be obtained by elementary means. It is easy to verify that a sufficient condition for the existence of a positive root of Equation (9) when $\mu =1$ is the fulfillment of the following inequalities:

$q_1-q_2{q}_{3,1}-\frac{1}{2}{q}_4{q}_5>0,q_1-q_2{\displaystyle \sum }_{k=0}^{\nu }{q}_{3,k}{q}_5^{\frac{\mu -k}{2}}\le q_5^{\frac{\mu +1}{2}}.$ (10)

Indeed, as $A\to +0$ , the asymptotic expansion holds:

$F\left(A\right)=\left(q_1-q_2{q}_{3,1}-\frac{1}{2}{q}_4{q}_5\right)A^2+O\left(A^4\right),$

from which, using the first inequality in (10), it follows that $F\left(A\right)>0$ as $A\to +0$ . Furthermore, using the second inequality in (10), it is easy to verify that $F\left(\frac{1}{\sqrt{q_5}}\right)\le 0$ , and therefore, by the Bolzano-Cauchy intermediate value theorem for continuous functions, there exists a positive root $0<A^{\text{*}}\le \frac{1}{\sqrt{q_5}}$ of Equation (9). It is also directly verified that ${\frac{\text{d}F}{\text{d}A}|}_{A=A^{\text{*}}}>0$ , i.e., the limit cycle ${\text{Γ}}_{\epsilon }\simeq {\ell }_{A^{\text{*}}}$ is unstable. Further, although numerical, analysis shows that when varying the coefficients of system (1), the function $F\left(A\right)$ acquires a second root. The complete picture of possible phase portrait bifurcations of system (1) is presented in the next section.

Remark. Note that the canonical form of the dynamic system (1) considered above includes, as a special case, the classical dynamic systems of the Rayleigh and van der Pol oscillators. Indeed, assuming $H\left(x\right)=0$ , $\phi \left(y\right)=\mu \left(\frac{1}{3}{y}^3-y\right)$ and $\alpha =\beta =1$ in (1), we obtain the dynamic system of the Rayleigh oscillator, which in turn is reduced to the van der Pol system by changing the variables as shown in Section 1.

3. Analysis of Local and Global Bifurcations

As the parameters of the considered dynamic system (1) vary, its phase portrait undergoes a series of both local and global bifurcations, particularly those associated with the emergence of a homoclinic curve [7]. Let us examine these bifurcations in more detail. First, note that system (1) has two equilibrium points: $O_1\left(\xi ,\eta \right)$ —a focus, and $O_2\left({\xi }_0,H\left({\xi }_0\right)\right)$ —a saddle point, where $x_0<{\xi }_0\le \xi$ .

We consider the wave resistance $Z$ as the bifurcation parameter. Wave resistance (impedance) is a physical property of a medium that characterizes the ratio between the amplitudes of the electric and magnetic fields (in electrodynamics), between voltage and current (in transmission lines), or between pressure and particle velocity (in acoustics), during wave propagation. Mathematically, it is defined differently depending on the type of wave and the medium. In nonlinear self-oscillatory systems, this impedance may be time-dependent, amplitude-sensitive, and nonlinear, reflecting the complex dynamics of energy transfer between the source and the system. In the context of the considered dynamical system (1), the wave resistance can be represented in the following form $Z=\sqrt{\frac{\beta }{\alpha }}$ [5]. For sufficiently large values of $Z$ , system (1) has two special trajectories (Figure 1(a)): the separatrix $K_1$ of the saddle $O_2$ , which is its stable invariant manifold, and the heteroclinic curve $K_2$ , which asymptotically approaches the focus $O_1$ as $t\to +\infty$ and the saddle $O_2$ as $t\to -\infty$ . In this case, system (1) has a single attractor - the stable focus $O_1$ .

Figure 1. Phase portraits of the system (1).

As the value of $Z$ decreases, at a certain first bifurcation value $Z=Z_1$ , a homoclinic curve $K_3$ appears, forming a loop of the separatrix of the saddle $O_2$ . At this point, the focus $O_1$ remains an attractor with a basin of attraction bounded by the homoclinic curve $K_3$ (Figure 1(b)). At the bifurcation point $Z=Z_1$ , the phase portrait is non-rough (structurally unstable), and with further decrease $Z<Z_1$ , the homoclinic curve breaks (Figure 1(c)), giving rise to an unstable limit cycle ${\text{Γ}}_1$ .

So far, the bifurcations have been global and are associated with qualitative topological changes in the entire phase portrait of the system (1). However, with further decrease in $Z$ , the following local bifurcations also occur. At a certain value $Z=Z_2$ , the focus $O_1$ becomes unstable and a stable limit cycle ${\text{Γ}}_2$ branches off from it (Figure 1(d)). This bifurcation is local and represents a supercritical Andronov-Hopf bifurcation [6] [8]. Indeed, the linearization matrix of system (3) at the origin has the following pair of complex conjugate eigenvalues:

${\lambda }_{1,2}\left(Z\right)=\frac{1}{2}\left(\frac{H^{\prime }\left(\xi \right)}{Z}-{\phi }^{\prime }\left(\eta \right)Z\right)\pm i\sqrt{1-\frac{1}{4}{\left(\frac{H^{\prime }\left(\xi \right)}{Z}+{\phi }^{\prime }\left(\eta \right)Z\right)}^2},$ (11)

from which, for the bifurcation value $Z=Z_2$ defined by the condition ${\text{Re}\left[{\lambda }_{1,2}\left(Z\right)\right]|}_{Z=Z_2}=0$ , we obtain the expression $Z_2=\sqrt{\frac{H^{\prime }\left(\xi \right)}{{\phi }^{\prime }\left(\eta \right)}}$ , and

$\frac{\partial }{\partial Z}{\mathrm{Re}\left[{\lambda }_{1,2}\left(Z\right)\right]|}_{Z=Z_2}=2\phi \left(\eta \right)\ne 0,$

${\mathrm{Im}\left[{\lambda }_{1,2}\left(Z\right)\right]|}_{Z=Z_2}=\pm \sqrt{1-H^{\prime }\left(\xi \right){\phi }^{\prime }\left(\eta \right)}\ne 0.$

The amplitude of the emerging limit cycle ${\text{Γ}}_2$ is determined in order by the square root of the supercriticality, i.e., the following asymptotic holds:

$\text{Amplitude}\left({\text{Γ}}_2\right)=\text{Const}{\left|Z-\sqrt{\frac{H^{\prime }\left(\xi \right)}{{\phi }^{\prime }\left(\eta \right)}}\right|}^{1/2}\left(1+o\left(1\right)\right)\text{as}Z\to Z_2.$

Near the bifurcation point (for $Z>Z_2$ ), the period of self-oscillations defined by the limit cycle ${\text{Γ}}_2$ does not depend on the magnitude of the supercriticality $\left|Z-Z_2\right|$ , and the following asymptotic expansion holds:

$T=\frac{2\pi }{\sqrt{1-H^{\prime }\left(\xi \right){\phi }^{\prime }\left(\eta \right)}}+o\left(1\right)\text{as}Z\to Z_2.$

With further decrease in $Z$ , the cycles ${\text{Γ}}_1$ and ${\text{Γ}}_2$ approach each other, and at a certain value $Z=Z_3$ , a bifurcation of merging of the stable ${\text{Γ}}_2$ and unstable ${\text{Γ}}_1$ cycles occurs. At the bifurcation point (i.e., at $Z=Z_3$ ), a double (semi-stable) limit cycle ${\text{Γ}}_3$ arises [8] (Figure 1(e)). For $Z<Z_3$ , the cycle ${\text{Γ}}_3$ breaks, the focus $O_1$ becomes a repeller, and system (1) has no attractors (Figure 1(f)).

4. Subcritical Andronov-Hopf Bifurcation

In this case, the parameter $\xi$ is considered as a bifurcation parameter and is defined by the stable focus $O_1$ of the dynamic system (1). As the value of $\xi$ decreases, the basin of attraction of the focus $O_1$ narrows. We show that in this case a reverse (subcritical) Andronov-Hopf bifurcation occurs [6] [8]. Indeed, the pair of complex conjugate eigenvalues (11) should now be considered as functions of the parameter $\xi$ . In this case,

$\mathrm{Re}\left[{\lambda }_{1,2}\right]=\frac{1}{2}\left(\frac{H^{\prime }\left(\xi \right)}{Z}-{\phi }^{\prime }\left(\eta \right)Z\right)=\frac{{\phi }^{\prime }\left(\eta \right)}{2Z}\left(\frac{H^{\prime }\left(\xi \right)}{{\phi }^{\prime }\left(\eta \right)}-Z^2\right).$

Furthermore, considering that ${\phi }^{\prime }\left(\eta \right)=\frac{a\xi }{\eta -H_0}>0$ , we obtain

$\begin{array}{l}\mathrm{sgn}\mathrm{Re}\left[{\lambda }_{1,2}\right]=\mathrm{sgn}\Psi \left(\xi \right),\text{where}\Psi \left(\xi \right)=\frac{H^{\prime }\left(\xi \right)}{{\phi }^{\prime }\left(\eta \right)}-Z^2\\ =\frac{\epsilon \left(\eta -H_0\right)}{a\xi }\left[\frac{k_1\epsilon }{{\left(\xi -x_0\right)}^{1+\epsilon }}-k_1\nu {\left(\xi -x_0\right)}^{\nu -1}+k_3\left(1-\mathrm{sgn}\left(\mu -1\right)\right)\right]-Z^2.\end{array}$

It is easy to verify that the function $\Psi \left(\xi \right)$ is monotonically decreasing and has a positive root ${\xi }^{\text{*}}>0$ such that $\Psi \left({\xi }^{*}\right)=0$ , and for $0<\xi <{\xi }^{\text{*}}$ the inequality $\Psi \left(\xi \right)>0$ holds, while for $\xi >{\xi }^{\text{*}}$ the reverse inequality $\Psi \left(\xi \right)<0$ holds. Thus, when the parameter $\xi$ crosses the bifurcation value ${\xi }^{\text{*}}$ from right to left (i.e., decreasing), the stable focus $O_1$ becomes unstable. Moreover, before the bifurcation point in the region $\xi >{\xi }^{\text{*}}$ , the function $F\left(A\right)$ has a positive root $A=A^{\text{*}}>0$ , for which the inequality ${\frac{\text{d}F}{\text{d}A}|}_{A=A^{\text{*}}}>0$ holds. According to Proposition 1, in this case the stable focus $O_1$ is surrounded by an unstable limit cycle ${\text{Γ}}_1\approx {\ell }_{A^{\text{*}}}$ . This cycle serves as the basin of attraction for the attractor $O_1$ , and as $\xi \to {\xi }^{\text{*}}+0$ , it shrinks and collapses to a point. At $\xi ={\xi }^{\text{*}}$ , a crisis (disappearance) of the attractor $O_1$ occurs, and for $\xi <{\xi }^{\text{*}}$ it becomes a repeller.

5. Conclusion

This paper investigates the dynamics of a generalized model that unifies the Van der Pol and Rayleigh oscillators. Conditions for the existence of limit cycles are derived, and their stability is analyzed. Special emphasis is placed on the subcritical Andronov-Hopf bifurcation and the regime of hard excitation of self-oscillations, characterized by hysteresis and multistability. The study also explores global bifurcations involving the formation of homoclinic loops and examines how the position of the equilibrium influences the structure of the phase portrait. It is important to note that the above analytical derivations and asymptotic analysis are valid primarily for small values of parameter $\epsilon \to 0$ . For larger values of this parameter, the (1) system’s behavior may deviate significantly from the predicted trends. In such cases, numerical approaches become essential for further investigation. Methods such as the Runge-Kutta algorithm or other finite-difference schemes for solving ordinary differential equations can be effectively employed to analyze the system’s (1) dynamics. The results contribute to a deeper understanding of nonlinear dynamics and enhance the analytical and control capabilities for self-oscillatory regimes in relevant hydrodynamic systems [5]. Let us also note that the considered dynamic system (1) is a mathematical model describing thermoacoustic self-oscillations and self-oscillations of vibration combustion in various heat and power systems [9]. Thus, understanding the qualitative structure of the phase portrait of this dynamic system, in particular, possible bifurcations, will improve the technological performance of the corresponding industrial units [10].

Acknowledgements

The authors acknowledge that the research is supported by the National Research Foundation of Ukraine under project No. 2025.06/0054, “Absorption and Prevention of Infrared Electromagnetic Radiation Propagation and Development of Means for Thermal Camouflage”.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Andronov, A.A., Vitt, A.A. and Khaikin, S.E. (1966) Theory of Oscillators. Princeton Univ. Press.
[2]	Dorodnitsyn, A.A. (1947) Asymptotic Solution of the Van der Pol Equation. Journal of Applied Mathematics and Mechanics (Prikladnaya Matematika i Mekhanika), 11, 313-328.
[3]	Cesari, L. (1998) Asymptotic Behavior and Stability Problems in Ordinary Differential Equations. Dover Publications.
[4]	Gotsulenko, V.V. (2014) Self-Oscillations in Implicit Singularly Perturbed Dynamical Systems on the Plane. Nelineinaya Dinamika, 10, 157-175.[CrossRef]
[5]	Gotsulenko, V.V. and Gotsulenko, V.N. (2012) Distinctive Feature of Self-Oscillations (Surging) of Impeller Pumps. Journal of Engineering Physics and Thermophysics, 85, 125-130.[CrossRef]
[6]	Arnold, V.I. (1973) Ordinary Differential Equations. MIT Press.
[7]	Kuznetsov, S.P. (2001) Dynamical Chaos. Fizmatlit. (In Russian)
[8]	Neimark, Yu.I. (1966) The Method of Point Mappings in the Theory of Nonlinear Oscillations. Nauka. (In Russian)
[9]	Basok, B.I. and Gotsulenko, V.V. (2015) Thermo-Hydrodynamic Instability of Heat Carrier Flow. KALITA Publ. (In Russian)
[10]	Basok, B.I., Gotsulenko, V.V. and Avramenko, A.A. (2019) Mechanisms of Thermophysical Instability of Heat Carrier Flows. Institute of Technical Thermo-Physics, NAS of Ukraine. (In Russian)