A Theoretical and Numerical Study on Fractional Order Biological Models with Caputo Fabrizio Derivative

Abstract

This article studies a biological population model in the context of a fractional Caputo-Fabrizio operator using double Laplace transform combined with the Adomian method. The conditions for the existence and uniqueness of solution of the problem under consideration is established with the use of the Banach principle and some theorems from fixed point theory. Furthermore, the convergence analysis is presented. For the accuracy and validation of the technique, some applications are presented. The numerical simulations present the obtained approximate solutions with a variety of fractional orders. From the numerical simulations, it is observed that when the fractional order is large, then the population density is also large; on the other hand, population density decreases with the decrease in the fractional order. The obtained results reveal that the considered technique is suitable and highly accurate in terms of the cost of computing, and can be used to analyze a wide range of complex non-linear fractional differential equations.

1. Introduction

For a long time, the theory of differential equations has been employed as a promising approach in several disciplines of science, engineering, and technology. Differential equations specifically, partial differential equations of various orders emerge in many realistic conditions, such as Brownian motion, fluid dynamic systems, populace increase, and numerous challenges to highways and gas dynamics [1,2]. The classical operators are local and do not have the property to preserve memory about the past. Due to this, these operators were not applicable to model certain real-world phenomena such as parabolic equations, groundwater flow equations, etc. On the other hand, non-integer order operators are non-local and memory preserving.

Over the last few decades, fractional calculus (FC) has been frequently considered in many areas, such as astronomy, viscoelasticity, mathematical biology, electrochemistry, signal processing, physics, economics, and social sciences [3,4,5,6]. It has been successfully used as a promising area to model real-world phenomena [7,8]. The fractional-order derivatives and integrals play a vital role in analyzing a specific problem in an interval, which is an advantage over classical derivatives and integrals. Various fractional operators, which include Riemann–Liouville, Hilfer, and Caputo, have been developed in FC to study real-world systems that are defined by convolution of the power-law with singular kernels [3,4]. This idea was further generalized to exponential decay type kernel which is the convolution of the exponential law with non-local kernels [9]. Since the Mittag–Leffler function is the generalized version of the exponential function, so further generalization of the fractional operators was made using the Mittag–Leffler kernel [10], which is non-singular and non-local in nature. Further study revealed a new type of fractional operator, combining Caputo and Proportional derivatives [11]. The aforesaid operators have a lot of advantages because due to the complexities of fractional-order nonlinear differential equations (FONDEs) integer-order operators cannot handle most of the problems to obtain explicit solutions [12]. These operators play a significant role in analyzing different types of FONDEs. Motivated by the above literature, we investigate the following biological model of fractional order with the CF operator [13,14].

$$\begin{array}{cc}\hfill & {}^{CF}{D}_t^{\alpha }u(x,y,t)=(D_x^2+D_y^2)u^2+k{u}^a-kr{u}^{(a+b)},\hfill \\ & t>0,x,y\in R,0\le \alpha \le 1\hfill \\ & u(x,y,0)=f_0.\hfill \end{array}$$

(1)

The authors in [15], presented the numerical solution of fractional order biological model using Caputo derivative in 2009. In 2012, D. Kumar et al. provided an algorithm of the Homotopy analysis method using Caputo fractional derivative [16].

In the fractional order derivatives, the most famous two kernels are singular and non-singular kernels. The prior one cannot model certain phenomena very effectively, that is why the new non-singular kernel was presented in [9]. This operator was found to be very efficient in predicting the dynamics of a biological model, for instance for the benefits of the aforesaid operator we can study the research works [17,18]. Inspired by the previous works, in this article we examine a general fractional-order biological population model in the sense of CF operator. We apply fixed point theorems to achieve existence and unique results. For the analytical solution of the considered model (1), we use the double Laplace decomposition method (DLDM), which is the combination of the double Laplace transform and decomposition method. We provide numerical examples for the accuracy and validation of our results.

Adomian proposed the Adomian decomposition method in 1980, which is an excellent approach for finding numerical and explicit solutions to a wide range of differential equations reflecting physical issues. This approach is effective for both initial and boundary value issues, PDEs and ODEs, linear, nonlinear equations, and also for stochastic systems. This approach does not require any perturbation or linearizations. ADM contributed significantly to the analytical solutions of nonlinear equations and the solution of FNDEs. We use the Double Laplace transform approach in this research work, which is a strong tool in engineering and mathematical analysis. We use this approach to transform fractions into algebraic equations, which could then be solved using DLADM (see [19,20]). The only limitation of the considered method is that it requires some space and analytical work which is a little bit time consuming.

The article is organized as follows. In Section 2 some definitions of the technique and fractional calculus are presented. Section 3 presents the existence and unique results of the proposed model. Section 4 presents the convergence analysis of the technique. In Section 5, the applications of the technique are presented while in Section 6 iterative examples are included. In Section 6.1, the numerical discussion is presented, while Section 7 concludes the article.

2. Preliminaries

Here we include some definitions regarding this article.

Definition 1

([9]). Consider that $\Theta \in {\mathrm{H}}^1[0,T]$ is a function, such that $0<T$, $\delta \in (0,1],$ then the definition of CFFD is

$${}^{CF}{D}^{\delta }(\Theta \left(t\right))=\frac{\mathfrak{M}\left(\delta \right)}{(1-\delta )}{\int }_0^t{\Theta }^{\prime }\left(t\right)exp\left[-\delta \frac{t-\rho }{1-\delta }\right]d\rho ,$$

