Exact Parametric and Semi-Analytical Solutions for the Rucklidge-Type Dynamical System

Abstract

The behavior of the Rucklidge-type dynamical system was investigated, providing some semi-analytical solutions, in this paper. This system was analytically investigated by means of the Optimal Auxiliary Functions Method (OAFM) for two cases. An exact parametric solution was obtained. The effect of the physical parameters was investigated on the asymptotic behaviors and damped oscillations of the solutions. Damped oscillations are essential for analyzing and designing various mechanical, biological, and electrical systems. Many of the applications involving these systems represent the main reason of this work. A comparison between the obtained results via the OAFM, the analytical solution obtained with the iterative method, and the corresponding numerical solution was performed. The accuracy of the analytical and corresponding numerical results is illustrated by graphical and tabular representations.

1. Introduction

In nature, there are many phenomena that involve the study of dynamical systems. The behavior of complex dynamical systems is described by the solutions of the equations of motion generated by numerical analysis.

The Rucklidge chaotic system was introduced in 1992 [1]. It exhibits extremely rich dynamical behavior and has diverse applications [1,2,3,4,5,6,7,8,9]. The model is used in the case of parameter regimes where chaotic solutions occur in a three-dimensional space. The numerical solution of the Rucklidge model can cause problems because the model is chaotic. The Rucklidge model is a family of three-dimensional systems of quadratic differential equations. Quadratic systems in ${\mathbb{R}}^3$ are the simplest systems after linear ones. The Rucklidge chaotic system is a dissipative system that tends to dissipate energy over time [4]. The dissipative nature of the system can also be seen in the dynamical behaviors of this model. The system exhibits a high sensitivity to initial conditions, which means that small changes in the initial configuration of the system can generate very different trajectories of the system. Thus, the basic property of such models is the extreme sensitivity of the solution to numerical inaccuracies and initial conditions [4]. Small changes in the initial values can lead to a large differences in time, and these are highlighted in certain dynamical properties of chaos, including their equilibria, stability, Lyapunov exponents, and the bifurcation diagram [2,3,5]. In the Rucklidge system, it is the physical parameters that affect the dynamical behaviors. For different values of these parameters, the model expresses different bifurcations. Bifurcation is one of the important tools that shows whether a map or a system is chaotic [5]. In the literature, numerous theoretical and numerical results investigated using the Rucklidge system have been given. Thus, this Rucklidge system was used in the fluid mechanics in the case of a Boussinesq fluid for the study of 2D convection in a horizontal layer [4,6]. It is known in the literature to describe and analyze the unstable periodic orbits of the Rucklidge system in order to explore the hidden topological properties in periodic orbits using a symbolic coding method. In [7], bifurcations of periodic orbits were explored, allowing for improved dynamics of the Rucklidge system. Unstable periodic orbits of turbulent flows have attracted considerable attention in fluid dynamics research as they provide the building blocks of disordered dynamics [7]. In a chaotic system, unstable periodic orbits also play important roles in the analysis of dynamical behavior. In [5], a three-dimensional discrete fractional-order Rucklidge system with complex state variables was studied. The dynamic nature and chaotic behavior exhibited by the Rucklidge system with real state variables and the higher-dimensional system derived from complex state variables were compared using various methods, such as bifurcation analysis and maximum Lyapunov exponents by the Jacobian matrix method. In [8], the modeling of the Rucklidge system and nonlinear fractional differential equations was carried out using variable-order differential operators using fractional calculus and Newton polynomial interpolation. In [7], unstable periodic orbits in the Rucklidge system were investigated by a variational method up to a certain topological length. The dynamics of the Rucklidge system were explored using Lyapunov exponents, phase portrait analysis, and Poincaré first-return maps. A chaotic oscillator was designed and simulated, and the synchronization and masking communication circuits of the Rucklidge attractor were studied [7]. Small-amplitude limit cycles in the Rucklidge system were considered, and the bifurcation analysis at equilibria was performed by calculating the second Lyapunov coefficient or the first Lyapunov value at the stability limit. Various techniques were used to control the chaotic behavior of the Rucklidge system, from which the constitutive elements of disordered dynamics, useful in turbulent motion in fluid mechanics, were found. The Rucklidge system is a model of a double convection process, which models convection in an imposed vertical magnetic field and in a uniformly rotating fluid layer. In [10], by analytical method and using the method of undetermined coefficients, the existence of heteroclinic and homoclinic orbits in the Rucklidge system was demonstrated and the explicit and uniformly convergent algebraic expressions of Si’lnikov-type orbits were given.

In [9], two widely used approaches to chaos synchronization were investigated: active control and backstepping control. Different initial conditions were considered for two identical chaotic systems (Master/Drive of Rucklidge Systems). In the concept of synchronization of two similar nonlinear chaotic systems, it is assumed that different initial conditions are initially assumed and that the Master system can follow the trajectories of the Drive system. An appropriate control law is applied, and it is assumed that the two chaotic systems cannot synchronize with each other. However, if the two systems share information regularly, they will be able to synchronize. For the synchronization of chaos, a multitude of methods can be used: the linear error feedback method, the delayed feedback method, the active control approach, the impulsive method, the backstepping approach, and other control methods. These methods have been applied to many practical systems [9], such as the Rikitake two-disk dynamo, Chua’s circuits, Van der Pol Duffing oscillators, nonlinear Bloch equations modeling nuclear magnetic resonance, electrical circuits modeling nonlinear “jerk” equations, nonlinear equations of acoustic gravitational waves, and other chaotic systems.

We propose, for this study, the Rucklidge-type dynamical system described in [2] with four dimensionless physical parameters (a, b, c, $d>0$):

