Semi-Analytical Closed-Form Solutions for the Rikitake-Type System through the Optimal Homotopy Perturbation Method

Abstract

The goal of this work is to build semi-analytical solutions of the Rikitake-type system by means of the optimal homotopy perturbation method (OHPM) using only two iterations. The chaotic behaviors are excepted. By taking into consideration the geometrical properties of the Rikitake-type system, the closed-form solutions can be established. The obtained solutions have a periodical behavior. These geometrical properties allow reducing the initial system to a second-order nonlinear differential equation. The latter equation is solved analytically using the OHPM procedure. The validation of the OHPM method is presented for three cases of the physical parameters. The advantages of the OHPM technique, such as the small number of iterations (the efficiency), the convergence control (in the sense that the semi-analytical solutions are approaching the exact solution), and the writing of the solutions in an effective form, are shown graphically and with tables. The accuracy of the results provides good agreement between the analytical and corresponding numerical results. Other dynamic systems with similar geometrical properties could be successfully solved using the same procedure.

1. Introduction

The first study of the Rikitake system [1] examined the behavior of two disk dynamos coupled one to another in relation to Earth’s magnetic field. These studies were continued by Steeb [2].

Some relevant dynamical and geometrical properties of the Rikitake system were examined in [3] by reporting the stability of the equilibria points, the existence of the symmetries [4], the Hamilton–Poisson realization of the Rikitake-type system with control, the existence of the periodic orbits, the Lax formulation [5,6,7], and the integrable deformations [8].

Other dynamical systems with geometrical properties are the Rabinovich system, studied in [9,10,11], the Kermack–McKendrick system, investigated in [12], and the Rossler system, explored in [13]. In addition, a Lotka–Volterra-type system was studied in [14], and a T system with an integrable deformation was examinated in [15].

Beyond the geometrical properties, many of the dynamical systems have chaotic behaviors and are found in practical applications in engineering (e.g., Chua’s circuits [16,17], the turbulence in the atmosphere and oceans [18], and plasma control systems [19]) or in medicine (e.g., multi-scale synthesis of Alzheimer’s pathogenesis [20], neuroethology [21], and a model of stress as a dynamical system [22]).

Other analytical methods for solving nonlinear differential problems were developed, namely the variational iteration method (VIM) [23,24], the multiple scales technique [25], the homotopy perturbation method (HPM) [26], the modified Laplace transform homotopy perturbation method (MLT-HPM) [27], the homotopy analysis method (HAM) [28], the least squares differential quadrature method [29], the polynomial least squares method [30], the function method [31], the optimal iteration parametrization method (OIPM) [32], the optimal variational iteration method (OVIM) [33], the optimal homotopy asymptotic method (OHAM) [34,35,36], the and optimal homotopy perturbation method (OHPM) [37,38]. The modified homotopy perturbation method (MHPM) is applied to solve sine-Gordon and coupled sine-Gordon equations [39], among other applications [40,41,42]. The modified homotopy analysis method (MHAM) is used to find appropriate solutions to Zakharov–Kuznetsov equations [43], among other uses [44]. To build some new exact solitary wave solutions for the sixth-order Boussinesq equation, which plays an important role in mathematical physics, engineering sciences and applied mathematics, the exp-function method is proposed [45], among other applications [46,47]. The modified exp-function method explores the strain wave equation in micro-structured solids, which is applied in solid physics [48], among other applications [49,50,51,52,53].

The structure of this paper is as follows. In Section 2, we present the geometrical properties and the closed-form solutions of the Rikitake-type system. Section 3 is dedicated to the basic ideas of the OHPM technique. Section 4 is focused on building semi-analytical solutions using the OHPM method. Section 5 provides the numerical results. The last section is dedicated to the conclusions.

2. The Rikitake-Type System

2.1. Global Analytic First Integrals and Hamilton–Poisson Realizations

The Rikitake-type system has the following form (see [5,6]):

