# Some Comparison of Solutions by Different Numerical Techniques on Mathematical Biology Problem

Authors: Susmita Paul; Sankar Prasad Mondal; Paritosh Bhattacharya; Kripasindhu Chaudhuri
Journal: International Journal of Differential Equations
 (2016)

---

## Abstract

We try to compare the solutions by some numerical techniques when we apply the methods on some mathematical biology problems. The Runge-Kutta-Fehlberg (RKF) method is a promising method to give an approximate solution of nonlinear ordinary differential equation systems, such as a model for insect population, one-species Lotka-Volterra model. The technique is described and illustrated by numerical examples. We modify the population models by taking the Holling type III functional response and intraspecific competition term and hence we solve it by this numerical technique and show that RKF method gives good results. We try to compare this method with the Laplace Adomian Decomposition Method (LADM) and with the exact solutions.


---

## Body

## 1. Introduction

### 1.1. Population Growth Model

Mathematical models of population growth have been formed to provide a significant angle of the real ecological situation. The meaning and importance of each parameter in the models have been defined biologically [1, 2]. In case of insect population, birth and death rate of a species typically are not constant; instead, they vary periodically with the passage of seasons, whereas the Lotka-Volterra equations demonstrate an arbitrary number of ecological competitors (or predator-prey) model which is dynamic in nature. The population growth model is very important in mathematical biology which is used basically to demonstrate a simple nonlinear control system in population growth.


### 1.2. Holling Type III Functional Response

For getting exact population Holling type III functional response plays an important role in population dynamics. Holling type III functional response should be taken into the predator-prey interactions, which is proposed by Holling [3, 4] based on the fundament of experiment. The Holling type III functional response in ecological model is the intake rate of a predator species as a function of prey density.


### 1.3. Numerical Techniques

In the field of science and technology, numerous significant physical phenomena are frequently modeled by nonlinear differential equations. Such equations are often very difficult to solve analytically. Yet, analytical approximate methods are very important for obtaining the accurate solutions which have gained much significance in recent years. There are various methods, undertaken to find out approximate solutions to nonlinear problems. Homotopy Perturbation Method (HPM), Homotopy Analysis Method (HAM), Differential Transform Method (DTM), Variational Iteration Method (VIM), Adomian Decomposition Method (ADM), Laplace Adomian Decomposition Method (LADM), and Runge-Kutta-Fehlberg (RKF) method are some very popular methods. The purpose of this paper is to bring out the numerical solution of various population models by using the approach, namely, Runge-Kutta-Fehlberg (RKF) method.Recently, different scientists use the numerical method in their different problems. Here we try to give some references showing the importance of quasi-numerical methods techniques in present time: Kumar and Baskar [5] consider B-spline quasi-interpolation based numerical methods for some Sobolev type equations. Analytical and numerical study of dirty bosons in a quasi-one-dimensional harmonic trap is discussed by Khellil et al. [6]. Numerical study of a quasi-hydrodynamic system of equations for flow computation at small Mach numbers is considered by Balashov and Savenkov [7]. Quasi-optimal complexity of adaptive finite element method for linear elasticity problems in two dimensions is discussed by Liu et al. [8]. Numerical solutions of quasi-two-dimensional models for laminar water hammer problems are considered by Zhao [9].


### 1.4. Motivation

The accurate solutions of population growth models may become a difficult task if the equations are highly nonlinear. To overcome the situation there we take the numerical simulation, until there are no particular numerical methods for solving such problems. So to fill up the gap, here we find the approximated solutions of some population models by applying such some reliable, efficient, and more comfortable numerical technique (e.g., LADM, RKF) and try to conclude which one is the best.


### 1.5. Novelties

The principal aim of this paper is to perform systematic analysis of the comparisons among exact solution and some reliable numerical techniques on the dynamics of the population model. The said numerical methods shall be more acceptable and reliable for solving such kind of problem. The issues which are addressed in the paper are described in the following:(i)
Adding Holling type III functional response in insect population model and Lotka-Volterra model for better formulation
(ii)
Applying two numerical methods LADM and RKF for solution of the models and find the solution
(iii)
Analysis of the comparisons among exact solution, Laplace Adomian Decomposition Method (LADM), and Runge-Kutta-Fehlberg (RKF) method on the models
(iv)
The solution of the models by different numerical techniques being illustrated numerically and graphically
(v)
The necessary algorithm for numerical techniques givenMoreover, we can say all these developments can help the researchers who engage with nonlinear differential equation, mathematical biology, and numerical techniques.


### 1.6. Structure of the Paper

The paper is organized as follows. The basic literature survey on population growth model, functional response, and numerical techniques is addressed in Section1. In Section 2 we discuss the numerical methods, namely, Runge-Kutta-Felhberg (RKF) method and Laplace Adomian Decomposition Method (LADM) for nonlinear equation. The necessary algorithm for finding the numerical results is also discussed in this section. Section 3 is followed by a numerical example in different mathematical biology models. The different biological model is formulated and we write the general expression of the numerical solution. Section 4 is illustrated by the comparison study between the solutions and error terms of these said methods. Finally conclusions and future research scope of this paper are drawn in the last section, Section 5.


## 1.1. Population Growth Model

Mathematical models of population growth have been formed to provide a significant angle of the real ecological situation. The meaning and importance of each parameter in the models have been defined biologically [1, 2]. In case of insect population, birth and death rate of a species typically are not constant; instead, they vary periodically with the passage of seasons, whereas the Lotka-Volterra equations demonstrate an arbitrary number of ecological competitors (or predator-prey) model which is dynamic in nature. The population growth model is very important in mathematical biology which is used basically to demonstrate a simple nonlinear control system in population growth.


## 1.2. Holling Type III Functional Response

For getting exact population Holling type III functional response plays an important role in population dynamics. Holling type III functional response should be taken into the predator-prey interactions, which is proposed by Holling [3, 4] based on the fundament of experiment. The Holling type III functional response in ecological model is the intake rate of a predator species as a function of prey density.


## 1.3. Numerical Techniques

In the field of science and technology, numerous significant physical phenomena are frequently modeled by nonlinear differential equations. Such equations are often very difficult to solve analytically. Yet, analytical approximate methods are very important for obtaining the accurate solutions which have gained much significance in recent years. There are various methods, undertaken to find out approximate solutions to nonlinear problems. Homotopy Perturbation Method (HPM), Homotopy Analysis Method (HAM), Differential Transform Method (DTM), Variational Iteration Method (VIM), Adomian Decomposition Method (ADM), Laplace Adomian Decomposition Method (LADM), and Runge-Kutta-Fehlberg (RKF) method are some very popular methods. The purpose of this paper is to bring out the numerical solution of various population models by using the approach, namely, Runge-Kutta-Fehlberg (RKF) method.Recently, different scientists use the numerical method in their different problems. Here we try to give some references showing the importance of quasi-numerical methods techniques in present time: Kumar and Baskar [5] consider B-spline quasi-interpolation based numerical methods for some Sobolev type equations. Analytical and numerical study of dirty bosons in a quasi-one-dimensional harmonic trap is discussed by Khellil et al. [6]. Numerical study of a quasi-hydrodynamic system of equations for flow computation at small Mach numbers is considered by Balashov and Savenkov [7]. Quasi-optimal complexity of adaptive finite element method for linear elasticity problems in two dimensions is discussed by Liu et al. [8]. Numerical solutions of quasi-two-dimensional models for laminar water hammer problems are considered by Zhao [9].


## 1.4. Motivation

The accurate solutions of population growth models may become a difficult task if the equations are highly nonlinear. To overcome the situation there we take the numerical simulation, until there are no particular numerical methods for solving such problems. So to fill up the gap, here we find the approximated solutions of some population models by applying such some reliable, efficient, and more comfortable numerical technique (e.g., LADM, RKF) and try to conclude which one is the best.


## 1.5. Novelties

The principal aim of this paper is to perform systematic analysis of the comparisons among exact solution and some reliable numerical techniques on the dynamics of the population model. The said numerical methods shall be more acceptable and reliable for solving such kind of problem. The issues which are addressed in the paper are described in the following:(i)
Adding Holling type III functional response in insect population model and Lotka-Volterra model for better formulation
(ii)
Applying two numerical methods LADM and RKF for solution of the models and find the solution
(iii)
Analysis of the comparisons among exact solution, Laplace Adomian Decomposition Method (LADM), and Runge-Kutta-Fehlberg (RKF) method on the models
(iv)
The solution of the models by different numerical techniques being illustrated numerically and graphically
(v)
The necessary algorithm for numerical techniques givenMoreover, we can say all these developments can help the researchers who engage with nonlinear differential equation, mathematical biology, and numerical techniques.


## 1.6. Structure of the Paper

The paper is organized as follows. The basic literature survey on population growth model, functional response, and numerical techniques is addressed in Section1. In Section 2 we discuss the numerical methods, namely, Runge-Kutta-Felhberg (RKF) method and Laplace Adomian Decomposition Method (LADM) for nonlinear equation. The necessary algorithm for finding the numerical results is also discussed in this section. Section 3 is followed by a numerical example in different mathematical biology models. The different biological model is formulated and we write the general expression of the numerical solution. Section 4 is illustrated by the comparison study between the solutions and error terms of these said methods. Finally conclusions and future research scope of this paper are drawn in the last section, Section 5.


## 2. Numerical Solution of Nonlinear Differential Equation

