Expedite homotopy perturbation method based on metaheuristic technique mimicked by the flashing behavior of fireflies

1.   Introduction

In the past few years, fractal calculus and fractional calculus has been one of the most rapidly growing areas of mathematical analysis, which are used to model various problems of real life. Fractal calculus is considered as a fruitful field of research in science and technology and is extensively used to elucidate the phenomena of hierarchical or porous media ^[1]. Owing to the simplicity and effectiveness of fractal derivative it has defined an alternative approach of fractional derivatives, to explain some of the most fundamental theories with considerable ease and elegance ^[2,3]. The fractional derivative and fractional integral operators are the valuable tools which are used in the modeling of various physical phenomena of science and engineering. Anomalous dynamics of numerous complex nonlinear systems are elegantly modeled by the assistance of fractional differential equations, which are acknowledged as the generalization of the classical differential equations of integer order ^[4,5,6].

The fractional order Ricatti equation (RE) accommodates an extremely significant class of nonlinear differential equations which includes: The nth order linear homogenous ODE, one-dimensional Schrödinger equation, wave solution of nonlinear partial differential equation (PDE), etc ^[7]. These REs has substantial importance in classical, as well as, modern science and engineering problems because of their multifarious applications in various fields. For instance, optimal control, variational calculus, random processes, quantum mechanics, thermodynamics, robust stabilization, stochastic realization theory and diffusion problems ^[8,9,10,11].

The mathematical formulation for the nonlinear Ricatti differential equations with the fractional derivative defined in the Caputo sense is read as,

with the initial condition defined as,

where, $\alpha$ is the order of equation, $\alpha > 0$ and $\alpha \in \Re$, with $T, {t_0}, \eta \in \Re$; $p\left(t \right), q\left(t \right), r\left(t \right)$ are the real continuous functions and $x\left(t \right)$is the solution of the equation. The behavior of the aforementioned Eq. (1) depends on the parameter $\alpha$, which can be varied to analyze the dynamical behavior of these equations. In order to study the nonlinear Eq. (1) various analytical and numerical techniques developed in the past were applied by many researchers and scientists to compute the solutions of these equations among which some of them are cited here: The homotopy perturbation method (HPM) ^[12], the fractional-order Legendre operational matrix method ^[13], the iterative reproducing kernel Hilbert space method ^[14], the optimal homotopy asymptotic method ^[15], the modified Laplace Adomian decomposition method ^[16], the B-spline operational method ^[17] and the Haar wavelet collocation method ^[18]. The HPM introduced by He has been successfully used to solve numerous linear and nonlinear problems of real life. Since, its development it has been modified by He and many other scientists to solve various ordinary and partial differential equations of integer and non-integer order ^[19,20,21].

The demand of global optimization technique is increasing day by day, which are utilized in the assessment of numerous nonlinear and multimodal problems of real life. The deterministic algorithm and the stochastic algorithm are the two types of optimization algorithm that are found in the literature ^[22,23]. The deterministic algorithm are often gradient-based whereas, the stochastic algorithm are further subcategorized as heuristic or metaheuristic algorithm. In recent trends, the popularity and demand of the nature inspired algorithms have increased extensively. Nature inspired metaheuristic algorithms which efficiently deals with the nonlinear optimization problems includes the genetic algorithm (GA) ^[24], simulated annealing (SA) algorithm ^[25], differential evolution (DE) algorithm ^[26], the ant colony optimization (ACO) algorithm ^[27], particle swarm optimization (PSO) algorithms ^[28], the shark smell algorithm ^[29] and the most powerful and demanding firefly (FA) algorithm ^[30,31,32].

The FA which mimics the flashing pattern and behavior of the fireflies was first developed by Yang in 2008 and since then it has been modified by various scientist and researcher to solve different types of challenging optimization problems. Some of the flashing characteristic of these unisex fireflies are idealized to develop the FA. The three idealized rule used by the firefly algorithm are as follows:

● A firefly is attracted to the other fireflies regardless of their gender.

● Attractiveness of fireflies is directly proportional to their brightness. The attractiveness and brightness of fireflies both increases as the distance between them deceases.

● Brightness of a firefly is determined by the objective function.

