Fractional Derivatives with the Power-Law and the Mittag–Leffler Kernel Applied to the Nonlinear Baggs–Freedman Model

Abstract

This paper considers the Freedman model using the Liouville–Caputo fractional-order derivative and the fractional-order derivative with Mittag–Leffler kernel in the Liouville–Caputo sense. Alternative solutions via Laplace transform, Sumudu–Picard and Adams–Moulton rules were obtained. We prove the uniqueness and existence of the solutions for the alternative model. Numerical simulations for the prediction and interaction between a unilingual and a bilingual population were obtained for different values of the fractional order.

1. Introduction

Interactions between groups that speak different languages are occurring continuously in several countries in the world due to globalization and cultural openness. Multilingualism is the use of more than one language, either by an individual speaker or by a community of speakers. A mathematical model portraying the interaction dynamics of a population considering bilingual components and a monolingual component was proposed in [1,2]. Baggs in [2] studied the condition under which the bilingual component could persist and conditions under which it could become extinct. The weakness of these models is that they do not take into account the degree of interest in time and also the memory of the interaction, meaning the recall of the original meeting or interaction or contact up to a particular period of time in the present. Fractional calculus is one of the most powerful mathematical tools used in recent decades to model real-world problems in many fields, such as science, technology and engineering. The Liouville–Caputo fractional derivative involves the power-law function. The Liouville–Caputo fractional-order derivative allows usual initial conditions when playing with the integral transform, for instance the Laplace transform [3,4,5]. Recently, Abdon Atangana and Dumitru Baleanu proposed two fractional-order operators involving the generalized Mittag–Leffler function. The generalized Mittag–Leffler function was introduced in the literature to improve the limitations posed by the power-law [6,7,8,9,10,11,12]. The two-parametric, three-parametric, four-parametric and multiple Mittag–Leffler functions were presented by Wiman, Prabhakar, Shukla and Srivastava in [13,14,15,16,17,18]. The kernel used in Atangana–Baleanu fractional differentiation appears naturally in several physical problems as generalized exponential decay and as a power-law asymptotic for a very large time [19,20,21,22,23,24]. The choice of this derivative is motivated by the fact that the interaction is not local, but global, and also, the trend observed in the field does not follow the power-law. The generalized Mittag–Leffler function completely induced the effect of memory, which is very important in the nonlinear Baggs–Freedman model.

Atangana and Koca in [25] studied the nonlinear Baggs and Freedman model. Starting from the integer-order Freedman model presented by [25], we have:

$$\begin{array}{c}{\mathcal{D}}_t{x}_1=(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right),\\ {\mathcal{D}}_t{x}_2=(A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right),\end{array}$$

(1)

where $0<A_i,M_i,\left|F_i\right|\le 1$ are the birth, death and emigration parameters for $i\in [1,2]$. $0<G_1\le 1$ is the infant language acquisition parameter, that proportion of births in the $x_2$ population raised unilingually. $0<\alpha \le 1$ is the non-infant language acquisition rate, the proportion of $x_1$ learning the $x_2$ language per unit time after infancy. The term $\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}$ describes that part of population $x_1$ lost to $x_2$ due to virtual predation on the part of $x_2$ [26].

The aim of this work is to obtain alternative representations of the Freedman model considering Liouville–Caputo and Atangana–Baleanu–Caputo fractional derivatives. The paper is organized as follows: Section 2 introduces the fractional operators. Alternative representations of the Freedman model are shown in Section 3. Finally, in Section 4, we conclude the manuscript.

2. Fractional Operators

The Liouville–Caputo fractional-order derivative of order $\gamma$ is defined by [27]:

$${}_a^C{\mathcal{D}}_t^{\gamma }\left\{f\left(t\right)\right\}=\frac{1}{\mathrm{\Gamma }(1-\gamma )}{\int }_a^t\dot{f}\left(\theta \right){(t-\theta )}^{-\gamma }d\theta ,0<\gamma \le 1.$$

(2)

The Laplace transform to Liouville–Caputo fractional-order derivative gives [27]:

$$\mathcal{L}{[}_a^C{\mathcal{D}}_t^{\gamma }f\left(t\right)]=S^{\gamma }F\left(S\right)-{\displaystyle \sum _{k=0}^{m-1}}{S}^{\gamma -k-1}{f}^{\left(k\right)}\left(0\right).$$

(3)

The Atangana–Baleanu–Caputo fractional-order derivative is defined as follows [19,20,21,22,24]:

$${}_a^{ABC}{\mathcal{D}}_t^{\gamma }\left\{f\left(t\right)\right\}=\frac{B\left(\gamma \right)}{1-\gamma }{\int }_a^t\dot{f}\left(\theta \right)E_{\gamma }\left[-\frac{\gamma }{1-\gamma }{(t-\theta )}^{\gamma }\right]d\theta ,0<\gamma \le 1,$$

(4)

where $B\left(\beta \right)=B\left(0\right)=B\left(1\right)=1$ is a normalization function and $E_{\gamma }$ is the Mittag–Leffler function [6,7,8,9,10,11,12]. The Mittag–Leffler kernel is a combination of both the exponential-law and power-law. For this fractional derivative, we have at the same time the power-law and the stretched exponential as the waiting time distribution.

The Laplace transform of Equation (4) is defined as follows:

$$\begin{array}{lll}\mathcal{L}{[}_a^{ABC}{\mathcal{D}}_t^{\gamma }f\left(t\right)]\left(s\right)& =& \frac{B\left(\gamma \right)}{1-\gamma }\mathcal{L}\left[{\int }_a^t\dot{f}\left(\theta \right)E_{\gamma }\left[-\gamma \frac{{(t-\theta )}^{\gamma }}{1-\gamma }\right]d\theta \right]\hfill \\ & =& \frac{B\left(\gamma \right)}{1-\gamma }\frac{s^{\gamma }\mathcal{L}\left[f\left(t\right)\right]\left(s\right)-s^{\gamma -1}f\left(0\right)}{s^{\gamma }+\frac{\gamma }{1-\gamma }}.\hfill \end{array}$$