### 2.1. Laplace Adomian Decomposition Method (LADM)

#### 2.1.1. LADM for Nonlinear Differential Equation

The Laplace Adomian Decomposition Method (LADM) [10, 11] was firstly introduced by Khuri and has been successfully used to find the solution of linear and nonlinear differential equations. The significant advantage of this method is that it is a combination of two powerful techniques, namely, Laplace transform and Adomian Decomposition Method [12–15]. The Laplace transform is an elementary but useful technique for solving linear ordinary differential equations that is widely used by scientists and engineers for tackling the linear models. The main impotence of this method is that the solution of this method is expressed in ∞-series which converges to the exact solution and will not take too much time for computing.Consider the following nonlinear differential equation:(1)
L
u
t
+
R
u
t
+
N
u
t
=
g
t
,where L is a linear operator of the highest-order derivative which is assumed to be invertible easily, R is the remaining linear operator of order less than L and N is a nonlinear operator, and g
(
t
) is a source term.Taking Laplace transform on both sides of the above equation, we get(2)
L
L
u
t
+
L
R
u
t
+
L
N
u
t
=
L
g
t
.Using the differential property of Laplace transform and using the initial condition, we get(3)
s
n
L
u
t
-
s
n
-
1
u
0
-
s
n
-
2
u
′
0
-
⋯
-
u
n
-
1
0
+
L
R
u
t
+
L
N
u
t
=
L
g
t
,or(4)
L
u
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
-
1
s
n
L
R
u
t
-
1
s
n
L
N
u
t
+
1
s
n
L
g
t
.Now we represent the unknown functions u
(
t
) by an infinite series of the form(5)
u
t
=
∑
n
=
0
∞
u
n
t
.Here the components u
n
(
t
) are usually determined recurrently and the nonlinear operator N
(
u
) can be decomposed into an infinite series of polynomials, given by(6)
N
u
=
∑
n
=
0
∞
A
n
,where A
n are Adomian polynomials of u
0
,
u
1
,
…
,
u
n defined by(7)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
i
=
0
∞
λ
i
u
i
λ
=
0
,
n
=
0,1
,
2
,
…
.Therefore,(8)
L
∑
n
=
0
∞
u
n
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
-
1
s
n
L
R
∑
n
=
1
∞
u
n
t
-
1
s
n
L
∑
n
=
1
∞
A
n
+
1
s
n
L
g
t
.In general, the recursive relation is given by(9)
L
u
0
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
+
1
s
n
L
g
t
,
L
u
n
+
1
t
=
-
1
s
n
L
R
u
n
t
-
1
s
n
L
A
n
.Applying the inverse Laplace transform to both sides of (9), we obtain  u
n, (
n
≥
0
), which is then substituted into (5).For numerical computation, we get the expression as(10)
ϕ
n
t
=
∑
k
=
0
n
u
k
t
,which is nth term approximation of u
(
t
) and the obtained series solution which converges to the exact solution.


#### 2.1.2. Algorithm for Finding the Numerical Solution by LADM

Step 1.
Split the given equation into two parts. The first part is(11)
-
R
u
t
+
N
u
tand the second part is(12)
L
u
t
=
g
t
.Step 2.
Apply the Laplace transform to second part, determine the coefficient ofL
u
t, and solve this equation for L
u
t. So, you can get L
u
0.Step 3.
Calculate the Adomian polynomials for the functionN
u
t. Apply the Laplace transform to these polynomials.Step 4.
Divide the first part to the coefficient ofL
u
t. In a loop, calculate L
u
n
+
1
t with substituting L
A
n and L
u
n
t into the first part.Step 5.
Construct the solution using inverse Laplace transform toL
u
n
t.Step 6.
End


### 2.2. Runge-Kutta-Fehlberg (RKF) Method

#### 2.2.1. RKF for Nonlinear Differential Equation

One of the most popular methods with a constant step size is the fourth-order Runge-Kutta (RK4) method. Reasonably the Runge-Kutta method can [16] obtain the accuracy of a Taylor Series approximation without the need of higher derivative calculations. This method can be considered as the basic form of other methods. However, in terms of error estimation, the one-step method with an adaptive step size like the Runge-Kutta-Fehlberg (RKF) method [17, 18] gives better error estimation than one-step method with a constant step size like the Runge-Kutta method. At each step the Runge-Kutta-Fehlberg method described the calculation of two Runge-Kutta methods of different order (RK4 and RK5). If the two answers are close enough then continue for the next step with the same step size. For a fixed accuracy the step size should be reduced. If the answers agree with more significant digits than required, the step size is increased. So we conclude that the one-step algorithm method associated with an adaptive step size automatically organizes the step size as a recomposition to the calculation truncation errors. This method has shown that it works in the case of nonlinear models and, hence, its application is found in wide range of deterministic and stochastic problems, linear and nonlinear problems, physics, biology, and chemical reactions problems, and so forth.Consider the initial value problem(13)
y
′
t
=
f
t
,
y
t
,
y
t
0
=
y
0
.The RKF is one way to try to resolve this problem.The problem is to solve the initial value problem in the above equation by means of Runge-Kutta methods [19, 20] of order 4 and order 5.First we need some definitions:(14)
k
1
=
h
f
t
i
,
y
i
,
k
2
=
h
f
t
i
+
1
4
h
,
y
i
+
1
4
k
1
,
k
3
=
h
f
t
i
+
3
8
h
,
y
i
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
f
t
i
+
12
13
h
,
y
i
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
f
t
i
+
h
,
y
i
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
f
t
i
+
h
,
y
i
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
.Then an approximation to the solution of initial value problem is made using Runge-Kutta method of order 4:(15)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(16)
z
i
+
1
=
y
i
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.The numerical approximation of the global (local) discretization error at the point t
i
+
1 is(17)
ϵ
=
y
i
+
1
-
z
i
+
1
.If ϵ
m
i
n
≤
ϵ
≤
ϵ
m
a
x, (where the maximum and minimum tolerance for local truncation error are ϵ
m
a
x and ϵ
m
i
n, respectively, which are prescribed at the beginning) then y
i
+
1 is an acceptable approximation of y
(
t
i
+
1
) and the step size is considerable for the next step. Otherwise, the new step size s
h can be determined by multiplying the scalar s times the step size h. The scalar s is given by(18)
s
=
ϵ
h
2
y
i
+
1
-
z
i
+
1
1
/
4
=
0.0840896
ϵ
h
y
i
+
1
-
z
i
+
1
1
/
4
,where ϵ is the specified error control tolerance.Note. RK4 requires four functions for finding the result and RK5 requires six; that is total of ten for RK4 and RK5. Fehlberg devised a method to get RK4 and RK5 results using only six function evaluations by using some of k values in both methods where k
=
∂
f
/
∂
y.


#### 2.2.2. Algorithm for Finding the Numerical Solution by RKF

Step 1.
F
(
t
,
y
1
)
← “Function to be supplied”.Step 2.
Readt
(
0
), y
1
(
0
), h, limit, ϵ
m
i
n
,
ϵ
m
a
x.Step 3.
Fori
=
0
(
1
) limit(19)
K
1
⟵
h
f
t
i
,
y
1
t
i
K
2
⟵
h
·
f
t
i
+
1
4
h
,
y
1
t
i
+
1
4
K
1
K
3
⟵
h
·
f
t
i
+
3
8
h
,
y
1
t
i
+
3
32
K
1
+
9
32
K
2
K
4
⟵
h
·
f
t
i
+
12
13
h
,
y
1
t
i
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
K
5
⟵
h
f
t
i
+
h
,
y
1
t
i
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
K
6
⟵
h
f
t
i
+
h
,
y
1
t
i
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
.Step 4.
(20)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.Step 5.
(21)
z
i
+
1
t
i
+
1
=
y
i
t
i
+
16
135
K
1
+
6656
12
,
825
K
3
+
28
,
561
56
,
430
K
4
-
9
50
K
5
+
2
55
K
6
.Step 6.
(22)
ϵ
=
y
i
+
1
-
z
i
+
1
.Step 7.
Considerϵ
m
i
n
≤
ϵ
≤
ϵ
m
a
x.Step 8.
t
i
+
1
=
t
i
+
h. Write z
i
t
i
+
1, t
i.Step 9.
Repeat for better approximation.Step 10.
End the programme.Step 11.
Otherwise,h
=
s
h with s
=
(
ϵ
h
/
2
|
y
i
+
1
-
z
i
+
1
|
)
1
/
4
=
0.0840896
(
ϵ
h
/
|
y
i
+
1
-
z
i
+
1
|
)
1
/
4.Step 12.
Repeat Steps8
–
10 for better accuracy; any one can stop here. For better accuracy the process can be run.


## 2.1. Laplace Adomian Decomposition Method (LADM)

### 2.1.1. LADM for Nonlinear Differential Equation

