Anali PAZU - Letnik 2, leto 2012, številka 1 
 
22 
 
Competitive Predator-Prey Systems with Time-Dependent 
Coefficients: A Multistage Homotopy Perturbation Analysis 
Tekmovalni sistemi plenilec-plen s časovno odvisnimi 
koeficienti: Večstopenjska homotopsko perturbacijska analiza 
Rudolf Pušenjak
1,
*, Maks Oblak
1 
and Martin Lipičnik
1
 
1
University of Maribor, Faculty of Logistics / Mariborska cesta 7, SI-3000 Celje 
E-Mails: rudi.pusenjak@teleing.com; maks.oblak@fl.uni-mb.si; martin.lipicnik@fl.uni-mb.si 
* Avtor za korespondenco; Full Prof. Dr. Rudolf  Pušenjak, Kidričeva 5, 9240 Ljutomer, Slovenia; Tel.: +386-02-582-1-
406 
 
Abstract: This paper treats competitive predator-prey systems, which growth of populations and their mutual 
interactions are time dependent. For solving a general Lotka-Volterra system of equations, the multistage 
homotopy perturbation (MH-P) method  is developed to predict the time evolution of the dynamical system 
and its properties, such as existence of stable periodic orbits. As the newest achievement, the efficiency of 
MH-P method is provedintreatment of almost-periodic variations of coefficients with incommensurate 
excitation frequencies.The periodic variations of coefficients are analyzed as special case by assuming that 
excitation frequencies are commensurate. By using MH-P method, the approximate analytic solutions are 
obtained, which are very accurate in the long term behaviour. Although an usefull convergence test of the 
computed solution is provided, the accuracy of  MH-P method is compared also by results of the numerical 
integration of Lotka-Volterra equations by using the Runge-Kutta method, where an excellent agreement is 
obtained. 
Key words: predator-prey system; almost-periodic coefficients; MH-P method 
Povzetek: Članek obravnava tekmovalne sisteme plenilcev-plen, katerih rasti populacij in njih medsebojni 
vplivi so časovno odvisni. Za reševanje splošnega sistema enačb tipa Lotka-Volterra je razvita večstopenjska 
homotopsko perturbacijska (VH-P) metoda, s katero določimo časovni razvoj dinamičnega sistema in njegove 
lastnosti, kot denimo obstoj stabilnih periodičnih orbit. Kot najnovejši dosežek je prikazana učinkovitost VH-P 
metode v obravnavi skoraj periodičnih variacij koeficientov z nekomenzurnimi vzbujevalnimi frekvencami. 
Periodične variacije koeficientov so analizirane kot posebni primer s komenzurnimi vzbujevalnimi 
frekvencami. Z uporabo VH-P metode dobimo približne analitične rešitve, ki so zelo natančne v daljšem 
časovnem obdobju. Čeprav metoda predvideva uporaben test konvergence izračunanih rešitev, je natančnost 
VH-P  metode preverjena tudi s primerjavo rezultatov, dobljenih z numerično integracijo enačb Lotka-Volterra 
ob uporabi metode Runge-Kutta, pri čemer je ugotavljeno odlično ujemanje rezultatov. 
Ključne besede: sistemi plenilci-plen; skoraj periodični koeficienti; VH-P metoda 
 
1. Introduction 
Dynamical systems of predator-prey type appear in 
diverse branches of science and engineering. Predator-prey 
systems are applied successfully in the combustion theory 
to model the evolution of chemical radicals formed during 
H
2
 and O
2
 combustion [1] and are widely used in 
chemistry, biology, economics, medicine and ecology [2]. 
Due to the wide applicability such systems represent a 
great research challenge.The interaction mechanisms in 
predator-prey ecosystems are nonlinear because harvest-, 
birth-, death-, and migration rates depend strongly on 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
23 
 
populations. Unfortunately, such models does not take into 
account the exogenous environmental factors, such as 
food, temperature, humidity, light, etc., which vary 
periodically or even almost-periodically in dependence on 
the season-,  daily-, month- and year cycles. The objective 
of this paper is the inclusion of such exogenous factors into 
dynamical system to obtain a more realistic model of the 
predator-prey system. 
Dynamical system describing the predator-prey 
ecosystem, which takes into account the time-dependent 
coefficients of populations growth rates as well as the 
time-dependent coefficients of mutual interactions between 
populations of species, is written in the following form of 
generalized Lotka-Volterra equations: 
   
   
1
0
,
0 , 1, ,
n
i i i ij i j
j
i i i
x t x a t x x
x x c i n



 



  

 (1) 
