# Exponential Time Differencing Method for Studying Prey-Predator Dynamic during Mating Period

Authors: Noufe H. Aljahdaly; H. A. Ashi
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

This paper addresses a first numerical simulation to the nonlinear dynamic system of equations that describes the prey-predator model at the predator mating period. Some male species accompany the females during the mating period. In this case, both male and female feed on the same prey. The presented work shows the numerical solution for this specific case of the prey-predator mathematical model via an exponential time differencing method. In addition, the paper provides the biological implication of the solution.


---

## Body

## 1. Introduction

Many applications and phenomena in physics, biology, the medical field, and other fields can be described as a mathematical model in order to understand or predict the dynamic of these phenomena or applications [1–3]. The prey-predator model is a famous biomathematical model that has been investigated widely to study the dynamic of the ecology system and an intersection between prey and predator [4]. The Lotka-Volterra model is known as the predator-prey model for the study population process [5–7]. The model has been studied when considering the influence of diseases [8], prey refuge [9], food chain [10], or competitive populations [11]. The model also was developed to study the intersection between two different populations of predators and one prey [12, 13]. Recently, the model was modified to study the intersection between two predators from the same species (female and male) and one prey. In some species, the male stays with the female all the time, and they feed together at the mating period. This model is given as follows [14]:
(1)ut=a1u−a2u2−a3uv−a4uv2,vt=D2Δv+a3uv−a5v−a6v2+a4uv2,where u and v are the total numbers of prey and predator, respectively. a1 is the prey’s growth rate. a2, a3, and a4 are the prey’s decay rate due to competition on food supply between the male and the female, between one predator and one prey, and between two predators and one prey, respectively. a5 is the mortality rate for the predator. a6 is the decay rate of the predator due to competition on the food supply between the male and the female. D2 is the predator’s diffusion term. In Reference [14], the system (1) was rescaled and became a dimensionless nonlinear system as follows:
(2)Ut=U−U2−UV−UV2,Vt=VXX−k2k1V−k1V2+k1VU+k1UV2,where k1=a4a1/a2=a3/a2 and k2=a5/a1k1.The mathematicians put forth more effort in the computational field to improve different techniques to be able to find the consistent solutions for realistic applications, for example, the geometric-qualitative method, semiapproximate method, qualitative method, numerical method, and analytical method. Analytical methods generate different solutions based on many factors such as used ansatz. There are many methods to obtain the analytical solutions such as theG′/G2-expansion approach [15] and the direct algebraic method [16–19]. They are able to obtain the solutions without an initial condition. Each method can introduce new solutions to the system. The solutions of system (2) were found by two analytical methods, the G′/G2-expansion approach and the direct algebraic method [15]. Since the analytical methods give different solutions, it is difficult to determine which solution predicts the applications. On the other hand, the numerical method for the initial value problem is the method that finds the approximate solution and converges to only one solution based on the initial condition that is obtained from the application, for example, the finite difference method [20–22], Adomian decomposition method (ADM) [23], Laplacian Adomian decomposition method (LADM), multistage differential transform method [24, 25], and exponential time differencing method [26]. The numerical solution for system (2) will be introduced in this paper via the exponential time differencing method since the method is fourth-order accurate and is not local convergent comparing to ADM or LADM. The results either determine which analytical solution that predicts the problem or find a new solution.The paper is organized as follows: the next section presents the numerical results and the scheme of the used method, the third section discusses the numerical solution from the biological perspective, the fourth section is a comparison with the previous studies for the aforementioned model, and the last section is the summary of our analysis and results.


## 2. Numerical Simulation

