EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 13, No. 5, 2020, 1176-1198
ISSN 1307-5543 – www.ejpam.com
Published by New York Business Global

Special Issue Dedicated to
Professor Hari M. Srivastava

On the Occasion of his 80th Birthday

Analysis of The Picard’s Iteration Method and Stability
for Ecological Initial Value Problems of Single Species

Models with Harvesting Factor

Mohammad Hossein Rahmani Doust1,∗,V. Lokesha2, Atena Ghasemabadi3

1 Department of Mathematics, Faculty of Sciences, University of Neyshabur, Neyshabur,
Khorasen-e-Razavi, Iran
2 Department of Studies in Mathematics, V. S. K. University, Vinayaka Nagara, Ballari,
India
3 Esfarayen University of Technology, Esfayen, North Khorasan, Iran

Abstract. In the recent decades, biology and ecology area and also computer and network sciences
are marched on at a rapid pace toward perfection by help of mathematical concepts such as stability,
bifurcation, chaos and etc. Because of no existing any interspecific interaction in the single species,
one is able to see that this is the simplest model. Meanwhile by adding some assumptions, we
see that it has so many practical applications in the nature and any branch of sciences. In this
article, some dynamical models of single species are studied. First, Picard’s iteration method for
exponential growth rate is analyzed. In continuation, some logistic models for both cases without
harvesting and having harvested factor which are constant or variable are studied. Indeed, the
solution and stability of equilibria for the said models are analyzed. Finally, in the section of
simulation analysis by help of Matlab software, we give some numerical simulations to support of
our mathematical conclusions which show the stability of the equilibria for I.V.Ps. of the logistic
equation developed.

2020 Mathematics Subject Classifications: 34D20, 92D40

Key Words and Phrases: Harvesting Factor, Logistic Equation, Picard’s Iteration Method,
Single Species, Stability

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v13i5.3713

Email addresses: mh.rahmanidoust@neyshabur.ac.ir (M.H. Rahmani Doust),
v.lokesha@gmail.com (V. Lokesha), ghasemabadi.math@gmail.com (A. Ghasemabadi)

https://www.ejpam.com 1176 c© 2020 EJPAM All rights reserved.



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1177

1. Introduction

Population ecologists use methods to model the population dynamics. An accurate
model should be able to describe the changes occurring in a population and predict future
changes. Population growth is the most topic in ecology and biology area. Two simple
models of population growth use deterministic equations to describe the rate of change in
the size of a population over time. The first one, exponential growth, describes theoretical
populations that increase in numbers without any limits to their growth. While the other
one, logistic growth, introduces limits to reproductive growth that become more intense
as the population size increases. Carrying capacity of the environment is so important to
study of logistic dynamics. Neither model adequately describes natural populations, but
they provide points of comparison.

As time goes, any researcher may see the application of mathematics in interdisci-
plinary sciences such as biomathematics and bioscience especially in the last century. One
of most important mathematical topics using in biology, ecology etc. is the solution find-
ing for system of differential equations. Most of time finding of solution for a differential
equation is impossible unfortunately. On other hand having an initial population for a
system of differential equations, an initial value problem may be formed. We should find
out the solution or analyze its behavior. Anybody knows the importance and effect of
single species in community and nature, and so we present a summary and brief literature
of mathematical models such as single species, exponential and logistic modeling used to
describe population dynamics, which may help us future model development and high-
lights the importance of population growth rate modeling in biology, ecology, and another
area even computer and network sciences.

The study of population change started by Leonardo of Pisa, who was only given
the nickname Fibonacci. He introduced in his arithmetic book of 1202 set, a modeling
exercise involving a hypothetical growing rabbit population [9]. Environmental variation
in ecological communities and inferences from single species data is studied by Abbott
et al [2]. They present a method for estimating dimension of community on a single
species. A dynamics of single species population growth is studied by Mueller et al in
[8]. The stability at carrying capacity of environment for 25 different populations of
Drosophila melanogaster are examined by them. The spatial dynamics single species can
contribute to long term rarity and commonness. A spatial population dynamics model and
regional commonness and regional extinction are investigated in [6]. Anybody is able to
find out the answer of question ”How great an effect does self-generated spatial structure
have on logistic population growth?” which is studied by Law et al in [7]. They showed
that population growth given by logistic model may differ greatly from that of the non-
spatial logistic equation. They moreover proved that populations may achieve asymptotic
densities greater than or less than the carrying capacity of logistic model where it is not
spatial, and can even tend towards extinction. Accounting for the local spatial processes
indeed brings the theory of single species population growth a step closer to the growth
of real spatially structured populations. The population growth of infection of virus and
worm in computer networks is analyzed by help of logistic modeling [1]. Pollutant and virus



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1178

induced disease effect of on single species animal population and its essential mathematical
features by help of the logistic modeling. Moreover, It is shown that the susceptible
population does not vanish when it is only under the effect of infection meanwhile in
the polluted environment, it can extinct [5]. An application of exponential and logistic
growths for single-cell models that are incapable or misleading for inferring population
dynamics as there is no any interactions between cells via metabolites or physical contact,
nor competition for limited resources such as nutrients or space is studied in [4]. Some
more text on applications of single species with exponential or logistic growth rate having
harvesting factor can be seen in [3], [9], [10], [11] and [13].