The literature of this modern, self-adaptive, highly efficient and truly intelligent algorithm has expanded dramatically ^{[30,31,32,33,34]}.

The fundamental aim of this paper is to provide a numerical technique for the assessment of nonlinear fractional differential model defined as (1) and (2). The concept of classical homotopy perturbation method is merged with the modern metaheuristic optimization technique for the development of an expedite homptopy perturbation method (EHPM). The developed method-EHPM transforms the fractional model into a system of algebraic equations leading to a fitness function determined by a particular fragment of the weighted series solution, which are trained by the using the powerful and reliable optimization technique-FA. The optimal values achieved by the FA are utilized to attain the accurate, convergent and reliable solutions. Comparative study is conducted by the comparing the EHPM computed results with the available exact solution and the solutions obtained by the modified homotopy perturbation (MHPM) ^[35], the residual power series method (RPSM) ^[36] and the Adam bashforth method (ABMA) ^[37]. Furthermore, the accuracy and competency of the EHPM is also ratified by computing results by the proposed design methodology in combination with accelerated particle swarm optimization (APSO), i.e., the fitness function determined by EHPM is optimized by using APSO. Various error measures are also carried out to validate the correctness and accuracy of the suggested scheme.

2.   Elementary tools

Some basic definitions of the fractional calculus which are going to be utilized in the further discussion are stated below:

Definition 2.1. Let $f\left(x \right)$ be a differential function with $\beta \in \left({0, 1} \right],$then the Caputo order fractional derivative ${}^CD_0^\beta f\left(x \right),$ is defined as ^[36],

with ${}^CD_0^\beta \omega = 0,$ for some constant $\omega \in {\mathfrak{R}}$.

Definition 2.2. For any function $f\left(x \right)$ with $x:\left[{0, \infty } \right) \to \Re,$ the Riemann fractional integral of order $\beta$ is given as ^[35],

where $\beta > 0$ and $\beta \in \Re.$

3.   Analysis of the method

This section comprehensively describes the procedure to determine the approximate solutions of the nonlinear fractional differential equations using the classical idea of HPM with the modern optimizing tool. A brief review of the learning solver FA, which is utilized for the development of the presented scheme (EHPM) is also demonstrated here.

3.1.   The expedite homotopy perturbation method (EHPM)

To exemplify the basic ideas of the proposed scheme consider the nonlinear FDE of the form,

with the initial condition defined as

where $g$ and $x$are function of $t$ and $\Omega$ is the boundary of the domain. Then $\rm{ \mathit{ ψ}}$ can be expressed as

Where L is the linear operator, N is the nonlinear operator and A is the known analytical function. By the homotopy procedure we rewrite Eq. (7) as

where ${{\mathit{\boldsymbol{x}}_0}\left(t \right)}$ is an initial approximation, and ${\bf{ {{ δ}} }}$ is an embedding parameter defined in some closed interval $\left[{0, 1} \right]$. Let the homotopy solution of equation Eq. (8) can be written as

where ${{\mathit{\boldsymbol{X}}_k}\left(t \right)}$ are attained by the $kth$-order homotopy derivative defined as

Now, we let the initial approximation of Eq. (8) of the form

where ${\xi _0}, {\xi _1}, {\xi _2}, \ldots, {\xi _n}$ are the unknown weights to be determined. By utilizing the Eq. (11), Eq. (9) and the basic ideas of MHPM ^[35] with ${\mathit{\boldsymbol{X}}_1}\left(t \right) = 0$ and ${\mathit{\boldsymbol{L}}^{\mathit{\boldsymbol{ - 1}}}}$ as the inverse operator of $L$ we obtain an algebraic system expressed as

and ${\mathit{\boldsymbol{X}}_2}\left(t \right) = {\mathit{\boldsymbol{X}}_3}\left(t \right) = \ldots = 0$. Setting ${\mathit{\boldsymbol{X}}_1}\left(t \right) = 0$ by MHPM, automatically leads the solution of Eq. (8) written as:

The aforementioned approximate solution (14) comprises of the unknown weights ${\xi _0}, {\xi _1}, {\xi _2}, \ldots, {\xi _n}$ which are determined by constructing a fitness function given as

The weights in Eq. (15) are learned by using a modern and powerful optimization technique FA. The graphical abstract of the above presented scheme EHPM is portrayed in Figure 1.