Functions  
i
t  represent the time-dependent 
coefficients of populations growth rates and functions 
 
ij
at  stand for time-dependent coefficients of mutual 
interactions between populations of species.Because the 
right hand sides of ODE's are explicitly dependent on the 
time, the dynamical system, described by Eqs. (1) is 
nonautonomous. It is required, that all  
i
t  and  
ij
at , 
respectively, are continuous functions of time, which for 
example may be polinomials, periodic or even almost-
periodic functions. In this paper, almost-periodic functions 
will be treated in details, because they represent the most 
relevant basis for modelling the realistic predator-prey 
ecosystems. When all coefficients of mutual interactions 
between populations of speciesare nonnegative,   0
ij
at  
for all times and for all indexes i,j, the ecosystem is 
competitive and when   0
ij i j
at

 , then ecosystem 
iscooperative [3]. To restrict yourself in this paper, only 
competitive ecosystems will be treated, although the 
present formulation of the MH-P method is not limited on 
the competitive systems only. By means of equations
0ii
xc  , the corresponding initial conditions in the time 
0 t  are prescribed. 
Homotopy perturbation (H-P) method is videly applied 
in a broad range of nonlinear problems describing various 
physical phenomena [4].Such a phenomena are governed 
by systems of nonlinear ordinary or partial differential 
equations. The strength of the H-P method comes from the 
fact, that solutions can be expressed in the form of power 
series, where the coefficients of power series are 
determined in a similar manner as in perturbation methods 
without the need for the introduction of a small expansion 
parameter [5],[6].Due to this, the H-P method is ideally 
suited to solve strongly nonlinear problems if the 
nonlinearities appear in the polinomial form. Besides this, 
the power series produce analytical solutions, which offer 
numerous advantages (the implementation of parametric 
studies, stability and bifurcation analysis).On the other 
side, there is a limiting factor in the H-P method, namely 
the length of the time intervalin which is desired to obtain 
an accurate solution. In dynamical systems exhibiting 
chaos, the length of the time interval as a rule must be very 
large in order to obtain the fully developed characteristics 
of chaos. The standard H-P method applied on dynamical 
systems in which chaos appears, fails on large time 
intervals.To overcome this deficiency, the multistage 
homotopy perturbation (MH-P) method is first proposed in 
[7] and then applied in predator-prey systems in [8]. The 
main idea behind the MH-P method is the division of the 
entire time interval on the finite number of subintervals 
with equal length and succesive application of the H-P 
algorithm on each subinterval. The only difficulty in such a 
procedure are unknown initial conditions on all 
subintervals except on the first. This difficulty can be 
easily overcomed by assuming that unknown initial 
conditions must be equal to the final values on the 
preceding interval in order to ensure the C
0
 continuity. 
MH-P method can be applied purely natural on Lotka-
Volterra Eqs. (1), because these equations contain 
nonlinearities in the polinomial form. In this paper, the 
general form of the MH-P method is presented, which is 
capable to solve Lotka-Volterra equations with time-
dependent coefficients. The solution obtained has an 
explicit form. Application of the method is shown in 
details for almost-periodic coefficients with 
incommensurate excitation frequencies including two 
special examples. The first special example is reduced on 
the standard Lotka-Volterra equations with constant 
coefficients giving an autonomous dynamical system. 
Nevertheless, interesting properties of the competitive 
predator-prey model are demonstrated in this case: the 
exinction of all populations of predators except the 
dominant one from ecosystem viewpoint and the absolute 
necessity of the multistage homotopy perturbation in the 
long term behaviour ensuring the desired accuracy. In the 
second special example, excitation frequencies are chosen 
commensurate, excluding the phenomenon of the 
extinction and leading to the stable periodic orbit in the 
long term behaviour. 
 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
24 
 
