A Numerical Study of Several Species Population Models

Francisco J. Sánchez-Bernabe¹, Maria del Rosario Escalona-Magdaleno²

¹Department of Mathematics, Universidad Autónoma Metropolitana Iztapalapa, Ciudad de México, Mexico.
²Instituto de Educación Media Superior de la Ciudad de México, Plantel Miravalle, Zapote y Mirador, Mexico.

DOI: 10.4236/jamp.2023.1112251                

Abstract

This work considers a special case of Lotka-Volterra equations, which means that in the system of two ordinary differential equations, we take the four parameters equal to one. The reason is that we want just to illustrate the procedure to reduce that system to only one ordinary differential equation, such that we know its analytical solution. This idea will be applied to study the relations between a system of three ordinary differential equations, and a couple of partial differential equations. Lotka-Volterra equations are solved numerically by a fourth-order predictor-corrector method, which is initialized by an improved Euler method with a rather small time step because it is only a second-order algorithm. Then, it is proposed a model with three species, defined by ordinary differential equations.

Keywords

Lotka-Volterra Equations, Adams-Bashfort Predictor, Adams-Moulton Corrector, Improved Euler Method, GeoGebra, Matlab

1. Introduction

Lotka-Volterra model [1] [2] is a well known system of ordinary differential equations that model interactions between two species, that it is numerically approximated by an Adams-Bashfort predictor and a corresponding Adams-Moulton corrector scheme [3] . Since an implicit exact solution for Lotka-Volterra equations is known, the numerical values may be compared with those defined for the exact solution to verify that predictor-corrector method produces good results [4] . Next, taking like a motivation Lotka-Volterra equations, we consider a couple of partial differential equations that have a similar structure. Then, a system of three ordinary differential equations [5] is related to that system of two coupled partial differential equations that is studied from a geometric point of view, with MATLAB [6] . It is shown that curves that are solutions of the system of ordinary differential equations are contained in the surface defined by a equation that it turns to be a solution of both partial differential equations.

2. Statement of the Problem

Given the system of ordinary differential equations

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x-xy\mathrm{,}\\ \frac{\text{d}y}{\text{d}t}=xy-y\end{array}$ (1)

where we consider the following initial conditions