(1)$\left\{\begin{array}{c}\dot{x}=-ax+by-cyz\hfill \\ \dot{y}=dx\hfill \\ \dot{z}=y^2-z\hfill \end{array}\right.,$

and with the initial conditions for the state variables $(x,y,z)\in {\mathbb{R}}^3$ given by

(2)$\begin{array}{c}x(0)=x_{i0},y(0)=y_{i0},z(0)=z_{i0}.\hfill \end{array}$

The Rucklidge chaotic system (1) and its equilibria, as well as their basic dynamical properties, were evaluated: Lyapunov exponents and bifurcation diagrams (as discussed in [2]).

For $c=d=1$, the existence of Darboux and analytic first integrals of System (1) were studied in [3]. Specifically, it was investigated that System (1) had no analytic first integral, but the linear part of System (1) had two independent first integrals in each of the cases: $a=0$, $b\ne 0$; $a\ne 0$, $b=0$; $ab\ne 0$, $b=-{\displaystyle \frac{a^2}{4}}$, and $ab\ne 0$, $b\ne -{\displaystyle \frac{a^2}{4}}$, respectively.

It is easy to see that the eigenvalues of the linear part of System (1) are ${\lambda }_1=-1<0$, ${\lambda }_2={\displaystyle \frac{-a+\sqrt{a^2+4bd}}{2}}>0$, ${\lambda }_3={\displaystyle \frac{-a-\sqrt{a^2+4bd}}{2}}<0$.

Based on these considerations, we built semi-analytical solutions of System (1) with their behaviors, asymptotic or damped oscillations, using the OAFM.

This work has the following structure: Section 1 introduces the general properties of the Rucklidge-type dynamical system. Section 2 describes the Optimal Auxiliary Functions Method (OAFM) and its application in deriving semi-analytical solutions. Section 3 is dedicated to the numerical results and validation of the applied method. Lastly, Section 4 highlights the main conclusions.

2. The Optimal Auxiliary Functions Method (OAFM)

2.1. Stepwise OAFM

The steps of the OAFM were presented in [11,12], where the differential equation in general form was considered:

(3)$\begin{array}{c}{\displaystyle \mathcal{L}\left[\overline{v}(t)\right]+g(t)+\mathcal{N}\left[\overline{v}(t)\right]=0,}\hfill \\ \mathcal{B}(\overline{v}(0),\dot{\overline{v}}(0))=0,\hfill \end{array}$

where $\mathcal{L}$ is a linear operator; g is a known function; and $\mathcal{N}$ is a nonlinear operator, where t is the independent variable, and the approximate solution $\overline{v}(t)$ is written with just two components in the following form:

(4)${\displaystyle \overline{v}(t)=v_0(t)+v_1(t,C_i),i=1,2,\dots ,s.}$

Remark 1. 

The linear operator $\mathcal{L}$ was chosen and the initial approximation $v_0(t)$ became an elementary function or combination of the elementary functions. These functions could be as follows: the exponential function $e^{-Kt}$ or the rational function ${\displaystyle \frac{1}{1+Kt}}$ in the case of the boundary value problems from fluid mechanics; the trigonometric functions $cos({\omega }_{0}t)$ and $sin({\omega }_{0}t)$, which describe the nonlinear vibrations with periodic behaviors; $e^{-Kt}cos({\omega }_{0}t)$, $e^{-Kt}sin({\omega }_{0}t)$, modeling the nonlinear vibrations with harmonic/anharmonic oscillations—i.e., the damping effect.

The next step was to obtain the initial approximation $v_0(t)$ and the first approximation $v_1(t,C_i)$.

Firstly, the initial approximation $v_0(t)$ was chosen as a solution of the following equation:

(5)$\begin{array}{c}{\displaystyle \mathcal{L}\left[v_0(t)\right]+g(t)=0,}\hfill \\ \mathcal{B}(v_0(0),{\dot{v}}_0(0))=0,\hfill \end{array}$

and Equation (3) became

(6)$\begin{array}{c}{\displaystyle \mathcal{L}\left[v_0(t)\right]+\mathcal{L}\left[v_1(t,C_i)\right]+g(t)+\mathcal{N}\left[v_0(t)+v_1(t,C_i)\right]=0.}\hfill \end{array}$

The nonlinear operator can be expanded in the form

(7)$\begin{array}{c}{\displaystyle \mathcal{N}\left[v_0(t)+v_1(t,C_i)\right]=\mathcal{N}\left[v_0(t)\right]+\sum _{k=1}^{\infty }\frac{v_1^k(t,C_i)}{k!}{\mathcal{N}}^{(k)}\left[v_0(t)\right].}\hfill \end{array}$

With two arbitrary auxiliary functions ${\mathcal{A}}_1$ and ${\mathcal{A}}_2$—which depend on the initial approximation $v_0(t)$ and unknown parameters $C_i$ and $C_j$, $i=1,2,\dots ,p$, $j=p+1,p+2,\dots ,s$—then, using Equations (6) and (7), the first approximation $v_1(t)$ is the solution of the problem:

(8)$\begin{array}{c}{\displaystyle \mathcal{L}\left[v_1(t,C_i)\right]+{\mathcal{A}}_1\left[v_0(t),C_i\right]\mathcal{N}\left[v_0(t)\right]+{\mathcal{A}}_2\left[v_0(t),C_j\right]=0,}\hfill \end{array}$

(9)${\displaystyle v_1(0,C_i)=0,{\dot{v}}_1(0,C_i)=0.}$

Finally, the semi-analytical solution was given using Equation (4) in the following form:

(10)${\displaystyle \overline{v}(t)=v_0(t)+v_1(t,C_i),i=1,2,\dots ,s,}$

where $v_0(t)$ and $v_1(t,C_i)$ are solutions of Equations (5) and (8), respectively.

In OAFM, the linear operator $\mathcal{L}$ is arbitrarily chosen, not the physical parameters. There are situations when the selection of the physical parameters leads to chaotic behavior.This happened in the case of choosing higher values for the damping factor or for exceeding the optimal resonance conditions, as well as in the case of arbitrary chosen of initial conditions.

Remark 2. 

The nonlinear operator ${\mathcal{N}}^{(k)}\left[v_0(t)\right]$ from Equation (7) has the form

(11)$\begin{array}{c}{\displaystyle {\mathcal{N}}^{(k)}\left[v_0(t)\right]=\sum _{i=1}^{n_s}{h}_i(t)g_i(t),k=1,2,\dots ,}\hfill \end{array}$

where $n_s$ is a positive integer, and $h_i(t)$ and $g_i(t)$ are known functions that depend on $u_0(t)$.

This justifies the form of Equation (8) as

(12)$\begin{array}{c}{\displaystyle \mathcal{L}\left[v_1(t,C_i)\right]+{\mathcal{A}}_1\left[v_0(t),C_i\right]h_1(t)+{\mathcal{A}}_2\left[v_0(t),C_j\right]h_2(t)+{\mathcal{A}}_3\left[v_0(t),C_k\right]h_3(t)+\dots =0,}\hfill \end{array}$

where ${\mathcal{A}}_1\left[v_0(t),C_i\right]$, ${\mathcal{A}}_2\left[v_0(t),C_j\right]$, ${\mathcal{A}}_3\left[v_0(t),C_k\right]$, … arbitrary auxiliary functions. These functions depend on the unknown parameters $C_i$, $C_j$, $C_k$, … being optimally computed via various methods: the collocation method, the least squares method, the Galerkin method, the weighted residual method, etc.

Using the linearly independent functions $h_1,h_2,\dots ,h_{n_s}$, we introduced some types of approximate solutions of Equation (3).

Definition 1. 

A sequence of functions ${\left\{s_m(t)\right\}}_{m\ge 1}$ of the form

(13)${\displaystyle s_m(t)=\sum _{i=1}^m{\alpha }_m^i·h_i(t),m\ge 1,{\alpha }_m^i\in \mathbb{R}}$

is called an OAFM sequence of Equation (3).

Functions of the OAFM sequences are called OAFM functions of Equation (3).

The OAFM sequences ${\left\{s_m(t)\right\}}_{m\ge 1}$ with the property

${\displaystyle \underset{m\to \infty }{\lim }\mathcal{R}(t,s_m(t))=0}$

are called convergent to the solution of Equation (3), where $\mathcal{R}(t,u(t))=\mathcal{L}[u(t)]+\mathcal{N}[u(t)]+g(t)$.

Definition 2. 

The OAFM functions $\tilde{F}$ satisfying the conditions

(14)${\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,\tilde{F}(t))dt\le \epsilon ,\mathcal{B}\left(\tilde{F}(t,C_i),\frac{d\tilde{F}(t,C_i)}{dt}\right)=0}$

are called weak ε-approximate OAFM solutions of Equation (3) on the real interval $(0,\infty )$.

Remark 3. 

An ε-approximate OAFM solution of Equation (3) is also a weak ε-approximate OAFM solution. It follows that the set of weak ε-approximate OAFM solutions of Equation (3) also contains the approximate OAFM solutions of Equation (3).

The existence of weak $\epsilon$-approximate OAFM solutions is built by the theorem presented above.

Theorem 1. 

Equation (3) admits a sequence of weak ε-approximate OAFM solutions.

Proof. 

This is similar to the theorem from [13].

Firstly, the OAFM sequences ${\left\{s_m\right\}}_{m\ge 1}$ are built by considering the approximate OAFM solutions of the following type:

(15)${\displaystyle \overline{F}(t)=\sum _{i=1}^n{C}_m^i·h_i(t),\mathrm{where}m\ge 1\mathrm{is}\mathrm{fixed}\mathrm{arbitrary}.}$

The unknown parameters $C_m^i$, $i\in \{1,2,\dots ,n\}$ will then be determined.

If the approximate solutions $\overline{F}$ is introduced in Equation (3), the following expression is yielded:

$\mathcal{R}(t,C_m^i):=\mathcal{R}(t,\overline{F})=\mathcal{L}\left[\overline{F}(t)\right]+\mathcal{N}\left[\overline{F}(t)\right]+g(t).$

By attaching to Equation (3) the following real functional

(16)${\displaystyle {\mathcal{J}}_1(C_m^i)=\underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,C_m^i)dt}$

and imposing the initial conditions, we can determine $l\in \mathbb{N}$, $l\le m$, such that $C_m^1$, $C_m^2$, …, $C_m^l$ are computed as $C_m^{l+1}$, $C_m^{l+2}$, …, $C_m^n$.

The values of ${\tilde{C}}_m^{l+1}$, ${\tilde{C}}_m^{l+2}$, …, ${\tilde{C}}_m^n$ are computed by replacing $C_m^1$, $C_m^2$, …, $C_m^l$ in Equation (16), which give the minimum of Functional (16).

By means of the initial conditions, the values ${\tilde{C}}_m^1$, ${\tilde{C}}_m^2$, …, ${\tilde{C}}_m^l$ as functions of ${\tilde{C}}_m^{l+1}$, ${\tilde{C}}_m^{l+2}$, …, ${\tilde{C}}_m^n$ are determined.

Using the constants ${\tilde{C}}_m^1$, ${\tilde{C}}_m^2$, …, ${\tilde{C}}_m^n$ thus determined, the following OAFM functions

(17)${\displaystyle s_m(t)=\sum _{i=1}^n{\tilde{C}}_m^i·h_i(t)}$

are constructed.

We propose to demonstrate that the OAFM functions $s_m(t)$ are weak $\epsilon$-approximate OAFM solutions of Equation (3).

By computing OAFM functions $s_m(t)$ and taking into account that the $\overline{F}$ given by (15) are OAFM functions for Equation (3), then the following result is obtained:

${\displaystyle 0\le \underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,s_m(t))dt\le \underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,\overline{F}(t))dt,\forall m\ge 1.}$

Thus,

${\displaystyle 0\le \underset{m\to \infty }{\lim }\underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,s_m(t))dt\le \underset{m\to \infty }{\lim }\underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,\overline{F}(t))dt.}$

Since $\overline{F}(t)$ is convergent to the solution of Equation (3), we obtain the following:

${\displaystyle \underset{m\to \infty }{\lim }\underset{0}{\overset{\infty }{\int }}{\mathcal{R}}^2(t,s_m(t))dt=0.}$