2. Multistage Homotopy Perturbation (MH-P) Method 
for Predator-Prey Systems with Time-Dependent 
Coefficients  
The multistage homotopy perturbation (MH-P) method 
is improvement of the standard homotopy perturbation 
(SH-P) method. The SH-P method is distinguished by high 
accuracy on short time intervals.Short time intervals 
represent a serious limitation in dynamical systems, 
exhibiting chaos. An example of such a dynamical system 
is the famous Lorenz system, where the chaos is fully 
developed on the relatively long time interval.In this paper 
will be shown, that short time intervals represent a serious 
limitation in respect of the required accuracy in the 
predator-prey systems, too. Fortunately, the desired high 
accuracyof the method can be maintained even on long 
time intervals, when the entire time interval is divided into 
finite number of subintervals with equal length and the 
homotopy perturbation method is applied on each 
subinterval.In such a way, the MH-P method is obtained. 
To implement MH-P method on the some subinterval, the 
initial values at the beginning of the subinterval must be 
known, however these values are actually known at the 
beginning of the first subinterval only. To overcome this 
difficulty, the final values of the preceding subintervals are 
computed and chosen as initial conditions on the next 
subinterval to ensure the C
0
 continuity (an exception of this 
procedure represents the  first subinterval, which does not 
have a preceding subinterval at all). 
Suppose, that we are interested in the solution of Eqs. 
(1) on the time interval  
0
, tt . According to the MH-P 
method, the time interval  
0
, tt is divided intos 
subintervals of equal length t  . The solution is sought 
iteratively on subintervals  
       
0 1 1 2 2 3 1
, , , , , , , ,
s
t t t t t t t t

, where initial values on 
each particular subinterval are unknown except initial 
values on the first subinterval in the time 
*
0
tt  . The 
homotopy perturbation method on a particular subinterval 
can be formulated in the simplest way, if the folowing 
homotopy: 
 
            
1
0, 1, ,
n
i i i i i ij i j
j
L x L v pL v p t x a t x x i n 


       


 (2)
is assigned to the system of ODE's (1), where 
d
dt
L  
denotes the linear operator,   0,1 p  is an embedding 
parameter and variables 
j
v denote initial approximations 
of the solution, which must satisfy initial conditions. 
Homotopy (2) is reduced on the linear system of equations, 
when 0 p  and is transformed during the so called 
deformation process into the original nonlinear system of 
equations, when embedding parameter approaches 1 p  . 
The approximate solution of Eqs. (1) is represented in the 
form of power series:  
 
              
23
0 1 2 3
0
, 1, ,
k
i i i i i ik
k
x t x t px t p x t p x t p x t i n


        , , (3)
where     , 1,2,
ik
x t k  denote unknown functions of 
time, which must be determined in the following 
perturbation procedure.By introduction of power series (3) 
into equation of homotopy (2) one obtains: 
 
 
         
       
1
00
1
1 0 0
0, 1, ,
kk
ik i i i ik
kk
n
kl
ij ik jl
j k l
L p x t L v pL v t p x t
a t p x t x t i n






  

   


     
 (4) 
By collecting coefficients at like powers of the embedding 
parameter p one obtain the following equations for 
0
p and 
1
p , respectively: 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
25 
 
    
0
00
: 0, ,
i i i i
p L x L v x v     (5) 
              
1
1 0 0 0
1
:0
n
i i i i ij i j
j
p L x t L v t x t a t x t x t 

     

, (6) 
and for each   ,1
r
pr  the equation: 
            
1
, 1 , 1
10
:0
nr
r
ir i i r ij ik j r k
jk
p L x t t x t a t x t x t 

  

    

. (7)
The system of equations (5-7) is obtained in an explicit 
form and can be solved by successive integrations of linear 
operators     , 1,2, ,
ir
L x t r m  

.Although Eqs. (5-7) 
are developed by the perturbation procedure, the small 
perturbation parameter, which is present in all standard 
perturbation methods, is not required here at all. In other 
words, functions  
i
t  and  
ij
at , respectively, do not 
contain such a small perturbation parameter. Consequently, 
the MH-P method can be applied in the analysis of the 
predator-prey ecosystems, which contain high polynomial 
nonlinearities and time-dependent coefficients. In the rest 
of the paper, we will apply the following explicit form of 
functions  
i
t  and  
ij
at :  
 
 
   
     
0
1
,0 , , , ,
1
cos sin ,
cos sin , , 1,2, ,
t
t
M
i i iq iq iq iq
q
M
ij ij ij q ij q ij q ij q
q
t t t
a t t t i j n
     
    


   
    
 (8) 
where 
0 i
 , 
,0 ij
 ,
iq
 , 
iq
 , 
, ij q
 , 
