Numerical Solution of Non-Linear Prey-Predator System using Finite Elements Method

113 Numerical Solution of Non-Linear Prey-Predator System using Finite Elements Method Saad A. Manaa Ahmed F. Qasem College of Computer sciences and Mathematics University of Mosul, Iraq Received on: 21/02/2007 Accepted on: 16/04/2007 ABSTRACT A non-linear prey-predator system solved numerically by Galerkin method, and we compare these results with the results of Pius Peter Nyaanga[6] who used finite difference methods, we found that Galerkin finite elements method is faster than finite difference method to reach equilibrium state where the density for the prey ) , ( t x u and the predator ) , ( t x v are equals for all the values for x and t , also we found that

) , ( t x u and the predator ) , ( t x v are equals for all the values for x and t , also we found that Galerkin method converges towards the steady state solutions faster than finite difference method with less steps in time. Keywords: prey-predator system, Galerkin method, finite difference methods.

Introduction:
The finite element method is one of the most flexible tools available for solving engneering problems of the kind involved in analyzing the deformation of solids, the transfer of heat, the flow of fluids, or electrical problems, it can be applied to systems with virtually any geometric configuration or boundary conditions. [1] The finite element method has developed simultaneously with the increasing use of high-speed electronic digital computers and with the growing emphasis on numerical methods for engineering analysis. Although the method was originally developed for structural analysis, the general nature of the theory on which it is based has also made possible its successful application for solutions of problems in other fields of engineering [3].
A basic model for studying the interraction of species in population biology is the prey-predator model, the prey-predator model is a planar system representing the behavior of a population of prey, as fish, and a population of predators, as sharks [6].
Du and Lou [4] studied some uniqueness and exact multiplicity results for a predator-prey model, they consider positive solutions of a predator-prey model with diffusion and under homogeneous dirichlet boundary conditions,it turns out that a certain parameter in this model plays a very important role.
Tyutyunov et.al [2] studied directed movement of predators and the emergence of density-dependence in predator-prey models. Numerical analysis shows that, on the spatially aggregated scale, the average predator density adversely affects the individual consumption, leading to a non-linear predator-dependent trophic function.
Meng and Wang [5] studied asymptotic behavior of a predator-prey diffusion system with time delays. They found that the global asymptotic convergence is established by the upper-lower solutions and iteration method in terms of the rate constants of the reaction function independent of the time delays and the effect of diffusion.
In this paper, we will study the numerical solution for prey-predator system by Galerkin method and we will compare these results with the results of [6].

2.Mathematical model:
A basic model for studying the interraction of species in population biology is the prey-predator model, the prey-predator model is [6]:  is the rate of growth of the prey, ( ) 2 D is diffusion coefficient of the predator, ( )  is positive constant and represents the death rate of the predator, 2 1 ,k k are positive constants and they represent the rate of interraction between the prey and the predator which facilitates the killing of prey by predator. With initial conditions:

3.Numerical solution:
Galerkin's method has been discussed by several authors [8]. It is a means of obtaining an approximate solution to a differential equation. It does this by requiring that the error between the approximate solution and the exact solution be orthogonal to the functions used in the approximation. If we start with a differential equation (L is a differential operator) and approximate the solution by where  is a residual or error. Our desire is to make  as small as possible.
One way of accomplishing this objective is to require the integral  =  R i 0 dx N for each basis function i N . This integral mathematically states that the basic function must be orthogonal to the error over the region R [7]. Let * u and * v have the normal finite element form: where the interpolation function are [1]: Also integral over the element must involve the integrals of the pyramid functions [1]: it has been shown that the general integration formula is [1]: Galerkin method uses the approximating functions as weighting functions: E is the total number of elements. for the (ith) node in element (e) in the equation (11), the residual is: also for the (jth) node in element (e) in the equation (11), the residual is: for the (ith) node in element (e) in the equation (12), the residual is: also for the (jth) node in element (e) in the equation (11), the residual is: Now, we use the integrating by parts and use the neumann boundary conditions(4), the equation (13) becomes: (6) and (7), we get: and by using equation (9), the above equation becomes: from the time term (14) in the (ith) node, we get: Similarly, the equation (15) for the (jth) node in element (e) becomes: by the same way as in residual (ith), we get: this completes all of the residuals of the equation (11).
For equation (12) and by the same mannar we have: by using the equations (5),(6), (7) and (9) , we get: also for the (jth) node in element (e), the equation (19) becomes: by the same way as in residual (ith), we get:

5.Element matrices:
It is convenient to organize the residauls in element matrix form because a set of simultaneous nodal equations must be solved at each time step. Then the normal assembly procedure can be used. The terms involving nodal values at time step (m+1) are unknown and appear in the element matrix. All other terms are known and go in the element column [1].
They have the form: the particular components can now be evaluated.
The element residuals are written as: (22) and (24) for the space residual yields: also the element column is:        we can solve the previous systems (37) and (38) by using the Gaussian elimination method.

6.Numerical results:
In this section, we have solved the systems (37) and (38) with the neumann boundary conditions [6]:         The main conclusion, which we can draw from this result is that Galerkin finite element method is faster than finite difference method to reach the equilibrium state and to reach the steady state solutions with less steps in time, from Fiqure (1) and table (1) it is obvious that Galerkin method is faster than finite difference method to reach the equilibrium state where the values for u(x,t) are equal for all the values (x) and (t), and v(x,t) is also equal for all the values (x) and (t). Also Galerkin method is faster than finite difference method to make the predator v(x,t) reach the steady state solution (v=1), while the prey u(x,t) decreases to steady state (u=0) in far time steps, given that the rate of interraction ( )  of the prey, Galerkin method faster than finite difference method to reach the prey u(x,t) and the predator v(x,t) to the equilibrium state and to make the prey reach a steady state (u=1) and the predator to steady state (v=0) in far time steps as shown in table (2) and figure (2).