A Comparative Study of the Fractional-Order Clock Chemical Model

Abstract

In this paper, a comparative study has been made between different algorithms to find the numerical solutions of the fractional-order clock chemical model (FOCCM). The spectral collocation method (SCM) with the shifted Legendre polynomials, the two-stage fractional Runge–Kutta method (TSFRK) and the four-stage fractional Runge–Kutta method (FSFRK) are used to approximate the numerical solutions of FOCCM. Our results are compared with the results obtained for the numerical solutions that are based upon the fundamental theorem of fractional calculus as well as the Lagrange polynomial interpolation (LPI). Firstly, the accuracy of the results is checked by computing the absolute error between the numerical solutions by using SCM, TSFRK, FSFRK, and LPI and the exact solution in the case of the fractional-order logistic equation (FOLE). The numerical results demonstrate the accuracy of the proposed method. It is observed that the FSFRK is better than those by SCM, TSFRK and LPI in the case of an integer order. However, the non-integer orders in the cases of the SCM and LPI are better than those obtained by using the TSFRK and FSFRK. Secondly, the absolute error between the numerical solutions of FOCCM based upon SCM, TSFFRK, FSFRK, and LPI for integer order and non-integer order has been computed. The absolute error in the case of the integer order by using the three methods of the third order is considered. For the non-integer order, the order of the absolute error in the case of SCM is found to be the best. Finally, these results are graphically illustrated by means of different figures.

1. Introduction, Definitions and Preliminaries

During the past several decades, the classical calculus has been generalized into fractional calculus. The generalization was the result of a simple question that Leibniz asked in his letter to L’Hôpital about the derivative of order $\frac{1}{2}$. In recent years, many researchers have caught the attention of modeling real-world phenomena by using fractional-order derivatives. The dynamics of many applied problems has been modeled and studied using the concept of fractional-order derivatives. Such problems appear in, for example, biology and physics, ecology, engineering, and various other fields of applied sciences, and so on (see, for details, [1,2,3]; also see the recently-published survey-cum-expository review article [4]).

The main goal of our study in this work is to find the numerical solutions of the fractional-order clock chemical model, which is, we have presented a numerical study. In fact, most of the problems related to the real world are not solved exactly. Hence, in order to verify an appropriate numerical method to find the approximate solution of the fractional-order clock chemical model, a comparative study is conducted between a number of applicable numerical methods. This is indeed useful for many researchers who may use these methods in order to solve other systems and models. Moreover, the fractional-order clock chemical model provides an application of a real-world problem and we need to contribute toward finding its numerical solution, as followed by many researchers in this field. To the best our knowledge, this fractional model has not been studied numerically so far, except possibly in one of our recent papers (see [3]) and, of course, in this work. This justifies the importance of our contribution.

One of the features of mathematical modeling of various biological phenomena is that the mathematical model is represented as a mathematical function in time and the involved parameters. Hence, in this case, we can find the numerical solution to the model and also the parameters that affect this model can be appropriately controlled (see [5,6,7,8,9,10,11,12,13,14,15,16,17,18]).

In this paper, we first contribute toward studying the effect of changes in the order of the fractional derivatives involved. We verify the accuracy of the approximate solution in the presence of a complete solution to the fractional logistic growth model. We then employ different numerical methods in finding approximate solutions for the fractional-order clock chemical model. We also numerically analyze and graphically illustrate these results by means of different figures.

In this section, we give some basic definitions and properties of fractional calculus theory (see, for example [1,2]).

Definition 1.

For $\nu >0,$ the Liouville–Caputo fractional derivative of order $\nu ,$ denoted by ${}_{a_{+}}^{\mathrm{LC}}{\mathcal{D}}_{\zeta }^{\nu },$ is defined by

$$\begin{array}{c}\hfill {}_{a_{+}}^{\mathrm{LC}}{\mathcal{D}}_{\zeta }^{\nu }\alpha \left(\zeta \right)=\frac{1}{\mathsf{\Gamma }(n-\nu )}{\int }_{a_{+}}^{\zeta }{(\zeta -\eta )}^{n-\nu -1}{\mathcal{D}}^n\alpha \left(\eta \right)d\eta ,\end{array}$$

