New analytical solution of the fractal anharmonic oscillator using an ancient Chinese algorithm: Investigating how plasma frequency changes with fractal parameter values

Abstract

This paper uses the two-scale fractal dimension transform and He’s formula derived from the ancient Chinese algorithm Ying Bu Zu Shu to find the approximate frequency–amplitude expression of the fractal and forced anharmonic oscillator that can be used to study the nonlinear oscillations produced by the plasma physics fractal structures. The results show how the electron frequency and wavelength change as a function of the plasma physics fractal structure. In fact, if the value of the fractal parameter is decreased, the wavelength increases, and consequently, the system frequency decreases. The introduced solution procedure sheds a bright light on the easy-to-follow steps to obtain an accurate steady-state analytical solution of fractal anharmonic nonlinear oscillators.

Introduction

Fractal calculus application to model several physical and engineering phenomena has gained popularity since introducing the two-scale fractal transform that unveils hidden phenomena beyond the conventional continuum mechanics.^1–6 The resulting equations of motion are usually described by nonlinear ordinary differential equations whose frequency-amplitude relationships are determined using, for instance, homotopy perturbation methods,^7–12 parameter expansion method (PEM), ¹³ He’s frequency formulation,^14–16 iteration perturbation method,^17,18 the fractal variational principle,^19–23 energy balance method, ²⁴ energy method, ²⁵ He’s frequency formulation based on energy conservation, ²⁶ frequency-amplitude’s formula for non-conservative oscillators, ²⁷ the enhanced homotopy perturbation method, ²⁸ Taylor series,^29,30 the equivalent power-form representation, ³¹ and an ancient Chinese algorithm ³² to name a few.

The use of small and large scales to describe phenomena that occur in nature has attracted the attention of several researchers since Ji-Huan He introduced the two-scale dimension transform. To illustrate its applicability, let us consider lotus leaves when water drops fall on their surface. Studies show super-hydrophobic properties of the lotus leaves’ surface due to their rough texture visible for the human eye (macro scale). However, when observing the surface morphology at the nano/micro scale using scanning electron microscope images (SEM), countless protrusions are evident on the leaf surface that can change the super-hydrophobic surface property to hydrophilic when interacting with water.^33,34 This physical behavior of the water drops on the lotus leaf surface can be captured if the mathematical model considers small- and large-scale dimensions as studied in Refs. 35–38.

Furthermore, this two-scale transform has been applied to study the mechanism of snow’s insulation properties, ³⁹ to investigate one-dimensional microgravity flows, ²⁰ to find the frequency-amplitude expression of nonlinear oscillators in fractal space,^40,41 to study convection-diffusion processes, ⁴² for electrochemical applications such as sensors, ⁴³ to model solvent evaporation during fabrication of porous fibers by electrospinning,^29,30,44,45 in porous electrodes, ⁴⁶ lasers, ³¹ and on water collector from air^47–49 to name a few.

Therefore, this article focuses on deriving the primary resonance frequency-amplitude expression of the following anharmonic oscillator (also known as the fractal Helmhotz-Duffing equation)

$$\frac{d}{d{S}^{\alpha }}\left(\frac{dy}{d{S}^{\alpha }}\right)+{\alpha }_{0}y+{\alpha }_{22}{y}^2+{\alpha }_1{y}^3=Q\cos {\omega }_{f}S,\hspace{1em}y\left(0^{\alpha }\right)=A,\hspace{1em}\frac{dy\left(0^{\alpha }\right)}{dS}=0$$

(1)

that can be used to investigate how the resonance frequency and wavelength change as a function of the fractal parameter values of α. Here, y represents the oscillations displacement; S is the time; α is the fractal dimensions of time; α₀ represents the system natural frequency; α₁ and α₂₂ are system parameters; Q and ω _f are the external excitation force and frequency, respectively; and A is the initial displacement.

Experimental evidence suggests that equation (1) can be used to study the nonlinear oscillation produced by the plasma physics fractal structures that arise when a metal target is hit with a sufficiently intense laser pulse producing aerosol particles, which are sources of fractal aggregates from laser plasma.^50–55 In other words, equation (1) can be used to characterize the dynamics phenomena that occur when high-density plasma interacts with high-frequency electromagnetic waves.^56–59

We shall next use the two-scale fractal transform to write equation (1) in equivalent form.

The two-scale fractal transform

In order to derive the analytical fractal frequency-amplitude relationship of equation (1) using an ancient Chinese mathematical algorithm, the fractal derivative definition

