Find the Equilibrium Solution of the Modified Predator-prey System

Prey Population

Solution of Differential Equations

George Lindfield , John Penny , in Numerical Methods (Fourth Edition), 2019

5.10.2 The Predator–Prey Problem

A system of differential equations which models the interaction of competing or predator–prey populations is based on the Volterra equations and may be written in the form

(5.28) $\begin{array}{c}dP/dt=K_{1}P-CPQ\hfill \\ dQ/dt=-K_{2}Q+DPQ\hfill \end{array}$

together with the initial conditions

$Q=Q_0\text{ and}P=P_0\text{ at time}t=0$

The variables P and Q give the size of the prey and predator populations, respectively, at time t. These two populations interact and compete. $K_1$ , $K_2$ , C, and D are positive constants. $K_1$ relates to the rate of growth of the prey population P, and $K_2$ relates to the rate of decay of the predator population Q. It seems reasonable to assume that the number of encounters of predator and prey is proportional to P multiplied by Q and that a proportion C of these encounters will be fatal to members of the prey population. Thus the term CPQ gives a measure of the decrease in the prey population and the unrestricted growth in this population, which could occur assuming ample food, must be modified by the subtraction of this term. Similarly the decrease in the population of the predator must be modified by the addition of the term DPQ since the predator population gains food from its encounters with its prey and therefore more of the predators survive.

The solution of the differential equation depends on the specific values of the constants and will often result in nature in a stable cyclic variation of the populations. This is because as the predators continue to eat the prey, the prey population will fall and become insufficient to support the predator population which itself then falls. However, as the predator population falls, more of the prey survive and consequently the prey population will then increase. This in turn leads to an increase in the predator population since it has more food and the cycle begins again. This cycle maintains the predator and prey populations between certain upper and lower limits. The Volterra differential equations can be solved directly but this solution does not provide a simple relation between the size of the predator and prey populations; therefore, numerical methods of solution should be applied. An interesting description of this problem is given by Simmons (1972).

We now use Matlab to study the behavior of a system of equations of the form (5.28) applied to the interaction of the lynx and its prey, the hare. The choice of the constants $K_1$ , $K_2$ , C, and D is not a simple matter if we wish to obtain a stable situation where the populations of the predator and prey never die out completely but oscillate between upper and lower limits. The Matlab script below uses $K_1=2$ , $K_2=10$ , $C=0.001$ , and $D=0.002$ , and considers the interaction of a population of lynxes and hares where it is assumed that this interaction is the crucial feature in determining the size of the two populations. With an initial population of 5000 hares and 100 lynxes, the script e4s505.m uses these values to produce the graph in Fig. 5.13.

Figure 5.13. Variation in the population of lynxes (dashed line) and hares (solid line) against time, beginning with 5000 hares and 100 lynxes. Accuracy 0.005.

% e4s505.m

% x(1) and x(2) are hare and lynx populations.

simtime = input('enter runtime ');

acc = input('enter accuracy value ');

fv = @(t,x) [2*x(1)-0.001*x(1)*x(2); -10*x(2)+0.002*x(1)*x(2)];

initx = [5000 100]';

options = odeset('RelTol',acc);

[t x] = ode23(fv,[0 simtime],initx,options);

plot(t,x(:,1),'k',t,x(:,2),'k--')

xlabel('Time'), ylabel('Population of hares and lynxes')

For these parameters, there is a remarkably wide variation in the populations of hares and lynxes. The lynx population, although periodically small, still recovers following a recovery of the hare population.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128122563000142

Applications of Systems of Ordinary Differential Equations

Martha L. Abell , James P. Braselton , in Introductory Differential Equations (Fourth Edition), 2014

Biological Systems: Predator-Prey Interaction

Let x(t) and y(t) represent the number of members at time t of the prey and predator populations, respectively. (Examples of such populations include fox-rabbit and shark-seal.) Suppose that the positive constant a is the birth rate of x(t), so that in the absence of the predator dx/dt  = ax and that c   >  0 is the death rate of y, which indicates that dy/dt  =   – cy in the absence of the prey-population. In addition to these factors, the number of interactions between predator and prey affects the number of members in the two populations. Note that an interaction increases the growth of the predator population and decreases the growth of the prey population, because an interaction between the two populations indicates that a predator overtakes a member of the prey population. To include these interactions in the model, we assume that the number of interactions is directly proportional to the product of x(t) and y(t). Therefore, the rate at which x(t) changes with respect to time is

$\mathrm{d}x/\mathrm{d}t=\mathit{ax}-\mathit{bxy},$

where b   >  0. Similarly, the rate at which y(t) changes with respect to time is

$\mathrm{d}y/\mathrm{d}t=-\mathit{cy}+\mathrm{d}\mathit{xy},$

where d   >  0. These equations and initial conditions form the Lotka-Volterra problem

(7.9) $\begin{array}{c}\hfill \mathrm{d}x/\mathrm{d}t=\mathit{ax}-\mathit{bxy}\hfill \\ \hfill \mathrm{d}y/\mathrm{d}t=-\mathit{cy}+\mathit{dx}{y}^{\prime }\hfill \end{array},x\left(0\right)=x_0,y\left(0\right)=y_0.$

Vito Volterra (May 3, 1860, Ancona, Italy-October 11, 1940, Rome, Italy) During World War I, Volterra was a member of the Italian Air Force. After the War, he returned to the University of Rome and his research interests moved from mathematical physics to mathematical biology.

Example 7.3.1

Find and classify the equilibrium points of the Lotka-Volterra system.

Solution

