# Approximating Solution of Fabrizio-Caputo Volterra’s Model for Population Growth in a Closed System by Homotopy Analysis Method

Authors: Tahereh Bashiri; S. Mansour Vaezpour; Juan J. Nieto
Journal: Journal of Function Spaces
 (2018)

---

## Abstract

Volterra’s model for population growth in a closed system consists in an integral term to indicate accumulated toxicity besides the usual terms of the logistic equation. Scudo in 1971 suggested the Volterra model for a populationu(t) of identical individuals to show crowding and sensitivity to “total metabolism”: du/dt=au(t)-bu2(t)-cu(t)∫0tu(s)ds. In this paper our target is studying the existence and uniqueness as well as approximating the following Caputo-Fabrizio Volterra’s model for population growth in a closed system: CFDαu(t)=au(t)-bu2(t)-cu(t)∫0tu(s)ds, α∈[0,1], subject to the initial condition u(0)=0. The mechanism for approximating the solution is Homotopy Analysis Method which is a semianalytical technique to solve nonlinear ordinary and partial differential equations. Furthermore, we use the same method to analyze a similar closed system by considering classical Caputo’s fractional derivative. Comparison between the results for these two factional derivatives is also included.


---

## Body

## 1. Introduction

Malthus was the first economist to propose a systematic theory of population [1] where he gathered experimental data to support his thesis. He proposed the principle that human populations grow exponentially. Consequently, a nonlinear growth equation was introduced into population dynamics by Verhulst [2] to solve the unbounded growth in human population proposed by Malthus. Verhulst introduced the nonlinear term into the rate equation and reached what afterwards became nominated as the logistic equation:(1)dutdt=kut1-ut.Volterra’s model for population growth in a closed system includes an integral term to indicate accumulated toxicity in addition to the usual terms of the logistic equation. Scudo, in 1971, indicates that Volterra proposed his following model for a populationu(t) of identical individuals which exhibits crowding and sensitivity to “total metabolism”:(2)dudt=aut-bu2t-cut∫0tusds.If the integral term on the right is removed the famous logistic equation with birth ratea and crowding coefficient b appears. The last term contains the integral that indicates the “total metabolism” or total amount of toxins produced since time zero.The individual death rate is corresponding to this integral, so the population death rate by virtue of toxicity must include a factoru. The existence of the toxic term as a result of the system being closed always causes the population level to fall to zero in the long run, as will be seen shortly. The relative size of the sensitivity to toxins, denoted by c, determines the manner in which the population thrives before its decay.The tool for describing the behaviour of the equations such as logistic equation must be nonlocal differential equations in time. With this purpose, in the last decades, the Fractional Calculus (FC) allows investigating the nonlocal response of multiple phenomena [3–9]. Fractional derivatives are memory operators which usually represent dissipative effects or damage.Some fundamental definitions in the context of FC are Erdelyi-Kober, Riesz, Riemann-Liouville, Hadamard, Graunwald-Letnikov, Weyl, Jumarie, or Caputo (see, e.g., [10–15] and references therein). The Riemann-Liouville definition brings about physically unacceptable initial conditions (fractional order initial conditions) [12]. In the Caputo concept, the initial conditions are expressed in terms of integer-order derivatives, so it has physical meaning [13]. These definitions have the disadvantage that their kernel has singularity; this kernel includes memory effects and therefore both definitions cannot precisely describe the full effect of the memory [16]. Due to this problem, Caputo and Fabrizio in [17] introduced a new definition of fractional derivative without singular kernel, the Caputo-Fabrizio (CF) fractional derivative; this derivative possesses very interesting properties which are reviewed in detail in [18]. Caputo and Fabrizio demonstrated that the new derivative involves further properties in comparison with the old version. They indicated, for instance, that it can interpret substance heterogeneities and configurations with different scales, which apparently cannot be investigated with the prominent local theories and also the old versions of fractional derivative. Another utilization is scrutinizing the macroscopic behaviours of some materials, identified with nonlocal communications between atoms. Other applications of the CF fractional derivative can be achieved, for example, in [19–21].The main aims of this paper are studying the existence and uniqueness as well as approximating the Caputo-Fabrizio fractional population growth:(3)CFDαut=aut-bu2t-cut∫0tusdsα∈0,1,subject to the initial condition(4)u0=0.Moreover, we shall compare the results with those obtained for the same problem by considering the classical Caputo fractional derivative in (3).The auxiliary method used here for approximating the solution is a semianalytical technique known as Homotopy Analysis Method to solve nonlinear ordinary and partial differential equations.Next, we introduce some basic definitions and notations that shall be used in next sections.Definition 1 (see [22]).
For at leastn-times continuously differentiable function g:[0,∞)→R, the Caputo fractional derivative of order σ>0 is defined as(5)cDσht=1Γn-σ∫0tt-wn-σ-1gnwdw,where n=[σ]+1.By changing the kernel(t-s)-a by the function exp(-(α/1-α)(t-s)) and 1/Γ(1-α) by M(α)/1-α, one obtains the new Caputo-Fabrizio fractional derivative of order 0<α<1, which has been recently introduced by Caputo and Fabrizio in [17]. More precisely, they introduced the following concept.Definition 2 (see [17]).
For0<α<1, h∈H1(0,b), b>0. The Caputo-Fabrizio fractional derivative is defined by (6)CFDαht=Mα1-α∫0te-α/1-αt-sh′sds,where H1(0,b) is the usual Sobolev space, the set of function in L2(0,b) such that the derivative, in the sense of distributions, is in L2(0,b). M(α) is a normalization constant (M(1)=M(0)=1). For simplicity, we shall consider M(α)=1.Lemma 3 (see [17]).
Let0<α<1, and let us consider the following fractional differential equation:(7)CFDαft=ut,t≥0.One can get(8)ft=1-αut-u0+α∫0tusds+f0,t≥0.


