An instrumental insight for a periodic solution of a fractal Mathieu–Duffing equation

Abstract

The primary goal of the present study is to investigate how to obtain a periodic solution for a fractal Mathieu–Duffing oscillator. To achieve this, the fractal oscillator in the fractal space has been transformed into a damping Mathieu–Duffing equation in the continuous space by employing a new modification of He’s definition of the fractal derivative. The required analytical periodic solution has been based on the rank upgrade technique (RUT) presented. The RUT successfully generates a periodic solution without sacrificing the damping coefficient by creating an alternate equation, aside from any challenges in managing the impact of the linear damping component. The homotopy perturbation method (HPM) has been used to find the required periodic solution for the alternate equation. A comparison of the numerical solutions of the original equation and the alternative equation showed good agreement. The stability behavior in the non-resonance case as well as in the sub-harmonic resonance case has also been discussed. Further, another method, “the non-perturbative approach”, that deals with the obtained equation has been introduced.

Introduction

The vibration of oscillatory systems is reduced or stopped by a damping action, which can have an internal or external impact. Friction can be a factor in damping, but in general, the definition of damping should not be confused with friction. This damping is divided into two types, the first of which does not depend on the loss of energy in oscillation, which is evident in biological systems and some mechanical systems such as the suspension system of bicycles. ¹ The second type depends mainly on the loss and dispersion of the energy stored in the oscillation, which appears clearly in physical systems such as the dispersion and absorption of light in optical oscillators and resistance in electronic oscillators, as well as some mechanical systems, such as viscous damping processes, where mechanical oscillations are slowed down due to the resistance resulting from fluid viscosity. ² The damping strength depends on a parameter called the damping ratio, which is an important measure that describes the rate at which vibrations decay from one feedback to another in the system after a perturbation occurs. If the damping ratio parameter is equal to zero, the system is called un-damped; if it is less than one, it is called un-damped; if it is equal to one, it is called critically damped; and if it is greater than one, it is called over-damped.^3,4 Oscillating systems have important applications in structural engineering, such as erecting a tall building in the wind, electrical engineering, as well as control engineering, such as controlling the speed of an electric motor. Therefore, applications of the damped nonlinear oscillator are of great importance in plasma physics and mechanical engineering.

In 1868, Émile Léonard Mathieu ⁵ introduced his equation while investigating the vibration of the elliptical drumheads. Since then, Mathieu’s equation has gained great importance because of its significant applications in a variety of areas, including general relativity, quantum mechanics, optics, applied mathematics, physics, and mechanical engineering. ⁶ An example of these applications is the analysis of dynamic systems that contain periodic forces in time or space, as well as the study of some boundary value problems in elliptical symmetry. ⁷ The stability zone analysis of the damped Mathieu model was studied by many researchers, such as Taylor and Narendra, ⁸ Gunderson et al., ⁹ Turyn, ¹⁰ Pradeep and Shrivastava, ¹¹ Allievi and Soudack, ¹² and Barrios et al. ¹³ The same study was also carried out in the presence of a time delay effect in Refs. ^14–17 Moreover, El-Dib and Elgazery ¹⁸ proposed an analytical solution for the damped Mathieu model by using a modulating homotopy perturbation technique. Recently, Li-He’s method developed as a novel modified homotopy perturbation technique to deal with this type of problems. ¹⁹ Moreover, He and El-Dib ²⁰ introduced another enhanced homotopy perturbation mechanism for solving the Klein–Gordon equation. Moreover, the rank upgrading mechanism is extended in Refs. ^21,22 for solving some nonlinear oscillations via homotopy perturbation technique. Also, for the Mathieu–Duffing oscillator with fractional order, the super-harmonic resonance case, and dynamical delayed response were studied by Niu et al. ²³ and Wen et al., ²⁴ respectively.

Theoretical applications of fractal analysis that have been of interest to various studies include the physical phenomena of hierarchical structure, thermal conductivity in micro-channels, diffusion, dynamics of the atmosphere and ocean, petroleum and nuclear engineering, meteorology, material design, and thermo-elasticity^25–33. Through a porous structure, Sheng et al. ³⁴ investigated a fractal model, and Miao et al. ³⁵ looked into another fractal model of shale. Moreover, for a forced damped fractal oscillator, Elías-Zúñiga et al. ³⁶ computed a precise solution during a steady state. He et al. ³⁷ investigated a rough solution for the Toda oscillator in fractal space employing the technique of homotopy perturbation. He et al.^38,39 also developed a contemporary fractal modification and used the Hamilton principle to determine a periodic solution for a fractal Duffing oscillator. Additionally, numerous publications have been written about the research of micro/nano-electromechanical fractal systems; for examples, see Refs. ^40–42 Further, Wang ⁴³ used the Fourier series to provide a novel method for a class of discontinuous fractal nonlinear oscillators in the context of microgravity space. El-Nabulsi et al. ⁴⁴ also developed a novel local derivative operator with fractals across Lagrangian/Hamiltonian functions for further research of the dynamical systems. The two-scale fractal theory was developed by Anjum and his school^45,46 to study the mathematical models of both population and tsunami waves. The same fractal theory was also employed by He and El-Dib, ⁴⁷ Feng, ⁴⁸ and Lu and Chen ⁴⁹ to determine the analytical solutions for the fractal Shabat-Zhiber, Duffing, and Yao-Cheng oscillators. Furthermore, the two-scale fractal theory has been utilized to explore several crucial models, including drinking of alcohol, some rheological issues, heat transfer, and porous media (for more details, see Refs ^50,51). On the other hand, a novel technique for computing the fractal dimensions of porous medium was presented by Xia et al. ⁵² Recently, He and Liu ⁵³ also discussed how to calculate fractal dimensions.

As mentioned in Refs^54,55, He’s frequency formulation could be extended to the nonlinear fractal oscillators, which constitutes the motivation for the work to obtain a periodic solution for the Mathieu–Duffing fractal oscillator by converting it into a damping Mathieu–Duffing equation in the continuous space via a new modification of He’s definition of the fractional derivative. A surrogate equation is produced using the rank upgrade technique (RUT), despite any difficulties in controlling the impact of the linear damping component. These were solved using two distinct approaches.

Parametric excitation in fractal space

The investigation of the Mathieu–Duffing equation has a great history, and it is mostly considered that most properties of its solutions have been perfectly elucidated. ⁶ In the present study, the properties of the solution to the fractal Mathieu–Duffing equation will be investigated. In the fractal space, it has the form