(1)

$$(n-1<\nu <n;n\in \mathbb{N}=\{1,2,3,\cdots \left\}\right),$$

$${}_{a_{+}}^{\mathrm{LC}}{\mathcal{D}}_{\zeta }^{\nu }=D^{\nu }(\nu \in \mathbb{N}),$$

where $D=\frac{d}{d\zeta }.$

The problem under study here is the chemical clock reaction model which is given by [3]. We now replace the ordinary derivative in this model by a fractional-order derivative in the Liouville-Caputo sense (see Definition 1 above), which we have chosen to use here for simplicity instead of Definition 1. The fractional-order chemical clock reaction model is then given by

$$D^{\nu }{\alpha }_1\left(\zeta \right)=\mu e^{-\epsilon \delta \zeta }-{\alpha }_1\left(\zeta \right){\alpha }_2\left(\zeta \right),$$

(2)

$$D^{\nu }{\alpha }_2\left(\zeta \right)={\alpha }_1\left(\zeta \right){\alpha }_2\left(\zeta \right)-\frac{{\alpha }_2\left(\zeta \right)}{\delta }\left({\alpha }_1\left(\zeta \right)+{\alpha }_2\left(\zeta \right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta \zeta }\right)+\lambda -1\right),$$

(3)

$${\alpha }_3\left(\zeta \right)={\alpha }_1\left(\zeta \right)+{\alpha }_2\left(\zeta \right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta \zeta }\right)+\lambda -1$$

(4)

and

$${\alpha }_1\left(0\right)=0,{\alpha }_2\left(0\right)=1\mathrm{and}{\alpha }_3\left(0\right)=\lambda ,\zeta \ge 0.$$

(5)

For more details about this model, one may see its mathematical formulation of in [19,20] to know the definitions for the chemical concentrations ${\alpha }_1\left(\zeta \right),{\alpha }_2\left(\zeta \right)$ and ${\alpha }_3\left(\zeta \right)$ and the parameters $\mu ,$$ϵ,$$\delta$, and $\lambda$.

The main structure of this paper is as follows. In Section 2, the fractional-order Logistic equation and is exact solution are presented. In Section 3, the implementation the spectral collocation and simulation are given. In Section 4, the generalized Runge–Kutta methods are proposed. In Section 5, the numerical results and a systematic discussion are introduced. Finally, in Section 6, the conclusion is presented.

2. The Fractional-Order Logistic Equation

In this section, we introduce the fractional-order Logistic growth model and their exact solution for $0<\nu \leqq 1.$ The dual aspects of complexity, fractional growth and nonlinearly-induced saturation emerge into the fractional-order logistic equation given by

$$D^{\nu }\alpha \left(\zeta \right)=K^{\nu }\alpha \left(\zeta \right)\left(1-\alpha \left(\zeta \right)\right)\left(\alpha \left(0\right)={\alpha }_0\right).$$

(6)

West [21] (see also [22]) constructed an exact solution to the FOLE by using the Carleman embedding technique. The exact solution of (6) is given by

$$\alpha \left(\zeta \right)=\sum _{n=0}^{\infty }{\left(\frac{{\alpha }_0-1}{{\alpha }_0}\right)}^n{E}_{\nu }(-n{K}^{\nu }{\zeta }^{\nu }),$$

(7)

where

$$E_{\nu }\left(\zeta \right)=\sum _{n=0}^{\infty }\frac{{\zeta }^n}{\mathsf{\Gamma }(n\nu +1)}$$

is the Mittag–Leffler function.

3. The Spectral Collocation Method (SCM)

In this section, we apply the SCM for evaluating the numerical solutions of the (2)–(4). In many scientific applications, orthogonal functions and spectral methods have been used to find approximate solutions for fractional differential equations. Using these methods, we can transform the fractional differential equations into a set of algebraic equations. Subsequently, the algebraic equations are solved by one or the other of several known numerical methods (see [23,24,25]).

