Series Solution of the Multispecies Lotka-Volterra Equations by Means of the Homotopy Analysis Method

The time evolution of the multispecies Lotka-Volterra system is investigated by the homotopy analysis method HAM . The continuous solution for the nonlinear system is given, which provides a convenient and straightforward approach to calculate the dynamics of the system. The HAM continuous solution generated by polynomial base functions is of comparable accuracy to the purely numerical fourth-order Runge-Kutta method. The convergence theorem for the three-dimensional case is also given.


Introduction
The simplest model of predator-prey interactions was developed independently by Lotka  where the a's, b's, and c's are constants.System 1.2 -1.3 has a wide applicability to a variety of different physical 5 , chemical 6 , and biological problems 7 .
In the study of nonlinear systems of differential equations such as the Lotka-Volterra equation, analytical solutions are usually unknown.In this case, in order to analyze the behavior of the system, one usually resorts to numerical integration techniques, such as the Runge-Kutta method 8 , or perturbation techniques 9 .The problem with purely numerical technique like the Runge-Kutta method is that it does not give a functional form of the solution to the problem at hand, which is often useful if we need to scrutinize the solution in detail.Perturbation techniques depend on the existence of small or large parameters in the nonlinear problems.
The homotopy analysis method HAM , initially proposed by Liao in his Ph.D. thesis 10 , is a powerful analytic method for nonlinear problems.A systematic and clear exposition on HAM is given in 11 .In recent years, this method has been successfully employed to solve many types of nonlinear problems in science and engineering 12-29 .More recently, Bataineh et al. 30-33 employed the standard HAM to solve some problems in engineering sciences.HAM yields rapidly convergent series solutions in most cases, usually only a few iterations leading to very accurate solutions.Very recently, Bataineh et al. 34,35 presented two modifications of HAM MHAM to solve systems of second-order BVPs and homogeneous or nonhomogeneous differential equations with constant or variable coefficients.HAM and its modifications contain a certain auxiliary parameter , which provides us with a simple way to adjust and control the convergence region and rate of convergence of the series solution.Moreover, by means of the so-called -curve, it is easy to find the valid regions of to gain a convergent series solution.Thus, through HAM, explicit analytic solutions of nonlinear problems are possible.Apart from providing us with a functional form of the solution to the nonlinear problem, another advantage of HAM is that the method is valid for equations without small or large parameters like the Lotka-Volterra equation.
In this paper, we will employ HAM to obtain series solutions to the multispecies Lotka-Volterra competition models which are governed by a system of nonlinear ordinary differential equations.The HAM gives continuous solution which is of comparable accuracy to purely numerical method like the classical fourth-order Runge-Kutta method RK4 .The convergence theorem for the three-dimensional case is also given.

HAM for system of ODEs
We consider the following system of differential equations: where N i are nonlinear operators, t denotes the independent variable, y i t are unknown functions, and g i t are known analytic functions representing the nonhomogeneous terms.If g i t 0, 2.1 reduces to the homogeneous equation.By means of generalizing the traditional homotopy method 11 , we construct the so-called zeroth-order deformation equation: where q ∈ 0, 1 is an embedding parameter, is a nonzero auxiliary function, L is an auxiliary linear operator, y i,0 t are the initial guesses of y i t , and φ i t; q are unknown functions.It is important to note that one has great freedom to choose the auxiliary objects such as and L in HAM.Obviously, when q 0 and q 1, both hold.Thus, as q increases from 0 to 1, the solutions φ i t; q vary from the initial guesses y i,0 t to the solutions y i t .Expanding φ i t; q in Taylor series with respect to q, one has where If the auxiliary linear operator, the initial guesses, the auxiliary parameters , and the auxiliary functions are so properly chosen, then the series 2.4 converges at q 1 and which must be one of the solutions of the original nonlinear equation, as proved by 11 .As −1, 2.2 becomes 1 − q L φ i t; q − y i,0 t q N i φ i t; q − g i t 0, 2.7 which is used mostly in the HPM 36 .According to 2.5 , the governing equations can be deduced from the zeroth-order deformation equations 2.2 .Define the vectors y i,n y i,0 t , y i,1 t , . . ., y i,n t .

2.8
Differentiating 2.2 m times with respect to the embedding parameter q, then setting q 0, and finally dividing them by m!, we have the so-called mth-order deformation equation where

2.11
It should be emphasized that y i,m t m ≥ 1 is governed by the linear equation 2.9 with the linear initial/boundary conditions that come from the original problem, which can be easily solved by symbolic computation softwares such as Maple and Mathematica.

Applications of HAM
We will next obtain series solutions to the one-, two-, and three-dimensional cases of 1.2 -1.3 by HAM separately.We assume, in this paper, that the solution to 1.2 -1.3 can be expressed by the set of polynomial base functions: where A m are coefficients to be determined.This provides us with the so-called rule of solution expression; that is, the solution of 1.2 must be expressed in the same form as 3.2 .