Lemma 3.1. In problem (5) if $x\left(t \right) = {\mathit{\boldsymbol{X}}_\mathit{\boldsymbol{0}}}$ is the solutions of Eq. (9), which comprises of the unknown weights ${\xi _i}\begin{array}{*{20}{c}}; \end{array}i = 0, 1, 2, ..., n,$ then the minimization problem will be simplified as:

Theorem 3.2. If the solution of Eq. (8) satisfies ${\mathit{\boldsymbol{X}}_1}\left(t \right) = 0,$ then Eq. (10) results ${\mathit{\boldsymbol{X}}_2}\left(t \right) = {\mathit{\boldsymbol{X}}_3}\left(t \right) = \ldots = 0$ and $x\left(t \right) = {\mathit{\boldsymbol{X}}_0}$ as the solution of Eq. (8).

Here we mention that if $g\left({r\left(t \right)} \right)$ and ${\mathit{\boldsymbol{X}}_0}$ are analytic at $t = {t_0}$, then their power series is defined as

where ${\xi _0}^*, {\xi _1}^*, {\xi _2}^*, \ldots$are known coefficients and ${\xi _0}, {\xi _1}, {\xi _2}, \ldots$are unknown ones which are to be determined. In order to obtain ${\xi _0}, {\xi _1}, {\xi _2}, \ldots$ we adapt a nature inspired strategy (FA) to find the approximate solution, by using the above Lemma 3.1. for constructing the fitness function formulated as:

Theorem 3.3. Let $x\left(t \right)$ be an integrable and continuous function defined in some domain $[a, T]$ where $a \in \Re,$ ${t^{i\alpha }}$ is bounded and continuous function defined in the same interval $[a, T]$ for some, $\varepsilon \in {\Re ^ + }$ and real bounded weights, ${\xi _i} \in \Re$ then by substituting the real bounded weights attained by the FA, the obtained series solution converges as $n$ approaches to infinity.

By simple calculations and assuming $\left| {{\xi _i}} \right| \leqslant \hat M$; $\hat M \in {\Re ^ + }$ then Eq. (14) can be written as:

As the number of terms in the initial approximation, $n \to \infty$ it leads$\left| {\Gamma \left({\left({n\alpha + 1} \right) + 1} \right)} \right| > 1$.

3.2.   Learning solver: Firefly algorithm

The FA is based on the flashing pattern and behavior of the fireflies. These flies are unisex, so they are attracted to the other fireflies regardless of their gender. The attractiveness of these tropical fireflies is proportional to their brightness. The main purpose of this attraction is to enable an algorithm to converge quickly, by allowing the swarming agents to interact and move towards the true global optimality. In FA the attractiveness between the fireflies at ${{\bf{ β}}_i}$ and ${{\bf{ β}}_j}$ are determined as

where the middle term is due to the attraction $\theta \left({{{\bf{ β}}_i}, {{\bf{ β}}_j}} \right)$ of the fireflies, which varies with the distance ${r_{ij}}$ between them and can be modeled as:

where

The attractiveness at distance zero is represented by ${\theta _0}$, $\tilde \alpha \in \left[{0, 1} \right]$ is a randomization parameter and $\mu$ is a random number generator that is uniformly distributed in the interval $\left[{0, 1} \right]$. The light absorption coefficient is denoted by $\gamma \in \left({0, \infty } \right]$, which is assumed to be very large but in practice it is determined by the characteristic distance $\Gamma$ of the system, over which the attractiveness varies from ${\theta _0}$ to ${\theta _0}{e^{ - 1}}$. The convergence of the algorithm can be further improved by varying the randomization parameter $\tilde \alpha$ that is decreased gradually as the optima is approached.

The nonlinear updating equation used by the FA yields a richer behavior and higher convergence than the other optimization algorithms with linear updating equating. The middle term of the updating equation becomes negligible by setting a very large value of $\gamma$, leading to the standard SA algorithm. For $\tilde \alpha = 0$ and $\gamma \to 0$ the exponential term in Eq. (21) tends to one, leading Eq. (20) to be a variant of the DE algorithm. Moreover, the APSO algorithm and the harmonic search (HS) algorithm are also the special cases of the FA. Therefore, we can essentially say that the FA is an amalgamation of all these four algorithm (SA, DE, APSO, HS) to a certain extent and can easily outperform other modern and powerful optimization algorithms ^[33,34].