, ij q
 areexcitation 
amplitudes and 
iq
 , 
, ij q
 ;   1,2, ,
t
qM  are 2
t
M
excitation frequencies. If all pairs of frequencies 
iq
 , 
, ij q

have an integer ratio, then frequencies 
iq
 , 
, ij q
 ;
  1,2, ,
t
qM  are commensurate and all excitations are 
periodic. If all pairs of frequencies 
iq
 , 
, ij q
 have an 
irational ratio, then frequenciesare incommensurate and 
functions   
i
t  and  
ij
at are almost-periodic. A special 
example of Eqs. (8) occurs, when there are 0
t
M  tones. 
In such a case, the predator-prey ecosystem is autonomous 
and contains the constant coefficients.  
When functions  
i
t  and  
ij
at are given in the form 
(8), then Eqs. (5-7) can be integrated analytically. 
Therefore, it is evident, that a very large class of predator-
prey ecosystems can be analytically solved by means of the 
presented MH-P method. Now look for the implementation 
of the MH-P method on each subinterval.  
Denote the starting time on an arbitrary subinterval by 
*
t and the corresponding final time of the subinterval by t. 
Then the entire time interval in SH-P method as well as the 
first subinterval in the MH-P method can be denoted by 
 
0
*, t t t  , the second subinterval of the MH-P method by 
 
1
*, t t t  and the last subinterval by  
1
*,
s
t t t

 , respectively. 
By using this convention, the entire time interval as well as 
all subintervals can be computed in the same way. As 
initial solutions on an arbitrary subinterval we choose 
initial conditions of that subinterval. Because initial 
conditions differ on subequent intervals, we use the same 
convention as above and write: 
  
   
 
* **
0
, 1, ,
i i i i
v t x t x t c i n     ,        (9) 
where 
*
i
c denote initial conditions on an arbitrary 
subinterval, which are equal to the computed values  
i
xt 
on the end of the preceding subinterval. From Eq. (9) it 
follows:  
   0
i
Lv  , (10) 
which in turn simplifies Eq. (6). As already mentioned, the 
system of Eqs. (5-7) is solved successively by using the 
inverse linear operator:   
 
   
*
1
d
t
t
Lt

  

. (11) 
Actually, integration (11) is first applied on Eq. (6), 
because Eq. (5) is fullfilled trivially. After this, the 
integration is continued on subsequent equations. 
Approximate solutions of the original Lotka-Volterra 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
26 
 
equations (1) are obtained by execution ofthe limit process 
1 p  on the power series (3), where only the first mterms 
are retained from the practical reasons:  
     
1
,
1
0
lim
m
i m i ik
p
k
t x t x t 



  (12) 
The algorithm of the MH-P method for predator-prey 
ecosystems with time-dependent coefficients is now 
completed. The algorithm is ideally suited for 
implementation in the programming environment 
Mathematica
®
. 
3. Results and discussion 
In this paper, three different competitive predator-prey 
ecosystems are analyzed by using the MH-P method. In the 
first example, the predator-prey system with three species 
is analyzed, where population growth rates and coefficients 
of mutual interactions between populationsare prescribed 
in the form of almost-periodic functions. The predator-prey 
system, which is studied in this example, is competitive, 
because in addition, coefficients of mutual interactions 
fullfill conditions   0
ij
at  .The governing Lotka-Volterra 
equations have the following form: 
 
 
           
           
               
1 1 1 1 1 2 3
2 2 2 1 2 2 3
3 3 3 1 2 3 3
5 4cos 3 cos
4 3sin 3 sin
1
3 2cos 2 4 cos 2
5
x t x t t t x t x t x t
x t x t t x t t x t x t
x t x t t x t x t t x t





      



      





     


 
 (13)
where excitation frequencies are incommensurate having 
values  
1
1 2 1 3 1
3
1, , 2         and where theinitial  
conditions
      