(5)

The Sumudu transform is derived from the classical Fourier integral [28]. The Sumudu transform of Equation (4) is defined as:

$$ST\left\{{}_a^{ABC}{\mathcal{D}}_t^{\gamma }f\left(t\right)\right\}=\frac{B\left(\gamma \right)}{1-\gamma }\left(\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)\right)\times [ST\left(f\left(t\right)\right)-f\left(0\right)].$$

(6)

The Atangana–Baleanu fractional integral of order $\gamma$ of a function $f\left(t\right)$ is defined as:

$${}_a^{AB}{I}_t^{\gamma }f\left(t\right)=\frac{1-\gamma }{B\left(\gamma \right)}f\left(t\right)+\frac{\gamma }{B\left(\gamma \right)\mathrm{\Gamma }\left(\gamma \right)}{\int }_0^{t}f\left(s\right){(t-s)}^{\gamma -1}ds.$$

(7)

3. Freedman Model

In this section, we obtain alternative representations of the Freedman model considering the Liouville–Caputo fractional derivative, and the special solution is obtained using a Laplace transform method.

3.1. Freedman Model with the Power-Law Kernel

Considering Equation (2), the modified Freedman model with the power-law kernel is given as:

$$\begin{array}{c}{}_0^C{\mathcal{D}}_t^{\gamma }{x}_1\left(t\right)=(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right),\\ {}_0^C{\mathcal{D}}_t^{\gamma }{x}_2\left(t\right)=(A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right),\end{array}$$

(8)

where $0<\gamma \le 1$ is the fractional order, $x_1\left(t\right)$ represents the interaction of the majority unilingual population and $x_2\left(t\right)$ is a bilingual population.

Applying the Laplace transform operator (3) and the inverse Laplace transform on both sides of Equation (8), we obtain:

$$\begin{array}{c}{x}_1\left(t\right)=x_1\left(0\right)+{\mathcal{L}}^{-1}\left\{\frac{1}{s^{\gamma }}\mathcal{L}\left[(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha \frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right)\right]\left(s\right)\right\}\left(t\right),\\ x_2\left(t\right)=x_2\left(0\right)+{\mathcal{L}}^{-1}\left\{\frac{1}{s^{\gamma }}\mathcal{L}\left[(A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha \frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right)\right]\left(s\right)\right\}\left(t\right).\end{array}$$

(9)

The following iterative formula is then proposed:

$$\begin{array}{ll}{x}_{1\left(n\right)}\left(t\right)& ={\mathcal{L}}^{-1}\left\{\frac{1}{s^{\gamma }}\mathcal{L}\right[(A_1-M_1-F_1)x_{1(n-1)}\left(t\right)-L_1{x}_{1(n-1)}^2\left(t\right)\\ & -\alpha \frac{x_{1(n-1)}\left(t\right)x_{2(n-1)}\left(t\right)}{1+x_{1(n-1)}\left(t\right)}+G_1{A}_2{x}_{2(n-1)}\left(t\right)\left]\left(s\right)\right\}\left(t\right),\\ x_{2\left(n\right)}\left(t\right)& ={\mathcal{L}}^{-1}\left\{\frac{1}{s^{\gamma }}\mathcal{L}\right[(A_2-M_2-F_2)x_{2(n-1)}\left(t\right)-L_2{x}_{1(n-1)}^2\left(t\right)\\ & +\alpha \frac{x_{1(n-1)}\left(t\right)x_{2(n-1)}\left(t\right)}{1+x_{1(n-1)}\left(t\right)}-G_1{A}_2{x}_{2(n-1)}\left(t\right)\left]\left(s\right)\right\}\left(t\right),\end{array}$$

(10)

where,

$$x_{1\left(0\right)}\left(t\right)=x_1\left(0\right);x_{2\left(0\right)}\left(t\right)=x_2\left(0\right),$$

(11)

where the approximate solution is assumed to be obtained as a limit when “n” tend to infinity:

$$x_1\left(t\right)=\underset{n\to \infty }{\lim }{x}_{1\left(n\right)}\left(t\right);x_2\left(t\right)=\underset{n\to \infty }{\lim }{x}_{2\left(n\right)}\left(t\right).$$

(12)

3.2. Stability Analysis of the Iteration Method

Theorem 1.

We demonstrate that the recursive method given by Equation (10) is stable.

Proof. 

It is possible to find two positive constants Y and Z such that, for all:

$$0\le t\le T\le \infty ,\left|\right|x_1\left(t\right)\left|\right|<Y\mathrm{and}\left|\right|x_1\left(t\right)\left|\right|<Z.$$

(13)

Now, we consider a subset of $C_2\left((a,b)(0,T)\right)$ defined by:

$$H=\left\{\phi :(a,b)(0,T)\to H,\frac{1}{\mathrm{\Gamma }\left(\gamma \right)}\int {(t-\phi )}^{\gamma -1}v\left(\phi \right)u\left(\phi \right)d\phi <\infty \right\},$$

(14)

we now consider the operator $\varphi$ defined as:

$$\varphi (x_1,x_2)=\left\{\begin{array}{c}(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right),\hfill \\ (A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right).\hfill \end{array}\right.$$

Then:

$$=\left\{\begin{array}{c}<\varphi (x_1,x_2)-\varphi (X_1,X_2),\varphi (x_1-X_1,x_2-X_2)>,\hfill \\ <(A_1-M_1-F_1)(x_1\left(t\right)-X_1\left(t\right))-L_1{(x_1\left(t\right)-X_1\left(t\right))}^2\hfill \\ -\alpha ·\frac{(x_1\left(t\right)-X_1\left(t\right))(x_2\left(t\right)-X_2\left(t\right))}{1+(x_1\left(t\right)-X_1\left(t\right))}+G_1{A}_2(x_2\left(t\right)-X_2\left(t\right))>,\hfill \\ <(A_2-M_2-F_2)(x_2\left(t\right)-X_2\left(t\right))-L_2{(x_1\left(t\right)-X_2\left(t\right))}^2\hfill \\ +\alpha ·\frac{(x_1\left(t\right)-X_1\left(t\right))(x_2\left(t\right)-X_2\left(t\right))}{1+(x_1\left(t\right)-X_1\left(t\right))}-G_1{A}_2(x_2\left(t\right)-X_2\left(t\right))>,\hfill \end{array}\right.$$

where,

$$x_1\left(t\right)\ne X_1\left(t\right),\mathrm{and}{x}_2\left(t\right)\ne X_2\left(t\right).$$

(15)

Applying the absolute value on both sides, we have:

$$\begin{array}{c}<\varphi (x_1,x_2)-\varphi (X_1,X_2),\varphi (x_1-X_1,x_2-X_2)>,\\ =\left\{\begin{array}{c}\{(A_1-M_1-F_1)-L_1\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|\hfill \\ -\alpha ·\frac{\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|}{1+\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2}+G_1{A}_2\frac{\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|}{\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|}\left\}\right||x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2,\hfill \\ \{(A_2-M_2-F_2)-L_2\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|\hfill \\ +\alpha ·\frac{\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|}{1+\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2}-G_1{A}_2\left\}\right||x_2\left(t\right)-X_2\left(t\right){\left|\right|}^2.\hfill \end{array}\right.\end{array}$$

(16)

Then,

$$\begin{array}{c}<\varphi (x_1,x_2)-\varphi (X_1,X_2),\varphi (x_1-X_1,x_2-X_2)>,\\ <\left\{\begin{array}{c}\{(A_1-M_1-F_1)+L_1\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|\hfill \\ +\alpha ·\frac{\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|}{1+\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2}+G_1{A}_2\frac{\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|}{\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|}\left\}\right||x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2,\hfill \\ \{(A_2-M_2-F_2)+L_2\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|\hfill \\ +\alpha ·\frac{\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|}{1+\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2}+G_1{A}_2\left\}\right||x_2\left(t\right)-X_2\left(t\right){\left|\right|}^2,\hfill \end{array}\right.\end{array}$$

(17)

where,

$$\begin{array}{c}<\varphi (x_1,x_2)-\varphi (X_1,X_2),\varphi (x_1-X_1,x_2-X_2)>,\\ <\left\{\begin{array}{c}M\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2,\hfill \\ N\left|\right|x_2\left(t\right)-X_2\left(t\right){\left|\right|}^2,\hfill \end{array}\right.\end{array}$$

(18)

with:

$$M=(A_1-M_1-F_1)+L_1\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|+\alpha ·\frac{\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|}{1+\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2}+G_1{A}_2\frac{\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|}{\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|},$$

and:

$$N=(A_2-M_2-F_2)+L_2\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|+\alpha ·\frac{\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|}{1+\left|\right|x_1\left(t\right)-X_1\left(t\right){\left|\right|}^2}+G_1{A}_2.$$

(19)

Furthermore, if we consider a given non-null vector $(x_1,x_2)$, then using the some routine as the above case, we obtain:

$$\begin{array}{c}<\varphi (x_1,x_2)-\varphi (X_1,X_2),\varphi (x_1-X_1,x_2-X_2)>,\\ <\left\{\begin{array}{c}M\left|\right|x_1\left(t\right)-X_1\left(t\right)\left|\right|\left|\right|x_1\left(t\right)\left|\right|,\hfill \\ N\left|\right|x_2\left(t\right)-X_2\left(t\right)\left|\right|\left|\right|x_2\left(t\right)\left|\right|,\hfill \end{array}\right.\end{array}$$

(20)

From the results obtained in Equations (18) and (20), we conclude that the iterative method used is stable. This complete the proof.  ☐

Now, we consider the Adams method [29,30] to solve the system given by Equation (8). The basic idea of the n-step Adams–Bashforth method is to use a polynomial interpolation for $f(t,y(t\left)\right)$ passing through n points: $(t_i,f_i)$, $(t_{i-1},f_{i-1})$, …, $(t_{i-n+1},f_{i-n+1})$. Correspondingly, the n-step Adams–Moulton method uses a polynomial interpolation for $f(t,y(t\left)\right)$ passing through $n+1$ points: $(t_{i+1},f_{i+1})$, $(t_i,f_i)$, …, $(t_{i-n},f_{i-n})$.

The fractional Adams method is derived as follows [30]:

$$\begin{array}{c}{f}_{k+1}^P={\displaystyle \sum _{j=0}^{n-1}}{\displaystyle \frac{t_{k+1}^j}{j!}}{f}_0^{\left(j\right)}+{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\displaystyle \sum _{j=0}^k}{b}_{j,k+1}g(t_j,f_j),\hfill \\ f_{k+1}={\displaystyle \sum _{j=0}^{n-1}}{\displaystyle \frac{t_{k+1}^j}{j!}}{f}_0^{\left(j\right)}+{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}\left({\displaystyle \sum _{j=0}^k}{a}_{j,k+1}g(t_j,f_j)+a_{k+1,k+1}g(t_{k+1},f_{k+1}^P)\right),\hfill \end{array}$$

(21)

where,

$$\begin{array}{cc}{a}_{j,k+1}& ={\displaystyle \frac{h^{\gamma }}{\gamma (\gamma +1)}}·\left\{\begin{array}{cc}(k^{\gamma +1}-(k-\gamma ){(k+1)}^{\gamma })\hfill & j=0,\hfill \\ ({(k-j+2)}^{\gamma +1}+{(k-j)}^{\gamma +1}-2{(k-j+1)}^{\gamma +1})\hfill & 1\le j\le k,\hfill \\ 1\hfill & j=k+1,\hfill \end{array}\right.\hfill \\ b_{j,k+1}& ={\displaystyle \frac{h^{\gamma }}{\gamma }}({(k+1-j)}^{\gamma }-{(k-j)}^{\gamma }),j=0,1,2,...,k.\hfill \end{array}$$

(22)

Following this procedure, we can propose a numerical solution for System (8) using the Adams method (21) as follows:

$$\begin{array}{c}{x}_1\left(t\right)={\displaystyle \sum _{k=0}^{n-1}}{x}_1{\left(0\right)}^{\left(k\right)}{\displaystyle \frac{t^k}{k!}}+{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\displaystyle \underset{0}{\overset{t}{\int }}}{(t-u)}^{\gamma -1}\left[(A_1-M_1-F_1)x_1\left(u\right)-L_1{x}_1^2\left(u\right)-\alpha ·\frac{x_1\left(u\right)x_2\left(u\right)}{1+x_1\left(u\right)}+G_1{A}_2{x}_2\left(u\right)\right]du,\\ x_2\left(t\right)={\displaystyle \sum _{k=0}^{n-1}}{x}_2{\left(0\right)}^{\left(k\right)}{\displaystyle \frac{t^k}{k!}}+{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\gamma \right)}}{\displaystyle \underset{0}{\overset{t}{\int }}}{(t-u)}^{\gamma -1}\left[(A_2-M_2-F_2)x_2\left(u\right)-L_2{x}_2^2\left(u\right)+\alpha ·\frac{x_1\left(u\right)x_2\left(u\right)}{1+x_1\left(u\right)}-G_1{A}_2{x}_2\left(u\right)\right]du.\end{array}$$

(23)

3.3. Freedman Model with the Mittag–Leffler Kernel

Considering Equation (4), the modified Freedman model is given as:

$$\begin{array}{c}{}_0^{ABC}{\mathcal{D}}_t^{\gamma }{x}_1\left(t\right)=(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right),\\ {}_0^{ABC}{\mathcal{D}}_t^{\gamma }{x}_2\left(t\right)=(A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha ·\frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right),\end{array}$$

(24)

where $0<\gamma \le 1$ is the fractional order, $x_1\left(t\right)$ represents the interaction of the majority unilingual population and $x_2\left(t\right)$ is a bilingual population.

Now, we obtain an alternative solution using an iterative scheme. The technique involves coupling the Sumudu transform and its inverse. The Sumudu transform is an integral transform similar to the Laplace transform, introduced by Watugala to solve differential equations [28,29,30,31,32]. Applying the Sumudu transform (6) and the inverse Sumudu transform on both sides of the system (24) yields:

$$\begin{array}{c}{x}_1\left(t\right)=x_1\left(0\right)+S{T}^{-1}\left\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}·ST\left[(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha \frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right)\right]\right\},\\ x_2\left(t\right)=x_2\left(0\right)+S{T}^{-1}\left\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}·ST\left[(A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha \frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right)\right]\right\}.\end{array}$$

(25)

The following recursive formula for Equation (25) is obtained:

$$\begin{array}{ll}{x}_{1(n+1)}\left(t\right)=& x_{1\left(n\right)}\left(0\right)+S{T}^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\hfill \\ & ·ST[(A_1-M_1-F_1)x_{1\left(n\right)}\left(t\right)-L_1{x}_{1\left(n\right)}^2\left(t\right)\hfill \\ & -\alpha \frac{x_{1\left(n\right)}\left(t\right)x_{2\left(n\right)}\left(t\right)}{1+x_{1\left(n\right)}\left(t\right)}+G_1{A}_2{x}_{2\left(n\right)}\left(t\right)\left]\right\},\hfill \\ x_{2(n+1)}\left(t\right)=& x_{1\left(n\right)}\left(0\right)+S{T}^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\hfill \\ & ·ST[(A_2-M_2-F_2)x_{2\left(n\right)}\left(t\right)-L_2{x}_{2\left(n\right)}^2\left(t\right)\hfill \\ & +\alpha \frac{x_{1\left(n\right)}\left(t\right)x_{2\left(n\right)}\left(t\right)}{1+x_{1\left(n\right)}\left(t\right)}-G_1{A}_2{x}_{2\left(n\right)}\left(t\right)\left]\right\},\hfill \end{array}$$