We begin by defining the shifted Legendre polynomials on the interval $[0,1]$ with the variable $z=2\zeta -1$. These polynomials have the following property:

$$P_s\left(\zeta \right)=P_s(2\zeta -1)=P_{2s}\left(\sqrt{\zeta }\right),$$

where the set $\{P_s\left(z\right):s=0,1,2,\cdots \}$ forms a family of orthogonal Legendre polynomials on the interval $[-1,1]$ (see, for details [26]).

The analytic form of the shifted Legendre polynomials of degree s is given by

$$\tilde{P}\left(\zeta \right)=\sum _{k=0}^s\frac{{(-1)}^{s+k}(s+k)!}{{(k!)}^2(s-k)!}{\zeta }^k({\tilde{P}}_0\left(\zeta \right)=1;{\tilde{P}}_1\left(\zeta \right)=2\zeta -1;s=2,3,4,\cdots ).$$

(8)

The function $\alpha \left(\zeta \right)\in {\mathcal{L}}_2[0,1]$ can be expressed and approximated as a linear combination of the first $(m+1)$ terms of ${\tilde{P}}_s\left(\zeta \right)$, as follows:

$$\alpha \left(\zeta \right)\simeq {\alpha }_m\left(\zeta \right)=\sum _{i=0}^m{a}_i{\tilde{P}}_i\left(\zeta \right),m=1,2,3,4,\cdots .$$

(9)

where the coefficients $a_i$ are given by

$$a_i=(2i+1){\int }_0^1\alpha \left(\zeta \right){\tilde{P}}_i\left(\zeta \right)d\zeta (i=0,1,2,\cdots ).$$

We now state the following useful theorem.

Theorem 1

([27]). Let $\alpha \left(\zeta \right)$ be approximated by the shifted Legendre polynomials in (9). Suppose also that $\nu >0$. Then

$$D^{\nu }\alpha \left(\zeta \right)=\sum _{i=\lceil \nu \rceil }^m\sum _{k=\lceil \nu \rceil }^i{a}_i{\chi }_{ik}^{\left(\nu \right)}{\zeta }^{k-\nu },$$

(10)

where

$${\chi }_{i,k}^{\left(\nu \right)}=\frac{(i+k)!\mathsf{\Gamma }(k+1)}{{(k!)}^2(i-k)!\mathsf{\Gamma }(k-\alpha +1)}.$$

We now implement the SCM for treating the given problem numerically and convert it into a non-linear system of algebraic equations. We will then solve the resulting system of algebraic equations using the Newton iteration method (see [28]).

In order to achieve this implementation of the SCM to solve numerically the suggested model given by (2)–(5), we approximate ${\alpha }_1\left(t\right),$ ${\alpha }_2\left(\zeta \right)$ and ${\alpha }_3\left(\zeta \right)$ by ${\alpha }_{1,m}\left(t\right),$${\alpha }_{2,m}\left(\zeta \right)$ and ${\alpha }_{3,m}\left(t\right)$, respectively, as follows:

$${\alpha }_{1,m}\left(\zeta \right)=\sum _{k=0}^m{a}_k{\tilde{P}}_k\left(\zeta \right),{\alpha }_{2,m}\left(\zeta \right)=\sum _{k=0}^m{b}_k{\tilde{P}}_k\left(\zeta \right),\mathrm{and}{\alpha }_{3,m}\left(\zeta \right)=\sum _{k=0}^m{c}_k{\tilde{P}}_k\left(\zeta \right).$$

(11)

The connection between the Equations (2)–(5), (10), and (11) leads us to the following consequences:

$$\sum _{i=\lceil \nu \rceil }^m\sum _{k=\lceil \nu \rceil }^i{a}_k{\chi }_{i,k}^{\left(\nu \right)}{\zeta }^{k-\nu }=\mu e^{-\epsilon \delta \zeta }-\left(\sum _{k=0}^m{a}_k{\tilde{P}}_k\left(\zeta \right)\right)\left(\sum _{k=0}^m{b}_k{\tilde{P}}_k\left(\zeta \right)\right),$$