10 1 1 20 2 2 30 3 3
0 0.4, 0 0.4, 0 0.4 x x c x x c x x c          (14)
are prescribed. The first equation of (13) belongs to the 
prey population, the second equation belongs to the 
superpredator population and the third belongs to the 
population of predators. The meaning of incommensurate 
frequenciesisthat grow rate cycles and cycles of mutual 
interactions between population of species are completely 
independent of each other. Because Eq. (5) is trivial, it can 
be omitted and only Eqs. (6,7) are written down: 
 
 
       
               
               
       
         
             
11 1 1 10
2
1 10 10 20 10 30 1
2
21 2 2 20 20 10 2 20
1
20 30 2
31 3 3 30
2
30 10 30 20 3 30
5 4cos
3 cos 0, 0
4 3sin 3 sin
:
0, 0
3 2cos 2
1
4 cos 2
5
L x t L v t x t
t x t x t x t x t x t L v
L x t L v t x t x t x t t x t
p
x t x t L v
L x t L v t x t
x t x t x t x t t x t





   

     
      

  
   

       
3
0, 0 Lv
















 (15) 
 
           
       
         
         
1
1 1 1, 1 1 1 1, 1
0
11
1 2, 1 1 3, 1
00
1
2 2 2, 1 2 1, 1
0
11
2 2 2, 1 2 3, 1
00
5 4cos 3 cos
0
4 3sin
:
3 sin
r
r r k r k
k
rr
k r k k r k
kk
r
r r k r k
k r
rr
k r k k r k
kk
L x t t x t t x t x t
x t x t x t x t
L x t t x t x t x t
p
t x t x t x t x t




  


   


  


   

     

   
    

    
           
           
1
3 3 3, 1 3 1, 1
0
11
3 2, 1 3 3 3, 1
00
0
3 2cos 2
1
4 cos 2 0
5
r
r r k r k
k
rr
k r k k r k
kk
L x t t x t x t x t
x t x t t x t x t



  


   















    




    


 (16) 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
27 
 
System of equations (15,16) is complete, because equation 
of an arbitrary order can be constructed in the explicit 
form. By successive integration of equations, analytical 
solutions        
1
, ; 1,2,3 , 2, ,
i ir
x t x t i r m  are 
obtained, which are then used in Eq. (12) to construct 
approximate analytical solutions  
, im
t  of the m-th 
order. As an example, the integration of Eq. (15) is shown, 
which leads on the result: 
 
 
 
 
   
 
       
 
 
   
 
 
*
*
2
* * * * * * *
11 11 1 1 1 1 1 2 1 3
* * * * * * *
1 1 1 1 1 2 3
1
2
* * * * * * *
21 21 2 2 2 1 2 2 2 3
**
2 2 2
2
5 4cos 3 cos d
1
4 sin sin 5 3 ,
4 3sin 3 sin d
1
3 cos cos
t
t
t
t
x t x t c c c c c c
c c t t c c c t t
x t x t c c c c c c
c c t
    


    



      




       



      



   
     
 
 
       
 
 
     
*
* * * * *
2 1 2 3
2
* * * * * * *
31 31 3 3 3 1 3 2 3 3
* * * * * * *
3 3 3 3 1 2 3
3
4 3 ,
1
3 2cos 2 4 cos 2 d
5
1 1 1
1 sin 2 sin 2 3 4
25
t
t
t c c c t t
x t x t c c c c c c
c c t t c c c t t
    















    






      
 





   

       
    

   
 
, (17) 
where is considered, that relations 
     
* * *
11 21 31
0 x t x t x t    always hold, because initial 
conditions are fulfilled already by values 
           
* * * * * * * * *
10 1 1 20 2 2 30 3 3
,, x t x t c x t x t c x t x t c      
. 
On the same way, the integration of subsequent equations 
of higher order is carried out. Results of computation of 
populations evolution by using the MH-P method and 
comparison with Runge-Kutta method are shown on 
Figures 1, 2 and 3, respectively, where an excellent 
agreement is established. The high accuracy of the MH-P 
method in the entire observation range is achieved with 
relatively small number (m=4) of terms in power series 
(12), but using a great number of subintervals s=500.  To 
evaluate the accuracy of results, one computes the integral 
of the absolute error between two subequent solutions 
 
 
1 s
i
xt

and 
 
    , 1, , 3
s
i
x t i n  over the entire time 
interval  
0
, tt , where the first solution contains s-1 
subintervals and the next solution has s subintervals, 
respectively and check, if the values of the computed 
integral are less than the prescribed tolerance tol: 
 
 
 
   
0
1
d , 1,2,3
t
ss
i i i
t
J x t x t t tol i

   

 (18) 
If this check fails, then the number of subintervals is too 
small. In the predator-prey system (13), the tolerance tolis 
prescribed to be equal 
4
10 tol

 . The number of 