4.   Design methodology

The step by step procedure of the proposed algorithm for Eq. (1) is as follows:

Step1: Construct homotopy $\mathit{\boldsymbol{X}}\left({r\left(t \right), {\bf{ {{ δ}} }}} \right):\Omega \times \left[{0, T} \right] \to \Re$ by using Eqs. (9) and (11).

Step2: Fix the number of terms in the initial approximation ${\mathit{\boldsymbol{x}}_0}\left(t \right)$.

Step3: Set ${\mathit{\boldsymbol{X}}_1}\left(t \right) = 0$ and the equidistant points for in the interval

Step4:Configure the error function $E\left({{\xi _i}} \right)$ specified in Eq. (15).

Step5: Utilize the FA with $\gamma \to \infty$, to optimize the derived error function and attain the minimum fitness value for an effective value of each unknown.

Step6: Substitute the optimal values of ${\xi _i}$ in Eq. (14) to acquire the approximate solution.

Step7: Calculate error norms to validate to correctness and efficiency of the suggested scheme.

5.   Problem simulation and discussion

In this section, the accuracy and competency of the proposed scheme is illustrated by considering several examples of the form (1). Comparison is made with the available exact solution and the results obtained by some former techniques such as, ABMA, APSO, MHPM and RPSM. Also, the accuracy of the algorithm are assessed in terms of the error norms mathematically formulated as

where $\hat n$ stands for the number of input grid points in the computational domain $\left[{{t_0}, T} \right]$, for $T \in {\Re ^ + }$ and $x_j^{rf}$ and $x_j^*$ represents the reference solution and the approximate solution, respectively.

Test problem 1

Consider a nonlinear fractional differential equation

with initial condition $x\left(0 \right) = 0.5$. The exact solution of above problem for $\alpha = 1$ is given as

Test problem 2

Consider a nonlinear factional Riccati differential equation with the fractional derivative defined in the Caputo sense expressed as

and the initial condition defined as $x\left(0 \right) = 0$. The exact solution of the above problem for $\alpha = 1$ is given as

Test problem 3

Now we consider another fractional differential equation

with initial condition $x\left(0 \right) = 0$ and the exact solution of Eq. (27) for $\alpha = 1$ is

Test problem 4

Consider another fractional differential equation

with initial condition $x\left(0 \right) = 0$. By the use of definition (2.1) the exact solution of Eq. (32) is given as

The proposed methodology EHPM accompanied by an efficient and powerful optimization technique is applied to the nonlinear problems defined above. The fitness function derived in the form (18) for the nonlinear problems (26), (28), (30) and (32) are optimized using the FA. The optimal weights achieved by the implementation of the developed technique at $\alpha = 1$, $n = 5$ and $x \in \left[{0, 1} \right]$ yielding a fitness value $3.177187 \times {10^{ - 17}}$, $4.316209 \times {10^{ - 12}}$, and $1.43145 \times {10^{ - 17}}$ for the test problems 1–4 are displayed in Figure 2, respectively. The numerical results accomplished by the utilization of these optimal weights for the respective test problems are presented in Tables 1–4. Comparative study is conducted by comparing the EHPM computed results with the respective exact solution and the results attained by ABMA with the step size $\bar h = 0.001$ and $\bar K = 1000$, $10$ term approximate solution by MHPM, the APSO computed results (for 800 iterations and swarm size equal to $40$) and the seventh term approximate solutions obtained by the RPSM. To illustrate the reliability and accuracy of the suggested scheme EHPM the values of absolute error attained by the considering the exact solution as the reference solution of the respective problems 1–4 are also presented in Tables 1–4, respectively. It is observed that generally, the EHPM computed results for all the problems coincides the exact solution within five to ten decimal places of accuracy. The error norms ${L_{rms}}$ and ${L_\infty }$ acquired for the above problems by consuming the developed scheme EHPM for the computational domain $x \in \left[{0, 1} \right]$ and at distinct values $n$ are tabulated in Table 5. These error estimates presented in Table 5 ensures about the convergence and accuracy of the deliberated technique which is seen to be increased as $n$ increases.