(26)

and the solution of Equation (26) is provided by:

$$x_1\left(t\right)=\underset{n\to \infty }{\lim }{x}_{1\left(n\right)}\left(t\right);x_2\left(t\right)=\underset{n\to \infty }{\lim }{x}_{2\left(n\right)}\left(t\right).$$

(27)

3.4. Stability Analysis of the Iteration Method

Now, we provide in detail the stability analysis of this method and show the uniqueness of the special solutions using the fixed point theory and properties of the inner product and the Hilbert space, respectively.

Let $(X,|·\left|\right)$ be a Banach space and Ha self-map of X. Let $z_{n+1}=g(H,z_n)$ be a particular recursive procedure. The following conditions must be satisfied for $z_{n+1}=H{z}_n$.

• The fixed point set of H has at least one element.
• $z_n$ converges to a point $P\in F\left(H\right)$.
• ${\lim }_{n\to \infty }{x}_n\left(t\right)=P.$

Property 1.

Let $(X,|·\left|\right)$ be a Banach space and Ha self-map of X satisfying:

$$\left|\right|H_x-H_z\left|\right|\le \eta \left|\right|X-H_x\left|\right|+\eta \left|\right|x-z\left|\right|,$$

(28)

for all $x,z\in X$, where $0\le \eta$, $0\le \eta <1$. Suppose that H is Picard H-stable.

Considering the following recursive formula, we have:

$$\begin{array}{ll}{x}_{1(n+1)}\left(t\right)=& x_{1\left(n\right)}\left(0\right)+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\hfill \\ & ·S[(A_1-M_1-F_1)x_{1\left(n\right)}\left(t\right)-L_1{x}_{1\left(n\right)}^2\left(t\right)\hfill \\ & -\alpha \frac{x_{1\left(n\right)}\left(t\right)x_{2\left(n\right)}\left(t\right)}{1+x_{1\left(n\right)}\left(t\right)}+G_1{A}_2{x}_{2\left(n\right)}\left(t\right)\left]\right\},\hfill \\ x_{2(n+1)}\left(t\right)=& x_{2\left(n\right)}\left(0\right)+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\hfill \\ & ·S[(A_2-M_2-F_2)x_{2\left(n\right)}\left(t\right)-L_2{x}_{2\left(n\right)}^2\left(t\right)\hfill \\ & +\alpha \frac{x_{1\left(n\right)}\left(t\right)x_{2\left(n\right)}\left(t\right)}{1+x_{1\left(n\right)}\left(t\right)}-G_1{A}_2{x}_{2\left(n\right)}\left(t\right)\left]\right\},\hfill \end{array}$$