$\mathrm{Solving}\left\{\begin{array}{l}\mathit{ax}-\mathit{bxy}=\mathit{x}\left(\mathit{a}-\mathit{by}\right)=0\\ -\mathit{cy}+\mathrm{d}\mathit{xy}=\mathit{y}\left(-\mathit{c}+\mathrm{d}\mathit{x}\right)=0\end{array}\right.\mathrm{for}\mathit{x}$

and y we have x   =  0 or y  = a/b from the first equation and y   =  0 or x  = c/d from the second equation. Thus, the equilibrium points are E ₀=  (0,0) and E_A  =  (c/d, a/b). The Jacobian matrix of the nonlinear system is $\mathbf{J}\left(x,y\right)=\left(\begin{array}{cc}\hfill a-\mathit{by}\hfill & \hfill -\mathit{bx}\hfill \\ \hfill \mathrm{d}y\hfill & \hfill -c+\mathrm{d}x\hfill \end{array}\right).$ Evaluated at E₀  =   (0,0), the Jacobian is J(E ₀) = $\left(\begin{array}{cc}\hfill a\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -c\hfill \end{array}\right)$ with eigenvalues λ₁  =   –c and λ₂  = a. Because these eigenvalues are real with opposite signs, we classify E ₀  =   (0,0) as a saddle. Similarly, evaluated at E _A   =   (c/d, a/b), the Jacobian is $\mathbf{J}\left(E_A\right)=\left(\begin{array}{cc}\hfill 0\hfill & \hfill -\mathit{bc}/d\hfill \\ \hfill \mathit{ad}/b\hfill & \hfill 0\hfill \end{array}\right),$ with eigenvalues ${\mathrm{\lambda }}_{1,2}=\pm \mathrm{i}\sqrt{\mathit{ac}},$ so the equilibrium point E _A   =  (c/d, a/b) is classified as a center.

In Figure 7.10(a), we show several solutions (which were found using different initial populations) parametrically in the phase plane of this system with a   =  2, b   =  1, c   =  3, and d   =  1. Notice that all of the solutions oscillate about the center. These solutions reveal the relationship between the two populations: prey, x(t), and predator, y(t ). As we follow one cycle counterclockwise beginning on the left extreme of the cycle, we notice that as the prey population, x(t), increases, the predator population, y(t), first slightly decreases (is that really possible?), and then increases until the predator becomes overpopulated. Then, because the prey population is too small to supply the predator population, the predator population decreases, which leads to an increase in the population of the prey. At this point, because the number of predators becomes too small to control the prey population, the number in the prey population becomes overpopulated and the cycle repeats itself, as illustrated in Figure 7.10(b).

How does the period of a solution with an initial point near the equilibrium point (c/d, a/b) compare to that of a solution with an initial point that is not located near this point?

Figure 7.10. (a) Typical solutions of the Lotka-Volterra system-x versus y. (b) A typical solution to the Lotka-Volterra system, x (in dark red; dark gray in print versions) and y as functions of t.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780124172197000077

Some applications of the Reduction theorem and the HKV method

Irina Kareva , Georgy Karev , in Modeling Evolution of Heterogenous Populations, 2020

Example 5.5 Inhomogeneous prey-predator Volterra model with three distributed parameters

The prey-predator Volterra model in its simplest form reads:

$\begin{array}{c}\frac{\mathit{dx}}{\mathit{dt}}=a_{1}x-a_2\mathit{xy}\\ \frac{\mathit{dy}}{\mathit{dt}}=-a_{3}y+a_4\mathit{xy}\end{array},$

where x(t) and y(t) denote prey and predator densities, respectively; a ₁ is the reproduction rate of the prey population; a ₂ is the per capita rate of the consumption of prey by the predator; a ₃ is the death rate of the predator; and a ₄/a ₂ is the fraction of prey biomass, which is converted into predator biomass.

Let us consider the inhomogeneous version of this classical model by assuming that parameters a ₁, a ₂, a ₃ are distributed and the ratio a ₄/a ₂ is fixed (and hence could be chosen equal to 1). Then we get a system

(5.13) $\begin{array}{c}\mathit{dx}/\mathit{dt}=a_{1}x-a_2\mathit{xy}\\ \mathit{dy}/\mathit{dt}=a_2\mathit{xy}-a_{3}y\end{array}.$

Let x(t; a ₁, a ₂), y(t; a ₃) be the densities of the prey and predator populations with respect to parameters a ₁, a ₂, and a ₃ correspondingly, and let

$X\left(t\right)=\underset{A}{\int }x\left(t;a_1,.a_2\right){\mathit{da}}_1{\mathit{da}}_2,Y\left(t\right)=\underset{A}{\int }y\left(t;a_3\right){\mathit{da}}_3$

be the total sizes of the two populations. The initial population sizes X(0), Y(0) and initial distributions P(0; a ₁, a ₂), P ₃(0; a ₃) are assumed to be given. We also assume that the reproduction and death rates are specific for each subpopulation, while consumption is driven by the interaction of the prey (predator) subpopulation with the entire predator (prey) population. Then, assuming the "proportional distribution" of prey among the predators, we can write the inhomogeneous version of Volterra' model in the form:

(5.14) $\begin{array}{l}\frac{\mathit{dx}\left(t;a_1,a_2\right)}{\mathit{dt}}=x\left(t;a_1,a_2\right)\left(a_1-a_2\text{Y}\left(t\right)\right)\\ \frac{\mathit{dy}\left(t;a_3\right)}{\mathit{dt}}=y\left(t;a_3\right)\left(G\left(t\right)-a_3\right)\end{array},$

where $G\left(t\right)=\underset{A}{\int }{a}_{2}x\left(t;a_1,a_2\right){da}_1{da}_2=X\left(t\right)E^t\left[a_2\right]$ .

Theorem 4.3 provides a method for studying this model. The key step is reduction of this model to an escort system of ODE. It is instructive to derive the escort system and the main results informally to clarify the main idea of the method in application to community models.

It is reasonable to assume that parameter a ₃ is stochastically independent of parameters a ₁, a ₂. Let M[λ ₁, λ ₂] be the mgf of the initial joint distribution of parameters a ₁, a ₂, and let M ³[λ ₃] be the mgf of the initial distribution of parameter a ₃.

Let us now formally introduce auxiliary keystone variables q ₁(t), q ₂(t) as solutions to the following Cauchy problem:

(5.15) $\begin{array}{l}\frac{{\mathit{dq}}_1}{\mathit{dt}}=Y\left(t\right)\\ \frac{{\mathit{dq}}_2}{\mathit{dt}}=G\left(t\right)=X\left(t\right)E^t\left[a_2\right],\\ q_1\left(0\right)=q_2\left(0\right)=0\end{array}$

Then system (5.14) can be written as

(5.16) $\begin{array}{l}\frac{\mathit{dx}\left(t;a_1,a_2\right)}{\mathit{dt}}=x\left(t;a_1,a_2\right)\left(a_1-a_2\frac{{\mathit{dq}}_1}{\mathit{dt}}\right)\\ \frac{\mathit{dy}\left(t;a_3\right)}{\mathit{dt}}=y\left(t;a_3\right)\left(\frac{{\mathit{dq}}_2}{\mathit{dt}}-a_3\right)\end{array}.$

Its solution is given by

(5.17) $\begin{array}{l}x\left(t;a_1,a_2\right)=x\left(0;a_1,a_2\right){\mathrm{e}}^{a_{1}t-a_2{q}_1\left(t\right)}\\ y\left(t;a_3\right)=y\left(0;a_3\right){\mathrm{e}}^{q_2\left(t\right)-a_{3}t}\end{array}.$

Now, we can express all values of interest with the help of the mgf-s of the initial distributions and the auxiliary keystone variables:

(5.18) $\begin{array}{l}X\left(t\right)=X\left(0\right)\underset{A}{\int }{\mathrm{e}}^{a_{1}t-a_2{q}_1\left(t\right)}{P}_1\left(0;a_1,a_2\right){da}_1{da}_2=X\left(0\right)M\left[t,,,-q_1\left(t\right)\right]\\ Y\left(t\right)=Y\left(0\right)\underset{A}{\int }{\mathrm{e}}^{q_2\left(t\right)-a_{3}t}{P}_3\left(0;a_3\right){da}_3=Y\left(0\right){\mathrm{e}}^{q_2\left(t\right)}{M}^3\left[-t\right]\end{array}\text{,}$

(5.19) $\begin{array}{l}P\left(t;a_1,a_2\right)=\frac{x\left(t;a_1,a_2\right)}{X\left(t\right)}=\frac{{\mathrm{e}}^{a_{1}t-a_2{q}_1\left(t\right)}}{M\left[t,-q_1\left(t\right)\right]}P\left(0;a_1,a_2\right)\\ P_3\left(t;a_3\right)=\frac{y\left(t;a_3\right)}{Y\left(t\right)}=\frac{{\mathrm{e}}^{-a_{3}t}}{M^3\left[-t\right]}{P}_3\left(0;a_3\right)\\ E^t\left[a_2\right]=\underset{A}{\int }\frac{a_2{\mathrm{e}}^{a_{1}t-a_2{q}_1\left(t\right)}}{M\left[t,-q_1\left(t\right)\right]}P\left(0;a_1,a_2\right){da}_1{da}_2=-\frac{\partial }{\partial q_1}lnM\left[t,,,-q_1\left(t\right)\right]\end{array}\text{.}$

Hence,

(5.20) $X\left(t\right)E^t\left[a_2\right]=-X\left(0\right)\frac{\partial }{\partial q_1}M\left[t,-q_1\right]\text{.}$

Finally, we obtain a closed escort system of nonautonomous equations:

(5.21) $\begin{array}{l}\frac{{dq}_1}{dt}=Y\left(0\right){\mathrm{e}}^{q_2\left(t\right)}{M}^3\left[-t\right],q_1\left(0\right)=0,\\ \frac{{dq}_2}{dt}=-X\left(0\right)\frac{\partial }{\partial q_1}M\left[t,,-q_1\right],q_2\left(0\right)=0\text{.}\end{array}$

Now that we have a solution to the Cauchy problem (5.21), we can obtain explicit formulas for population sizes defined in Eq. (5.18) and current distribution of the parameters defined in Eq. (5.19), which completely solve the problem. In particular, the current mean values of the parameters are given by

(5.22) $\begin{array}{l}{E}^t\left[a_1\right]=\frac{\partial }{\partial t}ln\left(M\left[t,-q_1\left(t\right)\right]\right),\\ E^t\left[a_2\right]=-\frac{\partial }{\partial q_1}ln\left(M\left[t,-q_1\left(t\right)\right]\right),\\ E^t\left[a_3\right]=-\frac{\partial }{\partial t}ln\left(M^3\left[-t\right]\right).\end{array}$

Integrating the equations of system (5.14) over the three parameters, we obtain the following system:

(5.23) $\begin{array}{l}\frac{\mathit{dX}}{\mathit{dt}}=X\left(E^t\left[a_1\right]-E^t\left[a_2\right]Y\right),\\ \frac{\mathit{dY}}{\mathit{dt}}=Y\left(E^t\left[a_2\right]X-E^t\left[a_3\right]\right).\end{array}$

These equations for total sizes of inhomogeneous populations have the same form as the initial Volterra system (5.13). The key difference is that now parameter values are not constants but vary over time according to Eq. (5.22). The phase-parameter portrait of the parametrically homogeneous Volterra' model is well known (see, e.g., Bazykin, 1998). The dynamics of system (5.23) are determined by the parametric point (E ^t [a ₁], E ^t [a ₂], E ^t [a ₃]), which moves across the phase-parameter portrait of model (5.13) over time as the system evolves. This phenomenon, which amounts to "traveling across the phase-parameter portrait of the homogeneous model," is a common feature of corresponding inhomogeneous models. It was well observed on the example of discrete-time models (see Chapter 15). A Volterra-type model of two inhomogeneous populations with logistic reproduction rates and ratio-dependent predator functional response will be studied in Chapter 12, together with additional examples of travel across the phase-parameter portrait.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128143681000059

Open Research Problems: Systems Dynamics, Complex Systems

Bernard P. Zeigler , ... Ernesto Kofman , in Theory of Modeling and Simulation (Third Edition), 2019

24.3.4 Example: Prey-Predator Model

A classical domain for SD models is population dynamics. A reference ecological model in this domain is the Predator-Prey type of system, with its simplest form being the Lotka-Volterra model (Lotka, 1956; Volterra, 1928).

A population of predators Pred and a population of preys Prey interact in the following way. Predators hunt for preys and, when encounters happen, a portion 0<k_eff<1 of the energy of preys is absorbed by the predators (defining their birth rate). The system is as follows:

$\begin{array}{c}\hfill d\mathbf{Prey}(t)/dt=\mathbf{-}prey\mathrm{\_}death\mathrm{\_}r⁎\mathbf{Pred}⁎\mathbf{Prey}\mathbf{+}prey\mathrm{\_}birth\mathrm{\_}r⁎\mathbf{Prey}\hfill \\ \hfill d\mathbf{Pred}(t)/dt=\mathbf{-}pred\mathrm{\_}death\mathrm{\_}r⁎\mathbf{Pred}\mathbf{+}prey\mathrm{\_}death\mathrm{\_}r⁎\mathbf{Pred}⁎\mathbf{Prey}⁎k\mathrm{\_}eff\hfill \end{array}$

where:

1.

$prey\mathrm{\_}birth\mathrm{\_}r$ is the natural growth rate of the prey population feeding on natural resources available in abundance

2.

$prey\mathrm{\_}death\mathrm{\_}r$ is the death rate of the prey population due to encounters with predators

3.

$pred\mathrm{\_}death\mathrm{\_}r$ is the natural death rate of the predator population

4.

$k\mathrm{\_}eff$ is the energy efficiency in growing predators from preys

Fig. 24.3 shows the schematic representations of the Lotka-Volterra model both in SD and DEVS. The DEVS model was obtained automatically using a translator from an XMILE source file into a DEVSML target file, applying the morphism relations proposed in the previous section.

Figure 24.3. Model of the Lotka-Volterra system. System Dynamics schematic (top) and DEVS coupled model (bottom).

Below are illustrative examples of excerpts taken from the XMILE and DEVSML specifications.

The XMILE specification for the Prey <stock> component is shown below.

The <eqn>100</eqn> tag sets the initial condition of the state variable to 100 units. The <inflow> and <outflow> tags define the name for the increment and decrement rates, namely BirthPrey and DeathPrey, respectively.

The DEVSML specification for the Prey stock component is the following:

We see the mapping of the initial condition, with value 100. One input and one output port are created, taking as input the derivative of the variable to be integrated and sending as output the quantized version of the integrated variable. This model must observe user-defined accuracy requirements for the integration procedure at each state variable. We could have left this choice to be set globally for all integrators (as is done in all SD tools) and remove these parameters from each QSS atomic model. ¹⁰

Below we focus on elements related to the Decrement Rate for the Prey Stock. An XMILE declaration involves the <flow> component for the <outflow>DeathPrey</outflow> property defined in the stock above. In turn it requires the declaration of a Converter to specify the prey_death_r parameter (set with the <aux> component). XMILE also requires to make explicit the dependencies of variables in the special section <isee:dependencies>.

A possible corresponding DEVSML specification is shown below.

The listing above presents DEVS atomic models of type F_decr and ConstGen. We also make explicit the Internal Couplings to wire the input ports of the F_decr model in the context of the overall DEVS coupled model. The full specification will of course require also the QSS model type for the Predator population, the F_der_Prey model to combine the decrement rate with the increment rate, and so on.

The atomic model ConstGen for the constant value generation could be avoided by hardcoding that value directly into the "function" of the F_decr model. Yet, the atomic model option enables desirable flexibilities such as easy replacement by more complex generators (e.g. a slowly varying modulated parameter, yielding a new variable) or easy rewiring to other submodels of the environment within which the population is evolving (e.g. a death rate could be a function of average temperature, in turn related to climate change effects)

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128133705000353

Applications of Systems of Ordinary Differential Equations

Martha L. Abell , James P. Braselton , in Introductory Differential Equations (Fifth Edition), 2018

7.3 Nonlinear Systems of Equations

•

Biological Systems: Predator-Prey Interaction

•

Physical Systems: Variable Damping

Several special equations and systems that arise in the study of many areas of applied mathematics can be solved using the techniques of Chapter 6 . These include the predator-prey population dynamics problem, the Van-der-Pol equation that models variable damping in a spring-mass system, and the Bonhoeffer-Van-der-Pol (BVP) oscillator.

Biological Systems: Predator-Prey Interaction

Let $x(t)$ and $y(t)$ represent the number of members at time t of the prey and predator populations, respectively. (Examples of such populations include fox-rabbit and shark-seal.) Suppose that the positive constant a is the birth rate of $x(t)$ , so that in the absence of the predator $dx/dt=ax$ and that $c>0$ is the death rate of y, which indicates that $dy/dt=-cy$ in the absence of the prey-population. In addition to these factors, the number of interactions between predator and prey affects the number of members in the two populations. Note that an interaction increases the growth of the predator population and decreases the growth of the prey population, because an interaction between the two populations indicates that a predator overtakes a member of the prey population. To include these interactions in the model, we assume that the number of interactions is directly proportional to the product of $x(t)$ and $y(t)$ . Therefore, the rate at which $x(t)$ changes with respect to time is

$dx/dt=ax-bxy,$

where $b>0$ . Similarly, the rate at which $y(t)$ changes with respect to time is

$dy/dt=-cy+dxy,$

where $d>0$ . These equations and initial conditions form the Lotka-Volterra problem

(7.9) $\begin{array}{cc}\hfill dx/dt& =ax-bxy\hfill \\ \hfill dy/dt& =-cy+dxy\hfill \end{array},x(0)=x_0,y(0)=y_0.$

Alfred James Lotka (1880–1949) American biophysicist, born in Ukraine; wrote first text on mathematical biology.

Example 7.8

Find and classify the equilibrium points of the Lotka-Volterra system.

Solution: Solving $\{\begin{array}{c}ax-bxy=x(a-by)=0\hfill \\ -cy+dxy=y(-c+dx)=0\hfill \end{array}$ for x and y we have $x=0$ or $y=a/b$ from the first equation and $y=0$ or $y=c/d$ from the second equation. Thus, the equilibrium points are $E_0=(0,0)$ and $E_A=(c/d,a/b)$ . The Jacobian matrix of the nonlinear system is $\mathbf{J}(x,y)=\left(\begin{array}{cc}\hfill a-by\hfill & \hfill -bx\hfill \\ \hfill dy\hfill & \hfill -c+dx\hfill \end{array}\right)$ . Evaluated at $E_0=(0,0)$ , the Jacobian is $\mathbf{J}(E_0)=\left(\begin{array}{cc}\hfill a\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -c\hfill \end{array}\right)$ with eigenvalues ${\lambda }_1=-c$ and ${\lambda }_2=a$ . Because these eigenvalues are real with opposite signs, we classify $E_0=(0,0)$ as a saddle. Similarly, evaluated at $E_A=(c/d,a/b)$ , the Jacobian is $\mathbf{J}(E_A)=\left(\begin{array}{cc}\hfill 0\hfill & \hfill -bc/d\hfill \\ \hfill ad/b\hfill & \hfill 0\hfill \end{array}\right)$ , with eigenvalues ${\lambda }_{1,2}=\pm i\sqrt{ac}$ , so the equilibrium point $E_A=(c/d,a/b)$ is classified as a center.

Vito Volterra (May 3, 1860, Ancona, Italy-October 11, 1940, Rome, Italy) During World War I, Volterra was a member of the Italian Air Force. After the War, he returned to the University of Rome and his research interests moved from mathematical physics to mathematical biology.

In Fig. 7.10A, we show several solutions (which were found using different initial populations) parametrically in the phase plane of this system with $a=2$ , $b=1$ , $c=3$ , and $d=1$ . Notice that all of the solutions oscillate about the center. These solutions reveal the relationship between the two populations: prey, $x(t)$ , and predator, $y(t)$ . As we follow one cycle counterclockwise beginning, we notice that as the prey population, $x(t)$ , increases, the predator population, $y(t)$ , first slightly decreases (is that really possible?), and then increases until the predator becomes overpopulated. Then, because the prey population is too small to supply the predator population, the predator population decreases, which leads to an increase in the population of the prey. At this point, because the number of predators becomes too small to control the prey population, the number in the prey population becomes overpopulated and the cycle repeats itself, as illustrated in Fig. 7.10B.

Figure 7.10. (A) Typical solutions of the Lotka-Volterra system-x versus y. (B) A typical solution to the Lotka-Volterra system, x (in dark red) and y as functions of t.

Interestingly, studies show that in addition to death, predation causes stress in the prey population due to the high rate of encounters, which causes reductions in reproduction rates and offspring viability as well as susceptibility to diseases, all of which are detrimental to the population. ¹ □

How does the period of a solution with an initial point near the equilibrium point $(c/d,a/b)$ compare to that of a solution with an initial point that is not located near this point?

Physical Systems: Variable Damping

In some physical systems, energy is fed into the system when there are small oscillations, and energy is taken from the system when there are large oscillations. This indicates that the system undergoes "negative damping" for small oscillations and "positive damping" for large oscillations. A differential equation that models this situation is Van-der-Pol's equation

(7.10) $x^{″}+\mu (x^2-1)x^{\prime }+x=0,$

where μ is a positive constant. We can transform this second order differential equation into a system of first order differential equations with the substitution $x^{\prime }=y$ . Hence $y^{\prime }=x^{″}=-x-\mu (x^2-1)x^{\prime }=-x-\mu (x^2-1)y$ , so the corresponding system of equations is

(7.11) $\begin{array}{cc}\hfill x^{\prime }& =y\hfill \\ \hfill y^{\prime }& =-x-\mu (x^2-1)y\hfill \end{array},$

which is solved using an initial position $x(0)=x_0$ and an initial velocity $y(0)=x^{\prime }(0)=y_0$ . Notice that $\mu (x^2-1)$ represents the damping coefficient. This system models variable damping because $\mu (x^2-1)<0$ when $-1<x<1$ and $\mu (x^2-1)>0$ when $|x|>1$ . Therefore, damping is negative for the small oscillations, $-1<x<1$ , and positive for the large oscillations, $|x|>1$ .

Balthasar Van-der-Pol (1889–1959) Dutch applied mathematician and engineer.

Example 7.9

Find and classify the equilibrium points of the system of differential equations that is equivalent to Van-der-Pol's equation.

Solution: We find the equilibrium points by solving

$\begin{array}{cc}\hfill y& =0\hfill \\ \hfill -x-\mu (x^2-1)y& =0\hfill \end{array}.$

From the first equation, we see that $y=0$ . Substitution of $y=0$ into the second equation yields $x=0$ as well. Therefore, the (only) equilibrium point is $(0,0)$ .

The Jacobian matrix for this system is

$\mathbf{J}(0,0)=\left(\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1-2\mu xy\hfill & \hfill -\mu (x^2-1)\hfill \end{array}\right).$

At $(0,0)$ , we have the matrix $\mathbf{J}(0,0)=\left(\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill \mu \hfill \end{array}\right)$ . We find the eigenvalues of $\mathbf{J}(0,0)$ by solving

$\left|\begin{array}{cc}\hfill -\lambda \hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill -\mu -\lambda \hfill \end{array}\right|={\lambda }^2-\mu \lambda +1=0,$

Graph several solution curves in the phase plane for the Van-der-Pol equation if $\mathit{\mu }=1/10$ and $\mathit{\mu }=1/1000$ . Compare these graphs to the corresponding solution curves for the equation $x^{″}+x=0$ . Are they similar? Why?

which has roots ${\lambda }_{1,2}=\frac{1}{2}(\mu \pm \sqrt{{\mu }^2-4})$ . Notice that if $\mu >2$ , then both eigenvalues are positive and real, so we classify $(0,0)$ as an unstable node. On the other hand, if $0<\lambda <2$ , the eigenvalues are a complex conjugate pair with a positive real part. Hence $(0,0)$ is an unstable spiral point. (We omit the case when $\mu =2$ because the eigenvalues are repeated.) In Fig. 7.11, we show several curves in the phase plane that begin at different points for various values of μ. In each figure, we see that all of the curves approach a curve called a limit cycle. Physically, the fact that the system has a limit cycle indicates that for all oscillations the motion eventually becomes periodic, which is represented by a closed curve in the phase plane.

Figure 7.11. From left to right, (A) μ = 1/2, (B) μ = 1, (C) μ = 3/2, (D) μ = 3.

On the other hand, in Fig. 7.12 we graph the solution that satisfies the initial conditions $x(0)=1$ and $y(0)=0$ parametrically and individually for various values of μ. Notice that for small values of μ the system more closely approximates that of the harmonic oscillator because the damping coefficient is small. The curves are more circular than for those for larger values of μ. □

Figure 7.12. From left to right, (A) μ = 1/4, (B) μ = 1/2, (C) μ = 1, (D) μ = 3/2, (E) μ = 2, (F) μ = 3.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128149485000070

PD Bifurcation and Chaos Behavior in a Predator-Prey Model With Allee Effect and Seasonal Perturbation

Afef Ben Saad , ... Zeraoulia Elhadj , in Recent Advances in Chaotic Systems and Synchronization, 2019

1 Introduction

Recently, dynamical relationship between species has been intensively studied and will continue to be one of the most important themes of ecology [1–5]. Such systems are generally depicted by nonlinear polynomial models which are based on nonlinear Lotka-Volterra model [6] and have the following basic form [7].

(1) $\begin{array}{l}\hfill \left\{\begin{array}{l}\frac{d{x}_1}{dt}=a{x}_1-f(x_1)x_2\\ \frac{d{x}_2}{dt}=c{x}_1{x}_2-cm{x}_2\\ \end{array}\right.\end{array}$

where x ₁ and x ₂ are the size of the prey population and the predator population, respectively. a is the prey's growth rate in absence of the predator. f(x ₁) is the functional response of the predator to prey density. c represents the predator's conversion efficiency, and m is the mortality rate of the predator depending on the predator's efficiency.

On the other hand, during the last few years, bifurcation analysis of predator-prey models which incorporate one or more extra factors has been considered as a challenging research topic. The most famous extra factors are time delay effect [8, 9], impulsive effect [10], seasonal effect [11], and Allee effects [12, 13].

In [14], a bifurcation analysis of the predator-prey Bazykin and Berezovskaya (BB) model affected by Allee effect considered in [15], and introduced noninteger polynomials, is carried out. To extend the set of information related to this model, the authors propose modeling the prey's growth rate and the predator's growth rate using polynomials with noninteger elements such as:

(2) $\begin{array}{l}\hfill a=x_1^{q_1-1}(x_1-l)(1-x_1),f(x_1)=x_1^{q_2}\end{array}$

where l is the Allee effect threshold , q ₁ and q ₂ are two noninteger elements. Using Eq. (2), system (1) is then given by

(3) $\begin{array}{l}\hfill \left\{\begin{array}{l}\frac{d{x}_1}{dt}=x_1^{q_1}(x_1-l)(1-x_1)-x_1^{q_2}{x}_2\\ \frac{d{x}_2}{dt}=c(x_1-m)x_2\\ \end{array}\right.\end{array}$

For the system described by Eq. (3), several bifurcation points are detected according to the Allee effect's type and the noninteger elements but with standard bifurcation types such as Hopf bifurcation (H) and branch point bifurcation (BP).

In this chapter, a modified structure of prey's growth rate and predator's efficiency based on noninteger elements is proposed. Detailed analyses of the nonperturbed and the perturbed system are investigated by considering the two cases study of Strong and Weak Allee effects.

This chapter is arranged as follows: Section 2 presents the predator-prey BB model with noninteger elements in the prey's growth rate and predator's efficiency structure. Preliminaries are given in Section 3. Stability and Hopf bifurcation analysis are investigated in Section 4. In Section 5, numerical simulations are depicted. Finally, the research work undertaken is discussed in Section 6.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B978012815838800011X

Systems of nonlinear differential equations

Henry J. Ricardo , in A Modern Introduction to Differential Equations (Third Edition), 2021

7.1 Equilibria of nonlinear systems

Recall that an equilibrium point of a differential equation or a system of differential equations is a constant solution. If we look at the two (somewhat similar) equations

$(1)y^{\prime }=-y(2)y^{\prime }=-y(1-y),$

we see some important differences between linear and nonlinear equations.

Eq. (1) is linear, but more fundamentally it is separable, so it is easy to find the general solution, $y=C{e}^{-t}$ , where C is an arbitrary constant. (We recognize that $C=y(0)$ , the initial state of the system being modeled by the equation.)

Now Eq. (2) is nonlinear and separable, and its general solution is $y=\frac{1}{1+C{e}^t}$ , where $C=1/y(0)-1$ . (Verify the solutions to both equations.)

Let's examine some typical solution curves for Eq. (1). Fig. 7.1 shows that there is only one equilibrium solution, $y\equiv 0$ , and this is a sink. (Review Section 2.6 if necessary.) If an object described by the equation starts off at zero (that is, if $C=0$ ), it remains at zero for all time. If the object's initial state is not zero, then the object will approach the solution $y\equiv 0$ as its asymptotically stable solution (or sink).

Figure 7.1. Solutions of y′ = −y; y(0)=5,3,1,−2,−4

On the other hand, Fig. 7.2 shows the same kind of information for Eq. (2). For such a nonlinear equation there can be more than one equilibrium solution, in this case $y\equiv 0$ and $y\equiv 1$ . Also note that some solutions of a nonlinear equation may become unbounded as t approaches some finite value. (If $y(0)>1$ , for what value of t does the denominator of the general solution to Eq. (2) vanish?)

Figure 7.2. Solutions of y′ = −y(1 −y)

In contrast, all solutions of a linear equation or a system of linear equations are defined for all values of the independent variable. Finally, looking closely at the behavior of solutions of Eq. (2) with different initial values, we see that the solutions starting off above 1 behave differently from those solutions with initial values less than 1. The equilibrium solution $y\equiv 0$ is a sink if $y(0)<1$ and $y\equiv 1$ is a source if $y(0)>1$ . Furthermore, for solutions with initial values C greater than 1, the line $t=\ln (-\frac{1}{C})$ is a vertical asymptote. We say that the solution "blows up in finite time." (Note that when $y(0)>1$ , $C=1/y(0)-1<0$ , so that the blow-up point t is the logarithm of a positive quantity.) The last three types of behavior cannot occur when we are dealing with a linear equation. You should expect that the situation with nonlinear systems is appropriately complicated.

Let's look at an example of a nonlinear system and its behavior near its equilibrium points.

Example 7.1.1 Stability of a Nonlinear System

The nonlinear system

$$\begin{gathered} x^{\prime }=x-x^2-xy \\ {y}^{\prime }=-y-y^2+2xy \end{gathered}$$

represents two populations interacting in a predator-prey relationship. This is essentially a Lotka–Volterra system (see Section 7.4) with "crowding" terms (the squared terms) added for both species.

To calculate the equilibrium points of this system, we solve the system $\{x^{\prime }=0,y^{\prime }=0\}$ , which is the same as the nonlinear algebraic system

$$\begin{gathered} \text{(A)}x(1-x-y)=0 \\ \text{(B)}y(-1-y+2x)=0. \end{gathered}$$

Clearly, the origin, $x=y=0$ , is an equilibrium point. Logically, there are only three other cases to examine: (1) $x=0$ , $y\ne 0$ ; (2) $x\ne 0$ , $y=0$ ; and (3) $x\ne 0$ , $y\ne 0$ . Assuming case 1, we can eliminate Eq. (A) and examine (B), which becomes $y(-1-y)=0$ . Because $y\ne 0$ , we conclude that $-1-y=0$ , or $y=-1$ . Thus, our second equilibrium point is $(0,-1)$ . Moving to case 2, we can ignore Eq. (B) and focus on (A), which now looks like $x(1-x)=0$ . Because we are assuming in case 2 that $x\ne 0$ , we can see that $x=1$ , which gives us the third equilibrium point (1, 0). Finally, if $x\ne 0$ and $y\ne 0$ , our system of algebraic equations becomes

$$\begin{gathered} \text{(A2)}x+y=1 \\ \text{(B2)}y-2x=-1. \end{gathered}$$

(We have divided out x and y in (A) and (B) and then rearranged the terms of each equation.) Subtracting (B2) from (A2) gives us $3x=2$ , or $x=\frac{2}{3}$ . Substituting this value of x in (A2) yields $y=\frac{1}{3}$ . Therefore, the last equilibrium point is $(\frac{2}{3},\frac{1}{3})$ .

In terms of a population problem, the only interesting equilibrium point is the last one we found. (Why is this so?) If we look at a slope field for the original system of nonlinear differential equations near the point $(\frac{2}{3},\frac{1}{3})$ , we see an interesting pattern (Fig. 7.3a).

Figure 7.3a. Slope field for x′ =x −x ² −xy, y′ = −y −y ² + 2xy near $(\frac{2}{3},\frac{1}{3})$

The apparent spiraling of solutions into the equilibrium point can be seen more clearly if we show some numerically generated solution curves (Fig. 7.3b). Fig. 7.3b represents a predator-prey population that is stabilizing. If the units are thousands of creatures, then the X population is heading for a steady population of about 667, whereas the Y population has 333 as its stable value.

Figure 7.3b. Phase portrait for x′ =x −x ² −xy, y′ = −y −y ² + 2xy near $(\frac{2}{3},\frac{1}{3})$ ; (x(0),y(0))=(0.2,1), (0.8,0.8), (0.8,0.5), (1,0.7), (1,0.2), (0.5,1)

Mathematically, however, we should look at the entire phase portrait to understand the complex behavior of nonlinear systems. We'll return for a detailed analysis of this system in Example 7.3.1.

Exercises 7.1

A

Find all the equilibrium points for each of the systems in Problems 1–14, using technology if necessary.

1.

$x^{\prime }=-x+xy$ , $y^{\prime }=-y+2xy$

2.

$x^{\prime }=x-xy$ , $y^{\prime }=y-xy$

3.

$x^{\prime }=x^2-y^2$ , $y^{\prime }=x-xy$

4.

$x^{\prime }=1-y^2$ , $y^{\prime }=1-x^2$

5.

$x^{\prime }=x+y+2xy$ , $y^{\prime }=-2x+y+y^3$

6.

$x^{\prime }=y(1-x^2)$ , $y^{\prime }=-x(1-y^2)$

7.

$x^{\prime }=x-x^2-xy$ , $y^{\prime }=3y-xy-2y^2$

8.

$x^{\prime }=1-y$ , $y^{\prime }=x^2-y^2$

9.

$x^{\prime }=(1+x)\sin y$ , $y^{\prime }=1-x-\cos y$

10.

$x^{\prime }=3y-e^x$ , $y^{\prime }=2x-y$ [Hint: There are two equilibrium points. Use your CAS to approximate these points.]

11.

$x^{\prime }=y^2-x^2$ , $y^{\prime }=x-1$

12.

$x^{\prime }=x^2-y$ , $y^{\prime }=y^2-x$

13.

$x^{\prime }=xy(1-x)$ , $y^{\prime }=y(1-\frac{y}{x})$

14.

$x^{\prime }=y$ , $y^{\prime }=-\sin x-3y$

B

15.

Use technology to find all the equilibrium points of the system

$x^{\prime }=-y,y^{\prime }=(x^4+4x^3-x^2-4x+y)/8.$

16.

A two-mode laser produces two different kinds of photons, whose numbers are $n_1$ and $n_2$ . The equations governing the rates of photon production are

$$\begin{gathered} {\dot{n}}_1=G_{1}N{n}_1-k_1{n}_1 \\ {\dot{n}}_2=G_{2}N{n}_2-k_2{n}_2, \end{gathered}$$

where $N(t)=N_0-a_1{n}_1-a_2{n}_2$ is the number of excited atoms. The parameters $G_1,G_2,k_1,k_2,a_1,a_2$ , and $N_0$ are all positive. Use a CAS "solve" command to find all the equilibrium points of the system.

17.

A chemostat is a device for growing and studying bacteria by supplying nutrients and maintaining convenient levels of the bacteria in a culture. One model of a chemostat is the nonlinear system

$$\begin{gathered} \frac{dN}{dt}=a_1\left(\frac{C}{1+C}\right)N-N \\ \frac{dC}{dt}=-\left(\frac{C}{1+C}\right)N-C+a_2, \end{gathered}$$

where $N=N(t)$ denotes the bacterial density at time t; $C=C(t)$ denotes the concentration of nutrient; and $a_1,a_2$ are positive parameters. Use technology to find all the equilibrium solutions $(N^{⁎},C^{⁎})$ of the system.

18.

In the absence of damping and any external force, the motion of a pendulum is described by the equation $\frac{d^2\theta }{d{t}^2}+\frac{g}{L}\sin \theta =0$ , where θ is the angle between the pendulum and the downward vertical, g is the acceleration due to gravity, and L is the length of the pendulum.

a.

Write this equation as a system of two first-order equations.

b.

Describe all the equilibrium points of the system.

C

19.

Use technology to find all the equilibrium solutions of the system

$x^{\prime }=x-y^2+a+bxy,y^{\prime }=0.2y-x+x^3,$

where $a=1.28$ and $b=1.4$ . (Round to the nearest thousandth.)

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128182178000142

Existence of periodic travelling waves solutions in predator prey model with diffusion

Radouane Yafia , M.A. Aziz-Alaoui , in Applied Mathematical Modelling, 2013

5 Numerical simulations and discussions

With Matlab software we illustrate our result by some numerical simulations. The method used to compute the travelling waves solutions is to use BVP solver with projection conditions and a phase condition with the following parameters values $a=0.5\text{,}b=0.25\text{,}c=2\text{,}{e}_1=2\text{,}{e}_2=2.5$ and $0<d<1$ From Fig. 1, we observe that for small density prey the density of predator decreases to reach the equilibrium state $u_{\ast }$ and the prey density of prey increases to reach the equilibrium $v_{\ast }$ for $d<d_0$ . But in Fig. 2, the two population coexist but they oscillate periodically around the nontrivial steady state $(u_{\ast }\text{,}{v}_{\ast })$ for $d=d_0$ and the two population can survive together for long time (because the prey may exists with a high density $({10}^{12})$ and there is no extinction of the two population). In Fig. 3, an illustration of the periodic travelling waves solution with respect to the time and space variables and Neumann boundary conditions. The solutions of the reaction diffusion system are represented by a surface, for each fixed x the solutions are represented by lines which correspond with travelling waves in dimension two. We observe that, the periodicity of the invasion of the predators imply the periodicity of the prey.

Fig. 1. Travelling wave solution connecting the steady states $(1\text{,}0)$ and $(u_{\ast }\text{,}{v}_{\ast })$ , predator (blue line) and prey (green line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. An illustration of the existence of periodic travelling waves connecting the steady states $(1\text{,}0)$ and $(u_{\ast }\text{,}{v}_{\ast })$ for $d=d_0$ predator (blue line) and prey (green line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Periodic travelling waves with respect to the time and space variables.

In this paper, we show the existence of periodic travelling waves solutions via Hopf bifurcation Theorem by considering the diffusion coefficient as a parameter of bifurcation. A numerical simulations are given to illustrate the theoretical study of the reaction diffusion system modelling the invasion of the prey by the predator species. The predation is an established case of cycling in prey species. Here, the ability of predation to explain periodic travelling waves in prey population (see Fig. 2), which have recently been found in a number of spatiotemporal field studies. The nature of periodic waves in this systems, and the way in which they can be generated by the invasion of predators into a prey population. A theoretical calculation that predicts, as a function of one parameter ratio, whether such an invasion will leads to a periodic travelling waves that would be observed. The result gives an insight into the types of predator prey systems in which one would expects to see periodic travelling waves following an invasion by predators.

5.1 Conclusion

In this work, a spatio-temporal system reaction–diffusion modelling predator–prey population with modified Leslie-Gower and Holling type-II functional response is studied. By using the semigroup theory, we showed the well-posedeness of the problem and the positivity of solutions. In Theorem 3.1, the conditions of change of stability of possible steady states are given. By considering the diffusion coefficient as a parameter of bifurcation and by applying the Hopf bifurcation Theorem, we prove that there exists a critical value of diffusion parameter at which a small amplitude periodic solution in ${\mathbb{R}}^4$ which corresponds to a small amplitude travelling wave solution connecting the two equilibrium points occurs (Theorem 4.1).

Read full article

URL:

https://www.sciencedirect.com/science/article/pii/S0307904X12004593

Controlling predator–prey discrete dynamics utilizing a threshold policy with hysteresis

Magno Enrique Mendoza Meza , Amit Bhaya , in Applied Mathematics and Computation, 2011

4 Rosenzweig–MacArthur model

The Rosenzweig–MacArthur model can be regarded as one of the simplest nontrivial paradigms that can be proposed by refining the more classical but biologically unrealistic Lotka–Volterra model. The Rosenzweig–MacArthur model is as follows:

(16) $\left\{\begin{array}{c}\dot{x}(t)=\mathit{rx}\left(1-\frac{x}{K}\right)-\frac{\mathit{xy}}{x+A}\text{,}\hfill \\ \dot{y}(t)=\frac{\mathit{sA}(x-J)}{(J+A)(x+A)}y\text{,}\hfill \end{array}\right.$

where r is the intrinsic growth rate of the prey, K is the carrying capacity of the environment, A is the half saturation constant, s conversion efficiency of predator, and J is the minimum prey population for which the predator can survive below the carrying capacity K; J  < K. The state variable x denotes the prey population density and the state variable y denotes the predator population density. The discretization of the Rosenzweig–MacArthur model (16) is carried out according to [2]. The following replacements are made in the first equation of (16): (i) $\frac{\mathit{dx}}{\mathit{dt}}\to \frac{x_{k+1}-x_k}{\varphi }$ , (ii) x  → x _k , (iii) x ²  → x _k+1 x _k , (iii) xy  → x _k+1 y _k ; and in the second equation of (16): (i) $\frac{\mathit{dy}}{\mathit{dt}}\to \frac{y_{k+1}-y_k}{\varphi }$ , (ii) y  → y _k , (iii) xy  → x _k y _k+1. The resulting Rosenzweig–MacArthur discrete model is as follows:

(17) $\left\{\begin{array}{c}{x}_{k+1}=\frac{(1+r{\varphi }_1(h))x_k}{1+\frac{r}{K}{\varphi }_1(h)x_k-\frac{{\varphi }_1(h)y_k}{x_k+A}}\text{,}\hfill \\ y_{k+1}=\frac{\left(1-\frac{\mathit{sAJ}{\varphi }_2(h)}{(J+A)x_k}\right)y_k}{1-\frac{\mathit{sA}{\varphi }_2(h)x_k}{(J+A)(x_k+A)}}\text{,}\hfill \end{array}\right.$

where ϕ ₁(h)   =   exp(h)   −   1 and ϕ ₂(h)   =   exp(h)   −   1. The Rosenzweig–MacArthur model subject to a TPH is as follows:

(18) $\left\{\begin{array}{c}{x}_{k+1}=\frac{(1+r{\varphi }_1(h))x_k}{1+\frac{r}{K}{\varphi }_1(h)x_k-\frac{{\varphi }_1(h)y_k}{x_k+A}}\text{,}\hfill \\ y_{k+1}=\frac{\left(1-\frac{\mathit{sAJ}{\varphi }_2(h)}{(J+A)x_k}\right)y_k}{1-\frac{\mathit{sA}{\varphi }_2(h)x_k}{(J+A)(x_k+A)}}-{\epsilon }_2{y}_k{\psi }_{\mathit{hys}}(\tau )\text{.}\hfill \end{array}\right.$

The equilibrium point of the system without harvesting (without control) is

$z^{\mathit{VNC}}=(x^{\mathit{VNC}}\text{,}{y}^{\mathit{VNC}})=\left(J\text{,}-\frac{r}{K}(x^{\mathit{VNC}}+A)(x^{\mathit{VNC}}-K)\right)\text{.}$

The equilibrium point of the system with harvesting (proportional control) is

$z^{\mathit{VC}}=(x^{\mathit{VC}}\text{,}{y}^{\mathit{VC}})=\left(\frac{A(\mathit{sJ}\varphi +{\epsilon }_2(J+A))}{-{\epsilon }_2(J+A-\mathit{sA}\varphi )+\mathit{sA}\varphi }\text{,}-\frac{r}{K}(x^{\mathit{VC}}+A)(x^{\mathit{VC}}-K)\right)\text{.}$

If ${\epsilon }_2<\frac{\mathit{sA}\varphi }{J+A-\mathit{sA}\varphi }$ the equilibrium point of the system is positive. We consider in this paper values of ${\epsilon }_2>\frac{\mathit{sA}\varphi }{J+A-\mathit{sA}\varphi }$ . We choose the threshold level as follows:

$\tau =y_k-y_{\mathit{th}}$

such that y ^VC   < y _th   < y ^VNC , i.e., z ^VNC and z ^VC are virtual. Simulation of the model (18) under TPH is shown in Fig. 7.

Read full article

URL:

https://www.sciencedirect.com/science/article/pii/S0096300311003055

Global asymptotic stability of almost periodic solution for a multispecies competition-predator system with time delays

Xiaojie Lin , ... Yansen Lv , in Applied Mathematics and Computation, 2013

1 Introduction

This paper deals with the almost periodic solution of a multispecies competition-predator system with Holling III functional response and time delays

(1.1) $\left\{\begin{array}{c}{\dot{x}}_i(t)=x_i(t)\left[b_i(t)-{\displaystyle \sum _{k=1}^n}{a}_{\mathit{ik}}(t)x_k(t-{\tau }_k(t))-{\displaystyle \sum _{k=1}^m}\frac{c_{\mathit{ik}}(t)x_i(t)y_k(t)}{x_i^2(t)+f_{\mathit{ik}}(t)}\right]\text{,}i=1\text{,}2\text{,}\dots \text{,}n\text{,}\hfill \\ {\dot{y}}_j(t)=y_j(t)\left[-r_j(t)+{\displaystyle \sum _{k=1}^n}\frac{d_{\mathit{kj}}(t)x_k^2(t-{\delta }_k(t))}{x_k^2(t-{\delta }_k(t))+f_{\mathit{kj}}(t)}-{\displaystyle \sum _{k=1}^m}{e}_{\mathit{jk}}(t)y_k(t-{\sigma }_k(t))\right]\text{,}j=1\text{,}2\text{,}\dots \text{,}m\text{,}\hfill \end{array}\right.$

with initial conditions

(1.2) $\begin{array}{c}{x}_i(\theta )={\varphi }_i(\theta )\text{,}\theta \in [-\tau \text{,}0]\text{,}{\varphi }_i(\theta )\in C([-\tau \text{,}0]\text{,}{R}_{+})\text{,}i=1\text{,}2\text{,}\dots \text{,}n\text{,}\hfill \\ y_j(\theta )={\psi }_j(\theta )\text{,}\theta \in [-\tau \text{,}0]\text{,}{\psi }_j(\theta )\in C([-\tau \text{,}0]\text{,}{R}_{+})\text{,}j=1\text{,}2\text{,}\dots \text{,}m\text{,}\hfill \end{array}$

where $x_i(t)\text{,}{y}_j(t)$ denote the size of prey and predator population at time $t\text{;}{b}_i(t)\text{,}{a}_{\mathit{il}}(t)\text{,}{c}_{\mathit{ik}}(t)\text{,}{d}_{\mathit{lj}}(t)\text{,}{f}_{\mathit{ij}}(t)\text{,}{r}_j(t)\text{,}{e}_{\mathit{jk}}(t)\text{,}(i\text{,}l=1\text{,}2\text{,}\dots \text{,}n\text{;}j\text{,}k=1\text{,}2\text{,}\dots \text{,}m)$ are all continuous positive almost periodic functions on R with the ecology meaning as follows: $r_i$ is the prey population grows in the absence of predators, $r_j$ is the predator population decays in the absence of preys, $a_{\mathit{ik}}$ is the prey population decays in the competition among the preys, $e_{\mathit{jk}}$ is the predator population decays in the competition among the predator, $c_{\mathit{ik}}$ is the prey is fed upon by the predators, $d_{\mathit{lj}}$ is the coefficient of transformation from preys to predators, $(i\text{,}l=1\text{,}2\text{,}\dots \text{,}n\text{,}j\text{,}k=1\text{,}2\text{,}\dots \text{,}m)$ . The constant $\tau$ is $\tau ={\max }_{t\in R}\{{\tau }_i(t)\text{,}{\delta }_j(t)\text{,}{\sigma }_k(t)\text{,}i\text{,}j=1\text{,}2\text{,}\dots \text{,}n\text{;}k=1\text{,}2\text{,}\dots \text{,}m\}$ , where ${\tau }_i(t)\text{,}{\delta }_j(t)\text{,}{\sigma }_k(t)$ are nonnegative and continuously differentiable almost periodic functions on R, and ${\min }_{t\in R}\{1-{\dot{\tau }}_i(t)\text{,}1-{\dot{\delta }}_j(t)\text{,}1-{\dot{\sigma }}_k(t)\}>0$ .

The study of predator–prey dynamics was originated in the 1920s in the works of Lotka and Volterra who showed for a one-predator-one-prey model (known as the standard Lotka–Volterra model) that the predator and prey permanently oscillate for any positive initial conditions. The dynamic relationship between predator and prey has long been and will continue to be one of the dominant themes in population dynamics due to its universal existence and importance. The predator–prey models have been extensively studied, see papers [1–13] and therein. For example, Zhou [1] discussed a discrete multispecies cooperation and competition predator–prey systems. For general nonautonomous case, the permanence and the global stability of the system were obtained; for periodic case, the existence results of a globally stable positive periodic solution of the system were established. Li, Zhao and Ye [2] considered a delay nonautonomous Lotka–Volterra type multispecies competitive system with harvesting terms. An existence theorem of $2n$ positive periodic solutions was established by using the coincidence degree theory. Zhang and Teng [3] were concerned with a class of nonautonomous N-species Lotka–Volterra type competitive system with time delays and impulsive perturbations. New criteria on coexistence and global attractivity for all species other than some inclining to extinction were established. Shi, Li and Chen [4] proposed a Lotka–Volterra competitive system with infinite delay and feedback controls. By using the method of multiple Lyapunov functionals and by developing a new analysis technique, the sufficient conditions which guarantee that some of the n species are driven to extinction were established. Liu and Wang [5] studied a general stochastic Lotka–Volterra system with infinite delays. Some sufficient criteria for global asymptotic stability of were obtained. Some simulation figures were introduced to support the analytical findings. Braza [6] discussed a predator–prey model in which a modified Lotka–Volterra interaction term is used as the functional response of the predator to the prey. The interaction term is proportional to the square root of the prey population, which appropriately models systems in which the prey exhibits strong herd structure implying that the predator generally interacts with the prey along the outer corridor of the herd. Hilhorst, Martin and Mimura [7] considered a competition-diffusion system for two competing species, under the situation where the two species spatially segregate as the interspecific competition rate becomes large, the authors showed that the resulting limit problem turns out to be a free boundary problem. Sun, Zhang and Song et al. [8] studied the spatiotemporal complexity in a predator–prey system with Holling–Tanner form. A theoretical analysis of emerging spatial pattern was presented and wavelength and pattern speed were calculated.

Cai, Huang and Chen [13] proposed the systems (1.1) and (1.2), they assume that the coefficients of the system (1.1)

$b_i(t)\text{,}{a}_{\mathit{il}}(t)\text{,}{c}_{\mathit{ik}}(t)\text{,}{d}_{\mathit{lj}}(t)\text{,}{f}_{\mathit{ij}}(t)\text{,}{r}_j(t)\text{,}{e}_{\mathit{jk}}(t)\text{,}(i\text{,}l=1\text{,}2\text{,}\dots \text{,}n\text{;}j\text{,}k=1\text{,}2\text{,}\dots \text{,}m)$

are all continuous positive $\omega$ -periodic functions on R. By applying Mawhin continuation theorem and constructing a suitable Lyapunov function, they obtained sufficient conditions which guarantee the existence of a unique globally asymptotic stable positive $\omega$ -periodic solution of the systems (1.1) and (1.2). However, there are little real periodic systems in ecosystem, more systems are heteronomy ones. It is well known that the assumption of almost periodicity of the coefficients in (1.1) is a way of incorporating the time-dependent variability of the environment, especially when the various components of the environment are periodic with not necessary commensurate periods (e.g. seasonal effects of weather, food supplies, mating habits, and harvesting). It is very important to discuss almost periodic systems. Furthermore, as we know, the method used to investigate the positive $\omega$ -periodic solution of the non-linear ecosystem could not be used to investigate the almost periodic solution of the systems (1.1) and (1.2).

Recently, there have been many nice works on the almost positive periodic solutions of competitive model with periodic coefficients [14–19]. Chen [14] discussed the almost periodic solution of the non-autonomous two-species competitive model with stage structure. Alzabut, Nieto and Stamov [15] studied the existence and exponential stability of positive almost periodic solutions for a model of hematopoiesis. Wang and Dai [16] considered the almost periodic solution for n-species Lotka–Volterra competitive system with delay and feedback controls.

As we all known, investigating the almost periodic solutions of population dynamics model with time delays has more extensively practical application value. Zhao, Wang and Zhao [17] studied the almost periodic solution for a neutral multi-species Logarithmic population model. By employing Banach's fixed point theorem and using differential inequality technique, the authors presented some sufficient conditions ensuring the existence, uniqueness and globally exponential stability of almost periodic solution for the model. The results obtained extended and improved the earlier publications. Li and Zhang [18] considered a discrete Hematopoiesis model with a time delay, some sufficient conditions for the existence of a unique uniformly asymptotically stable positive almost periodic solution were obtained.

Motivated by the above papers, the aim of this paper is to obtain sufficient conditions for the existence of a unique globally asymptotic stable almost periodic solution of the systems (1.1) and (1.2), by utilizing the comparison theorem of the differential equation and constructing a suitable Lyapunov functional and applying the analysis technique of papers [10,12,15,16,19–21].

The remaining part of this paper is organized as follows: In Section 2, we will state several definitions and lemmas which will be useful in the proving of main results of this paper. In Section 3, by applying the theory of differential inequality, we present the permanence results for the systems (1.1) and (1.2). In Section 4, by constructing a suitable Lyapunov function, a set of sufficient conditions which ensure the global asymptotic stability of the systems (1.1) and (1.2) are obtained. In Section 5, we present some sufficient conditions which guarantee existence and uniqueness of almost periodic solution of the systems (1.1) and (1.2). In the end, two suitable examples together with their numeric simulations are given to illustrate our results.

Read full article

URL:

https://www.sciencedirect.com/science/article/pii/S0096300312010995

Find the Equilibrium Solution of the Modified Predator-prey System

Source: https://www.sciencedirect.com/topics/mathematics/prey-population

Share :