Optimal Parametric Iteration Method for Solving Multispecies Lotka-Volterra Equations

We apply an analytical method called the Optimal Parametric Iteration Method OPIM to multispecies Lotka-Volterra equations. By using initial values, accurate explicit analytic solutions have been derived. The method does not depend upon small parameters and provides us with a convenient way to optimally control the convergence of the approximate solutions. An excellent agreement has been demonstrated between the obtained solutions and the numerical ones. This new approach, which can be easily applied to other strongly nonlinear problems, is very effective and yields very accurate results.


Introduction
A substantial amount of research work has been invested in the study of nonlinear systems of differential equations.Systems of nonlinear differential equations arise in many scientific models such as biological systems and are used in various fields as engineering, chemistry, and ecology.In 1925, Lotka 1 developed the motion of an evolutionary system based on two fundamental changes, those involving matter between components of a system and those involving exchanges of energy 2 .Unlike being grounded in chemistry, Lotka believed that these ideas could be applied to any biological system.In 1926, Volterra 3 developed the wellknown mathematical models of multispecies interaction.These models, the predator, prey and competition models are known today as Lotka-Volterra models.The environment of the species can be influenced by the effect of food availability, weather conditions, temperature, mating habits, contact with predators, and other resource or physical environmental quantities.Lotka-Volterra equations describe variations of population densities of few species that compete for the same resources.However, the ecological system is often affected by environmental changes and other human activities.In many practical situations, it is often the case that one of the species maybe suffers a significant loss or increase in density for some reason at some transitory time slots.These models, for instance, they can describe the competing fish species which are exploited by human activities, can also describe the dynamics of normal and tumour cells in a changing environment under the effects of the chemotherapy.These models are also applicable in case we are interested in the existence and stability of tumor-free solution and how treatment affects the interaction of tumour and normal cells.Although simplistic, these few models are still used as the foundation for mathematical models in biology 4 .These models can also describe the time history of a biological system and are used in various fields as engineering, chemistry, biology, or mathematics 5 .In fact the Lotka-Volterra model is one of the most popular ones to demonstrate a simple nonlinear control system.The accurate solutions of the Lotka-Volterra equations may become a difficult task either if the equations are stiff or when the number of species is large.
Nonlinear analytical techniques for solving nonlinear problems have been dominated by the perturbation methods, which have found wide applications in engineering 6 .But perturbation methods have their limitations: perturbation techniques are based on the existence of a small parameter.A majority of nonlinear problems have no small parameters at all.So, it is necessary to develop a kind of new nonlinear analytical methods which does not require small parameters at all.
There exist same alternative methods, such as the variational iteration method 7 , the Adomian decomposition method 8 , a modified Lindstedt-Poincare method 9 , the homotopy analysis method 10 , the homotopy perturbation method 11 , and the optimal homotopy asymptotic method 12 .
In recent years, a growing interest towards the application of iterative techniques in nonlinear problems has appeared in science and engineering.In 1987 Mickens 13 proposed an iterative scheme for nonlinear problems.In 2002 Lim and Wu 14 , in 2006 Hu 15 , and in 2008 Chen and Liu 16 proposed modified iteration procedures.Also He 17 proposed some iterative and asymptotic methods for nonlinear problems and later, Marinca and Heris ¸anu 18 proposed in 2006 a new iteration method by combining Mickens' and He's iteration methods.Ramos 19 investigated other iterative techniques for nonlinear differential equations.These methods are valid not only for a small parameter, but also for very large parameters and have been used to solve nonlinear both conservative and nonconservative oscillators.
In this paper we propose a new approach to find analytic approximate solution of multispecies Lotka-Volterra models using a new iterative procedure, namely, the Optimal Parametric Iteration Method OPIM .The efficiency of the proposed procedure is proved since an accurate solution is explicitly analytically achieved in an iterative manner after only one iteration.This new approach involves the presence of a finite number of initially unknown parameters, which are optimally determined, providing a rigorous way to control the convergence of the solutions 20 .Different methodologies have been applied to study Lotka-Volterra models 21-25 .Using this new approach, the obtained approximate solutions rapidly converge to the exact solution.Some alternative ways for mathematical biology problems are presented in 26-28 .

Formulation and Solution Approach
In order to introduce the OPIM framework for the most general form of systems of nonlinear ODEs, we consider the following system: where dot denotes derivative with respect to time t, f i , i 1, 2, . . ., s are nonlinear functions.
In the frame of OPIM we define the following iteration formulas:  A special case of this general system is the Lotka-Volterra model with three species, which is described by the following system of nonlinear differential equation 29, 30 : where α and β are known parameters not necessary small.Initial conditions are x 0 a, y 0 b, z 0 c.

2.4
In this special case, it is clear that s 3 and