where $\mathfrak{M}\left(\delta \right)$ is $\mathfrak{M}\left(\delta \right)=\frac{2}{2-\delta },0<\delta \le 1$. Further $\mathfrak{M}\left(0\right)=1.$ If Θ ∉ in ${\mathrm{H}}^1(0,T)$, then the CF operator is

$${}^{CF}{D}^{\delta }(\Theta \left(t\right))=\frac{\mathfrak{M}\left(\delta \right)}{(1-\delta )}{\int }_0^t(\Theta \left(t\right)-\Theta \left(\rho \right))exp\left[-\delta \frac{t-\rho }{1-\delta }\right]d\rho .$$

Definition 2

([9]). Consider $\delta \in (0,1],$ then CF integral with the order δ of Θ is

$${}^{CF}{I}^{\delta }[\Theta \left(t\right)]=G\Theta \left(t\right)+\overline{G}{\int }_0^t\Theta \left(\rho \right)d\rho ,t\ge 0.$$

When $\delta =1$, one can obtain classical integral of Θ, here

$$G=\frac{(1-\delta )}{\mathfrak{M}\left(\delta \right)},\overline{G}=\frac{\delta }{\mathfrak{M}\left(\delta \right)}.$$

Definition 3

([21]). Let $\varphi (x,t)$ be a function and for $x,t>0$ defined in the xt-plane, then the DLT of $\varphi (x,t)$ is presented as

$$\begin{array}{c}\hfill {\mathcal{L}}_x{\mathcal{L}}_t\left[\Theta (x,t)\right]={\int }_0^{\infty }{e}^{-px}{\int }_0^{\infty }{e}^{-st}\Theta (x,t)dtdx,\end{array}$$

where, p and s are complex numbers.

Definition 4.

The DLT of CF operator is

$${\mathcal{L}}_x{\mathcal{L}}_t\left\{{}^{CF}{D}_x^{\delta +n}\Theta (x,t)\right\}=\frac{M\left(\delta \right)}{p+(1-p)\delta }\left[p^{n+1}\overline{\Theta }(p,s)-\sum _{k=0}^n{p}^{n-k}{\mathcal{L}}_t\left\{\frac{{\partial }^k\Theta (0,t)}{\partial x^k}\right\}\right],$$

and

$${\mathcal{L}}_x{\mathcal{L}}_t\left\{{}^{CF}{D}_t^{\delta +m}\Theta (x,t)\right\}=\frac{M\left(\delta \right)}{s+(1-s)\delta }\left[s^{m+1}\overline{\Theta }(p,s)-\sum _{k=0}^m{s}^{m-k}{\mathcal{L}}_x\left\{\frac{{\partial }^k\Theta (x,0)}{\partial t^k}\right\}\right],$$

where, $n=\left[\delta \right]+1,m=\left[\delta \right]+1$.

Lemma 1

([22]).

$$\begin{array}{cc}& {}^{CF}{D}_t^{\delta }u(x,y,t)=(D_x^2+D_y^2)u^2+kr{u}^{(a+b)}k{u}^a,t>0,x,y\in R,0\le \delta \le 1\hfill \\ & u(x,y,0)=f_0\in R,\hfill \end{array}$$

in the terms of the following integral the above equation has a solution as

$$u\left(t\right)=u_0+\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}(\mathrm{E}(t,\omega \left(t\right))-\mathrm{E}(0,\omega \left(0\right)))+\frac{\left(2\delta \right)}{(2-\zeta )M\left(\delta \right)}{\int }_0^{t}f(\tau ,u\left(\tau \right))d\tau .$$

3. Existence Uniqueness of Results for Fractional Order Biological Model (1)

To present the results, we provide the following notions and lemma.

Consider $\vartheta =[0,1]$ and $\mathbb{C}\left(\vartheta \right)$ represents the space containing function which is continuous on $\vartheta$. Furthermore, consider set $B=\mathrm{Q}\left(t\right)/\mathrm{Q}\left(t\right)\in \mathbb{C}\left(\vartheta \right)$ with ${\left|\right|\mathrm{Q}\left(t\right)\left|\right|}_b\le ma{x}_{t\in \vartheta }\left|\mathrm{Q}\left(t\right)\right|$ represents Banach space. By Lemma 1, model (1) in integral form is

$$\mathrm{Q}\left(t\right)={\mathrm{Q}}_0+\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}(\mathrm{E}(t,\omega \left(t\right))-\mathrm{E}(0,\omega \left(0\right)))+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mathrm{E}(\xi ,\mathcal{Q}\left(\xi \right))d\xi .$$

Let suppose an operator $\mathfrak{T}:B\to B$ defined as

$$\mathfrak{T}\mathrm{Q}\left(t\right)={\mathrm{Q}}_0+\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}(\mathrm{E}(t,\omega \left(t\right))-\mathrm{E}(0,\omega \left(0\right)))+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mathrm{E}(\xi ,\mathcal{Q}\left(\xi \right))d\xi .$$

then operator $\mathfrak{T}$ has same fixed-point (FP) as (1).

Theorem 1.

Suppose that $f:\vartheta \times \mathcal{R}\to \mathcal{R}$ is continuous. Furthermore, consider in the following at-least one is satisfied.

$\left({\mathcal{H}}_1\right)$