(29)

where:

$$\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)},$$

(30)

correspond to the fractional Lagrange multiplier.

Theorem 2.

Let K be a self-map defined as:

$$\begin{array}{ll}K\left[x_{1(n+1)}\left(t\right)\right]=& x_{1(n+1)}\left(t\right)=x_{1\left(n\right)}\left(t\right)+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\hfill \\ & ·S\left[(A_1-M_1-F_1)x_{1\left(n\right)}\left(t\right)-L_1{x}_{1\left(n\right)}^2\left(t\right)-\alpha \frac{x_{1\left(n\right)}\left(t\right)x_{2\left(n\right)}\left(t\right)}{1+x_{1\left(n\right)}\left(t\right)}+G_1{A}_2{x}_{2\left(n\right)}\left(t\right)\right]\},\hfill \\ K\left[x_{2(n+1)}\left(t\right)\right]=& x_{2(n+1)}\left(t\right)=x_{2\left(n\right)}\left(t\right)+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\hfill \\ & ·S\left[(A_2-M_2-F_2)x_{2\left(n\right)}\left(t\right)-L_2{x}_{2\left(n\right)}^2\left(t\right)+\alpha \frac{x_{1\left(n\right)}\left(t\right)x_{2\left(n\right)}\left(t\right)}{1+x_{1\left(n\right)}\left(t\right)}-G_1{A}_2{x}_{2\left(n\right)}\left(t\right)\right]\},\hfill \end{array}$$

(31)

is K-stable in $L^1(a,b)$ if:

$$\begin{array}{c}1+(A_1-M_1-F_1)f\left(\gamma \right)-L_1\left|\right|x_{1\left(n\right)}\left(t\right)+x_{1\left(m\right)}\left(t\right)\left|\right|d\left(\gamma \right)-\alpha \frac{ϵ\beta +ϵ+\theta }{\eta \rho }\omega \left(\gamma \right)+G_1{A}_{2}k\left(\gamma \right)<1,\\ 1+(A_2-M_2-F_2)i\left(\gamma \right)-L_2\left|\right|x_{2\left(n\right)}\left(t\right)+x_{2\left(m\right)}\left(t\right)\left|\right|r\left(\gamma \right)+\alpha \frac{ϵ\beta +ϵ+\theta }{\eta \rho }s\left(\gamma \right)-G_1{A}_{2}o\left(\gamma \right)<1,\end{array}$$

(32)

where $f\left(\gamma \right),d\left(\gamma \right),\omega \left(\gamma \right),k\left(\gamma \right),i\left(\gamma \right),r\left(\gamma \right),s\left(\gamma \right)$ and $o\left(\gamma \right)$ are functions from:

$$S^{-1}\left\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}S\right\}.$$

(33)

Proof. 

The proof consists of showing that K has a fixed point. To achieve this, we consider:

$$\begin{array}{c}K\left[x_{1(n+1)}\left(t\right)\right]-K\left[x_{1(m+1)}\left(t\right)\right]=x_{1\left(n\right)}\left(0\right)-x_{1\left(m\right)}\left(0\right)+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\\ ·S[((A_1-M_1-F_1)x_{1\left(n\right)}\left(t\right)-(A_1-M_1-F_1)x_{1\left(m\right)}\left(t\right))-(L_1{x}_{1\left(n\right)}^2\left(t\right)+L_1{x}_{1\left(m\right)}^2\left(t\right))\\ -\left(\alpha \frac{(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))·(x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right))}{1+(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)}\right)+G_1{A}_2(x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right))\left]\right\},\end{array}$$

and:

$$\begin{array}{c}K\left[x_{2(n+1)}\left(t\right)\right]-K\left[x_{2(m+1)}\left(t\right)\right]=x_{2\left(n\right)}\left(0\right)-x_{2\left(m\right)}\left(0\right)+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\\ ·S[((A_2-M_2-F_2)x_{2\left(n\right)}\left(t\right)-(A_2-M_2-F_2)x_{2\left(m\right)}\left(t\right))-(L_2{x}_{2\left(n\right)}^2\left(t\right)+L_2{x}_{2\left(m\right)}^2\left(t\right))\\ +\left(\alpha \frac{(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))·(x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right))}{1+(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)}\right)-G_1{A}_2(x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right))\left]\right\}.\end{array}$$

(34)

Using the properties of the norm and considering the triangular inequality, we get:

$$\begin{array}{c}\left|\right|K\left[x_{1\left(n\right)}\left(t\right)\right]-K\left[x_{1\left(m\right)}\left(t\right)\right]\left|\right|\le \left|\right|x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)\left|\right|+S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\\ ·S[((A_1-M_1-F_1)x_{1\left(n\right)}\left(t\right)-(A_1-M_1-F_1)x_{1\left(m\right)}\left(t\right))-(L_1{x}_{1\left(n\right)}^2\left(t\right)+L_1{x}_{1\left(m\right)}^2\left(t\right))\\ +\left(\alpha \frac{(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))·(x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right))}{1+(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)}\right)-G_1{A}_2(x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right))\left]\right\}.\end{array}$$

(35)

we consider that the solutions play the some role, i.e., $\left|\right|x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right)\left|\right|\cong \left|\right|x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)\left|\right|$.

Using the linearity of the inverse Sumudu transform, we get:

$$\begin{array}{c}\left|\right|K\left[x_{1\left(n\right)}\left(t\right)\right]-K\left[x_{1\left(m\right)}\left(t\right)\right]\left|\right|\le \left|\right|x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)\left|\right|\\ +S^{-1}\left\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}S\left(\left|\right|(A_1-M_1-F_1)[x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)]\left|\right|\right)\right\}\\ +S^{-1}\left\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}S\left(\left|\right|-L_1[x_{1\left(n\right)}^2\left(t\right)-x_{1\left(m\right)}^2\left(t\right)]\left|\right|\right)\right\}\\ +S^{-1}\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}\\ ·S\left(\left|\right|-\alpha \frac{(x_{1\left(m\right)}\left(t\right)\left(x_{1\left(n\right)}\left(t\right)\right)((x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))}{1+(x_{1\left(n\right)}\left(t\right)(1+x_{1\left(m\right)}\left(t\right))}+\frac{x_{1\left(n\right)}\left(t\right)(x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))+x_{2\left(n\right)}\left(t\right)(x_{2\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))}{1+(x_{1\left(n\right)}\left(t\right)(1+x_{1\left(m\right)}\left(t\right))}\left|\right|\right)\}\\ +S^{-1}\left\{\frac{1-\gamma }{B\left(\gamma \right)\gamma \mathrm{\Gamma }(\gamma +1)E_{\gamma }\left(-\frac{1}{1-\gamma }{u}^{\gamma }\right)}S\left|\right|G_1{A}_2[x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right))\left]\right||\right\}.\end{array}$$

(36)

Since $x_{1\left(n\right)}\left(t\right)$ and $x_{1\left(m\right)}\left(t\right)$ are bounded, we can find the following positive constants, $ϵ$, $\beta$, $\theta$, $\eta$ and $\rho$ such that for all t:

$$\begin{array}{c}\left|\right|x_{1\left(n\right)}\left(t\right)\left|\right|\le ϵ,\left|\right|x_{1\left(m\right)}\left(t\right)\left|\right|\le \beta ,\left|\right|1+x_{1\left(n\right)}\left(t\right)\left|\right|\le \rho ,\\ \left|\right|x_{2\left(n\right)}\left(t\right)\left|\right|\le \theta ,\left|\right|1+x_{1\left(m\right)}\left(t\right)\left|\right|\le \eta ,(n,m)\in N\times N.\end{array}$$

(37)

Considering Equations (36) and (37), we obtain:

$$\begin{array}{c}\left|\right|K\left[x_{1\left(n\right)}\left(t\right)\right]-K\left[x_{1\left(m\right)}\left(t\right)\right]\left|\right|\le \left|\right|x_{1\left(n\right)}\left(t\right)-x_{1\left(m\right)}\left(t\right)\left|\right|\\ ·\left(1+(A_1-M_1-F_1)f\left(\gamma \right)-L_1\left|\right|x_{1\left(n\right)}\left(t\right)+x_{1\left(m\right)}\left(t\right)\left|\right|d\left(\gamma \right)-\alpha \frac{ϵ\beta +ϵ+\theta }{\eta \rho }\omega \left(\gamma \right)+G_1{A}_{2}k\left(\gamma \right)\right),\end{array}$$

(38)

and:

$$\begin{array}{c}\left|\right|K\left[x_{2\left(n\right)}\left(t\right)\right]-K\left[x_{2\left(m\right)}\left(t\right)\right]\left|\right|\le \left|\right|x_{2\left(n\right)}\left(t\right)-x_{2\left(m\right)}\left(t\right)\left|\right|\\ ·\left(1+(A_2-M_2-F_2)i\left(\gamma \right)-L_2\left|\right|x_{2\left(n\right)}\left(t\right)+x_{2\left(m\right)}\left(t\right)\left|\right|r\left(\gamma \right)+\alpha \frac{ϵ\beta +ϵ+\theta }{\eta \rho }s\left(\gamma \right)-G_1{A}_{2}o\left(\gamma \right)\right),\end{array}$$

(39)

where $f\left(\gamma \right),d\left(\gamma \right),\omega \left(\gamma \right),k\left(\gamma \right),i\left(\gamma \right),r\left(\gamma \right),s\left(\gamma \right)$ and $o\left(\gamma \right)$ are functions from (33).

We next show that K satisfies Property 1. Consider Equations (38) and (39), yielding:

$$\eta (0,0),\eta =\left\{\begin{array}{c}1+(A_1-M_1-F_1)f\left(\gamma \right)-L_1\left|\right|x_{1\left(n\right)}\left(t\right)+x_{1\left(m\right)}\left(t\right)\left|\right|d\left(\gamma \right)-\alpha \frac{ϵ\beta +ϵ+\theta }{\eta \rho }\omega \left(\gamma \right)+G_1{A}_{2}k\left(\gamma \right),\hfill \\ 1+(A_2-M_2-F_2)i\left(\gamma \right)-L_2\left|\right|x_{2\left(n\right)}\left(t\right)+x_{2\left(m\right)}\left(t\right)\left|\right|r\left(\gamma \right)+\alpha \frac{ϵ\beta +ϵ+\theta }{\eta \rho }s\left(\gamma \right)-G_1{A}_{2}o\left(\gamma \right),\hfill \end{array}\right.$$

We conclude that K is Picard K-stable. ☐

3.5. Uniqueness of the Special Solution

Theorem 3.

We consider the Hilbert space $H=L^2((a,b)\times (0,k))$ that can be defined as the set of those functions:

$$v:(a,b)\times [0,T]\to \mathbb{R},\int \int uvdudv<\infty .$$

(40)

We now consider the following operator:

$$\eta (0,0),\eta =\left\{\begin{array}{c}(A_1-M_1-F_1)x_1\left(t\right)-L_1{x}_1^2\left(t\right)-\alpha \frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}+G_1{A}_2{x}_2\left(t\right),\hfill \\ (A_2-M_2-F_2)x_2\left(t\right)-L_2{x}_2^2\left(t\right)+\alpha \frac{x_1\left(t\right)x_2\left(t\right)}{1+x_1\left(t\right)}-G_1{A}_2{x}_2\left(t\right),\hfill \end{array}\right.$$

Proof. 

We prove that the inner product of:

$$(T(x_{11}\left(t\right)-x_{12}\left(t\right),x_{21}\left(t\right)-x_{22}\left(t\right),({\omega }_1,{\omega }_2)),$$

(41)

where $(x_{11}\left(t\right)-x_{12}\left(t\right)),x_{21}\left(t\right)-x_{22}\left(t\right)$ are special solutions of the system. We can assume that $(x_{11}\left(t\right)-x_{12}\left(t\right))\cong x_{21}\left(t\right)-x_{22}\left(t\right)$.

Using the relationship between the norm and the inner function, we get:

$$\begin{array}{c}\left((A_1-M_1-F_1)(x_{11}\left(t\right)-x_{12}\left(t\right))-L_1{(x_{11}\left(t\right)-x_{12}\left(t\right))}^2-\alpha \frac{(x_{11}\left(t\right)-x_{12}\left(t\right))(x_{21}\left(t\right)-x_{22}\left(t\right))}{1+(x_{11}\left(t\right)-x_{12}\left(t\right))}+G_1{A}_2(x_{21}\left(t\right)-x_{22}\left(t\right)),{\omega }_1\right)\\ \le (A_1-M_1-F_1)\left|\right|x_{11}-x_{12}\left|\right|\left|\right|{\omega }_1\left|\right|+L_1\left|\right|x_{11}-x_{12}{\left|\right|}^2\left|\right|{\omega }_1\left|\right|+\alpha \frac{\left|\right|x_{11}\left(t\right)-x_{12}\left(t\right){\left|\right|}^2}{\left|\right|1+(x_{11}\left(t\right)-x_{12}\left(t\right))\left|\right|}\left|\right|{\omega }_1\left|\right|+G_1{A}_2\left|\right|x_{11}-x_{12}\left|\right|\left|\right|{\omega }_1\left|\right|\end{array}$$

(42)

and:

$$\begin{array}{c}\left((A_2-M_2-F_2)(x_{21}\left(t\right)-x_{22}\left(t\right))-L_2{(x_{21}\left(t\right)-x_{22}\left(t\right))}^2+\alpha \frac{(x_{11}\left(t\right)-x_{12}\left(t\right))(x_{21}\left(t\right)-x_{22}\left(t\right))}{1+(x_{11}\left(t\right)-x_{12}\left(t\right))}+G_1{A}_2(x_{22}\left(t\right)-x_{21}\left(t\right)),{\omega }_2\right)\\ \le (A_2+M_2-F_2)\left|\right|x_{21}-x_{22}\left|\right|\left|\right|{\omega }_2\left|\right|+L_2\left|\right|x_{22}-x_{21}{\left|\right|}^2\left|\right|{\omega }_2\left|\right|+\alpha \frac{\left|\right|x_{22}\left(t\right)-x_{21}\left(t\right){\left|\right|}^2}{\left|\right|1+(x_{22}\left(t\right)-x_{21}\left(t\right))\left|\right|}\left|\right|{\omega }_2\left|\right|+G_1{A}_2\left|\right|x_{22}-x_{21}\left|\right|\left|\right|{\omega }_2\left|\right|;\end{array}$$

(43)

for large number m and n, both solutions converge to the exact solution; if $\eta =\left|\right|X_1-X_{11}\left|\right|,\left|\right|X_1-X_{12}\left|\right|$ and $\nu =\left|\right|x_2-x_{21}\left|\right|,\left|\right|x_2-x_{22}\left|\right|$, we have:

$$\eta <\frac{{\lambda }_n}{2\left((A_1-M_1-F_1)+L_1\left|\right|x_{11}\left(t\right)-x_{12}\left(t\right)\left|\right|+\alpha \frac{\left|\right|x_{11}\left(t\right)-x_{12}\left(t\right)\left|\right|}{1+x_{11}\left(t\right)-x_{12}\left(t\right)}+G_1{A}_2\right)\left|\right|{\omega }_1\left|\right|},$$

(44)

and:

$$\nu <\frac{{\lambda }_m}{2\left((A_2+M_2-F_2)+L_2\left|\right|x_{22}\left(t\right)-x_{21}\left(t\right)\left|\right|+\alpha \frac{\left|\right|x_{22}\left(t\right)-x_{21}\left(t\right)\left|\right|}{1-(x_{22}\left(t\right)-x_{21}\left(t\right))}+G_1{A}_2\right)\left|\right|{\omega }_2\left|\right|},$$

(45)

where ${\lambda }_n$ and ${\lambda }_m$ are two very small positive parameters.

Using the topology concept, we conclude that ${\lambda }_n<0$ and ${\lambda }_m<0$, where:

$$\begin{array}{c}\left((A_1-M_1-F_1)+L_1\left|\right|x_{11}\left(t\right)-x_{12}\left(t\right)\left|\right|+\alpha \frac{\left|\right|x_{11}\left(t\right)-x_{12}\left(t\right)\left|\right|}{1+x_{11}\left(t\right)-x_{12}\left(t\right)}+G_1{A}_2\right)\ne 0,\\ \left((A_2+M_2-F_2)+L_2\left|\right|x_{22}\left(t\right)-x_{21}\left(t\right)\left|\right|+\alpha \frac{\left|\right|x_{22}\left(t\right)-x_{21}\left(t\right)\left|\right|}{1-(x_{22}\left(t\right)-x_{21}\left(t\right))}+G_1{A}_2\right)\ne 0.\end{array}$$

(46)

This completes the proof. ☐

Involving the Atangana–Baleanu fractional integral, we can propose a numerical solution using the Adams–Moulton rule:

$${}_0^{AB}{\mathcal{I}}_t^{\gamma }\left[f\left(t_{n+1}\right)\right]=\frac{1-\gamma }{B\left(\gamma \right)}\left[\frac{f\left(t_{n+1}\right)-f\left(t_n\right)}{2}\right]+\frac{\gamma }{\mathrm{\Gamma }\left(\gamma \right)}{\displaystyle \sum _{k=0}^{\infty }}\left[\frac{f\left(t_{k+1}\right)-f\left(t_k\right)}{2}\right]b_k^{\gamma },$$

(47)

where:

$$b_k^{\gamma }={(k+1)}^{1-\gamma }-{\left(k\right)}^{1-\gamma }.$$

(48)

Considering the above numerical scheme, we have:

$$\begin{array}{c}{x}_{1(n+1)}\left(t\right)-x_{1\left(n\right)}\left(t\right)=x_{0\left(1\right)}^n\left(t\right)+\left\{\frac{1-\gamma }{B\left(\gamma \right)}\right[(A_1-M_1-F_1)\frac{x_{1(n+1)}\left(t\right)-x_{1\left(n\right)}\left(t\right)}{2}\\ -L_1\left(\frac{x_{1(n+1)}^2\left(t\right)-x_{1\left(n\right)}^2\left(t\right)}{2}\right)-\alpha \frac{\frac{(x_{1(n+1)}\left(t\right)-x_{1\left(n\right)}\left(t\right))(x_{2(n+1)}\left(t\right)-x_{2\left(n\right)}\left(t\right))}{2}}{1+(x_{1(n+1)}\left(t\right)-x_{1\left(n\right)}\left(t\right))}+G_1{A}_2\left(\frac{x_{2(n+1)}-x_{2\left(n\right)}\left(t\right)}{2}\right)]\\ +\frac{\gamma }{B\left(\gamma \right)}{\displaystyle \sum _{k=0}^{\infty }}{b}_k^{\gamma }·[(A_1-M_1-F_1)\frac{x_{1(k+1)}\left(t\right)-x_{1\left(k\right)}\left(t\right)}{2}-L_1\left(\frac{x_{1(k+1)}^2\left(t\right)-x_{1\left(k\right)}^2\left(t\right)}{2}\right)\\ -\alpha \frac{\frac{(x_{1(k+1)}\left(t\right)-x_{1\left(k\right)}\left(t\right))(x_{2(k+1)}\left(t\right)-x_{2\left(k\right)}\left(t\right))}{2}}{1+(x_{1(k+1)}\left(t\right)-x_{1\left(k\right)}\left(t\right))}-G_1{A}_2\left(\frac{x_{2(k+1)}-x_{2\left(k\right)}\left(t\right)}{2}\right)],\end{array}$$

(49)

and:

$$\begin{array}{c}{x}_{2(n+1)}\left(t\right)-x_{2\left(n\right)}\left(t\right)=x_{0\left(2\right)}^n\left(t\right)+\left\{\frac{1-\gamma }{B\left(\gamma \right)}\right[(A_2-M_2-F_2)\frac{x_{2(n+1)}\left(t\right)-x_{2\left(n\right)}\left(t\right)}{2}-\\ -L_2\left(\frac{x_{2(n+1)}^2\left(t\right)-x_{2\left(n\right)}^2\left(t\right)}{2}\right)+\alpha \frac{\frac{(x_{1(n+1)}\left(t\right)-x_{1\left(n\right)}\left(t\right))(x_{2(n+1)}\left(t\right)-x_{2\left(n\right)}\left(t\right))}{2}}{1+(x_{1(n+1)}\left(t\right)-x_{1\left(n\right)}\left(t\right))}+G_1{A}_2\left(\frac{x_{2(n+1)}-x_{2\left(n\right)}\left(t\right)}{2}\right)]\\ +\frac{\gamma }{B\left(\gamma \right)}{\displaystyle \sum _{k=0}^{\infty }}{b}_k^{\gamma }·[(A_2-M_2-F_2)\frac{x_{2(k+1)}\left(t\right)-x_{2\left(k\right)}\left(t\right)}{2}-L_2\left(\frac{x_{2(k+1)}^2\left(t\right)-x_{2\left(k\right)}^2\left(t\right)}{2}\right)\\ -\alpha \frac{\frac{(x_{1(k+1)}\left(t\right)-x_{1\left(k\right)}\left(t\right))(x_{2(k+1)}\left(t\right)-x_{2\left(k\right)}\left(t\right))}{2}}{1+(x_{1(k+1)}\left(t\right)-x_{1\left(k\right)}\left(t\right))}-G_1{A}_2\left(\frac{x_{2(n+1)}-x_{2\left(n\right)}\left(t\right)}{2}\right)].\end{array}$$

(50)

The proof of existence is described in detail by Alkahtani in [20].

4. Numerical Results

Example 1.

We present numerical simulations of the special solution of our model using the Adams–Moulton rule given by Equation (23) and Equations (49) and (50) for different arbitrary values of fractional order γ. We consider $A_1=0.017$; $A_2=0.30$; $M_1=0.06$; $M_2=0.007$; $L_1=0.01$; $L_2=0.004$; $F_1=0.3$; $F_2=0.7$; $G_1=0.01$; $\alpha =0.05$; and initial conditions $x\left(0\right)=10$; $y\left(0\right)=10$, arbitrarily chosen. The simulation time is 10 s, and the step size used in evaluating the approximate solution was $h=0.001$. The numerical results given in Figure 1a–d, Figure 2a–d, Figure 3a–d, and Figure 4a–d show numerical simulations of the special solution of our model as a function of time for different values of γ.

The figures show the interaction dynamics between a bilingual component and a monolingual component of a population in a particular environment. These numerical results show the influence of fractional order $\gamma$ in the prediction between the unilingual and bilingual population.

5. Conclusions

A Freedman model was considered using the fractional derivatives of the Liouville–Caputo and Atangana–Baleanu–Caputo types. The solutions of the alternative models were obtained using an iterative scheme based on the Laplace transform and the Sumudu transform. Furthermore, we employed the fixed point theorem to study the stability analysis of the iterative methods, and using properties of the inner product and the Hilbert space, the uniqueness of the special solution was presented in detail. Additionally, special solutions via the Adams–Moulton rule were obtained for both fractional derivatives. The results obtained using the Liouville–Caputo and Atangana–Baleanu-Caputo derivatives are exactly the same as the ordinary case. However, as $\gamma$ takes values smaller than one, the results obtained become a little different, having a remarkable difference when $\gamma <0.9$. This is due to the kernel involved in the definitions of the fractional derivative. The computer used for obtaining the results in this paper is an Intel Core i7, 2.6-GHz processor, 16.0 GB RAM (MATLAB R.2013a).