Exponential growth is associated with the name of Thomas Robert Malthus (1766-
1834) who first realized that any species can potentially increase in numbers according to
a geometric series [12]. He published his book in 1798 stating that populations with abun-
dant natural resources grow very rapidly; however, they limit further growth by depleting
their resources. The early pattern of accelerating population size is called exponential
growth. Charles Darwin, in developing his theory of natural selection, was influenced by
Malthus.

Generally, in open population or open system, there is the following equality

Population Change=Births-Deaths+Immigrations-Emigrations

By assuming there is no any immigrations or emigrations in our population, we consider
a closed population or closed system. That is the population has changed only by the
occurrence of birth and death. Moreover, By considering x(t) denotes the size of population
at time t; b and d denote the birth rate and death rate respectively, in the time interval
[t, t + ∆t], then we get

x(t + ∆t)− x(t) = bx(t)∆t− dx(t)∆t.

⇒ dx

dt
= lim

∆t→0

x(t + ∆t)− x(t)

∆t

= lim
∆t→0

bx(t)∆t

∆t
− lim

∆t→0

dx(t)∆t

∆t

= bx(t)− dx(t),

and so, by taking r = b− d, we have

dx

dt
= rx(t). (1)

The value of parameter r as intrinsic growth rate of population can be positive, meaning
the population is increasing in size (the rate of change is positive); or negative, meaning the
population is decreasing in size; or zero, in which case the population size is unchanging
i.e. the population is constant. We now present an example of exponential growth for
bacteria population in the following:



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1179

Example 1. The famous example of exponential growth in organisms indeed is seen in
bacteria. Bacteria are prokaryotes that reproduce largely by binary fission. This division
takes about an hour for many bacterial species. If 1000 bacteria are placed in a large flask
with an abundant supply of nutrients, the number of bacteria will have doubled from 1000 to
2000 after just an hour. In another hour, each of the 2000 bacteria will divide, producing
4000 bacteria. After the third hour, there should be 8000 bacteria in the flask. The
important concept of exponential growth is that the growth rate-the number of organisms
added in each reproductive generation-is itself increasing; that is, the population size is
increasing at a greater and greater rate. After 24 of these cycles, the population would
have increased from 1000 to more than 16 billion bacteria. When the population size x(t)
is plotted over time, a J-shaped growth curve is produced.

The bacteria-in-a-flask example indeed is not truly representative of the real world
where resources are usually limited. However, when a species is introduced into a new
habitat that it finds suitable, it may show exponential growth for a while. In the case
of the bacteria in the flask, some bacteria will die during the experiment and thus not
reproduce; therefore, the growth rate is lowered from a maximal rate in which there is
no mortality. The growth rate of a population is largely determined by subtracting the
death rate (d) from the birth rate (b). The growth rate can be expressed in a simple
equation that combines the birth and death rates into a single factor r which is shown in
the previous formula.

It is clear that Malthus’s model is applicable only if the number of species is small.
Otherwise, the limited resources such as food and place limit the growth of population. The
above I.V.P. is valid for a short period and can’t go on forever. However, the exponential
growth applications in the nature is so much. Indeed, there is no an area of science
that no need to help of exponential growth. Some of them are microbiology (growth
of bacteria), conservation biology (restoration of disturbed populations), insect rearing
(prediction of yield), plant or insect quarantine (population growth of introduced species),
fishery (prediction of fish dynamics), biochemical(radioactivity).

2. Modeling and Discussion

We now in a position to construct and study some ecological I.V.Ps. of single species
models based on P. Verhulst(1838-1845) idea which is published for human population
growth rate [12]. The larger population means fewer resources such as food, place etc.
which are limited. And so, this implies a smaller rate of growth rate. In the simplest case,
he considered that the rate decreases linearly as a function of x.

2.1. The Exponential Model; Having Constant Harvesting Factor

