BULLETII~ OF MATHEMATICALBIOLOGY
VOLUME 37, 1975
APPLICATION OF PERTURBATION THEORY TO THE NONLINEAR VOLTERRA-GAUSE-WITT MODEL FOR P R E Y - P R E D A T O R INTERACTION
9 RANABIRDUTT, P. K. GHOSHAND B. B. KARM~AR Department of Physics, Visva-Bharati, Santiniketan-731235, W. Bengal, India
Krylov-Bogoliubov-Mitropolsky perturbation method was used to study the effect of nonlinearity in the V o l te r r a - G a u s e - W it t (VGW) model for a two species prey-predator system. The first order corrections to both the frequency of oscillation and the amplitude of the linearized system were computed. I t was found t h a t the basic qualitative features of the nonlinearity are exhibited by the first order result. We have also discussed the Lotka-Volterra problem which is a special case of VGW model.
One of the most interesting observations in the study of population dynamics is that the populations of a pair of species bound by the prey-predator relationship, fluctuate periodically. The mathematical investigation of this problem was initiated by Lotka and Volterra (Goel alld Maitra, 1971; Bartlett, 1960; iGause, 1971) who proposed a set of rate equations with quadratic nonlinearities which describe phenomenologieally the interaction of a prey species with its predator when both populations coexist in an ecological niche with finite resources. It is assumed that the prey (Species 1) would grow exponentially in the absence of the predator (Species 2) while the predator dies out exponentially in the absence of its prey. The Lotka-Volterra equations 1. Introduction.
are
dN~/dt = a l N ~ - fl~N1N2 dN2/dt = - a 2 N 2 + [3zN1N2
(1) 139
140
R A N A B I R D U T T , P. K. G H O S H A N D B. B. KA_RMAKAR
where N, (i = 1, 2) is the number of individuals of the species i at a given time, a t is the innate capacity for increase per individual (intra-specific parameter) and fl, is the coefficient of mutual interaction (interspecific or niche overlap parameter). I n this form of the equations describing a prey-predator system, these coefficients are all positive. A deficiency of the Lotka-Volterra (LV) model is the non-existence of a saturation level for the population of the prey species alone (in the absence of the predator) which, in general does not grow indefinitely in a given environment with limited space and resource. The saturation effect was later incorporated by Gause and Witt who introduced a self-interaction term in the Verhurlst form. The modified Volterra-Gause-Witt (VGW) model (Gause and Witt, 1935) is described by the instantaneous growth equations dN1/dt = alNl(1 - N1/O) - flIN1N2 dN2/dt = - a 2 N 2 + fi2N1N2
(2)
where 0 is the carrying capacity (self saturation level) of Species 1. I t is clear that the VGW model reproduces the LV system in the limit 0-+~. It is quite difficult to handle these nonlinear population models analytically because the exact solutions have not been obtained so far. In the case of the LV model there is an advantage: the equations (1) describe a conservative system and correspond to a simple closed loop in the (N1N2) phase plane. This certainly indicates that periodic solutions exist for N1 and Ng. with a single unique frequency. However, it is not at all clear whether the VGW model will simulate similar periodic fluctuations in these variables. When the exact solutions are not available, the standard procedure generally followed, is to linearize the nonlinear equations in the neighborhood of the equilibrium point in the phase plane, assuming a priori that the effect of the nonlinear terms would be small (Rosen, 1970; Minorsky, 1962). Essentially, this assumption motivates us to work out a perturbative calculation to determine the weak nonlinear effects. In Section 2, we shall give an outline of the prescription for linearization (Rosen, 1970) and obtain the analytic form of the amplitude and frequency of oscillation. I n Section 3, we shall discuss how the asymptotic perturbation theory of Krylov-Bogoliubov-Mitropolsky (KBM, cf. Minorsky, 1962) m a y be applied to our problem to estimate the first order correction, due to nonlinearity, to the linear frequency and to the amplitude. Some remarks on the first order results will be given in the concluding section. 2. Linearization of the Problem.
As a guideline, we shall utilize the lineariza-
A P P L I C A T I O N OF P E R T U R B A T I O N T H E O R Y
141
tion theorem (l~osen, 1970) giving sufficient conditions which, when satisfied will allow us to draw conclusions about the general feature of the nonlinear model as well as about the stability of the equilibrium point from the behavior of the linearized system. The nonzero equilibrium populations as obtained from (2) b y setting dN1,2/dt = O, are q~ %/~ q~ = ~ 1 / ~ ( 1 - ~//3~0).
(3)
=
For our purpose it is convenient to shift the origin of the phase plane to the point (ql, q2) and hence we substitute N~(t) = ql + x(t) N2(t) = q2 + y(t)
(4)
in (2) to obtain dx/dt = - alq---Ax - fllqlY + r 0 dy/dt = fi2q2x + r y)
y)
(5)
where (I)l(x, y) =
al x2 0
fllxy
(6)
~2(x, y) = fl2xy.
The equilibrium point (0, 0) in the xy-plane corresponds to the (ql, q2) point in the Zr12V2-plane. The nonlinear functions ~1, and (~2 satisfy the conditions r
o) = ~2(o, o) = o
Lt (b~(x,y) = Lt (P2(x,y) = 0 x, y-,o V ( x 2 + y2) x, v-,o V ( x2 + y2)
(7)
which are required for the validity of the linearization theorem. Thus in the neighborhood of the steady state populations, *he behavior of (5) would be represented approximately by the solution of the linear par~. For the linearized system, we then obtain the following equation from (5) ~ +--a?+
k2x = 0
or
(s) (zlqt..
ij + - - ~ y •
k~y = 0
in which k2 - fllfl2qlq2 = ala2 ( 1 -
a~20)"
(9)
142
RANABIRDUTT, P. K. GHOSH AND B. B. KARMAKAR
Equations (8) can be solved exactly and it is known from the solutions that x(t) and y(t) will have periodic fluctuation if the condition 4/~2 > (~qi/9) 2 is satisfied. For biological species, it is found that alql[O is, in general, a small quantity as compared to 2k and we shall consider that the system fluctuates effectively with the approximate frequency/r Thus the solution for x is written
in which ao and ~ are the arbitrary constants. Since x and y are related b y (5), the solution for y can be obtained from the knowledge of x and ~. Henceforth we consider the equation for x only. The reason we write the solutions with the approximate frequency/~ and not with the exact one, which is w = kV'[1 - 1/41c2(alql/O)2] would be clear in the next section where we shall discuss the perturbation method. However, we retain (alql[O) in the exponential since it distinguishes a damped system from an undamped one. From the topological consideration, the damping term in (10) indicates that the nonvanishing equilibrium point in the VGW model is a focus type singularity which, in our case, is asymptotically stable. The frequency and the amplitudes of oscillation in the LV model can be trivially obtained from (9) and (10) b y taking the limit (~--> ~ . With this limit, we find that the linearized LV model exhibits undamped oscillation with the frequency ]co = %/(ala2), which does not include the interspecific parameters fil or fie- Further, the singularity in this case becomes a center which is neutrally stable (i.e. stable b u t not asymptotically stable) (Rosen, 1970; Minorsky, 1962).
3. Nonlinear Corrections by the K B M Perturbation Technique.
Our task is now to obtain the corrections due to the nonlinear terms in r which have been neglected in the linearization procedure. For this purpose, the K B M asymptotic perturbation method seems to be an elegant prescription. This method has been applied successfully b y several authors (Ford, 1961; Ford and Waters, 1963; Mendelson, 1970) to nonlinear oscillation problems with or without damping. Following their technique, we reduce (5) and (6) to the canonical forms: + k2x = fl(x, (11) where fl(x,
=
(12} ~2 +q1+x+~
~2
X3 9
APPLICATION OF PERTURBATION
THEORY
143
Here ~ is a dimensionless small positive expansion parameter 9 I t is to be emphasized t h a t the perturbation procedure followed here is formal in the sense t h a t there is no w a y of estimating the size of the expansion parameter. W e note t h a t the linear damping terms have been included in f i in accordance with our starting assumption t h a t the coefficient alqi/8 is a small q u a n t i t y of the same order of E. This is done because we need u n d a m p e d periodic solution in the zero th order approximation (e = 0) in the K B M method. W h e n e r 0, (11) is assumed to have a series solution
x(t) = a cos ~b + EU
VOLUME 37, 1975
APPLICATION OF PERTURBATION THEORY TO THE NONLINEAR VOLTERRA-GAUSE-WITT MODEL FOR P R E Y - P R E D A T O R INTERACTION
9 RANABIRDUTT, P. K. GHOSHAND B. B. KARM~AR Department of Physics, Visva-Bharati, Santiniketan-731235, W. Bengal, India
Krylov-Bogoliubov-Mitropolsky perturbation method was used to study the effect of nonlinearity in the V o l te r r a - G a u s e - W it t (VGW) model for a two species prey-predator system. The first order corrections to both the frequency of oscillation and the amplitude of the linearized system were computed. I t was found t h a t the basic qualitative features of the nonlinearity are exhibited by the first order result. We have also discussed the Lotka-Volterra problem which is a special case of VGW model.
One of the most interesting observations in the study of population dynamics is that the populations of a pair of species bound by the prey-predator relationship, fluctuate periodically. The mathematical investigation of this problem was initiated by Lotka and Volterra (Goel alld Maitra, 1971; Bartlett, 1960; iGause, 1971) who proposed a set of rate equations with quadratic nonlinearities which describe phenomenologieally the interaction of a prey species with its predator when both populations coexist in an ecological niche with finite resources. It is assumed that the prey (Species 1) would grow exponentially in the absence of the predator (Species 2) while the predator dies out exponentially in the absence of its prey. The Lotka-Volterra equations 1. Introduction.
are
dN~/dt = a l N ~ - fl~N1N2 dN2/dt = - a 2 N 2 + [3zN1N2
(1) 139
140
R A N A B I R D U T T , P. K. G H O S H A N D B. B. KA_RMAKAR
where N, (i = 1, 2) is the number of individuals of the species i at a given time, a t is the innate capacity for increase per individual (intra-specific parameter) and fl, is the coefficient of mutual interaction (interspecific or niche overlap parameter). I n this form of the equations describing a prey-predator system, these coefficients are all positive. A deficiency of the Lotka-Volterra (LV) model is the non-existence of a saturation level for the population of the prey species alone (in the absence of the predator) which, in general does not grow indefinitely in a given environment with limited space and resource. The saturation effect was later incorporated by Gause and Witt who introduced a self-interaction term in the Verhurlst form. The modified Volterra-Gause-Witt (VGW) model (Gause and Witt, 1935) is described by the instantaneous growth equations dN1/dt = alNl(1 - N1/O) - flIN1N2 dN2/dt = - a 2 N 2 + fi2N1N2
(2)
where 0 is the carrying capacity (self saturation level) of Species 1. I t is clear that the VGW model reproduces the LV system in the limit 0-+~. It is quite difficult to handle these nonlinear population models analytically because the exact solutions have not been obtained so far. In the case of the LV model there is an advantage: the equations (1) describe a conservative system and correspond to a simple closed loop in the (N1N2) phase plane. This certainly indicates that periodic solutions exist for N1 and Ng. with a single unique frequency. However, it is not at all clear whether the VGW model will simulate similar periodic fluctuations in these variables. When the exact solutions are not available, the standard procedure generally followed, is to linearize the nonlinear equations in the neighborhood of the equilibrium point in the phase plane, assuming a priori that the effect of the nonlinear terms would be small (Rosen, 1970; Minorsky, 1962). Essentially, this assumption motivates us to work out a perturbative calculation to determine the weak nonlinear effects. In Section 2, we shall give an outline of the prescription for linearization (Rosen, 1970) and obtain the analytic form of the amplitude and frequency of oscillation. I n Section 3, we shall discuss how the asymptotic perturbation theory of Krylov-Bogoliubov-Mitropolsky (KBM, cf. Minorsky, 1962) m a y be applied to our problem to estimate the first order correction, due to nonlinearity, to the linear frequency and to the amplitude. Some remarks on the first order results will be given in the concluding section. 2. Linearization of the Problem.
As a guideline, we shall utilize the lineariza-
A P P L I C A T I O N OF P E R T U R B A T I O N T H E O R Y
141
tion theorem (l~osen, 1970) giving sufficient conditions which, when satisfied will allow us to draw conclusions about the general feature of the nonlinear model as well as about the stability of the equilibrium point from the behavior of the linearized system. The nonzero equilibrium populations as obtained from (2) b y setting dN1,2/dt = O, are q~ %/~ q~ = ~ 1 / ~ ( 1 - ~//3~0).
(3)
=
For our purpose it is convenient to shift the origin of the phase plane to the point (ql, q2) and hence we substitute N~(t) = ql + x(t) N2(t) = q2 + y(t)
(4)
in (2) to obtain dx/dt = - alq---Ax - fllqlY + r 0 dy/dt = fi2q2x + r y)
y)
(5)
where (I)l(x, y) =
al x2 0
fllxy
(6)
~2(x, y) = fl2xy.
The equilibrium point (0, 0) in the xy-plane corresponds to the (ql, q2) point in the Zr12V2-plane. The nonlinear functions ~1, and (~2 satisfy the conditions r
o) = ~2(o, o) = o
Lt (b~(x,y) = Lt (P2(x,y) = 0 x, y-,o V ( x 2 + y2) x, v-,o V ( x2 + y2)
(7)
which are required for the validity of the linearization theorem. Thus in the neighborhood of the steady state populations, *he behavior of (5) would be represented approximately by the solution of the linear par~. For the linearized system, we then obtain the following equation from (5) ~ +--a?+
k2x = 0
or
(s) (zlqt..
ij + - - ~ y •
k~y = 0
in which k2 - fllfl2qlq2 = ala2 ( 1 -
a~20)"
(9)
142
RANABIRDUTT, P. K. GHOSH AND B. B. KARMAKAR
Equations (8) can be solved exactly and it is known from the solutions that x(t) and y(t) will have periodic fluctuation if the condition 4/~2 > (~qi/9) 2 is satisfied. For biological species, it is found that alql[O is, in general, a small quantity as compared to 2k and we shall consider that the system fluctuates effectively with the approximate frequency/r Thus the solution for x is written
in which ao and ~ are the arbitrary constants. Since x and y are related b y (5), the solution for y can be obtained from the knowledge of x and ~. Henceforth we consider the equation for x only. The reason we write the solutions with the approximate frequency/~ and not with the exact one, which is w = kV'[1 - 1/41c2(alql/O)2] would be clear in the next section where we shall discuss the perturbation method. However, we retain (alql[O) in the exponential since it distinguishes a damped system from an undamped one. From the topological consideration, the damping term in (10) indicates that the nonvanishing equilibrium point in the VGW model is a focus type singularity which, in our case, is asymptotically stable. The frequency and the amplitudes of oscillation in the LV model can be trivially obtained from (9) and (10) b y taking the limit (~--> ~ . With this limit, we find that the linearized LV model exhibits undamped oscillation with the frequency ]co = %/(ala2), which does not include the interspecific parameters fil or fie- Further, the singularity in this case becomes a center which is neutrally stable (i.e. stable b u t not asymptotically stable) (Rosen, 1970; Minorsky, 1962).
3. Nonlinear Corrections by the K B M Perturbation Technique.
Our task is now to obtain the corrections due to the nonlinear terms in r which have been neglected in the linearization procedure. For this purpose, the K B M asymptotic perturbation method seems to be an elegant prescription. This method has been applied successfully b y several authors (Ford, 1961; Ford and Waters, 1963; Mendelson, 1970) to nonlinear oscillation problems with or without damping. Following their technique, we reduce (5) and (6) to the canonical forms: + k2x = fl(x, (11) where fl(x,
=
(12} ~2 +q1+x+~
~2
X3 9
APPLICATION OF PERTURBATION
THEORY
143
Here ~ is a dimensionless small positive expansion parameter 9 I t is to be emphasized t h a t the perturbation procedure followed here is formal in the sense t h a t there is no w a y of estimating the size of the expansion parameter. W e note t h a t the linear damping terms have been included in f i in accordance with our starting assumption t h a t the coefficient alqi/8 is a small q u a n t i t y of the same order of E. This is done because we need u n d a m p e d periodic solution in the zero th order approximation (e = 0) in the K B M method. W h e n e r 0, (11) is assumed to have a series solution
x(t) = a cos ~b + EU