$\frac{dy}{d{t}^{2\alpha }}+\left({\widehat{\omega }}_0^2+2{\delta }_0\mathit{cos}\Omega t\right)y+Q_0{y}^3=0;0<\alpha <1,$

(1)

where ${\widehat{\omega }}_0$ , ${\delta }_0$ , $\Omega$ , and $Q_0$ are the non-zero coefficients, oscillator natural frequency, amplitude, parametric force, and Duffing coefficient, respectively. This system is assumed to depend on the following initial conditions

$y\left(0\right)=A,\frac{d}{d{t}^{\alpha }}y\left(0\right)=0.$

(2)

The fractal operator obeys the modification of He’s definition of the fractal derivative ⁵¹ . This modification was established by El-Dib and Elgazery^56,57 and El-Dib et al. ⁵⁸ and was to transform the oscillator in the fractal space into the continuous space. It has the form

$\frac{dy}{d{t}^{\alpha }}=S^{\alpha }\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)y+S^{\alpha -1}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)\dot{y};0<\alpha <1,$

(3)

where the parameter $S$ represents the porosity of the medium and depends on the fractal order $\alpha$ . Under the modification, equation (3) can transform equation (1) into

$\ddot{y}+2\mu \dot{y}+\left({\omega }_0^2+2\delta \mathit{cos}\Omega t\right)y+Q{y}^3=0,$

(4)

One of the most important features of the transformation (3) is that it deals directly with the operator, unlike the previous works that deal with the variable ³⁶ .

The resulting equation (4) is a damping Mathieu–Duffing equation in continuous space, which represents the equivalent form of the original equation (1). This alternative equation has new coefficients defined as

$\mu =\frac{1}{2}S\mathit{cot}\left(\frac{1}{2}\pi \alpha \right),$

(5)

${\omega }_0^2=4{\mu }^2+\frac{{\widehat{\omega }}_0^2}{{\nu }^2};\nu =S^{\alpha -1}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right),$

(6)

$\delta =\frac{{\delta }_0}{{\nu }^2\left(\alpha \right)},Q=\frac{Q_0}{{\nu }^2\left(\alpha \right)}.$

(7)

Also, the fractal initial conditions (2) have been transformed to

$y\left(0\right)=A,\dot{y}\left(0\right)=-2\mu A.$

(8)

In the limit case ( $\alpha \to 1$ ), it should be emphasized that equation (4) will be reduced to equation (1) $\alpha =1$ . But when the fractal parameter $\alpha$ is less than one, equation (1) will produce a periodic solution in the fractal space, while equation (4) reads as having a damping solution in the continuous space.

The computed fractal parameter $S$ that appears in the transformation (3) can be done in the comparison between the original non-parametric with the alternative linear frequencies ${\widehat{\omega }}_0^2$ and ${\omega }_0^2$ , respectively, which gives

$\mu =\frac{{\widehat{\omega }}_0}{2\nu }\sqrt{\left({\nu }^2-1\right)}.$

(9)

Employing equations (5) and (6) into equation (9) leads to

$S=\frac{{\widehat{\omega }}_0}{\nu }\sqrt{\left({\nu }^2-1\right)}\mathit{tan}\left(\frac{1}{2}\pi \alpha \right).$

(10)

Since the original nonlinear Mathieu equation is defined without the linear damping coefficient $\mu$ , then its solutions behave as periodic ones. Because the damping force is explicitly included in the solution of the corresponding equation (4) in continuous space, damping behavior results are obtained. To return to the original case with the periodic behavior, we may follow El-Dib et al. ⁵⁹ to establish another new alternative equation in the continuous space. The concept of this method depends on deriving an alternative equation free of the linear damping part.

Alternating technique leads to a periodic solution

El-Dib et al.^21,22,59 succeeded in creating a new form that does not explicitly show damping forces. They used the rank-upgrading method to present equation (4) in its alternative form of a fourth-order system having parametric excitation forces. Then they proceeded to replace the first-order and the third-order derivatives with their equivalent expressions. The procedure may be followed as follows:

By differentiating equation (4) concerning the variable $t$ , one finds

$y^{\left(3\right)}+2\mu \ddot{y}+\left({\omega }_0^2+2\delta \mathit{cos}\Omega t+3Q{y}^2\right)\dot{y}-2\delta \Omega \mathit{sin}\Omega ty=0.$

(11)

A double differentiation of equation (4) concerning the variable $t$ gives

$y^{\left(4\right)}+2\mu y^{\left(3\right)}+\left({\omega }_0^2+2\delta \mathit{cos}\Omega t+3Q{y}^2\right)\ddot{y}-4\delta \Omega \mathit{sin}\Omega t\dot{y}+\left(6Q{\dot{y}}^2-2\delta {\Omega }^2\mathit{cos}\Omega t\right)y=0.$

(12)

By using the value of $y^{\left(3\right)}$ from equation (11) and the value of $\dot{y}$ from equation (4), thus equation (12) can be reduced to the next structure

$y^{\left(4\right)}+\left(P_0\left(t\right)+P_2\left(t\right)y^2+\frac{3Q^2}{{\mu }^2}{y}^4\right)\ddot{y}+\left(P_1\left(t\right)+P_3\left(t\right)y^2+\left(Q{P}_2\left(t\right)-1\right)y^4+\frac{3Q^3}{2{\mu }^2}{y}^6+\frac{3Q}{2{\mu }^2}{\ddot{y}}^2\right)y=0.$

(13)

Equation (13), a nonlinear fourth-order differential equation with a parametric coefficient, is free of the difficulties due to the linear damping term and depends on the following initial conditions

$y\left(0\right)=A,\dot{y}\left(0\right)=-2\mu A,\ddot{y}\left(0\right)=A\left(4{\mu }^2-{\omega }_0^2-2\delta -Q{A}^2\right),and\hspace{0.17em}{y}^{\left(3\right)}\left(0\right)=4\mu A\left({\omega }_0^2+2\delta +2Q{A}^2-2{\mu }^2\right).$

(14)

where $P_i\left(t\right);i=\mathrm{0,1,2,3}$ are periodic functions defined as

