Numerical Simulation for the Treatment of Nonlinear Predator–Prey Equations by Using the Finite Element Optimization Method

: This article aims to introduce an efﬁcient simulation to obtain the solution for a dynamical– biological system, which is called the Lotka–Volterra system, involving predator–prey equations. The ﬁnite element method (FEM) is employed to solve this problem. This technique is based mainly upon the appropriate conversion of the proposed model to a system of algebraic equations. The resulting system is then constructed as a constrained optimization problem and optimized in order to get the unknown coefﬁcients and, consequently, the solution itself. We call this combination of the two well-known methods the ﬁnite element optimization method (FEOM). We compare the obtained results with the solutions obtained by using the fourth-order Runge–Kutta method (RK4 method). The residual error function is evaluated, which supports the efﬁciency and the accuracy of the presented procedure. From the given results, we can say that the presented procedure provides an easy and efﬁcient tool to investigate the solution for such models as those investigated in this paper.


Introduction
Nonlinear differential equations have kept attracting the interest of many researchers. Some authors have realized the need for finding new mathematical models of many real-life problems in such different areas as, for example, fluid mechanics, viscoelasticity, biology, physics, and engineering. Most nonlinear differential equations do not have an exact solution. Therefore, approximate and numerical techniques must be used (see [1][2][3][4]).
For relevant further details about differential equations, together with their important properties and the tools and techniques used for them, one can see the works [5] (see also [6] for various developments involving ordinary differential equations as well as analogous fractional-order differential equations).
After the above-mentioned works, many researchers have successfully applied various numerical methods in this field (see, for example, [7]). Among these numerical methods, the finite element method (FEM) (see [8]) has some distinct advantages for handling this class of problems in which the coefficients for the solution can easily be found to exist after 1.
FEM allows for easier modeling of complex geometrical and irregular shapes.

2.
FEM can be adapted to obtain the solution of the problem with a high degree of accuracy to decrease the need for physical prototypes in the design process. 3.
FEM is remarkably useful for certain time-dependent simulations, such as crash simulations in which deformations in one area depend upon deformations in another area.
In this paper, we apply the finite element optimization method (FEOM) to obtain the solution for the given problem. The outline of the paper is as follows: • Section 2 is devoted to the presentation of the formulation of the Lotka-Volterra system; • Section 3 is devoted to the construction and prediction of the solution by FEOM; • Section 4 is devoted to the results and discussion; • Section 5 presents the conclusions of this investigation.

The Model of the Lotka-Volterra System (with Two Predators and One Prey)
Predator-prey equations are a pair of nonlinear differential equations of the first order that are frequently used to describe the dynamics of biological systems in which two species interact with each other; one is a predator and the other is its prey. The studied model was presented by Samardzija [9], who proposed a two-predator and one-prey concept of the Lotka-Volterra system.
In our study, the proposed model for the predator-prey equations is defined as follows: and together with the following initial conditions: x(0) = x 0 , y(0) = y 0 and z(0) = z 0 , where a prime denotes differentiation with respect to t, y(t) and z(t) are the numbers of the predators at time t, x(t) is the number of their prey at time t, and λ j (j = 1, . . . , 7) are the parameters that represent the interaction between the three species involved (see also [10,11]). In our present work, we solve this problem in the interval (0, t f ) = (0, 8). For more details about this model, see the earlier work [9].

Procedure of Solution by the Finite Element Method (FEM)
The FEM is a powerful technique for solving differential equations and integral equations, as well as their fractional-order analogs. The basic concept is that the whole domain is divided into smaller elements of finite dimensions, called Finite Elements. It is one of the most versatile numerical techniques in modern engineering analysis and is employed to study diverse problems in such areas as, for example, fluid mechanics [12], heat transfer [13], chemical processing [14], solid mechanics [15], electrical systems, and many other fields.
The analysis of the finite element method is implemented through the following steps (see, for details, [16]): 1.
Finite-element discretization: The whole domain is divided into a finite number of sub-domains, which is called the discretization of the domain. Each sub-domain is called an element. The collection of elements is called the finite-element mesh.

2.
Generation of the element equations: • From the mesh, a typical element is isolated and the variational formulation of the given problem over the typical element is then constructed. • An approximate solution of the variational problem is assumed and the element equations are formulated by substituting this solution into the above system. • The element matrix, which is called the stiffness matrix, is constructed by using the element interpolation functions.