The Laplace Adomian Decomposition Method (LADM) [10, 11] was firstly introduced by Khuri and has been successfully used to find the solution of linear and nonlinear differential equations. The significant advantage of this method is that it is a combination of two powerful techniques, namely, Laplace transform and Adomian Decomposition Method [12–15]. The Laplace transform is an elementary but useful technique for solving linear ordinary differential equations that is widely used by scientists and engineers for tackling the linear models. The main impotence of this method is that the solution of this method is expressed in ∞-series which converges to the exact solution and will not take too much time for computing.Consider the following nonlinear differential equation:(1)
L
u
t
+
R
u
t
+
N
u
t
=
g
t
,where L is a linear operator of the highest-order derivative which is assumed to be invertible easily, R is the remaining linear operator of order less than L and N is a nonlinear operator, and g
(
t
) is a source term.Taking Laplace transform on both sides of the above equation, we get(2)
L
L
u
t
+
L
R
u
t
+
L
N
u
t
=
L
g
t
.Using the differential property of Laplace transform and using the initial condition, we get(3)
s
n
L
u
t
-
s
n
-
1
u
0
-
s
n
-
2
u
′
0
-
⋯
-
u
n
-
1
0
+
L
R
u
t
+
L
N
u
t
=
L
g
t
,or(4)
L
u
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
-
1
s
n
L
R
u
t
-
1
s
n
L
N
u
t
+
1
s
n
L
g
t
.Now we represent the unknown functions u
(
t
) by an infinite series of the form(5)
u
t
=
∑
n
=
0
∞
u
n
t
.Here the components u
n
(
t
) are usually determined recurrently and the nonlinear operator N
(
u
) can be decomposed into an infinite series of polynomials, given by(6)
N
u
=
∑
n
=
0
∞
A
n
,where A
n are Adomian polynomials of u
0
,
u
1
,
…
,
u
n defined by(7)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
i
=
0
∞
λ
i
u
i
λ
=
0
,
n
=
0,1
,
2
,
…
.Therefore,(8)
L
∑
n
=
0
∞
u
n
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
-
1
s
n
L
R
∑
n
=
1
∞
u
n
t
-
1
s
n
L
∑
n
=
1
∞
A
n
+
1
s
n
L
g
t
.In general, the recursive relation is given by(9)
L
u
0
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
+
1
s
n
L
g
t
,
L
u
n
+
1
t
=
-
1
s
n
L
R
u
n
t
-
1
s
n
L
A
n
.Applying the inverse Laplace transform to both sides of (9), we obtain  u
n, (
n
≥
0
), which is then substituted into (5).For numerical computation, we get the expression as(10)
ϕ
n
t
=
∑
k
=
0
n
u
k
t
,which is nth term approximation of u
(
t
) and the obtained series solution which converges to the exact solution.


### 2.1.2. Algorithm for Finding the Numerical Solution by LADM

Step 1.
Split the given equation into two parts. The first part is(11)
-
R
u
t
+
N
u
tand the second part is(12)
L
u
t
=
g
t
.Step 2.
Apply the Laplace transform to second part, determine the coefficient ofL
u
t, and solve this equation for L
u
t. So, you can get L
u
0.Step 3.
Calculate the Adomian polynomials for the functionN
u
t. Apply the Laplace transform to these polynomials.Step 4.
Divide the first part to the coefficient ofL
u
t. In a loop, calculate L
u
n
+
1
t with substituting L
A
n and L
u
n
t into the first part.Step 5.
Construct the solution using inverse Laplace transform toL
u
n
t.Step 6.
End


## 2.1.1. LADM for Nonlinear Differential Equation

The Laplace Adomian Decomposition Method (LADM) [10, 11] was firstly introduced by Khuri and has been successfully used to find the solution of linear and nonlinear differential equations. The significant advantage of this method is that it is a combination of two powerful techniques, namely, Laplace transform and Adomian Decomposition Method [12–15]. The Laplace transform is an elementary but useful technique for solving linear ordinary differential equations that is widely used by scientists and engineers for tackling the linear models. The main impotence of this method is that the solution of this method is expressed in ∞-series which converges to the exact solution and will not take too much time for computing.Consider the following nonlinear differential equation:(1)
L
u
t
+
R
u
t
+
N
u
t
=
g
t
,where L is a linear operator of the highest-order derivative which is assumed to be invertible easily, R is the remaining linear operator of order less than L and N is a nonlinear operator, and g
(
t
) is a source term.Taking Laplace transform on both sides of the above equation, we get(2)
L
L
u
t
+
L
R
u
t
+
L
N
u
t
=
L
g
t
.Using the differential property of Laplace transform and using the initial condition, we get(3)
s
n
L
u
t
-
s
n
-
1
u
0
-
s
n
-
2
u
′
0
-
⋯
-
u
n
-
1
0
+
L
R
u
t
+
L
N
u
t
=
L
g
t
,or(4)
L
u
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
-
1
s
n
L
R
u
t
-
1
s
n
L
N
u
t
+
1
s
n
L
g
t
.Now we represent the unknown functions u
(
t
) by an infinite series of the form(5)
u
t
=
∑
n
=
0
∞
u
n
t
.Here the components u
n
(
t
) are usually determined recurrently and the nonlinear operator N
(
u
) can be decomposed into an infinite series of polynomials, given by(6)
N
u
=
∑
n
=
0
∞
A
n
,where A
n are Adomian polynomials of u
0
,
u
1
,
…
,
u
n defined by(7)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
i
=
0
∞
λ
i
u
i
λ
=
0
,
n
=
0,1
,
2
,
…
.Therefore,(8)
L
∑
n
=
0
∞
u
n
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
-
1
s
n
L
R
∑
n
=
1
∞
u
n
t
-
1
s
n
L
∑
n
=
1
∞
A
n
+
1
s
n
L
g
t
.In general, the recursive relation is given by(9)
L
u
0
t
=
u
0
s
+
u
′
0
s
2
+
⋯
+
u
n
-
1
0
s
n
+
1
s
n
L
g
t
,
L
u
n
+
1
t
=
-
1
s
n
L
R
u
n
t
-
1
s
n
L
A
n
.Applying the inverse Laplace transform to both sides of (9), we obtain  u
n, (
n
≥
0
), which is then substituted into (5).For numerical computation, we get the expression as(10)
ϕ
n
t
=
∑
k
=
0
n
u
k
t
,which is nth term approximation of u
(
t
) and the obtained series solution which converges to the exact solution.


## 2.1.2. Algorithm for Finding the Numerical Solution by LADM

Step 1.
Split the given equation into two parts. The first part is(11)
-
R
u
t
+
N
u
tand the second part is(12)
L
u
t
=
g
t
.Step 2.
Apply the Laplace transform to second part, determine the coefficient ofL
u
t, and solve this equation for L
u
t. So, you can get L
u
0.Step 3.
Calculate the Adomian polynomials for the functionN
u
t. Apply the Laplace transform to these polynomials.Step 4.
Divide the first part to the coefficient ofL
u
t. In a loop, calculate L
u
n
+
1
t with substituting L
A
n and L
u
n
t into the first part.Step 5.
Construct the solution using inverse Laplace transform toL
u
n
t.Step 6.
End


## 2.2. Runge-Kutta-Fehlberg (RKF) Method

### 2.2.1. RKF for Nonlinear Differential Equation

One of the most popular methods with a constant step size is the fourth-order Runge-Kutta (RK4) method. Reasonably the Runge-Kutta method can [16] obtain the accuracy of a Taylor Series approximation without the need of higher derivative calculations. This method can be considered as the basic form of other methods. However, in terms of error estimation, the one-step method with an adaptive step size like the Runge-Kutta-Fehlberg (RKF) method [17, 18] gives better error estimation than one-step method with a constant step size like the Runge-Kutta method. At each step the Runge-Kutta-Fehlberg method described the calculation of two Runge-Kutta methods of different order (RK4 and RK5). If the two answers are close enough then continue for the next step with the same step size. For a fixed accuracy the step size should be reduced. If the answers agree with more significant digits than required, the step size is increased. So we conclude that the one-step algorithm method associated with an adaptive step size automatically organizes the step size as a recomposition to the calculation truncation errors. This method has shown that it works in the case of nonlinear models and, hence, its application is found in wide range of deterministic and stochastic problems, linear and nonlinear problems, physics, biology, and chemical reactions problems, and so forth.Consider the initial value problem(13)
y
′
t
=
f
t
,
y
t
,
y
t
0
=
y
0
.The RKF is one way to try to resolve this problem.The problem is to solve the initial value problem in the above equation by means of Runge-Kutta methods [19, 20] of order 4 and order 5.First we need some definitions:(14)
k
1
=
h
f
t
i
,
y
i
,
k
2
=
h
f
t
i
+
1
4
h
,
y
i
+
1
4
k
1
,
k
3
=
h
f
t
i
+
3
8
h
,
y
i
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
f
t
i
+
12
13
h
,
y
i
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
f
t
i
+
h
,
y
i
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
f
t
i
+
h
,
y
i
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
.Then an approximation to the solution of initial value problem is made using Runge-Kutta method of order 4:(15)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(16)
z
i
+
1
=
y
i
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.The numerical approximation of the global (local) discretization error at the point t
i
+
1 is(17)
ϵ
=
y
i
+
1
-
z
i
+
1
.If ϵ
m
i
n
≤
ϵ
≤
ϵ
m
a
x, (where the maximum and minimum tolerance for local truncation error are ϵ
m
a
x and ϵ
m
i
n, respectively, which are prescribed at the beginning) then y
i
+
1 is an acceptable approximation of y
(
t
i
+
1
) and the step size is considerable for the next step. Otherwise, the new step size s
h can be determined by multiplying the scalar s times the step size h. The scalar s is given by(18)
s
=
ϵ
h
2
y
i
+
1
-
z
i
+
1
1
/
4
=
0.0840896
ϵ
h
y
i
+
1
-
z
i
+
1
1
/
4
,where ϵ is the specified error control tolerance.Note. RK4 requires four functions for finding the result and RK5 requires six; that is total of ten for RK4 and RK5. Fehlberg devised a method to get RK4 and RK5 results using only six function evaluations by using some of k values in both methods where k
=
∂
f
/
∂
y.


