# Series Solution of the Multispecies  Lotka-Volterra Equations by Means

Authors: A. Sami Bataineh; M. S. M. Noorani; I. Hashim
Journal: Differential Equations and Nonlinear Mechanics
 (2008)

---

## Abstract

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.


---

## Body

## 1. Introduction

The simplest  model of
predator-prey interactions was developed independently by Lotka [1] and Volterra [2]. The classical two-dimensional Lotka-Volterra equation
is given by(1.1)dy1dt=ay1−by1y2,dy2dt=−cy2+by1y2.System (1.1) has been one
of the most studied models for a two-dimensional dynamical system. The
generalized n-dimensional Lotka-Volterra equations are
given by (cf. [3, 4])(1.2)dyi(t)dt=yi(t)(bi+∑j=1naijyj(t)),i=1,2,…,n,subject to the initial conditions(1.3)yi(0)=ci,i=1,2,…,n,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.


## 2. HAM for System of ODEs

We consider the following system of differential
equations:(2.1)Ni[yi(t)]=gi(t),i=1,2,…,n,where Ni are nonlinear operators, t denotes the independent variable, yi(t) are unknown functions, and gi(t) are known analytic functions representing the nonhomogeneous terms. If gi(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: (2.2)(1−q)L[ϕi(t;q)−yi,0(t)]=qℏ{Ni[ϕi(t;q)]−gi(t)},where q∈[0,1] is an embedding parameter, ℏ is a nonzero auxiliary function, L is an auxiliary linear operator, yi,0(t) are the initial guesses of yi(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(2.3)ϕi(t;0)=yi,0(t),ϕi(t;1)=yi(t)hold. Thus, as q increases from 0 to 1, the solutions ϕi(t;q) vary from the initial guesses yi,0(t) to the solutions yi(t). Expanding ϕi(t;q) in Taylor series with respect to q,
one has(2.4)ϕi(t;q)=yi,0(t)+∑m=1+∞yi,m(t)qm,where(2.5)yi,m=1m!∂mϕi(t;q)∂qm|q=0.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(2.6)ϕi(t;1)=yi(t)=yi,0(t)+∑m=1+∞yi,m(t),which must be one of the
solutions of the original nonlinear equation, as proved by [11]. As ℏ=−1,
(2.2) becomes(2.7)(1−q)L[ϕi(t;q)−yi,0(t)]+q{Ni[ϕi(t;q)]−gi(t)}=0,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(2.8)y→i,n={yi,0(t),yi,1(t),…,yi,n(t)}.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(2.9)L[yi,m(t)−χmyi,m-1(t)]=ℏRi,m(y→i,m-1),where(2.10)Ri,m(y→i,m-1)=1(m−1)!∂m−1{Ni[ϕi(t;q)]−gi(t)}∂qm−1|q=0,(2.11)χm={0,m≤1,1,m>1.It should be emphasized that yi,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.


## 3. 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:(3.1){tm∣m=1,2,3,…},from which we
have(3.2)y(t)=∑m=0∞Amtm,where Am 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).


### 3.1. One-Dimensional Case

Consider the 1D case of
(1.2), known as the Verhulst equation,(3.3)dydt=by−ay2,where a and b are positive constants. The exact solution of
(3.3) can be found by direct integration and is given by(3.4)y(t)=bebt(b+ay(t))/y(0)−aebtforb≠0.For definiteness, we will assume
the following initial condition:(3.5)y(0)=0.1.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(3.6)y0(t)=0.1and the linear
operator(3.7)L[ϕ(t;q)]=∂ϕ(t;q)∂t,with the
property(3.8)L[c]=0,where c is an integral constant. Furthermore, (3.3)
suggests that we define the nonlinear operator as(3.9)N[ϕ(t;q)]=∂ϕ(t;q)∂t−bϕ(t;q)+aϕ2(t;q).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(3.10)ym(0)=0,where(3.11)Rm[y→m-1]=ym-1′(t)−bym-1(t)+a∑i=0m−1yi(t)ym-1-i(t).Now, the solution of the mth-order
deformation equation (2.9) becomes(3.12)ym(t)=χmym-1(t)+ℏ∫0tRm(y→m-1)dτ+c,where the integration of
constant c is determined by the initial condition (3.10).
We now successively obtain(3.13)y1(t)=−7100ℏt,y2(t)=−7100ℏt−7100ℏ2t+141000ℏ2t2⋮In general, the
analytic solution of (1.2) via the polynomial base functions is given
by(3.14)y(t)=∑m=1+∞dm(ℏ)tm.


### 3.2. Two-Dimensional Case

Now we apply HAM to solve the 2D version of
(1.2):(3.15)dy1(t)dt=y1(t)[b1+a11y1(t)+a12y2(t)],dy2(t)dt=y2(t)[b2+a21y1(t)+a22y2(t)],where a's and b's are constants and subject to the initial
conditions(3.16)y1(0)=4,y2(0)=10.According to HAM, the initial
approximations of system (3.15)-(3.16)
are(3.17)y1,0(t)=4,y2,0(0)=10,and the auxiliary linear
operators for i=1,2 are(3.18)L[ϕi(t;q)]=∂ϕi(t;q)∂t,i=1,2,with the
property(3.19)L[ci]=0and the nonlinear
operators(3.20)N1[ϕi(t;q)]=∂ϕ1(t;q)∂t−ϕ1(t;q)(b1+a11ϕ1(t;q)+a12ϕ2(t;q)),N2[ϕi(t;q)]=∂ϕ2(t;q)∂t−ϕ2(t,q)(b2+a21ϕ1(t;q)+a22ϕ2(t;q)).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(3.21)yi,m(0)=0,where(3.22)R1,m[y→i,m-1(t)]=y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y1,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t),R2,m[y→i,m-1(t)]=y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t),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(3.23)yi,m(t)=χmyi,m-1(t)+ℏ∫0tR1,m(y→i,m-1)dτ+ci,where the integration constants ci (i=1,2) are determined by the initial
condition (3.21). Thereafter, we successively
obtain(3.24)y1,1(t)=−0.3296ℏt,y1,2(t)=−0.3296ℏt−0.3296ℏ2t+0.011063ℏ2t2,y2,1(t)=−0.664ℏt,y2,2(t)=−0.664ℏt−0.664ℏ2t+0.0172416ℏ2t2and so forth. Thus, the analytic
solution of (3.15)-(3.16) has the general form(3.25)yi(t)=∑m=1+∞ai,m(ℏ)tm.