$$\left\{\begin{array}{c}\dot{x}=kyz+\alpha y\hfill \\ \dot{y}=xz-\alpha x\hfill \\ \dot{z}=-xy\hfill \end{array}\right.,k,\alpha \in \mathbb{R},$$

(1)

where k is the tuning parameter and the state variables x, y, and z depend on $t>0$ and $\dot{x}=\frac{dx}{dt}$.

For the physical system describing the behavior of two disk dynamos coupled to one another in relation to Earth’s magnetic field, in agreement with [54], the state variables x, y, and z are expressed as

$$\left\{\begin{array}{c}x=I_1·{\left(\frac{M}{G}\right)}^{\frac{1}{2}}\hfill \\ y=I_2·{\left(\frac{M}{G}\right)}^{\frac{1}{2}}\hfill \\ z={\Omega }_0·{\left(\frac{CM}{GL}\right)}^{\frac{1}{2}}\hfill \end{array}\right.,$$

(2)

where $I_1$ and $I_2$ are the currents, C is moment of inertia of a dynamo about its axis, M is the mutual inductance between the dynamo circuits, L is the self-inductance, G is the torque which describes the motion of the dynamos in the report of Ohmic losses in the coils and disks, and ${\Omega }_0$ is the angular velocity.

Remark 1.

Some geometrical properties of the system in Equation (1) are the following:

• This system admits a symmetry with respect to the $Oz$ axis (invariant to the transformation $(x,y,z)\to (-x,-y,z)$ for $\alpha \ne 0$, $k\in \mathbb{R}$) and symmetries with respect to all coordinate axes for $\alpha =0$ and $k\in \mathbb{R}$;
• The system has a Hamilton–Poisson realization, with the Hamiltonian and the Casimir given by $H(x,y,z)=\frac{1}{2}(x^2+y^2)+\frac{1+k}{2}{z}^2$ and $C(x,y,z)=\frac{1}{4}(y^2-x^2)-\alpha z+\frac{1-k}{4}{z}^2$, respectively, for $\alpha \ne 0$, $k\in \mathbf{R}$, $H(x,y,z)=\frac{1}{2}(x^2+y^2)+z^2$, and $C(x,y,z)=\frac{1}{2}(x^2-y^2)$ for the special case where $\alpha =0$ and $k=1$;
• The existence of the periodic orbits ensures that the solutions have periodic behaviors.

Remark 2.

Let the initial conditions be

$$x\left(0\right)=x_0,y\left(0\right)=y_0,z\left(0\right)=z_0.$$

(3)

Then, the phase curves of dynamics (1) are the intersections of the following surfaces

$$\left\{\begin{array}{c}\frac{1}{2}(x^2+y^2)+\frac{1+k}{2}{z}^2=\frac{1}{2}(x_0^2+y_0^2)+\frac{1+k}{2}{z}_0^2\hfill \\ \frac{1}{4}(y^2-x^2)-\alpha z+\frac{1-k}{4}{z}^2=\frac{1}{4}(y_0^2-x_0^2)-\alpha z_0+\frac{1-k}{4}{z}_0^2\hfill \end{array}\right.,\mathrm{for}\alpha \ne 0,k\in \mathbb{R},$$

(4)

and

$$\left\{\begin{array}{c}\frac{1}{2}(x^2+y^2)+z^2=\frac{1}{2}(x_0^2+y_0^2)+z_0^2\hfill \\ x^2-y^2=x_0^2-y_0^2\hfill \end{array}\right.,\mathrm{for}\alpha =0,k=1,$$

(5)

respectively.

2.2. Closed-Form Solutions

For three cases of the selected physical parameters $\alpha$, k, the closed-form solutions of the system in Equation (1) are obtained.

(i)

For $\alpha >0$, $k>0$,

The transformations

$$\left\{\begin{array}{c}y\left(t\right)=R·sin\left(u\right(t\left)\right)\hfill \\ z\left(t\right)=\alpha +R·cos\left(u\right(t\left)\right)\hfill \end{array}\right.,$$

(6)

where $R=\sqrt{y_0^2+{(z_0-\alpha )}^2}$, $u\left(t\right)$ is an unknown smooth function, giving the closed-form solutions.

The third equation from Equation (1) yields

$$x\left(t\right)=\dot{u}\left(t\right).$$

(7)

From Equation (1), we obtain

$$\ddot{u}\left(t\right)-\frac{k}{2}{R}^2·sin(2·u\left(t\right))-(k+1)R\alpha sin\left(u\left(t\right)\right)=0.$$

(8)

The initial conditions $u\left(0\right)$ and $\dot{u}\left(0\right)$ obtained from Equations (3), (6) and (7) are

$$u\left(0\right)=arctan\frac{y_0}{z_0-\alpha },\dot{u}\left(0\right)=x_0,$$

(9)

with $z_0\ne \alpha$.

Remark 3.

For $u\left(t\right)$ the exact solution of the problem in Equations (8) and (9), the closed-form solution of the system in Equation (1) results from Equations (6) and (7).

If $\overline{u}\left(t\right)$ is an semi-analytical solution of the problem Equations (8) and (9), then the Equations (6) and (7) provide semi-analytical closed-form solutions of the system in Equation (1) denoted by $\overline{x}\left(t\right)$, $\overline{y}\left(t\right)$ and $\overline{z}\left(t\right)$, respectively.

(ii)

In the case $\alpha \ne 0$, $k<0$,

The closed-form solutions can be written as

$$\left\{\begin{array}{c}y\left(t\right)=R·\frac{2·u\left(t\right)}{1+u^2\left(t\right)}\hfill \\ z\left(t\right)=\alpha +R·\frac{1-u^2\left(t\right)}{1+u^2\left(t\right)}\hfill \end{array}\right.,$$

(10)

where $R=\sqrt{{(z_0-\alpha )}^2+y_0^2}$.

Equation (1) yields

$$x\left(t\right)=\frac{2\dot{u}\left(t\right)}{1+u^2\left(t\right)}.$$

(11)

The following nonlinear problem gives the unknown function $u\left(t\right)$

$$\left\{\begin{array}{c}\ddot{u}\left(t\right)·(1+u^2\left(t\right))-2u\left(t\right)·{\left(\dot{u}\left(t\right)\right)}^2-k{R}^2·u\left(t\right)·(1-u^2\left(t\right))-\hfill \\ -(k+1)\alpha ·R·u\left(t\right)·(1+u^2\left(t\right))=0\hfill \\ u\left(0\right)=\frac{|y_0|}{|z_0-\alpha +R|},\dot{u}\left(0\right)=\frac{x_0}{2}·\left(1+u^2\left(0\right)\right),\hfill \end{array}\right.$$

(12)

with $z_0-\alpha +R\ne 0$.

(iii)

A remarkable case is $\alpha =0$, $k=1$.

The closed-form solutions could have the form

$$\left\{\begin{array}{c}x\left(t\right)=R·\frac{2·u\left(t\right)}{1+u^2\left(t\right)}\hfill \\ z\left(t\right)=-\alpha +R·\frac{1-u^2\left(t\right)}{1+u^2\left(t\right)}\hfill \end{array}\right.,$$

(13)

where $R=\sqrt{{(z_0+\alpha )}^2+y_0^2}$. Then the third equation from Equation (1) yields

$$y\left(t\right)=\frac{2\dot{u}\left(t\right)}{1+u^2\left(t\right)},$$

(14)

where $u\left(t\right)$ is the solution of the following nonlinear problem

$$\left\{\begin{array}{c}\ddot{u}\left(t\right)-\frac{2u\left(t\right)}{1+u^2\left(t\right)}·{\left(\dot{u}\left(t\right)\right)}^2-R^2·u\left(t\right)·\frac{1-u^2\left(t\right)}{1+u^2\left(t\right)}+2\alpha ·R·u\left(t\right)=0\hfill \\ u\left(0\right)=\frac{|x_0|}{|R+(z_0+\alpha \left)\right|},\dot{u}\left(0\right)=\frac{y_0}{2}·\left(1+u^2\left(0\right)\right).\hfill \end{array}\right.$$

(15)

Remark 4.

An analytical approximate solution of Equation (1) is said to be a semi-analytical solution and will be denoted by $\overline{x}\left(t\right)$, $\overline{y}\left(t\right)$, $\overline{z}\left(t\right)$ or ${\overline{x}}_{OHPM}\left(t\right)$, ${\overline{y}}_{OHPM}\left(t\right)$, ${\overline{z}}_{OHPM}\left(t\right)$ by the OHPM technique.

In this work, Equations (8), (9), (12) and (41) are semi-analytically solved using the optimal homotopy perturbation method (OHPM) technique.

A semi-analytical solution of the auxiliary problem given by Equations (8) or (9) or (12) or (41) will be denoted by $\overline{u}$ or ${\overline{u}}_{OHPM}$.

3. The Steps of the OHPM Technique

This section describes the ideas of the OHPM for the second-order nonlinear differential Equation [35,37,38]

$$\mathcal{L}\left[f,f^{\prime },f^{″},t\right]+g\left(t\right)+\mathcal{N}\left[f,f^{\prime },f^{″},t\right]=0,t\in D,$$

(16)

with the initial/boundary conditions

$$\mathcal{B}\left(f\left(t\right),\frac{df\left(t\right)}{dt}\right)=0,t\in \partial D,$$

(17)

where $\mathcal{L}$, $\mathcal{N}$ and $\mathcal{B}$ are the linear, nonlinear and boundary operators, respectively, and D is domain of interest and $g\left(t\right)$ is an arbitrary continuous function.

Let $f_0\left(t\right)$ be an initial approximation of the unknown function $f\left(t\right)$, such as

$$\mathcal{L}\left[f_0\left(t\right)\right]+g\left(t\right)=0,\mathcal{B}\left(f_0\left(t\right),\frac{d{f}_0\left(t\right)}{dt}\right)=0,$$

(18)

and $p\in \left[0,1\right]$ be an embedding parameter. In the following, we consider a function $\Phi (t,p,C_{k_1})$ in the particular form

$$\Phi (t,p,C_{k_1})=f_0\left(t\right)+p{f}_1(t,C_{k_1})+p^2{f}_2(t,C_{k_1}),$$

(19)

where $C_{k_1}$, $k_1=1,2,\dots$ are arbitrary parameters which will be determined later, $f_1$ is named the first approximation and $f_2$ is named the second approximation. Obviously, when $p=0$ and $p=1$ it holds that

$$\Phi (t,0,C_{k_1})=f_0\left(t\right),\Phi (t,1,C_{k_1})=f(t,C_{k_1}),$$

(20)

respectively. Thus, as p increases from 0 to 1, the solution $\Phi (t,p,C_{k_1})$ varies from $f_0\left(\eta \right)$ to the solution $f(t,C_{k_1})$. Consequently, the second-order approximate solution can be written as

$$\overline{f}(t,C_{k_1})=f_0\left(t\right)+f_1(t,C_{k_1})+f_2(t,C_{k_1}).$$

(21)

The initial/boundary conditions become

$$\mathcal{B}\left(\overline{f}(t,C_{k_1}),\frac{d\overline{f}(t,C_{k_1})}{dt}\right)=0.$$

(22)

For nonlinear Equation (16), the following family of equations is performed:

$$\begin{array}{c}\mathcal{H}\left(\Phi (t,p,C_{k_1})\right)=\mathcal{L}\left(f,f^{\prime },f^{″},t\right)+g\left(t\right)+p\mathcal{N}\left(f,f^{\prime },f^{″},t\right)=0,\hfill \end{array}$$

(23)

with the properties

$$\mathcal{H}\left(\Phi (t,0,C_{k_1})\right)=L\left[f_0\left(t\right)\right]+g\left(t\right)=0,$$

(24)

$$\mathcal{H}\left(\Phi (t,1,C_{k_1})\right)=-\left[L\left(\overline{f}(t,C_{k_1})\right)+\mathcal{N}\left(\overline{f}(t,C_{k_1})\right)\right]=0.$$

(25)

Series expansion of the nonlinear operator $\mathcal{N}$ with respect to the parameter p yields

$$\begin{array}{c}\hfill \begin{array}{c}\mathcal{N}\left(\overline{f},{\overline{f}}^{\prime },{\overline{f}}^{″},t\right)=\mathcal{N}\left(f_0,f_0^{\prime },f_0^{″},t\right)+p\left[f_1{\mathcal{N}}_{\overline{f}}\left(f_0,f_0^{\prime },f_0^{″},t\right)\right.\hfill \\ \left.+f_1^{\prime }{\mathcal{N}}_{{\overline{f}}^{\prime }}\left(f_0,f_0^{\prime },f_0^{″},t\right)+f_1^{″}{\mathcal{N}}_{{\overline{f}}^{″}}\left(f_0,f_0^{\prime },f_0^{″},t\right)\right]+p^2\left(f_2{\mathcal{N}}_{\overline{f}}+f_2^{\prime }{\mathcal{N}}_{{\overline{f}}^{\prime }}+\dots \right),\hfill \end{array}\end{array}$$

(26)

where $F_{\overline{u}}=\frac{\partial F}{\partial \overline{u}}$.

The homotopy relation given by Equation (23) becomes

$$\begin{array}{c}\mathcal{H}\left(\Phi (t,p,C_{k_1})\right)=\mathcal{L}\left(f,f^{\prime },f^{″},t\right)+g\left(t\right)+p·\left\{\mathcal{N}\left(f_0,f_0^{\prime },f_0^{″},t\right)+p\left[f_1{\mathcal{N}}_{\overline{f}}\left(f_0,f_0^{\prime },f_0^{″},t\right)\right.\right.\hfill \\ \left.+f_1^{\prime }{\mathcal{N}}_{{\overline{f}}^{\prime }}\left(f_0,f_0^{\prime },f_0^{″},t\right)+f_1^{″}{\mathcal{N}}_{{\overline{f}}^{″}}\left(f_0,f_0^{\prime },f_0^{″},t\right)\right]\left.+p^2\left(f_2{\mathcal{N}}_{\overline{f}}+f_2^{\prime }{\mathcal{N}}_{{\overline{f}}^{\prime }}+\dots \right)\right\}=0.\hfill \end{array}$$

(27)

Let $H_{i,j}(t,C_{k_1})$, $i=1,2,3$, $j=0,1,2$ be a number of unknown auxiliary functions that depend on the variable t and some parameters $C_{k_1}$, $k_1=1,2,\dots ,s$, where $s\ge 1$ is a given integer number. The following new homotopy relation is proposed:

$$\begin{array}{c}\mathcal{H}\left(\Phi (t,p,C_{k_1})\right)=\mathcal{L}\left[f_0\left(t\right)\right]+g\left(t\right)+p·H_{1,0}(t,C_{k_1})·\mathcal{N}\left[f_0\left(t\right)\right]\hfill \\ +p^2·\left[H_{2,0}(t,C_{k_1})·f_1·{\mathcal{N}}_{\overline{f}}\left[f_0\left(t\right)\right]+H_{2,1}(t,C_{k_1})·{f^{\prime }}_1·{\mathcal{N}}_{{\overline{f}}^{\prime }}\left[f_0\left(t\right)\right]+H_{2,2}(t,C_{k_1})·{f^{″}}_1·{\mathcal{N}}_{{\overline{f}}^{″}}\left[f_0\left(t\right)\right]\right]\hfill \\ +p^3·\left[H_{3,0}(t,C_{k_1})·f_2·{\mathcal{N}}_{\overline{f}}\left[f_0\left(t\right)\right]+H_{3,1}(t,C_{k_1})·{f^{\prime }}_2·{\mathcal{N}}_{{\overline{f}}^{\prime }}\left[f_0\left(t\right)\right]+H_{3,2}(t,C_{k_1})·{f^{″}}_2·{\mathcal{N}}_{{\overline{f}}^{″}}\left[f_0\left(t\right)\right]\right].\hfill \end{array}$$

(28)

In the following, the terms in $p^3$ are neglected.

From Equations (19) and (21), we obtain

$$\Phi (t,0,C_{k_1})=f_0\left(t\right),\Phi (t,1,C_{k_1})=\overline{f}(t,C_{k_1}).$$

(29)

Now, equating the coefficients of $p^0$, $p^1$ and $p^2$ into Equation (28), it holds that $f_0$ is given by Equation (18), and the first approximation $f_1(t,C_{k_1})$ is obtained from the following equation:

$$\mathcal{L}\left[f_1(t,C_{k_1})\right]+H_{1,0}(t,C_{k_1})·\mathcal{N}\left[f_0\left(t\right)\right]=0,B\left(f_1,\frac{d{f}_1}{dt}\right)=0.$$

(30)

and the second approximation $f_2(t,C_{k_1})$ is obtained from the following equation:

$$\begin{array}{c}\mathcal{L}\left[f_2(t,C_{k_1})\right]+H_{2,0}(t,C_{k_1})·f_1·{\mathcal{N}}_{\overline{f}}\left[f_0\left(t\right)\right]+H_{2,1}(t,C_{k_1})·{f^{\prime }}_1·{\mathcal{N}}_{{\overline{f}}^{\prime }}\left[f_0\left(t\right)\right]\hfill \\ +H_{2,2}(t,C_{k_1})·{f^{″}}_1·{\mathcal{N}}_{{\overline{f}}^{″}}\left[f_0\left(t\right)\right]=0,B\left(f_2,\frac{d{f}_2}{dt}\right)=0.\hfill \end{array}$$

(31)

The nonlinear operators $\mathcal{N}\left[f_0\left(t\right)\right]$, ${\mathcal{N}}_{\overline{f}}\left[f_0\left(t\right)\right]$, ${\mathcal{N}}_{{\overline{f}}^{\prime }}\left[f_0\left(t\right)\right]$ and ${\mathcal{N}}_{{\overline{f}}^{″}}\left[f_0\left(t\right)\right]$ which appear in Equations (30) and (31) may be written in the general form

$${\displaystyle \sum _{i=1}^m{h}_i\left(t\right)g_i\left(t\right)K_i\left(t\right),}$$

(32)

where the known functions $h_i\left(t\right)$, $g_i\left(t\right)$ and $K_i\left(t\right)$ depend on the function $f_0\left(t\right)$ and also on the nonlinear operator $\mathcal{N}$, m being a known integer. It is more convenient to consider the unknown function $f_1(t,C_{k_1})$ in the following forms

$${\displaystyle f_1(t,C_{k_1})=\sum _{i=1}^n{\tilde{H}}_{1,i}\left(t,h_j\left(t\right),C_{k_1}\right)g_i\left(t\right)K_i\left(t\right),k_1=1,2,\dots ,s,j=1,2,\dots ,}$$

(33)

or

$${\displaystyle f_1(t,C_{k_1})=\sum _{i=1}^n{\tilde{H}}_{1,i}\left(t,g_j\left(t\right),C_{k_1}\right)h_i\left(t\right)K_i\left(t\right),k_1=1,2,\dots ,s,j=1,2,\dots ,}$$

(34)

$$B\left(f_1(t,C_{k_1}),\frac{d{f}_1(t,C_{k_1})}{dt}\right)=0,$$

(35)

and the unknown function $f_2(t,C_{k_1})$ in the forms

$${\displaystyle f_2(t,C_{k_1})=\sum _{i=1}^n{\tilde{H}}_{2,i}\left(t,h_j\left(t\right),C_{k_1}\right)g_i\left(t\right)K_i\left(t\right),k_1=1,2,\dots ,s,j=1,2,\dots ,}$$

(36)

or

$${\displaystyle f_2(t,C_{k_1})=\sum _{i=1}^n{\tilde{H}}_{2,i}\left(t,g_j\left(t\right),C_{k_1}\right)h_i\left(t\right)K_i\left(t\right),k_1=1,2,\dots ,s,j=1,2,\dots ,}$$

(37)

$$B\left(f_2(t,C_{k_1}),\frac{d{f}_2(t,C_{k_1})}{dt}\right)=0.$$

(38)

The expression of the optimal auxiliary functions ${\tilde{H}}_{j,i}(t,h_j,C_{k_1})$ contains combinations of functions $h_j$ and several unknown parameters $C_{k_1}$, $k_1=1,2,\dots ,s$, and n is an arbitrary integer number. There are a lot a possibility to choose the value of the integer positive n and the optimal auxiliary functions ${\tilde{H}}_{j,i}$. The values of the parameters $C_{k_1}$ can be optimally computed via various methods: the collocation method, the least-square method, the Galerkin method, the Ritz method, and so on.

The second-order semi-analytical solutions given by Equation (21) are well-determined.

4. Approximate Analytic Solutions via OHPM

For the initial conditions $A=u\left(0\right)$ and $B=\dot{u}\left(0\right)$, the problem given by Equation (12) is semi-analytically solved by the OHPM technique.

Let ${\omega }_0>0$ be an unknown parameter. This problem can be rewritten as

$$\mathcal{L}\left[u\right(t\left)\right]+\mathcal{N}\left[u\right(t\left)\right]=0,$$

(39)

where the linear and nonlinear operators $\mathcal{L}$ and $\mathcal{N}$ could be chosen in the following form:

$$\begin{array}{c}\mathcal{L}\left[u\left(t\right)\right]=\ddot{u}\left(t\right)+{\omega }_0^{2}u\left(t\right),\hfill \\ \mathcal{N}\left[u\left(t\right)\right]=\left(-{\omega }_0^2-k{R}^2-(k+1)\alpha ·R\right)·u\left(t\right)+\ddot{u}\left(t\right)u^2\left(t\right)-2u\left(t\right){\left(u^{\prime }\left(t\right)\right)}^2\hfill \\ +\left(k{R}^2-(k+1)\alpha ·R\right)u^3\left(t\right).\hfill \end{array}$$

(40)

By using Equation (40) and the initial conditions given by Equation (12), the initial approximation $u_0\left(t\right)$ is the solution of the problem via Equation (18) ($g\left(t\right)$ = 0), namely

$${\displaystyle u_0\left(t\right)=A·cos{\omega }_{0}t+\frac{B}{{\omega }_0}·sin{\omega }_{0}t.}$$

(41)

Equation (30) becomes

$${\ddot{u}}_1\left(t\right)+{\omega }_0^2{u}_1\left(t\right)+H_{1,0}(t,C_{k_1})·\mathcal{N}\left[u_0\left(t\right)\right]=0,u_1\left(0\right)=0,{\dot{u}}_1\left(0\right)=0,$$

(42)

where the nonlinear operator $\mathcal{N}\left[u_0\left(t\right)\right]$ is

$$\mathcal{N}\left[u_0\left(t\right)\right]$$

$$=\left(-{\omega }_0^2-k{R}^2-(k+1)\alpha ·R\right)·u_0\left(t\right)+{\ddot{u}}_0\left(t\right)u_0^2\left(t\right)-2u_0\left(t\right){\left(u_0^{\prime }\left(t\right)\right)}^2+\left(k{R}^2-(k+1)\alpha ·R\right)u_0^3\left(t\right)$$

$$=(M_0+N_0)\left(A·cos{\omega }_{0}t+B·sin{\omega }_{0}t\right)+P_0·cos\left(3{\omega }_{0}t\right)+Q_0·sin\left(3{\omega }_{0}t\right)),$$

with

$$M_0=\left[-\frac{5}{4}{\omega }_0^2+\frac{3}{4}\left(k{R}^2+(k+1)\alpha ·R\right)\right](A^2+B^2),N_0=-{\omega }_0^2-k{R}^2-(k+1)\alpha ·R,$$

$$P_0=\left[\frac{1}{4}{\omega }_0^2+\frac{1}{4}\left(k{R}^2+(k+1)\alpha ·R\right)\right]A(A^2-3B^2),$$

$$Q_0=\left[\frac{1}{4}{\omega }_0^2+\frac{1}{4}\left(k{R}^2+(k+1)\alpha ·R\right)\right]B(3A^2-B^2).$$

Choosing the auxiliary function $H_{1,0}(t,C_{k_1})=cos\left(6{\omega }_{0}t\right)$, the solution of Equation (42) is:

$$\begin{array}{c}{u}_1\left(t\right)=\frac{1}{480{\omega }_0^2}·\left[(-15M_0-33P_0)cos{\omega }_{0}t+30P_{0}cos\left(3{\omega }_{0}t\right)+10M_{0}cos\left(5{\omega }_{0}t\right)\right.\hfill \\ +5M_{0}cos\left(7{\omega }_{0}t\right)+3P_{0}cos\left(9{\omega }_{0}t\right)+(15N_0+63Q_0)sin{\omega }_{0}t-30Q_{0}sin\left(3{\omega }_{0}t\right)\hfill \\ \left.-10N_{0}sin\left(5{\omega }_{0}t\right)+5N_{0}sin\left(7{\omega }_{0}t\right)+3Q_{0}sin\left(9{\omega }_{0}t\right)\right]\hfill \end{array}$$

(43)

Equation (31) becomes

$$\begin{array}{c}{\ddot{u}}_2\left(t\right)+{\omega }_0^2{u}_2\left(t\right)+H_{2,0}(t,C_{k_1})·u_1·{\mathcal{N}}_{\overline{u}}\left[u_0\left(t\right)\right]+H_{2,1}(t,C_{k_1})·{u^{\prime }}_1·{\mathcal{N}}_{{\overline{u}}^{\prime }}\left[u_0\left(t\right)\right]\hfill \\ +H_{2,2}(t,C_{k_1})·{u^{″}}_1·{\mathcal{N}}_{{\overline{u}}^{″}}\left[u_0\left(t\right)\right]=0,u_2\left(0\right)=0,{\dot{u}}_2\left(0\right)=0,\hfill \end{array}$$

(44)

where the nonlinear operators ${\mathcal{N}}_{\overline{u}}\left[u_0\left(t\right)\right]$, ${\mathcal{N}}_{{\overline{u}}^{\prime }}\left[u_0\left(t\right)\right]$ and ${\mathcal{N}}_{{\overline{u}}^{″}}\left[u_0\left(t\right)\right]$ are

$${\mathcal{N}}_{\overline{u}}\left[u_0\left(t\right)\right]=\left(-{\omega }_0^2-k{R}^2-(k+1)\alpha ·R\right)+2{\ddot{u}}_0\left(t\right)u_0\left(t\right)-2{\left(u_0^{\prime }\left(t\right)\right)}^2+3\left(k{R}^2-(k+1)\alpha ·R\right)u_0^2\left(t\right)$$

$$=M_1+N_{1}cos\left(2{\omega }_{0}t\right)+P_{1}sin\left(2{\omega }_{0}t\right),$$

$${\mathcal{N}}_{{\overline{u}}^{\prime }}\left[u_0\left(t\right)\right]=-4u_0{\dot{u}}_0=-4{\omega }_{0}ABcos\left(2{\omega }_{0}t\right)+2{\omega }_0(A^2-B^2)sin\left(2{\omega }_{0}t\right),$$

$${\mathcal{N}}_{{\overline{u}}^{″}}\left[u_0\left(t\right)\right]=u_0^2=\frac{A^2+B^2}{2}+\frac{A^2-B^2}{2}cos\left(2{\omega }_{0}t\right)+ABsin\left(2{\omega }_{0}t\right),$$

with

$$M_1=-{\omega }_0^2-k{R}^2-(k+1)\alpha ·R-2{\omega }_0^2(A^2+B^2)+\frac{3}{2}(A^2+B^2)\left(k{R}^2-(k+1)\alpha ·R\right),$$

$$N_1=\frac{3}{2}(A^2-B^2)\left(k{R}^2-(k+1)\alpha ·R\right),P_1=3AB\left(k{R}^2-(k+1)\alpha ·R\right).$$

If $n_1\ge 6$ is a fixed integer number, then choosing the auxiliary functions $H_{2,0}(t,C_{k_1})=0$, $H_{2,1}(t,C_{k_1})=0$ and ${\displaystyle H_{2,1}(t,C_{k_1})=\sum _{i=6}^{n_1}{b}_{i}cos\left((2i+1){\omega }_{0}t\right)+c_{i}cos\left((2i+1){\omega }_{0}t\right)}$, respectively, the solution of the Equation (44) has the following form:

$$\begin{array}{c}{\displaystyle u_2\left(t\right)=\sum _{i=6}^{N_{max}}{B}_{i}cos\left((2i+1){\omega }_{0}t\right)+C_{i}cos\left((2i+1){\omega }_{0}t\right),}\hfill \end{array}$$

(45)

where $N_{max}>1$ an arbitrary integer number depending on $n_1$ and the convergence-control parameters $B_i$, $C_i$, so will be optimally identified, such that ${\displaystyle {\left\{B_i\right\}}_{i=1}^{N_{max}}\cup {\left\{C_i\right\}}_{i=1}^{N_{max}}={\left\{C_{k_1}\right\}}_{k_1=1}^s}$.

Finally, the second-order approximate analytic solution is obtained from Equations (21), (41), (43) and (45) in the form

$${\displaystyle \overline{u}\left(t\right)=u_0\left(t\right)+u_1(t,C_{k_1})+u_2(t,C_{k_1}),k_1=1,2,\dots ,s.}$$

(46)

Therefore, by applying the same procedure, we obtain the second-order approximate analytic solution $\overline{u}\left(t\right)$ of nonlinear problem Equation (41) and has the form of Equation (46).

For $\alpha >0$, $k>0$, from Equation (8) are obtained the linear operator ${\displaystyle \mathcal{L}\left(u\left(t\right)\right)=\ddot{u}\left(t\right)+{\omega }_0^{2}u\left(t\right)}$ and the nonlinear operator as ${\displaystyle \mathcal{N}\left(u\left(t\right)\right)=-{\omega }_0^{2}u\left(t\right)-\frac{k}{2}{R}^2·sin(2·u\left(t\right))-(k+1)R\alpha sin\left(u\left(t\right)\right)}$. The nonlinear operators ${\mathcal{N}}_u$, ${\mathcal{N}}_{u^{\prime }}$ and ${\mathcal{N}}_{u^{″}}$, respectively, which appear in Equations (30) and (31), become ${\mathcal{N}}_u=-{\omega }_0^2-k{R}^2·\dot{u}\left(t\right)cos(2·u\left(t\right))-(k+1)R\alpha \dot{u}\left(t\right)cos\left(u\left(t\right)\right)$, ${\mathcal{N}}_{u^{\prime }}=0$ and ${\mathcal{N}}_{u^{″}}=0$.

The initial approximation $u_0\left(t\right)=u\left(0\right)·cos\left({\omega }_{0}t\right)+\frac{\dot{u}\left(0\right)}{{\omega }_0}·sin\left({\omega }_{0}t\right)$ satisfies ${\displaystyle \mathcal{L}\left(u\left(t\right)\right)=0}$, with the initial conditions given by Equation (9).

In Equation (8), we can use the approximate expansions

$$\begin{array}{c}{\displaystyle sin\left(u\right)\simeq \sum _{i=0}^{n_2}{(-1)}^i·\frac{u^{2i+1}}{(2i+1)!},}\hfill \\ {\displaystyle cos\left(u\right)\simeq \sum _{i=0}^{n_2}{(-1)}^i·\frac{u^{2i}}{\left(2i\right)!}.}\hfill \end{array}$$

(47)

where $n_2$ is an arbitrary integer number.

The expression ${\displaystyle \mathcal{N}\left(u_0\left(t\right)\right)}$ contains a combination of the elementary functions $cos\left((2i+1){\omega }_{0}t\right)$, $sin\left((2i+1){\omega }_{0}t\right)$, $i=1,2,\dots$. In Equation (30), choosing the auxiliary function $H_{1,0}(t,C_i)=cos\left(2{\omega }_{0}t\right)$, by a simple computation (neglecting the secular terms), the first approximation $u_1\left(t\right)$ has the form

$$\begin{array}{c}{\displaystyle u_1\left(t\right)=\sum _{k=1}^{N_1}{\tilde{a}}_k^{\left(1\right)}·cos\left((2k+1){\omega }_{0}t\right)+{\tilde{b}}_k^{\left(1\right)}·sin\left((2k+1){\omega }_{0}t\right),}\hfill \end{array}$$

(48)

where the convergence-control parameters ${\omega }_0$, ${\tilde{a}}_k^{\left(1\right)}$, ${\tilde{b}}_k^{\left(1\right)}$ will be optimally identified.

Otherwise, ${\displaystyle {\mathcal{N}}_u\left(u_0\left(t\right)\right)}$ contains a combination of the elementary functions $cos\left(\left(2i\right){\omega }_{0}t\right)$, $sin\left(\left(2i\right){\omega }_{0}t\right)$, $i=1,2,\dots$. In Equation (31), choosing the auxiliary functions $H_{2,0}(t,C_i)=cos\left({\omega }_{0}t\right)$, $H_{2,1}(t,C_i)=0$, $H_{2,2}(t,C_i)=0$, by a simple computation (neglecting the secular terms), the second approximation $u_2\left(t\right)$ has the form

$$\begin{array}{c}{\displaystyle u_2\left(t\right)=\sum _{k=1}^{N_2}{\tilde{a}}_k^{\left(2\right)}·cos\left((2k+1){\omega }_{0}t\right)+{\tilde{b}}_k^{\left(2\right)}·sin\left((2k+1){\omega }_{0}t\right),}\hfill \end{array}$$

(49)

where the convergence-control parameters ${\omega }_0$, ${\tilde{a}}_k^{\left(2\right)}$, ${\tilde{b}}_k^{\left(2\right)}$ will be optimally identified.

Therefore, the second-order analytic approximate solutions of the nonlinear problem in Equations (8) and (9) are the same form as in Equation (46).

In this way, the approximate analytic solutions of the nonlinear problems in Equations (12) and (41) can be constructed via the OHPM method.

5. Numerical Results and Discussions

By taking into account of the second-order semi-analytical solutions given by Equations (46), (48) and (49) with the index $N_{max}\in \{20,25,30,35\}$, the effective solutions are graphically presented below and are given in Appendix A. These solutions highlighted the accuracy of the OHPM method.

The precision of the obtained results is shown in Figure 1 and Figure 2 (for $k=0.5$, $\alpha =0.25$) and Figure 3 and Figure 4 (for $k=-0.5$, $\alpha =0.25$), respectively.

The semi-analytical solutions are in good agreement with the corresponding numerical integration results, computed by means of the fourth-order Runge–Kutta method, as shown in these figures.

On the other hand, the case where $k=1$, $\alpha =0$ is depicted in Figure 5 and Figure 6.

For different values of the parameter $N_{max}$, from Equations (46)–(48), the convergence-control parameters $C_1=a_0+v\left(0\right)$, $C_i=a_{k-1}^{\left(4\right)}$, $B_1=b_0+\frac{\dot{v}\left(0\right)}{{\omega }_0}$, $B_i=b_{k-1}^{\left(4\right)}$, $i=2,3,\dots ,N_{max}$, are optimally computed by the least square method.

The symmetry with respect to the $Oz$-axis for $k\in \mathbf{R}$, $\alpha \ne 0$ and the symmetry with respect to the all coordinate axes for $k=1$, $\alpha =0$ is highlighted in Figure 7, Figure 8 and Figure 9 (i.e., the periodic orbit of the Rikitake-type system).

The convergence-control parameters are presented in Appendix A.

Returning to the exact solutions given by Equations (4) and (5), the approximate residuals are as follows:

$$\left\{\begin{array}{c}{\mathcal{R}}_H\left(t\right)=\frac{1}{2}({\overline{x}}^2+{\overline{y}}^2)+\frac{1+k}{2}{\overline{z}}^2-\frac{1}{2}(x_0^2+y_0^2)-\frac{1+k}{2}{z}_0^2\hfill \\ {\mathcal{R}}_C\left(t\right)=\frac{1}{4}({\overline{y}}^2-{\overline{x}}^2)-\alpha \overline{z}+\frac{1-k}{4}{\overline{z}}^2-\frac{1}{4}(y_0^2-x_0^2)+\alpha z_0-\frac{1-k}{4}{z}_0^2\hfill \end{array}\right.,\mathrm{for}\alpha \ne 0,k\in \mathbb{R},$$

(50)

and

$$\left\{\begin{array}{c}{\mathcal{R}}_H\left(t\right)=\frac{1}{2}({\overline{x}}^2+{\overline{y}}^2)+{\overline{z}}^2-\frac{1}{2}(x_0^2+y_0^2)-z_0^2\hfill \\ {\mathcal{R}}_C\left(t\right)={\overline{x}}^2-{\overline{y}}^2-x_0^2+y_0^2\hfill \end{array}\right.,\mathrm{for}\alpha =0,k=1,$$

(51)

respectively, with $\overline{x}$, $\overline{y}$, $\overline{z}$ the semi-analytical solutions obtained via the OHPM method.

Figure 10 and Figure 11 show the profiles of the residuals ${\mathcal{R}}_H$ and ${\mathcal{R}}_C$, given by Equation (50) in the cases where $k=0.5$, $\alpha =0.25$ for $x_0=0.5$, $y_0=0.5$, $z_0=3.5$ and $N_{max}=30$. Figure 12 and Figure 13 show the profiles of the residuals ${\mathcal{R}}_H$ and ${\mathcal{R}}_C$, given by Equation (51) in the case $k=1$, $\alpha =0$ for $x_0=1.5$, $y_0=1.25$, $z_0=0.5$ and $N_{max}=35$. These figures emphasize the accuracy of the semi-analytical solutions obtained by the OHPM technique by an order of magnitude of $\simeq {10}^{-4}$ of errors. Moreover, the precision of the obtained solutions increases if the index number $N_{max}$ increases. This can be seen in

Table 1. Values of the absolute errors: ${ϵ}_u=|u_{numerical}-{\overline{u}}_{OHPM}|$ for $k=0.5$, $\alpha =0.25$, $x_0=0.5$, $y_0=0.5$, $z_0=3.5$ and different values of the index $N_{max}$; ${\overline{u}}_{OHPM}$ obtained from Equations (46) and (A1)–(A3).

t	${\mathit{N}}_{\mathit{max}}=10$	${\mathit{N}}_{\mathit{max}}=15$	${\mathit{N}}_{\mathit{max}}=20$	${\mathit{N}}_{\mathit{max}}=25$	${\mathit{N}}_{\mathit{max}}=30$
0	1.257327 $\times {10}^{-14}$	1.201336 $\times {10}^{-11}$	3.785027 $\times {10}^{-13}$	7.624456 $\times {10}^{-14}$	6.176531 $\times {10}^{-12}$
1/2	4.468917 $\times {10}^{-4}$	9.189849 $\times {10}^{-5}$	1.045288 $\times {10}^{-5}$	2.676522 $\times {10}^{-6}$	1.844087 $\times {10}^{-8}$
1	1.152456 $\times {10}^{-3}$	5.318678 $\times {10}^{-5}$	1.044831 $\times {10}^{-5}$	1.858868 $\times {10}^{-6}$	2.960618 $\times {10}^{-7}$
3/2	1.308409 $\times {10}^{-4}$	1.062711 $\times {10}^{-4}$	2.531494 $\times {10}^{-5}$	2.738590 $\times {10}^{-6}$	3.229869 $\times {10}^{-7}$
2	1.991649 $\times {10}^{-3}$	1.189236 $\times {10}^{-4}$	1.874768 $\times {10}^{-5}$	2.111766 $\times {10}^{-6}$	2.062485 $\times {10}^{-8}$
5/2	7.220934 $\times {10}^{-4}$	1.019888 $\times {10}^{-4}$	8.060565 $\times {10}^{-6}$	2.758030 $\times {10}^{-6}$	1.833048 $\times {10}^{-7}$
3	1.354596 $\times {10}^{-3}$	6.872820 $\times {10}^{-5}$	1.770506 $\times {10}^{-5}$	2.776405 $\times {10}^{-6}$	2.159764 $\times {10}^{-7}$
7/2	1.740249 $\times {10}^{-4}$	6.252692 $\times {10}^{-5}$	1.092423 $\times {10}^{-5}$	1.025640 $\times {10}^{-7}$	1.510173 $\times {10}^{-7}$
4	1.119755 $\times {10}^{-3}$	4.032325 $\times {10}^{-5}$	1.079516 $\times {10}^{-6}$	1.399745 $\times {10}^{-6}$	5.298769 $\times {10}^{-8}$
9/2	6.124672 $\times {10}^{-4}$	1.630336 $\times {10}^{-5}$	2.412921 $\times {10}^{-6}$	1.487217 $\times {10}^{-6}$	4.939037 $\times {10}^{-8}$
5	1.045304 $\times {10}^{-3}$	3.376512 $\times {10}^{-5}$	3.281965 $\times {10}^{-6}$	3.786340 $\times {10}^{-7}$	1.117548 $\times {10}^{-8}$

. Therefore, in each case, the semi-analytical solutions are approaching the exact solution.

Table 1,

Table 2. Values of the absolute errors: ${ϵ}_u=|u_{numerical}-{\overline{u}}_{OHPM}|$ for $k=-0.5$, $\alpha =0.25$, $x_0=0.5$, $y_0=0.5$, $z_0=3.5$ and different values of the index $N_{max}$; ${\overline{u}}_{OHPM}$ obtained from Equations (46) and (A4)–(A6).

t	${\mathit{N}}_{\mathit{max}}=20$	${\mathit{N}}_{\mathit{max}}=25$	${\mathit{N}}_{\mathit{max}}=30$
0	2.305100 $\times {10}^{-14}$	1.403738 $\times {10}^{-13}$	5.510869 $\times {10}^{-14}$
1/2	5.947881 $\times {10}^{-10}$	7.017700 $\times {10}^{-10}$	4.043440 $\times {10}^{-10}$
1	4.880902 $\times {10}^{-11}$	4.392379 $\times {10}^{-10}$	2.660502 $\times {10}^{-10}$
3/2	3.681659 $\times {10}^{-10}$	7.508715 $\times {10}^{-12}$	3.074690 $\times {10}^{-11}$
2	3.997419 $\times {10}^{-11}$	7.857658 $\times {10}^{-12}$	3.787098 $\times {10}^{-11}$
5/2	6.828394 $\times {10}^{-11}$	1.448007 $\times {10}^{-11}$	5.347465 $\times {10}^{-11}$
3	2.205485 $\times {10}^{-11}$	1.818600 $\times {10}^{-12}$	6.491557 $\times {10}^{-11}$
7/2	6.357189 $\times {10}^{-11}$	8.991882 $\times {10}^{-11}$	4.421671 $\times {10}^{-11}$
4	1.532726 $\times {10}^{-10}$	4.305329 $\times {10}^{-11}$	2.281703 $\times {10}^{-11}$
9/2	9.329380 $\times {10}^{-11}$	1.755498 $\times {10}^{-12}$	8.361002 $\times {10}^{-11}$
5	1.002945 $\times {10}^{-10}$	9.415189 $\times {10}^{-11}$	2.679244 $\times {10}^{-11}$

and

Table 3. Values of the absolute errors: ${ϵ}_u=|u_{numerical}-{\overline{u}}_{OHPM}|$ for $k=1$, $\alpha =0$, $x_0=1.5$, $y_0=1.25$, $z_0=0.5$ and different values of the index $N_{max}$; ${\overline{u}}_{OHPM}$ obtained from Equations (46) and (A7)–(A10).

t	${\mathit{N}}_{\mathit{max}}=20$	${\mathit{N}}_{\mathit{max}}=25$	${\mathit{N}}_{\mathit{max}}=30$	${\mathit{N}}_{\mathit{max}}=35$
0	2.575717 $\times {10}^{-14}$	1.154631 $\times {10}^{-14}$	9.769962 $\times {10}^{-15}$	1.776356 $\times {10}^{-15}$
1/2	8.589703 $\times {10}^{-5}$	2.034201 $\times {10}^{-5}$	3.926340 $\times {10}^{-6}$	7.512735 $\times {10}^{-6}$
1	1.072679 $\times {10}^{-4}$	1.670757 $\times {10}^{-5}$	3.924268 $\times {10}^{-6}$	7.510218 $\times {10}^{-6}$
3/2	8.272699 $\times {10}^{-5}$	2.130846 $\times {10}^{-5}$	4.292848 $\times {10}^{-6}$	7.505662 $\times {10}^{-6}$
2	7.439898 $\times {10}^{-5}$	2.030329 $\times {10}^{-5}$	3.818494 $\times {10}^{-6}$	7.512221 $\times {10}^{-6}$
5/2	7.643807 $\times {10}^{-5}$	1.820285 $\times {10}^{-5}$	4.163315 $\times {10}^{-6}$	7.525591 $\times {10}^{-6}$
3	7.848230 $\times {10}^{-5}$	1.887592 $\times {10}^{-5}$	3.981412 $\times {10}^{-6}$	7.511741 $\times {10}^{-6}$
7/2	7.910670 $\times {10}^{-5}$	2.118086 $\times {10}^{-5}$	4.118728 $\times {10}^{-6}$	7.517661 $\times {10}^{-6}$
4	8.160504 $\times {10}^{-5}$	2.078708 $\times {10}^{-5}$	4.070157 $\times {10}^{-6}$	7.532979 $\times {10}^{-6}$
9/2	8.540867 $\times {10}^{-5}$	1.944186 $\times {10}^{-5}$	4.058795 $\times {10}^{-6}$	7.515431 $\times {10}^{-6}$
5	8.781798 $\times {10}^{-5}$	2.054146 $\times {10}^{-5}$	4.001397 $\times {10}^{-6}$	7.502956 $\times {10}^{-6}$

show the influence of the index number $N_{max}$ on the values of the absolute errors. The values of the better semi-analytical solutions are shown in

Table 4. The semi-analytical solution ${\overline{x}}_{OHPM}$ given by Equations (7), (46) and (A3) by comparison with corresponding numerical solution for $k=0.5$, $\alpha =0.25$, $x_0=0.5$, $y_0=0.5$, $z_0=3.5$ and the index $N_{max}=30$; (absolute errors: ${ϵ}_x=|x_{numerical}-{\overline{x}}_{OHPM}|$).

t	${\mathit{x}}_{\mathit{numerical}}$	${\overline{\mathit{x}}}_{\mathit{OHPM}}$	${\mathit{ϵ}}_{\mathit{x}}$
0	0.5	0.5000000256	2.568413 $\times {10}^{-8}$
1/2	1.5214552317	1.5214806203	2.538851 $\times {10}^{-5}$
1	2.8716257775	2.8716245167	1.260769 $\times {10}^{-6}$
3/2	2.2488093942	2.2488113818	1.987665 $\times {10}^{-6}$
2	2.8435419487	2.8435503176	8.368986 $\times {10}^{-6}$
5/2	1.7711339212	1.7711378969	3.975644 $\times {10}^{-6}$
3	0.5827320936	0.5827293699	2.723671 $\times {10}^{-6}$
7/2	0.3106026991	0.3106006161	2.082953 $\times {10}^{-6}$
4	0.6229399329	0.6229416108	1.677919 $\times {10}^{-6}$
9/2	1.8801975655	1.8801969832	5.823633 $\times {10}^{-7}$
5	2.8117250468	2.8117289431	3.896274 $\times {10}^{-6}$

,

Table 5. The semi-analytical solution ${\overline{y}}_{OHPM}$ given by Equations (10), (46) and (A6) by comparison with corresponding numerical solution for $k=-0.5$, $\alpha =0.25$, $x_0=0.5$, $y_0=0.5$, $z_0=3.5$ and the index $N_{max}=30$; (absolute errors: ${ϵ}_y=|y_{numerical}-{\overline{y}}_{OHPM}|$).

t	${\mathit{y}}_{\mathit{numerical}}$	${\overline{\mathit{y}}}_{\mathit{OHPM}}$	${\mathit{ϵ}}_{\mathit{y}}$
0	0.5	0.4999999999	3.561040 $\times {10}^{-13}$
1/2	0.8847309585	0.8847308956	6.291253 $\times {10}^{-8}$
1	0.3219523370	0.3219522986	3.837036 $\times {10}^{-8}$
3/2	−0.6010213841	−0.6010213748	9.341504 $\times {10}^{-9}$
2	−0.8621287789	−0.8621288420	6.308938 $\times {10}^{-8}$
5/2	−0.1974129935	−0.1974130787	8.513892 $\times {10}^{-8}$
3	0.6890139943	0.6890140472	5.291821 $\times {10}^{-8}$
7/2	0.8221568424	0.8221570945	2.520653 $\times {10}^{-7}$
4	0.0683396173	0.0683395592	5.812823 $\times {10}^{-8}$
9/2	−0.7623835272	−0.7623835803	5.315807 $\times {10}^{-8}$
5	−0.7653990825	−0.7653994132	3.306396 $\times {10}^{-7}$

and

Table 6. The semi-analytical solution ${\overline{z}}_{OHPM}$ given by Equations (39), (46), (A9) by comparison with corresponding numerical solution for $k=1$, $\alpha =0$, $x_0=1.5$, $y_0=1.25$, $z_0=0.5$ and the index $N_{max}=35$; (absolute errors: ${ϵ}_z=|z_{numerical}-{\overline{z}}_{OHPM}|$).

t	${\mathit{z}}_{\mathit{numerical}}$	${\overline{\mathit{z}}}_{\mathit{OHPM}}$	${\mathit{ϵ}}_{\mathit{z}}$
0	0.5	0.4999999999	3.275157 $\times {10}^{-15}$
3/5	−0.6939047661	−0.6938988897	5.876456 $\times {10}^{-6}$
6/5	−1.2734211748	−1.2734197882	1.386637 $\times {10}^{-6}$
9/5	−1.3237111667	−1.3237099418	1.224830 $\times {10}^{-6}$
12/5	−0.9302874402	−0.9302835499	3.890361 $\times {10}^{-6}$
3	0.1717944434	0.1718080341	1.359061 $\times {10}^{-5}$
18/5	1.1016264922	1.1016412743	1.478212 $\times {10}^{-5}$
21/5	1.3448426458	1.3448541193	1.147350 $\times {10}^{-5}$
24/5	1.1874060275	1.1874195726	1.354508 $\times {10}^{-5}$
27/5	0.3999481041	0.3999621828	1.407868 $\times {10}^{-5}$
6	−0.7760295495	−0.7760248831	4.666312 $\times {10}^{-6}$

.

6. Conclusions

In the present work, the Hamilton–Poisson realization of the Rikitake-type system is emphasized and the semi-analytical solutions is established in a closed form. This 3D differential system describes the motion of the two disk dynamos coupled to one another in relation to Earth’s magnetic field. Only the periodic behaviors are taken into account, excluding chaotic ones.

The influence of the global analytic first integrals (the Hamiltonian and the Casimir functions, respectively) of the Rikitake-type system leads to a simplification to a second-order nonlinear differential equation for which the semi-analytical solutions are established.

The solution method is the optimal homotopy perturbation method (OHPM). Only two iterations are used. This method involves the presence of some auxiliary functions which depend on a set of convergence-control parameters. These parameters are optimally computed for some values of the physical parameters k and $\alpha$. The precision of the obtained solutions is shown by comparison with the corresponding numerical results (via the fourth-order Runge–Kutta method). This comparison highlights the effectiveness of the applied method.

The obtain semi-analytical solutions could be successfully used for the study of dynamical systems encountered in electrical engineering (the synchronization, secure communications or optimization of nonlinear system performance) or medicine.