## 2. Existence and Uniqueness

In this section, we aim to show the existence and uniqueness of solution of the fractional Volterra population growth model (3). By Lemma 3 and using initial conditions we obtain the following result.Theorem 4.
If(1-α)L+αTL<1, for t∈J=[0,T] and some L>0, then (3) has a unique solution which is zero solution.Proof.
First, we have(9)ut=1-αaut-bu2t-cut∫0tusds+α∫0taus-bu2s-cus∫0suxdxds.Let(10)Θt,u=aut-bu2t-cut∫0tusds.We have to find a positive real number L such that(11)Θt,u-Θt,v≤Lu-v.We obtain(12)Θt,u-Θt,v≤aut-vt+bv2t-u2t+cvt∫0tvsds-ut∫0tusds≤au-v+bu+vu-v+cvt∫0tvsds±vt∫0tusds-ut∫0tusds≤u-v+bu+vu-v+cv∫0tu-vds+cu-v∫0tuds.Due to the assumption that u and v are bounded, and t∈[0,T], there exists a positive constant η>0 such that v≤η and u≤η. Thus, we get(13)Θt,u-Θt,v≤1+2ηb+2cηTu-v.This shows the Lipschitz condition for Θ by taking L=1+2ηb+2cηT. Regarding (9) we can get(14)u0t=t,unt=1-αΘt,un-1+α∫0tΘt,un-1ds.Now, we set Fn(t)=un(t)-un-1(t). So we have(15)Fnt=1-αΘt,un-1-Θt,un-2+α∫0tΘs,un-1-Θs,un-2ds.Taking the norm of the latter equation gives(16)Fnt≤1-αΘt,un-1-Θt,un-2+α∫0tΘs,un-1-Θs,un-2ds≤1-αLun-1-un-2+αL∫0tun-1-un-2ds≤1-αLFn-1+αL∫0tFn-1ds.Then, by using the recursive principle it yields(17)Fnt≤1-αLn+αLtnu0t,which proves that the solution exists and is continuous. We show that u¯(t) defined as(18)u¯t=limn→∞unt,which is a solution for (3). Let Zn(t)=u¯(t)-un(t). We know Zn should tend to zero as n→∞. We can see(19)u¯t-un+1t≤1-αΘt,u¯-Θt,un+α∫0tΘs,u¯-Θs,unds≤1-αLu¯-un+αL∫0tu¯-unds≤1-αLu¯-un+αL∫0tu¯-unds≤1-αLZn+αL∫0tZnds. So, the right hand side gives(20)limn→∞un=u¯.We can take u¯ as a solution of (3) that is continuous. Furthermore, applying the Lipschitz condition for Θ, we have(21)u¯t-1-αΘt,u¯-α∫0tΘs,u¯ds=Zn+1-αΘt,u¯n+α∫0tΘs,u¯nds-1-αΘt,u¯-α∫0tΘs,u¯ds.So(22)u¯t-1-αΘt,u¯-α∫0tΘs,u¯ds≤Zn+1-αΘt,u¯n-Θt,u¯+α∫0tΘs,u¯n-Θs,u¯ds≤Zn+1-αLu¯n-u¯+αtLu¯n-u¯.Passing to the limit when n→∞ and regarding the initial condition, we obtain(23)u¯t=1-αΘt,u¯+α∫0tΘs,u¯ds.For uniqueness we consider u and v to be two different solutions of (3); then, the Lipschitz condition for Θ gives the following result:(24)u-v≤1-αLu-v+αTLu-v,rearranged to be(25)u-v1-1-αL-αTL≤0.Then, u-v=0 if(26)1-αL+αTL<1.On the other hand, the constant function equal to zero is solution of the equation. Thus, the proof is completed.