### 2.2.2. Algorithm for Finding the Numerical Solution by RKF

Step 1.
F
(
t
,
y
1
)
← “Function to be supplied”.Step 2.
Readt
(
0
), y
1
(
0
), h, limit, ϵ
m
i
n
,
ϵ
m
a
x.Step 3.
Fori
=
0
(
1
) limit(19)
K
1
⟵
h
f
t
i
,
y
1
t
i
K
2
⟵
h
·
f
t
i
+
1
4
h
,
y
1
t
i
+
1
4
K
1
K
3
⟵
h
·
f
t
i
+
3
8
h
,
y
1
t
i
+
3
32
K
1
+
9
32
K
2
K
4
⟵
h
·
f
t
i
+
12
13
h
,
y
1
t
i
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
K
5
⟵
h
f
t
i
+
h
,
y
1
t
i
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
K
6
⟵
h
f
t
i
+
h
,
y
1
t
i
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
.Step 4.
(20)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.Step 5.
(21)
z
i
+
1
t
i
+
1
=
y
i
t
i
+
16
135
K
1
+
6656
12
,
825
K
3
+
28
,
561
56
,
430
K
4
-
9
50
K
5
+
2
55
K
6
.Step 6.
(22)
ϵ
=
y
i
+
1
-
z
i
+
1
.Step 7.
Considerϵ
m
i
n
≤
ϵ
≤
ϵ
m
a
x.Step 8.
t
i
+
1
=
t
i
+
h. Write z
i
t
i
+
1, t
i.Step 9.
Repeat for better approximation.Step 10.
End the programme.Step 11.
Otherwise,h
=
s
h with s
=
(
ϵ
h
/
2
|
y
i
+
1
-
z
i
+
1
|
)
1
/
4
=
0.0840896
(
ϵ
h
/
|
y
i
+
1
-
z
i
+
1
|
)
1
/
4.Step 12.
Repeat Steps8
–
10 for better accuracy; any one can stop here. For better accuracy the process can be run.


## 2.2.1. RKF for Nonlinear Differential Equation

One of the most popular methods with a constant step size is the fourth-order Runge-Kutta (RK4) method. Reasonably the Runge-Kutta method can [16] obtain the accuracy of a Taylor Series approximation without the need of higher derivative calculations. This method can be considered as the basic form of other methods. However, in terms of error estimation, the one-step method with an adaptive step size like the Runge-Kutta-Fehlberg (RKF) method [17, 18] gives better error estimation than one-step method with a constant step size like the Runge-Kutta method. At each step the Runge-Kutta-Fehlberg method described the calculation of two Runge-Kutta methods of different order (RK4 and RK5). If the two answers are close enough then continue for the next step with the same step size. For a fixed accuracy the step size should be reduced. If the answers agree with more significant digits than required, the step size is increased. So we conclude that the one-step algorithm method associated with an adaptive step size automatically organizes the step size as a recomposition to the calculation truncation errors. This method has shown that it works in the case of nonlinear models and, hence, its application is found in wide range of deterministic and stochastic problems, linear and nonlinear problems, physics, biology, and chemical reactions problems, and so forth.Consider the initial value problem(13)
y
′
t
=
f
t
,
y
t
,
y
t
0
=
y
0
.The RKF is one way to try to resolve this problem.The problem is to solve the initial value problem in the above equation by means of Runge-Kutta methods [19, 20] of order 4 and order 5.First we need some definitions:(14)
k
1
=
h
f
t
i
,
y
i
,
k
2
=
h
f
t
i
+
1
4
h
,
y
i
+
1
4
k
1
,
k
3
=
h
f
t
i
+
3
8
h
,
y
i
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
f
t
i
+
12
13
h
,
y
i
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
f
t
i
+
h
,
y
i
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
f
t
i
+
h
,
y
i
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
.Then an approximation to the solution of initial value problem is made using Runge-Kutta method of order 4:(15)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(16)
z
i
+
1
=
y
i
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.The numerical approximation of the global (local) discretization error at the point t
i
+
1 is(17)
ϵ
=
y
i
+
1
-
z
i
+
1
.If ϵ
m
i
n
≤
ϵ
≤
ϵ
m
a
x, (where the maximum and minimum tolerance for local truncation error are ϵ
m
a
x and ϵ
m
i
n, respectively, which are prescribed at the beginning) then y
i
+
1 is an acceptable approximation of y
(
t
i
+
1
) and the step size is considerable for the next step. Otherwise, the new step size s
h can be determined by multiplying the scalar s times the step size h. The scalar s is given by(18)
s
=
ϵ
h
2
y
i
+
1
-
z
i
+
1
1
/
4
=
0.0840896
ϵ
h
y
i
+
1
-
z
i
+
1
1
/
4
,where ϵ is the specified error control tolerance.Note. RK4 requires four functions for finding the result and RK5 requires six; that is total of ten for RK4 and RK5. Fehlberg devised a method to get RK4 and RK5 results using only six function evaluations by using some of k values in both methods where k
=
∂
f
/
∂
y.


## 2.2.2. Algorithm for Finding the Numerical Solution by RKF

Step 1.
F
(
t
,
y
1
)
← “Function to be supplied”.Step 2.
Readt
(
0
), y
1
(
0
), h, limit, ϵ
m
i
n
,
ϵ
m
a
x.Step 3.
Fori
=
0
(
1
) limit(19)
K
1
⟵
h
f
t
i
,
y
1
t
i
K
2
⟵
h
·
f
t
i
+
1
4
h
,
y
1
t
i
+
1
4
K
1
K
3
⟵
h
·
f
t
i
+
3
8
h
,
y
1
t
i
+
3
32
K
1
+
9
32
K
2
K
4
⟵
h
·
f
t
i
+
12
13
h
,
y
1
t
i
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
K
5
⟵
h
f
t
i
+
h
,
y
1
t
i
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
K
6
⟵
h
f
t
i
+
h
,
y
1
t
i
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
.Step 4.
(20)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.Step 5.
(21)
z
i
+
1
t
i
+
1
=
y
i
t
i
+
16
135
K
1
+
6656
12
,
825
K
3
+
28
,
561
56
,
430
K
4
-
9
50
K
5
+
2
55
K
6
.Step 6.
(22)
ϵ
=
y
i
+
1
-
z
i
+
1
.Step 7.
Considerϵ
m
i
n
≤
ϵ
≤
ϵ
m
a
x.Step 8.
t
i
+
1
=
t
i
+
h. Write z
i
t
i
+
1, t
i.Step 9.
Repeat for better approximation.Step 10.
End the programme.Step 11.
Otherwise,h
=
s
h with s
=
(
ϵ
h
/
2
|
y
i
+
1
-
z
i
+
1
|
)
1
/
4
=
0.0840896
(
ϵ
h
/
|
y
i
+
1
-
z
i
+
1
|
)
1
/
4.Step 12.
Repeat Steps8
–
10 for better accuracy; any one can stop here. For better accuracy the process can be run.


## 3. Solution of Some Population Model by Numerical Techniques

### 3.1. Insect Population Model (Model I)

Suppose that an insect populationP shows the seasonal growth model which is discussed by Erbe et al. [21]. The differential equation of the insect population model is given by(23)
d
P
d
t
=
K
P
cos
⁡
λ
t
,where K and λ are positive constants.We now solve (23) by LADM with initial condition P
(
0
) using the Laplace transform and Adomian Decomposition Method.(24)
∑
n
=
0
∞
P
n
t
=
P
0
+
K
L
-
1
1
s
L
∑
n
=
0
∞
P
n
t
cos
⁡
λ
t
.Now, we have the following recursive algorithm by applying LADM,(25)
P
0
=
P
0
,
P
n
+
1
t
=
K
L
-
1
1
s
L
P
n
t
cos
⁡
λ
t
,
n
≥
0We now solve (23) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemP
′
t
=
f
t
,
P
t
=
K
P
cos
⁡
λ
t
,
P
t
0
=
P
0.First we define(26)
k
1
=
h
K
P
0
cos
⁡
λ
t
0
,
k
2
=
h
K
P
0
+
1
4
k
1
cos
⁡
λ
t
0
+
h
4
,
k
3
=
h
K
P
0
+
3
32
k
1
+
9
32
k
2
cos
⁡
λ
t
0
+
3
8
h
,
k
4
=
h
K
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
cos
⁡
λ
t
0
+
12
13
h
,
k
5
=
h
K
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
cos
⁡
λ
t
0
+
h
,
k
6
=
h
K
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
cos
⁡
λ
t
0
+
h
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(27)
P
n
+
1
=
P
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(28)
Z
n
+
1
=
P
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


### 3.2. Modified Insect Population Model (Model II)