One-dimensional case
Consider the 1D case of 1.2 , known as the Verhulst equation, where a and b are positive constants.The exact solution of 3.3 can be found by direct integration and is given by y t be bt b ay t /y 0 − ae bt for b / 0.

3.4
For definiteness, we will assume the following initial condition: To solve 3.3 by HAM with the initial condition 3.5 and with b 1 and a 3, we first choose the initial approximation and the linear operator L φ t; q ∂φ t; q ∂t , 3.7 with the property where c is an integral constant.Furthermore, 3.3 suggests that we define the nonlinear operator as N φ t; q ∂φ t; q ∂t − bφ t; q aφ 2 t; q .3.9 Using the above definition, we construct the zeroth-order deformation equation as in 2.2 , and the mth-order deformation equation for m ≥ 1 is as in 2.9 , subject to the initial condition where Now, the solution of the mth-order deformation equation 2.9 becomes where the integration of constant c is determined by the initial condition 3.10 .We now successively obtain . . .

3.13
In general, the analytic solution of 1.2 via the polynomial base functions is given by 3.14

3.20
Again, using the above definition, we construct the zeroth-order deformation equation as in 2.2 , and the mth-order deformation equation for m ≥ 1 is as in 2.9 , subject to the initial condition where

3.22
where the prime denotes differentiation with respect to the similarity variable t.Now, the solution of the mth-order deformation equation 2.9 for m ≥ 1 and i 1, 2 is given by and the auxiliary linear operators are as in 3.18 with the property 3.19 , where c i i 1, 2, 3 are constants of integrations.In a similar way as in the previous systems, we obtain the mthorder deformation equation 2.9 , where

3.29
subject to the initial condition y i,m 0 0.

3.30
Following similar procedure as in the previous section, we find the first two terms of the series solution:

3.31
Then, the solution expression via the polynomial base functions can be written as in 3.25 .
We will next give a convergence theorem for the HAM series solution of the 3D version of 1.2 .

Convergence theorem
As long as the series y i t y i,0 t ∞ m 1 y i,m t converges, where y i,m t is governed by 2.9 under the definitions 3.29 , 3.30 , and 2.11 , it must be the solution of 3.26 .
Proof.If the series is convergent, we can write for i 1, 2, 3 that From 2.9 and by using the definitions 2.11 and 3.18 , we then have On the other hand, substituting 3.29 , respectively, into the above expressions and simplifying them, we obtain y 3,j t S 3 t − S 3 t b 3 a 31 S 1 t a 32 S 2 t a 33 S 3 t 0.

3.36
From the initial conditions 3.28 and 3.30 , there hold y 3,m 0 0.5.

3.37
So, S i t satisfy 3.26 , and are therefore solutions of the 3D version of 1.2 with the initial condition 3.30 .This ends the proof.

Results and discussions
The series solutions of 1.2 given by HAM contain the auxiliary parameter .The validity of the method is based on such an assumption that the series 2.4 converges at q 1.It is the auxiliary parameter which ensures that this assumption can be satisfied.In general, by means of the so-called -curve, it is straightforward to choose a proper value of which ensures that the solution series is convergent.Figure 1 shows the -curves of 1D, 2D, and 3D obtained from the 20th-order HAM approximation solutions of 1.2 .From these figures, the valid regions of correspond to the line segments nearly parallel to the horizontal axis.In Figures 2, 3, and 4, it is demonstrated that the HAM solutions, taking −0.25, −0.8, and −0.6, agree very well with the solutions obtained by the classical fourth-order Runge-Kutta method at the step size Δt 0.001.

Conclusions
In this paper, the homotopy analysis method HAM was applied to solve the Lotka-Volterra equations.Polynomial base functions were found to give very good accuracy of HAM solutions for the Lotka-Volterra equations.The HAM gives continuous solution which is of comparable accuracy to purely numerical method like the classical fourth-order Runge-Kutta method RK4 .This is convenient for practical applications with minimum requirements on calculation and computation.The convergence theorem for the three-dimensional case is also given.We remark that the validity of the HAM series solutions can be enhanced by finding more terms and/or using the Padé technique.The functional form of the solution would be useful in the study of the stability of the system.
1 and Volterra 2 .The classical two-dimensional Lotka-Volterra equation is given by b 1 a 11 y 1 t a 12 y 2 t , dy 2 t dt y 2 t b 2 a 21 y 1 t a 22 y 2 t , 1 a 11 y 1 t a 12 y 2 t a 13 y 3 t , dy 2 t dt y 2 t b 2 a 21 y 1 t a 22 y 2 t a 23 y 3 t , dy 3 t dt y 3 t b 3 a 31 y 1 t a 32 y 2 t a 33 y 3 t ,