## 3. Homotopy Analysis Method

The Homotopy Analysis Method operates the theory of the homotopy from topology to generate a convergent series solution for nonlinear systems.The HAM was first constructed in 1992 by Liao in his Ph.D. Thesis [23] and later altered [24] in 1999 by proposing a nonzero auxiliary parameter, as the convergence-control parameter, h, to construct a homotopy on a differential system in general form [25]. The convergence-control parameter is a nonphysical variable that provides a simple way to certify convergence of a solution series. It is a series expansion method which is not directly dependent on small or large physical parameters. It also provides extreme flexibility to choose the basis functions of the desired solution and the corresponding auxiliary linear operator of the homotopy.Now, by using HAM to approximate the solution of the model (3), in view of properties of Caputo-Fabrizio fractional integrals, we assume that u(t) can be expressed by the functions(27)tαm±neαp±qt:m,n,p,q=0,1,2,….It is essential to know that we have a great freedom to choose auxiliary parameters in HAM. From base functions denoted by (27) and the initial condition given in (4), it is convenient to choose(28)u0t=t,as the initial approximation of u(t), as well as(29)Lφ=CFDαφ,as the auxiliary linear operator.Leth denote a nonzero auxiliary parameter. We prepare the HAM deformation equation as(30)1-qLΨt,q-u0t=hqNΨt,q,subject to the initial condition(31)Ψ0,q=0,where q∈[0,1] is an embedding parameter and NΨt,q is a nonlinear operator, given by(32)NΨt,q=CFDαΨt,q-aΨt,q+bΨ2t,q+cΨt,q∫0tΨs,qds.Evidently, when q=0, the solution of (30) is(33)Ψt,0=u0t.When q=1, (30) is the same as the original equation (3), provided(34)Ψt,1=ut.Thus, Ψ(t,q) alters from the initial approximation u0(t) to the exact solution u(t). Expanding Taylor’s series with respect to q, we have(35)Ψt,q=Ψt,0+∑k=1∞uktqk,where(36)ukt=1k!∂kΨt,q∂qkq=0.For briefness, one can define the vector(37)u→k=u0,u1,…,uk.Differentiating the HAM deformation equation (30) k times with respect to q, then setting q=0, and finally dividing them by k!, we attain the kth-order deformation equation(38)Lukt-χkuk-1t=hHtRku→k-1,subject to the initial condition(39)uk0=0,where(40)Rku→k-1=CFDαuk-1-auk-1+b∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-ids,χk=0k=1,1k>1,and H(t)=1.Next, our target is to reach a recursive equation to be implemented, for example, in MATLAB, in order to provide an approximation of the solution.From (38) we have(41)CFDαukt-χkuk-1t=hCFDαuk-1+b∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-ids-auk-1,which implies that(42)CFDαukt-χkCFDαuk-1t=hCFDαuk-1+hb∑i=0k-1uiuk-1-i+hc∑i=0k-1ui∫0tuk-1-ids-hauk-1;then,(43)CFDαukt=h+χkCFDαuk-1t+hb∑i=0k-1uiuk-1-i+hc∑i=0k-1ui∫0tuk-1-ids-hauk-1.Now, by Lemma 3 we have(44)ukt=1-αh+χkCFDαuk-1t+hb∑i=0k-1uiuk-1-i+hc∑i=0k-1ui∫0tuk-1-idx-hauk-1+α∫0th+χkCFDαuk-1s+hb∑i=0k-1uiuk-1-i+hc∑i=0k-1ui∫0suk-1-idx-hauk-1ds.Hence, we can conclude(45)ukt=1-αhb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-idx-auk-1+αh∫0tb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-idx-auk-1ds+h+χk1-αCFDαuk-1+α∫0tCFDαuk-1sds.By definition of the Caputo-Fabrizio fractional derivative, we get(46)ukt=1-αhb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-idx-auk-1+αh∫0tb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0suk-1-idx-auk-1ds+h+χke-α/1-αt∫0teα/1-αsuk-1′sds+α1-α∫0te-α/1-αs∫0seα/1-αwuk-1′wdwds.Now, using integration by parts as well as the initial condition we obtain(47)ukt=1-αhb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-idx-auk-1+αh∫0tb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0suk-1-idx-auk-1ds+h+χke-α/1-αteα/1-αtuk-1t-α1-α∫0teα/1-αsuk-1sds+α1-α∫0te-α/1-αseα/1-αsuk-1s-α1-α∫0weα/1-αwuk-1wdwds.So,(48)ukt=1-αhb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0tuk-1-idx-auk-1+αh∫0tb∑i=0k-1uiuk-1-i+c∑i=0k-1ui∫0suk-1-idx-auk-1ds+h+χkuk-1t-α1-αe-α/1-αt∫0teα/1-αsuk-1sds+α1-α∫0tuk-1s-α1-αe-α/1-αs∫0seα/1-αwuk-1wdwds.Then, we realize(49)ukt=1-αhb∑i=0k-1uiuk-1-i+1-αhc∑i=0k-1ui∫0tuk-1-idx-1-αhauk-1+αhb∫0t∑i=0k-1uiuk-1-ids+αhc∫0t∑i=0k-1ui∫0suk-1-idxds-αha∫0tuk-1ds+h+χkuk-1t-h+χkα1-αe-α/1-αt∫0teα/1-αsuk-1sds+h+χkα1-α∫0tuk-1sds-h+χkα1-α2∫0te-α/1-αs∫0seα/1-αwuk-1wdwds.Finally, we have the following equation:(50)ukt=1-αhb∑i=0k-1uiuk-1-i+h+χk-1-αhauk-1+1-αhc∑i=0k-1ui∫0tuk-1-idx+αhc∫0t∑i=0k-1ui∫0suk-1-idxds+αhb∫0t∑i=0k-1uiuk-1-ids+h+χkα1-α-αha∫0tuk-1ds-h+χkα1-αe-α/1-αt∫0teα/1-αsuk-1sds-h+χkα1-α2∫0te-α/1-αs∫0seα/1-αwuk-1wdwds.As suggested by Liao [25], the appropriate region for h is a horizontal line segment. We can investigate the influence of h on the convergence of u, by plotting the curve of it versus h, as shown in Figure 1. It seems that the best choice can be in the region [-1,1]. We can also see the results of the 8th order analytical approximations for u(t) given by HMA for Fabrizio-Caputo differential equation in Figure 2.Figure 1
Theh-curves obtained from the 6-order HAM approximate solution for Fabrizio-Caputo logistic equation.Figure 2
The 8th order analytical approximations foru(t) given by HMA for Fabrizio-Caputo logistic differential equation.Xu [26] has also studied the Caputo form of this equation for nonzero initial conditions. The result for this equation with zero initial conditions is as follows:(51)ukt=h+χkuk-1-h+χk∑i=0n-1uk-1i0tii!+hJαb∑i=0k-1uiuk-1-i-auk-1+c∑i=0k-1ui∫0tuk-1-ids.As shown in Figure 3, it seems that the best choice for h can be in the region [-1,1]. The 8th order analytical approximations for u(t) given by HMA for the differential equation by considering the classical Caputo fractional derivative is shown in Figure 4. This approximation is for a=b=c=10 and α=0.5. We also get here u0(t)=t.Figure 3
Theh-curves obtained from the 6th-order HAM approximate solution for Caputo logistic equation.Figure 4
The 8th order analytical approximations foru(t) given by HMA by considering the classical Caputo fractional derivative.As a result, we can see that the solutions are similar with more oscillation in the case of Caputo-Fabrizio type of the population growth.