2.5
With the notations 2.5 , the iteration formulas 2.2 can be written in the form , n e , and p e must follow the terms appearing in 2.6 .Therefore, we try to choose the auxiliary functions so that in 2.6 the products m e F x are of the same shape with the function F x , where F ∈ {f, g, h}.The constants C i , C j , C k which appear in the auxiliary functions m e , n e , and p e can be optimally determined via various methods, for example, putting the conditions that the residual functionals J 1 , J 2 , and J 3 given, respectively, by be minimum, that is: taking into account also the initial conditions.Other methods suitable to determine the constants C j are the Galerkin method, the collocation method, and so on.
On the other hand, the initial approximations x 0 , y 0 , and z 0 and also the auxiliary functions m e , n e , and p e are not unique.
In short, the basic idea of the proposed procedure OPIM are the construction a new iteration scheme 2.6 , the involvement of the auxiliary functions, and the convergencecontrol constants C 1 , C 2 , . . .which lead to an excellent agreement of the approximate solutions with the exact ones.

Results and Discussions
In this section we present an example to show the efficiency of the method described in the previous section for solving 2.3 .
Taking into account the initial conditions, we choose the initial approximations as functions of the form x 0 ae −t , y 0 be −t , z 0 ce −t .

3.1
For n 0 into 2.6 we can construct the following iteration formulas:

3.2
For the auxiliary functions m i , n i , p i , i 1, 2, 3 we choose the expressions Alternatively, we can consider either other expressions for the initial approximations and the auxiliary function, follows as: − βcC 14 t βcC 16 e −3t .

3.5
The solutions of 3.5 with the initial conditions can be written in the form

3.7
By substituting 3.7 into 2.3 it results in the residuals:

3.8
Making collocation in the arbitrary points t i , i 1, 2, . . ., 6 3.9 we obtain the optimal values of the constants C 1 , C 2 , . . ., C 18 and therefore the solution 3.7 in the first approximation is well determined.
In the case when a 0. It is easy to verify the accuracy of the obtained solution if we graphically compare the analytical solution with the numerical one.Figures 1, 2, and 3 show the comparison between the present solutions and the numerical integration results obtained by a fourth-order Runge-Kutta method.
It can be seen from Figures 1-3 that the solution obtained by OPIM is in very good agreement with numerical integration results.

Conclusions
In this work the Optimal Parametric Iteration Method OPIM is employed to propose new analytic solutions for some nonlinear differential equations.The validity of the method is   illustrated on the Lotka-Volterra model with three species.Our procedure is independent of the presence of small parameters.Our construction of iterations is different from other known iteration techniques.In short, the basic new ideas of our method are related to the construction of a new iteration scheme which involves some auxiliary functions m e , n e , and p e , e 1, 2, 3 whose parameters C 1 , C 2 , . . .viz. the convergence-control constants lead to an excellent agreement of the approximate solutions with the exact ones.Our procedure is very effective, explicit, and accurate for nonlinear approximations, rapidly converging to the exact solution after only one iteration.This procedure provides a simple but rigorous way to control and adjust the convergence of the solutions by optimally determining the parameters C 1 , C 2 , . ... Our method gives analytic solutions valid globally in time unlike other known methods, for instance.Adomian decomposition method, which unfortunately does not guarantee analytic solutions valid globally in time as proved by Rèpaci 31 .
y n , z n n 1 t, C i f y x n , y n , z n p 1 t, C i f z x n , y n , z n ẏn 1 g x n , y n , z n m 2 t, C j g x x n , y n , z n n 2 t, C j g y x n , y n , z n p 2 t, C j g z x n , y n , z n żn 1 h x n , y n , z n m 3 t, C k h x x n , y n , z n n 3 t, C k h y x n , y n , z n p 3 t, C k h z x n , y n , z n , x ∂f/∂x, n 0, 1, 2 . .., and C i , C j , C k , i, j, k 1, 2, . . .are unknown constants at this moment.There are many possibilities to choose the auxiliary functions m e t, C i , n e t, C j and p e t, C k , e 1, 2, 3. Basically, the shape of m e − αaC 4 e −t C 1 − βaC 5 t C 2 − βaC 6 − a 2 αab βac e −2t − 2a αb βc C 1 t C 2 e −3t ẏ1 C 9 t −βbC 7 t b − βbC 8 C 10 e −t − βa 2b αc C 9 t βab b 2 αbc βa 2b αc C 10 e −2t − αbC 11 t αbC 12 e −3t ż1 C 17 t C 18 c e −t − αa βc 2c C 17 αcC 13 t αac βbc c 2 αcC 14 αa βc 2c C 18 e −2t 2, b 0.3, c 0.5, α 0.1, β 0.1, t 1 0.3, t 2 0.7, t 3 1, t 4 2, t 5 3, t 6 5, from 3.9 we obtain