(12)

$$\begin{array}{cc}\hfill & \sum _{i=\lceil \nu \rceil }^m\sum _{k=\lceil \nu \rceil }^i{a}_k{\chi }_{i,k}^{\left(\nu \right)}{\zeta }^{k-\nu }=\left(\sum _{k=0}^m{a}_k{\tilde{P}}_k\left(\zeta \right)\right)\left(\sum _{k=0}^m{b}_k{\tilde{P}}_k\left(\zeta \right)\right)-\frac{1}{\delta }\left(\sum _{k=0}^m{b}_k{\tilde{P}}_k\left(\zeta \right)\right)\hfill \\ \hfill & ·\left(\sum _{k=0}^m{b}_k{\tilde{P}}_k\left(\zeta \right)+\sum _{k=0}^m{b}_k{\tilde{P}}_k\left(\zeta \right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta \zeta },\right)+\lambda -1\right)\hfill \end{array}$$

(13)

and

$$\sum _{k=0}^m{c}_k{\tilde{P}}_k\left(\zeta \right)=\sum _{k=0}^m{a}_k{\tilde{P}}_k\left(\zeta \right)+\sum _{k=0}^N{b}_k{\tilde{P}}_k\left(\zeta \right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta \zeta }\right)+\lambda -1.$$

(14)

Moreover, the collocation for the system given by (12)–(14) at m of the nodes ${\zeta }_p$, which is, the roots of ${\tilde{P}}_k\left(\zeta \right)$, yields

$$\sum _{i=\lceil \nu \rceil }^m\sum _{k=\lceil \nu \rceil }^i{a}_k{\chi }_{i,k}^{\left(\nu \right)}{\zeta }_p^{k-\nu }=\mu e^{-\epsilon \delta {\zeta }_p}-\left(\sum _{k=0}^m{a}_k{\tilde{P}}_k\left({\zeta }_p\right)\right)\left(\sum _{k=0}^m{b}_k{\tilde{P}}_k\left({\zeta }_p\right)\right),$$

(15)

$$\begin{array}{cc}\hfill & \sum _{i=\lceil \nu \rceil }^m\sum _{k=\lceil \nu \rceil }^i{b}_k{\chi }_{i,k}^{\left(\nu \right)}{\zeta }_p^{k-\nu }=\left(\sum _{k=0}^m{a}_k{\tilde{P}}_k\left({\zeta }_p\right)\right)\left(\sum _{k=0}^N{b}_k{\tilde{P}}_k\left({\zeta }_p\right)\right)-\frac{1}{\delta }\left(\sum _{k=0}^m{b}_k{\tilde{P}}_k\left({\zeta }_p\right)\right)\hfill \\ \hfill & ·\left(\sum _{k=0}^m{a}_k{\tilde{P}}_k\left({\zeta }_p\right)+\sum _{k=0}^m{b}_k{\tilde{P}}_k\left({\zeta }_p\right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta {\zeta }_p},\right)+\lambda -1\right)\hfill \end{array}$$

(16)

and

$$\sum _{k=0}^m{c}_k{\tilde{P}}_k\left({\zeta }_p\right)=\sum _{k=0}^m{a}_k{\tilde{P}}_k\left({\zeta }_p\right)+\sum _{k=0}^m{b}_k{\tilde{P}}_k\left({\zeta }_p\right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta {\zeta }_p}\right)+\lambda -1.$$

(17)

Now, upon substituting from Equation (11) into the initial conditions given by (5), we obtain the following three equations:

$$\sum _{k=0}^m{\Delta }_k{a}_k=0,$$

(18)

$$\sum _{k=0}^m{\Delta }_k{b}_k=1,$$

(19)

$$\sum _{k=0}^m{\Delta }_k{c}_k=\lambda ,$$

(20)

where