## 4. MATLAB Codes

Here, we prepare codes written by MATLAB related to plotting the figures you can see as the results of approximation by HAM.MATLAB code for Figure1:clc; clear all; close all;
symstuxshzq=0.5;a=10;b=10;c=10;forj=2:6ifj==2z(j)=0elsez(j)=1endendforj=2:6u(1)=t;g=0fori=1:(j-1)f=u(i)∗u(j-i)g=g+f;endl=0fori=1:(j-1)m=u(i)∗int(u(j-i),0,t)l=l+m;endu(j)=((h+z(j)-((1-q)∗h∗a))∗u(j-1)) + (((1-q)∗h∗a)∗l) + ((q∗h∗c)∗int(l,0,t)) + (((((h+z(j))∗q)/(1-q))-(q∗h∗a))∗int(u(j-1),0,t)) + ((q∗h∗b)∗int(g,0,t)) − ((((h+z(j))∗q)/(1-q))∗exp(-(q/(1-q))∗t)∗int(exp((q/(1-q))∗s)∗u(j-1),s,0,t)) + (((1-q)∗h∗b)∗g) − ((h+z(j))∗((q/(1-q))2)∗int(exp(-(q/(1-q))∗s)∗(int(exp((q/(1-q))∗x)∗u(j-1),x,0,s)),s,0,t))endusum=tforj=2:6usum=usum+u(j)endsymspop(t)=usumo(h)=p(1)ezplot(o(h),[-4,4])MATLAB code for Figure2:clc; clear all; close all;
symstuxshzq=0.5;a=10;b=10;c=10;forj=2:8ifj==2z(j)=0elsez(j)=1endendforj=2:8u(1)=t;g=0fori=1:(j-1)f=u(i)∗u(j-i)g=g+f;endl=0fori=1:(j-1)m=u(i)∗int(u(j-i),0,t)l=l+m;endu(j)=((h+z(j)-((1-q)∗h∗a))∗u(j-1)) + (((1-q)∗h∗c)∗l) + ((q∗h∗c)∗int(l,0,t)) + (((((h+z(j))∗q)/(1-q)) − (q∗h∗a))∗int(u(j-1),0,t)) + ((q∗h∗b)∗int(g,0,t)) − ((((h+z(j))∗q)/(1-q))∗exp(-(q/(1-q))∗t)∗int(exp((q/(1-q))∗s)∗u(j-1),s,0,t)) + (((1-q)∗h∗b)∗g) − ((h+z(j))∗((q/(1-q))2)∗int(exp(-(q/(1-q))∗s)∗(int(exp((q/(1-q))∗x)∗u(j-1),x,0,s)),s,0,t))endusum=tforj=2:8usum=usum+u(j)endsymspp(h)=usum;figure;hold on;p1=ezplot(p(-1),[0,1])p2=ezplot(p(-0.5),[0,1])p3=ezplot(p(0.5),[0,1])p4=ezplot(p(1),[0,1])set (p1, “color”, “r”, “linestyle”, “−−”)set (p2, “color”, “b”, “linestyle”, “−.”)set (p3, “color”, “k”, “linestyle”, “:”)set (p4, “color”, “g”, “linestyle”, “−”)legend (“h=-1”, “h=-0.5”, “h=0.5”, “h=1”)MATLAB code for Figure3:clc; clear all; close all;
symstuxshznq=0.5; a=10; n=1; b=10; c=10;u(1)=t;forj=2:4ifj==2z(j)=0elsez(j)=1endendforj=2:4g=0fori=1:(j-1)f=u(i)∗u(j-i)g=g+fendl=0fori=1:(j-1)m=u(i)∗int(u(j-i),0,t)l=l+mendw=0fori=1:nr(t)=diff(u(j-1),i-1)y(i)=(r(0)∗(t(i-1)))/(factorial(i-1))w=w+y(i)endu(j)=((h+z(j))∗u(j-1)) − ((h+z(j))∗w) + (((h∗b)/gamma(q))∗int((power((t-s),(q-1)))∗g,s,0,t)) − (((h∗a)/gamma(q))∗int((power((t-s),(q-1)))∗u(j-1),s,0,t)) + (((h∗c)/gamma(q))∗int((power((t-s),(q-1)))∗l,s,0,t))endv=tforj=2:4v=v+u(j)endsymspop(t)=vo(h)=p(1)ezplot(o(h),[-4,4])MATLAB code for Figure4:clc; clear all; close all;
symstuxshznq=0.5; a=10; n=1; b=10; c=10;u(1)=t;forj=2:8ifj==2z(j)=0elsez(j)=1endendforj=2:8g=0fori=1:(j-1)f=u(i)∗u(j-i)g=g+fendl=0fori=1:(j-1)m=u(i)∗int(u(j-i),0,t)l=l+mendw=0fori=1:nr(t)=diff(u(j-1),i-1)y(i)=(r(0)∗(t(i-1)))/(factorial(i-1))w=w+y(i)endu(j)=((h+z(j))∗u(j-1)) − ((h+z(j))∗w)+(((h∗b)/gamma(q))∗int((power((t-s),(q-1)))∗g,s,0,t)) − (((h∗a)/gamma(q))∗int((power((t-s),(q-1)))∗u(j-1),s,0,t)) + (((h∗c)/gamma(q))∗int((power((t-s),(q-1)))∗l,s,0,t))endv=tforj=2:8v=v+u(j)endsymspp(h)=v;figure;hold on;p1=ezplot(p(-1),[0,1])p2=ezplot(p(-0.5),[0,1])p3=ezplot(p(0.5),[0,1])p4=ezplot(p(1),[0,1])set (p1, “color”, “r”, “linestyle”, “−−”)set (p2, “color”, “b”, “linestyle”, “−.”)set (p3, “color”, “k”, “linestyle”, “:”)set (p4, “color”, “g”, “linestyle”, “−”)legend (“h=-1”, “h=-0.5”, “h=0.5”, “h=1”)



---

*Source: 3152502-2018-01-11.xml*