$$\frac{dy}{d{S}^{\alpha }}\left(S_0^{\alpha }\right)=\mathrm{\Gamma }\left(1+\alpha \right){\text{lim}}_{\frac{\left(S-S_0\right)\to \mathrm{\Delta }t}{\mathrm{\Delta }S\ne 0}}\frac{y\left(S\right)-y\left(S_0\right)}{{\left(S-S_0\right)}^{\alpha }}$$

(2)

and the two-scale fractal dimension transform

$$t=S^{\alpha }$$

(3)

are introduced.^1,6 Thus, equation (1) can be written as a nonlinear differential equation of the form

$$\frac{d^{2}y}{d{t}^2}+{\alpha }_{0}y+{\alpha }_{22}{y}^2+{\alpha }_1{y}^3=Q\cos {\omega }_f{t}^{1/\alpha },y\left(0^{1/\alpha }\right)=A,\hspace{1em}\frac{d}{d{t}^{1/\alpha }}y\left(0^{1/\alpha }\right)=0$$

(4)

Equation (4) models an anharmonic oscillator that describes nonlinear phenomena that appear in physics of plasma,^53,59 thin laminated plates, ⁶⁰ acoustics, ⁶¹ naval engineering, ⁶² eardrum oscillations, ⁶³ elasto-magnetic suspensions, ⁶⁴ graded beams, ⁶⁵ asymmetric oscillators, ⁶⁶ and electronics ⁶⁷ to name a few. Notice that equation (4) has an exact solution based on Jacobi elliptic functions when and only when Q = 0. ⁶⁸ For Q ≠ 0, the exact solution of equation (4) is unknown; therefore, in order to determine the angular frequency of equation (4), we use the ancient Chinese algorithm Ying Bu Zu Shu published from an ancient Chinese mathematics monograph—The Nine Chapters on the Mathematical Art. ⁶⁹ In this ancient Chinese method, the approximate solutions for equation (4) are assumed to be of the form y₁(t) and y₂(t). Substitution of y₁(t) and y₂(t) into equation (4) leads to the residuals R₁(t) and R₂(t) which can be used to obtain average trail residuals ${\tilde{R}}_1\left(t\right)$ and ${\tilde{R}}_2\left(t\right)$ needed to find the frequency-amplitude expression from Ji-Huan He’s formula. Thus, in accordance with the ancient Chinese algorithm,^69–73 the steady-state trial residual functions y₁ and y₂ are assumed to be of the form

$$y_i=A_i\cos \left({\omega }_i{t}^{1/\alpha }\right),i=\mathrm{1,2}\left(\text{no}\hspace{0.17em}\text{sum}\right)$$

(5)

where A _i and ω _i are, respectively, oscillation amplitudes and frequencies that need to be determined. Thus, the substitution of equation (5) into equation (4) yields the following trial residual functions R _i

$$R_i=\frac{d^2{y}_i}{d{t}^2}+{\alpha }_0{y}_i+{\alpha }_{22}{y}_i^2+{\alpha }_1{y}_i^3-Q\cos {\omega }_i{t}^{1/\alpha }$$

(6)

It is also assumed that ω₁ = ω _f and ω ₂ = ω _f .^69–75 Expanding equation (6) gives

$$\begin{array}{lll}{R}_i\hfill & =\hfill & A_i^2{\alpha }_{22}\cos {\left[t^{1/\alpha }{\omega }_i\right]}^2+1/2\cos \left[t^{1/\alpha }{\omega }_i\right](2A_i{\alpha }_0-2Q+\hfill \\ \hfill & \hfill & A_i^3{\alpha }_1-2A_i{t}^{2/\alpha -2}{\omega }_i^2/{\alpha }^2+A_i^3{\alpha }_1\cos \left[2t^{1/\alpha }{\omega }_i\right])+\hfill \\ \hfill & \hfill & A_i\left(\alpha -1\right)t^{1/\alpha -2}{\omega }_i\sin \left[t^{1/\alpha }{\omega }_i\right]/{\alpha }^2\hfill \end{array}$$

(7)

The next step in the ancient Chinese algorithm consists in calculating the average trail residuals ${\tilde{R}}_1$ and ${\tilde{R}}_2$ defined as

$${\tilde{R}}_i=\frac{4}{T}{\displaystyle {\int }_0^{T/4}{R}_i}\left(t\right)dt$$

(8)

where T = 2π/ω _i .^69–75 To simplify the computation of the integral (8) first, the transformation S = t^1/α is introduced. Hence, equations (7) and (8) become