$\{\begin{array}{l}x\left(0\right)\equiv x_0=e,\\ y\left(0\right)\equiv y_0=1\end{array}$ (2)

by applying chain rule

$\frac{\text{d}y}{\text{d}t}=\frac{\text{d}y}{\text{d}x}\frac{\text{d}x}{\text{d}t},$ (3)

then, solving for $\text{d}y/\text{d}x$

$\frac{\text{d}y}{\text{d}x}=\frac{\text{d}y/\text{d}t}{\text{d}x/\text{d}t},$ (4)

by substituting each Equation of (1) in (4)

$\frac{\text{d}y}{\text{d}x}=\frac{y}{x}\frac{x-1}{1-y}\mathrm{.}$ (5)

Next step is to introduce an equation that satisfies initial conditions (2), that is

$x+y=\ln xy+e\mathrm{.}$ (6)

Let us check that also satisfies the differential Equation (5). By taking derivatives with respect to x on both sides of (6)

$1+\frac{\text{d}y}{\text{d}x}=\frac{1}{xy}\left(x\frac{\text{d}y}{\text{d}x}+y\right),$ (7)

in other words

$1+\frac{\text{d}y}{\text{d}x}=\frac{1}{y}\frac{\text{d}y}{\text{d}x}+\frac{1}{x}\mathrm{,}$ (8)

that is

$\frac{1}{y}\frac{\text{d}y}{\text{d}x}-\frac{\text{d}y}{\text{d}x}=1-\frac{1}{x},$ (9)

in an equivalent way

$\frac{\text{d}y}{\text{d}x}\left(\frac{1}{y}-1\right)=\frac{x-1}{x},$ (10)

furthermore

$\frac{\text{d}y}{\text{d}x}\left(\frac{1-y}{y}\right)=\frac{x-1}{x}\mathrm{,}$ (11)

in consequence, we obtain the result by solving in (11) for the derivative.

On the other hand, system of diferential Equation (1) can be written in the following way

$\{\begin{array}{l}{x}^{\prime }\left(t\right)=f\left(t,x,y\right),\\ y^{\prime }\left(t\right)=g\left(t,x,y\right)\end{array}$ (12)

besides initial conditions (2).

Improved Euler method applied to each equation gives

$\{\begin{array}{l}{x}_{n+1}=x_n+\frac{h}{2}\left(f\left(t_n,x_n,y_n\right)+f\left(t_{n+1},x_n+hf\left(t_n,x_n,y_n\right),y_n\right)\right),\\ y_{n+1}=y_n+\frac{h}{2}\left(g\left(t_n,x_n,y_n\right)+g\left(t_{n+1},x_n,y_n+hg\left(t_n,x_n,y_n\right)\right)\right)\end{array}$ (13)

since

$\{\begin{array}{l}f\left(t_n,x_n,y_n\right)=x_n-x_n{y}_n,\\ g\left(t_n,x_n,y_n\right)=x_n{y}_n-y_n\end{array}$ (14)

by substituting the first Equation of (14)

$f\left(t_{n+1},x_n+hf\left(t_n,x_n,y_n\right),y_n\right)=f\left(t_{n+1},x_n+h\left(x_n-x_n{y}_n\right),y_n\right),$ (15)

that is

$f\left(t_{n+1}\mathrm{,}{x}_n+hf\left(t_n\mathrm{,}{x}_n\mathrm{,}{y}_n\right)\mathrm{,}{y}_n\right)=f\left(t_{n+1}\mathrm{,}\left(1+h\right)x_n-h{x}_n{y}_n\mathrm{,}{y}_n\right)\mathrm{,}$ (16)

again, by the first Equation of (14)

$f\left(t_{n+1},x_n+hf\left(t_n,x_n,y_n\right),y_n\right)=\left(\left(1+h\right)x_n-h{x}_n{y}_n\right)\left(1-y_n\right),$ (17)

in an equivalent way

$f\left(t_{n+1},x_n+hf\left(t_n,x_n,y_n\right),y_n\right)=x_n\left(1+h-h{y}_n\right)\left(1-y_n\right),$ (18)

by sustituting (18) into the first Equation of (13)

$x_{n+1}=x_n+\frac{h}{2}\left(x_n-x_n{y}_n+x_n\left(1+h-h{y}_n\right)\left(1-y_n\right)\right),$ (19)

that is

$x_{n+1}=x_n+x_n\frac{h}{2}\left(1-y_n+\left(1+h-h{y}_n\right)\left(1-y_n\right)\right)\mathrm{,}$ (20)

in consequence

$x_{n+1}=x_n+x_n\left(1-y_n\right)\frac{h}{2}\left(2+h\left(1-y_n\right)\right)\mathrm{.}$ (21)

Similarly, from the second Equation of (13)

$y_{n+1}=y_n+\left(x_n-1\right)y_n\frac{h}{2}\left(2+h\left(x_n-1\right)\right)\mathrm{.}$ (22)

From initial conditions (2), we choose $h=1/512$ , so that, with $n=0$ , we obtain

$\{\begin{array}{l}{x}_1=e\approx 2.718281828459,\\ y_1=1+\left(e-1\right)\frac{1}{1024}\left(2+\frac{1}{512}\left(e-1\right)\right)\approx 1.0033616\end{array}$ (23)

in the following step

$\{\begin{array}{l}{x}_2\approx \mathrm{2.71826398,}\\ y_2\approx 1.0067346\end{array}$ (24)

applying again, improved Euler method

$\{\begin{array}{l}{x}_3\approx \mathrm{2.71822822,}\\ y_3\approx 1.0101188\end{array}$ (25)

Calculated points and initial conditions iniciales, are denotated in the following way: $A=\left(x_0\mathrm{,}{y}_0\right)=\left(e\mathrm{,1}\right)$ , $B=\left(x_1\mathrm{,}{y}_1\right)$ , $C=\left(x_2\mathrm{,}{y}_2\right)$ , and $D=\left(x_3\mathrm{,}{y}_3\right)$ . These points are shown in Figure 1.

Up to now, we have considered a time step given by $h=1/512$ that is relatively small because improved Euler is a low order method. But, since we already have, beside the initial values, three ordered pairs of approximated values, it is possible to apply fourth order predictor scheme

$\begin{array}{c}{\tilde{x}}_{m+1}=x_m+\frac{h}{24}(55f\left(x_m,y_m\right)-59f\left(x_{m-1},y_{m-1}\right)\\ +37f\left(x_{m-2},y_{m-2}\right)-9f\left(x_{m-3},y_{m-3}\right)),\end{array}$ (26)

similarly

$\begin{array}{c}{\tilde{y}}_{m+1}=y_m+\frac{h}{24}(55g\left(x_m,y_m\right)-59g\left(x_{m-1},y_{m-1}\right)\\ +37g\left(x_{m-2},y_{m-2}\right)-9g\left(x_{m-3},y_{m-3}\right)),\end{array}$ (27)

then, a corrector method of fifth order is also considered

Figure 1. Exact solution (6), and points obtained by (21), (22), and (26), (27), (28), (29).

$\begin{array}{c}{x}_{m+1}=x_m+\frac{h}{720}(251f\left({\tilde{x}}_{m+1},{\tilde{y}}_{m+1}\right)+646f\left(x_m,y_m\right)-246f\left(x_{m-1},y_{m-1}\right)\\ +106f\left(x_{m-2},y_{m-2}\right)-19f\left(x_{m-3},y_{m-3}\right)),\end{array}$ (28)

besides

$\begin{array}{c}{y}_{m+1}=y_m+\frac{h}{720}(251g\left({\tilde{x}}_{m+1},{\tilde{y}}_{m+1}\right)+646g\left(x_m,y_m\right)-246g\left(x_{m-1},y_{m-1}\right)\\ +106g\left(x_{m-2},y_{m-2}\right)-19g\left(x_{m-3},y_{m-3}\right)).\end{array}$ (29)

these schemes, are known as Adams-Bashfort and, Adams-Moulton. In this way, we obtain from the four points $A\mathrm{,}B\mathrm{,}C$ y $D$ , with $h=1/512$ the following predictions

$\{\begin{array}{l}{\tilde{x}}_4\approx \mathrm{2.71816549,}\\ {\tilde{y}}_4\approx 1.01351437\end{array}$ (30)

and then, the corrections

$\{\begin{array}{l}{x}_4\approx \mathrm{2.71816460,}\\ y_4\approx 1.01359883\end{array}$ (31)

By defining $E_p=\left(x_4\mathrm{,}{y}_4\right)$ , from the points $B\mathrm{,}C\mathrm{,}D$ , and $E_p$ , we calculate the predictions

$\{\begin{array}{l}{\tilde{x}}_5\approx \mathrm{2.71808279,}\\ {\tilde{y}}_5\approx 1.01700624\end{array}$ (32)

and then, the corrections

$\{\begin{array}{l}{x}_5\approx \mathrm{2.71808191,}\\ y_5\approx 1.01643729\end{array}$ (33)

Next, we define $F_p=\left(x_5\mathrm{,}{y}_5\right)$ , and from the points $C\mathrm{,}D\mathrm{,}{E}_p$ , and $F_p$ , we obtain the predictions

$\{\begin{array}{l}{\tilde{x}}_6\approx \mathrm{2.71798206,}\\ {\tilde{y}}_6\approx 1.02050936\end{array}$ (34)

and then, the corrections

$\{\begin{array}{l}{x}_6\approx \mathrm{2.71798025,}\\ y_6\approx 1.02059439\end{array}$ (35)

Now, by defining $G_p=\left(x_6\mathrm{,}{y}_6\right)$ , from the points $A\mathrm{,}C\mathrm{,}{E}_p$ and $G_p$ , we are going to double the value of $h=1/256$ because we are using a higher order method than improved Euler. So, by applying Adams-Bashfort we obtain

$\{\begin{array}{l}{\tilde{x}}_7\approx \mathrm{2.71772389,}\\ {\tilde{y}}_7\approx 1.02746678\end{array}$ (36)

and then, we make a correction with Adams-Moulton

$\{\begin{array}{l}{x}_7\approx \mathrm{2.71772123,}\\ y_7\approx 1.02763628\end{array}$ (37)

Next step: we define $H_p=\left(x_7\mathrm{,}{y}_7\right)$ , we consider the points $C\mathrm{,}{E}_p\mathrm{,}{G}_p$ , and $H_p$ , we get the following predictions with Adams-Bashfort

$\{\begin{array}{l}{\tilde{x}}_8\approx \mathrm{2.71739061,}\\ {\tilde{y}}_8\approx 1.03455447\end{array}$ (38)

and then, by applying Adams-Moulton

$\{\begin{array}{l}{x}_8\approx \mathrm{2.71738546,}\\ y_8\approx 1.03472549\end{array}$ (39)

Finaly, we define $K_p=\left(x_8\mathrm{,}{y}_8\right)$

In Figure 1 it is shown the exact solution (6) of studied case of Lotka-Volterra equations, and the points obtained by numerical methods.

3. A Model for Three Species

We take as motivation, the Lotka-Volterra equations to propose the couple of partial differential equations

$\{\begin{array}{l}\frac{\partial z}{\partial x}=\frac{z}{x}\left(\frac{x-1}{1-z}\right),\\ \frac{\partial z}{\partial y}=\frac{z}{y}\left(\frac{y-1}{1-z}\right)\end{array}$ (40)

then, we define the following system of ordinary differential equations

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{x}{x-1},\\ \frac{\text{d}y}{\text{d}t}=\frac{y}{y-1},\\ \frac{\text{d}z}{\text{d}t}=\frac{2z}{1-z}\end{array}$ (41)

from the first Equation of (41)

$\frac{x-1}{x}\frac{\text{d}x}{\text{d}t}=1,$ (42)

that is

$\left(1-\frac{1}{x}\right)\frac{\text{d}x}{\text{d}t}=\mathrm{1,}$ (43)

then, we obtain

$\text{d}x-\frac{\text{d}x}{x}=\text{d}t\mathrm{,}$ (44)

by integrating

$x-\ln x=t+c_1\mathrm{,}$ (45)

similarly, from the second Equation of (41)

$y-\ln y=t+c_2\mathrm{,}$ (46)

next, from the third differential Equation of (41), we get

$\frac{1-z}{2z}\frac{\text{d}z}{\text{d}t}=\mathrm{1,}$ (47)

that is

$\left(\frac{1}{z}-1\right)\frac{\text{d}z}{\text{d}t}=2$

$\left(\frac{1}{z}-1\right)\frac{\text{d}z}{\text{d}t}=2,$ (48)

implying that

$\frac{\text{d}z}{z}-\text{d}z=2\text{d}t,$ (49)

by integrating

$\ln z-z=2t+{\tilde{c}}_3\mathrm{,}$ (50)

as a consequence, we have the following equations

$\{\begin{array}{l}x-\ln x=t+c_1,\\ y-\ln y=t+c_2,\\ z-\ln z=-2t+c_3\end{array}$ (51)

by adding each corresponding side of Equation of (51), we obtain

$x+y+z-\ln x-\ln y-\ln z=\tilde{c},$ (52)

that can be expressed in the following way

$x+y+z-\frac{1}{2}\ln x-\frac{1}{2}\ln y-\frac{1}{2}\ln x-\frac{1}{2}\ln z-\frac{1}{2}\ln y-\frac{1}{2}\ln z=\tilde{c}\mathrm{,}$ (53)

multiplicating by 2

$2x+2y+2z-\ln x-\ln y-\ln x-\ln z-\ln y-\ln z=c\mathrm{,}$ (54)

that is

$2x+2y+2z-\ln xy-\ln xz-\ln yz=c.$ (55)

Let us consider the value $c=6e-e$ , that is

$2x+2y+2z-\ln xy-\ln xz-\ln yz=6e-e.$ (56)

Then, in the system (41), we will have the following initial conditions

$\{\begin{array}{l}x\left(0\right)=e,\\ y\left(0\right)=e,\\ z\left(0\right)=e\end{array}$ (57)

now, let us consider the first Equation of (51)

$x-\ln x=t+c_1,$ (58)

then we substitute $x=e$ , if $t=0$ , to obtain

$e-\ln e=0+c_1\mathrm{,}$ (59)

that is

$e-1=0+c_1\mathrm{,}$ (60)

by solving for the constant

$c_1=e-\mathrm{1,}$ (61)

in consequence, from (58)

$x-\ln x=t+e-1.$ (62)

Similarly, the second Equation of (51), now becomes

$y-\ln y=t+e-1.$ (63)

Next, we consider the third Equation of (51)

$z-\ln z=-2t+c_3,$ (64)

by imposing the initial condition

$e-\ln e=-2\left(0\right)+c_3,$ (65)

where, the constant is given by

$c_3=e-\mathrm{1,}$ (66)

in consequence

$z-\ln z=-2t+e-1.$ (67)

Now, by adding the three solutions (62), (63), and (67), variable t is canceled

$x+y+z-\ln x-\ln y-\ln z=3e-3.$ (68)

that can be written in the following way

$x+y+z-\frac{1}{2}\ln x-\frac{1}{2}\ln y-\frac{1}{2}\ln x-\frac{1}{2}\ln z-\frac{1}{2}\ln y-\frac{1}{2}\ln z=3e-3.$ (69)

multiplicating by 2

$2x+2y+2z-\ln x-\ln y-\ln x-\ln z-\ln y-\ln z=6e-6.$ (70)

that coincides with (56) after applying logarithms properties (see Figure 2).

4. Main Results

Let us take partial derivatives with respect to x on both sides of (56)

$2+2\frac{\partial z}{\partial x}-\frac{y}{xy}-\frac{x{z}_x+z}{xz}-\frac{y{z}_x}{yz}=\mathrm{0,}$ (71)

that is

$2\frac{\partial z}{\partial x}-\frac{1}{x}-\frac{x{z}_x+z}{xz}-\frac{z_x}{z}=-\mathrm{2,}$ (72)

Figure 2. Surface defined by Equation (56).

multiplying both sides of (72) by xz

$2xz\frac{\partial z}{\partial x}-z-\left(x{z}_x+z\right)-x{z}_x=-2xz\mathrm{,}$ (73)

in other words

$2xz\frac{\partial z}{\partial x}-2z-2x\frac{\partial z}{\partial x}=-2xz\mathrm{,}$ (74)

from where, we obtain

$xz\frac{\partial z}{\partial x}-x\frac{\partial z}{\partial x}=z-xz\mathrm{,}$ (75)

that is

$x\frac{\partial z}{\partial x}\left(z-1\right)=z\left(1-x\right)\mathrm{,}$ (76)

in consequence, we recover the first Equation of system (40). Therefore, the surface defined by (56) satisfies the first partial differential Equation of system (40).

Similarly, it can be check that the Equation of the surface (56), also satisfies the second partial differential Equation of the system (40) of partial differential equations.

5. Conclusion and Suggestions

Lotka-Volterra equations are solved by the improved Euler method and fourth-order predictor-corrector scheme Adams-Bashfort and Adams-Moulton. Then a three-species model is defined by a couple of partial differential equations, which are transformed into a system of three ordinary differential equations, which are not coupled between them. It was shown that solutions of the system of three ordinary differential equations are contained in the surface defined by the solution of the couple of partial differential equations. Future work will be addressed to study systems of three ordinary differential equations [7] [8] , which are actually coupled. This means that each right-hand side of these differential equations depends on at least two independent variables.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to assistant editor Nancy Ho.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Evans, C. and Findley, G. (1999) Analytic Solutions to a Family of Lotka-Volterra Related Differential Equations. Journal of Mathematical Chemistry, 25, 181-189. https://doi.org/10.1023/A:1019184417025
[2]	Draper, P.W. (2017) Lotka-Volterra Predator-Prey Analytic and Numerical Models. Master’s Thesis, Laurentian University Sudbury, Ontario.
[3]	Sauer, T. (2012) Numerical Analysis. 2nd Edition, Pearson Education Inc., Boston.
[4]	Hall, J. and Lingefjärd, T. (2016) Mathematical Modeling: Applications with GeoGebra. J. Wiley & Sons, New Jersey.
[5]	Reddy, B.R. (2013) Stability Analysis of Mathematical Syn-Ecological Model Comprising of Prey-Predator, Host-Commensal, Mutualism and Neutral Pairs II (Three of the Four Species Are Washed out States). Journal of Engineering, Science and Mathematics, 2, 119-129.
[6]	Attaway, D.C. (2018) MATLAB: A Practical Introduction to Programming and Problem Solving. Butterworth-Heinemann, Oxford.
[7]	Escalona-Magdaleno, M.D.R. (2017) Modelos Poblacionales, Undergraduate Research Project, Universidad Autónoma Metropolitana Iztapalapa.
[8]	Fonseca-Riquelme, W.A. (2015) Modelos Matemáticos en Dinámica de Poblaciones. Universidad del Bío-Bío, Chile.