Let $g\left(t\right)\in L[0,1]$ be the exists function which is non-negative, ∋

$$\left|\mathrm{E}(t,x)\right|\le h\left(t\right)+c_0{\left|x\right|}^{\xi },here,c_0\ge 0,0<\xi <1.$$

$\left({\mathcal{H}}_2\right)$

The function $\mathrm{E}$ satisfies $\left|\mathrm{E}(t,x)\right|\le c_0{\left|x\right|}^{\xi }$, where $c_0>0,\xi >1$. Then model (1) has one solution.

Proof. 

To prove the results, we use the Schauder FP theorem. To achieve our goal consider that $\left({\mathcal{H}}_1\right)$ is satisfied. Let us consider $G=\left\{\mathrm{Q}\left(t\right)\right|\mathrm{Q}\left(t\right)\in B,\parallel \mathrm{Q}\left(t\right){\parallel }_B\le k,t\in \vartheta \}$, where $k\ge max{\left(2A{c}_0\right)}^{\frac{1}{1-\xi }},2l$ and $l={max}_{y\in \vartheta }\left({\mathrm{Q}}_0+\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)g\left(t\right)}+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\left|g\left(\mathrm{t}\right)\right|d\mathrm{t}\right)$. Clearly in B, G is a ball. Furthermore, we prove that $\mathfrak{T}:G\to G$.

∀$u\in G$ we obtain

$$\begin{array}{cc}& \left|\mathfrak{T}\mathrm{Q}\left(t\right)\right|=|{\mathrm{Q}}_0+\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}(\mathrm{E}(t,\omega \left(t\right))-\mathrm{E}(0,\omega \left(0\right)))+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mathrm{E}(\mathrm{t},\mathcal{Q}\left(\mathrm{t}\right))d\mathrm{t}|\hfill \\ & \le {\mathrm{Q}}_0+\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\left|\mathrm{E}(t,\omega \left(t\right))\right|+\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\left|\mathrm{E}(0,\omega \left(0\right))\right|+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mathrm{E}(\mathrm{t},\mathcal{Q}\left(\mathrm{t}\right))d\mathrm{t}\hfill \\ & \le {\mathrm{Q}}_0+\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)}\left(g\left(t\right)+c_0{k}^{\xi }\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t(g\left(t\right)+c_0{k}^{\xi })d\mathrm{t}\hfill \\ & \le {\mathrm{Q}}_0+\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)}\left(g\left(t\right)+c_0{k}^{\xi }\right)+\frac{2\delta c_0{k}^{\xi }t}{(2-\delta )M\left(\delta \right)}+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t(g\left(\mathrm{t}\right))d\mathrm{t}\hfill \\ & \le {\mathrm{Q}}_0+\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)}g\left(t\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\left(g\left(\mathrm{t}\right)\right)d\mathrm{t}+\left(\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)}+\frac{2\delta t}{(2-\delta )M\left(\delta \right)}\right)c_0{k}^{\xi }\hfill \\ & \le {\mathrm{Q}}_0+\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)}g\left(t\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\left(g\left(\mathrm{t}\right)\right)d\mathrm{t}+\left(\frac{4(1-\delta )}{(2-\delta )M\left(\delta \right)}+\frac{2c_0{k}^{\xi }}{M\left(\delta \right)}\right).\hfill \end{array}$$

Therefore

$$\begin{array}{ccc}& & {\left|\right|\mathfrak{T}\mathrm{Q}\left(t\right)\left|\right|}_B=\underset{t\in \vartheta }{max}\left|\mathfrak{T}\mathrm{Q}\left(t\right)\right|\le l+\frac{2c_0{k}^{\xi }}{M\left(\delta \right)}=l+A{c}_0{k}^{\xi }\le \frac{k}{2}+\frac{k}{2}=k.\hfill \end{array}$$

Hence $T\mathrm{Q}\left(t\right)$ is continuous on $\vartheta$.

Next, suppose that $\left({\mathcal{H}}_2\right)$ is also satisfied. Selecting $0\le k\le \left(\frac{1}{A{c}_0}\right)\left(\frac{1}{\xi -1}\right)$. Similarly, on repeating the procedure as used above, we obtain

$${\left|\right|\mathfrak{T}\mathrm{Q}\left(t\right)\left|\right|}_B\le A{c}_0{k}^{\xi }\le k.$$

As a result, we obtain $\mathfrak{T}:G\to G$, clearly, that the operator $\mathfrak{T}$ is continuous as a result of the continuity $\mathrm{E}$.

Next we show that the operator T is completely continuous. Let $R={max}_{t\in \vartheta }\left|\mathrm{E}(t,\omega \left(t\right))\right|$, for any $\vartheta \in G$. Let $t_1,t_2\in \vartheta$ such that $t_1<t_2$.

Furthermore, let $U_1=\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}$ and $U_2=\frac{2\delta }{(2-\delta )M\left(\delta \right)}$, we obtain