It follows that, for all $\epsilon >0$, there exists $m_0\ge 1$ such that, for all $m\ge 1$, $m>m_0$, the sequence $s_m(t)$ is a weak $\epsilon$-approximate OAFM solution of Equation (3). □

Remark 4. 

The proof of the above theorem gives us a way to determine a weak ε-approximate OAFM solution of Equation (3), $\overline{F}$. Moreover, taking into account Remark 3, if $|\mathcal{R}(t,\overline{F})|<\epsilon$, then $\overline{F}$ is also an ε-approximate OAFM solution of the considered equation.

In the following, a special case of System (1) when $b=-{\displaystyle \frac{a-\frac{1}{2}}{2d}}$ with $a<{\displaystyle \frac{1}{2}}$, is analytically approached using the OAFM. If $a>0.2$, then the solution of the system describes a damped oscillatory motion. If a approaches $1/2$, then the solution of the system has asymptotic behavior.

In the general case of the physical parameters a, b, c, $d>0$, the initial System (1) is reduced to a two-dimensional differential system whose solution is given by exact parametric form.

2.2. Semi-Analytical Solutions via the OAFM

Case 1: $b=-(a-{\displaystyle \frac{1}{2}})/(2d)$, with $a>0.2$.

Integrating this with System (1) will obtain the following:

