# Optimal Parametric Iteration Method for Solving Multispecies Lotka-Volterra Equations

Authors: Vasile Marinca; Nicolae Herişanu
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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 well-known 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 Heriş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].


## 2. 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:(2.1)ẋ1=f1(x1,x2,x3,…,xs)ẋ2=f2(x1,x2,x3,…,xs)⋮ẋs=fs(x1,x2,x3,…,xs)x1(0)=a1,x2(0)=a2,…,xs(0)=as,
where dot denotes derivative with respect to time t, fi,(i=1,2,…,s) are nonlinear functions.In the frame of OPIM we define the following iteration formulas:(2.2)ẋ1,n+1=f1(x1,n,x2,n,…,xs,n)+m1(t,Ci)f1,x1(x1,n,x2,n,xs,n)+n1(t,Ci)f1,x2(x1,n,x2,n,…,xs,n)+⋯+p1(t,Ci)f1,xs(x1,n,x2,n,…,xs,n)ẋ2,n+1=f2(x1,n,x2,n,…,xs,n)+m2(t,Cj)f2,x1(x1,n,x2,n,xs,n)+n2(t,Cj)f2,x2(x1,n,x2,n,…,xs,n)+⋯+p2(t,Cj)f2,xs(x1,n,x2,n,…,xs,n)⋮ẋs,n+1=fs(x1,n,x2,n,…,xs,n)+ms(t,Ck)fs,x1(x1,n,x2,n,xs,n)+ns(t,Ci)fs,x2(x1,n,x2,n,…,xs,n)+⋯+ps(t,Ck)fs,xs(x1,n,x2,n,…,xs,n),
where f1,xi=∂f1/∂xi, n=0,1,2… and Ci,Cj⋯Ck are unknown constants at this moment.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]:(2.3)ẋ=x(1-x-αy-βz),ẏ=y(1-βx-y-αz),ż=z(1-αx-βy-z),
where α and β are known parameters not necessary small. Initial conditions are(2.4)x(0)=a,y(0)=b,z(0)=c.In this special case, it is clear thats=3 and(2.5)x1=x,x2=y,x3=z,x4=x5=⋯=0f1=x(1-x-αy-βz),f2=y(1-βx-y-αz),f3=z(1-αx-βy-z),f4=f5=⋯=0.With the notations (2.5), the iteration formulas (2.2) can be written in the form(2.6)ẋn+1=f(xn,yn,zn)+m1(t,Ci)fx(xn,yn,zn)+n1(t,Ci)fy(xn,yn,zn)+p1(t,Ci)fz(xn,yn,zn)ẏn+1=g(xn,yn,zn)+m2(t,Cj)gx(xn,yn,zn)+n2(t,Cj)gy(xn,yn,zn)+p2(t,Cj)gz(xn,yn,zn)żn+1=h(xn,yn,zn)+m3(t,Ck)hx(xn,yn,zn)+n3(t,Ck)hy(xn,yn,zn)+p3(t,Ck)hz(xn,yn,zn),
where(2.7)f(x,y,z)=x(1-x-αy-βz),g(x,y,z)=y(1-βx-y-αz),h(x,y,z)=z(1-αx-βy-z),
and fx=∂f/∂x, n=0,1,2…, and Ci,Cj,Ck, i,j,k=1,2,… are unknown constants at this moment. There are many possibilities to choose the auxiliary functions me(t,Ci), ne(t,Cj) and pe(t,Ck), e=1,2,3. Basically, the shape of me, ne, and pe must follow the terms appearing in (2.6). Therefore, we try to choose the auxiliary functions so that in (2.6) the products meFx are of the same shape with the function Fx, where F∈{f,g,h}. The constants Ci,Cj,Ck which appear in the auxiliary functions me, ne, and pe can be optimally determined via various methods, for example, putting the conditions that the residual functionals J1, J2, and J3 given, respectively, by(2.8)J1=∫ab[ẋ-f(x,y,z)]2dt,J2=∫ab[ẏ-g(x,y,z)]2dt,J3=∫ab[ż-h(x,y,z)]2dt
be minimum, that is:(2.9)∂Ji∂Cj=0,i=1,2,3,j=1,2,…,s
taking into account also the initial conditions. Other methods suitable to determine the constants Cj are the Galerkin method, the collocation method, and so on.On the other hand, the initial approximationsx0, y0, and z0 and also the auxiliary functions me, ne, and pe 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 convergence-control constants C1,C2,… which lead to an excellent agreement of the approximate solutions with the exact ones.