subintervals s=500 is then sufficient to fulfill the criterion 
(18) on the entire time interval  
0
0, 50 tt  with 
computed values of integrals:
-5 -5
12
4.72776 10 , 1.53379 10 JJ     and
-5
3
8.23404 10 J . On the Figure 4, the parametric 
dependencies of evolutions of populations    
12
, x t x t 
and  
3
xt , which are computed by MH-P method and 
comparison of results with Runge-Kutta method are 
plotted, where again the excellent agreement is shown.
 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
28 
 
Figure 1. Evolution of populationx
1
(t) in the almost-periodic case.        Solution, computed by MH-P method,                
         Solution, computed by numerical integration using Runge-Kutta method. 
Figure 2. Evolution of populationx
2
(t) in the almost-periodic  case.        Solution, computed by MH-P 
method, - - -  Solution, computed by numerical integration using Runge-Kutta method. 
Figure 3.  Evolution of populationx
3
(t) in the almost-periodic case.       Solution, computed by MH-P 
method,- - -Solution,computed by numerical integration using Runge-Kutta method. 
Figure 4. Parametric dependencies of populationsx
1
(t), x
2
(t) in x
3
(t) in the almost-periodic case.               Solution, 
computed by MH-P method, - - - Solution, computed by numerical integration using Runge Kutta method. 
Time histories on Figures 1-3, as well as orbits on the 
Figure 4 reveal the essential nature of almost-periodic 
oscillations. Time histories consist from series of cycles, 
which do not repeat. This property is a direct consequence 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
29 
 
of incommensurate frequencies involved. On the other 
side, it can be recognized, that the time histories are not 
random or even chaotic. This property is evident also from 
the parametric plots, where numerous loops can be seen, 
but no random behaviour. 
Now look on the periodic counterpart of the predator-
prey ecosystem (13). To obtain periodicity, the general 
form of Eqs. (13) and corresponding values of initial 
conditions are retained, but excitation frequencies are now 
required to be commensurate. To obtain the 
commensurability of frequencies, we prescribe an 
additional condition in the form 
1 2 3
1       , which of 
course is restrictive and therefore it must be justified in 
reality. The results of computation by application of the 
MH-P method and comparison with results of numerical 
integration using Runge-Kutta method are plotted on 
Figures 5, 6 and 7, where again an excellent agreement can 
be found.The convergence test with prescribed tolerance 
4
10 tol

 results in values of integrals: 
-5 -5
12
3.05072 10 , 8.46863 10 JJ     and
-5
3
5.82898 10 J . The parametric dependencies of 
evolutions of populations    
12
, x t x t and  
3
xt , which 
are computed by MH-P method and comparison of results 
with Runge-Kutta method are shown on Figure 8, where 
again an excellent agreement is obtained. 
Unlike time histories of populations in almost-periodic 
case, oscillations on Figures 5-7 after the transient phase 
now consist from repeating cycles, which confirm the 
periodicity. The same conclusion follows also from the 
Figure 8, where trajectories approaches the periodic orbits 
after the transient phenomenon dies. 
Figure 5. The evolution of populationx
1
(t) in the periodic case.       Solution, computed by MH-P 
method,          Solution, computed by numerical integration using Runge-Kutta method. 
 
Figure 6. The evolution of populationx
2
(t) in the periodic case.          Solution, computed by MH-P method, - 
- -  Solution, computed by numerical integration using Runge-Kutta method. 
 
 
 
 
 
 
 
 
 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
30 
 
 
Figure 7. The evolution of populationx
3
(t) in the periodic case.            Solution, computed by MH-P method, 
- - - Solution,computed by numerical integration using Runge-Kutta method. 
 
Figure 8. Parametric dependencies of populationsx
1
(t), x
2
(t) in x
3
(t) in the periodic case.         Solution, 
computed by MH-P method, - - - Solution, computed by numerical integration using Runge-Kutta method
As third example, the competitive predator-prey 
ecosystem (13) is analyzed, which has constant coefficients 
only. Therefore, all time-dependent terms in Eqs. (13) are 
omitted and corresponding Lotka-Volterra equations are 
reduced on the form: 
         
         
         
1 1 1 2 3
2 2 1 2 3
3 3 1 2 3
53
43
1
34
5
x t x t x t x t x t
x t x t x t x t x t
x t x t x t x t x t


    



    





   


 
 (19) 
