Abstract

Translate

The present paper focuses on some classes of dynamical systems involving Hamilton–Poisson structures, while neglecting their chaotic behaviors. Based on this, the closed-form solutions are obtained. These solutions are derived using the Optimal Auxiliary Functions Method (OAFM). The impact of the physical parameters of the system is also investigated. Periodic orbits around the equilibrium points are performed. There are homoclinic or heteroclinic orbits and they are obtained in exact form. The dynamical system is reduced to a second-order nonlinear differential equation, which is analytically solved through the OAFM procedure. The influence of initial conditions on the system is explored, specifically regarding the presence of symmetries. A good agreement between the analytical and corresponding numerical results is demonstrated, reflecting the accuracy of the proposed method. A comparative analysis underlines the advantages of the OAFM compared with the iterative method. The results of this work encourage the study of dynamical systems with bi-Hamiltonian structure and similar properties as physical and biological problems.

Full text

Turn on search term navigation

Translate

1. Introduction

The study of phenomena in nature is closely related to the theory of dynamical systems. This theory is an important part of mathematics, used to describe the behavior of complex dynamical systems, usually through the use of differential equations. Using these equation, the nature of certain dynamical systems and the solutions of the equations of motion are studied, in the context of planetary orbits or the behavior of electronic circuits, as well as systems that occur in engineering, biology, medicine, economics, secure communications, the optimization of the performance of the nonlinear system, and in other important fields. A variety of complex dynamic behaviors are generated by numerical analysis.

Three-dimensional systems of common differential equations provide a simple and effective means of capturing the complicated interaction between multiple entities. These systems exhibit rich dynamical properties and are further investigated with a variety of methods. In many situations, these systems admit Hamilton–Poisson structures for some specific parameter values, such as Euler’s equations in rigid body dynamics [1], the Maxwell–Bloch system [2] in laser physics, Rikitake-type models in geophysics [3], and the Rössler system in chaos theory [4].

The Poisson structure that is obtained for various dynamical systems has applications in physics and biology. In atomic physics, the dynamical system could be studied as a model describing the change in the energy levels of an atom with only two levels as a function of the eigenvalue [5,6].

The bi-Hamiltonian structure of the SIR model of Kermack and McKendrick could be governing the spread of epidemics for a closed population. The basic structure of the Kermack–McKendrick equations governing the spread of epidemics allows considerable freedom for incorporating new effects without changing the underlying bi-Hamiltonian structure.

Tudoran et al. in 2009 [7,8], after studying the two-disk Rikitake dynamical system, established the connections between the dynamical elements of this systems and the semialgebraic stratifications associated with the image of the energy Casimir mapping for the two-disk Rikitake system. Thus, using the approach of non-standard Poisson geometry, they showed that there is a connection between the dynamical properties of a given system and the geometry of the image of a vector-valued constant of motion, showing that many dynamical elements (e.g., equilibria, periodic orbits, homoclinic and heteroclinic connections, and bifurcation phenomena) result from the image of this mapping. The identification of the semi-algebraic divisions of the image of the energy Casimir mapping gives a topological classification of the fibers of the energy-Casimir mapping. The relationships between the dynamical properties of 3D systems and the image of the energy-Casimir mappings are usually complicated. It has often been shown that some findings obtained from certain studied models have failed to generalize to a wider framework. Therefore, it is important to obtain more information about the dynamics of 3D systems through the energy-Casimir mappings, allowing us to find out how the dynamics of the system are influenced by the geometry of the image of the energy-Casimir mappings. Currently, much information is given for many specific systems for classes of dynamical systems [3,4,6,9,10,11] in the context of investigations into the dynamical properties of these systems.

Other geometric properties of dynamical systems, such as Hamiltonian realization, equilibrium points, and integrable deformations, have been analyzed in [1,2,3,4,6,7,8,9,10,11,12,13,14,15]. The interaction between laser light and a material sample composed of two-level atoms is described by Maxwell’s equations of the electric field [2] and the Schrödinger equations for the probability amplitudes of atomic levels. Many systems of first-order differential equations that model certain processes are three-dimensional systems. Some of them have two constants of motion. Accordingly, these are Hamiltonian Poisson systems [2], whose constants of motion are the energy or Hamiltonian function H and a Casimir C of the corresponding Lie algebra. The dynamics of such systems occurs at the intersection of the common level sets of the Hamiltonian and the Casimir [10]. The orbits of the system are included in the intersection of the level sets H = constant and C = constant. For many three-dimensional Hamilton–Poisson systems, connections have been reported between the associated $(H, C)$ energy-Casimir mapping and some of their dynamical properties [1,8]. These connections depend on the image of the energy-Casimir mapping. They also depend on the partition of this image given by the images of the equilibrium points.

The dual space of a Lie algebra admits a canonical Poisson structure, namely, the Lie–Poisson structure. Quadratic Hamilton–Poisson systems on these structures form a natural framework for a variety of dynamical systems (especially in mathematical physics); prominent examples are the classical Euler equations for the rigid body, as well as their extensions and generalizations (see, for example [8,9]). In particular, quadratic (positive) semidefinite Hamilton–Poisson systems (on Lie–Poisson spaces) arise in the study of optimal control problems invariant on Lie groups (see, for example [4,11]). Quadratic Hamilton–Poisson systems on three-dimensional Lie–Poisson spaces of lower dimensions have been investigated [16]. Lie–Poisson spaces admit a global Casimir function, the normal form of a system being essentially determined by its equilibrium states and their behavior. The stability of equilibria has been investigated using an (extended) Casimir-energy method, the Routh–Hurwitz criterion or the Lyapunov function method, while instability usually arises from spectral instability.

Recent achievements have been obtained in the case of Hamilton–Poisson jerk systems, which exhibit chaotic behavior, starting from the study of stability and demonstrating the existence of periodic orbits around nonlinearly stable equilibrium points. These systems correspond to the case of families of equations of the harmonic oscillator, the mathematical pendulum, the Duffing oscillator, and other anharmonic oscillators. Bifurcations in the dynamics of jerk systems are also analyzed [6].

The dynamics of a three-dimensional Hamilton–Poisson system takes place at the intersection of the level sets given by the two constants of motion. Thus, for certain constants of motion, the orbits of the system can be described. The advantages of studying the image of the energy-Casimir mapping are observed in more and more works, exploring the stability of equilibrium points, the existence of periodic orbits around nonlinearly stable equilibria and the existence of homoclinic or heteroclinic orbits, which lead to the description of the behavior of 3D dynamical systems [10].

A class of dynamical systems involving a Hamilton–Poisson structure is written as follows [9]:

(1) $\{\begin{matrix} \dot{x} = a y + b y z \\ \dot{y} = c x z + d x \\ \dot{z} = - k x y \end{matrix},$

subject to initial conditions

(2) $\begin{matrix} x (0) = x_{0}, y (0) = y_{0}, z (0) = z_{0}, \end{matrix}$

with

a, b, c, d, k \in R, k \neq 0

System (1) generalizes the Euler equations for a free rigid body [17] or the rigid body with a free spinning rotor [18], the real Maxwell–Bloch Equations [19]. This system can also describe a particular case of the following systems: the Rikitake system [20], Rabinovich system [21], May–Leonard system [5], and Chen–Lee system [22].

The structure of this paper is as follows: Section 1 introduces the topic, a class of dynamical systems, including two prime integrals (constants of motion), while Section 2 presents the equilibrium points, symmetries, periodic orbits, and homoclinic or heteroclinic orbits. In Section 3, we discuss the closed-form solutions of the studied dynamics. Section 4 describes the Optimal Auxiliary Functions Method (OAFM) and its application in deriving semi-analytical solutions. Section 5 includes numerical results and highlights the advantages of the applied method. Finally, the concluding remarks are presented in Section 6.

2. Symmetries, Equilibrium Points, Periodic Orbits,Homoclinic Orbits, Heteroclinic Orbits

There are two cases: $a c - b d \neq 0$ and $a c - b d = 0$ containing three subcases $b = 0$ , $c = 0$ and $b c \neq 0$ .

2.1. Case 1— $a c - b d \neq 0$

The Hamilton–Poisson structure is characterize by two prime integrals

(3) $\{\begin{matrix} H_{1} (x, y, z) = - \frac{d}{2} x^{2} + \frac{a}{2} y^{2} + \frac{a c - b d}{2 k} z^{2} \\ C_{1} (x, y, z) = \frac{c k}{2} x^{2} - \frac{b k}{2} y^{2} + (a c - b d) z \end{matrix},$

with

H_{1}

being the Hamiltonian and

C_{1}

the Casimir of the system (1).

For $b = 0$ or $c = 0$ or $b c \neq 0$ , the system (1) admits symmetries with respect to the $O z$ -axis, being invariant to transformation $(x, y, z) \mapsto (- x, - y, z)$ ;

The matrix of linear part of system (1) is

(4) $A = J (x, y, z) = (\begin{matrix} 0 & a + b z & b y \\ c z + d & 0 & c x \\ - k y & - k x & 0 \end{matrix}) .$

A. For $b = 0$ , system (1) becomes as follows:

(5) $\{\begin{matrix} \dot{x} = a y \\ \dot{y} = c x z + d x \\ \dot{z} = - k x y \end{matrix},$

with the following equilibrium points:

e_{1}^{M} = (M, 0, - \frac{d}{c})

M \neq 0

and

e_{2}^{M} = (0, 0, M)

The characteristic polynomial of $J (e_{1}^{M})$ is

(6) $\begin{matrix} P_{J (e_{1}^{M})} (λ) = - λ (λ^{2} + c k M^{2}), \end{matrix}$

whose roots are

(7) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | M | \sqrt{c k} \in C, \end{matrix}$

for

c k > 0

(i.e., purely imaginary).

There is one periodic orbit of system (5) around the equilibrium point $e_{1}^{M}$ whose period is closed to $\frac{2 π}{| M | \sqrt{c k}}$ .

The characteristic polynomial of $J (e_{2}^{M})$ is

(8) $\begin{matrix} P_{J (e_{2}^{M})} (λ) = - λ (λ^{2} - a (c M + d)), \end{matrix}$

whose roots are purely imaginary

(9) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i \sqrt{- a (c M + d)} \in C, \end{matrix}$

for

- a (c M + d) > 0

Therefore, there exists one periodic orbit of system (5) around the equilibrium point $e_{2}^{M}$ whose period is closed to $\frac{2 π}{\sqrt{- a (c M + d)}}$ .

For the equilibrium point $e_{2}^{M} = (0, 0, M)$ , when $- a (c M + d) < 0$ , we assume that $a > 0$ and $c M + d > 0$ , $c > 0$ , $k > 0$ . Then, $i \sqrt{c k} \in C$ . The periodic orbit does not exist, but there are homoclinic orbits given by the intersection between the level sets $H_{1} (x, y, z) = H_{1} (0, 0, M)$ and $C_{1} (x, y, z) = C_{1} (0, 0, M)$ as follows:

(10) $\{\begin{matrix} \frac{k}{c} y^{2} + {(z + \frac{d}{c})}^{2} = {(M + \frac{d}{c})}^{2} \\ \frac{k}{2 a} x^{2} + z = M \end{matrix} .$

The solution of the system (5) can be obtained by using the following transformation

(11) $\{\begin{matrix} x = \pm \sqrt{\frac{2 a}{k}} \sqrt{M - z} \\ y = \sqrt{\frac{c}{k}} (M + \frac{d}{c}) \frac{2 U}{1 + U^{2}} \\ z = - \frac{d}{c} + (M + \frac{d}{c}) \frac{1 - U^{2}}{1 + U^{2}} \end{matrix} .$

Substituting Equation (11) into the first relation from Equation (5) yields the following differential equation

$\frac{\dot{U}}{U \sqrt{1 + U^{2}}} = \pm \sqrt{a (c M + d)},$

whose solution is

(12) $\begin{matrix} U (t) = \pm \frac{1}{sinh (\sqrt{a (c M + d)} t)} . \end{matrix}$

Combining Equations (11) and (12) we obtain four solutions of system (5) $(\pm x (t),$ $\pm y (t), z (t))$ with $lim_{t \to \pm \infty} (\pm x (t), \pm y (t), z (t)) = (0, 0, M)$ . Therefore, we obtain four homoclinic orbits, which are graphically depicted in Figure 1.

Related to the equilibrium point $e_{1}^{M} = (M, 0, - \frac{c}{d})$ , when $- a (c M + d) \neq 0$ , we assume that $a > 0$ and $c < 0$ , $k > 0$ . Then, $i \sqrt{c k} \in R$ . The periodic orbit does not exists, but there exist heteroclinic orbits given by the intersection between the level sets $H_{1} (x, y, z) = H_{1} (M, 0, - \frac{d}{c})$ and $C_{1} (x, y, z) = C_{1} (M, 0, - \frac{d}{c})$ as follows:

(13) $\{\begin{matrix} \frac{k}{2} y^{2} + \frac{c}{2} {(z + \frac{d}{c})}^{2} = 0 \\ \frac{c k}{2} x^{2} + a c z = \frac{c k}{2} M^{2} - a d \end{matrix} .$

The solution of the system (5) can be obtained by solving the previous system by

(14) $\{\begin{matrix} x = \pm \sqrt{M^{2} - \frac{2 a}{k} U} \\ y = \pm \sqrt{\frac{- c}{k}} U \\ z = U - \frac{d}{c} \end{matrix} .$

Substituting Equation (14) into the first relation from Equation (5) yields the following differential equation

$\pm \frac{\dot{U}}{U \sqrt{M^{2} - \frac{2 a}{k} U}} = \mp \sqrt{- c k},$

whose solution is

(15) $\begin{matrix} | \frac{V (t) + | M |}{V (t) - | M |} | = e^{\pm | M | \sqrt{- c k} t + C}, \end{matrix}$

where

V^{2} (t) = M^{2} - \frac{2 a}{k} U (t)

. Then, we consider

(16) $\begin{matrix} V_{1} (t) = | M | \frac{1 - e^{- p_{1} (t)}}{1 + e^{- p_{1} (t)}}, V_{2} (t) = - | M | \frac{1 + e^{- p_{1} (t)}}{1 - e^{- p_{1} (t)}}, \end{matrix}$

where

p_{1} (t) = | M | \sqrt{- c k} t + C

. Combining Equations (14) and (16), we obtain the following solutions of system (5)

(17) $\begin{matrix} x_{1} (t) = V_{1} (t), y_{1} (t) = \frac{\sqrt{- c k}}{2 a} (V_{1}^{2} (t) - M^{2}), z_{1} (t) = - \frac{d}{c} - \frac{k}{2 a} (V_{1}^{2} (t) - M^{2}), \\ x_{2} (t) = V_{2} (t), y_{2} (t) = \frac{\sqrt{- c k}}{2 a} (V_{2}^{2} (t) - M^{2}), z_{2} (t) = - \frac{d}{c} - \frac{k}{2 a} (V_{2}^{2} (t) - M^{2}) . \end{matrix}$

These solutions yield four heteroclinic orbits as

$(x_{1} (t), y_{1} (t), z_{1} (t)), (- x_{1} (t), y_{1} (t), z_{1} (t)), (x_{1} (t), - y_{1} (t), z_{1} (t)), (- x_{1} (t), - y_{1} (t), z_{1} (t)),$

and four split-heteroclinic orbits as

$(x_{2} (t), y_{2} (t), z_{2} (t)), (- x_{2} (t), y_{2} (t), z_{2} (t)), (x_{2} (t), - y_{2} (t), z_{2} (t)), (- x_{2} (t), - y_{2} (t), z_{2} (t)),$

of system (5), respectively, with the following:

$\begin{matrix} lim_{t \to \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to \infty} (x_{1} (t), - y_{1} (t), z_{1} (t)) = (M, 0, - \frac{d}{c}), \\ lim_{t \to - \infty} (- x_{1} (t), - y_{1} (t), z_{1} (t)) = lim_{t \to - \infty} (- x_{1} (t), y_{1} (t), z_{1} (t)) = (- M, 0, - \frac{d}{c}), \\ lim_{t \to \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to \infty} (x_{2} (t), - y_{2} (t), z_{2} (t)) = (M, 0, - \frac{d}{c}), \\ lim_{t \to - \infty} (- x_{2} (t), - y_{2} (t), z_{2} (t)) = lim_{t \to - \infty} (- x_{2} (t), y_{2} (t), z_{2} (t)) = (- M, 0, - \frac{d}{c}), \end{matrix}$

respectively.

These eight heteroclinic orbits are graphically depicted in Figure 2.

B. For $c = 0$ , the system (1) becomes as follows:

(18) $\{\begin{matrix} \dot{x} = a y + b y z \\ \dot{y} = d x \\ \dot{z} = - k x y \end{matrix} .$

Then, the equilibrium points are $e_{3}^{M} = (0, M, - \frac{a}{b})$ , $M \neq 0$ , and $e_{4}^{M} = (0, 0, M)$ .

The characteristic polynomial of $J (e_{3}^{M})$ is

(19) $\begin{matrix} P_{J (e_{3}^{M})} (λ) = - λ (λ^{2} + b k M^{2}), \end{matrix}$

whose roots are

(20) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | M | \sqrt{b k} \in C, \end{matrix}$

for

b k > 0

(i.e., purely imaginary).

Therefore, there exists one periodic orbit of system (18) around the equilibrium point $e_{3}^{M}$ whose period is closed to $\frac{2 π}{| M | \sqrt{b k}}$ .

The characteristic polynomial of $J (e_{4}^{M})$ is

(21) $\begin{matrix} P_{J (e_{4}^{M})} (λ) = - λ (λ^{2} - d (a + b M)), \end{matrix}$

whose roots are purely imaginary

(22) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i \sqrt{- d (a + b M)} \in C, \end{matrix}$

for

- d (a + b M) > 0

Therefore, there exists one periodic orbit of system (18) around the equilibrium point $e_{4}^{M}$ whose period is closed to $\frac{2 π}{\sqrt{- d (a + b M)}}$ .

Related to the equilibrium point $e_{4}^{M} = (0, 0, M)$ , when $- d (b M + a) < 0$ , we assume that $d < 0$ and $b M + a < 0$ , $b > 0$ , $k > 0$ . Then, $i \sqrt{b k} \in C$ . The periodic orbit does not exists, but there exists homoclinic orbits given by the intersection between the level sets $H_{1} (x, y, z) = H_{1} (0, 0, M)$ and $C_{1} (x, y, z) = C_{1} (0, 0, M)$ as follows:

(23) $\{\begin{matrix} \frac{k}{b} x^{2} + {(z + \frac{a}{b})}^{2} = {(M + \frac{a}{b})}^{2} \\ \frac{k}{2 d} y^{2} + z = M \end{matrix} .$

The solution of system (18) can be obtained by using the following transformation

(24) $\{\begin{matrix} x = \sqrt{\frac{b}{k}} (M + \frac{a}{b}) \frac{2 U}{1 + U^{2}} \\ y = \pm \sqrt{\frac{2}{k}} \sqrt{(- d) (z - M)} \\ z = - \frac{a}{b} + (M + \frac{a}{b}) \frac{1 - U^{2}}{1 + U^{2}} \end{matrix} .$

Substituting Equation (24) into the second relation from Equation (18) yields the following differential equation

$\frac{\dot{U}}{U \sqrt{1 + U^{2}}} = \pm \sqrt{d (b M + a)},$

whose solution is

(25) $\begin{matrix} U (t) = \pm \frac{1}{sinh (\sqrt{d (b M + a)} t)} . \end{matrix}$

Combining Equations (24) and (25) we obtain four solutions of system (18) $(\pm x (t), \pm y (t),$ $z (t))$ with $lim_{t \to \pm \infty} (\pm x (t), \pm y (t), z (t)) = (0, 0, M)$ . Therefore, we obtain four homoclinic orbits which are graphically presented in Figure 3.

Related to the equilibrium point $e_{3}^{M} = (0, M, - \frac{a}{b})$ , when $- d (b M + a) \neq 0$ , we assume that $a > 0$ , $b < 0$ , $d < 0$ , $k > 0$ . Then, $i \sqrt{b k} \in R$ . The periodic orbit does not exists, but there exist heteroclinic orbits given by the intersection between the level sets $H_{1} (x, y, z) = H_{1} (0, M, - \frac{a}{b})$ and $C_{1} (x, y, z) = C_{1} (0, M, - \frac{a}{b})$ as follows:

(26) $\{\begin{matrix} \frac{k}{2} x^{2} + \frac{b}{2} {(z + \frac{a}{b})}^{2} = 0 \\ \frac{- b k}{2} y^{2} + (- b) d z = \frac{- b k}{2} M^{2} + a d \end{matrix} .$

The solution of system (18) can be obtained by solving the previous system by

(27) $\{\begin{matrix} x = \pm \sqrt{\frac{- b}{k}} U \\ y = \pm \sqrt{M^{2} + \frac{- 2 d}{k} U} \\ z = U - \frac{a}{b} \end{matrix} .$

Substituting Equation (27) into the third relation from Equation (18) yields the following differential equation

$\pm \frac{\dot{U}}{U \sqrt{M^{2} + \frac{- 2 d}{k} U}} = \mp \sqrt{- b k},$

whose solution is

(28) $\begin{matrix} | \frac{V (t) - | M |}{V (t) + | M |} | = e^{\pm | M | \sqrt{- b k} t + C}, \end{matrix}$

where

V^{2} (t) = M^{2} + \frac{- 2 d}{k} U (t)

. Then, we consider

(29) $\begin{matrix} V_{1} (t) = | M | \frac{1 - e^{- p_{3} (t)}}{1 + e^{- p_{3} (t)}}, V_{2} (t) = - | M | \frac{1 + e^{- p_{3} (t)}}{1 - e^{- p_{3} (t)}}, \end{matrix}$

where

p_{3} (t) = | M | \sqrt{- b k} t + C

. Combining Equations (27) and (29), we obtain the following solutions of system (18)

(30) $\begin{matrix} x_{1} (t) = \frac{\sqrt{- b k}}{- 2 d} (V_{1}^{2} (t) - M^{2}), y_{1} (t) = V_{1} (t), z_{1} (t) = - \frac{a}{b} + \frac{k}{- 2 d} (V_{1}^{2} (t) - M^{2}), \\ x_{2} (t) = \frac{\sqrt{- b k}}{- 2 d} (V_{2}^{2} (t) - M^{2}), y_{2} (t) = V_{2} (t), z_{2} (t) = - \frac{a}{b} + \frac{k}{- 2 d} (V_{2}^{2} (t) - M^{2}) . \end{matrix}$

These solutions yield four heteroclinic orbits of the system (18) as follows:

$(x_{1} (t), y_{1} (t), z_{1} (t)), (- x_{1} (t), y_{1} (t), z_{1} (t)), (x_{1} (t), - y_{1} (t), z_{1} (t)), (- x_{1} (t), - y_{1} (t), z_{1} (t)),$

and four split-heteroclinic orbits

$(x_{2} (t), y_{2} (t), z_{2} (t)), (- x_{2} (t), y_{2} (t), z_{2} (t)), (x_{2} (t), - y_{2} (t), z_{2} (t)), (- x_{2} (t), - y_{2} (t), z_{2} (t)),$

with the following:

$\begin{matrix} lim_{t \to \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to \infty} (- x_{1} (t), y_{1} (t), z_{1} (t)) = (0, M, - \frac{a}{b}), \\ lim_{t \to - \infty} (- x_{1} (t), - y_{1} (t), z_{1} (t)) = lim_{t \to - \infty} (x_{1} (t), - y_{1} (t), z_{1} (t)) = (0, - M, - \frac{a}{b}), \\ lim_{t \to \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to \infty} (- x_{2} (t), y_{2} (t), z_{2} (t)) = (0, M, - \frac{a}{b}), \\ lim_{t \to - \infty} (- x_{2} (t), - y_{2} (t), z_{2} (t)) = lim_{t \to - \infty} (x_{2} (t), - y_{2} (t), z_{2} (t)) = (0, - M, - \frac{a}{b}), \end{matrix}$

respectively.

These eight heteroclinic orbits are graphically depicted in Figure 4.

C. For $b c \neq 0$ , the equilibrium points are $e_{5}^{M} = (0, M, - \frac{a}{b})$ , $M \neq 0$ , $e_{6}^{M} = (M, 0, - \frac{d}{c})$ , $M \neq 0$ , and $e_{7}^{M} = (0, 0, M)$ .

The characteristic polynomial of $J (e_{5}^{M})$ is

(31) $\begin{matrix} P_{J (e_{5}^{M})} (λ) = - λ (λ^{2} + b k M^{2}), \end{matrix}$

whose roots are

(32) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | M | \sqrt{b k} \in C, \end{matrix}$

for

b k > 0

(i.e., purely imaginary).

Therefore, there exists one periodic orbit of system (1) around the equilibrium point $e_{5}^{M}$ whose period is closed to $\frac{2 π}{| M | \sqrt{b k}}$ .

The characteristic polynomial of $J (e_{6}^{M})$ is

(33) $\begin{matrix} P_{J (e_{6}^{M})} (λ) = - λ (λ^{2} + c k M^{2}), \end{matrix}$

whose roots are purely imaginary

(34) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | M | \sqrt{c k} \in C, \end{matrix}$

for

c k > 0

Therefore, there exists one periodic orbit of system (1) around the equilibrium point $e^{M}$ whose period is closed to $\frac{2 π}{| M | \sqrt{c k}}$ .

The characteristic polynomial of $J (e_{7}^{M})$ is

(35) $\begin{matrix} P_{J (e_{7}^{M})} (λ) = - λ (λ^{2} - (a + b M) (c M + d)), \end{matrix}$

whose roots are purely imaginary

(36) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i \sqrt{- (a + b M) (c M + d)} \in C, \end{matrix}$

for

- (a + b M) (c M + d) > 0

Therefore, there exists one periodic orbit of system (1) around the equilibrium point $e_{7}^{M}$ whose period is closed to $\frac{2 π}{\sqrt{- (a + b M) (c M + d)}}$ .

In the case $a c - b d > 0$ with $b c > 0$ , the condition $- (a + b M) (c M + d) > 0$ is true if and only if $M \in (- \frac{a}{b}, - \frac{d}{c}) \cup (- \frac{d}{c}, - \frac{a}{b})$ .

In the case $a c - b d > 0$ with $b c < 0$ , the condition $- (a + b M) (c M + d) > 0$ is true if and only if $M \in (- \infty, - \frac{a}{b}) \cup (- \frac{d}{c}, \infty)$ or $M \in (- \infty, - \frac{d}{c}) \cup (- \frac{a}{b}, \infty)$ .

The case $a c - b d < 0$ is similar to the case $a d - b c > 0$ .

In relation to the equilibrium point $e_{7}^{M} = (0, 0, M)$ , when $- (b M + a) (c M + d) < 0$ , we assume that $d < 0$ , $b > 0$ , $c > 0$ and $k > 0$ . Then, $i \sqrt{- (b M + a) (c M + d)} \in R$ . The periodic orbit does not exist, but there exists homoclinic orbits given by the intersection between the level sets $H_{1} (x, y, z) = H_{1} (0, 0, M)$ and $C_{1} (x, y, z) = C_{1} (0, 0, M)$ as follows:

(37) $\{\begin{matrix} \frac{k}{b} x^{2} + {(z + \frac{a}{b})}^{2} = {(M + \frac{a}{b})}^{2} \\ \frac{k}{c} y^{2} + {(z + \frac{d}{c})}^{2} = {(M + \frac{d}{c})}^{2} \end{matrix} .$

The solution of system (1) can be obtained by using the following transformation

(38) $\{\begin{matrix} x = \pm \sqrt{\frac{b}{k}} \sqrt{{(M + \frac{a}{b})}^{2} - {(z + \frac{a}{b})}^{2}} \\ y = \pm \sqrt{\frac{c}{k}} \sqrt{{(M + \frac{d}{c})}^{2} - {(z + \frac{d}{c})}^{2}} \\ z = | \frac{a}{b} - \frac{d}{c} | \frac{1 + U^{2}}{2 U} - M - (\frac{a}{b} + \frac{d}{c}) \end{matrix} .$

Substituting Equation (38) into the third relation from Equation (1) yields the following differential equation

$\frac{2 \dot{U}}{2 U (2 M + \frac{a}{b} + \frac{d}{c}) - (\frac{a}{b} - \frac{d}{c}) (1 + U^{2})} = \pm \sqrt{b c},$

whose solutions are

(39) $\begin{matrix} U_{1} (t) = A_{0} + B_{0} \frac{1 - e^{- p_{7} (t)}}{1 + e^{- p_{7} (t)}}, U_{2} (t) = A_{0} + B_{0} \frac{1 + e^{- p_{7} (t)}}{1 - e^{- p_{7} (t)}}, \end{matrix}$

where

p_{7} (t) = 2 \sqrt{(b M + a) (c M + d)} t + C

A_{0} = \frac{2 M + \frac{a}{b} + \frac{d}{c}}{| \frac{a}{b} - \frac{d}{c} |}

B_{0} = 2 \frac{\sqrt{(M + \frac{a}{b}) (M + \frac{d}{c})}}{| \frac{a}{b} - \frac{d}{c} |}

Combining Equations (38) and (39), for $U_{1} (t)$ , we obtain four solutions of system (1) $(x_{1} (t), y_{1} (t), z_{1} (t))$ , $(- x_{1} (t), - y_{1} (t), z_{1} (t))$ , $(- x_{1} (t), y_{1} (t), z_{1} (t))$ , $(x_{1} (t), - y_{1} (t), z_{1} (t))$ for $M > - d / c$ and $M < - 1$ , respectively, and two solutions of system (1) $(- x_{1} (t), - y_{1} (t), z_{1} (t))$ , $(x_{1} (t), - y_{1} (t), z_{1} (t))$ , for $- 1 < M < - a / b$ with $lim_{t \to \pm \infty} (\pm x_{1} (t), \pm y_{1} (t), z_{1} (t)) = (0, 0, M)$ . These homoclinic orbits are graphically shown in Figure 5 and Figure 6, respectively.

2.2. Case 2— $a c - b d = 0$

The system (1) becomes as follows:

(40) $\{\begin{matrix} \dot{x} = a y + b y z \\ \dot{y} = c x z + \frac{a c}{b} x \\ \dot{z} = - k x y \end{matrix}, b \neq 0 .$

In this case, the Hamiltonian $H_{2}$ and the Casimir $C_{2}$ of the system (40) are

(41) $\{\begin{matrix} H_{2} (x, y, z) = \frac{k}{2} x^{2} + \frac{b}{2} {(z + \frac{a}{b})}^{2} \\ C_{2} (x, y, z) = \frac{c}{2} x^{2} - \frac{b}{2} y^{2} \end{matrix} .$

By means of these prime integrals (constants of motion), the semi-analytical solutions in closed-form could be obtained using the OAFM procedure [23]. This solution will be called the OAFM-solution and is denoted by $\bar{u}$ or ${\bar{u}}_{O A F M}$ . The accuracy is validated in Section 4.2.

For $a c - b d = 0$ , the system (40) admits symmetries with respect to the Oz-axis, being invariant to transformation $(x, y, z) \to (- x, - y, z)$ .

Assume that $d = \frac{a c}{b}$ . Then, the equilibrium points are as follows: $e_{8}^{M} = (0, M, - \frac{a}{b})$ , $M \neq 0$ , $e_{9}^{M} = (M, 0, - \frac{a}{b})$ , $M \neq 0$ , and $e_{10}^{M} = (0, 0, M)$ .

The characteristic polynomial of $J (e_{8}^{M})$ is

(42) $\begin{matrix} P_{J (e_{8}^{M})} (λ) = - λ (λ^{2} + b k M^{2}), \end{matrix}$

whose roots are

(43) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | M | \sqrt{b k} \in C, \end{matrix}$

for

b k > 0

(i.e., purely imaginary).

Therefore, there exists one periodic orbit of system (40) around the equilibrium point $e_{8}^{M}$ whose period is closed to $\frac{2 π}{| M | \sqrt{b k}}$ .

The characteristic polynomial of $J (e_{9}^{M})$ is

(44) $\begin{matrix} P_{J (e_{9}^{M})} (λ) = - λ (λ^{2} + c k M^{2}), \end{matrix}$

whose roots are purely imaginary

(45) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | M | \sqrt{c k} \in C, \end{matrix}$

for

c k > 0

Therefore, there exists one periodic orbit of system (40) around the equilibrium point $e_{9}^{M}$ whose period is closed to $\frac{2 π}{| M | \sqrt{c k}}$ .

The characteristic polynomial of $J (e_{10}^{M})$ is

(46) $\begin{matrix} P_{J (e_{10}^{M})} (λ) = - λ (λ^{2} - \frac{c}{b} {(a + b M)}^{2}), \end{matrix}$

whose roots are purely imaginary

(47) $\begin{matrix} λ_{1} = 0, λ_{2, 3} = \pm i | a + b M | \sqrt{- \frac{c}{b}} \in C, \end{matrix}$

for

- \frac{c}{b} > 0

and

M \neq - \frac{a}{b}

Therefore, there exists one periodic orbit of system (40) around the equilibrium point $e_{10}^{M}$ whose period is closed to $\frac{2 π}{| a + b M | \sqrt{- \frac{c}{b}}}$ .

Related to the equilibrium point $e_{8}^{M} = (0, M, - a / b)$ , when $b < 0$ , we assume that $a > 0$ , $c > 0$ , and $k > 0$ . Then, $i \sqrt{- b k} \in C$ . The periodic orbit does not exist, but there exist eight heteroclinic orbits given by the intersection between the level sets $H_{2} (x, y, z) = H_{2} (0, M, - a / b)$ and $C_{2} (x, y, z) = C_{2} (0, M, - a / b)$ as follows:

(48) $\{\begin{matrix} \frac{k}{2} x^{2} - \frac{- b}{2} {(z + \frac{a}{b})}^{2} = 0 \\ \frac{c}{2} x^{2} - \frac{b}{2} y^{2} = - \frac{b}{2} M^{2} \end{matrix} .$

The solution of the system (40) can be obtained by solving the previous system by

(49) $\{\begin{matrix} x = \pm \sqrt{\frac{- b}{c}} M \frac{1 - U^{2}}{1 + U^{2}} \\ y = \pm M \frac{2 U}{1 + U^{2}} \\ z = - \frac{a}{b} + \sqrt{\frac{k}{c}} M \frac{1 - U^{2}}{1 + U^{2}} \end{matrix} .$

Substituting Equation (49) into the third relation from Equation (40) yields the following differential equation

$\frac{2 \dot{U}}{U^{2} - 1} = \pm \sqrt{- b k} M,$

whose solutions are

(50) $\begin{matrix} U_{1} (t) = \frac{1 - C e^{- p_{8} (t)}}{1 + C e^{- p_{8} (t)}}, U_{2} (t) = \frac{1 + C e^{- p_{8} (t)}}{1 - C e^{- p_{8} (t)}}, \end{matrix}$

where

p_{8} (t) = \sqrt{- b k} M t

C = | \frac{1 - U (0)}{1 + U (0)} |

, and

U (0) = - \sqrt{\frac{z (0) + \frac{a}{b} + M \sqrt{\frac{k}{c}}}{M \sqrt{\frac{k}{c}} - z (0) - \frac{a}{b}}}

. Combining Equations (49) and (50), we obtain the following solutions of system (40)

(51) $\begin{matrix} x_{1} (t) = \sqrt{\frac{- b}{c}} M \frac{1 - U_{1}^{2}}{1 + U_{1}^{2}}, y_{1} (t) = M \frac{2 U_{1}}{1 + U_{1}^{2}}, z_{1} (t) = - \frac{a}{b} + \sqrt{\frac{k}{c}} M \frac{1 - U_{1}^{2}}{1 + U_{1}^{2}}, \\ x_{2} (t) = \sqrt{\frac{- b}{c}} M \frac{1 - U_{2}^{2}}{1 + U_{2}^{2}}, y_{2} (t) = M \frac{2 U_{2}}{1 + U_{2}^{2}}, z_{2} (t) = - \frac{a}{b} + \sqrt{\frac{k}{c}} M \frac{1 - U_{2}^{2}}{1 + U_{2}^{2}} . \end{matrix}$

These solutions describe only four heteroclinic orbits of the system (40) as follows:

$(x_{1} (t), y_{1} (t), z_{1} (t)), (- x_{1} (t), y_{1} (t), z_{1} (t)), (x_{2} (t), y_{2} (t), z_{2} (t)), (- x_{2} (t), y_{2} (t), z_{2} (t)),$

with

$lim_{t \to \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to \infty} (- x_{1} (t), y_{1} (t), z_{1} (t)) = (0, M, - a / b),$

$lim_{t \to \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to \infty} (- x_{2} (t), y_{2} (t), z_{2} (t)) = (0, M, - a / b),$

$lim_{t \to - \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to - \infty} (- x_{1} (t), y_{1} (t), z_{1} (t)) = (0, - M, - a / b),$

$lim_{t \to - \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to - \infty} (- x_{2} (t), y_{2} (t), z_{2} (t)) = (0, - M, - a / b) .$

Therefore, we obtain only four heteroclinic orbits which are graphically depicted in Figure 7.

Related to the equilibrium point $e_{9}^{M} = (M, 0, - a / b)$ , when $c < 0$ , we assume that $a > 0$ and $b > 0$ , $k > 0$ . Then, $i \sqrt{- c k} \in C$ . The periodic orbit does not exist, but there exist two heteroclinic orbits given by the intersection between the level sets $H_{2} (x, y, z) = H_{2} (M, 0, - a / b)$ and $C_{2} (x, y, z) = C_{2} (M, 0, - a / b)$ as follows:

(52) $\{\begin{matrix} \frac{k}{2} x^{2} + \frac{b}{2} {(z + \frac{a}{b})}^{2} = \frac{b}{2} M^{2} \\ \frac{c}{2} x^{2} - \frac{b}{2} y^{2} = \frac{c}{2} M^{2} \end{matrix} .$

The solution of the system (40) can be obtained by solving the previous system by

(53) $\{\begin{matrix} x = \pm M \frac{2 U}{1 + U^{2}} \\ y = \pm \sqrt{\frac{- c}{b}} M \frac{1 - U^{2}}{1 + U^{2}} \\ z = - \frac{a}{b} + \sqrt{\frac{k}{b}} M \frac{1 - U^{2}}{1 + U^{2}} \end{matrix} .$

Substituting Equation (53) into the first relation from Equation (40) yields the following differential equation

$\frac{2 \dot{U}}{U^{2} - 1} = \pm \sqrt{- c k} M,$

whose solutions are

(54) $\begin{matrix} U_{1} (t) = \frac{1 - C e^{- p_{9} (t)}}{1 + C e^{- p_{9} (t)}}, U_{2} (t) = \frac{1 + C e^{- p_{9} (t)}}{1 - C e^{- p_{9} (t)}}, \end{matrix}$

where

p_{9} (t) = \sqrt{- c k} M t

C = | \frac{1 - U (0)}{1 + U (0)} |

, and

U (0) = - \sqrt{\frac{z (0) + \frac{a}{b} + M \sqrt{\frac{k}{b}}}{M \sqrt{\frac{k}{b}} - z (0) - \frac{a}{b}}}

Combining Equations (53) and (54), we obtain the following solutions of system (40)

(55) $\begin{matrix} x_{1} (t) = M \frac{2 U_{1}}{1 + U_{1}^{2}}, y_{1} (t) = \sqrt{\frac{- c}{b}} M \frac{1 - U_{1}^{2}}{1 + U_{1}^{2}}, z_{1} (t) = - \frac{a}{b} + \sqrt{\frac{k}{b}} M \frac{1 - U_{1}^{2}}{1 + U_{1}^{2}}, \\ x_{2} (t) = M \frac{2 U_{2}}{1 + U_{2}^{2}}, y_{2} (t) = \sqrt{\frac{- c}{b}} M \frac{1 - U_{2}^{2}}{1 + U_{2}^{2}}, z_{2} (t) = - \frac{a}{b} + \sqrt{\frac{k}{b}} M \frac{1 - U_{2}^{2}}{1 + U_{2}^{2}} . \end{matrix}$

These solutions describe only four heteroclinic orbits of the system (40) as follows:

$(x_{1} (t), y_{1} (t), z_{1} (t)), (x_{1} (t), - y_{1} (t), z_{1} (t)), (x_{2} (t), y_{2} (t), z_{2} (t)), (x_{2} (t), - y_{2} (t), z_{2} (t)),$

with

$lim_{t \to \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to \infty} (x_{1} (t), - y_{1} (t), z_{1} (t)) = (M, 0, - a / b),$

$lim_{t \to \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to \infty} (x_{2} (t), - y_{2} (t), z_{2} (t)) = (M, 0, - a / b),$

$lim_{t \to - \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to - \infty} (x_{1} (t), - y_{1} (t), z_{1} (t)) = (- M, 0, - a / b),$

$lim_{t \to - \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to - \infty} (x_{2} (t), - y_{2} (t), z_{2} (t)) = (- M, 0, - a / b) .$

Therefore, we obtain only four heteroclinic orbits, which are graphically depicted in Figure 8.

Related to the equilibrium point $e_{10}^{M} = (0, 0, M)$ , when $b M + a \neq 0$ and $- c / b < 0$ , we assume that $a > 0$ and $b > 0$ , $c > 0$ , $k > 0$ . Then, $i \sqrt{- c / b} \in R$ . The periodic orbit does not exist, but there exist four heteroclinic orbits given by the intersection between the level sets $H_{2} (x, y, z) = H_{2} (0, 0, M)$ and $C_{2} (x, y, z) = C_{2} (0, 0, M)$ as follows:

(56) $\{\begin{matrix} \frac{k}{2} x^{2} + \frac{b}{2} {(z + \frac{a}{b})}^{2} = \frac{b}{2} {(M + \frac{a}{b})}^{2} \\ \frac{c}{2} x^{2} - \frac{b}{2} y^{2} = 0 \end{matrix} .$

The solution of the system (40) can be obtained by solving the previous system by

(57) $\{\begin{matrix} x = \pm \sqrt{\frac{b}{k}} (M + \frac{a}{b}) \frac{2 U}{1 + U^{2}} \\ y = \pm \sqrt{\frac{c}{b}} x \\ z = - \frac{a}{b} + (M + \frac{a}{b}) \frac{1 - U^{2}}{1 + U^{2}} \end{matrix} .$

Substituting Equation (57) into the first relation from Equation (40) yields the following differential equation

$\frac{\dot{U}}{U} = \pm \sqrt{b c} (M + \frac{a}{b}),$

whose solutions are

(58) $\begin{matrix} U_{1} (t) = e^{- p_{10} (t)}, U_{2} (t) = e^{p_{10} (t)}, \end{matrix}$

where

p_{10} (t) = \sqrt{b c} (M + \frac{a}{b}) t + C

. Combining Equations (57) and (58), we obtain the following solutions of system (40)

(59) $\begin{matrix} x_{1} (t) = \sqrt{\frac{b}{k}} (M + \frac{a}{b}) \frac{2 U_{1}}{1 + U_{1}^{2}}, x_{2} (t) = \sqrt{\frac{b}{k}} (M + \frac{a}{b}) \frac{2 U_{2}}{1 + U_{2}^{2}}, \\ y_{1} (t) = \sqrt{\frac{c}{b}} x_{1}, y_{2} (t) = \sqrt{\frac{c}{b}} x_{2}, \\ z_{1} (t) = - \frac{a}{b} + (M + \frac{a}{b}) \frac{1 - U_{1}^{2}}{1 + U_{1}^{2}}, z_{2} (t) = - \frac{a}{b} + (M + \frac{a}{b}) \frac{1 - U_{2}^{2}}{1 + U_{2}^{2}} . \end{matrix}$

These solutions describe four heteroclinic orbits of the system (40) as follows:

$(x_{1} (t), y_{1} (t), z_{1} (t)), (- x_{1} (t), y_{1} (t), z_{1} (t)), (x_{2} (t), y_{2} (t), z_{2} (t)), (- x_{2} (t), y_{2} (t), z_{2} (t)),$

with

$lim_{t \to \infty} (x_{1} (t), y_{1} (t), z_{1} (t)) = lim_{t \to \infty} (- x_{1} (t), y_{1} (t), z_{1} (t)) = (0, 0, M),$

$lim_{t \to - \infty} (x_{2} (t), y_{2} (t), z_{2} (t)) = lim_{t \to - \infty} (- x_{2} (t), y_{2} (t), z_{2} (t)) = (0, 0, - M - \frac{2 a}{b}) .$

Therefore, we obtain four heteroclinic orbits, which are graphically depicted in Figure 9.

Remark 1.

For $b < 0$ , $c < 0$ , and $d = \frac{a c}{b}$ , the system (40) admits periodic solutions if and only if $k < 0$ .

Remark 2.

All the analyzed cases lead to the periodic behaviors of the system described by Equation (1). The periodic solutions will be numerically computed in Section 4.2 using the OAFM technique for several values of the physical parameters a, b, c, d, $k > 0$ and initial conditions $x_{0} > 0$ , $y_{0} > 0$ , $z_{0} > 0$ , respectively.

3. Closed-Form Solutions

In this section, closed-form solutions for the system (1) are obtained using the prime integrals presented in the previous section for each relevant case for our study.

In the case, $a c - b d \neq 0$ , $c > 0$ , $k > 0$ with $b = 0$ , considering the second prime integral (3) for $z (t)$ and the first equation from (5) for $y (t)$ .

The closed-form solutions are written as

(60) $\{\begin{matrix} x = u \\ y = \frac{1}{a} \dot{u} \\ z = - \frac{k}{2 a} u^{2} + (z_{0} + \frac{k}{2 a} x_{0}^{2}) \end{matrix},$

where the unknown smooth function

u (t)

is the solution of the following problem

(61) $\{\begin{matrix} \frac{1}{a} \ddot{u} + \frac{c k}{2 a} u^{3} - c (z_{0} + \frac{k}{2 a} x_{0}^{2}) u - d u = 0 \\ u (0) = x_{0}, \dot{u} (0) = a y_{0} \end{matrix} .$

Considering $a c - b d \neq 0$ , $b > 0$ , $k > 0$ for $c = 0$ and taking into account the second prime integral (3) for $z (t)$ , respectively, the second equation from (18) for $x (t)$ , the closed-form solutions are written as

(62) $\{\begin{matrix} x = \frac{1}{d} \dot{u} \\ y = u \\ z = - \frac{k}{2 d} u^{2} + (z_{0} + \frac{k}{2 d} y_{0}^{2}) \end{matrix},$

where the unknown smooth function

u (t)

is the solution of the following problem

(63) $\{\begin{matrix} \frac{1}{d} \ddot{u} + \frac{b k}{2 d} u^{3} - b (z_{0} + \frac{k}{2 d} y_{0}^{2}) u - a u = 0 \\ u (0) = y_{0}, \dot{u} (0) = d x_{0} \end{matrix} .$

In the case $a c - b d \neq 0$ , $k > 0$ , $b > 0$ and $c > 0$ , but $a < 0$ and $d < 0$ , using the second prime integral (3) for $y (t), z (t)$ , respectively, the last equation from (1) for $x (t)$ , the closed-form solutions are defined by the following:

(64) $\{\begin{matrix} x = \frac{1}{\sqrt{c k}} \frac{2 \dot{u}}{1 + u^{2}} \\ y = \frac{R}{\sqrt{k}} \frac{2 u}{1 + u^{2}} \\ z = - \frac{d}{c} + \frac{R}{\sqrt{c}} \frac{1 - u^{2}}{1 + u^{2}} \end{matrix},$

where the unknown smooth function

u (t)

is the solution of the following problem

(65) $\{\begin{matrix} \ddot{u} (1 + u^{2}) - 2 u {(\dot{u})}^{2} - \frac{R}{\sqrt{c}} (a c - b d) u (1 + u^{2}) - b R^{2} u (1 - u^{2}) = 0 \\ u (0) = \frac{\sqrt{k} | y_{0} |}{| R + \sqrt{c} (z_{0} + \frac{d}{c}) |}, \dot{u} (0) = \sqrt{c k} \frac{x_{0}}{2} (1 + u^{2} (0)) \end{matrix},$

with

R = \sqrt{k y_{0}^{2} + c {(z_{0} + \frac{d}{c})}^{2}}

(in the same manner for (

a < 0

d > 0

); (

a > 0

d > 0

) and all subcases (

b < 0

c > 0

) with

a, d \in R

For $a c - b d \neq 0$ , $k > 0$ , $b > 0$ , and $c > 0$ , but $a > 0$ and $d < 0$ , the closed-form solutions are obtained from the second prime integral (3) for $x (t), z (t)$ and the last equation from (1) for $y (t)$ , as follows:

(66) $\{\begin{matrix} x = \frac{R}{\sqrt{k}} \frac{2 u}{1 + u^{2}} \\ y = \frac{1}{\sqrt{b k}} \frac{2 \dot{u}}{1 + u^{2}} \\ z = - \frac{a}{b} + \frac{R}{\sqrt{b}} \frac{1 - u^{2}}{1 + u^{2}} \end{matrix},$

where the unknown smooth function

u (t)

is the solution of the following problem

(67) $\{\begin{matrix} \ddot{u} (1 + u^{2}) - 2 u {(\dot{u})}^{2} - \frac{R}{\sqrt{b}} (b d - a c) u (1 + u^{2}) - c R^{2} u (1 - u^{2}) = 0 \\ u (0) = \frac{\sqrt{k} | x_{0} |}{| R + \sqrt{b} (z_{0} + \frac{a}{b}) |}, \dot{u} (0) = \sqrt{b k} \frac{y_{0}}{2} (1 + u^{2} (0)) \end{matrix},$

with

R = \sqrt{k x_{0}^{2} + b {(z_{0} + \frac{a}{b})}^{2}}

(in the same manner for all subcases (

b > 0

c < 0

) with

a, d \in R

In the last case $a c - b d = 0$ , $k > 0$ , $a > 0$ , $b > 0$ , $c > 0$ with $d = \frac{a c}{b}$ , the closed-form solutions resulted from the second prime integral (41) for $y (t), z (t)$ , and the last equation from (40) for $x (t)$ as follows:

The closed-form solutions are written as

(68) $\{\begin{matrix} x = \frac{1}{\sqrt{c k}} \frac{2 \dot{u}}{1 + u^{2}} \\ y = \frac{R}{\sqrt{k}} \frac{2 u}{1 + u^{2}} \\ z = - \frac{a}{b} + \frac{R}{\sqrt{c}} \frac{1 - u^{2}}{1 + u^{2}} \end{matrix},$

where the unknown smooth function

u (t)

is the solution of the following problem

(69) $\{\begin{matrix} \ddot{u} (1 + u^{2}) - 2 u {(\dot{u})}^{2} - b R^{2} u (1 - u^{2}) = 0 \\ u (0) = \frac{\sqrt{k} | y_{0} |}{| R + \sqrt{c} (z_{0} + \frac{a}{b}) |}, \dot{u} (0) = \sqrt{c k} \frac{x_{0}}{2} (1 + u^{2} (0)) \end{matrix},$

with

R = \sqrt{k y_{0}^{2} + c {(z_{0} + \frac{a}{b})}^{2}}

In Section 4, the above Equations (61), (63), (65), (67), and (69) are semi-analytically solved using the OAFM technique.

4. The Optimal Auxiliary Functions Method (OAFM)

4.1. Main Steps of the OAFM

We introduce the basic ideas of the OAFM by considering the second order nonlinear differential equation in general form, as in [12,24,25], as follows:

(70) $\begin{matrix} L [\bar{v} (t)] + g (t) + N [\bar{v} (t)] = 0, \\ B (\bar{v} (0), \dot{\bar{v}} (0)) = 0, \end{matrix}$

where

L

is a linear operator, g is a known function, and

N

is a given nonlinear operator. t denotes the independent variable and the approximate solution

\bar{v} (t)

is written with just two components in the following form:

(71) $\bar{v} (t) = v_{0} (t) + v_{1} (t, C_{i}), i = 1, 2, \dots, s .$

Remark 3.

The linear operator $L$ is chosen such that the initial approximation $v_{0} (t)$ becomes an elementary function or combination of the elementary functions. These functions could be as follows: the exponential function $e^{- K t}$ or the rational function $\frac{1}{1 + K t}$ in the case of the boundary value problems from fluid mechanics; trigonometric functions $c o s (ω_{0} t)$ , $s i n (ω_{0} t)$ , describing the nonlinear vibrations with periodic behaviors; $e^{- K t} c o s (ω_{0} t)$ , $e^{- K t} s i n (ω_{0} t)$ modeling the nonlinear vibrations with harmonic/anharmonic oscillations—damping effect.

The initial approximation $v_{0} (t)$ and the first approximation $v_{1} (t, C_{i})$ will be determined as follows.

Firstly, the initial approximation $v_{0} (t)$ is chosen as solution of the following equation:

(72) $\begin{matrix} L [v_{0} (t)] + g (t) = 0, \\ B (v_{0} (0), {\dot{v}}_{0} (0)) = 0 . \end{matrix}$

Now, Equation (70) becomes

(73) $\begin{matrix} L [v_{0} (t)] + L [v_{1} (t, C_{i})] + g (t) + N [v_{0} (t) + v_{1} (t, C_{i})] = 0 . \end{matrix}$

The nonlinear operator can be expanded in the form

(74) $\begin{matrix} N [v_{0} (t) + v_{1} (t, C_{i})] = N [v_{0} (t)] + \sum_{k = 1}^{\infty} \frac{v_{1}^{k} (t, C_{i})}{k!} N^{(k)} [v_{0} (t)] . \end{matrix}$

Using Equations (73) and (74), the first approximation $v_{1} (t)$ is the solution of the following problem:

(75) $\begin{matrix} L [v_{1} (t, C_{i})] + A_{1} [v_{0} (t), C_{i}] N [v_{0} (t)] + A_{2} [v_{0} (t), C_{j}] = 0, \end{matrix}$

(76) $v_{1} (0, C_{i}) = 0, {\dot{v}}_{1} (0, C_{i}) = 0,$

where

A_{1}

and

A_{2}

are two arbitrary auxiliary functions depending on the initial approximation

v_{0} (t)

and several unknown parameters

C_{i}

and

C_{j}

i = 1, 2, \dots, p

j = p + 1, p + 2, \dots, s

Finally, the approximate analytic solution is obtained from the Equation (71) in the following form:

(77) $\bar{v} (t) = v_{0} (t) + v_{1} (t, C_{i}), i = 1, 2, \dots, s,$

with

v_{0} (t)

and

v_{1} (t, C_{i})

solutions of Equations (72) and (75), respectively.

In the OAFM, the linear operator $L$ is arbitrarily chosen, not the physical parameters. There are situations when the choice of the physical parameters gives rise to the chaotic behavior of the dynamic system. This happened in the case of choosing higher values of damping factor or when exceeding the optimal resonance conditions. Likewise, this occurred in the case of arbitrarily choosing the initial conditions.

Remark 4.

The nonlinear operator $N^{(k)} [v_{0} (t)]$ from Equation (74) has the form

(78) $\begin{matrix} N^{(k)} [v_{0} (t)] = \sum_{i = 1}^{n_{s}} h_{i} (t) g_{i} (t), k = 1, 2, \dots, \end{matrix}$

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 (75) as

(79) $\begin{matrix} L [v_{1} (t, C_{i})] + A_{1} [v_{0} (t), C_{i}] h_{1} (t) + A_{2} [v_{0} (t), C_{j}] h_{2} (t) + A_{3} [v_{0} (t), C_{k}] h_{3} (t) + \dots = 0, \end{matrix}$

with $A_{1} [v_{0} (t), C_{i}]$ , $A_{2} [v_{0} (t), C_{j}]$ , $A_{3} [v_{0} (t), C_{k}]$ , … arbitrary auxiliary functions depending on the unknown parameters $C_{i}$ , $C_{j}$ , $C_{k}$ , … being optimally computed via various methods, as follows: the collocation method, the least squares method, the Galerkin method, the weighted residual method, and so on.

Using the linearly independent functions $h_{1}, h_{2}, \dots, h_{n_{s}}$ , we introduce some types of approximate solutions of Equation (70).

Definition 1.

A sequence of functions ${s_{m} (t)}_{m \geq 1}$ of the form

(80) $s_{m} (t) = \sum_{i = 1}^{m} α_{m}^{i} \cdot h_{i} (t), m \geq 1, α_{m}^{i} \in R,$

is called an OAFM sequence of Equation (70).

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

The OAFM sequences ${s_{m} (t)}_{m \geq 1}$ with the property

$lim_{m \to \infty} R (t, s_{m} (t)) = 0$

is called convergent to the solution of Equation (70), where $R (t, u (t)) = L [u (t)] + N [u (t)] + g (t)$ .

Definition 2.

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

(81) $| R (t, \tilde{F} (t)) | < ε, B (\tilde{F} (t, C_{i}), \frac{d \tilde{F} (t, C_{i})}{d t}) = 0$

are called ε-approximate OAFM solutions of Equation (70).

Definition 3.

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

(82) $\int_{0}^{\infty} R^{2} (t, \tilde{F} (t)) d t \leq ε, B (\tilde{F} (t, C_{i}), \frac{d \tilde{F} (t, C_{i})}{d t}) = 0$

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

Remark 5.

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

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

Theorem 1.

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

Proof.

It is similar to the Theorem from [26].

Firstly, the OAFM sequences ${s_{m}}_{m \geq 1}$ are built by considering the approximate OAFM solutions of the following type:

(83) $\bar{F} (t) = \sum_{i = 1}^{n} C_{m}^{i} \cdot h_{i} (t), where m \geq 1 is fixed arbitrary .$

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

Introducing the approximate solutions $\bar{F}$ in Equation (70), the expression yields the following:

$R (t, C_{m}^{i}) : = R (t, \bar{F}) = L [\bar{F} (t)] + N [\bar{F} (t)] + g (t) .$

Attaching to Equation (70) the following real functional

(84) $J_{1} (C_{m}^{i}) = \int_{0}^{\infty} R^{2} (t, C_{m}^{i}) d t$

and imposing the initial conditions, we can determine

l \in N

l \leq 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 (84), which give the minimum of functional (84).

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

(85) $s_{m} (t) = \sum_{i = 1}^{n} {\tilde{C}}_{m}^{i} \cdot h_{i} (t)$

are constructed.

The next step is to show that the above OAFM functions $s_{m} (t)$ are weak $ε$ -approximate OAFM solutions of Equation (70).

Therefore, the OAFM functions $s_{m} (t)$ are computed, and taking into account that $\bar{F}$ given by (83) are OAFM functions for Equation (70), it follows that

$0 \leq \int_{0}^{\infty} R^{2} (t, s_{m} (t)) d t \leq \int_{0}^{\infty} R^{2} (t, \bar{F} (t)) d t, \forall m \geq 1 .$

Thus,

$0 \leq lim_{m \to \infty} \int_{0}^{\infty} R^{2} (t, s_{m} (t)) d t \leq lim_{m \to \infty} \int_{0}^{\infty} R^{2} (t, \bar{F} (t)) d t .$

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

$lim_{m \to \infty} \int_{0}^{\infty} R^{2} (t, s_{m} (t)) d t = 0 .$

It follows that for all $ε > 0$ , there exists $m_{0} \geq 1$ such that for all $m \geq 1$ and $m > m_{0}$ , the sequence $s_{m} (t)$ is a weak $ε$ -approximate OAFM solution of Equation (70). □

Remark 6.

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

4.2. Semi-Analytical Solutions via the OAFM Procedure

The nonlinear differential problems given by Equations (61) and (63) have the following form:

(86) $\{\begin{matrix} \ddot{u} + α u^{3} + β u = 0 \\ u (0) = γ, \dot{u} (0) = δ \end{matrix},$

where the coefficients

α

β

γ

, and

δ

could be identified in Equations (61) and (63).

Thus, the steps of the applicability of the OAFM technique presented in Section 3 become as follows:

(i). The linear and the corresponding nonlinear operators are as follows:

(87) $\begin{matrix} L [u (t)] = \ddot{u} + ω_{0}^{2} u, \\ N [u (t)] = α u^{3} + A_{0} u, \end{matrix}$

where

A_{0} = β - ω_{0}^{2}

;

(ii). Write the OAFM-solution $\bar{u} (t)$ as a sum of the initial approximate solution $u_{0} (t)$ and the first approximation $u_{1} (t)$ (see Equation (71));

(iii). Equation (72) (with $g (t) = 0$ ) is equivalent to

(88) $\begin{matrix} L [u_{0} (t)] = 0, \\ u_{0} (0) = γ, {\dot{u}}_{0} (0)) = δ, \end{matrix}$

whose solution is

(89) $\begin{matrix} u_{0} (t) = A cos (ω_{0} t) + B sin (ω_{0} t), \end{matrix}$

where

A = γ

and

B = \frac{δ}{ω_{0}}

;

(iv). The expression $N [u_{0} (t)]$ that appears in Equation (78) is as follows:

(90) $\begin{matrix} N [u_{0} (t)] = α u_{0}^{3} + A_{0} u_{0} = g_{1} (t) h_{1} (t) + g_{2} (t) h_{2} (t) + g_{3} (t) h_{3} (t) + g_{4} (t) h_{4} (t), \end{matrix}$

where

$\begin{matrix} g_{1} (t) = A A_{0} + \frac{1}{2} α A^{3} + \frac{1}{2} α A^{3} cos (2 ω_{0} t), & h_{1} (t) = cos (ω_{0} t), \\ g_{2} (t) = \frac{3}{2} α A B^{2} sin (2 ω_{0} t), & h_{2} (t) = sin (ω_{0} t), \\ g_{3} (t) = \frac{3}{2} α A^{2} B sin (2 ω_{0} t), & h_{3} (t) = cos (ω_{0} t), \\ g_{4} (t) = B A_{0} + \frac{1}{2} α B^{3} - \frac{1}{2} α B^{3} cos (2 ω_{0} t), & h_{4} (t) = sin (ω_{0} t); \end{matrix}$

(v). For the computation of the first approximation $u_{1} (t)$ (from Equation (79)), it is more convenient to chose the auxiliary functions defined as follows:

$\begin{matrix} A_{1} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{B}}_{j} cos (2 j ω_{0} t), & A_{2} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{B}}_{j} sin (2 j ω_{0} t), \\ A_{3} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{C}}_{j} sin (2 j ω_{0} t), & A_{4} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{C}}_{j} cos (2 j ω_{0} t); \end{matrix}$

(vi). Combining Equations (79) and (87) yields the following equation:

(91) $\begin{matrix} {\ddot{u}}_{1} + ω_{0}^{2} u_{1} = \sum_{j = 1}^{N_{m a x}} {\tilde{B}}_{j} cos ((2 j + 1) ω_{0} t) + {\tilde{C}}_{j} sin ((2 j + 1) ω_{0} t), \end{matrix}$

whose solution has the form

(92) $\begin{matrix} u_{1} = \sum_{j = 1}^{N_{m a x}} B_{j} cos ((2 j + 1) ω_{0} t) + C_{j} sin ((2 j + 1) ω_{0} t), \end{matrix}$

where the unknown parameters

B_{j}

and

C_{j}

depend on the

{\tilde{B}}_{j}

and

{\tilde{C}}_{j}

and will be optimally identified;

(vii). Finally, the OAFM-solution $\bar{u}$ , obtained via Equation (71), is well determined by

(93) $\begin{matrix} \bar{u} (t) = u_{0} (t) + u_{1} (t), \end{matrix}$

with

u_{0} (t)

and

u_{1} (t)

given by Equations (89) and (92), respectively.

The nonlinear differential problems given by Equations (65), (67) and (69) have the form:

(94) $\{\begin{matrix} \ddot{u} (1 + u^{2}) - 2 u {(\dot{u})}^{2} + α_{1} u (1 + u^{2}) + β_{1} u (1 - u^{2}) = 0 \\ u (0) = γ_{1}, \dot{u} (0) = δ_{1} \end{matrix},$

where the coefficients

α_{1}

β_{1}

γ_{1}

, and

δ_{1}

depend of the physical parameters a, b, c, d, and k from Equations (65), (67), and (69).

Thus, the steps of the applicability of the OAFM technique presented in Section 3 become as follows:

(i). The linear and the corresponding nonlinear operators are as follows:

(95) $\begin{matrix} L [u (t)] = \ddot{u} + ω_{0}^{2} u \\ N [u (t)] = {\ddot{u}}_{0} u_{0}^{2} - 2 u_{0} {({\dot{u}}_{0})}^{2} + α_{1} u_{0} (1 + u_{0}^{2}) + β_{1} u_{0} (1 - u_{0}^{2}) - ω_{0}^{2} u_{0}; \end{matrix}$

(ii). Write the OAFM-solution $\bar{u} (t)$ as a sum of the initial approximate solution $u_{0} (t)$ and the first approximation $u_{1} (t)$ (see Equation (71));

(iii). Equation (72) (with $g (t) = 0$ ) is equivalent to

(96) $\begin{matrix} L [u_{0} (t)] = 0, \\ u_{0} (0) = γ_{1}, {\dot{u}}_{0} (0)) = δ_{1}, \end{matrix}$

whose solution is

(97) $\begin{matrix} u_{0} (t) = A cos (ω_{0} t) + B sin (ω_{0} t), \end{matrix}$

where

A = γ_{1}

and

B = \frac{δ_{1}}{ω_{0}}

;

(iv). The expression $N [u_{0} (t)]$ that appears in Equation (78) is as follows:

(98) $\begin{matrix} N [u_{0} (t)] = {\ddot{u}}_{0} u_{0}^{2} - 2 u_{0} {({\dot{u}}_{0})}^{2} + α_{1} u_{0} (1 + u_{0}^{2}) + β_{1} u_{0} (1 - u_{0}^{2}) - ω_{0}^{2} u_{0} = \\ = g_{1} (t) h_{1} (t) + g_{2} (t) h_{2} (t) + g_{3} (t) h_{3} (t) + g_{4} (t) h_{4} (t), \end{matrix}$

where

$\begin{matrix} g_{1} (t) = A (α_{1} + β_{1} - ω_{0}^{2}) + [A^{3} (α_{1} - β_{1} - ω_{0}^{2}) - 2 A B^{2} ω_{0}^{2}] (\frac{1}{2} - \frac{1}{2} cos (2 ω_{0} t)), \\ g_{2} (t) = \frac{1}{2} (3 A B^{2} α_{1} - 3 A B^{2} β_{1} - 2 A^{3} ω_{0}^{2} + A B^{2} ω_{0}^{2}) sin (2 ω_{0} t), \\ g_{3} (t) = \frac{1}{2} (3 A^{2} B α_{1} - 3 A^{2} B β_{1} - 2 B^{3} ω_{0}^{2} + A^{2} B ω_{0}^{2}) sin (2 ω_{0} t), \\ g_{4} (t) = B (α_{1} + β_{1} - ω_{0}^{2}) + [B^{3} (α_{1} - β_{1} - ω_{0}^{2}) - 2 A^{2} B ω_{0}^{2}] (\frac{1}{2} - \frac{1}{2} cos (2 ω_{0} t)), \\ h_{1} (t) = cos (ω_{0} t), h_{2} (t) = sin (ω_{0} t), h_{3} (t) = cos (ω_{0} t), h_{4} (t) = sin (ω_{0} t); \end{matrix}$

(v). For the computation of the first approximation $u_{1} (t)$ (from Equation (79)), it is more convenient to chose the auxiliary functions defined by the following:

$\begin{matrix} A_{1} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{B}}_{j} cos (M j ω_{0} t), & A_{2} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{B}}_{j} sin (M j ω_{0} t), \\ A_{3} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{C}}_{j} sin (M j ω_{0} t), & A_{4} (t) = \sum_{j = 1}^{N_{m a x}} {\tilde{C}}_{j} cos (M j ω_{0} t), \end{matrix}$

where

M \in {2, 3, 4}

is an integer number used for the convergence control of the OAFM-solution;

(vi). Combining Equations (79) and (95) yields the following equation:

(99) $\begin{matrix} {\ddot{u}}_{1} + ω_{0}^{2} u_{1} = \sum_{j = 1}^{N_{m a x}} {\tilde{B}}_{j} cos ((M j + 1) ω_{0} t) + {\tilde{C}}_{j} sin ((M j + 1) ω_{0} t), \end{matrix}$

whose solution has the form

(100) $\begin{matrix} u_{1} = \sum_{j = 1}^{N_{m a x}} B_{j} cos ((M j + 1) ω_{0} t) + C_{j} sin ((M j + 1) ω_{0} t), \end{matrix}$

where the unknown parameters

B_{j}

and

C_{j}

depend on the

{\tilde{B}}_{j}

and

{\tilde{C}}_{j}

and will be optimally identified;

(vii). Finally, the OAFM-solution $\bar{u}$ , obtained via Equation (71), is well determined by

(101) $\begin{matrix} \bar{u} (t) = u_{0} (t) + u_{1} (t), \end{matrix}$

with

u_{0} (t)

and

u_{1} (t)

given by Equations (97) and (100), respectively.

5. Numerical Results and Validation

The existence of periodic orbits ensures the periodic behavior of OAFM-solutions and the corresponding numerical solutions obtained by the fourth-order Runge–Kutta method.

In this section, the OAFM solutions are validated for different numerical values of the physical parameters of the system (1). The comparison between OAFM solutions ${\bar{u}}_{O A F M}$ and the corresponding numerical solutions for the given initial conditions, some physical constants a, b, c, d, $k > 0$ , and index $N_{m a x}$ is presented in details in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11, Table 12, Table 13, Table 14, Table 15, Table 16, Table 17, Table 18, Table 19 and Table 20. Moreover, the absolute values $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ) are highlighted in these Tables to argue for the accuracy of the obtained results.

The performance of the OAFM-solutions ${\bar{u}}_{O A F M}$ and the states $(x (t), y (t), z (t))$ of the studied system, by comparison with the corresponding numerical ones, is qualitatively reflected by Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17 for different physical parameter values. Figure 18 shows periodic orbit $(x (t), y (t), z (t))$ for $a = 0.5$ , $b = 0.85$ , $c = 0.75$ , $d = - 0.45$ , and $k = 0.65$ .

Unlike the homotopy perturbation methods (requiring a large number of iterations), a first advantage of the OAFM technique is the convergence control with the index $N_{m a x}$ .

The convergence control is regarded as a qualitative property of the OAFM-solutions as close as possible to the exact solution implicitly defined by Equations (3) and (41).

Testing is done using a rigorous numerical method such as the fourth-order Runge–Kutta method. The absolute values denoted by $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ are presented in the following tables.

There are 3D dynamical systems that admit only a prime integral, so the exact solution is not known, but with the help of the prime integral and the OAFM method, one can construct semi-analytic solutions in closed form.

In the case of complex dynamical systems, even with several degrees of freedom, the convergence control is more difficult and new numerical/analytical methods are sought.

Comparison with the Iterative Method

In 2006, Daftardar-Gejji et al. [26] developed an iterative method for solving nonlinear functional equations. This method is validated for a class of nonlinear functional equations (such as the Voltera-type integral equation, the nonlinear 3D- differential system, and the fractional-order differential equation) whose exact solutions are known in analytical form.

The advantage of the OAFM method compared with the iterative method described in [27] (using eight iterations) is highlighted in this subsection. In Table 21 and Figure 19, respectively, the precision and efficiency of the OAFM method can be seen.

System (1) is integrated over the interval $[0, t]$ , and it obtains the following:

(102) $\begin{matrix} x (t) = x (0) + \int_{0}^{t} (a y (s) + b y (s) z (s)) d s \\ y (t) = y (0) + \int_{0}^{t} (c x (s) z (s) + d x (s)) d s \\ z (t) = z (0) + \int_{0}^{t} (- k x (s) y (s)) d s \end{matrix} .$

The iterative algorithm leads to the following:

(103) $\begin{matrix} x_{0} (t) = x (0), x_{1} (t) = N_{1} (x_{0}, y_{0}, z_{0}) = \int_{0}^{t} (a y (0) + b y (0) z (0)) d s, \\ y_{0} (t) = y (0), y_{1} (t) = N_{2} (x_{0}, y_{0}, z_{0}) = \int_{0}^{t} (c x (0) z (0) + d x (0)) d s, \\ z_{0} (t) = z (0), z_{1} (t) = N_{3} (x_{0}, y_{0}, z_{0}) = \int_{0}^{t} (- k x (0) y (0)) d s, \\ \dots \\ x_{m} (t) = N_{1} (\sum_{i = 0}^{m - 1} x_{i}, \sum_{i = 0}^{m - 1} y_{i}, \sum_{i = 0}^{m - 1} z_{i}) - N_{1} (\sum_{i = 0}^{m - 2} x_{i}, \sum_{i = 0}^{m - 2} y_{i}, \sum_{i = 0}^{m - 2} z_{i}), \\ y_{m} (t) = N_{2} (\sum_{i = 0}^{m - 1} x_{i}, \sum_{i = 0}^{m - 1} y_{i}, \sum_{i = 0}^{m - 1} z_{i}) - N_{2} (\sum_{i = 0}^{m - 2} x_{i}, \sum_{i = 0}^{m - 2} y_{i}, \sum_{i = 0}^{m - 2} z_{i}), \\ z_{m} (t) = N_{3} (\sum_{i = 0}^{m - 1} x_{i}, \sum_{i = 0}^{m - 1} y_{i}, \sum_{i = 0}^{m - 1} z_{i}) - N_{3} (\sum_{i = 0}^{m - 2} x_{i}, \sum_{i = 0}^{m - 2} y_{i}, \sum_{i = 0}^{m - 2} z_{i}), \\ m \geq 2 . \end{matrix}$

The iterative algorithm leads to the solution of the Equation (1) as follows:

$x_{i t e r} (t) = \sum_{m = 0}^{\infty} x_{m} (t), y_{i t e r} (t) = \sum_{m = 0}^{\infty} y_{m} (t), z_{i t e r} (t) = \sum_{m = 0}^{\infty} z_{m} (t),$

For the following initial conditions: $x (0) = 0.25$ , $y (0) = 0.5$ , $z (0) = 1.5$ and physical constants $a = 0.5 > 0$ , $b = 0$ , $c = 0.75$ , $d = - 0.45 < 0$ , $k = 0.65$ (presented in Table 21), the iterative solutions $x_{i t e r} (t)$ , after eight iterations using the algorithm (103), have the following form:

(104) $\begin{matrix} x_{i t e r} (t) = \sum_{m = 0}^{8} x_{m} (t) = 0.25 + 0.25 t + 0.0421875 t^{2} + 0.0127929687 t^{3} \\ + 0.0001272583 t^{4} - 0.0001868476 t^{5} - 0.0001127917 t^{6} - 0.0000280089 t^{7} \\ - 4.8371 \times 10^{- 6} t^{8} - 5.0463 \times 10^{- 7} t^{9} + 6.9139 \times 10^{- 8} t^{10} + 4.3810 \times 10^{- 8} t^{11} \\ + 1.2752 \times 10^{- 8} t^{12} + 2.4153 \times 10^{- 9} t^{13} + 2.8401 \times 10^{- 10} t^{14} - 4.4590 \times 10^{- 12} t^{15} \\ - 1.2718 \times 10^{- 11} t^{16} - 3.8729 \times 10^{- 12} t^{17} + O (10^{- 13}); \\ y_{i t e r} (t) = \sum_{m = 0}^{8} y_{m} (t) = 0.5 + 0.16875 t + 0.0767578125 t^{2} \\ + 0.0010180664 t^{3} - 0.0018684768 t^{4} - 0.0013535010 t^{5} - 0.0003921249 t^{6} \\ - 0.0000773944 t^{7} - 8.9987 \times 10^{- 6} t^{8} + 1.3144 \times 10^{- 6} t^{9} \\ + 9.5956 \times 10^{- 7} t^{10} + 3.1557 \times 10^{- 7} t^{11} + 6.7129 \times 10^{- 8} t^{12} + 8.9845 \times 10^{- 9} t^{13} \\ - 1.4733 \times 10^{- 10} t^{14} - 5.1401 \times 10^{- 10} t^{15} - 1.8139 \times 10^{- 10} t^{16} - 4.1114 \times 10^{- 11} t^{17} \\ - 6.1228 \times 10^{- 12} t^{18} + O (10^{- 13}); \\ z_{i t e r} (t) = \sum_{m = 0}^{8} z_{m} (t) = 1.5 - 0.081250 t - 0.0543359375 t^{2} \\ - 0.0178686523 t^{3} - 0.0053559341 t^{4} - 0.0006822478 t^{5} - 0.0000159755 t^{6} \\ + 0.0000538912 t^{7} + 0.0000199578 t^{8} + 5.1792 \times 10^{- 6} t^{9} \\ + 8.6714 \times 10^{- 7} t^{10} + 4.8166 \times 10^{- 8} t^{11} - 2.8151 \times 10^{- 8} t^{12} - 1.3599 \times 10^{- 8} t^{13} \\ - 3.6087 \times 10^{- 9} t^{14} - 6.5933 \times 10^{- 10} t^{15} - 6.1959 \times 10^{- 11} t^{16} + 9.3548 \times 10^{- 12} t^{17} \\ + 6.2900 \times 10^{- 12} t^{18} + 1.7947 \times 10^{- 12} t^{19} + O (10^{- 13}) . \end{matrix}$

Figure 19 shows the the OAFM solutions ${\bar{x}}_{O A F M}$ , ${\bar{y}}_{O A F M}$ , ${\bar{z}}_{O A F M}$ , compared with the corresponding iterative solutions given in Equation (104) and the corresponding numerical ones. This comparison reflects the accuracy of the OAFM approach using only one iteration.

6. Conclusions

The present paper is dedicated to developing the OAFM analytical approach for a class of dynamical system. This dynamical system admits two prime integrals (constants of motion), involving the existence of symmetry with respect to the $O z$ -axis and the existence of periodic, homoclinic, or heteroclinic orbits. Good agreement between the OAFM results and the corresponding numerical ones is demonstrated in this work. The comparative analysis reflects the accuracy of the method, as the obtained solutions are close to the exact solution. The advantages of the OAFM method compared with the iterative method are synthesized by the following:

The arbitrary choice of the linear operator L and auxiliary functions $A_{j}$ allows us to write the OAFM-solution in an effective form;

The convergence control (in the sense that the residual functions are smaller than 1): The auxiliary functions depend on the finite number of the unknown parameters which could be optimally computed to obtain the absolute values between the semi-analytical solutions and numerical ones smaller than 1;

The non-existence of small parameters: The OAFM-solutions are built for arbitrary values for the physical parameters a, b, c, d, k, namely, for $a c - b d \neq 0$ in all subcases $b c = 0$ , $b c > 0$ , and $b c < 0$ , respectively; A special case when $a c - b d = 0$ is analytically investigated. The OAFM-solution is obtained in the same manner.

The efficiency is highlighted by the comparison between the OAFM-solutions with the corresponding iterative solutions, more specifically, the iterative method is validated only when the exact solution is known; on the other hand, by the OAFM method, the exact solution is approximated by arbitrarily choosing the auxiliary functions and optimally finding the unknown parameters, proving the accuracy of the OAFM-results.

The achievement of these results encourages the study of dynamical systems with similar properties as models that can describe physical and biological problems as follows: the change in the energy levels of an atom with two levels as a function of the eigenvalue; or in biology, the bi-Hamiltonian structure of Kermack and McKendrick modeling for the spread of epidemics in a closed population.

Author Contributions

Conceptualization, N.P.; data curation, R.-D.E. and 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

Dataset available on request from the authors.

Conflicts of Interest

The authors declare no conflict 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 Equilibrium point $e_{2}^{M} (0, 0, M)$ (black point) and the homoclinic orbits $(x (t), y (t), z (t))$ (magenta), $(x (t), - y (t), z (t))$ (red), $(- x (t), y (t), z (t))$ (green), and $(x (t), y (t), z (t))$ (blue) using Equations (11) and (12) for $M = 5.15$ , $a = 0.5$ , $b = 0$ , $c = 0.75$ , $d = - 0.45$ , $k = 0.65$ ; initial data: $x_{0} = 0.01$ , $y_{0} = 0.002$ , $z_{0} = M$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 2 Equilibrium points $e_{1}^{M} (M, 0, - d / c)$ (black point) and $e_{1}^{- M} (- M, 0, - d / c)$ (green point); the heteroclinic orbits $(x_{1} (t), y_{1} (t), z_{1} (t))$ (red), $(- x_{1} (t), y_{1} (t), z_{1} (t))$ (green), $(x_{1} (t), - y_{1} (t), z_{1} (t))$ (brown), $(- x_{1} (t), - y_{1} (t), z_{1} (t))$ (magenta); and the split-heteroclinic orbits $(x_{2} (t), y_{2} (t), z_{2} (t))$ , $(- x_{2} (t), y_{2} (t), z_{2} (t))$ , $(x_{2} (t), - y_{2} (t), z_{2} (t))$ , $(- x_{2} (t), - y_{2} (t), z_{2} (t))$ (blue), using Equations (16) and (17) for $M = 3.25$ , $a = 0.5$ , $b = 0$ , $c = - 0.75$ , $d = - 0.45$ , $k = 0.65$ ; initial data: $x_{0} = 3.24$ , $y_{0} = 0.001$ , $z_{0} = - 0.59$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 3 Equilibrium point $e_{4}^{M} (0, 0, M)$ (black point) and the homoclinic orbits $(x (t), y (t), z (t))$ (magenta), $(- x (t), y (t), z (t))$ (red), $(x (t), - y (t), z (t))$ (green), and $(x (t), y (t), z (t))$ (blue), using Equations (24) and (25) for $M = - 1.15$ , $a = 0.5$ , $b = 0.75$ , $c = 0$ , $d = - 0.45$ , $k = 0.65$ ; initial data: $x_{0} = 0.01$ , $y_{0} = 0.02$ , $z_{0} = M + 0.01$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 4 Equilibrium points $e_{3}^{M} (0, M, - a / b)$ (black point) and $e_{3}^{- M} (0, - M, - d / c)$ (green point); heteroclinic orbits $(x_{1} (t), y_{1} (t), z_{1} (t))$ (red), $(- x_{1} (t), y_{1} (t), z_{1} (t))$ (green), $(x_{1} (t), - y_{1} (t), z_{1} (t))$ (brown), $(- x_{1} (t), - y_{1} (t), z_{1} (t))$ (magenta); and split-heteroclinic orbits $(x_{2} (t), y_{2} (t), z_{2} (t))$ , $(- x_{2} (t), y_{2} (t), z_{2} (t))$ , $(x_{2} (t), - y_{2} (t), z_{2} (t))$ , $(- x_{2} (t), - y_{2} (t), z_{2} (t))$ (blue), using Equations (29) and (30) for $M = 2.25$ , $a = 0.5$ , $b = - 0.75$ , $c = 0$ , $d = - 0.45$ , $k = 0.65$ ; initial data: $x_{0} = 0.001$ , $y_{0} = M + 0.001$ , $z_{0} = - a / b - 0.001$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 5 Equilibrium point $e_{7}^{M} (0, 0, M)$ (black point) and the homoclinic orbits $(x_{i} (t), y_{i} (t), z_{i} (t))$ , $i = 1, 2$ , using Equations (38) and (39) for $a = 0.5$ , $b = 0.85$ , $c = 0.75$ , $d = - 0.45$ , $k = 0.65$ ; initial data: $x_{0} = - 0.001$ , $y_{0} = - 0.002$ , $z_{0} = M + 0.001$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively: (a) $M = - 0.65$ ; (b) $M = - 2.65$ .

Figure 6 Equilibrium point $e_{7}^{M} (0, 0, M)$ (black point) and the homoclinic orbits $(x_{i} (t), y_{i} (t), z_{i} (t))$ , $i = 1, 2$ , using Equations (38) and (39) for $a = 0.5$ , $b = 0.85$ , $c = 0.75$ , $d = - 0.45$ , $k = 0.65$ ; initial data: $x_{0} = - 0.001$ , $y_{0} = - 0.002$ , $z_{0} = M + 0.001$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively. (a) $M = - 0.65$ ; (b) $M = - 2.65$ .

Figure 7 Equilibrium points $e_{8}^{M} (0, M, - a / b)$ (black point) and $e_{8}^{- M} (0, - M, - a / b)$ (orange point) and the heteroclinic orbits $(x_{1} (t), y_{1} (t), z_{1} (t))$ (red), $(- x_{1} (t), y_{1} (t), z_{1} (t))$ (green), $(x_{2} (t), y_{2} (t), z_{2} (t))$ (blue), $(- x_{2} (t), y_{2} (t), z_{2} (t))$ (orange) using Equations (50) and (51) for $M = 2.05$ , $a = 0.5$ , $b = - 0.25$ , $c = 0.75$ , $d = \frac{a c}{b}$ , $k = 0.65$ ; initial data: $x_{0} = 0.002$ , $y_{0} = M + 0.001$ , $z_{0} = - a / b + 0.001$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 8 Equilibrium points $e_{9}^{M} (M, 0, - a / b)$ (black point) and $e_{9}^{- M} (- M, 0, - a / b)$ (orange point) and the heteroclinic orbits $(x_{1} (t), y_{1} (t), z_{1} (t))$ (red), $(x_{1} (t), - y_{1} (t), z_{1} (t))$ (green), $(x_{2} (t), y_{2} (t), z_{2} (t))$ (blue), $(x_{2} (t), - y_{2} (t), z_{2} (t))$ (orange), using Equations (54) and (55) for $M = 2.05$ , $a = 0.5$ , $b = 0.25$ , $c = - 0.75$ , $d = \frac{a c}{b}$ , $k = 0.65$ ; initial data: $x_{0} = M + 0.001$ , $y_{0} = 0.002$ , $z_{0} = - a / b + 0.001$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 9 Equilibrium points $e_{10}^{M} (0, 0, M)$ (black point) and $e_{10}^{- M - 2 a / b} (0, 0, - M - \frac{2 a}{b})$ (green point) and the heteroclinic orbits $(x_{1} (t), y_{1} (t), z_{1} (t))$ (red), $(- x_{1} (t), y_{1} (t), z_{1} (t))$ (green), $(x_{2} (t), y_{2} (t), z_{2} (t))$ (blue), $(- x_{2} (t), y_{2} (t), z_{2} (t))$ (orange), using Equations (58) and (59) for $M = 2.05$ , $a = 0.5$ , $b = 0.25$ , $c = 0.75$ , $d = \frac{a c}{b}$ , $k = 0.65$ ; initial data: $x_{0} = 0.001$ , $y_{0} = 0.002$ , $z_{0} = M + 0.001$ —exact solution (dotted color curve) and numerical results (solid black curve), respectively.

Figure 10 Performance of the OAFM-solution ${\bar{u}}_{O A F M} (t)$ of Equation (61) (dotted black curve), compared with the corresponding numerical solutions (solid color curve), respectively.

Figure 11 Performance of the OAFM-solution ${\bar{u}}_{O A F M} (t)$ of Equation (63) (dotted black curve), compared with the corresponding numerical ones (solid color curve), respectively.

Figure 12 States of the system $({\bar{x}}_{O A F M}, {\bar{y}}_{O A F M}, {\bar{z}}_{O A F M})$ using Equations (101)–(A3) for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.85 > 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ compared with numerical integration: OAFM solutions (dashing black curve) and numerical results (solid color curve: green line for $x (t)$ , red curve for $y (t)$ , blue curve for $z (t)$ ), respectively.

Figure 13 Performance of the OAFM-solution ${\bar{u}}_{O A F M} (t)$ of Equation (65) (for $(a > 0, d > 0)$ , $(a < 0, d > 0)$ , $(a < 0, d < 0)$ ) and (67) (for $(a > 0, d < 0)$ ), respectively, (dotted black curve), compared with the corresponding numerical ones (solid color curve), respectively.

Figure 14 States of the system $({\bar{x}}_{O A F M}, {\bar{y}}_{O A F M}, {\bar{z}}_{O A F M})$ using Equations (101) and (A4) for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = - 0.85 < 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ compared with numerical integration: OAFM solutions (dashing black curve) and numerical results (solid color curve: green line for $x (t)$ , red curve for $y (t)$ , blue curve for $z (t)$ ), respectively.

Figure 15 Performance of the OAFM-solution ${\bar{u}}_{O A F M} (t)$ of Equation (65) (dotted black curve), compared with the corresponding numerical ones (solid color curve), respectively.

Figure 16 States of the system $({\bar{x}}_{O A F M}, {\bar{y}}_{O A F M}, {\bar{z}}_{O A F M})$ using Equations (101) and (A5) for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.85 > 0$ , $c = - 0.75 < 0$ , $d = - 0.45 < 0$ , $k = 0.65$ compared with numerical integration: OAFM solutions (dashing black curve) and numerical results (solid color curve: green line for $x (t)$ , red curve for $y (t)$ , blue curve for $z (t)$ ), respectively.

Figure 17 Performance of the OAFM-solution ${\bar{u}}_{O A F M} (t)$ of Equation (67) (dotted black curve), compared with the corresponding numerical ones (solid color curve), respectively.

Figure 18 Periodic orbit $(x (t), y (t), z (t))$ of Equation (67) for $a = 0.5$ , $b = 0.85$ , $c = 0.75$ , $d = - 0.45$ , $k = 0.65$ : OAFM-solution (dotted color curve), compared with the corresponding numerical ones (solid black curve), respectively.

Figure 19 Performance of the approximate solutions ${\bar{x}}_{O A F M} (t)$ , ${\bar{y}}_{O A F M} (t)$ , ${\bar{z}}_{O A F M} (t)$ of Equation (1) using Equation (A1) (dashing black curves), compared with the iterative solutions $x_{i t e r} (t)$ , $y_{i t e r} (t)$ , $z_{i t e r} (t)$ using Equation (104) (dotted red curves) and the corresponding numerical ones (solid color curves: green for $x (t)$ , magenta for $y (t)$ , blue for $z (t)$ ), respectively.

Table 1

Function ${\bar{u}}_{O A F M}$ of Equation (61) using Equation (93) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0$ , $c = 0.75$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 35$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.25	0.2499999332	6.679 $\times 10^{- 8}$
2	1.4558317036	1.4558314691	2.345 $\times 10^{- 7}$
4	2.1345855716	2.1345849344	6.372 $\times 10^{- 7}$
6	0.4436803612	0.4436798580	5.031 $\times 10^{- 7}$
8	0.0232122934	0.0232129759	6.824 $\times 10^{- 7}$
10	−0.3023772974	−0.3023767301	5.672 $\times 10^{- 7}$
12	−1.6718036804	−1.6718025094	1.171 $\times 10^{- 6}$
14	−1.9184774993	−1.9184801693	2.669 $\times 10^{- 6}$
16	−0.3706742998	−0.3706745418	2.420 $\times 10^{- 7}$
18	−0.0011567148	−0.0011568841	1.692 $\times 10^{- 7}$
20	0.3636280521	0.3636242064	3.845 $\times 10^{- 6}$

Table 2

Function ${\bar{u}}_{O A F M}$ of Equation (61) using Equation (93) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0$ , $c = 0.75$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 20$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.25	0.2500000099	9.999 $\times 10^{- 9}$
2	−0.3261410070	−0.3261415988	5.917 $\times 10^{- 7}$
4	−0.1103406696	−0.1103398391	8.305 $\times 10^{- 7}$
6	0.3740448500	0.3740444942	3.558 $\times 10^{- 7}$
8	−0.0490568640	−0.0490563370	5.269 $\times 10^{- 7}$
10	−0.3527164604	−0.3527165184	5.806 $\times 10^{- 8}$
12	0.1997022432	0.1997021360	1.072 $\times 10^{- 7}$
14	0.2663682565	0.2663686171	3.606 $\times 10^{- 7}$
16	−0.3143542225	−0.3143554040	1.181 $\times 10^{- 6}$
18	−0.1316087742	−0.1316083111	4.631 $\times 10^{- 7}$
20	0.3714483361	0.3714484021	6.598 $\times 10^{- 8}$

Table 3

Function ${\bar{u}}_{O A F M}$ of Equation (61) using Equations (93) and (A1) from Appendix A and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0$ , $c = 0.75$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.25	0.2500000099	9.977 $\times 10^{- 9}$
2	1.0050873071	1.0050873871	8.002 $\times 10^{- 8}$
4	1.7246494125	1.7246491591	2.534 $\times 10^{- 7}$
6	0.8378585353	0.8378582351	3.001 $\times 10^{- 7}$
8	0.1708461465	0.1708459727	1.737 $\times 10^{- 7}$
10	−0.2603347198	−0.2603347740	5.420 $\times 10^{- 8}$
12	−1.0266201727	−1.0266192992	8.735 $\times 10^{- 7}$
14	−1.7196398283	−1.7196399403	1.120 $\times 10^{- 7}$
16	−0.8183650472	−0.8183654069	3.596 $\times 10^{- 7}$
18	−0.1616133322	−0.1616137781	4.459 $\times 10^{- 7}$
20	0.2708170212	0.2708174195	3.982 $\times 10^{- 7}$

Table 4

Function ${\bar{u}}_{O A F M}$ of Equation (61) using Equation (93) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0$ , $c = 0.75$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 20$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.25	0.2500001000	1.000 $\times 10^{- 7}$
2	−0.2947941015	−0.2947941133	1.177 $\times 10^{- 8}$
4	−0.4642237912	−0.4642238389	4.764 $\times 10^{- 8}$
6	−0.0302489423	−0.0302489246	1.768 $\times 10^{- 8}$
8	0.4416713811	0.4416713733	7.815 $\times 10^{- 9}$
10	0.3413034018	0.3413033961	5.685 $\times 10^{- 9}$
12	−0.1957908959	−0.1957909204	2.447 $\times 10^{- 8}$
14	−0.4856459204	−0.4856460742	1.537 $\times 10^{- 7}$
16	−0.1423006483	−0.1423007238	7.547 $\times 10^{- 8}$
18	0.3801682672	0.3801682756	8.330 $\times 10^{- 9}$
20	0.4139574552	0.4139574204	3.478 $\times 10^{- 8}$

Table 5

Function ${\bar{u}}_{O A F M}$ of Equation (63) using Equation (93) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.75$ , $c = 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 20$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.5	0.5000000099	9.994 $\times 10^{- 9}$
2	0.0723411559	0.0723411567	7.387 $\times 10^{- 10}$
4	−0.4221694541	−0.4221694553	1.227 $\times 10^{- 9}$
6	−0.5006240790	−0.5006241583	7.934 $\times 10^{- 8}$
8	−0.0738573747	−0.0738573584	1.632 $\times 10^{- 8}$
10	0.4211685017	0.4211685458	4.411 $\times 10^{- 8}$
12	0.5012439117	0.5012439151	3.436 $\times 10^{- 9}$
14	0.0753730129	0.0753729844	2.855 $\times 10^{- 8}$
16	−0.4201640909	−0.4201642069	1.159 $\times 10^{- 7}$
18	−0.5018592500	−0.5018593009	5.090 $\times 10^{- 8}$
20	−0.0768881441	−0.0768881187	2.536 $\times 10^{- 8}$

Table 6

Function ${\bar{u}}_{O A F M}$ of Equation (63) using Equation (93) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.75$ , $c = 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.5	0.5000000094	9.480 $\times 10^{- 9}$
2	−0.1895452367	−0.1895452940	5.739 $\times 10^{- 8}$
4	−0.4513828282	−0.4513829176	8.931 $\times 10^{- 8}$
6	0.3021181602	0.3021180596	1.006 $\times 10^{- 7}$
8	0.3746120614	0.3746121816	1.202 $\times 10^{- 7}$
10	−0.3965318318	−0.3965320445	2.127 $\times 10^{- 7}$
12	−0.2749544847	−0.2749543939	9.081 $\times 10^{- 8}$
14	0.4665845178	0.4665846038	8.595 $\times 10^{- 8}$
16	0.1588574777	0.1588575485	7.085 $\times 10^{- 8}$
18	−0.5074046144	−0.5074045868	2.754 $\times 10^{- 8}$
20	−0.0334240359	−0.0334238797	1.562 $\times 10^{- 7}$

Table 7

Function ${\bar{u}}_{O A F M}$ of Equation (63) using Equations (93) and (A2) from Appendix A and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.75$ , $c = 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.5	0.5000000099	9.999 $\times 10^{- 9}$
2	1.0217586587	1.0217585727	8.606 $\times 10^{- 8}$
4	1.6122687506	1.6122687847	3.412 $\times 10^{- 8}$
6	0.9771755875	0.9771755460	4.148 $\times 10^{- 8}$
8	0.4881520246	0.4881522304	2.057 $\times 10^{- 7}$
10	0.5651825033	0.5651830287	5.253 $\times 10^{- 7}$
12	1.2073713225	1.2073719005	5.780 $\times 10^{- 7}$
14	1.5544874544	1.5544870665	3.878 $\times 10^{- 7}$
16	0.8151592588	0.8151585992	6.596 $\times 10^{- 7}$
18	0.4584677770	0.4584678287	5.165 $\times 10^{- 8}$
20	0.6606037837	0.6606045696	7.859 $\times 10^{- 7}$

Table 8

Function ${\bar{u}}_{O A F M}$ of Equation (63) using Equation (93) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.75$ , $c = 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.5	0.5000000099	9.978 $\times 10^{- 9}$
3/2	1.1527585884	1.1527597554	1.167 $\times 10^{- 6}$
3	2.4744494023	2.4744498101	4.078 $\times 10^{- 7}$
9/2	1.4748028496	1.4748031916	3.419 $\times 10^{- 7}$
6	0.5640614977	0.5640617186	2.209 $\times 10^{- 7}$
15/2	0.6087443039	0.6087453576	1.053 $\times 10^{- 6}$
9	1.6353223889	1.6353232912	9.023 $\times 10^{- 7}$
21/2	2.4024547180	2.4024570288	2.310 $\times 10^{- 6}$
12	1.0306913103	1.0306932177	1.907 $\times 10^{- 6}$
27/2	0.4867084231	0.4867092477	8.245 $\times 10^{- 7}$
15	0.8317655324	0.8317642484	1.284 $\times 10^{- 6}$

Table 9

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equations (101) and (A3) from Appendix A and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.85 > 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.2432905977	0.2432906175	1.979 $\times 10^{- 8}$
6/5	0.4846218181	0.4846218036	1.455 $\times 10^{- 8}$
12/5	1.1698729824	1.1698729786	3.855 $\times 10^{- 9}$
18/5	2.7879098198	2.7879098172	2.631 $\times 10^{- 9}$
24/5	5.1712051469	5.1712051326	1.427 $\times 10^{- 8}$
6	4.7825808563	4.7825808486	7.614 $\times 10^{- 9}$
36/5	2.3790273783	2.3790273599	1.835 $\times 10^{- 8}$
42/5	0.9840844666	0.9840844522	1.436 $\times 10^{- 8}$
48/5	0.4139115393	0.4139115307	8.646 $\times 10^{- 9}$
54/5	0.2257634423	0.2257634344	7.917 $\times 10^{- 9}$
12	0.2489790000	0.2489789925	7.488 $\times 10^{- 9}$

Table 10

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.85 > 0$ , $c = 0.75 > 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.1094984914	0.1094985994	1.079 $\times 10^{- 7}$
6/5	0.3568595307	0.3568595545	2.380 $\times 10^{- 8}$
12/5	0.8701398896	0.8701407808	8.912 $\times 10^{- 7}$
18/5	0.5130777811	0.5130785593	7.781 $\times 10^{- 7}$
24/5	0.1529333049	0.1529333257	2.078 $\times 10^{- 8}$
6	0.0767159741	0.0767155401	4.339 $\times 10^{- 7}$
36/5	0.1645666624	0.1645666889	2.646 $\times 10^{- 8}$
42/5	0.5489429386	0.5489443061	1.367 $\times 10^{- 6}$
48/5	0.8520534884	0.8520554740	1.985 $\times 10^{- 6}$
54/5	0.3305690525	0.3305705103	1.457 $\times 10^{- 6}$
12	0.1032733606	0.1032725673	7.932 $\times 10^{- 7}$

Table 11

Function ${\bar{u}}_{O A F M}$ of Equation (67) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.85 > 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.0522026653	0.0522026753	9.996 $\times 10^{- 9}$
6/5	0.3945359265	0.3945360394	1.128 $\times 10^{- 7}$
12/5	0.9175809646	0.9175809224	4.220 $\times 10^{- 8}$
18/5	0.4499470197	0.4499467854	2.343 $\times 10^{- 7}$
24/5	0.0752962496	0.0752958777	3.718 $\times 10^{- 7}$
6	−0.1584114594	−0.1584115081	4.877 $\times 10^{- 8}$
36/5	−0.6461206343	−0.6461211207	4.863 $\times 10^{- 7}$
42/5	−0.8115395235	−0.8115395858	6.225 $\times 10^{- 8}$
48/5	−0.2506525428	−0.2506524188	1.239 $\times 10^{- 7}$
54/5	0.0155266606	0.0155268233	1.627 $\times 10^{- 7}$
12	0.3120474137	0.3120477475	3.337 $\times 10^{- 7}$

Table 12

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.85 > 0$ , $c = 0.75 > 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.1094984914	0.1094985023	1.087 $\times 10^{- 8}$
0.65	0.2528697458	0.2528696148	1.310 $\times 10^{- 7}$
1.3	0.7350513344	0.7350512909	4.349 $\times 10^{- 8}$
1.95	2.1712602031	2.1712602039	7.611 $\times 10^{- 10}$
2.6	5.7155181284	5.7155183086	1.801 $\times 10^{- 7}$
3.25	8.0328720231	8.0328723947	3.715 $\times 10^{- 7}$
3.9	4.1422352619	4.1422355190	2.571 $\times 10^{- 7}$
4.55	1.4654146467	1.4654147873	1.406 $\times 10^{- 7}$
5.2	0.4940690502	0.4940691795	1.292 $\times 10^{- 7}$
5.85	0.1768058969	0.1768061185	2.216 $\times 10^{- 7}$
6.5	0.0969539633	0.0969537347	2.286 $\times 10^{- 7}$

Table 13

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equations (101) and (A4) from Appendix A and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = - 0.85 < 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.2432905977	0.2432906075	9.768 $\times 10^{- 9}$
6/5	0.1236797698	0.1236797694	3.230 $\times 10^{- 10}$
12/5	−0.1951295824	−0.1951295832	8.568 $\times 10^{- 10}$
18/5	−0.1988029060	−0.1988028917	1.424 $\times 10^{- 8}$
24/5	0.1186472859	0.1186472940	8.100 $\times 10^{- 9}$
6	0.2450325438	0.2450325428	1.012 $\times 10^{- 9}$
36/5	−0.0253309420	−0.0253309424	3.551 $\times 10^{- 10}$
42/5	−0.2549696312	−0.2549696381	6.876 $\times 10^{- 9}$
48/5	−0.0715331192	−0.0715331379	1.868 $\times 10^{- 8}$
54/5	0.2269712959	0.2269713040	8.103 $\times 10^{- 9}$
12	0.1583350656	0.1583350670	1.370 $\times 10^{- 9}$

Table 14

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = - 0.85 < 0$ , $c = 0.75 > 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 30$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.1094984914	0.1094985013	9.889 $\times 10^{- 9}$
6/5	−0.0010299734	−0.0010299733	3.567 $\times 10^{- 11}$
12/5	−0.1085945057	−0.1085945059	1.702 $\times 10^{- 10}$
18/5	0.0963577352	0.0963577342	1.047 $\times 10^{- 9}$
24/5	0.0240120471	0.0240120477	5.583 $\times 10^{- 10}$
6	−0.1174324085	−0.1174324087	2.006 $\times 10^{- 10}$
36/5	0.0790773281	0.0790773282	4.065 $\times 10^{- 11}$
42/5	0.0480220973	0.0480220973	5.149 $\times 10^{- 11}$
48/5	−0.1212261140	−0.1212261137	3.006 $\times 10^{- 10}$
54/5	0.0583992243	0.0583992240	3.417 $\times 10^{- 10}$
12	0.0699684587	0.0699684587	4.213 $\times 10^{- 11}$

Table 15

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = - 0.85 < 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.2432905977	0.2432906075	9.817 $\times 10^{- 9}$
6/5	0.2460836480	0.2460836510	2.972 $\times 10^{- 9}$
12/5	0.0594486339	0.0594486337	2.705 $\times 10^{- 10}$
18/5	−0.1731062884	−0.1731062883	8.856 $\times 10^{- 11}$
24/5	−0.2721618354	−0.2721618358	4.013 $\times 10^{- 10}$
6	−0.1618869684	−0.1618869695	1.072 $\times 10^{- 9}$
36/5	0.0732632563	0.0732632565	1.430 $\times 10^{- 10}$
42/5	0.2518290067	0.2518290069	1.173 $\times 10^{- 10}$
48/5	0.2365777280	0.2365777295	1.558 $\times 10^{- 9}$
54/5	0.0391650200	0.0391650206	6.471 $\times 10^{- 10}$
12	−0.1885034750	−0.1885034742	7.429 $\times 10^{- 10}$

Table 16

Function ${\bar{u}}_{O A F M}$ of Equation (65) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = - 0.85 < 0$ , $c = 0.75 > 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 30$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.1094984914	0.1094985013	9.816 $\times 10^{- 9}$
6/5	0.1041306918	0.1041306925	6.098 $\times 10^{- 10}$
12/5	−0.0600726892	−0.0600726887	4.096 $\times 10^{- 10}$
18/5	−0.1323845211	−0.1323845215	4.174 $\times 10^{- 10}$
24/5	−0.0033019941	−0.0033019941	8.815 $\times 10^{- 11}$
6	0.1308385249	0.1308385247	1.977 $\times 10^{- 10}$
36/5	0.0659039160	0.0659039162	2.064 $\times 10^{- 10}$
42/5	−0.0998133832	−0.0998133833	4.908 $\times 10^{- 11}$
48/5	−0.1132275458	−0.1132275460	2.798 $\times 10^{- 10}$
54/5	0.0459343591	0.0459343592	1.119 $\times 10^{- 10}$
12	0.1347910754	0.1347910757	2.867 $\times 10^{- 10}$

Table 17

Function ${\bar{u}}_{O A F M}$ of Equation (67) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.85 > 0$ , $c = - 0.75 < 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.1182123227	0.1182123331	1.034 $\times 10^{- 8}$
6/5	0.1953183019	0.1953182933	8.692 $\times 10^{- 9}$
12/5	−0.0214693501	−0.0214693536	3.485 $\times 10^{- 9}$
18/5	−0.2061591221	−0.2061592258	1.036 $\times 10^{- 7}$
24/5	−0.0804186595	−0.0804187439	8.442 $\times 10^{- 8}$
6	0.1656775930	0.1656775854	7.574 $\times 10^{- 9}$
36/5	0.1629273796	0.1629274019	2.236 $\times 10^{- 8}$
42/5	−0.0845011900	−0.0845011829	7.031 $\times 10^{- 9}$
48/5	−0.2054492568	−0.2054492960	3.916 $\times 10^{- 8}$
54/5	−0.0170509655	−0.0170510898	1.242 $\times 10^{- 7}$
12	0.1968386807	0.1968386429	3.783 $\times 10^{- 8}$

Table 18

Function ${\bar{u}}_{O A F M}$ of Equation (67) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.85 > 0$ , $c = - 0.75 < 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.1182123227	0.1182123327	9.945 $\times 10^{- 9}$
6/5	0.2824553972	0.2824554002	2.911 $\times 10^{- 9}$
12/5	0.2690054664	0.2690054664	3.635 $\times 10^{- 11}$
18/5	0.0861696160	0.0861696151	9.142 $\times 10^{- 10}$
24/5	−0.1513559476	−0.1513559445	3.034 $\times 10^{- 9}$
6	−0.2929996066	−0.2929996051	1.481 $\times 10^{- 9}$
36/5	−0.2504544150	−0.2504544151	1.462 $\times 10^{- 10}$
42/5	−0.0500180308	−0.0500180304	3.630 $\times 10^{- 10}$
48/5	0.1821824458	0.1821824445	1.328 $\times 10^{- 9}$
54/5	0.2991168963	0.2991168948	1.450 $\times 10^{- 9}$
12	0.2281008479	0.2281008477	1.368 $\times 10^{- 10}$

Table 19

Function ${\bar{u}}_{O A F M}$ of Equation (67) using Equations (101) and (A5) from Appendix A and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.85 > 0$ , $c = - 0.75 < 0$ , $d = - 0.45 < 0$ , $k = 0.65$ ; and index $N_{m a x} = 25$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.0522026653	0.0522026751	9.816 $\times 10^{- 9}$
13/10	0.0653671394	0.0653671395	1.382 $\times 10^{- 10}$
13/5	−0.1240143138	−0.1240143139	3.461 $\times 10^{- 11}$
39/10	0.0708838947	0.0708838946	7.741 $\times 10^{- 11}$
26/5	0.0461416083	0.0461416083	6.125 $\times 10^{- 11}$
13/2	−0.1215738696	−0.1215738699	2.653 $\times 10^{- 10}$
39/5	0.0874276503	0.0874276505	1.558 $\times 10^{- 10}$
91/10	0.0255245727	0.0255245727	3.568 $\times 10^{- 11}$
52/5	−0.1154680122	−0.1154680122	6.699 $\times 10^{- 12}$
117/10	0.1013351501	0.1013351500	1.002 $\times 10^{- 10}$
13	0.0041377795	0.0041377813	1.838 $\times 10^{- 9}$

Table 20

Function ${\bar{u}}_{O A F M}$ of Equation (67) using Equation (101) and numerical integration for $x_{0} = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.85 > 0$ , $c = - 0.75 < 0$ , $d = 0.45 > 0$ , $k = 0.65$ ; and index $N_{m a x} = 30$ (absolute values: $ϵ_{u} = | u_{n u m e r i c a l} - {\bar{u}}_{O A F M} |$ ).

t	$u_{numerical}$	${\bar{u}}_{OAFM}$	$ϵ_{u}$
0	0.0522026653	0.0522026751	9.852 $\times 10^{- 9}$
6/5	0.1862045108	0.1862045110	2.311 $\times 10^{- 10}$
12/5	0.0719839573	0.0719839575	1.071 $\times 10^{- 10}$
18/5	−0.1401037547	−0.1401037500	4.706 $\times 10^{- 9}$
24/5	−0.1634923955	−0.1634923801	1.534 $\times 10^{- 8}$
6	0.0322666420	0.0322666438	1.732 $\times 10^{- 9}$
36/5	0.1840423135	0.1840423143	8.054 $\times 10^{- 10}$
42/5	0.0903401684	0.0903401539	1.457 $\times 10^{- 8}$
48/5	−0.1259524257	−0.1259524557	3.003 $\times 10^{- 8}$
54/5	−0.1721822006	−0.1721822089	8.218 $\times 10^{- 9}$
12	0.0119219983	0.0119219978	5.307 $\times 10^{- 10}$

Table 21

Values of the OAFM-solution ${\bar{x}}_{O A F M} (t)$ (A1), the iterative solution $x_{i t e r} (t)$ , and the numerical integration.

t	$x_{numerical}$	${\bar{x}}_{OAFM}$	$x_{iter}$
0	0.25	0.2500000099	0.25
1/2	0.3871461205	0.3871462684	0.3871461097
1	0.5547748969	0.5547749864	0.5547748652
3/2	0.7604269013	0.7604266715	0.7604267852
2	1.0050873071	1.0050873871	1.0050857278
5/2	1.2742254393	1.2742254977	1.2742081489
3	1.5271576554	1.5271575136	1.5272332007
7/2	1.6978024263	1.6978026989	1.7009599136
4	1.7246494125	1.7246491591	1.7536803864
9/2	1.5967155859	1.5967157642	1.7306827097
5	1.3631319493	1.3631316176	1.7128681253

Appendix A

Example 1. ${\bar{u}}_{O A F M}$ is the OAFM-solution of Equations (61) and (72) for initial data $x 0 = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0$ , $c = 0.75$ , $d = - 0.45 < 0$ , $k = 0.65$ . The optimal parameters that appear in Equation (93) for different values for the index number $N_{m a x} = 25$ are as follows:

(A1)

\begin{matrix} ω_{0} = 0.1004121488, & B_{1} = 7.3119258362, & B_{2} = 42.6993026598, \\ B_{3} = 59.4213749369, & B_{4} = - 17.3757871188, & B_{5} = - 129.4179243917, \\ B_{6} = - 125.2333676742, & B_{7} = 19.4746605558, & B_{8} 142.5824124930, \\ B_{9} = 111.1428236138, & B_{10} = - 13.6338957406, & B_{11} = - 83.8801784433, \\ B_{12} = - 53.7620365001, & B_{13} = 5.9338443476, & B_{14} = 28.1051309583, \\ B_{15} = 14.5871505067, & B_{16} = - 1.5120260034, & B_{17} = - 5.1175424653, \\ B_{18} = - 2.0487152810, & B_{19} = 0.1980873163, & B_{20} = 0.4287108517, \\ B_{21} = 0.1174067649, & B_{22} = - 0.0098895124, & B_{23} = - 0.0102936430, \\ B_{24} = - 0.0012404443, & B_{25} = 0.0000663868, & C_{1} = 17.2240423624, \\ C_{2} = 12.2741255032, & C_{3} = - 50.8493401442, & C_{4} = - 110.7219385490, \\ C_{5} = - 52.9308854253, & C_{6} = 95.0719694550, & C_{7} = 161.0009809177, \\ C_{8} = 61.3272239797, & C_{9} = - 82.8291813381, & C_{10} = - 115.2112897092, \\ C_{11} = - 35.8193830285, & C_{12} = 40.5991539215, & C_{13} = 46.3589037273, \\ C_{14} = 11.5929461208, & C_{15} = - 11.3627188386, & C_{16} = - 10.3622108429, \\ C_{17} = - 1.9903649181, & C_{18} = 1.6740284016, & C_{19} = 1.1418639879, \\ C_{20} = 0.1515033205, & C_{21} = - 0.1034801174, & C_{22} = - 0.0446204683, \\ C_{23} = - 0.0029550038, & C_{24} = 0.0013201742, & C_{25} = 0.0001516801; \end{matrix}

Example 2. ${\bar{u}}_{O A F M}$ is the OAFM-solution of Equations (63) and (72) for initial data $x 0 = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.75$ , $c = 0$ , $d = 0.45 > 0$ , $k = 0.65$ The optimal parameters that appear in Equation (93) for different values for the index number $N_{m a x} = 25$ are as follows:

(A2)

\begin{matrix} ω_{0} = 0.1047139130, & B_{1} = 0.3086318952, & B_{2} = 0.5313403493, \\ B_{3} = 0.2993192208, & B_{4} = 0.3726412184, & B_{5} = - 0.6877516775, \\ B_{6} = - 1.5562906173, & B_{7} = - 1.2682734226, & B_{8} = 0.2192571336, \\ B_{9} = 1.4155345390, & B_{10} = 1.2923201865, & B_{11} = 0.2423126050, \\ B_{12} = - 0.6050546289, & B_{13} = - 0.6628999635, & B_{14} = - 0.2376657689, \\ B_{15} = 0.1150081078, & B_{16} = 0.1774048561, & B_{17} = 0.0802770323, \\ B_{18} = - 0.0028859700, & B_{19} = - 0.0222102008, & B_{20} = - 0.0110686269, \\ B_{21} = - 0.0013674245, & B_{22} = 0.0009211829, & B_{23} = 0.0004429820, \\ B_{24} = 0.0000601901, & B_{25} = - 3.1883 \times 10^{- 6}, & C_{1} = 0.2756822273, \\ C_{2} = 0.3797254820, & C_{3} = - 0.3435526898, & C_{4} = - 0.9666205940, \\ C_{5} = - 1.2332443309, & C_{6} = - 0.3082072428, & C_{7} = 1.2095840157, \\ C_{8} = 1.7317101440, & C_{9} = 0.8061022308, & C_{10} = - 0.5854077453, \\ C_{11} = - 1.1427897428, & C_{12} = - 0.6795263469, & C_{13} = 0.0616404935, \\ C_{14} = 0.3944407927, & C_{15} = 0.2768101108, & C_{16} = 0.0441247526, \\ C_{17} = - 0.0662765534, & C_{18} = - 0.0548639333, & C_{19} = - 0.0147978533, \\ C_{20} = 0.0039435754, & C_{21} = 0.0044008656, & C_{22} = 0.0012759316, \\ C_{23} = 0.0000114865, & C_{24} = - 0.0000729093, & C_{25} = - 0.0000113965; \end{matrix}

Example 3. ${\bar{u}}_{O A F M}$ is the OAFM-solution of Equations (65) and (72) for initial data $x 0 = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = 0.85 > 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ The optimal parameters that appear in Equation (101) for different values for the index number $N_{m a x} = 25$ are as follows:

(A3)

\begin{matrix} ω_{0} = 0.0837434277, & B_{1} = - 3.1515787737, & B_{2} = - 2.1385161706, \\ B_{3} = 1.5441444917, & B_{4} = 3.8932367709, & B_{5} = 3.8481666199, \\ B_{6} = 1.0761566007, & B_{7} = - 2.4309166849, & B_{8} = - 3.6161186440, \\ B_{9} = - 1.8543214418, & B_{10} = 0.6820900822, & B_{11} = 1.7571483521, \\ B_{12} = 1.1270525157, & B_{13} = 0.0338905245, & B_{14} = - 0.4780229448, \\ B_{15} = - 0.3541838012, & B_{16} = - 0.0677800531, & B_{17} = 0.0669893745, \\ B_{18} = 0.0563801632, & B_{19} = 0.0144682625, & B_{20} = - 0.0036900316, \\ B_{21} = - 0.0037377194, & B_{22} = - 0.0009404535, & B_{23} = 0.00002248507, \\ B_{24} = 0.0000533823, & B_{25} = 7.112 \times 10^{- 6}, & C_{1} = 2.2726517301, \\ C_{2} = 1.9751332671, & C_{3} = 3.3286075863, & C_{4} = 2.0229548103, \\ C_{5} = - 1.8249466878, & C_{6} = - 4.2698346367, & C_{7} = - 3.2137310100, \\ C_{8} = 0.0036033605, & C_{9} = 2.4211138073, & C_{10} = 2.3520437624, \\ C_{11} = 0.6494932824, & C_{12} = - 0.7527848517, & C_{13} = - 0.9315904606, \\ C_{14} = - 0.3700048736, & C_{15} = 0.1069448162, & C_{16} = 0.2011346233, \\ C_{17} = 0.0918006354, & C_{18} = - 0.0008263369, & C_{19} = - 0.0211728089, \\ C_{20} = - 0.0099105326, & C_{21} = - 0.0009913894, & C_{22} = 0.0007842430, \\ C_{23} = 0.0003178325, & C_{24} = 0.0000328774, & C_{25} = - 2.472 \times 10^{- 6}; \end{matrix}

Example 4. ${\bar{u}}_{O A F M}$ is the OAFM-solution of Equations (65) and (72) for initial data $x 0 = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = - 0.5 < 0$ , $b = - 0.85 < 0$ , $c = 0.75 > 0$ , $d = - 0.45 < 0$ , $k = 0.65$ The optimal parameters that appear in Equation (101) for different values for the index number $N_{m a x} = 25$ are as follows:

(A4)

\begin{matrix} ω_{0} = 0.0913911457, & B_{1} = - 226.8540043327, & B_{2} = 265.6605062478, \\ B_{3} = 931.7388690548, & B_{4} = - 896.2918959290, & B_{5} = - 516.7204243651, \\ B_{6} = - 160.0604662083, & B_{7} = 243.6332698017, & B_{8} = 522.5414937630, \\ B_{9} = 208.2575046181, & B_{10} = 108.1318190204, & B_{11} = - 428.8940069255, \\ B_{12} = - 185.5328461326, & B_{13} = - 335.7211183847, & B_{14} = 395.3573025962, \\ B_{15} = 98.5572914924, & B_{16} = 427.5090910475, & B_{17} = - 487.5210772261, \\ B_{18} = - 253.9743791380, & B_{19} = 383.5827206995, & B_{20} = - 53.7859131483, \\ B_{21} = - 58.3907261312, & B_{22} = 18.7014806950, & B_{23} = 0.6708570897, \\ B_{24} = - 0.6214611435, & B_{25} = 0.0261129487, & C_{1} = 9.9338943296, \\ C_{2} = 637.8121987193, & C_{3} = - 688.1561252406, & C_{4} = - 861.0433936404, \\ C_{5} = 562.0223108032, & C_{6} = 258.9221028658, & C_{7} = 655.4226597511, \\ C_{8} = - 172.0672563402, & C_{9} = - 85.1584803811, & C_{10} = - 568.5143050312, \\ C_{11} = - 71.0839617212, & C_{12} = - 65.2161590187, & C_{13} = 431.2983597185, \\ C_{14} = 218.7751370601, & C_{15} = 22.6481737014, & C_{16} = - 247.5848707888, \\ C_{17} = - 477.7488881146, & C_{18} = 539.7010052076, & C_{19} = 31.9744026607, \\ C_{20} = - 183.8288774652, & C_{21} = 45.2898700808, & C_{22} = 10.9746163940, \\ C_{23} = - 4.5615175082, & C_{24} = 0.1535303109, & C_{25} = 0.0355506432; \end{matrix}

Example 5. ${\bar{u}}_{O A F M}$ is the OAFM-solution of Equations (67) and (72) for initial data $x 0 = 0.25$ , $y_{0} = 0.5$ , $z_{0} = 1.5$ ; physical constants: $a = 0.5 > 0$ , $b = 0.85 > 0$ , $c = - 0.75 < 0$ , $d = - 0.45 < 0$ , $k = 0.65$ The optimal parameters that appear in Equation (101) for different values for the index number $N_{m a x} = 25$ are as follows:

(A5)

\begin{matrix} ω_{0} = 0.0522026653, & B_{1} = - 504.5389490182, & B_{2} = - 41.4720576342, \\ B_{3} = 858.8736751708, & B_{4} = - 12.2210727195, & B_{5} = 71.623966633, \\ B_{6} = - 362.5575056477, & B_{7} = - 215.7652450323, & B_{8} = 6.2003889342, \\ B_{9} = 24.7700733919, & B_{10} = 311.2777732871, & B_{11} = - 21.7436380405, \\ B_{12} = 22.9111383651, & B_{13} = - 178.1435782892, & B_{14} = - 151.0646749115, \\ B_{15} = 236.2925981990, & B_{16} = 34.4139955889, & B_{17} = - 102.2086336648, \\ B_{18} = 10.0220324566, & B_{19} = 18.0215765638, & B_{20} = - 3.8643195277, \\ B_{21} = - 1.1700305856, & B_{22} = 0.3332304667, & B_{23} = 0.0155027326, \\ B_{24} = - 0.0063196841, & B_{25} = 0.0000729754, & C_{1} = - 101.1656550811, \\ C_{2} = 913.0001926763, & C_{3} = - 63.6819904517, & C_{4} = - 407.5948692447, \\ C_{5} = - 200.6539424118, & C_{6} = - 279.4998547021, & C_{7} = 284.4709961065, \\ C_{8} = 80.2409804205, & C_{9} = 252.6380896874, & C_{10} = - 37.0805456848, \\ C_{11} = - 172.6372139042, & C_{12} = - 65.0999791648, & C_{13} = - 142.4634271405, \\ C_{14} = 245.2318200810, & C_{15} = 94.9140212750, & C_{16} = - 174.4664734121, \\ C_{17} = 0.2055439256, & C_{18} = 48.0298579498, & C_{19} = - 7.9468423486, \\ C_{20} = - 5.2949867586, & C_{21} = 1.3347963584, & C_{22} = 0.1792358943, \\ C_{23} = - 0.0579182983, & C_{24} = - 0.0001499925, & C_{25} = 0.0003270484 . \end{matrix}

References

1. Xu, M.; Shi, S.; Huang, K. The connection between the dynamical properties of 3D systems and the image of the energy-Casimir mapping. Discrete Contin. Dyn. Syst. AIMS; 2024; 44, pp. 791-807. [DOI: https://dx.doi.org/10.3934/dcds.2023126]

2. Lazureanu, C. On the Hamilton–Poisson realizations of the integrable deformations of the Maxwell-Bloch equations. C. R. Math.; 2017; 355, pp. 596-600. [DOI: https://dx.doi.org/10.1016/j.crma.2017.04.002]

3. Lazureanu, C. Hamilton–Poisson Realizations of the Integrable Deformations of the Rikitake System. Adv. Math. Phys.; 2017; 2017, 4596951.

4. Tudoran, R.M.; Girban, A. On the completely integrable case of the Rössler system. J. Math. Phys.; 2012; 53, 052701. [DOI: https://dx.doi.org/10.1063/1.4708621]

5. Gümral, H.; Nutku, Y. Poisson structure of dynamical systems with three degrees of freedom. J. Math. Phys.; 1993; 34, 5691. [DOI: https://dx.doi.org/10.1063/1.530278]

6. Lazureanu, C.; Cho, J. On a family of Hamilton–Poisson Jerk System. Mathematics; 2024; 12, 1260. [DOI: https://dx.doi.org/10.3390/math12081260]

7. Tudoran, R.M.; Aron, A.; Nicoara, S. On a Hamiltonian version of the Rikitake system. SIAM J. Appl. Dyn. Syst.; 2009; 8, pp. 454-479. [DOI: https://dx.doi.org/10.1137/080728822]

8. Tudoran, R.M.; Tudoran, R.A. On a large class of three-dimensional Hamiltonian systems. J. Math. Phys.; 2009; 50, 012703. [DOI: https://dx.doi.org/10.1063/1.3068405]

9. Lazuranu, C.; Binzar, T. Symmtries of some classes of dynamical systems. J. Nonlinear Math. Phy.; 2015; 22, pp. 265-274.

10. Lazureanu, C.; Hedrea, C.; Petrisor, C. On a deformed version of the two-disk dynamo system. Appl. Math.; 2021; 66, pp. 345-372. [DOI: https://dx.doi.org/10.21136/AM.2021.0303-19]

11. Giné, J.; Sinelshchikov, D.I. On the geometric and analytical properties of the anharmonic oscillator. Commun. Nonlinear Sci.; 2024; 131, 107875. [DOI: https://dx.doi.org/10.1016/j.cnsns.2024.107875]

12. 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.

13. Lazureanu, C.; Binzar, T. On the symmetries of a Rikitake type system. C. R. Math. Acad. Sci. Paris; 2012; 350, pp. 529-533.

14. Xu, M.; Shi, S.; Huang, K. On the integrable stretch-twist-fold flow: Bi-Hamiltonian structures and global dynamics. J. Math. Phys.; 2024; 65, 022704. [DOI: https://dx.doi.org/10.1063/5.0185673]

15. Huang, K.; Shi, S.; Xu, Z. Integrable deformations, bi-Hamiltonian structures and nonintegrability of a generalized Rikitake system. Int. J. Geom. Methods Mod. Phys.; 2019; 16, 1950059. [DOI: https://dx.doi.org/10.1142/s0219887819500592]

16. Biggs, R.; Remsing, C.C. Quadratic Hamilton–Poisson Systems in Three Dimensions: Equivalence, Stability, and Integration. Acta Appl. Math.; 2017; 148, pp. 1-59. [DOI: https://dx.doi.org/10.1007/s10440-016-0074-1]

17. Olver, P.J. Applications of Lie Groups to Differential Equations; Springer: New York, NY, USA, 1986.

18. Krishnaprasad, P.S. Lie–Poisson structures, dual-spin spacecraft and asymptotic stability. Nonlinear Anal. Theory Methods Appl.; 1985; 9, pp. 1011-1035.

19. David, D.; Holm, D.D. Multiple Lie–Poisson Structures, Reductions, and Geometric Phases for the Maxwell-Bloch Travelling Wave Equations. J. Nonlinear Sci.; 1992; 2, pp. 241-262.

20. Rikitake, T. Oscillations of a system of disk dynamos. Proc. Camb. Philos. Soc.; 1958; 54, pp. 89-105. [DOI: https://dx.doi.org/10.1017/S0305004100033223]

21. Pikovskii, A.S.; Rabinovich, M.I. Stochastic behavior of dissipative systems. Soc. Sci. Rev. C Math. Phys. Rev.; 1981; 2, pp. 165-208.

22. Chen, H.-K.; Lee, C.-I. Anti-control of chaos in rigid body motion. Chaos Soliton Fract.; 2004; 21, pp. 957-965. [DOI: https://dx.doi.org/10.1016/j.chaos.2003.12.034]

23. Ene, R.-D.; Pop, N. Semi-Analytical Solutions for the Qi–Type Dynamical System. Symmetry; 2024; 16, 1578. [DOI: https://dx.doi.org/10.3390/sym16121578]

24. 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.

25. Ene, R.-D.; Pop, N.; Lapadat, M. Approximate Closed-Form Solutions for the Rabinovich System via the Optimal Auxiliary Functions Method. Symmetry; 2022; 14, 2185. [DOI: https://dx.doi.org/10.3390/sym14102185]

26. 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]

27. Daftardar-Gejji, V.; Jafari, H. An iterative method for solving nonlinear functional equations. J. Math. Anal. Appl.; 2006; 316, pp. 753-763.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Symmetries and Closed-Form Solutions for Some Classes of Dynamical Systems

Content area

Abstract

Full text

1. Introduction

2. Symmetries, Equilibrium Points, Periodic Orbits,Homoclinic Orbits, Heteroclinic Orbits

2.1. Case 1— ac−bd≠0

2.2. Case 2— ac−bd=0

3. Closed-Form Solutions

4. The Optimal Auxiliary Functions Method (OAFM)

4.1. Main Steps of the OAFM

4.2. Semi-Analytical Solutions via the OAFM Procedure

5. Numerical Results and Validation

6. Conclusions

2.1. Case 1— $a c - b d \neq 0$

2.2. Case 2— $a c - b d = 0$