$\begin{array}{c}{P}_0\left(t\right)=\frac{2}{{\mu }^2}\left({\omega }_0^2{\mu }^2-2{\mu }^4+2\delta \Omega \mu \mathit{sin}\Omega t+2\delta {\mu }^2\mathit{cos}\Omega t\right),\\ P_1\left(t\right)=\frac{1}{\mu }\left[\mu \left(2{\delta }^2+{\omega }_0^4\right)+2\delta \Omega \left(2{\mu }^2+{\omega }_0^2\right)\mathit{sin}\Omega t+2\delta \mu \left(2{\omega }_0^2-{\Omega }^2\right)\mathit{cos}\Omega t+2{\delta }^2\Omega \mathit{sin}2\Omega t+2{\delta }^2\mu \mathit{cos}2\Omega t\right],\\ P_2\left(t\right)=\frac{3Q}{{\mu }^2}\left({\omega }_0^2+2{\mu }^2+2\delta \mathit{cos}\Omega t\right),\\ P_3\left(t\right)=\frac{1}{{\mu }^2}\left[\frac{Q}{2}\left(6{\delta }^2+8{\omega }_0^2{\mu }^2+3{\omega }_0^4\right)+2\delta Q\mu \Omega \mathit{sin}\Omega t+2\delta Q\left(4{\mu }^2+3{\omega }_0^2\right)\mathit{cos}\Omega t+3{\delta }^{2}Q\mathit{cos}2\Omega t\right].\end{array}\}.$

(15)

To ensure that equation (13) is indeed the alternative form of equation (4), the comparison of their numerical solutions can be done for $\alpha <1$ .

A comparison of the numerical solutions of the original and alternative equations is displayed in Figures 1–3 which shows good agreement. This prompts us to trust the possibility of dealing with the alternative equation, addressing how to solve it, and then using its results as a reliable indicator. Further, as we have seen in Figures 1–3, the fractal parameter $\alpha$ has a significant impact on the motion, and minimizing it has a damping effect on the oscillator movement.

Figure 1.

Comparison of the numerical solutions of both equations (4) and (13) for a system having.

$Q_0=0.1,{\delta }_0=0.1,{\widehat{\omega }}_0=1,A=1\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}\alpha =0.9$

Figure 2.

Comparison of the numerical solutions of both equations (4) and (13) for the same system specified in Figure 1 with $\alpha =0.7$ .

Figure 3.

Comparison of the numerical solutions of both equations (4) and (13) for the same system given in Figure 1 except that $\alpha =0.5$ .

The alternative fourth-order equation (13) has multiple forms of periodic coefficients, which represents a complicated problem to handle analytically. To calculate an analytical solution for this equation, we may reduce its rank to become a second-order one. This can be accomplished by using the HPM. This will be achieved in the following section.

Constructing a homotopy equation and its solution

In the present section, reducing the rank of equation (13) is well done. Also, the stability behavior will be discussed in the non-resonance and sub-harmonic resonance cases. The HPM should be used for these purposes.

$\left(-{\omega }_0^{4}y\right)$ is added to both sides of equation (13), or in other words, an auxiliary linear part $\left(-{\omega }_0^{4}y\right)$ is introduced and a homotopy equation is constructed in the next formula

$\left(D^{\left(4\right)}-{\omega }_0^4\right)y=-\rho \{P_0\left(t\right)\ddot{y}+\left(P_1\left(t\right)+{\omega }_0^4\right)y+P_3\left(t\right)y^3+P_2\left(t\right)y^2\ddot{y}+\frac{3Q}{2{\mu }^2}y{\ddot{y}}^2+\frac{3Q^2}{{\mu }^2}{y}^4\ddot{y}+\left(Q{P}_2\left(t\right)-1\right)y^5+\frac{3Q^3}{2{\mu }^2}{y}^7\}.$

(16)

Since the higher power in equation (16) is artificial, we need to reduce this power to a second-order power similar to the original one. This can be accomplished throughout integration, and then by replacing the left-hand side of equation (26) by $\left(D^{\left(2\right)}+{\omega }_0^2\right)\left(D^{\left(2\right)}-{\omega }_0^2\right)$ . Therefore, equation (16) is rewritten in the following form

$\left(D^2+{\omega }_0^2\right)y=-\frac{\rho }{D^2-{\omega }_0^2}\times \{P_0\left(t\right)\ddot{y}+\left(P_1\left(t\right)+{\omega }_0^4\right)y+P_3\left(t\right)y^3+P_2\left(t\right)y^2\ddot{y}+\frac{3Q}{2{\mu }^2}y{\ddot{y}}^2+\frac{3Q^2}{{\mu }^2}{y}^4\ddot{y}+\left(Q{P}_2\left(t\right)-1\right)y^5+\frac{3Q^3}{2{\mu }^2}{y}^7\}.$

(17)

Consider that the function $y\left(t;\rho \right)$ can be expanded in the form

$y\left(t;\rho \right)=y_0\left(t\right)+\rho y_1\left(t\right)+{\rho }^2{y}_2\left(t\right)+...$

(18)

Substituting from equation (18) into equation (17), the unknown function $y_0\left(t\right)$ is given by

$y_0\left(t\right)=A\mathit{cos}{\omega }_{0}t,$

(19)

while the function $y_1\left(t\right)$ is governed by the system of order ${\rho }^1$ as follows

$\left(D^2+{\omega }_0^2\right)y_1=-\frac{1}{D^2-{\omega }_0^2}\times \{P_0\left(t\right)\ddot{y}+\left(P_1\left(t\right)+{\omega }_0^4\right)y+P_3\left(t\right)y^3+P_2\left(t\right)y^2\ddot{y}+\frac{3Q}{2{\mu }^2}y{\ddot{y}}^2+\frac{3Q^2}{{\mu }^2}{y}^4\ddot{y}+\left(Q{P}_2\left(t\right)-1\right)y^5+\frac{3Q^3}{2{\mu }^2}{y}^7\}.$

(20)

The solution (19) is employed into equation (20) to obtain the first-order solution. The non-resonance case occurs, where the natural frequency ${\omega }_0$ and parametric frequency $\Omega$ do not interact^60,61. The removal of the secular terms, in the non-resonance case, requires

$4{\mu }^2\left(8{\mu }^2-3Q{A}^2\right){\omega }_0^2+2{\delta }^2\left(8{\mu }^2+9Q{A}^2\right)+\frac{15}{16}{Q}^2{A}^4\left(16{\mu }^2+7Q{A}^2\right)=0.$

(21)

The above equation shows that the stability of the non-resonance case requires