Now we modify the above model by introducing Holling type III functional response which is defined in the following:(29)
d
P
d
t
=
K
P
cos
⁡
λ
t
-
α
P
2
β
+
P
2
.The term α
P
2
/
(
β
+
P
2
) is a Holling type III functional response, where α and β are positive constants and α denotes the maximum capture rate of insects by a predator species.The feeding rate saturates at the maximal feeding rateα. This is intended to model that, for small densities, preys are able to evade predators by taking refuge. Such a functional response is known as sigmoidal. Biologically, a sigmoid functional response explains the fact that, at low densities of prey population, the effect of predation is low, but if the population size increases, predation is then more intensive. This phenomenon appears in a number of interactions in the real world and in this case it is said that the predator is generalist, because if the prey population size is low, it would then seek other food alternatives [22].We now solve (29) by LADM with initial condition N
(
0
) using the Laplace transform and Adomian Decomposition Method.Now applying LADM, we have the following recursive algorithm:(30)
P
0
=
P
0
-
α
t
,
P
n
+
1
t
=
K
L
-
1
1
s
L
P
n
cos
⁡
λ
t
+
α
β
L
-
1
1
s
L
A
n
,
n
≥
0
.Therefore P can be represented as a series P
=
∑
n
=
0
∞
P
n and the Adomian Polynomial is calculating by formula given below:(31)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
k
=
0
n
λ
k
P
k
λ
=
0
.We now solve (29) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemP
′
t
=
f
t
,
P
t
=
K
P
cos
⁡
λ
t
-
α
P
2
/
(
β
+
P
2
),  P
t
0
=
P
0.First we define(32)
k
1
=
h
K
P
0
cos
⁡
λ
t
0
-
α
P
0
2
β
+
P
0
2
,
k
2
=
h
K
P
0
+
1
4
k
1
cos
⁡
λ
t
0
+
h
4
-
α
P
0
+
1
/
4
k
1
2
β
+
P
0
+
1
/
4
k
1
2
,
k
3
=
h
K
P
0
+
3
32
k
1
+
9
32
k
2
cos
⁡
λ
t
0
+
3
8
h
-
α
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
β
+
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
,
k
4
=
h
K
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
cos
⁡
λ
t
0
+
12
13
h
-
α
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
β
+
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
,
k
5
=
h
K
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
cos
⁡
λ
t
0
+
h
-
α
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
β
+
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
,
k
6
=
h
K
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
cos
⁡
λ
t
0
+
h
-
α
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
β
+
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(33)
P
n
+
1
=
P
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(34)
Z
n
+
1
=
P
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


### 3.3. Introducing Interspecific Competition Term in Modified Insect Population Model (Model III)

Here we consider interspecific competition term in the insect population model. Interspecific competition is the most important aspect in the population growth dynamics model [23, 24]. It is assumed to be additional instantaneous deaths to the insect population and the increased death rate is proportional to the square of the predator density. This competition ennoble the death rate of the species; it is expected that the growth rate will considerably differ from the growth rate in a simple insect population model. Hence, we study the insect population model having interspecific competition in the population growth dynamics when either prey or predator or both the prey and predator populations are subject to harvest.(35)
d
P
d
t
=
K
P
cos
⁡
λ
t
-
α
P
2
β
+
P
2
-
μ
P
2
,where μ is the coefficient of interspecific competition term.Solving (35) by LADM yields the following recursive algorithm:(36)
P
0
=
P
0
-
α
t
,
P
n
+
1
t
=
K
L
-
1
1
s
L
P
n
cos
⁡
λ
t
-
μ
L
-
1
1
s
L
P
n
2
+
α
β
L
-
1
1
s
L
1
β
+
P
n
2
,
n
≥
0
.Therefore P can be represented as a series P
=
∑
n
=
0
∞
P
n and the Adomian Polynomial is calculating by formula given below:(37)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
k
=
0
n
λ
k
P
k
λ
=
0
.We now solve (35) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemP
′
t
=
f
t
,
P
t
=
K
P
cos
⁡
λ
t
-
α
P
2
/
(
β
+
P
2
)
-
μ
P
2,  P
t
0
=
P
0.First we define(38)
k
1
=
h
K
P
0
cos
⁡
λ
t
0
-
α
P
0
2
β
+
P
0
2
-
μ
P
0
2
,
k
2
=
h
K
P
0
+
1
4
k
1
cos
⁡
λ
t
0
+
h
4
-
α
P
0
+
1
/
4
k
1
2
β
+
P
0
+
1
/
4
k
1
2
-
μ
P
0
+
1
4
k
1
2
,
k
3
=
h
K
P
0
+
3
32
k
1
+
9
32
k
2
cos
⁡
λ
t
0
+
3
8
h
-
α
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
β
+
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
-
μ
P
0
+
3
32
k
1
+
9
32
k
2
2
,
k
4
=
h
K
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
cos
⁡
λ
t
0
+
12
13
h
-
α
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
β
+
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
-
μ
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
2
,
k
5
=
h
K
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
cos
⁡
λ
t
0
+
h
-
α
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
β
+
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
-
μ
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
2
,
k
6
=
h
K
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
cos
⁡
λ
t
0
+
h
-
α
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
β
+
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
-
μ
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
2
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(39)
P
n
+
1
=
P
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(40)
Z
n
+
1
=
P
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


### 3.4. Lotka-Volterra Model

We define the Lotka-Volterra model into different ways as follows.


#### 3.4.1. Multispecies Lotka-Volterra Model

Mathematical models of population growth have been formed to provide an inconceivable significant angle of true ecological situation. The meaning of each parameter in the models has been defined biologically. Forn species, we consider the following [25, 26] general Lotka-Volterra model: (41)
d
N
i
d
t
=
N
i
b
i
-
∑
j
=
1
n
a
i
j
N
j
,
i
=
1,2
,
…
,
n
.These equations may represent either predator-prey or competition cases.


#### 3.4.2. Single Species Lotka-Volterra Model (Model IV)

In case of one species, (41) is written for a given limited source of food:(42)
d
N
d
t
=
N
b
-
a
N
,
b
>
0
,
a
>
0
,
N
0
>
0
,where a and b are positive constants. This equation has an exact solution(43)
N
t
=
b
e
b
t
b
-
a
N
0
/
N
0
+
a
e
b
t
for
b
≠
0
,
=
N
0
t
1
+
a
t
N
0
t
for
b
=
0
,where N
(
0
) is the initial condition.Solving (42) by LADM yields the following recursive algorithm:(44)
N
0
=
N
0
,
N
n
+
1
t
=
b
L
-
1
1
s
L
N
n
t
-
a
L
-
1
1
s
L
A
1
,
n
,
n
≥
0
,where the Adomian Polynomials A
1
,
n are given by(45)
A
1
,
n
=
∑
i
=
0
n
N
i
N
n
-
i
.We now solve (42) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemN
′
t
=
f
t
,
N
t
=
N
b
-
a
N
,
N
t
0
=
N
0
=
0.1.First we define(46)
k
1
=
h
N
0
b
-
a
N
0
,
k
2
=
h
N
0
+
1
4
k
1
b
-
a
N
0
+
1
4
k
1
,
k
3
=
h
N
0
+
3
32
k
1
+
9
32
k
2
b
-
a
N
0
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
b
-
a
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
b
-
a
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
b
-
a
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(47)
N
n
+
1
=
N
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(48)
Z
n
+
1
=
N
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


### 3.5. Lotka-Volterra Equation with Holling Type III Functional Response (Model V)

Now we modify the Lotka-Volterra model by introducing Holling type III functional response which is defined as(49)
d
N
d
t
=
N
b
-
a
N
-
α
N
2
β
+
N
2
.The term α
N
2
/
(
β
+
N
2
) is a Holling type III functional, where α and β are positive constants and α denotes the predation rate and β denotes half-saturation constant.Solving (49) by LADM yields the following recursive algorithm:(50)
N
0
=
N
0
,
N
n
+
1
t
=
b
L
-
1
1
s
L
N
n
t
-
a
L
-
1
1
s
L
A
2,1
,
n
-
α
L
-
1
1
s
L
A
2,2
,
n
,
n
≥
0
,where the Adomian polynomials A
2
,
k
,
n are given by(51)
A
2
,
k
,
n
=
d
n
d
λ
n
N
∑
i
=
0
n
λ
i
N
i
λ
=
0
,
k
=
0,1
.We now solve (49) by Runge-Kutta-Fehlberg (RKF) method:Consider the initial value problem(52)
N
′
t
=
f
t
,
N
t
=
N
b
-
a
N
-
α
N
2
β
+
N
2
,
N
t
0
=
N
0
=
0.1
,where we know that(53)
k
1
=
h
N
0
b
-
a
N
0
-
α
N
0
2
β
+
N
0
2
,
k
2
=
h
N
0
+
1
4
k
1
b
-
a
N
0
+
1
4
k
1
-
α
N
0
+
1
/
4
k
1
2
β
+
N
0
+
1
/
4
k
1
2
,
k
3
=
h
N
0
+
3
32
k
1
+
9
32
k
2
b
-
a
N
0
+
3
32
k
1
+
9
32
k
2
-
α
N
0
+
3
/
32
k
1
+
9
/
32
k
2
2
β
+
N
0
+
3
/
32
k
1
+
9
/
32
k
2
2
,
k
4
=
h
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
b
-
a
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
-
α
N
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
β
+
N
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
,
k
5
=
h
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
b
-
a
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
-
α
N
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
β
+
N
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
,
k
6
=
h
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
b
-
a
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
-
α
N
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
β
+
N
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(54)
N
n
+
1
=
N
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(55)
Z
n
+
1
=
N
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 3.1. Insect Population Model (Model I)