3.
Assembly of the element equations: The algebraic equations so obtained are assembled by imposing the inter-element continuity conditions. This yields a large number of algebraic equations known as the global finite element model, which governs the whole domain. 4.
Imposition of the boundary conditions: The essential and natural boundary conditions are imposed on the assembled equations. 5.
Solution a system of algebraic equations: The obtained system of algebraic equations can be solved by any one of the appropriate numerical techniques.

Variational Formulation of the System
The variational form associated with Equations (1)-(3) over a typical linear element (t e , t e+1 ) is given by the following: where φ 1 , φ 2 and φ 3 are arbitrary test functions and may be viewed as the variation in x(t), y(t) and z(t), respectively.

Finite Element Formulation of the Model
The finite element model may be obtained from the above equations by substituting finite element approximations of the following form: In our computations, the shape functions of a typical element (t e , t e+1 ) are taken to be as follows: Linear element: and t e t t e+1 .
Quadratic element: and t e t t e+1 .
The finite element model of the equations thus formed is given by the following: where [K rs ] and [b r ] (r, s = 1, 2, 3) are defined as follows: in whichx Each element matrix is of the order 6 × 6. The entire flow domain is divided into a set of 800 line elements and, following the assembly of all of the element equations, a matrix of order 2403 × 2403 is generated. The resulting system of equations is strongly nonlinear and recourse must be made to a robust iterative scheme to solve it. The system is linearized by incorporating the functionsx,ȳ andz, which are assumed to be known. After applying the given boundary conditions, only a system of 2398 equations remains for finding the solution, which is performed by using a robust Gauss elimination method while maintaining an accuracy of 0.0001. A convergence criterion based on the relative difference between the current and previous iterations is employed. When these differences reach the desired accuracy, the solution is assumed to have converged and the iterative process is then terminated. The Gaussian quadrature is implemented for solving the integrations. However, the suitability of the shape functions varies from problem to problem. Due to the simple and efficient use in computations, linear, as well as quadratic, shape functions are used in the present problem. However, it is observed that the results do not vary considerably, indicating that both elements provide approximately the same accuracy.

Numerical Simulation
Here, in this section, we present a numerical test example to verify the accuracy and the quality of the presented scheme. With this aim in view, we consider the systems (1) to (4) with different values of the parameters λ j (j = 1, . . . , 7), as well as the initial conditions for x 0 , y 0 and z 0 . Figures 1-4 present the numerical simulation for the studied model by employing the technique that we have used in this paper.
In Figure 1, we present the behavior of the approximate solution by the proposed method (dotted line) and the fourth-order Runge-Kutta method (continuous line) in the interval (0, 8) with the initial conditions x 0 = 0.5, y 0 = 1 and z 0 = 1 and with the parameters given by the following: In Figure 2, we present the behavior of the approximate solution via distinct values of the initial conditions in the interval (0, 8) and with the parameters given by the following: In this situation, we take the following three cases: i. x 0 = y 0 = 1, z 0 = 2; (a) for the solution x(t); ii. x 0 = y 0 = 2, z 0 = 4; (b) for the solution y(t); iii. x 0 = y 0 = 3, z 0 = 6; (c) for the solution z(t).
In addition, in Figure 4, the residual error function (REF) is introduced at the initial conditions x 0 = y 0 = 1, and z 0 = 2, with the parameters given by the following: From these results, we can note that the behavior of the numerical solution depends on the values in the initial conditions and the included parameters λ j (j = 1, . . . , 7). This confirms that the used method has been implemented in a good way for solving the proposed problem.    In order to validate our numerical solutions, we present a comparison of the REF in Table 1 at the initial conditions x 0 = y 0 = 2, and z 0 = 3, and with the parameters given by the following: This comparison shows the thoroughness of the proposed method in this article.

Conclusions
The proposed model is solved numerically by employing the finite element optimization method (FEOM). The resulting numerical solutions for the studied problem confirm that the presented technique is quite suitable to study this problem effectively. We can confirm also that, if we add extra terms from the series of the approximate solution, the errors will decrease. Finally, the numerical solutions with different values of the order n of approximation and the residual error function are computed to illustrate the validity of the proposed technique. The derived results also show that the given procedure provides an efficient tool to investigate the numerical solution for such models. On the other hand, in capturing this numerical analysis, our work may provide stronger physical meanings for performing future theoretical and computational analysis on this subject. All numerical results in this paper were obtained by using the Mathematica software.