$${\Delta }_k={\tilde{P}}_k\left(0\right)={(-1)}^k(2k+1)(k=0,1,\cdots ,m).$$

Equations (15)–(17), together with the three equations in (18), provide a system of $3(m+1)$ algebraic equations. The involved unknowns $a_k,b_k,c_k(k=0,1,\cdots ,m)$ will be obtained by solving the system given by (15)–(18) with the help of the Newton iteration method.

4. Runge-Kutta Methods

In this section, we apply the fractional Runge–Kutta method of order 2 and of order 4, which was used in [29,30], respectively. To this end, we construct the schemes for fractional model given by (2)–(4) with these algorithms.

4.1. Two-Stage Fractional Runge-Kutta Method

The numerical solutions of Equations (2)–(4), which are based on TSFRK, are given by

$${\alpha }_1\left({\zeta }_{n+1}\right)={\alpha }_1\left({\zeta }_n\right)+\frac{h^{\nu }}{2\mathsf{\Gamma }(1+\nu )}(K_{2,1}+K_{2,2}),$$

(21)

$${\alpha }_2\left({\zeta }_{n+1}\right)={\alpha }_2\left({\zeta }_n\right)+\frac{h^{\nu }}{2\mathsf{\Gamma }(1+\nu )}(M_{2,1}+M_{2,2}),$$

(22)

$${\alpha }_3\left({\zeta }_{n+1}\right)={\alpha }_1\left({\zeta }_n\right)+{\alpha }_2\left({\zeta }_n\right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta {\zeta }_n}\right)+\lambda -1,$$

(23)

where

$$\begin{array}{c}\hfill K_{2,1}=f\left({\zeta }_n,{\alpha }_1\left({\zeta }_n\right),{\alpha }_2\left({\zeta }_n\right)\right),\end{array}$$

(24)

$$\begin{array}{c}\hfill M_{2,1}=g\left({\zeta }_n,{\alpha }_1\left({\zeta }_n\right),{\alpha }_2\left({\zeta }_n\right)\right),\end{array}$$

(25)

$$\begin{array}{cc}\hfill K_{2,2}& =f({\zeta }_n+\frac{2\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }(1+2\nu )}{h}^{\nu },{\alpha }_1\left({\zeta }_n\right)+\frac{2\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }(1+2\nu )}{h}^{\nu }{K}_{2,1},{\alpha }_2\left({\zeta }_n\right)\hfill \\ \hfill & +\frac{2\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }(1+2\nu )}{h}^{\nu }{M}_{2,1}),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill M_{2,2}& =g({\zeta }_n+\frac{2\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }(1+2\nu )}{h}^{\nu },{\alpha }_1\left({\zeta }_n\right)+\frac{2\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }(1+2\nu )}{h}^{\nu }{K}_{2,1},{\alpha }_2\left({\zeta }_n\right)\hfill \\ \hfill & +\frac{2\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }(1+2\nu )}{h}^{\nu }{M}_{2,1}),\hfill \end{array}$$

(27)

$$\begin{array}{c}\hfill f\left(\zeta ,{\alpha }_1\left(\zeta \right),{\alpha }_2\left(\zeta \right)\right)=\mu e^{-\epsilon \delta \zeta }-{\alpha }_1\left(\zeta \right){\alpha }_2\left(\zeta \right)\end{array}$$

(28)

and

$$\begin{array}{cc}\hfill & g\left(\zeta ,{\alpha }_1\left(\zeta \right),{\alpha }_2\left(\zeta \right)\right)={\alpha }_1\left(\zeta \right){\alpha }_2\left(\zeta \right)\hfill \\ \hfill & -\frac{{\alpha }_2\left(\zeta \right)}{\delta }\left({\alpha }_1\left(\zeta \right)+{\alpha }_2\left(\zeta \right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta \zeta }\right)+\lambda -1\right).\hfill \end{array}$$

(29)

For further details about the derivation of this algorithm, see [29].

4.2. Four-Stage Fractional Runge-Kutta

In this subsection, we employ the algorithm that was derived in [30] on the fractional model given by (2)–(4), as in the above section with TSFRK, such that