### 3.3. Three-Dimensional Case

Finally, we apply HAM to solve the 3D version of
(1.2):(3.26)dy1(t)dt=y1(t)(b1+a11y1(t)+a12y2(t)+a13y3(t)),dy2(t)dt=y2(t)(b2+a21y1(t)+a22y2(t)+a23y3(t)),dy3(t)dt=y3(t)(b3+a31y1(t)+a32y2(t)+a33y3(t)), with the initial
conditions(3.27)y1,0(0)=0.2,y2,0(0)=0.3,y3,0(0)=0.5,where the a's and b's are constants.According to HAM, the initial approximations of
(3.26) are(3.28)y1,0(t)=0.2,y2,0(t)=0.3,y3,0(t)=0.5,and the auxiliary linear
operators are as in (3.18) with the property (3.19), where ci (i=1,2,3) are constants of integrations. In a similar
way as in the previous systems, we obtain the mth-order deformation
equation (2.9), where(3.29)R1,m[y→i,m-1(t)]=y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y1,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t)−a13∑i=0m−1y1,i(t)y3,m-1-i(t),R2,m[y→i,m-1(t)]=y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t)−a23∑i=0m−1y2,i(t)y3,m-1-i(t),R3,m[y→i,m-1(t)]=y3,m-1′(t)−b3y3,m-1(t)−a31∑i=0m−1y3,i(t)y1,m-1-i(t)−a32∑i=0m−1y3,i(t)y2,m-1-i(t)−a33∑i=0m−1y3,i(t)y3,m-1-i(t), subject to the initial
condition(3.30)yi,m(0)=0.Following similar procedure as
in the previous section, we find the first two terms of the series solution:(3.31)y1,1(t)=−0.144ℏt,y1,2(t)=−0.144ℏt−0.144ℏ2t+0.0333ℏ2t2,y2,1(t)=−0.189ℏt,y2,2(t)=−0.189ℏt−0.189ℏ2t+0.02565ℏ2t2,y3,1(t)=−0.225ℏt,y3,2(t)=−0.225ℏt−0.225ℏ2t−0.01395ℏ2t2.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).