Suppose that an insect populationP shows the seasonal growth model which is discussed by Erbe et al. [21]. The differential equation of the insect population model is given by(23)
d
P
d
t
=
K
P
cos
⁡
λ
t
,where K and λ are positive constants.We now solve (23) by LADM with initial condition P
(
0
) using the Laplace transform and Adomian Decomposition Method.(24)
∑
n
=
0
∞
P
n
t
=
P
0
+
K
L
-
1
1
s
L
∑
n
=
0
∞
P
n
t
cos
⁡
λ
t
.Now, we have the following recursive algorithm by applying LADM,(25)
P
0
=
P
0
,
P
n
+
1
t
=
K
L
-
1
1
s
L
P
n
t
cos
⁡
λ
t
,
n
≥
0We now solve (23) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemP
′
t
=
f
t
,
P
t
=
K
P
cos
⁡
λ
t
,
P
t
0
=
P
0.First we define(26)
k
1
=
h
K
P
0
cos
⁡
λ
t
0
,
k
2
=
h
K
P
0
+
1
4
k
1
cos
⁡
λ
t
0
+
h
4
,
k
3
=
h
K
P
0
+
3
32
k
1
+
9
32
k
2
cos
⁡
λ
t
0
+
3
8
h
,
k
4
=
h
K
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
cos
⁡
λ
t
0
+
12
13
h
,
k
5
=
h
K
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
cos
⁡
λ
t
0
+
h
,
k
6
=
h
K
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
cos
⁡
λ
t
0
+
h
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(27)
P
n
+
1
=
P
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(28)
Z
n
+
1
=
P
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 3.2. Modified Insect Population Model (Model II)

Now we modify the above model by introducing Holling type III functional response which is defined in the following:(29)
d
P
d
t
=
K
P
cos
⁡
λ
t
-
α
P
2
β
+
P
2
.The term α
P
2
/
(
β
+
P
2
) is a Holling type III functional response, where α and β are positive constants and α denotes the maximum capture rate of insects by a predator species.The feeding rate saturates at the maximal feeding rateα. This is intended to model that, for small densities, preys are able to evade predators by taking refuge. Such a functional response is known as sigmoidal. Biologically, a sigmoid functional response explains the fact that, at low densities of prey population, the effect of predation is low, but if the population size increases, predation is then more intensive. This phenomenon appears in a number of interactions in the real world and in this case it is said that the predator is generalist, because if the prey population size is low, it would then seek other food alternatives [22].We now solve (29) by LADM with initial condition N
(
0
) using the Laplace transform and Adomian Decomposition Method.Now applying LADM, we have the following recursive algorithm:(30)
P
0
=
P
0
-
α
t
,
P
n
+
1
t
=
K
L
-
1
1
s
L
P
n
cos
⁡
λ
t
+
α
β
L
-
1
1
s
L
A
n
,
n
≥
0
.Therefore P can be represented as a series P
=
∑
n
=
0
∞
P
n and the Adomian Polynomial is calculating by formula given below:(31)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
k
=
0
n
λ
k
P
k
λ
=
0
.We now solve (29) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemP
′
t
=
f
t
,
P
t
=
K
P
cos
⁡
λ
t
-
α
P
2
/
(
β
+
P
2
),  P
t
0
=
P
0.First we define(32)
k
1
=
h
K
P
0
cos
⁡
λ
t
0
-
α
P
0
2
β
+
P
0
2
,
k
2
=
h
K
P
0
+
1
4
k
1
cos
⁡
λ
t
0
+
h
4
-
α
P
0
+
1
/
4
k
1
2
β
+
P
0
+
1
/
4
k
1
2
,
k
3
=
h
K
P
0
+
3
32
k
1
+
9
32
k
2
cos
⁡
λ
t
0
+
3
8
h
-
α
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
β
+
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
,
k
4
=
h
K
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
cos
⁡
λ
t
0
+
12
13
h
-
α
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
β
+
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
,
k
5
=
h
K
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
cos
⁡
λ
t
0
+
h
-
α
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
β
+
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
,
k
6
=
h
K
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
cos
⁡
λ
t
0
+
h
-
α
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
β
+
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(33)
P
n
+
1
=
P
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(34)
Z
n
+
1
=
P
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 3.3. Introducing Interspecific Competition Term in Modified Insect Population Model (Model III)

Here we consider interspecific competition term in the insect population model. Interspecific competition is the most important aspect in the population growth dynamics model [23, 24]. It is assumed to be additional instantaneous deaths to the insect population and the increased death rate is proportional to the square of the predator density. This competition ennoble the death rate of the species; it is expected that the growth rate will considerably differ from the growth rate in a simple insect population model. Hence, we study the insect population model having interspecific competition in the population growth dynamics when either prey or predator or both the prey and predator populations are subject to harvest.(35)
d
P
d
t
=
K
P
cos
⁡
λ
t
-
α
P
2
β
+
P
2
-
μ
P
2
,where μ is the coefficient of interspecific competition term.Solving (35) by LADM yields the following recursive algorithm:(36)
P
0
=
P
0
-
α
t
,
P
n
+
1
t
=
K
L
-
1
1
s
L
P
n
cos
⁡
λ
t
-
μ
L
-
1
1
s
L
P
n
2
+
α
β
L
-
1
1
s
L
1
β
+
P
n
2
,
n
≥
0
.Therefore P can be represented as a series P
=
∑
n
=
0
∞
P
n and the Adomian Polynomial is calculating by formula given below:(37)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
k
=
0
n
λ
k
P
k
λ
=
0
.We now solve (35) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemP
′
t
=
f
t
,
P
t
=
K
P
cos
⁡
λ
t
-
α
P
2
/
(
β
+
P
2
)
-
μ
P
2,  P
t
0
=
P
0.First we define(38)
k
1
=
h
K
P
0
cos
⁡
λ
t
0
-
α
P
0
2
β
+
P
0
2
-
μ
P
0
2
,
k
2
=
h
K
P
0
+
1
4
k
1
cos
⁡
λ
t
0
+
h
4
-
α
P
0
+
1
/
4
k
1
2
β
+
P
0
+
1
/
4
k
1
2
-
μ
P
0
+
1
4
k
1
2
,
k
3
=
h
K
P
0
+
3
32
k
1
+
9
32
k
2
cos
⁡
λ
t
0
+
3
8
h
-
α
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
β
+
P
0
+
3
/
32
k
1
+
9
/
32
k
2
2
-
μ
P
0
+
3
32
k
1
+
9
32
k
2
2
,
k
4
=
h
K
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
cos
⁡
λ
t
0
+
12
13
h
-
α
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
β
+
P
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
-
μ
P
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
2
,
k
5
=
h
K
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
cos
⁡
λ
t
0
+
h
-
α
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
β
+
P
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
-
μ
P
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
2
,
k
6
=
h
K
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
cos
⁡
λ
t
0
+
h
-
α
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
β
+
P
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
-
μ
P
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
2
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(39)
P
n
+
1
=
P
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(40)
Z
n
+
1
=
P
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 3.4. Lotka-Volterra Model

We define the Lotka-Volterra model into different ways as follows.


### 3.4.1. Multispecies Lotka-Volterra Model

Mathematical models of population growth have been formed to provide an inconceivable significant angle of true ecological situation. The meaning of each parameter in the models has been defined biologically. Forn species, we consider the following [25, 26] general Lotka-Volterra model: (41)
d
N
i
d
t
=
N
i
b
i
-
∑
j
=
1
n
a
i
j
N
j
,
i
=
1,2
,
…
,
n
.These equations may represent either predator-prey or competition cases.


### 3.4.2. Single Species Lotka-Volterra Model (Model IV)

In case of one species, (41) is written for a given limited source of food:(42)
d
N
d
t
=
N
b
-
a
N
,
b
>
0
,
a
>
0
,
N
0
>
0
,where a and b are positive constants. This equation has an exact solution(43)
N
t
=
b
e
b
t
b
-
a
N
0
/
N
0
+
a
e
b
t
for
b
≠
0
,
=
N
0
t
1
+
a
t
N
0
t
for
b
=
0
,where N
(
0
) is the initial condition.Solving (42) by LADM yields the following recursive algorithm:(44)
N
0
=
N
0
,
N
n
+
1
t
=
b
L
-
1
1
s
L
N
n
t
-
a
L
-
1
1
s
L
A
1
,
n
,
n
≥
0
,where the Adomian Polynomials A
1
,
n are given by(45)
A
1
,
n
=
∑
i
=
0
n
N
i
N
n
-
i
.We now solve (42) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemN
′
t
=
f
t
,
N
t
=
N
b
-
a
N
,
N
t
0
=
N
0
=
0.1.First we define(46)
k
1
=
h
N
0
b
-
a
N
0
,
k
2
=
h
N
0
+
1
4
k
1
b
-
a
N
0
+
1
4
k
1
,
k
3
=
h
N
0
+
3
32
k
1
+
9
32
k
2
b
-
a
N
0
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
b
-
a
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
b
-
a
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
b
-
a
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(47)
N
n
+
1
=
N
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(48)
Z
n
+
1
=
N
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 3.4.1. Multispecies Lotka-Volterra Model

Mathematical models of population growth have been formed to provide an inconceivable significant angle of true ecological situation. The meaning of each parameter in the models has been defined biologically. Forn species, we consider the following [25, 26] general Lotka-Volterra model: (41)
d
N
i
d
t
=
N
i
b
i
-
∑
j
=
1
n
a
i
j
N
j
,
i
=
1,2
,
…
,
n
.These equations may represent either predator-prey or competition cases.


## 3.4.2. Single Species Lotka-Volterra Model (Model IV)