$${\alpha }_1\left({\zeta }_{n+1}\right)={\alpha }_1\left({\zeta }_n\right)+\frac{1}{6}(K_{4,1}+2K_{4,2}+2K_{4,3}+K_{4,4}),$$

(30)

$${\alpha }_2\left({\zeta }_{n+1}\right)={\alpha }_2\left({\zeta }_n\right)+\frac{1}{6}(M_{4,1}+2M_{4,2}+2M_{4,3}+M_{4,4})$$

(31)

and

$${\alpha }_3\left({\zeta }_{n+1}\right)={\alpha }_1\left({\zeta }_n\right)+{\alpha }_2\left({\zeta }_n\right)-\frac{\mu }{\epsilon \delta }\left(1-e^{-\epsilon \delta {\zeta }_n}\right)+\lambda -1,$$

(32)

where

$$\begin{array}{c}\hfill K_{4,1}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}f\left({\zeta }_n,{\alpha }_1\left({\zeta }_n\right),{\alpha }_2\left({\zeta }_n\right))\right),\end{array}$$

(33)

$$\begin{array}{c}\hfill M_{4,1}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}g\left({\zeta }_n,{\alpha }_1\left({\zeta }_n\right),{\alpha }_2\left({\zeta }_n\right)\right),\end{array}$$

(34)

$$\begin{array}{c}\hfill K_{4,2}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}f\left({\zeta }_n+\frac{h^{\nu }}{2\mathsf{\Gamma }(1+\nu )},{\alpha }_1\left({\zeta }_n\right)+\frac{1}{2}{K}_{4,1},{\alpha }_2\left({\zeta }_n\right)+\frac{1}{2}{M}_{4,1}\right),\end{array}$$

(35)

$$\begin{array}{c}\hfill M_{4,2}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}g\left({\zeta }_n+\frac{h^{\nu }}{2\mathsf{\Gamma }(1+\nu )},{\alpha }_1\left({\zeta }_n\right)+\frac{1}{2}{K}_{4,1},{\alpha }_2\left({\zeta }_n\right)+\frac{1}{2}{M}_{4,1}\right),\end{array}$$

(36)

$$\begin{array}{c}\hfill K_{4,3}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}f\left({\zeta }_n+\frac{h^{\nu }}{2\mathsf{\Gamma }(1+\nu )},{\alpha }_1\left({\zeta }_n\right)+\frac{1}{2},K_{4,2},{\alpha }_2\left({\zeta }_n\right)+\frac{1}{2}{M}_{4,2}\right),\end{array}$$

(37)

$$\begin{array}{c}\hfill M_{4,3}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}g\left({\zeta }_n+\frac{h^{\nu }}{2\mathsf{\Gamma }(1+\nu )},{\alpha }_1\left({\zeta }_n\right)+\frac{1}{2}{K}_{4,2},{\alpha }_2\left({\zeta }_n\right)+\frac{1}{2}{M}_{4,2}\right),\end{array}$$

(38)

$$\begin{array}{c}\hfill K_{4,4}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}f\left({\zeta }_n+\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )},{\alpha }_1\left({\zeta }_n\right)+K_{4,3},{\alpha }_2\left({\zeta }_n\right)+M_{4,3}\right)\end{array}$$

(39)

and

$$\begin{array}{c}\hfill M_{4,4}=\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )}g\left({\zeta }_n+\frac{h^{\nu }}{\mathsf{\Gamma }(1+\nu )},{\alpha }_1\left({\zeta }_n\right)+K_{4,3},{\alpha }_2\left({\zeta }_n\right)+M_{4,3}\right).\end{array}$$

(40)

For further details on the four-stage fractional Runge–Kutta (FSFRK) method, see [30].

5. Numerical Results and Discussion