#### 3.3.1. Convergence Theorem

As long  as the seriesyi(t)=yi,0(t)+∑m=1+∞yi,m(t) converges, where yi,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 fori=1,2,3 that(3.32)Si=∑m=0+∞yi,m(t),and there hold(3.33)limn→+∞y1,n(t)=0,limn→+∞y2,n=0,limn→+∞y3,n=0.From (2.9) and by using the
definitions (2.11) and (3.18), we then have(3.34)ℏ∑m=1+∞Ri,m(y→i,m-1)=∑m=1+∞L[yi,m(t)−χmyi,m-1(t)]=limn→+∞∑m=1nL[yi,m(t)−χmyi,m-1(t)]=L[limn→+∞∑m=1nyi,m(t)−χmyi,m-1(t)]=L[limn→+∞yi,n(t)]=0,which gives, since ℏ≠0,(3.35)∑n=1+∞Ri,n(y→i,n-1)=0.On the other hand, substituting
(3.29), respectively, into the above expressions and simplifying
them, we obtain(3.36)∑m=1+∞R1,m(y→i,m-1)=∑m=1+∞[y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y2,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t)−a13∑i=0m−1y1,i(t)y3,m-1-i(t)]=∑m=1+∞y1,m−1′(t)−b1y1,m−1(t)−a11∑m=1+∞∑i=0m−1y1,i(t)y1,m−1−i(t)−a12∑m=1+∞∑i=0m−1y1,i(t)y2,m−1−i(t)−a13∑m=1+∞∑i=0m−1y1,i(t)y3,m−1−i(t)=∑m=1+∞y1,m−1′(t)−b1y1,m−1(t)−a11∑i=0+∞∑m=i+1+∞y1,i(t)y1,m−1−i(t)−a12∑i=0+∞∑m=i+1+∞y1,i(t)y2,m−1−i(t)−a13∑i=0+∞∑m=i+1+∞y1,i(t)y3,m−1−i(t)=∑m=0+∞y1,m′(t)−b1y1,m(t)−a11∑i=0+∞y1,i(t)∑j=0+∞y1,j(t)−a12∑i=0+∞y1,i(t)∑j=0+∞y2,j(t)−a13∑i=0+∞y1,i(t)∑j=0+∞y3,j(t)=S1′(t)−S1(t)[b1+a11S1(t)+a12S2(t)+a13S3(t)]=0,∑m=1+∞R2,m(y→i,m-1)=∑m=1+∞[y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t)−a23∑i=0m−1y2,i(t)y3,m-1-i(t)]=∑m=1+∞y2,m−1′(t)−b2y2,m−1(t)−a21∑m=1+∞∑i=0m−1y2,i(t)y1,m−1−i(t)−a22∑m=1+∞∑i=0m−1y2,i(t)y2,m−1−i(t)−a23∑m=1+∞∑i=0m−1y2,i(t)y3,m−1−i(t)=∑m=1+∞y2,m−1′(t)−b2y2,m−1(t)−a21∑i=0+∞∑m=i+1+∞y2,i(t)y1,m−1−i(t)−a22∑i=0+∞∑m=i+1+∞y2,i(t)y2,m−1−i(t)−a23∑i=0+∞∑m=i+1+∞y2,i(t)y3,m−1−i(t)=∑m=0+∞y2,m′(t)−b2y2,m(t)−a21∑i=0+∞y2,i(t)∑j=0+∞y1,j(t)−a22∑i=0+∞y2,i(t)∑j=0+∞y2,j(t)−a23∑i=0+∞y2,i(t)∑j=0+∞y3,j(t)=S2′(t)−S2(t)[b2+a21S1(t)+a12S2(t)+a13S3(t)]=0,∑m=1+∞R3,m(y→i,m-1)=∑m=1+∞[y3,m-1′(t)−b3y3,m−1(t)−a31∑i=0m−1y3,i(t)y1,m−1−i(t)−a32∑i=0m−1y3,i(t)y2,m−1−i(t)−a33∑i=0m−1y3,i(t)y3,m−1−i(t)]=∑m=1+∞y3,m-1′(t)−b3y3,m−1(t)−a31∑m=1+∞∑i=0m−1y3,i(t)y1,m−1−i(t)−a32∑m=1+∞∑i=0m−1y3,i(t)y2,m−1−i(t)−a33∑m=1+∞∑i=0m−1y3,i(t)y3,m−1−i(t)=∑m=1+∞y3,m-1′(t)−b3y3,m−1(t)−a31∑i=0+∞∑m=i+1+∞y3,i(t)y1,m−1−i(t)−a32∑i=0+∞∑m=i+1+∞y3,i(t)y2,m−1−i(t)−a33∑i=0+∞∑m=i+1+∞y3,i(t)y3,m−1−i(t)=∑m=0+∞y3,m′(t)−b3y3,m(t)−a31∑i=0+∞y3,i(t)∑j=0+∞y1,j(t)−a32∑i=0+∞y3,i(t)∑j=0+∞y2,j(t)−a33∑i=0+∞y3,i(t)∑j=0+∞y3,j(t)=S3′(t)−S3(t)[b3+a31S1(t)+a32S2(t)+a33S3(t)]=0.From the initial conditions (3.28)
and (3.30), there hold(3.37)S1(0)=y1,0(0)+∑m=1+∞y1,m(0)=0.2,S2(0)=y2,0(0)+∑m=1+∞y2,m(0)=0.3,S3(0)=y3,0(0)+∑m=1+∞y3,m(0)=0.5.So, Si(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.


## 3.1. One-Dimensional Case

Consider the 1D case of
(1.2), known as the Verhulst equation,(3.3)dydt=by−ay2,where a and b are positive constants. The exact solution of
(3.3) can be found by direct integration and is given by(3.4)y(t)=bebt(b+ay(t))/y(0)−aebtforb≠0.For definiteness, we will assume
the following initial condition:(3.5)y(0)=0.1.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(3.6)y0(t)=0.1and the linear
operator(3.7)L[ϕ(t;q)]=∂ϕ(t;q)∂t,with the
property(3.8)L[c]=0,where c is an integral constant. Furthermore, (3.3)
suggests that we define the nonlinear operator as(3.9)N[ϕ(t;q)]=∂ϕ(t;q)∂t−bϕ(t;q)+aϕ2(t;q).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(3.10)ym(0)=0,where(3.11)Rm[y→m-1]=ym-1′(t)−bym-1(t)+a∑i=0m−1yi(t)ym-1-i(t).Now, the solution of the mth-order
deformation equation (2.9) becomes(3.12)ym(t)=χmym-1(t)+ℏ∫0tRm(y→m-1)dτ+c,where the integration of
constant c is determined by the initial condition (3.10).
We now successively obtain(3.13)y1(t)=−7100ℏt,y2(t)=−7100ℏt−7100ℏ2t+141000ℏ2t2⋮In general, the
analytic solution of (1.2) via the polynomial base functions is given
by(3.14)y(t)=∑m=1+∞dm(ℏ)tm.


## 3.2. Two-Dimensional Case

Now we apply HAM to solve the 2D version of
(1.2):(3.15)dy1(t)dt=y1(t)[b1+a11y1(t)+a12y2(t)],dy2(t)dt=y2(t)[b2+a21y1(t)+a22y2(t)],where a's and b's are constants and subject to the initial
conditions(3.16)y1(0)=4,y2(0)=10.According to HAM, the initial
approximations of system (3.15)-(3.16)
are(3.17)y1,0(t)=4,y2,0(0)=10,and the auxiliary linear
operators for i=1,2 are(3.18)L[ϕi(t;q)]=∂ϕi(t;q)∂t,i=1,2,with the
property(3.19)L[ci]=0and the nonlinear
operators(3.20)N1[ϕi(t;q)]=∂ϕ1(t;q)∂t−ϕ1(t;q)(b1+a11ϕ1(t;q)+a12ϕ2(t;q)),N2[ϕi(t;q)]=∂ϕ2(t;q)∂t−ϕ2(t,q)(b2+a21ϕ1(t;q)+a22ϕ2(t;q)).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(3.21)yi,m(0)=0,where(3.22)R1,m[y→i,m-1(t)]=y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y1,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t),R2,m[y→i,m-1(t)]=y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t),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(3.23)yi,m(t)=χmyi,m-1(t)+ℏ∫0tR1,m(y→i,m-1)dτ+ci,where the integration constants ci (i=1,2) are determined by the initial
condition (3.21). Thereafter, we successively
obtain(3.24)y1,1(t)=−0.3296ℏt,y1,2(t)=−0.3296ℏt−0.3296ℏ2t+0.011063ℏ2t2,y2,1(t)=−0.664ℏt,y2,2(t)=−0.664ℏt−0.664ℏ2t+0.0172416ℏ2t2and so forth. Thus, the analytic
solution of (3.15)-(3.16) has the general form(3.25)yi(t)=∑m=1+∞ai,m(ℏ)tm.