(18)$\left\{\begin{array}{c}{\displaystyle y(t)=y(0)+d\underset{0}{\overset{t}{\int }}x(s)ds}\hfill \\ z(t)=e^{-t}\left({\displaystyle z(0)+\underset{0}{\overset{t}{\int }}{y}^2(s)e^{s}ds}\right)\hfill \end{array}\right..$

Based on this, it is more convenient to determine a semi-analytical solution $\overline{x}(t)$ for the unknown function $x(t)$.

Figure 1 highlights the choice of a linear operator to build the semi-analytical solutions with the OAFM. From this figure, it can be seen that the solution of System (1) describes a damped oscillatory motion for the values of the parameter $a>0.2$. However, the amplitude of the oscillation increased as the parameter $a>0$ decreased.

In this case, the initial System (1) could be written as follows:

(19)$\left\{\begin{array}{c}\dot{x}-{\phi }_1(t)+\left[{\phi }_1(t)+ax-by+cyz\right]=0\hfill \\ \dot{y}-{\phi }_2(t)+\left[{\phi }_2(t)-dx\right]=0\hfill \\ \dot{z}-{\phi }_3(t)+\left[{\phi }_3(t)-y^2+z\right]=0\hfill \end{array}\right.,$

where ${\phi }_i(t)=-K_i{w}_i+\left({\tilde{C}}_{i}cos({\omega }_{0}t)+{\tilde{D}}_{i}sin({\omega }_{0}t)\right)e^{-K_{i}t}$, $i=1,2,3$, with $w_1=x$, $w_2=y$, $w_3=z$, are the given functions depending on some unknown parameters ${\tilde{C}}_i$, ${\tilde{D}}_i$, $K_i$, ${\omega }_0$.

In this manner, the linear and nonlinear operators $\mathcal{L}$, $\mathcal{N}$ that—from Equation (3)—correspond to the unknown function $x(t)$, are $\mathcal{L}(x(t))=\dot{x}-{\phi }_1(t)$ and $\mathcal{N}(x(t))=-K_{1}x+\left({\tilde{C}}_{i}cos({\omega }_{0}t)+{\tilde{D}}_{i}sin({\omega }_{0}t)\right)e^{-K_{i}t}+ax-by+cyz$, respectively.

Using the OAFM, a semi-analytical solution $\overline{x}(t)$ (taking into account Equation (10)) is written as follows:

(20)$\begin{array}{c}\overline{x}(t)=x_0(t)+x_1(t).\hfill \end{array}$

Equation (5) then becomes the following:

(21)$\begin{array}{c}\mathcal{L}(x_0(t))\equiv {\dot{x}}_0+K_1{x}_0-\left({\tilde{C}}_{1}cos({\omega }_{0}t)+{\tilde{D}}_{1}sin({\omega }_{0}t)\right)e^{-K_{1}t}=0,\hfill \end{array}$

with the solution

(22)$\begin{array}{c}{x}_0(t)=M_0{e}^{-K_{1}t}+M_1{e}^{-K_{1}t}cos({\omega }_{0}t)+M_2{e}^{-K_{1}t}sin({\omega }_{0}t),\hfill \end{array}$

for $M_0={\displaystyle \frac{{\tilde{D}}_1+x(0){\omega }_0}{{\omega }_0}}$, $M_1=-{\displaystyle \frac{{\tilde{D}}_1}{{\omega }_0}}$, $M_2={\displaystyle \frac{{\tilde{C}}_1}{{\omega }_0}}$.

Analogously, the initial approximations $y_0(t)$ and $z_0(t)$ can be obtained from Equation (5) using the linear operators $\mathcal{L}(y_0(t))={\dot{y}}_0+K_2{y}_0-\left({\tilde{C}}_{2}cos({\omega }_{0}t)+{\tilde{D}}_{2}sin({\omega }_{0}t)\right)e^{-K_{2}t}$ and $\mathcal{L}(z_0(t))={\dot{z}}_0+K_3{z}_0-\left({\tilde{C}}_{3}cos({\omega }_{0}t)+{\tilde{D}}_{3}sin({\omega }_{0}t)\right)e^{-K_{3}t}$, which are linear combinations of the following functions set:$\{e^{-K_{i}t},e^{-K_{i}t}cos({\omega }_{0}t),e^{-K_{i}t}sin({\omega }_{0}t)\}$, $i=2,3$.

For obtaining the first approximation $x_1(t)$, it can be observed that the nonlinear operator $\mathcal{N}(x_0(t))$ becomes the following:

(23)$\begin{array}{c}\mathcal{N}(x_0(t))=\left({\tilde{C}}_{i}cos({\omega }_{0}t)+{\tilde{D}}_{i}sin({\omega }_{0}t)\right)e^{-K_{i}t}+a{x}_0-b{y}_0+c{y}_0{z}_0.\hfill \end{array}$

This expression is a linear combination of functions set:

(24)$\begin{array}{c}\left\{e^{-K_{1}t},e^{-K_{2}t},e^{-(K_2+K_3)t},e^{-K_{1}t}cos({\omega }_{0}t),e^{-K_{2}t}cos({\omega }_{0}t),e^{-K_{1}t}sin({\omega }_{0}t),\right.\hfill \\ e^{-K_{2}t}sin({\omega }_{0}t),e^{-(K_2+K_3)t}cos({\omega }_{0}t),e^{-(K_2+K_3)t}sin({\omega }_{0}t),e^{-(K_2+K_3)t}cos(2{\omega }_{0}t),\hfill \\ \left.e^{-(K_2+K_3)t}sin(2{\omega }_{0}t)\right\}.\hfill \end{array}$

Combining Equations (11) and (24), the linear independent functions are identified:

$\begin{array}{cccc}{h}_1(t)=e^{-K_{1}t},\hfill & g_1(1)=1,\hfill & h_2(t)=e^{-K_{2}t},\hfill & g_2(1)=1,\hfill \\ h_3(t)=e^{-(K_2+K_3)t},\hfill & g_3(t)=1,\hfill & h_4(t)=e^{-K_{1}t},\hfill & g_4(t)=cos({\omega }_{0}t),\hfill \\ h_5(t)=e^{-K_{2}t},\hfill & g_5(t)=cos({\omega }_{0}t),\hfill & h_6(t)=e^{-K_{1}t},\hfill & g_6(t)=sin({\omega }_{0}t),\hfill \\ h_7(t)=e^{-K_{2}t},\hfill & g_7(t)=sin({\omega }_{0}t),\hfill & h_8(t)=e^{-(K_2+K_3)t},\hfill & g_8(t)=cos({\omega }_{0}t),\hfill \\ h_9(t)=e^{-(K_2+K_3)t},\hfill & g_9(t)=sin({\omega }_{0}t),\hfill & h_{10}(t)=e^{-(K_2+K_3)t},\hfill & g_{10}(t)=cos(2{\omega }_{0}t),\hfill \\ h_{11}(t)=e^{-(K_2+K_3)t},\hfill & g_{11}(t)=sin(2{\omega }_{0}t).\hfill \end{array}$

At this moment, the first approximation $x_1(t)$ of the unknown function $x(t)$ can be computed using Equation (12), which becomes the following:

(25)$\begin{array}{c}{\displaystyle \mathcal{L}(x_1(t))+\sum _{i=1}^{10}{\mathcal{A}}_i(t)h_i(t)=0.}\hfill \end{array}$

There are more possibilities for choosing the auxiliary functions ${\mathcal{A}}_i(t)$, $i=\overline{1,10}$ as follows:

(26)$\begin{array}{c}{\displaystyle {\mathcal{A}}_1(t)=\sum _{j=1}^{N_{max}}\left({(-1)}^{j}j{\tilde{C}}_{j}cos(j{\omega }_{0}t)+{(-1)}^{j}j{\tilde{D}}_{j}sin(j{\omega }_{0}t)\right),}\hfill \\ {\mathcal{A}}_k(t)=0,k=\overline{1,10},\hfill \end{array}$

or

$\begin{array}{c}{\displaystyle {\mathcal{A}}_4(t)=\sum _{j=1}^{N_{max}}\left({(-1)}^{j}j{\tilde{C}}_{j}cos(j{\omega }_{0}t)+{(-1)}^{j}j{\tilde{D}}_{j}sin(j{\omega }_{0}t)\right),}\hfill \\ {\mathcal{A}}_k(t)=0,k=\overline{1,10},k\ne 4,\hfill \end{array}$

or

$\begin{array}{c}{\displaystyle {\mathcal{A}}_6(t)=\sum _{j=1}^{N_{max}}\left({(-1)}^{j}j{\tilde{C}}_{j}cos(j{\omega }_{0}t)+{(-1)}^{j}j{\tilde{D}}_{j}sin(j{\omega }_{0}t)\right),}\hfill \\ {\mathcal{A}}_k(t)=0,k=\overline{1,10},k\ne 6,\hfill \end{array}$

and so on, where ${\tilde{C}}_j$, ${\tilde{D}}_j$, $j=\overline{1,N_{max}}$ are unknown parameters, and $N_{max}$ is an arbitrary fixed integer number.

The first approximation $x_1(t)$ could be obtained using the auxiliary functions defined by Equation (26) and by integration of Equation (25) having the following form:

(27)$\begin{array}{c}{\displaystyle x_1(t)=\sum _{j=1}^{N_{max}}\left({(-1)}^{j}j{C}_{j}cos(j{\omega }_{0}t)+{(-1)}^{j}j{D}_{j}sin(j{\omega }_{0}t)\right)e^{-K_{1}t},}\hfill \end{array}$

where $N_{max}$ is an arbitrary fixed integer, and ${\omega }_0$, $K_1$, $C_j$, $D_j$, $j=\overline{1,N_{max}}$ are unknown convergence-control parameters depending on ${\tilde{C}}_j$, ${\tilde{D}}_j$, $j=\overline{1,N_{max}}$ that will be optimally identified at the end.

Therefore, the first-order approximate solution $\overline{x}(t)$ defined by Equation (20) is well determined by Equations (22) and (27) via the OAFM:

(28)$\begin{array}{c}\overline{x}(t)=M_0{e}^{-K_{1}t}+M_1{e}^{-K_{1}t}cos({\omega }_{0}t)+M_2{e}^{-K_{1}t}sin({\omega }_{0}t)+\hfill \\ {\displaystyle +\sum _{j=1}^{N_{max}}\left({(-1)}^{j}j{C}_{j}cos(j{\omega }_{0}t)+{(-1)}^{j}j{D}_{j}sin(j{\omega }_{0}t)\right)e^{-K_{1}t}.}\hfill \end{array}$

Combining Equations (18) and (20), the semi-analytical solutions for unknown functions $y(t)$ and $z(t)$ are obtained as follows:

(29)$\left\{\begin{array}{c}{\displaystyle \overline{y}(t)=y(0)+d\underset{0}{\overset{t}{\int }}\overline{x}(s)ds}\hfill \\ \overline{z}(t)=e^{-t}\left({\displaystyle z(0)+\underset{0}{\overset{t}{\int }}{\overline{y}}^2(s)e^{s}ds}\right)\hfill \end{array}\right.,$

where $\overline{x}(t)$ is defined by Equation (28).

Case 2: $b=-(a-{\displaystyle \frac{1}{2}})/(2d)$, with a closer to ${\displaystyle \frac{1}{2}}$.

As in the previous case, as shown in Figure 2, it can be seen that the solution of System (1) has asymptotic behavior.

In this case, the initial System (1) could be written as follows:

(30)$\left\{\begin{array}{c}\dot{x}+K_{1}x+\left[-K_{1}x+ax-by+cyz\right]=0\hfill \\ \dot{y}+K_{2}y+\left[-K_{2}y-dx\right]=0\hfill \\ \dot{z}+K_{3}z+\left[-K_{3}z-y^2+z\right]=0\hfill \end{array}\right.,$

where ${\tilde{C}}_i$, ${\tilde{D}}_i$, $K_i$, ${\omega }_0$ are unknown parameters at this moment.

In the same manner, the linear $\mathcal{L}$ and nonlinear operators $\mathcal{N}$ from Equation (3) that correspond to unknown function $x(t)$ are $\mathcal{L}(x(t))=\dot{x}+K_{1}x$ and $\mathcal{N}(x(t))=-K_{1}x+ax-by+cyz$, respectively.

By means of the OAFM, a semi-analytical solution $\overline{x}(t)$ has the same form as Equation (20).

Equation (5) becomes:

(31)$\begin{array}{c}\mathcal{L}(x_0(t))\equiv {\dot{x}}_0+K_1{x}_0=0,\hfill \end{array}$

with the solution

(32)$\begin{array}{c}{x}_0(t)=x(0)e^{-K_{1}t}.\hfill \end{array}$

Analogously, the initial approximations $y_0(t)$ and $z_0(t)$ can be obtained from Equation (5) using the linear operators $\mathcal{L}(y_0(t))={\dot{y}}_0+K_2{y}_0$ and $\mathcal{L}(z_0(t))={\dot{z}}_0+K_3{z}_0$, having the form $y_0(t)=y(0)e^{-K_{2}t}$, $z_0(t)=z(0)e^{-K_{3}t}$.

To build the first approximation $x_1(t)$, it can be observed that the nonlinear operator $\mathcal{N}(x_0(t))$ becomes the following:

(33)$\begin{array}{c}\mathcal{N}(x_0(t))=\left({\tilde{C}}_{i}cos({\omega }_{0}t)+{\tilde{D}}_{i}sin({\omega }_{0}t)\right)e^{-K_{i}t}+a{x}_0-b{y}_0+c{y}_0{z}_0=\hfill \\ (a-K_1)x(0)e^{-K_{1}t}-by(0)e^{-K_{2}t}+cy(0)z(0)e^{-(K_2+K_3)t}.\hfill \end{array}$

This expression is a linear combination of the following functions set:

(34)$\begin{array}{c}\{e^{-K_{1}t},e^{-K_{2}t},e^{-(K_2+K_3)t}\}.\hfill \end{array}$

When combining Equations (11) and (24), the linear independent functions are identified as follows:

$\begin{array}{cccc}{h}_1(t)=e^{-K_{1}t},\hfill & g_1(1)=1,\hfill & h_2(t)=e^{-K_{2}t},\hfill & g_2(1)=1,\hfill \\ h_3(t)=e^{-(K_2+K_3)t},\hfill & g_3(t)=1.\hfill \end{array}$

At this moment, the first approximation $x_1(t)$ of the unknown function $x(t)$ can be computed using Equation (12), which becomes the following:

(35)$\begin{array}{c}{\displaystyle \mathcal{L}(x_1(t))+\sum _{i=1}^3{\mathcal{A}}_i(t)h_i(t)=0.}\hfill \end{array}$

There are more possibilities to choose the auxiliary functions ${\mathcal{A}}_i(t)$, $i=\overline{1,3}$ as follows:

(36)$\begin{array}{c}{\displaystyle {\mathcal{A}}_1(t)=({\tilde{D}}_1+{\tilde{D}}_{2}t+{\tilde{D}}_3{t}^2+{\tilde{C}}_1{t}^3)e^{-K_{1}t},}\hfill \\ {\displaystyle {\mathcal{A}}_2(t)={\tilde{D}}_4+{\tilde{D}}_{5}t+{\tilde{D}}_6{t}^2+{\tilde{C}}_2{t}^3}\hfill \\ {\displaystyle {\mathcal{A}}_3(t)={\tilde{D}}_7+{\tilde{D}}_{8}t+{\tilde{D}}_9{t}^2+{\tilde{C}}_3{t}^3,}\hfill \end{array}$

or

$\begin{array}{c}{\displaystyle {\mathcal{A}}_1(t)=({\tilde{D}}_1+{\tilde{D}}_{2}t+{\tilde{D}}_3{t}^2)e^{-2K_{1}t},}\hfill \\ {\displaystyle {\mathcal{A}}_2(t)=({\tilde{D}}_4+{\tilde{D}}_{5}t+{\tilde{D}}_6{t}^2)e^{-K_{2}t}}\hfill \\ {\displaystyle {\mathcal{A}}_3(t)={\tilde{D}}_7+{\tilde{D}}_{8}t+{\tilde{D}}_9{t}^2+{\tilde{C}}_3{t}^3,}\hfill \end{array}$

or

$\begin{array}{c}{\displaystyle {\mathcal{A}}_1(t)=({\tilde{D}}_1+{\tilde{D}}_{2}t)e^{-K_{1}t},}\hfill \\ {\displaystyle {\mathcal{A}}_2(t)=({\tilde{D}}_4+{\tilde{D}}_{5}t)e^{-K_{2}t}}\hfill \\ {\displaystyle {\mathcal{A}}_3(t)={\tilde{D}}_7+{\tilde{D}}_{8}t+{\tilde{D}}_9{t}^2+{\tilde{C}}_3{t}^3,}\hfill \end{array}$

and so on, where ${\tilde{C}}_i$, ${\tilde{D}}_j$ are unknown parameters.

The first approximation $x_1(t)$ could be obtained using the auxiliary functions defined by Equation (36) and by integrating Equation (35). This solution has the following form:

(37)$\begin{array}{c}{\displaystyle x_1(t)=(D_1+D_{2}t+D_3{t}^2+C_1{t}^3)e^{-2K_{1}t}+(D_4+D_{5}t+D_6{t}^2+C_2{t}^3)e^{-K_{2}t}+}\hfill \\ +(D_7+D_{8}t+D_9{t}^2+C_3{t}^3)e^{-(K_2+K_3)t},\hfill \end{array}$

where $K_1$, $K_2$, $K_2$, $C_i$, $i=1,2,3$, $D_j$, $j=\overline{1,9}$ are unknown convergence-control parameters that depending on physical parameters $a,b,c,d$, initial conditions $x(0),y(0),z(0)$, and ${\tilde{C}}_i$, ${\tilde{D}}_j$. These parameters will be optimally identified at the end.

Therefore, the first-order approximate solution $\overline{x}(t)$ defined by Equation (20) is well determined by Equations (32) and (37) via the OAFM as follows:

(38)$\begin{array}{c}{\displaystyle \overline{x}(t)=x(0)e^{-K_{1}t}+(D_1+D_{2}t+D_3{t}^2+C_1{t}^3)e^{-2K_{1}t}+(D_4+D_{5}t+D_6{t}^2+C_2{t}^3)e^{-K_{2}t}+}\hfill \\ +(D_7+D_{8}t+D_9{t}^2+C_3{t}^3)e^{-(K_s)t}.\hfill \end{array}$

The semi-analytical solutions were built using Equations (18) and (38) for unknown functions $y(t)$ and $z(t)$.

Case 3:a, b, c, d are arbitrarily chosen (damped oscillatory motion).

In this case, a semi-analytical solution can be determined in parametric form. By changing the variables $x(t)=X(y(t))$ and $b-cz(t)=Z(y(t))$, the initial system becomes the following:

(39)$\left\{\begin{array}{c}dX{\displaystyle \frac{dX}{dy}}=-aX+yZ\hfill \\ dX{\displaystyle \frac{dZ}{dy}}=y^2+{\displaystyle \frac{1}{c}}Z-{\displaystyle \frac{b}{c}}\hfill \end{array}\right.$

subject to the initial conditions

(40)$\begin{array}{c}X(y(0))=x_{i0},Z(y(0))=z_{i0}.\hfill \end{array}$

It was assumed that the existence of a constant $K_0>0$ is such that $|y_{i0}-{\alpha }_0| <K_0$, where ${\alpha }_0=\sqrt{{\displaystyle \frac{b}{c}}}$. In this case the function $f(y)=y$ is an odd function and the function $g(y)=y^2-{\displaystyle \frac{b}{c}}$ is an even function. The expansion of this is as follows:

(41)$\begin{array}{c}{\displaystyle f(y)=y=(y-{\alpha }_0)+{\alpha }_0}\hfill \\ {\displaystyle g(y)=y^2-\frac{b}{c}={(y-{\alpha }_0)}^2+2{\alpha }_0(y-{\alpha }_0).}\hfill \end{array}$

This idea suggests building exact parametric solutions by power series expansion as follows:

(42)$\begin{array}{c}{\displaystyle X(y)=\sum _{n\ge 0}{(-1)}^n\frac{C_{2n+1}}{2n+1}{\left(\frac{y-{\alpha }_0}{K_0}\right)}^{2n+1},}\hfill \\ {\displaystyle Z(y)=D_0+\sum _{n\ge 1}{(-1)}^n\frac{D_{2n}}{2n}{\left(\frac{y-{\alpha }_0}{K_0}\right)}^{2n}.}\hfill \end{array}$

The power series from Equation (42) are absolutely convergent (via Abel’s Theorem from Mathematical Analysis).

The semi-analytical solutions $\overline{X}(y)$ and $\overline{Z}(y)$ can be obtained from Equation (42) as follows:

(43)$\begin{array}{c}{\displaystyle \overline{X}(y)=\sum _{n=0}^{N_{max}}{(-1)}^n\frac{C_{2n+1}}{2n+1}{\left(\frac{y-{\alpha }_0}{K_0}\right)}^{2n+1},}\hfill \\ {\displaystyle \overline{Z}(y)=D_0+\sum _{n=1}^{N_{max}}{(-1)}^n\frac{D_{2n}}{2n}{\left(\frac{y-{\alpha }_0}{K_0}\right)}^{2n},}\hfill \end{array}$

where $D_0$, $C_{2n+1}$, $D_{2n}$, $n=\overline{1,N_{max}}$ are unknown convergence-control parameters and will be optimally computed at the end.

3. Numerical Results and Validation

The OAFM solutions are presented in this section for two mentioned cases. The comparative analysis between the OAFM solutions ${\overline{u}}_{OAFM}$ and the corresponding numerical ones, for given initial conditions and physical constants a, b, c, d, are presented in detail in Table 1 and Table 2, respectively. The fourth-order Runge–Kutta method was used for testing the numerical method. The absolute values denoted by ${ϵ}_u=|u_{numerical}-{\overline{u}}_{OAFM}|$ dispute the accuracy of the obtained results in these tables.

Figure 3 and Figure 4 qualitatively reflect the agreement of the OAFM solutions ${\overline{u}}_{OAFM}$ and the states $(x(t),y(t),z(t))$ of the studied system with the corresponding numerical solutions.

The convergence-control parameters from Equations (28) and (38), respectively, are given in details in the Appendix A.

The very good agreement between the values for the function ${\overline{x}}_{OAFM}$, the numerical function $x_{numerical}$, for $x_{i0}=0.25$, $y_{i0}=0.5$, $z_{i0}=1.5$; the physical constants $a=0.27$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$, index $N_{max}=15$, results from the calculus of the absolute difference ${ϵ}_x=|x_{numerical}-{\overline{x}}_{OAFM}|$ and its order of magnitude (${10}^{-4}$–${10}^{-6}$).

The very good agreement between the values for the function ${\overline{x}}_{OAFM}$, as well as the numerical function $x_{numerical}$ for $x_{i0}=0.25$, $y_{i0}=0.5$, $z_{i0}=1.5$; the physical constants $a=0.49$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$, index $N_{max}=15$ results from the calculus of the absolute difference ${ϵ}_x=|x_{numerical}-{\overline{x}}_{OAFM}|$ and its order of magnitude (${10}^{-4}$–${10}^{-5}$).

For $a=0.56$, $b=0.84$, $c=0.275$, $d=2.165$, and the initial conditions $x_{i0}=0.25$, $y_{i0}=0.5$, and $z_{i0}=1.5$, we have $|y_{i0}-{\alpha }_0|=\left|y_{i0}-\sqrt{{\displaystyle \frac{b}{c}}}\right|=|-1.2477257950106058|<2=K_0$.

In this case, the semi-analytical solutions $\overline{X}(y)$ and $\overline{Z}(y)$ of System (39), taking into account Equation (43) for $N_{max}=10$ is as follows:

(44)$\begin{array}{c}{\displaystyle \overline{X}(y)=-4.2409358275\left(\frac{y-{\alpha }_0}{K_0}\right)+24.1508630465{\left(\frac{y-{\alpha }_0}{K_0}\right)}^3-}\hfill \\ 95.1251514294{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^5+254.2470586656{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^7-466.6522066628{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^9+\hfill \\ 594.9807405785{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{11}-525.8772313084{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{13}+315.8505259783{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{15}-\hfill \\ 122.9173308772{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{17}+27.9371537376{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{19}-2.8138962930{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{21},\hfill \\ \\ {\displaystyle \overline{Z}(y)=0.9217106664+0.8127141641{\left(\frac{y-{\alpha }_0}{K_0}\right)}^2-1.7309576021{\left(\frac{y-{\alpha }_0}{K_0}\right)}^4+}\hfill \\ 3.8071679067{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^6-6.0053357228{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^8+6.6628016497{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{10}-\hfill \\ 5.0246593808{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{12}+2.4616626473{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{14}-0.7123667470{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{16}+\hfill \\ 0.0965502691{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{18}-0.0017049572{\left({\displaystyle \frac{y-{\alpha }_0}{K_0}}\right)}^{20}.\hfill \end{array}$

The behavior of these solutions is graphically depicted in Figure 5. The variation of the functions ${\epsilon }_X=|{\overline{X}}_{OAFM}(y)-X_{numerical}(y)|$ and ${\epsilon }_Z=|{\overline{Z}}_{OAFM}(y)-Z_{numerical}(y)|$ is qualitatively presented in Figure 6 and Figure 7, respectively.

The OAFM via the Iterative Method

Recently, M. Farooq et al. [14] investigated the mechanical problem of the steady flow of couple stress fluid between two infinitely parallel inclined plates under the impact of MHD using the comparison between the techniques of HAM and OAFM. The performance of the OAFM was highlighted by the residual error of the OAFM and HAM solutions for the velocity.

The iterative method proposed by Daftardar-Gejji et al. [15] is used for solving nonlinear functional equations, such as the fractional-order differential equation, the nonlinear 3D-differential system, and the Voltera-type integral equation. This technique is only validated for a special class of nonlinear functional equations with known exact solutions.

By integration of System (1) over the interval $[0,t]$, we obtained the following results:

(45)$\begin{array}{c}{\displaystyle x(t)=x(0)+\underset{0}{\overset{t}{\int }}\left(-ax(s)+by(s)-cy(s)z(s)\right)ds}\hfill \\ {\displaystyle y(t)=y(0)+\underset{0}{\overset{t}{\int }}dx(s)ds}\hfill \\ {\displaystyle z(t)=z(0)+\underset{0}{\overset{t}{\int }}\left(y^2(s)-z(s)\right)ds}\hfill \end{array}.$

The iterative algorithm led to the following:

(46)$\begin{array}{c}{\displaystyle x_0(t)=x(0),x_1(t)=N_1(x_0,y_0,z_0)=\underset{0}{\overset{t}{\int }}\left(-ax(0)+by(0)-cy(0)z(0)\right)ds,}\hfill \\ {\displaystyle y_0(t)=y(0),y_1(t)=N_2(x_0,y_0,z_0)=\underset{0}{\overset{t}{\int }}dx(0)ds,}\hfill \\ {\displaystyle z_0(t)=z(0),z_1(t)=N_3(x_0,y_0,z_0)=\underset{0}{\overset{t}{\int }}\left(y^2(0)-z(0)\right)ds,}\hfill \\ \dots \hfill \\ x_m(t)=N_1\left({\displaystyle \sum _{i=0}^{m-1}{x}_i,\sum _{i=0}^{m-1}{y}_i,\sum _{i=0}^{m-1}{z}_i}\right)-N_1\left({\displaystyle \sum _{i=0}^{m-2}{x}_i,\sum _{i=0}^{m-2}{y}_i,\sum _{i=0}^{m-2}{z}_i}\right),\hfill \\ y_m(t)=N_2\left({\displaystyle \sum _{i=0}^{m-1}{x}_i,\sum _{i=0}^{m-1}{y}_i,\sum _{i=0}^{m-1}{z}_i}\right)-N_2\left({\displaystyle \sum _{i=0}^{m-2}{x}_i,\sum _{i=0}^{m-2}{y}_i,\sum _{i=0}^{m-2}{z}_i}\right),\hfill \\ z_m(t)=N_3\left({\displaystyle \sum _{i=0}^{m-1}{x}_i,\sum _{i=0}^{m-1}{y}_i,\sum _{i=0}^{m-1}{z}_i}\right)-N_3\left({\displaystyle \sum _{i=0}^{m-2}{x}_i,\sum _{i=0}^{m-2}{y}_i,\sum _{i=0}^{m-2}{z}_i}\right),\hfill \\ m\ge 2.\hfill \end{array}$

The iterative solution generated the solution of Equation (1) as

${\displaystyle x_{iter}(t)=\sum _{m=0}^{\infty }{x}_m(t),y_{iter}(t)=\sum _{m=0}^{\infty }{y}_m(t),z_{iter}(t)=\sum _{m=0}^{\infty }{z}_m(t).}$

For the initial conditions $x(0)=0.25$, $y(0)=0.5$, and $z(0)=1.5$; and the physical constants $a=0.27$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.650$ (presented in Table 3 and Table 4), the iterative solutions $x_{iter}(t)$ (after seven iterations using Algorithm (46)) have the following form:

(47)$\begin{array}{c}{\displaystyle x_{iter}(t)=\sum _{m=0}^7{x}_m(t)=0.25-0.5415384615t+0.2304514423t^2-}\hfill \\ 0.0026201169t^3-0.0476812213t^4+0.0324036061t^5-0.0123571245t^6+\hfill \\ 0.0022921155t^7+0.0004469592t^8-0.0005745121t^9+0.0002413778t^{10}-\hfill \\ 0.0000509306t^{11}-1.4110\times {10}^{-6}{t}^{12}+5.0074\times {10}^{-6}{t}^{13}-\hfill \\ 1.69095\times {10}^{-6}{t}^{14}+1.4761\times {10}^{-7}{t}^{15}+1.13293\times {10}^{-7}{t}^{16}-\hfill \\ 5.8624\times {10}^{-8}{t}^{17}+1.16251\times {10}^{-8}{t}^{18}+9.6537\times {10}^{-10}{t}^{19}-\hfill \\ 1.3883\times {10}^{-9}{t}^{20}+4.4314\times {10}^{-10}{t}^{21}-5.9803\times {10}^{-11}{t}^{22}-\hfill \\ 7.7661\times {10}^{-12}{t}^{23}+5.8424\times {10}^{-12}{t}^{24}-1.3419\times {10}^{-12}{t}^{25}+\hfill \\ 1.0590\times {10}^{-13}{t}^{26}+\mathcal{O}({10}^{-14}),\hfill \\ \\ {\displaystyle y_{iter}(t)=\sum _{m=0}^7{y}_m(t)=0.5+0.1625t-0.176t^2+0.0499311458t^3-}\hfill \\ 0.0004257690t^4-0.0061985587t^5+0.0035103906t^6-0.0011474472t^7+\hfill \\ 0.0001981113t^8+0.0000437522t^9-0.0000406537t^{10}+0.0000130644t^{11}-\hfill \\ 2.0380\times {10}^{-6}{t}^{12}-9.3442\times {10}^{-8}{t}^{13}+1.3163\times {10}^{-7}{t}^{14}-\hfill \\ 1.9260\times {10}^{-8}{t}^{15}-5.8684\times {10}^{-9}{t}^{16}+3.0713\times {10}^{-9}{t}^{17}-\hfill \\ 4.3939\times {10}^{-10}{t}^{18}-5.0669\times {10}^{-11}{t}^{19}+3.2836\times {10}^{-11}{t}^{20}-\hfill \\ 5.5704\times {10}^{-12}{t}^{21}+1.4442\times {10}^{-13}{t}^{22}+8.5337\times {10}^{-14}{t}^{23}+\mathcal{O}({10}^{-15}),\hfill \\ {\displaystyle z_{iter}(t)=\sum _{m=0}^7{z}_m(t)=1.5-1.2499999999t+0.7062499999t^2-}\hfill \\ 0.285281250t^3+0.0695030989t^4-0.004545049t^5-0.0032279413t^6+\hfill \\ 0.001052398t^7+0.0000477133t^8-0.000259671t^9+0.0001207900t^{10}-\hfill \\ 0.000022641t^{11}-3.5903\times {10}^{-6}{t}^{12}+3.8903\times {10}^{-6}{t}^{13}-\hfill \\ 1.3146\times {10}^{-6}{t}^{14}+2.0343\times {10}^{-7}{t}^{15}+2.5297\times {10}^{-8}{t}^{16}-\hfill \\ 2.5120\times {10}^{-8}{t}^{17}+7.3672\times {10}^{-9}{t}^{18}-8.9245\times {10}^{-10}{t}^{19}-\hfill \\ 1.6348\times {10}^{-10}{t}^{20}+9.8683\times {10}^{-11}{t}^{21}-1.8665\times {10}^{-11}{t}^{22}+\hfill \\ 2.8455\times {10}^{-13}{t}^{23}+6.3031\times {10}^{-13}{t}^{24}-1.2029\times {10}^{-13}{t}^{25}-\hfill \\ 6.3853\times {10}^{-16}{t}^{26}+3.1206\times {10}^{-15}{t}^{27}-2.7553\times {10}^{-16}{t}^{28}+\mathcal{O}({10}^{-17}).\hfill \end{array}$

Table 3 and Table 4, as well as Figure 8, show the accuracy of the OAFM by comparison with the iterative method (when using seven iterations).

4. Conclusions

This paper focused on an analytical approach to asymptotic behavior, but it also considered the damped oscillations of some solutions of the Rucklidge-type dynamical system using the OAFM. Damped harmonic oscillations appear naturally in many applications involving mechanical and electrical, as well as biological, systems. Understanding damped oscillations is essential for analyzing and designing these various systems. As is known, this dynamical system does not have Darboux polynomials, nor does it have primary polynomial integrals or invariant algebraic surfaces. Therefore, it does not have a Hamilton–Poisson structure. Thus, the exact solution can no longer be described as the intersection of two level surfaces $H(x,y,z)=constant$ (Hamiltonian function) and $C(x,y,z)=constant$ (independent Casimir function). The dynamics of the system is influenced by the values of the physical parameters. Two situations were analytically investigated, namely when the parameter a approaches 0 (damped oscillations) and when a approaches $1/2$ (asymptotic behaviors) for $b=-(a-1/2)/(2d)$. The case $b\ne -(a-1/2)/(2d)$ for which the solution of the system describes damped oscillations was also analytically approached. In this case, the solution was given in parametric form using trigonometric expansion series, and it was then presented as a semi-analytical solution.

The OAFM results were in good agreement with the corresponding numerical ones, as highlighted by the comparison between them in the presented figures and tables. The OAFM had some advantages to be applied. Some of them are mentioned below:

The OAFM solutions were written in effective form by arbitrary choice of the linear operator L and the auxiliary functions $A_j$;

Dependence of the auxiliary functions on the finite number of the unknown parameters were optimally computed to obtain the absolute values between the semi-analytical solutions and numerical ones smaller than one, which assured convergence control;

The OAFM solutions were obtained for arbitrary values of the physical parameters a, b, c, d, namely for $b=-(a-1/2)/(2d)$. In two subcases, a was closer to 0, and a approached $1/2$;

The OAFM is more efficient when compared with the iterative method by means of solutions, namely the iterative method is validated just when the exact solution is known.

Author Contributions

Conceptualization, N.P.; formal analysis, N.P.; investigation, R.-D.E. and R.B.; methodology, R.-D.E., R.B. and N.P.; software, R.-D.E. and R.B.; supervision, N.P.; validation, R.-D.E., R.B. and N.P.; visualization, R.-D.E., R.B. and N.P.; writing–original draft, R.-D.E., R.B. and N.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 The damped oscillatory behavior of the numerical solutions $x(t)$ of System (1) for specified values of the following parameters: a, $b=-(a-{\displaystyle \frac{1}{2}})/(2d)$, $c=0.75$, $d=0.65$.

Figure 2 The asymptotic behavior of the numerical solutions $(x(t),y(t),z(t))$ of System (1) for physical parameters $a=0.49$, $b=-(a-{\displaystyle \frac{1}{2}})/(2d)$, $c=0.75$, $d=0.65$ (asymptotic behaviors).

Figure 3 States of System (1) $({\overline{x}}_{OAFM},{\overline{y}}_{OAFM},{\overline{z}}_{OAFM})$ using Equations (28), (29) and (A1), for $x_{i0}=0.25$, $y_{i0}=0.5$, $z_{i0}=1.5$; for physical constants $a=0.27$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$ (dotted curves). This was achieved by comparison with the numerical integration results (solid curves of magenta for $x(t)$, red curve for $y(t)$, and blue for $z(t)$).

Figure 4 States of System (1) $({\overline{x}}_{OAFM},{\overline{y}}_{OAFM},{\overline{z}}_{OAFM})$ using Equations (38), (29) and (A2) for $x_{i0}=0.25$, $y_{i0}=0.5$, and $z_{i0}=1.5$; physical constants $a=0.49$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$ (dotted curves). This was achieved by comparison with the numerical integration results (solid curves for magenta for $x(t)$, red curve for $y(t)$, and blue for $z(t)$).

Figure 5 The OAFM solutions $({\overline{X}}_{OAFM}(y),{\overline{Z}}_{OAFM}(y))$ of System (39) for $x_{i0}=0.25$, $y_{i0}=0.5$, $z_{i0}=1.5$; for physical constants $a=0.56$, $b=0.84$, $c=0.275$, $d=2.165$, $N_{max}=10$. This was achieved by comparison with numerical integration of the OAFM solutions (dotted black curve) and numerical results (solid color curve: green line for $X(y)$ and red curve for $Z(y)$), respectively.

Figure 6 Variation in the function ${\epsilon }_X=|{\overline{X}}_{OAFM}(y)-X_{numerical}(y)|$ for $x_{i0}=0.25$, $y_{i0}=0.5$, $z_{i0}=1.5$; for physical constants $a=0.56$, $b=0.84$, $c=0.275$, $d=2.165$, $N_{max}=10$.

Figure 7 Variation in the function ${\epsilon }_Z=|{\overline{Z}}_{OAFM}(y)-Z_{numerical}(y)|$ for $x_{i0}=0.25$, $y_{i0}=0.5$, and $z_{i0}=1.5$; and physical constants $a=0.56$, $b=0.84$, $c=0.275$, $d=2.165$, and $N_{max}=10$.

Figure 8 The performance of the approximate solutions ${\overline{x}}_{OAFM}(t)$, ${\overline{y}}_{OAFM}(t)$, ${\overline{z}}_{OAFM}(t)$ of Equation (1) using Equation (A1) (dotted green curves). This was achieved by comparison with the iterative solutions $x_{iter}(t)$, $y_{iter}(t)$, $z_{iter}(t)$ using Equation (47) (dashing black curves) and the corresponding numerical ones (solid color curves for magenta for $x(t)$, red for $y(t)$, and blue for $z(t)$), respectively.

Table 1

The numerical function $x_{numerical}$ for $x_{i0}=0.25$ and $y_{i0}=0.5$, $z_{i0}=1.5$; physical constants $a=0.27$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$, index $N_{max}=15$; the ${\overline{x}}_{OAFM}$ solution from Equation (1) using Equations (28) and (A1); the absolute values ${ϵ}_x=|x_{numerical}-{\overline{x}}_{OAFM}|$.

t	${\mathit{x}}_{\mathit{numerical}}$	${\overline{\mathit{x}}}_{\mathit{OAFM}}$	${\mathit{ϵ}}_{\mathit{x}}$
0	0.25	0.25	2.8421 $\times {10}^{-14}$
6	−0.0282696503	−0.0283332544	6.3604 $\times {10}^{-5}$
12	0.0875446140	0.0876729147	1.2830 $\times {10}^{-4}$
18	−0.0080246353	−0.0081568880	1.3225 $\times {10}^{-4}$
24	−0.0258419393	−0.0257718485	7.0090 $\times {10}^{-5}$
30	0.0615665579	0.0615581458	8.4121 $\times {10}^{-6}$
36	−0.0546747900	−0.0547310459	5.6255 $\times {10}^{-5}$
42	0.0303556825	0.0303595022	3.8197 $\times {10}^{-6}$
48	−0.0091051729	−0.0090233939	8.1779 $\times {10}^{-5}$
54	−0.0064377130	−0.0064017731	3.5939 $\times {10}^{-5}$
60	0.0182054744	0.0180646680	1.4080 $\times {10}^{-4}$

Table 2

The numerical function $x_{numerical}$ for $x_{i0}=0.25$, $y_{i0}=0.5$, $z_{i0}=1.5$; for physical constants $a=0.49$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$, index $N_{max}=15$; and the ${\overline{x}}_{OAFM}$ solution from Equation (1) using Equations (38) and (A2), as well as the absolute values ${ϵ}_x=|x_{numerical}-{\overline{x}}_{OAFM}|$.

t	${\mathit{x}}_{\mathit{numerical}}$	${\overline{\mathit{x}}}_{\mathit{OAFM}}$	${\mathit{ϵ}}_{\mathit{x}}$
0	0.25	0.25	1.1102 $\times {10}^{-16}$
2	−0.2310125718	−0.2309915725	2.0999 $\times {10}^{-5}$
4	−0.1270907520	−0.1271941460	1.0339 $\times {10}^{-4}$
6	−0.0495725523	−0.0495150146	5.7537 $\times {10}^{-5}$
8	−0.0186868834	−0.0187247943	3.7910 $\times {10}^{-5}$
10	−0.0073059431	−0.0073250783	1.9135 $\times {10}^{-5}$
12	−0.0031036310	−0.0030636909	3.9940 $\times {10}^{-5}$
14	−0.0015389657	−0.0015240290	1.4936 $\times {10}^{-5}$
16	−0.0009541478	−0.0009876586	3.3510 $\times {10}^{-5}$
18	−0.0007359285	−0.0007665108	3.0582 $\times {10}^{-5}$
20	−0.0006552978	−0.0006176353	3.7662 $\times {10}^{-5}$

Table 3

Values of the OAFM solution ${\overline{y}}_{OAFM}(t)$ using (A1), the iterative solution $y_{iter}(t)$, and numerical integration.

t	${\mathit{y}}_{\mathit{numerical}}$	${\overline{\mathit{y}}}_{\mathit{OAFM}}$	${\mathit{y}}_{\mathit{iter}}$
0	0.5	0.5	0.5
0.35	0.5374226057	0.5374486596	0.5374226145
0.7	0.5438221770	0.5438941587	0.5438231287
1.05	0.5293352623	0.5294043272	0.5293606734
1.4	0.5017292155	0.5017640253	0.5019829926
1.75	0.4667123295	0.4667257564	0.4681816626
2.1	0.4283693903	0.4283893815	0.4344330377
2.45	0.3895448972	0.3895842146	0.4096806911
2.8	0.3521465461	0.3521957668	0.4109737727
3.15	0.3173816854	0.3174212857	0.4807833563
3.5	0.2859444082	0.2859602449	0.7463073937

Table 4

Values of the OAFM solution ${\overline{z}}_{OAFM}(t)$ using (A1), the iterative solution $z_{iter}(t)$, and numerical integration.

t	${\mathit{z}}_{\mathit{numerical}}$	${\overline{\mathit{z}}}_{\mathit{OAFM}}$	${\mathit{z}}_{\mathit{iter}}$
0	1.5	1.5	1.5
0.35	1.1377980798	1.1378008156	1.1377980375
0.7	0.8888423764	0.8888612686	0.8888369606
1.05	0.7116898089	0.7117270141	0.7115520959
1.4	0.5800303146	0.5800722098	0.5786930777
1.75	0.4778418539	0.4778775637	0.4702214720
2.1	0.3956568187	0.3956859302	0.3646682873
2.45	0.3279561047	0.3279838028	0.2279547171
2.8	0.2714770533	0.2715066870	−0.0023675260
3.15	0.2241809182	0.2242108232	−0.4447085587
3.5	0.1846601213	0.1846861393	−1.3396706314

Appendix A

Example A1. 

${\overline{x}}_{OAFM}$ is a semi-analytical solution for the system created by Equation (1) for initial conditions $x_{i0}=0.25$, $y_{i0}=0.5$, and $z_{i0}=1.5$; and physical constants $a=0.27$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$, and index $N_{max}=15$. The numerical values of the convergence-control parameters for ${\overline{x}}_{OAFM}$ were obtained from Equations (2), (22) and (27).

(A1) $$\begin{array}{ccc}{M}_0=-74.2469110791,\hfill & M_1=0.25,\hfill & M_2=-6.7436101551,\hfill \\ K_1=0.1848690496,\hfill & {\omega }_0=0.0803039394,\hfill & C_1=139.0054402824,\hfill \\ C_2=-22.8316572959,\hfill & C_3=-58.1056207469,\hfill & C_4=80.49563386601,\hfill \\ C_5=-2.6238448048,\hfill & C_6=-47.46330213287,\hfill & C_7=11.4210378291,\hfill \\ C_8=8.67200516746,\hfill & C_9=-6.4707783601,\hfill & C_{10}=0.88387227681,\hfill \\ C_{11}=1.7647286122,\hfill & C_{12}=-0.55303594910,\hfill & C_{13}=-0.2266744776,\hfill \\ C_{14}=0.08103920859,\hfill & C_{15}=0.0173679019,\hfill & D_1=-100.2318705179,\hfill \\ D_2=103.0314766276,\hfill & D_3=-68.2362545841,\hfill & D_4=-18.8861252680,\hfill \\ D_5=71.5941395127,\hfill & D_6=-9.3470883136,\hfill & D_7=-23.5202358023,\hfill \\ D_8=9.8452784472,\hfill & D_9=1.5769917873,\hfill & D_{10}=-3.5741473485,\hfill \\ D_{11}=0.9900318369,\hfill & D_{12}=0.6939046570,\hfill & D_{13}=-0.2438312600,\hfill \\ D_{14}=-0.0690954352,\hfill & D_{15}=0.0116149725.\hfill \end{array}$$ Example A2. 

${\overline{x}}_{OAFM}$ is a semi-analytical solution for the problem given by Equation (1) for the initial conditions $x_{i0}=0.25$, $y_{i0}=0.5$, and $z_{i0}=1.5$; and the physical constants $a=0.49$, $b=-(a-1/2)/(2d)$, $c=0.75$, $d=0.65$, and index $N_{max}=15$. The numerical values of the convergence-control parameters for ${\overline{x}}_{OAFM}$ were obtained from Equations (2), (32) and (37).

(A2) $$\begin{array}{ccc}{K}_1=0.4377235140,\hfill & K_2=0.4162288571,\hfill & K_s=0.4155127604,\hfill \\ D_1=1.1082468204,\hfill & D_2=-0.3028978305,\hfill & D_3=-0.1494396294,\hfill \\ C_1=0.0192987619,\hfill & D_4=-0.3344629683,\hfill & D_5=-0.0141957205,\hfill \\ D_6=2.4124142218,\hfill & C_2=0.1827701801,\hfill & D_7=-0.7737838521,\hfill \\ D_8=0.2548517753,\hfill & D_9=-2.4620720976,\hfill & C_3=-0.1767723559.\hfill \end{array}$$

References

1. Rucklidge, A.M. Chaos in models of double convection. J. Fluid Mech.; 1992; 237, pp. 209-229. [DOI: https://dx.doi.org/10.1017/S0022112092003392]

2. Chen, Z.; Wen, G.; Zhou, H.; Chen, J. A new M×N-grid double-scroll chaotic attractors from Rucklidge chaotic system. Optik; 2017; 136, pp. 27-35. [DOI: https://dx.doi.org/10.1016/j.ijleo.2017.01.088]

3. Lima, M.F.S.; Llibre, J.; Valls, C. Integrability of the Rucklidge system. Nonlinear Dyn.; 2014; 77, pp. 1441-1453. [DOI: https://dx.doi.org/10.1007/s11071-014-1389-y]

4. Hosen, M.Z.; Nurujjaman, M.; Tahura, S.; Ahmed, P. Controlling of chaotic Rucklidge system with dynamical behaviors. Seybold Report.; 2023; 18, pp. 849-865. [DOI: https://dx.doi.org/10.17605/OSF.IO/Ajuby]

5. Vignesh, D.; He, S.; Banerjee, S. Modelling discrete time fractional Rucklidge system with complex state variables and its synchronization. Appl. Math. Comput.; 2023; 455, 128111. [DOI: https://dx.doi.org/10.1016/j.amc.2023.128111]

6. Kocamaz, U.E.; Uyaroglu, Y. Controlling Rucklidge chaotic system with a single controller using linear feedback and passive control methods. Nonlinear Dyn.; 2014; 75, pp. 63-72. [DOI: https://dx.doi.org/10.1007/s11071-013-1049-7]

7. Dong, C.; Lian, J.; Qi, J.; Hantao, L. Symbolic Encoding of Periodic Orbits and Chaos in the Rucklidge System. Complexity; 2021; 2021, 4465151. [DOI: https://dx.doi.org/10.1155/2021/4465151]

8. Naveen, S.; Parthiban, V. Application of Newton’s polynomial interpolation scheme for variable order fractional derivative with power-law kernel. Sci. Rep.; 2024; 14, 16090. [DOI: https://dx.doi.org/10.1038/s41598-024-66494-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38997322]

9. Tarammim, A.; Akter, M. A Comparative Study of Synchronization Methods of Rucklidge Chaotic Systems with Design of Active Control and Backstepping Methods. Int. J. Mod. Theory Appl.; 2022; 11, pp. 31-51. [DOI: https://dx.doi.org/10.4236/ijmnta.2022.112003]

10. Wang, X. Si’lnikov chaos and Hopf bifurcation analysis of Rucklidge system. Chaos Solitons Fractals; 2009; 42, pp. 2208-2217. [DOI: https://dx.doi.org/10.1016/j.chaos.2009.03.137]

11. Marinca, V.; Herisanu, N. Approximate analytical solutions to Jerk equation. Springer Proceedings in Mathematics & Statistics, Proceedings of the Dynamical Systems: Theoretical and Experimental Analysis, Lodz, Poland, 7–10 December 2015; Springer: Cham, Switzerland, 2016; 182, pp. 169-176.

12. Marinca, V.; Ene, R.-D.; Marinca, V.B. Optimal Auxiliary Functions Method for viscous flow due to a stretching surface with partial slip. Open Eng.; 2018; 8, pp. 261-274. [DOI: https://dx.doi.org/10.1515/eng-2018-0028]

13. Ene, R.-D.; Pop, N.; Lapadat, M.; Dungan, L. Approximate closed-form solutions for the Maxwell-Bloch equations via the Optimal Homotopy Asymptotic Method. Mathematics; 2022; 10, 4118. [DOI: https://dx.doi.org/10.3390/math10214118]

14. Farooq, M.; Rahman, A.U.; Khan, A.; Ozsahin, I.; Uzun, B.; Ahmad, H. Comparative analysis of magnetohydrodynamic inclined Poiseulle flow of couple stress fluids. J. Comput. Appl. Mech.; 2025; 56, pp. 536-560.

15. Daftardar-Gejji, V.; Jafari, H. An iterative method for solving nonlinear functional equations. J. Math. Anal. Appl.; 2006; 316, pp. 753-763. [DOI: https://dx.doi.org/10.1016/j.jmaa.2005.05.009]