In case of one species, (41) is written for a given limited source of food:(42)
d
N
d
t
=
N
b
-
a
N
,
b
>
0
,
a
>
0
,
N
0
>
0
,where a and b are positive constants. This equation has an exact solution(43)
N
t
=
b
e
b
t
b
-
a
N
0
/
N
0
+
a
e
b
t
for
b
≠
0
,
=
N
0
t
1
+
a
t
N
0
t
for
b
=
0
,where N
(
0
) is the initial condition.Solving (42) by LADM yields the following recursive algorithm:(44)
N
0
=
N
0
,
N
n
+
1
t
=
b
L
-
1
1
s
L
N
n
t
-
a
L
-
1
1
s
L
A
1
,
n
,
n
≥
0
,where the Adomian Polynomials A
1
,
n are given by(45)
A
1
,
n
=
∑
i
=
0
n
N
i
N
n
-
i
.We now solve (42) by Runge-Kutta-Fehlberg (RKF) method.Consider the initial value problemN
′
t
=
f
t
,
N
t
=
N
b
-
a
N
,
N
t
0
=
N
0
=
0.1.First we define(46)
k
1
=
h
N
0
b
-
a
N
0
,
k
2
=
h
N
0
+
1
4
k
1
b
-
a
N
0
+
1
4
k
1
,
k
3
=
h
N
0
+
3
32
k
1
+
9
32
k
2
b
-
a
N
0
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
b
-
a
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
b
-
a
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
b
-
a
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(47)
N
n
+
1
=
N
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined by using a Runge-Kutta method of order 5:(48)
Z
n
+
1
=
N
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 3.5. Lotka-Volterra Equation with Holling Type III Functional Response (Model V)

Now we modify the Lotka-Volterra model by introducing Holling type III functional response which is defined as(49)
d
N
d
t
=
N
b
-
a
N
-
α
N
2
β
+
N
2
.The term α
N
2
/
(
β
+
N
2
) is a Holling type III functional, where α and β are positive constants and α denotes the predation rate and β denotes half-saturation constant.Solving (49) by LADM yields the following recursive algorithm:(50)
N
0
=
N
0
,
N
n
+
1
t
=
b
L
-
1
1
s
L
N
n
t
-
a
L
-
1
1
s
L
A
2,1
,
n
-
α
L
-
1
1
s
L
A
2,2
,
n
,
n
≥
0
,where the Adomian polynomials A
2
,
k
,
n are given by(51)
A
2
,
k
,
n
=
d
n
d
λ
n
N
∑
i
=
0
n
λ
i
N
i
λ
=
0
,
k
=
0,1
.We now solve (49) by Runge-Kutta-Fehlberg (RKF) method:Consider the initial value problem(52)
N
′
t
=
f
t
,
N
t
=
N
b
-
a
N
-
α
N
2
β
+
N
2
,
N
t
0
=
N
0
=
0.1
,where we know that(53)
k
1
=
h
N
0
b
-
a
N
0
-
α
N
0
2
β
+
N
0
2
,
k
2
=
h
N
0
+
1
4
k
1
b
-
a
N
0
+
1
4
k
1
-
α
N
0
+
1
/
4
k
1
2
β
+
N
0
+
1
/
4
k
1
2
,
k
3
=
h
N
0
+
3
32
k
1
+
9
32
k
2
b
-
a
N
0
+
3
32
k
1
+
9
32
k
2
-
α
N
0
+
3
/
32
k
1
+
9
/
32
k
2
2
β
+
N
0
+
3
/
32
k
1
+
9
/
32
k
2
2
,
k
4
=
h
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
b
-
a
N
0
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
-
α
N
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
β
+
N
0
+
1932
/
2197
k
1
-
7200
/
2197
k
2
+
7296
/
2197
k
3
2
,
k
5
=
h
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
b
-
a
N
0
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
-
α
N
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
β
+
N
0
+
439
/
216
k
1
-
8
k
2
+
3680
/
513
k
3
-
845
/
4104
k
4
2
,
k
6
=
h
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
b
-
a
N
0
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
-
α
N
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
β
+
N
0
-
8
/
27
k
1
+
2
k
2
-
3544
/
2565
k
3
+
1859
/
4104
k
4
-
11
/
40
k
5
2
.Then an approximation to the solution of initial value problem is made by using Runge-Kutta method of order 4:(54)
N
n
+
1
=
N
n
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(55)
Z
n
+
1
=
N
n
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.


## 4. Results and Discussion