## 3. 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(3.1)x0=ae-t,y0=be-t,z0=ce-t.Forn=0 into (2.6) we can construct the following iteration formulas:(3.2)ẋ1=x0-x02-αx0y0-βx0z0+m1(t,Ci)[1-(2x0+αy0+βz0)]+n1(t,Ci)(-αx0)+p1(t,Ci)(-βx0),ẏ1=y0-βx0y0-y02-αy0z0+m2(t,Cj)(-βy0)+n2(t,Cj)[1-(βx0+2y0+αz0)]+p2(t,Ci)(-αy0),ż1=z0-αx0z0-βy0z0-z02+m3(t,Ck)(-αz0)+n3(t,Ck)(-βz0)+p3(t,Ck)×[1-(αx0+βy0+2z0)].For the auxiliary functionsmi, ni, pi, i=1,2,3 we choose the expressions(3.3)m1(t,Ci)=(C1t+C2)e-t,n1(t,Ci)=C3t+C4,p1(t,Ci)=(C5t+C6)e-tm2(t,Cj)=(C7t+C8),n2(t,Cj)=(C9t+C10)e-t,p2(t,Cj)=(C11t+C12)e-2tm3(t,Ck)=(C13t+C14)e-t,n3(t,Ck)=(C15t+C16)e-2t,p3(t,Ck)=(C17t+C18)e-t.Alternatively, we can consider either other expressions for the initial approximations and the auxiliary function, follows as:(3.4)x0=ae-2t,y0=b,z0=ce-t,m1(t,Ci′)=C1′e-2t,n1(t,Ci′)=(C2′t+C3′)e-t,p1(t,Ci′)=(C4′t+C5)e-2tm2(t,Cj′)=(C6′t+C7′)e-t,n2(t,Cj′)=C8′e-2t,p2(t,Cj′)=(C9′t+C10′)e-tm3(t,Ck′)=C11′e-t,n3(t,Ck′)=C12′t+C13′,p3(t,Ck′)=(C14′t+C15′)e-2t,
and so on.Now, substituting (3.1) and (3.3) into (3.2), we obtain the equations in ẋ1, ẏ1, and ż1, respectively,(3.5)ẋ1=(-αaC3t+a-αaC4)e-t+[(C1-βaC5)t+C2-βaC6-(a2+αab+βac)]e-2t-(2a+αb+βc)(C1t+C2)e-3tẏ1=[(C9t-βbC7)t+b-βbC8+C10]e-t-[(βa+2b+αc)C9t+βab+b2+αbc+(βa+2b+αc)C10]e-2t-(αbC11t+αbC12)e-3tż1=(C17t+C18+c)e-t-{[(αa+βc+2c)C17+αcC13]t+αac+βbc+c2+αcC14+(αa+βc+2c)C18}e-2t-(βcC14t+βcC16)e-3t.The solutions of (3.5) with the initial conditions(3.6)x1(0)=a,y1(0)=b,z1(0)=c
can be written in the form(3.7)x1(t)=a+αaC3te-t+(a-αaC3-αaC4)(1-e-t)-12(C1-βaC5)te-2t+14[C1-βaC5+2(C2-βaC6)-2(a2+αab+βac)](1-e-2t)+13(2a+αb+βc)×C1te-3t-19(2a+αb+βc)(C1+3C2)(1-e-3t)y1(t)=b+(b-βbC8+C10)(1-e-t)+(βbC7-C9)te-t-14[(βa+2b+αc)C9+2(βab+b2+αbc)+2(βa+2b+αc)C10](1-e-2t)+12(βa+2b+αc)C9te-2t-19αb(C11+3C12)(1-e-3t)+13αbC11te-3tz1(t)=c+(C17+C18+c)(1-e-t)-C17te-t-14[2(αac+βbc+c2)+2αcC14+2(αa+βc+2c)C18+αcC1+(αa+βc+2c)C17]×(1-e-2t)+12[αcC13+(αa+βc+2c)C17]te-2t+13βcC15te-3t-14βc(C15+3C16)×(1-e-3t).
By substituting (3.7) into (2.3) it results in the residuals:(3.8)R1(t,Ci)=ẋ1-x1(1-x1-αy1-βz1),R2(t,Ci)=ẏ1-y1(1-βx1-y1-αz1),R3(t,Ci)=ż1-z1(1-αx1-βy1-z1).Making collocation in the arbitrary pointsti,i=1,2,…,6(3.9)R1(t1,Ci)=R1(t2,Ci)=⋯=R1(t6,Ci)=0R2(t1,Ci)=R2(t2,Ci)=⋯=R2(t6,Ci)=0R3(t1,Ci)=R3(t2,Ci)=⋯=R3(t6,Ci)=0,
we obtain the optimal values of the constants C1,C2,…,C18 and therefore the solution (3.7) in the first approximation is well determined.In the case whena=0.2, b=0.3, c=0.5, α=0.1, β=0.1, t1=0.3, t2=0.7, t3=1, t4=2, t5=3, t6=5, from (3.9) we obtain(3.10)C1=1.669841466,C2=-0.452081708,C3=0.899181851,C4=2.361528973,C5=-1.940901217,C6=0.582517076,C7=1033.008542,C8=-540.7564632,C9=30.90993616,C10=-14.59841284,C11=29.8078082,C12=94.01456893,C14=53.21346396,C15=-22.21139178,C16=-35.8538265,C17=-0.345181562,C18=0.993490199.By substituting (3.10) into (3.7) we obtain after only one iteration, the following expression for the first-order approximate solution of Lotka-Volterra equations in three species (2.3):(3.11)x1(t)=0.8281811-(0.452081708t+1.669841466)e-t+(2.361528973t-0.899181851)e-2t+(0.582517076t+1.940902217)e-3t,(3.12)y1(t)=0.83333001+(0.08032009t-1.924281055)e-t+(1.715501457t+0.342446025)e-2t+(0.298078082t+1.04850505)e-3t,(3.13)z1(t)=0.836102311+(0.345181562t-1.148308637)e-t+(0.001606241t+1.532769033)e-2t-(0.370189863t+0.720960396)e-3t.It is easy to verify the accuracy of the obtained solution if we graphically compare the analytical solution with the numerical one. Figures1, 2, and 3 show the comparison between the present solutions and the numerical integration results obtained by a fourth-order Runge-Kutta method.Figure 1
Comparison between the approximate solution (3.11) and numerical results obtained for (2.3) in case a=0.2,b=0.3,c=0.5,α=0.1,β=0.1: — numerical solution; - - - approximate solution.Figure 2
Comparison between the approximate solution (3.12) and numerical results obtained for (2.3) in case a=0.2, b=0.3, c=0.5, α=0.1, β=0.1: — numerical solution; - - - approximate solution.Figure 3
Comparison between the approximate solution (3.13) and numerical results obtained for (2.3) in case a=0.2, b=0.3, c=0.5, α=0.1, β=0.1: — numerical solution; - - - approximate solution.It can be seen from Figures1–3 that the solution obtained by OPIM is in very good agreement with numerical integration results.


## 4. 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 functionsme, ne, and pe, e=1,2,3 whose parameters C1,C2,… (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 C1,C2,…. 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].



---

*Source: 842121-2012-02-14.xml*