The computation by MH-P method is performed with the 
same initial conditions as in first two examples. The results 
of computation by application of the MH-P method are 
plotted on the Figure 9. The comparison with results of 
numerical integration using Runge-Kutta method is shown 
and an excellent agreement is obtaned again. The 
convergence test with prescribed tolerance 
4
10 tol

 
results in values of integrals: 
-5 -5
12
8.35791 10 , 3.40169 10 JJ     and
-6
3
2.45845 10 J . 
 
Figure 9. The evolution of populationx
3
(t) in the periodic case.            Solution, computed by MH-P method, 
- - -  Solution,computed by numerical integration using Runge-Kutta method. 
 
 
 
 
 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
31 
 
Time histories on the Figure 9 reveal, that first two species, 
which belong to the prey and predator, respectively, 
survive. Their populations grow according to the logistic 
curve, approaching the corresponding carrying capacities. 
On the contrary, the third species of superpredator is after 
some growth in a short initial phase driven to the 
extinction.  
4. Conclusions 
In this paper, a generalized predator-prey model with 
time-dependent coefficients of populations growth rates 
and time-dependent coefficients of mutual interactions 
between populations of speciesis treated, where 
competitive ecosystems with almost-periodic, periodic and 
constant coefficients are analyzed in details. To obtain 
approximate analytical solutions of generalized Lotka-
Volterra equations for almost-periodic, periodic and 
constant coefficients, respectively, the MH-P metod is 
developed in an explicit form, which is appropriate for 
analysis on long time intervals. To ensure the desired 
accuracy of the MH-P method on the entire time interval, 
the  useful convergence test of the subsequent solutions is 
provided.Nevertheless, the results computed by MH-P 
method are in all examples compared with results of the 
numerical integration using Runge-Kutta method, where an 
excellent agreement is always obtained. The analytical 
form of solutions has many advantages, which allow 
analysis of stability, parametric studies and bifurcations. 
Such possibilities as well as the study of the chaotic 
behaviour in predator-prey ecosystems [9] are challenging 
for the future research work. Besides this, the MH-P 
method is well suited for computation of highly nonlinear 
ecosystems. 
References 
1. Hoppensteadt, F.Predator-prey model. 
Scholarpedia1, 2006, pp.1563. 
2. Hoppensteadt, F. C.Mathematical methods of 
population biology. Courant Institute of 
Mathematical Sciences New York University, New 
York, 1977. 
3. Smith, H. L.Periodic orbits of competitive and 
cooperative systems.Journal of Differential 
Equations65, 1986, pp. 361-373. 
4. He, J. H. A coupling method of homotopy technique 
and perturbation technique for nonlinear 
problems,Int. J. Non-Linear Mech.35(1),2000, pp. 
37-43. 
5. Pušenjak, R. R., Oblak, M. M., Tičar, I. 
Nonstationary Vibration and Transition through 
Fundamental Resonance of Electromechanical 
Systems Forced by a Nonideal Energy Source. 
International Journal for Nonlinear Sciences and 
Numerical Simulation, vol. 10, no. 5, May 2009, pp. 
637-660. 
6. Pušenjak, R. R., Oblak, M. M., Tičar, I. Modified 
Lindstedt-Poincare method with multiple time scales 
for combination resonance of damped dynamical 
systems with strong non-linearities. Int. j. of 
nonlinear sci. & numer. simul., vol. 11, no. 3,2010, 
pp. 173-201. 
7. Chowdhury, M. S. H.,Hashim, I., Momani, S. The 
multistage homotopy-perturbation method: A 
powerful scheme for handling the Lorenz system. 
Chaos, Solitons and Fractals 40,2009, pp. 1929-
1937. 
8. Pušenjak, R., Oblak, M., Lipičnik, M.Evolution of a 
Predator-Prey System with Periodic Excitation by 
Using Multistage Homotopy Perturbation 
Method.Kuhljevi dnevi, 26-27. September 2012, 
Rogaška Slatina.(accepted for publication). 
9. Vano,J. A., Wildenberg, J. C., Anderson,M. B., 
Noel,J. K. and Sprott, J. C.Chaos in low-dimensional 
Lotka-Volterra models of 
competition.Nonlinearity19, 2006, pp. 2391-2404. 

Anali PAZU - Letnik 2, leto 2012, številka 1 
 
32