The numerical solutions obtained by using the RKF, LADM method and are compared with the exact solution for different population models.Remark 1.
From Table1 and Figure 1, we see the comparison among the RKF method, three-iterate LADM, and exact solutions for model I for particular numerical value of parameters and initial condition. The numerical results show that RKF method is of good accuracy.Table 1
Numerical comparison when initial condition isP
(
0
)
=
1000 and K
=
2, λ
=
π, h
=
0.1.
t
Exact solution
Solution by RKF method
Solution by LADM(three iterates)
ERKF
ELADM
0
1000.0000000
1000.0000000
1000.00000000
0.00
E
+
00
0.00
E
+
00
0.1
1217.4076457
1217.4447223
1216.07381402
3.71
E
-
02
1.33
E
+
00
0.2
1453.7920776
1453.8275070
1444.17894316
3.54
E
-
02
9.61
E
+
00
0.3
1673.5897680
1673.6084290
1647.56834151
1.87
E
-
02
2.60
E
+
01
0.4
1831.8331881
1831.8125036
1788.52160650
2.07
E
-
02
4.33
E
+
01
0.5
1889.5969624
1889.5159617
1838.84276708
8.10
E
-
02
5.08
E
+
01
0.6
1831.3776006
1831.2908389
1788.12230677
8.68
E
-
02
4.33
E
+
01
0.7
1672.7981525
1672.7809816
1646.85165910
1.72
E
-
02
2.59
E
+
01
0.8
1452.8457426
1452.9182247
1443.28432644
7.25
E
-
02
9.56
E
+
00
0.9
1216.4762073
1216.5807034
1215.15811427
1.04
E
-
01
1.32
E
+
00
1
999.1956876
999.2604905
999.19564898
6.48
E
-
02
3.86
E
-
05
ELADM→ Error term compared to the exact solution on LADM.
ERKF→ Error term compared to the exact solution on RKF.Figure 1
Comparison among the exact solution and the solutions obtained by using RKF method and LADM for model I.Remark 2.
In Table2 and Figure 2, we show the comparison among the RKF method, three-iterate LADM, and exact solutions for model II in the case when α
=
0.5, β
=
0.03, K
=
2, λ
=
π, h
=
0.05, and P
(
0
)
=
100. Again, the numerical results show that RKF method is of good accuracy.Table 2
Numerical comparison of solutions whenP
(
0
)
=
100, α
=
0.5, β
=
0.03, K
=
2,  λ
=
π
,
h
=
0
.
05.
t
Exact solution
Solution by RKF method
Solution by LADM (three iterates)
ERKF
ELADM
0
100.00000000
100.0000000
100.0000000
0.00
E
+
00
0.00
E
+
00
0.05
110.44538789
110.1789856
110.4297547
2.66
E
-
01
1.56
E
-
02
0.1
121.68557418
121.108106
121.5572303
5.77
E
-
01
1.28
E
-
01
0.15
133.42386083
132.497513
133.0020711
9.26
E
-
01
4.22
E
-
01
0.2
145.25872028
143.9584501
144.3172756
1.30
E
+
00
9.41
E
-
01
0.25
156.69395306
155.0138903
155.0163121
1.68
E
+
00
1.68
E
+
00
0.3
167.16677455
165.1252092
164.6053933
2.04
E
+
00
2.56
E
+
00
0.35
176.09213446
173.7346800
172.6181278
2.36
E
+
00
3.47
E
+
00
0.4
182.92083874
180.3198596
178.6494928
2.60
E
+
00
4.27
E
+
00
0.45
187.20173364
184.4522441
182.3860869
2.75
E
+
00
4.82
E
+
00
0.5
188.63838966
185.8501193
183.6299329
2.79
E
+
00
5.01
E
+
00
0.55
187.1284265
184.4153954
182.3136691
2.71
E
+
00
4.81
E
+
00
0.6
182.7773689
180.2468100
178.5057434
2.53
E
+
00
4.27
E
+
00
0.65
175.8845511
173.6266764
172.4051363
2.26
E
+
00
3.48
E
+
00
0.7
166.9030599
164.9840041
164.3260955
1.92
E
+
00
2.58
E
+
00
0.75
156.3834523
154.8416109
154.6742722
1.54
E
+
00
1.71
E
+
00
0.8
144.9109052
143.7574377
143.9164278
1.15
E
+
00
9.94
E
-
01
0.85
133.047794
132.2701488
132.5464437
7.78
E
-
01
5.01
E
-
01
0.9
121.2890397
120.8566363
121.0506779
4.32
E
-
01
2.38
E
-
01
0.95
110.0347113
109.90535847
109.8757207
1.29
E
-
01
1.59
E
-
01
1
99.57968623
99.70572240
99.4013312
1.26
E
-
01
1.78
E
-
01Figure 2
Evaluation among the exact solution and the solutions obtained by using LADM and RKF method for model II.Remark 3.
In Table3 and Figure 3, we show the comparison among the RKF method, three-iterate LADM, and exact solutions for model III in the case when α
=
0.5, β
=
0.03, K
=
2, λ
=
π, h
=
0.05, μ
=
0.00001, and P
(
0
)
=
100. We conclude that RKF method is of good accuracy.Table 3
Numerical comparison whenP
(
0
)
=
100, α
=
0.5, β
=
0.03, K
=
2, λ
=
π, h
=
0.05,  μ
=
0.00001.
t
Exact solution
Solution by RKF method
Solution by LADM (three iterates)
ERKF
ELADM
0
100.00000
100.00000
100.0000000
0.00
E
+
00
0.00
E
+
00
0.05
110.43960
110.43957
110.4247550
3.00
E
-
05
1.48
E
-
02
0.1
121.67210
121.67914
121.5472313
7.04
E
-
03
1.25
E
-
01
0.15
133.40060
133.41682
132.9870734
1.62
E
-
02
4.14
E
-
01
0.2
145.22320
145.25102
144.2972796
2.77
E
-
02
9.26
E
-
01
0.25
156.64390
156.68567
154.9913184
4.18
E
-
02
1.65
E
+
00
0.3
167.09980
167.15791
164.5754023
5.81
E
-
02
2.52
E
+
00
0.35
176.00650
176.08282
172.5831401
7.63
E
-
02
3.42
E
+
00
0.4
182.81540
182.91116
178.6095088
9.58
E
-
02
4.21
E
+
00
0.45
187.07650
187.19184
182.3411071
1.15
E
-
01
4.74
E
+
00
0.5
188.49450
188.62843
183.5799579
1.34
E
-
01
4.91
E
+
00
0.55
186.96800
187.11853
182.2586994
1.50
E
-
01
4.71
E
+
00
0.6
182.60380
182.76773
178.4457794
1.64
E
-
01
4.16
E
+
00
0.65
175.70170
175.87525
172.3401786
1.73
E
-
01
3.36
E
+
00
0.7
166.71530
166.89424
164.2561445
1.79
E
-
01
2.46
E
+
00
0.75
156.19490
156.37518
154.5993285
1.80
E
-
01
1.60
E
+
00
0.8
144.72520
144.90324
143.8364918
1.78
E
-
01
8.89
E
-
01
0.85
132.86810
133.04075
132.4615159
1.73
E
-
01
4.07
E
-
01
0.9
121.11750
121.28263
120.9607589
1.65
E
-
01
1.57
E
-
01
0.95
109.87270
110.02889
109.7808110
1.56
E
-
01
9.19
E
-
02
1
99.42780
99.57443
99.3014312
1.47
E
-
01
1.26
E
-
01Figure 3
Evaluation among the exact solution and the solutions obtained by using LADM and RKF method for model III.Remark 4.
From Table4 and Figure 4, it shows comparison among the RKF method, three-iterate LADM, and the exact solution for the single species Lotka-Volterra model in the case when b
=
1, a
=
3, and N
(
0
)
=
0.1, h
=
0.1. The results show error-free calculation between exact solution and RKF method whereas there are some amount of error in the calculation between exact solution and LADM.Table 4
Numerical comparison whenN
0
=
0.1, a
=
3, b
=
1, h
=
0.1.
t
Exact solution
Solution by RKF method
Solution by LADM method (three iterates)
ERKF
ELADM
0
0.10000000
0.10000000
0.10000000
0.00
E
+
00
0.00
E
+
00
0.1
0.10713679
0.10713679
0.10714000
0.00
E
+
00
0.00
E
+
00
0.2
0.11453291
0.11453291
0.11456000
0.00
E
+
00
0.00
E
+
00
0.3
0.12216385
0.12216385
0.12226000
0.00
E
+
00
9.61
E
-
05
0.4
0.13000114
0.13000114
0.13024000
0.00
E
+
00
2.39
E
-
04
0.5
0.13801261
0.13801261
0.13850000
0.00
E
+
00
4.87
E
-
04
0.6
0.14616290
0.14616290
0.14704000
0.00
E
+
00
8.77
E
-
04
0.7
0.15441399
0.15441399
0.15586000
0.00
E
+
00
1.45
E
-
03
0.8
0.16272591
0.16272591
0.16496000
0.00
E
+
00
2.23
E
-
03
0.9
0.17105750
0.17105750
0.17434000
0.00
E
+
00
3.28
E
-
03
1
0.17936718
0.17936718
0.17936718
0.00
E
+
00
4.63
E
-
03
1.1
0.18761383
0.18761383
0.19394000
0.00
E
+
00
6.33
E
-
03
1.2
0.19575756
0.19575756
0.20416000
1.00
E
-
10
8.40
E
-
03
1.3
0.20376050
0.20376050
0.21466000
0.00
E
+
00
1.09
E
-
02
1.4
0.21158743
0.21158743
0.22544000
1.00
E
-
10
1.39
E
-
02
1.5
0.21920638
0.21920638
0.23650000
0.00
E
+
00
1.73
E
-
02
1.6
0.22658907
0.22658907
0.24784000
0.00
E
+
00
2.13
E
-
02
1.7
0.23371122
0.23371122
0.25946000
0.00
E
+
00
2.57
E
-
02
1.8
0.24055276
0.24055276
0.27136000
0.00
E
+
00
3.08
E
-
02
1.9
0.24709782
0.24709782
0.28354000
0.00
E
+
00
3.64
E
-
02
2
0.25333471
0.25333471
0.29600000
0.00
E
+
00
4.27
E
-
02Figure 4
Evaluation between the exact solution and the solutions obtained by using LADM and RKF method for model IV.Remark 5.
From Table5 and Figure 5, we show the comparison between the RKF, two-iterate LADM, and exact solutions for model V in the case when N
(
0
)
=
0.1, a
=
3, b
=
1, α
=
0.5, β
=
1, and h
=
0.1. Again, the numerical results show that RKF is of high accuracy. The graphical representations of the 5 models reveal that the exact solution and RKF are overlapping each other whereas there is a little difference between exact solution and that of LADM. So from Figures 1–5, we can also say that RKF is a reliable numerical technique.Table 5
Numerical comparison whenN
(
0
)
=
0
.
1, a
=
3, b
=
1, α
=
0.5, β
=
1, h
=
0.1.
t
Exact solution
Solution by RKF method
Solution by LADM (two iterates)
ERKF
ELADM
0
0.1000000000
0.1000000000
0.1000000000
0.00
E
+
00
0.00
E
+
00
0.1
0.1065991220
0.1065991220
0.1065049505
0.00
E
+
00
9.42
E
-
05
0.2
0.1133693593
0.1133693590
0.1130099010
3.00
E
-
10
3.59
E
-
04
0.3
0.1202826611
0.1202826611
0.1195148515
0.00
E
+
00
7.68
E
-
04
0.4
0.1273083819
0.1273083815
0.1260198020
4.00
E
-
10
1.29
E
-
03
0.5
0.1344137893
0.1344137893
0.1325247525
0.00
E
+
00
1.89
E
-
03
0.6
0.1415646647
0.1415646643
0.1390297030
4.00
E
-
10
2.53
E
-
03
0.7
0.1487259514
0.1487259514
0.1455346535
0.00
E
+
00
3.19
E
-
03
0.8
0.1558624471
0.1558624468
0.1520396040
3.00
E
-
10
3.82
E
-
03
0.9
0.1629394898
0.1629394898
0.1585445545
0.00
E
+
00
4.39
E
-
03
1
0.1699236292
0.1699236290
0.1650495050
2.00
E
-
10
4.87
E
-
03
1.1
0.1767832382
0.1767832382
0.1715544554
0.00
E
+
00
5.23
E
-
03
1.2
0.1834890574
0.1834890572
0.1780594059
2.00
E
-
10
5.43
E
-
03
1.3
0.1900146405
0.1900146404
0.1845643564
1.00
E
-
10
5.45
E
-
03
1.4
0.1963367003
0.1963367001
0.1910693069
2.00
E
-
10
5.27
E
-
03
1.5
0.2024353397
0.2024353396
0.1975742574
1.00
E
-
10
4.86
E
-
03
1.6
0.2082941773
0.2082941768
0.2040792079
5.00
E
-
10
4.21
E
-
03
1.7
0.2139003623
0.2139003624
0.2105841584
1.00
E
-
10
3.32
E
-
03
1.8
0.2192445053
0.2192445047
0.2170891089
6.00
E
-
10
2.16
E
-
03
1.9
0.2243205144
0.2243205144
0.2235940594
0.00
E
+
00
7.26
E
-
04
2
0.2291253839
0.2291253840
0.2300990099
1.00
E
-
10
9.74
E
-
04Figure 5
Evaluation between the exact solution and the solutions obtained by using LADM and RKF method for model V.


## 5. Conclusion

In this paper, we describe the method for finding numerical solution of insect population model and Lotka-Volterra model. Here we apply two numerical methods called RKF and LADM for solutions of the said models. Here the numerical solutions obtained by using the RKF show high accuracy and these are compared with the LADM solution. So we can say that these numerical results show that the RKF method is an acceptable and reliable numerical technique for the solution of linear and nonlinear differential equation models on population models. It can be seen clearly from the graphical representations that RKF gives quite good results after a certain considerable time intervals. This is a very useful method, which will be undoubtedly found applicable in broad applications. The advantage of the RKF over the LADM is that there is no need for the evaluations of the Adomian polynomials and the advantage of RKF over RK4 is that it has a good accuracy using variable step size. Hence it provides an efficient numerical solution.



---

*Source: 8921710-2016-12-07.xml*