$$\begin{array}{cc}\hfill & \mid T\mathrm{E}\left(t_2\right)-T\mathrm{E}\left(t_1\right)\mid =\mid {\mathrm{Q}}_0+U_1[\mathrm{Q}(t_2,\mathrm{E}\left(t_2\right))-\mathrm{Q}(0,\omega \left(0\right))]+U_2{\int }_0^{t_2}f(t,u\left(t\right))dt\hfill \\ & -{\mathrm{Q}}_0-U_1[\mathrm{Q}(t_1,\mathrm{E}\left(t_1\right))-\mathrm{Q}(0,\omega \left(0\right))]+U_2{\int }_0^{t_2}\mathrm{Q}(t,u\left(t\right))dt\mid \hfill \\ & =\mid U_1[\mathrm{Q}(t_2,\mathrm{E}\left(t_2\right))-\mathrm{Q}(t_1,\mathrm{E}\left(t_1\right))]+U_2{\int }_{t_1}^{t_2}\mathrm{E}(t,u\left(t\right))dt\mid \hfill \\ & \le U_1\mid \mathrm{Q}(t_2,\mathrm{E}\left(t_2\right))\mid +U_1\mid \mathrm{Q}(t_1,\mathrm{E}\left(t_1\right))\mid +U_2{\int }_{t_1}^{t_2}\mid \mathrm{Q}(t,u\left(t\right))\mid dt\hfill \\ & \le 2R{U}_1+R{U}_2{\int }_{t_1}^{t_2}dt=R(2U_1+U_2(t_2-t_1)).\hfill \end{array}$$

According to uniform continuity of the function $(t_2-t_1)$ on interval $\vartheta$, we obtain that $\mathfrak{T}G$ is equi continuous set. It is also observed that the function is uniformly bounded and $\mathfrak{T}G\subseteq G$, thus $\mathfrak{T}$ is continuous completely. Therefore, using the Schauder FP theorem, ∃ a solution for Equation (1) in the set G. □

Corollary 1.

Consider a function $\mathrm{Q}$ is bounded continuous on $\vartheta \times \mathrm{R}$, then Equation (1) has a solution.

Proof. 

As $\mathrm{Q}$ is continuous as well as bounded on $\vartheta \times \mathrm{R}$, ∃$L>0$, satisfying $\left|\mathrm{Q}\right|<L$. Consider $h\left(t\right)=L,c_0=0$ in $\left({\mathcal{H}}_1\right)$ of Theorem 1, then the model (1) has a solution.

Next we use the Banach-contraction principle to establish uniqueness results for solutions to (1). □

Theorem 2.

Consider a function that $\mathrm{Q}:\vartheta \times \mathrm{R}\to \mathrm{R}$ is continuous, and also satisfies the following conditions.

$\left({\mathcal{H}}_3\right)$

Consider a positive function $h\left(t\right)\in L[0,1]$ exists, ∋

$$\mid \mathrm{Q}(t,x)-\mathrm{Q}(t,y)\mid \le h\mid x-y\mid ,t\in [0,1],$$

also function $\mathrm{Q}$ satisfying $\mathrm{Q}(t,0)=0$

$\left({\mathcal{H}}_4\right)$

Consider that $\xi ={max}_{t\in \vartheta }\mid \frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}h\left(t\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\left|h\left(\mathrm{t}\right)\right|d\mathrm{t}\mid <1$, then model (1) has one solution.

Proof. 

Let the operator $\mathfrak{T}$ can be represented as

$$\begin{array}{cc}\hfill & \mathfrak{T}\mathrm{Q}\left(t\right)=\varphi +\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mathrm{E}(t,\omega \left(t\right))+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\left|\mathrm{E}(t,u\left(t\right))\right|dt,\hfill \end{array}$$

where $\varphi ={\mathrm{Q}}_0-\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mathrm{E}(0,\omega \left(0\right))$. For $\mathrm{Q}\left(t\right)\in B$, we obtain

$$\begin{array}{cc}\hfill & \mid \mathfrak{T}\mathrm{Q}\left(t\right)\mid =\mid \varphi +\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mathrm{E}(t,\omega \left(t\right))+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mid \mathrm{E}(t,u\left(t\right))\mid dt,\hfill \\ & \le \mid \varphi \mid +\mid \frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mathrm{E}(t,\omega \left(t\right))-\mathrm{E}(t,0)\mid +\mid \frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mid \mathrm{E}(t,u\left(t\right))\mid dt\hfill \\ & \le \mid \varphi \mid +\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}h\left(t\right)\mid \mathrm{Q}\left(t\right)\mid +\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^{t}h\left(t\right)\mid u\left(t\right)dt\hfill \\ & \le \mid \varphi \mid +\left(\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}h\left(t\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^{t}h\left(t\right)dt\right)\left|\right|u\left|\right|,\hfill \end{array}$$

we have

$$\begin{array}{cc}\hfill & {\Vert \mathfrak{T}\mathrm{Q}\left(t\right)\Vert }_B\le \mid \varphi \mid +\left(\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}h\left(t\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^{t}h\left(t\right)dt\right)\left|\right|u\left|\right|\hfill \\ & \le \mid \varphi \mid +\xi \Vert u\Vert \le \Vert u\Vert .\hfill \end{array}$$

Let $\mathrm{Q}\left(t\right),\nu \left(t\right)\in B$ we have