In this section, we give the approximate numerical solution of system (2). First, the system is discretized in space utilizing the Fourier spectral method. Hence, we obtain the following system of ordinary differential equations (ODEs) that depends only on time:
(3a)du^ktdt=u^kt+F^ku^kt,v^kt,t,(3b)dv^ktdt=cv^kt+G^ku^kt,v^kt,t,where the parameter c=k2−k2k1 is associated to the linear part, k represents the wave numbers, both terms
(4)F^ku^kt,v^kt,t=fft−U2t−UtVt−UtV2t,G^ku^kt,v^kt,t=k1fft−V2t+UtVt+UtV2t,represent the nonlinear forcing term, and fft is the MATLAB command that represents the fast Fourier transform (FFT).Second, we integrate the obtained uncoupled system of ODEs (3) in time with the aid of the fourth-order exponential time differencing method ETD4RK [27–29]. Applying the ETD4RK method to (3b) takes the formula
(5)ak,n=v^k,necΔt/2+ecΔt/2−1G^k,nc,bk,n=v^k,necΔt/2+ecΔt/2−1G^ak,n,tn+Δt/2c,Ck,n=ak,necΔt/2+ecΔt/2−12G^bk,n,tn+Δt/2−G^k,nc,(6)v^k,n+1=v^k,necΔt+c2Δt2−3cΔt+4ecΔt−cΔt−4G^k,n+2cΔt−2ecΔt+cΔt+2G^ak,n,tn+Δt2+G^bk,n,tn+Δt2+−cΔt+4ecΔt−c2Δt2−3cΔt−4G^Ck,n,tn+Δtc3Δt2,where Δt represents the time step (a similar formula is obtained when applying the ETD4RK method to (3a) with c=1).The time integration starts with initial conditions (taken from [30]) given by
(7)u^0=635−2×10−7x−225x−675,v^0=116245−3×10−5x−450+180×10−4.The results of our experiment are shown in figures.The solution is stable at0≤k2≤1+5/2 based on the stability analysis in Reference [14]. A numerical solution is produced in Figure 1 for different value of k2. We realized that solutions are changed with k2. When k2=1.5,U increases while V decreases until they reach the steady state. For k2=0.1, the solution is periodic. The parameter k2=a2a5/a1a4a1=a5a2/a1a3 has direct relation with a5 and a2 and negative relation with a1,a3, and a4.Figure 1
Plot ofU,V for k1=1 and 0<k2<1+5/2.


## 3. Biological Relevance of the Solution

A new term in the studied model isa4v2u. Thus, the numerical solution is used to study the effect of this term on a traditional prey-predator model. The solution of the traditional model is an undamped oscillation, and Lotka explained this behavior because of the law of mass action [5], while the presented model (2) has an oscillation solution with small k2 at a bounded space and the solution converges to the steady state out of this space, because a4v2u is related to the interaction between prey and two predators at the mating period which disappears after this period and also happens in a specific location where the predators give birth.The large rate ofa4 indicates that the predators feed on a large amount of prey in this period, which lead to a large size of the predator population than the prey population. This is not an optimal dynamic of population to live normally (see Figure 2 when a4=10). Therefore, the predator population will decrease because of absence of a food source, while the prey population will grow again. The new term a4v2u has more effect than a3uv to decay the prey population as we see by comparing Figures 3 and 2.Figure 2
Plot ofu,v for a1=a2=1,a5=0.5 to capture the effect of the term a4uv2 on the model.
(a)
a4=10(b)
a4=5(c)
a4=1(d)
a4=0.4Figure 3
Plot ofu,v for a1=a2=1,a5=0.5 to capture the effect of the term a3uv on the model.
(a)
a3=5(b)
a3=3(c)
a3=1(d)
a3=0.4


## 4. Connection with the Previous Studies

The analytical solutions of Equations (2) was found by Aljahdaly and Alqudah [31] using G′/G and the generalized auxiliary equation [15]. The solutions are Ux,t and Vx,t which denote the population of the prey and predator, respectively, in a single patch. The solutions predict the two following behaviors


### 4.1. TheG′/G-Expansion Method (Solution (1))

As a result of two predators’ sexual interaction and the high prey consumption during the predator’s mating period,V increases and U decreases. Ultimately, V will become dominant in the patch.


### 4.2. Generalized Auxiliary Equation Method (Solution (2))

The predators graze in the mating period where the prey density (U) is large. Then, U decrease and V increase. At a certain time, U will be very low, and the predators will either die or move to different patches. Then, the prey population will recover and grow again in the same location.


### 4.3. Exponential Time Differencing Method

For a large value ofk2, we obtain the numerical solution that converges to solution (1), while for small k2, the numerical solution approximates solution (2).The numerical solution gives the important reasons behind the two different solutions by the analytical methods.


## 4.1. TheG′/G-Expansion Method (Solution (1))

As a result of two predators’ sexual interaction and the high prey consumption during the predator’s mating period,V increases and U decreases. Ultimately, V will become dominant in the patch.


## 4.2. Generalized Auxiliary Equation Method (Solution (2))

The predators graze in the mating period where the prey density (U) is large. Then, U decrease and V increase. At a certain time, U will be very low, and the predators will either die or move to different patches. Then, the prey population will recover and grow again in the same location.


## 4.3. Exponential Time Differencing Method

For a large value ofk2, we obtain the numerical solution that converges to solution (1), while for small k2, the numerical solution approximates solution (2).The numerical solution gives the important reasons behind the two different solutions by the analytical methods.


## 5. Conclusion

The prey-predator system with the extra terma4v2u (1) is a new nonlinear dynamic system of equations in the literature. It has been solved analytically with two different methods which give two different solutions. The importance of this paper is that the numerical solution gives connection between the obtained analytical solutions from the biological perspective. The numerical solution predicts the dynamic of the studied model as well as proves the effect of the new extra term on the model.



---

*Source: 2819145-2021-09-21.xml*