In this section, we investigate the numerical solutions of the fractional-order chemical clock reaction model by implementing the SCM, TSFRK, FSFRK, methods. We then compare the obtained results with these methods with the results of Lagrange polynomial interpolation (LPI) in [3]. We first conduct a test for all methods on fractional-order Logistic equation, which has an exact solution for integer as well as non-integer order. Figure 1a–d graphically illustrate the absolute error between the exact solution of FOLE and the corresponding numerical solutions by SCM, TSFRK, FSFRK, and LPI, respectively, for $\nu =1,$$\alpha \left(0\right)=0.75,h=0.03,m=21$ and $K=1.$ In these figures, we notice that the absolute error based upon the FSFRK are much better than those of the SCM, the TSFRK, and the LPI methods when the order of logistic equation is 1. In this case, the order of the error for FSFRK is ${10}^{-14}.$ In Figure 2a–d, we evaluate the absolute error for non-integer order of FOLE between the exact solution of FOLE and SCM, TSFRK, FSFRK, and LPI, respectively. In this case, we set $\nu =0.9,$$\alpha \left(0\right)=0.75,h=0.03,m=21$ and $K=1.$ In these figures, the orders of the SCM and the LPI are better than those of the TSFRK and FSFRK. The order of error for the SCM and the LPI is ${10}^{-2}$, whereas, in the case of the TSFRK and the FSFRK, is ${10}^{-1}$. Accordingly, based on this investigation, we can take the LPI [3] as a method to compare our numerical solution of the fractional-order clock chemical model based on the SCM, TSFRK and FSFRK. Additionally, we investigate in details the error in

Table 1. The absolute error between the exact solution (6) and the numerical solutions by spectral collocation method (SCM), two-stage fractional Runge–Kutta (TSFR), four-stage fractional Runge–Kutta (FSFRK), and Lagrange polynomial interpolation (LPI), respectively, for $\nu =1,\zeta =5,h=0.003,m=21$ with different values of K.

K	$|\mathbf{Exact}-\mathbf{SCM}|$	$|\mathbf{Exact}-\mathbf{TSFRK}|$	$|\mathbf{Exact}-\mathbf{FSFRK}|$	$|\mathbf{Exact}-\mathbf{LPI}|$
$0.2$	$1.53945\times {10}^{-17}$	$4.16541\times {10}^{-9}$	$3.33067\times {10}^{-16}$	$2.127988\times {10}^{-5}$
$0.4$	$6.72001\times {10}^{-15}$	$1.68964\times {10}^{-8}$	$1.88738\times {10}^{-15}$	$4.166692\times {10}^{-5}$
$0.6$	$6.696420\times {10}^{-13}$	$2.464570\times {10}^{-8}$	$3.663735\times {10}^{-15}$	$3.946389\times {10}^{-5}$
$0.8$	$6.34949\times {10}^{-11}$	$2.3406\times {10}^{-8}$	$7.43849\times {10}^{-15}$	$2.762099\times {10}^{-5}$
$1.0$	$1.281608\times {10}^{-10}$	$1.766755\times {10}^{-8}$	$9.325873\times {10}^{-15}$	$1.647932\times {10}^{-5}$

and

Table 2. The absolute error between the exact solution (6) and the numerical solutions by SCM, TSFR, FSFRK, and LPI, respectively, for $\nu =0.8,\zeta =5,h=0.003,m=21$ with different values of K.

K	$|\mathbf{Exact}-\mathbf{SCM}|$	$|\mathbf{Exact}-\mathbf{TSFRK}|$	$|\mathbf{Exact}-\mathbf{FSFRK}|$	$|\mathbf{Exact}-\mathbf{LPI}|$
$0.2$	$1.126692\times {10}^{-3}$	$1.08085\times {10}^{-1}$	$1.080896\times {10}^{-1}$	$1.184387\times {10}^{-3}$
$0.4$	$2.186600\times {10}^{-3}$	$6.576224\times {10}^{-2}$	$6.576269\times {10}^{-2}$	$2.257930\times {10}^{-3}$
$0.6$	$2.480864\times {10}^{-3}$	$4.453745\times {10}^{-2}$	$4.453749\times {10}^{-2}$	$2.553845\times {10}^{-3}$
$0.8$	$2.413834\times {10}^{-3}$	$3.293829\times {10}^{-2}$	$3.293830\times {10}^{-2}$	$2.486642\times {10}^{-3}$
$1.0$	$2.227433\times {10}^{-3}$	$2.591927\times {10}^{-2}$	$2.591927\times {10}^{-2}$	$2.301093\times {10}^{-3}$