$$\begin{array}{cc}\hfill & \mid \mathfrak{T}\mathrm{Q}\left(t\right)-\mathfrak{T}\nu \left(t\right)\mid =\mid \Phi +\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mathrm{E}(t,\omega \left(t\right))+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\left|\mathcal{E}(t,\mathcal{Q}\left(t\right))\right|dt\hfill \\ & -\Phi -\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}f(t,v\left(t\right))-\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t|\mathcal{E}(t,v\left(t\right))dt\mid \hfill \\ & \le \frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mid \mathrm{E}(t,\omega \left(t\right))-f(t,v\left(t\right))\mid +\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mid \mathcal{E}(t,\mathcal{Q}\left(\mathrm{t}\right))-(t,v\left(t\right))dt\mid \hfill \\ & \le \frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}\mid \mathrm{Q}\left(t\right)-v\left(t\right)\mid +\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mid \mathcal{Q}\left(\mathrm{t}\right)-v\left(\mathrm{t}\right)\mid d\mathrm{t}\hfill \\ & \le \left(\frac{2(1-\delta )}{(2-\delta )M\left(\delta \right)}h\left(t\right)+\frac{\left(2\delta \right)}{(2-\delta )M\left(\delta \right)}{\int }_0^t\mid h\left(t\right)dt\right)\mid \mathcal{Q}\left(t\right)-v\left(t\right)\mid \hfill \\ & \le \Phi \left|\right|\mathcal{Q}\left(t\right)-v\left(t\right)\left|\right|\le \left|\right|\mathcal{Q}\left(t\right)-v\left(t\right)\left|\right|\hfill \end{array}$$

Finally, in view of $\xi <1$, $\mathfrak{T}$ is contraction. As a result, $\mathfrak{T}$ has only one fixed point according to Banach contraction principle. □

4. Convergence of MDLDM for Considered System

In this part, we consider the convergence of the MDLDM for the given problem Equation (1). To do this, we use Equation (1) in the operator as

$$\begin{array}{c}\hfill \mathcal{T}\left(u\right)=D_t\left(u\right)=(D_x^2+D_y^2)u^2+k{u}^a-kr{u}^{(a+b)}.\end{array}$$

(2)

Let $\mathrm{H}\in {\mathrm{L}}^2\left(\mathcal{T}\right)$∀$u\in \mathrm{H}$, where $\mathrm{H}={\mathrm{L}}_u^2\left[(m,n)\times [0,T]\right]$ [23], such that

$$\begin{array}{c}\hfill u:=\left[(m,n)\times [0,T]\right]\to {\mathbb{R}}^3,\end{array}$$

with $m\ll 0$ and $B=\left[(m,n)\times [0,T]\right]$ where ${\left|\right|u\left|\right|}_{\mathrm{H}}^2={\int }_B{u}^{2}dxdydt,$ then

$$\begin{array}{c}\hfill {\mathcal{L}}_x^{-1}{\mathcal{L}}_t^{-1}\left\{{\mathcal{L}}_x{\mathcal{L}}_t\left\{u(x,y,t)\right\}\right\}<\infty .\end{array}$$

Now to show $\mathcal{T}$ be semi-continuous [23], we take the following assumption as

Assumption 1.

$H_1$ For $\sigma >0$, ∃ a constant $\beta >0$, and $\forall u,v\in \mathrm{H}$ with $k\left|\right|u+v\left|\right|\le \sigma$, we obtain $\parallel \mathcal{T}\left(u\right)-\mathcal{T}\left(v\right)\parallel \le \beta \left|\right|u-v\left|\right|,\forall u,vs.\in \mathrm{H}$.

Theorem 3

([24] (Convergence condition)). The considered problem is tested in Equation (1) without initial and boundary conditions converging to a particular solution.

Using the Assumption 1 for operator $\mathcal{T}\left(\varphi \right)$ in Equation (1), ∋

$$\begin{array}{ccc}\hfill \mathcal{T}\left(u\right)-\mathcal{T}\left(v\right)& =& (D_x^2+D_y^2)u^2+k{u}^a-kr{u}^{(a+b)}-[(D_x^2+D_y^2)v^2+k{v}^a-kr{v}^{(a+b)}],\hfill \\ & =& D_x^2(u^2-v^2)+D_y^2(u^2-v^2)+k(u^a-v^a)-kr(u^{a+b}-v^{a+b})\hfill \end{array}$$

On taking the norm

$$\parallel \mathcal{T}\left(u\right)-\mathcal{T}\left(v\right)\parallel \le \parallel D_x^2(u+v)(u-v)\parallel +\parallel D_y^2(u+v)(u-v)\parallel +k\parallel u^a-v^a\parallel -kr\parallel u^{a+b}-v^{a+b}\parallel ,$$

by using the conditions on the operators $D_x^2$ and $D_y^2$ in $\mathrm{H}$, ∋ for ${\eta }_1,{\eta }_2>0$, and if $a=b=1$ we can define

$$D_x^2(u+v)(u-v)\le {\eta }_1\left|\right|u-v\left|\right|,D_y^2(u+v)(u-v)\le {\eta }_2\left|\right|u-v\left|\right|.$$

Therefore,

$$\parallel \mathcal{T}\left(u\right)-\mathcal{T}\left(v\right)\parallel \le {\eta }_1\left|\right|u-v\left|\right|+{\eta }_2\left|\right|u-v\left|\right|+k\parallel u-v\parallel -\sigma r\parallel u-v\parallel ,$$

$$\parallel \mathcal{T}\left(u\right)-\mathcal{T}\left(v\right)\parallel \le ({\eta }_1+{\eta }_2+k-\sigma r)\parallel u-v\parallel ,$$

taking $\beta =({\eta }_1+{\eta }_2+k-\sigma r)>0$, we can write

$$\parallel \mathcal{T}\left(u\right)-\mathcal{T}\left(v\right)\parallel \le \beta \parallel u-v\parallel .$$