Numerical solutions of the above problems are also constructed by the discussed technique for the non-integer order. The solutions of the above nonlinear problems are attained by the learning of unknown weights in the derived fitness function at two distinct values of $\alpha$ i.e., $0.95$ and $0.85.$ The optimal weights achieved by the proposed scheme for the considered fractional values with $n = 2$ and $x \in \left[{0, 1} \right]$ for problem 1, with $n = 2$ and $x \in \left[{0, 4} \right]$ for problem 2, with $n = 5$ and $x \in \left[{0, 5} \right]$ for problem 3 and with $n = 5$ and $x \in \left[{0, 1} \right]$ for problem 4 are showcased in Figure 3, respectively. The results obtained by utilizing these real valued and bounded optimal weights for the respective problems 1–4 are presented in Tables 6–9, respectively. The results obtained by EHPM are compared with the available exact solution, MHPM computed results, seventh term approximate solutions obtained by the RPSM and the results attained by ABMA for the step $\bar h = 0.001$ and the values $\bar K$ equal to $4000$ and $5000$ for the test problems 1–3, respectively. In problem 1 by taking smaller value $n$ for a smaller domain $\left[{0, 1} \right]$ the results achieved by the EHPM and the RPSM shows a constructive agreement with the results obtained by ABMA. While, by taking a smaller value $n$ for a larger domain $\left[{0, 4} \right]$ in problem 2, the RPSM is seen to diverge whereas, the solutions obtained by the EHPM remains convergent. By taking a larger value $n$ for a larger domain $\left[{0, 5} \right]$ in problem 3 the RPSM computed results diverges whereas, the results accomplished by the proposed scheme shows a constructive agreement with the results attained by ABMA. The values of absolute error achieved by EHPM for the problems 1–3 with the results attained by ABMA as the reference solution are also depicted in Tables 6–8. In problem 4 by taking a larger value of $n$ for a smaller domain $\left[{0, 1} \right]$ the EHPM computed results are compared with available exact solution and the results obtained by MHPM, with the values of absolute error tabulated in Table 9. Furthermore, the correctness and reliability of the suggested scheme is validated by the minimize fitness function and the related graphical solution attained at distinct values of $\alpha$ $=$ 1, 0.95, 0.85 for a lager span, that are portrayed in Figures 4a–7a and Figures 4b–7b for the problems 1–4, respectively. One can infer that the accuracy of the presented scheme can be further enhanced but at the cost of more computation.

6.   Conclusion

In this work, the nonlinear fractional differential equations with the fractional derivative and integral operators defined in the Caputo sense are successfully solved by a simple, accurate and reliable numerical technique. The designed methodology EHPM is an amalgamation of the classical homotopy perturbation technique with the modern bio-inspired metaheuristic technique. The FA which is based on the flashing behavior of fireflies fast track the procedure to determine the approximate solution of the considered nonlinear FDE. Some instructive examples were undertaken, with varying span and different number of terms in the assumed initial approximation to expound the potential ability of the proposed technique. The validity, applicability and computational efficiency of the EHPM is exposed by comparing the fast track outcomes achieved with the available exact solution and the solutions obtained by ABMA, APSO, MHPM and the RPSM. High accuracy and consistent convergence were ascertained by the accomplishment of the optimal values of the error norms. Thus, the contribution of the discussed scheme is abridged as follows:

● The expedite design methodology is well suited to solve the nonlinear FDE.

● The classical approach transformed the nonlinear FDE into an algebraic system which leads to a weighted series solution.

● Fitness function determined by a trivial fragment of the weighted series solution was optimized using the FA.

● Fast track outcomes were achieved by the modern optimization technique (FA).

● Numerical experiments ratified the correctness and reliability of the EHPM.

● Converging solutions were attained even for a larger computational domain.

● High accuracy and consistent convergence achieved by EHPM can be further enhanced by increasing the number of term in the initial approximation.

As the proposed design methodology exhibits a merger of classical idea with the modern optimizing metaheuristic technique, which offers promising converging results. One may utilize the recommended fusion to solve numerous other mathematical models or may even attempt to determine a new combination of modern optimization technique with the former traditional schemes.

Conflict of interest

The authors declare no conflict of interest.