## 3.3. Three-Dimensional Case

Finally, we apply HAM to solve the 3D version of
(1.2):(3.26)dy1(t)dt=y1(t)(b1+a11y1(t)+a12y2(t)+a13y3(t)),dy2(t)dt=y2(t)(b2+a21y1(t)+a22y2(t)+a23y3(t)),dy3(t)dt=y3(t)(b3+a31y1(t)+a32y2(t)+a33y3(t)), with the initial
conditions(3.27)y1,0(0)=0.2,y2,0(0)=0.3,y3,0(0)=0.5,where the a's and b's are constants.According to HAM, the initial approximations of
(3.26) are(3.28)y1,0(t)=0.2,y2,0(t)=0.3,y3,0(t)=0.5,and the auxiliary linear
operators are as in (3.18) with the property (3.19), where ci (i=1,2,3) are constants of integrations. In a similar
way as in the previous systems, we obtain the mth-order deformation
equation (2.9), where(3.29)R1,m[y→i,m-1(t)]=y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y1,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t)−a13∑i=0m−1y1,i(t)y3,m-1-i(t),R2,m[y→i,m-1(t)]=y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t)−a23∑i=0m−1y2,i(t)y3,m-1-i(t),R3,m[y→i,m-1(t)]=y3,m-1′(t)−b3y3,m-1(t)−a31∑i=0m−1y3,i(t)y1,m-1-i(t)−a32∑i=0m−1y3,i(t)y2,m-1-i(t)−a33∑i=0m−1y3,i(t)y3,m-1-i(t), subject to the initial
condition(3.30)yi,m(0)=0.Following similar procedure as
in the previous section, we find the first two terms of the series solution:(3.31)y1,1(t)=−0.144ℏt,y1,2(t)=−0.144ℏt−0.144ℏ2t+0.0333ℏ2t2,y2,1(t)=−0.189ℏt,y2,2(t)=−0.189ℏt−0.189ℏ2t+0.02565ℏ2t2,y3,1(t)=−0.225ℏt,y3,2(t)=−0.225ℏt−0.225ℏ2t−0.01395ℏ2t2.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).