${\omega }_0^2>\frac{2{\delta }^2\left(8{\mu }^2+9Q{A}^2\right)+\frac{15}{16}{Q}^2{A}^4\left(16{\mu }^2+7Q{A}^2\right)}{4{\mu }^2\left(3Q{A}^2-8{\mu }^2\right)}.$

(22)

The first-order asymptotic periodic solution can be performed as

$y\left(t\right)=A\mathit{cos}{\omega }_{0}t+y_1\left(t\right),$

(23)

where the total solution $y_1\left(t\right)$ of equation (20) without secular terms is performed in the form

$y_1\left(t\right)=-a_1\left(\mathit{cos}3{\omega }_{0}t-\mathit{cos}{\omega }_{0}t\right)-a_2\left(\mathit{cos}5{\omega }_{0}t-\mathit{cos}{\omega }_{0}t\right)-a_3\left(\mathit{cos}7{\omega }_{0}t-\mathit{cos}{\omega }_{0}t\right)-a_4\left(\mathit{cos}\left({\omega }_0+\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-a_5\left(\mathit{cos}\left({\omega }_0-\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-b_1\left[{\omega }_0\mathit{sin}\left({\omega }_0+\Omega \right)t-\left({\omega }_0+\Omega \right)\mathit{sin}{\omega }_{0}t\right]-b_2\left[{\omega }_0\mathit{sin}\left({\omega }_0-\Omega \right)t-\left({\omega }_0-\Omega \right)\mathit{sin}{\omega }_{0}t\right]-a_6\left\{a_7\left(\mathit{cos}\left({\omega }_0+2\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-a_8\left(\mathit{cos}\left({\omega }_0-2\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)\right\}-b_3\left\{a_7\left({\omega }_0\mathit{sin}\left({\omega }_0+2\Omega \right)t-\left({\omega }_0+2\Omega \right)\mathit{sin}{\omega }_{0}t\right)-a_8\left({\omega }_0\mathit{sin}\left({\omega }_0-2\Omega \right)t-\left({\omega }_0-2\Omega \right)\mathit{sin}{\omega }_{0}t\right)\right\}-a_9\left(\mathit{cos}\left(3{\omega }_0+\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-a_{10}\left(\mathit{cos}\left(3{\omega }_0-\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-b_4\left({\omega }_0\mathit{sin}\left(3{\omega }_0+\Omega \right)t-\left(3{\omega }_0+\Omega \right)\mathit{sin}{\omega }_{0}t\right)-b_5\left({\omega }_0\mathit{sin}\left(3{\omega }_0-\Omega \right)t-\left(3{\omega }_0-\Omega \right)\mathit{sin}{\omega }_{0}t\right)-a_{11}\left(\mathit{cos}\left(3{\omega }_0+2\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-a_{12}\left(\mathit{cos}\left(3{\omega }_0-2\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-a_{13}\left(\mathit{cos}\left(5{\omega }_0+\Omega \right)t-\mathit{cos}{\omega }_{0}t\right)-a_{14}\left(\mathit{cos}\left(5{\omega }_0-\Omega \right)t-\mathit{cos}{\omega }_{0}t\right).$

(24)

The constant coefficients appearing in the above equation are listed below in the appendix.

Numerical estimations for the non-resonance case

To illustrate the stability picture of the non-resonance case, we estimate the stability condition (22). Besides, the periodic solution (23) should be plotted to show the behavior along with time in the non-resonant case. In Figure 4, we plot the natural frequency ${\omega }_0^2$ versus the amplitude of the parametric force for a system of $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=A=1$ . From this figure, it is clear that the system turns out to be destabilized for small values of the fractal order $\alpha$ . For some numerical data of the fractal order, damping, nonlinear, and amplitude of the periodically forced coefficients $\alpha$ , $\mu$ , $\delta$ , $Q$ are estimated and displayed in Table 1. The investigation of Table 1 shows that as $\alpha$ increases, the damping coefficient $\mu$ decreases, which means that the solution behaves as a periodic one for large values $\alpha$ . On the other side, increasing $\alpha$ leads to an increase in the amplitude of the parametric force $\delta$ . The same behavior is observed for the nonlinear coefficient. Further, the solution (23) that covers the non-resonance case is plotted in Figure 5. The analysis of the effect of the fractal order $\alpha$ is shown for a system having $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=A=\Omega =1$ . It is observed that there is no new conclusion as compared to the case of the non-parametric force. The inspection of this graph also reveals that reducing the parameter $\alpha$ leads to a reduction in the time cycle of the wave solution. Figure 6 also illustrates the impact of the coefficient of the parametric force $A$ on the periodic solution (23). The time cycle of the wave solution increases with the increase of $A$ . Further, a picture of the non-resonance periodic solution (23) profile is also displayed in Figure 7. This graph indicates that with increasing time, the influence of frequency $\Omega$ appears.

Figure 4.

Stability distribution in the plane ( ${\omega }_0^2-\delta$ ) for the condition (22) at $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=A=1$ and $\alpha =\mathrm{0.8,0.9}$ , and 0.95.

Table 1.

Some fractal order $\alpha$ values with the corresponding values of the damping, nonlinear, and amplitude of the periodically forced coefficients $\mu$ , $\delta$ , $Q$ by following equations (5) and (7), at $Q_0={\delta }_0=0.1$ and ${\widehat{\omega }}_0=1$ .

$\alpha$	$\mu$	$\delta$	$Q$
$0.1$	$0.308166$	$0.0620134$	$0.310067$
$0.2$	$0.280154$	$0.0686055$	$0.343028$
$0.3$	$0.255736$	$0.0738396$	$0.369198$
$0.4$	$0.231941$	$0.0784813$	$0.392406$
$0.6$	$0.179849$	$0.0870617$	$0.435308$
$0.7$	$0.14863$	$0.0911636$	$0.455818$
$0.8$	$0.111321$	$0.095043$	$0.475215$
$0.9$	$0.0643582$	$0.0983432$	$0.491716$

Figure 5.

The influence of the fractal order $\alpha$ on the wave solution (23) for a system having $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=A=\Omega =1$ .

Figure 6.

The variation of the coefficient of the parametric force $A$ for the system of $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=\Omega =1$ , and $\alpha =0.8$ .

Figure 7.

The time history of the periodic solution (23) for a system of $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=A=1$ , and $\alpha =0.8$ .

Sub-harmonic resonance case

In acoustics, resonance is defined as an interference modality between two sounds with somewhat different frequencies. It can also be viewed as a periodic alteration in magnitude whose rate is the distinction between the two frequencies. Three sub-harmonic resonance cases can be distinguished in the first-order problem (20), namely, $\Omega \approx 2{\omega }_0,$ $\Omega \approx 4{\omega }_0,$ and $\Omega \approx 6{\omega }_0.$ To study the stability in the sub-harmonic case $\Omega \approx 2{\omega }_0,$ as an example, we introduce a detaining parameter $\sigma$ such that

$\frac{1}{4}{\Omega }^2={\omega }_0^2+\rho \sigma ,$

(25)

where the detaining parameter will be determined later.

Considering the definition (25), and employing the expansion (18) into the homotopy equation (17), the analysis leads to the solution of the zero-order problem as

$y_0\left(t\right)=A\mathit{cos}\frac{1}{2}\Omega t.$

(26)

The first-order problem is converted to

$\left(D^2+\frac{1}{4}{\Omega }^2\right)y_1=-\frac{1}{\left(D^2-\frac{1}{4}{\Omega }^2\right)}\left\{{\stackrel{\sim }{P}}_0\left(t\right){\ddot{y}}_0+{\stackrel{\sim }{P}}_1\left(t\right)y_0+\frac{3Q}{2{\mu }^2}{y}_0{{\ddot{y}}_0}^2+{\stackrel{\sim }{P}}_3\left(t\right)y_0^3+{\stackrel{\sim }{P}}_2\left(t\right)y_0^2{\ddot{y}}_0+\frac{3Q^2}{{\mu }^2}{y}_0^4{\ddot{y}}_0+Q{\stackrel{\sim }{P}}_2\left(t\right)y_0^5+\frac{3Q^3}{2{\mu }^2}{y}_0^7\right\}$

(27)

where ${\stackrel{\sim }{P}}_i\left(t\right);i=\mathrm{0,1,2,3}$ are periodic functions defined as

$\begin{array}{l}{\stackrel{\sim }{P}}_0\left(t\right)=\frac{1}{{\mu }^2}\left[\frac{1}{2}{\Omega }^2{\mu }^2-4{\mu }^4+2\delta \Omega \mu \mathit{sin}\Omega t+4\delta {\mu }^2\mathit{cos}\Omega t\right],\\ {\stackrel{\sim }{P}}_1\left(t\right)=\frac{1}{\mu }\left[\mu \left(\frac{1}{8}{\Omega }^4+\frac{1}{2}{\Omega }^2\sigma +2{\delta }^2\right)+2\delta \Omega \left(2{\mu }^2+\frac{1}{4}{\Omega }^2\right)\mathit{sin}\Omega t-\delta \mu {\Omega }^2\mathit{cos}\Omega t+2{\delta }^2\Omega \mathit{sin}2\Omega t+2{\delta }^2\mu \mathit{cos}2\Omega t\right],\\ {\stackrel{\sim }{P}}_2\left(t\right)=\frac{3Q}{{\mu }^2}\left(\frac{1}{4}{\Omega }^2+2{\mu }^2+2\delta \mathit{cos}\Omega t\right),\\ {\stackrel{\sim }{P}}_3\left(t\right)=\frac{1}{{\mu }^2}\left[\frac{Q}{2}\left(6{\delta }^2+2{\Omega }^2{\mu }^2+\frac{3}{16}{\Omega }^4\right)+2\delta Q\mu \Omega \mathit{sin}\Omega t+2\delta Q\left(4{\mu }^2+\frac{3}{4}{\Omega }^2\right)\mathit{cos}\Omega t+3{\delta }^{2}Q\mathit{cos}2\Omega t\right].\end{array}\}.$

(28)

By utilizing the solution of (26) into the first-order equation (27), removing the secular terms requires

$Q^3{A}^6+\frac{8}{7}\left(3\delta +2{\mu }^2\right)Q^2{A}^4+\frac{16}{105}\left(21{\delta }^2+32\delta {\mu }^2-3{\Omega }^2{\mu }^2\right)Q{A}^2+\frac{64}{105}{\mu }^2\left(4{\delta }^2-2\delta {\Omega }^2+{\Omega }^2\left(\sigma +2{\mu }^2\right)\right)=0,$

(29)

and

${\Omega }^2=2\left(Q{A}^2+4{\mu }^2\right).$

(30)

Under the above conditions, the first-order solution arises in the form

$y_1=\frac{A}{256{\Omega }^2{\mu }^2}\left[\begin{array}{l}63Q^3{A}^6+24Q^2{A}^4\left(5{\mu }^2+11\delta \right)+16Q{A}^2\left(15{\delta }^2+24\delta {\mu }^2-{\Omega }^2{\mu }^2\right)\\ +128P{\mu }^2\left(\delta -{\Omega }^2\right)\end{array}\right]\left(\mathit{cos}\frac{3}{2}\Omega t-\mathit{cos}\frac{1}{2}\Omega t\right)+\frac{A\left(Q{A}^2+4\delta \right)}{768{\Omega }^2{\mu }^2}\left[21Q^2{A}^4+12Q{A}^2\left(2{\mu }^2+3\delta \right)+32\delta {\mu }^2\right]\left(\mathit{cos}\frac{5}{2}\Omega t-\mathit{cos}\frac{1}{2}\Omega t\right)+\frac{3Q{A}^3{\left(Q{A}^2+4\delta \right)}^2}{1536{\mu }^2{\Omega }^2}\left(\mathit{cos}\frac{7}{2}\Omega t-\mathit{cos}\frac{1}{2}\Omega t\right)+\frac{A\delta \Omega }{8\mu {\Omega }^2}\left(3Q{A}^2+4\delta +8{\mu }^2+{\Omega }^2\right)\left(\mathit{sin}\frac{3}{2}\Omega t-3\mathit{sin}\frac{1}{2}\Omega t\right)+\frac{A\delta \left(Q{A}^2+4\delta \right)}{24\mu {\Omega }^2}\left(\mathit{sin}\frac{5}{2}\Omega t-5\mathit{sin}\frac{1}{2}\Omega t\right).$

(31)

As mentioned before, the first-order periodic solution in the sub-harmonic resonance case $\Omega \approx 2{\omega }_0$ is carried out by using (26) and (31) into the expansion (18), letting $\rho \to 1,$ to get

$y\left(t\right)=A\mathit{cos}\frac{1}{2}\Omega t+\frac{A\left(\mathit{cos}\frac{3}{2}\Omega t-\mathit{cos}\frac{1}{2}\Omega t\right)}{256{\Omega }^2{\mu }^2}\left[\begin{array}{l}63Q^3{A}^6+24Q^2{A}^4\left(5{\mu }^2+11\delta \right)+16Q{A}^2\left(15{\delta }^2+24\delta {\mu }^2-{\Omega }^2{\mu }^2\right)\\ +128P{\mu }^2\left(\delta -{\Omega }^2\right)\end{array}\right]+\frac{A\left(Q{A}^2+4\delta \right)}{768{\Omega }^2{\mu }^2}\left[21Q^2{A}^4+12Q{A}^2\left(2{\mu }^2+3\delta \right)+32\delta {\mu }^2\right]\left(\mathit{cos}\frac{5}{2}\Omega t-\mathit{cos}\frac{1}{2}\Omega t\right)+\frac{3Q{A}^3{\left(Q{A}^2+4\delta \right)}^2}{1536{\mu }^2{\Omega }^2}\left(\mathit{cos}\frac{7}{2}\Omega t-\mathit{cos}\frac{1}{2}\Omega t\right)+\frac{A\delta \Omega }{8\mu {\Omega }^2}\left(3Q{A}^2+4\delta +8{\mu }^2+{\Omega }^2\right)\left(\mathit{sin}\frac{3}{2}\Omega t-3\mathit{sin}\frac{1}{2}\Omega t\right)+\frac{A\delta \left(Q{A}^2+4\delta \right)}{24\mu {\Omega }^2}\left(\mathit{sin}\frac{5}{2}\Omega t-5\mathit{sin}\frac{1}{2}\Omega t\right).$

(32)

Stability analysis in the sub-harmonic resonance case

It is noted that the behavior of the frequency $\Omega$ controls the stability in this case since the stability condition requires that ${\Omega }^2$ be real and have a positive value. We have a sufficient condition to satisfy the solvability conditions (30) due to Duffing’s coefficient; $Q>0.$ On the other hand, we need to solve the solvability condition (29) in the detaining parameter, which leads to

$\sigma =2\delta -{\mu }^2+\frac{3}{4}Q{A}^2-\frac{1}{64{\mu }^2{\Omega }^2}\left(105Q^3{A}^6+120Q^2{A}^4\left(3\delta +2{\mu }^2\right)+16Q{A}^2\delta \left(21\delta +32{\mu }^2\right)+256{\delta }^2{\mu }^2\right).$

(33)

To construct the frequency equation at the resonance case, we insert equation (33) into equation (25) and let $\rho \to 1,$ and the result gives

${\Omega }^4-{\stackrel{\sim }{a}}_1{\Omega }^2+{\stackrel{\sim }{a}}_0=0.$

(34)

This is a quadratic equation in the parameter ${\Omega }^2,$ where the constants ${\stackrel{\sim }{a}}_0$ and ${\stackrel{\sim }{a}}_1$ are listed below

$\begin{array}{l}{\stackrel{\sim }{a}}_1=\left(4{\omega }_0^2+8\delta -8{\mu }^2+3Q{A}^2\right),\\ {\stackrel{\sim }{a}}_0=16{\delta }^2+\frac{\delta Q{A}^2}{16{\mu }^2}\left[105Q^2{A}^4+120Q{A}^2\left(3\delta +2{\mu }^2\right)+16\left(21\delta +32{\mu }^2\right)\right].\end{array}\}$

(35)

The following conditions need to be satisfied, to ensure the presence of the positive real roots of the frequency equation (34)

${\stackrel{\sim }{a}}_1>0,{\stackrel{\sim }{a}}_0>0,\&\hspace{0.17em}{\stackrel{\sim }{a}}_1^2-4{\stackrel{\sim }{a}}_0>0.$

(36)

The above stability conditions (36) can be arranged in terms of ${\omega }_0^2$ as follows

${\omega }_0^2>2{\mu }^2-2\delta -\frac{3}{4}Q{A}^2,$

(37)

${\omega }_0^2>2{\mu }^2-2\delta -\frac{3}{4}Q{A}^2+\frac{1}{2}\sqrt{{\stackrel{\sim }{a}}_0},$

(38)

${\omega }_0^2<2{\mu }^2-2\delta -\frac{3}{4}Q{A}^2-\frac{1}{2}\sqrt{{\stackrel{\sim }{a}}_0}.$

(39)

The investigation of the above stability conditions (37)–(39) reveals that condition (39) lies between the two conditions (37) and (38). Therefore, the last two conditions are critical stability conditions. Their corresponding transition curves will surround the unstable region. It should be noted that the curves coincide when the parameter ${\omega }_0^2\to 0$ . The bifurcation will occur as the parameter ${\omega }_0^2$ divers from zero. It is important to note that the increase in the positive value of the nonlinear coefficient $Q$ will increase this bifurcation and, therefore, increase the instability behavior of the system.

Numerical estimations at the harmonic resonance case

The small detaining parameter $\sigma$ is defined in equation (25) and given by equation (33), measures the vibration of the resonance case, and is plotted in Figure 8. It is calculated for a system with a fixed value of $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=\Omega =1$ , and the subsequent values of the fractal order $\alpha =\mathrm{0.95,0.9,0.8,0.7}$ are collected together in this graph. Investigation of Figure 8 reveals that the increase in $\alpha$ leads to increasing the small detaining parameter $\sigma$ . The mechanism for decreasing $\sigma$ is faster for very small values in the amplitude $A$ . However, the vibration that occurs due to the resonance case will be suppressed by the decrease in the fractal order parameter $\alpha$ . Additionally, the numerical illustrations of the periodic solution in the sub-harmonic resonance that is given in (32) are plotted in Figures 9 and 10. The influence of the parametric force $A$ on the periodic solution (32) can be found in the graph of Figure 9. It is noticed that the increase in $A$ yields a greater boost in the periodic solution value. Finally, it is noticeable that the wave solution starts when $t=0$ at the point. After a slight increase in the variable $t$ , the amplitude has a relatively larger value. Small decreases in the fractal order $\alpha$ reduce the amplitude $A$ as shown in Figure 10. A slight rise in the coefficient of the parametric force leads to a greater rise in the value of the amplitude $A$ as observed in Figure 10.

Figure 8.

The change of the detaining parameter $\sigma$ versus the fractal-order $\alpha \hspace{0.17em}\mathrm{a}\mathrm{t}\hspace{0.17em}{Q}_0={\delta }_0=0.1,{\widehat{\omega }}_0=\mathrm{\Omega }=1$ .

Figure 9.

The change of the solution (32) with the parametric force $A$ for the system of $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=1$ , and $\alpha =0.6$ .

Figure 10.

The change of the solution (32) with the fractal-order $\alpha$ at $Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=A=1$ .

Advanced task for the analytical solution of the parametric fourth-order equation

In this section, another analytical periodic solution can be performed for equation (13) without using the perturbation method and without reducing its rank. This task can be achieved using the non-perturbative technique^59,62–65. As mentioned before, this system of a nonlinear oscillator having periodic coefficients represents a highly complicated equation. The plan for handling this equation has two steps. The first step deals with obtaining a linear representation that corresponds to the nonlinear equation (13). This aim can be accomplished using Caughey’s linearization method ⁶⁶ as proposed by El-Dib. ⁶³ The second step is concerned with calculating the solution of the parametric linear fourth order having multi-periodic coefficients. To provide the required solution, Galerkin’s method can be used to obtain the frequency for the parametric differential equation. The procedure is as follows:

It is noted that equation (13) can be viewed as

$y^{\left(4\right)}+g\left(y;t\right)\ddot{y}+f\left(y;t\right)y=0,$

(40)

where the parametric nonlinear coefficients $g\left(y;t\right)$ and $f\left(y;t\right)$ are listed below

$g\left(y;t\right)=P_0\left(t\right)+P_2\left(t\right)y^2+\frac{3Q}{{\mu }^2}{y}^4,$

(41)

$f\left(y;t\right)=P_1\left(t\right)+P_3\left(t\right)y^2+\left(Q{P}_2\left(t\right)-1\right)y^4+\frac{3Q^3}{2{\mu }^2}{y}^6+\frac{3Q}{2{\mu }^2}{\ddot{y}}^2.$

(42)

It is noted that the linearized form of equation (40) can be established by using the equivalent linearization technique given in Ref 66. Since both the functions $g\left(y;t\right)$ and $f\left(y;t\right)$ are even functions in the variable $y$ , then they are required to be rearranged in the form

$\widehat{g}\left(y;t\right)=\frac{P_0\left(t\right)y+P_2\left(t\right)y^3+\frac{3Q}{{\mu }^2}{y}^5}{y},\mathrm{a}\mathrm{n}\mathrm{d}$

(43)

$\widehat{f}\left(y;t\right)=\frac{P_1\left(t\right)y+P_3\left(t\right)y^3+Q{P}_2\left(t\right)y^5+\frac{3Q^3}{2{\mu }^2}{y}^7+\frac{3Q}{{\mu }^2}{\ddot{y}}^{2}y}{y}.$

(44)

Simply, the linearization will be available when equation (40) can be written as

$y^{\left(4\right)}+\beta \left(t\right)\ddot{y}+{\varpi }^4\left(t\right)y=0,$

(45)

where the parametric coefficients $\beta \left(t\right)$ are estimated by the simulation of the modification to He’s frequency formula proposed by El-Dib ⁶³ as

$\beta \left(t\right)={\frac{D_y\left(P_0\left(t\right)y+P_2\left(t\right)y^3\right)+\frac{3Q}{{\mu }^2}\frac{y{D}_y^2}{2!}\left(y^5\right)}{D_y\left(y\right)}|}_{y=\frac{1}{2}A}=P_0\left(t\right)+\frac{3}{4}{P}_2\left(t\right)A^2+\frac{15}{8}Q{A}^4,$

(46)

and

${\varpi }^4\left(t\right)={\frac{D_y\left(P_1\left(t\right)y+P_3\left(t\right)y^3\right)+\frac{\left(Q{P}_2\left(t\right)-1\right)y{D}_y^2{y}^5}{2!}+\frac{3Q^3}{2{\mu }^2}\frac{y^2}{3!}{D}_y^3{y}^7+\frac{3Q}{2{\mu }^2}{D}_y{\ddot{y}}^{2}y}{D_{y}y}|}_{y=\frac{1}{2}A,\ddot{y}=\frac{1}{2}A{\varpi }^2\left(t\right)}$

This leads to

${\varpi }^4\left(t\right)=\frac{8{\mu }^2}{8{\mu }^2-3Q{A}^2}\left(P_1\left(t\right)+\frac{3}{4}{P}_3\left(t\right)A^2+\frac{5}{8}\left(Q{P}_2\left(t\right)-1\right)A^4+\frac{105}{128}\frac{Q^3}{{\mu }^2}{A}^6\right),$

(47)

where $\omega$ is assumed to be the oscillation frequency. In terms of (46) and (47), equation (45) will represent a linear fourth-order equation having a parametric coefficient. To obtain its periodic solution using a non-perturbative approach, we may introduce a trial solution in the form

$\widehat{y}\left(t\right)=A\mathit{cos}\omega t.$

(48)

According to the suggested solution (48), equation (45) may have the form

$y^{\left(4\right)}-{\omega }^4\left(\Omega \right)y=0,$

(49)

where ${\omega }_0^4\left(\Omega \right)$ can be estimated using Galerkin’s method. The corresponding residual can be established by inserting (48) into equation (45), which yields $R\left(\omega \left(\Omega \right);t\right)=0$ . Employing the suggested solution (48) into equation (45) yields

$R\left({\omega }_0\left(\Omega \right);t\right)=\left({\omega }^4\left(\Omega \right)-\beta \left(t\right){\omega }^2\left(\Omega \right)+{\varpi }^4\left(t\right)\right)\widehat{y}\left(t\right)=0.$

(50)

Applying Galerkin’s method to (50) requires

$\underset{0}{\overset{T}{{\displaystyle \int }}}R\left(\omega \left(\Omega \right);t\right)\mathit{cos}\left(\omega t\right)dt=0;T=\frac{\pi }{2\omega }.$

(51)

To evaluate the above integration, the normalized form of the excitation frequency $\Omega$ may be established, so the frequency $\Omega$ of equation (50) may be replaced by $\left(\widehat{\Omega }/T\right)$ . Doing the integration in equation (51), the frequency equation results in the form

$A_4{\omega }^4+A_3{\omega }^3+A_2{\omega }^2+A_1\omega +A_0=0,$

(52)

where the coefficients $A^{\prime }s$ are listed below

$\begin{array}{l}{A}_4=16\pi \widehat{\Omega }\left({\pi }^4-5\pi {\widehat{\Omega }}^2+4{\widehat{\Omega }}^4\right)\left(3Q{A}^2+8{\mu }^2\right),\\ A_3=512\delta \mu \widehat{\Omega }\left({\pi }^2-4{\widehat{\Omega }}^2\right)\left[{\pi }^2\mathit{cos}\widehat{\Omega }-\left({\pi }^2-2{\widehat{\Omega }}^2\right)\right],\\ A_2=-16\widehat{\Omega }\pi \left({\pi }^4-5{\pi }^2{\widehat{\Omega }}^2+4{\widehat{\Omega }}^4\right)\left[15A^{4}Q+16{\mu }^2\left({\omega }_0^2-2{\mu }^2\right)+18A^{2}Q\left({\omega }_0^2+2{\mu }^2\right)\right]\\ -64\pi \delta \left({\pi }^2-4{\widehat{\Omega }}^2\right)\left[{\pi }^2\left(9A^{2}Q+8{\mu }^2\right)+16{\mu }^2{\widehat{\Omega }}^2\right]\mathit{sin}\widehat{\Omega },\\ A_1=128\delta \mu \widehat{\Omega }\left[\left({\pi }^4+8{\widehat{\Omega }}^4\right)\left(2{\omega }_0^2+\delta +4{\mu }^2\right)-6{\pi }^2{\widehat{\Omega }}^2\left(4{\omega }_0^2+3\delta +8{\mu }^2+3A^{2}Q\right)\right]\\ -128{\pi }^2\delta \mu \widehat{\Omega }\left(4{\omega }_0^2+8{\mu }^2+3A^{2}Q\right)\left({\pi }^2-4{\widehat{\Omega }}^2\right)\mathit{cos}\widehat{\Omega }-256{\pi }^2{\delta }^2\mu \widehat{\Omega }\left({\pi }^2-{\widehat{\Omega }}^2\right)\mathit{cos}2\widehat{\Omega },\\ A_0=\pi \widehat{\Omega }\left({\pi }^2-4{\widehat{\Omega }}^2\right)\left({\pi }^2-{\widehat{\Omega }}^2\right)[\left(105A^6{Q}^3+128({\omega }_0^4+2{\delta }^2\right){\mu }^2+240A^4{Q}^2\left({\omega }_0^2+{\mu }^2\right)\\ +48A^{2}Q\left(3{\omega }_0^4+8{\omega }_0^2{\mu }^2+6{\delta }^2\right)]+32{\pi }^3\delta \left({\pi }^2-4{\widehat{\Omega }}^2\right)\left[15A^4{Q}^2+16{\omega }_0^2{\mu }^2+6A^{2}Q\left(3{\omega }_0^2+4{\mu }^2\right)\right]\mathit{sin}\widehat{\Omega }\\ +16{\pi }^3{\delta }^2\left(9A^{2}Q+8{\mu }^2\right)\left({\pi }^2-{\widehat{\Omega }}^2\right)\mathit{sin}2\widehat{\Omega }.\end{array}\}.$

(53)

With the real roots of equation (52), the solution of equation (40) which simulates the suggested solution (48) is performed in the form

$y\left(t\right)=A\mathit{cos}\omega t.$

(54)

This solution represents a solution of the parametric fourth-order nonlinear equation (13). Rather, it represents a solution to the original fractal Mathieu equation (1). The comparison of the above analytical solution (54) with the numerical solution of equation (1) can be accomplished in the simulation as the fractal parameter $\alpha$ is very closer to unity, which proves the reliability of the method. As $\alpha \to 1$ , the parameter $\mu \to 0$ . Consequently, the frequency equation (52) is reduced to the form

$A_4{\widehat{\omega }}^4+A_2{\widehat{\omega }}^2+{\widehat{A}}_0=0,$

(55)

where ${\widehat{A}}_0=\underset{\mu \to 0}{\mathit{lim}}{A}_0$ and $\widehat{\omega }=\underset{\mu \to 0}{\mathit{lim}}\omega$ . Therefore, the frequency that controls the case $\alpha =1$ is given as

$\widehat{\omega }=\sqrt{-\frac{A_2}{2A_4}+\frac{1}{2A_4}\sqrt{{A_2}^2-4A_4{\widehat{A}}_0}}.$

(56)

This means that the standard nonlinear Mathieu equation takes the form

$\frac{d^{2}y}{d{t}^2}+\left({\varpi }^2+2{\delta }_0\mathit{cos}\Omega t\right)y+Q_0{y}^3=0.$

(57)

The comparison between the analytical solution (54) and the numerical solution of equation (57) appears in Figure 11. Finally, the analytical solution (54) for variation of the fractal order with a comparison with the numerical solution of equation (57) for the same system considered in Figure 11 has been represented in Figure 12. It is clear that with an increase in the value of the fractal parameter $\alpha$ , the oscillation increases. Also, the numerical results for the effect of the fractal parameter $\alpha$ show that the HPM solution in Figure 5 has the same behavior as the corresponding non-perturbative technique solution in Figure 12. It turns out that decreasing the value of the parameter $\alpha$ leads to a contraction in the wavelength of the solution.

Figure 11.

Comparison of the analytical solution (54) with the numerical solution of the original Mathieu equation (57) for a system having $\Omega =1,Q_0={\delta }_0=0.1,{\widehat{\omega }}_0=1,A=0.5$ , and $\alpha =1$ .

Figure 12.

Illustration of the analytical solution (54) for variation of the fractal order $\alpha$ with the comparison with the numerical solution of equation (57) for the same system considered in Figure 11.

Conclusions

In this article, the periodic solution to the fractal Mathieu–Duffing equation is discovered. The steps listed below have been used to do this.

1. The fractal oscillator in the fractal space has been converted into a damping Mathieu–Duffing one in the continuous space by employing a new modification of He’s definition of the fractal derivative.

2. Utilizing a novel method named the RUT can create a similar equation, aside from any challenges in managing the impact of the linear damping component.

3. Ensuring that the numerical solutions of the alternative equation agree with the solutions of the original equation.

Thus, a similar equation can be manipulated to obtain the desired solution. Two different techniques have been suggested. The homotopy perturbation technique has been used to find the required periodic solution without sacrificing the damping coefficient. The stability behavior in the non-resonance case as well as in the sub-harmonic resonance case has also been discussed. Further, another method called “the non-perturbative procedure” that deals with the obtained equation has been introduced. The current technique is very important in the study of electric resonators and their applications in engineering. Moreover, the current procedure is suitable to examine various oscillations and acoustic waves in plasma physics systems.