$$\begin{array}{lll}{R}_i\hfill & =\hfill & 1/2\cos \left[S{\omega }_i\right](2A_i{\alpha }_0-2Q+A_i^3{\alpha }_1+2A_i^2{\alpha }_{22}\cos \left[S{\omega }_i\right]+\hfill \\ \hfill & \hfill & A_i^3{\alpha }_1\cos \left[2S{\omega }_i\right])-\left(A_i{S}^{1-2\alpha }{\omega }_i\left(S{\omega }_i\cos \left[S{\omega }_i\right]-\left(\alpha -1\right)\sin \left[S{\omega }_i\right]\right)\right)/{\alpha }^2\hfill \end{array}$$

(9)

and

$${\tilde{R}}_i=\frac{4\alpha }{T}{\displaystyle {\int }_0^{{\left(T/4\right)}^{1/\alpha }}{R}_i}{S}^{\left(\alpha -1\right)}dS$$

(10)

Substitution of equation (9) into equation (10), and evaluating the integrals, we obtain the following expression for ${\tilde{R}}_i$

$$\begin{array}{l}{\tilde{R}}_i=1/T_i^2{4}^{-\left(1+2/\alpha \right)}(1/\left(\left(\alpha -2\right)\alpha \right)64A_i{T}_i^{2/\alpha }{\omega }_{ip}^2{F}_q\left[\left\{1-\alpha /2\right\},\left\{1/2,2-\alpha /2\right\},-4^{-1-2/\alpha }{T}_i^{2/\alpha }{\omega }_i^2\right]\\ -\left(1/\left(\left(\alpha -2\right)\alpha \right)\right)64A_i{T}_i^{2/\alpha }\left(\alpha -1\right){\omega }_{ip}^2{F}_q[\left\{1-\alpha /2\right\},\left\{3/2,2-\alpha /2\right\},-4^{-1-2/\alpha }{T}_i^{2/\alpha }{\omega }_i^2]\\ +{16}^{1+1/\alpha }{4}^{-2}{T}_i^2(A_i^3{\alpha }_{1p}{F}_q[\left\{\alpha /2\right\},\left\{1/2,1+\alpha /2\right\},\left(-9\right)4^{-1-2/\alpha }{T}_i^{2/\alpha }{\omega }_i^2]\\ +{\left(4A_i{\alpha }_0-4Q+3A_i^3{\alpha }_1\right)}_p{F}_q[\left\{\alpha /2\right\},\left\{1/2,1+\alpha /2\right\},-4^{-1-2/\alpha }{T}_i^{2/\alpha }{\omega }_i^2]\\ +2A_i^2{\alpha }_{22}\left(1{+}_p{F}_q\left[\left\{\alpha /2\right\},\left\{1/2,1+\alpha /2\right\},-{16}^{-1/\alpha }{T}_i^{2/\alpha }{\omega }_i^2\right]\right))\end{array}$$

(11)

where _p F _q is the generalized hypergeometric functions. ⁷⁶

The final step in the ancient Chinese algorithm consists in using He’s formula

$${\omega }_{AE}^2=\frac{{\omega }_1^2{\tilde{R}}_2-{\omega }_2^2{\tilde{R}}_1}{{\tilde{R}}_2-{\tilde{R}}_1}$$

(12)

to obtain an analytical fractal frequency-amplitude expression. As usual, ω₁ = 1 and ω₂ = 2.^69–75 Considering the nature of the Helmholtz-Duffing anharmonic oscillator, it is further assumed that A₁ = A and A₂ = βA where β is determined fitting the approximate backbone curve

$${\omega }_{AE}=\sqrt{\frac{12\left(4+\beta \left({\alpha }_0-4\right)-4{\alpha }_0\right)+8A^2\left({\beta }^3-4\right){\alpha }_1+3A\pi \left({\beta }^2-4\right){\alpha }_{22}}{12\left(1+\beta \left({\alpha }_0-4\right)-{\alpha }_0\right)+8A^2\left({\beta }^3-1\right){\alpha }_1+3A\pi \left({\beta }^2-1\right){\alpha }_{22}}}$$

(13)

derived from (12) with α = 1 and Q = 0, with the exact solution of the homogeneous equation (4) found in Ref. 68.

Results

This section elucidates the accuracy attained using the analytical expression of ω _AE derived using the ancient Chinese algorithm when compared to the exact solution (solid blue line) of equation (4) when Q = 0. ⁶⁸

Figure 1 shows the backbone curve obtained from the exact solution of equation (4) (solid blue line) ⁶⁸ and the one computed from equation (13) (dashed red line) for the system parameter values of α₀ = 1.014 82, α₁ = 0.02, and α₂₂ = 0.029 8 with β = 0.88. It is observed from Figure 1 that the theoretical backbone curve computed from equation (13) closely follows the exact backbone curve, which indicates the accuracy provided by our derived expression (13). Figure 2 illustrates the amplitude-time response curves obtained by using the exact solution of equation (4) derived in Ref. 68 and those computed from equation (13) considering the fractal parameter values of α = 0.9, 1, and 1.1. One can see from Figure 2 that the maximum oscillation amplitude is the same for all values of α; however, when α < 1 or α > 1, the oscillation frequency becomes lower or higher than that of α = 1, respectively. From Figures 1 and 2, one can conclude that the approximate frequency-amplitude and amplitude-time relationships derived using the ancient Chinese algorithm describe well the system dynamics when Q = 0.

Figure 1.

Frequency-amplitude backbone curve of the anharmonic equation for which α₀ = 1.014 8, α₁ = 0.02, α₂₂ = 0.0298, and β = 0.88. Notice that the backbone curve obtained from equation (13) closely follows the exact backbone curve computed from equation (4).

Figure 2.

Amplitude-time response curves for system parameter values of α₀ = 1.0148, α₁ = 0.02, α₂₂ = 0.0298, Q = 0, and β = 0.88. Here, the blue dashed lines describe the exact numerical solution of equation (4) obtained by using the solution derived in Ref. 68., while the color solid lines are the approximate solutions computed from equations (5) and (12).

We next plot the approximate frequency-amplitude response curves using equation (11) when an electric (restoring) force |Q| = 0.1 is acting with different fractal order values of α = 0.9, 1, and 1.1 with β = 0.88. One can see from Figure 3 that the frequency-amplitude curves shifted to the right of the curves computed with α = 1 when the fractal parameter values are bigger than 1 (α > 1), or to the left when the values of α are less than one (α < 1). Also, notice that an increase in α gives rise to a decrease in the wavelength since the oscillation frequency increases. Furthermore, if the value of the fractal parameter α is decreased, the wavelength increases, and consequently, the electron frequency decreases. ⁷⁷ Figure 4 illustrates the amplitude-time curves computed using equations (5) and (13) and those obtained from equation (4) using the fourth-order RungeKutta method provided by the symbolic program of Mathematica. Notice that our analytical solution predicts the qualitative and quantitative fractal dynamic behavior of the anharmonic oscillator well.

Figure 3.

Frequency-amplitude curve of the anharmonic fractal oscillator computed from equation (12) using the values of α₀ = 1.0148; α₁ = 0.02; α₂₂ = 0.0298; Q = 0.1; β = 0.88; and α = 0.9, 1, and 1.1. Notice that the fractal frequency-amplitude curves shifted to the right of the curves computed with α = 1 when the fractal dimension values are bigger than 1 (α > 1), or to the left when the values of α are less than one (α < 1).

Figure 4.

Amplitude-time response curves for system parameter values of α₀ = 1.0148; α₁ = 0.02; α₂₂ = 0.0298; Q = −0.1; β = 0.88; and α = 0.9, 1, and 1.1. Here, the blue dashed lines describe the exact numerical solution of equation (4) obtained using the fourth-order RungeKutta method provided by the symbolic program of Mathematica, while the color solid lines are the approximate solutions computed from equations (5) and (12).

Finally, when α₁ = 0 and α₂₂ = 0 in equation (1), it becomes the classical equation used to study fractional electromagnetic waves in plasma physics. ⁵⁹ In this case, the exact fractal frequency–amplitude response expression was determined by Elías-Zúñiga et al. in Ref. 78 using the Vakakis and Blanchard approach jointly with the two-scale fractal dimension transform. Following this approach, these researches avoided the complexity associated with fractional calculus providing a solution that describes physical observations well.

Conclusion

In this article, we have used the two-scale fractal calculus and He’s formula, which was developed from the ancient Chinese algorithm to derive an approximate analytical frequency-amplitude expression for the fractal Helmhotz-Duffing nonlinear differential equation that describes the nonlinear oscillations that occur in plasma physics. This analytical solution provides the fractal frequency-amplitude response curves from which it is possible to determine how the frequency and wavelength change as a function of the fractal dimension. Furthermore, the results show how the electron frequency and wavelength change as a function of the plasma physics fractal structures that arise when a metal target is hit with a sufficiently intense laser pulse. An increase in α decreases the wavelength with an increase in the electron frequency and vice versa.

The results obtained in this paper elucidate the applicability of the two-scale fractal dimension transform that unveils the critical role that the fractal parameter value α has in shifting to the left or to the right, of the backbone curve, the fractal frequency-amplitude response curves of anharmonic oscillators that model electromagnetic waves in plasma physics.