### 3.3.1. Convergence Theorem

As long  as the seriesyi(t)=yi,0(t)+∑m=1+∞yi,m(t) converges, where yi,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 fori=1,2,3 that(3.32)Si=∑m=0+∞yi,m(t),and there hold(3.33)limn→+∞y1,n(t)=0,limn→+∞y2,n=0,limn→+∞y3,n=0.From (2.9) and by using the
definitions (2.11) and (3.18), we then have(3.34)ℏ∑m=1+∞Ri,m(y→i,m-1)=∑m=1+∞L[yi,m(t)−χmyi,m-1(t)]=limn→+∞∑m=1nL[yi,m(t)−χmyi,m-1(t)]=L[limn→+∞∑m=1nyi,m(t)−χmyi,m-1(t)]=L[limn→+∞yi,n(t)]=0,which gives, since ℏ≠0,(3.35)∑n=1+∞Ri,n(y→i,n-1)=0.On the other hand, substituting
(3.29), respectively, into the above expressions and simplifying
them, we obtain(3.36)∑m=1+∞R1,m(y→i,m-1)=∑m=1+∞[y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y2,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t)−a13∑i=0m−1y1,i(t)y3,m-1-i(t)]=∑m=1+∞y1,m−1′(t)−b1y1,m−1(t)−a11∑m=1+∞∑i=0m−1y1,i(t)y1,m−1−i(t)−a12∑m=1+∞∑i=0m−1y1,i(t)y2,m−1−i(t)−a13∑m=1+∞∑i=0m−1y1,i(t)y3,m−1−i(t)=∑m=1+∞y1,m−1′(t)−b1y1,m−1(t)−a11∑i=0+∞∑m=i+1+∞y1,i(t)y1,m−1−i(t)−a12∑i=0+∞∑m=i+1+∞y1,i(t)y2,m−1−i(t)−a13∑i=0+∞∑m=i+1+∞y1,i(t)y3,m−1−i(t)=∑m=0+∞y1,m′(t)−b1y1,m(t)−a11∑i=0+∞y1,i(t)∑j=0+∞y1,j(t)−a12∑i=0+∞y1,i(t)∑j=0+∞y2,j(t)−a13∑i=0+∞y1,i(t)∑j=0+∞y3,j(t)=S1′(t)−S1(t)[b1+a11S1(t)+a12S2(t)+a13S3(t)]=0,∑m=1+∞R2,m(y→i,m-1)=∑m=1+∞[y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t)−a23∑i=0m−1y2,i(t)y3,m-1-i(t)]=∑m=1+∞y2,m−1′(t)−b2y2,m−1(t)−a21∑m=1+∞∑i=0m−1y2,i(t)y1,m−1−i(t)−a22∑m=1+∞∑i=0m−1y2,i(t)y2,m−1−i(t)−a23∑m=1+∞∑i=0m−1y2,i(t)y3,m−1−i(t)=∑m=1+∞y2,m−1′(t)−b2y2,m−1(t)−a21∑i=0+∞∑m=i+1+∞y2,i(t)y1,m−1−i(t)−a22∑i=0+∞∑m=i+1+∞y2,i(t)y2,m−1−i(t)−a23∑i=0+∞∑m=i+1+∞y2,i(t)y3,m−1−i(t)=∑m=0+∞y2,m′(t)−b2y2,m(t)−a21∑i=0+∞y2,i(t)∑j=0+∞y1,j(t)−a22∑i=0+∞y2,i(t)∑j=0+∞y2,j(t)−a23∑i=0+∞y2,i(t)∑j=0+∞y3,j(t)=S2′(t)−S2(t)[b2+a21S1(t)+a12S2(t)+a13S3(t)]=0,∑m=1+∞R3,m(y→i,m-1)=∑m=1+∞[y3,m-1′(t)−b3y3,m−1(t)−a31∑i=0m−1y3,i(t)y1,m−1−i(t)−a32∑i=0m−1y3,i(t)y2,m−1−i(t)−a33∑i=0m−1y3,i(t)y3,m−1−i(t)]=∑m=1+∞y3,m-1′(t)−b3y3,m−1(t)−a31∑m=1+∞∑i=0m−1y3,i(t)y1,m−1−i(t)−a32∑m=1+∞∑i=0m−1y3,i(t)y2,m−1−i(t)−a33∑m=1+∞∑i=0m−1y3,i(t)y3,m−1−i(t)=∑m=1+∞y3,m-1′(t)−b3y3,m−1(t)−a31∑i=0+∞∑m=i+1+∞y3,i(t)y1,m−1−i(t)−a32∑i=0+∞∑m=i+1+∞y3,i(t)y2,m−1−i(t)−a33∑i=0+∞∑m=i+1+∞y3,i(t)y3,m−1−i(t)=∑m=0+∞y3,m′(t)−b3y3,m(t)−a31∑i=0+∞y3,i(t)∑j=0+∞y1,j(t)−a32∑i=0+∞y3,i(t)∑j=0+∞y2,j(t)−a33∑i=0+∞y3,i(t)∑j=0+∞y3,j(t)=S3′(t)−S3(t)[b3+a31S1(t)+a32S2(t)+a33S3(t)]=0.From the initial conditions (3.28)
and (3.30), there hold(3.37)S1(0)=y1,0(0)+∑m=1+∞y1,m(0)=0.2,S2(0)=y2,0(0)+∑m=1+∞y2,m(0)=0.3,S3(0)=y3,0(0)+∑m=1+∞y3,m(0)=0.5.So, Si(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.


## 3.3.1. Convergence Theorem

As long  as the seriesyi(t)=yi,0(t)+∑m=1+∞yi,m(t) converges, where yi,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 fori=1,2,3 that(3.32)Si=∑m=0+∞yi,m(t),and there hold(3.33)limn→+∞y1,n(t)=0,limn→+∞y2,n=0,limn→+∞y3,n=0.From (2.9) and by using the
definitions (2.11) and (3.18), we then have(3.34)ℏ∑m=1+∞Ri,m(y→i,m-1)=∑m=1+∞L[yi,m(t)−χmyi,m-1(t)]=limn→+∞∑m=1nL[yi,m(t)−χmyi,m-1(t)]=L[limn→+∞∑m=1nyi,m(t)−χmyi,m-1(t)]=L[limn→+∞yi,n(t)]=0,which gives, since ℏ≠0,(3.35)∑n=1+∞Ri,n(y→i,n-1)=0.On the other hand, substituting
(3.29), respectively, into the above expressions and simplifying
them, we obtain(3.36)∑m=1+∞R1,m(y→i,m-1)=∑m=1+∞[y1,m-1′(t)−b1y1,m-1(t)−a11∑i=0m−1y1,i(t)y2,m-1-i(t)−a12∑i=0m−1y1,i(t)y2,m-1-i(t)−a13∑i=0m−1y1,i(t)y3,m-1-i(t)]=∑m=1+∞y1,m−1′(t)−b1y1,m−1(t)−a11∑m=1+∞∑i=0m−1y1,i(t)y1,m−1−i(t)−a12∑m=1+∞∑i=0m−1y1,i(t)y2,m−1−i(t)−a13∑m=1+∞∑i=0m−1y1,i(t)y3,m−1−i(t)=∑m=1+∞y1,m−1′(t)−b1y1,m−1(t)−a11∑i=0+∞∑m=i+1+∞y1,i(t)y1,m−1−i(t)−a12∑i=0+∞∑m=i+1+∞y1,i(t)y2,m−1−i(t)−a13∑i=0+∞∑m=i+1+∞y1,i(t)y3,m−1−i(t)=∑m=0+∞y1,m′(t)−b1y1,m(t)−a11∑i=0+∞y1,i(t)∑j=0+∞y1,j(t)−a12∑i=0+∞y1,i(t)∑j=0+∞y2,j(t)−a13∑i=0+∞y1,i(t)∑j=0+∞y3,j(t)=S1′(t)−S1(t)[b1+a11S1(t)+a12S2(t)+a13S3(t)]=0,∑m=1+∞R2,m(y→i,m-1)=∑m=1+∞[y2,m-1′(t)−b2y2,m-1(t)−a21∑i=0m−1y2,i(t)y1,m-1-i(t)−a22∑i=0m−1y2,i(t)y2,m-1-i(t)−a23∑i=0m−1y2,i(t)y3,m-1-i(t)]=∑m=1+∞y2,m−1′(t)−b2y2,m−1(t)−a21∑m=1+∞∑i=0m−1y2,i(t)y1,m−1−i(t)−a22∑m=1+∞∑i=0m−1y2,i(t)y2,m−1−i(t)−a23∑m=1+∞∑i=0m−1y2,i(t)y3,m−1−i(t)=∑m=1+∞y2,m−1′(t)−b2y2,m−1(t)−a21∑i=0+∞∑m=i+1+∞y2,i(t)y1,m−1−i(t)−a22∑i=0+∞∑m=i+1+∞y2,i(t)y2,m−1−i(t)−a23∑i=0+∞∑m=i+1+∞y2,i(t)y3,m−1−i(t)=∑m=0+∞y2,m′(t)−b2y2,m(t)−a21∑i=0+∞y2,i(t)∑j=0+∞y1,j(t)−a22∑i=0+∞y2,i(t)∑j=0+∞y2,j(t)−a23∑i=0+∞y2,i(t)∑j=0+∞y3,j(t)=S2′(t)−S2(t)[b2+a21S1(t)+a12S2(t)+a13S3(t)]=0,∑m=1+∞R3,m(y→i,m-1)=∑m=1+∞[y3,m-1′(t)−b3y3,m−1(t)−a31∑i=0m−1y3,i(t)y1,m−1−i(t)−a32∑i=0m−1y3,i(t)y2,m−1−i(t)−a33∑i=0m−1y3,i(t)y3,m−1−i(t)]=∑m=1+∞y3,m-1′(t)−b3y3,m−1(t)−a31∑m=1+∞∑i=0m−1y3,i(t)y1,m−1−i(t)−a32∑m=1+∞∑i=0m−1y3,i(t)y2,m−1−i(t)−a33∑m=1+∞∑i=0m−1y3,i(t)y3,m−1−i(t)=∑m=1+∞y3,m-1′(t)−b3y3,m−1(t)−a31∑i=0+∞∑m=i+1+∞y3,i(t)y1,m−1−i(t)−a32∑i=0+∞∑m=i+1+∞y3,i(t)y2,m−1−i(t)−a33∑i=0+∞∑m=i+1+∞y3,i(t)y3,m−1−i(t)=∑m=0+∞y3,m′(t)−b3y3,m(t)−a31∑i=0+∞y3,i(t)∑j=0+∞y1,j(t)−a32∑i=0+∞y3,i(t)∑j=0+∞y2,j(t)−a33∑i=0+∞y3,i(t)∑j=0+∞y3,j(t)=S3′(t)−S3(t)[b3+a31S1(t)+a32S2(t)+a33S3(t)]=0.From the initial conditions (3.28)
and (3.30), there hold(3.37)S1(0)=y1,0(0)+∑m=1+∞y1,m(0)=0.2,S2(0)=y2,0(0)+∑m=1+∞y2,m(0)=0.3,S3(0)=y3,0(0)+∑m=1+∞y3,m(0)=0.5.So, Si(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.


## 4. 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.Theℏ-curves obtained from the 20th-order HAM
approximation solutions of (a) (3.3); (b) (3.15); (c) (3.26).
(a)(b)(c)Figure 2
The 61th-order HAM solution of (3.3) with ℏ=−0.25 versus RK4 (Δt=0.001) solution for the 1D case.Figure 3
The 31th-order HAM solution of (3.15)
with ℏ=−0.8 versus RK4 (Δt=0.001) solution for the 2D case.Figure 4
The 41th-order HAM solution of (3.26) with ℏ=−0.6 versus RK4 (Δt=0.001) solution for the 3D case.


## 5. 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.



---

*Source: 816787-2008-07-28.xml*