Having considered x(0) = x0 initial population for equation (1), we obtain the following
the initial value problem named I.V.P. :{

dx
dt = rx
x(0) = x0

(2)



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1180

where its solution as follows:

x(t) = x0exp(rt). (3)

Adding constant harvesting factor h in I.V.P exponential growth rate (2), we obtain
the following I.V.P.: {

dx
dt = rx− h
x(0) = x0

(4)

Where the above I.V.P. solution is given by

x(t) = (x0 −
h

r
)exp(rt) +

h

r
. (5)

We are now going to analyze the solution of I.V.Ps. (2) and (4) by Picard’s iteration
method.

Theorem 1. The following statements for I.V.Ps. of exponential growth models are true:
(i) The series solution of I.V.P. (2) is

x(t) = x0

k=∞∑
k=0

(rt)k

k!
(6)

(ii) The series solution of I.V.P. (4) is

x(t) =
h

r
+ (x0 −

h

r
)
k=∞∑
k=0

(rt)k

k!
(7)

Proof. (i) First, we should consider the following iteration scheme:

xn(t) = x0 +

∫ s=t

s=0
rxn−1(s)ds. (8)

At the first step, we calculate x1(t)

x1(t) = x0 +

∫ s=t

s=0
rx0ds = x0 + rx0t = xo(1 + rt),

and so,

x1(t) = xo(1 + rt). (9)

At the second step, we calculate the x2(t)

x2(t) = x0 +

∫ s=t

s=0
rx1(s)ds

= x0 +

∫ s=t

s=0
rx0(1 + rs)ds



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1181

= x0 + x0[rs +
(rs)2

2
]s=t
s=0

= x0 + x0t(r +
(rt)2

2
).

Thus,

x2(t) = x0(1 + rt +
(rt)2

2
). (10)

At the third step, we calculate the x3(t)

x3(t) = x0 +

∫ s=t

s=0
rx2(s)ds

= x0 +

∫ s=t

s=0
rx0(1 + rs +

(rs)2

2
)ds

= x0 + x0[rs +
(rs)2

2
+

(rs)3

2× 3
]s=t
s=0

= x0(1 + rt +
(rt)2

2
+

(rt)3

2× 3
).

Hence,

x3(t) = x0(1 + rt +
(rt)2

2
+

(rt)3

2× 3
). (11)

At the forth step, we calculate the x4(t)

x4(t) = x0 +

∫ s=t

s=0
rx3(s)ds

= x0 +

∫ s=t

s=0
rx0(1 + rs +

(rs)2

2
+

(rs)3

2× 3
)ds

= x0 + x0[rs +
(rs)2

2
+

(rs)3

2× 3
]s=t
s=0

= x0(1 + rt +
(rt)2

2
+

(rt)3

2× 3
+

(rs)4

2× 3× 4
).

Then,

x4(t) = x0(1 + rt +
(rt)2

2
+

(rt)3

2× 3
+

(rt)4

2× 3× 4
). (12)

Paying attention to relations (9),(10),(11) and (12), one may guess the following series:

xn(t) = x0

k=n∑
k=0

(rt)k

k!
. (13)



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1182

Approaching n to infinity, we see that the series solution of I.V.P. (2) may be obtained
as follows:

x(t) = limn→∞xn(t) = x0

k=∞∑
k=0

(rt)k

k!
(14)

Therefore, the truth of relation (6) is proved.
(ii) Now, consider the following iteration scheme:

xn(t) = x0 +

∫ s=t

s=0
(rx(n− 1)(s)− h)ds (15)

At the first, step we calculate x1(t)

x1(t) = x0 +

∫ s=t

s=0
(rx0 − h)ds = x0 + (x0 −

h

r
)rt.

For simplifying take

A = x0 −
h

r
, (16)

which it implies that
x1(t) = x0 + Art. (17)

At the second step, we calculate the x2(t)

x2(t) = x0 +

∫ s=t

s=0
(rx1(s)− h)ds

= x0 +

∫ s=t

s=0
(r(x0 + Ars)− h)ds

= x0 + [rx0s + A
(rs)2

2
− hs]s=t

s=0

= x0 + rx0t + A
(rt)2

2
− ht.

Then,

x2(t) = x0 + A(rt +
(rt)2

2
). (18)

At the third step, we calculate the x3(t)

x3(t) = x0 +

∫ s=t

s=0
(rx2(s)− h)ds

= x0 +

∫ s=t

s=0
(r(x0 + Ars + A

(rs)2

2
)− h)ds



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1183

= x0 + [rx0 + Ar2s + A
r3

2
s2 − hs]s=t

s=0

= x0 + rx0t− ht + A
(rt)2

2
+ A

(rt)3

2
.

Thus,

x3(t) = x0 + A(rt +
(rt)2

2
+

(rt)3

2× 3
). (19)

At the forth step, we calculate the x4(t)

x4(t) = x0 +

∫ s=t

s=0
(rx3(s)− h)ds

= x0 +

∫ s=t

s=0
(r(x0 + Ars + A

(rs)2

2
+ A

(rs)3

2× 3
)− h)ds

= x0 + [rx0 + Ar2s + A
r3

2
s2 + A

r4

2× 3
s3 − hs]s=t

s=0

= x0 + rx0t− ht + A
(rt)2

2
+ A

(rt)3

2
+ A

(rt)4

2× 3× 4
.

Therefore,

x4(t) = x0 + A(rt +
(rt)2

2
+

(rt)2

2× 3
+

(rt)4

2× 3× 4
). (20)

Paying attention to relations (17), (18), (19) and (20), we guess the following series:

xn(t) = x0 + A
k=n∑
k=1

(rt)k

k!
. (21)

And so,

xn(t) = x0 −A
(rt)0

0!
+ A

k=n∑
k=0

(rt)k

k!

= x0 − (x0 −
h

r
) + A

k=n∑
k=0

(rt)k

k!
.

Therefore,

xn(t) =
h

r
+ A

k=n∑
k=0

(rt)k

k!
. (22)

Approaching n to infinity, we see that the series solution of I.V.P. (4) may be obtained
as follows:

x(t) = limn→∞xn(t) =
h

r
+ A

k=∞∑
k=0

(rt)k

k!



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1184

which implies that

x(t) =
h

r
+ A

k=∞∑
k=0

(rt)k

k!
(23)

Regarding the relation (16), we may work out the relation (7) which completes the
proof of (ii). Therefore, the proof of theorem is done.

2.2. The Logistic Model without Harvesting Factor

Extended exponential growth is possible only when infinite natural resources are avail-
able; this is not the case in the real world. Charles Darwin recognized this fact in his
description of the ”struggle for existence,” which states that individuals will compete
(with members of their own or other species) for limited resources. The successful ones
are more likely to survive and pass on the traits that made them successful to the next
generation at a greater rate (natural selection). To model the reality of limited resources,
population ecologists developed the logistic growth model.

In the real world, with its limited resources, exponential growth cannot continue in-
definitely. Exponential growth may occur in environments where there are few individuals
and plentiful resources, but when the number of individuals gets large enough, resources
will be depleted and the growth rate will slow down. Eventually, the growth rate will
plateau or level off. This population size, which is determined by the maximum popula-
tion size that a particular environment can sustain, is called the carrying capacity. In real
populations, a growing population often overshoots its carrying capacity, and the death
rate increases beyond the birth rate causing the population size to decline back to the
carrying capacity or below it. Most populations usually fluctuate around the carrying
capacity in an undulating fashion rather than existing right at it.

We are now in a position model the logistic growth. By this mean consider positive
parameters M and r as the carrying capacity of the environment and the rate of growth
for small population numbers, respectively. It is obvious that the factor rx represents

unhampered growth and reduced by the term
r

M
x2 that corresponds to competition for

limited resources within the population. Thus the form the following equation may be
obtained:

dx

dt
= rx(1− x

M
). (24)

This equation is known as the logistic equation. In precisely speaking, the above formula
used to calculate logistic growth adds the carrying capacity as a moderating force in the
growth rate. The expression M − x is equal to the number of individuals that may be
added to a population at a given time, and we saw that M − x is divided by M is the
fraction of the carrying capacity available for further growth. Thus, the exponential growth
model is restricted by this factor to generate the logistic growth equation. Therefore, we
can take a result as ”carrying capacity is the most and important subject in any logistic
modeling”. When the population size is equal to the carrying capacity, the quantity in



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1185

parentheses is equal to zero and growth is equal to zero. A graph of this equation yields the
S-shaped curve. It is a more realistic model of population growth than exponential growth.
There are three different sections to an S-shaped curve. Initially, growth is exponential
because there are few individuals and ample resources available. Then, as resources begin
to become limited, the growth rate decreases. Finally, the growth rate levels off at the
carrying capacity of the environment, with little change in population number over time.

Making assumption x(t0) = x0 as the initial population for time t0 in logistic equation
(24), we may obtain an ecological I.V.P. For simplifying, we denoted t0 by 0.{

dx
dt = rx(1− x

M
)

x(0) = x0.
(25)

By using integration part by part, the solution of above I.V.P. may be found as follows:

x(t) =
Mx0

x0 + (M − x0)exp(−rt)
(26)

It is clear that the following properties for the last I.V.P. may be obtained:

(a) It has two equilibria x = 0,M . That is the density does not change in these cases.
For 0 < x < M , it increases, and for x > M it decreases;

(b) If x0 < M , then the population grows and approaches to M asymptotically as
t→∞;

(c) If x0 > M , then the population decreases, again approaching M asymptotically as
t→∞;

(d) If x0 = M , then the population remains in time at x = M .

(e) Equilibrium pointx = M is globally stable; i.e.

lim
t→∞

x(t) = M ;

(f) The behavior of solution in the region between 0 and M is known as logistic growth.

The above discussion is verified by numerical simulation for some different values of r
and M in section 3 and shown in figure (1).

2.3. The Logistic Modeling with Constant Harvesting Factor

Having assumed the positive constant number h of population removed per each
duration, we can extend I.V.P. (25) as follows:{

dx
dt = rx(1− x

M
)− h

x(0) = x0.
(27)

The above I.V.P. exhibits the limited population.



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1186

Theorem 2. Ecological I.V.P. (27) has two equilibria x1 and x2 as follows:

x1 =
M

2
(1−

√
1− 4h

rM
) , x2 =

M

2
(1 +

√
1− 4h

rM
). (28)

Moreover, the above equilibrium points are unstable and asymptotic stable, respectively.

Proof. It is clear that by setting dx
dt , we get r

M x2 − rx + h = 0. The roots of last
equation are equilibria of I.V.P. which are in (28).

Paying attention to the above equilibria, we are able to see the following properties
immediately:

(a) Regarding x(t) denotes the density or number of population, we have h < rM
4 ;

(b) x1 is positive provided
√

1− 4h
rM < 1;

(c) Since h > 0, 0 < x1 < x2 < M .

Regarding relations (28), we get:

rx(1− x

M
)− h = a(x− x1)(x− x2). (29)

And so. I.V.P. (27) leads to the following I.V.P.:{
dx
dt = a(x− x1)(x− x2)
x(0) = x0.

(30)

where its solution is given by:

x(t) =
x2(x0 − x1)− x1(x0 − x2)exp(− r

M (x2 − x1)t)

(x0 − x1)− (x0 − x2)exp(− r
M (x2 − x1)t)

(31)

Since the following inequality is true

− r
M (x2 − x1)t < 0 , for all t > 0;

limn→∞exp(− r

M
(x2 − x1)t) = 0.

If x0 > x1, we get x(t) = x2.
Thus, x(t) = x2 is an equilibrium limiting solution. Therefore, the point of x = x2 is

stable point.
We now analyze the solution behavior of equation (31) in case of x0 < x1 .
By assuming

x0 − x1 = (x0 − x2)exp(− r

M
(x2 − x1)t),

we get



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1187

t =
1

− r
M (x2 − x1)

ln
x0 − x1

x0 − x2
.

Now, let t1 be as following constant value:

t1 =
1

− r
M (x2 − x1)

ln
x0 − x1

x0 − x2
. (32)

By setting the value of t1 in solution (31) we have:

x(t1) =
x2(x0 − x1)− x1(x0 − x2)exp[− r

M (x2 − x1) 1
− r

M
(x2−x1) ln

x0−x1
x0−x2

]

(x0 − x1)− (x0 − x2)exp[− r
M (x2 − x1) 1

− r
M

(x2−x1) ln
x0−x1
x0−x2

]
.

Since the numerator of the above fraction is negative and it’s denominator is zero, we get
x(t1) = −∞.

Therefore, the value of the population function x(t) at t = x1 is a threshold and
consequently, point x = x1 is a unstable equilibrium point. This means that the proof is
completed.

2.4. The Logistic Modeling with Variable Harvesting Factor

Now we are going to extend the logistic population harvesting into case of har-
vesting parameter h in I.V.P. (27) is not constant. Consider the situation that the pop-
ulation has logistic growth rate and harvesting coefficient is not constant. Therefore, the
general case of this model is as follows:

dx

dt
= rx(1− x

M
)− f(x) (33)

where f(x) is an arbitrary function. There is no general way to solve the above
equation; and so we should restrict our attention to a few special cases. In continuation,
we study two cases for variable harvesting factor which are cubic and fractional.

2.4.1. Case of Fractional Harvesting Factor

Let us make assumption that the harvesting factor be following function:

f(x) = h
x

1 + x
. (34)

And so, one of the extension of the logistic population harvesting factor will be ap-
peared as follows:

dx

dt
= rx(1− x

M
)− h

x

1 + x
(35)



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1188

The last equation exhibits that the harvesting coefficient for each term depends on it’s
population density. By adding initial population x(0) = x0 in the equation (35), we get
the following I.V.P.: {

dx
dt = rx(1− x

M
)− h x

x+1

x(0) = x0.
(36)

Since x
1+x < 1, the term of harvesting factor has inverse relation respect to density of

population. In the other word, the term of harvesting decreases as population increases.

Theorem 3. The following statements for logistic modeling I.V.P. (36) are true:
(i) It has three equilibria x1, x2 and x3 as follows:

x1 = 0, x2,3 =
M − 1

2
± M

2

√
(

1

M
− 1)2 +

4

M
(1− h

r
); (37)

(ii) The equilibria x2 and x3 are real number provided

h

r
≤ 1 +

M

4
(1− 1

M
)2;

(iii) It is impossible that both equilibria x2 and x3 are positive;
(iv) The solution of this I.V.P. is as follows:

xB1(x− x2)B2(x− x3)B3 = x0
B1(x0 − x2)B2(x0 − x3)B3exp(

−rt
M

) (38)

(v) In case of M = r = h = 1, I.V.P. has equilibria x = 0 which is stable point.

Proof.
(i)
We first set dx

dt = 0. by a simple calculation, we find equilibria which are given (37)
easily.

(ii)
paying attention to (37), we see that for equilibria x2 and x3 are real number provided

h

r
≤ 1 +

M

4
(1− 1

M
)2;

(iii)
Now, let us take they are positive. It implies that

1− 1

M
<

√
(1− 1

M
)2 − 4

M
(
h

r
− 1) <

1

M
− 1.

If we take M < 1, the following contradiction is followed:√
(1− 1

M
)2 − 4

M
(
h

r
− 1) < 0



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1189

which implies that the truth of (iii).
(iv)
If we take M = 1, we see that

x2,3 = ±
√

1− h

r

For simplifying we consider

dx

dt
=
−r
M

x

x + 1
(x− x2)(x− x3) (39)

Then, we get
1 + x

x(x− x2)(x− x3)
dx =

−r
M

dt.

By decomposing method, we may find the parameters B1, B2 and B3 which satisfy in
the following equation:

1 + x

x(x− x2)(x− x3)
=

B1

x
+

B2

x− x2
+

B3

x− x3
.

By multiply the above equation into the phrases x, x − x2 and x − x3; and setting
x = 0, x = x2 and x = x3 respectively in three steps, we have

B1 =
1

x2x3
, B2 =

1 + x2

x2(x2 − x3)
, B3 =

1 + x2

x3(x3 − x2)
.

Thus, ∫
1 + x

x(x− x2)(x− x3)
dx = Ln[xB1(x− x2)B2(x− x3)B3 ]

which it implies that

xB1(x− x2)B2(x− x3)B3 = Cexp(
−rt
M

), (40)

where C is constant.
By setting x(0) = x0 in (40), we see that the solution given by formula (38) is right.
(v) Taking M = r = h = 1, we see equation (35) leads to the following equation:

dx

dt
= x(1− x)− x

1 + x

⇒ x3 = 0

Thus x = 0 is equilibria of (35). As regarding x− x2 − x
x+1 = − x3

x+1 < 0; we have

x− x2 − x

x + 1
< 0



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1190

And so,x(t) is decreasing.
Therefore,

lim
t→∞

x(t) = 0.

The above population function is an implicit function and so, we can’t sketch the graph
of solution x(t) respect to t clearly.

Applying variation table for function g(x) = rx(1− x

M
)−h x

x+1 , one is able to determine

the monotonicity (increasing or decreasing) of the above implicit function.

In the next section, we describe the solution behavior by help of method of numerical
simulation.

2.4.2. Case of Cubic Harvesting Factor

In the second case, we make assumption that the harvesting function be as follows:

f(x) = x3, (41)

and so we get:
dx

dt
= rx(1− x

M
)− hx3 (42)

Like the previous case, by adding initial population x(0) = x0 in the equation (42), we
get the following I.V.P.: {

dx
dt = rx(1− x

M
)− hx3

x(0) = x0.
(43)

Theorem 4. The following statements for logistic modeling I.V.P. (43) are true:
(i) It has three equilibria x1, x2 and x3 as follows:

x1 = 0, x2,3 =
r ±∆

2hM
where ∆ =

√
r2 + 4hrM2; (44)

(ii) The solution of this I.V.P. is the following implicit function:

x
h
r (x− x2)

r−∆
−2∆ (x− x3)

r+∆
−2∆ = x

h
r
0 (x0 − x2)

r−∆
−2∆ (x0 − x3)

r+∆
−2∆ ; (45)

(iii) In case of M = r = h = 1, the solution of this I.V.P. is decreasing function.
Moreover, the equilibria x1 = 0 is asymptotically stable.

Proof.
(i)
By setting dx

dt = 0, then we see that the first equilibrium point is x1 = 0.
Another equilibria are the roots of following quadratic polynomial:

−hx2 − r

M
x + r,



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1191

⇒ x2,3 =
r ±
√
r2 + 4hrM

2hM
,

which implies the truth of (44).
(ii)
Then, we get

dx

dt
= x(x− r −∆

2hM
)(x− r + ∆

2hM
) (46)

.
After calculation, we get:∫

dt =

∫
(
h
r

x
+

r−∆
−2∆

x− x2
+

r+∆
−2∆

x− x3
)dx,

We therefore obtain the solution of (42) as follows:

x
h
r (x− x2)

r−∆
−2∆ (x− x3)

r+∆
−2∆ = Kexp(t). (47)

Therefore, considering initial population x(0) = x0, we see that the solution of I.V.P.
(43) may be found as follows:

x
h
r (x− x2)

r−∆
−2∆ (x− x3)

r+∆
−2∆ = x

h
r
0 (x0 − x2)

r−∆
−2∆ (x0 − x3)

r+∆
−2∆ exp(t) (48)

which is implicit function.
(iii)
Making assumption M = r = h = 1, implies that the equation (42) leads to the

following equation:

dx

dt
= x(1− x)− x3

As regarding x(t) describes the number of population which is positive integer. And
so it takes the number one at least. we see that

x < x2 < x3

⇒ dx

dt
= x− x2 − x3 < 0

Therefore, x(t) is decreasing function.

lim
t→∞

x(t) = 0.

Therefore, the equilibria x1 = 0 is asymptotically stable. Because of being implicit the
solution x(t), it is impossible to sketch the graph of this solution. Applying variation table,
one is able to determine the monotonicity (increasing or decreasing) of the above implicit
function. In the next section, we describe the solution behavior by help of simulation
method.

In the next section, we describe the solution behavior by help of method of numerical
simulation.



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1192

3. Numerical Simulation

In this section, we discuss the analytical results by making some simulations on stud-
ied I.V.P. describing logistic modeling having and without harvesting factor analyzed in
section 2. First, the behavior of solution (26) for I.V.P. (25) is simulated. An numerical
simulation for various values to the carrying capacity M and growth rate r is formulated
as follows:.

(i) x0 = M + i, where i = 20, 40, 60 and 80;

(ii) x0 = M
2 + j, where i = 10, 20, 30 and 40;

(iii) x0 = M
k , where i = 3, 4, 5 and 6.

Indeed, in each column of table (1) by considering carrying capacity (M), we worked out
formula for finding suitable initial populations. And so, for rows a, b, c, d, e, f; we have:

In case (i): Calculated initial populations is greater than the carrying capacity. In this
case, obtained graphs shows the being of solutions decreasing.

In case (ii): Calculated initial populations are located between M
2 and M . That is

they are greater than M
2 and less than the carrying capacity M .

In case (iii): Calculated initial populations is less than M
2 . In this case, obtained

graphs shows that all solutions have turning behavior at the point x = M
2 .

In all of above cases, the solutions tend to M asymptotically. This numerical argument
is shown in the table (1).

The graphs related to solution (26) are drawn in figure 1 (a, b, c and d) which verify
the presented mathematical discussion. Indeed, green graphs, which are started with value
greater than M , are decreasing and tend to M asymptotically. Blue graphs, which are lo-
cated betweenM

2 and M , are increasing and tend to M asymptotically. Red graphs, which
are started with value less than the M

2 , are increasing and tend to M asymptotically. It is
clear that it decrease slowly whenever growth rate r is less than 1, meanwhile it decreases
rapidly whenever growth rate r is greater than 1.

To verify the the results of theorem 2.2, we make simulation for solution (31) of
I.V.P.(27) describing logistic modeling with constant harvesting factor. By giving some
different values to carrying capacity M , growth rate r and harvesting factor h, we obtain
the related equilibria x1 and x2. This argument has been brought in Table 2.

Paying attention to drawn graphs in figure 2, we see that the equilibrium point x1 is
unstable and the equilibrium point x2 is asymptotically stable which are proved in theo-
rem 2.2. The details may be seen in figure 2. The solution graphs are shown in this figure
(a,b,c,d) clearly.

We now make simulation for solution (38) of I.V.P. (36). This I.V.P. describes the
logistic modeling variable harvesting factor. Indeed, in this case harvesting factor is a
fractional function f(x) = x

x+1 . The details of parameters: carrying capacity(M), growth



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1193

rate(r), harvesting factor(h) and initial populations(x0) are presented in table (2). We see
that for case of M = r = h = 1 all of equilibria for (36) are x = 0. The related graphs
verifying the results of theorem 2.3 are shown in figure (3). These graphs are simulated
in two scale 0 ≤ t < 30 and 0 ≤ t < 1000.

For the final case, we make simulation for solution of (45) of I.V.P. (43). This I.V.P.
describes logistic modeling having variable harvesting factor. It is assumed that in this
model harvesting factor is cubic function f(x) = x3. Considering assumption as M = r =
h = 1 implies that equilibria are given by

x1 = 0, x2 = −1.618and x3 = 6.18× 10.

If we consider M = ×106, r = 1, h = ×10−6, we have equilibria as follows:

x1 = 0, x2 = −1.0005× 103 and x3 = 9.995× 103.

In both of the above cases, the equilibria are negative, zero and positive which are
proved in theorem 2.4. The another result of solution behavior which are studied in said
theorem are shown in figure (4).

Table 1: Parameters, Equilibria and Initial Populations for I.V.P.(25)

a b c d

M 102 102 103 103

r 5 × 10 2 0.5 2
x1 0 0 0 0

x2 102 102 103 103

x01 1.667 × 10 1.666 × 102 50 1.666 × 102

x02 20 20 200 2 × 102

x03 25 25 250 2.5 × 102

x04 3.333 × 10 3.333 × 102 3.333 × 103 3.333 × 103

x05 6 × 10 6 × 10 5.1 × 102 5.1 × 102

x06 7 × 10 7 × 10 5.2 × 102 5.2 × 102

x07 8 × 10 8 × 10 5.3 × 102 5.3 × 102

x08 9 × 10 9 × 10 5.4 × 102 5.4 × 102

x09 1.2 × 10 1.2 × 102 1.2 × 103 1.2 × 103

x010 1.4 × 10 1.4 × 10 1.4 × 103 1.4 × 103

x011 1.6 × 10 1.6 × 10 1.6 × 103 1.6 × 103

x012 1.8 × 10 1.8 × 10 1.8 × 103 1.8 × 103

Table 2: Parameters, Equilibria and Initial Populations for I.V.P.(27)

a b c d

M 104 104 104 104

r .5 .5 10 10

h 102 102 103 103

x1 0 0 0 0

x2 2.041 × 102 2.041 × 102 1.001 × 10 1.001 × 10

x3 9.795 × 103 9.795 × 103 9.899 × 103 9.899 × 103

x01 2.2 × 102 1.8 × 102 50 10

x02 5 × 102 5 × 102 5 × 102 5 × 102

x03 2 × 103 2 × 103 2 × 103 2 × 103

x04 9 × 103 9 × 103 9 × 103 9 × 103

x05 1.2 × 103 1.2 × 103 1.3 × 103 1.3 × 103

x06 1.5 × 103 1.5 × 103 1.5 × 103 1.5 × 103



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1194

Table 3: Parameters, Equilibria and Initial Populations for I.V.P.(36)

a b

M 1 106

r 1 1

h 1 10−6

x1 0 0

x2 0 −9.9999 × 10−1

x3 0 9.9999 × 105

x01 105 106

x02 9 × 104 9 × 105

x03 8 × 104 8 × 105

x04 7 × 104 7 × 105

x05 6 × 104 6 × 105

x06 5 × 104 5 × 105

x07 4 × 104 4 × 105

x08 3 × 104 3 × 105

x09 2 × 104 2 × 105

x10 104 105

x011 6 × 103 6 × 104

x012 4 × 103 4 × 104

x013 2 × 103 2 × 104

x014 103 104

Table 4: Parameters, Equilibria and Initial Populations for I.V.P.(43)

a b

M 1 106

r 1 1

h 1 10−6

x1 0 0

x2 −1.618 −1.0005 × 103

x3 6.18 × 10−1 9.9950 × 103

x01 105 106

x02 9 × 104 9 × 105

x03 8 × 104 8 × 105

x04 7 × 104 7 × 105

x05 6 × 104 6 × 105

x06 5 × 104 5 × 105

x07 4 × 104 4 × 105

x08 3 × 104 3 × 105

x09 2 × 104 2 × 105

x10 104 105

x011 6 × 103 6 × 104

x012 4 × 103 4 × 104

x013 2 × 103 2 × 104

x014 104 104



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1195

Figure 1: I.V.P. (2.24)

Figure 2: I.V.P. (2.26)



M.H. R. Doust , V. Lokesha, A. Ghasemabadi / Eur. J. Pure Appl. Math, 13 (5) (2020), 1176-1198 1196

Figure 3: I.V.P. (2.35)

Figure 4: I.V.P. (2.42)



REFERENCES 1197

4. Conclusion

The importance and effects of single species in the nature and community is clear
for anybody. The discussed models in this work have so many practical application. In
deed, by help of these models, one may predict, check and defence to spread of viruses,
microbe and bacteria in a community. There are so many dangerous viruses for hu-
manity. Anybody gets involved with microbe viruses and such as Black Death, Spanish
Flu, HIV/AIDS, Swine Flu, Ebola virus, Zika virus, Corona viruses such as: SARS-Cov,
MERS-Cov, COVID-19. Especially, last one which is the most deadly and disastrous
viruses nowadays in the world. This viruses is spread in 2019, and all of the countries in
the world get involved by it. Some of them such as Corona viruses: COVID-19 are not
only epidemic but also are pandemic. As a consequence in this research, we worked out
the series solution for exponential modeling for both cases of having constant harvesting
factor or simple model without harvesting factor. Moreover, making some conditions for
single species of logistic modeling, we find out they have stable solutions. And also,their
asymptotical stability is obtained. Making some various simulations, we observed the
obtained results which are proved theorems are true. The important result is: Carrying
capacity of the environment, growth rate of population, harvesting factor are so important
to stability the equilibria.

References

[1] R. Aarthee, D. Ezhilmaran, A Logistic Model For The Population Of Virus Growth
In Local Area Network, International Journal of Pure and Applied Mathematics,
115(9):401-409, 2017.

[2] K.C. Abbott, J. Ripa, A.R. Ives, Environmental variation in ecological communities
and inferences from single species data, Ecology, (2009).

[3] J. Barlow , Nonlinear and Logistic Growth In Experimental Populations of Guppies,
Ecology, (1992).

[4] D.A. Charlebois, G.Balazsi, Modeling Cell Population Dynamics, In Silico Biology,
13:21-39, 2019.

[5] S. Chauhan, O.P. Misra, Modeling and Analysis of a Single Species Population with
Viral Infection in Polluted Environment , Applied Mathematics, 3: 662-672(2012).

[6] I. Hanski, Single Species Spatial Dynamics May Contribute to Long?Term Rarity and
Commonness, Ecology, (1985).

[7] R. Law, D.J. Murrell, U. Dieckmann, Population Growth In Space And Time: Spatial
Logistic Equations, Ecology, (2003).

[8] L.D. Mueller, F.J. Ayala, Dynamics of Single Species Population Growth: Stability
or Chaos?, Ecology, (1981).



REFERENCES 1198

[9] J.D. Murrray, Mathematical Biology I. An Introduction, Springer, New York, (2002).

[10] C.Pao, C.Hao Lin, A New Approach to The Logistic Modeling Population; Having
Harvesting Factor, Yugoslav Journal of Operations Research, 26(3):381-392,2016.

[11] M.H. Rahmani Doust, M. Saraj, The Logistic Modeling Population; Having Harvest-
ing Factor, Yugoslav Journal of Operations Research 25(1):107-115, 2013.

[12] S. Ruan, Delay Differential Equation In Single Species Dynamics, University of Mi-
ami, (2006).

[13] E. Sorouri, M. Eshaghi Gourji, R.Memarbashi, An Ajnalysis of a Model with Nonlin-
ear Harvesting function, Int. J. Nonlinear Anal Appl., 11(1):,37-80, 2020.