. Table 1 and Table 2 represent the absolute error between the SCM, the TSFR, the FSFRK, and the LPI methods, respectively, with the exact solution of fractional-order logistic equation. In these table, the effect of the changing of the order of the differential equation from the integer order to the non-integer order is the most important parameter in this work. This, in fact, is the aim of this work. The effect of the parameter K is also studied. The time and parameters h and m are fixed. Through Table 1 and Table 2, we make the following observations:

• In Table 1, we notice that the order of the errors are small in all of the methods used here. Additionally, in the cases of the SCM and the FSFRK methods are much better than those in the cases of the TSFR and the LPI methods.
• In Table 2, we see that the orders of the errors for all of the methods are close to each other. However, the SCM and LPI methods are better than the TSFRK and FSFRK methods. Therefore, we cannot say here that one of the methods is the best absolutely, but we can say that some of the methods are better than the other methods. In Figure 3, Figure 4 and Figure 5, we compute the absolute error between LPI and SCM, TSFRK, and FSFRK, respectively. In these figures, we set $\nu =1,h=0.03,m=21,\lambda =2,ϵ=1,\mu =2,\delta =1$ and the initial values are set as ${\alpha }_1\left(0\right)=0,{\alpha }_2\left(0\right)=1$ and ${\alpha }_3\left(0\right)=\lambda .$ From these figures, we observe that the order of the error for the SCM, TSFRK, and FSFRK methods is ${10}^{-3}$. In Figure 6, Figure 7 and Figure 8, we find the absolute error for the same caption as those of Figure 3, Figure 4 and Figure 5, but with $\nu =0.9.$ In the non-integer order, we have obtained the same results, as illustrated in Figure 2. The order of the error of the SCM is better than those of the TSFRK and FSFRK. From this investigation, we can say that the numerical solutions based upon the SCM and LPI methods are remarkably accurate and effective in the case of non-integer order.

6.  Conclusions

In this paper, we have presented three algorithms to find the numerical solutions of the chemical clock model. We applied each of these algorithms to the fractional-order chemical clock reaction model. In order to verify the efficiency of the results, which we have investigated in this paper, we introduced the fractional-order logistic equation with its exact solution. We compared the numerical solutions based on the SCM, TSFRK, FSFRK, and LPI, respectively, with the exact solution of the fractional-order logistic equation for integer and non integer orders. We have noted that the error for all of the methods were very small. However, the SCM and the FSFRK are better than TSFRK and LPI when the order of the logistic equation was an integer. But, if the order is non-integer, the SCM and LPI are better than the TSFRK and FSFRK. The order of error here, as is evident from the figures and the Tables, with respect to SCM and LPI is ${10}^{-2}$. In the case of the TSFRK and FSFRK, the order of error is found to be ${10}^{-1}$. For the fractional-order chemical clock reaction model, we have also compared the results with the corresponding previously-published results for integer order and non-integer order. We have thereby noticed that the error in the case of integer order for SCM, TSFRK and FSFRK are small for all methods. However, for the case of non-integer order, the SCM is remarkably better than the TSFRK and FSFRK. Finally, in this paper, we presented an analytical comparative study between four methods to find the numerical solutions to fractional differential equations, and the accuracy of the Runge-Kutta methods as compared to other methods was studied. As for the Runge-Kutta methods of finding numerical solutions to fractional differential equations, they are still new. Additionally, as far as we know, there are only two or three studies in this direction. Therefore, in our next work, we will use these methods to find and study the behavior of the numerical solutions for systems of fractional differential equations. The various results, which we have presented in this paper, have indeed been numerically analyzed and graphically illustrated by means of different figures and tables.