Hence assumption 1 is satisfied. Thus, the proposed method is convergent.

5. Applications

Here, we present an algorithm to study the analytical solution of a fractional-order Biological population model.

Double Laplace Adomian Decomposition Method

Taking double Laplace transform (DLT) on both side of the model (1) as,

$$\begin{array}{cc}\hfill & {\mathcal{L}}_x{\mathcal{L}}_t\left[D_t^{\delta }u(x,y,t)\right]={\mathcal{L}}_x{\mathcal{L}}_t[(D_x^2+D_y^2)u^2+kr{u}^{(a+b)}+k{u}^a],\hfill \end{array}$$

(3)

applying DLT and using initial value

$$\begin{array}{cc}& \frac{s}{s+\delta (1-s)}[{\mathcal{L}}_x{\mathcal{L}}_{t}u(x,y,t)-{\mathcal{L}}_{x}u(x,y,0)]={\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)u^2\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(u^{(a+b)}\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(u^a\right),\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left[u(x,y,t)\right]=u(x,y,0)+\frac{s+\delta (1-s)}{s}\left({\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)u^2\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(u^{(a+b)}\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(u^a\right)\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\mathcal{L}}_x{f}_0+\frac{s+\delta (1-s)}{s}\left({\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)u^2\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(u^{(a+b)}\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(u^a\right)\right)\hfill \end{array}$$

(4)

Since the considered model (1), contain non-linear terms $u^2,u^{(a+b)},(a+b)\ne 0,u^a$

$$\left\{\begin{array}{c}u(x,y,t)=\sum _{m=0}^{\infty }{u}_m(x,y,t),u^2(x,y,t)=\sum _{m=0}^{\infty }{P}_m,\hfill \\ u^{(a+b)}(x,y,t)=\sum _{m=0}^{\infty }{Q}_m.\hfill \end{array}\right.$$

(5)

After decomposing the non-linear terms, (4) can be written as

$$\begin{array}{cc}\hfill & {\mathcal{L}}_x{\mathcal{L}}_t\left[\sum _{m=0}^{\infty }{u}_m(x,y,t)\right]=f_0(x,y)+\frac{s+\delta (1-s)}{s}\left({\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)\sum _{m=0}^{\infty }{P}_m\right]\right.\hfill \\ & \left.+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(\sum _{m=0}^{\infty }{Q}_m\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(\sum _{m=0}^{\infty }{T}_m\right)\right),\hfill \end{array}$$

(6)

comparing terms both side, we obtain

$$\left\{\begin{array}{cc}\hfill & {\mathcal{L}}_x{\mathcal{L}}_t\left({\mathrm{Q}}_0(x,y,0)\right)={\mathcal{L}}_x{f}_0(x,y)\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left(u_1(x,y,t)\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)P_0\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(Q_0\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(T_0\right),\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left(u_2(x,y,t)\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)P_1\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(Q_1\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(T_1\right),\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left(u_3(x,y,t)\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)P_2\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(Q_2\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(T_2\right)\hfill \\ & ⋮\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left(u_{m+1}(x,y,t)\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t\left[(D_x^2+D_y^2)P_m\right]+kr{\mathcal{L}}_x{\mathcal{L}}_t\left(Q_m\right)+k{\mathcal{L}}_x{\mathcal{L}}_t\left(T_m\right).\hfill \end{array}\right.$$

(7)

Taking inverse DLT of (7), both side we obtain

$$\left\{\begin{array}{cc}\hfill & {\mathrm{Q}}_0(x,y,t)=f_0(x,y),\hfill \\ & u_1(x,y,t)=\left((D_x^2+D_y^2)P_0+kr{Q}_0+T_0\right)\left(1+\xi (t-1)\right)\hfill \\ & u_2(x,y,t)=\left((D_x^2+D_y^2)P_1+kr{Q}_1+T_1\right)\left(1+{\xi }^2(t-1)\right)\hfill \\ & u_3(x,y,t)=\left((D_x^2+D_y^2)P_2+kr{Q}_2+T_2\right)\left(1+{\xi }^3(t-1)\right)\hfill \\ & ⋮\hfill \\ & u_{m+1}(x,y,t)=\left((D_x^2+D_y^2)P_m+kr{Q}_m+T_m\right)\left(1+{\xi }^{m+1}(t-1)\right).\hfill \end{array}\right.$$

(8)

After simplification we obtain the values of ${\mathrm{Q}}_0,u_1,u_2,u_3$. Similarly in (8), we can obtain

$$\begin{array}{cc}& u(x,y,t)=u_0+u_1+u_2+u_3+\dots ,\hfill \end{array}$$

6. Iterative Examples

In this section, we use some examples for the validation of the proposed method.

Example 1

([15]). Our first example is

$$D_t^{\delta }u(x,y,t)=\left(D_x^2+D_y^2\right)u(x,y,t)+ku(x,y,t)$$

(9)

with initial condition $u(x,y,0)=\surd xy$

Applying DLT on both side, we obtain

$$\begin{array}{cc}& {\mathcal{L}}_x{\mathcal{L}}_t\left[D_t^{\delta }u(x,y,t)\right]={\mathcal{L}}_x{\mathcal{L}}_t[\left(D_x^2+D_y^2\right)u(x,y,t)+ku(x,y,t)]\hfill \end{array}$$

after calculation we obtain

$$\begin{array}{cc}& {\mathcal{L}}_x{\mathcal{L}}_t\left(\delta (x,y,t)\right)={\mathcal{L}}_x{f}_0(x,y)+\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t[(D_x^2+D_y^2)u+ku],\hfill \end{array}$$

Comparing both sides

$$\begin{array}{cc}& {\mathcal{L}}_x{\mathcal{L}}_t\left({\mathrm{Q}}_0(x,y,0)\right)={\mathcal{L}}_x{f}_0(x,y){\mathcal{L}}_x{\mathcal{L}}_t\left(u_1\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t[(D_x^2+D_y^2){\mathrm{Q}}_0+k{\mathrm{Q}}_0],\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left(u_2\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t[(D_x^2+D_y^2)u_1+k{u}_1],\hfill \\ & {\mathcal{L}}_x{\mathcal{L}}_t\left(u_3\right)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t[(D_x^2+D_y^2)u_2+k{u}_2]\hfill \end{array}$$

After applying inverse DLT, we obtain

$$\begin{array}{cc}\hfill & u(x,y,0)=\surd xy,\hfill \\ & u_1=k\surd xy(1+\xi (t-1))\hfill \\ & u_2=k^2\surd xy(1+{\xi }^2(t-1))\hfill \\ & u_3=k^3\surd xy(1+{\xi }^3(t-1))\hfill \end{array}$$

Thus, the series solution of (9),

$$\begin{array}{cc}\hfill & u(x,y,t)={\mathrm{Q}}_0+u_1+u_2+u_3+\dots +\hfill \\ & u(x,y,t)=\surd xy+k\surd xy(1+\xi (t-1))+k^2\surd xy(1+{\xi }^2(t-1))+\dots \hfill \end{array}$$

(10)

Example 2.

Let us suppose the following fractional order biological model

$$\begin{array}{cc}\hfill & D_t^{\delta }(x,y,t)=\left(D_x^2+D_y^2\right)u(x,y,t)-r{u}^2+u.\hfill \end{array}$$

(11)

with given subsidiary condition [15] $u(x,y,0)=exp\left(\frac{1}{2}\surd \left(\frac{r}{2}\right)(x+y)\right)$. Using DLT

$$\begin{array}{cc}\hfill & {\mathcal{L}}_x{\mathcal{L}}_t\left(D_t^{\delta }u(x,y,t)\right)={\mathcal{L}}_x{\mathcal{L}}_t[\left(D_x^2+D_y^2\right)u-r{u}^2+u],\hfill \\ & u(x,y,t)=\frac{s+\delta (1-s)}{s}{\mathcal{L}}_x{\mathcal{L}}_t[\left(D_x^2+D_y^2\right)u-r{u}^2+u],\hfill \\ & {\mathrm{Q}}_0=exp\left(\frac{1}{2}\surd \left(\frac{r}{2}\right)(x+y)\right),\hfill \\ & u_1=exp\left(\frac{1}{2}\surd \left(\frac{r}{2}\right)(x+y)\right)(1+\xi (t-1)),\hfill \\ & u_2=exp\left(\frac{1}{2}\surd \left(\frac{r}{2}\right)(x+y)\right)(1+{\xi }^2(t-1)),\hfill \\ & u_3=exp\left(\frac{1}{2}\surd \left(\frac{r}{2}\right)(x+y)\right)(1+{\xi }^3(t-1)).\hfill \end{array}$$

Thus, the series solution of (11),

$$\begin{array}{cc}\hfill & u(x,y,t)=u_0+u_1+u_2+u_3+\dots \hfill \end{array}$$

(12)

6.1. Numerical Plots and Comparison Tables

In this section, we present the numerical illustrations of the approximate solutions Equations (10) and (12) of the iterative examples considered on proposed model. For the graphical representations we consider the parameters as $k=3,r=50$ and the fractional orders are considered to be $\xi =1,0.9,0.8$ and 0.7. In each of the following figures the fractional orders are considered as Figure 1a–d. Figure 1 and Figure 2 shows the dynamics of the analytical approximate solution Equation (10). Similarly, Figure 3 depicts the dynamics in the behaviors of the solutions Equations (10) and (12) with $x=3$ and $x=4$ respectively and varying $y,t$. Furthermore, Figure 4 represents the changes that occur in the behaviors of the solutions Equation (12) respectively with $t=10$ and varying the variables $x,y$. Here we see that when fractional order is large then the population density is also large, on the other hand population density decreases with the decrease in the fractional order.

Figure 5 shows the dynamical behavior of the solution Equation (12) with fixed $y=4$ and $y=2$, respectively, and varying the other variables $x,t$. Furthermore, Figure 3 and Figure 6 depicts the dynamics in the behaviors of the solutions Equations (10) and (12) with $x=3$ and $x=4$, respectively, and varying $y,t$. We observed that the considered scheme is rapidly convergent and is highly accurate.

It is observed that at lower fractional orders, the population density decreases, as can be seen in the figure. The considered operator shows an astonishing impact on the dynamics of the considered population model. The simulations show that the double Laplace transform is a powerful method that can be used to study such complexities with ease, without perturbing or dealing with the long-lasting polynomials.

Next, we present the comparison table between the exact and the approximate series solutions of the considered iterative examples in the sense of the $CF$ operator. We considered $k=2$ for

Table 1. Comparison table between the exact [15] and the approximate series solution Equation (10).

(x, y, t)	Exact	Approximate	Error
(1, 1, 0)	1	1.0001	$9.8000\times {10}^{-5}$
(1, 1, 1)	31	31	0
(2, 1, 2)	86.2670	86.2670	$1.3859\times {10}^{-4}$
(2, 2, 2)	122	121.9998	$1.9600\times {10}^{-4}$
(3, 2, 1)	75.9342	75.9342	0
(5, 2, 0)	3.1623	3.1626	$3.0990\times {10}^{-4}$
(3, 5, 1)	120.0625	120.0625	0
(7, 5, 3)	538.3633	538.3621	0.0012
(1, 3, 0)	1.7321	1.7322	$1.6974\times {10}^{-4}$
(10, 10, 1)	310	310	0
(10, 10, 0)	10	10.0010	$9.8000\times {10}^{-4}$
(9, 2, 3)	386.0803	386.0795	$8.3156\times {10}^{-4}$
(4, 4, 1)	124	124	0
(4, 8, 4)	684.4794	684.4777	0.0017
(8, 9, 5)	$1.2813\times {10}^3$	$1.2813\times {10}^3$	0.0033
(9, 8, 1)	263.0437	263.0437	0
(0.5, 1, 1)	21.9203	21.9203	0
(0.2, 1, 0)	0.4472	0.4473	$4.3827\times {10}^{-5}$
(0.1, 0.1, 0)	0.1000	0.1000	$9.8000\times {10}^{-6}$

, fractional-order as 1 and different values of variables $x,y,t$. Similarly, we supposed $r=50$ for

Table 2. Comparison table between the exact [15] and the approximate series solution Equation (12).

(x, y, t)	Exact	Approximate	Error
(−10,−10, 0)	$3.7201\times {10}^{-44}$	$3.7201\times {10}^{-44}$	$3.7201\times {10}^{-52}$
(−10, −10, 1)	$1.8600\times {10}^{-43}$	$1.8600\times {10}^{-43}$	0
(−6, −6, 0)	$8.7565\times {10}^{-27}$	$8.7565\times {10}^{-27}$	$8.7565\times {10}^{-35}$
(−6, −6, 1)	$4.3783\times {10}^{-26}$	$4.3783\times {10}^{-26}$	0
(0, −2, 2)	$3.1780\times {10}^{-4}$	$3.1780\times {10}^{-4}$	$1.3620\times {10}^{-12}$
(0, 2, 2)	$1.9824\times {10}^5$	$1.9824\times {10}^5$	$2.2026\times {10}^{-4}$
(−10, 0, 10)	$7.9079\times {10}^{-21}$	$7.9079\times {10}^{-21}$	$1.7359\times {10}^{-29}$
(−10, 1, 10)	$1.1736\times {10}^{-18}$	$1.1736\times {10}^{-18}$	$2.5763\times {10}^{-27}$
(−8, 0, 8)	$1.4020\times {10}^{-16}$	$1.4020\times {10}^{-16}$	$2.9738\times {10}^{-25}$
(8, 1, 8)	$2.0807\times {10}^{-14}$	$2.0807\times {10}^{-14}$	$4.4136\times {10}^{-23}$
(−2, 0, 2)	$4.0860\times {10}^{-4}$	$4.0860\times {10}^{-4}$	$4.4948\times {10}^{-6}$
(2, 0, 2)	0.0606	0.0600	$6.6709\times {10}^{-4}$
(−1, −1, 0)	$4.5400\times {10}^{-5}$	$4.5400\times {10}^{-5}$	$4.4948\times {10}^{-6}$
(1, 1, 1)	$1.1013\times {10}^5$	$1.1013\times {10}^5$	0
(2, −2, 0)	1	1.0990	0.0990
(−3, −1, 2)	$1.8550\times {10}^{-8}$	$1.8550\times {10}^{-8}$	$2.0406\times {10}^{-10}$
(0, −2, 1)	$2.2700\times {10}^{-4}$	$2.2700\times {10}^{-4}$	0
(2, 0, 0)	$4.5400\times {10}^{-5}$	$4.5400\times {10}^{-5}$	$4.4948\times {10}^{-6}$

, fractional-order as 1 and various values of variables $x,y,t$. A good agreement is observed between the exact and approximate solutions which can be seen in the table.

7. Conclusions

In this manuscript, we proposed a fractional-order biological population system. Using fixed point theorems, we effectively developed proofs for the existence uniqueness of the considered model results. We obtained the numerical solutions of the fractional-order model (1) using DLADM. Furthermore, the convergence analysis of the proposed model with the considered method has been presented. From the numerical simulations, it is observed that when the fractional order is large then the population density is also large; on the other hand, population density decreases with the decrease in the fractional order. The error analysis shows that the technique approximate solution converges to the exact solution very quickly. As a result, we observed that the discussed methodology has significant advantages for obtaining analytical solutions for the FNDEs when compared to other methods. From the obtained results, we concluded that the DLADM is capable of reducing the volume of computational work as compared to other methods. The considered scheme does not need much analytical study and converges easily to the exact solution. In the future, we will suggest modifying the considered scheme of DLADM for the investigation of nonlinear partial differential equations and some advanced problems in fluid dynamics and elasticity dealing with integrals will be investigated in subsequent papers. Finally, we state that the considered method is indeed trustworthy and applicable to all nonlinear PDEs subject to initial condition.